AMBER Archive (2000)

Subject: Re: Plotting Hbond info from Carnal

From: samantha hughes (s.j.hughes_at_ic.ac.uk)
Date: Mon Oct 23 2000 - 13:41:10 CDT

Dear All (including Rebecca Perlow!),

Further to my recent request to this esteemed mailgroup re
plotting Carnal hydrogen bond information, I have written a short
Perl script which will extract the hydrogen bonds from the LIS
file generated by Carnal and create seperate files of
(timestep:hydrogen bond distance)data *for the most frequently
occupied hydrogen bonds* . These can include hydrogen bonds which
are not 100% occupied.

I find this script very useful for plotting hydrogen bonds over a
trajectory, without having to establish beforehand which specific
hydrogen bonds you want to monitor. The resulting files can be
easily read into standard graphical packages and plotted, unlike
the LIS, TAB or OUT files generated by Carnal (sorry Bill!)

I hope someone else will find it useful too.

Sam

> Hi All,
>
> Before I reinvent the wheel here, is there any easy/established
> way of plotting frequently occupied hydrogen bonds using a Carnal
> LIS, TAB or OUT file?
>
> I wish to have a series of plots of hydrogen bond length vs
> time, one for each hydrogen bond measured by Carnal (or better,
> only the most frequently occupied ones), preferably without lots
> of fiddling round with awk, grep or perl type programs.
>
> Any ideas/scripts would be appreciated
>
> Thanks for any input,
>
> Sam
>
> -----------------------------------------------------------
> Samantha Hughes
> Department of Chemistry
> Imperial College of Science, Technology and Medicine
> Exhibition Road
> London SW7 2AY
> -----------------------------------------------------------
>#!/usr/sbin/perl
# Program "extract.perl" by S.J.Hughes
# reads in the hbonds.lis file created by carnal and searches for the
# particular hbonds specified below

# prompt the user for information on the data

print "Enter minimal occupancy level as a fractional number\: " ;
$minocc = <STDIN>;
chomp ($minocc);

print "Enter the name of the hbond.lis file: ";
$lisfile = <STDIN>;

# Open data file and initialise hbond type counter

open (LISFILE, $lisfile) || die "Could not open $lisfile -- $!";

$i=1;

# For each line in .lis file, extract hydrogen bond info.
# If bond type is not in hbond array, create a new entry in hbond array, and
# create two new arrays, one for bond length, one for timestep.
# if bond is already in hbond array, index it and add data to relevant timestep and
# bond length arrays

while ( <LISFILE> ) {
  $string = $_;
  @a = split('\s+', $string);
  $locate = index( $string, "(");
  $modstring = substr($string,$locate, 34);
  $descr = " $modstring";
  $maxtime = $a[0];
 
  if ( $name{$descr}) {
    $bondno = $name{$descr} ;
    push(@{eval ts.$bondno}, $a[0]);
    push(@{eval hb.$bondno}, $a[10]);
  } else {
    $name{$descr} = $i ;
    push(@{eval ts.$i}, $a[0]);
    push(@{eval hb.$i}, $a[10]);
    $i++;
  }
}

close (LISFILE);

# Print the list of {timestep, bond length} for each hydrogen bond to a seperate file

while (($bondtyp, $i) = each(%name)) {

   $n=0 ;
   $arr_length = scalar @{eval ts.$i} ;

   if ( $arr_length >= ($minocc * $maxtime)) {
      open (OUTFILE, ">HBOND$i" ) || die "Could not open -- HBOND$i !";
      print OUTFILE "No. of occupied timesteps for hbond \#", $i, " is $arr_length \n" ;
      print OUTFILE "Time \thbond\#",$i, $bondtyp, "\n";

      do {
        print OUTFILE "${eval ts.$i}[$n] \t${eval hb.$i}[$n] \n" ;
        $n++;
      } until ${eval ts.$i}[$n] == undef ;
# print OUTFILE "\n";
        close(OUTFILE);
   }

}