source: utils/boxmaker.py@ bfbb62

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 bfbb62 was bfbb62, checked in by Frederik Heber <heber@…>, 13 years ago

Final version 1.0

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