new src file layout (warning: doesnt compile by now!)

[physik/posic.git] / potentials / lennard_jones.c

diff --git a/potentials/lennard_jones.c b/potentials/lennard_jones.c
new file mode 100644 (file)
index 0000000..7d45a17
--- /dev/null
+++ b/potentials/lennard_jones.c
@@ -0,0 +1,60 @@
+/*
+ * lennard_jones.c - lennard jones potential
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <math.h>
+
+#include "../moldyn.h"
+#inlcude "../math/math.h"
+//#include "lennard_jones.h"
+
+int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+       t_lj_params *params;
+       t_3dvec force,distance;
+       double d,h1,h2;
+       double eps,sig6,sig12;
+
+       params=moldyn->pot2b_params;
+       eps=params->epsilon4;
+       sig6=params->sigma6;
+       sig12=params->sigma12;
+
+       if(ai<aj) return 0;
+
+       v3_sub(&distance,&(aj->r),&(ai->r));
+       if(bc) check_per_bound(moldyn,&distance);
+       d=v3_absolute_square(&distance);        /* 1/r^2 */
+       if(d<=moldyn->cutoff_square) {
+               d=1.0/d;                        /* 1/r^2 */
+               h2=d*d;                         /* 1/r^4 */
+               h1=h2*h2;                       /* 1/r^12 */
+               moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc);
+               h2*=d;                          /* 1/r^8 */
+               h1*=d;                          /* 1/r^14 */
+               h2*=6*sig6;
+               h1*=12*sig12;
+               d=+h1-h2;
+               d*=eps;
+               v3_scale(&force,&distance,d);
+               v3_add(&(aj->f),&(aj->f),&force);
+               v3_scale(&force,&force,-1.0); /* f = - grad E */
+               v3_add(&(ai->f),&(ai->f),&force);
+               virial_calc(ai,&force,&distance);
+               virial_calc(aj,&force,&distance); /* f and d signe switched */
+       }
+
+       return 0;
+}
+