Coverage for gpaw/defects/__init__.py: 11%
231 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1import numpy as np
2import numbers
3from scipy.optimize import minimize
4from scipy.integrate import simpson
5from gpaw import GPAW, PW
6from ase.parallel import parprint
7import scipy.integrate as integrate
8from scipy.interpolate import InterpolatedUnivariateSpline
9from ase.units import Hartree as Ha
10from ase.units import Bohr
13class ElectrostaticCorrections():
14 """
15 Calculate the electrostatic corrections for charged defects.
16 """
17 def __init__(self, pristine, charged,
18 q=None, sigma=None, r0=None, dimensionality='3d'):
19 if isinstance(pristine, str):
20 pristine = GPAW(pristine, txt=None, parallel={'domain': 1})
21 if isinstance(charged, str):
22 charged = GPAW(charged, txt=None)
23 calc = GPAW(mode=PW(500, force_complex_dtype=True),
24 kpts={'size': (1, 1, 1),
25 'gamma': True},
26 parallel={'domain': 1},
27 symmetry='off',
28 txt=None)
29 atoms = pristine.atoms.copy()
30 calc.initialize(atoms)
32 self.pristine = pristine
33 self.charged = charged
34 self.calc = calc
35 self.q = q
36 self.sigma = sigma
37 self.dimensionality = dimensionality
39 self.pd = self.calc.wfs.pd
40 self.L = self.pd.gd.cell_cv[2, 2]
41 if r0 is not None:
42 self.r0 = np.array(r0) / Bohr
43 elif dimensionality == '2d':
44 self.r0 = np.array([0, 0, self.L / 2])
45 else:
46 self.r0 = np.array([0, 0, 0])
47 self.z0 = -self.r0[2]
48 self.G_Gv = self.pd.get_reciprocal_vectors(q=0, add_q=False)
49 self.G2_G = self.pd.G2_qG[0] # |\vec{G}|^2 in Bohr^-2
50 self.rho_G = self.calculate_gaussian_density()
51 self.Omega = abs(np.linalg.det(self.calc.density.gd.cell_cv))
52 self.data = None
53 self.El = None
55 # For the 2D case, we assume that the dielectric profile epsilon(z)
56 # follows the electronic density of the pristine system n(z).
57 density = self.pristine.get_pseudo_density(gridrefinement=2)
58 density_1d = density.mean(axis=(0, 1))
59 coarse_z = find_z(self.pristine.density.gd.refine())
60 fine_z = find_z(self.pd.gd.refine())
61 transformer = InterpolatedUnivariateSpline(coarse_z, density_1d)
62 self.density_1d = np.array([transformer(x) for x in fine_z])
64 # We potentially have two z axes -- one on the original calculation,
65 # and one on the calculator we just set up.
66 self.z_g = fine_z
67 self.density_z = coarse_z
69 self.G_z = find_G_z(self.G_Gv)
70 self.G_parallel = np.unique(self.G_Gv[:, :2], axis=0)
71 self.GG = np.outer(self.G_z, self.G_z) # G * Gprime
73 # We need the G_z vectors on the finely sampled grid, to calculate
74 # epsilon(G - G'), which goes up to 2G_max in magnitude.
75 self.G_z_fine = (np.fft.fftfreq(len(fine_z)) *
76 2 * np.pi / (self.L) * len(fine_z))
77 self.index_array = self.get_index_array()
79 self.epsilons_z = None
80 self.Eli = None
81 self.Elp = None
82 self.epsilon_GG = None
84 def calculate_gaussian_density(self):
85 # Fourier transformed gaussian:
86 return self.q * np.exp(-0.5 * self.G2_G * self.sigma ** 2)
88 def set_epsilons(self, epsilons, epsilon_bulk=1):
89 """Set the bulk dielectric constant of the system corresponding to the
90 atomic configuration of the pristine system. This is used to calculate
91 the screened coulomb potential of the gaussian model distribution that
92 we use. In 2D, the dielectric constant is poorly defined. However, when
93 we do DFT calculations with periodic systems, we are always working
94 with pseudo bulk-like structures, due to the periodic boundary
95 conditions. It is therefore possible to calculate some dielectric
96 constant. The screening of the system must come from where there is an
97 electronic density, and we thus set
99 epsilon(z) = k * n(z) + epsilon_bulk,
101 for some constant of proprotionality k. This is determined by the
102 constraint that the total screening int dz epsilon(z) gives the correct
103 value.
105 For 3d systems, the trick is to set k = 0. In that case, we have the
106 constant dielectric function that we require.
108 Parameters:
109 epsilons: A float, or a list of two floats. If a single
110 float is given, the screening in-plane and out-of-plane is assumed
111 to be the same. If two floats are given, the first represents the
112 in-plane response, and the second represents the out-of-plane
113 response.
114 epsilon_bulk: The value of the screening far away from the system,
115 given in the same format as epsilons. In most cases, this will be
116 1, corresponding to vacuum.
118 Returns:
119 epsilons_z: dictionary containnig the normalized dielectric
120 profiles in the in-plane and out-of-plane directions along the z
121 axis.
122 """
124 # Reset the state of the object; the corrections need to be
125 # recalculated
126 self.Elp = None
127 self.Eli = None
128 self.epsilon_GG = None
130 def normalize(value):
131 if isinstance(value, numbers.Number):
132 value = [value, value]
133 elif len(value) == 1:
134 value = [value[0]] * 2
135 return np.array(value)
137 epsilons = normalize(epsilons)
138 self.epsilons = epsilons
140 if self.dimensionality == '2d':
141 eb = normalize(epsilon_bulk)
142 elif self.dimensionality == '3d':
143 eb = normalize(epsilons)
144 self.eb = eb
146 z = self.z_g
147 L = self.L
148 density_1d = self.density_1d
149 N = simpson(density_1d, z)
150 epsilons_z = np.zeros((2,) + np.shape(density_1d))
151 epsilons_z += density_1d
153 # In-plane
154 epsilons_z[0] = epsilons_z[0] * (epsilons[0] - eb[0]) / N * L + eb[0]
156 # Out-of-plane
157 def objective_function(k):
158 k = k[0]
159 integral = simpson(1 / (k * density_1d + eb[1]), z) / L
160 return np.abs(integral - 1 / epsilons[1])
162 test = minimize(objective_function, [1], method='Nelder-Mead')
163 assert test.success, "Unable to normalize dielectric profile"
164 k = test.x[0]
165 epsilons_z[1] = epsilons_z[1] * k + eb[1]
166 self.epsilons_z = {'in-plane': epsilons_z[0],
167 'out-of-plane': epsilons_z[1]}
168 self.calculate_epsilon_GG()
170 def get_index_array(self):
171 """
172 Calculate the indices that map between the G vectors on the fine grid,
173 and the matrix defined by M(G, Gprime) = G - Gprime.
175 We use this to generate the matrix epsilon_GG = epsilon(G - Gprime)
176 based on knowledge of epsilon(G) on the fine grid.
177 """
178 G, Gprime = np.meshgrid(self.G_z, self.G_z)
179 difference_GG = G - Gprime
180 index_array = np.zeros(np.shape(difference_GG))
181 index_array[:] = np.nan
182 for idx, val in enumerate(self.G_z_fine):
183 mask = np.isclose(difference_GG, val)
184 index_array[mask] = idx
185 if np.isnan(index_array).any():
186 print("Missing indices found in mapping between plane wave "
187 "sets. Something is wrong with your G-vectors!")
188 print(np.isnan(index_array).all())
189 print(np.where(np.isnan(index_array)))
190 print(self.G_z_fine)
191 print(difference_GG)
192 assert False
193 return index_array.astype(int)
195 def calculate_epsilon_GG(self):
196 assert self.epsilons_z is not None, ("You provide the dielectric "
197 "constant first!")
198 N = len(self.density_1d)
199 epsilon_par_G = np.fft.fft(self.epsilons_z['in-plane']) / N
200 epsilon_perp_G = np.fft.fft(self.epsilons_z['out-of-plane']) / N
201 epsilon_par_GG = epsilon_par_G[self.index_array]
202 epsilon_perp_GG = epsilon_perp_G[self.index_array]
204 self.epsilon_GG = {'in-plane': epsilon_par_GG,
205 'out-of-plane': epsilon_perp_GG}
207 def calculate_periodic_correction(self):
208 G_z = self.G_z
209 if self.Elp is not None:
210 return self.Elp
211 if self.epsilon_GG is None:
212 self.calculate_epsilon_GG()
213 Elp = 0.0
215 for vector in self.G_parallel:
216 G2 = (np.dot(vector, vector) + G_z * G_z)
217 rho_G = self.q * np.exp(- G2 * self.sigma ** 2 / 2, dtype=complex)
218 if self.dimensionality == '2d':
219 rho_G *= np.exp(1j * self.z0 * G_z)
220 A_GG = (self.GG * self.epsilon_GG['out-of-plane']
221 + np.dot(vector, vector) * self.epsilon_GG['in-plane'])
222 if np.allclose(vector, 0):
223 A_GG[0, 0] = 1 # The d.c. potential is poorly defined
224 V_G = np.linalg.solve(A_GG, rho_G)
225 if np.allclose(vector, 0):
226 parprint('Skipping G^2=0 contribution to Elp')
227 V_G[0] = 0
228 Elp += (rho_G * V_G).sum().real
230 Elp *= 2.0 * np.pi * Ha / self.Omega
232 self.Elp = Elp
233 return Elp
235 def calculate_isolated_correction(self):
236 if self.Eli is not None:
237 return self.Eli
239 L = self.L
240 G_z = self.G_z
242 G, Gprime = np.meshgrid(G_z, G_z)
243 Delta_GG = G - Gprime
244 phase = Delta_GG * self.z0
245 gaussian = (Gprime ** 2 + G ** 2) * self.sigma * self.sigma / 2
247 prefactor = np.exp(1j * phase - gaussian)
249 if self.epsilon_GG is None:
250 self.calculate_epsilon_GG()
252 dE_GG_par = (self.epsilon_GG['in-plane']
253 - self.eb[0] * np.eye(len(self.G_z)))
254 dE_GG_perp = (self.epsilon_GG['out-of-plane']
255 - self.eb[1] * np.eye(len(self.G_z)))
257 def integrand(k):
258 K_inv_G = ((self.eb[0] * k ** 2 + self.eb[1] * G_z ** 2) /
259 (1 - np.exp(-k * L / 2) * np.cos(L * G_z / 2)))
260 K_inv_GG = np.diag(K_inv_G)
261 D_GG = L * (K_inv_GG + dE_GG_perp * self.GG + dE_GG_par * k ** 2)
262 return (k * np.exp(-k ** 2 * self.sigma ** 2) *
263 (prefactor * np.linalg.inv(D_GG)).sum()).real
265 I = integrate.quad(integrand, 0, np.inf, limit=500)
266 Eli = self.q * self.q * I[0] * Ha
267 self.Eli = Eli
268 return Eli
270 def calculate_potential_alignment(self):
271 V_neutral = -np.mean(self.pristine.get_electrostatic_potential(),
272 (0, 1))
273 V_charged = -np.mean(self.charged.get_electrostatic_potential(),
274 (0, 1))
276 z = self.density_z
277 z_model, V_model = self.calculate_z_avg_model_potential()
278 V_model = InterpolatedUnivariateSpline(z_model[:-1], V_model[:-1])
279 V_model = V_model(z)
280 self.V_model = V_model
281 Delta_V = self.average(V_model - V_charged + V_neutral, z)
282 return Delta_V
284 def calculate_z_avg_model_potential(self):
286 vox3 = self.calc.density.gd.cell_cv[2, :] / len(self.z_g)
288 # The grid is arranged with z increasing fastest, then y
289 # then x (like a cube file)
291 G_z = self.G_z[1:]
292 rho_Gz = self.q * np.exp(-0.5 * G_z * G_z * self.sigma * self.sigma)
294 zs = []
295 Vs = []
296 if self.dimensionality == '2d':
297 phase = np.exp(1j * (self.G_z * self.z0))
298 A_GG = (self.GG * self.epsilon_GG['out-of-plane'])
299 A_GG[0, 0] = 1
300 V_G = np.linalg.solve(A_GG,
301 phase * np.array([0] + list(rho_Gz)))[1:]
302 elif self.dimensionality == '3d':
303 phase = np.exp(1j * (G_z * self.z0))
304 V_G = phase * rho_Gz / self.eb[1] / G_z ** 2
306 for z in self.z_g:
307 phase_G = np.exp(1j * (G_z * z))
308 V = (np.sum(phase_G * V_G).real
309 * Ha * 4.0 * np.pi / (self.Omega))
310 Vs.append(V)
312 V = (np.sum(V_G.real) * Ha * 4.0 * np.pi / (self.Omega))
313 zs = list(self.z_g) + [vox3[2]]
314 Vs.append(V)
315 return np.array(zs), np.array(Vs)
317 def average(self, V, z):
318 N = len(V)
319 if self.dimensionality == '3d':
320 middle = np.argmin(np.abs(z + self.z0)) + N // 2
321 middle = middle % len(z)
322 points = range(middle - N // 8, middle + N // 8 + 1)
323 restricted = V[points]
324 elif self.dimensionality == '2d':
325 points = list(range(0, N // 8)) + list(range(7 * N // 8, N))
326 restricted = V[points]
327 V_mean = np.mean(restricted)
328 return V_mean
330 def calculate_corrected_formation_energy(self):
331 E_0 = self.pristine.get_potential_energy()
332 E_X = self.charged.get_potential_energy()
333 Eli = self.calculate_isolated_correction()
334 Elp = self.calculate_periodic_correction()
335 Delta_V = self.calculate_potential_alignment()
336 return E_X - E_0 - (Elp - Eli) + Delta_V * self.q
338 def calculate_uncorrected_formation_energy(self):
339 E_0 = self.pristine.get_potential_energy()
340 E_X = self.charged.get_potential_energy()
341 return E_X - E_0
343 def collect_electrostatic_data(self):
344 V_neutral = -np.mean(self.pristine.get_electrostatic_potential(),
345 (0, 1)),
346 V_charged = -np.mean(self.charged.get_electrostatic_potential(),
347 (0, 1)),
348 data = {'epsilon': self.eb[0],
349 'z': self.density_z,
350 'V_0': V_neutral,
351 'V_X': V_charged,
352 'Elc': self.Elp - self.Eli,
353 'D_V_mean': self.calculate_potential_alignment(),
354 'V_model': self.V_model,
355 'D_V': (self.V_model
356 + V_neutral
357 - V_charged)}
358 self.data = data
359 return data
362def find_G_z(G_Gv):
363 mask = (G_Gv[:, 0] == 0) & (G_Gv[:, 1] == 0)
364 G_z = G_Gv[mask][:, 2] # qG_z vectors in Bohr^{-1}
365 return G_z
368def find_z(gd):
369 r3_xyz = gd.get_grid_point_coordinates()
370 nrz = r3_xyz.shape[3]
371 return r3_xyz[2].flatten()[:nrz]