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?

  • 1
    $\begingroup$ +1. Welcome to our new community, and we hope to see much more of you in the future!!! Thank you for contributing your question here, and joining us. The "parameters" to which you refer, are called "integrals", so I've changed the question so that more of our community will know what you're seeking. $\endgroup$ Commented Feb 5, 2021 at 19:04

2 Answers 2


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')

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

  • 1
    $\begingroup$ If he specifically wants $V_{ijkl}$ then aosym='s1' is likely a better option to include in your answer. Additionally, there is no such intor named int1e, so that will raise an Exception - if it is the 'core' Hamiltonian (nuclear + kinetic), then use mol.get_hcore(). That is equivalent to mol.intor('int1e_nuc') + mol.intor('int1e_kin'). $\endgroup$
    – obackhouse
    Commented Jul 11, 2021 at 11:20

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()
  for j in range(10):
    mo = hf.stability()[0]
    if np.allclose(mo, hf.mo_coeff):
    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)

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .