+ /* dV_jk stuff | add force contribution on atom i immediately */
+ if(exchange->d_ij_between_rs) {
+ zeta=f_c*g;
+ v3_scale(&temp1,&temp2,f_c);
+ v3_scale(&temp2,&dist_ij,df_c*g);
+ v3_add(&temp2,&temp2,&temp1); /* -> dzeta_jk in temp2 */
+ }
+ else {
+ zeta=g;
+ // dzeta_jk is simply dg, which is stored in temp2
+ }
+ /* betajnj * zeta_jk ^ nj-1 */
+ tmp=exchange->betajnj*pow(zeta,(n-1.0));
+ tmp=-chi/2.0*pow((1+tmp*zeta),(-1.0/(2.0*n)-1))*tmp;
+ v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
+ v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
+ /* scaled with 0.5 ^ */
+
+ /* virial */
+ ai->virial.xx-=temp2.x*dist_jk.x;
+ ai->virial.yy-=temp2.y*dist_jk.y;
+ ai->virial.zz-=temp2.z*dist_jk.z;
+ ai->virial.xy-=temp2.x*dist_jk.y;
+ ai->virial.xz-=temp2.x*dist_jk.z;
+ ai->virial.yz-=temp2.y*dist_jk.z;
+
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVjk (3bp) contrib:\n");
+ printf("%f | %f\n",temp2.x,ai->f.x);
+ printf("%f | %f\n",temp2.y,ai->f.y);
+ printf("%f | %f\n",temp2.z,ai->f.z);
+}
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVjk (3bp) contrib:\n");
+ printf("%f | %f\n",temp2.x*dist_jk.x,ai->virial.xx);
+ printf("%f | %f\n",temp2.y*dist_jk.y,ai->virial.yy);
+ printf("%f | %f\n",temp2.z*dist_jk.z,ai->virial.zz);
+}
+#endif