Friday, July 11, 2014

Error Estimates for Solid-State Density-Functional Theory Predictions: An Overview by Means of the Ground-State Elemental Crystals

K. Lejaeghere, V. Van Speybroeck, G. Van Oost & S. Cottenier Critical Reviews in Solid State and Materials Sciences 2014, 39, 1-24
Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

The question of how to characterise the accuracy of a computer code is a difficult one, and I have touched on these issues before (here, here and here for instance). However, given the large number of codes available, it should be possible to compare them to each other and to experiment (or higher level calculations) to test them. This is a well-established process in the quantum chemistry community, where there are various test sets for different properties, including enthalpies of formation(G97/2)[1], weak bonding (S22)[2] (and a vast database of different properties[3]).

A recent paper[4] and associated web site[5] offers a first approach for solid state codes, with the comparison based on the differences between all-electron and pseudopotential calculations. The presumption here is that all-electron calculations are the touchstone (though the website notes that the all-electron results have been refined to use extremely accurate tolerances and small muffin-tin radii, so there is clearly room for improvement in any method).

The idea of the comparison is to calculate binding energy curves for most elements in the periodic table, and from these curves to derive a single number which characterises the deviation from all-electron results. The deviation per element can be viewed, as can the deviation from experiment for the all-electron calculations. The deviation, delta, is defined by an integral over all calculated values of volume, and thus includes implicitly both the lattice constant and the bulk modulus, though not the cohesive energy. Most results shown on the website are for plane-wave codes, which generally perform rather well (the old norm-conserving FHI pseudopotentials are less accurate, though, and should be treated with care).

The approach is a good one, though a little heavy on the tester, and scripts to perform the necessary calculations are made freely available. However, the choice of the elemental form makes the tests rather restricted: there is no way to examine different types of bonding or different oxidation state, for instance. It is quite easy to imagine developing a set of these test suites for different purposes in solid state codes, just as there are different test sets in quantum chemistry.

The accuracy tests seem to be most illuminating for the pseudopotentials, rather than the codes themselves, and I think that it would be of immense value to the community if the pseudopotential generation details were made available. This should not just include the core radii and the reference configurations, but also a clear description of the pseudopotential algorithm (or an appropriate reference along with details).

There is something of a danger of using codes and supplied pseudopotential libraries as black boxes: there is a need to test parameters, though it’s rare (and I am happy to acknowledge that I don’t do it as much as I should). This paper and associated developments should go some way to standardising plane wave codes, and giving quantitative information on their reliability.

Note As I was writing this, a new paper in Science has just been published which takes a different view, and compares results from different functionals; it’s worth a read[6]

[1] J. Chem. Phys. 106, 1063 (1997) DOI:10.1063/1.473182
[2] Phys. Chem. Chem. Phys. 8, 1985 (2006) DOI: 10.1039/B600027D
[3] and see arXiv
[4] Crit. Rev. Sol. Stat. Mat. Sci. 39, 1–24. DOI:10.1080/10408436.2013.772503
[6] Science 345, 197 (2014) DOI:10.1126/science.1253486

Synthesis, Characterization, and Properties of [4]Cyclo-2,7-pyrenylene: Effects of Cyclic Structure on the Electronic Properties of Pyrene Oligomers

Iwamoto, T.; Kayahara, E.; Yasuda, N.; Suzuki, T.; Yamago, S. Angew. Chem. Int. Ed. 2014, 53, 6430-6434
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Macrocycles composed of aromatic subunits, like polycycloparaphenylenes, are of interest as components of nanotubes and for possible interesting optical properties. Tremendous advances have occurred over the past decade in preparing these rings ; see for examples these posts. Yamago now reports on the synthesis, optical properties and structure of [4]cyclo-2,7-pyrenylene 1, made by joining four pyrene units together.1

B3LYP/6-31G(d) optimization of the structure of 1 reveals a D2 geometry (Figure 1). This structure shows a very distorted pyrene unit. The strain energy of 1 is estimated as 392 kJ mol-1 (though how this was arrived at is not mentioned!), which is much larger than the strain energy of [8]-cycloparaphenylene.

Figure 1. B3LYP/6-31G(d) optimized structure of 1
This is another molecule to be sure to click on and rotate using JMol.

The nature of the HOMO and LUMO of 1 is very different than that of linear tetra-2,7-pyrene. The degenerate HOMOs and degenerate LUMOs of the linear compound have a node at the 2 and 7 positions and are localized to the terminal and central pyrene units, respectively. The HOMO and LUMO of 1 are fully delocalized. The implications of this are seen in the spectroscopy and electrochemistry of 1.


