Dear Amber users,
This refers to all nmr group restraints, excluding distance restraints,
that have a zero energy area in the middle of the parabola (r2=r3). The
restraint forces seem to be not correctly set to zero in this area.
For example, see subroutine angnrg. In case (thet <= r3), the energy is
set to zero, but not the dfr force array. Then the values are returned
(omitting the later assignment of dfr):
nmr.F90 line 173:/
// if (thet < r1) then//
// dif1 = r1-r2//
// df = 2.0d0 * k2 * dif1//
// e = df * (thet-r1) + k2*dif1*dif1//
// else if (thet < r2) then//
// dif = thet - r2//
// df = 2.0d0 * k2 * dif//
// e = k2*dif*dif//
// else if (thet <= r3) then//
// e = 0.0d0//
// return//
// else if (thet < r4) then//
// dif = thet - r3//
// df = 2.0d0 * k3 * dif//
// e = k3*dif*dif//
// else//
// dif1 = r4-r3//
// df = 2.0d0 * k3 * dif1//
// e = df * (thet-r4) + k3*dif1*dif1//
// end if/
For group-restraints (usecom), after calling angnrg, the subroutine
nmrcmf calculates the new forces f using dfr which was not set to zero
properly by angnrg:
nmr.F90 line 3681:
/case (2) ! Angle restraint//
// if (nave(2) > 0 .and. nmrfty(i) == 0) iave = iavtyp(2)//
// if (usecom) then//
// call nmrcms(x,xcom,nmrat,nmrcom,rmstot,rimass,ricmmass,i)//
// call angnrg(xcom,f,dfr,rint,0,3,6, &//
// e,r1,r2,r3,r4,rk2,rk3,ntb,aave(2*i-1), &//
// aave0(i),nave(2),nexact(2),ipower(2), &//
// tauave(2),ravein(2),dt,navint(2),iave, &//
// incflg(2),ifirst,3)//
// call nmrcmf(f,dfr,nmrat,nmrcom,rmstot,rimass,i)//
// else//
// call angnrg(x,f,dfr,rint,nmrat(1,i),nmrat(2,i),nmrat(3,i), &//
// e,r1,r2,r3,r4,rk2,rk3,ntb,aave(2*i-1), &//
// aave0(i),nave(2),nexact(2),ipower(2), &//
// tauave(2),ravein(2),dt,navint(2),iave, &//
// incflg(2),ifirst,0)//
// end if/
Is this a known issue?
For distance restraints, the dfr array is properly set zero.
Regards,
Fabian
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jan 10 2013 - 06:30:04 PST
