Combining Explicit Quantum Solvent with a Polarizable Continuum Model

Toward a numerically efficient description of bulk-solvated anionic states

We investigate the vertical electron attachment energy (VAE) of 1-methyl-4-nitroimidazole, a model radiosensitizer, employing quantum mechanics/molecular mechanics (QM/MM) and QM/polarized continuum (QM/PCM) solvation models. We considered the solvent-excluded surface (QM/PCM-SES) and Van der Waals (QM/PCM-VDW) cavities within the PCM framework, the electrostatic embedding QM/MM (EE-QM/MM) model, and the self-consistent sequential QM/MM polarizable electrostatic embedding (scPEE-S-QM/MM) model. Due to slow VAE convergence concerning the number of QM solvent molecules, full QM calculations prove inefficient. Ensemble averages in these calculations do not align with VAEs computed for the representative solute-solvent configuration. QM/MM and QM/PCM calculations show agreement with each other for sufficiently large QM regions, although the QM/PCM-VDW model exhibits artifacts linked to the cavity. QM/MM models demonstrate good agreement between ensemble averages and VAEs calculated with the representative configuration. Notably, the VAE computed with the scPEE-S-QM/MM model achieves faster convergence concerning the number of QM water molecules compared to the EE-QM/MM model, attributed to enhanced efficiency from MM charge polarization in the scPEE-S-QM/MM approach. This emphasizes the importance of QM/classical models with accurate solute-solvent and solvent-solvent mutual polarization for obtaining converged VAEs at a reasonable computational cost. The full-QM approach is very inefficient, while the microsolvation model is inaccurate. Computational savings in QM/MM models result from electrostatic embedding and the representative configuration, with the scPEE-S-QM/MM approach emerging as an efficient tool for describing bulk-solvated anions within the QM/MM framework. Its potential extends to improving transient anion state descriptions in biomolecules and radiosensitizers, especially given the frequent employment of microsolvation models.

Role of Explicit Solvation in Computational Modeling of Chemical Reactions: Mechanism of Cu(I) Transfer Between Thiolate-Based Chelators in Water

Solvation plays important roles in controlling the thermodynamic and kinetic aspects of chemical reactions. The conventional approaches to treat solvation via electronic structure methods are likely to become inadequate, when the reacting solutes have strong electrostatic and hydrogen bonding interactions with the solvent and undergo significant structural changes during the course of the reaction. In this article, we present evidence of such solvent and structural effects in the computational study of the Cu(I) transfer reaction between thiolate-based chelators dithiobutylamine (DTBA) and dithiotheritol (DTT) in water, inspired from biological copper trafficking phenomena. We propose a general solution to the problem by combining classical molecular dynamics (MD) simulations of the bulk system and static quantum chemistry calculations. The fluctuating solvation shell was estimated from MD, and energetics was assessed by averaging QM energies of a series of molecular clusters constructed from the MD snapshots. Applying this approach, we propose a reaction pathway with estimates of relative intermediate stabilities and barriers, which suggest the overall reaction to be reversible in nature and likely to go through both two and three coordinated intermediates, confirming previous studies of similar protein analogues. An interesting fact that emerged from our study was the strong indication that the rate-determining step is the deprotonation of initial thiol bound Cu(I) complex, without involving any Cu(I)-S bonds. The proposed method will lead to a better treatment of solvations, and these mechanistic insights will aid our understanding of biological copper(I) trafficking.

Prediction of Binding Pose and Affinity of Nelfinavir, a SARS-CoV-2 Main Protease Repositioned Drug, by Combining Docking, Molecular Dynamics, and Fragment Molecular Orbital Calculations

A novel in silico drug design procedure is described targeting the Main protease (Mpro) of the SARS-CoV-2 virus. The procedure combines molecular docking, molecular dynamics (MD), and fragment molecular orbital (FMO) calculations. The binding structure and properties of Mpro were predicted for Nelfinavir (NFV), which had been identified as a candidate compound through drug repositioning, targeting Mpro. Several poses of the Mpro and NFV complexes were generated by docking, from which four docking poses were selected by scoring with FMO energy. Then, each pose was subjected to MD simulation, 100 snapshot structures were sampled from each of the generated MD trajectories, and the structures were evaluated by FMO calculations to rank the pose based on binding energy. Several residues were found to be important in ligand recognition, including Glu47, Asp48, Glu166, Asp187, and Gln189, all of which interacted strongly with NFV. Asn142 is presumably regarded to form hydrogen bonds or CH/π interaction with NFV; however, in the present calculation, their interactions were transient. Moreover, the tert-butyl group of NFV had no interaction with Mpro. Identifying such strong and weak interactions provides candidates for maintaining and substituting ligand functional groups and important suggestions for drug discovery using drug repositioning. Besides the interaction between NFV and the amino acid residues of Mpro, the desolvation effect of the binding pocket also affected the ranking order. A similar procedure of drug design was applied to Lopinavir, and the calculated interaction energy and experimental inhibitory activity value trends were consistent. Our approach provides a new guideline for structure-based drug design starting from a candidate compound whose complex crystal structure has not been obtained.

