Water molecule

This notebook demonstrates how to combine pyscf with the DMRG implementation in chemtensor for approximating the ground state of the water molecule.

import numpy as np
import chemtensor

# PySCF (https://pyscf.org/) defines a molecular orbital basis,
# computes overlap integrals and runs other computational methods for comparison.
import pyscf

# maximum number of OpenMP threads (0 indicates that OpenMP is not available)
chemtensor.get_max_openmp_threads()

16

Define the molecule and perform reference calculations

h2o_atoms = [
    ["O", ( 0.,   0.,   0.)],
    ["H", ( 0.75, 0.47, 0.)],
    ["H", (-0.75, 0.47, 0.)],
]
mol = pyscf.M(atom=h2o_atoms, basis="sto-3g")

# run Hartree-Fock
hf = mol.HF().run()

converged SCF energy = -74.9307084820999

# run coupled-cluster with single and double excitations (CCSD), for comparison
ccsd = pyscf.cc.CCSD(hf).run()
ccsd.e_tot

E(CCSD) = -74.97016403895393  E_corr = -0.03945555685402218

<class 'pyscf.cc.ccsd.CCSD'> does not have attributes  converged

-74.97016403895393

# run full configuration interaction (FCI)
fcisolver = pyscf.fci.FCI(hf)
en_fci, _ = fcisolver.kernel()
en_fci

-74.97027268959357

Electron overlap integrals

# overlap integrals in atomic basis
h1_ao  = mol.intor("int1e_kin") + mol.intor("int1e_nuc")
eri_ao = mol.intor("int2e")

# transform to molecular orbital basis
h1_mo = np.einsum("pi,pq,qj->ij", hf.mo_coeff, h1_ao, hf.mo_coeff)
eri_mo = pyscf.ao2mo.kernel(eri_ao, hf.mo_coeff)
print(h1_mo.shape)
print(eri_mo.shape)

(7, 7)
(7, 7, 7, 7)

# convert to physicists' convention
tkin = h1_mo
vint = np.transpose(eri_mo, (0, 2, 1, 3))

Construct Hamiltonian as MPO and run two-site DMRG

hamiltonian = chemtensor.construct_spin_molecular_hamiltonian_mpo(tkin, vint)
# virtual bond dimensions
hamiltonian.bond_dims

[1, 36, 58, 96, 96, 58, 36, 1]

# local physical quantum numbers (number of electrons and spin)
[chemtensor.decode_quantum_number_pair(qnum) for qnum in hamiltonian.qsite]

[(0, 0), (1, -1), (1, 1), (2, 0)]

# overall quantum number sector of quantum state (number of electrons and spin)
q_pnum = 10
q_spin = 0
qnum_sector = chemtensor.encode_quantum_number_pair(q_pnum, q_spin)

# run two-site DMRG
psi, en_sweeps, entropy = chemtensor.dmrg(hamiltonian, num_sweeps=6, maxiter_lanczos=25,
                                          tol_split=1e-9, qnum_sector=qnum_sector)

Evaluate results

# virtual bond dimensions of optimized MPS
psi.bond_dims

[1, 4, 13, 27, 19, 16, 4, 1]

# energy after each DMRG sweep
en_sweeps

[-84.88903443532217,
 -84.88903443537987,
 -84.88903443537988,
 -84.88903443537988,
 -84.8890344353799,
 -84.88903443537978]

# add nuclear repulsion energy
en_dmrg = en_sweeps[-1] + hf.energy_nuc()
en_dmrg

-74.97027263503881

# difference to CCSD energy
ccsd.e_tot - en_dmrg

0.00010859608488544836

# difference to FCI energy
en_dmrg - en_fci

5.455476070892473e-08