#include <math.h>
#include "../moldyn.h"
-#inlcude "../math/math.h"
-//#include "lennard_jones.h"
+#include "../math/math.h"
+#include "lennard_jones.h"
int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
t_3dvec force,distance;
double d,h1,h2;
double eps,sig6,sig12;
+ double energy;
- params=moldyn->pot2b_params;
+ params=moldyn->pot_params;
eps=params->epsilon4;
sig6=params->sigma6;
sig12=params->sigma12;
v3_sub(&distance,&(aj->r),&(ai->r));
if(bc) check_per_bound(moldyn,&distance);
- d=v3_absolute_square(&distance); /* 1/r^2 */
+ d=v3_absolute_square(&distance); /* r^2 */
if(d<=moldyn->cutoff_square) {
d=1.0/d; /* 1/r^2 */
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;