jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

add single options to create_bonds command

This commit is contained in:
Steve Plimpton 2017-06-30 11:30:43 -06:00
parent 3bf2c60276
commit 711afe5062
3 changed files with 418 additions and 69 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

179
doc/src/create_bonds.txt
View File

@ -10,53 +10,93 @@ create_bonds command :h3
[Syntax:] [Syntax:]
create_bonds group-ID group2-ID btype rmin rmax :pre create_bonds style args ... keyword value ... :pre
group-ID = ID of first group style = {many} or {single/bond} or {single/angle} or {single/dihedral} :ule,l
group2-ID = ID of second group, bonds will be between atoms in the 2 groups {many} args = group-ID group2-ID btype rmin rmax
btype = bond type of created bonds group-ID = ID of first group
rmin = minimum distance between pair of atoms to bond together group2-ID = ID of second group, bonds will be between atoms in the 2 groups
rmax = minimum distance between pair of atoms to bond together :ul 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:] [Examples:]
create_bonds all all 1 1.0 1.2 create_bonds many all all 1 1.0 1.2
create_bonds surf solvent 3 2.0 2.4 :pre 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:] [Description:]
Create bonds between pairs of atoms that meet specified distance Create bonds between pairs of atoms that meet a specified distance
criteria. The bond interactions can then be computed during a criteria. Or create a single bond, angle, or dihedral between 2, 3,
simulation by the bond potential defined by the or 4 specified atoms.
"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 The new bond (angle, diehdral) interactions will then be computed
several times to create different types of bond at different during a simulation by the bond (angle, dihedral) potential defined by
distances. But it cannot typically create all the bonds that would the "bond_style"_bond_style.html, "bond_coeff"_bond_coeff.html,
normally be defined in a complex system of molecules. Also note that "angle_style"_angle_style.html, "angle_coeff"_angle_coeff.html,
this command does not add any 3-body or 4-body interactions which, "dihedral_style"_dihedral_style.html,
depending on your model, may be induced by added bonds, "dihedral_coeff"_dihedral_coeff.html commands.
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 The {many} style is useful for adding bonds to a system, e.g. between
of the two specified groups, and J is in the other. The two groups nearest neighbors in a lattice of atoms, without having to enumerate
can be the same, e.g. group "all". The created bonds will be of bond all the bonds in the data file read by the "read_data"_read_data.html
type {btype}, where {btype} must be a value between 1 and the number command.
of bond types defined. This maximum value is set by the "bond types"
field in the header of the data file read by the The {single} styles are useful for adding bonds, angles, dihedrals
"read_data"_read_data.html command, or via the optional "bond/types" to a system incrementally, then continuing a simulation.
argument of the "create_box"_create_box.html command.
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 For a bond to be created, an I,J pair of atoms must be a distance D
apart such that {rmin} <= D <= {rmax}. apart such that {rmin} <= D <= {rmax}.
The following settings must have been made in an input The following settings must have been made in an input script before
script before this command is used: this style is used:
special_bonds weight for 1-2 interactions must be 0.0 special_bonds weight for 1-2 interactions must be 0.0
a "pair_style"_pair_style.html must be defined 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 bonds. The neighbor list for all atom type pairs must also extend to
a distance that encompasses the {rmax} for new bonds to create. a distance that encompasses the {rmax} for new bonds to create.
An additional requirement is that your system must be ready to perform An additional requirement for this style is that your system must be
a simulation. This means, for example, that all ready to perform a simulation. This means, for example, that all
"pair_style"_pair_style.html coefficients be set via the "pair_style"_pair_style.html coefficients be set via the
"pair_coeff"_pair_coeff.html command. A "bond_style"_bond_style.html "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 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 via the "kspace_style"_kspace_style.html command for your subsequent
simulation. simulation.
NOTE: If the system has no bonds to begin with, or if more bonds per :line
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 The {single/bond} style creates a signle bond of type {btype} between
set to large enough values. Otherwise an error may occur when too two atoms with IDs {batom1} and {batom2}. {Btype} must be a value
many bonds are added to an atom. If the "read_data"_read_data.html between 1 and the number of bond types defined.
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 The {single/angle} style creates a single angle of type {atype}
section of the data file. If the "create_box"_create_box.html command between three atoms with IDs {aatom1}, {aatom2}, and {aatom3}. The
is used to define the system, these 2 parameters can be set via its ordering of the atoms is the same as in the {Angles} section of a data
optional "bond/types" and "extra/bond/per/atom" arguments. See the file read by the "read_data"_read_data command. I.e. the 3 atoms are
doc pages for the 2 commands for details. 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 signle 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:] [Restrictions:]
@ -105,4 +186,6 @@ molecule template files via the "molecule"_molecule.html and
"create_atoms"_create_atoms.html, "delete_bonds"_delete_bonds.html "create_atoms"_create_atoms.html, "delete_bonds"_delete_bonds.html
[Default:] none [Default:]
The keyword default is special = yes.

288
src/create_bonds.cpp
View File

@ -11,6 +11,10 @@
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Mike Salerno (NRL) added single methods
------------------------------------------------------------------------- */
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include "create_bonds.h" #include "create_bonds.h"
@ -27,6 +31,8 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
enum{MANY,SBOND,SANGLE,SDIHEDRAL};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
CreateBonds::CreateBonds(LAMMPS *lmp) : Pointers(lmp) {} CreateBonds::CreateBonds(LAMMPS *lmp) : Pointers(lmp) {}
@ -42,25 +48,104 @@ void CreateBonds::command(int narg, char **arg)
if (atom->molecular != 1) if (atom->molecular != 1)
error->all(FLERR,"Cannot use create_bonds with non-molecular system"); 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 // parse args
int igroup = group->find(arg[0]); int style;
if (igroup == -1) error->all(FLERR,"Cannot find create_bonds group ID");
int group1bit = group->bitmask[igroup]; 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]); igroup = group->find(arg[1]);
if (igroup == -1) error->all(FLERR,"Cannot find create_bonds group ID"); if (igroup == -1) error->all(FLERR,"Cannot find create_bonds group ID");
int group2bit = group->bitmask[igroup]; 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");
int btype = force->inumeric(FLERR,arg[2]); // optional args
double rmin = force->numeric(FLERR,arg[3]);
double rmax = force->numeric(FLERR,arg[4]);
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) if (btype <= 0 || btype > atom->nbondtypes)
error->all(FLERR,"Invalid bond type in create_bonds command"); error->all(FLERR,"Invalid bond type in create_bonds command");
if (rmin > rmax) error->all(FLERR,"Illegal 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();
// trigger special list build
if (specialflag) {
Special special(lmp);
special.build();
}
}
/* ---------------------------------------------------------------------- */
void CreateBonds::many()
{
double rminsq = rmin*rmin; double rminsq = rmin*rmin;
double rmaxsq = rmax*rmax; double rmaxsq = rmax*rmax;
@ -221,9 +306,184 @@ void CreateBonds::command(int narg, char **arg)
nadd_bonds,atom->nbonds); nadd_bonds,atom->nbonds);
} }
} }
}
// re-trigger special list build
/* ---------------------------------------------------------------------- */
Special special(lmp);
special.build(); 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++;
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");
// 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]++;
}
} }

12
src/create_bonds.h
View File

@ -30,9 +30,15 @@ class CreateBonds : protected Pointers {
void command(int, char **); void command(int, char **);
private: private:
inline int sbmask(int j) const { int igroup,group1bit,group2bit;
return j >> SBBITS & 3; 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();
}; };
} }