Functional group corrections to the GFN2-xTB and PM6 semiempirical methods for noncovalent interactions in alkanes and alkenes.

Analytical corrections were developed to improve the accuracy of the PM6 and GFN2-xTB semiempirical quantum mechanical methods for the evaluation of noncovalent interaction energies in alkanes and alkenes. We followed the approach of functional group corrections, wherein the atom-atom pair corrections depend on the nature of the interacting functional groups. The training set includes 21 alkane and 13 alkene complexes taken from the Donchev et al.'s database [Sci. Data 8, 55 (2021)], with interaction energies calculated at the CCSD(T)/CBS level, and our own data obtained for medium-size complexes (of 100 and 112 atoms). In general, for the systems included in the training and validation sets, the errors obtained with the PM6-FGC and xTB-FGC methods are within the chemical accuracy.

Electrostatic penetration effects stand at the heart of aromatic π interactions.

The nature of the interaction in benzene-containing dimers has been analysed by means of Symmetry Adapted Perturbation Theory (SAPT). The total interaction energy and the preference for the dimers to adopt slipped structures are, apparently, consequence of the balance between repulsion and dispersion. However, our results indicate that this only holds when trends are analysed using fixed intermolecular distances. Employing the most favourable separations between rings it turns out that the changes on the total interaction energy are mostly controlled by electrostatics, while repulsion and dispersion cancel each other to a great extent. Most of the electrostatic contribution is accounted for by electrostatic penetration, so a description based on multipoles should not be employed to rationalise the interaction in benzene-containing dimers. The changes on the interaction energy in benzene-containing dimers are steered by electrostatic penetration which, though often overlooked, plays an essential role for the description of aromatic π interactions.

The PM6-FGC Method: Improved Corrections for Amines and Amides.

Recently, we reported a new approach to develop pairwise analytical corrections to improve the description of noncovalent interactions, by approximate methods of electronic structures, such as semiempirical quantum mechanical (SQM) methods. In particular, and as a proof of concept, we used the PM6 Hamiltonian and we named the method PM6-FGC, where the FGC acronym, corresponding to Functional Group Corrections, emphasizes the idea that the corrections work for specific functional groups rather than for individual atom pairs. The analytical corrections were derived from fits to B3LYP-D3/def2-TZVP (reference). PM6 interaction energy differences, evaluated for a reduced set of small bimolecular complexes, were chosen as representatives of saturated hydrocarbons, carboxylic, amine and, tentatively, amide functional groups. For the validation, the method was applied to several complexes of well-known databases, as well as to complexes of diglycine and dialanine, assuming the transferability of amine group corrections to amide groups. The PM6-FGC method showed great potential but revealed significant inaccuracies for the description of some interactions involving the -NH2 group in amines and amides, caused by the inadequate selection of the model compound used to represent these functional groups (an NH3 molecule). In this work, methylamine and acetamide are used as representatives of amine and amide groups, respectively. This new selection leads to significant improvements in the calculation of noncovalent interactions in the validation set.

AutoMeKin2021: An open-source program for automated reaction discovery.

AutoMeKin2021 is an updated version of tsscds2018, a program for the automated discovery of reaction mechanisms (J. Comput. Chem. 2018, 39, 1922). This release features a number of new capabilities: rare-event molecular dynamics simulations to enhance reaction discovery, extension of the original search algorithm to study van der Waals complexes, use of chemical knowledge, a new search algorithm based on bond-order time series analysis, statistics of the chemical reaction networks, a web application to submit jobs, and other features. The source code, manual, installation instructions and the website link are available at: https://rxnkin.usc.es/index.php/AutoMeKin.

The relative position of π-π interacting rings notably changes the nature of the substituent effect.

