- return 0;
-}
-
-/* tersoff 2 body part */
-int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
-
- t_tersoff_mult_params *params;
- t_tersoff_exchange *exchange;
- t_3dvec dist_ij,force;
- double d_ij;
- double A,B,R,S,lambda,mu;
- double f_r,df_r;
- double f_c,df_c;
- int num;
- double s_r;
- double arg;
- double scale;
-
- params=moldyn->pot2b_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]);
- 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);
- 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(-lambda*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=M_PI*(d_ij-R)/s_r;
- f_c=0.5+0.5*cos(arg);
- df_c=-0.5*sin(arg)*(M_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 temp,force;
- double R,S,s_r;
- double d_ij,d_ij2,d_ik,d_jk;
- double f_c,df_c,b_ij,f_a,df_a;
- double f_c_ik,df_c_ik,arg;
- double scale;
- double chi;
- double n,c,d,h,beta,betan;
- double c2,d2,c2d2;
- double numer,denom;
- double theta,cos_theta,sin_theta;
- double d_theta,d_theta1,d_theta2;
- double h_cos,h_cos2,d2_h_cos2;
- double frac1,bracket1,bracket2,bracket2_n_1,bracket2_n;
- double bracket3,bracket3_pow_1,bracket3_pow;
- 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;
-
- f_a=params->exchange.f_a;
- df_a=params->exchange.df_a;
-
- /* d_ij is <= S, as we didn't return so far! */
-
- /*
- * calc of b_ij (scalar) and db_ij (vector)
- *
- * - for b_ij: chi, beta, f_c_ik, w(=1), c, d, h, n, cos_theta
- *
- * - for db_ij: d_theta, sin_theta, cos_theta, f_c_ik, df_c_ik,
- * w_ik,
- *
- */
-
-
- v3_sub(&dist_ik,&(ai->r),&(ak->r));
- if(bc) check_per_bound(moldyn,&dist_ik);
- d_ik=v3_norm(&dist_ik);
-
- /* constants for f_c_ik calc */
- if(num==ak->bnum) {
- R=params->R[num];
- S=params->S[num];
- }
- else {
- R=params->Rmixed;
- S=params->Smixed;
- }