1#!/usr/bin/env python
2
3'''
4MP2 at an individual k-point
5'''
6
7from functools import reduce
8import numpy
9from pyscf.pbc import gto, scf, mp
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 KMP2 with 2x2x2 k-points
28#
29kpts = cell.make_kpts([2,2,2])
30kmf = scf.KRHF(cell)
31kmf.kpts = kpts
32ehf = kmf.kernel()
33
34mypt = mp.KMP2(kmf)
35mypt.kernel()
36print("KMP2 energy (per unit cell) =", mypt.e_tot)
37
38#
39# The KHF and KMP2 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
46mypt = mp.KMP2(kmf)
47mypt.kernel()
48print("KMP2 energy (per unit cell) =", mypt.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
62mypt = mp.RMP2(mf).run()
63print("RMP2 energy (per unit cell) at k-point =", mypt.e_tot)
64dm1 = mypt.make_rdm1()
65dm2 = mypt.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("RMP2 energy based on MP2 density matrices =", e_tot.real)
71
72
73mf = scf.addons.convert_to_uhf(mf)
74mypt = mp.UMP2(mf).run()
75print("UMP2 energy (per unit cell) at k-point =", mypt.e_tot)
76dm1a, dm1b = mypt.make_rdm1()
77dm2aa, dm2ab, dm2bb = mypt.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("UMP2 energy based on MP2 density matrices =", e_tot.real)
91
92
93mf = scf.addons.convert_to_ghf(mf)
94mypt = mp.GMP2(mf).run()
95print("GMP2 energy (per unit cell) at k-point =", mypt.e_tot)
96dm1 = mypt.make_rdm1()
97dm2 = mypt.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("GMP2 energy based on MP2 density matrices =", e_tot.real)
108
109