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