int brand;
double s_r;
double arg;
+ double energy;
printf("WARNING! - tersoff_mult_2bp is obsolete.\n");
printf("WARNING! - repulsive part handled in 3bp/j2 routine.\n");
virial_calc(ai,&force,&dist_ij);
/* energy 2bp contribution */
- moldyn->energy+=f_r*f_c;
+ energy=f_r*f_c;
+ moldyn->energy+=energy;
+ ai->e+=0.5*energy;
+ aj->e+=0.5*energy;
return 0;
}
exchange->zeta_ij+=f_c_ik*g;
}
+#ifdef DEBUG
+ if((ai==&(moldyn->atom[0]))|
+ (aj==&(moldyn->atom[864]))|
+ (ak==&(moldyn->atom[1003]))) {
+ printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos);
+ }
+#endif
+
/* store even more data for second k loop */
exchange->g[kcount]=g;
exchange->dg[kcount]=dg;
unsigned char brand;
double ni,tmp;
double S,R,s_r,arg;
+ double energy;
params=moldyn->pot_params;
exchange=&(params->exchange);
printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
if(aj==&(moldyn->atom[0]))
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
exchange->pre_dzeta=-0.5*f_a*f_c*db;
/* energy contribution */
- moldyn->energy+=0.5*f_c*(b*f_a+f_r);
+ energy=0.5*f_c*(b*f_a+f_r);
+ moldyn->energy+=energy;
+ ai->e+=energy;
/* reset k counter for second k loop */
exchange->kcount=0;