Protocol Sheet: ACETK.SMD.v1
Potential of mean force with Steered Molecular Dynamics
Authors: Ignasi Buch, Nathaniel Stanley, Gianni De Fabritiis (acemd@googlegroups.com)
Copyrights 2010-2011, University Pompeu Fabra. All rights reserved.
Software requirements: ACEMD, Perl and Matlab or Octave
Knowledge requirements: ACETK.EQ, Steered MD
Input requirements: This protocol assumes that you have already built and equilibrated a system. Here, we provide a solvated and equilibrated SH2 complex which you should substitute with yours. If you have not equilibrated your system yet, please refer to protocol ACETK.EQ.
Disclaimer. This protocol is based on a best practice. In no way, we guarantee that it is optimal for your system. You use it at your own risk.
Directory structure:
acetk.smd.v1
|-- analysis
| |-- get_work.m
| `-- parse_force.sh
`-- run
|-- functions.tcl
|-- input
|-- par_all27_prot_na.inp
|-- relax.coor
|-- relax.vel
|-- structure.pdb
|-- structure.psf
`-- structure.restrained.pdb
STEP0: Before You Begin
This section describes how the system was built, and how you could build it on your own. A built system is provided in the attached files if you would like to skip this part.
As with the ACETK.EQ and ACETK.BLDCHARMM protocols, we will use the SH2/peptide complex from the PDB 1LKK for this simulation. However, there are some important differences. You will have to start from the BLDCHARMM protocol again, but with two important changes. First, the axis along which we are going to pull the peptide should be situated along one of the basis axes of the system. This simplifies the procedure. Align the axis of the centers of mass of the protein and peptide along the z-axis (the midpoint between the two centers of mass should be at the origin). Make sure to rotate the entire system (waters, etc.), not just the protein and peptide. If you choose a different axis, you will need to make modifications to the input files accordingly. TCL scripts for VMD for finding the center of mass of a group of atoms can be found online. Reviewing vector and matrix operations in VMD will give you the tools you need to properly align the system along the pulling axis.
Second, the water box must have dimensions that accommodate our separation of the protein and peptide. In water box simulations, molecules that exit one side of the box enter on the other side. If the water box is not large enough, the peptide will be pulled through one side and appear again on the opposite side, potentially interacting again with the protein. To avoid this, we recommend a water box of dimension 64 A in both the x and y axis, and 96 A in the z axis. Ideally the midpoint between the center of mass of the protein and peptide should be roughly in the 1/3 of the way along the z-axis, though it is not critical. Using the solvate plug-in for VMD makes this process rather straightforward. Solvation should be done after the protein/peptide have been aligned properly along the pulling axis.
Once you've built the system, follow the procedure in the ACETK.EQ tutorial to relax and minimize the system. When you're done, continue to the next step.
STEP1: Execute the Steered MD
Run the SMD to sample the reaction coordinate. By default, pulling velocity and distance are 25 A/ns for 25 A.
- 1) Edit 'input' to modify the different variables for the pulling simulations (velocity, Kpull, Krestrain, etc). Default values have been tested and published on this same system. 2) Run the simulation. Save the output from STDOUT in a log file.
acemd input > run.log
STEP2: Compute Work associated with pulling
- 1) Parse the force applied to the ligand from the log file. It will create a file called 'F.dat'.
./parse_force.sh
- 2) Get cumulative work associated with the steering of the ligand along the reaction coordinate. Outputs 'W.dat'.
octave --eval get_work.m