Python XC.calculate方法代码示例

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate [as 别名]
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

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate [as 别名]
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)

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate [as 别名]
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)

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate [as 别名]
a = 1.0
gd = GridDescriptor((N, N, N), (a, a, a))

for name in ['LDA', 'PBE']:
    xc = XC(name)
    for nspins in [1, 2]:
        n = gd.empty(nspins)
        n.fill(0.03)
        z = np.arange(gd.beg_c[2], gd.end_c[2]) * a / N
        n[:] += 0.01 * np.sin(2 * pi * z / a)
        if nspins == 2:
            n[1] += 0.01 * np.cos(2 * pi * z / a)
        n /= nspins

        v = 0.0 * n
        E = xc.calculate(gd, n, v)

        here = (gd.beg_c[0] <= 1 < gd.end_c[0] and
                gd.beg_c[1] <= 2 < gd.end_c[1] and
                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(name, nspins, E, x, x2, x - x2)
            equal(x, x2, 1e-11)

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate [as 别名]
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

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

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate [as 别名]
    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)
    assert abs(E2 - E) < 1e-11

    for i in range(4):
        if here:
            x = v[i, 1, 2, 3] * gd.dv
            n2[i, 1, 2, 3] += 0.000001
        Ep = xc.calculate(gd, n2)

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate [as 别名]
class C_XC(Contribution):
    def __init__(self, nlfunc, weight, functional = 'LDA'):
        Contribution.__init__(self, nlfunc, weight)
        self.functional = functional

    def get_name(self):
        return 'XC'

    def get_desc(self):
        return "("+self.functional+")"
        
    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 initialize_1d(self):
        self.ae = self.nlfunc.ae
        self.xc = XC(self.functional) 
        self.v_g = np.zeros(self.ae.N)

    def calculate_spinpaired(self, e_g, n_g, v_g):
        self.e_g[:] = 0.0
        self.vt_sg[:] = 0.0
        self.xc.calculate(self.nlfunc.finegd, n_g[None, ...], self.vt_sg,
                          self.e_g)
        v_g += self.weight * self.vt_sg[0]
        e_g += self.weight * self.e_g

    def calculate_spinpolarized(self, e_g, n_sg, v_sg):
        self.e_g[:] = 0.0
        self.vt_sg[:] = 0.0
        self.xc.calculate(self.nlfunc.finegd, n_sg, self.vt_sg, self.e_g)
        #self.xc.get_energy_and_potential(na_g, self.vt_sg[0], nb_g, self.vt_sg[1], e_g=self.e_g)
        v_sg[0] += self.weight * self.vt_sg[0]
        v_sg[1] += self.weight * self.vt_sg[1]
        e_g += self.weight * self.e_g

    def calculate_energy_and_derivatives(self, setup, D_sp, H_sp, a, addcoredensity=True):
        E = self.xc.calculate_paw_correction(setup, D_sp, H_sp, True, a)
        E += setup.xc_correction.Exc0
        print("E", E)
        return E

    def add_xc_potential_and_energy_1d(self, v_g):
        self.v_g[:] = 0.0
        Exc = self.xc.calculate_spherical(self.ae.rgd,
                                          self.ae.n.reshape((1, -1)),
                                          self.v_g.reshape((1, -1)))
        v_g += self.weight * self.v_g
        return self.weight * Exc

    def add_smooth_xc_potential_and_energy_1d(self, vt_g):
        self.v_g[:] = 0.0
        Exc = self.xc.calculate_spherical(self.ae.rgd,
                                          self.ae.nt.reshape((1, -1)),
                                          self.v_g.reshape((1, -1)))
        vt_g += self.weight * self.v_g
        return self.weight * Exc

    def initialize_from_atomic_orbitals(self, basis_functions):
        # LDA needs only density, which is already initialized
        pass

    def add_extra_setup_data(self, dict):
        # LDA has not any special data
        pass

    def write(self, writer, natoms):
        # LDA has not any special data to be written
        pass

    def read(self, reader):
        # LDA has not any special data to be read
        pass

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate [as 别名]
    P_ni = 0.2 * ra.random((20, ni))
    P_ni[:, niAO:] = 0.0
    D_ii = np.dot(np.transpose(P_ni), P_ni)
    D_p = pack(D_ii)
    p = 0
    for i1 in range(niAO):
        for i2 in range(i1, niAO):
            n_g += D_p[p] * psit_ig[i1] * psit_ig[i2]
            p += 1
        p += ni - niAO

    p = create_localized_functions([s.nct], gd, (0.5, 0.5, 0.5))
    p.add(n_g[0], np.ones(1))
    e_g = gd.zeros()
    xc.calculate(gd, n_g, v_g, e_g)

    r2_g = np.sum((np.indices((n, n, n)) - n / 2)**2, axis=0)
    dv_g = gd.dv * np.less(r2_g, (rcut / a * n)**2)

    E2 = -np.dot(e_g.ravel(), dv_g.ravel())

    s.xc_correction.n_qg[:] = 0.0
    s.xc_correction.nc_g[:] = 0.0
    E1 = (s.xc_correction.calculate(xc, D_p.reshape(1, -1),
                                    H_p.reshape(1, -1))
          + s.xc_correction.Exc0)

    print name, E1, E2, E1 - E2
    equal(E1, E2, 0.0013)

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate [as 别名]
class C_XC(Contribution):
    def __init__(self, nlfunc, weight, functional = 'LDA'):
        Contribution.__init__(self, nlfunc, weight)
        self.functional = functional

    def get_name(self):
        return 'XC'

    def get_desc(self):
        return "("+self.functional+")"
        
    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 initialize_1d(self):
        self.ae = self.nlfunc.ae
        self.xc = XC(self.functional) 
        self.v_g = np.zeros(self.ae.N)

    def calculate_spinpaired(self, e_g, n_g, v_g):
        self.e_g[:] = 0.0
        self.vt_sg[:] = 0.0
        self.xc.calculate(self.nlfunc.finegd, n_g[None, ...], self.vt_sg,
                          self.e_g)
        v_g += self.weight * self.vt_sg[0]
        e_g += self.weight * self.e_g

    def calculate_spinpolarized(self, e_g, na_g, va_g, nb_g, vb_g):
        self.e_g[:] = 0.0
        self.vt_sg[:] = 0.0
        self.xc.get_energy_and_potential(na_g, self.vt_sg[0], nb_g, self.vt_sg[1], e_g=self.e_g)
        va_g += self.weight * self.vt_sg[0]
        vb_g += self.weight * self.vt_sg[1]
        e_g += (self.weight * self.e_g).ravel()

    def calculate_energy_and_derivatives(self, D_sp, H_sp, a):
        # Get the XC-correction instance
        c = self.nlfunc.setups[a].xc_correction

        assert self.nlfunc.nspins == 1
        D_p = D_sp[0]
        dEdD_p = H_sp[0][:]
        D_Lq = dot3(c.B_pqL.T, D_p)
        n_Lg = np.dot(D_Lq, c.n_qg)
        n_Lg[0] += c.nc_g * sqrt(4 * pi)
        nt_Lg = np.dot(D_Lq, c.nt_qg)
        nt_Lg[0] += c.nct_g * sqrt(4 * pi)
        dndr_Lg = np.zeros((c.Lmax, c.ng))
        dntdr_Lg = np.zeros((c.Lmax, c.ng))
        for L in range(c.Lmax):
            c.rgd.derivative(n_Lg[L], dndr_Lg[L])
            c.rgd.derivative(nt_Lg[L], dntdr_Lg[L])
                                                            
        E = 0
        vt_g = np.zeros(c.ng)
        v_g = np.zeros(c.ng)
        e_g = np.zeros(c.ng)
        y = 0
        for w, Y_L in zip(weight_n, c.Y_nL):
            A_Li = rnablaY_nLv[y, :c.Lmax]
            a1x_g = np.dot(A_Li[:, 0], n_Lg)
            a1y_g = np.dot(A_Li[:, 1], n_Lg)
            a1z_g = np.dot(A_Li[:, 2], n_Lg)
            a2_g = a1x_g**2 + a1y_g**2 + a1z_g**2
            a2_g[1:] /= c.rgd.r_g[1:]**2
            a2_g[0] = a2_g[1]
            a1_g = np.dot(Y_L, dndr_Lg)
            a2_g += a1_g**2
            deda2_g = np.zeros(c.ng)  

            v_g[:] = 0.0
            e_g[:] = 0.0
            n_g = np.dot(Y_L, n_Lg)
            self.xc.kernel.calculate(e_g, n_g.reshape((1, -1)),
                                     v_g.reshape((1, -1)),
                                     a2_g.reshape((1, -1)),
                                     deda2_g.reshape((1, -1)))
            
            E += w * np.dot(e_g, c.rgd.dv_g)
            x_g = -2.0 * deda2_g * c.rgd.dv_g * a1_g
            c.rgd.derivative2(x_g, x_g)
            x_g += v_g * c.rgd.dv_g
            dEdD_p += self.weight * w * np.dot(dot3(c.B_pqL, Y_L),
                                  np.dot(c.n_qg, x_g))
            x_g = 8.0 * pi * deda2_g * c.rgd.dr_g
            dEdD_p += w * np.dot(dot3(c.B_pqL,
                                       A_Li[:, 0]),
                                  np.dot(c.n_qg, x_g * a1x_g))
            dEdD_p += w * np.dot(dot3(c.B_pqL,
                                       A_Li[:, 1]),
                                  np.dot(c.n_qg, x_g * a1y_g))
            dEdD_p += w * np.dot(dot3(c.B_pqL,
                                       A_Li[:, 2]),
                                  np.dot(c.n_qg, x_g * a1z_g))

            n_g = np.dot(Y_L, nt_Lg)
            a1x_g = np.dot(A_Li[:, 0], nt_Lg)
            a1y_g = np.dot(A_Li[:, 1], nt_Lg)
