# Verifying periodicity of crystalline orbital modulus in PySCF

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]

• +1 and welcome to our new community! Thank you for your contribution, and we hope to see much more of you in the future!!! I hope you decide to use a more recognizable name in the future, but it's not (yet) required. Commented Sep 24, 2023 at 21:32
• Crystalline orbitals need not have the same periodicity as the cell - this is why we have to do k point sampling. What happens if you plot the total charge density? That should be periodic with periodicity the length of the cell. However I'm not totally sure this is the case - your 2.0/3.657 is sufficiently close to the bohr <-> Angstrom conversion factor that there could just be a unit conversion problem. To be honest I doubt it as it's not that close, but plot the charge density to rule it out. Commented Sep 25, 2023 at 5:46
• I think that the crystal orbitals should indeed have the periodicity of the cell (see Eq. 1 of arxiv.org/abs/1701.04832). Increasing the number of k-points to say 20 actually increases the precision of the periodicity (updated code attached). I checked that the Angstrom/Bohr conversion does not explain the discrepancy. @IanBush Commented Sep 25, 2023 at 13:29
• See mattermodeling.stackexchange.com/questions/7183/… and in particular "Note that a k-point component of the form 1/n, for any integer n, means that the wavefunction phase will be periodic with a period of nL. In other words, each rational k-point may be interpreted as an assumption that the wavefunction is periodic, but with period nL rather than L." Commented Sep 25, 2023 at 13:45
• I'm sorry, I don't really speak python (give me Fortran or C, or maybe Julia). It's just a fact from the maths that, as explained in the linked question, the orbitals need not have the same periodicity as the lattice. The fact it is fractional suggests a bug, but somone who knows the language will be a better person to take a look than I/ Commented Sep 25, 2023 at 14:06

Setting the unit='Bohr' flag in gto.M (as originally suggested in the comments) resolved the issue in the latest version of the code above. Presumably there is a discrepancy in the default values of units passed to numint functions and gto.M.