[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.7 $"
__date__ = "$Date: 2008-06-26 15:50:22 $"
__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
_haveNSoluteAtoms = 1
_haveNSolventAtoms = 1

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

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

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])
                        if 'pvqj' in L[1]:
                            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])
                        else:
                            quaternions1.append([0.0, 0.0, 0.0, 0.0])
                            angVels1.append([0.0, 0.0, 0.0])

                        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])
                        if 'pvqj' in L[1]:
                            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])
                        else:
                            quaternions1.append([0.0, 0.0, 0.0, 0.0])
                            angVels1.append([0.0, 0.0, 0.0])

                        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")


    newIndex = 0
    for i in range(len(indices1)):
        if (pvqj[i] == 'pv'):
            outputFile.write("%10d %7s %18.10g %18.10g %18.10g %14e %13e %13e\n" % (newIndex, pvqj1[i], positions1[i][0], positions1[i][1], positions1[i][2], velocities1[i][0], velocities1[i][1], velocities1[i][2]))
        elif(pvqj[i] == 'pvqj'):
            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 frange(start,stop,step=1.0):
    while start < stop:
        yield start
        start += step


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,nSolventAtoms,nSoluteAtoms):

    rcut2 = rcut*rcut
    nextMol = 0


    for i in range(0,len(indices2),nSolventAtoms):
        keepThisMolecule = 1
        for atom1 in range (i, (i+nSolventAtoms)):
            iPos = positions2[atom1]
        
            for j in range(0,len(indices1),nSoluteAtoms):
                for atom2 in range (j, (j+nSoluteAtoms)):
                    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(0,len(indices2),nSolventAtoms):
        if (keepers[i] == 1):
            for j in range (i, i+nSolventAtoms):
                if (pvqj2[j] == 'pv'):
                    positions1.append(positions2[j])
                    velocities1.append(velocities2[j])
                    quaternions1.append(quaternions2[j])
                    angVels1.append(angVels2[j])
                    indices1.append(indices2[j])

