NoteTaking_C6: Computational Organic Chemistry

CHAPTER 2: FUNDAMENTALS OF ORGANIC CHEMISTRY (2)

How to calculate ring strain energy (RSE)?

Take cyclobutane as an example, E(RSE)=E(cyclobutane)-E(CH₂CH₂CH₂CH₂).

Then what is E(CH₂CH₂CH₂CH₂)? 

E(CH₂CH₂CH₂CH₂)=4E(CH₂), where E(CH₂) is the energy of unstrained methylene group.

Using experimental gas-phase heats of formation:

E(CH₂)=$\Delta H_f(n{\rm{-hexane}})-\Delta H_f(n{\rm{-pentane}})$

Then E(RSE)=$\Delta H_f(\rm{cyclobutane})$-4E(CH$_2$)

Using quantum chemistry method:

"choose molecule(s) that can serve as reference(s) for the strained species"

(1) Isodesmic reaction:conserves the number of bonds of a given formal type.

(2) Homodesmotic reaction: conserves (1) the number of carbon (and other heavy) atoms and their state of hybridization and (2) the number of carbon (and other heavy) atoms bearing H atoms. 

(3) Group equivalent reaction: every equivalent group in the cyclic molecule must be paired with an equivalent group in a short acyclic molecule in the product.

Bachrach, S. M., “The Group Equivalent Reaction: An Improved Method for Determin- ing Ring Strain Energy,” J. Chem. Ed., 67, 907–908 (1990).

NoteTaking_C5: Computational Organic Chemistry

CHAPTER 2: FUNDAMENTALS OF ORGANIC CHEMISTRY (2)

1. Origin of the Acidity of Carboxylic Acids: Resonance vs Electrostatic/Inductive Effect

Electrostatic/Inductive Effect: Carbonyl group polarized the O-H bond, making the hydrogen more positively charged than the hydrogen of an alcohol.
Resonance Effect: The electrons and the nuclei reorganized to reproduce the ground state of anion.

Two theoretical studies supported the notion of importance of inductive effect:

(1). Rablen, P. R., “Is the Acetate Anion Stabilized by Resonance or Electrostatics? A Systematic Structural Comparison,” J. Am. Chem. Soc., 122, 357–368 (2000).

(2).  Holt, J. and Karty, J. M., “Origin of the Acidity Enhancement of Formic Acid over Methanol: Resonance versus Inductive Effects,” J. Am. Chem. Soc., 125, 2797–2803 (2003).

Both works proposed a clever way to partition the contribution of resonance energy and electrostatic energy.

2. RSE of Cyclopropane and Cyclobutane 

Ring strain energies of cycloporpane (27.6 kcal/mol) and cyclobutane (26.2 kcal/mol) are surprisingly similar.
The reason is not clear yet.

Various proposed contritions to RSE:
a. Baeyer strain (C-C-C)
b. Pitzer strain  (C-H)
c. Dunitz-Shomaker strain (1,3-CC repulsions)
d. re-hybridization (CH strengthening in cyclopropane)
e. $\sigma$-aromaticity (not well defined)

SCF Loop in Quantum Espresso

Summarized from Shobhana Narasimhan's presentation: PWscf: The SCF Loop and Some Relevant Input Parameters

Goal: Solve the Kohn-Sham equations:

$\left[\frac{1}{2}\nabla ^2 + V_{nuc}(r)+V_{H}[{n(r)}]+V_{xc}[n(r)] \right] \psi_i(r) = \varepsilon_i\psi_i(r) $

Since Hamiltonian, $H$, depends on solution $n(r)$, a self-consistent solution is required.

                                     Self-consistent Iterative Solution

=======================================================================
Q: Why using plane waves for periodic systems?
The most naive answer: as plane wave is periodic, it is nature to choose plane waves as the basis sets for periodic systems.

The wavefunction is a linear composition of plain waves,
$\psi_k(r)=\frac{1}{V}\sum_G c_{k,G}e^{i(k+G)r}$, where G is reciprocal lattice vector, k is the k points in BZ.
We can truncate the expansion at some value of |k+G|, as the contribution from higher Fourier components (large |k+G|) is usually small.
Related parameter: ecutwfc

