[ViewVC] Annotation of: OpenMD/branches/development/src/applications/utilities/md2md

#!/usr/bin/env python
"""MetaData file remapper

Takes a MetaData file and maps all StuntDoubles back to the periodic box.
Will optionally replicate the system in the three box directions, and
then writes out a new MetaData file.

Usage: md2md

Options:
  -h, --help              show this help
  -m, --meta-data=...     use specified meta-data (.md, .eor) file
  -o, --output-file=...   use specified output file
  -x, --repeatX=...       make the system repeat in the x direction
  -y, --repeatY=...       make the system repeat in the y direction
  -z, --repeatZ=...       make the system repeat in the z direction
  -d, --depthFirst        write out the replicated system depth-first


Example:
  md2md -m lipidSystem.md -o bigLipidSystem.md -x 2 -y 2 -z 1

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision: 1.3 $"
__date__ = "$Date: 2008-09-14 01:32:23 $"
__copyright__ = "Copyright (c) 2008 by the University of Notre Dame"
__license__ = "OOPSE"

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

_haveMDFileName = 0
_haveOutputFileName = 0

repeatX = 1
repeatY = 1
repeatZ = 1

metaData = []
frameData = []
indices = []
pvqj = []
p = []
v = []
q = []
j = []
Hmat = []
BoxInv = []

def usage():
    print __doc__

def readFile(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:
                metaData.append(line)
                line = mdFile.readline()
                if '</MetaData>' in line:
                    metaData.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:
                        frameData.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(','))
                            Hmat.append([Hxx, Hxy, Hxz])
                            Hmat.append([Hyx, Hyy, Hyz])
                            Hmat.append([Hzx, Hzy, Hzz])
                            BoxInv.append(1.0/Hxx)
                            BoxInv.append(1.0/Hyy)
                            BoxInv.append(1.0/Hzz)
                        line = mdFile.readline()
                        if '</FrameData>' in line:
                            frameData.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])
                        indices.append(myIndex)
                        pvqj.append(L[1])
                        x = float(L[2])
                        y = float(L[3])
                        z = float(L[4])
                        p.append([x, y, z])
                        vx = float(L[5])
                        vy = float(L[6])
                        vz = float(L[7])
                        v.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])
                            q.append([qw, qx, qy, qz])
                            jx = float(L[12])
                            jy = float(L[13])
                            jz = float(L[14])
                            j.append([jx, jy, jz])
                        else:
                            q.append([0.0, 0.0, 0.0, 0.0])
                            j.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,depthFirst):
    outputFile = open(outputFileName, 'w')

    outputFile.write("<OOPSE version=4>\n");

    print "writing MetaData"
    for metaline in metaData:
        if 'nMol' in metaline:
            metasplit = metaline.split()
            nMol = float(metasplit[2].strip(';'))
            newNmol = nMol * repeatX * repeatY * repeatZ
            outputFile.write('  nMol = %10d;\n' % (newNmol))
        else:
            outputFile.write(metaline)

    print "writing Snapshot"
    outputFile.write("  <Snapshot>\n")

    print "writing FrameData"
    for frameline in frameData:
        if 'Hmat:' in frameline:
            myH = []
            myH.append([repeatX * Hmat[0][0], repeatX * Hmat[0][1], repeatX * Hmat[0][2]])
            myH.append([repeatY * Hmat[1][0], repeatY * Hmat[1][1], repeatY * Hmat[1][2]])
            myH.append([repeatZ * Hmat[2][0], repeatZ * Hmat[2][1], repeatZ * Hmat[2][2]])
            outputFile.write("        Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n" % (myH[0][0], myH[0][1], myH[0][2], myH[1][0], myH[1][1], myH[1][2], myH[2][0], myH[2][1], myH[2][2]))
        else:
            outputFile.write(frameline)

    print "writing StuntDoubles"
    outputFile.write("    <StuntDoubles>\n")
    whichSD = 0

    print repeatX, repeatY, repeatZ


    if (depthFirst) :
        whichSD = 0
        for i in range(len(indices)):
            for ii in range(repeatX):
                for jj in range(repeatY):
                    for kk in range(repeatZ):
                        
                        myP = [] 
                        myP.append(p[i][0] + ii*Hmat[0][0] + jj*Hmat[1][0] + kk*Hmat[2][0])
                        myP.append(p[i][1] + ii*Hmat[0][1] + jj*Hmat[1][1] + kk*Hmat[2][1])
                        myP.append(p[i][2] + ii*Hmat[0][2] + jj*Hmat[1][2] + kk*Hmat[2][2])
                        
                        if (pvqj[i] == 'pv'):
                            outputFile.write("%10d %7s %18.10g %18.10g %18.10g %14e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], v[i][0], v[i][1], v[i][2]))
                        elif (pvqj[i] == 'pvqj'):
                            outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], v[i][0], v[i][1], v[i][2], q[i][0], q[i][1], q[i][2], q[i][3], j[i][0], j[i][1], j[i][2]))
                        
                        whichSD = whichSD + 1

    else:
        whichSD = 0
        for ii in range(repeatX):
            for jj in range(repeatY):
                for kk in range(repeatZ):
                    for i in range(len(indices)):
                        
                        myP = [] 
                        myP.append(p[i][0] + ii*Hmat[0][0] + jj*Hmat[1][0] + kk*Hmat[2][0])
                        myP.append(p[i][1] + ii*Hmat[0][1] + jj*Hmat[1][1] + kk*Hmat[2][1])
                        myP.append(p[i][2] + ii*Hmat[0][2] + jj*Hmat[1][2] + kk*Hmat[2][2])
                        
                        if (pvqj[i] == 'pv'):
                            outputFile.write("%10d %7s %18.10g %18.10g %18.10g %14e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], v[i][0], v[i][1], v[i][2]))
                        elif (pvqj[i] == 'pvqj'):
                            outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], v[i][0], v[i][1], v[i][2], q[i][0], q[i][1], q[i][2], q[i][3], j[i][0], j[i][1], j[i][2]))
                        
                        whichSD = whichSD + 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 wrapVector(myVect):
    scaled = [0.0, 0.0, 0.0]
    for i in range(3):
        scaled[i] = myVect[i] * BoxInv[i]
        scaled[i] = scaled[i] - roundMe(scaled[i])
        myVect[i] = scaled[i] * Hmat[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 mapToBox():
    for i in range(len(indices)):
        dpos = p[i]
        p[i] = wrapVector(dpos)


def main(argv):
    depthFirst = 0
    try:                                
        opts, args = getopt.getopt(argv, "hm:o:x:y:z:d", ["help", "meta-data=", "output-file=", "repeatX=", "repeatY=", "repeatZ=","depthFirst"]) 
    except getopt.GetoptError:           
        usage()                          
        sys.exit(2)                     
    for opt, arg in opts:                
        if opt in ("-h", "--help"):      
            usage()                     
            sys.exit()                  
        elif opt in ("-m", "--meta-data"): 
            mdFileName = arg
            global _haveMDFileName
            _haveMDFileName = 1
        elif opt in ("-o", "--output-file"): 
            outputFileName = arg
            global _haveOutputFileName
            _haveOutputFileName = 1
        elif opt in ("-x", "--repeatX"):
            global repeatX
            repeatX = int(arg)
        elif opt in ("-y", "--repeatY"): 
            global repeatY
            repeatY = int(arg)
        elif opt in ("-z", "--repeatZ"): 
            global repeatZ
            repeatZ = int(arg)
        elif opt in ("-d", "--depthFirst"):
            depthFirst = 1
    if (_haveMDFileName != 1):
        usage() 
        print "No meta-data file was specified"
        sys.exit()
    if (_haveOutputFileName != 1):
        usage()
        print "No output file was specified"
        sys.exit()

    readFile(mdFileName)
    mapToBox()    
    writeFile(outputFileName, depthFirst)

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1293
Committed:	Sun Sep 14 01:32:26 2008 UTC (16 years, 7 months ago) by chuckv
Original Path:	*trunk/src/applications/utilities/md2md*
File size:	11510 byte(s)
Log Message:	Added large quantities of code for convex hull and constant pressure langevin dynamics.
#	User	Rev	Content
1	gezelter	1227	#!/usr/bin/env python
2			"""MetaData file remapper
3
4			Takes a MetaData file and maps all StuntDoubles back to the periodic box.
5			Will optionally replicate the system in the three box directions, and
6			then writes out a new MetaData file.
7
8			Usage: md2md
9
10			Options:
11			-h, --help show this help
12			-m, --meta-data=... use specified meta-data (.md, .eor) file
13			-o, --output-file=... use specified output file
14			-x, --repeatX=... make the system repeat in the x direction
15			-y, --repeatY=... make the system repeat in the y direction
16			-z, --repeatZ=... make the system repeat in the z direction
17	chuckv	1293	-d, --depthFirst write out the replicated system depth-first
18	gezelter	1227
19
20			Example:
21			md2md -m lipidSystem.md -o bigLipidSystem.md -x 2 -y 2 -z 1
22
23			"""
24
25			__author__ = "Dan Gezelter (gezelter@nd.edu)"
26	chuckv	1293	__version__ = "$Revision: 1.3 $"
27			__date__ = "$Date: 2008-09-14 01:32:23 $"
28	gezelter	1227	__copyright__ = "Copyright (c) 2008 by the University of Notre Dame"
29			__license__ = "OOPSE"
30
31			import sys
32			import getopt
33			import string
34			import math
35			import random
36			from sets import *
37
38			_haveMDFileName = 0
39			_haveOutputFileName = 0
40
41			repeatX = 1
42			repeatY = 1
43			repeatZ = 1
44
45			metaData = []
46			frameData = []
47			indices = []
48			pvqj = []
49			p = []
50			v = []
51			q = []
52			j = []
53			Hmat = []
54			BoxInv = []
55
56			def usage():
57			print __doc__
58
59			def readFile(mdFileName):
60			mdFile = open(mdFileName, 'r')
61			# Find OOPSE version info first
62			line = mdFile.readline()
63			while 1:
64			if '<OOPSE version=' in line:
65			OOPSEversion = line
66			break
67			line = mdFile.readline()
68
69			# Rewind file and find start of MetaData block
70
71			mdFile.seek(0)
72			line = mdFile.readline()
73			print "reading MetaData"
74			while 1:
75			if '<MetaData>' in line:
76			while 2:
77			metaData.append(line)
78			line = mdFile.readline()
79			if '</MetaData>' in line:
80			metaData.append(line)
81			break
82			break
83			line = mdFile.readline()
84
85			mdFile.seek(0)
86			print "reading Snapshot"
87			line = mdFile.readline()
88			while 1:
89			if '<Snapshot>' in line:
90			line = mdFile.readline()
91			while 1:
92			print "reading FrameData"
93			if '<FrameData>' in line:
94			while 2:
95			frameData.append(line)
96			if 'Hmat:' in line:
97			L = line.split()
98			Hxx = float(L[2].strip(','))
99			Hxy = float(L[3].strip(','))
100			Hxz = float(L[4].strip(','))
101			Hyx = float(L[7].strip(','))
102			Hyy = float(L[8].strip(','))
103			Hyz = float(L[9].strip(','))
104			Hzx = float(L[12].strip(','))
105			Hzy = float(L[13].strip(','))
106			Hzz = float(L[14].strip(','))
107			Hmat.append([Hxx, Hxy, Hxz])
108			Hmat.append([Hyx, Hyy, Hyz])
109			Hmat.append([Hzx, Hzy, Hzz])
110			BoxInv.append(1.0/Hxx)
111			BoxInv.append(1.0/Hyy)
112			BoxInv.append(1.0/Hzz)
113			line = mdFile.readline()
114			if '</FrameData>' in line:
115			frameData.append(line)
116			break
117			break
118
119			line = mdFile.readline()
120			while 1:
121			if '<StuntDoubles>' in line:
122			line = mdFile.readline()
123			while 2:
124			L = line.split()
125			myIndex = int(L[0])
126			indices.append(myIndex)
127			pvqj.append(L[1])
128			x = float(L[2])
129			y = float(L[3])
130			z = float(L[4])
131			p.append([x, y, z])
132			vx = float(L[5])
133			vy = float(L[6])
134			vz = float(L[7])
135			v.append([vx, vy, vz])
136			if 'pvqj' in L[1]:
137			qw = float(L[8])
138			qx = float(L[9])
139			qy = float(L[10])
140			qz = float(L[11])
141			q.append([qw, qx, qy, qz])
142			jx = float(L[12])
143			jy = float(L[13])
144			jz = float(L[14])
145			j.append([jx, jy, jz])
146			else:
147			q.append([0.0, 0.0, 0.0, 0.0])
148			j.append([0.0, 0.0, 0.0])
149
150			line = mdFile.readline()
151			if '</StuntDoubles>' in line:
152			break
153			break
154			line = mdFile.readline()
155			if not line: break
156
157			mdFile.close()
158
159	chuckv	1293	def writeFile(outputFileName,depthFirst):
160	gezelter	1227	outputFile = open(outputFileName, 'w')
161
162			outputFile.write("<OOPSE version=4>\n");
163
164			print "writing MetaData"
165			for metaline in metaData:
166			if 'nMol' in metaline:
167			metasplit = metaline.split()
168			nMol = float(metasplit[2].strip(';'))
169			newNmol = nMol * repeatX * repeatY * repeatZ
170			outputFile.write(' nMol = %10d;\n' % (newNmol))
171			else:
172			outputFile.write(metaline)
173
174			print "writing Snapshot"
175			outputFile.write(" <Snapshot>\n")
176
177			print "writing FrameData"
178			for frameline in frameData:
179			if 'Hmat:' in frameline:
180			myH = []
181			myH.append([repeatX * Hmat[0][0], repeatX * Hmat[0][1], repeatX * Hmat[0][2]])
182			myH.append([repeatY * Hmat[1][0], repeatY * Hmat[1][1], repeatY * Hmat[1][2]])
183			myH.append([repeatZ * Hmat[2][0], repeatZ * Hmat[2][1], repeatZ * Hmat[2][2]])
184			outputFile.write(" Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n" % (myH[0][0], myH[0][1], myH[0][2], myH[1][0], myH[1][1], myH[1][2], myH[2][0], myH[2][1], myH[2][2]))
185			else:
186			outputFile.write(frameline)
187
188			print "writing StuntDoubles"
189			outputFile.write(" <StuntDoubles>\n")
190			whichSD = 0
191
192	xsun	1243	print repeatX, repeatY, repeatZ
193
194	chuckv	1293
195			if (depthFirst) :
196			whichSD = 0
197			for i in range(len(indices)):
198			for ii in range(repeatX):
199			for jj in range(repeatY):
200			for kk in range(repeatZ):
201
202			myP = []
203			myP.append(p[i][0] + iiHmat[0][0] + jjHmat[1][0] + kk*Hmat[2][0])
204			myP.append(p[i][1] + iiHmat[0][1] + jjHmat[1][1] + kk*Hmat[2][1])
205			myP.append(p[i][2] + iiHmat[0][2] + jjHmat[1][2] + kk*Hmat[2][2])
206
207			if (pvqj[i] == 'pv'):
208			outputFile.write("%10d %7s %18.10g %18.10g %18.10g %14e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], v[i][0], v[i][1], v[i][2]))
209			elif (pvqj[i] == 'pvqj'):
210			outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], v[i][0], v[i][1], v[i][2], q[i][0], q[i][1], q[i][2], q[i][3], j[i][0], j[i][1], j[i][2]))
211
212			whichSD = whichSD + 1
213
214			else:
215			whichSD = 0
216	gezelter	1227	for ii in range(repeatX):
217			for jj in range(repeatY):
218			for kk in range(repeatZ):
219	chuckv	1293	for i in range(len(indices)):
220	gezelter	1227
221	chuckv	1293	myP = []
222			myP.append(p[i][0] + iiHmat[0][0] + jjHmat[1][0] + kk*Hmat[2][0])
223			myP.append(p[i][1] + iiHmat[0][1] + jjHmat[1][1] + kk*Hmat[2][1])
224			myP.append(p[i][2] + iiHmat[0][2] + jjHmat[1][2] + kk*Hmat[2][2])
225
226			if (pvqj[i] == 'pv'):
227			outputFile.write("%10d %7s %18.10g %18.10g %18.10g %14e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], v[i][0], v[i][1], v[i][2]))
228			elif (pvqj[i] == 'pvqj'):
229			outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], v[i][0], v[i][1], v[i][2], q[i][0], q[i][1], q[i][2], q[i][3], j[i][0], j[i][1], j[i][2]))
230
231			whichSD = whichSD + 1
232	gezelter	1227
233	chuckv	1293
234	gezelter	1227	outputFile.write(" </StuntDoubles>\n")
235			outputFile.write(" </Snapshot>\n")
236			outputFile.write("</OOPSE>\n")
237			outputFile.close()
238
239			def roundMe(x):
240			if (x >= 0.0):
241			return math.floor(x + 0.5)
242			else:
243			return math.ceil(x - 0.5)
244
245			def wrapVector(myVect):
246			scaled = [0.0, 0.0, 0.0]
247			for i in range(3):
248			scaled[i] = myVect[i] * BoxInv[i]
249			scaled[i] = scaled[i] - roundMe(scaled[i])
250			myVect[i] = scaled[i] * Hmat[i][i]
251			return myVect
252
253			def dot(L1, L2):
254			myDot = 0.0
255			for i in range(len(L1)):
256			myDot = myDot + L1[i]*L2[i]
257			return myDot
258
259			def normalize(L1):
260			L2 = []
261			myLength = math.sqrt(dot(L1, L1))
262			for i in range(len(L1)):
263			L2.append(L1[i] / myLength)
264			return L2
265
266			def cross(L1, L2):
267			# don't call this with anything other than length 3 lists please
268			# or you'll be sorry
269			L3 = [0.0, 0.0, 0.0]
270			L3[0] = L1[1]L2[2] - L1[2]L2[1]
271			L3[1] = L1[2]L2[0] - L1[0]L2[2]
272			L3[2] = L1[0]L2[1] - L1[1]L2[0]
273			return L3
274
275			def mapToBox():
276			for i in range(len(indices)):
277			dpos = p[i]
278			p[i] = wrapVector(dpos)
279
280
281			def main(argv):
282	chuckv	1293	depthFirst = 0
283	gezelter	1227	try:
284	chuckv	1293	opts, args = getopt.getopt(argv, "hm:o:x:y:z:d", ["help", "meta-data=", "output-file=", "repeatX=", "repeatY=", "repeatZ=","depthFirst"])
285	gezelter	1227	except getopt.GetoptError:
286			usage()
287			sys.exit(2)
288			for opt, arg in opts:
289			if opt in ("-h", "--help"):
290			usage()
291			sys.exit()
292			elif opt in ("-m", "--meta-data"):
293			mdFileName = arg
294			global _haveMDFileName
295			_haveMDFileName = 1
296			elif opt in ("-o", "--output-file"):
297			outputFileName = arg
298			global _haveOutputFileName
299			_haveOutputFileName = 1
300	xsun	1243	elif opt in ("-x", "--repeatX"):
301			global repeatX
302	gezelter	1227	repeatX = int(arg)
303			elif opt in ("-y", "--repeatY"):
304	xsun	1243	global repeatY
305	gezelter	1227	repeatY = int(arg)
306			elif opt in ("-z", "--repeatZ"):
307	xsun	1243	global repeatZ
308	gezelter	1227	repeatZ = int(arg)
309	chuckv	1293	elif opt in ("-d", "--depthFirst"):
310			depthFirst = 1
311	gezelter	1227	if (_haveMDFileName != 1):
312			usage()
313			print "No meta-data file was specified"
314			sys.exit()
315			if (_haveOutputFileName != 1):
316			usage()
317			print "No output file was specified"
318			sys.exit()
319	chuckv	1293
320	gezelter	1227	readFile(mdFileName)
321			mapToBox()
322	chuckv	1293	writeFile(outputFileName, depthFirst)
323	gezelter	1227
324			if __name__ == "__main__":
325			if len(sys.argv) == 1:
326			usage()
327			sys.exit()
328			main(sys.argv[1:])