Python xc.XC类代码示例

def create_setup(symbol, xc='LDA', lmax=0,
                 type='paw', basis=None, setupdata=None,
                 filter=None, world=None):
    if isinstance(xc, str):
        xc = XC(xc)

    if isinstance(type, str) and ':' in type:
        # Parse DFT+U parameters from type-string:
        # Examples: "type:l,U" or "type:l,U,scale"
        type, lu = type.split(':')
        if type == '':
            type = 'paw'
        l = 'spdf'.find(lu[0])
        assert lu[1] == ','
        U = lu[2:]
        if ',' in U:
            U, scale = U.split(',')
        else:
            scale = True
        U = float(U) / units.Hartree
        scale = int(scale)
    else:
        U = None

    if setupdata is None:
        if type == 'hgh' or type == 'hgh.sc':
            lmax = 0
            from gpaw.hgh import HGHSetupData, setups, sc_setups
            if type == 'hgh.sc':
                table = sc_setups
            else:
                table = setups
            parameters = table[symbol]
            setupdata = HGHSetupData(parameters)
        elif type == 'ah':
            from gpaw.ah import AppelbaumHamann
            ah = AppelbaumHamann()
            ah.build(basis)
            return ah
        elif type == 'ae':
            from gpaw.ae import HydrogenAllElectronSetup
            assert symbol == 'H'
            ae = HydrogenAllElectronSetup()
            ae.build(basis)
            return ae
        elif type == 'ghost':
            from gpaw.lcao.bsse import GhostSetupData
            setupdata = GhostSetupData(symbol)
        else:
            setupdata = SetupData(symbol, xc.get_setup_name(),
                                  type, True,
                                  world=world)
    if hasattr(setupdata, 'build'):
        setup = LeanSetup(setupdata.build(xc, lmax, basis, filter))
        if U is not None:
            setup.set_hubbard_u(U, l, scale)
        return setup
    else:
        return setupdata

def vxc(paw, xc=None, coredensity=True):
    """Calculate XC-contribution to eigenvalues."""
    
    ham = paw.hamiltonian
    dens = paw.density
    wfs = paw.wfs

    if xc is None:
        xc = ham.xc
    elif isinstance(xc, str):
        xc = XC(xc)

    if dens.nt_sg is None:
        dens.interpolate_pseudo_density()

    thisisatest = not True
    
    if xc.orbital_dependent:
        paw.get_xc_difference(xc)

    # Calculate XC-potential:
    vxct_sg = ham.finegd.zeros(wfs.nspins)
    xc.calculate(dens.finegd, dens.nt_sg, vxct_sg)
    vxct_sG = ham.gd.empty(wfs.nspins)
    ham.restrict(vxct_sg, vxct_sG)
    if thisisatest:
        vxct_sG[:] = 1
        
    # ... and PAW corrections:
    dvxc_asii = {}
    for a, D_sp in dens.D_asp.items():
        dvxc_sp = np.zeros_like(D_sp)
        xc.calculate_paw_correction(wfs.setups[a], D_sp, dvxc_sp, a=a, 
                                    addcoredensity=coredensity)
        dvxc_asii[a] = [unpack(dvxc_p) for dvxc_p in dvxc_sp]
        if thisisatest:
            dvxc_asii[a] = [wfs.setups[a].dO_ii]

    vxc_un = np.empty((wfs.kd.mynks, wfs.bd.mynbands))
    for u, vxc_n in enumerate(vxc_un):
        kpt = wfs.kpt_u[u]
        vxct_G = vxct_sG[kpt.s]
        for n in range(wfs.bd.mynbands):
            psit_G = wfs._get_wave_function_array(u, n, realspace=True)
            vxc_n[n] = wfs.gd.integrate((psit_G * psit_G.conj()).real,
                                        vxct_G, global_integral=False)

        for a, dvxc_sii in dvxc_asii.items():
            P_ni = kpt.P_ani[a]
            vxc_n += (np.dot(P_ni, dvxc_sii[kpt.s]) *
                      P_ni.conj()).sum(1).real

    wfs.gd.comm.sum(vxc_un)
    vxc_skn = wfs.kd.collect(vxc_un)

    if xc.orbital_dependent:
        vxc_skn += xc.exx_skn

    return vxc_skn * Hartree

 def get_xc_difference(self, xc):
     if isinstance(xc, str):
         xc = XC(xc)
     xc.initialize(self.density, self.hamiltonian, self.wfs,
                   self.occupations)
     xc.set_positions(self.atoms.get_scaled_positions() % 1.0)
     if xc.orbital_dependent:
         self.converge_wave_functions()
     return self.hamiltonian.get_xc_difference(xc, self.density) * Hartree

 def linearize_to_xc(self, newxc):
     """Linearize Hamiltonian to difference XC functional.
     
     Used in real time TDDFT to perform calculations with various kernels.
     """
     if isinstance(newxc, str):
         newxc = XC(newxc)
     self.txt.write('Linearizing xc-hamiltonian to ' + str(newxc))
     newxc.initialize(self.density, self.hamiltonian, self.wfs,
                      self.occupations)
     self.hamiltonian.linearize_to_xc(newxc, self.density)

