-
- v3_sub(&dist_jk,&(aj->r),&(ak->r));
- if(bc) check_per_bound(moldyn,&dist_jk);
- d_jk=v3_norm(&dist_jk);
-
- beta=*(exchange->beta);
- betan=exchange->betan;
- n=*(exchange->n);
- c=*(exchange->c);
- d=*(exchange->d);
- h=*(exchange->h);
- c2=exchange->c2;
- d2=exchange->d2;
- c2d2=exchange->c2d2;
-
- numer=d_ij2+d_ik*d_ik-d_jk*d_jk;
- denom=2*d_ij*d_ik;
- cos_theta=numer/denom;
- sin_theta=sqrt(1.0-(cos_theta*cos_theta));
- theta=acos(cos_theta);
- d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom);
- d_theta1=2*denom-numer*2*d_ik/d_ij;
- d_theta2=2*denom-numer*2*d_ij/d_ik;
- d_theta1*=d_theta;
- d_theta2*=d_theta;
-
- h_cos=(h-cos_theta);
- h_cos2=h_cos*h_cos;
- d2_h_cos2=d2-h_cos2;
-
- /* some usefull expressions */
- frac1=c2/(d2-h_cos2);
- bracket1=1+c2d2-frac1;
- bracket2=f_c_ik*bracket1;
- bracket2_n_1=pow(bracket2,n-1.0);
- bracket2_n=bracket2_n_1*bracket2;
- bracket3=1+betan*bracket2_n;
- bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0);
- bracket3_pow=bracket3_pow_1*bracket3;
-
- /* now go on with calc of b_ij and derivation of b_ij */
- b_ij=chi*bracket3_pow;
-
- /* derivation of theta */
- v3_scale(&force,&dist_ij,d_theta1);
- v3_scale(&temp,&dist_ik,d_theta2);
- v3_add(&force,&force,&temp);
-
- /* part 1 of derivation of b_ij */
- v3_scale(&force,&force,sin_theta*2*h_cos*f_c_ik*frac1);
-
- /* part 2 of derivation of b_ij */
- v3_scale(&temp,&dist_ik,df_c_ik*bracket1);
-
- /* sum up and scale ... */
- v3_add(&temp,&temp,&force);
- scale=bracket2_n_1*n*betan*(1+betan*bracket3_pow_1)*chi*(1.0/(2.0*n));
- v3_scale(&temp,&temp,scale);
-
- /* now construct an energy and a force out of that */
- v3_scale(&temp,&temp,f_a);
- v3_scale(&force,&dist_ij,df_a*b_ij);
- v3_add(&temp,&temp,&force);
- v3_scale(&temp,&temp,f_c);
- v3_scale(&force,&dist_ij,df_c*b_ij*f_a);
- v3_add(&force,&force,&temp);
-
- /* add forces */
- v3_add(&(ai->f),&(ai->f),&force);
- /* energy is 0.5 f_r f_c, but we will sum it up twice ... */
- moldyn->energy+=(0.25*f_a*b_ij*f_c);
-