Python GridDescriptor.integrate方法代码示例

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import integrate [as 别名]
def f(n, p):
    N = 2 * n
    gd = GridDescriptor((N, N, N), (L, L, L))
    a = gd.zeros()
    print(a.shape)
    #p = PoissonSolver(nn=1, relax=relax)
    p.set_grid_descriptor(gd)
    p.initialize()
    cut = N / 2.0 * 0.9
    s = Spline(l=0, rmax=cut, f_g=np.array([1, 0.5, 0.0]))
    c = LFC(gd, [[s], [s]])
    c.set_positions([(0, 0, 0), (0.5, 0.5, 0.5)])
    c.add(a)

    I0 = gd.integrate(a)
    a -= gd.integrate(a) / L**3

    I = gd.integrate(a)
    b = gd.zeros()
    p.solve(b, a, charge=0)#, eps=1e-20)
    return gd.collect(b, broadcast=1)

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import integrate [as 别名]
nbands = int(sys.argv[2])
repeats = 10

gd = GridDescriptor([gpts, gpts, gpts], comm=world, parsize=size)

a = gd.empty(nbands)
b = gd.empty(nbands)

np.random.seed(10)
a[:] = np.random.random(a.shape)
b[:] = np.random.random(b.shape)
c = np.zeros((nbands, nbands))

# warm-up
for i in range(3):
    gd.integrate(a, b, hermitian=False, _transposed_result=c)
    gemm(1.0, a, c, 0.0, b)
# equal(np.sum(c), 3600.89536641, 1e-6)
t0 = time()
for i in range(repeats):
    gd.integrate(a, b, hermitian=False, _transposed_result=c)
    gemm(1.0, a, c, 0.0, b)
t1 = time()
if rank == 0:
    print "Check", np.sum(b), "Time", (t1 - t0) / repeats

a_mic = gd.empty(nbands, usemic=True)
b_mic = gd.empty(nbands, usemic=True)
c_mic = stream.bind(c)
np.random.seed(10)
a_mic.array[:] = np.random.random(a_mic.shape)

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import integrate [as 别名]
v = gauss.remove_moment(nH, 0)
m = gauss.get_moment(nH, 0)
print("\nZero'th moment of compensated Hydrogen density =", m)
equal(m, 0.0, 1e-7)


# /-------------------------------------------------\
# | Check if Gaussian potentials are made correctly |
# \-------------------------------------------------/

# Array for storing the potential
pot = gd.zeros(dtype=float, global_array=False)
for L in range(7):  # Angular index of gaussian
    # Get analytic functions
    ng = gauss.get_gauss(L)
    vg = gauss.get_gauss_pot(L)

    # Solve potential numerically
    niter = solve(pot, ng, charge=None, zero_initial_phi=True)

    # Determine residual
    residual = norm(pot - vg)
    residual = gd.integrate((pot - vg) ** 2) ** 0.5

    # print result
    print("L=%s, processor %s of %s: %s" % (L, gd.comm.rank + 1, gd.comm.size, residual))

    assert residual < 0.6

# mpirun -np 2 python gauss_func.py --gpaw-parallel --gpaw-debug

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import integrate [as 别名]
from gpaw.test import equal
from gpaw.grid_descriptor import GridDescriptor
from gpaw.spline import Spline
import gpaw.mpi as mpi
from gpaw.lfc import LocalizedFunctionsCollection as LFC

s = Spline(0, 1.0, [1.0, 0.5, 0.0])
n = 40
a = 8.0
gd = GridDescriptor((n, n, n), (a, a, a), comm=mpi.serial_comm)
c = LFC(gd, [[s], [s], [s]])
c.set_positions([(0.5, 0.5, 0.25 + 0.25 * i) for i in [0, 1, 2]])
b = gd.zeros()
c.add(b)
x = gd.integrate(b)

