<figure><img src="image_23.png" /><figcaption>Figure 14. Performance of QM/MM methods.¹⁴⁷ (A) Wall time in hours per 1-ns MD simulation for insulin receptor kinase versus the number of CPU cores, for the MM-only, SQUANTM, MNDO97 BOMD, and MNDO97 with DXL-BOMD methods. (B) Energy conservation in different SCF accelerator implementations for adenylate kinase, showing less than 1 kcal/mol deviation of the total energy during 100-ps NVE MD simulations. Simulation systems consist of (A) 76 QM atoms and 28,823 MM atoms and (B) 92 QM atoms and 47,201 MM atoms. In both cases, the QM region is treated with the AM1/d-PhoT SE QM model⁶⁰⁸ for the MNDO97 module and with the AM1 model for the SQUANTM module.</figcaption></figure>

with all QM atoms and thus directly polarizes the QM electron density.

**11.3. Recent Advances in SE-QM/MM Methods.** Main strengths of SE-QM/MM methods over AI/DFT-QM/MM methods are their efficiency and flexible functional forms that allow recalibration against target data.⁶⁰⁸,⁶³³,⁶³⁴ Extended MD simulations within reasonable computational time is thus possible while yielding accuracy tailored to individual systems. In enzyme mechanism studies, SE-QM/MM methods are frequently employed together with the RXNCOR module for US simulations and more recently, with the SM.¹³⁹,⁶⁴⁶,⁶⁴⁷ They can also be used in FEP simulations such as calculating solvation free energies of solutes and ligand pKₐ values, functionalities available in the QUANTUM, SQUANTM, and SCC-DFTB modules.⁵⁹⁸,⁶⁴⁵,⁶⁴⁸,⁶⁴⁹

CHARMM also has several acceleration algorithms for the SE-QM/MM methods, including the direct inversion of the iterative subspace (DIIS) extrapolation scheme for faster SCF convergence⁶⁵⁰,⁶⁵¹ and the pseudodiagonalization algorithm for the Fock matrix. In addition, the MNDO97 module has recently incorporated the MPI parallelization and several new SCF accelerators, achieving more than 10-fold speed up (Figure 14).¹⁴⁷,¹⁶⁵ The newly implemented SCF accelerators are as follows:

**11.3.1. Extended Lagrangian MD (ELMD).**¹⁶⁵ This method⁶⁵² performs MD simulations with the electron density of the QM subsystem propagated by the Lagrangian:

$$ \mathcal{L}_{\text{ELMD}} = \frac{1}{2} \sum_{A} M_{A} \dot{\mathbf{R}}_{A}^{2} + \frac{1}{2} \sum_{i,j} m_{i,j} \dot{P}_{ij}^{2} - E_{\text{tot}}(\mathbf{R}, \mathbf{P}) \\ - \text{Tr}[\mathbf{\Lambda}(\mathbf{P}^{2} - \mathbf{P})] \qquad (54) $$

Here, **R** denotes the QM and MM coordinate with mass *M*<sub>A</sub> (A: atom index). **P** is the electron density matrix with mass *m*<sub>ij</sub> for the corresponding element. **R**<sub>A</sub> and **P**<sub>ij</sub> are time derivatives. The Lagrange multiplier **Λ** enforces the idempotency constraint on **P**. Alternatively, **P** can be propagated using the curvy-steps unitary update algorithm.⁶⁵³ Both methods avoid time-consuming SCF iteration at the expense of a smaller integration time step for the density matrix propagation. This limitation can be alleviated by applying the multiple time step (MTS) approach where nuclear coordinates are propagated with a larger time step, typically 0.5 or 1 fs, as commonly used in QM/MM MD simulations.

**11.3.2. Extended Lagrangian Born–Oppenheimer MD with Dissipation (DXL-BOMD).**¹⁴⁷ The idempotency condition of Eq. 54 is modified into a harmonic restraint for the auxiliary density variable **D**, to oscillate around the converged SCF density:

$$ \mathcal{L}_{\text{DXL-BOMD}} = \frac{1}{2} \sum_{A} M_{A} \dot{\mathbf{R}}_{A}^{2} + \frac{m_{D}}{2} \sum_{i,j} \dot{\mathbf{D}}_{ij}^{2} - E_{\text{tot}}(\mathbf{R}, \mathbf{P}) \\ - \frac{\kappa_{D}}{2} (\mathbf{D} - \mathbf{P})^{2} \qquad (55) $$

The true (SCF) density **P** is approximated by **D** that serves as the initial guess in the SCF iteration, and κ<sub>D</sub> is the force constant. **D** is extrapolated based on a predetermined number of previous SCF densities.⁶⁵⁴,⁶⁵⁵ The SCF iteration is then performed for a given number of SCF steps. In practice, Eq. 55 is solved in the limit of vanishing *m*<sub>D</sub>, resulting in coupled equations of motion, one for the nuclear position and the other for **D**, thereby eliminating the dependence of results on *m*<sub>D</sub>. This method is also supported by the DFTB QM/MM module.⁶⁵⁵

**11.3.3. Fock Matrix Dynamics (FMD).**¹⁴⁷ The Fock matrix is directly extrapolated based on its elements determined from previous MD steps, followed by regular SCF iteration until convergence.⁶⁵⁶,⁶⁵⁷ Both DXL-BOMD and FMD methods significantly reduce the number of SCF iterations compared to conventional (BOMD) SCF calculations while maintaining energy conservation.¹⁴⁷

In addition to efficiency, SE-QM/MM methods in CHARMM are being developed to improve accuracy. Recognizing the importance of non-bonded interactions, the MNDO97 and DFTB modules have implemented Grimme's dispersion and hydrogen bond correction terms.⁶³⁸,⁶⁵⁸,⁶⁵⁹ To further improve the quality of the PES, a simple valence bond-like (SVB) term has been introduced in the NDDO-based SE-QM/MM modules.⁶⁴⁶,⁶⁶⁰ In addition, the SQUANTM module integrates the SE-QM/MM and GAMESS-UK AI/DFT-QM/MM methods, introducing a dual-level approach that interpolates the QM/MM energy to the AI/DFT-QM/MM level of theory.⁶⁴² This method is compatible with the MTS algorithm so that MD simulations are performed at the SE-QM/MM level while simultaneously correcting energies and gradients at the target AI/DFT-QM/MM level over a longer time step (Figure 15). These developments enable simulations of highly challenging systems with unprecedented accuracy and efficiency, pushing the boundaries of QM/MM methods.
<figure><img src="image_24.png" /><figcaption>Figure 15. PMF for the S$_{N}$2 reaction between CH$_{3}$Cl and Cl$^{-}$ in water.⁶⁴² Results from MTS simulations with varying number of MD steps N, for evaluating the AI-QM/MM correction term. The MTS AI-QM/MM simulations were carried out using the AM1 and HF/3-21G methods for the low- and high-level QM theories, respectively. “AM1 only” results are from the AM1 SE-QM/MM simulations. The inset compares the impact of low-level QM theory on the PMF, while the high-level theory remains at the HF/3-21G level.</figcaption></figure>

## 11.4. Path-Integral-Free Energy Perturbation for Nuclear Quantum Effects.
An important approach for treating nuclear quantum effects (NQE) is the Feynman path integral (PI) formalism⁶⁶¹ that describes the wave function by an ensemble of paths weighed by the classical action for each, to capture Schrödinger's delocalized wave behavior. It is readily generalizable to multiparticle systems, and it naturally accounts for thermal effects and sampling can be performed in MC and MD simulations.

The quantum transition state theory (QTST) rate constant may be computed with PI (PI-QTST).⁶⁶² Consider a system composed of a set of QM atoms with coordinates **r** embedded in a bath of classical atoms with coordinates **R** in thermal equilibrium. The QM partition function $Q_Q$ can be written as the trace of the thermal density matrix $\rho$

$$ 
\begin{aligned}
Q_Q &= \mathrm{Tr}(\rho) = \int d\mathbf{R} \int d\mathbf{r} \rho(\mathbf{r}, \mathbf{r}; \mathbf{R}, \beta) \\
&= \int d\mathbf{R} \int \Pi_{j=1}^p d\mathbf{r} \rho(\mathbf{r}_j, \mathbf{r}_{j+1}; \mathbf{R}, \tau)
\end{aligned}
\quad (56) $$

where P is the number of quasi-particles or beads and $\tau = \beta/P$ ($\beta = 1/k_B T$). For a closed ring-polymer chain, $\mathbf{r}_1 = \mathbf{r}_{p+1}$. In the high-T limit ($\tau \to 0$ and $P \to \infty$), the semiclassical primitive approximation⁶⁶³ can be used for $\rho$

$$ 
\begin{aligned}
\rho_{\mathrm{PA}}(\mathbf{r}_j, \mathbf{r}_{j+1}; \mathbf{R}, \tau) \\
&= \left( \frac{m}{2\pi\tau\hbar^2} \right)^{D/2} \\
&\quad \times e^{-(m/2\tau\hbar^2)(\mathbf{r}_{j+1}-\mathbf{r}_j)^2 + (\tau/2)[V(\mathbf{r}_{j+1}\mathbf{R}) + V(\mathbf{r}_j\mathbf{R})]}
\end{aligned}
\quad (57) $$

where V is the system potential (i.e., a QM/MM potential), m is the quantum particle mass, and D is the dimension of **r**. The description above is isomorphic to a classical ring of beads system connected via harmonic springs, and forms the basis for the implementation in CHARMM where the bead distribution is sampled using MC simulations.

The rate constant of PI-QTST⁶⁶² is defined as

$$ k_Q^{\mathrm{PI-QTST}} = \Omega_{\mathrm{FP}} Q_Q^{\frac{\pm}{2}} / Q_Q^{\mathrm{RS}} \quad (58) $$

where $\Omega_{\mathrm{FP}}$ is a free-particle (FP) prefactor, and $Q_Q^{\pm}$ and $Q_Q^{\mathrm{RS}}$ are the quantum partition functions for the transition and reactant states, respectively. Although one can compute the PI-QTST directly, when using expensive potentials like QM/MM, it is convenient to compute the correction to the classical TST due to NQE

$$ \delta^{\pm} = \frac{Q_Q^{\pm} / Q_C^{\pm}}{Q_Q^{\mathrm{RS}} / Q_C^{\mathrm{RS}}} \quad (59) $$

where the subscript C denotes the corresponding classical partition functions. The above includes both quantum vibration (zero-point energy) and tunneling effects. To calculate the quantum to classical ratio of the partition function, a double average scheme⁶⁶⁴,⁶⁶⁵ can be used

$$ \delta = \frac{Q_Q}{Q_C} = \langle \langle e^{-\tau \sum_{j=1}^p [V(\mathbf{r}_j \mathbf{R}) - V(\mathbf{r}_j \mathbf{R})]} \rangle_{\mathrm{FP}, \mathbf{r}_c} \rangle_{V(\mathbf{r}_j \mathbf{R})} \quad (60) $$

where $\mathbf{r}_c$ is the centroid coordinate. The outer (classical) average is obtained using standard simulation techniques, while the delocalized QM description comes from the inner FP average. This approach is practical since the classical and quantum simulations are performed separately. In CHARMM, the NQE atoms are defined as QM atoms. Currently, the PI method works with all SE-QM/MM modules in CHARMM as well as with AI/DFT-QM/MM using CHARMM and Q-Chem.

A well-known challenge with PI simulations is the difficulty with sampling the polymer ring due to the harmonic coupling between the beads (Eq. 57).⁶⁶³ In CHARMM, the inner average in Eq. 60 for the FP PI sampling is performed using the bisection algorithm⁶⁶⁶,⁶⁶⁷ extended to a ring of quasi-particles⁶⁶⁸,⁶⁶⁹ or using the staging algorithm.⁶⁷⁰,⁶⁷¹ Since each new configuration in the bisection or staging PI sampling is independent of the distribution of beads in previous configurations, the PI rapidly converges, which is essential for accurate calculation of the NQE and absolute rate constants.⁶³³,⁶⁷²⁻⁶⁷⁴ However, specialized techniques are required to precisely compute isotope effects due to minute differences in free energy. In CHARMM, a novel mass-perturbation technique termed PI-FEP, was developed to directly compute the free energy difference between isotopes.⁶⁷⁵ Since relative free energies between the distributions of heavy and light particles are determined by FEP (Figure 16), the precision of the computed kinetic isotope effect (KIE) is of experimental quality for both the primary and secondary KIEs.⁶⁷¹,⁶⁷⁶⁻⁶⁸³

To reduce the computational cost of PI simulations in CHARMM, higher-order factorization of the density matrix operator has been adopted,⁶⁸⁴⁻⁶⁸⁶ which converges with a considerably smaller number of beads at the expense of computing the potential gradient in addition to the potential.⁶⁷¹ This method can be combined with the PI-FEP approach to efficiently compute the KIE.⁶⁷⁹

Information about tunneling can be obtained by inspecting the particle momentum distribution computed using open-chain PI (i.e., $\mathbf{r}_1 \neq \mathbf{r}_{p+1}$).⁶⁸⁷⁻⁶⁸⁹ Whereas closed-chain PI only samples diagonal elements of $\rho$, open-chain PI simulations also sample off-diagonal elements, and it can sample both the anisotropic and isotropic momentum distribution of a transferring hydrogen ($H^+/H^-/H^\cdot$) during a reaction.
<figure><img src="image_25.png" /><figcaption>**Figure 16.** Hydride transfer reaction in dihydrofolate reductase.<sup>633</sup> (A) Reaction mechanism. (B) Classic (solid line) and quantum (dashed lines) PMFs for hydride and deuteride transfer. (C) Active site for the hydride transfer from NADPH to H<sub>2</sub> folate. The transferring hydride is described using PI with 32 beads.</figcaption></figure>

**11.5. Density Functional Tight Binding (DFTB) Module.** One versatile QM/MM module in CHARMM is based on the DFTB approach,<sup>692,693</sup> which is an approximation method that aims to strike the balance between computational efficiency and accuracy. Also referred to as SCC-DFTB,<sup>69</sup> it is an SE QM method in the sense that it employs a minimal basis and the two-center approximation for electron integrals, which are key approximations for NDDO-based SE methods<sup>694</sup> such as AM1, PM3 and OM2 (Section 11.1). On the other hand, most of the parameters in the DFTB approach are computed based on atomic or diatomic molecules, and the most empirical aspect of the parametrization concerns those used to derive the atomic/diatomic electronic properties and the pairwise repulsive potentials (e.g., confinement radius). The most popular approach for biomolecular applications is the DFTB3/3OB model,<sup>64,695</sup> which has been parametrized for elements commonly encountered in organic and biomolecular systems: O, N, C, H, S, P, Na, K, Mg, Ca, Zn, Cu and the halogens. For recent reviews of the development and application of the DFTB3 method for condensed phase applications, see refs 658 and 696.

The DFTB3 model is integrated with MM model in CHARMM through the standard electrostatic embedding scheme<sup>163</sup> where the DFTB3 atoms are represented as Mulliken charges; alternative DFTB3/MM electrostatic interaction models have also been implemented<sup>697</sup> that consider the finite spatial distributions of the DFTB3 and MM charges. In terms of boundary conditions, the DFTB3/MM model can be used together with either a GSBP<sup>644</sup> or the PBC with either Ewald summation<sup>649</sup> or PME<sup>641</sup> (Section 11.2). For localized reactions, the DFTB3/MM-GSBP approach is computationally most efficient and generally agrees with the more expensive DFTB3/MM-PME approach.<sup>141</sup> For systems in which the chemical reaction is coupled with considerable conformational rearrangements, the PBC-based approach is more appropriate. Another technical detail relevant to many QM/MM applications is the flexible inner region ensemble separator (FIRES) potential<sup>698</sup> available in CHARMM (Section 11.7), which prevents the exchange of QM and MM water molecules and thus particularly important for solution reactions<sup>699</sup> or for solvent-accessible active sites.<sup>700</sup>

The DFTB3/MM approach can be used together with many key functionalities in CHARMM, especially various types of

free energy simulations that are essential to quantitative analysis of chemical transformations. For chemical reactions, they include US with the RXNCOR module, an interface with PLUMED<sup>701</sup> for various metadynamics simulations,<sup>682,690,702–704</sup> and the SM available in the STRINGM module<sup>139,297,705</sup> (Section 7.3). Another useful approach for improved sampling is replica exchange US through the REPDSTR module.<sup>297</sup> For alchemical free energy simulations, DFTB3/MM works with the PERT module.<sup>297,706</sup> For applications such as redox potential<sup>310</sup> and pK<sub>a</sub> calculations,<sup>649,707</sup> the DFTB3/MM model also works with the BLOCK module.

In the following, the DFTB3/MM method is illustrated with two types of free energy simulations. For chemical reaction, the catalysis in Usb1, an exoribonuclease that shortens the oligouridine tail of U6 snRNA,<sup>690</sup> is used as an example. In the proposed catalytic mechanism, two active-site histidine residues serve as the catalytic base and acid, respectively (Figure 17A). It was studied with DFTB3/MM metadynamics

<figure><img src="image_26.png" /><figcaption>**Figure 17.** DFTB3/MM free energy simulations for the catalysis in Usb1.<sup>690</sup> (A) The putative catalytic mechanism of Usb1. (B) The 3D free energy surface from DFTB3/MM multiwalker metadynamics simulations. The three CVs describe proton transfers associated with the catalytic acid/base and the phosphoryl transfer (red arrows in panel A). (C) The transition state structure from the DFTB3/MM free energy simulations features the transfer of a single proton between H208 and the leaving group. This was subsequently confirmed with proton inventory experiments.<sup>690</sup> Panels A and C were adapted from ref 690, which was published by the Oxford University Press.</figcaption></figure>

simulations using three CVs that describe proton transfer involving the catalytic base (H120), the phosphoryl transfer reaction, and the proton transfer involving the catalytic acid (H208), respectively. Multiwalker metadynamics calculations were run with 300–500 walkers, each of which was sampled for 0.5–1 ns, leading to a cumulative sampling of 0.1–0.2 µs for constructing the 3D PMF (Figure 17B). This level of sampling for QM/MM simulations is possible only with SE-type QM methods such as DFTB3, highlighting the value of calibrated SE QM/MM simulations for complex biomolecular
processes.<sup>696,703</sup> The transition state structure captured in the DFTB3/MM simulations (Figure 17C) suggests that one proton is in flight, and the predicted feature was subsequently confirmed experimentally.<sup>690</sup>

For alchemical free energy simulations, the binding selectivity of Mg<sup>2+</sup> and Ca<sup>2+</sup> in the Ca<sup>2+</sup> binding protein carp parvalbumin (CP) and its mutant (D51A/E101D/F102W)<sup>691</sup> is considered (Figure 18A). Experimentally,<sup>708</sup> the WT CP was

<figure><img src="image_27.png" /><figcaption>**Figure 18.** DFTB3/MM FEP simulations for the Ca<sup>2+</sup>/Mg<sup>2+</sup> binding selectivity in the Ca<sup>2+</sup> binding protein, carp parvalbumin (CP), and its mutant (D51A/E101D/F102W). (A) The DFTB3/MM-GSBP setup that illustrates the QM regions for the WT and mutant CP simulations. (B) Thermodynamic cycle used to probe the free energy change for the Ca<sup>2+</sup>/Mg<sup>2+</sup> conversion in a given environment. In simulations, the horizontal transitions occur at the MM level, and in the vertical transitions, a metal binding site is converted between MM and QM treatments using the PERT module in CHARMM. Reproduced from ref 691. Copyright [2024] American Chemical Society.</figcaption></figure>

measured to bind more strongly to Ca²⁺ by 5.6 kcal/mol; in the mutant, Ca²⁺ is still preferred over Mg²⁺, although the selectivity is reduced to 1.6 kcal/mol. With a standard FF, the relative binding free energy of Mg²⁺ and Ca²⁺ to a protein is readily computed using alchemical free energy simulations; the metal ions are interconverted twice: once in the binding pocket and once in solution, and their free energy difference is the relative binding affinity. With non-polarizable FFs, such calculations did not yield the correct trend and the smaller Mg²⁺ was predicted to bind more strongly.<sup>691,709</sup> While the same set of alchemical free energy simulations can, in principle, be carried out at the DFTB3/MM level, we adopt an alternative thermodynamic cycle, which involves converting the description of the metal ion and its ligands between MM and DFTB3 levels<sup>297</sup> (Figure 18B). This has the advantage that structural changes during the MM/DFTB3 conversion are

expected to be small and therefore convergence of the free energy simulation is rapid, minimizing the required DFTB3/MM computations. Encouragingly, with DFTB3/MM simulations,<sup>691</sup> the calculated ΔΔG<sub>bind</sub> was ~6.2 kcal/mol, in good agreement with the experimental value. These results highlight the value of a QM description of the metal binding site and support the role of electronic polarization<sup>709</sup> and charge transfer<sup>710</sup> in metal binding to proteins. For the mutant, different binding site models led to considerable variations in the computed relative binding affinities. With a coordination number of seven for Ca²⁺, which was shown by DFTB3/MM metadynamics to be the dominant coordination number for the mutant, the calculated relative binding affinity was ~4.2 kcal/mol, also in fair agreement with the experimental value.

**11.6. Double Link Atom Method (DLAM).** Accurately modeling chemical reactions in condensed phases using QM/MM is challenging, especially when partitioning across a covalent bond. Introducing dummy or link atoms serves as a bridge across the severed bond, acting as both a connection between QM and MM interfaces and an electron density cap for partitions. Typically hydrogen, a link atom can resemble the electronic character or features lost during truncation and is attached to a host group (HG) via the host atom (HA). The link atom is subject to interface with the MM and QM HGs. In the standard single link atom (SLA) scheme,<sup>601,623</sup> the link atom is added to the QM HA to compensate or neutralize the charge of the QM fragment and cap the QM Hamiltonian. However, the addition of the QM link atom introduces a number of challenges at the QM/MM interface, primarily the treatment of electrostatics.

Various approaches have been developed to address QM/MM boundary effects in the QM fragment by adjusting the magnitude of polarization from the MM fragment, including the excluded group (EXGR) scheme,<sup>623</sup> the charge shift scheme (CHSH),<sup>711,712</sup> the divided frontier charge (DIV) scheme,<sup>624</sup> and the distributed Gaussian (DG) method.<sup>713,714</sup> For example, the SLA scheme treats MM atoms as point charges and excludes the MM HA from the QM/MM electrostatics, leaving an artificial partial charge at the interface on the MM HG. The EXGR scheme corrects for the added unrealistic partial charge at the interface by excluding all partial charges of the MM HG, whereas the DIV scheme corrects for the unrealistic charge at the interface by redistributing the excluded MM HA partial charge to the MM HG, and the CHSH scheme introduces a dipole to counterbalance for the charge shift. Alternatively, the DG method includes all electrostatic interactions of the MM HG; however, the MM partial charges are represented as Gaussian charge distributions and utilize a smoothing potential or blur width (σ) to smear the MM electron density where optimal σ values vary depending on the physical property of interest.<sup>626,714</sup> Although these SLA-based schemes provide a balance for the QM fragment and interfacial electrostatics, they neglect the MM fragment which results in unbalanced forces and unrealistic electrostatics.

DLAM<sup>626</sup> addresses the above issues by the addition of a link atom to cap the MM fragment but it has been infrequently used due to its challenging implementation. It can now be called directly in CHARMM to add link atoms to both the QM and MM host fragments. The DLAMadd command adds a QM link atom “QQ” to the QM HA, similar to the single link atom command (ADDLink), and adds an MM link atom “QM” to the MM HA. By default, the link atoms are placed 1.0
Å colinearly from the respective HA. The MM link, typically an MM hydrogen, bears a small partial charge that should preserve both the net charge and dipole of the MM HG. Balance can be achieved by shifting charge between the MM link and MM HA, but can be less straightforward in some systems. DLAM is used with the DG method and employs the same σ value for all MM atoms. As σ and MM link partial charge are free parameters in DLAM that can be tuned to balance electrostatics, identifying reliable parameters for complex systems can be a challenge. Currently, optimizations of σ and MM link partial charges for amino acids compatible with the CHARMM and Amber FFs are underway for ease of use.

## 11.7. Flexible Inner Region Ensemble Separator (FIRES).
A QM/MM methodology in which a solute and the nearest water molecules are represented at a high *ab initio* level, offers a powerful strategy to study the hydration structure around small ions in the aqueous phase. However, one challenge with solvent molecules in hybrid QM/MM simulation is that they are free to diffuse away from the region of interest, and be replaced by MM solvent molecules that provide presumably a less accurate model. To resolve this issue, FIRES was designed in which the ion and a fixed number of nearest water molecules form a dynamical and flexible inner region that is represented with a high level *ab initio* QM method, while the water molecules in the surrounding bulk form an outer region that is represented by a classical FF. Simulations with FIRES yield rigorously correct thermodynamic averages as long as the solvent molecules in the flexible inner and outer regions are not allowed to exchange. The method was used to study hydration structure around Na⁺ and K⁺,⁶⁹⁸ and Mg²⁺ and Zn²⁺.⁷¹⁵ To obtain a more efficient dynamical propagation algorithm, it is necessary to manage the computational cost of the QM part. To this end, a MTS dual-Hamiltonian propagation algorithm was designed by which the trajectory is propagated at every time step via a computationally inexpensive QM Hamiltonian, and then corrected less frequently using a more accurate and computationally expensive QM Hamiltonian.⁷¹⁶

