Coverage for gpaw/io/old.py: 9%
148 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-09 00:21 +0000
1"""Code for reading old gpw files."""
2import warnings
4import numpy as np
5import ase.io.ulm as ulm
6from ase import Atoms
7from ase.dft.kpoints import monkhorst_pack
8from ase.io.trajectory import write_atoms
9from ase.units import AUT, Bohr, Ha
10from gpaw.utilities import devnull
12from gpaw import PW
13from gpaw.occupations import FermiDirac
14from gpaw.io.tar import Reader
17def wrap_old_gpw_reader(filename):
18 warnings.warn('You are reading an old-style gpw-file. Please check '
19 'the results carefully!')
21 r = Reader(filename)
23 data = {'version': -1,
24 'gpaw_version': '1.0',
25 'ha': Ha,
26 'bohr': Bohr,
27 'scf.': {'converged': True},
28 'atoms.': {},
29 'wave_functions.': {}}
31 class DictBackend:
32 def write(self, **kwargs):
33 data['atoms.'].update(kwargs)
35 write_atoms(DictBackend(), read_atoms(r))
37 e_total_extrapolated = r.get('PotentialEnergy').item() * Ha
38 magmom_a = r.get('MagneticMoments')
39 data['results.'] = {
40 'energy': e_total_extrapolated,
41 'magmoms': magmom_a,
42 'magmom': magmom_a.sum()}
44 if r.has_array('CartesianForces'):
45 data['results.']['forces'] = r.get('CartesianForces') * Ha / Bohr
47 p = data['parameters.'] = {}
49 p['xc'] = r['XCFunctional']
50 p['nbands'] = r.dimension('nbands')
51 p['spinpol'] = (r.dimension('nspins') == 2)
53 bzk_kc = r.get('BZKPoints', broadcast=True)
54 if r.has_array('NBZKPoints'):
55 p['kpts'] = r.get('NBZKPoints', broadcast=True)
56 if r.has_array('MonkhorstPackOffset'):
57 offset_c = r.get('MonkhorstPackOffset', broadcast=True)
58 if offset_c.any():
59 p['kpts'] = monkhorst_pack(p['kpts']) + offset_c
60 else:
61 p['kpts'] = bzk_kc
63 if r['version'] < 4:
64 usesymm = r['UseSymmetry']
65 if usesymm is None:
66 p['symmetry'] = {'time_reversal': False, 'point_group': False}
67 elif usesymm:
68 p['symmetry'] = {'time_reversal': True, 'point_group': True}
69 else:
70 p['symmetry'] = {'time_reversal': True, 'point_group': False}
71 else:
72 p['symmetry'] = {'point_group': r['SymmetryOnSwitch'],
73 'symmorphic': r['SymmetrySymmorphicSwitch'],
74 'time_reversal': r['SymmetryTimeReversalSwitch'],
75 'tolerance': r['SymmetryToleranceCriterion']}
77 p['basis'] = r['BasisSet']
79 try:
80 h = r['GridSpacing']
81 except KeyError: # CMR can't handle None!
82 h = None
83 if h is not None:
84 p['h'] = Bohr * h
85 if r.has_array('GridPoints'):
86 p['gpts'] = r.get('GridPoints')
88 p['lmax'] = r['MaximumAngularMomentum']
89 p['setups'] = r['SetupTypes']
90 p['fixdensity'] = r['FixDensity']
91 nbtc = r['NumberOfBandsToConverge']
92 if not isinstance(nbtc, (int, str)):
93 # The string 'all' was eval'ed to the all() function!
94 nbtc = 'all'
95 p['convergence'] = {'density': r['DensityConvergenceCriterion'],
96 'energy': r['EnergyConvergenceCriterion'] * Ha,
97 'eigenstates': r['EigenstatesConvergenceCriterion'],
98 'bands': nbtc}
99 mixer = r['MixClass']
100 weight = r['MixWeight']
102 for key in ['basis', 'setups']:
103 dct = p[key]
104 if isinstance(dct, dict) and None in dct:
105 dct['default'] = dct.pop(None)
107 if mixer == 'Mixer':
108 from gpaw.mixer import Mixer
109 elif mixer == 'MixerSum':
110 from gpaw.mixer import MixerSum as Mixer
111 elif mixer == 'MixerSum2':
112 from gpaw.mixer import MixerSum2 as Mixer
113 elif mixer == 'MixerDif':
114 from gpaw.mixer import MixerDif as Mixer
115 elif mixer == 'DummyMixer':
116 from gpaw.mixer import DummyMixer as Mixer
117 else:
118 Mixer = None
120 if Mixer is None:
121 p['mixer'] = None
122 else:
123 p['mixer'] = Mixer(r['MixBeta'], r['MixOld'], weight)
125 p['stencils'] = (r['KohnShamStencil'],
126 r['InterpolationStencil'])
128 vt_sG = r.get('PseudoPotential') * Ha
129 ps = r['PoissonStencil']
130 if isinstance(ps, int) or ps == 'M':
131 poisson = {'name': 'fd'}
132 poisson['nn'] = ps
133 if data['atoms.']['pbc'] == [1, 1, 0]:
134 v1, v2 = vt_sG[0, :, :, [0, -1]].mean(axis=(1, 2))
135 if abs(v1 - v2) > 0.01:
136 warnings.warn('I am guessing that this calculation was done '
137 'with a dipole-layer correction?')
138 poisson['dipolelayer'] = 'xy'
139 p['poissonsolver'] = poisson
141 p['charge'] = r['Charge']
142 fixmom = r['FixMagneticMoment']
144 p['occupations'] = FermiDirac(r['FermiWidth'] * Ha,
145 fixmagmom=fixmom)
147 p['mode'] = r['Mode']
149 if p['mode'] == 'pw':
150 p['mode'] = PW(ecut=r['PlaneWaveCutoff'] * Ha)
152 if len(bzk_kc) == 1 and not bzk_kc[0].any():
153 # Gamma point only:
154 if r['DataType'] == 'Complex':
155 p['dtype'] = complex
157 data['occupations.'] = {
158 'fermilevel': r['FermiLevel'] * Ha,
159 'split': r.parameters.get('FermiSplit', 0) * Ha,
160 'homo': np.nan,
161 'lumo': np.nan}
163 data['density.'] = {
164 'density': r.get('PseudoElectronDensity') * Bohr**-3,
165 'atomic_density_matrices': r.get('AtomicDensityMatrices')}
167 data['hamiltonian.'] = {
168 'e_coulomb': r['Epot'] * Ha,
169 'e_entropy': -r['S'] * Ha,
170 'e_external': r['Eext'] * Ha,
171 'e_kinetic': r['Ekin'] * Ha,
172 'e_total_extrapolated': e_total_extrapolated,
173 'e_xc': r['Exc'] * Ha,
174 'e_zero': r['Ebar'] * Ha,
175 'potential': vt_sG,
176 'atomic_hamiltonian_matrices': r.get('NonLocalPartOfHamiltonian') * Ha}
177 data['hamiltonian.']['e_total_free'] = (
178 sum(data['hamiltonian.'][e] for e in ['e_coulomb', 'e_entropy',
179 'e_external', 'e_kinetic',
180 'e_xc', 'e_zero']))
182 if r.has_array('GLLBPseudoResponsePotential'):
183 data['hamiltonian.']['xc.'] = {
184 'gllb_pseudo_response_potential':
185 r.get('GLLBPseudoResponsePotential') * Ha,
186 'gllb_dxc_pseudo_response_potential':
187 r.get('GLLBDxcPseudoResponsePotential') * Ha / Bohr,
188 'gllb_atomic_density_matrices':
189 r.get('GLLBAtomicDensityMatrices'),
190 'gllb_atomic_response_matrices':
191 r.get('GLLBAtomicResponseMatrices'),
192 'gllb_dxc_atomic_density_matrices':
193 r.get('GLLBDxcAtomicDensityMatrices'),
194 'gllb_dxc_atomic_response_matrices':
195 r.get('GLLBDxcAtomicResponseMatrices')}
197 special = [('eigenvalues', 'Eigenvalues'),
198 ('occupations', 'OccupationNumbers'),
199 ('projections', 'Projections')]
201 if r['Mode'] == 'pw':
202 special.append(('coefficients', 'PseudoWaveFunctions'))
203 try:
204 data['wave_functions.']['indices'] = r.get('PlaneWaveIndices')
205 except KeyError:
206 pass
207 elif r['Mode'] == 'fd':
208 special.append(('values', 'PseudoWaveFunctions'))
209 else:
210 special.append(('coefficients', 'WaveFunctionCoefficients'))
212 for name, old in special:
213 try:
214 fd, shape, size, dtype = r.get_file_object(old, ())
215 except KeyError:
216 continue
217 offset = fd
218 data['wave_functions.'][name + '.'] = {
219 'ndarray': (shape, dtype.name, offset)}
221 new = ulm.Reader(devnull, data=data,
222 little_endian=r.byteswap ^ np.little_endian)
224 for ref in new._data['wave_functions']._data.values():
225 try:
226 ref.fd = ref.offset
227 except AttributeError:
228 continue
229 ref.offset = 0
231 return new
234def read_atoms(reader):
235 positions = reader.get('CartesianPositions', broadcast=True) * Bohr
236 numbers = reader.get('AtomicNumbers', broadcast=True)
237 cell = reader.get('UnitCell', broadcast=True) * Bohr
238 pbc = reader.get('BoundaryConditions', broadcast=True)
239 tags = reader.get('Tags', broadcast=True)
240 magmoms = reader.get('MagneticMoments', broadcast=True)
242 # Create instance of Atoms object, and set_tags and magnetic moments
243 atoms = Atoms(positions=positions,
244 numbers=numbers,
245 cell=cell,
246 pbc=pbc)
248 if tags.any():
249 atoms.set_tags(tags)
250 if magmoms.any():
251 atoms.set_initial_magnetic_moments(magmoms)
253 if reader.has_array('CartesianVelocities'):
254 velocities = reader.get('CartesianVelocities',
255 broadcast=True) * Bohr / AUT
256 atoms.set_velocities(velocities)
258 return atoms