1#!/usr/bin/env python 2# Copyright 2014-2019 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# Author: Qiming Sun <osirpt.sun@gmail.com> 17# 18 19''' 20RCCSD 21 22Ref: JCP 90, 1752 (1989); DOI:10.1063/1.456069 23''' 24 25from functools import reduce 26import numpy 27import scipy.linalg 28from pyscf import lib 29from pyscf import ao2mo 30from pyscf.cc import ccsd_rdm 31from pyscf.grad import ccsd as ccsd_grad 32 33def kernel(cc, t1, t2, l1, l2, eris=None): 34 if eris is None: 35 eris = _ERIS(cc, cc.mo_coeff) 36 mol = cc.mol 37 mo_coeff = cc.mo_coeff 38 mo_energy = cc._scf.mo_energy 39 nao, nmo = mo_coeff.shape 40 nocc = numpy.count_nonzero(cc.mo_occ > 0) 41 mo_e_o = mo_energy[:nocc] 42 mo_e_v = mo_energy[nocc:] 43 with_frozen = not ((cc.frozen is None) 44 or (isinstance(cc.frozen, (int, numpy.integer)) and cc.frozen == 0) 45 or (len(cc.frozen) == 0)) 46 47 d1 = _gamma1_intermediates(cc, t1, t2, l1, l2) 48 d2 = _gamma2_intermediates(cc, t1, t2, l1, l2) 49 50 dm2 = ccsd_rdm._make_rdm2(cc, d1, d2, with_dm1=False, with_frozen=False) 51 eri = ao2mo.restore(1, ao2mo.full(cc.mol, mo_coeff), nmo) 52 Imat = numpy.einsum('jqrs,iqrs->ij', dm2, eri) * -1 53 Ioo = Imat[:nocc,:nocc] 54 Ivv = Imat[nocc:,nocc:] 55 doo, dov, dvo, dvv = d1 56 if with_frozen: 57 OA, VA, OF, VF = index_frozen_active(cc) 58 doo[OF[:,None],OA] = Ioo[OF[:,None],OA] / lib.direct_sum('i-j->ij', mo_e_o[OF], mo_e_o[OA]) 59 doo[OA[:,None],OF] = Ioo[OA[:,None],OF] / lib.direct_sum('i-j->ij', mo_e_o[OA], mo_e_o[OF]) 60 dvv[VF[:,None],VA] = Ivv[VF[:,None],VA] / lib.direct_sum('a-b->ab', mo_e_v[VF], mo_e_v[VA]) 61 dvv[VA[:,None],VF] = Ivv[VA[:,None],VF] / lib.direct_sum('a-b->ab', mo_e_v[VA], mo_e_v[VF]) 62 dm1 = scipy.linalg.block_diag(doo+doo.T, dvv+dvv.T) 63 dm1ao = reduce(numpy.dot, (mo_coeff, dm1, mo_coeff.T)) 64 vj, vk = cc._scf.get_jk(cc.mol, dm1ao) 65 Xvo = reduce(numpy.dot, (mo_coeff[:,nocc:].T, vj*2-vk, mo_coeff[:,:nocc])) 66 Xvo += Imat[:nocc,nocc:].T - Imat[nocc:,:nocc] 67 68 dm1 += ccsd_grad._response_dm1(cc, Xvo, eris) 69 Imat[nocc:,:nocc] = Imat[:nocc,nocc:].T 70 71 h1 =-(mol.intor('int1e_ipkin', comp=3) + 72 mol.intor('int1e_ipnuc', comp=3)) 73 s1 =-mol.intor('int1e_ipovlp', comp=3) 74 #zeta = lib.direct_sum('i-j->ij', mo_energy, mo_energy) 75 eri1 = mol.intor('int2e_ip1', comp=3).reshape(3,nao,nao,nao,nao) 76 eri1 = numpy.einsum('xipkl,pj->xijkl', eri1, mo_coeff) 77 eri1 = numpy.einsum('xijpl,pk->xijkl', eri1, mo_coeff) 78 eri1 = numpy.einsum('xijkp,pl->xijkl', eri1, mo_coeff) 79 g0 = ao2mo.restore(1, ao2mo.full(mol, mo_coeff), nmo) 80 81 de = numpy.empty((mol.natm,3)) 82 for k,(sh0, sh1, p0, p1) in enumerate(mol.offset_nr_by_atom()): 83 mol.set_rinv_origin(mol.atom_coord(k)) 84 vrinv = -mol.atom_charge(k) * mol.intor('int1e_iprinv', comp=3) 85 86# 2e AO integrals dot 2pdm 87 de2 = numpy.zeros(3) 88 for i in range(3): 89 g1 = numpy.einsum('pjkl,pi->ijkl', eri1[i,p0:p1], mo_coeff[p0:p1]) 90 g1 = g1 + g1.transpose(1,0,2,3) 91 g1 = g1 + g1.transpose(2,3,0,1) 92 g1 *= -1 93 hx =(numpy.einsum('pq,pi,qj->ij', h1[i,p0:p1], mo_coeff[p0:p1], mo_coeff) + 94 reduce(numpy.dot, (mo_coeff.T, vrinv[i], mo_coeff))) 95 hx = hx + hx.T 96 sx = numpy.einsum('pq,pi,qj->ij', s1[i,p0:p1], mo_coeff[p0:p1], mo_coeff) 97 sx = sx + sx.T 98 99 fij =(hx[:nocc,:nocc] 100 - numpy.einsum('ij,j->ij', sx[:nocc,:nocc], mo_e_o) * .5 101 - numpy.einsum('ij,i->ij', sx[:nocc,:nocc], mo_e_o) * .5 102 - numpy.einsum('kl,ijlk->ij', sx[:nocc,:nocc], 103 g0[:nocc,:nocc,:nocc,:nocc]) * 2 104 + numpy.einsum('kl,iklj->ij', sx[:nocc,:nocc], 105 g0[:nocc,:nocc,:nocc,:nocc]) 106 + numpy.einsum('ijkk->ij', g1[:nocc,:nocc,:nocc,:nocc]) * 2 107 - numpy.einsum('ikkj->ij', g1[:nocc,:nocc,:nocc,:nocc])) 108 109 fab =(hx[nocc:,nocc:] 110 - numpy.einsum('ij,j->ij', sx[nocc:,nocc:], mo_e_v) * .5 111 - numpy.einsum('ij,i->ij', sx[nocc:,nocc:], mo_e_v) * .5 112 - numpy.einsum('kl,ijlk->ij', sx[:nocc,:nocc], 113 g0[nocc:,nocc:,:nocc,:nocc]) * 2 114 + numpy.einsum('kl,iklj->ij', sx[:nocc,:nocc], 115 g0[nocc:,:nocc,:nocc,nocc:]) 116 + numpy.einsum('ijkk->ij', g1[nocc:,nocc:,:nocc,:nocc]) * 2 117 - numpy.einsum('ikkj->ij', g1[nocc:,:nocc,:nocc,nocc:])) 118 119 if with_frozen: 120 fij[OA[:,None],OF] -= numpy.einsum('ij,j->ij', sx[OA[:,None],OF], mo_e_o[OF]) * .5 121 fij[OA[:,None],OF] += numpy.einsum('ij,i->ij', sx[OA[:,None],OF], mo_e_o[OA]) * .5 122 fij[OF[:,None],OA] -= numpy.einsum('ij,j->ij', sx[OF[:,None],OA], mo_e_o[OA]) * .5 123 fij[OF[:,None],OA] += numpy.einsum('ij,i->ij', sx[OF[:,None],OA], mo_e_o[OF]) * .5 124 fab[VA[:,None],VF] -= numpy.einsum('ij,j->ij', sx[VA[:,None],VF], mo_e_v[VF]) * .5 125 fab[VA[:,None],VF] += numpy.einsum('ij,i->ij', sx[VA[:,None],VF], mo_e_v[VA]) * .5 126 fab[VF[:,None],VA] -= numpy.einsum('ij,j->ij', sx[VF[:,None],VA], mo_e_v[VA]) * .5 127 fab[VF[:,None],VA] += numpy.einsum('ij,i->ij', sx[VF[:,None],VA], mo_e_v[VF]) * .5 128 129 fai =(hx[nocc:,:nocc] 130 - numpy.einsum('ai,i->ai', sx[nocc:,:nocc], mo_e_o) 131 - numpy.einsum('kl,ijlk->ij', sx[:nocc,:nocc], 132 g0[nocc:,:nocc,:nocc,:nocc]) * 2 133 + numpy.einsum('kl,iklj->ij', sx[:nocc,:nocc], 134 g0[nocc:,:nocc,:nocc,:nocc]) 135 + numpy.einsum('ijkk->ij', g1[nocc:,:nocc,:nocc,:nocc]) * 2 136 - numpy.einsum('ikkj->ij', g1[nocc:,:nocc,:nocc,:nocc])) 137 138 f1 = numpy.zeros((nmo,nmo)) 139 f1[:nocc,:nocc] = fij 140 f1[nocc:,nocc:] = fab 141 f1[nocc:,:nocc] = fai 142 f1[:nocc,nocc:] = fai.T 143 de2[i] += numpy.einsum('ij,ji', f1, dm1) 144 de2[i] += numpy.einsum('ij,ji', sx, Imat) 145 de2[i] += numpy.einsum('iajb,iajb', dm2, g1) * .5 146 147 de[k] = de2 148 149 return de 150 151 152class _ERIS: 153 def __init__(self, cc, mo_coeff): 154 nocc = numpy.count_nonzero(cc.mo_occ > 0) 155 eri0 = ao2mo.full(cc._scf._eri, mo_coeff) 156 eri0 = ao2mo.restore(1, eri0, mo_coeff.shape[1]) 157 eri0 = eri0.reshape((mo_coeff.shape[1],)*4) 158 self.oooo = eri0[:nocc,:nocc,:nocc,:nocc].copy() 159 self.ooov = eri0[:nocc,:nocc,:nocc,nocc:].copy() 160 self.ovoo = eri0[:nocc,nocc:,:nocc,:nocc].copy() 161 self.oovo = eri0[:nocc,:nocc,nocc:,:nocc].copy() 162 self.oovv = eri0[:nocc,:nocc,nocc:,nocc:].copy() 163 self.ovov = eri0[:nocc,nocc:,:nocc,nocc:].copy() 164 self.ovvo = eri0[:nocc,nocc:,nocc:,:nocc].copy() 165 self.ovvv = eri0[:nocc,nocc:,nocc:,nocc:].copy() 166 self.vvvv = eri0[nocc:,nocc:,nocc:,nocc:].copy() 167 self.vvvo = eri0[nocc:,nocc:,nocc:,:nocc].copy() 168 self.vovv = eri0[nocc:,:nocc,nocc:,nocc:].copy() 169 self.vvov = eri0[nocc:,nocc:,:nocc,nocc:].copy() 170 self.vvoo = eri0[nocc:,nocc:,:nocc,:nocc].copy() 171 self.voov = eri0[nocc:,:nocc,:nocc,nocc:].copy() 172 self.vooo = eri0[nocc:,:nocc,:nocc,:nocc].copy() 173 self.mo_coeff = mo_coeff 174 self.fock = numpy.diag(cc._scf.mo_energy) 175 176def index_frozen_active(cc): 177 nocc = numpy.count_nonzero(cc.mo_occ > 0) 178 moidx = cc.get_frozen_mask() 179 OA = numpy.where( moidx[:nocc])[0] # occupied active orbitals 180 OF = numpy.where(~moidx[:nocc])[0] # occupied frozen orbitals 181 VA = numpy.where( moidx[nocc:])[0] # virtual active orbitals 182 VF = numpy.where(~moidx[nocc:])[0] # virtual frozen orbitals 183 return OA, VA, OF, VF 184 185def _gamma1_intermediates(cc, t1, t2, l1, l2): 186 d1 = ccsd_rdm._gamma1_intermediates(cc, t1, t2, l1, l2) 187 if cc.frozen is None: 188 return d1 189 nocc = numpy.count_nonzero(cc.mo_occ>0) 190 nvir = cc.mo_occ.size - nocc 191 OA, VA, OF, VF = index_frozen_active(cc) 192 doo = numpy.zeros((nocc,nocc)) 193 dov = numpy.zeros((nocc,nvir)) 194 dvo = numpy.zeros((nvir,nocc)) 195 dvv = numpy.zeros((nvir,nvir)) 196 doo[OA[:,None],OA] = d1[0] 197 dov[OA[:,None],VA] = d1[1] 198 dvo[VA[:,None],OA] = d1[2] 199 dvv[VA[:,None],VA] = d1[3] 200 return doo, dov, dvo, dvv 201 202def _gamma2_intermediates(cc, t1, t2, l1, l2): 203 d2 = ccsd_rdm._gamma2_intermediates(cc, t1, t2, l1, l2) 204 nocc, nvir = t1.shape 205 if cc.frozen is None: 206 dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov = d2 207 dvvov = dovvv.transpose(2,3,0,1) 208 dvvvv = ao2mo.restore(1, d2[1], nvir) 209 return dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov 210 nocc0 = numpy.count_nonzero(cc.mo_occ>0) 211 nvir0 = cc.mo_occ.size - nocc0 212 OA, VA, OF, VF = index_frozen_active(cc) 213 dovov = numpy.zeros((nocc0,nvir0,nocc0,nvir0)) 214 dvvvv = numpy.zeros((nvir0,nvir0,nvir0,nvir0)) 215 doooo = numpy.zeros((nocc0,nocc0,nocc0,nocc0)) 216 doovv = numpy.zeros((nocc0,nocc0,nvir0,nvir0)) 217 dovvo = numpy.zeros((nocc0,nvir0,nvir0,nocc0)) 218 dovvv = numpy.zeros((nocc0,nvir0,nvir0,nvir0)) 219 dooov = numpy.zeros((nocc0,nocc0,nocc0,nvir0)) 220 dovov[OA[:,None,None,None],VA[:,None,None],OA[:,None],VA] = d2[0] 221 dvvvv[VA[:,None,None,None],VA[:,None,None],VA[:,None],VA] = ao2mo.restore(1, d2[1], nvir) 222 doooo[OA[:,None,None,None],OA[:,None,None],OA[:,None],OA] = d2[2] 223 doovv[OA[:,None,None,None],OA[:,None,None],VA[:,None],VA] = d2[3] 224 dovvo[OA[:,None,None,None],VA[:,None,None],VA[:,None],OA] = d2[4] 225 dovvv[OA[:,None,None,None],VA[:,None,None],VA[:,None],VA] = d2[6] 226 dooov[OA[:,None,None,None],OA[:,None,None],OA[:,None],VA] = d2[7] 227 dvvov = dovvv.transpose(2,3,0,1) 228 return dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov 229 230 231if __name__ == '__main__': 232 from pyscf import gto 233 from pyscf import scf 234 from pyscf.cc import ccsd 235 from pyscf import grad 236 237 mol = gto.M( 238 verbose = 0, 239 atom = [ 240 ["O" , (0. , 0. , 0.)], 241 [1 , (0. ,-0.757 , 0.587)], 242 [1 , (0. , 0.757 , 0.587)]], 243 basis = '631g' 244 ) 245 mf = scf.RHF(mol).run() 246 mycc = ccsd.CCSD(mf) 247 ecc, t1, t2 = mycc.kernel() 248 l1, l2 = mycc.solve_lambda() 249 g1 = kernel(mycc, t1, t2, l1, l2) 250 ghf = grad.RHF(mf).grad() 251 print('gcc') 252 print(ghf+g1) 253 print(lib.finger(g1) - -0.042511000925747583) 254#[[ 0 0 1.00950969e-02] 255# [ 0 2.28063353e-02 -5.04754844e-03] 256# [ 0 -2.28063353e-02 -5.04754844e-03]] 257 258 print('-----------------------------------') 259 mol = gto.M( 260 verbose = 0, 261 atom = [ 262 ["O" , (0. , 0. , 0.)], 263 [1 , (0. ,-0.757 , 0.587)], 264 [1 , (0. , 0.757 , 0.587)]], 265 basis = '631g' 266 ) 267 mf = scf.RHF(mol).run() 268 mycc = ccsd.CCSD(mf) 269 mycc.frozen = [0,1,10,11,12] 270 ecc, t1, t2 = mycc.kernel() 271 l1, l2 = mycc.solve_lambda() 272 g1 = kernel(mycc, t1, t2, l1, l2) 273 ghf = grad.RHF(mf).grad() 274 print('gcc') 275 print(ghf+g1) 276 print(lib.finger(g1) - 0.10048468674687236) 277#[[ -7.81105940e-17 3.81840540e-15 1.20415540e-02] 278# [ 1.73095055e-16 -7.94568837e-02 -6.02077699e-03] 279# [ -9.49844615e-17 7.94568837e-02 -6.02077699e-03]] 280 281 r = 1.76 282 mol = gto.M( 283 verbose = 0, 284 atom = '''H 0 0 0; H 0 0 %f''' % r, 285 basis = '631g', 286 unit = 'bohr') 287 mf = scf.RHF(mol) 288 mf.conv_tol = 1e-14 289 ehf0 = mf.scf() 290 ghf = grad.RHF(mf).grad() 291 mycc = ccsd.CCSD(mf) 292 ecc, t1, t2 = mycc.kernel() 293 l1, l2 = mycc.solve_lambda() 294 g1 = kernel(mycc, t1, t2, l1, l2) 295 ghf = grad.RHF(mf).grad() 296 print('ghf') 297 print(ghf) 298 print('gcc') 299 print(g1) # 0.015643667024 300 print('tot') 301 print(ghf+g1) # -0.0708003526454 302 303 mol = gto.M( 304 verbose = 0, 305 atom = '''H 0 0 0; H 0 0 %f''' % (r-.001), 306 basis = '631g', 307 unit = 'bohr') 308 mf = scf.RHF(mol) 309 mf.conv_tol = 1e-14 310 ehf0 = mf.scf() 311 mycc = ccsd.CCSD(mf) 312 ecc0 = mycc.kernel()[0] 313 314 mol = gto.M( 315 verbose = 0, 316 atom = '''H 0 0 0; H 0 0 %f''' % (r+.001), 317 basis = '631g', 318 unit = 'bohr') 319 mf = scf.RHF(mol) 320 mf.conv_tol = 1e-14 321 ehf1 = mf.scf() 322 mycc = ccsd.CCSD(mf) 323 ecc1 = mycc.kernel()[0] 324 print((ehf1-ehf0)*500 - ghf[1,2]) 325 print('decc', (ecc1-ecc0)*500 - g1[1,2]) 326 print('decc', (ehf1+ecc1-ehf0-ecc0)*500 - (ghf[1,2]+g1[1,2])) 327