+/*
+ * tersoff potential & force for 2 sorts of atoms
+ */
+
+/* create mixed terms from parameters and set them */
+int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
+
+ printf("[moldyn] tersoff parameter completion\n");
+ p->Smixed=sqrt(p->S[0]*p->S[1]);
+ p->Rmixed=sqrt(p->R[0]*p->R[1]);
+ p->Amixed=sqrt(p->A[0]*p->A[1]);
+ p->Bmixed=sqrt(p->B[0]*p->B[1]);
+ p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]);
+ p->mu_m=0.5*(p->mu[0]+p->mu[1]);
+
+ printf("[moldyn] tersoff mult parameter info:\n");
+ printf(" S (m) | %.12f | %.12f | %.12f\n",p->S[0],p->S[1],p->Smixed);
+ printf(" R (m) | %.12f | %.12f | %.12f\n",p->R[0],p->R[1],p->Rmixed);
+ printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
+ printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
+ printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
+ p->lambda_m);
+ printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
+ printf(" beta | %.10f | %.10f\n",p->beta[0],p->beta[1]);
+ printf(" n | %f | %f\n",p->n[0],p->n[1]);
+ printf(" c | %f | %f\n",p->c[0],p->c[1]);
+ printf(" d | %f | %f\n",p->d[0],p->d[1]);
+ printf(" h | %f | %f\n",p->h[0],p->h[1]);
+ printf(" chi | %f \n",p->chi);
+
+ return 0;
+}
+
+/* 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=&(params->beta[num]);
+ exchange->n=&(params->n[num]);
+ exchange->c=&(params->c[num]);
+ exchange->d=&(params->d[num]);
+ exchange->h=&(params->h[num]);
+
+ exchange->betan=pow(*(exchange->beta),*(exchange->n));
+ exchange->n_betan=*(exchange->n)*exchange->betan;
+ exchange->c2=params->c[num]*params->c[num];
+ exchange->d2=params->d[num]*params->d[num];
+ exchange->c2d2=exchange->c2/exchange->d2;
+
+ 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;
+ exchange->run2bp_post=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);
+
+ d_ij=v3_norm(&dist_ij);
+
+ /* save for use in 3bp */
+ exchange->d_ij=d_ij;
+ exchange->dist_ij=dist_ij;
+ exchange->d_ij2=d_ij*d_ij;
+
+ /* constants */
+ if(num==aj->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];
+ params->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)
+ return 0;
+
+ f_r=A*exp(-lambda*d_ij);
+ df_r=-lambda*f_r/d_ij;
+
+ /* f_a, df_a calc + save for later 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);