Energy calculation issue in dihedarl_charmm routine

Hi Steve,

While I was writing my own code and validating with lammps, I came across issue with dihedral_charmm.cpp. May be you can have a look into it.

I have shown a sample calculations for 4 atoms whose co-ordinates are given below. The discrepancy comes in sin theta values (-ve sign) and hence, leads to different energy of 36.069, whereas it should be 45.75.

I have also checked that you are using the same approach for dihedral_harmonics as well. This should create problem there as well.

Best,
Ateeque

I don’t know of any known bug/issue with the dihedral charmm implementation.

It has been checked (by others) against CHARMM itself. Unless the

geometry you are testing is an oddball case (e.g. extreme angle), then

i doubt that LAMMPS is doing it wrong. Others with more experience

with the CHARMM ff can comment.

Steve

Hi Steve,

Hi Steve,

No it is not an oddball configuration. The angle is somewhere around 68
deg.

​at any rate, the benchmark here is the CHARMM code. if you can reliably
prove that LAMMPS differs from CHARMM​, then somebody will likely feel
compelled to look into it. if it differs from your newly written code,
most likely not.

are you certain that in your code is following the same conventions for the
ordering of atoms in the definition of the dihedrals is consistent with the
force/energy computation? that has been the most common reason of
discrepancies in the past.

axel.

Hi Axel,

Thanks for your comment.

The values I have printed is from LAMMPS, not my code.

Here is an exhaustive print-out from LAMMPS dihedral_charmm.cpp. I have just added few fprintf(lines in the code). Please see that for some atoms, the energy calculated directly K*(1+cos(n*theta-d)) is same as calculated using LAMMPS. However, for the example I have taken its different.

The calculations are correct till cos (theta).

I believe the problem is with sin calculation. Once error comes in sin, then it gets amplified during iterative calculations of cos(n*theta).

The printing commands I have added to dihedral_charmm.cpp are given below.

Best,
Ateeque

21-atik
15 563 566 565 573
   1st atom: 563 4.268 -9.135 2.031
   2nd atom: 566 5.333 -10.041 2.683
   3rd atom: 565 5.199 -11.51 2.307
   4th atom: 573 5.169 -11.637 1.22
   1st bond: vb1x,vb1y,vb1z -1.065000 0.906000 -0.652000
   2nd bond: vb2x,vb2y,vb2z -0.134000 -1.469000 -0.376000
   2nd bond: vb1xm,vb1ym,vb2zm 0.134000 1.469000 0.376000
   3rd bond: vb3x,vb3y,vb3z -0.030000 -0.127000 -1.087000
   ax ay az 1.298444 0.313072 -1.685889
   bx by bz 1.549051 -0.134378 -0.027052
   rasq rbsq rgsq rg 4.626193 2.418348 2.317293 1.522266
   rginv ra2inv rb2inv rabinv 0.656916 0.216160 0.413505 0.298971
   c-cos s-sin acos(c) and sin(acos(c)) 0.602394 0.798199 0.924299 0.798199
   energy parameters k,n,d(rad) 30.000000 1 0.174533
   energy calculation directly using formulate k*(1+cos(n*theta-d)) 51.955442
   energy calculation from LAMMPS 51.955442
