Physical machine learning outperforms "human learning" in Quantum Chemistry

Anton V. Sinitskiy and Vijay S. Pande (2019)

Highlighted by Jan Jensen

Figure 1 from the paper

The paper presents a method to estimate DFT or CCSD(T) energies (computed using large basis sets) based only on HF/cc-pVDZ densities and energies. In order to avoid overfitting, the method must also estimate the corresponding DFT or CCSD densities. The method is trained and validated on a subset of the QM9 data set (i.e. on relative small molecules). It is first trained on DFT data (using 89K molecules) and then retrained on CC data for a smaller subset (3.6K molecules), both being subsets of the QM9 data set (i.e. relatively small molecules). The input density is evaluated in a 3D grid that is big enough to accommodate the largest molecule in the data set, so a new model would have to be trained for significantly larger molecules.

This “physical machine learning in Quantum Chemistry” (PML-QC) model reaches a mean absolute error in energies of molecules with up to eight non-hydrogen atoms as low as 0.9 kcal/mol relative to CCSD(T) values, which is quite impressive. In fact the authors speculate that 

With ML, it may become not required that an accurate quantum chemical method works fast enough for every new molecule that an end user may be interested in. Instead, the focus shifts to generating highly accurate results only for a finite dataset to be used for training, while the efficiency in practical applications is to be achieved via improvements in DNNs to make them faster and more accurate.

Much will depend how much the number of outliers can be reduced. For example, for PML-QC_DFT ca 5% of molecules have errors greater than 2.6 kcal/mol and this effect can be magnified for relative energies if the individual errors have opposite sign.

This work is licensed under a Creative Commons Attribution 4.0 International License.

Popular Integration Grids Can Result in Large Errors in DFT-Computed Free Energies

Andrea N. Bootsma and Steven E. Wheeler (2019)

Highlighted by Jan Jensen

 Figure 1B from the paper (CC BY-NC-ND 4.0)

 This paper has already been highlighted here and here, so I'll just briefly summarise.

The grid used for the numerical integration in DFT calculations is defined relative to the Cartesian axes, so rotating the molecule will change the integration grid and, hence, the energy. This has been known for some time and, f.eks. Gaussian09 uses a default grid size (Fine, 75,302) where the effect on the electronic energy variation is usually negligible.

Bootsma and Wheeler show that the vibrational entropy and, hence, the free energy is significantly more sensitive to grid size than the electronic energy. Using the Fine grid, the differences in relative free energy changes can be as large as 4 kcal/mol, which could significantly change conclusion regarding mechanisms, etc. The effect comes from the variation in low frequency vibrational modes and the effect can be reduced a little by scaling these frequencies. 

However, the errors really only become acceptable when using the UltraFine grid size, which is the default in Gaussian16, especially combined with frequency scaling (which one should do anyway to get consistent results). If you are using Gaussian09 or some other quantum program to compute relative free energies it is definitely a good idea to look at the default grid size and perform some tests.

Note that if you want to perform such tests yourself, you need to re-optimise the molecule after you rotate it because the gradient is also affected by the rotation.

This work is licensed under a Creative Commons Attribution 4.0 International License.

More DFT benchmarking

Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

Selecting the appropriate density functional for one’s molecular system at hand is often a very confounding problem, especially for non-expert or first-time users of computational chemistry. The DFT zoo is vast and confusing, and perhaps what makes the situation worse is that there is no lack of benchmarking studies. For example, I have made more than 30 posts on benchmark studies, and I made no attempt to be comprehensive over the past dozen years!

One such benchmark study that I missed was presented by Mardirossian and Head-Gordon in 2017.¹ They evaluated 200 density functional using the MGCDB84 database, a combination of data from a number of different groups. They make a series of recommendations for local GGA, local meta-GGA, hybrid GGA, and hybrid meta-GGA functionals. And when pressed to choose just one functional overall, they opt for ωB97M-V, a range-separated hybrid meta-GGA with VV10 nonlocal correlation.

Georigk and Mehta² just recently offer a review of the density functional zoo. Leaning heavily on benchmark studies using the GMTKN55³ database, they report a number of observations. Of no surprise to readers of this blog, their main conclusion is that accounting for London dispersion is essential, usually through some type of correction like those proposed by Grimme.

These authors also note the general disparity between the most accurate, best performing functional per the benchmark studies and the results of the DFT poll conducted for many years by Swart, Bickelhaupt and Duran. It is somewhat remarkable that PBE or PBE0 have topped the poll for many years, despite the fact that many newer functionals perform better. As always, when choosing a functional caveat emptor.

References