gd = GridDescriptor((n, n, n), (a, a, a), comm=mpi.serial_comm)
c = LFC(gd, [[s], [s], [s]])
c.set_positions([(0.5, 0.5, 0.25 + 0.25 * i) for i in [0, 1, 2]])
b = gd.zeros()
c.add(b)
y = gd.integrate(b)
equal(x, y, 1e-13)

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import integrate [as 别名]
from gpaw.transformers import Transformer
import numpy.random as ra
from gpaw.grid_descriptor import GridDescriptor

p = 0
n = 20
gd1 = GridDescriptor((n, n, n), (8.0, 8.0, 8.0), pbc_c=p)
a1 = gd1.zeros()
ra.seed(8)
a1[:] = ra.random(a1.shape)
gd2 = gd1.refine()
a2 = gd2.zeros()
i = Transformer(gd1, gd2).apply
i(a1, a2)
assert abs(gd1.integrate(a1) - gd2.integrate(a2)) < 1e-10
r = Transformer(gd2, gd1).apply
a2[0] = 0.0
a2[:, 0] = 0.0
a2[:, :, 0] = 0.0
a2[-1] = 0.0
a2[:, -1] = 0.0
a2[:, :, -1] = 0.0
r(a2, a1)
assert abs(gd1.integrate(a1) - gd2.integrate(a2)) < 1e-10

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import integrate [as 别名]
class UTDomainParallelSetup(TestCase):
    """
    Setup a simple domain parallel calculation."""

    # Number of bands
    nbands = 1

    # Spin-paired
    nspins = 1

    # Mean spacing and number of grid points per axis
    h = 0.2 / Bohr

    # Generic lattice constant for unit cell
    a = 5.0 / Bohr

    # Type of boundary conditions employed
    boundaries = None

    # Type of unit cell employed
    celltype = None

    def setUp(self):
        for virtvar in ['boundaries', 'celltype']:
            assert getattr(self,virtvar) is not None, 'Virtual "%s"!' % virtvar

        # Basic unit cell information:
        pbc_c = {'zero'    : (False,False,False), \
                 'periodic': (True,True,True), \
                 'mixed'   : (True, False, True)}[self.boundaries]
        a, b = self.a, 2**0.5*self.a
        cell_cv = {'general'   : np.array([[0,a,a],[a/2,0,a/2],[a/2,a/2,0]]),
                   'rotated'   : np.array([[0,0,b],[b/2,0,0],[0,b/2,0]]),
                   'inverted'   : np.array([[0,0,b],[0,b/2,0],[b/2,0,0]]),
                   'orthogonal': np.diag([b, b/2, b/2])}[self.celltype]
        cell_cv = np.array([(4-3*pbc)*c_v for pbc,c_v in zip(pbc_c, cell_cv)])

        # Decide how many kpoints to sample from the 1st Brillouin Zone
        kpts_c = np.ceil((10/Bohr)/np.sum(cell_cv**2,axis=1)**0.5).astype(int)
        kpts_c = tuple(kpts_c*pbc_c + 1-pbc_c)
        bzk_kc = kpts2ndarray(kpts_c)
        self.gamma = len(bzk_kc) == 1 and not bzk_kc[0].any()

        #p = InputParameters()
        #Z_a = self.atoms.get_atomic_numbers()
        #xcfunc = XC(p.xc)
        #setups = Setups(Z_a, p.setups, p.basis, p.lmax, xcfunc)
        #symmetry, weight_k, self.ibzk_kc = reduce_kpoints(self.atoms, bzk_kc,
        #                                                  setups, p.usesymm)

        self.ibzk_kc = bzk_kc.copy() # don't use symmetry reduction of kpoints
        self.nibzkpts = len(self.ibzk_kc)
        self.ibzk_kv = kpoint_convert(cell_cv, skpts_kc=self.ibzk_kc)

        # Parse parallelization parameters and create suitable communicators.
        #parsize_domain, parsize_bands = create_parsize_minbands(self.nbands, world.size)
        parsize_domain, parsize_bands = world.size//gcd(world.size, self.nibzkpts), 1
        assert self.nbands % np.prod(parsize_bands) == 0
        domain_comm, kpt_comm, band_comm = distribute_cpus(parsize_domain,
            parsize_bands, self.nspins, self.nibzkpts)

        # Set up band descriptor:
        self.bd = BandDescriptor(self.nbands, band_comm)

        # Set up grid descriptor:
        N_c = np.round(np.sum(cell_cv**2, axis=1)**0.5 / self.h)
        N_c += 4-N_c % 4 # makes domain decomposition easier
        self.gd = GridDescriptor(N_c, cell_cv, pbc_c, domain_comm, parsize_domain)
        self.assertEqual(self.gamma, np.all(~self.gd.pbc_c))

        # What to do about kpoints?
        self.kpt_comm = kpt_comm

        if debug and world.rank == 0:
            comm_sizes = tuple([comm.size for comm in [world, self.bd.comm, \
                                                   self.gd.comm, self.kpt_comm]])
            print '%d world, %d band, %d domain, %d kpt' % comm_sizes

    def tearDown(self):
        del self.ibzk_kc, self.ibzk_kv, self.bd, self.gd, self.kpt_comm

    # =================================

    def verify_comm_sizes(self):
        if world.size == 1:
            return
        comm_sizes = tuple([comm.size for comm in [world, self.bd.comm, \
                                                   self.gd.comm, self.kpt_comm]])
        self._parinfo =  '%d world, %d band, %d domain, %d kpt' % comm_sizes
        self.assertEqual(self.nbands % self.bd.comm.size, 0)
        self.assertEqual((self.nspins*self.nibzkpts) % self.kpt_comm.size, 0)

    def verify_grid_volume(self):
        gdvol = np.prod(self.gd.get_size_of_global_array())*self.gd.dv
        self.assertAlmostEqual(self.gd.integrate(1+self.gd.zeros()), gdvol, 10)

    def verify_grid_point(self):
        # Volume integral of cartesian coordinates of all available grid points
        gdvol = np.prod(self.gd.get_size_of_global_array())*self.gd.dv
        cmr_v = self.gd.integrate(self.gd.get_grid_point_coordinates()) / gdvol
