diff --git a/examples/cc/20-ip_ea_eom_ccsd.py b/examples/cc/20-ip_ea_eom_ccsd.py index 5d43d9ac49..e8a3c8b6cb 100644 --- a/examples/cc/20-ip_ea_eom_ccsd.py +++ b/examples/cc/20-ip_ea_eom_ccsd.py @@ -1,7 +1,7 @@ #!/usr/bin/env python ''' -Ground-state and IP/EA-EOM-CCSD for singlet (RHF) and triplet (UHF) O2. +Ground-state, EOM-EE-GCCSD and IP/EA-EOM-CCSD for singlet (RHF) and triplet (UHF) O2. ''' from pyscf import gto, scf, cc @@ -49,3 +49,10 @@ eip,cip = mycc.ipccsd(nroots=1) eea,cea = mycc.eaccsd(nroots=1) +# EOM-GCCSD +mf = mf.to_ghf() +mycc = mf.GCCSD() +ecc, t1, t2 = mycc.kernel() +e,v = mycc.ipccsd(nroots=6) +e,v = mycc.eaccsd(nroots=6) +e,v = mycc.eeccsd(nroots=6) diff --git a/pyscf/cc/ccsd.py b/pyscf/cc/ccsd.py index 12a658039d..30284b51e6 100644 --- a/pyscf/cc/ccsd.py +++ b/pyscf/cc/ccsd.py @@ -51,10 +51,11 @@ def kernel(mycc, eris=None, t1=None, t2=None, max_cycle=50, tol=1e-8, elif t2 is None: t2 = mycc.get_init_guess(eris)[1] + name = mycc.__class__.__name__ cput1 = cput0 = (logger.process_clock(), logger.perf_counter()) eold = 0 eccsd = mycc.energy(t1, t2, eris) - log.info('Init E_corr(CCSD) = %.15g', eccsd) + log.info('Init E_corr(%s) = %.15g', name, eccsd) if isinstance(mycc.diis, lib.diis.DIIS): adiis = mycc.diis @@ -82,13 +83,13 @@ def kernel(mycc, eris=None, t1=None, t2=None, max_cycle=50, tol=1e-8, t1new = t2new = None t1, t2 = mycc.run_diis(t1, t2, istep, normt, eccsd-eold, adiis) eold, eccsd = eccsd, mycc.energy(t1, t2, eris) - log.info('cycle = %d E_corr(CCSD) = %.15g dE = %.9g norm(t1,t2) = %.6g', - istep+1, eccsd, eccsd - eold, normt) - cput1 = log.timer('CCSD iter', *cput1) + log.info('cycle = %d E_corr(%s) = %.15g dE = %.9g norm(t1,t2) = %.6g', + istep+1, name, eccsd, eccsd - eold, normt) + cput1 = log.timer(f'{name} iter', *cput1) if abs(eccsd-eold) < tol and normt < tolnormt: conv = True break - log.timer('CCSD', *cput0) + log.timer(name, *cput0) return conv, eccsd, t1, t2 @@ -713,7 +714,8 @@ def energy(mycc, t1=None, t2=None, eris=None): e += 2 * numpy.einsum('ijab,iabj', tau, eris_ovvo) e -= numpy.einsum('jiab,iabj', tau, eris_ovvo) if abs(e.imag) > 1e-4: - logger.warn(mycc, 'Non-zero imaginary part found in CCSD energy %s', e) + logger.warn(mycc, 'Non-zero imaginary part found in %s energy %s', + mycc.__class__.__name__, e) return e.real def restore_from_diis_(mycc, diis_file, inplace=True): @@ -841,7 +843,7 @@ def __call__(self, mol_or_geom, **kwargs): return self.e_tot -class CCSD(lib.StreamObject): +class CCSDBase(lib.StreamObject): '''restricted CCSD Attributes: @@ -919,6 +921,7 @@ class CCSD(lib.StreamObject): async_io = getattr(__config__, 'cc_ccsd_CCSD_async_io', True) incore_complete = getattr(__config__, 'cc_ccsd_CCSD_incore_complete', False) cc2 = getattr(__config__, 'cc_ccsd_CCSD_cc2', False) + callback = None _keys = set(( 'max_cycle', 'conv_tol', 'iterative_damping', @@ -927,7 +930,7 @@ class CCSD(lib.StreamObject): 'async_io', 'incore_complete', 'cc2', 'callback', 'mol', 'verbose', 'stdout', 'frozen', 'level_shift', 'mo_coeff', 'mo_occ', 'converged', 'converged_lambda', 'emp2', 'e_hf', - 'e_corr', 't1', 't2', 'l1', 'l2', 'chkfile', 'callback', + 'e_corr', 't1', 't2', 'l1', 'l2', 'chkfile', )) def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None): @@ -977,7 +980,6 @@ def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None): self._nocc = None self._nmo = None self.chkfile = mf.chkfile - self.callback = None @property def ecc(self): @@ -1035,12 +1037,6 @@ def dump_flags(self, verbose=None): log.info('diis_start_energy_diff = %g', self.diis_start_energy_diff) log.info('max_memory %d MB (current use %d MB)', self.max_memory, lib.current_memory()[0]) - if (log.verbose >= logger.DEBUG1 and - self.__class__ == CCSD): - nocc = self.nocc - nvir = self.nmo - self.nocc - flops = _flops(nocc, nvir) - log.debug1('total FLOPs %s', flops) return self def get_init_guess(self, eris=None): @@ -1113,78 +1109,45 @@ def _finalize(self): as_scanner = as_scanner restore_from_diis_ = restore_from_diis_ - - def solve_lambda(self, t1=None, t2=None, l1=None, l2=None, - eris=None): - from pyscf.cc import ccsd_lambda - if t1 is None: t1 = self.t1 - if t2 is None: t2 = self.t2 - if eris is None: eris = self.ao2mo(self.mo_coeff) - self.converged_lambda, self.l1, self.l2 = \ - ccsd_lambda.kernel(self, eris, t1, t2, l1, l2, - max_cycle=self.max_cycle, - tol=self.conv_tol_normt, - verbose=self.verbose) - return self.l1, self.l2 + def solve_lambda(self, t1=None, t2=None, l1=None, l2=None, eris=None): + raise NotImplementedError def ccsd_t(self, t1=None, t2=None, eris=None): - from pyscf.cc import ccsd_t - if t1 is None: t1 = self.t1 - if t2 is None: t2 = self.t2 - if eris is None: eris = self.ao2mo(self.mo_coeff) - return ccsd_t.kernel(self, eris, t1, t2, self.verbose) + raise NotImplementedError def ipccsd(self, nroots=1, left=False, koopmans=False, guess=None, partition=None, eris=None): - from pyscf.cc import eom_rccsd - return eom_rccsd.EOMIP(self).kernel(nroots, left, koopmans, guess, - partition, eris) + raise NotImplementedError def eaccsd(self, nroots=1, left=False, koopmans=False, guess=None, partition=None, eris=None): - from pyscf.cc import eom_rccsd - return eom_rccsd.EOMEA(self).kernel(nroots, left, koopmans, guess, - partition, eris) + raise NotImplementedError def eeccsd(self, nroots=1, koopmans=False, guess=None, eris=None): - from pyscf.cc import eom_rccsd - return eom_rccsd.EOMEE(self).kernel(nroots, koopmans, guess, eris) + raise NotImplementedError def eomee_ccsd_singlet(self, nroots=1, koopmans=False, guess=None, eris=None): - from pyscf.cc import eom_rccsd - return eom_rccsd.EOMEESinglet(self).kernel(nroots, koopmans, guess, eris) + raise NotImplementedError def eomee_ccsd_triplet(self, nroots=1, koopmans=False, guess=None, eris=None): - from pyscf.cc import eom_rccsd - return eom_rccsd.EOMEETriplet(self).kernel(nroots, koopmans, guess, eris) + raise NotImplementedError def eomsf_ccsd(self, nroots=1, koopmans=False, guess=None, eris=None): - from pyscf.cc import eom_rccsd - return eom_rccsd.EOMEESpinFlip(self).kernel(nroots, koopmans, guess, eris) + raise NotImplementedError def eomip_method(self): - from pyscf.cc import eom_rccsd - return eom_rccsd.EOMIP(self) + raise NotImplementedError def eomea_method(self): - from pyscf.cc import eom_rccsd - return eom_rccsd.EOMEA(self) + raise NotImplementedError def eomee_method(self): - from pyscf.cc import eom_rccsd - return eom_rccsd.EOMEE(self) + raise NotImplementedError def make_rdm1(self, t1=None, t2=None, l1=None, l2=None, ao_repr=False, with_frozen=True, with_mf=True): '''Un-relaxed 1-particle density matrix in MO space''' - from pyscf.cc import ccsd_rdm - if t1 is None: t1 = self.t1 - if t2 is None: t2 = self.t2 - if l1 is None: l1 = self.l1 - if l2 is None: l2 = self.l2 - if l1 is None: l1, l2 = self.solve_lambda(t1, t2) - return ccsd_rdm.make_rdm1(self, t1, t2, l1, l2, ao_repr=ao_repr, - with_frozen=with_frozen, with_mf=with_mf) + raise NotImplementedError def make_rdm2(self, t1=None, t2=None, l1=None, l2=None, ao_repr=False, with_frozen=True, with_dm1=True): @@ -1193,14 +1156,7 @@ def make_rdm2(self, t1=None, t2=None, l1=None, l2=None, ao_repr=False, dm2[p,r,q,s] = ''' - from pyscf.cc import ccsd_rdm - if t1 is None: t1 = self.t1 - if t2 is None: t2 = self.t2 - if l1 is None: l1 = self.l1 - if l2 is None: l2 = self.l2 - if l1 is None: l1, l2 = self.solve_lambda(t1, t2) - return ccsd_rdm.make_rdm2(self, t1, t2, l1, l2, ao_repr=ao_repr, - with_frozen=with_frozen, with_dm1=with_dm1) + raise NotImplementedError def ao2mo(self, mo_coeff=None): # Pseudo code how eris are implemented: @@ -1288,6 +1244,111 @@ def dump_chk(self, t1_t2=None, frozen=None, mo_coeff=None, mo_occ=None): lib.chkfile.save(self.chkfile, 'ccsd', cc_chk) + def density_fit(self, auxbasis=None, with_df=None): + raise NotImplementedError + + def nuc_grad_method(self): + raise NotImplementedError + +class CCSD(CCSDBase): + __doc__ = CCSDBase.__doc__ + + def dump_flags(self, verbose=None): + CCSDBase.dump_flags(self, verbose) + if self.verbose >= logger.DEBUG1 and self.__class__ == CCSD: + nocc = self.nocc + nvir = self.nmo - self.nocc + flops = _flops(nocc, nvir) + logger.debug1(self, 'total FLOPs %s', flops) + return self + + def solve_lambda(self, t1=None, t2=None, l1=None, l2=None, eris=None): + from pyscf.cc import ccsd_lambda + if t1 is None: t1 = self.t1 + if t2 is None: t2 = self.t2 + if eris is None: eris = self.ao2mo(self.mo_coeff) + self.converged_lambda, self.l1, self.l2 = \ + ccsd_lambda.kernel(self, eris, t1, t2, l1, l2, + max_cycle=self.max_cycle, + tol=self.conv_tol_normt, + verbose=self.verbose) + return self.l1, self.l2 + + def ccsd_t(self, t1=None, t2=None, eris=None): + from pyscf.cc import ccsd_t + if t1 is None: t1 = self.t1 + if t2 is None: t2 = self.t2 + if eris is None: eris = self.ao2mo(self.mo_coeff) + return ccsd_t.kernel(self, eris, t1, t2, self.verbose) + + def ipccsd(self, nroots=1, left=False, koopmans=False, guess=None, + partition=None, eris=None): + from pyscf.cc import eom_rccsd + return eom_rccsd.EOMIP(self).kernel(nroots, left, koopmans, guess, + partition, eris) + + def eaccsd(self, nroots=1, left=False, koopmans=False, guess=None, + partition=None, eris=None): + from pyscf.cc import eom_rccsd + return eom_rccsd.EOMEA(self).kernel(nroots, left, koopmans, guess, + partition, eris) + + def eeccsd(self, nroots=1, koopmans=False, guess=None, eris=None): + from pyscf.cc import eom_rccsd + return eom_rccsd.EOMEE(self).kernel(nroots, koopmans, guess, eris) + + def eomee_ccsd_singlet(self, nroots=1, koopmans=False, guess=None, eris=None): + from pyscf.cc import eom_rccsd + return eom_rccsd.EOMEESinglet(self).kernel(nroots, koopmans, guess, eris) + + def eomee_ccsd_triplet(self, nroots=1, koopmans=False, guess=None, eris=None): + from pyscf.cc import eom_rccsd + return eom_rccsd.EOMEETriplet(self).kernel(nroots, koopmans, guess, eris) + + def eomsf_ccsd(self, nroots=1, koopmans=False, guess=None, eris=None): + from pyscf.cc import eom_rccsd + return eom_rccsd.EOMEESpinFlip(self).kernel(nroots, koopmans, guess, eris) + + def eomip_method(self): + from pyscf.cc import eom_rccsd + return eom_rccsd.EOMIP(self) + + def eomea_method(self): + from pyscf.cc import eom_rccsd + return eom_rccsd.EOMEA(self) + + def eomee_method(self): + from pyscf.cc import eom_rccsd + return eom_rccsd.EOMEE(self) + + def make_rdm1(self, t1=None, t2=None, l1=None, l2=None, ao_repr=False, + with_frozen=True, with_mf=True): + '''Un-relaxed 1-particle density matrix in MO space''' + from pyscf.cc import ccsd_rdm + if t1 is None: t1 = self.t1 + if t2 is None: t2 = self.t2 + if l1 is None: l1 = self.l1 + if l2 is None: l2 = self.l2 + if l1 is None: l1, l2 = self.solve_lambda(t1, t2) + return ccsd_rdm.make_rdm1(self, t1, t2, l1, l2, ao_repr=ao_repr, + with_frozen=with_frozen, with_mf=with_mf) + + def make_rdm2(self, t1=None, t2=None, l1=None, l2=None, ao_repr=False, + with_frozen=True, with_dm1=True): + '''2-particle density matrix in MO space. The density matrix is + stored as + + dm2[p,r,q,s] = + ''' + from pyscf.cc import ccsd_rdm + if t1 is None: t1 = self.t1 + if t2 is None: t2 = self.t2 + if l1 is None: l1 = self.l1 + if l2 is None: l2 = self.l2 + if l1 is None: l1, l2 = self.solve_lambda(t1, t2) + return ccsd_rdm.make_rdm2(self, t1, t2, l1, l2, ao_repr=ao_repr, + with_frozen=with_frozen, with_dm1=with_dm1) + def density_fit(self, auxbasis=None, with_df=None): from pyscf.cc import dfccsd mycc = dfccsd.RCCSD(self._scf, self.frozen, self.mo_coeff, self.mo_occ) @@ -1315,7 +1376,6 @@ def get_d2_diagnostic(self, t2=None): if t2 is None: t2 = self.t2 return get_d2_diagnostic(t2) - CC = RCCSD = CCSD diff --git a/pyscf/cc/ccsd_lambda.py b/pyscf/cc/ccsd_lambda.py index f55513e235..dd23b18636 100644 --- a/pyscf/cc/ccsd_lambda.py +++ b/pyscf/cc/ccsd_lambda.py @@ -57,7 +57,8 @@ def kernel(mycc, eris=None, t1=None, t2=None, l1=None, l2=None, adiis.space = mycc.diis_space else: adiis = None - cput0 = log.timer('CCSD lambda initialization', *cput0) + name = mycc.__class__.__name__ + cput0 = log.timer(f'{name} lambda initialization', *cput0) conv = False for istep in range(max_cycle): @@ -68,7 +69,7 @@ def kernel(mycc, eris=None, t1=None, t2=None, l1=None, l2=None, l1new = l2new = None l1, l2 = mycc.run_diis(l1, l2, istep, normt, 0, adiis) log.info('cycle = %d norm(lambda1,lambda2) = %.6g', istep+1, normt) - cput0 = log.timer('CCSD iter', *cput0) + cput0 = log.timer(f'{name} iter', *cput0) if normt < tol: conv = True break diff --git a/pyscf/cc/ccsd_t.py b/pyscf/cc/ccsd_t.py index 8009adb0fa..52e0b972c8 100644 --- a/pyscf/cc/ccsd_t.py +++ b/pyscf/cc/ccsd_t.py @@ -37,6 +37,7 @@ def kernel(mycc, eris, t1=None, t2=None, verbose=logger.NOTE): if t1 is None: t1 = mycc.t1 if t2 is None: t2 = mycc.t2 + name = mycc.__class__.__name__ nocc, nvir = t1.shape nmo = nocc + nvir @@ -52,7 +53,7 @@ def kernel(mycc, eris, t1=None, t2=None, verbose=logger.NOTE): mo_energy, t1T, t2T, vooo, fvo, restore_t2_inplace = \ _sort_t2_vooo_(mycc, orbsym, t1, t2, eris) - cpu1 = log.timer_debug1('CCSD(T) sort_eri', *cpu1) + cpu1 = log.timer_debug1(f'{name}(T) sort_eri', *cpu1) cpu2 = list(cpu1) orbsym = numpy.hstack((numpy.sort(orbsym[:nocc]),numpy.sort(orbsym[nocc:]))) @@ -124,11 +125,11 @@ def contract(a0, a1, b0, b1, cache): t2 = restore_t2_inplace(t2T) et_sum *= 2 if abs(et_sum[0].imag) > 1e-4: - logger.warn(mycc, 'Non-zero imaginary part of CCSD(T) energy was found %s', - et_sum[0]) + logger.warn(mycc, 'Non-zero imaginary part of %s(T) energy was found %s', + name, et_sum[0]) et = et_sum[0].real - log.timer('CCSD(T)', *cpu0) - log.note('CCSD(T) correction = %.15g', et) + log.timer(f'{name}(T)', *cpu0) + log.note('%s(T) correction = %.15g', name, et) return et def _sort_eri(mycc, eris, nocc, nvir, vvop, log): diff --git a/pyscf/cc/ccsd_t_rdm_slow.py b/pyscf/cc/ccsd_t_rdm_slow.py index df75d18c82..efe56a0f11 100644 --- a/pyscf/cc/ccsd_t_rdm_slow.py +++ b/pyscf/cc/ccsd_t_rdm_slow.py @@ -125,48 +125,3 @@ def r6(w): return (4 * w + w.transpose(0,1,2,4,5,3) + w.transpose(0,1,2,5,3,4) - 2 * w.transpose(0,1,2,5,4,3) - 2 * w.transpose(0,1,2,3,5,4) - 2 * w.transpose(0,1,2,4,3,5)) - - -if __name__ == '__main__': - from functools import reduce - from pyscf import gto - from pyscf import scf - from pyscf.cc import ccsd - from pyscf.cc import ccsd_t_slow as ccsd_t - from pyscf.cc import ccsd_t_lambda_slow as ccsd_t_lambda - from pyscf import ao2mo - - mol = gto.Mole() - #mol.verbose = 5 - #mol.output = 'out_h2o' - mol.atom = [ - [8 , (0. , 0. , 0.)], - [1 , (0. , -.957 , .587)], - [1 , (0.2, .757 , .487)]] - - #mol.basis = 'ccpvdz' - mol.basis = '631g' - mol.build() - mf = scf.RHF(mol) - mf.conv_tol = 1e-1 - mf.scf() - mcc = ccsd.CCSD(mf) - mcc.conv_tol = 1e-14 - ecc, t1, t2 = mcc.kernel() - eris = mcc.ao2mo() - e3ref = ccsd_t.kernel(mcc, eris, t1, t2) - l1, l2 = ccsd_t_lambda.kernel(mcc, eris, t1, t2)[1:] - print(ecc, e3ref) - - eri_mo = ao2mo.kernel(mf._eri, mf.mo_coeff, compact=False) - nmo = mf.mo_coeff.shape[1] - eri_mo = eri_mo.reshape(nmo,nmo,nmo,nmo) - dm1 = make_rdm1(mcc, t1, t2, l1, l2, eris=eris) - dm2 = make_rdm2(mcc, t1, t2, l1, l2, eris=eris) - print(lib.finger(dm1) - 1.289951975176953 ) - print(lib.finger(dm2) - 6.6184784979411164) - h1 = reduce(numpy.dot, (mf.mo_coeff.T, mf.get_hcore(), mf.mo_coeff)) - e3 =(numpy.einsum('ij,ji->', h1, dm1) + - numpy.einsum('ijkl,ijkl->', eri_mo, dm2)*.5 + mf.mol.energy_nuc()) - print(e3ref, e3-(mf.e_tot+ecc)) - diff --git a/pyscf/cc/gccsd.py b/pyscf/cc/gccsd.py index 0b66143c19..f9b4687727 100644 --- a/pyscf/cc/gccsd.py +++ b/pyscf/cc/gccsd.py @@ -111,14 +111,14 @@ def amplitudes_from_rccsd(t1, t2, orbspin=None): return spatial2spin(t1, orbspin), spatial2spin(t2, orbspin) -class GCCSD(ccsd.CCSD): +class GCCSD(ccsd.CCSDBase): conv_tol = getattr(__config__, 'cc_gccsd_GCCSD_conv_tol', 1e-7) conv_tol_normt = getattr(__config__, 'cc_gccsd_GCCSD_conv_tol_normt', 1e-6) def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None): assert (isinstance(mf, scf.ghf.GHF)) - ccsd.CCSD.__init__(self, mf, frozen, mo_coeff, mo_occ) + ccsd.CCSDBase.__init__(self, mf, frozen, mo_coeff, mo_occ) def init_amps(self, eris=None): if eris is None: @@ -160,7 +160,7 @@ def ccsd(self, t1=None, t2=None, eris=None, mbpt2=False): if not np.any(orbspin == -1): self.mo_coeff = lib.tag_array(self.mo_coeff, orbspin=orbspin) - e_corr, self.t1, self.t2 = ccsd.CCSD.ccsd(self, t1, t2, eris) + e_corr, self.t1, self.t2 = ccsd.CCSDBase.ccsd(self, t1, t2, eris) if getattr(eris, 'orbspin', None) is not None: self.t1 = lib.tag_array(self.t1, orbspin=eris.orbspin) self.t2 = lib.tag_array(self.t2, orbspin=eris.orbspin) @@ -290,27 +290,6 @@ def spin2spatial(self, tx, orbspin=None): orbspin = orbspin[self.get_frozen_mask()] return spin2spatial(tx, orbspin) - def density_fit(self): - raise NotImplementedError - - def nuc_grad_method(self): - raise NotImplementedError - - def get_t1_diagnostic(self, t1=None): - if t1 is None: t1 = self.t1 - raise NotImplementedError - #return get_t1_diagnostic(t1) - - def get_d1_diagnostic(self, t1=None): - if t1 is None: t1 = self.t1 - raise NotImplementedError - #return get_d1_diagnostic(t1) - - def get_d2_diagnostic(self, t2=None): - if t2 is None: t2 = self.t2 - raise NotImplementedError - #return get_d2_diagnostic(self.spin2spatial(t2)) - CCSD = GCCSD @@ -565,51 +544,3 @@ def _make_eris_outcore(mycc, mo_coeff=None): cput0 = log.timer_debug1('transforming vvvv', *cput0) return eris - - -if __name__ == '__main__': - from pyscf import gto - mol = gto.Mole() - mol.atom = [['O', (0., 0., 0.)], - ['O', (1.21, 0., 0.)]] - mol.basis = 'cc-pvdz' - mol.spin = 2 - mol.build() - mf = scf.UHF(mol).run() - mf = scf.addons.convert_to_ghf(mf) - - # Freeze 1s electrons - frozen = [0,1,2,3] - gcc = GCCSD(mf, frozen=frozen) - ecc, t1, t2 = gcc.kernel() - print(ecc - -0.3486987472235819) - - mol = gto.Mole() - mol.atom = [ - [8 , (0. , 0. , 0.)], - [1 , (0. , -0.757 , 0.587)], - [1 , (0. , 0.757 , 0.587)]] - mol.basis = 'cc-pvdz' - mol.spin = 0 - mol.build() - mf = scf.UHF(mol).run() - mf = scf.addons.convert_to_ghf(mf) - - mycc = GCCSD(mf) - ecc, t1, t2 = mycc.kernel() - print(ecc - -0.2133432712431435) - e,v = mycc.ipccsd(nroots=8) - print(e[0] - 0.4335604332073799) - print(e[2] - 0.5187659896045407) - print(e[4] - 0.6782876002229172) - - e,v = mycc.eaccsd(nroots=8) - print(e[0] - 0.16737886338859731) - print(e[2] - 0.24027613852009164) - print(e[4] - 0.51006797826488071) - - e,v = mycc.eeccsd(nroots=4) - print(e[0] - 0.2757159395886167) - print(e[1] - 0.2757159395886167) - print(e[2] - 0.2757159395886167) - print(e[3] - 0.3005716731825082) diff --git a/pyscf/cc/gccsd_t_lambda.py b/pyscf/cc/gccsd_t_lambda.py index 492f23d6cc..cb5a8b20f6 100644 --- a/pyscf/cc/gccsd_t_lambda.py +++ b/pyscf/cc/gccsd_t_lambda.py @@ -81,80 +81,3 @@ def update_lambda(mycc, t1, t2, l1, l2, eris=None, imds=None): l1 += imds.l1_t l2 += imds.l2_t return l1, l2 - - -if __name__ == '__main__': - from pyscf import gto - from pyscf import scf - from pyscf import cc - - mol = gto.Mole() - mol.atom = [ - [8 , (0. , 0. , 0.)], - [1 , (0. , -0.757 , 0.587)], - [1 , (0. , 0.757 , 0.587)]] - mol.basis = '631g' - mol.spin = 0 - mol.build() - mf0 = mf = scf.RHF(mol).run(conv_tol=1.) - mf = scf.addons.convert_to_ghf(mf) - - from pyscf.cc import ccsd_t_lambda_slow as ccsd_t_lambda - mycc0 = cc.CCSD(mf0) - eris0 = mycc0.ao2mo() - mycc0.kernel(eris=eris0) - t1 = mycc0.t1 - t2 = mycc0.t2 - imds = ccsd_t_lambda.make_intermediates(mycc0, t1, t2, eris0) - l1, l2 = ccsd_t_lambda.update_lambda(mycc0, t1, t2, t1, t2, eris0, imds) - l1ref, l2ref = ccsd_t_lambda.update_lambda(mycc0, t1, t2, l1, l2, eris0, imds) - mycc = cc.GCCSD(mf) - eris = mycc.ao2mo() - t1 = mycc.spatial2spin(t1, mycc.mo_coeff.orbspin) - t2 = mycc.spatial2spin(t2, mycc.mo_coeff.orbspin) - l1 = mycc.spatial2spin(l1, mycc.mo_coeff.orbspin) - l2 = mycc.spatial2spin(l2, mycc.mo_coeff.orbspin) - imds = make_intermediates(mycc, t1, t2, eris) - l1, l2 = update_lambda(mycc, t1, t2, l1, l2, eris, imds) - l1 = mycc.spin2spatial(l1, mycc.mo_coeff.orbspin) - l2 = mycc.spin2spatial(l2, mycc.mo_coeff.orbspin) - print(abs(l2[1]-l2[1].transpose(1,0,2,3)-l2[0]).max()) - print(abs(l2[1]-l2[1].transpose(0,1,3,2)-l2[0]).max()) - print(abs(l1[0]-l1ref).max()) - print(abs(l2[1]-l2ref).max()) - - mol = gto.Mole() - mol.atom = [ - [8 , (0. , 0. , 0.)], - [1 , (0. , -0.757 , 0.587)], - [1 , (0. , 0.757 , 0.587)]] - mol.basis = '631g' - mol.spin = 2 - mol.build() - mf0 = mf = scf.UHF(mol).run(conv_tol=1) - mf = scf.addons.convert_to_ghf(mf) - - from pyscf.cc import uccsd_t_lambda - mycc0 = cc.CCSD(mf0) - eris0 = mycc0.ao2mo() - mycc0.kernel(eris=eris0) - t1 = mycc0.t1 - t2 = mycc0.t2 - imds = uccsd_t_lambda.make_intermediates(mycc0, t1, t2, eris0) - l1, l2 = uccsd_t_lambda.update_lambda(mycc0, t1, t2, t1, t2, eris0, imds) - l1ref, l2ref = uccsd_t_lambda.update_lambda(mycc0, t1, t2, l1, l2, eris0, imds) - mycc = cc.GCCSD(mf) - eris = mycc.ao2mo() - t1 = mycc.spatial2spin(t1, mycc.mo_coeff.orbspin) - t2 = mycc.spatial2spin(t2, mycc.mo_coeff.orbspin) - l1 = mycc.spatial2spin(l1, mycc.mo_coeff.orbspin) - l2 = mycc.spatial2spin(l2, mycc.mo_coeff.orbspin) - imds = make_intermediates(mycc, t1, t2, eris) - l1, l2 = update_lambda(mycc, t1, t2, l1, l2, eris, imds) - l1 = mycc.spin2spatial(l1, mycc.mo_coeff.orbspin) - l2 = mycc.spin2spatial(l2, mycc.mo_coeff.orbspin) - print(abs(l1[0]-l1ref[0]).max()) - print(abs(l1[1]-l1ref[1]).max()) - print(abs(l2[0]-l2ref[0]).max()) - print(abs(l2[1]-l2ref[1]).max()) - print(abs(l2[2]-l2ref[2]).max()) diff --git a/pyscf/cc/qcisd.py b/pyscf/cc/qcisd.py index 765eca7713..becdeda3a5 100644 --- a/pyscf/cc/qcisd.py +++ b/pyscf/cc/qcisd.py @@ -28,7 +28,7 @@ from pyscf import gto from pyscf import lib from pyscf.lib import logger -from pyscf.cc.ccsd import CCSD, _add_vvvv, _flops, _ChemistsERIs +from pyscf.cc.ccsd import CCSDBase, as_scanner, _add_vvvv, _ChemistsERIs from pyscf import __config__ BLKMIN = getattr(__config__, 'cc_ccsd_blkmin', 4) @@ -305,7 +305,7 @@ def load_ovvv(buf, p0): return fswap['wVOov'], fswap['wVooV'] def energy(mycc, t1=None, t2=None, eris=None): - '''CCSD correlation energy''' + '''QCISD correlation energy''' if t1 is None: t1 = mycc.t1 if t2 is None: t2 = mycc.t2 if eris is None: eris = mycc.ao2mo() @@ -324,49 +324,8 @@ def energy(mycc, t1=None, t2=None, eris=None): logger.warn(mycc, 'Non-zero imaginary part found in QCISD energy %s', e) return e.real -def as_scanner(cc): - '''Generating a scanner/solver for QCISD PES. - The returned solver is a function. This function requires one argument - "mol" as input and returns total QCISD energy. - - ''' - if isinstance(cc, lib.SinglePointScanner): - return cc - - logger.info(cc, 'Set %s as a scanner', cc.__class__) - name = cc.__class__.__name__ + QCISD_Scanner.__name_mixin__ - return lib.set_class(QCISD_Scanner(cc), (QCISD_Scanner, cc.__class__), name) - -class QCISD_Scanner(lib.SinglePointScanner): - def __init__(self, cc): - self.__dict__.update(cc.__dict__) - self._scf = cc._scf.as_scanner() - - def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, gto.Mole): - mol = mol_or_geom - else: - mol = self.mol.set_geom_(mol_or_geom, inplace=False) - - if self.t2 is not None: - last_size = self.vector_size() - else: - last_size = 0 - - self.reset(mol) - - mf_scanner = self._scf - mf_scanner(mol) - self.mo_coeff = mf_scanner.mo_coeff - self.mo_occ = mf_scanner.mo_occ - if last_size != self.vector_size(): - self.t1 = self.t2 = None - self.kernel(self.t1, self.t2, **kwargs) - return self.e_tot - - -class QCISD(CCSD): +class QCISD(CCSDBase): '''restricted QCISD ''' @@ -387,40 +346,13 @@ def dump_flags(self, verbose=None): log.info('diis_start_energy_diff = %g', self.diis_start_energy_diff) log.info('max_memory %d MB (current use %d MB)', self.max_memory, lib.current_memory()[0]) - if (log.verbose >= logger.DEBUG1 and - self.__class__ == QCISD): - nocc = self.nocc - nvir = self.nmo - self.nocc - flops = _flops(nocc, nvir) - log.debug1('total FLOPs %s', flops) return self energy = energy _add_vvvv = _add_vvvv update_amps = update_amps - - def kernel(self, t1=None, t2=None, eris=None): - return self.qcisd(t1, t2, eris) - def qcisd(self, t1=None, t2=None, eris=None): - assert (self.mo_coeff is not None) - assert (self.mo_occ is not None) - - if self.verbose >= logger.WARN: - self.check_sanity() - self.dump_flags() - - self.e_hf = self.get_e_hf() - - if eris is None: - eris = self.ao2mo(self.mo_coeff) - - self.converged, self.e_corr, self.t1, self.t2 = \ - kernel(self, eris, t1, t2, max_cycle=self.max_cycle, - tol=self.conv_tol, tolnormt=self.conv_tol_normt, - verbose=self.verbose) - self._finalize() - return self.e_corr, self.t1, self.t2 - + kernel = CCSDBase.ccsd + ccsd = lib.invalid_method('ccsd') as_scanner = as_scanner def qcisd_t(self, t1=None, t2=None, eris=None): @@ -430,44 +362,4 @@ def qcisd_t(self, t1=None, t2=None, eris=None): if eris is None: eris = self.ao2mo(self.mo_coeff) return qcisd_t.kernel(self, eris, t1, t2, self.verbose) - RQCISD = QCISD - - -if __name__ == '__main__': - from pyscf import scf - - mol = gto.Mole() - #mol.atom = [ - # [8 , (0. , 0. , 0.)], - # [1 , (0. , -0.757 , 0.587)], - # [1 , (0. , 0.757 , 0.587)]] - mol.atom = [['Ne', (0,0,0)]] - mol.basis = 'cc-pvdz' - mol.verbose = 7 - mol.spin = 0 - mol.build() - mf = scf.RHF(mol).run(conv_tol=1e-14) - - mycc = QCISD(mf, frozen=1) - ecc, t1, t2 = mycc.kernel() - et = mycc.qcisd_t() - print("QCISD(T) =", mycc.e_tot+et) - - mol = gto.Mole() - mol.atom = """C 0.000 0.000 0.000 - H 0.637 0.637 0.637 - H -0.637 -0.637 0.637 - H -0.637 0.637 -0.637 - H 0.637 -0.637 -0.637""" - mol.basis = 'cc-pvdz' - mol.verbose = 7 - mol.spin = 0 - mol.build() - mf = scf.RHF(mol).run(conv_tol=1e-14) - - mycc = QCISD(mf, frozen=1) - ecc, t1, t2 = mycc.kernel() - print(mycc.e_tot - -40.383989) - et = mycc.qcisd_t() - print(mycc.e_tot+et - -40.387679) diff --git a/pyscf/cc/rccsd.py b/pyscf/cc/rccsd.py index f7ed87f980..6ed9bb7e64 100644 --- a/pyscf/cc/rccsd.py +++ b/pyscf/cc/rccsd.py @@ -186,7 +186,7 @@ def ccsd(self, t1=None, t2=None, eris=None, mbpt2=False): if eris is None: eris = self.ao2mo(self.mo_coeff) - return ccsd.CCSD.ccsd(self, t1, t2, eris) + return ccsd.CCSDBase.ccsd(self, t1, t2, eris) def ao2mo(self, mo_coeff=None): nmo = self.nmo @@ -226,12 +226,6 @@ def solve_lambda(self, t1=None, t2=None, l1=None, l2=None, verbose=self.verbose) return self.l1, self.l2 - def ccsd_t(self, t1=None, t2=None, eris=None): - return ccsd.CCSD.ccsd_t(self, t1, t2, eris) - - def density_fit(self, auxbasis=None, with_df=None): - raise NotImplementedError - class _ChemistsERIs(ccsd._ChemistsERIs): diff --git a/pyscf/cc/test/test_ccd.py b/pyscf/cc/test/test_ccd.py new file mode 100644 index 0000000000..0bd9d7d417 --- /dev/null +++ b/pyscf/cc/test/test_ccd.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python +# Copyright 2023 The PySCF Developers. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import unittest +import numpy as np +from pyscf import gto +from pyscf import scf +from pyscf import ao2mo +from pyscf.cc import ccd + +def setUpModule(): + global mol, mf, mycc + mol = gto.Mole() + mol.verbose = 7 + mol.output = '/dev/null' + mol.atom = [ + [8 , (0. , 0. , 0.)], + [1 , (0. , -0.757 , 0.587)], + [1 , (0. , 0.757 , 0.587)]] + + mol.basis = '631g' + mol.build() + mf = scf.RHF(mol) + mf.kernel() + mycc = ccd.CCD(mf) + mycc.kernel() + +def tearDownModule(): + global mol, mf, mycc + mol.stdout.close() + del mol, mf, mycc + + +class KnownValues(unittest.TestCase): + def test_ccd(self): + self.assertAlmostEqual(mycc.e_corr, -0.134712806, 7) + + def test_rdm(self): + dm1 = mycc.make_rdm1() + dm2 = mycc.make_rdm2() + h1 = mf.mo_coeff.T.dot(mf.get_hcore()).dot(mf.mo_coeff) + nmo = mf.mo_coeff.shape[1] + eri = ao2mo.restore(1, ao2mo.kernel(mf._eri, mf.mo_coeff), nmo) + e1 = np.einsum('ij,ji', h1, dm1) + e1+= np.einsum('ijkl,ijkl', eri, dm2) * .5 + e1+= mol.energy_nuc() + self.assertAlmostEqual(e1, mycc.e_tot, 6) + +if __name__ == "__main__": + print("Full Tests for CCD") + unittest.main() diff --git a/pyscf/cc/test/test_ccsd_t.py b/pyscf/cc/test/test_ccsd_t.py index c39b9d401a..79fe33e2b1 100644 --- a/pyscf/cc/test/test_ccsd_t.py +++ b/pyscf/cc/test/test_ccsd_t.py @@ -19,8 +19,11 @@ from pyscf import gto, scf, lib, symm from pyscf import cc +from pyscf import ao2mo from pyscf.cc import ccsd_t from pyscf.cc import gccsd, gccsd_t +from pyscf.cc import ccsd_t_lambda_slow as ccsd_t_lambda +from pyscf.cc import ccsd_t_rdm_slow as ccsd_t_rdm def setUpModule(): global mol, rhf, mcc @@ -191,6 +194,34 @@ def test_ccsd_t_complex(self): self.assertAlmostEqual(e0, e1.real, 9) self.assertAlmostEqual(e1, -0.98756910139720788-0.0019567929592079489j, 9) + def test_ccsd_t_rdm(self): + mol = gto.Mole() + mol.atom = [ + [8 , (0. , 0. , 0.)], + [1 , (0. , -.957 , .587)], + [1 , (0.2, .757 , .487)]] + + mol.basis = '631g' + mol.build() + mf = scf.RHF(mol) + mf.conv_tol = 1e-1 + mf.scf() + mcc = mf.CCSD() + ecc, t1, t2 = mcc.kernel() + eris = mcc.ao2mo() + e3ref = ccsd_t.kernel(mcc, eris, t1, t2) + l1, l2 = ccsd_t_lambda.kernel(mcc, eris, t1, t2)[1:] + + eri_mo = ao2mo.kernel(mf._eri, mf.mo_coeff, compact=False) + nmo = mf.mo_coeff.shape[1] + eri_mo = eri_mo.reshape(nmo,nmo,nmo,nmo) + dm1 = ccsd_t_rdm.make_rdm1(mcc, t1, t2, l1, l2, eris=eris) + dm2 = ccsd_t_rdm.make_rdm2(mcc, t1, t2, l1, l2, eris=eris) + h1 = reduce(numpy.dot, (mf.mo_coeff.T, mf.get_hcore(), mf.mo_coeff)) + e3 =(numpy.einsum('ij,ji->', h1, dm1) + + numpy.einsum('ijkl,ijkl->', eri_mo, dm2)*.5 + mf.mol.energy_nuc()) + self.assertAlmostEqual(e3ref, e3-(mf.e_tot+ecc), 7) + if __name__ == "__main__": print("Full Tests for CCSD(T)") diff --git a/pyscf/cc/test/test_gccsd.py b/pyscf/cc/test/test_gccsd.py index 24e3be681f..393e8bbced 100644 --- a/pyscf/cc/test/test_gccsd.py +++ b/pyscf/cc/test/test_gccsd.py @@ -54,6 +54,22 @@ class KnownValues(unittest.TestCase): def test_gccsd(self): self.assertAlmostEqual(gcc1.e_corr, -0.10805861695870976, 6) + def test_frozen(self): + mol = gto.Mole() + mol.atom = [['O', (0., 0., 0.)], + ['O', (1.21, 0., 0.)]] + mol.basis = 'cc-pvdz' + mol.spin = 2 + mol.build() + mf = scf.UHF(mol).run() + mf = scf.addons.convert_to_ghf(mf) + + # Freeze 1s electrons + frozen = [0,1,2,3] + gcc = gccsd.GCCSD(mf, frozen=frozen) + ecc, t1, t2 = gcc.kernel() + self.assertAlmostEqual(ecc, -0.3486987472235819, 6) + def test_ERIS(self): gcc = gccsd.GCCSD(mf, frozen=4) numpy.random.seed(9) diff --git a/pyscf/cc/test/test_gccsd_lambda.py b/pyscf/cc/test/test_gccsd_lambda.py index 5b74aca33c..6b0a55b205 100644 --- a/pyscf/cc/test/test_gccsd_lambda.py +++ b/pyscf/cc/test/test_gccsd_lambda.py @@ -25,6 +25,9 @@ from pyscf import cc from pyscf.cc import addons from pyscf.cc import gccsd_lambda +from pyscf.cc import ccsd_t_lambda_slow as ccsd_t_lambda +from pyscf.cc import gccsd_t_lambda +from pyscf.cc import uccsd_t_lambda def setUpModule(): global mol @@ -242,6 +245,76 @@ def test_update_amps(self): self.assertAlmostEqual(abs(l1-l1ref).max(), 0, 8) self.assertAlmostEqual(abs(l2-l2ref).max(), 0, 8) + def test_gccsd_t_lambda(self): + mol = gto.Mole() + mol.atom = [ + [8 , (0. , 0. , 0.)], + [1 , (0. , -0.757 , 0.587)], + [1 , (0. , 0.757 , 0.587)]] + mol.basis = '631g' + mol.spin = 0 + mol.build() + mf0 = mf = scf.RHF(mol).run(conv_tol=1.) + mf = scf.addons.convert_to_ghf(mf) + + mycc0 = cc.CCSD(mf0) + eris0 = mycc0.ao2mo() + mycc0.kernel(eris=eris0) + t1 = mycc0.t1 + t2 = mycc0.t2 + imds = ccsd_t_lambda.make_intermediates(mycc0, t1, t2, eris0) + l1, l2 = ccsd_t_lambda.update_lambda(mycc0, t1, t2, t1, t2, eris0, imds) + l1ref, l2ref = ccsd_t_lambda.update_lambda(mycc0, t1, t2, l1, l2, eris0, imds) + mycc = cc.GCCSD(mf) + eris = mycc.ao2mo() + t1 = mycc.spatial2spin(t1, mycc.mo_coeff.orbspin) + t2 = mycc.spatial2spin(t2, mycc.mo_coeff.orbspin) + l1 = mycc.spatial2spin(l1, mycc.mo_coeff.orbspin) + l2 = mycc.spatial2spin(l2, mycc.mo_coeff.orbspin) + imds = gccsd_t_lambda.make_intermediates(mycc, t1, t2, eris) + l1, l2 = gccsd_t_lambda.update_lambda(mycc, t1, t2, l1, l2, eris, imds) + l1 = mycc.spin2spatial(l1, mycc.mo_coeff.orbspin) + l2 = mycc.spin2spatial(l2, mycc.mo_coeff.orbspin) + self.assertAlmostEqual(abs(l2[1]-l2[1].transpose(1,0,2,3)-l2[0]).max(), 0, 12) + self.assertAlmostEqual(abs(l2[1]-l2[1].transpose(0,1,3,2)-l2[0]).max(), 0, 12) + self.assertAlmostEqual(abs(l1[0]-l1ref).max(), 0, 12) + self.assertAlmostEqual(abs(l2[1]-l2ref).max(), 0, 12) + + mol = gto.Mole() + mol.atom = [ + [8 , (0. , 0. , 0.)], + [1 , (0. , -0.757 , 0.587)], + [1 , (0. , 0.757 , 0.587)]] + mol.basis = '631g' + mol.spin = 2 + mol.build() + mf0 = mf = scf.UHF(mol).run(conv_tol=1) + mf = scf.addons.convert_to_ghf(mf) + + mycc0 = cc.CCSD(mf0) + eris0 = mycc0.ao2mo() + mycc0.kernel(eris=eris0) + t1 = mycc0.t1 + t2 = mycc0.t2 + imds = uccsd_t_lambda.make_intermediates(mycc0, t1, t2, eris0) + l1, l2 = uccsd_t_lambda.update_lambda(mycc0, t1, t2, t1, t2, eris0, imds) + l1ref, l2ref = uccsd_t_lambda.update_lambda(mycc0, t1, t2, l1, l2, eris0, imds) + mycc = cc.GCCSD(mf) + eris = mycc.ao2mo() + t1 = mycc.spatial2spin(t1, mycc.mo_coeff.orbspin) + t2 = mycc.spatial2spin(t2, mycc.mo_coeff.orbspin) + l1 = mycc.spatial2spin(l1, mycc.mo_coeff.orbspin) + l2 = mycc.spatial2spin(l2, mycc.mo_coeff.orbspin) + imds = gccsd_t_lambda.make_intermediates(mycc, t1, t2, eris) + l1, l2 = gccsd_t_lambda.update_lambda(mycc, t1, t2, l1, l2, eris, imds) + l1 = mycc.spin2spatial(l1, mycc.mo_coeff.orbspin) + l2 = mycc.spin2spatial(l2, mycc.mo_coeff.orbspin) + self.assertAlmostEqual(abs(l1[0]-l1ref[0]).max(), 0, 12) + self.assertAlmostEqual(abs(l1[1]-l1ref[1]).max(), 0, 12) + self.assertAlmostEqual(abs(l2[0]-l2ref[0]).max(), 0, 12) + self.assertAlmostEqual(abs(l2[1]-l2ref[1]).max(), 0, 12) + self.assertAlmostEqual(abs(l2[2]-l2ref[2]).max(), 0, 12) + if __name__ == "__main__": print("Full Tests for GCCSD lambda") diff --git a/pyscf/cc/test/test_qcisd.py b/pyscf/cc/test/test_qcisd.py new file mode 100644 index 0000000000..92514de318 --- /dev/null +++ b/pyscf/cc/test/test_qcisd.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python +# Copyright 2023 The PySCF Developers. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import unittest +import numpy as np +from pyscf import gto +from pyscf import scf +from pyscf import ao2mo +from pyscf.cc import qcisd + +def setUpModule(): + global mol, mf, mycc + mol = gto.Mole() + mol.verbose = 7 + mol.output = '/dev/null' + mol.atom = [ + [8 , (0. , 0. , 0.)], + [1 , (0. , -0.757 , 0.587)], + [1 , (0. , 0.757 , 0.587)]] + + mol.basis = '631g' + mol.build() + mf = scf.RHF(mol) + mf.kernel() + mycc = qcisd.QCISD(mf) + mycc.kernel() + +def tearDownModule(): + global mol, mf, mycc + mol.stdout.close() + del mol, mf, mycc + + +class KnownValues(unittest.TestCase): + def test_qcisd_t(self): + mol = gto.Mole() + mol.atom = """C 0.000 0.000 0.000 + H 0.637 0.637 0.637 + H -0.637 -0.637 0.637 + H -0.637 0.637 -0.637 + H 0.637 -0.637 -0.637""" + mol.basis = 'cc-pvdz' + mol.verbose = 7 + mol.output = '/dev/null' + mol.spin = 0 + mol.build() + mf = scf.RHF(mol).run() + + mycc = qcisd.QCISD(mf, frozen=1) + ecc, t1, t2 = mycc.kernel() + self.assertAlmostEqual(mycc.e_tot, -40.3839884, 6) + et = mycc.qcisd_t() + self.assertAlmostEqual(mycc.e_tot+et, -40.38767969, 6) + + def test_qcisd_t_frozen(self): + mol = gto.Mole() + mol.atom = [['Ne', (0,0,0)]] + mol.basis = 'cc-pvdz' + mol.verbose = 7 + mol.output = '/dev/null' + mol.spin = 0 + mol.build() + mf = scf.RHF(mol).run() + + mycc = qcisd.QCISD(mf, frozen=1).as_scanner() + mycc(mol) + et = mycc.qcisd_t() + self.assertAlmostEqual(mycc.e_tot+et, -128.6788843055109, 6) + +if __name__ == "__main__": + print("Full Tests for CCD") + unittest.main() diff --git a/pyscf/cc/test/test_uccsd_t.py b/pyscf/cc/test/test_uccsd_t.py index f5eb091b43..b7b7452c1c 100644 --- a/pyscf/cc/test/test_uccsd_t.py +++ b/pyscf/cc/test/test_uccsd_t.py @@ -19,9 +19,16 @@ from functools import reduce from pyscf import gto, scf, lib, symm +from pyscf import ao2mo from pyscf import cc from pyscf.cc import uccsd_t from pyscf.cc import gccsd, gccsd_t +from pyscf.cc import uccsd_t_slow +from pyscf.cc import uccsd_t_lambda +from pyscf.cc import uccsd_t_rdm +from pyscf.cc import ccsd_t_lambda_slow as ccsd_t_lambda +from pyscf.cc import ccsd_t_rdm_slow as ccsd_t_rdm + def setUpModule(): global mol, mol1, mf, mcc @@ -165,8 +172,63 @@ def test_uccsd_t_complex(self): self.assertAlmostEqual(e0, e1.real, 9) self.assertAlmostEqual(e1, -0.056092415718338388-0.011390417704868244j, 9) + def test_uccsd_t_rdm(self): + mol = gto.Mole() + mol.atom = [ + [8 , (0. , 0. , 0.)], + [1 , (0. , -.957 , .587)], + [1 , (0.2, .757 , .487)]] + mol.basis = '631g' + mol.build() + mf0 = mf = scf.RHF(mol).run(conv_tol=1.) + mf = scf.addons.convert_to_uhf(mf) + + mycc0 = cc.CCSD(mf0) + eris0 = mycc0.ao2mo() + mycc0.kernel(eris=eris0) + t1 = mycc0.t1 + t2 = mycc0.t2 + imds = ccsd_t_lambda.make_intermediates(mycc0, t1, t2, eris0) + l1, l2 = ccsd_t_lambda.update_lambda(mycc0, t1, t2, t1, t2, eris0, imds) + dm1ref = ccsd_t_rdm.make_rdm1(mycc0, t1, t2, l1, l2, eris0) + dm2ref = ccsd_t_rdm.make_rdm2(mycc0, t1, t2, l1, l2, eris0) + + t1 = (t1, t1) + t2aa = t2 - t2.transpose(1,0,2,3) + t2 = (t2aa, t2, t2aa) + l1 = (l1, l1) + l2aa = l2 - l2.transpose(1,0,2,3) + l2 = (l2aa, l2, l2aa) + mycc = cc.UCCSD(mf) + eris = mycc.ao2mo() + dm1 = uccsd_t_rdm.make_rdm1(mycc, t1, t2, l1, l2, eris) + dm2 = uccsd_t_rdm.make_rdm2(mycc, t1, t2, l1, l2, eris) + trdm1 = dm1[0] + dm1[1] + trdm2 = dm2[0] + dm2[1] + dm2[1].transpose(2,3,0,1) + dm2[2] + self.assertAlmostEqual(abs(trdm1 - dm1ref).max(), 0, 12) + self.assertAlmostEqual(abs(trdm2 - dm2ref).max(), 0, 12) + + ecc = mycc.kernel(eris=eris)[0] + l1, l2 = mycc.solve_lambda(eris=eris) + e3ref = mycc.e_tot + mycc.ccsd_t() + + nmoa, nmob = mycc.nmo + eri_aa = ao2mo.kernel(mf._eri, mf.mo_coeff[0], compact=False).reshape([nmoa]*4) + eri_bb = ao2mo.kernel(mf._eri, mf.mo_coeff[1], compact=False).reshape([nmob]*4) + eri_ab = ao2mo.kernel(mf._eri, [mf.mo_coeff[k] for k in (0,0,1,1)], + compact=False).reshape(nmoa,nmoa,nmob,nmob) + dm1 = uccsd_t_rdm.make_rdm1(mycc, t1, t2, l1, l2, eris=eris) + dm2 = uccsd_t_rdm.make_rdm2(mycc, t1, t2, l1, l2, eris=eris) + h1a = reduce(numpy.dot, (mf.mo_coeff[0].T, mf.get_hcore(), mf.mo_coeff[0])) + h1b = reduce(numpy.dot, (mf.mo_coeff[1].T, mf.get_hcore(), mf.mo_coeff[1])) + e3 =(numpy.einsum('ij,ji->', h1a, dm1[0]) + + numpy.einsum('ij,ji->', h1b, dm1[1]) + + numpy.einsum('ijkl,ijkl->', eri_aa, dm2[0])*.5 + + numpy.einsum('ijkl,ijkl->', eri_bb, dm2[2])*.5 + + numpy.einsum('ijkl,ijkl->', eri_ab, dm2[1]) + + mf.mol.energy_nuc()) + self.assertAlmostEqual(e3, e3ref, 9) if __name__ == "__main__": print("Full Tests for UCCSD(T)") unittest.main() - diff --git a/pyscf/cc/uccsd.py b/pyscf/cc/uccsd.py index 115ffec35d..b7ff3a57e3 100644 --- a/pyscf/cc/uccsd.py +++ b/pyscf/cc/uccsd.py @@ -536,7 +536,7 @@ def _add_vvvv(mycc, t1, t2, eris, out=None, with_ovvv=False, t2sym=None): return u2aa,u2ab,u2bb -class UCCSD(ccsd.CCSD): +class UCCSD(ccsd.CCSDBase): conv_tol = getattr(__config__, 'cc_uccsd_UCCSD_conv_tol', 1e-7) conv_tol_normt = getattr(__config__, 'cc_uccsd_UCCSD_conv_tol_normt', 1e-6) @@ -548,7 +548,7 @@ class UCCSD(ccsd.CCSD): # orbitals, second list is for beta orbitals def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None): assert isinstance(mf, scf.uhf.UHF) - ccsd.CCSD.__init__(self, mf, frozen, mo_coeff, mo_occ) + ccsd.CCSDBase.__init__(self, mf, frozen, mo_coeff, mo_occ) get_nocc = get_nocc get_nmo = get_nmo @@ -612,7 +612,7 @@ def ccsd(self, t1=None, t2=None, eris=None, mbpt2=False): self.t1 = (np.zeros((nocca,nvira)), np.zeros((noccb,nvirb))) return self.e_corr, self.t1, self.t2 - return ccsd.CCSD.ccsd(self, t1, t2, eris) + return ccsd.CCSDBase.ccsd(self, t1, t2, eris) def solve_lambda(self, t1=None, t2=None, l1=None, l2=None, eris=None): @@ -733,9 +733,6 @@ def eomee_method(self): from pyscf.cc import eom_uccsd return eom_uccsd.EOMEE(self) - def density_fit(self): - raise NotImplementedError - def nuc_grad_method(self): from pyscf.grad import uccsd return uccsd.Gradients(self) @@ -762,21 +759,6 @@ def vector_size(self, nmo=None, nocc=None): def amplitudes_from_rccsd(self, t1, t2): return amplitudes_from_rccsd(t1, t2) - def get_t1_diagnostic(self, t1=None): - if t1 is None: t1 = self.t1 - raise NotImplementedError - #return get_t1_diagnostic(t1) - - def get_d1_diagnostic(self, t1=None): - if t1 is None: t1 = self.t1 - raise NotImplementedError - #return get_d1_diagnostic(t1) - - def get_d2_diagnostic(self, t2=None): - if t2 is None: t2 = self.t2 - raise NotImplementedError - #return get_d2_diagnostic(t2) - CCSD = UCCSD diff --git a/pyscf/cc/uccsd_t_rdm.py b/pyscf/cc/uccsd_t_rdm.py index 90c64138bf..50c21cf471 100644 --- a/pyscf/cc/uccsd_t_rdm.py +++ b/pyscf/cc/uccsd_t_rdm.py @@ -292,71 +292,3 @@ def r4(w): w = w - w.transpose(0,2,1,3,4,5) w = w + w.transpose(0,2,1,3,5,4) return w - -if __name__ == '__main__': - from functools import reduce - from pyscf import gto - from pyscf import scf - from pyscf import ao2mo - from pyscf import cc - from pyscf.cc import uccsd_t_slow - from pyscf.cc import uccsd_t_lambda - - mol = gto.Mole() - mol.atom = [ - [8 , (0. , 0. , 0.)], - [1 , (0. , -.957 , .587)], - [1 , (0.2, .757 , .487)]] - mol.basis = '631g' - mol.build() - mf0 = mf = scf.RHF(mol).run(conv_tol=1.) - mf = scf.addons.convert_to_uhf(mf) - - from pyscf.cc import ccsd_t_lambda_slow as ccsd_t_lambda - from pyscf.cc import ccsd_t_rdm_slow as ccsd_t_rdm - mycc0 = cc.CCSD(mf0) - eris0 = mycc0.ao2mo() - mycc0.kernel(eris=eris0) - t1 = mycc0.t1 - t2 = mycc0.t2 - imds = ccsd_t_lambda.make_intermediates(mycc0, t1, t2, eris0) - l1, l2 = ccsd_t_lambda.update_lambda(mycc0, t1, t2, t1, t2, eris0, imds) - dm1ref = ccsd_t_rdm.make_rdm1(mycc0, t1, t2, l1, l2, eris0) - dm2ref = ccsd_t_rdm.make_rdm2(mycc0, t1, t2, l1, l2, eris0) - - t1 = (t1, t1) - t2aa = t2 - t2.transpose(1,0,2,3) - t2 = (t2aa, t2, t2aa) - l1 = (l1, l1) - l2aa = l2 - l2.transpose(1,0,2,3) - l2 = (l2aa, l2, l2aa) - mycc = cc.UCCSD(mf) - eris = mycc.ao2mo() - dm1 = make_rdm1(mycc, t1, t2, l1, l2, eris) - dm2 = make_rdm2(mycc, t1, t2, l1, l2, eris) - trdm1 = dm1[0] + dm1[1] - trdm2 = dm2[0] + dm2[1] + dm2[1].transpose(2,3,0,1) + dm2[2] - print(abs(trdm1 - dm1ref).max()) - print(abs(trdm2 - dm2ref).max()) - - ecc = mycc.kernel(eris=eris)[0] - l1, l2 = mycc.solve_lambda(eris=eris) - e3ref = mycc.e_tot + mycc.ccsd_t() - - nmoa, nmob = mycc.nmo - eri_aa = ao2mo.kernel(mf._eri, mf.mo_coeff[0], compact=False).reshape([nmoa]*4) - eri_bb = ao2mo.kernel(mf._eri, mf.mo_coeff[1], compact=False).reshape([nmob]*4) - eri_ab = ao2mo.kernel(mf._eri, [mf.mo_coeff[k] for k in (0,0,1,1)], - compact=False).reshape(nmoa,nmoa,nmob,nmob) - dm1 = make_rdm1(mycc, t1, t2, l1, l2, eris=eris) - dm2 = make_rdm2(mycc, t1, t2, l1, l2, eris=eris) - h1a = reduce(numpy.dot, (mf.mo_coeff[0].T, mf.get_hcore(), mf.mo_coeff[0])) - h1b = reduce(numpy.dot, (mf.mo_coeff[1].T, mf.get_hcore(), mf.mo_coeff[1])) - e3 =(numpy.einsum('ij,ji->', h1a, dm1[0]) + - numpy.einsum('ij,ji->', h1b, dm1[1]) + - numpy.einsum('ijkl,ijkl->', eri_aa, dm2[0])*.5 + - numpy.einsum('ijkl,ijkl->', eri_bb, dm2[2])*.5 + - numpy.einsum('ijkl,ijkl->', eri_ab, dm2[1]) + - mf.mol.energy_nuc()) - print(e3 - e3ref) - diff --git a/pyscf/ci/cisd.py b/pyscf/ci/cisd.py index a0c57c9a15..64875eeb35 100644 --- a/pyscf/ci/cisd.py +++ b/pyscf/ci/cisd.py @@ -860,7 +860,7 @@ class CISD(lib.StreamObject): direct = getattr(__config__, 'ci_cisd_CISD_direct', False) async_io = getattr(__config__, 'ci_cisd_CISD_async_io', True) - keys = set(( + _keys = set(( 'conv_tol', 'max_cycle', 'max_space', 'lindep', 'level_shift', 'direct', 'async_io', 'mol', 'max_memory', 'nroots', 'frozen', 'chkfile', 'converged', 'mo_coeff', 'mo_occ', diff --git a/pyscf/mcscf/casci.py b/pyscf/mcscf/casci.py index 510d16be8e..0bbea9f61d 100644 --- a/pyscf/mcscf/casci.py +++ b/pyscf/mcscf/casci.py @@ -927,11 +927,10 @@ def ao2mo(self, mo_coeff=None): ''' raise NotImplementedError - get_h1cas = h1e_for_cas = h1e_for_cas + def get_h1cas(self, mo_coeff=None, ncas=None, ncore=None): + raise DeprecationWarning - def get_h1eff(self, mo_coeff=None, ncas=None, ncore=None): - return self.h1e_for_cas(mo_coeff, ncas, ncore) - get_h1eff.__doc__ = h1e_for_cas.__doc__ + get_h1eff = h1e_for_cas def casci(self, mo_coeff=None, ci0=None, verbose=None): raise NotImplementedError @@ -1181,70 +1180,3 @@ def nuc_grad_method(self): scf.uhf.UHF.CASCI = None del (WITH_META_LOWDIN, LARGE_CI_TOL, PENALTY) - - -if __name__ == '__main__': - from pyscf import mcscf - mol = gto.Mole() - mol.verbose = 0 - mol.output = None#"out_h2o" - mol.atom = [ - ['O', ( 0., 0. , 0. )], - ['H', ( 0., -0.757, 0.587)], - ['H', ( 0., 0.757 , 0.587)],] - - mol.basis = {'H': 'sto-3g', - 'O': '6-31g',} - mol.build() - - m = scf.RHF(mol) - ehf = m.scf() - mc = mcscf.CASCI(m, 4, 4) - mc.fcisolver = fci.solver(mol) - mc.natorb = 1 - emc = mc.kernel()[0] - print(ehf, emc, emc-ehf) - #-75.9577817425 -75.9624554777 -0.00467373522233 - print(emc+75.9624554777) - -# mc = CASCI(m, 4, (3,1)) -# #mc.fcisolver = fci.direct_spin1 -# mc.fcisolver = fci.solver(mol, False) -# emc = mc.casci()[0] -# print(emc - -75.439016172976) -# -# mol = gto.Mole() -# mol.verbose = 0 -# mol.output = "out_casci" -# mol.atom = [ -# ["C", (-0.65830719, 0.61123287, -0.00800148)], -# ["C", ( 0.73685281, 0.61123287, -0.00800148)], -# ["C", ( 1.43439081, 1.81898387, -0.00800148)], -# ["C", ( 0.73673681, 3.02749287, -0.00920048)], -# ["C", (-0.65808819, 3.02741487, -0.00967948)], -# ["C", (-1.35568919, 1.81920887, -0.00868348)], -# ["H", (-1.20806619, -0.34108413, -0.00755148)], -# ["H", ( 1.28636081, -0.34128013, -0.00668648)], -# ["H", ( 2.53407081, 1.81906387, -0.00736748)], -# ["H", ( 1.28693681, 3.97963587, -0.00925948)], -# ["H", (-1.20821019, 3.97969587, -0.01063248)], -# ["H", (-2.45529319, 1.81939187, -0.00886348)],] -# -# mol.basis = {'H': 'sto-3g', -# 'C': 'sto-3g',} -# mol.build() -# -# m = scf.RHF(mol) -# ehf = m.scf() -# mc = CASCI(m, 9, 8) -# mc.fcisolver = fci.solver(mol) -# emc = mc.casci()[0] -# print(ehf, emc, emc-ehf) -# print(emc - -227.948912536) -# -# mc = CASCI(m, 9, (5,3)) -# #mc.fcisolver = fci.direct_spin1 -# mc.fcisolver = fci.solver(mol, False) -# mc.fcisolver.nroots = 3 -# emc = mc.casci()[0] -# print(emc[0] - -227.7674519720) diff --git a/pyscf/mcscf/test/test_h2o.py b/pyscf/mcscf/test/test_h2o.py index fe22680521..f11375cbd9 100644 --- a/pyscf/mcscf/test/test_h2o.py +++ b/pyscf/mcscf/test/test_h2o.py @@ -144,6 +144,30 @@ def test_spin_and_pointgroup_sa4_newton (self): for e1, e0 in zip (numpy.sort (mc.e_states), mc_ref.e_states): self.assertAlmostEqual (e1, e0, 5) + def test_casci(self): + mol = gto.Mole() + mol.verbose = 0 + mol.atom = [ + ['O', ( 0., 0. , 0. )], + ['H', ( 0., -0.757, 0.587)], + ['H', ( 0., 0.757 , 0.587)],] + + mol.basis = {'H': 'sto-3g', + 'O': '6-31g',} + mol.build() + + m = scf.RHF(mol).run() + mc = mcscf.CASCI(m, 4, 4) + mc.fcisolver = fci.solver(mol) + mc.natorb = 1 + emc = mc.kernel()[0] + self.assertAlmostEqual(emc, -75.9624554777, 9) + + mc = CASCI(m, 4, (3,1)) + mc.fcisolver = fci.solver(mol, False) + emc = mc.casci()[0] + self.assertAlmostEqual(emc, -75.439016172976, 9) + if __name__ == "__main__": print("Full Tests for H2O") unittest.main() diff --git a/pyscf/mcscf/test/test_n2_df.py b/pyscf/mcscf/test/test_n2_df.py index 7d7ab6ed98..376920c3cb 100644 --- a/pyscf/mcscf/test/test_n2_df.py +++ b/pyscf/mcscf/test/test_n2_df.py @@ -192,11 +192,11 @@ def test_init(self): from pyscf.mcscf import df mf = scf.RHF(mol) self.assertTrue(isinstance(mcscf.CASCI(mf, 2, 2), mcscf.casci.CASCI)) - self.assertTrue(isinstance(mcscf.CASCI(mf.density_fit(), 2, 2), df._DFCASSCF)) + self.assertTrue(isinstance(mcscf.CASCI(mf.density_fit(), 2, 2), df._DFCASCI)) self.assertTrue(isinstance(mcscf.CASCI(mf.newton(), 2, 2), mcscf.casci.CASCI)) - self.assertTrue(isinstance(mcscf.CASCI(mf.density_fit().newton(), 2, 2), df._DFCASSCF)) + self.assertTrue(isinstance(mcscf.CASCI(mf.density_fit().newton(), 2, 2), df._DFCASCI)) self.assertTrue(isinstance(mcscf.CASCI(mf.newton().density_fit(), 2, 2), mcscf.casci.CASCI)) - self.assertTrue(isinstance(mcscf.CASCI(mf.density_fit().newton().density_fit(), 2, 2), df._DFCASSCF)) + self.assertTrue(isinstance(mcscf.CASCI(mf.density_fit().newton().density_fit(), 2, 2), df._DFCASCI)) self.assertTrue(isinstance(mcscf.CASSCF(mf, 2, 2), mcscf.mc1step.CASSCF)) self.assertTrue(isinstance(mcscf.CASSCF(mf.density_fit(), 2, 2), df._DFCASSCF)) @@ -205,12 +205,12 @@ def test_init(self): self.assertTrue(isinstance(mcscf.CASSCF(mf.newton().density_fit(), 2, 2), mcscf.mc1step.CASSCF)) self.assertTrue(isinstance(mcscf.CASSCF(mf.density_fit().newton().density_fit(), 2, 2), df._DFCASSCF)) - self.assertTrue(isinstance(mcscf.DFCASCI(mf, 2, 2), df._DFCASSCF)) - self.assertTrue(isinstance(mcscf.DFCASCI(mf.density_fit(), 2, 2), df._DFCASSCF)) - self.assertTrue(isinstance(mcscf.DFCASCI(mf.newton(), 2, 2), df._DFCASSCF)) - self.assertTrue(isinstance(mcscf.DFCASCI(mf.density_fit().newton(), 2, 2), df._DFCASSCF)) - self.assertTrue(isinstance(mcscf.DFCASCI(mf.newton().density_fit(), 2, 2), df._DFCASSCF)) - self.assertTrue(isinstance(mcscf.DFCASCI(mf.density_fit().newton().density_fit(), 2, 2), df._DFCASSCF)) + self.assertTrue(isinstance(mcscf.DFCASCI(mf, 2, 2), df._DFCASCI)) + self.assertTrue(isinstance(mcscf.DFCASCI(mf.density_fit(), 2, 2), df._DFCASCI)) + self.assertTrue(isinstance(mcscf.DFCASCI(mf.newton(), 2, 2), df._DFCASCI)) + self.assertTrue(isinstance(mcscf.DFCASCI(mf.density_fit().newton(), 2, 2), df._DFCASCI)) + self.assertTrue(isinstance(mcscf.DFCASCI(mf.newton().density_fit(), 2, 2), df._DFCASCI)) + self.assertTrue(isinstance(mcscf.DFCASCI(mf.density_fit().newton().density_fit(), 2, 2), df._DFCASCI)) self.assertTrue(isinstance(mcscf.DFCASSCF(mf, 2, 2), df._DFCASSCF)) self.assertTrue(isinstance(mcscf.DFCASSCF(mf.density_fit(), 2, 2), df._DFCASSCF)) @@ -220,11 +220,11 @@ def test_init(self): self.assertTrue(isinstance(mcscf.DFCASSCF(mf.density_fit().newton().density_fit(), 2, 2), df._DFCASSCF)) self.assertTrue(isinstance(mcscf.CASCI(msym, 2, 2), mcscf.casci_symm.CASCI)) - self.assertTrue(isinstance(mcscf.CASCI(msym.density_fit(), 2, 2), df._DFCASSCF)) + self.assertTrue(isinstance(mcscf.CASCI(msym.density_fit(), 2, 2), df._DFCASCI)) self.assertTrue(isinstance(mcscf.CASCI(msym.newton(), 2, 2), mcscf.casci_symm.CASCI)) - self.assertTrue(isinstance(mcscf.CASCI(msym.density_fit().newton(), 2, 2), df._DFCASSCF)) + self.assertTrue(isinstance(mcscf.CASCI(msym.density_fit().newton(), 2, 2), df._DFCASCI)) self.assertTrue(isinstance(mcscf.CASCI(msym.newton().density_fit(), 2, 2), mcscf.casci_symm.CASCI)) - self.assertTrue(isinstance(mcscf.CASCI(msym.density_fit().newton().density_fit(), 2, 2), df._DFCASSCF)) + self.assertTrue(isinstance(mcscf.CASCI(msym.density_fit().newton().density_fit(), 2, 2), df._DFCASCI)) self.assertTrue(isinstance(mcscf.CASSCF(msym, 2, 2), mcscf.mc1step_symm.CASSCF)) self.assertTrue(isinstance(mcscf.CASSCF(msym.density_fit(), 2, 2), df._DFCASSCF)) @@ -244,4 +244,3 @@ def test_init(self): if __name__ == "__main__": print("Full Tests for density fitting N2") unittest.main() - diff --git a/pyscf/mcscf/ucasci.py b/pyscf/mcscf/ucasci.py index 8ba75ca1b8..7143bbb3c4 100644 --- a/pyscf/mcscf/ucasci.py +++ b/pyscf/mcscf/ucasci.py @@ -32,7 +32,7 @@ from pyscf import scf from pyscf import fci from pyscf.mcscf import addons -from pyscf.mcscf import casci +from pyscf.mcscf.casci import as_scanner, CASBase from pyscf import __config__ WITH_META_LOWDIN = getattr(__config__, 'mcscf_analyze_with_meta_lowdin', True) @@ -227,11 +227,7 @@ def get_veff(self, mol=None, dm=None, hermi=1): def _eig(self, h, *args): return scf.hf.eig(h, None) - get_h1cas = h1e_for_cas = h1e_for_cas - - def get_h1eff(self, mo_coeff=None, ncas=None, ncore=None): - return self.h1e_for_cas(mo_coeff, ncas, ncore) - get_h1eff.__doc__ = h1e_for_cas.__doc__ + get_h1eff = h1e_for_cas def _finalize(self): log = logger.Logger(self.stdout, self.verbose) @@ -443,6 +439,8 @@ def kernel(self, mo_coeff=None, ci0=None): self._finalize() return self.e_tot, self.e_cas, self.ci + as_scanner = as_scanner + CASCI = UCASCI del (WITH_META_LOWDIN, LARGE_CI_TOL) diff --git a/pyscf/mrpt/nevpt2.py b/pyscf/mrpt/nevpt2.py index 4f4f138577..484406c65e 100644 --- a/pyscf/mrpt/nevpt2.py +++ b/pyscf/mrpt/nevpt2.py @@ -618,7 +618,7 @@ class NEVPT(lib.StreamObject): _keys = set(( 'ncore', 'root', 'compressed_mps', 'e_corr', 'canonicalized', 'onerdm', - )).union(casci.CASCI._keys, mc1step.CASSCF._keys) + )).union(casci.CASBase._keys, mc1step.CASSCF._keys) def __init__(self, mc, root=0): self.__dict__.update(mc.__dict__) @@ -849,7 +849,7 @@ def sc_nevpt(mc, ci=None, verbose=None): # register NEVPT2 in MCSCF -casci.CASCI.NEVPT2 = NEVPT +casci.CASBase.NEVPT2 = NEVPT diff --git a/pyscf/pbc/cc/__init__.py b/pyscf/pbc/cc/__init__.py index 684b7c77b0..614cba498f 100644 --- a/pyscf/pbc/cc/__init__.py +++ b/pyscf/pbc/cc/__init__.py @@ -38,20 +38,17 @@ def GCCSD(mf, frozen=None, mo_coeff=None, mo_occ=None): return ccsd.GCCSD(mf, frozen, mo_coeff, mo_occ) def KGCCSD(mf, frozen=None, mo_coeff=None, mo_occ=None): - from pyscf.pbc.cc import kccsd mf = scf.addons.convert_to_ghf(mf) - return kccsd.GCCSD(mf, frozen, mo_coeff, mo_occ) + return kgccsd.GCCSD(mf, frozen, mo_coeff, mo_occ) def KRCCSD(mf, frozen=None, mo_coeff=None, mo_occ=None): - from pyscf.pbc.cc import kccsd_rhf if not isinstance(mf, scf.khf.KRHF): mf = scf.addons.convert_to_rhf(mf) - return kccsd_rhf.RCCSD(mf, frozen, mo_coeff, mo_occ) + return krccsd.RCCSD(mf, frozen, mo_coeff, mo_occ) KCCSD = KRCCSD def KUCCSD(mf, frozen=None, mo_coeff=None, mo_occ=None): - from pyscf.pbc.cc import kccsd_uhf if not isinstance(mf, scf.kuhf.KUHF): mf = scf.addons.convert_to_uhf(mf) - return kccsd_uhf.UCCSD(mf, frozen, mo_coeff, mo_occ) + return kuccsd.UCCSD(mf, frozen, mo_coeff, mo_occ) diff --git a/pyscf/pbc/cc/kccsd_rhf.py b/pyscf/pbc/cc/kccsd_rhf.py index 659bd7a824..692be9d98b 100644 --- a/pyscf/pbc/cc/kccsd_rhf.py +++ b/pyscf/pbc/cc/kccsd_rhf.py @@ -505,10 +505,8 @@ class RCCSD(pyscf.cc.ccsd.CCSD): def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None): assert (isinstance(mf, scf.khf.KSCF)) - if isinstance(mf, scf.khf_ksymm.KsymAdaptedKSCF): - # mf.to_khf converts mf to a non-symmetry object - mf = mf.to_khf() - pyscf.cc.ccsd.CCSD.__init__(self, mf, frozen, mo_coeff, mo_occ) + # mf.to_khf converts mf to a non-symmetry object + pyscf.cc.ccsd.CCSD.__init__(self, mf.to_khf(), frozen, mo_coeff, mo_occ) self.kpts = mf.kpts self.khelper = kpts_helper.KptsHelper(mf.cell, mf.kpts, init_symm_map=False) diff --git a/pyscf/pbc/cc/kccsd_rhf_ksymm.py b/pyscf/pbc/cc/kccsd_rhf_ksymm.py index acd6d57683..b83f2bdfbc 100644 --- a/pyscf/pbc/cc/kccsd_rhf_ksymm.py +++ b/pyscf/pbc/cc/kccsd_rhf_ksymm.py @@ -395,7 +395,7 @@ def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None): ''' # NOTE self._scf is a non-symmetry object, see RCCSD.__init__ RCCSD.__init__(self, mf, frozen, mo_coeff, mo_occ) - self.kqrts = KQuartets(self.kpts).build() + self.kqrts = KQuartets(mf.kpts).build() self.rmat = None self.ktensor_direct = False self.eris_outcore = False diff --git a/pyscf/pbc/cc/kccsd_uhf.py b/pyscf/pbc/cc/kccsd_uhf.py index 1c460bdc6f..1c6acb8a22 100644 --- a/pyscf/pbc/cc/kccsd_uhf.py +++ b/pyscf/pbc/cc/kccsd_uhf.py @@ -1112,183 +1112,3 @@ def _make_df_eris(cc, mo_coeff=None): scf.kuhf.KUHF.CCSD = lib.class_as_method(KUCCSD) - - -if __name__ == '__main__': - from pyscf.pbc import gto - from pyscf import lo - - cell = gto.Cell() - cell.atom=''' - He 0.000000000000 0.000000000000 0.000000000000 - He 1.685068664391 1.685068664391 1.685068664391 - ''' - #cell.basis = [[0, (1., 1.)], [1, (.5, 1.)]] - cell.basis = [[0, (1., 1.)], [0, (.5, 1.)]] - cell.a = ''' - 0.000000000, 3.370137329, 3.370137329 - 3.370137329, 0.000000000, 3.370137329 - 3.370137329, 3.370137329, 0.000000000''' - cell.unit = 'B' - cell.mesh = [13]*3 - cell.build() - - np.random.seed(2) - # Running HF and CCSD with 1x1x2 Monkhorst-Pack k-point mesh - kmf = scf.KUHF(cell, kpts=cell.make_kpts([1,1,3]), exxdiv=None) - nmo = cell.nao_nr() - kmf.mo_occ = np.zeros((2,3,nmo)) - kmf.mo_occ[0,:,:3] = 1 - kmf.mo_occ[1,:,:1] = 1 - kmf.mo_energy = np.arange(nmo) + np.random.random((2,3,nmo)) * .3 - kmf.mo_energy[kmf.mo_occ == 0] += 2 - - mo = (np.random.random((2,3,nmo,nmo)) + - np.random.random((2,3,nmo,nmo))*1j - .5-.5j) - s = kmf.get_ovlp() - kmf.mo_coeff = np.empty_like(mo) - nkpts = len(kmf.kpts) - for k in range(nkpts): - kmf.mo_coeff[0,k] = lo.orth.vec_lowdin(mo[0,k], s[k]) - kmf.mo_coeff[1,k] = lo.orth.vec_lowdin(mo[1,k], s[k]) - - def rand_t1_t2(mycc): - nkpts = mycc.nkpts - nocca, noccb = mycc.nocc - nmoa, nmob = mycc.nmo - nvira, nvirb = nmoa - nocca, nmob - noccb - np.random.seed(1) - t1a = (np.random.random((nkpts,nocca,nvira)) + - np.random.random((nkpts,nocca,nvira))*1j - .5-.5j) - t1b = (np.random.random((nkpts,noccb,nvirb)) + - np.random.random((nkpts,noccb,nvirb))*1j - .5-.5j) - t2aa = (np.random.random((nkpts,nkpts,nkpts,nocca,nocca,nvira,nvira)) + - np.random.random((nkpts,nkpts,nkpts,nocca,nocca,nvira,nvira))*1j - .5-.5j) - kconserv = kpts_helper.get_kconserv(kmf.cell, kmf.kpts) - t2aa = t2aa - t2aa.transpose(1,0,2,4,3,5,6) - tmp = t2aa.copy() - for ki, kj, kk in kpts_helper.loop_kkk(nkpts): - kl = kconserv[ki, kk, kj] - t2aa[ki,kj,kk] = t2aa[ki,kj,kk] - tmp[ki,kj,kl].transpose(0,1,3,2) - t2ab = (np.random.random((nkpts,nkpts,nkpts,nocca,noccb,nvira,nvirb)) + - np.random.random((nkpts,nkpts,nkpts,nocca,noccb,nvira,nvirb))*1j - .5-.5j) - t2bb = (np.random.random((nkpts,nkpts,nkpts,noccb,noccb,nvirb,nvirb)) + - np.random.random((nkpts,nkpts,nkpts,noccb,noccb,nvirb,nvirb))*1j - .5-.5j) - t2bb = t2bb - t2bb.transpose(1,0,2,4,3,5,6) - tmp = t2bb.copy() - for ki, kj, kk in kpts_helper.loop_kkk(nkpts): - kl = kconserv[ki, kk, kj] - t2bb[ki,kj,kk] = t2bb[ki,kj,kk] - tmp[ki,kj,kl].transpose(0,1,3,2) - - t1 = (t1a, t1b) - t2 = (t2aa, t2ab, t2bb) - return t1, t2 - - mycc = KUCCSD(kmf) - eris = mycc.ao2mo() - t1, t2 = rand_t1_t2(mycc) - Ht1, Ht2 = mycc.update_amps(t1, t2, eris) - print(lib.fp(Ht1[0]) - (2.2677885702176339-2.5150764056992041j)) - print(lib.fp(Ht1[1]) - (-51.643438947846086+526.58026126100458j)) - print(lib.fp(Ht2[0]) - (-29.490813482748258-8.7509143690136018j)) - print(lib.fp(Ht2[1]) - (2256.0440056839416-193.16480896707569j)) - print(lib.fp(Ht2[2]) - (-250.59447681063182-397.57189085666982j)) - - kmf.mo_occ[:] = 0 - kmf.mo_occ[:,:,:2] = 1 - mycc = KUCCSD(kmf) - eris = mycc.ao2mo() - t1, t2 = rand_t1_t2(mycc) - Ht1, Ht2 = mycc.update_amps(t1, t2, eris) - print(lib.fp(Ht1[0]) - (5.4622516572705662+1.990046725028729j)) - print(lib.fp(Ht1[1]) - (4.8801120611799043-5.9940463787453488j)) - print(lib.fp(Ht2[0]) - (-192.38864512375193+305.14191018543983j)) - print(lib.fp(Ht2[1]) - (23085.044505825954-11527.802302550244j)) - print(lib.fp(Ht2[2]) - (115.57932548288559-40.888597453928604j)) - - from pyscf.pbc.cc import kccsd - kgcc = kccsd.GCCSD(scf.addons.convert_to_ghf(kmf)) - kccsd_eris = kccsd._make_eris_incore(kgcc, kgcc._scf.mo_coeff) - r1 = kgcc.spatial2spin(t1) - r2 = kgcc.spatial2spin(t2) - ge = kccsd.energy(kgcc, r1, r2, kccsd_eris) - r1, r2 = kgcc.update_amps(r1, r2, kccsd_eris) - ue = energy(mycc, t1, t2, eris) - print(abs(ge - ue)) - print(abs(r1 - kgcc.spatial2spin(Ht1)).max()) - print(abs(r2 - kgcc.spatial2spin(Ht2)).max()) - - kmf = kmf.density_fit(auxbasis=[[0, (1., 1.)]]) - mycc = KUCCSD(kmf) - eris = _make_df_eris(mycc, mycc.mo_coeff) - t1, t2 = rand_t1_t2(mycc) - Ht1, Ht2 = mycc.update_amps(t1, t2, eris) - - print(lib.fp(Ht1[0]) - (6.9341372555790013+0.87313546297025901j)) - print(lib.fp(Ht1[1]) - (6.7538005829391992-0.95702422534126796j)) - print(lib.fp(Ht2[0]) - (-509.24544842179876+448.00925776269855j)) - print(lib.fp(Ht2[1]) - (107.5960392010511+40.869216223808067j) ) - print(lib.fp(Ht2[2]) - (-196.75910296082139+218.53005038057515j)) - kgcc = kccsd.GCCSD(scf.addons.convert_to_ghf(kmf)) - kccsd_eris = kccsd._make_eris_incore(kgcc, kgcc._scf.mo_coeff) - r1 = kgcc.spatial2spin(t1) - r2 = kgcc.spatial2spin(t2) - ge = kccsd.energy(kgcc, r1, r2, kccsd_eris) - r1, r2 = kgcc.update_amps(r1, r2, kccsd_eris) - print(abs(r1 - kgcc.spatial2spin(Ht1)).max()) - print(abs(r2 - kgcc.spatial2spin(Ht2)).max()) - - print(all([abs(lib.fp(eris.oooo) - (-0.18290712163391809-0.13839081039521306j) )<1e-8, - abs(lib.fp(eris.ooOO) - (-0.084752145202964035-0.28496525042110676j) )<1e-8, - #abs(lib.fp(eris.OOoo) - (0.43054922768629345-0.27990237216969871j) )<1e-8, - abs(lib.fp(eris.OOOO) - (-0.2941475969103261-0.047247498899840978j) )<1e-8, - abs(lib.fp(eris.ooov) - (0.23381463349517045-0.11703340936984277j) )<1e-8, - abs(lib.fp(eris.ooOV) - (-0.052655392703214066+0.69533309442418556j) )<1e-8, - abs(lib.fp(eris.OOov) - (-0.2111361247200903+0.85087916975274647j) )<1e-8, - abs(lib.fp(eris.OOOV) - (-0.36995992208047412-0.18887278030885621j) )<1e-8, - abs(lib.fp(eris.oovv) - (0.21107397525051516+0.0048714991438174871j) )<1e-8, - abs(lib.fp(eris.ooVV) - (-0.076411225687065987+0.11080438166425896j) )<1e-8, - abs(lib.fp(eris.OOvv) - (-0.17880337626095003-0.24174716216954206j) )<1e-8, - abs(lib.fp(eris.OOVV) - (0.059186286356424908+0.68433866387500164j) )<1e-8, - abs(lib.fp(eris.ovov) - (0.15402983765151051+0.064359681685222214j) )<1e-8, - abs(lib.fp(eris.ovOV) - (-0.10697649196044598+0.30351249676253234j) )<1e-8, - #abs(lib.fp(eris.OVov) - (-0.17619329728836752-0.56585020976035816j) )<1e-8, - abs(lib.fp(eris.OVOV) - (-0.63963235318492118+0.69863219317718828j) )<1e-8, - abs(lib.fp(eris.voov) - (-0.24137641647339092+0.18676684336011531j) )<1e-8, - abs(lib.fp(eris.voOV) - (0.19257709151227204+0.38929027819406414j) )<1e-8, - #abs(lib.fp(eris.VOov) - (0.07632606729926053-0.70350947950650355j) )<1e-8, - abs(lib.fp(eris.VOOV) - (-0.47970203195500816+0.46735207193861927j) )<1e-8, - abs(lib.fp(eris.vovv) - (-0.1342049915673903-0.23391327821719513j) )<1e-8, - abs(lib.fp(eris.voVV) - (-0.28989635223866056+0.9644368822688475j) )<1e-8, - abs(lib.fp(eris.VOvv) - (-0.32428269235420271+0.0029847254383674748j))<1e-8, - abs(lib.fp(eris.VOVV) - (0.45031779746222456-0.36858577475752041j) )<1e-8])) - - eris = _make_eris_outcore(mycc, mycc.mo_coeff) - print(all([abs(lib.fp(eris.oooo) - (-0.18290712163391809-0.13839081039521306j) )<1e-8, - abs(lib.fp(eris.ooOO) - (-0.084752145202964035-0.28496525042110676j) )<1e-8, - #abs(lib.fp(eris.OOoo) - (0.43054922768629345-0.27990237216969871j) )<1e-8, - abs(lib.fp(eris.OOOO) - (-0.2941475969103261-0.047247498899840978j) )<1e-8, - abs(lib.fp(eris.ooov) - (0.23381463349517045-0.11703340936984277j) )<1e-8, - abs(lib.fp(eris.ooOV) - (-0.052655392703214066+0.69533309442418556j) )<1e-8, - abs(lib.fp(eris.OOov) - (-0.2111361247200903+0.85087916975274647j) )<1e-8, - abs(lib.fp(eris.OOOV) - (-0.36995992208047412-0.18887278030885621j) )<1e-8, - abs(lib.fp(eris.oovv) - (0.21107397525051516+0.0048714991438174871j) )<1e-8, - abs(lib.fp(eris.ooVV) - (-0.076411225687065987+0.11080438166425896j) )<1e-8, - abs(lib.fp(eris.OOvv) - (-0.17880337626095003-0.24174716216954206j) )<1e-8, - abs(lib.fp(eris.OOVV) - (0.059186286356424908+0.68433866387500164j) )<1e-8, - abs(lib.fp(eris.ovov) - (0.15402983765151051+0.064359681685222214j) )<1e-8, - abs(lib.fp(eris.ovOV) - (-0.10697649196044598+0.30351249676253234j) )<1e-8, - #abs(lib.fp(eris.OVov) - (-0.17619329728836752-0.56585020976035816j) )<1e-8, - abs(lib.fp(eris.OVOV) - (-0.63963235318492118+0.69863219317718828j) )<1e-8, - abs(lib.fp(eris.voov) - (-0.24137641647339092+0.18676684336011531j) )<1e-8, - abs(lib.fp(eris.voOV) - (0.19257709151227204+0.38929027819406414j) )<1e-8, - #abs(lib.fp(eris.VOov) - (0.07632606729926053-0.70350947950650355j) )<1e-8, - abs(lib.fp(eris.VOOV) - (-0.47970203195500816+0.46735207193861927j) )<1e-8, - abs(lib.fp(eris.vovv) - (-0.1342049915673903-0.23391327821719513j) )<1e-8, - abs(lib.fp(eris.voVV) - (-0.28989635223866056+0.9644368822688475j) )<1e-8, - abs(lib.fp(eris.VOvv) - (-0.32428269235420271+0.0029847254383674748j))<1e-8, - abs(lib.fp(eris.VOVV) - (0.45031779746222456-0.36858577475752041j) )<1e-8, - abs(lib.fp(eris.vvvv) - (-0.080512851258903173-0.2868384266725581j) )<1e-8, - abs(lib.fp(eris.vvVV) - (-0.5137063762484736+1.1036785801263898j) )<1e-8, - #abs(lib.fp(eris.VVvv) - (0.16468487082491939+0.25730725586992997j) )<1e-8, - abs(lib.fp(eris.VVVV) - (-0.56714875196802295+0.058636785679170501j) )<1e-8])) diff --git a/pyscf/pbc/cc/kintermediates_uhf.py b/pyscf/pbc/cc/kintermediates_uhf.py index 485e2189ef..ecb332e092 100644 --- a/pyscf/pbc/cc/kintermediates_uhf.py +++ b/pyscf/pbc/cc/kintermediates_uhf.py @@ -214,9 +214,9 @@ def cc_Woooo(cc, t1, t2, eris): WOOOO[km,ki,kn] += tmp_bbbbJ[ki,kj] WooOO[km,ki,kn] += tmp_aabbJ[ki,kj] WooOO[kn,ki,km] -= tmp_baabJ[ki,kj].transpose(0,3,2,1,4) - Woooo[km,ki,kn] += 0.25*einsum('yxijef,xmenf->yminj', tau_aa[ki,kj], eris.ovov[km,ki,kn]) - WOOOO[km,ki,kn] += 0.25*einsum('yxijef,xmenf->yminj', tau_bb[ki,kj], eris.OVOV[km,ki,kn]) - WooOO[km,ki,kn] += 0.5*einsum('yxijef,xmenf->yminj', tau_ab[ki,kj], eris.ovOV[km,ki,kn]) + Woooo[km,ki,kn] += 0.25*einsum('yxijef,xmenf->yminj', tau_aa[ki,kj], eris.ovov[km,:,kn]) + WOOOO[km,ki,kn] += 0.25*einsum('yxijef,xmenf->yminj', tau_bb[ki,kj], eris.OVOV[km,:,kn]) + WooOO[km,ki,kn] += 0.5*einsum('yxijef,xmenf->yminj', tau_ab[ki,kj], eris.ovOV[km,:,kn]) Woooo = Woooo - Woooo.transpose(2,1,0,5,4,3,6) WOOOO = WOOOO - WOOOO.transpose(2,1,0,5,4,3,6) diff --git a/pyscf/pbc/cc/kuccsd_rdm.py b/pyscf/pbc/cc/kuccsd_rdm.py index 792050a7eb..a81325491e 100644 --- a/pyscf/pbc/cc/kuccsd_rdm.py +++ b/pyscf/pbc/cc/kuccsd_rdm.py @@ -155,74 +155,3 @@ def _make_rdm1(mycc, d1, with_frozen=True, ao_repr=False): dm1a = lib.einsum('xpi,xij,xqj->xpq', mo_a, dm1a, mo_a.conj()) dm1b = lib.einsum('xpi,xij,xqj->xpq', mo_b, dm1b, mo_b.conj()) return dm1a, dm1b - - - -if __name__ == '__main__': - import numpy as np - from pyscf.pbc import gto - from pyscf.pbc import scf - from pyscf import lib - from pyscf.pbc.cc import KUCCSD - - a0 = 3.0 - vac = 200 - nmp = [1,1,1] - dis = 10 - vec = [[ 3.0/2.0*a0*dis,np.sqrt(3.0)/2.0*a0*dis, 0], - [-3.0/2.0*a0*dis,np.sqrt(3.0)/2.0*a0*dis, 0], - [ 0, 0,vac]] - bas = 'cc-pvdz' - pos = [['H',(-a0/2.0,0,0)], - ['H',( a0/2.0,0,0)]] - - cell = gto.M(unit='B',a=vec,atom=pos,basis=bas,verbose=4) - cell.precision = 1e-11 - - kpts = cell.make_kpts(nmp) - kmf = scf.KUHF(cell,kpts,exxdiv=None).density_fit() - kmf.max_cycle=250 - - dm = kmf.get_init_guess() - aoind = cell.aoslice_by_atom() - idx0, idx1 = aoind - dm[0,:,idx0[2]:idx0[3], idx0[2]:idx0[3]] = 0 - dm[0,:,idx1[2]:idx1[3], idx1[2]:idx1[3]] = dm[0,:,idx1[2]:idx1[3], idx1[2]:idx1[3]] * 2 - dm[1,:,idx0[2]:idx0[3], idx0[2]:idx0[3]] = dm[1,:,idx0[2]:idx0[3], idx0[2]:idx0[3]] * 2 - dm[1,:,idx1[2]:idx1[3], idx1[2]:idx1[3]] = 0 - - ehf = kmf.kernel(dm) - - kcc = KUCCSD(kmf) - ecc, t1, t2 = kcc.kernel() - - dm1a,dm1b = make_rdm1(kcc, t1, t2, l1=None, l2=None) - - ((dooa, doob), (dova, dovb), (dvoa, dvob), (dvva, dvvb)) = _gamma1_intermediates(kcc, t1, t2, l1=t1, l2=t2) - - print(np.linalg.norm(dooa)) - print(np.linalg.norm(doob)) - print(np.linalg.norm(dova)) - print(np.linalg.norm(dovb)) - print(np.linalg.norm(dvoa)) - print(np.linalg.norm(dvob)) - print(np.linalg.norm(dvva)) - print(np.linalg.norm(dvvb)) - - t1a, t1b = t1 - t2aa, t2ab, t2bb = t2 - - print(np.linalg.norm(t1a)) - print(np.linalg.norm(t1b)) - print(np.linalg.norm(t2aa)) - print(np.linalg.norm(t2ab)) - print(np.linalg.norm(t2bb)) - - - #print(kmf.spin_square()) - print(np.linalg.norm(dm1a)) - print(np.linalg.norm(dm1b)) - - print(np.linalg.norm(dm1a - dm1b)) - print(dm1a[0]) - print(dm1b[0]) diff --git a/pyscf/pbc/cc/test/test_kccsd_ksymm.py b/pyscf/pbc/cc/test/test_kccsd_ksymm.py index b203d01043..24a404a565 100644 --- a/pyscf/pbc/cc/test/test_kccsd_ksymm.py +++ b/pyscf/pbc/cc/test/test_kccsd_ksymm.py @@ -18,13 +18,15 @@ import unittest import numpy as np +from pyscf import lib from pyscf.pbc import gto, scf, cc def setUpModule(): - global kccref, kmf + global kmf, kcc L = 2. He = gto.Cell() He.verbose = 5 + He.output = '/dev/null' He.a = np.eye(3)*L He.atom =[['He' , ( L/2+0., L/2+0., L/2+0.)],] He.basis = {'He': [[0, (4.0, 1.0)], [0, (1.0, 1.0)]]} @@ -32,31 +34,32 @@ def setUpModule(): He.build() nk = [2,2,2] - - kpts0 = He.make_kpts(nk) - kmf0 = scf.KRHF(He, kpts0, exxdiv=None).density_fit() - kmf0.kernel() - kccref = cc.KRCCSD(kmf0) - kccref.kernel() - kpts = He.make_kpts(nk,space_group_symmetry=True,time_reversal_symmetry=True) kmf = scf.KRHF(He, kpts, exxdiv=None).density_fit() kmf.kernel() + kcc = cc.KsymAdaptedRCCSD(kmf) + kcc.kernel() def tearDownModule(): - global kccref, kmf - del kccref, kmf + global kmf, kcc + del kmf, kcc class KnownValues(unittest.TestCase): - def test_krccsd(self): - kcc = cc.KsymAdaptedRCCSD(kmf) - kcc.kernel() - self.assertAlmostEqual(kcc.e_corr, kccref.e_corr, 8) + def test_krccsd_ksym(self): + self.assertAlmostEqual(kcc.e_corr, -0.0073791230287231875, 8) + + t1 = kcc.t1.todense() + self.assertAlmostEqual(lib.fp(t1).max(), 0.0009185224804818333, 7) + + def test_vs_krccsd(self): + kmf0 = kmf.to_khf() + kccref = cc.krccsd.KRCCSD(kmf0) + kccref.kernel() t1 = kcc.t1.todense() self.assertAlmostEqual(abs(kccref.t1 - t1).max(), 0, 6) - #t2 = kcc.t2.todense() - #self.assertAlmostEqual(abs(kccref.t2 - t2).max(), 0, 6) + t2 = kcc.t2.todense() + self.assertAlmostEqual(abs(kccref.t2 - t2).max(), 0, 6) if __name__ == '__main__': print("Full Tests for CCSD with k-point symmetry") diff --git a/pyscf/pbc/cc/test/test_kuccsd.py b/pyscf/pbc/cc/test/test_kuccsd.py index 5035beef05..5ed1132f02 100644 --- a/pyscf/pbc/cc/test/test_kuccsd.py +++ b/pyscf/pbc/cc/test/test_kuccsd.py @@ -16,15 +16,19 @@ import unittest import numpy +import numpy as np from pyscf import lib +from pyscf import lo from pyscf.pbc import gto from pyscf.pbc import scf from pyscf.pbc.cc import kccsd_uhf from pyscf.pbc.cc import kccsd +from pyscf.pbc.cc import kuccsd_rdm from pyscf.pbc.cc import eom_kccsd_ghf from pyscf.pbc.cc import eom_kccsd_uhf from pyscf.pbc.lib import kpts_helper + def setUpModule(): global cell cell = gto.Cell() @@ -119,6 +123,65 @@ def test_spatial2spin(self): self.assertAlmostEqual(abs(r2[1]-t2[1]).max(), 0, 12) self.assertAlmostEqual(abs(r2[2]-t2[2]).max(), 0, 12) + def test_update_amps(self): + np.random.seed(2) + # Running HF and CCSD with 1x1x2 Monkhorst-Pack k-point mesh + kmf = scf.KUHF(cell, kpts=cell.make_kpts([1,1,3]), exxdiv=None) + nmo = cell.nao_nr() + kmf.mo_occ = np.zeros((2,3,nmo)) + kmf.mo_occ[0,:,:3] = 1 + kmf.mo_occ[1,:,:1] = 1 + kmf.mo_energy = np.arange(nmo) + np.random.random((2,3,nmo)) * .3 + kmf.mo_energy[kmf.mo_occ == 0] += 2 + + mo = (np.random.random((2,3,nmo,nmo)) + + np.random.random((2,3,nmo,nmo))*1j - .5-.5j) + s = kmf.get_ovlp() + kmf.mo_coeff = np.empty_like(mo) + nkpts = len(kmf.kpts) + for k in range(nkpts): + kmf.mo_coeff[0,k] = lo.orth.vec_lowdin(mo[0,k], s[k]) + kmf.mo_coeff[1,k] = lo.orth.vec_lowdin(mo[1,k], s[k]) + + mycc = kccsd_uhf.KUCCSD(kmf) + eris = mycc.ao2mo() + np.random.seed(1) + nocc = mycc.nocc + nvir = nmo - nocc[0], nmo - nocc[1] + t1, t2 = rand_t1_t2(cell, kmf.kpts, nocc, nvir) + Ht1, Ht2 = mycc.update_amps(t1, t2, eris) + self.assertAlmostEqual(lib.fp(Ht1[0]), (2.2677885702176339-2.5150764056992041j ), 4) + self.assertAlmostEqual(lib.fp(Ht1[1]), (-51.643438947846086+526.58026126100458j), 4) + self.assertAlmostEqual(lib.fp(Ht2[0]), (-29.490813482748258-8.7509143690136018j), 4) + self.assertAlmostEqual(lib.fp(Ht2[1]), (2256.0435017189384-193.16481690849486j ), 4) + self.assertAlmostEqual(lib.fp(Ht2[2]), (-250.59447681063182-397.57189085666982j), 4) + + kmf.mo_occ[:] = 0 + kmf.mo_occ[:,:,:2] = 1 + mycc = kccsd_uhf.KUCCSD(kmf) + eris = mycc.ao2mo() + np.random.seed(1) + nocc = mycc.nocc + nvir = nmo - nocc[0], nmo - nocc[1] + t1, t2 = rand_t1_t2(cell, kmf.kpts, nocc, nvir) + Ht1, Ht2 = mycc.update_amps(t1, t2, eris) + self.assertAlmostEqual(lib.fp(Ht1[0]), (5.4622516572705662+1.990046725028729j ), 4) + self.assertAlmostEqual(lib.fp(Ht1[1]), (4.8801120611799043-5.9940463787453488j), 4) + self.assertAlmostEqual(lib.fp(Ht2[0]), (-192.38864512375193+305.1419101854398j), 4) + self.assertAlmostEqual(lib.fp(Ht2[1]), (23085.056910282714-11527.808856109988j), 4) + self.assertAlmostEqual(lib.fp(Ht2[2]), (115.57932548288559-40.888597453928604j), 4) + + kgcc = kccsd.GCCSD(scf.addons.convert_to_ghf(kmf)) + kccsd_eris = kccsd._make_eris_incore(kgcc, kgcc._scf.mo_coeff) + r1 = kgcc.spatial2spin(t1) + r2 = kgcc.spatial2spin(t2) + ge = kccsd.energy(kgcc, r1, r2, kccsd_eris) + r1, r2 = kgcc.update_amps(r1, r2, kccsd_eris) + ue = kccsd_uhf.energy(mycc, t1, t2, eris) + self.assertAlmostEqual(abs(ge - ue), 0, 9) + self.assertAlmostEqual(abs(r1 - kgcc.spatial2spin(Ht1)).max(), 0, 9) + self.assertAlmostEqual(abs(r2 - kgcc.spatial2spin(Ht2)).max(), 0, 5) + def test_eris(self): numpy.random.seed(1) kpts = cell.make_kpts([1,1,3]) @@ -172,6 +235,86 @@ def test_eris(self): self.assertAlmostEqual(lib.fp(eris.VOvv), -0.036061642075282056 +0.019284185856121634j , 10) #self.assertAlmostEqual(lib.fp(eris.VVvv), 0.13458896578260207 -0.11322854172459119j , 10) + def test_df_eris(self): + np.random.seed(2) + # Running HF and CCSD with 1x1x2 Monkhorst-Pack k-point mesh + kmf = scf.KUHF(cell, kpts=cell.make_kpts([1,1,3]), exxdiv=None) + nmo = cell.nao_nr() + kmf.mo_occ = np.zeros((2,3,nmo)) + kmf.mo_occ[0,:,:3] = 1 + kmf.mo_occ[1,:,:1] = 1 + kmf.mo_energy = np.arange(nmo) + np.random.random((2,3,nmo)) * .3 + kmf.mo_energy[kmf.mo_occ == 0] += 2 + + mo = (np.random.random((2,3,nmo,nmo)) + + np.random.random((2,3,nmo,nmo))*1j - .5-.5j) + s = kmf.get_ovlp() + kmf.mo_coeff = np.empty_like(mo) + nkpts = len(kmf.kpts) + for k in range(nkpts): + kmf.mo_coeff[0,k] = lo.orth.vec_lowdin(mo[0,k], s[k]) + kmf.mo_coeff[1,k] = lo.orth.vec_lowdin(mo[1,k], s[k]) + + kmf.mo_occ[:] = 0 + kmf.mo_occ[:,:,:2] = 1 + kmf = kmf.density_fit(auxbasis=[[0, (1., 1.)]]) + mycc = kccsd_uhf.KUCCSD(kmf) + eris = kccsd_uhf._make_df_eris(mycc, mycc.mo_coeff) + self.assertAlmostEqual(lib.fp(eris.oooo), (-0.18290712163391809-0.13839081039521306j) , 7) + self.assertAlmostEqual(lib.fp(eris.ooOO), (-0.084752145202964035-0.28496525042110676j), 7) + #self.assertAlmostEqual(lib.fp(eris.OOoo), (0.43054922768629345-0.27990237216969871j) , 7) + self.assertAlmostEqual(lib.fp(eris.OOOO), (-0.2941475969103261-0.047247498899840978j) , 7) + self.assertAlmostEqual(lib.fp(eris.ooov), (0.23381463349517045-0.11703340936984277j) , 7) + self.assertAlmostEqual(lib.fp(eris.ooOV), (-0.052655392703214066+0.69533309442418556j), 7) + self.assertAlmostEqual(lib.fp(eris.OOov), (-0.2111361247200903+0.85087916975274647j) , 7) + self.assertAlmostEqual(lib.fp(eris.OOOV), (-0.36995992208047412-0.18887278030885621j) , 7) + self.assertAlmostEqual(lib.fp(eris.oovv), (0.21107397525051516+0.0048714991438174871j), 7) + self.assertAlmostEqual(lib.fp(eris.ooVV), (-0.076411225687065987+0.11080438166425896j), 7) + self.assertAlmostEqual(lib.fp(eris.OOvv), (-0.17880337626095003-0.24174716216954206j) , 7) + self.assertAlmostEqual(lib.fp(eris.OOVV), (0.059186286356424908+0.68433866387500164j) , 7) + self.assertAlmostEqual(lib.fp(eris.ovov), (0.15402983765151051+0.064359681685222214j) , 7) + self.assertAlmostEqual(lib.fp(eris.ovOV), (-0.10697649196044598+0.30351249676253234j) , 7) + #self.assertAlmostEqual(lib.fp(eris.OVov), (-0.17619329728836752-0.56585020976035816j) , 7) + self.assertAlmostEqual(lib.fp(eris.OVOV), (-0.63963235318492118+0.69863219317718828j) , 7) + self.assertAlmostEqual(lib.fp(eris.voov), (-0.24137641647339092+0.18676684336011531j) , 7) + self.assertAlmostEqual(lib.fp(eris.voOV), (0.19257709151227204+0.38929027819406414j) , 7) + #self.assertAlmostEqual(lib.fp(eris.VOov), (0.07632606729926053-0.70350947950650355j) , 7) + self.assertAlmostEqual(lib.fp(eris.VOOV), (-0.47970203195500816+0.46735207193861927j) , 7) + self.assertAlmostEqual(lib.fp(eris.vovv), (-0.1342049915673903-0.23391327821719513j) , 7) + self.assertAlmostEqual(lib.fp(eris.voVV), (-0.28989635223866056+0.9644368822688475j) , 7) + self.assertAlmostEqual(lib.fp(eris.VOvv), (-0.32428269235420271+0.00298472543836747j) , 7) + self.assertAlmostEqual(lib.fp(eris.VOVV), (0.45031779746222456-0.36858577475752041j) , 7) + + eris = kccsd_uhf._make_eris_outcore(mycc, mycc.mo_coeff) + self.assertAlmostEqual(lib.fp(eris.oooo), (-0.18290712163391809-0.13839081039521306j) , 7) + self.assertAlmostEqual(lib.fp(eris.ooOO), (-0.084752145202964035-0.28496525042110676j), 7) + #self.assertAlmostEqual(lib.fp(eris.OOoo), (0.43054922768629345-0.27990237216969871j) , 7) + self.assertAlmostEqual(lib.fp(eris.OOOO), (-0.2941475969103261-0.047247498899840978j) , 7) + self.assertAlmostEqual(lib.fp(eris.ooov), (0.23381463349517045-0.11703340936984277j) , 7) + self.assertAlmostEqual(lib.fp(eris.ooOV), (-0.052655392703214066+0.69533309442418556j), 7) + self.assertAlmostEqual(lib.fp(eris.OOov), (-0.2111361247200903+0.85087916975274647j) , 7) + self.assertAlmostEqual(lib.fp(eris.OOOV), (-0.36995992208047412-0.18887278030885621j) , 7) + self.assertAlmostEqual(lib.fp(eris.oovv), (0.21107397525051516+0.0048714991438174871j), 7) + self.assertAlmostEqual(lib.fp(eris.ooVV), (-0.076411225687065987+0.11080438166425896j), 7) + self.assertAlmostEqual(lib.fp(eris.OOvv), (-0.17880337626095003-0.24174716216954206j) , 7) + self.assertAlmostEqual(lib.fp(eris.OOVV), (0.059186286356424908+0.68433866387500164j) , 7) + self.assertAlmostEqual(lib.fp(eris.ovov), (0.15402983765151051+0.064359681685222214j) , 7) + self.assertAlmostEqual(lib.fp(eris.ovOV), (-0.10697649196044598+0.30351249676253234j) , 7) + #self.assertAlmostEqual(lib.fp(eris.OVov), (-0.17619329728836752-0.56585020976035816j) , 7) + self.assertAlmostEqual(lib.fp(eris.OVOV), (-0.63963235318492118+0.69863219317718828j) , 7) + self.assertAlmostEqual(lib.fp(eris.voov), (-0.24137641647339092+0.18676684336011531j) , 7) + self.assertAlmostEqual(lib.fp(eris.voOV), (0.19257709151227204+0.38929027819406414j) , 7) + #self.assertAlmostEqual(lib.fp(eris.VOov), (0.07632606729926053-0.70350947950650355j) , 7) + self.assertAlmostEqual(lib.fp(eris.VOOV), (-0.47970203195500816+0.46735207193861927j) , 7) + self.assertAlmostEqual(lib.fp(eris.vovv), (-0.1342049915673903-0.23391327821719513j) , 7) + self.assertAlmostEqual(lib.fp(eris.voVV), (-0.28989635223866056+0.9644368822688475j) , 7) + self.assertAlmostEqual(lib.fp(eris.VOvv), (-0.32428269235420271+0.002984725438367474j), 7) + self.assertAlmostEqual(lib.fp(eris.VOVV), (0.45031779746222456-0.36858577475752041j) , 7) + self.assertAlmostEqual(lib.fp(eris.vvvv), (-0.080512851258903173-0.2868384266725581j) , 7) + self.assertAlmostEqual(lib.fp(eris.vvVV), (-0.5137063762484736+1.1036785801263898j) , 7) + #self.assertAlmostEqual(lib.fp(eris.VVvv), (0.16468487082491939+0.25730725586992997j) , 7) + self.assertAlmostEqual(lib.fp(eris.VVVV), (-0.56714875196802295+0.058636785679170501j), 7) + def test_spatial2spin_ip(self): numpy.random.seed(1) kpts = cell.make_kpts([1,2,3]) @@ -234,6 +377,58 @@ def test_spatial2spin_ea(self): self.assertAlmostEqual(abs(r1-spin_r1_ea).max(), 0, 12) self.assertAlmostEqual(abs(r2-spin_r2_ea).max(), 0, 12) + # TODO: more tests are needed for kuccsd_rdm + def test_kuccsd_rdm_slow(self): + a0 = 3.0 + nmp = [3,1,1] + vec = [[ 3.0/2.0*a0,np.sqrt(3.0)/2.0*a0, 0], + [-3.0/2.0*a0,np.sqrt(3.0)/2.0*a0, 0], + [ 0, 0,a0]] + bas = [[0, [2., 1.]], [0, [.4, 1.]]] + pos = [['H',(-a0/2.0,0,0)], + ['H',( a0/2.0,0,0)]] + + cell = gto.M(unit='B',a=vec,atom=pos,basis=bas,verbose=4,output='/dev/null') + + kpts = cell.make_kpts(nmp) + kmf = scf.KUHF(cell,kpts,exxdiv=None).density_fit() + + dm = kmf.get_init_guess() + aoind = cell.aoslice_by_atom() + idx0, idx1 = aoind + dm[0,:,idx0[2]:idx0[3], idx0[2]:idx0[3]] = 0 + dm[0,:,idx1[2]:idx1[3], idx1[2]:idx1[3]] = dm[0,:,idx1[2]:idx1[3], idx1[2]:idx1[3]] * 2 + dm[1,:,idx0[2]:idx0[3], idx0[2]:idx0[3]] = dm[1,:,idx0[2]:idx0[3], idx0[2]:idx0[3]] * 2 + dm[1,:,idx1[2]:idx1[3], idx1[2]:idx1[3]] = 0 + + ehf = kmf.kernel(dm) + + kcc = kmf.CCSD() + ecc, t1, t2 = kcc.kernel() + + dm1a,dm1b = kuccsd_rdm.make_rdm1(kcc, t1, t2, l1=None, l2=None) + + ((dooa, doob), (dova, dovb), (dvoa, dvob), (dvva, dvvb)) = kuccsd_rdm._gamma1_intermediates(kcc, t1, t2, l1=t1, l2=t2) + t1a, t1b = t1 + t2aa, t2ab, t2bb = t2 + + self.assertAlmostEqual(np.linalg.norm(dooa), 0.196350361437495, 7) + self.assertAlmostEqual(np.linalg.norm(doob), 0.192418260688681, 7) + self.assertAlmostEqual(np.linalg.norm(dova), 0.343229849478994, 7) + self.assertAlmostEqual(np.linalg.norm(dovb), 0.343226810305894, 7) + self.assertAlmostEqual(np.linalg.norm(dvoa), 0.38420973758367 , 7) + self.assertAlmostEqual(np.linalg.norm(dvob), 0.386371342174770, 7) + self.assertAlmostEqual(np.linalg.norm(dvva), 0.196174979640156, 7) + self.assertAlmostEqual(np.linalg.norm(dvvb), 0.202569918050194, 7) + self.assertAlmostEqual(np.linalg.norm(t1a), 0.343229849478994, 7) + self.assertAlmostEqual(np.linalg.norm(t1b), 0.343226810305894, 7) + self.assertAlmostEqual(np.linalg.norm(t2aa), 0.114653158783553, 7) + self.assertAlmostEqual(np.linalg.norm(t2ab), 0.484479199870810, 7) + self.assertAlmostEqual(np.linalg.norm(t2bb), 0.114653041996545, 7) + self.assertAlmostEqual(np.linalg.norm(dm1a), 1.622541608251682, 7) + self.assertAlmostEqual(np.linalg.norm(dm1b), 1.622541545604457, 7) + self.assertAlmostEqual(np.linalg.norm(dm1a - dm1b), 0.088686084829135, 7) + if __name__ == '__main__': print("KUCCSD tests") unittest.main() diff --git a/pyscf/pbc/ci/test/test_cisd.py b/pyscf/pbc/ci/test/test_cisd.py new file mode 100644 index 0000000000..3a3bfe16b9 --- /dev/null +++ b/pyscf/pbc/ci/test/test_cisd.py @@ -0,0 +1,41 @@ +import unittest +import numpy as np +from pyscf import lib +from pyscf.pbc import gto, scf +from pyscf.pbc.ci import cisd + +def setUpModule(): + global cell, kmf, kci, eris + cell = gto.Cell() + cell.a = np.eye(3) * 2.5 + cell.mesh = [11] * 3 + cell.atom = '''He 0. 2. 1.5 + He 1. 1. 1.''' + cell.basis = {'He': [(0, (1.5, 1)), (0, (1., 1))]} + cell.build() + +def tearDownModule(): + global cell + del cell + +class KnownValues(unittest.TestCase): + def test_cisd(self): + mf = scf.RHF(cell).run() + myci = mf.CISD().run() + self.assertTrue(isinstance(myci, cisd.RCISD)) + self.assertAlmostEqual(myci.e_corr, -0.0107155147353, 9) + + umf = mf.to_uhf() + myci = umf.CISD().run() + self.assertTrue(isinstance(myci, cisd.UCISD)) + self.assertAlmostEqual(myci.e_corr, -0.0107155147353, 9) + + gmf = mf.to_ghf() + myci = gmf.CISD().run() + self.assertTrue(isinstance(myci, cisd.GCISD)) + self.assertAlmostEqual(myci.e_corr, -0.0107155147353, 9) + + +if __name__ == "__main__": + print("Full Tests for gamma-point CISD") + unittest.main() diff --git a/pyscf/pbc/gw/krgw_ac.py b/pyscf/pbc/gw/krgw_ac.py index 5ecb167c75..42487058d9 100644 --- a/pyscf/pbc/gw/krgw_ac.py +++ b/pyscf/pbc/gw/krgw_ac.py @@ -86,6 +86,9 @@ def kernel(gw, mo_energy, mo_coeff, orbs=None, raise RuntimeError('Found incompatible integral scheme %s.' 'KGWAC can be only used with GDF integrals' % gw.with_df.__class__) + if rhf.with_df._j_only: + logger.debug(gw, 'Rebuild CDERI for exchange integrals') + rhf.with_df.build(j_only=False) vk = rhf.get_veff(gw.mol,dm_kpts=dm) - rhf.get_j(gw.mol,dm_kpts=dm) for k in range(nkpts): @@ -639,69 +642,3 @@ def kernel(self, mo_energy=None, mo_coeff=None, orbs=None, kptlist=None, nw=100) logger.warn(self, 'GW QP energies may not be sorted from min to max') logger.timer(self, 'GW', *cput0) return self.mo_energy - -if __name__ == '__main__': - from pyscf.pbc import gto - from pyscf.pbc.lib import chkfile - import os - # This test takes a few minutes - cell = gto.Cell() - cell.build(unit = 'angstrom', - a = ''' - 0.000000 1.783500 1.783500 - 1.783500 0.000000 1.783500 - 1.783500 1.783500 0.000000 - ''', - atom = 'C 1.337625 1.337625 1.337625; C 2.229375 2.229375 2.229375', - dimension = 3, - max_memory = 8000, - verbose = 4, - pseudo = 'gth-pade', - basis='gth-szv', - precision=1e-10) - - kpts = cell.make_kpts([3,1,1],scaled_center=[0,0,0]) - gdf = df.GDF(cell, kpts) - gdf_fname = 'gdf_ints_311.h5' - gdf._cderi_to_save = gdf_fname - if not os.path.isfile(gdf_fname): - gdf.build() - - chkfname = 'diamond_311.chk' - if os.path.isfile(chkfname): - kmf = dft.KRKS(cell, kpts) - kmf.xc = 'pbe' - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - data = chkfile.load(chkfname, 'scf') - kmf.__dict__.update(data) - else: - kmf = dft.KRKS(cell, kpts) - kmf.xc = 'pbe' - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - kmf.conv_tol = 1e-12 - kmf.chkfile = chkfname - kmf.kernel() - - gw = KRGWAC(kmf) - gw.linearized = False - gw.ac = 'pade' - # without finite size corrections - gw.fc = False - nocc = gw.nocc - gw.kernel(kptlist=[0,1,2],orbs=range(0,nocc+3)) - print(gw.mo_energy) - assert ((abs(gw.mo_energy[0][nocc-1]-0.62045797))<1e-5) - assert ((abs(gw.mo_energy[0][nocc]-0.96574324))<1e-5) - assert ((abs(gw.mo_energy[1][nocc-1]-0.52639137))<1e-5) - assert ((abs(gw.mo_energy[1][nocc]-1.07513258))<1e-5) - - # with finite size corrections - gw.fc = True - gw.kernel(kptlist=[0,1,2],orbs=range(0,nocc+3)) - print(gw.mo_energy) - assert ((abs(gw.mo_energy[0][nocc-1]-0.54277092))<1e-5) - assert ((abs(gw.mo_energy[0][nocc]-0.80148537))<1e-5) - assert ((abs(gw.mo_energy[1][nocc-1]-0.45073793))<1e-5) - assert ((abs(gw.mo_energy[1][nocc]-0.92910108))<1e-5) diff --git a/pyscf/pbc/gw/krgw_cd.py b/pyscf/pbc/gw/krgw_cd.py index a4501ed84b..47454471c6 100644 --- a/pyscf/pbc/gw/krgw_cd.py +++ b/pyscf/pbc/gw/krgw_cd.py @@ -84,6 +84,9 @@ def kernel(gw, mo_energy, mo_coeff, orbs=None, raise RuntimeError('Found incompatible integral scheme %s.' 'KGWCD can be only used with GDF integrals' % gw.with_df.__class__) + if rhf.with_df._j_only: + logger.debug(gw, 'Rebuild CDERI for exchange integrals') + rhf.with_df.build(j_only=False) vk = rhf.get_veff(gw.mol,dm_kpts=dm) - rhf.get_j(gw.mol,dm_kpts=dm) for k in range(nkpts): @@ -699,68 +702,3 @@ def kernel(self, mo_energy=None, mo_coeff=None, orbs=None, kptlist=None, nw=100) logger.warn(self, 'GW QP energies may not be sorted from min to max') logger.timer(self, 'GW', *cput0) return self.mo_energy - -if __name__ == '__main__': - from pyscf.pbc import gto - from pyscf.pbc.lib import chkfile - import os - # This test takes a few minutes - cell = gto.Cell() - cell.build(unit = 'angstrom', - a = ''' - 0.000000 1.783500 1.783500 - 1.783500 0.000000 1.783500 - 1.783500 1.783500 0.000000 - ''', - atom = 'C 1.337625 1.337625 1.337625; C 2.229375 2.229375 2.229375', - dimension = 3, - max_memory = 8000, - verbose = 4, - pseudo = 'gth-pade', - basis='gth-szv', - precision=1e-10) - - kpts = cell.make_kpts([3,1,1],scaled_center=[0,0,0]) - gdf = df.GDF(cell, kpts) - gdf_fname = 'gdf_ints_311.h5' - gdf._cderi_to_save = gdf_fname - if not os.path.isfile(gdf_fname): - gdf.build() - - chkfname = 'diamond_311.chk' - if os.path.isfile(chkfname): - kmf = dft.KRKS(cell, kpts) - kmf.xc = 'pbe' - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - data = chkfile.load(chkfname, 'scf') - kmf.__dict__.update(data) - else: - kmf = dft.KRKS(cell, kpts) - kmf.xc = 'pbe' - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - kmf.conv_tol = 1e-12 - kmf.chkfile = chkfname - kmf.kernel() - - gw = KRGWCD(kmf) - gw.linearized = False - # without finite size corrections - gw.fc = False - nocc = gw.nocc - gw.kernel(kptlist=[0,1,2],orbs=range(0,nocc+3)) - print(gw.mo_energy) - assert ((abs(gw.mo_energy[0][nocc-1]-0.62045796))<1e-5) - assert ((abs(gw.mo_energy[0][nocc]-0.96574426))<1e-5) - assert ((abs(gw.mo_energy[1][nocc-1]-0.52639129))<1e-5) - assert ((abs(gw.mo_energy[1][nocc]-1.07442235))<1e-5) - - # with finite size corrections - gw.fc = True - gw.kernel(kptlist=[0,1,2],orbs=range(0,nocc+3)) - print(gw.mo_energy) - assert ((abs(gw.mo_energy[0][nocc-1]-0.5427707))<1e-5) - assert ((abs(gw.mo_energy[0][nocc]-0.80148557))<1e-5) - assert ((abs(gw.mo_energy[1][nocc-1]-0.45073751))<1e-5) - assert ((abs(gw.mo_energy[1][nocc]-0.92910117))<1e-5) diff --git a/pyscf/pbc/gw/test/test_krgw.py b/pyscf/pbc/gw/test/test_krgw.py index 49e33e091b..f7d037c956 100644 --- a/pyscf/pbc/gw/test/test_krgw.py +++ b/pyscf/pbc/gw/test/test_krgw.py @@ -9,7 +9,7 @@ from pyscf.pbc.gw import krgw_cd def setUpModule(): - global cell + global cell, kpts, gdf cell = gto.Cell() cell.build( a = ''' @@ -24,20 +24,16 @@ def setUpModule(): basis='gth-szv', precision=1e-8) + kpts = cell.make_kpts([3,1,1],scaled_center=[0,0,0]) + gdf = df.GDF(cell, kpts) + def tearDownModule(): - global cell + global cell, kpts, gdf cell.stdout.close() - del cell + del cell, kpts, gdf class KnownValues(unittest.TestCase): def test_gwac_pade_high_cost(self): - kpts = cell.make_kpts([3,1,1],scaled_center=[0,0,0]) - gdf = df.GDF(cell, kpts) - gdf._cderi = gdf_fname = os.path.realpath(os.path.join(__file__, '..', 'gdf_ints_311.h5')) - if 1: - gdf._cderi_to_save = gdf_fname - gdf.build() - kmf = dft.KRKS(cell, kpts).density_fit(with_df=gdf) kmf.xc = 'pbe' kmf.kernel() @@ -64,13 +60,6 @@ def test_gwac_pade_high_cost(self): self.assertAlmostEqual(gw.mo_energy[1][nocc] , 0.92910108, 4) def test_gwcd_high_cost(self): - kpts = cell.make_kpts([3,1,1],scaled_center=[0,0,0]) - gdf = df.GDF(cell, kpts) - gdf._cderi = gdf_fname = os.path.realpath(os.path.join(__file__, '..', 'gdf_ints_311.h5')) - if 1: - gdf._cderi_to_save = gdf_fname - gdf.build() - kmf = dft.KRKS(cell, kpts).density_fit(with_df=gdf) kmf.xc = 'pbe' kmf.kernel() @@ -95,6 +84,42 @@ def test_gwcd_high_cost(self): self.assertAlmostEqual(gw.mo_energy[1][nocc-1], 0.45073751, 4) self.assertAlmostEqual(gw.mo_energy[1][nocc], 0.92910117, 4) + def test_gw(self): + cell = gto.Cell() + cell.build(a = ''' + 0.000000 1.783500 1.783500 + 1.783500 0.000000 1.783500 + 1.783500 1.783500 0.000000 + ''', + atom = 'H 1.337625 1.337625 1.337625; H 2.229375 2.229375 2.229375', + verbose = 4, + output = '/dev/null', + basis=[[0, [2., 1.]], [0, [.5, 1.]]]) + + kpts = cell.make_kpts([3,1,1],scaled_center=[0,0,0]) + kmf = dft.KRKS(cell, kpts).density_fit().run() + + gw = krgw_ac.KRGWAC(kmf) + gw.linearized = True + gw.ac = 'pade' + # without finite size corrections + gw.fc = False + nocc = gw.nocc + gw.kernel(kptlist=[0,1,2],orbs=range(0,nocc+3)) + self.assertAlmostEqual(gw.mo_energy[0][nocc-1], -0.257088388010083, 7) + self.assertAlmostEqual(gw.mo_energy[0][nocc] , 0.7377021147675703, 7) + self.assertAlmostEqual(gw.mo_energy[1][nocc-1], -0.121872186953884, 7) + self.assertAlmostEqual(gw.mo_energy[1][nocc] , 0.570710170186033 , 7) + + # with finite size corrections + gw.linearized = False + gw.fc = True + gw.kernel(kptlist=[0,1,2],orbs=range(0,nocc+3)) + self.assertAlmostEqual(gw.mo_energy[0][nocc-1], -0.464099926108335, 7) + self.assertAlmostEqual(gw.mo_energy[0][nocc] , 0.7105306664244474, 7) + self.assertAlmostEqual(gw.mo_energy[1][nocc-1], -0.347704595829313, 7) + self.assertAlmostEqual(gw.mo_energy[1][nocc] , 0.552136080110482 , 7) + if __name__ == '__main__': print('Full Tests for KRGW') unittest.main() diff --git a/pyscf/pbc/lib/kpts.py b/pyscf/pbc/lib/kpts.py index d36cddfd5f..959d174483 100644 --- a/pyscf/pbc/lib/kpts.py +++ b/pyscf/pbc/lib/kpts.py @@ -1129,6 +1129,7 @@ class KQuartets: The class holding the symmetry relations between k-quartets. ''' def __init__(self, kpts): + assert isinstance(kpts, KPoints) self.kpts = kpts self.kqrts_ibz = None self.weights_ibz = None diff --git a/pyscf/pbc/scf/addons.py b/pyscf/pbc/scf/addons.py index 02ca4d166c..5ad6f30a08 100644 --- a/pyscf/pbc/scf/addons.py +++ b/pyscf/pbc/scf/addons.py @@ -132,9 +132,10 @@ def get_occ(self, mo_energy_kpts=None, mo_coeff_kpts=None): mo_occs = f_occ(mu[1], mo_es[1], sigma) self.entropy = self._get_entropy(mo_es[0], mo_occs[0], mu[0]) self.entropy += self._get_entropy(mo_es[1], mo_occs[1], mu[1]) + self.entropy /= nkpts - fermi = (mol_addons._get_fermi_level(mo_es[0], nocc[0]), - mol_addons._get_fermi_level(mo_es[1], nocc[1])) + fermi = (mol_addons._get_fermi(mo_es[0], nocc[0]), + mol_addons._get_fermi(mo_es[1], nocc[1])) logger.debug(self, ' Alpha-spin Fermi level %g Sum mo_occ_kpts = %s should equal nelec = %s', fermi[0], mo_occs[0].sum(), nocc[0]) logger.debug(self, ' Beta-spin Fermi level %g Sum mo_occ_kpts = %s should equal nelec = %s', @@ -166,12 +167,12 @@ def get_occ(self, mo_energy_kpts=None, mo_coeff_kpts=None): # If mu0 is given, fix mu instead of electron number. XXX -Chong Sun mu = self.mu0 mo_occs = f_occ(mu, mo_es, sigma) - self.entropy = self._get_entropy(mo_es, mo_occs, mu) + self.entropy = self._get_entropy(mo_es, mo_occs, mu) / nkpts if is_rhf: mo_occs *= 2 self.entropy *= 2 - fermi = mol_addons._get_fermi_level(mo_es, nocc) + fermi = mol_addons._get_fermi(mo_es, nocc) logger.debug(self, ' Fermi level %g Sum mo_occ_kpts = %s should equal nelec = %s', fermi, mo_occs.sum(), nelectron) logger.info(self, ' sigma = %g Optimized mu = %.12g entropy = %.12g', @@ -483,14 +484,15 @@ def convert_to_kscf(mf, out=None): scf.ghf.GHF : scf.kghf.KGHF , } mf = mol_addons._object_without_soscf(mf, known_cls, False) - if isinstance(mf, scf.uhf.UHF): - mf.mo_occ = mf.mo_occ[:,numpy.newaxis] - mf.mo_coeff = mf.mo_coeff[:,numpy.newaxis] - mf.mo_energy = mf.mo_energy[:,numpy.newaxis] - else: - mf.mo_occ = mf.mo_occ[numpy.newaxis] - mf.mo_coeff = mf.mo_coeff[numpy.newaxis] - mf.mo_energy = mf.mo_energy[numpy.newaxis] + if mf.mo_energy is not None: + if isinstance(mf, scf.uhf.UHF): + mf.mo_occ = mf.mo_occ[:,numpy.newaxis] + mf.mo_coeff = mf.mo_coeff[:,numpy.newaxis] + mf.mo_energy = mf.mo_energy[:,numpy.newaxis] + else: + mf.mo_occ = mf.mo_occ[numpy.newaxis] + mf.mo_coeff = mf.mo_coeff[numpy.newaxis] + mf.mo_energy = mf.mo_energy[numpy.newaxis] if out is None: return mf diff --git a/pyscf/pbc/scf/khf_ksymm.py b/pyscf/pbc/scf/khf_ksymm.py index 8d25f8cac7..5b57cf3534 100644 --- a/pyscf/pbc/scf/khf_ksymm.py +++ b/pyscf/pbc/scf/khf_ksymm.py @@ -134,6 +134,8 @@ class KsymAdaptedKSCF(khf.KSCF): KRHF with k-point symmetry """ + _keys = set(['use_ao_symmetry']) + get_occ = get_occ get_rho = get_rho energy_elec = energy_elec diff --git a/pyscf/pbc/scf/test/test_addons.py b/pyscf/pbc/scf/test/test_addons.py index dd9e9a0964..a67f7aa51d 100644 --- a/pyscf/pbc/scf/test/test_addons.py +++ b/pyscf/pbc/scf/test/test_addons.py @@ -94,13 +94,13 @@ def test_kuhf_smearing1(self): cell.build() nks = [2,1,1] mf = pscf.KUHF(cell, cell.make_kpts(nks)).density_fit() - mf = smearing_(mf, .1) + mf = pscf.addons.smearing_(mf, .1) mf.kernel() self.assertAlmostEqual(mf.e_tot, -5.56769351866668, 6) - mf = smearing_(mf, .1, mu0=0.351195741757) + mf = pscf.addons.smearing_(mf, .1, mu0=0.351195741757) mf.kernel() self.assertAlmostEqual(mf.e_tot, -5.56769351866668, 6) - mf = smearing_(mf, .1, method='gauss') + mf = pscf.addons.smearing_(mf, .1, method='gauss') mf.kernel() self.assertAlmostEqual(mf.e_tot, -5.56785857886738, 6) diff --git a/pyscf/qmmm/itrf.py b/pyscf/qmmm/itrf.py index adc0e1e2aa..bd75036d3c 100644 --- a/pyscf/qmmm/itrf.py +++ b/pyscf/qmmm/itrf.py @@ -91,7 +91,7 @@ def qmmm_for_scf(method, mm_mol): Args: mm_mol : MM Mole object ''' - assert (isinstance(method, (scf.hf.SCF, mcscf.casci.CASCI))) + assert (isinstance(method, (scf.hf.SCF, mcscf.casci.CASBase))) if isinstance(method, scf.hf.SCF): # Avoid to initialize QMMM twice @@ -400,7 +400,7 @@ def grad_nuc_mm(self, mol=None): # Inject QMMM interface wrapper to other modules scf.hf.SCF.QMMM = mm_charge -mcscf.casci.CASCI.QMMM = mm_charge +mcscf.casci.CASBase.QMMM = mm_charge grad.rhf.Gradients.QMMM = mm_charge_grad if __name__ == '__main__': diff --git a/pyscf/scf/addons.py b/pyscf/scf/addons.py index 8100314846..6cb52a2042 100644 --- a/pyscf/scf/addons.py +++ b/pyscf/scf/addons.py @@ -74,7 +74,7 @@ def nelec_cost_fn(m, mo_es, nocc): mo_occ = f_occ(m, mo_es, sigma) return (mo_occ.sum() - nocc)**2 - fermi = _get_fermi_level(mo_es, nocc) + fermi = _get_fermi(mo_es, nocc) res = scipy.optimize.minimize( nelec_cost_fn, fermi, args=(mo_es, nocc), method='Powell', options={'xtol': 1e-5, 'ftol': 1e-5, 'maxiter': 10000}) @@ -82,7 +82,7 @@ def nelec_cost_fn(m, mo_es, nocc): mo_occs = f_occ(mu, mo_es, sigma) return mu, mo_occs -def _get_fermi_level(mo_energy, nocc): +def _get_fermi(mo_energy, nocc): mo_e_sorted = numpy.sort(mo_energy) return mo_e_sorted[nocc-1] @@ -151,8 +151,7 @@ def get_occ(self, mo_energy=None, mo_coeff=None): mo_occs = f_occ(mu[1], mo_es[1], sigma) self.entropy = self._get_entropy(mo_es[0], mo_occs[0], mu[0]) self.entropy += self._get_entropy(mo_es[1], mo_occs[1], mu[1]) - fermi = (_get_fermi_level(mo_es[0], nocc[0]), - _get_fermi_level(mo_es[1], nocc[1])) + fermi = (_get_fermi(mo_es[0], nocc[0]), _get_fermi(mo_es[1], nocc[1])) logger.debug(self, ' Alpha-spin Fermi level %g Sum mo_occ = %s should equal nelec = %s', fermi[0], mo_occs[0].sum(), nocc[0]) @@ -974,93 +973,99 @@ def convert_to_ghf(mf, out=None, remove_df=False): def _update_mo_to_uhf_(mf, mf1): from pyscf import scf - if mf.mo_energy is not None: - if isinstance(mf, scf.uhf.UHF): - mf1.mo_occ = mf.mo_occ - mf1.mo_coeff = mf.mo_coeff - mf1.mo_energy = mf.mo_energy - elif getattr(mf, 'kpts', None) is None: # RHF/ROHF - mf1.mo_occ = numpy.array((mf.mo_occ>0, mf.mo_occ==2), dtype=numpy.double) - # ROHF orbital energies, not canonical UHF orbital energies - mo_ea = getattr(mf.mo_energy, 'mo_ea', mf.mo_energy) - mo_eb = getattr(mf.mo_energy, 'mo_eb', mf.mo_energy) - mf1.mo_energy = numpy.array((mo_ea, mo_eb)) - mf1.mo_coeff = numpy.array((mf.mo_coeff, mf.mo_coeff)) - else: # This to handle KRHF object - mf1.mo_occ = numpy.array([ - [numpy.asarray(occ> 0, dtype=numpy.double) for occ in mf.mo_occ], - [numpy.asarray(occ==2, dtype=numpy.double) for occ in mf.mo_occ]]) - mo_ea = getattr(mf.mo_energy, 'mo_ea', mf.mo_energy) - mo_eb = getattr(mf.mo_energy, 'mo_eb', mf.mo_energy) - mf1.mo_energy = numpy.array((mo_ea, mo_eb)) - mf1.mo_coeff = numpy.array((mf.mo_coeff, mf.mo_coeff)) + if mf.mo_energy is None: + return mf1 + + if isinstance(mf, scf.uhf.UHF): + mf1.mo_occ = mf.mo_occ + mf1.mo_coeff = mf.mo_coeff + mf1.mo_energy = mf.mo_energy + elif getattr(mf, 'kpts', None) is None: # RHF/ROHF + mf1.mo_occ = numpy.array((mf.mo_occ>0, mf.mo_occ==2), dtype=numpy.double) + # ROHF orbital energies, not canonical UHF orbital energies + mo_ea = getattr(mf.mo_energy, 'mo_ea', mf.mo_energy) + mo_eb = getattr(mf.mo_energy, 'mo_eb', mf.mo_energy) + mf1.mo_energy = numpy.array((mo_ea, mo_eb)) + mf1.mo_coeff = numpy.array((mf.mo_coeff, mf.mo_coeff)) + else: # This to handle KRHF object + mf1.mo_occ = numpy.array([ + [numpy.asarray(occ> 0, dtype=numpy.double) for occ in mf.mo_occ], + [numpy.asarray(occ==2, dtype=numpy.double) for occ in mf.mo_occ]]) + mo_ea = getattr(mf.mo_energy, 'mo_ea', mf.mo_energy) + mo_eb = getattr(mf.mo_energy, 'mo_eb', mf.mo_energy) + mf1.mo_energy = numpy.array((mo_ea, mo_eb)) + mf1.mo_coeff = numpy.array((mf.mo_coeff, mf.mo_coeff)) return mf1 def _update_mo_to_rhf_(mf, mf1): from pyscf import scf - if mf.mo_energy is not None: - if isinstance(mf, scf.hf.RHF): # RHF/ROHF/KRHF/KROHF - mf1.mo_occ = mf.mo_occ - mf1.mo_coeff = mf.mo_coeff - mf1.mo_energy = mf.mo_energy - elif getattr(mf, 'kpts', None) is None: # UHF - mf1.mo_occ = mf.mo_occ[0] + mf.mo_occ[1] - mf1.mo_energy = mf.mo_energy[0] - mf1.mo_coeff = mf.mo_coeff[0] - if getattr(mf.mo_coeff[0], 'orbsym', None) is not None: - mf1.mo_coeff = lib.tag_array(mf1.mo_coeff, orbsym=mf.mo_coeff[0].orbsym) - mf1.converged = False - else: # KUHF - mf1.mo_occ = [occa+occb for occa, occb in zip(*mf.mo_occ)] - mf1.mo_energy = mf.mo_energy[0] - mf1.mo_coeff = mf.mo_coeff[0] - mf1.converged = False + if mf.mo_energy is None: + return mf1 + + if isinstance(mf, scf.hf.RHF): # RHF/ROHF/KRHF/KROHF + mf1.mo_occ = mf.mo_occ + mf1.mo_coeff = mf.mo_coeff + mf1.mo_energy = mf.mo_energy + elif getattr(mf, 'kpts', None) is None: # UHF + mf1.mo_occ = mf.mo_occ[0] + mf.mo_occ[1] + mf1.mo_energy = mf.mo_energy[0] + mf1.mo_coeff = mf.mo_coeff[0] + if getattr(mf.mo_coeff[0], 'orbsym', None) is not None: + mf1.mo_coeff = lib.tag_array(mf1.mo_coeff, orbsym=mf.mo_coeff[0].orbsym) + mf1.converged = False + else: # KUHF + mf1.mo_occ = [occa+occb for occa, occb in zip(*mf.mo_occ)] + mf1.mo_energy = mf.mo_energy[0] + mf1.mo_coeff = mf.mo_coeff[0] + mf1.converged = False return mf1 def _update_mo_to_ghf_(mf, mf1): from pyscf import scf - if mf.mo_energy is not None: - if isinstance(mf, scf.hf.RHF): # RHF - nao, nmo = mf.mo_coeff.shape - orbspin = get_ghf_orbspin(mf.mo_energy, mf.mo_occ, True) - - mf1.mo_energy = numpy.empty(nmo*2) - mf1.mo_energy[orbspin==0] = mf.mo_energy - mf1.mo_energy[orbspin==1] = mf.mo_energy - mf1.mo_occ = numpy.empty(nmo*2) - mf1.mo_occ[orbspin==0] = mf.mo_occ > 0 - mf1.mo_occ[orbspin==1] = mf.mo_occ == 2 - - mo_coeff = numpy.zeros((nao*2,nmo*2), dtype=mf.mo_coeff.dtype) - mo_coeff[:nao,orbspin==0] = mf.mo_coeff - mo_coeff[nao:,orbspin==1] = mf.mo_coeff - if getattr(mf.mo_coeff, 'orbsym', None) is not None: - orbsym = numpy.zeros_like(orbspin) - orbsym[orbspin==0] = mf.mo_coeff.orbsym - orbsym[orbspin==1] = mf.mo_coeff.orbsym - mo_coeff = lib.tag_array(mo_coeff, orbsym=orbsym) - mf1.mo_coeff = lib.tag_array(mo_coeff, orbspin=orbspin) - - else: # UHF - nao, nmo = mf.mo_coeff[0].shape - orbspin = get_ghf_orbspin(mf.mo_energy, mf.mo_occ, False) - - mf1.mo_energy = numpy.empty(nmo*2) - mf1.mo_energy[orbspin==0] = mf.mo_energy[0] - mf1.mo_energy[orbspin==1] = mf.mo_energy[1] - mf1.mo_occ = numpy.empty(nmo*2) - mf1.mo_occ[orbspin==0] = mf.mo_occ[0] - mf1.mo_occ[orbspin==1] = mf.mo_occ[1] - - mo_coeff = numpy.zeros((nao*2,nmo*2), dtype=mf.mo_coeff[0].dtype) - mo_coeff[:nao,orbspin==0] = mf.mo_coeff[0] - mo_coeff[nao:,orbspin==1] = mf.mo_coeff[1] - if getattr(mf.mo_coeff[0], 'orbsym', None) is not None: - orbsym = numpy.zeros_like(orbspin) - orbsym[orbspin==0] = mf.mo_coeff[0].orbsym - orbsym[orbspin==1] = mf.mo_coeff[1].orbsym - mo_coeff = lib.tag_array(mo_coeff, orbsym=orbsym) - mf1.mo_coeff = lib.tag_array(mo_coeff, orbspin=orbspin) + if mf.mo_energy is None: + return mf1 + + if isinstance(mf, scf.hf.RHF): # RHF + nao, nmo = mf.mo_coeff.shape + orbspin = get_ghf_orbspin(mf.mo_energy, mf.mo_occ, True) + + mf1.mo_energy = numpy.empty(nmo*2) + mf1.mo_energy[orbspin==0] = mf.mo_energy + mf1.mo_energy[orbspin==1] = mf.mo_energy + mf1.mo_occ = numpy.empty(nmo*2) + mf1.mo_occ[orbspin==0] = mf.mo_occ > 0 + mf1.mo_occ[orbspin==1] = mf.mo_occ == 2 + + mo_coeff = numpy.zeros((nao*2,nmo*2), dtype=mf.mo_coeff.dtype) + mo_coeff[:nao,orbspin==0] = mf.mo_coeff + mo_coeff[nao:,orbspin==1] = mf.mo_coeff + if getattr(mf.mo_coeff, 'orbsym', None) is not None: + orbsym = numpy.zeros_like(orbspin) + orbsym[orbspin==0] = mf.mo_coeff.orbsym + orbsym[orbspin==1] = mf.mo_coeff.orbsym + mo_coeff = lib.tag_array(mo_coeff, orbsym=orbsym) + mf1.mo_coeff = lib.tag_array(mo_coeff, orbspin=orbspin) + + else: # UHF + nao, nmo = mf.mo_coeff[0].shape + orbspin = get_ghf_orbspin(mf.mo_energy, mf.mo_occ, False) + + mf1.mo_energy = numpy.empty(nmo*2) + mf1.mo_energy[orbspin==0] = mf.mo_energy[0] + mf1.mo_energy[orbspin==1] = mf.mo_energy[1] + mf1.mo_occ = numpy.empty(nmo*2) + mf1.mo_occ[orbspin==0] = mf.mo_occ[0] + mf1.mo_occ[orbspin==1] = mf.mo_occ[1] + + mo_coeff = numpy.zeros((nao*2,nmo*2), dtype=mf.mo_coeff[0].dtype) + mo_coeff[:nao,orbspin==0] = mf.mo_coeff[0] + mo_coeff[nao:,orbspin==1] = mf.mo_coeff[1] + if getattr(mf.mo_coeff[0], 'orbsym', None) is not None: + orbsym = numpy.zeros_like(orbspin) + orbsym[orbspin==0] = mf.mo_coeff[0].orbsym + orbsym[orbspin==1] = mf.mo_coeff[1].orbsym + mo_coeff = lib.tag_array(mo_coeff, orbsym=orbsym) + mf1.mo_coeff = lib.tag_array(mo_coeff, orbspin=orbspin) return mf1 def get_ghf_orbspin(mo_energy, mo_occ, is_rhf=None): diff --git a/pyscf/scf/atom_hf.py b/pyscf/scf/atom_hf.py index 94d303f8ca..7657863b12 100644 --- a/pyscf/scf/atom_hf.py +++ b/pyscf/scf/atom_hf.py @@ -200,16 +200,3 @@ def _angular_momentum_for_each_ao(mol): p0, p1 = ao_loc[i], ao_loc[i+1] ao_ang[p0:p1] = mol.bas_angular(i) return ao_ang - - -if __name__ == '__main__': - mol = gto.Mole() - mol.verbose = 5 - mol.output = None - - mol.atom = [["N", (0. , 0., .5)], - ["N", (0. , 0.,-.5)] ] - - mol.basis = {"N": '6-31g'} - mol.build() - print(get_atm_nrhf(mol)) diff --git a/pyscf/scf/atom_ks.py b/pyscf/scf/atom_ks.py index 8d3785500e..77ad838b84 100644 --- a/pyscf/scf/atom_ks.py +++ b/pyscf/scf/atom_ks.py @@ -90,15 +90,3 @@ def __init__(self, mol, *args, **kwargs): get_grad = atom_hf.AtomSphericAverageRHF.get_grad AtomSphericAverageRKS = AtomSphAverageRKS - -if __name__ == '__main__': - mol = gto.Mole() - mol.verbose = 5 - mol.output = None - - mol.atom = [["N", (0. , 0., .5)], - ["N", (0. , 0.,-.5)] ] - - mol.basis = {"N": '6-31g'} - mol.build() - print(get_atm_nrks(mol)) diff --git a/pyscf/scf/test/test_he.py b/pyscf/scf/test/test_he.py index bc60431d22..84cb369345 100644 --- a/pyscf/scf/test/test_he.py +++ b/pyscf/scf/test/test_he.py @@ -21,6 +21,7 @@ from pyscf import scf from pyscf import gto from pyscf import lib +from pyscf.scf import atom_hf, atom_ks def setUpModule(): global mol @@ -71,6 +72,16 @@ def test_nr_uhf(self): # mol.nucmod = {1:1} # mol.build() + def test_atomic_ks(self): + mol = gto.M( + atom = [["N", (0. , 0., .5)], + ["N", (0. , 0.,-.5)] ], + basis = {"N": '6-31g'} + ) + result = atom_hf.get_atm_nrhf(mol) + self.assertAlmostEqual(result['N'][0], -53.823206125468346, 9) + result = atom_ks.get_atm_nrks(mol) + self.assertAlmostEqual(result['N'][0], -53.53518426665269, 9) if __name__ == "__main__": print("Full Tests for He") diff --git a/pyscf/sgx/sgx.py b/pyscf/sgx/sgx.py index 9dd68919c1..7f8153da3e 100644 --- a/pyscf/sgx/sgx.py +++ b/pyscf/sgx/sgx.py @@ -227,7 +227,7 @@ def method_not_implemented(self, *args, **kwargs): CASSCF = method_not_implemented scf.hf.SCF.COSX = sgx_fit -mcscf.casci.CASCI.COSX = sgx_fit +mcscf.casci.CASBase.COSX = sgx_fit def _make_opt(mol, pjs=False): diff --git a/pyscf/solvent/ddcosmo.py b/pyscf/solvent/ddcosmo.py index 93d95face6..861a9816b6 100644 --- a/pyscf/solvent/ddcosmo.py +++ b/pyscf/solvent/ddcosmo.py @@ -292,7 +292,7 @@ def ddcosmo_for_tdscf(method, solvent_obj=None, dm=None): scf.hf.SCF.ddCOSMO = scf.hf.SCF.DDCOSMO = ddcosmo_for_scf mp.mp2.MP2.ddCOSMO = mp.mp2.MP2.DDCOSMO = ddcosmo_for_post_scf ci.cisd.CISD.ddCOSMO = ci.cisd.CISD.DDCOSMO = ddcosmo_for_post_scf -cc.ccsd.CCSD.ddCOSMO = cc.ccsd.CCSD.DDCOSMO = ddcosmo_for_post_scf +cc.ccsd.CCSDBase.ddCOSMO = cc.ccsd.CCSDBase.DDCOSMO = ddcosmo_for_post_scf tdscf.rhf.TDBase.ddCOSMO = tdscf.rhf.TDBase.DDCOSMO = ddcosmo_for_tdscf mcscf.casci.CASCI.ddCOSMO = mcscf.casci.CASCI.DDCOSMO = ddcosmo_for_casci mcscf.mc1step.CASSCF.ddCOSMO = mcscf.mc1step.CASSCF.DDCOSMO = ddcosmo_for_casscf diff --git a/pyscf/solvent/ddpcm.py b/pyscf/solvent/ddpcm.py index c0fecff0e7..4f2607f61f 100644 --- a/pyscf/solvent/ddpcm.py +++ b/pyscf/solvent/ddpcm.py @@ -89,7 +89,7 @@ def ddpcm_for_tdscf(method, solvent_obj=None, dm=None): scf.hf.SCF.ddPCM = scf.hf.SCF.DDPCM = ddpcm_for_scf mp.mp2.MP2.ddPCM = mp.mp2.MP2.DDPCM = ddpcm_for_post_scf ci.cisd.CISD.ddPCM = ci.cisd.CISD.DDPCM = ddpcm_for_post_scf -cc.ccsd.CCSD.ddPCM = cc.ccsd.CCSD.DDPCM = ddpcm_for_post_scf +cc.ccsd.CCSDBase.ddPCM = cc.ccsd.CCSDBase.DDPCM = ddpcm_for_post_scf tdscf.rhf.TDBase.ddPCM = tdscf.rhf.TDBase.DDPCM = ddpcm_for_tdscf mcscf.casci.CASCI.ddPCM = mcscf.casci.CASCI.DDPCM = ddpcm_for_casci mcscf.mc1step.CASSCF.ddPCM = mcscf.mc1step.CASSCF.DDPCM = ddpcm_for_casscf