lambda=;
B=;
mu=;
+ chi=;
+ beta=;
+ betaN=;
if(d_ij<=R) {
f_c=1.0;
ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
neighbourk);
+ /* go for zeta - 3 body stuff! */
+ zeta=0.0;
+ d_ij2=d_ij*d_ij;
+
/* cell of btom */
thisk=&(neighbourk[0]);
list_reset(thisk);
/* 3 body stuff (1) */
- theta_ijk=;
- sin_theta=;
- cos_theta=;
- hi_cos=;
- hi_cos_square=;
+ 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);
+
+ }
+ else
+ continue;
/* end 3 body stuff (1) */