EDA Functionality for LASSI-PDFT

examples/laspdft/H2molecules.py

@@ -3,36 +3,50 @@
 from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF
 from mrh.my_dmet import localintegrals, dmet, fragments
 from mrh.my_pyscf import mcpdft
+from mrh.my_pyscf import lassi
 
 # Initialization
 mol = gto.M()
 mol.atom='''H -5.10574 2.01997 0.00000; H -4.29369 2.08633 0.00000; H -3.10185 2.22603 0.00000; H -2.29672 2.35095 0.00000''' 
-mol.basis='sto3g'
+mol.basis='6-31g'
+mol.verbose=4
 mol.build()
 
 # Mean field calculation
 mf = scf.ROHF(mol).newton().run()
 
-'''
-# Option-1: Perform LASSCF and then pass las object
 # LASSCF Calculations
 las = LASSCF(mf,(2, 2),(2, 2),spin_sub=(1, 1))
 frag_atom_list = ([0, 1], [2, 3]) 
 mo0 = las.localize_init_guess(frag_atom_list)
-elas = las.kernel(mo0)[0]
-
-# LAS-PDFT
-mc = mcpdft.LASSCF(las, 'tPBE', (2, 2), (2,2), grids_level=1)
-epdft = mc.kernel()[0]
-'''
-# Option-2: Feed the mean field object and fragment information to mcpdft.LASSCF
-mc = mcpdft.LASSCF(mf, 'tPBE', (2, 2), (2, 2), spin_sub=(1,1), grids_level=1)
-frag_atom_list = ([0, 1] , [2, 3])
-mo0 = mc.localize_init_guess (frag_atom_list)
-mc.kernel(mo0)
-
-elas = mc.e_mcscf[0]
-epdft = mc.e_tot
-
-print ("E(LASSCF) =", elas)
-print ("E(tPBE) =", epdft)
+las.max_cycle_macro=1
+las.kernel(mo0)
+
+las = lassi.spaces.all_single_excitations (las)
+las.lasci ()
+
+lsi = lassi.LASSI(las)
+lsi.kernel()
+
+# LASSI-PDFT
+mc = mcpdft.LASSI(lsi, 'tPBE',states=[0, 1])
+mc.kernel()
+
+# Note: 
+# mc.e_tot : LASSI-PDFT energy
+# mc.e_mcscf: LASSI energy for those roots
+# mc.e_roots: LASSI energies of all the roots
+# mc.e_states: Energy of LAS states (basis of LASSI)
+
+# Analyze: Takes the state no as input. By default it is for state=0
+mc.analyze(state=[0,1])
+
+# Energy decomposition analysis: Will be done on the same no of roots calculated 
+# during the LASSI-PDFT
+e_decomp = mc.get_energy_decomposition (split_x_c=False)
+print ("e_nuc =",e_decomp[0])
+print ("e_1e =",e_decomp[1])
+print ("e_Coul =",e_decomp[2])
+print ("e_OT =",e_decomp[3])
+print ("e_ncwfn (not included in total energy) =",e_decomp[4])
+

my_pyscf/mcpdft/laspdft.py

@@ -4,6 +4,8 @@
 from mrh.my_pyscf.lassi import lassi
 import h5py
 import tempfile
+from pyscf.mcpdft.otfnal import transfnal, get_transfnal
+from pyscf.mcpdft.mcpdft import _get_e_decomp
 
 try:
     from pyscf.mcpdft.mcpdft import _PDFT, _mcscf_env
@@ -35,6 +37,65 @@ def make_casdm2s(filename, i):
     return rdm2s
 
 