## 11.8. Combining QM/MM with Gaussian Process Machine Learning Potentials.
The pyCHARMM⁹ interface in CHARMM has facilitated advanced uses of QM/MM potentials in conjunction with Python-based ML potentials, including those described by neural networks⁵⁴ and Gaussian process regression (GPR).⁷¹⁷,⁷¹⁸ Built upon multivariate Gaussian distribution of latent functions, GPR is a non-parametric, kernel-based stochastic inference ML approach that maximizes the likelihood of training data observation.⁷¹⁹ In simulations, GPR has been employed to model the relationship between molecular descriptors and the PES (reviewed in ref 363). Recently, GPR has been utilized to develop delta-ML potentials to improve SE-QM/MM free energy simulations.⁷¹⁷,⁷¹⁸ By combining the AM1/MM potential in CHARMM and energy-based streaming sparse GPR (SSGPR) models, AI-QM/MM quality PES information can be learned along the string free energy paths in a data-efficient manner.⁷¹⁸ Using the extended-kernel GPR with derivative observations (GPRwDO), both energy and force matching can be employed to improve SE-QM/MM free energy simulations.⁷¹⁷ Figure 19A shows the PMFs for the Menshutkin reaction simulated at the AM1/MM level, before and after deploying the GPR correction model where the latter significantly alleviates overestimation of the free energy barrier

<figure><img src="image_28.png" /><figcaption>Figure 19. (A) PMFs of the Menshutkin reaction (NH₃ + CH₃Cl → NH₃CH₃⁺ + Cl⁻) simulated at the AM1/MM and AM1-GPR/MM levels.⁷¹⁷ (B) The scheme of combining GPR Python libraries and CHARMM through pyCHARMM.</figcaption></figure>

and the product free energy. In these QM-GPR/MM studies, GPR models trained using Python libraries such as GPflow²⁰ and GPyTorch²¹ are deployed on the fly during MD simulations (Figure 19B). A Colab-based tutorial is available to demonstrate the use of the related Python libraries to train basic ML models for reactive systems.²² A similar tutorial for training GPR models for QM/MM systems using CHARMM and pyCHARMM is currently under development.

## 11.9. Multistate Empirical Valence Bond (MS-EVB).
The MS-EVB module of CHARMM⁷²³,⁷²⁴ is an efficient method for representing reactive PES, e.g., in enzymes where a system moves from a reactant to a product topology (see ref 725 for implementation details). The most common approach to EVB involves a pseudo-Hamiltonian matrix $H(q)$ constructed from two diabatic reactant (R) and product (P) basis functions, yielding a $2 \times 2$ matrix

$$ \mathbf{H}(\mathbf{q}) = \begin{bmatrix} V_R + \epsilon_R & H_{1,2} \\ H_{1,2} & V_P + \epsilon_P \end{bmatrix} \qquad (61) $$

where $V_R$ and $V_P$ are the potential energies of the reactant and product diabatic states at a given geometry $q$, obtained from standard FF. $\epsilon_R$ and $\epsilon_P$ are constant diagonal energy shifts usually chosen to reproduce the known exo- or endothermicity of the reaction. The off-diagonal element $H_{1,2}$ couples the reactant and product basis functions, which is usually a simple function of atomic coordinates. In CHARMM, one can choose constants or 1D/2D Gaussians, which are functions of one or two distances between atoms. $H$ can be diagonalized into $D = U^T H U$, where the diagonal matrix $D$ contains the eigenvalues and $U$ consists of eigenvectors of $H$. Applying the Hellman-Feynman relation gives a matrix of Cartesian atomic forces

$$ \mathbf{F} = -\frac{d\mathbf{D}}{d\mathbf{q}} = -\mathbf{U}^T \frac{d\mathbf{H}}{d\mathbf{q}} \qquad (62) $$

which contains the gradient vector $F_i$ for each adiabatic state corresponding to the i-th eigenvalue of $D$ in increasing order of energy. $F_0$ contains forces corresponding to the lowest eigenvalue $\lambda_0$, and is used for dynamics propagation on the adiabatic ground state.

The CHARMM-EVB implementation utilizes MPI to parallelize the energy and force calculation for each topological replica at any given time step, achieving near-linear scaling with the number of topological replicas so that the number of topological replicas is limited only by the number of MPI threads running on the given hardware.
**11.10. Reactive MD.** Following chemical reactions in time and space is a central aspect of chemistry. For computer-based methods, ab initio MD methods at correlated levels are usually too prohibitive, in particular if statistically significant numbers of trajectories need to be run. Earlier and previous empirical efforts to describe bond breaking and formation include approaches based on bond order and bond strength.⁷²⁶⁻⁷²⁹ Alternatively, chemical reactivity can be modeled as a linear combination of empirical energy functions describing two or multiple atom connectivities (reactant and one or several products) and to mix these representations. This leads to reactive PESs as in multistate adiabatic reactive MD (MS-ARMD).⁷³⁰⁻⁷³² Here, the PESs are mixed according to

$$V_{\text{MS-ARMD}}(\mathbf{x}) = \sum_{i=1}^{n} w_i(\mathbf{x}) V_i(\mathbf{x}) \qquad (63)$$

The weights $w_i(\mathbf{x})$ are obtained by normalizing the Boltzmann distributed raw weights $w_{i,0}(\mathbf{x})$

$$w_i(\mathbf{x}) = \frac{w_{i,0}(\mathbf{x})}{\sum_{j=1}^{n} w_{j,0}(\mathbf{x})}, \quad w_{i,0}(\mathbf{x}) = e^{-(1/\Delta V)[V_i(\mathbf{x}) - V_{\min}(\mathbf{x})]} \qquad (64)$$

where $V_{\min}(\mathbf{x})$ is the minimal energy scale for a given configuration $\mathbf{x}$ and $\Delta V$ is a characteristic energy scale (switching parameter). Per construction (cf., Eq. 64), only surfaces within a few times of $\Delta V$ from $V_{\min}(\mathbf{x})$ will contribute to instantaneous configuration $\mathbf{x}$. ARMD mixes different PESs $V_i$ by using Gaussian and polynomial functions around the crossing points between states by fitting to reference data such as the MEP.⁷³² Because the mixed PES $V_{\text{MS-ARMD}}(\mathbf{x})$ depends on energies of different states through weights $w_i$ which in turn are analytical functions of the coordinates $\mathbf{x}$, energy-conserving MS simulations can be run using MS-ARMD.⁷³²

A more recent extension combines MS-ARMD⁷³³ with VALBOND, a FF that allows to describe the geometries and dynamics of metal complexes.⁷³⁴⁻⁷³⁶ The formulation is reminiscent of empirical valence bond theory⁶⁰⁰ where diagonal terms are VALBOND descriptions of the states involved and off-diagonal terms describe the orbital overlap. MS-ARMD can also be combined with MM with proton transfer (MMPT),⁷³⁷ to follow proton transfer in gas and condensed phases.⁷³⁸⁻⁷⁴¹

In the gas phase, MS-ARMD was used to study reactions such as hydrogen transfer in the photodissociation of $H_2SO_4 \rightarrow H_2O + SO_3$⁷⁴² and other atmospherically relevant molecules by following vibrational excitation of the OH stretch,⁷⁴³,⁷⁴⁴ the Claisen rearrangement reaction,⁷⁴⁵ or to investigate Diels-Alder reactions.⁷⁴⁶ Such studies provide insights into reaction mechanisms and relevant coordinates driving the process. As an example, for the Diels-Alder reaction between 2,3-dibromo-1,3-butadiene and maleic anhydride MS-ARMD emphasized the importance of rotations of the two reactants to reach the transition state.⁷⁴⁶

More recently, the unimolecular dissociation of vibrationally excited syn-CH₃CHOO to form OH and CH₂CHO was investigated (Figure 20).³⁸²,³⁸⁹ For the reactant and product states, MS-ARMD performs close to the chemical accuracy (~1 kcal/mol) whereas the MEP is described considerably more accurately (inset in Figure 20B). Atomistic simulations using the MS-ARMD PES are about 2 orders of magnitude more efficient than using a neural network-based PES and about 6 orders of magnitude faster than ab initio MD

<figure><img src="image_29.png" /><figcaption>**Figure 20.** (A) OH-elimination following vibrational excitation of syn-CH₃CHOO. (B) Performance of MS-ARMD. The RMSEs between reference ab initio energies and the fitted FF for reactant (blue) and product (green) are 1.1 and 1.2 kcal/mol, respectively. Inset: MEP calculated with MS-ARMD (red circles) compared with reference calculations (black line). (C) Distribution of the total kinetic energy release from several thousand trajectories following CH-excitation with ~2 quanta using the MS-ARMD PES with OO scission energies of 22, 25, and 27 kcal/mol (blue, red, green). Open symbols are experimental results.³⁸⁷ Panels A and B reproduced from ref 382. Copyright [2021] American Chemical Society.</figcaption></figure>

simulations at the MP2 level of theory at which the MS-ARMD PES was developed. In other words, MS-ARMD simulations can be run routinely with high quality and in statistically significant numbers, as exemplified in Figure 20C.

Finally, biological systems were also studied using a combination of RKHS-based PESs and empirical FFs,⁷⁴⁷,⁷⁴⁸ which allowed structural interpretation of metastable states in MbNO and a molecularly refined understanding of ligand exchange (NO vs O₂) at the heme-iron in truncated hemoglobin.

**11.11. Indirect QM/MM Free Energy Simulations.** The alchemical free energy functionality, specifically the PERT module of CHARMM, was described in detail in the 2009 paper.³ In addition to discrete intermediate states as a function of the coupling parameter $\lambda$, PERT supports slow-growth TI (SGTI) that changes $\lambda$ incrementally at each step of the MD simulation. SGTI suffers from the Hamiltonian lag problem,⁷⁴⁹ hence it is rarely used directly. A free energy difference obtained from SGTI should be treated as non-equilibrium work (NEW),⁷⁵⁰ and the equilibrium free energy difference can be obtained by applying the Jarzynski equality⁷⁵¹ or Crooks theorem⁷⁵² to a sufficient number of SGTI runs.⁷⁵³ Such calculations can be automated by CHARMM's scripting language, as illustrated below.

PERT fully supports CHARMM's multiscale capabilities (MSCALE module).²⁶² This makes it possible to compute free energy differences between two descriptions of a system, such as an MM description on one hand and a hybrid QM/MM description on the other hand. Most standard applications of alchemical free energy simulations, such as the calculation of relative binding free energies, employ equilibrium methods (TI,³²⁹ Bennett's acceptance ratio method, BAR,⁷⁵⁴ or its multistate extension MBAR¹²). NEW based methods are also used.¹⁶⁹,⁷⁵⁵ For such traditional applications, it is unclear whether equilibrium or non-equilibrium techniques are more efficient. The situation is different when one has to compute free energy differences between levels of theory, as is the case
for the so-called indirect cycle QM/MM alchemical free energy simulations.

The calculation of free energies when using a QM/MM description poses two challenges: (1) It is slow, making it difficult to achieve sufficient sampling; and (2) standard recipes to realize alchemical transformations, such as soft-core potentials do not work with QM/MM Hamiltonians.²⁹⁷,³⁷⁷,⁷⁵⁶ Indirect cycles can circumvent both issues. The basic idea is illustrated using the calculation of an absolute solvation free energy (Figure 21). Instead of computing ΔG<sub>solv</sub>(QM/MM)

<figure><img src="image_30.png" /><figcaption>Figure 21. Illustration of an indirect cycle to compute a free energy difference at the QM/MM level of theory.</figcaption></figure>

directly, one computes ΔG<sub>solv</sub>(MM) where soft-core potentials, etc., can be used without restrictions. In addition, free energy differences ΔG<sub>α</sub>(MM → QM/MM) (α: gas or aq) both in the gas and aqueous can be calculated, which yields

$$ \Delta G_{\text{solv}}(\text{QM/MM}) = -\Delta G_{\text{gas}}(\text{MM} \rightarrow \text{QM/MM}) + \Delta G_{\text{solv}}(\text{MM}) + \Delta G_{\text{aq}}(\text{MM} \rightarrow \text{QM/MM}) \quad (65) $$

The calculation of ΔG<sub>α</sub>(MM → QM/MM) is challenging,⁷⁵⁷ but NEW-based approaches have been shown to be reliable and efficient.⁷⁵⁸⁻⁷⁶¹

The combination of PERT and MSCALE enables the calculation of NEW values for transitioning from, e.g., a MM to a QM/MM description using solely CHARMM's scripting language. By inserting as few as two hundred of such work values obtained from trivially parallel simulations, into the Jarzynski equality, ΔG<sub>α</sub>(MM → QM/MM) can be calculated accurately and efficiently in most cases.⁷⁶⁰,⁷⁶¹ The NEW switches require equilibrium configurations sampled in the canonical ensemble, i.e., at the low level of theory. Restart files are saved at regular intervals, from which 2–5-ps long independent NEW switching simulations start in parallel.

A self-contained example illustrating this procedure is available at Zenodo.⁷⁶² While the example uses the SCC-DFTB method as the high level theory,¹⁶³ changes required for, e.g., a true DFT method are trivial. The key step is using MSCALE to employ one master (control) job and two slave jobs. The latter are responsible for computing energies and forces at the respective MM and QM/MM levels. The master process is primarily responsible for mixing forces/interactions and integrating the equations of motion. Care is needed to avoid double-counting if additional restraints or similar terms are used. The switching itself is realized as a SGTI calculation of the PERT module. By default, one switches linearly from λ = 0 to λ = 1 in 1,000–5,000 MD steps. The final result of each switch is the NEW value W, which can be saved or extracted automatically within the CHARMM script. Full automation

can be achieved by calling the relevant CHARMM jobs from e.g., the Unix shell, a Python script, etc.

If the convergence of results obtained by the Jarzynski equality⁷⁵¹ is in doubt, one can also use the Crooks theorem.⁷⁵² In this case, a QM/MM-level equilibrium simulation is needed, followed by switches in the QM/MM → MM direction. The latter are again trivially parallel. While the computational cost of generating the initial configurations from the equilibrium QM/MM simulation is high, this workflow is still significantly more efficient than equilibrium-based approaches, which would require adequate sampling at each intermediate state (typically ten or more).

## 12. BOUNDARY CONDITION, SYSTEM PREPARATION, AND TRAJECTORY ANALYSIS

In addition to performing simulation itself, an ability to prepare the simulation system in a desired initial state, impose appropriate boundary conditions or constraints, and analyze simulation trajectories are essential for making scientific discoveries. CHARMM has an extensive set of tools available for this purpose. Presently, at least 25% of more than 1.17 M lines of the CHARMM source code belong to this category. Example applications of the recently developed methods described below are given in references therein. A vast array of other existing tools can be found in the documentation as well as the example 'Testcase' input scripts provided in the CHARMM package.

### 12.1. Simulation and Analysis of Membrane Proteins.

By virtue of its considerable functional flexibility, CHARMM has been a tool of choice in many studies of ion channels and membrane proteins. A theory was developed and implemented to account for the membrane potential and its representation by a constant electric field in computer simulations.⁷⁶³ It was subsequently used in studies of the Kv1.2 potassium channel⁷⁶⁴ and the voltage sensing domain of the voltage-sensitive phosphatase from Ciona intestinalis.⁷⁶⁵ Ion permeation through various channels was characterized,⁷⁶⁶⁻⁷⁶⁹ and the fundamental principles governing ion selectivity were explored.⁷⁷⁰,⁷⁷¹ Computational methods with empirical energy restraints were developed to exploit information from low-resolution experimental data in structural refinement of membrane proteins. A particular attention was given to electron paramagnetic resonance (EPR) accessibility data,⁷⁷²⁻⁷⁷⁴ and double electron-electron resonance (DEER) technique that reports distance distribution between spin labels.⁷⁷⁵⁻⁷⁷⁷ The EPR/DEER methodology is also supported by CHARMM-GUI for easy setup of restrained simulations.⁹⁵ Using these methods, the structure of various ion channels and membrane transporters were refined on the basis of EPR experimental data.⁷⁷²,⁷⁷³,⁷⁷⁸⁻⁷⁸⁰ Energy restraints were also developed to exploit information from mutational cross-link data,⁷⁸¹ which resolved ambiguities about the conformation of the resting state of the voltage sensing domain of potassium channels.⁷⁸²,⁷⁸³

### 12.2. Coordinate Unwrapping for Diffusion Constants.
In CHARMM, fractional coordinates are used to unwrap trajectories, which yields the same diffusion constants as the more recent "toroidal view preserving" method.⁷⁸⁴ It also highlights the need to correct calculated diffusion constants for PBC artifacts, especially those in lipid bilayers. The translational diffusion constant D is typically calculated using the Einstein relation
$$2nD = \lim_{t \to \infty} \frac{\text{MSD}}{t} \qquad (66)$$

where n is the spatial dimension, MSD is the mean-squared displacement, and t is simulation time.⁷⁸⁵ In simulations with finite periodic box, it is standard to “wrap” or “image” positions of molecules such that a molecule crossing a unit cell boundary is translated to the opposite side. This effect must be removed via “unwrapping” for MSD calculation,⁷⁸⁶ which is done in either Cartesian or fractional coordinates.

Constant volume (NVT) or constant energy (NVE) ensembles are recommended for diffusion calculations, as they minimize perturbations to dynamical variables. In MD codes tracking atomic virals like CHARMM, another advantage of using NVT is that the system viscosity and particle diffusion constants can be computed from the same trajectory to compare with experiment.⁷⁸⁸ Unwrapping is straightforward under constant volume, where Cartesian coordinates can be used. However, it can be confounded by volume fluctuations in NPT simulations. For example, the heuristic unwrapping scheme (in which the position of a particle is unwrapped by comparing its current wrapped position to its unwrapped position at the previous time step) used by several MD software packages with Cartesian coordinates was shown to introduce cumulative errors in molecule’s position and calculated MSDs.⁷⁸⁹ A new method called the “toroidal view preserving scheme” was proposed by the same group to correctly unwrap such simulations.⁷⁸⁴

CHARMM avoids the preceding problem by always unwrapping coordinates in fractional space (i.e., the space where each unit cell vector is transformed into 3 orthonormal vectors and positions are mapped onto this lattice) before projecting the coordinates back in Cartesian space. In fractional coordinates, the box fluctuations are removed, and MSD vs t plots with the correct slope are produced. Noise at longer simulation times can be mitigated by projecting the coordinates back into Cartesian space using the average unit cell vectors; the keyword for this operation in CHARMM is XFLUC.

The MSD is often computed using multiple time origins as a difference correlation function, where any deviation is indicative of a problem.⁷⁸⁴ For comparison, Figure 22 shows the single time origin MSD vs t for 1340 TIP3P waters at 20 °C for NPT and NVT simulations. The NPT simulation was unwrapped in four ways: Cartesian coordinates (non-CHARMM), CHARMM fractional coordinates, CHARMM fractional coordinates with average box dimensions (XFLUC), and toroidal view preserving.⁷⁸⁴ The NPT simulation unwrapped with the Cartesian scheme (violet) shows a significant accumulation of error and deviation from linear behavior after 500 ns. The increase in noise with simulation time for fractional space unwrapping (black) is essentially eliminated by projecting back onto the average unit cell vectors (blue). The resulting MSD vs t is practically indistinguishable from the results obtained with NVT (orange) or the toroidal view preserving scheme (green). Diffusion constants obtained from any of the non-Cartesian methods are therefore statistically equal.

Note that even for Cartesian unwrapping, the first several hundred nanoseconds appear unaffected. Comparison of the wrapping frequency distribution for water and for self-diffusion in a DPPC bilayer shows stark differences (Figure 23). Wrapping events for water are 2 or 3 orders of magnitude more

frequent than lipids, and the distribution is approximately Gaussian. The distributions for lipids are Poisson-like with zero being the most frequent value.

Errors obtained from the Cartesian unwrapping scheme are related to the number of wrapping events for a given molecule. Evaluation of lipid diffusion in bilayers from published simulations⁵⁷¹ shows that for larger, slower moving molecules the unwrapping method had very little effect on the results. Table 3 shows that the standard deviation over five replicate simulations (bottom row) is an order of magnitude larger than that over the three different unwrapping methods (last column).

Diffusion constants are also affected by the PBC that causes underestimates of the infinite system-size value. PBC errors in diffusion constants in isotropic systems are relatively modest, approximately 10%, and can be corrected by the Yeh-Hummer formula.⁷⁹⁰ The PBC correction for diffusion in lipid bilayers requires the periodic Saffman–Delbrück model.⁷⁹¹,⁷⁹² It is quite large for lipid self-diffusion, approximately 3-fold for

<figure><img src="image_31.png" /><figcaption>Figure 22. MSD vs t for 1-μs NVT and NPT simulations of 1340 TIP3P waters. The plots are offset by intervals of 0.25 × 10⁶ Å² on the y-axis to better distinguish them. From top to bottom, violet: NPT with Cartesian unwrapping; green: NPT with toroidal view preserving scheme; orange: NVT with Cartesian unwrapping; blue: NPT with unwrapping in fractional space and using average unit cell vectors (CHARMM XFLUC); black: NPT with unwrapping in fractional space (CHARMM method). Trajectories were run with OpenMM and analyzed with CPPTRAJ version 6.19.3.⁷⁸⁷</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Figure 23. Probabilities of wrap counts for water molecules from a 1-μs NPT MD simulation of TIP3P water and for a bilayer containing 288 DPPC molecules from 5 replicate 400-ns MD simulations.⁵⁷¹ The inset omits the large peak for DPPC with zero wraps; the mean number of lipid wrap counts is 111.5.</figcaption></figure>
<table><caption>**Table 3.** *D*<sub>MSD</sub> (in Units of 10<sup>-7</sup> cm<sup>2</sup>/s) for NPT DPPC Bilayers with Drude2023 FF, with Cartesian (Cart.), Fractional (Frac.), and Toroidal (Tor.) Unwrapping</caption><thead><tr><th>Rep #</th><th>Cart.</th><th>Frac.</th><th>Tor.</th><th>Avg</th><th>Std</th></tr></thead><tbody><tr><td>1</td><td>0.71</td><td>0.71</td><td>0.70</td><td>0.71</td><td>0.001</td></tr><tr><td>2</td><td>0.73</td><td>0.73</td><td>0.73</td><td>0.73</td><td>0.001</td></tr><tr><td>3</td><td>0.74</td><td>0.74</td><td>0.75</td><td>0.74</td><td>0.001</td></tr><tr><td>4</td><td>0.67</td><td>0.67</td><td>0.67</td><td>0.67</td><td>0.003</td></tr><tr><td>5</td><td>0.69</td><td>0.69</td><td>0.69</td><td>0.69</td><td>0.003</td></tr><tr><td>Avg</td><td>0.71</td><td>0.71</td><td>0.71</td><td>0.71</td><td>0.001</td></tr><tr><td>Std</td><td>0.03</td><td>0.03</td><td>0.03</td><td></td><td></td></tr></tbody></table>

bilayers with 288 lipids, and should not be overlooked. Our overall recommendations for unwrapping are the following:

* The Cartesian-based method should not be used when the box dimensions can change during a simulation. Though for relatively short simulations of slowly diffusing particles the errors are not substantial (Table 3), accurate methods are readily available and should be used.

* The XFLUC method in CHARMM should only be used when the box dimension fluctuations are less than approximately 15%. This is its major limitation. Though XFLUC would still yield smooth plots, the slope (and hence D) could be incorrect if the aspect ratio changes significantly.

* Although the toroidal view preserving method provides a smoother result than with CHARMM fractional coordinates, both methods should give correct answers if the average slope is used.

**12.3. Calculation of Pressure Profiles in Lipid Bilayers.** The lateral pressure profile provides a detailed view of the forces within a planar lipid bilayer with respect to the average bilayer normal (typically z-axis). CHARMM provides capability for estimating the spontaneous curvatures in symmetric<sup>557</sup> and asymmetric bilayers<sup>44</sup> and the difference in leaflet surface tensions (or differential stress) in asymmetric bilayers.<sup>44,793</sup> For a simulation system Ω, the volume-averaged virial stress tensor, σ, is<sup>794</sup>

$$ \sigma = -\frac{1}{V} \left( \sum_{i \in \Omega} m_i \vec{v}_i \otimes \vec{v}_i + \sum_{\{i<j\} \in \Omega} \vec{r}_{ij} \otimes \vec{r}_{ij} \right) \quad (67) $$

where V is the volume of Ω and mᵢ and vᵢ are the mass and velocity of atom i. The symbol ⊗ represents the tensor product between two vectors, rᵢⱼ = rᵢ - rⱼ, and Fᵢⱼ is the force exerted from atom j on atom i. The first and second terms on the right-hand side of Eq. 67 are the kinetic and configurational contributions, respectively. In the configurational virial stress, the contributions from periodic images must be considered for non-bonded interactions. The pressure tensor (p) is defined as the negative of the stress tensor, i.e., p = -σ. For an isotropic system, the bulk pressure is P = -Tr(σ)/3. CHARMM calculates the virial stress (and thus the pressure) using Eq. 67. While constant pressure can be emulated using Monte Carlo barostat,<sup>34</sup> the virial is needed to calculate transport properties such as viscosity.<sup>788</sup>

Derivatives of the free energy with respect to virtual transformations<sup>795</sup> of the periodic box can be computed with the lateral pressure profile p<sub>L</sub>(z) typically along the z-direction, normal to the membrane surface

$$ p_L(z) = p_T(z) - p_N(z) \quad (68) $$

where p<sub>T</sub>(z) and p<sub>N</sub>(z) are the tangential and normal components of the pressure tensor, respectively. The zeroth moment of the profile is the tension (the derivative of the free energy with respect to area) while the first moment is the derivative of the free energy with respect to the curvature. Interpreted through the Helfrich/Canham Hamiltonian,<sup>796,797</sup> it provides a convenient route to calculate spontaneous curvature of the leaflet.<sup>798</sup>

For planar lipid bilayers, the slab geometry is a convenient choice for lateral pressure profile calculations

