+
+ /* reset force */
+ v3_zero(&(itom[i].f));
+
+ /* single particle potential/force */
+ if(itom[i].attr&ATOM_ATTR_1BP)
+ moldyn->func1b(moldyn,&(itom[i]));
+
+ /* 2 body pair potential/force */
+ if(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
+
+ link_cell_neighbour_index(moldyn,
+ (itom[i].r.x+moldyn->dim.x/2)/lc->x,
+ (itom[i].r.y+moldyn->dim.y/2)/lc->y,
+ (itom[i].r.z+moldyn->dim.z/2)/lc->z,
+ neighbour_i);
+
+ countn=lc->countn;
+ dnlc=lc->dnlc;
+
+ for(j=0;j<countn;j++) {
+
+ this=&(neighbour_i[j]);
+ list_reset(this);
+
+ if(this->start==NULL)
+ continue;
+
+ bc_ij=(j<dnlc)?0:1;
+
+ do {
+ jtom=this->current->data;
+
+ if(jtom==&(itom[i]))
+ continue;
+
+ if((jtom->attr&ATOM_ATTR_2BP)&
+ (itom[i].attr&ATOM_ATTR_2BP))
+ moldyn->func2b(moldyn,
+ &(itom[i]),
+ jtom,
+ bc_ij);
+
+ /* 3 body potential/force */
+
+ if(!(itom[i].attr&ATOM_ATTR_3BP)||
+ !(jtom->attr&ATOM_ATTR_3BP))
+ continue;
+
+ /*
+ * according to mr. nordlund, we dont need to take the
+ * sum over all atoms now, as 'this is centered' around
+ * atom i ...
+ * i am not quite sure though! there is a not vanishing
+ * part even if f_c_ik is zero ...
+ * this analytical potentials suck!
+ * switching from mc to md to dft soon!
+ */
+
+ // link_cell_neighbour_index(moldyn,
+ // (jtom->r.x+moldyn->dim.x/2)/lc->x,
+ // (jtom->r.y+moldyn->dim.y/2)/lc->y,
+ // (jtom->r.z+moldyn->dim.z/2)/lc->z,
+ // neighbour_j);
+
+// /* neighbours of j */
+// for(k=0;k<lc->countn;k++) {
+//
+// that=&(neighbour_j[k]);
+// list_reset(that);
+//
+// if(that->start==NULL)
+// continue;
+//
+// bc_ijk=(k<lc->dnlc)?0:1;
+//
+// do {
+//
+// ktom=that->current->data;
+//
+// if(!(ktom->attr&ATOM_ATTR_3BP))
+// continue;
+//
+// if(ktom==jtom)
+// continue;
+//
+// if(ktom==&(itom[i]))
+// continue;
+//
+// moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
+//
+/* } while(list_next(that)!=\ */
+// L_NO_NEXT_ELEMENT);
+//
+// }
+
+ /* copy the neighbour lists */
+ memcpy(neighbour_i2,neighbour_i,
+ 27*sizeof(t_list));
+
+ /* get neighbours of i */
+ for(k=0;k<countn;k++) {
+
+ that=&(neighbour_i2[k]);
+ list_reset(that);
+
+ if(that->start==NULL)
+ continue;
+
+ bc_ijk=(k<dnlc)?0:1;
+
+ do {
+
+ ktom=that->current->data;
+
+ if(!(ktom->attr&ATOM_ATTR_3BP))
+ continue;
+
+ if(ktom==jtom)
+ continue;
+
+ if(ktom==&(itom[i]))
+ continue;
+
+printf("Debug: atom %d before 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
+ moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
+printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
+
+ } while(list_next(that)!=\
+ L_NO_NEXT_ELEMENT);
+
+ }
+
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+ }