some bugfixes ... com'on

diff --git a/Makefile b/Makefile
index b037404..dabf04c 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,8 +1,8 @@
 CC=gcc
-CFLAGS=-Wall #-DSIMPLE_TESTING
+CFLAGS=-Wall -O3 #-DSIMPLE_TESTING
 LDFLAGS=-lm
 
-OBJS=init/init.o visual/visual.o math/math.o random/random.o list/list.o
+OBJS=init/init.o visual/visual.o random/random.o list/list.o
 OBJS+=moldyn.o
 
 all: sic

diff --git a/moldyn.c b/moldyn.c
index 716a087..2631d6a 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -873,15 +873,7 @@ int potential_force_calc(t_moldyn *moldyn) {
                                           !(jtom->attr&ATOM_ATTR_3BP))
                                                continue;
 
-                       /*
-                        * according to mr. nordlund, we dont need to take the 
-                        * sum over all atoms now, as 'this is centered' around
-                        * atom i ...
-                        * i am not quite sure though! there is a not vanishing
-                        * part even if f_c_ik is zero ...
-                        * this analytical potentials suck!
-                        * switching from mc to md to dft soon!
-                        */
+                       /* neighbourhood of atom j is not needed! */
 
                        //              link_cell_neighbour_index(moldyn,
                        //                 (jtom->r.x+moldyn->dim.x/2)/lc->x,
@@ -1156,6 +1148,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        exchange=&(params->exchange);
 
        exchange->run3bp=0;
+       exchange->run2bp_post=0;
        
        /*
         * we need: f_c, df_c, f_r, df_r
@@ -1200,7 +1193,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        f_r=A*exp(-lambda*d_ij);
        df_r=-lambda*f_r/d_ij;
 
-       /* f_a, df_a calc + save for 3bp use */
+       /* f_a, df_a calc + save for later use */
        exchange->f_a=-B*exp(-mu*d_ij);
        exchange->df_a=-mu*exchange->f_a/d_ij;
 
@@ -1228,8 +1221,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        exchange->f_c=f_c;
        exchange->df_c=df_c;
 
-       /* enable the run of 3bp function */
+       /* enable the run of 3bp function and 2bp post processing */
        exchange->run3bp=1;
+       exchange->run2bp_post=1;
 
        /* reset 3bp sums */
        exchange->sum1_3bp=0.0;
@@ -1259,6 +1253,10 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        params=moldyn->pot2b_params;
        exchange=&(params->exchange);
 
+       /* we do not run if f_c_ij was dtected to be 0! */
+       if(!(exchange->run2bp_post))
+               return 0;
+
        db_ij=&(exchange->db_ij);
        f_c=exchange->f_c;
        df_c=exchange->df_c;
@@ -1266,6 +1264,7 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        df_a=exchange->df_a;
        betan=exchange->betan;
        n=*(exchange->n);
+       chi=exchange->chi;
        dist_ij=&(exchange->dist_ij);
 
        db_ij_scale1=(1+betan*exchange->sum1_3bp);
@@ -1274,19 +1273,24 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        b_ij=chi*db_ij_scale1*help;
        db_ij_scale1=-chi/(2*n)*help;
 
+       /* db_ij part */
        v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2));
        v3_scale(db_ij,db_ij,f_a);
 
+       /* df_a part */
        v3_scale(&temp,dist_ij,b_ij*df_a);
 
+       /* db_ij + df_a part */
        v3_add(&force,&temp,db_ij);
        v3_scale(&force,&force,f_c);
 
+       /* df_c part */
        v3_scale(&temp,dist_ij,f_a*b_ij*df_c);
 
        /* add energy of 3bp sum */
        moldyn->energy+=(0.5*f_c*b_ij*f_a);
-       /* add force of 3bp calculation */      
+
+       /* add force of 3bp calculation (all three parts) */
        v3_add(&(ai->f),&temp,&force);
 
        return 0;
@@ -1415,13 +1419,13 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                bracket=0.0;
                bracket_n_1=0.0;
                bracket_n=0.0;
-               printf("Foo -> 0: ");
+               //printf("Foo -> 0: ");
        }
        else {
                bracket=f_c_ik*g;
                bracket_n_1=pow(bracket,n-1.0);
                bracket_n=bracket_n_1*bracket;
-               printf("Foo -> 1: ");
+               //printf("Foo -> 1: ");
        }
 //printf("%.15f %.15f %.15f\n",bracket_n_1,bracket_n,bracket);
 

diff --git a/moldyn.h b/moldyn.h
index 09471a1..1145f8a 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -197,6 +197,7 @@ typedef struct s_tersoff_exchange {
        double n_betan;
 
        u8 run3bp;
+       u8 run2bp_post;
 
        t_3dvec db_ij;
        double sum1_3bp;

diff --git a/sic.c b/sic.c
index 5f50bd8..26bb8c8 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -96,7 +96,7 @@ int main(int argc,char **argv) {
 
        /* set (initial) dimensions of simulation volume */
        printf("[sic] setting dimensions\n");
-       set_dim(&md,3*LC_SI,3*LC_SI,3*LC_SI,TRUE);
+       set_dim(&md,5*LC_SI,5*LC_SI,5*LC_SI,TRUE);
 
        /* set periodic boundary conditions in all directions */
        printf("[sic] setting periodic boundary conditions\n");
@@ -106,15 +106,17 @@ int main(int argc,char **argv) {
        printf("[sic] creating atoms\n");
        create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+                      //ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP,
                       //ATOM_ATTR_2BP|ATOM_ATTR_HB,
-                      0,3,3,3);
+                      0,5,5,5);
 
        /* setting a nearest neighbour distance for the moldyn checks */
        set_nn_dist(&md,sqrt(3.0)*LC_SI/4.0); /* diamond ! */
 
        /* set temperature */
        printf("[sic] setting temperature\n");
-       set_temperature(&md,0.0);
+       set_temperature(&md,273.0+450.0);
+       //set_temperature(&md,0.0);
 
        /* set p/t scaling */
        printf("[sic] set p/t scaling\n");
@@ -126,12 +128,12 @@ int main(int argc,char **argv) {
 
        /* create the simulation schedule */
        printf("[sic] adding schedule\n");
-       moldyn_add_schedule(&md,100,1.0e-15);
+       moldyn_add_schedule(&md,10000,1.0e-15);
 
        /* activate logging */
        printf("[sic] activate logging\n");
-       moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",1);
-       moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",1);
+       moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",100);
+       moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",100);
 
        /*
         * let's do the actual md algorithm now