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);