ParmEd

ParmEd (PARMtop EDitor) is a program I wrote for modifying topology files. This page is dedicated to describing why the program exists in the first place, and how to use it to do anything you might want to do.

There are two versions of ParmEd: (1) parmed.py which provides a text-based environment that operates very similarly to ptraj and cpptraj and (2) xparmed.py which is a Tkinter-based GUI version to give a (hopefully) intuitive and point-and-click interface to ParmEd functionality.

Why ParmEd?

Topology files are very complex objects, and there's very little 'locality' in the prmtop. That is, adding a bond is not as easy as just modifying the BOND sections — you need to adjust the appropriate POINTER variables, the exclusion list, etc. Amber did not have any programs that allowed you to directly modify topology files. For instance, if you wanted to change the default implicit solvent radii (set default PBradii in LEaP), you actually had to create a whole new topology file using LEaP.

ParmEd addresses this — it modifies topology files directly, and exposes the data so that you can interact with prmtop files while it takes care of many of the details (like updating pointers, etc.). The design principle of ParmEd was to get out of your way as much as possible — you can modify the topology files in ways that are impossible otherwise (except by hand, and I do not suggest that).

While this program gives you a ton of flexibility with topology files, that flexibility comes with a cost. That is, ParmEd will not prevent you from creating a completely ridiculous and stupid topology file that actually produces (meaningless) numbers. It is important to know what ParmEd is actually doing so you don't assume it's doing more than it is.

Using ParmEd

Using parmed.py

This is the graphical, cpptraj-like interface to the ParmEd functionality. Like cpptraj, parmed.py can be run using the syntax:

parmed.py [<prmtop>] [<script>]

If you do not provide a prmtop, you will be dropped into the interpreter with no loaded topology file. If you do not supply a script with ParmEd commands, it will read your commands from standard input, just like cpptraj.

If you want help at any point, just use the help [<action>] command. If you do not provide an action, you get a detailed overview of all available commands. If you do provide an action, you will get a detailed overview of how to use that particular action and what it does. Shown below is an example:

bash$ parmed.py -n

Reading input from STDIN...
> help

Documented commands (type help <topic>):
========================================
EOF                 checkValidity     parm            quit          writeFrcmod
addAtomicNumber     combineMolecules  parmout         scee          writeOFF   
addDihedral         defineSolvent     printAngles     scnb        
addExclusions       deleteDihedral    printBonds      setAngle    
addLJType           go                printDetails    setBond     
change              help              printDihedrals  setMolecules
changeLJ14Pair      listParms         printFlags      setOverwrite
changeLJPair        loadRestrt        printInfo       shell       
changeLJSingleType  ls                printLJMatrix   source      
changeProtState     netCharge         printLJTypes    strip       
changeRadii         outparm           printPointers   tiMerge     

> help addLJType
addLJType <mask> [radius <new_radius>] [epsilon <new_epsilon>] [radius_14 <new_radius14>] [epsilon_14 <new_epsilon14>] [parm <idx>|<name>]

   Turns given mask into a new LJ atom type. It uses the radius and Rmin from
   the first atom type in <mask> if new_radius or new_epsilon aren't provided

>

You can output a topology file at any time using outparm prmtop_name. There is also a parmout prmtop_name command that will execute after every command is read in or when you type go (it will not execute if you type quit). My suggestion is to just use outparm to print out prmtops when you intend to. The AmberTools manual has a detailed section about all of the available actions and what they do (and how to use them), so that is a good reference as well.

All command names are case-insensitive (but all arguments provided to these commands are case sensitive)

Using xparmed.py

The GUI version of ParmEd is a simple GUI that just provides every action as a button. To run it, use the command:

xparmed.py [<prmtop>]

Note that xparmed.py will not take a script. If you do not provide a prmtop, you will get a file selection box to pick one. If you do provide a prmtop, it will automatically be loaded.

After you've loaded a prmtop, you'll get a window that looks like the following:

xparmed_window.png

A couple notes about using xparmed.py:

  1. There is no setOverwrite command, since the file dialogs automatically ask you if you want to overwrite a file.
  2. It cannot be scripted — for scripting purposes use parmed.py

Adjusting van der Waals interactions

One of the places that ParmEd gives you flexibility that doesn't exist elsewhere is in VDW parameter manipulation. Thomas Steinbrecher has provided a very helpful document about the VDW equation and how Amber and LEaP implement it. It is available at http://ambermd.org/vdwequation.pdf.

How VDW equation is implemented in Amber

This document explains how each atom has a radius and well depth, and how they are combined with different atom pairs to give combined well depths and radii. Specifically, between pairs i and j, the well depth $\epsilon _ {i,j}$ is the geometric average and the radius $r_{i,j}$ is the sum

(1)
\begin{align} \epsilon_{i,j} = \sqrt {\epsilon _ i \epsilon _ j} \\ \\ r_{i,j} = r_i + r_j \end{align}

