security ci

diff --git a/moldyn.c b/moldyn.c
index 57f76c4..d74b391 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -694,14 +694,6 @@ printf("done\n");
 
 /* generic potential and force calculation */
 
-#define                CREATE_CELL_LIST(nb_list) \
-       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,\
-                                 nb_list);
-
-
 int potential_force_calc(t_moldyn *moldyn) {
 
        int i,count;
@@ -710,13 +702,15 @@ int potential_force_calc(t_moldyn *moldyn) {
        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 */
-       u=0.0;
+       moldyn->energy=0.0;
 
        for(i=0;i<count;i++) {
        
@@ -724,71 +718,99 @@ int potential_force_calc(t_moldyn *moldyn) {
                v3_zero(&(atom[i].f));
 
                /* single particle potential/force */
-               if(moldyn->status&MOLDYN_STAT_1BP)
+               if(atom[i].attr&ATOM_ATTR_1BP)
                        moldyn->pf_func1b(moldyn,&(atom[i]));
 
                /* 2 body pair potential/force */
-               if(moldyn->status&MOLDYN_STAT_2BP) {
+               if(atom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
                
-                       CREATE_CELL_LIST(neighbour);
+                       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);
 
-                       /*
-                        * processing cell of atom i
-                        * => no need to check for empty list
-                        *    (1 element at minimum)
-                        */
+                       countn=lc->countn;
+                       dnlc=lc->dnlc;
 
-                       this=&(neighbour[0]);
-                       list_reset(this);
-                       do {
-                               btom=this->current->data;
-                               if(btom!=&(atom[i]))
-                                       moldyn->pf_func2b(moldyn,
-                                                         &(atom[i]),btom);
-                       } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+                       for(j=0;j<countn;j++) {
 
-                       /*
-                        * direct neighbour cells
-                        * => no boundary condition check necessary
-                        */
-                       for(j=0;j<lc->dnlc;j++) {
                                this=&(neighbour[j]);
                                list_reset(this);
-                               if(this->start!=NULL) {
-                                       do {
-                                               btom=this->current->data;
-                                               moldyn->pf_func2b(moldyn,
-                                                                 &(atom[i]),
-                                                                 btom);
-                                       } while(list_next(this)!=\
-                                               L_NO_NEXT_ELEMENT);
-                       }
 
-                       /*
-                        * neighbour cells due to periodic bc
-                        * => check boundary conditions
-                        */
-                       for(j=lc->dnlc;j<lc->countn;j++) {
-                               this=&(neighbour[j]);
-                               list_reset(this);
-                               if(this->start!=NULL) {
-                                       do {
-                                               btom=this->current->data;
+                               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);
+                                                                 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
+ */
 
 /* harmonic oscillator potential and force */
 
-int harmonic_oscillator(t_moldyn *moldyn) {
+int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,unsigned char bc)) {
 
        t_ho_params *params;
        t_atom *atom,*btom;

diff --git a/moldyn.h b/moldyn.h
index 717388a..66f774e 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -15,6 +15,7 @@
 
 /* datatypes */
 
+typedef unsigned char u8;
 
 /* the atom of the md simulation */
 
@@ -30,6 +31,10 @@ typedef struct s_atom {
 #define ATOM_ATTR_FP   0x01    /* fixed position (bulk material) */
 #define ATOM_ATTR_HB   0x02    /* coupled to heat bath (velocity scaling) */
 
+#define ATOM_ATTR_1BP  0x10    /* single paricle potential */
+#define ATOM_ATTR_2BP  0x20    /* pair potential */
+#define ATOM_ATTR_3BP  0x40    /* 3 body potential */ 
+
 /* cell lists */
 
 typedef struct s_linkcell {
@@ -99,9 +104,6 @@ typedef struct s_moldyn {
 #define MOLDYN_LVSTAT_VISUAL           0x08
 #define MOLDYN_LVSTAT_INITIALIZED      0x10
 
-#define MOLDYN_STAT_1BP                        0x01
-#define MOLDYN_STAT_2BP                        0x02    /* define for pair potentials */
-#define MOLDYN_STAT_3BP                        0x04    /* define for 3 body pot */ 
 #define MOLDYN_STAT_PBX                        0x08    /* periodic boudaries in x */
 #define MOLDYN_STAT_PBY                        0x10    /* y */
 #define MOLDYN_STAT_PBZ                        0x20    /* and z direction */