I have a protein simulated. I want to measure the angle between two domains (domain II & IV) and use the third domain (domain III) as a hinge between them. First, I want to select center of mass of these domains (instead of using the closest atom to center of mass) then get the angle value throughout the trajectory (for every time step). Unfortunately, I could only get the center of mass coordinates of each domain. How can I make the script in VMD?
Hello Yunus Emre Dilek You should be able to define three separate selections to identify the center of masses of the three domains, and then perform the calculation over the trajectory. The main script might look something like this. You would need to iterate this over the entire trajectory.
set domain1 [atomselect top ]
set domain 2 [atomselect top ]
set domain 3 [atomselect top ]
set domain1c [measure center $domain1 weight mass]
set domain2c [measure center $domain2 weight mass]
set domain3c [measure center $domain3 weight mass]
set vec1 [vecsub $domain3c $domain1c]
set vec2 [vecsub $domain3c $domain2c]
set nvec1 [vecnorm $vec1]
set nvec2 [vecnorm $vec2]
set cost [vecdot $nvec1 $nvec2]
puts $cost
Hemanth Haridas, thank you for your answer. I made my own script using my knowledge. I am not sure if I did it correctly but I tried to use the regular formula for the angle between two vectors to get over it. Here's how it looks like without iteration:
set sel1 [atomselect top ]
set sel2 [atomselect top ]
set sel3 [atomselect top ]
set pi [expr {atan(1) * 4}]
set com1 [measure center $sel1 weight mass]
set com2 [measure center $sel2 weight mass]
set com3 [measure center $sel3 weight mass]
set vec3_2 [vecsub $com1 $com3]
set vec3_4 [vecsub $com2 $com3]
set angle_top [vecdot $vec3_2 $vec3_4]
set vec3_2_len [veclength $vec3_2]
set vec3_4_len [veclength $vec3_4]
set angle_bottom [vecdot $vec3_2_len $vec3_4_len]
set angle_cos [expr $angle_top/$angle_bottom]
set angle_rad [expr acos($angle_cos)]
set angle($i.r) [expr $angle_rad*(180/$pi)]
puts $angle($i.r)
This may be a bit of clutter. I also used your script to compare the results. But they differ from each other. How accurate is my script? To give an idea, a brief look at the results:
My script
========
55.98542800468799
56.80759281940429
56.677557130777366
56.25073782787841
57.09028138260412
57.07842772332671
57.49948688907287
57.11029026814625
57.637658875299834
57.41295716017668
58.11547140939341
Hemanth's script
==============
0.5594037339260591
0.5474523309306549
0.5493501632542666
0.5555595257021602
0.5433168600294368
0.5434905344141417
0.5373071612963759
0.5430236462527298
0.5352717267094204
0.5385802547443969
0.528209070153075
Yunus Emre Dilek What does the command "set angle_bottom [vecdot $vec3_2_len $vec3_4_len]" do? And if you are trying to find the normal vector, you should instead use vecnorm $vector. And just to add, my script would give you the cosine of the angle, and you would need to convert it to theta.
- Badges
- Science topic
- Similar topics
- Physical Sciences
- Atomic, Molecular and Optical Physics
- Atoms