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

modernize tinker2lmp scripts so they run with both python2.7 and python3.x

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-07-05 08:02:34 -04:00
parent 89e0989522
commit 212048d4cc
2 changed files with 130 additions and 129 deletions
Download Patch File Download Diff File Expand all files Collapse all files

64
tools/tinker/data.py
View File

@ -8,6 +8,7 @@
# data tool # data tool
from __future__ import print_function
oneline = "Read, write, manipulate LAMMPS data files" oneline = "Read, write, manipulate LAMMPS data files"
docstr = """ docstr = """
@ -79,7 +80,7 @@ except: PIZZA_GUNZIP = "gunzip"
# Class definition # Class definition
class data: class data(object):
# -------------------------------------------------------------------- # --------------------------------------------------------------------
@ -129,15 +130,14 @@ class data:
keyword,length = pair[0],pair[1] keyword,length = pair[0],pair[1]
if keyword == line: if keyword == line:
found = 1 found = 1
if not headers.has_key(length): if length not in headers:
raise StandardError, \ raise (Exception, "data section %s has no matching header value" % line)
"data section %s has no matching header value" % line
f.readline() f.readline()
list = [] list = []
for i in xrange(headers[length]): list.append(f.readline()) for i in range(headers[length]): list.append(f.readline())
sections[keyword] = list sections[keyword] = list
if not found: if not found:
raise StandardError,"invalid section %s in data file" % line raise (Exception,"invalid section %s in data file" % line)
f.readline() f.readline()
line = f.readline() line = f.readline()
if not line: if not line:
@ -153,7 +153,7 @@ class data:
def map(self,*pairs): def map(self,*pairs):
if len(pairs) % 2 != 0: if len(pairs) % 2 != 0:
raise StandardError, "data map() requires pairs of mappings" raise Exception("data map() requires pairs of mappings")
for i in range(0,len(pairs),2): for i in range(0,len(pairs),2):
j = i + 1 j = i + 1
self.names[pairs[j]] = pairs[i]-1 self.names[pairs[j]] = pairs[i]-1
@ -168,7 +168,7 @@ class data:
lines = self.sections[field] lines = self.sections[field]
for line in lines: for line in lines:
words = line.split() words = line.split()
values = map(float,words) values = list(map(float,words))
array.append(values) array.append(values)
return array return array
elif len(list) == 2: elif len(list) == 2:
@ -181,7 +181,7 @@ class data:
vec.append(float(words[n])) vec.append(float(words[n]))
return vec return vec
else: else:
raise StandardError, "invalid arguments for data.get()" raise Exception("invalid arguments for data.get()")
# -------------------------------------------------------------------- # --------------------------------------------------------------------
# reorder columns in a data file field # reorder columns in a data file field
@ -192,10 +192,10 @@ class data:
oldlines = self.sections[name] oldlines = self.sections[name]
newlines = natoms*[""] newlines = natoms*[""]
for index in order: for index in order:
for i in xrange(len(newlines)): for i in range(len(newlines)):
words = oldlines[i].split() words = oldlines[i].split()
newlines[i] += words[index-1] + " " newlines[i] += words[index-1] + " "
for i in xrange(len(newlines)): for i in range(len(newlines)):
newlines[i] += "\n" newlines[i] += "\n"
self.sections[name] = newlines self.sections[name] = newlines
@ -206,7 +206,7 @@ class data:
lines = self.sections[name] lines = self.sections[name]
newlines = [] newlines = []
j = icol - 1 j = icol - 1
for i in xrange(len(lines)): for i in range(len(lines)):
line = lines[i] line = lines[i]
words = line.split() words = line.split()
words[j] = str(vector[i]) words[j] = str(vector[i])
@ -229,7 +229,7 @@ class data:
self.replace("Atoms",self.names['y']+1,y) self.replace("Atoms",self.names['y']+1,y)
self.replace("Atoms",self.names['z']+1,z) self.replace("Atoms",self.names['z']+1,z)
if dm.names.has_key("ix") and self.names.has_key("ix"): if "ix" in dm.names and "ix" in self.names:
ix,iy,iz = dm.vecs(ntime,"ix","iy","iz") ix,iy,iz = dm.vecs(ntime,"ix","iy","iz")
self.replace("Atoms",self.names['ix']+1,ix) self.replace("Atoms",self.names['ix']+1,ix)
self.replace("Atoms",self.names['iy']+1,iy) self.replace("Atoms",self.names['iy']+1,iy)
@ -240,36 +240,36 @@ class data:
def delete(self,keyword): def delete(self,keyword):
if self.headers.has_key(keyword): del self.headers[keyword] if keyword in self.headers: del self.headers[keyword]
elif self.sections.has_key(keyword): del self.sections[keyword] elif keyword in self.sections: del self.sections[keyword]
else: raise StandardError, "keyword not found in data object" else: raise Exception("keyword not found in data object")
# -------------------------------------------------------------------- # --------------------------------------------------------------------
# write out a LAMMPS data file # write out a LAMMPS data file
def write(self,file): def write(self,file):
f = open(file,"w") f = open(file,"w")
print >>f,self.title print(self.title, file=f)
# write any keywords in standard list hkeywords # write any keywords in standard list hkeywords
# in the order they are in hkeywords # in the order they are in hkeywords
# then write any extra keywords at end of header section # then write any extra keywords at end of header section
for keyword in hkeywords: for keyword in hkeywords:
if self.headers.has_key(keyword): if keyword in self.headers:
if keyword == "xlo xhi" or keyword == "ylo yhi" or \ if keyword == "xlo xhi" or keyword == "ylo yhi" or \
keyword == "zlo zhi": keyword == "zlo zhi":
pair = self.headers[keyword] pair = self.headers[keyword]
print >>f,pair[0],pair[1],keyword print(pair[0],pair[1],keyword, file=f)
elif keyword == "xy xz yz": elif keyword == "xy xz yz":
triple = self.headers[keyword] triple = self.headers[keyword]
print >>f,triple[0],triple[1],triple[2],keyword print(triple[0],triple[1],triple[2],keyword, file=f)
else: else:
print >>f,self.headers[keyword],keyword print(self.headers[keyword],keyword, file=f)
for keyword in self.headers.keys(): for keyword in list(self.headers.keys()):
if keyword not in hkeywords: if keyword not in hkeywords:
print >>f,self.headers[keyword],keyword print(self.headers[keyword],keyword, file=f)
# write any sections in standard list skeywords # write any sections in standard list skeywords
# in the order they are in skeywords # in the order they are in skeywords
@ -277,18 +277,18 @@ class data:
for pair in skeywords: for pair in skeywords:
keyword = pair[0] keyword = pair[0]
if self.sections.has_key(keyword): if keyword in self.sections:
print >>f,"\n%s\n" % keyword print("\n%s\n" % keyword, file=f)
for line in self.sections[keyword]: for line in self.sections[keyword]:
print >>f,line, print(line, end='', file=f)
skeyfirst = [pair[0] for pair in skeywords] skeyfirst = [pair[0] for pair in skeywords]
for keyword in self.sections.keys(): for keyword in list(self.sections.keys()):
if keyword not in skeyfirst: if keyword not in skeyfirst:
print >>f,"\n%s\n" % keyword print("\n%s\n" % keyword, file=f)
for line in self.sections[keyword]: for line in self.sections[keyword]:
print >>f,line, print(line, end='', file=f)
f.close() f.close()
@ -304,13 +304,13 @@ class data:
def findtime(self,n): def findtime(self,n):
if n == 0: return 0 if n == 0: return 0
raise StandardError, "no step %d exists" % (n) raise(Exception, "no step %d exists" % (n))
# -------------------------------------------------------------------- # --------------------------------------------------------------------
# return list of atoms and bonds to viz for data object # return list of atoms and bonds to viz for data object
def viz(self,isnap): def viz(self,isnap):
if isnap: raise StandardError, "cannot call data.viz() with isnap != 0" if isnap: raise Exception("cannot call data.viz() with isnap != 0")
id = self.names["id"] id = self.names["id"]
type = self.names["type"] type = self.names["type"]
@ -336,7 +336,7 @@ class data:
# assumes atoms are sorted so can lookup up the 2 atoms in each bond # assumes atoms are sorted so can lookup up the 2 atoms in each bond
bonds = [] bonds = []
if self.sections.has_key("Bonds"): if "Bonds" in self.sections:
bondlines = self.sections["Bonds"] bondlines = self.sections["Bonds"]
for line in bondlines: for line in bondlines:
words = line.split() words = line.split()