#.........这里部分代码省略.........

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate [as 别名]
class C_GLLBScr(Contribution):
    def __init__(self, nlfunc, weight, functional='GGA_X_B88'):
        Contribution.__init__(self, nlfunc, weight)
        self.functional = functional
        self.old_coeffs = None
        self.iter = 0
        
    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)
        self.vt_sg = self.nlfunc.finegd.empty(self.nlfunc.nspins)
        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_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):
        if kpt_u[0].psit_nG is None or isinstance(kpt_u[0].psit_nG,
                                                  TarFileReference): 
            return None

        if homolumo == None:
            e_ref, e_ref_lumo = self.occupations.get_homo_lumo(self.nlfunc.wfs)
        else:
            e_ref, e_ref_lumo = homolumo

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

        if lumo_perturbation:
            return [np.array([
                f * K_G * (self.f( np.where(e_ref_lumo - e>ee, e_ref_lumo-e,0))
                         -self.f( np.where(e_ref      - e>ee, e_ref-e,0)))
                     for e, f in zip(kpt.eps_n, kpt.f_n) ])
                     for kpt in kpt_u ]
            
            
        else:
            coeff = [ np.array([ f * K_G * self.f( np.where(e_ref - e>ee, e_ref-e,0))
                     for e, f in zip(kpt.eps_n, kpt.f_n) ])
                     for kpt in kpt_u ]
            if self.old_coeffs is None:
                self.old_coeffs = coeff
            else:
                # Mix the coefficients with 25%
                mix = 0.25
                self.old_coeffs = [ (1-mix) * old + mix * new for new, old in zip(coeff, self.old_coeffs) ]
            return self.old_coeffs
        

    def calculate_spinpaired(self, e_g, n_g, v_g):
        self.e_g[:] = 0.0
        self.vt_sg[:] = 0.0
        self.xc.calculate(self.nlfunc.finegd, n_g[None, ...], self.vt_sg,
#.........这里部分代码省略.........