The substituent effect in monosubstituted benzene dimers mostly follows changes on electrostatics mainly controlled by the direct interaction of the substituent and the other phenyl ring, whereas the contribution from the interacting rings is smaller. As the substituent is located further away the two contributions become of similar magnitude, so the global result is a combination of both effects. These trends are confirmed in larger systems containing a contact between phenyl rings; at closer distances the interaction of the substituent and the other ring clearly dominates over changes associated with the substituted ring, but as the substituent is located further away its contribution decreases and the contribution from the ring becomes more relevant. Care should be taken in larger systems because the observed energy change can also be affected by interactions with other regions of the molecule not directly involved in the π-π interaction.

Influence of Multiple Conformations and Paths on Rate Constants and Product Branching Ratios. Thermal Decomposition of 1-Propanol Radicals.

The potential energy surface involved in the thermal decomposition of 1-propanol radicals was investigated in detail using automated codes (tsscds2018 and Q2DTor). From the predicted elementary reactions, a relevant reaction network was constructed to study the decomposition at temperatures in the range 1000-2000 K. Specifically, this relevant network comprises 18 conformational reaction channels (CRCs), which in general exhibit a large wealth of conformers of reactants and transition states. Rate constants for all the CRCs were calculated using two approaches within the formulation of variational transition-state theory (VTST), as incorporated in the TheRa program. The simplest, one-well (1W) approach considers only the most stable conformer of the reactant and that of the transition state. In the second, more accurate approach, contributions from all the reactant and transition-state conformers are taken into account using the multipath (MP) formulation of VTST. In addition, kinetic Monte Carlo (KMC) simulations were performed to compute product branching ratios. The results show significant differences between the values of the rate constants calculated with the two VTST approaches. In addition, the KMC simulations carried out with the two sets of rate constants indicate that, depending on the radical considered as reactant, the 1W and the MP approaches may display different qualitative pictures of the whole decomposition process.

A Trajectory-Based Method to Explore Reaction Mechanisms.

The tsscds method, recently developed in our group, discovers chemical reaction mechanisms with minimal human intervention. It employs accelerated molecular dynamics, spectral graph theory, statistical rate theory and stochastic simulations to uncover chemical reaction paths and to solve the kinetics at the experimental conditions. In the present review, its application to solve mechanistic/kinetics problems in different research areas will be presented. Examples will be given of reactions involved in photodissociation dynamics, mass spectrometry, combustion chemistry and organometallic catalysis. Some planned improvements will also be described.

Photodissociation of acryloyl chloride at 193 nm: interpretation of the product energy distributions, and new elimination pathways.

The ground electronic state potential energy surface of acryloyl chloride, CH2CHC(O)Cl, has been mapped using an automated transition state search procedure. A total of 174 minima, 527 TSs, and 20 different dissociation channels have been found. Among others, three novel HCl elimination pathways, namely, a five-center mechanism and two three-body dissociations (leading to CO + HCl + HCCH) have been discovered. While the bimodal character of the experimental HCl rotational distributions was previously attributed to the presence of two competing channels, our dynamics simulations show that a single channel, the four-center HCl elimination of CH2ClCHCO following a 1,3-Cl-shift of CH2CHC(O)Cl, displays a bimodal distribution in nearly prefect agreement with the experiment. Overall, our simulation results suggest that, as far as molecular elimination is concerned, this channel dominates in the 193 nm photodissociation of the molecule. The simulations also show evidence of non-IRC dynamics for this channel.

HCN elimination from vinyl cyanide: product energy partitioning, the role of hydrogen-deuterium exchange reactions and a new pathway.

The different HCN elimination pathways from vinyl cyanide (VCN) are studied in this paper using RRKM, Kinetic Monte Carlo (KMC), and quasi-classical trajectory (QCT) calculations. A new HCN elimination pathway proves to be very competitive with the traditional 3-center and 4-center mechanisms, particularly at low excitation energies. However, low excitation energies have never been experimentally explored, and the high and low excitation regions are dynamically different. The KMC simulations carried out using singly deuterated VCN (CH2=CD-CN) at 148 kcal mol(-1) show the importance of hydrogen-deuterium exchange reactions: both DCN and HCN will be produced in any of the 1,1 and 1,2 elimination pathways. The QCT simulation results obtained for the 3-center pathway are in agreement with the available experimental results, with the 4-center results showing much more excitation of the products. In general, our results seem to be consistent with a photodissociation mechanism at 193 nm, where the molecule dissociates (at least the HCN elimination pathways) in the ground electronic state. However, our simulations assume that internal conversion is a fully statistical process, i.e., the HCN elimination channels proceed on the ground electronic state according to RRKM theory, which might not be the case. In future studies it would be of interest to include the photo-prepared electronically excited state(s) in the dynamics simulations.

