How are 1-electron and 2-electron integrals done in PySCF?

Question

I am right now using the PySCF for some quantum chemistry problems. I need to get the parameters, which are labeled $h_{ij}$ and $V_{ijkl}$. How can I get this done with PySCF?

Answer 1

While this earlier question on our site was asking why fewer 2e- integrals were printed by PySCF than expected, and your question is simply asking how to calculate the 1e- and 2e- integrals in the first place (so the two questions are not duplicated), that question gave an excellent example for how to print the 2e- integrals for the He atom in a 6-31g basis set:

from pyscf import gto
mol = gto.M(atom='He 0 0 0', basis='6-31g')
eri = mol.intor('int2e', aosym='s8')
print(eri)

Similarly you can get the 1e- integrals by changing int2e to int1e.

Answer 2

Here is the code that I currently used. i hope it could be helpful for other's research. To obtain the integrals, the first step is to run a Hartree Fock approximation and check if it has converged.

m = M(atom=atoms, basis="sto3g", spin=spin, charge=charge)
hf = m.RHF()
hf.kernel()
  for j in range(10):
    mo = hf.stability()[0]
    if np.allclose(mo, hf.mo_coeff):
       break
    dm = hf.make_rdm1(mo, hf.mo_occ)
    hf = hf.run(dm)

Then, one can extract the one- h1e and two-body h2e integral from the molecular class.

h1e = m.intor("int1e_kin") + m.intor("int1e_nuc")
h2e = m.intor("int2e")
scf_c = hf.mo_coeff
nuclear_repulsion = m.energy_nuc()
constant = nuclear_repulsion

If the molecular system that one considers is too large, the complete active space method can also be used similarly:

mycas = hf.CASCI(len(active_indices), n_elec).run()
h1e, constant = mycas.get_h1cas()
h2e = mycas.get_h2cas()
h2e = pyscf.ao2mo.addons.restore('1', h2e, len(active_indices))
h2e = h2e.transpose(0, 2, 3, 1)

The last step is to rotate the basis to the computed Hartree Fock basis as follows.

# Get the one and two electron integral in the Hatree Fock basis
h1e = scf_c.T @ h1e @ scf_c
for i in range(4):
   h2e = np.tensordot(h2e, scf_c, axes=1).transpose(3, 0, 1, 2)
h2e = h2e.transpose(0, 2, 3, 1)