$\frac{\hbar^2|k+G|^2}{2m}\le E_{cut}$

=======================================================================

0. Define the System

Related parameter: ibrav, celldm(i), nat, ntyp, 
Field: Atomic_Positions

1. Obtaining $V_{nuc}$

Given that (1) core wavefunctions sharply peaked near nucleus and (2) valence wavefunctions have lots of wiggles near nucleus, plane-wave basis sets with very high energy cutoff (|k+G|) are required to correctly describe the total wavefunctions. The problem is solved (computational expense is reduced) via using pseudopotentials. 
Related parameter: ecutwfc, ecutrho
Note:

"ecutwfc and ecutrho are the kinetic energy cutoff of the plane wave basis used to describe the wavefunctions and the charge density, respectively. The choice of the cutoffs depends only on the pseudopotentials (PPs) you are using and on which quantities you are interested in (for instance, magnetic properties are harder to converge wrt total energy). Once you've chosen your PPs you can test convergence of the energy wrt to cutoffs in the cheapest system (usually in the bulk). The threshold of convergence is somehow your choice: personally I don't go below 1mRy since this is the expected order of magnitude of the errors due to PPs."
For normconserving PPs you don't need to specify ecutrho, since the default (4*ecutwf) is good. For ultrasoft PPs you will need higher ecutrho due to the augmentation terms in the charge. "

2. Initial Guess for $n(r)$

Related parameter: startingwfc
a. superpositions of atomic densities, atomic;
b. converged n(r) from previous calculation, file;
c. random numbers, random.

3. $V_H$ and $V_{xc}$

Exchange-Correlation functions
LDA: pz
GGA: pw91, pbe

4. Potential & Hamiltonian

Kohn-Sham equations in plane wave basis

5. Diagonalization

Exact diagonaization can give $N_{PW}$ eigenvalues and eigenfunctions,  but the process is highly expensive as 
a. $N_{PW}$ is usually large
b. $N_{PW}^2$ elements to sore 
c. $N_{PW}^3$ operations are required to do the diagonalization
d. We usually only care limited number of bands 

Iterative diagonalization: recast diagonalization as a minimization process
Conjugate Gradient and Davidson
Related parameter: diagonalization, nbnd

6. New charge density: $n(r)=\sum_i|\psi_i(r)|^2$

Sum over all the k points in Brillouin Zone

k-point meshes:
1. Special points [Chadi & Cohen]
2. Monkhorst-Pack

Field: K_POINTS
Related parameter : nk1, nk2, nk2, k1, k2, k3

7. Test for Convergence 

a. compare new and old wavefunctions, charge densities, and total energies
b. compare with energy estimated with Harris-Foulkes functional
Related parameter: conv_thr (1.0e-8), electron_maxstep

8. Mixing to generate new $n(r)$

Take some combinations of input and output densities
 Related parameter: mixing_mode, mixing_beta

9. Output Quantities

$E_{tot}$: sum over eigenvalues, correct for double-counting of Hartree & XC terms, add ion-ion interactions
Force: Hellmann-Feynman theoreom 

Properties of BiFeO3

   BiFeO$_3$ has draw extensive studies in past years. I will try to summarize some of its properties in this post. This post will keep update.

A nice review: 

Physics and Applications of Bismuth Ferrite,  Gustau Catalan and James F. Scott, Advanced Materials, 2009, 21, 2463-2485. 

1. Crystal Structure

Room-temperature phase of BiFeO$_3$ is classed as rhombohedral ($R3C$).

The perovskite-type unit cell has a lattice parameter, $a_{\rm rh}$, of 3.965  and a rhombohedral angle, $\alpha_{\rm rh}$, of ca. 89.3–89.48 at room temperature with ferroelectric polarization along [111]. The unit cell can also be described in a hexagonal frame of reference, with the hexagonal $c$-axis parallel to the diagonals of the perovskite cube. 

The Goldschmid tolerance factor $t=(r_{\rm Bi} +r_{\rm O})/ \sqrt{2} (r_{\rm Fe}+r_{\rm O})=0.88$, which is much smaller than 1. Therefore, the oxgen cage has to rotate so as to fit into a cell that is too small. The oxgen cage is rotated around the polar [111] axis by 11-14$^\circ$, with the directly related Fe-O-Fe angle 154-156$^\circ$. 

2. Temperature-induced phase transition for bulk BiFeO$_3$  ($\alpha \rightarrow \beta$ transition)

$T_c$ = 825 ˚C

$R3C$ $\rightarrow$  $Pbnm$

$T_c$ = 931 ˚C

$Pbnm$ $\rightarrow$  $Pm\bar{3}m$

Atomic displacements in BiFeO3 as a function of temperature: neutron diffraction study,

A. Palewicz, R. Przeniosło,I. Sosnowska and A. W. Hewat, Acta Cryst.  B63, 537–544 (2007).

Ferroelectric-Paraelectric Transition in BiFeO3: Crystal Structure of the Orthorhombic Phase, Donna C. Arnold, Kevin S. Knight, Finlay D. Morrison, and Philip Lightfoot, Phys. Rev. Lett. 102, 027602 (2009). 

BiFeO3 Crystal Structure at Low Temperatures, A. Palewicz, I. Sosnowska, R. Przeniosło and A.W. Hewat, Acta Physica Polonica A, 117, 296 (2010).

3. Pressure-induced phase transition 

Structural Properties of Multiferroic BiFeO3 under Hydrostatic Pressure, Alexei, A. Belik,
Hitoshi Yusa,  Naohisa Hirao,  Yasuo Ohishi,  and Eiji Takayama-Muromachi, Chem. Mater. 2009, 21, 3400–3405

They found three new orthorhombic phases of BiFeO3 at room temperature below 9.7 GPa. Two transitions occurs at ~4 GPa and ~7 GPa, respectively.  

$R3C$ $\rightarrow$ OI $\rightarrow$  OII $\rightarrow$ OII

Effect of high pressure on multiferroic BiFeO3, R. Haumont, P. Bouvier, A. Pashkin, K. Rabia, S. Frank, B. Dkhil, W. A. Crichton, C. A. Kuntscher, and J. Kreisel, Phys. Rev. B, 79, 184110 (2009)

Multiple high-pressure phase transitions in BiFeO3, Mael Guennou, Pierre Bouvier, Grace S. Chen, Brahim Dkhil, Raphae Haumont, Gaston Garbarino, and Jens Kreise, Phys. Rev. B,  84, 174107 (2011)

They found six phase transitions in the 0-60-GPa range. At low pressure, four transitions occurred at 4, 5, 7, and 11 GPa.  The nonpolar $Pnma$ phase remains stable over a large pressure range between 11 and 38 GPa. Two high-pressure phase transitions at 38 and 48 GPa are marked by the occurrence of larger unit cells and increase of the distortion away from the cubic parent perovskite cell.

4. Thin film BiFeO$_3$

A mechanism for the 150 μC cm$^{−2}$ polarization of BiFeO3 films based on first-principles calculations and new structural data, Dan Ricinschi, Kwi-Young Yun and Masanori Okuyama,  J. Phys.: Condens. Matter 18 (2006) L97–L105

A Strain-Driven Morphotropic Phase Boundary in BiFeO3, R. J. Zeches et. al, Science 326, 977 (2009);

Stress-induced R-MA-MC-T symmetry changes in BiFeO3 films, H. M. Christen, J. H. Nam, H. S. Kim, A. J. Hatt, and N. A. Spaldin, Phys. Rev. B 83, 144107 (2011)

With compressive strain the structural proression is rhombohedral $\rightarrow$ (R-like)  monoclinic $\rightarrow$ (T-like) monoclinic $\rightarrow$ tetragonal. When the compressive strain exceeds 4.5%,  $c/a$ ratio is near 1.25.

A phase transition close to room temperature in BiFeO3 thin films, J Kreisel, P Jadhav, O Chaix-Pluchery, M Varela, N Dix, F Sanchez and J Fontcuberta, J. Phys.: Condens. Matter 23 (2011) 342202

