Python GridDescriptor.get_grid_point_coordinates方法代码示例

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import get_grid_point_coordinates [as 别名]
    def get_combined_data(self, qmdata=None, cldata=None, spacing=None):
        
        if qmdata is None:
            qmdata = self.density.rhot_g
        
        if cldata is None:
            cldata = self.classical_material.charge_density
        
        if spacing is None:
            spacing = self.cl.gd.h_cv[0, 0]
        
        spacing_au = spacing / Bohr  # from Angstroms to a.u.
        
        # Collect data from different processes
        cln = self.cl.gd.collect(cldata)
        qmn = self.qm.gd.collect(qmdata)

        clgd = GridDescriptor(self.cl.gd.N_c,
                              self.cl.cell,
                              False,
                              serial_comm,
                              None)

        if world.rank == 0:
            cln *= self.classical_material.sign
            # refine classical part
            while clgd.h_cv[0, 0] > spacing_au * 1.50:  # 45:
                cln = Transformer(clgd, clgd.refine()).apply(cln)
                clgd = clgd.refine()
                
            # refine quantum part
            qmgd = GridDescriptor(self.qm.gd.N_c,
                                  self.qm.cell,
                                  False,
                                  serial_comm,
                                  None)                           
            while qmgd.h_cv[0, 0] < clgd.h_cv[0, 0] * 0.95:
                qmn = Transformer(qmgd, qmgd.coarsen()).apply(qmn)
                qmgd = qmgd.coarsen()
            
            assert np.all(qmgd.h_cv == clgd.h_cv), " Spacings %.8f (qm) and %.8f (cl) Angstroms" % (qmgd.h_cv[0][0] * Bohr, clgd.h_cv[0][0] * Bohr)
            
            # now find the corners
            r_gv_cl = clgd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
            cind = self.qm.corner1 / np.diag(clgd.h_cv) - 1
            
            n = qmn.shape

            # print 'Corner points:     ', self.qm.corner1*Bohr,      ' - ', self.qm.corner2*Bohr
            # print 'Calculated points: ', r_gv_cl[tuple(cind)]*Bohr, ' - ', r_gv_cl[tuple(cind+n+1)]*Bohr
                        
            cln[cind[0] + 1:cind[0] + n[0] + 1,
                cind[1] + 1:cind[1] + n[1] + 1,
                cind[2] + 1:cind[2] + n[2] + 1] += qmn
        
        world.barrier()
        return cln, clgd

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import get_grid_point_coordinates [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 get_grid_point_coordinates [as 别名]
if size == 1:
    for name, D, cell in cells:
        if name == 'Jelver':
            # Strange one!
            continue
            
        print('------------------')
        print(name, D)
        print(cell[0])
        print(cell[1])
        print(cell[2])
        for n in range(1, 6):
            N = 2 * n + 2
            gd = GridDescriptor((N, N, N), cell)
            b_g = gd.zeros()
            r_gv = gd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
            c_v = gd.cell_cv.sum(0) / 2
            r_gv -= c_v
            lap = Laplace(gd, n=n)
            grad_v = [Gradient(gd, v, n=n) for v in range(3)]
            assert lap.npoints == D * 2 * n + 1
            for m in range(0, 2 * n + 1):
                for ix in range(m + 1):
                    for iy in range(m - ix + 1):
                        iz = m - ix - iy
                        a_g = (r_gv**(ix, iy, iz)).prod(3)
                        if ix + iy + iz == 2 and max(ix, iy, iz) == 2:
                            r = 2.0
                        else:
                            r = 0.0
                        lap.apply(a_g, b_g)

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import get_grid_point_coordinates [as 别名]
    def get_eigenmodes(self,filename = None, chi0 = None, calc = None, dm = None, 
                       xc = 'RPA', sum = None, vcut = None, checkphase = False, 
                       return_full = False):
        """
        Calculate the plasmonic eigenmodes as eigenvectors of the dielectric matrix.  

        Parameters:

        filename:  pckl file
                   output from response calculation.
         
        chi0:  gpw file
               chi0_wGG from response calculation.

        calc:  gpaw calculator instance
               ground state calculator used in response calculation.
               Wavefunctions only needed if chi0 is calculated from scratch

        dm:  gpw file
             dielectric matrix from response calculation

        xc:  str 'RPA'or 'ALDA' XC- Kernel
        
        sum:  str
              '2D': sum in the x and y directions
              '1D': To be implemented

        vcut:  str '0D','1D' or '2D'
               Cut the Coulomb potential 

        checkphase:   Bool
                      if True, the eigenfunctions id rotated in the complex
                      plane, to be made as real as posible

        return_full:  Bool
                      if True, the eigenvectors in reciprocal space is also
                      returned. 
           
        """
        self.read(filename)
        self.pbc = [1,1,1]
        #self.calc.atoms.pbc = [1,1,1]
        npw = self.npw
        self.w_w = np.linspace(0, self.dw * (self.Nw - 1)*Hartree, self.Nw)
        self.vcut = vcut
        dm_wGG = self.get_dielectric_matrix(xc=xc,
                                            symmetric=False,
                                            chi0_wGG=chi0,
                                            calc=calc,
                                            vcut=vcut)
    
        q = self.q_c
        

        # get grid on which the eigenmodes are calculated
        #gd = self.calc.wfs.gd
        #r = gd.get_grid_point_coordinates()
        #rrr = r*Bohr 
        from gpaw.utilities.gpts import get_number_of_grid_points
        from gpaw.grid_descriptor import GridDescriptor
        grid_size = [1,1,1]
        h=0.2
        cell_cv = self.acell_cv*np.diag(grid_size)
        mode = 'fd'
        realspace = True
        h /= Bohr
        N_c = get_number_of_grid_points(cell_cv, h, mode, realspace)
        gd = GridDescriptor(N_c, cell_cv, self.pbc) 
        #gd = self.calc.wfs.gd
        r = gd.get_grid_point_coordinates()
        rrr = r*Bohr
        
        eig_0 = np.array([], dtype = complex)
        eig_left = np.array([], dtype = complex)
        eig_right = np.array([], dtype = complex)
        vec_modes = np.zeros([1, self.npw], dtype = complex)
        vec_modes_dual = np.zeros([1, self.npw], dtype = complex)
        vec_modes_density = np.zeros([1, self.npw], dtype = complex)
        vec_modes_norm = np.zeros([1, self.npw], dtype = complex)
        eig_all = np.zeros([1, self.npw], dtype = complex)
        eig_dummy = np.zeros([1, self.npw], dtype = complex)
        v_dummy = np.zeros([1, self.npw], dtype = complex)
        vec_dummy = np.zeros([1, self.npw], dtype = complex)
        vec_dummy2 = np.zeros([1, self.npw], dtype = complex)
        w_0 = np.array([]) 
        w_left = np.array([])
        w_right = np.array([])
     
        if sum == '2D':
            v_ind = np.zeros([1, r.shape[-1]], dtype = complex)
            n_ind = np.zeros([1, r.shape[-1]], dtype = complex)
        elif sum == '1D':            
            self.printtxt('1D sum not implemented')
            return 
        else:
            v_ind = np.zeros([1, r.shape[1], r.shape[2], r.shape[3]], dtype = complex)
            n_ind = np.zeros([1, r.shape[1], r.shape[2], r.shape[3]], dtype = complex)

        eps_GG_plus = dm_wGG[0]
        eig_plus, vec_plus = np.linalg.eig(eps_GG_plus)  # find eigenvalues and eigenvectors
#.........这里部分代码省略.........

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

# 需要导入模块: from gpaw.grid_descriptor import GridDescriptor [as 别名]
# 或者: from gpaw.grid_descriptor.GridDescriptor import get_grid_point_coordinates [as 别名]
    def initialize(self):
                        
        self.eta /= Hartree
        self.ecut /= Hartree

        calc = self.calc
        
        # kpoint init
        self.kd = kd = calc.wfs.kd
        self.bzk_kc = kd.bzk_kc
        self.ibzk_kc = kd.ibzk_kc
        self.nkpt = kd.nbzkpts
        self.ftol /= self.nkpt

        # band init
        if self.nbands is None:
            self.nbands = calc.wfs.nbands
        self.nvalence = calc.wfs.nvalence

        # cell init
        self.acell_cv = calc.atoms.cell / Bohr
        self.bcell_cv, self.vol, self.BZvol = get_primitive_cell(self.acell_cv)

        # grid init
        self.nG = calc.get_number_of_grid_points()
        self.nG0 = self.nG[0] * self.nG[1] * self.nG[2]
        gd = GridDescriptor(self.nG, calc.wfs.gd.cell_cv, pbc_c=True, comm=serial_comm)
        self.gd = gd        
        self.h_cv = gd.h_cv

        # obtain eigenvalues, occupations
        nibzkpt = kd.nibzkpts
        kweight_k = kd.weight_k

        try:
            self.e_kn
        except:
            self.printtxt('Use eigenvalues from the calculator.')
            self.e_kn = np.array([calc.get_eigenvalues(kpt=k)
                    for k in range(nibzkpt)]) / Hartree
            self.printtxt('Eigenvalues(k=0) are:')
            print  >> self.txt, self.e_kn[0] * Hartree
        self.f_kn = np.array([calc.get_occupation_numbers(kpt=k) / kweight_k[k]
                    for k in range(nibzkpt)]) / self.nkpt


        # k + q init
        assert self.q_c is not None
        self.qq_v = np.dot(self.q_c, self.bcell_cv) # summation over c

        if self.optical_limit:
            kq_k = np.arange(self.nkpt)
            self.expqr_g = 1.
        else:
            r_vg = gd.get_grid_point_coordinates() # (3, nG)
            qr_g = gemmdot(self.qq_v, r_vg, beta=0.0)
            self.expqr_g = np.exp(-1j * qr_g)
            del r_vg, qr_g
            kq_k = kd.find_k_plus_q(self.q_c)
        self.kq_k = kq_k

        # Plane wave init
        self.npw, self.Gvec_Gc, self.Gindex_G = set_Gvectors(self.acell_cv, self.bcell_cv, self.nG, self.ecut)

        # Projectors init
        setups = calc.wfs.setups
        pt = LFC(gd, [setup.pt_j for setup in setups],
                 dtype=calc.wfs.dtype, forces=True)
        spos_ac = calc.atoms.get_scaled_positions()
        pt.set_k_points(self.bzk_kc)
        pt.set_positions(spos_ac)
        self.pt = pt

        # Printing calculation information
        self.print_stuff()

        return