JoonhoLee-Group / ipie

ipie stands for Intelligent Python-based Imaginary-time Evolution with a focus on simplicity and speed.
Apache License 2.0
43 stars 22 forks source link

Building Hamiltonian for some molecules results in an error. #313

Open Koberard opened 1 week ago

Koberard commented 1 week ago

Describe the bug When running generate_hamiltonian() I get that the modified Cholesky calculation goes out of bounds.

To Reproduce Perform UHF in pyscf. Then build AFQMC. import numpy as np from pyscf import gto, scf from ipie.analysis.autocorr import reblock_by_autocorr from ipie.config import MPI

We can extract the qmc data as as a pandas data frame like so

from ipie.analysis.extraction import extract_observable from ipie.estimators.energy import EnergyEstimator, local_energy_batch from ipie.qmc.afqmc import AFQMC from ipie.systems.generic import Generic from ipie.trial_wavefunction.single_det import SingleDet from ipie.utils.backend import arraylib as xp from ipie.utils.from_pyscf import generate_hamiltonian, generate_wavefunction_from_mo_coeff from ipie.qmc.calc import build_afqmc_driver from ipie.analysis.extraction import extract_observable from ipie.estimators.estimator_base import EstimatorBase from ipie.estimators.greens_function import greens_function from ipie.estimators.local_energy import local_energy_G import os Dir = os.getcwd()+'/Mols/' Mol = os.listdir(Dir) print(Mol) for i in Mol: molecule = Dir + i

print(molecule)
mol = gto.M(
    atom=molecule,
    basis="sto-3g",
    verbose=4,
    unit="Angstrom",

)
mf = scf.UHF(mol)
mf.max_cycle=200
mf.chkfile = "scf.chk"
mf.kernel()
comm = MPI.COMM_WORLD
# Now let's build our custom AFQMC algorithm
num_walkers = 100
num_steps_per_block = 25
num_blocks = 10
timestep = 0.000005
# 1. Build Hamiltonian
ham = generate_hamiltonian(
    mol,
    mf.mo_coeff,
    mf.get_hcore(),
    mf.mo_coeff,  # should be optional
)
# 2. Build trial wavefunction
orbs = generate_wavefunction_from_mo_coeff(
    mf.mo_coeff,
    mf.mo_occ,
    mf.mo_coeff,  # Make optional argument
    mol.nelec,
)
num_basis = mf.mo_coeff[0].shape[-1]
trial = NoisySingleDet(
    np.hstack([orbs, orbs]),
    mol.nelec,
    num_basis,
    noise_level=0,
)
# TODO: would be nice if we didn't need to do this.
trial.build()
trial.half_rotate(ham)
afqmc = AFQMC.build(
    mol.nelec,
    ham,
    trial,
    num_walkers=num_walkers,
    num_steps_per_block=num_steps_per_block,
    num_blocks=num_blocks,
    timestep=timestep,
    seed=59306159,
)

afqmc.run(additional_estimators=add_est)

Expected behavior The expected behavior is the Hamiltonian should be built.

Screenshots Output can be seen in Jupiter notebook

Other information:

[e.g. macOS Monterey]

Additional context Add any other context about the problem here. Bug_1.zip

fdmalone commented 1 week ago

could you specify the molecule in a comment like so:

mol = gto.M(
    atom=molecule, # <<<< we would need this to debug
    basis="sto-3g",
    verbose=4,
    unit="Angstrom",
)
Koberard commented 1 week ago

The bug only works with specific molecules for some reason. I provided the XYZ file of the molecule in the Mols folder which is read in via

import os Dir = os.getcwd()+'/Mols/' Mol = os.listdir(Dir) print(Mol) for i in Mol: molecule = Dir + i

print(molecule)
mol = gto.M(
    atom=molecule, ...

The contents of the XYZ file is, 18 p-xylene[B3LYP/aug-cc-pVTZ] C -0.006126 0.695514 1.193928 C 0.006126 -0.695514 1.193928 C -0.005956 0.693932 -1.191741 C 0.005956 -0.693932 -1.191741 C -0.010665 1.417282 0.002438 C 0.010665 -1.417282 0.002438 C 0.005956 2.923660 -0.002240 C -0.005956 -2.923660 -0.002240 H -0.011287 1.226083 2.138293 H 0.011287 -1.226083 2.138293 H -0.010880 1.224834 -2.136140 H 0.010880 -1.224834 -2.136140 H -0.642816 3.326662 -0.780885 H -0.324391 3.327444 0.954402 H 1.012311 3.306699 -0.189982 H 0.642816 -3.326662 -0.780885 H 0.324391 -3.327444 0.954402 H -1.012311 -3.306699 -0.189982 .

Koberard commented 1 week ago

My software versions are the following.

OS: Ubuntu 22.04.3 LTS python 3.11.7 Version -e git+https://github.com/linusjoonho/ipie.git@80a2444e19b846516fbf62962b7953e3a5de6137#egg=ipie

fdmalone commented 6 days ago

From #314 it looks like this is an issue with passing UHF mo coeffs in:

ham = generate_hamiltonian(
    mol,
    mf.mo_coeff,
    mf.get_hcore(),
    mf.mo_coeff,  # <<< should raise an error if these are UHF

@linusjoonho could someone look at this example and 1) confirm that it's the issue and 2) put in some assertions / raise and error if UHF integrals are found.

More generally we should try to make it as easy as possible to generate from a pyscf molecule from python, one option would be to add another factory method, or just some function which builds the objects from pyscf (trial, hamiltonian etc,) to remove some steps.