From: Floris Buelens (floris_buelens_at_yahoo.com)
Date: Wed Dec 20 2006 - 04:53:12 CST
Hi, I have a question about NAMD's constraints implementation. I'm trying to replicate exactly the forces that are generated by NAMD's constraints using TCLForces, so that both the forces applied and the energies reported can be compared directly. I found the NAMD code that does constraints and have re-implemented it in TCL, but there are some things I'm not clear about. Specifically, with reference to this piece of the NAMD (2.6) code (where k = constraint force constant, r is distance of atom from reference, and r2 is distance squared):
value=k;
// Loop through and multiple k by r consExp times.
// i.e. calculate kr^e
// I know, I could use pow(), but I don't trust it.
for (int k=0; k<consExp; ++k)
{
value *= r;
}
// Add to the energy total
energy += value;
// Now calculate the force, which is ekr^(e-1). Also, divide
// by another factor of r to normalize the vector before we
// multiple it by the magnitude
value *= consExp;
value /= r2;
Rij *= value;
//iout << iINFO << "restraining force" << Rij << "\n";
f[localID] += Rij;
netForce += Rij;
virial += outer(Rij,vpos);
In English: take force constant k, and multiply by distance squared (assuming constraint exponent of 2). This gives the constraint energy: add this to the running total.
The applied force is then derived by taking that same energy, multiplying by constraint exponent and dividing by r^2, so the distance vector Rij is scaled by 2 * k (assuming default exponent of 2) to give the force vector.
I've blindly implemented this in my TCLForces code but I'd be grateful if someone could tell me I'm thinking along the right lines - in particular I was surprised that the force applied is in fact linearly proportional to the distance from reference co-ordinates, whereas the energy scales by r^2... or am I missing something? And most importantly for my purposes, will the reported energy be comparable, (in the same units etc), with the rest of the system energy?
Thanks a lot!
Floris Buelens
Department of Crystallography, Birkbeck College, London
___________________________________________________________
All New Yahoo! Mail – Tired of Vi_at_gr@! come-ons? Let our SpamGuard protect you. http://uk.docs.yahoo.com/nowyoucan.html
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:44:16 CST