Refactor. Fix bug in atom initialization.

include_CLANG.mk

@@ -7,7 +7,7 @@ ANSI_CFLAGS += -std=c99
 ANSI_CFLAGS += -pedantic
 ANSI_CFLAGS += -Wextra
 
-CFLAGS   = -O3 #-g  $(ANSI_CFLAGS)
+CFLAGS   = -O3 $(ANSI_CFLAGS) #-g
 ASFLAGS  = -masm=intel
 CXXFLAGS = $(CFLAGS)
 FCFLAGS  =

src/allocate.c

@@ -64,8 +64,7 @@ void* reallocate (
         size_t newBytesize,
         size_t oldBytesize)
 {
-    void* newarray;
+    void* newarray =  allocate(alignment, newBytesize);
-    newarray =  allocate(alignment, newBytesize);
 
     if(ptr != NULL) {
         memcpy(newarray, ptr, oldBytesize);

src/atom.c

@@ -24,15 +24,27 @@
  *
  * =======================================================================================
  */
+#include <stdlib.h>
+#include <stdio.h>
 #include <math.h>
 
-#include <parameter.h>
-#include <allocate.h>
 #include <atom.h>
+#include <allocate.h>
 #include <util.h>
 
 #define DELTA 20000
 
+void initAtom(Atom *atom)
+{
+    atom->x  = NULL; atom->y  = NULL; atom->z  = NULL;
+    atom->vx = NULL; atom->vy = NULL; atom->vz = NULL;
+    atom->fx = NULL; atom->fy = NULL; atom->fz = NULL;
+    atom->Natoms = 0;
+    atom->Nlocal = 0;
+    atom->Nghost = 0;
+    atom->Nmax   = 0;
+}
+
 void createAtom(Atom *atom, Parameter *param)
 {
     double xlo = 0.0; double xhi = param->xprd;

src/includes/atom.h

@@ -36,6 +36,7 @@ typedef struct {
     double *fx, *fy, *fz;
 } Atom;
 
+extern void initAtom(Atom*);
 extern void createAtom(Atom*, Parameter*);
 extern void growAtom(Atom*);
 #endif

src/includes/neighbor.h

@@ -39,7 +39,7 @@ typedef struct {
 } Neighbor;
 
 extern void initNeighbor(Neighbor*, Parameter*);
-extern void setupNeighbor(Parameter*);
+extern void setupNeighbor();
 extern void binatoms(Atom*);
 extern void buildNeighbor(Atom*, Neighbor*);
 #endif

src/main.c

@@ -33,6 +33,7 @@
 #include <math.h>
 #include <float.h>
 
+#include <timing.h>
 #include <allocate.h>
 #include <neighbor.h>
 #include <parameter.h>
@@ -42,16 +43,16 @@
 
 #define HLINE "----------------------------------------------------------------------------\n"
 
-#define FACTOR 0.999
+typedef enum {
-#define SMALL 1.0e-6
+    TOTAL = 0,
-#define DELTA 20000
+    NEIGH,
+    FORCE,
+    NUMTIMER
+} timertype;
 
-void init(Atom *atom, Parameter *param)
+
+void init(Parameter *param)
 {
-    atom->x = NULL; atom->y = NULL; atom->z = NULL;
-    atom->vx = NULL; atom->vy = NULL; atom->vz = NULL;
-    atom->fx = NULL; atom->fy = NULL; atom->fz = NULL;
-
     param->epsilon = 1.0;
     param->sigma6 = 1.0;
     param->rho = 0.8442;
@@ -73,6 +74,19 @@ void init(Atom *atom, Parameter *param)
     param->zprd = param->nz * lattice;
 }
 
+void setup(Parameter *param, Atom *atom, Neighbor *neighbor){
+    initAtom(atom);
+    initNeighbor(neighbor, param);
+    initPbc();
+    setupNeighbor();
+    createAtom(atom, param);
+    setupThermo(param, atom->Natoms);
+    adjustThermo(param, atom);
+    setupPbc(atom, param);
+    updatePbc(atom, param);
+    buildNeighbor(atom, neighbor);
+}
+
 void initialIntegrate(Parameter *param, Atom *atom)
 {
     double* x = atom->x; double* y = atom->y; double* z = atom->z;
@@ -101,7 +115,7 @@ void finalIntegrate(Parameter *param, Atom *atom)
     }
 }
 
-void computeForce(Neighbor *neighbor, Atom *atom, Parameter *param)
+void computeForce(Parameter *param, Atom *atom, Neighbor *neighbor)
 {
     int Nlocal = atom->Natoms;
     int* neighs;
@@ -164,26 +178,52 @@ void printAtomState(Atom *atom)
     }
 }
 
-
 int main (int argc, char** argv)
 {
+    double timer[NUMTIMER];
+    double S;
     Atom atom;
     Neighbor neighbor;
     Parameter param;
 
-    init(&atom, &param);
+    init(&param);
-    initNeighbor(&neighbor, &param);
+
-    initPbc();
+    for(int i = 0; i < argc; i++)
-    setupNeighbor(&param);
+    {
-    createAtom(&atom, &param);
+        if((strcmp(argv[i], "-n") == 0) || (strcmp(argv[i], "--nsteps") == 0))
-    setupThermo(&param, atom.Natoms);
+        {
-    adjustThermo(&param, &atom);
+            param.ntimes = atoi(argv[++i]);
-    /* printAtomState(&atom); */
+            continue;
-    setupPbc(&atom, &param);
+        }
-    updatePbc(&atom, &param);
+        if((strcmp(argv[i], "-nx") == 0))
-    buildNeighbor(&atom, &neighbor);
+        {
+            param.nx = atoi(argv[++i]);
+            continue;
+        }
+        if((strcmp(argv[i], "-ny") == 0))
+        {
+            param.ny = atoi(argv[++i]);
+            continue;
+        }
+        if((strcmp(argv[i], "-nz") == 0))
+        {
+            param.nz = atoi(argv[++i]);
+            continue;
+        }
+        if((strcmp(argv[i], "-h") == 0) || (strcmp(argv[i], "--help") == 0))
+        {
+            printf("MD Bench: A minimalistic re-implementation of miniMD\n");
+            printf(HLINE);
+            printf("-n / --nsteps <int>:    set number of timesteps for simulation\n");
+            printf("-nx/-ny/-nz <int>:      set linear dimension of systembox in x/y/z direction\n");
+            printf(HLINE);
+            exit(EXIT_SUCCESS);
+        }
+    }
+
+    setup(&param, &atom, &neighbor);
     computeThermo(0, &param, &atom);
-    computeForce(&neighbor, &atom, &param);
+    computeForce(&param, &atom, &neighbor);
 
     for(int n = 0; n < param.ntimes; n++) {
 
@@ -199,7 +239,7 @@ int main (int argc, char** argv)
             buildNeighbor(&atom, &neighbor);
         }
 
-        computeForce(&neighbor, &atom, &param);
+        computeForce(&param, &atom, &neighbor);
         finalIntegrate(&param, &atom);
 
         if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) {

src/neighbor.c

@@ -76,7 +76,7 @@ void initNeighbor(Neighbor *neighbor, Parameter *param)
     neighbor->neighbors = NULL;
 }
 
-void setupNeighbor(Parameter *param)
+void setupNeighbor()
 {
     double coord;
     int mbinxhi, mbinyhi, mbinzhi;

src/pbc.c

@@ -39,17 +39,7 @@ static int *PBCx, *PBCy, *PBCz;
 
 static void growPbc();
 
-void growPbc()
+/* exported subroutines */
-{
-    int nold = NmaxGhost;
-    NmaxGhost += DELTA;
-
-    BorderMap = (int*) reallocate(BorderMap, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
-    PBCx = (int*) reallocate(PBCx, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
-    PBCy = (int*) reallocate(PBCy, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
-    PBCz = (int*) reallocate(PBCz, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
-}
-
 void initPbc()
 {
     NmaxGhost = 0;
@@ -57,6 +47,8 @@ void initPbc()
     PBCx = NULL; PBCy = NULL; PBCz = NULL;
 }
 
+/* update coordinates of ghost atoms */
+/* uses mapping created in setupPbc */
 void updatePbc(Atom *atom, Parameter *param)
 {
     int nlocal = atom->Nlocal;
@@ -74,6 +66,8 @@ void updatePbc(Atom *atom, Parameter *param)
     }
 }
 
+/* relocate atoms that have left domain according
+ * to periodic boundary conditions */
 void updateAtomsPbc(Atom *atom, Parameter *param)
 {
     double* x = atom->x;
@@ -105,6 +99,10 @@ void updateAtomsPbc(Atom *atom, Parameter *param)
     }
 }
 
+/* setup periodic boundary conditions by
+ * defining ghost atoms around domain
+ * only creates mapping and coordinate corrections
+ * that are then enforced in updatePbc */
 #define ADDGHOST(dx,dy,dz) Nghost++; BorderMap[Nghost] = i; PBCx[Nghost] = dx; PBCy[Nghost] = dy; PBCz[Nghost] = dz;
 void setupPbc(Atom *atom, Parameter *param)
 {
@@ -159,3 +157,15 @@ void setupPbc(Atom *atom, Parameter *param)
     // increase by one to make it the ghost atom count
     atom->Nghost = Nghost + 1;
 }
+
+/* internal subroutines */
+void growPbc()
+{
+    int nold = NmaxGhost;
+    NmaxGhost += DELTA;
+
+    BorderMap = (int*) reallocate(BorderMap, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
+    PBCx = (int*) reallocate(PBCx, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
+    PBCy = (int*) reallocate(PBCy, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
+    PBCz = (int*) reallocate(PBCz, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
+}

src/thermo.c

@@ -44,6 +44,7 @@ static double p_act;
 static double e_act;
 static int mstat;
 
+/* exported subroutines */
 void setupThermo(Parameter *param, int natoms)
 {
     int maxstat = param->ntimes / param->nstat + 2;

util/mdBench.c

@@ -325,7 +325,7 @@ void init(Neighbor *neighbor, Parameter *param)
     param->dt = 0.005;
     param->nx = 32;
     param->ny = 32;
-    param->nz = 32;
+    param->nz = 64;
     param->cutforce = 2.5;
     param->temp = 1.44;
     param->nstat = 100;