[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


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

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision: 1.2 $"
__date__ = "$Date: 2008-05-14 20:49:17 $"
__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):
    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

    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

    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):
    try:                                
        opts, args = getopt.getopt(argv, "hm:o:x:y:z:", ["help", "meta-data=", "output-file=", "repeatX=", "repeatY=", "repeatZ="]) 
    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)
    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)

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