[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, or 
translate every object in the box and before writing 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
  -t,  --translateX=...   translate all x coordinates by some amount 
  -u,  --translateY=...   translate all y coordinates by some amount
  -v,  --translateZ=...   translate all z coordinates by some amount

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

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision$"
__date__ = "$Date$"
__copyright__ = "Copyright (c) 2009 by the University of Notre Dame"
__license__ = "OpenMD"

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

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

def usage():
    print __doc__

def readFile(mdFileName):
    mdFile = open(mdFileName, 'r')        
    # Find OpenMD version info first
    line = mdFile.readline()
    while 1:
        if '<OOPSE version=' in line or '<OpenMD 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:
                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,repeatX,repeatY,repeatZ):
    outputFile = open(outputFileName, 'w')

    outputFile.write("<OpenMD version=1>\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")

    print repeatX, repeatY, repeatZ

    deco = [ (index, i) for i, index in enumerate(indices) ]
    deco.sort()
    whichSD = 0
    for ii in range(repeatX):
        for jj in range(repeatY):
            for kk in range(repeatZ):
                for index in range(len(deco)):
                    (index,i) = deco[index]
                    print i
                    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("</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 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(translateX, translateY, translateZ):
    print "translating by %f\t%f\t%f" % (translateX, translateY, translateZ)
    for i in range(len(indices)):
        dpos = []
        dpos.append(p[i][0] + translateX)
        dpos.append(p[i][1] + translateY)
        dpos.append(p[i][2] + translateZ)
        p[i] = wrapVector(dpos)

def main(argv):
    repeatX = 1
    repeatY = 1
    repeatZ = 1
    translateX = 0.0
    translateY = 0.0
    translateZ = 0.0
    _haveMDFileName = 0
    _haveOutputFileName = 0
    try:                                
        opts, args = getopt.getopt(argv, "hm:o:x:y:z:t:u:v:", ["help", "meta-data=", "output-file=", "repeatX=", "repeatY=", "repeatZ=","translateX=","translateY=","translateZ="]) 
    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
            _haveMDFileName = 1
        elif opt in ("-o", "--output-file"): 
            outputFileName = arg
            _haveOutputFileName = 1
        elif opt in ("-x", "--repeatX"):
            repeatX = int(arg)
        elif opt in ("-y", "--repeatY"): 
            repeatY = int(arg)
        elif opt in ("-z", "--repeatZ"): 
            repeatZ = int(arg)
        elif opt in ("-t", "--translateX"):
            translateX = float(arg)
        elif opt in ("-u", "--translateY"): 
            translateY = float(arg)
        elif opt in ("-v", "--translateZ"): 
            translateZ = float(arg)
    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(translateX, translateY, translateZ)    
    writeFile(outputFileName, repeatX, repeatY, repeatZ)

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