-/* 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,betan;
- double help;
- double n;
-
- 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;
- betan=exchange->betan;
- n=*(exchange->n);
- chi=exchange->chi;
- dist_ij=&(exchange->dist_ij);
-
- db_ij_scale1=(1+betan*exchange->sum1_3bp);
- db_ij_scale2=(exchange->n_betan*exchange->sum2_3bp);
- help=pow(db_ij_scale1,-1.0/(2*n)-1);
- b_ij=chi*db_ij_scale1*help;
- db_ij_scale1=-chi/(2*n)*help;
-printf("debug: %.20f %.20f %.20f\n",db_ij->x,exchange->sum1_3bp,exchange->sum2_3bp);
-
- /* 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);
-
- return 0;
-}
-
-/* tersoff 3 body part */
-
-int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
-
- t_tersoff_mult_params *params;
- t_tersoff_exchange *exchange;
- t_3dvec dist_ij,dist_ik,dist_jk;
- t_3dvec temp,force;
- double R,S,s_r;
- double d_ij,d_ij2,d_ik,d_jk;
- double f_c,df_c,f_a,df_a;
- double f_c_ik,df_c_ik,arg;
- double n,c,d,h;
- double c2,d2,c2d2;
- double numer,denom;
- double theta,cos_theta,sin_theta;
- double d_theta,d_theta1,d_theta2;
- double h_cos,d2_h_cos2;
- double frac,bracket,bracket_n_1,bracket_n;
- double g;
- int num;
-
- params=moldyn->pot3b_params;
- num=ai->bnum;
- exchange=&(params->exchange);
-
- if(!(exchange->run3bp))
- return 0;
-
- /*
- * we need: f_c, d_fc, b_ij, db_ij, f_a, df_a
- *
- * we got f_c, df_c, f_a, df_a from 2bp calculation
- */
-
- d_ij=exchange->d_ij;
- d_ij2=exchange->d_ij2;
- dist_ij=exchange->dist_ij;
-
- f_a=params->exchange.f_a;
- df_a=params->exchange.df_a;