ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add improper amoeba class

Browse Source
This commit is contained in:
Steve Plimpton
2021-12-10 12:52:11 -07:00
parent 657fcfa30d
commit 161fdec540
3 changed files with 427 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

249
src/AMOEBA/improper_amoeba.cpp Normal file
View File

@ -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 <cmath>
#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]);
}

59
src/AMOEBA/improper_amoeba.h Normal file
View File

@ -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.
*/

121
tools/tinker2lmp.py
View File

@ -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")
@ -319,6 +324,28 @@ class PRMfile:
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
for iline,line in enumerate(self.lines):
@ -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)