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);
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;
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 */
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 */
- 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);
-#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
+ /* 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;