+
+/* tersoff 3 body potential function (second k loop) */
+int tersoff_mult_3bp_k2(t_moldyn *moldyn,
+ t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+
+ t_tersoff_mult_params *params;
+ t_tersoff_exchange *exchange;
+ int kcount;
+ t_3dvec dist_ik,dist_ij;
+ double d_ik2,d_ik,d_ij2,d_ij;
+ unsigned char brand;
+ double S2;
+ double g,dg,cos_theta;
+ double pre_dzeta;
+ double f_c_ik,df_c_ik;
+ double dijdik_inv,fcdg,dfcg;
+ t_3dvec dcosdri,dcosdrj,dcosdrk;
+ t_3dvec force,tmp;
+
+ params=moldyn->pot_params;
+ exchange=&(params->exchange);
+ kcount=exchange->kcount;
+
+ if(kcount>TERSOFF_MAXN)
+ printf("FATAL: neighbours!\n");
+
+ /* d_ik2 */
+ d_ik2=exchange->d_ik2[kcount];
+
+ brand=ak->brand;
+ if(brand==ai->brand)
+ S2=params->S2[brand];
+ else
+ S2=params->S2mixed;
+
+ /* return if d_ik > S */
+ if(d_ik2>S2) {
+ exchange->kcount++;
+ return 0;
+ }
+
+ /* prefactor dzeta */
+ pre_dzeta=exchange->pre_dzeta;
+
+ /* dist_ik, d_ik */
+ dist_ik=exchange->dist_ik[kcount];
+ d_ik=exchange->d_ik[kcount];
+
+ /* f_c_ik, df_c_ik */
+ f_c_ik=exchange->f_c_ik[kcount];
+ df_c_ik=exchange->df_c_ik[kcount];
+
+ /* dist_ij, d_ij2, d_ij */
+ dist_ij=exchange->dist_ij;
+ d_ij2=exchange->d_ij2;
+ d_ij=exchange->d_ij;
+
+ /* g, dg, cos_theta */
+ g=exchange->g[kcount];
+ dg=exchange->dg[kcount];
+ cos_theta=exchange->cos_theta[kcount];
+
+ /* cos_theta derivatives wrt i,j,k */
+ dijdik_inv=1.0/(d_ij*d_ik);
+ v3_scale(&dcosdrj,&dist_ik,dijdik_inv);
+ v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
+ v3_add(&dcosdrj,&dcosdrj,&tmp);
+ v3_scale(&dcosdrk,&dist_ij,dijdik_inv);
+ v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
+ v3_add(&dcosdrk,&dcosdrk,&tmp);
+ v3_add(&dcosdri,&dcosdrj,&dcosdrk);
+ v3_scale(&dcosdri,&dcosdri,-1.0);
+
+ /* f_c_ik * dg, df_c_ik * g */
+ fcdg=f_c_ik*dg;
+ dfcg=df_c_ik*g;
+
+ /* derivative wrt i */
+ v3_scale(&force,&dist_ik,dfcg);
+ v3_scale(&tmp,&dcosdri,fcdg);
+ v3_add(&force,&force,&tmp);
+ v3_scale(&force,&force,pre_dzeta);
+
+ /* force contribution */
+ v3_add(&(ai->f),&(ai->f),&force);
+
+#ifdef DEBUG
+ if(ai==&(moldyn->atom[0])) {
+ printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
+ printf("adding %f %f %f\n",force.x,force.y,force.z);
+ printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
+ }
+#endif
+
+ /* derivative wrt j */
+ v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
+
+ /* force contribution */
+ v3_add(&(aj->f),&(aj->f),&force);
+
+#ifdef DEBUG
+ if(aj==&(moldyn->atom[0])) {
+ printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
+ printf("adding %f %f %f\n",force.x,force.y,force.z);
+ printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
+ }