VMD Distance Measurements of a Simulation in NAMD
ฝัง
- เผยแพร่เมื่อ 11 ก.ย. 2024
- This video shows you how to a load a simulation of a protein and calculate distances using the software VMD from the University of Illinois (www.ks.uiuc.edu/). The distance.tcl procedure (below) will create several files including the distance between two selections and a histogram of those distances that we analyze in Excel.
Parameter (topology) files necessary for structures with ATP:
toppar_all27_na_nad_ppi.str
This is a set of videos to help use VMD and NAMD for science class or lab projects. The videos step you through visualizing and simulating a protein to better understand how structure and function relate. For more information please see our article
"Molecular mechanics and dynamics characterization of an in silico mutated protein: A stand-alone lab module or support activity for in vivo and in vitro analyses of targeted proteins," Harry Chiang, Lucy C. Robinson, Cynthia J. Brame, and Troy C. Messina, Biochemistry and Molecular Biology Education, 9 Nov 2013.
onlinelibrary.w...
The distance.tcl script is shown below. You will need to search for "lt" and replace with a less than symbol and "gt" with a greater than symbol. These symbols are not allowed in TH-cam descriptions.
###################### distance.tcl ########################
A script for calculating the distance between the center ##
of mass of two selections for every frame that is ##
currently in the top structure of VMD. The script will ##
save the distance for every frame to the file "f_r_out". ##
It will also save a histogram with N_d bins to the file ##
"f_d_out". ##
#################################################################
#################################################################
proc distance {seltext1 seltext2 N_d f_r_out f_d_out} {
set sel1 [atomselect top "$seltext1"]
set sel2 [atomselect top "$seltext2"]
set nf [molinfo top get numframes]
set outfile [open $f_r_out w]
for {set i 0} {$i "lt" $nf} {incr i} {
puts "frame $i of $nf"
$sel1 frame $i
$sel2 frame $i
set com1 [measure center $sel1 weight mass]
set com2 [measure center $sel2 weight mass]
set simdata($i.r) [veclength [vecsub $com1 $com2]]
puts $outfile "$i , $simdata($i.r)"
}
close $outfile
set r_min $simdata(0.r)
set r_max $simdata(0.r)
for {set i 0} {$i "lt" $nf} {incr i} {
set r_tmp $simdata($i.r)
if {$r_tmp "lt" $r_min} {set r_min $r_tmp}
if {$r_tmp "gt" $r_max} {set r_max $r_tmp}
}
set dr [expr ($r_max - $r_min) /($N_d - 1)]
for {set k 0} {$k "lt" $N_d} {incr k} {
set distribution($k) 0
}
for {set i 0} {$i "lt" $nf} {incr i} {
set k [expr int(($simdata($i.r) - $r_min) / $dr)]
incr distribution($k)
}
set outfile [open $f_d_out w]
for {set k 0} {$k "lt" $N_d} {incr k} {
puts $outfile "[expr $r_min + $k*$dr] , $distribution($k)"
}
close $outfile
}
Thanx a lot! :D I have been thinking how to run an MD simulation for about a year and now i got it!!!
I'm glad you found it useful. Thanks!
Thank you so much sir
You are welcome.
it's really very interesting, but i just have one question: between which points of the two amino acids that you use it to measure the distances?
Abdulkarim Najjar The script calculates the center of mass for the two amino acids. It uses a formula like this (upload.wikimedia.org/math/b/e/2/be208f63480ba3b456ac669bea6319cb.png) where R is the center of mass vector, m_i is the mass of each atom at vector position r_i. The distance is the vector length between the two (see below).
If you give the script two single atoms to calculate the distance between, the centers of mass are at each atom, and it will use the atom positions (x, y, z) to calculate a distance = sqrt((x1-x2)^2+(y1-y2)^2+z1-z2)^2).
Thanks a lot for this helpful video and the script to calculate the distance
I have a question. I am only using a single gro file, how is the script needed to modify?
Also, does this works for a uploading a gro file and an xtc file to for the trajectories?
Thanks.
How do you alter the script to measure the min and max distance of two atoms or resid over time?
You can add variables to the script that update on each distance calculation to store min and max, or in Excel open the result from the script as it is and calculate min and max with Excel functions.
how do i get the topology files for the protein?
Hi, Troy. Great video, very helpful. One critical question, how to calculate the shortest distance across the periodic boundary condition (PBC) ?
Open your structure as the top molecule. Then use the commands
set all [atomselect top all]
measure minmax $all
This will give the minimum and maximum x, y, z values. Calculate the differences to find which is the smallest. Here is a discussion that includes a script to calculate the periodic cell. www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2011-2012/2735.html
Thanks a lot. Tutorial was very helpful. How to calculate angle between two helices?
You can calculate the angle made by three atoms. I think the command is measure angle {index1, index2, index3} (www.ks.uiuc.edu/Research/vmd/vmd-1.8.6/ug/node124.html). Or you can create vectors from two atoms on each helix and calculate the dot product -- for example, a = [(h1x-h2x), (h1y-h2y), (h1z-h2z)], where h1x is the x coordinate of atom 1 at one end of helix 1 and h2x is the x coordinate of atom 2 at the opposite end of helix 1. Repeat for b for helix 2. The command is then vecdot a b. (www.ks.uiuc.edu/Research/vmd/vmd-1.7.1/ug/node163.html). The dot product of two vectors is a . b = |a||b| cos(angle), where |a| = sqrt[(h1x-h2x)^2 + (h1y-h2y)^2 + (h1z-h2z)^2].
Hey Troy! Tutorial was very helpful, I have one doubt. How do we compute simulation time ( in nano seconds )? Could you clarify the question?
I assume you mean the simulation time with respect to the distance measurement. It will depend on the settings you use in your simulation. The timestep variable is in femtoseconds (0.0000001 nanoseconds). The frequency of saving to the dcd simulation file is in steps (dcdfreq). The run variable is a total number of steps. Each frame of the simulation that is saved in the dcd file corresponds to t = timestep * dcd, which is units of time/step * step. If you do this correctly, the last time value you calculate is equal to the number of run steps * timestep.
Hi troy,liked your videos please upload a video on how to move a protein using VMD.
Hi Surya, I am not sure what you mean by "move". Do you want to position two molecules? Do you want to simulate the motion of a molecule?
Troy Messina thanks for your reply !I meant to say how to use move script to position two molecules in VMD
If you want to visually align two molecules, you can use the move command (www.ks.uiuc.edu/Training/Tutorials/vmd/tutorial-html/node4.html#SECTION00042300000000000000). If you want to align two similar molecules by best fit, you can use measure-fit commands (www.ks.uiuc.edu/Training/Tutorials/vmd/tutorial-html/node5.html#SECTION00052000000000000000). If you want to do docking studies, you will need another program such as Autodock Vina (vina.scripps.edu/). On Mac and Linux Autodock can be installed as a plugin to VMD.
hi Troy,
I ran the simulation for 5 ns with the dcd frequency 10fs . This means that i have 500000 frames. when i m trying to load this file in vmd . It doesnt take more than 30000 frames. Why is like that. and what can we do to rectify that, Thanks.
It's probably a memory issue. You may only be able to load 30,000 at a time. You can select which frames to load when you open the dcd file. Typically, the dcd file is analyzed through tcl scripting rather than visual inspection.
If I load two files into VMD and I want to find the distance between two different residues on two different protein how do I change the distance command to do that. for example this is what I tried to run:
Main< (post1(3j59and3j5a)(4V75)Middle) 10 % distance "Chain C and resid 20" "Chain d and resid 11" 1 distance.csv histogram.csv
I would also like to do this for a number of different structure that each have their own pdb file.
+KMAMK123 The distance between two different files could result in something completely arbitrary. You will need to align the structures somehow. If they are the same structure from different simulations this is straightforward (www.ks.uiuc.edu/Training/Tutorials/vmd/tutorial-html/node5.html). It looks like you are comparing two different chains of a dimer/trimer/tetramer? If so, you will have to find a way to align them so that they are arranged as they would be when interacting. I'm not sure of a quick answer to that problem. If you can provide more information, I might be able to help more.
Hello Troy, first of all, thank u very much, your videos are very useful! But Troy let's me ask u something, do you have any idea of why the simulation ends, but the dcd file is not produced? I've been having this issue :( Thanks in advance (but this just happen to me with my own files, if I run a file from namd tutorials this works properly)
Hi Diana,
Have you checked your output log to verify there are no errors? Other possibilities are that NAMD is trying to save to a different directory than you expect. Are you running from a terminal or the GUI in VMD. If using the GUI, be sure to use TkConsole to set VMD to the directory where your simulation files are, for example, cd Documents/NAMD/.
Wow Troy, what a fast answer! Thanks again, I'll take a look about what you said. Regards.
I have used gromacs for MD simulation of my protein. So can distance.tcl still be used ?
I don't think so. You might need to use gmx distance (manual.gromacs.org/documentation/2018/onlinehelp/gmx-distance.html)
Sir I need a tutorial video on the execution of MM/PBSA - MM/GBSA on NAMD to calculate the average binding energies of ligand-Protin
I have never used those tools. I am unable to help.
hey troy, its really a great video. its really helped for my studies. i want to calculate hydrogen bonds between protin and ligand in a trajectory. kindly help me in this case.
Hi, I'm not sure what you mean by calculate hydrogen bonds. Do you mean you want to calculate distances where hydrogen bonds may exist? This video shows you how to calculate these distances in a trajectory. You only need to decide which two points for the script to use for distance calculations.
How can we set cell1, cell2 cell3 and cellorigin in configuration file?
achchirangeebi@gmail.com
@@physicssansar802 You can set them in VMD using the NAMD GUI setup. th-cam.com/video/2hEur0xMbiw/w-d-xo.html If you do not follow this method, I would need to know more about your configuration files.
How would this work if I want to plot end-to-end distances?
+FK A You should select the first and last amino acid numbers ("resid 1" for example) . Everything else is the same.
how do i download or form the namd files containing the topology files before starting the simulation?
+devadrita dey sarkar mackerell.umaryland.edu/charmm_ff.shtml#charmm
Ya i opened the link but i am not getting as to where and what should i download for the particular protein as for now i only have the pdb file of my sample protein 2fum but how do i get its topology files from the link?
+Troy Messina also like in your video you have made a separate folder for casein kinase so,if i take my protein 2fum so do i have to make different folder named after my protein and then save the topolgy files over there in the folder like in 1-1-sphere and all just like it is given in the namd tutorial?
I am really new to NAMD and vmd so i am clueless as to how do i start
+devadrita dey sarkar Download the first link from the mackerell site."toppar_c36_feb16.tgz" It is an archive of files that describe the most common atoms and molecules for simulating in Charmm-based programs. You will need to research how to unpack tar and gunzip files on your computer operating system.
There are no rules for how to name and structure folders. It is important that you maintain a structure and naming scheme that is easy for you to identify where files are. The namd tutorials have examples where the topology files go into a folder titled "common". I put my "topology" folder in a folder called "namd simulations". For each protein or mutant I create a new folder in "namd simulations" to store simulations for that specific biomolecule. For example, my folders are "Documents -> namd-simulations->topology" and "Documents -> namd-simulations->casein-kinase". You may want to find a resource for learning about folder and file hierarchy. This will be important for many operations in namd/vmd. Also, do not use spaces in the names of files and folders. It causes a huge headache in namd/vmd.
Thank you so much
Yes, it will measure the distance between any two selections.
@@TroyMessina Amazing. I would like to ask you one more question. What if I want to know which atoms are within a distance x of another atom? It would have to make some modification to the script or would have to start another from scratch.
Thank you for the help
@@andresamessias9310 Finding nearby atoms is a different calculation. The result is a list of atoms and their positions similar to what is in a pdb file. If you do this for a trajectory, you might have thousands of pdb files as a result, or maybe they could be stored in a dcd trajectory. I have never done this calculation on a trajectory. The tcl file would be something like
set sel [atomselect top "within x of resid 10"]
set n [molinfo top get numframes]
for { set i 0 } { $i < $n } { incr i } {
$sel frame $i
$sel update
$sel writepdb withinxof10_$i.pdb
}
Call it in TkConsole with source filename.tcl, and it will run automatically on the top molecule. With more advanced scripting, one could calculate the probability of atoms moving within a distance x of another atom.
Hi Troy, great job and thank you. Is it possible to do it taking atom to atom distance? if yes, can u please describe its syntax, thanks.
Hi Sinha, yes you can calculate distances between atoms. Replace "resid 179" with "index 1234". Of course, 1234 should be the true number of the atom you want to select and subtract 1 because the index begins with 0. If you want to be certain you have chosen the correct atom, you can use boolean logic. "index 1234 and name CA" for example.
Hi Troy, yeah I figured that out after posing that question... very nice tool
i am using charmm force field
Very well-presented video, thanks. However, I have been trying to load the command "source wat_sphere.tcl" for some time now and the Tk reply is "can't read molname: no such variable". I had no problem with the NAMD ubiquitin tutorial, but it doesn't work for my molecular structure. Also, although autoPSF build seems to work OK, the "Add solvation box" facility won't display my molecule in the water box: they sit side by side in VMD. Sorry to bother you, but I'm going round in circles with this. Appreciate your comments.
Hi Simon. Do you want a water box or a water sphere? I prefer the sphere since it is more natural and doesn't require the same kind of boundary conditions as a box that can make the simulation appear as though it is misbehaving. The water sphere requires you to do your own scripting for the simulations though. VMD cannot currently handle water spheres as far as I know.
Are you using the original wat_sphere.tcl or my edited version? The original version is looking for a file called ubq.pdb and ubq.psf in your current Tk directory. It will need to be modified for your particular molecule. My modified version will let you enter the filename. Source it (source wat_sphere.tcl). Then, use the command "addsphere filename". The filename is the name without extensions, i.e., do not put .psf or .pdb. You will need to have psf and pdb files with the same filename. drive.google.com/open?id=0B9hkwleIkrN0bmRaaURYblVZbU0
For the water box, VMD may be behaving correctly. You may need to close all structures and reload only the solvated structure that is saved from the VMD solvation process.
Many thanks for your prompt reply. I am using the original wat-sphere.tcl script, adapted for my molecular structure. I will try your suggestion tomorrow.
It has been very frustrating trying to run the simulation with my molecule...the tutorial went reasonably well, and I was optimistic NAMD would be well-suited for my MSc dissertation. But it really doesn't want to run... I wonder if it's because I'm running an oligonucleotide rather than a protein? I changed the top file but that had no effect.
I can solvate my molecule in a water box as required by the tutorial, but only via the "solvation box" facility. Command line doesn't work, so after that point, the simulation stalls.
Kind regards Simon
PS...it's good to see a video guide that is paced properly. I have made several for my university, and taking a procedure slowly is far more effective than being flash and speedy!
Hello Troy...I'm afraid that didn't work. I can add a water box to my molecule, but the simulation won't go any further. I'm wondering if the problem is loading an unorthodox nucleic acid structure rather than a protein? I have adapted all the relevant scripts, but they are consistently 'not found' in Tk console.
"not found' sounds like a path problem rather than a molecule problem. If you have a na molecule in a water box, you should be able to run from VMD using the NAMD GUI. You'll need to point to all of the relevant topology and molecule files as well as indicate thermodynamic and boundary information. This video shows usage of the GUI. th-cam.com/video/2hEur0xMbiw/w-d-xo.html
nice content
I'm new in VMD and MD. Which file contein the ResID and ResName?
Jahzeel Rosado Thanks!
The resid and resname are in the simulation file (.dcd), but this file is not easily readable because it is binary format. The resid and resname are also in the pdb files. You can open pdb files in a text editor. You can also save a pdb file from frames of a simulation (.dcd). See the thread linked below.
www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2007-2008/3487.html
Great thanks a lot!
Hi
Thank your for the vedio and the great explaination, but I have a problem in the frame loading, I have about 2000000 frames but I can only load 606 frames also my RAM is 32G
Hi, you can load fewer frames and do the analysis on a small subset of frames. You can combine the data of all of the subsets in your program of choice (Excel, Igor, Matlab, etc.). To load fewer frames, see this site www.ks.uiuc.edu/Research/vmd/current/ug/node42.html .
yes this is helpful.. thnx a lot .. but i have question regarding tcl script...the tcl script u have given, this is only for simpler structure of protein .. but i'm having structure of protein with mRNA and i dont have any idea how will i write the script. please help me on this
Troy Messina
Hi Sakshi, I do not work with DNA or RNA. Therefore, I am not certain. I assume you have been able to load and view structures in VMD. You can make distance measurements between any keywords that are indexed in your structure (pdb?) file. For example, if you know the atom "index", you can use that instead of "resid". VMD has a variety of selection keywords. You can see the list of keywords in the Graphics -> Graphical Representations interface. Use the "Selections" tab to view the keywords.
okay thank you Troy
but i'm having 17 conformations, and each conformation contains 3000 atoms (excluded water sphere), so how am gonna measure disance???
The number of conformations doesn't matter. The script goes through "frames", which are the conformations You will get a list of distances for the pair of selections you make for every conformation. You only need to tell the script what selections. You will need to decide which atoms have a distance between them that is important.
Hay,
You said we can scale the steps into fs. How we can do that?
In this video th-cam.com/video/2hEur0xMbiw/w-d-xo.html go to 15:14. Adjust the value "TimeStep (femtoseconds)".