[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