From 3028b6f34c9db9e0345cc62237da01bec7f856e7 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 6 May 2024 19:16:06 -0600 Subject: [PATCH] clean up of docs and code --- doc/src/replicate.rst | 20 +- examples/replicate/README | 9 + src/replicate.cpp | 621 +++++--------------------------------- src/replicate.h | 4 +- 4 files changed, 96 insertions(+), 558 deletions(-) diff --git a/doc/src/replicate.rst b/doc/src/replicate.rst index 64d87e525a..eb62fbfa21 100644 --- a/doc/src/replicate.rst +++ b/doc/src/replicate.rst @@ -8,7 +8,7 @@ Syntax .. code-block:: LAMMPS - replicate nx ny nz *keyword* ... + replicate nx ny nz keyword ... nx,ny,nz = replication factors in each dimension @@ -17,8 +17,8 @@ nx,ny,nz = replication factors in each dimension .. parsed-literal:: - *bbox* = only check atoms in replicas that overlap with a processor's subdomain - *bond/periodic* = use a different algorithm that correctly replicates periodic bond loops + *bbox* = use a bounding-box algorithm which is faster for large proc counts + *bond/periodic* = use an algorithm that correctly replicates periodic bond loops Examples """""""" @@ -56,7 +56,7 @@ are created between pairs of new atoms as well as between old and new atoms. .. note:: - + The bond discussion which follows only refers to models with permanent covalent bonds typically defined in LAMMPS via a data 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 all image flags are zero (in each dimension). --- note: +.. note:: LAMMPS does not check for image flag consistency before performing 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 there is a periodic bond loop. See the next note. --- note: +.. note:: 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 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 (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, before it is replicated. --- note: +.. note:: The algorithm used by the *bond/periodic* keyword builds on the 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. diff --git a/examples/replicate/README b/examples/replicate/README index 1363158b37..e33739672c 100644 --- a/examples/replicate/README +++ b/examples/replicate/README @@ -7,6 +7,8 @@ cross the periodic boundary to close the loop. To run these scripts, LAMMPS should be built with the MOLECULE and CLASS2 packages. The latter is only needed for the CNT example. +-------- + These scripts are tiny examples which illustrate both kinds of 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 @@ -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.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 periodic bond loops in all 3 dimensions xyz. diff --git a/src/replicate.cpp b/src/replicate.cpp index 3b6704f734..e07c7d9a26 100644 --- a/src/replicate.cpp +++ b/src/replicate.cpp @@ -13,7 +13,7 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing authors: + Contributing authors: Chris Knight (ANL) for bbox 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); // optional keywords - + bbox_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); maxmol = maxmol_all; } - + // reset image flags to zero for bond/periodic option - + if (bond_flag) for (i=0; inlocal; ++i) atom->image[i] = ((imageint) IMGMAX << IMG2BITS) | @@ -178,7 +178,7 @@ void Replicate::command(int narg, char **arg) for (i = 0; i < atom->nlocal; i++) domain->unmap(atom->x[i],atom->image[i]); - + // communication buffer for all my atom's info // max_size = largest buffer needed by any proc // 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) { replicate_by_proc(nx,ny,nz,sublo,subhi,buf); } else { 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 (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 (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 (buf[m]); - } - } - } - } - } - } // if (bbox_flag || bond_flag) - - */ - - - - - // free communication buffer and old atom class 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(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 // x = new replicated position, remapped into new simulation box // if atom is within my new subdomain, unpack it into new atom class // adjust tag, mol #, coord, topology info as needed - + m = 0; while (m < n) { image = ((imageint) IMGMAX << IMG2BITS) | @@ -983,25 +512,25 @@ void Replicate::replicate_by_proc(int nx, int ny, int nz, 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; @@ -1072,14 +601,14 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz, 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", @@ -1089,25 +618,25 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz, (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] + @@ -1115,7 +644,7 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz, _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; @@ -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[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; - + // if bond/periodic option // store old_x and old_tag for the entire original system - + 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]; @@ -1156,40 +685,40 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz, } // replication loop - + 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 can skip testing of any individual atoms - + int xoverlap = 1; int yoverlap = 1; int zoverlap = 1; @@ -1198,7 +727,7 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz, 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) @@ -1213,42 +742,42 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz, 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) { @@ -1256,15 +785,15 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz, _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) { @@ -1272,16 +801,16 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz, _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) { @@ -1289,25 +818,25 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz, _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) | @@ -1326,27 +855,27 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz, 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; @@ -1405,7 +934,7 @@ void Replicate::replicate_by_bbox(int nx, int ny, int nz, } } } - + memory->destroy(size_buf_rnk); memory->destroy(disp_buf_rnk); 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_tag); } - + int sum = 0; MPI_Reduce(&num_replicas_added, &sum, 1, MPI_INT, MPI_SUM, 0, world); double avg = (double) sum / nprocs; diff --git a/src/replicate.h b/src/replicate.h index 3877d2d720..defb35d1c6 100644 --- a/src/replicate.h +++ b/src/replicate.h @@ -33,7 +33,7 @@ class Replicate : public Command { private: int bbox_flag, bond_flag; - + class Atom *old; 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_bbox(int, int, int, double *, double *, double *); - + void newtag(tagint, tagint &); };