<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step1:  Forte Tutorial 3.00
Step2:  Forte objects (integrals, rdms, ...)
Here we use the utility function forte.utils.prepare_ints_rdms to generate the MO integrals, the reference energy, and RDMs (which are not used in this tutorial)
Step3:  MP2 energy expression
Please follow the notes. We use spin-integrated expression here.
\begin{equation}
E_\text{MP2} = \frac{1}{4} v_{i^\alpha j^\alpha}^{a^\alpha b^\alpha} t_{a^\alpha b^\alpha}^{i^\alpha j^\alpha}
+ \frac{1}{4} v_{i^\beta j^\beta}^{a^\beta b^\beta} t_{a^\beta b^\beta}^{i^\beta j^\beta}
+ v_{i^\alpha j^\beta}^{a^\alpha b^\beta} t_{a^\alpha b^\beta}^{i^\alpha j^\beta}.
\end{equation}
Since $v$ is Hermitian, we will put both $v$ and $t$ to shape $O \times O \times V \times V$.
Step6:  Build Fock matrix
The general expression of spin-orbital Fock matrix is given by
\begin{equation}
f^p_q = h^p_q + \sum_{i} \langle pi || qi \rangle.
\end{equation}
We only need to build the $OO$, $VV$, and $OV$ blocks of the Fock matrix.
Step9:  Build MP2 amplitudes
The spin-orbital MP2 amplitudes read as
\begin{align}
t^{ij}{ab} = \frac{v^{ij}{ab}}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b}.
\end{align}
For convenience, the energy denominator is denoated as
\begin{equation}
\Delta^{ij}_{ab} \equiv \epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b.
\end{equation}
Step13:  DSRG-PT2 energy expression
In DSRG-PT2, we introduce the flow parameter $s \in [0, \infty)$.
The DSRG-PT2 energy expression is given by
\begin{equation}
E_\text{DSRG-PT2} = \frac{1}{4} \tilde{v}{i^\alpha j^\alpha}^{a^\alpha b^\alpha} t{a^\alpha b^\alpha}^{i^\alpha j^\alpha}
+ \frac{1}{4} \tilde{v}{i^\beta j^\beta}^{a^\beta b^\beta} t{a^\beta b^\beta}^{i^\beta j^\beta}
+ \tilde{v}{i^\alpha j^\beta}^{a^\alpha b^\beta} t{a^\alpha b^\beta}^{i^\alpha j^\beta}.
\end{equation}
The spin-orbital renormalized 2e-integrals are
\begin{equation}
\tilde{v}{ij}^{ab} = v{ij}^{ab} \left[ 1 + e^{-s (\Delta^{ij}{ab})^2} \right],
\end{equation}
and the spin-orbital amplitudes are expressed as
\begin{equation}
t^{ij}{ab} = v^{ij}{ab} \frac{1 - e^{-s (\Delta^{ij}{ab})^2}}{\Delta^{ij}_{ab}}.
\end{equation}
<ASSISTANT_TASK:>
Python Code:
import psi4
import forte
import forte.utils
# water molecule
geom = 0 1
O
H 1 1.2
H 1 1.2 2 120.0
basis = '6-31g'
Escf, wfn = forte.utils.psi4_scf(geom, basis, 'rhf')
print(f"RHF energy: {Escf:.8f} Eh")
# Run Psi4 MP2
psi4.set_options({'mp2_type': 'conv',
                  'freeze_core': False})
Emp2_psi4 = psi4.energy('mp2')
print(f"All-electron MP2 energy from Psi4: {Emp2_psi4:.12f} Eh")
psi4.set_options({'mp2_type': 'conv',
                  'freeze_core': True})
Emp2_psi4 = psi4.energy('mp2')
print(f"Frozen-core MP2 energy from Psi4:  {Emp2_psi4:.12f} Eh")
Explanation: Forte Tutorial 3.00: A python implementation of the single-reference DSRG-PT2
In this tutorial we are going to explore how to generate orbitals and visualize them in Jupyter notebooks using the Python API.
End of explanation
# occupations for H2O we will use for MP2 or DSRG-PT2
mo_spaces = {'FROZEN_DOCC': [1,0,0,0],
             'RESTRICTED_DOCC': [2,0,1,1],
             'ACTIVE': [0,0,0,0]}
objs = forte.utils.prepare_ints_rdms(wfn, mo_spaces)
mo_space_info = objs['mo_space_info']
ints = objs['ints']
Eref = objs['reference_energy']
Explanation: Forte objects (integrals, rdms, ...)
Here we use the utility function forte.utils.prepare_ints_rdms to generate the MO integrals, the reference energy, and RDMs (which are not used in this tutorial)
End of explanation
import numpy as np
Explanation: MP2 energy expression
Please follow the notes. We use spin-integrated expression here.
\begin{equation}
E_\text{MP2} = \frac{1}{4} v_{i^\alpha j^\alpha}^{a^\alpha b^\alpha} t_{a^\alpha b^\alpha}^{i^\alpha j^\alpha}
+ \frac{1}{4} v_{i^\beta j^\beta}^{a^\beta b^\beta} t_{a^\beta b^\beta}^{i^\beta j^\beta}
+ v_{i^\alpha j^\beta}^{a^\alpha b^\beta} t_{a^\alpha b^\beta}^{i^\alpha j^\beta}.
\end{equation}
Since $v$ is Hermitian, we will put both $v$ and $t$ to shape $O \times O \times V \times V$.
End of explanation
# figure out occupied and virtual MO indices
core_mos = mo_space_info.corr_absolute_mo("RESTRICTED_DOCC")
ncore = len(core_mos)
print(core_mos)
virt_mos = mo_space_info.corr_absolute_mo("RESTRICTED_UOCC")
nvirt = len(virt_mos)
print(virt_mos)
# build alpha Fock
def build_fock_a(core_mos, virt_mos):
    
    Build alpha Fock matrix.
    :param core_mos: a list of core MO indices
    :param virt_mos: a list of virtual MO indices
    :return: the alpha Fock matrix in format {'oo': mat0, 'ov': mat1, 'vv': mat2}
    
    ncore = len(core_mos)
    nvirt = len(virt_mos)
    
    Fa = {'oo': ints.oei_a_block(core_mos, core_mos),
          'vv': ints.oei_a_block(virt_mos, virt_mos),
          'ov': ints.oei_a_block(core_mos, virt_mos)}
    
    # OO block
    v = ints.tei_aa_block(core_mos, core_mos, core_mos, core_mos)
    Fa['oo'] += np.einsum('piqi->pq', v)
    
    v = ints.tei_ab_block(core_mos, core_mos, core_mos, core_mos)
    Fa['oo'] += np.einsum('piqi->pq', v)
    
    # VV block
    v = ints.tei_aa_block(virt_mos, core_mos, virt_mos, core_mos)
    Fa['vv'] += np.einsum('piqi->pq', v)
    
    v = ints.tei_ab_block(virt_mos, core_mos, virt_mos, core_mos)
    Fa['vv'] += np.einsum('piqi->pq', v)
    
    # OV block
    v = ints.tei_aa_block(core_mos, core_mos, virt_mos, core_mos)
    Fa['ov'] += np.einsum('piqi->pq', v)
    
    v = ints.tei_ab_block(core_mos, core_mos, virt_mos, core_mos)
    Fa['ov'] += np.einsum('piqi->pq', v)
    return Fa
Fa = build_fock_a(core_mos, virt_mos)
Fa
# build beta Fock
def build_fock_b(core_mos, virt_mos):
    
    Build beta Fock matrix.
    :param core_mos: a list of core MO indices
    :param virt_mos: a list of virtual MO indices
    :return: the beta Fock matrix in format {'oo': mat0, 'ov': mat1, 'vv': mat2}
    
    ncore = len(core_mos)
    nvirt = len(virt_mos)
    
    Fb = {'oo': ints.oei_b_block(core_mos, core_mos),
          'vv': ints.oei_b_block(virt_mos, virt_mos),
          'ov': ints.oei_b_block(core_mos, virt_mos)}
    
    # OO block
    v = ints.tei_bb_block(core_mos, core_mos, core_mos, core_mos)
    Fb['oo'] += np.einsum('piqi->pq', v)
    
    v = ints.tei_ab_block(core_mos, core_mos, core_mos, core_mos)
    Fb['oo'] += np.einsum('ipiq->pq', v)
    
    # VV block
    v = ints.tei_bb_block(virt_mos, core_mos, virt_mos, core_mos)
    Fb['vv'] += np.einsum('piqi->pq', v)
    
    v = ints.tei_ab_block(core_mos, virt_mos, core_mos, virt_mos)
    Fb['vv'] += np.einsum('ipiq->pq', v)
    
    # OV block
    v = ints.tei_bb_block(core_mos, core_mos, virt_mos, core_mos)
    Fb['ov'] += np.einsum('piqi->pq', v)
    
    v = ints.tei_ab_block(core_mos, core_mos, core_mos, virt_mos)
    Fb['ov'] += np.einsum('ipiq->pq', v)
    return Fb
