Coverage for /builds/ase/ase/ase/calculators/openmx/reader.py : 62.91%
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
29def read_openmx(filename=None, debug=False):
30 from ase.calculators.openmx import OpenMX
31 from ase import Atoms
32 """
33 Read results from typical OpenMX output files and returns the atom object
34 In default mode, it reads every implementd properties we could get from
35 the files. Unlike previous version, we read the information based on file.
36 previous results will be eraised unless previous results are written in the
37 next calculation results.
39 Read the 'LABEL.log' file seems redundant. Because all the
40 information should already be written in '.out' file. However, in the
41 version 3.8.3, stress tensor are not written in the '.out' file. It only
42 contained in the '.log' file. So... I implented reading '.log' file method
43 """
44 log_data = read_file(get_file_name('.log', filename), debug=debug)
45 restart_data = read_file(get_file_name('.dat#', filename), debug=debug)
46 dat_data = read_file(get_file_name('.dat', filename), debug=debug)
47 out_data = read_file(get_file_name('.out', filename), debug=debug)
48 scfout_data = read_scfout_file(get_file_name('.scfout', filename))
49 band_data = read_band_file(get_file_name('.Band', filename))
50 # dos_data = read_dos_file(get_file_name('.Dos.val', filename))
51 """
52 First, we get every data we could get from the all results files. And then,
53 reform the data to fit to data structure of Atom object. While doing this,
54 Fix the unit to ASE format.
55 """
56 parameters = get_parameters(out_data=out_data, log_data=log_data,
57 restart_data=restart_data, dat_data=dat_data,
58 scfout_data=scfout_data, band_data=band_data)
59 atomic_formula = get_atomic_formula(out_data=out_data, log_data=log_data,
60 restart_data=restart_data,
61 scfout_data=scfout_data,
62 dat_data=dat_data)
63 results = get_results(out_data=out_data, log_data=log_data,
64 restart_data=restart_data, scfout_data=scfout_data,
65 dat_data=dat_data, band_data=band_data)
67 atoms = Atoms(**atomic_formula)
68 atoms.calc = OpenMX(**parameters)
69 atoms.calc.results = results
70 return atoms
73def read_file(filename, debug=False):
74 """
75 Read the 'LABEL.out' file. Using 'parameters.py', we read every 'allowed_
76 dat' dictionory. while reading a file, if one find the key matcheds That
77 'patters', which indicates the property we want is written, it will returns
78 the pair value of that key. For example,
79 example will be written later
80 """
81 from ase.calculators.openmx import parameters as param
82 if not os.path.isfile(filename):
83 return {}
84 param_keys = ['integer_keys', 'float_keys', 'string_keys', 'bool_keys',
85 'list_int_keys', 'list_float_keys', 'list_bool_keys',
86 'tuple_integer_keys', 'tuple_float_keys', 'tuple_float_keys']
87 patterns = {
88 'Stress tensor': ('stress', read_stress_tensor),
89 'Dipole moment': ('dipole', read_dipole),
90 'Fractional coordinates of': ('scaled_positions', read_scaled_positions),
91 'Utot.': ('energy', read_energy),
92 'energies in': ('energies', read_energies),
93 'Chemical Potential': ('chemical_potential', read_chemical_potential),
94 '<coordinates.forces': ('forces', read_forces),
95 'Eigenvalues (Hartree)': ('eigenvalues', read_eigenvalues)}
96 special_patterns = {
97 'Total spin moment': (('magmoms', 'total_magmom'),
98 read_magmoms_and_total_magmom),
99 }
100 out_data = {}
101 line = '\n'
102 if(debug):
103 print('Read results from %s' % filename)
104 with open(filename, 'r') as fd:
105 '''
106 Read output file line by line. When the `line` matches the pattern
107 of certain keywords in `param.[dtype]_keys`, for example,
109 if line in param.string_keys:
110 out_data[key] = read_string(line)
112 parse that line and store it to `out_data` in specified data type.
113 To cover all `dtype` parameters, for loop was used,
115 for [dtype] in parameters_keys:
116 if line in param.[dtype]_keys:
117 out_data[key] = read_[dtype](line)
119 After found matched pattern, escape the for loop using `continue`.
120 '''
121 while line != '':
122 pattern_matched = False
123 line = fd.readline()
124 try:
125 _line = line.split()[0]
126 except IndexError:
127 continue
128 for dtype_key in param_keys:
129 dtype = dtype_key.rsplit('_', 1)[0]
130 read_dtype = globals()['read_' + dtype]
131 for key in param.__dict__[dtype_key]:
132 if key in _line:
133 out_data[get_standard_key(key)] = read_dtype(line)
134 pattern_matched = True
135 continue
136 if pattern_matched:
137 continue
139 for key in param.matrix_keys:
140 if '<'+key in line:
141 out_data[get_standard_key(key)] = read_matrix(line, key, fd)
142 pattern_matched = True
143 continue
144 if pattern_matched:
145 continue
146 for key in patterns.keys():
147 if key in line:
148 out_data[patterns[key][0]] = patterns[key][1](
149 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, fd):
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(fd.read(8*TNO2)))
288 return myOLP
290 def readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, fd):
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(fd.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 = fd.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 = fd.readline().split()
366 band_data['band_kpath_unitcell'] = [line[:3], line[3:6], line[6:9]]
367 line = fd.readline().split()
368 band_data['band_nkpath'] = int(line[0])
369 for i in range(band_data['band_nkpath']):
370 line = fd.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 = fd.readline()
380 kpts[j] = np.array(line.split(), dtype=float)[1:]
381 line = fd.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 = fd.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, fd):
461 matrix = []
462 line = fd.readline()
463 while key not in line:
464 matrix.append(line.split())
465 line = fd.readline()
466 return matrix
469def read_stress_tensor(line, fd, debug=None):
470 fd.readline() # passing empty line
471 fd.readline()
472 line = fd.readline()
473 xx, xy, xz = read_tuple_float(line)
474 line = fd.readline()
475 yx, yy, yz = read_tuple_float(line)
476 line = fd.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, fd, debug=None):
483 total_magmom = read_float(line)
484 fd.readline() # Skip empty lines
485 fd.readline()
486 line = fd.readline()
487 magmoms = []
488 while not(line == '' or line.isspace()):
489 magmoms.append(read_float(line))
490 line = fd.readline()
491 return magmoms, total_magmom
494def read_energy(line, fd, debug=None):
495 # It has Hartree unit yet
496 return read_float(line)
499def read_energies(line, fd, debug=None):
500 line = fd.readline()
501 if '***' in line:
502 point = 7 # Version 3.8
503 else:
504 point = 16 # Version 3.9
505 for i in range(point):
506 fd.readline()
507 line = fd.readline()
508 energies = []
509 while not(line == '' or line.isspace()):
510 energies.append(float(line.split()[2]))
511 line = fd.readline()
512 return energies
515def read_eigenvalues(line, fd, debug=False):
516 """
517 Read the Eigenvalues in the `.out` file and returns the eigenvalue
518 First, it assumes system have two spins and start reading until it reaches
519 the end('*****...').
521 eigenvalues[spin][kpoint][nbands]
523 For symmetry reason, `.out` file prints the eigenvalues at the half of the
524 K points. Thus, we have to fill up the rest of the half.
525 However, if the calculation was conducted only on the gamma point, it will
526 raise the 'gamma_flag' as true and it will returns the original samples.
527 """
528 def prind(*line, end='\n'):
529 if debug:
530 print(*line, end=end)
531 prind("Read eigenvalues output")
532 current_line = fd.tell()
533 fd.seek(0) # Seek for the kgrid information
534 while line != '':
535 line = fd.readline().lower()
536 if 'scf.kgrid' in line:
537 break
538 fd.seek(current_line) # Retrun to the original position
540 kgrid = read_tuple_integer(line)
542 if kgrid != ():
543 prind('Non-Gamma point calculation')
544 prind('scf.Kgrid is %d, %d, %d' % kgrid)
545 gamma_flag = False
546 # fd.seek(f.tell()+57)
547 else:
548 prind('Gamma point calculation')
549 gamma_flag = True
550 line = fd.readline()
551 line = fd.readline()
553 eigenvalues = []
554 eigenvalues.append([])
555 eigenvalues.append([]) # Assume two spins
556 i = 0
557 while True:
558 # Go to eigenvalues line
559 while line != '':
560 line = fd.readline()
561 prind(line)
562 ll = line.split()
563 if line.isspace():
564 continue
565 elif len(ll) > 1:
566 if ll[0] == '1':
567 break
568 elif "*****" in line:
569 break
571 # Break if it reaches the end or next parameters
572 if "*****" in line or line == '':
573 break
575 # Check Number and format is valid
576 try:
577 # Suppose to be looks like
578 # 1 -2.33424746491277 -2.33424746917880
579 ll = line.split()
580 # Check if it reaches the end of the file
581 assert line != ''
582 assert len(ll) == 3
583 float(ll[1])
584 float(ll[2])
585 except (AssertionError, ValueError):
586 raise ParseError("Cannot read eigenvalues")
588 # Read eigenvalues
589 eigenvalues[0].append([])
590 eigenvalues[1].append([])
591 while not (line == '' or line.isspace()):
592 eigenvalues[0][i].append(float(rn(line, 2)))
593 eigenvalues[1][i].append(float(rn(line, 1)))
594 line = fd.readline()
595 prind(line, end='')
596 i += 1
597 prind(line)
598 if gamma_flag:
599 return np.asarray(eigenvalues)
600 eigen_half = np.asarray(eigenvalues)
601 prind(eigen_half)
602 # Fill up the half
603 spin, half_kpts, bands = eigen_half.shape
604 even_odd = np.array(kgrid).prod() % 2
605 eigen_values = np.zeros((spin, half_kpts*2-even_odd, bands))
606 for i in range(half_kpts):
607 eigen_values[0, i] = eigen_half[0, i, :]
608 eigen_values[1, i] = eigen_half[1, i, :]
609 eigen_values[0, 2*half_kpts-1-i-even_odd] = eigen_half[0, i, :]
610 eigen_values[1, 2*half_kpts-1-i-even_odd] = eigen_half[1, i, :]
611 return eigen_values
614def read_forces(line, fd, debug=None):
615 # It has Hartree per Bohr unit yet
616 forces = []
617 fd.readline() # Skip Empty line
618 line = fd.readline()
619 while 'coordinates.forces>' not in line:
620 forces.append(read_tuple_float(line))
621 line = fd.readline()
622 return np.array(forces)
625def read_dipole(line, fd, debug=None):
626 dipole = []
627 while 'Total' not in line:
628 line = fd.readline()
629 dipole.append(read_tuple_float(line))
630 return dipole
633def read_scaled_positions(line, fd, debug=None):
634 scaled_positions = []
635 fd.readline() # Skip Empty lines
636 fd.readline()
637 fd.readline()
638 line = fd.readline()
639 while not(line == '' or line.isspace()): # Detect empty line
640 scaled_positions.append(read_tuple_float(line))
641 line = fd.readline()
642 return scaled_positions
645def read_chemical_potential(line, fd, debug=None):
646 return read_float(line)
649def get_parameters(out_data=None, log_data=None, restart_data=None,
650 scfout_data=None, dat_data=None, band_data=None):
651 """
652 From the given data sets, construct the dictionary 'parameters'. If data
653 is in the paramerters, it will save it.
654 """
655 from ase.calculators.openmx import parameters as param
656 scaned_data = [dat_data, out_data, log_data, restart_data, scfout_data,
657 band_data]
658 openmx_keywords = [param.tuple_integer_keys, param.tuple_float_keys,
659 param.tuple_bool_keys, param.integer_keys,
660 param.float_keys, param.string_keys, param.bool_keys,
661 param.list_int_keys, param.list_bool_keys,
662 param.list_float_keys, param.matrix_keys]
663 parameters = {}
664 for scaned_datum in scaned_data:
665 for scaned_key in scaned_datum.keys():
666 for openmx_keyword in openmx_keywords:
667 if scaned_key in get_standard_key(openmx_keyword):
668 parameters[scaned_key] = scaned_datum[scaned_key]
669 continue
670 translated_parameters = get_standard_parameters(parameters)
671 parameters.update(translated_parameters)
672 return {k: v for k, v in parameters.items() if v is not None}
675def get_standard_key(key):
676 """
677 Standard ASE parameter format is to USE unerbar(_) instead of dot(.). Also,
678 It is recommended to use lower case alphabet letter. Not Upper. Thus, we
679 change the key to standard key
680 For example:
681 'scf.XcType' -> 'scf_xctype'
682 """
683 if isinstance(key, str):
684 return key.lower().replace('.', '_')
685 elif isinstance(key, list):
686 return [k.lower().replace('.', '_') for k in key]
687 else:
688 return [k.lower().replace('.', '_') for k in key]
691def get_standard_parameters(parameters):
692 """
693 Translate the OpenMX parameters to standard ASE parameters. For example,
695 scf.XcType -> xc
696 scf.maxIter -> maxiter
697 scf.energycutoff -> energy_cutoff
698 scf.Kgrid -> kpts
699 scf.EigenvalueSolver -> eigensolver
700 scf.SpinPolarization -> spinpol
701 scf.criterion -> convergence
702 scf.Electric.Field -> external
703 scf.Mixing.Type -> mixer
704 scf.system.charge -> charge
706 We followed GPAW schem.
707 """
708 from ase.calculators.openmx import parameters as param
709 from ase.units import Bohr, Ha, Ry, fs, m, s
710 units = param.unit_dat_keywords
711 standard_parameters = {}
712 standard_units = {'eV': 1, 'Ha': Ha, 'Ry': Ry, 'Bohr': Bohr, 'fs': fs,
713 'K': 1, 'GV / m': 1e9/1.6e-19 / m, 'Ha/Bohr': Ha/Bohr,
714 'm/s': m/s, '_amu': 1, 'Tesla': 1}
715 translated_parameters = {
716 'scf.XcType': 'xc',
717 'scf.maxIter': 'maxiter',
718 'scf.energycutoff': 'energy_cutoff',
719 'scf.Kgrid': 'kpts',
720 'scf.EigenvalueSolver': 'eigensolver',
721 'scf.SpinPolarization': 'spinpol',
722 'scf.criterion': 'convergence',
723 'scf.Electric.Field': 'external',
724 'scf.Mixing.Type': 'mixer',
725 'scf.system.charge': 'charge'
726 }
728 for key in parameters.keys():
729 for openmx_key in translated_parameters.keys():
730 if key == get_standard_key(openmx_key):
731 standard_key = translated_parameters[openmx_key]
732 unit = standard_units.get(units.get(openmx_key), 1)
733 standard_parameters[standard_key] = parameters[key] * unit
734 standard_parameters['spinpol'] = parameters.get('scf_spinpolarization')
735 return standard_parameters
738def get_atomic_formula(out_data=None, log_data=None, restart_data=None,
739 scfout_data=None, dat_data=None):
740 """_formula'.
741 OpenMX results gives following information. Since, we should pick one
742 between position/scaled_position, scaled_positions are suppressed by
743 default. We use input value of position. Not the position after
744 calculation. It is temporal.
746 Atoms.SpeciesAndCoordinate -> symbols
747 Atoms.SpeciesAndCoordinate -> positions
748 Atoms.UnitVectors -> cell
749 scaled_positions -> scaled_positions
750 If `positions` and `scaled_positions` are both given, this key deleted
751 magmoms -> magmoms, Single value for each atom or three numbers for each
752 atom for non-collinear calculations.
753 """
754 atomic_formula = {}
755 parameters = {'symbols': list, 'positions': list, 'scaled_positions': list,
756 'magmoms': list, 'cell': list}
757 datas = [out_data, log_data, restart_data, scfout_data, dat_data]
758 atoms_unitvectors = None
759 atoms_spncrd_unit = 'ang'
760 atoms_unitvectors_unit = 'ang'
761 for data in datas:
762 # positions unit save
763 if 'atoms_speciesandcoordinates_unit' in data:
764 atoms_spncrd_unit = data['atoms_speciesandcoordinates_unit']
765 # cell unit save
766 if 'atoms_unitvectors_unit' in data:
767 atoms_unitvectors_unit = data['atoms_unitvectors_unit']
768 # symbols, positions or scaled_positions
769 if 'atoms_speciesandcoordinates' in data:
770 atoms_spncrd = data['atoms_speciesandcoordinates']
771 # cell
772 if 'atoms_unitvectors' in data:
773 atoms_unitvectors = data['atoms_unitvectors']
774 # pbc
775 if 'scf_eigenvaluesolver' in data:
776 scf_eigenvaluesolver = data['scf_eigenvaluesolver']
777 # ???
778 for openmx_keyword in data.keys():
779 for standard_keyword in parameters.keys():
780 if openmx_keyword == standard_keyword:
781 atomic_formula[standard_keyword] = data[openmx_keyword]
783 atomic_formula['symbols'] = [i[1] for i in atoms_spncrd]
785 openmx_spncrd_keyword = [[i[2], i[3], i[4]] for i in atoms_spncrd]
786 # Positions
787 positions_unit = atoms_spncrd_unit.lower()
788 positions = np.array(openmx_spncrd_keyword, dtype=float)
789 if positions_unit == 'ang':
790 atomic_formula['positions'] = positions
791 elif positions_unit == 'frac':
792 scaled_positions = np.array(openmx_spncrd_keyword, dtype=float)
793 atomic_formula['scaled_positions'] = scaled_positions
794 elif positions_unit == 'au':
795 positions = np.array(openmx_spncrd_keyword, dtype=float) * Bohr
796 atomic_formula['positions'] = positions
798 # If Cluster, pbc is False, else it is True
799 atomic_formula['pbc'] = scf_eigenvaluesolver.lower() != 'cluster'
801 # Cell Handling
802 if atoms_unitvectors is not None:
803 openmx_cell_keyword = atoms_unitvectors
804 cell = np.array(openmx_cell_keyword, dtype=float)
805 if atoms_unitvectors_unit.lower() == 'ang':
806 atomic_formula['cell'] = openmx_cell_keyword
807 elif atoms_unitvectors_unit.lower() == 'au':
808 atomic_formula['cell'] = cell * Bohr
810 # If `positions` and `scaled_positions` are both given, delete `scaled_..`
811 if atomic_formula.get('scaled_positions') is not None and \
812 atomic_formula.get('positions') is not None:
813 del atomic_formula['scaled_positions']
814 return atomic_formula
817def get_results(out_data=None, log_data=None, restart_data=None,
818 scfout_data=None, dat_data=None, band_data=None):
819 """
820 From the gien data sets, construct the dictionary 'results' and return it'
821 OpenMX version 3.8 can yield following properties
822 free_energy, Ha # Same value with energy
823 energy, Ha
824 energies, Ha
825 forces, Ha/Bohr
826 stress(after 3.8 only) Ha/Bohr**3
827 dipole Debye
828 read_chemical_potential Ha
829 magmoms muB ?? set to 1
830 magmom muB ?? set to 1
831 """
832 from numpy import array as arr
833 results = {}
834 implemented_properties = {'free_energy': Ha, 'energy': Ha, 'energies': Ha,
835 'forces': Ha/Bohr, 'stress': Ha/Bohr**3,
836 'dipole': Debye, 'chemical_potential': Ha,
837 'magmom': 1, 'magmoms': 1, 'eigenvalues': Ha}
838 data = [out_data, log_data, restart_data, scfout_data, dat_data, band_data]
839 for datum in data:
840 for key in datum.keys():
841 for property in implemented_properties.keys():
842 if key == property:
843 results[key] = arr(datum[key])*implemented_properties[key]
844 return results
847def get_file_name(extension='.out', filename=None, absolute_directory=True):
848 directory, prefix = os.path.split(filename)
849 if directory == '':
850 directory = os.curdir
851 if absolute_directory:
852 return os.path.abspath(directory + '/' + prefix + extension)
853 else:
854 return os.path.basename(directory + '/' + prefix + extension)