From dece53fc37f9ebb52b33c9743333c213be2d6f26 Mon Sep 17 00:00:00 2001 From: hackbard Date: Thu, 30 Nov 2006 15:55:32 +0000 Subject: [PATCH] into the tersoff potential --- moldyn.c | 31 +++++++++++++++++++++------- moldyn.h | 28 ++++++++++++++++++++++++++ sic.c | 61 +++++++++++++++++++++++++++++++++++++------------------- 3 files changed, 92 insertions(+), 28 deletions(-) diff --git a/moldyn.c b/moldyn.c index d52fcd3..ed831ae 100644 --- a/moldyn.c +++ b/moldyn.c @@ -231,13 +231,15 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, printf("[moldyn] created lattice with %d atoms\n",count); while(count) { - moldyn->atom[count-1].element=element; - moldyn->atom[count-1].mass=mass; - moldyn->atom[count-1].attr=attr; - moldyn->atom[count-1].bnum=bnum; count-=1; + moldyn->atom[count].element=element; + moldyn->atom[count].mass=mass; + moldyn->atom[count].attr=attr; + moldyn->atom[count].bnum=bnum; + check_per_bound(moldyn,&(moldyn->atom[count].r)); } + return ret; } @@ -529,7 +531,7 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { } } - lc->dnlc=count2; + lc->dnlc=count1; lc->countn=27; return count2; @@ -618,7 +620,6 @@ int moldyn_integrate(t_moldyn *moldyn) { /* sqaure of some variables */ moldyn->tau_square=moldyn->tau*moldyn->tau; moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; - /* calculate initial forces */ potential_force_calc(moldyn); @@ -774,7 +775,7 @@ int potential_force_calc(t_moldyn *moldyn) { moldyn->energy=0.0; for(i=0;iattr&ATOM_ATTR_3BP)) continue; +printf("DEBUG: problem exists here ...\n"); link_cell_neighbour_index(moldyn, (btom->r.x+moldyn->dim.x/2)/lc->x, (btom->r.y+moldyn->dim.y/2)/lc->y, (btom->r.z+moldyn->dim.z/2)/lc->z, neighbourk); +printf("DEBUG: as you won't see that!\n"); for(k=0;kcountn;k++) { @@ -971,6 +974,20 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { * tersoff potential & force for 2 sorts of atoms */ +/* create mixed terms from parameters and set them */ +int tersoff_mult_complete_params(t_tersoff_mult_params *p) { + + printf("[moldyn] tersoff parameter completion\n"); + p->Smixed=sqrt(p->S[0]*p->S[1]); + 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]); + p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]); + p->mu_m=0.5*(p->mu[0]+p->mu[1]); + + return 0; +} + /* tersoff 1 body part */ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) { diff --git a/moldyn.h b/moldyn.h index 63dd14f..2d4b563 100644 --- a/moldyn.h +++ b/moldyn.h @@ -254,6 +254,7 @@ typedef struct s_tersoff_mult_params { #define K_BOLTZMANN 1.3807e-27 /* Nm/K */ #define AMU 1.660540e-27 /* kg */ +#define EV 1.60217733e-19 /* Nm */ #define FCC 0x01 #define DIAMOND 0x02 @@ -267,6 +268,32 @@ typedef struct s_tersoff_mult_params { #define LJ_SIGMA_SI ((0.25*sqrt(3.0)*LC_SI)/1.122462) /* m */ #define LJ_EPSILON_SI (2.1678*1.60e-19) /* Nm */ +#define TM_R_SI 2.7e-10 /* m */ +#define TM_S_SI 3.0e-10 /* m */ +#define TM_A_SI (1830.8*EV) /* Nm */ +#define TM_B_SI (471.18*EV) /* Nm */ +#define TM_LAMBDA_SI 2.4799e10 /* 1/m */ +#define TM_MU_SI 1.7322e10 /* 1/m */ +#define TM_BETA_SI 1.1000e-6 +#define TM_N_SI 0.78734 +#define TM_C_SI 1.0039e5 +#define TM_D_SI 1.62170 +#define TM_H_SI (-0.59825) + +#define TM_R_C 1.8e-10 /* m */ +#define TM_S_C 2.1e-10 /* m */ +#define TM_A_C (1393.6*EV) /* Nm */ +#define TM_B_C (346.7*EV) /* Nm */ +#define TM_LAMBDA_C 3.4879e10 /* 1/m */ +#define TM_MU_C 2.2119e10 /* 1/m */ +#define TM_BETA_C 1.5724e-7 +#define TM_N_C 0.72751 +#define TM_C_C 3.8049e4 +#define TM_D_C 4.384 +#define TM_H_C (-0.57058) + +#define TM_CHI_SIC 0.9776 + /* * @@ -325,6 +352,7 @@ int potential_force_calc(t_moldyn *moldyn); int check_per_bound(t_moldyn *moldyn,t_3dvec *a); int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); +int tersoff_mult_complete_params(t_tersoff_mult_params *p); int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai); int tersoff_mult_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); diff --git a/sic.c b/sic.c index 7c40a7e..2cd32eb 100644 --- a/sic.c +++ b/sic.c @@ -15,7 +15,7 @@ #include "posic.h" int main(int argc,char **argv) { - +printf("%d\n",sizeof(t_atom)); /* main moldyn structure */ t_moldyn md; @@ -24,12 +24,6 @@ int main(int argc,char **argv) { t_ho_params ho; t_tersoff_mult_params tp; - /* misc variables, mainly to initialize stuff */ - t_3dvec r,v; - - /* temperature */ - double t; - /* initialize moldyn */ printf("[sic] moldyn init\n"); moldyn_init(&md,argc,argv); @@ -40,10 +34,10 @@ int main(int argc,char **argv) { /* choose potential */ printf("[sic] selecting potential\n"); - //set_potential1b(&md,tersoff_mult_1bp,&tp); - //set_potential2b(&md,tersoff_mult_2bp,&tp); - //set_potential3b(&md,tersoff_mult_3bp,&tp); - set_potential2b(&md,lennard_jones,&lj); + set_potential1b(&md,tersoff_mult_1bp,&tp); + set_potential2b(&md,tersoff_mult_2bp,&tp); + set_potential3b(&md,tersoff_mult_3bp,&tp); + //set_potential2b(&md,lennard_jones,&lj); /* * potential parameters @@ -59,9 +53,39 @@ int main(int argc,char **argv) { ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI; ho.spring_constant=1; + /* + * tersoff mult potential parameters for SiC + */ + tp.S[0]=TM_S_SI; + tp.R[0]=TM_R_SI; + tp.A[0]=TM_A_SI; + tp.B[0]=TM_B_SI; + tp.lambda[0]=TM_LAMBDA_SI; + tp.mu[0]=TM_MU_SI; + tp.beta[0]=TM_BETA_SI; + tp.n[0]=TM_N_SI; + tp.c[0]=TM_C_SI; + tp.d[0]=TM_D_SI; + tp.h[0]=TM_H_SI; + + tp.S[1]=TM_S_C; + tp.R[1]=TM_R_C; + tp.A[1]=TM_A_C; + tp.B[1]=TM_B_C; + tp.lambda[1]=TM_LAMBDA_C; + tp.mu[1]=TM_MU_C; + tp.beta[1]=TM_BETA_C; + tp.n[1]=TM_N_C; + tp.c[1]=TM_C_C; + tp.d[1]=TM_D_C; + + tp.chi=TM_CHI_SIC; + + tersoff_mult_complete_params(&tp); + /* cutoff radius */ printf("[sic] setting cutoff radius\n"); - set_cutoff(&md,1*LC_SI); + set_cutoff(&md,TM_S_SI); /* set (initial) dimensions of simulation volume */ printf("[sic] setting dimensions\n"); @@ -73,21 +97,16 @@ int main(int argc,char **argv) { /* create the lattice / place atoms */ printf("[sic] creating atoms\n"); - //memset(&v,0,sizeof(t_3dvec)); - //r.y=0; - //r.z=0; - //r.x=0.23*sqrt(3.0)*LC_SI/2.0; - //add_atom(&md,SI,M_SI,0,ATOM_ATTR_2BP,&r,&v); - //r.x=-r.x; - //add_atom(&md,SI,M_SI,0,ATOM_ATTR_2BP,&r,&v); - create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,ATOM_ATTR_2BP,0,4,4,4); + create_lattice(&md,DIAMOND,LC_SI,SI,M_SI, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP, + 0,4,4,4); /* setting a nearest neighbour distance for the moldyn checks */ set_nn_dist(&md,sqrt(3.0)*LC_SI/4.0); /* diamond ! */ /* set temperature */ printf("[sic] setting temperature\n"); - set_temperature(&md,273.0); + set_temperature(&md,273.0+450.0); /* initial thermal fluctuations of particles */ printf("[sic] thermal init\n"); -- 2.39.2