Python XC.calculate_paw_correction方法代码示例

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


x = 0.000001
ra.seed(8)
for xc in ['LDA', 'PBE']:
    print(xc)
    xc = XC(xc)
    s = create_setup('N', xc)
    ni = s.ni
    nii = ni * (ni + 1) // 2
    D_p = 0.1 * ra.random(nii) + 0.2
    H_p = np.zeros(nii)

    E = xc.calculate_paw_correction(s, D_p.reshape(1, -1), H_p.reshape(1, -1))
    dD_p = x * ra.random(nii)
    dE = np.dot(H_p, dD_p) / x
    D_p += dD_p
    Ep = xc.calculate_paw_correction(s, D_p.reshape(1, -1))
    D_p -= 2 * dD_p
    Em = xc.calculate_paw_correction(s, D_p.reshape(1, -1))
    print(dE, dE - 0.5 * (Ep - Em) / x)
    equal(dE, 0.5 * (Ep - Em) / x, 1e-6)

    Ems = xc.calculate_paw_correction(s, np.array([0.5 * D_p, 0.5 * D_p]))
    print(Em - Ems)
    equal(Em, Ems, 1.0e-12)

    D_sp = 0.1 * ra.random((2, nii)) + 0.2
    H_sp = np.zeros((2, nii))

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate_paw_correction [as 别名]
    'GGA_X_PBE+GGA_C_PBE', 'GGA_X_PBE_R+GGA_C_PBE',
    'GGA_X_B88+GGA_C_P86', 'GGA_X_B88+GGA_C_LYP',
    'GGA_X_FT97_A+GGA_C_LYP'
    ]

x = 0.000001
for xcname in libxc_set:
    ra.seed(8)
    xc = XC(xcname)
    s = create_setup('N', xc)
    ni = s.ni
    nii = ni * (ni + 1) // 2
    D_p = 0.1 * ra.random(nii) + 0.4
    H_p = np.zeros(nii)

    E1 = xc.calculate_paw_correction(s, D_p.reshape(1, -1), H_p.reshape(1, -1))
    dD_p = x * ra.random(nii)
    D_p += dD_p
    dE = np.dot(H_p, dD_p) / x
    E2 = xc.calculate_paw_correction(s, D_p.reshape(1, -1))
    print xcname, dE, (E2 - E1) / x
    equal(dE, (E2 - E1) / x, 0.003)

    E2s = xc.calculate_paw_correction(s,
        np.array([0.5 * D_p, 0.5 * D_p]), np.array([H_p, H_p]))
    print E2, E2s
    equal(E2, E2s, 1.0e-12)

    if xcname in reference: # compare with old gpaw
        print 'A:', E2, reference[xcname]
        equal(E2, reference[xcname], tolerance)

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

x = 0.0000001
ra.seed(8)
for xc in ['LDA']:#, 'PBE']:
    print(xc)
    xc = XC(xc)
    s = create_setup('N', xc)
    ni = s.ni
    D_sp = np.array([pack(np.outer(P_i, P_i))
                     for P_i in 0.2 * ra.random((2, ni)) - 0.1])
    D_sp[1] += D_sp[0]

    nii = ni * (ni + 1) // 2

    E = xc.calculate_paw_correction(s, D_sp)

    xc = NonCollinearFunctional(xc)

    Dnc_sp = np.zeros((4, nii))
    Dnc_sp[0] = D_sp.sum(0)
    Dnc_sp[3] = D_sp[0] - D_sp[1]
    Enc = xc.calculate_paw_correction(s, Dnc_sp)
    print(E, E-Enc)
    assert abs(E - Enc) < 1e-11

    Dnc_sp[1] = 0.1 * Dnc_sp[3]
    Dnc_sp[2] = 0.2 * Dnc_sp[3]
    Dnc_sp[3] *= (1 - 0.1**2 - 0.2**2)**0.5
    H_sp = 0 * Dnc_sp
    Enc = xc.calculate_paw_correction(s, Dnc_sp, H_sp)

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate_paw_correction [as 别名]

