OpenMD 3.0
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
15
16Example:
17 waterReplacer -x basal.xyz -o basal.omd
18
19"""
20
21__author__ = "Dan Gezelter (gezelter@nd.edu)"
22__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
23__license__ = "OpenMD"
24
25import sys
26import getopt
27import string
28import math
29import random
30import numpy
31
32_haveXYZFileName = 0
33_haveOutputFileName = 0
34
35atypes = []
36positions = []
37metaData = []
38frameData = []
39WaterPos = []
40WaterQuats = []
41indices = []
42Hmat = []
43BoxInv = []
44Eliminate = []
45
46#Hmat = zeros([3,3],Float)
47#BoxInv = zeros([3],Float)
48
49def usage():
50 print(__doc__)
51
52def readFile(XYZFileName):
53 print("reading XYZ file")
54
55 XYZFile = open(XYZFileName, 'r')
56 # Find number of atoms first
57 line = XYZFile.readline()
58 L = line.split()
59 nAtoms = int(L[0])
60 # skip comment line
61 line = XYZFile.readline()
62 for i in range(nAtoms):
63 line = XYZFile.readline()
64 L = line.split()
65 myIndex = i
66 indices.append(myIndex)
67 atomType = L[0]
68 atypes.append(atomType)
69 x = float(L[1])
70 y = float(L[2])
71 z = float(L[3])
72 positions.append([x, y, z])
73 XYZFile.close()
74
75def findWaters():
76 print("finding water molecules")
77 # simpler since we only have to find H atoms within a few
78 # angstroms of each water:
79 H = []
80 hCovRad = 0.32
81 oCovRad = 0.73
82 covTol = 0.45
83 OHbond = hCovRad + oCovRad + covTol
84 Hmass = 1.0079
85 Omass = 15.9994
86
87 for i in range(len(indices)):
88 if (atypes[i] == "O"):
89 H.clear()
90 COM = [0.0, 0.0, 0.0]
91 opos = positions[i]
92 for j in range(len(indices)):
93 if (atypes[j] == "H"):
94 hpos = positions[j]
95 dx = opos[0] - hpos[0]
96 dy = opos[1] - hpos[1]
97 dz = opos[2] - hpos[2]
98 dist = math.sqrt(dx*dx + dy*dy + dz*dz)
99 if (dist < OHbond):
100 if (len(H) >= 2):
101 print("oxygen %d had too many hydrogens" % (i))
102 H.append(j)
103 if (len(H) != 2):
104 print("oxygen %d had %d hydrogens, skipping" % (i, len(H)))
105 if (len(H) == 2):
106 Xcom = Omass * opos[0] + Hmass*(positions[H[0]][0] + positions[H[1]][0])
107 Ycom = Omass * opos[1] + Hmass*(positions[H[0]][1] + positions[H[1]][1])
108 Zcom = Omass * opos[2] + Hmass*(positions[H[0]][2] + positions[H[1]][2])
109
110 totalMass = Omass + 2.0*Hmass
111 Xcom = Xcom / totalMass
112 Ycom = Ycom / totalMass
113 Zcom = Zcom / totalMass
114 COM = [Xcom, Ycom, Zcom]
115 WaterPos.append(COM)
116 bisector = [0.0, 0.0, 0.0]
117 ux = [0.0, 0.0, 0.0]
118 uy = [0.0, 0.0, 0.0]
119 uz = [0.0, 0.0, 0.0]
120 RotMat = numpy.zeros((3, 3), numpy.float)
121
122 for j in range(3):
123 bisector[j] = 0.5*(positions[H[0]][j] + positions[H[1]][j])
124 uz[j] = bisector[j] - opos[j]
125 uy[j] = positions[H[0]][j] - positions[H[1]][j]
126
127 uz = normalize(uz)
128 uy = normalize(uy)
129 ux = cross(uy, uz)
130 ux = normalize(ux)
131
132 q = [0.0, 0.0, 0.0, 0.0]
133
134 # RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code:
135
136 RotMat[0] = ux
137 RotMat[1] = uy
138 RotMat[2] = uz
139
140 t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
141
142 if( t > 1e-6 ):
143 s = 0.5 / math.sqrt( t )
144 q[0] = 0.25 / s
145 q[1] = (RotMat[1][2] - RotMat[2][1]) * s
146 q[2] = (RotMat[2][0] - RotMat[0][2]) * s
147 q[3] = (RotMat[0][1] - RotMat[1][0]) * s
148 else:
149 ad1 = RotMat[0][0]
150 ad2 = RotMat[1][1]
151 ad3 = RotMat[2][2]
152
153 if( ad1 >= ad2 and ad1 >= ad3 ):
154 s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
155 q[0] = (RotMat[1][2] - RotMat[2][1]) * s
156 q[1] = 0.25 / s
157 q[2] = (RotMat[0][1] + RotMat[1][0]) * s
158 q[3] = (RotMat[0][2] + RotMat[2][0]) * s
159 elif ( ad2 >= ad1 and ad2 >= ad3 ):
160 s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
161 q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
162 q[1] = (RotMat[0][1] + RotMat[1][0]) * s
163 q[2] = 0.25 / s
164 q[3] = (RotMat[1][2] + RotMat[2][1]) * s
165 else:
166 s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
167 q[0] = (RotMat[0][1] - RotMat[1][0]) * s
168 q[1] = (RotMat[0][2] + RotMat[2][0]) * s
169 q[2] = (RotMat[1][2] + RotMat[2][1]) * s
170 q[3] = 0.25 / s
171
172 WaterQuats.append(q)
173 Eliminate.append(i)
174 Eliminate.append(H[0])
175 Eliminate.append(H[1])
176
177
178def writeFile(outputFileName):
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 = 0
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:", ["help", "xyz-file=", "output-file="])
237 except getopt.GetoptError:
238 usage()
239 sys.exit(2)
240 for opt, arg in opts:
241 if opt in ("-h", "--help"):
242 usage()
243 sys.exit()
244 elif opt in ("-x", "--xyz-file"):
245 xyzFileName = arg
246 global _haveXYZFileName
247 _haveXYZFileName = 1
248 elif opt in ("-o", "--output-file"):
249 outputFileName = arg
250 global _haveOutputFileName
251 _haveOutputFileName = 1
252 if (_haveXYZFileName != 1):
253 usage()
254 print("No XYZ file was specified")
255 sys.exit()
256 if (_haveOutputFileName != 1):
257 usage()
258 print("No output file was specified")
259 sys.exit()
260 readFile(xyzFileName)
261 findWaters()
262 writeFile(outputFileName)
263
264if __name__ == "__main__":
265 if len(sys.argv) == 1:
266 usage()
267 sys.exit()
268 main(sys.argv[1:])