Journal of Chemical Theory and Computation最新文献

In the pursuit of informing experimental techniques with in silico optimizations, we propose a pip deployable Python package, pothos, to easily determine polymer crystallites within molecular dynamic melts and the chain orientation parameters of atomistic and coarse-grained simulations. The package supports the commonly used ⟨P₂⟩, ⟨P₄⟩, and ⟨P₆⟩ order parameters based on the chain chord vector and utilizes a modified DBSCAN algorithm to determine crystalline regions. The results of analysis are written to text and LAMMPS dump files for visualization and analysis.

Nonadiabatic dynamics simulations complement time-resolved experiments by revealing ultrafast excited-state mechanistic information in photochemical reactions. Understanding the relaxation mechanisms of photoexcited molecules finds application in energy, material, and medicinal research. However, with substantial computational costs, the nonadiabatic dynamics simulations have been restricted to ultrafast timescales, typically less than a few picoseconds, thus neglecting a wide range of photoactivated processes occurring in much longer timescales. Before developing new methodologies, we must ask: How well do the popular nonadiabatic dynamics methods perform in a long timescale simulation? In this study, we employ the multiconfiguration time-dependent Hartree (MCTDH) and its multilayer variants (ML-MCTDH), ab initio multiple spawning (AIMS), and fewest-switches surface hopping (FSSH) methodologies to simulate the excited-states dynamics of a weakly coupled multidimensional Spin-Boson model Hamiltonian designed for a long timescale decay behavior. Our study assures that despite having very different theoretical backgrounds, all the above methods deliver qualitatively similar results. While quantum dynamics would be very costly for long timescale simulations, the trajectory-based approaches are paving the way for future advancements.

High-harmonic generation (HHG) is a nonlinear process in which a material sample is irradiated by intense laser pulses, causing the emission of high harmonics of incident light. HHG has historically been explained by theories employing a classical electromagnetic field, successfully capturing its spectral and temporal characteristics. However, recent research indicates that quantum-optical effects naturally exist or can be artificially induced in HHG, such as entanglement between emitted harmonics. Even though the fundamental equations of motion for quantum electrodynamics (QED) are well-known, a unifying framework for solving them to explore HHG is missing. So far, numerical solutions have employed a wide range of basis-sets, methods, and untested approximations. Based on methods originally developed for cavity polaritonics, here we formulate a numerically accurate QED model consisting of a single active electron and a single quantized photon mode. Our framework can, in principle, be extended to higher electronic dimensions and multiple photon modes to be employed in ab initio codes for realistic physical systems. We employ it as a model of an atom interacting with a photon mode and predict a characteristic minimum structure in the HHG yield vs phase-squeezing. We find that this phenomenon, which can be used for novel ultrafast quantum spectroscopies, is partially captured by a multitrajectory Ehrenfest dynamics approach, with the exact minima position sensitive to the level of theory. On the one hand, this motivates using multitrajectory approaches as an alternative for costly exact calculations. On the other hand, it suggests an inherent limitation of the multitrajectory formalism, indicating the presence of entanglement and true quantum effects (especially prominent for atomic and molecular resonances). Our work creates a roadmap for a universal formalism of QED-HHG that can be employed for benchmarking approximate theories, predicting novel phenomena for advancing quantum applications, and for the measurements of entanglement and entropy.

In this work, we develop a novel Bayesian approach to study the adsorption and desorption of CO onto a Pd(111) surface, a process of great importance in natural sciences. The motivation for this work comes from the recent availability of time-resolved infrared spectroscopy data and the need for model interpretability and uncertainty quantification in chemical processes. The objective is to learn the relevant parameters that characterize the process: coverage with time, rate constants, activation energies, and pre-exponential factors. Our approach consists of three main schemes: (i) a problem design and probabilistic model for the whole system, (ii) a particle Markov chain Monte Carlo sampler to learn the hidden coverages and rate constant parameters, and (iii) two Bayesian formulations to infer the activation energies and pre-exponential factors. The flexibility of the Bayesian framework allows for uncertainty quantification where possible and integration of mathematical constraints in the model to reflect the system physically. We found that our results for the activation energies and pre-exponential factor are in agreement with those reported in the experimental literature, independently, and we provide discussions on the advantages and disadvantages as well as applicability to other systems.

