]> hackdaworld.org Git - physik/posic.git/commitdiff
first init of new albe implementation (still warnings!)
authorhackbard <hackbard@sage.physik.uni-augsburg.de>
Thu, 10 Jul 2008 20:13:33 +0000 (22:13 +0200)
committerhackbard <hackbard@sage.physik.uni-augsburg.de>
Thu, 10 Jul 2008 20:13:33 +0000 (22:13 +0200)
potentials/albe.c
potentials/albe.h

index c01ab15106d51834980fed6f9c0836c9336efcb1..9dfba13ed9993857b30dfef146e422147c6a44d7 100644 (file)
@@ -113,26 +113,47 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
        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) {
+/* first i loop */
+int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) {
 
+       int i;
        t_albe_mult_params *params;
        t_albe_exchange *exchange;
-       unsigned char brand;
-       double S2;
-       t_3dvec dist_ij;
-       double d_ij2,d_ij;
-
+       
        params=moldyn->pot_params;
        exchange=&(params->exchange);
 
-       /* reset zeta sum */
-       exchange->zeta_ij=0.0;
+       /* zero exchange values */
+       memset(exchange->dist,0,ALBE_MAXN*sizeof(t_3dvec));
+       memset(exchange->d2,0,ALBE_MAXN*sizeof(double));
+       memset(exchange->d,0,ALBE_MAXN*sizeof(double));
+       memset(exchange->zeta,0,ALBE_MAXN*sizeof(double));
+       for(i=0;i<ALBE_MAXN;i++)
+               memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(double));
+       memset(exchange->skip,0,ALBE_MAXN*sizeof(u8));
+       exchange->jcnt=0;
+       exchange->j2cnt=0;
+
+       return 0;
+}
+
+/* first j loop within first i loop */
+int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+       t_albe_mult_params *params;
+       t_albe_exchange *exchange;
+       double S2,S,R,d2,d,s_r,arg;
+       t_3dvec dist;
+       int j;
+       u8 brand;
+       
+       params=moldyn->pot_params;
+       exchange=&(params->exchange);
 
-       /*
-        * set ij depending values
-        */
+       /* get j counter */
+       j=exchange->jcnt;
 
+       /* set ij depending values */
        brand=ai->brand;
        if(brand==aj->brand) {
                S2=params->S2[brand];
@@ -142,161 +163,207 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        }
 
        /* dist_ij, d_ij2 */
-       v3_sub(&dist_ij,&(aj->r),&(ai->r));
-       if(bc) check_per_bound(moldyn,&dist_ij);
-       d_ij2=v3_absolute_square(&dist_ij);
+       v3_sub(&dist,&(aj->r),&(ai->r));
+       exchange->dist[j]=dist;
+       if(bc) check_per_bound(moldyn,&dist);
+       d2=v3_absolute_square(&dist);
+       exchange->d2[j]=d2;
 
        /* if d_ij2 > S2 => no force & potential energy contribution */
-       if(d_ij2>S2) {
+       if(d2>S2) {
                moldyn->run3bp=0;
+               exchange->skip[j]=1;
+               exchange->jcnt+=1;
                return 0;
        }
 
+       /* more ij depending values */
+       if(brand==aj->brand) {
+               R=params->R[brand];
+               S=params->S[brand];
+               /* set albe needs i,(j/k) depending c,d,h and gamma values */
+               exchange->gamma_[j]=&(params->gamma[brand]);
+               exchange->c_[j]=&(params->c[brand]);
+               exchange->d_[j]=&(params->d[brand]);
+               exchange->h_[j]=&(params->h[brand]);
+       }
+       else {
+               R=params->Rmixed;
+               S=params->Smixed;
+               /* albe needs i,(j/k) depending c,d,h and gamma values */
+               exchange->gamma_[j]=&(params->gamma_m);
+               exchange->c_[j]=&(params->c_mixed);
+               exchange->d_[j]=&(params->d_mixed);
+               exchange->h_[j]=&(params->h_mixed);
+       }
+       exchange->c2_[j]=*exchange->c_[j]**exchange->c_[j];
+       exchange->d2_[j]=*exchange->d_[j]**exchange->d_[j];
+       exchange->c2d2_[j]=exchange->c2_[j]/exchange->d2_[j];
+
        /* d_ij */
-       d_ij=sqrt(d_ij2);
+       d=sqrt(exchange->d2[j]);
+       exchange->d[j]=d;
+       
+       /* f_c, df_c */
+       if(d<R) {
+               exchange->f_c[j]=1.0;
+               exchange->df_c[j]=0.0;
+       }
+       else {
+               s_r=S-R;
+               arg=M_PI*(d-R)/s_r;
+               exchange->f_c[j]=0.5+0.5*cos(arg);
+               exchange->df_c[j]=0.5*sin(arg)*(M_PI/(s_r*d));
+       }
 
-       /* store values */
-       exchange->dist_ij=dist_ij;
-       exchange->d_ij2=d_ij2;
-       exchange->d_ij=d_ij;
+       /* reset k counter */
+       exchange->kcnt=0;
 
-       /* reset k counter for first k loop */
-       exchange->kcount=0;
-               
        return 0;
 }
 
