safety checkin

diff --git a/math/math.h b/math/math.h
index 8ff1bec..147f4c0 100644 (file)
--- a/math/math.h
+++ b/math/math.h
@@ -80,7 +80,12 @@ static inline int v3_set(t_3dvec *vec,double *ptr) {
 
 static inline int v3_copy(t_3dvec *trg,t_3dvec *src) {
 
-       memcpy(trg,src,sizeof(t_3dvec));
+       //memcpy(trg,src,sizeof(t_3dvec));
+
+       /* faster this way? */
+       trg->x=src->x;
+       trg->y=src->y;
+       trg->z=src->z;
 
        return 0;
 }

diff --git a/moldyn.c b/moldyn.c
index 44b53b6..73becde 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -802,14 +802,16 @@ int velocity_verlet(t_moldyn *moldyn) {
                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);
+//if(i==5) printf("v: %f %f %f\n",atom[i].r.x,(atom[i].r.x+moldyn->dim.x/2)/moldyn->lc.x,2*atom[i].r.x/moldyn->dim.x);
                check_per_bound(moldyn,&(atom[i].r));
+//if(i==5) printf("n: %f %f %f\n",atom[i].r.x,(atom[i].r.x+moldyn->dim.x/2)/moldyn->lc.x,2*atom[i].r.x/moldyn->dim.x);
 
                /* velocities */
                v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
                v3_add(&(atom[i].v),&(atom[i].v),&delta);
        }
 
-moldyn_bc_check(moldyn);
+//moldyn_bc_check(moldyn);
        /* neighbour list update */
        link_cell_update(moldyn);
 
@@ -841,7 +843,6 @@ int potential_force_calc(t_moldyn *moldyn) {
        t_linkcell *lc;
        t_list neighbour_i[27];
        t_list neighbour_i2[27];
-       //t_list neighbour_j[27];
        t_list *this,*that;
        u8 bc_ij,bc_ik;
        int dnlc;
@@ -988,7 +989,6 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
                if(a->z>=z) a->z-=dim->z;
                else if(-a->z>z) a->z+=dim->z;
        }
-printf("%f %f %f\n",a->x,x,a->x/x);
 
        return 0;
 }
@@ -1236,7 +1236,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                s_r=S-R;
                arg=M_PI*(d_ij-R)/s_r;
                f_c=0.5+0.5*cos(arg);
-               df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij));
+               //df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij)); /* MARK! */
+               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 */
@@ -1507,7 +1508,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                        s_r=S-R;
                        arg=M_PI*(d_ik-R)/s_r;
                        f_c_ik=0.5+0.5*cos(arg);
-                       df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik));
+                       //df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik)); /* MARK */
+                       df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 
                        /* zeta_ij */
                        exchange->zeta_ij+=f_c_ik*g;
@@ -1557,13 +1559,13 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                c2d2=exchange->cj2dj2;
 
                /* cosine of theta_jik by scalaproduct */
-               rr=v3_scalar_product(&dist_ij,&dist_jk); /* times -1 */
+               rr=-v3_scalar_product(&dist_ij,&dist_jk); /* -1, as ij -> ji */
                dd=d_ij*d_jk;
                cos_theta=rr/dd;
 
                /* d_costheta */
-               d_costheta1=1.0/(d_jk*d_ij);
-               d_costheta2=cos_theta/(d_ij*d_ij); /* in fact -cos(), but ^ */
+               d_costheta1=1.0/dd;
+               d_costheta2=cos_theta/(d_ij*d_ij);
 
                /* some usefull values */
                h_cos=(h-cos_theta);
@@ -1579,6 +1581,9 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                v3_add(&temp1,&temp1,&temp2);
                v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
 
+               /* store dg in temp2 and use it for dVjk later */
+               v3_copy(&temp2,&temp1);
+
                /* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */
                dzeta=&(exchange->dzeta_ji);
                if(d_jk<R) {
@@ -1600,7 +1605,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                        /* zeta_ji */
                        exchange->zeta_ji+=f_c_jk*g;
 
-                       /* dzeta_ij */
+                       /* dzeta_ji */
                        v3_scale(&temp1,&temp1,f_c_jk);
                        v3_add(dzeta,dzeta,&temp1);
                }
@@ -1608,19 +1613,19 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                /* dV_jk stuff | add force contribution on atom i immediately */
                if(exchange->d_ij_between_rs) {
                        zeta=f_c*g;
-                       v3_scale(&temp1,&temp1,f_c);
-                       v3_scale(&temp2,&dist_ij,df_c);
-                       v3_add(&temp1,&temp1,&temp2);
+                       v3_scale(&temp1,&temp2,f_c);
+                       v3_scale(&temp2,&dist_ij,df_c*g);
+                       v3_add(&temp2,&temp2,&temp1); /* -> dzeta_jk in temp2 */
                }
                else {
                        zeta=g;
-                       // dzeta_jk is simply dg, which is temp1
+                       // dzeta_jk is simply dg, which is stored in temp2
                }
                /* 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*0.5);
-               v3_add(&(ai->f),&(ai->f),&temp1); /* -1 skipped in f_a calc ^ */
+               tmp=-chi/2.0*pow((1+tmp*zeta),(-1.0/(2.0*n)-1))*tmp;
+               v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
+               v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
                                                  /* scaled with 0.5 ^ */
        }
 
@@ -1637,14 +1642,37 @@ int moldyn_bc_check(t_moldyn *moldyn) {
        t_atom *atom;
        t_3dvec *dim;
        int i;
+double x;
+u8 byte;
+int j,k;
 
        atom=moldyn->atom;
        dim=&(moldyn->dim);
+x=dim->x/2;
 
        for(i=0;i<moldyn->count;i++) {
-               if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2)
+               if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) {
                        printf("FATAL: atom %d: x: %.20f (%.20f)\n",
                               i,atom[i].r.x,dim->x/2);
+                       printf("diagnostic:\n");
+                       for(j=0;j<8;j++) {
+                               memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
+                               for(k=0;k<8;k++)
+                                       printf("%d%c",
+                                       ((byte)&(1<<k))?1:0,
+                                       (k==7)?'\n':'|');
+                       }
+                       printf("---------------\n");
+                       for(j=0;j<8;j++) {
+                               memcpy(&byte,(u8 *)(&x)+j,1);
+                               for(k=0;k<8;k++)
+                                       printf("%d%c",
+                                       ((byte)&(1<<k))?1:0,
+                                       (k==7)?'\n':'|');
+                       }
+                       if(atom[i].r.x==x) printf("hier gleich!\n");
+                       else printf("hier NICHT gleich!\n");
+               }
                if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
                        printf("FATAL: atom %d: y: %.20f (%.20f)\n",
                               i,atom[i].r.y,dim->y/2);

diff --git a/sic.c b/sic.c
index 5edb699..bc10e31 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -132,7 +132,8 @@ int main(int argc,char **argv) {
        /* set temperature */
        printf("[sic] setting temperature\n");
        //set_temperature(&md,273.0+450.0);
-       set_temperature(&md,0.0);
+       set_temperature(&md,1.0);
+       //set_temperature(&md,273.0);
 
        /* set p/t scaling */
        printf("[sic] set p/t scaling\n");