forcefield-tutorial/forcefield-files
From now on, we will call this directory the working directory.
When you look at the contents of the directory, you should have the following files:
Topology file: top_all27_prot_lipid_prfar_cyg_cyn.inp
Parameter file: par_all27_prot_na_lipids_full.inp
NAMD config file: hisHmini.namd
PDB file: hisH_cyg_mod.pdb
The parameter file was generated with the same protocol used in the
previous section of this tutorial. However, rather than
performing a semi-empirical calculation, the parameters in
par_all27_prot_na_lipids_full.inp
were calculated by full ab-initio quantum mechanics using the 6-31G* basis set, a
process that takes several hours on a fast machine.
First we have to generate NAMD2-compatible PSF and PDB files.
We will do this with PSFGEN. PSFGEN will add hydrogens to the atoms;
in the original pdb file we supplied you, the correct protonation state has
already been specified for each of the residues.
As we did previously, we will be running PSFGEN directly from the terminal. This time we
will run PSFGEN exactly as we had before, only we will use our new topology file. To do this,
open the PSFGEN input file called build_covsub.pgn in a text editor.
You can see that the topology file we specified is the old one, not the new one. We need to update
the PSFGEN input file so that it uses our new topology file.
Change the first line of the PSFGEN input file from:
topology top_all27_prot_na.inp
to:
topology top_all27_prot_lipid_prfar_cyg_cyn.inp
Save the new file with the same name (build_covsub.pgn) in the working directory.
Try running PSFGEN again by typing the following at the prompt:
psfgen < build_covsub.pgn > myoutput2
Again, PSFGEN will take a few seconds to run. When it is finished, it will return the prompt. Look at the output file myoutput2. This time it worked! You can see it has generated a PDB file containing all coordinates including hydrogens (called hisH_covsub.pdb), a PSF file containing all the information about connectivity, mass and charge of each individual atom in the structure (called hisH_covsub.psf), and a logfile of the PSFGEN program's activities in the file called myoutput2.
Now we are going to solvate the system so we can minimize it. To solvate the system, launch VMD from a
terminal window. Be sure you are already in the working directory
before you start VMD. If you open VMD from the start bar or another directory, open a TK Console window
(select Extensions: TK Console from the VMD Main menu) and cd into the directory where your
new pdb and psf files are.
Open the TK Console in VMD (select Extensions: TK Console from the main VMD menu). At the
prompt, type:
package require solvate
solvate hisH_covsub.psf hisH_covsub.pdb -t 9 -o hisH_solv
This will run SOLVATE, a program which creates a box or water molecules around the protein/object of interest. Notice that it created two new PDB and PSF files: hisH_solv.pdb and hisH_solv.psf.
You just created a rather large and well solvated globular protein system. Let's take a look at the newly solvated system you just created. Load the solvated system by selecting File: Load New Molecule. Browse for your new PDB and PSF file and load them both into VMD.
Now you should be looking at your solvated system. Before we minimize the system, we need to be sure that
the system has an overall charge of zero. To do this, in the TK Console window, type:
> set sel [atomselect top "all"]
> eval vecadd [$sel get charge]
The last command should return a number- this number is the overall charge. It should be something very close
to zero. This is good, because it means that our system is already neutral and we will not need to add any ions
to counter the charge. We will run the minimization on the solvated system using the configuration
file hisHmini.namd, which we supply you.
We will do the minimization in the NPT ensemble with a timestep of 1 femtosecond, a uniform dielectric constant of 1 and
periodic boundary conditions. The electrostatic interactions will be calculated with the Particle Mesh Ewald Method.
Go back to a normal terminal window and be sure you are once again in the working directory.
At the unix command prompt in a terminal shell, type:
> namd2 hisHmini.namd > hisHmini.namd.out
This command begins the minimization, which should take about 10 minutes to run. You can use the tail command
to track its progress. To do this, type this at the unix prompt:
> tail hisHmini.namd.out
While it is running, open
the NAMD2 configuration file and briefly explain what
each command does. The commands are defined in the users guide available
online at http://www.ks.uiuc.edu/Research/namd/current/ug/. Look under
the Chapter on Basic Simulation Parameters and the NAMD configuration
parameters. What minimization method is being used? What is temperature
during the minimization run? How many steps are we minimizing for?
The output of your NAMD job can be found in
HisHmini.namd.out. You can use namdplot to look at
energy changes during the minimization. The minimization job you just ran was for
only 100 steps. Minimizing for 10,000 steps
would likely require over 2 hours of compute time for a single processor! We have provided you with the full
minimization run results for a much longer run. All the minimization files are located in the directory
Long_min. To get there from the working directory, at the unix prompt, type:
> cd Long_min
Use namdplot to look at the energy changes over the longer minimization run we provide you:
hisH_solv_min_long.namd.out. To so this, at a terminal window in the Long_min directory, type:
> namdplot TOTAL hisH_solv_min_long.namd.out
We will need to load the DCD file into VMD in order to monitor the movement of the CYG residue over
the course of the minimization run. To do this, in the VMD Main window, click on the hisH_solv molecule, and
then select File: Load Data Into Molecule.
Browse into the Long_min directory and load the hisH_solv_min_long.dcd file into
the hisH_solv.pdb or hisH_solv.psf file you already have loaded in VMD. Once you've done that, you will see that
there are 100 frames.
Highlight the
CYG residue (VMD selected atoms phrase ''resname CYG") and play the trajectory. Watch how the residue moves. You
might also want to highlight the other residues in the catalytic triad (VMD selected atoms phrase ''protein and resid
84 178 180") and watch how the 3 residues move with respect to eachother.
You can determine how much the structure of the covalently bound substrate intermediate has changed during the minimization by using the RMSDscript you learned about in the NAMD tutorial. Compare the first frame and the last frame to see how the novel residue has changed.