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
+
+    if isinstance(mc_or_mf_or_mol, lassi.LASSI):
+        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 # Not Sure acc. to issue #34 of mrh
 
 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

@@ -1,36 +1,48 @@
-import sys
-from pyscf import gto, scf, tools, dft, lib
-from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF
-from mrh.my_dmet import localintegrals, dmet, fragments
-from mrh.my_pyscf import mcpdft
 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.states import all_single_excitations
+from mrh.my_pyscf.mcscf.lasci import get_space_info
 
 class KnownValues(unittest.TestCase):
     def test_ethene (self):
-        mol = gto.M()
-        mol.atom = '''C  -0.662958  0.000000  0.000000;
-                            C   0.662958  0.000000  0.000000;
-                            H  -1.256559  -0.924026  0.000000;
-                            H  1.256559  -0.924026  0.000000;
-                            H  -1.256559  0.924026  0.000000;
-                            H  1.256559  0.924026  0.000000'''
-        mol.basis='sto3g'
-        mol.verbose = 0
-        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)
-        elassi = mc.e_mcscf[0]
-        epdft = mc.e_tot[0]
-        self.assertAlmostEqual (elassi , -77.1154672717181, 7) # Reference values of CASSCF and CAS-PDFT
-        self.assertAlmostEqual (epdft , -77.49805221093968, 7)
+        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 ()
+
+        # 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)
+        las1.dump_spaces ()
+        lsi = LASSI (las1)
+        lsi.kernel (opt=0)
+
+        from mrh.my_pyscf import mcpdft
+        lsipdft = mcpdft.LASSI(lsi, 'tPBE', (2,2), (2,2))
+        lsipdft.kernel()
+
+        # CASCI limit
+        from pyscf import mcpdft
+        mc = mcpdft.CASCI (mf, 'tPBE', 4, 4).run ()
+
+        self.assertAlmostEqual (lsi.e_roots[0], mc.e_mcscf, 7)
+        self.assertAlmostEqual (lsipdft.e_tot[0], mc.e_tot, 7)
 
 if __name__ == "__main__":
     print("Full Tests for LASSI-PDFT")
     unittest.main()
 
-