diff --git a/doc/src/Errors_messages.rst b/doc/src/Errors_messages.rst index ebd189faaa..252a17b012 100644 --- a/doc/src/Errors_messages.rst +++ b/doc/src/Errors_messages.rst @@ -4724,6 +4724,12 @@ Doc page with :doc:`WARNING messages ` *Invalid Masses section in molecule file* Self-explanatory. +*Invalid molecule ID in molecule file* + Molecule ID must be a non-zero positive integer. + +*Invalid Molecules section in molecule file* + Self-explanatory. + *Invalid REAX atom type* There is a mis-match between LAMMPS atom types and the elements listed in the ReaxFF force field file. @@ -4790,6 +4796,9 @@ Doc page with :doc:`WARNING messages ` Atom IDs must be positive integers and within range of defined atoms. +*Invalid atom ID in Fragments section of molecule file* + Self-explanatory. + *Invalid atom ID in Impropers section of data file* Atom IDs must be positive integers and within range of defined atoms. diff --git a/doc/src/molecule.rst b/doc/src/molecule.rst index 32e6d846ae..930b6cc5ce 100644 --- a/doc/src/molecule.rst +++ b/doc/src/molecule.rst @@ -60,6 +60,7 @@ templates include: * :doc:`fix rigid/small ` * :doc:`fix shake ` * :doc:`fix gcmc ` +* :doc:`fix bond/react ` * :doc:`create_atoms ` * :doc:`atom_style template ` @@ -144,6 +145,7 @@ appear if the value(s) are different than the default. * Na *angles* = # of angles Na in molecule, default = 0 * Nd *dihedrals* = # of dihedrals Nd in molecule, default = 0 * Ni *impropers* = # of impropers Ni in molecule, default = 0 +* Nf *fragments* = # of fragments in molecule, default = 0 * Mtotal *mass* = total mass of molecule * Xc Yc Zc *com* = coordinates of center-of-mass of molecule * Ixx Iyy Izz Ixy Ixz Iyz *inertia* = 6 components of inertia tensor of molecule @@ -166,7 +168,7 @@ internally. These are the allowed section keywords for the body of the file. -* *Coords, Types, Charges, Diameters, Masses* = atom-property sections +* *Coords, Types, Molecules, Fragments, Charges, Diameters, Masses* = atom-property sections * *Bonds, Angles, Dihedrals, Impropers* = molecular topology sections * *Special Bond Counts, Special Bonds* = special neighbor info * *Shake Flags, Shake Atoms, Shake Bond Types* = SHAKE info @@ -223,6 +225,26 @@ listed in order from 1 to Nlines, but LAMMPS does not check for this. ---------- +*Molecules* section: + +* one line per atom +* line syntax: ID molecule-ID +* molecule-ID = molecule ID of atom + +---------- + +*Fragments* section: + +* one line per fragment +* line syntax: ID a b c d ... +* a,b,c,d,... = IDs of atoms in fragment + +The ID of a fragment can only contain alphanumeric characters and +underscores. The atom IDs should be values from 1 to Natoms, where +Natoms = # of atoms in the molecule. + +---------- + *Charges* section: * one line per atom diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index 2255f64eb2..25806a1ac8 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -632,7 +632,13 @@ void FixPour::pre_exchange() int n = atom->nlocal - 1; atom->tag[n] = maxtag_all + m+1; if (mode == MOLECULE) { - if (atom->molecule_flag) atom->molecule[n] = maxmol_all+1; + if (atom->molecule_flag) { + if (onemols[imol]->moleculeflag) { + atom->molecule[n] = maxmol_all + onemols[imol]->molecule[m]; + } else { + atom->molecule[n] = maxmol_all+1; + } + } if (atom->molecular == 2) { atom->molindex[n] = 0; atom->molatom[n] = m; @@ -666,7 +672,13 @@ void FixPour::pre_exchange() fixshake->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat); maxtag_all += natom; - if (mode == MOLECULE && atom->molecule_flag) maxmol_all++; + if (mode == MOLECULE && atom->molecule_flag) { + if (onemols[imol]->moleculeflag) { + maxmol_all += onemols[imol]->nmolecules; + } else { + maxmol_all++; + } + } } // warn if not successful with all insertions b/c too many attempts diff --git a/src/MISC/fix_deposit.cpp b/src/MISC/fix_deposit.cpp index 90b1063924..a15a359881 100644 --- a/src/MISC/fix_deposit.cpp +++ b/src/MISC/fix_deposit.cpp @@ -539,7 +539,13 @@ void FixDeposit::pre_exchange() n = atom->nlocal - 1; atom->tag[n] = maxtag_all + m+1; if (mode == MOLECULE) { - if (atom->molecule_flag) atom->molecule[n] = maxmol_all+1; + if (atom->molecule_flag) { + if (onemols[imol]->moleculeflag) { + atom->molecule[n] = maxmol_all + onemols[imol]->molecule[m]; + } else { + atom->molecule[n] = maxmol_all+1; + } + } if (atom->molecular == 2) { atom->molindex[n] = 0; atom->molatom[n] = m; @@ -603,7 +609,13 @@ void FixDeposit::pre_exchange() maxtag_all += natom; if (maxtag_all >= MAXTAGINT) error->all(FLERR,"New atom IDs exceed maximum allowed ID"); - if (mode == MOLECULE && atom->molecule_flag) maxmol_all++; + if (mode == MOLECULE && atom->molecule_flag) { + if (onemols[imol]->moleculeflag) { + maxmol_all += onemols[imol]->nmolecules; + } else { + maxmol_all++; + } + } if (atom->map_style) { atom->map_init(); atom->map_set(); diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 212a70f80a..5325b8b9b1 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -490,7 +490,13 @@ void CreateAtoms::command(int narg, char **arg) for (int i = 0; i < molcreate; i++) { if (tag) offset = tag[ilocal]-1; for (int m = 0; m < natoms; m++) { - if (molecule_flag) molecule[ilocal] = moloffset + i+1; + if (molecule_flag) { + if (onemol->moleculeflag) { + molecule[ilocal] = moloffset + onemol->molecule[m]; + } else { + molecule[ilocal] = moloffset + 1; + } + } if (molecular == 2) { atom->molindex[ilocal] = 0; atom->molatom[ilocal] = m; @@ -524,6 +530,13 @@ void CreateAtoms::command(int narg, char **arg) } ilocal++; } + if (molecule_flag) { + if (onemol->moleculeflag) { + moloffset += onemol->nmolecules; + } else { + moloffset++; + } + } } // perform irregular comm to migrate atoms to new owning procs diff --git a/src/molecule.cpp b/src/molecule.cpp index 38887ebb0c..a8bf469a4a 100644 --- a/src/molecule.cpp +++ b/src/molecule.cpp @@ -37,7 +37,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : - Pointers(lmp), id(NULL), x(NULL), type(NULL), q(NULL), radius(NULL), + Pointers(lmp), id(NULL), x(NULL), type(NULL), molecule(NULL), q(NULL), radius(NULL), rmass(NULL), num_bond(NULL), bond_type(NULL), bond_atom(NULL), num_angle(NULL), angle_type(NULL), angle_atom1(NULL), angle_atom2(NULL), angle_atom3(NULL), num_dihedral(NULL), dihedral_type(NULL), dihedral_atom1(NULL), @@ -46,7 +46,7 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : improper_atom3(NULL), improper_atom4(NULL), nspecial(NULL), special(NULL), shake_flag(NULL), shake_atom(NULL), shake_type(NULL), avec_body(NULL), ibodyparams(NULL), dbodyparams(NULL), dx(NULL), dxcom(NULL), dxbody(NULL), quat_external(NULL), - fp(NULL), count(NULL) + fp(NULL), count(NULL), fragmentmask(NULL) { me = comm->me; @@ -146,19 +146,19 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : if (me == 0) { if (screen) - fprintf(screen,"Read molecule %s:\n" + fprintf(screen,"Read molecule template %s:\n %d molecules\n" " %d atoms with max type %d\n %d bonds with max type %d\n" " %d angles with max type %d\n %d dihedrals with max type %d\n" " %d impropers with max type %d\n", - id,natoms,ntypes, + id,nmolecules,natoms,ntypes, nbonds,nbondtypes,nangles,nangletypes, ndihedrals,ndihedraltypes,nimpropers,nimpropertypes); if (logfile) - fprintf(logfile,"Read molecule %s:\n" + fprintf(logfile,"Read molecule template %s:\n %d molecules\n" " %d atoms with max type %d\n %d bonds with max type %d\n" " %d angles with max type %d\n %d dihedrals with max type %d\n" " %d impropers with max type %d\n", - id,natoms,ntypes, + id,nmolecules,natoms,ntypes, nbonds,nbondtypes,nangles,nangletypes, ndihedrals,ndihedraltypes,nimpropers,nimpropertypes); } @@ -447,6 +447,9 @@ void Molecule::read(int flag) } else if (strstr(line,"impropers")) { nmatch = sscanf(line,"%d",&nimpropers); nwant = 1; + } else if (strstr(line,"fragments")) { + nmatch = sscanf(line,"%d",&nfragments); + nwant = 1; } else if (strstr(line,"mass")) { massflag = 1; nmatch = sscanf(line,"%lg",&masstotal); @@ -519,6 +522,14 @@ void Molecule::read(int flag) typeflag = 1; if (flag) types(line); else skip_lines(natoms,line); + } else if (strcmp(keyword,"Molecules") == 0) { + moleculeflag = 1; + if (flag) molecules(line); + else skip_lines(natoms,line); + } else if (strcmp(keyword,"Fragments") == 0) { + fragmentflag = 1; + if (flag) fragments(line); + else skip_lines(nfragments,line); } else if (strcmp(keyword,"Charges") == 0) { qflag = 1; if (flag) charges(line); @@ -694,6 +705,58 @@ void Molecule::types(char *line) ntypes = MAX(ntypes,type[i]); } +/* ---------------------------------------------------------------------- + read molecules from file + set nmolecules = max of any molecule type +------------------------------------------------------------------------- */ + +void Molecule::molecules(char *line) +{ + int tmp; + for (int i = 0; i < natoms; i++) { + readline(line); + if (2 != sscanf(line,"%d %d",&tmp,&molecule[i])) + error->all(FLERR,"Invalid Molecules section in molecule file"); + // molecule[i] += moffset; // placeholder for possible molecule offset + } + + for (int i = 0; i < natoms; i++) + if (molecule[i] <= 0) + error->all(FLERR,"Invalid molecule ID in molecule file"); + + for (int i = 0; i < natoms; i++) + nmolecules = MAX(nmolecules,molecule[i]); +} + +/* ---------------------------------------------------------------------- + read fragments from file +------------------------------------------------------------------------- */ + +void Molecule::fragments(char *line) +{ + int n,m,atomID,nwords; + char **words = new char*[natoms+1]; + + for (int i = 0; i < nfragments; i++) { + readline(line); + nwords = parse(line,words,natoms+1); + if (nwords > natoms+1) + error->all(FLERR,"Invalid atom ID in Fragments section of molecule file"); + n = strlen(words[0]) + 1; + fragmentnames[i] = new char[n]; + strcpy(fragmentnames[i],words[0]); + + for (m = 1; m < nwords; m++) { + atomID = atoi(words[m]); + if (atomID <= 0 || atomID > natoms) + error->all(FLERR,"Invalid atom ID in Fragments section of molecule file"); + fragmentmask[i][atomID-1] = 1; + } + } + + delete [] words; +} + /* ---------------------------------------------------------------------- read charges from file ------------------------------------------------------------------------- */ @@ -1359,6 +1422,17 @@ void Molecule::body(int flag, int pflag, char *line) } } +/* ---------------------------------------------------------------------- + return fragment index if name matches existing fragment, -1 if no such fragment +------------------------------------------------------------------------- */ + +int Molecule::findfragment(const char *name) +{ + for (int i = 0; i < nfragments; i++) + if (fragmentnames[i] && strcmp(name,fragmentnames[i]) == 0) return i; + return -1; +} + /* ---------------------------------------------------------------------- error check molecule attributes and topology against system settings flag = 0, just check this molecule @@ -1432,13 +1506,14 @@ void Molecule::initialize() natoms = 0; nbonds = nangles = ndihedrals = nimpropers = 0; ntypes = 0; + nmolecules = 1; nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0; nibody = ndbody = 0; bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0; maxspecial = 0; - xflag = typeflag = qflag = radiusflag = rmassflag = 0; + xflag = typeflag = moleculeflag = fragmentflag = qflag = radiusflag = rmassflag = 0; bondflag = angleflag = dihedralflag = improperflag = 0; nspecialflag = specialflag = 0; shakeflag = shakeflagflag = shakeatomflag = shaketypeflag = 0; @@ -1493,6 +1568,11 @@ void Molecule::allocate() { if (xflag) memory->create(x,natoms,3,"molecule:x"); if (typeflag) memory->create(type,natoms,"molecule:type"); + if (moleculeflag) memory->create(molecule,natoms,"molecule:molecule"); + if (fragmentflag) fragmentnames = new char*[nfragments]; + if (fragmentflag) memory->create(fragmentmask,nfragments,natoms,"molecule:fragmentmask"); + for (int i = 0; i < nfragments; i++) + for (int j = 0; j < natoms; j++) fragmentmask[i][j] = 0; if (qflag) memory->create(q,natoms,"molecule:q"); if (radiusflag) memory->create(radius,natoms,"molecule:radius"); if (rmassflag) memory->create(rmass,natoms,"molecule:rmass"); @@ -1580,10 +1660,17 @@ void Molecule::deallocate() { memory->destroy(x); memory->destroy(type); + memory->destroy(molecule); memory->destroy(q); memory->destroy(radius); memory->destroy(rmass); + memory->destroy(fragmentmask); + if (fragmentflag) { + for (int i = 0; i < nfragments; i++) delete [] fragmentnames[i]; + delete [] fragmentnames; + } + memory->destroy(num_bond); memory->destroy(bond_type); memory->destroy(bond_atom); diff --git a/src/molecule.h b/src/molecule.h index 064e58aabd..757eb3415c 100644 --- a/src/molecule.h +++ b/src/molecule.h @@ -30,10 +30,14 @@ class Molecule : protected Pointers { int natoms; int nbonds,nangles,ndihedrals,nimpropers; - int ntypes; + int ntypes,nmolecules,nfragments; int nbondtypes,nangletypes,ndihedraltypes,nimpropertypes; int nibody,ndbody; + // fragment info + char **fragmentnames; + int **fragmentmask; // nfragments by natoms + // max bond,angle,etc per atom int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom; @@ -41,7 +45,7 @@ class Molecule : protected Pointers { // 1 if attribute defined in file, 0 if not - int xflag,typeflag,qflag,radiusflag,rmassflag; + int xflag,typeflag,moleculeflag,fragmentflag,qflag,radiusflag,rmassflag; int bondflag,angleflag,dihedralflag,improperflag; int nspecialflag,specialflag; int shakeflag,shakeflagflag,shakeatomflag,shaketypeflag; @@ -59,6 +63,7 @@ class Molecule : protected Pointers { double **x; // displacement of each atom from origin int *type; // type of each atom + tagint *molecule; // molecule of each atom double *q; // charge on each atom double *radius; // radius of each atom double *rmass; // mass of each atom @@ -118,6 +123,7 @@ class Molecule : protected Pointers { void compute_mass(); void compute_com(); void compute_inertia(); + int findfragment(const char *); void check_attributes(int); private: @@ -131,6 +137,8 @@ class Molecule : protected Pointers { void read(int); void coords(char *); void types(char *); + void molecules(char *); + void fragments(char *); void charges(char *); void diameters(char *); void masses(char *); @@ -363,6 +371,18 @@ E: Invalid improper type in impropers section of molecule file Self-explanatory. +E: Invalid molecule ID in molecule file + +Molecule ID must be a non-zero positive integer. + +E: Invalid Molecules section in molecule file + +Self-explanatory. + +E: Invalid atom ID in Fragments section of molecule file + +Self-explanatory. + E: Invalid Special Bond Counts section in molecule file Self-explanatory.