[MMTK] Force constants in CalphaFF
Konrad Hinsen
research at khinsen.fastmail.net
Wed Jul 3 07:54:13 UTC 2013
Lars Skjærven writes:
> I suspect my error might be related to the force constant calculation. For the CalphaFF
> I use equation 17 in your paper "Harmonicity in slow protein dynamics":
> -----
> if( r < .4)
> k = (8.6*(10**5)*r) - (2.39*(10**5))
> else
> k = 128 * r**(-6)
> -----
> Is this also the implementation in MMTK? (sorry - couldn't find it in the code).
Yes, that's the formula used in MMTK. You can find it in Src/MMTK_deformation.c
starting from line 570:
if (r_sq < 0.16) {
double r = sqrt(r_sq);
deriv2 = (860000.*r-239000.)*self->param[1];
}
else {
deriv2 = 128.*self->param[1]/(r_sq*r_sq*r_sq);
}
> Is there a way I can print these forceconstants for all pairs of ca-atoms (for my own
> debugging purpose)? I know the energyAndForceConstants(), but this gives the entire
> Hessian matrix..
True, but you can then just print the entries you need:
_, fc = universe.energyAndForceConstants()
for a1 in calphas:
for a2 in calphas:
print a1
print a2
print fc[a1, a2]
Or you can just make up a two- or three-atom system for testing and
look at the complete Hessian.
> For the record I'm comparing my modes with EnergeticModes() and EnergeticMode.rawMode()
> to avoid temperature scaling and mass-weighting:
That looks like the right approach.
Another useful approach for testing is to define your own ENM using
HarmonicDistanceRestraint terms:
from MMTK.ForceFields.ForceField import CompoundForceField
from MMTK.ForceFields.Restraints import HarmonicDistanceRestraint
def enm_forcefield(params):
terms = []
for a1, a2 in it.combinations(universe.atomIterator(), 2):
d = universe.distance(a1, a2)
if d < 2.5:
if d < 0.4:
k = params[0]
elif d < 2.5:
d3 = d*d*d
d6 = d3*d3
k = params[1]/d6
terms.append(HarmonicDistanceRestraint(a1, a2, d, k))
return CompoundForceField(*terms)
...
universe.setForceField(forcefield(params))
Note that this is not exactly the CalphaForceField, but a variant I
have used for evaluation purposes. But you can quickly make any ENM you
like using this approach.
Konrad.
--
---------------------------------------------------------------------
Konrad Hinsen
Centre de Biophysique Moléculaire, CNRS Orléans
Synchrotron Soleil - Division Expériences
Saint Aubin - BP 48
91192 Gif sur Yvette Cedex, France
Tel. +33-1 69 35 97 15
E-Mail: research AT khinsen DOT fastmail DOT net
http://dirac.cnrs-orleans.fr/~hinsen/
ORCID: http://orcid.org/0000-0003-0330-9428
Twitter: @khinsen
---------------------------------------------------------------------
More information about the mmtk
mailing list