|
Model Preparation: Reactant
compounds were modeled by a combination of Classical molecular mechanics (MM)
and quantum mechanical (QM) methods available in Hyperchem Version 8.10
(Hypercube, Inc., Gainsville, Florida, USA). Following initial construction and
rough 3D rendering, molecules were geometry optimized to ~0.1 kcal/mol-Å convergence
via Hyperchem’s MM+ default-parameter force field, with bond dipoles activated
to handle electrostatic interactions. Partial atomic charges were subsequently
assigned by performing a groundstate (zero-charge system) self-consistent field
(SCF), unrestricted Hartree-Fock (UHF) QM calculation using an RM1
semiempirical method (Rocha et al., 2006).
Models were then re-optimized (~0.1 kcal/mol-Å) by an alternating series of
MM/QM calculations using AMBER99 and RM1 methods, respectively. In certain
cases, models were further refined and optimized by carrying out ab initio
UHF calculations using the 6-31G** basis set and a spin multiplicity of 1. Ab
initio atomic charges and structures established in those instances,
however, reverted to those yielded by the RM1 method, which was used to control
the dynamic reaction simulations described below. Reactants and TIP3P water
models (Jorgensen et al., 1983)
were introduced into periodic micro-solvation cells (MSC) as described below in
preparation for reaction modeling. Dynamic Reaction Simulations:
Molecular (reactive) collisions were managed according to a novel “Stepped
Forced Molecular Dynamics” (SFMD) algorithm (see Ridgway
et al., 2017) that was scripted using Tcl/Tk (https://www.tcl.tk/software/tcltk/).
Reactants were placed in a cubic micro-solvation cell (MSC) containing
~75-100 TIP3P water molecules adjusted
to a density of ~1.0g/cc. Reactants were separated from one another by
approximately 10-15Å and independently randomly rotated along their XYZ axes to
scramble their initial pre-collision orientations. Following system (water +
reactants) optimization to ~0.1 kcal/mol-Å using the AMBER99 force field with a
constant non-scaled dielectric (Ɛ = 80), the oxidant molecule was accelerated toward the
center-of-mass of the target organic compound. To help avoid large repulsive
energy gradients over the reaction coordinate and reduce the likelihood of
collapsing the SCF calculations, the trajectory was broken into 4 to 6 discrete
“quantum runs” (QR) controlled by the RM1 method using a spin multiplicity of 3
and Pulay mixing, i.e., direct inversion of the iterative subspace; Pulay,
1980; 1982 to improve energy convergence. Each QR corresponded
to one segment of the total collision trajectory. At the end of each QR, the
reacting species were “frozen” and the water phase was re-optimized by a
suitable Newtonian force field (typically AMBER99) to permit adaptation of
water to the updated nuclear positions of the reactants. Following water
re-optimization, the aqueous phase was “frozen” and the reactants were allowed
to continue along their collision trajectory and adapt to the immobilized water
background, which effectively served as a perturbation on the RM1 electronic (i.e.,
clamped nucleus) Hamiltonian (Helec). The collision velocity (V1…Vn)
and duration (D1…Dn) of each QR were adjustable parameters and typically
gradually decremented (by about 10-15% per QR until a specified minimum “floor”
was reached) as the reactants approached one another, ostensibly to allow time
for reactant reorientation and reduction of repulsive non-bonded (Lennard-Jones
+ Coulomb) interactions. Standard initial velocity and QR duration were usually
set at ~75-100 Å/ps and 0.1 ps, respectively (using a 0.5 fs MD time step),
resulting in a total reaction time of approximately 0.2-0.5 ps; but this range
could vary depending on reactant chemistries (see below). Reactant “impact”, i.e.,
the point where molecular orbitals and covalent bond distances and angles began
to be noticeably affected, generally commenced at the end of QR1, or early in
the QR2 cycle. This QM/MM cycling process was continued until the last QR was
completed, at which time a custom subroutine (also coded in Tcl/Tk) was
implemented to “re-bond” the molecular system based on the degree of departure
of bond lengths from thermodynamically ideal values. Bonding and non-bonding
distance criteria were adjustable and based on prior RM1 studies (data not
presented). Initial QR velocities and total reaction duration depended to some
extent on the chemical nature of the reacting species; thus, numerous SFMD test
runs were conducted for each reactant pair to establish the initial parameters.
Whereas the SFMD algorithm incorporated provisions to allow water to enter into
the reaction coordinate (also via RM1), this feature was generally not used in
this study. Based on their mass, charge, elemental composition and atomic
coordination, the re-bonding subroutine automatically enumerated and
categorized all reaction products regardless of their chemical stability, e.g.,
hydrogen-abstracted, halogen substituted species, and so forth. The combined (i.e.,
the SFMD + rebonding) algorithms were iterated over multiple independent
collision events (typically ~1000) for each reactant pair to generate a
statistical distribution of stable (e.g.,
a formaldehyde molecule) and unstable (e.g., a carbocation) reaction
products, which were then converted to SMILES notation for further analysis and
processing. Product/fragment categories:
SP = stable products; C+ =
carbocations; H- = hydrogen
abstracted products; A/S = addition/substitution products.
References…
- William L.
Jorgensen, Jayaraman Chandrasekhar, and Jeffry D. Madura, Roger W. Impey and
Michael L. Klein (1983). “Comparison of simple potential functions for
simulating liquid water”. The Journal of Chemical Physics 79, 926; doi:
http://dx.doi.org/10.1063/1.445869
- Pulay, Péter (1980). "Convergence acceleration of
iterative sequences: The case of SCF iteration". Chemical Physics Letters.
73 (2): 393–398. doi:10.1016/0009-2614(80)80396-4.
- Pulay,
Péter (1982). "Improved SCF Convergence Acceleration". Journal of Computational
Chemistry. 3 (4): 556–560.
- H.F.
Ridgway, B. Mohan, X. Cui, K.J. Chua, and M.R. Islam (2017). “Molecular
dynamics simulation of gas-phase ozone reactions with sabinene and benzene”.
Journal of Molecular Graphics and Modelling 74 (2017) 241–250;
doi.org/10.1016/j.jmgm.2017.04.020
- Gerd Bruno Rocha, Ricardo Oliveira Freire, Alfredo
Mayall Simas*, and James J. P. Stewart (2006). “RM1: A reparameterization of
AM1 for H, C, N, O, P, S, F, Cl, Br, and I”. Journal of Computational Chemistry
27(10), 1101-1111; doi:10.1002/jcc.20425
|
