Difference between revisions of "2008 AMBER tutorial"

Latest revision as of 10:01, 31 January 2011

For additional Rizzo Lab tutorials see AMBER Tutorials.

This tutorial is almost entirely adapted from the tutorial found at the Amber webpage at amber.scripps.edu. The Amber site also contains users manuals available for download.

This wiki tutorial makes modifications to the original Amber tutorial to address the specific software available at SUNY Stony Brook.

The following programs are central to the AMBER package:

Leap - For preparation of the parameter and coordinate files.

tLeap - Terminal version (no GUI)

Antechamber - special program for preparation of parameter sets (“prep” files) for small molecules (drugs/ligands, lipids).

Sander - The number cruncher (aka the “brains” of the Amber MD package) that performs the computations.

Ptraj - Trajectory analysis tool.

Research Problem: We will study the x-ray crystal structure of the oxytocin-neurophysin complex PDB:1NPO. We will focus on oxytocin itself. Remember, the x-ray crystal structure has no hydrogen atoms or lone pairs. The AMBER Leap program will take care of adding hydrogens, etc.

A very important information item in a PDB file of a protein is the SSBOND records, which designate all of the disulfide bonds found in the structure. You will need this information for processing your file in Leap.

Download the protein coordinates

1. Make a project directory and use it for this exercise.

mkdir myproj

cd myproj

2. Download 1NPO.pdb from the protein data bank PDB:1NPO Review the PDB file in a text editor (VI)

vi 1NPO.pdb

You will notice that the structure is a dimer of two complexes. You should also notice that in the REMARK section we find that REMARK 6 informs us that the first several residues and the last several residues (86-95) in neurophysin were not found or resolved in the structure. We will ignore these for now as they are not important. However, if you know that you will need missing residues, you may rebuild them using the SEQRES information in your PDB file and the jackal suite of programs (another tutorial, another day!).

3. Using your text editor, delete all header information from 1NPO.pdb. Also, remove chains A, C, and D along with any water molecules and ions from the PDB file. We only need one complex! Save your new file as oxyt.pdb. You can now use VMD to analyze the structure.

Structure preparation with tleap

4. tLeap - The principal function of this programs is to prepare the AMBER coordinate (inpcrd or crd for short) and topology (prmtop or just top) files.

tleap -s -f leaprc.ff03

The leaprc.ff03 file loads the parameters for the AMBER03 force field. In the tleap prompt type the following commands one at a time:

oxy = loadPdb oxyt.pdb

The bond command creates bonds between atoms. The syntax is: unit_name.residue_number.atom_name. Create the disulfide bond and check the structure using

bond oxy.1.SG oxy.6.SG
check oxy

After using the check command you will notice several warnings about close contacts and possibly a warning about the overall charge of the unit if the unit is not neutral in overall formal charge. The overall charge is important, as we must neutralize this charge for particle-mesh ewald electrostatics later on.

*NOTE*: If an error other than that described above occurs, open your oxyt.pdb file and ensure that all residue types CYS have been changed to CYX then re-run tleap commands.

5. The structure should be fine at this point. You should notice that tLeap has added hydrogen atoms to the peptide structure. Save a topology and coordinate file for in vacuo runs:

saveAmberParm oxy oxy_vac.top oxy_vac.crd

We will use the solvateOct command to solvate the structure in TIP3P water using a truncated octahedron periodic box.

solvateOct oxy TIP3PBOX 9.0

You have told Leap to solvate the unit in a truncated octahedral box using a spacing distance of 9.0 angstroms around the molecule. Ideally, you should set the spacing at no less than 8.5 Å (~ 3 water layers) to avoid periodicity artifacts. For particle-mesh ewald electrostatics, our box side length must be > 2 X NB cutoff. We will use a 10.0 Å cutoff in our solvated system; therefore, our box side must be > 20 Å. Our box side length will be (2 X 9) + peptide dimension, which should easily be greater than 20 Å. The system must be neutral in terms of overall formal charge. Fortunately, our system is neutral as is. If this had not been the case then we would have used the addIons command to neutralize the charge (use Na+ to counter a negative charge or Clto counter a positive charge).

The saveAmberParm command saves the parameter file (prmtop or top file) and the initial coordinate file (inpcrd or crd file). We did this before solvating the system so we could perform an in vacuo dynamics simulation for comparision to the solvated system.

saveAmberParm oxy oxy.top oxy.crd
quit

6. Terminal Leap can be used to automate this process: Create a new file called oxy.scr and paste the following code into it:

source leaprc.ff03
oxy = loadPdb oxyt.pdb
bond oxy.1.SG oxy.6.SG
check oxy
saveAmberParm oxy oxy_vac.top oxy_vac.crd
solvateOct oxy TIP3PBOX 9.0
saveAmberParm oxy oxy.top oxy.crd
Quit

Tleap is very useful in situations where you might need to process a very large structure that requires a large solvent box and ions. Use the following command format for running tLeap.

tleap -s -f oxy.scr

As you can see the multiple steps of the previous portion of the tutorial has been conveniently combined into a single step with one command.

