+ /* 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 */
+ dzjj=&(exchange->dzeta[j][j]);
+ dzkk=&(exchange->dzeta[k][k]);
+ dzjk=&(exchange->dzeta[j][k]);
+ dzkj=&(exchange->dzeta[k][j]);
+ v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk);
+ v3_add(dzjj,dzjj,&tmp); // j j
+ v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj);
+ v3_add(dzkk,dzkk,&tmp); // k k
+ v3_scale(&tmp,&distk,-exchange->df_c[k]*gk); // dk rik = - di rik
+ v3_add(dzjk,dzjk,&tmp);
+ v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk);
+ v3_add(dzjk,dzjk,&tmp); // j k
+ v3_scale(&tmp,&distj,-exchange->df_c[j]*gj); // dj rij = - di rij
+ v3_add(dzkj,dzkj,&tmp);
+ v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj);
+ v3_add(dzkj,dzkj,&tmp); // k j