-/* tersoff 1 body part */
-int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
-
- int num;
- t_tersoff_mult_params *params;
- t_tersoff_exchange *exchange;
-
- num=ai->bnum;
- params=moldyn->pot1b_params;
- exchange=&(params->exchange);
-
- /*
- * simple: point constant parameters only depending on atom i to
- * their right values
- */
-
- exchange->beta_i=&(params->beta[num]);
- exchange->n_i=&(params->n[num]);
- exchange->c_i=&(params->c[num]);
- exchange->d_i=&(params->d[num]);
- exchange->h_i=&(params->h[num]);
-
- exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i));
- exchange->ci2=params->c[num]*params->c[num];
- exchange->di2=params->d[num]*params->d[num];
- exchange->ci2di2=exchange->ci2/exchange->di2;
-
- 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;
-
- params=moldyn->pot2b_params;
- num=aj->bnum;
- exchange=&(params->exchange);
-
- /* clear 3bp and 2bp post run */
- exchange->run3bp=0;
- exchange->run2bp_post=0;
-
- /* reset S > r > R mark */
- exchange->d_ij_between_rs=0;
-
- /*
- * calc of 2bp contribution of V_ij and dV_ij/ji
- *
- * for Vij and dV_ij we need:
- * - f_c_ij, df_c_ij
- * - f_r_ij, df_r_ij
- *
- * for dV_ji we need:
- * - f_c_ji = f_c_ij, df_c_ji = df_c_ij
- * - f_r_ji = f_r_ij; df_r_ji = df_r_ij
- *
- */
-
- /* dist_ij, d_ij */
- v3_sub(&dist_ij,&(aj->r),&(ai->r));
- if(bc) check_per_bound(moldyn,&dist_ij);
- d_ij=v3_norm(&dist_ij);
-
- /* save for use in 3bp */
- exchange->d_ij=d_ij;
- exchange->dist_ij=dist_ij;
-
- /* constants */
- if(num==ai->bnum) {
- S=params->S[num];
- R=params->R[num];
- A=params->A[num];
- B=params->B[num];
- lambda=params->lambda[num];
- mu=params->mu[num];
- exchange->chi=1.0;
- }
- else {
- S=params->Smixed;
- R=params->Rmixed;
- A=params->Amixed;
- B=params->Bmixed;
- lambda=params->lambda_m;
- mu=params->mu_m;
- params->exchange.chi=params->chi;
- }
-
- /* if d_ij > S => no force & potential energy contribution */
- if(d_ij>S)
- return 0;
-
- /* more constants */
- exchange->beta_j=&(params->beta[num]);
- exchange->n_j=&(params->n[num]);
- exchange->c_j=&(params->c[num]);
- exchange->d_j=&(params->d[num]);
- exchange->h_j=&(params->h[num]);
- if(num==ai->bnum) {
- exchange->betajnj=exchange->betaini;
- exchange->cj2=exchange->ci2;
- exchange->dj2=exchange->di2;
- exchange->cj2dj2=exchange->ci2di2;
- }
- else {
- exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j));
- exchange->cj2=params->c[num]*params->c[num];
- exchange->dj2=params->d[num]*params->d[num];
- exchange->cj2dj2=exchange->cj2/exchange->dj2;
- }