eval_gto_precision should be a little smarter

[lkwagner]

If you have a small cell with funny angles, the default eval_gto_precision might not be enough. This shows it:

It should be possible to automatically set this or at least have a test suite to make sure your variables are set correctly.

import pyqmc.api as pyq
import helper
import numpy as np
import matplotlib.pyplot as plt
import jax
import pyqmc.jax.jastrowspin as jax_jastrow
import pyqmc.jastrowspin as pyq_jastrow
import h5py
import pyqmc.jax.slater
from pyscf.scf import addons

cpu = True
if cpu:
    jax.config.update("jax_platform_name", "cpu")
    jax.config.update("jax_enable_x64", True)
else:
    pass


cell, mf = pyq.recover_pyscf(f"pbc_hf.hdf5", cancel_outputs=False)
cell_cart, mf_cart = pyq.recover_pyscf(f"unc-pbc_hf.hdf5", cancel_outputs=False)


wfs = {}
for precision in [1e-12, 1e-6, 1e-3, 1e-2]:
    wf, _ = pyq.generate_slater(cell, mf, eval_gto_precision=precision)
    wfs[precision] = wf
coords = pyq.initial_guess(cell, 1)

# For the same coordinates in my plots posted on slack
x = np.array(
     [[[1.4830396,  5.31432205, 3.92031976],
   [9.2106652,  4.41300995, 4.82653279],
   [2.23978882, 1.54766513, 3.21768169],
   [2.94263674, 3.18792307, 0.51652064],
   [7.23572943, 7.290782,   9.30377819],
   [0.81286747, 0.59370649, 0.87913851],
   [2.87981611, 3.65634976, 2.50512875],
   [3.06987531, 1.55142514, 2.02394672]]]
 )

xpbc, _ = pyqmc.pbc.enforce_pbc(cell.lattice_vectors(), x)
coords.configs = np.tile(xpbc, (coords.configs.shape[0], 1, 1))
lattice_vecs = cell.lattice_vectors()


npoints = 1000
all_epos = np.array([np.linspace(0, lattice_vecs[2, 0], npoints), 
                 np.linspace(0, lattice_vecs[2, 1], npoints), 
                 np.linspace(0, lattice_vecs[2, 2], npoints)]).T

wfvals = {}
for nm in wfs.keys():
    wfvals[nm] = np.zeros(npoints)
    wfs[nm].recompute(coords)

e = 5
for i in range(npoints):
    epos = np.tile(all_epos[i], (coords.configs.shape[0], 1))
    epos = coords.make_irreducible(e, epos)
    for nm in wfs.keys():
        wfvals[nm][i] = wfs[nm].testvalue(e, epos)[0][0].real

plt.figure()
for nm, item in wfvals.items():
    plt.plot(item, label=nm, alpha=0.75)
plt.legend()
plt.ylabel(r"$\psi$")
plt.tight_layout()
plt.savefig("pbc_wfval.png")
plt.show()

[lkwagner]

pbc_hf.hdf5.zip