Dear all,
I would like to calculate pairwise non-bonded interaction energies between amino acid residues in a protein over the length of MD simulation trajectories produced by the NAMD 2.9 software using the CHARMM 27 force-field.
Currently I use the namdenergy.tcl script for this purpose. I apply a minimum distance cutoff between alpha carbon positions to reduce the number of combinations as well as a stride for frames to be read. I have also written a Python program which divides the list of combinations into chunks (to take advantage of the multiple processors) and calls vmd (and subsequently the namdenergy) using an another external tcl script in parallel processes via the python multiprocessing module.
I want to include as many frames as possible in my calculation to avoid missing possibly interesting changes in the interaction energies since I also calculate the correlations between them in a later analysis step. However, when I reduce the stride value, the time required to finish the computation increases rather exponentially.
Namdenergy creates a temporary DCD file prior to each calculation. Also, if I move the files to a solid-state drive and run the calculation there, it speeds up by a factor of two, leading me to think that the computation is I/O bound.
So is there a more efficient/faster way to do this kind of computation (e.g. another computational tool) or is this just a speed limit that I have to live with?
In my experience, if one already has a DCD trajectory file from a previous NAMD calculation and in VMD (1) loads the PSF file, (2) loads onto it the DCD file, and (3) selects the "Extensions, Analysis, NAMDEnergy" tool, the job you comment is done very fast.
I've used that tool with several chain/residue selections in proteins and protein-ligand complexes of disparate sizes and with different selections of energy terms for loaded DCD files having about 2000-2500 frames and data (which you can directly plot or save in a separate text file with the tool) were rapidly obtained in a PC running Windows 7.
What are the numbers involved (protein size, number of frames, ...) in your case?
Thank you Dr. Pacios for your answer. In my case, the protein complex has 385 residues. I have a total of 50 MD simulations. Each simulation is 50 ns long with 50000 frames (I remove solvent molecules prior to analysis). With a stride value of 100, I analyze every 100th frame from the simulation which corresponds to 500 frames per simulation.
If I were interested in only the interactions of all residues with e.g. a pre-defined binding pocket, this speed limitation wouldn't be much of a problem. However, I am interested in investigating the changes in the whole complex, therefore I need to calculate the interaction energies for all possible logical combinations of residue pairs.
For 385 residues, the number of all unique pair combinations is about 74111. It is true that namdenergy.tcl is quite fast, but for this many combinations, calculation takes really long time as the script utilizes only one CPU. So I reduce this number by appyling a distance cutoff and a stride in addition to a workaround with simultaneous parallel calls to namdenergy.tcl. At the end the number is greatly reduced to around 3000-3500 (depends on the trajectory) and I am actually able to perform the computation in reasonable time using a 24 cores workstation (in about 4 hours, even faster if I increase the stride value). However, there seems to be a bottleneck in speed due to the creation of temporary dcd file, which prevents me from reducing the stride to include more frames in the analysis.
You have to be smart about this. Why do you want an exhaustive comparison -- is it cause you want to see which of the residues are important? for this you could try windowing the number of residues to say 10. See which windows interact. Then resolve the windows further to identify which residues are of interest.
there can be other ways --- i made this up right now -- I am sure you can as well.
Best,
/A
Hi all, I tend to disagree with Ashar's opinion, as in some cases it is possible that for example the repulsive interaction of a pair is compensated (summed up) by the attraction energy of another pair and you might see a small interaction energy as a result of window screening. However, I don't totally disregard it. With some modification e.g applying your knowledge about your protein sites etc you may be able to use it effectively. What algorithm do these docking programs use? Can you adopt that and apply it for your purpose?
Sorry Navid, although you made an astute observation however I don't see how that would be a problem. The summation error would imply that similar charged residues are within 10 residues on the backbone and contribute identical amount of interactions which is not true. The heatmap would be able to show this in a better way. Moreover I wouldn't be surprised if you can pull this off with window sizes larger than 10. Again your choice.
Dear Onur, for this type of NatomXNatom problem, without selecting specific interacting groups, you could write a fortran or C code which will be much more efficient than VMD. You must specify in a separate file atom types with charges as well as with vdW parameters, which will be read by your code. This will be the most time consuming part of the work. You will then first calculate pairwise atom interactions and then you can regroup back to amino acids.
http://www.ks.uiuc.edu/Research/folding/
You can see this link
It may be useful
- Badges
- Science topic