Recently, we showed the basics of quantum chemistry is, and how to do it using Python. This will be a continuation of it. In this article, we will discuss how to write a python script for solving mp2 equations. Make sure you check out the previous article to learn the very basics.
The theory behind it requires years of studying quantum mechanics, mathematics, chemistry, and sometimes even some biology, and combining these fields together to study nature using computer models. This is not what we are going to do here! However, I can give you some hints about how a complicated problem can be made much easier using Python script. We are going to use the PySCF package to perform some calculations.
In the examples in the previous article, you can perform calculations of the energy of a molecule using the Hartree-Fock theory and using Density Functional Theory. These methods are very important and useful, though in some cases, one wants to get to the next level in terms of accuracy.
Table of Contents
What is correlation energy?
After multiple approximations, it is possible to recover about 99% of the exact energy of a molecule. While it may sound excellent, the whole of chemistry is in the remaining 1%, which is also called correlation energy. Among the methods improving on this starting point (i.e., Hartree-Fock theory) is perturbation theory.
The name sounds funny and weird at the same time. You can think of something that disturbs the molecule a bit. But not too much. Let’s say the molecule is placed in a magnetic field. The interaction between them affects the energy and the properties of the molecule. They can be modelled by defining new parameters in the Schrödinger equation.
How to calculate the MP2 correlation energy?
If you read the previous article, you know that the electrons of atoms and molecules live in what we call orbitals. Within Hartree-Fock theory, we can model the orbitals using a set of mathematical functions called a basis set. There are coefficients telling what the weight of each basis function is. These coefficients are written in matrix form, which is called the Fock matrix.
After lots of mathematical operations, one obtains a set of equations that have to be solved. The solution gives the total energy of the molecule and optimises the orbital coefficients that give the best energy. However, these orbitals can be used as input in the perturbation theory approach. The most common method within it carries the names of two scientists – Møller–Plesset perturbation theory. It is usually truncated to the second order and abbreviated as MP2 for conciseness.
Here is an example of how to do it in practice using Python. We will use water as an example. First, we will define the coordinates, then specify the basis set and that the molecule is symmetric. Next, we have to run the Hartree-Fock calculation. Finally, we can also perform the MP2 calculation using the optimized wave function that was just calculated.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
from pyscf import gto, scf, mp mol = gto.Mole() # defining the molecule mol.atom = ''' O -1.82473 1.24653 -0.01294 H -0.85684 1.28190 0.01714 H -2.10358 1.95337 0.58838''' # specifying the basis set mol.basis = 'ccpvdz' mol.symmetry = 1 mol.build() # Computational method # Hartree-Fock theory mf = mol.RHF().run() # MP2 mf_mp2 = mp.MP2(mf).run() |
The code above prints the Hartree-Fock energy, as well as the MP2 energy. Their difference is the MP2 correlation energy.
|
1 2 |
converged SCF energy = -76.0259881597798 E(MP2) = -76.2308195860744 E_corr = -0.204831426294622 |
A note about the units. Probably these numbers do not speak much to you. They are in what we call “atomic units”, also known as hartree units. It is a bit hard to give you a perspective about its magnitude in terms of everyday objects or phenomena, but I will try. If you are into tech, you might know that modern high-end CPUs run at less than 1.5 volts. 1 hartree is equal to about 27 electron-volts. Check out Wikipedia if you want to know more about the units of energy.
What do we mean by the computational cost of calculations?
The example for water finishes immediately. However, some quantum chemical calculations can take hours and even days. Improving the results using second-order perturbation theory involves using the orbitals obtained in a Hartree-Fock calculation. This approach is more accurate, but it is computationally heavy. It depends on the number of integrals that the program needs to calculate in order to reach the final result.
Why are quantum chemistry calculations so computationally heavy?
There are mathematical techniques that introduce some approximations which do not affect the accuracy significantly, but they speed up the calculations. They carry names such as density fitting and resolution of identity. The hardest part of the equations is solving two-electron integrals. Essentially, they depict the electrostatic repulsion between two electrons occupying two possible orbitals. We place electron 1 in orbital 1 and electron 2 in orbital 2. Then each of them is squared in the integral expression. The other set of two-electron integrals describes the so-called exchange interaction, which has no analogue in the macroscopic world. The two electrons are placed each in two different orbitals – electron 1 in orbital 1 and electron 2 in orbital 2. Then, their places are flipped – electron 1 goes to orbital 2, and electron 2 goes to orbital 1. Sounds simple, but the math behind it is very demanding.
How to simplify and speed up the calculations?
The two-electron integrals can be expanded using different functions called an auxiliary basis set. They give little penalty to the accuracy, but the speedup is quite good. Still, it is impossible to calculate a whole protein with millions of atoms using correlated quantum chemistry methods. They are usually performed using classical Newtonian mechanics, where the nuclei are moving in an electrostatic potential created by the electrons. But that is a story for another day.
How big molecules can be studied using MP2?
One little challenge would be to check how the number of atoms affects the speed of the calculation. Or more precisely, the number of electrons. A hydrogen atom only has one electron, but iron has 26. You can guess which one takes longer to calculate. You can test on your own how many seconds it takes to calculate water versus a big porphyrin molecule. Here are the coordinates of its atoms
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
mol.atom = ''' H 3.209599995 3.172309777 0.000000000 N 0.000000000 -2.113600613 0.000000000 N 2.031813914 0.000000000 0.000000000 N 0.000000000 2.113600613 0.000000000 N -2.031813914 0.000000000 0.000000000 C 1.126850920 -2.888627221 0.000000000 C -1.126850920 -2.888627221 0.000000000 C 1.126850920 2.888627221 0.000000000 C -1.126850920 2.888627221 0.000000000 C 2.850160984 -1.083939750 0.000000000 C 2.850160984 1.083939750 0.000000000 C -2.850160984 -1.083939750 0.000000000 C -2.850160984 1.083939750 0.000000000 C 0.683827200 -4.249372023 0.000000000 C -0.683827200 -4.249372023 0.000000000 C 0.683827200 4.249372023 0.000000000 C -0.683827200 4.249372023 0.000000000 C 4.248378966 -0.675836785 0.000000000 C 4.248378966 0.675836785 0.000000000 C -4.248378966 -0.675836785 0.000000000 C -4.248378966 0.675836785 0.000000000 C 2.433989580 -2.416542873 0.000000000 C -2.433989580 -2.416542873 0.000000000 C -2.433989580 2.416542873 0.000000000 C 2.433989580 2.416542873 0.000000000 H 1.341416542 -5.103998738 0.000000000 H -1.341416542 -5.103998738 0.000000000 H 1.341416542 5.103998738 0.000000000 H -1.341416542 5.103998738 0.000000000 H 5.096138011 -1.344322017 0.000000000 H 5.096138011 1.344322017 0.000000000 H -5.096138011 -1.344322017 0.000000000 H -5.096138011 1.344322017 0.000000000 H 3.209599995 -3.172309777 0.000000000 H -3.209599995 -3.172309777 0.000000000 H -3.209599995 3.172309777 0.000000000 H 0.000000000 -1.102285542 0.000000000 H 0.000000000 1.102285542 0.000000000''' |
This is what the molecule looks like if you open it in a molecular viewer such as Avogadro.
What have we learned in this article?
Do you feel more of a quantum chemist after this article? No worries if you don’t. You can always learn more. There is an excellent implementation of a wide variety of quantum chemistry methods available in PySCF that you can test on your own, and you will learn more about it. How about you try some fancier molecule? You can do that in the PyScripter IDE. You can find installation instructions for PySCF inside PyScripter in the previous article about quantum chemistry.
PyScripter is available for download for free at https://www.embarcadero.com/free-tools/pyscripter/free-download
