CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed

Wonmuk Hwang, Steven L Austin, Arnaud Blondel, Eric D Boittier, Stefan Boresch, Matthias Buck, Joshua Buckner, Amedeo Caflisch, Hao-Ting Chang, Xi Cheng, et al.

► To cite this version:

Wonmuk Hwang, Steven L Austin, Arnaud Blondel, Eric D Boittier, Stefan Boresch, et al. CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed. Journal of Physical Chemistry B, 2024, 128 (41), pp.9976-10042. 10.1021/acs.jpcb.4c04100. pasteur-05110290

HAL Id: pasteur-05110290

https://pasteur.hal.science/pasteur-05110290v1

Submitted on 12 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed

Published as part of *The Journal of Physical Chemistry B* special issue "Recent Advances in Simulation Software and Force Fields".

Wonmuk Hwang,* Steven L. Austin, Arnaud Blondel, Eric D. Boittier, Stefan Boresch, Matthias Buck, Joshua Buckner, Amedeo Caflisch, Hao-Ting Chang, Xi Cheng, Yeol Kyo Choi, Jhih-Wei Chu, Michael F. Crowley, Qiang Cui, Ana Damjanovic, Yuqing Deng, Mike Devereux, Xinqiang Ding, Michael F. Feig, Jiali Gao, David R. Glowacki, James E. Gonzales, II, Mehdi Bagerhi Hamaneh, Edward D. Harder, Ryan L. Hayes, Jing Huang, Yandong Huang, Phillip S. Hudson, Wonpil Im, Shahidul M. Islam, Wei Jiang, Michael R. Jones, Silvan Käser, Fiona L. Kearns, Nathan R. Kern, Jeffery B. Klauda, Themis Lazaridis, Jinhyuk Lee, Justin A. Lemkul, Xiaorong Liu, Yun Luo, Alexander D. MacKerell, Jr., Dan T. Major, Markus Meuwly, Kwangho Nam, Lennart Nilsson, Victor Ovchinnikov, Emanuele Paci, Soohyung Park, Richard W. Pastor, Amanda R. Pittman, Carol Beth Post, Samarjeet Prasad, Jingzhi Pu, Yifei Qi, Thenmalarchelvi Rathinavelan, Daniel R. Roe, Benoit Roux, Christopher N. Rowley, Jana Shen, Andrew C. Simmonett, Alexander J. Sodt, Kai Töpfer, Meenu Upadhyay, Arjan van der Vaart, Luis Itza Vazquez-Salazar, Richard M. Venable, Luke C. Warrensford, H. Lee Woodcock, Yujin Wu, Charles L. Brooks, III,* Bernard R. Brooks,* and Martin Karplus

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

Cite This: *J. Phys. Chem. B* 2024, 128, 9976–10042

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

Read Online

## ACCESS

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

Metrics & More

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

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

Article Recommendations

**ABSTRACT:** Since its inception nearly a half century ago, CHARMM has been playing a central role in computational biochemistry and biophysics. Commensurate with the developments in experimental research and advances in computer hardware, the range of methods and applicability of CHARMM have also grown. This review summarizes major developments that occurred after 2009 when the last review of CHARMM was published. They include the following: new faster simulation engines, accessible user interfaces for convenient workflows, and a vast array of simulation and analysis methods that encompass quantum mechanical, atomistic, and coarse-grained levels, as well as extensive coverage of force fields. In addition to providing the current snapshot of the CHARMM development, this review may

serve as a starting point for exploring relevant theories and computational methods for tackling contemporary and emerging problems in biomolecular systems. CHARMM is freely available for academic and nonprofit research at https://academiccharmm.org/program.

## CHARMM

Force Field
System Preparation
Structure Manipulation

API
Accelerator
MD Protocol
QM/MM
Analysis

Constant pH
Docking
Free Energy
Sampling

1979–1983

2009

2024

## CONTENTS

1. Introduction 9978

2. CHARMM Accelerator Engines 9978

2.1. The OpenMM API 9979

2.2. DOMDEC Parallel-Scalable Platform 9979

2.3. BLaDE 9979

Received: June 20, 2024

Revised: August 15, 2024

Accepted: August 22, 2024

Published: September 20, 2024

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

<figure><img src="image_9.png" /></figure>
<table>
  <tbody>
    <tr>
      <td>2.4. apoCHARMM, Embracing CHARMM-Centric Functionality</td>
      <td>9980</td>
      <td>9.2. SSNMR Module</td>
      <td>9999</td>
    </tr>
    <tr>
      <td>3. Streamlining CHARMM Workflows</td>
      <td>9980</td>
      <td>9.3. Residual Dipolar Coupling (RDC) NMR Orientational Restraint</td>
      <td>10000</td>
    </tr>
    <tr>
      <td>3.1. pyCHARMM</td>
      <td>9980</td>
      <td>9.4. Torque Application</td>
      <td>10001</td>
    </tr>
    <tr>
      <td>3.2. crimm</td>
      <td>9981</td>
      <td>10. CHARMM Force Field Development</td>
      <td>10001</td>
    </tr>
    <tr>
      <td>3.3. CHARMM-GUI</td>
      <td>9981</td>
      <td>10.1. Water, Ions, and Polar Solvents</td>
      <td>10001</td>
    </tr>
    <tr>
      <td>4. Docking Methods</td>
      <td>9982</td>
      <td>10.2. Proteins</td>
      <td>10001</td>
    </tr>
    <tr>
      <td>4.1. CDOCKER</td>
      <td>9982</td>
      <td>10.3. Nucleic Acids</td>
      <td>10002</td>
    </tr>
    <tr>
      <td>4.2. EnzyDock</td>
      <td>9982</td>
      <td>10.4. Lipids</td>
      <td>10002</td>
    </tr>
    <tr>
      <td>4.3. CIFDock</td>
      <td>9983</td>
      <td>10.5. Carbohydrates</td>
      <td>10003</td>
    </tr>
    <tr>
      <td>5. Free Energy Methods</td>
      <td>9983</td>
      <td>10.6. Small Molecules</td>
      <td>10003</td>
    </tr>
    <tr>
      <td>5.1. λ-Dynamics, Multisite λ-Dynamics, and Constant pH MD</td>
      <td>9983</td>
      <td>11. Mixed Quantum Mechanics/Molecular Mechanics (QM/MM) Methods</td>
      <td>10004</td>
    </tr>
    <tr>
      <td>5.2. Hybrid Sampling and Free Energy Algorithms</td>
      <td>9985</td>
      <td>11.1. Background</td>
      <td>10004</td>
    </tr>
    <tr>
      <td>5.3. Optimal Variance Alchemical Path for Free Energy Calculation</td>
      <td>9985</td>
      <td>11.2. QM/MM Potentials and Practical Considerations</td>
      <td>10004</td>
    </tr>
    <tr>
      <td>6. Constant pH Methods</td>
      <td>9986</td>
      <td>11.3. Recent Advances in SE-QM/MM Methods</td>
      <td>10006</td>
    </tr>
    <tr>
      <td>6.1. Hybrid-Solvent and All-Atom Continuous Constant pH Methods</td>
      <td>9986</td>
      <td>11.3.1. Extended Lagrangian MD (ELMD).</td>
      <td>10006</td>
    </tr>
    <tr>
      <td>6.2. Constant pH MD with Discrete Protonation States</td>
      <td>9987</td>
      <td>11.3.2. Extended Lagrangian Born-Oppenheimer MD with Dissipation (DXL-BOMD).</td>
      <td>10006</td>
    </tr>
    <tr>
      <td>6.3. Proton Hopping Simulations</td>
      <td>9987</td>
      <td>11.3.3. Fock Matrix Dynamics (FMD).</td>
      <td>10006</td>
    </tr>
    <tr>
      <td>7. Enhanced Sampling and Transition Path Methods</td>
      <td>9988</td>
      <td>11.4. Path-Integral-Free Energy Perturbation for Nuclear Quantum Effects</td>
      <td>10007</td>
    </tr>
    <tr>
      <td>7.1. Replica Exchange MD (REMD)</td>
      <td>9988</td>
      <td>11.5. Density Functional Tight Binding (DFTB) Module</td>
      <td>10008</td>
    </tr>
    <tr>
      <td>7.2. Biasing Methods</td>
      <td>9988</td>
      <td>11.6. Double Link Atom Method (DLAM)</td>
      <td>10009</td>
    </tr>
    <tr>
      <td>7.2.1. Targeted MD (TMD)</td>
      <td>9988</td>
      <td>11.7. Flexible Inner Region Ensemble Separator (FIRES)</td>
      <td>10010</td>
    </tr>
    <tr>
      <td>7.2.2. Related Conformational Free Energy Sampling</td>
      <td>9989</td>
      <td>11.8. Combining QM/MM with Gaussian Process Machine Learning Potentials</td>
      <td>10010</td>
    </tr>
    <tr>
      <td>7.3. String Method (SM) for Conformational Transitions</td>
      <td>9989</td>
      <td>11.9. Multistate Empirical Valence Bond (MS-EVB)</td>
      <td>10010</td>
    </tr>
    <tr>
      <td>7.3.1. Zero-Temperature SM</td>
      <td>9989</td>
      <td>11.10. Reactive MD</td>
      <td>10011</td>
    </tr>
    <tr>
      <td>7.3.2. Finite-Temperature SM</td>
      <td>9989</td>
      <td>11.11. Indirect QM/MM Free Energy Simulations</td>
      <td>10011</td>
    </tr>
    <tr>
      <td>7.3.3. String in Collective Variables</td>
      <td>9990</td>
      <td>12. Boundary Condition, System Preparation, and Trajectory Analysis</td>
      <td>10012</td>
    </tr>
    <tr>
      <td>7.3.4. Reparameterization</td>
      <td>9990</td>
      <td>12.1. Simulation and Analysis of Membrane Proteins</td>
      <td>10012</td>
    </tr>
    <tr>
      <td>7.3.5. String with Swarms-of-Trajectories</td>
      <td>9990</td>
      <td>12.2. Coordinate Unwrapping for Diffusion Constants</td>
      <td>10012</td>
    </tr>
    <tr>
      <td>7.3.6. Script-Based SM Approach and Structure Building</td>
      <td>9990</td>
      <td>12.3. Calculation of Pressure Profiles in Lipid Bilayers</td>
      <td>10014</td>
    </tr>
    <tr>
      <td>7.4. Adaptively Biased Path Optimization (ABPO) for Transition Path Sampling</td>
      <td>9991</td>
      <td>12.4. P2<sub>1</sub> Periodic Boundary Condition</td>
      <td>10015</td>
    </tr>
    <tr>
      <td>7.5. Reaction Path Optimization with Holonomic Constraints</td>
      <td>9991</td>
      <td>12.5. Primary Hydration Shell (PHS) Model</td>
      <td>10016</td>
    </tr>
    <tr>
      <td>7.6. Boxed MD (BXD)</td>
      <td>9992</td>
      <td>12.6. Hydration Map</td>
      <td>10016</td>
    </tr>
    <tr>
      <td>7.7. Extended Adaptive Biasing Force (eABF) Method</td>
      <td>9992</td>
      <td>12.7. Conformational Entropy</td>
      <td>10017</td>
    </tr>
    <tr>
      <td>8. Advanced Energy Functions, Coarse Graining, and Implicit Models</td>
      <td>9993</td>
      <td>12.8. Identifying Non-Polar Contacts</td>
      <td>10017</td>
    </tr>
    <tr>
      <td>8.1. Multipolar Electrostatics</td>
      <td>9993</td>
      <td>12.9. Vectorial Analysis of Long-Range Concerted Motions in MD Trajectories</td>
      <td>10017</td>
    </tr>
    <tr>
      <td>8.2. Machine-Learning-Based Energy Functions</td>
      <td>9994</td>
      <td>12.10. Other Updated Preparation and Analysis Features</td>
      <td>10017</td>
    </tr>
    <tr>
      <td>8.3. Multipole and Point-Induced Dipole (MPID)</td>
      <td>9994</td>
      <td>12.10.1. System Generation</td>
      <td>10017</td>
    </tr>
    <tr>
      <td>8.4. Polarizable Intermolecular Potential Functions (PIPFs)</td>
      <td>9994</td>
      <td>12.10.2. Trajectory Handling</td>
      <td>10018</td>
    </tr>
    <tr>
      <td>8.5. Long-Range Lennard-Jones Interactions</td>
      <td>9995</td>
      <td>12.10.3. Time Series Analysis in the CORREL Module</td>
      <td>10018</td>
    </tr>
    <tr>
      <td>8.6. FACTS Implicit Solvent</td>
      <td>9996</td>
      <td>12.10.4. Similarity Analysis of Snapshots from Trajectory Files</td>
      <td>10018</td>
    </tr>
    <tr>
      <td>8.7. Implicit Modeling of Membranes</td>
      <td>9997</td>
      <td>13. Concluding Discussion</td>
      <td>10018</td>
    </tr>
    <tr>
      <td>8.7.1. Membrane Pores</td>
      <td>9997</td>
      <td>Author Information</td>
      <td>10018</td>
    </tr>
    <tr>
      <td>8.7.2. Curved Membranes</td>
      <td>9997</td>
      <td>Corresponding Authors</td>
      <td>10018</td>
    </tr>
    <tr>
      <td>8.7.3. Mean-Field Modeling of Deformable Membrane Bilayers via HDGB</td>
      <td>9997</td>
      <td>Authors</td>
      <td>10018</td>
    </tr>
    <tr>
      <td>8.8. Transferable Coarse-Graining via PRIMO</td>
      <td>9998</td>
      <td>Author Contributions</td>
      <td>10021</td>
    </tr>
    <tr>
      <td>9. Specialized Restraint Methods</td>
      <td>9999</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>9.1. CONSHELIX Module</td>
      <td>9999</td>
      <td></td>
      <td></td>
    </tr>
  </tbody>
</table>
<table><tr><td>Notes</td><td>10021</td></tr><tr><td>Acknowledgments</td><td>10021</td></tr><tr><td>References</td><td>10021</td></tr></table>

# 1. INTRODUCTION

CHARMM, the program for simulation and modeling of biomolecular systems, is now more than 40 years in continuous development and use. It is freely available to academic and government laboratory users, and is fast, as we describe below due to recent advances in GPU acceleration. From its earliest incarnations as Gandalf, renamed to HARMM (Harvard Macromolecular Mechanics), and finally CHARMM (Chemistry at Harvard Macromolecular Mechanics), it has provided a working framework for the exploration of biomolecular structure–function–dynamics relationships, beginning with the first molecular dynamics (MD) simulation of the small protein pancreatic trypsin inhibitor.¹ For an interesting perspective on the evolution of simulation methods and the applications of statistical mechanics to the study of biological molecules, we point readers to a recent review.²

CHARMM provides many practical and functional features that distinguish it from other programs that have evolved for similar purposes. Probably the most significant is that CHARMM was created as, and has remained, a repository for many of the trend-setting methods and models. It now comprises over 1,170,000 lines of code (modular Fortran 90, C, C++, CUDA, OpenCL, and Python) encapsulating its extensive functionality. An equally important and distinctive feature that has been integral to the software since its earliest days, is an interpreted language as its command parser, enabling “programs” to be written in “CHARMM-language” in contrast to other programs in this area whereby one prepares an input “script” describing the particular (one-pass) calculation one wishes to run. This feature greatly facilitates testing and prototyping of many of the statistical mechanical methods and techniques that have been integrated into CHARMM. While much of this has been described in earlier works and will not be further elaborated here, it is important to recognize this fundamental differentiator of CHARMM and other programs utilized in the field.

In this review, we provide an update to the developments that have occurred in CHARMM since it was last described in 2009.³ We refer readers to the previous two papers describing CHARMM for basic organization of the program and other functionality incorporated prior to 2009.⁵,⁴ Among the key new developments during the past 15 years, the most significant are the establishment of GPU-accelerated kernels to perform many of the unique calculations available in CHARMM. We describe the three published CHARMM accelerator engines in Sections 2.1 (CHARMM/OpenMM API), 2.2 (CHARMM/DOMDEC), 2.3 (CHARMM/BLaDE), and 2.4 (apoCHARMM). Each of the former three platforms are fully integrated with CHARMM and support a significant range of CHARMM functionality, and thereby provide powerful platforms for establishing complex simulation workflows utilizing CHARMM scripting language. apoCHARMM is a new GPU accelerator that is currently being developed. Except for DOMDEC, these engines provide performance comparable to any existing GPU-based biomolecular simulation code. Thus, CHARMM performance is on par with other codes while the accessibility of methods is typically richer.

Beyond the advanced simulation engines, significant efforts to integrate CHARMM into modern workflows have advanced, as described in Sections 3.1 (pyCHARMM), 3.2 (crimm), and 3.3 (CHARMM-GUI). CHARMM-GUI continues to mature and is a vital service to the community through its system and simulation setup facilities, providing input scripts for simulations for a number of current biomolecular MD simulation packages, including CHARMM, OpenMM, Amber, and GROMACS (Section 3.3). Section 3.1 on pyCHARMM describes recent efforts to use CHARMM functionality within the context of Python language. This enables the straightforward integration of CHARMM and the immense base of Python modules. Finally, crimm is a Python-based package that integrates many simulation preparation tasks, e.g., building of missing residues and loops, choosing protonation states appropriate for a given pH, patching to represent disulfide and other system modifications and solvation. These tools and methods of CHARMM provide an essential platform for modern simulation and modeling workflows.

New docking methods and procedures are described in Section 4, followed by an update of the free energy (Section 5) and constant pH (Section 6) methods. We then discuss the new enhanced sampling and transition path methods (Section 7). CHARMM has supported a range of reactive, implicit and coarse-grained models for simulation of biorelated systems throughout its history and new developments are discussed in Section 8. System-specific and specialized restraint methods that have been developed and implemented in CHARMM recently are in Section 9.

CHARMM’s force fields (FFs) are corner posts of molecular simulations throughout the field. Section 10 gives updates to the CHARMM fixed-charge and polarizable FFs. Following is a discussion of the recent approaches and methods for quantum mechanics and molecular mechanics (QM/MM) simulations that are integrated into CHARMM (see Section 11). Lastly, we describe new methods, procedures and analysis tools that have been integrated into CHARMM in Section 12.

Through this review, we hope to convey the immense base of methods and models that are available and supported in CHARMM. Additionally, we hope that readers will use this review as an entry to the growing online repositories for tutorials, advanced simulation methods, and templates.

# 2. CHARMM ACCELERATOR ENGINES

Since the previous review of CHARMM functionality,³ there has been a significant growth in utilizing new generation of processors, especially GPUs, to accelerate MD and related modeling tasks within biomolecular simulation methods. The CHARMM development community has embraced this effort by developing specialized GPU accelerator and highly parallel kernels, or adaptor APIs that support a range of the extensive functionality available in CHARMM, from free energy methods and constant pH simulation techniques, to implicit solvent models such as generalized Born (GB) models and a host of others. These interfaces provide a straightforward means to set up systems, manipulate, patch or otherwise prepare them for simulation and then simulate using GPU or parallel CPU execution without creating extraneous intermediate files. This seamless interface also enables straightforward analysis and visualization of results from within the same pyCHARMM/CHARMM script, thereby providing an integrated framework for modeling, dynamics, and analysis. We
<table><caption>Table 1. CHARMM Accelerated Engine Benchmarks on GPUs (ns/day; avg ± sd over 5 runs)ᵃ</caption><thead><tr><th></th><th>System</th><th>DHFR</th><th>APOA1</th><th>DMPG</th><th>T4L (N<sub>s</sub> = 6)</th><th>HSP90 (N<sub>s</sub> = 11)</th></tr></thead><tbody><tr><td>GPU</td><td>Atoms</td><td>23558</td><td>92224</td><td>291168</td><td>37723</td><td>26685</td></tr><tr><td rowspan="3">Ada RTX6000</td><td>BLaDE</td><td>601.6 ± 2.4</td><td>255.6 ± 1.2</td><td>82.1 ± 0.4</td><td>-</td><td>-</td></tr><tr><td>OpenMM</td><td>1052.0 ± 3.5</td><td>257.8 ± 2.8</td><td>78.0 ± 0.9</td><td>-</td><td>-</td></tr><tr><td>apoCHARMM</td><td>423.5 ± 1.3</td><td>289.2 ± 1.4</td><td>95.3 ± 2.5</td><td>-</td><td>-</td></tr><tr><td rowspan="3">Ada RTX5000</td><td>BLaDE</td><td>514.3 ± 6.0</td><td>191.7 ± 0.6</td><td>63.6 ± 0.2</td><td>-</td><td>-</td></tr><tr><td>OpenMM</td><td>897.3 ± 1.1</td><td>197.7 ± 2.3</td><td>62.7 ± 1.3</td><td>-</td><td>-</td></tr><tr><td>apoCHARMM</td><td>323.6 ± 11.5</td><td>218.3 ± 0.7</td><td>72.3 ± 0.3</td><td>-</td><td>-</td></tr><tr><td rowspan="3">Ada RTX4500</td><td>BLaDE</td><td>436.7 ± 0.45</td><td>142.9 ± 0.2</td><td>45.6 ± 0.04</td><td>-</td><td>-</td></tr><tr><td>OpenMM</td><td>730.0 ± 0.3</td><td>138.9 ± 1.2</td><td>43.6 ± 0.6</td><td>-</td><td>-</td></tr><tr><td>apoCHARMM</td><td>388.6 ± 0.7</td><td>152.2 ± 0.7</td><td>51.3 ± 0.2</td><td>-</td><td>-</td></tr><tr><td rowspan="3">RTX A6000</td><td>BLaDE</td><td>355.9 ± 2.7</td><td>124.5 ± 0.2</td><td>40.3 ± 0.4</td><td>-</td><td>-</td></tr><tr><td>OpenMM</td><td>656.0 ± 6.5</td><td>129.8 ± 1.6</td><td>39.5 ± 0.5</td><td>-</td><td>-</td></tr><tr><td>apoCHARMM</td><td>295.8 ± 0.7</td><td>133.8 ± 1.1</td><td>42.6 ± 0.3</td><td>-</td><td>-</td></tr><tr><td rowspan="3">RTX A5500</td><td>BLaDE</td><td>333.5 ± 1.0</td><td>111.8 ± 4.7</td><td>37.9 ± 0.3</td><td>-</td><td>-</td></tr><tr><td>OpenMM</td><td>594.3 ± 2.7</td><td>113.4 ± 0.6</td><td>31.2 ± 0.1</td><td>-</td><td>-</td></tr><tr><td>apoCHARMM</td><td>277.2 ± 5.6</td><td>135.3 ± 2.4</td><td>38.5 ± 0.2</td><td>-</td><td>-</td></tr><tr><td rowspan="4">A100</td><td>BLaDE</td><td>402.9 ± 4.7</td><td>149.4 ± 0.9</td><td>54.0 ± 0.2</td><td>208.8 ± 0.3</td><td>246.8 ± 0.4</td></tr><tr><td>DOMDEC</td><td>106.3 ± 0.7</td><td>-</td><td>9.07 ± 0.06</td><td>33.5 ± 0.0</td><td>46.4 ± 0.3</td></tr><tr><td>OpenMM</td><td>558.5 ± 8.6</td><td>121.7 ± 2.3</td><td>36.2 ± 0.6</td><td>-</td><td>-</td></tr><tr><td>apoCHARMM</td><td>276.1 ± 1.5</td><td>172.6 ± 1.3</td><td>63.6 ± 0.3</td><td>-</td><td>-</td></tr></tbody></table>

ᵃDHFR, APOA1, and DMPG are for NVT simulations with 9-Å real-space non-bonded cutoff distance. T4L and HSP90 are for λ-dynamics in NPT ensemble and with 10-Å real-space non-bonded cutoff distance. N<sub>s</sub> in parentheses is the number of λ variables used (cf., Eq. 4). PME was used to account for long-range electrostatic interactions.²⁴ In all simulations, the integration time step was 2 fs.