-/* albe 3 body potential function (first k loop) */
-int albe_mult_3bp_k1(t_moldyn *moldyn,
-                     t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+/* first k loop within first j loop within first i loop */
+int albe_mult_i0_j0_k0(t_moldyn *moldyn,
+                       t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
        t_albe_mult_params *params;
        t_albe_exchange *exchange;
-       unsigned char brand;
-       double R,S,S2;
-       t_3dvec dist_ij,dist_ik;
-       double d_ik2,d_ik,d_ij;
-       double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
-       double f_c_ik,df_c_ik;
-       int kcount;
+       int j,k;
+       t_3dvec distj,distk;
+       double dj,dk,djdk_inv,cos_theta;
+       double gj,dgj,h_cos_j,d2_h_cos2_j,frac_j;
+       double gk,dgk,h_cos_k,d2_h_cos2_k,frac_k;
+       t_3dvec dcosdrj,dcosdrk,tmp,**dzeta;
 
        params=moldyn->pot_params;
        exchange=&(params->exchange);
-       kcount=exchange->kcount;
 
-       if(kcount>ALBE_MAXN) {
-               printf("FATAL: neighbours = %d\n",kcount);
-               printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
-       }
-
-       /* ik constants */
-       brand=ai->brand;
-       if(brand==ak->brand) {
-               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]);
+       /* k<j & check whether to run k */
+       j=exchange->jcnt;
+       k=exchange->kcnt;
+       if(k>=ALBE_MAXN) {
+               printf("FATAL: too many neighbours! (%d)\n",k);
+               printf("  atom i:%d | j:%d | k:%d\n",ai->tag,aj->tag,ak->tag);
        }
-       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));
-       if(bc) check_per_bound(moldyn,&dist_ik);
-       d_ik2=v3_absolute_square(&dist_ik);
-
-       /* store data for second k loop */
-       exchange->dist_ik[kcount]=dist_ik;
-       exchange->d_ik2[kcount]=d_ik2;
-
-       /* return if not within cutoff */
-       if(d_ik2>S2) {
-               exchange->kcount++;
+       if((k>=j)|(exchange->skip[k])) {
+               exchange->kcnt+=1;
                return 0;
        }
-
-       /* d_ik */
-       d_ik=sqrt(d_ik2);
-
-       /* dist_ij, d_ij */
-       dist_ij=exchange->dist_ij;
-       d_ij=exchange->d_ij;
+       
+       /* distances */
+       distj=exchange->dist[j];
+       distk=exchange->dist[k];
+       dj=exchange->d[j];
+       dk=exchange->d[k];
+       djdk_inv=1.0/(dj*dk);
 
        /* cos theta */
-       cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
-
-       /* g_ijk */
-       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; // + in albe f..
-
-       /* zeta sum += f_c_ik * g_ijk */
-       if(d_ik<=R) {
-               exchange->zeta_ij+=g;
-               f_c_ik=1.0;
-               df_c_ik=0.0;
+       cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
+
+       /* g(cos(theta)) - in albe: ik-depending values! */
+       h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism
+       d2_h_cos2_j=exchange->d2_[j]+(h_cos_j*h_cos_j);
+       frac_j=exchange->c2_[j]/d2_h_cos2_j;
+       gj=1.0+exchange->c2d2_[j]-frac_j;
+       gj*=*(exchange->gamma_[j]);
+       dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe 
+       if(ak->brand==aj->brand) {
+               gk=gj;
+               dgk=dgj;
        }
        else {
-               s_r=S-R;
-               arg=M_PI*(d_ik-R)/s_r;
-               f_c_ik=0.5+0.5*cos(arg);
-               df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
-               exchange->zeta_ij+=f_c_ik*g;
+               h_cos_k=*(exchange->h_[k])+cos_theta;
+               d2_h_cos2_k=exchange->d2_[k]+(h_cos_k*h_cos_k);
+               frac_k=exchange->c2_[k]/d2_h_cos2_k;
+               gk=1.0+exchange->c2d2_[k]-frac_k;
+               gk*=*(exchange->gamma_[k]);
+               dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
        }
 
-       /* store even more data for second k loop */
-       exchange->g[kcount]=g;
-       exchange->dg[kcount]=dg;
-       exchange->d_ik[kcount]=d_ik;
-       exchange->cos_theta[kcount]=cos_theta;
-       exchange->f_c_ik[kcount]=f_c_ik;
-       exchange->df_c_ik[kcount]=df_c_ik;
+       /* zeta */
+       exchange->zeta[j]+=(exchange->f_c[k]*gk);
+       exchange->zeta[k]+=(exchange->f_c[j]*gj);
+
+       /* cos theta derivatives */
+       v3_scale(&dcosdrj,&distk,djdk_inv);             // j
+       v3_scale(&tmp,&distj,-cos_theta/exchange->d2[j]);
+       v3_add(&dcosdrj,&dcosdrj,&tmp);
+       v3_scale(&dcosdrk,&distj,djdk_inv);             // k
+       v3_scale(&tmp,&distk,-cos_theta/exchange->d2[k]);
+       v3_add(&dcosdrk,&dcosdrk,&tmp);
+
+       /* zeta derivatives */
+       dzeta=exchange->dzeta;
+       v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk);
+       v3_add(&dzeta[j][j],&dzeta[j][j],&tmp);         // j j
+       v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj);
+       v3_add(&dzeta[k][k],&dzeta[k][k],&tmp);         // k k
+       v3_scale(&tmp,&distk,exchange->df_c[k]*gk/dk);
+       v3_add(&dzeta[j][k],&dzeta[j][k],&tmp);
+       v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk);
+       v3_add(&dzeta[j][k],&dzeta[j][k],&tmp);         // j k
+       v3_scale(&tmp,&distj,exchange->df_c[j]*gj/dj);
+       v3_add(&dzeta[k][j],&dzeta[k][j],&tmp);
+       v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj);
+       v3_add(&dzeta[k][j],&dzeta[k][j],&tmp);         // k j
 
        /* increase k counter */
