From 3961d57b84198e336085fd79263fec40837066a0 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Thu, 30 Mar 2006 23:56:11 +0000
Subject: [PATCH] basic integration method and functions added

---
 Makefile        |  4 +--
 moldyn.c        | 71 +++++++++++++++++++++++++++++++++++++++++++++++++
 moldyn.h        |  8 ++++++
 posic.c         | 16 ++++++++---
 posic.h         | 12 +--------
 random/random.h |  2 +-
 6 files changed, 95 insertions(+), 18 deletions(-)

diff --git a/Makefile b/Makefile
index 854b5e6..e171354 100644
--- a/Makefile
+++ b/Makefile
@@ -3,9 +3,9 @@ CFLAGS=-Wall
 
 OBJS=init/init.o visual/visual.o math/math.o random/random.o moldyn.o
 
-all: moldyn.o posic
+all: posic
 
-posic: moldyn.o $(OBJS)
+posic: $(OBJS) moldyn.o
 	$(CC) $(CFLAGS) -lm -o $@ $(OBJS) $(LIBS) posic.c
 
 clean:
diff --git a/moldyn.c b/moldyn.c
index 57ae62f..08901c7 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -14,6 +14,7 @@
 #include "math/math.h"
 #include "init/init.h"
 #include "random/random.h"
+#include "visual/visual.h"
 
 
 int create_lattice(unsigned char type,int element,double mass,double lc,
@@ -177,6 +178,76 @@ t_3dvec get_total_p(t_atom *atom, int count) {
 }
 
 
+/*
+ *
+ * 'integration of newtons equation' - algorithms
+ *
+ */
+
+/* start the integration */
+
+int moldyn_integrate(t_moldyn *moldyn) {
+
+	int i;
+
+	/* calculate initial forces */
+	moldyn->force(moldyn);
+
+	for(i=0;i<moldyn->time_steps;i++) {
+		/* integration step */
+		moldyn->integrate(moldyn);
+
+		/* check for visualiziation */
+		// to be continued ...
+		if(!(i%100)) 
+			visual_atoms(moldyn->visual,i*moldyn->tau,
+			             moldyn->atom,moldyn->count);
+	}
+
+	return 0;
+}
+
+/* velocity verlet */
+
+int velocity_verlet(t_moldyn *moldyn) {
+
+	int i,count;
+	double tau,tau_square;
+	t_3dvec delta;
+	t_atom *atom;
+
+	atom=moldyn->atom;
+	count=moldyn->count;
+	tau=moldyn->tau;
+
+	tau_square=tau*tau;
+
+	for(i=0;i<count;i++) {
+		/* new positions */
+		v3_scale(&delta,&(atom[i].v),tau);
+		v3_add(&(atom[i].r),&(atom[i].r),&delta);
+		v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
+		v3_add(&(atom[i].r),&(atom[i].r),&delta);
+		v3_per_bound(&(atom[i].r),&(moldyn->dim));
+
+		/* velocities */
+		v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
+		v3_add(&(atom[i].r),&(atom[i].r),&delta);
+	}
+
+	/* forces depending on chosen potential */
+	moldyn->force(moldyn);
+
+	for(i=0;i<count;i++) {
+		/* again velocities */
+		v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
+		v3_add(&(atom[i].v),&(atom[i].v),&delta);
+	}
+
+	return 0;
+}
+
+
 /*
  *
  * potentials & corresponding forces
diff --git a/moldyn.h b/moldyn.h
index f515f41..d94e193 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -10,6 +10,7 @@
 
 #include "math/math.h"
 #include "random/random.h"
+//#include "visual/visual.h"
 
 /* datatypes */
 
@@ -29,6 +30,10 @@ typedef struct s_moldyn {
 	int (*force)(struct s_moldyn *moldyn);
 	double cutoff_square;
 	t_3dvec dim;
+	int (*integrate)(struct s_moldyn *moldyn);
+	int time_steps;
+	double tau;
+	void *visual;
 	unsigned char status;
 } t_moldyn;
 
@@ -74,6 +79,9 @@ double get_e_pot(t_moldyn *moldyn);
 double get_total_energy(t_moldyn *moldyn);
 t_3dvec get_total_p(t_atom *atom,int count);
 
+int moldyn_integrate(t_moldyn *moldyn);
+int velocity_verlet(t_moldyn *moldyn);
+
 double potential_lennard_jones(t_moldyn *moldyn);
 int force_lennard_jones(t_moldyn *moldyn);
 
diff --git a/posic.c b/posic.c
index 79654cd..04d9b7f 100644
--- a/posic.c
+++ b/posic.c
@@ -57,9 +57,6 @@ int main(int argc,char **argv) {
 	printf("setting thermal fluctuations\n");
 	thermal_init(si,&random,count,t);
 
-	/* visualize */
-
-	visual_atoms(&vis,0.0,si,count);
 
 	/* check kinetic energy */
 
@@ -78,8 +75,11 @@ int main(int argc,char **argv) {
 	md.force=force_lennard_jones;
 	md.cutoff_square=((LC_SI/4.0)*(LC_SI/4.0));
 	md.pot_params=&lj;
-	md.force=NULL;
+	md.integrate=velocity_verlet;
+	md.time_steps=RUNS;
+	md.tau=TAU;
 	md.status=0;
+	md.visual=&vis;
 
 	lj.sigma6=3.0/16.0*LC_SI*LC_SI;
 	help=lj.sigma6*lj.sigma6;
@@ -103,6 +103,14 @@ int main(int argc,char **argv) {
 	 * integration of newtons equations
 	 */
 
+	/* visualize */
+	//visual_atoms(&vis,0.0,si,count);
+	
+
+	moldyn_integrate(&md);
+
+	printf("total energy (after integration): %f\n",get_total_energy(&md));
+
 	/* close */
 
 	visual_tini(&vis);
diff --git a/posic.h b/posic.h
index 8f306ec..854dfb8 100644
--- a/posic.h
+++ b/posic.h
@@ -18,24 +18,14 @@
 #define POSIC_H
 
 #define RUNS 15000
-#define TAU 0.001
+#define TAU 0.000001
 
 #define TEMPERATURE 273.0 
 
-#define SI_M 1
-#define SI_LC 5.43105
-#define LJ_SIGMA SI_LC
-#define LJ_SIGMA_02 (LJ_SIGMA*LJ_SIGMA)
-#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 1
 #define LEN_Y 1
 #define LEN_Z 1
 
 #define R_CUTOFF 20
-#define R2_CUTOFF (R_CUTOFF*R_CUTOFF)
-
-#define AMOUNT_SI ((LEN_X/SI_LC)*(LEN_Y/SI_LC)*(LEN_Z/SI_LC)*2)
 
 #endif
diff --git a/random/random.h b/random/random.h
index 35d5706..5c0d7b7 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) /* 4 k byte */
+#define RAND_BUFSIZE 			(1024*1024) /* 4 mega byte */
 
 #define URAND_MAX			0xffffffff
 #define URAND_MAX_PLUS_ONE		0x100000000
-- 
2.20.1