+ /* 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(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
+ v3_add(&(ai->f),&(ai->f),&temp1); /* -1 skipped in f_a calc ^ */
+ /* scaled with 0.5 ^ */
+ }
+
+ return 0;
+}
+
+
+/*
+ * debugging / critical check functions
+ */
+
+int moldyn_bc_check(t_moldyn *moldyn) {
+
+ t_atom *atom;
+ t_3dvec *dim;
+ int i;
+
+ atom=moldyn->atom;
+ dim=&(moldyn->dim);
+
+ for(i=0;i<moldyn->count;i++) {
+ if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2)
+ printf("FATAL: atom %d: x: %.20f (%.20f)\n",
+ i,atom[i].r.x*1e10,dim->x/2*1e10);
+ if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
+ printf("FATAL: atom %d: y: %.20f (%.20f)\n",
+ i,atom[i].r.y*1e10,dim->y/2*1e10);
+ if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
+ printf("FATAL: atom %d: z: %.20f (%.20f)\n",
+ i,atom[i].r.z*1e10,dim->z/2*1e10);