+ v3_sub(&dist_ik,ktom,&(atom[i]));
+ d_ik=v3_norm(&dist_ik);
+ if(d_ik<=Sik) {
+
+ Rik=;
+ Sik=;
+ Aik=;
+ lambda_ik=;
+ Bik=;
+ mu_ik=;
+ omega_ik=;
+ c_i=;
+ d_i=;
+ h_i=;
+
+
+ if(d_ik<=Rik) {
+ f_cik=1.0;
+ df_cik=0.0;
+ }
+ else {
+ sik_rik=Sik-Rik;
+ arg1ik=PI*(d_ik-Rik)/sik_rik;
+ f_cik=0.5+0.5*cos(arg1ik);
+ df_cik=-0.5*sin(arg1ik)* \
+ (PI/(sik_rik*d_ik));
+ f_rik=Aik*exp(-lambda_ik*d_ik);
+ f_aik=-Bik*exp(-mu_ik*d_ik);
+ }
+
+ v3_sub(&distance_jk,ktom,btom);
+ cos_theta=(d_ij2+d_ik*d_ik-d_jk*d_jk)/\
+ (2*d_ij*d_ik);
+ theta=arccos(cos_theta);