These combined radii and depths are then combined into A coefficients and B coefficients using the equations:

(2)
\begin{align} ACOEF _ {i,j} = \epsilon_{i,j} r_{i,j} ^ {12} \\ \\ BCOEF _ {i,j} = 2 \epsilon_{i,j} r_{i,j} ^ {6} \end{align}

Equations 1 and 2 are evaluated in LEaP, and $\epsilon _ i$ and $r_i$ are provided as input in the parameter files. Because there are more ACOEF and BCOEF parameters than there are input parameters, the way LEaP handles VDW parameters restricts some flexibility. You can think of the A and B coefficients as a matrix of pairwise combined terms (see Eq. 2) in which you just get to specify the diagonal terms. You cannot independently set how each type will interact with each other type independently. For instance, suppose you want a bunch of Oxygen molecules, but for some reason you don't them to interact with each other (but you want them to interact with everybody else) via VDW terms. This is impossible to do with LEaP, since you cannot selectively zero-out an off-diagonal A or B coefficient.

We switch gears a little here and talk about how LEaP compresses the parameter count in topology files. Because the number of ACOEF and BCOEF parameters scale as $N_{types} ^ 2$, LEaP compresses them as much as possible by assigning all atoms that have the same starting $r_i$ and $\epsilon _ i$ to the same Lennard Jones (LJ), or nonbonded type. You cannot assume that if two atoms have different types (e.g. OH for hydroxyl oxygen vs. O for carbonyl oxygen), that they will also be different types for the VDW sections.

ParmEd provides a command, printLJTypes that will print out every atom that shares the same LJ type. This tells you which atoms will be affected if you change how a particular LJ atom type behaves. You should always use this command to make sure you're not changing atoms you don't want to.

If it turns out that a particular LJ type applies to more atoms than you want to modify, you can use the command addLJType to turn the selected atoms into a new LJ type, so you can modify exactly what you intend to modify.

The next sections will cover what ParmEd will allow you to do, along with some of the caveats.

Changing Lennard-Jones Parameters

ParmEd will let you change the default $r_i$ and $\epsilon _ i$ for any Lennard Jones (LJ) type. Do this using the changeLJSingleType command in parmed.py or xparmed.py

Here, you select an atom and provide a new input $r_i$ and $\epsilon _ i$ and ParmEd will apply the standard LEaP combining rules described in Eq.(1) with all of the other types to regenerate a new A and B coefficient arrays. Note, this will overwrite any off-diagonal terms that you changed before you ran this command.

To make sure you only adjust what you want to adjust, make sure to use the printLJTypes and addLJType commands described above as appropriate.

Changing off-diagonal Lennard-Jones Parameters

NOTE: pmemd.cuda from Amber 12 and older does not use the A/B Coefficient arrays as stored in the topology file. It recomputes the radius and well depth and recombines them. Therefore, modifications made in this section will not work for the GPU implementation of pmemd unless you upgrade to Amber 14 or later.

ParmEd will let you change off-diagonal Lennard Jones (LJ) interactions. That is, you can change how atom type 1 interacts with atom type 2 without changing how atom type 1 interacts with any other atom type. Do this with the changeLJPair function. This action takes pre-combined $r_{i,j}$ and $\epsilon _{i,j}$ values, so you have to combine them yourself (using Eq. 1 if you want).

To make sure you only adjust what you want to adjust, make sure to use the printLJTypes and addLJType commands described above as appropriate.