(1) Iwamoto, T.; Kayahara, E.; Yasuda, N.; Suzuki, T.; Yamago, S. "Synthesis, Characterization, and Properties of [4]Cyclo-2,7-pyrenylene: Effects of Cyclic Structure on the Electronic Properties of Pyrene Oligomers," Angew. Chem. Int. Ed. 201453, 6430-6434, DOI:


1: InChI=1S/C64H32/c1-2-34-18-50-20-36-4-3-35-19-49(17-33(1)57(35)58(34)36)51-21-37-5-7-41-25-53(26-42-8-6-38(22-51)59(37)61(41)42)55-29-45-13-15-47-31-56(32-48-16-14-46(30-55)63(45)64(47)48)54-27-43-11-9-39-23-52(50)24-40-10-12-44(28-54)62(43)60(39)40/h1-32H/b51-49-,52-50-,55-53-,56-54-

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

Monday, June 30, 2014

Band offsets of lattice-matched semiconductor heterojunctions through hybrid functionals and G0W0

Karim Steiner, Wei Chen, and Alfredo Pasquarello Phys. Rev. B 89, 205309 
Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

Fixing the well-known problems with DFT band gaps can be done, but there is no universal solution.  Hybrid functionals have long been used in the chemistry community (where B3LYP is very popular) while solid state approaches tend to use the PBE0 functional (which mixes 25% exact exchange into the PBE functional in its original form).  There are also many-body perturbation theory approaches such as GW, which tend to be complex but are generally very good at getting band structures correct.
The paper I’m writing about[1] considers how to calculate offsets between both the valence and conduction bands in different bulk semiconductors.  These have obvious applications in areas such as photovoltaics and electronics, and a consistent way to calculate them would be very valuable.  This isn’t the only paper which has been written about this problem; however, this paper is really very careful, and goes to some lengths to be open.  They give details of all the relevant parameters used in the calculations, even giving details of the pseudopotentials used.  This is unusual, and highly laudable: it makes the science they do open and reproducible.
They consider non-polar interfaces, and a range of group IV, III-V and II-VI semiconductors.  What do they find ? They use the common one-shot G0W0 method (which calculates corrections to DFT bands without any form of self-consistency) and find excellent agreement with experiment.  For PBE0, the hybrid functional, they find that the original mixing parameter of 0.25 gives poor results; they optimise it by fitting either to experimental gaps or to GW gaps.  When modelling the interfaces, they find that the experimental gap fitting process (taking an average of the values for the two semiconductors) is much better.  Interestingly, they find that the valence band offsets are modelled extremely well by PBE alone, though they have applied a similar method to semiconductor/insulator interfaces where PBE fared rather poorly[2].
Overall, the conclusions add to our understanding of hybrids and GW.  A single hybrid mixing parameter for all systems performs fairly poorly; however, fitting the parameter to (say) bulk properties can be successful.  GW is generally reliable, though it cannot calculate electron-hole effects, and can depend on input bands (in particular, small gap systems in DFT can form a poor starting point for GW calculations – see, for instance [3]).  But I really like their approach, constructing pseudopotentials very carefully and publishing their protocols.
[1] Phys. Rev. B 89, 205309 (2014) doi:10.1103/PhysRevB.89.205309
[2] Phys. Rev. Lett. 101, 106802 doi:10.1103/PhysRevLett.101.106802
[3] New J. Phys. 7 126 (2005) doi:10.1088/1367-2630/7/1/126

Wednesday, June 25, 2014

Total synthesis and isolation of citrinalin and cyclopiamine congeners

Mercado-Marin, E. V.; Garcia-Reynaga, P.; Romminger, S.; Pimenta, E. F.; Romney, D. K.; Lodewyk, M. W.; Williams, D. E.; Andersen, R. J.; Miller, S. J.; Tantillo, D. J.; Berlinck, R. G. S.; Sarpong, R. Nature 2014, 509, 318-324
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Here is another nice example of the partnership between experiment and computation in ascertaining molecular structure. The Sarpong, Tantillo and Miller groups collaborated on the synthesis, characterization and biosynthesis of some metabolites from Penniculium strains.1 I will focus here on just the structural identification component of this paper; the synthesis and the biosynthesis are very interesting too!

Cyclopiamine A 1 and cyclopiamine B 2 interconvert through an intermediate that allows for the epimerization at carbon bearing the nitro group.2


Citrinalin A 3 might also seem to undergo the same type of ring opening-ring closing reaction to produce citrinalin B. However, the original proposed structure3 of citrinalin B 4 implies an epimerization at a different carbon (at the ring fusion to the terminal 5 member ring). These authors suggested that perhaps the proper structure of citrinalin B is 5, which differs from citrinalin A only at the carbon bearing the nitro group, analogous to the relationship between 1 and 2.



The low energy conformations of both 4 and 5 (actually the trifluoroacetic acid salts) were optimized at B3LYP/6-31+G(d,p) and the chemical shifts for both 1H and 13C were computed, Boltzmann-weighted and scaled, and then compared with the NMR spectra of authentic citrinalin B. (The lowest energy conformations of 4 and 5 are shown in Figure 1.) The corrected mean absolute deviations for the 1H and13C chemical shift for the original structure 4 are 0.45 ppm and 2.0 ppm, respectively (with the largest outliers of 2.3 ppm for H and 9.6 ppm for C). These errors are about twice what is observed in comparing the experimental and computed 1H and 13C chemical shifts of 3. The agreement between the computed and experimental values using 5 are much improved, with mean deviations of 0.12 and 1.6ppm, and largest deviations of 0.38 ppm for 1H and 4.4 ppm for 13C. Use of Goodman’s DP4 method indicates a 100% probability that the structure of citrinalin B is 5. This prediction is confirmed by the x-ray structure.


Figure 1. B3LYP/6-31+G(d,p) optimized lowest energy conformers of 4 and 5.


(1) Mercado-Marin, E. V.; Garcia-Reynaga, P.; Romminger, S.; Pimenta, E. F.; Romney, D. K.; Lodewyk, M. W.; Williams, D. E.; Andersen, R. J.; Miller, S. J.; Tantillo, D. J.; Berlinck, R. G. S.; Sarpong, R. "Total synthesis and isolation of citrinalin and cyclopiamine congeners," Nature 2014509, 318-324, DOI:10.1038/nature13273.
(2) Bond, R. F.; Boeyens, J. C. A.; Holzapfel, C. W.; Steyn, P. S. "Cyclopiamines A and B, novel oxindole metabolites of Penicillium cyclopium westling," J. Chem. Soc., Perkin Trans I 1979, 1751-1761, DOI:10.1039/P19790001751.
(3) Pimenta, E. F.; Vita-Marques, A. M.; Tininis, A.; Seleghim, M. H. R.; Sette, L. D.; Veloso, K.; Ferreira, A. G.; Williams, D. E.; Patrick, B. O.; Dalisay, D. S.; Andersen, R. J.; Berlinck, R. G. S. "Use of Experimental Design for the Optimization of the Production of New Secondary Metabolites by Two Penicillium Species,"J. Nat. Prod. 201073, 1821-1832, DOI: 10.1021/np100470h.


1: InChI=1S/C26H33N3O5/c1-23(2)12-17(30)20-18(34-5)9-8-16-21(20)28(23)22(31)26(16)13-25(29(32)33)14-27-10-6-7-15(27)11-19(25)24(26,3)4/h8-9,15,19H,6-7,10-14H2,1-5H3/t15-,19+,25+,26-/m1/s1
2: InChI=1S/C26H33N3O5/c1-23(2)12-17(30)20-18(34-5)9-8-16-21(20)28(23)22(31)26(16)13-25(29(32)33)14-27-10-6-7-15(27)11-19(25)24(26,3)4/h8-9,15,19H,6-7,10-14H2,1-5H3/t15-,19+,25-,26-/m1/s1
3: InChI=1S/C25H31N3O5/c1-22(2)11-16(29)19-17(33-22)8-7-15-20(19)26-21(30)25(15)12-24(28(31)32)13-27-9-5-6-14(27)10-18(24)23(25,3)4/h7-8,14,18H,5-6,9-13H2,1-4H3,(H,26,30)/t14-,18+,24+,25-/m0/s1
4: InChI=1S/C25H31N3O5/c1-22(2)11-16(29)19-17(33-22)8-7-15-20(19)26-21(30)25(15)12-24(28(31)32)13-27-9-5-6-14(27)10-18(24)23(25,3)4/h7-8,14,18H,5-6,9-13H2,1-4H3,(H,26,30)/t14-,18-,24-,25+/m1/s1
5: InChI=1S/C25H31N3O5/c1-22(2)11-16(29)19-17(33-22)8-7-15-20(19)26-21(30)25(15)12-24(28(31)32)13-27-9-5-6-14(27)10-18(24)23(25,3)4/h7-8,14,18H,5-6,9-13H2,1-4H3,(H,26,30)/t14-,18+,24-,25-/m0/s1

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