def main(argv):                         
    try:                                
        opts, args = getopt.getopt(argv, "hu:v:n:p:r:o:", ["help", "solute=","solvent=","nSoluteAtoms=","nSolventAtoms=", "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 ("-n", "--nSoluteAtoms"): 
            nSoluteAtoms = int(arg)
            global _haveNSoluteAtoms
            _haveNSoluteAtoms = 1
        elif opt in ("-p", "--nSolventAtoms"): 
            nSolventAtoms = int(arg)
            global _haveNSolventAtoms
            _haveNSolventAtoms = 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

    if (_haveNSoluteAtoms != 1):
        nSoluteAtoms = 1

    if (_haveNSolventAtoms != 1):
        nSolventAtoms = 1

    readFile1(mdFileName1)
    readFile2(mdFileName2)
    removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms)
    writeFile(outputFileName)

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1265
Committed:	Thu Jun 26 15:50:22 2008 UTC (17 years ago) by chuckv
Original Path:	*trunk/src/applications/utilities/md-solvator*
File size:	14938 byte(s)
Log Message:	More progress...
#	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	chuckv	1265	__version__ = "$Revision: 1.7 $"
25			__date__ = "$Date: 2008-06-26 15:50:22 $"
26	kfletch2	1262	__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	chuckv	1265	_haveNSoluteAtoms = 1
42			_haveNSolventAtoms = 1
43	gezelter	1116
44	kfletch2	1262	metaData1 = []
45			frameData1 = []
46			positions1 = []
47			velocities1 = []
48			quaternions1 = []
49			angVels1 = []
50			indices1 = []
51			Hmat1 = []
52			BoxInv1 = []
53	chuckv	1263	pvqj1 = []
54	gezelter	1116
55	kfletch2	1262	metaData2 = []
56			frameData2 = []
57			positions2 = []
58			velocities2 = []
59			quaternions2 = []
60			angVels2 = []
61			indices2 = []
62			Hmat2 = []
63			BoxInv2 = []
64	chuckv	1263	pvqj2 = []
65	gezelter	1116
66	kfletch2	1262	keepers = []
67	gezelter	1116
68	kfletch2	1262	def usage():
69			print __doc__
70	gezelter	1116
71	kfletch2	1262	def readFile1(mdFileName):
72			mdFile = open(mdFileName, 'r')
73			# Find OOPSE version info first
74			line = mdFile.readline()
75			while 1:
76			if '<OOPSE version=' in line:
77			OOPSEversion = line
78			break
79			line = mdFile.readline()
80
81			# Rewind file and find start of MetaData block
82
83			mdFile.seek(0)
84			line = mdFile.readline()
85			print "reading MetaData"
86			while 1:
87			if '<MetaData>' in line:
88			while 2:
89			metaData1.append(line)
90			line = mdFile.readline()
91			if '</MetaData>' in line:
92			metaData1.append(line)
93			break
94			break
95			line = mdFile.readline()
96	gezelter	1116
97	kfletch2	1262	mdFile.seek(0)
98			print "reading Snapshot"
99			line = mdFile.readline()
100			while 1:
101			if '<Snapshot>' in line:
102			line = mdFile.readline()
103			while 1:
104			print "reading FrameData"
105			if '<FrameData>' in line:
106			while 2:
107			frameData1.append(line)
108			if 'Hmat:' in line:
109			L = line.split()
110			Hxx = float(L[2].strip(','))
111			Hxy = float(L[3].strip(','))
112			Hxz = float(L[4].strip(','))
113			Hyx = float(L[7].strip(','))
114			Hyy = float(L[8].strip(','))
115			Hyz = float(L[9].strip(','))
116			Hzx = float(L[12].strip(','))
117			Hzy = float(L[13].strip(','))
118			Hzz = float(L[14].strip(','))
119			Hmat1.append([Hxx, Hxy, Hxz])
120			Hmat1.append([Hyx, Hyy, Hyz])
121			Hmat1.append([Hzx, Hzy, Hzz])
122			print Hmat1
123			BoxInv1.append(1.0/Hxx)
124			BoxInv1.append(1.0/Hyy)
125			BoxInv1.append(1.0/Hzz)
126			print BoxInv1
127			line = mdFile.readline()
128			if '</FrameData>' in line:
129			frameData1.append(line)
130			break
131			break
132	gezelter	1116
133	kfletch2	1262	line = mdFile.readline()
134			while 1:
135			if '<StuntDoubles>' in line:
136			line = mdFile.readline()
137			while 2:
138			L = line.split()
139			myIndex = int(L[0])
140			indices1.append(myIndex)
141			pvqj1.append(L[1])
142			x = float(L[2])
143			y = float(L[3])
144			z = float(L[4])
145			positions1.append([x, y, z])
146			vx = float(L[5])
147			vy = float(L[6])
148			vz = float(L[7])
149			velocities1.append([vx, vy, vz])
150	chuckv	1263	if 'pvqj' in L[1]:
151			qw = float(L[8])
152			qx = float(L[9])
153			qy = float(L[10])
154			qz = float(L[11])
155			quaternions1.append([qw, qx, qy, qz])
156			jx = float(L[12])
157			jy = float(L[13])
158			jz = float(L[14])
159			angVels1.append([jx, jy, jz])
160			else:
161			quaternions1.append([0.0, 0.0, 0.0, 0.0])
162			angVels1.append([0.0, 0.0, 0.0])
163
164	kfletch2	1262	line = mdFile.readline()
165			if '</StuntDoubles>' in line:
166			break
167			break
168			line = mdFile.readline()
169			if not line: break
170
171			mdFile.close()
172	gezelter	1116
173	kfletch2	1262	def readFile2(mdFileName):
174			mdFile = open(mdFileName, 'r')
175			# Find OOPSE version info first
176			line = mdFile.readline()
177			while 1:
178			if '<OOPSE version=' in line:
179			OOPSEversion = line
180			break
181			line = mdFile.readline()
182
183			# Rewind file and find start of MetaData block
184
185			mdFile.seek(0)
186			line = mdFile.readline()
187			print "reading MetaData"
188			while 1:
189			if '<MetaData>' in line:
190			while 2:
191			metaData2.append(line)
192			line = mdFile.readline()
193			if '</MetaData>' in line:
194			metaData2.append(line)
195			break
196			break
197			line = mdFile.readline()
198	gezelter	1116
199	kfletch2	1262	mdFile.seek(0)
200			print "reading Snapshot"
201			line = mdFile.readline()
202			while 1:
203			if '<Snapshot>' in line:
204			line = mdFile.readline()
205			while 1:
206			print "reading FrameData"
207			if '<FrameData>' in line:
208			while 2:
209			frameData2.append(line)
210			if 'Hmat:' in line:
211			L = line.split()
212			Hxx = float(L[2].strip(','))
213			Hxy = float(L[3].strip(','))
214			Hxz = float(L[4].strip(','))
215			Hyx = float(L[7].strip(','))
216			Hyy = float(L[8].strip(','))
217			Hyz = float(L[9].strip(','))
218			Hzx = float(L[12].strip(','))
219			Hzy = float(L[13].strip(','))
220			Hzz = float(L[14].strip(','))
221			Hmat2.append([Hxx, Hxy, Hxz])
222			Hmat2.append([Hyx, Hyy, Hyz])
223			Hmat2.append([Hzx, Hzy, Hzz])
224			print Hmat2
225			BoxInv2.append(1.0/Hxx)
226			BoxInv2.append(1.0/Hyy)
227			BoxInv2.append(1.0/Hzz)
228			print BoxInv2
229			line = mdFile.readline()
230			if '</FrameData>' in line:
231			frameData2.append(line)
232			break
233			break
234	gezelter	1116
235	kfletch2	1262	line = mdFile.readline()
236			while 1:
237			if '<StuntDoubles>' in line:
238			line = mdFile.readline()
239			while 2:
240			L = line.split()
241			myIndex = int(L[0])
242			indices2.append(myIndex)
243			pvqj2.append(L[1])
244			x = float(L[2])
245			y = float(L[3])
246			z = float(L[4])
247			positions2.append([x, y, z])
248			vx = float(L[5])
249			vy = float(L[6])
250			vz = float(L[7])
251			velocities2.append([vx, vy, vz])
252	chuckv	1263	if 'pvqj' in L[1]:
253			qw = float(L[8])
254			qx = float(L[9])
255			qy = float(L[10])
256			qz = float(L[11])
257			quaternions2.append([qw, qx, qy, qz])
258			jx = float(L[12])
259			jy = float(L[13])
260			jz = float(L[14])
261			angVels2.append([jx, jy, jz])
262			else:
263			quaternions1.append([0.0, 0.0, 0.0, 0.0])
264			angVels1.append([0.0, 0.0, 0.0])
265
266	kfletch2	1262	line = mdFile.readline()
267			if '</StuntDoubles>' in line:
268			break
269			break
270			line = mdFile.readline()
271			if not line: break
272
273			mdFile.close()
274	xsun	1214
275	kfletch2	1262	def writeFile(outputFileName):
276			outputFile = open(outputFileName, 'w')
277	gezelter	1116
278	kfletch2	1262	outputFile.write("<OOPSE version=4>\n");
279
280			for metaline in metaData1:
281			outputFile.write(metaline)
282	gezelter	1116
283	kfletch2	1262	outputFile.write(" <Snapshot>\n")
284	gezelter	1116
285	kfletch2	1262	for frameline in frameData1:
286			outputFile.write(frameline)
287
288			outputFile.write(" <StuntDoubles>\n")
289	gezelter	1116
290	chuckv	1263
291	kfletch2	1262	newIndex = 0
292			for i in range(len(indices1)):
293	chuckv	1263	if (pvqj[i] == 'pv'):
294			outputFile.write("%10d %7s %18.10g %18.10g %18.10g %14e %13e %13e\n" % (newIndex, pvqj1[i], positions1[i][0], positions1[i][1], positions1[i][2], velocities1[i][0], velocities1[i][1], velocities1[i][2]))
295			elif(pvqj[i] == 'pvqj'):
296			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]))
297
298	kfletch2	1262	newIndex = newIndex + 1
299	gezelter	1116
300	kfletch2	1262	outputFile.write(" </StuntDoubles>\n")
301			outputFile.write(" </Snapshot>\n")
302			outputFile.write("</OOPSE>\n")
303			outputFile.close()
304
305			def roundMe(x):
306			if (x >= 0.0):
307			return math.floor(x + 0.5)
308			else:
309			return math.ceil(x - 0.5)
310	gezelter	1116
311	chuckv	1263	def frange(start,stop,step=1.0):
312			while start < stop:
313			yield start
314			start += step
315
316
317	kfletch2	1262	def wrapVector(myVect):
318			scaled = [0.0, 0.0, 0.0]
319			for i in range(3):
320			scaled[i] = myVect[i] * BoxInv1[i]
321			scaled[i] = scaled[i] - roundMe(scaled[i])
322			myVect[i] = scaled[i] * Hmat1[i][i]
323			return myVect
324	gezelter	1116
325	kfletch2	1262	def dot(L1, L2):
326			myDot = 0.0
327			for i in range(len(L1)):
328			myDot = myDot + L1[i]*L2[i]
329			return myDot
330	gezelter	1116
331	kfletch2	1262	def normalize(L1):
332			L2 = []
333			myLength = math.sqrt(dot(L1, L1))
334			for i in range(len(L1)):
335			L2.append(L1[i] / myLength)
336			return L2
337	gezelter	1116
338	kfletch2	1262	def cross(L1, L2):
339			# don't call this with anything other than length 3 lists please
340			# or you'll be sorry
341			L3 = [0.0, 0.0, 0.0]
342			L3[0] = L1[1]L2[2] - L1[2]L2[1]
343			L3[1] = L1[2]L2[0] - L1[0]L2[2]
344			L3[2] = L1[0]L2[1] - L1[1]L2[0]
345			return L3
346	gezelter	1116
347	chuckv	1265	def removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms):
348	gezelter	1116
349	kfletch2	1262	rcut2 = rcut*rcut
350	chuckv	1263	nextMol = 0
351
352
353			for i in range(0,len(indices2),nSolventAtoms):
354	kfletch2	1262	keepThisMolecule = 1
355	chuckv	1264	for atom1 in range (i, (i+nSolventAtoms)):
356	chuckv	1263	iPos = positions2[atom1]
357	kfletch2	1262
358	chuckv	1263	for j in range(0,len(indices1),nSoluteAtoms):
359	chuckv	1264	for atom2 in range (j, (j+nSoluteAtoms)):
360	chuckv	1263	jPos = positions1[j]
361			dpos = [jPos[0]-iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
362			dpos = wrapVector(dpos)
363			dist2 = dot(dpos,dpos)
364	gezelter	1116
365	chuckv	1263	if (dist2 < rcut2):
366			keepThisMolecule = 0
367			break
368
369	kfletch2	1262	keepers.append(keepThisMolecule)
370	gezelter	1116
371	chuckv	1263	for i in range(0,len(indices2),nSolventAtoms):
372	kfletch2	1262	if (keepers[i] == 1):
373	chuckv	1263	for j in range (i, i+nSolventAtoms):
374	chuckv	1265	if (pvqj2[j] == 'pv'):
375			positions1.append(positions2[j])
376			velocities1.append(velocities2[j])
377			quaternions1.append(quaternions2[j])
378			angVels1.append(angVels2[j])
379			indices1.append(indices2[j])
380	gezelter	1116
381	kfletch2	1262	def main(argv):
382			try:
383	chuckv	1265	opts, args = getopt.getopt(argv, "hu:v:n:p:r:o:", ["help", "solute=","solvent=","nSoluteAtoms=","nSolventAtoms=", "rcut=" "output-file="])
384	kfletch2	1262	except getopt.GetoptError:
385			usage()
386			sys.exit(2)
387			for opt, arg in opts:
388			if opt in ("-h", "--help"):
389			usage()
390			sys.exit()
391			elif opt in ("-u", "--solute"):
392			mdFileName1 = arg
393			global _haveMDFileName1
394			_haveMDFileName1 = 1
395			elif opt in ("-v", "--solvent"):
396			mdFileName2 = arg
397			global _haveMDFileName2
398			_haveMDFileName2 = 1
399	chuckv	1263	elif opt in ("-n", "--nSoluteAtoms"):
400	chuckv	1265	nSoluteAtoms = int(arg)
401	chuckv	1263	global _haveNSoluteAtoms
402			_haveNSoluteAtoms = 1
403			elif opt in ("-p", "--nSolventAtoms"):
404	chuckv	1265	nSolventAtoms = int(arg)
405	chuckv	1263	global _haveNSolventAtoms
406			_haveNSolventAtoms = 1
407	kfletch2	1262	elif opt in ("-r", "--rcut"):
408			rcut = float(arg)
409			global _haveRcut
410			_haveRcut = 1
411			elif opt in ("-o", "--output-file"):
412			outputFileName = arg
413			global _haveOutputFileName
414			_haveOutputFileName = 1
415	gezelter	1116
416	kfletch2	1262	if (_haveMDFileName1 != 1):
417			usage()
418			print "No meta-data file was specified for the solute"
419			sys.exit()
420	gezelter	1116
421	kfletch2	1262	if (_haveMDFileName2 != 1):
422			usage()
423			print "No meta-data file was specified for the solvent"
424			sys.exit()
425	gezelter	1116
426	kfletch2	1262	if (_haveOutputFileName != 1):
427			usage()
428			print "No output file was specified"
429			sys.exit()
430	gezelter	1116
431	kfletch2	1262	if (_haveRcut != 1):
432			print "No cutoff radius was specified, using 4 angstroms"
433			rcut = 4.0
434	gezelter	1116
435	chuckv	1263	if (_haveNSoluteAtoms != 1):
436			nSoluteAtoms = 1
437
438			if (_haveNSolventAtoms != 1):
439			nSolventAtoms = 1
440
441	kfletch2	1262	readFile1(mdFileName1)
442			readFile2(mdFileName2)
443	chuckv	1263	removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms)
444	kfletch2	1262	writeFile(outputFileName)
445	gezelter	1116
446	kfletch2	1262	if __name__ == "__main__":
447			if len(sys.argv) == 1:
448			usage()
449			sys.exit()
450			main(sys.argv[1:])