1.  Mardirossian, N.; Head-Gordon, M., “Thirty years of density functional theory in computational chemistry: an overview and extensive assessment of 200 density functionals.” Mol. Phys. 2017, 115, 2315-2372, DOI: 10.1080/00268976.2017.1333644.

2. Goerigk, L.; Mehta, N., “A Trip to the Density Functional Theory Zoo: Warnings and Recommendations for the User.” Aust. J. Chem. 2019, ASAP, DOI: 10.1071/CH19023.

3. Goerigk, L.; Hansen, A.; Bauer, C.; Ehrlich, S.; Najibi, A.; Grimme, S., “A look at the density functional theory zoo with the advanced GMTKN55 database for general main group thermochemistry, kinetics and noncovalent interactions.” Phys. Chem. Chem. Phys. 2017, 19, 32184-32215, DOI: 10.1039/C7CP04913G.

'
This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Triplet-Tuning: A Novel Non-Empirical Construction Scheme of Exchange Functionals

Zhou Lin and Troy Van Voorhis (2018)

Highlighted by Jan Jensen

The absolute error of the optical band gap computed with the CC-PVDZ basis set

The difference between the lowest triplet and singlet energy (E_T) can be calculated both by UDFT and TDDFT but give different results because we don't know the exact density functional. So the authors suggest that functional can be improved by minimising this difference, i.e. without comparison to experimental data. 

Indeed, such a triplet tuned (TT) functional "provide more accurate predictions for key observables in photochemical measurements, including but not limited to E_T, optical band gaps (E_S), singlet–triplet gaps (∆E_ST), and ionization potentials (I)" for a set of 100 organic molecules. Two parameters in the PBE exchange functional, the fraction of short-range HF exchange and the range-separation parameter, were adjusted.

One thing that is not clear to me is if the equivalence of UDFT and TDDFT is valid only for the exact density. If so, this would only be valid in the complete basis set limit. Either way, the results clearly improve using the CC-PVDZ basis set.


This work is licensed under a Creative Commons Attribution 4.0 International License.

Rational Density Functional Selection Using Game Theory

Suzanne McAnanama-Brereton and Mark P. Waller (2018)
Highlighted by Jan Jensen

Copyright 2018 American Chemical Society

This paper describes a method for finding the optimum combination of DFT functional and basis set for a particular calculation by using game theory, more specifically by finding the corresponding Nash equilibria. As the authors write:

A Nash equilibrium (NE) is an optimal strategy in game theory, given that all players have differing payoffs for different strategies. It is a strategy that, when all players have full knowledge of their own and their opponents’ payoffs for each strategy, every rational player will choose.

For this application the three "players" are 1) the complexity or computational cost of the method, 2) the accuracy of the method for a benchmark set, and 3) the similarities of the system of interest with the molecules in the benchmark set.

Scores for each player are entered into a payoff matrix and then fed to the Gambit program that computes the Nash equilibria, which are then converted into a score for each combination of functional and basis set. The method, called Decider, is made available as via a web-interface.

The method is for the benzene dimerization energy using the S22 benchmark set (where the benzene dimer has been removed). The top 81 functional/basis set combinations are given and the top five, middle five, and bottom five combinations are tested, i.e. the benzene dimerization energy is computed and compared to the reference value. As expected, the top five combinations are among both the fastest and most accurate of the 15 combinations tested.

Because the molecules in the benchmark set are all quite small, Decider invariably suggested the def2-QZVP basis because it affordable for these kinds of systems and provides accurate results (in fact this happens for every molecule pair I have tried to far).

As far as I can tell the version available via the web interface is based on the S22 benchmark, so it can only be used for interaction energies at present. Every molecule I have tried to far has led to def2-QZVP as the only basis set suggestion.  Unfortunately, it doesn't appear that the source code is available.  The source code is available on Github.

The authors have also set up a Turing test of the method.


This work is licensed under a Creative Commons Attribution 4.0 International License.

Is the Accuracy of Density Functional Theory for Atomization Energies and Densities in Bonding Regions Correlated?

Brorsen, K. R.; Yang, Y.; Pak, M. V.; Hammes-Schiffer, S., J. Phys. Chem. Lett. 2017, 8, 2076-2081
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

I recently blogged about a paper arguing that modern density functional development has strayed from the path of improving density description, in favor of improved energetics. The Medvedev paper¹ was met with a number of criticisms. A potential “out” from the conclusions of the work was that perhaps molecular densities do not fare so poorly with more modern functionals, following the argument that better energies might reflect better densities in bonding regions.

The Hammes-Schiffer group have now examined 14 diatomic molecules with the goal of testing just this hypothesis.² They subjected both homonuclear diatomics, like N₂, Cl₂, and Li₂, and heteronuclear diatomics, like HF, LiF, and SC, to 90 different density functionals using the very large aug-cc-pCVQZ basis set. Using the CCSD density as a reference, they examined the differences in the densities predicted by the various functional both along the internuclear axis and perpendicular to it.

