Merge pull request #121 from JangidBhavnesh/dev

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
@@ -59,53 +61,6 @@ def get_h2eff(self, mo_coeff=None):
                              max_memory=self.max_memory)
         return eri
 
-    def compute_pdft_energy_(self, mo_coeff=None, ci=None, ot=None, otxc=None,
-                             grids_level=None, grids_attr=None, **kwargs):
-        '''Compute the MC-PDFT energy(ies) (and update stored data)
-        with the MC-SCF wave function fixed. '''
-        '''
-        Instead of finding the energies of all the states, this can allow
-        to take state number for which you want to add the PDFT corrections
-        '''
-        if mo_coeff is not None: self.mo_coeff = mo_coeff
-        if ci is not None: self.ci = ci
-        if ot is not None: self.otfnal = ot
-        if otxc is not None: self.otxc = otxc
-        if grids_attr is None: grids_attr = {}
-        if grids_level is not None: grids_attr['level'] = grids_level
-        if len(grids_attr): self.grids.__dict__.update(**grids_attr)
-        nroots = getattr(self.fcisolver, 'nroots', 1)
-        if isinstance(nroots, list):
-            epdft = [self.energy_tot(mo_coeff=self.mo_coeff, ci=self.ci, state=ix,
-                                     logger_tag='MC-PDFT state {}'.format(ix))
-                     for ix in nroots]
-        else:
-            epdft = [self.energy_tot(mo_coeff=self.mo_coeff, ci=self.ci, state=ix,
-                                     logger_tag='MC-PDFT state {}'.format(ix))
-                     for ix in range(nroots)]
-
-        self.e_ot = [e_ot for e_tot, e_ot in epdft]
-
-        if isinstance(self, StateAverageMCSCFSolver):
-            e_states = [e_tot for e_tot, e_ot in epdft]
-            try:
-                self.e_states = e_states
-            except AttributeError as e:
-                self.fcisolver.e_states = e_states
-                assert (self.e_states is e_states), str(e)
-            # TODO: redesign this. MC-SCF e_states is stapled to
-            # fcisolver.e_states, but I don't want MS-PDFT to be
-            # because that makes no sense
-            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]
-            e_states = self.e_tot
-        else:  # nroots==1 not StateAverage class
-            self.e_tot, self.e_ot = epdft[0]
-            e_states = [self.e_tot]
-        return self.e_tot, self.e_ot, e_states
-
     def multi_state(self, method='Lin'):
         if method.upper() == "LIN":
             from mrh.my_pyscf.mcpdft._lpdft import linear_multi_state
@@ -130,11 +85,12 @@ def get_h2eff(self, mo_coeff=None):
                 return mc.__class__.get_h2eff(self, mo_coeff=mo_coeff)
             else:
                 return _LASPDFT.get_h2eff(self, mo_coeff=mo_coeff)
-
+
+        # Have to pass this due to dump_chk, which won't work for LAS.
         def compute_pdft_energy_(self, mo_coeff=None, ci=None, ot=None, otxc=None,
-                                 grids_level=None, grids_attr=None, states=states, **kwargs):
+                                 grids_level=None, grids_attr=None, dunp_chk=False, **kwargs):
             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)
+                    grids_level=grids_level, grids_attr=grids_attr, dump_chk=False, **kwargs)
 
         def multi_state(self, **kwargs):
             """
@@ -149,13 +105,25 @@ 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)
+
         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 +160,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}'
@@ -203,11 +176,11 @@ def _store_rdms(self):
 
             def make_one_casdm1s(self, ci=None, state=0, **kwargs):
                 rdmstmpfile = self.rdmstmpfile
-                return make_casdm1s(rdmstmpfile, state)
+                return make_casdm1s(rdmstmpfile, self.states[state])
 
             def make_one_casdm2(self, ci=None, state=0, **kwargs):
                 rdmstmpfile = self.rdmstmpfile
-                return make_casdm2s(rdmstmpfile, state).sum((0, 3))
+                return make_casdm2s(rdmstmpfile, self.states[state]).sum((0, 3))
 
         else:
             make_one_casdm1s = mc.__class__.state_make_casdm1s
@@ -221,22 +194,26 @@ 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
                     if self.states is None:
-                        self.fcisolver.nroots = len(self.e_roots)
                         self.states = list(range(len(self.e_roots)))
+                        self.fcisolver.nroots = len(self.e_roots)
                     else:
-                        self.fcisolver.nroots = self.states
+                        self.fcisolver.nroots = len(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()
+
+