+/* generic potential and force calculation */
+
+int potential_force_calc(t_moldyn *moldyn) {
+
+ int i,count;
+ t_atom *atom;
+ t_linkcell *lc;
+ t_list neighbour[27];
+ t_list *this;
+ double u;
+ unsigned char bc,bc3;
+ int countn,dnlc;
+
+ count=moldyn->count;
+ atom=moldyn->atom;
+ lc=&(moldyn->lc);
+
+ /* reset energy */
+ moldyn->energy=0.0;
+
+ for(i=0;i<count;i++) {
+
+ /* reset force */
+ v3_zero(&(atom[i].f));
+
+ /* single particle potential/force */
+ if(atom[i].attr&ATOM_ATTR_1BP)
+ moldyn->pf_func1b(moldyn,&(atom[i]));
+
+ /* 2 body pair potential/force */
+ if(atom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
+
+ link_cell_neighbour_index(moldyn,
+ (atom[i].r.x+moldyn->dim.x/2)/lc->x,
+ (atom[i].r.y+moldyn->dim.y/2)/lc->y,
+ (atom[i].r.z+moldyn->dim.z/2)/lc->z,
+ neighbour);
+
+ countn=lc->countn;
+ dnlc=lc->dnlc;
+
+ for(j=0;j<countn;j++) {
+
+ this=&(neighbour[j]);
+ list_reset(this);
+
+ if(this->start==NULL)
+ continue;
+
+ bc=(j<dnlc)?0:1;
+
+ do {
+ btom=this->current->data;
+
+ if(btom==&(atom[i]))
+ continue;
+
+ if((btom->attr&ATOM_ATTR_2BP)&
+ (atom[i].attr&ATOM_ATTR_2BP))
+ moldyn->pf_func2b(moldyn,
+ &(atom[i]),
+ btom,
+ bc);
+
+ /* 3 body potential/force */
+
+ if(!(atom[i].attr&ATOM_ATTR_3BP)||
+ !(btom->attr&ATOM_ATTR_3BP))
+ continue;
+
+ link_cell_neighbour_index(moldyn,
+ (btom->r.x+moldyn->dim.x/2)/lc->x,
+ (btom->r.y+moldyn->dim.y/2)/lc->y,
+ (btom->r.z+moldyn->dim.z/2)/lc->z,
+ neighbourk);
+
+ for(k=0;k<lc->countn;k++) {
+
+ thisk=&(neighbourk[k]);
+ list_reset(thisk);
+
+ if(thisk->start==NULL)
+ continue;
+
+ bck=(k<lc->dnlc)?0:1;
+
+ do {
+
+ ktom=thisk->current->data;
+
+ if(!(ktom->attr&ATOM_ATTR_3BP))
+ continue;
+
+ if(ktom==btom)
+ continue;
+
+ if(ktom==&(atom[i]))
+ continue;
+
+ moldyn->pf_func3b(moldyn,&(atom[i]),btom,ktom,bck);
+
+ } while(list_next(thisk)!=\
+ L_NO_NEXT_ELEMENT);
+
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+ }
+ }
+ }
+
+ return 0;
+}
+
+/*
+ * example potentials
+ */
+