Fb = build_fock_b(core_mos, virt_mos)
Fb
# test if alpha and beta Fock matrices are indentical
assert np.all(np.abs(Fa['oo'] - Fb['oo']) < 1.0e-10)
assert np.all(np.abs(Fa['ov'] - Fb['ov']) < 1.0e-10)
assert np.all(np.abs(Fa['vv'] - Fb['vv']) < 1.0e-10)
Explanation: Build Fock matrix
The general expression of spin-orbital Fock matrix is given by
\begin{equation}
f^p_q = h^p_q + \sum_{i} \langle pi || qi \rangle.
\end{equation}
We only need to build the $OO$, $VV$, and $OV$ blocks of the Fock matrix.
End of explanation
# build orbital energies
Fdiag = {'o': np.diag(Fa['oo']), 'v': np.diag(Fa['vv'])}
Fdiag
# grab 2e-integrals
V = {'aa': ints.tei_aa_block, 'ab': ints.tei_ab_block, 'bb': ints.tei_bb_block}
Vaa = V['aa'](core_mos, core_mos, virt_mos, virt_mos)
Vab = V['ab'](core_mos, core_mos, virt_mos, virt_mos)
Vbb = V['bb'](core_mos, core_mos, virt_mos, virt_mos)
# build MP2 amplitudes
def build_amplitudes_mp2(Fdiag, V_spin):
    
    Build MP2 amplitudes for spin cases corresponding to V.
    :param Fdiag: a map for orbital energies {'o': occupied epsilon, 'v': virtual epsilon}
    :param V: spin cases for antisymmetrized 2e-integrals
    :return: the MP2 amplitudes
    
    T = V[V_spin](core_mos, core_mos, virt_mos, virt_mos).copy()
    for ii, i in enumerate(Fdiag['o']):
        for ij, j in enumerate(Fdiag['o']):
            for ia, a in enumerate(Fdiag['v']):
                for ib, b in enumerate(Fdiag['v']):
                    T[ii][ij][ia][ib] /= i + j - a - b
    return T
Taa = build_amplitudes_mp2(Fdiag, 'aa')
Tab = build_amplitudes_mp2(Fdiag, 'ab')
Tbb = build_amplitudes_mp2(Fdiag, 'bb')
# contract 2e-integrals with 2-amplitudes
def contraction(v, t, factor):
    
    Fully contract v and t.
    :param v: 2e-integrals
    :param t: 2-amplitudes
    :param factor: scaling factor
    :return: the contracted value
    
    return factor * np.einsum('ijab,ijab->', v, t)
# compute MP2 energy
Emp2_corr = contraction(Vaa, Taa, 0.25) + contraction(Vbb, Tbb, 0.25) \
          + contraction(Vab, Tab, 1.0)
Emp2 = Eref + Emp2_corr
print(f"MP2 correlation energy:     {Emp2_corr:18.12f}")
print(f"MP2 total energy tutorial:  {Emp2:18.12f}")
print(f"MP2 total energy from Psi4: {Emp2_psi4:18.12f}")
assert(abs(Emp2 - Emp2_psi4) < 1.0e-8)
Explanation: Build MP2 amplitudes
The spin-orbital MP2 amplitudes read as
\begin{align}
t^{ij}{ab} = \frac{v^{ij}{ab}}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b}.
\end{align}
For convenience, the energy denominator is denoated as
\begin{equation}
\Delta^{ij}_{ab} \equiv \epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b.
\end{equation}
End of explanation
flow_param = 1.0
# renormalize 2e-integrals
def renorm_aptei(Fdiag, V_spin, s):
    
    Renormalize antisymmetrized 2e-integrals for DSRG-PT2.
    :param Fdiag: a map for orbital energies {'o': occupied epsilon, 'v': virtual epsilon}
    :param V_spin: spin case for antisymmetrized 2e-integrals
    :param s: the DSRG flow parameter
    :return: the renormalized 2e-integrals
    
    R = V[V_spin](core_mos, core_mos, virt_mos, virt_mos).copy()
    for ii, i in enumerate(Fdiag['o']):
        for ij, j in enumerate(Fdiag['o']):
            for ia, a in enumerate(Fdiag['v']):
                for ib, b in enumerate(Fdiag['v']):
                    delta = i + j - a - b
                    R[ii][ij][ia][ib] *= 1.0 + np.exp(-s * delta ** 2)
    return R