-       exchange->kcount++;
+       exchange->kcnt+=1;
+               
+       return 0;
+}
+
+/* first j loop within first i loop */
+int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+       t_albe_mult_params *params;
+       t_albe_exchange *exchange;
+
+       params=moldyn->pot_params;
+       exchange=&(params->exchange);
+
+       /* increase j counter */
+       exchange->jcnt+=1;
 
        return 0;
 }
 
-int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+/* second j loop within first i loop */
+int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        t_albe_mult_params *params;
        t_albe_exchange *exchange;
-       t_3dvec force;
-       double f_a,df_a,b,db,f_c,df_c;
-       double f_r,df_r;
-       double scale;
-       double mu,B;
-       double lambda,A;
-       double d_ij,r0;
-       unsigned char brand;
-       double S,R,s_r,arg;
+
+       int j;
+       double d,f_a,df_a,f_r,df_r,f_c,df_c,b,db;
+       double A,B,mu,lambda,r0;
        double energy;
+       t_3dvec *dist,force;
+       double scale;
+       u8 brand;
 
        params=moldyn->pot_params;
        exchange=&(params->exchange);
 
+       /* get j counter */
+       j=exchange->j2cnt;
+
+       /* skip if j not within cutoff */
+       if(exchange->skip[j]) {
+               moldyn->run3bp=0;
+               exchange->j2cnt+=1;
+               return 0;
+       }
+
+       /* distance */
+       d=exchange->d[j];
+       dist=&(exchange->dist[j]);
+       f_c=exchange->f_c[j];
+       df_c=exchange->df_c[j];
+
+       /* determine parameters to calculate fa, dfa, fr, dfr */
        brand=aj->brand;
        if(brand==ai->brand) {
-               S=params->S[brand];
-               R=params->R[brand];
                B=params->B[brand];
                A=params->A[brand];
                r0=params->r0[brand];
@@ -304,8 +371,6 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                lambda=params->lambda[brand];
        }
        else {
-               S=params->Smixed;
-               R=params->Rmixed;
                B=params->Bmixed;
                A=params->Amixed;
                r0=params->r0_mixed;
@@ -313,212 +378,101 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                lambda=params->lambda_m;
        }
 