+def get_energy_decomposition(mc, mo_coeff=None, ci=None, ot=None, otxc=None,
+                             grids_level=None, grids_attr=None,
+                             split_x_c=None, verbose=None):
+    """
+    I have to include this function to do the EDA for selected roots only.
+    """
+    if verbose is None: verbose = mc.verbose
+    log = lib.logger.new_logger(mc, verbose)
+    if mo_coeff is None: mo_coeff = mc.mo_coeff
+    if ci is None: ci = mc.ci
+    if grids_attr is None: grids_attr = {}
+    if grids_level is not None: grids_attr['level'] = grids_level
+    if len(grids_attr) or (otxc is not None):
+        old_ot = ot if (ot is not None) else mc.otfnal
+        old_grids = old_ot.grids
+        if otxc is None: otxc = old_ot.otxc
+        new_ot = get_transfnal(mc.mol, otxc)
+        new_ot.grids.__dict__.update(old_grids.__dict__)
+        new_ot.grids.__dict__.update(**grids_attr)
+        ot = new_ot
+    if ot is None: ot = mc.otfnal
+    if split_x_c is None:
+        split_x_c = True
+        log.warn('Currently, split_x_c in get_energy_decomposition defaults to '
+                 'True.\nThis default will change to False in the near future.')
+    hyb_x, hyb_c = ot._numint.hybrid_coeff(ot.otxc)
+    if hyb_x > 1e-10 or hyb_c > 1e-10:
+        raise NotImplementedError("Decomp for hybrid PDFT fnals")
+    if not isinstance(ot, transfnal):
+        raise NotImplementedError("Decomp for non-translated PDFT fnals")
+
+    if split_x_c:
+        ot = list(ot.split_x_c())
+    else:
+        ot = [ot, ]
+    e_nuc = mc._scf.energy_nuc()
+
+    nroots = getattr(mc.fcisolver, 'nroots', 1)
+    if not isinstance(nroots, list):
+        nroots = list(range(nroots))
+
+    if isinstance(nroots, list) and len(nroots) > 1:
+        e_1e, e_coul, e_otxc, e_ncwfn = [], [], [], []
+        for state in nroots:
+            row = _get_e_decomp(mc, mo_coeff=mo_coeff,ci=mc.ci,ot=ot,state=state)
+            e_1e.append(row[0])
+            e_coul.append(row[1])
+            e_otxc.append(row[2])
+            e_ncwfn.append(row[3])
+        e_otxc = [[e[i] for e in e_otxc] for i in range(len(e_otxc[0]))]
+    else:
+        e_1e, e_coul, e_otxc, e_ncwfn = _get_e_decomp(mc,mo_coeff=mo_coeff,ci=mc.ci,ot=ot,state=nroots[0])
+    if split_x_c:
+        e_otx, e_otc = e_otxc
+        return e_nuc, e_1e, e_coul, e_otx, e_otc, e_ncwfn
+    else:
+        return e_nuc, e_1e, e_coul, e_otxc[0], e_ncwfn
+
+
 class _LASPDFT(_PDFT):
     'MC-PDFT energy for a LASSCF wavefunction'
 
@@ -84,7 +145,7 @@ def compute_pdft_energy_(self, mo_coeff=None, ci=None, ot=None, otxc=None,
                                      logger_tag='MC-PDFT state {}'.format(ix))
                      for ix in range(nroots)]
 
-        self.e_ot = [e_ot for e_tot, e_ot in epdft]
+        self.e_ot = np.array([e_ot for e_tot, e_ot in epdft])
 
         if isinstance(self, StateAverageMCSCFSolver):
             e_states = [e_tot for e_tot, e_ot in epdft]
