Python GridDescriptor.empty方法代码示例

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import empty [as 别名]
def go(comm, ngpts, repeat, narrays, out, prec):
    N_c = np.array((ngpts, ngpts, ngpts))
    a = 10.0
    gd = GridDescriptor(N_c, (a, a, a), comm=comm))
    gdcoarse = gd.coarsen()
    gdfine = gd.refine()
    kin1 = Laplace(gd, -0.5, 1).apply
    laplace = Laplace(gd, -0.5, 2)
    kin2 = laplace.apply
    restrict = Transformer(gd, gdcoarse, 1).apply
    interpolate = Transformer(gd, gdfine, 1).apply
    precondition = Preconditioner(gd, laplace, np_float)
    a1 = gd.empty(narrays)
    a1[:] = 1.0
    a2 = gd.empty(narrays)
    c = gdcoarse.empty(narrays)
    f = gdfine.empty(narrays)

    T = [0, 0, 0, 0, 0]
    for i in range(repeat):
        comm.barrier()
        kin1(a1, a2)
        comm.barrier()
        t0a = time()
        kin1(a1, a2)
        t0b = time()
        comm.barrier()
        t1a = time()
        kin2(a1, a2)
        t1b = time()
        comm.barrier()
        t2a = time()
        for A, C in zip(a1,c):
            restrict(A, C)
        t2b = time()
        comm.barrier()
        t3a = time()
        for A, F in zip(a1,f):
            interpolate(A, F)
        t3b = time()
        comm.barrier()
        if prec:
            t4a = time()
            for A in a1:
                precondition(A, None, None, None)
            t4b = time()
            comm.barrier()

        T[0] += t0b - t0a
        T[1] += t1b - t1a
        T[2] += t2b - t2a
        T[3] += t3b - t3a
        if prec:
            T[4] += t4b - t4a

    if mpi.rank == 0:
        out.write(' %2d %2d %2d' % tuple(gd.parsize_c))
        out.write(' %12.6f %12.6f %12.6f %12.6f %12.6f\n' %
                  tuple([t / repeat / narrays for t in T]))
        out.flush()

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import empty [as 别名]
class WignerSeitzTruncatedCoulomb:
    def __init__(self, cell_cv, nk_c, txt=sys.stdout):
        self.nk_c = nk_c
        bigcell_cv = cell_cv * nk_c[:, np.newaxis]
        L_c = (np.linalg.inv(bigcell_cv)**2).sum(0)**-0.5
        
        rc = 0.5 * L_c.min()
        prnt('Inner radius for %dx%dx%d Wigner-Seitz cell: %.3f Ang' %
             (tuple(nk_c) + (rc * Bohr,)), file=txt)
        
        self.a = 5 / rc
        prnt('Range-separation parameter: %.3f Ang^-1' % (self.a / Bohr),
             file=txt)
        
#        nr_c = [get_efficient_fft_size(2 * int(L * self.a * 1.5))
        nr_c = [get_efficient_fft_size(2 * int(L * self.a * 3.0))
                for L in L_c]
        prnt('FFT size for calculating truncated Coulomb: %dx%dx%d' %
             tuple(nr_c), file=txt)
        
        self.gd = GridDescriptor(nr_c, bigcell_cv, comm=mpi.serial_comm)
        v_R = self.gd.empty()
        v_i = v_R.ravel()
        
        pos_iv = self.gd.get_grid_point_coordinates().reshape((3, -1)).T
        corner_jv = np.dot(np.indices((2, 2, 2)).reshape((3, 8)).T, bigcell_cv)
        for i, pos_v in enumerate(pos_iv):
            r = ((pos_v - corner_jv)**2).sum(axis=1).min()**0.5
            if r == 0:
                v_i[i] = 2 * self.a / pi**0.5
            else:
                v_i[i] = erf(self.a * r) / r
                
        self.K_Q = np.fft.fftn(v_R) * self.gd.dv
        
    def get_potential(self, pd):
        q_c = pd.kd.bzk_kc[0]
        shift_c = (q_c * self.nk_c).round().astype(int)
        max_c = self.gd.N_c // 2
        K_G = pd.zeros()
        N_c = pd.gd.N_c
        for G, Q in enumerate(pd.Q_qG[0]):
            Q_c = (np.unravel_index(Q, N_c) + N_c // 2) % N_c - N_c // 2
            Q_c = Q_c * self.nk_c + shift_c
            if (abs(Q_c) < max_c).all():
                K_G[G] = self.K_Q[tuple(Q_c)]

        G2_G = pd.G2_qG[0]
        a = self.a
        if pd.kd.gamma:
            K_G[0] += pi / a**2
        else:
            K_G[0] += 4 * pi * (1 - np.exp(-G2_G[0] / (4 * a**2))) / G2_G[0]
        K_G[1:] += 4 * pi * (1 - np.exp(-G2_G[1:] / (4 * a**2))) / G2_G[1:]
        assert pd.dtype == complex
        return K_G

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import empty [as 别名]
from gpaw.test import equal
import numpy as np
import time
import timeit
import os

import pyMIC as mic

offload_enabled = os.environ.get("GPAW_OFFLOAD");

device = mic.devices[0]

gpts = 64
gd = GridDescriptor([gpts, gpts, gpts])

a = gd.empty() # source
b = gd.empty() # result

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

op = Laplace(gd, 1.0, 3)

if offload_enabled:
    offl_a = device.bind(a)
    offl_b = device.bind(b)

print "--------------------------------------------------------------"
print "Starting relaxation step"
relax_start = time.time()

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

N = 20
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

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import empty [as 别名]
from gpaw.mpi import world, size, rank
from gpaw.utilities.blas import gemm
from gpaw.mic.micblas import gemm as mic_gemm
import pymic
import sys

device = pymic.devices[0]
stream = device.get_default_stream()

gpts = int(sys.argv[1])
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)

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import empty [as 别名]
                        parsize_bands=B,
                        nspins=1, nibzkpts=1)
