From dfbe75140f5a113f898809c529a882034507f6d1 Mon Sep 17 00:00:00 2001 From: hackbard Date: Mon, 4 Dec 2006 14:55:02 +0000 Subject: [PATCH] added 2bp post function .. this is getting sick! --- moldyn.c | 87 ++++++++++++++++++++++++++++++++++++++++++++++++++++---- moldyn.h | 16 ++++++++--- sic.c | 7 +++-- 3 files changed, 98 insertions(+), 12 deletions(-) diff --git a/moldyn.c b/moldyn.c index f504e68..52a716e 100644 --- a/moldyn.c +++ b/moldyn.c @@ -950,6 +950,13 @@ printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[ } } while(list_next(this)!=L_NO_NEXT_ELEMENT); + + /* 2bp post function */ + if(moldyn->func2b_post) + mlodyn->func2b_post(moldyn, + &(itom[i]), + jtom,bc_ij); + } } } @@ -1072,6 +1079,21 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) { p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]); p->mu_m=0.5*(p->mu[0]+p->mu[1]); + printf("[moldyn] tersoff mult parameter info:\n"); + printf(" S (m) | %.12f | %.12f | %.12f\n",p->S[0],p->S[1],p->Smixed); + printf(" R (m) | %.12f | %.12f | %.12f\n",p->R[0],p->R[1],p->Rmixed); + printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV); + printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV); + printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1], + p->lambda_m); + printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m); + printf(" beta | %.10f | %.10f\n",p->beta[0],p->beta[1]); + printf(" n | %f | %f\n",p->n[0],p->n[1]); + printf(" c | %f | %f\n",p->c[0],p->c[1]); + printf(" d | %f | %f\n",p->d[0],p->d[1]); + printf(" h | %f | %f\n",p->h[0],p->h[1]); + printf(" chi | %f \n",p->chi); + return 0; } @@ -1200,6 +1222,57 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* enable the run of 3bp function */ exchange->run3bp=1; + /* reset 3bp sums */ + exchange->3bp_sum1=0.0; + exchange->3bp_sum2=0.0; + + return 0; +} + +/* tersoff 2 body post part */ + +int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { + + /* here we have to allow for the 3bp sums */ + + t_tersoff_mult_params *params; + t_tersoff_exchange *exchange; + + t_3dvec force,temp,*db_ij; + double db_ij_scale1,db_ij_scale2; + double b_ij; + double f_c,df_c,f_a,df_a; + + params=moldyn->pot2b_params; + exchange=&(moldyn->exchange); + + db_ij=&(exchange->db_ij); + f_c=exchange->f_c; + df_c=exchange->df_c; + f_a=exchange->f_a; + df_a=exchange->df_a; + + db_ij_scale1=(1+betan*3bp_sum1); + db_ij_scale2=(n*betan*3bp_sum2); + help=pow(db_ij_scale1,-1.0/(2*n)-1); + b_ij=chi*db_ij_scale1*help; + db_ij_scale1=-chi/(2*n)*help; + + v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2)); + v3_scale(db_ij,db_ij,f_a); + + v3_scale(&temp,dist_ij,b_ij*df_a); + + v3_add(&force,&temp,db_ij); + v3_scale(&force,&force,f_c); + + v3_scale(&temp,&dist_ij,f_a*b_ij*df_c); + + /* add energy of 3bp sum */ + moldyn->energy+=(0.5*f_c*b_ij*f_a); + /* add force of 3bp calculation */ + v3_add(&(ai->f),&temp,&force); + return 0; } @@ -1242,9 +1315,13 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { d_ij=exchange->d_ij; d_ij2=exchange->d_ij2; + dist_ij=exchange->dist_ij; f_a=params->exchange.f_a; df_a=params->exchange.df_a; + + f_c=exchange->f_c; + df_c=exchange->df_c; /* d_ij is <= S, as we didn't return so far! */ @@ -1305,8 +1382,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { numer=d_ij2+d_ik*d_ik-d_jk*d_jk; denom=2*d_ij*d_ik; - //cos_theta=numer/denom; - cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); + cos_theta=numer/denom; + //cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); sin_theta=sqrt(1.0-(cos_theta*cos_theta)); theta=acos(cos_theta); d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom); @@ -1317,10 +1394,10 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { h_cos=(h-cos_theta); h_cos2=h_cos*h_cos; - d2_h_cos2=d2-h_cos2; + d2_h_cos2=d2+h_cos2; /* some usefull expressions */ - frac1=c2/(d2-h_cos2); + frac1=c2/(d2_h_cos2); bracket1=1+c2d2-frac1; if(f_c_ik==0.0) { bracket2=0.0; @@ -1337,8 +1414,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { printf("Foo -> 1: "); } bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0); -printf("THETA: %.15f %.15f\n",cos_theta,theta*180/(2*M_PI)); bracket3_pow=bracket3_pow_1*bracket3; +printf("%.15f %.15f %.15f\n",bracket2_n_1,bracket2_n); /* now go on with calc of b_ij and derivation of b_ij */ b_ij=chi*bracket3_pow; diff --git a/moldyn.h b/moldyn.h index 813af68..41b711a 100644 --- a/moldyn.h +++ b/moldyn.h @@ -36,9 +36,9 @@ typedef struct s_atom { #define ATOM_ATTR_FP 0x01 /* fixed position (bulk material) */ #define ATOM_ATTR_HB 0x02 /* coupled to heat bath (velocity scaling) */ -#define ATOM_ATTR_1BP 0x10 /* single paricle potential */ -#define ATOM_ATTR_2BP 0x20 /* pair potential */ -#define ATOM_ATTR_3BP 0x40 /* 3 body potential */ +#define ATOM_ATTR_1BP 0x10 /* single paricle potential */ +#define ATOM_ATTR_2BP 0x20 /* pair potential */ +#define ATOM_ATTR_3BP 0x40 /* 3 body potential */ /* cell lists */ typedef struct s_linkcell { @@ -73,6 +73,8 @@ typedef struct s_moldyn { int (*func1b)(struct s_moldyn *moldyn,t_atom *ai); void *pot1b_params; int (*func2b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + void (*func2b_post)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj, + u8 bc); void *pot2b_params; int (*func3b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak, u8 bck); @@ -194,6 +196,10 @@ typedef struct s_tersoff_exchange { double betan; u8 run3bp; + + t_3dvec db_ij; + double 3bp_sum1; + double 3bp_sum2; } t_tersoff_exchange; /* tersoff multi (2!) potential parameters */ @@ -285,7 +291,7 @@ typedef struct s_tersoff_mult_params { #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_D_SI 16.217 #define TM_H_SI (-0.59825) #define TM_R_C 1.8e-10 /* m */ @@ -311,6 +317,7 @@ typedef struct s_tersoff_mult_params { typedef int (*pf_func1b)(t_moldyn *,t_atom *ai); typedef int (*pf_func2b)(t_moldyn *,t_atom *,t_atom *,u8 bc); +typedef int (*pf_func2b_post)(t_moldyn *,t_atom *,t_atom *,u8 bc); typedef int (*pf_func3b)(t_moldyn *,t_atom *,t_atom *,t_atom *,u8 bc); int moldyn_init(t_moldyn *moldyn,int argc,char **argv); @@ -325,6 +332,7 @@ 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 moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer); diff --git a/sic.c b/sic.c index b10c681..3408895 100644 --- a/sic.c +++ b/sic.c @@ -83,6 +83,7 @@ int main(int argc,char **argv) { tp.n[1]=TM_N_C; tp.c[1]=TM_C_C; tp.d[1]=TM_D_C; + tp.h[1]=TM_H_C; tp.chi=TM_CHI_SIC; @@ -94,7 +95,7 @@ int main(int argc,char **argv) { /* set (initial) dimensions of simulation volume */ printf("[sic] setting dimensions\n"); - set_dim(&md,4*LC_SI,4*LC_SI,4*LC_SI,TRUE); + set_dim(&md,3*LC_SI,3*LC_SI,3*LC_SI,TRUE); /* set periodic boundary conditions in all directions */ printf("[sic] setting periodic boundary conditions\n"); @@ -105,7 +106,7 @@ int main(int argc,char **argv) { 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,4,4,4); + 0,3,3,3); /* setting a nearest neighbour distance for the moldyn checks */ set_nn_dist(&md,sqrt(3.0)*LC_SI/4.0); /* diamond ! */ @@ -124,7 +125,7 @@ int main(int argc,char **argv) { /* create the simulation schedule */ printf("[sic] adding schedule\n"); - moldyn_add_schedule(&md,100,1.0e-15); + moldyn_add_schedule(&md,10,1.0e-15); /* activate logging */ printf("[sic] activate logging\n"); -- 2.20.1