I am new to the md simulation community. I have a peptide which is helical in a bilayer. I want to calculate the angle of the peptide relative to the bilayer with gromacs. For that, I created an index with the peptide residue at the N terminal and the C terminal, so that I can then create a vector with the gmx gangle command. The problem is that for all my trials the computer says me : Inconsistency in user input: Number of positions in selection 2 in the first group not divisible by 3 and I don't understand what It means, I tried to find the explanation on internet but I didn't understand the explanation either. Can anyone of you kindly suggest how can I use this gmx gangle command to compute the angle made by the helical's axis with the normal to the plane? Thank you very much in advance.
9
-
1$\begingroup$ Can you please provide additional information on how is the axis of the helix defined? If yes, I can write an answer with the script that you can use. $\endgroup$– Hemanth HaridasCommented Apr 19, 2023 at 16:19
-
$\begingroup$ I defined the helical axis from the initial two carbon alpha atoms, and two final carbon alpha atoms. Thank you very much for your help. $\endgroup$– Anna VDPCommented Apr 20, 2023 at 8:30
-
$\begingroup$ Why two carbon atoms? $\endgroup$– Hemanth HaridasCommented Apr 20, 2023 at 11:29
-
$\begingroup$ I don't have a specific reason for it, but is there any other way around? Shall I use more Carbon alpha, or all carbon alpha? $\endgroup$– Anna VDPCommented Apr 20, 2023 at 11:45
-
$\begingroup$ Can you please add an image for the system that you are studying? $\endgroup$– Hemanth HaridasCommented Apr 20, 2023 at 12:48
|
Show 4 more comments
1 Answer
$\begingroup$
$\endgroup$
set traj name of trajectory file
set gro name of corresponding gro file
set output [open filename.dat w]
mol load type trj $traj waitfor all
mol addfile gro $gro
set numframes [molinfo top get numframes]
for {set frameC 0} {$frameC < $numframes} {incr frameC} {
molinfo top set frame $frameC
set sel1 [selection string for helix first atom]
set sel2 [selection string for helix second atom]
set vec [vecsub [$sel1 get pos] [$sel2 get pos]]
set nvec [vecnorm $vec]
set theta [lindex $nvec 2]
puts $output "$frameC $theta"
}
I wrote this script for analyzing a NAMD DCD file and have not tested this for a gromacs trajectory. However, this should give you something to start with.
