Extracting High-Frequency Motions From the NMFF Pathway
  Integrating the low-resolution models with atomically detailed simulations and subsequently determining thermodynamic quantities (e.g., free energy) along the pathway  is of extreme interest in the field of computational biophysics. Here, we present a way to integrate information from low-resolution models with the atomically detailed models. 
  Specifically, we will perform  minimization with restraints starting from each structure along the pathway generated using NMFF. We will use CHARMM for this purpose. The example CHARMM input file (kinase_bead1.inp) that acomplish's this task is as follows:
  
  * This tutorial example illustrates how to perform minimization with restraints 
    * in CHARMM.
    * Note first we will build the missing atoms in the starting structure.
    * In an accompanying example we will do all this with mmtsb command.
    *
  ! Read in the parameter and topology files.
    ! We'll use the $CHARMMDATADIR as our source for top/param files
  open unit 1 read form name "$CHARMMDATA/top_all22_prot_cmap.inp"
    read rtf card unit 1
    close unit 1
  open unit 1 read form name "$CHARMMDATA/par_all22_prot_cmap.inp"
    read param card unit 1
    close unit 1
  ! Read the sequence explicitly 
    !read sequ card unit 5
    !* Sequence for 214 residue ADK
    !* 
    !*
    read sequ card
    * Adenylate Kinase
    * 
  214
  MET ARG ILE ILE LEU LEU GLY ALA PRO GLY ALA GLY LYS 
  GLY THR GLN ALA GLN PHE ILE MET GLU LYS TYR GLY ILE 
  PRO GLN ILE SER THR GLY ASP MET LEU ARG ALA ALA VAL 
  LYS SER GLY SER GLU LEU GLY LYS GLN ALA LYS ASP ILE 
  MET ASP ALA GLY LYS LEU VAL THR ASP GLU LEU VAL ILE 
  ALA LEU VAL LYS GLU ARG ILE ALA GLN GLU ASP CYS ARG 
  ASN GLY PHE LEU LEU ASP GLY PHE PRO ARG THR ILE PRO 
  GLN ALA ASP ALA MET LYS GLU ALA GLY ILE ASN VAL ASP 
  TYR VAL LEU GLU PHE ASP VAL PRO ASP GLU LEU ILE VAL 
  ASP ARG ILE VAL GLY ARG ARG VAL HSD ALA PRO SER GLY 
  ARG VAL TYR HSD VAL LYS PHE ASN PRO PRO LYS VAL GLU 
  GLY LYS ASP ASP VAL THR GLY GLU GLU LEU THR THR ARG 
  LYS ASP ASP GLN GLU GLU THR VAL ARG LYS ARG LEU VAL 
  GLU TYR HSD GLN MET THR ALA PRO LEU ILE GLY TYR TYR 
  SER LYS GLU ALA GLU ALA GLY ASN THR LYS TYR ALA LYS 
  VAL ASP GLY THR LYS PRO VAL ALA GLU VAL ARG ALA ASP 
  LEU GLU LYS ILE LEU GLY 
  generate PRO0 first NTER last CTER setup
  open unit 9 read card name adk_mov1.pdb ! Read coordinates of the first structure 
    read coor pdb unit 9
    close unit 9
  ! Since the structures from the NMFF path do not have hydrogens atoms
    ! we add them here. Use ic build and then hbuild.
    ic param
    ic build
    coor init select hydrogen end
  ! For consistency with defaults in the mmtsb tool minCHARMM.pl, we use
    ! this hbuild set-up
    hbuild atom cdie eps 80.0 cutnb 10.0 ctofnb 7.5 ctonnb 6.5 shift vshift bygr
    update atom rdie eps 1 cutnb 12 ctofnb 11 ctonnb 8 switch vswitch bygr
  ! Now we will do a short minimization with harmonic restraints on CA atoms 
          !Restraint CA atoms with high force constant
  ! Copy current coordinates to comparison set 
      coor copy compare
  cons harm force 100 mass select ( ( resid 1:214 .and. segid PRO0 ) ) .and. ( atom * * ca ) end
  mini sd nprint 10 step 0.02 nstep 500 ! Steepest-Descent minimization
    mini abnr nprint 10 step 0.005 nstep 1000 tolenr 1e-05 ! ABNR minimization scheme
  open unit 10 write form name "movmin1.pdb" ! minimized coordinates
    write coor pdb unit 10
    *
    close unit 10
  STOP
  Use the MMTSB tool convpdb.pl to "convert" from generic resiude naming convention to "charmm22" convention for residues and atoms.
  convpdb.pl -out charm22 mov1.pdb > adk_mov1.pdb
  The ADK conformation is minimized using CHARMM with the command:
      $CHARMMEXEC < kinase_bead1.inp > kinase_bead1.out
  Preferably, use  minCHARMM.pl routine in the MMTSB Tool Set to perform the same task.
  minCHARMM.pl -par minsteps=0,sdsteps=500,sdstepsz=0.02 \
  -par minsteps=1000,minmode=abnr \
  -par trunc=switch,dielec=rdie,cutnb=12,cuton=8,cutoff=11 \
  -par param=22x,cmap\
  -par xtop=top_all22_prot_cmap.inp \
  -par xpar=par_all22_prot_cmap.inp \
  -cons ca self 1:214_100 -cmd mmtsbmin.inp -log mmtsbmin$i.log \
  mov$i.pdb > movmin$i.pdb
  Now minimize other intermediate states along the pathway using the same protocol. You have a choice of either submitting jobs manually after making changes to CHARMM input file (kinase_bead1.inp) to read the next structure along the pathway or execute this simple shell script (runcharmm.sh) to perform the task automatically (it will take approximately 1 hour to minimize all the intermediate conformations). If more than one CPU is available this task can also be performed in parallel. This is particularly attractive feature of the method and makes it computationally tractable to investigate conformational transitions of large biomolecules.
  After minimization jobs are finished, load in VMD all the minimized intermediate conformations along the ADK pathway. You can do this by issuing the following commands in the "tkcon" console of VMD:
  source animatepbs.tcl 
    animatepdbs 1 30 "movmin%0d.pdb"
  Turn on the simultaneously the representation of the molecule  to "Cartoon" for the entire protein and to "Licorice" for only few selected residues, numbered 36, 88, 123, 156, 16 and 13 in the PDB files. Now play the trajectory and   visualize the atomically detailed residue motions that were missing in your initial NMFF pathway. 
  In principle, it is also feasible to extract thermodynamically relevant information from the NMFF pathways. Although, to put this into practice much more computationally intensive calculations are required and is beyond the scope of this workshop. If time permits, interested participants may work on the suggested exercise.