#.........这里部分代码省略.........

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import integrate [as 别名]
        [(2, 4, 4), (4, 8, 8)],
        [(2, 4, 2), (4, 8, 4)]
        ]:
        print(size1, size2)
        gd1 = GridDescriptor(size1, size1)
        gd2 = GridDescriptor(size2, size1)
        pd1 = PWDescriptor(1, gd1, complex)
        pd2 = PWDescriptor(1, gd2, complex)
        pd1r = PWDescriptor(1, gd1)
        pd2r = PWDescriptor(1, gd2)
        for R1, R2 in [[(0,0,0), (0,0,0)],
                       [(0,0,0), (0,0,1)]]:
            x = test(gd1, gd2, pd1, pd2, R1, R2)
            y = test(gd1, gd2, pd1r, pd2r ,R1, R2)
            equal(x, y, 1e-9)
        
        a1 = np.random.random(size1)
        a2 = pd1r.interpolate(a1, pd2r)[0]
        c2 = pd1.interpolate(a1 + 0.0j, pd2)[0]
        d2 = pd1.interpolate(a1 * 1.0j, pd2)[0]
        equal(gd1.integrate(a1), gd2.integrate(a2), 1e-13)
        equal(abs(c2 - a2).max(), 0, 1e-14)
        equal(abs(d2 - a2 * 1.0j).max(), 0, 1e-14)

        a1 = pd2r.restrict(a2, pd1r)[0]
        c1 = pd2.restrict(a2 + 0.0j, pd1)[0]
        d1 = pd2.restrict(a2 * 1.0j, pd1)[0]
        equal(gd1.integrate(a1), gd2.integrate(a2), 1e-13)
        equal(abs(c1 - a1).max(), 0, 1e-14)
        equal(abs(d1 - a1 * 1.0j).max(), 0, 1e-14)

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import integrate [as 别名]
from gpaw.grid_descriptor import GridDescriptor
gd = GridDescriptor([4, 4, 4])
a = gd.empty(dtype=complex)
a[:] = 1.0
assert gd.integrate(a.real, a.real) == 1.0

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import integrate [as 别名]
from gpaw.grid_descriptor import GridDescriptor
from gpaw.transformers import Transformer


