How to prevent large imaginary part in fock matrix/ folded supercell?

[efertitta]

I occasionally get fatal errors due to imaginary part of folded Fock matrix above tolerance in the form, such as

[ERROR]   |      Imaginary part in folded Fock matrix: L(2)= 3.32e-03 L(inf)= 2.92e-04 !!!
[WARNING] |      Error in folded MO energies: L(2)= 3.53e-06 L(inf)= 1.39e-06 !
[FATAL]   |      Imaginary part of supercell integrals: 2.68e-04 (tolerance= 1.00e-06)

Sometimes I am able to solve the issue by using tighter parameters for the mean-field, but it is not guaranteed to solve it

What is the best strategy to tackle this?

[abhishekkhedkar09]

Hi @efertitta Could you please share the input for the mean-field calculation? Passing time_reversal_symmetry=True to make_kpts should fix this problem...

[efertitta]

Hi @abhishekkhedkar09 I get this error in different circumstances, sometimes changing the k-mesh or the ke_cutoff and precision or cell volume might fix it, but not always This is an example where it happens

import os
import numpy as np
import pyscf
from pyscf.pbc import gto, scf, df
import vayesta
import vayesta.ewf
from vayesta.mpi import mpi

here = os.path.dirname(__file__)

cell = gto.Cell()

cell.a = np.array([[ 1.40562849,  0.81153999,  4.63648548],[-1.40562849,  0.81153999,  4.63648548],[ 0.,         -1.62307997,  4.63648548]])
cell.atom = [['Li', [0.0, 0.0, 0.0]], ['Co', [0.0, 1.62307997235, 2.31824273875]], ['O', [-6.492169233150939e-19, 2.400065585912509e-11, 3.338364128103741]], ['O', [1.4056284884, 0.8115399861759994, 1.2981213493962593]]]

cell.spin = 0
cell.unit = 'A'
cell.verbose = 4
cell.basis= 'gth-dzvp-molopt-sr'
cell.pseudo='gth-hf-rev'
cell.ke_cutoff = 200.0
cell.precision = 1.e-8
cell.build()

kpts = cell.make_kpts([3]*3)

print('Running HF')
mf = scf.KRHF(cell, kpts)
mf.with_df = df.GDF(cell, kpts)
mf.with_df.auxbasis = 'def2-qzvpp-ri'
mf.level_shift = 0.05
mf.max_cycle = 700
mf.chkfile = os.path.join(here,'hf.chk')
mf.kernel()

with open(os.path.join(here,'energies.txt'), 'a') as f:
        f.write('%6s %16.8f\n' % ('HF', mf.e_tot))

etas = [1.e-3, 1.e-4, 1.e-5, 1.e-6, 1.e-7, 1.e-8, 1.e-9, 1.e-10]
for eta in etas:

  emb = vayesta.ewf.EWF(mf, solver='CCSD', bath_options=dict(bathtype='mp2', threshold=eta))

  with emb.iao_fragmentation() as f:
      f.add_all_atomic_fragments()

  emb.kernel()
  e_wf = emb.get_wf_energy()
  e_dm_incluster_cumulant = emb.get_dm_energy()
  if mpi.is_master:
    with open(os.path.join(here,'energies.txt'), 'a') as f:
      f.write('%.2e %16.8f %16.8f\n' % (
               eta , e_wf, e_dm_incluster_cumulant))

[efertitta]

I tried to use kpts = cell.make_kpts([3]*3, time_reversal_symmetry=True) in the above example but pyscf crashed with the following error. Maybe a problem in version 2.3.0?

Traceback (most recent call last):
  File "/cluster/projects/nn8040k/SEP-QUEST/BIGMAP/LiCoO2/KPOINTS/FORGEORGE/Imag_Fmatrix/Time_reversal/Embedding.py", line 34, in <module>
    mf.kernel()
  File "<string>", line 2, in kernel
  File "/cluster/home/edoardof/.local/lib/python3.11/site-packages/pyscf/scf/hf.py", line 1685, in scf
    kernel(self, self.conv_tol, self.conv_tol_grad,
  File "/cluster/home/edoardof/.local/lib/python3.11/site-packages/pyscf/scf/hf.py", line 128, in kernel
    vhf = mf.get_veff(mol, dm)
          ^^^^^^^^^^^^^^^^^^^^
  File "/cluster/home/edoardof/.local/lib/python3.11/site-packages/pyscf/pbc/scf/khf.py", line 629, in get_veff
    vj, vk = self.get_jk(cell, dm_kpts, hermi, kpts, kpts_band)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/cluster/home/edoardof/.local/lib/python3.11/site-packages/pyscf/pbc/scf/khf_ksymm.py", line 248, in get_jk
    vj, vk = self.with_df.get_jk(dm_kpts, hermi, kpts.kpts, kpts_band,
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/cluster/home/edoardof/.local/lib/python3.11/site-packages/pyscf/pbc/df/df.py", line 434, in get_jk
    vk = df_jk.get_k_kpts(self, dm, hermi, kpts, kpts_band, exxdiv)
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/cluster/home/edoardof/.local/lib/python3.11/site-packages/pyscf/pbc/df/df_jk.py", line 178, in get_k_kpts
    mydf.build(kpts_band=kpts_band)
  File "/cluster/home/edoardof/.local/lib/python3.11/site-packages/pyscf/pbc/df/df.py", line 259, in build
    self._make_j3c(self.cell, self.auxcell, None, cderi)
  File "/cluster/home/edoardof/.local/lib/python3.11/site-packages/pyscf/pbc/df/df.py", line 275, in _make_j3c
    kpts_union = unique(numpy.vstack([self.kpts, self.kpts_band]))[0]
                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/cluster/home/edoardof/.local/lib/python3.11/site-packages/numpy/core/shape_base.py", line 289, in vstack
    return _nx.concatenate(arrs, 0, dtype=dtype, casting=casting)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

[efertitta]

Hi @abhishekkhedkar09 I digged into it and it turns out that the above is a bug in pyscf. In pyscf/pbc/df/df.py at line 275 kpts_union = unique(numpy.vstack([self.kpts, self.kpts_band]))[0] is called which is supposed to stack the two arrays self.kpts and self.kpts_band However if time_reversal_symmetry is called upon in cell.make_kpts, self.kpts is created via a different library and it is returned as a KPoints class object rather than an array I added one line just above line 275 in pyscf/pbc/df/df.py: self.kpts = self.kpts.kpts This calls for the kpts array stored in the object and everything runs smoothly

As a bonus, time reversal symmetry now works and the issue with the imaginary component is solved as you suggested!

Since this appears to be a bug in pyscf rather than a problem in Vayesta, I will open a case in pyscf repo and I think you can close this issue

[efertitta]

The error is now reported here https://github.com/pyscf/pyscf/issues/2155#issue-2224671702

[abhishekkhedkar09]

@efertitta That is great... Thank you for hunting down the bug! hopefully PySCF will get this sorted soon as it seems straightforward to rectify it.

[ghb24]

This appears to now be fixed in pyscf master.