lc=&(moldyn->lc);
count=moldyn->count;
- /* reset enrgy counter */
+ /* reset energy counter */
u=0.0;
for(i=0;i<count;i++) {
/* 2 body stuff */
+ v3_sub(&dist_ij,btom,&(atom[i]));
+ d_ij=v3_norm(&dist_ij);
+ if(d_ij<=S) {
+
+ S=;
+ R=;
+ A=;
+ lambda=;
+ B=;
+ mu=;
+ chi=;
+ beta=;
+ betaN=;
+
+ if(d_ij<=R) {
+ f_c=1.0;
+ df_c=0.0;
+ }
+ else {
+ s_r=S-R;
+ arg1=PI*(d_ij-R)/s_r;
+ f_c=0.5+0.5*cos(arg1);
+ df_c=-0.5*sin(arg1)*(PI/(s_r*d_ij));
+ f_r=A*exp(-lambda*d_ij);
+ f_a=-B*exp(-mu*d_ij);
+ }
+ }
+ else
+ continue;
+
+
+ /* end 2 body stuff */
+
/* determine cell neighbours of btom */
ki=(btom->r.x+(moldyn->dim.x/2))/lc->x;
kj=(btom->r.y+(moldyn->dim.y/2))/lc->y;
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);
continue;
/* 3 body stuff (1) */
+
+ 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) */
+
} while(list_next(thisk)!=L_NO_NEXT_ELEMENT);