1#!/usr/bin/env python 2 3''' 4CCSD with k-point sampling or at an individual k-point 5''' 6 7from functools import reduce 8import numpy 9from pyscf.pbc import gto, scf, cc 10 11cell = gto.Cell() 12cell.atom=''' 13C 0.000000000000 0.000000000000 0.000000000000 14C 1.685068664391 1.685068664391 1.685068664391 15''' 16cell.basis = 'gth-szv' 17cell.pseudo = 'gth-pade' 18cell.a = ''' 190.000000000, 3.370137329, 3.370137329 203.370137329, 0.000000000, 3.370137329 213.370137329, 3.370137329, 0.000000000''' 22cell.unit = 'B' 23cell.verbose = 5 24cell.build() 25 26# 27# KHF and KCCSD with 2x2x2 k-points 28# 29kpts = cell.make_kpts([2,2,2]) 30kmf = scf.KRHF(cell) 31kmf.kpts = kpts 32ehf = kmf.kernel() 33 34mycc = cc.KCCSD(kmf) 35mycc.kernel() 36print("KRCCSD energy (per unit cell) =", mycc.e_tot) 37 38# 39# The KHF and KCCSD for single k-point calculation. 40# 41kpts = cell.get_abs_kpts([0.25, 0.25, 0.25]) 42kmf = scf.KRHF(cell) 43kmf.kpts = kpts 44ehf = kmf.kernel() 45 46mycc = cc.KRCCSD(kmf) 47mycc.kernel() 48print("KRCCSD energy (per unit cell) =", mycc.e_tot) 49 50 51# 52# The PBC module provides an separated implementation specified for the single 53# k-point calculations. They are more efficient than the general implementation 54# with k-point sampling. For gamma point, integrals and orbitals are all real 55# in this implementation. They can be mixed with other post-HF methods that 56# were provided in the molecular program. 57# 58kpt = cell.get_abs_kpts([0.25, 0.25, 0.25]) 59mf = scf.RHF(cell, kpt=kpt) 60ehf = mf.kernel() 61 62mycc = cc.RCCSD(mf).run() 63print("RCCSD energy (per unit cell) at k-point =", mycc.e_tot) 64dm1 = mycc.make_rdm1() 65dm2 = mycc.make_rdm2() 66nmo = mf.mo_coeff.shape[1] 67eri_mo = mf.with_df.ao2mo(mf.mo_coeff, kpts=kpt).reshape([nmo]*4) 68h1 = reduce(numpy.dot, (mf.mo_coeff.conj().T, mf.get_hcore(), mf.mo_coeff)) 69e_tot = numpy.einsum('ij,ji', h1, dm1) + numpy.einsum('ijkl,jilk', eri_mo, dm2)*.5 + mf.energy_nuc() 70print("RCCSD energy based on CCSD density matrices =", e_tot.real) 71 72 73mf = scf.addons.convert_to_uhf(mf) 74mycc = cc.UCCSD(mf).run() 75print("UCCSD energy (per unit cell) at k-point =", mycc.e_tot) 76dm1a, dm1b = mycc.make_rdm1() 77dm2aa, dm2ab, dm2bb = mycc.make_rdm2() 78nmo = dm1a.shape[0] 79eri_aa = mf.with_df.ao2mo(mf.mo_coeff[0], kpts=kpt).reshape([nmo]*4) 80eri_bb = mf.with_df.ao2mo(mf.mo_coeff[1], kpts=kpt).reshape([nmo]*4) 81eri_ab = mf.with_df.ao2mo((mf.mo_coeff[0],mf.mo_coeff[0],mf.mo_coeff[1],mf.mo_coeff[1]), kpts=kpt).reshape([nmo]*4) 82hcore = mf.get_hcore() 83h1a = reduce(numpy.dot, (mf.mo_coeff[0].conj().T, hcore, mf.mo_coeff[0])) 84h1b = reduce(numpy.dot, (mf.mo_coeff[1].conj().T, hcore, mf.mo_coeff[1])) 85e_tot = (numpy.einsum('ij,ji', h1a, dm1a) + 86 numpy.einsum('ij,ji', h1b, dm1b) + 87 numpy.einsum('ijkl,jilk', eri_aa, dm2aa)*.5 + 88 numpy.einsum('ijkl,jilk', eri_ab, dm2ab) + 89 numpy.einsum('ijkl,jilk', eri_bb, dm2bb)*.5 + mf.energy_nuc()) 90print("UCCSD energy based on CCSD density matrices =", e_tot.real) 91 92 93mf = scf.addons.convert_to_ghf(mf) 94mycc = cc.GCCSD(mf).run() 95print("GCCSD energy (per unit cell) at k-point =", mycc.e_tot) 96dm1 = mycc.make_rdm1() 97dm2 = mycc.make_rdm2() 98nao = cell.nao_nr() 99nmo = mf.mo_coeff.shape[1] 100mo = mf.mo_coeff[:nao] + mf.mo_coeff[nao:] 101eri_mo = mf.with_df.ao2mo(mo, kpts=kpt).reshape([nmo]*4) 102orbspin = mf.mo_coeff.orbspin 103eri_mo[orbspin[:,None]!=orbspin] = 0 104eri_mo[:,:,orbspin[:,None]!=orbspin] = 0 105h1 = reduce(numpy.dot, (mf.mo_coeff.conj().T, mf.get_hcore(), mf.mo_coeff)) 106e_tot = numpy.einsum('ij,ji', h1, dm1) + numpy.einsum('ijkl,jilk', eri_mo, dm2)*.5 + mf.energy_nuc() 107print("GCCSD energy based on CCSD density matrices =", e_tot.real) 108