1import warnings 2from math import sqrt, pi 3import numpy as np 4 5from ase.units import Ha 6from gpaw import BadParallelization 7from gpaw.mpi import world 8from gpaw.density import redistribute_array, redistribute_atomic_matrices 9from gpaw.sphere.lebedev import weight_n 10from gpaw.utilities import pack, pack_atomic_matrices, unpack_atomic_matrices 11from gpaw.xc.gllb import safe_sqr 12from gpaw.xc.gllb.contribution import Contribution 13 14# XXX Work in process 15debug = False 16 17 18def d(*args): 19 if debug: 20 print(args) 21 22 23class ResponsePotential: 24 """Container for response potential""" 25 def __init__(self, response, vt_sG, vt_sg, D_asp, Dresp_asp): 26 self.response = response 27 self.vt_sG = vt_sG 28 self.vt_sg = vt_sg 29 self.D_asp = D_asp 30 self.Dresp_asp = Dresp_asp 31 32 def redistribute(self, new_response): 33 old_response = self.response 34 new_wfs = new_response.wfs 35 new_density = new_response.density 36 self.vt_sG = redistribute_array(self.vt_sG, 37 old_response.density.gd, 38 new_density.gd, 39 new_wfs.nspins, 40 new_wfs.kptband_comm) 41 if self.vt_sg is not None: 42 self.vt_sg = redistribute_array(self.vt_sg, 43 old_response.density.finegd, 44 new_density.finegd, 45 new_wfs.nspins, 46 new_wfs.kptband_comm) 47 48 def redist_asp(D_asp): 49 return redistribute_atomic_matrices(D_asp, 50 new_density.gd, 51 new_wfs.nspins, 52 new_density.setups, 53 new_density.atom_partition, 54 new_wfs.kptband_comm) 55 56 self.D_asp = redist_asp(self.D_asp) 57 self.Dresp_asp = redist_asp(self.Dresp_asp) 58 self.response = new_response 59 60 61class C_Response(Contribution): 62 def __init__(self, weight, coefficients, damp=1e-10): 63 Contribution.__init__(self, weight) 64 d('In c_Response __init__', self) 65 self.coefficients = coefficients 66 self.vt_sg = None 67 self.vt_sG = None 68 self.D_asp = None 69 self.Dresp_asp = None 70 self.Ddist_asp = None 71 self.Drespdist_asp = None 72 self.damp = damp 73 self.fix_potential = False 74 75 def set_damp(self, damp): 76 self.damp = damp 77 78 def get_name(self): 79 return 'RESPONSE' 80 81 def get_desc(self): 82 return self.coefficients.get_description() 83 84 def initialize(self, density, hamiltonian, wfs): 85 Contribution.initialize(self, density, hamiltonian, wfs) 86 self.coefficients.initialize(wfs) 87 if self.Dresp_asp is None: 88 assert self.density.D_asp is None 89 # XXX With the deprecated `fixdensity=True` option this 90 # initialize() function is called both before AND after read()! 91 # Thus, the second call of initialize() would discard the read 92 # potential unless we check if it was already allocated. 93 # Remove this if statement once `fixdensity=True` option has 94 # been removed from the calculator. 95 if self.vt_sg is None: 96 self.vt_sG = self.gd.empty(self.nspins) 97 self.vt_sg = self.finegd.empty(self.nspins) 98 99 def initialize_1d(self, ae): 100 Contribution.initialize_1d(self, ae) 101 self.coefficients.initialize_1d(ae) 102 103 def initialize_from_other_response(self, response): 104 pot = ResponsePotential(response, 105 response.vt_sG, 106 response.vt_sg, 107 response.D_asp, 108 response.Dresp_asp) 109 pot.redistribute(self) 110 self.vt_sG = pot.vt_sG 111 self.vt_sg = pot.vt_sg 112 self.D_asp = pot.D_asp 113 self.Dresp_asp = pot.Dresp_asp 114 self.Ddist_asp = self.distribute_D_asp(pot.D_asp) 115 self.Drespdist_asp = self.distribute_D_asp(pot.Dresp_asp) 116 117 # Calcualte the GLLB potential and energy 1d 118 def add_xc_potential_and_energy_1d(self, v_g): 119 w_i = self.coefficients.get_coefficients_1d() 120 u2_j = safe_sqr(self.ae.u_j) 121 v_g += self.weight * np.dot(w_i, u2_j) / ( 122 np.dot(self.ae.f_j, u2_j) + self.damp) 123 return 0.0 # Response part does not contribute to energy 124 125 def update_potentials(self): 126 d('In update response potential') 127 if self.fix_potential: 128 # Skip the evaluation of the potential so that 129 # the existing potential is used. This is relevant for 130 # band structure calculation so that the potential 131 # does not get updated with the other k-point sampling. 132 return 133 134 if self.wfs.kpt_u[0].eps_n is None: 135 # This skips the evaluation of the potential so that 136 # the existing one is used. 137 # This happens typically after reading when occupations 138 # haven't been calculated yet and the potential read earlier 139 # is used. 140 return 141 142 w_kn = self.coefficients.get_coefficients(self.wfs.kpt_u) 143 f_kn = [kpt.f_n for kpt in self.wfs.kpt_u] 144 if w_kn is not None: 145 self.vt_sG[:] = 0.0 146 nt_sG = self.gd.zeros(self.nspins) 147 148 for kpt, w_n in zip(self.wfs.kpt_u, w_kn): 149 self.wfs.add_to_density_from_k_point_with_occupation( 150 self.vt_sG, kpt, w_n) 151 self.wfs.add_to_density_from_k_point(nt_sG, kpt) 152 153 self.wfs.kptband_comm.sum(nt_sG) 154 self.wfs.kptband_comm.sum(self.vt_sG) 155 156 if self.wfs.kd.symmetry: 157 for nt_G, vt_G in zip(nt_sG, self.vt_sG): 158 self.wfs.kd.symmetry.symmetrize(nt_G, self.gd) 159 self.wfs.kd.symmetry.symmetrize(vt_G, self.gd) 160 161 d('response update D_asp', world.rank, self.Dresp_asp.keys(), 162 self.D_asp.keys()) 163 self.wfs.calculate_atomic_density_matrices_with_occupation( 164 self.Dresp_asp, w_kn) 165 self.Drespdist_asp = self.distribute_D_asp(self.Dresp_asp) 166 d('response update Drespdist_asp', world.rank, 167 self.Dresp_asp.keys(), self.D_asp.keys()) 168 self.wfs.calculate_atomic_density_matrices_with_occupation( 169 self.D_asp, f_kn) 170 self.Ddist_asp = self.distribute_D_asp(self.D_asp) 171 172 self.vt_sG /= nt_sG + self.damp 173 174 self.density.distribute_and_interpolate(self.vt_sG, self.vt_sg) 175 176 def calculate(self, e_g, n_sg, v_sg): 177 self.update_potentials() 178 v_sg += self.weight * self.vt_sg 179 180 def distribute_D_asp(self, D_asp): 181 return self.density.atomdist.to_work(D_asp) 182 183 def calculate_energy_and_derivatives(self, setup, D_sp, H_sp, a, 184 addcoredensity=True): 185 # Get the XC-correction instance 186 c = setup.xc_correction 187 ncresp_g = setup.extra_xc_data['core_response'] / self.nspins 188 if not addcoredensity: 189 ncresp_g[:] = 0.0 190 191 for D_p, dEdD_p, Dresp_p in zip(self.Ddist_asp[a], H_sp, 192 self.Drespdist_asp[a]): 193 D_Lq = np.dot(c.B_pqL.T, D_p) 194 n_Lg = np.dot(D_Lq, c.n_qg) # Construct density 195 if addcoredensity: 196 n_Lg[0] += c.nc_g * sqrt(4 * pi) / self.nspins 197 nt_Lg = np.dot( 198 D_Lq, c.nt_qg 199 ) # Construct smooth density (_without smooth core density_) 200 201 Dresp_Lq = np.dot(c.B_pqL.T, Dresp_p) 202 nresp_Lg = np.dot(Dresp_Lq, c.n_qg) # Construct 'response density' 203 nrespt_Lg = np.dot( 204 Dresp_Lq, c.nt_qg 205 ) # Construct smooth 'response density' (w/o smooth core) 206 207 for w, Y_L in zip(weight_n, c.Y_nL): 208 nt_g = np.dot(Y_L, nt_Lg) 209 nrespt_g = np.dot(Y_L, nrespt_Lg) 210 x_g = nrespt_g / (nt_g + self.damp) 211 dEdD_p -= self.weight * w * np.dot( 212 np.dot(c.B_pqL, Y_L), np.dot(c.nt_qg, x_g * c.rgd.dv_g)) 213 214 n_g = np.dot(Y_L, n_Lg) 215 nresp_g = np.dot(Y_L, nresp_Lg) 216 x_g = (nresp_g + ncresp_g) / (n_g + self.damp) 217 218 dEdD_p += self.weight * w * np.dot( 219 np.dot(c.B_pqL, Y_L), np.dot(c.n_qg, x_g * c.rgd.dv_g)) 220 221 return 0.0 222 223 def integrate_sphere(self, a, Dresp_sp, D_sp, Dwf_p, spin=0): 224 c = self.wfs.setups[a].xc_correction 225 Dresp_p, D_p = Dresp_sp[spin], D_sp[spin] 226 D_Lq = np.dot(c.B_pqL.T, D_p) 227 n_Lg = np.dot(D_Lq, c.n_qg) # Construct density 228 n_Lg[0] += c.nc_g * sqrt(4 * pi) / len(D_sp) 229 nt_Lg = np.dot(D_Lq, c.nt_qg 230 ) # Construct smooth density (without smooth core) 231 Dresp_Lq = np.dot(c.B_pqL.T, Dresp_p) # Construct response 232 nresp_Lg = np.dot(Dresp_Lq, c.n_qg) # Construct 'response density' 233 nrespt_Lg = np.dot( 234 Dresp_Lq, c.nt_qg 235 ) # Construct smooth 'response density' (w/o smooth core) 236 Dwf_Lq = np.dot(c.B_pqL.T, Dwf_p) # Construct lumo wf 237 nwf_Lg = np.dot(Dwf_Lq, c.n_qg) 238 nwft_Lg = np.dot(Dwf_Lq, c.nt_qg) 239 E = 0.0 240 for w, Y_L in zip(weight_n, c.Y_nL): 241 v = np.dot(Y_L, nwft_Lg) * np.dot(Y_L, nrespt_Lg) / ( 242 np.dot(Y_L, nt_Lg) + self.damp) 243 E -= self.weight * w * np.dot(v, c.rgd.dv_g) 244 v = np.dot(Y_L, nwf_Lg) * np.dot(Y_L, nresp_Lg) / ( 245 np.dot(Y_L, n_Lg) + self.damp) 246 E += self.weight * w * np.dot(v, c.rgd.dv_g) 247 return E 248 249 def add_smooth_xc_potential_and_energy_1d(self, vt_g): 250 w_ln = self.coefficients.get_coefficients_1d(smooth=True) 251 v_g = np.zeros(self.ae.N) 252 n_g = np.zeros(self.ae.N) 253 for w_n, f_n, u_n in zip(w_ln, self.ae.f_ln, 254 self.ae.s_ln): # For each angular momentum 255 u2_n = safe_sqr(u_n) 256 v_g += np.dot(w_n, u2_n) 257 n_g += np.dot(f_n, u2_n) 258 259 vt_g += self.weight * v_g / (n_g + self.damp) 260 return 0.0 # Response part does not contribute to energy 261 262 def calculate_delta_xc(self, homolumo=None): 263 warnings.warn( 264 'The function calculate_delta_xc() is deprecated. ' 265 'Please use calculate_discontinuity_potential() instead. ' 266 'See documentation on calculating band gap with GLLBSC.', 267 DeprecationWarning) 268 269 if homolumo is None: 270 # Calculate band gap 271 272 # This always happens, so we don't warn. 273 # We should perhaps print it as an ordinary message, 274 # but we do not have a file here to which to print. 275 # import warnings 276 # warnings.warn('Calculating KS-gap directly from the k-points, ' 277 # 'can be inaccurate.') 278 homolumo = self.wfs.get_homo_lumo() 279 homo, lumo = homolumo 280 dxc_pot = self.calculate_discontinuity_potential(homo * Ha, lumo * Ha) 281 self.Dxc_vt_sG = dxc_pot.vt_sG 282 self.Dxc_D_asp = dxc_pot.D_asp 283 self.Dxc_Dresp_asp = dxc_pot.Dresp_asp 284 285 def calculate_discontinuity_potential(self, homo, lumo): 286 r"""Calculate the derivative discontinuity potential. 287 288 This function calculates $`\Delta_{xc}(r)`$ as given by 289 Eq. (24) of https://doi.org/10.1103/PhysRevB.82.115106 290 291 Parameters: 292 293 homo: 294 homo energy (or energies in spin-polarized case) in eV 295 lumo: 296 lumo energy (or energies in spin-polarized case) in eV 297 298 Returns: 299 300 dxc_pot: Discontinuity potential 301 """ 302 homolumo = np.asarray((homo, lumo)) / Ha 303 304 dxc_Dresp_asp = self.empty_atomic_matrix() 305 dxc_D_asp = self.empty_atomic_matrix() 306 for a in self.density.D_asp: 307 ni = self.wfs.setups[a].ni 308 dxc_Dresp_asp[a] = np.zeros((self.nspins, ni * (ni + 1) // 2)) 309 dxc_D_asp[a] = np.zeros((self.nspins, ni * (ni + 1) // 2)) 310 311 # Calculate new response potential with LUMO reference 312 w_kn = self.coefficients.get_coefficients_for_lumo_perturbation( 313 self.wfs.kpt_u, 314 homolumo=homolumo) 315 f_kn = [kpt.f_n for kpt in self.wfs.kpt_u] 316 317 dxc_vt_sG = self.gd.zeros(self.nspins) 318 nt_sG = self.gd.zeros(self.nspins) 319 for kpt, w_n in zip(self.wfs.kpt_u, w_kn): 320 self.wfs.add_to_density_from_k_point_with_occupation(dxc_vt_sG, 321 kpt, w_n) 322 self.wfs.add_to_density_from_k_point(nt_sG, kpt) 323 324 self.wfs.kptband_comm.sum(nt_sG) 325 self.wfs.kptband_comm.sum(dxc_vt_sG) 326 327 if self.wfs.kd.symmetry: 328 for nt_G, dxc_vt_G in zip(nt_sG, dxc_vt_sG): 329 self.wfs.kd.symmetry.symmetrize(nt_G, self.gd) 330 self.wfs.kd.symmetry.symmetrize(dxc_vt_G, self.gd) 331 332 dxc_vt_sG /= nt_sG + self.damp 333 334 self.wfs.calculate_atomic_density_matrices_with_occupation( 335 dxc_Dresp_asp, w_kn) 336 self.wfs.calculate_atomic_density_matrices_with_occupation( 337 dxc_D_asp, f_kn) 338 dxc_pot = ResponsePotential(self, dxc_vt_sG, None, 339 dxc_D_asp, dxc_Dresp_asp) 340 return dxc_pot 341 342 def calculate_delta_xc_perturbation(self): 343 warnings.warn( 344 'The function calculate_delta_xc_perturbation() is deprecated. ' 345 'Please use calculate_discontinuity() instead. ' 346 'See documentation on calculating band gap with GLLBSC.', 347 DeprecationWarning) 348 dxc_pot = ResponsePotential(self, self.Dxc_vt_sG, None, self.Dxc_D_asp, 349 self.Dxc_Dresp_asp) 350 if self.nspins != 1: 351 ret = [] 352 for spin in range(self.nspins): 353 ret.append(self.calculate_discontinuity(dxc_pot, spin=spin)) 354 return ret 355 return self.calculate_discontinuity(dxc_pot) 356 357 def calculate_delta_xc_perturbation_spin(self, s=0): 358 warnings.warn( 359 'The function calculate_delta_xc_perturbation_spin() is ' 360 'deprecated. Please use calculate_discontinuity_spin() instead. ' 361 'See documentation on calculating band gap with GLLBSC.', 362 DeprecationWarning) 363 dxc_pot = ResponsePotential(self, self.Dxc_vt_sG, None, self.Dxc_D_asp, 364 self.Dxc_Dresp_asp) 365 return self.calculate_discontinuity(dxc_pot, spin=s) 366 367 def calculate_discontinuity(self, 368 dxc_pot: ResponsePotential, 369 spin: int = None): 370 r"""Calculate the derivative discontinuity. 371 372 This function evaluates the perturbation theory expression 373 for the derivative discontinuity value as given by 374 Eq. (25) of https://doi.org/10.1103/PhysRevB.82.115106 375 376 Parameters: 377 378 dxc_pot: 379 Discontinuity potential. 380 spin: 381 Spin value. 382 383 Returns: 384 385 KSgap: 386 Kohn-Sham gap in eV 387 dxc: 388 Derivative discontinuity in eV 389 """ 390 if spin is None: 391 if self.nspins != 1: 392 raise ValueError('Specify spin for the discontinuity.') 393 spin = 0 394 395 # Redistribute discontinuity potential 396 if dxc_pot.response is not self: 397 dxc_pot.redistribute(self) 398 399 homo, lumo = self.wfs.get_homo_lumo(spin) 400 KSgap = lumo - homo 401 402 # Find the lumo-orbital of this spin 403 if self.wfs.bd.comm.size != 1: 404 raise BadParallelization( 405 'Band parallelization is not supported by ' 406 'calculate_discontinuity().') 407 eps_n = self.wfs.kpt_qs[0][spin].eps_n 408 lumo_n = (eps_n < self.wfs.fermi_level).sum() 409 410 # Calculate the expectation value for all k points, and keep 411 # the smallest energy value 412 nt_G = self.gd.empty() 413 min_energy = np.inf 414 for u, kpt in enumerate(self.wfs.kpt_u): 415 if kpt.s == spin: 416 nt_G[:] = 0.0 417 self.wfs.add_orbital_density(nt_G, kpt, lumo_n) 418 E = 0.0 419 for a in dxc_pot.D_asp: 420 D_sp = dxc_pot.D_asp[a] 421 Dresp_sp = dxc_pot.Dresp_asp[a] 422 P_ni = kpt.P_ani[a] 423 Dwf_p = pack(np.outer(P_ni[lumo_n].T.conj(), 424 P_ni[lumo_n]).real) 425 E += self.integrate_sphere(a, Dresp_sp, D_sp, Dwf_p, 426 spin=spin) 427 E = self.gd.comm.sum(E) 428 E += self.gd.integrate(nt_G * dxc_pot.vt_sG[spin]) 429 E += kpt.eps_n[lumo_n] - lumo 430 min_energy = min(min_energy, E) 431 432 # Take the smallest value over all distributed k points 433 dxc = self.wfs.kd.comm.min(min_energy) 434 return KSgap * Ha, dxc * Ha 435 436 def calculate_discontinuity_from_average(self, 437 dxc_pot: ResponsePotential, 438 spin: int, 439 testing: bool = False): 440 msg = ('This function estimates discontinuity by calculating ' 441 'the average of the discontinuity potential. ' 442 'Use only for testing and debugging.') 443 if not testing: 444 raise RuntimeError(msg) 445 else: 446 warnings.warn(msg) 447 assert self.wfs.world.size == 1 448 449 # Calculate average of lumo reference response potential 450 dxc = np.average(dxc_pot.vt_sG[spin]) 451 return dxc * Ha 452 453 def initialize_from_atomic_orbitals(self, basis_functions): 454 # Initialize 'response-density' and density-matrices 455 # print('Initializing from atomic orbitals') 456 self.Dresp_asp = self.empty_atomic_matrix() 457 setups = self.wfs.setups 458 459 for a in self.density.D_asp.keys(): 460 ni = setups[a].ni 461 self.Dresp_asp[a] = np.zeros((self.nspins, ni * (ni + 1) // 2)) 462 463 self.D_asp = self.empty_atomic_matrix() 464 f_asi = {} 465 w_asi = {} 466 467 for a in basis_functions.atom_indices: 468 if setups[a].type == 'ghost': 469 w_j = [0.0] 470 else: 471 w_j = setups[a].extra_xc_data['w_j'] 472 473 # Basis function coefficients based of response weights 474 w_si = setups[a].calculate_initial_occupation_numbers( 475 0, False, 476 charge=0, 477 nspins=self.nspins, 478 f_j=w_j) 479 # Basis function coefficients based on density 480 f_si = setups[a].calculate_initial_occupation_numbers( 481 0, False, 482 charge=0, 483 nspins=self.nspins) 484 485 if a in basis_functions.my_atom_indices: 486 self.Dresp_asp[a] = setups[a].initialize_density_matrix(w_si) 487 self.D_asp[a] = setups[a].initialize_density_matrix(f_si) 488 489 f_asi[a] = f_si 490 w_asi[a] = w_si 491 492 self.Drespdist_asp = self.distribute_D_asp(self.Dresp_asp) 493 self.Ddist_asp = self.distribute_D_asp(self.D_asp) 494 nt_sG = self.gd.zeros(self.nspins) 495 basis_functions.add_to_density(nt_sG, f_asi) 496 self.vt_sG = self.gd.zeros(self.nspins) 497 basis_functions.add_to_density(self.vt_sG, w_asi) 498 # Update vt_sG to correspond atomic response potential. This will be 499 # used until occupations and eigenvalues are available. 500 self.vt_sG /= nt_sG + self.damp 501 self.vt_sg = self.finegd.zeros(self.nspins) 502 self.density.distribute_and_interpolate(self.vt_sG, self.vt_sg) 503 504 def get_extra_setup_data(self, extra_data): 505 ae = self.ae 506 njcore = ae.njcore 507 w_ln = self.coefficients.get_coefficients_1d(smooth=True) 508 w_j = [] 509 for w_n in w_ln: 510 for w in w_n: 511 w_j.append(w) 512 extra_data['w_j'] = w_j 513 514 w_j = self.coefficients.get_coefficients_1d() 515 x_g = np.dot(w_j[:njcore], safe_sqr(ae.u_j[:njcore])) 516 x_g[1:] /= ae.r[1:] ** 2 * 4 * np.pi 517 x_g[0] = x_g[1] 518 extra_data['core_response'] = x_g 519 520 # For debugging purposes 521 w_j = self.coefficients.get_coefficients_1d() 522 u2_j = safe_sqr(self.ae.u_j) 523 v_g = self.weight * np.dot(w_j, u2_j) / ( 524 np.dot(self.ae.f_j, u2_j) + self.damp) 525 v_g[0] = v_g[1] 526 extra_data['all_electron_response'] = v_g 527 528 # Calculate Hardness of spherical atom, for debugging purposes 529 l = [np.where(f < 1e-3, e, 1000) 530 for f, e in zip(self.ae.f_j, self.ae.e_j)] 531 h = [np.where(f > 1e-3, e, -1000) 532 for f, e in zip(self.ae.f_j, self.ae.e_j)] 533 lumo_e = min(l) 534 homo_e = max(h) 535 if lumo_e < 999: # If there is unoccpied orbital 536 w_j = self.coefficients.get_coefficients_1d_for_lumo_perturbation() 537 v_g = self.weight * np.dot(w_j, u2_j) / ( 538 np.dot(self.ae.f_j, u2_j) + self.damp) 539 e2 = [e + np.dot(u2 * v_g, self.ae.dr) 540 for u2, e in zip(u2_j, self.ae.e_j)] 541 lumo_2 = min([np.where(f < 1e-3, e, 1000) 542 for f, e in zip(self.ae.f_j, e2)]) 543 # print('New lumo eigenvalue:', lumo_2 * 27.2107) 544 self.hardness = lumo_2 - homo_e 545 # print('Hardness predicted: %10.3f eV' % 546 # (self.hardness * 27.2107)) 547 548 def write(self, writer): 549 """Writes response specific data.""" 550 wfs = self.wfs 551 kpt_comm = wfs.kd.comm 552 gd = wfs.gd 553 554 # Write the pseudodensity on the coarse grid: 555 shape = (wfs.nspins,) + tuple(gd.get_size_of_global_array()) 556 writer.add_array('gllb_pseudo_response_potential', shape) 557 if kpt_comm.rank == 0: 558 for vt_G in self.vt_sG: 559 writer.fill(gd.collect(vt_G) * Ha) 560 561 def pack(X0_asp): 562 X_asp = wfs.setups.empty_atomic_matrix( 563 wfs.nspins, wfs.atom_partition) 564 # XXX some of the provided X0_asp contain strangely duplicated 565 # elements. Take only the minimal set: 566 for a in X_asp: 567 X_asp[a][:] = X0_asp[a] 568 return pack_atomic_matrices(X_asp) 569 570 writer.write(gllb_atomic_density_matrices=pack(self.D_asp)) 571 writer.write(gllb_atomic_response_matrices=pack(self.Dresp_asp)) 572 573 def empty_atomic_matrix(self): 574 assert self.wfs.atom_partition is self.density.atom_partition 575 return self.wfs.setups.empty_atomic_matrix(self.wfs.nspins, 576 self.wfs.atom_partition) 577 578 def read(self, reader): 579 r = reader.hamiltonian.xc 580 wfs = self.wfs 581 582 d('Reading vt_sG') 583 self.gd.distribute(r.gllb_pseudo_response_potential / reader.ha, 584 self.vt_sG) 585 self.density.distribute_and_interpolate(self.vt_sG, self.vt_sg) 586 587 def unpack(D_sP): 588 return unpack_atomic_matrices(D_sP, wfs.setups) 589 590 # Read atomic density matrices and non-local part of hamiltonian: 591 D_asp = unpack(r.gllb_atomic_density_matrices) 592 Dresp_asp = unpack(r.gllb_atomic_response_matrices) 593 594 # All density matrices are loaded to all cores 595 # First distribute them to match density.D_asp 596 self.D_asp = self.empty_atomic_matrix() 597 self.Dresp_asp = self.empty_atomic_matrix() 598 for a in self.density.D_asp: 599 self.D_asp[a][:] = D_asp[a] 600 self.Dresp_asp[a][:] = Dresp_asp[a] 601 assert len(self.Dresp_asp) == len(self.density.D_asp) 602 603 # Further distributes the density matricies for xc-corrections 604 self.Ddist_asp = self.distribute_D_asp(self.D_asp) 605 self.Drespdist_asp = self.distribute_D_asp(self.Dresp_asp) 606 607 def heeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeelp(self, olddens): 608 # XXX This function should be removed once the deprecated 609 # `fixdensity=True` option is removed. 610 from gpaw.density import redistribute_array 611 self.vt_sg = redistribute_array(self.vt_sg, 612 olddens.finegd, self.finegd, 613 self.wfs.nspins, self.wfs.kptband_comm) 614 self.vt_sG = redistribute_array(self.vt_sG, 615 olddens.gd, self.gd, 616 self.wfs.nspins, self.wfs.kptband_comm) 617