From 4e6604346459a4a8619bfc45f2a4b7dbd925ae86 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Thu, 30 Mar 2006 19:00:00 +0000
Subject: [PATCH] lennard jones force

---
 math/math.c |  2 +-
 math/math.h |  2 +-
 moldyn.c    | 13 +++++++++++++
 moldyn.h    |  1 +
 posic.c     |  2 ++
 5 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/math/math.c b/math/math.c
index 9217c88..6488248 100644
--- a/math/math.c
+++ b/math/math.c
@@ -76,7 +76,7 @@ double v3_norm(t_3dvec *a) {
 	return(sqrt(v3_absolute_square(a)));
 }
 
-int v3_per_bounds(t_3dvec *a,t_3dvec *dim) {
+int v3_per_bound(t_3dvec *a,t_3dvec *dim) {
 
 	double x,y,z;
 
diff --git a/math/math.h b/math/math.h
index 5ef7163..77100f3 100644
--- a/math/math.h
+++ b/math/math.h
@@ -27,6 +27,6 @@ int v3_copy(t_3dvec *trg,t_3dvec *src);
 int v3_cmp(t_3dvec *a,t_3dvec *b);
 double v3_absolute_square(t_3dvec *a);
 double v3_norm(t_3dvec *a);
-int v3_per_bounds(t_3dvec *a,t_3dvec *dim);
+int v3_per_bound(t_3dvec *a,t_3dvec *dim);
 
 #endif
diff --git a/moldyn.c b/moldyn.c
index 82749b2..57ae62f 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -226,10 +226,16 @@ int force_lennard_jones(t_moldyn *moldyn) {
 	t_3dvec distance;
 	t_3dvec force;
 	double d,h1,h2;
+	double eps,sig6,sig12;
 
 	atom=moldyn->atom;	
 	count=moldyn->count;
 	params=moldyn->pot_params;
+	eps=params->epsilon;
+	sig6=params->sigma6;
+	sig12=params->sigma12;
+
+	for(i=0;i<count;i++) v3_zero(&(atom[i].f));
 
 	for(i=0;i<count;i++) {
 		for(j=0;j<i;j++) {
@@ -242,6 +248,13 @@ int force_lennard_jones(t_moldyn *moldyn) {
 				h2=d*d;				/* 1/r^8 */
 				h1*=d;				/* 1/r^6 */
 				h1*=h2;				/* 1/r^14 */
+				h1*=sig12;
+				h2*=sig6;
+				d=-12.0*h1+6.0*h2;
+				d*=eps;
+				v3_scale(&force,&distance,d);
+				v3_add(&(atom[j].f),&(atom[j].f),&force);
+				v3_sub(&(atom[i].f),&(atom[i].f),&force);
 			}
 		}
 	}
diff --git a/moldyn.h b/moldyn.h
index b43ed8e..f515f41 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -75,5 +75,6 @@ double get_total_energy(t_moldyn *moldyn);
 t_3dvec get_total_p(t_atom *atom,int count);
 
 double potential_lennard_jones(t_moldyn *moldyn);
+int force_lennard_jones(t_moldyn *moldyn);
 
 #endif
diff --git a/posic.c b/posic.c
index 6d901b4..79654cd 100644
--- a/posic.c
+++ b/posic.c
@@ -75,6 +75,8 @@ int main(int argc,char **argv) {
 	md.count=count;
 	md.atom=si;
 	md.potential=potential_lennard_jones;
+	md.force=force_lennard_jones;
+	md.cutoff_square=((LC_SI/4.0)*(LC_SI/4.0));
 	md.pot_params=&lj;
 	md.force=NULL;
 	md.status=0;
-- 
2.20.1