+ params=moldyn->pot_params;
+ num=ai->bnum;
+ exchange=&(params->exchange);
+
+ exchange->run3bp=0;
+
+ /*
+ * we need: f_c, df_c, f_r, df_r
+ *
+ * therefore we need: R, S, A, lambda
+ */
+
+ v3_sub(&dist_ij,&(ai->r),&(aj->r));
+
+ if(bc) check_per_bound(moldyn,&dist_ij);
+
+ /* save for use in 3bp */ /* REALLY ?!?!?! */
+ exchange->dist_ij=dist_ij;
+
+ /* constants */
+ if(num==aj->bnum) {
+ S=params->S[num];
+ R=params->R[num];
+ A=params->A[num];
+ lambda=params->lambda[num];
+ /* more constants depending of atoms i and j, needed in 3bp */
+ params->exchange.B=&(params->B[num]);
+ params->exchange.mu=params->mu[num];
+ params->exchange.chi=1.0;
+ }
+ else {
+ S=params->Smixed;
+ R=params->Rmixed;
+ A=params->Amixed;
+ lambda=params->lambda_m;
+ /* more constants depending of atoms i and j, needed in 3bp */
+ params->exchange.B=&(params->Bmixed);
+ params->exchange.mu=&(params->mu_m);
+ params->exchange.chi=params->chi;
+ }
+
+ d_ij=v3_norm(&dist_ij);
+
+ /* save for use in 3bp */
+ exchange->d_ij=d_ij;
+
+ if(d_ij>S)
+ return 0;
+
+ f_r=A*exp(-lamda*d_ij);
+ df_r=-lambda*f_r/d_ij;
+
+ /* f_a, df_a calc + save for 3bp use */
+ exchange->f_a=-B*exp(-mu*d_ij);
+ exchange->df_a=-mu*exchange->f_a/d_ij;
+
+ if(d_ij<R) {
+ /* f_c = 1, df_c = 0 */
+ f_c=1.0;
+ df_c=0.0;
+ v3_scale(&force,&dist_ij,df_r);
+ }
+ else {
+ s_r=S-R;
+ arg=PI*(d_ij-R)/s_r;
+ f_c=0.5+0.5*cos(arg);
+ df_c=-0.5*sin(arg)*(PI/(s_r*d_ij));
+ scale=df_c*f_r+df_r*f_c;
+ v3_scale(&force,&dist_ij,scale);
+ }
+
+ /* add forces */
+ v3_add(&(ai->f),&(ai->f),&force);
+ /* energy is 0.5 f_r f_c, but we will sum it up twice ... */
+ moldyn->energy+=(0.25*f_r*f_c);
+
+ /* save for use in 3bp */
+ exchange->f_c=f_c;
+ exchange->df_c=df_c;
+
+ /* enable the run of 3bp function */
+ exchange->run3bp=1;
+
+ 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 db_ij,temp,force;
+ double R,S,s_r;
+ double d_ij,d_ik,d_jk;
+ double f_c,df_c,b_ij,f_a,df_a;
+ double B,mu;
+ int num;