- /* 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);
+
+#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
+
+ /* virial */
+ virial_calc(ai,&force,dist);
+
+ 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