present below in chronological order, the accelerated performance platforms that are integrated into CHARMM. As indicated in Table 1, the performance of CHARMM through the GPU-accelerated APIs compares well with those observed in other GPU-accelerated packages such as pmemd.cuda,⁵ GROMACS,⁶ and NAMD.⁷ CHARMM/pyCHARMM (Section 3.1) also offers a range of developed and developing interfaces to meet needs across a breadth of the methodological application areas.

**2.1. The OpenMM API.** The first GPU-accelerated engine coupled to CHARMM comprised a FORTRAN-90 API that takes advantage of the significant developments of OpenMM.⁸ This interface provides direct calls to OpenMM functionality for MD, energy minimization, and free energy methods. A host of restraints existing in CHARMM are also implemented using OpenMM's custom forces.⁸ This effort, spearheaded by Michael Garrahan and Charles Brooks, first appeared in CHARMM release c37b1. A key advantage of using OpenMM through the CHARMM/OpenMM API is that applications that need to move between system setup and preparation, processing, and analysis can occur within a single workflow using the CHARMM interpreted command language or directly through Python with pyCHARMM.⁹ Aside from the many simulation environments, restraints, heavy atom-hydrogen constraints, NVT and NPT (using either an isotropic Monte-Carlo barostat or anisotropic pressure coupling often used in membrane simulations), particle-mesh Ewald (PME), CHARMM shifting and switching methods for van der Waals and/or electrostatic interactions, this interface also supports free energy perturbation (FEP) methods utilizing fixed windowing approaches, e.g., FEP with analysis by MBAR (Section 11.11).¹⁰⁻¹² The CHARMM/OpenMM API provides a robust platform for integrated modeling tasks and workflows and has been utilized in the extension of the CDOCKER approaches¹³⁻¹⁵ to parallel simulated annealing with GPU acceleration with any of the CHARMM-compatible physically

based FFs, including CHARMM,¹⁶ CHARMM General FF (CGenFF),¹⁷,¹⁸ AMBER,¹⁹ GAFF,²⁰ OPLS,²¹ and the LigParGen OPLS extension for small molecules.²² At present, not all of the CHARMM functionalities are fully implemented through the API, and a key missing element is full access to the Drude polarizable FF that has been rapidly developing over the past several years.²³ Plans are underway to provide availability through the CHARMM/OpenMM API in the near future.

**2.2. DOMDEC Parallel-Scalable Platform.** In 2014, a new DOMain DEComposition (DOMDEC) MD engine was introduced into CHARMM by Hynninen and Crowley.²⁵ It was faster both in the execution on serial and parallel CPU platforms. Serial performance was approximately two times higher than in the previous versions of CHARMM with its “fast” CPU-based options. The parallel version enabled efficient utilization up to hundreds of CPU cores.

The DOMDEC module of CHARMM served as an early platform for the development of the multisite λ-dynamics (MSλD)²⁶⁻²⁸ and explicit-solvent constant pH MD (CpHMD)²⁹⁻³¹ methods. Also implemented as part of this effort were GPU-resident kernels that accelerated components of the computation and enabled partitioning of the system being studied across both CPU and GPU cores for scaling and acceleration. Finally, a GPU implementation that handles all energy calculations except for the SHAKE constraint and position propagation was implemented in the “GPU only” functionality of DOMDEC. All of CHARMM’s MSλD and CpHMD functions were integrated into the “GPU-only” kernel, which was the platform on which NAMD’s GPU-accelerated kernel (gpu-offload) was based as well as more recent faster engines including the newly developing apoCHARMM (Section 2.4).

**2.3. BLaDE.** The BLaDE (Basic LAmbda Dynamics Engine) module of CHARMM²⁴ was developed to optimize the speed of λ-dynamics simulations on GPUs, but it also provides a robust and accelerated platform for conventional
MD simulations. Previously, the DOMDEC module²⁵ was the fastest implementation of λ-dynamics as noted above, but the SHAKE constraint and position propagation being handled by the CPU were rate limiting. Also, DOMDEC performed suboptimally on smaller systems on a single GPU. BLaDE optimizes these tasks and achieves 5- to 6-fold speedup over DOMDEC. Although the CHARMM/OpenMM API (discussed above) exhibits similar performance on standard MD (Table 1), it is less suited for λ-dynamics.

Table 1 shows benchmarks for BLaDE, DOMDEC, and OpenMM through their CHARMM interfaces, and a stand-alone version of apoCHARMM. We have focused on previously established benchmarks as presented earlier²⁴ and elsewhere. DHFR is a small globular protein, ApoA1 is a solvated lipid nanodisc, DMPG is a larger lipid bilayer. T4L is a protein mutation calculation, and HSP90 is a ligand perturbation calculation. Benchmarks were repeated 5 times and run on single NVIDIA GPUs as noted in Table 1. These benchmarks demonstrate that BLaDE scales well, especially for λ-dynamics.

Since BLaDE is designed to be simple and fast, not all features present in CHARMM are available in BLaDE. For example, while much of CHARMM uses the Langevin Piston barostat,³² BLaDE uses the Monte Carlo barostat³³,³⁴ in constant pressure simulations, removing the overhead of computing the virial and rectifying SHAKE constraints after coordinate update. Similarly, BLaDE only includes a Langevin thermostat, and is not yet implemented in energy minimization routines, though energy calls can be made directly to BLaDE. New features continue to be added, including support for harmonic and nuclear Overhauser effect (NOE) restraints, support for non-orthogonal boxes, PME³⁵ or force/energy switching electrostatics,³⁶ and support for AMBER FF with different 1–4 scaling and improper torsion potentials.

## 2.4. apoCHARMM, Embracing CHARMM-Centric Functionality.
apoCHARMM is a developing open source package designed specifically to support some of the distinctive methods of CHARMM absent in the CHARMM/OpenMM or CHARMM/BLaDE APIs (see Sections 2.1 and 2.3), and at the speeds provided by modern GPU architectures (Table 1). In particular, apoCHARMM in its current development supports, or plans to support:

* A complete analytic virial tensor.

* Multiple PSFs (protein structure files) simultaneously (upper limit set by the hardware resource).

* Uncommon crystal symmetries such as $P2_1$ (Section 12.4).

By accounting for the complete virial tensor, an implementation of the Langevin piston algorithm³² for constant pressure or constant surface tension ensembles is enabled. Simultaneous support for multiple PSFs allows free energy methods modeled after the CHARMM PERT approach to be run. It also allows the enveloping distribution sampling (EDS) based method for free energies,³⁷⁻³⁹ and for state-based CpHMD.⁴⁰,⁴¹ Support for $P2_1$ crystal symmetry⁴²,⁴³ allows lipid bilayer systems to be simulated without chemical potential mismatch between upper and lower leaflets, which can be very useful when making membrane insertions⁴⁴ (Section 12.4).

A number of different integrators have been implemented for different ensembles: Velocity-Verlet and leapfrog integrators for the microcanonical ensemble, and Langevin thermostat

and Nosé-Hoover integrators for the canonical ensemble. Holonomic constraints are handled using SHAKE and SETTLE algorithms. Since the virial is calculated during force calculation, the isobaric ensemble can be sampled using the Langevin piston method.³² This is an extended ensemble method with additional degrees of freedom corresponding to pistons that are used to control the pressure. Thus, a number of ensembles are available including constant area (NPAT) and constant surface tension (NPγT).

Several methods for free energy difference calculations are implemented in apoCHARMM.⁴³ A unifying scheme, which depends on a variant of energy interpolation, is implemented using a composite design pattern, where forces and energies of the end states after being separately calculated are interpolated. Additionally, soft-core formulation of the van der Waals interaction⁴⁶ to calculate λ-specific energy is available. The double exponential method⁴⁷ has been implemented as well. Although it is slightly slower than the van der Waals formulation, it provides a base version of soft-core.

apoCHARMM is derived from the erstwhile GitHub package by Antti-Pekka Hynninen.⁴⁸ It is written in CUDA and modern C++ to leverage the full potential of NVIDIA GPU architectures. Additionally, it features a Pybind11-based Python interface, ensuring convenience for end-users. The codebase adheres to test driven development (TDD) principles and incorporates Catch2-based unit tests with extensive code coverage. Notably, apoCHARMM is designed as a GPU-exclusive implementation, with all aspects of MD including energy and force calculations, restraints, constraints, and integration, executed entirely on GPU. Minimizing host-GPU memory transfers, the system only necessitates such transfers during logging or trajectory saving operations. One of the modular design patterns employed is the mediator pattern that reduces dependencies between different components by mandating communication through a central mediator object. In a similar vein, loggers and integrators leverage a publisher-subscriber design pattern, facilitating the versatile reuse of different loggers with distinct integrators. Overall, apoCHARMM performance is comparable to or better than other GPU based MD engines (Table 1). Since it is optimized for larger systems, its performance is not as good for the smaller DHFR.

## 3. STREAMLINING CHARMM WORKFLOWS

The rich methodology and broad functionality of CHARMM, including its unique scripting language, has enabled many complex workflows to be created and tested prior to committing them to code in FORTRAN 90/C/C++.³ These scripting capabilities have led to extensive libraries of CHARMM scripts in various forums and repositories as well as seeded the establishment of the web-based CHARMM-GUI⁴⁹ and a range of other modeling tasks, e.g., MCSS,⁵⁰ early stages of XPLOR,⁵¹ and SILCS.⁵² CHARMM scripting language, although extremely powerful, is not naturally integrable with other workflows that are convenient and widely used in the modeling of biomolecules. This realization has provided impetus for establishing a complete Python interface to CHARMM’s full range of functionality and efforts to facilitate the utilization of CHARMM and pyCHARMM, namely crimm and CHARMM-GUI.

### 3.1. pyCHARMM.
Efforts were initiated in the group of Charles Brooks to develop CHARMM callable functionality though a Python interface and APIs, called pyCHARMM.⁹ It
enables CHARMM variables and data structures to be explored and used, and in some instances manipulated at the Python level, providing the means of creating complex workflows that integrate and extend tools built in Python for numerical and graphical tasks. Native Python functions and modules complement and extend the already rich landscape of CHARMM functionalities. Examples include a framework that enables novel energy functions to be integrated with CHARMM's modeling tools through Python callable routines available in Python, CUDA, and OpenCL, as well as utilize machine learned functions such as TORCH-ANI⁵³ and PhysNet⁵⁴ for energy and force calculation. Analogous 'hooks' are built into the CHARMM dynamics engine. Graphical engines are also readily integrated into pyCHARMM for rapid visualization of simulation models and results. Loosely coupling tasks across many processors too is straightforward within pyCHARMM workflows using MPI frameworks such as MPI4PY and this has facilitated free energy calculations using multiscale Bennett's acceptance ratio (MBAR) and thermodynamic integration (TI) approaches, or high-throughput MS/MS free energy methods (Section 5.1), string path optimization calculations (Section 7.3), replica exchange (Section 7.1), and fully automated docking workflows employing CDOCKER (Section 4.1).¹³,¹⁴,⁵⁵,⁵⁶ pyCHARMM is integrated with the accelerated platform kernels and APIs of CHARMM/OpenMM (Section 2.1) and CHARMM/BLADE (Section 2.3).

Since its release in early 2023,⁹ two workshops have been given on integrating modeling tasks using pyCHARMM. The first focus was on general modeling tasks and methods with examples provided as Jupyter Notebooks and Python scripts (July 2022),⁵⁷ followed by an advanced workshop held in July 2023, focused on utilizing pyCHARMM for high-accuracy, high-throughput free energy calculations.⁵⁸ The Jupyter Notebooks and scripts associated with both workshops are also available through the GitHub page of Charles Brooks' lab: https://github.com/BrooksResearchGroup-UM. The release and ongoing development of pyCHARMM represent an important milestone for integrating biomolecular modeling, FFs, advanced simulation, sampling, and docking protocols, into the widely used Python programming language.

### 3.2. crimm.
Despite the best efforts from the developers of major biomolecular modeling and simulation softwares, there exists a substantial barrier when researchers first start to learn these tools. Frustrations often arise from the unfamiliarity of command scripts and in system preparation protocols that involve multiple steps to process macromolecules, small molecule ligands, water, and ions. A number of structural preparation tools have been developed for this purpose. For example, CHARMM-GUI (Section 3.3), originally designed to prepare simulation systems for CHARMM, provides a web interface to assist users with building a system for simulation.⁴⁹,⁵⁹⁻⁶² Despite convenience, it lacks scriptability and integratability. Other tools such as PDBfixer from OpenMM provide Python APIs for scripting and possible integration with other tools. However, they rely on structure files such as PDB as an intermediary to pass structural information to other software packages. The limitation of the PDB file format has rendered it insufficient in keeping complete information of a macromolecular system. Python packages such as Biotite⁶³ and BioPython⁶⁴ offer adequate APIs for structure manipulation and protocols for integration with other computational tools, but they are limited in utility

for structural preparation for MD simulations. Functions such as building missing loop regions, building missing atoms, adding hydrogens, solvation, etc., are currently absent. To address issues of scriptability and integratability in structure preparation, crimm (Chemistry with ReInvented Macro-molecular Mechanics)⁶⁵ is being designed with the following software principles and aims:

1. Accurate, consistent, and complete structure information and annotations for biomolecules maintained throughout the structure preparation pipeline.

2. Intuitive object design that organizes structural entities (e.g., model, chain, residue, atom) for retrieving information and manipulating structures, thus providing greater flexibility for programming.

3. Abstraction on routines (e.g., protonation, solvation, loop building) to create high-level APIs provided in Python for a convenient and intuitive scripting on structure preparation.

4. Clear protocols and reference implementations to Adaptors (interfaces to convert between Python classes in memory) to pass structural information to other software library or platforms, where accurate and efficient transfer of data can be guaranteed.

5. Visualizations in an interactive programming environment, i.e., Jupyter Notebook, to aid examination of structures.

6. Ease of installation, free and open source, and support for all major hardware platforms to encourage adaptation.

Crimm is directly built on the BioPython library and adopts the SMCRA model (Structure, Model, Chain, Residue, and Atom) for representing structures.⁶⁴ A BioPython-based object class provides optimal classification of macromolecular chain entities (protein, RNA, DNA, oligosaccharide, etc.). Importantly, functions of BioPython can be directly called with crimm structural object as an argument. All structural objects can be directly visualized using NGLView⁶⁶ in a Jupyter Notebook/Lab.

Structure preparation in crimm begins by fetching structures from the RCSB⁶⁷ or AlphaFold DataBase⁶⁸ as mmCIF format files for the complete and consistent organization of information.⁶⁹ In the context of CHARMM, the topology generation functions with CHARMM naming conventions and the CHARMM C36 FF is used. Currently available routines to process initial structures from the RCSB include automated missing loop/residue and disulfide bond assignment based on the data in the mmCIF file, patching of titratable residues with protonation state assigned by using interface to PropKa,⁷⁰ and topology generation. A solvation module is under development. Adapters to pyCHARMM⁹ have been implemented and crimm structures can be operated on or one can run simulations with CHARMM functions via pyCHARMM. An Adapter to RDKit⁷¹ has also been implemented for small molecule ligands integral to PDB entries. These are created as mol objects in RDKit to guarantee they maintain the correct bond orders. Other features are currently being developed to aid structure preparation and will address interfaces and adaptors to packages such as OpenMM,⁸ OpenFF,⁷² and Autodock Vina.⁷³