The 20 functionals that do the best job in mimicking the CCSD density are all hybrid GGA functionals, along with the sole double hybrid functional included in the study (B2PLYP). These functionals date from 1993 to 2012. The 20 functionals that do the poorest job include functionals from all rung-types, and date from 1980-2012. A very slight upward trend can be observed in the density error increasing with development year, while the error in the dissociation energy clearly is decreasing over time.

They note that six functionals of the Minnesota-type, those that are highly parameterized and of recent vintage, perform very poorly at predicting atomic densities, but do well with the diatomic densities.

Hammes-Schiffer concludes that their diatomic results support the general trend noted by Medvedev’s atomic results, that density description is lagging in more recently developed functionals. I’d add that this trend is not as dramatic for the diatomics as for atoms.

They pose what is really the key question: “Is the purpose to approximate the exact functional or simply to provide chemists with a useful tool for exploring chemical systems?” Since, as they note, the modern highly parameterized functionals have worked so well for predicting energies and geometries, “the observation that many modern functionals produce incorrect densities could be of no great consequence for many studies”. Nonetheless, “the ultimate goal is still to obtain both accurate densities and accurate energies”.

References

1) Medvedev, M. G.; Bushmarinov, I. S.; Sun, J.; Perdew, J. P.; Lyssenko, K. A., "Density functional theory is straying from the path toward the exact functional." Science 2017, 355, 49-52, DOI: 10.1126/science.aah5975.

2) Brorsen, K. R.; Yang, Y.; Pak, M. V.; Hammes-Schiffer, S., "Is the Accuracy of Density Functional Theory for Atomization Energies and Densities in Bonding Regions Correlated?" J. Phys. Chem. Lett. 2017, 8, 2076-2081, DOI: 10.1021/acs.jpclett.7b00774.

'
This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Progress in DFT development and the density they predict

Medvedev, M. G.; Bushmarinov, I. S.; Sun, J.; Perdew, J. P.; Lyssenko, K. A., "Density functional theory is straying from the path toward the exact functional." Science 2017, 355, 49-52
Hammes-Schiffer, S., "A conundrum for density functional theory." Science 2017, 355, 28-29
Korth, M., "Density Functional Theory: Not Quite the Right Answer for the Right Reason Yet." Angew. Chem. Int. Ed. 2017, 56, 5396-5398
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

“Getting the right answer for the right reason” – how important is this principle when it comes to computational chemistry? Medvedev and co-workers argue that when it comes to DFT, trends in functional development have overlooked this maxim in favor of utility.¹ Specifically, they note that

There exists an exact functional that yields the exact energy of a system from its exact density.

Over the past two decades a great deal of effort has gone into functional development, mostly in an empirical way done usually to improve energy prediction. This approach has a problem:

[It], however, overlooks the fact that the reproduction of exact energy is not a feature of the exact functional, unless the input electron density is exact as well.

So, these authors have studied functional performance with regards to obtaining proper electron densities. Using CCSD/aug-cc-pwCV5Z as the benchmark, they computed the electron density for a number of neutral and cationic atoms having 2, 4, or 10 electrons. Then, they computed the densities with 128 different functionals of all of the rungs of Jacob’s ladder. They find that accuracy was increasing as new functionals were developed from the 1970s to the early 2000s. Since then, however, newer functionals have tended towards poorer electron densities, even though energy prediction has continued to improve. Medvedev et al argue that the recent trend in DFT development has been towards functionals that are highly parameterized to fit energies with no consideration given to other aspects including the density or constraints of the exact functional.

In the same issue of Science, Hammes-Schiffer comments about this paper.² She notes some technical issues, most importantly that the benchmark study is for atoms and that molecular densities might be a different issue. But more philosophically (and practically), she points out that for many chemical and biological systems, the energy and structure are of more interest than the density. Depending on where the errors in density occur, these errors may not be of particular relevance in understanding reactivity; i.e., if the errors are largely near the nuclei but the valence region is well described then reactions (transition states) might be treated reasonably well. She proposes that future development of functionals, likely still to be driven by empirical fitting, might include other data to fit to that may better reflect the density, such as dipole moments. This seems like a quite logical and rational step to take next.

A commentary by Korth³ summarizes a number of additional concerns regarding the Medvedev paper. The last concern is the one I find most striking:

Even if there really are (new) problems, it is as unclear as before how they can be overcome…With this in mind, it does not seem unreasonable to compromise on the quality of the atomic densities to improve the description of more relevant properties, such as the energetics of molecules.