Direct and indirect hydrogen abstraction in Cl + alkene reactions.

Reactions between Cl atoms and propene can lead to HCl formation either by direct H abstraction or through a chloropropyl addition complex. Barring stabilizing collisions, the chloropropyl radical will either decompose to reactants or form HCl and allyl products. Using velocity-map imaging to measure the quantum state and velocity of the HCl products provides a view into the reaction dynamics, which show signs of both direct and indirect reaction mechanisms. Simulated trajectories of the reaction highlight the role of the direct H-abstraction pathways, and the resultant simulated scattering images show reasonable agreement with measurement. The simulations also show the importance of large excursions of the Cl atom far from equilibrium geometries within the chloropropyl complex, and these large-amplitude motions are the ultimate drivers toward HCl + allyl fragmentation. Gas-phase measurements of larger alkenes, 2-methylpropene and 2,3-dimethylbut-2-ene, show slightly different product distributions but still feature similar reaction dynamics. The current suite of experiments offers ready extensions to liquid-phase bimolecular reactions.

Collision-induced dissociation mechanisms of [Li(uracil)]+.

The collision-induced dissociation (CID) of the [Li(uracil)](+) complex with Xe is studied by means of quasi-classical trajectory calculations. The potential energy surface is obtained "on the fly" from AM1 semiempirical calculations, supplemented with two-body analytical potentials to model the intermolecular interactions. The simulations show that Li(+) production is the primary channel, in agreement with a previous experimental study [M. T. Rodgers and P. B. Armentrout, J. Am. Chem. Soc., 2000, 122, 8548]. Collision-induced isomerization of [Li(uracil)](+) was found to be very important as well in the 2.5-10 eV collision energy range. Three minor channels are also identified: complex formation between Xe and [Li(uracil)](+), ligand exchange to form LiXe(+), and fragmentations of the uracil ring, which are strongly nonstatistical. Additional quasi-classical trajectory calculations carried out to investigate in more detail the fragmentations of the uracil ring reveal the presence of bifurcations in the potential energy surface, as trajectories starting from a transition state give rise to four different product channels. The integral cross sections for Li(+) production calculated in this work agree well with those obtained in the experiments only for the lowest collision energies, being ∼20 times greater than the experimental values for a collision energy of 10 eV. Finally, the initial translational energy is transferred preferentially to the [Li(uracil)](+) vibrational degrees of freedom, with energy transfer to rotation being modest. The amount of energy transfer to the different degrees of freedom as a function of the collision energy follows quite nicely a model recently proposed by our group.

Ab initio and RRKM study of the HCN/HNC elimination channels from vinyl cyanide.

Ab initio CCSD and CCSD(T) calculations with the 6-311+G(2d,2p) and the 6-311++G(3df,3pd) basis sets were carried out to characterize the vinyl cyanide (C(3)H(3)N) dissociation channels leading to hydrogen cyanide (HCN) and its isomer hydrogen isocyanide (HNC). Our computations predict three elimination channels giving rise to HCN and another four channels leading to HNC formation. The relative HCN/HNC branching ratios as a function of internal energy of vinyl cyanide were computed using RRKM theory and the kinetic Monte Carlo method. At low internal energies (120 kcal/mol), the total HCN/HNC ratio is about 14, but at 148 kcal/mol (193 nm) this ratio becomes 1.9, in contrast with the value 124 obtained in a previous ab initio/RRKM study at 193 nm (Derecskei-Kovacs, A.; North, S. W. J. Chem. Phys.1999, 110, 2862). Moreover, our theoretical results predict a ratio of rovibrationally excited acetylene over total acetylene of 3.3, in perfect agreement with very recent experimental measurements (Wilhelm, M. J.; Nikow, M.; Letendre, L.; Dai, H.-L. J. Chem. Phys.2009, 130, 044307).

