2. Energy Minimization
After we have created the initial model, it is a good idea to visually inspect the model to make sure everything looks reasonable. You can use the command "ambpdb" to convert the topology and restart file into the pdb format:
| New format PARM file being parsed.
| Version = 1.000 Date = 06/14/03 Time = 14:53:33
$
You can then inspect the pdb file using your favorite visualization software. Again, I use the SwissPDB Viewer here as an example. You should see something like this when you load the pdb file into your viewer.
In the SwissPDB Viewer, there is an option called "aa making clashes" under the "Select" menu. This command will highlight the amino acid residues that are making bad contacts with its neighbors. To visualize the clashes, choose "by selection" from the "Color" menu. You should now see something like this:
The region colored by blue are the residues that are making bad contacts. Not bad, only three of them.
Before we start a molecular dynamics simulation, we need to remove these bad contacts. The reason is that if we start the molecular dynamics with these bad contacts, the energy in that region will be unrealistically high and that can either crash the simulation or cause the trajectory to proceed in an unrealistic direction.
We remove these bad contacts by performing energy minimization. Even if there are no obvious bad contacts, it is still a good idea to run a short energy minimization to relax the structure a bit.
We will perform the energy minimization in 2 stages. In the first stage, we'll only minimize the water molecules and hold the protein and Mg2+ fixed. Following is the control input file I used for this calculation:
&cntrl
imin=1, maxcyc=200,
ntpr=5,
ntr=1,
&end
Group input for restrained atoms
100.0
RES 1 155
END
END
Save these control information in a file called "min.in" feed it to sander like this:
Notice the "-ref" option at the end of the command line. This option is normally not required. In this case, we need to have it because we are doing a restrained minimization and this information is required for sander to be able to map the residue selections we specified in our "min.in". The file "wt1mg.rst" is the same file as the original restart file we generated at the end of tleap preparation. Simply copy the original file wt1mg.crd into wt1mg.rst. (The reason we make another copy of the same file is only for bookkeeping. You can use the same wt1mg.crd if you want.)
This command asks sander to take "min.in" as the control input, wt1mg.parm7 as the parameter file for the system, wt1mg.crd as the input coordinate. It will write out wt1mg.rst as the restart file, and wt1mg_min_water.out as the output file.
After the calculation is finished, check your output to see if it looks something like this:
The calculation may take a couple of minutes to complete, so this is a good time to take a short break.