Nonadiabatic molecular dynamics simulations aim to describe the coupled electron- nuclear dynamics of molecules in excited electronic states, beyond the celebrated Born-Oppenheimer approximation. These simulations have been applied to understand a plethora of photochemical and photophysical processes and to support the interpretation of ultrafast spectroscopy experiments at advanced light sources. As a result, the number of nonadiabatic dynamics simulations has been growing significantly over the past decade. Yet, the field remains in its infancy, and a potential user may find it difficult to approach this type of simulation, given their complexity and the number of elements that should be considered for a (hopefully) successful nonadiabatic dynamics simulation. Nonadiabatic molecular dynamics relies on several key steps: finding a level of electronic-structure theory to describe the molecule in its Franck-Condon region and beyond, describing the photoexcitation process, selecting a method to perform the nonadiabatic dynamics, and analyzing the final results before calculating observables for a more direct comparison with experiment. This Best Practices guide aims to provide a general guide for the user of nonadiabatic molecular dynamics by (i) discussing the fundamentals of nonadiabatic molecular dynamics and the various trajectory-based methods developed for molecular systems, (ii) introducing the different electronic-structure methods and concepts – adiabatic/diabatic representation, conical intersections – that can be used with nonadiabatic molecular dynamics (or for benchmarking), (iii) providing details on the various steps required to perform a nonadiabatic dynamics simulation and their practical use, as well as guided examples and a discussion on the calculation of observables, (iv) proposing a FAQ with the typical questions a user may have when performing nonadiabatic dynamics, and (v) sketching a checklist for the key practical steps when performing a (trajectory-based) nonadiabatic molecular dynamics. Each section is self-contained, but we endeavor to provide additional key references for each concept discussed, making this Guide a starting point for the interested reader to dig further into the field of nonadiabatic dynamics.
Parameterizing modified nucleic acids is a difficult but necessary task for expanding the simulated space of oligonucleotides, including both naturally occurring structures and those with pharmaceutical relevance. In lieu of expensive and difficult chemical synthesis in the laboratory, computer simulations are often performed to make predictions for sequence and structure effects, as well as downstream critical quality attributes. To enable these simulations, modifications have to be parameterized to faithfully represent their effect on nucleotides. This is a non-trivial process, complicated by the fact that it may be the first thing researchers must figure out before they can build their structures and start their initial simulations. To enable these research projects, we created modXNA, a code that assembles pre-parameterized modules of the base, backbone, and sugar, to create bespoke combinations of modifications. In the following tutorial, we provide background on force field parameterization in the Amber software ecosystem and detail the steps necessary to perform parameterization of modified nucleic acids using modXNA.
This tutorial describes the practical use of some recent methodological advances implemented in the GROMOS software for biomolecular simulations. It is envisioned as a living document, with additional tutorials being added in the course of time. Currently, it consists of six distinct tutorials. The first tutorial describes the use of time-averaged restraints to enforce agreement with order parameters derived from NMR experiments. The second tutorial describes the use of extended thermodynamic integration in the double-decoupling method to compute the affinity of a small molecule to a protein. The molecule involved bears a negative charge, necessitating the application of post-simulation corrections. The third tutorial is based on the same molecular system, but computes the binding free energy from a path-sampling method with distance-field distance restraints and Hamiltonian replica exchange simulations. The fourth tutorial describes the use of Gaussian accelerated MD, an enhanced sampling technique. The fifth tutorial is about the use of Accelerated Enveloping Distribution Sampling (AEDS) an enhanced version of EDS. The sixth tutorial illustrates how to run QM/MM simulations with the Buffer Region Neural Network (BuRNN) method. The tutorials are written for users with some experience in the application of molecular dynamics simulations.
This tutorial equips natural scientists with the essential knowledge needed to utilize and comprehend molecular dynamics simulations effectively. Beginning with the stability of integration algorithms and the conservation of energy, this article proceeds with the description of atomistic forcefields, thermodynamic ensembles, long-range electrostatics treatment, and free energy calculations within molecular dynamics simulations. It then extends to the simulation of proteins and drug discovery applications. This comprehensive overview includes numerous references to relevant publications and tackles real-world problems. The tutorial is based on a 10-week master’s course typically undertaken by students in the first or second semester of their master’s program.
The availability of open-source molecular simulation software packages allows scientists and engineers to focus on running and analyzing simulations without having to write, parallelize, and validate their own simulation software. While molecular simulations thus become accessible to a larger audience, the ‘‘black box’’ nature of such software packages and wide array of options and features can make it challenging to use them correctly, particularly for beginners in the topic of simulations. LAMMPS is one such versatile molecular simulation code, designed for modeling particle-based systems across a broad range of materials science and computational chemistry applications, including atomistic, coarse-grained, mesoscale, grid-free continuum, and discrete element models. LAMMPS is capable of efficiently running simulations of varying sizes from small desktop computers to large-scale supercomputing environments. Its flexibility and extensibility make it ideal for complex and extensive simulations of atomic and molecular systems, and beyond. This article introduces a suite of tutorials designed to make learning LAMMPS more accessible to new users. The first four tutorials cover the basics of running molecular simulations in LAMMPS with systems of varying complexities. The second four tutorials address more advanced molecular simulation techniques, specifically the use of a reactive force field, grand canonical Monte Carlo, enhanced sampling, and the REACTER protocol. In addition, we introduce LAMMPS–GUI, an enhanced cross-platform graphical text editor specifically designed for use with LAMMPS and able to run LAMMPS directly on the edited input. LAMMPS–GUI is used as the primary tool in the tutorials to edit inputs, run LAMMPS, extract data, and visualize the simulated systems.
This review article provides an overview of structurally oriented experimental datasets that can be used to benchmark protein force fields, focusing on data generated by nuclear magnetic resonance (NMR) spectroscopy and room temperature (RT) protein crystallography. We discuss what the observables are, what they tell us about structure and dynamics, what makes them useful for assessing force field accuracy, and how they can be connected to molecular dynamics simulations carried out using the force field one wishes to benchmark. We also touch on best practices for setup and analysis of benchmark simulations. We hope this article will be particularly useful to computational researchers and trainees who develop, benchmark, or use protein force fields or machine learning models that generate protein ensembles.
Grid Inhomogeneous Solvation Theory (GIST) is a method to compute the free energy of solvation of a solute molecule on a three-dimensional grid based on sampling from molecular dynamics (MD) simulations. The high spatial resolution of the GIST output, as well as the decomposition into energy and entropy contributions, allow for highly detailed analyses of solvation around both proteins and small molecules. However, this versatility also comes with a significant entry barrier for new users. In this tutorial, we aim to guide the reader through the most common steps involved in a GIST analysis using the streptavidin-biotin complex as a demonstrative system. To this end, Jupyter notebooks and a Python package (gisttools) are provided to simplify the analysis. Furthermore, we discuss the theory of GIST with a focus on practical aspects. We highlight potential pitfalls and provide strategies to avoid technical difficulties. This tutorial assumes familiarity with molecular dynamics simulations and the AmberTools package.
Recent advances in machine learning (ML) are reshaping drug discovery. Structure-based ML methods use physically-inspired models to predict binding affinities from protein:ligand complexes. These methods promise to enable the integration of data for many related targets, which addresses issues related to data scarcity for single targets and could enable generalizable predictions for a broad range of targets, including mutations. In this work, we report our experiences in building KinoML, a novel framework for ML in target-based small molecule drug discovery with an emphasis on structure-enabled methods. KinoML focuses currently on kinases as the relative structural conservation of this protein superfamily, particularly in the kinase domain, means it is possible to leverage data from the entire superfamily to make structure-informed predictions about binding affinities, selectivities, and drug resistance. Some key lessons learned in building KinoML include the importance of reproducible data collection and deposition, the harmonization of molecular data and featurization, and the selection of the right data format to ensure reusability and reproducibility of ML models. As a result, KinoML allows users to easily achieve three tasks: accessing and curating molecular data; featurizing this data with representations suitable for ML applications; and running reproducible ML experiments that require access to ligand, protein, and assay information to predict ligand affinity. Despite KinoML focusing on kinases, this framework can be applied to other proteins. The lessons reported here can help guide the development of platforms for structure-enabled ML in other areas of drug discovery.
Gaussian-accelerated molecular dynamics (GaMD) simulations are an advanced technique that enhances the sampling of configurational space by applying biasing potentials that reduce energy barriers, enabling faster exploration of the free energy landscape. This tutorial demonstrates the application of GaMD to the alanine dipeptide, serving as an accessible model system, and guides users through all GaMD simulation stages: conventional MD, GaMD equilibration, GaMD production, and reweighting. Users will gain practical insights into the preparation of input files, monitoring of GaMD convergence, and analysis of free energy profiles using PyReweighting. We make a particular effort to connect the underlying theory with the GaMD workflow. This tutorial is intended for users with prior molecular dynamics experience, Linux and command-line navigation, and with basic Python knowledge. The step-by-step instructions and accompanying scripts aim to streamline the GaMD workflow, making it accessible for the broader research community to explore enhanced sampling for a range of biomolecular systems.
Although Monte Carlo (MC) is a very powerful molecular simulation method in statistical mechanics, the development and application of novel MC trials to optimize sampling in complex systems is hindered by the difficulty in deriving their acceptance probabilities. We present a checklist approach to deriving acceptance probabilities, and apply this approach to a variety of trials in the canonical, isothermal-isobaric, grand-canonical, semi-grand canonical and Gibbs ensembles. The ideal gas is then shown to be a useful test case to compare the results of simulations with those from theoretical expectations, providing a computational benchmark that can be easily and rapidly implemented for determining if the acceptance criteria were derived correctly. More complex models and trials are also considered with this checklist approach, including configurational bias, cavity bias, energy bias, aggregation volume bias, dual-cut configurational bias and rigid cluster moves. The result is a framework designed to help researchers implement and test specialized MC trials that expand the model complexity and length scales currently available in open-source MC molecular simulation software. Sample code is also provided in the GitHub repository [https://github.com/usnistgov/best-practices-mc].
This review article provides an overview of structurally oriented experimental datasets that can be used to benchmark protein force fields, focusing on data generated by nuclear magnetic resonance (NMR) spectroscopy and room temperature (RT) protein crystallography. We discuss what the observables are, what they tell us about structure and dynamics, what makes them useful for assessing force field accuracy, and how they can be connected to molecular dynamics simulations carried out using the force field one wishes to benchmark. We also touch on statistical issues that arise when comparing simulations with experiment. We hope this article will be particularly useful to computational researchers and trainees who develop, benchmark, or use protein force fields for molecular simulations.
SEEKR2 (Simulation enabled estimation of kinetic rates v. 2) is a powerful and versatile software tool designed to computationally estimate the kinetics and thermodynamics of complex molecular processes, particularly emphasizing the process of receptor-ligand binding and unbinding. We present a suite of tutorials for the SEEKR2 (Simulation enabled estimation of kinetic rates v. 2) multiscale milestoning software. This tutorial presents a comprehensive guide for users offering the best practices for preparing, executing, and analyzing molecular dynamics (MD) and Brownian dynamics (BD) simulations using SEEKR2. This tutorial highlights the advancements presented in SEEKR2 - the latest iteration within the SEEKR programs, including significant improvements in speed and capabilities compared to its earlier versions. SEEKR2 now supports both NAMD and OpenMM simulation engines, providing users with more flexibility in their simulation setups. Additionally, the BD component has been upgraded to the Browndye2 engine, enhancing the accuracy and efficiency of simulations. This tutorial aims to guide users to install SEEKR2, run MD and BD simulations within the framework of the SEEKR2 program, and analyze and interpret the kinetics and thermodynamics of binding and unbinding of model host-guest systems, thereby demonstrating its ease of usability and extensible features that allow for future expansions of the method. This tutorial equips users with the necessary knowledge to effectively prepare, execute, and analyze simulations using SEEKR2. By following the best practices outlined in the tutorial, users can leverage the power of the SEEKR2 program to gain insights into complex molecular processes and accelerate their understanding of key biomolecular interactions.
The weighted ensemble (WE) strategy has been demonstrated to be highly efficient in generating pathways and rate constants for rare events such as protein folding and protein binding using atomistic molecular dynamics simulations. Here we present two sets of tutorials instructing users in the best practices for preparing, carrying out, and analyzing WE simulations for various applications using the WESTPA software. The first set of more basic tutorials describes a range of simulation types, from a molecular association process in explicit solvent to more complex processes such as host-guest association, peptide conformational sampling, and protein folding. The second set ecompasses six advanced tutorials instructing users in the best practices of using key new features and plugins/extensions of the WESTPA 2.0 software package, which consists of major upgrades for larger systems and/or slower processes. The advanced tutorials demonstrate the use of the following key features: (i) a generalized resampler module for the creation of "binless" schemes, (ii) a minimal adaptive binning scheme for more efficient surmounting of free energy barriers, (iii) streamlined handling of large simulation datasets using an HDF5 framework, (iv) two different schemes for more efficient rate-constant estimation, (v) a Python API for simplified analysis of WE simulations, and (vi) plugins/extensions for Markovian Weighted Ensemble Milestoning and WE rule-based modeling for systems biology models. Applications of the advanced tutorials include atomistic and non-spatial models, and consist of complex processes such as protein folding and the membrane permeability of a drug-like molecule. Users are expected to already have significant experience with running conventional molecular dynamics or systems biology simulations.
Free Energy Perturbation (FEP) is a powerful but challenging computational technique for estimating differences in free energy between two or more states. This document is intended both as a tutorial and as an adaptable protocol for computing free energies of binding using free energy perturbations in NAMD. We present the Streamlined Alchemical Free Energy Perturbation (SAFEP) framework. SAFEP shifts the computational frame of reference from the ligand to the binding site itself. This both simplifies the thermodynamic cycle and makes the approach more broadly applicable to superficial sites and other less common geometries. As a practical example, we give instructions for calculating the absolute binding free energy of phenol to lysozyme. We assume familiarity with standard procedures for setting up, running, and analyzing molecular dynamics simulations using NAMD and VMD. While simulation times will vary, the human tasks should take no more than 3 to 4 hours for a reader without previous training in free energy calculations or experience with the VMD Colvars Dashboard. Sample data are provided for all key calculations both for comparison and readers’ convenience.
The weighted ensemble (WE) strategy has been demonstrated to be highly efficient in generating pathways and rate constants for rare events such as protein folding and protein binding using atomistic molecular dynamics simulations. Here we present two sets of tutorials instructing users in the best practices for preparing, carrying out, and analyzing WE simulations for various applications using the WESTPA software. The first set of more basic tutorials describes a range of simulation types, from a molecular association process in explicit solvent to more complex processes such as host-guest association, peptide conformational sampling, and protein folding. The second set ecompasses six advanced tutorials instructing users in the best practices of using key new features and plugins/extensions of the WESTPA 2.0 software package, which consists of major upgrades for larger systems and/or slower processes. The advanced tutorials demonstrate the use of the following key features: (i) a generalized resampler module for the creation of "binless" schemes, (ii) a minimal adaptive binning scheme for more efficient surmounting of free energy barriers, (iii) streamlined handling of large simulation datasets using an HDF5 framework, (iv) two different schemes for more efficient rate-constant estimation, (v) a Python API for simplified analysis of WE simulations, and (vi) plugins/extensions for Markovian Weighted Ensemble Milestoning and WE rule-based modeling for systems biology models. Applications of the advanced tutorials include atomistic and non-spatial models, and consist of complex processes such as protein folding and the membrane permeability of a drug-like molecule. Users are expected to already have significant experience with running conventional molecular dynamics or systems biology simulations.
Pysimm is a framework for molecular simulations of polymers and polymer-based nanostructures, which enables their direct chemical synthesis and preparation. Pysimm facilitates the understanding of novel, amorphous, processable materials for a broad range of applications, including heterogeneous catalysts, adsorbents and gas storage materials, as well as protein-polymer conjugates. This tutorial provides a detailed guide on the construction of atomistic and united-atom models of polymers using Pysimm: an open-source Python Application Programming Interface for atomistic molecular simulations. The API complements and simplifies the work of widely known molecular simulation software, such as LAMMPS, CASSANDRA, NAMD and Amber. Readers should be familiar with the basic concepts of atomistic molecular simulations, as well as the basic knowledge of Python programming language, before attempting to follow this tutorial. This work is separated into 3 main sections. First, the process of building an atomic-level model of a polymer chain from its repetitive units is described. The second section shows how to work with existing forcefields, and how Pysimm can automatically read, recognize, and assign appropriate Force Field parameters to a molecule. The final section discusses how to use Pysimm to construct polymer chains with pre-specified tacticity. The section is also available in the form of an interactive Jupyter notebook tutorial outlining simple guidelines to construct polymer models.
Hydrogen-deuterium exchange (HDX) is a comprehensive yet detailed probe of protein structure and dynamics and, coupled to mass spectrometry, has become a powerful tool for investigating an increasingly large array of systems. Computer simulations are often used to help rationalize experimental observations of exchange, but interpretations have frequently been limited to simple, subjective correlations between microscopic dynamical fluctuations and the observed macroscopic exchange behavior. With this in mind, we previously developed the HDX ensemble reweighting approach and associated software, HDXer, to aid the objective interpretation of HDX data using molecular simulations. HDXer has two main functions; first, to compute H-D exchange rates that describe each structure in a candidate ensemble of protein structures, for example from molecular simulations, and second, to objectively reweight the conformational populations present in a candidate ensemble to conform to experimental exchange data. In this article, we first describe the HDXer approach, theory, and implementation. We then guide users through a suite of tutorials that demonstrate the practical aspects of preparing experimental data, computing HDX levels from molecular simulations, and performing ensemble reweighting analyses. Finally we provide a practical discussion of the capabilities and limitations of the HDXer methods including recommendations for a user's own analyses. Overall, this article is intended to provide an up-to-date, pedagogical counterpart to the software, which is freely available at https://github.com/Lucy-Forrest-Lab/HDXer.