Korth concludes with

In the meantime, while theoreticians should not rest until they have the right answer for the right reason, computational chemists and experimentalists will most likely continue to be happy with helpful answers for good reasons.

I do really think this is the correct take-away message: DFT does appear to provide good predictions of a variety of chemical and physical properties, and it will remain a widely utilized tool even if the density that underpins the theory is incorrect. Functional development must continue, and Medvedev et al. remind us of this need.

References

1) Medvedev, M. G.; Bushmarinov, I. S.; Sun, J.; Perdew, J. P.; Lyssenko, K. A., "Density functional theory is straying from the path toward the exact functional." Science 2017, 355, 49-52, DOI: 10.1126/science.aah5975.

2) Hammes-Schiffer, S., "A conundrum for density functional theory." Science 2017, 355, 28-29, DOI: 10.1126/science.aal3442.

3) Korth, M., "Density Functional Theory: Not Quite the Right Answer for the Right Reason Yet." Angew. Chem. Int. Ed. 2017, 56, 5396-5398, DOI: 10.1002/anie.201701894.

'
This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Local Fitting of the Kohn−Sham Density in a Gaussian and Plane Waves Scheme for Large-Scale Density Functional Theory Simulations

Dorothea Golze, Marcella Iannuzzi, and Jürg Hutter (2017)

Highlighted by Michael Banck

Reprinted (adapted) with permission from Dorothea Golze, Marcella Iannuzzi, and Jürg Hutter. Journal of Chemical Theory and Computation, 2017 ASAP,  Copyright 2017 American Chemical Society.

Hutter et al. have published their LRIGPW (local resolution-of-the-identity gaussian and plane waves method) paper in JCTC. The image above taken from that paper highlights that much of the total runtime for conventional GPW (the main method implemented in the CP2K package) is spent on the description of the total charge density on real-space grids ("GPW grid", dark blue). Can you spot the orange bars for the same work done in the LRIGPW approach? This takes CP2K another big step forward.

Effect of Complex-Valued Optimal Orbitals on Atomization Energies with the Perdew–Zunger Self-Interaction Correction to Density Functional Theory

Susi Lehtola, Elvar Ö. Jónsson, and Hannes Jónsson J. Chem. Theory Comput. in press (2016)

Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

One of the biggest problems facing DFT is that of self-interaction: each electron effectively interacts with itself, because the potential derives from the total charge density of the system. This is not an issue for the exact (unknown) density functional, or for Hartree-Fock, but is the cause of significant error in many DFT functionals. Approaches such as DFT+U[1],[2],[3] and hybrid functionals (far too many to reference !) are aimed in part at fixing this problem.

Probably the earliest attempt to remove this error is the self-interaction correction of Perdew and Zunger[4] which corrects the potential for each Kohn-Sham orbital, complicating the calculation considerably over a standard DFT calculation. (Ironically, this paper, which has over 11,000 citations, is best known for its appendix C, where a parameterisation of the LDA XC energy is given.) However, this process is notoriously slow to converge and is not widely used.

A recent paper[5] showed that, even for isolated molecules, complex orbitals were required to achieve convergence, and this approach has now been tested for atomisation energies of a standard set of 140 molecules[6]. The tests compare the new complex SIC implementation against the standard, real implementation, as well as various GGAs, hybrid functionals and meta-GGAs. The complex SIC, when coupled with the PBEsol functional[7], gives good results (though ironically the PBEsol functional was developed to improve PBE for solids). Not surprisingly, the best results are from hybrids, but meta-GGA improves the energies almost as well.

This study highlights the problem with DFT at the moment: there are many different approaches, which often work well for specific problems. SIC is cheaper than hybrid calculations, and can be important for charge transfer problems (and Rydberg states). The results for convergence and complex orbitals are interesting, but based on these results, I would use meta-GGA for atomisation energies, as a good compromise between accuracy and cost (almost the same as GGA).

[1] Phys. Rev. B 52, R5467 (1995) DOI:10.1103/PhysRevB.52.R5467

[2] Phys. Rev. B 57, 1505 (1998) DOI:10.1103/PhysRevB.57.1505

[3] Int. J. Quantum Chemistry 114, 14 (2014) DOI:10.1002/qua.24521

[4] Phys. Rev. B. 23, 5048 (1981) DOI:10.1103/PhysRevB.23.5048

[5] J. Chem. Theory Comput. 12, 3195 (2016) DOI:10.1021/acs.jctc.6b00347

[6] J. Chem. Theory Comput. in press (2016) DOI:10.1021/acs.jctc.6b00622

[7] Phys. Rev. Lett. 100, 136406 (2008) DOI:10.1103/PhysRevLett.100.136406

