Source code for pyscf.pbc.cc.kccsd_t_rhf_slow

#!/usr/bin/env python
# Copyright 2017-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.
#
# Authors: James D. McClain <jmcclain@princeton.edu>
#
"""Module for running restricted closed-shell k-point ccsd(t)"""

import h5py
import itertools
import numpy as np
import pyscf.pbc.cc.kccsd_rhf


from itertools import product
from pyscf import lib
#from pyscf import _ccsd
from pyscf.lib import logger
from pyscf.lib.misc import tril_product
from pyscf.lib.misc import flatten
from pyscf.lib.numpy_helper import cartesian_prod
from pyscf.lib.numpy_helper import pack_tril
from pyscf.lib.parameters import LARGE_DENOM
from pyscf.pbc.lib import kpts_helper
from pyscf.pbc.mp.kmp2 import (get_frozen_mask, get_nocc, get_nmo,
                               padded_mo_coeff, padding_k_idx)
from pyscf import __config__

#einsum = np.einsum
einsum = lib.einsum

# CCSD(T) equations taken from Scuseria, JCP (94), 1991
#
# NOTE: As pointed out in cc/ccsd_t_slow.py, there is an error in this paper
#     and the equation should read [ia] >= [jb] >= [kc] (since the only
#     symmetry in spin-less operators is the exchange of a column of excitation
#     ooperators).
[docs] def kernel(mycc, eris, t1=None, t2=None, max_memory=2000, verbose=logger.INFO): '''Returns the CCSD(T) for restricted closed-shell systems with k-points. Note: Returns real part of the CCSD(T) energy, raises warning if there is a complex part. Args: mycc (:class:`RCCSD`): Coupled-cluster object storing results of a coupled-cluster calculation. eris (:class:`_ERIS`): Integral object holding the relevant electron- repulsion integrals and Fock matrix elements t1 (:obj:`ndarray`): t1 coupled-cluster amplitudes t2 (:obj:`ndarray`): t2 coupled-cluster amplitudes max_memory (float): Maximum memory used in calculation (NOT USED) verbose (int, :class:`Logger`): verbosity of calculation Returns: energy_t (float): The real-part of the k-point CCSD(T) energy. ''' assert isinstance(mycc, pyscf.pbc.cc.kccsd_rhf.RCCSD) cpu1 = cpu0 = (logger.process_clock(), logger.perf_counter()) if isinstance(verbose, logger.Logger): log = verbose else: log = logger.Logger(mycc.stdout, verbose) if t1 is None: t1 = mycc.t1 if t2 is None: t2 = mycc.t2 if eris is None: raise TypeError('Electron repulsion integrals, `eris`, must be passed in ' 'to the CCSD(T) kernel or created in the cc object for ' 'the k-point CCSD(T) to run!') if t1 is None or t2 is None: raise TypeError('Must pass in t1/t2 amplitudes to k-point CCSD(T)! (Maybe ' 'need to run `.ccsd()` on the ccsd object?)') cell = mycc._scf.cell kpts = mycc.kpts # The dtype of any local arrays that will be created dtype = t1.dtype nkpts, nocc, nvir = t1.shape mo_energy_occ = [eris.mo_energy[ki][:nocc] for ki in range(nkpts)] mo_energy_vir = [eris.mo_energy[ki][nocc:] for ki in range(nkpts)] fov = eris.fock[:, :nocc, nocc:] # Set up class for k-point conservation kconserv = kpts_helper.get_kconserv(cell, kpts) cpu1 = log.timer_debug1('CCSD(T) tmp eri creation', *cpu1) def get_w(ki, kj, kk, ka, kb, kc, a0, a1, b0, b1, c0, c1, out=None): '''Wijkabc intermediate as described in Scuseria paper before Pijkabc acts''' km = kconserv[ki, ka, kj] kf = kconserv[kk, kc, kj] ret = einsum('kjcf,fiba->abcijk', t2[kk,kj,kc,:,:,c0:c1,:], eris.vovv[kf,ki,kb,:,:,b0:b1,a0:a1].conj()) ret = ret - einsum('mkbc,jima->abcijk', t2[km,kk,kb,:,:,b0:b1,c0:c1], eris.ooov[kj,ki,km,:,:,:,a0:a1].conj()) return ret def get_permuted_w(ki, kj, kk, ka, kb, kc, orb_indices): '''Pijkabc operating on Wijkabc intermediate as described in Scuseria paper''' a0, a1, b0, b1, c0, c1 = orb_indices out = get_w(ki, kj, kk, ka, kb, kc, a0, a1, b0, b1, c0, c1) out = out + get_w(kj, kk, ki, kb, kc, ka, b0, b1, c0, c1, a0, a1).transpose(2,0,1,5,3,4) out = out + get_w(kk, ki, kj, kc, ka, kb, c0, c1, a0, a1, b0, b1).transpose(1,2,0,4,5,3) out = out + get_w(ki, kk, kj, ka, kc, kb, a0, a1, c0, c1, b0, b1).transpose(0,2,1,3,5,4) out = out + get_w(kk, kj, ki, kc, kb, ka, c0, c1, b0, b1, a0, a1).transpose(2,1,0,5,4,3) out = out + get_w(kj, ki, kk, kb, ka, kc, b0, b1, a0, a1, c0, c1).transpose(1,0,2,4,3,5) return out def get_rw(ki, kj, kk, ka, kb, kc, orb_indices): '''R operating on Wijkabc intermediate as described in Scuseria paper''' a0, a1, b0, b1, c0, c1 = orb_indices ret = (4. * get_permuted_w(ki,kj,kk,ka,kb,kc,orb_indices) + 1. * get_permuted_w(kj,kk,ki,ka,kb,kc,orb_indices).transpose(0,1,2,5,3,4) + 1. * get_permuted_w(kk,ki,kj,ka,kb,kc,orb_indices).transpose(0,1,2,4,5,3) - 2. * get_permuted_w(ki,kk,kj,ka,kb,kc,orb_indices).transpose(0,1,2,3,5,4) - 2. * get_permuted_w(kk,kj,ki,ka,kb,kc,orb_indices).transpose(0,1,2,5,4,3) - 2. * get_permuted_w(kj,ki,kk,ka,kb,kc,orb_indices).transpose(0,1,2,4,3,5)) return ret def get_v(ki, kj, kk, ka, kb, kc, a0, a1, b0, b1, c0, c1): '''Vijkabc intermediate as described in Scuseria paper''' out = np.zeros((a1-a0,b1-b0,c1-c0) + (nocc,)*3, dtype=dtype) if kk == kc: out = out + einsum('kc,ijab->abcijk', t1[kk,:,c0:c1], eris.oovv[ki,kj,ka,:,:,a0:a1,b0:b1].conj()) out = out + einsum('kc,ijab->abcijk', fov[kk,:,c0:c1], t2[ki,kj,ka,:,:,a0:a1,b0:b1]) return out def get_permuted_v(ki, kj, kk, ka, kb, kc, orb_indices): '''Pijkabc operating on Vijkabc intermediate as described in Scuseria paper''' a0, a1, b0, b1, c0, c1 = orb_indices ret = get_v(ki, kj, kk, ka, kb, kc, a0, a1, b0, b1, c0, c1) ret = ret + get_v(kj, kk, ki, kb, kc, ka, b0, b1, c0, c1, a0, a1).transpose(2,0,1,5,3,4) ret = ret + get_v(kk, ki, kj, kc, ka, kb, c0, c1, a0, a1, b0, b1).transpose(1,2,0,4,5,3) ret = ret + get_v(ki, kk, kj, ka, kc, kb, a0, a1, c0, c1, b0, b1).transpose(0,2,1,3,5,4) ret = ret + get_v(kk, kj, ki, kc, kb, ka, c0, c1, b0, b1, a0, a1).transpose(2,1,0,5,4,3) ret = ret + get_v(kj, ki, kk, kb, ka, kc, b0, b1, a0, a1, c0, c1).transpose(1,0,2,4,3,5) return ret energy_t = 0.0 # Get location of padded elements in occupied and virtual space nonzero_opadding, nonzero_vpadding = padding_k_idx(mycc, kind="split") mem_now = lib.current_memory()[0] max_memory = max(0, mycc.max_memory - mem_now) blkmin = 4 # temporary t3 array is size: blksize**3 * nocc**3 * 16 vir_blksize = min(nvir, max(blkmin, int((max_memory*.9e6/16/nocc**3)**(1./3)))) tasks = [] log.debug('max_memory %d MB (%d MB in use)', max_memory, mem_now) log.debug('virtual blksize = %d (nvir = %d)', nvir, vir_blksize) for a0, a1 in lib.prange(0, nvir, vir_blksize): for b0, b1 in lib.prange(0, nvir, vir_blksize): for c0, c1 in lib.prange(0, nvir, vir_blksize): tasks.append((a0,a1,b0,b1,c0,c1)) for ka in range(nkpts): for kb in range(ka+1): for ki, kj, kk in product(range(nkpts), repeat=3): # eigenvalue denominator: e(i) + e(j) + e(k) eijk = LARGE_DENOM * np.ones((nocc,)*3, dtype=mo_energy_occ[0].dtype) n0_ovp_ijk = np.ix_(nonzero_opadding[ki], nonzero_opadding[kj], nonzero_opadding[kk]) eijk[n0_ovp_ijk] = lib.direct_sum('i,j,k->ijk', mo_energy_occ[ki], mo_energy_occ[kj], mo_energy_occ[kk])[n0_ovp_ijk] # Find momentum conservation condition for triples # amplitude t3ijkabc kc = kpts_helper.get_kconserv3(cell, kpts, [ki, kj, kk, ka, kb]) if not (ka >= kb and kb >= kc): continue if ka == kb and kb == kc: symm_kpt = 1. elif ka == kb or kb == kc: symm_kpt = 3. else: symm_kpt = 6. eabc = LARGE_DENOM * np.ones((nvir,)*3, dtype=mo_energy_occ[0].dtype) n0_ovp_abc = np.ix_(nonzero_vpadding[ka], nonzero_vpadding[kb], nonzero_vpadding[kc]) eabc[n0_ovp_abc] = lib.direct_sum('i,j,k->ijk', mo_energy_vir[ka], mo_energy_vir[kb], mo_energy_vir[kc])[n0_ovp_abc] for task_id, task in enumerate(tasks): eijkabc = (eijk[None,None,None,:,:,:] - eabc[a0:a1,b0:b1,c0:c1,None,None,None]) pwijk = (get_permuted_w(ki,kj,kk,ka,kb,kc,task) + get_permuted_v(ki,kj,kk,ka,kb,kc,task) * 0.5) rwijk = get_rw(ki,kj,kk,ka,kb,kc,task) / eijkabc energy_t += symm_kpt * einsum('abcijk,abcijk', pwijk, rwijk.conj()) energy_t *= (1. / 3) energy_t /= nkpts if abs(energy_t.imag) > 1e-4: log.warn('Non-zero imaginary part of CCSD(T) energy was found %s', energy_t.imag) log.timer('CCSD(T)', *cpu0) log.note('CCSD(T) correction per cell = %.15g', energy_t.real) log.note('CCSD(T) correction per cell (imag) = %.15g', energy_t.imag) return energy_t.real
if __name__ == '__main__': from pyscf.pbc import gto from pyscf.pbc import scf from pyscf.pbc import cc cell = gto.Cell() cell.atom = ''' C 0.000000000000 0.000000000000 0.000000000000 C 1.685068664391 1.685068664391 1.685068664391 ''' cell.basis = 'gth-szv' cell.pseudo = 'gth-pade' cell.a = ''' 0.000000000, 3.370137329, 3.370137329 3.370137329, 0.000000000, 3.370137329 3.370137329, 3.370137329, 0.000000000''' cell.conv_tol = 1e-12 cell.conv_tol_grad = 1e-12 cell.direct_scf_tol = 1e-16 cell.unit = 'B' cell.verbose = 4 cell.mesh = [24, 24, 24] cell.build() nmp = [1,1,4] kpts = cell.make_kpts(nmp) kpts -= kpts[0] kmf = scf.KRHF(cell, kpts=kpts, exxdiv=None) kmf.conv_tol = 1e-12 kmf.conv_tol_grad = 1e-12 kmf.direct_scf_tol = 1e-16 ehf = kmf.kernel() mycc = cc.KRCCSD(kmf) eris = mycc.ao2mo() ecc, t1, t2 = mycc.kernel(eris=eris) energy_t = kernel(mycc, eris=eris, verbose=9) # Start of supercell calculations from pyscf.pbc.tools.pbc import super_cell supcell = super_cell(cell, nmp) supcell.build() kmf = scf.RHF(supcell, exxdiv=None) kmf.conv_tol = 1e-12 kmf.conv_tol_grad = 1e-12 kmf.direct_scf_tol = 1e-16 sup_ehf = kmf.kernel() myscc = cc.RCCSD(kmf) eris = myscc.ao2mo() sup_ecc, t1, t2 = myscc.kernel(eris=eris) sup_energy_t = myscc.ccsd_t(eris=eris) print("Kpoint CCSD: %20.16f" % ecc) print("Supercell CCSD: %20.16f" % (sup_ecc/np.prod(nmp))) print("Kpoint CCSD(T): %20.16f" % energy_t) print("Supercell CCSD(T): %20.16f" % (sup_energy_t/np.prod(nmp)))