4Finds atomistic waters in an xyz file and generates an OpenMD (omd)
5file with center of mass and orientational coordinates for rigid body
11 -h, --help show this help
12 -x, use the specified input (.xyz) file
13 -o, --output-file=... use specified output (.omd) file
17 waterReplacer -x basal.xyz -o basal.omd
21__author__ = "Dan Gezelter (gezelter@nd.edu)"
22__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
33_haveOutputFileName = 0
46#Hmat = zeros([3,3],Float)
47#BoxInv = zeros([3],Float)
52def readFile(XYZFileName):
53 print("reading XYZ file")
55 XYZFile = open(XYZFileName, 'r')
56 # Find number of atoms first
57 line = XYZFile.readline()
61 line = XYZFile.readline()
62 for i in range(nAtoms):
63 line = XYZFile.readline()
66 indices.append(myIndex)
68 atypes.append(atomType)
72 positions.append([x, y, z])
76 print("finding water molecules")
77 # simpler since we only have to find H atoms within a few
78 # angstroms of each water:
83 OHbond = hCovRad + oCovRad + covTol
87 for i in range(len(indices)):
88 if (atypes[i] == "O"):
92 for j in range(len(indices)):
93 if (atypes[j] == "H"):
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)
101 print("oxygen %d had too many hydrogens" % (i))
104 print("oxygen %d had %d hydrogens, skipping" % (i, len(H)))
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])
110 totalMass = Omass + 2.0*Hmass
111 Xcom = Xcom / totalMass
112 Ycom = Ycom / totalMass
113 Zcom = Zcom / totalMass
114 COM = [Xcom, Ycom, Zcom]
116 bisector = [0.0, 0.0, 0.0]
120 RotMat = numpy.zeros((3, 3), numpy.float)
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]
132 q = [0.0, 0.0, 0.0, 0.0]
134 # RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code:
140 t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
143 s = 0.5 / math.sqrt( t )
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
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
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
164 q[3] = (RotMat[1][2] + RotMat[2][1]) * s
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
174 Eliminate.append(H[0])
175 Eliminate.append(H[1])
178def writeFile(outputFileName):
179 outputFile = open(outputFileName, 'w')
181 outputFile.write("<OpenMD version=1>\n");
183 for metaline in metaData:
184 outputFile.write(metaline)
186 outputFile.write(" <Snapshot>\n")
188 for frameline in frameData:
189 outputFile.write(frameline)
191 outputFile.write(" <StuntDoubles>\n")
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))
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))
207 outputFile.write(" </StuntDoubles>\n")
208 outputFile.write(" </Snapshot>\n")
209 outputFile.write("</OpenMD>\n")
214 for i in range(len(L1)):
215 myDot = myDot + L1[i]*L2[i]
220 myLength = math.sqrt(dot(L1, L1))
221 for i in range(len(L1)):
222 L2.append(L1[i] / myLength)
226 # don't call this with anything other than length 3 lists please
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]
236 opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz-file=", "output-file="])
237 except getopt.GetoptError:
240 for opt, arg in opts:
241 if opt in ("-h", "--help"):
244 elif opt in ("-x", "--xyz-file"):
246 global _haveXYZFileName
248 elif opt in ("-o", "--output-file"):
250 global _haveOutputFileName
251 _haveOutputFileName = 1
252 if (_haveXYZFileName != 1):
254 print("No XYZ file was specified")
256 if (_haveOutputFileName != 1):
258 print("No output file was specified")
260 readFile(xyzFileName)
262 writeFile(outputFileName)
264if __name__ == "__main__":
265 if len(sys.argv) == 1: