of 8
Supporting Information for
Chen
,
Cusumano,
Flesch
, Strong, Goddard,
and Stoltz
S
1
Supporting Information for
Molecular Dynamics Investigations of Dienolate [4+2]
Reactions
Peng
-
Jui Chen
a,
, Alexander Q. Cusumano
a,
,
Kaylin N. Flesch
a
, Christian Santiago Strong
a
,
William A. Goddard III
b,
*, Brian M. Stoltz
a,
*
a
The Warren and
Katharine Schlinger Laboratory for Chemistry and Chemical Engineering,
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena,
California 91125, United States
b
Materials and Process Simulation Center, Beckman Institute, California Institute of
Technology, Pasadena, California 91125, United States
These authors contributed equally.
stoltz@caltech.edu
wag@caltech.edu
Table of Contents:
General Notes
Quantum Mechanics Calculations
................................
................................
................................
.....
2
General Notes
Molecular Dynamics Simulations
................................
................................
................................
......
3
Configuration and Parameters for Progdyn
................................
................................
................................
..................
4
Comparison of Computational Results from Different Functionals
................................
................................
......
6
References
................................
................................
................................
................................
................................
................
7
Supporting Information for
Chen
,
Cusumano,
Flesch
, Strong, Goddard,
and Stoltz
S
2
General Notes
Quantum Mechanics Calculations
All quantum mechanic
s
calculations were carried out with the ORCA program.
1
Geometry
optimizations
,
harmonic frequency calculations
, and single
-
point energy evaluations
were carried
out with density functional theory (DFT).
While PBE0
-
D4, the PBE0 hybrid functional
2
paired
with Becke
Johnson damped D4 dispersion corrections
3
, proved to be a robust method for Pd
-
catalyzed reactions in our prior studies,
4
generalized gradient approximation (GGA) functional
PBE paired with D4 corrections was utilized in this study as a compromise between robustness
and computational efficiency due to the size of the systems of interest.
For geometry optimization
and harmonic frequency calculations
,
Pd
is described by the def2
-
TZVP all electron basis se
t
5
and
the ECP28MWB small
-
core (
18
explicit valence electrons) quasi
-
relativistic pseudopotential
,
6
while
C, H, N, and P
are assigned the def2
-
SV
(
P
)
basis. Diffuse functions are added to
oxyge
n
(ma
-
def2
-
SVP
)
.
Herein, we refer to this composite basis set as BS1.
Geometry o
ptimization and
harmonic
frequency calculations were carried out with the CPCM implicit solvation model for
respective solvents and their corresponding dielectric constants
.
For all calculations employing
CPCM, surface charges are described by the improved Gaussian charge scheme of Neese and
coworkers with a scaled Van der Waals cavity (
a
= 1.2).
7
All Hessians were computed analytically.
Stationary points are characterized by the correct number of imaginary vibrational modes (zero for
minima and one for saddle points). Cartesian coordinates of all optimized structures are included
as “.xyz” files
,
available online in a compressed zip file format.
Final Gibbs free energies were obtained by applying thermodynamic corrections to
the
calculated
electronic energies. Thermodynamic corrections from harmonic frequency calculations
employ the quasi
-
ridged rotor harmonic oscillator approach to correct for the breakdown of the
harmonic oscillator approximation at low vibrational frequencies.
8
Note that free energies are
adjusted to a 1 M standard state. The translational (S
trans
) and rotational entropy (S
rot
) contributions
to the Gibbs free energy calculated for a complex in condensed phase are
ca.
40
60% of the values
obtained assuming an ideal gas.
9
As suggested in the literature, S
trans
and S
rot
obtained by ideal gas
treatment are scaled by a factor of 0.5 to obtain the final condensed phase values.
10
Hence, the
final Gibbs free energy at
333
.15 K is calculated as:
!"#$
=
&#
,
!"#$
()
*
+
푍푃퐸
+
+,-.!
+
,"+
+
$/0
+
0
*
&#
+
$/0
+
1
2
+,-.!
+
1
2
,"+
.
+
1
Supporting Information for
Chen
,
Cusumano,
Flesch
, Strong, Goddard,
and Stoltz
S
3
The resolution of identity (RI) and Chain
-
of
-
Spheres (COS) approximations are employed
for efficient evaluation of Coulomb and exchange integrals, respectively.
11
The def2/J auxiliary
basis
12
is employed for all atoms except
oxygen
, for which a suitable auxiliary was obtained via
the automatic generation algorithm in the ORCA program (keyword:
AutoAux
).
13
Very fine grid
settings are employed in all calculations (optimization/frequency calculations:
DefGrid2
).
Transition states of the intramolecular Li dienolate [4+2] reaction and conjugated Pd
enolate [4+2] reaction were computed through a relaxed surface scan and subsequent geometry
optimization to a saddle point with a single imaginary frequency (keyword:
optTS
). The observed
bond distances in these transition states were then exploited in the transition state search for the
intermolecular Li dienolate and siloxy diene [4+2] reactions.
T
he reactive bond distances in these
two systems were constrained to the bo
nd distances of the intramolecular Li dienolate [4+2]
transition state
,
and the resulting structures
were
directly
optimized to
the respective transition
states.
Conformer searching was carried out for
the
stationary point
s of siloxy diene and
conjugated Pd enolate reactions
using the meta
-
dynamics
-
based CREST program
(using GNF
-
FF)
from the Grimme group.
14
Duplicate conformers were removed, low energy conformers were
subsequently
optimized
,
and energies
were
evaluated
at
the level of theory
mentioned prior.
The
lowest energy conformer of the two reactions were subsequently used as starting points in
molecular dynamics simulations.
General Notes
Molecular Dynamics Simulations
Quasiclassical trajectories were initialized by normal mode sampling at the saddle points.
Starting geometries
were derived
that correspond to
a
Boltzmann
distribution (corresponding to
thermal energy available at the respective reaction temperatures)
of geometries
along each normal
mode
.
The result is a weighted distribution of starting structures along the dividing surface at the
saddle point.
Normal modes are given vibrational zero
-
point energy and a randomized phase.
Trajectories were propagated in the
forward and reverse directions with a velocity
-
Verlet
algorithm in 1.0 fs timesteps, as implemented in Sing
l
eton’s Progdyn program
.
15
Trajectories were
propagated
until either formation of cycloadduct (forming both C
C bond lengths < 1.5 Å),
formation of
reactant
(the shorter C
C bond length > 3 Å and the longer > 4Å), or
of a maximum
simulation time of 500 fs was reached
.
Supporting Information for
Chen
,
Cusumano,
Flesch
, Strong, Goddard,
and Stoltz
S
4
Considering the size of the system,
100 randomized
trajectories were sampled to maintain
tractability, and we
employed the efficient
generalized gradient approximation (GGA) functional
PBE for the
QM energies and forces
. To investigate the influence of DFT functional choice on
simulation results, a suite of control experiments was performed on a model reaction using hybrid,
meta
-
hybrid, and range separated functionals. While simulations with different functionals yield
s
lightly variable results, all simulat
ion results described in this report were performed with the
same functionals for comparison
.
Most importantly,
the same qualitative results are obtained
between polarized and non
-
polarized systems
(vide infra)
.
The initial configurations for Progdyn
are described below and samples of the remaining
supporting files are provided online.
Configuration and
Parameters for Progdyn
A sample
progdyn.conf
file for the Pd enolate [4+2] cycloaddition is provided below for
reference.
Supporting Information for
Chen
,
Cusumano,
Flesch
, Strong, Goddard,
and Stoltz
S
5
Note that the value of
etolerance
must be increased for the large conformationally flexible
system. A value of 10 still leads to several initial trajectories being rejected. This value is adjusted
accordingly for each system. The
methodfile
contains additional input information, such as basis
set assignments. For the system above as an example:
method PBE
basis def2
-
SV(P)
charge 0
multiplicity 1
processors 16
memory 2000
diagnostics 3
title G
initialdis 2
timestep 1E
-
15
scaling 1.0
temperature 33
3
.15
method2 D4
method3 Def2/J
method4 CPCM(Toluene)
methodfile 6
numimag 1
searchdir positive
classical 0
keepevery 1
highlevel 999
boxon 0
boxsize 7.5
etolerance 10
damping 1
reversetraj true
#JAN 2020 1.0 version, based on 2010 version progdyn
Supporting Information for
Chen
,
Cusumano,
Flesch
, Strong, Goddard,
and Stoltz
S
6
The value of 6 for the
methodfile
input in the
progdyn.conf
file corresponds to the six lines
in the
methodfile
file.
Comparison of
Computational
Results from
D
ifferent Functionals
Limited by the size of the system of interest,
generalized gradient approximation (GGA)
functional
PBE
-
D4 was utilized
to maintain computational efficiency
.
To investigate the effect of
the lack of terms describing Hartree
Fock exchange energy
, a suite of control experiments
, each
consisted of 100 trajectories initiated from the transition state for molecular dynamics simulation,
w
as
performed on
non
-
polarized and polarized
model system
s
(
Figure S1 and S2
)
using different
functionals, including hybrid functionals PBE0,
2
B3LYP,
16
and M06
-
2X.
17
In addition, range
-
separated functional
w
B97x
-
D3
was also sampled
(keyword
in ORCA
:
wB97X
-
D3BJ
)
.
18
In this
study, PBE0 and B3LYP were paired with
Becke
Johnson damped D4 dispersion corrections
.
For
comparison, calculations in vacuum at different temperatures were also conducted
with
the PBE
-
D4
functional
.
The
computed transition structures, PES features, and molecular dynamics
simulation results
are
summarized in Figure
s S1
-
S2
.
It was observed that
the absence of the
Hartree
Fock exchange term generally led to more asynchronous transition structures and longer
average ∆t
. In addition,
flatter potential energy surfaces with imaginary vibrational frequencies of
smaller magnitudes
at the transition state w
ere
observed in the absence
o
f Hartree
Fock exchange
term
.
19
Performing the calculations in vacuum led to slightly
less flat PES and
shorter ∆t, but the
effect was not significant.
%basis
NewGTO Pd "def2
-
TZVP" end
NewECP Pd "def2
-
ECP" end
NewGTO O "ma
-
def2
-
SV(P)" end
NewAuxJGTO O "AutoAux" end
end
Supporting Information for
Chen
,
Cusumano,
Flesch
, Strong, Goddard,
and Stoltz
S
7
Figure S
1
.
Comparison of transition struct
u
res, PES features, and molecular dynamics simulation results
of non
-
polarized [4+2] reaction
from different functionals.
a
Imaginary vibrational frequencies were multiplied by the
corresponding scaling factors at the level of theory utilized.
20
b
Bond distance cutoff of 1.6 Å used for fully
formed C
C bond.
Figure S
2
.
Comparison of transition structures, PES features, and molecular dynamics simulation results
of
polarized [4+2] reaction
from different functionals.
a
Imaginary vibrational frequencies were multiplied by the
corresponding scaling factors at the level of theory utilized.
20
b
Bond distance cutoff of 1.6 Å used for fully
formed C
C bond.
References
1
(a) Neese, F. Software Update: The ORCA Program System, Version 4.0.
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2018,
8
, No. e1327. (b) Neese, F. The ORCA Program System.
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2012
,
2
,
73−78.
2
Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0
Model.
J. Chem. Phys.
1999
,
110
, 6158−6170.
3
(a) Caldeweyher, E.; Ehlert, S.; Hansen, A.; Neugebauer, H.; Spicher, S.; Bannwarth, C.; Grimme, S. A Generally
Applicable Atomic
-
Charge Dependent London Dispersion Correction.
J. Chem. Phys.
2019
,
150
(15), 154122. (b)
+
1
2
3
TS1
Method
Solvent
Temperature (K)
d
1
(Å)
d
2
(Å)
Imaginary
Vibrational
Frequency (cm
–1
)
a
Average
t (fs)
b
Median
t (fs)
Stepwise
Trajectories (%)
Concerted
Trajectories (%)
333.15
298.15
333.15
333.15
333.15
333.15
PhMe
gas
gas
PhMe
PhMe
PhMe
PBE-D4
PBE-D4
PBE-D4
PBE0-D4
B3LYP-D4
M06-2X
2.35
2.35
2.35
2.31
2.28
2.27
2.35
2.35
2.35
2.31
2.28
2.27
–410.52
–408.58
–408.58
–477.82
–518.23
–497.06
4.5
4.1
4.8
4.0
3.7
2.8
4
3
4
3
4
3
0
0
0
0
0
0
100
100
100
100
100
100
MeO
CN
MeO
CN
MeO
CN
+
d
2
d
1
S1
S2
S3
TS-S1
Method
Solvent
Temperature (K)
d
1
(Å)
d
2
(Å)
Imaginary
Vibrational
Frequency (cm
–1
)
a
Average
t (fs)
b
Median
t (fs)
Stepwise
Trajectories (%)
Concerted
Trajectories (%)
333.15
298.15
333.15
333.15
333.15
333.15
333.15
PhMe
gas
gas
PhMe
PhMe
PhMe
PhMe
PBE-D4
PBE-D4
PBE-D4
PBE0-D4
B3LYP-D4
M06-2X
ω
B97x-D3
2.06
2.06
2.06
2.03
1.99
2.04
2.00
2.90
2.85
2.85
2.74
2.72
2.50
2.59
–354.86
–369.77
–369.77
–422.32
–442.18
–543.29
–544.25
45.9
40.6
37.0
25.3
26.2
10.2
14.1
38.0
36.0
34.0
24.0
23.0
10.0
12.5
19
14
11
3
1
0
0
81
86
89
97
99
100
100
Supporting Information for
Chen
,
Cusumano,
Flesch
, Strong, Goddard,
and Stoltz
S
8
Caldeweyher, E.; Bannwarth, C.; Grimme, S. Extension of the D3 Dispersion Coefficient Model.
J. Chem. Phys.
2017
,
147
, 034112.
4
(A) Cusumano, A. Q.; Goddard, W. A. I.; Stoltz, B. M. The Transition Metal Catalyzed [π2s + π2s +
s
2s +
s
2s]
Pericyclic Reaction: Woodward
Hoffmann Rules, Aromaticity, and Electron Flow.
J. Am. Chem. Soc.
2020
,
142
,
19033
19039.
(B) Cusumano, A. Q.; Stoltz, B. M.; Goddard, W. A. Reaction Mechanism, Origins of
Enantioselectivity, and Reactivity Trends in Asymmetric Allylic Alkylation: A Comprehensive Quantum Mechanics
Investigation of a C(sp
3
)
C(sp
3
) Cross
-
Coupling.
J. Am. Chem. Soc.
2020
,
142
, 13917
13933.
5
Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence
Quality for H to Rn: Design and Assessment of Accuracy.
Phys. Chem. Chem. Phys.
2005
,
7
, 3297
3305.
6
Peterson, K. A.; Figgen, D.; Goll, E.; Stoll, H.; Dolg, M. Systematically Convergent Basis Sets with Relativistic
Pseudopotentials. II. Small
-
Core Pseudopotentials and Correlation Consistent Basis Sets for the Post
-
d
Group 16
18
Elements.
J
.
Chem
.
Phys
.
2003
,
119
, 11113
11123.
7
Garcia‐Ratés, M.; Neese, F. Effect of the Solute Cavity on the Solvation Energy and Its Derivatives within the
Framework of the Gaussian Charge Scheme.
J
.
Comput
.
Chem
.
2020
,
41
, 922
939.
8
Grimme, S. Supramolecular Binding Thermodynamics by Dispersion
-
Corrected Density Functional Theory.
Chem.
Eur. J.
2012
,
18
, 9955
9964.
9
Izato, Y.; Matsugi, A.; Koshi, M.; Miyake, A. A Simple Heuristic Approach to Estimate the Thermochemistry of
Condensed
-
Phase Molecules Based on the Polarizable Continuum Model.
Phys. Chem. Chem. Phys.
2019
,
21
, 18920
18929.
10
Finkelstein, A. V.; Janin, J. The Price of Lost Freedom: Entropy of Bimolecular Complex Formation.
Protein
Engineering, Design and Selection
1989
,
3
, 1
3.
11
Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Efficient, Approximate and Parallel Hartree
-
Fock and Hybrid
DFT Calculations. A ‘Chain
-
of
-
Spheres’ Algorithm for the Hartree
-
Fock Exchange.
Chemical Physics
2009
,
356
, 98
109.
12
Weigend, F. Accurate Coulomb
-
Fitting Basis Sets for H to Rn.
Phys. Chem. Chem. Phys.
2006
,
8
, 1057
1065.
13
Stoychev, G. L.; Auer, A. A.; Neese, F. Automatic Generation of Auxiliary Basis Sets.
J. Chem. Theory Comput.
2017
,
13
, 554
562.
14
Grimme, S. Exploration of Chemical Compound, Conformer, and Reaction Space with Meta
-
Dynamics Simulations
Based on Tight
-
Binding Quantum Chemical Calculations.
J. Chem. Theory Comput.
2019
,
15
, 2847
2862.
15
Aziz, H. R.; Singleton, D. A. Concert along the Edge: Dynamics and the Nature of the Border between General and
Specific Acid
Base Catalysis.
J. Am. Chem. Soc.
2017
,
139
, 5965
5972.
16
Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and
Circular Dichroism Spectra Using Density Functional Force Fields.
J. Phys. Chem.
1994
,
98
, 11623
11627.
17
Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical
Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic
Testing of Four M06
-
Class Functiona
ls and 12 Other Functionals.
Theor Chem Account
2008
,
120
, 215
241.
18
Najibi, A.; Goerigk, L. The Nonlocal Kernel in van Der Waals Density Functionals as an Additive Correction: An
Extensive Analysis with Special Emphasis on the B97M
-
V and ωB97M
-
V Approaches.
J. Chem. Theory Comput.
2018
,
14
, 5725
5738.
19
See: Yepes, D.; Valenzuela, J.; Martínez
-
Araya, J. I.; Pérez, P.; Jaque, P. Effect of the Exchange
Correlation
Functional on the Synchronicity/Nonsynchronicity in Bond Formation in Diels
Alder Reactions: A Reaction Force
Constant Analysis.
Phys. Chem. Chem. Phys.
2019
,
21
, 7412
7428.
20
Kesharwani, M. K.; Brauer, B.; Martin, J. M. L. Frequency and Zero
-
Point Vibrational Energy Scale Factors
for Double
-
Hybrid Density Functionals (and Other Selected Methods): Can Anharmonic Force Fields Be
Avoided?
J. Phys. Chem. A
2015
,
119
,
1701
1714.