[ViewVC] Annotation of: OpenMD/trunk/src/applications/utilities/pack2md

#!@PYTHON_EXECUTABLE@
"""Packmol RigidBody Replacer

Finds atomistic rigid bodies in a packmol-generated xyz file and
generates a meta-data file with center of mass and orientational
coordinates for rigid bodies.

Usage: pack2md

Options:
  -h, --help              show this help
  -x,                     use the specified packmol (.xyz) file
  -r, --rigid-body=...    use this xyz structure as the rigid body
  -o, --output-file=...   use specified output (.md) file


Example:
   pack2md -x tolueneBox.xyz -r singleToluene.xyz -o tolueneBox.md

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision: 1646 $"
__date__ = "$Date: 2011-09-26 09:30:00 -0400 (Mon, 26 Sep 2011) $"
__copyright__ = "Copyright (c) 2011 by the University of Notre Dame"
__license__ = "OpenMD"

import sys
import getopt
import string
import math
import random
import numpy

_haveXYZFileName = 0
_haveRBFileName = 0
_haveOutputFileName = 0

atypes = []
positions = []
metaData = []
frameData = []
RBPos = []
RBQuats = []
indices = []
Hmat = []
BoxInv = []
H = []
Eliminate = []

mass_table = {
    'H': 1.00794,
    'C': 12.0107,
    'Cl': 35.453,
    'O': 15.999,
    'N': 14.007,
    'S': 32.0655,
    'Au': 196.466569,
    'CH4': 16.05,           
    'CH3': 15.04,           
    'CH2': 14.03,           
    'CH':  13.02,
    'CHar': 13.02,
    'CHa': 13.02,
    'RCHar': 12.0107,
    'RCH': 12.0107,    
    'CH3S': 15.04,
    'CH2S': 14.03,
    'CHS': 13.02,
    'CS': 12.0107,
    'SYZ': 32.0655,
    'SH': 32.0655,
    'HS': 1.0079,
    'S': 32.0655,
    'SZ': 32.0655,
    'SS': 32.0655, 
    'SP': 32.0655,
    'CS': 12.0107,
    'SAu': 228.9807,
    'SSD': 18.0153,         
    'SSD1': 18.0153,
    'SSD_E': 18.0153,
    'SSD_RF': 18.0153,                 
    'O_TIP3P': 15.9994,                 
    'O_TIP4P': 15.9994,
    'O_TIP4P-Ew': 15.9994,
    'O_TIP5P': 15.9994,
    'O_TIP5P-E': 15.9994, 
    'O_SPCE': 15.9994,
    'O_SPC': 15.9994,
    'H_TIP3P':    1.0079,
    'H_TIP4P':    1.0079,
    'H_TIP4P-Ew': 1.0079,
    'H_TIP5P':    1.0079,
    'H_SPCE':     1.0079,
    'H_SPC':      1.0079,
    'EP_TIP4P':   0.0,
    'EP_TIP4P-Ew':0.0,
    'EP_TIP5P':   0.0,
    'Ni':      58.710,
    'Cu':      63.550,
    'Rh':      102.90550,
    'Pd':      106.42,
    'Ag':      107.8682,
    'Ir':      192.217,
    'Pt':      195.09
}

#Hmat = zeros([3,3],Float)
#BoxInv = zeros([3],Float)

def usage():
    print __doc__

def readFile(XYZFileName):
    print "reading XYZ file"

    XYZFile = open(XYZFileName, 'r')        
    # Find number of atoms first
    line = XYZFile.readline()
    L = line.split()
    nAtoms = int(L[0])
    # skip comment line
    line = XYZFile.readline()
    for i in range(nAtoms):
        line = XYZFile.readline()
        L = line.split()
        myIndex = i
        indices.append(myIndex)
        atomType = L[0]
        atypes.append(atomType)
        x = float(L[1])
        y = float(L[2])
        z = float(L[3])
        positions.append([x, y, z])
    XYZFile.close()

def readRBFile(RBFileName):
    print "reading Rigid Body file"

    RBFile = open(RBFileName, 'r')        
    # Find number of atoms first
    line = RBFile.readline()
    L = line.split()
    nAtoms = int(L[0])
    # skip comment line
    line = RBFile.readline()
    for i in range(nAtoms):
        line = RBFile.readline()
        L = line.split()
        myIndex = i
        RBindices.append(myIndex)
        atomType = L[0]
        RBatypes.append(atomType)
        x = float(L[1])
        y = float(L[2])
        z = float(L[3])
        RBpositions.append([x, y, z])
    RBFile.close()

def findCOM():
    #find center of mass
    Xcom = 0.0
    Ycom = 0.0
    Zcom = 0.0
    totalMass = 0.0

    for i in range(0,len(RBindices)):
        myMass = mass_table[RBatypes[i]]
        
        Xcom = Xcom + myMass * RBpositions[i][0]
        Ycom = Ycom + myMass * RBpositions[i][1]
        Zcom = Zcom + myMass * RBpositions[i][2]
        totalMass = totalMass + myMass

    Xcom = Xcom / totalMass
    Ycom = Ycom / totalMass
    Zcom = Zcom / totalMass

    COM = [Xcom, Ycom, Zcom]

    return COM

def calcMoments():

    COM = findCOM()
    
    #find inertia tensor matrix elements

    I = numpy.zeros((3,3), numpy.float)
    
    for i in range(0,len(RBindices)):
        myMass = mass_table[RBatypes[i]]

        # move the origin of the reference coordinate system to the COM
        RBpositions[i][0] -= COM[0]
        RBpositions[i][1] -= COM[1]
        RBpositions[i][2] -= COM[2]

        dx = RBpositions[i][0]
        dy = RBpositions[i][1]
        dz = RBpositions[i][2]

        I[0,0] = I[0,0] + myMass * ( dy * dy + dz * dz )
        I[1,1] = I[1,1] + myMass * ( dx * dx + dz * dz )
        I[2,2] = I[2,2] + myMass * ( dx * dx + dy * dy )

        I[0,1] = I[0,1] - myMass * ( dx * dy )
        I[0,2] = I[0,2] - myMass * ( dx * dz )

        I[1,2] = I[1,2] - myMass * ( dy * dz )

        I[1,0] = I[0,1]
        I[2,0] = I[0,2]
        I[2,1] = I[1,2]

    print "Inertia Tensor:"
    print I
    print

    (evals, evects) = numpy.linalg.eig(I)
    print "evals:"
    print evals
    print 
    print "evects:"
    print evects
    print

    return (COM, evals, evects)

def findRBs():

    for i in range(len(RBindices)):
        ref_[i] = numpy.array([RBpositions[i][0], RBpositions[i][1], RBpositions[i][2]], numpy.float)

    print "finding rigid bodies (assuming strict packmol ordering)"
    xyzIndex =  0    

    for j in range(nBodies):
        mov_com = numpy.zeros(3, numpy.float)
        totalMass = 0.0
        
        for i in range(len(RBindices)):
            mov[i] = numpy.array([positions[xyzIndex][0], positions[xyzIndex][1], positions[xyzIndex][2]], numpy.float)
            myMass = mass_table[RBatypes[i]]
            mov_com = mov_com + myMass*mov[i]
            totalMass = totalMass + myMass
            xyzIndex = xyzIndex + 1

        mov_com = mov_com / totalMass

        RBpos.append(mov_com)

        for i in range(len(RBindices)):
            mov[i] = mov[i] - mov_com

        R = numpy.zeros((3,3), numpy.float)
        E0 = 0.0

        for n in range(len(RBindices)):
    
            # correlation matrix R:   
            #   R(i,j) = sum(over n): y(n,i) * x(n,j)  
            #   where x(n) and y(n) are two vector sets   

            R = R + numpy.linalg.outer(mov[n], ref_[n])

        v, s, w = numpy.linalg.svd(R, full_matrices = True)

        if (numpy.linalg.det(v) * numpy.linalg.det(w) < 0.0):
            is_reflection = true
        else
            is_reflection = false
            
        if (is_reflection):
            s[2] = -s[2]

        RotMat = numpy.zeros((3,3), numpy.float)
        RotMat = v * w
    
        q = numpy.array([0.0, 0.0, 0.0, 0.0], numpy.float)

        t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
        
        if( t > 1e-6 ):
            s = 0.5 / math.sqrt( t )
            q[0] = 0.25 / s
            q[1] = (RotMat[1][2] - RotMat[2][1]) * s
            q[2] = (RotMat[2][0] - RotMat[0][2]) * s
            q[3] = (RotMat[0][1] - RotMat[1][0]) * s
        else:
            ad1 = RotMat[0][0]
            ad2 = RotMat[1][1]
            ad3 = RotMat[2][2]

            if( ad1 >= ad2 and ad1 >= ad3 ):
                s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
                q[0] = (RotMat[1][2] - RotMat[2][1]) * s
                q[1] = 0.25 / s
                q[2] = (RotMat[0][1] + RotMat[1][0]) * s
                q[3] = (RotMat[0][2] + RotMat[2][0]) * s
            elif ( ad2 >= ad1 and ad2 >= ad3 ):
                s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
                q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
                q[1] = (RotMat[0][1] + RotMat[1][0]) * s
                q[2] = 0.25 / s
                q[3] = (RotMat[1][2] + RotMat[2][1]) * s
            else:
                s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
                q[0] = (RotMat[0][1] - RotMat[1][0]) * s
                q[1] = (RotMat[0][2] + RotMat[2][0]) * s
                q[2] = (RotMat[1][2] + RotMat[2][1]) * s
                q[3] = 0.25 / s

                RBQuats.append(q)


def writeFile(outputFileName):
    outputFile = open(outputFileName, 'w')

    outputFile.write("<OpenMD version=1>\n");
    
    for metaline in metaData:
        outputFile.write(metaline)

    outputFile.write("  <Snapshot>\n")

    for frameline in frameData:
        outputFile.write(frameline)
        
    outputFile.write("    <StuntDoubles>\n")

    sdFormat = 'pvqj'

    index = 0
    for i in range(len(RBPos)):
        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))
        index = index + 1


    sdFormat = 'pv'
    for i in range(len(indices)):
        if  i not in Eliminate:
            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))

    outputFile.write("    </StuntDoubles>\n")
    outputFile.write("  </Snapshot>\n")
    outputFile.write("</OpenMD>\n")
    outputFile.close()

def main(argv):                         
    try:                                
        opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz-file=", "output-file="]) 
    except getopt.GetoptError:           
        usage()                          
        sys.exit(2)                     
    for opt, arg in opts:                
        if opt in ("-h", "--help"):      
            usage()                     
            sys.exit()                  
        elif opt in ("-x", "--xyz-file"): 
            xyzFileName = arg
            global _haveXYZFileName
            _haveXYZFileName = 1
        elif opt in ("-r", "--rigid-body"): 
            rbFileName = arg
            global _haveRBFileName
            _haveRBFileName = 1
        elif opt in ("-o", "--output-file"): 
            outputFileName = arg
            global _haveOutputFileName
            _haveOutputFileName = 1
    if (_haveXYZFileName != 1):
        usage() 
        print "No input packmol (xyz) file was specified"
        sys.exit()
    if (_haveRBFileName != 1):
        usage() 
        print "No Rigid Body file (xyz) was specified"
        sys.exit()
    if (_haveOutputFileName != 1):
        usage()
        print "No output file was specified"
        sys.exit()
    readRBFile(rbFileName)
    
    readFile(xyzFileName)

    findRBs()
    writeFile(outputFileName)

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1888
Committed:	Tue Jun 18 17:52:37 2013 UTC (12 years, 4 months ago) by gezelter
File size:	10828 byte(s)
Log Message:	Added initial work on packmol->OpenMD converter, fixed a bug in RNEMDStats
#	User	Rev	Content
1	gezelter	1888	#!@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:])