1#!/usr/bin/env python 2# Copyright 2017-2021 The PySCF Developers. All Rights Reserved. 3# 4# Licensed under the Apache License, Version 2.0 (the "License"); 5# you may not use this file except in compliance with the License. 6# You may obtain a copy of the License at 7# 8# http://www.apache.org/licenses/LICENSE-2.0 9# 10# Unless required by applicable law or agreed to in writing, software 11# distributed under the License is distributed on an "AS IS" BASIS, 12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13# See the License for the specific language governing permissions and 14# limitations under the License. 15# 16# Authors: James D. McClain <jmcclain@princeton.edu> 17# 18"""Module for running k-point ccsd(t)""" 19 20import numpy as np 21import pyscf.pbc.cc.kccsd 22 23from pyscf import lib 24from pyscf.lib import logger 25from pyscf.lib.misc import tril_product 26from pyscf.pbc import scf 27from pyscf.pbc.lib import kpts_helper 28from pyscf.lib.numpy_helper import cartesian_prod 29from pyscf.lib.parameters import LOOSE_ZERO_TOL, LARGE_DENOM # noqa 30from pyscf.pbc.cc.kccsd_t_rhf import _get_epqr 31from pyscf.pbc.mp.kmp2 import (get_frozen_mask, get_nocc, get_nmo, 32 padded_mo_coeff, padding_k_idx) # noqa 33 34#einsum = np.einsum 35einsum = lib.einsum 36 37 38# CCSD(T) equations taken from Tu, Yang, Wang, and Guo JPC (135), 2011 39# 40# There are some complex conjugates not included in the equations 41# by Watts, Gauss, Bartlett JCP (98), 1993 42def kernel(mycc, eris, t1=None, t2=None, max_memory=2000, verbose=logger.INFO): 43 '''Returns the CCSD(T) for general spin-orbital integrals with k-points. 44 45 Note: 46 Returns real part of the CCSD(T) energy, raises warning if there is 47 a complex part. 48 49 Args: 50 mycc (:class:`GCCSD`): Coupled-cluster object storing results of 51 a coupled-cluster calculation. 52 eris (:class:`_ERIS`): Integral object holding the relevant electron- 53 repulsion integrals and Fock matrix elements 54 t1 (:obj:`ndarray`): t1 coupled-cluster amplitudes 55 t2 (:obj:`ndarray`): t2 coupled-cluster amplitudes 56 max_memory (float): Maximum memory used in calculation (NOT USED) 57 verbose (int, :class:`Logger`): verbosity of calculation 58 59 Returns: 60 energy_t (float): The real-part of the k-point CCSD(T) energy. 61 ''' 62 assert isinstance(mycc, pyscf.pbc.cc.kccsd.GCCSD) 63 if isinstance(verbose, logger.Logger): 64 log = verbose 65 else: 66 log = logger.Logger(mycc.stdout, verbose) 67 68 if t1 is None: t1 = mycc.t1 69 if t2 is None: t2 = mycc.t2 70 71 if t1 is None or t2 is None: 72 raise TypeError('Must pass in t1/t2 amplitudes to k-point CCSD(T)! (Maybe ' 73 'need to run `.ccsd()` on the ccsd object?)') 74 75 cell = mycc._scf.cell 76 kpts = mycc.kpts 77 78 # The dtype of any local arrays that will be created 79 dtype = t1.dtype 80 81 nkpts, nocc, nvir = t1.shape 82 83 mo_energy_occ = [eris.mo_energy[ki][:nocc] for ki in range(nkpts)] 84 mo_energy_vir = [eris.mo_energy[ki][nocc:] for ki in range(nkpts)] 85 fov = eris.fock[:, :nocc, nocc:] 86 87 mo_e_o = mo_energy_occ 88 mo_e_v = mo_energy_vir 89 90 # Set up class for k-point conservation 91 kconserv = kpts_helper.get_kconserv(cell, kpts) 92 93 # Get location of padded elements in occupied and virtual space 94 nonzero_opadding, nonzero_vpadding = padding_k_idx(mycc, kind="split") 95 96 energy_t = 0.0 97 98 for ki in range(nkpts): 99 for kj in range(ki + 1): 100 for kk in range(kj + 1): 101 # eigenvalue denominator: e(i) + e(j) + e(k) 102 eijk = _get_epqr([0,nocc,ki,mo_e_o,nonzero_opadding], 103 [0,nocc,kj,mo_e_o,nonzero_opadding], 104 [0,nocc,kk,mo_e_o,nonzero_opadding]) 105 106 # Factors to include for permutational symmetry among k-points for occupied space 107 if ki == kj and kj == kk: 108 symm_ijk_kpt = 1. # only one degeneracy 109 elif ki == kj or kj == kk: 110 symm_ijk_kpt = 3. # 3 degeneracies when only one k-point is unique 111 else: 112 symm_ijk_kpt = 6. # 3! combinations of arranging 3 distinct k-points 113 114 for ka in range(nkpts): 115 for kb in range(ka + 1): 116 117 # Find momentum conservation condition for triples 118 # amplitude t3ijkabc 119 kc = kpts_helper.get_kconserv3(cell, kpts, [ki, kj, kk, ka, kb]) 120 if kc not in range(kb + 1): 121 continue 122 123 # -1.0 so the LARGE_DENOM does not cancel with the one from eijk 124 eabc = _get_epqr([0,nvir,ka,mo_e_v,nonzero_vpadding], 125 [0,nvir,kb,mo_e_v,nonzero_vpadding], 126 [0,nvir,kc,mo_e_v,nonzero_vpadding], 127 fac=[-1.,-1.,-1.]) 128 129 # Factors to include for permutational symmetry among k-points for virtual space 130 if ka == kb and kb == kc: 131 symm_abc_kpt = 1. # one unique combination of (ka, kb, kc) 132 elif ka == kb or kb == kc: 133 symm_abc_kpt = 3. # 3 unique combinations of (ka, kb, kc) 134 else: 135 symm_abc_kpt = 6. # 6 unique combinations of (ka, kb, kc) 136 137 # Determine the a, b, c indices we will loop over as 138 # determined by the k-point symmetry. 139 abc_indices = cartesian_prod([range(nvir)] * 3) 140 symm_3d = symm_2d_ab = symm_2d_bc = False 141 if ka == kc: # == kb from lower triangular summation 142 abc_indices = tril_product(range(nvir), repeat=3, tril_idx=[0, 1, 2]) # loop a >= b >= c 143 symm_3d = True 144 elif ka == kb: 145 abc_indices = tril_product(range(nvir), repeat=3, tril_idx=[0, 1]) # loop a >= b 146 symm_2d_ab = True 147 elif kb == kc: 148 abc_indices = tril_product(range(nvir), repeat=3, tril_idx=[1, 2]) # loop b >= c 149 symm_2d_bc = True 150 151 for a, b, c in abc_indices: 152 # See symm_3d and abc_indices above for description of factors 153 symm_abc = 1. 154 if symm_3d: 155 if a == b == c: 156 symm_abc = 1. 157 elif a == b or b == c: 158 symm_abc = 3. 159 else: 160 symm_abc = 6. 161 elif symm_2d_ab: 162 if a == b: 163 symm_abc = 1. 164 else: 165 symm_abc = 2. 166 elif symm_2d_bc: 167 if b == c: 168 symm_abc = 1. 169 else: 170 symm_abc = 2. 171 172 # Form energy denominator 173 eijkabc = (eijk[:,:,:] + eabc[a,b,c]) 174 175 # Form connected triple excitation amplitude 176 t3c = np.zeros((nocc, nocc, nocc), dtype=dtype) 177 178 # First term: 1 - p(ij) - p(ik) 179 ke = kconserv[kj, ka, kk] 180 t3c = t3c + einsum('jke,ie->ijk', t2[kj, kk, ka, :, :, a, :], -eris.ovvv[ki, ke, kc, :, :, c, b].conj()) 181 ke = kconserv[ki, ka, kk] 182 t3c = t3c - einsum('ike,je->ijk', t2[ki, kk, ka, :, :, a, :], -eris.ovvv[kj, ke, kc, :, :, c, b].conj()) 183 ke = kconserv[kj, ka, ki] 184 t3c = t3c - einsum('jie,ke->ijk', t2[kj, ki, ka, :, :, a, :], -eris.ovvv[kk, ke, kc, :, :, c, b].conj()) 185 186 km = kconserv[kb, ki, kc] 187 t3c = t3c - einsum('mi,jkm->ijk', t2[km, ki, kb, :, :, b, c], eris.ooov[kj, kk, km, :, :, :, a].conj()) 188 km = kconserv[kb, kj, kc] 189 t3c = t3c + einsum('mj,ikm->ijk', t2[km, kj, kb, :, :, b, c], eris.ooov[ki, kk, km, :, :, :, a].conj()) 190 km = kconserv[kb, kk, kc] 191 t3c = t3c + einsum('mk,jim->ijk', t2[km, kk, kb, :, :, b, c], eris.ooov[kj, ki, km, :, :, :, a].conj()) 192 193 # Second term: - p(ab) + p(ab) p(ij) + p(ab) p(ik) 194 ke = kconserv[kj, kb, kk] 195 t3c = t3c - einsum('jke,ie->ijk', t2[kj, kk, kb, :, :, b, :], -eris.ovvv[ki, ke, kc, :, :, c, a].conj()) 196 ke = kconserv[ki, kb, kk] 197 t3c = t3c + einsum('ike,je->ijk', t2[ki, kk, kb, :, :, b, :], -eris.ovvv[kj, ke, kc, :, :, c, a].conj()) 198 ke = kconserv[kj, kb, ki] 199 t3c = t3c + einsum('jie,ke->ijk', t2[kj, ki, kb, :, :, b, :], -eris.ovvv[kk, ke, kc, :, :, c, a].conj()) 200 201 km = kconserv[ka, ki, kc] 202 t3c = t3c + einsum('mi,jkm->ijk', t2[km, ki, ka, :, :, a, c], eris.ooov[kj, kk, km, :, :, :, b].conj()) 203 km = kconserv[ka, kj, kc] 204 t3c = t3c - einsum('mj,ikm->ijk', t2[km, kj, ka, :, :, a, c], eris.ooov[ki, kk, km, :, :, :, b].conj()) 205 km = kconserv[ka, kk, kc] 206 t3c = t3c - einsum('mk,jim->ijk', t2[km, kk, ka, :, :, a, c], eris.ooov[kj, ki, km, :, :, :, b].conj()) 207 208 # Third term: - p(ac) + p(ac) p(ij) + p(ac) p(ik) 209 ke = kconserv[kj, kc, kk] 210 t3c = t3c - einsum('jke,ie->ijk', t2[kj, kk, kc, :, :, c, :], -eris.ovvv[ki, ke, ka, :, :, a, b].conj()) 211 ke = kconserv[ki, kc, kk] 212 t3c = t3c + einsum('ike,je->ijk', t2[ki, kk, kc, :, :, c, :], -eris.ovvv[kj, ke, ka, :, :, a, b].conj()) 213 ke = kconserv[kj, kc, ki] 214 t3c = t3c + einsum('jie,ke->ijk', t2[kj, ki, kc, :, :, c, :], -eris.ovvv[kk, ke, ka, :, :, a, b].conj()) 215 216 km = kconserv[kb, ki, ka] 217 t3c = t3c + einsum('mi,jkm->ijk', t2[km, ki, kb, :, :, b, a], eris.ooov[kj, kk, km, :, :, :, c].conj()) 218 km = kconserv[kb, kj, ka] 219 t3c = t3c - einsum('mj,ikm->ijk', t2[km, kj, kb, :, :, b, a], eris.ooov[ki, kk, km, :, :, :, c].conj()) 220 km = kconserv[kb, kk, ka] 221 t3c = t3c - einsum('mk,jim->ijk', t2[km, kk, kb, :, :, b, a], eris.ooov[kj, ki, km, :, :, :, c].conj()) 222 223 # Form disconnected triple excitation amplitude contribution 224 t3d = np.zeros((nocc, nocc, nocc), dtype=dtype) 225 226 # First term: 1 - p(ij) - p(ik) 227 if ki == ka: 228 t3d = t3d + einsum('i,jk->ijk', t1[ki, :, a], -eris.oovv[kj, kk, kb, :, :, b, c].conj()) 229 t3d = t3d + einsum('i,jk->ijk',-fov[ki, :, a], t2[kj, kk, kb, :, :, b, c]) 230 231 if kj == ka: 232 t3d = t3d - einsum('j,ik->ijk', t1[kj, :, a], -eris.oovv[ki, kk, kb, :, :, b, c].conj()) 233 t3d = t3d - einsum('j,ik->ijk',-fov[kj, :, a], t2[ki, kk, kb, :, :, b, c]) 234 235 if kk == ka: 236 t3d = t3d - einsum('k,ji->ijk', t1[kk, :, a], -eris.oovv[kj, ki, kb, :, :, b, c].conj()) 237 t3d = t3d - einsum('k,ji->ijk',-fov[kk, :, a], t2[kj, ki, kb, :, :, b, c]) 238 239 # Second term: - p(ab) + p(ab) p(ij) + p(ab) p(ik) 240 if ki == kb: 241 t3d = t3d - einsum('i,jk->ijk', t1[ki, :, b], -eris.oovv[kj, kk, ka, :, :, a, c].conj()) 242 t3d = t3d - einsum('i,jk->ijk',-fov[ki, :, b], t2[kj, kk, ka, :, :, a, c]) 243 244 if kj == kb: 245 t3d = t3d + einsum('j,ik->ijk', t1[kj, :, b], -eris.oovv[ki, kk, ka, :, :, a, c].conj()) 246 t3d = t3d + einsum('j,ik->ijk',-fov[kj, :, b], t2[ki, kk, ka, :, :, a, c]) 247 248 if kk == kb: 249 t3d = t3d + einsum('k,ji->ijk', t1[kk, :, b], -eris.oovv[kj, ki, ka, :, :, a, c].conj()) 250 t3d = t3d + einsum('k,ji->ijk',-fov[kk, :, b], t2[kj, ki, ka, :, :, a, c]) 251 252 # Third term: - p(ac) + p(ac) p(ij) + p(ac) p(ik) 253 if ki == kc: 254 t3d = t3d - einsum('i,jk->ijk', t1[ki, :, c], -eris.oovv[kj, kk, kb, :, :, b, a].conj()) 255 t3d = t3d - einsum('i,jk->ijk',-fov[ki, :, c], t2[kj, kk, kb, :, :, b, a]) 256 257 if kj == kc: 258 t3d = t3d + einsum('j,ik->ijk', t1[kj, :, c], -eris.oovv[ki, kk, kb, :, :, b, a].conj()) 259 t3d = t3d + einsum('j,ik->ijk',-fov[kj, :, c], t2[ki, kk, kb, :, :, b, a]) 260 261 if kk == kc: 262 t3d = t3d + einsum('k,ji->ijk', t1[kk, :, c], -eris.oovv[kj, ki, kb, :, :, b, a].conj()) 263 t3d = t3d + einsum('k,ji->ijk',-fov[kk, :, c], t2[kj, ki, kb, :, :, b, a]) 264 265 t3c_plus_d = t3c + t3d 266 t3c_plus_d /= eijkabc 267 268 energy_t += symm_abc_kpt * symm_ijk_kpt * symm_abc * einsum('ijk,ijk', t3c, t3c_plus_d.conj()) 269 270 energy_t = (1. / 36) * energy_t / nkpts 271 272 if abs(energy_t.imag) > 1e-4: 273 log.warn('Non-zero imaginary part of CCSD(T) energy was found %s', 274 energy_t.imag) 275 log.note('CCSD(T) correction per cell = %.15g', energy_t.real) 276 return energy_t.real 277 278 279# Gamma point calculation 280# 281# Parameters 282# ---------- 283# mesh : [24, 24, 24] 284# kpt : [1, 1, 2] 285# Returns 286# ------- 287# SCF : -8.65192329453 Hartree per cell 288# CCSD : -0.15529836941 Hartree per cell 289# CCSD(T) : -0.00191451068 Hartree per cell 290 291if __name__ == '__main__': 292 from pyscf.pbc import gto 293 from pyscf.pbc import cc 294 295 cell = gto.Cell() 296 cell.atom = ''' 297 C 0.000000000000 0.000000000000 0.000000000000 298 C 1.685068664391 1.685068664391 1.685068664391 299 ''' 300 cell.basis = 'gth-szv' 301 cell.pseudo = 'gth-pade' 302 cell.a = ''' 303 0.000000000, 3.370137329, 3.370137329 304 3.370137329, 0.000000000, 3.370137329 305 3.370137329, 3.370137329, 0.000000000''' 306 cell.unit = 'B' 307 cell.verbose = 5 308 cell.mesh = [24, 24, 24] 309 cell.build() 310 311 kpts = cell.make_kpts([1, 1, 3]) 312 kpts -= kpts[0] 313 kmf = scf.KRHF(cell, kpts=kpts, exxdiv=None) 314 ehf = kmf.kernel() 315 316 mycc = cc.KGCCSD(kmf) 317 ecc, t1, t2 = mycc.kernel() 318 319 energy_t = kernel(mycc) 320