From: hackbard Date: Sun, 17 Dec 2006 01:29:39 +0000 (+0000) Subject: testing foo X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=150c5f804a5088a01ef2bb1e934b5e65d10ac23e;p=physik%2Fposic.git testing foo --- diff --git a/moldyn.c b/moldyn.c index 5aa6d87..db91859 100644 --- a/moldyn.c +++ b/moldyn.c @@ -519,7 +519,8 @@ int scale_volume(t_moldyn *moldyn) { t_atom *atom; t_3dvec *dim,*vdim; - double virial,scale; + double scale,v; + t_virial virial; t_linkcell *lc; int i; @@ -528,12 +529,23 @@ int scale_volume(t_moldyn *moldyn) { vdim=&(moldyn->vis.dim); lc=&(moldyn->lc); - for(i=0;icount;i++) - virial+=v3_norm(&(atom[i].virial)); + memset(&virial,0,sizeof(t_virial)); + + for(i=0;icount;i++) { + virial.xx+=atom[i].virial.xx; + virial.yy+=atom[i].virial.yy; + virial.zz+=atom[i].virial.zz; + virial.xy+=atom[i].virial.xy; + virial.xz+=atom[i].virial.xz; + virial.yz+=atom[i].virial.yz; + } + + /* just a guess so far ... */ + v=sqrt(virial.xx*virial.xx+virial.yy*virial.yy+virial.zz+virial.zz); -printf("%f\n",virial); +printf("%f\n",v); /* get pressure from virial */ - moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*virial; + moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*v; moldyn->p/=moldyn->volume; printf("%f\n",moldyn->p/(ATM)); @@ -984,6 +996,7 @@ int potential_force_calc(t_moldyn *moldyn) { int i,j,k,count; t_atom *itom,*jtom,*ktom; + t_virial *virial; t_linkcell *lc; t_list neighbour_i[27]; t_list neighbour_i2[27]; @@ -1005,7 +1018,13 @@ int potential_force_calc(t_moldyn *moldyn) { v3_zero(&(itom[i].f)); /* reset viral of atom i */ - v3_zero(&(itom[i].virial)); + virial=&(itom[i].virial); + virial->xx=0.0; + virial->yy=0.0; + virial->zz=0.0; + virial->xy=0.0; + virial->xz=0.0; + virial->yz=0.0; /* reset site energy */ itom[i].e=0.0; @@ -1224,7 +1243,10 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { int tersoff_mult_complete_params(t_tersoff_mult_params *p) { printf("[moldyn] tersoff parameter completion\n"); + p->S2[0]=p->S[0]*p->S[0]; + p->S2[1]=p->S[1]*p->S[1]; p->Smixed=sqrt(p->S[0]*p->S[1]); + p->S2mixed=p->Smixed*p->Smixed; p->Rmixed=sqrt(p->R[0]*p->R[1]); p->Amixed=sqrt(p->A[0]*p->A[1]); p->Bmixed=sqrt(p->B[0]*p->B[1]); @@ -1285,8 +1307,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_tersoff_mult_params *params; t_tersoff_exchange *exchange; t_3dvec dist_ij,force; - double d_ij; - double A,B,R,S,lambda,mu; + double d_ij,d_ij2; + double A,B,R,S,S2,lambda,mu; double f_r,df_r; double f_c,df_c; int brand; @@ -1317,18 +1339,10 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { * */ - /* dist_ij, d_ij */ - v3_sub(&dist_ij,&(aj->r),&(ai->r)); - if(bc) check_per_bound(moldyn,&dist_ij); - d_ij=v3_norm(&dist_ij); - - /* save for use in 3bp */ - exchange->d_ij=d_ij; - exchange->dist_ij=dist_ij; - /* constants */ if(brand==ai->brand) { S=params->S[brand]; + S2=params->S2[brand]; R=params->R[brand]; A=params->A[brand]; B=params->B[brand]; @@ -1338,6 +1352,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { } else { S=params->Smixed; + S2=params->S2mixed; R=params->Rmixed; A=params->Amixed; B=params->Bmixed; @@ -1346,10 +1361,23 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { params->exchange.chi=params->chi; } - /* if d_ij > S => no force & potential energy contribution */ - if(d_ij>S) + /* dist_ij, d_ij */ + v3_sub(&dist_ij,&(aj->r),&(ai->r)); + if(bc) check_per_bound(moldyn,&dist_ij); + d_ij2=v3_absolute_square(&dist_ij); + + /* if d_ij2 > S2 => no force & potential energy contribution */ + if(d_ij2>S2) return 0; + /* now we will need the distance */ + //d_ij=v3_norm(&dist_ij); + d_ij=sqrt(d_ij2); + + /* save for use in 3bp */ + exchange->d_ij=d_ij; + exchange->dist_ij=dist_ij; + /* more constants */ exchange->beta_j=&(params->beta[brand]); exchange->n_j=&(params->n[brand]); @@ -1389,7 +1417,6 @@ 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)); /* 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); @@ -1685,7 +1712,6 @@ 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)); /* MARK */ df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik)); /* zeta_ij */ diff --git a/moldyn.h b/moldyn.h index 38e4999..dd94b0e 100644 --- a/moldyn.h +++ b/moldyn.h @@ -22,12 +22,22 @@ /* general */ typedef unsigned char u8; +/* virial */ +typedef struct s_virial { + double xx; /* | xx xy xz | */ + double yy; /* V = | yx yy yz | */ + double zz; /* | zx zy zz | */ + double xy; /* */ + double xz; /* with: xy=yx, xz=zx, yz=zy */ + double yz; /* */ +} t_virial; + /* the atom of the md simulation */ typedef struct s_atom { t_3dvec r; /* position */ t_3dvec v; /* velocity */ t_3dvec f; /* force */ - t_3dvec virial; /* virial (v_xx, v_yy, v_zz) */ + t_virial virial; /* virial */ double e; /* site energy */ int element; /* number of element in pse */ double mass; /* atom mass */ @@ -222,8 +232,10 @@ typedef struct s_tersoff_exchange { /* tersoff multi (2!) potential parameters */ typedef struct s_tersoff_mult_params { double S[2]; /* tersoff cutoff radii */ + double S2[2]; /* tersoff cutoff radii squared */ double R[2]; /* tersoff cutoff radii */ double Smixed; /* mixed S radius */ + double S2mixed; /* mixed S radius squared */ double Rmixed; /* mixed R radius */ double A[2]; /* factor of tersoff attractive part */ double B[2]; /* factor of tersoff repulsive part */ @@ -310,6 +322,7 @@ typedef struct s_tersoff_mult_params { #define SI 0x0e #define LC_SI (0.543105e-9*METER) /* A */ #define M_SI 28.08553 /* amu */ + #define LJ_SIGMA_SI ((0.25*sqrt(3.0)*LC_SI)/1.122462) /* A */ #define LJ_EPSILON_SI (2.1678*EV) /* NA */ diff --git a/sic.c b/sic.c index 35931e0..c4355b3 100644 --- a/sic.c +++ b/sic.c @@ -112,24 +112,26 @@ 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, - 0,5,5,5); - moldyn_bc_check(&md); + //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); + //moldyn_bc_check(&md); /* testing configuration */ - //r.x=2.95/2; 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=-2.95/2; 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=2.8/2; 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, + // ATOM_ATTR_2BP, + &r,&v); + r.x=-2.8/2; 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, + // ATOM_ATTR_2BP, + &r,&v); /* setting a nearest neighbour distance for the moldyn checks */ set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */ @@ -151,25 +153,21 @@ int main(int argc,char **argv) { printf("[sic] set p/t scaling\n"); //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0, // T_SCALE_BERENDSEN,100.0); - set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0); + //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0); /* initial thermal fluctuations of particles (in equilibrium) */ printf("[sic] thermal init\n"); - thermal_init(&md,TRUE); + //thermal_init(&md,TRUE); /* create the simulation schedule */ printf("[sic] adding schedule\n"); - moldyn_add_schedule(&md,20000,.1); - moldyn_add_schedule(&md,10000,.2); - moldyn_add_schedule(&md,6667,.3); - moldyn_add_schedule(&md,5000,.4); - moldyn_add_schedule(&md,4001,.5); + moldyn_add_schedule(&md,10000,.1); /* activate logging */ printf("[sic] activate logging\n"); moldyn_set_log_dir(&md,argv[1]); - moldyn_set_log(&md,LOG_TOTAL_ENERGY,100); - moldyn_set_log(&md,VISUAL_STEP,100); + moldyn_set_log(&md,LOG_TOTAL_ENERGY,1); + moldyn_set_log(&md,VISUAL_STEP,20); /* * let's do the actual md algorithm now