Python XC.calculate_spherical方法代码示例

# 需要导入模块: from gpaw.xc import XC [as 别名]
# 或者: from gpaw.xc.XC import calculate_spherical [as 别名]
from gpaw.atom.radialgd import EquidistantRadialGridDescriptor
from gpaw.grid_descriptor import GridDescriptor
from gpaw.xc import XC
import numpy as np
from gpaw.test import equal

rgd = EquidistantRadialGridDescriptor(0.01, 100)

for name in ['LDA', 'PBE']:
    xc = XC(name)
    for nspins in [1, 2]:
        n = rgd.zeros(nspins)
        v = rgd.zeros(nspins)
        n[:] = np.exp(-rgd.r_g**2)
        n[-1] *= 2
        E = xc.calculate_spherical(rgd, n, v)
        i = 23
        x = v[-1, i] * rgd.dv_g[i]
        n[-1, i] += 0.000001
        Ep = xc.calculate_spherical(rgd, n, v)
        n[-1, i] -= 0.000002
        Em = xc.calculate_spherical(rgd, n, v)
        x2 = (Ep - Em) / 0.000002
        print(name, nspins, E, x, x2, x - x2)
        equal(x, x2, 1e-9)
        n[-1, i] += 0.000001
        if nspins == 1:
            ns = rgd.empty(2)
            ns[:] = n / 2
            Es = xc.calculate_spherical(rgd, ns, 0 * ns)
            equal(E, Es, 1e-13)

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

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

        self.log('Basis functions: %s (%s)' %
                 (', '.join([str(nb) for nb in nb_l]),
                  ', '.join('spdf'[:lmax + 1])))

        self.vr_sg = self.rgd.zeros(self.nspins)
        self.vr_sg[:] = -self.Z

    def solve(self):
        """Diagonalize Schrödinger equation."""
        self.eeig = 0.0
        for channel in self.channels:
            if self.method == 'Gaussian basis-set':
                channel.solve(self.vr_sg[channel.s])
            else:
                channel.solve2(self.vr_sg[channel.s], self.scalar_relativistic)
            self.eeig += channel.get_eigenvalue_sum()

    def calculate_density(self):
        """Calculate elctron density and kinetic energy."""
        self.n_sg = self.rgd.zeros(self.nspins)
        for channel in self.channels:
            self.n_sg[channel.s] += channel.calculate_density()

    def calculate_electrostatic_potential(self):
        """Calculate electrostatic potential and energy."""
        n_g = self.n_sg.sum(0)
        self.vHr_g = self.rgd.poisson(n_g)        
        self.eH = 0.5 * self.rgd.integrate(n_g * self.vHr_g, -1)
        self.eZ = -self.Z * self.rgd.integrate(n_g, -1)
        
    def calculate_xc_potential(self):
        self.vxc_sg = self.rgd.zeros(self.nspins)
        self.exc = self.xc.calculate_spherical(self.rgd, self.n_sg, self.vxc_sg)

    def step(self):
        self.solve()
        self.calculate_density()
        self.calculate_electrostatic_potential()
        self.calculate_xc_potential()
        self.vr_sg = self.vxc_sg * self.rgd.r_g
        self.vr_sg += self.vHr_g
        self.vr_sg -= self.Z
        self.ekin = (self.eeig -
                     self.rgd.integrate((self.vr_sg * self.n_sg).sum(0), -1))
        
    def run(self, mix=0.4, maxiter=117, dnmax=1e-9):
        if self.channels is None:
            self.initialize()

        if self.dirac:
            equation = 'Dirac'
        elif self.scalar_relativistic:
            equation = 'scalar-relativistic Schrödinger'
        else:
            equation = 'non-relativistic Schrödinger'
        self.log('\nSolving %s equation using %s:' % (equation, self.method))

        dn = self.Z
        
        for iter in range(maxiter):
            self.log('.', end='')
            self.fd.flush()
            if iter > 0:
                self.vr_sg *= mix
                self.vr_sg += (1 - mix) * vr_old_sg

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