113
tools/tinker/tinker2lmp.py
View File

@ -15,6 +15,7 @@
# Author: Steve Plimpton # Author: Steve Plimpton
from __future__ import print_function
import sys,os,math import sys,os,math
from data import data from data import data
@ -29,20 +30,20 @@ DELTA = 0.001 # delta on LAMMPS shrink-wrap box size, in Angstroms
def error(txt=""): def error(txt=""):
if not txt: if not txt:
print "Syntax: tinker2lmp.py -switch args ..." print("Syntax: tinker2lmp.py -switch args ...")
print " -xyz file" print(" -xyz file")
print " -amoeba file" print(" -amoeba file")
print " -hippo file" print(" -hippo file")
print " -data file" print(" -data file")
print " -bitorsion file" print(" -bitorsion file")
print " -nopbc" print(" -nopbc")
print " -pbc xhi yhi zhi" print(" -pbc xhi yhi zhi")
else: print "ERROR:",txt else: print("ERROR:",txt)
#sys.exit() #sys.exit()
# read and store values from a Tinker xyz file # read and store values from a Tinker xyz file
class XYZfile: class XYZfile(object):
def __init__(self,file): def __init__(self,file):
lines = open(file,'r').readlines() lines = open(file,'r').readlines()
header = lines[0] header = lines[0]
@ -212,7 +213,7 @@ class XYZfile:
def output(self,outfile): def output(self,outfile):
fp = open(outfile,'w') fp = open(outfile,'w')
words = self.header.split() words = self.header.split()
print >>fp,self.natoms,"replicated",' '.join(words[1:]) print(self.natoms,"replicated",' '.join(words[1:]), file=fp)
id = self.id id = self.id
label = self.label label = self.label
@ -225,9 +226,9 @@ class XYZfile:
# 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 >>fp,i+1,label[i],x[i],y[i],z[i],type[i], print(i+1,label[i],x[i],y[i],z[i],type[i], end=' ', file=fp)
for j in bonds[i]: print >>fp,j, for j in bonds[i]: print(j, end=' ', file=fp)
print >>fp print(file=fp)
fp.close() fp.close()
@ -255,7 +256,7 @@ class XYZfile:
# scalar force field params in Force Field Definition section # scalar force field params in Force Field Definition section
# bond, angle, dihedral coeffs indexed by Tinker classes # bond, angle, dihedral coeffs indexed by Tinker classes
class PRMfile: class PRMfile(object):
def __init__(self,file): def __init__(self,file):
lines = open(file,'r').readlines() lines = open(file,'r').readlines()
self.nlines = len(lines) self.nlines = len(lines)
@ -519,7 +520,7 @@ class PRMfile:
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))
mfourier = (len(words)-5) / 3 mfourier = int((len(words)-5)/3)
oneparams = [class1,class2,class3,class4,mfourier] oneparams = [class1,class2,class3,class4,mfourier]
for iset in range(mfourier): for iset in range(mfourier):
@ -743,7 +744,7 @@ if pbcflag:
else: else:
xlo = ylo = zlo = BIG xlo = ylo = zlo = BIG
xhi = yhi = zhi = -BIG xhi = yhi = zhi = -BIG
for i in xrange(natoms): for i in range(natoms):
xlo = min(xlo,float(x[i])) xlo = min(xlo,float(x[i]))
ylo = min(ylo,float(y[i])) ylo = min(ylo,float(y[i]))
zlo = min(zlo,float(z[i])) zlo = min(zlo,float(z[i]))
@ -1097,11 +1098,11 @@ for i,one in enumerate(alist):
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 #sys.exit() NOTE: allow this for now
if hcount == 0: which = 1 if hcount == 0: which = 1
@ -1109,22 +1110,22 @@ for i,one in enumerate(alist):
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 #sys.exit() NOTE: allow this for now
if hcount == 0: which = 1 if hcount == 0: which = 1
@ -1170,7 +1171,7 @@ 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) 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))
@ -1600,17 +1601,17 @@ if nbitorsions:
nbitorsions) nbitorsions)
fp = open(bitorsionfile,'w') fp = open(bitorsionfile,'w')
print >>fp,"Tinker BiTorsion parameter file for fix bitorsion\n" print("Tinker BiTorsion parameter file for fix bitorsion\n", file=fp)
print >>fp,"%d bitorsion types" % len(bitorsionparams) print("%d bitorsion types" % len(bitorsionparams), file=fp)
itype = 0 itype = 0
for nx,ny,array in bitorsionparams: for nx,ny,array in bitorsionparams:
itype += 1 itype += 1
print >>fp print(file=fp)
print >>fp,itype,nx,ny print(itype,nx,ny, file=fp)
for ix in range(nx): for ix in range(nx):
for iy in range(ny): for iy in range(ny):
xgrid,ygrid,value = array[ix][iy] xgrid,ygrid,value = array[ix][iy]
print >>fp," ",xgrid,ygrid,value print(" ",xgrid,ygrid,value, file=fp)
fp.close() fp.close()
lines = [] lines = []
@ -1624,21 +1625,21 @@ d.write(datafile)
# print stats to screen # print stats to screen
print "Natoms =",natoms 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)
print "Ndihedrals =",ndihedrals print("Ndihedrals =",ndihedrals)
print "Nimpropers =",nimpropers print("Nimpropers =",nimpropers)
print "Npitorsions =",npitorsions print("Npitorsions =",npitorsions)
print "Nbitorsions =",nbitorsions print("Nbitorsions =",nbitorsions)
print "Nbondtypes =",len(bparams) print("Nbondtypes =",len(bparams))
print "Nangletypes =",len(aparams) print("Nangletypes =",len(aparams))
print "Ndihedraltypes =",len(dparams) print("Ndihedraltypes =",len(dparams))
print "Nimpropertypes =",len(oparams) print("Nimpropertypes =",len(oparams))
print "Npitorsiontypes =",len(pitorsionparams) print("Npitorsiontypes =",len(pitorsionparams))
print "Nbitorsiontypes =",len(bitorsionparams) print("Nbitorsiontypes =",len(bitorsionparams))