else {
tmp=betaini*pow(zeta,ni-1.0); /* beta^n * zeta^n-1 */
b=(1+zeta*tmp); /* 1 + beta^n zeta^n */
- db=chi*pow(b,-1.0/(2*ni)-1); /* x(...)^(-1/2n - 1) */
+ db=chi*pow(b,-1.0/(2.0*ni)-1); /* x(...)^(-1/2n - 1) */
b=db*b; /* b_ij */
db*=-0.5*tmp; /* db_ij */
v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
v3_sub(&temp1,&temp1,&temp2); /* there is a minus! */
v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
- /* store dg in temp2 and use it for dVjk later */
+ /* store dg for use in dV_jk */
v3_copy(&temp2,&temp1);
/* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */
// dzeta_jk is simply dg, which is stored in temp2
}
/* betajnj * zeta_jk ^ nj-1 */
+printf("FATAL db_jk calc!\n");
+printf("(z1 + z2 + z3 ...)^n != z1^n + z2^n + z3^n + ...\n");
+printf("st00pid me => tersoff_orig is obsolete!\n");
tmp=exchange->betajnj*pow(zeta,(n-1.0));
tmp=-chi/2.0*pow((1+tmp*zeta),(-1.0/(2.0*n)-1))*tmp;
+#ifdef DEBUG
+ if((ai->tag==0)&(aj->tag==864)) { // &(ak->tag==23)) {
+ printf("\n\n");
+ printf("db: ni zeta = %f %f\n",
+ n,zeta);
+ printf("\n\n");
+ }
+#endif
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 ^ */