{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.8.8"}},"nbformat_minor":5,"nbformat":4,"cells":[{"cell_type":"code","source":"\"\"\"Møller-Plesset perturbation theory to second order\"\"\"\n\n__author__ = \"Roberto Di Remigio\"\n__credit__ = [\"Roberto Di Remigio\", \"Xin Li\"]\n\n__copyright__ = \"(c) 2021, ENCSS and PDC\"\n__license__ = \"MIT\"\n__date__ = \"2021-04-14\"","metadata":{"tags":["remove_cell"],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n \n
","metadata":{"tags":["remove_cell"]}},{"cell_type":"markdown","source":"# Møller-Plesset perturbation theory to second order","metadata":{"tags":[]}},{"cell_type":"markdown","source":"
\n \n \n
\n To run the selected code cell, hit 
Shift + Enter
\n
\n
","metadata":{"tags":["remove_cell"]}},{"cell_type":"markdown","source":"## Theory refresher\n\n\n### Formal Rayleigh-Schrödinger perturbation theory\n\nWe want to solve the time-independent Schrödinger equation: \n\n$$\nH \\Psi = E \\Psi\n$$\n\nIn perturbation theory, the Hamiltonian is partitioned into *zeroth-order* and *perturbation* terms:\n\n$$\n(H_0 + \\lambda V) \\Psi = E \\Psi\n$$\n\nwith $\\lambda$ the coupling strength of the perturbation. We assume that the complete spectrum of the zeroth-order Hamiltonian is known:\n\n$$\nH_0 \\Psi^{(0)}_{n} = E_{n}^{(0)} \\Psi^{(0)}_{n}\n$$\n\nThe energy and wavefunction are expanded in a formal power series in $\\lambda$ and terms to the same order are collected on the left- and right-hand sides:\n\n$$\n(H_0 + \\lambda V) \\left( \\sum_{k} \\lambda^{k} \\Psi_{n}^{(k)} \\right) = \\left( \\sum_{k} \\lambda^{k} E_{n}^{(k)} \\right) \\left(\\sum_{j} \\lambda^{j} \\Psi_{n}^{(j)}\\right).\n$$\n\nThus:\n\n$$\n(E_{n}^{(0)} - H_{0})\\Psi_{n}^{(m)} = V\\Psi_{n}^{(m-1)} - \\sum_{l = 0}^{m-1} E_{n}^{(m-l)}\\Psi_{n}^{(l)}.\n$$\n\n\nEnergy corrections and perturbative expansion coefficients, $a_{kn}^{(m)} \\equiv \\langle \\Psi_{k}^{(0)} | \\Psi_{n}^{(m)} \\rangle$ , for the wavefunction are obtained by projection onto the set of known zeroth-order eigenfunctions:\n\n$$\nE_{n}^{(m)} = \\left\\langle \\Psi_{n}^{(0)} \\left| V \\right| \\Psi_{n}^{(m-1)} \\right\\rangle,\n$$\n\nand:\n\n$$\n[E_{n}^{(0)} - E_{n}^{(0)}]a_{kn}^{(m)} = \n\\sum_{j} \\langle \\Psi_{k}^{(0)} | V | \\Psi_{j}^{(0)} \\rangle a_{jn}^{(m-1)} - \n\\sum_{l = 0}^{m-1} E_{n}^{(m-l)} a_{kn}^{(l)},\\quad \na_{kn}^{(0)} = \\delta_{kn},\\,\\,a_{nn}^{(m)} = \\delta_{m0}\n$$","metadata":{}},{"cell_type":"markdown","source":"### Møller-Plesset partitioning and second order energy\n\nIn general, the partitioning of the Hamiltonian can be achieved in a number of ways. In molecular electronic structure theory, $H$ is the Born-Oppenheimer, many-body, molecular electronic Hamiltonian. The Møller-Plesset partitioning is probably the most common, since it follows quite naturally when first approximating the desired ground state in terms of a *single reference determinant*. The Hamiltonian is in fact rewritten as:\n\n$$\nH = F + \\Phi\n$$\n\nwhere $F$ is the Fock operator and $\\Phi$ is the so-called *fluctuation potential*. $F$ is a one-body operator whose spectrum consists of all determinants that can be built from excitation of the reference $| 0 \\rangle$:\n\n$$\nF | 0 \\rangle = E_{0}^{(0)}| 0 \\rangle, \\quad\nF\\left|_{ij\\ldots} ^{ab\\ldots} \\right\\rangle =\n\\left(E_{0}^{(0)} + \\varepsilon^{ab\\cdots}_{ij\\cdots} \\right)\n\\left|_{ij\\ldots} ^{ab\\ldots} \\right\\rangle,\n$$\n\nwhere the *zeroth-order energy* and the *orbital energy denominators* are:\n\n$$\nE_{0}^{(0)} = \\sum_{i} \\varepsilon_{i},\\quad\n\\varepsilon^{ab\\cdots}_{ij\\cdots} = \\varepsilon_{a} + \\varepsilon_{b} + \\cdots - \\varepsilon_{i} - \\varepsilon_{j} - \\cdots\n$$\n\nThe fluctuation potential is a two-body operator. With this partinioning, the zeroth-order energy is the sum of orbital energies. The first-order correction is:\n\n$$\nE_0^{(1)} = \\left\\langle 0 \\left| \\Phi \\right| 0 \\right\\rangle = -\\frac{1}{2} \\sum_{ij} \\langle ij \\| ij \\rangle\n$$\n\nthat is, the energy of the reference single determinant is *correct* throught first order in the perturbative series: $E_{\\mathrm{ref}} = E_{0}^{(0)} + E_0^{(1)}$.\n\nThe first-order wavefunction is obtained from the general RSPT expression. In a basis of molecular spin-orbitals:\n\n$$\n| \\Psi^{(1)} \\rangle = -\\frac{1}{4}\\sum_{ij}^{N_{\\mathrm{O}}} \\sum_{ab}^{N_{\\mathrm{V}}} \n\\frac{\\langle ab \\| ij \\rangle}{\\varepsilon_{ij}^{ab}} |_{ij}^{ab}\\rangle,\n$$\n\nwhere the *orbital energy denominator* is: $\\varepsilon_{ij}^{ab} = \\varepsilon_{i} + \\varepsilon_{j} -\\varepsilon_{a} - \\varepsilon_{b} $. The second order energy correction follows:\n\n$$\nE_{0}^{(2)} \\equiv \nE_{\\mathrm{MP2}} = -\n\\frac{1}{4} \n\\sum_{ij}^{N_{\\mathrm{O}}} \\sum_{ab}^{N_{\\mathrm{V}}} \n\\frac{\\langle ij \\| ab \\rangle \\langle ab \\| ij \\rangle}{\\varepsilon_{ij}^{ab}}.\n$$\n\nFor a closed-shell, restricted reference using real MOs:\n\n$$\nE_{\\mathrm{MP2}} = -\n\\sum_{ij}^{N_{\\mathrm{O}}} \\sum_{ab}^{N_{\\mathrm{V}}} \n\\frac{\\langle ij | ab \\rangle}{\\varepsilon_{ij}^{ab}}\n[ 2 \\langle ij | ab \\rangle - \\langle ij | ba \\rangle ],\n$$\n\nwhich we can further rearrange into to two terms, *opposite-spin* and *same-spin*:\n\n$$\nE_{\\mathrm{MP2}} = -\n\\sum_{ij}^{N_{\\mathrm{O}}} \\sum_{ab}^{N_{\\mathrm{V}}} \n\\frac{\\langle ij | ab \\rangle\\langle ij | ab \\rangle}{\\varepsilon_{ij}^{ab}} -\n\\sum_{ij}^{N_{\\mathrm{O}}} \\sum_{ab}^{N_{\\mathrm{V}}} \n\\frac{\\langle ij | ab \\rangle[ \\langle ij | ab \\rangle - \\langle ij | ba \\rangle ]}{\\varepsilon_{ij}^{ab}} = \nE_{\\mathrm{MP2}}^{\\mathrm{OS}} + E_{\\mathrm{MP2}}^{\\mathrm{SS}}.\n$$","metadata":{}},{"cell_type":"markdown","source":"## Implementation\n\nTo compute $E_{\\mathrm{MP2}}$ we need to:\n\n1. Obtain the reference closed-shell determinant from a Hartree-Fock calculation.\n2. Transform the AO basis ERI tensor to MO basis.\n3. Assemble the energy denominators.\n4. Combine the results of steps 2 and 3 to form the perturbative correction.\n\n![Obtaining the MP2 energy correction](../img/mp2.svg)\n\nWe start with the declaration of the usual water molecule and its basis set. We also perform the SCF calculation with the `ScfRestrictedDriver`.","metadata":{}},{"cell_type":"code","source":"import veloxchem as vlx\n\nh2o_xyz = \"\"\"3\nwater \nO 0.000000000000 0.000000000000 0.000000000000 \nH 0.000000000000 0.740848095288 0.582094932012 \nH 0.000000000000 -0.740848095288 0.582094932012\n\"\"\"\n\nmol = vlx.Molecule.from_xyz_string(h2o_xyz)\n\nbasis = vlx.MolecularBasis.read(mol, \"cc-pvdz\")\n\nscfdrv = vlx.ScfRestrictedDriver()\nscfdrv.compute(mol, basis)","metadata":{"tags":["hide_output"]},"execution_count":1,"outputs":[{"name":"stdout","text":"* Warning * Environment variable OMP_NUM_THREADS not set.\n* Warning * Setting OMP_NUM_THREADS to 8.\n \n Self Consistent Field Driver Setup \n ==================================== \n \n Wave Function Model : Spin-Restricted Hartree-Fock \n Initial Guess Model : Superposition of Atomic Densities \n Convergence Accelerator : Two Level Direct Inversion of Iterative Subspace \n Max. Number of Iterations : 50 \n Max. Number of Error Vectors : 10 \n Convergence Threshold : 1.0e-06 \n ERI Screening Scheme : Cauchy Schwarz + Density \n ERI Screening Mode : Dynamic \n ERI Screening Threshold : 1.0e-12 \n Linear Dependence Threshold : 1.0e-06 \n \n* Info * Nuclear repulsion energy: 9.3436381580 au \n \n* Info * Overlap matrix computed in 0.05 sec. \n \n* Info * Kinetic energy matrix computed in 0.03 sec. \n \n* Info * Nuclear potential matrix computed in 0.01 sec. \n \n* Info * Orthogonalization matrix computed in 0.10 sec. \n \n* Info * SAD initial guess computed in 0.06 sec. \n \n* Info * Starting Reduced Basis SCF calculation... \n* Info * ...done. SCF energy in reduced basis set: -75.979046359568 au. Time: 1.66 sec. \n \n* Info * Overlap matrix computed in 0.08 sec. \n \n* Info * Kinetic energy matrix computed in 0.00 sec. \n \n* Info * Nuclear potential matrix computed in 0.01 sec. \n \n* Info * Orthogonalization matrix computed in 0.00 sec. \n \n \n Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change \n -------------------------------------------------------------------------------------------- \n 1 -76.025869744699 0.0000000000 0.09892066 0.01121867 0.00000000 \n 2 -76.026932827150 -0.0010630825 0.01920269 0.00294570 0.02594960 \n 3 -76.026981056596 -0.0000482294 0.00412658 0.00076652 0.00507499 \n 4 -76.026983601528 -0.0000025449 0.00247905 0.00043111 0.00166621 \n 5 -76.026984175646 -0.0000005741 0.00021766 0.00003639 0.00052311 \n 6 -76.026984187023 -0.0000000114 0.00003206 0.00000540 0.00012084 \n 7 -76.026984187254 -0.0000000002 0.00000208 0.00000029 0.00001663 \n 8 -76.026984187255 -0.0000000000 0.00000053 0.00000008 0.00000087 \n \n *** SCF converged in 8 iterations. Time: 4.48 sec. \n \n Spin-Restricted Hartree-Fock: \n ----------------------------- \n Total Energy : -76.0269841873 au \n Electronic Energy : -85.3706223452 au \n Nuclear Repulsion Energy : 9.3436381580 au \n ------------------------------------ \n Gradient Norm : 0.0000005316 au \n \n \n Ground State Information \n ------------------------ \n Charge of Molecule : 0.0 \n Multiplicity (2S+1) : 1.0 \n Magnetic Quantum Number (M_S) : 0.0 \n \n \n Spin Restricted Orbitals \n ------------------------ \n \n Molecular Orbital No. 1: \n -------------------------- \n Occupation: 2.0 Energy: -20.54819 au \n ( 1 O 1s : 1.00) \n \n Molecular Orbital No. 2: \n -------------------------- \n Occupation: 2.0 Energy: -1.34520 au \n ( 1 O 2s : -0.44) ( 1 O 3s : -0.36) ( 2 H 1s : -0.20) \n ( 3 H 1s : -0.20) \n \n Molecular Orbital No. 3: \n -------------------------- \n Occupation: 2.0 Energy: -0.70585 au \n ( 1 O 1p-1: -0.49) ( 1 O 2p-1: -0.22) ( 2 H 1s : -0.33) \n ( 3 H 1s : 0.33) \n \n Molecular Orbital No. 4: \n -------------------------- \n Occupation: 2.0 Energy: -0.57109 au \n ( 1 O 2s : 0.15) ( 1 O 3s : 0.36) ( 1 O 1p0 : -0.55) \n ( 1 O 2p0 : -0.36) ( 2 H 1s : -0.21) ( 3 H 1s : -0.21) \n \n Molecular Orbital No. 5: \n -------------------------- \n Occupation: 2.0 Energy: -0.49457 au \n ( 1 O 1p+1: -0.63) ( 1 O 2p+1: -0.49) \n \n Molecular Orbital No. 6: \n -------------------------- \n Occupation: 0.0 Energy: 0.18787 au \n ( 1 O 3s : 1.02) ( 1 O 1p0 : 0.19) ( 1 O 2p0 : 0.33) \n ( 2 H 2s : -0.83) ( 3 H 2s : -0.83) \n \n Molecular Orbital No. 7: \n -------------------------- \n Occupation: 0.0 Energy: 0.25852 au \n ( 1 O 1p-1: -0.28) ( 1 O 2p-1: -0.67) ( 2 H 2s : 1.48) \n ( 3 H 2s : -1.48) \n \n Molecular Orbital No. 8: \n -------------------------- \n Occupation: 0.0 Energy: 0.79749 au \n ( 1 O 1p-1: 0.26) ( 1 O 2p-1: 0.49) ( 2 H 1s : -0.95) \n ( 2 H 2s : 0.66) ( 2 H 1p0 : -0.16) ( 3 H 1s : 0.95) \n ( 3 H 2s : -0.66) ( 3 H 1p0 : 0.16) \n \n Molecular Orbital No. 9: \n -------------------------- \n Occupation: 0.0 Energy: 0.87271 au \n ( 1 O 2s : 0.25) ( 1 O 3s : -0.34) ( 1 O 1p0 : 0.34) \n ( 2 H 1s : -0.77) ( 2 H 2s : 0.54) ( 2 H 1p-1: -0.31) \n ( 3 H 1s : -0.77) ( 3 H 2s : 0.54) ( 3 H 1p-1: 0.31) \n \n Molecular Orbital No. 10: \n -------------------------- \n Occupation: 0.0 Energy: 1.16315 au \n ( 1 O 3s : -0.78) ( 1 O 1p0 : 0.74) ( 1 O 2p0 : -1.31) \n ( 2 H 1s : 0.59) ( 2 H 1p0 : -0.24) ( 3 H 1s : 0.59) \n ( 3 H 1p0 : -0.24) \n \n","output_type":"stream"}]},{"cell_type":"markdown","source":"We can now access orbital energies and MO coefficients from the driver:","metadata":{}},{"cell_type":"code","source":"epsilon = scfdrv.scf_tensors[\"E\"]\nC = scfdrv.scf_tensors[\"C\"]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## The Integral transformation\n\nWe compute the MP2 energy correction with the ERI expressed in MO basis: we need to transform the ERI tensor from AO basis.\nThe transformation reads:\n\n$$\n\\langle pq | rs \\rangle = \\sum_{\\mu\\nu\\kappa\\lambda} C_{\\mu p}C_{\\nu r} (\\mu\\nu|\\kappa\\lambda) C_{\\kappa q} C_{\\lambda s},\n$$\n\nwith the MO integrals in [**physicists' notation**](http://vergil.chemistry.gatech.edu/notes/permsymm/permsymm.html). The transformation requires $O(N^{8})$ operation count. \nHowever, we can perform it more efficiently as a stepwise contraction:\n\n$$\n\\langle pq | rs \\rangle = \\sum_{\\mu} C_{\\mu p} \\left(\\sum_{\\nu} C_{\\nu r} \\left (\\sum_{\\kappa} \\left(\\sum_{\\lambda} (\\mu\\nu|\\kappa\\lambda) C_{\\lambda s} \\right) C_{\\kappa q} \\right)\\right).\n$$\n\nWe should also note that we do **not** need the full ERI tensor in MO basis, but rather the *OOVV* class of integrals, which involve two occupied and two virtual MO indices:\n\n$$\n\\langle ij | ab \\rangle = \n\\sum_{\\mu} C_{\\mu i} \n\\left(\\sum_{\\nu} C_{\\nu j} \n\\left(\\sum_{\\kappa} \n\\left(\\sum_{\\lambda} (\\mu\\kappa|\\nu\\lambda) C_{\\lambda b} \\right)\nC_{\\kappa a}\\right)\\right).\n$$","metadata":{}},{"cell_type":"code","source":"eridrv = vlx.ElectronRepulsionIntegralsDriver()\nmknl = eridrv.compute_in_mem(mol, basis)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\n\nN_O = mol.number_of_electrons() // 2\nN_V = scfdrv.mol_orbs.number_mos() - N_O\n\n# lambda -> b transformation\nmknb = np.einsum(\"mknl,lB->mknB\", mknl, C[:, N_O:])\nprint(f\"{mknb.shape=}\")\n# TODO kappa -> a transformation\nmnab =\nprint(f\"{mnab.shape=}\")\n# TODO nu -> j transformation\nmjab =\nprint(f\"{mjab.shape=}\")\n# TODO mu -> i transformation\nijab =\nprint(f\"{ijab.shape=}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's compare our *OOVV* ERI tensor with the one computed by using VeloxChem's own `MOIntegralsDriver`:","metadata":{}},{"cell_type":"code","source":"moeridrv = vlx.MOIntegralsDriver()\nmoeri = moeridrv.compute_in_mem(mol, basis, mol_orbs=scfdrv.mol_orbs, mints_type=\"OOVV\")\n\nnp.testing.assert_allclose(ijab, moeri, atol=1.e-10)","metadata":{"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## The MP2 energy correction\n\nWe now have all the ingrediente to compute the *opposite-spin* and *same-spin* components of the MP2 energy correction:\n\n$$\nE_{\\mathrm{MP2}} = -\n\\sum_{ij}^{N_{\\mathrm{O}}} \\sum_{ab}^{N_{\\mathrm{V}}} \n\\frac{\\langle ij | ab \\rangle\\langle ij | ab \\rangle}{\\varepsilon_{ij}^{ab}} - \n\\sum_{ij}^{N_{\\mathrm{O}}} \\sum_{ab}^{N_{\\mathrm{V}}} \n\\frac{\\langle ij | ab \\rangle[ \\langle ij | ab \\rangle - \\langle ij | ba \\rangle ]}{\\varepsilon_{ij}^{ab}} = \nE_{\\mathrm{MP2}}^{\\mathrm{OS}} + E_{\\mathrm{MP2}}^{\\mathrm{SS}}.\n$$","metadata":{}},{"cell_type":"code","source":"e_mp2_ss = 0.0\ne_mp2_os = 0.0\n\n# TODO extract the occupied subset of the orbital energies using NumPy slicing\ne_ij =\n# TODO extract the virtual subset of the orbital energies using NumPy slicing\ne_ab =\n\n# loop over occupied orbitals\nfor i in range(N_O):\n # loop over occupied orbitals\n for j in range(N_O):\n # loop over virtual orbitals\n for a in range(N_V):\n # loop over virtual orbitals\n for b in range(N_V):\n # TODO form enegy denominator for given i, j, a, b quadruplet of indices\n e_ijab =\n \n # TODO update opposite-spin component of the energy\n e_mp2_os -=\n \n # TODO update same-spin component of the energy\n e_mp2_ss -=","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Opposite-spin MP2 energy: {e_mp2_os:20.12f}\")\nprint(f\"Same-spin MP2 energy: {e_mp2_ss:20.12f}\")\nprint(f\"MP2 energy: {e_mp2_os + e_mp2_ss:20.12f}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"VeloxChem has its own implementation of the MP2 energy correction. We can check our result against it.","metadata":{}},{"cell_type":"code","source":"mp2drv = vlx.Mp2Driver()\nmp2drv.compute_conventional(mol, basis, scfdrv.mol_orbs)\n\nnp.testing.assert_allclose(e_mp2_os + e_mp2_ss, mp2drv.e_mp2, atol=1e-9)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## An alternative implmentation using `np.einsum`\n\nWe can compute the MP2 energy terms using `np.einsum` to handle the tensor contractions. First, we need to define the energy denominators as a 4-index tensor with `np.reshape`.","metadata":{}},{"cell_type":"code","source":"e_ijab = 1 / (e_ab.reshape(1, 1, N_V, 1) + e_ab.reshape(1, 1, 1, N_V) - e_ij.reshape(N_O, 1, 1, 1) - e_ij.reshape(1, N_O, 1, 1))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The same-spin and opposite-spin are the contraction of the ERI and denominators tensors:","metadata":{}},{"cell_type":"code","source":"mp2_os = - np.einsum(\"ijab,ijab,ijab->\", ijab, ijab, e_ijab)\nmp2_ss = - np.einsum(\"ijab,ijab,ijab->\", ijab, ijab - ijab.swapaxes(2, 3), e_ijab)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The results are unchanged with respect to the quadruple-loop approach:","metadata":{}},{"cell_type":"code","source":"print(f\"Opposite-spin MP2 energy: {mp2_os:20.12f}\")\nprint(f\"Same-spin MP2 energy: {mp2_ss:20.12f}\")\nprint(f\"MP2 energy: {mp2_os + mp2_ss:20.12f}\")\n\n# compare with MP2 energy components computed with a quadruple loop\nnp.testing.assert_allclose(mp2_os, e_mp2_os, atol=1e-9)\nnp.testing.assert_allclose(mp2_ss, e_mp2_ss, atol=1e-9)\n# compare with MP2 energy computed by VeloxChem\nnp.testing.assert_allclose(mp2_os + mp2_ss, mp2drv.e_mp2, atol=1e-9)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Bonus: the MP2 one-particle density matrix","metadata":{}},{"cell_type":"markdown","source":"The computation of expectation values of one-electron observables, such as the electric dipole moment, requires the one-particle density matrix (OPDM). In restricted Hartree-Fock theory, the MO basis OPDM is quite simple:\n\n$$\n\\mathbf{D}_{\\mathrm{RHF}}^{\\mathrm{MO}} = \n\\left(\n\\begin{array}{c | c} \n \\begin{array}{c c c} \n 2 & 0 & 0 \\\\ \n 0 & \\ddots & 0 \\\\ \n 0 & 0 & 2 \n \\end{array} & \\mathbf{0} \\\\ \n \\hline \n \\mathbf{0} & \\mathbf{0}\n \\end{array} \n \\right)\n$$\n\nLet us compute the electronic component of the electric dipole moment. First, we obtain the necessary integrals in AO basis and transform them to MO basis. Note that the dipole moment has a 3 Cartesian components: the `ElectricDipoleIntegralsDriver` will compute all of them and the `to_numpy` method thus returns a *list* of AO basis integral matrices.","metadata":{}},{"cell_type":"code","source":"mudrv = vlx.ElectricDipoleIntegralsDriver()\nmu_ao = mudrv.compute(mol, basis).to_numpy()\n# TODO obtain the MO basis electric dipole moment integrals\nmu_mo =","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we construct the RHF density matrix and compute the dipole moment.","metadata":{}},{"cell_type":"code","source":"# TODO form the OPDM of the RHF reference determinant\nDpq = \n\n# TODO compute the RHF electronic component of the electric dipole moment\nmu_rhf =\nprint(mu_rhf)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Lagrangian, amplitudes, and multipliers\n\nComputing the MP2 correction to the density matrix requires not only the MP1 *amplitudes* (wavefunction parameters), but also the MP2 *multipliers*. The MP2 Lagrangian is:\n\n$$\n\\mathcal{L}^{(2)} =\n\\sum_{ijab} \\varepsilon_{ij}^{ab} t_{ij}^{ab\\,(1)}\\lambda^{ij\\,(1)}_{ab} + \n\\sum_{ijab} t_{ij}^{ab\\,(1)} \\langle 0 | \\Phi | _{ij} ^{ab} \\rangle + \n\\sum_{ijab} \\lambda^{ij\\,(1)}_{ab} \\langle _{ij} ^{ab} | \\Phi | 0 \\rangle\n$$\n\nAt its stationary points, the Lagrangian is the MP2 energy. The amplitudes and multipliers that make the Lagrangian stationary *are* \nthe MP1 amplitudes, $t_{ij}^{ab\\,(1)}$, and multipliers, $\\lambda^{ij\\,(1)}_{ab}$:\n\n$$\n\\frac{\\partial \\mathcal{L}^{(2)}}{\\partial \\lambda^{ij\\,(1)}_{ab}} = 0 \\Longleftrightarrow \nt_{ij}^{ab\\,(1)} = - \\frac{\\langle _{ij} ^{ab} | \\Phi | 0 \\rangle}{\\varepsilon_{ij}^{ab}},\\quad\n\\frac{\\partial \\mathcal{L}^{(2)}}{\\partial t_{ij}^{ab\\,(1)}} = 0 \\Longleftrightarrow \n\\lambda^{ij\\,(1)}_{ab} = - \\frac{\\langle 0 | \\Phi | _{ij} ^{ab} \\rangle}{\\varepsilon_{ij}^{ab}}\n$$\n\nFor a closed-shell system, the amplitudes and multipliers are explicitly given as:\n\n$$\nt_{ij}^{ab\\,(1)} = - \\frac{\\langle ab | ij \\rangle}{\\varepsilon_{ij}^{ab}},\\quad\n\\lambda^{ij\\,(1)}_{ab} = - 2\\frac{[2\\langle ij | ab \\rangle - \\langle ij | ba \\rangle]}{\\varepsilon_{ij}^{ab}}\n$$","metadata":{}},{"cell_type":"markdown","source":"### The *unrelaxed* MP2 density matrix\n\nThe expectation value of a one-electron operator $X$ is then:\n\n$$\n\\langle X \\rangle = \\langle 0 | X | 0 \\rangle +\n\\sum_{ijab}\\sum_{klcd} \\lambda^{ij\\,(1)}_{ab} t_{kl}^{cd\\,(1)}\n\\langle 0 | [X, \\tau_{kl}^{cd}] | 0 \\rangle = \n\\langle 0 | X | 0 \\rangle + \n\\sum_{ij} X_{ij} D^{(2)}_{ji} + \\sum_{ab} X_{ab} D^{(2)}_{ba} \n$$\n\nwhere we identifyied, after some algebra, the MP2 density matrix as:\n\n$$\nD^{(2)}_{ij} = - \\frac{1}{2} \\sum_{abk} \\lambda^{ki\\,(1)}_{ab} t_{kj}^{ab\\,(1)},\\quad\nD^{(2)}_{ab} = \\frac{1}{2} \\sum_{ijc} \\lambda^{ij\\,(1)}_{bc} t_{ij}^{ac\\,(1)}\n$$\n\nNote that:\n1. These expressions are *not manifestly symmetric*. In practice, we will symmetrize the MP2 density matrix before computing expectation values.\n2. This is the so-called *unrelaxed* density matrix. We are using *fixed* molecular orbitals, unaffected by electron correlation.\n\nLet us use `np.einsum` to form all these quantities and compute the *unrelaxed* MP2 electronic component of the electric dipole moment.","metadata":{}},{"cell_type":"code","source":"# TODO form the doubles amplitudes, i.e. the MP1 wavefunction coefficients, from the ERI and denominator tensors\nt2 =\n \n# TODO form the doubles multiplitiers, from the ERI and denominator tensors\nl2 =","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We now form the OO and VV blocks of the MP2 OPDM, following the equations above, and update the reference MO density matrix with these contributions:","metadata":{}},{"cell_type":"code","source":"# TODO form the OO block of the MP2 OPDM, don't forget to symmetrize it!\nDij =\n\n# TODO form VV block of the MP2 OPDM, don't forget to symmetrize it!\nDab =\n\n# add the MP2 contributions to the OO and VV blocks\nDpq[:N_O, :N_O] += Dij\nDpq[N_O:, N_O:] += Dab\n\n# test that the total density matrix is symmetric\nnp.testing.assert_allclose(Dpq, np.transpose(Dpq))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, the dipole moment:","metadata":{}},{"cell_type":"code","source":"# TODO compute the MP2 unrelaxed electronic component of the electric dipole moment\nmu_mp2 =\nprint(mu_mp2)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Natural occupations and natural orbitals\n\nThe eigenvectors and eigenvalues of the one-particle density matrix provide the *natural orbitals* and the *natural-orbital occupation numbers*. These are useful, for example, in setting up the active space in multireference calculations.","metadata":{}},{"cell_type":"code","source":"# TODO diagonalize the density matrix\nomega, NOs =\n\n# eigenvalues and eigenvectors are in *ascending* order, so we resort them in *descending* order\nidx = omega.argsort()[::-1] \nomega = omega[idx]\nNOs = NOs[:,idx]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"One way to check that our MP2 unrelaxed OPDM is correct is to sum the natural occupation numbers: this should give us the number of electrons in the system.","metadata":{}},{"cell_type":"code","source":"print(sum(omega))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can also double-check against reference values for the occupations from another code.","metadata":{}},{"cell_type":"code","source":"# reference natural orbital occupation numbers from DALTON\nref_omega = np.array([\n 1.99990540, 1.98720752, 1.97426785, 1.97108868, 1.96924405,\n 0.02241866, 0.02020351, 0.01713431, 0.01024357, 0.00551830,\n 0.00517755, 0.00472951, 0.00414944, 0.00404548, 0.00090056,\n 0.00086293, 0.00060545, 0.00051955, 0.00046726, 0.00045319,\n 0.00039740, 0.00037924, 0.00004262, 0.00003795 \n])\n\nnp.testing.assert_allclose(omega, ref_omega, atol=1e-3)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## References\n\n- Shavitt, I.; Bartlett, R. J. *Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory* Cambridge Molecular Science; Cambridge University Press, 2009.\n- Helgaker, T.; Jørgensen, P.; Olsen, J. *Molecular Electronic-Structure Theory* Wiley, 2000.","metadata":{}}]}