McCallum Group

General Notes

QwikMD creates a narrative file describing what it has done; this is filename.infoMD. The information here is *nearly* good enough to include in a journal article.

QwikMD creates a "setup" folder into which all necessary startup files are copied: parameter and topology files, constructed PDB and PSF, etc.

QwikMD creates a "run" folder which is the run-time input and output folder.

Current SMD issues w/ QwikMD:

1. Error message when an SMD stage is run — some sort of selection issue. This does not prevent the stage from running, but it probably causes the second issue.
2. Run stages will not proceed automatically — Qwikmd does not seem to recognize when a stage has finished. SOLUTION: watch the log file ('tail -f run/Minimization.log', for example). When it is complete, re-load the .qwikmd file, and when the dialog box pops up, include the last stage. Then you should be able to start the next stage. Better solution: Once all the stages are prepared, create a shell script which will run all the configuration files sequentially. Here is simple example (it's safest to specify the FULL PATH when calling files. Use 'pwd' to find the path to your current directory):

#!/bin/csh
cd /Users/mccallum/Res/Liang/mutated_SMD/run
namd3 Minimization.conf > min.log
namd3 /Users/mccallum/Res/Liang/mutated_SMD/run/Annealing.conf > ann.log
namd3 /Users/mccallum/Res/Liang/mutated_SMD/run/Equilibration.conf > eq.log

SMD with QwikMD

You may have to test different values of SMDVel together with 'run nsteps' in order to observe enough pulling in your system. Ideally, SMDVel is small, and you can generate the necessary total 'distance' by increasing nsteps. However, this may require very long run times. For example, with triplex DNA, the provided
SMDVel 5e-6
.
.
.
run 500000

Does not generate significant movement in the Chain C strand.

To start with, a value of
SMDVel 5e-5
coupled with 500000 steps seems to generate about the right amount of movement. Theoretically, using
SMDVel 1e-5
with 2500000 should be equivalent, yet generate less noise in the output data. You should probably decrease the output (DCD and restart, etc) frequency by the same factor.

In order to create multiple SMD.conf files (SMD.0.conf, SMD.1.conf, etc), you may need to simply copy the one that QwikMD generates, and make some changes to the OUTPUT part so you generate unique sets of data.

Here is a sample movie:

ChainC_5e-5

Simple analysis

You are able to plot the force (in pN) versus tilmestep (time) by using the script ft.tcl. First, create an analysis directory (to hold any analyses you might create):
mkdir analysis

To create the force versus time file, simply source the file in your tkcon (NOT terminal):
source ft.tcl

You will be prompted for the force vector values n_x, n_y, and n_z. Because QwikMD aligns the molecule along the force pulling direction, the vector is (0, 0, 1), that is, the three values for n_x, n_y, and n_z are 0, 0, and 1 (this may be found in the configuration file under 'SMDDir 0.0 0.0 1.0'). After these values are entered, the file ft.dat is created in the analysis directory. Here is a plot (from gnuplot). Do you think you can tell when a hydrogen bond has broken from this data?


Further ideas to explore:

• order of residue (hydrogen bond) removal?
• are residues shifting from one location to another as an interim step before full detachment?
• is there a zipper effect?

Advanced analysis — averaging

The SMD force/work results only become meaningful after we combine many simulations and we can get some averaged behavior. This is not as obvious as simply summing up all the individual forces at each step for each simulation.

Advanced Analysis — hydrogen bonding

You may display hydrogen bonds in the triplex by creating a new representation (Style -> Special -> HBonds). To display all of them, enter "all" in the selection entry box. You will probably have to play with color (Color: ColorID) and line thickness (Draw style tab) to see them better. It is desirable to show only the HBonds that involve chain C; however It is difficult to figure out how to make that selection in VMD itself. If you know the three atoms involved (say, through their 'index'), you can simply type out that list such as

index 252 760 254 or index 219 728 220 or index…

(see the Timeline discussion, below)

You can isolate them better in the Hydrogen Bond analysis plugin. This provides a graphical analysis. In the main VMD window, choose Analysis -> Hydrogen Bonds. In the new window that pops up, make the selections

• Selection 1: chain C
• Selection 2: chain A or chain B

(The main VMD selection window does not allow such a bicameral selection, not sure why!) Now at the bottom, hit the "Find hydrogen bonds!" button and a plot will pop up. The data is very noisy, and benefits from multiple runs being combined in some way (and perhaps using a running average in your plotting program).

After that, you may play around with the settings, and choose "Write output to files?" so you may save this data and use a better plotting program.

In this particular case, you should see the number of hydrogen bonds decreasing as chain C is pulled of the duplex. If you go back to the selections, and make Selection 1 "all" (and make Selection 2 "all" or leave it blank), you should see more total hydrogen bonds, with a decrease as chain C is pulled off.

Next we will learn how to drill down and examine specific hydrogen bonds, to see if we can detect any preferential order in how they come apart.

Using Timeline

Timeline has a (feature or bug?!) whereby the calculation becomes slower and slower as the simulation steps are processed. This makes it unusable unless you carefully plan your selection(s) before you choose "Calculate"! Here is how we avoid this slowdown:

1. First we have to choose the hydrogen bonds we are really interested in. We can select these by index (best), or residue id. In the main VMD selection window, make sure the HBond representation selection is "all".
2. Use a CPK representation, and find the hydrogen bonds which involve chain C. Use the atom picker (Mouse -> Query, and then use the VMD terminal; or Mouse -> label, and look in the display) to find these atoms. Write down the three atoms (index) involved (typically an H, N, and O). Remember that a hydrogen bond is defined through N–H...O or the N…H–O.
3. When you have written down all these triplets of atoms, use them as the selections in the Timeline entry (under "molecule" on the left-hand side of the Timeline window). For example "index 252 760 254 or index 219 728 220".

You have to be careful when the Timeline window is open, as it takes over VMD and assumes any mouse action is selecting one of the Timeline selections in the VMD display window.

You can filter by atom types too — "type CA" for example. VMD selections allow wildcards, but they can be nonsensical. For example, nitrogens could be type NN1, NN2U, NN3A, etc. To include all types beginning with "NN.." the correct wildcard is "NN.*" (double quotes INCLUDED. This link is useful to keep handy.

Another way to isolate your selections (index s) is to use the Timeline to highlight triplets in the VMD display window, to make sure these involve chain C. Perhaps the most efficient way to do this is to shorten your trajectory (long enough to make sure all the hydrogen bonds show up somewhere, but short enough to use Timeline efficiently). Because the indexes are sequential chain-by-chain, you may be able to use that to help isolate triplets that involve chain C. For example, (252 760 254) can be seen to involve a single atom on chain C, so chain C indexes are in the 700s. This allows me to find (186 696 188) as a chain C-involved H Bond. (If you have to, you can refer to the PDB file to make sure a particular index belongs to chain C.) Using the previous section's example (Find hydrogen bonds!), you will know about how many total hydrogen bonds you need to identify. In this case, 5 or 6.

Perhaps the simplest way is to look at the HBond Timeline, and see which HBond index go from formed (white) to broken (black) sometime during the trajectory. Then you can go back and put only these indexes in the selection, to isolate them and enable you to save the data (File -> Write Data File).

Here is an example after the selection has been culled down. Note the threshold plot at the bottom. The first slider is set at "0" and the second is set at "1" in order to count the density of Hbonds across these selected triplets.