diff --git a/src/create_bonds.cpp b/src/create_bonds.cpp index 5b7c354595..437f5959db 100644 --- a/src/create_bonds.cpp +++ b/src/create_bonds.cpp @@ -12,7 +12,9 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing authors: Mike Salerno (NRL) added single methods + Contributing authors: + Mike Salerno (NRL) added single methods + Thomas Farmer (ISIS) added single/improper ------------------------------------------------------------------------- */ #include "create_bonds.h" @@ -31,7 +33,7 @@ using namespace LAMMPS_NS; -enum{MANY,SBOND,SANGLE,SDIHEDRAL}; +enum{MANY,SBOND,SANGLE,SDIHEDRAL,SIMPROPER}; /* ---------------------------------------------------------------------- */ @@ -100,6 +102,18 @@ void CreateBonds::command(int narg, char **arg) (datom2 == datom3) || (datom2 == datom4) || (datom3 == datom4)) error->all(FLERR,"Illegal create_bonds command"); iarg = 6; + } else if (strcmp(arg[0],"single/improper") == 0) { + style = SIMPROPER; + 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]); + if ((datom1 == datom2) || (datom1 == datom3) || (datom1 == datom4) || + (datom2 == datom3) || (datom2 == datom4) || (datom3 == datom4)) + error->all(FLERR,"Illegal create_bonds command"); + iarg = 6; } else error->all(FLERR,"Illegal create_bonds command"); // optional args @@ -132,6 +146,9 @@ void CreateBonds::command(int narg, char **arg) } else if (style == SDIHEDRAL) { if (dtype <= 0 || dtype > atom->ndihedraltypes) error->all(FLERR,"Invalid dihedral type in create_bonds command"); + } else if (style == SIMPROPER) { + if (dtype <= 0 || dtype > atom->nimpropertypes) + error->all(FLERR,"Invalid improper type in create_bonds command"); } // invoke creation method @@ -140,6 +157,7 @@ void CreateBonds::command(int narg, char **arg) else if (style == SBOND) single_bond(); else if (style == SANGLE) single_angle(); else if (style == SDIHEDRAL) single_dihedral(); + else if (style == SIMPROPER) single_improper(); // trigger special list build @@ -512,3 +530,89 @@ void CreateBonds::single_dihedral() num_dihedral[m]++; } } + +/* ---------------------------------------------------------------------- */ + +void CreateBonds::single_improper() +{ + int m; + + // check that 4 atoms exist + + const int nlocal = atom->nlocal; + const int idx1 = atom->map(datom1); + const int idx2 = atom->map(datom2); + const int idx3 = atom->map(datom3); + const int idx4 = atom->map(datom4); + + int count = 0; + if ((idx1 >= 0) && (idx1 < nlocal)) count++; + if ((idx2 >= 0) && (idx2 < nlocal)) count++; + if ((idx3 >= 0) && (idx3 < nlocal)) count++; + if ((idx4 >= 0) && (idx4 < nlocal)) count++; + + int allcount; + MPI_Allreduce(&count,&allcount,1,MPI_INT,MPI_SUM,world); + if (allcount != 4) + error->all(FLERR,"Create_bonds single/improper atoms do not exist"); + + // create bond once or 4x if newton_bond set + + int *num_improper = atom->num_improper; + int **improper_type = atom->improper_type; + tagint **improper_atom1 = atom->improper_atom1; + tagint **improper_atom2 = atom->improper_atom2; + tagint **improper_atom3 = atom->improper_atom3; + tagint **improper_atom4 = atom->improper_atom4; + + if ((m = idx2) >= 0) { + if (num_improper[m] == atom->improper_per_atom) + error->one(FLERR, + "New improper exceeded impropers per atom in create_bonds"); + improper_type[m][num_improper[m]] = dtype; + improper_atom1[m][num_improper[m]] = datom1; + improper_atom2[m][num_improper[m]] = datom2; + improper_atom3[m][num_improper[m]] = datom3; + improper_atom4[m][num_improper[m]] = datom4; + num_improper[m]++; + } + atom->nimpropers++; + + if (force->newton_bond) return; + + if ((m = idx1) >= 0) { + if (num_improper[m] == atom->improper_per_atom) + error->one(FLERR, + "New improper exceeded impropers per atom in create_bonds"); + improper_type[m][num_improper[m]] = dtype; + improper_atom1[m][num_improper[m]] = datom1; + improper_atom2[m][num_improper[m]] = datom2; + improper_atom3[m][num_improper[m]] = datom3; + improper_atom4[m][num_improper[m]] = datom4; + num_improper[m]++; + } + + if ((m = idx3) >= 0) { + if (num_improper[m] == atom->improper_per_atom) + error->one(FLERR, + "New improper exceeded impropers per atom in create_bonds"); + improper_type[m][num_improper[m]] = dtype; + improper_atom1[m][num_improper[m]] = datom1; + improper_atom2[m][num_improper[m]] = datom2; + improper_atom3[m][num_improper[m]] = datom3; + improper_atom4[m][num_improper[m]] = datom4; + num_improper[m]++; + } + + if ((m = idx4) >= 0) { + if (num_improper[m] == atom->improper_per_atom) + error->one(FLERR, + "New improper exceeded impropers per atom in create_bonds"); + improper_type[m][num_improper[m]] = dtype; + improper_atom1[m][num_improper[m]] = datom1; + improper_atom2[m][num_improper[m]] = datom2; + improper_atom3[m][num_improper[m]] = datom3; + improper_atom4[m][num_improper[m]] = datom4; + num_improper[m]++; + } +} diff --git a/src/create_bonds.h b/src/create_bonds.h index 0c71242ed9..eea99b0113 100644 --- a/src/create_bonds.h +++ b/src/create_bonds.h @@ -39,6 +39,7 @@ class CreateBonds : protected Pointers { void single_bond(); void single_angle(); void single_dihedral(); + void single_improper(); }; } @@ -87,6 +88,10 @@ E: Invalid dihedral type in create_bonds command UNDOCUMENTED +E: Invalid improper type in create_bonds command + +UNDOCUMENTED + E: Create_bonds requires a pair style be defined Self-explanatory. @@ -135,4 +140,12 @@ E: New dihedral exceeded dihedrals per atom in create_bonds UNDOCUMENTED +E: Create_bonds single/improper atoms do not exist + +UNDOCUMENTED + +E: New improper exceeded impropers per atom in create_bonds + +UNDOCUMENTED + */