<jupyter_start><jupyter_text>The TB Hamiltonian represents the Fock operator (kinetic + elec-ion + Hartree + exchange (or exchange-corrrelation in RKS)), while the int_1e represents only the kinetic + elec-ion. You may check that the int_1e in the AO basis does not depend on what method you use.<jupyter_code># Obtain the 1e TB model Hamiltonian in LO
ham_tb_lo = np.einsum('i,ai,bi->ab', mf.mo_energy, R_mo_lo, R_mo_lo)
print('TB Ham:\n', ham_tb_lo)

# check that the TB model gives the same mo energies
eigvals, eigvecs = np.linalg.eigh(ham_tb_lo)
print('TB eigenvals:', eigvals)
print('HF mo energies:', mf.mo_energy)
# compute the 1e integrals
ham_1e_kin = pyscf.gto.getints("int1e_kin_sph", mol._atm, mol._bas, mol._env)
ham_1e_nuc = pyscf.gto.getints("int1e_nuc_sph", mol._atm, mol._bas, mol._env)

ham_1e = ham_1e_kin + ham_1e_nuc

# transform into MO basis
ham_1e_mo = np.einsum('kl,ki,lj', ham_1e, mf.mo_coeff, mf.mo_coeff)
print('direct 1e integrals in MO:\n', ham_1e_mo)

print('direct 1e integrals in AO basis:\n', ham_1e)
# transform into the localized basis
ham_1e_lo = np.einsum('kl,ki,lj', ham_1e, a, a)
print('direct 1e integrals in LO (orthogonalized AO):\n', ham_1e_lo)
# compute the 2e integrals
int_2e_ao = pyscf.gto.getints("int2e_sph", mol._atm, mol._bas, mol._env)
int_2e_mo = np.einsum('klmn,ka,lb,mc,nd->abcd', int_2e_ao, mf.mo_coeff, mf.mo_coeff, mf.mo_coeff, mf.mo_coeff)

print('direct 2e integrals in MO:\n',int_2e_mo)
# The density-density interactions betweem MO 0 and 0, 0 and 1 are close in value
print(int_2e_mo[0,0,0,0], int_2e_mo[1,1,1,1])
# The exchange interactions
print(int_2e_mo[0,1,1,0], int_2e_mo[1,0,0,1])
# transform to the localized orbital basis
int_2e_lo = np.einsum('klmn,ka,lb,mc,nd->abcd', int_2e_ao, a, a, a, a)

print('direct 2e integrals in LO:\n',int_2e_lo)<jupyter_output>direct 2e integrals in LO:
 [[[[ 0.70634442 -0.00491842]
   [-0.00491842  0.42596012]]

  [[-0.00491842  0.00987475]
   [ 0.00987475 -0.00491842]]]


 [[[-0.00491842  0.00987475]
   [ 0.00987475 -0.00491842]]

  [[ 0.42596012 -0.00491842]
   [-0.00491842  0.70634442]]]]<jupyter_text>Next, we compare the 1b and 2b energy contributions from the ab-initio (RKS) with those obtained using exact diagonalization of the model Hamiltonian* RKS: 1b and 2b contributions<jupyter_code>rdm1_lo = np.einsum('ij,ai,bj->ab',mf.make_rdm1(), R_ao_lo, R_ao_lo)
e_1b = np.einsum('ij,ij->',ham_1e_lo, rdm1_lo)

rdm2_mo = np.zeros([2,2,2,2])
rdm2_mo[0,0,0,0] = 2
rdm2_lo = np.einsum('ijkl,ai,bj,ck,dl->abcd',rdm2_mo, R_mo_lo, R_mo_lo, R_mo_lo, R_mo_lo)
e_2b = 0.5*np.einsum('ijkl,ijkl->',int_2e_lo, rdm2_lo)

print('RKS (LDA) energy contribution from 1b and 2b:', e_1b, e_2b)
print('RKS total energy - E_nuc:', e_1b+e_2b)<jupyter_output>RKS (LDA) energy contribution from 1b and 2b: -2.3710421147641005 0.5661901887808973
RKS total energy - E_nuc: -1.8048519259832032<jupyter_text>* Exact diagonalization:<jupyter_code>from pyscf import fci

def solver(n,h1,h2):
  nelec=(int(n/2),int(n/2))
  e, fcivec = fci.direct_spin1.kernel(h1, h2, n, nelec, verbose=5)
  dm1, dm2 = fci.direct_spin1.make_rdm12(fcivec, n, nelec)
  return {'e':e, 'dm1':dm1, 'dm2':dm2}

def get_energy_components(n,h1,h2):
  data = solver(n,h1,h2)
  int_expectation = 0.5*np.einsum('ijkl,ijkl->', data['dm2'], h2)
  oneb_expectation = np.einsum('ij,ij->', data['dm1'], h1)

  return {'total_energy_model':data['e'],  
          '1b_expectation':oneb_expectation,
          '2b_expectation':int_expectation} 
# Note that this total energy does not have the contribution from the nuclei interaction
get_energy_components(2,ham_1e_lo,int_2e_lo)
# dmd model from optimized vmc
h1 = np.zeros([2,2])
h1[0,1] = -0.2894
h1[1,0] = h1[0,1]
h2 = np.zeros([2,2,2,2])
h2[0,0,0,0] = 0.1985
h2[1,1,1,1] = 0.1985
e0 = -0.6661

d = get_energy_components(2,h1,h2)
print('The ground state energy from ED of the DMD model:', d['total_energy_model'] + e0)

# This is to be compared with the total ground state energy
# Note that it is doing much better than the RKS<jupyter_output>The ground state energy from ED of the DMD model: -1.154097820345039
{'fraction_non_alphanumeric': 0.1440389294403893, 'fraction_numerical': 0.09416058394160584, 'mean_word_length': 5.334360554699538, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '163424', 'license': 'CC0-1.0', 'n_tokens_mistral': 2024, 'n_tokens_neox': 1833, 'n_words': 409, 'path': '4_15_bare_interaction/h2.ipynb', 'repo_name': 'WagnerGroup/downfolding_school', 'size': '4110'}