@@ -99,7 +160,7 @@ def compute_pdft_energy_(self, mo_coeff=None, ci=None, ot=None, otxc=None,
             self.e_tot = np.dot(e_states, self.weights)
             e_states = self.e_states
         elif (len(nroots) > 1 if isinstance(nroots, list) else nroots > 1):
-            self.e_tot = [e_tot for e_tot, e_ot in epdft]
+            self.e_tot = np.array([e_tot for e_tot, e_ot in epdft])
             e_states = self.e_tot
         else:  # nroots==1 not StateAverage class
             self.e_tot, self.e_ot = epdft[0]
@@ -136,6 +197,9 @@ def compute_pdft_energy_(self, mo_coeff=None, ci=None, ot=None, otxc=None,
             return _LASPDFT.compute_pdft_energy_(self, mo_coeff=mo_coeff, ci=ci, ot=ot, otxc=otxc,
                                                  grids_level=grids_level, grids_attr=grids_attr, **kwargs)
 
+        def get_energy_decomposition(self, **kwargs):
+            raise NotImplementedError('EDA is not yet defined for LAS-PDFT')
+
         def multi_state(self, **kwargs):
             """
             In future will have to change this to consider the modal space selection, weights...
@@ -149,13 +213,30 @@ def multi_state(self, **kwargs):
             _mc_class.DoLASSI = True
             _mc_class.rdmstmpfile = tempfile.NamedTemporaryFile(dir=lib.param.TMPDIR)
 
+            def analyze(self, state=0, **kwargs):
+                log = lib.logger.new_logger(self, self.verbose)
+                log.warn("Analyze function is not yet defined for LAS-PDFT. Turning on the analyze function of LASSI")
+                from mrh.my_pyscf.lassi.sitools import analyze
+                return analyze(self, self.si, state=state, **kwargs)
+
+            def get_energy_decomposition(self, mo_coeff=None, ci=None, ot=None, otxc=None, grids_level=None,
+                                         grids_attr=None, split_x_c=None, verbose=None):
+                return get_energy_decomposition(self, mo_coeff=mo_coeff, ci=ci, ot=ot, otxc=otxc,
+                                                grids_level=grids_level, grids_attr=grids_attr,
+                                                split_x_c=split_x_c, verbose=verbose)
         else:
             _mc_class.DoLASSI = False
             if mc.ci is not None:
                 mc.fcisolver.nroots = mc.nroots
             else:
                 mc.fcisolver.nroots = 1
 
+            def analyze(self):
+                raise NotImplementedError('Analyze function is not yet implemented for LAS-PDFT')
+
+            def get_energy_decomposition(self, **kwargs):
+                raise NotImplementedError('EDA is not yet implemented for LAS-PDFT')
+
         if states is not None: _mc_class.states = states
 
         if _mc_class.DoLASSI:
@@ -192,6 +273,11 @@ def _store_rdms(self):
 
                         rdm1s, rdm2s = lassi.root_make_rdm12s(self, self.ci, self.si,
                                                               state=self.states[i:j])
+
+                        if len(self.states[i:j]) == 1:
+                            rdm1s = [rdm1s]
+                            rdm2s = [rdm2s]
+
                         for k in range(i, j):
                             stateno = self.states[k]
                             rdm1s_dname = f'rdm1s_{stateno}'
@@ -221,6 +307,7 @@ def make_one_casdm2(self, ci=None, state=0, **kwargs):
         def optimize_mcscf_(self, mo_coeff=None, ci0=None, **kwargs):
             '''Optimize the MC-SCF wave function underlying an MC-PDFT calculation.
             Has the same calling signature as the parent kernel method. '''
+
             with _mcscf_env(self):
                 if self.DoLASSI:
                     self.statlis = [x for x in range(len(self.e_roots))]  # LASSI-LPDFT
@@ -229,14 +316,17 @@ def optimize_mcscf_(self, mo_coeff=None, ci0=None, **kwargs):
                         self.states = list(range(len(self.e_roots)))
                     else:
                         self.fcisolver.nroots = self.states
+
                     self._store_rdms()
                 else:
                     self.e_mcscf, self.e_cas, self.ci, self.mo_coeff, self.mo_energy = \
                         self._mc_class.kernel(self, mo_coeff, ci0=ci0, **kwargs)[:-2]
                     self.fcisolver.nroots = self.nroots
 
-    pdft = PDFT(mc._scf, mc.ncas_sub, mc.nelecas_sub, my_ot=ot, **kwargs)
+            if self.DoLASSI:
+                self.e_mcscf = self.e_roots[self.states]  # To be consistent with PySCF
 
+    pdft = PDFT(mc._scf, mc.ncas_sub, mc.nelecas_sub, my_ot=ot, **kwargs)
     _keys = pdft._keys.copy()
     pdft.__dict__.update(mc.__dict__)
     pdft._keys = pdft._keys.union(_keys)

tests/mcpdft/test_edalaspdft.py

@@ -0,0 +1,55 @@
+import unittest
+from pyscf import lib, gto, scf
+from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF
+from mrh.my_pyscf.lassi import LASSI
+from mrh.my_pyscf.lassi.spaces import all_single_excitations
+from mrh.my_pyscf.mcscf.lasci import get_space_info
+
+class KnownValues(unittest.TestCase):
+    def test_casci_limit (self):
+        xyz='''H 0 0 0
+               H 1 0 0
+               H 3 0 0
+               H 4 0 0'''
+
+        mol = gto.M (atom=xyz, basis='sto3g', symmetry=False, verbose=0, output='/dev/null')
+        mol.stdout.close()
+        mf = scf.RHF (mol).run ()
+        # LASSCF and LASSI
+        las = LASSCF (mf, (2,2), (2,2), spin_sub=(1,1))
+        las.lasci ()
+        las1 = las
+        for i in range (2): las1 = all_single_excitations (las1)
+        charges, spins, smults, wfnsyms = get_space_info (las1)
+        lroots = 4 - smults
+        idx = (charges!=0) & (lroots==3)
+        lroots[idx] = 1
+        las1.conv_tol_grad = las.conv_tol_self = 9e99
+        las1.lasci (lroots=lroots.T)
+        # CASCI limit
+        from pyscf import mcpdft
+        mc = mcpdft.CASCI (mf, 'tPBE', 4, 4).run ()
+        ref_e_decomp1 = mc.get_energy_decomposition (split_x_c=False)
+        ref_e_decomp2 = mc.get_energy_decomposition (split_x_c=True)
+        for opt in range (2):
+            with self.subTest (opt=opt):
+                lsi = LASSI (las1)
+                lsi.opt = opt
+                lsi.kernel ()
+                self.assertAlmostEqual (lsi.e_roots[0], mc.e_mcscf, 7)
+                from mrh.my_pyscf import mcpdft
+                lsipdft = mcpdft.LASSI (lsi, 'tPBE')
+                lsipdft.opt = opt
+                lsipdft.kernel()
+                e_decomp_ = lsipdft.get_energy_decomposition (split_x_c=False)
+                e_decomp1 = [e_decomp_[0], *[e[0] for e in e_decomp_[1:5]]]
+                e_decomp_ = lsipdft.get_energy_decomposition (split_x_c=True)
+                e_decomp2 = [e_decomp_[0], *[e[0] for e in e_decomp_[1:6]]]
+                [self.assertAlmostEqual(e1, ref_e1, 7) for e1, ref_e1 in zip(e_decomp1, ref_e_decomp1)]
+                [self.assertAlmostEqual(e2, ref_e2, 7) for e2, ref_e2 in zip(e_decomp2, ref_e_decomp2)]
+
+if __name__ == "__main__":
+    print("Full Tests for Energy Decomposition of LASSI-PDFT")
+    unittest.main()
+
+