- /*
- * here we have to allow for the 3bp sums
- *
- * that is:
- * - zeta_ij, dzeta_ij
- * - zeta_ji, dzeta_ji
- *
- * to compute the 3bp contribution to:
- * - Vij, dVij
- * - dVji
- *
- */
-
- t_tersoff_mult_params *params;
- t_tersoff_exchange *exchange;
-
- t_3dvec force,temp;
- t_3dvec *dist_ij;
- double b,db,tmp;
- double f_c,df_c,f_a,df_a;
- double chi,ni,betaini,nj,betajnj;
- double zeta;
-
- params=moldyn->pot2b_params;
- exchange=&(params->exchange);
-
- /* we do not run if f_c_ij was detected to be 0! */
- if(!(exchange->run2bp_post))
- return 0;
-
- f_c=exchange->f_c;
- df_c=exchange->df_c;
- f_a=exchange->f_a;
- df_a=exchange->df_a;
- betaini=exchange->betaini;
- betajnj=exchange->betajnj;
- ni=*(exchange->n_i);
- nj=*(exchange->n_j);
- chi=exchange->chi;
- dist_ij=&(exchange->dist_ij);
-
- /* Vij and dVij */
- zeta=exchange->zeta_ij;
- 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); /* chi * (...)^(-1/2n - 1) */
- b=db*b; /* b_ij */
- db*=-0.5*tmp; /* db_ij */
- v3_scale(&force,&(exchange->dzeta_ij),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 energy of 3bp sum */
- 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);
-
- /* add force (sub, as F = - dVji) */
- v3_sub(&(ai->f),&(ai->f),&force);