$$ p_T = (p^{xx} + p^{yy})/2, \quad p_N = p^{zz} \quad (69) $$

Here, p<sup>xx</sup>, p<sup>yy</sup>, and p<sup>zz</sup> are the diagonal elements of p(z). Since a typical lipid bilayer cannot support in-plane shear strain, there is no off-diagonal coupling of x and y.

Denoting the z-dimension of the simulation box as L<sub>z</sub>, the bilayer and leaflet surface tensions are given as zeroth-moments of p<sub>L</sub>(z):

$$ \begin{gathered} \gamma = -\int_{-L_z/2}^{+L_z/2} dz p_L(z), \\ \gamma_1 = -\int_0^{+L_z/2} dz p_L(z), \quad \gamma_2 = -\int_{-L_z/2}^0 dz p_L(z) \end{gathered} \quad (70) $$

Without any external force, the bilayer tension must vanish (γ = 0). While leaflets are tensionless (γ₁ = γ₂ = 0) in symmetric bilayers, leaflet tensions and their difference (Δ = γ₁ - γ₂) in asymmetric bilayers do not necessarily vanish.<sup>799</sup>

The leaflet spontaneous curvature c₀ of a (planar) symmetric bilayer can be calculated from the first moment of p<sub>L</sub>(z)

$$ k_c c_0 = \int_0^{L_z/2} dz p_L(z) \quad (71) $$

where k<sub>c</sub> is the bending modulus of the leaflet. Spontaneous curvatures of asymmetric bilayers can be calculated with a generalization of Eq. 71.<sup>44,800</sup>

While the bulk pressure tensor can be readily calculated using Eq. 67, the calculation of the local pressure tensor (including the lateral pressure profile) is complicated by the need to assign each contribution to the virial locally in space (here, z). Briefly, a vector contour is integrated from one force center i to the other center j, creating a 3D function whose gradient is a Dirac-δ function at each end point multiplied by the force. With this requirement met, the pressure tensor contains spatial correlations of force such that the work required to reshape the box by a virtual deformation can be computed. Many choices of contour satisfy this requirement, leading to the inherent ambiguity of the profile unless complemented by a virtual deformation yielding an observable that resolves the ambiguity.<sup>795</sup> Because of this complexity, the pressure profile calculation has not been supported in most other simulation programs.

The LOPR module in CHARMM calculates p<sub>T</sub> in Eq. 69 where their profiles from the full electrostatics can be obtained by either the Ewald sum<sup>801</sup> or by the PME method.<sup>176</sup> Presently there are only two programs in addition to CHARMM that support the pressure profile from full electrostatics: NAMD (by Ewald sum in post analyses) and GROMACS (by PME<sup>802</sup> in a branch version). Due to the limitation of the Harasima contour<sup>801</sup> employed for the virial
calculation, the normal component $p_N$ (Eq. 69) is not calculated but can be evaluated in post analysis. Furthermore, a planar bilayer cannot support heterogeneous normal pressure. For a bilayer without any external force ($\gamma = 0$), $p_N$ is calculated from Eqs. 68 and 70 as

$$p_N = \frac{1}{L_z} \int dz p_T(z) \qquad (72)$$

The pressure profiles $p^{xx}(z)$ and $p^{yy}(z)$ are typically calculated by binning the z-dimension.⁸⁰³ For accurate profiles, the bin size is set typically to ~1 Å. Alternatively, $p^{xx}$ and $p^{yy}$ can be calculated in a binless manner using Fourier series