How Aqueous Solvation Impacts the Frequencies and Intensities of Infrared Absorption Bands in Flavin: The Quest for a Suitable Solvent Model

FTIR spectroscopy accompanied by quantum chemical simulations can reveal important information about molecular structure and intermolecular interactions in the condensed phase. Simulations typically account for the solvent either through cluster quantum mechanical (QM) models, polarizable continuum models (PCM), or hybrid quantum mechanical/molecular mechanical (QM/MM) models. Recently, we studied the effect of aqueous solvent interactions on the vibrational frequencies of lumiflavin, a minimal flavin model, using cluster QM and PCM models. Those models successfully reproduced the relative frequencies of four prominent stretching modes of flavin's isoalloxazine ring in the diagnostic 1450-1750 cm-1 range but poorly reproduced the relative band intensities. Here, we extend our studies on this system and account for solvation through a series of increasingly sophisticated models. Only by combining elements of QM clusters, QM/MM, and PCM approaches do we obtain an improved agreement with the experiment. The study sheds light more generally on factors that can impact the computed frequencies and intensities of IR bands in solution.

Ab initio molecular dynamics free energy study of enhanced copper (II) dimerization on mineral surfaces

Understanding the adsorption of isolated metal cations from water on to mineral surfaces is critical for toxic waste retention and cleanup in the environment. Heterogeneous nucleation of metal oxyhydroxides and other minerals on material surfaces is key to crystal growth and dissolution. The link connecting these two areas, namely cation dimerization and polymerization, is far less understood. In this work we apply ab initio molecular dynamics calculations to examine the coordination structure of hydroxide-bridged Cu(II) dimers, and the free energy changes associated with Cu(II) dimerization on silica surfaces. The dimer dissociation pathway involves sequential breaking of two Cu²⁺-OH^- bonds, yielding three local minima in the free energy profiles associated with 0-2 OH^- bridges between the metal cations, and requires the design of a (to our knowledge) novel reaction coordinate for the simulations. Cu(II) adsorbed on silica surfaces are found to exhibit stronger tendency towards dimerization than when residing in water. Cluster-plus-implicit-solvent methods yield incorrect trends if OH^- hydration is not correctly depicted. The predicted free energy landscapes are consistent with fast equilibrium times (seconds) among adsorbed structures, and favor Cu²⁺ dimer formation on silica surfaces over monomer adsorption.

Progress on Exploring the Luminescent Properties of Organic Molecular Aggregates by Multiscale Modeling

Luminescent molecular aggregates have attracted worldwide attention because of their potential applications in many fields. The luminescent properties of organic aggregates are complicated and highly morphology-dependent, unraveling the intrinsic mechanism behind is urgent. This review summarizes recent works on investigating the structure-property relationships of organic molecular aggregates at different environments, including crystal, cocrystal, amorphous aggregate, and doped systems by multiscale modeling protocol. We aim to explore the influence of intermolecular non-covalent interactions on molecular packing and their photophysical properties and then pave the effective way to design, synthesize, and develop advanced organic luminescent materials.

Single-ion solvation free energy: A new cluster-continuum approach based on the cluster expansion method

Accurate calculation of the solvation free energy of single ions remains an important goal, involving development in the dielectric continuum solvation models, and statistical mechanics with explicit solvent and hybrid discrete-continuum methods. In the last case, many of the research studies involve a quasi-chemical approach using the monomer cycle or the cluster cycle to calculate the solvation free energy of single ions. In this work, a different cluster-continuum approach based on the cluster expansion method was tested for solvation of 16 cations and 32 anions in aqueous solution. The SMD model was used for the dielectric continuum part and three explicit water molecules were introduced in the region of the solute with the highest interaction energy. Harmonic frequency calculations and molecular dynamics sampling of configurations are not required. An empirical γ_N parameter for cations and another for anions is introduced. The method produces a substantial improvement of the SMD model with a mean absolute deviation of 2.3 kcal mol^-1 for cations and 2.9 kcal mol^-1 for anions. The analysis of the correlation between theoretical and experimental data produces a linear regression line with a slope of 1.09 for cations and 1.01 for anions. The good results of this approximated cluster expansion approach suggest that the method could be further improved by including more solvent molecules and sampling the configurations.

Hydration and Charge-Transfer Effects of Alkaline Earth Metal Ions Binding to a Carboxylate Anion, Phosphate Anion, and Guanine Nucleobase

