Source code for pyscf.cc.ccsd

#!/usr/bin/env python
# Copyright 2014-2021 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.
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

'''
RCCSD for real integrals
8-fold permutation symmetry has been used
(ij|kl) = (ji|kl) = (kl|ij) = ...
'''


import ctypes
from functools import reduce
import numpy
from pyscf import gto
from pyscf import lib
from pyscf.lib import logger
from pyscf import ao2mo
from pyscf.ao2mo import _ao2mo
from pyscf.cc import _ccsd
from pyscf.mp.mp2 import get_nocc, get_nmo, get_frozen_mask, get_e_hf, _mo_without_core
from pyscf import __config__

BLKMIN = getattr(__config__, 'cc_ccsd_blkmin', 4)
MEMORYMIN = getattr(__config__, 'cc_ccsd_memorymin', 2000)


# t1: ia
# t2: ijab
[docs] def kernel(mycc, eris=None, t1=None, t2=None, max_cycle=50, tol=1e-8, tolnormt=1e-6, verbose=None, callback=None): log = logger.new_logger(mycc, verbose) if eris is None: eris = mycc.ao2mo(mycc.mo_coeff) if t1 is None and t2 is None: t1, t2 = mycc.get_init_guess(eris) 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(%s) = %.15g', name, eccsd) if isinstance(mycc.diis, lib.diis.DIIS): adiis = mycc.diis elif mycc.diis: adiis = lib.diis.DIIS(mycc, mycc.diis_file, incore=mycc.incore_complete) adiis.space = mycc.diis_space else: adiis = None conv = False for istep in range(max_cycle): t1new, t2new = mycc.update_amps(t1, t2, eris) if callback is not None: callback(locals()) tmpvec = mycc.amplitudes_to_vector(t1new, t2new) tmpvec -= mycc.amplitudes_to_vector(t1, t2) normt = numpy.linalg.norm(tmpvec) tmpvec = None if mycc.iterative_damping < 1.0: alpha = mycc.iterative_damping t1new = (1-alpha) * t1 + alpha * t1new t2new *= alpha t2new += (1-alpha) * t2 t1, t2 = t1new, t2new 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(%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(name, *cput0) return conv, eccsd, t1, t2
[docs] def update_amps(mycc, t1, t2, eris): if mycc.cc2: raise NotImplementedError assert (isinstance(eris, _ChemistsERIs)) time0 = logger.process_clock(), logger.perf_counter() log = logger.Logger(mycc.stdout, mycc.verbose) nocc, nvir = t1.shape fock = eris.fock mo_e_o = eris.mo_energy[:nocc] mo_e_v = eris.mo_energy[nocc:] + mycc.level_shift t1new = numpy.zeros_like(t1) t2new = mycc._add_vvvv(t1, t2, eris, t2sym='jiba') t2new *= .5 # *.5 because t2+t2.transpose(1,0,3,2) in the end time1 = log.timer_debug1('vvvv', *time0) #** make_inter_F fov = fock[:nocc,nocc:].copy() t1new += fov foo = fock[:nocc,:nocc] - numpy.diag(mo_e_o) foo += .5 * numpy.einsum('ia,ja->ij', fock[:nocc,nocc:], t1) fvv = fock[nocc:,nocc:] - numpy.diag(mo_e_v) fvv -= .5 * numpy.einsum('ia,ib->ab', t1, fock[:nocc,nocc:]) if mycc.incore_complete: fswap = None else: fswap = lib.H5TmpFile() fwVOov, fwVooV = _add_ovvv_(mycc, t1, t2, eris, fvv, t1new, t2new, fswap) time1 = log.timer_debug1('ovvv', *time1) woooo = numpy.asarray(eris.oooo).transpose(0,2,1,3).copy() unit = nocc**2*nvir*7 + nocc**3 + nocc*nvir**2 mem_now = lib.current_memory()[0] max_memory = max(0, mycc.max_memory - mem_now) blksize = min(nvir, max(BLKMIN, int((max_memory*.9e6/8-nocc**4)/unit))) log.debug1('max_memory %d MB, nocc,nvir = %d,%d blksize = %d', max_memory, nocc, nvir, blksize) for p0, p1 in lib.prange(0, nvir, blksize): wVOov = fwVOov[p0:p1] wVooV = fwVooV[p0:p1] eris_ovoo = eris.ovoo[:,p0:p1] eris_oovv = numpy.empty((nocc,nocc,p1-p0,nvir)) def load_oovv(p0, p1): eris_oovv[:] = eris.oovv[:,:,p0:p1] with lib.call_in_background(load_oovv, sync=not mycc.async_io) as prefetch_oovv: #:eris_oovv = eris.oovv[:,:,p0:p1] prefetch_oovv(p0, p1) foo += numpy.einsum('kc,kcji->ij', 2*t1[:,p0:p1], eris_ovoo) foo += numpy.einsum('kc,icjk->ij', -t1[:,p0:p1], eris_ovoo) tmp = lib.einsum('la,jaik->lkji', t1[:,p0:p1], eris_ovoo) woooo += tmp + tmp.transpose(1,0,3,2) tmp = None wVOov -= lib.einsum('jbik,ka->bjia', eris_ovoo, t1) t2new[:,:,p0:p1] += wVOov.transpose(1,2,0,3) wVooV += lib.einsum('kbij,ka->bija', eris_ovoo, t1) eris_ovoo = None load_oovv = prefetch_oovv = None eris_ovvo = numpy.empty((nocc,p1-p0,nvir,nocc)) def load_ovvo(p0, p1): eris_ovvo[:] = eris.ovvo[:,p0:p1] with lib.call_in_background(load_ovvo, sync=not mycc.async_io) as prefetch_ovvo: #:eris_ovvo = eris.ovvo[:,p0:p1] prefetch_ovvo(p0, p1) t1new[:,p0:p1] -= numpy.einsum('jb,jiab->ia', t1, eris_oovv) wVooV -= eris_oovv.transpose(2,0,1,3) wVOov += wVooV*.5 #: bjia + bija*.5 load_ovvo = prefetch_ovvo = None t2new[:,:,p0:p1] += (eris_ovvo*0.5).transpose(0,3,1,2) eris_voov = eris_ovvo.conj().transpose(1,0,3,2) t1new[:,p0:p1] += 2*numpy.einsum('jb,aijb->ia', t1, eris_voov) eris_ovvo = None tmp = lib.einsum('ic,kjbc->ibkj', t1, eris_oovv) tmp += lib.einsum('bjkc,ic->jbki', eris_voov, t1) t2new[:,:,p0:p1] -= lib.einsum('ka,jbki->jiba', t1, tmp) eris_oovv = tmp = None fov[:,p0:p1] += numpy.einsum('kc,aikc->ia', t1, eris_voov) * 2 fov[:,p0:p1] -= numpy.einsum('kc,akic->ia', t1, eris_voov) tau = numpy.einsum('ia,jb->ijab', t1[:,p0:p1]*.5, t1) tau += t2[:,:,p0:p1] theta = tau.transpose(1,0,2,3) * 2 theta -= tau fvv -= lib.einsum('cjia,cjib->ab', theta.transpose(2,1,0,3), eris_voov) foo += lib.einsum('aikb,kjab->ij', eris_voov, theta) tau = theta = None tau = t2[:,:,p0:p1] + numpy.einsum('ia,jb->ijab', t1[:,p0:p1], t1) woooo += lib.einsum('ijab,aklb->ijkl', tau, eris_voov) tau = None def update_wVooV(q0, q1, tau): wVooV[:] += lib.einsum('bkic,jkca->bija', eris_voov[:,:,:,q0:q1], tau) with lib.call_in_background(update_wVooV, sync=not mycc.async_io) as update_wVooV: for q0, q1 in lib.prange(0, nvir, blksize): tau = t2[:,:,q0:q1] * .5 tau += numpy.einsum('ia,jb->ijab', t1[:,q0:q1], t1) #:wVooV += lib.einsum('bkic,jkca->bija', eris_voov[:,:,:,q0:q1], tau) update_wVooV(q0, q1, tau) tau = update_wVooV = None def update_t2(q0, q1, tmp): t2new[:,:,q0:q1] += tmp.transpose(2,0,1,3) tmp *= .5 t2new[:,:,q0:q1] += tmp.transpose(0,2,1,3) with lib.call_in_background(update_t2, sync=not mycc.async_io) as update_t2: for q0, q1 in lib.prange(0, nvir, blksize): tmp = lib.einsum('jkca,ckib->jaib', t2[:,:,p0:p1,q0:q1], wVooV) #:t2new[:,:,q0:q1] += tmp.transpose(2,0,1,3) #:tmp *= .5 #:t2new[:,:,q0:q1] += tmp.transpose(0,2,1,3) update_t2(q0, q1, tmp) tmp = None wVOov += eris_voov eris_VOov = -.5 * eris_voov.transpose(0,2,1,3) eris_VOov += eris_voov eris_voov = None def update_wVOov(q0, q1, tau): wVOov[:,:,:,q0:q1] += .5 * lib.einsum('aikc,kcjb->aijb', eris_VOov, tau) with lib.call_in_background(update_wVOov, sync=not mycc.async_io) as update_wVOov: for q0, q1 in lib.prange(0, nvir, blksize): tau = t2[:,:,q0:q1].transpose(1,3,0,2) * 2 tau -= t2[:,:,q0:q1].transpose(0,3,1,2) tau -= numpy.einsum('ia,jb->ibja', t1[:,q0:q1]*2, t1) #:wVOov[:,:,:,q0:q1] += .5 * lib.einsum('aikc,kcjb->aijb', eris_VOov, tau) update_wVOov(q0, q1, tau) tau = None def update_t2(q0, q1, theta): t2new[:,:,q0:q1] += lib.einsum('kica,ckjb->ijab', theta, wVOov) with lib.call_in_background(update_t2, sync=not mycc.async_io) as update_t2: for q0, q1 in lib.prange(0, nvir, blksize): theta = t2[:,:,p0:p1,q0:q1] * 2 theta -= t2[:,:,p0:p1,q0:q1].transpose(1,0,2,3) #:t2new[:,:,q0:q1] += lib.einsum('kica,ckjb->ijab', theta, wVOov) update_t2(q0, q1, theta) theta = None eris_VOov = wVOov = wVooV = update_wVOov = None time1 = log.timer_debug1('voov [%d:%d]'%(p0, p1), *time1) fwVOov = fwVooV = fswap = None for p0, p1 in lib.prange(0, nvir, blksize): theta = t2[:,:,p0:p1].transpose(1,0,2,3) * 2 - t2[:,:,p0:p1] t1new += numpy.einsum('jb,ijba->ia', fov[:,p0:p1], theta) t1new -= lib.einsum('jbki,kjba->ia', eris.ovoo[:,p0:p1], theta) tau = numpy.einsum('ia,jb->ijab', t1[:,p0:p1], t1) tau += t2[:,:,p0:p1] t2new[:,:,p0:p1] += .5 * lib.einsum('ijkl,klab->ijab', woooo, tau) theta = tau = None ft_ij = foo + numpy.einsum('ja,ia->ij', .5*t1, fov) ft_ab = fvv - numpy.einsum('ia,ib->ab', .5*t1, fov) t2new += lib.einsum('ijac,bc->ijab', t2, ft_ab) t2new -= lib.einsum('ki,kjab->ijab', ft_ij, t2) eia = mo_e_o[:,None] - mo_e_v t1new += numpy.einsum('ib,ab->ia', t1, fvv) t1new -= numpy.einsum('ja,ji->ia', t1, foo) t1new /= eia #: t2new = t2new + t2new.transpose(1,0,3,2) for i in range(nocc): if i > 0: t2new[i,:i] += t2new[:i,i].transpose(0,2,1) t2new[i,:i] /= lib.direct_sum('a,jb->jab', eia[i], eia[:i]) t2new[:i,i] = t2new[i,:i].transpose(0,2,1) t2new[i,i] = t2new[i,i] + t2new[i,i].T t2new[i,i] /= lib.direct_sum('a,b->ab', eia[i], eia[i]) time0 = log.timer_debug1('update t1 t2', *time0) return t1new, t2new
def _add_ovvv_(mycc, t1, t2, eris, fvv, t1new, t2new, fswap): time1 = logger.process_clock(), logger.perf_counter() log = logger.Logger(mycc.stdout, mycc.verbose) nocc, nvir = t1.shape nvir_pair = nvir * (nvir+1) // 2 if fswap is None: wVOov = numpy.zeros((nvir,nocc,nocc,nvir)) else: wVOov = fswap.create_dataset('wVOov', (nvir,nocc,nocc,nvir), 'f8') wooVV = numpy.zeros((nocc,nocc*nvir_pair)) max_memory = mycc.max_memory - lib.current_memory()[0] unit = nocc*nvir**2*3 + nocc**2*nvir + 2 blksize = min(nvir, max(BLKMIN, int((max_memory*.95e6/8-wooVV.size)/unit))) if not mycc.direct: unit = nocc*nvir**2*3 + nocc**2*nvir + 2 + nocc*nvir**2 + nocc*nvir blksize = min(nvir, max(BLKMIN, int((max_memory*.95e6/8-wooVV.size-nocc**2*nvir)/unit))) log.debug1('max_memory %d MB, nocc,nvir = %d,%d blksize = %d', max_memory, nocc, nvir, blksize) def load_ovvv(buf, p0): if p0 < nvir: p1 = min(nvir, p0+blksize) buf[:p1-p0] = eris.ovvv[:,p0:p1].transpose(1,0,2) with lib.call_in_background(load_ovvv, sync=not mycc.async_io) as prefetch: buf = numpy.empty((blksize,nocc,nvir_pair)) buf_prefetch = numpy.empty((blksize,nocc,nvir_pair)) load_ovvv(buf_prefetch, 0) for p0, p1 in lib.prange(0, nvir, blksize): buf, buf_prefetch = buf_prefetch, buf prefetch(buf_prefetch, p1) eris_vovv = buf[:p1-p0] #:wooVV -= numpy.einsum('jc,ciba->jiba', t1[:,p0:p1], eris_vovv) lib.ddot(numpy.asarray(t1[:,p0:p1], order='C'), eris_vovv.reshape(p1-p0,-1), -1, wooVV, 1) eris_vovv = lib.unpack_tril(eris_vovv.reshape((p1-p0)*nocc,nvir_pair)) eris_vovv = eris_vovv.reshape(p1-p0,nocc,nvir,nvir) fvv += 2*numpy.einsum('kc,ckab->ab', t1[:,p0:p1], eris_vovv) fvv[:,p0:p1] -= numpy.einsum('kc,bkca->ab', t1, eris_vovv) if not mycc.direct: vvvo = eris_vovv.transpose(0,2,3,1).copy() for i in range(nocc): tau = t2[i,:,p0:p1] + numpy.einsum('a,jb->jab', t1[i,p0:p1], t1) tmp = lib.einsum('jcd,cdbk->jbk', tau, vvvo) t2new[i] -= lib.einsum('ka,jbk->jab', t1, tmp) tau = tmp = None wVOov[p0:p1] = lib.einsum('biac,jc->bija', eris_vovv, t1) theta = t2[:,:,p0:p1].transpose(1,2,0,3) * 2 theta -= t2[:,:,p0:p1].transpose(0,2,1,3) t1new += lib.einsum('icjb,cjba->ia', theta, eris_vovv) theta = None time1 = log.timer_debug1('vovv [%d:%d]'%(p0, p1), *time1) if fswap is None: wooVV = lib.unpack_tril(wooVV.reshape(nocc**2,nvir_pair)) return wVOov, wooVV.reshape(nocc,nocc,nvir,nvir).transpose(2,1,0,3) else: fswap.create_dataset('wVooV', (nvir,nocc,nocc,nvir), 'f8') wooVV = wooVV.reshape(nocc,nocc,nvir_pair) tril2sq = lib.square_mat_in_trilu_indices(nvir) for p0, p1 in lib.prange(0, nvir, blksize): fswap['wVooV'][p0:p1] = wooVV[:,:,tril2sq[p0:p1]].transpose(2,1,0,3) return fswap['wVOov'], fswap['wVooV'] def _add_vvvv(mycc, t1, t2, eris, out=None, with_ovvv=None, t2sym=None): '''t2sym: whether t2 has the symmetry t2[ijab]==t2[jiba] or t2[ijab]==-t2[jiab] or t2[ijab]==-t2[jiba] ''' #TODO: Guess the symmetry of t2 amplitudes #if t2sym is None: # if t2.shape[0] != t2.shape[1]: # t2sym = '' # elif abs(t2-t2.transpose(1,0,3,2)).max() < 1e-12: # t2sym = 'jiba' # elif abs(t2+t2.transpose(1,0,2,3)).max() < 1e-12: # t2sym = '-jiab' # elif abs(t2+t2.transpose(1,0,3,2)).max() < 1e-12: # t2sym = '-jiba' if t2sym in ('jiba', '-jiba', '-jiab'): Ht2tril = _add_vvvv_tril(mycc, t1, t2, eris, with_ovvv=with_ovvv) nocc, nvir = t2.shape[1:3] Ht2 = _unpack_t2_tril(Ht2tril, nocc, nvir, out, t2sym) else: Ht2 = _add_vvvv_full(mycc, t1, t2, eris, out, with_ovvv) return Ht2 def _add_vvvv_tril(mycc, t1, t2, eris, out=None, with_ovvv=None): '''Ht2 = numpy.einsum('ijcd,acdb->ijab', t2, vvvv) Using symmetry t2[ijab] = t2[jiba] and Ht2[ijab] = Ht2[jiba], compute the lower triangular part of Ht2 ''' time0 = logger.process_clock(), logger.perf_counter() log = logger.Logger(mycc.stdout, mycc.verbose) if with_ovvv is None: with_ovvv = mycc.direct nocc, nvir = t2.shape[1:3] nocc2 = nocc*(nocc+1)//2 if t1 is None: tau = t2[numpy.tril_indices(nocc)] else: tau = numpy.empty((nocc2,nvir,nvir), dtype=t2.dtype) p1 = 0 for i in range(nocc): p0, p1 = p1, p1 + i+1 tau[p0:p1] = numpy.einsum('a,jb->jab', t1[i], t1[:i+1]) tau[p0:p1] += t2[i,:i+1] if mycc.direct: # AO-direct CCSD mo = getattr(eris, 'mo_coeff', None) if mo is None: # If eris does not have the attribute mo_coeff mo = _mo_without_core(mycc, mycc.mo_coeff) nao, nmo = mo.shape aos = numpy.asarray(mo[:,nocc:].T, order='F') tau = _ao2mo.nr_e2(tau.reshape(nocc2,nvir**2), aos, (0,nao,0,nao), 's1', 's1') tau = tau.reshape(nocc2,nao,nao) time0 = log.timer_debug1('vvvv-tau', *time0) buf = eris._contract_vvvv_t2(mycc, tau, mycc.direct, out, log) buf = buf.reshape(nocc2,nao,nao) Ht2tril = _ao2mo.nr_e2(buf, mo.conj(), (nocc,nmo,nocc,nmo), 's1', 's1') Ht2tril = Ht2tril.reshape(nocc2,nvir,nvir) if with_ovvv: #: tmp = numpy.einsum('ijcd,ka,kdcb->ijba', tau, t1, eris.ovvv) #: t2new -= tmp + tmp.transpose(1,0,3,2) tmp = _ao2mo.nr_e2(buf, mo.conj(), (nocc,nmo,0,nocc), 's1', 's1') Ht2tril -= lib.ddot(tmp.reshape(nocc2*nvir,nocc), t1).reshape(nocc2,nvir,nvir) tmp = _ao2mo.nr_e2(buf, mo.conj(), (0,nocc,nocc,nmo), 's1', 's1') #: Ht2tril -= numpy.einsum('xkb,ka->xab', tmp.reshape(-1,nocc,nvir), t1) tmp = lib.transpose(tmp.reshape(nocc2,nocc,nvir), axes=(0,2,1), out=buf) tmp = lib.ddot(tmp.reshape(nocc2*nvir,nocc), t1, 1, numpy.ndarray((nocc2*nvir,nvir), buffer=tau), 0) tmp = lib.transpose(tmp.reshape(nocc2,nvir,nvir), axes=(0,2,1), out=buf) Ht2tril -= tmp.reshape(nocc2,nvir,nvir) else: assert (not with_ovvv) Ht2tril = eris._contract_vvvv_t2(mycc, tau, mycc.direct, out, log) return Ht2tril def _add_vvvv_full(mycc, t1, t2, eris, out=None, with_ovvv=False): '''Ht2 = numpy.einsum('ijcd,acdb->ijab', t2, vvvv) without using symmetry t2[ijab] = t2[jiba] in t2 or Ht2 ''' time0 = logger.process_clock(), logger.perf_counter() log = logger.Logger(mycc.stdout, mycc.verbose) if t1 is None: tau = t2 else: tau = numpy.einsum('ia,jb->ijab', t1, t1) tau += t2 if mycc.direct: # AO-direct CCSD if with_ovvv: raise NotImplementedError mo = getattr(eris, 'mo_coeff', None) if mo is None: # If eris does not have the attribute mo_coeff mo = _mo_without_core(mycc, mycc.mo_coeff) nocc, nvir = t2.shape[1:3] nao, nmo = mo.shape aos = numpy.asarray(mo[:,nocc:].T, order='F') tau = _ao2mo.nr_e2(tau.reshape(nocc**2,nvir,nvir), aos, (0,nao,0,nao), 's1', 's1') tau = tau.reshape(nocc,nocc,nao,nao) time0 = log.timer_debug1('vvvv-tau mo2ao', *time0) buf = eris._contract_vvvv_t2(mycc, tau, mycc.direct, out, log) buf = buf.reshape(nocc**2,nao,nao) Ht2 = _ao2mo.nr_e2(buf, mo.conj(), (nocc,nmo,nocc,nmo), 's1', 's1') else: assert (not with_ovvv) Ht2 = eris._contract_vvvv_t2(mycc, tau, mycc.direct, out, log) return Ht2.reshape(t2.shape) def _contract_vvvv_t2(mycc, mol, vvvv, t2, out=None, verbose=None): '''Ht2 = numpy.einsum('ijcd,acbd->ijab', t2, vvvv) Args: vvvv : None or integral object if vvvv is None, contract t2 to AO-integrals using AO-direct algorithm ''' if vvvv is None or len(vvvv.shape) == 2: # AO-direct or vvvv in 4-fold symmetry return _contract_s4vvvv_t2(mycc, mol, vvvv, t2, out, verbose) else: return _contract_s1vvvv_t2(mycc, mol, vvvv, t2, out, verbose) def _contract_s4vvvv_t2(mycc, mol, vvvv, t2, out=None, verbose=None): '''Ht2 = numpy.einsum('ijcd,acbd->ijab', t2, vvvv) where vvvv has to be real and has the 4-fold permutation symmetry Args: vvvv : None or integral object if vvvv is None, contract t2 to AO-integrals using AO-direct algorithm ''' assert (t2.dtype == numpy.double) if t2.size == 0: return numpy.zeros_like(t2) _dgemm = lib.numpy_helper._dgemm time0 = logger.process_clock(), logger.perf_counter() log = logger.new_logger(mycc, verbose) nvira, nvirb = t2.shape[-2:] x2 = t2.reshape(-1,nvira,nvirb) nocc2 = x2.shape[0] nvir2 = nvira * nvirb Ht2 = numpy.ndarray(x2.shape, dtype=x2.dtype, buffer=out) Ht2[:] = 0 def contract_blk_(eri, i0, i1, j0, j1): ic = i1 - i0 jc = j1 - j0 #:Ht2[:,j0:j1] += numpy.einsum('xef,efab->xab', x2[:,i0:i1], eri) _dgemm('N', 'N', nocc2, jc*nvirb, ic*nvirb, x2.reshape(-1,nvir2), eri.reshape(-1,jc*nvirb), Ht2.reshape(-1,nvir2), 1, 1, i0*nvirb, 0, j0*nvirb) if i0 > j0: #:Ht2[:,i0:i1] += numpy.einsum('xef,abef->xab', x2[:,j0:j1], eri) _dgemm('N', 'T', nocc2, ic*nvirb, jc*nvirb, x2.reshape(-1,nvir2), eri.reshape(-1,jc*nvirb), Ht2.reshape(-1,nvir2), 1, 1, j0*nvirb, 0, i0*nvirb) max_memory = max(MEMORYMIN, mycc.max_memory - lib.current_memory()[0]) if vvvv is None: # AO-direct CCSD ao_loc = mol.ao_loc_nr() assert (nvira == nvirb == ao_loc[-1]) intor = mol._add_suffix('int2e') ao2mopt = _ao2mo.AO2MOpt(mol, intor, 'CVHFnr_schwarz_cond', 'CVHFsetnr_direct_scf') blksize = max(BLKMIN, numpy.sqrt(max_memory*.9e6/8/nvirb**2/2.5)) blksize = int(min((nvira+3)/4, blksize)) sh_ranges = ao2mo.outcore.balance_partition(ao_loc, blksize) blksize = max(x[2] for x in sh_ranges) eribuf = numpy.empty((blksize,blksize,nvirb,nvirb)) loadbuf = numpy.empty((blksize,blksize,nvirb,nvirb)) fint = gto.moleintor.getints4c for ip, (ish0, ish1, ni) in enumerate(sh_ranges): for jsh0, jsh1, nj in sh_ranges[:ip]: eri = fint(intor, mol._atm, mol._bas, mol._env, shls_slice=(ish0,ish1,jsh0,jsh1), aosym='s2kl', ao_loc=ao_loc, cintopt=ao2mopt._cintopt, out=eribuf) i0, i1 = ao_loc[ish0], ao_loc[ish1] j0, j1 = ao_loc[jsh0], ao_loc[jsh1] tmp = numpy.ndarray((i1-i0,nvirb,j1-j0,nvirb), buffer=loadbuf) _ccsd.libcc.CCload_eri(tmp.ctypes.data_as(ctypes.c_void_p), eri.ctypes.data_as(ctypes.c_void_p), (ctypes.c_int*4)(i0, i1, j0, j1), ctypes.c_int(nvirb)) contract_blk_(tmp, i0, i1, j0, j1) time0 = log.timer_debug1('AO-vvvv [%d:%d,%d:%d]' % (ish0,ish1,jsh0,jsh1), *time0) eri = fint(intor, mol._atm, mol._bas, mol._env, shls_slice=(ish0,ish1,ish0,ish1), aosym='s4', ao_loc=ao_loc, cintopt=ao2mopt._cintopt, out=eribuf) i0, i1 = ao_loc[ish0], ao_loc[ish1] eri = lib.unpack_tril(eri, axis=0) tmp = numpy.ndarray((i1-i0,nvirb,i1-i0,nvirb), buffer=loadbuf) _ccsd.libcc.CCload_eri(tmp.ctypes.data_as(ctypes.c_void_p), eri.ctypes.data_as(ctypes.c_void_p), (ctypes.c_int*4)(i0, i1, i0, i1), ctypes.c_int(nvirb)) eri = None contract_blk_(tmp, i0, i1, i0, i1) time0 = log.timer_debug1('AO-vvvv [%d:%d,%d:%d]' % (ish0,ish1,ish0,ish1), *time0) else: nvir_pair = nvirb * (nvirb+1) // 2 unit = nvira*nvir_pair*2 + nvirb**2*nvira/4 + 1 if mycc.async_io: fmap = lib.map_with_prefetch unit += nvira*nvir_pair else: fmap = map blksize = numpy.sqrt(max(BLKMIN**2, max_memory*.95e6/8/unit)) blksize = int(min((nvira+3)/4, blksize)) def load(v_slice): i0, i1 = v_slice off0 = i0*(i0+1)//2 off1 = i1*(i1+1)//2 return numpy.asarray(vvvv[off0:off1], order='C') tril2sq = lib.square_mat_in_trilu_indices(nvira) loadbuf = numpy.empty((blksize,blksize,nvirb,nvirb)) slices = [(i0, i1) for i0, i1 in lib.prange(0, nvira, blksize)] for istep, wwbuf in enumerate(fmap(load, lib.prange(0, nvira, blksize))): i0, i1 = slices[istep] off0 = i0*(i0+1)//2 for j0, j1 in lib.prange(0, i1, blksize): eri = wwbuf[tril2sq[i0:i1,j0:j1]-off0] tmp = numpy.ndarray((i1-i0,nvirb,j1-j0,nvirb), buffer=loadbuf) _ccsd.libcc.CCload_eri(tmp.ctypes.data_as(ctypes.c_void_p), eri.ctypes.data_as(ctypes.c_void_p), (ctypes.c_int*4)(i0, i1, j0, j1), ctypes.c_int(nvirb)) contract_blk_(tmp, i0, i1, j0, j1) wwbuf = None time0 = log.timer_debug1('vvvv [%d:%d]'%(i0,i1), *time0) return Ht2.reshape(t2.shape) def _contract_s1vvvv_t2(mycc, mol, vvvv, t2, out=None, verbose=None): '''Ht2 = numpy.einsum('ijcd,acdb->ijab', t2, vvvv) where vvvv can be real or complex and no permutation symmetry is available in vvvv. Args: vvvv : None or integral object if vvvv is None, contract t2 to AO-integrals using AO-direct algorithm ''' # vvvv == None means AO-direct CCSD. It should redirect to # _contract_s4vvvv_t2(mycc, mol, vvvv, t2, out, verbose) assert (vvvv is not None) time0 = logger.process_clock(), logger.perf_counter() log = logger.new_logger(mycc, verbose) nvira, nvirb = t2.shape[-2:] x2 = t2.reshape(-1,nvira,nvirb) nocc2 = x2.shape[0] dtype = numpy.result_type(t2, vvvv) Ht2 = numpy.ndarray(x2.shape, dtype=dtype, buffer=out) max_memory = mycc.max_memory - lib.current_memory()[0] unit = nvirb**2*nvira*2 + nocc2*nvirb + 1 blksize = min(nvira, max(BLKMIN, int(max_memory*1e6/8/unit))) for p0,p1 in lib.prange(0, nvira, blksize): Ht2[:,p0:p1] = lib.einsum('xcd,acbd->xab', x2, vvvv[p0:p1]) time0 = log.timer_debug1('vvvv [%d:%d]' % (p0,p1), *time0) return Ht2.reshape(t2.shape) def _unpack_t2_tril(t2tril, nocc, nvir, out=None, t2sym='jiba'): t2 = numpy.ndarray((nocc,nocc,nvir,nvir), dtype=t2tril.dtype, buffer=out) idx,idy = numpy.tril_indices(nocc) if t2sym == 'jiba': t2[idy,idx] = t2tril.transpose(0,2,1) t2[idx,idy] = t2tril elif t2sym == '-jiba': t2[idy,idx] = -t2tril.transpose(0,2,1) t2[idx,idy] = t2tril elif t2sym == '-jiab': t2[idy,idx] =-t2tril t2[idx,idy] = t2tril t2[numpy.diag_indices(nocc)] = 0 return t2 def _unpack_4fold(c2vec, nocc, nvir, anti_symm=True): t2 = numpy.zeros((nocc**2,nvir**2), dtype=c2vec.dtype) if nocc > 1 and nvir > 1: t2tril = c2vec.reshape(nocc*(nocc-1)//2,nvir*(nvir-1)//2) otril = numpy.tril_indices(nocc, k=-1) vtril = numpy.tril_indices(nvir, k=-1) lib.takebak_2d(t2, t2tril, otril[0]*nocc+otril[1], vtril[0]*nvir+vtril[1]) lib.takebak_2d(t2, t2tril, otril[1]*nocc+otril[0], vtril[1]*nvir+vtril[0]) if anti_symm: # anti-symmetry when exchanging two particle indices t2tril = -t2tril lib.takebak_2d(t2, t2tril, otril[0]*nocc+otril[1], vtril[1]*nvir+vtril[0]) lib.takebak_2d(t2, t2tril, otril[1]*nocc+otril[0], vtril[0]*nvir+vtril[1]) return t2.reshape(nocc,nocc,nvir,nvir)
[docs] def amplitudes_to_vector(t1, t2, out=None): nocc, nvir = t1.shape nov = nocc * nvir size = nov + nov*(nov+1)//2 vector = numpy.ndarray(size, t1.dtype, buffer=out) vector[:nov] = t1.ravel() lib.pack_tril(t2.transpose(0,2,1,3).reshape(nov,nov), out=vector[nov:]) return vector
[docs] def vector_to_amplitudes(vector, nmo, nocc): nvir = nmo - nocc nov = nocc * nvir t1 = vector[:nov].copy().reshape((nocc,nvir)) # filltriu=lib.SYMMETRIC because t2[iajb] == t2[jbia] t2 = lib.unpack_tril(vector[nov:], filltriu=lib.SYMMETRIC) t2 = t2.reshape(nocc,nvir,nocc,nvir).transpose(0,2,1,3) return t1, numpy.asarray(t2, order='C')
[docs] def amplitudes_to_vector_s4(t1, t2, out=None): nocc, nvir = t1.shape nov = nocc * nvir size = nov + nocc*(nocc-1)//2*nvir*(nvir-1)//2 vector = numpy.ndarray(size, t1.dtype, buffer=out) vector[:nov] = t1.ravel() otril = numpy.tril_indices(nocc, k=-1) vtril = numpy.tril_indices(nvir, k=-1) lib.take_2d(t2.reshape(nocc**2,nvir**2), otril[0]*nocc+otril[1], vtril[0]*nvir+vtril[1], out=vector[nov:]) return vector
[docs] def vector_to_amplitudes_s4(vector, nmo, nocc): nvir = nmo - nocc nov = nocc * nvir size = nov + nocc*(nocc-1)//2*nvir*(nvir-1)//2 t1 = vector[:nov].copy().reshape(nocc,nvir) t2 = numpy.zeros((nocc,nocc,nvir,nvir), dtype=vector.dtype) t2 = _unpack_4fold(vector[nov:size], nocc, nvir) return t1, t2
[docs] def energy(mycc, t1=None, t2=None, eris=None): '''CCSD correlation energy''' if t1 is None: t1 = mycc.t1 if t2 is None: t2 = mycc.t2 if eris is None: eris = mycc.ao2mo() nocc, nvir = t1.shape fock = eris.fock e = numpy.einsum('ia,ia', fock[:nocc,nocc:], t1) * 2 max_memory = mycc.max_memory - lib.current_memory()[0] blksize = int(min(nvir, max(BLKMIN, max_memory*.3e6/8/(nocc**2*nvir+1)))) for p0, p1 in lib.prange(0, nvir, blksize): eris_ovvo = eris.ovvo[:,p0:p1] tau = t2[:,:,p0:p1] + numpy.einsum('ia,jb->ijab', t1[:,p0:p1], t1) 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 %s energy %s', mycc.__class__.__name__, e) return e.real
[docs] def restore_from_diis_(mycc, diis_file, inplace=True): '''Reuse an existed DIIS object in the CCSD calculation. The CCSD amplitudes will be restored from the DIIS object to generate t1 and t2 amplitudes. The t1/t2 amplitudes of the CCSD object will be overwritten by the generated t1 and t2 amplitudes. The amplitudes vector and error vector will be reused in the CCSD calculation. ''' adiis = lib.diis.DIIS(mycc, mycc.diis_file, incore=mycc.incore_complete) adiis.restore(diis_file, inplace=inplace) ccvec = adiis.extrapolate() mycc.t1, mycc.t2 = mycc.vector_to_amplitudes(ccvec) if inplace: mycc.diis = adiis return mycc
[docs] def get_t1_diagnostic(t1): '''Returns the t1 amplitude norm, normalized by number of correlated electrons.''' nelectron = 2 * t1.shape[0] return numpy.sqrt(numpy.linalg.norm(t1)**2 / nelectron)
[docs] def get_d1_diagnostic(t1): '''D1 diagnostic given in Janssen, et. al Chem. Phys. Lett. 290 (1998) 423 ''' f = lambda x: numpy.sqrt(numpy.sort(numpy.abs(x[0])))[-1] d1norm_ij = f(numpy.linalg.eigh(numpy.einsum('ia,ja->ij',t1,t1))) d1norm_ab = f(numpy.linalg.eigh(numpy.einsum('ia,ib->ab',t1,t1))) d1norm = max(d1norm_ij, d1norm_ab) return d1norm
[docs] def get_d2_diagnostic(t2): '''D2 diagnostic given in Nielsen, et. al Chem. Phys. Lett. 310 (1999) 568 Note: This is currently only defined in the literature for restricted closed-shell systems. ''' f = lambda x: numpy.sqrt(numpy.sort(numpy.abs(x[0])))[-1] d2norm_ij = f(numpy.linalg.eigh(numpy.einsum('ikab,jkab->ij',t2,t2))) d2norm_ab = f(numpy.linalg.eigh(numpy.einsum('ijac,ijbc->ab',t2,t2))) d2norm = max(d2norm_ij, d2norm_ab) return d2norm
[docs] def set_frozen(mycc, method='auto', window=(-1000.0, 1000.0), is_gcc=False): if method == 'auto': from pyscf.data import elements mycc.frozen = elements.chemcore(mycc.mol, spinorb=is_gcc) elif method == 'window': emin, emax = window mo_e = numpy.asarray(mycc._scf.mo_energy) if mo_e.ndim == 1: fr1 = list(numpy.flatnonzero(mo_e < emin)) fr2 = list(numpy.flatnonzero(mo_e > emax)) frozen = fr1 + fr2 elif mo_e.ndim == 2: fr1a = list(numpy.flatnonzero(mo_e[0] < emin)) fr2a = list(numpy.flatnonzero(mo_e[0] > emax)) fr1b = list(numpy.flatnonzero(mo_e[1] < emin)) fr2b = list(numpy.flatnonzero(mo_e[1] > emax)) frozen = [fr1a+fr2a, fr1b+fr2b] mycc.frozen = frozen return mycc
[docs] def as_scanner(cc): '''Generating a scanner/solver for CCSD PES. The returned solver is a function. This function requires one argument "mol" as input and returns total CCSD energy. The solver will automatically use the results of last calculation as the initial guess of the new calculation. All parameters assigned in the CCSD and the underlying SCF objects (conv_tol, max_memory etc) are automatically applied in the solver. Note scanner has side effects. It may change many underlying objects (_scf, with_df, with_x2c, ...) during calculation. Examples:: >>> from pyscf import gto, scf, cc >>> mol = gto.M(atom='H 0 0 0; F 0 0 1') >>> cc_scanner = cc.CCSD(scf.RHF(mol)).as_scanner() >>> e_tot = cc_scanner(gto.M(atom='H 0 0 0; F 0 0 1.1')) >>> e_tot = cc_scanner(gto.M(atom='H 0 0 0; F 0 0 1.5')) ''' if isinstance(cc, lib.SinglePointScanner): return cc logger.info(cc, 'Set %s as a scanner', cc.__class__) name = cc.__class__.__name__ + CCSD_Scanner.__name_mixin__ return lib.set_class(CCSD_Scanner(cc), (CCSD_Scanner, cc.__class__), name)
[docs] class CCSD_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.MoleBase): 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
[docs] class CCSDBase(lib.StreamObject): '''restricted CCSD Attributes: verbose : int Print level. Default value equals to :class:`Mole.verbose` max_memory : float or int Allowed memory in MB. Default value equals to :class:`Mole.max_memory` conv_tol : float converge threshold. Default is 1e-7. conv_tol_normt : float converge threshold for norm(t1,t2). Default is 1e-5. max_cycle : int max number of iterations. Default is 50. diis_space : int DIIS space size. Default is 6. diis_start_cycle : int The step to start DIIS. Default is 0. iterative_damping : float The self consistent damping parameter. direct : bool AO-direct CCSD. Default is False. async_io : bool Allow for asynchronous function execution. Default is True. incore_complete : bool Avoid all I/O (also for DIIS). Default is False. level_shift : float A shift on virtual orbital energies to stablize the CCSD iteration frozen : int or list If integer is given, the inner-most orbitals are frozen from CC amplitudes. Given the orbital indices (0-based) in a list, both occupied and virtual orbitals can be frozen in CC calculation. >>> mol = gto.M(atom = 'H 0 0 0; F 0 0 1.1', basis = 'ccpvdz') >>> mf = scf.RHF(mol).run() >>> # freeze 2 core orbitals >>> mycc = cc.CCSD(mf).set(frozen = 2).run() >>> # auto-generate the number of core orbitals to be frozen (1 in this case) >>> mycc = cc.CCSD(mf).set_frozen().run() >>> # freeze 2 core orbitals and 3 high lying unoccupied orbitals >>> mycc.set(frozen = [0,1,16,17,18]).run() callback : function(envs_dict) => None callback function takes one dict as the argument which is generated by the builtin function :func:`locals`, so that the callback function can access all local variables in the current environment. Saved results converged : bool CCSD converged or not e_corr : float CCSD correlation correction e_tot : float Total CCSD energy (HF + correlation) t1, t2 : T amplitudes t1[i,a], t2[i,j,a,b] (i,j in occ, a,b in virt) l1, l2 : Lambda amplitudes l1[i,a], l2[i,j,a,b] (i,j in occ, a,b in virt) ''' max_cycle = getattr(__config__, 'cc_ccsd_CCSD_max_cycle', 50) conv_tol = getattr(__config__, 'cc_ccsd_CCSD_conv_tol', 1e-7) iterative_damping = getattr(__config__, 'cc_ccsd_CCSD_iterative_damping', 1.0) conv_tol_normt = getattr(__config__, 'cc_ccsd_CCSD_conv_tol_normt', 1e-5) diis = getattr(__config__, 'cc_ccsd_CCSD_diis', True) diis_space = getattr(__config__, 'cc_ccsd_CCSD_diis_space', 6) diis_file = None diis_start_cycle = getattr(__config__, 'cc_ccsd_CCSD_diis_start_cycle', 0) # FIXME: Should we avoid DIIS starting early? diis_start_energy_diff = getattr(__config__, 'cc_ccsd_CCSD_diis_start_energy_diff', 1e9) direct = getattr(__config__, 'cc_ccsd_CCSD_direct', False) 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', 'conv_tol_normt', 'diis', 'diis_space', 'diis_file', 'diis_start_cycle', 'diis_start_energy_diff', 'direct', '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', )) def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None): from pyscf.scf import hf if isinstance(mf, hf.KohnShamDFT): raise RuntimeError( 'CCSD Warning: The first argument mf is a DFT object. ' 'CCSD calculation should be initialized with HF object.\n' 'DFT can be converted to HF object with the mf.to_hf() method\n') if mo_coeff is None: mo_coeff = mf.mo_coeff if mo_occ is None: mo_occ = mf.mo_occ self.mol = mf.mol self._scf = mf self.verbose = self.mol.verbose self.stdout = self.mol.stdout self.max_memory = mf.max_memory self.frozen = frozen self.incore_complete = self.incore_complete or self.mol.incore_anyway self.level_shift = 0 ################################################## # don't modify the following attributes, they are not input options self.mo_coeff = mo_coeff self.mo_occ = mo_occ self.converged = False self.converged_lambda = False self.emp2 = None self.e_hf = None self.e_corr = None self.t1 = None self.t2 = None self.l1 = None self.l2 = None self._nocc = None self._nmo = None self.chkfile = mf.chkfile @property def ecc(self): return self.e_corr @property def e_tot(self): return self.e_hf + self.e_corr @property def nocc(self): return self.get_nocc() @nocc.setter def nocc(self, n): self._nocc = n @property def nmo(self): return self.get_nmo() @nmo.setter def nmo(self, n): self._nmo = n
[docs] def reset(self, mol=None): if mol is not None: self.mol = mol self._scf.reset(mol) return self
get_nocc = get_nocc get_nmo = get_nmo get_frozen_mask = get_frozen_mask get_e_hf = get_e_hf
[docs] def set_frozen(self, method='auto', window=(-1000.0, 1000.0)): from pyscf import cc is_gcc = isinstance(self, cc.gccsd.GCCSD) return set_frozen(self, method=method, window=window, is_gcc=is_gcc)
[docs] def dump_flags(self, verbose=None): log = logger.new_logger(self, verbose) log.info('') log.info('******** %s ********', self.__class__) log.info('CC2 = %g', self.cc2) log.info('CCSD nocc = %s, nmo = %s', self.nocc, self.nmo) if self.frozen is not None: log.info('frozen orbitals %s', self.frozen) log.info('max_cycle = %d', self.max_cycle) log.info('direct = %d', self.direct) log.info('conv_tol = %g', self.conv_tol) log.info('conv_tol_normt = %s', self.conv_tol_normt) log.info('diis_space = %d', self.diis_space) #log.info('diis_file = %s', self.diis_file) log.info('diis_start_cycle = %d', self.diis_start_cycle) 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]) return self
[docs] def get_init_guess(self, eris=None): return self.init_amps(eris)[1:]
[docs] def init_amps(self, eris=None): time0 = logger.process_clock(), logger.perf_counter() if eris is None: eris = self.ao2mo(self.mo_coeff) e_hf = self.e_hf if e_hf is None: e_hf = self.get_e_hf(mo_coeff=self.mo_coeff) mo_e = eris.mo_energy nocc = self.nocc nvir = mo_e.size - nocc eia = mo_e[:nocc,None] - mo_e[None,nocc:] t1 = eris.fock[:nocc,nocc:] / eia t2 = numpy.empty((nocc,nocc,nvir,nvir), dtype=eris.ovov.dtype) max_memory = self.max_memory - lib.current_memory()[0] blksize = int(min(nvir, max(BLKMIN, max_memory*.3e6/8/(nocc**2*nvir+1)))) emp2 = 0 for p0, p1 in lib.prange(0, nvir, blksize): eris_ovov = eris.ovov[:,p0:p1] t2[:,:,p0:p1] = (eris_ovov.transpose(0,2,1,3).conj() / lib.direct_sum('ia,jb->ijab', eia[:,p0:p1], eia)) emp2 += 2 * numpy.einsum('ijab,iajb', t2[:,:,p0:p1], eris_ovov) emp2 -= numpy.einsum('jiab,iajb', t2[:,:,p0:p1], eris_ovov) self.emp2 = emp2.real logger.info(self, 'Init t2, MP2 energy = %.15g E_corr(MP2) %.15g', e_hf + self.emp2, self.emp2) logger.timer(self, 'init mp2', *time0) return self.emp2, t1, t2
energy = energy _add_vvvv = _add_vvvv update_amps = update_amps
[docs] def kernel(self, t1=None, t2=None, eris=None): return self.ccsd(t1, t2, eris)
[docs] def ccsd(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, callback=self.callback) self._finalize() return self.e_corr, self.t1, self.t2
def _finalize(self): '''Hook for dumping results and clearing up the object.''' if self.converged: logger.info(self, '%s converged', self.__class__.__name__) else: logger.note(self, '%s not converged', self.__class__.__name__) logger.note(self, 'E(%s) = %.16g E_corr = %.16g', self.__class__.__name__, self.e_tot, self.e_corr) return self as_scanner = as_scanner restore_from_diis_ = restore_from_diis_
[docs] def solve_lambda(self, t1=None, t2=None, l1=None, l2=None, eris=None): raise NotImplementedError
[docs] def ccsd_t(self, t1=None, t2=None, eris=None): raise NotImplementedError
[docs] def ipccsd(self, nroots=1, left=False, koopmans=False, guess=None, partition=None, eris=None): raise NotImplementedError
[docs] def eaccsd(self, nroots=1, left=False, koopmans=False, guess=None, partition=None, eris=None): raise NotImplementedError
[docs] def eeccsd(self, nroots=1, koopmans=False, guess=None, eris=None): raise NotImplementedError
[docs] def eomee_ccsd_singlet(self, nroots=1, koopmans=False, guess=None, eris=None): raise NotImplementedError
[docs] def eomee_ccsd_triplet(self, nroots=1, koopmans=False, guess=None, eris=None): raise NotImplementedError
[docs] def eomsf_ccsd(self, nroots=1, koopmans=False, guess=None, eris=None): raise NotImplementedError
[docs] def eomip_method(self): raise NotImplementedError
[docs] def eomea_method(self): raise NotImplementedError
[docs] def eomee_method(self): raise NotImplementedError
[docs] 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''' raise NotImplementedError
[docs] 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] = <p^+ q^+ s r> ''' raise NotImplementedError
[docs] def ao2mo(self, mo_coeff=None): # Pseudo code how eris are implemented: # nocc = self.nocc # nmo = self.nmo # nvir = nmo - nocc # eris = _ChemistsERIs() # eri = ao2mo.incore.full(self._scf._eri, mo_coeff) # eri = ao2mo.restore(1, eri, nmo) # eris.oooo = eri[:nocc,:nocc,:nocc,:nocc].copy() # eris.ovoo = eri[:nocc,nocc:,:nocc,:nocc].copy() # eris.ovvo = eri[nocc:,:nocc,nocc:,:nocc].copy() # eris.ovov = eri[nocc:,:nocc,:nocc,nocc:].copy() # eris.oovv = eri[:nocc,:nocc,nocc:,nocc:].copy() # ovvv = eri[:nocc,nocc:,nocc:,nocc:].copy() # eris.ovvv = lib.pack_tril(ovvv.reshape(-1,nvir,nvir)) # eris.vvvv = ao2mo.restore(4, eri[nocc:,nocc:,nocc:,nocc:], nvir) # eris.fock = numpy.diag(self._scf.mo_energy) # return eris nmo = self.nmo nao = self.mo_coeff.shape[0] nmo_pair = nmo * (nmo+1) // 2 nao_pair = nao * (nao+1) // 2 mem_incore = (max(nao_pair**2, nmo**4) + nmo_pair**2) * 8/1e6 mem_now = lib.current_memory()[0] if (self._scf._eri is not None and (mem_incore+mem_now < self.max_memory or self.incore_complete)): return _make_eris_incore(self, mo_coeff) elif getattr(self._scf, 'with_df', None): logger.warn(self, 'CCSD detected DF being used in the HF object. ' 'MO integrals are computed based on the DF 3-index tensors.\n' 'It\'s recommended to use dfccsd.CCSD for the ' 'DF-CCSD calculations') return _make_df_eris_outcore(self, mo_coeff) else: return _make_eris_outcore(self, mo_coeff)
[docs] def run_diis(self, t1, t2, istep, normt, de, adiis): if (adiis and istep >= self.diis_start_cycle and abs(de) < self.diis_start_energy_diff): vec = self.amplitudes_to_vector(t1, t2) t1, t2 = self.vector_to_amplitudes(adiis.update(vec)) logger.debug1(self, 'DIIS for step %d', istep) return t1, t2
[docs] def amplitudes_to_vector(self, t1, t2, out=None): return amplitudes_to_vector(t1, t2, out)
[docs] def vector_to_amplitudes(self, vec, nmo=None, nocc=None): if nocc is None: nocc = self.nocc if nmo is None: nmo = self.nmo return vector_to_amplitudes(vec, nmo, nocc)
[docs] def vector_size(self, nmo=None, nocc=None): if nocc is None: nocc = self.nocc if nmo is None: nmo = self.nmo nvir = nmo - nocc nov = nocc * nvir return nov + nov*(nov+1)//2
[docs] def dump_chk(self, t1_t2=None, frozen=None, mo_coeff=None, mo_occ=None): if not self.chkfile: return self if t1_t2 is None: t1_t2 = self.t1, self.t2 t1, t2 = t1_t2 if frozen is None: frozen = self.frozen # "None" cannot be serialized by the chkfile module if frozen is None: frozen = 0 cc_chk = {'e_corr': self.e_corr, 't1': t1, 't2': t2, 'frozen': frozen} if mo_coeff is not None: cc_chk['mo_coeff'] = mo_coeff if mo_occ is not None: cc_chk['mo_occ'] = mo_occ if self._nmo is not None: cc_chk['_nmo'] = self._nmo if self._nocc is not None: cc_chk['_nocc'] = self._nocc lib.chkfile.save(self.chkfile, 'ccsd', cc_chk)
[docs] def density_fit(self, auxbasis=None, with_df=None): raise NotImplementedError
[docs] def nuc_grad_method(self): raise NotImplementedError
[docs] class CCSD(CCSDBase): __doc__ = CCSDBase.__doc__
[docs] 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
[docs] 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
[docs] 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)
[docs] 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)
[docs] 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)
[docs] 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)
[docs] 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)
[docs] 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)
[docs] 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)
[docs] def eomip_method(self): from pyscf.cc import eom_rccsd return eom_rccsd.EOMIP(self)
[docs] def eomea_method(self): from pyscf.cc import eom_rccsd return eom_rccsd.EOMEA(self)
[docs] def eomee_method(self): from pyscf.cc import eom_rccsd return eom_rccsd.EOMEE(self)
[docs] 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)
[docs] 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] = <p^+ q^+ s r> ''' 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)
[docs] 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) if with_df is not None: mycc.with_df = with_df if mycc.with_df.auxbasis != auxbasis: mycc.with_df = mycc.with_df.copy() mycc.with_df.auxbasis = auxbasis return mycc
[docs] def nuc_grad_method(self): from pyscf.grad import ccsd return ccsd.Gradients(self)
[docs] def get_t1_diagnostic(self, t1=None): if t1 is None: t1 = self.t1 return get_t1_diagnostic(t1)
[docs] def get_d1_diagnostic(self, t1=None): if t1 is None: t1 = self.t1 return get_d1_diagnostic(t1)
[docs] def get_d2_diagnostic(self, t2=None): if t2 is None: t2 = self.t2 return get_d2_diagnostic(t2)
CC = RCCSD = CCSD class _ChemistsERIs: '''(pq|rs)''' def __init__(self, mol=None): self.mol = mol self.mo_coeff = None self.nocc = None self.fock = None self.oooo = None self.ovoo = None self.oovv = None self.ovvo = None self.ovov = None self.ovvv = None self.vvvv = None def _common_init_(self, mycc, mo_coeff=None): if mo_coeff is None: mo_coeff = mycc.mo_coeff self.mo_coeff = mo_coeff = _mo_without_core(mycc, mo_coeff) # Note: Recomputed fock matrix and HF energy since SCF may not be fully converged. dm = mycc._scf.make_rdm1(mycc.mo_coeff, mycc.mo_occ) vhf = mycc._scf.get_veff(mycc.mol, dm) fockao = mycc._scf.get_fock(vhf=vhf, dm=dm) self.fock = reduce(numpy.dot, (mo_coeff.conj().T, fockao, mo_coeff)) nocc = self.nocc = mycc.nocc self.mol = mycc.mol # Note self.mo_energy can be different to fock.diagonal(). # self.mo_energy is used in the initial guess function (to generate # MP2 amplitudes) and CCSD update_amps preconditioner. # fock.diagonal() should only be used to compute the expectation value # of Slater determinants. mo_e = self.mo_energy = self.fock.diagonal().real try: gap = abs(mo_e[:nocc,None] - mo_e[None,nocc:]).min() if gap < 1e-5: logger.warn(mycc, 'HOMO-LUMO gap %s too small for CCSD.\n' 'CCSD may be difficult to converge. Increasing ' 'CCSD Attribute level_shift may improve ' 'convergence.', gap) except ValueError: # gap.size == 0 pass return self def get_ovvv(self, *slices): '''To access a subblock of ovvv tensor''' ovw = numpy.asarray(self.ovvv[slices]) nocc, nvir, nvir_pair = ovw.shape ovvv = lib.unpack_tril(ovw.reshape(nocc*nvir,nvir_pair)) nvir1 = ovvv.shape[2] return ovvv.reshape(nocc,nvir,nvir1,nvir1) def _contract_vvvv_t2(self, mycc, t2, vvvv_or_direct=False, out=None, verbose=None): if isinstance(vvvv_or_direct, numpy.ndarray): vvvv = vvvv_or_direct elif vvvv_or_direct: # AO-direct contraction vvvv = None else: vvvv = self.vvvv return _contract_vvvv_t2(mycc, self.mol, vvvv, t2, out, verbose) def _contract_vvvv_oov(self, mycc, r2, out=None): raise NotImplementedError def _contract_vvvv_ovv(self, mycc, r2, out=None): raise NotImplementedError def _make_eris_incore(mycc, mo_coeff=None): cput0 = (logger.process_clock(), logger.perf_counter()) eris = _ChemistsERIs() eris._common_init_(mycc, mo_coeff) nocc = eris.nocc nmo = eris.fock.shape[0] nvir = nmo - nocc eri1 = ao2mo.incore.full(mycc._scf._eri, eris.mo_coeff) #:eri1 = ao2mo.restore(1, eri1, nmo) #:eris.oooo = eri1[:nocc,:nocc,:nocc,:nocc].copy() #:eris.ovoo = eri1[:nocc,nocc:,:nocc,:nocc].copy() #:eris.ovvo = eri1[:nocc,nocc:,nocc:,:nocc].copy() #:eris.ovov = eri1[:nocc,nocc:,:nocc,nocc:].copy() #:eris.oovv = eri1[:nocc,:nocc,nocc:,nocc:].copy() #:ovvv = eri1[:nocc,nocc:,nocc:,nocc:].copy() #:eris.ovvv = lib.pack_tril(ovvv.reshape(-1,nvir,nvir)).reshape(nocc,nvir,-1) #:eris.vvvv = ao2mo.restore(4, eri1[nocc:,nocc:,nocc:,nocc:], nvir) if eri1.ndim == 4: eri1 = ao2mo.restore(4, eri1, nmo) nvir_pair = nvir * (nvir+1) // 2 eris.oooo = numpy.empty((nocc,nocc,nocc,nocc)) eris.ovoo = numpy.empty((nocc,nvir,nocc,nocc)) eris.ovvo = numpy.empty((nocc,nvir,nvir,nocc)) eris.ovov = numpy.empty((nocc,nvir,nocc,nvir)) eris.ovvv = numpy.empty((nocc,nvir,nvir_pair)) eris.vvvv = numpy.empty((nvir_pair,nvir_pair)) ij = 0 outbuf = numpy.empty((nmo,nmo,nmo)) oovv = numpy.empty((nocc,nocc,nvir,nvir)) for i in range(nocc): buf = lib.unpack_tril(eri1[ij:ij+i+1], out=outbuf[:i+1]) for j in range(i+1): eris.oooo[i,j] = eris.oooo[j,i] = buf[j,:nocc,:nocc] oovv[i,j] = oovv[j,i] = buf[j,nocc:,nocc:] ij += i + 1 eris.oovv = oovv oovv = None ij1 = 0 for i in range(nocc,nmo): buf = lib.unpack_tril(eri1[ij:ij+i+1], out=outbuf[:i+1]) eris.ovoo[:,i-nocc] = buf[:nocc,:nocc,:nocc] eris.ovvo[:,i-nocc] = buf[:nocc,nocc:,:nocc] eris.ovov[:,i-nocc] = buf[:nocc,:nocc,nocc:] eris.ovvv[:,i-nocc] = lib.pack_tril(buf[:nocc,nocc:,nocc:]) dij = i - nocc + 1 lib.pack_tril(buf[nocc:i+1,nocc:,nocc:], out=eris.vvvv[ij1:ij1+dij]) ij += i + 1 ij1 += dij logger.timer(mycc, 'CCSD integral transformation', *cput0) return eris def _make_eris_outcore(mycc, mo_coeff=None): cput0 = (logger.process_clock(), logger.perf_counter()) log = logger.Logger(mycc.stdout, mycc.verbose) eris = _ChemistsERIs() eris._common_init_(mycc, mo_coeff) mol = mycc.mol mo_coeff = numpy.asarray(eris.mo_coeff, order='F') nocc = eris.nocc nao, nmo = mo_coeff.shape nvir = nmo - nocc orbo = mo_coeff[:,:nocc] orbv = mo_coeff[:,nocc:] nvpair = nvir * (nvir+1) // 2 eris.feri1 = lib.H5TmpFile() eris.oooo = eris.feri1.create_dataset('oooo', (nocc,nocc,nocc,nocc), 'f8') eris.oovv = eris.feri1.create_dataset('oovv', (nocc,nocc,nvir,nvir), 'f8', chunks=(nocc,nocc,1,nvir)) eris.ovoo = eris.feri1.create_dataset('ovoo', (nocc,nvir,nocc,nocc), 'f8', chunks=(nocc,1,nocc,nocc)) eris.ovvo = eris.feri1.create_dataset('ovvo', (nocc,nvir,nvir,nocc), 'f8', chunks=(nocc,1,nvir,nocc)) eris.ovov = eris.feri1.create_dataset('ovov', (nocc,nvir,nocc,nvir), 'f8', chunks=(nocc,1,nocc,nvir)) eris.ovvv = eris.feri1.create_dataset('ovvv', (nocc,nvir,nvpair), 'f8') def save_occ_frac(p0, p1, eri): eri = eri.reshape(p1-p0,nocc,nmo,nmo) eris.oooo[p0:p1] = eri[:,:,:nocc,:nocc] eris.oovv[p0:p1] = eri[:,:,nocc:,nocc:] def save_vir_frac(p0, p1, eri): eri = eri.reshape(p1-p0,nocc,nmo,nmo) eris.ovoo[:,p0:p1] = eri[:,:,:nocc,:nocc].transpose(1,0,2,3) eris.ovvo[:,p0:p1] = eri[:,:,nocc:,:nocc].transpose(1,0,2,3) eris.ovov[:,p0:p1] = eri[:,:,:nocc,nocc:].transpose(1,0,2,3) vvv = lib.pack_tril(eri[:,:,nocc:,nocc:].reshape((p1-p0)*nocc,nvir,nvir)) eris.ovvv[:,p0:p1] = vvv.reshape(p1-p0,nocc,nvpair).transpose(1,0,2) cput1 = logger.process_clock(), logger.perf_counter() if not mycc.direct: max_memory = max(MEMORYMIN, mycc.max_memory-lib.current_memory()[0]) eris.feri2 = lib.H5TmpFile() ao2mo.full(mol, orbv, eris.feri2, max_memory=max_memory, verbose=log) eris.vvvv = eris.feri2['eri_mo'] cput1 = log.timer_debug1('transforming vvvv', *cput1) fswap = lib.H5TmpFile() max_memory = max(MEMORYMIN, mycc.max_memory-lib.current_memory()[0]) int2e = mol._add_suffix('int2e') ao2mo.outcore.half_e1(mol, (mo_coeff,orbo), fswap, int2e, 's4', 1, max_memory, verbose=log) ao_loc = mol.ao_loc_nr() nao_pair = nao * (nao+1) // 2 blksize = int(min(8e9,max_memory*.5e6)/8/(nao_pair+nmo**2)/nocc) blksize = min(nmo, max(BLKMIN, blksize)) log.debug1('blksize %d', blksize) cput2 = cput1 fload = ao2mo.outcore._load_from_h5g buf = numpy.empty((blksize*nocc,nao_pair)) buf_prefetch = numpy.empty_like(buf) def load(buf_prefetch, p0, rowmax): if p0 < rowmax: p1 = min(rowmax, p0+blksize) fload(fswap['0'], p0*nocc, p1*nocc, buf_prefetch) outbuf = numpy.empty((blksize*nocc,nmo**2)) with lib.call_in_background(load, sync=not mycc.async_io) as prefetch: prefetch(buf_prefetch, 0, nocc) for p0, p1 in lib.prange(0, nocc, blksize): buf, buf_prefetch = buf_prefetch, buf prefetch(buf_prefetch, p1, nocc) nrow = (p1 - p0) * nocc dat = ao2mo._ao2mo.nr_e2(buf[:nrow], mo_coeff, (0,nmo,0,nmo), 's4', 's1', out=outbuf, ao_loc=ao_loc) save_occ_frac(p0, p1, dat) cput2 = log.timer_debug1('transforming oopp', *cput2) prefetch(buf_prefetch, nocc, nmo) for p0, p1 in lib.prange(0, nvir, blksize): buf, buf_prefetch = buf_prefetch, buf prefetch(buf_prefetch, nocc+p1, nmo) nrow = (p1 - p0) * nocc dat = ao2mo._ao2mo.nr_e2(buf[:nrow], mo_coeff, (0,nmo,0,nmo), 's4', 's1', out=outbuf, ao_loc=ao_loc) save_vir_frac(p0, p1, dat) cput2 = log.timer_debug1('transforming ovpp [%d:%d]'%(p0,p1), *cput2) cput1 = log.timer_debug1('transforming oppp', *cput1) log.timer('CCSD integral transformation', *cput0) return eris def _make_df_eris_outcore(mycc, mo_coeff=None): cput0 = (logger.process_clock(), logger.perf_counter()) log = logger.Logger(mycc.stdout, mycc.verbose) eris = _ChemistsERIs() eris._common_init_(mycc, mo_coeff) mo_coeff = numpy.asarray(eris.mo_coeff, order='F') nocc = eris.nocc nao, nmo = mo_coeff.shape nvir = nmo - nocc nvir_pair = nvir*(nvir+1)//2 naux = mycc._scf.with_df.get_naoaux() Loo = numpy.empty((naux,nocc,nocc)) Lov = numpy.empty((naux,nocc,nvir)) Lvo = numpy.empty((naux,nvir,nocc)) Lvv = numpy.empty((naux,nvir_pair)) ijslice = (0, nmo, 0, nmo) Lpq = None p1 = 0 for eri1 in mycc._scf.with_df.loop(): Lpq = _ao2mo.nr_e2(eri1, mo_coeff, ijslice, aosym='s2', out=Lpq).reshape(-1,nmo,nmo) p0, p1 = p1, p1 + Lpq.shape[0] Loo[p0:p1] = Lpq[:,:nocc,:nocc] Lov[p0:p1] = Lpq[:,:nocc,nocc:] Lvo[p0:p1] = Lpq[:,nocc:,:nocc] Lvv[p0:p1] = lib.pack_tril(Lpq[:,nocc:,nocc:].reshape(-1,nvir,nvir)) Loo = Loo.reshape(naux,nocc*nocc) Lov = Lov.reshape(naux,nocc*nvir) Lvo = Lvo.reshape(naux,nocc*nvir) eris.feri1 = lib.H5TmpFile() eris.oooo = eris.feri1.create_dataset('oooo', (nocc,nocc,nocc,nocc), 'f8') eris.oovv = eris.feri1.create_dataset('oovv', (nocc,nocc,nvir,nvir), 'f8', chunks=(nocc,nocc,1,nvir)) eris.ovoo = eris.feri1.create_dataset('ovoo', (nocc,nvir,nocc,nocc), 'f8', chunks=(nocc,1,nocc,nocc)) eris.ovvo = eris.feri1.create_dataset('ovvo', (nocc,nvir,nvir,nocc), 'f8', chunks=(nocc,1,nvir,nocc)) eris.ovov = eris.feri1.create_dataset('ovov', (nocc,nvir,nocc,nvir), 'f8', chunks=(nocc,1,nocc,nvir)) eris.ovvv = eris.feri1.create_dataset('ovvv', (nocc,nvir,nvir_pair), 'f8') eris.vvvv = eris.feri1.create_dataset('vvvv', (nvir_pair,nvir_pair), 'f8') eris.oooo[:] = lib.ddot(Loo.T, Loo).reshape(nocc,nocc,nocc,nocc) eris.ovoo[:] = lib.ddot(Lov.T, Loo).reshape(nocc,nvir,nocc,nocc) eris.oovv[:] = lib.unpack_tril(lib.ddot(Loo.T, Lvv)).reshape(nocc,nocc,nvir,nvir) eris.ovvo[:] = lib.ddot(Lov.T, Lvo).reshape(nocc,nvir,nvir,nocc) eris.ovov[:] = lib.ddot(Lov.T, Lov).reshape(nocc,nvir,nocc,nvir) eris.ovvv[:] = lib.ddot(Lov.T, Lvv).reshape(nocc,nvir,nvir_pair) eris.vvvv[:] = lib.ddot(Lvv.T, Lvv) log.timer('CCSD integral transformation', *cput0) return eris def _flops(nocc, nvir): '''Total float points''' return (nocc**3*nvir**2*2 + nocc**2*nvir**3*2 + # Ftilde nocc**4*nvir*2 * 2 + nocc**4*nvir**2*2 + # Wijkl nocc*nvir**4*2 * 2 + # Wabcd nocc**2*nvir**3*2 + nocc**3*nvir**2*2 + nocc**3*nvir**3*2 + nocc**3*nvir**3*2 + nocc**2*nvir**3*2 + nocc**3*nvir**2*2 + # Wiabj nocc**2*nvir**3*2 + nocc**3*nvir**2*2 + # t1 nocc**3*nvir**2*2 * 2 + nocc**4*nvir**2*2 + nocc*(nocc+1)/2*nvir**4*2 + # vvvv nocc**2*nvir**3*2 * 2 + nocc**3*nvir**2*2 * 2 + # t2 nocc**3*nvir**3*2 + nocc**3*nvir**3*2 * 2 + nocc**3*nvir**2*2 * 4) # Wiabj 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.basis = {'H': 'cc-pvdz', 'O': 'cc-pvdz',} mol.build() rhf = scf.RHF(mol) rhf.scf() # -76.0267656731 mf = rhf.density_fit(auxbasis='weigend') mf._eri = None mcc = CCSD(mf) eris = mcc.ao2mo() emp2, t1, t2 = mcc.init_amps(eris) print(abs(t2).sum() - 4.9318753386922278) print(emp2 - -0.20401737899811551) t1, t2 = update_amps(mcc, t1, t2, eris) print(abs(t1).sum() - 0.046961325647584914) print(abs(t2).sum() - 5.378260578551683 ) mcc = CCSD(rhf) eris = mcc.ao2mo() emp2, t1, t2 = mcc.init_amps(eris) print(abs(t2).sum() - 4.9556571218177) print(emp2 - -0.2040199672883385) t1, t2 = update_amps(mcc, t1, t2, eris) print(abs(t1).sum()-0.0475038989126) print(abs(t2).sum()-5.401823846018721) print(energy(mcc, t1, t2, eris) - -0.208967840546667) t1, t2 = update_amps(mcc, t1, t2, eris) print(energy(mcc, t1, t2, eris) - -0.212173678670510) print(abs(t1).sum() - 0.05470123093500083) print(abs(t2).sum() - 5.5605208391876539) mcc.ccsd() print(mcc.ecc - -0.213343234198275) print(abs(mcc.t2).sum() - 5.63970304662375) mcc.max_memory = 1 mcc.direct = True mcc.ccsd() print(mcc.ecc - -0.213343234198275) print(abs(mcc.t2).sum() - 5.63970304662375) e, v = mcc.ipccsd(nroots=3) print(e[0] - 0.43356041409195489) print(e[1] - 0.51876598058509493) print(e[2] - 0.6782879569941862 ) e, v = mcc.eeccsd(nroots=4) print(e[0] - 0.2757159395886167) print(e[1] - 0.2757159395886167) print(e[2] - 0.2757159395886167) print(e[3] - 0.3005716731825082)