Systematic Error Estimation for Chemical Reaction Energies

Gregor N. Simm and Markus Reiher (2016)
Contributed by Jan Jensen

Simm and Reiher present an approach to estimate the error of density functionals parameterized for a specific system.  System-specific parameterization can yield very accurate predictions for the training set but the applicability of the resulting functional to any other systems, including closely related ones, is not guaranteed precisely because the training set is so specific.  However, by Simm and Reiher provide confidence intervals for each result to assess whether the results are reliable. 

Four functional parameters are optimized to best fit 23 CCSD(T) reaction energies of a model system. As one would expect the functional outperforms several other standard functionals for that particular model system. The real question is whether this performance translates to the real system of interest where CCSD(T) benchmarks aren’t available.  To answer this Simm and Reiher generate a set of new values for one of the parameters that lead to a range of predicted reaction energies from which the uncertainty (standard deviations) in the prediction can be computed.  

For the model system the error is within two standard deviations for 21 of the 23 reaction energies. When the functional is then applied to a slightly larger system, for which CCSD(T) reference values also could be computed, the error is within two standard deviations for 17 of the 20 reaction energies. Thus for most of the reaction energies this approach can be used to reliably gauge the accuracy of the results computed using the system dependent density functional.

It should be noted that the current approach only applies to linear parameters, which is the reason the analysis what only performed on one of the four optimized parameters. But it should be possible to extend it to other parameter types, in which case one also needs to address the large number of parameter-combinations that need to be checked.  

This work is licensed under a Creative Commons Attribution 4.0

Different approaches to creating (meta-GGA) DFT functionals

Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

Within the DFT community, John Perdew’s idea of the Jacob’s ladder of accuracy[1] starts with LDA, moves to GGAs with the inclusion of the gradient of the electron density, and is further extended (to rungs three and four in the ladder analogy) with meta-GGA, where kinetic energy density is added, and hybrids, where exact exchange plays a role. Although meta-GGAs have been around for ten to fifteen years, they are only starting to become widely used. I will compare two examples which both seem promising, but also encapsulate the two most common approaches to functional creation.

Meta-GGAs are promising because the kinetic energy density allows some discrimination between areas with one or two electrons (in this way, they are similar in some ways to the electron localisation function, or ELF, which can be used to analyse bonding[2]). This gives some hope that they may be able to fit both strong and weak bonding as well as possibly mitigating the self-interaction error that plagues DFT.

The Minnesota meta-GGA MN15-L[3] takes a very complex functional form with 58 parameters and a non-separable form for the exchange-correlation (adding an extra correlation functional to this) and fits it to a portion of a very large database (I estimate that there are at least 900 entries in the database, covering many different chemical properties, specifically including solid state properties, transition barriers and weak interactions). The resulting functional is local (no hybrid or non-local van der Waals terms were included) and produces extremely small errors when compared to those parts of the database which were not used in fitting. Notably, it out-performs functionals with exchange and van der Waals included.

By contrast, the SCAN functional[4] uses only seven parameters (close to, if not at, the minimal number for a meta-GGA) and includes various physically motivated norms and constraints for the electron gas. The early progress in GGAs was made by satisfying important constraints, so this is seen as a good route to reliability. This functional was also tested on various databases, particularly focussing on solid-state and weak interactions. It has excellent agreement with these, though was published before the MN15-L functional so is not compared directly. In a follow-up paper[5] the performance of the SCAN functional for band gaps is found to be good, though it does not calculate Kohn-Sham gaps, but gaps within a generalized Kohn-Sham theory (a distinction which I don’t have time to discuss here; I may write another blog on this, as it is relevant to hybrid functionals among other things.)

What can we learn from this ? Both functionals perform well in the tests which are published. The functional that you choose will depend in part, as always, on which system you wish to study; however, both of these functionals show some promise in being widely applicable. Your choice will also depend on your attitude to fitting[6]: is a reasonable functional form with many parameters something that you trust, or do you prefer to be more prescriptive, and deal with fewer parameters ? Fifty years since its inception, DFT is still developing, communities are still somewhat divided in the approaches that they take to functional development, but there are an increasing number of ways to achieve efficiency and accuracy.

[1] My colleague, Mike Gillan, reckons that we should instead talk about wrestling Jacob when considering how to improve DFT functionals.
[2] J. Chem. Phys. 92, 5397 (1990) DOI: 10.1063/1.458517
[3] J. Chem. Theor. Comput. 12, 1280 (2016) DOI: 10.1021/acs.jctc.5b01082
[4] Phys. Rev. Lett. 115, 036402 (2015) DOI: 10.1103/PhysRevLett.115.036402
[5] Phys. Rev. B 93, 205205 (2016) DOI: 10.1103/PhysRevB.93.205205
[6] The oft-quoted maxim about fitting an elephant with four parameters has been put into practice (see the original paper here and a nice write-up and a python implementation here)

