+ /* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */
+ if(d_jk<S) {
+
+ /* constants_j from exchange data */
+ n=*(exchange->n_j);
+ c=*(exchange->c_j);
+ d=*(exchange->d_j);
+ h=*(exchange->h_j);
+ c2=exchange->cj2;
+ d2=exchange->dj2;
+ c2d2=exchange->cj2dj2;
+
+ /* cosine of theta_jik by scalaproduct */
+ rxxryy=v3_scalar_product(&dist_ij,&dist_jk); /* times -1 */
+ dxxdyy=d_ij*d_jk;
+ cos_theta=rxxryy/dxxdyy;
+
+ /* d_costheta */
+ d_costheta1=1.0/(d_jk*d_ij);
+ d_costheta2=cos_theta/(d_ij*d_ij); /* in fact -cos(), but ^ */
+
+ /* some usefull values */
+ h_cos=(h-cos_theta);
+ d2_h_cos2=d2+(h_cos*h_cos);
+ frac=c2/(d2_h_cos2);
+
+ /* g(cos_theta) */
+ g=1.0+c2d2-frac;
+
+ /* d_costheta_ij and dg(cos_theta) - needed in any case! */
+ v3_scale(&temp1,&dist_jk,d_costheta1);
+ v3_scale(&temp2,&dist_ij,-d_costheta2); /* ji -> ij => -1 */
+ v3_add(&temp1,&temp1,&temp2);
+ v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
+
+ /* dV_jk stuff | add force contribution on atom i immediately */
+ if(exchange->d_ij_between_rs) {
+ tmp=pow(f_c_ij*g,(n_j-1.0)); /* zeta_jk ^ n_j-1 */
+ v3_scale(&temp2,&temp1,f_c_ij)
+ v3_scale(&temp3,&dist_ij,df_c_ij);
+ v3_add(&temp3,&temp3,&temp2); /* dzeta_jk */
+ }
+ else {
+ /* f_c_ij = 1, df_c_ij = 0 */
+ tmp=pow(g,(n_j-1.0)); /* zeta_jk ^ n_j-1 */
+ tmp
+ /* dzeta_jk in temp1 */
+ /* HIER WEITER !!! */
+ }