From: hackbard <hackbard>
Date: Thu, 30 Mar 2006 13:49:53 +0000 (+0000)
Subject: started potential and force stuff
X-Git-Url: https://www.hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=commitdiff_plain;h=85abe46fecc79a3d7885ff6124186b2be2ca96ac

started potential and force stuff
---

diff --git a/math/math.c b/math/math.c
index 9cfaff6..9217c88 100644
--- a/math/math.c
+++ b/math/math.c
@@ -76,3 +76,20 @@ double v3_norm(t_3dvec *a) {
 	return(sqrt(v3_absolute_square(a)));
 }
 
+int v3_per_bounds(t_3dvec *a,t_3dvec *dim) {
+
+	double x,y,z;
+
+	x=0.5*dim->x;
+	y=0.5*dim->y;
+	z=0.5*dim->z;
+
+	if(a->x>x) a->x-=dim->x;
+	else if(-a->x>x) a->x+=dim->x;
+	if(a->y>y) a->y-=dim->y;
+	else if(-a->y>y) a->y+=dim->y;
+	if(a->z>z) a->z-=dim->z;
+	else if(-a->z>z) a->z+=dim->z;
+
+	return 0;
+}
diff --git a/math/math.h b/math/math.h
index 53e2114..5ef7163 100644
--- a/math/math.h
+++ b/math/math.h
@@ -27,5 +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);
 
 #endif
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;
+}
 
diff --git a/moldyn.h b/moldyn.h
index 30c4bf8..b43ed8e 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -21,12 +21,14 @@ typedef struct s_atom {
 	double mass;	/* atom mass */
 } t_atom;
 
-
 typedef struct s_moldyn {
 	int count;
 	t_atom *atom;
-	double (*potential)(void *ptr);
-	int (*force)(struct s_moldyn *moldyn,void *ptr);
+	double (*potential)(struct s_moldyn *moldyn);
+	void *pot_params;
+	int (*force)(struct s_moldyn *moldyn);
+	double cutoff_square;
+	t_3dvec dim;
 	unsigned char status;
 } t_moldyn;
 
@@ -68,8 +70,10 @@ int destroy_lattice(t_atom *atom);
 int thermal_init(t_atom *atom,t_random *random,int count,double t);
 int scale_velocity(t_atom *atom,int count,double t);
 double get_e_kin(t_atom *atom,int count);
+double get_e_pot(t_moldyn *moldyn);
+double get_total_energy(t_moldyn *moldyn);
 t_3dvec get_total_p(t_atom *atom,int count);
 
-double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr);
+double potential_lennard_jones(t_moldyn *moldyn);
 
 #endif
diff --git a/posic.c b/posic.c
index 4eb2d04..6d901b4 100644
--- a/posic.c
+++ b/posic.c
@@ -14,6 +14,8 @@
 
 int main(int argc,char **argv) {
 
+	t_moldyn md;
+
 	t_atom *si;
 
 	t_visual vis;
@@ -21,10 +23,13 @@ int main(int argc,char **argv) {
 	t_random random;
 
 	int a,b,c;
-	double t,e;
+	double t,e,u;
+	double help;
 	t_3dvec p;
 	int count;
 
+	t_lj_params lj;
+
 	char fb[32]="saves/fcc_test";
 
 	/* init */
@@ -66,6 +71,30 @@ int main(int argc,char **argv) {
 	p=get_total_p(si,count);
 	printf("total momentum: %f\n",v3_norm(&p));
 
+	/* check potential energy */
+	md.count=count;
+	md.atom=si;
+	md.potential=potential_lennard_jones;
+	md.pot_params=&lj;
+	md.force=NULL;
+	md.status=0;
+
+	lj.sigma6=3.0/16.0*LC_SI*LC_SI;
+	help=lj.sigma6*lj.sigma6;
+	lj.sigma6*=help;
+	lj.sigma12=lj.sigma6*lj.sigma6;
+	lj.epsilon=1;
+
+	u=get_e_pot(&md);
+
+	printf("potential energy: %f\n",u);
+	printf("total energy (1): %f\n",e+u);
+	printf("total energy (2): %f\n",get_total_energy(&md));
+
+	md.dim.x=a*LC_SI;
+	md.dim.y=b*LC_SI;
+	md.dim.z=c*LC_SI;
+
 	/*
 	 * let's do the actual md algorithm now
 	 *