Dynamic Reaction Simulations...
("Stepped Forced Molecular Dynamics")

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.



  • 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
Return to Home Page...