From 92ef07d77a4c879527180224acea73a3f6564497 Mon Sep 17 00:00:00 2001 From: hackbard Date: Wed, 30 May 2007 15:43:57 +0000 Subject: [PATCH] albe force calc bug fixed, sic mods, mean virial --- moldyn.c | 84 ++++++++++++++++++++++---------------------- moldyn.h | 25 +++++++++---- potentials/albe.c | 72 +++++++++++++++++-------------------- potentials/albe.h | 6 ++-- potentials/tersoff.c | 11 ++++++ sic.c | 79 +++++++++++++++++------------------------ 6 files changed, 139 insertions(+), 138 deletions(-) diff --git a/moldyn.c b/moldyn.c index 070ce69..e554a2f 100644 --- a/moldyn.c +++ b/moldyn.c @@ -803,8 +803,12 @@ double pressure_calc(t_moldyn *moldyn) { v+=(virial->xx+virial->yy+virial->zz); } + /* virial sum and mean virial */ + moldyn->virial_sum+=v; + moldyn->mean_v=moldyn->virial_sum/moldyn->total_steps; + /* assume up to date kinetic energy */ - moldyn->p=2.0*moldyn->ekin+v; + moldyn->p=2.0*moldyn->ekin+moldyn->mean_v; moldyn->p/=(3.0*moldyn->volume); moldyn->p_sum+=moldyn->p; moldyn->mean_p=moldyn->p_sum/moldyn->total_steps; @@ -823,8 +827,8 @@ double pressure_calc(t_moldyn *moldyn) { double thermodynamic_pressure_calc(t_moldyn *moldyn) { t_3dvec dim,*tp; - double u,p; - double scale,dv; + double u_up,u_down,dv; + double scale,p; t_atom *store; /* @@ -832,13 +836,11 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { * * => p = - dU/dV * - * dV: dx,y,z = 0.001 x,y,z */ - scale=1.00000000000001; -printf("\n\nP-DEBUG:\n"); + scale=0.00001; + dv=8*scale*scale*scale*moldyn->volume; - tp=&(moldyn->tp); store=malloc(moldyn->count*sizeof(t_atom)); if(store==NULL) { printf("[moldyn] allocating store mem failed\n"); @@ -846,59 +848,44 @@ printf("\n\nP-DEBUG:\n"); } /* save unscaled potential energy + atom/dim configuration */ - u=moldyn->energy; memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom)); dim=moldyn->dim; - /* derivative with respect to x direction */ - scale_dim(moldyn,scale,TRUE,0,0); - scale_atoms(moldyn,scale,TRUE,0,0); - dv=0.00000000000001*moldyn->dim.x*moldyn->dim.y*moldyn->dim.z; - link_cell_shutdown(moldyn); - link_cell_init(moldyn,QUIET); - potential_force_calc(moldyn); - tp->x=(moldyn->energy-u)/dv; - p=tp->x*tp->x; - - /* restore atomic configuration + dim */ - memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); - moldyn->dim=dim; - - /* derivative with respect to y direction */ - scale_dim(moldyn,scale,0,TRUE,0); - scale_atoms(moldyn,scale,0,TRUE,0); - dv=0.00000000000001*moldyn->dim.y*moldyn->dim.x*moldyn->dim.z; + /* scale up dimension and atom positions */ + scale_dim(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE); link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - tp->y=(moldyn->energy-u)/dv; - p+=tp->y*tp->y; + u_up=moldyn->energy; /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; - /* derivative with respect to z direction */ - scale_dim(moldyn,scale,0,0,TRUE); - scale_atoms(moldyn,scale,0,0,TRUE); - dv=0.00000000000001*moldyn->dim.z*moldyn->dim.x*moldyn->dim.y; + /* scale down dimension and atom positions */ + scale_dim(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE); link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - tp->z=(moldyn->energy-u)/dv; - p+=tp->z*tp->z; + u_down=moldyn->energy; + + /* calculate pressure */ + p=-(u_up-u_down)/dv; +printf("-------> %.10f %.10f %f\n",u_up/EV/moldyn->count,u_down/EV/moldyn->count,p/BAR); /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; /* restore energy */ - moldyn->energy=u; + potential_force_calc(moldyn); link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); - return sqrt(p); + return p; } double get_pressure(t_moldyn *moldyn) { @@ -907,12 +894,18 @@ double get_pressure(t_moldyn *moldyn) { } -int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) { +int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) { t_3dvec *dim; dim=&(moldyn->dim); + if(dir==SCALE_UP) + scale=1.0+scale; + + if(dir==SCALE_DOWN) + scale=1.0-scale; + if(x) dim->x*=scale; if(y) dim->y*=scale; if(z) dim->z*=scale; @@ -920,11 +913,17 @@ int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) { return 0; } -int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) { +int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) { int i; t_3dvec *r; + if(dir==SCALE_UP) + scale=1.0+scale; + + if(dir==SCALE_DOWN) + scale=1.0-scale; + for(i=0;icount;i++) { r=&(moldyn->atom[i].r); if(x) r->x*=scale; @@ -956,8 +955,8 @@ int scale_volume(t_moldyn *moldyn) { moldyn->debug=scale; /* scale the atoms and dimensions */ - scale_atoms(moldyn,scale,TRUE,TRUE,TRUE); - scale_dim(moldyn,scale,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); + scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); /* visualize dimensions */ if(vdim->x!=0) { @@ -1306,7 +1305,7 @@ return 0; temperature_calc(moldyn); pressure_calc(moldyn); //tp=thermodynamic_pressure_calc(moldyn); -//printf("thermodynamic p: %f %f %f - %f\n",moldyn->tp.x/BAR,moldyn->tp.y/BAR,moldyn->tp.z/BAR,tp/BAR); +//printf("thermodynamic p: %f\n",thermodynamic_pressure_calc(moldyn)/BAR); /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) @@ -1498,7 +1497,8 @@ int potential_force_calc(t_moldyn *moldyn) { /* single particle potential/force */ if(itom[i].attr&ATOM_ATTR_1BP) - moldyn->func1b(moldyn,&(itom[i])); + if(moldyn->func1b) + moldyn->func1b(moldyn,&(itom[i])); if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP))) continue; diff --git a/moldyn.h b/moldyn.h index 86444f4..5160668 100644 --- a/moldyn.h +++ b/moldyn.h @@ -111,6 +111,9 @@ typedef struct s_moldyn { double gp_sum; /* sum over all gp */ double mean_gp; /* mean value of gp */ + double mean_v; /* mean of virial */ + double virial_sum; /* sum over all calculated virials */ + double p_ref; /* reference pressure */ double p; /* actual pressure (computed by virial) */ double p_sum; /* sum over all p */ @@ -234,6 +237,10 @@ typedef struct s_moldyn { #define VERBOSE 1 #define QUIET 0 +#define SCALE_UP 'u' +#define SCALE_DOWN 'd' +#define SCALE_DIRECT 'D' + /* * potential related phsical values / constants * @@ -242,8 +249,7 @@ typedef struct s_moldyn { #define ONE_THIRD (1.0/3.0) #define C 0x06 -//#define LC_C 3.567 /* A */ -#define LC_C 3.560 /* A */ +#define LC_C 3.567 /* A */ #define M_C 12.011 /* amu */ #define SI 0x0e @@ -254,7 +260,7 @@ typedef struct s_moldyn { #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.7 /* A */ @@ -297,6 +303,8 @@ typedef struct s_moldyn { #define ALBE_D_SI 0.81472 #define ALBE_H_SI 0.259 +#define LC_SI_ALBE 5.429 + #define ALBE_R_C (2.00-0.15) #define ALBE_S_C (2.00+0.15) #define ALBE_A_C (6.00*EV/1.167) @@ -309,6 +317,8 @@ typedef struct s_moldyn { #define ALBE_D_C 6.28433 #define ALBE_H_C 0.5556 +#define LC_C_ALBE 3.566 + #define ALBE_R_SIC (2.40-0.20) #define ALBE_S_SIC (2.40+0.10) #define ALBE_A_SIC (4.36*EV/0.847) @@ -321,10 +331,11 @@ typedef struct s_moldyn { #define ALBE_D_SIC 180.314 #define ALBE_H_SIC 0.68 -#define ALBE_CHI_SIC 1.0 +#define LC_SIC_ALBE 4.359 + /* - * lattice constants + * lattice types */ #define CUBIC 0x01 @@ -384,8 +395,8 @@ double pressure_calc(t_moldyn *moldyn); double thermodynamic_pressure_calc(t_moldyn *moldyn); double get_pressure(t_moldyn *moldyn); int scale_volume(t_moldyn *moldyn); -int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z); -int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z); +int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z); +int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z); double e_kin_calc(t_moldyn *moldyn); double get_total_energy(t_moldyn *moldyn); diff --git a/potentials/albe.c b/potentials/albe.c index c2acf2f..a79f58d 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -43,34 +43,6 @@ int albe_mult_complete_params(t_albe_mult_params *p) { return 0; } -/* albe 1 body part */ -int albe_mult_1bp(t_moldyn *moldyn,t_atom *ai) { - - int brand; - t_albe_mult_params *params; - t_albe_exchange *exchange; - - brand=ai->brand; - params=moldyn->pot_params; - exchange=&(params->exchange); - - /* - * simple: point constant parameters only depending on atom i to - * their right values - */ - - exchange->gamma_i=&(params->gamma[brand]); - exchange->c_i=&(params->c[brand]); - exchange->d_i=&(params->d[brand]); - exchange->h_i=&(params->h[brand]); - - exchange->ci2=params->c[brand]*params->c[brand]; - exchange->di2=params->d[brand]*params->d[brand]; - exchange->ci2di2=exchange->ci2/exchange->di2; - - return 0; -} - /* albe 3 body potential function (first ij loop) */ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { @@ -92,11 +64,12 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { */ brand=ai->brand; - - if(brand==aj->brand) + if(brand==aj->brand) { S2=params->S2[brand]; - else + } + else { S2=params->S2mixed; + } /* dist_ij, d_ij2 */ v3_sub(&dist_ij,&(aj->r),&(ai->r)); @@ -152,12 +125,25 @@ int albe_mult_3bp_k1(t_moldyn *moldyn, R=params->R[brand]; S=params->S[brand]; S2=params->S2[brand]; + /* albe needs i,k depending c,d,h and gamma values */ + exchange->gamma_i=&(params->gamma[brand]); + exchange->c_i=&(params->c[brand]); + exchange->d_i=&(params->d[brand]); + exchange->h_i=&(params->h[brand]); } else { R=params->Rmixed; S=params->Smixed; S2=params->S2mixed; + /* albe needs i,k depending c,d,h and gamma values */ + exchange->gamma_i=&(params->gamma_m); + exchange->c_i=&(params->c_mixed); + exchange->d_i=&(params->d_mixed); + exchange->h_i=&(params->h_mixed); } + exchange->ci2=*(exchange->c_i)**(exchange->c_i); + exchange->di2=*(exchange->d_i)**(exchange->d_i); + exchange->ci2di2=exchange->ci2/exchange->di2; /* dist_ik, d_ik2 */ v3_sub(&dist_ik,&(ak->r),&(ai->r)); @@ -185,11 +171,11 @@ int albe_mult_3bp_k1(t_moldyn *moldyn, cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); /* g_ijk */ - h_cos=*(exchange->h_i)-cos_theta; + h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism d2_h_cos2=exchange->di2+(h_cos*h_cos); frac=exchange->ci2/d2_h_cos2; g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac); - dg=-2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; + dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f.. /* zeta sum += f_c_ik * g_ijk */ if(d_ik<=R) { @@ -205,6 +191,14 @@ int albe_mult_3bp_k1(t_moldyn *moldyn, exchange->zeta_ij+=f_c_ik*g; } +#ifdef DEBUG + if((ai==&(moldyn->atom[0]))| + (aj==&(moldyn->atom[864]))| + (ak==&(moldyn->atom[1003]))) { + printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos); + } +#endif + /* store even more data for second k loop */ exchange->g[kcount]=g; exchange->dg[kcount]=dg; @@ -289,7 +283,7 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { } /* force contribution */ - scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*df_a)); + scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*f_a)); v3_scale(&force,&(exchange->dist_ij),scale); v3_add(&(ai->f),&(ai->f),&force); v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij @@ -304,12 +298,12 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r), f_c,b,f_a,f_r); + printf(" %f %f %f\n",exchange->zeta_ij,.0,.0); } #endif /* virial */ - if(ajdist_ij)); + virial_calc(ai,&force,&(exchange->dist_ij)); /* dzeta prefactor = - 0.5 f_c f_a db */ exchange->pre_dzeta=-0.5*f_a*f_c*db; @@ -435,8 +429,7 @@ int albe_mult_3bp_k2(t_moldyn *moldyn, /* virial */ //v3_scale(&force,&force,-1.0); - if(ajkcount++; diff --git a/potentials/albe.h b/potentials/albe.h index 35a0762..09b96b1 100644 --- a/potentials/albe.h +++ b/potentials/albe.h @@ -63,17 +63,19 @@ typedef struct s_albe_mult_params { double mu_m; /* mixed mu */ double gamma[2]; + double gamma_m; double c[2]; + double c_mixed; double d[2]; + double d_mixed; double h[2]; + double h_mixed; t_albe_exchange exchange; /* exchange between 2bp and 3bp calc */ } t_albe_mult_params; /* function prototypes */ int albe_mult_complete_params(t_albe_mult_params *p); -int albe_mult_1bp(t_moldyn *moldyn,t_atom *ai); -int albe_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int albe_mult_3bp_k1(t_moldyn *moldyn, t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); diff --git a/potentials/tersoff.c b/potentials/tersoff.c index da92ea5..6bfcf30 100644 --- a/potentials/tersoff.c +++ b/potentials/tersoff.c @@ -312,6 +312,14 @@ int tersoff_mult_3bp_k1(t_moldyn *moldyn, exchange->zeta_ij+=f_c_ik*g; } +#ifdef DEBUG + if((ai==&(moldyn->atom[0]))| + (aj==&(moldyn->atom[864]))| + (ak==&(moldyn->atom[1003]))) { + printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos); + } +#endif + /* store even more data for second k loop */ exchange->g[kcount]=g; exchange->dg[kcount]=dg; @@ -414,6 +422,9 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); if(aj==&(moldyn->atom[0])) printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); + printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r), + f_c,b,f_a,f_r); + printf(" %f %f %f\n",exchange->zeta_ij,.0,.0); } #endif diff --git a/sic.c b/sic.c index 73b2935..cea2644 100644 --- a/sic.c +++ b/sic.c @@ -82,8 +82,6 @@ int main(int argc,char **argv) { t_moldyn md; /* potential parameters */ - t_lj_params lj; - t_ho_params ho; t_tersoff_mult_params tp; t_albe_mult_params ap; @@ -102,56 +100,36 @@ int main(int argc,char **argv) { set_int_alg(&md,MOLDYN_INTEGRATE_VERLET); /* choose potential */ - set_potential1b(&md,tersoff_mult_1bp); -#ifdef TERSOFF_ORIG - set_potential3b_j1(&md,tersoff_mult_2bp); - set_potential3b_k1(&md,tersoff_mult_3bp); - set_potential3b_j2(&md,tersoff_mult_post_2bp); -#elif ALBE - set_potential1b(&md,albe_mult_1bp); +#ifdef ALBE set_potential3b_j1(&md,albe_mult_3bp_j1); set_potential3b_k1(&md,albe_mult_3bp_k1); set_potential3b_j2(&md,albe_mult_3bp_j2); set_potential3b_k2(&md,albe_mult_3bp_k2); #else + set_potential1b(&md,tersoff_mult_1bp); set_potential3b_j1(&md,tersoff_mult_3bp_j1); set_potential3b_k1(&md,tersoff_mult_3bp_k1); set_potential3b_j2(&md,tersoff_mult_3bp_j2); set_potential3b_k2(&md,tersoff_mult_3bp_k2); #endif - //set_potential2b(&md,lennard_jones); - //set_potential2b(&md,harmonic_oscillator); #ifdef ALBE set_potential_params(&md,&ap); #else set_potential_params(&md,&tp); #endif - //set_potential_params(&md,&lj); - //set_potential_params(&md,&ho); /* cutoff radius */ +#ifdef ALBE set_cutoff(&md,ALBE_S_SI); - //set_cutoff(&md,TM_S_SI); - //set_cutoff(&md,LC_SI*sqrt(3.0)); - //set_cutoff(&md,2.0*LC_SI); +#else + set_cutoff(&md,TM_S_SI); +#endif /* * potential parameters */ - /* lennard jones */ - lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI*LJ_SIGMA_SI; - lj.sigma6*=lj.sigma6; - lj.sigma12=lj.sigma6*lj.sigma6; - lj.epsilon4=4.0*LJ_EPSILON_SI; - 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.spring_constant=LJ_EPSILON_SI; - /* * tersoff mult potential parameters for SiC */ @@ -217,36 +195,43 @@ int main(int argc,char **argv) { ap.r0_mixed=ALBE_R0_SIC; ap.lambda_m=ALBE_LAMBDA_SIC; ap.mu_m=ALBE_MU_SIC; + ap.gamma_m=ALBE_GAMMA_SIC; + ap.c_mixed=ALBE_C_SIC; + ap.d_mixed=ALBE_D_SIC; + ap.h_mixed=ALBE_H_SIC; albe_mult_complete_params(&ap); /* set (initial) dimensions of simulation volume */ - set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE); + //set_dim(&md,6*LC_SI_ALBE,6*LC_SI_ALBE,6*LC_SI_ALBE,TRUE); + //set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE); + //set_dim(&md,6*LC_C_ALBE,6*LC_C_ALBE,6*LC_C_ALBE,TRUE); //set_dim(&md,6*LC_C,6*LC_C,6*LC_C,TRUE); - //set_dim(&md,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,TRUE); + set_dim(&md,6*LC_SIC_ALBE,6*LC_SIC_ALBE,6*LC_SIC_ALBE,TRUE); /* set periodic boundary conditions in all directions */ set_pbc(&md,TRUE,TRUE,TRUE); /* 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, - //create_lattice(&md,DIAMOND,LC_C,C,M_C, - ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI, + //create_lattice(&md,DIAMOND,LC_C_ALBE,C,M_C, + //create_lattice(&md,DIAMOND,LC_SI_ALBE,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,NULL); + // 0,6,6,6,NULL); // 1,6,6,6,NULL); /* create centered zinc blende lattice */ - //r.x=0.5*0.25*TM_LC_3C_SIC; r.y=r.x; r.z=r.x; - //create_lattice(&md,FCC,TM_LC_3C_SIC,SI,M_SI, - // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, - // 0,6,6,6,&r); - //r.x+=0.25*TM_LC_3C_SIC; r.y=r.x; r.z=r.x; - //create_lattice(&md,FCC,TM_LC_3C_SIC,C,M_C, - // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, - // 1,6,6,6,&r); + /**/ + r.x=0.5*0.25*LC_SIC_ALBE; r.y=r.x; r.z=r.x; + create_lattice(&md,FCC,LC_SIC_ALBE,SI,M_SI, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + 0,6,6,6,&r); + r.x+=0.25*LC_SIC_ALBE; r.y=r.x; r.z=r.x; + create_lattice(&md,FCC,LC_SIC_ALBE,C,M_C, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + 1,6,6,6,&r); + /**/ moldyn_bc_check(&md); @@ -297,7 +282,7 @@ int main(int argc,char **argv) { /* create the simulation schedule */ /* initial configuration */ - moldyn_add_schedule(&md,100,1.0); + moldyn_add_schedule(&md,10000,1.0); /* adding atoms */ //for(inject=0;inject