diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md index d1cb2d480978b2becfa2dbf94f9e509e959b2438..af539f7a1039e3f4c7a9a02637f57f299ae12dab 100644 --- a/.github/PULL_REQUEST_TEMPLATE.md +++ b/.github/PULL_REQUEST_TEMPLATE.md @@ -4,7 +4,11 @@ _Briefly describe the new feature(s), enhancement(s), or bugfix(es) included in ## Implementation Notes -_Provide any relevant details about how the changes are implemented, how correctness was verified, how other features, if any, in LAMMPS are affected_ +_Provide any relevant details about how the changes are implemented, how correctness was verified, how other features - if any - in LAMMPS are affected_ + +## Backward Compatibility + +_Please state whether any changes in the pull request break backward compatibility for inputs, and - if yes - explain what has been changed and why_ ## Post Submission Checklist diff --git a/doc/src/create_bonds.txt b/doc/src/create_bonds.txt index f5a4c1f505df8d473f1f18544c2e99674f4c3d9c..5a878521693c968f09bbbbefbbca130a91832fa6 100644 --- a/doc/src/create_bonds.txt +++ b/doc/src/create_bonds.txt @@ -10,53 +10,93 @@ create_bonds command :h3 [Syntax:] -create_bonds group-ID group2-ID btype rmin rmax :pre - -group-ID = ID of first group -group2-ID = ID of second group, bonds will be between atoms in the 2 groups -btype = bond type of created bonds -rmin = minimum distance between pair of atoms to bond together -rmax = minimum distance between pair of atoms to bond together :ul +create_bonds style args ... keyword value ... :pre + +style = {many} or {single/bond} or {single/angle} or {single/dihedral} :ule,l + {many} args = group-ID group2-ID btype rmin rmax + group-ID = ID of first group + group2-ID = ID of second group, bonds will be between atoms in the 2 groups + btype = bond type of created bonds + rmin = minimum distance between pair of atoms to bond together + rmax = minimum distance between pair of atoms to bond together + {single/bond} args = btype batom1 batom2 + btype = bond type of new bond + batom1,batom2 = atom IDs for two atoms in bond + {single/angle} args = atype aatom1 aatom2 aatom3 + atype = bond type of new angle + aatom1,aatom2,aatom3 = atom IDs for three atoms in angle + {single/dihedral} args = dtype datom1 datom2 datom3 datom4 + dtype = bond type of new dihedral + datom1,datom2,datom3,datom4 = atom IDs for four atoms in dihedral :pre +zero or more keyword/value pairs may be appended :l +keyword = {special} :l + {special} value = {yes} or {no} :pre +:ule [Examples:] -create_bonds all all 1 1.0 1.2 -create_bonds surf solvent 3 2.0 2.4 :pre +create_bonds many all all 1 1.0 1.2 +create_bonds many surf solvent 3 2.0 2.4 +create_bond single/bond 1 1 2 +create_bond single/angle 5 52 98 107 special no :pre [Description:] -Create bonds between pairs of atoms that meet specified distance -criteria. The bond interactions can then be computed during a -simulation by the bond potential defined by the -"bond_style"_bond_style.html and "bond_coeff"_bond_coeff.html -commands. This command is useful for adding bonds to a system, -e.g. between nearest neighbors in a lattice of atoms, without having -to enumerate all the bonds in the data file read by the -"read_data"_read_data.html command. - -Note that the flexibility of this command is limited. It can be used -several times to create different types of bond at different -distances. But it cannot typically create all the bonds that would -normally be defined in a complex system of molecules. Also note that -this command does not add any 3-body or 4-body interactions which, -depending on your model, may be induced by added bonds, -e.g. "angle"_angle_style.html, "dihedral"_dihedral_style.html, or -"improper"_improper_style.html interactions. - -All created bonds will be between pairs of atoms I,J where I is in one -of the two specified groups, and J is in the other. The two groups -can be the same, e.g. group "all". The created bonds will be of bond -type {btype}, where {btype} must be a value between 1 and the number -of bond types defined. This maximum value is set by the "bond types" -field in the header of the data file read by the -"read_data"_read_data.html command, or via the optional "bond/types" -argument of the "create_box"_create_box.html command. +Create bonds between pairs of atoms that meet a specified distance +criteria. Or create a single bond, angle, or dihedral between 2, 3, +or 4 specified atoms. + +The new bond (angle, dihedral) interactions will then be computed +during a simulation by the bond (angle, dihedral) potential defined by +the "bond_style"_bond_style.html, "bond_coeff"_bond_coeff.html, +"angle_style"_angle_style.html, "angle_coeff"_angle_coeff.html, +"dihedral_style"_dihedral_style.html, +"dihedral_coeff"_dihedral_coeff.html commands. + +The {many} style is useful for adding bonds to a system, e.g. between +nearest neighbors in a lattice of atoms, without having to enumerate +all the bonds in the data file read by the "read_data"_read_data.html +command. + +The {single} styles are useful for adding bonds, angles, dihedrals +to a system incrementally, then continuing a simulation. + +Note that this command does not auto-create any angle or dihedral +interactions when a bond is added. Nor does it auto-create any bonds +when an angle or dihedral is added. Or auto-create any angles when a +dihedral is added. Thus the flexibility of this command is limited. +It can be used several times to create different types of bond at +different distances. But it cannot typically auto-create all the +bonds or angles or dihedral that would normally be defined in a data +file for a complex system of molecules. + +NOTE: If the system has no bonds (angles, dihedrals) to begin with, or +if more bonds per atom are being added than currently exist, then you +must insure that the number of bond types and the maximum number of +bonds per atom are set to large enough values. And similarly for +angles and dihedrals. Otherwise an error may occur when too many +bonds (angles, dihedrals) are added to an atom. If the +"read_data"_read_data.html command is used to define the system, these +parameters can be set via the "bond types" and "extra bond per atom" +fields in the header section of the data file. If the +"create_box"_create_box.html command is used to define the system, +these 2 parameters can be set via its optional "bond/types" and +"extra/bond/per/atom" arguments. And similarly for angles and +dihedrals. See the doc pages for these 2 commands for details. + +:line + +The {many} style will create bonds between pairs of atoms I,J where I +is in one of the two specified groups, and J is in the other. The two +groups can be the same, e.g. group "all". The created bonds will be +of bond type {btype}, where {btype} must be a value between 1 and the +number of bond types defined. For a bond to be created, an I,J pair of atoms must be a distance D apart such that {rmin} <= D <= {rmax}. -The following settings must have been made in an input -script before this command is used: +The following settings must have been made in an input script before +this style is used: special_bonds weight for 1-2 interactions must be 0.0 a "pair_style"_pair_style.html must be defined @@ -69,8 +109,8 @@ cannot appear in the neighbor list, to avoid creation of duplicate bonds. The neighbor list for all atom type pairs must also extend to a distance that encompasses the {rmax} for new bonds to create. -An additional requirement is that your system must be ready to perform -a simulation. This means, for example, that all +An additional requirement for this style is that your system must be +ready to perform a simulation. This means, for example, that all "pair_style"_pair_style.html coefficients be set via the "pair_coeff"_pair_coeff.html command. A "bond_style"_bond_style.html command and all bond coefficients must also be set, even if no bonds @@ -83,17 +123,58 @@ executes, e.g. if you wish to use long-range Coulombic interactions via the "kspace_style"_kspace_style.html command for your subsequent simulation. -NOTE: If the system has no bonds to begin with, or if more bonds per -atom are being added than currently exist, then you must insure that -the number of bond types and the maximum number of bonds per atom are -set to large enough values. Otherwise an error may occur when too -many bonds are added to an atom. If the "read_data"_read_data.html -command is used to define the system, these 2 parameters can be set -via the "bond types" and "extra bond per atom" fields in the header -section of the data file. If the "create_box"_create_box.html command -is used to define the system, these 2 parameters can be set via its -optional "bond/types" and "extra/bond/per/atom" arguments. See the -doc pages for the 2 commands for details. +:line + +The {single/bond} style creates a single bond of type {btype} between +two atoms with IDs {batom1} and {batom2}. {Btype} must be a value +between 1 and the number of bond types defined. + +The {single/angle} style creates a single angle of type {atype} +between three atoms with IDs {aatom1}, {aatom2}, and {aatom3}. The +ordering of the atoms is the same as in the {Angles} section of a data +file read by the "read_data"_read_data command. I.e. the 3 atoms are +ordered linearly within the angle; the central atom is {aatom2}. +{Atype} must be a value between 1 and the number of angle types +defined. + +The {single/dihedral} style creates a single dihedral of type {btype} +between two atoms with IDs {batom1} and {batom2}. The ordering of the +atoms is the same as in the {Dihedrals} section of a data file read by +the "read_data"_read_data command. I.e. the 4 atoms are ordered +linearly within the dihedral. {Dtype} must be a value between 1 and +the number of dihedral types defined. + +:line + +The keyword {special} controls whether an internal list of special +bonds is created after one or more bonds, or a single angle or +dihedral is added to the system. + +The default value is {yes}. A value of {no} cannot be used +with the {many} style. + +This is an expensive operation since the bond topology for the system +must be walked to find all 1-2, 1-3, 1-4 interactions to store in an +internal list, which is used when pairwise interactions are weighted; +see the "special_bonds"_special_bonds.html command for details. + +Thus if you are adding a few bonds or a large list of angles all at +the same time, by using this command repeatedly, it is more efficient +to only trigger the internal list to be created once, after the last +bond (or angle, or dihedral) is added: + +create_bonds single/bond 5 52 98 special no +create_bonds single/bond 5 73 74 special no +... +create_bonds single/bond 5 17 386 special no +create_bonds single/bond 4 112 183 special yes :pre + +Note that you MUST insure the internal list is re-built after the last +bond (angle, dihedral) is added, before performing a simulation. +Otherwise pairwise interactions will not be properly excluded or +weighted. LAMMPS does NOT check that you have done this correctly. + +:line [Restrictions:] @@ -105,4 +186,6 @@ molecule template files via the "molecule"_molecule.html and "create_atoms"_create_atoms.html, "delete_bonds"_delete_bonds.html -[Default:] none +[Default:] + +The keyword default is special = yes. diff --git a/src/create_bonds.cpp b/src/create_bonds.cpp index 7b3ca714a638dc7a36cf83f229589cb9929adc66..c62c6df91c5f1ddfbb197aff8b7d12cc7edd215f 100644 --- a/src/create_bonds.cpp +++ b/src/create_bonds.cpp @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Mike Salerno (NRL) added single methods +------------------------------------------------------------------------- */ + #include <stdlib.h> #include <string.h> #include "create_bonds.h" @@ -27,6 +31,8 @@ using namespace LAMMPS_NS; +enum{MANY,SBOND,SANGLE,SDIHEDRAL}; + /* ---------------------------------------------------------------------- */ CreateBonds::CreateBonds(LAMMPS *lmp) : Pointers(lmp) {} @@ -42,25 +48,104 @@ void CreateBonds::command(int narg, char **arg) if (atom->molecular != 1) error->all(FLERR,"Cannot use create_bonds with non-molecular system"); - if (narg != 5) error->all(FLERR,"Illegal create_bonds command"); + if (narg < 4) error->all(FLERR,"Illegal create_bonds command"); // parse args - int igroup = group->find(arg[0]); - if (igroup == -1) error->all(FLERR,"Cannot find create_bonds group ID"); - int group1bit = group->bitmask[igroup]; - igroup = group->find(arg[1]); - if (igroup == -1) error->all(FLERR,"Cannot find create_bonds group ID"); - int group2bit = group->bitmask[igroup]; + int style; + + int iarg; + if (strcmp(arg[0],"many") == 0) { + style = MANY; + if (narg != 6) error->all(FLERR,"Illegal create_bonds command"); + igroup = group->find(arg[1]); + if (igroup == -1) error->all(FLERR,"Cannot find create_bonds group ID"); + group1bit = group->bitmask[igroup]; + igroup = group->find(arg[2]); + if (igroup == -1) error->all(FLERR,"Cannot find create_bonds group ID"); + group2bit = group->bitmask[igroup]; + btype = force->inumeric(FLERR,arg[3]); + rmin = force->numeric(FLERR,arg[4]); + rmax = force->numeric(FLERR,arg[5]); + if (rmin > rmax) error->all(FLERR,"Illegal create_bonds command"); + iarg = 6; + } else if (strcmp(arg[0],"single/bond") == 0) { + style = SBOND; + if (narg < 4) error->all(FLERR,"Illegal create_bonds command"); + btype = force->inumeric(FLERR,arg[1]); + batom1 = force->tnumeric(FLERR,arg[2]); + batom2 = force->tnumeric(FLERR,arg[3]); + iarg = 4; + } else if (strcmp(arg[0],"single/angle") == 0) { + style = SANGLE; + if (narg < 5) error->all(FLERR,"Illegal create_bonds command"); + atype = force->inumeric(FLERR,arg[1]); + aatom1 = force->tnumeric(FLERR,arg[2]); + aatom2 = force->tnumeric(FLERR,arg[3]); + aatom3 = force->tnumeric(FLERR,arg[4]); + iarg = 5; + } else if (strcmp(arg[0],"single/dihedral") == 0) { + style = SDIHEDRAL; + if (narg < 6) error->all(FLERR,"Illegal create_bonds command"); + dtype = force->inumeric(FLERR,arg[1]); + datom1 = force->tnumeric(FLERR,arg[2]); + datom2 = force->tnumeric(FLERR,arg[3]); + datom3 = force->tnumeric(FLERR,arg[4]); + datom4 = force->tnumeric(FLERR,arg[5]); + iarg = 6; + } else error->all(FLERR,"Illegal create_bonds command"); + + // optional args + + int specialflag = 1; + + while (iarg < narg) { + if (strcmp(arg[iarg],"special") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal create_bonds command"); + if (strcmp(arg[iarg+1],"yes") == 0) specialflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) specialflag = 0; + else error->all(FLERR,"Illegal create_bonds command"); + iarg += 2; + } else error->all(FLERR,"Illegal create_bonds command"); + } + + // error checks + + if (style == MANY) { + if (btype <= 0 || btype > atom->nbondtypes) + error->all(FLERR,"Invalid bond type in create_bonds command"); + if (specialflag == 0) + error->all(FLERR,"Cannot use special no with create_bonds many"); + } else if (style == SBOND) { + if (btype <= 0 || btype > atom->nbondtypes) + error->all(FLERR,"Invalid bond type in create_bonds command"); + } else if (style == SANGLE) { + if (atype <= 0 || atype > atom->nangletypes) + error->all(FLERR,"Invalid angle type in create_bonds command"); + } else if (style == SDIHEDRAL) { + if (dtype <= 0 || dtype > atom->ndihedraltypes) + error->all(FLERR,"Invalid dihedral type in create_bonds command"); + } + + // invoke creation method + + if (style == MANY) many(); + else if (style == SBOND) single_bond(); + else if (style == SANGLE) single_angle(); + else if (style == SDIHEDRAL) single_dihedral(); - int btype = force->inumeric(FLERR,arg[2]); - double rmin = force->numeric(FLERR,arg[3]); - double rmax = force->numeric(FLERR,arg[4]); + // trigger special list build - if (btype <= 0 || btype > atom->nbondtypes) - error->all(FLERR,"Invalid bond type in create_bonds command"); - if (rmin > rmax) error->all(FLERR,"Illegal create_bonds command"); + if (specialflag) { + Special special(lmp); + special.build(); + } +} +/* ---------------------------------------------------------------------- */ + +void CreateBonds::many() +{ double rminsq = rmin*rmin; double rmaxsq = rmax*rmax; @@ -221,9 +306,184 @@ void CreateBonds::command(int narg, char **arg) nadd_bonds,atom->nbonds); } } +} + +/* ---------------------------------------------------------------------- */ + +void CreateBonds::single_bond() +{ + int m; + + // check that 2 atoms exist + + int count = 0; + if (atom->map(batom1) >= 0) count++; + if (atom->map(batom2) >= 0) count++; - // re-trigger special list build + int allcount; + MPI_Allreduce(&count,&allcount,1,MPI_INT,MPI_SUM,world); + if (allcount != 2) + error->all(FLERR,"Create_bonds single/bond atoms do not exist"); - Special special(lmp); - special.build(); + // create bond once or 2x if newton_bond set + + int *num_bond = atom->num_bond; + int **bond_type = atom->bond_type; + tagint **bond_atom = atom->bond_atom; + + if ((m = atom->map(batom1)) >= 0) { + if (num_bond[m] == atom->bond_per_atom) + error->one(FLERR,"New bond exceeded bonds per atom in create_bonds"); + bond_type[m][num_bond[m]] = btype; + bond_atom[m][num_bond[m]] = batom2; + num_bond[m]++; + } + + if (force->newton_bond) return; + + if ((m = atom->map(batom2)) >= 0) { + if (num_bond[m] == atom->bond_per_atom) + error->one(FLERR,"New bond exceeded bonds per atom in create_bonds"); + bond_type[m][num_bond[m]] = btype; + bond_atom[m][num_bond[m]] = batom1; + num_bond[m]++; + } +} + +/* ---------------------------------------------------------------------- */ + +void CreateBonds::single_angle() +{ + int m; + + // check that 3 atoms exist + + int count = 0; + if (atom->map(aatom1) >= 0) count++; + if (atom->map(aatom2) >= 0) count++; + if (atom->map(aatom3) >= 0) count++; + + int allcount; + MPI_Allreduce(&count,&allcount,1,MPI_INT,MPI_SUM,world); + if (allcount != 3) + error->all(FLERR,"Create_bonds single/angle atoms do not exist"); + + // create angle once or 3x if newton_bond set + + int *num_angle = atom->num_angle; + int **angle_type = atom->angle_type; + tagint **angle_atom1 = atom->angle_atom1; + tagint **angle_atom2 = atom->angle_atom2; + tagint **angle_atom3 = atom->angle_atom3; + + if ((m = atom->map(aatom2)) >= 0) { + if (num_angle[m] == atom->angle_per_atom) + error->one(FLERR,"New angle exceeded angles per atom in create_bonds"); + angle_type[m][num_angle[m]] = atype; + angle_atom1[m][num_angle[m]] = aatom1; + angle_atom2[m][num_angle[m]] = aatom2; + angle_atom3[m][num_angle[m]] = aatom3; + num_angle[m]++; + } + + if (force->newton_bond) return; + + if ((m = atom->map(aatom1)) >= 0) { + if (num_angle[m] == atom->angle_per_atom) + error->one(FLERR,"New angle exceeded angles per atom in create_bonds"); + angle_type[m][num_angle[m]] = atype; + angle_atom1[m][num_angle[m]] = aatom1; + angle_atom2[m][num_angle[m]] = aatom2; + angle_atom3[m][num_angle[m]] = aatom3; + num_angle[m]++; + } + + if ((m = atom->map(aatom3)) >= 0) { + if (num_angle[m] == atom->angle_per_atom) + error->one(FLERR,"New angle exceeded angles per atom in create_bonds"); + angle_type[m][num_angle[m]] = atype; + angle_atom1[m][num_angle[m]] = aatom1; + angle_atom2[m][num_angle[m]] = aatom2; + angle_atom3[m][num_angle[m]] = aatom3; + num_angle[m]++; + } +} + +/* ---------------------------------------------------------------------- */ + +void CreateBonds::single_dihedral() +{ + int m; + + // check that 4 atoms exist + + int count = 0; + if (atom->map(datom1) >= 0) count++; + if (atom->map(datom2) >= 0) count++; + if (atom->map(datom3) >= 0) count++; + if (atom->map(datom4) >= 0) count++; + + int allcount; + MPI_Allreduce(&count,&allcount,1,MPI_INT,MPI_SUM,world); + if (allcount != 4) + error->all(FLERR,"Create_bonds single/dihedral atoms do not exist"); + + // create bond once or 4x if newton_bond set + + int *num_dihedral = atom->num_dihedral; + int **dihedral_type = atom->dihedral_type; + tagint **dihedral_atom1 = atom->dihedral_atom1; + tagint **dihedral_atom2 = atom->dihedral_atom2; + tagint **dihedral_atom3 = atom->dihedral_atom3; + tagint **dihedral_atom4 = atom->dihedral_atom4; + + if ((m = atom->map(datom2)) >= 0) { + if (num_dihedral[m] == atom->dihedral_per_atom) + error->one(FLERR, + "New dihedral exceeded dihedrals per atom in create_bonds"); + dihedral_type[m][num_dihedral[m]] = dtype; + dihedral_atom1[m][num_dihedral[m]] = datom1; + dihedral_atom2[m][num_dihedral[m]] = datom2; + dihedral_atom3[m][num_dihedral[m]] = datom3; + dihedral_atom4[m][num_dihedral[m]] = datom4; + num_dihedral[m]++; + } + + if (force->newton_bond) return; + + if ((m = atom->map(datom1)) >= 0) { + if (num_dihedral[m] == atom->dihedral_per_atom) + error->one(FLERR, + "New dihedral exceeded dihedrals per atom in create_bonds"); + dihedral_type[m][num_dihedral[m]] = dtype; + dihedral_atom1[m][num_dihedral[m]] = datom1; + dihedral_atom2[m][num_dihedral[m]] = datom2; + dihedral_atom3[m][num_dihedral[m]] = datom3; + dihedral_atom4[m][num_dihedral[m]] = datom4; + num_dihedral[m]++; + } + + if ((m = atom->map(datom3)) >= 0) { + if (num_dihedral[m] == atom->dihedral_per_atom) + error->one(FLERR, + "New dihedral exceeded dihedrals per atom in create_bonds"); + dihedral_type[m][num_dihedral[m]] = dtype; + dihedral_atom1[m][num_dihedral[m]] = datom1; + dihedral_atom2[m][num_dihedral[m]] = datom2; + dihedral_atom3[m][num_dihedral[m]] = datom3; + dihedral_atom4[m][num_dihedral[m]] = datom4; + num_dihedral[m]++; + } + + if ((m = atom->map(datom4)) >= 0) { + if (num_dihedral[m] == atom->dihedral_per_atom) + error->one(FLERR, + "New dihedral exceeded dihedrals per atom in create_bonds"); + dihedral_type[m][num_dihedral[m]] = dtype; + dihedral_atom1[m][num_dihedral[m]] = datom1; + dihedral_atom2[m][num_dihedral[m]] = datom2; + dihedral_atom3[m][num_dihedral[m]] = datom3; + dihedral_atom4[m][num_dihedral[m]] = datom4; + num_dihedral[m]++; + } } diff --git a/src/create_bonds.h b/src/create_bonds.h index 24b1596e37ad68b8d6d5115a2a5e124a922e5a7f..175c74f5bcb8ac5e40b80e82a699cb32110300fd 100644 --- a/src/create_bonds.h +++ b/src/create_bonds.h @@ -30,9 +30,15 @@ class CreateBonds : protected Pointers { void command(int, char **); private: - inline int sbmask(int j) const { - return j >> SBBITS & 3; - } + int igroup,group1bit,group2bit; + int btype,atype,dtype; + tagint batom1,batom2,aatom1,aatom2,aatom3,datom1,datom2,datom3,datom4; + double rmin,rmax; + + void many(); + void single_bond(); + void single_angle(); + void single_dihedral(); }; }