Methane Tutorial

Overview

• Leap:
  • Build residue
  • Solvate
  • Save 'parm' files for dynamics
  • Save 'parm' files for perturbation with and without bond shrinkage
• Sander:
  • Energy minimize to clean up box
  • Run dynamics to equilibrate
• Interface:
  • Run alternative equilibration with sander
  • Run perturbation forward/back with gibbs

Leap

  Build residue

	% xleap &
	> edit MTH

    Select 'C', draw a carbon, then 'H' and draw the four hydrogens out 
    from the carbon. Use 'Erase' mode if you make a mistake. Don't worry
    about the geometry - that will be cleaned up later.

	


    Now pull down 'Unit / Build' and the geometry is complete. Even simpler
    would have been to just draw the carbon and choose 'Add H & Build', but
    it's fun to draw the hydrogens.

    Hold down the two right buttons and push forward/back to zoom in/out.
    The middle button alone rotates, the right button translates, and
    the spacebar recenters the molecule.

    To complete the model, we must add atom types and charges.  While we are
    at it, we'll add the perturbed parameters for 'disappearing' the molecule
    as well. Choose 'Select' mode and snap a box around the molecule. Then pull 
    down 'Edit / Edit selected atoms' and fill in the table. Clicking on a cell
    'activates' it, then you can type in; you can also paste into a cell (using
    the middle mouse button) without activating it. When done, pull down
    'Operations / Check Table'.

	

    The 'PERT.name' is not important; the PERT.type will be used for selecting
    the end-state VDW, bond, and angle parameters.  Now pull down 'Table / Save
    and quit', and deselect the atoms by holding down the SHift key while
    clicking the left mouse button anywhere in the blackness away from the
    atoms (this is more for cosmetic purposes). Now pull down 'Unit / Close'
    and save this model for future reference from the main window:

	> saveoff MTH methane.lib

    In a later session, this molecule can be loaded by

	> loadoff methane.lib


  Solvate:

	> x = copy MTH
	> edit x
	> solvatebox x WATBOX216 10

	


  Save 'parm' files for dynamics:

	> saveamberparm x mebox.top mebox.crd


  Save 'parm' files for perturbation:

    With no bond shrinkage:

	> fmod1 = loadamberparams me_1.frcmod
	> saveamberparmpert x me_p1.top xxx

    With bond shrinkage:

	> fmod2 = loadamberparams me_2.frcmod
	> saveamberparmpert x me_p2.top xxx

    (The coordinate file has the throwaway name xxx because the one
    resulting from the dynamics equilibration will be used for each
    perturbation.)  Since this is a simple system, we can immediately
    spot the relevant difference in the perturbed bond length:

	% diff me_p1.top me_p2.top
	781c781
	<   1.09000000E+00  1.09000000E+00  9.57200000E-01  1.51360000E+00
	---
	>   1.09000000E+00  2.00000000E-01  9.57200000E-01  1.51360000E+00

Sander

  Energy minimize:

     This is done to clean up any steric problems at the box boundary. 
     These can occur because the initial, replicated periodic box (WATBOX216) 
     is trimmed at the cutoff we specified (10A) without considering the 
     waters on the opposite sides of the box:

	% sander_classic -O \
		-i min.in \
		-p mebox.top \
		-c mebox.crd \
		-o min.out \
		-r min.rst


  Equilibrate the system: note

	% sander_classic -O \
		-i md0.in \
		-p mebox.top \
		-c min.rst \
		-o md0.out \
		-e md0.en \
		-x md0.crd \
		-r md0.rst

  The protocol is very simple because this is a conformationally
  limited system, i.e. the methane can't adopt any 'wrong'
  conformation. It turns out that the final density could be
  improved (the initial one is low owing to the superimpose-and-
  subtract algorithm used to place waters), so we run some more 
  dynamics, adjusting the .in to use the velocities and box from 
  the previous run (IREST=1, NTX=7):

	% sander_classic -O \
		-i md1.in \
		-p mebox.top \
		-c md0.rst \
		-o md1.out \
		-e md1.en \
		-x md1.crd \
		-r md1.rst

  Finally, another run with the pressure coupling relaxed. This
  could probably have been done for the previous run (or possibly
  the equilibration could have been considered complete at md0 -
  see note).

	% sander_classic -O \
		-i md2.in \
		-p mebox.top \
		-c md1.rst \
		-o md2.out \
		-e md2.en \
		-x md2.crd \
		-r md2.rst

Interface

  An alternative equilibration:

    Simple minimization
    Dynamics with 5 temperature schemes

  Perturbation

    Simple Gibbs perturbation
    Gibbs perturbation with electrostatic decoupling

  Equilibration + perturbation together

    Run everything: min, dynamics, 
	  perturbation 1->0, dynamics, perturbation 0->1