Consistent structures and interactions by density functional theory with small atomic orbital basis sets

Stefan Grimme, Jan Gerit Brandenburg, Christoph Bannwarth, and Andreas Hansen (2015)
Contributed by Jan Jensen

B3LYP/6-31G* is still the de facto default level of theory for geometry optimizations of large-ish (ca 50-200 atoms) molecules and this paper introduces a cheaper, more accurate replacement called PBEh-3c. PBEh-3c is a basically PBE0/def2-SV(P) with added dispersion and BSSE corrections, where both functional and basis set has been modified slightly (for B-Ne in the case of the basis set).

As the authors write

Most striking is the roughly “MP2-quality” (or slightly better) obtained for the non-covalent complexes in the S22/S66 sets and equilibrium structures ($B_e$ values) for medium-sized organic molecules in the ROT34 set.

For example, the S22 and S66 geometries can be reproduced with an RMSD of 0.08 and 0.05 Å, respectively.  But this method is aimed at the entire periodic table and bond lengths between heavier atoms are also tested and well reproduced.

Though not specifically designed for it the method is also tested for intermolecular interaction energies, reaction energies, barrier heights.  Here PBEh-3c doesn't always outperform B3LYP/def2-SV(P) or M06-2X/def2-SV(P), but when it does it's typically a very significant improvement.  For example, the S30L (which includes host-guest complexes with multiple hydrogen bonds and/or charged systems) interaction energies are reproduced with an MAD of 3.4 kcal/mol, compared to 7.4 and 25.9 kcal/mol for M06-2X/def2-SV(P) and B3LYP/def2-SV(P), respectively.  3.3 kcal/mol may still sound like a lot but, for comparison, the corresponding MAD for PW6B95-D3/def2-QZVP(D) is 2.5 kcal/mol.

This work is licensed under a Creative Commons Attribution 4.0

Metal oxidation states in the oxygen-evolving complex: A computational challenge

V. Krewald, M. Retegan, N. Cox, J. Messinger, W. Lubitz, S. DeBeer, F. Neese and D.A. Pantazis,
Chemical Science 2015, early-view (Open-Access)
Contributed by Marcel Swart

The determination of oxidation states of metals in biological systems is not an easy task and can lead to surprises or controversies. This was for instance shown last year when an unusual molybdenum(III) oxidation state in the nitrogenase enzyme was reported [1], or the year before [2] for a Sc-capped iron-oxygen complex where the crystal structure shows a iron(III) oxidation state unlike the expected iron(IV) state based on solution data [3].
The situation is infinitely more complex when more than one metal ion is present, as is the case in the oxygen evolving complex (OEC) of photosystem II that is studied here through a combination of computational chemistry and spectroscopy. The OEC contains 4 manganese and 1 calcium as the active species to convert water into oxygen. Many studies have been performed on the five stages along the catalytic cycle, often with contradictory conclusions about the geometry, electronic structure, oxidation and protonation states involved.

Summary of the catalytic cycle with the five stages (reproduced with permission from Chem. Sci., 2015, Advance Article, DOI: 10.1039/C4SC03720K; Published by The Royal Society of Chemistry) 

A breakthrough was a recent atomic resolution crystal structure [4] that largely confirmed the theoretical predictions by Siegbahn [5] about the geometrical positioning of the manganese ions and the oxygens. Still, a detailed understanding of the oxidation states was lacking, which is being explored in the paper by Pantazis and co-workers through a combination of a variety of experimental and theoretical techniques.

Typical model system used in the calculations (reproduced with permission from Chem. Sci., 2015, Advance Article, DOI: 10.1039/C4SC03720K; Published by The Royal Society of Chemistry) 

It is this combination of both theory and experiment, and the systematic exploration of protonation states, open- vs. closed cubane structure, distribution of +3/+4 oxidation states over the different manganese ions that makes this a complete and convincing story for the assignment of high-valent (HV) oxidation states along the catalytic cycle.

One example of the different possible distributions of oxidation states (III or IV) and spin states (1/2, 5/2 or 13/2), for state 2 models (reproduced with permission from Chem. Sci., 2015, Advance Article, DOI: 10.1039/C4SC03720K; Published by The Royal Society of Chemistry) 

The HV assignments made are consistent with the experimental data obtained so far, including a very recent “radiation-damage-free” crystal structure [6], for the first three stages of the catalytic cycle (S1-S3). What remains to be explored is how the enzyme produces the dioxygen molecule in the S4 stage and returns back to the resting state S0. Without any doubt, further unexpected findings will be coming along for this intriguing catalytic cycle.

