]> hackdaworld.org Git - physik/posic.git/commitdiff
new albe potential + new force calc routine (old potentials need to be
authorhackbard <hackbard@sage.physik.uni-augsburg.de>
Thu, 24 Jul 2008 11:55:41 +0000 (13:55 +0200)
committerhackbard <hackbard@sage.physik.uni-augsburg.de>
Thu, 24 Jul 2008 11:55:41 +0000 (13:55 +0200)
adapted!)

moldyn.c
moldyn.h
potentials/albe.c
potentials/albe.h

index bf7b9f0cde519d988f962b6576a6c58d947142c8..1aff5a8df207c2b5d22077c70ccbc31c128afdd6 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -215,11 +215,11 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 
        switch(type) {
                case MOLDYN_POTENTIAL_TM:
-                       moldyn->func1b=tersoff_mult_1bp;
-                       moldyn->func3b_j1=tersoff_mult_3bp_j1;
-                       moldyn->func3b_k1=tersoff_mult_3bp_k1;
-                       moldyn->func3b_j2=tersoff_mult_3bp_j2;
-                       moldyn->func3b_k2=tersoff_mult_3bp_k2;
+                       moldyn->func_i0=tersoff_mult_1bp;
+                       moldyn->func_j1=tersoff_mult_3bp_j1;
+                       moldyn->func_j1_k0=tersoff_mult_3bp_k1;
+                       moldyn->func_j1c=tersoff_mult_3bp_j2;
+                       moldyn->func_j1_k1=tersoff_mult_3bp_k2;
                        moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
                        break;
                /*
@@ -232,21 +232,21 @@ int set_potential(t_moldyn *moldyn,u8 type) {
                        break;
                */
                case MOLDYN_POTENTIAL_AM:
-                       moldyn->func1b=albe_mult_i0;
-                       moldyn->func2b=albe_mult_i0_j0;
-                       moldyn->func3b_0=albe_mult_i0_j0_k0;
-                       moldyn->func3b_1=albe_mult_i0_j1;
-                       moldyn->func3b_j1=albe_mult_i0_j2;
-                       moldyn->func3b_k2=albe_mult_i0_j2_k0;
-                       moldyn->func3b_j2=albe_mult_i0_j3;
+                       moldyn->func_i0=albe_mult_i0;
+                       moldyn->func_j0=albe_mult_i0_j0;
+                       moldyn->func_j0_k0=albe_mult_i0_j0_k0;
+                       moldyn->func_j0e=albe_mult_i0_j1;
+                       moldyn->func_j1=albe_mult_i0_j2;
+                       moldyn->func_j1_k0=albe_mult_i0_j2_k0;
+                       moldyn->func_j1c=albe_mult_i0_j3;
                        moldyn->check_2b_bond=albe_mult_check_2b_bond;
                        break;
                case MOLDYN_POTENTIAL_HO:
-                       moldyn->func2b=harmonic_oscillator;
+                       moldyn->func_j0=harmonic_oscillator;
                        moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond;
                        break;
                case MOLDYN_POTENTIAL_LJ:
-                       moldyn->func2b=lennard_jones;
+                       moldyn->func_j0=lennard_jones;
                        moldyn->check_2b_bond=lennard_jones_check_2b_bond;
                        break;
                default:
@@ -1860,8 +1860,8 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                /* single particle potential/force */
                if(itom[i].attr&ATOM_ATTR_1BP)
-                       if(moldyn->func1b)
-                               moldyn->func1b(moldyn,&(itom[i]));
+                       if(moldyn->func_i0)
+                               moldyn->func_i0(moldyn,&(itom[i]));
 
                if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
                        continue;
@@ -1876,8 +1876,14 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                dnlc=lc->dnlc;
 
+#ifndef STATIC_LISTS
+               /* check for later 3 body interaction */
+               if(itom[i].attr&ATOM_ATTR_3BP)
+                       memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
+#endif
+
                /* first loop over atoms j */
-               if(moldyn->func2b) {
+               if(moldyn->func_j0) {
                        for(j=0;j<27;j++) {
 
                                bc_ij=(j<dnlc)?0:1;
@@ -1888,18 +1894,6 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                        jtom=&(atom[neighbour_i[j][p]]);
                                        p++;
-
-                                       if(jtom==&(itom[i]))
-                                               continue;
-
-                                       if((jtom->attr&ATOM_ATTR_2BP)&
-                                          (itom[i].attr&ATOM_ATTR_2BP)) {
-                                               moldyn->func2b(moldyn,
-                                                              &(itom[i]),
-                                                              jtom,
-                                                              bc_ij);
-                                       }
-                               }
 #else
                                this=&(neighbour_i[j]);
                                list_reset_f(this);
@@ -1909,6 +1903,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                do {
                                        jtom=this->current->data;
+#endif
 
                                        if(jtom==&(itom[i]))
                                                continue;
@@ -1918,13 +1913,12 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                        if((jtom->attr&ATOM_ATTR_2BP)&
                                           (itom[i].attr&ATOM_ATTR_2BP)) {
-                                               moldyn->func2b(moldyn,
-                                                              &(itom[i]),
-                                                              jtom,
-                                                              bc_ij);
+                                               moldyn->func_j0(moldyn,
+                                                               &(itom[i]),
+                                                               jtom,
+                                                               bc_ij);
                                        }
 
-// REWRITE ONCE WORKING!!!
                                        /* 3 body potential/force */
 
                                        /* in j loop, 3bp run can be skipped */
@@ -1937,9 +1931,10 @@ int potential_force_calc(t_moldyn *moldyn) {
                                        if(!(jtom->attr&ATOM_ATTR_3BP))
                                                continue;
 
-                                       if(moldyn->func3b_0==NULL)
+                                       if(moldyn->func_j0_k0==NULL)
                                                continue;
 
+                               /* first loop over atoms k in first j loop */
                                for(k=0;k<27;k++) {
 
                                        bc_ik=(k<dnlc)?0:1;
@@ -1964,17 +1959,17 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                if(!(ktom->attr&ATOM_ATTR_3BP))
                                                        continue;
 
-                                               if(ktom==jtom)
-                                                       continue;
+                                               //if(ktom==jtom)
+                                               //      continue;
 
                                                if(ktom==&(itom[i]))
                                                        continue;
 
-                                               moldyn->func3b_0(moldyn,
-                                                                &(itom[i]),
-                                                                jtom,
-                                                                ktom,
-                                                                bc_ik|bc_ij);
+                                               moldyn->func_j0_k0(moldyn,
+                                                                  &(itom[i]),
+                                                                  jtom,
+                                                                  ktom,
+                                                                  bc_ik|bc_ij);
 #ifdef STATIC_LISTS
                                        }
 #else
@@ -1984,35 +1979,27 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                }
 
-
-                               if(moldyn->func3b_1)
-                                       moldyn->func3b_1(moldyn,
+                               /* finish of first j loop after first k loop */
+                               if(moldyn->func_j0e)
+                                       moldyn->func_j0e(moldyn,
                                                         &(itom[i]),
                                                         jtom,
                                                         bc_ij);
 
-// UNTIL HERE + change function names!
-
-
-
+#ifdef STATIC_LISTS
+                               }
+#else
                                } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 #endif
 
                        }
                }
 
-               /* 3 body potential/force */
+               /* continued 3 body potential/force */
 
                if(!(itom[i].attr&ATOM_ATTR_3BP))
                        continue;
 
-               /* copy the neighbour lists */
-#ifdef STATIC_LISTS
-               /* no copy needed for static lists */
-#else
-               memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
-#endif
-
                /* second loop over atoms j */
                for(j=0;j<27;j++) {
 
@@ -2045,18 +2032,18 @@ int potential_force_calc(t_moldyn *moldyn) {
                                /* reset 3bp run */
                                moldyn->run3bp=1;
 
-                               if(moldyn->func3b_j1)
-                                       moldyn->func3b_j1(moldyn,
-                                                         &(itom[i]),
-                                                         jtom,
-                                                         bc_ij);
+                               if(moldyn->func_j1)
+                                       moldyn->func_j1(moldyn,
+                                                       &(itom[i]),
+                                                       jtom,
+                                                       bc_ij);
 
-                               /* in first j loop, 3bp run can be skipped */
+                               /* in j loop, 3bp run can be skipped */
                                if(!(moldyn->run3bp))
                                        continue;
                        
-                               /* first loop over atoms k */
-                               if(moldyn->func3b_k1) {
+                               /* first loop over atoms k in second j loop */
+                               if(moldyn->func_j1_k0) {
 
                                for(k=0;k<27;k++) {
 
@@ -2082,17 +2069,17 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                if(!(ktom->attr&ATOM_ATTR_3BP))
                                                        continue;
 
-                                               if(ktom==jtom)
-                                                       continue;
+                                               //if(ktom==jtom)
+                                               //      continue;
 
                                                if(ktom==&(itom[i]))
                                                        continue;
 
-                                               moldyn->func3b_k1(moldyn,
-                                                                 &(itom[i]),
-                                                                 jtom,
-                                                                 ktom,
-                                                                 bc_ik|bc_ij);
+                                               moldyn->func_j1_k0(moldyn,
+                                                                  &(itom[i]),
+                                                                  jtom,
+                                                                  ktom,
+                                                                  bc_ik|bc_ij);
 #ifdef STATIC_LISTS
                                        }
 #else
@@ -2104,14 +2091,15 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                }
 
-                               if(moldyn->func3b_j2)
-                                       moldyn->func3b_j2(moldyn,
-                                                         &(itom[i]),
-                                                         jtom,
-                                                         bc_ij);
+                               /* continued j loop after first k loop */
+                               if(moldyn->func_j1c)
+                                       moldyn->func_j1c(moldyn,
+                                                        &(itom[i]),
+                                                        jtom,
+                                                        bc_ij);
 
                                /* second loop over atoms k */
-                               if(moldyn->func3b_k2) {
+                               if(moldyn->func_j1_k1) {
 
                                for(k=0;k<27;k++) {
 
@@ -2137,17 +2125,17 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                if(!(ktom->attr&ATOM_ATTR_3BP))
                                                        continue;
 
-                                               if(ktom==jtom)
-                                                       continue;
+                                               //if(ktom==jtom)
+                                               //      continue;
 
                                                if(ktom==&(itom[i]))
                                                        continue;
 
-                                               moldyn->func3b_k2(moldyn,
-                                                                 &(itom[i]),
-                                                                 jtom,
-                                                                 ktom,
-                                                                 bc_ik|bc_ij);
+                                               moldyn->func_j1_k1(moldyn,
+                                                                  &(itom[i]),
+                                                                  jtom,
+                                                                  ktom,
+                                                                  bc_ik|bc_ij);
 
 #ifdef STATIC_LISTS
                                        }
@@ -2160,11 +2148,11 @@ int potential_force_calc(t_moldyn *moldyn) {
                                
                                }
 
-                               /* 2bp post function */
-                               if(moldyn->func3b_j3) {
-                                       moldyn->func3b_j3(moldyn,
-                                                         &(itom[i]),
-                                                         jtom,bc_ij);
+                               /* finish of j loop after second k loop */
+                               if(moldyn->func_j1e) {
+                                       moldyn->func_j1e(moldyn,
+                                                        &(itom[i]),
+                                                        jtom,bc_ij);
                                }
 #ifdef STATIC_LISTS
                        }
@@ -2173,7 +2161,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 #endif
                
                }
-               
+
 #ifdef DEBUG
        //printf("\n\n");
 #endif
index 6345f3faf1cfb75cce7394d05cdf67d670e6cc15..066ad3642cfa9768809c5c17f47001ebaa3aca71 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -104,18 +104,18 @@ typedef struct s_moldyn {
        double volume;          /* volume of sim cell (dim.x*dim.y*dim.z) */
 
        /* potential force function and parameter pointers */
-       int (*func1b)(struct s_moldyn *moldyn,t_atom *ai);
-       int (*func2b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
-       int (*func3b_0)(struct s_moldyn *moldyn,
-                       t_atom *ai,t_atom *aj,t_atom *ak,u8 bck);
-       int (*func3b_1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
-       int (*func3b_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
-       int (*func3b_j2)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
-       int (*func3b_j3)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
-       int (*func3b_k1)(struct s_moldyn *moldyn,
-                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bck);
-       int (*func3b_k2)(struct s_moldyn *moldyn,
-                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bck);
+       int (*func_i0)(struct s_moldyn *moldyn,t_atom *ai);
+       int (*func_j0)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+       int (*func_j0_k0)(struct s_moldyn *moldyn,
+                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bck);
+       int (*func_j0e)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+       int (*func_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+       int (*func_j1_k0)(struct s_moldyn *moldyn,
+                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bck);
+       int (*func_j1c)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+       int (*func_j1_k1)(struct s_moldyn *moldyn,
+                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bck);
+       int (*func_j1e)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
        void *pot_params;
        unsigned char run3bp;
 
index 9dfba13ed9993857b30dfef146e422147c6a44d7..e8aee74d4be12cf9abc3821acb5abb064b12fbc2 100644 (file)
@@ -53,7 +53,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
                        p->mu[0]=ALBE_MU_SI;
                        p->gamma[0]=ALBE_GAMMA_SI;
                        p->c[0]=ALBE_C_SI;
+                       p->c2[0]=p->c[0]*p->c[0];
                        p->d[0]=ALBE_D_SI;
+                       p->d2[0]=p->d[0]*p->d[0];
+                       p->c2d2[0]=p->c2[0]/p->d2[0];
                        p->h[0]=ALBE_H_SI;
                        switch(element2) {
                                case C:
@@ -67,7 +70,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
                                        p->mu[1]=ALBE_MU_C;
                                        p->gamma[1]=ALBE_GAMMA_C;
                                        p->c[1]=ALBE_C_C;
+                                       p->c2[1]=p->c[1]*p->c[1];
                                        p->d[1]=ALBE_D_C;
+                                       p->d2[1]=p->d[1]*p->d[1];
+                                       p->c2d2[1]=p->c2[1]/p->d2[1];
                                        p->h[1]=ALBE_H_C;
                                        /* mixed type: silicon carbide */
                                        p->Smixed=ALBE_S_SIC;
@@ -79,7 +85,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
                                        p->mu_m=ALBE_MU_SIC;
                                        p->gamma_m=ALBE_GAMMA_SIC;
                                        p->c_mixed=ALBE_C_SIC;
+                                       p->c2_mixed=p->c_mixed*p->c_mixed;
                                        p->d_mixed=ALBE_D_SIC;
+                                       p->d2_mixed=p->d_mixed*p->d_mixed;
+                                       p->c2d2_m=p->c2_mixed/p->d2_mixed;
                                        p->h_mixed=ALBE_H_SIC;
                                        break;
                                default:
@@ -105,10 +114,13 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
        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("  gamma  | %f | %f\n",p->gamma[0],p->gamma[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("  gamma  | %f | %f | %f\n",p->gamma[0],p->gamma[1],p->gamma_m);
+       printf("  c      | %f | %f | %f\n",p->c[0],p->c[1],p->c_mixed);
+       printf("  d      | %f | %f | %f\n",p->d[0],p->d[1],p->d_mixed);
+       printf("  c2     | %f | %f | %f\n",p->c2[0],p->c2[1],p->c2_mixed);
+       printf("  d2     | %f | %f | %f\n",p->d2[0],p->d2[1],p->d2_mixed);
+       printf("  c2d2   | %f | %f | %f\n",p->c2d2[0],p->c2d2[1],p->c2d2_m);
+       printf("  h      | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed);
 
        return 0;
 }
@@ -116,21 +128,18 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 /* first i loop */
 int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) {
 
-       int i;
        t_albe_mult_params *params;
        t_albe_exchange *exchange;
        
+       int i;
+
        params=moldyn->pot_params;
        exchange=&(params->exchange);
 
        /* 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));
+               memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(t_3dvec));
        exchange->jcnt=0;
        exchange->j2cnt=0;
 
@@ -142,6 +151,7 @@ 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;
@@ -164,8 +174,8 @@ int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        /* dist_ij, d_ij2 */
        v3_sub(&dist,&(aj->r),&(ai->r));
-       exchange->dist[j]=dist;
        if(bc) check_per_bound(moldyn,&dist);
+       exchange->dist[j]=dist;
        d2=v3_absolute_square(&dist);
        exchange->d2[j]=d2;
 
@@ -176,16 +186,20 @@ int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                exchange->jcnt+=1;
                return 0;
        }
+       exchange->skip[j]=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 */
+               /* 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]);
+               exchange->c2_[j]=&(params->c2[brand]);
+               exchange->d2_[j]=&(params->d2[brand]);
+               exchange->c2d2_[j]=&(params->c2d2[brand]);
        }
        else {
                R=params->Rmixed;
@@ -195,13 +209,13 @@ int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                exchange->c_[j]=&(params->c_mixed);
                exchange->d_[j]=&(params->d_mixed);
                exchange->h_[j]=&(params->h_mixed);
+               exchange->c2_[j]=&(params->c2_mixed);
+               exchange->d2_[j]=&(params->d2_mixed);
+               exchange->c2d2_[j]=&(params->c2d2_m);
        }
-       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=sqrt(exchange->d2[j]);
+       d=sqrt(d2);
        exchange->d[j]=d;
        
        /* f_c, df_c */
@@ -228,16 +242,23 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn,
 
        t_albe_mult_params *params;
        t_albe_exchange *exchange;
+
        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;
+       t_3dvec dcosdrj,dcosdrk,tmp;
+       t_3dvec *dzjj,*dzkk,*dzjk,*dzkj;
 
        params=moldyn->pot_params;
        exchange=&(params->exchange);
 
+       if(aj==ak) {
+               exchange->kcnt+=1;
+               return 0;
+       }
+
        /* k<j & check whether to run k */
        j=exchange->jcnt;
        k=exchange->kcnt;
@@ -249,7 +270,7 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn,
                exchange->kcnt+=1;
                return 0;
        }
-       
+
        /* distances */
        distj=exchange->dist[j];
        distk=exchange->dist[k];
@@ -260,11 +281,11 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn,
        /* cos theta */
        cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
 
-       /* g(cos(theta)) - in albe: ik-depending values! */
+       /* g(cos(theta)) ij and ik 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;
+       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) {
@@ -273,14 +294,19 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn,
        }
        else {
                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;
+               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;
        }
 
-       /* zeta */
+       /* zeta - for albe: ik depending g function */
+//if(ai->tag==0) {
+//     printf("------> %.15f %.15f\n",dj,dk);
+//     printf("------> %.15f %.15f\n",dj,dk);
+//}
+
        exchange->zeta[j]+=(exchange->f_c[k]*gk);
        exchange->zeta[k]+=(exchange->f_c[j]*gj);
 
@@ -293,19 +319,22 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn,
        v3_add(&dcosdrk,&dcosdrk,&tmp);
 
        /* zeta derivatives */
-       dzeta=exchange->dzeta;
+       dzjj=&(exchange->dzeta[j][j]);
+       dzkk=&(exchange->dzeta[k][k]);
+       dzjk=&(exchange->dzeta[j][k]);
+       dzkj=&(exchange->dzeta[k][j]);
        v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk);
-       v3_add(&dzeta[j][j],&dzeta[j][j],&tmp);         // j j
+       v3_add(dzjj,dzjj,&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_add(dzkk,dzkk,&tmp);                         // k k
+       v3_scale(&tmp,&distk,-exchange->df_c[k]*gk);    // dk rik = - di rik
+       v3_add(dzjk,dzjk,&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_add(dzjk,dzjk,&tmp);                         // j k
+       v3_scale(&tmp,&distj,-exchange->df_c[j]*gj);    // dj rij = - di rij
+       v3_add(dzkj,dzkj,&tmp);
        v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj);
-       v3_add(&dzeta[k][j],&dzeta[k][j],&tmp);         // k j
+       v3_add(dzkj,dzkj,&tmp);                         // k j
 
        /* increase k counter */
        exchange->kcnt+=1;
@@ -400,22 +429,51 @@ int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        v3_scale(&force,dist,scale);
        v3_add(&(ai->f),&(ai->f),&force);
 
+#ifdef NDEBUG
+if(ai->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d (ij) %.15f\n",force.x,force.y,force.z,ai->tag,aj->tag,exchange->zeta[j]);
+printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
+}
+#endif
+
        /* 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);
 
+#ifdef NDEBUG
+if(aj->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d (ji) %.15f\n",force.x,force.y,force.z,aj->tag,ai->tag,exchange->zeta[j]);
+printf("    t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
+}
+#endif
+
        /* virial */
        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*exchange->f_c[j]*db;
+       exchange->pre_dzeta=0.5*f_a*f_c*db;
 
        /* force contribution (drj derivative) */
        v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta);
        v3_add(&(aj->f),&(aj->f),&force);
+
+#ifdef NDEBUG
+if(aj->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d (j der)\n",force.x,force.y,force.z,aj->tag,ai->tag);
+printf("    t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
+}
+#endif
+
        v3_scale(&force,&force,-1.0);
        v3_add(&(ai->f),&(ai->f),&force);
 
+#ifdef NDEBUG
+if(ai->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d (i contr j der)\n",force.x,force.y,force.z,ai->tag,aj->tag);
+printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
+}
+#endif
+
        /* virial */
        virial_calc(ai,&force,dist);
 
@@ -438,6 +496,11 @@ int albe_mult_i0_j2_k0(t_moldyn *moldyn,
        params=moldyn->pot_params;
        exchange=&(params->exchange);
 
+       if(aj==ak) {
+               exchange->kcnt+=1;
+               return 0;
+       }
+
        /* k!=j & check whether to run k */
        k=exchange->kcnt;
        j=exchange->j2cnt;
@@ -449,9 +512,25 @@ int albe_mult_i0_j2_k0(t_moldyn *moldyn,
        /* force contribution (drk derivative) */
        v3_scale(&force,&(exchange->dzeta[j][k]),exchange->pre_dzeta);
        v3_add(&(ak->f),&(ak->f),&force);
+
+#ifdef NDEBUG
+if(ak->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d %d (k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag);
+printf("    t: %.15f %.15f %.15f\n",ak->f.x,ak->f.y,ak->f.z);
+}
+#endif
+
        v3_scale(&force,&force,-1.0);
        v3_add(&(ai->f),&(ai->f),&force);
 
+#ifdef NDEBUG
+if(ai->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d %d -- %d(i contr k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag,k);
+printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
+printf("  ## %f\n",exchange->d[k]);
+}
+#endif
+
        /* virial */
        virial_calc(ai,&force,&(exchange->dist[k]));
 
index 0e4d4eece731f44decced154116a1e594a769140..39c4df870261ad68b06139ec502ee3058e806c3f 100644 (file)
@@ -8,7 +8,7 @@
 #ifndef ALBE_H
 #define ALBE_H
 
-#define ALBE_MAXN      16*27
+#define ALBE_MAXN      (4*27)
 
 /* albe exchange type */
 typedef struct s_albe_exchange {
@@ -28,9 +28,9 @@ typedef struct s_albe_exchange {
        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 *c2_[ALBE_MAXN];
+       double *d2_[ALBE_MAXN];
+       double *c2d2_[ALBE_MAXN];
        double *h_[ALBE_MAXN];
 
        double pre_dzeta;
@@ -62,9 +62,15 @@ typedef struct s_albe_mult_params {
        double gamma[2];
        double gamma_m;
        double c[2];
+       double c2[2];
        double c_mixed;
+       double c2_mixed;
        double d[2];
+       double d2[2];
        double d_mixed;
+       double d2_mixed;
+       double c2d2[2];
+       double c2d2_m;
        double h[2];
        double h_mixed;