7.You’re done with file preparation (hooray!). You will notice a file called “leap.log” in your working directory. Leap records every command and result (as well as commands issued by Leap under the hood) in a log file. Some of the information in this log file is required for publication (e.g. periodic box size; # of water molecules; # of ions if any, etc).

Make a separate directory for your in vacuo dynamics runs and move your oxy_vac.top and oxy_vac.crd files there.

mkdir in_vacuo
mv oxy_vac.* in_vacuo

Molecular Dynamics with Sander

The in vacuo Model. This computation will be performed in two steps using an NVT ensemble. An energy minimization will be done followed by a production dynamics run.

Step 1. Energy Minimization of the System.

We perform a steepest descents energy minimization to relieve bad steric interactions that would otherwise cause problems with or dynamics runs followed up with conjugate gradient minimization to get closer to an energy minimum. Using the Sander program.

The SANDER program is the number crunching juggernaut of the AMBER software package. SANDER will perform energy minimization, dynamics and NMR refinement calculations. You must specify an input file to tell SANDER what computations you want to perform and how you would like to perform those computations. Study the input file for energy minimization below.

create a new file min_vac.in and paste the following code into the file:

oxytocin: initial minimization prior to MD
&cntrl
imin = 1,
maxcyc = 500,
ncyc = 250,
ntb = 0,
igb = 0,
cut = 12
/

General format of the sander command and its flags (The following is just an example, do not try to run it as it will crash!):

sander -O -i in -o out -p prmtop -c inpcrd -r restrt [-ref refc -x mdcrd -v mdvel -e mden -inf mdinfo]

Issue the following command to run your in vacuo energy minimization.

nohup sander -O -i min_vac.in -o min_vac.out -p oxy_vac.top -c oxy_vac.crd -r oxy_vacmin.rst &

We use “nohup” to prevent hang-ups. The ampersand at the end of the command runs the job in background. Monitor the progress using “tail -20 min_vac.out” command. The in vacuo minimization should proceed very rapidly (It only took 1 second on our 2.2 GHz linux workstation!). It is important to fully understand what the flags in the above command mean and what the following files being called are. In the above example -O means "overwrite" files if they exist, -i specifies the "input" file name, -o specifies the "output" file name, -p specifies the "parameter" file name, -c the "coordinate" file name, and -r the "restart" file name. The Amber manual should be consulted for a more complete description.

Step 2: Molecular Dynamics

Make a new file md_vac.in with the following code:

oxytocin MD in-vacuo, 12 angstrom cut off, 250 ps
&cntrl
imin = 0, ntb = 0,
igb = 0, ntpr = 100, ntwx = 500,
ntt = 3, gamma_ln = 1.0,
tempi = 300.0, temp0 = 300.0,
nstlim = 125000, dt = 0.002,
cut = 12.0
/

Issue the following command:

nohup sander -O -i md_vac.in -o md_vac.out -p oxy_vac.top -c
oxy_vacmin.rst -r oxy_vacmd.rst -x oxy_vacmd.mdcrd -ref
oxy_vacmin.rst -inf mdvac.info &

This computation took about 3.8 minutes on a 2.2 GHz linux workstation. The stuff in the input files in a nutshell. There are many more settings/flags above than we can possibly explain here. So, for more in-depth info please consult the AMBER User’s Manual. The settings important for minimization are highlighted and explained below.

&cntrl and / - Most if not all of your instructions must appear in the “control” block (hence &cntrl).

cut = nonbonded cutoff in angstroms.

ntr = Flag used to perform position restraints (1 = on, 0 = off)

imin = Flag to run energy minimization (if = 1 then perform minimization; if = 0 then perform molecular dynamics).

macyc = maximum # of cycles

ncyc = After ncyc cycles the minimization method will switch from steepest descents to conjugate gradient.

ntmin = Flag for minimization method. (if = 0 then perform full conjugate gradient min with the first 10 cycles being steepest descent and every nonbonded pairlist update; if = 1 for ncyc cycles steepest descent is used then conjugate gradient is switched on [default]; if = 2 then only use steepest descent)

dx0 = The initial step length

dxm = The maximum step allowed

drms = gradient convergence criterion 1.0E-4 kcal/mol•Å is the default

MOLECULAR DYNAMICS IN A WATER BOX

This job will be accomplished in 4 steps:

Step 1. Restrained Minimization - relieve bad vander Waals contacts in the surrounding solvent while keeping the solute (protein) restrained.

Step 2. Unrestrained Minimization - Relieve bad contacts in the entire system.

Step 3. Restrained Dynamics - Relax the solvent layers around the solute while gradually bringing the system temperature from 0 K to 300 K.

Step 4. Production Run - Run the production dynamics at 300 K and 1 bar pressure.

Step 1: Restrained Minimization

Create a file called min1.in containing the following code:

oxytocin: initial minimisation solvent + ions
&cntrl
imin = 1,
maxcyc = 1000,
ncyc = 500,
ntb = 1,
ntr = 1,
cut = 10
/
Hold the protein fixed
500.0
RES 1 9
END
END

Hold the peptide fixed 500.0 (This is the force in kcal/mol used to restrain the atom positions.) RES 1 9 (Tells AMBER to apply this force to residue #’s 1 to 9).

Issue the following command:

nohup sander -O -i min1.in -o min1.out -p oxy.top -c oxy.crd -r
oxy_min1.rst -ref oxy.crd &

Step 2: Unrestrained Minimization

Create a file min2.in with the following code:

oxytocin: initial minimisation whole system
&cntrl
imin = 1,
maxcyc = 2500,
ncyc = 1000,
ntb = 1,
ntr = 0,
cut = 10
/

Issue the following command:

nohup sander -O -i min2.in -o min2.out -p oxy.top -c
oxy_min1.rst -r oxy_min2.rst &

Step 3. Position-Restrained Dynamics

NOTE: IF IN AMS 535 CLASSROOM, CONSIDER LOWERING MAX STEPS TO SAVE TIME

This initial dynamics run is performed to relax the positions of the solvent molecules. In this dynamics run, we will keep the macromolecule atom positions restrained (not fixed, however). In a position-restrained run, we apply a force to the specified atoms to minimize their movement during the dynamics. The solvent we are using in our system, water, has a relaxation time of 10 ps; therefore we need to perform at least > 10 ps of position restrained dynamics to relax the water in our periodic box.

Create a file md1.in with the following code:

oxytocin: 20ps MD with res on protein
&cntrl
imin = 0,
irest = 0,
ntx = 1,
ntb = 1,
cut = 10,
ntr = 1,
ntc = 2,
ntf = 2,
tempi = 0.0,
temp0 = 300.0,
ntt = 3,
gamma_ln = 1.0,
nstlim = 10000, dt = 0.002,
ntpr = 100, ntwx = 500, ntwr = 1000
/
Keep protein fixed with weak restraints
10.0
RES 1 9
END
END

ntb = 1 Constant volume dynamics.

imin = 0 Switch to indicate that we are running a dynamics.

nstlim = # of steps limit.

dt = 0.002 time step in ps (2 fs)

temp0 = 300 reference temp (in degrees K) at which system is to be kept.

tempi = 100 initial temperature (in degrees K)

gamma_ln = 1 collision frequency in ps-1 when ntt = 3 (see Amber 8 manual).

ntt = 3 temperature scaling switch (3 = use langevin dynamics)

tautp = 0.1 Time constant for the heat bath (default = 1.0) smaller constant gives tighter coupling.

vlimit = 20.0 used to avoid occasional instability in dynamics runs (velocity limit; 20.0 is the default). If any velocity component is > vlimit, then the component will be reduced to vlimit.

comp = 44.6 unit of compressibility for the solvent (H2O)

ntc = 2 Flag for the Shake algorithm (1 - No Shake is performed; 2 - bonds to hydrogen are constrained; 3 - all bonds are constrained).

tol = #.##### relative geometric tolerance for coordinate resetting in shake.

You will note that we used a smaller restraint force (10 kcal/mol). For dynamics, one only needs to use 5 to 10 kcal/mol restraint force when ntr = 1 (uses a harmonic potential to restrain coordinates to a reference frame; hence, the need to include reference coordinates with the -ref flag.). Larger restraint forces lead to instability in the shake algorithm with a 2 fs time step. Larger restraint force constants lead to higher frequency vibrations, which in turn lead to the instability. Excess motion away from the reference coordinates is not possible due to the steepness of the harmonic potential. Therefore, large restraint force constants are not necessary.

Issue the following command:

nohup sander -O -i md1.in -o md1.out -p oxy.top -c oxy_min2.rst
-r oxy_md1.rst -x oxy_md1.mdcrd -ref oxy_min2.rst -inf md1.info
&

Monitor the progress using “tail -40 md1.out”. Our run took approximately 32 minutes with no other processes running in background on a 2.2 GHz linux workstation.

Step 4: The Production Run

NOTE IF IN AMS535 CLASSROOM CONSIDER LOWERING THE SIMULATION TIME

Create a file md2.in with the following code:

oxytocin: 250ps MD
&cntrl
imin = 0, irest = 1, ntx = 7,
ntb = 2, pres0 = 1.0, ntp = 1,
taup = 2.0,
cut = 10, ntr = 0,
ntc = 2, ntf = 2,
tempi = 300.0, temp0 = 300.0,
ntt = 3, gamma_ln = 1.0,
nstlim = 125000, dt = 0.002,
ntpr = 100, ntwx = 500, ntwr = 1000
/

ntb = 2 Constant pressure dynamics.

ntp = 1 md with isotropic position scaling.

taup = 2.0 pressure relaxation time in ps

pres0 = 1 reference pressure in bar.

Issue the following command to initiate your molecular dynamics simulation:

nohup sander -O -i md2.in -o md2.out -p oxy.top -c oxy_md1.rst -
r oxy_md2.rst -x oxy_md2.mdcrd -ref oxy_md1.rst -inf md2.info &

Monitor the progress using “tail -40 md2.out”. Our run took approximately 7.5 hours to run on a 2.2 GHz linux workstation.

THAT IS THE END OF THE PRODUCTION TUTORIAL, NOW ON TO ANALYSIS

ANALYSIS