[1] Chem. Sci. 2014, 5, 3096-3103, DOI: 10.1039/C4SC00337C 
[2] Chem. Commun. 2013, 49, 6650-6652, DOI: 10.1039/c3cc42200c
[3] Nature Chem. 2010, 2, 756-759, DOI: 10.1038/nchem.731
[4] Nature 2011, 473, 55-60, DOI: 10.1038/nature09913
[5] Acc. Chem. Res. 2009, 42, 1871-1880, DOI: 10.1021/ar900117k
[6] Nature 2015, 517, 99-103, DOI: 10.1038/nature13991


This work is licensed under a Creative Commons Attribution 3.0 Unported License. 

Orbital-free density functional theory implementation with the projector augmented-wave method

Lehtomäki, I. Makkonen, M. A. Caro, A. Harju, and O. Lopez-Acevedo, J. Chem. Phys. 141, 234102 (2014)
Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

Orbital-free DFT[1,2] is a linear-scaling approach[3] to DFT which returns to the original spirit of DFT and Thomas-Fermi theory, in seeking only the charge density of the system.  The potential computational saving is immense, as we only have to consider one three dimensional variable rather than enough variables to describe the orbitals of the system. As there is no diagonalisation or orthogonalisation step, we achieve linear scaling.  However, there are two considerable problems with the method: first, the kinetic energy functional, which is not known; second, the potential due to the ions.  It is this second problem that a recent paper[4] seeks to address, by using the projector-augmented wave (PAW) method[5].

There has been considerable research on possible kinetic energy functionals for use in OFDFT, with two key approximations: the Thomas-Fermi (which uses a simple power of the density); and the von Weizsacker (which considers the gradient of the density as well as the density).  Most OFDFT approaches mix these at some level, and research into new functionals is on-going.

The second problem, with potentials due to the ions, is significant.  Since the only variable considered by the OFDFT method is the density, which is a local function, the potentials used must also be local. (In solid-state approaches, pseudopotentials are almost always non-local, that is functions of two points in space, which is possible as we are acting on different orbitals with different energies; the non-local form makes the pseudopotentials significantly more transferrable.)  There has been work on developing transferrable local potentials[6] for OFDFT, but this still presents a significant challenge.

Given this problem, an approach that uses the bare Coulomb potential is attractive.  Most OFDFT approaches use a grid-based approach, which makes the use of the bare Coulomb potential challenging, and the proposal to introduce PAWs for OFDFT is attractive.  But PAWs have the same requirement that non-local pseudopotentials do: they use multiple projector functions in multiple angular momentum channels, and are generally non-local.

The authors do not address this question clearly: there is a passing statement that all quantities  are generated as for a 1s orbital (i.e. spherically symmetric), and a single note that only one partial wave is used.  This assumption really deserves more discussion, and investigation: since the projectors are operating on the density (or more properly its square root), it is not clear whether the PAW method can be used faithfully in this context.

The results presented lead to further questions.  While isolated binary molecules are represented reasonably well, this is not a significant test of the method.  The bulk results are more concerning: the lattice constant and bulk modulus for lithium and carbon are very poorly described, and the authors introduce a new parameter to the kinetic energy, which they adjust to match the lattice constant.

While the idea of combining PAWs and OFDFT is an interesting one, the present paper leaves many unanswered questions, most relating to the implementation.  It is far from clear that the use of a single partial wave for PAW offers any advantage over the local pseudopotential approaches already proposed for OFDFT, and while the implementation in an existing code is simple, much more work is needed to test the reliability of the method.

[1] L. W. Wang and M. P. Teter, Phys. Rev. B 45, 13196 (1992) DOI:10.1103/PhysRevB.45.13196

[2] M. Pearson, E. Smargiassi, and P. A. Madden, J. Phys.: Condens. Matter 5, 3221 (1993) DOI:10.1088/0953-8984/5/19/019

[3] D. R. Bowler and T. Miyazaki, Rep. Prog. Phys. 75, 036503 (2012) DOI:10.1088/0034-4885/75/3/036503

[4] J. Lehtomäki, I. Makkonen, M. A. Caro, A. Harju, and O. Lopez-Acevedo, J. Chem. Phys. 141, 234102 (2014) DOI:10.1063/1.4903450

[5] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994) DOI:10.1103/PhysRevB.50.17953

[6] S. Watson, B. J. Jesson, E. A. Carter, and P. A. Madden, Europhys. Lett. 41, 37 (1998) DOI:10.1209/epl/i1998-00112-5

Assessment of theoretical procedures for a diverse set of isomerization reactions involving double-bond migration in conjugated dienes