22-atik
15 563 566 565 574
   1st atom: 563 4.268 -9.135 2.031
   2nd atom: 566 5.333 -10.041 2.683
   3rd atom: 565 5.199 -11.51 2.307
   4th atom: 574 4.295 -11.953 2.735
   1st bond: vb1x,vb1y,vb1z -1.065000 0.906000 -0.652000
   2nd bond: vb2x,vb2y,vb2z -0.134000 -1.469000 -0.376000
   2nd bond: vb1xm,vb1ym,vb2zm 0.134000 1.469000 0.376000
   3rd bond: vb3x,vb3y,vb3z -0.904000 -0.443000 0.428000
   ax ay az 1.298444 0.313072 -1.685889
   bx by bz -0.795300 0.397256 -1.268614
   rasq rbsq rgsq rg 4.626193 2.399696 2.317293 1.522266
   rginv ra2inv rb2inv rabinv 0.656916 0.216160 0.416719 0.300130
   c-cos s-sin acos(c) and sin(acos(c)) 0.369298 -0.929311 1.192542 0.929311
   energy parameters k,n,d(rad) 30.000000 1 0.174533
   energy calculation directly using formulate k*(1+cos(n*theta-d)) 45.751831
   energy calculation from LAMMPS 36.069443

      fprintf(screen,"%d-atik\n",n);
      fprintf(screen,"%d %d %d %d %d\n",type,i1,i2,i3,i4);

         fprintf(screen," 1st atom: %d %g %g %g\n",
                 i1,x[i1][0],x[i1][1],x[i1][2]);
         fprintf(screen," 2nd atom: %d %g %g %g\n",
                 i2,x[i2][0],x[i2][1],x[i2][2]);
         fprintf(screen," 3rd atom: %d %g %g %g\n",
                 i3,x[i3][0],x[i3][1],x[i3][2]);
         fprintf(screen," 4th atom: %d %g %g %g\n",
                 i4,x[i4][0],x[i4][1],x[i4][2]);

         fprintf(screen," 1st bond: vb1x,vb1y,vb1z %f %f %f\n",
                 vb1x,vb1y,vb1z);
         fprintf(screen," 2nd bond: vb2x,vb2y,vb2z %f %f %f\n",
                 vb2x,vb2y,vb2z);
         fprintf(screen," 2nd bond: vb1xm,vb1ym,vb2zm %f %f %f\n",
                 vb2xm,vb2ym,vb2zm);
         fprintf(screen," 3rd bond: vb3x,vb3y,vb3z %f %f %f\n",
                 vb3x,vb3y,vb3z);
         fprintf(screen," ax ay az %f %f %f\n",
                 ax,ay,az);
         fprintf(screen," bx by bz %f %f %f\n",
                 bx,by,bz);
         fprintf(screen," rasq rbsq rgsq rg %f %f %f %f\n",
                 rasq,rbsq,rgsq,rg);
         fprintf(screen," rginv ra2inv rb2inv rabinv %f %f %f %f\n",
                 rginv,ra2inv,rb2inv,rabinv);

         fprintf(screen," c-cos s-sin acos(c) and sin(acos(c)) %f %f %f %f \n",
                    c,s,acos(c),sin(acos(c)));

     fprintf(screen," energy parameters k,n,d(rad) %f %d %f \n",
                 k[type],m,acos(cos_shift[type]) );
     fprintf(screen," energy calculation directly using formulate k*(1+cos(n*theta-d)) %f\n",
                 k[type]*( 1+cos( acos(c)*m-acos(cos_shift[type]) ) ) );
     fprintf(screen," energy calculation from LAMMPS %f\n",k[type]*p);

Dear Axel and Steve,

Please see my previous email, where I have printed various variables by including few lines in dihedral_charmm.cpp file of LAMMPS. I have also printed the energy K*(1+cos(n*theta-d)) in lammps dihedral_chamm.cpp file itself.

I show the data for 2 dihedral case. Interestingly, for 1-case the calculated value is same as calculated using direct formula. However, for the second case, the two answers are different, and the reason is -ve sign of sin(theta), which is further used iteratively to calculate cos(n*theta). You can check the calculations simply as well.

Best,
Ateeque

Hello Ateeque,

I’m not very familiar with CHARMM, but isn’t the dihedral angle defined over the whole circle? In this case the angle is signed and only its magnitude is given by acos(cos(theta)). Assuming the angle is defined in the plane normal to the vector vb2, the sign of dot(vb2,cross(a,b)) would determine the sign of the angle.

-David

Hi David,

Thank you for your comment. That does make sense.

Apologies to Steve and Axel.

Best,
Ateeque