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