Yu, L.-J.; Karton, A. Chem. Phys. 2014, 441, 166-177
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

With the proliferation of density functionals, selecting the functional to use in your particular application requires some care. That is why there have been quite a number of benchmark studies (see these posts for some examples). Yu and Karton have now added to our benchmark catalog with a study of π-conjugation.¹

They looked at a set of 60 reactions which involve a reactant with π-conjugation and a product which lacks conjugation. A few examples, showing examples involving linear and cyclic systems, are shown in Scheme 1.

Scheme 1.

The reaction energies were evaluated at W2-F12, which should have an error of a fraction of a kcal mol^-1. Three of the reactions can be compared with experimental values, and difference in the experimental and computed values are well within the error bars of the experiment. It is too bad that the authors did not also examine 1,3-cyclohexadiene → 1,4-cyclohexadiene, a reaction that is both of broader interest than many of the ones included in the test set and can also be compared with experiment.

These 60 reactions were then evaluated with a slew of functionals from every rung of Jacob’s ladder. The highlights of this benchmark study are that most GGA and meta-GGA and hybrid functionals (like B3LYP) have errors that exceed chemical accuracy (about 1 kcal mol^-1). However, the range-separated functionals give very good energies, including ωB97X-D. The best results are provided with double hybrid functionals. Lastly, the D3 dispersion correction does generally improve energies by 10-20%. On the wavefunction side, SCS-MPs gives excellent results, and may be one of the best choices when considering computational resources.

References

(1) Yu, L.-J.; Karton, A. "Assessment of theoretical procedures for a diverse set of isomerization reactions involving double-bond migration in conjugated dienes," Chem. Phys. 2014, 441, 166-177, DOI:10.1016/j.chemphys.2014.07.015.

This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Energy-conserving ab initio molecular dynamics

Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

Energy conservation in molecular dynamics is an important test of the propagators and parameters being used.  In first principles Born-Oppenheimer molecular dynamics (BOMD) this is complicated because of the need for self-consistent iteration to a consistent ground state charge density and potential.  The self-consistency cycle introduces errors since it is not possible to achieve a perfect ground state.  Moreover, to improve efficiency, the starting charge density is often predicted from previous charge densities; while this can results in impressive speed gains, it breaks the time-reversibility of the dynamics and will introduce an unavoidable energy drift.

This blog will discuss a set of papers[1-4] from Anders Niklasson in Los Alamos, which seek to remove this problem, and which have also shown a way to avoid the self-consistency cycle entirely.  The work is perhaps less well-known that it should be, and offers an excellent route to efficiency.  The recent linear scaling DFT molecular dynamics that we have published (going to over 32,000 atoms with DFT)[5] used the time-reversible formulation to conserve energy.

The time reversibility is built on an extended Lagrangian (and is known as XL-BOMD) which introduces a new dynamical variable, or set of variables, for the electronic degrees of freedom (whether charge density or density matrix) alongside the nuclear degrees of freedom, in a trick reminiscent of Car-Parrinello MD[6].  This degree of freedom, written below as P, evolves harmonically in time, depending on a frequency (set by the user) and the difference to the ground state charge density (or whatever is being used).

This degree of freedom, P, can be propagated in a time-reversible manner with an accuracy that only depends on the timestep.  The BOMD can thus be made time-reversible by starting the self-consistency cycle from P – resulting in stability and energy conservation in MD over long timescales.  On its own, this is a very valuable advance.

However, in a set of papers[2-4], Niklasson takes the idea further.  He suggests that, given the right functional for the energy, the self-consistency cycle can be removed.  This involves performing a single diagonalisation/minimisation without self-consistent update, and using P as the input charge density.  The scheme has been developed, with the latest paper[4] generalising the earlier work, and showing that the approach is applicable to systems regardless of gap, provided that an adiabatic separation between electronic and nuclear degrees of freedom can be made.

This work provides techniques to accelerate ab initio MD significantly; we have found that, in combination with linear scaling DFT, we can access picoseconds of MD for tens of thousands of atoms even with a simple initial implementation.  The two approaches together promise to make large-scale, long time accurate MD significantly more affordable.

[1] Phys. Rev. Lett. 100, 123004 (2008) DOI:10.1103/PhysRevLett.100.123004

[2] J. Chem. Phys. 139, 214102 (2013) DOI:10.1063/1.4834015

[3] J. Chem. Phys. 140, 044117 (2014) DOI:10.1063/1.4862907

[4] J. Chem. Phys. 140, 164123 (2014) DOI:10.1063/1.4898803

[5] arxiv 1409.6085 and J. Chem. Theor. Comput. in press (2014) DOI:10.1021/ct500847y