rVaa = renorm_aptei(Fdiag, 'aa', flow_param)
rVab = renorm_aptei(Fdiag, 'ab', flow_param)
rVbb = renorm_aptei(Fdiag, 'bb', flow_param)
# build DSRG-PT2 amplitudes
def build_amplitudes_pt2(Fdiag, V_spin, s):
    
    Build DSRG-PT2 amplitudes for spin cases corresponding to V.
    :param Fdiag: a map for orbital energies {'o': occupied epsilon, 'v': virtual epsilon}
    :param V_spin: spin case for antisymmetrized 2e-integrals
    :param s: the DSRG flow parameter
    :return: the DSR-PT2 amplitudes
    
    T = V[V_spin](core_mos, core_mos, virt_mos, virt_mos).copy()
    for ii, i in enumerate(Fdiag['o']):
        for ij, j in enumerate(Fdiag['o']):
            for ia, a in enumerate(Fdiag['v']):
                for ib, b in enumerate(Fdiag['v']):
                    delta = i + j - a - b
                    if abs(delta) < 1.0e-8:
                        value = 0.0
                    else:
                        value = (1.0 - np.exp(-s * delta ** 2)) / delta
                    T[ii][ij][ia][ib] *= value
    return T
rTaa = build_amplitudes_pt2(Fdiag, 'aa', flow_param)
rTab = build_amplitudes_pt2(Fdiag, 'ab', flow_param)
rTbb = build_amplitudes_pt2(Fdiag, 'bb', flow_param)
from math import isclose
# reference value (s = 1.0)
Ept2_corr_ref = -0.148255452641
# compute DSRG-PT2 energy
Ept2_corr = contraction(rVaa, rTaa, 0.25) + contraction(rVbb, rTbb, 0.25) \
          + contraction(rVab, rTab, 1.0)
Ept2 = Eref + Ept2_corr
print(f"DSRG-PT2 correlation energy:     {Ept2_corr:18.12f}")
print(f"DSRG-PT2 total energy tutorial:  {Ept2:18.12f}")
print(f"MP2 total energy from Psi4:      {Emp2_psi4:18.12f}")
# test against the reference value
assert isclose(Ept2_corr, Ept2_corr_ref, rel_tol=0, abs_tol=1e-8)
# play around with flow parameter
def compute_dsrgpt2_energy(s):
    
    Compute the DSRG-PT2 energy.
    :param s: the DSRG flow parameter
    :return: the DSRG-PT2 correlation energy
    
    rV = renorm_aptei(Fdiag, 'aa', s)
    rT = build_amplitudes_pt2(Fdiag, 'aa', s)
    Eout = contraction(rV, rT, 0.25)
    
    rV = renorm_aptei(Fdiag, 'bb', s)
    rT = build_amplitudes_pt2(Fdiag, 'bb', s)
    Eout += contraction(rV, rT, 0.25)
    
    rV = renorm_aptei(Fdiag, 'ab', s)
    rT = build_amplitudes_pt2(Fdiag, 'ab', s)
    Eout += contraction(rV, rT, 1.0)
    
    return Eout
dsrg_s = [0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 5.0, 7.5, 10.0, 25.0, 50.0, 100.0, 1000.0]
energies = [compute_dsrgpt2_energy(s) for s in dsrg_s]
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib notebook
mpl.rcParams['pdf.fonttype'] = 42
fig = plt.figure(figsize=(6,4), dpi=150)
ax = plt.subplot(111)
ax.semilogx(dsrg_s, energies)
ax.axhline(Emp2_corr, lw=1.0, ls=':', color='k')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
Explanation: DSRG-PT2 energy expression
In DSRG-PT2, we introduce the flow parameter $s \in [0, \infty)$.
The DSRG-PT2 energy expression is given by
\begin{equation}
E_\text{DSRG-PT2} = \frac{1}{4} \tilde{v}{i^\alpha j^\alpha}^{a^\alpha b^\alpha} t{a^\alpha b^\alpha}^{i^\alpha j^\alpha}
+ \frac{1}{4} \tilde{v}{i^\beta j^\beta}^{a^\beta b^\beta} t{a^\beta b^\beta}^{i^\beta j^\beta}
+ \tilde{v}{i^\alpha j^\beta}^{a^\alpha b^\beta} t{a^\alpha b^\beta}^{i^\alpha j^\beta}.
\end{equation}
The spin-orbital renormalized 2e-integrals are
\begin{equation}
\tilde{v}{ij}^{ab} = v{ij}^{ab} \left[ 1 + e^{-s (\Delta^{ij}{ab})^2} \right],
\end{equation}
and the spin-orbital amplitudes are expressed as
\begin{equation}
t^{ij}{ab} = v^{ij}{ab} \frac{1 - e^{-s (\Delta^{ij}{ab})^2}}{\Delta^{ij}_{ab}}.
\end{equation}
End of explanation
<END_TASK>