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 Haridas– Hemanth Haridas2023-04-19 16:19:01 +00:00Commented 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– Anna VDP2023-04-20 08:30:32 +00:00Commented Apr 20, 2023 at 8:30
- $\begingroup$ Why two carbon atoms? $\endgroup$Hemanth Haridas– Hemanth Haridas2023-04-20 11:29:22 +00:00Commented 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– Anna VDP2023-04-20 11:45:42 +00:00Commented Apr 20, 2023 at 11:45
- $\begingroup$ Can you please add an image for the system that you are studying? $\endgroup$Hemanth Haridas– Hemanth Haridas2023-04-20 12:48:18 +00:00Commented 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.
