Python GridDescriptor.zeros方法代码示例

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import zeros [as 别名]
 def __init__(self, gd, spline_j, spos_c, index=None):
     rcut = max([spline.get_cutoff() for spline in spline_j])
     corner_c = np.ceil(spos_c * gd.N_c - rcut / gd.h_c).astype(int)
     size_c = np.ceil(spos_c * gd.N_c + rcut / gd.h_c).astype(int) - corner_c
     smallgd = GridDescriptor(N_c=size_c + 1, cell_cv=gd.h_c * (size_c + 1), pbc_c=False, comm=mpi.serial_comm)
     lfc = LFC(smallgd, [spline_j])
     lfc.set_positions((spos_c[np.newaxis, :] * gd.N_c - corner_c + 1) / smallgd.N_c)
     ni = lfc.Mmax
     f_iG = smallgd.zeros(ni)
     lfc.add(f_iG, {0: np.eye(ni)})
     LocalizedFunctions.__init__(self, gd, f_iG, corner_c, index=index)

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import zeros [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 zeros [as 别名]
def interpolate_2d(mat):
    from gpaw.grid_descriptor import GridDescriptor
    from gpaw.transformers import Transformer
    nn = 10
    N_c = np.zeros([3], dtype=int)
    N_c[1:] = mat.shape[:2]
    N_c[0] = nn
    bmat = np.resize(mat, N_c)
    gd = GridDescriptor(N_c, N_c)
    finegd = GridDescriptor(N_c * 2, N_c)
    interpolator = Transformer(gd, finegd, 3)
    fine_bmat = finegd.zeros()
    interpolator.apply(bmat, fine_bmat)
    return fine_bmat[0]

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import zeros [as 别名]
def make_dummy_kpt_reference(l, function, k_c,
                             rcut=6., a=10., n=60, dtype=complex):
    r = np.linspace(0., rcut, 300)
    mcount = 2*l + 1
    fcount = 1
    kcount = 1
    gd = GridDescriptor((n, n, n), (a, a, a), (True, True, True))
    kpt = KPoint([], gd, 1., 0, 0, 0, k_c, dtype)
    spline = Spline(l, r[-1], function(r))
    center = (.5, .5, .5)
    lf = create_localized_functions([spline], gd, center, dtype=dtype)
    lf.set_phase_factors([kpt.k_c])
    psit_nG = gd.zeros(mcount, dtype=dtype)
    coef_xi = np.identity(mcount * fcount, dtype=dtype)
    lf.add(psit_nG, coef_xi, k=0)
    kpt.psit_nG = psit_nG
    print 'Number of boxes', len(lf.box_b)
    print 'Phase kb factors shape', lf.phase_kb.shape
    return gd, kpt, center

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import zeros [as 别名]
def test():
    from gpaw.grid_descriptor import GridDescriptor

    ngpts = 40
    h = 1.0 / ngpts
    N_c = (ngpts, ngpts, ngpts)
    a = h * ngpts
    gd = GridDescriptor(N_c, (a, a, a))
    
    from gpaw.spline import Spline
    a = np.array([1, 0.9, 0.8, 0.0])
    s = Spline(0, 0.2, a)
    x = LocalizedFunctionsCollection(gd, [[s], [s]])
    x.set_positions([(0.5, 0.45, 0.5), (0.5, 0.55, 0.5)])
    n_G = gd.zeros()
    x.add(n_G)
    import pylab as plt
    plt.contourf(n_G[20, :, :])
    plt.axis('equal')
    plt.show()

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import zeros [as 别名]
def make_dummy_reference(l, function=None, rcut=6., a=12., n=60,
                         dtype=float):
    """Make a mock reference wave function using a made-up radial function
    as reference"""
    #print 'Dummy reference: l=%d, rcut=%.02f, alpha=%.02f' % (l, rcut, alpha)
    r = np.arange(0., rcut, .01)

    if function is None:
        function = QuasiGaussian(4., rcut)

    norm = get_norm(r, function(r), l)
    function.renormalize(norm)
    #g = QuasiGaussian(alpha, rcut)
    
    mcount = 2*l + 1
    fcount = 1
    gd = GridDescriptor((n, n, n), (a, a, a), (False, False, False))
    spline = Spline(l, r[-1], function(r), points=50)
    center = (.5, .5, .5)
    lf = create_localized_functions([spline], gd, center, dtype=dtype)
    psit_k = gd.zeros(mcount, dtype=dtype)
    coef_xi = np.identity(mcount * fcount, dtype=dtype)
    lf.add(psit_k, coef_xi)
    return gd, psit_k, center, function

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import zeros [as 别名]
from gpaw.fd_operators import Gradient
import numpy as np
from gpaw.grid_descriptor import GridDescriptor
from gpaw.mpi import world

if world.size > 4:
    # Grid is so small that domain decomposition cannot exceed 4 domains
    assert world.size % 4 == 0
    group, other = divmod(world.rank, 4)
    ranks = np.arange(4*group, 4*(group+1))
    domain_comm = world.new_communicator(ranks)
else:
    domain_comm = world

gd = GridDescriptor((8, 1, 1), (8.0, 1.0, 1.0), comm=domain_comm)
a = gd.zeros()
dadx = gd.zeros()
a[:, 0, 0] = np.arange(gd.beg_c[0], gd.end_c[0])
gradx = Gradient(gd, v=0)
print a.itemsize, a.dtype, a.shape
print dadx.itemsize, dadx.dtype, dadx.shape
gradx.apply(a, dadx)

#   a = [ 0.  1.  2.  3.  4.  5.  6.  7.]
#
#   da
#   -- = [-2.5  1.   1.   1.   1.   1.  1.  -2.5]
#   dx

dadx = gd.collect(dadx, broadcast=True)
assert dadx[3, 0, 0] == 1.0 and np.sum(dadx[:, 0, 0]) == 0.0

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import zeros [as 别名]
from time import time
from gpaw.grid_descriptor import GridDescriptor
from gpaw.fd_operators import Laplace
import gpaw.mpi as mpi

n = 96
h = 0.1
L = n * h
gd = GridDescriptor((n, n, n), [L, L, L])

# Allocate arrays:
a = gd.zeros(100) + 1.2
b = gd.empty(100)

o = Laplace(gd, 2).apply

t0 = time()
for r in range(10):
    o(a, b)

if mpi.rank == 0:
    print time() - t0, a.shape
    

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import zeros [as 别名]
from gpaw.grid_descriptor import GridDescriptor
from gpaw.spline import Spline
a = 4.0
gd = GridDescriptor(N_c=[16, 20, 20], cell_cv=[a, a + 1, a + 2],
                    pbc_c=(0, 1, 1))
spos_ac = np.array([[0.25, 0.15, 0.35], [0.5, 0.5, 0.5]])
kpts_kc = None
s = Spline(l=0, rmax=2.0, f_g=np.array([1, 0.9, 0.1, 0.0]))
p = Spline(l=1, rmax=2.0, f_g=np.array([1, 0.9, 0.1, 0.0]))
spline_aj = [[s], [s, p]]
c = LFC(gd, spline_aj, cut=True, forces=True)
c.set_positions(spos_ac)
C_ani = c.dict(3, zero=True)
if 1 in C_ani:
    C_ani[1][:, 1:] = np.eye(3)
psi = gd.zeros(3)
c.add(psi, C_ani)
c.integrate(psi, C_ani)
if 1 in C_ani:
    d = C_ani[1][:, 1:].diagonal()
    assert d.ptp() < 4e-6
    C_ani[1][:, 1:] -= np.diag(d)
    assert abs(C_ani[1]).max() < 5e-17
d_aniv = c.dict(3, derivative=True)
c.derivative(psi, d_aniv)
if 1 in d_aniv:
    for v in range(3):
        assert abs(d_aniv[1][v - 1, 0, v] + 0.2144) < 5e-5
        d_aniv[1][v - 1, 0, v] = 0
    assert abs(d_aniv[1]).max() < 3e-16
eps = 0.0001

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import zeros [as 别名]
class DF(CHI):
    """This class defines dielectric function related physical quantities."""

    def __init__(self,
                 calc=None,
                 nbands=None,
                 w=None,
                 q=None,
                 eshift=None,
                 ecut=10.,
                 density_cut=None,
                 G_plus_q=False,
                 eta=0.2,
                 rpad=None,
                 vcut=None,
                 ftol=1e-7,
                 txt=None,
                 xc='ALDA',
                 print_xc_scf=False,
                 hilbert_trans=True,
                 time_ordered=False,
                 optical_limit=False,
                 comm=None,
                 kcommsize=None):

        CHI.__init__(self, calc=calc, nbands=nbands, w=w, q=q, eshift=eshift,
                     ecut=ecut, density_cut=density_cut,
                     G_plus_q=G_plus_q, eta=eta, rpad=rpad, vcut=vcut,
                     ftol=ftol, txt=txt, xc=xc, hilbert_trans=hilbert_trans,
                     time_ordered=time_ordered, optical_limit=optical_limit,
                     comm=comm, kcommsize=kcommsize)

        self.df_flag = False
        self.print_bootstrap = print_xc_scf
        self.df1_w = None  # NLF RPA
        self.df2_w = None  # LF RPA
        self.df3_w = None  # NLF ALDA
        self.df4_w = None  # LF ALDA

    def get_dielectric_matrix(self,
                              xc='RPA',
                              overwritechi0=False,
                              symmetric=True,
                              chi0_wGG=None,
                              calc=None,
                              vcut=None,
                              dir=None):
	if self.chi0_wGG is None and chi0_wGG is None:
            self.initialize()
            self.calculate()
        elif self.chi0_wGG is None and chi0_wGG is not None:
            #Read from file and reinitialize 
            self.xc = xc

            from gpaw.response.parallel import par_read 
            self.chi0_wGG = par_read(chi0_wGG, 'chi0_wGG')
            self.nvalbands = self.nbands
            #self.parallel_init() # parallelization not yet implemented
            self.Nw_local = self.Nw  # parallelization not yet implemented
            if self.calc is None:
                from gpaw import GPAW
                self.calc = GPAW(calc,txt=None)
            if self.xc == 'ALDA' or self.xc == 'ALDA_X':
                from gpaw.response.kernel import calculate_Kxc
                from gpaw.grid_descriptor import GridDescriptor
                from gpaw.mpi import world, rank, size, serial_comm
                    
                self.pbc = self.calc.atoms.pbc
                self.gd = GridDescriptor(self.calc.wfs.gd.N_c*self.rpad, self.acell_cv,
                                         pbc_c=True, comm=serial_comm)

                R_av = self.calc.atoms.positions / Bohr
                nt_sg = self.calc.density.nt_sG
                    
                if (self.rpad > 1).any() or (self.pbc - True).any():
                    nt_sG = self.gd.zeros(self.nspins)
                    #nt_sG = np.zeros([self.nspins, self.nG[0], self.nG[1], self.nG[2]])
                    for s in range(self.nspins):
                        nt_G = self.pad(nt_sg[s])
                        nt_sG[s] = nt_G
                else:
                    nt_sG = nt_sg
                        
                self.Kxc_sGG = calculate_Kxc(self.gd, 
                                             nt_sG,
                                             self.npw, self.Gvec_Gc,
                                             self.gd.N_c, self.vol,
                                             self.bcell_cv, R_av,
                                             self.calc.wfs.setups,
                                             self.calc.density.D_asp,
                                             functional=self.xc,
                                             density_cut=self.density_cut)

        if overwritechi0:
            dm_wGG = self.chi0_wGG
        else:
            dm_wGG = np.zeros_like(self.chi0_wGG)

        if dir is None:
            q_c = self.q_c
#.........这里部分代码省略.........

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import zeros [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 zeros [as 别名]
from gpaw.response.math_func import two_phi_planewave_integrals

# Initialize s, p, d (9 in total) wave and put them on grid
rc = 2.0
a = 2.5 * rc
n = 64
lmax = 2
b = 8.0
m = (lmax + 1)**2
gd = GridDescriptor([n, n, n], [a, a, a])
r = np.linspace(0, rc, 200)
g = np.exp(-(r / rc * b)**2)
splines = [Spline(l=l, rmax=rc, f_g=g) for l in range(lmax + 1)]
c = LFC(gd, [splines])
c.set_positions([(0.5, 0.5, 0.5)])
psi = gd.zeros(m)
d0 = c.dict(m)
if 0 in d0:
    d0[0] = np.identity(m)
c.add(psi, d0)

# Calculate on 3d-grid < phi_i | e**(-ik.r) | phi_j >
R_a = np.array([a/2,a/2,a/2])
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)

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

    for i in range(4):
        if here:

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