Multiferroic phase transition near room temperature in BiFeO3 films,  I.C. Infante et, al.
Recent experiments reveals that the T-like phase exhibits both structural, magnetic and FE phase transitions in a narrow temperature range around 360 K

5. First-principle studies

First-principles study of spontaneous polarization in multiferroic BiFeO3,

J. B. Neaton,1, C. Ederer, U. V. Waghmare, N. A. Spaldin, and K. M. Rabe, Phys. Rev. B,  71, 014113 (2005)

Theoretical investigation of magnetoelectric behavior in BiFeO3, P. Ravindran, R. Vidya, A. Kjekshus, and H. Fjellvåg, Phys. Rev. B 74, 224412 (2006)

First-principles study of ferroelectric domain walls in multiferroic bismuth ferrite

Axel Lubk, S. Gemming,  and N. A. Spaldin, Phys. Rev. B 80, 104110 (2009)

Phase Transitions in Epitaxial (-110) BiFeO3 Films from First Principles,  S. Prosandeev, Igor A. Kornev, and L. Bellaiche, Phys. Rev. Lett. 107, 117602 (2011)

6. Classical Force Field

Development of a bond-valence based interatomic potential for BiFeO3 for accurate molecular dynamics simulations, S. Liu, I. Grinberg, and A. M. Rappe, J. Phys.: Condens. Matter 25 102202 (2013)

NoteTaking_B3: Crystal & Crystal Structures

CH2 TWO-DIMENSIONAL PATTERNS 

A point group describes the symmetry of an isolated shape

A plane group describes the symmetry of a two-dimensional repeating pattern. 

1. Rotation axes

monad, diad, triad,  tetrad, pentad,  hexad,  center of symmetry (inversion center)

Only rotation axes that can occur in a lattice are 1, 2, 3, 4, 6. 

The five-fold rotation axis and all axes with $n$ higher than 6, cannot occur in a plain lattice.

A unit cell can not show overall five-fold rotation symmetry does not mean that a unit cell cannot obtain a pentagonal arrangement of atoms. 

2. Five plane lattices

Oblique (mp)

Rectangular primitive (op)

Rectangular centered (oc)

Square (tp)

Hexagonal primitive (hp)

3. Ten plane crystallographic point symmetry groups

rotation axes + mirror planes

1, 2, 3, 4, 6, $m$, 2$mm$, 3$m$, 4$mm$, 6$mm$

See figures here.

4. Symmetry of patterns: the 17 plane groups (Wallpaper group)

1) 17 plane groups = 5 plane lattices + 10 point symmetry groups

2) glide line = translation + reflection

     The glide displacement is restricted to one half of the lattice repeat. 

3) Name of plane groups

Taken from Wikipedia:

For wallpaper groups the full notation begins with either p or c, for a primitive cell or a face-centred cell.  This is followed by a digit, n, indicating the highest order of rotational symmetry: 1-fold (none), 2-fold, 3-fold, 4-fold, or 6-fold. The next two symbols indicate symmetries relative to one translation axis of the pattern, referred to as the "main" one; if there is a mirror perpendicular to a translation axis we choose that axis as the main one (or if there are two, one of them). The symbols are either m, g, or 1, for mirror, glide reflection, or none. The axis of the mirror or glide reflection is perpendicular to the main axis for the first letter, and either parallel or tilted 180°/n (when n > 2) for the second letter. Many groups include other symmetries implied by the given ones. The short notation drops digits or an m that can be deduced, so long as that leaves no confusion with another group.

The first letter gives the lattice type. The number in the second place gives the rotation axis

Mirrors and glide lines are placed last. 

5. General and special positions

general position: position does not fall on a symmetry element, high multiplicity, low site symmetry

special position: position falls on a symmetry element, low multiplicity, high site symmetry

NoteTaking_C4: Computational Organic Chemistry

CHAPTER 2: FUNDAMENTALS OF ORGANIC CHEMISTRY (1)

1. Bond Dissociation Enthalpy (BDE)

AB → A⋅ + B⋅
BDE=E(A⋅) + E(B⋅) - E(A-B)

