update dump_mo

README.md

@@ -12,7 +12,8 @@ This program aims to do:
 
 ## Features
 * UHF -> UNO -> CASSCF
-* RHF -> PM LMO -> CASSCF
+* RHF -> vir MO projection -> PM LMO -> CASSCF
+
 UHF, RHF can be auto-detected.
 
 ## Utilities
@@ -22,6 +23,7 @@ UHF, RHF can be auto-detected.
   + from_fch
 * UHF stability
 * dump CASCI coefficients
+* dump active orbital compositions
 
 ## Quick Start
 ```

autocas.py

@@ -27,7 +27,7 @@ def get_uno(mf, st='st2'):
     elif st=='st2':
         dm = mf.make_rdm1()
         unos, noon = get_uno_st2(dm[0] + dm[1], S)
-        nacto = check_uno(noon)
+        nacto, ndb, nex = check_uno(noon)
     print('UNO ON:', noon)
     #ndb, nocc, nopen = idx
     #nacto = nocc - ndb
@@ -36,13 +36,13 @@ def get_uno(mf, st='st2'):
     print('nacto, nacta, nactb: %d %d %d' % (nacto, nacta, nactb))
     mf = mf.to_rhf()
     mf.mo_coeff = unos
-    return mf, unos, noon, nacto, (nacta, nactb)
+    return mf, unos, noon, nacto, (nacta, nactb), ndb, nex
 
 def check_uno(noon, thresh=1.98):
-    n1 = np.count_nonzero(noon < thresh)
-    n2 = np.count_nonzero(noon > (2.0-thresh))
-    nacto = n1 + n2 - len(noon)
-    return nacto
+    ndb = np.count_nonzero(noon > thresh)
+    nex = np.count_nonzero(noon < (2.0-thresh))
+    nacto = len(noon) - ndb - nex
+    return nacto, ndb, nex
 
 def get_uno_st2(dm, S):
     A = reduce(np.dot, (S, dm, S))
@@ -141,7 +141,7 @@ def get_locorb(mf, localize='pm1', pair=True):
         mf.mo_coeff = alpha_coeff.copy()
         print('MOs after pairing')
         dump_mat.dump_mo(mf.mol,mf.mo_coeff[:,idx1:idx3], ncol=10)
-    return mf, alpha_coeff, npair
+    return mf, alpha_coeff, npair, ncore
 
 def check_uhf(mf):
     dm = mf.make_rdm1()
@@ -162,9 +162,9 @@ def check_uhf(mf):
 def cas(mf, crazywfn=False, max_memory=2000):
     is_uhf, mf = check_uhf(mf)
     if is_uhf:
-        mf, unos, unoon, nacto, (nacta, nactb) = get_uno(mf)
+        mf, unos, unoon, nacto, (nacta, nactb), ndb, nex = get_uno(mf)
     else:
-        mf, lmos, npair = get_locorb(mf)
+        mf, lmos, npair, ndb = get_locorb(mf)
         nacto = npair*2
         nacta = nactb = npair
     nopen = nacta - nactb
@@ -182,8 +182,11 @@ def cas(mf, crazywfn=False, max_memory=2000):
     else:
         mc.fcisolver.max_cycle = 100
     mc.natorb = True
-    mc.verbose = 5
+    mc.verbose = 4
     mc.kernel()
+    #mc.analyze(with_meta_lowdin=False)
+    print('Natrual Orbs')
+    dump_mat.dump_mo(mf.mol,mf.mo_coeff[:,ndb:ndb+nacto], ncol=10)
     return mc
 
 def nevpt2(mc):

cidump.py

@@ -30,6 +30,7 @@ def dump(mc, thresh=1e-2):
     Gaussian format: 11a00
     ORCA format    : 22100
     '''
+    print('***** CI components ******')
     ci = mc.ci
     ncas = mc.ncas
     na, nb = mc.nelecas
@@ -59,5 +60,6 @@ def dump(mc, thresh=1e-2):
     print(" c**2     ORCA-type vector")
     for k,v in dump_o:
         print("{: .6f}  {:7s}".format(v,k))
+
     return dump_g, dump_o
 

dump_mat.py

@@ -1,7 +1,7 @@
 import sys
 
 def dump_mo(mol, c, label2=None,
-             ncol=6, digits=5, start=0):
+             ncol=6, digits=3, start=0):
     ''' Print an array in rectangular format
     Args:
         mol: `gto.Mole` object
@@ -40,8 +40,114 @@ def dump_mo(mol, c, label2=None,
             coefs = ''
             for j in v:
                 if abs(j) < 0.1:
-                    coefs += '          '
+                    coefs += '        '
                 else:
                     coefs += fmt % j
+            if len(coefs.strip()) == 0:
+                continue
             stdout.write(('%12s' % label[k]) + coefs + '\n')
-    stdout.flush()
+    stdout.flush()
+
+'''
+def analyze(casscf, mo_coeff=None, ci=None, verbose=None,
+            large_ci_tol=LARGE_CI_TOL, with_meta_lowdin=WITH_META_LOWDIN,
+            **kwargs):
+    from pyscf.lo import orth
+    from pyscf.tools import dump_mat
+    from pyscf.mcscf import addons
+    log = logger.new_logger(casscf, verbose)
+
+    if mo_coeff is None: mo_coeff = casscf.mo_coeff
+    if ci is None: ci = casscf.ci
+    nelecas = casscf.nelecas
+    ncas = casscf.ncas
+    ncore = casscf.ncore
+    nocc = ncore + ncas
+    mocore = mo_coeff[:,:ncore]
+    mocas = mo_coeff[:,ncore:nocc]
+
+    label = casscf.mol.ao_labels()
+    if (isinstance(ci, (list, tuple, RANGE_TYPE)) and
+        not isinstance(casscf.fcisolver, addons.StateAverageFCISolver)):
+        log.warn('Mulitple states found in CASCI/CASSCF solver. Density '
+                 'matrix of the first state is generated in .analyze() function.')
+        civec = ci[0]
+    else:
+        civec = ci
+    if getattr(casscf.fcisolver, 'make_rdm1s', None):
+        casdm1a, casdm1b = casscf.fcisolver.make_rdm1s(civec, ncas, nelecas)
+        casdm1 = casdm1a + casdm1b
+        dm1b = numpy.dot(mocore, mocore.conj().T)
+        dm1a = dm1b + reduce(numpy.dot, (mocas, casdm1a, mocas.conj().T))
+        dm1b += reduce(numpy.dot, (mocas, casdm1b, mocas.conj().T))
+        dm1 = dm1a + dm1b
+        if log.verbose >= logger.DEBUG2:
+            log.info('alpha density matrix (on AO)')
+            dump_mat.dump_tri(log.stdout, dm1a, label, **kwargs)
+            log.info('beta density matrix (on AO)')
+            dump_mat.dump_tri(log.stdout, dm1b, label, **kwargs)
+    else:
+        casdm1 = casscf.fcisolver.make_rdm1(civec, ncas, nelecas)
+        dm1a = (numpy.dot(mocore, mocore.conj().T) * 2 +
+                reduce(numpy.dot, (mocas, casdm1, mocas.conj().T)))
+        dm1b = None
+        dm1 = dm1a
+
+    if log.verbose >= logger.INFO:
+        ovlp_ao = casscf._scf.get_ovlp()
+        # note the last two args of ._eig for mc1step_symm
+        occ, ucas = casscf._eig(-casdm1, ncore, nocc)
+        log.info('Natural occ %s', str(-occ))
+        mocas = numpy.dot(mocas, ucas)
+        if with_meta_lowdin:
+            log.info('Natural orbital (expansion on meta-Lowdin AOs) in CAS space')
+            orth_coeff = orth.orth_ao(casscf.mol, 'meta_lowdin', s=ovlp_ao)
+            mocas = reduce(numpy.dot, (orth_coeff.conj().T, ovlp_ao, mocas))
+        else:
+            log.info('Natural orbital (expansion on AOs) in CAS space')
+        dump_mat.dump_rec(log.stdout, mocas, label, start=1, **kwargs)
+        if log.verbose >= logger.DEBUG2:
+            if not casscf.natorb:
+                log.debug2('NOTE: mc.mo_coeff in active space is different to '
+                           'the natural orbital coefficients printed in above.')
+            if with_meta_lowdin:
+                c = reduce(numpy.dot, (orth_coeff.conj().T, ovlp_ao, mo_coeff))
+                log.debug2('MCSCF orbital (expansion on meta-Lowdin AOs)')
+            else:
+                c = mo_coeff
+                log.debug2('MCSCF orbital (expansion on AOs)')
+            dump_mat.dump_rec(log.stdout, c, label, start=1, **kwargs)
+
+        if casscf._scf.mo_coeff is not None:
+            addons.map2hf(casscf, casscf._scf.mo_coeff)
+
+        if (ci is not None and
+            (getattr(casscf.fcisolver, 'large_ci', None) or
+             getattr(casscf.fcisolver, 'states_large_ci', None))):
+            log.info('** Largest CI components **')
+            if isinstance(ci, (list, tuple, RANGE_TYPE)):
+                if hasattr(casscf.fcisolver, 'states_large_ci'):
+                    # defined in state_average_mix_ mcscf object
+                    res = casscf.fcisolver.states_large_ci(ci, casscf.ncas, casscf.nelecas,
+                                                           large_ci_tol, return_strs=False)
+                else:
+                    res = [casscf.fcisolver.large_ci(civec, casscf.ncas, casscf.nelecas,
+                                                     large_ci_tol, return_strs=False)
+                           for civec in ci]
+                for i, civec in enumerate(ci):
+                    log.info('  [alpha occ-orbitals] [beta occ-orbitals]  state %-3d CI coefficient', i)
+                    for c,ia,ib in res[i]:
+                        log.info('  %-20s %-30s %.12f', ia, ib, c)
+            else:
+                log.info('  [alpha occ-orbitals] [beta occ-orbitals]            CI coefficient')
+                res = casscf.fcisolver.large_ci(ci, casscf.ncas, casscf.nelecas,
+                                                large_ci_tol, return_strs=False)
+                for c,ia,ib in res:
+                    log.info('  %-20s %-30s %.12f', ia, ib, c)
+
+        if with_meta_lowdin:
+            casscf._scf.mulliken_meta(casscf.mol, dm1, s=ovlp_ao, verbose=log)
+        else:
+            casscf._scf.mulliken_pop(casscf.mol, dm1, s=ovlp_ao, verbose=log)
+    return dm1a, dm1b
+    '''