From eceebe3ee412aa8cea3e6a7f0038883707f78460 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Wed, 29 Mar 2006 16:44:28 +0000
Subject: [PATCH] started potentials

---
 Makefile        |  2 +-
 moldyn.c        | 76 ++++++++++++++++++++++++++++++++++++++++---------
 moldyn.h        | 31 ++++++++++++++++++--
 posic.c         | 20 +++++++++++--
 posic.h         |  8 +++---
 random/random.c | 15 +++++++---
 random/random.h |  2 +-
 7 files changed, 127 insertions(+), 27 deletions(-)

diff --git a/Makefile b/Makefile
index ff3fbe3..854b5e6 100644
--- a/Makefile
+++ b/Makefile
@@ -9,4 +9,4 @@ posic: moldyn.o $(OBJS)
 	$(CC) $(CFLAGS) -lm -o $@ $(OBJS) $(LIBS) posic.c
 
 clean:
-	rm -f *.o posic
+	rm -f *.o posic */*.o
diff --git a/moldyn.c b/moldyn.c
index b507d65..ecba8ba 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -87,25 +87,26 @@ int thermal_init(t_atom *atom,t_random *random,int count,double t) {
 	/* gaussian distribution of velocities */
 	v3_zero(&p_total);
 	for(i=0;i<count;i++) {
-		sigma=sqrt(2.0*K_BOLTZMANN*t/atom[count].mass);
+		sigma=sqrt(2.0*K_BOLTZMANN*t/atom[i].mass);
 		/* x direction */
 		v=sigma*rand_get_gauss(random);
-		atom[count].v.x=v;
-		p_total.x+=atom[count].mass*v;
+		atom[i].v.x=v;
+		p_total.x+=atom[i].mass*v;
 		/* y direction */
 		v=sigma*rand_get_gauss(random);
-		atom[count].v.y=v;
-		p_total.x+=atom[count].mass*v;
+		atom[i].v.y=v;
+		p_total.y+=atom[i].mass*v;
 		/* z direction */
 		v=sigma*rand_get_gauss(random);
-		atom[count].v.z=v;
-		p_total.x+=atom[count].mass*v;
+		atom[i].v.z=v;
+		p_total.z+=atom[i].mass*v;
 	}
 
 	/* zero total momentum */
+	v3_scale(&p_total,&p_total,1.0/count);
 	for(i=0;i<count;i++) {
-		v3_scale(&delta,&p_total,1.0/atom[count].mass);
-		v3_sub(&(atom[count].v),&(atom[count].v),&delta);
+		v3_scale(&delta,&p_total,1.0/atom[i].mass);
+		v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 	}
 
 	/* velocity scaling */
@@ -124,10 +125,10 @@ int scale_velocity(t_atom *atom,int count,double t) {
 	 */
 	e=0.0;
 	for(i=0;i<count;i++)
-		e+=0.5*atom[count].mass*v3_absolute_square(&(atom[count].v));
+		e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 	c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*t));
 	for(i=0;i<count;i++)
-		v3_scale(&(atom[count].v),&(atom[count].v),(1.0/c));
+		v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
 
 	return 0;
 }
@@ -139,9 +140,58 @@ double get_e_kin(t_atom *atom,int count) {
 
 	e=0.0;
 
-	for(i=0;i<count;i++)
-		e+=0.5*atom[count].mass*v3_absolute_square(&(atom[count].v));
+	for(i=0;i<count;i++) {
+		e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
+	}
 
 	return e;
 }
 
+t_3dvec get_total_p(t_atom *atom, int count) {
+
+	t_3dvec p,p_total;
+	int i;
+
+	v3_zero(&p_total);
+	for(i=0;i<count;i++) {
+		v3_scale(&p,&(atom[i].v),atom[i].mass);
+		v3_add(&p_total,&p_total,&p);
+	}
+
+	return p_total;
+}
+
+double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr) {
+
+	t_lj_params *params;
+	t_atom *atom;
+	int i,j;
+	int count;
+	t_3dvec distance;
+	double d,help;
+	double u;
+	double eps,sig6,sig12;
+
+	params=ptr;
+	atom=moldyn->atom;
+	count=moldyn->count;
+	eps=params->epsilon;
+	sig6=params->sigma6;
+	sig12=params->sigma12;
+
+	u=0.0;
+	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 */
+			help=d*d; 				/* 1/r^4 */
+			help*=d; 				/* 1/r^6 */
+			d=help*help; 				/* 1/r^12 */
+			u+=eps*(sig12*d-sig6*help);
+		}
+	}
+	
+	return u;
+}
+
+
diff --git a/moldyn.h b/moldyn.h
index 65f0d8d..30c4bf8 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -22,9 +22,33 @@ typedef struct s_atom {
 } t_atom;
 
 
-/* defines */
+typedef struct s_moldyn {
+	int count;
+	t_atom *atom;
+	double (*potential)(void *ptr);
+	int (*force)(struct s_moldyn *moldyn,void *ptr);
+	unsigned char status;
+} t_moldyn;
 
-#define K_BOLTZMANN		1.3807E-23
+typedef struct s_lj_params {
+	double sigma6;
+	double sigma12;
+	double epsilon;
+} t_lj_params;
+
+/*
+ *  defines
+ */
+
+/* general defines */
+
+#define MOLDYN_STAT_POTENTIAL		0x01
+#define MOLDYN_STAT_FORCE		0x02
+
+/* phsical values */
+
+//#define K_BOLTZMANN		1.3807E-23
+#define K_BOLTZMANN		1.0
 
 #define FCC			0x01
 #define DIAMOND			0x02
