t_3dvec force,distance;
double d,h1,h2;
double eps,sig6,sig12;
+ double energy;
params=moldyn->pot_params;
eps=params->epsilon4;
h2=d*d; /* 1/r^4 */
h2*=d; /* 1/r^6 */
h1=h2*h2; /* 1/r^12 */
- moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc);
+ energy=(eps*(sig12*h1-sig6*h2)-params->uc);
+ moldyn->energy+=energy; /* total energy */
+ ai->e+=0.5*energy; /* site energy */
+ aj->e+=0.5*energy;
h2*=d; /* 1/r^8 */
h1*=d; /* 1/r^14 */
h2*=6*sig6;