1#!/usr/bin/env python
2
3'''
4Showing use of the parallelized CCSD with K-point sampling.
5'''
6
7import numpy as np
8from pyscf import cc
9from pyscf.pbc import cc as pbccc
10from pyscf.pbc import scf as pbchf
11from pyscf.pbc import gto
12from pyscf.pbc.tools.pbc import super_cell
13
14nmp = [1, 1, 2]
15cell = gto.M(
16    unit='B',
17    a=[[0., 3.37013733, 3.37013733],
18       [3.37013733, 0., 3.37013733],
19       [3.37013733, 3.37013733, 0.]],
20    mesh=[24,]*3,
21    atom='''C 0 0 0
22              C 1.68506866 1.68506866 1.68506866''',
23    basis='gth-szv',
24    pseudo='gth-pade',
25    verbose=4
26)
27
28# We build a supercell composed of 'nmp' replicated units and run
29# our usual molecular Hartree-Fock program, but using integrals
30# between periodic gaussians.
31#cell = build_cell(ase_atom, ke=50., basis=basis)
32supcell = super_cell(cell, nmp)
33mf = pbchf.RHF(supcell)
34mf.kernel()
35supcell_energy = mf.energy_tot() / np.prod(nmp)
36
37mycc = cc.RCCSD(mf)
38gccsd_energy = mycc.ccsd()[0] / np.prod(nmp)
39eip, wip = mycc.ipccsd(nroots=2)
40eea, wea = mycc.eaccsd(nroots=2)
41
42# We now begin our k-point calculations for the same system, making
43# sure we shift the k-points to be gamma-centered.
44kpts = cell.make_kpts(nmp)
45kpts -= kpts[0]
46kmf = pbchf.KRHF(cell, kpts)
47kpoint_energy = kmf.kernel()
48
49mykcc = pbccc.KRCCSD(kmf)
50kccsd_energy = mykcc.ccsd()[0]
51ekcc = mykcc.ecc
52# We look at a gamma-point transition for IP/EA
53ekip, wkip = mykcc.ipccsd(nroots=2, kptlist=[0])
54ekea, wkea = mykcc.eaccsd(nroots=2, kptlist=[0])
55
56print('Difference between gamma/k-point mean-field calculation = %.15g' % (
57    abs(supcell_energy-kpoint_energy)))
58print('Difference between gamma/k-point ccsd calculation = %.15g' % (
59    abs(gccsd_energy - kccsd_energy)))
60print('Difference between gamma/k-point ip-eomccsd calculation = %.15g' % (
61    np.linalg.norm(np.array(eip) - np.array(ekip))))
62print('Difference between gamma/k-point ea-eomccsd calculation = %.15g' % (
63    np.linalg.norm(np.array(eea) - np.array(ekea))))
64