0

I am trying to create a package to do a modified version of FCI. I am trying to learn how to index determinants, solve equations using slater rules, etc. I am currently stuck on the 2 electron integrals. I need to get the set of integrals for my basis set that will be used in solving for slater-condon rules for my system. I currently am running a few tests. Progress below:

import pyscf as ps
import pyscf.mcscf as mc
import pyscf.tools.molden as md
import pyscf.lo as lo
import scipy
import pyscf.scf.ghf as ghf
import pyscf.ao2mo as ao2mo
import numpy as np

print("yes")

mol_h2o = ps.gto.M(atom = 'H 0 1 0; H 0 0 1', basis = 'ccpvdz', symmetry=False)

uhf_h2o = ps.scf.UHF(mol_h2o)
uhf_h2o.kernel()
print(mol_h2o.intor('int1e_kin').shape)
print(mol_h2o.intor('int1e_nuc').shape)
print(mol_h2o.intor('int1e_ovlp').shape)
print(mol_h2o.intor('int2e').shape)
print(uhf_h2o.mo_coeff.shape)
print(uhf_h2o.mo_occ.shape)

eri1 = mol_h2o.ao2mo(uhf_h2o.mo_coeff[0], uhf_h2o.mo_coeff[0], aosym = 's1', compact=False)
print(eri1.shape)

H=np.zeros([2,2])

N = 10
sum=0
for i in np.arange(0,N):
    if(uhf_h2o.mo_occ[0,i]==1):
        for j in np.arange(0,N):
            for k in np.arange(0,N):
                #sum = sum+ uhf_h2o.mo_coeff[0,i,j]*uhf_h2o.mo_coeff[0,i,k]*(mol_h2o.intor('int1e_kin')[j,k]+mol_h2o.intor('int1e_nuc')[j,k])
                sum = sum+ uhf_h2o.mo_coeff[0,j,i]*uhf_h2o.mo_coeff[0,k,i]*(mol_h2o.intor('int1e_kin')[j,k]+mol_h2o.intor('int1e_nuc')[j,k])
    if(uhf_h2o.mo_occ[1,i]==1):
        for j in np.arange(0,N):
            for k in np.arange(0,N):
                #sum = sum+ uhf_h2o.mo_coeff[0,i,j]*uhf_h2o.mo_coeff[0,i,k]*(mol_h2o.intor('int1e_kin')[j,k]+mol_h2o.intor('int1e_nuc')[j,k])
                sum = sum+ uhf_h2o.mo_coeff[1,j,i]*uhf_h2o.mo_coeff[1,k,i]*(mol_h2o.intor('int1e_kin')[j,k]+mol_h2o.intor('int1e_nuc')[j,k])
```
yes
converged SCF energy = -1.01844561914196  <S^2> = 1.5343282e-13  2S+1 = 1
(10, 10)
(10, 10)
(10, 10)
(10, 10, 10, 10)
(2, 10, 10)
(2, 10)
(100, 100)
'''`

I think I was able to properly get the 1 electron integrals. I am confused why er1.shape is returning (100,100) and how I can get a certain <ij||kl> integral from this 2D nd.array? 

Thanks in advance.
Elijah123
  • 1
  • 1
    Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. – itprorh66 Nov 16 '22 at 00:47

1 Answers1

0

PySCF returns the two-body electron repulsion integrals as a two-dimensional array. To convert it to a four-dimensional array indexed as <ij||kl> you need to reshape the returned array. In your code this corresponds to eri1 = eri1.reshape(10, 10, 10, 10).

Schoyen
  • 7
  • 3