From 47bc224440051c430cefe4b82a7fb320ba76eba4 Mon Sep 17 00:00:00 2001 From: hackbard Date: Thu, 23 Nov 2006 15:20:46 +0000 Subject: [PATCH] nearly finished tersoff potential ... --- moldyn.c | 496 +++++++++++++++++++------------------------------------ moldyn.h | 60 +++++-- 2 files changed, 223 insertions(+), 333 deletions(-) diff --git a/moldyn.c b/moldyn.c index 9f88211..830b8f0 100644 --- a/moldyn.c +++ b/moldyn.c @@ -903,380 +903,230 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { * tersoff potential & force for 2 sorts of atoms */ -/* tersoff 2 body part */ +/* tersoff 1 body part */ +int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) { + + int num; + t_tersoff_mult_params *params; + t_tersoff_exchange *exchange; + + num=ai->bnum; + params=moldyn->pot1b_params; + exchange=&(params->exchange); -int tersoff_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + /* + * simple: point constant parameters only depending on atom i to + * their right values + */ - t_tersoff_params *params; + exchange->beta=&(params->beta[num]); + exchange->n=&(params->n[num]); + exchange->c=&(params->c[num]); + exchange->d=&(params->d[num]); + exchange->h=&(params->h[num]); + + exchange->c2=params->c[num]*params->c[num]; + exchange->d2=params->d[num]*params->d[num]; + exchange->c2d2=exchange->c2/exchange->d2; + + return 0; +} + +/* tersoff 2 body part */ +int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_tersoff_mult_params *params; + t_tersoff_exchange *exchange; t_3dvec dist_ij; double d_ij; + double A,B,R,S,lambda; + int num; params=moldyn->pot_params; + num=ai->bnum; + exchange=&(params->exchange); + + exchange->run3bp=0; /* * we need: f_c, df_c, f_r, df_r * - * therefore we need: R, S + * therefore we need: R, S, A, lambda */ v3_sub(&dist_ij,&(ai->r),&(aj->r)); if(bc) check_per_bound(moldyn,&dist_ij); - if(ai->bnum==aj->bnum) { - S=params->S[ai->bnum]; - R=params->R[ai->bnum]; + /* save for use in 3bp */ /* REALLY ?!?!?! */ + exchange->dist_ij=dist_ij; + + /* constants */ + if(num==aj->bnum) { + S=params->S[num]; + R=params->R[num]; + A=params->A[num]; + lambda=params->lambda[num]; + /* more constants depending of atoms i and j, needed in 3bp */ + params->exchange.B=&(params->B[num]); + params->exchange.mu=params->mu[num]; + params->exchange.chi=1.0; } else { S=params->Smixed; R=params->Rmixed; + A=params->Amixed; + lambda=params->lambda_m; + /* more constants depending of atoms i and j, needed in 3bp */ + params->exchange.B=&(params->Bmixed); + params->exchange.mu=&(params->mu_m); + params->exchange.chi=params->chi; } d_ij=v3_norm(&dist_ij); - if(d_ij<=S) { - f_r=A*exp(-lamda*d_ij); - df_r=-lambda*f_r/d_ij; - if(d_ijf),&(ai->f),&force); - } - moldyn->energy+=(f_r*f_c); - } - - return 0; -} - -/* tersoff 3 body part */ - -int tersoff(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc,u8 bck) { - - t_tersoff_params *params; - t_3dvec dist_ij; - double d_ij; - - params=moldyn->pot_params; - - /* 2 body part of the tersoff potential */ - - v3_sub(&dist_ij,&(ai->r),&(aj->r)); - if(bc) check_per_bound(moldyn,&dist_ij); - d_ij=v3_norm(&dist_ij); - if(d_ij<=S) { - - /* determine the tersoff parameters */ - if(atom[i].element!=btom->element) { - S=sqrt(TERSOFF_S[e1]*TERSOFF_S[e2]); - R=R_m; - A=; - lambda=; - B=; - mu=; - chi=; - beta=; - betaN=; - - if(d_ij<=R) { - df_r=-lambda*A*exp(-lambda*d_ij)/d_ij; - v3_scale(&force,&dist_ij,df_r); - v3_add(&(atom[i].f),&(atom[i].f), - &force); - } - else { - s_r=S-R; - arg1=PI*(d_ij-R)/s_r; - f_c=0.5+0.5*cos(arg1); - df_c=-0.5*sin(arg1)*(PI/(s_r*d_ij)); - f_r=A*exp(-lambda*d_ij); - df_r=-lambda*f_r/d_ij; - scale=df_c*f_r+df_r*f_c; - v3_scale(&force,&dist_ij,scale); - v3_add(&(atom[i].f),&(atom[i].f), - &force); - } - } - else - continue; - - - /* end 2 body stuff */ - - /* determine cell neighbours of btom */ - ki=(btom->r.x+(moldyn->dim.x/2))/lc->x; - kj=(btom->r.y+(moldyn->dim.y/2))/lc->y; - kk=(btom->r.z+(moldyn->dim.z/2))/lc->z; - ck=link_cell_neighbour_index(moldyn,ki,kj,kk, - neighbourk); - - /* go for zeta - 3 body stuff! */ - zeta=0.0; - d_ij2=d_ij*d_ij; - - /* cell of btom */ - thisk=&(neighbourk[0]); - list_reset(thisk); - do { - ktom=thisk->current->data; - if(ktom==btom) - continue; - if(ktom==&(atom[i])) - continue; - - /* 3 body stuff (1) */ - - v3_sub(&dist_ik,ktom,&(atom[i])); - d_ik=v3_norm(&dist_ik); - if(d_ik<=Sik) { - - Rik=; - Sik=; - Aik=; - lambda_ik=; - Bik=; - mu_ik=; - omega_ik=; - c_i=; - d_i=; - h_i=; - - - if(d_ik<=Rik) { - f_cik=1.0; - df_cik=0.0; - } - else { - sik_rik=Sik-Rik; - arg1ik=PI*(d_ik-Rik)/sik_rik; - f_cik=0.5+0.5*cos(arg1ik); - df_cik=-0.5*sin(arg1ik)* \ - (PI/(sik_rik*d_ik)); - f_rik=Aik*exp(-lambda_ik*d_ik); - f_aik=-Bik*exp(-mu_ik*d_ik); - } - - v3_sub(&distance_jk,ktom,btom); - cos_theta=(d_ij2+d_ik*d_ik-d_jk*d_jk)/\ - (2*d_ij*d_ik); - sin_theta=sqrt(1.0/\ - (cos_theta*cos_theta)); - theta=arccos(cos_theta); - - - } - else - continue; - - /* end 3 body stuff (1) */ - + /* save for use in 3bp */ + exchange->d_ij=d_ij; - } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); - - /* direct neighbours of btom cell */ - for(k=1;kstart!=NULL) { - - do { - ktom=thisk->current->data; - if(ktom==&(atom[i])) - continue; + if(d_ij>S) + return 0; - /* 3 body stuff (2) */ + f_r=A*exp(-lamda*d_ij); + df_r=-lambda*f_r/d_ij; - } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); - - } - } - - /* indirect neighbours of btom cell */ - for(k=ck;k<27;k++) { - thisk=&(neighbourk[k]); - list_reset(thisk); - if(thisk->start!=NULL) { - - do { - ktom=thisk->current->data; - if(ktom==&(atom[i])) - continue; - - /* 3 body stuff */ - - } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); - - } - } + /* f_a, df_a calc + save for 3bp use */ + exchange->f_a=-B*exp(-mu*d_ij); + exchange->df_a=-mu*exchange->f_a/d_ij; + if(d_ijf),&(ai->f),&force); + /* energy is 0.5 f_r f_c, but we will sum it up twice ... */ + moldyn->energy+=(0.25*f_r*f_c); - /* - * direct neighbour cells of atom i - */ - for(j=1;jstart!=NULL) { + /* save for use in 3bp */ + exchange->f_c=f_c; + exchange->df_c=df_c; - do { - btom=this->current->data; + /* enable the run of 3bp function */ + exchange->run3bp=1; - /* 2 body stuff */ + return 0; +} +/* tersoff 3 body part */ - /* determine cell neighbours of btom */ - ki=(btom->r.x+(moldyn->dim.x/2))/lc->x; - kj=(btom->r.y+(moldyn->dim.y/2))/lc->y; - kk=(btom->r.z+(moldyn->dim.z/2))/lc->z; - ck=link_cell_neighbour_index(moldyn,ki,kj,kk, - neighbourk); +int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { - /* cell of btom */ - thisk=&(neighbourk[0]); - list_reset(thisk); - do { - ktom=thisk->current->data; - if(ktom==btom) - continue; - if(ktom==&(atom[i])) - continue; - - /* 3 body stuff (1) */ + t_tersoff_mult_params *params; + t_tersoff_exchange *exchange; + t_3dvec dist_ij,dist_ik,dist_jk; + t_3dvec db_ij,temp,force; + double R,S,s_r; + double d_ij,d_ik,d_jk; + double f_c,df_c,b_ij,f_a,df_a; + double B,mu; + int num; - } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + params=moldyn->pot_params; + num=ai->bnum; + exchange=params->exchange; - /* direct neighbours of btom cell */ - for(k=1;kstart!=NULL) { + if(!(exchange->run3bp)) + return 0; - do { - ktom=thisk->current->data; - if(ktom==&(atom[i])) - continue; + /* + * we need: f_c, d_fc, b_ij, db_ij, f_a, df_a + * + * we got f_c, df_c, f_a, df_a from 2bp calculation + */ - /* 3 body stuff (2) */ + d_ij=exchange->d_ij; + d_ij2=exchange->d_ij2; - } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + B=*(params->exchange.B); + mu=*(params->exchange.mu); - } - } + f_a=params->exchange.f_a; + df_a=params->exchange.df_a; + + /* d_ij is <= S, as we didn't return so far! */ - /* indirect neighbours of btom cell */ - for(k=ck;k<27;k++) { - thisk=&(neighbourk[k]); - list_reset(thisk); - if(thisk->start!=NULL) { + /* + * calc of b_ij (scalar) and db_ij (vector) + * + * - for b_ij: chi, beta, f_c_ik, w(=1), c, d, h, n, cos_theta + * + * - for db_ij: d_theta, sin_theta, cos_theta, f_c_ik, df_c_ik, + * w_ik, + * + */ - do { - ktom=thisk->current->data; - if(ktom==&(atom[i])) - continue; + + v3_sub(&dist_ik,&(aj->i),&(ak->r)); + if(bc) check_per_bound(moldyn,&dist_ik); + d_ik=v3_norm(&dist_ik); + + /* constants for f_c_ik calc */ + if(num==ak->bnum) { + R=params->R[num]; + S=params->S[num]; + } + else { + R=params->Rmixed; + S=params->Smixed; + } - /* 3 body stuff (3) */ + /* calc of f_c_ik */ + if(d_ik>S) + return 0; - } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + if(d_ikr),&(ak->r)); + if(bc) check_per_bound(moldyn,&dist_jk); + d_jk=v3_norm(&dist_jk); - } - } + // GO ON HERE !!! - } while(list_next(this)!=L_NO_NEXT_ELEMENT); + cos_theta=(d_ij2+d_ik*d_ik-d_jk*d_jk)/(2*d_ij*d_ik); + sin_theta=sqrt(1.0-(cos_theta*cos_theta)); + theta=arccos(cos_theta); - } - } + h_cos=(h-cos_theta); + h_cos2=h_cos*h_cos; + d2_h_cos2=d2-h_cos2; - /* - * indirect neighbour cells of atom i - */ - for(j=c;j<27;j++) { - this=&(neighbour[j]); - list_reset(this); - if(this->start!=NULL) { - - do { - btom=this->current->data; - - /* 2 body stuff */ - - - /* determine cell neighbours of btom */ - ki=(btom->r.x+(moldyn->dim.x/2))/lc->x; - kj=(btom->r.y+(moldyn->dim.y/2))/lc->y; - kk=(btom->r.z+(moldyn->dim.z/2))/lc->z; - ck=link_cell_neighbour_index(moldyn,ki,kj,kk, - neighbourk); - - /* cell of btom */ - thisk=&(neighbourk[0]); - list_reset(thisk); - do { - ktom=thisk->current->data; - if(ktom==btom) - continue; - if(ktom==&(atom[i])) - continue; + /* add forces */ + v3_add(&(ai->f),&(ai->f),&force); + /* energy is 0.5 f_r f_c, but we will sum it up twice ... */ + moldyn->energy+=(0.25*f_a*b_ij*f_c); - /* 3 body stuff (1) */ - - } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); - - /* direct neighbours of btom cell */ - for(k=1;kstart!=NULL) { - - do { - ktom=thisk->current->data; - if(ktom==&(atom[i])) - continue; - - /* 3 body stuff (2) */ - - } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); - - } - } - - /* indirect neighbours of btom cell */ - for(k=ck;k<27;k++) { - thisk=&(neighbourk[k]); - list_reset(thisk); - if(thisk->start!=NULL) { - - do { - ktom=thisk->current->data; - if(ktom==&(atom[i])) - continue; - - /* 3 body stuff (3) */ - - } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); - - } - } - - - } while(list_next(this)!=L_NO_NEXT_ELEMENT); - - } - } - - } - - moldyn->energy=0.5*u; - return 0; } diff --git a/moldyn.h b/moldyn.h index d98b428..9c63307 100644 --- a/moldyn.h +++ b/moldyn.h @@ -122,18 +122,58 @@ typedef struct s_lj_params { double epsilon4; } t_lj_params; -typedef struct s_tersoff_params { +/* tersoff exchange structure to exchange 2bp and 3bp calculated values */ + +typedef struct s_tersoff_exchange { + double f_c,df_c; + + t_3dvec dist_ij; + double d_ij; + double d_ij2; + + double chi; + + double *B; + double *mu; + + double *beta; + double *n; + double *c; + double *d; + double *h; + + double c2; + double d2; + double c2d2; + + u8 run3bp; +} t_tersoff_exchange; + +/* tersoff multi (2!) potential parameters */ + +typedef struct s_tersoff_mult_params { double S[2]; /* tersoff cutoff radii */ double R[2]; /* tersoff cutoff radii */ - double Smixed /* mixed S radius */ - double Rmixed /* mixed R radius */ - - double l_1,l_2; - double m_1,m_2; - double a_1,a_2; - double b_1,b_2; - double r_1,r_2; - double s_1,s_2; + double Smixed; /* mixed S radius */ + double Rmixed; /* mixed R radius */ + double A[2]; /* factor of tersoff attractive part */ + double B[2]; /* factor of tersoff repulsive part */ + double Amixed; /* mixed A factor */ + double Bmixed; /* mixed B factor */ + double lambda[2]; /* tersoff lambda */ + double lambda_m; /* mixed lambda */ + double mu[2]; /* tersoff mu */ + double mu_m; /* mixed mu */ + + double chi; + + double beta[2]; + double n[2]; + double c[2]; + double d[2]; + double h[2]; + + t_tersoff_exchange exchange; /* exchange between 2bp and 3bp calc */ } t_tersoff_params; /* -- 2.39.2