Applying RMSD restraints in Sander using the Colvars library

Timothy J. Giese 1, and Darrin M. York1
1Laboratory for Biomolecular Simulation Research, Institute for Quantitative Biomedicine and Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, NJ 08854, USA

Learning objectives

  • Use colvars to apply coordinate root mean square deviation (RMSD) restraints within Sander simulations

Relevant literature

Tutorial

In this Tutorial, you will learn how to activate and run restrained dynamics in Sander using the colvars library.

Introduction

The colvars library implements many novel collective variable definitions, including RMSD restraints The colvars implementation of RMSD restraints allows one to rotate and translate a reference structure onto the current configuration using a selection of atoms, and then calculate the RMSD value from an independent selection of atoms.

Use colvars restraints within Sander

The official version of AmberTools does not have the colvars library. It is only available within internally modified York Group versions of AmberTools. There are 2 versions available on Amarel.

[user@amarel] module use /projects/community/modulefiles
[user@amarel] module use /projectsp/f_lbsr/YorkGroup/software/modulefiles
[user@amarel] module use /projectsp/f_lbsr/YorkGroup/software/23jun25/modulefiles
#
# Now load the version of AmberTools with colvars
#
[user@amarel] module load 23jun25/ambertools/libcolvars
#
# It is also available in the HFDF version of AmberTools, but
# this should only be used with ab initio qm/mm simulations
#
[user@amarel] module load 23jun25/ambertools/hfdf

The standard RMSD restraint functionality is specified in the mdin file by setting ntr=1 and restraint_wt=100 in the &cntrl section, and providing a reference structure on the command line with -refc. In constrast, the colvars implementation requires one to set the following options in the mdin file, and the -refc option is not used.

Title: Umbrella sampling example
&cntrl
   [...simulation options...]
   !
   ! Do not activate NMR restraints
   !
   nmropt = 0
/

&colvars
   icolvars=1
   output="restraints.cv.dumpave"
   configfile="restraints.cv.disang"
/
  • icolvars: integer, default=0; Activates colvar restraints if set to 1

  • configfile: string, default=””; File that defines the variables and restraints

  • output: string, default=””; The output file containing the observed collective variable values.

Defining the RMSD restraint within the colvars CONFIG file

An example colvars config file is shown below. It defines a RMSD coordinate and a biasing potential. The sander/colvars interface has been modified to read reference coordinates in the amber-native restart file format. The filename is specified using the refPositionsFile The biasing potential is an effective half-harmonic.

# The colvarsTrajFrequency option controls how often
# the collective variables are written to the dumpave
# file

colvarsTrajFrequency 100

# This section defines the RMSD coordinate

colvar {

   name rmsd
   width 1.0

   rmsd {

      # The reference coordinates

      refPositionsFile ref.rst7

      atoms {
         name ligand

         # The atoms used to calculate the RMSD value once
         # the coordinates have been fit. In this example,
         # we are pretending that the first 33 atoms are a
         # protein-bound ligand

         atomNumbersRange 1-33

         centerToReference

         rotateToReference

         fittingGroup {
            name protein

            # The atoms used to calculate the translation and rotation
            # In this example, we are pretending the selected atoms
            # are part of the protein near the binding pocket.

            atomNumbers { 838 840 842 845 848 851 853 854 857 860 861 }
         }

         # The reference coordinates (again)

         refPositionsFile ref.rst7
      }
   }
}

# This creates a biasing potential that uses the
# rmsd coordinate.

harmonicWalls {

   # This effectively creates a half-harmonic that is
   # flat from rmsd = [0.0,0.1]

   # The force constants are for a potential 1/2 k dx^2
   # In other words, they appear to be twice as large
   # as what one would type into the mdin file or
   # NMR restraint file.

   name bias_rmsd
   colvars rmsd
   lowerWalls 0.0
   upperWalls 0.1
   lowerWallConstant   0.00000001
   upperWallConstant 200.00000000
}