From: hackbard Date: Wed, 7 Mar 2007 16:23:21 +0000 (+0000) Subject: tersoff MAYBE fixed now! X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e25ff194682ff5fac86c60701343103e74973bed;p=physik%2Fposic.git tersoff MAYBE fixed now! --- diff --git a/moldyn.c b/moldyn.c index dae975c..ca6f84a 100644 --- a/moldyn.c +++ b/moldyn.c @@ -162,28 +162,28 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) { return 0; } -int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params) { +int set_potential1b(t_moldyn *moldyn,pf_func1b func) { moldyn->func1b=func; return 0; } -int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) { +int set_potential2b(t_moldyn *moldyn,pf_func2b func) { moldyn->func2b=func; return 0; } -int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) { +int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func) { moldyn->func2b_post=func; return 0; } -int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) { +int set_potential3b(t_moldyn *moldyn,pf_func3b func) { moldyn->func3b=func; diff --git a/moldyn.h b/moldyn.h index 1300998..b140372 100644 --- a/moldyn.h +++ b/moldyn.h @@ -230,9 +230,9 @@ typedef struct s_moldyn { #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_SIGMA_SI ((0.25*sqrt(3.0)*LC_SI)/1.122462) /* A */ //#define LJ_SIGMA_SI (LC_SI/1.122462) /* A */ -#define LJ_SIGMA_SI (0.5*sqrt(2.0)*LC_SI/1.122462) /* A */ +//#define LJ_SIGMA_SI (0.5*sqrt(2.0)*LC_SI/1.122462) /* A */ #define LJ_EPSILON_SI (2.1678*EV) /* NA */ #define TM_R_SI (2.7e-10*METER) /* A */ @@ -292,10 +292,11 @@ int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc); int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize); int set_nn_dist(t_moldyn *moldyn,double dist); int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z); -int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params); -int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params); -int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params); -int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params); +int set_potential1b(t_moldyn *moldyn,pf_func1b func); +int set_potential2b(t_moldyn *moldyn,pf_func2b func); +int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func); +int set_potential3b(t_moldyn *moldyn,pf_func3b func); +int set_potential_params(t_moldyn *moldyn,void *params); int moldyn_set_log_dir(t_moldyn *moldyn,char *dir); int moldyn_set_report(t_moldyn *moldyn,char *author,char *title); diff --git a/potentials/lennard_jones.c b/potentials/lennard_jones.c index 1738ea5..1277c46 100644 --- a/potentials/lennard_jones.c +++ b/potentials/lennard_jones.c @@ -35,10 +35,11 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { v3_sub(&distance,&(aj->r),&(ai->r)); if(bc) check_per_bound(moldyn,&distance); - d=v3_absolute_square(&distance); /* 1/r^2 */ + d=v3_absolute_square(&distance); /* r^2 */ if(d<=moldyn->cutoff_square) { d=1.0/d; /* 1/r^2 */ h2=d*d; /* 1/r^4 */ + h2*=d; /* 1/r^6 */ h1=h2*h2; /* 1/r^12 */ moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc); h2*=d; /* 1/r^8 */ diff --git a/potentials/tersoff.c b/potentials/tersoff.c index 972075d..01364ea 100644 --- a/potentials/tersoff.c +++ b/potentials/tersoff.c @@ -95,6 +95,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { double s_r; double arg; + /* use newtons third law */ + //if(aipot_params; brand=aj->brand; exchange=&(params->exchange); @@ -119,27 +122,11 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { * */ - /* constants */ - if(brand==ai->brand) { - S=params->S[brand]; + /* determine cutoff square */ + if(brand==ai->brand) S2=params->S2[brand]; - R=params->R[brand]; - A=params->A[brand]; - B=params->B[brand]; - lambda=params->lambda[brand]; - mu=params->mu[brand]; - exchange->chi=1.0; - } - else { - S=params->Smixed; + else S2=params->S2mixed; - R=params->Rmixed; - A=params->Amixed; - B=params->Bmixed; - lambda=params->lambda_m; - mu=params->mu_m; - params->exchange.chi=params->chi; - } /* dist_ij, d_ij */ v3_sub(&dist_ij,&(aj->r),&(ai->r)); @@ -166,12 +153,26 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { exchange->d_j=&(params->d[brand]); exchange->h_j=&(params->h[brand]); if(brand==ai->brand) { + S=params->S[brand]; + R=params->R[brand]; + A=params->A[brand]; + B=params->B[brand]; + lambda=params->lambda[brand]; + mu=params->mu[brand]; + exchange->chi=1.0; exchange->betajnj=exchange->betaini; exchange->cj2=exchange->ci2; exchange->dj2=exchange->di2; exchange->cj2dj2=exchange->ci2di2; } else { + S=params->Smixed; + R=params->Rmixed; + A=params->Amixed; + B=params->Bmixed; + lambda=params->lambda_m; + mu=params->mu_m; + params->exchange.chi=params->chi; exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j)); exchange->cj2=params->c[brand]*params->c[brand]; exchange->dj2=params->d[brand]*params->d[brand]; diff --git a/sic.c b/sic.c index e148709..167735f 100644 --- a/sic.c +++ b/sic.c @@ -51,6 +51,8 @@ int main(int argc,char **argv) { /* testing location & velocity vector */ t_3dvec r,v; + memset(&r,0,sizeof(t_3dvec)); + memset(&v,0,sizeof(t_3dvec)); /* initialize moldyn */ moldyn_init(&md,argc,argv); @@ -59,17 +61,20 @@ int main(int argc,char **argv) { set_int_alg(&md,MOLDYN_INTEGRATE_VERLET); /* choose potential */ - //set_potential1b(&md,tersoff_mult_1bp,&tp); - //set_potential2b(&md,tersoff_mult_2bp,&tp); - //set_potential2b_post(&md,tersoff_mult_post_2bp,&tp); - //set_potential3b(&md,tersoff_mult_3bp,&tp); - set_potential2b(&md,lennard_jones,&lj); - //set_potential2b(&md,harmonic_oscillator,&ho); + set_potential1b(&md,tersoff_mult_1bp); + set_potential2b(&md,tersoff_mult_2bp); + set_potential2b_post(&md,tersoff_mult_post_2bp); + set_potential3b(&md,tersoff_mult_3bp); + //set_potential2b(&md,lennard_jones); + //set_potential2b(&md,harmonic_oscillator); + set_potential_params(&md,&tp); + //set_potential_params(&md,&lj); + //set_potential_params(&md,&ho); /* cutoff radius */ - //set_cutoff(&md,TM_S_SI); - //set_cutoff(&md,2*LC_SI*0.5*sqrt(1.5)); - set_cutoff(&md,2.0*LC_SI); + set_cutoff(&md,TM_S_SI); + //set_cutoff(&md,LC_SI*sqrt(3.0)); + //set_cutoff(&md,2.0*LC_SI); /* * potential parameters @@ -83,8 +88,8 @@ int main(int argc,char **argv) { lj.uc=lj.epsilon4*(lj.sigma12/pow(md.cutoff,12.0)-lj.sigma6/pow(md.cutoff,6)); /* harmonic oscillator */ - //ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI; - ho.equilibrium_distance=LC_SI; + ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI; + //ho.equilibrium_distance=LC_SI; ho.spring_constant=LJ_EPSILON_SI; /* @@ -126,16 +131,16 @@ int main(int argc,char **argv) { /* create the lattice / place atoms */ //create_lattice(&md,CUBIC,LC_SI,SI,M_SI, - create_lattice(&md,FCC,LC_SI,SI,M_SI, - //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, + //create_lattice(&md,FCC,LC_SI,SI,M_SI, + 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,6,6,6); moldyn_bc_check(&md); /* testing configuration */ - //r.x=0.28*sqrt(3)*LC_SI/2; v.x=0; - //r.x=1.75*LC_SI; v.x=-0.01; + //r.x=0.27*sqrt(3.0)*LC_SI/2.0; v.x=0; + //r.x=(TM_S_SI+TM_R_SI)/4.0; v.x=0; //r.y=0; v.y=0; //r.z=0; v.z=0; //add_atom(&md,SI,M_SI,0, @@ -149,6 +154,13 @@ int main(int argc,char **argv) { // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, // ATOM_ATTR_2BP|ATOM_ATTR_HB, // &r,&v); + //r.x=0; v.x=0; + //r.y=0; v.y=0; + //r.z=sin(M_PI*60.0/180.0)*(TM_S_SI+TM_R_SI)/4.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|ATOM_ATTR_HB, + // &r,&v); /* set temperature & pressure */ set_temperature(&md,atof(argv[2])+273.0); @@ -164,7 +176,7 @@ int main(int argc,char **argv) { thermal_init(&md,TRUE); /* create the simulation schedule */ - moldyn_add_schedule(&md,101,1.0); + moldyn_add_schedule(&md,10001,1.0); //moldyn_add_schedule(&md,501,1.0); //moldyn_add_schedule(&md,501,1.0); @@ -175,7 +187,7 @@ int main(int argc,char **argv) { moldyn_set_log_dir(&md,argv[1]); moldyn_set_report(&md,"Frank Zirkelbach","Test 1"); moldyn_set_log(&md,LOG_TOTAL_ENERGY,1); - moldyn_set_log(&md,VISUAL_STEP,1); + moldyn_set_log(&md,VISUAL_STEP,100); moldyn_set_log(&md,CREATE_REPORT,0); /*