<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 pyscf
from pyscf import gto, scf
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 = pyscf.gto.M(
    atom="O 0.0000000 0.0000000 0.0000000; H 0.7569685 0.0000000 -0.5858752; H -0.7569685 0.0000000 -0.5858752",
    basis='sto-3g',
    unit="Ang",
    verbose=0,
    symmetry=False,
    spin=0,
    charge=0
)
# get number of atomic orbitals
num_ao = mol.nao_nr()
# get number of electrons
num_elec_alpha, num_elec_beta = mol.nelec
num_elec = num_elec_alpha + num_elec_beta
# get nuclear repulsion energy
E_nuc = mol.energy_nuc()<jupyter_output><empty_output><jupyter_text>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 = mol.intor('cint1e_ovlp_sph')
# calculate kinetic energy integrals
T = mol.intor('cint1e_kin_sph')
# calculate nuclear attraction integrals
V = mol.intor('cint1e_nuc_sph')
# form core Hamiltonian
H = T + V
# calculate two electron integrals
eri = mol.intor('cint2e_sph', aosym='s8')
# since we are using the 8 fold symmetry of the 2 electron integrals
# the functions below will help us when accessing elements
__idx2_cache = {}


def idx2(i, j):
    if (i, j) in __idx2_cache:
        return __idx2_cache[i, j]
    elif i >= j:
        __idx2_cache[i, j] = int(i*(i+1)/2+j)
    else:
        __idx2_cache[i, j] = int(j*(j+1)/2+i)
    return __idx2_cache[i, j]


def idx4(i, j, k, l):
    return idx2(idx2(i, j), idx2(k, l))


print(np.shape(eri))<jupyter_output>(406,)<jupyter_text>Perform Hartree-Fock SCF<jupyter_code># 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_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_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_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[idx4(i, j, k, l)])) -
                         (eri[idx4(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>Iter      Time(s)      RMSC DM      delta E       E_elec            
           ****      *******      *******      *******       ******            
              1     0.012981  2.69561E+00  1.27367E+02  -127.366748            
              2     0.011632  1.84626E+00  4.69671E+01   -80.399634            
              3     0.012944  1.84892E-01  4.07021E+00   -84.469843            
              4     0.012603  3.65179E-02  3.36584E-01   -84.133260            
              5     0.016098  1.41819E-02  2.77766E-02   -84.161036            
              6     0.011819  5.65413E-03  2.81184E-03   -84.158224            
              7     0.011283  2.37192E-03  1.60719E-04   -84.158064            
              8     0.015524  1.00945E-03  1.58481E-04   -84.157905            
              9     0.011104  4.33408E-04  6.31723E-05   -84.157842            
             10     0.014282  1.86875E-04  2.79308E-05   -84.157814            
             11     0.013606  8.07615E-0[...]<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.$$This is implemented in the cell block below. There are a few ways to implement this, below is by far the worst. The algorithm coded below is the naive approach known as the Noddy algorithm. This algorithm scales as $N^8$, although MP2 is formally known to scale as $N^5$; however. The noddy algorithm is a great starting point.<jupyter_code>eri_mo = np.zeros((num_ao, num_ao, num_ao, num_ao))
for p in range(num_ao):
    for q in range(num_ao):
        for r in range(num_ao):
            for s in range(num_ao):
                for mu in range(num_ao):
                    for nu in range(num_ao):
                        for lmda in range(num_ao):
                            for sigma in range(num_ao):
                                eri_mo[p, q, r, s] += C[mu, p]*C[nu, q]*C[lmda,
                                                                          r]*C[sigma, s]*eri[idx4(mu, nu, lmda, sigma)]<jupyter_output><empty_output><jupyter_text>Compute the MP2 Energy<jupyter_code>E_corr_mp2 = 0
for i in range(num_elec_alpha):
    for j in range(num_elec_alpha):
        for a in range(num_elec_alpha, num_ao):
            for b in range(num_elec_alpha, num_ao):
                temp = eri_mo[i, a, j, b] * \
                    (2*eri_mo[i, a, j, b] - eri_mo[i, b, j, a])
                temp /= (E_orbitals[i] + E_orbitals[j] - E_orbitals[a] - E_orbitals[b])
                E_corr_mp2 += temp
print("{:^79}".format("Total MP2 correlation energy : {:>11f}".format(E_corr_mp2)))<jupyter_output><empty_output><jupyter_text>The correlation energy is very small compared to the total energy, which is generally the case. However, this correlation energy can be very important to describing properties such as dispersion.   A comparison with PySCF<jupyter_code>from pyscf import mp
m = scf.RHF(mol)
print('E(HF) = %g' % m.kernel())
mp2 = mp.MP2(m)
E_corr_mp2_pyscf = mp2.kernel()[0]
print('E(MP2) = {:.9g}'.format(E_corr_mp2_pyscf))
#comparison from pyscf
E_diff = E_corr_mp2_pyscf - E_corr_mp2
print(E_diff)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13903928118880313, 'fraction_numerical': 0.0643935030526437, 'mean_word_length': 2.660202360876897, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 45, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 5, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '109911', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3559, 'n_tokens_neox': 3225, 'n_words': 745, 'path': '03_MP2/basics/mp2_pyscf_solution.ipynb', 'repo_name': 'b-v-henders/QMMM_study_group', 'size': '8681'}