Introduction to Molecular Modelling: Part 9 (Partial Charges)
I’m back with another part of the series of blog posts that I started last year! This is part 9, here I will cover the calculation of partial charges. As with my previous posts, I will be using GAMESS as the QM software, and Avogadro & wxMacMolPlt as GUI for GAMESS.
What are partial charges?
In chemistry, there are ions and neutral molecules. Ions have a total positive or negative charge, while neutral molecules are often polar. The polarity is due to the differences in electronegativities of different elements. For ionic molecules the formal charge is often drawn on one atom — perhaps the central atom, even though it would be unlikely that all of the charge is concentrated on one atom. The bottom line of all of this is that the electron density of a chemical species (ion, molecule) is generally distorted in some way, which allows us to talk about partial charges. For example, an electronegative atom will have a negative partial charge (as it is drawing electron density toward itself) and vice versa. In some way, the partial charge is therefore, a representation of the distortion of electron density in a chemical species. Partial charges are generally non-integer i.e. fractional.
Partial charges are not real in the sense that they cannot be directly measured through an experiment. There are some experimental techniques like NMR, X-ray scattering, dipole moment measurement which can be used to indirectly calculate partial charges, but as we cannot define exactly where an atom is present, the values given always depend on how we define atom boundaries.
QM softwares like GAMESS or Gaussian can calculate the electron density of the molecule. Using that data, it is possible to assign partial charges to atoms. Here again, the exact values of partial charges depend on the method used to extract the charges from the electron density.
So where are partial charges useful? Calculated partial charges can be often useful for a qualitative understanding of chemical phenomena. For example, imagine a reaction where an alkoxide (RO⁻) is a leaving group. When studying the reaction path, we can choose to look at the partial charge on the oxygen atom, and when it reaches close to 1, we can consider that the leaving group has completely left.
Partial charges are tremendously useful for molecular dynamics. In molecular dynamics, the atoms are considered as a discrete particles which follow Newton’s laws of motion. (In QM the electrons are considered as waves which exist around nuclei.) To model the electrostatic attraction between positively polarized parts and negatively polarized parts of molecules, partial charges are assigned to each atom and the interaction is calculated using Coulomb’s law by summing over atom pairs. This is a simple approach, however, it can give very accurate results if the parameters for molecular dynamics are fitted well.
Calculating partial charges from a QM software
There are a huge number of schemes for calculating partial charges computationally (see this). There are pros and cons of every charge scheme and they all give varying results. I do not have extensive knowledge regarding the different schemes, so I will not go into detail about the various types.
GAMESS can calculate both population-analysis based charges (Mulliken, Löwdin) and electrostatic potential i.e. ESP based charges (CHELPG, Merz-Kollman)
Population analysis based charges take the total wavefunction, and attempts to partition the total number of electrons among different atoms using some scheme. This means the calculation returns the fractional number of electrons which is considered to be present on one atom, i.e. the population. This can then be added to the nuclear charge to get the partial atomic charges. These charges are easy to calculate, however, they can change a lot depending on the method/basis set used, so they are less reproducible and reliable.
ESP based charges calculate the electrostatic potential in a 3D grid around the molecule, from the wavefunction. Then the atomic partial charges are fitted to reproduce the electrostatic potential (and also the total charge). The partial charges can be fitted to more parameters such as the dipole moment or quadrupole moment if required. Sometimes, such fitting can make the partial charges more chemically meaningful. There are different schemes for ESP charges, all of which use different ways to generate the 3D grid. The good thing about these charges are that they are more stable with respect to method and basis set, and they are used for many molecular dynamics force-fields. The only problem with these charges is that they are not reliable for atoms which are away from the grid points.
Population based charges with GAMESS
Mulliken and Lowdin charges are calculated and shown automatically in the output after every single point energy calculation in GAMESS. (Only for methods which have density matrix, which means CCSD(T) won’t show partial charges)
ESP based charges with GAMESS
ESP charges are requested in GAMESS through the $ELPOT group. The keywords to perform ESP charge calculation are IEPOT=1 WHERE=PDC in this group. The first keyword turns on electrostatic potential calculation, and the second keyword requests the use of grids for the calculation.
The grid options are present in $PDC group. The two main keywords in this group are PTSEL and CONSTR. The first keyword gives the choice of the grid. The defaults are PTSEL=GEODESIC, which uses an algorithm by Mark Spackman, and CONSTR=CHARGE, which means the partial charges are fitted to the total charge of the system too.
The Merz-Kollman (or Merz-Singh-Kollman) scheme for partial charges is quite popular. For the original method of Merz, Singh and Kollman, PTSEL=CONNOLLY has to be used. (In Gaussian this is requested with pop=mk). The default grid choice (geodesic) also gives similar values to the Merz-Kollman method, but it is less rotation dependent.
Let’s try to calculate the partial charges with the original Merz-Singh-Kollman method for the HCl molecule. First, optimize an H-Cl molecule with B3LYP/aug-pcseg-1, and get the strucure:
Note that this is a linear molecule, so the automatic generation of internal coordinates (with $ZMAT DLC=.T. AUTO=.T. $END) does not work. Either optimize in Cartesian coordinates or define the one internal coordinate (the bond) with IZMAT(1)=1,1,2 in $ZMAT group.
Next, open the log file from optimization with wxMacMolPlt or Avogadro and write an input file. Edit the file with a text editor:
! File created by the GAMESS Input Deck Generator Plugin for Avogadro
$BASIS GBASIS=APCseg-1 $END
$CONTRL SCFTYP=RHF RUNTYP=ENERGY ISPHER=1 DFTTYP=B3LYP $END
$ELPOT IEPOT=1 WHERE=PDC OUTPUT=NONE $END
$PDC PTSEL=CONNOLLY $END$DATA
Title
C1
Cl 17.0 -7.46498 2.90098 0.00000
H 1.0 -6.16880 2.90098 0.00000
$END
Notice that there is an extra keyword OUTPUT=NONE in $ELPOT group. This means the data of the grid points is not printed in the output, only the partial charges on each atom are printed in the end. Also note that the partial charge calculation is running after a single point calculation here.
Run the input file, and look near the end of the output file:
The partial charge on Cl is -0.2142e and the partial charge on H is +0.2142e. (Here e is the charge of an electron). As you can see, the partial charges match the intuitive notion of electronegativity, because the Cl atom is more electronegative, so it draws the electron density towards it.
Let’s try the same calculation with the acetate anion. An intuitive chemical analysis would suggest that the negative charge in the acetate ion is spread equally between the two oxgyen atoms.
The partial charges at B3LYP/aug-pcseg-1 level are:
As you can see the partial charges on each oxygen atom is about -0.86e. The charge is more negative than -0.5e (if we assume the negative charge being shared equally) is because the oxygen atoms also withdraw some electron density from the carbon atom they are attached to, which leads to a higher partial negative charge. I had mentioned before, that even though we draw the negative charge in structures on one or more atoms, the calculated partial charges do not always match this exactly, even though the trend is generally the same.
Other softwares
There are other softwares such as Multiwfn, which can use a GAMESS single-point calculation log-file as an input and then calculate partial charges. This program has a large number of partial charge algorithms available (such as Hirshfield, Voronoi, CM5, RESP, Gasteiger, modified Mulliken etc.).
Lastly, I want to reiterate the point that it is probably not wise to take any particular computational algorithm for generating partial charges as sacrosanct. All models attempt to capture the unequal distribution of electron density, with the aim that the partial charge values are somehow related to the chemical behaviour that we can see experimentally.
GAMESS: https://www.msg.chem.iastate.edu/gamess/download.html
Avogadro: https://avogadro.cc/
wxMacMolPlt: https://brettbode.github.io/wxmacmolplt/
Multiwfn: http://sobereva.com/multiwfn/
