OpenMD 3.1
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
16
17Example:
18 waterReplacer -x basal.xyz -o basal.omd
19
20"""
21
22__author__ = "Dan Gezelter (gezelter@nd.edu)"
23__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
24__license__ = "OpenMD"
25
26import sys
27import getopt
28import string
29import math
30import random
31import numpy
32
33_haveXYZFileName = 0
34_haveOutputFileName = 0
35
36atypes = []
37positions = []
38metaData = []
39frameData = []
40WaterPos = []
41WaterQuats = []
42indices = []
43Hmat = []
44BoxInv = []
45Eliminate = []
46
47#Hmat = zeros([3,3],Float)
48#BoxInv = zeros([3],Float)
49
50def usage():
51 print(__doc__)
52
53def readFile(XYZFileName):
54 print("reading XYZ file")
55
56 XYZFile = open(XYZFileName, 'r')
57 # Find number of atoms first
58 line = XYZFile.readline()
59 L = line.split()
60 nAtoms = int(L[0])
61 # skip comment line
62 line = XYZFile.readline()
63 for i in range(nAtoms):
64 line = XYZFile.readline()
65 L = line.split()
66 myIndex = i
67 indices.append(myIndex)
68 atomType = L[0]
69 atypes.append(atomType)
70 x = float(L[1])
71 y = float(L[2])
72 z = float(L[3])
73 positions.append([x, y, z])
74 XYZFile.close()
75
76def findWaters():
77 print("finding water molecules")
78 # simpler since we only have to find H atoms within a few
79 # angstroms of each water:
80 H = []
81 hCovRad = 0.32
82 oCovRad = 0.73
83 covTol = 0.45
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.float_)
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
178
179def writeFile(outputFileName, startingIndex):
180 outputFile = open(outputFileName, 'w')
181
182 outputFile.write("<OpenMD version=1>\n");
183
184 for metaline in metaData:
185 outputFile.write(metaline)
186
187 outputFile.write(" <Snapshot>\n")
188
189 for frameline in frameData:
190 outputFile.write(frameline)
191
192 outputFile.write(" <StuntDoubles>\n")
193
194 sdFormat = 'pvqj'
195
196 index = startingIndex
197 for i in range(len(WaterPos)):
198 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))
199 index = index + 1
200
201
202 sdFormat = 'pv'
203 for i in range(len(indices)):
204 if i not in Eliminate:
205 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))
206 index = index + 1
207
208 outputFile.write(" </StuntDoubles>\n")
209 outputFile.write(" </Snapshot>\n")
210 outputFile.write("</OpenMD>\n")
211 outputFile.close()
212
213def dot(L1, L2):
214 myDot = 0.0
215 for i in range(len(L1)):
216 myDot = myDot + L1[i]*L2[i]
217 return myDot
218
219def normalize(L1):
220 L2 = []
221 myLength = math.sqrt(dot(L1, L1))
222 for i in range(len(L1)):
223 L2.append(L1[i] / myLength)
224 return L2
225
226def cross(L1, L2):
227 # don't call this with anything other than length 3 lists please
228 # or you'll be sorry
229 L3 = [0.0, 0.0, 0.0]
230 L3[0] = L1[1]*L2[2] - L1[2]*L2[1]
231 L3[1] = L1[2]*L2[0] - L1[0]*L2[2]
232 L3[2] = L1[0]*L2[1] - L1[1]*L2[0]
233 return L3
234
235def main(argv):
236 try:
237 opts, args = getopt.getopt(argv, "hx:o:s:", ["help", "xyz-file=", "output-file=", "starting-index="])
238 except getopt.GetoptError:
239 usage()
240 sys.exit(2)
241 startingIndex = 0
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 if (_haveXYZFileName != 1):
257 usage()
258 print("No XYZ file was specified")
259 sys.exit()
260 if (_haveOutputFileName != 1):
261 usage()
262 print("No output file was specified")
263 sys.exit()
264 readFile(xyzFileName)
265 findWaters()
266 writeFile(outputFileName, startingIndex)
267
268if __name__ == "__main__":
269 if len(sys.argv) == 1:
270 usage()
271 sys.exit()
272 main(sys.argv[1:])