dihedral restraint force constant not converted from degrees
The function of the dihedral restraints is currently a bit odd, and probably should be considered a bug:
phi 0 and dphi are read from the topology file in degrees and converted internally to radians, but the
force constant (kfac) is used as read. In other words, to get the expected behavior, one would need to
specify phi and dphi in the topology file in units of degrees, but the force constant in units of kJ/rad^2
(rather than the expected kJ/degree^2).
Probably this should be fixed so that expected input units are consistent. The relevant lines are in
phi0 = ip[type].dihres.phi*d2r;
dphi = ip[type].dihres.dphi*d2r;
kfac = ip[type].dihres.kfac*fc;
Sorry, I should have caught this before. I finally figured it out because I realized I was using units of kJ/
rad^2 for my input force constants but the input angles were in degrees. I had adjusted the force
constant to get the restraint energy I was expecting, so I didn't notice before.
#1 Updated by David van der Spoel about 14 years ago
This is by design, and is the same for all other angle and dihedral forces.
Hence it is not a bug. I agree that it is not very elegant, it would be very
messy to change it however, because it would imply reading interpreting the
force constants in another way. I added a small remark in the manual in the
#2 Updated by David Mobley about 14 years ago
OK, great. I realized after I submitted it that the angle force constant has the same units, so I can see why
it is a bit confusing. Fixing the documentation should solve the problem, since essentially the problem
was just that the 3.2 manual doesn't even mention dihedral restraints, so it's not clear what units one
should expect things to be in.