Support linear dependent overlap matrix for k point sampling (#572)

* Linear dependent overlap matrix roughly supported for k point sampling

* Forgot to add one file

* Bug fix with previous merge and previous functionality

* Another round of cleanup

* Tests done!

Author: henryw7

gpu4pyscf/df/df_jk.py

@@ -134,7 +134,8 @@ def get_jk(self, mol=None, dm=None, hermi=1, with_j=True, with_k=True,
             vj, vk = self.with_df.get_jk(dm, hermi, with_j, with_k,
                                          self.direct_scf_tol, omega)
         else:
-            vj, vk = super().get_jk(mol, dm, hermi, with_j, with_k, omega)
+            raise ValueError(f"with_df field not found in a df object (type = {type(self)}) during a get_jk() call.")
+            # vj, vk = super().get_jk(mol, dm, hermi, with_j, with_k, omega)
         return vj, vk
 
     def nuc_grad_method(self):

gpu4pyscf/lib/cupy_helper.py

@@ -836,7 +836,7 @@ def pinv(a, lindep=1e-10):
     j2c = cupy.dot(v1/w[mask], v1.conj().T)
     return j2c
 
-def cond(a, sympos=False):
+def cond(a, sympos=False, verbose=logger.WARN):
     """
     Calculate the condition number of a matrix.
 
@@ -847,27 +847,33 @@ def cond(a, sympos=False):
     Returns:
     float: The condition number of the matrix.
     """
+    if isinstance(verbose, logger.Logger):
+        log = verbose
+    else:
+        log = logger.Logger(sys.stdout, verbose)
+
     if a.shape[0] > cusolver.MAX_EIGH_DIM:
         if not SCIPY_EIGH_FOR_LARGE_ARRAYS:
             raise RuntimeError(
                 f'Array size exceeds the maximum size {cusolver.MAX_EIGH_DIM}.')
         a = a.get()
         if sympos:
             s = scipy.linalg.eigvalsh(a)
-            if s[0] <= 0:
-                raise RuntimeError('matrix is not positive definite')
-            return s[-1] / s[0]
-        else:
-            _, s, _ = scipy.linalg.svd(a)
-            cond_number = s[0] / s[-1]
-            return cond_number
-
-    if sympos:
-        s = cupy.linalg.eigvalsh(a)
-        if s[0] <= 0:
-            raise RuntimeError('matrix is not positive definite')
-        return s[-1] / s[0]
+            if s[0] > 0:
+                return s[-1] / s[0]
+            else:
+                log.warn(f'In condition number calculation, matrix is assumed to be positive definite, but it is not (minimal eigenvalue = {s[0]:e})')
+        _, s, _ = scipy.linalg.svd(a)
+        cond_number = s[0] / s[-1]
+        return cond_number
+
     else:
+        if sympos:
+            s = cupy.linalg.eigvalsh(a)
+            if s[0] > 0:
+                return s[-1] / s[0]
+            else:
+                log.warn(f'In condition number calculation, matrix is assumed to be positive definite, but it is not (minimal eigenvalue = {s[0]:e})')
         _, s, _ = cupy.linalg.svd(a)
         cond_number = s[0] / s[-1]
         return cond_number

gpu4pyscf/pbc/dft/krks.py

@@ -59,6 +59,7 @@ def get_veff(ks, cell=None, dm=None, dm_last=None, vhf_last=None, hermi=1,
         n, exc, vxc = ni.nr_rks(cell, ks.grids, ks.xc, dm, 0, hermi, kpts, kpts_band)
         log.debug('nelec by numeric integration = %s', n)
         if ks.do_nlc():
+            raise NotImplementedError("VV10 not implemented for periodic system")
             if ni.libxc.is_nlc(ks.xc):
                 xc = ks.xc
             else:
@@ -73,12 +74,12 @@ def get_veff(ks, cell=None, dm=None, dm_last=None, vhf_last=None, hermi=1,
     vj, vk = _get_jk(ks, cell, dm, hermi, kpts, kpts_band, not j_in_xc,
                      dm_last, vhf_last)
     if not j_in_xc:
-        vxc += vj
+        vxc = vxc + vj
         ecoul = None
         if ground_state:
             ecoul = float(cp.einsum('Kij,Kji->', dm, vj).real.get()) * .5 * weight
     if hybrid:
-        vxc -= .5 * vk
+        vxc = vxc - .5 * vk
         if ground_state:
             exc -= float(cp.einsum('Kij,Kji->', dm, vk).real.get()) * .25 * weight
     vxc = tag_array(vxc, ecoul=ecoul, exc=exc, vj=vj, vk=vk)

gpu4pyscf/pbc/dft/kuks.py

@@ -60,6 +60,7 @@ def get_veff(ks, cell=None, dm=None, dm_last=None, vhf_last=None, hermi=1,
         ks.initialize_grids(cell, dm, kpts)
         n, exc, vxc = ni.nr_uks(cell, ks.grids, ks.xc, dm, 0, hermi, kpts, kpts_band)
         if ks.do_nlc():
+            raise NotImplementedError("VV10 not implemented for periodic system")
             if ni.libxc.is_nlc(ks.xc):
                 xc = ks.xc
             else:
@@ -75,12 +76,12 @@ def get_veff(ks, cell=None, dm=None, dm_last=None, vhf_last=None, hermi=1,
     vj, vk = krks._get_jk(ks, cell, dm, hermi, kpts, kpts_band, not j_in_xc,
                           dm_last, vhf_last)
     if not j_in_xc:
-        vxc += vj[0] + vj[1]
+        vxc = vxc + vj[0] + vj[1]
         ecoul = None
         if ground_state:
             ecoul = float(cp.einsum('nKij,mKji->', dm, vj).real.get()) * .5 * weight
     if hybrid:
-        vxc -= vk
+        vxc = vxc - vk
         if ground_state:
             exc -= float(cp.einsum('nKij,nKji->', dm, vk).real.get()) * .5 * weight
     vxc = tag_array(vxc, ecoul=ecoul, exc=exc, vj=vj, vk=vk)

gpu4pyscf/pbc/dft/rks.py

@@ -81,6 +81,13 @@ def get_veff(ks, cell=None, dm=None, dm_last=None, vhf_last=None, hermi=1,
         n, exc, vxc = ni.nr_rks(cell, ks.grids, ks.xc, dm, 0, hermi, kpt, kpts_band)
         log.debug('nelec by numeric integration = %s', n)
         if ks.do_nlc():
+            warning_message = "ATTENTION!!! VV10 is only valid for open boundary, and it is incorrect for actual periodic system! " \
+                              "Lattice summation is not performed for the double integration. " \
+                              "Please use only under open boundary, i.e. neighbor images are well separated, and " \
+                              "all atoms belonging to one image is placed in the same image in the input."
+            log.warn(warning_message)
+            print(warning_message) # This is an important warning, so print even if verbose == 0.
+
             if ni.libxc.is_nlc(ks.xc):
                 xc = ks.xc
             else:

gpu4pyscf/pbc/dft/uks.py

@@ -68,6 +68,13 @@ def get_veff(ks, cell=None, dm=None, dm_last=None, vhf_last=None, hermi=1,
         ks.initialize_grids(cell, dm, kpt)
         n, exc, vxc = ni.nr_uks(cell, ks.grids, ks.xc, dm, 0, hermi, kpt, kpts_band)
         if ks.do_nlc():
+            warning_message = "ATTENTION!!! VV10 is only valid for open boundary, and it is incorrect for actual periodic system! " \
+                              "Lattice summation is not performed for the double integration. " \
+                              "Please use only under open boundary, i.e. neighbor images are well separated, and " \
+                              "all atoms belonging to one image is placed in the same image in the input."
+            log.warn(warning_message)
+            print(warning_message) # This is an important warning, so print even if verbose == 0.
+
             if ni.libxc.is_nlc(ks.xc):
                 xc = ks.xc
             else:

gpu4pyscf/pbc/scf/hf.py

@@ -125,8 +125,7 @@ def check_sanity(self):
             isinstance(self.with_df, df.DF)):
             logger.warn(self, 'exxdiv %s is not supported in DF', self.exxdiv)
 
-        if self.verbose >= logger.DEBUG:
-            mol_hf.SCF.check_sanity(self)
+        mol_hf.SCF.check_sanity(self)
         return self
 
     @property

gpu4pyscf/pbc/scf/khf.py

@@ -27,7 +27,7 @@
 from pyscf.pbc import tools
 from gpu4pyscf.lib import logger, utils
 from gpu4pyscf.lib.cupy_helper import (
-    return_cupy_array, contract, tag_array, sandwich_dot)
+    return_cupy_array, contract, tag_array, sandwich_dot, eigh)
 from gpu4pyscf.scf import hf as mol_hf
 from gpu4pyscf.pbc.scf import hf as pbchf
 from gpu4pyscf.pbc import df
@@ -168,6 +168,9 @@ def energy_elec(mf, dm_kpts=None, h1e_kpts=None, vhf_kpts=None):
     return (e1+e_coul).real, e_coul.real
 
 def canonicalize(mf, mo_coeff_kpts, mo_occ_kpts, fock=None):
+    if hasattr(mf, 'overlap_canonical_decomposed_x') and mf.overlap_canonical_decomposed_x is not None:
+        raise NotImplementedError("Overlap matrix canonical decomposition (removing linear dependency for diffused orbitals) "
+                                  "not supported for canonicalize() function with k-point sampling")
     if fock is None:
         dm = mf.make_rdm1(mo_coeff_kpts, mo_occ_kpts)
         fock = mf.get_fock(dm=dm)
@@ -389,10 +392,26 @@ def eig(self, h_kpts, s_kpts):
         nkpts, nao = h_kpts.shape[:2]
         eig_kpts = cp.empty((nkpts, nao))
         mo_coeff_kpts = cp.empty((nkpts, nao, nao), dtype=h_kpts.dtype)
-        for k in range(nkpts):
-            e, c = self._eigh(h_kpts[k], s_kpts[k])
-            eig_kpts[k] = e
-            mo_coeff_kpts[k] = c
+
+        x_kpts = None
+        if hasattr(self, 'overlap_canonical_decomposed_x') and self.overlap_canonical_decomposed_x is not None:
+            x_kpts = [cp.asarray(x) for x in self.overlap_canonical_decomposed_x]
+
+        if x_kpts is None:
+            for k in range(nkpts):
+                e, c = eigh(h_kpts[k], s_kpts[k])
+                eig_kpts[k] = e
+                mo_coeff_kpts[k] = c
+        else:
+            for k in range(nkpts):
+                xk = x_kpts[k]
+                ek, ck = cp.linalg.eigh(xk.T.conj() @ h_kpts[k] @ xk)
+                ck = xk @ ck
+                _, nmo_k = xk.shape
+                eig_kpts[k, :nmo_k] = ek
+                eig_kpts[k, nmo_k:] = float(cp.max(cp.abs(ek))) * 2 + 1e5
+                mo_coeff_kpts[k, :, :nmo_k] = ck
+                mo_coeff_kpts[k, :, nmo_k:] = 0
         return eig_kpts, mo_coeff_kpts
 
     def make_rdm1(self, mo_coeff_kpts=None, mo_occ_kpts=None, **kwargs):

gpu4pyscf/pbc/scf/kuhf.py

@@ -174,6 +174,9 @@ def canonicalize(mf, mo_coeff_kpts, mo_occ_kpts, fock=None):
     '''Canonicalization diagonalizes the UHF Fock matrix within occupied,
     virtual subspaces separatedly (without change occupancy).
     '''
+    if hasattr(mf, 'overlap_canonical_decomposed_x') and mf.overlap_canonical_decomposed_x is not None:
+        raise NotImplementedError("Overlap matrix canonical decomposition (removing linear dependency for diffused orbitals) "
+                                  "not supported for canonicalize() function with k-point sampling")
     if fock is None:
         dm = mf.make_rdm1(mo_coeff_kpts, mo_occ_kpts)
         fock = mf.get_fock(dm=dm)