Since I am currently looking at the solvation dynamics of a series of isomers of a metal oxide cluster I quickly got tired of working out the bonds and bond constraints for my 101 atom clusters by hand.
So here's a python (2.7.x) script that generates parameters suitable for a GROMACS .top or .itp file.
Let's call the script genparam. You can call it as follows:
genparam molecule.xyz 0.21 3
where molecule.xyz is the molecule in xyz coordinates. 0.21 (nm) is used to generate bonds -- atoms that are 0.21 nm or closer to each other are bonded. The final number should be 1, 2 or 3. If it is 1, then only bonds (and constraints for bonds that are 0.098 nm or shorter -- such as O-H bonds. It's specific for my use, so you can edit the code.) are printed. If it is 2, then bonds and angles are printed. If it is 3, then bonds, angles and dihedrals are printed.
A simple and understandable example using methanol looks like this:
./genparams.py methanol.xyz 0.15 3[bonds] [bonds] 1 2 6 0.143 50000.00 ;C OH 1 3 6 0.108 50000.00 ;C H 1 4 6 0.108 50000.00 ;C H 1 5 6 0.108 50000.00 ;C H [constraints] 2 6 2 0.096 ;OH HO [ angles ] 2 1 3 1 109.487 50000.00 ;OH C H 2 1 4 1 109.487 50000.00 ;OH C H 2 1 5 1 109.466 50000.00 ;OH C H 1 2 6 1 109.578 50000.00 ;C OH HO 3 1 4 1 109.542 50000.00 ;H C H 3 1 5 1 109.534 50000.00 ;H C H 4 1 5 1 109.534 50000.00 ;H C H [ dihedrals ] 6 2 1 3 1 60.01 50000.00 2 ;HO OH C H 6 2 1 4 1 60.01 50000.00 2 ;HO OH C H 6 2 1 5 1 180.00 50000.00 2 ;HO OH C H
where methanol.xyz looks like this:
6 Methanol C -0.000000010000 0.138569980000 0.355570700000 OH -0.000000010000 0.187935770000 -1.074466460000 H 0.882876920000 -0.383123830000 0.697839450000 H -0.882876940000 -0.383123830000 0.697839450000 H -0.000000010000 1.145042790000 0.750208830000 HO -0.000000010000 -0.705300580000 -1.426986340000
Note that the bond length value may need to be tuned for each molecule -- e.g. 0.21 nm is fine for my metal oxides which have long M-O bonds, while 0.15 nm is better for my methanol.xyz. In the end, you'll probably have to do some manual editing to remove excess bonds, angles and dihedrals.
Also note that there's a switch (prevent_hhbonds) which prevent the formation of bonds between atoms that start with H -- this is to prevent bonds between e.g. H in CH3 units (177 pm). You can change this in the code in the __main__ section.
Finally, note that the function types, multiplicities and forces may need to be edited. I suggest not doing that in the code, but in the output as it will depend on each molecule.
This script won't do science for you, but it'll take care of some of the drudgery of providing geometric descriptors.
Anyway, having said all that, here's the code:
#!/usr/bin/python
import sys
from math import sqrt, pi
from itertools import permutations
from math import acos,atan2
# see table 5.5 (p 132, v 4.6.3) in the gromacs manual
# for the different function types
#from
#http://stackoverflow.com/questions/1984799/cross-product-of-2-different-\
#vectors-in-python
def crossproduct(a, b):
c = [a[1]*b[2] - a[2]*b[1],
a[2]*b[0] - a[0]*b[2],
a[0]*b[1] - a[1]*b[0]]
return c
#end of code snippet
# mostly from
# http://www.emoticode.net/python/calculate-angle-between-two-vectors.html
def dotproduct(a,b):
return sum([a[i]*b[i] for i in range(len(a))])
def veclength(a):
length=sqrt(a[0]**2+a[1]**2+a[2]**2)
return length
def ange(a,b,la,lb,angle_unit):
dp=dotproduct(a,b)
costheta=dp/(la*lb)
angle=acos(costheta)
if angle_unit=='deg':
angle=angle*180/pi
return angle
# end of code snippet
def diheder(a,b,c,angle_unit):
dihedral=atan2(veclength(crossproduct(crossproduct(a,b),\
crossproduct(b,c))), dotproduct(crossproduct(a,b),crossproduct(b,c)))
if angle_unit=='deg':
dihedral=dihedral*180/pi
return dihedral
def readatoms(infile):
positions=[]
f=open(infile,'r')
atomno=-2
for line in f:
atomno+=1
if atomno >=1:
position=filter(None,line.rstrip('\n').split(' '))
if len(position)>3:
positions+=[[position[0],int(atomno),\
float(position[1]),float(position[2]),\
float(position[3])]]
return positions
def makebonds(positions,rcutoff,prevent_hhbonds):
bonds=[]
for firstatom in positions:
for secondatom in positions:
distance=round(sqrt((firstatom[2]-secondatom[2])**2\
+(firstatom[3]-secondatom[3])**2+(firstatom[4]-secondatom[4])**2)/10.0,3)
xyz=[[firstatom[2],firstatom[3],firstatom[4]],\
[secondatom[2],secondatom[3],secondatom[4]]]
if distance<=rcutoff and firstatom[1]!=secondatom[1]:
if prevent_hhbonds and (firstatom[0][0:1]=='H' and\
secondatom[0][0:1]=='H'):
pass
elif firstatom[1]<secondatom[1]:
bonds+=[[firstatom[1],secondatom[1],\
distance,firstatom[0],secondatom[0],xyz[0],xyz[1]]]
else:
bonds+=[[secondatom[1],firstatom[1],\
distance,firstatom[0],secondatom[0],xyz[1],xyz[0]]]
return bonds
def dedupe_bonds(bonds):
newbonds=[]
for olditem in bonds:
dupe=False
for newitem in newbonds:
if newitem[0]==olditem[0] and newitem[1]==olditem[1]:
dupe=True
break;
if dupe==False:
newbonds+=[olditem]
return(newbonds)
def genvec(a,b):
vec=[b[0]-a[0],b[1]-a[1],b[2]-a[2]]
return vec
def findangles(bonds,angle_unit):
# for atoms 1,2,3 we can have the following situations
# 1-2, 1-3
# 1-2, 2-3
# 1-3, 2-3
# The indices are sorted so that the lower number is always first
angles=[]
for firstbond in bonds:
for secondbond in bonds:
if firstbond[0]==secondbond[0] and not (firstbond[1]==\
secondbond[1]): # 1-2, 1-3
vec=[genvec(firstbond[6],firstbond[5])]
vec+=[genvec(secondbond[6],secondbond[5])]
angle=ange(vec[0],vec[1],firstbond[2]*10,\
secondbond[2]*10,angle_unit)
angles+=[[firstbond[1],firstbond[0],\
secondbond[1],angle,firstbond[4],firstbond[3],secondbond[4],firstbond[6],\
firstbond[5],secondbond[6]]]
if firstbond[0]==secondbond[1] and not (firstbond[1]==\
secondbond[1]): # 1-2, 3-1
#this should never be relevant since we've sorted the atom numbers
pass
if firstbond[1]==secondbond[0] and not (firstbond[0]==\
secondbond[1]): # 1-2, 2-3
vec=[genvec(firstbond[5],firstbond[6])]
vec+=[genvec(secondbond[6],secondbond[5])]
angle=ange(vec[0],vec[1],firstbond[2]*10,\
secondbond[2]*10,angle_unit)
angles+=[[firstbond[0],firstbond[1],\
secondbond[1],angle,firstbond[3],firstbond[4],secondbond[4],firstbond[5],\
firstbond[6],secondbond[6]]]
if firstbond[1]==secondbond[1] and not (firstbond[0]==\
secondbond[0]): # 1-3, 2-3
vec=[genvec(firstbond[6],firstbond[5])]
vec+=[genvec(secondbond[6],secondbond[5])]
angle=ange(vec[0],vec[1],firstbond[2]*10,\
secondbond[2]*10,angle_unit)
angles+=[[firstbond[0],firstbond[1],\
secondbond[0],angle,firstbond[3],firstbond[4],secondbond[3],firstbond[5],\
firstbond[6],secondbond[5]]]
return angles
def dedupe_angles(angles):
dupe=False
newangles=[]
for item in angles:
dupe=False
for anotheritem in newangles:
if item[0]==anotheritem[2] and (item[2]==anotheritem[0]\
and item[1]==anotheritem[1]):
dupe=True
break
if dupe==False:
newangles+=[item]
newerangles=[]
dupe=False
for item in newangles:
dupe=False
for anotheritem in newerangles:
if item[2]==anotheritem[2] and (item[0]==anotheritem[1]\
and item[1]==anotheritem[0]):
dupe_item=anotheritem
dupe=True
break
if dupe==False:
newerangles+=[item]
elif dupe==True:
if dupe_item[3]>item[3]:
pass;
else:
newerangles[len(newerangles)-1]=item
newestangles=[]
dupe=False
for item in newerangles:
dupe=False
for anotheritem in newestangles:
if (sorted(item[0:3]) == sorted (anotheritem[0:3])):
dupe_item=anotheritem
dupe=True
break
if dupe==False:
newestangles+=[item]
elif dupe==True:
if dupe_item[3]>item[3]:
pass;
else:
newestangles[len(newestangles)-1]=item
return newestangles
def finddihedrals(angles,bonds,angle_unit):
dihedrals=[]
for item in angles:
for anotheritem in bonds:
if item[2]==anotheritem[0]:
vec=[genvec(item[7],item[8])]
vec+=[genvec(item[8],item[9])]
vec+=[genvec(item[9],anotheritem[6])]
dihedral=diheder(vec[0],vec[1],vec[2],\
angle_unit)
dihedrals+=[ [ item[0],item[1],item[2],\
anotheritem[1], dihedral, item[4],item[5],item[6], anotheritem[4] ] ]
if item[0]==anotheritem[0] and not item[1]==anotheritem[1]:
vec=[genvec(anotheritem[6],item[7])]
vec+=[genvec(item[7],item[8])]
vec+=[genvec(item[8],item[9])]
dihedral=diheder(vec[0],vec[1],vec[2],\
angle_unit)
dihedrals+=[ [ anotheritem[1],item[0],item[1],\
item[2], dihedral, anotheritem[4],item[4],item[5], item[6] ] ]
return dihedrals
def dedup_dihedrals(dihedrals):
newdihedrals=[]
for item in dihedrals:
dupe=False
for anotheritem in newdihedrals:
rev=anotheritem[0:4]
rev.reverse()
if item[0:4] == rev:
dupe=True
if not dupe:
newdihedrals+=[item]
return newdihedrals
def print_bonds(bonds,func):
constraints=""
funcconstr='2'
print '[bonds]'
for item in bonds:
if item[2]<=0.098:
constraints+=(str(item[0])+'\t'+str(item[1])+'\t'+\
funcconstr+'\t'+str("%1.3f"% item[2])+'\t;'+str(item[3])+'\t'+str(item[4])+'\n')
else:
print str(item[0])+'\t'+str(item[1])+'\t'+func+'\t'+\
str("%1.3f"% item[2])+'\t'+'50000.00'+'\t;'+str(item[3])+'\t'+str(item[4])
print '[constraints]'
print constraints
return 0
def print_angles(angles,func):
print '[ angles ]'
for item in angles:
print str(item[0])+'\t'+str(item[1])+'\t'+str(item[2])+'\t'+\
func+'\t'+str("%3.3f" % item[3])+'\t'+'50000.00'+'\t;'\
+str(item[4])+'\t'+str(item[5])+'\t'+str(item[6])
return 0
def print_dihedrals(dihedrals,func):
print "[ dihedrals ]"
force='50000.00'
mult='2'
for item in dihedrals:
print str(item[0])+'\t'+str(item[1])+'\t'+str(item[2])+'\t'+\
str(item[3])+\
'\t'+func+'\t'+str("%3.2f" % item[4])+'\t'+force+'\t'+mult+\
'\t;'+str(item[5])+'\t'+str(item[6])+\
'\t'+str(item[7])+'\t'+str(item[8])
return 0
if __name__ == "__main__":
infile=sys.argv[1]
rcutoff=float(sys.argv[2]) # in nm
itemstoprint=int(sys.argv[3])
angle_unit='deg' #{'rad'|'deg'}
prevent_hhbonds=True # False is safer -- it prevents bonds between
# atoms whose names start with H
positions=readatoms(infile)
bonds=makebonds(positions,rcutoff,prevent_hhbonds)
bonds=dedupe_bonds(bonds)
print_bonds(bonds,'6')
angles=findangles(bonds,angle_unit)
angles=dedupe_angles(angles)
if itemstoprint>=2:
print_angles(angles,'1')
dihedrals=finddihedrals(angles,bonds,angle_unit)
dihedrals=dedup_dihedrals(dihedrals)
if itemstoprint>=3:
print_dihedrals(dihedrals,'1')