Self-diffusion coefficients, D*, are routinely estimated from molecular dynamics simulations by fitting a linear model to the observed mean squared displacements (MSDs) of mobile species. MSDs derived from simulations exhibit statistical noise that causes uncertainty in the resulting estimate of D*. An optimal scheme for estimating D* minimizes this uncertainty, i.e., it will have high statistical efficiency, and also gives an accurate estimate of the uncertainty itself. We present a scheme for estimating D* from a single simulation trajectory with a high statistical efficiency and accurately estimating the uncertainty in the predicted value. The statistical distribution of MSDs observable from a given simulation is modeled as a multivariate normal distribution using an analytical covariance matrix for an equivalent system of freely diffusing particles, which we parametrize from the available simulation data. We use Bayesian regression to sample the distribution of linear models that are compatible with this multivariate normal distribution to obtain a statistically efficient estimate of D* and an accurate estimate of the associated statistical uncertainty.

In this work, we describe a computational tool designed to determine the local dielectric constants (ε) of charge-neutral heterogeneous systems by analyzing dipole moment fluctuations from molecular dynamics (MD) trajectories. Unlike conventional methods, our tool can calculate dielectric constants for dynamically evolving selections of molecules within a defined region of space, rather than for fixed sets of molecules. We validated our approach by computing the dielectric constants of TIP3P water nanospheres, achieving results consistent with literature values for bulk water. We then applied our tool to more complex systems, the water slabs around solvated phospholipid bilayers, where we observed a lower dielectric constant of water near the bilayer headgroups (ε = 20-50) compared to nanospheres of bulk water (ε = 58-62) with the same number of molecules. Our tool also enabled us to compute the dielectric constants of water in more heterogeneous systems, where water surrounding asymmetrically distributed phospholipids on single-walled carbon nanotubes also exhibited lower dielectric constants than in bulk water nanospheres. Addition of positively charged peptides that bind to phospholipid-nanotube conjugates further lowered the dielectric constants of water in the immediate vicinity of these conjugates. Moreover, we estimated dielectric constants for lipids in symmetric bilayers, where values are well-documented, and for asymmetric phospholipid-wrapped nanotube systems, which were previously unexplored, and found that dielectric constants of phospholipids depend on their arrangement in the assembled aggregate. The results align with the literature for bilayers and provide new insights for phospholipid-nanotube systems. The ability of our tool to provide local dielectric constants for both well-studied and novel systems advances our understanding of molecular environments and interactions.

The successful simulation of proteins by molecular dynamics (MD) critically depends on the accuracy of the applied force field. Here, we modify the AMBER-family ff99SBnmr2 force field through improvements to the side-chain χ₁ dihedral angle potentials in a residue-specific manner using conformational dihedral angle distributions from an experimental coil library as targets. Based on significant deviations observed for the parent force field with respect to the coil library, the χ₁ dihedral angle potentials of seven amino acids were modified, namely, Val, Ser, His, Asn, Trp, Tyr, and Phe. The new force field, named ff99SBnmr2Chi1, was benchmarked against NMR-derived χ₁ rotamer populations of denatured proteins, overall resulting in much better agreement and without any noticeable adverse consequences on the quality of the simulation of folded proteins. The new force field should allow more realistic modeling of protein side-chain properties by MD of both folded and unfolded protein systems, such as for the better in-silico characterization of protein-protein and protein-ligand interactions.

