[MMTK] getting hessian matrix (and equation) and force vector of a protein conformaion
Konrad Hinsen
research at khinsen.fastmail.net
Wed Jan 11 11:36:19 UTC 2012
Hyuntae Na writes:
> My first thought about the second derivative of the (likely say) van der waals
> (vdW) potential energy was like the follows:
> $$
> d^2E / d dx^2 = eps_ij * (12*13*r_ij^12/ r^12 - 2*6*7*r_ij^6/r^6) * (dx/r)^
> 2,
> d^2E / d dx dy = eps_ij * (12*13*r_ij^12/ r^12 - 2*6*7*r_ij^6/r^6) * (dx/r)*
> (dy/r),
> $$
> where eps_ij = sqrt(eps_i * eps_j), r_ij = vdW_rad_i + vdW_rad_j, dx (and dy)
> is the difference of x(and y)-coordinate of atoms i and j, and r=sqrt(dx^2+dy^
> 2+dz^2).
I don't have the exact expressions at hand, but this looks wrong. The diagonal
terms (xx, yy, zz) are sums of two contributions with different powers of r.
> However, after running a symbolic operation tool (Mathematica) to check its
> validity, the program returns a more complex equations:
> $$
> d^2E / d dx^2 = eps_ij * (168*r_ij^12/ r^12 - 96*r_ij^6/r^6) * (dx/r)^2 +
> (-12*r_ij^12/r^12 + 12*r_ij^6/r^6)*(1/r)^2,
> d^2E / d dx dy = eps_ij * (168*r_ij^12/ r^12 - 96*r_ij^6/r^6) * (dx/r)*(dy/
> r).
> $$
This looks wrong as well: the two terms in the sum for d^2E / d dx^2 have
different dimensions.
Note also that the elements of the Hessian are derivatives with
respect to the individual coordinates x, y, z of each atom, not the
relative coordinates that you call dx, dy, dz.
I suggest to start by deriving the general form of the second derivatives
for a generic pair potential of the form f(r). That will be a rather compact
expression in terms of the first and second derivative of f(r). You can then
use it for vdW, Coulomb, and bond-length interactions. Angles and Dihedrals
are much messier, unfortunately.
> Additionally, would you tell me how to get the hessian matrix and forces of
> each atoms of a protein using MMTK?
See the last message from this list; it contains all you need.
You can also have a look at MMTK's source code. The code that adds
the contribution from one pair potential to the total Hessian is at
line 462 of MMTK_forcefield.c. The individiual force field terms
call that code with the values of the first and second derivative
of the potential function. For the vdW term, they are calculated
in the code at line 483 of nonbonded.c.
Good luck,
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/
---------------------------------------------------------------------
More information about the mmtk
mailing list