Tutorial: Using custom CHARMM parameters with Desmond in Schrödinger Suite 2013-2

by Dmitry Lupyan (Updated on 10/9/2013)

Introduction

Dictionary

Before you start

Disclaimer

Input files for this tutorial

Custom CHARMM parameter files

Generating CHARMM Ibuprofen parameters with SwissParam

Converting Ibuprofen viparr template

Using the IBP viparr template

Preparing the Protein-Ligand complex for simulation

Solvating and neutralizing the protein-ligand complex

Assigning CHARMM force fields to the molecular system

Adding constraints

Running an MD simulation with Desmond

Introduction

The Desmond MD package includes the latest CHARMM force fields for proteins, DNA and RNA, lipids, and sugars, which can be readily used for molecular dynamics simulations.  Often one would like to run an MD simulation with custom ligand parameters that are not part of the standard CHARMM-family parameter set.  This tutorial covers the basics on how to convert CHARMM parameters to a viparr template.  After the template is prepared, we will use it to assign CHARMM parameters, with viparr, to a molecular system for MD simulation.

Dictionary

• Maestro -- Schrödinger’s 3D visualizer software
• Desmond -- Molecular Dynamics software by D.E. Shaw Research, which is distributed by Schrödinger Inc.  Please visit the DESRES site to learn more
• Viparr Utility  -- Desmond’s tool to manage force fields.  More information can be found in Schrödinger’s Desmond User Manual

Before you start

Make sure you have the $SCHRODINGER variable set to your Schrodinger software installation.

Disclaimer

The converter tool described in this document is not supported by either Schrödinger Inc. or D.E. Shaw Research.  For help please consult other Desmond users via Google Desmond Group.

Input files for this tutorial

⚛All the input files needed to follow this tutorial can be downloaded from Google Drive.

The link contains the following files:

Custom CHARMM parameter files

To simulate a molecule that is not covered by the standard CHARMM parameter set, we first need to obtain a topology and parameter file for that molecule.  The source of the parameters can be one of the following:

• Search the CHARMM message board to see if anyone else is willing to share their parameters | CHARMM forums
• Adopt parameters from the literature
• Use other computational chemistry tools to generate parameters, see the following:

• A. D. MacKerrel’s slide deck | PDF link
• N. Foloppe and A. D. MacKerell (2000). JCC 21 (2) | DOI link

• Use the CHARMM General Force Field (CGenFF) facility | website | DOI link

• If you’re using CGenFF please make sure you read the Note attached to one of the sections below.

• Use the Paratool plugin for VMD | website
• Use the SwissParam server | website | DOI link

Any of the above methods can be used to generate CHARMM-compatible parameters for the viparr converter tool.

Note: to ensure that the parameters you intend to use are the ones that are being assigned, it is best if the topology and parameter files are self-contained and do not depend on other parameter files.

For the sake of simplicity and reproducibility we will use the SwissParam server to generate topology and parameter files for this tutorial.  Please refer to the original SwissParam documentation for more details on parameterization methods.

Generating CHARMM Ibuprofen parameters with SwissParam

Suppose we would like to run a simulation of a protein in complex with an Ibuprofen (IBP) ligand, Figure 1.  The standard CHARMM parameter set lacks parameters for Ibuprofen.  To generate CHARMM parameters for IBP, we first need its structure.  The ligand structure can be obtained or generated by:

• Drawing the ligand using the 2d Sketcher (in Maestro choose Edit → 2D Sketcher)
• Build the molecule with “Fragments” toolbar in Maestro
• Generate the ligand 2D structure from a SMILES string (“[O-]C(=O)[C@@H](C)c1ccc(cc1)CC(C)C”) by pasting it directly into Maestro, and then use LigPrep tool to convert it to 3D.
• Extract it from a PDB file (e.g. 3p6h)

• Use the ibp.pdb included with this tutorial

Note: When converting from 2D to 3D structures, please check you have the correct chirality.

SwissParam accepts only mol2 file format as an input, so after you have IBP 3D structure loaded in Maestro, with proper bond orders and all hydrogens, export it as ibp.mol2.

Upload the ibp.mol2 file to the SwissParam website. The SwissParam server prepares a ibp.zip file, which is ready to download after a short time.  Follow the instructions on the site to download the file.  The zip archive contains two files that we need; please extract them for further use.

• CHARMM topology file: ibp.rdf
• CHARMM parameter file: ibp.par

Converting Ibuprofen viparr template

Now that we have a custom parameter and topology file for IBP, we can convert it to a viparr-formatted template.  To do so, we will run the converter script in the terminal:

The error is encountered because the converter script interprets the atom type ‘CR’ as chromium (chemical symbol ‘Cr’) and atom type ‘CO2M’ is believed to be cobalt (chemical symbol ‘Co’).  To correct this, we need to change the atom types of both CO2M and CR to an atom type that is interpreted correctly by the viparr conversion tool.

Note: A safe name would include a number character in the second position of the atom type.  Also we have to make sure that the new atom type is not already used in the standard CHARMM parameter set.

