-
- /* 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 j,k */
- dijdik_inv=1.0/(d_ij*d_ik);
- v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
- v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
- v3_add(&dcosdrj,&dcosdrj,&tmp);
- v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
- v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
- v3_add(&dcosdrk,&dcosdrk,&tmp);
-
- /* f_c_ik * dg, df_c_ik * g */
- fcdg=f_c_ik*dg;
- dfcg=df_c_ik*g;
-
- /* derivative wrt j */
- v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
-
- /* force contribution */
- v3_add(&(aj->f),&(aj->f),&force);
-
-#ifdef DEBUG
-if(moldyn->time>DSTART&&moldyn->time<DEND) {
- if(aj==&(moldyn->atom[DATOM])) {
- 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);
- printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
- printf(" d ij ik = %f %f\n",d_ij,d_ik);
- }
-}
-#endif
-
- /* force contribution to atom i */