-       d_ij=exchange->d_ij;
-
-       /* f_c, df_c */
-       if(d_ij<R) {
-               f_c=1.0;
-               df_c=0.0;
-       }
-       else {
-               s_r=S-R;
-               arg=M_PI*(d_ij-R)/s_r;
-               f_c=0.5+0.5*cos(arg);
-               df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
-       }
-
        /* f_a, df_a */
-       f_a=-B*exp(-mu*(d_ij-r0));
-       df_a=mu*f_a/d_ij;
+       f_a=-B*exp(-mu*(d-r0));
+       df_a=mu*f_a/d;
 
        /* f_r, df_r */
-       f_r=A*exp(-lambda*(d_ij-r0));
-       df_r=lambda*f_r/d_ij;
+       f_r=A*exp(-lambda*(d-r0));
+       df_r=lambda*f_r/d;
 
        /* b, db */
-       if(exchange->zeta_ij==0.0) {
-               b=1.0;
-               db=0.0;
-       }
-       else {
-               b=1.0/sqrt(1.0+exchange->zeta_ij);
-               db=-0.5*b/(1.0+exchange->zeta_ij);
-       }
+       b=1.0/sqrt(1.0+exchange->zeta[j]);
+       db=-0.5*b/(1.0+exchange->zeta[j]);
+
+       /* energy contribution */
+       energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
+       moldyn->energy+=energy;
+       ai->e+=energy;
 
-       /* force contribution for atom i */
+       /* force contribution for atom i due to ij bond */
        scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
-       v3_scale(&force,&(exchange->dist_ij),scale);
+       v3_scale(&force,dist,scale);
        v3_add(&(ai->f),&(ai->f),&force);
 
-       /* force contribution for atom j */
+       /* force contribution for atom j due to ij bond */
        v3_scale(&force,&force,-1.0); // dri rij = - drj rij
        v3_add(&(aj->f),&(aj->f),&force);
 
        /* virial */
-       virial_calc(aj,&force,&(exchange->dist_ij));
-
-#ifdef DEBUG
-if(moldyn->time>DSTART&&moldyn->time<DEND) {
-       if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
-               printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
-               printf("  adding %f %f %f\n",force.x,force.y,force.z);
-               if(ai==&(moldyn->atom[0]))
-                       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
+       virial_calc(aj,&force,dist);
 
        /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
-       exchange->pre_dzeta=0.5*f_a*f_c*db;
+       exchange->pre_dzeta=0.5*f_a*exchange->f_c[j]*db;
 
-       /* energy contribution */
-       energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
-       moldyn->energy+=energy;
-       ai->e+=energy;
+       /* force contribution (drj derivative) */
+       v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta);
+       v3_add(&(aj->f),&(aj->f),&force);
+       v3_scale(&force,&force,-1.0);
+       v3_add(&(ai->f),&(ai->f),&force);
+
+       /* virial */
+       virial_calc(ai,&force,dist);
 
        /* reset k counter for second k loop */
