- /* store even more data for second k loop */
- exchange->g[kcount]=g;
- exchange->dg[kcount]=dg;
- exchange->d_ik[kcount]=d_ik;
- exchange->cos_theta[kcount]=cos_theta;
- exchange->f_c_ik[kcount]=f_c_ik;
- exchange->df_c_ik[kcount]=df_c_ik;
+ /* zeta */
+ exchange->zeta[j]+=(exchange->f_c[k]*gk);
+ exchange->zeta[k]+=(exchange->f_c[j]*gj);
+
+ /* cos theta derivatives */
+ v3_scale(&dcosdrj,&distk,djdk_inv); // j
+ v3_scale(&tmp,&distj,-cos_theta/exchange->d2[j]);
+ v3_add(&dcosdrj,&dcosdrj,&tmp);
+ v3_scale(&dcosdrk,&distj,djdk_inv); // k
+ v3_scale(&tmp,&distk,-cos_theta/exchange->d2[k]);
+ v3_add(&dcosdrk,&dcosdrk,&tmp);
+
+ /* zeta derivatives */
+ dzeta=exchange->dzeta;
+ v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk);
+ v3_add(&dzeta[j][j],&dzeta[j][j],&tmp); // j j
+ v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj);
+ v3_add(&dzeta[k][k],&dzeta[k][k],&tmp); // k k
+ v3_scale(&tmp,&distk,exchange->df_c[k]*gk/dk);
+ v3_add(&dzeta[j][k],&dzeta[j][k],&tmp);
+ v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk);
+ v3_add(&dzeta[j][k],&dzeta[j][k],&tmp); // j k
+ v3_scale(&tmp,&distj,exchange->df_c[j]*gj/dj);
+ v3_add(&dzeta[k][j],&dzeta[k][j],&tmp);
+ v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj);
+ v3_add(&dzeta[k][j],&dzeta[k][j],&tmp); // k j