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

#!@PYTHON_EXECUTABLE@
"""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 new combined md file.  The md file must be
edited to run properly in OpenMD.  Note that the two boxes must have
identical box geometries (specified on the Hmat line).

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 solute.md -v solvent.md -n 3 -p 3 -r 4.0 -o combined.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 solute 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 solute Snapshot"
    line = mdFile.readline()
    while 1:
        if '<Snapshot>' in line:
            line = mdFile.readline()
            while 1:
                print "reading solute 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])
                            BoxInv1.append(1.0/Hxx)
                            BoxInv1.append(1.0/Hyy)
                            BoxInv1.append(1.0/Hzz)
                        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 solvent 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 solvent Snapshot"
    line = mdFile.readline()
    while 1:
        if '<Snapshot>' in line:
            line = mdFile.readline()
            while 1:
                print "reading solvent 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])
                            BoxInv2.append(1.0/Hxx)
                            BoxInv2.append(1.0/Hyy)
                            BoxInv2.append(1.0/Hzz)
                        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 checkBoxes():
    boxTolerance = 1.0e-3
    maxDiff = 0.0
    for i in range(3):
        for j in range(3):
            diff = math.fabs( Hmat1[i][j] - Hmat2[i][j])
            if (diff > maxDiff):
               maxDiff = diff
    if (maxDiff > boxTolerance):
       print "The solute and solvent boxes have different geometries:"
       print "                     Solute           |                   Solvent"
       print " -------------------------------------|------------------------------------"
       for i in range(3):
           print( "|  %10.4g %10.4g %10.4g   |  %10.4g %10.4g %10.4g  |" % (Hmat1[i][0], Hmat1[i][1], Hmat1[i][2], Hmat2[i][0], Hmat2[i][1], Hmat2[i][2]))

       print " -------------------------------------|------------------------------------"
       sys.exit()

    
def roundMe(x):
    if (x >= 0.0):
        return math.floor(x + 0.5)
    else:
        return math.ceil(x - 0.5)

def frange(start,stop,step=1.0):
    while start < stop:
        yield start
        start += step


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

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

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

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

def removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms):

    rcut2 = rcut*rcut
    nextMol = 0
    for i in range(0,len(indices2),nSolventAtoms):
        keepThisMolecule = 1
        for atom1 in range (i, (i+nSolventAtoms)):
           
            iPos = positions2[atom1]
            for j in range(0,len(indices1)):
                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)
    checkBoxes()
    removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms)
    writeFile(outputFileName)

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