Professor Tom Markland focuses on problems at the interface of quantum mechanics and statistical mechanics, with applications ranging from chemistry and biology to geology and materials science. Markland Group research frequently explores theories of hydrogen bonding, the interplay between structure and dynamics, systems with multiple time and length-scales, and quantum mechanical effects. Particular current interests include proton and electron transfer in materials and enzymatic systems, atmospheric isotope separation, and the control of catalytic chemical reactivity in heterogeneous environments.

Thomas E. Markland studied chemistry at Balliol College, University of Oxford (MChem 2006), where as a Brackenbury Scholar he performed thesis work in the area of non-adiabatic dynamics. He continued at Oxford (D.Phil. 2009), working in quantum dynamics under the supervision of Professor David Manolopoulos. Together, the two developed an approach to allow quantum effects of nuclei to be included in condensed phase simulation at near classical computational cost, as well as elucidating isotope effects observed in liquids. Next, during postdoctoral work with Bruce Berne at Columbia University, Professor Markland focused on structure and dynamics in classical and quantum biophysical systems. He moved to Stanford in 2011 as an Assistant Professor in the Department of Chemistry and was promoted to Associate Professor with tenure in 2018. He has received recognition in a number of awards, including a Research Corporation Cottrell Scholarship, Alfred P. Sloan Research Fellowship, Terman Fellowship, Hellman Faculty Scholarship, the ACS OpenEye Outstanding Junior Faculty Award, the NSF CAREER award, the Camille Dreyfus Teacher-Scholar award, the H&S Dean's Award for Distinguished Teaching, the Kavli Emerging Leader in Chemistry Lectureship, and the ACS Early Career Award in Theoretical Chemistry.

Research in the Markland Group lies in the application and development of theoretical methods to model condensed phase systems, with a particular emphasis on the role of quantum mechanical effects. Treatment of these problems requires a range of theoretical approaches as well as molecular mechanics and ab initio simulations. The group is particularly interested in developing and applying methods based on the path integral formulation of quantum mechanics to include quantum fluctuations such as zero-point energy and tunneling in the dynamics of reactive condensed phase systems. The group has also developed methods to treat non-equilibrium excited state dynamics by exploiting the combination of quantum-classical theory and quantum master equation approaches.

Work in the Markland Group has already provided insights into several systems, including reactions in liquids and enzymes, and the quantum liquid–glass transition. Group members have also introduced methods to perform path integral calculations at near classical computational cost, expanding the ability to treat large-scale condensed phase systems.

Please visit the Markland Group website to learn more.

Academic Appointments

Administrative Appointments

  • Member, Stanford PULSE Institute, SLAC National Accelerator Laboratory (2014 - Present)

Honors & Awards

  • ACS Early Career Award in Theoretical Chemistry, American Chemical Society (2021)
  • Kavli Emerging Leader in Chemistry Lectureship, Kavli Foundation (2019)
  • Camille Dreyfus Teacher-Scholar Award, Camille and Henry Dreyfus Foundation (2017)
  • Terman Fellow, Stanford University (2017)
  • NSF CAREER award, National Science Foundation (2016)
  • Cottrell Scholar, Research Corporation for Science Advancement (2015)
  • Dean's Award for Distinguished Teaching, Stanford University School of Humanities & Sciences (2015)
  • Hellman Faculty Scholar, Stanford University (2014)
  • OpenEye Outstanding Junior Faculty Award, American Chemical Society (2014)
  • Sloan Research Fellowship, Alfred P. Sloan Foundation (2014)
  • Terman Fellow, Stanford University (2012)
  • Coulson Prize, Royal Society of Chemistry (2009)

Boards, Advisory Committees, Professional Organizations

  • General Member, Telluride Science Research Center (2015 - Present)

Professional Education

  • Postdoc, Columbia University, Theoretical Chemistry (2010)
  • DPhil, University of Oxford, Chemistry (2009)
  • MChem, University of Oxford, Chemistry (2006)

Current Research and Scholarly Interests

Our research centers on problems at the interface of quantum and statistical mechanics. Particular themes that occur frequently in our research are hydrogen bonding, the interplay between structure and dynamics, systems with multiple time and length-scales and quantum mechanical effects. The applications of our methods are diverse, ranging from chemistry to biology to geology and materials science. Particular current interests include proton and electron transfer in fuel cells and enzymatic systems, atmospheric isotope separation and the control of catalytic chemical reactivity using electric fields.

Treatment of these problems requires a range of analytic techniques as well as molecular mechanics and ab initio simulations. We are particularly interested in developing and applying methods based on the path integral formulation of quantum mechanics to include quantum fluctuations such as zero-point energy and tunneling in the dynamics of liquids and glasses. This formalism, in which a quantum mechanical particle is mapped onto a classical "ring polymer," provides an accurate and physically insightful way to calculate reaction rates, diffusion coefficients and spectra in systems containing light atoms. Our work has already provided intriguing insights in systems ranging from diffusion controlled reactions in liquids to the quantum liquid-glass transition as well as introducing methods to perform path integral calculations at near classical computational cost, expanding our ability to treat large-scale condensed phase systems.

2023-24 Courses

Stanford Advisees

