OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
pack2omd
1#!@Python3_EXECUTABLE@
2"""Packmol RigidBody Replacer
3
4Finds atomistic rigid bodies in a packmol-generated xyz file and
5generates an OpenMD (omd) file with center of mass and orientational
6coordinates for rigid bodies.
7
8Usage: pack2omd
9
10Options:
11 -h, --help show this help
12 -x, use the specified packmol (.xyz) file
13 -r, --rigid-body=... use this xyz structure as the rigid body
14 -o, --output-file=... use specified output (.omd) file
15
16
17Example:
18 pack2omd -x tolueneBox.xyz -r singleToluene.xyz -o tolueneBox.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_haveRBFileName = 0
35_haveOutputFileName = 0
36
37atypes = []
38positions = []
39metaData = []
40frameData = []
41RBpos = []
42RBQuats = []
43RBindices = []
44RBatypes = []
45RBpositions = []
46indices = []
47Hmat = []
48BoxInv = []
49H = []
50Eliminate = []
51
52mass_table = {
53 'H': 1.00794,
54 'C': 12.0107,
55 'Cl': 35.453,
56 'O': 15.999,
57 'N': 14.007,
58 'S': 32.0655,
59 'Au': 196.466569,
60 'CH4': 16.05,
61 'CH3': 15.04,
62 'CH2': 14.03,
63 'CH': 13.02,
64 'CHar': 13.02,
65 'CHa': 13.02,
66 'RCHar': 12.0107,
67 'RCH': 12.0107,
68 'CH3S': 15.04,
69 'CH2S': 14.03,
70 'CHS': 13.02,
71 'CS': 12.0107,
72 'SYZ': 32.0655,
73 'SH': 32.0655,
74 'HS': 1.0079,
75 'S': 32.0655,
76 'SZ': 32.0655,
77 'SS': 32.0655,
78 'SP': 32.0655,
79 'CS': 12.0107,
80 'SAu': 228.9807,
81 'SSD': 18.0153,
82 'SSD1': 18.0153,
83 'SSD_E': 18.0153,
84 'SSD_RF': 18.0153,
85 'O_TIP3P': 15.9994,
86 'O_TIP4P': 15.9994,
87 'O_TIP4P-Ew': 15.9994,
88 'O_TIP5P': 15.9994,
89 'O_TIP5P-E': 15.9994,
90 'O_SPCE': 15.9994,
91 'O_SPC': 15.9994,
92 'H_TIP3P': 1.0079,
93 'H_TIP4P': 1.0079,
94 'H_TIP4P-Ew': 1.0079,
95 'H_TIP5P': 1.0079,
96 'H_SPCE': 1.0079,
97 'H_SPC': 1.0079,
98 'EP_TIP4P': 0.0,
99 'EP_TIP4P-Ew':0.0,
100 'EP_TIP5P': 0.0,
101 'Ni': 58.710,
102 'Cu': 63.550,
103 'Rh': 102.90550,
104 'Pd': 106.42,
105 'Ag': 107.8682,
106 'Ir': 192.217,
107 'Pt': 195.09
108}
109
110#Hmat = zeros([3,3],Float)
111#BoxInv = zeros([3],Float)
112
113def usage():
114 print(__doc__)
115
116def readFile(XYZFileName):
117 print("reading XYZ file")
118
119 XYZFile = open(XYZFileName, 'r')
120 # Find number of atoms first
121 line = XYZFile.readline()
122 L = line.split()
123 nAtoms = int(L[0])
124 # skip comment line
125 line = XYZFile.readline()
126 for i in range(nAtoms):
127 line = XYZFile.readline()
128 L = line.split()
129 myIndex = i
130 indices.append(myIndex)
131 atomType = L[0]
132 atypes.append(atomType)
133 x = float(L[1])
134 y = float(L[2])
135 z = float(L[3])
136 positions.append([x, y, z])
137 XYZFile.close()
138
139def readRBFile(RBFileName):
140 print("reading Rigid Body file")
141
142 RBFile = open(RBFileName, 'r')
143 # Find number of atoms first
144 line = RBFile.readline()
145 L = line.split()
146 nAtoms = int(L[0])
147 # skip comment line
148 line = RBFile.readline()
149 for i in range(nAtoms):
150 line = RBFile.readline()
151 L = line.split()
152 myIndex = i
153 RBindices.append(myIndex)
154 atomType = L[0]
155 RBatypes.append(atomType)
156 x = float(L[1])
157 y = float(L[2])
158 z = float(L[3])
159 RBpositions.append([x, y, z])
160 RBFile.close()
161
162def findCOM():
163 #find center of mass
164 Xcom = 0.0
165 Ycom = 0.0
166 Zcom = 0.0
167 totalMass = 0.0
168
169 for i in range(0, len(RBindices)):
170 myMass = mass_table[RBatypes[i]]
171
172 Xcom = Xcom + myMass * RBpositions[i][0]
173 Ycom = Ycom + myMass * RBpositions[i][1]
174 Zcom = Zcom + myMass * RBpositions[i][2]
175 totalMass = totalMass + myMass
176
177 Xcom = Xcom / totalMass
178 Ycom = Ycom / totalMass
179 Zcom = Zcom / totalMass
180
181 COM = [Xcom, Ycom, Zcom]
182
183 return COM
184
185def calcMoments():
186
187 COM = findCOM()
188
189 #find inertia tensor matrix elements
190
191 I = numpy.zeros((3, 3), numpy.float)
192
193 for i in range(0, len(RBindices)):
194 myMass = mass_table[RBatypes[i]]
195
196 # move the origin of the reference coordinate system to the COM
197 RBpositions[i][0] -= COM[0]
198 RBpositions[i][1] -= COM[1]
199 RBpositions[i][2] -= COM[2]
200
201 dx = RBpositions[i][0]
202 dy = RBpositions[i][1]
203 dz = RBpositions[i][2]
204
205 I[0, 0] = I[0, 0] + myMass * ( dy * dy + dz * dz )
206 I[1, 1] = I[1, 1] + myMass * ( dx * dx + dz * dz )
207 I[2, 2] = I[2, 2] + myMass * ( dx * dx + dy * dy )
208
209 I[0, 1] = I[0, 1] - myMass * ( dx * dy )
210 I[0, 2] = I[0, 2] - myMass * ( dx * dz )
211
212 I[1, 2] = I[1, 2] - myMass * ( dy * dz )
213
214 I[1, 0] = I[0, 1]
215 I[2, 0] = I[0, 2]
216 I[2, 1] = I[1, 2]
217
218 print("Inertia Tensor:")
219 print(I)
220 print()
221
222 (evals, evects) = numpy.linalg.eig(I)
223 print("evals:")
224 print(evals)
225 print()
226 print("evects:")
227 print(evects)
228 print()
229
230 return (COM, evals, evects)
231
232def findRBs():
233 ref_ = []
234 mov = []
235
236 for i in range(len(RBindices)):
237 ref_.append(numpy.array([RBpositions[i][0], RBpositions[i][1], RBpositions[i][2]], numpy.float))
238 mov.append(numpy.array([0, 0, 0], numpy.float))
239
240 print("finding rigid bodies (assuming strict packmol ordering)")
241
242 nBodies = int( len(indices) / len(RBindices))
243 print(nBodies)
244 xyzIndex = 0
245
246 for j in range(nBodies):
247 mov_com = numpy.zeros(3, numpy.float)
248 totalMass = 0.0
249
250 for i in range(len(RBindices)):
251 mov[i] = numpy.array([positions[xyzIndex][0], positions[xyzIndex][1], positions[xyzIndex][2]], numpy.float)
252 myMass = mass_table[RBatypes[i]]
253 mov_com = mov_com + myMass*mov[i]
254 totalMass = totalMass + myMass
255 Eliminate.append(xyzIndex)
256 xyzIndex = xyzIndex + 1
257
258 mov_com = mov_com / totalMass
259
260 RBpos.append(mov_com)
261
262 for i in range(len(RBindices)):
263 mov[i] = mov[i] - mov_com
264
265 R = numpy.zeros((3, 3), numpy.float)
266 E0 = 0.0
267
268 for n in range(len(RBindices)):
269
270 # correlation matrix R:
271 # R(i,j) = sum(over n): y(n,i) * x(n,j)
272 # where x(n) and y(n) are two vector sets
273
274 R = R + numpy.outer(mov[n], ref_[n])
275
276 v, s, w = numpy.linalg.svd(R, full_matrices = True)
277
278 if (numpy.linalg.det(v) * numpy.linalg.det(w) < 0.0):
279 is_reflection = True
280 else:
281 is_reflection = False
282
283 if (is_reflection):
284 s[2] = -s[2]
285
286 RotMat = numpy.zeros((3, 3), numpy.float)
287 RotMat = v * w
288
289 q = numpy.array([0.0, 0.0, 0.0, 0.0], numpy.float)
290
291 t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
292
293 if( t > 1e-6 ):
294 s = 0.5 / math.sqrt( t )
295 q[0] = 0.25 / s
296 q[1] = (RotMat[1][2] - RotMat[2][1]) * s
297 q[2] = (RotMat[2][0] - RotMat[0][2]) * s
298 q[3] = (RotMat[0][1] - RotMat[1][0]) * s
299 else:
300 ad1 = RotMat[0][0]
301 ad2 = RotMat[1][1]
302 ad3 = RotMat[2][2]
303
304 if( ad1 >= ad2 and ad1 >= ad3 ):
305 s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
306 q[0] = (RotMat[1][2] - RotMat[2][1]) * s
307 q[1] = 0.25 / s
308 q[2] = (RotMat[0][1] + RotMat[1][0]) * s
309 q[3] = (RotMat[0][2] + RotMat[2][0]) * s
310 elif ( ad2 >= ad1 and ad2 >= ad3 ):
311 s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
312 q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
313 q[1] = (RotMat[0][1] + RotMat[1][0]) * s
314 q[2] = 0.25 / s
315 q[3] = (RotMat[1][2] + RotMat[2][1]) * s
316 else:
317 s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
318 q[0] = (RotMat[0][1] - RotMat[1][0]) * s
319 q[1] = (RotMat[0][2] + RotMat[2][0]) * s
320 q[2] = (RotMat[1][2] + RotMat[2][1]) * s
321 q[3] = 0.25 / s
322
323 RBQuats.append(q)
324
325
326
327def writeFile(outputFileName):
328 outputFile = open(outputFileName, 'w')
329
330 outputFile.write("<OpenMD version=1>\n");
331
332 for metaline in metaData:
333 outputFile.write(metaline)
334
335 outputFile.write(" <Snapshot>\n")
336
337 for frameline in frameData:
338 outputFile.write(frameline)
339
340 outputFile.write(" <StuntDoubles>\n")
341
342 sdFormat = 'pvqj'
343
344 index = 0
345 for i in range(len(RBpos)):
346 outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (index, sdFormat, RBpos[i][0], RBpos[i][1], RBpos[i][2], 0.0, 0.0, 0.0, RBQuats[i][0], RBQuats[i][1], RBQuats[i][2], RBQuats[i][3], 0.0, 0.0, 0.0))
347 index = index + 1
348
349
350 sdFormat = 'pv'
351 for i in range(len(indices)):
352 if i not in Eliminate:
353 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))
354
355 outputFile.write(" </StuntDoubles>\n")
356 outputFile.write(" </Snapshot>\n")
357 outputFile.write("</OpenMD>\n")
358 outputFile.close()
359
360def main(argv):
361 try:
362 opts, args = getopt.getopt(argv, "hr:x:o:", ["help", "rigid-body=", "xyz-file=", "output-file="])
363 except getopt.GetoptError:
364 usage()
365 sys.exit(2)
366 for opt, arg in opts:
367 if opt in ("-h", "--help"):
368 usage()
369 sys.exit()
370 elif opt in ("-x", "--xyz-file"):
371 xyzFileName = arg
372 global _haveXYZFileName
373 _haveXYZFileName = 1
374 elif opt in ("-r", "--rigid-body"):
375 rbFileName = arg
376 global _haveRBFileName
377 _haveRBFileName = 1
378 elif opt in ("-o", "--output-file"):
379 outputFileName = arg
380 global _haveOutputFileName
381 _haveOutputFileName = 1
382 if (_haveXYZFileName != 1):
383 usage()
384 print("No input packmol (xyz) file was specified")
385 sys.exit()
386 if (_haveRBFileName != 1):
387 usage()
388 print("No Rigid Body file (xyz) was specified")
389 sys.exit()
390 if (_haveOutputFileName != 1):
391 usage()
392 print("No output file was specified")
393 sys.exit()
394
395 readRBFile(rbFileName)
396 readFile(xyzFileName)
397 findRBs()
398 writeFile(outputFileName)
399
400if __name__ == "__main__":
401 if len(sys.argv) == 1:
402 usage()
403 sys.exit()
404 main(sys.argv[1:])