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