+ moldyn->energy+=(0.5*f_c*b*f_a);
+
+ /* add force (sub, as F = - dVij) */
+ v3_sub(&(ai->f),&(ai->f),&force);
+
+ /* dVji */
+ zeta=exchange->zeta_ji;
+ tmp=betajnj*pow(zeta,nj-1.0); /* beta^n * zeta^n-1 */
+ b=(1+zeta*tmp); /* 1 + beta^n * zeta^n */
+ db=chi*pow(b,-1.0/(2*nj)-1); /* chi * (...)^(-1/2n - 1) */
+ b=db*b; /* b_ij */
+ db*=-0.5*tmp; /* db_ij */
+ v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
+ v3_scale(&temp,dist_ij,df_a*b);
+ v3_add(&force,&force,&temp);
+ v3_scale(&force,&force,f_c);
+ v3_scale(&temp,dist_ij,df_c*b*f_a);
+ v3_add(&force,&force,&temp);