1 | #!@PYTHON@
|
---|
2 |
|
---|
3 | # Boxmaker 1.0
|
---|
4 | # Creates tremolo-datafiles with arbitrary size and specific density from a single input molecule,
|
---|
5 | # supporting numerous pre- and postprocessing features such as unit conversion.
|
---|
6 | # Gregor Bollerhey - bollerhe@ins.uni-bonn.de
|
---|
7 |
|
---|
8 |
|
---|
9 | import re, os, os.path, sys, operator
|
---|
10 |
|
---|
11 | avogadro = 6.022143e23
|
---|
12 |
|
---|
13 | class c_opt():
|
---|
14 | basename = None
|
---|
15 | tremofiledir = './'
|
---|
16 | potentialsfiledir = './'
|
---|
17 | outfilename = 'out'
|
---|
18 |
|
---|
19 | source = None
|
---|
20 | molarmass = None
|
---|
21 | density = None
|
---|
22 | temp = None
|
---|
23 |
|
---|
24 | number = '1000'
|
---|
25 |
|
---|
26 | cubicdomain = 'on'
|
---|
27 | cubiccell = 'off'
|
---|
28 | autorotate = 'off'
|
---|
29 | autodim = 'on'
|
---|
30 | postprocess = 'on'
|
---|
31 | automass = 'on'
|
---|
32 |
|
---|
33 | def update(self, name, value):
|
---|
34 | shortcuts = {'tf': 'temofiledir', 'pf': 'potentialsfiledir', 'o': 'outfilename',
|
---|
35 | 'i': 'source', 'm': 'molarmass', 'rho': 'density',
|
---|
36 | 't': 'temp', 'n': 'number', 'cd': 'cubicdomain',
|
---|
37 | 'cc': 'cubiccell', 'ar': 'autorotate', 'ad': 'autodim',
|
---|
38 | 'pp': 'postprocess', 'am': 'automass'}
|
---|
39 |
|
---|
40 | if name in shortcuts:
|
---|
41 | name = shortcuts[name]
|
---|
42 |
|
---|
43 | if name in dir(self):
|
---|
44 | exec('self.%s = "%s"' % (name, value))
|
---|
45 | else:
|
---|
46 | print 'Warning: Unknown option:', name
|
---|
47 |
|
---|
48 |
|
---|
49 | def ReadSettings(opt):
|
---|
50 | # Obtain basename
|
---|
51 | if len(sys.argv) >= 2:
|
---|
52 | opt.basename = sys.argv[1]
|
---|
53 | else:
|
---|
54 | print 'Usage: boxmaker.py <basename> [options]'
|
---|
55 | exit()
|
---|
56 |
|
---|
57 | # Read settings file
|
---|
58 | try:
|
---|
59 | with open('boxmaker.' + opt.basename) as f:
|
---|
60 | for line in f:
|
---|
61 | if len(line) > 0 and line[0] != '#':
|
---|
62 | L, S, R = line.partition('=')
|
---|
63 | opt.update(L.strip(), R.strip())
|
---|
64 | except IOError:
|
---|
65 | print 'Warning: Configuration file not readable, CLI only'
|
---|
66 |
|
---|
67 | # Parse parameters
|
---|
68 | i = 2
|
---|
69 | while i < len(sys.argv):
|
---|
70 | L = sys.argv[i]
|
---|
71 |
|
---|
72 | if L[0] in '+-':
|
---|
73 | LN = L[1:]
|
---|
74 |
|
---|
75 | if L[0] == '+':
|
---|
76 | R = 'on'
|
---|
77 | else:
|
---|
78 | R = 'off'
|
---|
79 | else:
|
---|
80 | LN = L
|
---|
81 | i += 1
|
---|
82 | R = sys.argv[i]
|
---|
83 |
|
---|
84 | opt.update(LN, R)
|
---|
85 | i += 1
|
---|
86 |
|
---|
87 |
|
---|
88 | def ReadUnits(opt):
|
---|
89 | lines = [] # The file needs to be processed twice, so we save the lines in the first run
|
---|
90 |
|
---|
91 | with open(opt.tremofiledir + opt.basename + '.tremolo') as f:
|
---|
92 | for line in f:
|
---|
93 | if len(line) > 0 and line[0] != '#':
|
---|
94 | line = line.strip()
|
---|
95 | lines.append(line)
|
---|
96 |
|
---|
97 | if 'systemofunits' in line:
|
---|
98 | L, S, SOU = line.partition('=')
|
---|
99 | SOU = SOU.strip()[:-1] # Remove semicolon
|
---|
100 |
|
---|
101 | if SOU == 'custom':
|
---|
102 | units = {}
|
---|
103 | quantities = ['length', 'mass', 'temperature']
|
---|
104 |
|
---|
105 | for quantity in quantities:
|
---|
106 | units[quantity] = [None, None] # Init with scaling factor and unit 'None'.
|
---|
107 |
|
---|
108 | for line in lines:
|
---|
109 | for quantity in quantities:
|
---|
110 | if quantity in line:
|
---|
111 | L, S, R = line.partition('=')
|
---|
112 | R = R.strip()[:-1] # Remove semicolon
|
---|
113 |
|
---|
114 | if 'scalingfactor' in line:
|
---|
115 | units[quantity][0] = float(R)
|
---|
116 | else:
|
---|
117 | units[quantity][1] = R
|
---|
118 |
|
---|
119 | elif SOU == 'kcalpermole':
|
---|
120 | units = {'length': [1.0, 'angstrom'], 'mass': [1.0, 'u'], 'temperature': [503.556, 'K']}
|
---|
121 |
|
---|
122 | elif SOU == 'evolt':
|
---|
123 | units = {'length': [1.0, 'angstrom'], 'mass': [1.0, 'u'], 'temperature': [11604.0, 'K']}
|
---|
124 |
|
---|
125 | else: # SI
|
---|
126 | units = {'length': [1.0, 'm'], 'mass': [1.0, 'kg'], 'temperature': [1.0, 'K']}
|
---|
127 |
|
---|
128 | return units
|
---|
129 |
|
---|
130 |
|
---|
131 | def ConvertUnits(have, want):
|
---|
132 | if have[0] == '!':
|
---|
133 | return float(have[1:])
|
---|
134 |
|
---|
135 | # Redo with pipes?
|
---|
136 | ret = os.system("units '%s' '%s' > temp_units_output" % (have, want))
|
---|
137 |
|
---|
138 | if ret == 0:
|
---|
139 | with open('temp_units_output') as f:
|
---|
140 | line = f.readline()
|
---|
141 |
|
---|
142 | os.system('rm temp_units_output')
|
---|
143 |
|
---|
144 | return float(line[3:-1])
|
---|
145 | else:
|
---|
146 | raise NameError('UnitError')
|
---|
147 |
|
---|
148 |
|
---|
149 | def GetSourceMolareMass(opt):
|
---|
150 | with open(opt.potentialsfiledir+opt.basename+'.potentials') as f:
|
---|
151 | potfile = f.read()
|
---|
152 |
|
---|
153 | elementmasses = dict(re.findall(r'element_name=(\w{1,2}).*?mass=([0-9.]*)', potfile))
|
---|
154 |
|
---|
155 | for key in elementmasses:
|
---|
156 | elementmasses[key] = float(elementmasses[key])
|
---|
157 |
|
---|
158 | mass_sum = 0.0
|
---|
159 |
|
---|
160 | with open('temp_source.xyz') as f:
|
---|
161 | N = int(f.readline())
|
---|
162 | comment = f.readline()
|
---|
163 |
|
---|
164 | for i in range(N):
|
---|
165 | elem = f.readline().split(None, 1)[0].strip()
|
---|
166 | mass_sum += elementmasses[elem]
|
---|
167 |
|
---|
168 | return mass_sum*avogadro
|
---|
169 |
|
---|
170 |
|
---|
171 | def UpdateSettingsAndSource(opt):
|
---|
172 | # Map boolean values
|
---|
173 | boolmap = {'on': True, 'off': False}
|
---|
174 |
|
---|
175 | for name in ['cubicdomain', 'cubiccell', 'autorotate', 'autodim', 'postprocess', 'automass']:
|
---|
176 | value = eval('opt.' + name)
|
---|
177 |
|
---|
178 | if value in boolmap:
|
---|
179 | value = boolmap[value]
|
---|
180 | else:
|
---|
181 | print 'Not a boolean value:', value
|
---|
182 | exit()
|
---|
183 |
|
---|
184 | exec('opt.' + name + '= value')
|
---|
185 |
|
---|
186 | # Convert dimensions
|
---|
187 | if opt.autodim:
|
---|
188 | units = ReadUnits(opt)
|
---|
189 |
|
---|
190 | if not opt.automass:
|
---|
191 | have = opt.molarmass
|
---|
192 | want = '%f*%s / mol' % tuple(units['mass'])
|
---|
193 | opt.molarmass = ConvertUnits(have, want)
|
---|
194 |
|
---|
195 | have = opt.density
|
---|
196 | want = '(%f*%s) ' % tuple(units['mass']) + '/ (%f*%s)**3' % tuple(units['length'])
|
---|
197 | opt.density = ConvertUnits(have, want)
|
---|
198 |
|
---|
199 | if opt.temp:
|
---|
200 | have = opt.temp
|
---|
201 | want = '%f*%s' % tuple(units['temperature'])
|
---|
202 | opt.temp = ConvertUnits(have, want)
|
---|
203 | else:
|
---|
204 | if not opt.automass:
|
---|
205 | opt.molarmass = float(opt.molarmass)
|
---|
206 |
|
---|
207 | opt.density = float(opt.density)
|
---|
208 |
|
---|
209 | if opt.temp:
|
---|
210 | opt.temp = float(opt.temp)
|
---|
211 |
|
---|
212 | # Number might be an integer or a 3-vector
|
---|
213 | nvec = opt.number.split()
|
---|
214 | if len(nvec) == 3:
|
---|
215 | opt.number = [0]*3
|
---|
216 |
|
---|
217 | for i in range(3):
|
---|
218 | opt.number[i] = int(nvec[i])
|
---|
219 | else:
|
---|
220 | opt.number = int(opt.number)
|
---|
221 |
|
---|
222 | UpdateSource(opt)
|
---|
223 |
|
---|
224 | # Automatic source mass
|
---|
225 | if opt.automass:
|
---|
226 | opt.molarmass = GetSourceMolareMass(opt)
|
---|
227 | print '======== MOLAR MASS:', opt.molarmass
|
---|
228 |
|
---|
229 |
|
---|
230 | def FindBestCube(opt):
|
---|
231 | newroot = int( round(opt.number**(1./3)) )
|
---|
232 | newnumber = newroot**3
|
---|
233 |
|
---|
234 | if newnumber != opt.number:
|
---|
235 | print 'Warning: Number changed to %d.' % newnumber
|
---|
236 |
|
---|
237 | return [newroot] * 3
|
---|
238 |
|
---|
239 |
|
---|
240 | def FindBestCuboid(opt):
|
---|
241 | n = opt.number
|
---|
242 |
|
---|
243 | # Prime factors of n
|
---|
244 | factors = []
|
---|
245 |
|
---|
246 | for i in [2, 3]:
|
---|
247 | while n % i == 0:
|
---|
248 | factors.append(i)
|
---|
249 | n /= 2
|
---|
250 |
|
---|
251 | t = 5
|
---|
252 | diff = 2
|
---|
253 |
|
---|
254 | while t*t <= n:
|
---|
255 | while n % t == 0:
|
---|
256 | factors.append(t)
|
---|
257 | n /= t
|
---|
258 |
|
---|
259 | t = t + diff
|
---|
260 | diff = 6 - diff
|
---|
261 |
|
---|
262 | if n > 1:
|
---|
263 | factors.append(n)
|
---|
264 |
|
---|
265 | # Even distribution of current biggest prime to each vector -> similar sizes
|
---|
266 | if len(factors) < 3:
|
---|
267 | print 'Warning: Not enough prime factors - falling back to cubic placement'
|
---|
268 | return FindBestCube(opt)
|
---|
269 |
|
---|
270 | factors.sort()
|
---|
271 | distri = [[],[],[]]
|
---|
272 | current = 0
|
---|
273 |
|
---|
274 | for factor in factors:
|
---|
275 | distri[current].append(factor)
|
---|
276 | current += 1
|
---|
277 | if current == 3:
|
---|
278 | current = 0
|
---|
279 |
|
---|
280 | result = [0]*3
|
---|
281 |
|
---|
282 | print '======== CUBOID USED:',
|
---|
283 |
|
---|
284 | for i in range(3):
|
---|
285 | result[i] = int( reduce(operator.mul, distri[i]) )
|
---|
286 |
|
---|
287 | print result
|
---|
288 | return result
|
---|
289 |
|
---|
290 |
|
---|
291 | def GetSourceBBabs(opt):
|
---|
292 | bbmax = [0.0]*3
|
---|
293 | bbmin = [float('inf')]*3
|
---|
294 |
|
---|
295 | s_name_ext = os.path.basename(opt.source).rsplit('.', 1)
|
---|
296 | s_namepart = s_name_ext[0]
|
---|
297 |
|
---|
298 | if len(s_name_ext) > 1:
|
---|
299 | s_ext = s_name_ext[1]
|
---|
300 | else:
|
---|
301 | s_ext = ''
|
---|
302 |
|
---|
303 | # Calculate bounding box from xyz-file
|
---|
304 | with open('temp_source.xyz') as f:
|
---|
305 | N = int(f.readline())
|
---|
306 | comment = f.readline()
|
---|
307 |
|
---|
308 | for i in xrange(N):
|
---|
309 | buf = f.readline()
|
---|
310 | xyz = buf.split()[1:]
|
---|
311 |
|
---|
312 | for i in range(3):
|
---|
313 | bbmax[i] = max(bbmax[i], float(xyz[i]))
|
---|
314 | bbmin[i] = min(bbmin[i], float(xyz[i]))
|
---|
315 |
|
---|
316 | bb = [0.0]*3
|
---|
317 |
|
---|
318 | for i in range(3):
|
---|
319 | bb[i] = abs(bbmax[i] - bbmin[i])
|
---|
320 |
|
---|
321 | return bb
|
---|
322 |
|
---|
323 |
|
---|
324 | def UpdateSource(opt):
|
---|
325 | potfilepath = opt.potentialsfiledir + opt.basename + '.potentials'
|
---|
326 |
|
---|
327 | cmd = 'molecuilder -o xyz tremolo --parse-tremolo-potentials %s -i temp_source.xyz -l %s' % (potfilepath, opt.source)
|
---|
328 |
|
---|
329 | if opt.autorotate:
|
---|
330 | cmd += ' --select-all-atoms --rotate-to-principal-axis-system "0, 1, 0"'
|
---|
331 |
|
---|
332 | os.system(cmd)
|
---|
333 |
|
---|
334 | opt.source = 'temp_source.data'
|
---|
335 |
|
---|
336 |
|
---|
337 | # Global options with sensible default parameters
|
---|
338 | opt = c_opt()
|
---|
339 |
|
---|
340 | ReadSettings(opt)
|
---|
341 | UpdateSettingsAndSource(opt)
|
---|
342 |
|
---|
343 | if type(opt.number) == type([]):
|
---|
344 | # Number is a vector - use it without any modification
|
---|
345 | nbox = opt.number
|
---|
346 | else:
|
---|
347 | if opt.cubicdomain:
|
---|
348 | nbox = FindBestCube(opt)
|
---|
349 | else:
|
---|
350 | nbox = FindBestCuboid(opt)
|
---|
351 |
|
---|
352 | VolumePerMolecule = opt.molarmass / (avogadro * opt.density)
|
---|
353 | cell = [VolumePerMolecule**(1./3)] * 3
|
---|
354 |
|
---|
355 | if not opt.cubiccell:
|
---|
356 | try:
|
---|
357 | bb = GetSourceBBabs(opt)
|
---|
358 | print '======== BBOX:', bb
|
---|
359 | # Scaling factor - the molecules bounding box is scaled to fit the volume suiting the density
|
---|
360 | s = (VolumePerMolecule / (bb[0]*bb[1]*bb[2])) ** (1./3)
|
---|
361 |
|
---|
362 | if s < 1:
|
---|
363 | print 'Warning: Molecular cells will overlap.'
|
---|
364 |
|
---|
365 | for i in range(3):
|
---|
366 | cell[i] = bb[i]*s
|
---|
367 | except ZeroDivisionError:
|
---|
368 | print 'Warning: Singularity in bounding box, falling back to cubic cell.'
|
---|
369 |
|
---|
370 |
|
---|
371 | print '======== CELL: ', cell
|
---|
372 |
|
---|
373 | import pyMoleCuilder as mol
|
---|
374 | mol.CommandVerbose('0')
|
---|
375 | mol.ParserParseTremoloPotentials(opt.potentialsfiledir + opt.basename + '.potentials')
|
---|
376 | mol.WorldInput(opt.source)
|
---|
377 | mol.WorldCenterInBox('%f 0 0 %f 0 %f' % tuple(cell))
|
---|
378 | mol.WorldRepeatBox('%d %d %d' % tuple(nbox))
|
---|
379 | mol.WorldOutput(opt.outfilename + '.data')
|
---|
380 | mol.WorldOutput(opt.outfilename + '.xyz')
|
---|
381 |
|
---|
382 | domain = [0.0]*3
|
---|
383 |
|
---|
384 | for i in range(3):
|
---|
385 | domain[i] = cell[i]*nbox[i]
|
---|
386 |
|
---|
387 | print '======== DOMAIN: ', domain
|
---|
388 |
|
---|
389 | # Postprocessing
|
---|
390 |
|
---|
391 | if opt.postprocess:
|
---|
392 | with open(opt.outfilename + '.data') as f:
|
---|
393 | ofile = f.read()
|
---|
394 |
|
---|
395 | with open(opt.outfilename + '.data', 'w') as f:
|
---|
396 | f.write('# INPUTCONV shift center\n')
|
---|
397 |
|
---|
398 | if opt.temp:
|
---|
399 | f.write('# INPUTCONV temp %.4f\n' % opt.temp)
|
---|
400 |
|
---|
401 | f.write(ofile)
|
---|
402 |
|
---|
403 | os.system('rm temp_source.data temp_source.xyz')
|
---|