1#!/usr/bin/env python
2
3'''
4Looks at a hydrogen metallic lattice and looks at using the level shift in
5k-point ccsd.  While for most systems the level shift will not affect results,
6this is one instance where the system will converge on a different ccsd solution
7depending on the initial guess and whether one is using a level shift.
8'''
9
10import numpy as np
11
12from pyscf.lib import finger
13from pyscf.pbc import gto as pbcgto
14from pyscf.pbc import scf as pbcscf
15
16import pyscf.cc
17import pyscf.pbc.cc as pbcc
18from pyscf.pbc.lib import kpts_helper
19import pyscf.pbc.cc.kccsd_t_rhf as kccsd_t_rhf
20
21cell = pbcgto.Cell()
22cell.atom = [['H', (0.000000000, 0.000000000, 0.000000000)],
23             ['H', (0.000000000, 0.500000000, 0.250000000)],
24             ['H', (0.500000000, 0.500000000, 0.500000000)],
25             ['H', (0.500000000, 0.000000000, 0.750000000)]]
26cell.unit = 'Bohr'
27cell.a = [[1.,0.,0.],[0.,1.,0],[0,0,2.2]]
28cell.verbose = 3
29cell.spin = 0
30cell.charge = 0
31cell.basis = 'gth-szv'
32cell.pseudo = 'gth-pade'
33for i in range(len(cell.atom)):
34    cell.atom[i][1] = tuple(np.dot(np.array(cell.atom[i][1]),np.array(cell.a)))
35cell.build()
36
37nmp = [2, 1, 1]
38
39kmf = pbcscf.KRHF(cell)
40kmf.kpts = cell.make_kpts(nmp, scaled_center=[0.0,0.0,0.0])
41e = kmf.kernel()  # 2.30510338236481
42
43mycc = pbcc.KCCSD(kmf)
44eris = mycc.ao2mo(kmf.mo_coeff)
45eris.mo_energy = [eris.fock[k].diagonal() for k in range(mycc.nkpts)]
46print('\nCCSD energy w/o level shift and MP2 initial guess:')  # 0.02417522810234485
47ekccsd, t1, t2 = mycc.kernel(eris=eris)
48
49# Use a level shift with a level shift equal to the Madelung
50# constant for this system.  Using the previous t1/t2 as an initial
51# guess, we see that these amplitudes still solve the CCSD amplitude
52# equations.
53
54def _adjust_occ(mo_energy, nocc, shift):
55    '''Modify occupied orbital energy'''
56    mo_energy = mo_energy.copy()
57    mo_energy[:nocc] += shift
58    return mo_energy
59
60madelung = 1.36540204381
61eris.mo_energy = [_adjust_occ(mo_e, mycc.nocc, madelung) for mo_e in eris.mo_energy]
62print('\nCCSD energy w/o level shift and previous t1/t2 as initial guess:')  # 0.02417522810234485
63ekccsd, _, _ = mycc.kernel(t1=t1, t2=t2, eris=eris)
64
65# Use level shift with an MP2 guess.  Here the results will differ from
66# those before.
67
68print('\nCCSD energy w/ level shift and MP2 initial guess:')  # -0.11122802032348603
69ekccsd, t1, t2 = mycc.kernel(eris=eris)
70
71# Check to see it satisfies the CCSD amplitude equations.
72
73print('\nCCSD energy w/ level shift and previous t1/t2 as initial guess:')  # -0.11122802032348603
74ekccsd, _, _ = mycc.kernel(t1=t1, t2=t2, eris=eris)
75