]> hackdaworld.org Git - physik/posic.git/commitdiff
nearly finished tersoff potential ...
authorhackbard <hackbard>
Thu, 23 Nov 2006 15:20:46 +0000 (15:20 +0000)
committerhackbard <hackbard>
Thu, 23 Nov 2006 15:20:46 +0000 (15:20 +0000)
moldyn.c
moldyn.h

index 9f8821164dd9ae3b0eb406ec21d8c3b5cbb12709..830b8f044cbab93357598b57f59180b17f4c5d45 100644 (file)
--- 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_ij<R) {
-                       /* f_c = 1, df_c = 0 */
-                       v3_scale(&force,&dist_ij,df_r);
-               }
-               else {
-                       s_r=S-R;
-                       arg=PI*(d_ij-R)/s_r;
-                       f_c=0.5+0.5*cos(arg);
-                       df_c=-0.5*sin(arg)*(PI/(s_r*d_ij));
-                       scale=df_c*f_r+df_r*f_c;
-                       v3_scale(&force,&dist_ij,scale);
-                       v3_add(&(ai->f),&(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;k<ck;k++) {
-                               thisk=&(neighbourk[k]);
-                               list_reset(thisk);
-                               if(thisk->start!=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_ij<R) {
+               /* f_c = 1, df_c = 0 */
+               f_c=1.0;
+               df_c=0.0;
+               v3_scale(&force,&dist_ij,df_r);
+       }
+       else {
+               s_r=S-R;
+               arg=PI*(d_ij-R)/s_r;
+               f_c=0.5+0.5*cos(arg);
+               df_c=-0.5*sin(arg)*(PI/(s_r*d_ij));
+               scale=df_c*f_r+df_r*f_c;
+               v3_scale(&force,&dist_ij,scale);
+       }
 
-               } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+       /* 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_r*f_c);
 
-               /*
-                * direct neighbour cells of atom i
-                */
-               for(j=1;j<c;j++) {
-                       this=&(neighbour[j]);
-                       list_reset(this);
-                       if(this->start!=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;k<ck;k++) {
-                               thisk=&(neighbourk[k]);
-                               list_reset(thisk);
-                               if(thisk->start!=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_ik<R) {
+               /* f_c_ik = 1, df_c_ik = 0 */
+               f_c_ik=1.0;
+               df_c_ik=0.0;
+       }
+       else {
+               s_r=S-R;
+               arg=PI*(d_ik-R)/s_r;
+               f_c_ik=0.5+0.5*cos(arg);
+               df_c_ik=-0.5*sin(arg)*(PI/(s_r*d_ik));
+       }
+       
+       v3_sub(&dist_jk,&(aj->r),&(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;k<ck;k++) {
-                               thisk=&(neighbourk[k]);
-                               list_reset(thisk);
-                               if(thisk->start!=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;
 }
 
index d98b4286268ce6dfaebace9bb1fb9d234d43630f..9c633076270e6a3df40aa2df6992e118a2114dc7 100644 (file)
--- 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;
 
 /*