OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
waterReplacer
1#!@Python3_EXECUTABLE@
2"""Water Replacer
3
4Finds atomistic waters in an xyz file and generates an OpenMD (omd)
5file with center of mass and orientational coordinates for rigid body
6waters.
7
8Usage: waterReplacer
9
10Options:
11 -h, --help show this help
12 -x, use the specified input (.xyz) file
13 -o, --output-file=... use specified output (.omd) file
14 -s, --starting-index=. start the water objects with an index
15 -t, --tolerance=...... tolerance to add to ideal O-H bond length
16
17
18Example:
19 waterReplacer -x basal.xyz -o basal.omd
20
21"""
22
23__author__ = "Dan Gezelter (gezelter@nd.edu)"
24__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
25__license__ = "OpenMD"
26
27import sys
28import getopt
29import string
30import math
31import random
32import numpy
33
34_haveXYZFileName = 0
35_haveOutputFileName = 0
36
37atypes = []
38positions = []
39metaData = []
40frameData = []
41WaterPos = []
42WaterQuats = []
43indices = []
44Hmat = []
45BoxInv = []
46Eliminate = []
47
48#Hmat = zeros([3,3],Float)
49#BoxInv = zeros([3],Float)
50
51def usage():
52 print(__doc__)
53
54def readFile(XYZFileName):
55 print("reading XYZ file")
56
57 XYZFile = open(XYZFileName, 'r')
58 # Find number of atoms first
59 line = XYZFile.readline()
60 L = line.split()
61 nAtoms = int(L[0])
62 # skip comment line
63 line = XYZFile.readline()
64 for i in range(nAtoms):
65 line = XYZFile.readline()
66 L = line.split()
67 myIndex = i
68 indices.append(myIndex)
69 atomType = L[0]
70 atypes.append(atomType)
71 x = float(L[1])
72 y = float(L[2])
73 z = float(L[3])
74 positions.append([x, y, z])
75 XYZFile.close()
76
77def findWaters(covTol):
78 print("finding water molecules")
79 # simpler since we only have to find H atoms within a few
80 # angstroms of each water:
81 H = []
82 hCovRad = 0.32
83 oCovRad = 0.73
84 OHbond = hCovRad + oCovRad + covTol
85 Hmass = 1.0079
86 Omass = 15.9994
87
88 for i in range(len(indices)):
89 if (atypes[i] == "O"):
90 H.clear()
91 COM = [0.0, 0.0, 0.0]
92 opos = positions[i]
93 for j in range(len(indices)):
94 if (atypes[j] == "H"):
95 hpos = positions[j]
96 dx = opos[0] - hpos[0]
97 dy = opos[1] - hpos[1]
98 dz = opos[2] - hpos[2]
99 dist = math.sqrt(dx*dx + dy*dy + dz*dz)
100 if (dist < OHbond):
101 if (len(H) >= 2):
102 print("oxygen %d had too many hydrogens" % (i))
103 H.append(j)
104 if (len(H) != 2):
105 print("oxygen %d had %d hydrogens, skipping" % (i, len(H)))
106 if (len(H) == 2):
107 Xcom = Omass * opos[0] + Hmass*(positions[H[0]][0] + positions[H[1]][0])
108 Ycom = Omass * opos[1] + Hmass*(positions[H[0]][1] + positions[H[1]][1])
109 Zcom = Omass * opos[2] + Hmass*(positions[H[0]][2] + positions[H[1]][2])
110
111 totalMass = Omass + 2.0*Hmass
112 Xcom = Xcom / totalMass
113 Ycom = Ycom / totalMass
114 Zcom = Zcom / totalMass
115 COM = [Xcom, Ycom, Zcom]
116 WaterPos.append(COM)
117 bisector = [0.0, 0.0, 0.0]
118 ux = [0.0, 0.0, 0.0]
119 uy = [0.0, 0.0, 0.0]
120 uz = [0.0, 0.0, 0.0]
121 RotMat = numpy.zeros((3, 3), numpy.float64)
122
123 for j in range(3):
124 bisector[j] = 0.5*(positions[H[0]][j] + positions[H[1]][j])
125 uz[j] = bisector[j] - opos[j]
126 uy[j] = positions[H[0]][j] - positions[H[1]][j]
127
128 uz = normalize(uz)
129 uy = normalize(uy)
130 ux = cross(uy, uz)
131 ux = normalize(ux)
132
133 q = [0.0, 0.0, 0.0, 0.0]
134
135 # RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code:
136
137 RotMat[0] = ux
138 RotMat[1] = uy
139 RotMat[2] = uz
140
141 t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
142
143 if( t > 1e-6 ):
144 s = 0.5 / math.sqrt( t )
145 q[0] = 0.25 / s
146 q[1] = (RotMat[1][2] - RotMat[2][1]) * s
147 q[2] = (RotMat[2][0] - RotMat[0][2]) * s
148 q[3] = (RotMat[0][1] - RotMat[1][0]) * s
149 else:
150 ad1 = RotMat[0][0]
151 ad2 = RotMat[1][1]
152 ad3 = RotMat[2][2]
153
154 if( ad1 >= ad2 and ad1 >= ad3 ):
155 s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
156 q[0] = (RotMat[1][2] - RotMat[2][1]) * s
157 q[1] = 0.25 / s
158 q[2] = (RotMat[0][1] + RotMat[1][0]) * s
159 q[3] = (RotMat[0][2] + RotMat[2][0]) * s
160 elif ( ad2 >= ad1 and ad2 >= ad3 ):
161 s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
162 q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
163 q[1] = (RotMat[0][1] + RotMat[1][0]) * s
164 q[2] = 0.25 / s
165 q[3] = (RotMat[1][2] + RotMat[2][1]) * s
166 else:
167 s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
168 q[0] = (RotMat[0][1] - RotMat[1][0]) * s
169 q[1] = (RotMat[0][2] + RotMat[2][0]) * s
170 q[2] = (RotMat[1][2] + RotMat[2][1]) * s
171 q[3] = 0.25 / s
172
173 WaterQuats.append(q)
174 Eliminate.append(i)
175 Eliminate.append(H[0])
176 Eliminate.append(H[1])
177
178def writeFile(outputFileName, startingIndex):
179 outputFile = open(outputFileName, 'w')
180
181 outputFile.write("<OpenMD version=1>\n");
182
183 for metaline in metaData:
184 outputFile.write(metaline)
185
186 outputFile.write(" <Snapshot>\n")
187
188 for frameline in frameData:
189 outputFile.write(frameline)
190
191 outputFile.write(" <StuntDoubles>\n")
192
193 sdFormat = 'pvqj'
194
195 index = startingIndex
196 for i in range(len(WaterPos)):
197 outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (index, sdFormat, WaterPos[i][0], WaterPos[i][1], WaterPos[i][2], 0.0, 0.0, 0.0, WaterQuats[i][0], WaterQuats[i][1], WaterQuats[i][2], WaterQuats[i][3], 0.0, 0.0, 0.0))
198 index = index + 1
199
200
201 sdFormat = 'pv'
202 for i in range(len(indices)):
203 if i not in Eliminate:
204 outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e \n" % (index, sdFormat, positions[i][0], positions[i][1], positions[i][2], 0.0, 0.0, 0.0))
205 index = index + 1
206
207 outputFile.write(" </StuntDoubles>\n")
208 outputFile.write(" </Snapshot>\n")
209 outputFile.write("</OpenMD>\n")
210 outputFile.close()
211
212def dot(L1, L2):
213 myDot = 0.0
214 for i in range(len(L1)):
215 myDot = myDot + L1[i]*L2[i]
216 return myDot
217
218def normalize(L1):
219 L2 = []
220 myLength = math.sqrt(dot(L1, L1))
221 for i in range(len(L1)):
222 L2.append(L1[i] / myLength)
223 return L2
224
225def cross(L1, L2):
226 # don't call this with anything other than length 3 lists please
227 # or you'll be sorry
228 L3 = [0.0, 0.0, 0.0]
229 L3[0] = L1[1]*L2[2] - L1[2]*L2[1]
230 L3[1] = L1[2]*L2[0] - L1[0]*L2[2]
231 L3[2] = L1[0]*L2[1] - L1[1]*L2[0]
232 return L3
233
234def main(argv):
235 try:
236 opts, args = getopt.getopt(argv, "hx:o:s:t:", ["help", "xyz-file=", "output-file=", "starting-index=", "tolerance="])
237 except getopt.GetoptError:
238 usage()
239 sys.exit(2)
240 startingIndex = 0
241 tolerance = 0.45
242 for opt, arg in opts:
243 if opt in ("-h", "--help"):
244 usage()
245 sys.exit()
246 elif opt in ("-x", "--xyz-file"):
247 xyzFileName = arg
248 global _haveXYZFileName
249 _haveXYZFileName = 1
250 elif opt in ("-o", "--output-file"):
251 outputFileName = arg
252 global _haveOutputFileName
253 _haveOutputFileName = 1
254 elif opt in ("-s", "--starting-index"):
255 startingIndex = int(arg)
256 elif opt in ("-t", "--tolerance"):
257 tolerance = float(arg)
258 if (_haveXYZFileName != 1):
259 usage()
260 print("No XYZ file was specified")
261 sys.exit()
262 if (_haveOutputFileName != 1):
263 usage()
264 print("No output file was specified")
265 sys.exit()
266 readFile(xyzFileName)
267 findWaters(tolerance)
268 writeFile(outputFileName, startingIndex)
269
270if __name__ == "__main__":
271 if len(sys.argv) == 1:
272 usage()
273 sys.exit()
274 main(sys.argv[1:])