HF significantly underestimated the BDE.
MP2 is better by taking into account some correlation energy.
B3LYP can work quite well for some types of bond cleavages, but fails completely for some molecules. B3LYP may be acceptable for small molecules but its error increases with the size of the molecule. Be cautious when using B3LYP!
The composite methods are the best and are always chosen as the benchmark method. But the computational  expense becomes unaffordable for large molecules.

Case 1: BDF of R-X with different α-substitution (J. Phys. Chem. A, 2005, 109, 7558)
BDE of the C-H bond and C-C decreases with increasing α-substitution. The opposite trend is observed for the C-F bond and C-OH bonds. And there is a turn over for C-OMe (first increases then decreases). B3LYP predicted the correct trend of decreasing C-H and C-C BDEs with increasing α-substitution, and increasing C-F BDEs with increasing α-substitution. But B3LYP fails to predict the trend of C-OH, and the turn over in C-OMe is too easy.

α-substitution increases the electron-donating capability of R.   When α-substitution increases, the R-X bond becomes more polar. R-X can be considered as a highbrid of three resonant structures: R:X ↔ R^- X⁺ ↔ R⁺X^-.  When X=F and OH, the ionic structure R⁺X^- becomes very important, and a more polared R-X bond will have high strength. When X=H or C, the stability of alkyl radical dominates the trend in BDEs.  When X=O, the two trends work in opposite directions.

Case2: BDE of Oximes

R₁R₂C=NOH → R₁R₂C=NO⋅+ H⋅

DFT and CBS-QB3 help to calirfy the controversy of the BDE of oximes.

2. Acidity 

AH → A^- + H⁺

DPE: DeProtonation Energy
The effect of diffuse functions is very important to get accruate DPE.
The role of diffuse functions is to preperly describe the tail behavior of orbitals, especially in electron rich systems where electron-electron repulsion leads to a more extensive distribution than usual.
MP2/6-31+G(d) is a reasonable accurate method.

B3LYP:

" For a set of 45 acids, the average absolute error for the B3LYP/6-31+G(d)-predicted DPEs is 4.7 kcal/mol. Improvement of the basis set leads to smaller errors. With the B3LYP/6-311+G(2d,p) method, even the DPE of acetic acid is predicted with an error less than 2.5 kcal/mol. For the set of 45 acids, the average error is only 2.1 kcal mol21 with the B3LYP/6-311++G(3df,2pd) level. "

B3PW91 is also a good functional. 

NoteTaking_C3: Computational Organic Chemistry

CHAPTER 1: QUANTUM MECHANICS FOR ORGANIC CHEMISTRY (3)

Density Functional Theory

1) Hohenberg-Kohn existence theorem: The exact electronic energy is a unique functional of electron density.

The electron density obeys the variational theorem.

2) Functional relates a function to a scalar quantity.

3) Kohn-Sham equation

A key point in Kohn-Sham equation is that the kinetic energy is defined as the kinetic energy of noninteracting electrons whose density is the same as the density of the real interacting electrons. The exchange-correlation functional takes care of all the other aspects of a true system.
The Kohn-Sham procedure is then to solve for the orbitals that minimize the energy.

The Kohn–Sham orbitals are separable by definition (the electrons they describe are noninteracting), analogous to the HF MOs. Therefore, the Kohn-Sham equation can be solved in a self-consistent way using similar steps as was used in the Hatree-Fork-Roothann method. 

Therefore, for a similar cost as the HF methods, DFT reproduces the energy of a molecule that includes the electron correlation.  But, there is no guidance as to the form of the density functional. 

Paraphrasing Cramer’s description of the contrast between HF and DFT, HF and the various post-HF electron correlation methods provide an exact solution to an approximate theory, but DFT provides an exact theory with an approximate solution.

4)  Exchange-Correlation Functionals

                                            E_{xc}=E_x + E_c

Local density Approximation (LDA) and Generalized Gradient Approximation (GGA)

Exchange functional:

"B" Becke exchange 

Correlational functional:

"LYP" Lee, Yang and Parr, second derivative of density

"PW" Perdew and Wang , first derivative of density

Hybrid functional: add some HF exchange