X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=82749b2f10abf978c9636e213585e21c8393741a;hb=85abe46fecc79a3d7885ff6124186b2be2ca96ac;hp=ecba8baffc40a07c641eee9bfed8838c7c8116cc;hpb=eceebe3ee412aa8cea3e6a7f0038883707f78460;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index ecba8ba..82749b2 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -147,6 +147,21 @@ double get_e_kin(t_atom *atom,int count) {
 	return e;
 }
 
+double get_e_pot(t_moldyn *moldyn) {
+
+	return(moldyn->potential(moldyn));
+}
+
+double get_total_energy(t_moldyn *moldyn) {
+
+	double e;
+
+	e=get_e_kin(moldyn->atom,moldyn->count);
+	e+=get_e_pot(moldyn);
+
+	return e;
+}
+
 t_3dvec get_total_p(t_atom *atom, int count) {
 
 	t_3dvec p,p_total;
@@ -161,7 +176,16 @@ t_3dvec get_total_p(t_atom *atom, int count) {
 	return p_total;
 }
 
-double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr) {
+
+/*
+ *
+ * potentials & corresponding forces
+ * 
+ */
+
+/* lennard jones potential & force for one sort of atoms */
+ 
+double potential_lennard_jones(t_moldyn *moldyn) {
 
 	t_lj_params *params;
 	t_atom *atom;
@@ -172,7 +196,7 @@ double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr) {
 	double u;
 	double eps,sig6,sig12;
 
-	params=ptr;
+	params=moldyn->pot_params;
 	atom=moldyn->atom;
 	count=moldyn->count;
 	eps=params->epsilon;
@@ -183,7 +207,7 @@ double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr) {
 	for(i=0;i<count;i++) {
 		for(j=0;j<i;j++) {
 			v3_sub(&distance,&(atom[j].r),&(atom[i].r));
-			d=v3_absolute_square(&distance); 	/* 1/r^2 */
+			d=1.0/v3_absolute_square(&distance); 	/* 1/r^2 */
 			help=d*d; 				/* 1/r^4 */
 			help*=d; 				/* 1/r^6 */
 			d=help*help; 				/* 1/r^12 */
@@ -194,4 +218,34 @@ double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr) {
 	return u;
 }
 
+int force_lennard_jones(t_moldyn *moldyn) {
+
+	t_lj_params *params;
+	int i,j,count;
+	t_atom *atom;
+	t_3dvec distance;
+	t_3dvec force;
+	double d,h1,h2;
+
+	atom=moldyn->atom;	
+	count=moldyn->count;
+	params=moldyn->pot_params;
+
+	for(i=0;i<count;i++) {
+		for(j=0;j<i;j++) {
+			v3_sub(&distance,&(atom[j].r),&(atom[i].r));
+			v3_per_bound(&distance,&(moldyn->dim));
+			d=v3_absolute_square(&distance);
+			if(d<=moldyn->cutoff_square) {
+				h1=1.0/d;			/* 1/r^2 */
+				d=h1*h1;			/* 1/r^4 */
+				h2=d*d;				/* 1/r^8 */
+				h1*=d;				/* 1/r^6 */
+				h1*=h2;				/* 1/r^14 */
+			}
+		}
+	}
+
+	return 0;
+}