From 161fdec54057e5000bd8bfdf8493caa9045f673f Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 10 Dec 2021 12:52:11 -0700 Subject: [PATCH] add improper amoeba class --- src/AMOEBA/improper_amoeba.cpp | 249 +++++++++++++++++++++++++++++++++ src/AMOEBA/improper_amoeba.h | 59 ++++++++ tools/tinker2lmp.py | 121 +++++++++++++++- 3 files changed, 427 insertions(+), 2 deletions(-) create mode 100644 src/AMOEBA/improper_amoeba.cpp create mode 100644 src/AMOEBA/improper_amoeba.h diff --git a/src/AMOEBA/improper_amoeba.cpp b/src/AMOEBA/improper_amoeba.cpp new file mode 100644 index 0000000000..d32b23ae2a --- /dev/null +++ b/src/AMOEBA/improper_amoeba.cpp @@ -0,0 +1,249 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "improper_amoeba.h" + +#include +#include "atom.h" +#include "comm.h" +#include "neighbor.h" +#include "force.h" +#include "update.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define TOLERANCE 0.05 +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +ImproperAmoeba::ImproperAmoeba(LAMMPS *lmp) : Improper(lmp) +{ + writedata = 1; +} + +/* ---------------------------------------------------------------------- */ + +ImproperAmoeba::~ImproperAmoeba() +{ + if (allocated && !copymode) { + memory->destroy(setflag); + memory->destroy(k); + } +} + +/* ---------------------------------------------------------------------- */ + +void ImproperAmoeba::compute(int eflag, int vflag) +{ + int ia,ib,ic,id,n,type; + double xia,yia,zia,xib,yib,zib,xic,yic,zic,xid,yid,zid; + double xab,yab,zab,xcb,ycb,zcb,xdb,ydb,zdb,xad,yad,zad,xcd,ycd,zcd; + double rad2,rcd2,rdb2,dot,cc,ee; + double sine,angle; + double eimproper,f1[3],f2[3],f3[3],f4[3]; + + eimproper = 0.0; + ev_init(eflag,vflag); + + double **x = atom->x; + double **f = atom->f; + int **improperlist = neighbor->improperlist; + int nimproperlist = neighbor->nimproperlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nimproperlist; n++) { + + // in Tinker code, atom1 = D, atom2 = B, atom3 = A, atom4 = C + // for Alligner angle: + // atoms A,C,D form a plane, B is out-of-plane + // angle is between plane and the vector from D to B + + id = improperlist[n][0]; + ib = improperlist[n][1]; + ia = improperlist[n][2]; + ic = improperlist[n][3]; + type = improperlist[n][4]; + + // coordinates of the atoms at trigonal center + + xia = x[ia][0]; + yia = x[ia][1]; + zia = x[ia][2]; + xib = x[ib][0]; + yib = x[ib][1]; + zib = x[ib][2]; + xic = x[ic][0]; + yic = x[ic][1]; + zic = x[ic][2]; + xid = x[id][0]; + yid = x[id][1]; + zid = x[id][2]; + + // compute the out-of-plane bending angle + + xab = xia - xib; + yab = yia - yib; + zab = zia - zib; + xcb = xic - xib; + ycb = yic - yib; + zcb = zic - zib; + xdb = xid - xib; + ydb = yid - yib; + zdb = zid - zib; + xad = xia - xid; + yad = yia - yid; + zad = zia - zid; + xcd = xic - xid; + ycd = yic - yid; + zcd = zic - zid; + + // Allinger angle between A-C-D plane and D-B vector for D-B < AC + + rad2 = xad*xad + yad*yad + zad*zad; + rcd2 = xcd*xcd + ycd*ycd + zcd*zcd; + dot = xad*xcd + yad*ycd + zad*zcd; + cc = rad2*rcd2 - dot*dot; + + // find the out-of-plane angle bending energy + + ee = xdb*(yab*zcb-zab*ycb) + ydb*(zab*xcb-xab*zcb) + zdb*(xab*ycb-yab*xcb); + rdb2 = xdb*xdb + ydb*ydb + zdb*zdb; + if (rdb2 == 0.0 || cc == 0.0) continue; + + sine = fabs(ee) / sqrt(cc*rdb2); + sine = MIN(1.0,sine); + angle = asin(sine) * MY_PI/180.0; + + //dt = angle; + //dt2 = dt * dt; + //dt3 = dt2 * dt; + //dt4 = dt2 * dt2; + //e = opbunit * force * dt2 * (1.0d0+copb*dt+qopb*dt2+popb*dt3+sopb*dt4); + //deddt = opbunit * force * dt * radian * + // (2.0d0 + 3.0d0*copb*dt + 4.0d0*qopb*dt2 + 5.0d0*popb*dt3 + 6.0d0*sopb*dt4); + //dedcos = -deddt * sign(1.0d0,ee) / sqrt(cc*rdb2-ee*ee); + + /* + // apply force to each of 4 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += f1[0]; + f[i1][1] += f1[1]; + f[i1][2] += f1[2]; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] += f2[0]; + f[i2][1] += f2[1]; + f[i2][2] += f2[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (newton_bond || i4 < nlocal) { + f[i4][0] += f4[0]; + f[i4][1] += f4[1]; + f[i4][2] += f4[2]; + } + + if (evflag) + ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4, + vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); + */ + } +} + +/* ---------------------------------------------------------------------- */ + +void ImproperAmoeba::allocate() +{ + allocated = 1; + int n = atom->nimpropertypes; + + memory->create(k,n+1,"improper:k"); + + memory->create(setflag,n+1,"improper:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void ImproperAmoeba::coeff(int narg, char **arg) +{ + if (narg != 2) error->all(FLERR,"Incorrect args for improper coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + utils::bounds(FLERR,arg[0],1,atom->nimpropertypes,ilo,ihi,error); + + double k_one = utils::numeric(FLERR,arg[1],false,lmp); + + // convert chi from degrees to radians + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for improper coefficients"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void ImproperAmoeba::write_restart(FILE *fp) +{ + fwrite(&k[1],sizeof(double),atom->nimpropertypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void ImproperAmoeba::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR,&k[1],sizeof(double),atom->nimpropertypes,fp,nullptr,error); + } + MPI_Bcast(&k[1],atom->nimpropertypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void ImproperAmoeba::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nimpropertypes; i++) + fprintf(fp,"%d %g\n",i,k[i]); +} diff --git a/src/AMOEBA/improper_amoeba.h b/src/AMOEBA/improper_amoeba.h new file mode 100644 index 0000000000..439a6c39fb --- /dev/null +++ b/src/AMOEBA/improper_amoeba.h @@ -0,0 +1,59 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef IMPROPER_CLASS +// clang-format off +ImproperStyle(amoeba,ImproperAmoeba); +// clang-format on +#else + +#ifndef LMP_IMPROPER_AMOEBA_H +#define LMP_IMPROPER_AMOEBA_H + +#include "improper.h" + +namespace LAMMPS_NS { + +class ImproperAmoeba : public Improper { + public: + ImproperAmoeba(class LAMMPS *); + ~ImproperAmoeba(); + void compute(int, int); + void coeff(int, char **); + void write_restart(FILE *); + void read_restart(FILE *); + void write_data(FILE *); + + protected: + double *k; + + virtual void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif + +/* ERROR/WARNING messages: + +W: Improper problem: %d %ld %d %d %d %d + +Conformation of the 4 listed improper atoms is extreme; you may want +to check your simulation geometry. + +E: Incorrect args for improper coefficients + +Self-explanatory. Check the input script or data file. + +*/ diff --git a/tools/tinker2lmp.py b/tools/tinker2lmp.py index a78d0d263c..63a5b6d154 100644 --- a/tools/tinker2lmp.py +++ b/tools/tinker2lmp.py @@ -10,6 +10,8 @@ # -nopbc = non-periodic system (default) # -pbc xhi yhi zhi = periodic system from 0 to hi in each dimension (optional) +# Author: Steve Plimpton + import sys,os,math path = os.environ["LAMMPS_PYTHON_TOOLS"] sys.path.append(path) @@ -105,6 +107,7 @@ class PRMfile: self.angleparams = self.angles() self.bondangleparams = self.bondangles() self.torsionparams = self.torsions() + self.opbendparams = self.opbend() self.ntypes = len(self.masses) def force_field_definition(self): @@ -280,6 +283,8 @@ class PRMfile: iline += 1 return params + # Dihedral interactions + def torsions(self): params = [] iline = self.find_section("Torsional Parameters") @@ -318,6 +323,28 @@ class PRMfile: params.append(oneparams) iline += 1 return params + + # Improper or out-of-plane bend interactions + + def opbend(self): + params = [] + iline = self.find_section("Out-of-Plane Bend Parameters") + if iline < 0: return params + iline += 3 + while iline < self.nlines: + words = self.lines[iline].split() + if len(words): + if words[0].startswith("###########"): break + if words[0] == "opbend": + class1 = int(words[1]) + class2 = int(words[2]) + class3 = int(words[3]) + class4 = int(words[4]) + value1 = float(words[5]) + lmp1 = value1 + params.append((class1,class2,class3,class4,lmp1)) + iline += 1 + return params def find_section(self,txt): txt = "## %s ##" % txt @@ -522,8 +549,30 @@ for atom2 in id: dlist.append((atom1,atom2,atom3,atom4)) ddict[(atom1,atom2,atom3,atom4)] = 1 +# create olist = list of out-of-plane impropers +# generate topology by triple loop over bonds of center atom2 +# atom2 must have 3 or more bonds to be part of an improper +# avoid double counting by requiring atom3 < atom4 +# this is since in Tinker the final 2 atoms in the improper are interchangeable + +id = xyz.id +type = xyz.type +bonds = xyz.bonds + +olist = [] + +for atom2 in id: + if len(bonds[atom2-1]) < 3: continue + for atom1 in bonds[atom2-1]: + for atom3 in bonds[atom2-1]: + for atom4 in bonds[atom2-1]: + if atom1 == atom3: continue + if atom1 == atom4: continue + if atom3 >= atom4: continue + olist.append((atom1,atom2,atom3,atom4)) + # ---------------------------------------- -# create lists of bond/angle/dihedral types +# create lists of bond/angle/dihedral/improper types # ---------------------------------------- # generate btype = LAMMPS type of each bond @@ -697,7 +746,7 @@ for atom1,atom2,atom3 in alist: # flags[i] = which LAMMPS dihedral type (1-N) the Tinker FF file dihedral I is # 0 = none # convert prm.torsionparams to a dictionary for efficient searching -# key = (class1,class2) +# key = (class1,class2,class3,class4) # value = (M,params) where M is index into prm.torsionparams id = xyz.id @@ -734,6 +783,54 @@ for atom1,atom2,atom3,atom4 in dlist: flags[m] = len(dparams) dtype.append(flags[m]) +# generate otype = LAMMPS type of each out-of-plane improper +# generate oparams = LAMMPS params for each improper type +# flags[i] = which LAMMPS improper type (1-N) the Tinker FF file improper I is +# 0 = none +# convert prm.opbendparams to a dictionary for efficient searching +# key = (class1,class2) +# value = (M,params) where M is index into prm.opbendparams + +id = xyz.id +type = xyz.type +classes = prm.classes + +odict = {} +for m,params in enumerate(prm.opbendparams): + odict[(params[0],params[1])] = (m,params) + +flags = len(prm.opbendparams)*[0] +otype = [] +oparams = [] +olist_reduced = [] + +for atom1,atom2,atom3,atom4 in olist: + type1 = type[atom1-1] + type2 = type[atom2-1] + type3 = type[atom3-1] + type4 = type[atom4-1] + c1 = classes[type1-1] + c2 = classes[type2-1] + c3 = classes[type3-1] + c4 = classes[type4-1] + + # 4-tuple is only an improper if matches an entry in PRM file + # olist_reduced = list of just these 4-tuples + + if (c1,c2) in odict: + m,params = odict[(c1,c2)] + olist_reduced.append((atom1,atom2,atom3,atom4)) + + if not flags[m]: + oneparams = params[4:] + oparams.append(oneparams) + flags[m] = len(oparams) + otype.append(flags[m]) + +# replace original olist with reduced version + +olist = olist_reduced + # ---------------------------------------- # assign each atom to a Tinker group # NOTE: doing this inside LAMMPS now @@ -787,6 +884,7 @@ ttype = xyz.type nbonds = len(blist) nangles = len(alist) ndihedrals = len(dlist) +nimpropers = len(olist) # data file header values @@ -885,6 +983,23 @@ if ndihedrals: lines.append(line+'\n') d.sections["Dihedrals"] = lines +if nimpropers: + d.headers["impropers"] = len(olist) + d.headers["improper types"] = len(oparams) + + lines = [] + for i,one in enumerate(oparams): + strone = [str(single) for single in one] + line = "%d %s" % (i+1,' '.join(strone)) + lines.append(line+'\n') + d.sections["Improper Coeffs"] = lines + + lines = [] + for i,one in enumerate(olist): + line = "%d %d %d %d %d %d" % (i+1,otype[i],one[0],one[1],one[2],one[3]) + lines.append(line+'\n') + d.sections["Impropers"] = lines + d.write(datafile) # print stats to screen @@ -898,6 +1013,8 @@ print "Nmol =",nmol print "Nbonds =",len(blist) print "Nangles =",len(alist) print "Ndihedrals =",len(dlist) +print "Nimpropers =",len(olist) print "Nbondtypes =",len(bparams) print "Nangletypes =",len(aparams) print "Ndihedraltypes =",len(dparams) +print "Nimpropertypes =",len(oparams)