n = 20
gd = GridDescriptor((n,n,n))
np.random.seed(8)
a = gd.empty()
a[:] = np.random.random(a.shape)

gd2 = gd.refine()
b = gd2.zeros()
for k in [2, 4, 6, 8]:
    inter = Transformer(gd, gd2, k // 2).apply
    inter(a, b)
    assert abs(gd.integrate(a) - gd2.integrate(b)) < 1e-14

gd2 = gd.coarsen()
b = gd2.zeros()
for k in [2, 4, 6, 8]:
    restr = Transformer(gd, gd2, k // 2).apply
    restr(a, b)
    assert abs(gd.integrate(a) - gd2.integrate(b)) < 1e-14

# complex versions
a = gd.empty(dtype=complex)
a.real = np.random.random(a.shape)
a.imag = np.random.random(a.shape)

phase = np.ones((3, 2), complex)

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import integrate [as 别名]
gd = GridDescriptor((n, n, n), (a, a, a))
c_LL = np.identity(9, float)
a_Lg = gd.zeros(9)
nspins = 2
xc = XC('LDA')
for soft in [False]:
    s = create_setup('Cu', xc, lmax=2)
    ghat_l = s.ghat_l
    ghat_Lg = create_localized_functions(ghat_l, gd, (0.54321, 0.5432, 0.543))
    a_Lg[:] = 0.0
    ghat_Lg.add(a_Lg, c_LL)
    for l in range(3):
        for m in range(2 * l + 1):
            L = l**2 + m
            a_g = a_Lg[L]
            Q0 = gd.integrate(a_g) / sqrt(4 * pi)
            Q1_m = -gd.calculate_dipole_moment(a_g) / sqrt(4 * pi / 3)
            print Q0
            if l == 0:
                Q0 -= 1.0
                Q1_m[:] = 0.0
            elif l == 1:
                Q1_m[(m + 1) % 3] -= 1.0
            print soft, l, m, Q0, Q1_m
            assert abs(Q0) < 2e-6
            assert np.alltrue(abs(Q1_m) < 3e-5)
    b_Lg = np.reshape(a_Lg, (9, -1))
    S_LL = np.inner(b_Lg, b_Lg)
    gd.comm.sum(S_LL)
    S_LL.ravel()[::10] = 0.0
    print max(abs(S_LL).ravel())

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import integrate [as 别名]
rr = gd.get_grid_point_coordinates()
for dim in range(3):
    rr[dim] -= R_a[dim]

k_G = np.array([[11.,0.2,0.1],[10., 0., 10.]])
nkpt = k_G.shape[0]

d0 = np.zeros((nkpt,m,m), dtype=complex)
for i in range(m):
    for j in range(m):
        for ik in range(nkpt):
            k = k_G[ik]
            kk = np.sqrt(np.inner(k,k))
            kr = np.inner(k,rr.T).T
            expkr = np.exp(-1j * kr)
            d0[ik, i,j] = gd.integrate(psi[i] * psi[j] * expkr)

# Calculate on 1d-grid < phi_i | e**(-ik.r) | phi_j >
rgd = RadialGridDescriptor(r, np.ones_like(r) * r[1])
g = [np.exp(-(r / rc * b)**2)*r**l for l in range(lmax + 1)]
l_j = range(lmax + 1)
d1 = two_phi_planewave_integrals(k_G, rgd=rgd, phi_jg=g,
                            phit_jg=np.zeros_like(g), l_j=l_j)

d1 = d1.reshape(nkpt, m, m)

for i in range(m):
    for j in range(m):
        for ik in range(nkpt):
            if np.abs(d0[ik,i,j] - d1[ik,i,j]) > 1e-10:
                print i, j, d0[ik,i,j]- d1[ik,i,j]