From: hackbard Date: Mon, 11 Dec 2006 10:19:02 +0000 (+0000) Subject: security checkin X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=1d32be0a3b47add316f22f506b8595d56ac09a34;p=physik%2Fposic.git security checkin --- diff --git a/moldyn.c b/moldyn.c index 6736e45..1f6391f 100644 --- 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 --- 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 --- 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");