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

#!/usr/bin/env python
"""MD Solvator

Opens two md files, one with a solute structure and one with a solvent structure.
Deletes any solvent molecules that overlap with solute molecules and produces a md file.
The md file is bare and must be edited to run properly in OpenMD.

Usage: md-solvator

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
  -n,  --nSoluteAtoms=...  Number of atoms in solute molecule, default is 1 atom.
  -p,  --nSolventAtoms=... Number of atoms in solvent molecule, default is 1 atom. 

Example:
   md-solvator -u frosty.md -v tepid.md --nSoluteAtoms=3 --nSolventAtoms=3 -r 4.0 -o lukewarm.md

"""

__author__ = "Charles Vardeman (cvardema@nd.edu)"
__version__ = "$Revision$"
__date__ = "$Date$"
__copyright__ = "Copyright (c) 2008 by the University of Notre Dame"
__license__ = "OpenMD"

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 = 0
_haveNSolventAtoms = 0

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

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

keepers = []

soluteTypeLine = str()
solventTypeLine = str()
soluteMolLine = str()
nSolvents = 0

def usage():
    print __doc__

def readFile1(mdFileName):
    mdFile = open(mdFileName, 'r')        
    # Find OpenMD version info first
    line = mdFile.readline()
    while 1:
        if '<OpenMD version=' in line or '<OOPSE version=' in line:
            OpenMDversion = 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 'type' in line:
                    global soluteTypeLine
                    soluteTypeLine = line
                if 'nMol' in line:
                    global soluteMolLine
                    soluteMolLine = line
                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 OpenMD version info first
    line = mdFile.readline()
    while 1:
        if '<OpenMD version=' in line or '<OOPSE version=':
            OpenMDversion = 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:
                if 'type' in line:
                    global solventTypeLine
                    solventTypeLine = line
                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("<OpenMD version=1>\n")
    
