LASSI-PDFT updated the framework & example

examples/laspdft/c2h4n4_si_laspdft.py

@@ -1,7 +1,7 @@
 from pyscf import gto, scf, lib, mcscf
 from mrh.my_pyscf.fci import csf_solver
 from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF
-from mrh.my_pyscf.mcscf import lassi
+from mrh.my_pyscf import lassi
 from mrh.my_pyscf import mcpdft
 from mrh.my_pyscf.tools.molden import from_lasscf
 from c2h4n4_struct import structure as struct
@@ -22,16 +22,27 @@
 guess_mo = las.sort_mo([16,18,22,23,24,26])
 mo0 = las.localize_init_guess((list(range (5)), list(range (5,10))), guess_mo)
 las.kernel(mo0)
-
+
+las = lassi.states.all_single_excitations (las)
+las.lasci ()
+
+lsi = lassi.LASSI(las)
+lsi.kernel()
+
 # LASSI-PDFT
-mc = mcpdft.LASSI(las, 'tPBE', (3, 3), ((2,1),(1,2)))
-mc = all_single_excitations(mc) # Level of charge transfer
-mc.kernel(las.mo_coeff) # SA-LAS orbitals
+mc = mcpdft.LASSI(lsi, 'tPBE', (3, 3), ((2,1),(1,2)))
+mc.kernel() 
+
+# CASCI-PDFT in las orbitals
+from pyscf import mcpdft
+mc_ci = mcpdft.CASCI(mf, 'tPBE', 6, 6)
+mc_ci.kernel(las.mo_coeff)
 
 # Results
 print("\n----Results-------\n")
-#print("State",' \t',  "LASSCF Energy",'\t\t',"LASSI Energy",'\t\t', "LASSI-PDFT Energy") 
-#[print(sn,'\t',x,'\t', y,'\t', z) for sn, x, y, z in zip(list(range(mc.nroots)), mc.e_mcscf, mc.e_lassi, mc.e_tot)]
-print("LASSI state-0 =", mc.e_mcscf[0])
+print("CASCI state-0 =", mc_ci.e_mcscf)
+print("LASSI state-0 =", lsi.e_roots[0])
+
+print("CASCI-PDFT state-0 =", mc_ci.e_tot)
 print("LASSI-PDFT state-0 =", mc.e_tot[0])
 

my_pyscf/mcpdft/__init__.py

@@ -1,5 +1,6 @@
 # Lahh dee dah
 import copy
+import numpy as np
 from pyscf.mcpdft.mcpdft import get_mcpdft_child_class
 from pyscf import mcscf, gto
 from pyscf.lib import logger
@@ -48,7 +49,6 @@ def CASCIPDFT (mc_or_mf_or_mol, ot, ncas, nelecas, ncore=None, **kwargs):
 CASSCF=CASSCFPDFT
 CASCI=CASCIPDFT
 
-
 # LAS-PDFT
 def _laspdftEnergy(mc_class, mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, DoLASSI=False, ncore=None, spin_sub=None,
                    frozen=None, **kwargs):