def paired():
    xc = XC('vdW-DF')
    n = 0.3 * np.ones((1, N, N, N))
    n += 0.01 * np.cos(np.arange(N) * 2 * pi / N)
    v = 0.0 * n
    xc.calculate(gd, n, v)
    n2 = 1.0 * n
    i = 1
    n2[0, i, i, i] += 0.00002
    x = v[0, i, i, i] * gd.dv
    E2 = xc.calculate(gd, n2, v)
    n2[0, i, i, i] -= 0.00004
    E2 -= xc.calculate(gd, n2, v)
    x2 = E2 / 0.00004
    print(i, x, x2, x - x2, x / x2)
    equal(x, x2, 2e-11)

    def initialize(self):
        self.occupations = self.nlfunc.occupations
        self.xc = XC(self.functional)

        # Always 1 spin, no matter what calculation nspins is
        self.vt_sg = self.nlfunc.finegd.empty(1) 
        self.e_g = self.nlfunc.finegd.empty()#.ravel()

def polarized():
    xc = XC('vdW-DF')
    n = 0.04 * np.ones((2, N, N, N))
    n[1] = 0.3
    n[0] += 0.02 * np.sin(np.arange(N) * 2 * pi / N)
    n[1] += 0.2 * np.cos(np.arange(N) * 2 * pi / N)
    v = 0.0 * n
    xc.calculate(gd, n, v)
    n2 = 1.0 * n
    i = 1
    n2[0, i, i, i] += 0.00002
    x = v[0, i, i, i] * gd.dv
    E2 = xc.calculate(gd, n2, v)
    n2[0, i, i, i] -= 0.00004
    E2 -= xc.calculate(gd, n2, v)
    x2 = E2 / 0.00004
    print(i, x, x2, x - x2, x / x2)
    equal(x, x2, 1e-10)

    def __init__(self, symbol, xc='LDA', spinpol=False, dirac=False,
                 log=sys.stdout):
        """All-electron calculation for spherically symmetric atom.

        symbol: str (or int)
            Chemical symbol (or atomic number).
        xc: str
            Name of XC-functional.
        spinpol: bool
            If true, do spin-polarized calculation.  Default is spin-paired.
        dirac: bool
            Solve Dirac equation instead of Schrödinger equation.
        log: stream
            Text output."""

        if isinstance(symbol, int):
            symbol = chemical_symbols[symbol]
        self.symbol = symbol
        self.Z = atomic_numbers[symbol]

        self.nspins = 1 + int(bool(spinpol))

        self.dirac = bool(dirac)
        
        self.scalar_relativistic = False

        if isinstance(xc, str):
            self.xc = XC(xc)
        else:
            self.xc = xc

        if log is None:
            log = devnull
        self.fd = log

        self.vr_sg = None  # potential * r
        self.n_sg = 0.0    # density
        self.rgd = None     # radial grid descriptor

        # Energies:
        self.ekin = None
        self.eeig = None
        self.eH = None
        self.eZ = None

        self.channels = None

        self.initialize_configuration()

        self.log('Z:              ', self.Z)
        self.log('Name:           ', atomic_names[self.Z])
        self.log('Symbol:         ', symbol)
        self.log('XC-functional:  ', self.xc.name)
        self.log('Equation:       ', ['Schrödinger', 'Dirac'][self.dirac])

        self.method = 'Gaussian basis-set'

 def initialize(self):
     self.xc = XC(self.functional)
     self.vt_sg = self.nlfunc.finegd.empty(self.nlfunc.nspins)
     self.e_g = self.nlfunc.finegd.empty()

