basic integration method and functions added

diff --git a/Makefile b/Makefile
index 854b5e6..e171354 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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