source: utils/boxmaker.py.in@ bcfb77

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since bcfb77 was bcfb77, checked in by Frederik Heber <heber@…>, 13 years ago

Added correct shebang to boxmaker.py.

  • also configure fills in @PYTHON@.
  • Property mode set to 100755
File size: 10.3 KB
Line 
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
9import re, os, os.path, sys, operator
10
11avogadro = 6.022143e23
12
13class 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
49def 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
88def 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
131def 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
149def 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
171def 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
230def 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
240def 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
291def 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
324def 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
338opt = c_opt()
339
340ReadSettings(opt)
341UpdateSettingsAndSource(opt)
342
343if type(opt.number) == type([]):
344 # Number is a vector - use it without any modification
345 nbox = opt.number
346else:
347 if opt.cubicdomain:
348 nbox = FindBestCube(opt)
349 else:
350 nbox = FindBestCuboid(opt)
351
352VolumePerMolecule = opt.molarmass / (avogadro * opt.density)
353cell = [VolumePerMolecule**(1./3)] * 3
354
355if 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
371print '======== CELL: ', cell
372
373import pyMoleCuilder as mol
374mol.CommandVerbose('0')
375mol.ParserParseTremoloPotentials(opt.potentialsfiledir + opt.basename + '.potentials')
376mol.WorldInput(opt.source)
377mol.WorldCenterInBox('%f 0 0 %f 0 %f' % tuple(cell))
378mol.WorldRepeatBox('%d %d %d' % tuple(nbox))
379mol.WorldOutput(opt.outfilename + '.data')
380mol.WorldOutput(opt.outfilename + '.xyz')
381
382domain = [0.0]*3
383
384for i in range(3):
385 domain[i] = cell[i]*nbox[i]
386
387print '======== DOMAIN: ', domain
388
389# Postprocessing
390
391if 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
403os.system('rm temp_source.data temp_source.xyz')
Note: See TracBrowser for help on using the repository browser.