Merge branch 'issue_86_errhandle' into dev

my_pyscf/mcscf/productstate.py

@@ -1,6 +1,7 @@
 import numpy as np
 from scipy import linalg
 from pyscf import lib
+from pyscf.scf.addons import canonical_orth_
 from pyscf.fci import cistring
 from pyscf.mcscf.addons import state_average as state_average_mcscf
 from mrh.my_pyscf.fci.csf import CSFFCISolver
@@ -30,10 +31,12 @@ def kernel (self, h1, h2, norb_f, nelec_f, ecore=0, ci0=None, orbsym=None,
         converged = False
         e_sigma = 0
         e = [0 for n in norb_f]
-        ci1 = self.get_init_guess (ci0, norb_f, nelec_f, h1, h2)
+        ci1 = ci0
         log.info ('Entering product-state fixed-point CI iteration')
         for it in range (max_cycle_macro):
-            ci0 = ci1
+            ci0 = self.get_init_guess (ci1, norb_f, nelec_f, h1, h2)
+            # Issue #86: put get_init_guess INSIDE the iteration in case _1shot below encounters
+            # linear dependencies and can't populate all CI vectors for all fragments.
             h1eff, h0eff, ci0 = self.project_hfrag (h1, h2, ci0, norb_f, nelec_f,
                 ecore=ecore, **kwargs)
             grad = self._get_grad (h1eff, h2, ci0, norb_f, nelec_f, **kwargs)
@@ -55,12 +58,35 @@ def kernel (self, h1, h2, norb_f, nelec_f, ecore=0, ci0=None, orbsym=None,
         conv_str = ['NOT converged','converged'][int (converged)]
         log.info (('Product_state fixed-point CI iteration {} after {} '
                    'cycles').format (conv_str, it))
-        if not converged: self._debug_csfs (log, ci0, ci1, norb_f, nelec_f, grad)
+        if not converged:
+            ci1 = self.get_init_guess (ci1, norb_f, nelec_f, h1, h2)
+            # Issue #86: see above, same problem
+            self._debug_csfs (log, ci0, ci1, norb_f, nelec_f, grad)
         energy_elec = self.energy_elec (h1, h2, ci1, norb_f, nelec_f,
             ecore=ecore, **kwargs)
         return converged, energy_elec, ci1
 
     def get_init_guess (self, ci0, norb_f, nelec_f, h1, h2):
+        '''Generate CI guess vectors for all fragments.
+
+        Args:
+            ci0: list of length nfrag or None
+                Contains either None or ndarrays of guess CI vectors. Any new guess CI vectors
+                constructed by this function are constrained to be orthogonal to those already
+                provided here, if any.
+            norb_f: list of length nfrag of integers
+                Number of orbitals in each fragment
+            nelec_f: list of length nfrag of integers
+                Number of electrons (in reference state) in each fragment
+            h1: ndarray of shape (ncas,ncas) or (2,ncas,ncas)
+                One-electron Hamiltonian amplitudes
+            h2: ndarray of shape (ncas,ncas,ncas,ncas)
+                Two-electron Hamiltonian amplitudes
+
+        Returns:
+            ci1: list of length nfrag of ndarrays
+                Orthonormal guess CI vectors. Any vectors present in ci0 are preserved unaltered.
+        '''
         if ci0 is None: ci0 = [None for i in range (len (norb_f))]
         ci1 = [c for c in ci0] # reference safety
         if h1.ndim < 3: h1 = np.stack ([h1, h1], axis=0)
@@ -78,13 +104,21 @@ def get_init_guess (self, ci0, norb_f, nelec_f, h1, h2):
             elif np.asarray (ci1[ix]).reshape (-1,na*nb).shape[0] < solver.nroots:
                 ci1_inp = np.asarray (ci1[ix]).reshape (-1,na*nb)
                 ci1_guess = np.asarray (ci1_guess).reshape (-1,na*nb)
-                ovlp = ci1_inp.conj () @ ci1_guess.T
-                ci1_guess -= ovlp.T @ ci1_inp
-                ovlp = ci1_guess.conj () @ ci1_guess.T
-                Q, R = linalg.qr (ovlp)
-                ci1_guess = Q.T @ ci1_guess
-                ci1[ix] = np.append (ci1_inp, ci1_guess, axis=0)[:solver.nroots].reshape (
-                    solver.nroots, na, nb)
+                x = np.append (ci1_inp, ci1_guess, axis=0)
+                x = canonical_orth_(x.conj () @ x.T).T @ x
+                # ^ an orthonormal basis
+                assert (x.shape[0] >= solver.nroots)
+                x2inp = x.conj () @ ci1_inp.T
+                ninp = ci1_inp.shape[0]
+                u, svals, vh = linalg.svd (x2inp, full_matrices=True)
+                u[:,:ninp] = u[:,:ninp] @ vh
+                ci1_new = u.T @ x
+                nnew = ci1_new.shape[0]
+                ovlp = ci1_new.conj () @ ci1_inp.T
+                assert (np.all (np.abs (ovlp[:ninp,:ninp] - np.eye (ninp)) < 1e-6)), '{}'.format (ovlp)
+                ovlp = (ci1_new.conj () @ ci1_new.T)
+                assert (np.all (np.abs (ovlp - np.eye (nnew)) < 1e-6)), '{}'.format (ovlp)
+                ci1[ix] = ci1_new[:solver.nroots].reshape (solver.nroots, na, nb)
         return self._check_init_guess (ci1, norb_f, nelec_f)
 
     def _check_init_guess (self, ci0, norb_f, nelec_f):