This is similar to the FIXNB functionality in CHARMM (I've never used CHARMM myself, but I heard about this feature).

CHAMBER note

The CHARMM force field does not scale 1-4 non-bonded VDW terms. Instead, it has a different set of A and B coefficients for these interactions. To also adjust how two atoms will interact in a 1-4 interaction, use the changeLJ14Pair command.

Example

This example will demonstrate you can pick two atoms in a given topology file and eliminate their VDW interactions without affecting how they interact with all other atoms. As an example, let's use atom numbers 9 and 35.

To do this, follow these steps:

  1. Make atom 9 a new LJ type
  2. Make atom 35 a new LJ type
  3. Change how atom 9 and atom 35 interact

To implement this, use the following parmed.py script:

addLJType @9
addLJType @35
changeLJPair @9 @35 0.0 0.0

This will eliminate the VDW interaction between (only) atoms 9 and 35. Because I did not specify a well depth or radius for either of the addLJType commands, it uses the original depth and radius from the original atom type.

Setting up for Hamiltonian Replica Exchange Simulations

ParmEd is also helpful for setting up for H-REMD simulations in many instances. For instance, if you're looking at a deprotonation event, or disappearing a particular atom, you can create a series of topology files in which the partial charges and VDW parameters for a particular atom (or set of atoms) change according to some alchemical ($\lambda$) parameter.

To change the charges, use the command

change charge <mask> <charge>

This will change the charge of every atom in mask to charge

You can also slowly disappear the VDW radius if desired by creating a new type with a shallower well and shorter distance using addLJType <radius> <epsilon>.

Changing Parameters

ParmEd gives users the flexibility to directly manipulate internal potential parameters (e.g., bond, angle, and dihedral parameters).

NOTE: You need ParmEd from AmberTools 14 or higher (or the standalone ParmEd from Github version 1.0 or higher) in order to have these commands work for chamber-style topology files.

Changing Bonds

Bonds can be created or changed using setBond. If the bond does not exist, it will be created with the given force constant and equilibrium distance. If it does exist, the force constant and equilibrium distance will be changed. Only the bond specified will be changed — it will not change the bond type of bonds that share the same type. It will create a new bond type if necessary for that bond.

Note: If you add a bond using setBond, it will NOT add any angle or dihedral parameters for the angles/dihedrals that may have been created by your bond. You will have to add those yourself if you want them. The exclusion list is automatically updated so these atoms will be excluded from interacting with each other in non-bonded interactions.

Changing Angles

Angles can be created or changed using setAngle. If the angle does not exist, it will be created with the given force constant and equilibrium angle. If it does exist, the force constant and equilibrium angle will be changed. Only the angle specified will be changed — it will not change the angle type of angles that share the same type. It will create a new angle type if necessary for that angle.

Note: If you add an angle using setAngle, it will NOT add any bond or dihedral parameters for the bonds/dihedrals that may have been created by your angle. You will have to add those yourself if you want them. The exclusion list is automatically updated so these atoms will be excluded from interacting with each other in non-bonded interactions.

Changing Dihedrals

Changing dihedrals is more difficult than changing bonds and angles, since dihedrals can be composed of multiple Fourier terms. As a result, "changing" an existing dihedral is ambiguous (do you change one term? all of them? get rid of all and replace it with your new term(s)?)

As a result, there is no setDihedral like there is a setAngle and setBond. There are 2 other commands instead, which you will have to combine in order to effect a change in a given dihedral parameter:

addDihedral

This adds a dihedral. You can add an improper dihedral, or any number of terms to a multiterm dihedral.

When adding a multiterm dihedral, make sure the last term in that dihedral is a normal type (the rest must be multiterm). This affects how the 1-4 nonbonded interactions are calculated. They are calculated only for normal-type dihedrals — since you don't want to calculate them multiple times for a dihedral that has multiple terms (so it's only calculated for the 'last' term). Note, if you are adding a dihedral to a ring system, check to see if the 1-4 interactions should be calculated. If not, make this a multiterm dihedral. If you are unsure, don't use this function — let LEaP do it.

Note, you also have to supply a 1-4 EEL scaling factor (SCEE) and 1-4 VDW scaling factor (SCNB). This is the value by which the EEL and VDW interactions between the 2 ends of a dihedral are divided. By default, SCEE=1.2 for Amber force fields (1.0 for GLYCAM) and SCNB=2.0 for Amber force fields (1.0 for GLYCAM).

deleteDihedral

This will remove all dihedral terms between the specified atoms.

Changing Implicit Solvent Radii

Finally! If you're doing MM/PBSA analysis or something, or you just want to switch GB models, and you need a new set of implicit solvent radii for it, you don't have to go through LEaP anymore! You can change it directly using changeRadii <radius_set> where radius_set is any of the radius sets you can set in LEaP (e.g., mbondi, mbondi2, amber6, mbondi3, and bondi)

Generating Stripped Topologies

You can strip a set of atoms from your topology files as well. It will print out a new topology file with all of the remaining parameters (although they will be reordered, so "diffing" them won't tell you anything helpful).

You can strip atoms that are involved in bonds with atoms that aren't stripped — those bonds will simply be deleted (same with Angles and Dihedrals). Use the strip command to access this functionality.

ParmEd Extensions and Scripts

I have tried to make ParmEd easily extendable and more useful in three ways:

  1. I have exposed a ParmEd API in Python so you can write your own Python scripts and use the ParmEd machinery that is already present.
  2. parmed.py has a built-in Python interpreter (limited for security) so you can perform tasks that ParmEd does not explicitly support.
  3. There is a source command that will load a ParmEd script and execute every command inside it. This is useful for creating "actions" that can be performed as a combination of existing actions.

List of available potentially useful ParmEd scripts

This is a section with ParmEd scripts (typically used via the source command) that are not full actions, but perform potentially useful tasks. Feel free to send me a message on Wikidot if you have any scripts you would like added (with appropriate attribution) to this list.

  • water_to_tip3p.parmed: This ParmEd script will convert all :WAT residues with extra points (e.g., TIP4P, TIP5P, etc.) into TIP3P waters. It works on the currently active topology file followed by sourcing it (source water_to_tip3p.parmed). You should load a restart file with the topology file so you can write out a coordinate file to accompany the new topology file.
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License