+ /* reset force */
+ v3_zero(&(atom[i].f));
+
+ /* determin cell neighbours */
+ ni=(atom[i].r.x+(moldyn->dim.x/2))/lc->x;
+ nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y;
+ nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z;
+ c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
+
+ /*
+ * processing cell of atom i
+ * => no need to check for empty list (1 element at minimum)
+ */
+ this=&(neighbour[0]);
+ list_reset(this);
+ do {
+ btom=this->current->data;
+ if(btom==&(atom[i]))
+ continue;
+
+ /* 2 body stuff */
+
+ /* we need: f_c, df_c, f_r, df_r */
+
+ v3_sub(&dist_ij,btom,&(atom[i]));
+ d_ij=v3_norm(&dist_ij);
+ if(d_ij<=S) {
+
+ /* determine the tersoff parameters */
+ if(atom[i].element!=btom->element) {
+ S=sqrt(TERSOFF_S[e1]*TERSOFF_S[e2]);
+ R=R_m;
+ A=;
+ lambda=;
+ B=;
+ mu=;
+ chi=;
+ beta=;
+ betaN=;
+
+ if(d_ij<=R) {
+ df_r=-lambda*A*exp(-lambda*d_ij)/d_ij;
+ v3_scale(&force,&dist_ij,df_r);
+ v3_add(&(atom[i].f),&(atom[i].f),
+ &force);
+ }
+ 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);
+ df_r=-lambda*f_r/d_ij;
+ scale=df_c*f_r+df_r*f_c;
+ v3_scale(&force,&dist_ij,scale);
+ v3_add(&(atom[i].f),&(atom[i].f),
+ &force);
+ }
+ }
+ 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;
+ kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
+ 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);
+ do {
+ ktom=thisk->current->data;
+ if(ktom==btom)
+ continue;
+ if(ktom==&(atom[i]))
+ 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);
+ sin_theta=sqrt(1.0/\
+ (cos_theta*cos_theta));
+ theta=arccos(cos_theta);
+
+
+ }
+ else
+ continue;
+
+ /* end 3 body stuff (1) */
+
+
+ } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+
+ /* direct neighbours of btom cell */
+ for(k=1;k<ck;k++) {
+ thisk=&(neighbourk[k]);
+ list_reset(thisk);
+ if(thisk->start!=NULL) {
+
+ do {
+ ktom=thisk->current->data;
+ if(ktom==&(atom[i]))
+ continue;
+
+ /* 3 body stuff (2) */
+
+ } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+
+ }
+ }
+
+ /* indirect neighbours of btom cell */
+ for(k=ck;k<27;k++) {
+ thisk=&(neighbourk[k]);
+ list_reset(thisk);
+ if(thisk->start!=NULL) {
+
+ do {
+ ktom=thisk->current->data;
+ if(ktom==&(atom[i]))
+ continue;
+
+ /* 3 body stuff */
+
+ } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+
+ }
+ }
+
+
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+
+ /*
+ * direct neighbour cells of atom i
+ */
+ for(j=1;j<c;j++) {
+ this=&(neighbour[j]);
+ list_reset(this);
+ if(this->start!=NULL) {
+
+ do {
+ btom=this->current->data;
+
+ /* 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;
+ kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
+ ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
+ neighbourk);
+
+ /* cell of btom */
+ thisk=&(neighbourk[0]);
+ list_reset(thisk);
+ do {
+ ktom=thisk->current->data;
+ if(ktom==btom)
+ continue;
+ if(ktom==&(atom[i]))
+ continue;
+
+ /* 3 body stuff (1) */
+
+ } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+
+ /* direct neighbours of btom cell */
+ for(k=1;k<ck;k++) {
+ thisk=&(neighbourk[k]);
+ list_reset(thisk);
+ if(thisk->start!=NULL) {
+
+ do {
+ ktom=thisk->current->data;
+ if(ktom==&(atom[i]))
+ continue;
+
+ /* 3 body stuff (2) */
+
+ } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+
+ }
+ }
+
+ /* indirect neighbours of btom cell */
+ for(k=ck;k<27;k++) {
+ thisk=&(neighbourk[k]);
+ list_reset(thisk);
+ if(thisk->start!=NULL) {
+
+ do {
+ ktom=thisk->current->data;
+ if(ktom==&(atom[i]))
+ continue;
+
+ /* 3 body stuff (3) */
+
+ } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+
+ }
+ }
+
+
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+
+ }
+ }
+
+ /*
+ * indirect neighbour cells of atom i
+ */
+ for(j=c;j<27;j++) {
+ this=&(neighbour[j]);
+ list_reset(this);
+ if(this->start!=NULL) {
+
+ do {
+ btom=this->current->data;
+
+ /* 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;
+ kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
+ ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
+ neighbourk);
+
+ /* cell of btom */
+ thisk=&(neighbourk[0]);
+ list_reset(thisk);
+ do {
+ ktom=thisk->current->data;
+ if(ktom==btom)
+ continue;
+ if(ktom==&(atom[i]))
+ continue;
+
+ /* 3 body stuff (1) */
+
+ } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+
+ /* direct neighbours of btom cell */
+ for(k=1;k<ck;k++) {
+ thisk=&(neighbourk[k]);
+ list_reset(thisk);
+ if(thisk->start!=NULL) {
+
+ do {
+ ktom=thisk->current->data;
+ if(ktom==&(atom[i]))
+ continue;
+
+ /* 3 body stuff (2) */
+
+ } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+
+ }
+ }
+
+ /* indirect neighbours of btom cell */
+ for(k=ck;k<27;k++) {
+ thisk=&(neighbourk[k]);
+ list_reset(thisk);
+ if(thisk->start!=NULL) {
+
+ do {
+ ktom=thisk->current->data;
+ if(ktom==&(atom[i]))
+ continue;
+
+ /* 3 body stuff (3) */
+
+ } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+
+ }
+ }
+
+
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+