#    for metaline in metaData1:
#        outputFile.write(metaline)
    outputFile.write("   <MetaData>\n")
    outputFile.write("\n\n")
    outputFile.write("component{\n")
    outputFile.write(soluteTypeLine)
    outputFile.write(soluteMolLine)
    outputFile.write("}\n")

    outputFile.write("component{\n")
    outputFile.write(solventTypeLine)
    outputFile.write("nMol = %d;\n" % (nSolvents))
    outputFile.write("}\n")
    outputFile.write("\n\n")
    outputFile.write("    </MetaData>\n")
    outputFile.write("  <Snapshot>\n")

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


    newIndex = 0
    for i in range(len(indices1)):
        if (pvqj1[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(pvqj1[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("</OpenMD>\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
    print len(indices1),len(indices2)
    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)):
                for atom2 in range (j, (j+nSoluteAtoms), nSoluteAtoms):
                    jPos = positions1[atom2]
                    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
            if (keepThisMolecule == 0):
                break

        keepers.append(keepThisMolecule)

       
    global nSolvents
    myIndex = len(indices2) - 1
    for i in range(0,len(keepers)):
        
        if (keepers[i] == 1):
            nSolvents = nSolvents + 1
            atomStartIndex = i * nSolventAtoms
            for j in range (atomStartIndex, (atomStartIndex+nSolventAtoms)):
                indices1.append(myIndex)
                pvqj1.append(pvqj2[j])
                if (pvqj2[j] == 'pv'):
                    positions1.append(positions2[j])
                    velocities1.append(velocities2[j])
                    quaternions1.append([0.0, 0.0, 0.0, 0.0])
                    angVels1.append([0.0, 0.0, 0.0])
                else:
                    positions1.append(positions2[j])
                    velocities1.append(velocities2[j])
                    quaternions1.append(quaternions2[j])
                    angVels1.append(angVels2[j])
 #                   indices1.append(indices2[j])
                myIndex = myIndex +1
    
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):
        print "Number of solute atoms was not specified. Using 1 atom."
        nSoluteAtoms = 1

    if (_haveNSolventAtoms != 1):
        print "Number of solute atoms was not specified. Using 1 atom."
        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:	1442
Committed:	Mon May 10 17:28:26 2010 UTC (15 years, 5 months ago) by gezelter
File size:	16922 byte(s)
Log Message:	Adding property set to svn entries
#	User	Rev	Content
1	kfletch2	1262	#!/usr/bin/env python
2	chuckv	1268	"""MD Solvator
3	gezelter	1116
4	chuckv	1268	Opens two md files, one with a solute structure and one with a solvent structure.
5			Deletes any solvent molecules that overlap with solute molecules and produces a md file.
6	gezelter	1390	The md file is bare and must be edited to run properly in OpenMD.
7	gezelter	1116
8	chuckv	1268	Usage: md-solvator
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	chuckv	1268	-n, --nSoluteAtoms=... Number of atoms in solute molecule, default is 1 atom.
17			-p, --nSolventAtoms=... Number of atoms in solvent molecule, default is 1 atom.
18	gezelter	1116
19	kfletch2	1262	Example:
20	chuckv	1268	md-solvator -u frosty.md -v tepid.md --nSoluteAtoms=3 --nSolventAtoms=3 -r 4.0 -o lukewarm.md
21	gezelter	1116
22	kfletch2	1262	"""
23	gezelter	1116
24	chuckv	1268	__author__ = "Charles Vardeman (cvardema@nd.edu)"
25	gezelter	1442	__version__ = "$Revision$"
26			__date__ = "$Date$"
27	chuckv	1268	__copyright__ = "Copyright (c) 2008 by the University of Notre Dame"
28	gezelter	1390	__license__ = "OpenMD"
29	gezelter	1116
30	kfletch2	1262	import sys
31			import getopt
32			import string
33			import math
34			import random
35			from sets import *
36			#from Numeric import *
37	gezelter	1116
38	kfletch2	1262	_haveMDFileName1 = 0
39			_haveMDFileName2 = 0
40			_haveRcut = 0
41			_haveOutputFileName = 0
42	chuckv	1268	_haveNSoluteAtoms = 0
43			_haveNSolventAtoms = 0
44	gezelter	1116
45	kfletch2	1262	metaData1 = []
46			frameData1 = []
47			positions1 = []
48			velocities1 = []
49			quaternions1 = []
50			angVels1 = []
51			indices1 = []
52			Hmat1 = []
53			BoxInv1 = []
54	chuckv	1263	pvqj1 = []
55	gezelter	1116
56	kfletch2	1262	metaData2 = []
57			frameData2 = []
58			positions2 = []
59			velocities2 = []
60			quaternions2 = []
61			angVels2 = []
62			indices2 = []
63			Hmat2 = []
64			BoxInv2 = []
65	chuckv	1263	pvqj2 = []
66	gezelter	1116
67	kfletch2	1262	keepers = []
68	gezelter	1116
69	chuckv	1266	soluteTypeLine = str()
70			solventTypeLine = str()
71			soluteMolLine = str()
72			nSolvents = 0
73
74	kfletch2	1262	def usage():
75			print __doc__
76	gezelter	1116
77	kfletch2	1262	def readFile1(mdFileName):
78			mdFile = open(mdFileName, 'r')
79	gezelter	1390	# Find OpenMD version info first
80	kfletch2	1262	line = mdFile.readline()
81			while 1:
82	gezelter	1390	if '<OpenMD version=' in line or '<OOPSE version=' in line:
83			OpenMDversion = line
84	kfletch2	1262	break
85			line = mdFile.readline()
86
87			# Rewind file and find start of MetaData block
88
89			mdFile.seek(0)
90			line = mdFile.readline()
91	chuckv	1266
92	kfletch2	1262	print "reading MetaData"
93			while 1:
94			if '<MetaData>' in line:
95			while 2:
96			metaData1.append(line)
97			line = mdFile.readline()
98	chuckv	1266	if 'type' in line:
99			global soluteTypeLine
100			soluteTypeLine = line
101			if 'nMol' in line:
102			global soluteMolLine
103			soluteMolLine = line
104	kfletch2	1262	if '</MetaData>' in line:
105			metaData1.append(line)
106			break
107			break
108			line = mdFile.readline()
109	gezelter	1116
110	kfletch2	1262	mdFile.seek(0)
111			print "reading Snapshot"
112			line = mdFile.readline()
113			while 1:
114			if '<Snapshot>' in line:
115			line = mdFile.readline()
116			while 1:
117			print "reading FrameData"
118			if '<FrameData>' in line:
119			while 2:
120			frameData1.append(line)
121			if 'Hmat:' in line:
122			L = line.split()
123			Hxx = float(L[2].strip(','))
124			Hxy = float(L[3].strip(','))
125			Hxz = float(L[4].strip(','))
126			Hyx = float(L[7].strip(','))
127			Hyy = float(L[8].strip(','))
128			Hyz = float(L[9].strip(','))
129			Hzx = float(L[12].strip(','))
130			Hzy = float(L[13].strip(','))
131			Hzz = float(L[14].strip(','))
132			Hmat1.append([Hxx, Hxy, Hxz])
133			Hmat1.append([Hyx, Hyy, Hyz])
134			Hmat1.append([Hzx, Hzy, Hzz])
135			print Hmat1
136			BoxInv1.append(1.0/Hxx)
137			BoxInv1.append(1.0/Hyy)
138			BoxInv1.append(1.0/Hzz)
139			print BoxInv1
140			line = mdFile.readline()
141			if '</FrameData>' in line:
142			frameData1.append(line)
143			break
144			break
145	gezelter	1116
146	kfletch2	1262	line = mdFile.readline()
147			while 1:
148			if '<StuntDoubles>' in line:
149			line = mdFile.readline()
150			while 2:
151			L = line.split()
152			myIndex = int(L[0])
153			indices1.append(myIndex)
154			pvqj1.append(L[1])
155			x = float(L[2])
156			y = float(L[3])
157			z = float(L[4])
158			positions1.append([x, y, z])
159			vx = float(L[5])
160			vy = float(L[6])
161			vz = float(L[7])
162			velocities1.append([vx, vy, vz])
163	chuckv	1263	if 'pvqj' in L[1]:
164			qw = float(L[8])
165			qx = float(L[9])
166			qy = float(L[10])
167			qz = float(L[11])
168			quaternions1.append([qw, qx, qy, qz])
169			jx = float(L[12])
170			jy = float(L[13])
171			jz = float(L[14])
172			angVels1.append([jx, jy, jz])
173			else:
174			quaternions1.append([0.0, 0.0, 0.0, 0.0])
175			angVels1.append([0.0, 0.0, 0.0])
176
177	kfletch2	1262	line = mdFile.readline()
178			if '</StuntDoubles>' in line:
179			break
180			break
181			line = mdFile.readline()
182			if not line: break
183
184			mdFile.close()
185	gezelter	1116
186	kfletch2	1262	def readFile2(mdFileName):
187			mdFile = open(mdFileName, 'r')
188	gezelter	1390	# Find OpenMD version info first
189	kfletch2	1262	line = mdFile.readline()
190			while 1:
191	gezelter	1390	if '<OpenMD version=' in line or '<OOPSE version=':
192			OpenMDversion = line
193	kfletch2	1262	break
194			line = mdFile.readline()
195
196			# Rewind file and find start of MetaData block
197
198			mdFile.seek(0)
199			line = mdFile.readline()
200			print "reading MetaData"
201			while 1:
202	chuckv	1266	if '<MetaData>' in line:
203	kfletch2	1262	while 2:
204	chuckv	1266	if 'type' in line:
205			global solventTypeLine
206			solventTypeLine = line
207	kfletch2	1262	metaData2.append(line)
208			line = mdFile.readline()
209			if '</MetaData>' in line:
210			metaData2.append(line)
211			break
212			break
213			line = mdFile.readline()
214	gezelter	1116
215	kfletch2	1262	mdFile.seek(0)
216			print "reading Snapshot"
217			line = mdFile.readline()
218			while 1:
219			if '<Snapshot>' in line:
220			line = mdFile.readline()
221			while 1:
222			print "reading FrameData"
223			if '<FrameData>' in line:
224			while 2:
225			frameData2.append(line)
226			if 'Hmat:' in line:
227			L = line.split()
228			Hxx = float(L[2].strip(','))
229			Hxy = float(L[3].strip(','))
230			Hxz = float(L[4].strip(','))
231			Hyx = float(L[7].strip(','))
232			Hyy = float(L[8].strip(','))
233			Hyz = float(L[9].strip(','))
234			Hzx = float(L[12].strip(','))
235			Hzy = float(L[13].strip(','))
236			Hzz = float(L[14].strip(','))
237			Hmat2.append([Hxx, Hxy, Hxz])
238			Hmat2.append([Hyx, Hyy, Hyz])
239			Hmat2.append([Hzx, Hzy, Hzz])
240			print Hmat2
241			BoxInv2.append(1.0/Hxx)
242			BoxInv2.append(1.0/Hyy)
243			BoxInv2.append(1.0/Hzz)
244			print BoxInv2
245			line = mdFile.readline()
246			if '</FrameData>' in line:
247			frameData2.append(line)
248			break
249			break
250	gezelter	1116
251	kfletch2	1262	line = mdFile.readline()
252			while 1:
253			if '<StuntDoubles>' in line:
254			line = mdFile.readline()
255			while 2:
256			L = line.split()
257			myIndex = int(L[0])
258			indices2.append(myIndex)
259			pvqj2.append(L[1])
260			x = float(L[2])
261			y = float(L[3])
262			z = float(L[4])
263			positions2.append([x, y, z])
264			vx = float(L[5])
265			vy = float(L[6])
266			vz = float(L[7])
267			velocities2.append([vx, vy, vz])
268	chuckv	1263	if 'pvqj' in L[1]:
269			qw = float(L[8])
270			qx = float(L[9])
271			qy = float(L[10])
272			qz = float(L[11])
273			quaternions2.append([qw, qx, qy, qz])
274			jx = float(L[12])
275			jy = float(L[13])
276			jz = float(L[14])
277			angVels2.append([jx, jy, jz])
278			else:
279			quaternions1.append([0.0, 0.0, 0.0, 0.0])
280			angVels1.append([0.0, 0.0, 0.0])
281
282	kfletch2	1262	line = mdFile.readline()
283			if '</StuntDoubles>' in line:
284			break
285			break
286			line = mdFile.readline()
287			if not line: break
288
289			mdFile.close()
290	xsun	1214
291	kfletch2	1262	def writeFile(outputFileName):
292			outputFile = open(outputFileName, 'w')
293	gezelter	1116
294	gezelter	1390	outputFile.write("<OpenMD version=1>\n")
295	kfletch2	1262
296	chuckv	1266	# for metaline in metaData1:
297			# outputFile.write(metaline)
298			outputFile.write(" <MetaData>\n")
299			outputFile.write("\n\n")
300			outputFile.write("component{\n")
301			outputFile.write(soluteTypeLine)
302			outputFile.write(soluteMolLine)
303			outputFile.write("}\n")
304	gezelter	1116
305	chuckv	1266	outputFile.write("component{\n")
306			outputFile.write(solventTypeLine)
307	chuckv	1268	outputFile.write("nMol = %d;\n" % (nSolvents))
308	chuckv	1266	outputFile.write("}\n")
309			outputFile.write("\n\n")
310			outputFile.write(" </MetaData>\n")
311	kfletch2	1262	outputFile.write(" <Snapshot>\n")
312	gezelter	1116
313	kfletch2	1262	for frameline in frameData1:
314			outputFile.write(frameline)
315
316			outputFile.write(" <StuntDoubles>\n")
317	gezelter	1116
318	chuckv	1263
319	kfletch2	1262	newIndex = 0
320			for i in range(len(indices1)):
321	chuckv	1266	if (pvqj1[i] == 'pv'):
322	chuckv	1263	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]))
323	chuckv	1266	elif(pvqj1[i] == 'pvqj'):
324	chuckv	1263	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]))
325
326	kfletch2	1262	newIndex = newIndex + 1
327	gezelter	1116
328	kfletch2	1262	outputFile.write(" </StuntDoubles>\n")
329			outputFile.write(" </Snapshot>\n")
330	gezelter	1390	outputFile.write("</OpenMD>\n")
331	kfletch2	1262	outputFile.close()
332
333			def roundMe(x):
334			if (x >= 0.0):
335			return math.floor(x + 0.5)
336			else:
337			return math.ceil(x - 0.5)
338	gezelter	1116
339	chuckv	1263	def frange(start,stop,step=1.0):
340			while start < stop:
341			yield start
342			start += step
343
344
345	kfletch2	1262	def wrapVector(myVect):
346			scaled = [0.0, 0.0, 0.0]
347			for i in range(3):
348			scaled[i] = myVect[i] * BoxInv1[i]
349			scaled[i] = scaled[i] - roundMe(scaled[i])
350			myVect[i] = scaled[i] * Hmat1[i][i]
351			return myVect
352	gezelter	1116
353	kfletch2	1262	def dot(L1, L2):
354			myDot = 0.0
355			for i in range(len(L1)):
356			myDot = myDot + L1[i]*L2[i]
357			return myDot
358	gezelter	1116
359	kfletch2	1262	def normalize(L1):
360			L2 = []
361			myLength = math.sqrt(dot(L1, L1))
362			for i in range(len(L1)):
363			L2.append(L1[i] / myLength)
364			return L2
365	gezelter	1116
366	kfletch2	1262	def cross(L1, L2):
367			# don't call this with anything other than length 3 lists please
368			# or you'll be sorry
369			L3 = [0.0, 0.0, 0.0]
370			L3[0] = L1[1]L2[2] - L1[2]L2[1]
371			L3[1] = L1[2]L2[0] - L1[0]L2[2]
372			L3[2] = L1[0]L2[1] - L1[1]L2[0]
373			return L3
374	gezelter	1116
375	chuckv	1265	def removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms):
376	gezelter	1116
377	kfletch2	1262	rcut2 = rcut*rcut
378	chuckv	1263	nextMol = 0
379	chuckv	1293	print len(indices1),len(indices2)
380	chuckv	1263	for i in range(0,len(indices2),nSolventAtoms):
381	kfletch2	1262	keepThisMolecule = 1
382	chuckv	1264	for atom1 in range (i, (i+nSolventAtoms)):
383	chuckv	1293
384			iPos = positions2[atom1]
385	chuckv	1272	for j in range(0,len(indices1)):
386			for atom2 in range (j, (j+nSoluteAtoms), nSoluteAtoms):
387	chuckv	1268	jPos = positions1[atom2]
388	chuckv	1263	dpos = [jPos[0]-iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
389			dpos = wrapVector(dpos)
390			dist2 = dot(dpos,dpos)
391	chuckv	1266
392	chuckv	1272
393	chuckv	1263	if (dist2 < rcut2):
394			keepThisMolecule = 0
395	chuckv	1272	break
396			if (keepThisMolecule == 0):
397			break
398	chuckv	1268
399	kfletch2	1262	keepers.append(keepThisMolecule)
400	gezelter	1116
401	chuckv	1272
402	chuckv	1266	global nSolvents
403			myIndex = len(indices2) - 1
404			for i in range(0,len(keepers)):
405
406	kfletch2	1262	if (keepers[i] == 1):
407	chuckv	1266	nSolvents = nSolvents + 1
408	chuckv	1272	atomStartIndex = i * nSolventAtoms
409			for j in range (atomStartIndex, (atomStartIndex+nSolventAtoms)):
410	chuckv	1266	indices1.append(myIndex)
411			pvqj1.append(pvqj2[j])
412	chuckv	1265	if (pvqj2[j] == 'pv'):
413			positions1.append(positions2[j])
414			velocities1.append(velocities2[j])
415	chuckv	1266	quaternions1.append([0.0, 0.0, 0.0, 0.0])
416			angVels1.append([0.0, 0.0, 0.0])
417			else:
418			positions1.append(positions2[j])
419			velocities1.append(velocities2[j])
420	chuckv	1265	quaternions1.append(quaternions2[j])
421			angVels1.append(angVels2[j])
422	chuckv	1266	# indices1.append(indices2[j])
423			myIndex = myIndex +1
424
425	kfletch2	1262	def main(argv):
426			try:
427	chuckv	1265	opts, args = getopt.getopt(argv, "hu:v:n:p:r:o:", ["help", "solute=","solvent=","nSoluteAtoms=","nSolventAtoms=", "rcut=" "output-file="])
428	kfletch2	1262	except getopt.GetoptError:
429			usage()
430			sys.exit(2)
431			for opt, arg in opts:
432			if opt in ("-h", "--help"):
433			usage()
434			sys.exit()
435			elif opt in ("-u", "--solute"):
436			mdFileName1 = arg
437			global _haveMDFileName1
438			_haveMDFileName1 = 1
439			elif opt in ("-v", "--solvent"):
440			mdFileName2 = arg
441			global _haveMDFileName2
442			_haveMDFileName2 = 1
443	chuckv	1268	elif opt in ("-n", "--nSoluteAtoms"):
444	chuckv	1265	nSoluteAtoms = int(arg)
445	chuckv	1263	global _haveNSoluteAtoms
446			_haveNSoluteAtoms = 1
447			elif opt in ("-p", "--nSolventAtoms"):
448	chuckv	1265	nSolventAtoms = int(arg)
449	chuckv	1263	global _haveNSolventAtoms
450			_haveNSolventAtoms = 1
451	kfletch2	1262	elif opt in ("-r", "--rcut"):
452			rcut = float(arg)
453			global _haveRcut
454			_haveRcut = 1
455			elif opt in ("-o", "--output-file"):
456			outputFileName = arg
457			global _haveOutputFileName
458			_haveOutputFileName = 1
459	gezelter	1116
460	kfletch2	1262	if (_haveMDFileName1 != 1):
461			usage()
462			print "No meta-data file was specified for the solute"
463			sys.exit()
464	gezelter	1116
465	kfletch2	1262	if (_haveMDFileName2 != 1):
466			usage()
467			print "No meta-data file was specified for the solvent"
468			sys.exit()
469	gezelter	1116
470	kfletch2	1262	if (_haveOutputFileName != 1):
471			usage()
472			print "No output file was specified"
473			sys.exit()
474	gezelter	1116
475	kfletch2	1262	if (_haveRcut != 1):
476			print "No cutoff radius was specified, using 4 angstroms"
477	chuckv	1266	rcut =4.0
478	gezelter	1116
479	chuckv	1263	if (_haveNSoluteAtoms != 1):
480	chuckv	1268	print "Number of solute atoms was not specified. Using 1 atom."
481	chuckv	1263	nSoluteAtoms = 1
482
483			if (_haveNSolventAtoms != 1):
484	chuckv	1268	print "Number of solute atoms was not specified. Using 1 atom."
485	chuckv	1263	nSolventAtoms = 1
486
487	kfletch2	1262	readFile1(mdFileName1)
488			readFile2(mdFileName2)
489	chuckv	1263	removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms)
490	kfletch2	1262	writeFile(outputFileName)
491	gezelter	1116
492	kfletch2	1262	if __name__ == "__main__":
493			if len(sys.argv) == 1:
494			usage()
495			sys.exit()
496			main(sys.argv[1:])