Merge branch 'leadoff'

[physik/posic.git] / potentials / albe.c

diff --git a/potentials/albe.c b/potentials/albe.c
index 9d8336e..9f98282 100644 (file)
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -105,6 +105,15 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
        p->S2[0]=p->S[0]*p->S[0];
        p->S2[1]=p->S[1]*p->S[1];
        p->S2mixed=p->Smixed*p->Smixed;
        p->S2[0]=p->S[0]*p->S[0];
        p->S2[1]=p->S[1]*p->S[1];
        p->S2mixed=p->Smixed*p->Smixed;

+       p->c2[0]=p->c[0]*p->c[0];
+       p->c2[1]=p->c[1]*p->c[1];
+       p->c2_mixed=p->c_mixed*p->c_mixed;
+       p->d2[0]=p->d[0]*p->d[0];
+       p->d2[1]=p->d[1]*p->d[1];
+       p->d2_mixed=p->d_mixed*p->d_mixed;
+       p->c2d2[0]=p->c2[0]/p->d2[0];
+       p->c2d2[1]=p->c2[1]/p->d2[1];
+       p->c2d2_m=p->c2_mixed/p->d2_mixed;

 
        printf("[albe] mult parameter info:\n");
        printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
 
        printf("[albe] mult parameter info:\n");
        printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);


@@ -301,40 +310,18 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn,
                dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
        }
 
                dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
        }
 

-       /* 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);
-
-       /* 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);
+#ifdef DEBUG
+       if(ai==&(moldyn->atom[DATOM])) 
+               printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik);
+#endif

 
 

-       /* zeta derivatives */
-       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(dzjj,dzjj,&tmp);                         // j j
-       v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj);
-       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(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(dzkj,dzkj,&tmp);                         // k j
+       /* 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;

 
        /* increase k counter */
        exchange->kcnt+=1;
 
        /* increase k counter */
        exchange->kcnt+=1;


@@ -440,11 +427,21 @@ printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
        v3_scale(&force,&force,-1.0); // dri rij = - drj rij
        v3_add(&(aj->f),&(aj->f),&force);
 
        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);
-}
+       /* virial */
+       virial_calc(ai,&force,&(exchange->dist_ij));
+
+#ifdef DEBUG
+       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[DATOM]))
+                       printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
+               if(aj==&(moldyn->atom[DATOM]))
+                       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 */
 #endif
 
        /* virial */


@@ -501,39 +498,89 @@ int albe_mult_i0_j2_k0(t_moldyn *moldyn,
                return 0;
        }
 
                return 0;
        }
 

-       /* k!=j & check whether to run k */
-       k=exchange->kcnt;
-       j=exchange->j2cnt;
-       if((k==j)|(exchange->skip[k])) {
-               exchange->kcnt+=1;
-               return 0;
-       }
-       
-       /* force contribution (drk derivative) */
-       v3_scale(&force,&(exchange->dzeta[j][k]),exchange->pre_dzeta);
-       v3_add(&(ak->f),&(ak->f),&force);
+       /* prefactor dzeta */
+       pre_dzeta=exchange->pre_dzeta;

 
 

-#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);
-}
+       /* 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(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
 
        /* virial */
 #endif
 
        /* virial */

-       virial_calc(ai,&force,&(exchange->dist[k]));
+       virial_calc(ai,&force,&dist_ij);

 
 

+       /* force contribution to atom i */

        v3_scale(&force,&force,-1.0);
        v3_add(&(ai->f),&(ai->f),&force);
 
        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]);
-}
+       /* 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);
+
+       v3_scale(&force,&force,-1.0);
+       v3_add(&(ai->f),&(ai->f),&force);
+
+#ifdef DEBUG
+       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);
+       }

 #endif
 
 #endif
 

+       /* virial */
+       virial_calc(ai,&force,&dist_ik);
+       
+       /* force contribution to atom i */
+       v3_scale(&force,&force,-1.0);
+       v3_add(&(ai->f),&(ai->f),&force);
+

        /* increase k counter */
        exchange->kcnt+=1;
 
        /* increase k counter */
        exchange->kcnt+=1;