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..6a7faf8519 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(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/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_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_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()