#.........这里部分代码省略.........
                #              stored on each nucleus
                timer2.start('init v grids')
                vp_s=np.zeros(nt_s.shape,nt_s.dtype.char)
                vm_s=np.zeros(nt_s.shape,nt_s.dtype.char)
                if kss.npspins == 2: # spin polarised
                    nv_s = nt_s.copy()
                    nv_s[kss[ij].pspin] += ns * kss[ij].get(fg)
                    xc.calculate(gd, nv_s, vp_s)
                    nv_s = nt_s.copy()
                    nv_s[kss[ij].pspin] -= ns * kss[ij].get(fg)
                    xc.calculate(gd, nv_s, vm_s)
                else: # spin unpolarised
                    nv = nt_s + ns * kss[ij].get(fg)
                    xc.calculate(gd, nv, vp_s)
                    nv = nt_s - ns * kss[ij].get(fg)
                    xc.calculate(gd, nv, vm_s)
                vvt_s = (0.5 / ns) * (vp_s - vm_s)
                timer2.stop()

                # initialize the correction matrices
                timer2.start('init v corrections')
                I_asp = {}
                for a, P_ni in wfs.kpt_u[kss[ij].spin].P_ani.items():
                    # create the modified density matrix
                    Pi_i = P_ni[kss[ij].i]
                    Pj_i = P_ni[kss[ij].j]
                    P_ii = np.outer(Pi_i, Pj_i)
                    # we need the symmetric form, hence we can pack
                    P_p = pack(P_ii)
                    D_sp = self.paw.density.D_asp[a].copy()
                    D_sp[kss[ij].pspin] -= ns * P_p
                    setup = wfs.setups[a]
                    I_sp = np.zeros_like(D_sp)
                    self.xc.calculate_paw_correction(setup, D_sp, I_sp)
                    I_sp *= -1.0
                    D_sp = self.paw.density.D_asp[a].copy()
                    D_sp[kss[ij].pspin] += ns * P_p
                    self.xc.calculate_paw_correction(setup, D_sp, I_sp)
                    I_sp /= 2.0 * ns
                    I_asp[a] = I_sp
                timer2.stop()
                    
            timer.stop()
            t0 = timer.get_time('init')
            timer.start(ij)
            
            for kq in range(ij,nij):
                weight = self.weight_Kijkq(ij, kq)
                
                if self.derivativeLevel == 0:
                    # only Exc is available
                    
                    if kss.npspins==2: # spin polarised
                        nv_g = nt_sg.copy()
                        nv_g[kss[ij].pspin] += kss[ij].get(fg)
                        nv_g[kss[kq].pspin] += kss[kq].get(fg)
                        Excpp = xc.get_energy_and_potential(
                            nv_g[0], v_g, nv_g[1], v_g)
                        nv_g = nt_sg.copy()
                        nv_g[kss[ij].pspin] += kss[ij].get(fg)
                        nv_g[kss[kq].pspin] -= kss[kq].get(fg)
                        Excpm = xc.get_energy_and_potential(\
                                            nv_g[0],v_g,nv_g[1],v_g)
                        nv_g = nt_sg.copy()
                        nv_g[kss[ij].pspin] -=\
                                        kss[ij].get(fg)

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate_paw_correction [as 别名]
    v_g = np.zeros((1, n, n, n))

    P_ni = 0.2 * ra.random((20, ni))
    P_ni[:, nao:] = 0.0
    D_ii = np.dot(np.transpose(P_ni), P_ni)
    D_p = pack(D_ii)
    p = 0
    for i1 in range(nao):
        for i2 in range(i1, nao):
            n_g += D_p[p] * psit_ig[i1] * psit_ig[i2]
            p += 1
        p += ni - nao

    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 = (xc.calculate_paw_correction(s, D_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_paw_correction [as 别名]
    n_g = np.zeros((1, n, n, n))
    v_g = np.zeros((1, n, n, n))

    P_ni = 0.2 * ra.random((20, ni))
    P_ni[:, nao:] = 0.0
    D_ii = np.dot(np.transpose(P_ni), P_ni)
    D_p = pack(D_ii)
    p = 0
    for i1 in range(nao):
        for i2 in range(i1, nao):
            n_g += D_p[p] * psit_ig[i1] * psit_ig[i2]
            p += 1
        p += ni - nao

    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 = xc.calculate_paw_correction(s, D_p.reshape(1, -1)) + s.xc_correction.Exc0

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