ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

clean up of docs and code

Browse Source
This commit is contained in:
Steve Plimpton
2024-05-06 19:16:06 -06:00
parent 1c11de8a64
commit 3028b6f34c
4 changed files with 96 additions and 558 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
doc/src/replicate.rst
View File

@ -8,7 +8,7 @@ Syntax
.. code-block:: LAMMPS .. code-block:: LAMMPS
replicate nx ny nz *keyword* ... replicate nx ny nz keyword ...
nx,ny,nz = replication factors in each dimension nx,ny,nz = replication factors in each dimension
@ -17,8 +17,8 @@ nx,ny,nz = replication factors in each dimension
.. parsed-literal:: .. parsed-literal::
*bbox* = only check atoms in replicas that overlap with a processor's subdomain *bbox* = use a bounding-box algorithm which is faster for large proc counts
*bond/periodic* = use a different algorithm that correctly replicates periodic bond loops *bond/periodic* = use an algorithm that correctly replicates periodic bond loops
Examples Examples
"""""""" """"""""
@ -56,7 +56,7 @@ are created between pairs of new atoms as well as between old and new
atoms. atoms.
.. note:: .. note::
The bond discussion which follows only refers to models with The bond discussion which follows only refers to models with
permanent covalent bonds typically defined in LAMMPS via a data permanent covalent bonds typically defined in LAMMPS via a data
file. It is not relevant to sytems modeled with many-body file. It is not relevant to sytems modeled with many-body
@ -83,7 +83,7 @@ atoms and uses a different algorithm to find new (nearby) bond
neighbors in the replicated system. In the final replicated system neighbors in the replicated system. In the final replicated system
all image flags are zero (in each dimension). all image flags are zero (in each dimension).
-- note: .. note::
LAMMPS does not check for image flag consistency before performing LAMMPS does not check for image flag consistency before performing
the replication (it does issue a warning about this before a the replication (it does issue a warning about this before a
@ -92,13 +92,13 @@ all image flags are zero (in each dimension).
will otherwise be correctly replicated. This is NOT the case if will otherwise be correctly replicated. This is NOT the case if
there is a periodic bond loop. See the next note. there is a periodic bond loop. See the next note.
-- note: .. note::
LAMMPS does not check for periodic bond loops. If you use the LAMMPS does not check for periodic bond loops. If you use the
*bond/periodic* option for a system without periodic bond loops, *bond/periodic* keyword for a system without periodic bond loops,
the system will be correctly replicated, but image flag information the system will be correctly replicated, but image flag information
will be lost (which may or may not be important to your model). If will be lost (which may or may not be important to your model). If
you do not use the *bond/periodic* option for a system with you do not use the *bond/periodic* keyword for a system with
periodic bond loops, the replicated system will have invalid bonds periodic bond loops, the replicated system will have invalid bonds
(typically very long), resulting in bad dynamics. (typically very long), resulting in bad dynamics.
@ -112,7 +112,7 @@ requires a temporary use of more memory. Each processor must be able
to store all atoms (and their per-atom data) in the original system, to store all atoms (and their per-atom data) in the original system,
before it is replicated. before it is replicated.
-- note: .. note::
The algorithm used by the *bond/periodic* keyword builds on the The algorithm used by the *bond/periodic* keyword builds on the
algorithm used by the *bbox* keyword and thus has the same memory algorithm used by the *bbox* keyword and thus has the same memory
@ -121,7 +121,7 @@ before it is replicated.
---------- ----------
Restrictions Restrictions
"""""""""""" """"""""""""
A 2d simulation cannot be replicated in the z dimension. A 2d simulation cannot be replicated in the z dimension.

9
examples/replicate/README
View File

@ -7,6 +7,8 @@ cross the periodic boundary to close the loop.
To run these scripts, LAMMPS should be built with the MOLECULE and To run these scripts, LAMMPS should be built with the MOLECULE and
CLASS2 packages. The latter is only needed for the CNT example. CLASS2 packages. The latter is only needed for the CNT example.
--------
These scripts are tiny examples which illustrate both kinds of These scripts are tiny examples which illustrate both kinds of
systems. Each produces a series of images which can be visualized. systems. Each produces a series of images which can be visualized.
If the 3 lines for a dump movie command are uncommented, a MPG movie If the 3 lines for a dump movie command are uncommented, a MPG movie
@ -17,6 +19,13 @@ in.replcate.bond.x.y # 2d grid of bonded atoms, bond loops in x and y
in.replicate.bond.xy # linear chains in diagonal direction, bond loop in x and y in.replicate.bond.xy # linear chains in diagonal direction, bond loop in x and y
in.replicate.bond.noloop # linear chains in x direction, no bond loop in.replicate.bond.noloop # linear chains in x direction, no bond loop
If you do not use the bond/periodic keyword with the replicate command
in the first 3 of these scripts (which have periodic bond loops), and
visualize the dynamics of hee simulation, you will see how the
replication creates a bogus system.
--------
This script is for a complex system of 3 orthogonal CNTs which has This script is for a complex system of 3 orthogonal CNTs which has
periodic bond loops in all 3 dimensions xyz. periodic bond loops in all 3 dimensions xyz.

621
src/replicate.cpp
View File

@ -13,7 +13,7 @@
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
Contributing authors: Contributing authors:
Chris Knight (ANL) for bbox option Chris Knight (ANL) for bbox option
Jake Gissinger (Stevens Institute of Technology) for bond/periodic option Jake Gissinger (Stevens Institute of Technology) for bond/periodic option
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -77,7 +77,7 @@ void Replicate::command(int narg, char **arg)
nx, ny, nz, nrep); nx, ny, nz, nrep);
// optional keywords // optional keywords
bbox_flag = 0; bbox_flag = 0;
bond_flag = 0; bond_flag = 0;
@ -134,9 +134,9 @@ void Replicate::command(int narg, char **arg)
MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world); MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
maxmol = maxmol_all; maxmol = maxmol_all;
} }
// reset image flags to zero for bond/periodic option // reset image flags to zero for bond/periodic option
if (bond_flag) if (bond_flag)
for (i=0; i<atom->nlocal; ++i) for (i=0; i<atom->nlocal; ++i)
atom->image[i] = ((imageint) IMGMAX << IMG2BITS) | atom->image[i] = ((imageint) IMGMAX << IMG2BITS) |
@ -178,7 +178,7 @@ void Replicate::command(int narg, char **arg)
for (i = 0; i < atom->nlocal; i++) for (i = 0; i < atom->nlocal; i++)
domain->unmap(atom->x[i],atom->image[i]); domain->unmap(atom->x[i],atom->image[i]);
// communication buffer for all my atom's info // communication buffer for all my atom's info
// max_size = largest buffer needed by any proc // max_size = largest buffer needed by any proc
// must do before new Atom class created, since size_restart() uses atom->nlocal // must do before new Atom class created, since size_restart() uses atom->nlocal
@ -378,485 +378,14 @@ void Replicate::command(int narg, char **arg)
} }
} }
// use // use one of two algorithms for replication
if (!bbox_flag) { if (!bbox_flag) {
replicate_by_proc(nx,ny,nz,sublo,subhi,buf); replicate_by_proc(nx,ny,nz,sublo,subhi,buf);
} else { } else {
replicate_by_bbox(nx,ny,nz,sublo,subhi,buf); replicate_by_bbox(nx,ny,nz,sublo,subhi,buf);
} }
/*
AtomVec *old_avec = old->avec;
AtomVec *avec = atom->avec;
int ix,iy,iz;
tagint atom_offset,mol_offset,atom0tag;
imageint image;
double x[3],lamda[3];
double *coord;
int tag_enable = atom->tag_enable;
if (bbox_flag || bond_flag) {
// allgather size of buf on each proc
n = 0;
for (i = 0; i < old->nlocal; i++) n += old_avec->pack_restart(i,&buf[n]);
int * size_buf_rnk;
memory->create(size_buf_rnk, nprocs, "replicate:size_buf_rnk");
MPI_Allgather(&n, 1, MPI_INT, size_buf_rnk, 1, MPI_INT, world);
// size of buf_all
int size_buf_all = 0;
MPI_Allreduce(&n, &size_buf_all, 1, MPI_INT, MPI_SUM, world);
if (me == 0) {
auto mesg = fmt::format(" bounding box image = ({} {} {}) "
"to ({} {} {})\n",
_imagelo[0],_imagelo[1],_imagelo[2],
_imagehi[0],_imagehi[1],_imagehi[2]);
mesg += fmt::format(" bounding box extra memory = {:.2f} MB\n",
(double)size_buf_all*sizeof(double)/1024/1024);
utils::logmesg(lmp,mesg);
}
// rnk offsets
int *disp_buf_rnk;
memory->create(disp_buf_rnk, nprocs, "replicate:disp_buf_rnk");
disp_buf_rnk[0] = 0;
for (i = 1; i < nprocs; i++)
disp_buf_rnk[i] = disp_buf_rnk[i-1] + size_buf_rnk[i-1];
// allgather buf_all
double * buf_all;
memory->create(buf_all, size_buf_all, "replicate:buf_all");
MPI_Allgatherv(buf,n,MPI_DOUBLE,buf_all,size_buf_rnk,disp_buf_rnk,
MPI_DOUBLE,world);
// bounding box of original unwrapped system
double _orig_lo[3], _orig_hi[3];
if (triclinic) {
_orig_lo[0] = domain->boxlo[0] +
_imagelo[0] * old_xprd + _imagelo[1] * old_xy + _imagelo[2] * old_xz;
_orig_lo[1] = domain->boxlo[1] +
_imagelo[1] * old_yprd + _imagelo[2] * old_yz;
_orig_lo[2] = domain->boxlo[2] + _imagelo[2] * old_zprd;
_orig_hi[0] = domain->boxlo[0] +
(_imagehi[0]+1) * old_xprd +
(_imagehi[1]+1) * old_xy + (_imagehi[2]+1) * old_xz;
_orig_hi[1] = domain->boxlo[1] +
(_imagehi[1]+1) * old_yprd + (_imagehi[2]+1) * old_yz;
_orig_hi[2] = domain->boxlo[2] + (_imagehi[2]+1) * old_zprd;
} else {
_orig_lo[0] = domain->boxlo[0] + _imagelo[0] * old_xprd;
_orig_lo[1] = domain->boxlo[1] + _imagelo[1] * old_yprd;
_orig_lo[2] = domain->boxlo[2] + _imagelo[2] * old_zprd;
_orig_hi[0] = domain->boxlo[0] + (_imagehi[0]+1) * old_xprd;
_orig_hi[1] = domain->boxlo[1] + (_imagehi[1]+1) * old_yprd;
_orig_hi[2] = domain->boxlo[2] + (_imagehi[2]+1) * old_zprd;
}
double _lo[3], _hi[3];
int num_replicas_added = 0;
// store x and tag for the whole system (before replication)
if (bond_flag) {
memory->create(old_x,old->natoms,3,"replicate:old_x");
memory->create(old_tag,old->natoms,"replicate:old_tag");
i = m = 0;
while (m < size_buf_all) {
old_x[i][0] = buf_all[m+1];
old_x[i][1] = buf_all[m+2];
old_x[i][2] = buf_all[m+3];
old_tag[i] = (tagint) ubuf(buf_all[m+4]).i;
old_map.insert({old_tag[i],i});
m += static_cast<int> (buf_all[m]);
i++;
}
}
for (ix = 0; ix < nx; ix++) {
for (iy = 0; iy < ny; iy++) {
for (iz = 0; iz < nz; iz++) {
thisrep[0] = ix;
thisrep[1] = iy;
thisrep[2] = iz;
// domain->remap() overwrites coordinates, so always recompute here
if (triclinic) {
_lo[0] = _orig_lo[0] + ix * old_xprd + iy * old_xy + iz * old_xz;
_hi[0] = _orig_hi[0] + ix * old_xprd + iy * old_xy + iz * old_xz;
_lo[1] = _orig_lo[1] + iy * old_yprd + iz * old_yz;
_hi[1] = _orig_hi[1] + iy * old_yprd + iz * old_yz;
_lo[2] = _orig_lo[2] + iz * old_zprd;
_hi[2] = _orig_hi[2] + iz * old_zprd;
} else {
_lo[0] = _orig_lo[0] + ix * old_xprd;
_hi[0] = _orig_hi[0] + ix * old_xprd;
_lo[1] = _orig_lo[1] + iy * old_yprd;
_hi[1] = _orig_hi[1] + iy * old_yprd;
_lo[2] = _orig_lo[2] + iz * old_zprd;
_hi[2] = _orig_hi[2] + iz * old_zprd;
}
// test if bounding box of shifted replica overlaps sub-domain of proc
// if not, then skip testing atoms
int xoverlap = 1;
int yoverlap = 1;
int zoverlap = 1;
if (triclinic) {
double _llo[3];
domain->x2lamda(_lo,_llo);
double _lhi[3];
domain->x2lamda(_hi,_lhi);
if (_llo[0] > (subhi[0] - EPSILON)
|| _lhi[0] < (sublo[0] + EPSILON) ) xoverlap = 0;
if (_llo[1] > (subhi[1] - EPSILON)
|| _lhi[1] < (sublo[1] + EPSILON) ) yoverlap = 0;
if (_llo[2] > (subhi[2] - EPSILON)
|| _lhi[2] < (sublo[2] + EPSILON) ) zoverlap = 0;
} else {
if (_lo[0] > (subhi[0] - EPSILON)
|| _hi[0] < (sublo[0] + EPSILON) ) xoverlap = 0;
if (_lo[1] > (subhi[1] - EPSILON)
|| _hi[1] < (sublo[1] + EPSILON) ) yoverlap = 0;
if (_lo[2] > (subhi[2] - EPSILON)
|| _hi[2] < (sublo[2] + EPSILON) ) zoverlap = 0;
}
int overlap = 0;
if (xoverlap && yoverlap && zoverlap) overlap = 1;
// if no overlap, test if bounding box wrapped back into new system
if (!overlap) {
// wrap back into cell
imageint imagelo = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
domain->remap(&(_lo[0]), imagelo);
int xboxlo = (imagelo & IMGMASK) - IMGMAX;
int yboxlo = (imagelo >> IMGBITS & IMGMASK) - IMGMAX;
int zboxlo = (imagelo >> IMG2BITS) - IMGMAX;
imageint imagehi = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
domain->remap(&(_hi[0]), imagehi);
int xboxhi = (imagehi & IMGMASK) - IMGMAX;
int yboxhi = (imagehi >> IMGBITS & IMGMASK) - IMGMAX;
int zboxhi = (imagehi >> IMG2BITS) - IMGMAX;
if (triclinic) {
double _llo[3];
_llo[0] = _lo[0]; _llo[1] = _lo[1]; _llo[2] = _lo[2];
domain->x2lamda(_llo,_lo);
double _lhi[3];
_lhi[0] = _hi[0]; _lhi[1] = _hi[1]; _lhi[2] = _hi[2];
domain->x2lamda(_lhi,_hi);
}
// test all fragments for any overlap; ok to include false positives
int _xoverlap1 = 0;
int _xoverlap2 = 0;
if (!xoverlap) {
if (xboxlo < 0) {
_xoverlap1 = 1;
if (_lo[0] > (subhi[0] - EPSILON)) _xoverlap1 = 0;
}
if (xboxhi > 0) {
_xoverlap2 = 1;
if (_hi[0] < (sublo[0] + EPSILON)) _xoverlap2 = 0;
}
if (_xoverlap1 || _xoverlap2) xoverlap = 1;
}
int _yoverlap1 = 0;
int _yoverlap2 = 0;
if (!yoverlap) {
if (yboxlo < 0) {
_yoverlap1 = 1;
if (_lo[1] > (subhi[1] - EPSILON)) _yoverlap1 = 0;
}
if (yboxhi > 0) {
_yoverlap2 = 1;
if (_hi[1] < (sublo[1] + EPSILON)) _yoverlap2 = 0;
}
if (_yoverlap1 || _yoverlap2) yoverlap = 1;
}
int _zoverlap1 = 0;
int _zoverlap2 = 0;
if (!zoverlap) {
if (zboxlo < 0) {
_zoverlap1 = 1;
if (_lo[2] > (subhi[2] - EPSILON)) _zoverlap1 = 0;
}
if (zboxhi > 0) {
_zoverlap2 = 1;
if (_hi[2] < (sublo[2] + EPSILON)) _zoverlap2 = 0;
}
if (_zoverlap1 || _zoverlap2) zoverlap = 1;
}
// does either fragment overlap w/ sub-domain
if (xoverlap && yoverlap && zoverlap) overlap = 1;
}
// while loop over one proc's atom list
if (overlap) {
num_replicas_added++;
m = 0;
while (m < size_buf_all) {
image = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
if (triclinic == 0) {
x[0] = buf_all[m+1] + ix*old_xprd;
x[1] = buf_all[m+2] + iy*old_yprd;
x[2] = buf_all[m+3] + iz*old_zprd;
} else {
x[0] = buf_all[m+1] + ix*old_xprd + iy*old_xy + iz*old_xz;
x[1] = buf_all[m+2] + iy*old_yprd + iz*old_yz;
x[2] = buf_all[m+3] + iz*old_zprd;
}
domain->remap(x,image);
if (triclinic) {
domain->x2lamda(x,lamda);
coord = lamda;
} else coord = x;
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) {
m += avec->unpack_restart(&buf_all[m]);
i = atom->nlocal - 1;
if (tag_enable)
atom_offset = iz*ny*nx*maxtag + iy*nx*maxtag + ix*maxtag;
else atom_offset = 0;
mol_offset = iz*ny*nx*maxmol + iy*nx*maxmol + ix*maxmol;
atom->x[i][0] = x[0];
atom->x[i][1] = x[1];
atom->x[i][2] = x[2];
atom0tag = atom->tag[i];
atom->tag[i] += atom_offset;
atom->image[i] = image;
if (atom->molecular != Atom::ATOMIC) {
if (atom->molecule[i] > 0)
atom->molecule[i] += mol_offset;
if (atom->molecular == Atom::MOLECULAR) {
if (atom->avec->bonds_allow)
for (j = 0; j < atom->num_bond[i]; j++) {
if (bond_flag)
newtag(atom0tag,atom->bond_atom[i][j]);
else atom->bond_atom[i][j] += atom_offset;
}
if (atom->avec->angles_allow)
for (j = 0; j < atom->num_angle[i]; j++) {
if (bond_flag) {
newtag(atom0tag,atom->angle_atom1[i][j]);
newtag(atom0tag,atom->angle_atom2[i][j]);
newtag(atom0tag,atom->angle_atom3[i][j]);
} else {
atom->angle_atom1[i][j] += atom_offset;
atom->angle_atom2[i][j] += atom_offset;
atom->angle_atom3[i][j] += atom_offset;
}
}
if (atom->avec->dihedrals_allow)
for (j = 0; j < atom->num_dihedral[i]; j++) {
if (bond_flag) {
newtag(atom0tag,atom->dihedral_atom1[i][j]);
newtag(atom0tag,atom->dihedral_atom2[i][j]);
newtag(atom0tag,atom->dihedral_atom3[i][j]);
newtag(atom0tag,atom->dihedral_atom4[i][j]);
} else {
atom->dihedral_atom1[i][j] += atom_offset;
atom->dihedral_atom2[i][j] += atom_offset;
atom->dihedral_atom3[i][j] += atom_offset;
atom->dihedral_atom4[i][j] += atom_offset;
}
}
if (atom->avec->impropers_allow)
for (j = 0; j < atom->num_improper[i]; j++) {
if (bond_flag) {
newtag(atom0tag,atom->improper_atom1[i][j]);
newtag(atom0tag,atom->improper_atom2[i][j]);
newtag(atom0tag,atom->improper_atom3[i][j]);
newtag(atom0tag,atom->improper_atom4[i][j]);
} else {
atom->improper_atom1[i][j] += atom_offset;
atom->improper_atom2[i][j] += atom_offset;
atom->improper_atom3[i][j] += atom_offset;
atom->improper_atom4[i][j] += atom_offset;
}
}
}
}
} else m += static_cast<int> (buf_all[m]);
}
} // if (overlap)
}
}
}
memory->destroy(size_buf_rnk);
memory->destroy(disp_buf_rnk);
memory->destroy(buf_all);
if (bond_flag) {
memory->destroy(old_x);
memory->destroy(old_tag);
}
int sum = 0;
MPI_Reduce(&num_replicas_added, &sum, 1, MPI_INT, MPI_SUM, 0, world);
double avg = (double) sum / nprocs;
if (me == 0)
utils::logmesg(lmp," average # of replicas added to proc = {:.2f} out "
"of {} ({:.2f}%)\n",avg,nx*ny*nz,avg/(nx*ny*nz)*100.0);
} else {
for (int iproc = 0; iproc < nprocs; iproc++) {
if (me == iproc) {
n = 0;
for (i = 0; i < old->nlocal; i++) n += old_avec->pack_restart(i,&buf[n]);
}
MPI_Bcast(&n,1,MPI_INT,iproc,world);
MPI_Bcast(buf,n,MPI_DOUBLE,iproc,world);
for (ix = 0; ix < nx; ix++) {
for (iy = 0; iy < ny; iy++) {
for (iz = 0; iz < nz; iz++) {
// while loop over one proc's atom list
m = 0;
while (m < n) {
image = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
if (triclinic == 0) {
x[0] = buf[m+1] + ix*old_xprd;
x[1] = buf[m+2] + iy*old_yprd;
x[2] = buf[m+3] + iz*old_zprd;
} else {
x[0] = buf[m+1] + ix*old_xprd + iy*old_xy + iz*old_xz;
x[1] = buf[m+2] + iy*old_yprd + iz*old_yz;
x[2] = buf[m+3] + iz*old_zprd;
}
domain->remap(x,image);
if (triclinic) {
domain->x2lamda(x,lamda);
coord = lamda;
} else coord = x;
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) {
m += avec->unpack_restart(&buf[m]);
i = atom->nlocal - 1;
if (tag_enable)
atom_offset = iz*ny*nx*maxtag + iy*nx*maxtag + ix*maxtag;
else atom_offset = 0;
mol_offset = iz*ny*nx*maxmol + iy*nx*maxmol + ix*maxmol;
atom->x[i][0] = x[0];
atom->x[i][1] = x[1];
atom->x[i][2] = x[2];
atom->tag[i] += atom_offset;
atom->image[i] = image;
if (atom->molecular != Atom::ATOMIC) {
if (atom->molecule[i] > 0)
atom->molecule[i] += mol_offset;
if (atom->molecular == Atom::MOLECULAR) {
if (atom->avec->bonds_allow)
for (j = 0; j < atom->num_bond[i]; j++)
atom->bond_atom[i][j] += atom_offset;
if (atom->avec->angles_allow)
for (j = 0; j < atom->num_angle[i]; j++) {
atom->angle_atom1[i][j] += atom_offset;
atom->angle_atom2[i][j] += atom_offset;
atom->angle_atom3[i][j] += atom_offset;
}
if (atom->avec->dihedrals_allow)
for (j = 0; j < atom->num_dihedral[i]; j++) {
atom->dihedral_atom1[i][j] += atom_offset;
atom->dihedral_atom2[i][j] += atom_offset;
atom->dihedral_atom3[i][j] += atom_offset;
atom->dihedral_atom4[i][j] += atom_offset;
}
if (atom->avec->impropers_allow)
for (j = 0; j < atom->num_improper[i]; j++) {
atom->improper_atom1[i][j] += atom_offset;
atom->improper_atom2[i][j] += atom_offset;
atom->improper_atom3[i][j] += atom_offset;
atom->improper_atom4[i][j] += atom_offset;
}
}
}
} else m += static_cast<int> (buf[m]);
}
}
}
}
}
} // if (bbox_flag || bond_flag)
*/
// free communication buffer and old atom class // free communication buffer and old atom class
memory->destroy(buf); memory->destroy(buf);
@ -955,16 +484,16 @@ void Replicate::replicate_by_proc(int nx, int ny, int nz,
} }
MPI_Bcast(&n,1,MPI_INT,iproc,world); MPI_Bcast(&n,1,MPI_INT,iproc,world);
MPI_Bcast(buf,n,MPI_DOUBLE,iproc,world); MPI_Bcast(buf,n,MPI_DOUBLE,iproc,world);
for (ix = 0; ix < nx; ix++) { for (ix = 0; ix < nx; ix++) {
for (iy = 0; iy < ny; iy++) { for (iy = 0; iy < ny; iy++) {
for (iz = 0; iz < nz; iz++) { for (iz = 0; iz < nz; iz++) {
// while loop over one proc's atom list // while loop over one proc's atom list
// x = new replicated position, remapped into new simulation box // x = new replicated position, remapped into new simulation box
// if atom is within my new subdomain, unpack it into new atom class // if atom is within my new subdomain, unpack it into new atom class
// adjust tag, mol #, coord, topology info as needed // adjust tag, mol #, coord, topology info as needed
m = 0; m = 0;
while (m < n) { while (m < n) {
image = ((imageint) IMGMAX << IMG2BITS) | image = ((imageint) IMGMAX << IMG2BITS) |
@ -983,25 +512,25 @@ void Replicate::replicate_by_proc(int nx, int ny, int nz,
domain->x2lamda(x,lamda); domain->x2lamda(x,lamda);
coord = lamda; coord = lamda;
} else coord = x; } else coord = x;
if (coord[0] >= sublo[0] && coord[0] < subhi[0] && if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] && coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) { coord[2] >= sublo[2] && coord[2] < subhi[2]) {
m += avec->unpack_restart(&buf[m]); m += avec->unpack_restart(&buf[m]);
i = atom->nlocal - 1; i = atom->nlocal - 1;
if (tag_enable) atom_offset = iz*ny*nx*maxtag + iy*nx*maxtag + ix*maxtag; if (tag_enable) atom_offset = iz*ny*nx*maxtag + iy*nx*maxtag + ix*maxtag;
else atom_offset = 0; else atom_offset = 0;
mol_offset = iz*ny*nx*maxmol + iy*nx*maxmol + ix*maxmol; mol_offset = iz*ny*nx*maxmol + iy*nx*maxmol + ix*maxmol;
atom->x[i][0] = x[0]; atom->x[i][0] = x[0];
atom->x[i][1] = x[1]; atom->x[i][1] = x[1];
atom->x[i][2] = x[2]; atom->x[i][2] = x[2];
atom->tag[i] += atom_offset; atom->tag[i] += atom_offset;
atom->image[i] = image; atom->image[i] = image;
if (atom->molecular != Atom::ATOMIC) { if (atom->molecular != Atom::ATOMIC) {
if (atom->molecule[i] > 0) if (atom->molecule[i] > 0)
atom->molecule[i] += mol_offset; atom->molecule[i] += mol_offset;
@ -1072,14 +601,14 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
int * size_buf_rnk; int * size_buf_rnk;
memory->create(size_buf_rnk, nprocs, "replicate:size_buf_rnk"); memory->create(size_buf_rnk, nprocs, "replicate:size_buf_rnk");
MPI_Allgather(&n, 1, MPI_INT, size_buf_rnk, 1, MPI_INT, world); MPI_Allgather(&n, 1, MPI_INT, size_buf_rnk, 1, MPI_INT, world);
// size of buf_all // size of buf_all
int size_buf_all = 0; int size_buf_all = 0;
MPI_Allreduce(&n, &size_buf_all, 1, MPI_INT, MPI_SUM, world); MPI_Allreduce(&n, &size_buf_all, 1, MPI_INT, MPI_SUM, world);
if (me == 0) { if (me == 0) {
auto mesg = fmt::format(" bounding box image = ({} {} {}) " auto mesg = fmt::format(" bounding box image = ({} {} {}) "
"to ({} {} {})\n", "to ({} {} {})\n",
@ -1089,25 +618,25 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
(double)size_buf_all*sizeof(double)/1024/1024); (double)size_buf_all*sizeof(double)/1024/1024);
utils::logmesg(lmp,mesg); utils::logmesg(lmp,mesg);
} }
// rnk offsets // rnk offsets
int *disp_buf_rnk; int *disp_buf_rnk;
memory->create(disp_buf_rnk, nprocs, "replicate:disp_buf_rnk"); memory->create(disp_buf_rnk, nprocs, "replicate:disp_buf_rnk");
disp_buf_rnk[0] = 0; disp_buf_rnk[0] = 0;
for (i = 1; i < nprocs; i++) for (i = 1; i < nprocs; i++)
disp_buf_rnk[i] = disp_buf_rnk[i-1] + size_buf_rnk[i-1]; disp_buf_rnk[i] = disp_buf_rnk[i-1] + size_buf_rnk[i-1];
// allgather buf_all // allgather buf_all
double *buf_all; double *buf_all;
memory->create(buf_all, size_buf_all, "replicate:buf_all"); memory->create(buf_all, size_buf_all, "replicate:buf_all");
MPI_Allgatherv(buf,n,MPI_DOUBLE,buf_all,size_buf_rnk,disp_buf_rnk, MPI_Allgatherv(buf,n,MPI_DOUBLE,buf_all,size_buf_rnk,disp_buf_rnk,
MPI_DOUBLE,world); MPI_DOUBLE,world);
// bounding box of original unwrapped system // bounding box of original unwrapped system
double _orig_lo[3], _orig_hi[3]; double _orig_lo[3], _orig_hi[3];
if (triclinic) { if (triclinic) {
_orig_lo[0] = domain->boxlo[0] + _orig_lo[0] = domain->boxlo[0] +
@ -1115,7 +644,7 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
_orig_lo[1] = domain->boxlo[1] + _orig_lo[1] = domain->boxlo[1] +
_imagelo[1] * old_yprd + _imagelo[2] * old_yz; _imagelo[1] * old_yprd + _imagelo[2] * old_yz;
_orig_lo[2] = domain->boxlo[2] + _imagelo[2] * old_zprd; _orig_lo[2] = domain->boxlo[2] + _imagelo[2] * old_zprd;
_orig_hi[0] = domain->boxlo[0] + _orig_hi[0] = domain->boxlo[0] +
(_imagehi[0]+1) * old_xprd + (_imagehi[0]+1) * old_xprd +
(_imagehi[1]+1) * old_xy + (_imagehi[2]+1) * old_xz; (_imagehi[1]+1) * old_xy + (_imagehi[2]+1) * old_xz;
@ -1126,23 +655,23 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
_orig_lo[0] = domain->boxlo[0] + _imagelo[0] * old_xprd; _orig_lo[0] = domain->boxlo[0] + _imagelo[0] * old_xprd;
_orig_lo[1] = domain->boxlo[1] + _imagelo[1] * old_yprd; _orig_lo[1] = domain->boxlo[1] + _imagelo[1] * old_yprd;
_orig_lo[2] = domain->boxlo[2] + _imagelo[2] * old_zprd; _orig_lo[2] = domain->boxlo[2] + _imagelo[2] * old_zprd;
_orig_hi[0] = domain->boxlo[0] + (_imagehi[0]+1) * old_xprd; _orig_hi[0] = domain->boxlo[0] + (_imagehi[0]+1) * old_xprd;
_orig_hi[1] = domain->boxlo[1] + (_imagehi[1]+1) * old_yprd; _orig_hi[1] = domain->boxlo[1] + (_imagehi[1]+1) * old_yprd;
_orig_hi[2] = domain->boxlo[2] + (_imagehi[2]+1) * old_zprd; _orig_hi[2] = domain->boxlo[2] + (_imagehi[2]+1) * old_zprd;
} }
double _lo[3], _hi[3]; double _lo[3], _hi[3];
int num_replicas_added = 0; int num_replicas_added = 0;
// if bond/periodic option // if bond/periodic option
// store old_x and old_tag for the entire original system // store old_x and old_tag for the entire original system
if (bond_flag) { if (bond_flag) {
memory->create(old_x,old->natoms,3,"replicate:old_x"); memory->create(old_x,old->natoms,3,"replicate:old_x");
memory->create(old_tag,old->natoms,"replicate:old_tag"); memory->create(old_tag,old->natoms,"replicate:old_tag");
i = m = 0; i = m = 0;
while (m < size_buf_all) { while (m < size_buf_all) {
old_x[i][0] = buf_all[m+1]; old_x[i][0] = buf_all[m+1];
@ -1156,40 +685,40 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
} }
// replication loop // replication loop
for (ix = 0; ix < nx; ix++) { for (ix = 0; ix < nx; ix++) {
for (iy = 0; iy < ny; iy++) { for (iy = 0; iy < ny; iy++) {
for (iz = 0; iz < nz; iz++) { for (iz = 0; iz < nz; iz++) {
thisrep[0] = ix; thisrep[0] = ix;
thisrep[1] = iy; thisrep[1] = iy;
thisrep[2] = iz; thisrep[2] = iz;
// domain->remap() overwrites coordinates, so always recompute here // domain->remap() overwrites coordinates, so always recompute here
if (triclinic) { if (triclinic) {
_lo[0] = _orig_lo[0] + ix * old_xprd + iy * old_xy + iz * old_xz; _lo[0] = _orig_lo[0] + ix * old_xprd + iy * old_xy + iz * old_xz;
_hi[0] = _orig_hi[0] + ix * old_xprd + iy * old_xy + iz * old_xz; _hi[0] = _orig_hi[0] + ix * old_xprd + iy * old_xy + iz * old_xz;
_lo[1] = _orig_lo[1] + iy * old_yprd + iz * old_yz; _lo[1] = _orig_lo[1] + iy * old_yprd + iz * old_yz;
_hi[1] = _orig_hi[1] + iy * old_yprd + iz * old_yz; _hi[1] = _orig_hi[1] + iy * old_yprd + iz * old_yz;
_lo[2] = _orig_lo[2] + iz * old_zprd; _lo[2] = _orig_lo[2] + iz * old_zprd;
_hi[2] = _orig_hi[2] + iz * old_zprd; _hi[2] = _orig_hi[2] + iz * old_zprd;
} else { } else {
_lo[0] = _orig_lo[0] + ix * old_xprd; _lo[0] = _orig_lo[0] + ix * old_xprd;
_hi[0] = _orig_hi[0] + ix * old_xprd; _hi[0] = _orig_hi[0] + ix * old_xprd;
_lo[1] = _orig_lo[1] + iy * old_yprd; _lo[1] = _orig_lo[1] + iy * old_yprd;
_hi[1] = _orig_hi[1] + iy * old_yprd; _hi[1] = _orig_hi[1] + iy * old_yprd;
_lo[2] = _orig_lo[2] + iz * old_zprd; _lo[2] = _orig_lo[2] + iz * old_zprd;
_hi[2] = _orig_hi[2] + iz * old_zprd; _hi[2] = _orig_hi[2] + iz * old_zprd;
} }
// test if bounding box of shifted replica overlaps sub-domain of proc // test if bounding box of shifted replica overlaps sub-domain of proc
// if not, then can skip testing of any individual atoms // if not, then can skip testing of any individual atoms
int xoverlap = 1; int xoverlap = 1;
int yoverlap = 1; int yoverlap = 1;
int zoverlap = 1; int zoverlap = 1;
@ -1198,7 +727,7 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
domain->x2lamda(_lo,_llo); domain->x2lamda(_lo,_llo);
double _lhi[3]; double _lhi[3];
domain->x2lamda(_hi,_lhi); domain->x2lamda(_hi,_lhi);
if (_llo[0] > (subhi[0] - EPSILON) if (_llo[0] > (subhi[0] - EPSILON)
|| _lhi[0] < (sublo[0] + EPSILON) ) xoverlap = 0; || _lhi[0] < (sublo[0] + EPSILON) ) xoverlap = 0;
if (_llo[1] > (subhi[1] - EPSILON) if (_llo[1] > (subhi[1] - EPSILON)
@ -1213,42 +742,42 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
if (_lo[2] > (subhi[2] - EPSILON) if (_lo[2] > (subhi[2] - EPSILON)
|| _hi[2] < (sublo[2] + EPSILON) ) zoverlap = 0; || _hi[2] < (sublo[2] + EPSILON) ) zoverlap = 0;
} }
int overlap = 0; int overlap = 0;
if (xoverlap && yoverlap && zoverlap) overlap = 1; if (xoverlap && yoverlap && zoverlap) overlap = 1;
// if no overlap, test if bounding box wrapped back into new system // if no overlap, test if bounding box wrapped back into new system
if (!overlap) { if (!overlap) {
// wrap back into cell // wrap back into cell
imageint imagelo = ((imageint) IMGMAX << IMG2BITS) | imageint imagelo = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX; ((imageint) IMGMAX << IMGBITS) | IMGMAX;
domain->remap(&(_lo[0]), imagelo); domain->remap(&(_lo[0]), imagelo);
int xboxlo = (imagelo & IMGMASK) - IMGMAX; int xboxlo = (imagelo & IMGMASK) - IMGMAX;
int yboxlo = (imagelo >> IMGBITS & IMGMASK) - IMGMAX; int yboxlo = (imagelo >> IMGBITS & IMGMASK) - IMGMAX;
int zboxlo = (imagelo >> IMG2BITS) - IMGMAX; int zboxlo = (imagelo >> IMG2BITS) - IMGMAX;
imageint imagehi = ((imageint) IMGMAX << IMG2BITS) | imageint imagehi = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX; ((imageint) IMGMAX << IMGBITS) | IMGMAX;
domain->remap(&(_hi[0]), imagehi); domain->remap(&(_hi[0]), imagehi);
int xboxhi = (imagehi & IMGMASK) - IMGMAX; int xboxhi = (imagehi & IMGMASK) - IMGMAX;
int yboxhi = (imagehi >> IMGBITS & IMGMASK) - IMGMAX; int yboxhi = (imagehi >> IMGBITS & IMGMASK) - IMGMAX;
int zboxhi = (imagehi >> IMG2BITS) - IMGMAX; int zboxhi = (imagehi >> IMG2BITS) - IMGMAX;
if (triclinic) { if (triclinic) {
double _llo[3]; double _llo[3];
_llo[0] = _lo[0]; _llo[1] = _lo[1]; _llo[2] = _lo[2]; _llo[0] = _lo[0]; _llo[1] = _lo[1]; _llo[2] = _lo[2];
domain->x2lamda(_llo,_lo); domain->x2lamda(_llo,_lo);
double _lhi[3]; double _lhi[3];
_lhi[0] = _hi[0]; _lhi[1] = _hi[1]; _lhi[2] = _hi[2]; _lhi[0] = _hi[0]; _lhi[1] = _hi[1]; _lhi[2] = _hi[2];
domain->x2lamda(_lhi,_hi); domain->x2lamda(_lhi,_hi);
} }
// test all fragments for any overlap; ok to include false positives // test all fragments for any overlap; ok to include false positives
int _xoverlap1 = 0; int _xoverlap1 = 0;
int _xoverlap2 = 0; int _xoverlap2 = 0;
if (!xoverlap) { if (!xoverlap) {
@ -1256,15 +785,15 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
_xoverlap1 = 1; _xoverlap1 = 1;
if (_lo[0] > (subhi[0] - EPSILON)) _xoverlap1 = 0; if (_lo[0] > (subhi[0] - EPSILON)) _xoverlap1 = 0;
} }
if (xboxhi > 0) { if (xboxhi > 0) {
_xoverlap2 = 1; _xoverlap2 = 1;
if (_hi[0] < (sublo[0] + EPSILON)) _xoverlap2 = 0; if (_hi[0] < (sublo[0] + EPSILON)) _xoverlap2 = 0;
} }
if (_xoverlap1 || _xoverlap2) xoverlap = 1; if (_xoverlap1 || _xoverlap2) xoverlap = 1;
} }
int _yoverlap1 = 0; int _yoverlap1 = 0;
int _yoverlap2 = 0; int _yoverlap2 = 0;
if (!yoverlap) { if (!yoverlap) {
@ -1272,16 +801,16 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
_yoverlap1 = 1; _yoverlap1 = 1;
if (_lo[1] > (subhi[1] - EPSILON)) _yoverlap1 = 0; if (_lo[1] > (subhi[1] - EPSILON)) _yoverlap1 = 0;
} }
if (yboxhi > 0) { if (yboxhi > 0) {
_yoverlap2 = 1; _yoverlap2 = 1;
if (_hi[1] < (sublo[1] + EPSILON)) _yoverlap2 = 0; if (_hi[1] < (sublo[1] + EPSILON)) _yoverlap2 = 0;
} }
if (_yoverlap1 || _yoverlap2) yoverlap = 1; if (_yoverlap1 || _yoverlap2) yoverlap = 1;
} }
int _zoverlap1 = 0; int _zoverlap1 = 0;
int _zoverlap2 = 0; int _zoverlap2 = 0;
if (!zoverlap) { if (!zoverlap) {
@ -1289,25 +818,25 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
_zoverlap1 = 1; _zoverlap1 = 1;
if (_lo[2] > (subhi[2] - EPSILON)) _zoverlap1 = 0; if (_lo[2] > (subhi[2] - EPSILON)) _zoverlap1 = 0;
} }
if (zboxhi > 0) { if (zboxhi > 0) {
_zoverlap2 = 1; _zoverlap2 = 1;
if (_hi[2] < (sublo[2] + EPSILON)) _zoverlap2 = 0; if (_hi[2] < (sublo[2] + EPSILON)) _zoverlap2 = 0;
} }
if (_zoverlap1 || _zoverlap2) zoverlap = 1; if (_zoverlap1 || _zoverlap2) zoverlap = 1;
} }
// does either fragment overlap w/ sub-domain // does either fragment overlap w/ sub-domain
if (xoverlap && yoverlap && zoverlap) overlap = 1; if (xoverlap && yoverlap && zoverlap) overlap = 1;
} }
// while loop over one proc's atom list // while loop over one proc's atom list
if (overlap) { if (overlap) {
num_replicas_added++; num_replicas_added++;
m = 0; m = 0;
while (m < size_buf_all) { while (m < size_buf_all) {
image = ((imageint) IMGMAX << IMG2BITS) | image = ((imageint) IMGMAX << IMG2BITS) |
@ -1326,27 +855,27 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
domain->x2lamda(x,lamda); domain->x2lamda(x,lamda);
coord = lamda; coord = lamda;
} else coord = x; } else coord = x;
if (coord[0] >= sublo[0] && coord[0] < subhi[0] && if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] && coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) { coord[2] >= sublo[2] && coord[2] < subhi[2]) {
m += avec->unpack_restart(&buf_all[m]); m += avec->unpack_restart(&buf_all[m]);
i = atom->nlocal - 1; i = atom->nlocal - 1;
if (tag_enable) if (tag_enable)
atom_offset = iz*ny*nx*maxtag + iy*nx*maxtag + ix*maxtag; atom_offset = iz*ny*nx*maxtag + iy*nx*maxtag + ix*maxtag;
else atom_offset = 0; else atom_offset = 0;
mol_offset = iz*ny*nx*maxmol + iy*nx*maxmol + ix*maxmol; mol_offset = iz*ny*nx*maxmol + iy*nx*maxmol + ix*maxmol;
atom->x[i][0] = x[0]; atom->x[i][0] = x[0];
atom->x[i][1] = x[1]; atom->x[i][1] = x[1];
atom->x[i][2] = x[2]; atom->x[i][2] = x[2];
atom0tag = atom->tag[i]; atom0tag = atom->tag[i];
atom->tag[i] += atom_offset; atom->tag[i] += atom_offset;
atom->image[i] = image; atom->image[i] = image;
if (atom->molecular != Atom::ATOMIC) { if (atom->molecular != Atom::ATOMIC) {
if (atom->molecule[i] > 0) if (atom->molecule[i] > 0)
atom->molecule[i] += mol_offset; atom->molecule[i] += mol_offset;
@ -1405,7 +934,7 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
} }
} }
} }
memory->destroy(size_buf_rnk); memory->destroy(size_buf_rnk);
memory->destroy(disp_buf_rnk); memory->destroy(disp_buf_rnk);
memory->destroy(buf_all); memory->destroy(buf_all);
@ -1413,7 +942,7 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz,
memory->destroy(old_x); memory->destroy(old_x);
memory->destroy(old_tag); memory->destroy(old_tag);
} }
int sum = 0; int sum = 0;
MPI_Reduce(&num_replicas_added, &sum, 1, MPI_INT, MPI_SUM, 0, world); MPI_Reduce(&num_replicas_added, &sum, 1, MPI_INT, MPI_SUM, 0, world);
double avg = (double) sum / nprocs; double avg = (double) sum / nprocs;

4
src/replicate.h
View File

@ -33,7 +33,7 @@ class Replicate : public Command {
private: private:
int bbox_flag, bond_flag; int bbox_flag, bond_flag;
class Atom *old; class Atom *old;
double old_xprd, old_yprd, old_zprd; double old_xprd, old_yprd, old_zprd;
@ -51,7 +51,7 @@ private:
void replicate_by_proc(int, int, int, double *, double *, double *); void replicate_by_proc(int, int, int, double *, double *, double *);
void replicate_by_bbox(int, int, int, double *, double *, double *); void replicate_by_bbox(int, int, int, double *, double *, double *);
void newtag(tagint, tagint &); void newtag(tagint, tagint &);
}; };