Vibrational Energy Relaxation of Deuterium Fluoride in d-Dichloromethane: Insights from Different Potentials.

Vibrationally excited deuterium fluoride (DF) formed by fluorine atom reaction with a solvent was found (Science, 2015, 347, 530) to relax rapidly (less than 10 ps) in acetonitrile-d3 (CD3CN) and dichloromethane-d2 (CD2Cl2). However, insights into how CD2Cl2 facilitates this energy relaxation have so far been lacking, given the weak interaction between DF and a single CD2Cl2. In this work, we report the results of reactive simulations with a two-state reactive empirical valence bond (EVB) potential to study the energy deposited into nascent DF after transition-state passage and of nonequilibrium molecular dynamics simulations using multiple different potential energy functions to model the relaxation dynamics. For these second simulations, we used the standard Merck molecular force field (MMFF) potential, an MMFF-based covalent-ionic empirical valence bond (EVB) potential (EVBCI), a newly developed potential [referred to as MMFF(rDF)] which extends upon the MMFF potential by making the DF/CD2Cl2 interaction depend on the value of the D-F bond stretching coordinate and by taking the anisotropic charge distribution of the solvent molecules into account, the polarizable atomic multipole optimized energetics for biomolecular applications (AMOEBA) potential, and the quantum mechanics/molecular mechanics (QM/MM) potential. The relaxation is revealed to be highly sensitive to the potential used. Neither standard MMFF nor EVBCI reproduces the experimentally observed rapid relaxation dynamics, and they also fail to provide a good description of the interaction potential between DF and CD2Cl2 as calculated using CCSD(T)-F12. This is attributed to the use of a point-charge model for the solute and to failing to model the anisotropic electrostatic properties of CD2Cl2. The MMFF(rDF), AMOEBA, and QM/MM potentials all reproduce the CCSD(T)-F12 two-body DF---CD2Cl2 interaction potential rather well but only with the QM/MM approach is fast vibrational relaxation obtained (lifetimes of ∼288, ∼186, and ∼8 ps, respectively), which we attribute to differences in the solute-solvent local structure. With QM/MM, a unique "many-body" interaction pattern in which DF is in close contact with two solvent Cl atoms and more than three solvent D atoms is found, but this structure is not seen with other potentials. The QM/MM dynamics also display enhanced solute-solvent interactions with vibrationally excited DF that induce a DF band redshift and hence a resonant overlap with solvent C-D modes, which facilitate the intermolecular energy transfer. Our work also suggests that potentials used to model energy relaxation need to capture the fine structure of solute-solvent interactions and not just the two-body part.

New Approach for Correcting Noncovalent Interactions in Semiempirical Quantum Mechanical Methods: The Importance of Multiple-Orientation Sampling.

A new approach is presented to improve the performance of semiempirical quantum mechanical (SQM) methods in the description of noncovalent interactions. To show the strategy, the PM6 Hamiltonian was selected, although, in general, the procedure can be applied to other semiempirical Hamiltonians and to different methodologies. A set of small molecules were selected as representative of various functional groups, and intermolecular potential energy curves (IPECs) were evaluated for the most relevant orientations of interacting molecular pairs. Then, analytical corrections to PM6 were derived from fits to B3LYP-D3/def2-TZVP reference-PM6 interaction energy differences. IPECs provided by the B3LYP-D3/def2-TZVP combination of the electronic structure method and basis set were chosen as the reference because they are in excellent agreement with CCSD(T)/aug-cc-pVTZ curves for the studied systems. The resulting method, called PM6-FGC (from functional group corrections), significantly improves the performance of PM6 and shows the importance of including a sufficient number of orientations of the interacting molecules in the reference data set in order to obtain well-balanced descriptions.

