From 72df64eacc634e315f2205fc7bd2406223f92bf3 Mon Sep 17 00:00:00 2001 From: hackbard Date: Wed, 13 Dec 2006 00:37:10 +0000 Subject: [PATCH] safety checkin --- math/math.h | 7 +++++- moldyn.c | 62 ++++++++++++++++++++++++++++++++++++++--------------- sic.c | 3 ++- 3 files changed, 53 insertions(+), 19 deletions(-) diff --git a/math/math.h b/math/math.h index 8ff1bec..147f4c0 100644 --- 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 --- 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_jkzeta_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;icount;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<=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 --- 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"); -- 2.20.1