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

import updated version of tinker2lmp from @sjplimp. fix whitespace.

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-07-29 19:11:25 -04:00
parent 3a7a941cd5
commit 0dc9b31620
1 changed files with 133 additions and 132 deletions
Download Patch File Download Diff File Expand all files Collapse all files

265
tools/tinker/tinker2lmp.py
View File

@ -28,18 +28,16 @@ DELTA = 0.001 # delta on LAMMPS shrink-wrap box size, in Angstroms
# print an error message and exit # print an error message and exit
def error(txt=""): def error(txt="""
if not txt: Syntax: tinker2lmp.py -switch args ...
print("Syntax: tinker2lmp.py -switch args ...") -xyz file
print(" -xyz file") -amoeba file
print(" -amoeba file") -hippo file
print(" -hippo file") -data file
print(" -data file") -bitorsion file
print(" -bitorsion file") -nopbc
print(" -nopbc") -pbc xhi yhi zhi"""):
print(" -pbc xhi yhi zhi") sys.exit("ERROR: " + txt)
else: print("ERROR:",txt)
#sys.exit()
# read and store values from a Tinker xyz file # read and store values from a Tinker xyz file
@ -55,7 +53,7 @@ class XYZfile(object):
y = [] y = []
z = [] z = []
bonds = [] bonds = []
for line in lines[1:natoms+1]: for line in lines[1:natoms+1]:
words = line.split() words = line.split()
id.append(int(words[0])) id.append(int(words[0]))
@ -89,19 +87,19 @@ class XYZfile(object):
xhalfsq = 0.25 * boxhi[0]*boxhi[0] xhalfsq = 0.25 * boxhi[0]*boxhi[0]
yhalfsq = 0.25 * boxhi[1]*boxhi[1] yhalfsq = 0.25 * boxhi[1]*boxhi[1]
zhalfsq = 0.25 * boxhi[2]*boxhi[2] zhalfsq = 0.25 * boxhi[2]*boxhi[2]
x = self.x x = self.x
y = self.y y = self.y
z = self.z z = self.z
bonds = self.bonds bonds = self.bonds
imageflags = [] imageflags = []
for i in range(self.natoms): for i in range(self.natoms):
xi = float(x[i]) xi = float(x[i])
yi = float(y[i]) yi = float(y[i])
zi = float(z[i]) zi = float(z[i])
oneflags = [] oneflags = []
for j in bonds[i]: for j in bonds[i]:
xj = float(x[j-1]) xj = float(x[j-1])
yj = float(y[j-1]) yj = float(y[j-1])
@ -120,7 +118,7 @@ class XYZfile(object):
elif zi < zj: zimage = -1 elif zi < zj: zimage = -1
elif zi > zj: zimage = 1 elif zi > zj: zimage = 1
oneflags.append((ximage,yimage,zimage)) oneflags.append((ximage,yimage,zimage))
imageflags.append(oneflags) imageflags.append(oneflags)
self.imageflags = imageflags self.imageflags = imageflags
@ -128,7 +126,7 @@ class XYZfile(object):
# replicate system by Nx and Ny and Nz in each dim # replicate system by Nx and Ny and Nz in each dim
# rebuild data structs (except imageflags) in this class # rebuild data structs (except imageflags) in this class
# also alter global boxhi # also alter global boxhi
def replicate(self,nx,ny,nz): def replicate(self,nx,ny,nz):
natoms = self.natoms natoms = self.natoms
@ -140,11 +138,11 @@ class XYZfile(object):
z = self.z z = self.z
bonds = self.bonds bonds = self.bonds
imageflags = self.imageflags imageflags = self.imageflags
xbox = boxhi[0] xbox = boxhi[0]
ybox = boxhi[1] ybox = boxhi[1]
zbox = boxhi[2] zbox = boxhi[2]
idnew = [] idnew = []
labelnew = [] labelnew = []
typenew = [] typenew = []
@ -168,23 +166,23 @@ class XYZfile(object):
oneflags = imageflags[m] oneflags = imageflags[m]
onebonds = [] onebonds = []
for n,mbond in enumerate(bonds[m]): for n,mbond in enumerate(bonds[m]):
# ijk new = which replicated box the mbond atom is in # ijk new = which replicated box the mbond atom is in
if oneflags[n][0] == 0: inew = i if oneflags[n][0] == 0: inew = i
elif oneflags[n][0] == -1: inew = i - 1 elif oneflags[n][0] == -1: inew = i - 1
elif oneflags[n][0] == 1: inew = i + 1 elif oneflags[n][0] == 1: inew = i + 1
if inew < 0: inew = nx-1 if inew < 0: inew = nx-1
if inew == nx: inew = 0 if inew == nx: inew = 0
if oneflags[n][1] == 0: jnew = j if oneflags[n][1] == 0: jnew = j
elif oneflags[n][1] == -1: jnew = j - 1 elif oneflags[n][1] == -1: jnew = j - 1
elif oneflags[n][1] == 1: jnew = j + 1 elif oneflags[n][1] == 1: jnew = j + 1
if jnew < 0: jnew = ny-1 if jnew < 0: jnew = ny-1
if jnew == ny: jnew = 0 if jnew == ny: jnew = 0
if oneflags[n][2] == 0: knew = k if oneflags[n][2] == 0: knew = k
elif oneflags[n][2] == -1: knew = k - 1 elif oneflags[n][2] == -1: knew = k - 1
elif oneflags[n][2] == 1: knew = k + 1 elif oneflags[n][2] == 1: knew = k + 1
@ -193,12 +191,12 @@ class XYZfile(object):
# bondnew = replicated atom ID of bond partner # bondnew = replicated atom ID of bond partner
# based on replication box (inew,jnew,knew) that owns it # based on replication box (inew,jnew,knew) that owns it
bondnew = mbond + natoms*(knew*ny*nx + jnew*nx + inew) bondnew = mbond + natoms*(knew*ny*nx + jnew*nx + inew)
onebonds.append(bondnew) onebonds.append(bondnew)
bondsnew.append(onebonds) bondsnew.append(onebonds)
self.natoms = natoms * nx*ny*nz self.natoms = natoms * nx*ny*nz
self.id = idnew self.id = idnew
self.label = labelnew self.label = labelnew
@ -209,7 +207,7 @@ class XYZfile(object):
self.bonds = bondsnew self.bonds = bondsnew
# write out new xyzfile for replicated system # write out new xyzfile for replicated system
def output(self,outfile): def output(self,outfile):
fp = open(outfile,'w') fp = open(outfile,'w')
words = self.header.split() words = self.header.split()
@ -224,32 +222,32 @@ class XYZfile(object):
bonds = self.bonds bonds = self.bonds
# NOTE: worry about formatting of line # NOTE: worry about formatting of line
for i in range(self.natoms): for i in range(self.natoms):
print(i+1,label[i],x[i],y[i],z[i],type[i], end=' ', file=fp) print(i+1,label[i],x[i],y[i],z[i],type[i], end=' ', file=fp)
for j in bonds[i]: print(j, end=' ', file=fp) for j in bonds[i]: print(j, end=' ', file=fp)
print(file=fp) print(file=fp)
fp.close() fp.close()
# triplet of atoms in an angle = atom 1,2,3 # triplet of atoms in an angle = atom 1,2,3
# atom2 is center atom # atom2 is center atom
# nbonds = number of atoms which atom2 is bonded to # nbonds = number of atoms which atom2 is bonded to
# hcount = # of H atoms which atom2 is bonded to, excluding atom1 and atom3 # hcount = # of H atoms which atom2 is bonded to, excluding atom1 and atom3
def angle_hbond_count(self,atom1,atom2,atom3,lmptype,lmpmass): def angle_hbond_count(self,atom1,atom2,atom3,lmptype,lmpmass):
bondlist = self.bonds[atom2-1] bondlist = self.bonds[atom2-1]
nbonds = len(bondlist) nbonds = len(bondlist)
hcount = 0 hcount = 0
for bondID in bondlist: for bondID in bondlist:
if atom1 == bondID: continue if atom1 == bondID: continue
if atom3 == bondID: continue if atom3 == bondID: continue
massone = lmpmass[lmptype[bondID-1]-1] massone = lmpmass[lmptype[bondID-1]-1]
if int(massone+0.5) == 1: hcount += 1 if int(massone+0.5) == 1: hcount += 1
return nbonds,hcount return nbonds,hcount
# read and store select values from a Tinker force field PRM file # read and store select values from a Tinker force field PRM file
# ntypes = # of Tinker types # ntypes = # of Tinker types
# per-type values: class and mass # per-type values: class and mass
@ -276,7 +274,7 @@ class PRMfile(object):
self.ntypes = len(self.masses) self.ntypes = len(self.masses)
# find a section in the PRM file # find a section in the PRM file
def find_section(self,txt): def find_section(self,txt):
txt = "## %s ##" % txt txt = "## %s ##" % txt
for iline,line in enumerate(self.lines): for iline,line in enumerate(self.lines):
@ -284,7 +282,7 @@ class PRMfile(object):
return -1 return -1
# scalar params # scalar params
def force_field_definition(self): def force_field_definition(self):
iline = self.find_section("Force Field Definition") iline = self.find_section("Force Field Definition")
iline += 3 iline += 3
@ -301,7 +299,7 @@ class PRMfile(object):
iline += 1 iline += 1
# atom classes and masses # atom classes and masses
def peratom(self): def peratom(self):
classes = [] classes = []
masses = [] masses = []
@ -320,12 +318,12 @@ class PRMfile(object):
return classes,masses return classes,masses
# replicate per-atom data: classes and masses # replicate per-atom data: classes and masses
def replicate_peratom(self,nx,ny,nz): def replicate_peratom(self,nx,ny,nz):
classes = self.classes classes = self.classes
masses = self.masses masses = self.masses
natoms = len(self.masses) natoms = len(self.masses)
classes_new = [] classes_new = []
masses_new = [] masses_new = []
@ -338,9 +336,9 @@ class PRMfile(object):
self.classes = classes_new self.classes = classes_new
self.masses = masses_new self.masses = masses_new
# polarization groups # polarization groups
def polarize(self): def polarize(self):
polgroup = [] polgroup = []
iline = self.find_section("Dipole Polarizability Parameters") iline = self.find_section("Dipole Polarizability Parameters")
@ -362,7 +360,7 @@ class PRMfile(object):
return polgroup return polgroup
# convert PRMfile params to LAMMPS bond_style class2 params # convert PRMfile params to LAMMPS bond_style class2 params
def bonds(self): def bonds(self):
params = [] params = []
iline = self.find_section("Bond Stretching Parameters") iline = self.find_section("Bond Stretching Parameters")
@ -384,11 +382,11 @@ class PRMfile(object):
params.append((class1,class2,lmp1,lmp2,lmp3,lmp4)) params.append((class1,class2,lmp1,lmp2,lmp3,lmp4))
iline += 1 iline += 1
return params return params
# convert PRMfile params to LAMMPS angle_style class2/p6 params # convert PRMfile params to LAMMPS angle_style class2/p6 params
# line may have prefactor plus 1,2,3 angle0 params # line may have prefactor plus 1,2,3 angle0 params
# save prefactor/angle0 pairs as option 1,2,3 # save prefactor/angle0 pairs as option 1,2,3
def angles(self): def angles(self):
r2d = 180.0 / math.pi r2d = 180.0 / math.pi
ubflag = 0 ubflag = 0
@ -405,7 +403,7 @@ class PRMfile(object):
if words[0] == "angle" or words[0] == "anglep": if words[0] == "angle" or words[0] == "anglep":
if words[0] == "angle": pflag = 0 if words[0] == "angle": pflag = 0
if words[0] == "anglep": pflag = 1 if words[0] == "anglep": pflag = 1
class1 = int(words[1]) class1 = int(words[1])
class2 = int(words[2]) class2 = int(words[2])
class3 = int(words[3]) class3 = int(words[3])
@ -421,14 +419,14 @@ class PRMfile(object):
lmp4 = self.angle_quartic * value1 * r2d*r2d lmp4 = self.angle_quartic * value1 * r2d*r2d
lmp5 = self.angle_pentic * value1 * r2d*r2d*r2d lmp5 = self.angle_pentic * value1 * r2d*r2d*r2d
lmp6 = self.angle_sextic * value1 * r2d*r2d*r2d*r2d lmp6 = self.angle_sextic * value1 * r2d*r2d*r2d*r2d
option1 = (pflag,ubflag,lmp1,lmp2,lmp3,lmp4,lmp5,lmp6) option1 = (pflag,ubflag,lmp1,lmp2,lmp3,lmp4,lmp5,lmp6)
if len(words) >= 7: if len(words) >= 7:
value3 = float(words[6]) value3 = float(words[6])
lmp1 = value3 lmp1 = value3
option2 = (pflag,ubflag,lmp1,lmp2,lmp3,lmp4,lmp5,lmp6) option2 = (pflag,ubflag,lmp1,lmp2,lmp3,lmp4,lmp5,lmp6)
if len(words) == 8: if len(words) == 8:
value4 = float(words[7]) value4 = float(words[7])
lmp1 = value4 lmp1 = value4
@ -440,7 +438,7 @@ class PRMfile(object):
params.append((class1,class2,class3,[option1,option2])) params.append((class1,class2,class3,[option1,option2]))
else: else:
params.append((class1,class2,class3,[option1,option2,option3])) params.append((class1,class2,class3,[option1,option2,option3]))
iline += 1 iline += 1
return params return params
@ -448,7 +446,7 @@ class PRMfile(object):
# lmp3,lmp4 = equilibrium bond lengths for 2 bonds in angle # lmp3,lmp4 = equilibrium bond lengths for 2 bonds in angle
# need to find these values in self.bondparams # need to find these values in self.bondparams
# put them in a dictionary for efficient searching # put them in a dictionary for efficient searching
def bondangles(self): def bondangles(self):
params = [] params = []
iline = self.find_section("Stretch-Bend Parameters") iline = self.find_section("Stretch-Bend Parameters")
@ -471,15 +469,15 @@ class PRMfile(object):
value2 = float(words[5]) value2 = float(words[5])
lmp1 = value1 lmp1 = value1
lmp2 = value2 lmp2 = value2
if (class1,class2) in bdict: if (class1,class2) in bdict:
lmp3 = bdict[(class1,class2)] lmp3 = bdict[(class1,class2)]
elif (class2,class1) in bdict: elif (class2,class1) in bdict:
lmp3 = bdict[(class2,class1)] lmp3 = bdict[(class2,class1)]
else: else:
error("1st bond in BondAngle term not found: %d %d %d" % \
(class1,class2,class3))
# NOTE: just for debugging # NOTE: just for debugging
#error("1st bond in BondAngle term not found: %d %d %d" % \
# (class1,class2,class3))
lmp3 = 0.0 lmp3 = 0.0
if (class2,class3) in bdict: if (class2,class3) in bdict:
@ -487,17 +485,17 @@ class PRMfile(object):
elif (class3,class2) in bdict: elif (class3,class2) in bdict:
lmp4 = bdict[(class3,class2)] lmp4 = bdict[(class3,class2)]
else: else:
error("2nd bond in BondAngle term not found: %d %d %d" % \
(class1,class2,class3))
# NOTE: just for debugging # NOTE: just for debugging
#error("2nd bond in BondAngle term not found: %d %d %d" % \
# (class1,class2,class3))
lmp4 = 0.0 lmp4 = 0.0
params.append((class1,class2,class3,lmp1,lmp2,lmp3,lmp4)) params.append((class1,class2,class3,lmp1,lmp2,lmp3,lmp4))
iline += 1 iline += 1
return params return params
# dihedral interactions # dihedral interactions
def torsions(self): def torsions(self):
params = [] params = []
iline = self.find_section("Torsional Parameters") iline = self.find_section("Torsional Parameters")
@ -512,11 +510,11 @@ class PRMfile(object):
class2 = int(words[2]) class2 = int(words[2])
class3 = int(words[3]) class3 = int(words[3])
class4 = int(words[4]) class4 = int(words[4])
if len(words) <= 5: if len(words) <= 5:
error("torsion has no params: %d %d %d %d" % \ error("torsion has no params: %d %d %d %d" % \
(class1,class2,class3,class4)) (class1,class2,class3,class4))
if (len(words)-5) % 3: if (len(words)-5) % 3:
error("torsion does not have triplets of params: %d %d %d %d" % \ error("torsion does not have triplets of params: %d %d %d %d" % \
(class1,class2,class3,class4)) (class1,class2,class3,class4))
@ -531,13 +529,13 @@ class PRMfile(object):
lmp2 = value3 lmp2 = value3
lmp3 = value2 lmp3 = value2
oneparams += [lmp1,lmp2,lmp3] oneparams += [lmp1,lmp2,lmp3]
params.append(oneparams) params.append(oneparams)
iline += 1 iline += 1
return params return params
# improper or out-of-plane bend interactions # improper or out-of-plane bend interactions
def opbend(self): def opbend(self):
params = [] params = []
iline = self.find_section("Out-of-Plane Bend Parameters") iline = self.find_section("Out-of-Plane Bend Parameters")
@ -557,10 +555,10 @@ class PRMfile(object):
params.append((class1,class2,class3,class4,lmp1)) params.append((class1,class2,class3,class4,lmp1))
iline += 1 iline += 1
return params return params
# convert PRMfile params to LAMMPS angle_style amoeba UB params # convert PRMfile params to LAMMPS angle_style amoeba UB params
# coeffs for K2,K3 = 0.0 since Urey-Bradley is simple harmonic # coeffs for K2,K3 = 0.0 since Urey-Bradley is simple harmonic
def ureybonds(self): def ureybonds(self):
params = [] params = []
iline = self.find_section("Urey-Bradley Parameters") iline = self.find_section("Urey-Bradley Parameters")
@ -579,7 +577,7 @@ class PRMfile(object):
value2 = float(words[5]) value2 = float(words[5])
lmp1 = value1 lmp1 = value1
lmp2 = value2 lmp2 = value2
params.append((class1,class2,class3,lmp1,lmp2)) params.append((class1,class2,class3,lmp1,lmp2))
iline += 1 iline += 1
return params return params
@ -601,7 +599,7 @@ class PRMfile(object):
class2 = int(words[2]) class2 = int(words[2])
value1 = float(words[3]) value1 = float(words[3])
lmp1 = value1 lmp1 = value1
params.append((class1,class2,lmp1)) params.append((class1,class2,lmp1))
iline += 1 iline += 1
return params return params
@ -646,7 +644,7 @@ class PRMfile(object):
# ---------------------------------------- # ----------------------------------------
# main program # main program
# ---------------------------------------- # ----------------------------------------
args = sys.argv[1:] args = sys.argv[1:]
narg = len(args) narg = len(args)
@ -714,7 +712,7 @@ if not datafile: error("datafile not specified")
if not pbcflag and repflag: error("cannot replicate non-periodic system") if not pbcflag and repflag: error("cannot replicate non-periodic system")
if repflag and (nxrep <= 0 or nyrep <= 0 or nzrep <= 0): if repflag and (nxrep <= 0 or nyrep <= 0 or nzrep <= 0):
error("replication factors <= 0 not allowed") error("replication factors <= 0 not allowed")
# read Tinker xyz and prm files # read Tinker xyz and prm files
xyz = XYZfile(xyzfile) xyz = XYZfile(xyzfile)
@ -731,7 +729,7 @@ if repflag:
boxhi[0] *= nxrep boxhi[0] *= nxrep
boxhi[1] *= nyrep boxhi[1] *= nyrep
boxhi[2] *= nzrep boxhi[2] *= nzrep
# create LAMMPS box bounds based on pbcflag # create LAMMPS box bounds based on pbcflag
natoms = xyz.natoms natoms = xyz.natoms
@ -750,7 +748,7 @@ else:
zlo = min(zlo,float(z[i])) zlo = min(zlo,float(z[i]))
xhi = max(xhi,float(x[i])) xhi = max(xhi,float(x[i]))
yhi = max(yhi,float(y[i])) yhi = max(yhi,float(y[i]))
zhi = max(zhi,float(z[i])) zhi = max(zhi,float(z[i]))
boxlo = [xlo-DELTA,ylo-DELTA,zlo-DELTA] boxlo = [xlo-DELTA,ylo-DELTA,zlo-DELTA]
boxhi = [xhi+DELTA,yhi+DELTA,zhi+DELTA] boxhi = [xhi+DELTA,yhi+DELTA,zhi+DELTA]
@ -941,7 +939,7 @@ for atom1 in id:
if atom1 < atom2: if atom1 < atom2:
if len(bonds[atom1-1]) != 3: continue if len(bonds[atom1-1]) != 3: continue
if len(bonds[atom2-1]) != 3: continue if len(bonds[atom2-1]) != 3: continue
if bonds[atom1-1][0] == atom2: if bonds[atom1-1][0] == atom2:
atom3 = bonds[atom1-1][1] atom3 = bonds[atom1-1][1]
atom4 = bonds[atom1-1][2] atom4 = bonds[atom1-1][2]
@ -1006,7 +1004,7 @@ classes = prm.classes
bdict = {} bdict = {}
for m,params in enumerate(prm.bondparams): for m,params in enumerate(prm.bondparams):
bdict[(params[0],params[1])] = (m,params) bdict[(params[0],params[1])] = (m,params)
flags = len(prm.bondparams)*[0] flags = len(prm.bondparams)*[0]
btype = [] btype = []
bparams = [] bparams = []
@ -1016,17 +1014,17 @@ for atom1,atom2 in blist:
type2 = type[atom2-1] type2 = type[atom2-1]
c1 = classes[type1-1] c1 = classes[type1-1]
c2 = classes[type2-1] c2 = classes[type2-1]
if (c1,c2) in bdict: m,params = bdict[(c1,c2)] if (c1,c2) in bdict: m,params = bdict[(c1,c2)]
elif (c2,c1) in bdict: m,params = bdict[(c2,c1)] elif (c2,c1) in bdict: m,params = bdict[(c2,c1)]
else: error("bond not found: %d %d: %d %d" % (atom1,atom2,c1,c2)) else: error("bond not found: %d %d: %d %d" % (atom1,atom2,c1,c2))
if not flags[m]: if not flags[m]:
v1,v2,v3,v4 = params[2:] v1,v2,v3,v4 = params[2:]
bparams.append((v1,v2,v3,v4)) bparams.append((v1,v2,v3,v4))
flags[m] = len(bparams) flags[m] = len(bparams)
btype.append(flags[m]) btype.append(flags[m])
# generate atype = LAMMPS type of each angle # generate atype = LAMMPS type of each angle
# generate aparams = LAMMPS params for each angle type # generate aparams = LAMMPS params for each angle type
# flags[i] = which LAMMPS angle type (1-N) the Tinker FF file angle I is # flags[i] = which LAMMPS angle type (1-N) the Tinker FF file angle I is
@ -1047,7 +1045,7 @@ for m,params in enumerate(prm.angleparams):
adict[(params[0],params[1],params[2])] = (noptions,params) adict[(params[0],params[1],params[2])] = (noptions,params)
n = len(params[3]) n = len(params[3])
noptions += n noptions += n
flags = noptions*[0] flags = noptions*[0]
atype = [] atype = []
aparams = [] aparams = []
@ -1072,11 +1070,11 @@ for i,one in enumerate(alist):
# from Bond Stretching section of PRM file # from Bond Stretching section of PRM file
# since in general r1 != r2, the LAMMPS AngleAmoeba class requires # since in general r1 != r2, the LAMMPS AngleAmoeba class requires
# the 3 atoms in the angle be in the order that matches r1 and r2 # the 3 atoms in the angle be in the order that matches r1 and r2
if c1 != c3 and (c3,c2,c1) in adict: if c1 != c3 and (c3,c2,c1) in adict:
m,params = adict[(c3,c2,c1)] m,params = adict[(c3,c2,c1)]
alist[i] = (atom3,atom2,atom1) alist[i] = (atom3,atom2,atom1)
# params is a sequence of 1 or 2 or 3 options # params is a sequence of 1 or 2 or 3 options
# which = which of 1,2,3 options this atom triplet matches # which = which of 1,2,3 options this atom triplet matches
# for which = 2 or 3, increment m to index correct position in flags # for which = 2 or 3, increment m to index correct position in flags
@ -1093,41 +1091,43 @@ for i,one in enumerate(alist):
if len(params[3]) == 1: if len(params[3]) == 1:
which = 1 which = 1
elif len(params[3]) == 2: elif len(params[3]) == 2:
nbonds,hcount = xyz.angle_hbond_count(atom1,atom2,atom3,lmptype,lmpmass) nbonds,hcount = xyz.angle_hbond_count(atom1,atom2,atom3,lmptype,lmpmass)
if nbonds != 3: #if nbonds != 3:
print("Center angle atom has wrong bond count") #print("Center angle atom has wrong bond count")
print(" angle atom IDs:",atom1,atom2,atom3) #print(" angle atom IDs:",atom1,atom2,atom3)
print(" angle atom classes:",c1,c2,c3) #print(" angle atom classes:",c1,c2,c3)
print(" Tinker FF file param options:",len(params[3])) #print(" Tinker FF file param options:",len(params[3]))
print(" Nbonds and hydrogen count:",nbonds,hcount) #print(" Nbonds and hydrogen count:",nbonds,hcount)
#sys.exit() NOTE: allow this for now # NOTE: allow this for now
#sys.exit()
if hcount == 0: which = 1 if hcount == 0: which = 1
elif hcount == 1: elif hcount == 1:
which = 2 which = 2
m += 1 m += 1
print("3-bond angle") #print("3-bond angle")
print(" angle atom IDs:",atom1,atom2,atom3) #print(" angle atom IDs:",atom1,atom2,atom3)
print(" angle atom classes:",c1,c2,c3) #print(" angle atom classes:",c1,c2,c3)
print(" Tinker FF file param options:",len(params[3])) #print(" Tinker FF file param options:",len(params[3]))
print(" Nbonds and hydrogen count:",nbonds,hcount) #print(" Nbonds and hydrogen count:",nbonds,hcount)
print(" which:",which,m) #print(" which:",which,m)
elif len(params[3]) == 3: elif len(params[3]) == 3:
nbonds,hcount = xyz.angle_hbond_count(atom1,atom2,atom3,lmptype,lmpmass) nbonds,hcount = xyz.angle_hbond_count(atom1,atom2,atom3,lmptype,lmpmass)
if nbonds != 4: #if nbonds != 4:
print("Center angle atom has wrong bond count") #print("Center angle atom has wrong bond count")
print(" angle atom IDs:",atom1,atom2,atom3) #print(" angle atom IDs:",atom1,atom2,atom3)
print(" angle atom classes:",c1,c2,c3) #print(" angle atom classes:",c1,c2,c3)
print(" Tinker FF file param options:",len(params[3])) #print(" Tinker FF file param options:",len(params[3]))
print(" Nbonds and hydrogen count:",nbonds,hcount) #print(" Nbonds and hydrogen count:",nbonds,hcount)
#sys.exit() NOTE: allow this for now # NOTE: allow this for now
#sys.exit()
if hcount == 0: which = 1 if hcount == 0: which = 1
elif hcount == 1: elif hcount == 1:
which = 2 which = 2
@ -1138,7 +1138,7 @@ for i,one in enumerate(alist):
else: else:
error("angle not found: %d %d %d: %d %d %d" % (atom1,atom2,atom3,c1,c2,c3)) error("angle not found: %d %d %d: %d %d %d" % (atom1,atom2,atom3,c1,c2,c3))
if not flags[m]: if not flags[m]:
pflag,ubflag,v1,v2,v3,v4,v5,v6 = params[3][which-1] pflag,ubflag,v1,v2,v3,v4,v5,v6 = params[3][which-1]
aparams.append((pflag,ubflag,v1,v2,v3,v4,v5,v6)) aparams.append((pflag,ubflag,v1,v2,v3,v4,v5,v6))
@ -1157,7 +1157,7 @@ for v1,v2,v3,v4,v5,v6,v7 in prm.bondangleparams:
baparams = [] baparams = []
for itype in range(len(aparams)): for itype in range(len(aparams)):
iangle = atype.index(itype+1) iangle = atype.index(itype+1)
atom1,atom2,atom3 = alist[iangle] atom1,atom2,atom3 = alist[iangle]
type1 = type[atom1-1] type1 = type[atom1-1]
type2 = type[atom2-1] type2 = type[atom2-1]
@ -1171,9 +1171,10 @@ for itype in range(len(aparams)):
elif (c3,c2,c1) in badict: elif (c3,c2,c1) in badict:
n1,n2,r1,r2 = badict[(c3,c2,c1)] n1,n2,r1,r2 = badict[(c3,c2,c1)]
else: else:
print("Bond-stretch angle triplet not found: %d %d %d" % (c1,c2,c3)) # NOTE: just for debugging
#print("Bond-stretch angle triplet not found: %d %d %d" % (c1,c2,c3))
n1,n2,r1,r2 = 4*[0.0] n1,n2,r1,r2 = 4*[0.0]
baparams.append((n1,n2,r1,r2)) baparams.append((n1,n2,r1,r2))
# augment the aparams with Urey_Bradley terms from ureyparams # augment the aparams with Urey_Bradley terms from ureyparams
@ -1188,7 +1189,7 @@ for v1,v2,v3,v4,v5 in prm.ureyparams:
ubparams = [] ubparams = []
for itype in range(len(aparams)): for itype in range(len(aparams)):
iangle = atype.index(itype+1) iangle = atype.index(itype+1)
atom1,atom2,atom3 = alist[iangle] atom1,atom2,atom3 = alist[iangle]
type1 = type[atom1-1] type1 = type[atom1-1]
type2 = type[atom2-1] type2 = type[atom2-1]
@ -1198,7 +1199,7 @@ for itype in range(len(aparams)):
c3 = classes[type3-1] c3 = classes[type3-1]
# if UB settings exist for this angle type, set ubflag in aparams to 1 # if UB settings exist for this angle type, set ubflag in aparams to 1
if (c1,c2,c3) in ubdict: if (c1,c2,c3) in ubdict:
r1,r2 = ubdict[(c1,c2,c3)] r1,r2 = ubdict[(c1,c2,c3)]
pflag,ubflag,v1,v2,v3,v4,v5,v6 = aparams[itype] pflag,ubflag,v1,v2,v3,v4,v5,v6 = aparams[itype]
@ -1211,7 +1212,7 @@ for itype in range(len(aparams)):
aparams[itype] = (pflag,ubflag,v1,v2,v3,v4,v5,v6) aparams[itype] = (pflag,ubflag,v1,v2,v3,v4,v5,v6)
else: else:
r1,r2 = 2*[0.0] r1,r2 = 2*[0.0]
ubparams.append((r1,r2)) ubparams.append((r1,r2))
# generate dtype = LAMMPS type of each dihedral # generate dtype = LAMMPS type of each dihedral
@ -1243,7 +1244,7 @@ for atom1,atom2,atom3,atom4 in dlist:
c2 = classes[type2-1] c2 = classes[type2-1]
c3 = classes[type3-1] c3 = classes[type3-1]
c4 = classes[type4-1] c4 = classes[type4-1]
if (c1,c2,c3,c4) in ddict: m,params = ddict[(c1,c2,c3,c4)] if (c1,c2,c3,c4) in ddict: m,params = ddict[(c1,c2,c3,c4)]
elif (c4,c3,c2,c1) in ddict: m,params = ddict[(c4,c3,c2,c1)] elif (c4,c3,c2,c1) in ddict: m,params = ddict[(c4,c3,c2,c1)]
else: else:
@ -1289,11 +1290,11 @@ for atom1,atom2,atom3,atom4 in olist:
# 4-tuple is only an improper if matches an entry in PRM file # 4-tuple is only an improper if matches an entry in PRM file
# olist_reduced = list of just these 4-tuples # olist_reduced = list of just these 4-tuples
if (c1,c2) in odict: if (c1,c2) in odict:
m,params = odict[(c1,c2)] m,params = odict[(c1,c2)]
olist_reduced.append((atom1,atom2,atom3,atom4)) olist_reduced.append((atom1,atom2,atom3,atom4))
if not flags[m]: if not flags[m]:
oneparams = params[4:] oneparams = params[4:]
oparams.append(oneparams) oparams.append(oneparams)
@ -1319,7 +1320,7 @@ classes = prm.classes
pitdict = {} pitdict = {}
for m,params in enumerate(prm.pitorsionparams): for m,params in enumerate(prm.pitorsionparams):
pitdict[(params[0],params[1])] = (m,params) pitdict[(params[0],params[1])] = (m,params)
flags = len(prm.pitorsionparams)*[0] flags = len(prm.pitorsionparams)*[0]
pitorsiontype = [] pitorsiontype = []
pitorsionparams = [] pitorsionparams = []
@ -1333,7 +1334,7 @@ for tmp1,tmp2,atom1,atom2,tmp3,tmp4 in pitorsionlist:
# 6-tuple is only a PiTorsion if central 2 atoms match an entry in PRM file # 6-tuple is only a PiTorsion if central 2 atoms match an entry in PRM file
# pitorsionlist_reduced = list of just these 6-tuples # pitorsionlist_reduced = list of just these 6-tuples
if (c1,c2) in pitdict or (c2,c1) in pitdict: if (c1,c2) in pitdict or (c2,c1) in pitdict:
if (c1,c2) in pitdict: m,params = pitdict[(c1,c2)] if (c1,c2) in pitdict: m,params = pitdict[(c1,c2)]
else: m,params = pitdict[(c2,c1)] else: m,params = pitdict[(c2,c1)]
@ -1344,7 +1345,7 @@ for tmp1,tmp2,atom1,atom2,tmp3,tmp4 in pitorsionlist:
pitorsionparams.append(v1) pitorsionparams.append(v1)
flags[m] = len(pitorsionparams) flags[m] = len(pitorsionparams)
pitorsiontype.append(flags[m]) pitorsiontype.append(flags[m])
# replace original pitorsionlist with reduced version # replace original pitorsionlist with reduced version
pitorsionlist = pitorsionlist_reduced pitorsionlist = pitorsionlist_reduced
@ -1384,7 +1385,7 @@ for atom1,atom2,atom3,atom4,atom5 in bitorsionlist:
# 5-tuple is only a BiTorsion if 5 atoms match an entry in PRM file # 5-tuple is only a BiTorsion if 5 atoms match an entry in PRM file
# bitorsionlist_reduced = list of just these 5-tuples # bitorsionlist_reduced = list of just these 5-tuples
if (c1,c2,c3,c4,c5) in bitdict or (c5,c4,c3,c2,c1) in bitdict: if (c1,c2,c3,c4,c5) in bitdict or (c5,c4,c3,c2,c1) in bitdict:
if (c1,c2,c3,c4,c5) in bitdict: m,params = bitdict[(c1,c2,c3,c4,c5)] if (c1,c2,c3,c4,c5) in bitdict: m,params = bitdict[(c1,c2,c3,c4,c5)]
else: m,params = bitdict[(c5,c4,c3,c2,c1)] else: m,params = bitdict[(c5,c4,c3,c2,c1)]
@ -1395,7 +1396,7 @@ for atom1,atom2,atom3,atom4,atom5 in bitorsionlist:
bitorsionparams.append(v1) bitorsionparams.append(v1)
flags[m] = len(bitorsionparams) flags[m] = len(bitorsionparams)
bitorsiontype.append(flags[m]) bitorsiontype.append(flags[m])
# replace original bitorsionlist with reduced version # replace original bitorsionlist with reduced version
bitorsionlist = bitorsionlist_reduced bitorsionlist = bitorsionlist_reduced
@ -1436,7 +1437,7 @@ bitorsionlist = bitorsionlist_reduced
# if tgroup[j-1] > 0: continue # if tgroup[j-1] > 0: continue
# if ttype[j-1] not in polgroup[ttype[m-1]-1]: continue # if ttype[j-1] not in polgroup[ttype[m-1]-1]: continue
# stack.append(j) # stack.append(j)
# ---------------------------------------- # ----------------------------------------
# write LAMMPS data file via Pizza.py data class # write LAMMPS data file via Pizza.py data class
# ---------------------------------------- # ----------------------------------------
@ -1494,7 +1495,7 @@ d.sections["Tinker Types"] = lines
if nbonds: if nbonds:
d.headers["bonds"] = len(blist) d.headers["bonds"] = len(blist)
d.headers["bond types"] = len(bparams) d.headers["bond types"] = len(bparams)
lines = [] lines = []
for i,one in enumerate(bparams): for i,one in enumerate(bparams):
strone = [str(single) for single in one] strone = [str(single) for single in one]
@ -1502,7 +1503,7 @@ if nbonds:
lines.append(line+'\n') lines.append(line+'\n')
d.sections["Bond Coeffs"] = lines d.sections["Bond Coeffs"] = lines
lines = [] lines = []
for i,one in enumerate(blist): for i,one in enumerate(blist):
line = "%d %d %d %d" % (i+1,btype[i],one[0],one[1]) line = "%d %d %d %d" % (i+1,btype[i],one[0],one[1])
lines.append(line+'\n') lines.append(line+'\n')
@ -1511,7 +1512,7 @@ if nbonds:
if nangles: if nangles:
d.headers["angles"] = len(alist) d.headers["angles"] = len(alist)
d.headers["angle types"] = len(aparams) d.headers["angle types"] = len(aparams)
lines = [] lines = []
for i,one in enumerate(aparams): for i,one in enumerate(aparams):
strone = [str(single) for single in one] strone = [str(single) for single in one]
@ -1533,7 +1534,7 @@ if nangles:
lines.append(line+'\n') lines.append(line+'\n')
d.sections["UreyBradley Coeffs"] = lines d.sections["UreyBradley Coeffs"] = lines
lines = [] lines = []
for i,one in enumerate(alist): for i,one in enumerate(alist):
line = "%d %d %d %d %d" % (i+1,atype[i],one[0],one[1],one[2]) line = "%d %d %d %d %d" % (i+1,atype[i],one[0],one[1],one[2])
lines.append(line+'\n') lines.append(line+'\n')
@ -1542,7 +1543,7 @@ if nangles:
if ndihedrals: if ndihedrals:
d.headers["dihedrals"] = len(dlist) d.headers["dihedrals"] = len(dlist)
d.headers["dihedral types"] = len(dparams) d.headers["dihedral types"] = len(dparams)
lines = [] lines = []
for i,one in enumerate(dparams): for i,one in enumerate(dparams):
strone = [str(single) for single in one] strone = [str(single) for single in one]
@ -1550,7 +1551,7 @@ if ndihedrals:
lines.append(line+'\n') lines.append(line+'\n')
d.sections["Dihedral Coeffs"] = lines d.sections["Dihedral Coeffs"] = lines
lines = [] lines = []
for i,one in enumerate(dlist): for i,one in enumerate(dlist):
line = "%d %d %d %d %d %d" % (i+1,dtype[i],one[0],one[1],one[2],one[3]) line = "%d %d %d %d %d %d" % (i+1,dtype[i],one[0],one[1],one[2],one[3])
lines.append(line+'\n') lines.append(line+'\n')
@ -1559,7 +1560,7 @@ if ndihedrals:
if nimpropers: if nimpropers:
d.headers["impropers"] = len(olist) d.headers["impropers"] = len(olist)
d.headers["improper types"] = len(oparams) d.headers["improper types"] = len(oparams)
lines = [] lines = []
for i,one in enumerate(oparams): for i,one in enumerate(oparams):
strone = [str(single) for single in one] strone = [str(single) for single in one]
@ -1567,7 +1568,7 @@ if nimpropers:
lines.append(line+'\n') lines.append(line+'\n')
d.sections["Improper Coeffs"] = lines d.sections["Improper Coeffs"] = lines
lines = [] lines = []
for i,one in enumerate(olist): 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]) line = "%d %d %d %d %d %d" % (i+1,otype[i],one[0],one[1],one[2],one[3])
lines.append(line+'\n') lines.append(line+'\n')
@ -1576,7 +1577,7 @@ if nimpropers:
if npitorsions: if npitorsions:
d.headers["pitorsions"] = len(pitorsionlist) d.headers["pitorsions"] = len(pitorsionlist)
d.headers["pitorsion types"] = len(pitorsionparams) d.headers["pitorsion types"] = len(pitorsionparams)
lines = [] lines = []
for i,one in enumerate(pitorsionparams): for i,one in enumerate(pitorsionparams):
strone = [str(single) for single in one] strone = [str(single) for single in one]
@ -1584,7 +1585,7 @@ if npitorsions:
lines.append(line+'\n') lines.append(line+'\n')
d.sections["PiTorsion Coeffs"] = lines d.sections["PiTorsion Coeffs"] = lines
lines = [] lines = []
for i,one in enumerate(pitorsionlist): for i,one in enumerate(pitorsionlist):
line = "%d %d %d %d %d %d %d %d" % \ line = "%d %d %d %d %d %d %d %d" % \
(i+1,pitorsiontype[i],one[0],one[1],one[2],one[3],one[4],one[5]) (i+1,pitorsiontype[i],one[0],one[1],one[2],one[3],one[4],one[5])
@ -1595,7 +1596,7 @@ if nbitorsions:
d.headers["bitorsions"] = len(bitorsionlist) d.headers["bitorsions"] = len(bitorsionlist)
# if there are bitorsions, then -bitorsion file must have been specified # if there are bitorsions, then -bitorsion file must have been specified
if not bitorsionfile: if not bitorsionfile:
error("no -bitorsion file was specified, but %d bitorsions exist" % \ error("no -bitorsion file was specified, but %d bitorsions exist" % \
nbitorsions) nbitorsions)
@ -1613,8 +1614,8 @@ if nbitorsions:
xgrid,ygrid,value = array[ix][iy] xgrid,ygrid,value = array[ix][iy]
print(" ",xgrid,ygrid,value, file=fp) print(" ",xgrid,ygrid,value, file=fp)
fp.close() fp.close()
lines = [] lines = []
for i,one in enumerate(bitorsionlist): for i,one in enumerate(bitorsionlist):
line = "%d %d %d %d %d %d %d" % \ line = "%d %d %d %d %d %d %d" % \
(i+1,bitorsiontype[i],one[0],one[1],one[2],one[3],one[4]) (i+1,bitorsiontype[i],one[0],one[1],one[2],one[3],one[4])
@ -1629,7 +1630,7 @@ print("Natoms =",natoms)
print("Ntypes =",ntypes) print("Ntypes =",ntypes)
print("Tinker XYZ types =",len(tink2lmp)) print("Tinker XYZ types =",len(tink2lmp))
print("Tinker PRM types =",prm.ntypes) print("Tinker PRM types =",prm.ntypes)
#print "Tinker groups =",ngroups #print("Tinker groups =",ngroups)
print("Nmol =",nmol) print("Nmol =",nmol)
print("Nbonds =",nbonds) print("Nbonds =",nbonds)
print("Nangles =",nangles) print("Nangles =",nangles)