Monday, June 23, 2014

Hydrogen Bonding Switches the Spin State of Diphenylcarbene from Triplet to Singlet

Costa, P.; Sander, W. Angew. Chem. Int. Ed. 2014, 53, 5122-5125
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Carbenes remain an active area of interest for computational chemists, as seen in Chapter 5 of my book. For many carbenes, the triplet is the ground state, and that is true of diphenylcarbene 1. Can solvent play a role in the stability of carbene spin states? Surprisingly, the answer, provided in a recent paper by Sander,1 is yes!

In the gas phase, the singlet-triplet gap of 1 is computed to be 5.62 kcal mol-1 at (U)B3LYP/6-311++G(d,p) (and this reduces to 5.06 at (U)B3LYP+D3/6-311++G(d,p)) with the ground state as a triplet. If a single methanol molecules is allowed to approach 1, then the complex involving the singlet has a short hydrogen bond distance of 1.97 Å but the triplet has a much longer distance of 2.33 Å. These structures are shown in Figure 1. This manifests in a dramatic change in the relative stability, with the singlet complex now 0.26 kcal mol-1 (0.44 with the D3 correction) lower in energy than the triplet.


Figure 1. (U)B3LYP/6-311++G(d,p) optimized geometries of the compelxes of methanol with singlet or triplet 1.

IR spectroscopy of 1 in an argon matrix doped with a small amount of methanol confirms the presence of the singlet carbene, and detailed description of the difference in the reactivities of the singlet and triplet are provided.


(1) Costa, P.; Sander, W. "Hydrogen Bonding Switches the Spin State of Diphenylcarbene from Triplet to Singlet," Angew. Chem. Int. Ed. 201453, 5122-5125, DOI: 10.1002/anie.201400176.


1: InChI=1S/C13H10/c1-3-7-12(8-4-1)11-13-9-5-2-6-10-13/h1-10H

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

Tuesday, June 17, 2014

Coherence penalty functional: A simple method for adding decoherence in Ehrenfest dynamics

A. V. Akimov, R. Long and O. V. Prezhdo, J. Chem. Phys. 140, 194107 (2014)
Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

As I discussed in my last blog post, probably the most common approach to non-adiabatic methods is called surface hopping.  This stochastic method has to explore configuration space, and relies on multiple trajectories.  There must be multiple electronic trajectories to sample the effects of the randomly determined jumps between electronic energy surfaces.  However, to sample the phase space properly, there should also be multiple MD trajectories, and potentially different starting energy surfaces.  This makes the cost rather large, and often it is the cost that dictates the choice of number of trajectories (though cost is often the determining factor in any ab initio simulations).

Ehrenfest dynamics, the approach where electronic degrees of freedom are propagated according to the time-dependent Schrodinger equation while nuclear degrees of freedom obey classical forces calculated from the mean-field potential energy surface, is simpler but does not allow for decoherence between energy surfaces.  The simplest way to understand this is to think about what electronic states would be found for each set of nuclear positions in a standard approach: these states define the different adiabatic surfaces.  Ehrenfest dynamics allows the electronic states to evolve in time, so they diverge from the adiabatic surfaces.  Physically, there are mechanisms that tend to collapse electrons onto specific surfaces: this is decoherence, and is missing in Ehrenfest dynamics.  Ehrenfest dynamics is simple: only a single trajectory needs to be followed; and the nuclear motion is determined by a mean-field potential energy.  However, this simplicity is offset by the lack of decoherence between states.

A recent paper[1] presents an alternative approach, designed to introduce decoherence into Ehrenfest dynamics.  But to understand the papers discussed in this post, we must first define some notation. We index the surfaces found above with labels $i, j$ etc.  Then if we represent the electronic structure in terms of these surfaces, the weight for each surface is $c_i, c_j$ etc.  (Note that these weights apply to the entire electronic system – the different many-body states.) Coherence in the simulation is often characterised in terms of the density matrix elements $c^*_i c_j$ – a measure of how the contributions of different surfaces interfere.  Ehrenfest is well-known to preserve coherence more than it should.
The authors rewrite the time-dependent Schrodinger equation (following an approach known as MMST in the chemistry literature) as classical Hamilton dynamics (the real part of the $c_i$ coefficients become the positions and the imaginary parts the momenta).  Having done this, they add a penalty functional to the Hamiltonian to increase the energy of system when there is coherence.  The form of the functional is not specified by they theory, and they choose a simple form: the sum over states of  the square of the matrix element defined above, each scaled by a multiplier.  These multipliers are related to the decoherence rates – which have to be calculated.  This is not a trivial problem, but the authors suggest a series of possible approximate routes to find the rates.