Dynamics of CO2 scattering off a perfluorinated self-assembled monolayer. Influence of the incident collision energy, mass effects, and use of different surface models.

The dynamics of collisions of CO2 with a perfluorinated alkanethiol self-assembled monolayer (F-SAM) on gold were investigated by classical trajectory calculations using explicit atom (EA) and united atom (UA) models to represent the F-SAM surface. The CO2 molecule was directed perpendicularly to the surface at initial collision energies of 1.6, 4.7, 7.7, and 10.6 kcal/mol. Rotational distributions of the scattered CO2 molecules are in agreement with experimental distributions determined for collisions of CO2 with liquid surfaces of perfluoropolyether. The agreement is especially good for the EA model. The role of the mass in the efficiency of the energy transfer was investigated in separate simulations in which the mass of the F atoms was replaced by either that of hydrogen or chlorine, while keeping the potential energy function unchanged. The calculations predict the observed trend that less energy is transferred to the surface as the mass of the alkyl chains increases. Significant discrepancies were found between results obtained with the EA and UA models. The UA surface leads to an enhancement of the energy transfer efficiency in comparison with the EA surface. The reason for this is in the softer structure of the UA surface, which facilitates transfer from translation to interchain vibrational modes.

Internal energy of HCl upon photolysis of 2-chloropropene at 193 nm investigated with time-resolved Fourier-transform spectroscopy and quasiclassical trajectories.

Following photodissociation of 2-chloropropene (H(2)CCClCH(3)) at 193 nm, vibration-rotationally resolved emission spectra of HCl (upsilon < or = 6) in the spectral region of 1900-2900 cm(-1) were recorded with a step-scan time-resolved Fourier-transform spectrometer. All vibrational levels show a small low-J component corresponding to approximately 400 K and a major high-J component corresponding to 7100-18,700 K with average rotational energy of 39+/-(3)(11) kJ mol(-1). The vibrational population of HCl is inverted at upsilon = 2, and the average vibrational energy is 86+/-5 kJ mol(-1). Two possible channels of molecular elimination producing HCl + propyne or HCl + allene cannot be distinguished positively based on the observed internal energy distribution of HCl. The observed rotational distributions fit qualitatively with the distributions of both channels obtained with quasiclassical trajectories (QCTs), but the QCT calculations predict negligible populations for states at small J. The observed vibrational distribution agrees satisfactorily with the total QCT distribution obtained as a weighted sum of contributions from both four-center elimination channels. Internal energy distributions of HCl from 2-chloropropene and vinyl chloride are compared.

Intermolecular potential for binding of protonated peptide ions with perfluorinated hydrocarbon surfaces.

An analytic potential energy function was developed to model both short-range and long-range interactions between protonated peptide ions and perfluorinated hydrocarbon chains. The potential function is defined as a sum of two-body potentials of the Buckingham form. The parameters of the two-body potentials were obtained by fits to intermolecular potential energy curves (IPECs) calculated for CF4, which represents the F and C atoms of a perfluoroalkane chain, interacting with small molecules chosen as representatives of the main functional groups and atoms present in protonated peptide ions: specifically, CH4, NH3, NH4(+), and HCOOH. The IPECs were calculated at the MP2/aug-cc-pVTZ level of theory, with basis set superposition error (BSSE) corrections. Good fits were obtained for an energy range extending up to about 400 kcal/mol. It is shown that the pair potentials derived from the NH3/CF4 and HCOOH/CF4 fits reproduce acceptably well the intermolecular interactions in HCONH2/CF4, which indicates that the parameters obtained for the amine and carbonyl atoms may be transferable to the corresponding atoms of the amide group. The derived potential energy function may be used in chemical dynamics simulations of collisions of peptide-H(+) ions with perfluorinated hydrocarbon surfaces.

Interaction and dimerization energies in methyl-blocked alpha,gamma-peptide nanotube segments.

