From 0e07089de689663a06c9754621e17a306e67ba18 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Wed, 7 Jun 2023 22:04:28 -0400 Subject: [PATCH 1/4] create_atoms:overlap_keyword_w_atomic_molecules --- src/create_atoms.cpp | 103 ++++++++++++++++++++++++++++++++----------- src/create_atoms.h | 5 ++- 2 files changed, 81 insertions(+), 27 deletions(-) diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 32be85e647..146eb23915 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -308,6 +308,7 @@ void CreateAtoms::command(int narg, char **arg) onemol->compute_center(); ranmol = new RanMars(lmp, molseed + comm->me); + memory->create(xmol, onemol->natoms, 3, "create_atoms:xmol"); } if (style == MESH) { @@ -634,6 +635,7 @@ void CreateAtoms::command(int narg, char **arg) delete[] xstr; delete[] ystr; delete[] zstr; + if (mode == MOLECULE) memory->destroy(xmol); // for MOLECULE mode: // create special bond lists for molecular systems, @@ -697,7 +699,8 @@ void CreateAtoms::add_single() if (mode == ATOM) { atom->avec->create_atom(ntype, xone); } else { - add_molecule(xone); + get_xmol(xone); + add_molecule(); } } } @@ -708,15 +711,17 @@ void CreateAtoms::add_single() void CreateAtoms::add_random() { + int incr; double xlo, ylo, zlo, xhi, yhi, zhi, zmid; double delx, dely, delz, distsq, odistsq; double lamda[3], *coord; - double *boxlo, *boxhi; + double *boxlo, *boxhi, *xmolbuf; if (overlapflag) { double odist = overlap; - if (mode == MOLECULE) odist += onemol->molradius; odistsq = odist * odist; + if (mode == MOLECULE) + memory->create(xmolbuf, onemol->natoms*3, "create_atoms:xmolbuf"); } // random number generator, same for all procs @@ -799,25 +804,54 @@ void CreateAtoms::add_random() // check for overlap of new atom/mol with all other atoms // including prior insertions // minimum_image() needed to account for distances across PBC - // new molecule only checks its center pt against others - // odistsq is expanded for mode=MOLECULE to account for molecule size if (overlapflag) { double **x = atom->x; int nlocal = atom->nlocal; int reject = 0; - for (int i = 0; i < nlocal; i++) { - delx = xone[0] - x[i][0]; - dely = xone[1] - x[i][1]; - delz = xone[2] - x[i][2]; - domain->minimum_image(delx, dely, delz); - distsq = delx * delx + dely * dely + delz * delz; - if (distsq < odistsq) { - reject = 1; - break; + if (mode == ATOM) { + for (int i = 0; i < nlocal; i++) { + delx = xone[0] - x[i][0]; + dely = xone[1] - x[i][1]; + delz = xone[2] - x[i][2]; + domain->minimum_image(delx, dely, delz); + distsq = delx * delx + dely * dely + delz * delz; + if (distsq < odistsq) { + reject = 1; + break; + } + } + } else { + int incr; + if (comm->me == 0) { + get_xmol(xone); + incr = 0; + for (int i = 0; i < onemol->natoms; i++) + for (int j = 0; j < 3; j++) + xmolbuf[incr++] = xmol[i][j]; + } + MPI_Bcast(xmolbuf, onemol->natoms*3, MPI_DOUBLE, 0, world); + incr = 0; + for (int i = 0; i < onemol->natoms; i++) + for (int j = 0; j < 3; j++) + xmol[i][j] = xmolbuf[incr++]; + + for (int i = 0; i < nlocal; i++) { + for (int j = 0; j < onemol->natoms; j++) { + delx = xmol[j][0] - x[i][0]; + dely = xmol[j][1] - x[i][1]; + delz = xmol[j][2] - x[i][2]; + domain->minimum_image(delx, dely, delz); + distsq = delx * delx + dely * dely + delz * delz; + if (distsq < odistsq) { + reject = 1; + break; + } + } } } + int reject_any; MPI_Allreduce(&reject, &reject_any, 1, MPI_INT, MPI_MAX, world); if (reject_any) continue; @@ -844,7 +878,11 @@ void CreateAtoms::add_random() if (mode == ATOM) { atom->avec->create_atom(ntype, xone); } else { - add_molecule(xone); + + // atomic coordinates calculated above for overlap check + + if (!overlapflag) get_xmol(xone); + add_molecule(); } } } @@ -857,6 +895,7 @@ void CreateAtoms::add_random() // clean-up delete random; + if (overlapflag && mode == MOLECULE) memory->destroy(xmolbuf); } /* ---------------------------------------------------------------------- @@ -1346,7 +1385,8 @@ void CreateAtoms::loop_lattice(int action) if (mode == ATOM) { atom->avec->create_atom(basistype[m], x); } else { - add_molecule(x); + get_xmol(x); + add_molecule(); } } else if (action == COUNT) { if (nlatt == MAXSMALLINT) nlatt_overflow = 1; @@ -1354,7 +1394,8 @@ void CreateAtoms::loop_lattice(int action) if (mode == ATOM) { atom->avec->create_atom(basistype[m], x); } else { - add_molecule(x); + get_xmol(x); + add_molecule(); } } @@ -1369,7 +1410,7 @@ void CreateAtoms::loop_lattice(int action) add a molecule with its center at center ------------------------------------------------------------------------- */ -void CreateAtoms::add_molecule(double *center) +void CreateAtoms::get_xmol(double *center) { double r[3], rotmat[3][3]; @@ -1392,10 +1433,6 @@ void CreateAtoms::add_molecule(double *center) MathExtra::quat_to_mat(quatone, rotmat); - // create atoms in molecule with atom ID = 0 and mol ID = 0 - // IDs are reset in caller after all molecules created by all procs - // pass add_molecule_atom an offset of 0 since don't know - // max tag of atoms in previous molecules at this point // onemol->quat_external is used by atom->add_moleclue_atom() onemol->quat_external = quatone; @@ -1405,9 +1442,25 @@ void CreateAtoms::add_molecule(double *center) for (int m = 0; m < natoms; m++) { MathExtra::matvec(rotmat, onemol->dx[m], xnew); MathExtra::add3(xnew, center, xnew); - atom->avec->create_atom(ntype + onemol->type[m], xnew); - n = atom->nlocal - 1; - atom->add_molecule_atom(onemol, m, n, 0); + for (int i = 0; i < 3; i++) + xmol[m][i] = xnew[i]; + } +} + +/* ---------------------------------------------------------------------- + add a molecule with atom coordinates from xmol +------------------------------------------------------------------------- */ + +void CreateAtoms::add_molecule() +{ + // create atoms in molecule with atom ID = 0 and mol ID = 0 + // IDs are reset in caller after all molecules created by all procs + // pass add_molecule_atom an offset of 0 since don't know + // max tag of atoms in previous molecules at this point + + for (int m = 0; m < onemol->natoms; m++) { + atom->avec->create_atom(ntype + onemol->type[m], xmol[m]); + atom->add_molecule_atom(onemol, m, atom->nlocal - 1, 0); } } diff --git a/src/create_atoms.h b/src/create_atoms.h index ae6f1b9d33..f839e3f0df 100644 --- a/src/create_atoms.h +++ b/src/create_atoms.h @@ -40,7 +40,7 @@ class CreateAtoms : public Command { bigint nsubset; double subsetfrac; int *basistype; - double xone[3], quatone[4]; + double xone[3], quatone[4], **xmol; double radthresh, radscale, mesh_density; int varflag, vvar, xvar, yvar, zvar; @@ -71,7 +71,8 @@ class CreateAtoms : public Command { int add_quasirandom(const double[3][3], tagint); void add_lattice(); void loop_lattice(int); - void add_molecule(double *); + void get_xmol(double *); + void add_molecule(); int vartest(double *); // evaluate a variable with new atom position }; From 3e2d5098c029a40edc7968167f468bc3c0b9c288 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Wed, 7 Jun 2023 22:07:57 -0400 Subject: [PATCH 2/4] Update create_atoms.rst --- doc/src/create_atoms.rst | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/doc/src/create_atoms.rst b/doc/src/create_atoms.rst index 5d1e7c872c..8748df9def 100644 --- a/doc/src/create_atoms.rst +++ b/doc/src/create_atoms.rst @@ -401,11 +401,9 @@ to. The *overlap* keyword only applies to the *random* style. It prevents newly created particles from being created closer than the specified -*Doverlap* distance from any other particle. When the particles being -created are molecules, the radius of the molecule (from its geometric -center) is added to *Doverlap*. If particles have finite size (see -:doc:`atom_style sphere ` for example) *Doverlap* should -be specified large enough to include the particle size in the +*Doverlap* distance from any other particle. If particles have finite +size (see :doc:`atom_style sphere ` for example) *Doverlap* +should be specified large enough to include the particle size in the non-overlapping criterion. .. note:: From 05438d23579c86d3abe4b38dbc1ee045c6436824 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Wed, 8 May 2024 18:50:02 -0400 Subject: [PATCH 3/4] Update create_atoms.rst --- doc/src/create_atoms.rst | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/doc/src/create_atoms.rst b/doc/src/create_atoms.rst index 1e2ea38bbf..1b7bcecd13 100644 --- a/doc/src/create_atoms.rst +++ b/doc/src/create_atoms.rst @@ -471,7 +471,10 @@ newly created particles from being created closer than the specified *Doverlap* distance from any other particle. If particles have finite size (see :doc:`atom_style sphere ` for example) *Doverlap* should be specified large enough to include the particle size in the -non-overlapping criterion. +non-overlapping criterion. If molecules are being randomly inserted, then +an insertion is only accepted if each particle in the molecule meets the +overlap criterion with respect to other particles (not including particles +in the molecule itself). .. note:: From 46b4c090832412c77d93777976c8dee73bd5e468 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 10 May 2024 00:15:21 -0400 Subject: [PATCH 4/4] simplify xmol comm --- src/create_atoms.cpp | 20 +++----------------- 1 file changed, 3 insertions(+), 17 deletions(-) diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index a0da871e97..2912701159 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -731,13 +731,11 @@ void CreateAtoms::add_random() double xlo, ylo, zlo, xhi, yhi, zhi; double delx, dely, delz, distsq, odistsq; double lamda[3], *coord; - double *boxlo, *boxhi, *xmolbuf; + double *boxlo, *boxhi; if (overlapflag) { double odist = overlap; odistsq = odist * odist; - if (mode == MOLECULE) - memory->create(xmolbuf, onemol->natoms*3, "create_atoms:xmolbuf"); } // random number generator, same for all procs @@ -836,19 +834,8 @@ void CreateAtoms::add_random() } } } else { - int incr; - if (comm->me == 0) { - get_xmol(xone); - incr = 0; - for (int i = 0; i < onemol->natoms; i++) - for (int j = 0; j < 3; j++) - xmolbuf[incr++] = xmol[i][j]; - } - MPI_Bcast(xmolbuf, onemol->natoms*3, MPI_DOUBLE, 0, world); - incr = 0; - for (int i = 0; i < onemol->natoms; i++) - for (int j = 0; j < 3; j++) - xmol[i][j] = xmolbuf[incr++]; + if (comm->me == 0) get_xmol(xone); + MPI_Bcast(&xmol[0][0], onemol->natoms*3, MPI_DOUBLE, 0, world); for (int i = 0; i < nlocal; i++) { for (int j = 0; j < onemol->natoms; j++) { @@ -908,7 +895,6 @@ void CreateAtoms::add_random() // clean-up delete random; - if (overlapflag && mode == MOLECULE) memory->destroy(xmolbuf); } /* ----------------------------------------------------------------------