To investigate the ability of alkaline earth metal ions to tune ion-mediated DNA adsorption, hydrated Mg²⁺, Ca²⁺, Sr²⁺, and Ba²⁺ ions bound to a carboxylate anion, phosphate anion, and guanine nucleobase were modeled using density functional theory (DFT) and a combined explicit and continuum solvent model. The large first solvation shell of Ba²⁺ requires a larger solute cavity defined by a solvent-accessible surface, which is used to model all hydrated ions. Alkaline earth metal ions bind indirectly or directly to each binding site. DFT binding energies decrease with increasing ion size, which is likely due to ion size and hydration structure, rather than quantum effects such as charge transfer. However, charge transfer explains weaker ion binding to guanine compared to phosphate or carboxylate. Overall, carboxylate and phosphate anions are expected to compete equally for hydrated Mg²⁺, Ca²⁺, Sr²⁺, and Ba²⁺ ions and larger alkaline earth metal ions may induce weaker ion-mediated adsorption. The ion size and hydration structure of alkaline earth metal ions may effectively tune ion-mediated adsorption processes, such as DNA adsorption to functionalized surfaces.

Accurate pKa Evaluations for Complex Bio-Organic Molecules in Aqueous Media

Despite the numerous computational efforts on estimating acid dissociation constant (pK_a's), an accurate estimation of pK_a's of bio-organic molecules in the aqueous medium is still a challenge. The major difficulty lies in the accurate description of the aqueous environment and the cost and accuracy of quantum mechanical (QM) methods. Herein, we report a well-defined quantum chemical protocol for accurately calculating pK_a's of a wide range of bio-organic molecules in aqueous media. The performance of our method has been assessed using test sets containing molecules with a range of sizes and a variety of functional groups, including alcohols, phenols, amines, and carboxylic acids, and obtained an impressive mean absolute accuracy of 0.5 pK_a units. For the smaller molecules, we use a high-level QM method (CBS-QB3) and a calibrated explicit-implicit solvation model that yields accurate pK_a values for a range of functional groups. For the larger molecules, we combine this approach with an efficient error-cancellation scheme that eliminates the systematic errors in different density functional methods to yield accurate pK_a values for simple to complex molecular systems. Our protocol is efficient, applicable to large molecules, covers all the common functional groups present in bio-organic molecules, and should find widespread applications in diverse research areas including drug-protein binding, catalysis, and chemical synthesis.

Tuning the photoreactivity of Z-hexatriene photoswitches by substituents - a non-adiabatic molecular dynamics study

To understand how substituents can be used to increase the quantum yield of photochemical electrocyclic ring-closing of the Z-hexa-1,3,5-triene (HT) photoswitch forming cyclohexadiene (CHD), we investigate the S1 photo dynamics of HT and its derivatives 2,5-dimethyl-HT (DMHT), 2-isopropyl-5-methyl-HT (2,5-IMHT), 1-isopropyl-4-methyl-HT (1,4-IMHT), and 2,5-diisopropyl-HT (DIHT) using time-dependent density functional theory surface hopping dynamics. We report detailed photoproduct distributions, formation mechanisms, branching ratios, and wavelength-dependent product quantum yields. Most products have been confirmed experimentally and include all-trans HT derivatives, cyclopropanes, cyclobutenes, cyclopentene, cyclohexadienes, and bicyclic compounds. Regarding CHD formation, we find that for the 2,5-substituted derivatives DMHT, 2,5-IMHT, and DIHT, the branching ratios increase with increasing size of the substituents. In contrast the branching ratios of the E/Z-isomerization decrease with increasing size of the substituents. Due to steric interactions, increasing the size of the substituents increases the amount of gZg rotamers in the ground state, which are prone to CHD formation and have lower E/Z-isomerization probability. Furthermore, we find [1,4], [1,5], and [1,6]-sigmatropic hydrogen shift reactions occurring at large percentages (5% to 15%); for sterical reasons these reactions stem from tZg conformers. DIHT shows the lowest percentage of side product formation among the 2,5-substituted molecules and highest CHD branching ratio; its CHD quantum yield can be increased up to more than 64%, by excitation of DIHT on the red tail of its absorption spectrum, whereas the Z/E-isomerization is reduced below 5% and side reactions practically vanish. This makes DIHT the best candidate for applications in molecular switches.

First-Principles Prediction of Wavelength-Dependent Product Quantum Yields

We present a method to predict wavelength-dependent product quantum yields (PQYs) for photochemical reactions and applied it to Z/E-isomerization and several ring-closing reactions of Z-2,5-dimethyl-1,3,5-hexatriene and truncated previtamin D. Using branching ratios from surface hopping molecular dynamics, individual trajectories are correlated with the absorption spectra of their initial structures. The wavelength-dependent PQYs are computed by dividing the average spectrum of the initial structures of the product-forming trajectories by the average spectrum of all initial structures. Accurate absorption spectra are calculated using the correlated ADC(2) method with an implicit solvent. Calculations reproduce the experimentally found trend of increasing six-ring formation and decreasing Z/E-isomerization on the red side of the spectrum. Over all seven reactions studied, the mean absolute error (MAE) between experimental and calculated PQYs (MAE) amounts to 8.1%, and the largest MAE is 18.6%. For four reactions, predicted values agree quantitatively with experiments within 5.6%.