+ /* reset 3bp sums */
+ exchange->3bp_sum1=0.0;
+ exchange->3bp_sum2=0.0;
+
+ return 0;
+}
+
+/* tersoff 2 body post part */
+
+int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+
+ /* here we have to allow for the 3bp sums */
+
+ t_tersoff_mult_params *params;
+ t_tersoff_exchange *exchange;
+
+ t_3dvec force,temp,*db_ij;
+ double db_ij_scale1,db_ij_scale2;
+ double b_ij;
+ double f_c,df_c,f_a,df_a;
+
+ params=moldyn->pot2b_params;
+ exchange=&(moldyn->exchange);
+
+ db_ij=&(exchange->db_ij);
+ f_c=exchange->f_c;
+ df_c=exchange->df_c;
+ f_a=exchange->f_a;
+ df_a=exchange->df_a;
+
+ db_ij_scale1=(1+betan*3bp_sum1);
+ db_ij_scale2=(n*betan*3bp_sum2);
+ help=pow(db_ij_scale1,-1.0/(2*n)-1);
+ b_ij=chi*db_ij_scale1*help;
+ db_ij_scale1=-chi/(2*n)*help;
+
+ v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2));
+ v3_scale(db_ij,db_ij,f_a);
+
+ v3_scale(&temp,dist_ij,b_ij*df_a);
+
+ v3_add(&force,&temp,db_ij);
+ v3_scale(&force,&force,f_c);
+
+ v3_scale(&temp,&dist_ij,f_a*b_ij*df_c);
+
+ /* add energy of 3bp sum */
+ moldyn->energy+=(0.5*f_c*b_ij*f_a);
+ /* add force of 3bp calculation */
+ v3_add(&(ai->f),&temp,&force);
+