diff --git a/doc/src/create_atoms.rst b/doc/src/create_atoms.rst index 7935c676ef..1b7bcecd13 100644 --- a/doc/src/create_atoms.rst +++ b/doc/src/create_atoms.rst @@ -468,12 +468,13 @@ 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 -non-overlapping criterion. +*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. 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:: diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index fd1d535792..2912701159 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -316,6 +316,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) { @@ -651,6 +652,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, @@ -714,7 +716,8 @@ void CreateAtoms::add_single() if (mode == ATOM) { atom->avec->create_atom(ntype, xone); } else { - add_molecule(xone); + get_xmol(xone); + add_molecule(); } } } @@ -732,7 +735,6 @@ void CreateAtoms::add_random() if (overlapflag) { double odist = overlap; - if (mode == MOLECULE) odist += onemol->molradius; odistsq = odist * odist; } @@ -813,25 +815,43 @@ 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 { + 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++) { + 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; @@ -858,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(); } } } @@ -1429,7 +1453,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; @@ -1437,7 +1462,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(); } } @@ -1452,7 +1478,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]; @@ -1475,10 +1501,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; @@ -1488,9 +1510,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 };