<html>
<head>
</head>
<body class='hmmessage'><div dir='ltr'>


<style><!--
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Tahoma
}
--></style>
<div dir="ltr">Hi.<br><br>Sorry about my inaccurate explanation in my previous email.<br><br>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:<br>* dx = xi - xj<br>* dy = yi - yj<br>* dz = zi - zj<br>* eps_i = vdW epsilon for atom i<br>* eps_j = vdW epsilon for atom j<br>* esp_ij = sqrt(eps_i * eps_j)<br>* r_i = vdW radii for atom i<br>* r_j = vdW radii for atom j<br>* r_ij = (r_i + r_j)<br>* r = sqrt(dx^2 + dy^2 + dz^2)<br><br>I used LJ12 for the van der waals potential energy, and its energy equation E is as follows:<br>E = eps_ij * [ (r_ij / r)^12 - 2* (r_ij / r)^6 ].<br><br>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.<br><br>Its first derivative by xi (=dx) is<br>dE/dxi<br>= (dE / dr) * (dr / d dx)<br>=  eps_ij * [ -12*(r_ij / r)^13 + 2*6* (r_ij / r)^7 ]     *     dx/sqrt(dx^2 + dy^2 + dz^2)<br>=  eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ] * dx<br>= dEr * dx<br><br>The second derivative by xi and yi is<br>d(dE/dxi) / dxi<br>= d(
dEr * dx
) / d dx <br>= [d(dEr / d dx ] * dx<br>   + dEr
 * [d(dx) / d dx]<br>= [d( eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ] * dx
) / d dx ] * dx<br>    + [eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ]] * 1<br>= [ eps_ij * [ 12*14*(r_ij / r)^15 - 2*6*8* (r_ij / r)^9 ] * dx/r ] * dx<br>    + (eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ])<br>= [ eps_ij * [ 168*(r_ij / r)^16 - 96*(r_ij / r)^10 ] * dx ] * dx<br>
    + (eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ])
<br>
<BR>= eps_ij * [ 168*(r_ij / r)^14 - 96*(r_ij / r)^8 ] * dx ] * (dx/r)^2<br>    + (eps_ij * [ -12*(r_ij / r)^14 + 2*6* (r_ij / r)^8 ]),<br>
<BR>//<br>d(dE/dxi) / dyi<br>= d(
dEr * dx
) / d dy <br>...<BR>
= eps_ij * [ 168*(r_ij / r)^14 - 96*(r_ij / r)^8 ] * dx ] * (dx/r) * (dy/r)<br><br>This is similar to the second equation set in my previous email (There was some typos).<br><br>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.<br><br>Does my above approach correct?<BR><br><BR><br><BR><br><BR>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?<BR><br><BR>Thank you very much.<BR><br><BR>Best regards,<BR>-- Hyuntae<br><br><br><BR><br> 
 <br><br><div>> Date: Wed, 11 Jan 2012 12:36:19 +0100<br>> To: htna@hotmail.com<br>> CC: mmtk@starship.python.net<br>> Subject: [MMTK] getting hessian matrix (and equation) and force vector of a protein conformaion<br>> From: research@khinsen.fastmail.net<br>> <br>> Hyuntae Na writes:<br>> <br>>  > My first thought about the second derivative of the (likely say) van der waals<br>>  > (vdW) potential energy was like the follows:<br>>  > $$<br>>  >  d^2E / d dx^2 =  eps_ij * (12*13*r_ij^12/ r^12 - 2*6*7*r_ij^6/r^6) * (dx/r)^<br>>  > 2,<br>>  >  d^2E / d dx dy =  eps_ij * (12*13*r_ij^12/ r^12 - 2*6*7*r_ij^6/r^6) * (dx/r)*<br>>  > (dy/r),<br>>  > $$<br>>  > where eps_ij = sqrt(eps_i * eps_j), r_ij = vdW_rad_i + vdW_rad_j, dx (and dy)<br>>  > is the difference of x(and y)-coordinate of atoms i and j, and r=sqrt(dx^2+dy^<br>>  > 2+dz^2).<br>> <br>> I don't have the exact expressions at hand, but this looks wrong. The diagonal<br>> terms (xx, yy, zz) are sums of two contributions with different powers of r.<br>> <br>>  > However, after running a symbolic operation tool (Mathematica) to check its<br>>  > validity, the program returns a more complex equations:<br>>  > $$<br>>  >   d^2E / d dx^2 =  eps_ij * (168*r_ij^12/ r^12 - 96*r_ij^6/r^6) * (dx/r)^2   +<br>>  >   (-12*r_ij^12/r^12 + 12*r_ij^6/r^6)*(1/r)^2,<br>>  >   d^2E / d dx dy =  eps_ij * (168*r_ij^12/ r^12 - 96*r_ij^6/r^6) * (dx/r)*(dy/<br>>  > r).<br>>  > $$<br>> <br>> This looks wrong as well: the two terms in the sum for d^2E / d dx^2 have<br>> different dimensions.<br>> <br>> Note also that the elements of the Hessian are derivatives with<br>> respect to the individual coordinates x, y, z of each atom, not the<br>> relative coordinates that you call dx, dy, dz.<br>> <br>> I suggest to start by deriving the general form of the second derivatives<br>> for a generic pair potential of the form f(r). That will be a rather compact<br>> expression in terms of the first and second derivative of f(r). You can then<br>> use it for vdW, Coulomb, and bond-length interactions. Angles and Dihedrals<br>> are much messier, unfortunately.<br>> <br>>  > Additionally, would you tell me how to get the hessian matrix and forces of<br>>  > each atoms of a protein using MMTK?<br>> <br>> See the last message from this list; it contains all you need.<br>> <br>> You can also have a look at MMTK's source code. The code that adds<br>> the contribution from one pair potential to the total Hessian is at<br>> line 462 of MMTK_forcefield.c. The individiual force field terms<br>> call that code with the values of the first and second derivative<br>> of the potential function. For the vdW term, they are calculated<br>> in the code at line 483 of nonbonded.c.<br>> <br>> Good luck,<br>>   Konrad.<br>> -- <br>> ---------------------------------------------------------------------<br>> Konrad Hinsen<br>> Centre de Biophysique Moléculaire, CNRS Orléans<br>> Synchrotron Soleil - Division Expériences<br>> Saint Aubin - BP 48<br>> 91192 Gif sur Yvette Cedex, France<br>> Tel. +33-1 69 35 97 15<br>> E-Mail: research AT khinsen DOT fastmail DOT net<br>> http://dirac.cnrs-orleans.fr/~hinsen/<br>> ---------------------------------------------------------------------<br></div></div>
                                          </div></body>
</html>