From acdd9a6aa72f3692ccd03cc0df20e3d60555f555 Mon Sep 17 00:00:00 2001 From: hackbard Date: Fri, 18 May 2007 11:54:07 +0000 Subject: [PATCH] first sic tests --- Makefile | 3 ++- moldyn.c | 32 +++++++++++++++++++--------- moldyn.h | 7 +++++- potentials/tersoff.c | 9 +++++--- potentials/tersoff_orig.c | 13 ++++++++++- sic.c | 45 ++++++++++++++++++++++++++++----------- 6 files changed, 81 insertions(+), 28 deletions(-) diff --git a/Makefile b/Makefile index 70fcb8b..ebad7e3 100644 --- a/Makefile +++ b/Makefile @@ -6,12 +6,13 @@ CFLAGS+=-ffloat-store CFLAGS+=-g #CFLAGS+=-DDEBUG #CFLAGS+=-DVDEBUG +#CFLAGS+=-DTERSOFF_ORIG LDFLAGS=-lm SOURCE=moldyn.c visual/visual.c random/random.c list/list.c -POT_SRC=potentials/tersoff.c #POT_SRC=potentials/tersoff_orig.c +POT_SRC=potentials/tersoff.c POT_SRC+= potentials/lennard_jones.c potentials/harmonic_oscillator.c all: sic diff --git a/moldyn.c b/moldyn.c index 705a28c..070ce69 100644 --- a/moldyn.c +++ b/moldyn.c @@ -421,11 +421,11 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { */ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, - u8 attr,u8 brand,int a,int b,int c) { + u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin) { int new,count; int ret; - t_3dvec origin; + t_3dvec orig; void *ptr; t_atom *atom; @@ -447,24 +447,33 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, atom=&(moldyn->atom[count]); /* no atoms on the boundaries (only reason: it looks better!) */ - origin.x=0.5*lc; - origin.y=0.5*lc; - origin.z=0.5*lc; + if(!origin) { + orig.x=0.5*lc; + orig.y=0.5*lc; + orig.z=0.5*lc; + } + else { + orig.x=origin->x; + orig.y=origin->y; + orig.z=origin->z; + } switch(type) { case CUBIC: set_nn_dist(moldyn,lc); - ret=cubic_init(a,b,c,lc,atom,&origin); + ret=cubic_init(a,b,c,lc,atom,&orig); break; case FCC: - v3_scale(&origin,&origin,0.5); + if(!origin) + v3_scale(&orig,&orig,0.5); set_nn_dist(moldyn,0.5*sqrt(2.0)*lc); - ret=fcc_init(a,b,c,lc,atom,&origin); + ret=fcc_init(a,b,c,lc,atom,&orig); break; case DIAMOND: - v3_scale(&origin,&origin,0.25); + if(!origin) + v3_scale(&orig,&orig,0.25); set_nn_dist(moldyn,0.25*sqrt(3.0)*lc); - ret=diamond_init(a,b,c,lc,atom,&origin); + ret=diamond_init(a,b,c,lc,atom,&orig); break; default: printf("unknown lattice type (%02x)\n",type); @@ -1252,6 +1261,9 @@ int moldyn_integrate(t_moldyn *moldyn) { /* calculate initial forces */ potential_force_calc(moldyn); +#ifdef DEBUG +return 0; +#endif /* some stupid checks before we actually start calculating bullshit */ if(moldyn->cutoff>0.5*moldyn->dim.x) diff --git a/moldyn.h b/moldyn.h index 5b12048..3bc40c0 100644 --- a/moldyn.h +++ b/moldyn.h @@ -242,12 +242,15 @@ typedef struct s_moldyn { #define ONE_THIRD (1.0/3.0) #define C 0x06 +#define LC_C (0.3567e-9*METER) /* A */ #define M_C 12.011 /* amu */ #define SI 0x0e #define LC_SI (0.543105e-9*METER) /* A */ #define M_SI 28.08553 /* amu */ +#define LC_3C_SIC (0.43596e-9*METER) /* 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 */ @@ -279,6 +282,8 @@ typedef struct s_moldyn { #define TM_CHI_SIC 0.9776 +#define TM_LC_3C_SIC (0.432e-9*METER) /* A */ + /* * lattice constants */ @@ -324,7 +329,7 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer); int moldyn_log_shutdown(t_moldyn *moldyn); int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, - u8 attr,u8 brand,int a,int b,int c); + u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin); int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin); int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin); int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin); diff --git a/potentials/tersoff.c b/potentials/tersoff.c index a4a87e3..da92ea5 100644 --- a/potentials/tersoff.c +++ b/potentials/tersoff.c @@ -418,7 +418,8 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { #endif /* virial */ - virial_calc(ai,&force,&(exchange->dist_ij)); + if(ajdist_ij)); /* dzeta prefactor = - 0.5 f_c f_a db */ exchange->pre_dzeta=-0.5*f_a*f_c*db; @@ -544,7 +545,8 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn, /* virial */ //v3_scale(&force,&force,-1.0); - virial_calc(ai,&force,&dist_ij); + if(ajkcount++; diff --git a/potentials/tersoff_orig.c b/potentials/tersoff_orig.c index 69b67ac..431cfeb 100644 --- a/potentials/tersoff_orig.c +++ b/potentials/tersoff_orig.c @@ -629,7 +629,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { v3_sub(&temp1,&temp1,&temp2); /* there is a minus! */ v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */ - /* store dg in temp2 and use it for dVjk later */ + /* store dg for use in dV_jk */ v3_copy(&temp2,&temp1); /* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */ @@ -670,8 +670,19 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { // dzeta_jk is simply dg, which is stored in temp2 } /* betajnj * zeta_jk ^ nj-1 */ +printf("FATAL db_jk calc!\n"); +printf("(z1 + z2 + z3 ...)^n != z1^n + z2^n + z3^n + ...\n"); +printf("st00pid me => tersoff_orig is obsolete!\n"); tmp=exchange->betajnj*pow(zeta,(n-1.0)); tmp=-chi/2.0*pow((1+tmp*zeta),(-1.0/(2.0*n)-1))*tmp; +#ifdef DEBUG + if((ai->tag==0)&(aj->tag==864)) { // &(ak->tag==23)) { + printf("\n\n"); + printf("db: ni zeta = %f %f\n", + n,zeta); + printf("\n\n"); + } +#endif v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5); v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */ /* scaled with 0.5 ^ */ diff --git a/sic.c b/sic.c index 2221025..3c815ba 100644 --- a/sic.c +++ b/sic.c @@ -13,8 +13,12 @@ /* potential */ #include "potentials/harmonic_oscillator.h" #include "potentials/lennard_jones.h" + +#ifdef TERSOFF_ORIG +#include "potentials/tersoff_orig.h" +#else #include "potentials/tersoff.h" -//#include "potentials/tersoff_orig.h" +#endif #define INJECT 1 #define NR_ATOMS 4 @@ -97,14 +101,16 @@ int main(int argc,char **argv) { /* choose potential */ set_potential1b(&md,tersoff_mult_1bp); - //set_potential2b(&md,tersoff_mult_2bp); - //set_potential3b_j1(&md,tersoff_mult_2bp); - //set_potential3b_k1(&md,tersoff_mult_3bp); - //set_potential3b_j3(&md,tersoff_mult_post_2bp); +#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); +#else 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); set_potential_params(&md,&tp); @@ -164,7 +170,8 @@ int main(int argc,char **argv) { tersoff_mult_complete_params(&tp); /* 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,6*LC_SI,6*LC_SI,TRUE); + set_dim(&md,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,TRUE); /* set periodic boundary conditions in all directions */ set_pbc(&md,TRUE,TRUE,TRUE); @@ -172,10 +179,21 @@ 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, + //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, // ATOM_ATTR_2BP|ATOM_ATTR_HB, - 0,6,6,6); + // 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); moldyn_bc_check(&md); /* testing configuration */ @@ -225,7 +243,7 @@ int main(int argc,char **argv) { /* create the simulation schedule */ /* initial configuration */ - moldyn_add_schedule(&md,10000,1.0); + moldyn_add_schedule(&md,1000,1.0); /* adding atoms */ //for(inject=0;inject