Imaginary vaules when unfolding AO integrals.

[basilib]

Silicon calculation with def2-svp crashes due to imaginary values in unfolded AO integrals. Traceback and script below:

Traceback (most recent call last): File "lithium.py", line 73, in edft, erhf, eemb = run(a) File "lithium.py", line 53, in run emb = vayesta.ewf.EWF(kmf, solver="CCSD", bno_threshold=1e-6) File "/users/git/Vayesta/vayesta/ewf/init.py", line 17, in EWF return REWF(mf, *args, kwargs) File "/users/git/Vayesta/vayesta/ewf/ewf.py", line 98, in init super().init(mf, options=options, log=log, kwargs) File "/users/git/Vayesta/vayesta/core/qemb/qemb.py", line 152, in init self._veff_orig = self.get_init_veff() File "/users/git/Vayesta/vayesta/core/qemb/qemb.py", line 272, in get_init_veff return self.mf.get_veff() File "/users/git/Vayesta/vayesta/core/foldscf.py", line 85, in get_veff veff = k2bvk_2d(vk, self.kphase) File "/users/git/Vayesta/vayesta/core/foldscf.py", line 208, in k2bvk_2d assert abs(ag.imag).max() <= imag_tol AssertionError

import numpy as np
import pyscf
import pyscf.pbc
import pyscf.pbc.scf
import pyscf.pbc.cc
import pyscf.pbc.tools
import vayesta
import vayesta.ewf

erhfs = []
edfts = []
eembs = []
def run(a):
    print(a/0.529177)
    cell = pyscf.pbc.gto.Cell()
    cell.atom = ['Si 0.0 0.0 0.0', 'Si %f %f %f' % (a/4, a/4, a/4)]
    cell.a = np.asarray([
        [a/2, a/2, 0],
        [0, a/2, a/2],
        [a/2, 0, a/2]])
    cell.basis = 'def2-svp'
    cell.output = 'si.log'
    cell.ke_cutoff = 80
    cell.exp_to_discard=0.1
    cell.build()

    kmesh = [3,3,3]
    kpts = cell.make_kpts(kmesh)

    print("Running RHF")
    kmf = pyscf.pbc.scf.KRHF(cell, [kpts])

    kmf = kmf.density_fit()

    kmf.kernel()
    print(kmf.e_tot)

    print("Running Vayesta")
    vayesta.new_log("vayesta-si-%.1f.log" % a)
    emb = vayesta.ewf.EWF(kmf, solver="CCSD", bno_threshold=1e-6)
    emb.iao_fragmentation()
    emb.add_atomic_fragment(0, sym_factor=2) # 2 Si-atoms per unit cell
    emb.kernel()
    print(emb.e_tot)

    return kmf.e_tot, emb.e_tot

rng = np.arange(9.6, 11.3, 0.1)
for a in rng*0.529177:
    erhf, eemb = run(a)
    erhfs.append(erhf)
    eembs.append(eemb)

[maxnus]

I am not sure what is happening here... I think this would also occur for PySCF, using their k2gamma tools, so maybe we should consider making a PySCF minimal example and post it there

[maxnus]

Not much progress on this. The critical value of lattice constant seems to be between 9.6 and 9.7 B - for 9.7 the make real works, for 9.6 it doesn't. I think the mean-field solution plays a role here; it seem related to the HOMO-LUMO gap, for the 9.6 B lattice constant the gap is only ~0.04 Htr, whereas it's more like 0.5 Htr for 9.7 B. Additionally, there might also be some crossing of MO energies between these solutions. Nevertheless, I would have thought that a gap of 40 mHtr is still plenty large enough to cleanly separate it from the virtual space...

[maxnus]

This is a consequence of the MF breaking k -- (-k) conjugation symmetry, in which case the solution cannot be represented by real orbitals, see https://github.com/pyscf/pyscf/issues/1161

If this solution is desired, complex orbitals are needed, which will require quite some work (but we will want comlex orbitals in the long run anyways).

Alternatively, one can try to instead find a MF solution with k -- (-k) conjugation symmetry. This, for example, is the default in VASP (see https://www.vasp.at/wiki/index.php/ISYM). This will also save computation cost, as k and -k do not need to be treated independently. Unfortunately, this symmetry is not implemented in PySCF. It may be possible to force such a solution by symmetrizing between k and -k each SCF iteration, although this would not reduce the computational cost.

[basilib]

Importing from pyscf.pbc.lib.kpts_helper import conj_mapping and from pyscf.pbc.scf.addons import kconj_symmetry_ and adding kmf = kconj_symmetry_(kmf) before the density fitting line in the above will result in a gapless non-converged SCF and vayesta crashes with

Traceback (most recent call last): File "minimal.py", line 99, in erhf, eemb = run(a) File "minimal.py", line 89, in run emb = vayesta.ewf.EWF(kmf, solver="CCSD", bno_threshold=1e-6) File "/users/k21127360/git/Vayesta/vayesta/ewf/init.py", line 17, in EWF return REWF(mf, args, kwargs) File "/users/k21127360/git/Vayesta/vayesta/ewf/ewf.py", line 105, in init super().init(mf, options=options, log=log, kwargs) File "/users/k21127360/git/Vayesta/vayesta/core/qemb/qemb.py", line 143, in init self.init_mf(mf) File "/users/k21127360/git/Vayesta/vayesta/core/qemb/qemb.py", line 188, in init_mf mf = fold_scf(mf) File "/users/k21127360/git/Vayesta/vayesta/core/foldscf.py", line 30, in fold_scf return FoldedRHF(kmf, args, **kwargs) File "/users/k21127360/git/Vayesta/vayesta/core/foldscf.py", line 138, in init fold_mos(kmf.mo_energy, kmf.mo_coeff, kmf.mo_occ, self.kphase, ovlp, hcore) File "/users/k21127360/git/Vayesta/vayesta/core/foldscf.py", line 189, in fold_mos mo_energy = np.hstack(kmo_energy) File "<__array_function__ internals>", line 5, in hstack TypeError: dispatcher for __array_function__ did not return an iterable