@@ -64,16 +64,43 @@ def _laspdftEnergy(mc_class, mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, DoLASSI
 
     if frozen is not None: mc1 = mc_class(mf_or_mol, ncas_sub, nelecas_sub, ncore=ncore, spin_sub = spin_sub, frozen=frozen)
     else:  mc1 = mc_class(mf_or_mol, ncas_sub, nelecas_sub, ncore=ncore, spin_sub = spin_sub)
-    
+
     from mrh.my_pyscf.mcpdft.laspdft import get_mcpdft_child_class
     mc2 = get_mcpdft_child_class(mc1, ot, DoLASSI=DoLASSI, **kwargs)
 
     if mc0 is not None:
-        mc2.mo_coeff = mc_or_mf_or_mol.mo_coeff.copy ()    
+        mc2.mo_coeff = mc_or_mf_or_mol.mo_coeff.copy ()
         mc2.ci = copy.deepcopy (mc_or_mf_or_mol.ci)
         mc2.converged = mc0.converged
     return mc2
-
+
+
+def _lassipdftEnergy(mc_class, mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, DoLASSI=False, ncore=None, spin_sub=None,
+                   frozen=None, **kwargs):
+
+    from mrh.my_pyscf.lassi import lassi, lassis
+
+    if isinstance(mc_or_mf_or_mol, (lassi.LASSI, lassis.LASSIS)):
+        mc0 = mc_or_mf_or_mol._las
+        mf_or_mol = mc_or_mf_or_mol._las._scf
+    else:
+        raise "Requires lassi instance"
+
+    mc1 = mc_class(mf_or_mol, ncas_sub, nelecas_sub, ncore=ncore, spin_sub=spin_sub)
+
+    from mrh.my_pyscf.mcpdft.laspdft import get_mcpdft_child_class
+    mc2 = get_mcpdft_child_class(mc1, ot, DoLASSI=DoLASSI, **kwargs)
+
+    if mc0 is not None:
+        mc2.mo_coeff = mc_or_mf_or_mol.mo_coeff.copy()
+        mc2.ci = copy.deepcopy(mc_or_mf_or_mol.ci)
+        mc2.converged = mc0.converged
+        _keys = mc_or_mf_or_mol._keys.copy()
+        mc2.__dict__.update(mc_or_mf_or_mol.__dict__)
+        mc2._keys = mc2._keys.union(_keys)
+        mc2.e_mcscf = np.average(mc_or_mf_or_mol.e_roots).copy()
+    return mc2
+
 def LASSCFPDFT(mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub,  ncore=None, spin_sub=None, frozen=None,
                **kwargs):
     from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF
@@ -83,11 +110,13 @@ def LASSCFPDFT(mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub,  ncore=None, spin_sub
 def LASSIPDFT(mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, ncore=None, spin_sub=None, frozen=None,
                **kwargs):
     from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF
-    return _laspdftEnergy(LASSCF,  mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, DoLASSI=True, ncore=ncore,
+    return _lassipdftEnergy(LASSCF,  mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, DoLASSI=True, ncore=ncore,
                           spin_sub=spin_sub, frozen=frozen, **kwargs)
 
+
 LASSCF = LASSCFPDFT
 LASSI = LASSIPDFT
+LASSIS = LASSIPDFT
 
 def CIMCPDFT (*args, **kwargs):
     from mrh.my_pyscf.mcpdft.var_mcpdft import CIMCPDFT as fn

my_pyscf/mcpdft/laspdft.py

@@ -68,13 +68,13 @@ 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):
-                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
                 if self.DoLASSI:
-                    self.e_states, self.si = self.lassi()
-                return self.e_mcscf, self.e_cas, self.ci, self.mo_coeff, self.mo_energy
-
+                    self.fcisolver.nroots = len(self.e_states)
+                    self.e_states = self.e_roots
+                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)
 
     _keys = pdft._keys.copy()

tests/mcpdft/test_lassipdft.py

     mc = mcscf.CASCI (mf, 4, 4).run () 
     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') 
     mf = scf.RHF (mol).run () 
     # Random Hamiltonian 
     rng = np.random.default_rng (424) 
     mf._eri = rng.random (mf._eri.shape) 
     hcore = rng.random ((4,4)) 
     hcore = hcore + hcore.T 
     mf.get_hcore = lambda *args: hcore 
     # LASSCF with CASCI-limit model space 
     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) 
     las1.dump_spaces () 
     lsi = LASSI (las1) 
     lsi.kernel (opt=0) 
     # CASCI limit 
     mc = mcscf.CASCI (mf, 4, 4).run () 
@@ -3,6 +3,7 @@
 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
 import unittest
 
 class KnownValues(unittest.TestCase):
@@ -19,10 +20,15 @@ def test_ethene (self):
         mol.output = '/dev/null'
         mol.build()
         mf = scf.ROHF(mol).newton().run()
-        mc = mcpdft.LASSI(mf, 'tPBE', (2, ), (2, ), grid_level=1).state_average(
-            [1, 0], spins=[[0,], [2, ]], smults=[[1, ], [3, ]], charges=[[0, ],[0, ]])
-        mo0 = mc.localize_init_guess(([0, 1],), mc.sort_mo([8, 9]))
-        mc.kernel(mo0)
+
+        las = LASSCF(mf, (2,),(2,))
+        las = las.state_average([1, 0], spins=[[0,], [2, ]], smults=[[1, ], [3, ]], charges=[[0, ],[0, ]])
+        mo0 = las.localize_init_guess(([0, 1],), las.sort_mo([8, 9]))
+        las.kernel(mo0)
+        lsi = lassi.LASSI(las)
+        lsi.kernel()
+        mc = mcpdft.LASSI(lsi, 'tPBE', (2, ), (2, ))
+        mc.kernel()
         elassi = mc.e_mcscf[0]
         epdft = mc.e_tot[0]
         self.assertAlmostEqual (elassi , -77.1154672717181, 7) # Reference values of CASSCF and CAS-PDFT
@@ -32,5 +38,4 @@ def test_ethene (self):
     print("Full Tests for LASSI-PDFT")
     unittest.main()
 
-