diff --git a/examples/laspdft/lassi_lpdft.py b/examples/laspdft/lassi_lpdft.py new file mode 100644 index 00000000..11438d69 --- /dev/null +++ b/examples/laspdft/lassi_lpdft.py @@ -0,0 +1,62 @@ +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 import lassi +from mrh.my_pyscf import mcpdft +from c2h4n4_struct import structure as struct +from mrh.my_pyscf.lassi.spaces import all_single_excitations + +# Mean field calculation +mol = struct(0, 0, '6-31g') +mol.verbose = lib.logger.INFO +mol.build() +mf = scf.RHF(mol).run() + +# SA-LASSCF: For orbitals +las = LASSCF(mf, (3,3), ((2,1),(1,2))) +las = las.state_average([0.5,0.5], spins=[[1,-1],[-1,1]], smults=[[2,2],[2,2]], charges=[[0,0],[0,0]],wfnsyms=[[1,1],[1,1]]) +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.spaces.all_single_excitations (las) +las.lasci () + +lsi = lassi.LASSI(las) +lsi.kernel() + +# LASSI-PDFT +mc = mcpdft.LASSI(lsi, 'tPBE') +mc = mc.multi_state() +mc.kernel() + +# CASCI-PDFT in las orbitals +from pyscf import mcpdft +mc_sci = mcpdft.CASCI(mf, 'tPBE', 6, 6) +mc_sci.fcisolver = csf_solver(mol, smult=1) +mc_sci.kernel(las.mo_coeff) + +mc_tci = mcpdft.CASCI(mf, 'tPBE', 6, (4, 2)) +mc_tci.fcisolver = csf_solver(mol, smult=3) +mc_tci.kernel(las.mo_coeff) + + +# CASSCF-PDFT in las orbitals +from pyscf import mcpdft +mcas_sci = mcpdft.CASSCF(mf, 'tPBE', 6, 6) +mcas_sci.fcisolver = csf_solver(mol, smult=1) +mcas_sci.kernel(las.mo_coeff) + +mcas_tci = mcpdft.CASSCF(mf, 'tPBE', 6, (4, 2)) +mcas_tci.fcisolver = csf_solver(mol, smult=3) +mcas_tci.kernel(las.mo_coeff) + +# Results Singlet-Triplet Gap +print("\n----Results-------\n") +print("CASSCF S-T =", 27.21139*(mcas_tci.e_mcscf - mcas_sci.e_mcscf)) +print("CASCI S-T =", 27.21139*(mc_tci.e_mcscf - mc_sci.e_mcscf)) +print("LASSI S-T =", 27.21139*(lsi.e_roots[1]-lsi.e_roots[0])) +print("CASSCF-LPDFT S-T =", 27.21139*(mcas_tci.e_tot - mcas_sci.e_tot)) +print("CASCI-LPDFT S-T =", 27.21139*(mc_tci.e_tot - mc_sci.e_tot)) +print("LASSI-LPDFT S-T =", 27.21139*(mc.e_states[1]-mc.e_states[0])) + diff --git a/my_pyscf/mcpdft/_lpdft.py b/my_pyscf/mcpdft/_lpdft.py new file mode 100644 index 00000000..2218da69 --- /dev/null +++ b/my_pyscf/mcpdft/_lpdft.py @@ -0,0 +1,304 @@ +import numpy as np +from scipy import linalg +from pyscf.lib import logger +from pyscf import mcpdft, lib +from pyscf.mcpdft import _dms +from mrh.my_pyscf.lassi.lassi import las_symm_tuple, iterate_subspace_blocks +from mrh.my_pyscf.lassi import op_o1 +from pyscf.mcpdft import lpdft as lpdft_fns +''' +This file is taken from pyscf-forge and adopted for the LAS wavefunctions. +''' + + +def weighted_average_densities(mc): + """ + Compute the weighted average 1- and 2-electron LAS densities + in the selected modal space + """ + casdm1s = [mc.make_one_casdm1s(mc.ci, state=state) for state in mc.statlis] + casdm2 = [mc.make_one_casdm2(mc.ci, state=state) for state in mc.statlis] + weights = [1 / len(mc.statlis), ] * len(mc.statlis) + return (np.tensordot(weights, casdm1s, axes=1)), (np.tensordot(weights, casdm2, axes=1)) + + +# Importing functions from the PySCF-forge +get_lpdft_hconst = lpdft_fns.get_lpdft_hconst +transformed_h1e_for_cas = lpdft_fns.transformed_h1e_for_cas +get_transformed_h2eff_for_cas = lpdft_fns.get_transformed_h2eff_for_cas + +def make_lpdft_ham_(mc, mo_coeff=None, ci=None, ot=None): + '''Compute the L-PDFT Hamiltonian + + Args: + mo_coeff : ndarray of shape (nao, nmo) + A full set of molecular orbital coefficients. Taken from self if + not provided. + + ci : list of ndarrays of length nroots + CI vectors should be from a converged CASSCF/CASCI calculation + + ot : an instance of on-top functional class - see otfnal.py + + Returns: + lpdft_ham : ndarray of shape (nroots, nroots) or (nirreps, nroots, nroots) + Linear approximation to the MC-PDFT energy expressed as a + hamiltonian in the basis provided by the CI vectors. If + StateAverageMix, then returns the block diagonal of the lpdft + hamiltonian for each irrep. + ''' + + if mo_coeff is None: mo_coeff = mc.mo_coeff + if ci is None: ci = mc.ci + if ot is None: ot = mc.otfnal + + ot.reset(mol=mc.mol) + + spin = abs(mc.nelecas[0] - mc.nelecas[1]) + omega, _, hyb = ot._numint.rsh_and_hybrid_coeff(ot.otxc, spin=spin) + if abs(omega) > 1e-11: + raise NotImplementedError("range-separated on-top functionals") + if abs(hyb[0] - hyb[1]) > 1e-11: + raise NotImplementedError( + "hybrid functionals with different exchange, correlations components") + + cas_hyb = hyb[0] + + ncas = mc.ncas + casdm1s_0, casdm2_0 = mc.get_casdm12_0() + + mc.veff1, mc.veff2, E_ot = mc.get_pdft_veff(mo=mo_coeff, casdm1s=casdm1s_0, + casdm2=casdm2_0, drop_mcwfn=True, incl_energy=True) + + # This is all standard procedure for generating the hamiltonian in PySCF + h1, h0 = mc.get_h1lpdft(E_ot, casdm1s_0, casdm2_0, hyb=1.0 - cas_hyb) + h2 = mc.get_h2lpdft() + + statesym, s2_states = las_symm_tuple(mc, verbose=0) + + # Initialize matrices + e_roots = [] + s2_roots = [] + rootsym = [] + si = [] + s2_mat = [] + idx_allprods = [] + # Loop over symmetry blocks + qn_lbls = ['neleca', 'nelecb', 'irrep'] + for it, (las1, sym, indices, indexed) in enumerate(iterate_subspace_blocks(mc, ci, statesym)): + idx_space, idx_prod = indices + ci_blk, nelec_blk = indexed + idx_allprods.extend(list(np.where(idx_prod)[0])) + lib.logger.info(mc, 'Build + diag H matrix L-PDFT-LASSI symmetry block %d\n' + + '{} = {}\n'.format(qn_lbls, sym) + + '(%d rootspaces; %d states)', it, + np.count_nonzero(idx_space), + np.count_nonzero(idx_prod)) + + ham_blk, s2_blk, ovlp_blk = op_o1.ham(mc, h1, h2, ci_blk, nelec_blk) + diag_idx = np.diag_indices_from(ham_blk) + ham_blk[diag_idx] += h0 + cas_hyb * mc.e_roots + + try: + e, c = linalg.eigh(ham_blk, b=ovlp_blk) + except linalg.LinAlgError as err: + ovlp_det = linalg.det(ovlp_blk) + lc = 'checking if L-PDFT-LASSI basis has lindeps: |ovlp| = {:.6e}'.format(ovlp_det) + lib.logger.info(las, 'Caught error %s, %s', str(err), lc) + if ovlp_det < LINDEP_THRESH: + x = canonical_orth_(ovlp_blk, thr=LINDEP_THRESH) + lib.logger.info(las, '%d/%d linearly independent model states', + x.shape[1], x.shape[0]) + xhx = x.conj().T @ ham_blk @ x + e, c = linalg.eigh(xhx) + c = x @ c + else: + raise (err) from None + + s2_mat.append(s2_blk) + si.append(c) + s2_blk = c.conj().T @ s2_blk @ c + lib.logger.debug2(mc, 'Block S**2 in adiabat basis:') + lib.logger.debug2(mc, '{}'.format(s2_blk)) + e_roots.extend(list(e)) + s2_roots.extend(list(np.diag(s2_blk))) + rootsym.extend([sym, ] * c.shape[1]) + + idx_allprods = np.argsort(idx_allprods) + si = linalg.block_diag(*si)[idx_allprods, :] + s2_mat = linalg.block_diag(*s2_mat)[np.ix_(idx_allprods, idx_allprods)] + idx = np.argsort(e_roots) + rootsym = np.asarray(rootsym)[idx] + s2_roots = np.asarray(s2_roots)[idx] + si = si[:, idx] + return ham_blk, e, si, rootsym, s2_roots, s2_mat + + +def kernel(mc, mo_coeff=None, ot=None, **kwargs): + if ot is None: ot = mc.otfnal + if mo_coeff is None: mo_coeff = mc.mo_coeff + mc.optimize_mcscf_(mo_coeff=mo_coeff, **kwargs) + mc.lpdft_ham, mc.e_states, mc.si_pdft, mc.rootsym, mc.s2_roots, s2_mat = mc.make_lpdft_ham_(ot=ot) + logger.debug(mc, f"L-PDFT Hamiltonian in LASSI Basis:\n{mc.get_lpdft_ham()}") + + logger.debug(mc, f"L-PDFT SI:\n{mc.si_pdft}") + + mc.e_tot = mc.e_states[0] + logger.info(mc, 'LASSI-LPDFT eigenvalues (%d total):', len(mc.e_states)) + + fmt_str = '{:2s} {:>16s} {:2s} ' + col_lbls = ['ix', 'Energy', ''] + logger.info(mc, fmt_str.format(*col_lbls)) + fmt_str = '{:2d} {:16.10f} {:6.3f} ' + for ix, (er, s2r) in enumerate(zip(mc.e_states, mc.s2_roots)): + row = [ix, er, s2r] + logger.info(mc, fmt_str.format(*row)) + if ix >= 99 and mc.verbose < lib.logger.DEBUG: + lib.logger.info(mc, ('Remaining %d eigenvalues truncated; ' + 'increase verbosity to print them all'), len(mc.e_states) - 100) + break + + return mc.e_tot, mc.e_states, mc.si_pdft, mc.s2_roots + + +class _LPDFT(mcpdft.MultiStateMCPDFTSolver): + '''Linerized PDFT + + Saved Results + + e_tot : float + Weighted-average L-PDFT final energy + e_states : ndarray of shape (nroots) + L-PDFT final energies of the adiabatic states + ci : list of length (nroots) of ndarrays + CI vectors in the optimized adiabatic basis of MC-SCF. Related to + the L-PDFT adiabat CI vectors by the expansion coefficients + ``si_pdft''. + si_pdft : ndarray of shape (nroots, nroots) + Expansion coefficients of the L-PDFT adiabats in terms of the + optimized + MC-SCF adiabats + e_mcscf : ndarray of shape (nroots) + Energies of the MC-SCF adiabatic states + lpdft_ham : ndarray of shape (nroots, nroots) + L-PDFT Hamiltonian in the MC-SCF adiabatic basis + veff1 : ndarray of shape (nao, nao) + 1-body effective potential in the AO basis computed using the + zeroth-order densities. + veff2 : pyscf.mcscf.mc_ao2mo._ERIS instance + Relevant 2-body effective potential in the MO basis. + ''' + + def __init__(self, mc): + self.__dict__.update(mc.__dict__) + keys = set(('lpdft_ham', 'si_pdft', 'veff1', 'veff2')) + self.lpdft_ham = None + self.si_pdft = None + self.veff1 = None + self.veff2 = None + self._e_states = None + self._keys = set(self.__dict__.keys()).union(keys) + + make_lpdft_ham_ = make_lpdft_ham_ + make_lpdft_ham_.__doc__ = make_lpdft_ham_.__doc__ + + get_lpdft_hconst = get_lpdft_hconst + get_lpdft_hconst.__doc__ = get_lpdft_hconst.__doc__ + + get_h1lpdft = transformed_h1e_for_cas + get_h1lpdft.__doc__ = transformed_h1e_for_cas.__doc__ + + get_h2lpdft = get_transformed_h2eff_for_cas + get_h2lpdft.__doc__ = get_transformed_h2eff_for_cas.__doc__ + + get_casdm12_0 = weighted_average_densities + get_casdm12_0.__doc__ = weighted_average_densities.__doc__ + + def get_lpdft_ham(self): + '''The L-PDFT effective Hamiltonian matrix + + Returns: + lpdft_ham : ndarray of shape (nroots, nroots) + Contains L-PDFT Hamiltonian elements on the off-diagonals + and PDFT approx energies on the diagonals + ''' + return self.lpdft_ham + + def kernel(self, mo_coeff=None, ci0=None, ot=None, verbose=None): + ''' + Returns: + 6 elements, they are + total energy, + the MCSCF energies, + the active space CI energy, + the active space FCI wave function coefficients, + the MCSCF canonical orbital coefficients, + the MCSCF canonical orbital energies + + They are attributes of the QLPDFT object, which can be accessed by + .e_tot, .e_mcscf, .e_cas, .ci, .mo_coeff, .mo_energy + ''' + + if ot is None: ot = self.otfnal + ot.reset(mol=self.mol) # scanner mode safety + + if mo_coeff is None: + mo_coeff = self.mo_coeff + else: + self.mo_coeff = mo_coeff + + log = logger.new_logger(self, verbose) + + kernel(self, mo_coeff, ot=ot, verbose=log) + + return ( + self.e_tot, self.e_mcscf, self.e_cas, self.ci, + self.mo_coeff, self.mo_energy) + + def get_lpdft_hcore_only(self, casdm1s_0, hyb=1.0): + ''' + Returns the lpdft hcore AO integrals weighted by the + hybridization factor. Excludes the MC-SCF (wfn) component. + ''' + + dm1s = _dms.casdm1s_to_dm1s(self, casdm1s=casdm1s_0) + dm1 = dm1s[0] + dm1s[1] + v_j = self._scf.get_j(dm=dm1) + return hyb * self.get_hcore() + self.veff1 + hyb * v_j + + def get_lpdft_hcore(self, casdm1s_0=None): + ''' + Returns the full lpdft hcore AO integrals. Includes the MC-SCF + (wfn) component for hybrid functionals. + ''' + if casdm1s_0 is None: + casdm1s_0 = self.get_casdm12_0()[0] + + spin = abs(self.nelecas[0] - self.nelecas[1]) + cas_hyb = self.otfnal._numint.rsh_and_hybrid_coeff(self.otfnal.otxc, spin=spin)[2] + hyb = 1.0 - cas_hyb[0] + + return cas_hyb[0] * self.get_hcore() + self.get_lpdft_hcore_only(casdm1s_0, hyb=hyb) + + +def linear_multi_state(mc, **kwargs): + ''' Build linearized multi-state MC-PDFT method object + + Args: + mc : instance of class _PDFT + + Kwargs: + weights : sequence of floats + + Returns: + si : instance of class _LPDFT + ''' + + base_name = mc.__class__.__name__ + mcbase_class = mc.__class__ + + class LPDFT(_LPDFT, mcbase_class): + pass + + LPDFT.__name__ = "LIN" + base_name + return LPDFT(mc) diff --git a/my_pyscf/mcpdft/laspdft.py b/my_pyscf/mcpdft/laspdft.py index 568648ea..061a2707 100644 --- a/my_pyscf/mcpdft/laspdft.py +++ b/my_pyscf/mcpdft/laspdft.py @@ -13,9 +13,10 @@ try: from pyscf.mcpdft.mcpdft import _PDFT, _mcscf_env except ImportError: - msg = "For performing LASPDFT, you will require pyscf-forge.\n" +\ - "pyscf-forge can be found at : https://github.com/pyscf/pyscf-forge" - raise ImportError(msg) + msg = "For performing LASPDFT, you will require pyscf-forge.\n" + \ + "pyscf-forge can be found at : https://github.com/pyscf/pyscf-forge" + raise ImportError(msg) + def make_casdm1s(filename, i): ''' @@ -27,6 +28,7 @@ def make_casdm1s(filename, i): rdm1s = np.array(rdm1s) return rdm1s + def make_casdm2s(filename, i): ''' This function stores the rdm2s for the given state 'i' in a tempfile @@ -37,23 +39,26 @@ def make_casdm2s(filename, i): rdm2s = np.array(rdm2s) return rdm2s + class _LASPDFT(_PDFT): 'MC-PDFT energy for a LASSCF wavefunction' - + def get_h2eff(self, mo_coeff=None): 'Compute the active space two-particle Hamiltonian.' ncore = self.ncore ncas = self.ncas nocc = ncore + ncas - if mo_coeff is None: mo_coeff = self.mo_coeff[:,ncore:nocc] - elif mo_coeff.shape[1] != ncas: mo_coeff = mo_coeff[:,ncore:nocc] + if mo_coeff is None: + mo_coeff = self.mo_coeff[:, ncore:nocc] + elif mo_coeff.shape[1] != ncas: + mo_coeff = mo_coeff[:, ncore:nocc] - if getattr (self._scf, '_eri', None) is not None: + if getattr(self._scf, '_eri', None) is not None: eri = ao2mo.full(self._scf._eri, mo_coeff, - max_memory=self.max_memory) + max_memory=self.max_memory) else: eri = ao2mo.full(self.mol, mo_coeff, verbose=self.verbose, - max_memory=self.max_memory) + max_memory=self.max_memory) return eri def compute_pdft_energy_(self, mo_coeff=None, ci=None, ot=None, otxc=None, @@ -74,15 +79,15 @@ def compute_pdft_energy_(self, mo_coeff=None, ci=None, ot=None, otxc=None, 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] + 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)] + 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: @@ -103,41 +108,64 @@ def compute_pdft_energy_(self, mo_coeff=None, ci=None, ot=None, otxc=None, e_states = [self.e_tot] return self.e_tot, self.e_ot, e_states -def get_mcpdft_child_class(mc, ot, DoLASSI=False,states=None,**kwargs): + def multi_state(self, method='Lin'): + if method.upper() == "LIN": + from mrh.my_pyscf.mcpdft._lpdft import linear_multi_state + return linear_multi_state(self) + else: + raise NotImplementedError(f"StateAverageMix not available for {method}") + + +def get_mcpdft_child_class(mc, ot, DoLASSI=False, states=None, **kwargs): mc_doc = (mc.__class__.__doc__ or 'No docstring for MC-SCF parent method') - + class PDFT(_LASPDFT, mc.__class__): - __doc__= mc_doc + '\n\n' + _LASPDFT.__doc__ + __doc__ = mc_doc + '\n\n' + _LASPDFT.__doc__ _mc_class = mc.__class__ setattr(_mc_class, 'DoLASSI', None) setattr(_mc_class, 'states', None) + setattr(_mc_class, 'statlis', None) setattr(_mc_class, 'rdmstmpfile', None) - + def get_h2eff(self, mo_coeff=None): - if self._in_mcscf_env: return mc.__class__.get_h2eff(self, mo_coeff=mo_coeff) - else: return _LASPDFT.get_h2eff(self, mo_coeff=mo_coeff) - + if self._in_mcscf_env: + return mc.__class__.get_h2eff(self, mo_coeff=mo_coeff) + else: + return _LASPDFT.get_h2eff(self, mo_coeff=mo_coeff) + 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, states=states, **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, **kwargs) + + def multi_state(self, **kwargs): + ''' + In future will have to change this to consider the modal space selection, weights... + ''' + assert self.DoLASSI, "multi_state is only defined for post LAS methods" + return _LASPDFT.multi_state(self, **kwargs) - if DoLASSI: + multi_state_mix = multi_state + + if DoLASSI: _mc_class.DoLASSI = True _mc_class.rdmstmpfile = tempfile.NamedTemporaryFile(dir=lib.param.TMPDIR) - - else: _mc_class.DoLASSI = False - - if states is not None: _mc_class.states=states + + else: + _mc_class.DoLASSI = False + + if states is not None: _mc_class.states = states + if _mc_class.DoLASSI: - + ''' The cost of the RDM build is similar to LASSI diagonalization step. Therefore, calling it 2n time for n-states becomes prohibitively expensive. One alternative can be just call it once and store all the generated casdm1 and casdm2 and later on just call a reader function which will read the rdms from this temp file. ''' + def _store_rdms(self): # MRH: I made it loop over blocks of states to handle the O(N^5) memory cost # If there's enough memory it'll still do them all at once @@ -145,48 +173,52 @@ def _store_rdms(self): safety_factor = 1.3 mem_per_state = safety_factor*8*(2*(self.ncas**2) + 4*(self.ncas**4)) / 1e6 current_mem = lib.current_memory ()[0] + if current_mem > self.max_memory: - log.warn ("Current memory usage (%d MB) exceeds maximum memory (%d MB)", - current_mem, self.max_memory) + log.warn("Current memory usage (%d MB) exceeds maximum memory (%d MB)", + current_mem, self.max_memory) nblk = 1 else: nblk = max (1, int ((self.max_memory - current_mem) / mem_per_state)-1) + log.debug ('_store_rdms: looping over %d states at a time of %d total', nblk, len (self.e_states)) + rdmstmpfile = self.rdmstmpfile with h5py.File(rdmstmpfile, 'w') as f: - for i in range (0, len (self.e_states), nblk): - j = min (i+nblk, len (self.e_states)) + for i in range(0, len(self.e_states), nblk): + j = min(i + nblk, len(self.e_states)) rdm1s, rdm2s = lassi.root_make_rdm12s(self, self.ci, self.si, - state=list(range(i,j))) - for k in range (i, j): + state=list(range(i, j))) + for k in range(i, j): rdm1s_dname = f'rdm1s_{k}' f.create_dataset(rdm1s_dname, data=rdm1s[k]) rdm2s_dname = f'rdm2s_{k}' f.create_dataset(rdm2s_dname, data=rdm2s[k]) - rdm1s = rdm2s = None + rdm1s = rdm2s = None + + # # This code doesn't seem efficent, have to calculate the casdm1 and casdm2 in different functions. - # # This code doesn't seem efficent, have to calculate the casdm1 and casdm2 in different functions. # def make_one_casdm1s(self, ci=None, state=0, **kwargs): - # with lib.temporary_env (self, verbose=2): - # casdm1s = lassi.root_make_rdm12s (self, ci=ci, si=self.si, state=state)[0] - # return casdm1s + # with lib.temporary_env (self, verbose=2): + # casdm1s = lassi.root_make_rdm12s (self, ci=ci, si=self.si, state=state)[0] + # return casdm1s # def make_one_casdm2(self, ci=None, state=0, **kwargs): - # with lib.temporary_env (self, verbose=2): - # casdm2s = lassi.root_make_rdm12s (self, ci=ci, si=self.si, state=state)[1] - # return casdm2s.sum ((0,3)) - + # with lib.temporary_env (self, verbose=2): + # casdm2s = lassi.root_make_rdm12s (self, ci=ci, si=self.si, state=state)[1] + # return casdm2s.sum ((0,3)) + def make_one_casdm1s(self, ci=None, state=0, **kwargs): rdmstmpfile = self.rdmstmpfile return make_casdm1s(rdmstmpfile, 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, state).sum((0, 3)) + else: - make_one_casdm1s=mc.__class__.state_make_casdm1s - make_one_casdm2=mc.__class__.state_make_casdm2 + make_one_casdm1s = mc.__class__.state_make_casdm1s + make_one_casdm2 = mc.__class__.state_make_casdm2 # TODO: in pyscf-forge/pyscf/mcpdft/mcpdft.py::optimize_mcscf_, generalize the number # of return arguments. Then the redefinition below will be unnecessary. @@ -195,18 +227,18 @@ def optimize_mcscf_(self, mo_coeff=None, ci0=None, **kwargs): Has the same calling signature as the parent kernel method. ''' with _mcscf_env(self): if self.DoLASSI: - self._store_rdms() - self.fcisolver.nroots = len(self.e_states) if self.states is None else self.states self.e_states = self.e_roots + self.statlis = [x for x in range(len(self.e_states))] + self.fcisolver.nroots = len(self.e_states) if self.states is None else 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) _keys = pdft._keys.copy() - pdft.__dict__.update (mc.__dict__) + pdft.__dict__.update(mc.__dict__) pdft._keys = pdft._keys.union(_keys) return pdft - -