@@ -44,5 +68,8 @@ 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);
+t_3dvec get_total_p(t_atom *atom,int count);
+
+double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr);
 
 #endif
diff --git a/posic.c b/posic.c
index d7a659a..4eb2d04 100644
--- a/posic.c
+++ b/posic.c
@@ -22,6 +22,7 @@ int main(int argc,char **argv) {
 
 	int a,b,c;
 	double t,e;
+	t_3dvec p;
 	int count;
 
 	char fb[32]="saves/fcc_test";
@@ -31,6 +32,11 @@ int main(int argc,char **argv) {
 	rand_init(&random,NULL,1);
 	random.status|=RAND_STAT_VERBOSE;
 
+	/* testing random numbers */
+	//for(a=0;a<1000000;a++)
+	//	printf("%f %f\n",rand_get_gauss(&random),
+	//	                 rand_get_gauss(&random));
+
 	visual_init(&vis,fb);
 
 	a=LEN_X;
@@ -39,13 +45,13 @@ int main(int argc,char **argv) {
 
 	t=TEMPERATURE;
 
-	printf("placing silicon atoms\n");
+	printf("placing silicon atoms ... ");
 	count=create_lattice(DIAMOND,Si,M_SI,LC_SI,a,b,c,&si);
+	printf("(%d) ok!\n",count);
 
 	printf("setting thermal fluctuations\n");
 	thermal_init(si,&random,count,t);
 
-
 	/* visualize */
 
 	visual_atoms(&vis,0.0,si,count);
@@ -56,6 +62,16 @@ int main(int argc,char **argv) {
 	printf("kinetic energy: %f\n",e);
 	printf("3/2 N k T = %f\n",1.5*count*K_BOLTZMANN*t);
 
+	/* check total momentum */
+	p=get_total_p(si,count);
+	printf("total momentum: %f\n",v3_norm(&p));
+
+	/*
+	 * let's do the actual md algorithm now
+	 *
+	 * integration of newtons equations
+	 */
+
 	/* close */
 
 	visual_tini(&vis);
diff --git a/posic.h b/posic.h
index de0106d..8f306ec 100644
--- a/posic.h
+++ b/posic.h
@@ -20,7 +20,7 @@
 #define RUNS 15000
 #define TAU 0.001
 
-#define TEMPERATURE 0.0
+#define TEMPERATURE 273.0 
 
 #define SI_M 1
 #define SI_LC 5.43105
@@ -29,9 +29,9 @@
 #define LJ_SIGMA_06 (LJ_SIGMA_02*LJ_SIGMA_02*LJ_SIGMA_02)
 #define LJ_SIGMA_12 (LJ_SIGMA_06*LJ_SIGMA_06)
 
-#define LEN_X 2
-#define LEN_Y 2
-#define LEN_Z 2
+#define LEN_X 1
+#define LEN_Y 1
+#define LEN_Z 1
 
 #define R_CUTOFF 20
 #define R2_CUTOFF (R_CUTOFF*R_CUTOFF)
diff --git a/random/random.c b/random/random.c
index c06b98a..ddbd7b5 100644
--- a/random/random.c
+++ b/random/random.c
@@ -62,21 +62,28 @@ int rand_close(t_random *random) {
 
 unsigned int rand_get(t_random *random) {
 
+	int left;
+
 	if(random->b_ptr==random->buffer+RAND_BUFSIZE) {
 		if(random->status&RAND_STAT_VERBOSE)
 			dprintf(random->logfd,
 			        "[random] getting new random numbers\n");
-		random->b_ptr=random->buffer;
 		if(!(random->status&RAND_STAT_UDEV)) {
 			lseek(random->fd,0,SEEK_SET);
 			dprintf(random->logfd,
 			        "[random] warning, rereading random file\n");
 		}
-		read(random->fd,random->b_ptr,
-		     RAND_BUFSIZE*sizeof(unsigned int));
+		left=RAND_BUFSIZE*sizeof(unsigned int);
+		while(left) {
+			left-=read(random->fd,
+			           random->buffer+\
+			           RAND_BUFSIZE*sizeof(unsigned int)-left,
+			           left);
+		}
 		if(random->status&RAND_STAT_VERBOSE)
 			dprintf(random->logfd,
 			        "[random] got new random numbers\n");
+		random->b_ptr=random->buffer;
 	}
 
 	return(*(random->b_ptr++));
@@ -84,7 +91,7 @@ unsigned int rand_get(t_random *random) {
 
 double rand_get_double(t_random *random) {
 
-	return(1.0*rand_get(random)/(double)(URAND_MAX+1));
+	return(1.0*rand_get(random)/((long long unsigned int)URAND_MAX+1));
 }
 
 double rand_get_double_linear(t_random *random,double a,double b) {
diff --git a/random/random.h b/random/random.h
index 8a21276..35d5706 100644
--- a/random/random.h
+++ b/random/random.h
@@ -26,7 +26,7 @@ typedef struct s_random {
 #define RAND_STAT_UDEV			0x04
 #define RAND_STAT_GAUSS			0x08
 
-#define RAND_BUFSIZE 			(1024) /* 16 mega byte */
+#define RAND_BUFSIZE 			(1024) /* 4 k byte */
 
 #define URAND_MAX			0xffffffff
 #define URAND_MAX_PLUS_ONE		0x100000000
-- 
2.20.1