[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