[MMTK] getting hessian matrix (and equation) and force vector of a protein conformaion
Hyuntae Na
htna at hotmail.com
Wed Jan 11 19:51:22 UTC 2012
Hi.
Sorry about my inaccurate explanation in my previous email.
Assuming that the coordinate of atom i and j are (xi, yi, zi) and (xj, yj, zj), respectively. Followings are the values I used for the explain:
* dx = xi - xj
* dy = yi - yj
* dz = zi - zj
* eps_i = vdW epsilon for atom i
* eps_j = vdW epsilon for atom j
* esp_ij = sqrt(eps_i * eps_j)
* r_i = vdW radii for atom i
* r_j = vdW radii for atom j
* r_ij = (r_i + r_j)
* r = sqrt(dx^2 + dy^2 + dz^2)
I used LJ12 for the van der waals potential energy, and its energy equation E is as follows:
E = eps_ij * [ (r_ij / r)^12 - 2* (r_ij / r)^6 ].
Here, I look for the 1st and 2nd derivatives of the energy by the coordinate change of atom i. Considering the relative coordinate of atom j, (xj,yj,zj) becomes the origin, and (dx,dy,dz) = (xi,yi,zi). This simplify the overall equations.
Its first derivative by xi (=dx) is
dE/dxi
= (dE / dr) * (dr / d dx)
= eps_ij * [ -12*(r_ij / r)^13 + 2*6* (r_ij / r)^7 ] * dx/sqrt(dx^2 + dy^2 + dz^2)
= eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ] * dx
= dEr * dx
The second derivative by xi and yi is
d(dE/dxi) / dxi
= d(
dEr * dx
) / d dx
= [d(dEr / d dx ] * dx
+ dEr
* [d(dx) / d dx]
= [d( eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ] * dx
) / d dx ] * dx
+ [eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ]] * 1
= [ eps_ij * [ 12*14*(r_ij / r)^15 - 2*6*8* (r_ij / r)^9 ] * dx/r ] * dx
+ (eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ])
= [ eps_ij * [ 168*(r_ij / r)^16 - 96*(r_ij / r)^10 ] * dx ] * dx
+ (eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ])
= eps_ij * [ 168*(r_ij / r)^14 - 96*(r_ij / r)^8 ] * dx ] * (dx/r)^2
+ (eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ]),
//
d(dE/dxi) / dyi
= d(
dEr * dx
) / d dy
...
= eps_ij * [ 168*(r_ij / r)^14 - 96*(r_ij / r)^8 ] * dx ] * (dx/r) * (dy/r)
This is similar to the second equation set in my previous email (There was some typos).
The above equations are the effect between only atom-i and atom-j (relative effect of atom-i by atom-j), so it is only about the off-diagonal. I agree that, for the diagonal of the hessian matrix, the summation of the 2nd derivatives of potential energies should be used.
Does my above approach correct?
Additionally, I am interested in adding my custom force-field into MMTK. It is excellent tool because it adds hydrogen to protein when it is missed in PDB file, and minimizes the conformational energy using Conjugate-Gradient and Steepest-Descent. Can I add my own force-field with giving few constraints to several few atoms without touching c source (but within only python part)? Can you give me a guide to do it?
Thank you very much.
Best regards,
-- Hyuntae
> Date: Wed, 11 Jan 2012 12:36:19 +0100
> To: htna at hotmail.com
> CC: mmtk at starship.python.net
> Subject: [MMTK] getting hessian matrix (and equation) and force vector of a protein conformaion
> From: research at khinsen.fastmail.net
>
> 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/
> ---------------------------------------------------------------------
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://starship.python.net/pipermail/mmtk/attachments/20120111/7533de88/attachment.html>
More information about the mmtk
mailing list