Free energy calculations for protein-ligand complexes have become widespread in recent years owing to several conceptual, methodological and technological advances. Central among these is the use of ensemble methods which permits accurate, precise and reproducible predictions and is necessary for uncertainty quantification. Absolute binding free energies (ABFEs) are challenging to predict using alchemical methods and their routine application in drug discovery has remained out of reach until now. Here, we apply ensemble alchemical ABFE methods to a large data set comprising 219 ligand-protein complexes and obtain statistically robust results with high accuracy (<1 kcal/mol). We compare equilibrium and nonequilibrium methods for ABFE predictions at large scale and provide a systematic critical assessment of each method. The equilibrium method is more accurate, precise, faster, computationally more cost-effective and requires a much simpler protocol, making it preferable for large scale and blind applications. We find that the calculated free energy distributions are non-normal and discuss the consequences. We recommend a definitive protocol to perform ABFE calculations optimally. Using this protocol, it is possible to perform thousands of ABFE calculations within a few hours on modern exascale machines.

The evolution of photosynthetic reaction centers (RCs) from anoxygenic bacteria to higher-order oxygenic cynobacteria and plants highlights a remarkable journey of structural and functional diversification as an adaptation to environmental conditions. The role of chirality in these centers is important, influencing the arrangement and function of key molecules involved in photosynthesis. Investigating the role of chirality may provide a deeper understanding of photosynthesis and the evolutionary history of life on Earth. In this study, we explore chirality-related energy transfer in two evolutionarily distinct RCs: one from the anoxygenic purple sulfur bacterium Thermochromatium tepidum (BRC) and the other from the oxygenic cyanobacterium Thermosynechococcus vulcanus (PSII RC), utilizing two-dimensional electronic spectroscopy (2DES). By employing circularly polarized laser pulses, we can extract transient chiral dynamics within these RCs, offering a detailed view of their chiral contribution to energy transfer processes. We also compute traditional 2DES and compare these results with spectra related to circular dichroism. Our findings indicate that two-dimensional circular dichroism spectroscopy effectively reveals chiral dynamics, emphasizing the structural symmetries of pigments and their interactions with associated proteins. Despite having similar pigment-protein architectures, the BRC and PSII RC exhibit significantly different chiral dynamics on an ultrafast time scale. In the BRC, the complex contributions of pigments such as BCh_M, BPh_L, BCh, and P_M to key excitonic states lead to more pronounced chiral features and dynamic behavior. In contrast, the PSII RC, although significantly influenced by Chl_D1 and Chl_D2, shows less complex chiral effects and more subdued chiral dynamics. Notably, the PSII RC demonstrates a faster decay of coherence to localized excitonic populations compared to the BRC, which may represent an adaptive mechanism to minimize oxidative stress in oxygenic photosystems. By examining and comparing the chiral excitonic interactions and dynamics of BRC and PSII RC, this study offers valuable insights into the mechanisms of photosynthetic complexes. These findings could contribute to understanding how the functional optimization of photosynthetic proteins in ultrafast time scales is linked to biological evolution.

Classical density functional theory (DFT) provides a versatile framework to study the polymers with complex topological structure. Generally, a classical DFT describes the excess Helmholtz free energy of nonbonded chain connectivity due to excluded-volume effects and electrostatic correlations using the first-order thermodynamic perturbation theory (referred to as DFT-TPT1). Beyond first-order perturbation, the second-order TPT (TPT2) captures not only the correlations between neighboring monomers but also the interactions within three consecutive monomers, playing a crucial role in describing the polymer topology. However, the numerical implementation of TPT2 is limited by the lack of an effective triple correlation function (CF), especially for charged systems. Here, we propose an effective triple CF and incorporate it into DFT using TPT2 (referred to as DFT-eTPT2) to describe the nonbonded chain connectivity due to excluded-volume effects and electrostatic correlations. Using the data from molecular dynamics simulation as a benchmark, DFT-eTPT2 shows a clear improvement over DFT-TPT1 in predicting the density profiles of both neutral and charged branched polymer brushes, accurately capturing key structural features, such as the significant peaks near the branching point in the density profiles. In short, this work provides a precise and efficient theoretical tool for revealing molecular-level insights into branched polymers and their brushes.