Coverage for gpaw/response/bse.py: 90%
813 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1from dataclasses import dataclass
2from datetime import timedelta
3from functools import cached_property
4from time import time, ctime
6from ase.units import Hartree, Bohr
7from ase.dft import monkhorst_pack
8import numpy as np
9from scipy.linalg import eigh
11from gpaw.blacs import BlacsGrid, Redistributor, BlacsDescriptor
12from gpaw.kpt_descriptor import KPointDescriptor
13from gpaw.mpi import world, serial_comm
14from gpaw.response import ResponseContext
15from gpaw.response.groundstate import CellDescriptor
16from gpaw.response.chi0 import Chi0Calculator, get_frequency_descriptor
17from gpaw.response.context import timer
18from gpaw.response.coulomb_kernels import CoulombKernel
19from gpaw.response.df import Chi0DysonEquations
20from gpaw.response.df import write_response_function
21from gpaw.response.frequencies import FrequencyDescriptor
22from gpaw.response.pair import KPointPairFactory, get_gs_and_context
23from gpaw.response.qpd import SingleQPWDescriptor
24from gpaw.response.screened_interaction import (initialize_w_calculator,
25 GammaIntegrationMode)
26from gpaw.utilities.elpa import LibElpa
27from gpaw.response.pw_parallelization import Blocks1D
30def decide_whether_tammdancoff(val_m, con_m):
31 for n in val_m:
32 if n in con_m:
33 return False
34 return True
37@dataclass
38class BSEMatrix:
39 df_S: np.ndarray
40 H_sS: np.ndarray
41 deps_S: np.ndarray
42 deps_max: float
44 def diagonalize_nontammdancoff(self, bse, deps_max=None):
45 df_S = self.df_S
46 H_sS = self.H_sS
47 if deps_max is None:
48 deps_max = self.deps_max
49 excludef_S = np.where(np.abs(df_S) < 0.001)[0]
50 excludedeps_S = np.where(np.abs(self.deps_S) > deps_max)[0]
51 exclude_S = np.unique(np.concatenate((excludef_S, excludedeps_S)))
52 bse.context.print(' Using numpy.linalg.eig...')
53 bse.context.print(' Eliminated %s pair orbitals' % len(
54 exclude_S))
55 H_SS = bse.collect_A_SS(H_sS)
56 w_T = np.zeros(bse.nS - len(exclude_S), complex)
57 v_ST = None
58 if world.rank == 0:
59 H_SS = np.delete(H_SS, exclude_S, axis=0)
60 H_SS = np.delete(H_SS, exclude_S, axis=1)
61 w_T, v_ST = np.linalg.eig(H_SS)
62 else:
63 v_ST = None
64 world.broadcast(w_T, 0)
65 return w_T, v_ST, exclude_S
67 def diagonalize_tammdancoff(self, bse, deps_max=None, elpa=False):
68 if deps_max is None:
69 deps_max = self.deps_max
70 exclude_S = np.where(np.abs(self.deps_S) > deps_max)[0]
71 H_rr, new_grid_desc = self.exclude_states(bse, exclude_S)
72 if world.size == 1:
73 bse.context.print(' Using lapack...')
74 w_T, v_Rt = eigh(H_rr)
75 return w_T, v_Rt, exclude_S
76 nR = bse.nS - len(exclude_S)
77 w_T = np.empty(nR)
78 v_rt = new_grid_desc.empty(dtype=complex)
79 if elpa:
80 bse.context.print('Using elpa...')
81 elpa = LibElpa(new_grid_desc)
82 elpa.diagonalize(H_rr, v_rt, w_T)
83 else:
84 bse.context.print('Using scalapack...')
85 new_grid_desc.diagonalize_dc(H_rr, v_rt, w_T)
87 # redistribute eigenvectors
88 # we want them to be parallelized over the last index only
90 grid_tR = BlacsGrid(world, world.size, 1)
91 nt = -((-nR) // world.size)
92 desc_tR = grid_tR.new_descriptor(nR, nR, nt, nR)
93 v_tR = desc_tR.zeros(dtype=complex)
94 Redistributor(world, new_grid_desc,
95 desc_tR).redistribute(v_rt, v_tR)
96 v_Rt = v_tR.conj().T
97 return w_T, v_Rt, exclude_S
99 def exclude_states(self, bse, exclude_S):
100 """
101 Removes pairs from the BSE Hamiltonian.
102 A pair is removed if the absolute value of the
103 transition energy eps_c - eps_v is greater than deps_max
104 """
105 H_sS = self.H_sS
106 grid = BlacsGrid(world, world.size, 1)
107 nS = bse.nS
108 ns = bse.ns
109 desc = grid.new_descriptor(nS, nS, ns, nS)
110 H_rr, new_desc = parallel_delete(H_sS, exclude_S, desc)
111 bse.context.print(' Eliminated %s pair orbitals' % len(
112 exclude_S))
113 return H_rr, new_desc
116def parallel_delete(A_nn: np.ndarray,
117 deleteN: np.ndarray,
118 grid_desc: BlacsDescriptor,
119 new_desc=None):
120 """
121 Removes rows and columns from the distributed square matrix A_nn.
122 This is done by redistributing the matrix to first make the second index
123 global (A_nn -> A_mN), and then deleting along the second (global) index;
124 then repeating the same procedure for the first index.
125 --------------
126 Parameters:
127 A_nn: distributed matrix
128 deleteN : list of (global) indices to delete
129 grid_desc: BlacsDescriptor for A_nn
130 new_grid: BlacsGrid on which A_nn will be returned; optional.
131 If None, the output grid will be determined automatically
132 by the gpaw.matrix.suggest_blocking function.
133 -------------------
134 Returns:
135 A_rs: np.ndarray
136 new_desc: BlacsDescriptor
138 A_rs is the (new) distributed matrix after rows and columns
139 have been deleted.
140 new_grid is the grid on which it is distributed.
141 """
142 from gpaw.matrix import suggest_blocking
143 N = grid_desc.N
144 assert N == grid_desc.M, 'Matrix must be square'
146 R = N - len(deleteN) # global matrix dimension after deletion
147 dtype = A_nn.dtype
148 # redistribute matrix, so it is distributed over first index only.
149 # then we can safely delete entries from the second index.
150 grid_mN = BlacsGrid(world, world.size, 1)
151 m = -((-N) // world.size)
152 desc_mN = grid_mN.new_descriptor(N, N, m, N)
153 A_mN = desc_mN.zeros(dtype=dtype)
154 Redistributor(world, grid_desc,
155 desc_mN).redistribute(A_nn, A_mN)
157 # delete, and ensure that array is still contiguous in memory
158 A_mR = np.delete(A_mN, deleteN, axis=1)
159 A_mR = np.ascontiguousarray(A_mR)
160 desc_mR = grid_mN.new_descriptor(N, R, m, R)
162 # now distribute over second index, so we can delete entries from 1st
163 r = -((-R) // world.size)
164 grid_Nr = BlacsGrid(world, 1, world.size)
165 desc_Nr = grid_Nr.new_descriptor(N, R, N, r)
166 A_Nr = desc_Nr.zeros(dtype=dtype)
167 Redistributor(world, desc_mR,
168 desc_Nr).redistribute(A_mR, A_Nr)
169 A_Rr = np.delete(A_Nr, deleteN, axis=0)
170 A_Rr = np.ascontiguousarray(A_Rr)
171 desc_Rr = grid_Nr.new_descriptor(R, R, R, r)
173 # Redistribute to final grid.
174 # If this is not specified by the user, we try to find the most
175 # efficient grid using the suggest_blocking function.
176 if new_desc is None:
177 nrows, ncols, blocksize = suggest_blocking(R, world.size)
178 new_grid = BlacsGrid(world, nrows, ncols)
179 new_desc = new_grid.new_descriptor(R, R, blocksize, blocksize)
180 A_rr = new_desc.zeros(dtype=dtype)
181 Redistributor(world, desc_Rr,
182 new_desc).redistribute(A_Rr, A_rr)
184 return A_rr, new_desc
187@dataclass
188class ScreenedPotential:
189 pawcorr_q: list
190 W_qGG: list
191 qpd_q: list
194class SpinorData:
195 def __init__(self, con_m, val_m, e_km, f_km, v_kmn, soc_tol):
196 self.e_km = e_km
197 self.f_km = f_km
198 self.v0_kmn = v_kmn[:, :, ::2]
199 self.v1_kmn = v_kmn[:, :, 1::2]
201 mi = val_m[0]
202 mf = con_m[-1] + 1
203 tmp_n = np.argwhere(np.abs(v_kmn[:, mi:mf])**2 > soc_tol)[:, 2]
204 self.ni = np.min(tmp_n) // 2
205 self.nf = np.max(tmp_n) // 2 + 1
207 self.vslice_m = slice(val_m[0], val_m[-1] + 1)
208 self.cslice_m = slice(con_m[0], con_m[-1] + 1)
210 def _transform_rho(self, K1, K2, slice1, slice2,
211 rho0_nnG, rho1_nnG, susc_component='00'):
212 slice_n = slice(self.ni, self.nf)
213 vec0k1_mn = self.v0_kmn[K1, slice1, slice_n]
214 vec0k2_mn = self.v0_kmn[K2, slice2, slice_n]
215 vec1k1_mn = self.v1_kmn[K1, slice1, slice_n]
216 vec1k2_mn = self.v1_kmn[K2, slice2, slice_n]
217 if susc_component == '00':
218 if rho1_nnG is None:
219 rho1_nnG = rho0_nnG
220 rho0_mmG = np.dot(vec0k1_mn.conj(), np.dot(vec0k2_mn, rho0_nnG))
221 rho1_mmG = np.dot(vec1k1_mn.conj(), np.dot(vec1k2_mn, rho1_nnG))
222 return rho0_mmG + rho1_mmG
223 elif susc_component == '+-':
224 assert rho1_nnG is None
225 return np.dot(vec0k1_mn.conj(), np.dot(vec1k2_mn, rho0_nnG))
226 elif susc_component == '-+':
227 assert rho1_nnG is None
228 return np.dot(vec1k1_mn.conj(), np.dot(vec0k2_mn, rho0_nnG))
229 else:
230 raise NotImplementedError('Susceptibility component not '
231 'implemented. Please choose between '
232 '"00", "+-" or "-+"')
234 def rho_valence_valence(self, K1, K2, rho0_nnG, rho1_nnG=None):
235 return self._transform_rho(K1, K2, self.vslice_m,
236 self.vslice_m, rho0_nnG, rho1_nnG)
238 def rho_conduction_conduction(self, K1, K2, rho0_nnG, rho1_nnG=None):
239 return self._transform_rho(K1, K2, self.cslice_m, self.cslice_m,
240 rho0_nnG, rho1_nnG)
242 def rho_valence_conduction(self, K1, K2, rho0_nnG,
243 rho1_nnG=None, susc_component='00'):
244 return self._transform_rho(K1, K2, self.vslice_m, self.cslice_m,
245 rho0_nnG, rho1_nnG,
246 susc_component=susc_component)
248 def get_deps(self, K1, K2):
249 epsv_m = self.e_km[K1, self.vslice_m]
250 epsc_m = self.e_km[K2, self.cslice_m]
251 return -(epsv_m[:, np.newaxis] - epsc_m)
253 def get_df(self, K1, K2):
254 fv_m = self.f_km[K1, self.vslice_m]
255 fc_m = self.f_km[K2, self.cslice_m]
256 return fv_m[:, np.newaxis] - fc_m
259class BSEBackend:
260 def __init__(self, *, gs, context,
261 valence_bands,
262 conduction_bands,
263 deps_max=None,
264 add_soc=False,
265 soc_tol=0.0001,
266 ecut=10.,
267 scale=1.0,
268 nbands=None,
269 eshift=None,
270 gw_kn=None,
271 truncation=None,
272 integrate_gamma='reciprocal',
273 mode='BSE',
274 q_c=[0.0, 0.0, 0.0],
275 direction=0):
277 integrate_gamma = GammaIntegrationMode(integrate_gamma)
279 self.gs = gs
280 self.q_c = q_c
281 self.direction = direction
282 self.context = context
283 self.add_soc = add_soc
284 self.scale = scale
286 assert mode in ['RPA', 'BSE']
288 if deps_max is None:
289 self.deps_max = np.inf
290 else:
291 self.deps_max = deps_max / Hartree
292 self.ecut = ecut / Hartree
293 self.nbands = nbands
294 self.mode = mode
296 if integrate_gamma.is_analytical and truncation is not None:
297 self.context.print('***WARNING*** Analytical Coulomb integration' +
298 ' is not expected to work with Coulomb ' +
299 'truncation. ' +
300 'Use integrate_gamma=\'reciprocal\'')
301 self.integrate_gamma = integrate_gamma
303 # Find q-vectors and weights in the IBZ:
304 self.kd = self.gs.kd
305 if -1 in self.kd.bz2bz_ks:
306 self.context.print('***WARNING*** Symmetries may not be right. '
307 'Use gamma-centered grid to be sure')
308 offset_c = 0.5 * ((self.kd.N_c + 1) % 2) / self.kd.N_c
309 bzq_qc = monkhorst_pack(self.kd.N_c) + offset_c
310 self.qd = KPointDescriptor(bzq_qc)
311 self.qd.set_symmetry(self.gs.atoms, self.kd.symmetry)
313 # By default calculate the density-density response
314 self.susc_component = '00'
316 # Bands and spin
317 self.nspins = self.gs.nspins
318 self.val_m = self.parse_bands(valence_bands,
319 band_type='valence')
320 self.con_m = self.parse_bands(conduction_bands,
321 band_type='conduction')
323 self.use_tammdancoff = decide_whether_tammdancoff(self.val_m,
324 self.con_m)
326 self.nK = self.kd.nbzkpts
327 self.nv = len(self.val_m)
328 self.nc = len(self.con_m)
329 if eshift is not None:
330 eshift /= Hartree
331 if gw_kn is not None:
332 assert self.nv + self.nc == len(gw_kn[0])
333 assert self.kd.nibzkpts == len(gw_kn)
334 gw_kn = gw_kn[self.kd.bz2ibz_k]
335 gw_kn /= Hartree
336 self.gw_kn = gw_kn
337 self.eshift = eshift
339 self.coulomb = CoulombKernel.from_gs(self.gs, truncation=truncation)
341 # Distribution of kpoints
342 self.myKrange, self.myKsize = self.parallelisation_kpoints()
344 # Number of global and local pair orbitals. Note that self.ns is the
345 # the same everywhere and adds up to a value that is larger that nS.
346 # This is required for BlacsGrids in the ScalaPack diagonalization.
347 self.nS = self.nK * self.nv * self.nc
348 self.ns = -(-self.nK // world.size) * self.nv * self.nc
350 # Print all the details
351 self.print_initialization(self.use_tammdancoff, self.eshift,
352 self.gw_kn)
354 # Treat spin-polarized states as spinors without soc
355 if self.nspins == 2 and not self.add_soc:
356 self.add_soc = True
357 self.scale = 0.0
359 # Setup bands
360 if self.add_soc:
361 self.spinors_data = self._spinordata(soc_tol)
362 # Get a wide range of pair densities to allow for SOC mixing.
363 # The number of no-SOC states included are determined by soc_tol
364 # such that components of the soc eigenstates are included if
365 # their norm square is above soc_tol.
366 # The no-SOC pair densities are then transformed to the SOC pair
367 # densities. vi, vf, ci, cf here determines initial and final
368 # indices of no-SOC valence ond conduction states included.
369 self.vi = self.spinors_data.ni
370 self.vf = self.spinors_data.nf
371 self.ci = self.spinors_data.ni
372 self.cf = self.spinors_data.nf
373 else:
374 # Here we just need the pair densities of the specified bands
375 self.vi, self.vf = self.val_m[0], self.val_m[-1] + 1
376 self.ci, self.cf = self.con_m[0], self.con_m[-1] + 1
378 def parse_bands(self, bands, band_type='valence'):
379 """Helper function that checks whether bands are correctly specified,
380 and brings them to the format used later in the code.
382 Either integers (numbers of desired bands) or lists of band indices
383 must be provided. For spin-polarized calculations all
384 valence/condiction bands must be specified as a single list (one
385 for each) - regardless of spin. Same as if one includes SOC.
387 band_type is an optional parameter that is relevant when a desired
388 number of bands is given (rather than a list) to help figure out the
389 correct band indices.
390 """
391 if hasattr(bands, '__iter__'):
392 if not isinstance(bands[0], int):
393 raise ValueError('The bands must be specified as a single '
394 'list or an integer (number of bands).')
395 return bands
397 n_fully_occupied_bands, n_partially_occupied_bands = \
398 self.gs.count_occupied_bands()
400 if self.nspins == 2:
401 n_fully_occupied_bands += n_partially_occupied_bands
402 elif self.add_soc:
403 n_fully_occupied_bands *= 2
405 if band_type == 'valence':
406 bands_m = range(n_fully_occupied_bands - bands,
407 n_fully_occupied_bands)
408 elif band_type == 'conduction':
409 bands_m = range(n_fully_occupied_bands,
410 n_fully_occupied_bands + bands)
411 else:
412 raise ValueError(f'Invalid band type: {band_type}')
414 return bands_m
416 def _spinordata(self, soc_tol):
417 self.context.print('Diagonalizing spin-orbit Hamiltonian')
418 # We dont need ALL the screening bands to mix in the SOC here.
419 # This will give us bands up to two times the highest conduction band.
420 n2 = np.min([self.con_m[-1], self.nbands])
421 soc = self.gs.soc_eigenstates(n2=n2, scale=self.scale)
422 f_km = np.array([wf.f_m for wf in soc])
423 e_km = soc.eigenvalues()
424 e_km /= Hartree
425 v_kmn = soc.eigenvectors()
426 return SpinorData(self.con_m, self.val_m, e_km, f_km, v_kmn, soc_tol)
428 @timer('BSE calculate')
429 def calculate(self, optical, irreducible=False):
430 """Calculate the BSE Hamiltonian. This includes setting up all
431 machinery for pair densities, KS eignevalues and occupation factors.
432 At the end the direct and indirect interaction are included through
433 calls to separate functions.
435 The indices are indicative such that all capital letters imply
436 global indices and lower case letters imply local ("my CPU")
437 indices. For example, K, S and k, s are global and local k-point
438 and pair state indices respectively.
440 In addition the KS state indices are used such that n represents
441 states without SOC and m represents states with SOC. This is not always
442 possible though - the Hamiltonian, for example, is always
443 denoted H_kmmKmm - also for calculations without SOC. G is reciprocal
444 lattice index.
446 The parameter 'irreducible' puts V=0 such that the BSE kernel only
447 contians W. It is used for BSE+ calculations.
448 """
449 qpd0 = SingleQPWDescriptor.from_q(self.q_c, self.ecut, self.gs.gd)
451 self.ikq_k = self.kd.find_k_plus_q(self.q_c)
452 if irreducible:
453 self.v_G = np.zeros(qpd0.ng_q[0])
454 else:
455 self.v_G = self.coulomb.V(qpd=qpd0, q_v=None)
457 if optical:
458 self.v_G[0] = 0.0
460 context = ResponseContext(txt='pair.txt', timer=self.context.timer,
461 comm=serial_comm)
462 kptpair_factory = KPointPairFactory(gs=self.gs, context=context)
463 pair_calc = kptpair_factory.pair_calculator()
464 pawcorr = self.gs.pair_density_paw_corrections(qpd0)
466 if self.mode != 'RPA':
467 screened_potential = self.calculate_screened_potential()
468 else:
469 screened_potential = None
471 # Calculate pair densities, eigenvalues and occupations
472 self.context.timer.start('Pair densities')
473 if self.susc_component != '00':
474 assert self.nspins == 2
475 rhomag_KmmG = np.zeros((self.nK, self.nv,
476 self.nc, len(self.v_G)), complex)
477 rhoex_KmmG = np.zeros((self.nK, self.nv,
478 self.nc, len(self.v_G)), complex)
479 df_Kmm = np.zeros((self.nK, self.nv,
480 self.nc), float) # (fc - fv)
481 deps_kmm = np.zeros((self.myKsize, self.nv,
482 self.nc), float) # (ec - ev)
483 deps_Kmm = np.zeros((self.nK, self.nv, self.nc), float)
484 optical_limit = np.allclose(self.q_c, 0.0)
486 get_pair = kptpair_factory.get_kpoint_pair
487 get_pair_density = pair_calc.get_pair_density
489 # Calculate all properties diagonal in k-point
490 # These include the indirect (exchange) kernel,
491 # pseudo-energies, and occupation numbers
492 for ik, iK in enumerate(self.myKrange):
493 pair0 = get_pair(qpd0, 0, iK, self.vi, self.vf,
494 self.ci, self.cf)
495 v_n = np.arange(self.vi, self.vf)
496 c_n = np.arange(self.ci, self.cf)
497 iKq = self.gs.kd.find_k_plus_q(self.q_c, [iK])[0]
499 # Energies
500 if self.gw_kn is not None:
501 epsv_m = self.gw_kn[iK, :self.nv]
502 epsc_m = self.gw_kn[iKq, self.nv:]
503 deps_kmm[ik] = -(epsv_m[:, np.newaxis] - epsc_m)
504 elif self.add_soc:
505 deps_kmm[ik] = self.spinors_data.get_deps(iK, iKq)
506 else:
507 deps_kmm[ik] = -pair0.get_transition_energies()
508 if optical_limit:
509 deps_kmm[np.where(deps_kmm == 0)] = 1.0e-9
511 # Occupation factors
512 if self.add_soc:
513 df_Kmm[iK] = self.spinors_data.get_df(iK, iKq)
514 else:
515 df_Kmm[iK] = pair0.get_occupation_differences()
517 # Pair densities
518 rho0_nnG = get_pair_density(qpd0, pair0, v_n, c_n,
519 pawcorr=pawcorr)
520 if optical_limit:
521 n0_nnv = pair_calc.get_optical_pair_density_head(
522 qpd0, pair0, v_n, c_n)
523 rho0_nnG[:, :, 0] = n0_nnv[:, :, self.direction]
524 if self.nspins == 2:
525 pair1 = get_pair(qpd0, 1, iK, self.vi, self.vf,
526 self.ci, self.cf)
527 rho1_nnG = get_pair_density(qpd0, pair1, v_n, c_n,
528 pawcorr=pawcorr)
529 if optical_limit:
530 n1_nnv = pair_calc.get_optical_pair_density_head(
531 qpd0, pair1, v_n, c_n)
532 rho1_nnG[:, :, 0] = n1_nnv[:, :, self.direction]
533 deps1_nn = -pair1.get_transition_energies()
534 rho1_nnG[:, :, 0] *= deps1_nn
535 else:
536 rho1_nnG = None
538 # Generate the pair density matrix (with soc) used below
539 if self.add_soc:
540 if optical_limit:
541 deps0_nn = -pair0.get_transition_energies()
542 rho0_nnG[:, :, 0] *= deps0_nn
543 rhoex_KmmG[iK] = \
544 self.spinors_data.rho_valence_conduction(
545 iK, iKq, rho0_nnG, rho1_nnG)
546 if optical_limit:
547 rhoex_KmmG[iK, :, :, 0] /= deps_kmm[ik]
548 else:
549 rhoex_KmmG[iK] = rho0_nnG
551 # Generate the magnetic spin flip pair density for magnons
552 if self.susc_component != '00':
553 if self.susc_component == '+-':
554 s = 0
555 elif self.susc_component == '-+':
556 s = 1
557 pairflip = get_pair(qpd0, s, iK, self.vi, self.vf,
558 self.ci, self.cf, flipspin=True)
559 rhoflip_nnG = get_pair_density(qpd0, pairflip, v_n, c_n,
560 pawcorr=pawcorr)
561 rhomag_KmmG[iK] = self.spinors_data.rho_valence_conduction(
562 iK, iKq, rhoflip_nnG, susc_component=self.susc_component)
564 # Scissors operator shift
565 if self.eshift is not None:
566 deps_kmm[np.where(df_Kmm[self.myKrange] > 1e-3)] += self.eshift
567 deps_kmm[np.where(df_Kmm[self.myKrange] < -1e-3)] -= self.eshift
568 deps_Kmm[self.myKrange] = deps_kmm
570 world.sum(deps_Kmm)
571 world.sum(df_Kmm)
572 world.sum(rhoex_KmmG)
574 self.rhoG0_S = np.reshape(rhoex_KmmG[:, :, :, 0], -1)
575 self.rho_SG = np.reshape(rhoex_KmmG, (len(self.rhoG0_S), -1))
576 if self.susc_component != '00':
577 world.sum(rhomag_KmmG)
578 self.rhomag_SG = np.reshape(rhomag_KmmG, (self.nS, -1))
579 G_Gv = qpd0.get_reciprocal_vectors(add_q=False)
580 self.G_Gc = np.dot(G_Gv, qpd0.gd.cell_cv.T / (2 * np.pi))
581 self.context.timer.stop('Pair densities')
583 # Calculate Hamiltonian
584 self.context.timer.start('Calculate Hamiltonian')
585 t0 = time()
587 def update_progress(iK1):
588 dt = time() - t0
589 tleft = dt * self.myKsize / (iK1 + 1) - dt
591 self.context.print(
592 ' Finished %s pair orbitals in %s - Estimated %s left'
593 % ((iK1 + 1) * self.nv * self.nc * world.size,
594 timedelta(seconds=round(dt)),
595 timedelta(seconds=round(tleft))))
597 self.context.print('Calculating {} matrix elements at q_c = {}'.format(
598 self.mode, self.q_c))
600 # Hamiltonian buffer array
601 H_kmmKmm = np.zeros((self.myKsize, self.nv, self.nc,
602 self.nK, self.nv, self.nc),
603 complex)
605 # Add kernels to buffer array
606 self.add_indirect_kernel(kptpair_factory, rhoex_KmmG, H_kmmKmm)
607 if self.mode != 'RPA':
608 self.add_direct_kernel(kptpair_factory, pair_calc,
609 screened_potential, update_progress,
610 H_kmmKmm)
611 H_kmmKmm /= self.gs.volume
612 self.context.timer.stop('Calculate Hamiltonian')
614 if self.myKsize > 0:
615 iS0 = self.myKrange[0] * self.nv * self.nc
617 # multiply by 2 when spin-paired and no SOC
618 df_Kmm *= 2.0 / self.nK / (self.add_soc + 1)
619 df_S = np.reshape(df_Kmm, -1)
620 self.df_S = df_S
622 deps_S = np.reshape(deps_Kmm, -1)
623 deps_s = np.reshape(deps_kmm, -1)
625 mySsize = self.myKsize * self.nv * self.nc
626 H_sS = np.reshape(H_kmmKmm, (mySsize, self.nS))
627 for iS in range(mySsize):
628 # Multiply by occupations
629 H_sS[iS] *= df_S[iS0 + iS]
630 # add bare transition energies
631 H_sS[iS, iS0 + iS] += deps_s[iS]
633 return BSEMatrix(df_S, H_sS, deps_S, self.deps_max)
635 @timer('add_direct_kernel')
636 def add_direct_kernel(self, kptpair_factory, pair_calc, screened_potential,
637 update_progress, H_kmmKmm):
638 kpf = kptpair_factory
639 for ik1, iK1 in enumerate(self.myKrange):
640 kptv1_s = [kpf.get_k_point(s, iK1, self.vi, self.vf)
641 for s in range(self.nspins)]
642 kptc1_s = [kpf.get_k_point(s, self.ikq_k[iK1], self.ci, self.cf)
643 for s in range(self.nspins)]
644 for Q_c in self.qd.bzk_kc:
645 iK2 = self.kd.find_k_plus_q(Q_c, [kptv1_s[0].K])[0]
646 kptv2_s = [kptpair_factory.get_k_point(s, iK2, self.vi,
647 self.vf)
648 for s in range(self.nspins)]
649 kptc2_s = [kptpair_factory.get_k_point(s, self.ikq_k[iK2],
650 self.ci, self.cf)
651 for s in range(self.nspins)]
653 rho3_nnG, iq = self.get_density_matrix(
654 pair_calc, screened_potential, kptv1_s[0], kptv2_s[0])
656 rho4_nnG, iq = self.get_density_matrix(
657 pair_calc, screened_potential, kptc1_s[0], kptc2_s[0])
659 if self.nspins == 2:
660 rho3s1_nnG, iq = self.get_density_matrix(
661 pair_calc, screened_potential, kptv1_s[1], kptv2_s[1])
663 rho4s1_nnG, iq = self.get_density_matrix(
664 pair_calc, screened_potential, kptc1_s[1], kptc2_s[1])
665 else:
666 rho3s1_nnG = None
667 rho4s1_nnG = None
669 # Here we use n instead of m for the soc indices to save memory
670 if self.add_soc:
671 rho3_nnG = self.spinors_data.rho_valence_valence(
672 kptv1_s[0].K, kptv2_s[0].K, rho3_nnG, rho3s1_nnG)
674 rho4_nnG = self.spinors_data.rho_conduction_conduction(
675 kptc1_s[0].K, kptc2_s[0].K, rho4_nnG, rho4s1_nnG)
677 self.context.timer.start('Screened exchange')
678 W_mmmm = np.einsum(
679 'ijk,km,pqm->ipjq',
680 rho3_nnG.conj(),
681 screened_potential.W_qGG[iq],
682 rho4_nnG,
683 optimize='optimal')
684 # Only include 0.5*W for spinpaired calculations without soc
685 H_kmmKmm[ik1, :, :, iK2] -= W_mmmm * (self.add_soc + 1) / 2
686 self.context.timer.stop('Screened exchange')
688 if iK1 % (self.myKsize // 5 + 1) == 0:
689 update_progress(iK1=iK1)
691 @timer('add_indirect_kernel')
692 def add_indirect_kernel(self, kptpair_factory, rhoex_KmmG, H_kmmKmm):
693 for ik1, iK1 in enumerate(self.myKrange):
694 kptv1 = kptpair_factory.get_k_point(
695 0, iK1, self.vi, self.vf)
696 rho1V_mmG = rhoex_KmmG.conj()[iK1, :, :] * self.v_G
697 for Q_c in self.qd.bzk_kc:
698 iK2 = self.kd.find_k_plus_q(Q_c, [kptv1.K])[0]
699 rho2_mmG = rhoex_KmmG[iK2]
700 self.context.timer.start('Coulomb')
701 H_kmmKmm[ik1, :, :, iK2, :, :] += np.einsum(
702 'ijG,mnG->ijmn', rho1V_mmG, rho2_mmG,
703 optimize='optimal')
704 self.context.timer.stop('Coulomb')
706 @timer('get_density_matrix')
707 def get_density_matrix(self, pair_calc, screened_potential, kpt1, kpt2):
708 self.context.timer.start('Symop')
709 from gpaw.response.g0w0 import QSymmetryOp, get_nmG
710 symop, iq = QSymmetryOp.get_symop_from_kpair(self.kd, self.qd,
711 kpt1, kpt2)
712 qpd = screened_potential.qpd_q[iq]
713 nG = qpd.ngmax
714 pawcorr0 = screened_potential.pawcorr_q[iq]
715 pawcorr, I_G = symop.apply_symop_q(qpd, pawcorr0, kpt1, kpt2)
716 self.context.timer.stop('Symop')
718 rho_nnG = np.zeros((len(kpt1.eps_n), len(kpt2.eps_n), nG), complex)
719 for n in range(len(rho_nnG)):
720 rho_nnG[n] = get_nmG(kpt1, kpt2, pawcorr, n, qpd, I_G,
721 pair_calc, timer=self.context.timer)
723 return rho_nnG, iq
725 @cached_property
726 def _chi0calc(self):
727 return Chi0Calculator(
728 self.gs, self.context.with_txt('chi0.txt'),
729 wd=FrequencyDescriptor([0.0]),
730 eta=0.001,
731 ecut=self.ecut * Hartree,
732 intraband=False,
733 hilbert=False,
734 nbands=self.nbands)
736 @cached_property
737 def blockcomm(self):
738 return self._chi0calc.chi0_body_calc.blockcomm
740 @cached_property
741 def wcontext(self):
742 return ResponseContext(txt='w.txt', comm=world)
744 @cached_property
745 def _wcalc(self):
746 return initialize_w_calculator(
747 self._chi0calc, self.wcontext,
748 coulomb=self.coulomb,
749 integrate_gamma=self.integrate_gamma)
751 @timer('calculate_screened_potential')
752 def calculate_screened_potential(self):
753 """Calculate W_GG(q)."""
755 pawcorr_q = []
756 W_qGG = []
757 qpd_q = []
759 t0 = time()
760 self.context.print('Calculating screened potential')
761 for iq, q_c in enumerate(self.qd.ibzk_kc):
762 chi0 = self._chi0calc.calculate(q_c)
763 W_wGG = self._wcalc.calculate_W_wGG(chi0)
764 W_GG = W_wGG[0]
765 # This is such a terrible way to access the paw
766 # corrections. Attributes should not be groped like
767 # this... Change in the future! XXX
768 pawcorr_q.append(self._chi0calc.chi0_body_calc.pawcorr)
769 qpd_q.append(chi0.qpd)
770 W_qGG.append(W_GG)
772 if iq % (self.qd.nibzkpts // 5 + 1) == 2:
773 dt = time() - t0
774 tleft = dt * self.qd.nibzkpts / (iq + 1) - dt
775 self.context.print(
776 ' Finished {} q-points in {} - Estimated {} left'.format(
777 iq + 1, timedelta(seconds=round(dt)), timedelta(
778 seconds=round(tleft))))
780 return ScreenedPotential(pawcorr_q, W_qGG, qpd_q)
782 @timer('diagonalize')
783 def diagonalize_bse_matrix(self, bsematrix):
784 self.context.print('Diagonalizing Hamiltonian')
785 if self.use_tammdancoff:
786 return bsematrix.diagonalize_tammdancoff(self)
787 else:
788 return bsematrix.diagonalize_nontammdancoff(self)
790 @timer('get_bse_matrix')
791 def get_bse_matrix(self, optical=True, irreducible=False):
792 """Calculate BSE matrix."""
793 return self.calculate(optical=optical, irreducible=irreducible)
795 @timer('get_spectral_weights')
796 def get_spectral_weights(self, eig_data, df_S, mode_c):
797 if mode_c is None:
798 rho_S = self.rhoG0_S
799 else:
800 G_Gc = self.G_Gc
801 index = np.where(np.all(np.round(G_Gc) == mode_c, axis=1))[0][0]
802 rho_S = self.rhomag_SG[:, index]
804 w_T, v_St = eig_data[0], eig_data[1]
805 exclude_S = eig_data[2]
806 nS = self.nS - len(exclude_S)
807 ns = -(-nS // world.size)
808 dft_S = np.delete(df_S, exclude_S)
809 rhot_S = np.delete(rho_S, exclude_S)
810 C_T = np.zeros(nS, complex)
811 # Calculate the spectral weights C_T
812 if self.use_tammdancoff:
813 A_t = np.dot(rhot_S, v_St)
814 B_t = np.dot(rhot_S * dft_S, v_St)
815 if world.size == 1:
816 C_T = B_t.conj() * A_t
817 else:
818 grid = BlacsGrid(world, world.size, 1)
819 desc = grid.new_descriptor(nS, 1, ns, 1)
820 C_t = desc.empty(dtype=complex)
821 C_t[:, 0] = B_t.conj() * A_t
822 C_T = desc.collect_on_master(C_t)[:, 0]
823 if world.rank != 0:
824 C_T = np.empty(nS, dtype=complex)
825 world.broadcast(C_T, 0)
826 else:
827 if world.rank == 0:
828 A_T = np.dot(rhot_S, v_St)
829 B_T = np.dot(rhot_S * dft_S, v_St)
830 tmp = np.dot(v_St.conj().T, v_St)
831 overlap_TT = np.linalg.inv(tmp)
832 C_T = np.dot(B_T.conj(), overlap_TT.T) * A_T
833 world.broadcast(C_T, 0)
835 return w_T, C_T
837 def _cache_eig_data(self, irreducible, optical, w_w):
838 if (not hasattr(self, 'eig_data')
839 or self.eig_data_irreducible != irreducible
840 or self.eig_data_optical != optical):
841 bsematrix = self.get_bse_matrix(optical=optical,
842 irreducible=irreducible)
843 self.context.print('Calculating response function at %s frequency '
844 'points' % len(w_w))
845 self.eig_data = self.diagonalize_bse_matrix(bsematrix)
846 self.eig_data_irreducible = irreducible
847 self.eig_data_optical = optical
849 @timer('get_vchi')
850 def get_vchi(self, w_w=None, eta=0.1, optical=True, write_eig=None,
851 mode_c=None, irreducible=False):
852 """Returns v * chi where v is the bare Coulomb interaction"""
854 vchi_w = np.zeros(len(w_w), dtype=complex)
856 self._cache_eig_data(irreducible, optical, w_w)
858 w_T, C_T = self.get_spectral_weights(self.eig_data,
859 self.df_S, mode_c)
861 if write_eig is not None:
862 assert isinstance(write_eig, str)
863 filename = write_eig
864 if world.rank == 0:
865 write_bse_eigenvalues(filename, self.mode,
866 w_T * Hartree, C_T)
868 eta /= Hartree
869 for iw, w in enumerate(w_w / Hartree):
870 tmp_T = 1. / (w - w_T + 1j * eta)
871 vchi_w[iw] += np.dot(tmp_T, C_T)
872 vchi_w *= 4 * np.pi / self.gs.volume
874 if not np.allclose(self.q_c, 0.0):
875 cell_cv = self.gs.gd.cell_cv
876 B_cv = 2 * np.pi * np.linalg.inv(cell_cv).T
877 q_v = np.dot(self.q_c, B_cv)
878 vchi_w /= np.dot(q_v, q_v)
880 # Check f-sum rule
881 nv = self.gs.nvalence
882 dw_w = (w_w[1:] - w_w[:-1]) / Hartree
883 wvchi_w = (w_w[1:] * vchi_w[1:] + w_w[:-1] * vchi_w[:-1]) / Hartree / 2
884 N = -np.dot(dw_w, wvchi_w.imag) * self.gs.volume / (2 * np.pi**2)
885 Nt = 2 * np.dot(w_T, C_T).real
886 self.context.print('', flush=False)
887 self.context.print('Checking f-sum rule', flush=False)
888 self.context.print(f' Valence electrons : {nv}', flush=False)
889 self.context.print(f' Frequency integral: {N:f}', flush=False)
890 self.context.print(f' Sum of weights : {Nt:f}', flush=False)
891 self.context.print('')
893 return vchi_w
895 @timer('get_chi_wGG')
896 def get_chi_wGG(self, w_w=None, eta=0.1, readfile=None, optical=True,
897 irreducible=False):
898 """Returns chi_wGG'"""
900 self._cache_eig_data(irreducible, optical, w_w)
902 w_T, v_Rt, exclude_S = \
903 self.eig_data[0], self.eig_data[1], self.eig_data[2]
904 rho_SG = self.rho_SG
905 df_S = self.df_S
906 df_R = np.delete(df_S, exclude_S)
907 rho_RG = np.delete(rho_SG, exclude_S, axis=0)
909 nG = rho_RG.shape[-1]
910 nR = self.nS - len(exclude_S)
911 nr = -(-nR // world.size)
912 # nr is the local size of the array
914 self.context.print('Calculating response function at %s frequency '
915 'points' % len(w_w))
916 self.blocks = Blocks1D(world, len(w_T))
917 w_t = w_T[self.blocks.myslice]
919 if not self.use_tammdancoff:
920 if world.rank == 0:
921 v_RT = v_Rt
922 A_GT = rho_RG.T @ v_RT
923 B_GT = rho_RG.T * df_R[np.newaxis] @ v_RT
924 tmp = v_RT.conj().T @ v_RT
925 overlap_tt = np.linalg.inv(tmp)
926 C_tGG = ((B_GT.conj() @ overlap_tt.T).T)[..., np.newaxis] *\
927 A_GT.T[:, np.newaxis]
928 C_tGG = C_tGG[:nR].reshape((nR, nG, nG))
929 flat_C_tGG = C_tGG.ravel()
930 else:
931 flat_C_tGG = np.empty(nR * nG * nG, dtype=complex)
932 world.broadcast(flat_C_tGG, 0)
933 C_tGG = flat_C_tGG.reshape((nR, nG, nG))[self.blocks.myslice]
934 C_tGG1 = None
935 else:
936 A_Gt = rho_RG.T @ v_Rt
937 B_Gt = (rho_RG.T * df_R[np.newaxis]) @ v_Rt
938 '''The following computes
939 C_tGG1 = A_Gt.T.conj()[..., np.newaxis] * B_Gt.T[:, np.newaxis]
940 C_tGG = B_Gt.T.conj()[..., np.newaxis] * A_Gt.T[:, np.newaxis]
941 '''
942 grid = BlacsGrid(world, world.size, 1)
943 desc = grid.new_descriptor(nR, nG * nG, nr, nG * nG)
944 C_tGG = desc.empty(dtype=complex)
945 np.einsum('Gt,Ht->tGH', B_Gt.conj(), A_Gt,
946 out=C_tGG.reshape((-1, nG, nG)))
947 desc1 = grid.new_descriptor(nR, nG * nG, nr, nG * nG)
948 C_tGG1 = desc1.empty(dtype=complex)
949 np.einsum('Gt,Ht->tGH', A_Gt.conj(), B_Gt,
950 out=C_tGG1.reshape((-1, nG, nG)))
951 print(f'shape is {C_tGG.shape}')
952 C_tGG = C_tGG[:C_tGG.shape[0]].reshape((C_tGG.shape[0], nG, nG))
953 C_tGG1 = C_tGG1[:C_tGG1.shape[0]].reshape(
954 (C_tGG1.shape[0], nG, nG))
956 eta /= Hartree
958 if C_tGG is not None:
959 tmp_tw = 1 / (w_w[None, :] / Hartree - w_t[:, None] + 1j * eta)
960 chi_wGG_local = np.einsum('tw,tAB->wAB', tmp_tw, C_tGG)
962 if C_tGG1 is not None:
963 n_tmp_tw = - 1 / (w_w[None, :] / Hartree
964 + w_t[:, None] + 1j * eta)
965 chi_wGG_local += np.einsum('tw,tAB->wAB', n_tmp_tw, C_tGG1)
967 chi_wGG_local *= 1 / self.gs.volume
969 world.sum(chi_wGG_local)
970 chi_wGG = chi_wGG_local
972 return np.swapaxes(chi_wGG, -1, -2)
974 def get_dielectric_function(self, *args, filename='df_bse.csv', **kwargs):
975 vchi = self.vchi(*args, optical=True, **kwargs)
976 return vchi.dielectric_function(filename=filename)
978 def get_eels_spectrum(self, *args, filename='df_bse.csv', **kwargs):
979 vchi = self.vchi(*args, optical=False, **kwargs)
980 return vchi.eels_spectrum(filename=filename)
982 def get_polarizability(self, *args, filename='pol_bse.csv', **kwargs):
983 vchi = self.vchi(*args, optical=True, **kwargs)
984 return vchi.polarizability(filename=filename)
986 def get_magnetic_susceptibility(self, *args, modes_Gc=[[0, 0, 0]],
987 susc_component='+-',
988 write_eig='eig',
989 filename='susc_+-_bse_',
990 **kwargs):
991 """Returns and writes real and imaginary part of the magnetic
992 susceptibility.
994 susc_componenet: str
995 Component of the susceptibility tensor. '+-' and '-+'
996 are supported.
997 modes_Gc: list
998 List of reciprocal lattice vectors in reduced on which the
999 susceptibility is calculated. Default is the [0, 0, 0]
1000 component which gives the response over the Brillouin zone,
1001 but for optical magnons other components are needed.
1002 """
1003 self.susc_component = susc_component
1004 assert susc_component in ['+-', '-+']
1005 chi_Gw = []
1006 for mode_c in modes_Gc:
1007 assert len(mode_c) == 3
1008 assert all(isinstance(x, int) for x in mode_c)
1009 file_G = filename + ''.join(str(x) for x in mode_c) + '.csv'
1010 eig = write_eig + ''.join(str(x) for x in mode_c) + '.dat'
1011 vchi = self.vchi(*args, optical=False, mode_c=mode_c,
1012 write_eig=eig, **kwargs)
1013 chi_Gw.append(vchi.magnetic_susceptibility(filename=file_G))
1014 return chi_Gw
1016 def vchi(self, w_w=None, eta=0.1, write_eig='eig.dat',
1017 optical=None, mode_c=None):
1018 vchi_w = self.get_vchi(w_w=w_w, eta=eta, optical=optical,
1019 write_eig=write_eig, mode_c=mode_c)
1020 return VChi(self.gs.cd, self.context, w_w, vchi_w, optical=optical)
1022 def collect_A_SS(self, A_sS):
1023 if world.rank == 0:
1024 A_SS = np.zeros((self.nS, self.nS), dtype=complex)
1025 A_SS[:len(A_sS)] = A_sS
1026 Ntot = len(A_sS)
1027 for rank in range(1, world.size):
1028 buf = np.empty((self.ns, self.nS), dtype=complex)
1029 world.receive(buf, rank, tag=123)
1030 A_SS[Ntot:Ntot + self.ns] = buf
1031 Ntot += self.ns
1032 else:
1033 world.send(A_sS, 0, tag=123)
1034 world.barrier()
1035 if world.rank == 0:
1036 return A_SS
1038 def parallelisation_kpoints(self, rank=None):
1039 if rank is None:
1040 rank = world.rank
1041 nK = self.kd.nbzkpts
1042 myKsize = -(-nK // world.size)
1043 myKrange = range(rank * myKsize,
1044 min((rank + 1) * myKsize, nK))
1045 myKsize = len(myKrange)
1046 return myKrange, myKsize
1048 def print_initialization(self, td, eshift, gw_kn):
1049 isl = ['----------------------------------------------------------',
1050 f'{self.mode} Hamiltonian',
1051 '----------------------------------------------------------',
1052 f'Started at: {ctime()}', '',
1053 'Atoms : '
1054 f'{self.gs.atoms.get_chemical_formula(mode="hill")}',
1055 f'Ground state XC functional : {self.gs.xcname}',
1056 f'Valence electrons : {self.gs.nvalence}',
1057 f'Spinor calculations : {self.add_soc}',
1058 f'Number of bands : {self.gs.bd.nbands}',
1059 f'Number of spins : {self.gs.nspins}',
1060 f'Number of k-points : {self.kd.nbzkpts}',
1061 f'Number of irreducible k-points : {self.kd.nibzkpts}',
1062 f'Number of q-points : {self.qd.nbzkpts}',
1063 f'Number of irreducible q-points : {self.qd.nibzkpts}', '']
1065 for q in self.qd.ibzk_kc:
1066 isl.append(f' q: [{q[0]:1.4f} {q[1]:1.4f} {q[2]:1.4f}]')
1067 isl.append('')
1068 if gw_kn is not None:
1069 isl.append('User specified BSE bands')
1070 isl.extend([f'Response PW cutoff : {self.ecut * Hartree} '
1071 f'eV',
1072 f'Screening bands included : {self.nbands}'])
1073 isl.extend([f'Valence bands : {self.val_m}',
1074 f'Conduction bands : {self.con_m}'])
1075 if eshift is not None:
1076 isl.append(f'Scissors operator : {eshift * Hartree}'
1077 f'eV')
1078 isl.extend([
1079 f'Tamm-Dancoff approximation : {td}',
1080 f'Number of pair orbitals : {self.nS}',
1081 '',
1082 f'Truncation of Coulomb kernel : {self.coulomb.truncation}'])
1083 integrate_gamma = self.integrate_gamma.type
1084 if self.integrate_gamma.reduced:
1085 integrate_gamma += '2D'
1086 isl.append(
1087 f'Coulomb integration scheme : {integrate_gamma}')
1088 isl.extend([
1089 '',
1090 '----------------------------------------------------------',
1091 '----------------------------------------------------------',
1092 '',
1093 f'Parallelization - Total number of CPUs : {world.size}',
1094 ' Screened potential',
1095 f' K-point/band decomposition : {world.size}',
1096 ' Hamiltonian',
1097 f' Pair orbital decomposition : {world.size}'])
1098 self.context.print('\n'.join(isl))
1101class BSE(BSEBackend):
1102 def __init__(self, calc=None, timer=None, txt='-', **kwargs):
1103 """Creates the BSE object
1105 calc: str or calculator object
1106 The string should refer to the .gpw file contaning KS orbitals
1107 ecut: float
1108 Plane wave cutoff energy (eV)
1109 nbands: int
1110 Number of bands used for the screened interaction
1111 valence_bands: list or integer
1112 Valence bands used in the BSE Hamiltonian
1113 conduction_bands: list or integer
1114 Conduction bands used in the BSE Hamiltonian
1115 deps_max: float or None
1116 Maximum absolute value of transition energy for pair
1117 to be included in the BSE Hamiltonian
1118 add_soc: bool
1119 If True the calculation will included non-selfconsitent SOC.
1120 All band indices m refers to spinors, while n indices refer to
1121 states without SOC.
1122 scale: float
1123 Scaling of the SOC. A value of scale=1.0 yields a proper SOC
1124 calculation (id add_soc=True), whereas soc_scale=0 is equivalent
1125 to having add_soc=False.
1126 soc_tol: float
1127 Tolerance for how many non-SOC states are included when the SOC
1128 states are constructed (if add_soc=True). The SOC pair densities
1129 are constructed as linear combinations of pair densities without
1130 SOC. We include all states that contribute by more than soc_tol
1131 to the corresponding soc eigenstates.
1132 eshift: float
1133 Scissors operator opening the gap (eV)
1134 q_c: list of three floats
1135 Wavevector in reduced units on which the response is calculated
1136 direction: int
1137 if q_c = [0, 0, 0] this gives the direction in cartesian
1138 coordinates - 0=x, 1=y, 2=z
1139 gw_kn: list / array
1140 List or array defining the gw quasiparticle energies in eV
1141 used in the BSE Hamiltonian. Should match k-points and
1142 valence + conduction bands
1143 truncation: str or None
1144 Coulomb truncation scheme. Can be None or 2D.
1145 integrate_gamma: dict
1146 txt: str
1147 txt output
1148 mode: str
1149 Theory level used. can be RPA TDHF or BSE. Only BSE is screened.
1150 """
1151 gs, context = get_gs_and_context(
1152 calc, txt, world=world, timer=timer)
1154 super().__init__(gs=gs, context=context, **kwargs)
1157def write_bse_eigenvalues(filename, mode, w_w, C_w):
1158 with open(filename, 'w') as fd:
1159 print('# %s eigenvalues (in eV) and weights' % mode, file=fd)
1160 print('# Number eig weight', file=fd)
1161 for iw, (w, C) in enumerate(zip(w_w, C_w)):
1162 print('%8d %12.6f %12.16f' % (iw, w.real, C.real),
1163 file=fd)
1166def read_bse_eigenvalues(filename):
1167 _, w_w, C_w = np.loadtxt(filename, unpack=True)
1168 return w_w, C_w
1171def write_spectrum(filename, w_w, A_w):
1172 with open(filename, 'w') as fd:
1173 for w, A in zip(w_w, A_w):
1174 print(f'{w:.9f}, {A:.9f}', file=fd)
1177def read_spectrum(filename):
1178 w_w, A_w = np.loadtxt(filename, delimiter=',',
1179 unpack=True)
1180 return w_w, A_w
1183@dataclass
1184class VChi:
1185 cd: CellDescriptor
1186 context: ResponseContext
1187 w_w: np.ndarray
1188 vchi_w: np.ndarray
1189 optical: bool
1191 def epsilon(self):
1192 assert self.optical
1193 return -self.vchi_w + 1.0
1195 def eels(self):
1196 assert not self.optical
1197 return -self.vchi_w.imag
1199 def alpha(self):
1200 assert self.optical
1201 L = self.cd.nonperiodic_hypervolume
1202 return -L * self.vchi_w / (4 * np.pi)
1204 def susceptibility(self):
1205 assert not self.optical
1206 return self.vchi_w
1208 def dielectric_function(self, filename='df_bse.csv'):
1209 """Returns and writes real and imaginary part of the dielectric
1210 function.
1212 w_w: list of frequencies (eV)
1213 Dielectric function is calculated at these frequencies
1214 eta: float
1215 Lorentzian broadening of the spectrum (eV)
1216 filename: str
1217 data file on which frequencies, real and imaginary part of
1218 dielectric function is written
1219 write_eig: str
1220 File on which the BSE eigenvalues are written
1221 """
1223 return self._hackywrite(self.epsilon(), filename)
1225 # XXX The default filename clashes with that of dielectric function!
1226 def eels_spectrum(self, filename='df_bse.csv'):
1227 """Returns and writes real and imaginary part of the dielectric
1228 function.
1230 w_w: list of frequencies (eV)
1231 Dielectric function is calculated at these frequencies
1232 eta: float
1233 Lorentzian broadening of the spectrum (eV)
1234 filename: str
1235 data file on which frequencies, real and imaginary part of
1236 dielectric function is written
1237 write_eig: str
1238 File on which the BSE eigenvalues are written
1239 """
1240 return self._hackywrite(self.eels(), filename)
1242 def polarizability(self, filename='pol_bse.csv'):
1243 r"""Calculate the polarizability alpha.
1244 In 3D the imaginary part of the polarizability is related to the
1245 dielectric function by Im(eps_M) = 4 pi * Im(alpha). In systems
1246 with reduced dimensionality the converged value of alpha is
1247 independent of the cell volume. This is not the case for eps_M,
1248 which is ill defined. A truncated Coulomb kernel will always give
1249 eps_M = 1.0, whereas the polarizability maintains its structure.
1250 pbs should be a list of booleans giving the periodic directions.
1252 By default, generate a file 'pol_bse.csv'. The three colomns are:
1253 frequency (eV), Real(alpha), Imag(alpha). The dimension of alpha
1254 is \AA to the power of non-periodic directions.
1255 """
1256 return self._hackywrite(self.alpha(), filename)
1258 def magnetic_susceptibility(self, filename='susc_+-_0_bse.csv'):
1259 return self._hackywrite(self.susceptibility(), filename)[1]
1261 def _hackywrite(self, array, filename):
1262 if world.rank == 0 and filename is not None:
1263 if array.dtype == complex:
1264 write_response_function(filename, self.w_w, array.real,
1265 array.imag)
1266 else:
1267 assert array.dtype == float
1268 write_spectrum(filename, self.w_w, array)
1270 world.barrier()
1272 self.context.print('Calculation completed at:', ctime(), flush=False)
1273 self.context.print('')
1275 return self.w_w, array
1278class BSEPlus:
1280 def create_chi0_full_calculator(self):
1281 chi0calc_full = Chi0Calculator(self.gs, self.context,
1282 wd=self.wd,
1283 nbands=self.rpa_nbands,
1284 intraband=False,
1285 hilbert=False,
1286 eta=self.eta,
1287 ecut=self.ecut,
1288 eshift=self.eshift)
1290 return chi0calc_full
1292 def create_bse_calculator(self):
1293 bse = BSE(self.bse_gpw,
1294 ecut=self.ecut,
1295 valence_bands=self.bse_valence_bands,
1296 conduction_bands=self.bse_conduction_bands,
1297 nbands=self.bse_nbands,
1298 eshift=self.eshift,
1299 mode='BSE',
1300 truncation=self.truncation,
1301 q_c=self.q_c,
1302 direction=self.bse_direction,
1303 add_soc=self.bse_add_soc,
1304 txt='bse_calculation.txt')
1306 return bse
1308 def create_chi0_limited_calculator(self):
1309 self.gs, self.context = get_gs_and_context(
1310 self.rpa_gpw, txt=None, world=world, timer=None)
1311 self.wd = get_frequency_descriptor(
1312 self.w_w, gs=self.gs, nbands=self.rpa_nbands)
1313 chi0calc_limited = Chi0Calculator(self.gs, self.context,
1314 wd=self.wd,
1315 nbands=slice(
1316 self.n1_chi0, self.m2_chi0),
1317 intraband=False,
1318 hilbert=False,
1319 eta=self.eta,
1320 ecut=self.ecut,
1321 eshift=self.eshift)
1322 return chi0calc_limited
1324 def get_chi_RPA(self, chi0calc, q_c, coulomb_kernel, xc_kernel,
1325 CellDescriptor, direction):
1326 self.chi0_data = chi0calc.calculate(q_c)
1327 dyson_eqs = Chi0DysonEquations(self.chi0_data, coulomb_kernel,
1328 xc_kernel, CellDescriptor)
1329 self.v_G = coulomb_kernel.V(self.chi0_data.qpd)
1330 chi0_wGG = dyson_eqs.get_chi0_wGG(direction)
1331 chi0_WGG = dyson_eqs.wblocks.all_gather(chi0_wGG)
1332 del chi0calc, dyson_eqs, chi0_wGG
1333 return chi0_WGG
1335 def __init__(self,
1336 bse_gpw,
1337 bse_valence_bands,
1338 bse_conduction_bands,
1339 bse_nbands,
1340 rpa_gpw,
1341 rpa_nbands,
1342 w_w,
1343 eshift=0.0,
1344 bse_add_soc=False,
1345 eta=0.1,
1346 q_c=[0.0, 0.0, 0.0],
1347 direction=0,
1348 truncation=None,
1349 ecut=10):
1351 """ BSE+ calculation of chi. BSE+ offers a way to improve
1352 the convergence of the BSE by including transitions outside
1353 the active BSE electron-hole subspace at the RPA level in
1354 the irreducible polarizability. It saves the chi matrix
1355 calculated with the BSE+, BSE and RPA as npy-files.
1357 Parameters
1358 ----------
1359 bse_gpw: Path or str
1360 Name of the calculator that the BSE calculation should be
1361 made from (typically a fixed density calculator with less
1362 kpts)
1363 bse_valence_bands: range or list of integers
1364 Number of valence bands to be included in the bse calculation
1365 bse_conduction_bands: range or list of integers
1366 Number of conduction bands to be included in the bse calculation
1367 bse_nbands: integer
1368 Number of bands used for the screened interaction
1369 rpa_nbands: integer
1370 Number of bands to be included in the RPA calculation.
1371 w_w: list of floats
1372 Dielectric function is calculated at these frequencies (eV)
1373 eshift: float
1374 Scissors operator opening the gap (eV)
1375 bse_add_soc: bool
1376 If True the calculation will included non-self-consitent SOC in the
1377 underlying BSE calculation. SOC is not implemented in the RPA code.
1378 eta: float
1379 Lorentzian broadening of the spectrum (eV)
1380 q_c: list of three floats
1381 Wavevector in reduced units on which the response is calculated
1382 direction: int
1383 If q_c = [0, 0, 0] this gives the direction in cartesian
1384 coordinates - 0=x, 1=y, 2=z
1385 truncation: str or None
1386 Coulomb truncation scheme. Can be None or 2D.
1387 ecut: float
1388 Plane wave cutoff energy (eV)
1389 """
1391 self.bse_gpw = bse_gpw
1392 self.bse_valence_bands = bse_valence_bands
1393 self.bse_conduction_bands = bse_conduction_bands
1394 self.bse_nbands = bse_nbands
1395 self.rpa_gpw = rpa_gpw
1396 self.rpa_nbands = rpa_nbands
1397 self.w_w = w_w
1398 self.eshift = eshift
1399 self.bse_add_soc = bse_add_soc
1400 self.eta = eta
1401 self.q_c = q_c
1402 self.bse_direction = direction
1403 self.rpa_direction = ('x', 'y', 'z')[direction]
1404 self.truncation = truncation
1405 self.ecut = ecut
1407 self.n1_BSE = self.bse_valence_bands[0]
1408 self.m2_BSE = self.bse_conduction_bands[-1]
1409 self.m2_chi0_full = rpa_nbands - 1
1411 if bse_add_soc:
1412 self.n1_chi0 = int(self.n1_BSE / 2)
1413 self.m2_chi0 = int((self.m2_BSE + 1) / 2)
1414 else:
1415 self.n1_chi0 = self.n1_BSE
1416 self.m2_chi0 = self.m2_BSE + 1
1418 assert truncation in [None, '2D']
1420 assert self.m2_chi0 < self.m2_chi0_full, \
1421 'Large chi0 calculation should contain more ' \
1422 'bands than the BSE calculation'
1424 def calculate_chi_wGG(self, optical=True, xc_kernel=None, comm=None,
1425 bsep_name='chi_BSEPlus',
1426 save_chi_BSE=False, save_chi_RPA=False):
1428 # irreducibale bse chi
1429 bse = self.create_bse_calculator()
1430 chi_irr_BSE_WGG = bse.get_chi_wGG(
1431 eta=self.eta,
1432 optical=optical,
1433 irreducible=True,
1434 w_w=self.w_w)
1435 del bse
1437 # chi0 calculation with the same bands as in the bse
1438 chi0calc_limited = self.create_chi0_limited_calculator()
1439 coulomb_kernel = CoulombKernel.from_gs(self.gs,
1440 truncation=self.truncation)
1441 chi0_limited_WGG = self.get_chi_RPA(chi0calc_limited, self.q_c,
1442 coulomb_kernel, xc_kernel,
1443 self.gs.cd, self.rpa_direction)
1445 # chi0 fully converged
1446 chi0calc_full = self.create_chi0_full_calculator()
1447 chi0_full_WGG = self.get_chi_RPA(chi0calc_full, self.q_c,
1448 coulomb_kernel, xc_kernel,
1449 self.gs.cd, self.rpa_direction)
1451 if self.truncation == '2D':
1452 pbc_c = self.gs.pbc
1453 assert sum(pbc_c) == 2
1454 coulomb_kernel_bare = CoulombKernel.from_gs(
1455 self.gs, truncation=None)
1456 v_G_bare = coulomb_kernel_bare.V(self.chi0_data.qpd, q_v=None)
1457 self.v_G = self.v_G / v_G_bare
1458 if optical:
1459 v_G_bare[0] = 0.0
1460 chi0_limited_WGG = chi0_limited_WGG * \
1461 v_G_bare[np.newaxis, np.newaxis, :]
1462 chi0_full_WGG = chi0_full_WGG * v_G_bare[np.newaxis, np.newaxis, :]
1463 chi_irr_BSE_WGG = chi_irr_BSE_WGG * \
1464 v_G_bare[np.newaxis, np.newaxis, :]
1465 cell_cv = self.gs.gd.cell_cv
1466 V = np.abs(np.linalg.det(cell_cv[~pbc_c][:, ~pbc_c]))
1467 V *= Bohr
1468 elif self.truncation is None and optical:
1469 self.v_G[0] = 0.0
1471 del self.chi0_data
1473 if comm is None:
1474 comm = world
1475 nR = len(self.w_w)
1476 self.blocks = Blocks1D(comm, nR)
1478 chi_irr_BSE_wGG = chi_irr_BSE_WGG[self.blocks.myslice]
1479 chi0_full_wGG = chi0_full_WGG[self.blocks.myslice]
1480 chi0_limited_wGG = chi0_limited_WGG[self.blocks.myslice]
1482 chi_irr_BSEPlus_wGG = \
1483 chi_irr_BSE_wGG - chi0_limited_wGG + chi0_full_wGG
1484 eye = np.eye(chi_irr_BSEPlus_wGG.shape[1])
1486 chi_BSEPlus_wGG = \
1487 np.linalg.solve(eye - chi_irr_BSEPlus_wGG @ np.diag(self.v_G),
1488 chi_irr_BSEPlus_wGG)
1490 if self.truncation == '2D':
1491 chi_BSEPlus_wGG *= V / (4 * np.pi)
1493 chi_BSEPlus_WGG = self.blocks.gather(chi_BSEPlus_wGG, 0)
1495 if world.rank == 0:
1496 np.save(bsep_name + '.npy', chi_BSEPlus_WGG)
1497 del chi_BSEPlus_WGG
1498 del chi_BSEPlus_wGG, chi0_limited_wGG
1500 if save_chi_BSE:
1501 chi_BSE_wGG = \
1502 np.linalg.solve(eye - chi_irr_BSE_wGG @ np.diag(self.v_G),
1503 chi_irr_BSE_wGG)
1505 if self.truncation == '2D':
1506 chi_BSE_wGG *= V / (4 * np.pi)
1508 chi_BSE_WGG = self.blocks.gather(chi_BSE_wGG, 0)
1510 if world.rank == 0:
1511 np.save('chi_BSE.npy' if save_chi_BSE is True else
1512 save_chi_BSE, chi_BSE_WGG)
1513 del chi_BSE_WGG
1514 del chi_BSE_wGG
1516 if save_chi_RPA:
1517 chi_full_wGG = \
1518 np.linalg.solve(eye - chi0_full_wGG @ np.diag(self.v_G),
1519 chi0_full_wGG)
1521 if self.truncation == '2D':
1522 chi_full_wGG *= V / (4 * np.pi)
1524 chi_full_WGG = self.blocks.gather(chi_full_wGG, 0)
1526 if world.rank == 0:
1527 np.save('chi_RPA.npy' if save_chi_RPA is True else
1528 save_chi_RPA, chi_full_WGG)
1529 del chi_full_WGG
1530 del chi_full_wGG