The tests that the authors apply, ranging from simple tests suggested by Tully in the context of the surface hopping method to large scale quantum dots on TiO$_2$, show significant improvement over Ehrenfest, and, for the large systems, sufficient accuracy to suggest that it produces reliable results.  In particular the times for electron transfers match well with experimental results.  This offers a new route to non-adiabatic calculations for large systems, which will be important for sufficiently accurate functionals.

In some ways the most interesting idea in the paper is the analogy between their penalty functional and DFT.  The authors note that Ehrenfest dynamics is a mean-field theory, which has problems correctly describing correlations between the nuclear motion and electronic states (this is the key source of excessive coherence in Ehrenfest).  Just as DFT lacks a correct description of electron correlations, Ehrenfest lacks a correct description of electron-nuclear correlations.  They suggest that their penalty functional might offer a similarly fruitful approach to Ehrenfest dynamics as DFT has been to electronic structure.  I think that the possibility to improve Ehrenfest dynamics in this way is very exciting.
This is not the first approach to augmenting Ehrenfest with decoherence rates.  The Schwartz group produced a method known as mean-field dynamics with stochastic decoherence[2].  They derive decoherence rates by attaching frozen Gaussians to nuclei, and at each Ehrenfest step calculate the probability for decoherence to a single energy surface.  The method requires one external parameter, the non-adiabatic interaction width which in turn determines the nuclear size. Subotnik presented a similar method to define decoherence rates[3], based on a partial Wigner transform.  The method is not simple, but shows that decoherence rates based on nuclear-electronic correlations are an excellent route to improving Ehrenfest dynamics.  He also shows that the approach gives a method to monitor the accuracy of an Ehrenfest calculation.

Non-adiabatic calculations are technically demanding, and there are many different methods available, normally from the quantum chemistry literature.  The recent movement to augmenting Ehrenfest or mean-field calculations is encouraging, and offers the hope of large-scale calculations on these important processes, in the same way that DFT has allowed large-scale calculations on electronic structure.

[1] A. V. Akimov, R. Long and O. V. Prezhdo, J. Chem. Phys. 140, 194107 (2014)10.1063/1.4875702
[2] M. J. Bedard-Hearn, R. E. Larsen and B. J. Schwartz, J. Chem. Phys. 123, 234106 (2005)10.1063/1.2131056
[3] J. E. Subotnik, J. Chem. Phys. 132, 134112 (2010) 10.1063/1.3314248

Tuesday, June 10, 2014

Gas-Phase Structure Determination of Dihydroxycarbene, One of the Smallest Stable Singlet Carbenes

Womack, C. C.; Crabtree, K. N.; McCaslin, L.; Martinez, O.; Field, R. W.; Stanton, J. F.; McCarthy, M. C.  Angew. Chem. Int. Ed. 2014, 53, 4089
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Dihdroxycarbene was the subject of a post a few years ago relating to how this carbene does not undergo tunneling,1 while related hydroxycarbene do undergo a tunneling rearrangement.

Now we have a gas-phase microwave determination of the trans,cis isomer of dihydroxycarbene.2 The computed CCSD(T)/cc-pCVQZ structure is shown in Figure 1. What is truly remarkable here is the amazing agreement between the experimental and computed structure – as seen in Table 1.The bond distance are in agreement within 0.001 Å and the bond angles agree within 0.3°! Just further evidence of the quality one can expect from high-level computations. And computing this structure was certainly far easier than the experiments!

Figure 1. CCSD(T)/cc-pCVQZ optimized geometry of dihydroxycarbene.

Table 1. Experimental and computed (CCSD(T)/cc-pCVQZ) geometric parameters of dihydroxycarbene.a


aDistances in Å and angles in deg.


(1) Schreiner, P. R.; Reisenauer, H. P. "Spectroscopic Identification of Dihydroxycarbene," Angew. Chem. Int. Ed. 200847, 7071-7074, DOI: 10.1002/anie.200802105.
(2) Womack, C. C.; Crabtree, K. N.; McCaslin, L.; Martinez, O.; Field, R. W.; Stanton, J. F.; McCarthy, M. C. "Gas-Phase Structure Determination of Dihydroxycarbene, One of the Smallest Stable Singlet Carbenes,"Angew. Chem. Int. Ed. 201453, 4089-4092, DOI: 10.1002/anie.201311082.


Dihydroxycarbene: InChI=1S/CH2O2/c2-1-3/h2-3H

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