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.

  • 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$ Commented 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 VDP
    Commented Apr 20, 2023 at 8:30
  • $\begingroup$ Why two carbon atoms? $\endgroup$ Commented 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 VDP
    Commented Apr 20, 2023 at 11:45
  • $\begingroup$ Can you please add an image for the system that you are studying? $\endgroup$ Commented Apr 20, 2023 at 12:48

1 Answer 1

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.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .