security checkin

diff --git a/moldyn.c b/moldyn.c
index 6736e45..1f6391f 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -681,7 +681,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
        /* zero absolute time */
        moldyn->time=0.0;
 
-       /* debugging, ignre */
+       /* debugging, ignore */
        moldyn->debug=0;
 
        /* executing the schedule */
@@ -741,7 +741,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
                        if(!(i%v)) {
                                visual_atoms(&(moldyn->vis),moldyn->time,
                                             moldyn->atom,moldyn->count);
-                               printf("\rsched: %d, steps: %d, theta: %d",
+                               printf("\rsched: %d, steps: %d, debug: %d",
                                       sched,i,moldyn->debug);
                                fflush(stdout);
                        }
@@ -920,15 +920,14 @@ int potential_force_calc(t_moldyn *moldyn) {
                
                                /* 2bp post function */
                                if(moldyn->func2b_post) {
-printf("DEBUG: vor 2bp post\n");
                                        moldyn->func2b_post(moldyn,
                                                            &(itom[i]),
                                                            jtom,bc_ij);
-printf("DEBUG: nach 2bp post\n");
                                }
 
                        }
                }
+//printf("DEBUG: %.15f \n",itom[i].f.x);
        }
 
        return 0;
@@ -1119,6 +1118,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        /* clear 3bp and 2bp post run */
        exchange->run3bp=0;
        exchange->run2bp_post=0;
+
+       /* reset S > r > R mark */
+       exchange->d_ij_between_rs=0;
        
        /*
         * calc of 2bp contribution of V_ij and dV_ij/ji
@@ -1208,6 +1210,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij));
                /* two body contribution (ij, ji) */
                v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
+               /* tell 3bp that S > r > R */
+               exchange->d_ij_between_rs=1;
        }
 
        /* add forces of 2bp (ij, ji) contribution
@@ -1281,40 +1285,55 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        
        /* Vij and dVij */
        zeta=exchange->zeta_ij;
-       tmp=betaini*pow(zeta,ni-1.0);           /* beta^n * zeta^n-1 */
-       b=(1+zeta*tmp);                         /* 1 + beta^n * zeta^n */
-       db=chi*pow(b,-1.0/(2*ni)-1);            /* chi * (...)^(-1/2n - 1) */
-       b=db*b;                                 /* b_ij */
-       db*=-0.5*tmp;                           /* db_ij */
-       v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
-       v3_scale(&temp,dist_ij,df_a*b);
-       v3_add(&force,&force,&temp);
-       v3_scale(&force,&force,f_c);
-       v3_scale(&temp,dist_ij,df_c*b*f_a);
-       v3_add(&force,&force,&temp);
+       if(zeta==0.0) {
+               moldyn->debug++;                /* just for debugging ... */
+               db=0.0;
+               b=chi;
+       }
+       else {
+               tmp=betaini*pow(zeta,ni-1.0);           /* beta^n * zeta^n-1 */
+               b=(1+zeta*tmp);                         /* 1 + beta^n zeta^n */
+               db=chi*pow(b,-1.0/(2*ni)-1);            /* x(...)^(-1/2n - 1) */
+               b=db*b;                                 /* b_ij */
+               db*=-0.5*tmp;                           /* db_ij */
+               v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
+               v3_scale(&temp,dist_ij,df_a*b);
+               v3_add(&force,&force,&temp);
+               v3_scale(&force,&force,f_c);
+               v3_scale(&temp,dist_ij,df_c*b*f_a);
+               v3_add(&force,&force,&temp);
+               v3_scale(&force,&force,-0.5);
+
+               /* add force */
+               v3_add(&(ai->f),&(ai->f),&force);
+       }
 
        /* add energy of 3bp sum */
        moldyn->energy+=(0.5*f_c*b*f_a);
 
-       /* add force (sub, as F = - dVij) */
-       v3_sub(&(ai->f),&(ai->f),&force);
-
        /* dVji */
        zeta=exchange->zeta_ji;
-       tmp=betajnj*pow(zeta,nj-1.0);           /* beta^n * zeta^n-1 */
-       b=(1+zeta*tmp);                         /* 1 + beta^n * zeta^n */
-       db=chi*pow(b,-1.0/(2*nj)-1);            /* chi * (...)^(-1/2n - 1) */
-       b=db*b;                                 /* b_ij */
-       db*=-0.5*tmp;                           /* db_ij */
-       v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
-       v3_scale(&temp,dist_ij,df_a*b);
-       v3_add(&force,&force,&temp);
-       v3_scale(&force,&force,f_c);
-       v3_scale(&temp,dist_ij,df_c*b*f_a);
-       v3_add(&force,&force,&temp);
-
-       /* add force (sub, as F = - dVji) */
-       v3_sub(&(ai->f),&(ai->f),&force);
+       if(zeta==0.0) {
+               // do nothing ... (db=0, b=chi)
+               moldyn->debug++;
+       }
+       else {
+               tmp=betajnj*pow(zeta,nj-1.0);           /* beta^n * zeta^n-1 */
+               b=(1+zeta*tmp);                         /* 1 + beta^n zeta^n */
+               db=chi*pow(b,-1.0/(2*nj)-1);            /* x(...)^(-1/2n - 1) */
+               b=db*b;                                 /* b_ij */
+               db*=-0.5*tmp;                           /* db_ij */
+               v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
+               v3_scale(&temp,dist_ij,df_a*b);
+               v3_add(&force,&force,&temp);
+               v3_scale(&force,&force,f_c);
+               v3_scale(&temp,dist_ij,0.5*df_c*b*f_a);
+               v3_add(&force,&force,&temp);
+               v3_scale(&force,&force,-0.5);
+
+               /* add force */
+               v3_sub(&(ai->f),&(ai->f),&force);
+       }
 
        return 0;
 }
@@ -1569,8 +1588,9 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                /* betajnj * zeta_jk ^ nj-1 */
                tmp=exchange->betajnj*pow(zeta,(n-1.0));
                tmp=-chi/2.0*pow(1+tmp*zeta,-1.0/(2.0*n)-1)*tmp;
-               v3_scale(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk);
+               v3_scale(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
                v3_add(&(ai->f),&(ai->f),&temp1); /* -1 skipped in f_a calc ^ */
+                                                 /* scaled with 0.5 ^ */
        }
 
        return 0;

diff --git a/sic.c b/sic.c
index d141782..cbb81a3 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -26,6 +26,9 @@ int main(int argc,char **argv) {
        /* misc parameters */
        double tau;
 
+       /* testing location & velocity vector */
+       t_3dvec r,v;
+
        /* values */
        tau=1.0e-15;    /* delta t = 1 fs */
 
@@ -93,10 +96,11 @@ int main(int argc,char **argv) {
        /* cutoff radius */
        printf("[sic] setting cutoff radius\n");
        set_cutoff(&md,TM_S_SI);
+       //set_cutoff(&md,2*LC_SI);
 
        /* set (initial) dimensions of simulation volume */
        printf("[sic] setting dimensions\n");
-       set_dim(&md,2*LC_SI,2*LC_SI,2*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");
@@ -104,10 +108,20 @@ int main(int argc,char **argv) {
 
        /* create the lattice / place atoms */
        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_2BP|ATOM_ATTR_HB,
-                      0,2,2,2);
+       //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+        //              ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+         //             //ATOM_ATTR_2BP|ATOM_ATTR_HB,
+          //            0,5,5,5);
+
+       /* testing configuration */
+       r.x=-0.55*0.25*sqrt(3.0)*LC_SI; v.x=0;
+       r.y=0;                  v.y=0;
+       r.z=0;                  v.z=0;
+       add_atom(&md,SI,M_SI,0,ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,&r,&v);
+       r.x=+0.55*0.25*sqrt(3.0)*LC_SI; v.x=0;
+       r.y=0;                  v.y=0;
+       r.z=0;                  v.z=0;
+       add_atom(&md,SI,M_SI,0,ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,&r,&v);
 
        /* setting a nearest neighbour distance for the moldyn checks */
        set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */
@@ -127,12 +141,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,1000,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",10);
+       moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",10);
 
        /*
         * let's do the actual md algorithm now

diff --git a/visual/visual.c b/visual/visual.c
index 88bc6e6..687c561 100644 (file)
--- a/visual/visual.c
+++ b/visual/visual.c
@@ -93,7 +93,7 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
        dprintf(v->fd,"set ambient 20\n");
        dprintf(v->fd,"set specular on\n");
        dprintf(v->fd,"set boundbox on\n");
-       //dprintf(v->fd,"label\n");
+       dprintf(v->fd,"label\n");
        sprintf(file,"%s-%.15f.ppm",v->fb,time);
        dprintf(v->fd,"write ppm %s\n",file);
        dprintf(v->fd,"zap\n");