-       exchange->kcount=0;
+       exchange->kcnt=0;
                
        return 0;
 }
 
-/* albe 3 body potential function (second k loop) */
-int albe_mult_3bp_k2(t_moldyn *moldyn,
-                     t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+/* second k loop within second j loop within first i loop */
+int albe_mult_i0_j2_k0(t_moldyn *moldyn,
+                       t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
        t_albe_mult_params *params;
        t_albe_exchange *exchange;
-       int kcount;
-       t_3dvec dist_ik,dist_ij;
-       double d_ik2,d_ik,d_ij2,d_ij;
-       unsigned char brand;
-       double S2;
-       double g,dg,cos_theta;
-       double pre_dzeta;
-       double f_c_ik,df_c_ik;
-       double dijdik_inv,fcdg,dfcg;
-       t_3dvec dcosdrj,dcosdrk;
-       t_3dvec force,tmp;
+
+       int j,k;
+       t_3dvec force;
 
        params=moldyn->pot_params;
        exchange=&(params->exchange);
-       kcount=exchange->kcount;
-
-       if(kcount>ALBE_MAXN)
-               printf("FATAL: neighbours!\n");
-
-       /* d_ik2 */
-       d_ik2=exchange->d_ik2[kcount];
-
-       brand=ak->brand;
-       if(brand==ai->brand)
-               S2=params->S2[brand];
-       else
-               S2=params->S2mixed;
 
-       /* return if d_ik > S */
-       if(d_ik2>S2) {
-               exchange->kcount++;
+       /* k!=j & check whether to run k */
+       k=exchange->kcnt;
+       j=exchange->j2cnt;
+       if((k==j)|(exchange->skip[k])) {
+               exchange->kcnt+=1;
                return 0;
        }
-
-       /* prefactor dzeta */
-       pre_dzeta=exchange->pre_dzeta;
-
-       /* dist_ik, d_ik */
-       dist_ik=exchange->dist_ik[kcount];
-       d_ik=exchange->d_ik[kcount];
-
-       /* f_c_ik, df_c_ik */
-       f_c_ik=exchange->f_c_ik[kcount];
-       df_c_ik=exchange->df_c_ik[kcount];
-
-       /* dist_ij, d_ij2, d_ij */
-       dist_ij=exchange->dist_ij;
-       d_ij2=exchange->d_ij2;
-       d_ij=exchange->d_ij;
-
-       /* g, dg, cos_theta */
-       g=exchange->g[kcount];
-       dg=exchange->dg[kcount];
-       cos_theta=exchange->cos_theta[kcount];
-
-       /* cos_theta derivatives wrt j,k */
-       dijdik_inv=1.0/(d_ij*d_ik);
-       v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
-       v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
-       v3_add(&dcosdrj,&dcosdrj,&tmp);
-       v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
-       v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
-       v3_add(&dcosdrk,&dcosdrk,&tmp);
-
-       /* f_c_ik * dg, df_c_ik * g */
-       fcdg=f_c_ik*dg;
-       dfcg=df_c_ik*g;
-
-       /* derivative wrt j */
-       v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
-
-       /* force contribution */
-       v3_add(&(aj->f),&(aj->f),&force);
-
-#ifdef DEBUG
-if(moldyn->time>DSTART&&moldyn->time<DEND) {
-       if(aj==&(moldyn->atom[DATOM])) {
-               printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
-               printf("  adding %f %f %f\n",force.x,force.y,force.z);
-               printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
-               printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
-               printf("    d ij ik = %f %f\n",d_ij,d_ik);
-       }
-}
-#endif
-
-       /* force contribution to atom i */
+       
+       /* force contribution (drk derivative) */
+       v3_scale(&force,&(exchange->dzeta[j][k]),exchange->pre_dzeta);
+       v3_add(&(ak->f),&(ak->f),&force);
        v3_scale(&force,&force,-1.0);
        v3_add(&(ai->f),&(ai->f),&force);
 
        /* virial */
-       virial_calc(ai,&force,&dist_ij);
-
-       /* derivative wrt k */
-       v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
-       v3_scale(&tmp,&dcosdrk,fcdg);
-       v3_add(&force,&force,&tmp);
-       v3_scale(&force,&force,pre_dzeta);
+       virial_calc(ai,&force,&(exchange->dist[k]));
 
-       /* force contribution */
-       v3_add(&(ak->f),&(ak->f),&force);
+       /* increase k counter */
+       exchange->kcnt+=1;
 
-#ifdef DEBUG
-if(moldyn->time>DSTART&&moldyn->time<DEND) {
-       if(ak==&(moldyn->atom[DATOM])) {
-               printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
-               printf("  adding %f %f %f\n",force.x,force.y,force.z);
-               printf("  total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
-               printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
-               printf("    d ij ik = %f %f\n",d_ij,d_ik);
-       }
+       return 0;
 }
-#endif
 
-       /* force contribution to atom i */
-       v3_scale(&force,&force,-1.0);
-       v3_add(&(ai->f),&(ai->f),&force);
+int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
-       /* virial */
-       virial_calc(ai,&force,&dist_ik);
-       
-       /* increase k counter */
-       exchange->kcount++;     
+       t_albe_mult_params *params;
+       t_albe_exchange *exchange;
 
-       return 0;
+       params=moldyn->pot_params;
+       exchange=&(params->exchange);
 
+       /* increase j counter */
+       exchange->j2cnt+=1;
+
+       return 0;
 }
 
 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
index 15386e76e5392f997e26e7e2639108cc24a985a6..0e4d4eece731f44decced154116a1e594a769140 100644 (file)
 /* albe exchange type */
 typedef struct s_albe_exchange {
 
-       t_3dvec dist_ij;
-       double d_ij2;
-       double d_ij;
+       t_3dvec dist[ALBE_MAXN];
+       double d2[ALBE_MAXN];
+       double d[ALBE_MAXN];
 
-       t_3dvec dist_ik[ALBE_MAXN];
-       double d_ik2[ALBE_MAXN];
-       double d_ik[ALBE_MAXN];
+       double f_c[ALBE_MAXN];
+       double df_c[ALBE_MAXN];
 
-       double f_c_ik[ALBE_MAXN];
-       double df_c_ik[ALBE_MAXN];
+       double zeta[ALBE_MAXN];
+       t_3dvec dzeta[ALBE_MAXN][ALBE_MAXN];
 
-       double g[ALBE_MAXN];
-       double dg[ALBE_MAXN];
-       double cos_theta[ALBE_MAXN];
+       u8 skip[ALBE_MAXN];
 
-       double *gamma_i;
-       double *c_i;
-       double *d_i;
-       double *h_i;
+       double *gamma_[ALBE_MAXN];
+       double *c_[ALBE_MAXN];
+       double *d_[ALBE_MAXN];
+       double c2_[ALBE_MAXN];
+       double d2_[ALBE_MAXN];
+       double c2d2_[ALBE_MAXN];
+       double *h_[ALBE_MAXN];
 
-       double ci2;
-       double di2;
-       double ci2di2;
-
-       double zeta_ij;
        double pre_dzeta;
 
-       int kcount;
+       int jcnt;
+       int j2cnt;
+       int kcnt;
 } t_albe_exchange;
 
 /* albe mult (2!) potential parameters */
@@ -76,12 +73,15 @@ typedef struct s_albe_mult_params {
 
 /* function prototypes */
 int albe_mult_set_params(t_moldyn *moldyn,int element1,int elemnt2);
-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);
-int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
-int albe_mult_3bp_k2(t_moldyn *moldyn,
-                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
+int albe_mult_i0(t_moldyn *moldyn,t_atom *ai);
+int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int albe_mult_i0_j0_k0(t_moldyn *moldyn,
+                       t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
+int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int albe_mult_i0_j2_k0(t_moldyn *moldyn,
+                       t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
+int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc);
 
 /* albe potential parameter defines */