- 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);
-
- /* add force (sub, as F = - dVji) */
- v3_sub(&(ai->f),&(ai->f),&force);
+ if(zeta==0.0) {
+ // do nothing ... (db=0, b=chi)
+ moldyn->debug++;
+ }
+ else {
+ 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); /* x(...)^(-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,0.5*df_c*b*f_a);
+ v3_add(&force,&force,&temp);
+ v3_scale(&force,&force,-0.5);
+
+ /* add force */
+ v3_sub(&(ai->f),&(ai->f),&force);
+ }