<jupyter_start><jupyter_text>MP2  Some useful resources: - [original paper](https://journals.aps.org/pr/abstract/10.1103/PhysRev.46.618) - Levine Chapter 16 - [psi4numpy tutorial](https://github.com/psi4/psi4numpy/blob/master/Tutorials/05_Moller-Plesset/5a_conventional-mp2.ipynb) - [Crawdad programming notes](http://sirius.chem.vt.edu/wiki/doku.php?id=crawdad:programming:project4)  MP2 algorithm1. The starting point will be the Hartree-Fock wavefunction  Imports<jupyter_code>import numpy as np
import scipy.linalg as spla
import psi4
import matplotlib.pyplot as plt
import time
%matplotlib notebook<jupyter_output><empty_output><jupyter_text>Specify the molecule<jupyter_code># start timer
start_time = time.time()
# define molecule
mol = psi4.geometry("""
O 0.0000000 0.0000000 0.0000000
H 0.7569685 0.0000000 -0.5858752
H -0.7569685 0.0000000 -0.5858752
symmetry c1
""")
psi4.set_options({'basis': 'sto-3g'})
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option('BASIS'))
mints = psi4.core.MintsHelper(wfn.basisset())
# get number of electrons
num_elec_alpha =  wfn.nalpha()
num_elec_beta =  wfn.nbeta()
num_elec = num_elec_alpha + num_elec_beta
# get nuclear repulsion energy
E_nuc =  mol.nuclear_repulsion_energy()<jupyter_output><empty_output><jupyter_text>STEP 2 : Calculate molecular integrals Overlap $$ S_{\mu\nu} = (\mu|\nu) = \int dr \phi^*_{\mu}(r) \phi_{\nu}(r) $$Kinetic$$ T_{\mu\nu} = (\mu\left|-\frac{\nabla}{2}\right|\nu) = \int dr \phi^*_{\mu}(r) \left(-\frac{\nabla}{2}\right) \phi_{\nu}(r) $$Nuclear Attraction$$ V_{\mu\nu} = (\mu|r^{-1}|\nu) = \int dr \phi^*_{\mu}(r) r^{-1} \phi_{\nu}(r) $$Form Core Hamiltonian$$ H = T + V $$Two electron integrals$$ (\mu\nu|\lambda\sigma) = \int dr_1 dr_2 \phi^*_{\mu}(r_1) \phi_{\nu}(r_1) r_{12}^{-1} \phi_{\lambda}(r_2) \phi_{\sigma}(r_2) $$<jupyter_code># calculate overlap integrals
S = np.asarray(mints.ao_overlap())
# calculate kinetic energy integrals
T = np.asarray(mints.ao_kinetic())
# calculate nuclear attraction integrals
V = np.asarray(mints.ao_potential())
# form core Hamiltonian
H = T + V
# calculate two electron integrals
eri = np.asarray(mints.ao_eri())
# get number of atomic orbitals
num_ao = np.shape(S)[0]



print(np.shape(eri))
# set inital density matrix to zero
D = np.zeros((num_ao,num_ao))

# 2 helper functions for printing during SCF
def print_start_iterations():
    print("{:^79}".format("{:>4}  {:>11}  {:>11}  {:>11}  {:>11}".format("Iter", "Time(s)", "RMSC DM", "delta E", "E_scf_elec")))
    print("{:^79}".format("{:>4}  {:>11}  {:>11}  {:>11}  {:>11}".format("****", "*******", "*******", "*******", "******")))
def print_iteration(iteration_num, iteration_start_time, iteration_end_time, iteration_rmsc_dm, iteration_E_diff, E_scf_elec):
    print("{:^79}".format("{:>4d}  {:>11f}  {:>.5E}  {:>.5E}  {:>11f}".format(iteration_num, iteration_end_time - iteration_start_time, iteration_rmsc_dm, iteration_E_diff, E_scf_elec)))