domain_comm, kpt_comm, band_comm, block_comm = \
    [comms[name] for name in ['d', 'k', 'b', 'K']]

assert kpt_comm.size == 1
if world.rank == 0:
    print('MPI: %d domains, %d band groups' % (domain_comm.size, band_comm.size))

# Set up band and grid descriptors:
bd = BandDescriptor(N, band_comm, False)
gd = GridDescriptor((G, G, G), (a, a, a), True, domain_comm, parsize=D)
ksl = BandLayouts(gd, bd, block_comm, float)

# Random wave functions:
psit_mG = gd.empty(M)
for m in range(M):
    np.random.seed(world.rank * M + m)
    psit_mG[m] = np.random.uniform(-0.5, 0.5, tuple(gd.n_c))
if world.rank == 0:
    print('Size of wave function array:', psit_mG.shape)
P_ani = {0: psit_mG[:, :2, 0, 0].copy(),
         1: psit_mG[:, -1, -1, -3:].copy()}

kin = Laplace(gd, -0.5, 2).apply
vt_G = gd.empty()
vt_G.fill(0.567)

def run(psit_mG):
    overlap = MatrixOperator(ksl, J)
    def H(psit_xG):

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

ngpts = 80
N_c = (ngpts, ngpts, ngpts)
a = 10.0
gd = GridDescriptor(N_c, (a, a, a))
gdfine = gd.refine()
interpolate = Transformer(gd, gdfine, 3).apply
a1 = gd.empty()
a1[:] = 1.0
f = gdfine.empty()
ta = time()
r = 600
for i in range(r):
    interpolate(a1, f)
tb = time()
n = 8 * (1 + 2 + 4) * ngpts**3
print '%.3f GFlops' % (r * n / (tb - ta) * 1e-9)

"""
python: 0.330 GFlops
mpirun -np 2 gpaw-python: 0.500 GFlops
gpaw-python + OMP: 0.432 GFlops
"""

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import empty [as 别名]
from math import pi
from gpaw.grid_descriptor import GridDescriptor
from gpaw.xc import XC
from gpaw.xc.noncollinear import NonCollinearLDAKernel, NonCollinearFunctional
import numpy as np
from gpaw.test import equal

N = 20
a = 1.0
gd = GridDescriptor((N, N, N), (a, a, a))

