Update examples

Author: sunqm

examples/pbc/31-pbc_0D_as_mol.py

@@ -45,16 +45,20 @@
 print('E(HF) of molecular RHF with cell %s' % mf.e_tot)
 
 #
-# Convert mol back to cell.
+# Convert mol back to cell using the to_cell function.
+# Lattice vectors "a" are not available in the mole object. Specify a "box" to
+# place the mole in the cell.
 #
-# The mol ojbect above contains all information of the pbc system which was
-# initialized at the beginning. Using the "view" method to convert mol back to
-# the cell object, all information can be transfer to the resultant cell
-# object. Lattice vectors "a" are not available in the mole object. It needs
-# to be specified in the cell.
-#
-cell_0D = mol.view(pbcgto.Cell)
-cell_0D.a = numpy.eye(3)
-cell_0D.dimension = 0
+cell_0D = mol.to_cell(box=a, dimension=0)
 mf = pbcscf.RHF(cell).density_fit().run()
 print('E(HF) with 0D PBC RHF calculation %s' % mf.e_tot)
+
+#
+# By transforming Mole to Cell instance, we can apply the MultiGrid integral
+# algorithm for DFT calculations to fastly evaluate the Coulomb and XC functional.
+#
+cell_0D.verbose = 5
+mf = cell_0D.RKS(xc='pbe')
+mf.run() # This calls the standard numint module
+mf = mf.multigrid_numint()
+mf.run() # This calls the MultiGridNumInt algorithm

pyscf/gto/mole.py

@@ -3849,14 +3849,13 @@ def to_cell(self, box=None, dimension=None, margin=None):
             dimension : integer
                 PBC dimensions
             margin: float
-                The distance from the edge of the molecule to the edge of the box. 
+                The distance from the edge of the molecule to the edge of the box.
                 If not provided, a default margin will be estimated, to ensure
                 that the electron density decays to approximately 1e-7 outside
                 of the box.
         '''
         from pyscf.pbc.gto import Cell, rcut_by_shells
-        cell = Cell()
-        cell.__dict__.update(self.__dict__)
+        cell = mol.view(Cell)
         if box is None:
             # Place molecule in a big box for dimension=0
             dimension = 0

pyscf/pbc/gto/cell.py

@@ -746,25 +746,22 @@ def ewald(cell, ew_eta=None, ew_cut=None):
     # where
     #   ZS_I(G) = \sum_a Z_a exp (i G.R_a)
 
-    Gv, Gvbase, weights = cell.get_Gv_weights(mesh)
-    absG2 = np.einsum('gi,gi->g', Gv, Gv)
-    absG2[absG2==0] = 1e200
+    if cell.dimension == 3 or cell.low_dim_ft_type == 'inf_vacuum':
+        Gv, Gvbase, weights = cell.get_Gv_weights(mesh)
+        absG2 = np.einsum('gi,gi->g', Gv, Gv)
+        absG2[absG2==0] = 1e200
 
-    if cell.dimension != 2 or cell.low_dim_ft_type == 'inf_vacuum':
-        # TODO: truncated Coulomb for 0D. The non-uniform grids for inf-vacuum
-        # have relatively large error
         coulG = 4*np.pi / absG2
         coulG *= weights
 
         #:ZSI = np.einsum('i,ij->j', chargs, cell.get_SI(Gv))
-        ngrids = len(Gv)
-        ZSI = np.empty((ngrids,), dtype=np.complex128)
-        mem_avail = cell.max_memory - lib.current_memory()[0]
-        blksize = int((mem_avail*1e6 - cell.natm*24)/((3+cell.natm*2)*8))
-        blksize = min(ngrids, max(mesh[2], blksize))
-        for ig0, ig1 in lib.prange(0, ngrids, blksize):
-            np.einsum('i,ij->j', chargs, cell.get_SI(Gv[ig0:ig1]), out=ZSI[ig0:ig1])
-
+        basex, basey, basez = cell.get_Gv_weights(mesh)[1]
+        b = cell.reciprocal_vectors()
+        rb = np.dot(coords, b.T)
+        SIx = np.exp(-1j*np.einsum('z,g->zg', rb[:,0], basex))
+        SIy = np.exp(-1j*np.einsum('z,g->zg', rb[:,1], basey))
+        SIz = np.exp(-1j*np.einsum('z,g->zg', rb[:,2], basez))
+        ZSI = np.einsum('i,ix,iy,iz->xyz', chargs, SIx, SIy, SIz).ravel()
         ZexpG2 = ZSI * np.exp(-absG2/(4*ew_eta**2))
         ewg = .5 * np.einsum('i,i,i', ZSI.conj(), ZexpG2, coulG).real
 
@@ -789,6 +786,8 @@ def gn0(eta,z):
         inv_area = np.linalg.norm(np.cross(b[0], b[1]))/(2*np.pi)**2
         # Perform the reciprocal space summation over  all reciprocal vectors
         # within the x,y plane.
+        Gv, Gvbase, weights = cell.get_Gv_weights(mesh)
+        absG2 = np.einsum('gi,gi->g', Gv, Gv)
         planarG2_idx = np.logical_and(Gv[:,2] == 0, absG2 > 0.0)
         Gv = Gv[planarG2_idx]
         absG2 = absG2[planarG2_idx]