I'm trying to verify that the modulus of the crystalline atomic/molecular orbitals computed by PySCF satisfy the correct periodicity condition for the supercell. I'm considering a hydrogen chain oriented along the z-axis with H-H separation 1.0 and with two H atoms per cell, giving a cell size of 2.0.
The following code evaluates the modulus of the atomic/molecular orbitals for all k-points at positions $(0,0,0)$ and $(0,0,2)$.
One can see from the output of the code that they are far from equality. How can one understand this bizarre behavior?
import numpy as np
import matplotlib.pyplot as plt
from pyscf.pbc import gto, dft
from pyscf.pbc.dft import gen_grid, numint
num_kpts = 20
cell = gto.M(atom=[['H', (0.0, 0.0, 1.0)], ['H', (0.0, 0.0, 2.0)]],
a=[[40.0, 0.0, 0.0], [0., 40.0, 0.0], [0.0, 0.0, 2.0]],
verbose=1)
kpts = cell.make_kpts([1, 1, num_kpts])
mf = dft.KRKS(cell)
mf.kpts = kpts
mf.xc = 'lda,vwn'
mf = mf.density_fit()
energy = mf.kernel()
dm = cell.make_rdm1()
orbital_index = 0 # select an orbital for evaluation
points = [(0.0, 0.0, 0.0), (0.0, 0.0, 2.0)] # modulus should be equal at these points
print('orbital: ' + str(orbital_index))
for kpt_index in range(num_kpts):
print('k-point: ' + str(kpt_index))
kpt = kpts[kpt_index]
ao = numint.eval_ao(cell, points, kpt=kpt) # (2, nao)
rho = numint.eval_rho(cell, ao, dm).real # (2,)
print('rho: ' + str(rho))
ao_density = np.abs(ao[:, orbital_index])**2
print('ao_density: ' + str(ao_density))
mo_coeff = mf.mo_coeff # (num_kpts, nao, nao)
mo = np.einsum('xi,kij->xkj', ao, mo_coeff) # (2, num_kpts, nao)
mo_density = np.abs(mo[:, kpt_index, orbital_index])**2
print('mo_density: ' + str(mo_density))
# orbital: 0
# k-point: 0
# rho: [0.27688535 0.26870818]
# ao_density: [0.02243983 0.39680968]
# mo_density: [0.13844267 0.13435409]
# k-point: 1
# rho: [0.27400583 0.26613968]
# ao_density: [0.02186969 0.39587295]
# mo_density: [0.13927244 0.13514004]
# k-point: 2
# rho: [0.26569365 0.25872769]
# ao_density: [0.0202211 0.39316138]
# mo_density: [0.14179107 0.13752496]
# k-point: 3
# rho: [0.25287852 0.24730755]
# ao_density: [0.01767116 0.38895845]
# mo_density: [0.14608699 0.14159039]
# k-point: 4
# rho: [0.23695731 0.23313186]
# ao_density: [0.01448881 0.38369784]
# mo_density: [0.15230891 0.14747331]
# k-point: 5
# rho: [0.21960232 0.21769588]
# ao_density: [0.01100108 0.37791244]
# mo_density: [0.16066695 0.15536683]
# k-point: 6
# rho: [0.20255418 0.20255011]
# ao_density: [0.00755514 0.37217532]
# mo_density: [0.17143259 0.16551978]
# k-point: 7
# rho: [0.18743636 0.18913423]
# ao_density: [0.00448221 0.36704113]
# mo_density: [0.18493519 0.17823305]
# k-point: 8
# rho: [0.17561427 0.17865325]
# ao_density: [0.00206752 0.36299446]
# mo_density: [0.20155235 0.19385034]
# k-point: 9
# rho: [0.16810548 0.17200107]
# ao_density: [0.00052834 0.3604093 ]
# mo_density: [0.22168984 0.21274038]
# k-point: 10
# rho: [0.16553284 0.1697228 ]
# ao_density: [1.51068664e-28 3.59520862e-01]
# mo_density: [0.24574325 0.25196794]
# k-point: 11
# rho: [0.16810548 0.17200107]
# ao_density: [0.00052834 0.3604093 ]
# mo_density: [0.22168984 0.21274038]
# k-point: 12
# rho: [0.17561427 0.17865325]
# ao_density: [0.00206752 0.36299446]
# mo_density: [0.20155235 0.19385034]
# k-point: 13
# rho: [0.18743636 0.18913423]
# ao_density: [0.00448221 0.36704113]
# mo_density: [0.18493519 0.17823305]
# k-point: 14
# rho: [0.20255418 0.20255011]
# ao_density: [0.00755514 0.37217532]
# mo_density: [0.17143259 0.16551978]
# k-point: 15
# rho: [0.21960232 0.21769588]
# ao_density: [0.01100108 0.37791244]
# mo_density: [0.16066695 0.15536683]
# k-point: 16
# rho: [0.23695731 0.23313186]
# ao_density: [0.01448881 0.38369784]
# mo_density: [0.15230891 0.14747331]
# k-point: 17
# rho: [0.25287852 0.24730755]
# ao_density: [0.01767116 0.38895845]
# mo_density: [0.14608699 0.14159039]
# k-point: 18
# rho: [0.26569365 0.25872769]
# ao_density: [0.0202211 0.39316138]
# mo_density: [0.14179107 0.13752496]
# k-point: 19
# rho: [0.27400583 0.26613968]
# ao_density: [0.02186969 0.39587295]
# mo_density: [0.13927244 0.13514004]