def calculate_Kxc(pd, nt_sG, R_av, setups, D_asp, functional='ALDA',
                  density_cut=None):
    """ALDA kernel"""

    gd = pd.gd
    npw = pd.ngmax
    nG = pd.gd.N_c
    vol = pd.gd.volume
    bcell_cv = np.linalg.inv(pd.gd.cell_cv)
    G_Gv = pd.get_reciprocal_vectors()
    
    # The soft part
    #assert np.abs(nt_sG[0].shape - nG).sum() == 0
    if functional == 'ALDA_X':
        x_only = True
        A_x = -3. / 4. * (3. / np.pi)**(1. / 3.)
        nspins = len(nt_sG)
        assert nspins in [1, 2]
        fxc_sg = nspins**(1. / 3.) * 4. / 9. * A_x * nt_sG**(-2. / 3.)
    else:
        assert len(nt_sG) == 1
        x_only = False
        fxc_sg = np.zeros_like(nt_sG)
        xc = XC(functional[1:])
        xc.calculate_fxc(gd, nt_sG, fxc_sg)

    if density_cut is not None:
        fxc_sg[np.where(nt_sG * len(nt_sG) < density_cut)] = 0.0

    # FFT fxc(r)
    nG0 = nG[0] * nG[1] * nG[2]
    tmp_sg = [np.fft.fftn(fxc_sg[s]) * vol / nG0 for s in range(len(nt_sG))]

    r_vg = gd.get_grid_point_coordinates()
    Kxc_sGG = np.zeros((len(fxc_sg), npw, npw), dtype=complex)
    for s in range(len(fxc_sg)):
        for iG, iQ in enumerate(pd.Q_qG[0]):
            iQ_c = (np.unravel_index(iQ, nG) + nG // 2) % nG - nG // 2
            for jG, jQ in enumerate(pd.Q_qG[0]):
                jQ_c = (np.unravel_index(jQ, nG) + nG // 2) % nG - nG // 2
                ijQ_c = (iQ_c - jQ_c)
                if (abs(ijQ_c) < nG // 2).all():
                    Kxc_sGG[s, iG, jG] = tmp_sg[s][tuple(ijQ_c)]

    # The PAW part
    KxcPAW_sGG = np.zeros_like(Kxc_sGG)
    dG_GGv = np.zeros((npw, npw, 3))
    for v in range(3):
        dG_GGv[:, :, v] = np.subtract.outer(G_Gv[:, v], G_Gv[:, v])

    for a, setup in enumerate(setups):
        if rank == a % size:
            rgd = setup.xc_correction.rgd
            n_qg = setup.xc_correction.n_qg
            nt_qg = setup.xc_correction.nt_qg
            nc_g = setup.xc_correction.nc_g
            nct_g = setup.xc_correction.nct_g
            Y_nL = setup.xc_correction.Y_nL
            dv_g = rgd.dv_g

            D_sp = D_asp[a]
            B_pqL = setup.xc_correction.B_pqL
            D_sLq = np.inner(D_sp, B_pqL.T)
            nspins = len(D_sp)

            f_sg = rgd.empty(nspins)
            ft_sg = rgd.empty(nspins)

            n_sLg = np.dot(D_sLq, n_qg)
            nt_sLg = np.dot(D_sLq, nt_qg)

            # Add core density
            n_sLg[:, 0] += np.sqrt(4. * np.pi) / nspins * nc_g
            nt_sLg[:, 0] += np.sqrt(4. * np.pi) / nspins * nct_g

            coefatoms_GG = np.exp(-1j * np.inner(dG_GGv, R_av[a]))
            for n, Y_L in enumerate(Y_nL):
                w = weight_n[n]
                f_sg[:] = 0.0
                n_sg = np.dot(Y_L, n_sLg)
                if x_only:
                    f_sg = nspins * (4 / 9.) * A_x * (nspins * n_sg)**(-2 / 3.)
                else:
                    xc.calculate_fxc(rgd, n_sg, f_sg)

                ft_sg[:] = 0.0
                nt_sg = np.dot(Y_L, nt_sLg)
                if x_only:
                    ft_sg = nspins * (4 / 9.) * (A_x
                                                 * (nspins * nt_sg)**(-2 / 3.))
                else:
                    xc.calculate_fxc(rgd, nt_sg, ft_sg)
                for i in range(len(rgd.r_g)):
                    coef_GG = np.exp(-1j * np.inner(dG_GGv, R_nv[n])
                                     * rgd.r_g[i])
                    for s in range(len(f_sg)):
                        KxcPAW_sGG[s] += w * np.dot(coef_GG,
                                                    (f_sg[s, i] -
                                                     ft_sg[s, i])
                                                    * dv_g[i]) * coefatoms_GG
#.........这里部分代码省略.........

    def run(self, use_restart_file=True):
        #     beta g
        # r = ------, g = 0, 1, ..., N - 1
        #     N - g
        #
        #        rN
        # g = --------
        #     beta + r

        t = self.text
        N = self.N
        beta = self.beta
        t(N, 'radial gridpoints.')
        self.rgd = AERadialGridDescriptor(beta / N, 1.0 / N, N)
        g = np.arange(N, dtype=float)
        self.r = self.rgd.r_g
        self.dr = self.rgd.dr_g
        self.d2gdr2 = self.rgd.d2gdr2()

        # Number of orbitals:
        nj = len(self.n_j)

        # Radial wave functions multiplied by radius:
        self.u_j = np.zeros((nj, self.N))

        # Effective potential multiplied by radius:
        self.vr = np.zeros(N)

        # Electron density:
        self.n = np.zeros(N)

        # Always spinpaired nspins=1
        self.xc = XC(self.xcname)

        # Initialize for non-local functionals
        if self.xc.type == 'GLLB':
            self.xc.pass_stuff_1d(self)
            self.xc.initialize_1d()
            
        n_j = self.n_j
        l_j = self.l_j
        f_j = self.f_j
        e_j = self.e_j
        
        Z = self.Z    # nuclear charge
        r = self.r    # radial coordinate
        dr = self.dr  # dr/dg
        n = self.n    # electron density
        vr = self.vr  # effective potential multiplied by r

        vHr = np.zeros(self.N)
        self.vXC = np.zeros(self.N)

        restartfile = '%s/%s.restart' % (tempdir, self.symbol)
        if self.xc.type == 'GLLB' or not use_restart_file:
            # Do not start from initial guess when doing
            # non local XC!
            # This is because we need wavefunctions as well
            # on the first iteration.
            fd = None
        else:
            try:
                fd = open(restartfile, 'r')
            except IOError:
                fd = None
            else:
                try:
                    n[:] = pickle.load(fd)
                except (ValueError, IndexError):
                    fd = None
                else:
                    norm = np.dot(n * r**2, dr) * 4 * pi
                    if abs(norm - sum(f_j)) > 0.01:
                        fd = None
                    else:
                        t('Using old density for initial guess.')
                        n *= sum(f_j) / norm

        if fd is None:
            self.initialize_wave_functions()
            n[:] = self.calculate_density()

        bar = '|------------------------------------------------|'
        t(bar)
        
        niter = 0
        nitermax = 117
        qOK = log(1e-10)
        mix = 0.4
        
        # orbital_free needs more iterations and coefficient
        if self.orbital_free:
            #qOK = log(1e-14)
            e_j[0] /= self.tf_coeff
            mix = 0.01
            nitermax = 1000
            
        vrold = None
        
        while True:
#.........这里部分代码省略.........

class C_GLLBScr(Contribution):
    def __init__(self, nlfunc, weight, functional='GGA_X_B88', metallic=False):
        Contribution.__init__(self, nlfunc, weight)
        self.functional = functional
        self.old_coeffs = None
        self.iter = 0
        self.metallic = metallic
        
    def get_name(self):
        return 'SCREENING'

    def get_desc(self):
        return '(' + self.functional + ')'
        
    # Initialize GLLBScr functional
    def initialize_1d(self):
        self.ae = self.nlfunc.ae
        self.xc = XC(self.functional)
        self.v_g = np.zeros(self.ae.N)
        self.e_g = np.zeros(self.ae.N)

    # Calcualte the GLLB potential and energy 1d
    def add_xc_potential_and_energy_1d(self, v_g):
        self.v_g[:] = 0.0
        self.e_g[:] = 0.0
        self.xc.calculate_spherical(self.ae.rgd, self.ae.n.reshape((1, -1)),
                                    self.v_g.reshape((1, -1)), self.e_g)
        v_g += 2 * self.weight * self.e_g / (self.ae.n + 1e-10)
        Exc = self.weight * np.sum(self.e_g * self.ae.rgd.dv_g)
        return Exc

    def initialize(self):
        self.occupations = self.nlfunc.occupations
        self.xc = XC(self.functional)

        # Always 1 spin, no matter what calculation nspins is
        self.vt_sg = self.nlfunc.finegd.empty(1) 
        self.e_g = self.nlfunc.finegd.empty()#.ravel()

    def get_coefficient_calculator(self):
        return self

    def f(self, f):
        return sqrt(f)
    
    def get_coefficients(self, e_j, f_j):
        homo_e = max( [ np.where(f>1e-3, e, -1000) for f,e in zip(f_j, e_j)] ) 
        return [ f * K_G * self.f( max(0, homo_e - e)) for e,f in zip(e_j, f_j) ]

    def get_coefficients_1d(self, smooth=False, lumo_perturbation = False):
        homo_e = max( [ np.where(f>1e-3, e, -1000) for f,e in zip(self.ae.f_j, self.ae.e_j)]) 
        if not smooth:
            if lumo_perturbation:
                lumo_e = min( [ np.where(f<1e-3, e, 1000) for f,e in zip(self.ae.f_j, self.ae.e_j)])
                return np.array([ f * K_G * (self.f( max(0, lumo_e - e)) - self.f(max(0, homo_e -e)))
                                        for e,f in zip(self.ae.e_j, self.ae.f_j) ])
            else:
                return np.array([ f * K_G * (self.f( max(0, homo_e - e)))
                                   for e,f in zip(self.ae.e_j, self.ae.f_j) ])
        else:
            return [ [ f * K_G * self.f( max(0, homo_e - e))
                    for e,f in zip(e_n, f_n) ]
                     for e_n, f_n in zip(self.ae.e_ln, self.ae.f_ln) ]
        

    def get_coefficients_by_kpt(self, kpt_u, lumo_perturbation=False, homolumo=None, nspins=1):
        if not hasattr(kpt_u[0],'orbitals_ready'):
            kpt_u[0].orbitals_ready = True
            return None
        #if kpt_u[0].psit_nG is None or isinstance(kpt_u[0].psit_nG,
        #                                          TarFileReference): 
        #    if kpt_u[0].C_nM==None:
        #        return None

        if homolumo == None:
            if self.metallic:
                # For metallic systems, the calculated fermi level represents 
                # the most accurate estimate for reference-energy
                eref_lumo_s = eref_s = nspins * [ self.occupations.get_fermi_level() ]
            else:
                # Find homo and lumo levels for each spin
                eref_s = []
                eref_lumo_s = []
                for s in range(nspins):
                    homo, lumo = self.occupations.get_homo_lumo_by_spin(self.nlfunc.wfs, s)
                    eref_s.append(homo)
                    eref_lumo_s.append(lumo)
        else:
            eref_s, eref_lumo_s = homolumo
            if not isinstance(eref_s, (list, tuple)):
                eref_s = [ eref_s ]
                eref_lumo_s = [ eref_lumo_s ]

        # The parameter ee might sometimes be set to small thereshold value to
        # achieve convergence on small systems with degenerate HOMO.
        if len(kpt_u) > nspins:
            ee = 0.0
        else:
            ee = 0.05 / 27.21

#.........这里部分代码省略.........

            gd.beg_c[2] <= 3 < gd.end_c[2])
    if here:
        x = v[-1, 1, 2, 3] * gd.dv
        n[-1, 1, 2, 3] += 0.000001
    Ep = xc.calculate(gd, n, v)
    if here:
        n[-1, 1, 2, 3] -= 0.000002
    Em = xc.calculate(gd, n, v)
    x2 = (Ep - Em) / 0.000002
    if here:
        print(xc.name, E, x, x2, x - x2)
        equal(x, x2, 1e-11)
        n[-1, 1, 2, 3] += 0.000001

    if 0:#xc.type == 'LDA':
        xc = XC(NonCollinearLDAKernel())
    else:
        xc = NonCollinearFunctional(xc)

    n2 = gd.zeros(4)
    n2[0] = n.sum(0)
    n2[3] = n[0] - n[1]
    E2 = xc.calculate(gd, n2)
    print(E, E2-E)
    assert abs(E2 - E) < 1e-11
    n2[1] = 0.1 * n2[3]
    n2[2] = 0.2 * n2[3]
    n2[3] *= (1 - 0.1**2 - 0.2**2)**0.5
    v = n2 * 0
    E2 = xc.calculate(gd, n2, v)
    print(E, E2-E)

def calculate_Kxc(gd,
                  nt_sG,
                  npw,
                  Gvec_Gc,
                  nG,
                  vol,
                  bcell_cv,
                  R_av,
                  setups,
                  D_asp,
                  functional='ALDA',
                  density_cut=None):
    """ALDA kernel"""

    # The soft part
    #assert np.abs(nt_sG[0].shape - nG).sum() == 0
    if functional == 'ALDA_X':
        x_only = True
        A_x = -3. / 4. * (3. / np.pi)**(1. / 3.)
        nspins = len(nt_sG)
        assert nspins in [1, 2]
        fxc_sg = nspins**(1. / 3.) * 4. / 9. * A_x * nt_sG**(-2. / 3.)
    else:
        assert len(nt_sG) == 1
        x_only = False
        fxc_sg = np.zeros_like(nt_sG)
        xc = XC(functional[1:])
        xc.calculate_fxc(gd, nt_sG, fxc_sg)

    if density_cut is not None:
        fxc_sg[np.where(nt_sG * len(nt_sG) < density_cut)] = 0.0

    # FFT fxc(r)
    nG0 = nG[0] * nG[1] * nG[2]
    tmp_sg = [np.fft.fftn(fxc_sg[s]) * vol / nG0 for s in range(len(nt_sG))]

    r_vg = gd.get_grid_point_coordinates()
    Kxc_sGG = np.zeros((len(fxc_sg), npw, npw), dtype=complex)
    for s in range(len(fxc_sg)):
        for iG in range(npw):
            for jG in range(npw):
                dG_c = Gvec_Gc[iG] - Gvec_Gc[jG]
                if (nG / 2 - np.abs(dG_c) > 0).all():
                    index = (dG_c + nG) % nG
                    Kxc_sGG[s, iG, jG] = tmp_sg[s][index[0],
                                                   index[1],
                                                   index[2]]
                else:  # not in the fft index
                    dG_v = np.dot(dG_c, bcell_cv)
                    dGr_g = gemmdot(dG_v, r_vg, beta=0.0)
                    Kxc_sGG[s, iG, jG] = gd.integrate(np.exp(-1j * dGr_g)
                                                      * fxc_sg[s])

    # The PAW part
    KxcPAW_sGG = np.zeros_like(Kxc_sGG)
    dG_GGv = np.zeros((npw, npw, 3))
    for iG in range(npw):
        for jG in range(npw):
            dG_c = Gvec_Gc[iG] - Gvec_Gc[jG]
            dG_GGv[iG, jG] = np.dot(dG_c, bcell_cv)

    for a, setup in enumerate(setups):
        if rank == a % size:
            rgd = setup.xc_correction.rgd
            n_qg = setup.xc_correction.n_qg
            nt_qg = setup.xc_correction.nt_qg
            nc_g = setup.xc_correction.nc_g
            nct_g = setup.xc_correction.nct_g
            Y_nL = setup.xc_correction.Y_nL
            dv_g = rgd.dv_g

            D_sp = D_asp[a]
            B_pqL = setup.xc_correction.B_pqL
            D_sLq = np.inner(D_sp, B_pqL.T)
            nspins = len(D_sp)

            f_sg = rgd.empty(nspins)
            ft_sg = rgd.empty(nspins)

            n_sLg = np.dot(D_sLq, n_qg)
            nt_sLg = np.dot(D_sLq, nt_qg)

            # Add core density
            n_sLg[:, 0] += np.sqrt(4. * np.pi) / nspins * nc_g
            nt_sLg[:, 0] += np.sqrt(4. * np.pi) / nspins * nct_g

            coefatoms_GG = np.exp(-1j * np.inner(dG_GGv, R_av[a]))
            for n, Y_L in enumerate(Y_nL):
                w = weight_n[n]
                f_sg[:] = 0.0
                n_sg = np.dot(Y_L, n_sLg)
                if x_only:
                    f_sg = nspins * (4 / 9.) * A_x * (nspins * n_sg)**(-2 / 3.)
                else:
                    xc.calculate_fxc(rgd, n_sg, f_sg)

                ft_sg[:] = 0.0
                nt_sg = np.dot(Y_L, nt_sLg)
                if x_only:
                    ft_sg = nspins * (4 / 9.) * (A_x
#.........这里部分代码省略.........

    def __init__(self,
                 calculator=None,
                 kss=None,
                 xc=None,
                 derivativeLevel=None,
                 numscale=0.001,
                 filehandle=None,
                 txt=None,
                 finegrid=2,
                 eh_comm=None,
                 ):
        
        if not txt and calculator:
            txt = calculator.txt
        self.txt, firsttime = initialize_text_stream(txt, mpi.rank)

        if eh_comm == None:
            eh_comm = mpi.serial_comm

        self.eh_comm = eh_comm

        if filehandle is not None:
            self.kss = kss
            self.read(fh=filehandle)
            return None

        self.fullkss = kss
        self.finegrid = finegrid

        if calculator is None:
            return

        self.paw = calculator
        wfs = self.paw.wfs
        
        # handle different grid possibilities
        self.restrict = None
        self.poisson = PoissonSolver(nn=self.paw.hamiltonian.poisson.nn)
        if finegrid:
            self.poisson.set_grid_descriptor(self.paw.density.finegd)
            self.poisson.initialize()
            
            self.gd = self.paw.density.finegd
            if finegrid == 1:
                self.gd = wfs.gd
        else:
            self.poisson.set_grid_descriptor(wfs.gd)
            self.poisson.initialize()
            self.gd = wfs.gd
        self.restrict = Transformer(self.paw.density.finegd, wfs.gd,
                                    self.paw.input_parameters.stencils[1]
                                    ).apply

        if xc == 'RPA': 
            xc = None # enable RPA as keyword
        if xc is not None:
            self.xc = XC(xc)
            self.xc.initialize(self.paw.density, self.paw.hamiltonian,
                               wfs, self.paw.occupations)

            # check derivativeLevel
            if derivativeLevel is None:
                derivativeLevel= \
                    self.xc.get_functional().get_max_derivative_level()
            self.derivativeLevel = derivativeLevel
            # change the setup xc functional if needed
            # the ground state calculation may have used another xc
            if kss.npspins > kss.nvspins:
                spin_increased = True
            else:
                spin_increased = False
        else:
            self.xc = None

        self.numscale = numscale
    
        self.singletsinglet = False
        if kss.nvspins<2 and kss.npspins<2:
             # this will be a singlet to singlet calculation only
             self.singletsinglet=True

        nij = len(kss)
        self.Om = np.zeros((nij,nij))
        self.get_full()