### 3.3. CHARMM-GUI.
Since its original development in 2006,⁴⁹,⁵⁹⁻⁶² CHARMM-GUI has proven to be an ideal web-based platform (https://www.charmm-gui.org) to interactively
<figure><img src="image_10.png" /><figcaption>**Figure 1.** Main concepts of EnzyDock applied to the mechanism in the Diels–Alderase enzyme, LepI.<sup>133,134</sup> (A) Similar (mapped) atoms are marked in green. (B) EnzyDock docking with the transition state as a template ("seed") for docking the remaining states.</figcaption></figure>

build complex systems and prepare their inputs with well-established and reproducible simulation protocols for widely used simulation packages such as CHARMM, AMBER,<sup>74</sup> Desmond,<sup>75</sup> GENESIS,<sup>76</sup> GROMACS,<sup>77</sup> LAMMPS,<sup>78</sup> NAMD,<sup>79</sup> OpenMM,<sup>8</sup> and Tinker.<sup>80</sup> CHARMM-GUI has been widely adopted for various purposes and it now contains more than 20 modules designed to set up a broad range of molecular simulation systems.<sup>81-83</sup> CHARMM-GUI also provides educational resources including online lecture materials, an online user forum, and workshops. Its archives support scientific reproducibility by providing the lipid conformation library<sup>59,61,62</sup> used in membrane generation, prebuilt COVID-19 systems,<sup>81,84,85</sup> prebuilt membrane complexes,<sup>82,83</sup> and a searchable CHARMM small molecule library (CSML). Many original modules were developed as an in-house effort, but close collaborations with the developers of CHARMM and other simulation packages have been established for adding newer modules.<sup>86-88</sup>

The philosophy in CHARMM-GUI development is less about providing the nuts and bolts of molecular modeling, but instead it focuses on helping users to achieve a task, including building membrane systems,<sup>61,62,89-94</sup> modifying and solvating proteins,<sup>95,96</sup> characterizing protein–ligand interactions,<sup>97-104</sup> or modeling complex carbohydrates<sup>105-107</sup> via a streamlined interface.<sup>108-116</sup> This makes CHARMM-GUI broadly accessible to users with little experience in modeling tools while remaining useful to experts, especially for batch generation of systems.

CHARMM-GUI development is not only guided by requests from general users and experts, but also in response to an emerging need for a unified platform to prepare and execute various advanced simulation approaches developed in diverse simulation communities and packages.<sup>113,117-120</sup> In addition to building complex molecular systems, CHARMM-GUI also assists with preparing input files for both general and advanced modeling and simulation tasks.

## 4. DOCKING METHODS

**4.1. CDOCKER.** First introduced in 2003,<sup>55</sup> CDOCKER provides an integrated CHARMM-based scripting framework for small molecule-receptor docking studies. It employs a numerical grid-based representation for the van der Waals and electrostatic interactions utilizing a fully molecular mechanics (MM)-based FF representation of the interactions.<sup>55</sup> CDOCKER utilizes conformational search based on simulated annealing, and it is also compatible with enhanced sampling

search approaches such as self-guided Langevin dynamics.<sup>58</sup> It has been used in a broad range of applications, including early efforts in community-based docking.<sup>121-125</sup> In this capacity it served as a platform to explore a range of docking and scoring approaches, including some of the early flexible receptor<sup>126</sup> and covalent docking methods.<sup>56</sup>

In the past few years, CDOCKER has been significantly updated to utilize accelerated platforms such as GPUs.<sup>13-15</sup> While the basic philosophy has remained centered on sampling via (accelerated) simulated annealing and structured around the state-of-the-art small molecule and biomacromolecular FFs, new fast Fourier transform (FFT)-based approaches have been introduced for binding pocket and ligandable-site discovery via functional probe docking,<sup>127</sup> representing important hydrogen bonding by use of hydrogen-bond-specific donor–acceptor grids,<sup>14</sup> and hybrid sampling methods that combine simulated annealing with genetic algorithm moves.<sup>14</sup>

CDOCKER has also been implemented as a package within pyCHARMM (Section 3.1). In addition to providing full access to the methods available within CDOCKER, pyCHARMM greatly simplifies the workflow through use of ‘best practice’ parameter choices and a single callable pyCHARMM command. This enables large-scale virtual screening via a single script that integrates ligand building via RDKit and SMILES strings, parametrization of ligands with small molecule FF parameter estimators such as CGenFF,<sup>17,18</sup> GAFF,<sup>20</sup> LigParGen,<sup>22</sup> and OpenFF,<sup>128</sup> automated protein grid generation, parallel docking, clustering of results and ranking, including reranking with implicit solvent models such as GBSA/GBSW/GBMV, and FACTS.<sup>129-132</sup> In summary, CDOCKER is a fast, flexible and accurate GPU-accelerated molecular docking engine that can handle cases from high-throughput small probe docking to flexible receptor–ligand docking.

**4.2. EnzyDock.** Modeling enzyme reactions requires a carefully designed computational protocol that relies on well-established theoretical foundation. The starting point is reliable 3-dimensional (3D) structures of the substrate, product, intermediates, or transition states bound to the enzyme. EnzyDock<sup>135</sup> is a CHARMM-based docking program like the well-known CDOCKER<sup>56,136,137</sup> (Section 4.1), with emphasis on enzymes. Its main feature is mechanism-based multistate consensus docking that allows the docking of reaction substrate, intermediates, transition states, and products in a mechanistically consistent and induced-fit manner (Figure 1). EnzyDock is written as a series of CHARMM scripts (>10,000
lines of script code), Python codes (~3,000 lines), and shell scripts. EnzyDock is a docking-tool and it does not compute free energy profiles that can be obtained using other methods in CHARMM such as umbrella sampling (US),¹³⁸ string-based methods¹³⁹ (Section 7.3), or metadynamics.¹⁴⁰,¹⁴¹ Consensus docking in EnzyDock is achieved by applying geometric restraints implemented via NOE restraints on reaction states relative to a predetermined “seed” state, such that all states are docked with similar poses under a given user-defined threshold (Figure 1). For instance, the seed state could be a tightly bound transition state or a known inhibitor-bound state. Conversely, if unrestrained multistate docking is performed, a reaction pathfinder module identifies all matching poses along a reaction path.¹⁴² Additional restraints such as on dihedral angles can enforce specific stereo- and regio-chemistry during docking, while positional harmonic and NOE restraints can be employed to include chemical information such as the initial cleavage site, nucleophilic attack, or ligand positions relative to key active-site residues or cofactors. Different protonation states of enzyme and cofactors during docking of different states is facilitated via CHARMM patching. Sampling of configurational space is performed using MD³ or Monte Carlo (MC)¹⁴³ simulated annealing on a grid representing the enzyme,⁵⁵ and poses are scored using the C36¹⁴⁴ and CGenFF¹⁴⁵ FFs. Flexible residues, cofactors, and waters are treated as explicit atoms on the grid. Following ligand pose clustering, final energy minimization and scoring is performed using all-atom description of the entire system and optional refinement using a QM/MM approach¹³⁵,¹⁴⁶ with a range of QM methods, e.g., semiempirical (SE)¹⁴⁷ or density functional theory (DFT)¹⁴⁸ (Section 11). Bulk solvation is modeled using an implicit solvation model (e.g., GB).³

EnzyDock has been applied to diverse systems such as terpene synthases, racemases, Diels-Alderases, phosphotriesterase,¹⁴⁹ and covalently bound ligands.¹³⁵,¹⁵⁰ It has also been used along with other docking programs in a benchmark study on ligand binding in the main protease in SARS-CoV-2.¹⁵¹

From a user perspective, the enzyme must be provided as a PDB file or CHARMM PSF and CRD (coordinate) files, ligand states as PDB files or SMILES strings, and atom mapping between similar states along a reaction path must be provided by the user. Additional restraints can be provided by the user. EnzyDock is available via GitHub and has recently been implemented in CHARMM-GUI⁴⁹ (Section 3.3).

**4.3. CIFDock.** Accurately modeling protein and ligand flexibility is vital when using molecular docking to elucidate binding modes and predict binding affinity. Binding events may rely on “induced fit” where the ligand induces conformational changes in the protein binding site.¹⁵² Accounting for induced fit has been shown to be critical for accurate modeling of the complex.¹⁵³,¹⁵⁴ To this end, we developed a novel CHARMM-based induced fit docking protocol (CIFDock)¹⁵⁵ that employs all-atom FFs and enhanced sampling MD.

The CIFDock protocol begins with processing the protein structure through CHARMM-GUI⁹⁶ to fix bond orders, add hydrogens, and correct protonation states of the protein residues. The resulting PDB file is then fed into a series of CHARMM scripts which separate the protein, ligand, ion, and water molecules into CHARMM-compatible structure files, and they will be combined during subsequent steps. Key to the CIFDock protocol is the definition of active site residues that are mutated to alanine, to allow for a more “open” active site that can accommodate larger ligands and facilitate greater

ligand conformational searching. In the final preparation step, the Confab module of OpenBabel¹⁵⁶ is used to generate a ligand conformational ensemble to seed initial binding pose searching.

The main docking procedure begins with the initial placement of the ligand in the active site of the protein in a random orientation. The ligands are then sampled using a 20-ps self-guided Langevin dynamics (SGLD)¹⁵⁷ simulation. Following this step, pairwise root-mean-square deviation (RMSD) clustering of ligand conformations is performed using the CORREL module to avoid further sampling of overlapping conformations. Each cluster is saved as a trajectory file which consists of conformations that were within a predefined cutoff radius of the cluster center. Each of the resulting protein-ligand complexes is then “backmutated” (i.e., the residues mutated to Ala are mutated back to their original residues) and a random dihedral angle-based rotamer library¹⁵⁸ is generated, and side chains are relaxed by a short energy minimization and SGLD simulation. Explicit water molecules and ions saved in the preparation stage are added back, and a second SGLD simulation is conducted on the active site complex.

The resulting “docked” poses are scored and ranked using a set of custom scoring functions that are based on the well-validated SWISSDOCK¹⁵⁹ scoring function. They are linear combinations of energy terms calculated by CHARMM, which include FF-based energies and the GMBV II implicit solvent model for solvation energy.¹⁶⁰ CIFDock was validated by cross-docking studies on a set of 21 pharmaceutically relevant proteins. Results obtained were comparable to, or in some cases improved upon, commercial docking programs. This can be attributed to the treatment of the ligand, active site, and explicit waters as fully flexible components during the docking procedure. Additionally, because CIFDock is based on short classical MD simulations, its computational cost is minimal.

To handle the formation of covalent bonds and allow covalent inhibitors to be screened, we integrated both MNDO and SCC-DFTB¹⁶¹⁻¹⁶⁵ minimizations (Section 11) into the CIFDock workflow. These minimizations together with additional dynamics simulations using positional restraints ensure adequate protein-ligand complex sampling pre- and postreaction. The covalent-based CIFDock (CovCIFDock) has been validated on a cross-dock and self-dock test set,¹⁶⁶,¹⁶⁷ with an average RMSD of 1.91 and 1.89 Å, respectively, and a 76% success rate. This compares favorably with commercial covalent docking programs such as Schrödinger's CovDock-LO (Lead Optimization) that has a 74% success rate on the same test sets. The hybrid QM/MM minimizations also add little computational overhead to the docking procedure.

## 5. FREE ENERGY METHODS

**5.1. λ-Dynamics, Multisite λ-Dynamics, and Constant pH MD.** Alchemical free energy simulations are an important class of statistical mechanical methods used in computing free energy values and differences in small molecule design and refinement,²⁸,¹⁶⁸,¹⁶⁹ as well as protein design¹⁷⁰,¹⁷¹ and CpHMD simulations.⁵⁰,¹⁷² Alchemical methods determine free energy differences by simulating chemical transformations along a non-physical pathway, often using a chemical progress variable λ. λ-dynamics is a particularly efficient and scalable alchemical method that takes advantage of natural fluctuations in the systems being studied to “drive” the chemical coordinate between the desired end points, and is generalizable to
multidimensional chemical spaces, allowing exploration of many substituents at a site or even at multiple sites (MSλD) in a single simulation.²⁶,¹⁷³

For two states A and B of a molecular species (e.g., a protein or a side chain) and the environment (e.g., solvent and/or the receptor pocket), the alchemical hybrid Hamiltonian (or Lagrangian) for the λ-dynamics is

$$ 
\begin{aligned}
\mathcal{H}(\mathbf{r}_A, \mathbf{r}_B, \mathbf{R}_e, \lambda) &= (1 - \lambda)U_A(\mathbf{r}_A, \mathbf{R}_e) + \lambda U_B(\mathbf{r}_B, \mathbf{R}_e) \\
&\quad + U_e(\mathbf{R}_e) + U_{\text{bias}}(\lambda) + K_A + K_B + K_e \\
&\quad + K_{\lambda} \\
&= \mathcal{U}(\mathbf{r}_A, \mathbf{r}_B, \mathbf{R}_e, \lambda) + U_{\text{bias}}(\lambda) + K_A + K_B \\
&\quad + K_e + K_{\lambda}
\end{aligned}
\quad (1) $$

where $U_{A/B}$ represent potential energy of A/B interacting with themselves and the environment and $U_e$ is the potential energy of the environment itself. Terms involving these three potential energies are denoted together as $\mathcal{U}(\mathbf{r}_A, \mathbf{r}_B, \mathbf{R}_e, \lambda)$. $U_{\text{bias}}$ is a biasing (umbrella) potential to facilitate sampling in the chemical coordinate $\lambda$. $K_p$ ($p \in \{A, B, e, \lambda\}$) is the corresponding kinetic energy term for the conformational or chemical variable. From Eq. 1, one can derive coupled equations of motion for the atomic coordinates $\mathbf{r}_{A/B}$ and $\mathbf{R}_e$ and the chemical coordinate $\lambda$ with a suitably assigned mass. Integrating the equations of motion subject to a holonomic constraint on $\lambda \in [0, 1]$ allows sampling of the “extended system” in a statistical ensemble of choice.¹⁷³ In the canonical ensemble, the partition function is

$$ Z(\lambda, T) = \frac{\int \delta(\lambda - \lambda') e^{-\beta[\mathcal{U}(\mathbf{r}_A, \mathbf{r}_B, \mathbf{R}_e, \lambda') + U_{\text{bias}}(\lambda')] \mathbf{d}\mathbf{r}_A \mathbf{d}\mathbf{r}_B \mathbf{d}\mathbf{R}_e \mathbf{d}\lambda'}{\int e^{-\beta[\mathcal{U}(\mathbf{r}_A, \mathbf{r}_B, \mathbf{R}_e, \lambda') + U_{\text{bias}}(\lambda')] \mathbf{d}\mathbf{r}_A \mathbf{d}\mathbf{r}_B \mathbf{d}\mathbf{R}_e \mathbf{d}\lambda'} \quad (2) $$

where $\delta(\lambda - \lambda')$ is the Dirac-δ function and $\beta = 1/k_B T$ is the inverse temperature ($k_B$: Boltzmann constant, $T$: temperature). It follows that $\Delta G_{AB}$ is given by

$$ e^{-\beta \Delta G_{AB}(T)} = \frac{Z(\lambda = 1, T)}{Z(\lambda = 0, T)} \quad (3) $$

Extension to multiple sites at one or multiple positions of a scaffold is generalized from the terms in $\mathcal{U}(\mathbf{r}_A, \mathbf{r}_B, \mathbf{R}_e, \lambda)$ above to

$$ 
\begin{aligned}
\mathcal{U}(\mathbf{r}_{\{s\}_{\{i\}}}, \mathbf{R}_e, \lambda) = & U_e(\mathbf{R}_e) + \sum_{s}^{M} \sum_{i}^{N_s} \lambda_{si} U(\mathbf{R}_e, \mathbf{r}_{si}) \\
& + \sum_{s}^{M} \sum_{t>s} \sum_{i}^{N_t} \sum_{j}^{N_t} U(\mathbf{r}_{si}, \mathbf{r}_{tj}) + U_{\text{bias}}(\{\lambda\})
\end{aligned}
\quad (4) $$

Thus, each substituent i of the $N_s$ substituents at each site s of the M total sites gets its own $\lambda_{si}$. Interactions of a substituent with itself and the environment $U(\mathbf{R}_e, \mathbf{r}_{si})$ are scaled by $\lambda_{si}$ while interactions between sites $U(\mathbf{r}_{si}, \mathbf{r}_{tj})$ are scaled by the product of $\lambda_s \lambda_{tj}$, and all remaining interactions, $U_e(\mathbf{R}_e)$, are unscaled. Although λ-dynamics has been primarily implemented in CHARMM, it can be implemented in OpenMM using custom non-bonded forces.⁸ However, for large chemical spaces the

computational efficiency is poor. CpHMD methods based on λ-dynamics have also been implemented in GROMACS,⁶ Amber,⁵ and AMOEBA.¹⁷⁴

In CHARMM, λ-dynamics is implemented through the BLOCK module, where many new features have been introduced to improve the accuracy, robustness, scope, and sampling. It is computationally expedient to ensure that λ remains between 0 and 1 (boundary constraint), and all λ values at a particular site add up to 1 (normalization). While these criteria can be maintained approximately or exactly with restraints or constraints, respectively, as was done in the earliest implementations,¹⁷³ it is more convenient to maintain them with implicit constraints through change of variables. This provides an alternative set of alchemical variables θ that map back to λ such that the constraints and normalization are satisfied by construction.²⁶

Soft-core interactions that remove non-bonded singularities near the alchemical end points of 0 or 1 are important for convergence. They are especially critical for the accuracy and reproducibility of λ-dynamics because free energy is estimated by binning together states near alchemical end points, where hard cores can lead to very sharp changes in the free energy. The BLOCK module contains a special set of soft core functions for λ-dynamics that enables van der Waals and electrostatic interactions to be turned off concurrently.¹⁷⁵ For λ-dynamics, the PME electrostatics³⁵,¹⁷⁶ gives better results than force switching electrostatics,³⁶ especially for longer simulations.¹⁷⁷,¹⁷⁸ The BLOCK module includes commands to enable a generalization of PME for λ dynamics.¹⁷⁹ It is also worth noting that the MSLD command ('L' for λ) in the BLOCK module accepts an FFIX option that will run otherwise identical simulations, but with fixed values of λ for FEP validation or discrete λ sampling.¹⁰,¹⁸⁰⁻¹⁸² Several additions to BLOCK allow broader applicability of λ-dynamics to more unusual perturbations. Protein mutations to proline, and ligand calculations involving ring changes, core hopping, or macrocyclization require special considerations to ensure that when a substituent is non-interacting at λ = 0, the dummy atoms in the substituent are only bonded to one environment atom so they do not exert a net force on the rest of the system.¹⁸³ To satisfy these considerations, the BLOCK RMLA command allows removing λ scaling for classes of interactions, and it is recommended to only leave bond and angle interactions unscaled and to scale dihedrals. For finer granularity, soft bonds are implemented in BLOCK to, for example, break the proline ring at λ = 0 and allow free rotation of other amino acids at the same site around their λ backbone angle.¹⁸⁴,¹⁸⁵ If significant portions of a molecule are similar but cannot be incorporated into the common core of a hybrid topology model due to differing charge or atom types, they may be harmonically restrained together with their bonded interactions scaled with the CATS command in BLOCK,¹⁸⁵ analogous to a similar process in NAMD.⁷

Another set of features crucial for sampling of the chemical space is adaptive landscape flattening (ALF) where a biasing potential in the λ space is iteratively developed to flatten the chemical landscape for enhanced sampling.¹⁷⁵,¹⁷⁷,¹⁷⁸ These biases are implemented by the LDIN and LDBV commands and are typically tuned by an external ALF python package.¹⁷⁵ Sampling can also be improved with Hamiltonian replica-exchange MD (REMD) through the REPD module.²⁷,³¹ More rapid sampling can be achieved with the BLADE module²⁴ (Section 2.3).
The above developments enabled sampling of massive chemical spaces spanning 512 HIV reverse transcriptase inhibitors,²⁷ 240 T4 lysozyme mutants,¹⁷⁷ and 32768 ribonuclease H variants,¹⁷¹ as well as challenging perturbations of both ligands¹⁸⁶ and proteins.¹⁸⁵

## 5.2. Hybrid Sampling and Free Energy Algorithms.

The calculation of solvation free energy and binding affinity of small molecules to macromolecules are among the most important practical applications of MD simulations, especially with the potential impact on drug discovery efforts. A wide range of methodological advances were implemented in CHARMM to improve the statistical convergence and physical accuracy of free energy calculations. Conceptual advances in free energy methodologies implemented in CHARMM were reviewed in ref 187. For example, a version of λ-dynamics was introduced via a MC multicanonical REMD (FEP/REMD).¹⁸⁸,¹⁸⁹ Specifically, the FEP/REMD helps resolve the poor convergence of the free energy estimates as a function of λ near the end points (λ = 0 and 1), which is often reflected as hysteresis between the forward (0 → 1) and backward (1 → 0) calculations from traditional FEP calculations based on single trajectories.

Applications to the calculation of the binding free energy of different kinase inhibitors demonstrated that the FEP/REMD algorithm was critical for tackling complex ligands accurately.¹⁹⁰⁻¹⁹² A similar general strategy improved the convergence of multidimensional US calculations by swapping configurations from different windows via Hamiltonian REMD (US/H-REMD)¹⁹³ (Section 7.1). Another issue concerns the sampling of solvent configurations. The binding of a ligand to a receptor frequently involves the displacement of a certain number of bound water molecules. This is not an issue if the binding site is in direct contact with the bulk solution. However, the convergence and accuracy in FEP/MD calculations can be severely compromised when a binding site is deeply buried and is inaccessible to bulk water. In this case, simple MD does not guarantee a complete sampling of the solvent during the FEP calculation. As an illustration, the binding of camphor to a deeply buried pocket in cytochrome P450cam causes about 7 water molecules to be expelled.¹⁹⁴ To address this, standard MD was coupled with the grand canonical MC (GCMC) algorithm to allow the number of water to fluctuate in any chosen region during an alchemical FEP calculation.¹⁹⁴ GCMC helps better sample the solvent configurations in the binding pocket that are poorly accessible to bulk solvent. It is also powerful by introducing fluctuations in the number of solvent molecules in FEP calculations carried out with a reduced model where only the region surrounding the binding site is explicitly considered while the effect of the surrounding solvent and protein is mimicked implicitly with the generalized solvent boundary potential (GSBP).¹⁹⁵ Such a strategy made it possible to calculate the standard binding free energy of antibiotics to the peptidyl-transferase P-site of the bacterial ribosome.¹⁹⁶,¹⁹⁷

Over the years, increasing efforts were made to streamline free energy calculations, enabling automated calculation of the absolute solvation free energy of a large number of small drug-like molecules using explicit solvent.¹⁹⁸ Moreover, collaborative efforts were made to test the accuracy and reproducibility of free energy calculations across different software packages.¹⁹⁹ One of the principal advantage of CHARMM is that different methodologies can be naturally integrated within a single job. For example, a US formulation of equilibrium

binding²⁰⁰ was used to characterize the binding specificity of a large number of SH2 domains²⁰¹ with the generalized Born with a simple switching (GBSW) implicit solvent model.¹³¹ As another example, the PBEQ continuum electrostatics module of CHARMM¹⁰⁸,²⁰² conveniently allows one to directly access and read MD trajectory snapshots, and then combine its MM potential energy together with the solvation contribution based on the Poisson-Boltzmann and surface area approximation (PBSA). This MM/PBSA strategy, seamlessly integrated within CHARMM, has been used, for example, to process a large number of protein complexes to assess the binding specificity within a family of synaptic surface receptors.²⁰³

A growing family of hybrid sampling methods combining the strength of MD and Metropolis MC were tested and implemented, benefiting from the flexibility of the control flow from the native CHARMM scripting command language at the level of the input file.²⁰⁴,²⁰⁵ These algorithms typically consider new configurations generated by driving the system via a non-equilibrium MD (NEMD) trajectory that are subsequently treated as putative candidates for MC acceptance or rejection.²⁰⁴,²⁰⁵ The hybrid NEMD/MC algorithms can be exploited in a variety of context and offer a promising avenue to sample the configurations of complex systems. For example, the discrete ionization state of titratable residues can be sampled, effectively as a constant-pH simulation.²⁰⁶ Another example is to consider new configurations of an all-atom system generated by driving it via NEMD toward a configuration that originated from a CG simulation. It was shown that the CG-guided hybrid NEMD/MC algorithm can enhance the sampling of solvated peptides even with fairly rudimentary CG models as a guide.²⁰⁷

## 5.3. Optimal Variance Alchemical Path for Free Energy Calculation.
Despite continuous development of free energy calculation methods,²⁸,⁹⁹,²⁰⁸ practical challenges impede their precision and possibly reliability.⁴⁵,¹⁹⁹ Options for improvement include enhanced sampling,⁹⁹,¹⁸⁹,²⁰⁹⁻²¹¹ careful design of alchemical cycles,¹⁹⁹,²¹¹,²¹² variational and integration approaches,²¹¹,²¹³⁻²¹⁶ and the design of the alchemical path itself,⁹⁹,²¹⁵,²¹⁷⁻²²⁰ the latter being the focus of this section.

The hybrid Hamiltonian method relies on the ‘optimal alchemical path’ theory.²¹⁷ To overcome barriers between reactant and product phase spaces, it is implemented at the interaction pair level, treating each pair separately though in parallel. Denoting abolished (A) interacting pairs as $p_i \in P_A$ and created (B) ones as $p_i \in P_B$, the corresponding Hamiltonian contributions are

$$ 
\begin{aligned}
H_{Oab}(p_i, \lambda) &= -\frac{2}{\beta} \ln([1 - \lambda] + \lambda e^{-\beta(H_A(p_i) - G_{A,0,i})/2}) \\
H_{Ocr}(p_i, \lambda) &= -\frac{2}{\beta} \ln([1 - \lambda]e^{-\beta(H_B(p_i) - G_{B,0,i})/2} + \lambda)
\end{aligned}
\qquad (5) $$

where $H_A$, $G_{A,0,i}$, $H_B$, and $G_{B,0,i}$ are the energy functions and estimates for the free energy of abolishment or creation of each of the pairs, $p_i$. Approximations for Eq. 5, denoted by $H_{cr}$ and $H_{ab}$, follow from Eq. 25 of ref 217, and the hybrid Hamiltonian is given by

$$ H_{HH}(\lambda) = \sum_{p_i \in P_A} H_{ab}(p_i; \lambda) + \sum_{p_i \in P_B} H_{cr}(p_i; \lambda) + H_C(\text{terms } \notin P_A \cup P_B) \qquad (6) $$
where $H_C$ is for all other terms unaffected by the transformation. The improper and proper dihedral angle fluctuations being modest, simple multiplication factors are used for created and abolished terms, respectively (isomorphous to Eq. 25 in ref 217). For Ewald sum, a linear scheme for charge, $q = q_C + (1 - \lambda)q_A + \lambda q_B$, is used.

As a result, the derivative with respect to $\lambda$ can be intertwined as an additional dimension to that of the system spatial coordinates, $\mathbf{r}$, extending Eqs. 4.6 and 4.9 of ref 176 as

$$ \frac{\partial E_{\text{Wt}}}{\partial \{\mathbf{r}, \lambda\}} = \left\{ \frac{\partial Q(q)}{\partial \{\mathbf{r}\}}, Q\left(\frac{dq}{d\lambda}\right) \right\} \cdot (\theta_{\text{rec}} * Q(q)) \quad (7) $$

where * indicates convolution, $Q$ the charge mesh, and $\theta_{\text{rec}}$ the reciprocal factor mesh defined in ref 176. An application of the method was on the R67 DHFR system that is a pseudohomotetramer, a dimer of dimers.221–223 To simulate the mutation process, the two subunits of one dimer had an A-hybrid residue at position 59 and the two subunits of the other dimer had a B-hybrid residue at position 62 (Figure 2).

<figure><img src="image_11.png" /><figcaption>**Figure 2.** Thermodynamic cycle. Inset: general cycle design; horizontal arrows: measured affinities in kcal/mol;222,223 vertical arrows: computations for tetramer (left) and the two types of dimers (right). Graphical panels: local molecular surface at the interfaces with the mutated residues displayed as spheres for WT (S59: red, H62: blue), S59A/H62L (A59: red, L62: yellow), and S59A/H62F (A59: red F62: green). For WT, a 0.41-kcal/mol entropic term is added to account for higher symmetry.222 Computed (C:), measured (M:) differences, and discrepancies (D:) are given. Global discrepancy (ΔΔΔG) for the 3 cycles provides a self-consistency check. Average standard deviation (StdDev$_{calc}$) and error (⟨error⟩$_{calc}$) were computed using autocorrelation functions224 considering $\lambda$ windows as independent. The average error ⟨error⟩$_{obs}$ and the maximum observed error $Max_{error-obs}$ that compare experimental results with calculations are also reported.</figcaption></figure>

Simulations were run sequentially for 10 discretized values of $\lambda$ from 0 to 1. Hybrid residues were also subjected to TI in their isolated acetylated and aminated form as a control.

Branches of the various thermodynamic cycles in Figure 2 are further analyzed in Figure 3. Individual curves are bell-shaped, mirroring the quadratic form of the partition function of the optimal path as function of $\lambda$ (Eq. 11 of ref 217) and yielding a linear integrand. Due to differences in the position of the maxima for the different branches of a same cycle, the

<figure><img src="image_12.png" /><figcaption>**Figure 3.** Integration along the thermodynamic cycles in Figure 2. Cumulative error estimates24 are also shown for the tetramer, each of the dimers, the sum of the dimers, and the global cycle for (A) WT to S59A/H62L, (B) WT to S59A/H62F, and (C) S59A/H62L to S59A/H62F for which only one dimer is involved since residue 59 remains as Ala. (D) Integrand for the transformation of the tetramer from WT to S59A/H62L shown in black as an example. The average for each $\lambda$ window is marked by a stepwise white line. Linear regression along the whole trace is shown as a light gray line to appraise the linearity of the integrand with respect to $\lambda$. Dashed lines mark ±100 kcal/mol. (E) Integrand for isolated hybrid residues (acetylated and aminated) in a vacuum to evaluate the intrinsic energy contributions due to the FF energy difference of the original residues.</figcaption></figure>

global cycle profiles are sinusoidal rather than quadratic, nonetheless very tempered.

As previously reported, S59A/H62L is favorable despite loss of a hydrogen bond and the formation of a small hydrophobic cavity.222 In comparison, S59A/H62F is less favorable despite good shape complementarity and creation of new hydrophobic contacts.223 Despite modest calculation effort only intended to illustrate the optimal alchemical path integrand properties, simulations reproduce those unexpected results. Interestingly, TI on isolated hybrid residues revealed the predominance of the amino acids intrinsic FF potential differences on the integrand, suggesting that reducing those differences could reduce the difficulty to reach accurate results. The linearity of the integrand with respect to $\lambda$ for the method presented here facilitates integration, hence it is a desired property. It avoids the need for evolved integration schemes that can amplify errors, but are required to treat irregularity or singularity found, for example in conventional van der Waals creation.214

# 6. CONSTANT PH METHODS

## 6.1. Hybrid-Solvent and All-Atom Continuous Constant pH Methods.
Describing protonation state changes due to a change in solution pH or conformational environment was first enabled in CHARMM through the GB CpHMD methods.225,226 In these methods, an auxiliary set of ($\lambda$) coordinates representing the evolution of protonation states are propagated based on the idea of $\lambda$-dynamics173 (Section 5.1). Since 2010, the CpHMD framework was further developed to be carried out in explicit-solvent MD simulations (see, e.g., Section 5.1). An example is the hybrid-solvent CpHMD227 that samples solute conformation in explicit solvent but leverages the GBSW implicit solvent model131 for propagating protonation states. The pH REMD method
was also developed to accelerate sampling of the coupled conformation and protonation states.²²⁷ The hybrid-solvent CpHMD was later extended for transmembrane protein simulations²²⁸ by including the implicit membrane GBSW model²²⁹ with a water cylinder to account for water molecules in the pore of a channel or a transporter. To remove the dependence on the GB models which limits the accuracy, the all-atom CpHMD methods with generalized reaction field²³⁰ or PME for long-range electrostatics¹⁷⁹ have been developed. To enforce net charge neutrality in all-atom CpHMD, an approach based on cotitrating ions²³⁰ or water²³¹ has been developed. The hybrid-solvent and all-atom CpHMD have enabled not only new lines of inquiries, e.g., pH-dependent self-assembly mechanism of chitosan in which a total of 160 glucosamine units were allowed to titrate,²³² but they also provided fresh perspectives to resolve old questions where, e.g., the hybrid-solvent CpHMD simulations revealed the formation of proton-coupled hydrogen bonds as a major determinant for acid/base.²³³

## 6.2. Constant pH MD with Discrete Protonation States.
There are two main classes of constant pH simulations depending on whether the protonation states vary discretely (either deprotonated or protonated)²⁰⁷,²³⁴⁻²⁴¹ or continuously.¹⁷²,²²⁵⁻²²⁷ Two types of the former class are implemented in CHARMM. The first is based on the MD/ MC constant pH method²³⁷,²³⁸ that is available only in implicit solvent. The second type is based on the EDS method,³⁷ and is available for explicit solvent. Constant pH simulations with continuous protonation states, also in CHARMM, are available for implicit,²²⁵ combined implicit and explicit,²²⁷ as well as explicit-only solvent.¹⁷²

The MD/MC method,²³⁷ originally available for Amber, has been implemented in CHARMM and further extended to include constant pH REMD.²³⁸ During MD simulation, attempts to change the protonation state according to the Metropolis criterion are made at a user-defined interval. The deprotonated state is modeled with the proton present following the charge distribution for the deprotonated state. For replica *i*, let the positions and momenta of atoms be $q_i$ and $p_i$, respectively, $N_i^p$ be the number of titratable residues that are protonated, and $p_H$ its pH. Similarly define $q_H$, $p_H$, $N_H^p$ and $p_H^m$ for replica *j*. Denoting $X_i^l \equiv (q_i, p_i, N_i^p, p_H)$ and $X_i^m \equiv (q_i, p_i, N_i^p, p_H^m)$, the probability of exchange between replica *i* and replica *j* is

$$w(X_i^l, X_j^m \rightarrow X_i^m, X_j^l) = \begin{cases} 1 & (\Delta \le 0) \\ e^{-\Delta} & (\Delta > 0) \end{cases} \quad (8)$$

$$\Delta \equiv \ln 10 \cdot (pH_m - pH_i)(N_i^p - N_j^p)$$

In addition to the constant pH REMD that greatly improves sampling of the protonation state,²³⁸ the reservoir constant pH REMD method was developed to better sample conformational states.²⁴³ It relies on pregenerated reservoirs of conformations with fixed protonation states. The reservoirs can be generated either by long MD simulations, or with an enhanced sampling method, so that conformations with a given protonation state follow the Boltzmann distribution. Then an attempt to replace the current conformation with a random reservoir structure is made after a given number of steps. The attempt is accepted if protonation states of all ionizable residues match with those of the reservoir structure, and

rejected if not. In this way, the system can sample conformations from the Boltzmann ensemble of the reservoir.

Another method implemented in CHARMM is the EDS with Hamiltonian REMD (EDS-HREM).⁴⁰ In the EDS approach,³⁷ a hybrid Hamiltonian enveloping both states is defined such that the corresponding partition function is the sum of partition functions for individual Hamiltonians. In addition, a smoothness parameter can be introduced to facilitate conformational transitions between states with high energy barrier. In its constant pH implementation,⁴⁰ the two states are protonated and deprotonated, and a pH-dependent energy offset between the two states is introduced.

$$E_{\text{EDS}}(\mathbf{x}, s, \text{pH}) = -\frac{1}{\beta s} \ln \left( \sum_{i=1}^{N} e^{-\beta s[E_i(\mathbf{x}) - E_i^{\text{offset}}(\text{pH})]} \right) \quad (9)$$

Here, $E_i(\mathbf{x})$ is the potential energy of state *i* with coordinate $\mathbf{x}$, *s* is the smoothness parameter, and $E_i^{\text{offset}}(\text{pH})$ is the pH-dependent energy offset calculated ahead of the simulation via thermodynamic cycling. The Hamiltonian can be extended to several titrating groups, conveniently describing clusters of coupled residues. For example, it has been used to calculate the $pK_a$ values of four glutamic acid residues in the selectivity filter of a sodium channel.²⁴⁴

Different replicas have different values of *s*, which allows for replicas with low *s* (very smoothed) to cross energy barriers, while replica with *s* = 1 yields the conformational ensemble identical to the semigrand canonical ensemble at convergence. As a follow-up, a 2-dimensional (2D) replica exchange pH method was added in CHARMM, where the second dimension is pH.⁴¹

FEP methods have also been used for protein $pK_a$ calculation in both implicit and explicit solvent,²⁴⁵ as well as in QM/MM settings.²⁴⁶ In a recent study, $pK_a$ calculations from the 2D EDS-HREM in explicit solvent have been found to agree well with FEP results for a complicated system consisting of four selectivity filter glutamate residues of an ion channel with bound ions.²⁴⁷ Additional FEP simulations led to a new proposed mechanism of selectivity in this ion channel, based on the shift of the $pK_a$ value in the presence of different ions.²⁴⁷

## 6.3. Proton Hopping Simulations.
Classical biomolecular MD simulations normally do not allow changes in covalent bonding. This is an issue in systems involving proton transfer, as e.g. in proton diffusion in water where a proton breaks a bond with one water and forms a new one with a neighboring water molecule. The MOBHY (for “mobile hydrogen”) module in CHARMM allows proton mobility by interspersing discrete proton moves during a dynamics trajectory.²⁴⁸ After a given number of MD steps, an attempt is made to move a titratable proton to an eligible alternative location, i.e., a potential acceptor to which the titratable proton is hydrogen bonded. Upon the hop attempt, the molecular geometries and FF parameters of protonated and deprotonated species are changed accordingly. The missing protons are represented by dummy atoms (no charge and no interactions with surroundings). The initial protein structure is generated with all potential protons present, i.e., all specified titratable residues should be fully protonated in the PSF (whether they are truly protonated is selected by the user). Thus, no actual changes in bonding take place during a proton hopping simulation; only the atom types and charges change. The excess proton is represented as a classical hydronium ion. Acceptance of a
proton move is based on a Metropolis-like criterion that employs an empirical threshold for the energy change upon proton hopping. The threshold is chosen to reproduce the experimental proton diffusion coefficient in water. Similar empirical thresholds are used for proton hopping between water and protein side chains, while the true rates can be obtained by more elaborate methods.²⁴⁹ This method has been applied to proton conduction by gramicidin A₇,²⁴⁸ investigation of the asymmetry of proton conduction in the influenza M2 proton channel,²⁵⁰ and evaluation of models for the human voltage gated proton channel.²⁵¹

## 7. ENHANCED SAMPLING AND TRANSITION PATH METHODS

**7.1. Replica Exchange MD (REMD).** In REMD, N independent copies (or replicas) of a system are run in parallel and are periodically swapped (i.e., exchanged) to enhance the crossing of potential energy barriers.²⁵²,²⁵³ REMD is useful in systems where energy barriers lead to poor sampling and slow convergence in conventional MD, hindering accurate calculation of thermodynamic quantities.¹⁸⁸,¹⁸⁹,²⁵⁴

When the system volume does not change, the probability of observing a system in a configuration represented by coordinates X and Hamiltonian a with energy $E_a \equiv E(X)$ is

$$P(X, E_a) = \frac{e^{-E(X)\beta_a}}{Z_a} \qquad (10)$$

where $\beta_a = 1/k_B T_a$, with Boltzmann constant $k_B$ and temperature $T_a$. $Z_a$ is the partition function. For N non-interacting replicas, the probability of observing the system in a particular state is the product of the probabilities for individual replicas:

$$P = \prod_{n=1}^{N} P_n \qquad (11)$$

The enhancement of sampling in REMD comes from periodic swapping of the coordinates and velocities between two replicas. By imposing detailed balance, the ratio of the forward and backward transition rates between replicas a and b in exchanging their coordinates X and Y is given by²⁵³

$$\frac{P_{\text{forward}}}{P_{\text{backward}}} = e^{(\beta_b - \beta_a)(E(Y) - E(X))} \equiv e^{\Delta_T} \qquad (12)$$

where $e^{\Delta_T}$ is the temperature replica exchange probability. In Hamiltonian REMD, $e^{\Delta_H}$ can be similarly defined. They are used to accept or reject the exchange using the Metropolis criterion:

$$P_{\text{exchange}} = \begin{cases} 1 & (\Delta_{\{H,T\}} \ge 0) \\ e^{\Delta_{\{H,T\}}} & (\Delta_{\{H,T\}} < 0) \end{cases} \qquad (13)$$

When the system's volume changes (NPT ensemble), Eq. 10 changes to²⁵⁵

$$P(X, E_a) = \frac{e^{-\beta_a [E(X) - P_a V(X)]}}{Z_a} \qquad (14)$$

where $P_a$ is the external pressure at Hamiltonian a and $V(X)$ is the volume of the coordinates X. The exponents of $P_{\text{exchange}}$ then become

$$\begin{aligned} \Delta_{\{H,T\},P} &= \Delta_{\{H,T\}} + \Delta_{\text{correction}} \\ \Delta_{\text{correction}} &= (\beta_b P_b - \beta_a P_a)(V(Y) - V(X)) \end{aligned} \qquad (15)$$

REMD in CHARMM is handled through the REPD (REPlica Distributed) command. It requires MPI parallelism with one or more MPI processes per replica. Exchanges are attempted at a user-specified interval, typically on the order of 1 ps. The exchange direction alternates between “up” and “down” in the replica space. While it is not strictly necessary to attempt exchanges only between neighboring replicas, acceptance of an exchange between two replicas requires overlap between their potential energy distributions that is typically highest for neighboring replicas. Exchanges are accomplished by swapping coordinates and velocities between MPI processes, so that each MPI process yields a “replica” trajectory, i.e., coordinate frames corresponding to a single temperature or Hamiltonian.

CHARMM supports REMD for temperature, general Hamiltonian, self-guided Langevin,²⁵⁶ and CpHMD using either discrete²³⁷,²³⁸ or continuous¹⁷⁹,²²⁵,²²⁷ protonation states. CHARMM also supports coupling of the top and/or bottom replicas (i.e., the highest and lowest in replica space) to pregenerated structure reservoirs. Exchanges with the reservoir can further accelerate conformational sampling,²⁵⁷ and can be done assuming either Boltzmann (recommended)²⁵⁸ or non-Boltzmann²⁵⁹ weighting. In constant pH REMD, CHARMM supports exchanges with reservoirs that have fixed protonation states, where exchanges with a structure in the reservoir can only be accepted if the protonation state of all ionizable residues matches the structure to be exchanged (Section 6.2).²⁴³

CHARMM also supports multidimensional REMD,²⁶⁰,²⁶¹ with the only restriction being that a dimension aside from the general Hamiltonian may only be used once (e.g., one temperature dimension and one self-guided Langevin dimension is permitted, but not two temperature dimensions). The combination is multiplicative: For example, a setup with 4 temperatures and 2 Hamiltonians will use 8 replicas in total. To simplify scripting, CHARMM sets up several user-accessible variables, such as ?NREP and ?MYREP, which refer to the total number of replicas and the global replica index respectively. For multidimensional REMD, ?NREPD⟨X⟩ and ?MYREPD⟨X⟩ refer to the total number of replicas and replica index in dimension ⟨X⟩, respectively (for example, ?NREPD1 is the number of replicas in the first replica dimension).

REMD in CHARMM can be combined with other ensemble methods such as EDS³⁷ via the MSCALE module.²⁶² Earlier, a constant pH method in explicit solvent with discrete protonation states was developed based on a combination of EDS and a 1D REMD.⁴⁰ A more recent version features EDS with a 2D REMD (the second dimension being pH), which significantly accelerates the convergence of constant pH simulations.⁴¹

**7.2. Biasing Methods.** **7.2.1. Targeted MD (TMD).** Conformational transition pathways can be simulated with a number of TMD methods. The original implementation²⁶³ introduces a holonomic constraint that reduces the RMSD from the target coordinates with a preset value at each MD step. While this guarantees to reach the target conformation, generated pathways are generally irreversible²⁶⁴ and they can cross large free energy barriers.²⁶⁵ By using a perturbation of a
fixed magnitude that minimizes the RMSD with the target at every step, the restricted perturbation TMD (RPTMD) method<sup>265</sup> generates low free energy pathways along which potential of mean force (PMF) profiles can be readily calculated.<sup>266</sup> The RMSD can also be decreased by a restraint potential (RTMD) that can be symmetrized to yield more reversible paths.<sup>264</sup> Due to the use of global best-fit rotations, these TMD methods tend to favor large scale motion before small conformational changes,<sup>265,267</sup> which is subdued in locally restrained TMD (LRTMD) by applying a number of TMD restraints on subsets of atoms.<sup>267</sup>

### 7.2.2. Related Conformational Free Energy Sampling.
CHARMM supports a number of enhanced sampling techniques to evaluate conformational free energy differences. US<sup>138</sup> and adaptive US<sup>268,269</sup> of distances, angles, torsions, RMSD, and more complex geometrical order parameters are supported by the CONS, RXNCOR, and ADUMB modules. US is typically performed through the use of harmonic restraints that bias the system toward a desired target. CHARMM also supports best-fit positional restraints in which the reference coordinates are first rotated and translated to minimize the restraint energy. These best-fit restraints are key to the efficiency of confinement methods<sup>270-274</sup> that calculate conformational free energy differences by transforming (part of) the system to the desolvated harmonic oscillator state. The Gaussian-mixture US (GAMUS) method allows enhanced sampling of multidimensional order parameters (3-6 dimensions).<sup>275,276</sup> Like adaptive US, GAMUS uses the negative of the calculated free energy as the biasing potential, which is updated periodically while taking all sampled data into account. GAMUS constructs its biasing potential from a Gaussian-mixture model that fits the probability distribution using fully optimized Gaussian functions. By foregoing grids, GAMUS can sample higher dimensional spaces than traditional adaptive US. CHARMM also supports Tsallis-based biasing potentials<sup>277</sup> that increase sampling by reducing the force near energy barriers. In CHARMM, Tsallis-based sampling can also be coupled to replica exchange with solute tempering<sup>278,279</sup> for faster sampling.<sup>280,281</sup>

### 7.3. String Method (SM) for Conformational Transitions.
If a process of a system with positions x is described by the reaction coordinate $q(x) \in \mathbb{R}$, the free energy $\mathcal{F}$ of a conformational state $q(x) = q_0$ is

$$e^{-\beta \mathcal{F}(q_0)} = \int e^{-\beta E} \delta(q(\mathbf{x}) - q_0) \, d\mathbf{x} \qquad (16)$$

One often wishes to follow the progress of an actual chemical or physical reaction as $q_0$ is varied from the initial (reactants) to the final (products) value. Below, we focus on a set of methods in which the reaction coordinate is optimized from an initial pathway or a set of intermediate configurations.<sup>282-287</sup>

The essential idea of SM<sup>288-290</sup> is to assume that the optimized path is everywhere tangent (possibly up to a constant multiplicative tensor) to the reaction coordinate gradient without needing to specify an analytical form for it (Figure 4). Three versions of SMs implemented in CHARMM are described below.

#### 7.3.1. Zero-Temperature SM.
The zero-temperature SM (ZTSM) computes a minimal-energy path (MEP) which is a curve in the space of $N_a$ atom coordinates defined as $C = \{x(\alpha) \in \mathbb{R}^{3 \times N_a}, \alpha \in [0, 1]\}$ that, for any $\alpha \in (0, 1)$, satisfies

<figure><img src="image_13.png" /><figcaption>Figure 4. Illustration of the SM on the 2D Mueller potential. An MEP (dotted white curve) and a finite temperature string (solid white curve) connect the reactant (R; $q=0$) and product states (P; $q=1$) enclosed within red ellipses. Black contours represent isocommittor surfaces obtained from a 2nd order finite difference solution of the backward Kolmogorov equation for overdamped Langevin dynamics. White straight lines are planar approximations to the isocommittor surfaces, which also partition the configurational space into a Voronoi tessellation with nodes (red bullets). Gray dots are simulation coordinates from overdamped Langevin dynamics restrained to reaction coordinate planes and collectively define a transition tube.</figcaption></figure>

$$\mathbf{x}'(\alpha) = \frac{d\mathbf{x}}{d\alpha} \| \nabla_{\mathbf{x}} E(\mathbf{x}) \qquad (17)$$

where $\mathbf{x}(0)$ and $\mathbf{x}(1)$ correspond to the reactant and product state, respectively. The ZTSM evolves an initially assigned guess to the MEP using the steepest descent (SD) minimization while enforcing uniform parametrization by arc length, $|d\mathbf{x}/d\alpha| = \text{constant}$. In the CHARMM implementation, the continuous string is discretized into N replicas or 'images,' each assigned to a separate group of processors for parallel execution

$$\mathbf{x}_i = \mathbf{x}(\alpha_i), \quad i = 1, \dots, N \qquad (18)$$

with N typically determined by the available computing processors. The string evolves to the MEP as

$$\hat{\mathbf{x}}_i(\tau + \Delta\tau) = \mathbf{x}_i(\tau) - \Delta\tau\gamma_0^{-1}\nabla_{\mathbf{x}}E(\mathbf{x}_i) \qquad (19)$$

$$\mathbf{x}_i(\tau + \Delta\tau) = [R\hat{\mathbf{x}}(\tau + \Delta\tau)]_i \qquad (20)$$

where $\Delta\tau\gamma_0^{-1}$ controls the speed of SD evolution ($\Delta\tau$ is an artificial time step, and $\gamma_0$ is a friction constant that ensures dimensional consistency).

Eq. 19 is advanced independently for each image, and R is the reparameterization operator that corrects the provisional coordinates $\hat{\mathbf{x}}$ so that $\ln(\alpha)'1$ is constant along the string. R is common to the SMs in CHARMM (Section 7.3.4). The evolution step ($\Delta\tau\gamma_0^{-1}$) and convergence criteria can be set manually, or automatically by the SD minimizer of CHARMM. While SD is the default minimizer for ZTSM, other minimizers in CHARMM can also be used.

#### 7.3.2. Finite-Temperature SM.
The finite-temperature SM (FTSM) can be derived from the backward Kolmogorov equation (BKE)<sup>291</sup> corresponding to overdamped Langevin dynamics.<sup>292-295</sup> In FTSM, the desired reaction coordinate $q$ is assumed to be the committor function that solves the BKE,<sup>291</sup>
and the committor isosurface $q = q_0$ is approximated by a hyperplane (see Figure 4)

$$|\nabla q(\mathbf{x})|\delta(q(\mathbf{x}) - q_0) = \delta(\nu(q_0) \cdot [\mathbf{x} - \mathbf{\phi}(q_0)]) \quad (21)$$

The Jacobian $|\nabla q(\mathbf{x})|$ preserves the volume and $\nu(q_0)$ is the unit normal to the hyperplane $P_{q_0}$ that approximates the isosurface $q(x) = q_0$. $\mathbf{\phi}$ is constrained by

$$\begin{aligned} \mathbf{\phi}(q_0) &= Z(q_0)^{-1} \int \mathbf{x} e^{-\beta E} \delta(\nu(q_0) \cdot [\mathbf{x} - \mathbf{\phi}(q_0)]) \, d\mathbf{x} \\ &= \langle \mathbf{x} \rangle_{P_{q_0}} \end{aligned} \quad (22)$$

where $Z(q_0) = \int e^{-\beta E} \delta(\nu(q_0) \cdot [\mathbf{x} - \mathbf{\phi}(q_0)]) \, d\mathbf{x}$ is the partition function of the hyperplane. In analogy with an MEP, we can parametrize a continuous curve $\mathbf{\phi}(q_0(\alpha))$ having $d_{q_0}/d\alpha > 0$, and identify it with the average reaction path. Provided that the transition 'tube' (Figure 4), as measured by the variance of $\mathbf{x} - \mathbf{\phi}(q_0)$, is not too large, $\mathbf{\phi}$ also represents the dominant reaction path. From Eqs. 21 and 22, it can be shown that$^{292}$

$$\nu(q_0) \| \frac{d\mathbf{\phi}(q_0)}{d\alpha} \quad (23)$$

i.e., the reaction coordinate hyperplanes are locally perpendicular to the reaction path (string). In FTSM, Eqs. 22 and 23 are iteratively solved.$^{285,293,294}$ From an approximation to the string at iteration $n$ ($\mathbf{\phi}^n$), one obtains $\nu^n$ using Eq. 23, which permits computing $\mathbf{\phi}^{n+1}$ using Eq. 22. This is repeated until $\mathbf{\phi}^n$ does not change (up to thermal noise). The free energy can then be obtained by TI of the free energy derivatives sampled on the hyperplanes,$^{285,294}$ or by sampling a Voronoi tessellation (Figure 4).$^{290}$

The FTSM in CHARMM can optionally use Hamiltonian REMD to accelerate sampling, and an upper bound on the transition tube width can be set to limit sampling near a predefined path. In a parallel implementation,$^{296}$ FTSM starts from an initial string discretized into $N$ images $\mathbf{\phi}_i^0$, $i \in \{1, ..., N\}$, which can be obtained from, e.g., an MEP or a biased dynamics trajectory. To each image $\mathbf{\phi}_i$, one assigns a separate CPU group and a complete all-atom MD simulation system denoted by $\mathbf{x}_i$, to be used for sampling each reaction coordinate hypersurface. Each CPU group receives the neighbor images $\mathbf{\phi}_{i \pm 1}$ in addition to $\mathbf{\phi}_i$, which are required to compute $\nu(q_0)$ in Eq. 23, and samples the hyperplanes independently of the other groups.

### 7.3.3. String in Collective Variables.
There are cases when variables other than Cartesian coordinates, e.g., distances,$^{297}$ are more suitable for the reaction coordinate. Following the steps in ref 289, SM in CHARMM has been reformulated in a coarse-grained (CG) space of collective variables (CVs).$^{139}$ Assume that the reaction coordinate is determined by a set of CVs $\theta_i(\mathbf{x})$ ($j = 1, ..., K$) via some function $f$ (which does not need to be specified explicitly): $q(\mathbf{x}) = f(\theta_1(\mathbf{x}), \theta_2(\mathbf{x}), ..., \theta_K(\mathbf{x}))$. The coarse-graining leads to a $K$-dimensional free energy landscape as a function of CV coordinates denoted by $\mathbf{z}$:

$$\begin{aligned} F(\mathbf{z}) &= -\frac{1}{\beta} \ln \left( \frac{1}{Z} \int_{\mathbb{R}^{3 \times N_d}} e^{-\beta E(\mathbf{x})} \prod_{j=1}^{K} \delta(z_j - \theta_j(\mathbf{x})) \right) \\ &= -\frac{1}{\beta} \ln \langle \delta(\mathbf{z} - \mathbf{\theta}(\mathbf{x})) \rangle \end{aligned} \quad (24)$$

and a metric tensor $M(\mathbf{z})$ defined by

$$\frac{M_{jk}(\mathbf{z})}{e^{\beta F(\mathbf{z})}} = \sum_{l=1}^{N_d} \frac{1}{m_l} \langle \nabla_{\mathbf{x}_l} \theta_j(\mathbf{x}) \cdot \nabla_{\mathbf{x}_l} \theta_k(\mathbf{x}) \delta(\mathbf{z} - \mathbf{\theta}(\mathbf{x})) \rangle \quad (25)$$

where $m_l$ is the mass of atom $l$. Using $F$ and $M$, it is possible to write down Langevin equations governing the evolution of $\mathbf{z}$.289 Further, assume that the reaction proceeds via a localized reaction channel that contains a minimum free energy pathway (MFEP) on the CV landscape

$$\mathbf{z}'(\alpha) \| M \nabla_z F \quad (26)$$

with parameter $\alpha \in [0, 1]$ and $|z'(\alpha)|' = 0$ (equal arc length) in analogy with Eq. 17 for the MEP.

The SM in collective variables is an iterative algorithm for computing the MFEP using local averaging of the force $\nabla_z F$ and metric tensor $M$ obtained from restrained MD simulations.$^{139,289}$ After the string converges to the MFEP, two types of free energy profiles can be computed, $F[z(\alpha)]$ in the $K$-dimensional space of the CVs, and a 1-dimensional profile $F(\alpha)$ associated with the reaction coordinate hyperplanes on $\mathbf{x}$. An approximate calculation of $F(\alpha)$ in CHARMM is implemented using Voronoi tessellation, which also allows computation of the mean first passage time along the reaction coordinate using the Markov state model.$^{290}$

### 7.3.4. Reparameterization.
The SMs described here involve optimization of continuous curves (strings) specified by a parameter, e.g., $\{\mathbf{\phi}(\alpha), \alpha \in [0, 1]\}$. In numerical implementation, a set of discrete points along a string are used instead. To maintain uniform string resolution, parametrization by arc length is used, i.e.,

$$\alpha = \int_{0}^{\alpha} |\mathbf{\phi}'| d\alpha / \int_{0}^{1} |\mathbf{\phi}'| d\alpha \quad (27)$$

which implies that $|\mathbf{\phi}'|$ is constant along the string, or that $|\mathbf{\phi}_i - \mathbf{\phi}_{i-1}|$ is constant for all images $i > 0$. Because the string deforms as it evolves, points (images) along the curve $\mathbf{\phi}(\alpha)$ must be periodically reassigned to satisfy equidistance. This reparameterization operation (R in Eq. 20), is implemented by interpolating the string onto a refined parameter grid, i.e., $\alpha_j, j = 1, ..., N_j$ with $N_j = 5 \times N$, computing arc length on this grid normalized to the unit interval, and interpolating onto the original uniform parameter grid. Linear interpolation is the default and recommended method. Others such as B-splines and cubic splines can also be used.

### 7.3.5. String with Swarms-of-Trajectories.
Rather than refining the string in the multidimensional space of CVs by estimating the average force and metric tensor from restrained trajectories via Eqs. 24 and 25 as described above, an alternative approach considers the average dynamic drift of those variables determined on-the-fly via ensemble of short unbiased trajectories starting at different points along the string.$^{298}$ One advantage of this so-called "SM with swarms-of-trajectories" over the traditional procedure is that the computational task can be naturally distributed over many computer nodes with negligible interprocessor communication. The formal equivalence between the two approaches in the limit of very short trajectories was established,$^{299,300}$ and their respective significance has been clarified.$^{300}$

### 7.3.6. Script-Based SM Approach and Structure Building.
In a first application of the SM with swarms-of-trajectories to an all-atom solvated protein,$^{298}$ the activation pathway of Hck kinase and the inactivating DFG-flip were determined.$^{301,302}$ It
bears emphasizing that the SM could be scripted directly in the input file of CHARMM, and required no new source code. The powerful scripting facilities within CHARMM, especially the ability to modify the bonding topology of the system on the fly using the Patch Residue (PRES) facility, made it possible to generate all-atom models of the polymerized FT-30 membranes, which are widely used in reverse osmosis operations.<sup>303,304</sup>

### 7.4. Adaptively Biased Path Optimization (ABPO) for Transition Path Sampling.
Algorithms to compute the energetics and conformations associated with protein conformational transitions are most often based on path-restrained sampling using a chain-of-states defined at specified intervals along the path. The ABPO method<sup>305</sup> is an alternative approach that does not require the protein system be restrained to the path. ABPO is implemented in CHARMM through the ENSEMBLE module with options for defining CVs (also called reduced variables, RVs) and path optimization parameters. An adaptive biasing potential, V<sub>b</sub>, is utilized to enhance sampling of the path without restraining the system to specific points on the path<sup>306</sup>

$$V_b(\lambda, t) = k_B T \frac{b}{1-b} \ln[c(1-b)h(\lambda, t) + 1] \quad (28)$$

where b is the fraction of the free energy flattened by the bias, c has an inverse time unit and controls how the bias couples to the dynamics. V<sub>b</sub> adapts from the sampling histograms h(λ, t) that counts visits to the region of the path around λ over time t. The PMF is a direct result of the adaptive bias potential obtained for the optimal path.

A second distinction of ABPO compared to path-restrained methods is that evolution of the ABPO path begins by initiating multiple trajectories from an equilibrium ensemble simulated at each end state (Figure 5A). As such, the

<figure><img src="image_14.png" /><figcaption>**Figure 5.** Illustration of ABPO. (A) Energy landscape in the CV space at the initial stage of path optimization. Multiple trajectories are launched from the two end-state energy wells (blue) and sample freely along an arbitrary initial path (red line) enhanced by V<sub>b</sub> and within a tube centered on the path (white transparent rectangle) by a tube potential. (B) Trajectory visits to hyperplanes (small gray rectangles) perpendicular to the path tangent are counted. After sufficient sampling, the mean position in each hyperplane from counts over all replicate trajectories (blue X's) is determined, and the path and tube center are updated to these new values in CV space. The process is repeated to move the path incrementally (dashed arrows) until convergence to the optimal one (white curve).</figcaption></figure>

generation of unphysical structures at specified intervals along the initial chain-of-states path is avoided when starting the ABPO calculation. As trajectories move out of the end-state basin, their proximity to the path is retained with a tube potential of specified radius and centered on the path. An

advantage of free sampling within the tube is to reduce frustration in sampling a rugged free energy landscape.

Formulation of the ABPO path follows that of the finite-temperature SM<sup>307</sup> (cf., Section 7.3.2). The path is specified by CVs, the definition of which is key for the computation of the PMF. Sampling of the path is counted in terms of hits to hyperplanes orthogonal to the tangent at each path index point, and statistics over multiple trajectories in a time period are used to update the ABP. The path is evolved by computing the mean position of trajectory hits in the hyperplanes and updating the path variables to coincide with those of the mean (Figure 5B). A redistribution of the updated path index points is needed for smoothing and rescaling using a mollifier.<sup>305</sup> The optimum path is reached when the distance between the last and penultimate curves falls below a specified threshold.

The PMF A(λ, t) along the path parametrized with λ (within an additive constant) is computed from the histograms obtained from exhaustive sampling of the optimized path over time t,

$$A(\lambda, t) = -k_B T \frac{1}{1-b} \ln \left[ \frac{h(\lambda, t)}{\max[h(\lambda, t)]} \right] \quad (29)$$

As a directed approach, ABPO readily affords an atomistic description of a transition process in a reasonable simulation time depending on the choice of the selected CVs. Further, ABPO samples in a tube region surrounding the path and thereby generates a range of conformations orthogonal to the path that would not be obtained with path-restricted methods. The algorithm also provides a convenient way to assess the choice of CVs as well as the convergence of the path by following the time-course of individual CVs as a function of λ, so-called CV plots.<sup>308,309</sup> ABPO has the potential limitation of insufficient sampling in regions of high free energy, whereas path-restrained methods by nature ensure sampling all parts of the defined path.

### 7.5. Reaction Path Optimization with Holonomic Constraints.
When studying protein conformational changes, a chain of intermediate replicas of the system resolve the transition between the initial and final states. To find the most probable pathway, an objective function such as the total energy or free energy of replicas is defined and minimized.<sup>282,307</sup> Success of reaction path optimization depends on auxiliary schemes to ensure proper distribution of replicas for capturing kinetic bottlenecks.<sup>310</sup> In general, it is desirable to maintain equal distances between neighboring replicas while the distance is free to change since the actual reaction path is not known *a priori*. A folded-back path should also be avoided as replicas are placed to take forward steps in crossing kinetic barriers rather than going back and forth in a basin. In this regard, the angles between three consecutive replicas are often restrained<sup>310</sup> to prevent drastic changes in the tangent vectors along the path that are represented by the position vector differences of replicas. A key challenge of reaction path optimization is the auxiliary scheme of managing path quality interfering with the optimization of the objective function. Keeping equal distance between replicas, for example, tends to conflict with the forces along the path in energy minimization. Although the tangential component of the force can be removed,<sup>311</sup> the non-conservative projected force makes the application of fast-converging gradient-based optimization methods difficult.<sup>312</sup> The robustness and efficiency in capturing low-energy kinetic barriers are thus limited, especially with a
large number of degrees of freedom and a rugged potential energy surface (PES).

The RCONS module in CHARMM overcomes this by treating equal distance between replicas as holonomic constraints.313 Built on top of the REPLICA module, the reaction path optimization with RCONS is entirely gradient-based, readily allowing quasi-Newtonian methods and other optimization schemes assuming conservative forces. With Lagrange multipliers in constraint optimization, ad hoc numerical procedures such as rearranging atomic positions or force projections are not needed.313 Furthermore, the distance between replicas can be defined by using a non-commutative RMS best-fit procedure312 that is particularly useful for modeling transitions of macromolecules. Convergence of reaction path optimization provides a way to analyze if a sufficient number of replicas are used by testing whether the accumulated work along the optimized path agrees with the potential energy difference.313 Since the tangent vectors in this work-energy analysis are based on positional differences between replicas, the energy or free energy difference along a path can be decomposed into contributions from different atoms to deduce the kinetic bottleneck.314 It was also found that the straightness over replicas can be formulated as a kinetic energy potential and a temperature scale can be used to characterize the restraints regulating curvatures along the path.313

In principle, any potential energy function can be used to describe the energetics of replicas, and using RCONS with MSCALE provides a versatile framework for the general applications of reaction path optimization. Each replica along the path is treated as a subsystem for using a CHARMM potential energy function or in programs supported by MSCALE such as those providing a QM or QM/MM PES. For complex reactions involving conformational changes, implicit solvent model can be used to obtain an initial MEP followed by explicit-solvent MD simulations to obtain MFEP.315,316 RCONS can also be used to constrain the sampling of MD simulation over perpendicular directions to compute the PMF along a path. In this case, the chain defined by replicas is used as a 1-dimensional order parameter for the PMF calculation.315,316 Coupled with trajectory analysis, MD simulations constrained on the hyperplanes along a reaction path provide information about mechanistic details of a transition pathway. For example, the VIBRAN facility in CHARMM317 can be used to compute the scale-free mechanical coupling network in proteins and nucleic acids.318–321

**7.6. Boxed MD (BXD).** BXD322,323 is a simple technique to estimate rates and PMF G(ρ) along a CV ρ in a single MD simulation. BXD falls within a class of sampling methods such as milestones290,324 where molecular configuration space is divided into a set of boundaries (or hypersurfaces). ρ is kept within a perfectly reflecting “box” for a time interval sufficiently long to reach convergence. This is done by reversing the velocity of the particles involved in the definition of ρ. After a given number of collisions with the boundaries, ρ is allowed to increase or decrease so that a neighboring box can be sampled. From the number of collisions with the boundaries, G(ρ) over the whole range of ρ can be reconstructed, as well as the absolute rate of entering or exiting a specific “box” (Figure 6). Velocity inversion is carried out at each of the boundaries. Assuming that at a certain time the trajectory is in box m, i.e., ρ<sub>m−1</sub> < ρ(r) < ρ<sub>m</sub>, the transition rate from box m to box m + 1

<figure><img src="image_15.png" /><figcaption>Figure 6. In BXD, the range of values assumed by the CV ρ is partitioned in boxes separated by reflective boundaries.</figcaption></figure>

is: $k_{m,m+1} = h_{m,m+1}/t_m$ where $t_m$ is the time the trajectory spends in box m, and $h_{m,m+1}$ is the number of hits (i.e., velocity inversions) at the boundary between $\rho_m$ and $\rho_{m+1}$. After the forward and reverse transition rates are determined, the equilibrium constants between the neighboring boxes m and m + 1 is

$$K_{m,m+1} = \frac{k_{m,m+1}}{k_{m+1,m}} = e^{-\beta \Delta G_{m,m+1}} \quad (30)$$

The free energy $G_m$ can be determined by setting e.g., $G_1 = 0$, and the probability of finding ρ in box m is

$$p_m = \frac{\exp(-\beta G_m)}{\sum_n \exp(-\beta G_n)} \quad (31)$$

which can be multiplied by the normalized probability $P_m(\rho)$ estimated from the histograms within boxes to obtain the probability distribution function $P(\rho) = p_m P_m(\rho)$, or equivalently, $G(\rho) = -k_B T \ln P(\rho)$. In practice, the user sets the position of the boundaries and the number of times the trajectory hits a boundary before it is let into the adjacent one. Both affect the convergence, which can be assessed by performing a single simulation spanning multiple times in both directions over the range of $\rho(r)$ (from the lowest value to the largest, and vice versa). BXD is generalizable to multidimensional CVs using a general velocity-reflection procedure that conserves energy.325

**7.7. Extended Adaptive Biasing Force (eABF) Method.** Adaptive Biasing Force (ABF) is based on estimating the average force acting along a chosen CV, $\xi$, in order to construct and apply a biasing potential $f_m(\xi)$ that augments fluctuations of targeted dynamics.326–328 The classical ABF method is based on TI329 of the average force estimates which are computed in bins along the CV,

$$\Delta A = - \int \langle F_{\xi}(k) \rangle_{\xi} d\xi \\ = \int \left( \left\langle \frac{\partial V(x)}{\partial \xi} \right\rangle_{\xi} - \frac{1}{\beta} \left\langle \frac{\partial \ln|J|}{\partial \xi} \right\rangle_{\xi} \right) d\xi \quad (32)$$

where ΔA is the free energy difference, $\langle F_{\xi}(k) \rangle_{\xi}$ is the average force along $\xi$ in the k-th bin, $V(x)$ is the potential energy function, and $|J|$ is the determinant of the Jacobian.327,330

The biasing force $-\langle F_{\xi}(k) \rangle_{\xi}$ effectively flattens curvatures in the potential energy surface encountered along $\xi$, allowing for extensive sampling of transitions along $\xi$. The biasing potential is adaptive because it is updated by the current estimate of the average force along $\xi$ until convergence.331
Estimating $\langle F_{\xi}(k) \rangle_{\xi}$ brings about complications that hinder the utility of ABF.330 For example, calculating $\partial \ln|f| / \partial \xi$ in Eq. 32 can be challenging.332 The extended ABF method (eABF) was developed to overcome these limitations by introducing an extended potential energy function

$$V_m(x, \lambda) = V_0(x) + \frac{k_{\lambda}}{2}(\lambda - \xi)^2 + f_m(\lambda) \quad (33)$$

where $V_m(x, \lambda)$ is the extended potential energy, $\lambda$ is a virtual particle, and $k_{\lambda}$ is the associated spring constant.333–336 The key distinction of eABF from ABF (Eq. 33) is the extension of the system via the $\lambda$ particle; force estimates are now calculated via Hooke's law and the biasing potential is applied to $\lambda$, which augments transitions in $\xi$ via the harmonic coupling. Since the force estimates come from the harmonic restraint between $\lambda$ and $\xi$, the recovered PMF (along $\lambda$) may deviate from that of the physical system (along $\xi$) depending on the coupling strength. Several estimators have been developed to recover PMF.331,336

As an illustration, a simulation of gas-phase deca-alanine was performed where $\xi$ was defined as the end-to-end distance between the terminal $C_{\alpha}$ atoms. In 500 ns, a number of transitions between the helical state and extended states are realized, with an accompanying PMF along the distance consistent with previous studies (Figure 7).331,337

<figure><img src="image_16.png" /><figcaption>Figure 7. End-to-end distance ($\xi$) over time and the PMF of deca-alanine obtained using eABF.</figcaption></figure>

# 8. ADVANCED ENERGY FUNCTIONS, COARSE GRAINING, AND IMPLICIT MODELS

## 8.1. Multipolar Electrostatics.
Anisotropic charge distributions can be conveniently represented as a superposition of atom-centered multipoles.338 Halogen modifications are a noteworthy example which lead to a $\sigma$-hole on the halogen atom. Such features can be represented by using multipole expansions, often up to quadrupoles.339–343 Multipole-based electrostatics requires introducing local axes to define the orientation of higher-order multipole moments relative to the molecular geometry.

Multipolar interactions have been considered early on in molecular recognition.344 Compared to the spherically symmetric field around a single point charge, atomic multipoles can better capture anisotropic interactions. An example is carbon monoxide that cannot be modeled well with only atom-centered point charges located at nuclear positions of the two atoms because the total charge ($Q=0$) and the total molecular dipole $\mu = 0.048 ea_0$ ($e = 1.6 \times 10^{-19}$ C, the charge of an electron, and $1a_0 = 0.53$ Å, the Bohr atomic length) lead to two opposite partial charges that are small in magnitude. In order to describe its substantial quadrupole moment345–347 ($\Theta = -1.58 ea_0^2$) either a third interaction site halfway between the two atoms is included348 or the two atoms are described by a distributed multipole expansion.338,349–351

The electrostatic potential (ESP) around a molecule can be represented in general as an expansion in multipole moments where the zeroth order contribution arises from atom-centered point charges. Capturing strongly anisotropic and/or directional features, e.g., lone pairs, hydrogen bonding, $\pi$-electron density or $\sigma$-holes352–354 requires a description beyond a single partial charge at each nuclear position. The ESP $\Phi(\mathbf{r})$ is related to the electron charge density $\rho(\mathbf{r})$ through355

$$\begin{aligned} 4\pi\epsilon_0\Phi(\mathbf{r}) &= \int \frac{\rho(\mathbf{r}') d\mathbf{r}'}{|\mathbf{r} - \mathbf{r}'|} \\ &= \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \frac{Q_{lm}}{r^{l+1}} \sqrt{\frac{4\pi}{2l+1}} Y_{lm}(\theta, \phi) \end{aligned} \quad (34)$$

where $\mathbf{r}$ and $\mathbf{r}'$ are spatial variables and $1/|\mathbf{r} - \mathbf{r}'|$ was expanded in powers of $r'/r < 1$ to represent the ESP as a sum over spherical harmonics $Y_{lm}(\theta, \phi)$ from which the spherical multipole moment $Q_{lm}$ is defined as

$$Q_{lm} = \int d\mathbf{r}' \rho(\mathbf{r}') (r')^l \sqrt{\frac{4\pi}{2l+1}} Y_{lm}^*(\theta', \phi') \quad (35)$$

The above can be integrated to yield a compact atom-centered representation of the ESP around a molecule and are used together with the MTPL module of CHARMM.

Alternatively, multipoles of a given order can be represented by fixed charge arrangements, as is done in the distributed charge model (DCM).356,357 It replaces the evaluation of multipole–multipole interactions with the same number of charge–charge terms at the expense of introducing additional charge sites (Figure 8A). The magnitude $q_i$ and position of the

<figure><img src="image_17.png" /><figcaption>Figure 8. (A) 9-charge symmetry-constrained MDCM for $CCl_4$. Red and blue points respectively correspond to negative and positive charge positions. A charge at the C-atom nuclear position is hidden. (B) DFT reference ESP (kcal/[mol·e]) mapped onto the 0.001 au molecular isodensity surface (left); fitted multipolar model truncated at quadrupole (middle); 9-charge MDCM model (right, root-mean-square error (RMSE) of 0.29 kcal/mol over the grid used for fitting).</figcaption></figure>

DCM charges are defined with respect to a reference atom. During fitting it may be useful to constrain the maximum displacement of the DCM charges. Reducing the number of interaction sites can be accomplished using differential evolution to optimize charge positions and magnitudes, and find arrangements that achieve a desired accuracy using a minimal number of charges (MDCM).357 Figure 8B center and right show the ESP of a multipole representation and the
corresponding 9-charge MDCM model for $CCl_4$. More recently, the positions of the MDCM charges were explicitly coupled to the molecular geometry which leads to flexible MDCM (f-MDCM),³⁵⁸ available through the DCM module in CHARMM. It can effectively capture intramolecular polarization, and meaningful atomistic simulations can be carried out for condensed-phase systems. The latest development is kernel-based MDCM (kMDCM) that uses intramolecular separations as the features in a Gaussian-kernel to describe the charge displacements depending on molecular geometry.³⁵⁹

**8.2. Machine-Learning-Based Energy Functions.** Over the past few years, machine learning (ML)-based approaches have flourished for constructing PESs for molecular simulations.³⁶⁰⁻³⁶³ Typical approaches include permutationally invariant polynomials,³⁶⁴,³⁶⁵ neural networks (NNs),⁵⁴,³⁶⁶,³⁷⁰ or kernel-based methods.³⁶⁸⁻³⁷¹ The resulting PESs have been used for gas- and condensed-phase simulations to compute observables including spectroscopic properties and reaction rates.

PES representations based on reproducing kernel Hilbert space (RKHS) have long been used for small molecules.³⁷⁰,³⁷² The TRIAKERN module in CHARMM provides the functionality to use such representations in MD simulations. The kernel coefficients $\alpha$ required for the RKHS-based PES are determined through a versatile external utility.³⁷¹ Typical applications include reactive atom plus diatom collision systems,³⁷³,³⁷⁴ but the method has also been extended to spectroscopic investigations of larger molecules.³⁷⁵

The MLpot module³⁷⁶ in pyCHARMM⁹ allows running of mixed ML/MM simulations. It follows more established hybrid QM/MM strategies in which a usually smaller part of the system is treated with a QM method whereas the larger remainder is represented using an empirical energy function.³⁷⁷⁻³⁸⁰ In ML/MM, a ML representation, for example PhysNet,⁵⁴ is combined via mechanical embedding with an empirical FF such as CGenFF available in CHARMM.³,¹⁷,³⁷⁹

The NN-PES computes the total ML energy and forces together with electrostatic interactions between the predicted fluctuating point charges of the ML-atoms and the static atomic charges of the empirical MM atoms. CGenFF¹⁷ handles energies and forces for the remaining MM atoms and van der Waals interactions between MM and ML atoms. Therefore, a set of van der Waals parameters must be assigned to the ML atoms. In PhysNet, charges of the ML atoms fluctuate depending on solute structure (intramolecular charge redistribution). This is akin to the fluctuating MDCM approach³⁵⁸ where geometry-dependent point charges reproduce the molecular electrostatic potential and describe intramolecular polarization. The advantage of mechanical embedding is the direct application of ML-based models of atomic systems trained in the gas-phase for condensed-phase simulations without additional training. Environment-dependent electrostatics can be included at the training stage by including solvent-surrounded solutes in the training set.

The input file to run pyCHARMM (Section 3.1)⁹ together with MLPot is a Python script. The MLPot module initializes an external model potential and evaluates potential energy and forces for the subset of ML atoms together with the CHARMM FF energy. By adapting the MLPot module in the source code, it is possible to link different model potentials such as ANI³⁸¹ or SchNet.³⁶⁶ If the ML-based PES does not predict atomic charges, the electrostatic contribution between

assigned static point charges of the ML and MM atoms are computed by the empirical energy function.

For chemical reactions that was long a domain of ab initio MD simulation, ML-based energy functions now provide means to run statistically significant numbers of trajectories (~10³ or more) which was previously not possible.³⁸²,³⁸³ More recent examples include malonaldehyde in the gas phase,³⁸⁴,³⁸⁵ double proton transfer in hydrated formic acid dimer,³⁸⁶ or for atmospherically relevant reactions using permutationally invariant polynomials and NN-based energy functions.³⁸²,³⁸⁷⁻³⁸⁹

**8.3. Multipole and Point-Induced Dipole (MPID).** CHARMM now supports advanced electrostatic interactions with multipole expansion up to hexadecapole via the developmental MPOLe module. Different real-space cutoffs can be applied for PME calculations, enabling the use of smaller cutoffs for higher-order multipoles. Additionally, it allows selective exclusion of specific multipole–multipole interactions, a feature utilized for the development of a water model.³⁹⁰ The module also facilitates the calculation of point induced dipoles. Several Thole damping functions, including those used in the AMOEBA³⁹¹ and MPID³⁹² FFs are supported. Moreover, the anisotropic atomic polarizability as utilized in the MPID model has been implemented. Induced dipole moments can be calculated either via full self-consistent field (SCF) relaxation or with an extrapolation scheme that uses weighted average of dipole moments after each of several cycles,³⁹⁵ where the third-order extrapolation with empirically optimized weighting coefficients (OPT3)³⁹⁴ is the default recommended method. See Section 10 for further explanation about CHARMM FFs.

**8.4. Polarizable Intermolecular Potential Functions (PIPFs).** The point-dipole representation of the electronic response of a molecular system to an external field offers a systematic description of polarization effects within a classical framework, where several approaches can be equivalently derived.³³⁸ One example is the PIPF³⁹⁵,³⁹⁶ module of CHARMM.³⁹⁷,³⁹⁸ It complements other polarizable treatments available in CHARMM such as the fluctuating charge³⁹⁹,⁴⁰⁰ and Drude oscillator⁴⁰¹,⁴⁰² models.

In PIPF, each interaction site carries a fixed point charge and an inducible point dipole whose magnitude is determined by the total electric field due to all other point charges and induced point dipoles of the system. Assuming linear response, the induced dipole moment $\mu_i^{ind}$ at center i is proportional to the total electric field ($E_i^{tot}$) typically with a scalar isotropic polarizability $\alpha_i$,

$$ \mu_i^{ind} = \alpha_i E_i^{tot} = \alpha_i [E_i^0 + \sum_{i \neq j} T_{(2)}^{ij} \mu_j^{ind}] \quad (36) $$

The second part of Eq. 36 includes two contributions, the permanent electric field ($E_i^0$) due to fixed-point charges and the induced electric field due to induced dipole moments $\mu_j^{ind}$ at other sites. $T_{(n)}^{ij}$ is the rank-n polarization tensor.³³⁸,⁴⁰⁰

Thole’s interaction dipole (TID) model⁴⁰⁴ is employed in the PIPF model. Although an isotropic polarizability is used (a reasonable approximation for a non-interacting atom), the overall molecular polarizability is anisotropic. It has been shown that with the use of only a single parameter for each atom, the computed molecular polarizabilities using the TID model agree well with experimental data for molecules with a wide range of functional groups.⁴⁰⁴,⁴⁰⁵
Due to the interdependence of $\mu_i^{\text{ind}}$ in Eq. 36, a self-consistent iterative procedure is used in simulations to find them within a given threshold.⁴⁰⁶ The PIPF module in CHARMM includes three complementary approaches. The first is solving Eq. 36 by direct matrix inversion.⁴⁰⁶,⁴⁰⁷ Despite its $O(N^3)$ scaling behavior, exact results provided by direct matrix inversion are important for validating the convergence threshold for the faster iterative approach with $O(N^2)$ scaling. Direct matrix inversion is also needed for handling the intramolecular polarization part in the coupled polarization-matrix inversion and iteration (CPII) method explained below.

The next is to propagate induced dipoles dynamically via an extended Lagrangian.⁴⁰⁸ A Verlet integrator has been implemented to couple the fictitious dipole degrees of freedom to a low-temperature bath using the Nosé-Hoover thermostat.⁴⁰⁹⁻⁴¹² The low temperature dynamics makes the dipole fluctuations close to the true converged results. The extended Lagrangian method accelerates the self-consistent iteration scheme by nearly 2-fold.

The CPII approach involves a preconditioning algorithm. At convergence, the total polarization energy is

$$E^{\text{pol}} = -\frac{1}{2} \sum_{i} \mu_{i}^{\text{ind}} E_{i}^{0} \quad (37)$$

The nuclear gradients can be obtained by differentiating Eq. 37 with respect to atomic Cartesian coordinates.⁴⁰³,⁴⁰⁸ On the other hand, when direct dipole dynamics is used, evaluating energy derivatives invokes additional terms containing $T_{(3)}$⁴⁰⁸ since induced dipoles are not at the variational minimum; those third order terms can be obtained in a compact form in CHARMM.⁴¹³

To evaluate the permanent electric field, standard non-bonded list in CHARMM is used, where intramolecular atom pairs up to 1–3 bonded pairs (connected via a single atom) are excluded and 1–4 bonded pairs (connected via two atoms) are included. To alleviate spurious polarization interactions at short-range, Thole’s second damping function in the TID model⁴⁰⁴ is used by default to determine $E_i^0$ and $T_{(2)}^{ij}$.

Although excluding intramolecular polarization is practical for converging induced dipoles in MD simulations, works based on the TID model suggest that both intramolecular and intermolecular polarization should be included and damped in the same way to obtain consistent molecular polarizability. Unfortunately, including intramolecular polarization between bonded pairs can introduce numerical instability in the convergence of dipoles. While it can be avoided by matrix inversion, it is computationally too expensive for the entire system. The CPII method addresses this³⁹⁸ where iterative convergence of induced dipoles under intermolecular polarization is preconditioned using atom-distributed molecular polarizability tensor obtained from the TID model via matrix inversion of individual molecules. Consequently, Eq. 36 is modified to

$$\mu_K^{\text{ind}} = A_K[\mathbf{E}_K^0 + \sum_{L \neq K}^{M} \mathbf{T}_{(2)}^{KL} \mu_L^{\text{ind}}] \quad (K = 1, \dots, M) \quad (38)$$

where M is the number of molecules, K and L are indices for the corresponding atomic quantities grouped by molecules, and $A_K = [\alpha_K^{-1} - \mathbf{T}_{(2)}^{KK}]^{-1}$ is the atom-distributed molecular polarizability tensor⁴⁰⁴,⁴¹⁴ which is in units of Å³ (see Eqs. 14 and 17 of ref 398). The CPII preconditioning algorithm has

been implemented in CHARMM and it accelerates dipole convergence in liquid simulations of amide and polypeptide systems.

In addition to gradients, second derivatives of PIPFs can be calculated in CHARMM with the point dipole formalism.⁴⁰³ With analytical Hessian available, PIPFs can be used in conjunction with the VIBRan module in CHARMM for vibrational FF analysis.

The PIPF model has been employed in MD simulations to examine polarization effects in a series of organic liquids including alkanes, alcohols, and amides.³⁹⁵⁻³⁹⁷ The results obtained with the classical point-dipole model were found to be in good agreement with those from combined QM/MM simulations in which polarization effects are described quantum mechanically.

Recently, the PIPF model has been employed in the doubly polarized QM/MM (dp-QM/MM) method to enhance the accuracy of SE-QM/MM (SE: semiempirical) free energy simulations.⁴¹⁵ A well-known limitation of SE-QM methods is their tendency to underestimate molecular polarizability compared with experiments and AI/DFT-QM (AI: ab initio) benchmarks, leading to significant errors in free energy profile determined at SE-QM/MM levels. The dp-QM/MM method addresses this by improving the response properties of SE-QM/MM methods through high-level molecular polarizability fitting. Specifically, additional induced point dipoles are introduced on the QM atoms through a set of corrective polarizabilities (“chaperone polarizabilities”), whose magnitudes are determined from ML to reproduce the condensed-phase AI-DFT molecular polarizability along the MEP. These chaperone polarizabilities are then used in PIPF calculations in conjunction with QM/MM to compensate for the polarization energy underestimate in conventional SE-QM/MM simulations. Applied to the Menshutkin reaction in water, the dp-QM/MM method brought the computed and experimental free energy results into closer agreement.⁴¹⁵

### 8.5. Long-Range Lennard-Jones Interactions.
While the PME method³⁵,¹⁷⁶ for evaluating long-range electrostatic interactions was added to CHARMM in 1995,⁴¹⁶ implementing the long-range Lennard-Jones (LJ) interaction lagged considerably. Long-range effects of dispersion are important for accurate calculation of free energies and interfacial properties of liquids and surfactants.⁴¹⁷⁻⁴¹⁹ Ignoring them leads to inconsistencies in the surface tension of lipid bilayers and monolayers.⁴²⁰ To address long-range LJ interactions, a lattice-based method termed LJ-PME⁴²¹ was implemented into CHARMM.⁴²² The name LJ-PME is arguably a misnomer, in that only the C6 ($r^{-6}$) dispersion is calculated with an Ewald summation. The C12 ($r^{-12}$) term continues to be truncated with a standard switching function.⁴ Dispersion-PME has been used interchangeably and is perhaps a better name.

Key to the efficiency of electrostatic PME is the simple multiplicatively separable functional form. While the analogous PME method for dispersion has been available for decades,¹⁷⁶,⁴²³ its adoption has likely been hindered by the fact that many FFs, including CHARMM, use the Lorentz rule

$$\sigma_{ij} = \frac{\sigma_i + \sigma_j}{2} \quad (39)$$

to combine the LJ geometric parameters $\sigma_i$ and $\sigma_j$ for distinct atom types i and j. Since this does not yield a multiplicatively separable form, binomial expansion can be used:
$$ (\sigma_i + \sigma_j)^6 = \sum_{k=0}^{6} \binom{6}{k} \sigma_i^k \sigma_j^{6-k} \quad (40) $$

However, this approach requires several PME evaluations in addition to those used in electrostatic PME. An elegant solution⁴²¹ is to instead assume geometric mean combination rule

$$ \sigma_{ij} = \sqrt{\sigma_i \sigma_j} \quad (41) $$

and apply PME to the dispersion part of the LJ potential, requiring similar effort to that employed in the electrostatic term. This allows the C6 term to be calculated with a single PME calculation involving a simple multiplicative form, $C6_{ij} = \sqrt{C6_i C6_j}$.

To correct for this approximation, the geometric-mean term is analytically subtracted (similar to how 1–2 and 1–3 exclusions are handled in electrostatics) for all pairs within a cutoff distance and substituted with the correct form, which may possibly include NBFIX pair-specific corrections (Section 10). The net effect is that the LJ potential is exact up to the chosen cutoff, beyond which the geometric-mean functional form is used. For an ~8-Å cutoff, the geometric- and arithmetic-mean potentials coincide very closely and the overall approximation is excellent. To account for small but abrupt changes in energy as two schemes handoff at the cutoff, the original formulation in CHARMM applied a shift term to ensure continuity in energy. However this approach does not guarantee continuity in the first derivative. For this reason, a sigmoidal switch function has been implemented to seamlessly transition between the two regimes.⁴²²

The C6 terms are generally computed using the van der Waals parameters in CHARMM. To achieve greater flexibility, these values may be input directly regardless of atom types or short-range C6 values. In this way, specific interactions can be “fixed” to simplify parameter fitting and optimization by treating short-range and long-range dispersion terms independently.

The LJ-PME algorithm provided the impetus to revisit the use of cutoffs in the CHARMM lipid FF (Section 10.4). It has long been known that the high anisotropy present in such systems renders the standard isotropic corrections inappropriate. Because PME makes no assumption about the isotropy of the system, it is well suited to general systems including lipids. Previous iterations of the CHARMM lipid FF were parametrized with a given cutoff, with the cutoff errors implicitly absorbed into the parametrization. While successful, this strategy makes the resulting FF sensitive to the cutoff used at runtime; a choice other than that used for parametrization can yield erroneous results.

Long-range LJ terms can also be evaluated in CHARMM using the Isotropic Periodic Sum (IPS)⁴²⁴ and extended IPS methods.⁴¹⁷,⁴¹⁸ However, LJ-PME provides higher efficiency and transferability of FF among different simulation programs.

## 8.6. FACTS Implicit Solvent.
The Fast Analytical Continuum Treatment of Solvation (FACTS) model is an efficient GB-based implicit solvent method for calculating the solvation free energy of proteins, protein complexes, and protein-ligand interactions.¹³² FACTS is based on the analytical evaluation of the volume and spatial symmetry of the solvent that is displaced from a solute atom by nearby atoms. For each solute atom, these two measures of solvent displacement are combined into an empirical sigmoidal

equation for the calculation of the atomic (or self) electrostatic solvation energy and the solvent accessible surface area (SASA). The former is used to calculate the Born radii in the GB equation. The SASA is used to evaluate the non-polar contribution to solvation.

FACTS is fully analytical and because of its speed, it is useful for MD simulations. It has two main advantages over other implementations of the GB model. First, FACTS does not use the so-called Coulomb field approximation where the electric displacement field of a solute atom is calculated by assuming that the solute-solvent dielectric boundary is spherical with the atom at the center of the sphere. This assumption breaks down particularly for solute molecules with substantial aspherical volume and/or charges located close to the solute-solvent boundary. Second, FACTS does not require setting the dielectric discontinuity surface. Importantly, the CHARMM energy calculation with the FACTS model is only about four times slower than the vacuum energy, and FACTS scales linearly with system size (see Figure 11 in ref 132).

FACTS is versatile as parameters for new or unknown atom types (e.g., non-proteinaceous atoms in organic compounds) can be generated automatically by interpolation from the existing FACTS parameters by the FACTS keyword TAVW. The effect of salt (ionic strength) is treated by the linearized Debye-Hückel approximation.⁴²⁵ The FACTS energy and force terms have been parallelized for multiple CPUs which provides substantial speed-up particularly for large systems. Furthermore, FACTS is compatible with the IMAGE module for periodic boundary condition (PBC) and the BLOCK module for energy decomposition. The INTE command of CHARMM can be used to evaluate the FACTS energy between two groups of solute atoms, e.g., a protein and a small-molecule ligand.

Since the original publication in 2008, FACTS has been employed in many simulation studies of (small) protein folding, (poly)peptide amyloid aggregation, and binding of ligands to proteins. One interesting example is the MD study of the interactions of the toxic Alzheimer's Aβ₁⁻₄₂ peptide with carnosine, the endogenous brain dipeptide β-Ala-His, which revealed salt bridges with charged side chains, and van der Waals contacts with residues in and around the ¹⁷LVFFA²¹ central hydrophobic cluster of Aβ₁⁻₄₂.⁴²⁶ In 2014, the FACTS model was extended to lipid bilayer (membrane) environment by using a position-dependent dielectric constant and an empirical surface tension parameter. It was shown to reproduce the self-energy and pairwise interaction energies in solution calculated by the finite-difference Poisson method.⁴²⁷ However, the FACTS model for the membrane has not been implemented into the official version of CHARMM yet.

The last sentence of the original FACTS paper mentioned potential applications beyond MD: “The accuracy and efficiency of FACTS suggest that it could also be used for protein structure prediction and docking”.¹³² In fact, FACTS has been employed in several docking programs. A recent example is the FASTDock pipeline for the efficient scoring of poses in ortho- and allosteric pockets generated by MD.¹²⁷ Another example is the docking protocol called Attractive Cavities which uses energy minimization and a smoothed potential energy for guiding small molecules into protein cavities. In a successive refinement step, the binding energy is calculated as the sum based on the CHARMM FF and the FACTS model.⁴²⁸,⁴²⁹ Concerning structure prediction, a Python tool allows for the evaluation of the FACTS total energy and its
contribution as descriptors for ML models of protein-protein interactions.⁴³⁰

### 8.7. Implicit Modeling of Membranes.
**8.7.1. Membrane Pores.** The IMM1 implicit membrane model,⁴³¹,⁴³² an extension of the EEF1 effective energy function for soluble proteins,⁴³³ has been adapted to account for aqueous pores.⁴³⁴,⁴³⁵ IMM1 uses two sets of solvation parameters, one for water and one for the non-polar membrane interior. A continuous switching function *f* describes the transition from one environment to the other. Modeling of pores is accomplished by using a switching function *F* dependent on the vertical (*z*) and radial coordinate (the distance *r* from the z-axis). Denoting the thickness of the non-polar part of the membrane as *T* and the pore radius as *R* (Figure 9), the

<figure><img src="image_18.png" /><figcaption>**Figure 9.** Illustration of pore shapes in IMM1. (A) Constant radius. (B) Radius depending on the z-coordinate.</figcaption></figure>

switching functions can be expressed using dimensionless variables:

$$ 
\begin{aligned}
F(z', r') &= f(z') + b(r') - f(z')b(r'), \\
f(z') &= \frac{z'^n}{1 + z'^n}, \quad b(r') = 1 - f(r'), \quad z' = |z|^{\frac{2}{T}}, \quad r' = \frac{r}{R}
\end{aligned}
\quad (42) 
$$

Different pore shapes can be modeled by making *R* dependent on *z*, e.g., $R = R_0 + kz'^2$ (Figure 9B). Because the Gouy-Chapman formulas are invalid for pores in anionic membranes, an alternative approach is based on numerical solution of the Poisson-Boltzmann equation.⁴³⁶ It has been used to investigate the pore forming activity of antimicrobial peptides.⁴³⁵,⁴³⁷⁻⁴⁴² More recently, this model was used for initial evaluation of putative structures of *β* barrel membrane pores formed by fibril-forming peptides and proteins, such as amyloid *β*, IAPP, and *α*-synuclein.⁴⁴³⁻⁴⁴⁶ A similar energy function has been implemented in the Rosetta protein design package.⁴⁴⁷

### 8.7.2. Curved Membranes.
The IMM1 implicit membrane model has been extended to spherical and cylindrical membranes (vesicles and tubes) by changing the definition of the relative depth *z'* from $|z|/(T/2)$ to $|r - R|/(T/2)$ where *R* is the radius of the vesicle or tube and *r* the radial position of an atom.⁴⁴⁸ The model can also account for changes in lateral pressure profile as the membrane bends.⁴⁴⁹ It has been used to study ESCRT-III snf7,⁴⁵⁰ the mechanism of negative curvature generation by IBAR domains,⁴⁵¹ and the interaction of caveolin oligomers with membranes.⁴⁵²

### 8.7.3. Mean-Field Modeling of Deformable Membrane Bilayers via HDGB.
The Heterogeneous Dielectric Generalized Born (HDGB) model⁴⁵³ is an extension of the GBMV implicit solvent model¹⁶⁰ to capture the interaction of biomolecules with biological membranes. The most straightforward approach for implicitly modeling membrane-water interfaces is via a two-dielectric system where the membrane is modeled as a low-dielectric slab ($ε = 1-2$) embedded in a high-

dielectric region ($ε = 80$). This idea was implemented in earlier implicit membrane models.²²⁹,⁴³¹,⁴⁵⁸ The HDGB model refines it by introducing a continuously varying dielectric profile across the membrane-water interface to better describe the actual dielectric profile of membrane bilayers.⁴⁵⁹ The variable dielectric profile is then used in a modified GB equation

$$ 
\Delta G_{\text{solv}}^{\text{GB}} = -166 \sum_{i \neq j}^{N} \left[ 1 - \frac{2}{\epsilon_i + \epsilon_j} \right] \times \frac{q_i q_j}{(r_{ij}^2 + \alpha(\epsilon_i)\alpha(\epsilon_j)e^{-r_{ij}^2/[F\alpha(\epsilon_i)\alpha(\epsilon_j)]})^{1/2}} \quad (43) 
$$

where $166 = (8πε_0)^{-1} (ε_0: permittivity of vacuum)$ is a factor arising from units used in CHARMM, *N* is the number of atoms, $ε_i$ is the variable dielectric profile at the position of atom *i* with charge $q_i$, typically along the *z*-direction coinciding with the membrane normal,⁴⁵³,⁴⁶⁰ $r_{ij}$ is the distance between atoms *i* and *j*, $α(ε_i)$ are Born radii dependent on $ε_i$ and *F* = 8 is a dimensionless parameter. In addition to a variable dielectric profile, the HDGB model introduces a SASA-dependent non-polar contribution:

$$ \Delta G_{\text{cavity}} = \sum_{i=1}^{N} \text{SASA}_i \gamma S(z) \quad (44) $$

where SASAᵢ is the SASA of atom *i*, and *γ* is the surface tension reflecting the strength of the non-polar term, and *S(z)* is an optimizable profile along *z* implemented as a spline-interpolated function. In principle, HDGB can be applied to any heterogeneous dielectric environment including those with non-slab geometries. For example, spherical micelles can be modeled by varying the dielectric and non-polar profiles as a function of the radial position from the center. It is also possible to exclude select atoms from the variable dielectric and non-polar profiles, which allows modeling of membrane channels where membrane-facing atoms would be in contact with water or ions instead of the lipid bilayer.

HDGB variants available in CHARMM are summarized in Table 2. In the original HDGB,⁴⁵³ the dielectric profile was

<table><caption>**Table 2.** HDGB Model Variants Implemented in CHARMM<sup>a</sup></caption><thead><tr><th>Model</th><th>Dielectric Profile</th><th>Non-Polar Profile</th><th>VDW Term</th><th>Deformable</th></tr></thead><tbody><tr><td>HDGB<sup>453</sup></td><td>PB calculations for multilayer dielectric</td><td>Insertion profiles of O<sub>2</sub> and water</td><td>No</td><td>No</td></tr><tr><td>HDGBv2<sup>454</sup></td><td>ΔG<sub>insert</sub> for amino acid analogues</td><td></td><td>No</td><td>No</td></tr><tr><td>HDGBv3<sup>455</sup></td><td>Intramembrane interactions and ΔG<sub>insert</sub> for amino acid analogues</td><td></td><td>No</td><td>No</td></tr><tr><td>HDGBvdW<sup>456</sup></td><td>Same as HDGBv3</td><td>Reoptimized as for HDGBv3</td><td>Yes</td><td>No</td></tr><tr><td>DHDBG<sup>457</sup></td><td>Any</td><td>Any</td><td>No</td><td>Yes</td></tr></tbody></table>

<sup>a</sup>PB, Poisson-Boltzmann; ΔG<sub>insert</sub> insertion free energy of amino acid analogues; VDW, van der Waals.

taken from Poisson-Boltzmann calculations for a probe sphere at different locations in a multilayer dielectric system and the non-polar profile was adjusted to match insertion free energies of model molecules. In HDGBv2,⁴⁵⁴ both dielectric and non-polar profiles were optimized further to improve insertion free energies of amino acid analogues. In HDGBv3,⁴⁵⁵ side-chain
interactions within the membrane from all-atom simulations were taken into account. HDGBvdW⁴⁵⁶ includes implicit van der Waals interactions inside and outside of the membrane separate from the cavity non-polar term. This further improved intramembrane interactions. Furthermore, DHDGB⁴⁵⁷ adds dynamically fluctuating membrane deformations by coupling HDGB to membrane deformation energies from elasticity theory. This allows for a more realistic modeling of charged and polar compounds near and inside the membrane bilayer via membrane deformations. The DHDGB model can be combined with any other HDGB models in Table 2.

With HDGB, it is possible to study a variety of peptide–membrane interactions. A significant advantage of an implicit membrane model is that slow bilayer reorganization kinetics can be avoided. This is especially relevant for peptide–membrane insertion where atomistic simulations may be too slow to converge. HDGB has been used successfully to study the insertion of viral fusion peptides.⁴⁶¹,⁴⁶² Another advantage is that the width of the bilayer can be easily varied simply by scaling the dielectric and non-polar profiles. This allowed a comparison of phospholamban conformational sampling in different physiologically relevant bilayers.⁴⁵⁴ It also led to a method for estimating the optimal membrane width based on the structure of a given integral membrane protein.⁴⁶³

HDGB was also used to estimate the membrane permeability of drug-like molecules,⁴⁶⁴ and as a scoring function for membrane protein structures,⁴⁶⁵ which led to a MD-based structure refinement protocol for integral membrane proteins.⁴⁶⁶ HDGB can also be used for simply simulating the dynamics of membrane-embedded integral proteins.⁴⁶⁷ However, for larger systems, the computational advantage of HDGB over explicit lipid simulations is not as significant.⁴⁶⁸

### 8.8. Transferable Coarse-Graining via PRIMO.
The Protein Intermediate Resolution MOdel (PRIMO) is a CG model for proteins and nucleic acids with resolution intermediate between atomistic and residue levels.⁴⁶⁹ In PRIMO, protein backbones are represented with three particles: C<sub>α</sub>, N, and CO (at the midpoint of the carbonyl group). Non-glycine side chains are represented with one to five beads depending on the size of the side chain (e.g., Figure 10A). Nucleic acids are represented with a similar level of coarse-graining.⁴⁶⁹

The CG sites in PRIMO were chosen to allow an analytical reconstruction of atomistic detail with minimal loss of accuracy by applying known standard bond geometries.⁴⁶⁹ On average, all-atom reconstructions from PRIMO deviate by only 0.1 Å

from the original all-atom models.⁴⁶⁹ For comparison, at the time PRIMO was developed, all-atom reconstructions from C<sub>α</sub>-only models deviated by 1.7 Å⁴⁶⁹ on average and even when C<sub>α</sub> sites were combined with a site at the side chain center, the reconstruction error remained at 0.9 Å.⁴⁶⁹ More accurate all-atom reconstructions from residue-level CG models up to 0.5 Å are now possible with advanced ML models,⁴⁷⁰ but PRIMO has an advantage by maintaining very close connection to atomistic models with significantly reduced number of interaction sites. The near-exact mapping between CG and atomistic levels in PRIMO can be used to compress all-atom trajectory data.⁴⁷¹

The PRIMO FF incorporates a combination of bonded and non-bonded terms in all-atom FF.⁴⁷² It augments standard bonded terms with spline-based functions because interactions between many of the CG sites have multiple minima that are not approximated well with single-well harmonic terms. Moreover, some bonded terms operate on virtual sites that are reconstructed on-the-fly from the CG sites, which is possible because of the computationally efficient analytical mapping from PRIMO to atomistic sites. The virtual site approach is also used for an explicit hydrogen-bonding potential using 2D spline interpolations of PMF as a function of hydrogen bond angle and distance similar to the CMAP torsion potential in CHARMM.⁴⁷³ Finally, PRIMO captures solvation effects via implicit solvent terms. A GBMV-based model¹⁶⁰ is used to capture electrostatic contributions to the solvation free energy. It is complemented with a per-residue SASA term that adds non-polar contributions and compensates for incomplete electrostatic solvation contributions with the GB model because of less polarized PRIMO interaction sites. By replacing GBMV with the HDGB implicit membrane model,⁴⁵³ PRIMO can be extended to simulate protein–membrane interactions.⁴⁷⁴ The PRIMO FF were parametrized primarily by matching energies from CHARMM’s all-atom FFs, in particular CHARMM22/CMAP⁴⁷⁵ and CHARMM36⁴⁷⁶ with further adjustments made based on simulations of test peptides.⁴⁷²

The PRIMO FF is fully transferable to other systems.⁴⁷⁷ It is possible to run stable MD simulations of arbitrary protein systems without restraints⁴⁷²,⁴⁷⁴ and PRIMO can be used in combination with enhanced sampling techniques in peptide folding simulations,⁴⁷² or to study the insertion of peptides into membranes.⁴⁷⁴ PRIMO is also useful for extensive conformational sampling in protein structure refinement.⁴⁷⁸

Because PRIMO and atomistic interaction potentials are similar and compatible with each other, it is possible to run hybrid multiscale simulations where parts of a system are represented by PRIMO whereas other parts are represented in atomistic detail. One example is the simulation of a peptide in atomistic detail surrounded by crowder molecules represented at the CG level using PRIMO,⁴⁷⁹ where the coupling between the CG and atomistic levels involves only the non-bonded and solvation terms. It is also possible to run multiscale simulations where different parts of the same molecule are represented either atomistically or via PRIMO.⁴⁸⁰ The coupling between the all-atom and CG levels extends to the bonded terms by maintaining dual resolution across the interface between CG and atomistic regions. PRIMO resolution is trivially obtained from the atomistic level whereas atomistic sites are reconstructed analytically from the PRIMO model (Figure 10B).⁴⁸⁰ Using this approach, it is possible, for example, to efficiently sample dynamic regions of a given system at the CG

<figure><img src="image_19.png" /><figcaption>**Figure 10.** Illustration of the PRIMO CG model. (A) PRIMO interaction sites (red spheres) for asparagine as an example. (B) Hybrid all-atom/CG model of a protein with coupling between a PRIMO region and an atomistic region.</figcaption></figure>
level while applying atomistic detail to maintain accuracy of more conserved structural elements.⁴⁸⁰

PRIMO is unique because of its close coupling to the CHARMM all-atom FF. Its advantage is a high degree of transferability compared to other CG models and a suitability for multiscale simulation approaches where a given system can be represented simultaneously at different levels of resolution. The CG nature of PRIMO improves computational efficiency over comparable all-atom simulations with the same GBMV implicit solvent model by about a factor of 10 in wall-time.⁴⁷² There are additional gains in efficiency due to the smoother energy landscape at the CG level. However, the use of the relatively expensive GBMV model limits PRIMO's overall performance, especially for larger systems where the advantage of implicit solvent over explicit solvent diminishes.⁴⁶⁸ Combining PRIMO with other less expensive implicit solvent models in CHARMM is an option for potentially overcoming these limitations.

## 9. SPECIALIZED RESTRAINT METHODS

Restraint energy functions apply biases on particular degrees of freedom, and can drive the system into conformational states that may be otherwise inaccessible, thereby improving sampling around those states.

**9.1. CONSHELIX Module.** Most restraint potentials control reaction coordinates between atoms, such as CONS HARM, NOE, and RESD for atoms, and CONS DIHE for dihedral angles in CHARMM. By comparison, the restraint potentials in the CONSHELIX module are applied to molecular-level reaction coordinates such as helices and hairpins (Figure 11). These energy functions are especially

<figure><img src="image_20.png" /><figcaption>**Figure 11.** Helical reaction coordinates.⁴⁸¹,⁴⁸² The tilt angle $\tau$ is relative to the z-axis; the rotation angle $\rho$ is measured for a designated atom about the helical axis; the bend angle $\theta$ is between two helical axes; the minimum distance $D$ and crossing angle $\Omega$ are between two neighboring helices. Corresponding restraint energy functions are handled by the CONSHELIX module in CHARMM. For a hairpin, tilt, rotation, and distance restraint potentials are also available.</figcaption></figure>

useful for controlling motions of transmembrane domains within the lipid bilayer,⁴⁸³,⁴⁸⁴ which are critical for, e.g., signal transduction,⁴⁸⁵,⁴⁸⁶ transport of ions and small molecules, antimicrobial activity, and transmembrane responses.⁴⁸⁷

The CONSHELIX restraint potential $U_{\xi}(\mathbf{R})$ (Eq. 45) takes a quadratic form for non-periodic variables ($\xi = \tau, \theta, D$, or $\Omega$; defined in Figure 11) and a cosine function for periodic variables ($\xi = \rho$), where $\mathbf{R}$ represents coordinates of atoms selected to define the helix/hairpin principal axis. Denoting the force constant and the target value as $k_{\xi}$ and $\xi_0$, respectively,

$$ U_{\xi}(\mathbf{R}) = \begin{cases} \frac{1}{2} k_{\xi} [\xi(\mathbf{R}) - \xi_0]^2 & (\text{non-periodic}) \\ k_{\xi} [1 - \cos(\xi(\mathbf{R}) - \xi_0)] & (\text{periodic}) \end{cases} \quad (45) $$

An example application is the study of the mismatch in the thickness of the hydrophobic region between the protein and the lipid bilayer,⁴⁸⁷ which leads to changes in lipid length, tilting of transmembrane proteins, and association of multiple transmembrane proteins. To explore the resultant tilting motion of the WALP19 model helix peptide (sequence: GWW(LA)₆LWWA) in a DMPC bilayer, umbrella sampling was performed.⁴⁸⁸,⁴⁸⁹ The weighted histogram analysis method (WHAM) and TI were utilized to obtain the PMF as a function of the tilt angle $\tau$, where 'precession entropy' was proposed as the driving force for tilting. As $\tau$ increases, the accessible volume of the helix conformation also increases (precession entropy), which stabilizes the helix in a tilted orientation.

Another application is the association free energy between two helices using the distance $D$ in a DMPC bilayer.⁴⁹⁰ The asparagine residue at the center of the helix drives helix-helix association in a membrane by forming bifurcated hydrogen bonding. However, the interaction between the helix and the lipid promotes dissociation of the two helices. The free energy cost for lipid depletion between helices is relatively small, around 7.6 cal/[mol·Å³], compared to that for cavity formation in water, 24–33 cal/[mol·Å³]. Through TI, various energy terms within a residue can be decomposed, which inform residues that favor helix-helix or helix-lipid interactions.

**9.2. SSNMR Module.** Solid-state NMR (SSNMR) spectroscopy is used to determine membrane protein structures in a native-like membrane environment. It utilizes 2D $^{15}$N-$^{1}$H NMR polarization inversion spin exchange at magic angle (PISEMA) spectrum experiments to obtain orientational information from dipolar coupling and chemical shift. To harness this experimental data effectively, the SSNMR module has been implemented in CHARMM.⁴⁹¹

Experimental observables ($O_{ssnmr}$) obtained through the SSNMR spectroscopy are dipolar coupling ($\nu$) for $^{15}$N-$^{1}$H pair and chemical shift ($\sigma$) for $^{15}$N atom:

$$ \begin{aligned} \nu &= \frac{\nu_0}{2}(3 \cos^2 \theta - 1) \\ \sigma &= \left\langle \sum_{i=1}^{3} [\hat{n} \cdot \hat{e}_i(t)] \sigma_{ii}(t) [\hat{e}_i(t) \cdot \hat{n}] \right\rangle \end{aligned} \quad (46) $$

Here, $\nu_0 \equiv \gamma_N \gamma_H h \mu_0 / (8 \pi r^3)$ is the dipolar coupling constant, where $\gamma_N = -2.71 \times 10^7/(\text{T} \cdot \text{sec})$, $\gamma_H = 2.675 \times 10^8/(\text{T} \cdot \text{sec})$, $h$ is the Planck constant, $\mu_0$ is the vacuum permeability, and $r$ is the N-H bond length. $\theta$ is the angle between the N → H bond vector and the direction $\hat{n}$ of the magnetic field that is assumed to be normal to the membrane. $\hat{e}_i$ ($i = 1, 2, 3$) are the instantaneous basis vectors for the chemical shift tensor where $\sigma_{ii}$ is its diagonal element. $\hat{e}_1$ and $\hat{e}_3$ are on the plane spanned by C, H, and N atoms of the peptide bond, and $\hat{e}_2 = \hat{e}_3 \times \hat{e}_1$ 491–493

Due to $\cos^2 \theta$ in Eq. 46, for one dipolar coupling experimental value, up to four helix orientations are possible, potentially with multiple helical structural conformers in membrane ('structural ambiguity'). The simplest form of restraint energy $U$ is a quadratic function that minimizes experimental values
$$U = \begin{cases} \sum_{i}^{M} k_{dc} (\nu_{i} - \nu_{exp,i})^{2} & (\text{dipolar coupling}) \\ \sum_{i}^{M} k_{cs} (\sigma_{i} - \sigma_{exp,i})^{2} & (\text{chemical shift}) \end{cases} \quad (47)$$

where M is the number of restraint potentials and the subscript 'exp' denotes experimental value. Eq. 47 can be minimized through various methods such as MD, MC, genetic algorithm optimization, and simulated annealing. A structure that satisfies experimental data can thereby be obtained. An example is modeling the structure of the fd-coat, the major pVIII coat protein of the fd filamentous bacteriophage (Figure 12). It has

<figure><img src="image_21.png" /><figcaption>Figure 12. Structural changes of the fd-coat membrane protein through MD simulations with the SSNMR restraint potential. (A) The fd-coat system. IP (blue): In-plane helix; TM (red): Trans-membrane helix. (B) Distribution of chemical shift (x-axis) and dipolar coupling (y-axis) calculated for the initial state structure. (C) Structure and calculated distribution in an intermediate stage of the simulation. Note changes in the distributions and conformations. (D) The final stage of the simulation. IP is at the lipid-water interface, while TM is positioned within the membrane.</figcaption></figure>

two types of α-helices within the lipid membrane: the amphipathic in-plane (IP) helix at the water-lipid interface, and the longer transmembrane (TM) helix. Solution and solid-state NMR structures are available (PDB ID: 1FDM and 1MZT).⁴⁹⁴ Calculations were performed in vacuum while considering a virtual lipid environment where Eq. 47 was used with experimental values as minima. Initially, the structure is perpendicular to the lipid membrane. Subsequently the IP region adheres to the lipid membrane surface, and the TM region adopts a tilted structure to match the hydrophobic length of TM and the membrane thickness (Figure 12), which satisfy constraints based on experimental data.

Other applications of the SSNMR module on trans-membrane proteins include MerF (a mercuric ion transporter), M2 (transmembrane domain from influenza A virus), and Vpu (viral protein u from HIV-1).⁴⁹¹ In another application, the SSNMR module was combined with the ensemble dynamics (ED) technique⁴⁹⁵⁻⁴⁹⁸ in an explicit membrane system to address discrepancies between semistatic and dynamic fitting models of SSNMR observables. Compared to these two fitting approaches, the main advantage of the SSNMR ED is its ability to generate an ensemble of structures (e.g., TM helix orientational distribution) that satisfy experimental observables within a reasonable physical model, without prior knowledge about the underlying distribution or motion.

### 9.3. Residual Dipolar Coupling (RDC) NMR Orientational Restraint.
The RDC module in CHARMM leverages experimental time-averaged RDC orientational NMR restraints (Figure 13).⁴⁹⁹ RDC informs about the orientation of each

<figure><img src="image_22.png" /><figcaption>Figure 13. From RDC ($D_{PQ}$) measurement to structure calculation using CHARMM illustrated by using the N-H internuclear vector of protein G (PDB 1P7E).⁵⁰⁰ The protein solution in the NMR tube is indicted in green. Oval-shaped slabs indicate the alignment media. $B_0$ is the applied static magnetic field, and $A_i$ and $M_i$ ($i = 1, 2, 3$) represent the 3 principal axes for the alignment tensor and the inertia tensor of the molecule, respectively.</figcaption></figure>

internuclear vector $r_{PQ} = r_P - r_Q$ formed by a pair of NMR active nuclei P and Q in a molecule with respect to the static magnetic field $B_0$. CHARMM uses decoupled RDC orientational information⁵⁰⁰ that consists of (1) the angle $\psi_i$ between $B_0$ and the i-th principal axis of the inertia tensor of the molecule, $M_i$, and (2) the angle $\theta_i$ formed between $r_{PQ}$ and $M_i$. The RDC between the two nuclei is

$$D_{PQ} = \frac{2 D_{\text{const}}}{3 r_{PQ}} \sum_{i,j=1}^{3} \left\langle \frac{3}{2} \cos \psi_i \cos \psi_j - \frac{1}{2} \right\rangle \times \left\langle \frac{3}{2} \cos \theta_i \cos \theta_j - \frac{1}{2} \right\rangle \quad (48)$$

where $D_{\text{const}} = -S\mu_0\gamma_p\gamma_Qh/8\pi^3$, where $\mu_0$ is the magnetic permeability of vacuum, $\gamma_p$ and $\gamma_Q$ are gyromagnetic ratios of the nuclei P and Q, h is the Planck constant, and S is the generalized order parameter that describes the internal motion of the internuclear vector. The angular brackets indicate time average. Eq. 48 can be expressed in terms of the alignment tensor or Saupe order matrix A

$$D_{PQ} = \frac{D_{\text{const}}}{r_{PQ}^5} \text{Tr}[AO^T R_{PQ} O] \quad (49)$$

where O consists of the three principal axes (eigenvectors) of $M_i$ and the 3 × 3 matrix $R_{PQ} = r_{PQ} \otimes r_{PQ}$. Since A is traceless, only 5 components are independent, which can be determined by using singular value decomposition (SVD) with the aid of experimental RDC,⁵⁰¹ $D_{n \times 1}^{\text{exp}}$ ($n$: number of experimental values) and $M_i^{500}$

$$\hat{A}_{5 \times 1} = V_{5 \times 5} [1/W]_{5 \times 5} U_{5 \times n}^T D_{n \times 1}^{\text{exp}} \quad (50)$$

where V, 1/W, U and D are matrices arising from SVD.

CHARMM supports simultaneous use of RDCs measured between different NMR active nuclei (e.g., N-H, Cα-C, and C-N). Since RDC provides only the orientational information
tion, it alone cannot fully determine the structure of a
biomolecule. Thus, RDC and other NMR restraints such as
NOE (via the NOE keyword) are used together. Alternatively,
one can use a fragment-based approach⁵⁰² to determine
molecular structure using multiple RDCs collected in multiple
alignment media. It is recommended to “reset” the RDC
module before using it for structure determination.

The module requires the name of the input file(s) containing experimental RDC values, their upper and lower bounds, and the corresponding internuclear atoms. Since the CHARMM RDC module uses the harmonic potential with soft asymptotic behavior, for each RDC, one can specify the force constant for the harmonic potential and the slope for the asymptotic function (default = 1), the value for the exponential function used in the soft asymptote if it is other than 1, and the cutoff length for the harmonic function (default: 1 Å). The input files may contain RDCs collected in different alignment media, stored in the CHARMM (default), BMRB, or XPLOR format. Other command options include the maximum RDC restraints, whether the RDCs (other than for N–H) need to be scaled with respect to the N–H RDC, and whether the principal axes are calculated only with respect to the RDC atoms. For validation, users can compare experimental versus calculated RDC data and back-calculate other types of RDC (e.g., C–N and/or Cα–C) for the structure based on the input RDC used (e.g., N–H).

**9.4. Torque Application.** Application of a torque on selected atoms about a user specified axis can be achieved by the PULL TORQue command.⁵⁰³

10. CHARMM FORCE FIELD DEVELOPMENT

Since the publication of the second CHARMM paper in 2009,³
significant progress has been made in the CHARMM-related
FFs. While the additive CHARMM36 (C36) FF was mature at
that time, additional enhancements and refinements, including
creation of CGenFF, were performed. With respect to
polarizable FFs, major developments in the classical Drude
oscillator model were made including coverage of all the major
classes of biomolecules and progress was made toward a Drude
General FF (DGenFF). A detailed description of the potential
energy functions for the additive and Drude FFs is presented in
ref 23.

In addition, a polarizable model based on MPID (Section 8.3) onto which Drude FF-based parameters can be mapped was presented.³⁹² To avoid the costly SCF procedure (Section 8.3), a small mass is assigned to the Drude particles, which are then propagated as dynamic variables during simulations via a dual-thermostat extended Lagrangian algorithm, with a “cold” temperature imposed on the degrees of freedom corresponding to the induced dipole. The statistical mechanical validity of this procedure was clarified.⁵⁰⁴ In addition to CHARMM, the additive and Drude FFs are available in OpenMM,⁵⁰⁵ facilitating the CHARMM/OpenMM API (Section 2.1), as well as in GROMACS⁵⁰⁶ and NAMD,⁵⁰⁷ though the implementation of the Drude FF in GROMACS is currently limited, and NAMD currently does not include LJ-PME capabilities (Section 8.5).⁴²¹ Notable is the ability to setup and generate inputs for complex molecules for a range of programs using CHARMM-GUI⁸⁷ (Section 3.3), including for the Drude polarizable FF.¹¹⁵ An advantage of the Drude FF over other polarizable FFs is the computational efficiency, i.e., it is only about 4 time slower to run compared to additive FFs.

**10.1. Water, ions, and Polar Solvents.** A number of extensions of the Drude FF, including additions to the potential energy function, have been made with respect to water and ions since 2009. Extensions to the additive FF include the use of alternative LJ interactions between the CHARMM TIP3P water model and proteins¹⁴⁴ (Section 10.2), and revised LJ parameters for Na⁺ and Ca²⁺ with lipids.⁵⁰⁸,⁵⁰⁹ The default water model with the Drude FF is the SWM4-NDP (simple water model with 4 sites and negative Drude polarization) model,⁵¹⁰ though a 6-point model SWM6 with an improved condensed phase hydrogen bonding properties is available with an increased computational cost.⁵¹¹ The energy function was expanded to account for polarization anisotropy to describe the dielectric constant of liquid amides more accurately.⁵¹² Substantial work was undertaken on the ions in the Drude FF, including a range of monatomic ions⁵¹³,⁵¹⁴ and later for molecular ions, including a number of ions uncommon in biomolecular systems.⁵¹⁵⁻⁵¹⁷

Developing the parametrization for all charged moieties made it possible to adopt a consistent absolute solvation scale for monatomic and molecular ions.⁵¹³,⁵¹⁵ Notable with the optimization of Mg²⁺ was the use of an LJ repulsion between the SWM4 Drude particle and the Mg²⁺ ion.⁵¹⁷ This enabled steric repulsion between water and ion as the water polarizes in the electric field of the ion. This avoids overbinding and yields a model that can reproduce both the experimental thermodynamic and kinetic properties, a capability not attained in any other FF to our knowledge. It highlights the importance of explicitly including electronic polarizability in a FF as well as the advantage of using a Drude oscillator particle in representing the electronic degrees of freedom. It also captures the cation-π interactions for aromatic side chains more accurately.⁵¹⁸ As parametrization of ions and charged molecules is typically carried out to reproduce experimental data in the infinite dilution limit, accounting for the osmotic pressure made it possible to extend the models to concentrated solutions via pair-specific LJ (NBFIX) and Thole electric shielding (NBThole).⁵¹⁹ A similar philosophy was exploited to optimize the parametrization of amide solutions.⁵²⁰ Early efforts led to a preliminary set of optimized parameters for ion-protein interactions,⁵²¹ although additional tests revealed a number of issues that are currently being addressed.

**10.2. Proteins.** Developments in the C36 additive protein FF have involved two iterations building on the C27 FF, also known as C22/CMAP.⁴⁷³ The revisions primarily involved optimization of the bonded parameters with only minimal changes in non-bonded terms.⁴⁷⁶ The first iteration published in 2012 yielding C36 focused on the CMAP term targeting NMR solution data for non-Gly, non-Pro amino acids, the CMAP terms for Gly and Pro residues targeting high-level QM data and the χ₁ and χ₂ side-chain dihedral parameters, with the latter targeting condensed phase data from simulations of (Ala)₄-X-(Ala)₄ model peptides.⁵²² Subsequent focus was on the CMAP term to account for oversampling of the αL conformation with C36 and improvements in the interaction between Arg and carboxylate groups using an off-diagonal LJ term (NBFIX in CHARMM nomenclature), yielding C36m.¹⁴⁴ These additional optimizations lead to improved treatment of intrinsically disordered peptides (IDPs) while maintaining accurate treatment of folded proteins. C36m is considered the default FF for additive protein simulations.
An interesting outcome of that study regarded the role of LJ
interactions between the water molecule and the protein in the
sampling of folded versus unfolded states of IDPs. While
changing the magnitude of the water–protein interactions
could improve the equilibrium between folded and unfolded
states for a specific protein, a general solution that universally
treats all IDPs in the context of an additive FF may not be
accessible. Additional works included improved treatment of
cation–π interactions,⁵²³,⁵²⁴ and halogen–protein interactions
important for ligand–protein simulations.⁵²⁵ Furthermore, the
additive FF was extended to over 100 non-standard amino
acids⁵²⁶ and to α-methyl amino acids.⁵²⁷

Advances in the Drude FF have been substantial as prior to
2009 only water, ion, and small molecule parameters had been
published. Parameters have been released for many other
biological molecules, including significant work on the protein
portion of the FF. Small molecule developments included
heteroaromatics,⁵²⁸ sulfur-containing compounds,⁵²⁹ and
ethers⁵³⁰ with the Drude FF shown to yield accurate hydration
free energies facilitated by the use of atom-pair-specific LJ
parameters (i.e., NBFIX).⁵³¹ Building on the foundation of the
small molecule parameters the first generation of the Drude
protein FF, Drude-2013, was presented. It overcame challenges
with moving from individual molecules to a polymer in a
polarizable FF, where unexpected overpolarization was avoided
by accounting for the conformational properties of the
polypeptide backbone.⁵³²

Application to polypeptide simulations revealed importance
of explicit polarization in both peptide folding<sup>533</sup> and
unfolding,<sup>534</sup> by capturing cooperativity inaccessible to
additive FFs. A number of issues found in Drude-2013
including the stability of β-sheet structures, led to additional
optimization yielding Drude-2019.<sup>535</sup> Improvements involve
both the polypeptide backbone and side-chain conformational
properties including optimization of the electrostatic parame-
ters of the atoms linking side chains to the backbone, as well as
the treatment of cation-π and anion-π interactions.<sup>536</sup>
Drude-2019 shows systematic improvements as compared to
C36m and Drude-2013, and it allows stable simulations of
proteins on the microsecond time scale. While the Drude FF
has largely been developed assuming an explicit solvent model,
a Poisson-Boltzmann (PB) implicit solvation model has been
developed<sup>537</sup> and subsequently used to predict pK<sub>a</sub> values of
proteins.<sup>538</sup> With the Drude PB model, pK<sub>a</sub>'s calculated for 8
proteins were insensitive to the assigned dielectric constant, in
contrast to the need for a value of 4 with C36m. This indicates
a potential advantage of the polarizable model in implicit
solvent approaches.

**10.3. Nucleic Acids.** Both additive and polarizable FFs for nucleic acids have been updated. Updates to the C27 nucleic acid FF include improved treatment of RNA, largely to account for contributions of the 2'OH group to conformational heterogeneity⁵³⁹ as reported in a combined QM/bioinformatics study.⁵⁴⁰ Work on DNA focused on the equilibrium between the BI and BII conformations in duplex structures.⁵⁴¹ In both cases, adjustments to the C36 FF only involved select dihedral parameters, suggesting that minimal improvements within the non-polarizable additive approximation were necessary. Beyond canonical DNA and RNA, the C36 FF was extended to a range of naturally modified ribonucleotides as required for the ever increasing list of non-coding RNAs.⁵⁴²

Given the high charge density of polyanionic nucleic acids, a
polarizable model is of particular interest. Development of the
Drude FF was based on carbohydrate, ion, and heteroaromatic
parameters. The first step involved optimization of Drude
parameters for nucleic acid bases targeting a range of QM data
for interactions with water and for base-base interactions, as
well as experimental data including base crystal geometries and
heats of sublimation.⁵⁴³ They were then combined with initial
parameters for the phosphodiester linkage yielding the Drude-
2013 DNA model³⁴⁴,⁵⁴⁵ that was iteratively optimized with
particular emphasis on dihedral parameters associated with
phosphodiester, sugar, and sugar-base glycosidic linkages. The
optimization involved comparison with experiments including
crystal data, to ensure suitability in simulations of duplexes in
the condensed phase such as the equilibrium between A- and
B-DNA and the BI/BII forms. This initial model improved
agreement with experimental data regarding base flipping,⁵⁴⁶
and yielded insights into the distribution and competition
between ions around duplex DNA.⁵⁴⁷ Impact of ions on DNA
conformation including the minor groove width could also be
addressed.⁵⁴⁸ Such results cannot be captured well by additive
FFs, again emphasizing the utility of the polarizable model for
studying charged species.

Subsequent optimization of the Drude-2013 DNA FF
focused on the underestimation of base stacking in duplexes
and unwinding of Z-DNA. It involved additional QM
calculations on Z conformations and application of higher-
level model chemistries for other QM data.⁵⁴⁹ The resulting
model reproduced both crystal and solution scattering data
over a range of duplexes in microsecond simulations.⁵⁵⁰ The
FF was also extended to RNA⁵⁵¹ which focused on the role of
the 2'OH group on the conformational properties using QM
data on RNA-specific model compounds. Condensed phase
testing involved stem-loop structures, adenine riboswitch, and
canonical duplexes, showing good agreement with crystallo-
graphic and NMR data.

The combination of the DNA and RNA FF, termed Drude-
2017, was applied successfully to a number of systems
including quadruplexes where the ions in the G tetramer are
stabilized by the explicit inclusion of electronic polar-
izability.⁵⁵² However, a tendency of the Drude-2017 FF to
overpolarize the Drude particle during MD simulations was
noted. While this was addressed by using the Drude hardwall
constraint⁵⁵³ it represents a non-adiabatic condition. The
electrostatic parameters were subsequently adjusted to yield a
model that was successfully used in simulations of RNA
hairpins.⁵⁵⁴

**10.4. Lipids.** Significant advances were made to the
additive and polarizable lipid FFs. Revised parameters for 6
neutral lipids were introduced, yielding the C36 lipid FF.⁵⁵⁵ As
background, previous CHARMM lipid FF⁵⁵⁵ required an
applied surface tension to avoid unphysical bilayer surface area
contraction in NPT simulations. Adjustments to charges and
torsion angles in the headgroup region in C36: (1) reduced the
surface tension to zero at the observed experimental surface
area per lipid for free-standing dipalmitoylphosphatidylcholine
(DPPC) bilayers; (2) increased area compressibility to
experimental ranges; and (3) captured the experimentally
observed splitting in deuterium order parameters for carbons
in glycerol and carbon 2 of the chain 2. The C36 FF was
further validated by the agreement with experimental bending
constants⁵⁵⁶ and spontaneous curvatures⁵⁵⁷,⁵⁵⁸ (those for
bilayers required new code for pressure profiles described in
**Section 12.3).** Extensions to new lipids are ongoing. To date, common phospholipids are parametrized, including 13 variants of inositol lipids, sphingolipids, 5 hydroxylations for ceramide lipids, ether lipids, glycolipids, and acyl chain variants (saturated, monounsaturated, polyunsaturated, branched, and cyclic). Excluding the nearly unlimited variations in glycolipids and lipopolysaccharides, over 300 lipids have been parametrized for C36⁵⁵⁹⁻⁵⁶¹ and they are readily available in CHARMM-GUI. A united-atom representation wherein hydrogen atoms are combined with their bonded heavy atom, has also been formulated and tested for most common lipids (C36UAr).⁵⁶² It is currently being extended to other lipid head groups and chain types such as sphingolipids.⁵⁶³

Despite its extensive refinements and wide usage, C36 has two fundamental limitations: sensitivity of the truncation method used for the LJ interactions, and lack of polarizability. The former manifests as inconsistent bilayer and monolayer surface tensions, which is because bilayers are parametrized to agree with their experimental surface area at particular temperatures. The acyl chain-air interface of monolayers is highly sensitive to truncation of the LJ potentials, causing underestimation of surface tension when using the same parameters that otherwise yield accurate results for bilayers. Conversely, the surface tension of the C36 DPPC monolayer agrees well with the experimental value when long-range LJ terms are included, but the bilayer contracts.⁴²⁰ It was thus necessary to parametrize the bilayer and monolayer consistently. While it can in principle be carried out with a truncated LJ potential, it is physically more reasonable to parametrize both with long-range interactions, which also avoids the sensitive dependence of bilayer properties such as for phase changes, on user-specified cutoff values. While a computationally efficient way of including long-range LJ terms in anisotropic systems such as bilayers and monolayers was unavailable when C36 was developed, subsequent incorporation of LJ-PME⁴²¹,⁴²² (Section 8.5) led to reparameterizing C36 to C36/LJ-PME.⁵⁶⁴,⁵⁶⁵ Consistency of bilayer and monolayer surface tensions for DPPC in C36/LJ-PME was obtained without compromising the overall quality of C36 for bilayers. Yet, monolayer isotherms at very large surface area where the surface tension of water-air is important are not well-described in C36/LJ-PME because the water-air surface tension of TIP3P (the default water in the additive CHARMM FF) is substantially lower than the experimental value.⁴¹⁷

The second limitation of C36, the lack of polarizability, manifests as water permeability in saturated lipids being 5-fold lower than experimental values.⁵⁶⁶ This is because the transfer free energy of water into hexadecane (a good model for the interior of a bilayer) is overestimated by 1 kcal/mol. Also, the dipole of the additive water cannot readjust when it is in the lipid environment. This motivated the development of the CHARMM Drude polarizable FF. Early versions of the lipid Drude FF⁵⁵³,⁵⁶⁷,⁵⁶⁸ provided insight into membrane dipole potentials,⁵⁶⁷ the mechanism of permeation of arginine as a function of membrane thickness,⁵⁶⁹ as well as ion conduction along the narrow gramicidin A channel.⁵⁷⁰ While these studies demonstrated the importance of a polarizable FF for membranes, the initial parametrization of phospholipid molecules had a number of shortcomings, including overestimated bilayer area compressibility. Furthermore, it was optimized with a truncated LJ potential without accounting for long-range dispersion. Bilayer surface areas and compressibility of the revised Drude-2023 FF⁵⁷¹ agree much better with

experiments. More importantly, Drude-2023 yields more accurate dipole potentials, water permeability, monolayer isotherms, and lipid diffusion constants compared to C36 or C36/LJ-PME. Efforts toward a more comprehensive collection of lipids in the context of Drude FF, including charged lipids and ceramide-based lipids, are ongoing.

**10.5. Carbohydrates.** Building FF for carbohydrates poses a particular challenge given the wide range of monosaccharides, including both furanoses and pyranoses, the large number of chemical functional groups beyond hydroxyls, and various glycosidic linkages in poly- and oligosaccharides. Additive carbohydrate FF developments since 2009 included acyclic species⁵⁷² and furanoses⁵⁷³ along with the required glycosidic linkage.⁵⁷⁴ Significantly increasing the coverage of the FF was the inclusion of a variety of chemical groups along with testing on polysaccharides and glycan-protein interactions.⁵⁷⁵,⁵⁷⁶ They together represent the carbohydrate portion of the C36 FF that has been widely used for carbohydrates, glycolipids and glycoproteins.

A similar path was taken with the Drude polarizable FF. Extensive non-bonded parameter optimization was undertaken on acyclic polyalcohols,⁵⁷⁷ aldehydes, and ketones,⁵⁷⁸ as required for the treatment of linear alcohols. FF for furanose⁵⁷⁹ and pyranose⁵⁸⁰ monosaccharides were completed and subsequent adjustments were made to the LJ parameters of pyranoses to improve the treatment of stacking interactions that led to better diffusion behaviors of glucose.⁵⁸¹ This was followed by parametrization of glycosidic linkages involving furanoses and pyranoses⁵⁸² and extension to N-acetyl groups⁵⁸³ and both N- and O-linkages for glycoproteins.⁵⁸⁴ Application of both the C36 and Drude FF to mannose disaccharides showed good agreement with NMR observables.⁵⁸⁵ As with the rest of the Drude FF, the nomenclature, with few exceptions, has been designed to be identical to that of C36 for ready access.

**10.6. Small Molecules.** CGenFF and DGenFF have been developed to greatly broaden the coverage of FFs by rapidly generating topologies and parameters for a wide range of molecules including those of medicinal chemistry and ionic liquids. CGenFF initially leveraged the wide collection of topologies and parameters of C36. Its coverage then extended to drug-like molecules by applying an optimization protocol that maintains compatibility with C36. While the initial CGenFF paper focused on the general philosophy of the model and details of parameter optimization,¹⁴⁵ the machinery for rapid application to small molecules was already in place. This included bond perception and atom typing algorithms along with the charge assignment protocol compatible with C36.¹⁷,¹⁸ In addition, the CGenFF program outputs penalties associated with charges and parameters not in the existing CGenFF parameter set, where the penalty is assessed based on the similarity between the algorithmically derived parameter and those available in the FF. As the penalty is not a direct measure of the “quality” of a given parameter, in many cases, parameters with relatively high penalties are often appropriate for modeling and simulation. CGenFF has been extended to include sulfonyl- and halogen-containing compounds.⁵⁸⁶,⁵⁸⁷ Treatment of halogens made use of lone pairs on aromatic Cl, Br and I atoms, allowing for modeling of halogen bonds involving weak favorable interactions with hydrogen bond acceptors along the C-X bond. Another major extension included parametrization of non-standard amino acids
mentioned above,⁵²⁶ which were treated with CGenFF combined with the C36 protein FF.

The small-molecule DGenFF was designed to be analogous to CGenFF with some important differences. Notably, DGenFF takes advantage of the original CGenFF program in which the bond perception, atom typing, parameter assignment and bonded penalty assignment algorithms were based on a rule-based approach that creates a rules file specific for the DGenFF while using the same CGenFF program. Assigning electrostatic parameters for the Drude FF requires partial atomic charges, atomic polarizabilities, and Thole scaling factors, where a deep neural network (DNN) was developed for each term,⁵⁸⁸ building upon an earlier DNN model.⁵⁸⁹ Features were based on atom connectivity up to 1–5 bonded atoms along with local through-space atom type–atom type pairs with training targeting QM data on nearly 40,000 small model compounds (<200 Da). This approach rapidly assigns electrostatic parameters along with penalties based on populations of different atom types and their connectivity in the DNN training set. Additional validation against QM dipole moments and molecular polarizabilities on 900 FDA-approved compounds showed excellent agreement, indicating that the method is appropriate for drug-like molecules 200 to 700 Da in size. Note that including distance features in the DNN to model asymmetric electrostatic parameters on lone pairs requires that molecules have approximately correct 3D geometries as well as correct ionization and tautomer states.

Efforts are ongoing to extend the coverage of DGenFF comparable to that of CGenFF. To date, extension to halogens has been completed.⁵⁹⁰ It includes the presence of lone pairs to accurately treat halogen bonds along with careful optimization of the anisotropic atomic polarizability on the halogens Cl, Br, and I and inclusion of LJ parameters on Drude particles, as was done for water–Mg²⁺ interactions discussed above. These latter terms allow for accurate modeling of out-of-plane interactions with hydrogen bond donors where the halogen atom serves as a hydrogen bond acceptor, which is more favorable than halogen bonds and is present in a large number of ligand–protein complexes.⁵⁹¹ Further extension of DGenFF will be facilitated by a DNN-based workflow to optimize LJ parameters of new atom types.⁵⁹² Upon completing full coverage, global optimization of LJ parameters will be undertaken, as done for a subset of atom types in the Drude FF.⁵⁹³ The resulting FF parameters yield good agreement with both pure solvent properties and hydration free energies for a large collection of small molecules. To our knowledge, agreement with both classes of condensed phase properties has not been attained with any additive FF, even when a similar global optimization protocol was used.⁵⁹⁴

A final issue concerning FFs in general is the validity of implementation.⁵⁹⁵ CGenFF is based on a specific algorithm to assign atom types yielding model compound topologies used for parameter assignment and optimization. When other algorithms are used based on analogy, the resulting charges and parameters are inconsistent with those on which the FF was optimized. For correct implementation of CGenFF, individual molecules can be uploaded online by users from educational institutions at https://cgenff.silcsbio.com/ and the CGenFF program can be obtained at no charge for users from educational institutions from Silcsbio LLC (silcsbio.com).

## 11. MIXED QUANTUM MECHANICS/MOLECULAR MECHANICS (QM/MM) METHODS

### 11.1. Background.
QM/MM methods are practical and efficient approaches for simulating chemical reactions in condensed phase including enzyme catalysis.⁵⁹⁶⁻⁵⁹⁹ A QM method is necessary for modeling changes in electronic structure such as bond formation and cleavage, photochemical reactions, and electron transfer in redox catalysis by metalloenzymes. However, it is neither practical nor necessary to treat an entire substrate–enzyme complex and the surrounding solvent at the QM level. This is further complicated by the need to sample multiple protein and solvent configurations to determine the free energy change along a reaction pathway. A combined QM/MM approach addresses this challenge by selectively applying the QM treatment to a region of the system involved in the reaction, such as the substrate, cofactors and key amino acid residues directly participating in the chemical event. This ‘QM subsystem’ is embedded in the rest of the system represented by an MM FF.⁵⁹⁶,⁶⁰⁰,⁶⁰¹ Because of its effectiveness and simplicity, QM/MM methods have become the de facto choice for simulating enzyme reactions.⁵⁹⁷⁻⁵⁹⁹,⁶⁰²,⁶⁰³

A QM/MM approach was first implemented in CHARMM by Field and Bash in 1987,⁶⁰¹,⁶⁰⁴ employing the SE ‘neglect of diatomic differential overlap’ (NDDO) method along with the CHARMM FF. QM/MM approaches have since continuously embraced diverse methods, including ab initio (AI), DFT, and SE QM alternatives. The SE-QM methods encompass both the NDDO-based models¹⁶¹,⁶⁰⁵⁻⁶⁰⁸ and the density functional tight binding (DFTB) methods, also referred to as the self-consistent-charge DFTB (SCC-DFTB) methods.⁶⁰⁹ Both of them are incorporated into CHARMM.

The QUANTUM module was the first SE-QM/MM method implemented in CHARMM, which was based on the MOPAC program (version 4.0).⁶¹⁰ Subsequently, two new NDDO-based SQUANTM³ and MNDO97¹⁴⁷,⁶¹¹ modules were added. The SQUANTM module was based on an implementation in the AMBER program,⁶¹² and the MNDO97 module was derived from a stand-alone MNDO97 program.⁶¹³ The latter has recently been rewritten for computational speed and parallelization.¹⁴⁷ The DFTB method was similarly implemented in CHARMM.¹⁶³ Due to their computational efficiency, these SE-QM/MM methods are frequently used in conjunction with other free energy simulation techniques, including US,¹³⁸ SM,¹³⁹,⁶¹⁴ reaction path,¹⁴⁸,⁶¹⁵ and FEP.⁶⁰⁴

For AI and DFT-based QM/MM methods, CHARMM provides robust interfaces to external softwares including Q-Chem,⁶¹⁶ GAMESS-US,⁶¹⁷,⁶¹⁸ GAMESS-UK,⁶¹⁹ CADPAC,⁶²⁰ and Gaussian16.⁶²¹ In addition, the MSCALE module provides a flexible means for accessing other QM programs, such as MOLPRO.⁶²² Unlike the SE-QM/MM modules of CHARMM, other packages for AI and DFT calculations must be obtained separately. Except for GAMESS-US and GAMESS-UK that can be compiled as a single executable within CHARMM, other packages should be installed separately.

### 11.2. QM/MM Potentials and Practical Considerations.
The effective QM/MM Hamiltonian operator is

$$ \hat{H}_{\text{eff}} = \hat{H}_{0} + \hat{H}_{\text{MM}} + \hat{H}_{\text{QM/MM}} \quad (51) $$

where $\hat{H}_0$ and $\hat{H}_{MM}$ describe the QM and MM subsystems, respectively, and $\hat{H}_{QM/MM}$ describes the interaction between
the two. The latter is further decomposed into electrostatic, van der Waals, and the QM-MM boundary terms:

$$ \hat{H}_{QM/MM} = \hat{H}_{QM/MM}^{\text{elec}} + \hat{H}_{QM/MM}^{\text{vdW}} + \hat{H}_{MM}^{\text{boundary}} \quad (52) $$

In CHARMM, $\hat{H}_{QM/MM}^{\text{elec}}$ is solved self-consistently with $\hat{H}_0$, while $\hat{H}_{QM/MM}^{\text{vdW}}$ is modeled with a LJ potential. The total energy of the system is

$$ E_{\text{tot}} = \langle \Psi_{\text{el}} | \hat{H}^0 + \hat{H}_{QM/MM}^{\text{elec}} | \Psi_{\text{el}} \rangle + E_{QM/MM}^{\text{vdW}} + E_{MM}^{\text{boundary}} + E_{MM} \quad (53) $$

where $\Psi_{\text{el}}$ represents the Hartree-Fock wave function for electrons and nuclei of the QM region. For evaluation, the first three terms in Eq. 53 are determined in the QM/MM module while $E_{MM}$ uses MM energy routines. For performing MD simulation or energy minimization, the QM electron density matrix from the previous MD/minimization step can be used as the initial guess for the next SCF calculation in Eq. 53, with optional addition of small random perturbations to reduce hysteresis and accelerate the convergence of the QM and QM/MM energies. This is implemented in all SE-QM/MM modules as well as in AI-QM/MM modules supporting GAMESS-US, GAMESS-UK and Q-Chem.

In Eq. 53, $E_{MM}^{\text{boundary}}$ addresses cases where the QM and MM division occurs across covalent bonds leaving the QM region with unsaturated dangling bonds. This commonly occurs in enzyme simulations where specific side chains are included in the QM region, for which three methods are available in CHARMM. The first is the hydrogen link (H-link) atom approach, where a hydrogen atom is added within the QM region to saturate and cap the dangling bond.⁶⁰¹,⁶²³,⁶²⁴ The H-link method is conceptually simple, so it has been widely adopted in various packages. In CHARMM, users can introduce an MM angle term to keep the H-link atom aligned with the replaced QM-MM bond. Also, charges on nearby MM atoms can be reassigned to minimize artificial polarization around the QM-MM boundary. The second model is the double-link atom method (DLAM) where an additional H-link atom is introduced at the MM atom site of the QM-MM covalent bond to achieve bond saturation at both of the loose QM and MM ends⁶²⁵ (Section 11.6). This method can be used together with delocalized Gaussian MM (DGMM) charges to mimic the delocalization of charge densities on MM atoms.⁶²⁶ The third is the generalized hybrid orbital (GHO) method.⁶²⁷⁻⁶²⁹ While the first two methods introduce additional degrees of freedom via the link atoms and alter local electrostatic potential by adjusting partial charges, the GHO method treats the QM boundary atom as a special sp³-hybridized carbon and also as an MM atom connected to nearby MM atoms. Its three sp³-hybrid orbitals pointing toward the connected MM atoms called auxiliary orbitals, are fixed with their electron densities assigned based on their MM charges. The remaining sp³-hybrid orbital called the active orbital, is optimized during the SCF iteration. In addition, MM FFs are applied to the GHO atoms to maintain the surrounding geometry.

The above boundary methods, particularly the GHO method, are specifically designed for covalent bonds between two sp³-hybridized carbon atoms, and they are not recommended for arbitrary covalent boundaries. This ensures that the covalent boundary does not perturb the geometry

around the QM-MM bond and minimizes artificial polarization of the QM electron density. The GHO method implemented in CHARMM supports all four SE-QM/MM modules (QUANTUM, SQUANTM, MNDO97 and SCC-DFTB) and AI/DFT QM/MM methods through the GAMESS-US interface. The H-link atom approach is available for all QM/MM modules of CHARMM.

In QM/MM simulations, one must decide on: (1) the QM model, (2) atoms for the QM region, and (3) representation for the QM/MM covalent boundary, and (4) the boundary condition of the whole system.⁶³⁰,⁶³¹ For (1), among AI-QM, DFT and SE-QM, computational errors contributing to the final results should be considered. While AI-QM and DFT methods are more accurate, their high computational cost limits routine use in extensive MD and free energy simulations. Thus, selection of the QM theory level and the basis set should be tailored to individual problems. Also, the high accuracy of computationally expensive QM methods such as a coupled-cluster model may be overshadowed by the statistical noise itself. In such cases, an SE-QM method would be more suitable for lengthy QM/MM simulations. However, they require calibration against AI-QM/DFT levels for the reaction under consideration.⁶⁰⁸,⁶³²⁻⁶³⁴ Among the three NDDO-based SE modules, QUANTUM and MNDO97 in CHARMM now have the option to read non-standard parameters without modifying the source code, obviating the need to rebuild the executable for reaction-specific parametrization.

About the choice of the QM region, there are recent debates about the minimum QM region size required for convergence.³⁷⁷,⁶³⁵⁻⁶⁴⁰ It affects the computational cost and the extent of sampling needed. Advances in efficient algorithms and specialized hardware enable systematic exploration of the QM size for the desired accuracy. However, such an investigation is currently feasible only for relatively small systems, and chemical intuition still remains crucial for selecting the QM region. In any event, we note that QM/MM methods are fundamentally empirical approximations. There is no reason to expect that an arbitrary combination of QM and MM models will reproduce the full quantum results. One should carefully optimize parameters for separating a full QM system into two distinct QM and MM models, to determine the minimum size of the QM region.

CHARMM can perform the QM/MM calculations using both PBC and the solvent boundary condition.⁶³⁰ When using PBC, the QM/MM-Ewald⁶⁴¹ and QM/MM-PME methods⁶¹²,⁶⁴² can be used. In this case, the $E_{MM}$ term in Eq. 53 includes MM-MM interactions with all images as for the regular PME and the $\hat{H}_{QM/MM}^{\text{elec}}$ term includes long-range electrostatic interactions of all MM and QM periodic images with the QM charges. The QM/MM-Ewald method is available in all SE-QM/MM methods, while the QM/MM-PME method is available in the SQUANTM and MNDO97 QM modules; the QM/MM-Ewald method is also supported by the AI/DFT-QM/MM method employing the QChem package.⁶⁴³ For the DFTB method, the GSBP method is also available as an alternative way to incorporate long-range electrostatic interactions into the QM/MM framework.⁶⁴⁴,⁶⁴⁵ Otherwise, it is recommended not to use any cutoff scheme for non-bonded interactions for balanced interactions between the QM-MM and MM-MM pairs.⁶⁴¹ This is because in Eq. 53, any MM atom included in the QM/MM interactions, e.g., those within the cutoff distance of any QM atom, interacts