Please I have a problem: I concatenated 3 trajectories and got one long trajectory on which I ran the perl script it gives me lets say 60% for hbond 1.
But when I look in each of the 3 trajectories alone, this hbond doesn't even exist in 2 of them. Do you have an idea where the problem might be because I checked if it's well concatenated and it is. And I can't figure out why I have this error. Thanks, Carla On Tue, Oct 26, 2010 at 2:08 PM, Justin A. Lemkul <jalem...@vt.edu> wrote: > > > Carla Jamous wrote: > >> Hi Justin, >> >> Please can you explain this particular comment in your perl script: "# >> There should now be $nres lines left in the file >> # The HB map for the last index is written first (top-down in .xpm file) >> # * Element 0 is index $nres, element 1 is $nres-1, etc." >> >> Is the perl script reading the index file beginning from the bottom? If >> so, why is that? Because in the manual, it says that in the matrix, the >> ordering is identical to that in the index file. >> >> > The order is the same, and the indices are mapped exactly as you would > expect them. The .xpm file is written, however, from top down. So the > final index in hbond.ndx is the top line in the matrix. As such, even > though they are in the same "order," per se, the two files (.xpm and .ndx) > have to be read in opposite directions. > > -Justin > > Thank you. >> >> Carla >> >> >> On Mon, Oct 11, 2010 at 4:31 PM, Justin A. Lemkul <jalem...@vt.edu<mailto: >> jalem...@vt.edu>> wrote: >> >> >> I wrote a Perl script to do a similar task (appended below). >> Perhaps it will be useful to you. I hope it works; I had to hack >> out some things that were specific to my needs and have done only >> limited testing. >> >> -Justin >> >> >> >> #!/usr/bin/perl >> >> # >> # plot_hbmap.pl <http://plot_hbmap.pl> - plot the probability of >> >> finding a particular hydrogen bond >> # based on several input files: >> # 1. coordinate file (for atom naming) - MUST be a .pdb file with >> NO CHAIN IDENTIFIERS >> # 2. hbmap.xpm >> # 3. hbond.ndx (modified to contain only the atom numbers in the >> [hbonds...] section, nothing else) >> # >> >> use strict; >> >> unless(@ARGV) { >> die "Usage: perl $0 -s structure.pdb -map hbmap.xpm -index >> hbond.ndx\n"; >> } >> >> # define input hash >> my %args = @ARGV; >> >> # input variables >> my $map_in; >> my $struct; >> my $ndx; >> >> if (exists($args{"-map"})) { >> $map_in = $args{"-map"}; >> } else { >> die "No -map specified!\n"; >> } >> >> if (exists($args{"-s"})) { >> $struct = $args{"-s"}; >> } else { >> die "No -s specified!\n"; >> } >> >> if (exists($args{"-index"})) { >> $ndx = $args{"-index"}; >> } else { >> die "No -index specified!\n"; >> } >> >> # open the input >> open(MAP, "<$map_in") || die "Cannot open input map file!\n"; >> my @map = <MAP>; >> close(MAP); >> >> open(STRUCT, "<$struct") || die "Cannot open input coordinate file!\n"; >> my @coord = <STRUCT>; >> close(STRUCT); >> >> open(NDX, "<$ndx") || die "Cannot open input index file!\n"; >> my @index = <NDX>; >> close(NDX); >> >> # determine number of HB indices and frames >> my $nres = 0; >> my $nframes = 0; >> for (my $i=0; $i<scalar(@map); $i++) { >> if ($map[$i] =~ /static char/) { >> my $res_line = $map[$i+1]; >> my @info = split(" ", $res_line); >> >> $nframes = $info[0]; >> my @nframes = split('', $nframes); >> shift(@nframes); # get rid of the " >> $nframes = join('', @nframes); >> >> $nres = $info[1]; >> } >> } >> >> print "Processing the map file...\n"; >> print "There are $nres HB indices.\n"; >> print "There are $nframes frames.\n"; >> >> # initialize hashes for later output writing >> # counter $a holds the HB index from hbond.ndx >> my %hbonds; >> for (my $a=0; $a<$nres; $a++) { >> $hbonds{$a+1} = 0; >> } >> >> # donor/acceptor hashes for bookkeeping purposes >> my %donors; >> for (my $b=1; $b<=$nres; $b++) { >> $donors{$b} = 0; >> } >> >> my %acceptors; >> for (my $c=1; $c<=$nres; $c++) { >> $acceptors{$c} = 0; >> } >> >> # clean up the output - up to 18 lines of comments, etc. >> splice(@map, 0, 18); >> >> # remove any "x-axis" or "y-axis" lines >> for (my $n=0; $n<scalar(@map); $n++) { >> if (($map[$n] =~ /x-axis/) || ($map[$n] =~ /y-axis/)) { >> shift(@map); >> $n--; >> } >> } >> >> # There should now be $nres lines left in the file >> # The HB map for the last index is written first (top-down in .xpm >> file) >> # * Element 0 is index $nres, element 1 is $nres-1, etc. >> >> for (my $i=$nres; $i>=1; $i--) { >> # There will be $nframes+2 elements in @line (extra two are " at >> beginning >> # and end of the line) >> # Establish a conversion factor and split the input lines >> my $j = $nres - $i; >> my @line = split('', $map[$j]); >> >> # for each index, write to hash >> for (my $k=1; $k<=($nframes+1); $k++) { >> if ($line[$k] =~ /o/) { >> $hbonds{$i}++; >> } >> } >> } >> >> print "Processing the index file...\n"; >> >> # Open up the index file and work with it >> for (my $n=0; $n<$nres; $n++) { >> my @line = split(" ", $index[$n]); >> $donors{$n+1} = $line[0]; >> $acceptors{$n+1} = $line[2]; >> } >> >> >> # some arrays for donor and acceptor atom names >> my @donor_names; >> my @donor_resn; >> my @acceptor_names; >> my @acceptor_resn; >> >> # Open up the structure file and work with it >> print "Processing coordinate file...\n"; >> foreach $_ (@coord) { >> my @line = split(" ", $_); >> my $natom = $line[1]; >> my $name = $line[2]; >> my $resn = $line[3]; >> my $resnum = $line[4]; >> >> if ($line[0] =~ /ATOM/) { >> unless ($resn =~ /SOL/) { >> for (my $z=1; $z<=$nres; $z++) { >> if ($donors{$z} == $natom) { >> $donor_names[$z] = $name; >> $donor_resn[$z] = join('', $resn, $resnum); >> } elsif ($acceptors{$z} == $natom) { >> $acceptor_names[$z] = $name; >> $acceptor_resn[$z] = join('', $resn, $resnum); >> } >> } >> } >> } >> } >> >> # open a single output file for writing >> open(OUT, ">>summary_HBmap.dat") || die "Cannot open output file!\n"; >> printf(OUT "%10s\t%10s\t%10s\t%10s\%10s\n", "# Donor", " ", >> "Acceptor", " ", "% Exist."); >> >> for (my $o=1; $o<=$nres; $o++) { >> printf(OUT "%10s\t%10s\t%10s\t%10s\t%10.3f\n", $donor_resn[$o], >> $donor_names[$o], $acceptor_resn[$o], $acceptor_names[$o], >> (($hbonds{$o}/$nframes)*100)); >> } >> >> close(OUT); >> >> exit; >> >> >> >> >> Carla Jamous wrote: >> >> Hi everyone, >> >> I tried to analyze the H-bonds in my trajectory with g-hbond and >> I analysed the xpm and ndx file. But now I need to know the >> percentage of existence of each hbond during my trajectory. Is >> there a way to do it with a command line? Or is there a program >> (someone told me there are python programs for analysis of >> gromacs trajectories) to extract this information from the .xpm >> file? >> >> Thank you. >> >> Cheers, >> Carla >> >> >> -- ======================================== >> >> Justin A. Lemkul >> Ph.D. Candidate >> ICTAS Doctoral Scholar >> MILES-IGERT Trainee >> Department of Biochemistry >> Virginia Tech >> Blacksburg, VA >> jalemkul[at]vt.edu <http://vt.edu> | (540) 231-9080 >> >> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin >> >> ======================================== >> -- gmx-users mailing list gmx-users@gromacs.org >> <mailto:gmx-users@gromacs.org> >> >> http://lists.gromacs.org/mailman/listinfo/gmx-users >> Please search the archive at >> http://www.gromacs.org/Support/Mailing_Lists/Search before posting! >> Please don't post (un)subscribe requests to the list. Use the www >> interface or send it to gmx-users-requ...@gromacs.org >> <mailto:gmx-users-requ...@gromacs.org>. >> >> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >> >> >> > -- > ======================================== > > Justin A. Lemkul > Ph.D. Candidate > ICTAS Doctoral Scholar > MILES-IGERT Trainee > Department of Biochemistry > Virginia Tech > Blacksburg, VA > jalemkul[at]vt.edu | (540) 231-9080 > http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin > > ======================================== > -- > gmx-users mailing list gmx-users@gromacs.org > http://lists.gromacs.org/mailman/listinfo/gmx-users > Please search the archive at > http://www.gromacs.org/Support/Mailing_Lists/Search before posting! > Please don't post (un)subscribe requests to the list. Use the www interface > or send it to gmx-users-requ...@gromacs.org. > Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >
-- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists