Coverage for /builds/kinetik161/ase/ase/calculators/singlepoint.py: 79.61%
206 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1import numpy as np
3from ase.calculators.calculator import (Calculator,
4 PropertyNotImplementedError,
5 PropertyNotPresent, all_properties)
6from ase.outputs import Properties
7from ase.utils import lazyproperty
10class SinglePointCalculator(Calculator):
11 """Special calculator for a single configuration.
13 Used to remember the energy, force and stress for a given
14 configuration. If the positions, atomic numbers, unit cell, or
15 boundary conditions are changed, then asking for
16 energy/forces/stress will raise an exception."""
18 name = 'unknown'
20 def __init__(self, atoms, **results):
21 """Save energy, forces, stress, ... for the current configuration."""
22 Calculator.__init__(self)
23 self.results = {}
24 for property, value in results.items():
25 assert property in all_properties, property
26 if value is None:
27 continue
28 if property in ['energy', 'magmom', 'free_energy']:
29 self.results[property] = value
30 else:
31 self.results[property] = np.array(value, float)
32 self.atoms = atoms.copy()
34 def __str__(self):
35 tokens = []
36 for key, val in sorted(self.results.items()):
37 if np.isscalar(val):
38 txt = f'{key}={val}'
39 else:
40 txt = f'{key}=...'
41 tokens.append(txt)
42 return '{}({})'.format(self.__class__.__name__, ', '.join(tokens))
44 def get_property(self, name, atoms=None, allow_calculation=True):
45 if atoms is None:
46 atoms = self.atoms
47 if name not in self.results or self.check_state(atoms):
48 if allow_calculation:
49 raise PropertyNotImplementedError(
50 f'The property "{name}" is not available.')
51 return None
53 result = self.results[name]
54 if isinstance(result, np.ndarray):
55 result = result.copy()
56 return result
59class SinglePointKPoint:
60 def __init__(self, weight, s, k, eps_n=None, f_n=None):
61 self.weight = weight
62 self.s = s # spin index
63 self.k = k # k-point index
64 if eps_n is None:
65 eps_n = []
66 self.eps_n = eps_n
67 if f_n is None:
68 f_n = []
69 self.f_n = f_n
72def arrays_to_kpoints(eigenvalues, occupations, weights):
73 """Helper function for building SinglePointKPoints.
75 Convert eigenvalue, occupation, and weight arrays to list of
76 SinglePointKPoint objects."""
77 nspins, nkpts, nbands = eigenvalues.shape
78 assert eigenvalues.shape == occupations.shape
79 assert len(weights) == nkpts
80 kpts = []
81 for s in range(nspins):
82 for k in range(nkpts):
83 kpt = SinglePointKPoint(
84 weight=weights[k], s=s, k=k,
85 eps_n=eigenvalues[s, k], f_n=occupations[s, k])
86 kpts.append(kpt)
87 return kpts
90class SinglePointDFTCalculator(SinglePointCalculator):
91 def __init__(self, atoms,
92 efermi=None, bzkpts=None, ibzkpts=None, bz2ibz=None,
93 kpts=None,
94 **results):
95 self.bz_kpts = bzkpts
96 self.ibz_kpts = ibzkpts
97 self.bz2ibz = bz2ibz
98 self.eFermi = efermi
100 SinglePointCalculator.__init__(self, atoms, **results)
101 self.kpts = kpts
103 def get_fermi_level(self):
104 """Return the Fermi-level(s)."""
105 return self.eFermi
107 def get_bz_to_ibz_map(self):
108 return self.bz2ibz
110 def get_bz_k_points(self):
111 """Return the k-points."""
112 return self.bz_kpts
114 def get_number_of_spins(self):
115 """Return the number of spins in the calculation.
117 Spin-paired calculations: 1, spin-polarized calculation: 2."""
118 if self.kpts is not None:
119 nspin = set()
120 for kpt in self.kpts:
121 nspin.add(kpt.s)
122 return len(nspin)
123 return None
125 def get_number_of_bands(self):
126 values = {len(kpt.eps_n) for kpt in self.kpts}
127 if not values:
128 return None
129 elif len(values) == 1:
130 return values.pop()
131 else:
132 raise RuntimeError('Multiple array sizes')
134 def get_spin_polarized(self):
135 """Is it a spin-polarized calculation?"""
136 nos = self.get_number_of_spins()
137 if nos is not None:
138 return nos == 2
139 return None
141 def get_ibz_k_points(self):
142 """Return k-points in the irreducible part of the Brillouin zone."""
143 return self.ibz_kpts
145 def get_kpt(self, kpt=0, spin=0):
146 if self.kpts is not None:
147 counter = 0
148 for kpoint in self.kpts:
149 if kpoint.s == spin:
150 if kpt == counter:
151 return kpoint
152 counter += 1
153 return None
155 def get_k_point_weights(self):
156 """ Retunrs the weights of the k points """
157 if self.kpts is not None:
158 weights = []
159 for kpoint in self.kpts:
160 if kpoint.s == 0:
161 weights.append(kpoint.weight)
162 return np.array(weights)
163 return None
165 def get_occupation_numbers(self, kpt=0, spin=0):
166 """Return occupation number array."""
167 kpoint = self.get_kpt(kpt, spin)
168 if kpoint is not None:
169 if len(kpoint.f_n):
170 return kpoint.f_n
171 return None
173 def get_eigenvalues(self, kpt=0, spin=0):
174 """Return eigenvalue array."""
175 kpoint = self.get_kpt(kpt, spin)
176 if kpoint is not None:
177 return kpoint.eps_n
178 return None
180 def get_homo_lumo(self):
181 """Return HOMO and LUMO energies."""
182 if self.kpts is None:
183 raise RuntimeError('No kpts')
184 eH = -np.inf
185 eL = np.inf
186 for spin in range(self.get_number_of_spins()):
187 homo, lumo = self.get_homo_lumo_by_spin(spin)
188 eH = max(eH, homo)
189 eL = min(eL, lumo)
190 return eH, eL
192 def get_homo_lumo_by_spin(self, spin=0):
193 """Return HOMO and LUMO energies for a given spin."""
194 if self.kpts is None:
195 raise RuntimeError('No kpts')
196 for kpt in self.kpts:
197 if kpt.s == spin:
198 break
199 else:
200 raise RuntimeError(f'No k-point with spin {spin}')
201 if self.eFermi is None:
202 raise RuntimeError('Fermi level is not available')
203 eH = -1.e32
204 eL = 1.e32
205 for kpt in self.kpts:
206 if kpt.s == spin:
207 for e in kpt.eps_n:
208 if e <= self.eFermi:
209 eH = max(eH, e)
210 else:
211 eL = min(eL, e)
212 return eH, eL
214 def properties(self) -> Properties:
215 return OutputPropertyWrapper(self).properties()
218def propertygetter(func):
219 from functools import wraps
221 @wraps(func)
222 def getter(self):
223 value = func(self)
224 if value is None:
225 raise PropertyNotPresent(func.__name__)
226 return value
227 return lazyproperty(getter)
230class OutputPropertyWrapper:
231 def __init__(self, calc):
232 self.calc = calc
234 @propertygetter
235 def nspins(self):
236 return self.calc.get_number_of_spins()
238 @propertygetter
239 def nbands(self):
240 return self.calc.get_number_of_bands()
242 @propertygetter
243 def nkpts(self):
244 return len(self.calc.kpts) // self.nspins
246 def _build_eig_occ_array(self, getter):
247 arr = np.empty((self.nspins, self.nkpts, self.nbands))
248 for s in range(self.nspins):
249 for k in range(self.nkpts):
250 value = getter(spin=s, kpt=k)
251 if value is None:
252 return None
253 arr[s, k, :] = value
254 return arr
256 @propertygetter
257 def eigenvalues(self):
258 return self._build_eig_occ_array(self.calc.get_eigenvalues)
260 @propertygetter
261 def occupations(self):
262 return self._build_eig_occ_array(self.calc.get_occupation_numbers)
264 @propertygetter
265 def fermi_level(self):
266 return self.calc.get_fermi_level()
268 @propertygetter
269 def kpoint_weights(self):
270 return self.calc.get_k_point_weights()
272 @propertygetter
273 def ibz_kpoints(self):
274 return self.calc.get_ibz_k_points()
276 def properties(self) -> Properties:
277 dct = {}
278 for name in ['eigenvalues', 'occupations', 'fermi_level',
279 'kpoint_weights', 'ibz_kpoints']:
280 try:
281 value = getattr(self, name)
282 except PropertyNotPresent:
283 pass
284 else:
285 dct[name] = value
287 for name, value in self.calc.results.items():
288 dct[name] = value
290 return Properties(dct)