for xc in [XC('LDA'), XC('PBE')]:
    n = gd.empty(2)
    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)
    n[1] += 0.02 + 0.01 * np.cos(2 * pi * z / a)
    n /= 2

    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

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import empty [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 empty [as 别名]
# Set up communicators:
r = world.rank // D * D
domain_comm = world.new_communicator(np.arange(r, r + D))
band_comm = world.new_communicator(np.arange(world.rank % D, world.size, D))

# Set up grid descriptor:
gd = GridDescriptor((G, G, G), (a, a, a), True, domain_comm, parsize)

# Random wave functions:
np.random.seed(world.rank)
psit_mG = np.random.uniform(-0.5, 0.5, size=(M,) + tuple(gd.n_c))
if world.rank == 0:
    print 'Size of wave function array:', psit_mG.shape

# Send and receive buffers:
send_mG = gd.empty(M)
recv_mG = gd.empty(M)

def run():
    S_nn = overlap(psit_mG, send_mG, recv_mG)

    t1 = time()
    if world.rank == 0:
        inverse_cholesky(S_nn)
        C_nn = S_nn
    else:
        C_nn = np.empty((N, N))
    t2 = time()

    if world.rank == 0:
        print 'Cholesky Time %f' % (t2-t1)

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import empty [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 empty [as 别名]
device = mic.devices[0]

gpts = 48
nbands = 256
if len(sys.argv) > 1:
    nbands = int(sys.argv[1])

gd = GridDescriptor([gpts, gpts, gpts])

alpha = 1.0
beta = 0.0

repeats = 1

a = gd.empty(nbands)
# a = a.reshape((nbands, -1))
b = gd.empty(nbands)
# b = a.reshape((nbands, -1))
# c = np.zeros((a.shape[0], b.shape[0]))

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

t0 = time()
for r in range(repeats):
    gemm(alpha, a, b, beta, c, "c")
t1 = time()
print "Host time", t1 - t0

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import empty [as 别名]
tip.positions[:, 2] = np.arange(5) * d + a / 2
tip.positions[:, :2] = a / 2
calc = GPAW(mode='lcao', basis='sz', width=0.1)
tip.set_calculator(calc)
tip.get_potential_energy()
calc.write('tip')

tip = GPAW('tip')
srf = GPAW('tip')
stm = STM(tip, srf)
stm.initialize(0, dmin=4.0)

from gpaw.spline import Spline
import numpy as np
a = np.array([1, 0.9, 0.1, 0.0])
s = Spline(0, 2.0, a)
p = Spline(1, 2.0, a)
from gpaw.transport.stm import AtomCenteredFunctions as ACF
from gpaw.grid_descriptor import GridDescriptor
a = 10.0
N = 32
gd = GridDescriptor((N, N, N), (a, a, a), (1, 1, 0))

a = ACF(gd, [s], np.array((0.5, 0.5, 0.5)))
b = ACF(gd, [p], np.array((0.6, 0.4, 0.54)))
print a.overlap(a)
print (a|a), (a|b), (b|b)
v = gd.empty()
v[:] = 1.0
print (a|v|a), (a|v|b), (b|v|a), (b|v|b)

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import empty [as 别名]
from gpaw.grid_descriptor import GridDescriptor
from gpaw.spline import Spline
import gpaw.mpi as mpi
from gpaw.wavefunctions.pw import PWDescriptor, PWLFC


x = 2.0
rc = 3.5
r = np.linspace(0, rc, 100)

n = 40
a = 8.0
cell_cv = np.array([[a, 0.5, -1], [0, a, 2], [-1, 0, a + 1]])
gd = GridDescriptor((n, n, n), cell_cv, comm=mpi.serial_comm)

a_R = gd.empty()
z = np.linspace(0, n, n, endpoint=False)
a_R[:] = 2 + np.sin(2 * np.pi * z / n)

spos_ac = np.array([(0.15, 0.45, 0.95)])

pd = PWDescriptor(45, gd)
a_G = pd.fft(a_R)

s = Spline(0, rc, 2 * x**1.5 / np.pi * np.exp(-x * r**2))
p = Spline(1, rc, 2 * x**1.5 / np.pi * np.exp(-x * r**2))
d = Spline(2, rc, 2 * x**1.5 / np.pi * np.exp(-x * r**2))

lfc = PWLFC([[s, p, d]], pd)
lfc.set_positions(spos_ac)
b_LG = pd.zeros(9)

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import empty [as 别名]
# Set up communicators:
domain_comm, kpt_comm, band_comm = distribute_cpus(parsize_domain=D,
                                                   parsize_bands=B,
                                                   nspins=1, nibzkpts=1)
assert kpt_comm.size == 1
if world.rank == 0:
    print 'MPI: %d domains, %d band groups' % (domain_comm.size, band_comm.size)

# Set up band and grid descriptors:
bd = BandDescriptor(N, band_comm, False)
gd = GridDescriptor((G, G, G), (a, a, a), True, domain_comm, parsize=D)
ksl = BandLayouts(gd, bd, float)

# Random wave functions:
psit_mG = gd.empty(M)
for m in range(M):
    np.random.seed(world.rank * M + m)
    psit_mG[m] = np.random.uniform(-0.5, 0.5, tuple(gd.n_c))
if world.rank == 0:
    print 'Size of wave function array:', psit_mG.shape
P_ani = {0: psit_mG[:, :2, 0, 0].copy(),
         1: psit_mG[:, -1, -1, -3:].copy()}
if 0:
    X = M // K
    assert K * X == M
    if G**3 // D // K * K < G**3 // D:
        X += 1
    print X
    work1_xG = gd.empty(X)
    work2_xG = gd.empty(X)