Coverage for /builds/debichem-team/python-ase/ase/calculators/openmx/reader.py : 70.22%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# flake8: noqa
2"""
3The ASE Calculator for OpenMX <http://www.openmx-square.org>: Python interface
4to the software package for nano-scale material simulations based on density
5functional theories.
6 Copyright (C) 2018 JaeHwan Shim and JaeJun Yu
8 This program is free software: you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation, either version 2.1 of the License, or
11 (at your option) any later version.
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU Lesser General Public License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with ASE. If not, see <http://www.gnu.org/licenses/>.
20"""
21# from ase.calculators import SinglePointDFTCalculator
22import os
23import struct
24import numpy as np
25from ase.units import Ha, Bohr, Debye
26from ase.io import ParseError
30def read_openmx(filename=None, debug=False):
31 from ase.calculators.openmx import OpenMX
32 from ase import Atoms
33 """
34 Read results from typical OpenMX output files and returns the atom object
35 In default mode, it reads every implementd properties we could get from
36 the files. Unlike previous version, we read the information based on file.
37 previous results will be eraised unless previous results are written in the
38 next calculation results.
40 Read the 'LABEL.log' file seems redundant. Because all the
41 information should already be written in '.out' file. However, in the
42 version 3.8.3, stress tensor are not written in the '.out' file. It only
43 contained in the '.log' file. So... I implented reading '.log' file method
44 """
45 log_data = read_file(get_file_name('.log', filename), debug=debug)
46 restart_data = read_file(get_file_name('.dat#', filename), debug=debug)
47 dat_data = read_file(get_file_name('.dat', filename), debug=debug)
48 out_data = read_file(get_file_name('.out', filename), debug=debug)
49 scfout_data = read_scfout_file(get_file_name('.scfout', filename))
50 band_data = read_band_file(get_file_name('.Band', filename))
51 # dos_data = read_dos_file(get_file_name('.Dos.val', filename))
52 """
53 First, we get every data we could get from the all results files. And then,
54 reform the data to fit to data structure of Atom object. While doing this,
55 Fix the unit to ASE format.
56 """
57 parameters = get_parameters(out_data=out_data, log_data=log_data,
58 restart_data=restart_data, dat_data=dat_data,
59 scfout_data=scfout_data, band_data=band_data)
60 atomic_formula = get_atomic_formula(out_data=out_data, log_data=log_data,
61 restart_data=restart_data,
62 scfout_data=scfout_data,
63 dat_data=dat_data)
64 results = get_results(out_data=out_data, log_data=log_data,
65 restart_data=restart_data, scfout_data=scfout_data,
66 dat_data=dat_data, band_data=band_data)
68 atoms = Atoms(**atomic_formula)
69 atoms.calc = OpenMX(**parameters)
70 atoms.calc.results = results
71 return atoms
74def read_file(filename, debug=False):
75 """
76 Read the 'LABEL.out' file. Using 'parameters.py', we read every 'allowed_
77 dat' dictionory. while reading a file, if one find the key matcheds That
78 'patters', which indicates the property we want is written, it will returns
79 the pair value of that key. For example,
80 example will be written later
81 """
82 from ase.calculators.openmx import parameters as param
83 if not os.path.isfile(filename):
84 return {}
85 param_keys = ['integer_keys', 'float_keys', 'string_keys', 'bool_keys',
86 'list_int_keys', 'list_float_keys', 'list_bool_keys',
87 'tuple_integer_keys', 'tuple_float_keys', 'tuple_float_keys']
88 patterns = {
89 'Stress tensor': ('stress', read_stress_tensor),
90 'Dipole moment': ('dipole', read_dipole),
91 'Fractional coordinates of': ('scaled_positions', read_scaled_positions),
92 'Utot.': ('energy', read_energy),
93 'energies in': ('energies', read_energies),
94 'Chemical Potential': ('chemical_potential', read_chemical_potential),
95 '<coordinates.forces': ('forces', read_forces),
96 'Eigenvalues (Hartree)': ('eigenvalues', read_eigenvalues)}
97 special_patterns = {
98 'Total spin moment': (('magmoms', 'total_magmom'),
99 read_magmoms_and_total_magmom),
100 }
101 out_data = {}
102 line = '\n'
103 if(debug):
104 print('Read results from %s' % filename)
105 with open(filename, 'r') as fd:
106 '''
107 Read output file line by line. When the `line` matches the pattern
108 of certain keywords in `param.[dtype]_keys`, for example,
110 if line in param.string_keys:
111 out_data[key] = read_string(line)
113 parse that line and store it to `out_data` in specified data type.
114 To cover all `dtype` parameters, for loop was used,
116 for [dtype] in parameters_keys:
117 if line in param.[dtype]_keys:
118 out_data[key] = read_[dtype](line)
120 After found matched pattern, escape the for loop using `continue`.
121 '''
122 while line != '':
123 pattern_matched = False
124 line = fd.readline()
125 try:
126 _line = line.split()[0]
127 except IndexError:
128 continue
129 for dtype_key in param_keys:
130 dtype = dtype_key.rsplit('_', 1)[0]
131 read_dtype = globals()['read_' + dtype]
132 for key in param.__dict__[dtype_key]:
133 if key in _line:
134 out_data[get_standard_key(key)] = read_dtype(line)
135 pattern_matched = True
136 continue
137 if pattern_matched:
138 continue
140 for key in param.matrix_keys:
141 if '<'+key in line:
142 out_data[get_standard_key(key)] = read_matrix(line, key, fd)
143 pattern_matched = True
144 continue
145 if pattern_matched:
146 continue
147 for key in patterns.keys():
148 if key in line:
149 out_data[patterns[key][0]] = patterns[key][1](line, fd, debug=debug)
150 pattern_matched = True
151 continue
152 if pattern_matched:
153 continue
154 for key in special_patterns.keys():
155 if key in line:
156 a, b = special_patterns[key][1](line, fd)
157 out_data[special_patterns[key][0][0]] = a
158 out_data[special_patterns[key][0][1]] = b
159 pattern_matched = True
160 continue
161 if pattern_matched:
162 continue
163 return out_data
166def read_scfout_file(filename=None):
167 """
168 Read the Developer output '.scfout' files. It Behaves like read_scfout.c,
169 OpenMX module, but written in python. Note that some array are begin with
170 1, not 0
172 atomnum: the number of total atoms
173 Catomnum: the number of atoms in the central region
174 Latomnum: the number of atoms in the left lead
175 Ratomnum: the number of atoms in the left lead
176 SpinP_switch:
177 0: non-spin polarized
178 1: spin polarized
179 TCpyCell: the total number of periodic cells
180 Solver: method for solving eigenvalue problem
181 ChemP: chemical potential
182 Valence_Electrons: total number of valence electrons
183 Total_SpinS: total value of Spin (2*Total_SpinS = muB)
184 E_Temp: electronic temperature
185 Total_NumOrbs: the number of atomic orbitals in each atom
186 size: Total_NumOrbs[atomnum+1]
187 FNAN: the number of first neighboring atoms of each atom
188 size: FNAN[atomnum+1]
189 natn: global index of neighboring atoms of an atom ct_AN
190 size: natn[atomnum+1][FNAN[ct_AN]+1]
191 ncn: global index for cell of neighboring atoms of an atom ct_AN
192 size: ncn[atomnum+1][FNAN[ct_AN]+1]
193 atv: x,y,and z-components of translation vector of periodically copied cell
194 size: atv[TCpyCell+1][4]:
195 atv_ijk: i,j,and j number of periodically copied cells
196 size: atv_ijk[TCpyCell+1][4]:
197 tv[4][4]: unit cell vectors in Bohr
198 rtv[4][4]: reciprocal unit cell vectors in Bohr^{-1}
199 note:
200 tv_i dot rtv_j = 2PI * Kronecker's delta_{ij}
201 Gxyz[atomnum+1][60]: atomic coordinates in Bohr
202 Hks: Kohn-Sham matrix elements of basis orbitals
203 size: Hks[SpinP_switch+1]
204 [atomnum+1]
205 [FNAN[ct_AN]+1]
206 [Total_NumOrbs[ct_AN]]
207 [Total_NumOrbs[h_AN]]
208 iHks:
209 imaginary Kohn-Sham matrix elements of basis orbitals
210 for alpha-alpha, beta-beta, and alpha-beta spin matrices
211 of which contributions come from spin-orbit coupling
212 and Hubbard U effective potential.
213 size: iHks[3]
214 [atomnum+1]
215 [FNAN[ct_AN]+1]
216 [Total_NumOrbs[ct_AN]]
217 [Total_NumOrbs[h_AN]]
218 OLP: overlap matrix
219 size: OLP[atomnum+1]
220 [FNAN[ct_AN]+1]
221 [Total_NumOrbs[ct_AN]]
222 [Total_NumOrbs[h_AN]]
223 OLPpox: overlap matrix with position operator x
224 size: OLPpox[atomnum+1]
225 [FNAN[ct_AN]+1]
226 [Total_NumOrbs[ct_AN]]
227 [Total_NumOrbs[h_AN]]
228 OLPpoy: overlap matrix with position operator y
229 size: OLPpoy[atomnum+1]
230 [FNAN[ct_AN]+1]
231 [Total_NumOrbs[ct_AN]]
232 [Total_NumOrbs[h_AN]]
233 OLPpoz: overlap matrix with position operator z
234 size: OLPpoz[atomnum+1]
235 [FNAN[ct_AN]+1]
236 [Total_NumOrbs[ct_AN]]
237 [Total_NumOrbs[h_AN]]
238 DM: overlap matrix
239 size: DM[SpinP_switch+1]
240 [atomnum+1]
241 [FNAN[ct_AN]+1]
242 [Total_NumOrbs[ct_AN]]
243 [Total_NumOrbs[h_AN]]
244 dipole_moment_core[4]:
245 dipole_moment_background[4]:
246 """
247 from numpy import insert as ins
248 from numpy import cumsum as cum
249 from numpy import split as spl
250 from numpy import sum, zeros
251 if not os.path.isfile(filename):
252 return {}
254 def easyReader(byte, data_type, shape):
255 data_size = {'d': 8, 'i': 4}
256 data_struct = {'d': float, 'i': int}
257 dt = data_type
258 ds = data_size[data_type]
259 unpack = struct.unpack
260 if len(byte) == ds:
261 if dt == 'i':
262 return data_struct[dt].from_bytes(byte, byteorder='little')
263 elif dt == 'd':
264 return np.array(unpack(dt*(len(byte)//ds), byte))[0]
265 elif shape is not None:
266 return np.array(unpack(dt*(len(byte)//ds), byte)).reshape(shape)
267 else:
268 return np.array(unpack(dt*(len(byte)//ds), byte))
270 def inte(byte, shape=None):
271 return easyReader(byte, 'i', shape)
273 def floa(byte, shape=None):
274 return easyReader(byte, 'd', shape)
276 def readOverlap(atomnum, Total_NumOrbs, FNAN, natn, f):
277 myOLP = []
278 myOLP.append([])
279 for ct_AN in range(1, atomnum + 1):
280 myOLP.append([])
281 TNO1 = Total_NumOrbs[ct_AN]
282 for h_AN in range(FNAN[ct_AN] + 1):
283 myOLP[ct_AN].append([])
284 Gh_AN = natn[ct_AN][h_AN]
285 TNO2 = Total_NumOrbs[Gh_AN]
286 for i in range(TNO1):
287 myOLP[ct_AN][h_AN].append(floa(f.read(8*TNO2)))
288 return myOLP
290 def readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, f):
291 Hks = []
292 for spin in range(SpinP_switch + 1):
293 Hks.append([])
294 Hks[spin].append([np.zeros(FNAN[0] + 1)])
295 for ct_AN in range(1, atomnum + 1):
296 Hks[spin].append([])
297 TNO1 = Total_NumOrbs[ct_AN]
298 for h_AN in range(FNAN[ct_AN] + 1):
299 Hks[spin][ct_AN].append([])
300 Gh_AN = natn[ct_AN][h_AN]
301 TNO2 = Total_NumOrbs[Gh_AN]
302 for i in range(TNO1):
303 Hks[spin][ct_AN][h_AN].append(floa(f.read(8*TNO2)))
304 return Hks
306 fd = open(filename, mode='rb')
307 atomnum, SpinP_switch = inte(fd.read(8))
308 Catomnum, Latomnum, Ratomnum, TCpyCell = inte(fd.read(16))
309 atv = floa(fd.read(8*4*(TCpyCell+1)), shape=(TCpyCell+1, 4))
310 atv_ijk = inte(fd.read(4*4*(TCpyCell+1)), shape=(TCpyCell+1, 4))
311 Total_NumOrbs = np.insert(inte(fd.read(4*(atomnum))), 0, 1, axis=0)
312 FNAN = np.insert(inte(fd.read(4*(atomnum))), 0, 0, axis=0)
313 natn = ins(spl(inte(fd.read(4*sum(FNAN[1:] + 1))), cum(FNAN[1:] + 1)),
314 0, zeros(FNAN[0] + 1), axis=0)[:-1]
315 ncn = ins(spl(inte(fd.read(4*np.sum(FNAN[1:] + 1))), cum(FNAN[1:] + 1)),
316 0, np.zeros(FNAN[0] + 1), axis=0)[:-1]
317 tv = ins(floa(fd.read(8*3*4), shape=(3, 4)), 0, [0, 0, 0, 0], axis=0)
318 rtv = ins(floa(fd.read(8*3*4), shape=(3, 4)), 0, [0, 0, 0, 0], axis=0)
319 Gxyz = ins(floa(fd.read(8*(atomnum)*4), shape=(atomnum, 4)), 0,
320 [0., 0., 0., 0.], axis=0)
321 Hks = readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, fd)
322 iHks = []
323 if SpinP_switch == 3:
324 iHks = readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, fd)
325 OLP = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd)
326 OLPpox = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd)
327 OLPpoy = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd)
328 OLPpoz = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd)
329 DM = readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, fd)
330 Solver = inte(fd.read(4))
331 ChemP, E_Temp = floa(fd.read(8*2))
332 dipole_moment_core = floa(fd.read(8*3))
333 dipole_moment_background = floa(fd.read(8*3))
334 Valence_Electrons, Total_SpinS = floa(fd.read(8*2))
336 fd.close()
337 scf_out = {'atomnum': atomnum, 'SpinP_switch': SpinP_switch,
338 'Catomnum': Catomnum, 'Latomnum': Latomnum, 'Hks': Hks,
339 'Ratomnum': Ratomnum, 'TCpyCell': TCpyCell, 'atv': atv,
340 'Total_NumOrbs': Total_NumOrbs, 'FNAN': FNAN, 'natn': natn,
341 'ncn': ncn, 'tv': tv, 'rtv': rtv, 'Gxyz': Gxyz, 'OLP': OLP,
342 'OLPpox': OLPpox, 'OLPpoy': OLPpoy, 'OLPpoz': OLPpoz,
343 'Solver': Solver, 'ChemP': ChemP, 'E_Temp': E_Temp,
344 'dipole_moment_core': dipole_moment_core, 'iHks': iHks,
345 'dipole_moment_background': dipole_moment_background,
346 'Valence_Electrons': Valence_Electrons, 'atv_ijk': atv_ijk,
347 'Total_SpinS': Total_SpinS, 'DM': DM
348 }
349 return scf_out
352def read_band_file(filename=None):
353 band_data = {}
354 if not os.path.isfile(filename):
355 return {}
356 band_kpath = []
357 eigen_bands = []
358 with open(filename, 'r') as fd:
359 line = f.readline().split()
360 nkpts = 0
361 nband = int(line[0])
362 nspin = int(line[1]) + 1
363 band_data['nband'] = nband
364 band_data['nspin'] = nspin
365 line = f.readline().split()
366 band_data['band_kpath_unitcell'] = [line[:3], line[3:6], line[6:9]]
367 line = f.readline().split()
368 band_data['band_nkpath'] = int(line[0])
369 for i in range(band_data['band_nkpath']):
370 line = f.readline().split()
371 band_kpath.append(line)
372 nkpts += int(line[0])
373 band_data['nkpts'] = nkpts
374 band_data['band_kpath'] = band_kpath
375 kpts = np.zeros((nkpts, 3))
376 eigen_bands = np.zeros((nspin, nkpts, nband))
377 for i in range(nspin):
378 for j in range(nkpts):
379 line = f.readline()
380 kpts[j] = np.array(line.split(), dtype=float)[1:]
381 line = f.readline()
382 eigen_bands[i, j] = np.array(line.split(), dtype=float)[:]
383 band_data['eigenvalues'] = eigen_bands
384 band_data['band_kpts'] = kpts
385 return band_data
388def read_electron_valency(filename='H_CA13'):
389 array = []
390 with open(filename, 'r') as fd:
391 array = fd.readlines()
392 fd.close()
393 required_line = ''
394 for line in array:
395 if 'valence.electron' in line:
396 required_line = line
397 return rn(required_line)
400def rn(line='\n', n=1):
401 """
402 Read n'th to last value.
403 For example:
404 ...
405 scf.XcType LDA
406 scf.Kgrid 4 4 4
407 ...
408 In Python,
409 >>> str(rn(line, 1))
410 LDA
411 >>> line = f.readline()
412 >>> int(rn(line, 3))
413 4
414 """
415 return line.split()[-n]
418def read_tuple_integer(line):
419 return tuple([int(x) for x in line.split()[-3:]])
422def read_tuple_float(line):
423 return tuple([float(x) for x in line.split()[-3:]])
426def read_integer(line):
427 return int(rn(line))
430def read_float(line):
431 return float(rn(line))
434def read_string(line):
435 return str(rn(line))
438def read_bool(line):
439 bool = str(rn(line)).lower()
440 if bool == 'on':
441 return True
442 elif bool == 'off':
443 return False
444 else:
445 return None
448def read_list_int(line):
449 return [int(x) for x in line.split()[1:]]
452def read_list_float(line):
453 return [float(x) for x in line.split()[1:]]
456def read_list_bool(line):
457 return [read_bool(x) for x in line.split()[1:]]
460def read_matrix(line, key, f):
461 matrix = []
462 line = f.readline()
463 while key not in line:
464 matrix.append(line.split())
465 line = f.readline()
466 return matrix
469def read_stress_tensor(line, f, debug=None):
470 f.readline() # passing empty line
471 f.readline()
472 line = f.readline()
473 xx, xy, xz = read_tuple_float(line)
474 line = f.readline()
475 yx, yy, yz = read_tuple_float(line)
476 line = f.readline()
477 zx, zy, zz = read_tuple_float(line)
478 stress = [xx, yy, zz, (zy + yz)/2, (zx + xz)/2, (yx + xy)/2]
479 return stress
482def read_magmoms_and_total_magmom(line, f, debug=None):
483 total_magmom = read_float(line)
484 f.readline() # Skip empty lines
485 f.readline()
486 line = f.readline()
487 magmoms = []
488 while not(line == '' or line.isspace()):
489 magmoms.append(read_float(line))
490 line = f.readline()
491 return magmoms, total_magmom
494def read_energy(line, f, debug=None):
495 # It has Hartree unit yet
496 return read_float(line)
498def read_energies(line, f, debug=None):
499 line = f.readline()
500 if '***' in line:
501 point = 7 # Version 3.8
502 else:
503 point = 16 # Version 3.9
504 for i in range(point):
505 f.readline()
506 line = f.readline()
507 energies = []
508 while not(line == '' or line.isspace()):
509 energies.append(float(line.split()[2]))
510 line = f.readline()
511 return energies
513def read_eigenvalues(line, f, debug=False):
514 """
515 Read the Eigenvalues in the `.out` file and returns the eigenvalue
516 First, it assumes system have two spins and start reading until it reaches
517 the end('*****...').
519 eigenvalues[spin][kpoint][nbands]
521 For symmetry reason, `.out` file prints the eigenvalues at the half of the
522 K points. Thus, we have to fill up the rest of the half.
523 However, if the calculation was conducted only on the gamma point, it will
524 raise the 'gamma_flag' as true and it will returns the original samples.
525 """
526 def prind(*line, end='\n'):
527 if debug:
528 print(*line, end=end)
529 prind("Read eigenvalues output")
530 current_line = f.tell()
531 f.seek(0) # Seek for the kgrid information
532 while line != '':
533 line = f.readline().lower()
534 if 'scf.kgrid' in line:
535 break
536 f.seek(current_line) # Retrun to the original position
538 kgrid = read_tuple_integer(line)
540 if kgrid != ():
541 prind('Non-Gamma point calculation')
542 prind('scf.Kgrid is %d, %d, %d' % kgrid)
543 gamma_flag = False
544 # f.seek(f.tell()+57)
545 else:
546 prind('Gamma point calculation')
547 gamma_flag = True
548 line = f.readline()
549 line = f.readline()
551 eigenvalues = []
552 eigenvalues.append([])
553 eigenvalues.append([]) # Assume two spins
554 i = 0
555 while True:
556 # Go to eigenvalues line
557 while line != '':
558 line = f.readline()
559 prind(line)
560 ll = line.split()
561 if line.isspace():
562 continue
563 elif len(ll) > 1:
564 if ll[0] == '1':
565 break
566 elif "*****" in line:
567 break
569 # Break if it reaches the end or next parameters
570 if "*****" in line or line == '':
571 break
573 # Check Number and format is valid
574 try:
575 # Suppose to be looks like
576 # 1 -2.33424746491277 -2.33424746917880
577 ll = line.split()
578 # Check if it reaches the end of the file
579 assert line != ''
580 assert len(ll) == 3
581 float(ll[1]); float(ll[2])
582 except (AssertionError, ValueError):
583 raise ParseError("Cannot read eigenvalues")
585 # Read eigenvalues
586 eigenvalues[0].append([])
587 eigenvalues[1].append([])
588 while not (line == '' or line.isspace()):
589 eigenvalues[0][i].append(float(rn(line, 2)))
590 eigenvalues[1][i].append(float(rn(line, 1)))
591 line = f.readline()
592 prind(line, end='')
593 i += 1
594 prind(line)
595 if gamma_flag:
596 return np.asarray(eigenvalues)
597 eigen_half = np.asarray(eigenvalues)
598 prind(eigen_half)
599 # Fill up the half
600 spin, half_kpts, bands = eigen_half.shape
601 even_odd = np.array(kgrid).prod() % 2
602 eigen_values = np.zeros((spin, half_kpts*2-even_odd, bands))
603 for i in range(half_kpts):
604 eigen_values[0, i] = eigen_half[0, i, :]
605 eigen_values[1, i] = eigen_half[1, i, :]
606 eigen_values[0, 2*half_kpts-1-i-even_odd] = eigen_half[0, i, :]
607 eigen_values[1, 2*half_kpts-1-i-even_odd] = eigen_half[1, i, :]
608 return eigen_values
611def read_forces(line, f, debug=None):
612 # It has Hartree per Bohr unit yet
613 forces = []
614 f.readline() # Skip Empty line
615 line = f.readline()
616 while 'coordinates.forces>' not in line:
617 forces.append(read_tuple_float(line))
618 line = f.readline()
619 return np.array(forces)
622def read_dipole(line, f, debug=None):
623 dipole = []
624 while 'Total' not in line:
625 line = f.readline()
626 dipole.append(read_tuple_float(line))
627 return dipole
630def read_scaled_positions(line, f, debug=None):
631 scaled_positions = []
632 f.readline() # Skip Empty lines
633 f.readline()
634 f.readline()
635 line = f.readline()
636 while not(line == '' or line.isspace()): # Detect empty line
637 scaled_positions.append(read_tuple_float(line))
638 line = f.readline()
639 return scaled_positions
642def read_chemical_potential(line, f, debug=None):
643 return read_float(line)
646def get_parameters(out_data=None, log_data=None, restart_data=None,
647 scfout_data=None, dat_data=None, band_data=None):
648 """
649 From the given data sets, construct the dictionary 'parameters'. If data
650 is in the paramerters, it will save it.
651 """
652 from ase.calculators.openmx import parameters as param
653 scaned_data = [dat_data, out_data, log_data, restart_data, scfout_data,
654 band_data]
655 openmx_keywords = [param.tuple_integer_keys, param.tuple_float_keys,
656 param.tuple_bool_keys, param.integer_keys,
657 param.float_keys, param.string_keys, param.bool_keys,
658 param.list_int_keys, param.list_bool_keys,
659 param.list_float_keys, param.matrix_keys]
660 parameters = {}
661 for scaned_datum in scaned_data:
662 for scaned_key in scaned_datum.keys():
663 for openmx_keyword in openmx_keywords:
664 if scaned_key in get_standard_key(openmx_keyword):
665 parameters[scaned_key] = scaned_datum[scaned_key]
666 continue
667 translated_parameters = get_standard_parameters(parameters)
668 parameters.update(translated_parameters)
669 return {k: v for k, v in parameters.items() if v is not None}
672def get_standard_key(key):
673 """
674 Standard ASE parameter format is to USE unerbar(_) instead of dot(.). Also,
675 It is recommended to use lower case alphabet letter. Not Upper. Thus, we
676 change the key to standard key
677 For example:
678 'scf.XcType' -> 'scf_xctype'
679 """
680 if isinstance(key, str):
681 return key.lower().replace('.', '_')
682 elif isinstance(key, list):
683 return [k.lower().replace('.', '_') for k in key]
684 else:
685 return [k.lower().replace('.', '_') for k in key]
688def get_standard_parameters(parameters):
689 """
690 Translate the OpenMX parameters to standard ASE parameters. For example,
692 scf.XcType -> xc
693 scf.maxIter -> maxiter
694 scf.energycutoff -> energy_cutoff
695 scf.Kgrid -> kpts
696 scf.EigenvalueSolver -> eigensolver
697 scf.SpinPolarization -> spinpol
698 scf.criterion -> convergence
699 scf.Electric.Field -> external
700 scf.Mixing.Type -> mixer
701 scf.system.charge -> charge
703 We followed GPAW schem.
704 """
705 from ase.calculators.openmx import parameters as param
706 from ase.units import Bohr, Ha, Ry, fs, m, s
707 units = param.unit_dat_keywords
708 standard_parameters = {}
709 standard_units = {'eV': 1, 'Ha': Ha, 'Ry': Ry, 'Bohr': Bohr, 'fs': fs,
710 'K': 1, 'GV / m': 1e9/1.6e-19 / m, 'Ha/Bohr': Ha/Bohr,
711 'm/s': m/s, '_amu': 1, 'Tesla': 1}
712 translated_parameters = {
713 'scf.XcType': 'xc',
714 'scf.maxIter': 'maxiter',
715 'scf.energycutoff': 'energy_cutoff',
716 'scf.Kgrid': 'kpts',
717 'scf.EigenvalueSolver': 'eigensolver',
718 'scf.SpinPolarization': 'spinpol',
719 'scf.criterion': 'convergence',
720 'scf.Electric.Field': 'external',
721 'scf.Mixing.Type': 'mixer',
722 'scf.system.charge': 'charge'
723 }
725 for key in parameters.keys():
726 for openmx_key in translated_parameters.keys():
727 if key == get_standard_key(openmx_key):
728 standard_key = translated_parameters[openmx_key]
729 unit = standard_units.get(units.get(openmx_key), 1)
730 standard_parameters[standard_key] = parameters[key] * unit
731 standard_parameters['spinpol'] = parameters.get('scf_spinpolarization')
732 return standard_parameters
735def get_atomic_formula(out_data=None, log_data=None, restart_data=None,
736 scfout_data=None, dat_data=None):
737 """_formula'.
738 OpenMX results gives following information. Since, we should pick one
739 between position/scaled_position, scaled_positions are suppressed by
740 default. We use input value of position. Not the position after
741 calculation. It is temporal.
743 Atoms.SpeciesAndCoordinate -> symbols
744 Atoms.SpeciesAndCoordinate -> positions
745 Atoms.UnitVectors -> cell
746 scaled_positions -> scaled_positions
747 If `positions` and `scaled_positions` are both given, this key deleted
748 magmoms -> magmoms, Single value for each atom or three numbers for each
749 atom for non-collinear calculations.
750 """
751 atomic_formula = {}
752 parameters = {'symbols': list, 'positions': list, 'scaled_positions': list,
753 'magmoms': list, 'cell': list}
754 datas = [out_data, log_data, restart_data, scfout_data, dat_data]
755 atoms_unitvectors = None
756 atoms_spncrd_unit = 'ang'
757 atoms_unitvectors_unit = 'ang'
758 for data in datas:
759 # positions unit save
760 if 'atoms_speciesandcoordinates_unit' in data:
761 atoms_spncrd_unit = data['atoms_speciesandcoordinates_unit']
762 # cell unit save
763 if 'atoms_unitvectors_unit' in data:
764 atoms_unitvectors_unit = data['atoms_unitvectors_unit']
765 # symbols, positions or scaled_positions
766 if 'atoms_speciesandcoordinates' in data:
767 atoms_spncrd = data['atoms_speciesandcoordinates']
768 # cell
769 if 'atoms_unitvectors' in data:
770 atoms_unitvectors = data['atoms_unitvectors']
771 # pbc
772 if 'scf_eigenvaluesolver' in data:
773 scf_eigenvaluesolver = data['scf_eigenvaluesolver']
774 # ???
775 for openmx_keyword in data.keys():
776 for standard_keyword in parameters.keys():
777 if openmx_keyword == standard_keyword:
778 atomic_formula[standard_keyword] = data[openmx_keyword]
780 atomic_formula['symbols'] = [i[1] for i in atoms_spncrd]
782 openmx_spncrd_keyword = [[i[2], i[3], i[4]] for i in atoms_spncrd]
783 # Positions
784 positions_unit = atoms_spncrd_unit.lower()
785 positions = np.array(openmx_spncrd_keyword, dtype=float)
786 if positions_unit == 'ang':
787 atomic_formula['positions'] = positions
788 elif positions_unit == 'frac':
789 scaled_positions = np.array(openmx_spncrd_keyword, dtype=float)
790 atomic_formula['scaled_positions'] = scaled_positions
791 elif positions_unit == 'au':
792 positions = np.array(openmx_spncrd_keyword, dtype=float) * Bohr
793 atomic_formula['positions'] = positions
795 # If Cluster, pbc is False, else it is True
796 atomic_formula['pbc'] = scf_eigenvaluesolver.lower() != 'cluster'
798 # Cell Handling
799 if atoms_unitvectors is not None:
800 openmx_cell_keyword = atoms_unitvectors
801 cell = np.array(openmx_cell_keyword, dtype=float)
802 if atoms_unitvectors_unit.lower() == 'ang':
803 atomic_formula['cell'] = openmx_cell_keyword
804 elif atoms_unitvectors_unit.lower() == 'au':
805 atomic_formula['cell'] = cell * Bohr
807 # If `positions` and `scaled_positions` are both given, delete `scaled_..`
808 if atomic_formula.get('scaled_positions') is not None and \
809 atomic_formula.get('positions') is not None:
810 del atomic_formula['scaled_positions']
811 return atomic_formula
814def get_results(out_data=None, log_data=None, restart_data=None,
815 scfout_data=None, dat_data=None, band_data=None):
816 """
817 From the gien data sets, construct the dictionary 'results' and return it'
818 OpenMX version 3.8 can yield following properties
819 free_energy, Ha # Same value with energy
820 energy, Ha
821 energies, Ha
822 forces, Ha/Bohr
823 stress(after 3.8 only) Ha/Bohr**3
824 dipole Debye
825 read_chemical_potential Ha
826 magmoms muB ?? set to 1
827 magmom muB ?? set to 1
828 """
829 from numpy import array as arr
830 results = {}
831 implemented_properties = {'free_energy': Ha, 'energy': Ha, 'energies': Ha,
832 'forces': Ha/Bohr, 'stress': Ha/Bohr**3,
833 'dipole': Debye, 'chemical_potential': Ha,
834 'magmom': 1, 'magmoms': 1, 'eigenvalues': Ha}
835 data = [out_data, log_data, restart_data, scfout_data, dat_data, band_data]
836 for datum in data:
837 for key in datum.keys():
838 for property in implemented_properties.keys():
839 if key == property:
840 results[key] = arr(datum[key])*implemented_properties[key]
841 return results
844def get_file_name(extension='.out', filename=None, absolute_directory=True):
845 directory, prefix = os.path.split(filename)
846 if directory == '':
847 directory = os.curdir
848 if absolute_directory:
849 return os.path.abspath(directory + '/' + prefix + extension)
850 else:
851 return os.path.basename(directory + '/' + prefix + extension)