From 814f561993377b31628e68b9c4bd59ea50da3fa3 Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Fri, 8 Jul 2022 02:30:03 +0200 Subject: [PATCH] Allow PBC in just some directions Signed-off-by: Rafael Ravedutti --- data/dem/mdbench_params.conf | 3 ++ lammps/atom.c | 54 +++++++++++----------- lammps/includes/parameter.h | 1 + lammps/parameter.c | 7 +++ lammps/pbc.c | 86 +++++++++++++++++++++--------------- 5 files changed, 91 insertions(+), 60 deletions(-) diff --git a/data/dem/mdbench_params.conf b/data/dem/mdbench_params.conf index 622c19f..d4401ad 100644 --- a/data/dem/mdbench_params.conf +++ b/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 diff --git a/lammps/atom.c b/lammps/atom.c index 27199fb..60de5df 100644 --- a/lammps/atom.c +++ b/lammps/atom.c @@ -448,32 +448,37 @@ int readAtom_in(Atom* atom, Parameter* param) { return -1; } + fgets(line, MAXLINE, fp); + natoms = atoi(strtok(line, " ")); + atom->Natoms = natoms; + atom->Nlocal = natoms; atom->ntypes = 1; - while(!feof(fp)) { + + while(atom->Nlocal >= atom->Nmax) { + growAtom(atom); + } + + for(int i = 0; i < natoms; i++) { fgets(line, MAXLINE, fp); - natoms = atoi(line); - for(int i = 0; i < natoms; i++) { - fgets(line, MAXLINE, fp); - // TODO: store mass per atom - char *s_mass = strtok(line, " "); - if(strncmp(s_mass, "inf", 3) == 0) { - // Set atom's mass to INFINITY - } else { - param->mass = atof(s_mass); - } - - atom->radius[atom_id] = atof(strtok(NULL, " ")); - atom_x(atom_id) = atof(strtok(NULL, " ")); - atom_y(atom_id) = atof(strtok(NULL, " ")); - atom_z(atom_id) = atof(strtok(NULL, " ")); - atom_vx(atom_id) = atof(strtok(NULL, " ")); - atom_vy(atom_id) = atof(strtok(NULL, " ")); - atom_vz(atom_id) = atof(strtok(NULL, " ")); - atom->type[atom_id] = 0; - atom->ntypes = MAX(atom->type[atom_id], atom->ntypes); - atom_id++; + // TODO: store mass per atom + char *s_mass = strtok(line, " "); + if(strncmp(s_mass, "inf", 3) == 0) { + // Set atom's mass to INFINITY + } else { + param->mass = atof(s_mass); } + + atom->radius[atom_id] = atof(strtok(NULL, " ")); + atom_x(atom_id) = atof(strtok(NULL, " ")); + atom_y(atom_id) = atof(strtok(NULL, " ")); + atom_z(atom_id) = atof(strtok(NULL, " ")); + atom_vx(atom_id) = atof(strtok(NULL, " ")); + atom_vy(atom_id) = atof(strtok(NULL, " ")); + atom_vz(atom_id) = atof(strtok(NULL, " ")); + atom->type[atom_id] = 0; + atom->ntypes = MAX(atom->type[atom_id], atom->ntypes); + atom_id++; } if(!natoms) { @@ -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 - atom->radius = (int *) reallocate(atom->radius, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int)); + // DEM + 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 } diff --git a/lammps/includes/parameter.h b/lammps/includes/parameter.h index 30c03bd..9a1add6 100644 --- a/lammps/includes/parameter.h +++ b/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; diff --git a/lammps/parameter.c b/lammps/parameter.c index 51c0ade..faeb7ba 100644 --- a/lammps/parameter.c +++ b/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); diff --git a/lammps/pbc.c b/lammps/pbc.c index 3610701..934dff9 100644 --- a/lammps/pbc.c +++ b/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,34 +122,54 @@ void setupPbc(Atom *atom, Parameter *param) /* Setup ghost atoms */ /* 6 planes */ - if (x < Cutneigh) { ADDGHOST(+1,0,0); } - if (x >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); } - if (y < Cutneigh) { ADDGHOST(0,+1,0); } - if (y >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); } - if (z < Cutneigh) { ADDGHOST(0,0,+1); } - if (z >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); } + 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 (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); } - if (x < Cutneigh && y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); } - if (x >= (xprd-Cutneigh) && y < Cutneigh && z < Cutneigh) { ADDGHOST(-1,+1,+1); } - 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); } + 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); } + if (x < Cutneigh && y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); } + if (x >= (xprd-Cutneigh) && y < Cutneigh && z < Cutneigh) { ADDGHOST(-1,+1,+1); } + 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 (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 (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 (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); } + 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;