+ exchange->run2bp_post=1;
+
+ /* reset 3bp sums */
+ exchange->zeta=0.0;
+ v3_zero(&(exchange->db_ij));
+
+ return 0;
+}
+
+/* tersoff 2 body post part */
+
+int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,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,*dist_ij;
+ double db_ij_scale1,db_ij_scale2;
+ double b_ij;
+ double f_c,df_c,f_a,df_a;
+ double chi,n,n_betan;
+ 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;
+
+ 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;
+ n_betan=exchange->n_betan;
+ n=*(exchange->n);
+ chi=exchange->chi;
+ dist_ij=&(exchange->dist_ij);
+ zeta=exchange->zeta;
+
+ db_ij_scale2=pow(zeta,n-1.0);
+printf("DEBUG: %.15f %.15f\n",zeta,db_ij_scale2);
+ db_ij_scale1=db_ij_scale2*zeta;
+ db_ij_scale2*=n_betan;
+ db_ij_scale1=pow((1.0+n_betan*db_ij_scale1),-1.0/(2*n)-1);
+ b_ij=chi*db_ij_scale1*(1.0+n_betan*db_ij_scale1);
+ db_ij_scale1*=(-1.0*chi/(2*n));
+
+ /* db_ij part */
+ v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2));
+ v3_scale(db_ij,db_ij,f_a);
+
+ /* df_a part */
+ v3_scale(&temp,dist_ij,b_ij*df_a);
+
+ /* db_ij + df_a part */
+ v3_add(&force,&temp,db_ij);
+ v3_scale(&force,&force,f_c);
+
+ /* df_c part */
+ 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 (all three parts) */
+ v3_add(&(ai->f),&temp,&force);