So make the following change: ‘CR’→‘C2R’ and ‘CO2M’→‘C2OM’ by either editing the ibp.par and ibp.rtf files manually with a text editor, or running the following two SED command in the terminal.

Running the two commands creates two new files: ibp_new.par and ibp_new.rtf with the atom types changed as described.

Note: Please manually inspect the new files to ensure that unintended substitutions were not made.

Now rerun the above command with the updated parameter and topology files as an input.

The converter script now generates the IBP directory which contains the viparr-formatted CHARMM template for Ibuprofen.  During the generation of the original parameter files, the SwissParam server changed the ligand’s residue name to ‘LIG.’  This is OK since viparr uses structure matching to assign force fields parameters, so the residue name is unimportant.

Using the IBP viparr template

In this section we will prepare a protein in complex with IBP for MD simulations with CHARMM force fields. To do this, we will first build the system with Maestro and Desmond tools.  The default force field used by Schrödinger Inc. is OPLS_2005, so we will later reassign CHARMM parameters to the system.

Preparing the Protein-Ligand complex for simulation

To use the IBP viparr template we generated in the previous section, we will use a PDB entry of a protein that is already complexed with IBP.  We will use PDB ID: 3p6h in this section.  To do open Maestro and:

1. Choose  Workflows → Protein Preparation Wizard or click the PrepWiz toolbar button
2. In the Protein Preparation Wizard panel enter ‘3p6h’ in the PDB text box, and click the Import button
3. After the PDB has been imported into the 3D workspace, we want to preprocess the protein and the ligand. In the  Protein Preparation Wizard panel,  make sure the following options are selected:

• Assign bond orders
• Add hydrogens
• Create disulfide bonds
• Cap termini
• Delete waters beyond 5 Å from het groups

4. Click the Preprocess button.
5. A warning window that lists protein preparation problems opens.  You can just click ‘OK’ to close it.

Note: This window just alerts you that there are missing atoms on the terminal residues, which is not problematic since we’re capping these residues anyway.  Also this window alerts you of overlapping atoms and alternative residue positions for some residues in the PDB file.  All of these are just warnings and can be ignored for this system.

Completing these steps produces a Protein-Ligand complex in the Workspace where the partial and formal charges, protonation states, bond orders of protein residues and ligand groups have been determined and assigned.  Hydrogens are also added to the protein, the ligand and the crystallographic water molecules.

Solvating and neutralizing the protein-ligand complex

Now that we have a protein-ligand complex loaded in the Workspace, we can solvate and neutralize the complex.  To do this in Maestro:

1. Open the System Builder by choosing Applications → Desmond → System Builder
2. Click the Run button to start building the system with the default settings.

After a short wait, the solvated system is displayed in the Workspace.  This system is solvated with SPC waters, and neutralized with counterions.  The completed system is also written out to the file desmond_setup-out.cms, which contains the coordinates and the OPLS_2005 parameters of the system.

Assigning CHARMM force fields to the molecular system

The CMS file generated in the previous step is ready to run Desmond, but OPLS_2005 force fields have been used.  In order to use CHARMM parameters, we need to reassign force fields to CHARMM.

The viparr utility comes prepackaged with standard sets of protein parameters.  In an earlier step we generated a viparr template for the IBP ligand with CHARMM parameters (which should be in the IBP directory), so we have all the ingredients to reassign force fields in our system.   To do this, run the following command in a Linux terminal.

This command assigned CHARMM27 parameters to the protein, ions, and waters (tip3p).  We also assigned our custom IBP parameters to the Ibuprofen ligand by specifying the  -d IBP flag.  The output  charmm-out.cms file contains the system’s structural and CHARMM force field information.

Note: If your custom parameters depend on the predefined parameters set, instead of “-d IBP” flag you should use a merge flag “-m IBP”.  This is often the case with CGenFF, where parameters in the prm file are incomplete and depend on other parameters included in the generalized CGenFF set.  In this case your viparr command should also include ‘-f cgenff_base_v2b7’.

Adding constraints

Now that the system force field parameters have been assigned in the charmm-out.cms file,   we need to add hydrogen constraints.  To do this, run the following command:

The resulting file, charmm-start.cms is now ready for use in a Desmond simulation.

Running an MD simulation with Desmond

Now that the system is prepared with CHARMM parameters, Desmond can minimize, equilibrate and run an MD simulation of the system.  In Maestro, do the following:

1. Open the Molecular Dynamics panel from  Applications → Desmond
2. In the Model system section, choose Import from file, and click Browse to locate and load charmm-start.cms
3. Ensure the ‘Relax model system before simulation’ checkbox is selected
4. Enter a Job name and make settings, then click Run to launch the simulation

Note: To test the described method and parameters, we ran a 10 ns simulation, please see the Simulation Interactions Diagram (SID, available in 2013-2 Schrodinger suite release) PDF report here.

For further information please consult the Desmond User Manual.

Acknowledgements

Great thanks to Vincent Zoete and Olivier Michielin, from Swiss Institute of Bioinformatics and coauthors of SwissParam utilities.