From bb3e6f4ab3dfea9083e0d5b7d403ee6197f4041d Mon Sep 17 00:00:00 2001 From: hackbard Date: Mon, 11 Dec 2006 16:03:33 +0000 Subject: [PATCH] bc check probs --- moldyn.c | 62 +++++++++++++++++++++++++++++++++++++++++++------------- moldyn.h | 2 ++ sic.c | 29 +++++++++++++------------- 3 files changed, 64 insertions(+), 29 deletions(-) diff --git a/moldyn.c b/moldyn.c index 1f6391f..4ddc42e 100644 --- a/moldyn.c +++ b/moldyn.c @@ -781,12 +781,15 @@ 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==0) printf("und dann %.20f\n",atom[i].r.x*1e10); check_per_bound(moldyn,&(atom[i].r)); +//if(i==0) printf("und danach %.20f\n",atom[i].r.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); /* neighbour list update */ link_cell_update(moldyn); @@ -1189,7 +1192,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* f_r_ij = f_r_ji, df_r_ij = df_r_ji */ f_r=A*exp(-lambda*d_ij); - df_r=-lambda*f_r/d_ij; + df_r=lambda*f_r/d_ij; /* f_a, df_a calc (again, same for ij and ji) | save for later use! */ exchange->f_a=-B*exp(-mu*d_ij); @@ -1289,6 +1292,7 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { moldyn->debug++; /* just for debugging ... */ db=0.0; b=chi; + v3_scale(&force,dist_ij,df_a*b*f_c); } else { tmp=betaini*pow(zeta,ni-1.0); /* beta^n * zeta^n-1 */ @@ -1300,13 +1304,13 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { 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); } + 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); @@ -1314,8 +1318,9 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* dVji */ zeta=exchange->zeta_ji; if(zeta==0.0) { - // do nothing ... (db=0, b=chi) moldyn->debug++; + b=chi; + v3_scale(&force,dist_ij,df_a*b*f_c); } else { tmp=betajnj*pow(zeta,nj-1.0); /* beta^n * zeta^n-1 */ @@ -1327,13 +1332,13 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { 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); } + v3_scale(&temp,dist_ij,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; } @@ -1596,3 +1601,32 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { return 0; } + +/* + * debugging / critical check functions + */ + +int moldyn_bc_check(t_moldyn *moldyn) { + + t_atom *atom; + t_3dvec *dim; + int i; + + atom=moldyn->atom; + dim=&(moldyn->dim); + + for(i=0;icount;i++) { + 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*1e10,dim->x/2*1e10); + 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*1e10,dim->y/2*1e10); + if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2) + printf("FATAL: atom %d: z: %.20f (%.20f)\n", + i,atom[i].r.z*1e10,dim->z/2*1e10); + } + + return 0; +} + diff --git a/moldyn.h b/moldyn.h index d62ec47..534287b 100644 --- a/moldyn.h +++ b/moldyn.h @@ -388,4 +388,6 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); +int moldyn_bc_check(t_moldyn *moldyn); + #endif diff --git a/sic.c b/sic.c index cbb81a3..0d3b798 100644 --- a/sic.c +++ b/sic.c @@ -108,20 +108,19 @@ 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,5,5,5); + create_lattice(&md,DIAMOND,LC_SI,SI,M_SI, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|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); + //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 ! */ @@ -141,12 +140,12 @@ int main(int argc,char **argv) { /* create the simulation schedule */ printf("[sic] adding schedule\n"); - moldyn_add_schedule(&md,1000,1.0e-15); + moldyn_add_schedule(&md,100,1.0e-15); /* activate logging */ printf("[sic] activate logging\n"); - moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",10); - moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",10); + moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",1); + moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",1); /* * let's do the actual md algorithm now -- 2.20.1