ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/applications/utilities/affineScale

Comparing trunk/src/applications/utilities/affineScale (file contents):
Revision 1063 by chrisfen, Tue Oct 10 14:08:10 2006 UTC vs.
Revision 1919 by jmarr, Wed Jul 31 18:03:35 2013 UTC

#	Line 1 | Line 1
1	<	#!/usr/bin/env perl
1	>	#!@PYTHON_EXECUTABLE@
2	>	"""MetaData affine scaling transform
3
4	<	# program that scales an OOPSE .in or .eor file to a new volume
4	>	Takes a MetaData file and scales both the periodic box and the
5	>	coordinates of all StuntDoubles in the system by the same amount.
6
7	<	# author    = "Chris Fennell
8	<	# version   = "$Revision: 1.1 $"
7	<	# date      = "$Date: 2006-10-10 14:08:10 $"
8	<	# copyright = "Copyright (c) 2006 by the University of Notre Dame"
9	<	# license   = "OOPSE"
7	>	You can either specify a new volume scaling for isotropic scaling,
8	>	or specify one (or more) of the coordinates for non-isotropic scaling.
9
10	+	Usage: affineScale
11
12	<	use Getopt::Std;
12	>	Options:
13	>	-h,  --help             show this help
14	>	-m,  --meta-data=...    use specified meta-data (.md, .eor) file
15	>	-o,  --output-file=...  use specified output file
16	>	-x,  --newX=...         scale the system to a new x dimension
17	>	-y,  --newY=...         scale the system to a new y dimension
18	>	-z,  --newZ=...         scale the system to a new z dimension
19	>	-v,  --newV=...         scale the system to a new volume
20
21	<	# some variables to get things going
22	<	$startSnap = 0;
16	<	$startFrame = 0;
17	<	$startStunts = 0;
21	>	Example:
22	>	affineScale -m lipidSystem.md -o scaledSystem.md -v 77000.0
23
24	<	# get our options
20	<	getopts('hv', \%opts);
24	>	"""
25
26	<	# if we don't have a filename, drop to -h
27	<	$opts{h} = 'true' if $#ARGV != 1;
26	>	__author__ = "Dan Gezelter (gezelter@nd.edu)"
27	>	__version__ = "$Revision$"
28	>	__date__ = "$Date$"
29	>	__copyright__ = "Copyright (c) 2009 by the University of Notre Dame"
30	>	__license__ = "OpenMD"
31
32	<	# our option output
33	<	if ($opts{h}){
34	<	print "affineScale: performs an affine transform on an OOPSE .in or .eor\n";
35	<	print "\tand writes to a new file named [file name].scale\n\n";
36	<	print "usage: affineScale [-hv] [file name] [new volume (Ang^3)]\n\n";
37	<	print "  -h : show this message\n";
31	<	die   "  -v : more verbose output\n";
32	<	}
32	>	import sys
33	>	import getopt
34	>	import string
35	>	import math
36	>	import random
37	>	from sets import *
38
39	<	# set some variables to be used in the code
40	<	$fileName = $ARGV[0];
41	<	$newVolume = $ARGV[1];
39	>	metaData = []
40	>	frameData = []
41	>	indices = []
42	>	pvqj = []
43	>	p = []
44	>	wrap = []
45	>	v = []
46	>	q = []
47	>	j = []
48	>	Hmat = []
49	>	BoxInv = []
50
51	<	# some crazy input checking
52	<	if ($newVolume =~ /^[0-9]/){
40	<	} else {
41	<	die "\t[new volume] value ($newVolume) is not a valid number\n\tPlease choose a non-negative numeric value for [new volume]\n";
42	<	}
51	>	def usage():
52	>	print __doc__
53
54	<	# split up the name to obtain an output filename
55	<	@names = split('\.', $fileName);
56	<	$names[$#names] = 'scale';
57	<	$outName = join('.',@names);
54	>	def readFile(mdFileName):
55	>	mdFile = open(mdFileName, 'r')
56	>	# Find OpenMD version info first
57	>	line = mdFile.readline()
58	>	while 1:
59	>	if '<OOPSE version=' in line or '<OpenMD version=' in line:
60	>	OpenMDversion = line
61	>	break
62	>	line = mdFile.readline()
63	>
64	>	# Rewind file and find start of MetaData block
65	>
66	>	mdFile.seek(0)
67	>	line = mdFile.readline()
68	>	print "reading MetaData"
69	>	while 1:
70	>	if '<MetaData>' in line:
71	>	while 2:
72	>	metaData.append(line)
73	>	line = mdFile.readline()
74	>	if '</MetaData>' in line:
75	>	metaData.append(line)
76	>	break
77	>	break
78	>	line = mdFile.readline()
79
80	<	open(STATFILE, "./$fileName") || die "\tError: can't find file $fileName\n";
81	<	open (SCALEFILE, ">./$outName") || die "\tError: can't open $outName";
80	>	mdFile.seek(0)
81	>	print "reading Snapshot"
82	>	line = mdFile.readline()
83	>	while 1:
84	>	if '<Snapshot>' in line:
85	>	line = mdFile.readline()
86	>	while 1:
87	>	print "reading FrameData"
88	>	if '<FrameData>' in line:
89	>	while 2:
90	>	frameData.append(line)
91	>	if 'Hmat:' in line:
92	>	L = line.split()
93	>	Hxx = float(L[2].strip(','))
94	>	Hxy = float(L[3].strip(','))
95	>	Hxz = float(L[4].strip(','))
96	>	Hyx = float(L[7].strip(','))
97	>	Hyy = float(L[8].strip(','))
98	>	Hyz = float(L[9].strip(','))
99	>	Hzx = float(L[12].strip(','))
100	>	Hzy = float(L[13].strip(','))
101	>	Hzz = float(L[14].strip(','))
102	>	Hmat.append([Hxx, Hxy, Hxz])
103	>	Hmat.append([Hyx, Hyy, Hyz])
104	>	Hmat.append([Hzx, Hzy, Hzz])
105	>	BoxInv.append(1.0/Hxx)
106	>	BoxInv.append(1.0/Hyy)
107	>	BoxInv.append(1.0/Hzz)
108	>	line = mdFile.readline()
109	>	if '</FrameData>' in line:
110	>	frameData.append(line)
111	>	break
112	>	break
113
114	<	while (<STATFILE>){
115	<	$startSnap = 0 if /\/Snapshot/;
116	<	$startFrame = 0 if /\/FrameData/;
117	<	$startStunts = 0 if /\/StuntDoubles/;
114	>	line = mdFile.readline()
115	>	while 1:
116	>	if '<StuntDoubles>' in line:
117	>	line = mdFile.readline()
118	>	while 2:
119	>	L = line.split()
120	>	myIndex = int(L[0])
121	>	indices.append(myIndex)
122	>	pvqj.append(L[1])
123	>	x = float(L[2])
124	>	y = float(L[3])
125	>	z = float(L[4])
126	>	p.append([x, y, z])
127	>	wrap.append([0,0,0])
128	>	vx = float(L[5])
129	>	vy = float(L[6])
130	>	vz = float(L[7])
131	>	v.append([vx, vy, vz])
132	>	if 'pvqj' in L[1]:
133	>	qw = float(L[8])
134	>	qx = float(L[9])
135	>	qy = float(L[10])
136	>	qz = float(L[11])
137	>	q.append([qw, qx, qy, qz])
138	>	jx = float(L[12])
139	>	jy = float(L[13])
140	>	jz = float(L[14])
141	>	j.append([jx, jy, jz])
142	>	else:
143	>	q.append([0.0, 0.0, 0.0, 0.0])
144	>	j.append([0.0, 0.0, 0.0])
145	>
146	>	line = mdFile.readline()
147	>	if '</StuntDoubles>' in line:
148	>	break
149	>	break
150	>	line = mdFile.readline()
151	>	if not line: break
152	>
153	>	mdFile.close()
154
155	<	if ($startSnap == 1){
156	<	if ($startFrame == 1){
59	<	if (/Hmat/){
60	<	@line = split;
61	<
62	<	chop $line[2];
63	<	chop $line[8];
64	<	$oldVolume = $line[2]*$line[8]*$line[14];
65	<
66	<	$scale = ($newVolume/$oldVolume)**0.333333333333333;
67	<	# scale the hmat vectors (only orthorhombic)
68	<	$line[2] *= $scale;
69	<	$line[8] *= $scale;
70	<	$line[14] *= $scale;
71	<	$volume = $line[2]*$line[8]*$line[14];
72	<
73	<	print SCALEFILE "\t$line[0]";
74	<	for ($i=1; $i<=$#line; $i++) {
75	<	print SCALEFILE " $line[$i]";
76	<	if ($i == 2 || $i == 8) {print SCALEFILE ",";}
77	<	}
78	<	print SCALEFILE "\n";
79	<	} else {
80	<	print SCALEFILE;
81	<	}
82	<	}
83	<	if ($startStunts == 1){
84	<	@line = split;
85	<	$line[2] *= $scale;
86	<	$line[3] *= $scale;
87	<	$line[4] *= $scale;
88	<	print SCALEFILE "$line[0]\t$line[1]\t";
89	<	for ($i=2; $i<=$#line; $i++) {
90	<	print SCALEFILE "$line[$i] ";
91	<	}
92	<	print SCALEFILE "\n";
93	<	}
94	<	if ($startFrame == 0 && $startStunts == 0) {
95	<	print SCALEFILE;
96	<	}
97	<	} else {
98	<	print SCALEFILE;
99	<	}
100	<	$startSnap = 1 if /Snapshot/;
101	<	$startFrame = 1 if /FrameData/;
102	<	$startStunts = 1 if /StuntDoubles/;
103	<	# check again since "/value" contains "value"
104	<	$startSnap = 0 if /\/Snapshot/;
105	<	$startFrame = 0 if /\/FrameData/;
106	<	$startStunts = 0 if /\/StuntDoubles/;
107	<	}
155	>	def writeFile(outputFileName,repeatX,repeatY,repeatZ):
156	>	outputFile = open(outputFileName, 'w')
157
158	<	print "Affine transformed configuration written to $outName\n";
158	>	outputFile.write("<OpenMD version=1>\n");
159
160	<	if (defined($opts{v})){
161	<	print "New Volume   : $volume\n";
162	<	print "Old Volume   : $oldVolume\n";
163	<	print "Scale Factor : $scale\n";
164	<	}
160	>	print "writing MetaData"
161	>	for metaline in metaData:
162	>	if 'nMol' in metaline:
163	>	metasplit = metaline.split()
164	>	nMol = float(metasplit[2].strip(';'))
165	>	newNmol = nMol * repeatX * repeatY * repeatZ
166	>	outputFile.write('  nMol = %10d;\n' % (newNmol))
167	>	else:
168	>	outputFile.write(metaline)
169
170	+	print "writing Snapshot"
171	+	outputFile.write("  <Snapshot>\n")
172	+
173	+	print "writing FrameData"
174	+	for frameline in frameData:
175	+	if 'Hmat:' in frameline:
176	+	myH = []
177	+	myH.append([repeatX * Hmat[0][0], repeatX * Hmat[0][1], repeatX * Hmat[0][2]])
178	+	myH.append([repeatY * Hmat[1][0], repeatY * Hmat[1][1], repeatY * Hmat[1][2]])
179	+	myH.append([repeatZ * Hmat[2][0], repeatZ * Hmat[2][1], repeatZ * Hmat[2][2]])
180	+	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]))
181	+	else:
182	+	outputFile.write(frameline)
183	+
184	+	print "writing StuntDoubles"
185	+	outputFile.write("    <StuntDoubles>\n")
186	+
187	+	print repeatX, repeatY, repeatZ
188	+
189	+	deco = [ (index, i) for i, index in enumerate(indices) ]
190	+	deco.sort()
191	+	whichSD = 0
192	+	for index in range(len(deco)):
193	+	(index,i) = deco[index]
194	+	print i
195	+	for ii in range(repeatX):
196	+	for jj in range(repeatY):
197	+	for kk in range(repeatZ):
198	+	myP = []
199	+	myP.append(p[i][0] + ii*Hmat[0][0] + jj*Hmat[1][0] + kk*Hmat[2][0])
200	+	myP.append(p[i][1] + ii*Hmat[0][1] + jj*Hmat[1][1] + kk*Hmat[2][1])
201	+	myP.append(p[i][2] + ii*Hmat[0][2] + jj*Hmat[1][2] + kk*Hmat[2][2])
202	+
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	+	elif (pvqj[i] == 'pvqj'):
206	+	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	+	whichSD = whichSD + 1
208	+
209	+	outputFile.write("    </StuntDoubles>\n")
210	+	outputFile.write("  </Snapshot>\n")
211	+	outputFile.write("</OpenMD>\n")
212	+	outputFile.close()
213	+
214	+	def roundMe(x):
215	+	if (x >= 0.0):
216	+	return math.floor(x + 0.5)
217	+	else:
218	+	return math.ceil(x - 0.5)
219	+
220	+	def wrapVector(myVect):
221	+	scaled = [0.0, 0.0, 0.0]
222	+	wrappingNumber = [0,0,0]
223	+	for i in range(3):
224	+	scaled[i] = myVect[i] * BoxInv[i]
225	+	wrappingNumber[i] = roundMe(scaled[i])
226	+	scaled[i] = scaled[i] - wrappingNumber[i]
227	+	myVect[i] = scaled[i] * Hmat[i][i]
228	+	return myVect, wrappingNumber
229	+
230	+	def dot(L1, L2):
231	+	myDot = 0.0
232	+	for i in range(len(L1)):
233	+	myDot = myDot + L1[i]*L2[i]
234	+	return myDot
235	+
236	+	def normalize(L1):
237	+	L2 = []
238	+	myLength = math.sqrt(dot(L1, L1))
239	+	for i in range(len(L1)):
240	+	L2.append(L1[i] / myLength)
241	+	return L2
242	+
243	+	def cross(L1, L2):
244	+	# don't call this with anything other than length 3 lists please
245	+	# or you'll be sorry
246	+	L3 = [0.0, 0.0, 0.0]
247	+	L3[0] = L1[1]*L2[2] - L1[2]*L2[1]
248	+	L3[1] = L1[2]*L2[0] - L1[0]*L2[2]
249	+	L3[2] = L1[0]*L2[1] - L1[1]*L2[0]
250	+	return L3
251	+
252	+	def mapToBox():
253	+	for i in range(len(indices)):
254	+	dpos = []
255	+	dpos.append(p[i][0])
256	+	dpos.append(p[i][1])
257	+	dpos.append(p[i][2])
258	+	p[i], wrap[i] = wrapVector(dpos)
259	+
260	+	def scaleBox(scaleX, scaleY, scaleZ):
261	+	for i in range(3):
262	+	Hmat[0][i] = scaleX * Hmat[0][i]
263	+	for i in range(3):
264	+	Hmat[1][i] = scaleY * Hmat[1][i]
265	+	for i in range(3):
266	+	Hmat[2][i] = scaleZ * Hmat[2][i]
267	+	for i in range(3):
268	+	BoxInv[i] = 1.0/Hmat[i][i]
269	+
270	+	def scaleCoordinates(scaleX, scaleY, scaleZ):
271	+	for i in range(len(indices)):
272	+	p[i][0] = p[i][0]*scaleX + wrap[i][0] * Hmat[0][0]
273	+	p[i][1] = p[i][1]*scaleY + wrap[i][1] * Hmat[1][1]
274	+	p[i][2] = p[i][2]*scaleZ + wrap[i][2] * Hmat[2][2]
275	+
276	+	def main(argv):
277	+	_haveMDFileName = 0
278	+	_haveOutputFileName = 0
279	+	try:
280	+	opts, args = getopt.getopt(argv, "hm:o:x:y:z:v:", ["help", "meta-data=", "output-file=", "newX=", "newY=", "newZ=","newV="])
281	+	except getopt.GetoptError:
282	+	usage()
283	+	sys.exit(2)
284	+	doV = 0
285	+	doX = 0
286	+	doY = 0
287	+	doZ = 0
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	+	_haveMDFileName = 1
295	+	elif opt in ("-o", "--output-file"):
296	+	outputFileName = arg
297	+	_haveOutputFileName = 1
298	+	if opt in ("-x", "--newX"):
299	+	newX = float(arg)
300	+	doX = 1
301	+	elif opt in ("-y", "--newY"):
302	+	newY = float(arg)
303	+	doY = 1
304	+	elif opt in ("-z", "--newZ"):
305	+	newZ = float(arg)
306	+	doZ = 1
307	+	elif opt in ("-v", "--newV"):
308	+	newV = float(arg)
309	+	doV = 1
310	+
311	+	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	+
320	+	if not (doV or doX or doY or doZ):
321	+	usage()
322	+	print "no scaling options given.  Nothing to do!"
323	+	sys.exit()
324	+
325	+	if doV and (doX or doY or doZ):
326	+	usage()
327	+	print "-v is mutually exclusive with any of the -x, -y, and -z options"
328	+	sys.exit()
329	+
330	+	readFile(mdFileName)
331	+
332	+	scaleX = 1.0
333	+	scaleY = 1.0
334	+	scaleZ = 1.0
335	+
336	+	if doX:
337	+	scaleX = newX / Hmat[0][0]
338	+	if doY:
339	+	scaleY = newY / Hmat[1][1]
340	+	if doZ:
341	+	scaleZ = newZ / Hmat[2][2]
342	+
343	+	if doV:
344	+	oldV = Hmat[0][0] * Hmat[1][1] * Hmat[2][2]
345	+	scaleX = pow(newV/oldV, 1.0/3.0)
346	+	scaleY = scaleX
347	+	scaleZ = scaleX
348	+
349	+	mapToBox()
350	+	scaleBox(scaleX, scaleY, scaleZ)
351	+	scaleCoordinates(scaleX, scaleY, scaleZ)
352	+	writeFile(outputFileName, 1, 1, 1)
353	+
354	+	if __name__ == "__main__":
355	+	if len(sys.argv) == 1:
356	+	usage()
357	+	sys.exit()
358	+	main(sys.argv[1:])

Comparing trunk/src/applications/utilities/affineScale (property svn:keywords):
Revision 1063 by chrisfen, Tue Oct 10 14:08:10 2006 UTC vs.
Revision 1919 by jmarr, Wed Jul 31 18:03:35 2013 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines