[ViewVC] Annotation of: OpenMD/branches/development/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: 1.12 $"
__date__ = "$Date: 2009-11-25 20:01:59 $"
__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("#include \"thiols.md\"\n")
    outputFile.write("#include \"metals.md\"\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("forceField = \"MnM\";\n\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:	1390
Committed:	Wed Nov 25 20:02:06 2009 UTC (15 years, 7 months ago) by gezelter
Original Path:	*trunk/src/applications/utilities/md-solvator*
File size:	17157 byte(s)
Log Message:	Almost all of the changes necessary to create OpenMD out of our old project (OOPSE-4)
#	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	1390	__version__ = "$Revision: 1.12 $"
26			__date__ = "$Date: 2009-11-25 20:01:59 $"
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("#include \"thiols.md\"\n")
301			outputFile.write("#include \"metals.md\"\n")
302			outputFile.write("\n\n")
303			outputFile.write("component{\n")
304			outputFile.write(soluteTypeLine)
305			outputFile.write(soluteMolLine)
306			outputFile.write("}\n")
307	gezelter	1116
308	chuckv	1266	outputFile.write("component{\n")
309			outputFile.write(solventTypeLine)
310	chuckv	1268	outputFile.write("nMol = %d;\n" % (nSolvents))
311	chuckv	1266	outputFile.write("}\n")
312			outputFile.write("\n\n")
313			outputFile.write("forceField = \"MnM\";\n\n")
314			outputFile.write("\n\n")
315			outputFile.write(" </MetaData>\n")
316	kfletch2	1262	outputFile.write(" <Snapshot>\n")
317	gezelter	1116
318	kfletch2	1262	for frameline in frameData1:
319			outputFile.write(frameline)
320
321			outputFile.write(" <StuntDoubles>\n")
322	gezelter	1116
323	chuckv	1263
324	kfletch2	1262	newIndex = 0
325			for i in range(len(indices1)):
326	chuckv	1266	if (pvqj1[i] == 'pv'):
327	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]))
328	chuckv	1266	elif(pvqj1[i] == 'pvqj'):
329	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]))
330
331	kfletch2	1262	newIndex = newIndex + 1
332	gezelter	1116
333	kfletch2	1262	outputFile.write(" </StuntDoubles>\n")
334			outputFile.write(" </Snapshot>\n")
335	gezelter	1390	outputFile.write("</OpenMD>\n")
336	kfletch2	1262	outputFile.close()
337
338			def roundMe(x):
339			if (x >= 0.0):
340			return math.floor(x + 0.5)
341			else:
342			return math.ceil(x - 0.5)
343	gezelter	1116
344	chuckv	1263	def frange(start,stop,step=1.0):
345			while start < stop:
346			yield start
347			start += step
348
349
350	kfletch2	1262	def wrapVector(myVect):
351			scaled = [0.0, 0.0, 0.0]
352			for i in range(3):
353			scaled[i] = myVect[i] * BoxInv1[i]
354			scaled[i] = scaled[i] - roundMe(scaled[i])
355			myVect[i] = scaled[i] * Hmat1[i][i]
356			return myVect
357	gezelter	1116
358	kfletch2	1262	def dot(L1, L2):
359			myDot = 0.0
360			for i in range(len(L1)):
361			myDot = myDot + L1[i]*L2[i]
362			return myDot
363	gezelter	1116
364	kfletch2	1262	def normalize(L1):
365			L2 = []
366			myLength = math.sqrt(dot(L1, L1))
367			for i in range(len(L1)):
368			L2.append(L1[i] / myLength)
369			return L2
370	gezelter	1116
371	kfletch2	1262	def cross(L1, L2):
372			# don't call this with anything other than length 3 lists please
373			# or you'll be sorry
374			L3 = [0.0, 0.0, 0.0]
375			L3[0] = L1[1]L2[2] - L1[2]L2[1]
376			L3[1] = L1[2]L2[0] - L1[0]L2[2]
377			L3[2] = L1[0]L2[1] - L1[1]L2[0]
378			return L3
379	gezelter	1116
380	chuckv	1265	def removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms):
381	gezelter	1116
382	kfletch2	1262	rcut2 = rcut*rcut
383	chuckv	1263	nextMol = 0
384	chuckv	1293	print len(indices1),len(indices2)
385	chuckv	1263	for i in range(0,len(indices2),nSolventAtoms):
386	kfletch2	1262	keepThisMolecule = 1
387	chuckv	1264	for atom1 in range (i, (i+nSolventAtoms)):
388	chuckv	1293
389			iPos = positions2[atom1]
390	chuckv	1272	for j in range(0,len(indices1)):
391			for atom2 in range (j, (j+nSoluteAtoms), nSoluteAtoms):
392	chuckv	1268	jPos = positions1[atom2]
393	chuckv	1263	dpos = [jPos[0]-iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
394			dpos = wrapVector(dpos)
395			dist2 = dot(dpos,dpos)
396	chuckv	1266
397	chuckv	1272
398	chuckv	1263	if (dist2 < rcut2):
399			keepThisMolecule = 0
400	chuckv	1272	break
401			if (keepThisMolecule == 0):
402			break
403	chuckv	1268
404	kfletch2	1262	keepers.append(keepThisMolecule)
405	gezelter	1116
406	chuckv	1272
407	chuckv	1266	global nSolvents
408			myIndex = len(indices2) - 1
409			for i in range(0,len(keepers)):
410
411	kfletch2	1262	if (keepers[i] == 1):
412	chuckv	1266	nSolvents = nSolvents + 1
413	chuckv	1272	atomStartIndex = i * nSolventAtoms
414			for j in range (atomStartIndex, (atomStartIndex+nSolventAtoms)):
415	chuckv	1266	indices1.append(myIndex)
416			pvqj1.append(pvqj2[j])
417	chuckv	1265	if (pvqj2[j] == 'pv'):
418			positions1.append(positions2[j])
419			velocities1.append(velocities2[j])
420	chuckv	1266	quaternions1.append([0.0, 0.0, 0.0, 0.0])
421			angVels1.append([0.0, 0.0, 0.0])
422			else:
423			positions1.append(positions2[j])
424			velocities1.append(velocities2[j])
425	chuckv	1265	quaternions1.append(quaternions2[j])
426			angVels1.append(angVels2[j])
427	chuckv	1266	# indices1.append(indices2[j])
428			myIndex = myIndex +1
429
430	kfletch2	1262	def main(argv):
431			try:
432	chuckv	1265	opts, args = getopt.getopt(argv, "hu:v:n:p:r:o:", ["help", "solute=","solvent=","nSoluteAtoms=","nSolventAtoms=", "rcut=" "output-file="])
433	kfletch2	1262	except getopt.GetoptError:
434			usage()
435			sys.exit(2)
436			for opt, arg in opts:
437			if opt in ("-h", "--help"):
438			usage()
439			sys.exit()
440			elif opt in ("-u", "--solute"):
441			mdFileName1 = arg
442			global _haveMDFileName1
443			_haveMDFileName1 = 1
444			elif opt in ("-v", "--solvent"):
445			mdFileName2 = arg
446			global _haveMDFileName2
447			_haveMDFileName2 = 1
448	chuckv	1268	elif opt in ("-n", "--nSoluteAtoms"):
449	chuckv	1265	nSoluteAtoms = int(arg)
450	chuckv	1263	global _haveNSoluteAtoms
451			_haveNSoluteAtoms = 1
452			elif opt in ("-p", "--nSolventAtoms"):
453	chuckv	1265	nSolventAtoms = int(arg)
454	chuckv	1263	global _haveNSolventAtoms
455			_haveNSolventAtoms = 1
456	kfletch2	1262	elif opt in ("-r", "--rcut"):
457			rcut = float(arg)
458			global _haveRcut
459			_haveRcut = 1
460			elif opt in ("-o", "--output-file"):
461			outputFileName = arg
462			global _haveOutputFileName
463			_haveOutputFileName = 1
464	gezelter	1116
465	kfletch2	1262	if (_haveMDFileName1 != 1):
466			usage()
467			print "No meta-data file was specified for the solute"
468			sys.exit()
469	gezelter	1116
470	kfletch2	1262	if (_haveMDFileName2 != 1):
471			usage()
472			print "No meta-data file was specified for the solvent"
473			sys.exit()
474	gezelter	1116
475	kfletch2	1262	if (_haveOutputFileName != 1):
476			usage()
477			print "No output file was specified"
478			sys.exit()
479	gezelter	1116
480	kfletch2	1262	if (_haveRcut != 1):
481			print "No cutoff radius was specified, using 4 angstroms"
482	chuckv	1266	rcut =4.0
483	gezelter	1116
484	chuckv	1263	if (_haveNSoluteAtoms != 1):
485	chuckv	1268	print "Number of solute atoms was not specified. Using 1 atom."
486	chuckv	1263	nSoluteAtoms = 1
487
488			if (_haveNSolventAtoms != 1):
489	chuckv	1268	print "Number of solute atoms was not specified. Using 1 atom."
490	chuckv	1263	nSolventAtoms = 1
491
492	kfletch2	1262	readFile1(mdFileName1)
493			readFile2(mdFileName2)
494	chuckv	1263	removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms)
495	kfletch2	1262	writeFile(outputFileName)
496	gezelter	1116
497	kfletch2	1262	if __name__ == "__main__":
498			if len(sys.argv) == 1:
499			usage()
500			sys.exit()
501			main(sys.argv[1:])