The building blocks of a promising class of peptide nanotubes composed of alternating D-alpha-amino acids and (1R,3S)-3-aminocyclohexane (or cyclopentane) carboxylic acid (D-gamma-Ach or D-gamma-Acp) were explored by computational methods. Specifically, density functional theory (DFT) calculations on monomers and dimers of gamma-Ach-based and gamma-Acp-based alpha,gamma-cyclo-hexapeptides and cyclo-octapeptides were carried out to investigate the experimentally observed preference for alpha-alpha over gamma-gamma dimerization, associated with the two types of stacking patterns present in these peptide nanotubes, as well as the preference for heterodimerization versus homodimerization. Full geometry optimizations were performed at the B3LYP/6-31G(d) level, and single point calculations were subsequently carried out with the B3LYP and M05-2X functionals and the 6-31+G(d,p) basis set. The calculations predict that the interaction energies in the alpha-alpha species are quite similar to those in the gamma-gamma dimers. However, a comparison of dimerization energies (i.e., interaction energies plus deformation energies of monomers) shows that alpha-alpha dimerization is energetically favored over gamma-gamma dimerization. The calculations strongly suggest that the preference for alpha-alpha binding is governed by differences between the deformation energies in the alpha and gamma monomers, rather than by differences between the relative strengths of the alpha-alpha and gamma-gamma hydrogen-bonding patterns. Calculations based on local properties of the electron density support the previous suggestion that the H-N bonds of the alpha-amino acids are more polarized than those of the gamma-amino acids.

Hydrogen transfer vs proton transfer in 7-hydroxy-quinoline.(NH3)3: a CASSCF/CASPT2 study.

Multiconfigurational CASSCF and CASPT2 calculations were performed to investigate the enol --> keto tautomerization in the lowest singlet excited state of the 7-hydroxyquinoline.(NH3)3 cluster. Two different reaction mechanisms were explored. The first one corresponds to that proposed previously by Tanner et al. (Science 2003, 302, 1736) on the basis of experimental observations and CASSCF optimizations under Cs-symmetry constraints. This mechanism comprises four consecutive steps and involves nonadiabatic transitions between the valence 1pipi* state and a pisigma* Rydberg-type state, resulting in hydrogen-atom transfer. Single-point CASPT2 calculations corroborate that for Cs-symmetry pathways hydrogen-atom transfer is clearly preferred over proton transfer. The second mechanism, predicted by CASSCF optimizations without constraints, implies proton transfer along a pathway on the 1pipi* surface in which one or more ammonia molecules depart significantly from the molecular plane defined by the hydroxyquinoline ring. The results suggest that both mechanisms may be competitive with proton transfer being somewhat favorable over hydrogen-atom transfer.

Direct-dynamics VTST study of the [1,7] hydrogen shift in 7-methylocta-1,3(Z),5(Z)-triene. A model system for the hydrogen transfer reaction in previtamin D3.

Direct-dynamics canonical variational transition-state theory calculations with microcanonically optimized multidimensional transmission coefficient (CVT/muOMT) for tunneling were carried out at the MPWB1K/6-31+G(d,p) level to study the [1,7] sigmatropic hydrogen rearrangement in 7-methylocta-1,3(Z),5(Z)-triene. This compound has seven conformers, of which only one leads to products, although all of them have to be included in the theoretical treatment. The calculated CVT/muOMT rate constants are in good agreement with the available experimental data. To try to understand the role of tunneling in the hydrogen shift reaction, we have also calculated the thermal rate constants for the monodeuterated compound in the interval T = 333.2-388.2 K. This allowed us to evaluate primary kinetic isotope effects (KIEs) and make a direct comparison with the experiment. Our calculations show that both the large measured KIE and the large measured difference in the activation energies between the deuterated and root compounds are due to the quantum tunneling. The tunneling contribution to the KIE becomes noticeable only when the coupling between the reaction coordinate and the transverse modes is taken into account. Our results confirm previous experimental and theoretical works, which guessed that the obtained kinetic parameters pointed to a reaction with an important contribution due to tunneling. The above conclusion would be essentially valid for the case of the [1,7] hydrogen shift in previtamin D3 because of the similarity to the studied model system.