Allow PBC in just some directions

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>

data/dem/mdbench_params.conf

@@ -1,3 +1,6 @@
+pbc_x 0
+pbc_y 0
+pbc_z 0
 ntimes 200
 dt 0.001
 reneigh_every 20

lammps/atom.c

@@ -448,10 +448,16 @@ int readAtom_in(Atom* atom, Parameter* param) {
         return -1;
     }
 
-    atom->ntypes = 1;
-    while(!feof(fp)) {
     fgets(line, MAXLINE, fp);
-        natoms = atoi(line);
+    natoms = atoi(strtok(line, " "));
+    atom->Natoms = natoms;
+    atom->Nlocal = natoms;
+    atom->ntypes = 1;
+
+    while(atom->Nlocal >= atom->Nmax) {
+        growAtom(atom);
+    }
+
     for(int i = 0; i < natoms; i++) {
         fgets(line, MAXLINE, fp);
 
@@ -474,7 +480,6 @@ int readAtom_in(Atom* atom, Parameter* param) {
         atom->ntypes = MAX(atom->type[atom_id], atom->ntypes);
         atom_id++;
     }
-    }
 
     if(!natoms) {
         fprintf(stderr, "Input error: atom data was not read!\n");
@@ -517,9 +522,8 @@ void growAtom(Atom *atom) {
     atom->fz = (MD_FLOAT*) reallocate(atom->fz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
     #endif
     atom->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int));
-    #ifdef DEM
+    // DEM
-    atom->radius = (int *) reallocate(atom->radius, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int));
+    atom->radius = (MD_FLOAT *) reallocate(atom->radius, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
     atom->av = (MD_FLOAT*) reallocate(atom->av, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 3, nold * sizeof(MD_FLOAT) * 3);
     atom->r  = (MD_FLOAT*) reallocate(atom->r,  ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 4, nold * sizeof(MD_FLOAT) * 4);
-    #endif
 }

lammps/includes/parameter.h

@@ -53,6 +53,7 @@ typedef struct {
     MD_FLOAT cutforce;
     MD_FLOAT cutneigh;
     int nx, ny, nz;
+    int pbc_x, pbc_y, pbc_z;
     MD_FLOAT lattice;
     MD_FLOAT xlo, xhi, ylo, yhi, zlo, zhi;
     MD_FLOAT xprd, yprd, zprd;

lammps/parameter.c

@@ -42,6 +42,9 @@ void initParameter(Parameter *param) {
     param->nx = 32;
     param->ny = 32;
     param->nz = 32;
+    param->pbc_x = 1;
+    param->pbc_y = 1;
+    param->pbc_z = 1;
     param->cutforce = 2.5;
     param->skin = 0.3;
     param->cutneigh = param->cutforce + param->skin;
@@ -116,6 +119,9 @@ void readParameter(Parameter *param, const char *filename) {
             PARSE_INT(nx);
             PARSE_INT(ny);
             PARSE_INT(nz);
+            PARSE_INT(pbc_x);
+            PARSE_INT(pbc_y);
+            PARSE_INT(pbc_z);
             PARSE_INT(nstat);
             PARSE_INT(reneigh_every);
             PARSE_INT(x_out_every);
@@ -147,6 +153,7 @@ void printParameter(Parameter *param) {
     printf("\tForce field: %s\n", ff2str(param->force_field));
     printf("\tUnit cells (nx, ny, nz): %d, %d, %d\n", param->nx, param->ny, param->nz);
     printf("\tDomain box sizes (x, y, z): %e, %e, %e\n", param->xprd, param->yprd, param->zprd);
+    printf("\tPeriodic (x, y, z): %d, %d, %d\n", param->pbc_x, param->pbc_y, param->pbc_z);
     printf("\tLattice size: %e\n", param->lattice);
     printf("\tEpsilon: %e\n", param->epsilon);
     printf("\tSigma: %e\n", param->sigma);

lammps/pbc.c

@@ -35,8 +35,7 @@ static int *PBCx, *PBCy, *PBCz;
 static void growPbc(Atom*);
 
 /* exported subroutines */
-void initPbc(Atom* atom)
+void initPbc(Atom* atom) {
-{
     NmaxGhost = 0;
     atom->border_map = NULL;
     PBCx = NULL; PBCy = NULL; PBCz = NULL;
@@ -44,8 +43,7 @@ void initPbc(Atom* atom)
 
 /* update coordinates of ghost atoms */
 /* uses mapping created in setupPbc */
-void updatePbc(Atom *atom, Parameter *param)
+void updatePbc(Atom *atom, Parameter *param) {
-{
     int *border_map = atom->border_map;
     int nlocal = atom->Nlocal;
     MD_FLOAT xprd = param->xprd;
@@ -61,8 +59,7 @@ void updatePbc(Atom *atom, Parameter *param)
 
 /* relocate atoms that have left domain according
  * to periodic boundary conditions */
-void updateAtomsPbc(Atom *atom, Parameter *param)
+void updateAtomsPbc(Atom *atom, Parameter *param) {
-{
     MD_FLOAT xprd = param->xprd;
     MD_FLOAT yprd = param->yprd;
     MD_FLOAT zprd = param->zprd;
@@ -101,8 +98,7 @@ void updateAtomsPbc(Atom *atom, Parameter *param)
     PBCz[Nghost] = dz;                                  \
     atom->type[atom->Nlocal + Nghost] = atom->type[i]
 
-void setupPbc(Atom *atom, Parameter *param)
+void setupPbc(Atom *atom, Parameter *param) {
-{
     int *border_map = atom->border_map;
     MD_FLOAT xprd = param->xprd;
     MD_FLOAT yprd = param->yprd;
@@ -111,10 +107,10 @@ void setupPbc(Atom *atom, Parameter *param)
     int Nghost = -1;
 
     for(int i = 0; i < atom->Nlocal; i++) {
-
         if (atom->Nlocal + Nghost + 7 >= atom->Nmax) {
             growAtom(atom);
         }
+
         if (Nghost + 7 >= NmaxGhost) {
             growPbc(atom);
             border_map = atom->border_map;
@@ -126,13 +122,23 @@ void setupPbc(Atom *atom, Parameter *param)
 
         /* Setup ghost atoms */
         /* 6 planes */
+        if(param->pbc_x != 0) {
             if (x < Cutneigh)         { ADDGHOST(+1,0,0); }
             if (x >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); }
+        }
+
+        if(param->pbc_y != 0) {
             if (y < Cutneigh)         { ADDGHOST(0,+1,0); }
             if (y >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); }
+        }
+
+        if(param->pbc_z != 0) {
             if (z < Cutneigh)         { ADDGHOST(0,0,+1); }
             if (z >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); }
+        }
+
         /* 8 corners */
+        if(param->pbc_x != 0 && param->pbc_y != 0 && param->pbc_z != 0) {
             if (x < Cutneigh         && y < Cutneigh         && z < Cutneigh)         { ADDGHOST(+1,+1,+1); }
             if (x < Cutneigh         && y >= (yprd-Cutneigh) && z < Cutneigh)         { ADDGHOST(+1,-1,+1); }
             if (x < Cutneigh         && y >= Cutneigh        && z >= (zprd-Cutneigh)) { ADDGHOST(+1,+1,-1); }
@@ -141,20 +147,30 @@ void setupPbc(Atom *atom, Parameter *param)
             if (x >= (xprd-Cutneigh) && y >= (yprd-Cutneigh) && z < Cutneigh)         { ADDGHOST(-1,-1,+1); }
             if (x >= (xprd-Cutneigh) && y < Cutneigh         && z >= (zprd-Cutneigh)) { ADDGHOST(-1,+1,-1); }
             if (x >= (xprd-Cutneigh) && y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(-1,-1,-1); }
+        }
+
         /* 12 edges */
+        if(param->pbc_x != 0 && param->pbc_z != 0) {
             if (x < Cutneigh         && z < Cutneigh)         { ADDGHOST(+1,0,+1); }
             if (x < Cutneigh         && z >= (zprd-Cutneigh)) { ADDGHOST(+1,0,-1); }
             if (x >= (xprd-Cutneigh) && z < Cutneigh)         { ADDGHOST(-1,0,+1); }
             if (x >= (xprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(-1,0,-1); }
+        }
+
+        if(param->pbc_y != 0 && param->pbc_z != 0) {
             if (y < Cutneigh         && z < Cutneigh)         { ADDGHOST(0,+1,+1); }
             if (y < Cutneigh         && z >= (zprd-Cutneigh)) { ADDGHOST(0,+1,-1); }
             if (y >= (yprd-Cutneigh) && z < Cutneigh)         { ADDGHOST(0,-1,+1); }
             if (y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(0,-1,-1); }
+        }
+
+        if(param->pbc_x != 0 && param->pbc_y != 0) {
             if (y < Cutneigh         && x < Cutneigh)         { ADDGHOST(+1,+1,0); }
             if (y < Cutneigh         && x >= (xprd-Cutneigh)) { ADDGHOST(-1,+1,0); }
             if (y >= (yprd-Cutneigh) && x < Cutneigh)         { ADDGHOST(+1,-1,0); }
             if (y >= (yprd-Cutneigh) && x >= (xprd-Cutneigh)) { ADDGHOST(-1,-1,0); }
         }
+    }
     // increase by one to make it the ghost atom count
     atom->Nghost = Nghost + 1;
 }