Applying RMSD restraints in Sander using the Colvars library
================================================================================

| Timothy J. Giese :sup:`1`, and Darrin M.
  York\ :sup:`1`
| :sup:`1`\ Laboratory 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
-------------------

- `Collective variables library for molecular simulation and analysis programs <https://colvars.github.io/>`__
- `Colvars reference documentation related to RMSD coordinates <https://colvars.github.io/master/colvars-refman-gromacs.html#sec:cvc_rmsd>`__


Tutorial
--------

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

.. contents::
   :local:
   :depth: 4

.. start-tutorial

  

Introduction
------------

The colvars library implements many novel collective variable definitions, including `RMSD restraints <https://colvars.github.io/master/colvars-refman-gromacs.html#sec:cvc_rmsd>`__
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.

.. code-block:: bash

      [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.

.. code-block::

   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.

.. code-block::

   # 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
   }