[ViewVC] Annotation of: OpenMD/branches/development/src/applications/utilities/md-solvator

#!/usr/bin/env python
"""Ice Cube Solvator

Opens two md files, one with water in an ice structure,
and one with water in a liquid phase.   Deletes any overlapping
liquid molecules and merges the two md files.

Usage: waterRotator

Options:
  -h,  --help              show this help
  -u,  --solute=...        use specified meta-data (.md) file as the solute
  -v,  --solvent=...       use specified meta-data (.md) file as the solvent
  -r,  --rcut=...          specify the cutoff radius for deleting solvent
  -o,  --output-file=...   use specified output (.md) file


Example:
   iceCubeSolvator -u frosty.md -v tepid.md -r 4.0 -o lukewarm.md

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision: 1.4 $"
__date__ = "$Date: 2008-06-25 18:22:08 $"
__copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
__license__ = "OOPSE"

import sys
import getopt
import string
import math
import random
from sets import *
#from Numeric import *

_haveMDFileName1 = 0
_haveMDFileName2 = 0
_haveRcut = 0
_haveOutputFileName = 0

metaData1 = []
frameData1 = []
positions1 = []
velocities1 = []
quaternions1 = []
angVels1 = []
indices1 = []
Hmat1 = []
BoxInv1 = []

metaData2 = []
frameData2 = []
positions2 = []
velocities2 = []
quaternions2 = []
angVels2 = []
indices2 = []
Hmat2 = []
BoxInv2 = []

keepers = []

def usage():
    print __doc__

def readFile1(mdFileName):
    mdFile = open(mdFileName, 'r')        
    # Find OOPSE version info first
    line = mdFile.readline()
    while 1:
        if '<OOPSE version=' in line:
            OOPSEversion = line
            break
        line = mdFile.readline()
        
        # Rewind file and find start of MetaData block
        
    mdFile.seek(0)
    line = mdFile.readline()
    print "reading MetaData"
    while 1:
        if '<MetaData>' in line:
            while 2:
                metaData1.append(line)
                line = mdFile.readline()
                if '</MetaData>' in line:
                    metaData1.append(line)
                    break
            break
        line = mdFile.readline()

    mdFile.seek(0)
    print "reading Snapshot"
    line = mdFile.readline()
    while 1:
        if '<Snapshot>' in line:
            line = mdFile.readline()
            while 1:
                print "reading FrameData"
                if '<FrameData>' in line:
                    while 2:
                        frameData1.append(line)
                        if 'Hmat:' in line:
                            L = line.split()
                            Hxx = float(L[2].strip(','))
                            Hxy = float(L[3].strip(','))
                            Hxz = float(L[4].strip(','))
                            Hyx = float(L[7].strip(','))
                            Hyy = float(L[8].strip(','))
                            Hyz = float(L[9].strip(','))
                            Hzx = float(L[12].strip(','))
                            Hzy = float(L[13].strip(','))
                            Hzz = float(L[14].strip(','))
                            Hmat1.append([Hxx, Hxy, Hxz])
                            Hmat1.append([Hyx, Hyy, Hyz])
                            Hmat1.append([Hzx, Hzy, Hzz])
                            print Hmat1
                            BoxInv1.append(1.0/Hxx)
                            BoxInv1.append(1.0/Hyy)
                            BoxInv1.append(1.0/Hzz)
                            print BoxInv1
                        line = mdFile.readline()
                        if '</FrameData>' in line:
                            frameData1.append(line)
                            break
                    break

            line = mdFile.readline()
            while 1:
                if '<StuntDoubles>' in line:
                    line = mdFile.readline()
                    while 2:
                        L = line.split()
                        myIndex = int(L[0])
                        indices1.append(myIndex)
                        pvqj1.append(L[1])
                        x = float(L[2])
                        y = float(L[3])
                        z = float(L[4])
                        positions1.append([x, y, z])
                        vx = float(L[5])
                        vy = float(L[6])
                        vz = float(L[7])
                        velocities1.append([vx, vy, vz])
                        qw = float(L[8])
                        qx = float(L[9])
                        qy = float(L[10])
                        qz = float(L[11])
                        quaternions1.append([qw, qx, qy, qz])
                        jx = float(L[12])
                        jy = float(L[13])
                        jz = float(L[14])
                        angVels1.append([jx, jy, jz])
                        line = mdFile.readline()
                        if '</StuntDoubles>' in line:
                            break
                    break
        line = mdFile.readline()
        if not line: break
    
    mdFile.close()

def readFile2(mdFileName):
    mdFile = open(mdFileName, 'r')        
    # Find OOPSE version info first
    line = mdFile.readline()
    while 1:
        if '<OOPSE version=' in line:
            OOPSEversion = line
            break
        line = mdFile.readline()
        
        # Rewind file and find start of MetaData block
        
    mdFile.seek(0)
    line = mdFile.readline()
    print "reading MetaData"
    while 1:
        if '<MetaData>' in line:
            while 2:
                metaData2.append(line)
                line = mdFile.readline()
                if '</MetaData>' in line:
                    metaData2.append(line)
                    break
            break
        line = mdFile.readline()

    mdFile.seek(0)
    print "reading Snapshot"
    line = mdFile.readline()
    while 1:
        if '<Snapshot>' in line:
            line = mdFile.readline()
            while 1:
                print "reading FrameData"
                if '<FrameData>' in line:
                    while 2:
                        frameData2.append(line)
                        if 'Hmat:' in line:
                            L = line.split()
                            Hxx = float(L[2].strip(','))
                            Hxy = float(L[3].strip(','))
                            Hxz = float(L[4].strip(','))
                            Hyx = float(L[7].strip(','))
                            Hyy = float(L[8].strip(','))
                            Hyz = float(L[9].strip(','))
                            Hzx = float(L[12].strip(','))
                            Hzy = float(L[13].strip(','))
                            Hzz = float(L[14].strip(','))
                            Hmat2.append([Hxx, Hxy, Hxz])
                            Hmat2.append([Hyx, Hyy, Hyz])
                            Hmat2.append([Hzx, Hzy, Hzz])
                            print Hmat2
                            BoxInv2.append(1.0/Hxx)
                            BoxInv2.append(1.0/Hyy)
                            BoxInv2.append(1.0/Hzz)
                            print BoxInv2
                        line = mdFile.readline()
                        if '</FrameData>' in line:
                            frameData2.append(line)
                            break
                    break

            line = mdFile.readline()
            while 1:
                if '<StuntDoubles>' in line:
                    line = mdFile.readline()
                    while 2:
                        L = line.split()
                        myIndex = int(L[0])
                        indices2.append(myIndex)
                        pvqj2.append(L[1])
                        x = float(L[2])
                        y = float(L[3])
                        z = float(L[4])
                        positions2.append([x, y, z])
                        vx = float(L[5])
                        vy = float(L[6])
                        vz = float(L[7])
                        velocities2.append([vx, vy, vz])
                        qw = float(L[8])
                        qx = float(L[9])
                        qy = float(L[10])
                        qz = float(L[11])
                        quaternions2.append([qw, qx, qy, qz])
                        jx = float(L[12])
                        jy = float(L[13])
                        jz = float(L[14])
                        angVels2.append([jx, jy, jz])
                        line = mdFile.readline()
                        if '</StuntDoubles>' in line:
                            break
                    break
        line = mdFile.readline()
        if not line: break
    
    mdFile.close()

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

    outputFile.write("<OOPSE version=4>\n");
    
    for metaline in metaData1:
        outputFile.write(metaline)

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

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

    sdFormat = 'pvqj'
    newIndex = 0
    for i in range(len(indices1)):
        
        outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (newIndex, pvqj1[i], positions1[i][0], positions1[i][1], positions1[i][2], velocities1[i][0], velocities1[i][1], velocities1[i][2], quaternions1[i][0], quaternions1[i][1], quaternions1[i][2], quaternions1[i][3], angVels1[i][0], angVels1[i][1], angVels1[i][2]))
        newIndex = newIndex + 1

    outputFile.write("    </StuntDoubles>\n")
    outputFile.write("  </Snapshot>\n")
    outputFile.write("</OOPSE>\n")
    outputFile.close()
    
def roundMe(x):
    if (x >= 0.0):
        return math.floor(x + 0.5)
    else:
        return math.ceil(x - 0.5)

def wrapVector(myVect):
    scaled = [0.0, 0.0, 0.0]
    for i in range(3):
        scaled[i] = myVect[i] * BoxInv1[i]
        scaled[i] = scaled[i] - roundMe(scaled[i])
        myVect[i] = scaled[i] * Hmat1[i][i]
    return myVect

def dot(L1, L2):
    myDot = 0.0
    for i in range(len(L1)):
        myDot = myDot + L1[i]*L2[i]
    return myDot

def normalize(L1):
    L2 = []
    myLength = math.sqrt(dot(L1, L1))
    for i in range(len(L1)):
        L2.append(L1[i] / myLength)
    return L2

def cross(L1, L2):
    # don't call this with anything other than length 3 lists please
    # or you'll be sorry
    L3 = [0.0, 0.0, 0.0]
    L3[0] = L1[1]*L2[2] - L1[2]*L2[1]
    L3[1] = L1[2]*L2[0] - L1[0]*L2[2]
    L3[2] = L1[0]*L2[1] - L1[1]*L2[0]
    return L3

def removeOverlaps(rcut):

    rcut2 = rcut*rcut
    for i in range(len(indices2)):
        keepThisMolecule = 1

        iPos = positions2[i]
        
        for j in range(len(indices1)):
            jPos = positions1[j]
            dpos = [jPos[0]-iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
            dpos = wrapVector(dpos)
            dist2 = dot(dpos,dpos)

            if (dist2 < rcut2):
                keepThisMolecule = 0
                break

        keepers.append(keepThisMolecule)

    for i in range(len(indices2)):
        if (keepers[i] == 1):
            positions1.append(positions2[i])
            velocities1.append(velocities2[i])
            quaternions1.append(quaternions2[i])
            angVels1.append(angVels2[i])
            indices1.append(indices2[i])

def main(argv):                         
    try:                                
        opts, args = getopt.getopt(argv, "hu:v:r:o:", ["help", "solute=","solvent=", "rcut" "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 ("-u", "--solute"): 
            mdFileName1 = arg
            global _haveMDFileName1
            _haveMDFileName1 = 1
        elif opt in ("-v", "--solvent"): 
            mdFileName2 = arg
            global _haveMDFileName2
            _haveMDFileName2 = 1
        elif opt in ("-r", "--rcut"): 
            rcut = float(arg)
            global _haveRcut
            _haveRcut = 1
        elif opt in ("-o", "--output-file"): 
            outputFileName = arg
            global _haveOutputFileName
            _haveOutputFileName = 1

    if (_haveMDFileName1 != 1):
        usage() 
        print "No meta-data file was specified for the solute"
        sys.exit()

    if (_haveMDFileName2 != 1):
        usage() 
        print "No meta-data file was specified for the solvent"
        sys.exit()

    if (_haveOutputFileName != 1):
        usage()
        print "No output file was specified"
        sys.exit()

    if (_haveRcut != 1):
        print "No cutoff radius was specified, using 4 angstroms"
        rcut = 4.0

    readFile1(mdFileName1)
    readFile2(mdFileName2)
    removeOverlaps(rcut)
    writeFile(outputFileName)

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1262
Committed:	Wed Jun 25 18:22:08 2008 UTC (17 years ago) by kfletch2
Original Path:	*trunk/src/applications/utilities/md-solvator*
File size:	13121 byte(s)
Log Message:	* empty log message *
#	User	Rev	Content
1	kfletch2	1262	#!/usr/bin/env python
2			"""Ice Cube Solvator
3	gezelter	1116
4	kfletch2	1262	Opens two md files, one with water in an ice structure,
5			and one with water in a liquid phase. Deletes any overlapping
6			liquid molecules and merges the two md files.
7	gezelter	1116
8	kfletch2	1262	Usage: waterRotator
9	gezelter	1116
10	kfletch2	1262	Options:
11			-h, --help show this help
12			-u, --solute=... use specified meta-data (.md) file as the solute
13			-v, --solvent=... use specified meta-data (.md) file as the solvent
14			-r, --rcut=... specify the cutoff radius for deleting solvent
15			-o, --output-file=... use specified output (.md) file
16	gezelter	1116
17
18	kfletch2	1262	Example:
19			iceCubeSolvator -u frosty.md -v tepid.md -r 4.0 -o lukewarm.md
20	gezelter	1116
21	kfletch2	1262	"""
22	gezelter	1116
23	kfletch2	1262	__author__ = "Dan Gezelter (gezelter@nd.edu)"
24			__version__ = "$Revision: 1.4 $"
25			__date__ = "$Date: 2008-06-25 18:22:08 $"
26			__copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
27			__license__ = "OOPSE"
28	gezelter	1116
29	kfletch2	1262	import sys
30			import getopt
31			import string
32			import math
33			import random
34			from sets import *
35			#from Numeric import *
36	gezelter	1116
37	kfletch2	1262	_haveMDFileName1 = 0
38			_haveMDFileName2 = 0
39			_haveRcut = 0
40			_haveOutputFileName = 0
41	gezelter	1116
42	kfletch2	1262	metaData1 = []
43			frameData1 = []
44			positions1 = []
45			velocities1 = []
46			quaternions1 = []
47			angVels1 = []
48			indices1 = []
49			Hmat1 = []
50			BoxInv1 = []
51	gezelter	1116
52	kfletch2	1262	metaData2 = []
53			frameData2 = []
54			positions2 = []
55			velocities2 = []
56			quaternions2 = []
57			angVels2 = []
58			indices2 = []
59			Hmat2 = []
60			BoxInv2 = []
61	gezelter	1116
62	kfletch2	1262	keepers = []
63	gezelter	1116
64	kfletch2	1262	def usage():
65			print __doc__
66	gezelter	1116
67	kfletch2	1262	def readFile1(mdFileName):
68			mdFile = open(mdFileName, 'r')
69			# Find OOPSE version info first
70			line = mdFile.readline()
71			while 1:
72			if '<OOPSE version=' in line:
73			OOPSEversion = line
74			break
75			line = mdFile.readline()
76
77			# Rewind file and find start of MetaData block
78
79			mdFile.seek(0)
80			line = mdFile.readline()
81			print "reading MetaData"
82			while 1:
83			if '<MetaData>' in line:
84			while 2:
85			metaData1.append(line)
86			line = mdFile.readline()
87			if '</MetaData>' in line:
88			metaData1.append(line)
89			break
90			break
91			line = mdFile.readline()
92	gezelter	1116
93	kfletch2	1262	mdFile.seek(0)
94			print "reading Snapshot"
95			line = mdFile.readline()
96			while 1:
97			if '<Snapshot>' in line:
98			line = mdFile.readline()
99			while 1:
100			print "reading FrameData"
101			if '<FrameData>' in line:
102			while 2:
103			frameData1.append(line)
104			if 'Hmat:' in line:
105			L = line.split()
106			Hxx = float(L[2].strip(','))
107			Hxy = float(L[3].strip(','))
108			Hxz = float(L[4].strip(','))
109			Hyx = float(L[7].strip(','))
110			Hyy = float(L[8].strip(','))
111			Hyz = float(L[9].strip(','))
112			Hzx = float(L[12].strip(','))
113			Hzy = float(L[13].strip(','))
114			Hzz = float(L[14].strip(','))
115			Hmat1.append([Hxx, Hxy, Hxz])
116			Hmat1.append([Hyx, Hyy, Hyz])
117			Hmat1.append([Hzx, Hzy, Hzz])
118			print Hmat1
119			BoxInv1.append(1.0/Hxx)
120			BoxInv1.append(1.0/Hyy)
121			BoxInv1.append(1.0/Hzz)
122			print BoxInv1
123			line = mdFile.readline()
124			if '</FrameData>' in line:
125			frameData1.append(line)
126			break
127			break
128	gezelter	1116
129	kfletch2	1262	line = mdFile.readline()
130			while 1:
131			if '<StuntDoubles>' in line:
132			line = mdFile.readline()
133			while 2:
134			L = line.split()
135			myIndex = int(L[0])
136			indices1.append(myIndex)
137			pvqj1.append(L[1])
138			x = float(L[2])
139			y = float(L[3])
140			z = float(L[4])
141			positions1.append([x, y, z])
142			vx = float(L[5])
143			vy = float(L[6])
144			vz = float(L[7])
145			velocities1.append([vx, vy, vz])
146			qw = float(L[8])
147			qx = float(L[9])
148			qy = float(L[10])
149			qz = float(L[11])
150			quaternions1.append([qw, qx, qy, qz])
151			jx = float(L[12])
152			jy = float(L[13])
153			jz = float(L[14])
154			angVels1.append([jx, jy, jz])
155			line = mdFile.readline()
156			if '</StuntDoubles>' in line:
157			break
158			break
159			line = mdFile.readline()
160			if not line: break
161
162			mdFile.close()
163	gezelter	1116
164	kfletch2	1262	def readFile2(mdFileName):
165			mdFile = open(mdFileName, 'r')
166			# Find OOPSE version info first
167			line = mdFile.readline()
168			while 1:
169			if '<OOPSE version=' in line:
170			OOPSEversion = line
171			break
172			line = mdFile.readline()
173
174			# Rewind file and find start of MetaData block
175
176			mdFile.seek(0)
177			line = mdFile.readline()
178			print "reading MetaData"
179			while 1:
180			if '<MetaData>' in line:
181			while 2:
182			metaData2.append(line)
183			line = mdFile.readline()
184			if '</MetaData>' in line:
185			metaData2.append(line)
186			break
187			break
188			line = mdFile.readline()
189	gezelter	1116
190	kfletch2	1262	mdFile.seek(0)
191			print "reading Snapshot"
192			line = mdFile.readline()
193			while 1:
194			if '<Snapshot>' in line:
195			line = mdFile.readline()
196			while 1:
197			print "reading FrameData"
198			if '<FrameData>' in line:
199			while 2:
200			frameData2.append(line)
201			if 'Hmat:' in line:
202			L = line.split()
203			Hxx = float(L[2].strip(','))
204			Hxy = float(L[3].strip(','))
205			Hxz = float(L[4].strip(','))
206			Hyx = float(L[7].strip(','))
207			Hyy = float(L[8].strip(','))
208			Hyz = float(L[9].strip(','))
209			Hzx = float(L[12].strip(','))
210			Hzy = float(L[13].strip(','))
211			Hzz = float(L[14].strip(','))
212			Hmat2.append([Hxx, Hxy, Hxz])
213			Hmat2.append([Hyx, Hyy, Hyz])
214			Hmat2.append([Hzx, Hzy, Hzz])
215			print Hmat2
216			BoxInv2.append(1.0/Hxx)
217			BoxInv2.append(1.0/Hyy)
218			BoxInv2.append(1.0/Hzz)
219			print BoxInv2
220			line = mdFile.readline()
221			if '</FrameData>' in line:
222			frameData2.append(line)
223			break
224			break
225	gezelter	1116
226	kfletch2	1262	line = mdFile.readline()
227			while 1:
228			if '<StuntDoubles>' in line:
229			line = mdFile.readline()
230			while 2:
231			L = line.split()
232			myIndex = int(L[0])
233			indices2.append(myIndex)
234			pvqj2.append(L[1])
235			x = float(L[2])
236			y = float(L[3])
237			z = float(L[4])
238			positions2.append([x, y, z])
239			vx = float(L[5])
240			vy = float(L[6])
241			vz = float(L[7])
242			velocities2.append([vx, vy, vz])
243			qw = float(L[8])
244			qx = float(L[9])
245			qy = float(L[10])
246			qz = float(L[11])
247			quaternions2.append([qw, qx, qy, qz])
248			jx = float(L[12])
249			jy = float(L[13])
250			jz = float(L[14])
251			angVels2.append([jx, jy, jz])
252			line = mdFile.readline()
253			if '</StuntDoubles>' in line:
254			break
255			break
256			line = mdFile.readline()
257			if not line: break
258
259			mdFile.close()
260	xsun	1214
261	kfletch2	1262	def writeFile(outputFileName):
262			outputFile = open(outputFileName, 'w')
263	gezelter	1116
264	kfletch2	1262	outputFile.write("<OOPSE version=4>\n");
265
266			for metaline in metaData1:
267			outputFile.write(metaline)
268	gezelter	1116
269	kfletch2	1262	outputFile.write(" <Snapshot>\n")
270	gezelter	1116
271	kfletch2	1262	for frameline in frameData1:
272			outputFile.write(frameline)
273
274			outputFile.write(" <StuntDoubles>\n")
275	gezelter	1116
276	kfletch2	1262	sdFormat = 'pvqj'
277			newIndex = 0
278			for i in range(len(indices1)):
279
280			outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (newIndex, pvqj1[i], positions1[i][0], positions1[i][1], positions1[i][2], velocities1[i][0], velocities1[i][1], velocities1[i][2], quaternions1[i][0], quaternions1[i][1], quaternions1[i][2], quaternions1[i][3], angVels1[i][0], angVels1[i][1], angVels1[i][2]))
281			newIndex = newIndex + 1
282	gezelter	1116
283	kfletch2	1262	outputFile.write(" </StuntDoubles>\n")
284			outputFile.write(" </Snapshot>\n")
285			outputFile.write("</OOPSE>\n")
286			outputFile.close()
287
288			def roundMe(x):
289			if (x >= 0.0):
290			return math.floor(x + 0.5)
291			else:
292			return math.ceil(x - 0.5)
293	gezelter	1116
294	kfletch2	1262	def wrapVector(myVect):
295			scaled = [0.0, 0.0, 0.0]
296			for i in range(3):
297			scaled[i] = myVect[i] * BoxInv1[i]
298			scaled[i] = scaled[i] - roundMe(scaled[i])
299			myVect[i] = scaled[i] * Hmat1[i][i]
300			return myVect
301	gezelter	1116
302	kfletch2	1262	def dot(L1, L2):
303			myDot = 0.0
304			for i in range(len(L1)):
305			myDot = myDot + L1[i]*L2[i]
306			return myDot
307	gezelter	1116
308	kfletch2	1262	def normalize(L1):
309			L2 = []
310			myLength = math.sqrt(dot(L1, L1))
311			for i in range(len(L1)):
312			L2.append(L1[i] / myLength)
313			return L2
314	gezelter	1116
315	kfletch2	1262	def cross(L1, L2):
316			# don't call this with anything other than length 3 lists please
317			# or you'll be sorry
318			L3 = [0.0, 0.0, 0.0]
319			L3[0] = L1[1]L2[2] - L1[2]L2[1]
320			L3[1] = L1[2]L2[0] - L1[0]L2[2]
321			L3[2] = L1[0]L2[1] - L1[1]L2[0]
322			return L3
323	gezelter	1116
324	kfletch2	1262	def removeOverlaps(rcut):
325	gezelter	1116
326	kfletch2	1262	rcut2 = rcut*rcut
327			for i in range(len(indices2)):
328			keepThisMolecule = 1
329	gezelter	1116
330	kfletch2	1262	iPos = positions2[i]
331
332			for j in range(len(indices1)):
333			jPos = positions1[j]
334			dpos = [jPos[0]-iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
335			dpos = wrapVector(dpos)
336			dist2 = dot(dpos,dpos)
337	gezelter	1116
338	kfletch2	1262	if (dist2 < rcut2):
339			keepThisMolecule = 0
340			break
341	gezelter	1116
342	kfletch2	1262	keepers.append(keepThisMolecule)
343	gezelter	1116
344	kfletch2	1262	for i in range(len(indices2)):
345			if (keepers[i] == 1):
346			positions1.append(positions2[i])
347			velocities1.append(velocities2[i])
348			quaternions1.append(quaternions2[i])
349			angVels1.append(angVels2[i])
350			indices1.append(indices2[i])
351	gezelter	1116
352	kfletch2	1262	def main(argv):
353			try:
354			opts, args = getopt.getopt(argv, "hu:v:r:o:", ["help", "solute=","solvent=", "rcut" "output-file="])
355			except getopt.GetoptError:
356			usage()
357			sys.exit(2)
358			for opt, arg in opts:
359			if opt in ("-h", "--help"):
360			usage()
361			sys.exit()
362			elif opt in ("-u", "--solute"):
363			mdFileName1 = arg
364			global _haveMDFileName1
365			_haveMDFileName1 = 1
366			elif opt in ("-v", "--solvent"):
367			mdFileName2 = arg
368			global _haveMDFileName2
369			_haveMDFileName2 = 1
370			elif opt in ("-r", "--rcut"):
371			rcut = float(arg)
372			global _haveRcut
373			_haveRcut = 1
374			elif opt in ("-o", "--output-file"):
375			outputFileName = arg
376			global _haveOutputFileName
377			_haveOutputFileName = 1
378	gezelter	1116
379	kfletch2	1262	if (_haveMDFileName1 != 1):
380			usage()
381			print "No meta-data file was specified for the solute"
382			sys.exit()
383	gezelter	1116
384	kfletch2	1262	if (_haveMDFileName2 != 1):
385			usage()
386			print "No meta-data file was specified for the solvent"
387			sys.exit()
388	gezelter	1116
389	kfletch2	1262	if (_haveOutputFileName != 1):
390			usage()
391			print "No output file was specified"
392			sys.exit()
393	gezelter	1116
394	kfletch2	1262	if (_haveRcut != 1):
395			print "No cutoff radius was specified, using 4 angstroms"
396			rcut = 4.0
397	gezelter	1116
398	kfletch2	1262	readFile1(mdFileName1)
399			readFile2(mdFileName2)
400			removeOverlaps(rcut)
401			writeFile(outputFileName)
402	gezelter	1116
403	kfletch2	1262	if __name__ == "__main__":
404			if len(sys.argv) == 1:
405			usage()
406			sys.exit()
407			main(sys.argv[1:])