pyCHARMM: Embedding CHARMM Functionality in a Python Framework

Abstract

CHARMM is rich in methodology and functionality as one of the first programs addressing problems of molecular dynamics and modeling of biological macromolecules and their partners, e.g., small molecule ligands. When combined with the highly developed CHARMM parameters for proteins, nucleic acids, small molecules, lipids, sugars, and other biologically relevant building blocks, and the versatile CHARMM scripting language, CHARMM has been a trendsetting platform for modeling studies of biological macromolecules. To further enhance the utility of accessing and using CHARMM functionality in increasingly complex workflows associated with modeling biological systems, we introduce pyCHARMM, Python bindings, functions, and modules to complement and extend the extensive set of modeling tools and methods already available in CHARMM. These include access to CHARMM function-generated variables associated with the system (psf), coordinates, velocities and forces, atom selection variables, and force field related parameters. The ability to augment CHARMM forces and energies with energy terms or methods derived from machine learning or other sources, written in Python, CUDA, or OpenCL and expressed as Python callable routines is introduced together with analogous functions callable during dynamics calculations. Integration of Python-based graphical engines for visualization of simulation models and results is also accessible. Loosely coupled parallelism is available for workflows such as free energy calculations, using MBAR/TI approaches or high-throughput multisite λ-dynamics (MSλD) free energy methods, string path optimization calculations, replica exchange, and molecular docking with a new Python-based CDOCKER module. CHARMM accelerated platform kernels through the CHARMM/OpenMM API, CHARMM/DOMDEC, and CHARMM/BLaDE API are also readily integrated into this Python framework. We anticipate that pyCHARMM will be a robust platform for the development of comprehensive and complex workflows utilizing Python and its extensive functionality as well as an optimal platform for users to learn molecular modeling methods and practices within a Python-friendly environment such as Jupyter Notebooks.

Figure 1:

pyCHARMM architecture. (Top) The end user writes a Python script to execute CHARMM commands by their Python expressions after importing the necessary Python modules in pyCHARMM. (Middle) The functions available to the user are the Python modules in the Python API, which in many cases represent bindings to the Fortran routines in the pyCHARMM API (or the Fortran module). The pyCHARMM API is a collection of Fortran modules that is included with the CHARMM source code.

Figure 2:

Building an alanine dipeptide. An illustration of the correspondence between CHARMM scripting language commands and their expression as Python functions. Shown are the specific tasks of building the structure of Ala-dipeptide, assigning coordinates to it, and setting up a list of non-bonded interaction parameters for energy minimization calculations.

Figure 3:

Configuring and performing a molecular dynamics run using the ‘DynamicsScript’ class in pyCHARMM. The correspondence with the DYNAmic command in CHARMM in terms of the use of command keywords as function arguments to instantiate an object of the class is shown and the presence of different categories of run parameters is illustrated. In Figure S1, the use of Python dictionaries for the same purpose is shown.

Figure 4:

(a) Illustration of the workflow in the absolute solvation free energy calculation pipeline. Based on the thermodynamic cycle, the absolute solvation free energy is obtained from the difference between ΔG_vaccum and ΔG_water, i.e., the free energy change of molecule annihilation in vacuum and water, respectively. To compute ΔG_vaccum and ΔG_water, multiple molecular dynamics simulations are performed at different λ values and potential energies are computed for each trajectory under perturbed thermodynamic states. These simulations and energy calculations can be executed in parallel to reduce the overall wall-time. Finally, FastMBAR is utilized to calculate the free energy difference between different alchemical states. (b) The absolute solvation free energies of 206 representative molecules from the FreeSolv database obtained using a protocol written in pyCHARMM based on the examples noted above is compared with their experimental values. The dashed line denotes the y=x line. Representative molecules were obtained by sampling from 50 functional fingerprint clusters into which all 642 molecules of the FreeSolv database was partitioned. Individual molecules were parameterized with CGenFF force field and program.^, A representative example of this workflow is included in the GitHub repository as noted above. A detailed spreadsheet containing the molecules SMILES strings, the computed free energies and the experimental free energies taken from the FreeSolv database are included in the SI.

Figure 5:

Illustration of pyCHARMM CDOCKER. These complex workflows have been compiled into a single Python module (pycharmm.cdocker) and standard docking calculation can be performed through one-line of Python code. In this example, a ligand has been docked at a binding site after CDOCKER searched a cubic volume of edge length of 25.762 Å, centered at (12.33 Å, 33.48 Å, 19.70 Å). See the GitHub site pyCHARMM-Workshop for complete examples of rigid receptor - flexible ligand CDOCKER and flexible CDOCKER.

Figure 6:

Demonstration of the string method applied to determine the minimum energy path between two metastable states of an alanine dipeptide in vacuum using the harmonic Fourier beads method. The initial path is shown in red and the path after 2000 iterations is shown in black. See the GitHub site pyCHARMM-Workshop for a complete example.