#.........这里部分代码省略.........
            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:
            # calculate hartree potential
            hartree(0, n * r * dr, r, vHr)

            # add potential from nuclear point charge (v = -Z / r)
            vHr -= Z

            # calculated exchange correlation potential and energy
            self.vXC[:] = 0.0

            if self.xc.type == 'GLLB':
                # Update the potential to self.vXC an the energy to self.Exc
                Exc = self.xc.get_xc_potential_and_energy_1d(self.vXC)
            else:
                Exc = self.xc.calculate_spherical(self.rgd,
                                                  n.reshape((1, -1)),
                                                  self.vXC.reshape((1, -1)))

            # calculate new total Kohn-Sham effective potential and
            # admix with old version
            vr[:] = (vHr + self.vXC * r) / self.tf_coeff

            if niter > 0:
                vr[:] = mix * vr + (1 - mix) * vrold
            vrold = vr.copy()

            # solve Kohn-Sham equation and determine the density change
            self.solve()
            dn = self.calculate_density() - n
            n += dn

            # estimate error from the square of the density change integrated
            q = log(np.sum((r * dn)**2))

            # print progress bar
            if niter == 0:
                q0 = q
                b0 = 0
            else:
                b = int((q0 - q) / (q0 - qOK) * 50)
                if b > b0:
                    self.txt.write(bar[b0:min(b, 50)])
                    self.txt.flush()
                    b0 = b

            # check if converged and break loop if so
            if q < qOK:

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

#.........这里部分代码省略.........
                    j = abs(k) - 0.5
                    f_n = (2 * j + 1) / (4 * l + 2) * np.array(f_n)
                    self.channels.append(DiracChannel(k, f_n, basis))

        self.log('Basis functions: %s (%s)' %
                 (', '.join([str(nb) for nb in nb_l]),
                  ', '.join('spdf'[:lmax + 1])))

        self.vr_sg = self.gd.zeros(self.nspins)
        self.vr_sg[:] = -self.Z

    def solve(self):
        """Diagonalize Schrödinger equation."""
        self.eeig = 0.0
        for channel in self.channels:
            channel.solve(self.vr_sg[channel.s])
            self.eeig += channel.get_eigenvalue_sum()

    def calculate_density(self):
        """Calculate elctron density and kinetic energy."""
        self.n_sg = self.gd.zeros(self.nspins)
        for channel in self.channels:
            self.n_sg[channel.s] += channel.calculate_density()

    def calculate_electrostatic_potential(self):
        """Calculate electrostatic potential and energy."""
        n_g = self.n_sg.sum(0)
        self.vHr_g = self.gd.poisson(n_g)        
        self.eH = 0.5 * self.gd.integrate(n_g * self.vHr_g, -1)
        self.eZ = -self.Z * self.gd.integrate(n_g, -1)
        
    def calculate_xc_potential(self):
        self.vxc_sg = self.gd.zeros(self.nspins)
        self.exc = self.xc.calculate_spherical(self.gd, self.n_sg, self.vxc_sg)

    def step(self):
        self.solve()
        self.calculate_density()
        self.calculate_electrostatic_potential()
        self.calculate_xc_potential()
        self.vr_sg = self.vxc_sg * self.gd.r_g
        self.vr_sg += self.vHr_g
        self.vr_sg -= self.Z
        self.ekin = (self.eeig -
                     self.gd.integrate((self.vr_sg * self.n_sg).sum(0), -1))
        
    def run(self, mix=0.4, maxiter=117, dnmax=1e-9):
        if self.channels is None:
            self.initialize()

        dn = self.Z
        pb = ProgressBar(log(dnmax / dn), 0, 53, self.fd)
        self.log()
        
        for iter in range(maxiter):
            if iter > 1:
                self.vr_sg *= mix
                self.vr_sg += (1 - mix) * vr_old_sg
                dn = self.gd.integrate(abs(self.n_sg - n_old_sg).sum(0))
                pb(log(dnmax / dn))
                if dn <= dnmax:
                    break

            vr_old_sg = self.vr_sg
            n_old_sg = self.n_sg
            self.step()

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

#.........这里部分代码省略.........
            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)
            a1z_g = np.dot(A_Li[:, 2], nt_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, dntdr_Lg)
            a2_g += a1_g**2
            v_g = np.zeros(c.ng)
            e_g = np.zeros(c.ng)
            deda2_g = np.zeros(c.ng)

            v_g[:] = 0.0
            e_g[:] = 0.0
            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.dv_g)
            x_g = -2.0 * deda2_g * c.dv_g * a1_g
            c.rgd.derivative2(x_g, x_g)
            x_g += v_g * c.dv_g

            B_Lqp = c.B_pqL.T
            dEdD_p -= w * np.dot(dot3(c.B_pqL, Y_L),
                                  np.dot(c.nt_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.nt_qg, x_g * a1x_g))
            dEdD_p -= w * np.dot(dot3(c.B_pqL,
                                       A_Li[:, 1]),
                                  np.dot(c.nt_qg, x_g * a1y_g))
            
            dEdD_p -= w * np.dot(dot3(c.B_pqL,
                                       A_Li[:, 2]),
                                  np.dot(c.nt_qg, x_g * a1z_g))
            
            y += 1
        
        return (E) * self.weight

    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_spherical [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,
#.........这里部分代码省略.........