All Publications

  • TorchMD-Net 2.0: Fast Neural Network Potentials for Molecular Simulations. ArXiv Pelaez, R. P., Simeon, G., Galvelis, R., Mirarchi, A., Eastman, P., Doerr, S., Thölke, P., Markland, T. E., De Fabritiis, G. 2024


    Achieving a balance between computational speed, prediction accuracy, and universal applicability in molecular simulations has been a persistent challenge. This paper presents substantial advancements in the TorchMD-Net software, a pivotal step forward in the shift from conventional force fields to neural network-based potentials. The evolution of TorchMD-Net into a more comprehensive and versatile framework is highlighted, incorporating cutting-edge architectures such as TensorNet. This transformation is achieved through a modular design approach, encouraging customized applications within the scientific community. The most notable enhancement is a significant improvement in computational efficiency, achieving a very remarkable acceleration in the computation of energy and forces for TensorNet models, with performance gains ranging from 2-fold to 10-fold over previous iterations. Other enhancements include highly optimized neighbor search algorithms that support periodic boundary conditions and the smooth integration with existing molecular dynamics frameworks. Additionally, the updated version introduces the capability to integrate physical priors, further enriching its application spectrum and utility in research. The software is available at

    View details for DOI 10.5281/zenodo.4774216

    View details for PubMedID 38463504

    View details for PubMedCentralID PMC10925388

  • Enhancing Protein-Ligand Binding Affinity Predictions Using Neural Network Potentials. Journal of chemical information and modeling Sabanés Zariquiey, F., Galvelis, R., Gallicchio, E., Chodera, J. D., Markland, T. E., De Fabritiis, G. 2024


    This letter gives results on improving protein-ligand binding affinity predictions based on molecular dynamics simulations using machine learning potentials with a hybrid neural network potential and molecular mechanics methodology (NNP/MM). We compute relative binding free energies with the Alchemical Transfer Method and validate its performance against established benchmarks and find significant enhancements compared with conventional MM force fields like GAFF2.

    View details for DOI 10.1021/acs.jcim.3c02031

    View details for PubMedID 38376463

  • OpenMM 8: Molecular Dynamics Simulation with Machine Learning Potentials. The journal of physical chemistry. B Eastman, P., Galvelis, R., Peláez, R. P., Abreu, C. R., Farr, S. E., Gallicchio, E., Gorenko, A., Henry, M. M., Hu, F., Huang, J., Krämer, A., Michel, J., Mitchell, J. A., Pande, V. S., Rodrigues, J. P., Rodriguez-Guerra, J., Simmonett, A. C., Singh, S., Swails, J., Turner, P., Wang, Y., Zhang, I., Chodera, J. D., De Fabritiis, G., Markland, T. E. 2023


    Machine learning plays an important and growing role in molecular simulation. The newest version of the OpenMM molecular dynamics toolkit introduces new features to support the use of machine learning potentials. Arbitrary PyTorch models can be added to a simulation and used to compute forces and energy. A higher-level interface allows users to easily model their molecules of interest with general purpose, pretrained potential functions. A collection of optimized CUDA kernels and custom PyTorch operations greatly improves the speed of simulations. We demonstrate these features in simulations of cyclin-dependent kinase 8 (CDK8) and the green fluorescent protein chromophore in water. Taken together, these features make it practical to use machine learning to improve the accuracy of simulations with only a modest increase in cost.

    View details for DOI 10.1021/acs.jpcb.3c06662

    View details for PubMedID 38154096

  • NNP/MM: Accelerating Molecular Dynamics Simulations with Machine Learning Potentials and Molecular Mechanics. Journal of chemical information and modeling Galvelis, R., Varela-Rial, A., Doerr, S., Fino, R., Eastman, P., Markland, T. E., Chodera, J. D., De Fabritiis, G. 2023


    Machine learning potentials have emerged as a means to enhance the accuracy of biomolecular simulations. However, their application is constrained by the significant computational cost arising from the vast number of parameters compared with traditional molecular mechanics. To tackle this issue, we introduce an optimized implementation of the hybrid method (NNP/MM), which combines a neural network potential (NNP) and molecular mechanics (MM). This approach models a portion of the system, such as a small molecule, using NNP while employing MM for the remaining system to boost efficiency. By conducting molecular dynamics (MD) simulations on various protein-ligand complexes and metadynamics (MTD) simulations on a ligand, we showcase the capabilities of our implementation of NNP/MM. It has enabled us to increase the simulation speed by ∼5 times and achieve a combined sampling of 1 μs for each complex, marking the longest simulations ever reported for this class of simulations.

    View details for DOI 10.1021/acs.jcim.3c00773

    View details for PubMedID 37694852

  • Developing machine-learned potentials to simultaneously capture the dynamics of excess protons and hydroxide ions in classical and path integral simulations. The Journal of chemical physics Atsango, A. O., Morawietz, T., Marsalek, O., Markland, T. E. 2023; 159 (7)


    The transport of excess protons and hydroxide ions in water underlies numerous important chemical and biological processes. Accurately simulating the associated transport mechanisms ideally requires utilizing abinitio molecular dynamics simulations to model the bond breaking and formation involved in proton transfer and path-integral simulations to model the nuclear quantum effects relevant to light hydrogen atoms. These requirements result in a prohibitive computational cost, especially at the time and length scales needed to converge proton transport properties. Here, we present machine-learned potentials (MLPs) that can model both excess protons and hydroxide ions at the generalized gradient approximation and hybrid density functional theory levels of accuracy and use them to perform multiple nanoseconds of both classical and path-integral proton defect simulations at a fraction of the cost of the corresponding abinitio simulations. We show that the MLPs are able to reproduce abinitio trends and converge properties such as the diffusion coefficients of both excess protons and hydroxide ions. We use our multi-nanosecond simulations, which allow us to monitor large numbers of proton transfer events, to analyze the role of hypercoordination in the transport mechanism of the hydroxide ion and provide further evidence for the asymmetry in diffusion between excess protons and hydroxide ions.

    View details for DOI 10.1063/5.0162066

    View details for PubMedID 37581418

  • Elucidating the Role of Hydrogen Bonding in the Optical Spectroscopy of the Solvated Green Fluorescent Protein Chromophore: Using Machine Learning to Establish the Importance of High-Level Electronic Structure. The journal of physical chemistry letters Chen, M. S., Mao, Y., Snider, A., Gupta, P., Montoya-Castillo, A., Zuehlsdorff, T. J., Isborn, C. M., Markland, T. E. 2023: 6610-6619


    Hydrogen bonding interactions with chromophores in chemical and biological environments play a key role in determining their electronic absorption and relaxation processes, which are manifested in their linear and multidimensional optical spectra. For chromophores in the condensed phase, the large number of atoms needed to simulate the environment has traditionally prohibited the use of high-level excited-state electronic structure methods. By leveraging transfer learning, we show how to construct machine-learned models to accurately predict the high-level excitation energies of a chromophore in solution from only 400 high-level calculations. We show that when the electronic excitations of the green fluorescent protein chromophore in water are treated using EOM-CCSD embedded in a DFT description of the solvent the optical spectrum is correctly captured and that this improvement arises from correctly treating the coupling of the electronic transition to electric fields, which leads to a larger response upon hydrogen bonding between the chromophore and water.

    View details for DOI 10.1021/acs.jpclett.3c01444

    View details for PubMedID 37459252

  • Electron transfer at electrode interfaces via a straightforward quasiclassical fermionic mapping approach. The Journal of chemical physics Jung, K. A., Kelly, J., Markland, T. E. 2023; 159 (1)


    Electron transfer at electrode interfaces to molecules in solution or at the electrode surface plays a vital role in numerous technological processes. However, treating these processes requires a unified and accurate treatment of the fermionic states of the electrode and their coupling to the molecule being oxidized or reduced in the electrochemical processes and, in turn, the way the molecular energy levels are modulated by the bosonic nuclear modes of the molecule and solvent. Here we present a physically transparent quasiclassical scheme to treat these electrochemical electron transfer processes in the presence of molecular vibrations by using an appropriately chosen mapping of the fermionic variables. We demonstrate that this approach, which is exact in the limit of non-interacting fermions in the absence of coupling to vibrations, is able to accurately capture the electron transfer dynamics from the electrode even when the process is coupled to vibrational motions in the regimes of weak coupling. This approach, thus, provides a scalable strategy to explicitly treat electron transfer from electrode interfaces in condensed-phase molecular systems.

    View details for DOI 10.1063/5.0156136

    View details for PubMedID 37409707

  • Building insightful, memory-enriched models to capture long-time biochemical processes from short-time simulations. Proceedings of the National Academy of Sciences of the United States of America Dominic, A. J., Sayer, T., Cao, S., Markland, T. E., Huang, X., Montoya-Castillo, A. 2023; 120 (12): e2221048120


    The ability to predict and understand complex molecular motions occurring over diverse timescales ranging from picoseconds to seconds and even hours in biological systems remains one of the largest challenges to chemical theory. Markov state models (MSMs), which provide a memoryless description of the transitions between different states of a biochemical system, have provided numerous important physically transparent insights into biological function. However, constructing these models often necessitates performing extremely long molecular simulations to converge the rates. Here, we show that by incorporating memory via the time-convolutionless generalized master equation (TCL-GME) one can build a theoretically transparent and physically intuitive memory-enriched model of biochemical processes with up to a three order of magnitude reduction in the simulation data required while also providing a higher temporal resolution. We derive the conditions under which the TCL-GME provides a more efficient means to capture slow dynamics than MSMs and rigorously prove when the two provide equally valid and efficient descriptions of the slow configurational dynamics. We further introduce a simple averaging procedure that enables our TCL-GME approach to quickly converge and accurately predict long-time dynamics even when parameterized with noisy reference data arising from short trajectories. We illustrate the advantages of the TCL-GME using alanine dipeptide, the human argonaute complex, and FiP35 WW domain.

    View details for DOI 10.1073/pnas.2221048120

    View details for PubMedID 36920924

  • A derivation of the conditions under which bosonic operators exactly capture fermionic structure and dynamics. The Journal of chemical physics Montoya-Castillo, A., Markland, T. E. 2023; 158 (9): 094112


    The dynamics of many-body fermionic systems are important in problems ranging from catalytic reactions at electrochemical surfaces to transport through nanojunctions and offer a prime target for quantum computing applications. Here, we derive the set of conditions under which fermionic operators can be exactly replaced by bosonic operators that render the problem amenable to a large toolbox of dynamical methods while still capturing the correct dynamics of n-body operators. Importantly, our analysis offers a simple guide on how one can exploit these simple maps to calculate nonequilibrium and equilibrium single- and multi-time correlation functions essential in describing transport and spectroscopy. We use this to rigorously analyze and delineate the applicability of simple yet effective Cartesian maps that have been shown to correctly capture the correct fermionic dynamics in select models of nanoscopic transport. We illustrate our analytical results with exact simulations of the resonant level model. Our work provides new insights as to when one can leverage the simplicity of bosonic maps to simulate the dynamics of many-electron systems, especially those where an atomistic representation of nuclear interactions becomes essential.

    View details for DOI 10.1063/5.0138664

    View details for PubMedID 36889969

  • An accurate and efficient Ehrenfest dynamics approach for calculating linear and nonlinear electronic spectra. The Journal of chemical physics Atsango, A. O., Montoya-Castillo, A., Markland, T. E. 2023; 158 (7): 074107


    Linear and nonlinear electronic spectra provide an important tool to probe the absorption and transfer of electronic energy. Here, we introduce a pure state Ehrenfest approach to obtain accurate linear and nonlinear spectra that is applicable to systems with large numbers of excited states and complex chemical environments. We achieve this by representing the initial conditions as sums of pure states and unfolding multi-time correlation functions into the Schrodinger picture. By doing this, we show that one can obtain significant improvements in accuracy over the previously used projected Ehrenfest approach and that these benefits are particularly pronounced in cases where the initial condition is a coherence between excited states. While such initial conditions do not arise when calculating linear electronic spectra, they play a vital role in capturing multidimensional spectroscopies. We demonstrate the performance of our method by showing that it is able to quantitatively capture the exact linear, 2D electronic spectroscopy, and pump-probe spectra for a Frenkel exciton model in slow bath regimes and is even able to reproduce the main spectral features in fast bath regimes.

    View details for DOI 10.1063/5.0138671

    View details for PubMedID 36813724

  • Data-Efficient Machine Learning Potentials from Transfer Learning of Periodic Correlated Electronic Structure Methods: Liquid Water at AFQMC, CCSD, and CCSD(T) Accuracy. Journal of chemical theory and computation Chen, M. S., Lee, J., Ye, H. Z., Berkelbach, T. C., Reichman, D. R., Markland, T. E. 2023


    Obtaining the atomistic structure and dynamics of disordered condensed-phase systems from first-principles remains one of the forefront challenges of chemical theory. Here we exploit recent advances in periodic electronic structure and provide a data-efficient approach to obtain machine-learned condensed-phase potential energy surfaces using AFQMC, CCSD, and CCSD(T) from a very small number (≤200) of energies by leveraging a transfer learning scheme starting from lower-tier electronic structure methods. We demonstrate the effectiveness of this approach for liquid water by performing both classical and path integral molecular dynamics simulations on these machine-learned potential energy surfaces. By doing this, we uncover the interplay of dynamical electron correlation and nuclear quantum effects across the entire liquid range of water while providing a general strategy for efficiently utilizing periodic correlated electronic structure methods to explore disordered condensed-phase systems.

    View details for DOI 10.1021/acs.jctc.2c01203

    View details for PubMedID 36730728

  • SPICE, A Dataset of Drug-like Molecules and Peptides for Training Machine Learning Potentials. Scientific data Eastman, P., Behara, P. K., Dotson, D. L., Galvelis, R., Herr, J. E., Horton, J. T., Mao, Y., Chodera, J. D., Pritchard, B. P., Wang, Y., De Fabritiis, G., Markland, T. E. 2023; 10 (1): 11


    Machine learning potentials are an important tool for molecular simulation, but their development is held back by a shortage of high quality datasets to train them on. We describe the SPICE dataset, a new quantum chemistry dataset for training potentials relevant to simulating drug-like small molecules interacting with proteins. It contains over 1.1 million conformations for a diverse set of small molecules, dimers, dipeptides, and solvated amino acids. It includes 15 elements, charged and uncharged molecules, and a wide range of covalent and non-covalent interactions. It provides both forces and energies calculated at the ωB97M-D3(BJ)/def2-TZVPPD level of theory, along with other useful quantities such as multipole moments and bond orders. We train a set of machine learning potentials on it and demonstrate that they can achieve chemical accuracy across a broad region of chemical space. It can serve as a valuable resource for the creation of transferable, ready to use potential functions for use in molecular simulations.

    View details for DOI 10.1038/s41597-022-01882-6

    View details for PubMedID 36599873

  • 2D spectroscopies from condensed phase dynamics: Accessing third-order response properties from equilibrium multi-time correlation functions. The Journal of chemical physics Jung, K. A., Markland, T. E. 2022; 157 (9): 094111


    The third-order response lies at the heart of simulating and interpreting nonlinear spectroscopies ranging from two-dimensional infrared (2D-IR) to 2D electronic (2D-ES), and 2D sum frequency generation (2D-SFG). The extra time and frequency dimensions in these spectroscopic techniques provide access to rich information on the electronic and vibrational states present, the coupling between them, and the resulting rates at which they exchange energy that are obscured in linear spectroscopy, particularly for condensed phase systems that usually contain many overlapping features. While the exact quantum expression for the third-order response is well established, it is incompatible with the methods that are practical for calculating the atomistic dynamics of large condensed phase systems. These methods, which include both classical mechanics and quantum dynamics methods that retain quantum statistical properties while obeying the symmetries of classical dynamics, such as LSC-IVR, centroid molecular dynamics, and Ring Polymer Molecular Dynamics (RPMD), naturally provide short-time approximations to the multi-time symmetrized Kubo transformed correlation function. Here, we show how the third-order response can be formulated in terms of equilibrium symmetrized Kubo transformed correlation functions. We demonstrate the utility and accuracy of our approach by showing how it can be used to obtain the third-order response of a series of model systems using both classical dynamics and RPMD. In particular, we show that this approach captures features such as anharmonically induced vertical splittings and peak shifts while providing a physically transparent framework for understanding multidimensional spectroscopies.

    View details for DOI 10.1063/5.0107087

    View details for PubMedID 36075710

  • Optically Induced Anisotropy in Time-Resolved Scattering: Imaging Molecular-Scale Structure and Dynamics in Disordered Media with Experiment and Theory. Physical review letters Montoya-Castillo, A., Chen, M. S., Raj, S. L., Jung, K. A., Kjaer, K. S., Morawietz, T., Gaffney, K. J., van Driel, T. B., Markland, T. E. 2022; 129 (5): 056001


    Time-resolved scattering experiments enable imaging of materials at the molecular scale with femtosecond time resolution. However, in disordered media they provide access to just one radial dimension thus limiting the study of orientational structure and dynamics. Here we introduce a rigorous and practical theoretical framework for predicting and interpreting experiments combining optically induced anisotropy and time-resolved scattering. Using impulsive nuclear Raman and ultrafast x-ray scattering experiments of chloroform and simulations, we demonstrate that this framework can accurately predict and elucidate both the spatial and temporal features of these experiments.

    View details for DOI 10.1103/PhysRevLett.129.056001

    View details for PubMedID 35960558

  • Solvent Organization and Electrostatics Tuned by Solute Electronic Structure: Amide versus Non-Amide Carbonyls. The journal of physical chemistry. B Fried, S. D., Zheng, C., Mao, Y., Markland, T. E., Boxer, S. G. 2022


    The ability to exploit carbonyl groups to measure electric fields in enzymes and other complex reactive environments by using the vibrational Stark effect has inspired growing interest in how these fields can be measured, tuned, and ultimately designed. Previous studies have concentrated on the role of the solvent in tuning the fields exerted on the solute. Here, we explore instead the role of the solute electronic structure in modifying the local solvent organization and electric field exerted on the solute. By measuring the infrared absorption spectra of amide-containing molecules, as prototypical peptides, and contrasting them with non-amide carbonyls in a wide range of solvents, we show that these solutes experience notable differences in their frequency shifts in polar solvents. Using vibrational Stark spectroscopy and molecular dynamics simulations, we demonstrate that while some of these differences can be rationalized by using the distinct intrinsic Stark tuning rates of the solutes, the larger frequency shifts for amides and dimethylurea primarily result from the larger solvent electric fields experienced by their carbonyl groups. These larger fields arise due to their stronger p-π conjugation, which results in larger C═O bond dipole moments that further induce substantial solvent organization. Using electronic structure calculations, we decompose the electric fields into contributions from solvent molecules that are in the first solvation shell and those from the bulk and show that both of these contributions are significant and become larger with enhanced conjugation in solutes. These results show that structural modifications of a solute can be used to tune both the solvent organization and electrostatic environment, indicating the importance of a solute-centric paradigm in modulating and designing the electrostatic environment in condensed-phase chemical processes.

    View details for DOI 10.1021/acs.jpcb.2c03095

    View details for PubMedID 35901512

  • A two-directional vibrational probe reveals different electric field orientations in solution and an enzyme active site. Nature chemistry Zheng, C., Mao, Y., Kozuch, J., Atsango, A. O., Ji, Z., Markland, T. E., Boxer, S. G. 2022


    The catalytic power of an electric field depends on its magnitude and orientation with respect to the reactive chemical species. Understanding and designing new catalysts for electrostatic catalysis thus requires methods to measure the electric field orientation and magnitude at the molecular scale. We demonstrate that electric field orientations can be extracted using a two-directional vibrational probe by exploiting the vibrational Stark effect of both the C=O and C-D stretches of a deuterated aldehyde. Combining spectroscopy with molecular dynamics and electronic structure partitioning methods, we demonstrate that, despite distinct polarities, solvents act similarly in their preference for electrostatically stabilizing large bond dipoles at the expense of destabilizing small ones. In contrast, we find that for an active-site aldehyde inhibitor of liver alcohol dehydrogenase, the electric field orientation deviates markedly from that found in solvents, which provides direct evidence for the fundamental difference between the electrostatic environment of solvents and that of a preorganized enzyme active site.

    View details for DOI 10.1038/s41557-022-00937-w

    View details for PubMedID 35513508

  • A framework for automated structure elucidation from routine NMR spectra. Chemical science Huang, Z., Chen, M. S., Woroch, C. P., Markland, T. E., Kanan, M. W. 2021; 12 (46): 15329-15338


    Methods to automate structure elucidation that can be applied broadly across chemical structure space have the potential to greatly accelerate chemical discovery. NMR spectroscopy is the most widely used and arguably the most powerful method for elucidating structures of organic molecules. Here we introduce a machine learning (ML) framework that provides a quantitative probabilistic ranking of the most likely structural connectivity of an unknown compound when given routine, experimental one dimensional 1H and/or 13C NMR spectra. In particular, our ML-based algorithm takes input NMR spectra and (i) predicts the presence of specific substructures out of hundreds of substructures it has learned to identify; (ii) annotates the spectrum to label peaks with predicted substructures; and (iii) uses the substructures to construct candidate constitutional isomers and assign to them a probabilistic ranking. Using experimental spectra and molecular formulae for molecules containing up to 10 non-hydrogen atoms, the correct constitutional isomer was the highest-ranking prediction made by our model in 67.4% of the cases and one of the top-ten predictions in 95.8% of the cases. This advance will aid in solving the structure of unknown compounds, and thus further the development of automated structure elucidation tools that could enable the creation of fully autonomous reaction discovery platforms.

    View details for DOI 10.1039/d1sc04105c

    View details for PubMedID 34976353

    View details for PubMedCentralID PMC8635205

  • A framework for automated structure elucidation from routine NMR spectra CHEMICAL SCIENCE Huang, Z., Chen, M. S., Woroch, C. P., Markland, T. E., Kanan, M. W. 2021

    View details for DOI 10.1039/d1sc04105c

    View details for Web of Science ID 000718969800001

  • Characterizing and Contrasting Structural Proton Transport Mechanisms in Azole Hydrogen Bond Networks Using Ab Initio Molecular Dynamics. The journal of physical chemistry letters Atsango, A. O., Tuckerman, M. E., Markland, T. E. 2021: 8749-8756


    Imidazole and 1,2,3-triazole are promising hydrogen-bonded heterocycles that conduct protons via a structural mechanism and whose derivatives are present in systems ranging from biological proton channels to proton exchange membrane fuel cells. Here, we leverage multiple time-stepping to perform ab initio molecular dynamics of imidazole and 1,2,3-triazole at the nanosecond time scale. We show that despite the close structural similarities of these compounds, their proton diffusion constants vary by over an order of magnitude. Our simulations reveal the reasons for these differences in diffusion constants, which range from the degree of hydrogen-bonded chain linearity to the effect of the central nitrogen atom in 1,2,3-triazole on proton transport. In particular, we uncover evidence of two "blocking" mechanisms in 1,2,3-triazole, where covalent and hydrogen bonds formed by the central nitrogen atom limit the mobility of protons. Our simulations thus provide insights into the origins of the experimentally observed 10-fold difference in proton conductivity.

    View details for DOI 10.1021/acs.jpclett.1c02266

    View details for PubMedID 34478302

  • AENET-LAMMPS and AENET-TINKER: Interfaces for accurate and efficient molecular dynamics simulations with machine learning potentials JOURNAL OF CHEMICAL PHYSICS Chen, M. S., Morawietz, T., Mori, H., Markland, T. E., Artrith, N. 2021; 155 (7): 074801


    Machine-learning potentials (MLPs) trained on data from quantum-mechanics based first-principles methods can approach the accuracy of the reference method at a fraction of the computational cost. To facilitate efficient MLP-based molecular dynamics and Monte Carlo simulations, an integration of the MLPs with sampling software is needed. Here, we develop two interfaces that link the atomic energy network (ænet) MLP package with the popular sampling packages TINKER and LAMMPS. The three packages, ænet, TINKER, and LAMMPS, are free and open-source software that enable, in combination, accurate simulations of large and complex systems with low computational cost that scales linearly with the number of atoms. Scaling tests show that the parallel efficiency of the ænet-TINKER interface is nearly optimal but is limited to shared-memory systems. The ænet-LAMMPS interface achieves excellent parallel efficiency on highly parallel distributed-memory systems and benefits from the highly optimized neighbor list implemented in LAMMPS. We demonstrate the utility of the two MLP interfaces for two relevant example applications: the investigation of diffusion phenomena in liquid water and the equilibration of nanostructured amorphous battery materials.

    View details for DOI 10.1063/5.0063880

    View details for Web of Science ID 000685163500002

    View details for PubMedID 34418919

  • Persistent Homology Metrics Reveal Quantum Fluctuations and Reactive Atoms in Path Integral Dynamics. Frontiers in chemistry Hu, Y., Ounkham, P., Marsalek, O., Markland, T. E., Krishmoorthy, B., Clark, A. E. 2021; 9: 624937


    Nuclear quantum effects (NQEs) are known to impact a number of features associated with chemical reactivity and physicochemical properties, particularly for light atoms and at low temperatures. In the imaginary time path integral formalism, each atom is mapped onto a "ring polymer" whose spread is related to the quantum mechanical uncertainty in the particle's position, i.e., its thermal wavelength. A number of metrics have previously been used to investigate and characterize this spread and explain effects arising from quantum delocalization, zero-point energy, and tunneling. Many of these shape metrics consider just the instantaneous structure of the ring polymers. However, given the significant interest in methods such as centroid molecular dynamics and ring polymer molecular dynamics that link the molecular dynamics of these ring polymers to real time properties, there exists significant opportunity to exploit metrics that also allow for the study of the fluctuations of the atom delocalization in time. Here we consider the ring polymer delocalization from the perspective of computational topology, specifically persistent homology, which describes the 3-dimensional arrangement of point cloud data, (i.e. atomic positions). We employ the Betti sequence probability distribution to define the ensemble of shapes adopted by the ring polymer. The Wasserstein distances of Betti sequences adjacent in time are used to characterize fluctuations in shape, where the Fourier transform and associated principal components provides added information differentiating atoms with different NQEs based on their dynamic properties. We demonstrate this methodology on two representative systems, a glassy system consisting of two atom types with dramatically different de Broglie thermal wavelengths, and ab initio molecular dynamics simulation of an aqueous 4M HCl solution where the H-atoms are differentiated based on their participation in proton transfer reactions.

    View details for DOI 10.3389/fchem.2021.624937

    View details for PubMedID 33748074

  • Exploiting Machine Learning to Efficiently Predict Multidimensional Optical Spectra in Complex Environments. The journal of physical chemistry letters Chen, M. S., Zuehlsdorff, T. J., Morawietz, T., Isborn, C. M., Markland, T. E. 2020: 7559–68


    The excited-state dynamics of chromophores in complex environments determine a range of vital biological and energy capture processes. Time-resolved, multidimensional optical spectroscopies provide a key tool to investigate these processes. Although theory has the potential to decode these spectra in terms of the electronic and atomistic dynamics, the need for large numbers of excited-state electronic structure calculations severely limits first-principles predictions of multidimensional optical spectra for chromophores in the condensed phase. Here, we leverage the locality of chromophore excitations to develop machine learning models to predict the excited-state energy gap of chromophores in complex environments for efficiently constructing linear and multidimensional optical spectra. By analyzing the performance of these models, which span a hierarchy of physical approximations, across a range of chromophore-environment interaction strengths, we provide strategies for the construction of machine learning models that greatly accelerate the calculation of multidimensional optical spectra from first principles.

    View details for DOI 10.1021/acs.jpclett.0c02168

    View details for PubMedID 32808797

  • On the advantages of exploiting memory in Markov state models for biomolecular dynamics. The Journal of chemical physics Cao, S., Montoya-Castillo, A., Wang, W., Markland, T. E., Huang, X. 2020; 153 (1): 014105


    Biomolecular dynamics play an important role in numerous biological processes. Markov State Models (MSMs) provide a powerful approach to study these dynamic processes by predicting long time scale dynamics based on many short molecular dynamics (MD) simulations. In an MSM, protein dynamics are modeled as a kinetic process consisting of a series of Markovian transitions between different conformational states at discrete time intervals (called "lag time"). To achieve this, a master equation must be constructed with a sufficiently long lag time to allow interstate transitions to become truly Markovian. This imposes a major challenge for MSM studies of proteins since the lag time is bound by the length of relatively short MD simulations available to estimate the frequency of transitions. Here, we show how one can employ the generalized master equation formalism to obtain an exact description of protein conformational dynamics both at short and long time scales without the time resolution restrictions imposed by the MSM lag time. Using a simple kinetic model, alanine dipeptide, and WW domain, we demonstrate that it is possible to construct these quasi-Markov State Models (qMSMs) using MD simulations that are 5-10 times shorter than those required by MSMs. These qMSMs only contain a handful of metastable states and, thus, can greatly facilitate the interpretation of mechanisms associated with protein dynamics. A qMSM opens the door to the study of conformational changes of complex biomolecules where a Markovian model with a few states is often difficult to construct due to the limited length of available MD simulations.

    View details for DOI 10.1063/5.0010787

    View details for PubMedID 32640825

  • Quantum kinetic energy and isotope fractionation in aqueous ionic solutions. Physical chemistry chemical physics : PCCP Wang, L. n., Ceriotti, M. n., Markland, T. E. 2020


    At room temperature, the quantum contribution to the kinetic energy of a water molecule exceeds the classical contribution by an order of magnitude. The quantum kinetic energy (QKE) of a water molecule is modulated by its local chemical environment and leads to uneven partitioning of isotopes between different phases in thermal equilibrium, which would not occur if the nuclei behaved classically. In this work, we use ab initio path integral simulations to show that QKEs of the water molecules and the equilibrium isotope fractionation ratios of the oxygen and hydrogen isotopes are sensitive probes of the hydrogen bonding structures in aqueous ionic solutions. In particular, we demonstrate how the QKE of water molecules in path integral simulations can be decomposed into translational, rotational and vibrational degrees of freedom, and use them to determine the impact of solvation on different molecular motions. By analyzing the QKEs and isotope fractionation ratios, we show how the addition of the Na+, Cl- and HPO42- ions perturbs the competition between quantum effects in liquid water and impacts their local solvation structures.

    View details for DOI 10.1039/c9cp06483d

    View details for PubMedID 31942581

  • Elucidating the Proton Transport Pathways in Liquid Imidazole with First-Principles Molecular Dynamics. The journal of physical chemistry letters Long, Z. n., Atsango, A. O., Napoli, J. A., Markland, T. E., Tuckerman, M. E. 2020: 6156–63


    Imidazole is a promising anhydrous proton conductor with a high conductivity comparable to that of water at a similar temperature relative to its melting point. Previous theoretical studies of the mechanism of proton transport in imidazole have relied either on empirical models or on ab initio trajectories that have been too short to draw significant conclusions. Here, we present the results of multiple time-step ab initio molecular dynamics simulations of an excess proton in liquid imidazole reaching 1 ns in total simulation time. We find that the proton transport is dominated by structural diffusion, with the diffusion constant of the proton defect being ∼8 times higher than that of self-diffusion of the imidazole molecules. By using correlation function analysis, we decompose the mechanism for proton transport into a series of first-order processes and show that the proton transport mechanism occurs over three distinct time and length scales. Although the mechanism at intermediate times is dominated by hopping along pseudo-one-dimensional chains, at longer times the overall rate of diffusion is limited by the re-formation of these chains. These results provide a more complete picture of the traditional idealized Grotthuss structural diffusion mechanism.

    View details for DOI 10.1021/acs.jpclett.0c01744

    View details for PubMedID 32633523

  • Excited state diabatization on the cheap using DFT: Photoinduced electron and hole transfer. The Journal of chemical physics Mao, Y. n., Montoya-Castillo, A. n., Markland, T. E. 2020; 153 (24): 244111


    Excited state electron and hole transfer underpin fundamental steps in processes such as exciton dissociation at photovoltaic heterojunctions, photoinduced charge transfer at electrodes, and electron transfer in photosynthetic reaction centers. Diabatic states corresponding to charge or excitation localized species, such as locally excited and charge transfer states, provide a physically intuitive framework to simulate and understand these processes. However, obtaining accurate diabatic states and their couplings from adiabatic electronic states generally leads to inaccurate results when combined with low-tier electronic structure methods, such as time-dependent density functional theory, and exorbitant computational cost when combined with high-level wavefunction-based methods. Here, we introduce a density functional theory (DFT)-based diabatization scheme that directly constructs the diabatic states using absolutely localized molecular orbitals (ALMOs), which we denote as Δ-ALMO(MSDFT2). We demonstrate that our method, which combines ALMO calculations with the ΔSCF technique to construct electronically excited diabatic states and obtains their couplings with charge-transfer states using our MSDFT2 scheme, gives accurate results for excited state electron and hole transfer in both charged and uncharged systems that underlie DNA repair, charge separation in donor-acceptor dyads, chromophore-to-solvent electron transfer, and singlet fission. This framework for the accurate and efficient construction of excited state diabats and evaluation of their couplings directly from DFT thus offers a route to simulate and elucidate photoinduced electron and hole transfer in large disordered systems, such as those encountered in the condensed phase.

    View details for DOI 10.1063/5.0035593

    View details for PubMedID 33380087

  • Resolving Heterogeneous Dynamics of Excess Protons in Aqueous Solution with Rate Theory. The journal of physical chemistry. B Roy, S. n., Schenter, G. K., Napoli, J. A., Baer, M. D., Markland, T. E., Mundy, C. J. 2020


    Rate theories have found great utility across the chemical sciences by providing a physically transparent way to analyze dynamical processes. Here we demonstrate the benefits of using transition state theory and Marcus theory to study the rate of proton transfer in HCl solutions. By using long ab initio molecular dynamics simulations, we show that good agreement is obtained between these two different formulations of rate theory and how they can be used to study the pathways and lifetime of proton transfer in aqueous solution. Since both rate theory formulations utilize identical sets of molecular data, this provides a self-consistent theoretical picture of the rates of aqueous phase proton transfer. Specifically, we isolate and quantify the rates of proton transfer, ion-pair dissociation, and solvent exchange in concentrated HCl solutions. Our analysis predicts a concentration dependence to both proton transfer and ion-pairing. Moreover, our estimate of the lifetime for the Zundel species is 0.8 and 1.3 ps for 2 M and 8 M HCl, respectively. We demonstrate that concentration effects can indeed be quantified through the combination of state-of-the-art simulation and theory and provides a picture of the important correlations between the cation (hydronium) and the counterion in acid solutions.

    View details for DOI 10.1021/acs.jpcb.0c02649

    View details for PubMedID 32482074

  • Accurate and efficient DFT-based diabatization for hole and electron transfer using absolutely localized molecular orbitals. The Journal of chemical physics Mao, Y., Montoya-Castillo, A., Markland, T. E. 2019; 151 (16): 164114


    Diabatic states and the couplings between them are important for quantifying, elucidating, and predicting the rates and mechanisms of many chemical and biochemical processes. Here, we propose and investigate approaches to accurately compute diabatic couplings from density functional theory (DFT) using absolutely localized molecular orbitals (ALMOs). ALMOs provide an appealing approach to generate variationally optimized diabatic states and obtain their associated forces, which allows for the relaxation of the donor and acceptor orbitals in a way that is internally consistent in how the method treats both the donor and acceptor states. Here, we show that one can obtain more accurate electronic couplings between ALMO-based diabats by employing the symmetrized transition density matrix to evaluate the exchange-correlation contribution. We demonstrate that this approach yields accurate results in comparison to other commonly used DFT-based diabatization methods across a wide array of electron and hole transfer processes occurring in systems ranging from conjugated organic molecules, such as thiophene and pentacene, to DNA base pairs. We also show that this approach yields accurate diabatic couplings even when combined with lower tiers of the DFT hierarchy, opening the door to combining it with quantum dynamics approaches to provide an ab initio treatment of nonadiabatic processes in the condensed phase.

    View details for DOI 10.1063/1.5125275

    View details for PubMedID 31675855

  • Tracking Aqueous Proton Transfer by Two-Dimensional Infrared Spectroscopy and ab Initio Molecular Dynamics Simulations ACS CENTRAL SCIENCE Yuan, R., Napoli, J. A., Yan, C., Marsalek, O., Markland, T. E., Fayer, M. D. 2019; 5 (7): 1269–77
  • Tracking Aqueous Proton Transfer by Two-Dimensional Infrared Spectroscopy and ab Initio Molecular Dynamics Simulations. ACS central science Yuan, R., Napoli, J. A., Yan, C., Marsalek, O., Markland, T. E., Fayer, M. D. 2019; 5 (7): 1269-1277


    Proton transfer in water is ubiquitous and a critical elementary event that, via proton hopping between water molecules, enables protons to diffuse much faster than other ions. The problem of the anomalous nature of proton transport in water was first identified by Grotthuss over 200 years ago. In spite of a vast amount of modern research effort, there are still many unanswered questions about proton transport in water. An experimental determination of the proton hopping time has remained elusive due to its ultrafast nature and the lack of direct experimental observables. Here, we use two-dimensional infrared spectroscopy to extract the chemical exchange rates between hydronium and water in acid solutions using a vibrational probe, methyl thiocyanate. Ab initio molecular dynamics (AIMD) simulations demonstrate that the chemical exchange is dominated by proton hopping. The observed experimental and simulated acid concentration dependence then allow us to extrapolate the measured single step proton hopping time to the dilute limit, which, within error, gives the same value as inferred from measurements of the proton mobility and NMR line width analysis. In addition to obtaining the proton hopping time in the dilute limit from direct measurements and AIMD simulations, the results indicate that proton hopping in dilute acid solutions is induced by the concerted multi-water molecule hydrogen bond rearrangement that occurs in pure water. This proposition on the dynamics that drive proton hopping is confirmed by a combination of experimental results from the literature.

    View details for DOI 10.1021/acscentsci.9b00447

    View details for PubMedID 31403075

    View details for PubMedCentralID PMC6661862

  • i-PI 2.0: A universal force engine for advanced molecular simulations COMPUTER PHYSICS COMMUNICATIONS Kapil, V., Rossi, M., Marsalek, O., Petraglia, R., Litman, Y., Spura, T., Cheng, B., Cuzzocrea, A., Meissner, R. H., Wilkins, D. M., Helfrecht, B. A., Juda, P., Bienvenue, S. P., Fang, W., Kessler, J., Poltavsky, I., Vandenbrande, S., Wieme, J., Corminboeuf, C., Kuehne, T. D., Manolopoulos, D. E., Markland, T. E., Richardson, J. O., Tkatchenko, A., Tribello, G. A., Van Speybroeck, V., Ceriotti, M. 2019; 236: 214–23
  • Efficient construction of generalized master equation memory kernels for multi-state systems from nonadiabatic quantum-classical dynamics. The Journal of chemical physics Pfalzgraff, W. C., Montoya-Castillo, A. n., Kelly, A. n., Markland, T. E. 2019; 150 (24): 244109


    Methods derived from the generalized quantum master equation (GQME) framework have provided the basis for elucidating energy and charge transfer in systems ranging from molecular solids to photosynthetic complexes. Recently, the nonperturbative combination of the GQME with quantum-classical methods has resulted in approaches whose accuracy and efficiency exceed those of the original quantum-classical schemes while offering significant accuracy improvements over perturbative expansions of the GQME. Here, we show that, while the non-Markovian memory kernel required to propagate the GQME scales quartically with the number of subsystem states, the number of trajectories required scales at most quadratically when using quantum-classical methods to construct the kernel. We then present an algorithm that allows further acceleration of the quantum-classical GQME by providing a way to selectively sample the kernel matrix elements that are most important to the process of interest. We demonstrate the utility of these advances by applying the combination of Ehrenfest mean field theory with the GQME (MF-GQME) to models of the Fenna-Matthews-Olson (FMO) complex and the light harvesting complex II (LHCII), with 7 and 14 states, respectively. This allows us to show that the MF-GQME is able to accurately capture all the relevant dynamical time scales in LHCII: the initial nonequilibrium population transfer on the femtosecond time scale, the steady state-type trapping on the picosecond time scale, and the long time population relaxation. Remarkably, all of these physical effects spanning tens of picoseconds can be encoded in a memory kernel that decays only after ∼65 fs.

    View details for DOI 10.1063/1.5095715

    View details for PubMedID 31255061

  • Hiding in the Crowd: Spectral Signatures of Overcoordinated Hydrogen-Bond Environments. The journal of physical chemistry letters Morawietz, T. n., Urbina, A. S., Wise, P. K., Wu, X. n., Lu, W. n., Ben-Amotz, D. n., Markland, T. E. 2019: 6067–73


    Molecules with an excess number of hydrogen-bonding partners play a crucial role in fundamental chemical processes, ranging from anomalous diffusion in supercooled water to transport of aqueous proton defects and ordering of water around hydrophobic solutes. Here we show that overcoordinated hydrogen-bond environments can be identified in both the ambient and supercooled regimes of liquid water by combining experimental Raman multivariate curve resolution measurements and machine learning accelerated quantum simulations. In particular, we find that OH groups appearing in spectral regions usually associated with non-hydrogen-bonded species actually correspond to hydrogen bonds formed in overcoordinated environments. We further show that only these species exhibit a turnover in population as a function of temperature, which is robust and persists under both constant pressure and density conditions. This work thus provides a new tool to identify, interpret, and elucidate the spectral signatures of crowded hydrogen-bond networks.

    View details for DOI 10.1021/acs.jpclett.9b01781

    View details for PubMedID 31549833

  • Optical spectra in the condensed phase: Capturing anharmonic and vibronic features using dynamic and static approaches. The Journal of chemical physics Zuehlsdorff, T. J., Montoya-Castillo, A. n., Napoli, J. A., Markland, T. E., Isborn, C. M. 2019; 151 (7): 074111


    Simulating optical spectra in the condensed phase remains a challenge for theory due to the need to capture spectral signatures arising from anharmonicity and dynamical effects, such as vibronic progressions and asymmetry. As such, numerous simulation methods have been developed that invoke different approximations and vary in their ability to capture different physical regimes. Here, we use several models of chromophores in the condensed phase and ab initio molecular dynamics simulations to rigorously assess the applicability of methods to simulate optical absorption spectra. Specifically, we focus on the ensemble scheme, which can address anharmonic potential energy surfaces but relies on the applicability of extreme nuclear-electronic time scale separation; the Franck-Condon method, which includes dynamical effects but generally only at the harmonic level; and the recently introduced ensemble zero-temperature Franck-Condon approach, which straddles these limits. We also devote particular attention to the performance of methods derived from a cumulant expansion of the energy gap fluctuations and test the ability to approximate the requisite time correlation functions using classical dynamics with quantum correction factors. These results provide insights as to when these methods are applicable and able to capture the features of condensed phase spectra qualitatively and, in some cases, quantitatively across a range of regimes.

    View details for DOI 10.1063/1.5114818

    View details for PubMedID 31438704

  • Beyond Badger's Rule: The Origins and Generality of the Structure-Spectra Relationship of Aqueous Hydrogen Bonds. The journal of physical chemistry letters Boyer, M. A., Marsalek, O. n., Heindel, J. P., Markland, T. E., McCoy, A. B., Xantheas, S. S. 2019: 918–24


    The structure of hydrogen bonded networks is intimately intertwined with their dynamics. Despite the incredibly wide range of hydrogen bond strengths encountered in water clusters, ion-water clusters, and liquid water, we demonstrate that the previously reported correlation between the change in the equilibrium bond length of the hydrogen bonded OH covalent bond and the corresponding shift in its harmonic frequency in water clusters is much more broadly applicable. Surprisingly, this correlation describes the ratios for both the equilibrium OH bond length/harmonic frequency and the vibrationally averaged bond length/anharmonic frequency in water, hydronium water, and halide water clusters. Consideration of harmonic and anaharmonic data leads to a correlation of -19 ± 1 cm-1/0.001 Å. The fundamental nature of this correlation is further confirmed through the analysis of ab initio Molecular Dynamics (AIMD) trajectories for liquid water. We demonstrate that this simple correlation for both harmonic and anharmonic systems can be modeled by the response of an OH bond to an external field. Treating the OH bond as a Morse oscillator, we develop analytic expressions, which relate the ratio of the shift in the vibrational frequency of the hydrogen-bonded OH bond to the shift in OH bond length, to parameters in the Morse potential and the ratio of the first and second derivatives of the field-dependent projection of the dipole moment of water onto the hydrogen-bonded OH bond. Based on our analysis, we develop a protocol for reconstructing the AIMD spectra of liquid water from the sampled distribution of the OH bond lengths. Our findings elucidate the origins of the relationship between the molecular structure of the fleeting hydrogen-bonded network and the ensuing dynamics, which can be probed by vibrational spectroscopy.

    View details for PubMedID 30735052

  • The Quest for Accurate Liquid Water Properties from First Principles JOURNAL OF PHYSICAL CHEMISTRY LETTERS Pestana, L., Marsalek, O., Markland, T. E., Head-Gordon, T. 2018; 9 (17): 5009–16


    Developing accurate ab initio molecular dynamics (AIMD) models that capture both electronic reorganization and nuclear quantum effects associated with hydrogen bonding is key to quantitative understanding of bulk water and its anomalies as well as its role as a universal solvent. For condensed phase simulations, AIMD has typically relied on the generalized gradient approximation (GGA) of density functional theory (DFT) as the underlying model chemistry for the potential energy surface, with nuclear quantum effects (NQEs) sometimes modeled by performing classical molecular dynamics simulations at elevated temperatures. Here we show that the properties of liquid water obtained from the meta-GGA B97M-rV functional, when evaluated using accelerated path integral molecular dynamics simulations, display accuracy comparable to a computationally expensive dispersion-corrected hybrid functional, revPBE0-D3. We show that the meta-GGA DFT functional reproduces bulk water properties including radial distribution functions, self-diffusion coefficients, and infrared spectra with comparable accuracy of a much more expensive hybrid functional. This work demonstrates that the underlying quality of a good DFT functional requires evaluation with quantum nuclei and that high-temperature simulations are a poor proxy for properly treating NQEs.

    View details for PubMedID 30118601

  • On the exact continuous mapping of fermions. Scientific reports Montoya-Castillo, A., Markland, T. E. 2018; 8 (1): 12929


    We derive a rigorous, quantum mechanical map of fermionic creation and annihilation operators to continuous Cartesian variables that exactly reproduces the matrix structure of the many-fermion problem. We show how our scheme can be used to map a general many-fermion Hamiltonian and then consider two specific models that encode the fundamental physics of many fermionic systems, the Anderson impurity and Hubbard models. We use these models to demonstrate how efficient mappings of these Hamiltonians can be constructed using a judicious choice of index ordering of the fermions. This development provides an alternative exact route to calculate the static and dynamical properties of fermionic systems and sets the stage to exploit the quantum-classical and semiclassical hierarchies to systematically derive methods offering a range of accuracies, thus enabling the study of problems where the fermionic degrees of freedom are coupled to complex anharmonic nuclear motion and spins which lie beyond the reach of most currently available methods.

    View details for PubMedID 30154503

  • Unraveling electronic absorption spectra using nuclear quantum effects: Photoactive yellow protein and green fluorescent protein chromophores in water JOURNAL OF CHEMICAL PHYSICS Zuehlsdorff, T. J., Napoli, J. A., Milanese, J. M., Markland, T. E., Isborn, C. M. 2018; 149 (2): 024107


    Many physical phenomena must be accounted for to accurately model solution-phase optical spectral line shapes, from the sampling of chromophore-solvent configurations to the electronic-vibrational transitions leading to vibronic fine structure. Here we thoroughly explore the role of nuclear quantum effects, direct and indirect solvent effects, and vibronic effects in the computation of the optical spectrum of the aqueously solvated anionic chromophores of green fluorescent protein and photoactive yellow protein. By analyzing the chromophore and solvent configurations, the distributions of vertical excitation energies, the absorption spectra computed within the ensemble approach, and the absorption spectra computed within the ensemble plus zero-temperature Franck-Condon approach, we show how solvent, nuclear quantum effects, and vibronic transitions alter the optical absorption spectra. We find that including nuclear quantum effects in the sampling of chromophore-solvent configurations using ab initio path integral molecular dynamics simulations leads to improved spectral shapes through three mechanisms. The three mechanisms that lead to line shape broadening and a better description of the high-energy tail are softening of heavy atom bonds in the chromophore that couple to the optically bright state, widening the distribution of vertical excitation energies from more diverse solvation environments, and redistributing spectral weight from the 0-0 vibronic transition to higher energy vibronic transitions when computing the Franck-Condon spectrum in a frozen solvent pocket. The absorption spectra computed using the combined ensemble plus zero-temperature Franck-Condon approach yield significant improvements in spectral shape and width compared to the spectra computed with the ensemble approach. Using the combined approach with configurations sampled from path integral molecular dynamics trajectories presents a significant step forward in accurately modeling the absorption spectra of aqueously solvated chromophores.

    View details for PubMedID 30007372

  • Decoding the spectroscopic features and time scales of aqueous proton defects. The Journal of chemical physics Napoli, J. A., Marsalek, O., Markland, T. E. 2018; 148 (22): 222833


    Acid solutions exhibit a variety of complex structural and dynamical features arising from the presence of multiple interacting reactive proton defects and counterions. However, disentangling the transient structural motifs of proton defects in the water hydrogen bond network and the mechanisms for their interconversion remains a formidable challenge. Here, we use simulations treating the quantum nature of both the electrons and nuclei to show how the experimentally observed spectroscopic features and relaxation time scales can be elucidated using a physically transparent coordinate that encodes the overall asymmetry of the solvation environment of the proton defect. We demonstrate that this coordinate can be used both to discriminate the extremities of the features observed in the linear vibrational spectrum and to explain the molecular motions that give rise to the interconversion time scales observed in recent nonlinear experiments. This analysis provides a unified condensed-phase picture of the proton structure and dynamics that, at its extrema, encompasses proton sharing and spectroscopic features resembling the limiting Eigen [H3O(H2O)3]+ and Zundel [H(H2O)2]+ gas-phase structures, while also describing the rich variety of interconverting environments in the liquid phase.

    View details for PubMedID 29907063

  • Nuclear quantum effects enter the mainstream NATURE REVIEWS CHEMISTRY Markland, T. E., Ceriotti, M. 2018; 2 (3)
  • The Interplay of Structure and Dynamics in the Raman Spectrum of Liquid Water over the Full Frequency and Temperature Range JOURNAL OF PHYSICAL CHEMISTRY LETTERS Morawietz, T., Marsalek, O., Pattenaude, S. R., Streacker, L. M., Ben-Amotz, D., Markland, T. E. 2018; 9 (4): 851–57


    While many vibrational Raman spectroscopy studies of liquid water have investigated the temperature dependence of the high-frequency O-H stretching region, few have analyzed the changes in the Raman spectrum as a function of temperature over the entire spectral range. Here, we obtain the Raman spectra of water from its melting to boiling point, both experimentally and from simulations using an ab initio-trained machine learning potential. We use these to assign the Raman bands and show that the entire spectrum can be well described as a combination of two temperature-independent spectra. We then assess which spectral regions exhibit strong dependence on the local tetrahedral order in the liquid. Further, this work demonstrates that changes in this structural parameter can be used to elucidate the temperature dependence of the Raman spectrum of liquid water and provides a guide to the Raman features that signal water ordering in more complex aqueous systems.

    View details for PubMedID 29394069

  • Proton Network Flexibility Enables Robustness and Large Electric Fields in the Ketosteroid Isomerase Active Site JOURNAL OF PHYSICAL CHEMISTRY B Wang, L., Fried, S. D., Markland, T. E. 2017; 121 (42): 9807–15


    Hydrogen-bond networks play vital roles in biological functions ranging from protein folding to enzyme catalysis. Here we combine electronic structure calculations and ab initio path integral molecular dynamics simulations, which incorporate both nuclear and electronic quantum effects, to show why the network of short hydrogen bonds in the active site of ketosteroid isomerase is remarkably robust to mutations along the network and how this gives rise to large local electric fields. We demonstrate that these properties arise from the network's ability to respond to a perturbation by shifting proton positions and redistributing electronic charge density. This flexibility leads to small changes in properties such as the partial ionization of residues and pKa isotope effects upon mutation of the residues, consistent with recent experiments. This proton flexibility is further enhanced when an extended hydrogen-bond network forms in the presence of an intermediate analogue, which allows us to explain the chemical origins of the large electric fields in the enzyme's active site observed in recent experiments.

    View details for PubMedID 28915043

  • Unravelling the influence of quantum proton delocalization on electronic charge transfer through the hydrogen bond CHEMICAL PHYSICS LETTERS Schran, C., Marsalek, O., Markland, T. E. 2017; 678: 289–95
  • Quantum Dynamics and Spectroscopy of Ab Initio Liquid Water: The Interplay of Nuclear and Electronic Quantum Effects JOURNAL OF PHYSICAL CHEMISTRY LETTERS Marsalek, O., Markland, T. E. 2017; 8 (7): 1545-1551


    Understanding the reactivity and spectroscopy of aqueous solutions at the atomistic level is crucial for the elucidation and design of chemical processes. However, the simulation of these systems requires addressing the formidable challenges of treating the quantum nature of both the electrons and nuclei. Exploiting our recently developed methods that provide acceleration by up to 2 orders of magnitude, we combine path integral simulations with on-the-fly evaluation of the electronic structure at the hybrid density functional theory level to capture the interplay between nuclear quantum effects and the electronic surface. Here we show that this combination provides accurate structure and dynamics, including the full infrared and Raman spectra of liquid water. This allows us to demonstrate and explain the failings of lower-level density functionals for dynamics and vibrational spectroscopy when the nuclei are treated quantum mechanically. These insights thus provide a foundation for the reliable investigation of spectroscopy and reactivity in aqueous environments.

    View details for DOI 10.1021/acs.jpclett.7b00391

    View details for Web of Science ID 000398884800035

    View details for PubMedID 28296422

  • Electrostatic Control of Regioselectivity in Au(I)-Catalyzed Hydroarylation JOURNAL OF THE AMERICAN CHEMICAL SOCIETY Lau, V. M., Pfalzgraff, W. C., Markland, T. E., Kanan, M. W. 2017; 139 (11): 4035-4041


    Competing pathways in catalytic reactions often involve transition states with very different charge distributions, but this difference is rarely exploited to control selectivity. The proximity of a counterion to a charged catalyst in an ion paired complex gives rise to strong electrostatic interactions that could be used to energetically differentiate transition states. Here we investigate the effects of ion pairing on the regioselectivity of the hydroarylation of 3-substituted phenyl propargyl ethers catalyzed by cationic Au(I) complexes, which forms a mixture of 5- and 7-substituted 2H-chromenes. We show that changing the solvent dielectric to enforce ion pairing to a SbF6(-) counterion changes the regioselectivity by up to a factor of 12 depending on the substrate structure. Density functional theory (DFT) is used to calculate the energy difference between the putative product-determining isomeric transition states (ΔΔE(‡)) in both the presence and absence of the counterion. The change in ΔΔE(‡) upon switching from the unpaired transition states in high solvent dielectric to ion paired transition states in low solvent dielectric (Δ(ΔΔE(‡))) was found to be in good agreement with the experimentally observed selectivity changes across several substrates. Our calculations indicate that the origin of Δ(ΔΔE(‡)) lies in the preferential electrostatic stabilization of the transition state with greater charge separation by the counterion in the ion paired case. By performing calculations at multiple different values of the solvent dielectric, we show that the role of the solvent in changing selectivity is not solely to enforce ion pairing, but rather that interactions between the ion paired complex and the solvent also contribute to Δ(ΔΔE(‡)). Our results provide a foundation for exploiting electrostatic control of selectivity in other ion paired systems.

    View details for DOI 10.1021/jacs.6b11971

    View details for Web of Science ID 000397477700021

    View details for PubMedID 28225605

  • Nuclear Quantum Effects in Water and Aqueous Systems: Experiment, Theory, and Current Challenges CHEMICAL REVIEWS Ceriotti, M., Fang, W., Kusalik, P. G., McKenzie, R. H., Michaelides, A., Morales, M. A., Markland, T. E. 2016; 116 (13): 7529-7550


    Nuclear quantum effects influence the structure and dynamics of hydrogen-bonded systems, such as water, which impacts their observed properties with widely varying magnitudes. This review highlights the recent significant developments in the experiment, theory, and simulation of nuclear quantum effects in water. Novel experimental techniques, such as deep inelastic neutron scattering, now provide a detailed view of the role of nuclear quantum effects in water's properties. These have been combined with theoretical developments such as the introduction of the principle of competing quantum effects that allows the subtle interplay of water's quantum effects and their manifestation in experimental observables to be explained. We discuss how this principle has recently been used to explain the apparent dichotomy in water's isotope effects, which can range from very large to almost nonexistent depending on the property and conditions. We then review the latest major developments in simulation algorithms and theory that have enabled the efficient inclusion of nuclear quantum effects in molecular simulations, permitting their combination with on-the-fly evaluation of the potential energy surface using electronic structure theory. Finally, we identify current challenges and future opportunities in this area of research.

    View details for DOI 10.1021/acs.chemrev.5b00674

    View details for Web of Science ID 000379794000004

    View details for PubMedID 27049513

  • Generalized quantum master equations in and out of equilibrium: When can one win? JOURNAL OF CHEMICAL PHYSICS Kelly, A., Montoya-Castillo, A., Wang, L., Markland, T. E. 2016; 144 (18)


    Generalized quantum master equations (GQMEs) are an important tool in modeling chemical and physical processes. For a large number of problems, it has been shown that exact and approximate quantum dynamics methods can be made dramatically more efficient, and in the latter case more accurate, by proceeding via the GQME formalism. However, there are many situations where utilizing the GQME approach with an approximate method has been observed to return the same dynamics as using that method directly. Here, for systems both in and out of equilibrium, we provide a more detailed understanding of the conditions under which using an approximate method can yield benefits when combined with the GQME formalism. In particular, we demonstrate the necessary manipulations, which are satisfied by exact quantum dynamics, that are required to recast the memory kernel in a form that can be analytically shown to yield the same result as a direct application of the dynamics regardless of the approximation used. By considering the connections between these forms of the kernel, we derive the conditions that approximate methods must satisfy if they are to offer different results when used in conjunction with the GQME formalism. These analytical results thus provide new insights as to when proceeding via the GQME approach can be used to improve the accuracy of simulations.

    View details for DOI 10.1063/1.4948612

    View details for Web of Science ID 000377711900006

    View details for PubMedID 27179469

  • Unraveling the dynamics and structure of functionalized self-assembled monolayers on gold using 2D IR spectroscopy and MD simulations PROCEEDINGS OF THE NATIONAL ACADEMY OF SCIENCES OF THE UNITED STATES OF AMERICA Yan, C., Yuan, R., Pfalzgraff, W. C., Nishida, J., Wang, L., Markland, T. E., Fayer, M. D. 2016; 113 (18): 4929-4934


    Functionalized self-assembled monolayers (SAMs) are the focus of ongoing investigations because they can be chemically tuned to control their structure and dynamics for a wide variety of applications, including electrochemistry, catalysis, and as models of biological interfaces. Here we combine reflection 2D infrared vibrational echo spectroscopy (R-2D IR) and molecular dynamics simulations to determine the relationship between the structures of functionalized alkanethiol SAMs on gold surfaces and their underlying molecular motions on timescales of tens to hundreds of picoseconds. We find that at higher head group density, the monolayers have more disorder in the alkyl chain packing and faster dynamics. The dynamics of alkanethiol SAMs on gold are much slower than the dynamics of alkylsiloxane SAMs on silica. Using the simulations, we assess how the different molecular motions of the alkyl chain monolayers give rise to the dynamics observed in the experiments.

    View details for DOI 10.1073/pnas.1603080113

    View details for Web of Science ID 000375395700026

    View details for PubMedID 27044113

    View details for PubMedCentralID PMC4983838

  • Ab initio molecular dynamics with nuclear quantum effects at classical cost: Ring polymer contraction for density functional theory JOURNAL OF CHEMICAL PHYSICS Marsalek, O., Markland, T. E. 2016; 144 (5)


    Path integral molecular dynamics simulations, combined with an ab initio evaluation of interactions using electronic structure theory, incorporate the quantum mechanical nature of both the electrons and nuclei, which are essential to accurately describe systems containing light nuclei. However, path integral simulations have traditionally required a computational cost around two orders of magnitude greater than treating the nuclei classically, making them prohibitively costly for most applications. Here we show that the cost of path integral simulations can be dramatically reduced by extending our ring polymer contraction approach to ab initio molecular dynamics simulations. By using density functional tight binding as a reference system, we show that our ring polymer contraction scheme gives rapid and systematic convergence to the full path integral density functional theory result. We demonstrate the efficiency of this approach in ab initio simulations of liquid water and the reactive protonated and deprotonated water dimer systems. We find that the vast majority of the nuclear quantum effects are accurately captured using contraction to just the ring polymer centroid, which requires the same number of density functional theory calculations as a classical simulation. Combined with a multiple time step scheme using the same reference system, which allows the time step to be increased, this approach is as fast as a typical classical ab initio molecular dynamics simulation and 35× faster than a full path integral calculation, while still exactly including the quantum sampling of nuclei. This development thus offers a route to routinely include nuclear quantum effects in ab initio molecular dynamics simulations at negligible computational cost.

    View details for DOI 10.1063/1.4941093

    View details for Web of Science ID 000369893900014

    View details for PubMedID 26851913

  • Simulating Nuclear and Electronic Quantum Effects in Enzymes. Methods in enzymology Wang, L., Isborn, C. M., Markland, T. E. 2016; 577: 389-418


    An accurate treatment of the structures and dynamics that lead to enhanced chemical reactivity in enzymes requires explicit treatment of both electronic and nuclear quantum effects. The former can be captured in ab initio molecular dynamics (AIMD) simulations, while the latter can be included by performing ab initio path integral molecular dynamics (AI-PIMD) simulations. Both AIMD and AI-PIMD simulations have traditionally been computationally prohibitive for large enzymatic systems. Recent developments in streaming computer architectures and new algorithms to accelerate path integral simulations now make these simulations practical for biological systems, allowing elucidation of enzymatic reactions in unprecedented detail. In this chapter, we summarize these recent developments and discuss practical considerations for applying AIMD and AI-PIMD simulations to enzymes.

    View details for DOI 10.1016/bs.mie.2016.05.047

    View details for PubMedID 27498646

  • Nonadiabatic Dynamics in Atomistic Environments: Harnessing Quantum-Classical Theory with Generalized Quantum Master Equations JOURNAL OF PHYSICAL CHEMISTRY LETTERS Pfalzgraff, W. C., Kelly, A., Markland, T. E. 2015; 6 (23): 4743-4748
  • Nonadiabatic Dynamics in Atomistic Environments: Harnessing Quantum-Classical Theory with Generalized Quantum Master Equations. The journal of physical chemistry letters Pfalzgraff, W. C., Kelly, A., Markland, T. E. 2015; 6 (23): 4743-8


    The development of methods that can efficiently and accurately treat nonadiabatic dynamics in quantum systems coupled to arbitrary atomistic environments remains a significant challenge in problems ranging from exciton transport in photovoltaic materials to electron and proton transfer in catalysis. Here we show that our recently introduced MF-GQME approach, which combines Ehrenfest mean field theory with the generalized quantum master equation framework, is able to yield quantitative accuracy over a wide range of charge-transfer regimes in fully atomistic environments. This is accompanied by computational speed-ups of up to 3 orders of magnitude over a direct application of Ehrenfest theory. This development offers the opportunity to efficiently investigate the atomistic details of nonadiabatic quantum relaxation processes in regimes where obtaining accurate results has previously been elusive.

    View details for DOI 10.1021/acs.jpclett.5b02131

    View details for PubMedID 26563917

  • Accurate nonadiabatic quantum dynamics on the cheap: making the most of mean field theory with master equations. journal of chemical physics Kelly, A., Brackbill, N., Markland, T. E. 2015; 142 (9): 094110-?


    In this article, we show how Ehrenfest mean field theory can be made both a more accurate and efficient method to treat nonadiabatic quantum dynamics by combining it with the generalized quantum master equation framework. The resulting mean field generalized quantum master equation (MF-GQME) approach is a non-perturbative and non-Markovian theory to treat open quantum systems without any restrictions on the form of the Hamiltonian that it can be applied to. By studying relaxation dynamics in a wide range of dynamical regimes, typical of charge and energy transfer, we show that MF-GQME provides a much higher accuracy than a direct application of mean field theory. In addition, these increases in accuracy are accompanied by computational speed-ups of between one and two orders of magnitude that become larger as the system becomes more nonadiabatic. This combination of quantum-classical theory and master equation techniques thus makes it possible to obtain the accuracy of much more computationally expensive approaches at a cost lower than even mean field dynamics, providing the ability to treat the quantum dynamics of atomistic condensed phase systems for long times.

    View details for DOI 10.1063/1.4913686

    View details for PubMedID 25747064

  • Accurate nonadiabatic quantum dynamics on the cheap: Making the most of mean field theory with master equations. journal of chemical physics Kelly, A., Brackbill, N., Markland, T. E. 2015; 142 (9): 094110-?

    View details for DOI 10.1063/1.4913686

    View details for PubMedID 25747064

  • Quantum delocalization of protons in the hydrogen-bond network of an enzyme active site. Proceedings of the National Academy of Sciences of the United States of America Wang, L., Fried, S. D., Boxer, S. G., Markland, T. E. 2014; 111 (52): 18454-18459


    Enzymes use protein architectures to create highly specialized structural motifs that can greatly enhance the rates of complex chemical transformations. Here, we use experiments, combined with ab initio simulations that exactly include nuclear quantum effects, to show that a triad of strongly hydrogen-bonded tyrosine residues within the active site of the enzyme ketosteroid isomerase (KSI) facilitates quantum proton delocalization. This delocalization dramatically stabilizes the deprotonation of an active-site tyrosine residue, resulting in a very large isotope effect on its acidity. When an intermediate analog is docked, it is incorporated into the hydrogen-bond network, giving rise to extended quantum proton delocalization in the active site. These results shed light on the role of nuclear quantum effects in the hydrogen-bond network that stabilizes the reactive intermediate of KSI, and the behavior of protons in biological systems containing strong hydrogen bonds.

    View details for DOI 10.1073/pnas.1417923111

    View details for PubMedID 25503367

    View details for PubMedCentralID PMC4284547

  • Quantum fluctuations and isotope effects in ab initio descriptions of water JOURNAL OF CHEMICAL PHYSICS Wang, L., Ceriotti, M., Markland, T. E. 2014; 141 (10)


    Isotope substitution is extensively used to investigate the microscopic behavior of hydrogen bonded systems such as liquid water. The changes in structure and stability of these systems upon isotope substitution arise entirely from the quantum mechanical nature of the nuclei. Here, we provide a fully ab initio determination of the isotope exchange free energy and fractionation ratio of hydrogen and deuterium in water treating exactly nuclear quantum effects and explicitly modeling the quantum nature of the electrons. This allows us to assess how quantum effects in water manifest as isotope effects, and unravel how the interplay between electronic exchange and correlation and nuclear quantum fluctuations determine the structure of the hydrogen bond in water.

    View details for DOI 10.1063/1.4894287

    View details for Web of Science ID 000342209400040

  • Quantum fluctuations and isotope effects in ab initio descriptions of water. journal of chemical physics Wang, L., Ceriotti, M., Markland, T. E. 2014; 141 (10): 104502-?


    Isotope substitution is extensively used to investigate the microscopic behavior of hydrogen bonded systems such as liquid water. The changes in structure and stability of these systems upon isotope substitution arise entirely from the quantum mechanical nature of the nuclei. Here, we provide a fully ab initio determination of the isotope exchange free energy and fractionation ratio of hydrogen and deuterium in water treating exactly nuclear quantum effects and explicitly modeling the quantum nature of the electrons. This allows us to assess how quantum effects in water manifest as isotope effects, and unravel how the interplay between electronic exchange and correlation and nuclear quantum fluctuations determine the structure of the hydrogen bond in water.

    View details for DOI 10.1063/1.4894287

    View details for PubMedID 25217932

  • Multiple time step integrators in ab initio molecular dynamics. journal of chemical physics Luehr, N., Markland, T. E., Martínez, T. J. 2014; 140 (8): 084116-?


    Multiple time-scale algorithms exploit the natural separation of time-scales in chemical systems to greatly accelerate the efficiency of molecular dynamics simulations. Although the utility of these methods in systems where the interactions are described by empirical potentials is now well established, their application to ab initio molecular dynamics calculations has been limited by difficulties associated with splitting the ab initio potential into fast and slowly varying components. Here we present two schemes that enable efficient time-scale separation in ab initio calculations: one based on fragment decomposition and the other on range separation of the Coulomb operator in the electronic Hamiltonian. We demonstrate for both water clusters and a solvated hydroxide ion that multiple time-scale molecular dynamics allows for outer time steps of 2.5 fs, which are as large as those obtained when such schemes are applied to empirical potentials, while still allowing for bonds to be broken and reformed throughout the dynamics. This permits computational speedups of up to 4.4x, compared to standard Born-Oppenheimer ab initio molecular dynamics with a 0.5 fs time step, while maintaining the same energy conservation and accuracy.

    View details for DOI 10.1063/1.4866176

    View details for PubMedID 24588157

  • Interface-Limited Growth of Heterogeneously Nucleated Ice in Supercooled Water JOURNAL OF PHYSICAL CHEMISTRY B Nistor, R. A., Markland, T. E., Berne, B. J. 2014; 118 (3): 752-760


    Heterogeneous ice growth exhibits a maximum in freezing rate arising from the competition between kinetics and the thermodynamic driving force between the solid and liquid states. Here, we use molecular dynamics simulations to elucidate the atomistic details of this competition, focusing on water properties in the interfacial region along the secondary prismatic direction. The crystal growth velocity is maximized when the efficiency of converting interfacial water molecules to ice, collectively known as the attachment kinetics, is greatest. We find water molecules that contact the intermediate ice layer in concave regions along the atomistically roughened surface are more likely to freeze directly. An increased roughening of the solid surface at large undercoolings consequently plays an important limiting role in the rate of ice growth, as water molecules are unable to integrate into increasingly deeper surface pockets. These results provide insight into the molecular mechanisms for self-assembly of solid phases that are important in many biological and atmospheric processes.

    View details for DOI 10.1021/jp408832b

    View details for Web of Science ID 000330252700012

    View details for PubMedID 24393086

  • Efficient and accurate surface hopping for long time nonadiabatic quantum dynamics JOURNAL OF CHEMICAL PHYSICS Kelly, A., Markland, T. E. 2013; 139 (1)


    The quantum-classical Liouville equation offers a rigorous approach to nonadiabatic quantum dynamics based on surface hopping type trajectories. However, in practice the applicability of this approach has been limited to short times owing to unfavorable numerical scaling. In this paper we show that this problem can be alleviated by combining it with a formally exact generalized quantum master equation treatment. This allows dramatic improvements in the efficiency of the approach in nonadiabatic regimes, making it computationally tractable to treat the quantum dynamics of complex systems for long times. We demonstrate our approach by applying it to a model of condensed phase charge transfer where our method is shown to be numerically exact in regimes where fewest-switches surface hopping and mean field approaches fail to obtain either the correct rates or long-time populations.

    View details for DOI 10.1063/1.4812355

    View details for Web of Science ID 000321716400006

    View details for PubMedID 23822290

  • Efficient methods and practical guidelines for simulating isotope effects JOURNAL OF CHEMICAL PHYSICS Ceriotti, M., Markland, T. E. 2013; 138 (1)


    The shift in chemical equilibria due to isotope substitution is frequently exploited to obtain insight into a wide variety of chemical and physical processes. It is a purely quantum mechanical effect, which can be computed exactly using simulations based on the path integral formalism. Here we discuss how these techniques can be made dramatically more efficient, and how they ultimately outperform quasi-harmonic approximations to treat quantum liquids not only in terms of accuracy, but also in terms of computational cost. To achieve this goal we introduce path integral quantum mechanics estimators based on free energy perturbation, which enable the evaluation of isotope effects using only a single path integral molecular dynamics trajectory of the naturally abundant isotope. We use as an example the calculation of the free energy change associated with H/D and (16)O/(18)O substitutions in liquid water, and of the fractionation of those isotopes between the liquid and the vapor phase. In doing so, we demonstrate and discuss quantitatively the relative benefits of each approach, thereby providing a set of guidelines that should facilitate the choice of the most appropriate method in different, commonly encountered scenarios. The efficiency of the estimators we introduce and the analysis that we perform should in particular facilitate accurate ab initio calculation of isotope effects in condensed phase systems.

    View details for DOI 10.1063/1.4772676

    View details for Web of Science ID 000313330000013

    View details for PubMedID 23298033

  • Ring-Polymer Molecular Dynamics: Quantum Effects in Chemical Dynamics from Classical Trajectories in an Extended Phase Space ANNUAL REVIEW OF PHYSICAL CHEMISTRY, VOL 64 Habershon, S., Manolopoulos, D. E., Markland, T. E., Miller, T. F. 2013; 64: 387-413


    This article reviews the ring-polymer molecular dynamics model for condensed-phase quantum dynamics. This model, which involves classical evolution in an extended ring-polymer phase space, provides a practical approach to approximating the effects of quantum fluctuations on the dynamics of condensed-phase systems. The review covers the theory, implementation, applications, and limitations of the approximation.

    View details for DOI 10.1146/annurev-physchem-040412-110122

    View details for Web of Science ID 000321771600018

    View details for PubMedID 23298242

  • Isotope effects in water as investigated by neutron diffraction and path integral molecular dynamics JOURNAL OF PHYSICS-CONDENSED MATTER Zeidler, A., Salmon, P. S., Fischer, H. E., Neuefeind, J. C., Simonson, J. M., Markland, T. E. 2012; 24 (28)


    The structures of heavy and light water at 300 K were investigated by using a joint approach in which the method of neutron diffraction with oxygen isotope substitution was complemented by path integral molecular dynamics simulations. The diffraction results, which give intra-molecular O-D and O-H bond distances of 0.985(5) and 0.990(5) Å, were found to be in best agreement with those obtained by using the flexible anharmonic TTM3-F water model. Both techniques show a difference of  ≃ 0.5% between the O-D and O-H intra-molecular bond lengths, and the results support a competing quantum effects model for water in which its structural and dynamical properties are governed by an offset between intra-molecular and inter-molecular quantum contributions. Further consideration of the O-O correlations is needed in order to improve agreement with experiment.

    View details for DOI 10.1088/0953-8984/24/28/284126

    View details for Web of Science ID 000305786400028

    View details for PubMedID 22738936

  • Zeidler et al. Reply PHYSICAL REVIEW LETTERS Zeidler, A., Salmon, P. S., Fischer, H. E., Neuefeind, J. C., Simonson, J. M., Lemmel, H., Rauch, H., Markland, T. E. 2012; 108 (25)
  • Growing Point-to-Set Length Scale Correlates with Growing Relaxation Times in Model Supercooled Liquids PHYSICAL REVIEW LETTERS Hocky, G. M., Markland, T. E., Reichman, D. R. 2012; 108 (22)


    It has been demonstrated recently that supercooled liquids sharing simple structural features (e.g. pair distribution functions) may exhibit strikingly distinct dynamical behavior. Here we show that a more subtle structural feature correlates with relaxation times in three simulated systems that have nearly identical radial distribution functions but starkly different dynamical behavior. In particular, for the first time we determine the thermodynamic "point-to-set" length scale in several canonical model systems and demonstrate the quantitative connection between this length scale and the growth of relaxation times. Our results provide clues necessary for distinguishing competing theories of the glass transition.

    View details for DOI 10.1103/PhysRevLett.108.225506

    View details for Web of Science ID 000304695900010

    View details for PubMedID 23003622

  • Unraveling quantum mechanical effects in water using isotopic fractionation PROCEEDINGS OF THE NATIONAL ACADEMY OF SCIENCES OF THE UNITED STATES OF AMERICA Markland, T. E., Berne, B. J. 2012; 109 (21): 7988-7991


    When two phases of water are at equilibrium, the ratio of hydrogen isotopes in each is slightly altered because of their different phase affinities. This isotopic fractionation process can be utilized to analyze water's movement in the world's climate. Here we show that equilibrium fractionation ratios, an entirely quantum mechanical property, also provide a sensitive probe to assess the magnitude of nuclear quantum fluctuations in water. By comparing the predictions of a series of water models, we show that those describing the OH chemical bond as rigid or harmonic greatly overpredict the magnitude of isotope fractionation. Models that account for anharmonicity in this coordinate are shown to provide much more accurate results because of their ability to give partial cancellation between inter- and intramolecular quantum effects. These results give evidence of the existence of competing quantum effects in water and allow us to identify how this cancellation varies across a wide-range of temperatures. In addition, this work demonstrates that simulation can provide accurate predictions and insights into hydrogen fractionation.

    View details for DOI 10.1073/pnas.1203365109

    View details for Web of Science ID 000304445800021

    View details for PubMedID 22566650

    View details for PubMedCentralID PMC3361433

  • Reduced density matrix hybrid approach: Application to electronic energy transfer JOURNAL OF CHEMICAL PHYSICS Berkelbach, T. C., Markland, T. E., Reichman, D. R. 2012; 136 (8)


    Electronic energy transfer in the condensed phase, such as that occurring in photosynthetic complexes, frequently occurs in regimes where the energy scales of the system and environment are similar. This situation provides a challenge to theoretical investigation since most approaches are accurate only when a certain energetic parameter is small compared to others in the problem. Here we show that in these difficult regimes, the Ehrenfest approach provides a good starting point for a dynamical description of the energy transfer process due to its ability to accurately treat coupling to slow environmental modes. To further improve on the accuracy of the Ehrenfest approach, we use our reduced density matrix hybrid framework to treat the faster environmental modes quantum mechanically, at the level of a perturbative master equation. This combined approach is shown to provide an efficient and quantitative description of electronic energy transfer in a model dimer and the Fenna-Matthews-Olson complex and is used to investigate the effect of environmental preparation on the resulting dynamics.

    View details for DOI 10.1063/1.3687342

    View details for Web of Science ID 000300944000006

    View details for PubMedID 22380029

  • Theory and simulations of quantum glass forming liquids JOURNAL OF CHEMICAL PHYSICS Markland, T. E., Morrone, J. A., Miyazaki, K., Berne, B. J., Reichman, D. R., Rabani, E. 2012; 136 (7)


    A comprehensive microscopic dynamical theory is presented for the description of quantum fluids as they transform into glasses. The theory is based on a quantum extension of mode-coupling theory. Novel effects are predicted, such as reentrant behavior of dynamical relaxation times. These predictions are supported by path integral ring polymer molecular dynamics simulations. The simulations provide detailed insight into the factors that govern slow dynamics in glassy quantum fluids. Connection to other recent work on both quantum glasses as well as quantum optimization problems is presented.

    View details for DOI 10.1063/1.3684881

    View details for Web of Science ID 000300551000025

    View details for PubMedID 22360252

  • Reduced density matrix hybrid approach: An efficient and accurate method for adiabatic and non-adiabatic quantum dynamics JOURNAL OF CHEMICAL PHYSICS Berkelbach, T. C., Reichman, D. R., Markland, T. E. 2012; 136 (3)


    We present a new approach to calculate real-time quantum dynamics in complex systems. The formalism is based on the partitioning of a system's environment into "core" and "reservoir" modes with the former to be treated quantum mechanically and the latter classically. The presented method only requires the calculation of the system's reduced density matrix averaged over the quantum core degrees of freedom which is then coupled to a classically evolved reservoir to treat the remaining modes. We demonstrate our approach by applying it to the spin-boson problem using the noninteracting blip approximation to treat the system and core, and Ehrenfest dynamics to treat the reservoir. The resulting hybrid methodology is accurate for both fast and slow baths, since it naturally reduces to its composite methods in their respective regimes of validity. In addition, our combined method is shown to yield good results in intermediate regimes where neither approximation alone is accurate and to perform equally well for both strong and weak system-bath coupling. Our approach therefore provides an accurate and efficient methodology for calculating quantum dynamics in complex systems.

    View details for DOI 10.1063/1.3671372

    View details for Web of Science ID 000299387700015

    View details for PubMedID 22280750

  • Oxygen as a Site Specific Probe of the Structure of Water and Oxide Materials PHYSICAL REVIEW LETTERS Zeidler, A., Salmon, P. S., Fischer, H. E., Neuefeind, J. C., Simonson, J. M., Lemmel, H., Rauch, H., Markland, T. E. 2011; 107 (14)


    The method of oxygen isotope substitution in neutron diffraction is introduced as a site specific structural probe. It is employed to measure the structure of light versus heavy water, thus circumventing the assumption of isomorphism between H and D as used in more traditional neutron diffraction methods. The intramolecular and intermolecular O-H and O-D pair correlations are in excellent agreement with path integral molecular dynamics simulations, both techniques showing a difference of ≃0.5% between the O-H and O-D intramolecular bond distances. The results support the validity of a competing quantum effects model for water in which its structural and dynamical properties are governed by an offset between intramolecular and intermolecular quantum contributions.

    View details for DOI 10.1103/PhysRevLett.107.145501

    View details for Web of Science ID 000296285800014

    View details for PubMedID 22107211

  • Quantum fluctuations can promote or inhibit glass formation NATURE PHYSICS Markland, T. E., Morrone, J. A., Berne, B. J., Miyazaki, K., Rabani, E., Reichman, D. R. 2011; 7 (2): 134-137

    View details for DOI 10.1038/NPHYS1865

    View details for Web of Science ID 000286807000015

  • Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances JOURNAL OF CHEMICAL PHYSICS Morrone, J. A., Markland, T. E., Ceriotti, M., Berne, B. J. 2011; 134 (1)


    Multiple time scale molecular dynamics enhances computational efficiency by updating slow motions less frequently than fast motions. However, in practice, the largest outer time step possible is limited not by the physical forces but by resonances between the fast and slow modes. In this paper we show that this problem can be alleviated by using a simple colored noise thermostatting scheme which selectively targets the high frequency modes in the system. For two sample problems, flexible water and solvated alanine dipeptide, we demonstrate that this allows the use of large outer time steps while still obtaining accurate sampling and minimizing the perturbation of the dynamics. Furthermore, this approach is shown to be comparable to constraining fast motions, thus providing an alternative to molecular dynamics with constraints.

    View details for DOI 10.1063/1.3518369

    View details for Web of Science ID 000286010600006

    View details for PubMedID 21218993

  • Efficient stochastic thermostatting of path integral molecular dynamics JOURNAL OF CHEMICAL PHYSICS Ceriotti, M., Parrinello, M., Markland, T. E., Manolopoulos, D. E. 2010; 133 (12)


    The path integral molecular dynamics (PIMD) method provides a convenient way to compute the quantum mechanical structural and thermodynamic properties of condensed phase systems at the expense of introducing an additional set of high frequency normal modes on top of the physical vibrations of the system. Efficiently sampling such a wide range of frequencies provides a considerable thermostatting challenge. Here we introduce a simple stochastic path integral Langevin equation (PILE) thermostat which exploits an analytic knowledge of the free path integral normal mode frequencies. We also apply a recently developed colored noise thermostat based on a generalized Langevin equation (GLE), which automatically achieves a similar, frequency-optimized sampling. The sampling efficiencies of these thermostats are compared with that of the more conventional Nosé-Hoover chain (NHC) thermostat for a number of physically relevant properties of the liquid water and hydrogen-in-palladium systems. In nearly every case, the new PILE thermostat is found to perform just as well as the NHC thermostat while allowing for a computationally more efficient implementation. The GLE thermostat also proves to be very robust delivering a near-optimum sampling efficiency in all of the cases considered. We suspect that these simple stochastic thermostats will therefore find useful application in many future PIMD simulations.

    View details for DOI 10.1063/1.3489925

    View details for Web of Science ID 000282648000007

    View details for PubMedID 20886921

  • A fast path integral method for polarizable force fields JOURNAL OF CHEMICAL PHYSICS Fanourgakis, G. S., Markland, T. E., Manolopoulos, D. E. 2009; 131 (9)


    A quantum simulation of an imaginary time path integral typically requires around n times more computational effort than the corresponding classical simulation, where n is the number of ring polymer beads (or imaginary time slices) used in the calculation. It is however possible to improve on this estimate by decomposing the potential into a sum of slowly and rapidly varying contributions. If the slowly varying contribution changes only slightly over the length scale of the ring polymer, it can be evaluated on a contracted ring polymer with fewer than the full n beads (or equivalently on a lower order Fourier decomposition of the imaginary time path). Here we develop and test this idea for systems with polarizable force fields. The development consists of iterating the induction on the contracted ring polymer and applying an appropriate transformation to obtain the forces on the original n beads. In combination with a splitting of the Coulomb potential into its short- and long-range parts, this results in a method with little more than classical computational effort in the limit of large system size. The method is illustrated with simulations of liquid water at 300 K and hexagonal ice at 100 K using a recently developed flexible and polarizable Thole-type potential energy model.

    View details for DOI 10.1063/1.3216520

    View details for Web of Science ID 000269625400003

    View details for PubMedID 19739844

  • Competing quantum effects in the dynamics of a flexible water model JOURNAL OF CHEMICAL PHYSICS Habershon, S., Markland, T. E., Manolopoulos, D. E. 2009; 131 (2)


    Numerous studies have identified large quantum mechanical effects in the dynamics of liquid water. In this paper, we suggest that these effects may have been overestimated due to the use of rigid water models and flexible models in which the intramolecular interactions were described using simple harmonic functions. To demonstrate this, we introduce a new simple point charge model for liquid water, q-TIP4P/F, in which the O-H stretches are described by Morse-type functions. We have parametrized this model to give the correct liquid structure, diffusion coefficient, and infrared absorption frequencies in quantum (path integral-based) simulations. The model also reproduces the experimental temperature variation of the liquid density and affords reasonable agreement with the experimental melting temperature of hexagonal ice at atmospheric pressure. By comparing classical and quantum simulations of the liquid, we find that quantum mechanical fluctuations increase the rates of translational diffusion and orientational relaxation in our model by a factor of around 1.15. This effect is much smaller than that observed in all previous simulations of empirical water models, which have found a quantum effect of at least 1.4 regardless of the quantum simulation method or the water model employed. The small quantum effect in our model is a result of two competing phenomena. Intermolecular zero point energy and tunneling effects destabilize the hydrogen-bonding network, leading to a less viscous liquid with a larger diffusion coefficient. However, this is offset by intramolecular zero point motion, which changes the average water monomer geometry resulting in a larger dipole moment, stronger intermolecular interactions, and a slower diffusion. We end by suggesting, on the basis of simulations of other potential energy models, that the small quantum effect we find in the diffusion coefficient is associated with the ability of our model to produce a single broad O-H stretching band in the infrared absorption spectrum.

    View details for DOI 10.1063/1.3167790

    View details for Web of Science ID 000267983100037

    View details for PubMedID 19603998

  • A refined ring polymer contraction scheme for systems with electrostatic interactions CHEMICAL PHYSICS LETTERS Markland, T. E., Manolopoulos, D. E. 2008; 464 (4-6): 256-261
  • An efficient ring polymer contraction scheme for imaginary time path integral simulations JOURNAL OF CHEMICAL PHYSICS Markland, T. E., Manolopoulos, D. E. 2008; 129 (2)


    A quantum simulation of an imaginary time path integral typically requires around n times more computational effort than the corresponding classical simulation, where n is the number of ring polymer beads (or imaginary time slices) used in the calculation. However, this estimate neglects the fact that the potential energies of many systems can be decomposed into a sum of rapidly varying short-range and slowly varying long-range contributions. For such systems, the computational effort of the path integral simulation can be reduced considerably by evaluating the long-range forces on a contracted ring polymer with fewer beads than are needed to evaluate the short-range forces. This idea is developed and then illustrated with an application to a flexible model of liquid water in which the intramolecular forces are evaluated with 32 beads, the oxygen-oxygen Lennard-Jones forces with seven, and the intermolecular electrostatic forces with just five. The resulting static and dynamic properties are within a few percent of those of a full 32-bead calculation, and yet they are obtained with a computational effort less than six times (rather than 32 times) that of a classical simulation. We hope that this development will encourage future studies of quantum mechanical fluctuations in liquid water and aqueous solutions and in many other systems with similar interaction potentials.

    View details for DOI 10.1063/1.2953308

    View details for Web of Science ID 000257629100006

    View details for PubMedID 18624514

  • Quantum diffusion of hydrogen and muonium atoms in liquid water and hexagonal ice JOURNAL OF CHEMICAL PHYSICS Markland, T. E., Habershon, S., Manolopoulos, D. E. 2008; 128 (19)


    We have used the ring polymer molecular dynamics method to study the diffusion of muonium, hydrogen, and deuterium atoms in liquid water and hexagonal ice over a wide temperature range (8-361 K). Quantum effects are found to dramatically reduce the diffusion of muonium in water relative to that predicted by classical simulation. This leads to a simple explanation for the lack of any significant isotope effect in the observed diffusion coefficients of these species in the room temperature liquid. Our results indicate that the mechanism of the diffusion in liquid water is similar to the intercavity hopping mechanism observed in ice, supplemented by the diffusion of the cavities in the liquid. Within the same model, we have also been able to simulate the observed crossover in the c-axis diffusion coefficients of hydrogen and deuterium in hexagonal ice. Finally, we have been able to obtain good agreement with experimental data on the diffusion of muonium in hexagonal ice at 8 K, where the process is entirely quantum mechanical.

    View details for DOI 10.1063/1.2925792

    View details for Web of Science ID 000256205200034

    View details for PubMedID 18500879