$$\begin{aligned} p^{\alpha\alpha}(z') &= \frac{w_{00}}{2} + \sum_{n=1}^{N} \sum_{l=0}^{1} w_{nl} \cos(2\pi nz' - \phi_l) \\ w_{00} &= 2 \int_{-1/2}^{+1/2} dz' p^{\alpha\alpha}(z') \\ w_{nl} &= 2 \int_{-1/2}^{+1/2} dz' p^{\alpha\alpha}(z') \cos(2\pi nz' - \phi_l) \end{aligned} \qquad (73)$$

where $\alpha$ represents x or y and $z' = (z - z_{cm})/L_z$ is the fractional z-coordinate with respect to the bilayer center $z_{cm}$. N is the order of the Fourier series, $w_{nl}$ are Fourier coefficients, and $\phi_0 = 0$ and $\phi_1 = \pi/2$ are phase shifts for even and odd series. The LOPR module supports both methods for the lateral pressure profile calculation. The Fourier series method is currently supported only in CHARMM, where accurate pressure profiles can be obtained with a moderate number of coefficients, $N \sim 20$ for typical bilayers with $L_z \sim 80$ Å. For simulation systems with larger $L_z$, larger N is required.

The LOPR module supports both on-the-fly and post analysis calculation of $p_T(z)$ including the full electrostatics via the PME method. This allows efficient resampling where sparsely sampled coordinate and velocity trajectories from various programs including CHARMM, NAMD, and OpenMM (with a customized Velocity Reporter for trajectory generation)⁷⁹³ can be utilized (Figure 24A,B). In the resampling approach, multiple short CHARMM simulations for chosen frames can be run simultaneously, which can be easily realized in typical computational resources. If there are sufficient samples from long simulation times, one can also calculate the lateral pressure profile in a post analysis using a single-step dynamics for each frame with an integration time step shorter than the one used for the original simulations. The post analysis method is faster than the resampling method, and yields sufficiently accurate results from 500-ns trajectories saved at every 10 ps (Figure 24B). For an asymmetric bilayer, the pressure profile is also asymmetric (Figure 24C) from which the non-vanishing leaflet tensions are calculated from Eq. 70.

The LOPR module currently does not support LJ-PME⁵⁶⁴ (Section 8.5) and the polarizable Drude model⁵⁷¹ which will be supported in future updates. Additionally, it is implemented only for CPU calculations without DOMDEC²⁵ (Section 2.2), which results in poor scalability over multiple nodes. Thus, the current LOPR module is suitable practically for only single node-jobs, and parallelization of the code will greatly improve its performance.

### 12.4. P2₁ Periodic Boundary Condition.
A novel approach utilizing the P2₁ PBC has been implemented to relax differential stress between the leaflets during MD

<figure><img src="image_33.png" /><figcaption>**Figure 24.** (A) Comparison of pressure profiles from a 100-ns CHARMM simulation of palmitoylsphingomyelin bilayer with complete sampling (Direct) or via the 10% resampling (Resamp) of 100-ps intervals spaced 1-ns apart, using restart files from the fully sampled simulation. The values of the first moment and their standard errors are 0.190 ± 0.009 and 0.192 ± 0.004, respectively. (B) Pressure profiles of a bilayer composed of 72 1,2-dipalmitoyl-sn-glycero-3-phosphocholine at 1 atm and 323 K. Data from resampling were obtained by the binning method with 100 bins from the previous CHARMM simulations (black), where 100-ps resampling was done for every 1 ns.⁵⁵⁸ Post analysis using single-step dynamics was tested using 500-ns OpenMM NPT trajectories of the same bilayer, where the pressure profiles were calculated for each frame by both the binning (blue) and Fourier series (red) methods. The numbers of bins and Fourier coefficients were set to 80 and N = 20 (Eq. 73), respectively. Standard errors from five 100-ns blocks are shown in cyan and pink areas, respectively, which are smaller than the line thickness except near dips. (C) Pressure profiles for a 1,2-dinervonoyl-sn-glycero-3-phosphocholine bilayer with a model peptide of 9 monomeric units of gramicidin A in the upper leaflet at a peptide area fraction $\phi \sim 0.40$ before (red) and after P2₁ equilibration (black). The P2₁ PBC in CHARMM allows lipid translocation between bilayer leaflets, which reduces area stress and alters the lateral pressure profile (see Section 12.4). The bilayer midplane was set to z = 0. Pressure profiles from the last 300-ns trajectories of 5 replicas from OpenMM were averaged for each asymmetric bilayer, whose standard errors are shown as pink and gray areas (typically smaller than the line thickness). Pressure profiles were calculated by single-step dynamics for each frame of trajectories with 200 bins along the z-direction. Panel A reproduced from ref 804. Copyright [2014] Cell Press. Panel C reproduced from ref 44. Copyright [2023] Wiley.</figcaption></figure>

simulations of lipid bilayers. The inherent difficulty in accurately estimating the number of lipids in each layer *a priori* gives rise to the differential stress. The widely used P1 PBC involves tessellating the simulation space with translated images of the box. As an atom exits the simulation box, it is replaced by its image located on the opposite face. Unlike P1, the P2₁ PBC introduces a half-screw symmetry between the images. In this scheme, the image of the simulation box is not a mere translated copy but a 180°-rotated image translated along the same screw axis. In the CHARMM non-DOMDEC version (DOMDEC is explained in Section 2.2), the screw axis could be oriented in any direction. However, with the Extended Eighth Shell (EES) method⁴³ in DOMDEC, the screw axis is constrained to the x-axis for enhanced performance. The P2₁ symmetry operation is denoted as $(x + 1/2, -y, -z)$, representing a half-unit cell length translation along the x-axis and reflection along the y- and z-axes. Reflection along two perpendicular axes is equivalent to 180° rotation along the
screw axis where lipids departing from the top layer along the x-axis (yz-faces) reenter the cell in the bottom layer, and vice versa. However, since only the x-axis is allowed as the screw axis, lipids leaving the cell along the xz-faces reenter along the same leaflet. This can be visualized as a torus-shaped structure along the screw axis where the top layer transitions to the bottom layer in the neighboring image cell. A bilayer simulation starting with 108 lipids on the top layer and 92 lipids in the bottom equilibrates to 100 lipids in both layers using this method.<sup>43</sup>

The scaling performance of the EES method for P2₁ simulations is similar to that for DOMDEC in P1 simulations. Compared to the latter, there is slightly larger import volume during the message transfer among nodes during the direct space calculations. However, by restricting the screw axis to the x-axis, there is a minimal impact of the larger message size on the overall performance. The reciprocal space calculations are done by distributing the charge on the full unit cell. However, forces are calculated only in the asymmetric unit, and extra bookkeeping is also performed to rotate the forces and velocities as the images are 180°-rotated along the yz faces.

While it may first appear that the exchange of lipids between the layers would allow only symmetric bilayer simulation, the P2₁ PBC is specifically useful for setting up asymmetric bilayers by restraining specific lipids to their original leaflet and redistributing others between the leaflets through P2₁ PBC. The differential stress or the difference in surface tension between the two layers is a consequence of the intrinsic bending and the asymmetric lipid packing. Methods for simulating asymmetric bilayers can be categorized as lipid-based, leaflet-based, and bilayer-based.<sup>793</sup> In the lipid-based approach termed APL, surface areas of the top and bottom are matched using the area per lipid from homogeneous lipid bilayers. This is the simplest and the most commonly used approach that disregards any coupling between the bilayers and assumes ideal mixing of lipids. The leaflet-based approach termed SA minimizes the differential area strain between the leaflets. It first equilibrates symmetric bilayers corresponding to each leaflet composition and then combines one leaflet from each bilayer. While this removes the differential strain, it also disregards coupling between the two leaflets. The bilayer-based approach termed 0-DS removes the differential stress by adjusting the number of lipids. When these approaches are followed by P2₁, agreement in mechanical properties significantly improves among APL/P2₁, SA/P2₁, and 0-DS/ P2₁.<sup>793</sup> These findings align with a theoretical framework emphasizing the intricate interplay between bending and asymmetric lipid packing. Torque balance conditions and stress indices provide theoretical support, showcasing promising results for P2₁ simulations in capturing lipid asymmetry observed in biological membranes.

As another example of the importance of P2₁ for setting up bilayer simulations, the curvature induced by a peptide in the so-called “peptide-asymmetric bilayer” was studied.<sup>44</sup> In these simulations, while both layers contained the same types of lipids, the asymmetry was induced by the presence of the peptide in the cis-leaflet. A series of gramicidin A (gA)-based peptides were simulated: a single monomer, a fused tetramer, a fused nonamer, nine gA monomers, and a fused tetramer of a gA mutant whose Trp residues were replaced by Gln. These assemblies were used to investigate effects of the size of chemically similar peptides spanning a single leaflet. Utilizing the APL method mentioned above, systems were created at

three distinct peptide area fractions. Subsequently, equilibration employing the P2₁ PBC was performed. The P1 simulation showed significant condensation of lipids in the trans leaflet, resulting in large differential stress and bilayer bending moment, which relaxed via the exchange of lipids between leaflets in P2₁ simulation (Figure 24C).

**12.5. Primary Hydration Shell (PHS) Model.** While the representation of solvent surrounding macromolecules should closely approximate physical reality in MD simulations, the common use of a sizable volume of solvent with PBC is computationally expensive. A PHS consisting of 2–3 layers of explicit water molecules around a protein may be sufficient to maintain the conformational stability and dynamics of macromolecules. The initial work<sup>631</sup> on the PHS model has been refined in two stages. First, the method was tested with hen egg lysozyme, where good agreement with the protein and solvent behavior was observed, including Lipari–Szabo order parameters for N–H main-chain and N–H<sub>2</sub> side-chain motions on the ps–ns time scale in simulations of 25–150 ns, compared to full PBC treatments.<sup>805</sup> The original PHS method has a modest half-harmonic restraint of waters to their nearest protein atom, should they become more distant than a threshold value (5.8 Å by default). As a part of this work a simpler GEO restraint has been implemented which saves computer time as it is calculated relative to three perpendicular principal axes that follow the protein frame. The GEO approach follows global conformational changes and is similarly good for the lysozyme tested.

The PHS method was subsequently refined to overcome issues when applied to larger systems.<sup>806</sup> A neighbor list was implemented to efficiently track the nearest protein atoms to the water oxygen atom, to avoid calculating distances at each step. Also, an asymmetric harmonic potential instead of a half-harmonic one was used to ensure correct water density close to the boundary. In addition, pressure control was implemented, and the confining potential was scaled to keep waters near hydrophobic residues. This approach showed a 14-fold reduction in computing time for a 82-kDa protein. Future developments of the PHS model should handle situations with extensive structural changes<sup>807</sup> and association between proteins where treatment of long-range force is important.<sup>808</sup>

**12.6. Hydration Map.** Surface hydration of proteins and nucleic acids are important for their biological function,<sup>809</sup> which has long been a subject of computer simulation.<sup>810</sup> For measuring location-dependent average behavior of water molecules near biomolecular surfaces, the COORDinates SMAP (Solvation MAP) command has been implemented in CHARMM.<sup>811,812</sup> It divides the simulation water box into a grid of cubic cells (default size: 0.7 Å<sup>3</sup>, half the radius of a water molecule), and locally calculates time-averaged properties of water within each cell. To account for protein motion during simulation, coordinate frames are aligned to a reference structure so that the calculated map is relative to the surface. In the current implementation, local water density (SDENSity), translational diffusion coefficient (DIFTrans), and the average number of hydrogen bonds formed by water molecules (HBOND) can be calculated. By default, water oxygen atoms are used for calculation. Other atoms such as ions can be selected as ‘solvent’ atoms, to build the corresponding maps. This capability has been used to find the preferred location of sodium ions around a double-stranded DNA as the center of the minor groove, between the ‘double water spines.’<sup>812</sup>
The solvation map can be saved as a data file for further analysis, or an MRC format electron density map for visualization, e.g., by using the UCSF ChimeraX.813 Since coordinate frames are oriented to a reference structure, calculated water densities around flexible loops may become low. To examine hydration around a moving loop, a separate solvation map can be built by selecting only the loop as the orientational reference.

The water density map can be used to calculate the solvation free energy.814 For N cells surrounding the protein under consideration, if the water density of cell i during the simulation is $ρ_i$ and denoting the bulk water density as $ρ_b$, the free energy of water for N cells is

$$G_{\text{solv}} = -k_{B}T \sum_{i}^{N} \ln(\rho_{i}/\rho_{b}) \qquad (74)$$

Since the size of a cell is smaller than that of a water molecule (0.7 Å by default), a high water density at cell i means that the cell is visited more frequently rather than water molecules are packed more tightly. The underlying idea for Eq. 74 is that a more frequently visited cell has a lower free energy (a favorable location) compared to the bulk water. To use Eq. 74 in practice, cells corresponding to the first hydration shell are selected, for which a distance cutoff of 4.5 Å from heavy atoms of the protein, and water density cutoff of 0.034 Å⁻³ (cf., $ρ_b = 0.0333$ Å⁻³) are used.814 Compared to the popular grid inhomogeneous solvation theory (GIST),815,816 the above method does not require the protein to be constrained during simulation, where constraint on proteins drastically alters surface hydration. And the measured solvation energy values are in physically more reasonable range compared to those from GIST.814 Calculation of the density-based solvation free energy is being implemented in CHARMM as the COORDinates SMAP SLVE command.

### 12.7. Conformational Entropy.
Conformational entropy is an essential component of the conformational free energy of a biomolecule. A widely used class of entropy calculation methods rely on quasi-harmonic approximation where frequencies of different vibrational modes are used to estimate conformational entropy.817,818 However, they cannot account for transitions between states. With increases in computational power, more direct evaluation of entropy from distributions of degrees of freedom (DOF), in particular, backbone and side- side chain dihedral angles has become possible. If all N DOFs are treated independently, the cost of entropy calculation scales linearly with N. However, since DOFs can be mutually correlated (e.g., by correlated motion of side chains forming contacts), higher order corrections should be made, which amounts to calculating multivariate histograms. The maximal information spanning tree (MIST) approach systematically handles higher order corrections in such a way that the estimated entropy monotonically approaches its asymptotic value.819,820 It should be noted that, an accurate calculation of higher order terms requires a greater number of coordinate frames since for n-th order correlation, the total number of bins for the histogram scales with $N^n$. With a limited number of coordinate frames, bins will be sparsely populated, leading to an increased statistical uncertainty. In practice, MIST calculation up to the second order is good for most purposes since it provides statistically reliable result, fast to calculate, and correlations beyond the second order do not usually contribute significantly.

By limiting DOFs to an amino acid side chain, MIST can be used to estimate the side-chain entropy of individual residues. In this case, since N is small (e.g., total number of dihedral angles of a side chain), calculation up to the third order MIST can be done for higher accuracy. This method has been used to calculate changes in the side-chain entropy for binding of proline-rich ligands to an SH3 domain.821 In this study, “entropy hotspots” were identified where the side-chain entropy of remote residues in the SH3 domain increases upon ligand binding. This arises from the rearrangement of contacts across the protein’s surface that makes the side chains of entropy hotspots to become more mobile upon ligand binding. While initially developed as a separate code,821 the COORDinates MIST command is currently being implemented in CHARMM.

### 12.8. Identifying Non-Polar Contacts.
In CHARMM, hydrogen bonds can be readily identified by the COORDinates HBONd command, which works either for a single structure or across coordinate frames. In comparison, identifying non-polar contacts has been less established. It is often desirable to determine non-polar contacts at the level of individual residues rather than between pairs of atoms. Measuring distances between Cα atoms or centers of mass between side chains of non-polar residues and applying an ad hoc cutoff distance does not provide an accurate picture of non-polar contacts. CHARMM now has the COORDinates DISTance RESIdue command. Without the RESIdue keyword, all pairwise distances between two groups of selected atoms are reported. With the RESIdue keyword, pairwise minimum distances between residues in the two groups (e.g., between two domains) are reported. To identify non-polar contacts, selecting atoms with the absolute value of charges less than $0.3e$ ($e = 1.6 \times 10^{-19}$ C) and distance cutoff of 3 Å can be used. These are based on charges of non-polar hydrogen atoms and their van der Waals radii (1.32 Å). In this way, physical contacts between non-polar residues can be identified and further processed to analyze their dynamics, i.e., occupancy, formation, and breakage.822,823

### 12.9. Vectorial Analysis of Long-Range Concerted Motions in MD Trajectories.
Long-range concerted motions in proteins and other biomolecules is best captured with a correlation coefficient (CC) based on covariance using Euclidean distances between entries of a position-vector time series. This coefficient, DCOR, is a vector equivalent of Pearson’s CC or a generalized CC.824 The relative accuracy of DCOR is established by an assessment conducted using vector displacements generated with a known CC.825 DCOR is least sensitive to angular variation between two vectors compared to Pearson’s CC or a vector CC. Nor is DCOR as sensitive to large variations between vector components compared to the (scalar) generalized CC, which was found to give inflated CCs relative to the actual values when only one of the vector components is highly correlated.

The DCOR value between any two vector time series can be computed using the CORREL module. The vector dimensions need not be equal. For each time series, a matrix of intravector Euclidean distances between all pairs of time points in the series is used to calculate the covariance. DCOR reflects both linear and non-linear correlations,825 and can detect long-distance concerted motions that neither Pearson’s CCs nor the generalized CCs can reveal.826

### 12.10. Other Updated Preparation and Analysis Features. 12.10.1. System Generation.
The residue sequence
to be used in the generation of a segment in the PSF can be
read from the ATOM and/or HETATM records in a PDB file
with a specified chain or segment ID to the READ SEQU PDB
command. If there are residues in the PDB file that do not exist
in the RTF (Residue Topology File), or with names that differ
from those used in the RTF, it is possible to either skip those
residues or map the name in the PDB to the name in the RTF.
The sequence can also be read from the SEQRES records,
which is useful when there are missing residues in the PDB-file.
These enhancements allow the generation of a PSF directly
from a PDB-file without editing it.

**12.10.2. Trajectory Handling.** CHARMM can read trajectories that are contained in multiple files, and normally applies a number of checks to ensure that the set of files constitutes a valid contiguous trajectory with no overlaps or gaps. Sometimes it is desirable to override these checks, for instance, to analyze a set of independent replicate trajectories together to obtain overall statistics. These checks can be disabled, making it possible to mix files that are not contiguous or differ in a other ways (e.g., time step or coordinate saving frequency), as long as they use the same PSF.

Binary trajectories can be read automatically, irrespective of
big-endian or little-endian format of the trajectory (and of the
executable). The CHARMM file OPEN command also allows
endianness to be specified with a keyword.

**12.10.3. Time Series Analysis in the CORREL Module.**
Time series data can be mapped to specified interval, which is useful, e.g., to avoid spurious jumps in dihedral angles (MANTIM command). Time series of protein secondary structure content can also be extracted using CORREL.

**12.10.4. Similarity Analysis of Snapshots from Trajectory Files.** The RMSDYN command now allows different numbers of frames in two trajectories to be compared. Two new metrics have been added: the interatomic average coordinate difference without superposition (DIFF option), and the RMS distance (DRMS option), which compares interatomic distances in one structure with the corresponding distances in another structure, obviating the need for structural superposition. The latter is also available for single coordinate sets, and for time-series analysis in the CORREL module.

13. CONCLUDING DISCUSSION

Since its first publication in 1983,⁴ CHARMM has been
continuously developing as the need and demand for
computational biophysics and biochemistry grew. The present
review of the major developments since 2009³ highlight
improvements as well as new capability, many of which are
uniquely available in CHARMM. This review may thereby
serve as a guide for exploring new methods in addition to
providing a broad overview of the current state of the art.

The advances also reflect changes in the research landscape
at large. Faster simulation engines are needed as the molecular
systems to study are becoming larger and also as the computer
hardware continues to develop. With its modularity and
flexibility, CHARMM now employs a number of engines either
within the program or through APIs for external engines,
which include DOMDEC, BLADE, CHARMM/OpenMM, and
the newly developed apoCHARMM (Section 2). Thus, the
multicore/multithread scaling and speed of CHARMM should
be comparable to those of other fast simulation engines
presently available, while maintaining highest accuracy.

Accessibility is another practical issue, for which CHARMM
is now readily available for academic and non-profit

laboratories.⁸²⁷ While the powerful CHARMM scripting
language enables sophisticated tasks, a potential downside is
the steep learning curve and difficulty in programming. This is
being addressed through the development of pyCHARMM
(Section 3.1). Its Python-based workflow also allows leveraging
the capability of the Python language. Additionally,
pyCHARMM is beginning to serve as a teaching platform
from which the general principles and ideas of molecular
modeling and biomolecular simulation can be taught.
Preparation of the simulation system and the CHARMM
script can also be done through CHARMM-GUI (Section 3.3).

Among other significant features of CHARMM are a wide
range of docking and sampling methods (Section 4—Section 7)
and various energy functions including implicit solvent and
membranes, coarse graining, as well as a host of constraint
capabilities (Section 8—Section 9). The QM/MM methods
described in Section 11 are uniquely available in CHARMM.
Likewise, CHARMM has distinct capabilities in system
preparation, structure manipulation, and coordinate/trajectory
analysis (Section 12). Finally, the ever-expanding CHARMM
FF (Section 10) is becoming the *de facto* standard that is
widely adopted in other simulation packages.

The extensive capabilities of CHARMM enable simulations
and quantitative analyses of systems ranging from small
molecules to large biomolecular assemblies and membrane
systems at both atomistic and CG levels. Beyond studying
small model systems, CHARMM is now increasingly used to
tackle problems of practical importance that involve larger
sizes, longer simulation times, and more extensive sampling,
which will continue to grow with advances in computer
hardware and methodologies. For the latter, CHARMM has
been the testbed for new computational methods, thereby it
stays on the forefront of biomolecular modeling and simulation
with a fertile link to its developers and users. We anticipate
CHARMM will continue to play an essential role for
addressing current problems and also for opening new avenues
of research in biomolecular systems.

<figure><img src="image_34.png" /></figure>

AUTHOR INFORMATION

**Corresponding Authors**

Wonmuk Hwang – Department of Biomedical Engineering,
Texas A&M University, College Station, Texas 77843,
United States; Department of Materials Science and
Engineering and Department of Physics and Astronomy,
Texas A&M University, College Station, Texas 77843,
United States; Center for AI and Natural Sciences, Korea
Institute for Advanced Study, Seoul 02455, Republic of
Korea; orcid.org/0000-0001-7514-3186; Email: hwm@
tamu.edu

Charles L. Brooks, III – Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109, United States; orcid.org/0000-0002-8149-5417; Email: brookscl@umich.edu

Bernard R. Brooks – Laboratory of Computational Biology,
National Heart Lung and Blood Institute, National Institutes
of Health, Bethesda, Maryland 20892, United States;
orcid.org/0000-0002-3586-2730; Email: brb@
nhlbi.nih.gov

Authors

Steven L. Austin – Department of Chemistry, University of South Florida, Tampa, Florida 33620, United States
Arnaud Blondel – Institut Pasteur, Université Paris Cité, CNRS UMR3825, Structural Bioinformatics Unit, F-75015 Paris, France; ● orcid.org/0000-0001-9586-3025

Eric D. Boittier – Department of Chemistry, University of Basel, CH-4056 Basel, Switzerland

Stefan Boresch – Faculty of Chemistry, Department of Computational Biological Chemistry, University of Vienna, 1090 Vienna, Austria; ● orcid.org/0000-0002-2793-6656

Matthias Buck – Department of Physiology and Biophysics, Case Western Reserve University, School of Medicine, Cleveland, Ohio 44106, United States

Joshua Buckner – Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109, United States

Amedeo Caflisch – Department of Biochemistry, University of Zürich, CH-8057 Zürich, Switzerland; ● orcid.org/0000-0002-2317-6792

Hao-Ting Chang – Institute of Bioinformatics and Systems Biology, National Yang Ming Chiao Tung University, Hsinchu 30010, Taiwan, ROC

Xi Cheng – Shanghai Institute of Materia Medica, Chinese Academy of Sciences, Shanghai 201203, China; ● orcid.org/0000-0003-3735-645X

Yeol Kyo Choi – Department of Biological Sciences, Lehigh University, Bethlehem, Pennsylvania 18015, United States; ● orcid.org/0000-0002-4218-7139

Jhih-Wei Chu – Institute of Bioinformatics and Systems Biology, Department of Biological Science and Technology, Institute of Molecular Medicine and Bioengineering, and Center for Intelligent Drug Systems and Smart Bio-devices (IDS²B), National Yang Ming Chiao Tung University, Hsinchu 30010, Taiwan, ROC; ● orcid.org/0000-0003-3842-2893

Michael F. Crowley – Renewable Resources and Enabling Sciences Center, National Renewable Energy Laboratory, Golden, Colorado 80401, United States; ● orcid.org/0000-0001-5163-9398

Qiang Cui – Department of Chemistry, Boston University, Boston, Massachusetts 02215, United States; Department of Physics, Boston University, Boston, Massachusetts 02215, United States; Department of Biomedical Engineering, Boston University, Boston, Massachusetts 02215, United States; ● orcid.org/0000-0001-6214-5211

Ana Damjanovic – Department of Biophysics, Johns Hopkins University, Baltimore, Maryland 21218, United States; Department of Physics and Astronomy, Johns Hopkins University, Baltimore, Maryland 21218, United States; Laboratory of Computational Biology, National Heart Lung and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892, United States; ● orcid.org/0000-0001-7635-2705

Yuqing Deng – Shanghai R&D Center, DP Technology, Ltd., Shanghai 201210, China

Mike Devereux – Department of Chemistry, University of Basel, CH-4056 Basel, Switzerland; ● orcid.org/0000-0002-1561-1635

Xinqiang Ding – Department of Chemistry, Tufts University, Medford, Massachusetts 02155, United States; ● orcid.org/0000-0002-4598-8732

Michael F. Feig – Department of Biochemistry and Molecular Biology, Michigan State University, East Lansing, Michigan 48824, United States; ● orcid.org/0000-0001-9380-6422

Jiali Gao – School of Chemical Biology & Biotechnology, Peking University Shenzhen Graduate School, Shenzhen,

Guangdong 518055, China; Institute of Systems and Physical Biology, Shenzhen Bay Laboratory, Shenzhen, Guangdong 518055, China; Department of Chemistry and Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55455, United States; ● orcid.org/0000-0003-0106-7154

David R. Glowacki – CiTIUS Centro Singular de Investigación en Tecnoloxías Intelixentes da USC, 15705 Santiago de Compostela, Spain; ● orcid.org/0000-0002-9608-3845

James E. Gonzales, II – Department of Biomedical Engineering, Texas A&M University, College Station, Texas 77843, United States; Laboratory of Computational Biology, National Heart Lung and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892, United States; ● orcid.org/0009-0009-7696-7561

Mehdi Bagerhi Hamaneh – Department of Physiology and Biophysics, Case Western Reserve University, School of Medicine, Cleveland, Ohio 44106, United States; Present Address: National Library of Medicine, National Institutes of Health, Bethesda, Maryland 20894, United States; ● orcid.org/0009-0000-5222-889X

Edward D. Harder – Schrodinger, Inc., New York, New York 10036, United States; ● orcid.org/0000-0001-9213-787X

Ryan L. Hayes – Department of Chemical and Biomolecular Engineering and Department of Pharmaceutical Sciences, University of California, Irvine, Irvine, California 92697, United States; ● orcid.org/0000-0003-1052-6391

Jing Huang – Key Laboratory of Structural Biology of Zhejiang Province, School of Life Sciences, Westlake University, Hangzhou, Zhejiang 310024, China; ● orcid.org/0000-0001-9639-2907

Yandong Huang – College of Computer Engineering, Jimei University, Xiamen 361021, China; ● orcid.org/0000-0002-1452-6383

Phillip S. Hudson – Department of Chemistry, University of South Florida, Tampa, Florida 33620, United States; Medicine Design, Pfizer Inc., Cambridge, Massachusetts 02139, United States

Wonpil Im – Department of Biological Sciences, Lehigh University, Bethlehem, Pennsylvania 18015, United States; ● orcid.org/0000-0001-5642-6041

Shahidul M. Islam – Department of Chemistry, Delaware State University, Dover, Delaware 19901, United States; ● orcid.org/0000-0001-5769-6844

Wei Jiang – Computational Science Division, Argonne National Laboratory, Argonne, Illinois 60439, United States; ● orcid.org/0000-0001-7707-6247

Michael R. Jones – Laboratory of Computational Biology, National Heart Lung and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892, United States; ● orcid.org/0000-0002-0542-3796

Silvan Käser – Department of Chemistry, University of Basel, CH-4056 Basel, Switzerland; ● orcid.org/0000-0002-3641-8519

Fiona L. Kearns – Department of Chemistry, University of South Florida, Tampa, Florida 33620, United States; ● orcid.org/0000-0002-5469-9035

Nathan R. Kern – Department of Biological Sciences, Lehigh University, Bethlehem, Pennsylvania 18015, United States; ● orcid.org/0000-0003-0936-7274

Jeffery B. Klauda – Department of Chemical and Biomolecular Engineering, Institute for Physical Science and
Technology, Biophysics Program, University of Maryland, College Park, Maryland 20742, United States; ● orcid.org/0000-0001-8725-1870

**Themis Lazaridis** — Department of Chemistry, City College of New York, New York, New York 10031, United States; ● orcid.org/0000-0003-4218-7590

**Jinhyuk Lee** — Disease Target Structure Research Center, Korea Research Institute of Bioscience and Biotechnology, Daejeon 34141, Republic of Korea; Department of Bioinformatics, KRIBB School of Bioscience, University of Science and Technology, Daejeon 34141, Republic of Korea

**Justin A. Lemkul** — Department of Biochemistry, Virginia Polytechnic Institute and State University, Blacksburg, Virginia 24061, United States; ● orcid.org/0000-0001-6661-8653

**Xiaorong Liu** — Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109, United States; ● orcid.org/0000-0002-7687-3604

**Yun Luo** — Department of Biotechnology and Pharmaceutical Sciences, College of Pharmacy, Western University of Health Sciences, Pomona, California 91766, United States; ● orcid.org/0000-0003-3581-754X

**Alexander D. MacKerell, Jr.** — Department of Pharmaceutical Sciences, University of Maryland School of Pharmacy, Baltimore, Maryland 21201, United States; ● orcid.org/0000-0001-8287-6804

**Dan T. Major** — Department of Chemistry and Institute for Nanotechnology & Advanced Materials, Bar-Ilan University, Ramat-Gan 52900, Israel; ● orcid.org/0000-0002-9231-0676

**Markus Meuwly** — Department of Chemistry, University of Basel, CH-4056 Basel, Switzerland; Department of Chemistry, Brown University, Providence, Rhode Island 02912, United States; ● orcid.org/0000-0001-7930-8806

**Kwangho Nam** — Department of Chemistry and Biochemistry, University of Texas at Arlington, Arlington, Texas 76019, United States; ● orcid.org/0000-0003-0723-7839

**Lennart Nilsson** — Karolinska Institutet, Department of Biosciences and Nutrition, SE-14183 Huddinge, Sweden; ● orcid.org/0000-0002-5067-6397

**Victor Ovchinnikov** — Harvard University, Department of Chemistry and Chemical Biology, Cambridge, Massachusetts 02138, United States; ● orcid.org/0000-0002-1793-2352

**Emanuele Paci** — Dipartimento di Fisica e Astronomia, Università di Bologna, Bologna 40127, Italy; ● orcid.org/0000-0002-4891-2768

**Soohyung Park** — Department of Biological Sciences, Lehigh University, Bethlehem, Pennsylvania 18015, United States; ● orcid.org/0000-0002-4883-3031

**Richard W. Pastor** — Laboratory of Computational Biology, National Heart Lung and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892, United States; ● orcid.org/0000-0002-2454-5131

**Amanda R. Pittman** — Department of Chemistry, University of South Florida, Tampa, Florida 33620, United States

**Carol Beth Post** — Borch Department of Medicinal Chemistry and Molecular Pharmacology, Purdue University, West Lafayette, Indiana 47907, United States; ● orcid.org/0000-0002-5233-899X

**Samarjeet Prasad** — Laboratory of Computational Biology, National Heart Lung and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892, United States

**Jingzhi Pu** — Department of Chemistry and Chemical Biology, Indiana University Indianapolis, Indianapolis, Indiana 46202, United States; ● orcid.org/0000-0002-3042-335X

**Yifei Qi** — School of Pharmacy, Fudan University, Shanghai 201203, China; ● orcid.org/0000-0003-2853-7910

**Thenmalarchelvi Rathinavelan** — Department of Biotechnology, Indian Institute of Technology Hyderabad, Kandi, Telangana State 502284, India; ● orcid.org/0000-0002-1142-0583

**Daniel R. Roe** — Laboratory of Computational Biology, National Heart Lung and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892, United States; ● orcid.org/0000-0002-5834-2447

**Benoit Roux** — Department of Chemistry, University of Chicago, Chicago, Illinois 60637, United States; ● orcid.org/0000-0002-5254-2712

**Christopher N. Rowley** — Department of Chemistry, Carleton University, Ottawa, Ontario K1S 5B6, Canada; ● orcid.org/0000-0002-0205-952X

**Jana Shen** — Department of Pharmaceutical Sciences, University of Maryland School of Pharmacy, Baltimore, Maryland 21201, United States; ● orcid.org/0000-0002-3234-0769

**Andrew C. Simmonett** — Laboratory of Computational Biology, National Heart Lung and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892, United States; ● orcid.org/0000-0002-5921-9272

**Alexander J. Sodt** — Eunice Kennedy Shriver National Institute of Child Health and Human Development, National Institutes of Health, Bethesda, Maryland 20892, United States; ● orcid.org/0000-0002-5570-8212

**Kai Töpfer** — Department of Chemistry, University of Basel, CH-4056 Basel, Switzerland; ● orcid.org/0000-0002-4650-9641

**Meenu Upadhyay** — Department of Chemistry, University of Basel, CH-4056 Basel, Switzerland; ● orcid.org/0000-0003-3802-0835

**Arjan van der Vaart** — Department of Chemistry, University of South Florida, Tampa, Florida 33620, United States; ● orcid.org/0000-0002-8950-1850

**Luis Itza Vazquez-Salazar** — Department of Chemistry, University of Basel, CH-4056 Basel, Switzerland; ● orcid.org/0000-0001-6347-5108

**Richard M. Venable** — Laboratory of Computational Biology, National Heart Lung and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892, United States

**Luke C. Warrensford** — Department of Chemistry, University of South Florida, Tampa, Florida 33620, United States; ● orcid.org/0000-0002-9206-5121

**H. Lee Woodcock** — Department of Chemistry, University of South Florida, Tampa, Florida 33620, United States; ● orcid.org/0000-0003-3539-273X

**Yujin Wu** — Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109, United States

**Martin Karplus** — Harvard University, Department of Chemistry and Chemical Biology, Cambridge, Massachusetts 02138, United States; Laboratoire de Chimie Biophysique, ISIS, Université de Strasbourg, 67000 Strasbourg, France

Complete contact information is available at:
https://pubs.acs.org/10.1021/acs.jpcb.4c04100
## Author Contributions

M. Karplus provided the overall guidance of the CHARMM Development Project. W. Hwang assembled the manuscript; W. Hwang, C.L. Brooks, III, and B.R. Brooks edited the manuscript; J.E. Gonzales, II and W. Hwang edited the references. Sections ¹ and ¹³: W. Hwang, C.L. Brooks, III, and B.R. Brooks. Sections ²,¹,²,²,³,³,³,³,²,⁴,¹, and ⁵,¹: R.L. Hayes, J. Buckner, X. Liu, Y. Wu, X. Ding, M.F. Crowley, and C.L. Brooks, III. Sections ²,⁴ and ¹²,⁴: S. Prasad and B.R. Brooks. Section ³,³: Y. Qi, X. Cheng, Y. Choi, N.R. Kern, and W. Im. Sections ⁴,²,⁸,⁴,¹¹,¹¹,¹¹,²,¹¹,³,¹¹,⁴, and ¹¹,⁸: K. Nam, J. Pu, D.T. Major, and J. Gao. Section ⁴,³: L.C. Warrensford, A.R. Pittman, F.L. Kearns, and H.L. Woodcock. Section ⁵,²: B. Roux, W. Jiang, and Y. Deng. Section ⁵,³: A. Blondel. Section ⁶,¹: J. Shen and Y. Huang. Section ⁶,²: A. Damjanovic and B.R. Brooks. Sections ⁶,³,⁸,⁷,¹ and ⁸,⁷,²: T. Lazaridis. Section ⁷,¹: D.R. Roe and B.R. Brooks. Section ⁷,²: A. van der Vaart. Sections ⁷,³ (introduction), ⁷,³,¹ ⁷,³,² ⁷,³,³, and ⁷,³,⁴: V. Ovchinnikov and M. Karplus. Sections ⁷,³,⁵ and ⁷,³,⁶: B. Roux, Y. Luo, and E.D. Harder. Sections ⁷,⁴ and ¹²,⁹: C. Post. Section ⁷,⁵: H.T. Chang and J.W. Chu. Sections ⁷,⁶ and ¹¹,⁹: E. Paci and D.R. Glowacki. Section ⁷,⁷: S.L. Austin and H.L. Woodcock. Section ⁸,¹: M. Devereux, E.D. Boittier, K. Töpfer, and M. Meuwly. Section ⁸,²: K. Töpfer, S. Käser, L.I. Vazquez-Salazar, and M. Meuwly. Section ⁸,³: J. Huang, A. MacKerell, Jr., and B.R. Brooks. Section ⁸,⁵: A.C. Simmonett, B.R. Brooks, and R.W. Pastor. Section ⁸,⁶: A. Caflisch. Sections ⁸,⁷,³ and ⁸,⁸: M. Feig. Sections ⁹,¹ and ⁹,²: W. Im and J. Lee. Section ⁹,³: W. Im and T. Rathinavelan. Sections ⁹,⁴ and ¹²,¹⁰: L. Nilsson. Section ¹⁰: A.D. MacKerell, Jr., J. Huang, J.B. Klauda, J.A. Lemkul, L. Nilsson, R.W. Pastor, and B. Roux. Section ¹¹,⁵: Q. Cui. Section ¹¹,⁶: M.R. Jones and B.R. Brooks. Section ¹¹,⁷: B. Roux and C.N. Rowley. Section ¹¹,¹⁰: M. Upadhyay, and M. Meuwly. Section ¹¹,¹¹: S. Boresch, H.L. Woodcock, P.S. Hudson, and F.L. Kearns. Section ¹²,¹¹: B. Roux and S.M. Islam. Section ¹²,²: D.R. Roe, R.M. Venable, R.W. Pastor, and B.R. Brooks. Section ¹²,³: S. Park, A.J. Sodt, W. Im, and R.W. Pastor. Section ¹²,⁵: M. Buck and M.B. Hamaneh. Sections ¹²,⁶,¹²,⁷, and ¹²,⁸: W. Hwang and J.E. Gonzales, II.

## Notes

The authors declare the following competing financial interest(s): Alexander D. MacKerell, Jr. is Co-founder and CSO of SilcsBio LLC.

## ACKNOWLEDGMENTS

This work was supported by various funding agencies including the U.S. National Institutes of Health (NIH) and the U.S. National Science Foundation (NSF). The following separate acknowledgments are in response to requests from the authors: BR Brooks and RW Pastor: Intramural Research Program of the NIH, National Heart, Lung, and Blood Institute (BR Brooks: ZIA HL001051 and ZIA HL001052; RW Pastor: ZIA HL000340). CL Brooks: NIH R35GM137058. A Caflisch: Swiss National Science Foundation (grant number 310030-212195). JW Chu: National Science and Technology Council of Taiwan (112-2113-M-A49-018-MY3). Q Cui: NIH R35GM141930. M Feig: NIH R35GM126948. J Gao: NIH R01GM046736. D Glowacki: European Research Council Consolidator Grant NANOVR 866559. MB Hamaneh: Intramural program of the National Library of Medicine at the NIH. W Im: NSF MCB-211172 and NIH R01GM138472. W Jiang: The Office of Science, U.S. Department of Energy

(DE-AC02-06CH11357). T Lazaridis: NSF MCB-1855942. JA Lemkul: R35GM133754. AD MacKerell, Jr.: NIH R35GM131710. M Meuwly: The Swiss National Science Foundation (grants 215088 and 219779), the NCCR MUST, and the University of Basel. K Nam and DT Major: NIH R01GM132481 (K Nam) and NIH R21GM148895 (K Nam and DT Major). CB Post: NIH R01GM039478. J Pu: NIH R01GM135392. J Shen: NIH R35GM148261. A van der Vaart: NSF MCB-1919096. HL Woodcock: NIH R01GM129519, R01GM143810, and R01GM140316.

## REFERENCES

(1) McCammon, J. A.; Gelin, B. R.; Karplus, M. Dynamics of folded proteins. *Nature* **1977**, *267*, 585–90.

(2) Macuglia, D.; Roux, B.; Ciccotti, G. The emergence of protein dynamics simulations: how computational statistical mechanics met biochemistry. *Euro. Phys. J. H* **2022**, *47*, 13.

(3) Brooks, B. R.; Brooks, C. L., III; MacKerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: the biomolecular simulation program. *J. Comput. Chem.* **2009**, *30*, 1545–1614.

(4) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM - a program for macro-molecular energy, minimization, and dynamics calculations. *J. Comput. Chem.* **1983**, *4*, 187–217.

(5) Harris, J. A.; Liu, R.; de Oliveira, V. M.; Vázquez-Montelongo, E. A.; Henderson, J. A.; Shen, J. GPU-accelerated all-atom particle-mesh Ewald continuous constant pH molecular dynamics in Amber. *J. Chem. Theory Comp.* **2022**, *18*, 7510–7527.

(6) Aho, N.; Buslaev, P.; Jansen, A.; Bauer, P.; Groenhof, G.; Hess, B. Scalable constant pH molecular dynamics in GROMACS. *J. Chem. Theory Comp.* **2022**, *18*, 6148–6160.

(7) Jiang, W.; Chipot, C.; Roux, B. Computing relative binding affinity of ligands to receptor: an effective hybrid single-dual-topology free-energy perturbation approach in NAMD. *J. Chem. Inf. and Mod.* **2019**, *59*, 3794–3802.

(8) Eastman, P.; Swails, J.; Chodera, J. D.; McGibbon, R. T.; Zhao, Y.; Beauchamp, K. A.; Wang, L. P.; Simmonett, A. C.; Harrigan, M. P.; Stern, C. D.; et al. OpenMM 7: rapid development of high performance algorithms for molecular dynamics. *PLoS Comput. Biol.* **2017**, *13*, No. e1005659.

(9) Buckner, J.; Liu, X.; Chakravorty, A.; Wu, Y.; Cervantes, L. F.; Lai, T. T.; Brooks, C. L., III. pyCHARMM: embedding CHARMM functionality in a python framework. *J. Chem. Theory Comp.* **2023**, *19*, 3752–3762.

(10) Chakravorty, A.; Hussain, A.; Cervantes, L. F.; Lai, T.; Brooks, C. L., III. Exploring the limits of small molecule force fields through predictions of solvation free energy of small molecules. *J. Chem. Inf. and Mod.* **2024**, *64*, 4089–4101.

(11) Ding, X.; Vilseck, J. Z.; Brooks, C. L., III. Fast solver for large scale multistate Bennett acceptance ratio equations. *J. Chem. Theory Comp.* **2019**, *15*, 799–802.

(12) Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. *J. Chem. Phys.* **2008**, *129*, 124105.

(13) Ding, X.; Wu, Y.; Wang, Y.; Vilseck, J. Z.; Brooks, C. L., III. Accelerated CDOCKER with GPUs, parallel simulated annealing, and fast Fourier transforms. *J. Chem. Theory Comp.* **2020**, *16*, 3910–3919.

(14) Wu, Y.; Brooks, C. L., III. Flexible CDOCKER: hybrid searching algorithm and scoring function with side chain conformational entropy. *J. Chem. Inf. and Mod.* **2021**, *61*, 5535–5549.

(15) Wu, Y.; Brooks, C. L., III. Covalent docking in CDOCKER. *J. Comput. Aided. Mol. Des.* **2022**, *36*, 563–574.

(16) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-atom empirical potential for molecular modeling and dynamics studies of proteins. *J. Phys. Chem. B* **1998**, *102*, 3586–616.
(17) Vanommeslaeghe, K.; MacKerell, A. D., Jr. Automation of the CHARMM general force field (CGenFF) I: bond perception and atom typing. *J. Chem. Inf. and Mod.* **2012**, *52*, 3144–3154.

(18) Vanommeslaeghe, K.; Raman, E. P.; MacKerell, A. D., Jr. Automation of the CHARMM general force field (CGenFF) II: assignment of bonded parameters and partial atomic charges. *J. Chem. Inf. and Mod.* **2012**, *52*, 3155–3168.

(19) Tian, C.; Kasavajhala, K.; Belfon, K. A. A.; Raguette, L.; Huang, H.; Migues, A. N.; Bickel, J.; Wang, Y.; Pincay, J.; Wu, Q.; et al. ff19SB: amino-acid-specific protein backbone parameters trained against quantum mechanics energy surfaces in solution. *J. Chem. Theory Comp.* **2020**, *16*, 528–552.

(20) He, X.; Man, V. H.; Yang, W.; Lee, T.-S.; Wang, J. A fast and high-quality charge model for the next generation general AMBER force field. *J. Chem. Phys.* **2020**, *153*, 114502.

(21) Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. *J. Am. Chem. Soc.* **1996**, *118*, 11225–11236.

(22) Dodda, L. S.; Cabeza de Vaca, I.; Tirado-Rives, J.; Jorgensen, W. L. LigParGen web server: an automatic OPLS-AA parameter generator for organic ligands. *Nucleic Acids Res.* **2017**, *45*, W331–W336.

(23) Lemkul, J. A.; Huang, J.; Roux, B.; MacKerell, A. D., Jr. An empirical polarizable force field based on the classical Drude oscillator model: development history and recent applications. *Chem. Rev.* **2016**, *116*, 4983–5013.

(24) Hayes, R. L.; Buckner, J.; Brooks, C. L., III. BLaDE: a basic lambda dynamics engine for GPU accelerated molecular dynamics free energy calculations. *J. Chem. Theory Comp.* **2021**, *17*, 6799–6807.

(25) Hynninen, A. P.; Crowley, M. F. New faster CHARMM molecular dynamics engine. *J. Comput. Chem.* **2014**, *35*, 406–413.

(26) Knight, J. L.; Brooks, C. L., III. Multisite $\lambda$ dynamics for simulated structure-activity relationship studies. *J. Chem. Theory Comp.* **2011**, *7*, 2728–2739.

(27) Vilseck, J. Z.; Armacost, K. A.; Hayes, R. L.; Goh, G. B.; Brooks, C. L., III. Predicting binding free energies in a large combinatorial chemical space using multisite $\lambda$ dynamics. *J. Phys. Chem. Lett.* **2018**, *9*, 3328–3332.

(28) Raman, E. P.; Paul, T. J.; Hayes, R. L.; Brooks, C. L., III. Automated, accurate, and scalable relative protein-ligand binding free-energy calculations using lambda dynamics. *J. Chem. Theory Comp.* **2020**, *16*, 7895–7914.

(29) Goh, G. B.; Knight, J. L.; Brooks, C. L., III. Constant pH molecular dynamics simulations of nucleic acids in explicit solvent. *J. Chem. Theory Comp.* **2012**, *8*, 36–46.

(30) Goh, G. B.; Knight, J. L.; Brooks, C. L., III. Towards accurate prediction of protonation equilibrium of nucleic acids. *J. Phys. Chem. Lett.* **2013**, *4*, 760–766.

(31) Armacost, K. A.; Goh, G. B.; Brooks, C. L., III. Biasing potential replica exchange multisite $\lambda$-dynamics for efficient free energy calculations. *J. Chem. Theory Comp.* **2015**, *11*, 1267–1277.

(32) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant pressure molecular dynamics simulation: the Langevin piston method. *J. Chem. Phys.* **1995**, *103*, 4613–4621.

(33) Chow, K.-H.; Ferguson, D. M. Isothermal-isobaric molecular dynamics simulations with Monte Carlo volume sampling. *Comput. Phys. Commun.* **1995**, *91*, 283–289.

(34) Åqvist, J.; Wennerström, P.; Nervall, M.; Bjelic, S.; Brandsdal, B. O. Molecular dynamics simulations of water and biomolecules with a Monte Carlo constant pressure algorithm. *Chem. Phys. Lett.* **2004**, *384*, 288–294.

(35) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald - an NLog(N) method for Ewald sums in large systems. *J. Chem. Phys.* **1993**, *98*, 10089–10092.

(36) Steinbach, P. J.; Brooks, B. R. New spherical-cutoff methods for long-range forces in macromolecular simulation. *J. Comput. Chem.* **1994**, *15*, 667–683.

(37) Christ, C. D.; van Gunsteren, W. F. Enveloping distribution sampling: a method to calculate free energy differences from a single simulation. *J. Chem. Phys.* **2007**, *126*, 184110.

(38) Christ, C. D.; van Gunsteren, W. F. Multiple free energies from a single simulation: extending enveloping distribution sampling to nonoverlapping phase-space distributions. *J. Chem. Phys.* **2008**, *128*, 174112.

(39) Christ, C. D.; Van, W. F. Comparison of three enveloping distribution sampling Hamiltonians for the estimation of multiple free energy differences from a single simulation. *J. Comput. Chem.* **2009**, *30*, 1664–1679.

(40) Lee, J.; Miller, B. T.; Damjanović, A.; Brooks, B. R. Constant pH molecular dynamics in explicit solvent with enveloping distribution sampling and Hamiltonian exchange. *J. Chem. Theory Comp.* **2014**, *10*, 2738–2750.

(41) Lee, J.; Miller, B. T.; Damjanović, A.; Brooks, B. R. Enhancing constant-pH simulation in explicit solvent with a two-dimensional replica exchange method. *J. Chem. Theory Comp.* **2015**, *11*, 2560–2574.

(42) Dolan, E. A.; Venable, R. M.; Pastor, R. W.; Brooks, B. R. Simulations of membranes and other interfacial systems using P21 and Pc periodic boundary conditions. *Biophys. J.* **2002**, *82*, 2317–2325.

(43) Prasad, S.; Simmonett, A. C.; Meana-Paneda, R.; Brooks, B. R. The extended eighth-shell method for periodic boundary conditions with rotational symmetry. *J. Comput. Chem.* **2021**, *42*, 1373–1383.

(44) Park, S.; Rice, A.; Im, W.; Pastor, R. W. Spontaneous curvature generation by peptides in asymmetric bilayers. *J. Comput. Chem.* **2023**, *45*, 512–522.

(45) Mey, A. S. J. S.; Allen, B. K.; Bruce McDonald, H. E.; Chodera, J. D.; Hahn, D. F.; Kuhn, M.; Michel, J.; Mobley, D. L.; Naden, L. N.; Prasad, S.; et al. Best practices for alchemical free energy calculations. *Living Journal of Computational Molecular Science* **2020**, *2*, 18378.

(46) Zacharias, M.; Straatsma, T. P.; McCammon, J. A. Separation-shifted scaling, a new scaling method for Lennard-Jones interactions in thermodynamic integration. *J. Chem. Phys.* **1994**, *100*, 9025–9031.

(47) Wu, X.; Brooks, B. R. A double exponential potential for van der Waals interaction. *AIP Adv.* **2019**, *9*, 065304.

(48) Hynninen, A.-P. ap-hynninen/GSE. 2017.

(49) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: a web-based graphical user interface for CHARMM. *J. Comput. Chem.* **2008**, *29*, 1859–1865.

(50) Martins de Oliveira, V.; Liu, R.; Shen, J. Constant pH molecular dynamics simulations: current status and recent applications. *Curr. Op. Struct. Biol.* **2022**, *77*, 102498.

(51) Brunger, A. T.; Campbell, R. L.; Clore, G. M.; Gronenborn, A. M.; Karplus, M.; Petsko, G. A.; Teeter, M. M. Solution of a protein crystal structure with a model obtained from NMR interproton distance restraints. *Science* **1987**, *235*, 1049–1053.

(52) Raman, E. P.; Yu, W.; Guvench, O.; Mackerell, A. D., Jr. Reproducing crystal binding modes of ligand functional groups using site-identification by ligand competitive saturation (SILCS) simulations. *J. Chem. Inf. and Mod.* **2011**, *51*, 877–896.

(53) Gao, X.; Ramezanghorbani, F.; Isayev, O.; Smith, J. S.; Roitberg, A. E. TorchANI: a free and open source PyTorch-based deep learning implementation of the ANI neural network potentials. *J. Chem. Inf. and Mod.* **2020**, *60*, 3408–3415.

(54) Unke, O. T.; Meuwly, M. PhysNet: a neural network for predicting energies, forces, dipole moments, and partial charges. *J. Chem. Theory Comp.* **2019**, *15*, 3678–3693.

(55) Wu, G.; Robertson, D. H.; Brooks, C. L., III; Vieth, M. Detailed analysis of grid-based molecular docking: a case study of CDOCKER-a CHARMM-based Md docking algorithm. *J. Comput. Chem.* **2003**, *24*, 1549–1562.

(56) Gagnon, J. K.; Law, S. M.; Brooks, C. L., III. Flexible CDOCKER: development and application of a pseudo-explicit structure-based docking method within CHARMM. *J. Comput. Chem.* **2016**, *37*, 753–762.
(57) Buckner, J.; Liu, X.; Chakravorty, A.; Wu, Y.; Cervantes, L. F.; Lai, T. T.; Brooks, C. L., III. pyCHARMM-Workshop (accessed August 5, 2024); https://github.com/BrooksResearchGroup-UM/pyCHARMM-Workshop.

(58) Hayes, L.; Ryan, Vilseck, J. Z., Buckner, J., Liu, X., Chakravorty, A., Cervantes, L. F., Lai, T. T., Brooks, C. L., III. MSLD-Workshop 2024, https://github.com/BrooksResearchGroup-UM/MSLD-Workshop.

(59) Jo, S.; Kim, T.; Im, W. Automated builder and database of protein/membrane complexes for molecular dynamics simulations. *PLoS One* **2007**, 2, No. e880.

(60) Jo, S.; Cheng, X.; Lee, J.; Kim, S.; Park, S. J.; Patel, D. S.; Beaven, A. H.; Lee, K. I.; Rui, H.; Park, S.; et al. CHARMM-GUI 10 years for biomolecular modeling and simulation. *J. Comput. Chem.* **2017**, 38, 1114–1124.

(61) Jo, S.; Lim, J. B.; Klauda, J. B.; Im, W. CHARMM-GUI membrane builder for mixed bilayers and its application to yeast membranes. *Biophys. J.* **2009**, *97*, 50–58.

(62) Wu, E. L.; Cheng, X.; Jo, S.; Rui, H.; Song, K. C.; Davila-Contreras, E. M.; Qi, Y.; Lee, J.; Monje-Galvan, V.; Venable, R. M.; et al. CHARMM-GUI membrane builder toward realistic biological membrane simulations. *J. Comput. Chem.* **2014**, *35*, 1997–2004.

(63) Kunzmann, P.; Muller, T. D.; Greil, M.; Krumbach, J. H.; Anter, J. M.; Bauer, D.; Islam, F.; Hamacher, K. Biotite: new tools for a versatile Python bioinformatics library. *BMC Bioinf* **2023**, *24*, 236.

(64) Cock, P. J.; Antao, T.; Chang, J. T.; Chapman, B. A.; Cox, C. J.; Dalke, A.; Friedberg, I.; Hamelryck, T.; Kauff, F.; Wilczynski, B.; et al. Biopython: freely available Python tools for computational molecular biology and bioinformatics. *Bioinf* **2009**, *25*, 1422–3.

(65) Xu, L., Zigqiao; Brooks, C. L., III. crimm (accessed August 5, 2024); https://github.com/BrooksResearchGroup-UM/crimm.

(66) Rose, A. S.; Bradley, A. R.; Valasatava, Y.; Duarte, J. M.; Prlic, A.; Rose, P. W. NGL viewer: web-based molecular graphics for large complexes. *Bioinf* **2018**, *34*, 3755–3758.

(67) Burley, S. K.; Bhikadiya, C.; Bi, C.; Bittrich, S.; Chao, H.; Chen, L.; Craig, P. A.; Crichlow, G. V.; Dalenberg, K.; Duarte, J. M.; et al. RCSB Protein Data Bank (RCSB.org): delivery of experimentally-determined PDB structures alongside one million computed structure models of proteins from artificial intelligence/machine learning. *Nucleic Acids Res.* **2023**, *51*, D488–D508.

(68) Tunyasuvunakool, K.; Adler, J.; Wu, Z.; Green, T.; Zielinski, M.; Zidek, A.; Bridgland, A.; Cowie, A.; Meyer, C.; Laydon, A.; et al. Highly accurate protein structure prediction for the human proteome. *Nature* **2021**, *596*, 590–596.

(69) Westbrook, J. D.; Young, J. Y.; Shao, C.; Feng, Z.; Guranovic, V.; Lawson, C. L.; Vallat, B.; Adams, P. D.; Berrisford, J. M.; Bricogne, G.; et al. PDBx/mmCIF ecosystem: foundational semantic tools for structural biology. *J. Mol. Biol.* **2022**, *434*, 167599.

(70) Sondergaard, C. R.; Olsson, M. H.; Rostkowski, M.; Jensen, J. H. Improved treatment of ligands and coupling effects in empirical calculation and rationalization of pKa values. *J. Chem. Theory Comp.* **2011**, *7*, 2284–95.

(71) Landrum, G.; Tosco, P.; Kelley, B.; Rodriguez, R.; Cosgrove, D.; Vianello, R.; Schneider, N.; Kawashima, E.; Jones, G.; Nealschneider, D.; et al. RDKit (accessed August 5, 2024);

(72) Boothroyd, S.; Behara, P. K.; Madin, O. C.; Hahn, D. F.; Jang, H.; Gapsys, V.; Wagner, J. R.; Horton, J. T.; Dotson, D. L.; Thompson, M. W.; et al. Development and benchmarking of open force field 2.0.0: the sage small molecule force field. *J. Chem. Theory Comp.* **2023**, *19*, 3251–3275.

(73) Eberhardt, J.; Santos-Martins, D.; Tillack, A. F.; Forli, S. AutoDock Vina 1.2.0: new docking methods, expanded force field, and python bindings. *J. Chem. Inf. and Mod.* **2021**, *61*, 3891–3898.

(74) Case, D. A.; Cheatham III, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M., J.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber biomolecular simulation programs. *J. Comput. Chem.* **2005**, *26*, 1668–1688.

(75) Bowers, K. J.; Chow, E.; Xu, H.; Dror, R. O.; Eastwood, M. P.; Gregersen, B. A.; Klepeis, J. L.; Kolossvary, I.; Moraes, M. A;

Sacerdoti, F. D.; et al. Scalable algorithms for molecular dynamics simulations on commodity clusters. *Proceedings of the 2006 ACM/ IEEE Conference on Supercomputing. 2006; p. 84.

(76) Jung, J.; Mori, T.; Kobayashi, C.; Matsunaga, Y.; Yoda, T.; Feig, M.; Sugita, Y. GENESIS: a hybrid-parallel and multi-scale molecular dynamics simulator with enhanced sampling algorithms for biomolecular and cellular simulations. *Wiley Interdiscip. Rev. Comput. Mol. Sci.* **2015**, *5*, 310–323.

(77) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: high performance molecular simulations through multi-level parallelism from laptops to super-computers. *SoftwareX* **2015**, *1–2*, 19–25.

(78) Thompson, A. P.; Aktulga, H. M.; Berger, R.; Bolintineanu, D. S.; Brown, W. M.; Crozier, P. S.; in 't Veld, P. J.; Kohlmeyer, A.; Moore, S. G.; Nguyen, T. D.; et al. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. *Comput. Phys. Commun.* **2022**, *271*, 108171.

(79) Phillips, J. C.; Hardy, D. J.; Maia, J. D. C.; Stone, J. E.; Ribeiro, J. V.; Bernardi, R. C.; Buch, R.; Fiorin, G.; Henin, J.; Jiang, W.; et al. Scalable molecular dynamics on CPU and GPU architectures with NAMD. *J. Chem. Phys.* **2020**, *153*, 044130.

(80) Rackers, J. A.; Wang, Z.; Lu, C.; Laury, M. L.; Lagardere, L.; Schnieders, M. J.; Piquemal, J. P.; Ren, P.; Ponder, J. W. Tinker 8: software tools for molecular design. *J. Chem. Theory Comp.* **2018**, *14*, 5273–5289.

(81) Woo, H.; Park, S. J.; Choi, Y. K.; Park, T.; Tanveer, M.; Cao, Y.; Kern, N. R.; Lee, J.; Yeom, M. S.; Croll, T. I.; et al. Developing a fully glycosylated full-length SARS-CoV-2 spike protein model in a viral membrane. *J. Phys. Chem. B* **2020**, *124*, 7128–7137.

(82) Pogozheva, I. D.; Armstrong, G. A.; Kong, L.; Hartnagel, T. J.; Carpino, C. A.; Gee, S. E.; Picarello, D. M.; Rubin, A. S.; Lee, J.; Park, S.; et al. Comparative molecular dynamics simulation studies of realistic eukaryotic, prokaryotic, and archaeal membranes. *J. Chem. Inf. and Mod.* **2022**, *62*, 1036–1051.

(83) Pogozheva, I. D.; Cherepanov, S.; Park, S. J.; Raghavan, M.; Im, W.; Lomize, A. L. Structural modeling of cytokine-receptor-JAK2 signaling complexes using AlphaFold multimer. *J. Chem. Inf. and Mod.* **2023**, *63*, 5874–5895.

(84) Cao, Y.; Choi, Y. K.; Frank, M.; Woo, H.; Park, S. J.; Yeom, M. S.; Seok, C.; Im, W. Dynamic interactions of fully glycosylated SARS-CoV-2 spike protein with various antibodies. *J. Chem. Theory Comp.* **2017**, *17*, 6559–6569.

(85) Choi, Y. K.; Cao, Y.; Frank, M.; Woo, H.; Park, S. J.; Yeom, M. S.; Croll, T. I.; Seok, C.; Im, W. Structure, dynamics, receptor binding, and antibody binding of the fully glycosylated full-length SARS-CoV-2 spike protein in a viral membrane. *J. Chem. Theory Comp.* **2021**, *17*, 2479–2487.

(86) Qi, Y.; Cheng, X.; Han, W.; Jo, S.; Schulten, K.; Im, W. CHARMM-GUI PACE CG Builder for solution, micelle, and bilayer coarse-grained simulations. *J. Chem. Inf. and Mod.* **2014**, *54*, 1003–1009.

(87) Lee, J.; Cheng, X.; Swails, J. M.; Yeom, M. S.; Eastman, P. K.; Lemkul, J. A.; Wei, S.; Buckner, J.; Jeong, J. C.; Qi, Y.; et al. CHARMM-GUI input generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 additive force field. *J. Chem. Theory Comp.* **2016**, *12*, 405–413.

(88) Lee, J.; Hitzenberger, M.; Rieger, M.; Kern, N. R.; Zacharias, M.; Im, W. CHARMM-GUI supports the Amber force fields. *J. Chem. Phys.* **2020**, *153*, 035103.

(89) Cheng, X.; Jo, S.; Lee, H. S.; Klauda, J. B.; Im, W. CHARMM-GUI micelle builder for pure/mixed micelle and protein/micelle complex systems. *J. Chem. Inf. and Mod.* **2013**, *53*, 2171–2180.

(90) Qi, Y.; Cheng, X.; Lee, J.; Vermaas, J. V.; Pogorelov, T. V.; Tajkhorshid, E.; Park, S.; Klauda, J. B.; Im, W. CHARMM-GUI HMMM builder for membrane simulations with the highly mobile membrane-mimetic model. *Biophys. J.* **2015**, *109*, 2012–2022.

(91) Lee, J.; Patel, D. S.; Stahle, J.; Park, S. J.; Kern, N. R.; Kim, S.; Lee, J.; Cheng, X.; Valvano, M. A.; Holst, O.; et al. CHARMM-GUI
membrane builder for complex biological membrane simulations with glycolipids and lipoglycans. *J. Chem. Theory Comp.* **2019**, *15*, 775–786.

(92) Qi, Y.; Lee, J.; Klauda, J. B.; Im, W. CHARMM-GUI nanodisc builder for modeling and simulation of various nanodisc systems. *J. Comput. Chem.* **2019**, *40*, 893–899.

(93) Park, S.; Choi, Y. K.; Kim, S.; Lee, J.; Im, W. CHARMM-GUI membrane builder for lipid nanoparticles with ionizable cationic lipids and PEGylated lipids. *J. Chem. Inf. and Mod.* **2021**, *61*, 5192–5202.

(94) Feng, S.; Park, S.; Choi, Y. K.; Im, W. CHARMM-GUI membrane builder: past, current, and future developments and applications. *J. Chem. Theory Comp.* **2023**, *19*, 2161–2185.

(95) Jo, S.; Cheng, X.; Islam, S. M.; Huang, L.; Rui, H.; Zhu, A.; Lee, H. S.; Qi, Y.; Han, W.; Vanommeslaeghe, K.; et al. CHARMM-GUI PDB manipulator for advanced modeling and simulations of proteins containing nonstandard residues. *Adv. Prot. Chem. Struct. Biol.* **2014**, *96*, 235–265.

(96) Park, S. J.; Kern, N.; Brown, T.; Lee, J.; Im, W. CHARMM-GUI PDB manipulator: various PDB structural modifications for biomolecular modeling and simulation. *J. Mol. Biol.* **2023**, *435*, 167995.

(97) Jo, S.; Jiang, W.; Lee, H. S.; Roux, B.; Im, W. CHARMM-GUI ligand binder for absolute binding free energy calculations and its application. *J. Chem. Inf. and Mod.* **2013**, *53*, 267–277.

(98) Kim, S.; Lee, J.; Jo, S.; Brooks, C. L., III; Lee, H. S.; Im, W. CHARMM-GUI ligand reader and modeler for CHARMM force field generation of small molecules. *J. Comput. Chem.* **2017**, *38*, 1879–1886.

(99) Kim, S.; Oshima, H.; Zhang, H.; Kern, N. R.; Re, S.; Lee, J.; Roux, B.; Sugita, Y.; Jiang, W.; Im, W. CHARMM-GUI free energy calculator for absolute and relative ligand solvation and binding free energy simulations. *J. Chem. Theory Comp.* **2020**, *16*, 7207–7218.

(100) Guterres, H.; Park, S. J.; Zhang, H.; Im, W. CHARMM-GUI LBS finder & refiner for ligand binding site prediction and refinement. *J. Chem. Inf. and Mod.* **2021**, *61*, 3744–3751.

(101) Guterres, H.; Park, S. J.; Cao, Y.; Im, W. CHARMM-GUI ligand designer for template-based virtual ligand design in a binding site. *J. Chem. Inf. and Mod.* **2021**, *61*, 5336–5342.

(102) Zhang, H.; Kim, S.; Giese, T. J.; Lee, T. S.; Lee, J.; York, D. M.; Im, W. CHARMM-GUI free energy calculator for practical ligand binding free energy simulations with AMBER. *J. Chem. Inf. and Mod.* **2021**, *61*, 4145–4151.

(103) Guterres, H.; Park, S.-J.; Zhang, H.; Perone, T.; Kim, J.; Im, W. CHARMM-GUI high-throughput simulator for efficient evaluation of protein-ligand interactions with different force fields. *Protein Sci.* **2022**, *31*, No. e4413.

(104) Guterres, H.; Im, W. CHARMM-GUI-based induced fit docking workflow to generate reliable protein-ligand binding modes. *J. Chem. Inf. and Mod.* **2023**, *63*, 4772–4779.

(105) Jo, S.; Song, K. C.; Desaire, H.; MacKerell, A. D., Jr.; Im, W. Glycan Reader: automated sugar identification and simulation preparation for carbohydrates and glycoproteins. *J. Comput. Chem.* **2011**, *32*, 3135–3141.

(106) Park, S. J.; Lee, J.; Patel, D. S.; Ma, H.; Lee, H. S.; Jo, S.; Im, W. Glycan reader is improved to recognize most sugar types and chemical modifications in the protein data bank. *Bioinf* **2017**, *33*, 3051–3057.

(107) Park, S. J.; Lee, J.; Qi, Y.; Kern, N. R.; Lee, H. S.; Jo, S.; Joung, L.; Joo, K.; Lee, J.; Im, W. CHARMM-GUI glycan modeler for modeling and simulation of carbohydrates and glycoconjugates. *Glycobiol* **2019**, *29*, 320–331.

(108) Jo, S.; Vargyas, M.; Vasko-Szedlar, J.; Roux, B.; Im, W. PBEQ-solver for online visualization of electrostatic potential of biomolecules. *Nucleic Acids Res.* **2008**, *36*, W270–W275.

(109) Qi, Y.; Ingolfsson, H. I.; Cheng, X.; Lee, J.; Marrink, S. J.; Im, W. CHARMM-GUI martini maker for coarse-grained simulations with the martini force field. *J. Chem. Theory Comp.* **2015**, *11*, 4486–4494.

(110) Hsu, P. C.; Bruininks, B. M. H.; Jefferies, D.; Cesar Telles de Souza, P.; Lee, J.; Patel, D. S.; Marrink, S. J.; Qi, Y.; Khalid, S.; Im, W. CHARMM-GUI martini maker for modeling and simulation of complex bacterial membranes with lipopolysaccharides. *J. Comput. Chem.* **2017**, *38*, 2354–2363.

(111) Qi, Y.; Lee, J.; Cheng, X.; Shen, R.; Islam, S. M.; Roux, B.; Im, W. CHARMM-GUI DEER facilitator for spin-pair distance distribution calculations and preparation of restrained-ensemble molecular dynamics simulations. *J. Comput. Chem.* **2020**, *41*, 415–420.

(112) Choi, Y. K.; Park, S. J.; Park, S.; Kim, S.; Kern, N. R.; Lee, J.; Im, W. CHARMM-GUI polymer builder for modeling and simulation of synthetic polymers. *J. Chem. Theory Comp.* **2021**, *17*, 2431–2443.

(113) Gao, Y.; Lee, J.; Smith, I. P. S.; Lee, H.; Kim, S.; Qi, Y.; Klauda, J. B.; Widmalm, G.; Khalid, S.; Im, W. CHARMM-GUI supports hydrogen mass repartitioning and different protonation states of phosphates in lipopolysaccharides. *J. Chem. Inf. and Mod.* **2021**, *61*, 831–839.

(114) Choi, Y. K.; Kern, N. R.; Kim, S.; Kanhaiya, K.; Afshar, Y.; Jeon, S. H.; Jo, S.; Brooks, B. R.; Lee, J.; Tadmor, E. B.; et al. CHARMM-GUI nanomaterial modeler for modeling and simulation of nanomaterial systems. *J. Chem. Theory Comp.* **2022**, *18*, 479–493.

(115) Kognole, A. A.; Lee, J.; Park, S. J.; Jo, S.; Chatterjee, P.; Lemkul, J. A.; Huang, J.; MacKerell, A. D., Jr.; Im, W. CHARMM-GUI Drude prepper for molecular dynamics simulation using the classical Drude polarizable force field. *J. Comput. Chem.* **2022**, *43*, 359–375.

(116) Zhang, H.; Kim, S.; Im, W. Practical guidance for consensus scoring and force field selection in protein-ligand binding free energy simulations. *J. Chem. Inf. and Mod.* **2022**, *62*, 6084–6093.

(117) Lee, K. I.; Jo, S.; Rui, H.; Egwolf, B.; Roux, B.; Pastor, R. W.; Im, W. Web interface for Brownian dynamics simulation of ion transport and its applications to beta-barrel pores. *J. Comput. Chem.* **2012**, *33*, 331–339.

(118) Qi, Y.; Lee, J.; Singharoy, A.; McGreevy, R.; Schulten, K.; Im, W. CHARMM-GUI MDFF/xMDFF Utilizer for molecular dynamics flexible fitting simulations in various environments. *J. Phys. Chem. B* **2017**, *121*, 3718–3723.

(119) Suh, D.; Feng, S.; Lee, H.; Zhang, H.; Park, S. J.; Kim, S.; Lee, J.; Choi, S.; Im, W. CHARMM-GUI enhanced sampler for various collective variables and enhanced sampling methods. *Protein Sci.* **2022**, *31*, No. e4446.

(120) Wang, K. W.; Lee, J.; Zhang, H.; Suh, D.; Im, W. CHARMM-GUI implicit solvent modeler for various generalized Born models in different simulation programs. *J. Phys. Chem. B* **2022**, *126*, 7354–7364.

(121) Kondo, D.; Taufer, M.; Brooks, C. L., III; Casanova, H.; Chien, A. A. Characterizing and evaluating desktop grids: An empirical study. *18th International Parallel and Distributed Processing Symposium*, 2004. Proceedings. 2004, p 26.

(122) Taufer, M.; Kerstens, A.; Estrada, T.; Flores, D. A.; Zamudio, R.; Teller, P. J.; Armen, R. S.; Brooks, C. L., III. Moving volunteer computing towards knowledge-constructed, dynamically-adaptive modeling and scheduling. *International Parallel and Distributed Processing Symposium* 2007, 1–8.

(123) Taufer, M.; An, C.; Kerstens, A.; Brooks, C. L., III. Predictor@ Home: a “protein structure prediction supercomputer” based on global computing. *IEEE Transac. Paral. Distrib. Sys.* **2006**, *17*, 786–796.

(124) Taufer, M.; Armen, R.; Chen, J.; Teller, P.; Brooks, C. L., III. Computational multiscale modeling in protein–ligand docking. *IEEE Eng. Med. Biol. Mag.* **2009**, *28*, 58–69.

(125) Rahaman, O.; Estrada, T. P.; Doren, D. J.; Taufer, M.; Brooks, C. L., III; Armen, R. S. Evaluation of several two-step scoring functions based on linear interaction energy, effective ligand size, and empirical pair potentials for prediction of protein-ligand binding geometry and free energy. *J. Chem. Inf. and Mod.* **2011**, *51*, 2047–2065.
(126) Armen, R. S.; Chen, J.; Brooks, C. L., III. An evaluation of explicit receptor flexibility in molecular docking using molecular dynamics and torsion angle molecular dynamics. *J. Chem. Theory Comp.* **2009**, *5*, 2909–2923.

(127) Ahmed, F.; Brooks, C. L., III. FASTDock: a pipeline for allosteric drug discovery. *J. Chem. Inf. and Mod.* **2023**, *63*, 7219–7227.

(128) Horton, J. T.; Boothroyd, S.; Wagner, J.; Mitchell, J. A.; Gokey, T.; Dotson, D. L.; Behara, P. K.; Ramaswamy, V. K.; Mackey, M.; Chodera, J. D.; et al. Open Force Field BespokeFit: automating bespoke torsion parametrization at scale. *J. Chem. Inf. and Mod.* **2022**, *62*, 5622–5633.

(129) Dominy, B. N.; Brooks, C. L., III. Development of a generalized Born model parametrization for proteins and nucleic acids. *J. Phys. Chem. B* **1999**, *103*, 3765–3773.

(130) Lee, M. S.; Salsbury, F. R., Jr.; Brooks, C. L., III. Novel generalized Born methods. *J. Chem. Phys.* **2002**, *116*, 10606–10614.

(131) Im, W.; Lee, M. S.; Brooks, C. L., III. Generalized Born model with a simple smoothing function. *J. Comput. Chem.* **2003**, *24*, 1691–1702.

(132) Haberthür, U.; Caflisch, A. FACTS: fast analytical continuum treatment of solvation. *J. Comput. Chem.* **2008**, *29*, 701–715.

(133) Ohashi, M.; Liu, F.; Hai, Y.; Chen, M.; Tang, M.; Yang, Z.; Sato, M.; Watanabe, K.; Houk, K. N.; Tang, Y. SAM-dependent enzyme-catalysed pericyclic reactions in natural product biosynthesis. *Nature* **2017**, *549*, 502–506.

(134) Chang, Z.; Ansbacher, T.; Zhang, L.; Yang, Y.; Ko, T.-P.; Zhang, G.; Liu, W.; Huang, J.-W.; Dai, L.; Guo, R.-T.; et al. Crystal structure of Lepi, a multifunctional SAM-dependent enzyme which catalyzes pericyclic reactions in leporin biosynthesis. *Org. Biomol. Chem.* **2019**, *17*, 2070–2076.

(135) Das, S.; Shimshi, M.; Raz, K.; Nitoker Eliaz, N.; Mhashal, A. R.; Ansbacher, T.; Major, D. T. EnzyDock: protein-ligand docking of multiple reactive states along a reaction coordinate in enzymes. *J. Chem. Theory Comp.* **2019**, *15*, 5116–5134.

(136) Vieth, M.; Hirst, J. D.; Dominy, B. N.; Daigler, H.; Brooks, C. L., III. Assessing search strategies for flexible docking. *J. Comput. Chem.* **1998**, *19*, 1623–1631.

(137) Vieth, M.; Hirst, J. D.; Kolinski, A.; Brooks, C. L., III. Assessing energy functions for flexible docking. *J. Comput. Chem.* **1998**, *19*, 1612–1622.

(138) Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: umbrella sampling. *J. Chem. Phys.* **1977**, *23*, 187–199.

(139) Ovchinnikov, V.; Karplus, M.; Vanden-Eijnden, E. Free energy of conformational transition paths in biomolecules: the string method and its application to myosin VI. *J. Chem. Phys.* **2011**, *134*, 085103.

(140) Bussi, G.; Laio, A. Using metadynamics to explore complex free-energy landscapes. *Nat. Rev. Phys.* **2020**, *2*, 200–212.

(141) Demapan, D.; Kussman, J.; Ochsenfeld, C.; Cui, Q. Factors that determine the variation of equilibrium and kinetic properties of QM/MM enzyme simulations: QM region, conformation and boundary condition. *J. Chem. Theory Comp.* **2022**, *18*, 2530–2542.

(142) Raz, K.; Driller, R.; Dimos, N.; Ringel, M.; Brück, T.; Loll, B.; Major, D. T. The impression of a nonexisting catalytic effect: the role of CotB2 in guiding the complex biosynthesis of Cyclooctat-9-en-7-ol. *J. Am. Chem. Soc.* **2020**, *142*, 21562–21574.

(143) Hu, J.; Ma, A.; Dinner, A. R. Monte Carlo simulations of biomolecules: the mc module in CHARMM. *J. Comput. Chem.* **2006**, *27*, 203–216.

(144) Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; de Groot, B. L.; Grubmüller, H.; MacKerell, A. D., Jr. CHARMM36: an improved force field for folded and intrinsically disordered proteins. *Nat. Meth.* **2017**, *14*, 71–73.

(145) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. *J. Comput. Chem.* **2010**, *31*, 671–690.

(146) Hudson, P. S.; Aviat, F.; Meana-Pañeda, R.; Warrensford, L.; Pollard, B. C.; Prasad, S.; Jones, M. R.; Woodcock III, H. L.; Brooks, B. R. Obtaining QM/MM binding free energies in the SAMPL8 drugs of abuse challenge: indirect approaches. *J. Comput. Aided Mol. Des.* **2022**, *36*, 263–277.

(147) Ojeda-May, P.; Nam, K. Acceleration of semiempirical QM/MM methods through message passing interface (MPI), hybrid MPI/Open Mmltiprocessing, and self-consistent field accelerator implementations. *J. Chem. Theory Comp.* **2017**, *13*, 3525–3536.

(148) Woodcock III, H. L.; Hodošček, M.; Gilbert, A. T. B.; Gill, P. M. W.; Schaefer, H. F., III; Brooks, B. R. Interfacing Q-Chem and CHARMM to perform QM/MM reaction path calculations. *J. Comput. Chem.* **2007**, *28*, 1485–1502.

(149) Gupta, P. K.; Tarannam, N.; Zev, S.; Major, D. T. Multistate multiscale docking study of the hydrolysis of toxic nerve agents by phosphotriesterase. *Electron. Struct.* **2023**, *5*, 035003.

(150) Raz, K.; Levi, S.; Gupta, P. K.; Major, D. T. Enzymatic control of product distribution in terpene syntheses: insights from multiscale simulations. *Curr. Opin. Biotechnol.* **2020**, *65*, 248–258.

(151) Zev, S.; Raz, K.; Schwartz, R.; Tarabeh, R.; Gupta, P. K.; Major, D. T. Benchmarking the ability of common docking programs to correctly reproduce and score binding modes in SARS-CoV-2 protease Mpro. *J. Chem. Inf. and Mod.* **2021**, *61*, 2957–2966.

(152) Koshland, D. Application of a theory of enzyme specificity to protein synthesis. *Proc. Natl. Acad. Sci. USA* **1958**, *44*, 98–104.

(153) Cavasotto, C. N.; Abagyan, R. A. Protein Flexibility in Ligand Docking and Virtual Screening to Protein Kinases. *J. Mol. Biol.* **2004**, *337*, 209–25.

(154) Sherman, W.; Day, T.; Jacobson, M. P.; Friesner, R. A.; Farid, R. Novel procedure for modeling ligand/receptor induced fit effects. *J. Med. Chem.* **2006**, *49*, 534–553.

(155) Vankayala, S.; Warrensford, L. C.; Pittman, A.; Pollard, B.; Kearns, F.; Larkin, J.; Woodcock, H. L. CIFDock: A novel CHARMM-based flexible receptor-flexible ligand docking protocol. *J. Comput. Chem.* **2022**, *43*, 84–95.

(156) O'Boyle, N. M.; Banck, M.; James, C. A.; Morely, C.; Vandermeersch, T.; Hutchison, G. R. Open Babel: An open chemical toolbox. *J. Cheminform.* **2011**, *3*, 33.

(157) Wu, X.; Brooks, B. R. Self-guided Langevin dynamics simulation method. *Chem. Phys. Lett.* **2003**, *381*, 512–518.

(158) Tuffery, P.; Etchebest, C.; Hazout, S.; Lavery, R. A new approach to the rapid determination of protein side chain conformations. *J. Biomol. Struct. Dyn.* **1991**, *8*, 1267–89.

(159) Grosdidier, A.; Zoete, V.; Michielin, O. SwissDock, a protein-small molecule docking web service based on EADock DSS. *Nucleic Acids Res.* **2011**, *39*, W270–W277.

(160) Lee, M. S.; Feig, M.; Salsbury Jr, F. R.; Brooks, C. L., III. New analytic approximation to the standard molecular volume definition and its application to generalized Born calculations. *J. Comput. Chem.* **2003**, *24*, 1348–1356.

(161) Stewart, J. J. P. Optimization of parameters for semiempirical methods I. method. *J. Comput. Chem.* **1989**, *10*, 209–220.

(162) Elstner, M. The SCC-DFTB method and its application to biological systems. *Theor. Chem. Acc.* **2006**, *116*, 316–325.

(163) Cui, Q.; Elstner, M.; Kaxiras, E.; Frauenheim, T.; Karplus, M. A QM/MM implementation of the self-consistent charge density functional tight binding (SCC-DFTB) method. *J. Phys. Chem. B* **2001**, *105*, 569–585.

(164) Gaus, M.; Cui, Q.; Elstner, M. DFTB-3rd: extension of the self-consistent-charge density-functional tight-binding method SCC-DFTB. *J. Chem. Theory Comp.* **2011**, *7*, 931–948.

(165) Nam, K. Acceleration of semiempirical quantum mechanical calculations by extended Lagrangian molecular dynamics approach. *J. Chem. Theory Comp.* **2013**, *9*, 3393–3403.

(166) Toledo Warshaviak, D.; Golan, G.; Borrelli, K. W.; Zhu, K.; Kalid, O. Structure-based virtual screening approach for discovery of covalently bound ligands. *J. Chem. Inf. and Mod.* **2014**, *54*, 1941–1950.
(167) Zhu, K.; Borrelli, K. W.; Greenwood, J. R.; Day, T.; Abel, R.; Farid, R. S.; Harder, E. Docking covalent inhibitors: A parameter free approach to pose prediction and scoring. *J. Chem. Inf. and Mod. 2014*, *54*, 1932–1940.

(168) Wang, L.; Chambers, J.; Abel, R. *Biomolecular Simulations*; Meth. Mol. Biol. 2019, 202, 201–232.

(169) Gapsys, V.; Pérez-Benito, L.; Aldeghi, M.; Seeliger, D.; van Vlijmen, H.; Tresadern, G.; de Groot, B. L. Large scale relative protein ligand binding affinities using non-equilibrium alchemy. *Chem. Sci.* 2020, *11*, 1140–1152.

(170) Gapsys, V.; Michielssens, S.; Seeliger, D.; de Groot, B. L. Accurate and rigorous prediction of the changes in protein free energies in a large-scale mutation scan. *Angew. Chem. 2016*, *55*, 7364–7368.

(171) Hayes, R. L.; Nixon, C. F.; Marqusee, S.; Brooks, C. L., III. Selection pressures on evolution of ribonuclease H explored with rigorous free-energy-based design. *Proc. Natl. Acad. Sci. USA 2024*, *121*, No. e2312029121.

(172) Goh, G. B.; Hulbert, B. S.; Zhou, H.; Brooks, C. L., III. Constant pH molecular dynamics of proteins in explicit solvent with proton tautomerism. *PROTEINS: Struct., Func. and Bioinf.* 2014, *82*, 1319–1331.

(173) Kong, X.; Brooks, C. L., III. Lambda-Dynamics: a new approach to free energy calculations. *J. Chem. Phys. 1996*, *105*, 2414–2423.

(174) Thiel, A. C.; Speranza, M. J.; Jadhav, S.; Stevens, L. L.; Unruh, D. K.; Ren, P.; Ponder, J. W.; Shen, J.; Schnieders, M. J. Constant-pH Simulations with the Polarizable Atomic Multipole AMOEBA Force Field. *J. Chem. Theory Comp. 2024*, *20*, 2921–2933.

(175) Hayes, R. L.; Armacost, K. A.; Vilseck, J. Z.; Brooks, C. L., III. Adaptive landscape flattening accelerates sampling of alchemical space in multisite $\lambda$ dynamics. *J. Phys. Chem. B 2017*, *121*, 3626–3635.

(176) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. *J. Chem. Phys.* 1995, *103*, 8577–8593.

(177) Hayes, R. L.; Vilseck, J. Z.; Brooks, C. L., III. Approaching protein design with multisite $\lambda$ dynamics: accurate and scalable mutational folding free energies in T4 lysozyme. *Protein Sci.* 2018, *27*, 1910–1922.

(178) Hayes, R. L.; Vilseck, J. Z.; Brooks, C. L., III. Addressing intersite coupling unlocks large combinatorial chemical spaces for alchemical free energy methods. *J. Chem. Theory Comp. 2022*, *18*, 2114–2123.

(179) Huang, Y.; Chen, W.; Wallace, J. A.; Shen, J. K. All-atom continuous constant pH molecular dynamics with particle mesh Ewald and titratable water. *J. Chem. Theory Comp. 2016*, *12*, 5411–5421.

(180) Ding, X.; Vilseck, J. Z.; Hayes, R. L.; Brooks, C. L., III. Gibbs sampler-based $\lambda$-dynamics and Rao-Blackwell estimator for alchemical free energy calculation. *J. Chem. Theory Comp. 2017*, *13*, 2501–2510.

(181) Vilseck, J. Z.; Ding, X.; Hayes, R. L.; Brooks, C. L., III. Generalizing the discrete Gibbs sampler-based $\lambda$-dynamics approach for multisite sampling of many ligands. *J. Chem. Theory Comp. 2021*, *17*, 3895–3907.

(182) Robo, M. T.; Hayes, R. L.; Ding, X.; Pulawski, B.; Vilseck, J. Z. Fast free energy estimates from $\lambda$-dynamics with bias-updated Gibbs sampling. *Nat. Commun.* 2023, *14*, 8515.

(183) Liu, S.; Wang, L.; Mobley, D. L. Is ring breaking feasible in relative binding free energy calculations? *J. Chem. Inf. and Mod. 2015*, *55*, 727–735.

(184) Wang, L.; Deng, Y.; Wu, Y.; Kim, B.; LeBard, D. N.; Wandschneider, D.; Beachy, M.; Friesner, R. A.; Abel, R. Accurate modeling of scaffold hopping transformations in drug discovery. *J. Chem. Theory Comp. 2017*, *13*, 42–54.

(185) Hayes, R. L.; Brooks, C. L., III. A strategy for proline and glycine mutations to proteins with alchemical free energy calculations. *J. Comput. Chem. 2021*, *42*, 1088–1094.

(186) Vilseck, J. Z.; Sohail, N.; Hayes, R. L.; Brooks, C. L., III. Overcoming challenging substituent perturbations with multisite $\lambda$-

dynamics: a case study targeting $\beta$-Secretase 1. *J. Phys. Chem. Lett.* 2019, *10*, 4875–4880.

(187) Deng, Y. Q.; Roux, B. Computations of standard binding free energies with molecular dynamics simulations. *J. Chem. Theory Comp.* 2009, *113*, 2234–2246.

(188) Jiang, W.; Hodoscek, M.; Roux, B. Computation of absolute hydration and binding free energy with free energy perturbation distributed replica-exchange molecular dynamics. *J. Chem. Theory Comp.* 2009, *5*, 2583–2588.

(189) Jiang, W.; Roux, B. Free energy perturbation hamiltonian replica-exchange molecular dynamics (FEP/H-REMD) for absolute ligand binding free energy calculations. *J. Chem. Theory Comp.* 2010, *6*, 2559–2565.

(190) Lin, Y. L.; Meng, Y.; Jiang, W.; Roux, B. Explaining why Gleevec is a specific and potent inhibitor of Abl kinase. *Proc. Natl. Acad. Sci. USA 2013*, *110*, 1664–1669.

(191) Lin, Y. L.; Roux, B. Computational analysis of the binding specificity of Gleevec to Abl, c-Kit, Lck, and c-Src tyrosine kinases. *J. Am. Chem. Soc.* 2013, *135*, 14741–14753.

(192) Lin, Y. L.; Meng, Y.; Huang, L.; Roux, B. Computational study of Gleevec and G6G reveals molecular determinants of kinase inhibitor selectivity. *J. Am. Chem. Soc.* 2014, *136*, 14753–14762.

(193) Jiang, W.; Luo, Y.; Maragliano, L.; Roux, B. Calculation of free energy landscape in multi-dimensions with Hamiltonian-exchange umbrella sampling on petascale supercomputer. *J. Chem. Theory Comp.* 2012, *8*, 4672–4680.

(194) Deng, Y. Q.; Roux, B. Computation of binding free energy with molecular dynamics and grand canonical Monte Carlo simulations. *J. Chem. Phys.* 2008, *128*, 115103.

(195) Im, W.; Bernèche, S.; Roux, B. Generalized solvent boundary potentials for computer simulations. *J. Chem. Phys.* 2001, *114*, 2924–2937.

(196) Ge, X.; Roux, B. Calculation of the standard binding free energy of sparsomycin to the ribosomal peptidyl-transferase P-site using molecular dynamics simulations with restraining potentials. *J. Mol. Recognit.* 2010, *23*, 128–141.

(197) Ge, X.; Roux, B. Absolute binding free energy calculations of sparsomycin analogs to the bacterial ribosome. *J. Phys. Chem. B 2010*, *114*, 9525–9539.

(198) Shivakumar, D.; Deng, Y.; Roux, B. Computations of absolute solvation free energies of small molecules using explicit and implicit solvent model. *J. Chem. Theory Comp.* 2009, *5*, 919–930.

(199) Loeffler, H. H.; Bosisio, S.; Duarte Ramos Matos, G.; Suh, D.; Roux, B.; Mobley, D. L.; Michel, J. Reproducibility of free energy calculations across different molecular simulation software packages. *J. Chem. Theory Comp.* 2018, *14*, 5567–5582.

(200) Woo, H.; Roux, B. Calculation of absolute protein-ligand binding free energy from computer simulations. *Proc. Natl. Acad. Sci. USA 2005*, *102*, 6825–6830.

(201) Gan, W.; Roux, B. Binding specificity of SH2 domains: insight from free energy simulations. *Proteins* 2009, *74*, 996–1007.

(202) Im, W.; Beglov, D.; Roux, B. Continuum solvation model: electrostatic forces from numerical solutions to the Poisson-Bolzmann equation. *Comput. Phys. Commun.* 1998, *111*, 59–75.

(203) Nandigrami, P.; Szczepaniak, F.; Boughter, C. T.; Dehez, F.; Chipot, C.; Roux, B. Computational assessment of protein-protein binding specificity within a family of synaptic surface receptors. *J. Phys. Chem. B 2022*, *126*, 7510–7527.

(204) Chen, Y.; Roux, B. Efficient hybrid non-equilibrium molecular dynamics - Monte Carlo simulations with symmetric momentum reversal. *J. Chem. Phys.* 2014, *141*, 114107.

(205) Chen, Y.; Roux, B. Generalized metropolis acceptance criterion for hybrid non-equilibrium molecular dynamics - Monte Carlo simulations. *J. Chem. Phys.* 2015, *142*, 024101.

(206) Chen, Y.; Roux, B. Enhanced sampling of an atomic model with hybrid nonequilibrium molecular dynamics: Monte Carlo simulations guided by a coarse-grained model. *J. Chem. Theory Comp.* 2015, *11*, 3572–3583.
(207) Chen, Y.; Roux, B. Constant-pH hybrid nonequilibrium molecular dynamics-Monte Carlo simulation method. *J. Chem. Theory Comp.* **2015**, *11*, 3919-3931.

(208) Suh, D.; Jo, S.; Jiang, W.; Chipot, C.; Roux, B. String method for protein-protein binding free-energy calculations. *J. Chem. Theory Comp.* **2019**, *15*, 5829-5844.

(209) Chen, H.; Chipot, C. Enhancing sampling with free-energy calculations. *Curr. Op. Struct. Biol.* **2022**, *77*, 102497.

(210) Liu, X.; Brooks, C. L., III. Enhanced sampling of buried charges in free energy calculations using replica exchange with charge tempering. *J. Chem. Theory Comp.* **2024**, *20*, 1051-1061.

(211) Lee, J.; Miller, B. T.; Brooks, B. R. Computational scheme for pH-dependent binding free energy calculation with explicit solvent. *Protein Sci.* **2016**, *25*, 231-243.

(212) Allen, T. W.; Andersen, O. S.; Roux, B. Ion permeation through a narrow channel: using gramicidin to ascertain all-atom molecular dynamics potential of mean force methodology and biomolecular force fields. *Biophys. J.* **2006**, *90*, 3447-3468.

(213) Wieder, M.; Fleck, M.; Braunsfeld, B.; Boresch, S. Alchemical free energy simulations without speed limits. A generic framework to calculate free energy differences independent of the underlying molecular dynamics program. *J. Comput. Chem.* **2022**, *43*, 1151-1160.

(214) Bruckner, S.; Boresch, S. Efficiency of alchemical free energy simulations. II. improvements for thermodynamic integration. *J. Comput. Chem.* **2011**, *32*, 1320-1333.

(215) Reinhardt, M.; Grubmüller, H. Determining free-energy differences through variationally derived intermediates. *J. Chem. Theory Comp.* **2020**, *16*, 3504-3512.

(216) König, G.; Brooks, B. R.; Thiel, W.; York, D. M. On the convergence of multi-scale free energy simulations. *Mol. Sim.* **2018**, *44*, 1062-1081.

(217) Blondel, A. Ensemble variance in free energy calculations by thermodynamic integration: theory, optimal “alchemical” path, and practical solutions. *J. Comput. Chem.* **2004**, *25*, 985-993.

(218) Pham, T. T.; Shirts, M. R. Optimal pairwise and non-pairwise alchemical pathways for free energy calculations of molecular transformation in solution phase. *J. Chem. Phys.* **2012**, *136*, 124120.

(219) König, G.; Ries, B.; Hünenberger, P. H.; Riniker, S. Efficient alchemical intermediate states in free energy calculations using λ-enveloping distribution sampling. *J. Chem. Theory Comp.* **2021**, *17*, 5805-5815.

(220) König, G.; Glaser, N.; Schroeder, B.; Kubincová, A.; Hünenberger, P. H.; Riniker, S. An alternative to conventional λ-intermediate states in alchemical free energy calculations: λ-enveloping distribution sampling. *J. Chem. Inf. and Mod.* **2020**, *60*, 5407-5423.

(221) Narayana, N.; Matthews, D. A.; Howell, E. E.; Xuong, N.-h. A plasmid-encoded dihydrofolate reductase from trimethoprim-resistant bacteria has a novel D2-symmetric active site. *Nat. Struct. Biol.* **1995**, *2*, 1018-1025.

(222) Dam, J.; Rose, T.; Goldberg, M. E.; Blondel, A. Complementation between dimeric mutants as a probe of dimer-dimer interactions in tetrameric dihydrofolate reductase encoded by R67 plasmid of E. coli. *J. Mol. Biol.* **2000**, *302*, 235-250.

(223) Dam, J.; Blondel, A. Effect of multiple symmetries on the association of R67 DHFR subunits bearing interfacial complementing mutations. *Protein Sci.* **2004**, *13*, 1-14.

(224) Straatsma, T. P.; Berendsen, H. J. C.; Stam, A. J. Estimation of statistical errors in molecular simulation calculations. *Mol. Phys.* **1986**, *57*, 89-95.

(225) Lee, M. S.; Salsbury, F. R.; Brooks, C. L., III. Constant-pH molecular dynamics using continuous titration coordinates. *PROTEINS: Struct., Func. and Bioinf.* **2004**, *56*, 738-752.

(226) Khandogin, J.; Brooks, C. L., III. Constant pH molecular dynamics with proton tautomerism. *Biophys. J.* **2005**, *89*, 141-157.

(227) Wallace, J. A.; Shen, J. K. Continuous constant pH molecular dynamics in explicit solvent with pH-based replica exchange. *J. Chem. Theory Comp.* **2011**, *7*, 2617-2629.

(228) Huang, Y.; Chen, W.; Dotson, D. L.; Beckstein, O.; Shen, J. Mechanism of pH-dependent activation of the sodium-proton antiporter NhaA. *Nat. Commun.* **2016**, *7*, 12940.

(229) Im, W.; Feig, M.; Brooks, C. L., III. An implicit membrane generalized Born theory for the study of structure, stability, and interactions of membrane proteins. *Biophys. J.* **2003**, *85*, 2900-2918.

(230) Wallace, J. A.; Shen, J. K. Charge-leveling and proper treatment of long-range electrostatics in all-atom molecular dynamics at constant pH. *J. Chem. Phys.* **2012**, *137*, 184105.

(231) Chen, W.; Wallace, J. A.; Yue, Z.; Shen, J. K. Introducing titratable water to all-atom molecular dynamics at constant pH. *Biophys. J.* **2013**, *105*, L15-L17.

(232) Morrow, B. H.; Payne, G. F.; Shen, J. pH-responsive self-assembly of polysaccharide through a rugged energy landscape. *J. Am. Chem. Soc.* **2015**, *137*, 13024-13030.

(233) Huang, Y.; Yue, Z.; Tsai, C.-C.; Henderson, J. A.; Shen, J. Predicting catalytic proton donors and nucleophiles in enzymes: how adding dynamics helps elucidate the structure-function relationships. *J. Phys. Chem. Lett.* **2018**, *9*, 1179-1184.

(234) Baptista, A. M.; Teixeira, V. H.; Soares, C. M. Constant-pH molecular dynamics using stochastic titration. *J. Chem. Phys.* **2002**, *117*, 4184-4200.

(235) Bürgi, R.; Kollman, P. A.; Van Gunsteren, W. F. Simulating proteins at constant pH: An approach combining molecular dynamics and Monte Carlo simulation. *PROTEINS: Struct., Func. and Bioinf.* **2002**, *47*, 469-480.

(236) Dlugosz, M.; Antosiewicz, J. M. Constant-pH molecular dynamics simulations: a test case of succinic acid. *Chem. Phys.* **2004**, *302*, 161-170.

(237) Mongan, J.; Case, D. A.; McCammon, J. A. Constant pH molecular dynamics in generalized Born implicit solvent. *J. Comput. Chem.* **2004**, *25*, 2038-2048.

(238) Itoh, S. G.; Damjanović, A.; Brooks, B. R. pH replica-exchange method based on discrete protonation states. *PROTEINS: Struct., Func. and Bioinf.* **2011**, *79*, 3420-3436.

(239) Swails, J. M.; York, D. M.; Roitberg, A. E. Constant pH replica exchange molecular dynamics in explicit solvent using discrete protonation states: implementation, testing, and validation. *J. Chem. Theory Comp.* **2014**, *10*, 1341-1352.

(240) Stern, H. A. Molecular simulation with variable protonation states at constant pH. *J. Chem. Phys.* **2007**, *126*, 164112.

(241) Radak, B. K.; Chipot, C.; Suh, D.; Jo, S.; Jiang, W.; Phillips, J. C.; Schulten, K.; Roux, B. Constant-pH molecular dynamics simulations for large biomolecular systems. *J. Chem. Theory Comp.* **2017**, *13*, 5933-5944.

(242) Donnini, S.; Tegeler, F.; Groenhof, G.; Grubmüller, H. Constant pH molecular dynamics in explicit solvent with λ-dynamics. *J. Chem. Theory Comp.* **2011**, *7*, 1962-1978.

(243) Damjanović, A.; Miller, B. T.; Okur, A.; Brooks, B. R. Reservoir pH replica exchange. *J. Chem. Phys.* **2018**, *149*, 072321.

(244) Damjanovic, A.; Chen, A. Y.; Rosenberg, R. L.; Roe, D. R.; Wu, X.; Brooks, B. R. Protonation state of the selectivity filter of bacterial voltage-gated sodium channels is modulated by ions. *PROTEINS: Struct., Func. and Bioinf.* **2020**, *88*, 527-539.

(245) Simonson, T.; Carlsson, J.; Case, D. A. Proton binding to proteins: p K a calculations with explicit and implicit solvent models. *J. Am. Chem. Soc.* **2004**, *126*, 4167-4180.

(246) Ghosh, N.; Cui, Q. PKa of Residue 66 in Staphylococal nuclease. I. Insights from QM/MM Simulations with Conventional Sampling. *J. Phys. Chem. B* **2008**, *112*, 8387-8397.

(247) Chen, A. Y.; Brooks, B. R.; Damjanovic, A. Ion channel selectivity through ion-modulated changes of selectivity filter PKa values. *Proc. Natl. Acad. Sci. USA* **2023**, *120*, No. e2220343120.

(248) Lazaridis, T.; Hummer, G. Classical molecular dynamics with mobile protons. *J. Chem. Inf. and Mod.* **2017**, *57*, 2833-2845.

(249) Lazaridis, T.; Sepehri, A. Amino acid deprotonation rates from classical force fields. *J. Chem. Phys.* **2022**, *157*, 085101.

(250) Lazaridis, T. Molecular origins of asymmetric proton conduction in the influenza M2 channel. *Biophys. J.* **2023**, *122*, 90-98.
(251) Lazaridis, T. Proton paths in models of the Hv1 proton channel. *J. Phys. Chem. B* **2023**, *127*, 7937–7945.

(252) Hansmann, U. H. E. Parallel tempering algorithm for conformational studies of biological molecules. *Chem. Phys. Lett.* **1997**, *281*, 140–150.

(253) Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. *Chem. Phys. Lett.* **1999**, *314*, 141–151.

(254) Bernardi, R. C.; Melo, M. C. R.; Schulten, K. Enhanced sampling techniques in molecular dynamics simulations of biological systems. *Biochimica et Biophysica Acta - General Subjects* **2015**, *1850*, 872–877.

(255) Paschek, D.; García, A. E. Reversible temperature and pressure denaturation of a protein fragment: a replica exchange molecular dynamics simulation study. *Phys. Rev. Lett.* **2004**, *93*, 10–13.

(256) Wu, X.; Hodoscek, M.; Brooks, B. R. Replica exchanging self-guided Langevin dynamics for efficient and accurate conformational sampling. *J. Chem. Phys.* **2012**, *137*, 1–13.

(257) Okur, A.; Miller, B. T.; Joo, K.; Lee, J.; Brooks, B. R. Generating reservoir conformations for replica exchange through the use of the conformational space annealing method. *J. Chem. Theory Comp.* **2013**, *9*, 1115–1124.

(258) Okur, A.; Roe, D. R.; Cui, G.; Hornak, V.; Simmerling, C. Improving convergence of replica-exchange simulations through coupling to a high-temperature structure reservoir. *J. Chem. Theory Comp.* **2007**, *3*, 557–568.

(259) Roitberg, A. E.; Okur, A.; Simmerling, C. Coupling of replica exchange simulations to a non-Boltzmann structure reservoir. *J. Phys. Chem. B* **2007**, *111*, 2415–2418.

(260) Sugita, Y.; Kitao, A.; Okamoto, Y. Multidimensional replica-exchange method for free-energy calculations. *J. Chem. Phys.* **2000**, *113*, 6042–6051.

(261) Bergonzo, C.; Henriksen, N. M.; Roe, D. R.; Swails, J. M.; Roitberg, A. E.; Cheatham, T. E.; Cheatham, T. E., III. Multidimensional replica exchange molecular dynamics yields a converged ensemble of an RNA tetranucleotide. *J. Chem. Theory Comp.* **2014**, *10*, 492–499.

(262) Woodcock III, H. L.; Miller, B. T.; Hodoscek, M.; Okur, A.; Larkin, J. D.; Ponder, J. W.; Brooks, B. R. MSCALE: a general utility for multiscale modeling. *J. Chem. Theory Comp.* **2011**, *7*, 1208–1219.

(263) Schlitter, J.; Engels, M.; Krüger, P.; Jacoby, E.; Wollmer, A. Targeted molecular-dynamics simulation of conformational change - application to the T ↔ R transition in insulin. *Mol. Sim.* **1993**, *10*, 291–308.

(264) Apostolakis, J.; Ferrara, P.; Caflisch, A. Calculation of conformational transitions and barriers in solvated systems: application to the alanine dipeptide in water. *J. Chem. Phys.* **1999**, *110*, 2099–2108.

(265) van der Vaart, A.; Karplus, M. Simulation of conformational transitions by the restricted perturbation-targeted molecular dynamics method. *J. Chem. Phys.* **2005**, *122*, 114903.

(266) van der Vaart, A.; Karplus, M. Minimum free energy pathways and free energy profiles for conformational transitions based on atomistic molecular dynamics simulations. *J. Chem. Phys.* **2007**, *126*, 164106.

(267) Ovchinnikov, V.; Karplus, M. Analysis and elimination of a bias in targeted molecular dynamics simulations of conformational transitions: application to calmodulin. *J. Phys. Chem. B* **2012**, *116*, 8584–8603.

(268) Hooft, R.; Vaneijk, B.; Kroon, J. An adaptive umbrella sampling procedure in conformational-analysis using molecular-dynamics and its application to glycol. *J. Chem. Phys.* **1992**, *97*, 6690–6694.

(269) Bartels, C.; Karplus, M. Multidimensional adaptive umbrella sampling: applications to main chain and side chain peptide conformations. *J. Comput. Chem.* **1997**, *18*, 1450–1462.

(270) Ovchinnikov, V.; Cecchini, M.; Karplus, M. A simplified confinement method for calculating absolute free energies and free energy and entropy differences. *J. Phys. Chem. B* **2013**, *117*, 750–762.

(271) Esque, J.; Cecchini, M. Accurate calculation of conformational free energy differences in explicit water: the confinement-solvation free energy approach. *J. Phys. Chem. B* **2015**, *119*, 5194–5207.

(272) Villemot, F.; Peguero-Tejada, A.; van der Vaart, A. Calculation of conformational free energies by confinement simulations in explicit water with implicit desolvation. *Mol. Sim.* **2018**, *44*, 1082–1089.

(273) van der Vaart, A.; Orndorff, P. B.; Le Phan, S. T. Calculation of conformational free energies with the focused confinement method. *J. Chem. Theory Comp.* **2019**, *15*, 6760–6768.

(274) Orndorff, P. B.; Le Phan, S. T.; Li, K. H.; van der Vaart, A. Conformational free-energy differences of large solvated systems with the focused confinement method. *J. Chem. Theory Comp.* **2020**, *16*, 5163–5173.

(275) Maragakis, P.; van der Vaart, A.; Karplus, M. Gaussian-mixture umbrella sampling. *J. Phys. Chem. B* **2009**, *113*, 4664–4673.

(276) Spiriti, J.; Binder, J. K.; Levitus, M.; van der Vaart, A. Cy3-DNA stacking interactions strongly depend on the identity of the terminal basepair. *Biophys. J.* **2011**, *100*, 1049–1057.

(277) Andricioaei, I.; Straub, J. E. On Monte Carlo and molecular dynamics methods inspired by Tsallis statistics: methodology, optimization, and application to atomic clusters. *J. Chem. Phys.* **1997**, *107*, 9117–9124.

(278) Liu, P.; Kim, B.; Friesner, R. A.; Berne, B. J. Replica exchange with solute tempering: a method for sampling biological systems in explicit water. *Proc. Natl. Acad. Sci. USA* **2005**, *102*, 13749–13754.

(279) Wang, L.; Friesner, R. A.; Berne, B. J. Replica exchange with solute scaling: a more efficient version of replica exchange with solute tempering (REST2). *J. Phys. Chem. B* **2011**, *115*, 9431–9438.

(280) Kamberaj, H.; van der Vaart, A. Multiple scaling replica exchange for the conformational sampling of biomolecules in explicit water. *J. Chem. Phys.* **2007**, *127*, 234102.

(281) Gray, G. M.; Thiessen, B.; van der Vaart, A. Secondary structure of peptides mimicking the Gly-rich regions of major ampullate spidroin protein 1 and 2. *Biophys. Chem.* **2022**, *284*, 106783.

(282) Elber, R.; Karplus, M. A method for determining reaction paths in large biomolecules: application to myoglobin. *Chem. Phys. Lett.* **1987**, *139*, 375–380.

(283) Jónsson, G. M.; Jacobsen, K. W. In *Classical and Quantum Dynamics in Condensed Phase Simulations*; Berne, B., Ciccotti, G., Coker, D., Eds.; World Scientific: Singapore, 1998; pp 385–404.

(284) Neria, E.; Fischer, S.; Karplus, M. Simulation of activation free energies in molecular systems. *J. Chem. Phys.* **1996**, *105*, 1902–1921.

(285) Ren, W.; Vanden-Eijnden, E.; Maragakis, P.; E, W. Transition pathways in complex systems: application of the finite-temperature string method to the alanine dipeptide. *J. Chem. Phys.* **2005**, *123*, 134109.

(286) Huo, S.; Straub, J. E. The MaxFlux algorithm for calculating variationally optimized reaction paths for conformational transitions in many body systems at finite temperature. *J. Chem. Phys.* **1997**, *107*, 5000–5006.

(287) Branduardi, D.; Gervasio, F. L.; Parrinello, M. From a to b in free energy space. *J. Chem. Phys.* **2007**, *126*, 054103.

(288) E, W.; Ren, W.; Vanden-Eijnden, E. Simplified and improved string method for computing the minimum energy paths in barrier-crossing events. *J. Chem. Phys.* **2007**, *126*, 164103.

(289) Maragliano, L.; Fischer, A.; Vanden-Eijnden, E.; Ciccotti, G. String method in collective variables: minimum free energy paths and isocommittor surfaces. *J. Chem. Phys.* **2006**, *125*, 024106.

(290) Vanden-Eijnden, E.; Venturoli, M. Markovian milestoning with Voronoi tessellations. *J. Chem. Phys.* **2009**, *130*, 194101.

(291) Gardiner, C. W. *Handbook of stochastic methods*, 3rd ed.; Springer-Verlag: Berlin, 2003.

(292) E, W.; Vanden-Eijnden, E. In *Multiscale Modeling and Simulation*; Attinger, S., Koumoutsakos, P., Eds.; Springer-Verlag: Berlin, 2004.

(293) E, W.; Ren, W.; Vanden-Eijnden, E. String method for the study of rare events. *Phys. Rev. B* **2002**, *66*, 052301.
(294) E, W.; Ren, W.; Vanden-Eijnden, E. Transition pathways in complex systems: reaction coordinates, isocommittor surfaces, and transition tubes. *Chem. Phys. Lett.* **2005**, *413*, 242–247.

(295) E, W.; Vanden-Eijnden, E. Transition-path theory and path-finding algorithms for the study of rare events. *Annu. Rev. Phys. Chem.* **2010**, *61*, 391–420.

(296) Ovchinnikov, V.; Karplus, M. Investigations of α-helix ↔ β-sheet transition pathways in a miniprotein using the finite-temperature string method. *J. Chem. Phys.* **2014**, *140*, 175103.

(297) Lu, X.; Fang, D.; Ito, S.; Okamoto, Y.; Ovchinnikov, V.; Cui, Q; QM/MM free energy simulations: recent progress and challenges. *Mol. Sim.* **2016**, *42*, 1056–1078.

(298) Pan, A.; Sezer, D.; Roux, B. Finding transition pathways using the string method with swarms of trajectories. *J. Phys. Chem. B* **2008**, *112*, 3432–3440.

(299) Maragliano, L.; Roux, B.; Vanden-Eijnden, E. Comparison between mean forces and swarms-of-trajectories string methods. *J. Chem. Theory Comp.* **2014**, *10*, 524–533.

(300) Roux, B. String method with swarms-of-trajectories, mean drifts, lag time, and committor. *J. Phys. Chem. A* **2021**, *125*, 7558–7571.

(301) Gan, W.; Yang, S.; Roux, B. Atomistic view of the conformational activation of Src kinase using the string method with swarms-of-trajectories. *Biophys. J.* **2009**, *97*, L8–L10.

(302) Meng, Y.; Lin, Y.-l.; Roux, B. Computational study of the “DFG-flip” conformational transition in c-Abl and c-Src tyrosine kinases. *J. Phys. Chem. B* **2015**, *119*, 1443–1456.

(303) Harder, E.; Walters, D. E.; Bodnar, Y. D.; Faibish, R. S.; Roux, B. Molecular dynamics study of a polymeric reverse osmosis membrane. *J. Phys. Chem. B* **2009**, *113*, 10177–10182.

(304) Luo, Y.; Harder, E.; Faibish, R. S.; Roux, B. Computer simulations of water flux and salt permeability of the reverse osmosis FT-30 aromatic polyamide membrane. *J. Membr. Sci.* **2011**, *384*, 1–9.

(305) Dickson, B. M.; Huang, H.; Post, C. B. Unrestrained computation of free energy along a path. *J. Phys. Chem. B* **2012**, *116*, 11046–11055.

(306) Dickson, B. M. Approaching a parameter-free metadynamics. *Phys. Rev. E* **2011**, *84*, 037701.

(307) E, W.; Ren, W.; Vanden-Eijnden, E. Finite temperature string method for the study of rare events. *J. Phys. Chem. B* **2005**, *109*, 6688–6693.

(308) Wu, H.; Post, C. B. Protein conformational transitions from all-atom adaptively biased path optimization. *J. Chem. Theory Comp.* **2018**, *14*, 5372–5382.

(309) Wu, H.; Huang, H.; Post, C. B. All-atom adaptively biased path optimization of Src kinase conformational inactivation: switched electrostatic network in the concerted motion of αC helix and the activation loop. *J. Chem. Phys.* **2020**, *153*, 175101.

(310) Lee Woodcock, H.; Hodoscek, M.; Sherwood, P.; Lee, Y. S.; Schaefer, H. F., III; Brooks, B. R. Exploring the quantum mechanical/molecular mechanical replica path method: a pathway optimization of the chorismate to prephenate Claisen rearrangement catalyzed by chorismate mutase. *Theor. Chem. Acc.* **2003**, *109*, 140–148.

(311) Henkelman, G.; Jónsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. *J. Chem. Phys.* **2000**, *113*, 9978–9985.

(312) Chu, J.-W.; Trout, B. L.; Brooks, B. R. A super-linear minimization scheme for the nudged elastic band method. *J. Chem. Phys.* **2003**, *119*, 12708.

(313) Brokaw, J. B.; Haas, K. R.; Chu, J.-W. Reaction path optimization with holonomic constraints and kinetic energy potentials. *J. Chem. Theory Comp.* **2009**, *5*, 2050–2061.

(314) Haas, K.; Chu, J.-W. Decomposition of energy and free energy changes by following the flow of work along reaction path. *J. Chem. Phys.* **2009**, *131*, 144105.

(315) Cho, H. M.; Gross, A. S.; Chu, J.-W. Dissecting force interactions in cellulose deconstruction reveals the required solvent versatility for overcoming biomass recalcitrance. *J. Am. Chem. Soc.* **2011**, *133*, 14033–14041.

(316) Lin, Y.; Beckham, G. T.; Himmel, M. E.; Crowley, M. F.; Chu, J.-W. Endoglucanase peripheral loops facilitate complexation of glucan chains on cellulose via adaptive coupling to the emergent substrate structures. *J. Phys. Chem. B* **2013**, *117*, 10750–10758.

(317) Brooks, B. R.; Janežič, D.; Karplus, M. Harmonic analysis of large systems. *I. methodology*. *J. Comput. Chem.* **1995**, *16*, 1522–1542.

(318) Chen, Y.-T.; Yang, H.; Chu, J.-W. Structure-mechanics statistical learning unravels the linkage between local rigidity and global flexibility in nucleic acids. *Chem. Sci.* **2020**, *11*, 4969–4979.

(319) Raj, N.; Click, T.; Yang, H.; Chu, J.-W. Mechanical couplings of protein backbone and side chains exhibit scale-free network properties and specific hotspots for function. *Comput. Struct. Biotechnol. J.* **2021**, *19*, 5309–5320.

(320) Raj, N.; Click, T. H.; Yang, H.; Chu, J.-W. Structure-mechanics statistical learning uncovers mechanical relay in proteins. *Chem. Sci.* **2022**, *13*, 3688–3696.

(321) Chen, Y.-T.; Yang, H.; Chu, J.-W. Mechanical codes of chemical-scale specificity in DNA motifs. *Chem. Sci.* **2023**, *14*, 10155–10166.

(322) Glowacki, D. R.; Paci, E.; Shalashilin, D. V. Boxed molecular dynamics: a simple and general technique for accelerating rare event kinetics and mapping free energy in large molecular systems. *J. Phys. Chem. B* **2009**, *113*, 16603–16611.

(323) Glowacki, D. R.; Paci, E.; Shalashilin, D. V. Boxed molecular dynamics: decorrelation time scales and the kinetic master equation. *J. Chem. Theory Comp.* **2011**, *7*, 1244–1252.

(324) Faradjian, A. K.; Elber, R. Computing time scales from reaction coordinates by milestoneing. *J. Chem. Phys.* **2004**, *120*, 10880–10889.

(325) O'Connor, M.; Paci, E.; McIntosh-Smith, S.; Glowacki, D. R. Adaptive free energy sampling in multidimensional collective variable space using boxed molecular dynamics. *Faraday Discuss* **2016**, *195*, 395–419.

(326) Darve, E.; Pohorille, A. Calculating free energies using average force. *J. Chem. Phys.* **2001**, *115*, 9169–9183.

(327) Darve, E.; Wilson, M. A.; Pohorille, A. Calculating free energies using a scaled-force molecular dynamics algorithm. *Molecular Simulation* **2002**, *28*, 113–144.

(328) Hénin, J.; Chipot, C. Overcoming free energy barriers using unconstrained molecular dynamics simulations. *J. Chem. Phys.* **2004**, *121*, 2904–2914.

(329) Kirkwood, J. G. Statistical mechanics of fluid mixtures. *J. Chem. Phys.* **1935**, *3*, 300–313.

(330) Fu, H.; Shao, X.; Chipot, C.; Cai, W. Extended adaptive biasing force algorithm. An on-the-fly implementation for accurate free-energy calculations. *J. Chem. Theory Comp.* **2016**, *12*, 3506–3513.

(331) Lesage, A.; Leliévre, T.; Stoltz, G.; Hénin, J. Smoothed biasing forces yield unbiased free energies with the extended-system adaptive biasing force method. *J. Phys. Chem. B* **2017**, *121*, 3676–3685.

(332) Darve, E.; Rodriguez-Gómez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. *J. Chem. Phys.* **2008**, *128*, 144120.

(333) Leliévre, T.; Rousset, M.; Stoltz, G. *Free Energy Computations*; Imperial College Press: 2010.

(334) Leliévre, T.; Rousset, M.; Stoltz, G. Computation of free energy profiles with parallel adaptive dynamics. *J. Chem. Phys.* **2007**, *126*, 134111.

(335) Zheng, L.; Chen, M.; Yang, W. Random walk in orthogonal space to achieve efficient free-energy simulation of complex systems. *Proc. Natl. Acad. Sci. USA* **2008**, *105*, 20227–20232.

(336) Zheng, L.; Yang, W. Practically efficient and robust free energy calculations: double-integration orthogonal space tempering. *J. Chem. Theory Comp.* **2012**, *8*, 810–823.

(337) Hénin, J.; Fiorin, G.; Chipot, C.; Klein, M. L. Exploring multidimensional free energy landscapes using time-dependent biases on collective variables. *J. Chem. Theory Comp.* **2010**, *6*, 35–47.

(338) Stone, A. J. *The theory of intermolecular forces*; Clarendon Press: Oxford, U.K., 1996; Vol. 32.
(339) Stone, A. J. Distributed multipole analysis: stability for large basis sets. *J. Chem. Theory Comp.* **2005**, *1*, 1128-1132.

(340) Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A., Jr.; et al. Current status of the AMOEBA polarizable force field. *J. Phys. Chem. B* **2010**, *114*, 2549-2564.

(341) Piquemal, J.-P.; Gresh, N.; Giessner-Prettre, C. Improved formulas for the calculation of the electrostatic contribution to the intermolecular interaction energy from multipolar expansion of the electronic distribution. *J. Phys. Chem. A* **2003**, *107*, 10353-10359.

(342) Bereau, T.; Kramer, C.; Meuwly, M. Leveraging symmetries of static atomic multipole electrostatics in molecular dynamics simulations. *J. Chem. Theory Comp.* **2013**, *9*, 5450-5459.

(343) Cardamone, S.; Popelier, P. L. A. Prediction of conformationally dependent atomic multipole moments in carbohydrates. *J. Comput. Chem.* **2015**, *36*, 2361-2373.

(344) Rein, R. On physical properties and interactions of polyatomic molecules with application to molecular recognition in biology. *Adv. Quantum Chem.* **1973**, *7*, 335-396.

(345) Maroulis, G. Electric polarizability and hyperpolarizability of carbon monoxide. *J. Phys. Chem.* **1996**, *100*, 13466-13473.

(346) Roco, J. M. M.; Calvo Hernández, A.; Velasco, S. Far-infrared permanent and induced dipole absorption of diatomic molecules in rare-gas fluids. I. spectral theory. *J. Chem. Phys.* **1995**, *103*, 9161-9174.

(347) Roco, J. M. M.; Medina, A.; Calvo Hernández, A.; Velasco, S. Far-infrared permanent and induced dipole absorption of diatomic molecules in rare-gas fluids. II. application to the CO–Ar system. *J. Chem. Phys.* **1995**, *103*, 9175-9186.

(348) Straub, J. E.; Karplus, M. Molecular dynamics study of the photodissociation of carbon monoxide from myoglobin: ligand dynamics in the 1st 10 ps. *Chem. Phys.* **1991**, *158*, 221-248.

(349) Nutt, D.; Meuwly, M. Theoretical investigation of infrared spectra and pocket dynamics of photodissociated carbonmonoxy Myoglobin. *Biophys. J.* **2003**, *85*, 3612-3623.

(350) Nutt, D.; Meuwly, M. CO migration in native and mutant myoglobin: atomistic simulations for the understanding of protein function. *Proc. Natl. Acad. Sci. USA* **2004**, *101*, 5998-6002.

(351) Plattner, N.; Meuwly, M. The role of higher CO-multipole moments in understanding the dynamics of photodissociated carbonmonoxide in myoglobin. *Biophys. J.* **2008**, *94*, 2505-2515.

(352) Clark, T.; Hennemann, M.; Murray, J. S.; Politzer, P. Halogen bonding: the σ - hole. *J. Mol. Model.* **2007**, *13*, 291-296.

(353) Murray, J. S.; Lane, P.; Politzer, P. Expansion of the σ - hole concept. *J. Mol. Model.* **2009**, *15*, 723-729.

(354) Politzer, P.; Murray, J. S.; Clark, T. Halogen bonding: an electrostatically-driven highly directional noncovalent interaction. *Phys. Chem. Chem. Phys.* **2010**, *12*, 7748-7757.

(355) Jackson, J. D. *Classical electrodynamics*; John Wiley & Sons: New York, 1998.

(356) Devereux, M.; Raghunathan, S.; Fedorov, D. G.; Meuwly, M. A novel, computationally efficient multipolar model employing distributed charges for molecular dynamics simulations. *J. Chem. Theory Comp.* **2014**, *10*, 4229-4241.

(357) Unke, O. T.; Devereux, M.; Meuwly, M. Minimal distributed charges: Multipolar quality at the cost of point charge electrostatics. *J. Chem. Phys.* **2017**, *147*, 161712.

(358) Boittier, E. D.; Devereux, M.; Meuwly, M. Molecular dynamics with conformationally dependent, distributed charges. *J. Chem. Theory Comp.* **2022**, *18*, 7544-7554.

(359) Boittier, E. D.; Devereux, M.; Töpfer, K.; Meuwly, M. Kernel-based minimally distributed charges: a conformationally dependent ESP-model for molecular simulations. *J. Chem. Theory Comp.* **2024**, in press. DOI: 10.48550/arXiv.2406.00513.

(360) Unke, O. T.; Chmiela, S.; Sauceda, H. E.; Gastegger, M.; Poltavsky, I.; Schütt, K. T.; Tkatchenko, A.; Müller, K.-R. Machine learning force fields. *Chem. Rev.* **2021**, *121*, 10142-10186.

(361) Manzhos, S.; Carrington, T., Jr Neural network potential energy surfaces for small molecules and reactions. *Chem. Rev.* **2021**, *121*, 10187-10217.

(362) Meuwly, M. Machine learning for chemical reactions. *Chem. Rev.* **2021**, *121*, 10218-10239.

(363) Deringer, V. L.; Bartók, A. P.; Bernstein, N.; Wilkins, D. M.; Ceriotti, M.; Csányi, G. Gaussian process regression for materials and molecules. *Chem. Rev.* **2021**, *121*, 10073-10141.

(364) Houston, P. L.; Qu, C.; Yu, Q.; Conte, R.; Nandi, A.; Li, J. K.; Bowman, J. M. PESPIP: software to fit complex molecular and many-body potential energy surfaces with permutationally invariant polynomials. *J. Chem. Phys.* **2023**, *158*, 044109.

(365) Braams, B. J.; Bowman, J. M. Permutationally invariant potential energy surfaces in high dimensionality. *Intern. Rev. Phys. Chem.* **2009**, *28*, 577-606.

(366) Schuett, K. T.; Sauceda, H. E.; Kindermans, P. J.; Tkatchenko, A.; Mueller, K. R. SchNet - a deep learning architecture for molecules and materials. *J. Chem. Phys.* **2018**, *148*, 241722.

(367) Zhang, L.; Han, J.; Wang, H.; Saidi, W. A.; Car, R.; E, W. End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems. *Advances in Neural Information Processing Systems.* **2018**.

(368) Chmiela, S.; Tkatchenko, A.; Sauceda, H. E.; Poltavsky, I.; Schütt, K. T.; Müller, K.-R. Machine learning of accurate energy-conserving molecular force fields. *Sci. Adv.* **2017**, *3*, No. e1603015.

(369) Sauceda, H. E.; Chmiela, S.; Poltavsky, I.; Müller, K.-R.; Tkatchenko, A. Molecular force fields with gradient-domain machine learning: construction and application to dynamics of small molecules with coupled cluster forces. *J. Chem. Phys.* **2019**, *150*, 114102.

(370) Ho, T.-S.; Rabitz, H. A general method for constructing multidimensional molecular potential energy surfaces from ab initio calculations. *J. Chem. Phys.* **1996**, *104*, 2584-2597.

(371) Unke, O. T.; Meuwly, M. Toolkit for the construction of reproducing kernel-based representations of data: application to multidimensional potential energy surfaces. *J. Chem. Inf. and Mod.* **2017**, *57*, 1923-1931.

(372) Meuwly, M.; Hutson, J. M. The potential energy surface and near-dissociation states of He-H₂⁺. *J. Chem. Phys.* **1999**, *110*, 3418-3427.

(373) Koner, D.; Bemish, R. J.; Meuwly, M. The C(³P) + NO(X²Π) → O(³P) + CN(X²Σ⁺), N(²D)/N(⁴S) + CO(X¹Σ⁺) reaction: rates, branching ratios, and final states from 15 to 20000 K. *J. Chem. Phys.* **2018**, *149*, 094305.

(374) San Vicente Veliz, J. C.; Koner, D.; Schwilk, M.; Bemish, R. J.; Meuwly, M. The N(⁴S) + O₂(X³Σ) ↔ O(³P) + NO(X²Π) reaction: thermal and vibrational relaxation rates for the ²A', ⁴A' and ²A'' states. *Phys. Chem. Chem. Phys.* **2020**, *22*, 3927-3939.

(375) Koner, D.; Meuwly, M. Permutationally invariant, reproducing kernel-based potential energy surfaces for polyatomic molecules: from formaldehyde to acetone. *J. Chem. Theory Comp.* **2020**, *16*, 5474-5484.

(376) Song, K.; Käser, S.; Töpfer, K.; Vazquez-Salazar, L. I.; Meuwly, M. PhysNet meets CHARMM: A framework for routine machine learning/molecular mechanics simulations. *J. Chem. Phys.* **2023**, *159*, 024125.

(377) Cui, Q.; Pal, T.; Xie, L. Biomolecular QM/MM simulations: what are some of the “burning issues”. *J. Phys. Chem. B* **2021**, *125*, 689-702.

(378) Gastegger, M.; Schütt, K. T.; Müller, K.-R. Machine learning of solvent effects on molecular spectra and reactions. *Chem. Sci.* **2021**, *12*, 11473-11483.

(379) Töpfer, K.; Käser, S.; Meuwly, M. Double proton transfer in hydrated formic acid dimer: interplay of spatial symmetry and solvent-generated force on reactivity. *Phys. Chem. Chem. Phys.* **2022**, *24*, 13869-13882.

(380) Lier, B.; Poliak, P.; Marquetand, P.; Westermayr, J.; Oostenbrink, C. BuRNN: buffer region neural network approach for polarizable-embedding neural network/molecular mechanics simulations. *J. Phys. Chem. Lett.* **2022**, *13*, 3812-3818.
(381) Smith, J. S.; Isayev, O.; Roitberg, A. E. ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost. *Chem. Sci.* **2017**, *8*, 3192–3203.

(382) Upadhyay, M.; Meuwly, M. Thermal and vibrationally activated decomposition of the Syn-CH₃CHOO Criegee intermediate. *ACS Earth Space Chem.* **2021**, *5*, 3396–3406.

(383) Hänninen, V.; Murdachaew, G.; Nathanson, G. M.; Gerber, R. B.; Halonen, L. Ab initio molecular dynamics studies of formic acid dimer colliding with liquid water. *Phys. Chem. Chem. Phys.* **2018**, *20*, 23717–23725.

(384) Schütt, K. T.; Gastegger, M.; Tkatchenko, A.; Müller, K. R.; Maurer, R. J. Unifying machine learning and quantum chemistry with a deep neural network for molecular wavefunctions. *Nat. Commun.* **2019**, *10*, 5024.

(385) Käser, S.; Richardson, J. O.; Meuwly, M. Transfer learning for affordable and high-quality tunneling splittings from instanton calculations. *J. Chem. Theory Comp.* **2022**, *18*, 6840–6850.

(386) Käser, S.; Meuwly, M. Transfer learned potential energy surfaces: accurate anharmonic vibrational dynamics and dissociation energies for the formic acid monomer and dimer. *Phys. Chem. Chem. Phys.* **2022**, *24*, 5269–5281.

(387) Kidwell, N. M.; Li, H.; Wang, X.; Bowman, J. M.; Lester, M. I. Unimolecular dissociation dynamics of vibrationally activated CH₃CHOO Criegee intermediates to OH radical products. *Nat. Chem.* **2016**, *8*, 509–514.

(388) Käser, S.; Unke, O. T.; Meuwly, M. Isomerization and decomposition reactions of acetaldehyde relevant to atmospheric processes from dynamics simulations on neural network-based potential energy surfaces. *J. Chem. Phys.* **2020**, *152*, 214304.

(389) Upadhyay, M.; Töpfer, K.; Meuwly, M. Molecular simulation for atmospheric reactions: non-equilibrium dynamics, roaming, and glycolaldehyde formation following photoinduced decomposition of syn-acetaldehyde oxide. *J. Phys. Chem. Lett.* **2024**, *15*, 90–96.

(390) Tan, M.-L.; Tran, K. N.; Pickard IV, F. C.; Simmonett, A. C.; Brooks, B. R.; Ichiye, T. Molecular multipole potential energy functions for water. *J. Phys. Chem. B* **2016**, *120*, 1833–1842.

(391) Ren, P.; Ponder, J. W. Polarizable atomic multipole water model for molecular mechanics simulation. *J. Phys. Chem. B* **2003**, *107*, 5933–5947.

(392) Huang, J.; Simmonett, A. C.; Pickard, F. C.; MacKerell, A. D., Jr.; Brooks, B. R. Mapping the Drude polarizable force field onto a multipole and induced dipole model. *J. Chem. Phys.* **2017**, *147*, 161702.

(393) Simmonett, A. C.; Pickard, F. C.; Shao, Y.; Cheatham, T. E.; Brooks, B. R. Efficient treatment of induced dipoles. *J. Chem. Phys.* **2015**, *143*, 074115.

(394) Simmonett, A. C.; Pickard, F. C.; Ponder, J. W.; Brooks, B. R. An empirical extrapolation scheme for efficient treatment of induced dipoles. *J. Chem. Phys.* **2016**, *145*, 164101.

(395) Gao, J.; Habibollazadeh, D.; Shao, L. A polarizable intermolecular potential function for simulation of liquid alcohols. *J. Phys. Chem.* **1995**, *99*, 16460–16467.

(396) Gao, J.; Pavelites, J. J.; Habibollazadeh, D. Simulation of liquid amides using a polarizable intermolecular potential function. *J. Phys. Chem.* **1996**, *100*, 2689–2697.

(397) Xie, W.; Pu, J.; MacKerell, A. D., Jr.; Gao, J. Development of a polarizable intermolecular potential function (PIPF) for liquid amides and alkanes. *J. Chem. Theory Comp.* **2007**, *3*, 1878–1889.

(398) Xie, W.; Pu, J.; Gao, J. A coupled polarization-matrix inversion and iteration approach for accelerating the dipole convergence in a polarizable potential function. *J. Phys. Chem. A* **2009**, *113*, 2109–2116.

(399) Patel, S.; Mackerell, A. D., Jr.; Brooks, C. L., III. CHARMM fluctuating charge force field for proteins: II protein/solvent properties from molecular dynamics simulations using a nonadditive electrostatic model. *J. Comput. Chem.* **2004**, *25*, 1504–1514.

(400) Patel, S.; Brooks, C. L., III. CHARMM fluctuating charge force field for proteins: I parametrization and application to bulk organic liquid simulations. *J. Comput. Chem.* **2004**, *25*, 1–16.

(401) Lamoureux, G.; MacKerell, A. D., Jr.; Roux, B. A simple polarizable model of water based on classical Drude oscillator. *J. Chem. Phys.* **2003**, *119*, 5185–5197.

(402) Lamoureux, G.; Roux, B. Modeling induced polarization with classical Drude oscillators: theory and molecular dynamics simulation algorithm. *J. Chem. Phys.* **2003**, *119*, 3025–3039.

(403) Angyan, J. G.; Colonna-Cesari, F.; Tapia, O. Analytical first and second energy derivatives in the polarization model. *Chem. Phys. Lett.* **1990**, *166*, 180–188.

(404) Thole, B. T. Molecular polarizabilities calculated with a modified dipole interaction. *Chem. Phys.* **1981**, *59*, 341–350.

(405) van Duijnen, P. T.; Swart, M. Molecular and atomic polarizabilities: Thole's model revisited. *J. Phys. Chem. A* **1998**, *102*, 2399–2407.

(406) Vesely, F. J. N-particle dynamics of polarizable Stockmayer-type molecules. *J. Chem. Phys.* **1977**, *24*, 361–371.

(407) Ding, Y.; Bernardo, D. N.; Krogh-Jespersen, K.; Levy, R. M. Solvation free energies of small amides and amines from molecular dynamics/free energy perturbation simulations using pairwise additive and many-body polarizable potentials. *J. Phys. Chem.* **1995**, *99*, 11575–11583.

(408) Van Belle, D.; Froeyen, M.; Lippens, G.; Wodak, S. J. Molecular dynamics simulation of polarizable water by an extended Lagrangian method. *Mol. Phys.* **1992**, *77*, 239–255.

(409) Hoover, W. Canonical dynamics: equilibrium phase-space distributions. *Phys. Rev. A* **1985**, *31*, 1695–1697.

(410) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. *J. Chem. Phys.* **1984**, *81*, 511–519.

(411) Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. *Mol. Phys.* **1984**, *52*, 255–268.

(412) Winkler, R. G.; Kraus, V.; Reineker, P. Time reversible and phase-space conserving molecular dynamics at constant temperature. *J. Chem. Phys.* **1995**, *102*, 9018–9025.

(413) Ahlstrom, P.; Wallqvist, A.; Engstrom, S.; Jonsson, B. A molecular dynamics study of polarizable water. *Mol. Phys.* **2006**, *68*, 563–581.

(414) Applequist, J.; Carl, J. R.; Fung, K.-K. An atom dipole interaction model for molecular polarizability. application to polyatomic molecules and determination of atom polarizabilities. *J. Am. Chem. Soc.* **1972**, *94*, 2952–2960.

(415) Kim, B.; Shao, Y.; Pu, J. Doubly polarized QM/MM with machine learning chaperone polarizability. *J. Chem. Theory Comp.* **2021**, *17*, 7682–7695.

(416) Feller, S. E.; Pastor, R. W.; Rojnuckarin, A.; Bogusz, S.; Brooks, B. R. Effect of electrostatic force truncation on interfacial and transport properties of water. *J. Phys. Chem.* **1996**, *100*, 17011–17020.

(417) Klauda, J. B.; Wu, X. W.; Pastor, R. W.; Brooks, B. R. Long-range Lennard-Jones and electrostatic interactions in interfaces: application of the isotropic periodic sum method. *J. Phys. Chem. B* **2007**, *111*, 4393–4400.

(418) Venable, R. M.; Chen, L. E.; Pastor, R. W. Comparison of the extended isotropic periodic sum and particle mesh Ewald methods for simulations of lipid bilayers and monolayers. *J. Phys. Chem. B* **2009**, *113*, 5855–5862.

(419) Wennberg, C. L.; Murtola, T.; Hess, B.; Lindahl, E. Lennard-Jones lattice summation in bilayer simulations has critical effects on surface tension and lipid properties. *J. Chem. Theory Comp.* **2013**, *9*, 3527–3537.

(420) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O'Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D., Jr.; Pastor, R. W. Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. *J. Phys. Chem. B* **2010**, *114*, 7830–7843.

(421) Wennberg, C. L.; Murtola, T.; Pall, S.; Abraham, M. J.; Hess, B.; Lindahl, E. Direct-space corrections enable fast and accurate Lorentz-Berthelot combination rule Lennard-Jones lattice summation. *J. Chem. Theory Comp.* **2015**, *11*, 5737–5746.

(422) Leonard, A. N.; Simmonett, A. C.; Pickard, F. C.; Huang, J.; Venable, R. M.; Klauda, J. B.; Brooks, B. R.; Pastor, R. W. Comparison
of additive and polarizable models with explicit treatment of long-range Lennard-Jones interactions using alkane simulations. *J. Chem. Theory Comp.* **2018**, *14*, 948–958.

(423) Perram, J. W.; Petersen, H. G.; De Leeuw, S. W. An algorithm for the simulation of condensed matter which grows as the 3/2 power of the number of particles. *Mol. Phys.* **1988**, *65*, 875–893.

(424) Wu, X. W.; Brooks, B. R. Isotropic periodic sum: a method for the calculation of long-range interactions. *J. Chem. Phys.* **2005**, *122*, 044107.

(425) Srinivasan, J.; Trevathan, M. W.; Beroza, P.; Case, D. A. Application of a pairwise generalized Born model to proteins and nucleic acids: inclusion of salt effects. *Theor. Chem. Acc.* **1999**, *101*, 426–434.

(426) Attanasio, F.; Convertino, M.; Magno, A.; Caflisch, A.; Corazza, A.; Haridas, H.; Esposito, G.; Cataldo, S.; Pignataro, B.; Milardi, D.; et al. Carnosine inhibits Aβ42 aggregation by perturbing the h-bond network in and around the central hydrophobic cluster. *ChemBioChem* **2013**, *14*, 583–592.

(427) Carballo-Pacheco, M.; Vancea, I.; Strodel, B. Extension of the FACTS implicit solvation model to membranes. *J. Chem. Theory Comp.* **2014**, *10*, 3163–3176.

(428) Zoete, V.; Schuepbach, T.; Bovigny, C.; Chaskar, P.; Daina, A.; Röhrig, U. F.; Michielin, O. Attracting cavities for docking. Replacing the rough energy landscape of the protein by a smooth attracting landscape. *J. Comput. Chem.* **2016**, *37*, 437–447.

(429) Röhrig, U. F.; Goullieux, M.; Bugnon, M.; Zoete, V. Attracting cavities 2.0: improving the flexibility and robustness for small-molecule docking. *J. Chem. Inf. and Mod.* **2023**, *63*, 3925–3940.

(430) Conti, S.; Ovchinnikov, V.; Karplus, M. ppdx: automated modeling of protein–protein interaction descriptors for use with machine learning. *J. Comput. Chem.* **2022**, *43*, 1747–1757.

(431) Lazaridis, T. Effective energy function for proteins in lipid membranes. *PROTEINS: Struct., Func. and Bioinf.* **2003**, *52*, 176–192.

(432) Lazaridis, T. Implicit solvent simulations of peptide interactions with anionic lipid membranes. *PROTEINS: Struct., Func. and Bioinf.* **2005**, *58*, 518–527.

(433) Lazaridis, T.; Karplus, M. Effective energy function for proteins in solution. *Proteins: Structure, Function, and Bioinformatics* **1999**, *35*, 133–152.

(434) Lazaridis, T. Structural determinants of transmembrane beta-barrels. *J. Chem. Theory Comp.* **2005**, *1*, 716–722.

(435) Mihajlovic, M.; Lazaridis, T. Antimicrobial peptides bind more strongly to membrane pores. *Biochim. et Biophys. Acta Biomembranes* **2010**, *1798*, 1494–1502.

(436) He, Y.; Prieto, L.; Lazaridis, T. Modeling peptide binding to anionic membrane pores. *J. Comput. Chem.* **2013**, *34*, 1463–1475.

(437) Lazaridis, T.; He, Y.; Prieto, L. Membrane interactions and pore formation by the antimicrobial peptide protegrin. *Biophys. J.* **2013**, *104*, 633–642.

(438) Prieto, L.; He, Y.; Lazaridis, T. Protein arcs may form stable pores in lipid membranes. *Biophys. J.* **2014**, *106*, 154–161.

(439) He, Y.; Lazaridis, T. Activity determinants of helical antimicrobial peptides: a large-scale computational study. *PLOS One* **2013**, *8*, No. e66440.

(440) Rahaman, A.; Lazaridis, T. Alamethicin pore formation: a thermodynamic approach. *Biochim. et Biophys. Acta Biomembranes* **2014**, *1838*, 1440–1447.

(441) Lipkin, R. B.; Lazaridis, T. Implicit membrane investigation of the stability of antimicrobial peptide β-barrels and arcs. *J. Membr. Biol.* **2015**, *248*, 469–486.

(442) Lipkin, R.; Lazaridis, T. Computational prediction of the optimal oligomeric state for membrane-inserted β-barrels of protegrin-1 and related mutants. *J. Pep. Sci.* **2017**, *23*, 334–345.

(443) Sepehri, A.; Nepal, B.; Lazaridis, T. Distinct modes of action of IAPP oligomers on membranes. *J. Chem. Inf. and Mod.* **2021**, *61*, 4645–4655.

(444) Sepehri, A.; Lazaridis, T. Putative structures of membrane-embedded amyloid β oligomers. *ACS Chem. Neurosci.* **2023**, *14*, 99–110.

(445) Dutta, A.; Sepehri, A.; Lazaridis, T. Putative pore structures of amyloid β 25–35 in lipid bilayers. *Biochemistry* **2023**, *62*, 2549–2558.

(446) Maurer, M.; Lazaridis, T. Transmembrane β-barrel models of α-synuclein oligomers. *J. Chem. Inf. and Mod.* **2023**, *63*, 7171–7179.

(447) Alford, R. F.; Fleming, P. J.; Fleming, K. G.; Gray, J. J. Protein structure prediction and design in a biologically realistic implicit membrane. *Biophys. J.* **2020**, *118*, 2042–2055.

(448) Nepal, B.; Leveritt, J.; Lazaridis, T. Membrane curvature sensing by amphipathic helices: insights from implicit membrane modeling. *Biophys. J.* **2018**, *114*, 2128–2141.

(449) Zhan, H.; Lazaridis, T. Inclusion of lateral pressure/curvature stress effects in implicit membrane models. *Biophys. J.* **2013**, *104*, 643–654.

(450) Nepal, B.; Sepehri, A.; Lazaridis, T. Mechanisms of negative membrane curvature sensing and generation by ESCRT III subunit Snf7. *Protein Sci.* **2020**, *29*, 1473–1485.

(451) Nepal, B.; Sepehri, A.; Lazaridis, T. Mechanism of negative membrane curvature generation by I-BAR domains. *Structure* **2021**, *29*, 1440–1452.

(452) Vasquez Rodriguez, S. Y.; Lazaridis, T. Simulations suggest a scaffolding mechanism of membrane deformation by the caveolin 8S complex. *Biophys. J.* **2023**, *122*, 4082–4090.

(453) Tanizaki, S.; Feig, M. A generalized Born formalism for heterogeneous dielectric environments: application to the implicit modeling of biological membranes. *J. Chem. Phys.* **2005**, *122*, 124706.

(454) Sayadi, M.; Tanizaki, S.; Feig, M. Effect of membrane thickness on conformational sampling of phospholamban from computer simulations. *Biophys. J.* **2010**, *98*, 805–814.

(455) Mirjalili, V.; Feig, M. Interactions of amino acid side-chain analogs within membrane environments. *J. Phys. Chem. B* **2015**, *119*, 2877–2885.

(456) Dutagaci, B.; Sayadi, M.; Feig, M. Heterogeneous dielectric generalized Born model with a van der Waals term provides improved association energetics of membrane-embedded transmembrane helices. *J. Comput. Chem.* **2017**, *38*, 1308–1320.

(457) Panahi, A.; Feig, M. Dynamic heterogeneous dielectric generalized Born (DHDGB): an implicit membrane model with a dynamically varying bilayer thickness. *J. Chem. Theory Comp.* **2013**, *9*, 1709–1719.

(458) Spassov, V. Z.; Yan, L.; Szalma, S. Introducing an implicit membrane in generalized Born/solvent accessibility continuum solvent models. *J. Phys. Chem. B* **2002**, *106*, 8726–8738.

(459) Stern, H. A.; Feller, S. E. Calculation of the dielectric permittivity profile for a nonuniform system: application to a lipid bilayer simulation. *J. Chem. Phys.* **2003**, *118*, 3401–3412.

(460) Feig, M.; Im, W.; Brooks, C. L., III. Implicit solvation based on generalized Born theory in different dielectric environments. *J. Chem. Phys.* **2004**, *120*, 903–911.

(461) Jaskierny, A. J.; Panahi, A.; Feig, M. Effect of flanking residues on the conformational sampling of the internal fusion peptide from Ebola virus. *PROTEINS: Struct., Func. and Bioinf.* **2011**, *79*, 1109–1117.

(462) Panahi, A.; Feig, M. Conformational sampling of influenza fusion peptide in membrane bilayers as a function of termini and protonation states. *J. Phys. Chem. B* **2010**, *114*, 1407–1416.

(463) Dutagaci, B.; Feig, M. Determination of hydrophobic lengths of membrane proteins with the hdgb implicit membrane model. *J. Chem. Inf. and Mod.* **2017**, *57*, 3032–3042.

(464) Brocke, S. A.; Degen, A.; MacKerell, A. D., Jr.; Dutagaci, B.; Feig, M. Prediction of membrane permeation of drug molecules by combining an implicit membrane model with machine learning. *J. Chem. Inf. and Mod.* **2019**, *59*, 1147–1162.

(465) Dutagaci, B.; Wittayanarakul, K.; Mori, T.; Feig, M. Discrimination of native-like states of membrane proteins with implicit membrane-based scoring functions. *J. Chem. Theory Comp.* **2017**, *13*, 3049–3059.

(466) Dutagaci, B.; Heo, L.; Feig, M. Structure refinement of membrane proteins via molecular dynamics simulations. *PROTEINS: Struct., Func. and Bioinf.* **2018**, *86*, 738–750.