Minimization with new parameters

In this unit, you'll use our state of the art molecular dynamics program NAMD2 with the CHARMM27 force field and the parameters we developed. You should be in the following directory (if not, go there now):

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.