# set stopping criteria
iteration_max = 100
convergence_E = 1e-9
convergence_DM = 1e-5
# loop variables
iteration_num = 0
E_scf_total = 0
E_scf_elec = 0.0
iteration_E_diff = 0.0
iteration_rmsc_dm = 0.0
converged = False
exceeded_iterations = False
print_start_iterations()
while (not converged and not exceeded_iterations):
    # store last iteration and increment counters
    iteration_start_time = time.time()
    iteration_num += 1
    E_elec_last = E_scf_elec
    D_last = np.copy(D)
    # form G matrix
    G = np.zeros((num_ao,num_ao))
    for i in range(num_ao):
        for j in range(num_ao):
            for k in range(num_ao):
                for l in range(num_ao):
                    G[i,j] += D[k,l] * ((2.0*(eri[i,j,k,l])) - (eri[i,k,j,l]))
    # build fock matrix
    F  = H + G
    # solve the generalized eigenvalue problem
    E_orbitals, C = spla.eigh(F,S)
    # compute new density matrix
    D = np.zeros((num_ao,num_ao))
    for i in range(num_ao):
        for j in range(num_ao):
            for k in range(num_elec_alpha):
                D[i,j] +=  C[i,k] * C[j,k]
    # calculate electronic energy
    E_scf_elec = np.sum(np.multiply(D , (H +  F)))
    # calculate energy change of iteration
    iteration_E_diff = np.abs(E_scf_elec - E_elec_last)
    # rms change of density matrix
    iteration_rmsc_dm = np.sqrt(np.sum((D - D_last)**2))
    iteration_end_time = time.time()
    print_iteration(iteration_num, iteration_start_time, iteration_end_time, iteration_rmsc_dm, iteration_E_diff, E_scf_elec)
    if(np.abs(iteration_E_diff) < convergence_E and iteration_rmsc_dm < convergence_DM): 
        converged = True
    if(iteration_num == iteration_max):
        exceeded_iterations = True

# calculate total energy
E_scf_total = E_scf_elec + E_nuc
print("{:^79}".format("Total HF energy : {:>11f}".format(E_scf_total)))<jupyter_output><empty_output><jupyter_text>Perform MP2 calculation  Convert the two-electron integrals from AO basis to the MO basis$$(pq|rs) = \sum_\mu \sum_\nu \sum_\lambda \sum_\sigma C_\mu^p C_\nu^q(\mu \nu|\lambda \sigma) C_\lambda^r C_\sigma^s.$$Attempt to code this conversion below, remember that the electron repulsion integrals above are stored as vector `eri` that is of the shape (num_ao,num_ao,num_ao,num_ao). Here the num_ao's for sto-3g water is 7. The resulting tensor will have the same shape as `eri`.<jupyter_code>## place code for two-electron integral conversion here.<jupyter_output><empty_output><jupyter_text>Compute the MP2 EnergyNow we can calculate the MP2 estimation of the correlation energy. $$E_{\mathrm{corr(MP2)}}\ =\ \frac{( ia \mid jb ) [ 2 (ia \mid jb ) - ( ib \mid ja )]}{\epsilon_i + \epsilon_j + \epsilon_a - \epsilon_b}$$Here $i$ and $j$ represent all occupied orbitals, where as $a$ and $b$ will be unoccupied orbitals. Remember during this coding step that we are basing our MP2 correction on an RHF calculation and thus there are the same amount of $\alpha$ and $\beta$ electrons.<jupyter_code>#initialize the variable forthe mp2 correlation energy
E_corr_mp2 = 0
# code the equation above and adjust the value of E_corr_mp2


#this will print your E_corr mp2
print("{:^79}".format("Total MP2 correlation energy : {:>11f}".format(E_corr_mp2)))<jupyter_output><empty_output><jupyter_text>The correlation energy should be very small compared to the total energy (-0.035493 Ha), which is generally the case. However, this correlation energy can be very important to describing properties such as dispersion.   A comparison with Psi4<jupyter_code># Get the SCF wavefunction & energies# Get t 
scf_e, scf_wfn = psi4.energy('scf', return_wfn=True)
mp2_e = psi4.energy('mp2')
print(mp2_e)

E_diff = (mp2_e - (E_total + E_corr_mp2)) 
print(E_diff)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14547622596512147, 'fraction_numerical': 0.031003130123714414, 'mean_word_length': 4.601001669449082, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 42, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 5, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '42825', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2631, 'n_tokens_neox': 2463, 'n_words': 610, 'path': '03_MP2/basics/mp2_psi4.ipynb', 'repo_name': 'kcg12/QMMM_study_group', 'size': '6709'}