diff --git a/pyscf/ci/cisd.py b/pyscf/ci/cisd.py index 7c9fd08e27..c63f216a69 100644 --- a/pyscf/ci/cisd.py +++ b/pyscf/ci/cisd.py @@ -1143,7 +1143,6 @@ def _cp(a): if __name__ == '__main__': - from pyscf import gto from pyscf import ao2mo mol = gto.Mole() diff --git a/pyscf/df/grad/casscf.py b/pyscf/df/grad/casscf.py index 9ed69841eb..a5bce4404f 100644 --- a/pyscf/df/grad/casscf.py +++ b/pyscf/df/grad/casscf.py @@ -229,7 +229,6 @@ def _finalize(self): if __name__ == '__main__': - from pyscf import gto from pyscf import scf from pyscf import mcscf from pyscf import df diff --git a/pyscf/fci/direct_spin1_symm.py b/pyscf/fci/direct_spin1_symm.py index c59b4f096f..134916fbbf 100644 --- a/pyscf/fci/direct_spin1_symm.py +++ b/pyscf/fci/direct_spin1_symm.py @@ -615,6 +615,8 @@ def sym_allowed_indices(nelec, orbsym, wfnsym): class FCISolver(direct_spin1.FCISolver): + _keys = {'wfnsym', 'sym_allowed_idx'} + pspace_size = getattr(__config__, 'fci_direct_spin1_symm_FCI_pspace_size', 400) def __init__(self, mol=None, **kwargs): diff --git a/pyscf/gto/basis/__init__.py b/pyscf/gto/basis/__init__.py index 9d0706d359..b54517e582 100644 --- a/pyscf/gto/basis/__init__.py +++ b/pyscf/gto/basis/__init__.py @@ -669,20 +669,20 @@ def load_ecp(filename_or_basisname, symb): try: return parse_nwchem_ecp.parse(filename_or_basisname, symb) - except IndexError: - raise BasisNotFoundError(filename_or_basisname) except BasisNotFoundError as basis_err: pass + except Exception: + raise BasisNotFoundError(filename_or_basisname) try: return parse_nwchem_ecp.parse(filename_or_basisname) - except IndexError: - raise BasisNotFoundError(f'Invalid ECP {filename_or_basisname}') except BasisNotFoundError: pass + except Exception: + raise BasisNotFoundError(f'Invalid ECP {filename_or_basisname}') # Last, a trial to access Basis Set Exchange database - from pyscf.basis import bse + from pyscf.gto.basis import bse if bse.basis_set_exchange is not None: try: bse_obj = bse.basis_set_exchange.api.get_basis( @@ -708,7 +708,9 @@ def load_pseudo(filename_or_basisname, symb): try: return parse_cp2k_pp.parse(filename_or_basisname) - except IndexError: + except BasisNotFoundError: + raise + except Exception: raise BasisNotFoundError(f'Invalid PP {filename_or_basisname}') def _load_external(module, filename_or_basisname, symb, **kwargs): diff --git a/pyscf/gto/basis/parse_molpro.py b/pyscf/gto/basis/parse_molpro.py index 13e1fa2995..75e9aea8b8 100644 --- a/pyscf/gto/basis/parse_molpro.py +++ b/pyscf/gto/basis/parse_molpro.py @@ -24,7 +24,6 @@ import re import numpy -import numpy as np try: from pyscf.gto.basis.parse_nwchem import optimize_contraction @@ -155,10 +154,7 @@ def parse_terms(terms): try: coef = [float(x) for x in line[1:]] except ValueError: - if DISABLE_EVAL: - raise ValueError('Failed to parse ecp %s' % line) - else: - coef = list(eval(','.join(line[1:]))) + raise ValueError('Failed to parse ecp %s' % line) r_orders[order].append(coef) return r_orders diff --git a/pyscf/mcscf/addons.py b/pyscf/mcscf/addons.py index 91b04de061..c2a676cb22 100644 --- a/pyscf/mcscf/addons.py +++ b/pyscf/mcscf/addons.py @@ -840,10 +840,6 @@ def spin_square(casscf, mo_coeff=None, ci=None, ovlp=None): return ss, s*2+1 # A tag to label the derived FCI class -class StateAverageFCISolver: - pass -class StateAverageMixFCISolver(StateAverageFCISolver): - pass class StateSpecificFCISolver: pass # A tag to label the derived MCSCF class @@ -873,238 +869,241 @@ def state_average(casscf, weights=(0.5,0.5), wfnsym=None): of the SA-CASSCF method; see: Mol. Phys. 99, 103 (2001). ''' assert (abs(sum(weights)-1) < 1e-3) - fcibase_class = casscf.fcisolver.__class__ - has_spin_square = getattr(casscf.fcisolver, 'spin_square', None) - if wfnsym is None: - wfnsym = casscf.fcisolver.wfnsym - - class FakeCISolver(fcibase_class, StateAverageFCISolver): - _keys = set (('weights','e_states','_base_class')) - - def __init__(self, fcibase): - self.__dict__.update (fcibase.__dict__) - self.nroots = len(weights) - self.weights = weights - self.wfnsym = wfnsym - self.e_states = [None] - # MRH 09/09/2022: I turned the _base_class property into an - # attribute to prevent conflict with fix_spin_ dynamic class - self._base_class = fcibase_class - - def dump_flags(self, verbose=None): - if hasattr(fcibase_class, 'dump_flags'): - fcibase_class.dump_flags(self, verbose) - log = logger.new_logger(self, verbose) - log.info('State-average over %d states with weights %s', - len(self.weights), self.weights) - return self - - def kernel(self, h1, h2, norb, nelec, ci0=None, **kwargs): - if 'nroots' not in kwargs: - kwargs['nroots'] = self.nroots - - # call fcibase_class.kernel function because the attribute orbsym - # is available in self but undefined in fcibase object - e, c = fcibase_class.kernel(self, h1, h2, norb, nelec, ci0=ci0, - wfnsym=self.wfnsym, **kwargs) - self.e_states = e - - log = logger.new_logger(self, kwargs.get('verbose')) - if log.verbose >= logger.DEBUG: - if has_spin_square: - ss = self.states_spin_square(c, norb, nelec)[0] - for i, ei in enumerate(e): - log.debug('state %d E = %.15g S^2 = %.7f', i, ei, ss[i]) - else: - for i, ei in enumerate(e): - log.debug('state %d E = %.15g', i, ei) - return numpy.einsum('i,i->', e, self.weights), c - - def approx_kernel(self, h1, h2, norb, nelec, ci0=None, **kwargs): - if hasattr(fcibase_class, 'approx_kernel'): - e, c = fcibase_class.approx_kernel(self, h1, h2, norb, nelec, - ci0=ci0, nroots=self.nroots, - wfnsym=self.wfnsym, - **kwargs) - else: - e, c = fcibase_class.kernel(self, h1, h2, norb, nelec, ci0=ci0, - nroots=self.nroots, - wfnsym=self.wfnsym, **kwargs) - return numpy.einsum('i,i->', e, self.weights), c - - def states_make_rdm1(self, ci0, norb, nelec, *args, **kwargs): - fcibase = super() - dm1 = [fcibase.make_rdm1(c, norb, nelec, *args, **kwargs) for c in ci0] - return dm1 - - def make_rdm1(self, ci0, norb, nelec, *args, **kwargs): - return sum ([w * dm for w, dm in zip(self.weights, - self.states_make_rdm1(ci0, norb, nelec, *args, **kwargs))]) - - def states_make_rdm1s(self, ci0, norb, nelec, *args, **kwargs): - fcibase = super() - dm1a = [] - dm1b = [] - for c in ci0: - dm1s = fcibase.make_rdm1s(c, norb, nelec, *args, **kwargs) - dm1a.append (dm1s[0]) - dm1b.append (dm1s[1]) - return dm1a, dm1b - - def make_rdm1s(self, ci0, norb, nelec, *args, **kwargs): - dm1s = self.states_make_rdm1s(ci0, norb, nelec, *args, **kwargs) - dm1s = numpy.einsum ('r,srpq->spq', self.weights, dm1s) - return dm1s[0], dm1s[1] - - def states_make_rdm12(self, ci0, norb, nelec, *args, **kwargs): - fcibase = super() - rdm1 = [] - rdm2 = [] - for c in ci0: - dm1, dm2 = fcibase.make_rdm12(c, norb, nelec, *args, **kwargs) - rdm1.append (dm1) - rdm2.append (dm2) - return rdm1, rdm2 - - def make_rdm12(self, ci0, norb, nelec, *args, **kwargs): - rdm1, rdm2 = self.states_make_rdm12(ci0, norb, nelec, *args, **kwargs) - rdm1 = numpy.einsum ('r,rpq->pq', self.weights, rdm1) - rdm2 = numpy.einsum ('r,rpqst->pqst', self.weights, rdm2) - return rdm1, rdm2 - - def states_make_rdm12s(self, ci0, norb, nelec, *args, **kwargs): - fcibase = super() - dm1a, dm1b = [], [] - dm2aa, dm2ab, dm2bb = [], [], [] - for c in ci0: - dm1s, dm2s = fcibase.make_rdm12s(c, norb, nelec, *args, **kwargs) - dm1a.append(dm1s[0]) - dm1b.append(dm1s[1]) - dm2aa.append(dm2s[0]) - dm2ab.append(dm2s[1]) - dm2bb.append(dm2s[2]) - return (dm1a, dm1b), (dm2aa, dm2ab, dm2bb) - - def make_rdm12s(self, ci0, norb, nelec, *args, **kwargs): - rdm1s, rdm2s = self.states_make_rdm12s(ci0, norb, nelec, *args, **kwargs) - rdm1s = numpy.einsum ('r,srpq->spq', self.weights, rdm1s) - rdm2s = numpy.einsum ('r,srpqtu->spqtu', self.weights, rdm2s) - return rdm1s, rdm2s - - def states_trans_rdm12 (self, ci1, ci0, norb, nelec, *args, **kwargs): - fcibase = super() - tdm1 = [] - tdm2 = [] - for c1, c0 in zip (ci1, ci0): - dm1, dm2 = fcibase.trans_rdm12 (c1, c0, norb, nelec) - tdm1.append (dm1) - tdm2.append (dm2) - return tdm1, tdm2 - - def trans_rdm12 (self, ci1, ci0, norb, nelec, *args, **kwargs): - tdm1, tdm2 = self.states_trans_rdm12 (ci1, ci0, norb, nelec, *args, **kwargs) - tdm1 = numpy.einsum ('r,rpq->pq', self.weights, tdm1) - tdm2 = numpy.einsum ('r,rpqst->pqst', self.weights, tdm2) - return tdm1, tdm2 - - if has_spin_square: - def spin_square(self, ci0, norb, nelec, *args, **kwargs): - ss, multip = self.states_spin_square(ci0, norb, nelec, *args, **kwargs) - weights = self.weights - return numpy.dot(ss, weights), numpy.dot(multip, weights) - - def states_spin_square(self, ci0, norb, nelec, *args, **kwargs): - fcibase = super() - s = [fcibase.spin_square(ci0[i], norb, nelec, *args, **kwargs) - for i, wi in enumerate(self.weights)] - return [x[0] for x in s], [x[1] for x in s] - - # No recursion of FakeCISolver is allowed! - if isinstance (casscf.fcisolver, StateAverageFCISolver): - fcisolver = casscf.fcisolver + fcisolver = casscf.fcisolver + + # No recursion is allowed! + if isinstance (fcisolver, StateAverageFCISolver): fcisolver.nroots = len(weights) fcisolver.weights = weights else: - fcisolver = FakeCISolver(casscf.fcisolver) + def spin_square(self, ci0, norb, nelec, *args, **kwargs): + ss, multip = self.states_spin_square(ci0, norb, nelec, *args, **kwargs) + weights = self.weights + return numpy.dot(ss, weights), numpy.dot(multip, weights) + + def states_spin_square(self, ci0, norb, nelec, *args, **kwargs): + fcibase = super(StateAverageFCISolver, self) + s = [fcibase.spin_square(ci0[i], norb, nelec, *args, **kwargs) + for i, wi in enumerate(self.weights)] + return [x[0] for x in s], [x[1] for x in s] + + methods = {} + if hasattr(fcisolver, 'spin_square'): + methods['spin_square'] = spin_square + methods['states_spin_square'] = states_spin_square + + fcisolver = lib.set_class(StateAverageFCISolver(fcisolver, weights, wfnsym), + (StateAverageFCISolver, fcisolver.__class__), + attrs=methods) return _state_average_mcscf_solver(casscf, fcisolver) +class StateAverageFCISolver: + __name_mixin__ = 'StateAverage' + + _keys = {'weights', 'e_states'} + + def __init__(self, fcibase, weights, wfnsym): + self.__dict__.update (fcibase.__dict__) + self.nroots = len(weights) + self.weights = weights + if wfnsym is not None: + self.wfnsym = wfnsym + self.e_states = [None] + + def undo_state_average(self): + obj = lib.view(self, lib.drop_class(self.__class__, StateAverageFCISolver)) + del obj.weights + del obj.e_states + return obj + + def dump_flags(self, verbose=None): + super().dump_flags(verbose) + log = logger.new_logger(self, verbose) + log.info('State-average over %d states with weights %s', + len(self.weights), self.weights) + return self + + def kernel(self, h1, h2, norb, nelec, ci0=None, **kwargs): + if 'nroots' not in kwargs: + kwargs['nroots'] = self.nroots + + fcibase = super() + # call fcibase_class.kernel function because the attribute orbsym + # is available in self but undefined in fcibase object + e, c = fcibase.kernel(h1, h2, norb, nelec, ci0=ci0, + wfnsym=self.wfnsym, **kwargs) + self.e_states = e + + log = logger.new_logger(self, kwargs.get('verbose')) + if log.verbose >= logger.DEBUG: + if hasattr(fcibase, 'spin_square'): + ss = self.states_spin_square(c, norb, nelec)[0] + for i, ei in enumerate(e): + log.debug('state %d E = %.15g S^2 = %.7f', i, ei, ss[i]) + else: + for i, ei in enumerate(e): + log.debug('state %d E = %.15g', i, ei) + return numpy.einsum('i,i->', e, self.weights), c + + def approx_kernel(self, h1, h2, norb, nelec, ci0=None, **kwargs): + fcibase = super() + if hasattr(fcibase, 'approx_kernel'): + e, c = fcibase.approx_kernel(h1, h2, norb, nelec, + ci0=ci0, nroots=self.nroots, + wfnsym=self.wfnsym, **kwargs) + else: + e, c = fcibase.kernel(h1, h2, norb, nelec, ci0=ci0, + nroots=self.nroots, wfnsym=self.wfnsym, **kwargs) + return numpy.einsum('i,i->', e, self.weights), c + + def states_make_rdm1(self, ci0, norb, nelec, *args, **kwargs): + fcibase = super() + dm1 = [fcibase.make_rdm1(c, norb, nelec, *args, **kwargs) for c in ci0] + return dm1 + + def make_rdm1(self, ci0, norb, nelec, *args, **kwargs): + return sum ([w * dm for w, dm in zip(self.weights, + self.states_make_rdm1(ci0, norb, nelec, *args, **kwargs))]) + + def states_make_rdm1s(self, ci0, norb, nelec, *args, **kwargs): + fcibase = super() + dm1a = [] + dm1b = [] + for c in ci0: + dm1s = fcibase.make_rdm1s(c, norb, nelec, *args, **kwargs) + dm1a.append (dm1s[0]) + dm1b.append (dm1s[1]) + return dm1a, dm1b + + def make_rdm1s(self, ci0, norb, nelec, *args, **kwargs): + dm1s = self.states_make_rdm1s(ci0, norb, nelec, *args, **kwargs) + dm1s = numpy.einsum ('r,srpq->spq', self.weights, dm1s) + return dm1s[0], dm1s[1] + + def states_make_rdm12(self, ci0, norb, nelec, *args, **kwargs): + fcibase = super() + rdm1 = [] + rdm2 = [] + for c in ci0: + dm1, dm2 = fcibase.make_rdm12(c, norb, nelec, *args, **kwargs) + rdm1.append (dm1) + rdm2.append (dm2) + return rdm1, rdm2 + + def make_rdm12(self, ci0, norb, nelec, *args, **kwargs): + rdm1, rdm2 = self.states_make_rdm12(ci0, norb, nelec, *args, **kwargs) + rdm1 = numpy.einsum ('r,rpq->pq', self.weights, rdm1) + rdm2 = numpy.einsum ('r,rpqst->pqst', self.weights, rdm2) + return rdm1, rdm2 + + def states_make_rdm12s(self, ci0, norb, nelec, *args, **kwargs): + fcibase = super() + dm1a, dm1b = [], [] + dm2aa, dm2ab, dm2bb = [], [], [] + for c in ci0: + dm1s, dm2s = fcibase.make_rdm12s(c, norb, nelec, *args, **kwargs) + dm1a.append(dm1s[0]) + dm1b.append(dm1s[1]) + dm2aa.append(dm2s[0]) + dm2ab.append(dm2s[1]) + dm2bb.append(dm2s[2]) + return (dm1a, dm1b), (dm2aa, dm2ab, dm2bb) + + def make_rdm12s(self, ci0, norb, nelec, *args, **kwargs): + rdm1s, rdm2s = self.states_make_rdm12s(ci0, norb, nelec, *args, **kwargs) + rdm1s = numpy.einsum ('r,srpq->spq', self.weights, rdm1s) + rdm2s = numpy.einsum ('r,srpqtu->spqtu', self.weights, rdm2s) + return rdm1s, rdm2s + + def states_trans_rdm12 (self, ci1, ci0, norb, nelec, *args, **kwargs): + fcibase = super() + tdm1 = [] + tdm2 = [] + for c1, c0 in zip (ci1, ci0): + dm1, dm2 = fcibase.trans_rdm12 (c1, c0, norb, nelec) + tdm1.append (dm1) + tdm2.append (dm2) + return tdm1, tdm2 + + def trans_rdm12 (self, ci1, ci0, norb, nelec, *args, **kwargs): + tdm1, tdm2 = self.states_trans_rdm12 (ci1, ci0, norb, nelec, *args, **kwargs) + tdm1 = numpy.einsum ('r,rpq->pq', self.weights, tdm1) + tdm2 = numpy.einsum ('r,rpqst->pqst', self.weights, tdm2) + return tdm1, tdm2 + def _state_average_mcscf_solver(casscf, fcisolver): '''A common routine for function state_average and state_average_mix to generate state-average MCSCF solver. ''' - mcscfbase_class = casscf.__class__ if isinstance (casscf, StateAverageMCSCFSolver): - raise TypeError('mc is not base MCSCF solver\n' - 'state_average function cannot work with decorated ' - 'MCSCF solver %s.\nYou can restore the base MCSCF ' - 'then call state_average function, eg\n' - ' mc = %s.%s(mc._scf, %s, %s)\n' - ' mc.state_average_()\n' % - (mcscfbase_class, mcscfbase_class.__base__.__module__, - mcscfbase_class.__base__.__name__, casscf.ncas, casscf.nelecas)) - - has_spin_square = getattr(fcisolver, 'spin_square', None) - - print(mcscfbase_class.__mro__) - class StateAverageMCSCF(mcscfbase_class, StateAverageMCSCFSolver): - _keys = set (('weights', '_base_class')) - - def __init__(self, my_mc): - self.__dict__.update (my_mc.__dict__) - self.fcisolver = fcisolver - - @property - def _base_class (self): - ''' for convenience; this is equal to mcscfbase_class ''' - return self.__class__.__bases__[0] - - @property - def weights (self): - ''' I want these to be accessible but not separable from fcisolver.weights ''' - return self.fcisolver.weights - - @weights.setter - def weights (self, x): - self.fcisolver.weights = x - return self.fcisolver.weights - - @property - def e_average(self): - return numpy.dot(self.fcisolver.weights, self.fcisolver.e_states) - - @property - def e_states(self): - return self.fcisolver.e_states - - def _finalize(self): - mcscfbase_class._finalize(self) - # Do not overwrite self.e_tot because .e_tot needs to be used in - # geometry optimization. self.e_states can be used to access the - # energy of each state - #self.e_tot = self.fcisolver.e_states - logger.note(self, 'CASCI state-averaged energy = %.15g', self.e_average) - logger.note(self, 'CASCI energy for each state') - if has_spin_square: - ss = self.fcisolver.states_spin_square(self.ci, self.ncas, - self.nelecas)[0] - for i, ei in enumerate(self.e_states): - logger.note(self, ' State %d weight %g E = %.15g S^2 = %.7f', - i, self.weights[i], ei, ss[i]) - else: - for i, ei in enumerate(self.e_states): - logger.note(self, ' State %d weight %g E = %.15g', - i, self.weights[i], ei) - return self - - def nuc_grad_method (self, state=None): - if callable (getattr (self, '_state_average_nuc_grad_method', None)): - return self._state_average_nuc_grad_method (state=state) - else: # Avoid messing up state-average CASCI - return self._base_class.nuc_grad_method (self) + casscf = casscf.undo_state_average() + + return lib.set_class(StateAverageMCSCF(casscf, fcisolver), + (StateAverageMCSCF, casscf.__class__)) + +class StateAverageMCSCF(StateAverageMCSCFSolver): + __name_mixin__ = 'StateAverage' + + def __init__(self, my_mc, fcisolver): + self.__dict__.update (my_mc.__dict__) + self.fcisolver = fcisolver + + def undo_state_average(self): + obj = lib.view(self, lib.drop_class(self.__class__, StateAverageMCSCF)) + if isinstance(self.fcisolver, StateAverageFCISolver): + obj.fcisolver = self.fcisolver.undo_state_average() + return obj + + @property + def _base_class (self): + ''' for convenience; this is equal to mcscfbase_class ''' + return self.__class__.__bases__[0] + + @property + def weights (self): + ''' I want these to be accessible but not separable from fcisolver.weights ''' + return self.fcisolver.weights + + @weights.setter + def weights (self, x): + self.fcisolver.weights = x + return self.fcisolver.weights + + @property + def e_average(self): + return numpy.dot(self.fcisolver.weights, self.fcisolver.e_states) + + @property + def e_states(self): + return self.fcisolver.e_states + + def _finalize(self): + super()._finalize() + # Do not overwrite self.e_tot because .e_tot needs to be used in + # geometry optimization. self.e_states can be used to access the + # energy of each state + #self.e_tot = self.fcisolver.e_states + logger.note(self, 'CASCI state-averaged energy = %.15g', self.e_average) + logger.note(self, 'CASCI energy for each state') + if hasattr(self, 'spin_square'): + ss = self.fcisolver.states_spin_square(self.ci, self.ncas, + self.nelecas)[0] + for i, ei in enumerate(self.e_states): + logger.note(self, ' State %d weight %g E = %.15g S^2 = %.7f', + i, self.weights[i], ei, ss[i]) + else: + for i, ei in enumerate(self.e_states): + logger.note(self, ' State %d weight %g E = %.15g', + i, self.weights[i], ei) + return self - Gradients = nuc_grad_method + def nuc_grad_method (self, state=None): + if callable (getattr (self, '_state_average_nuc_grad_method', None)): + return self._state_average_nuc_grad_method (state=state) + else: # Avoid messing up state-average CASCI + return self._base_class.nuc_grad_method (self) - return StateAverageMCSCF(casscf) + Gradients = nuc_grad_method def state_average_(casscf, weights=(0.5,0.5), wfnsym=None): ''' Inplace version of state_average ''' @@ -1121,74 +1120,79 @@ def state_specific_(casscf, state=1, wfnsym=None): state : int 0 for ground state; 1 for first excited state. ''' - fcibase_class = casscf.fcisolver.__class__ - if fcibase_class.__name__ == 'FakeCISolver': - raise TypeError('mc.fcisolver is not base FCI solver\n' - 'state_specific function cannot work with decorated ' - 'fcisolver %s.\nYou can restore the base fcisolver ' - 'then call state_specific function, eg\n' - ' mc.fcisolver = %s.%s(mc.mol)\n' - ' mc.state_specific_()\n' % - (casscf.fcisolver, fcibase_class.__base__.__module__, - fcibase_class.__base__.__name__)) - if wfnsym is None: - wfnsym = casscf.fcisolver.wfnsym - - class FakeCISolver(fcibase_class, StateSpecificFCISolver): - def __init__(self): - self.nroots = state+1 - self._civec = None + fcisolver = casscf.fcisolver + if isinstance(fcisolver, StateSpecificFCISolver): + fcisolver.nroots = state+1 + fcisolver._civec = None + if wfnsym is not None: + fcisolver.wfnsym = wfnsym + elif isinstance(fcisolver, StateAverageFCISolver): + fcisolver = fcisolver.undo_state_average() + else: + fcisolver = lib.set_class(StateSpecificFCISolver(fcisolver, state, wfnsym), + (StateSpecificFCISolver, fcisolver.__class__)) + casscf.fcisolver = fcisolver + return casscf +state_specific = state_specific_ + +class StateSpecificFCISolver: + __name_mixin__ = 'StateSpecific' + + _keys = {'state', 'nroots'} + + def __init__(self, fcibase, state, wfnsym): + self.__dict__.update(fcibase.__dict__) + self.state = state + self.nroots = state+1 + self._civec = None + if wfnsym is not None: self.wfnsym = wfnsym - def kernel(self, h1, h2, norb, nelec, ci0=None, **kwargs): - if self._civec is not None: - ci0 = self._civec - e, c = fcibase_class.kernel(self, h1, h2, norb, nelec, ci0=ci0, - nroots=self.nroots, wfnsym=self.wfnsym, - **kwargs) - if state == 0: - e = [e] - c = [c] + def undo_state_specific(self): + obj = lib.view(self, lib.drop_class(self.__class__, StateSpecificFCISolver)) + del obj._civec + return obj + + def kernel(self, h1, h2, norb, nelec, ci0=None, **kwargs): + if self._civec is not None: + ci0 = self._civec + fcibase = super() + e, c = fcibase.kernel(h1, h2, norb, nelec, ci0=ci0, + nroots=self.nroots, wfnsym=self.wfnsym, **kwargs) + state = self.state + if state == 0: + e = [e] + c = [c] + self._civec = c + log = logger.new_logger(self, kwargs.get('verbose')) + if log.verbose >= logger.DEBUG: + if hasattr(fcibase, 'spin_square'): + ss = fcibase.spin_square(c[state], norb, nelec) + log.debug('state %d E = %.15g S^2 = %.7f', + state, e[state], ss[0]) + else: + log.debug('state %d E = %.15g', state, e[state]) + return e[state], c[state] + + def approx_kernel(self, h1, h2, norb, nelec, ci0=None, **kwargs): + if self._civec is not None: + ci0 = self._civec + fcibase = super() + if hasattr(fcibase, 'approx_kernel'): + e, c = fcibase.approx_kernel(h1, h2, norb, nelec, + ci0=ci0, nroots=self.nroots, + wfnsym=self.wfnsym, **kwargs) + else: + e, c = fcibase.kernel(h1, h2, norb, nelec, ci0=ci0, + nroots=self.nroots, wfnsym=self.wfnsym, **kwargs) + state = self.state + if state == 0: + self._civec = [c] + return e, c + else: self._civec = c - log = logger.new_logger(self, kwargs.get('verbose')) - if log.verbose >= logger.DEBUG: - if getattr(fcibase_class, 'spin_square', None): - fcibase = super() - ss = fcibase.spin_square(c[state], norb, nelec) - log.debug('state %d E = %.15g S^2 = %.7f', - state, e[state], ss[0]) - else: - log.debug('state %d E = %.15g', state, e[state]) return e[state], c[state] - def approx_kernel(self, h1, h2, norb, nelec, ci0=None, **kwargs): - if self._civec is not None: - ci0 = self._civec - if hasattr(fcibase_class, 'approx_kernel'): - e, c = fcibase_class.approx_kernel(self, h1, h2, norb, nelec, - ci0=ci0, nroots=self.nroots, - wfnsym=self.wfnsym, - **kwargs) - else: - e, c = fcibase_class.kernel(self, h1, h2, norb, nelec, ci0=ci0, - nroots=self.nroots, - wfnsym=self.wfnsym, **kwargs) - if state == 0: - self._civec = [c] - return e, c - else: - self._civec = c - return e[state], c[state] - - fcisolver = FakeCISolver() - fcisolver.__dict__.update(casscf.fcisolver.__dict__) - fcisolver.nroots = state + 1 - if wfnsym is not None: - fcisolver.wfnsym = wfnsym - casscf.fcisolver = fcisolver - return casscf -state_specific = state_specific_ - class StateAverageMixFCISolver_solver_args: def __init__(self, data): self._data = data @@ -1201,266 +1205,298 @@ def __getitem__(self, key): # Handle data = None class StateAverageMixFCISolver_state_args (StateAverageMixFCISolver_solver_args): pass -def state_average_mix(casscf, fcisolvers, weights=(0.5,0.5)): - '''State-average CASSCF over multiple FCI solvers. - ''' - fcibase_class = fcisolvers[0].__class__ - nroots = sum(solver.nroots for solver in fcisolvers) - assert (nroots == len(weights)) - has_spin_square = all(getattr(solver, 'spin_square', None) - for solver in fcisolvers) - has_large_ci = all(getattr(solver, 'large_ci', None) - for solver in fcisolvers) - has_transform_ci = all(getattr(solver, 'transform_ci_for_orbital_rotation', None) - for solver in fcisolvers) +class StateAverageMixFCISolver(StateAverageFCISolver): + __name_mixin__ = 'StateAverageMix' + + _keys = set (('weights','e_states','fcisolvers')) + _solver_args = StateAverageMixFCISolver_solver_args _state_args = StateAverageMixFCISolver_state_args - class FakeCISolver(fcibase_class, StateAverageMixFCISolver): - _keys = set (('weights','e_states','_base_class','fcisolvers')) - - def __init__(self, mol): - fcibase_class.__init__(self, mol) - self.nroots = len(weights) - self.weights = weights - self.e_states = [None] - self.fcisolvers = fcisolvers - - @property - def _base_class (self): - ''' for convenience; this is equal to mcscfbase_class ''' - return self.__class__.__bases__[0] - - # MRH 06/24/2020: I need these functions in newton_casscf! - # TODO: handle things like linkstr somehow (variables that - # have to be different for different solvers or ci vecs) - def _loop_solver(self, *args, **kwargs): - p0 = 0 - for ix, solver in enumerate (self.fcisolvers): + + def __init__(self, fcisolvers, weights): + self.__dict__.update(fcisolvers[0].__dict__) + self.nroots = len(weights) + self.weights = weights + self.e_states = [None] + self.fcisolvers = fcisolvers + + def undo_smearing(self): + obj = lib.view(self, lib.drop_class(self.__class__, StateAverageMixFCISolver)) + del obj.weights + del obj.e_states + del obj.fcisolvers + return obj + + @property + def _base_class (self): + return self.fcisolvers[0].__base__ + + # MRH 06/24/2020: I need these functions in newton_casscf! + # TODO: handle things like linkstr somehow (variables that + # have to be different for different solvers or ci vecs) + def _loop_solver(self, *args, **kwargs): + _solver_args = self._solver_args + _state_args = self._state_args + p0 = 0 + for ix, solver in enumerate (self.fcisolvers): + my_args = [] + for arg in args: + if isinstance (arg, _state_args): + my_arg = arg[p0:p0+solver.nroots] + if solver.nroots == 1 and my_arg is not None: my_arg = my_arg[0] + my_args.append (my_arg) + elif isinstance (arg, _solver_args): + my_args.append (arg[ix]) + else: + my_args.append (arg) + my_kwargs = {} + for key, item in kwargs.items (): + if isinstance (item, _state_args): + my_arg = item[p0:p0+solver.nroots] + if solver.nroots == 1 and my_arg is not None: my_arg = my_arg[0] + my_kwargs[key] = my_arg + elif isinstance (item, _solver_args): + my_kwargs[key] = item[ix] + else: + my_kwargs[key] = item + yield solver, my_args, my_kwargs + p0 += solver.nroots + + def _loop_civecs(self, *args, **kwargs): + _solver_args = self._solver_args + _state_args = self._state_args + p0 = 0 + for i, solver in enumerate (self.fcisolvers): + for j in range(p0, p0+solver.nroots): my_args = [] for arg in args: if isinstance (arg, _state_args): - my_arg = arg[p0:p0+solver.nroots] - if solver.nroots == 1 and my_arg is not None: my_arg = my_arg[0] - my_args.append (my_arg) + my_args.append (arg[j]) elif isinstance (arg, _solver_args): - my_args.append (arg[ix]) + my_args.append (arg[i]) else: my_args.append (arg) my_kwargs = {} for key, item in kwargs.items (): if isinstance (item, _state_args): - my_arg = item[p0:p0+solver.nroots] - if solver.nroots == 1 and my_arg is not None: my_arg = my_arg[0] - my_kwargs[key] = my_arg + my_kwargs[key] = item[j] elif isinstance (item, _solver_args): - my_kwargs[key] = item[ix] + my_kwargs[key] = item[i] else: my_kwargs[key] = item yield solver, my_args, my_kwargs - p0 += solver.nroots - - def _loop_civecs(self, *args, **kwargs): - p0 = 0 - for i, solver in enumerate (self.fcisolvers): - for j in range(p0, p0+solver.nroots): - my_args = [] - for arg in args: - if isinstance (arg, _state_args): - my_args.append (arg[j]) - elif isinstance (arg, _solver_args): - my_args.append (arg[i]) - else: - my_args.append (arg) - my_kwargs = {} - for key, item in kwargs.items (): - if isinstance (item, _state_args): - my_kwargs[key] = item[j] - elif isinstance (item, _solver_args): - my_kwargs[key] = item[i] - else: - my_kwargs[key] = item - yield solver, my_args, my_kwargs - p0 += solver.nroots - - def _get_nelec(self, solver, nelec): - # FCISolver does not need this function. Some external solver may not - # have the function to handle nelec and spin - # MRH 06/24/2020: Yes, FCISolver DOES need this function! - if solver.spin is not None: - nelec = numpy.sum(nelec) - nelec = (nelec+solver.spin)//2, (nelec-solver.spin)//2 - return nelec - - def _collect(self, fname, *args, **kwargs): - for solver, args, kwargs in self._loop_civecs(*args, **kwargs): - fn = getattr(solver, fname) - yield fn(*args, **kwargs) - - def kernel(self, h1, h2, norb, nelec, ci0=None, verbose=0, **kwargs): - # Note self.orbsym is initialized lazily in mc1step_symm.kernel function - log = logger.new_logger(self, verbose) - es = [] - cs = [] - for solver, my_args, my_kwargs in self._loop_solver(_state_args (ci0)): - c0 = my_args[0] + p0 += solver.nroots + + def _get_nelec(self, solver, nelec): + # FCISolver does not need this function. Some external solver may not + # have the function to handle nelec and spin + # MRH 06/24/2020: Yes, FCISolver DOES need this function! + if solver.spin is not None: + nelec = numpy.sum(nelec) + nelec = (nelec+solver.spin)//2, (nelec-solver.spin)//2 + return nelec + + def _collect(self, fname, *args, **kwargs): + for solver, args, kwargs in self._loop_civecs(*args, **kwargs): + fn = getattr(solver, fname) + yield fn(*args, **kwargs) + + def kernel(self, h1, h2, norb, nelec, ci0=None, verbose=0, **kwargs): + # Note self.orbsym is initialized lazily in mc1step_symm.kernel function + _state_args = self._state_args + log = logger.new_logger(self, verbose) + es = [] + cs = [] + for solver, my_args, my_kwargs in self._loop_solver(_state_args (ci0)): + c0 = my_args[0] + e, c = solver.kernel(h1, h2, norb, self._get_nelec(solver, nelec), ci0=c0, + orbsym=self.orbsym, verbose=log, **kwargs) + if solver.nroots == 1: + es.append(e) + cs.append(c) + else: + es.extend(e) + cs.extend(c) + self.e_states = es + self.converged = numpy.all(getattr(sol, 'converged', True) + for sol in self.fcisolvers) + + if log.verbose >= logger.DEBUG: + if all(hasattr(solver, 'spin_square') for solver in self.fcisolvers): + ss = self.states_spin_square(cs, norb, nelec)[0] + for i, ei in enumerate(es): + log.debug('state %d E = %.15g S^2 = %.7f', i, ei, ss[i]) + else: + for i, ei in enumerate(es): + log.debug('state %d E = %.15g', i, ei) + return numpy.einsum('i,i', numpy.array(es), self.weights), cs + + def approx_kernel(self, h1, h2, norb, nelec, ci0=None, **kwargs): + _state_args = self._state_args + es = [] + cs = [] + for ix, (solver, my_args, my_kwargs) in enumerate (self._loop_solver(_state_args (ci0))): + c0 = my_args[0] + if hasattr(solver, 'approx_kernel'): + e, c = solver.approx_kernel(h1, h2, norb, self._get_nelec(solver, nelec), ci0=c0, + orbsym=self.orbsym, **kwargs) + else: e, c = solver.kernel(h1, h2, norb, self._get_nelec(solver, nelec), ci0=c0, - orbsym=self.orbsym, verbose=log, **kwargs) - if solver.nroots == 1: - es.append(e) - cs.append(c) - else: - es.extend(e) - cs.extend(c) - self.e_states = es - self.converged = numpy.all(getattr(sol, 'converged', True) - for sol in fcisolvers) - - if log.verbose >= logger.DEBUG: - if has_spin_square: - ss = self.states_spin_square(cs, norb, nelec)[0] - for i, ei in enumerate(es): - log.debug('state %d E = %.15g S^2 = %.7f', i, ei, ss[i]) - else: - for i, ei in enumerate(es): - log.debug('state %d E = %.15g', i, ei) - return numpy.einsum('i,i', numpy.array(es), weights), cs - - def approx_kernel(self, h1, h2, norb, nelec, ci0=None, **kwargs): - es = [] - cs = [] - for ix, (solver, my_args, my_kwargs) in enumerate (self._loop_solver(_state_args (ci0))): - c0 = my_args[0] - if hasattr(solver, 'approx_kernel'): - e, c = solver.approx_kernel(h1, h2, norb, self._get_nelec(solver, nelec), ci0=c0, - orbsym=self.orbsym, **kwargs) - else: - e, c = solver.kernel(h1, h2, norb, self._get_nelec(solver, nelec), ci0=c0, - orbsym=self.orbsym, **kwargs) - if solver.nroots == 1: - es.append(e) - cs.append(c) - else: - es.extend(e) - cs.extend(c) - return numpy.einsum('i,i->', es, weights), cs + orbsym=self.orbsym, **kwargs) + if solver.nroots == 1: + es.append(e) + cs.append(c) + else: + es.extend(e) + cs.extend(c) + return numpy.einsum('i,i->', es, self.weights), cs + + def states_make_rdm1 (self, ci0, norb, nelec, link_index=None, **kwargs): + _solver_args = self._solver_args + _state_args = self._state_args + ci0 = _state_args (ci0) + link_index = _solver_args (link_index) + nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) + return [dm for dm in self._collect ('make_rdm1', ci0, norb, nelec, link_index=link_index, **kwargs)] + + def make_rdm1(self, ci0, norb, nelec, link_index=None, **kwargs): + dm1 = self.states_make_rdm1 (ci0, norb, nelec, link_index=link_index, **kwargs) + return numpy.einsum ('r,rpq->pq', self.weights, dm1) + + def states_make_rdm1s (self, ci0, norb, nelec, link_index=None, **kwargs): + _solver_args = self._solver_args + _state_args = self._state_args + ci0 = _state_args (ci0) + link_index = _solver_args (link_index) + nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) + dm1a = [] + dm1b = [] + for dm1s in self._collect ('make_rdm1s', ci0, norb, nelec, link_index=link_index, **kwargs): + dm1a.append (dm1s[0]) + dm1b.append (dm1s[1]) + return dm1a, dm1b + + def make_rdm1s(self, ci0, norb, nelec, link_index=None, **kwargs): + dm1a, dm1b = self.states_make_rdm1s (ci0, norb, nelec, link_index=link_index, **kwargs) + dm1s = numpy.einsum ('r,srpq->spq', self.weights, [dm1a, dm1b]) + return dm1s[0], dm1s[1] + + def states_make_rdm12 (self, ci0, norb, nelec, link_index=None, **kwargs): + _solver_args = self._solver_args + _state_args = self._state_args + ci0 = _state_args (ci0) + link_index = _solver_args (link_index) + nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) + rdm1 = [] + rdm2 = [] + for dm1, dm2 in self._collect ('make_rdm12', ci0, norb, nelec, link_index=link_index, **kwargs): + rdm1.append (dm1) + rdm2.append (dm2) + return rdm1, rdm2 + + def make_rdm12(self, ci0, norb, nelec, link_index=None, **kwargs): + rdm1, rdm2 = self.states_make_rdm12 (ci0, norb, nelec, link_index=link_index, **kwargs) + rdm1 = numpy.einsum ('r,rpq->pq', self.weights, rdm1) + rdm2 = numpy.einsum ('r,rpqst->pqst', self.weights, rdm2) + return rdm1, rdm2 + + def states_make_rdm12s(self, ci0, norb, nelec, link_index=None, **kwargs): + _solver_args = self._solver_args + _state_args = self._state_args + ci0 = _state_args (ci0) + link_index = _solver_args (link_index) + nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) + dm1a, dm1b = [], [] + dm2aa, dm2ab, dm2bb = [], [], [] + for dm1s, dm2s in self._collect ('make_rdm12s', ci0, norb, nelec, link_index=link_index, **kwargs): + dm1a.append(dm1s[0]) + dm1b.append(dm1s[1]) + dm2aa.append(dm2s[0]) + dm2ab.append(dm2s[1]) + dm2bb.append(dm2s[2]) + return (dm1a, dm1b), (dm2aa, dm2ab, dm2bb) + + def make_rdm12s(self, ci0, norb, nelec, link_index=None, **kwargs): + rdm1s, rdm2s = self.states_make_rdm12s(ci0, norb, nelec, link_index=link_index, **kwargs) + rdm1s = numpy.einsum ('r,srpq->spq', self.weights, rdm1s) + rdm2s = numpy.einsum ('r,srpqtu->spqtu', self.weights, rdm2s) + return rdm1s, rdm2s + + # TODO: linkstr support + def states_trans_rdm12 (self, ci1, ci0, norb, nelec, link_index=None, **kwargs): + _solver_args = self._solver_args + _state_args = self._state_args + ci1 = _state_args (ci1) + ci0 = _state_args (ci0) + link_index = _solver_args (link_index) + nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) + tdm1 = [] + tdm2 = [] + for dm1, dm2 in self._collect ('trans_rdm12', ci1, ci0, norb, nelec, link_index=link_index, **kwargs): + tdm1.append (dm1) + tdm2.append (dm2) + return tdm1, tdm2 + + def trans_rdm12 (self, ci1, ci0, norb, nelec, link_index=None, **kwargs): + tdm1, tdm2 = self.states_trans_rdm12 (ci1, ci0, norb, nelec, link_index=link_index, **kwargs) + tdm1 = numpy.einsum ('r,rpq->pq', self.weights, tdm1) + tdm2 = numpy.einsum ('r,rpqst->pqst', self.weights, tdm2) + return tdm1, tdm2 - def states_make_rdm1 (self, ci0, norb, nelec, link_index=None, **kwargs): - ci0 = _state_args (ci0) - link_index = _solver_args (link_index) - nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) - return [dm for dm in self._collect ('make_rdm1', ci0, norb, nelec, link_index=link_index, **kwargs)] +def state_average_mix(casscf, fcisolvers, weights=(0.5,0.5)): + '''State-average CASSCF over multiple FCI solvers. + ''' + nroots = sum(solver.nroots for solver in fcisolvers) + assert (nroots == len(weights)) - def make_rdm1(self, ci0, norb, nelec, link_index=None, **kwargs): - dm1 = self.states_make_rdm1 (ci0, norb, nelec, link_index=link_index, **kwargs) - return numpy.einsum ('r,rpq->pq', self.weights, dm1) + methods = {} + if all(hasattr(solver, 'spin_square') for solver in fcisolvers): + def spin_square(self, ci0, norb, nelec, *args, **kwargs): + ss, multip = self.states_spin_square(ci0, norb, nelec, *args, **kwargs) + weights = self.weights + return numpy.dot(ss, weights), numpy.dot(multip, weights) - def states_make_rdm1s (self, ci0, norb, nelec, link_index=None, **kwargs): + def states_spin_square(self, ci0, norb, nelec, *args, **kwargs): + _solver_args = self._solver_args + _state_args = self._state_args ci0 = _state_args (ci0) - link_index = _solver_args (link_index) nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) - dm1a = [] - dm1b = [] - for dm1s in self._collect ('make_rdm1s', ci0, norb, nelec, link_index=link_index, **kwargs): - dm1a.append (dm1s[0]) - dm1b.append (dm1s[1]) - return dm1a, dm1b - - def make_rdm1s(self, ci0, norb, nelec, link_index=None, **kwargs): - dm1a, dm1b = self.states_make_rdm1s (ci0, norb, nelec, link_index=link_index, **kwargs) - dm1s = numpy.einsum ('r,srpq->spq', self.weights, [dm1a, dm1b]) - return dm1s[0], dm1s[1] - - def states_make_rdm12 (self, ci0, norb, nelec, link_index=None, **kwargs): - ci0 = _state_args (ci0) - link_index = _solver_args (link_index) - nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) - rdm1 = [] - rdm2 = [] - for dm1, dm2 in self._collect ('make_rdm12', ci0, norb, nelec, link_index=link_index, **kwargs): - rdm1.append (dm1) - rdm2.append (dm2) - return rdm1, rdm2 - - def make_rdm12(self, ci0, norb, nelec, link_index=None, **kwargs): - rdm1, rdm2 = self.states_make_rdm12 (ci0, norb, nelec, link_index=link_index, **kwargs) - rdm1 = numpy.einsum ('r,rpq->pq', self.weights, rdm1) - rdm2 = numpy.einsum ('r,rpqst->pqst', self.weights, rdm2) - return rdm1, rdm2 - - def states_make_rdm12s(self, ci0, norb, nelec, link_index=None, **kwargs): - ci0 = _state_args (ci0) - link_index = _solver_args (link_index) + res = list(self._collect('spin_square', ci0, norb, nelec, *args, **kwargs)) + ss = [x[0] for x in res] + multip = [x[1] for x in res] + return ss, multip + methods['spin_square'] = spin_square + methods['states_spin_square'] = states_spin_square + else: + methods['spin_square'] = None + + if all(hasattr(solver, 'large_ci') for solver in fcisolvers): + def states_large_ci(self, fcivec, norb, nelec, *args, **kwargs): + _solver_args = self._solver_args + _state_args = self._state_args + fcivec = _state_args (fcivec) nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) - dm1a, dm1b = [], [] - dm2aa, dm2ab, dm2bb = [], [], [] - for dm1s, dm2s in self._collect ('make_rdm12s', ci0, norb, nelec, link_index=link_index, **kwargs): - dm1a.append(dm1s[0]) - dm1b.append(dm1s[1]) - dm2aa.append(dm2s[0]) - dm2ab.append(dm2s[1]) - dm2bb.append(dm2s[2]) - return (dm1a, dm1b), (dm2aa, dm2ab, dm2bb) - - def make_rdm12s(self, ci0, norb, nelec, link_index=None, **kwargs): - rdm1s, rdm2s = self.states_make_rdm12s(ci0, norb, nelec, link_index=link_index, **kwargs) - rdm1s = numpy.einsum ('r,srpq->spq', self.weights, rdm1s) - rdm2s = numpy.einsum ('r,srpqtu->spqtu', self.weights, rdm2s) - return rdm1s, rdm2s - - # TODO: linkstr support - def states_trans_rdm12 (self, ci1, ci0, norb, nelec, link_index=None, **kwargs): - ci1 = _state_args (ci1) - ci0 = _state_args (ci0) - link_index = _solver_args (link_index) + return list(self._collect('large_ci', fcivec, norb, nelec, *args, **kwargs)) + methods['states_large_ci'] = states_large_ci + else: + methods['large_ci'] = None + + if all(hasattr(solver, 'transform_ci_for_orbital_rotation') for solver in fcisolvers): + def states_transform_ci_for_orbital_rotation(self, fcivec, norb, nelec, + *args, **kwargs): + _solver_args = self._solver_args + _state_args = self._state_args + fcivec = _state_args (fcivec) nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) - tdm1 = [] - tdm2 = [] - for dm1, dm2 in self._collect ('trans_rdm12', ci1, ci0, norb, nelec, link_index=link_index, **kwargs): - tdm1.append (dm1) - tdm2.append (dm2) - return tdm1, tdm2 - - def trans_rdm12 (self, ci1, ci0, norb, nelec, link_index=None, **kwargs): - tdm1, tdm2 = self.states_trans_rdm12 (ci1, ci0, norb, nelec, link_index=link_index, **kwargs) - tdm1 = numpy.einsum ('r,rpq->pq', self.weights, tdm1) - tdm2 = numpy.einsum ('r,rpqst->pqst', self.weights, tdm2) - return tdm1, tdm2 - - if has_spin_square: - def spin_square(self, ci0, norb, nelec, *args, **kwargs): - ss, multip = self.states_spin_square(ci0, norb, nelec, *args, **kwargs) - weights = self.weights - return numpy.dot(ss, weights), numpy.dot(multip, weights) - - def states_spin_square(self, ci0, norb, nelec, *args, **kwargs): - ci0 = _state_args (ci0) - nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) - res = list(self._collect('spin_square', ci0, norb, nelec, *args, **kwargs)) - ss = [x[0] for x in res] - multip = [x[1] for x in res] - return ss, multip - else: - spin_square = None - - large_ci = None - if has_large_ci: - def states_large_ci(self, fcivec, norb, nelec, *args, **kwargs): - fcivec = _state_args (fcivec) - nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) - return list(self._collect('large_ci', fcivec, norb, nelec, *args, **kwargs)) - - transform_ci_for_orbital_rotation = None - if has_transform_ci: - def states_transform_ci_for_orbital_rotation(self, fcivec, norb, nelec, - *args, **kwargs): - fcivec = _state_args (fcivec) - nelec = _solver_args ([self._get_nelec (solver, nelec) for solver in self.fcisolvers]) - return list(self._collect('transform_ci_for_orbital_rotation', - fcivec, norb, nelec, *args, **kwargs)) - - fcisolver = FakeCISolver(casscf.mol) - fcisolver.__dict__.update(casscf.fcisolver.__dict__) - fcisolver.nroots = nroots + return list(self._collect('transform_ci_for_orbital_rotation', + fcivec, norb, nelec, *args, **kwargs)) + methods['states_transform_ci_for_orbital_rotation'] = states_transform_ci_for_orbital_rotation + else: + methods['transform_ci_for_orbital_rotation'] = None + + fcisolver = lib.set_class(StateAverageMixFCISolver(fcisolvers, weights), + (StateAverageMixFCISolver, fcisolvers[0].__class__), + attrs=methods) mc = _state_average_mcscf_solver(casscf, fcisolver) return mc @@ -1473,70 +1509,3 @@ def state_average_mix_(casscf, fcisolvers, weights=(0.5,0.5)): del (BASE, MAP2HF_TOL) - - -if __name__ == '__main__': - from pyscf import mcscf - from pyscf.tools import ring - - mol = gto.M(verbose=0, - output=None, - atom=[['H', c] for c in ring.make(6, 1.2)], - basis='6-31g') - - m = scf.RHF(mol) - ehf = m.scf() - - mc = mcscf.CASSCF(m, 6, 6) - mc.verbose = 4 - emc, e_ci, fcivec, mo, mo_energy = mc.mc1step() - print(ehf, emc, emc-ehf) - print(emc - -3.272089958) - - rdm1 = make_rdm1(mc, mo, fcivec) - rdm1, rdm2 = make_rdm12(mc, mo, fcivec) - print(rdm1) - - mo1 = cas_natorb(mc)[0] - numpy.set_printoptions(2) - print(reduce(numpy.dot, (mo1[:,:6].T, mol.intor('int1e_ovlp_sph'), - mo[:,:6]))) - -# state average - mol.atom = [ - ['O', ( 0., 0. , 0. )], - ['H', ( 0., -0.757, 0.587)], - ['H', ( 0., 0.757 , 0.587)],] - mol.basis = '6-31g' - mol.symmetry = 1 - mol.build() - - m = scf.RHF(mol) - ehf = m.scf() - - - mc = mcscf.CASSCF(m, 4, 4) - mc.verbose = 4 - mc.fcisolver = fci.solver(mol, False) # to mix the singlet and triplet - mc = state_average_(mc, (.64,.36)) - emc, e_ci, fcivec, mo, mo_energy = mc.mc1step()[:5] - mc = mcscf.CASCI(m, 4, 4) - emc = mc.casci(mo)[0] - print(ehf, emc, emc-ehf) - print(emc - -76.003352190262532) - - mc = mcscf.CASSCF(m, 4, 4) - mc.verbose = 4 - mc = state_average_(mc, (.64,.36)) - emc, e_ci, fcivec, mo, mo_energy = mc.mc1step()[:5] - mc = mcscf.CASCI(m, 4, 4) - emc = mc.casci(mo)[0] - print(ehf, emc, emc-ehf) - print(emc - -75.982520334896776) - - - mc = mcscf.CASSCF(m, 4, 4) - mc.verbose = 4 - mc = state_specific_(mc, 2) - emc = mc.kernel()[0] - diff --git a/pyscf/mcscf/mc1step.py b/pyscf/mcscf/mc1step.py index 9ce98e9317..59040526dd 100644 --- a/pyscf/mcscf/mc1step.py +++ b/pyscf/mcscf/mc1step.py @@ -1331,7 +1331,7 @@ def expmat(a): from pyscf import fci from pyscf.mcscf import addons - mol = Mole() + mol = gto.Mole() mol.verbose = 0 mol.output = None#"out_h2o" mol.atom = [ diff --git a/pyscf/mcscf/test/test_h2o.py b/pyscf/mcscf/test/test_h2o.py index f11375cbd9..970683c286 100644 --- a/pyscf/mcscf/test/test_h2o.py +++ b/pyscf/mcscf/test/test_h2o.py @@ -161,12 +161,27 @@ def test_casci(self): mc.fcisolver = fci.solver(mol) mc.natorb = 1 emc = mc.kernel()[0] - self.assertAlmostEqual(emc, -75.9624554777, 9) + self.assertAlmostEqual(emc, -75.9624554777, 7) - mc = CASCI(m, 4, (3,1)) + mc = mcscf.CASCI(m, 4, (3,1)) mc.fcisolver = fci.solver(mol, False) emc = mc.casci()[0] - self.assertAlmostEqual(emc, -75.439016172976, 9) + self.assertAlmostEqual(emc, -75.439016172976, 6) + + def test_addons(self): + mc = mcscf.CASSCF(msym, 4, 4) + mc.fcisolver = fci.solver(molsym, False) # to mix the singlet and triplet + mc = mc.state_average_((.64,.36)) + emc, e_ci, fcivec, mo, mo_energy = mc.mc1step()[:5] + self.assertAlmostEqual(emc, -75.85387884606675, 8) + mc = mcscf.CASCI(msym, 4, 4) + emc = mc.casci(mo)[0] + self.assertAlmostEqual(emc, -75.98341123168858, 8) + + mc = mcscf.CASSCF(msym, 4, 4) + mc = mc.state_specific_(2) + emc = mc.kernel()[0] + self.assertAlmostEqual(emc, -75.59353002290788, 8) if __name__ == "__main__": print("Full Tests for H2O")