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

#!/usr/bin/env python
"""Ice Cube Solvator

Opens two md files, one with water in an ice structure,
and one with water in a liquid phase.   Deletes any overlapping
liquid molecules and merges the two md files.

Usage: waterRotator

Options:
  -h,  --help              show this help
  -u,  --solute=...        use specified meta-data (.md) file as the solute
  -v,  --solvent=...       use specified meta-data (.md) file as the solvent
  -r,  --rcut=...          specify the cutoff radius for deleting solvent
  -o,  --output-file=...   use specified output (.md) file


Example:
   iceCubeSolvator -u frosty.md -v tepid.md -r 4.0 -o lukewarm.md

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision: 1.8 $"
__date__ = "$Date: 2008-06-27 15:47:06 $"
__copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
__license__ = "OOPSE"

import sys
import getopt
import string
import math
import random
from sets import *
#from Numeric import *

_haveMDFileName1 = 0
_haveMDFileName2 = 0
_haveRcut = 0
_haveOutputFileName = 0
_haveNSoluteAtoms = 1
_haveNSolventAtoms = 1

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

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

keepers = []

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

def usage():
    print __doc__

def readFile1(mdFileName):
    mdFile = open(mdFileName, 'r')        
    # Find OOPSE version info first
    line = mdFile.readline()
    while 1:
        if '<OOPSE version=' in line:
            OOPSEversion = line
            break
        line = mdFile.readline()
        
        # Rewind file and find start of MetaData block
        
    mdFile.seek(0)
    line = mdFile.readline()

    print "reading MetaData"
    while 1:
        if '<MetaData>' in line:
            while 2:
                metaData1.append(line)
                line = mdFile.readline()
                if '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 OOPSE version info first
    line = mdFile.readline()
    while 1:
        if '<OOPSE version=' in line:
            OOPSEversion = line
            break
        line = mdFile.readline()
        
        # Rewind file and find start of MetaData block
        
    mdFile.seek(0)
    line = mdFile.readline()
    print "reading MetaData"
    while 1:
        if '<MetaData>' in line:                            
            while 2:
                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("<OOPSE version=4>\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;" % (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("</OOPSE>\n")
    outputFile.close()
    
def roundMe(x):
    if (x >= 0.0):
        return math.floor(x + 0.5)
    else:
        return math.ceil(x - 0.5)

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


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

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

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

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

def removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms):

    rcut2 = rcut*rcut
    nextMol = 0

    
    for i in range(0,len(indices2),nSolventAtoms):
        keepThisMolecule = 1
        for atom1 in range (i, (i+nSolventAtoms)):
            iPos = positions2[atom1]
        
            for j in range(0,len(indices1),nSoluteAtoms):
                for atom2 in range (j, (j+nSoluteAtoms)):
                    jPos = positions1[j]
                    dpos = [jPos[0]-iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
                    dpos = wrapVector(dpos)
                    dist2 = dot(dpos,dpos)
                    
                    if (dist2 < rcut2):
                        keepThisMolecule = 0
 #                       break
        
        keepers.append(keepThisMolecule)

    global nSolvents
    myIndex = len(indices2) - 1
    for i in range(0,len(keepers)):
        
        if (keepers[i] == 1):
            nSolvents = nSolvents + 1
            for j in range (i, i+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):
        nSoluteAtoms = 1

    if (_haveNSolventAtoms != 1):
        nSolventAtoms = 1

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

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1266
Committed:	Fri Jun 27 15:47:06 2008 UTC (17 years, 1 month ago) by chuckv
Original Path:	*trunk/src/applications/utilities/md-solvator*
File size:	16511 byte(s)
Log Message:	Md solvator writes out md files now. Index problem.
#	User	Rev	Content
1	kfletch2	1262	#!/usr/bin/env python
2			"""Ice Cube Solvator
3	gezelter	1116
4	kfletch2	1262	Opens two md files, one with water in an ice structure,
5			and one with water in a liquid phase. Deletes any overlapping
6			liquid molecules and merges the two md files.
7	gezelter	1116
8	kfletch2	1262	Usage: waterRotator
9	gezelter	1116
10	kfletch2	1262	Options:
11			-h, --help show this help
12			-u, --solute=... use specified meta-data (.md) file as the solute
13			-v, --solvent=... use specified meta-data (.md) file as the solvent
14			-r, --rcut=... specify the cutoff radius for deleting solvent
15			-o, --output-file=... use specified output (.md) file
16	gezelter	1116
17
18	kfletch2	1262	Example:
19			iceCubeSolvator -u frosty.md -v tepid.md -r 4.0 -o lukewarm.md
20	gezelter	1116
21	kfletch2	1262	"""
22	gezelter	1116
23	kfletch2	1262	__author__ = "Dan Gezelter (gezelter@nd.edu)"
24	chuckv	1266	__version__ = "$Revision: 1.8 $"
25			__date__ = "$Date: 2008-06-27 15:47:06 $"
26	kfletch2	1262	__copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
27			__license__ = "OOPSE"
28	gezelter	1116
29	kfletch2	1262	import sys
30			import getopt
31			import string
32			import math
33			import random
34			from sets import *
35			#from Numeric import *
36	gezelter	1116
37	kfletch2	1262	_haveMDFileName1 = 0
38			_haveMDFileName2 = 0
39			_haveRcut = 0
40			_haveOutputFileName = 0
41	chuckv	1265	_haveNSoluteAtoms = 1
42			_haveNSolventAtoms = 1
43	gezelter	1116
44	kfletch2	1262	metaData1 = []
45			frameData1 = []
46			positions1 = []
47			velocities1 = []
48			quaternions1 = []
49			angVels1 = []
50			indices1 = []
51			Hmat1 = []
52			BoxInv1 = []
53	chuckv	1263	pvqj1 = []
54	gezelter	1116
55	kfletch2	1262	metaData2 = []
56			frameData2 = []
57			positions2 = []
58			velocities2 = []
59			quaternions2 = []
60			angVels2 = []
61			indices2 = []
62			Hmat2 = []
63			BoxInv2 = []
64	chuckv	1263	pvqj2 = []
65	gezelter	1116
66	kfletch2	1262	keepers = []
67	gezelter	1116
68	chuckv	1266	soluteTypeLine = str()
69			solventTypeLine = str()
70			soluteMolLine = str()
71			nSolvents = 0
72
73	kfletch2	1262	def usage():
74			print __doc__
75	gezelter	1116
76	kfletch2	1262	def readFile1(mdFileName):
77			mdFile = open(mdFileName, 'r')
78			# Find OOPSE version info first
79			line = mdFile.readline()
80			while 1:
81			if '<OOPSE version=' in line:
82			OOPSEversion = line
83			break
84			line = mdFile.readline()
85
86			# Rewind file and find start of MetaData block
87
88			mdFile.seek(0)
89			line = mdFile.readline()
90	chuckv	1266
91	kfletch2	1262	print "reading MetaData"
92			while 1:
93			if '<MetaData>' in line:
94			while 2:
95			metaData1.append(line)
96			line = mdFile.readline()
97	chuckv	1266	if 'type' in line:
98			global soluteTypeLine
99			soluteTypeLine = line
100			if 'nMol' in line:
101			global soluteMolLine
102			soluteMolLine = line
103	kfletch2	1262	if '</MetaData>' in line:
104			metaData1.append(line)
105			break
106			break
107			line = mdFile.readline()
108	gezelter	1116
109	kfletch2	1262	mdFile.seek(0)
110			print "reading Snapshot"
111			line = mdFile.readline()
112			while 1:
113			if '<Snapshot>' in line:
114			line = mdFile.readline()
115			while 1:
116			print "reading FrameData"
117			if '<FrameData>' in line:
118			while 2:
119			frameData1.append(line)
120			if 'Hmat:' in line:
121			L = line.split()
122			Hxx = float(L[2].strip(','))
123			Hxy = float(L[3].strip(','))
124			Hxz = float(L[4].strip(','))
125			Hyx = float(L[7].strip(','))
126			Hyy = float(L[8].strip(','))
127			Hyz = float(L[9].strip(','))
128			Hzx = float(L[12].strip(','))
129			Hzy = float(L[13].strip(','))
130			Hzz = float(L[14].strip(','))
131			Hmat1.append([Hxx, Hxy, Hxz])
132			Hmat1.append([Hyx, Hyy, Hyz])
133			Hmat1.append([Hzx, Hzy, Hzz])
134			print Hmat1
135			BoxInv1.append(1.0/Hxx)
136			BoxInv1.append(1.0/Hyy)
137			BoxInv1.append(1.0/Hzz)
138			print BoxInv1
139			line = mdFile.readline()
140			if '</FrameData>' in line:
141			frameData1.append(line)
142			break
143			break
144	gezelter	1116
145	kfletch2	1262	line = mdFile.readline()
146			while 1:
147			if '<StuntDoubles>' in line:
148			line = mdFile.readline()
149			while 2:
150			L = line.split()
151			myIndex = int(L[0])
152			indices1.append(myIndex)
153			pvqj1.append(L[1])
154			x = float(L[2])
155			y = float(L[3])
156			z = float(L[4])
157			positions1.append([x, y, z])
158			vx = float(L[5])
159			vy = float(L[6])
160			vz = float(L[7])
161			velocities1.append([vx, vy, vz])
162	chuckv	1263	if 'pvqj' in L[1]:
163			qw = float(L[8])
164			qx = float(L[9])
165			qy = float(L[10])
166			qz = float(L[11])
167			quaternions1.append([qw, qx, qy, qz])
168			jx = float(L[12])
169			jy = float(L[13])
170			jz = float(L[14])
171			angVels1.append([jx, jy, jz])
172			else:
173			quaternions1.append([0.0, 0.0, 0.0, 0.0])
174			angVels1.append([0.0, 0.0, 0.0])
175
176	kfletch2	1262	line = mdFile.readline()
177			if '</StuntDoubles>' in line:
178			break
179			break
180			line = mdFile.readline()
181			if not line: break
182
183			mdFile.close()
184	gezelter	1116
185	kfletch2	1262	def readFile2(mdFileName):
186			mdFile = open(mdFileName, 'r')
187			# Find OOPSE version info first
188			line = mdFile.readline()
189			while 1:
190			if '<OOPSE version=' in line:
191			OOPSEversion = line
192			break
193			line = mdFile.readline()
194
195			# Rewind file and find start of MetaData block
196
197			mdFile.seek(0)
198			line = mdFile.readline()
199			print "reading MetaData"
200			while 1:
201	chuckv	1266	if '<MetaData>' in line:
202	kfletch2	1262	while 2:
203	chuckv	1266	if 'type' in line:
204			global solventTypeLine
205			solventTypeLine = line
206	kfletch2	1262	metaData2.append(line)
207			line = mdFile.readline()
208			if '</MetaData>' in line:
209			metaData2.append(line)
210			break
211			break
212			line = mdFile.readline()
213	gezelter	1116
214	kfletch2	1262	mdFile.seek(0)
215			print "reading Snapshot"
216			line = mdFile.readline()
217			while 1:
218			if '<Snapshot>' in line:
219			line = mdFile.readline()
220			while 1:
221			print "reading FrameData"
222			if '<FrameData>' in line:
223			while 2:
224			frameData2.append(line)
225			if 'Hmat:' in line:
226			L = line.split()
227			Hxx = float(L[2].strip(','))
228			Hxy = float(L[3].strip(','))
229			Hxz = float(L[4].strip(','))
230			Hyx = float(L[7].strip(','))
231			Hyy = float(L[8].strip(','))
232			Hyz = float(L[9].strip(','))
233			Hzx = float(L[12].strip(','))
234			Hzy = float(L[13].strip(','))
235			Hzz = float(L[14].strip(','))
236			Hmat2.append([Hxx, Hxy, Hxz])
237			Hmat2.append([Hyx, Hyy, Hyz])
238			Hmat2.append([Hzx, Hzy, Hzz])
239			print Hmat2
240			BoxInv2.append(1.0/Hxx)
241			BoxInv2.append(1.0/Hyy)
242			BoxInv2.append(1.0/Hzz)
243			print BoxInv2
244			line = mdFile.readline()
245			if '</FrameData>' in line:
246			frameData2.append(line)
247			break
248			break
249	gezelter	1116
250	kfletch2	1262	line = mdFile.readline()
251			while 1:
252			if '<StuntDoubles>' in line:
253			line = mdFile.readline()
254			while 2:
255			L = line.split()
256			myIndex = int(L[0])
257			indices2.append(myIndex)
258			pvqj2.append(L[1])
259			x = float(L[2])
260			y = float(L[3])
261			z = float(L[4])
262			positions2.append([x, y, z])
263			vx = float(L[5])
264			vy = float(L[6])
265			vz = float(L[7])
266			velocities2.append([vx, vy, vz])
267	chuckv	1263	if 'pvqj' in L[1]:
268			qw = float(L[8])
269			qx = float(L[9])
270			qy = float(L[10])
271			qz = float(L[11])
272			quaternions2.append([qw, qx, qy, qz])
273			jx = float(L[12])
274			jy = float(L[13])
275			jz = float(L[14])
276			angVels2.append([jx, jy, jz])
277			else:
278			quaternions1.append([0.0, 0.0, 0.0, 0.0])
279			angVels1.append([0.0, 0.0, 0.0])
280
281	kfletch2	1262	line = mdFile.readline()
282			if '</StuntDoubles>' in line:
283			break
284			break
285			line = mdFile.readline()
286			if not line: break
287
288			mdFile.close()
289	xsun	1214
290	kfletch2	1262	def writeFile(outputFileName):
291			outputFile = open(outputFileName, 'w')
292	gezelter	1116
293	chuckv	1266	outputFile.write("<OOPSE version=4>\n")
294	kfletch2	1262
295	chuckv	1266	# for metaline in metaData1:
296			# outputFile.write(metaline)
297			outputFile.write(" <MetaData>\n")
298			outputFile.write("\n\n")
299			outputFile.write("#include \"thiols.md\"\n")
300			outputFile.write("#include \"metals.md\"\n")
301			outputFile.write("\n\n")
302			outputFile.write("component{\n")
303			outputFile.write(soluteTypeLine)
304			outputFile.write(soluteMolLine)
305			outputFile.write("}\n")
306	gezelter	1116
307	chuckv	1266	outputFile.write("component{\n")
308			outputFile.write(solventTypeLine)
309			outputFile.write("nMol = %d;" % (nSolvents))
310			outputFile.write("}\n")
311			outputFile.write("\n\n")
312			outputFile.write("forceField = \"MnM\";\n\n")
313			outputFile.write("\n\n")
314			outputFile.write(" </MetaData>\n")
315	kfletch2	1262	outputFile.write(" <Snapshot>\n")
316	gezelter	1116
317	kfletch2	1262	for frameline in frameData1:
318			outputFile.write(frameline)
319
320			outputFile.write(" <StuntDoubles>\n")
321	gezelter	1116
322	chuckv	1263
323	kfletch2	1262	newIndex = 0
324			for i in range(len(indices1)):
325	chuckv	1266	if (pvqj1[i] == 'pv'):
326	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]))
327	chuckv	1266	elif(pvqj1[i] == 'pvqj'):
328	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]))
329
330	kfletch2	1262	newIndex = newIndex + 1
331	gezelter	1116
332	kfletch2	1262	outputFile.write(" </StuntDoubles>\n")
333			outputFile.write(" </Snapshot>\n")
334			outputFile.write("</OOPSE>\n")
335			outputFile.close()
336
337			def roundMe(x):
338			if (x >= 0.0):
339			return math.floor(x + 0.5)
340			else:
341			return math.ceil(x - 0.5)
342	gezelter	1116
343	chuckv	1263	def frange(start,stop,step=1.0):
344			while start < stop:
345			yield start
346			start += step
347
348
349	kfletch2	1262	def wrapVector(myVect):
350			scaled = [0.0, 0.0, 0.0]
351			for i in range(3):
352			scaled[i] = myVect[i] * BoxInv1[i]
353			scaled[i] = scaled[i] - roundMe(scaled[i])
354			myVect[i] = scaled[i] * Hmat1[i][i]
355			return myVect
356	gezelter	1116
357	kfletch2	1262	def dot(L1, L2):
358			myDot = 0.0
359			for i in range(len(L1)):
360			myDot = myDot + L1[i]*L2[i]
361			return myDot
362	gezelter	1116
363	kfletch2	1262	def normalize(L1):
364			L2 = []
365			myLength = math.sqrt(dot(L1, L1))
366			for i in range(len(L1)):
367			L2.append(L1[i] / myLength)
368			return L2
369	gezelter	1116
370	kfletch2	1262	def cross(L1, L2):
371			# don't call this with anything other than length 3 lists please
372			# or you'll be sorry
373			L3 = [0.0, 0.0, 0.0]
374			L3[0] = L1[1]L2[2] - L1[2]L2[1]
375			L3[1] = L1[2]L2[0] - L1[0]L2[2]
376			L3[2] = L1[0]L2[1] - L1[1]L2[0]
377			return L3
378	gezelter	1116
379	chuckv	1265	def removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms):
380	gezelter	1116
381	kfletch2	1262	rcut2 = rcut*rcut
382	chuckv	1263	nextMol = 0
383
384	chuckv	1266
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	1263	iPos = positions2[atom1]
389	kfletch2	1262
390	chuckv	1263	for j in range(0,len(indices1),nSoluteAtoms):
391	chuckv	1264	for atom2 in range (j, (j+nSoluteAtoms)):
392	chuckv	1263	jPos = positions1[j]
393			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	1263	if (dist2 < rcut2):
398			keepThisMolecule = 0
399	chuckv	1266	# break
400
401	kfletch2	1262	keepers.append(keepThisMolecule)
402	gezelter	1116
403	chuckv	1266	global nSolvents
404			myIndex = len(indices2) - 1
405			for i in range(0,len(keepers)):
406
407	kfletch2	1262	if (keepers[i] == 1):
408	chuckv	1266	nSolvents = nSolvents + 1
409	chuckv	1263	for j in range (i, i+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	1263	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			nSoluteAtoms = 1
481
482			if (_haveNSolventAtoms != 1):
483			nSolventAtoms = 1
484
485	kfletch2	1262	readFile1(mdFileName1)
486			readFile2(mdFileName2)
487	chuckv	1263	removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms)
488	kfletch2	1262	writeFile(outputFileName)
489	gezelter	1116
490	kfletch2	1262	if __name__ == "__main__":
491			if len(sys.argv) == 1:
492			usage()
493			sys.exit()
494			main(sys.argv[1:])