+/* tersoff potential & force for 2 sorts of atoms */
+
+int tersoff(t_moldyn *moldyn) {
+
+ t_tersoff_params *params;
+ t_atom *atom,*btom,*ktom;
+ t_linkcell *lc;
+ t_list *this,*thisk,neighbour[27],neighbourk[27];
+ int i,j,k,c,ck;
+ int count;
+ double u;
+ int ni,nj,nk;
+ int ki,kj,kk;
+
+
+ params=moldyn->pot_params;
+ atom=moldyn->atom;
+ lc=&(moldyn->lc);
+ count=moldyn->count;
+
+ /* reset enrgy counter */
+ u=0.0;
+
+ for(i=0;i<count;i++) {
+ /* 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 */
+
+ /* 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 */
+
+ } 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);
+
+ }
+ }
+
+ }
+
+ moldyn->energy=0.5*u;
+
+ return 0;
+}
+