#!/usr/bin/perl

#$n_resachanger 222
##########################################################################
# These values should be adapted to the system under investigation
# ------------------------------------------------------------------------
$offset=4000;      # length of trajectory excluded from calculations (in ps)
$length=6000;    # length of trajectory included in calculations (in ps)
$dyna_step=1000; # Number of dynamics step required to make 1 ps
$n_res=155;       # Number of residues included in calculations
$ref=xtal;   # a nickname for the ref structure MUST MATCH REF IN inp/calc_mean2.inp
##########################################################################

# Do not forget to adapt calc_corre_std.inp to your system.
# i.e. change trajectory, topology, parameters and PSF name, and command :
# CORREL MAXSERIES 300 MAXTIMESTEPS 2000 MAXATOMS 3500
# Also adpt SKIP value in trajectory reading

# length of averaging
$win=@ARGV[0]; # windows length (in ps)

$charmm="/net/bccs/cbu/linux/Charmm/c33b1/exec/gnu/charmm";


$nwin=int($length/$win);
$end_dir="win".$win;

if (-e $end_dir) {
  die "Directory $end_dir already exists. Remove or rename, and restart.\n";
}


$end_mean_dir=$end_dir."/mean_windows";
$mean_out=$end_mean_dir."/out.dat";


print "Calculating fluct $nwin windows of $win ps\n";

mkdir $end_dir, 0755;
mkdir $end_mean_dir, 0755;

for $n (1..$nwin) {
  $dir_name="window_".$n;
  mkdir $dir_name, 0755;
 
  $perl_ini=$dyna_step*($offset+($n-1)*$win);
  $perl_end=$dyna_step*($offset+$n*$win);
  
  print "   Calculating fluctuations around aver. for window $n \n";
  system "sed -e 's/perl-ini/$perl_ini/g' -e 's/perl-end/$perl_end/g' -e 's/perl-rep/$dir_name/g' -e 's/perl-res/$n_res/g' rmsfluct.inp > tmp$win.inp";
  system "$charmm  < tmp$win.inp > tmp$win.out";

  chdir $dir_name;

  print "   Calculating averages  for window $n \n";
# this subroutine will take the data from CA file, look at column that contains
# fluct  for backbone (x column) and keep track of info
  analyse_series();
# this subroutine will take the data from CB file, look at column that contains
# fluct  for side chain (y column) and keep track of info
  analyse_series_b();

# Data from CA is temporarily in out.dat
  open(IN,"out.dat");
  foreach $line (<IN>) {
    @w=split(" ",$line);
    $IC[$w[0]]+=$w[1];
  }
  close(IN);

# Data from CB is temporarily in out.b.dat
  open(IN,"out.b.dat");
  foreach $line (<IN>) {
    @w=split(" ",$line);
    $ICB[$w[0]]+=$w[1];
    $ICX[$w[0]]+=$w[2];
  }
  close(IN);
  
  chdir "..";

  system "mv $dir_name $end_dir";

}

#Finally, an average of the FLUCT is done over all time windows.
# this goes into the mean_window dir, in file out.dat
# there is no header to this file, for ease of plotting.
# the order of data is : residue number, average fluct around mean bb, average fluct
# around mean for side chain. NOT THAT THE OVERALL fluct  is not kept.

open(OUT,">$mean_out");
for $i (1..$n_res) {
    $IC=$IC[$i]/$nwin;
    $ICB=$ICB[$i]/$nwin;
    $ICX=$ICX[$i]/$nwin;
    print OUT "$i $IC $ICB $ICX \n";
}
close(OUT);

print "Calculations completed\n\n";











sub analyse_series {


  $in_name="../fluct.ca.".$perl_ini.".crd"; 
  $last_line=`tail -1 $in_name`;
  @w=split(" ",$last_line);
  $num_line=$w[1];

    open(IN,"$in_name") or die "Could not open file $in_name\n";
    foreach $line (<IN>) {
      @w=split(" ",$line);
      $rms[1][$w[1]]=$w[4];
      print "$w[1] $rms[1][$w[1]] $w[4] \n";
    }
    close(IN);
  
  open(OUT,">out.dat");

  
  for $i (1..$num_line) {
      $sum_sq_i=0; 
	$sum_sq_i+=$rms[$i][$n][1];
      print OUT "$i $rms[1][$i]\n";
      print "$i $rms[1][$i]\n";
  }
  
  close(OUT);

}



sub analyse_series_b {


  $in_name="../fluct.cb.".$perl_ini.".crd"; 
  $last_line=`tail -1 $in_name`;
  @w=split(" ",$last_line);
  $num_line=$w[1];

    open(IN,"$in_name") or die "Could not open file $in_name\n";
    foreach $line (<IN>) {
      @w=split(" ",$line);
      $rms[2][$w[1]]=$w[5];
      $rms[3][$w[1]]=$w[6];
      print "$w[1] $rms[2][$w[1]] $w[5] $rms[3][$w[1]] $w[6] \n";
    }
    close(IN);
#  }
  
  open(OUT,">out.b.dat");
  
  for $i (1..$num_line) {
      $sum_sq_i=0; 
#      for $n (1..$num_line) {
#          print $rms[$i][$n][1];
#          $sum_sq_i+=$rms[$i][$n][1];
#          print "$sum_sq_i $rms[$i][$n][1] \n";
#      }
#      $C[$i][$j]=$prod/sqrt($sum_sq_i*$sum_sq_j);
#       print "$sum_sq_i $num_line $rms[1][$w[0]]  \n";
#      $C[$i]=$sum_sq_i/$num_line;
      print OUT "$i $rms[2][$i] $rms[3][$i]\n";
      print "$i $rms[2][$i] $rms[3][$i]\n";
#    }
  }
  
  close(OUT);

}


sub makedata_xmgr {
  
  open(IN,"out.dat");

  open(OUTp025,">plus025.dat");
  open(OUTp045,">plus045.dat");
  open(OUTp065,">plus065.dat");
  open(OUTp085,">plus085.dat");
  
  open(OUTm025,">minus025.dat");
  open(OUTm045,">minus045.dat");
  open(OUTm065,">minus065.dat");
  open(OUTm085,">minus085.dat");
  
  
  foreach $line (<IN>) {
    @w=split(" ",$line);
    $x=$w[0]-0.5;
    $y=$w[1]-0.5;
    if ($w[2]>0.25 and $w[2]<=0.45) {
      print OUTp025 "$y $x $w[2]\n";
    }
    if ($w[2]>0.45 and $w[2]<=0.65) {
      print OUTp045 "$y $x $w[2]\n";
    }
    if ($w[2]>0.65 and $w[2]<=0.85) {
      print OUTp065 "$y $x $w[2]\n";
    } 
    if ($w[2]>0.85 and $w[2]<=1) {
      print OUTp085 "$y $x $w[2]\n";
    }
    
    
    if ($w[2]<-0.25 and $w[2]>=-0.45) {
      print OUTm025 "$x $y $w[2]\n";
    }
    if ($w[2]<-0.45 and $w[2]>=-0.65) {
      print OUTm045 "$x $y $w[2]\n";
    }
    if ($w[2]<-0.65 and $w[2]>=-0.85) {
      print OUTm065 "$x $y $w[2]\n";
    } 
    if ($w[2]<-0.85 and $w[2]>=-1) {
      print OUTm085 "$x $y $w[2]\n";
    }
    
  }
  
  close(OUTp025);
  close(OUTp045);
  close(OUTp065);
  close(OUTp085);

  close(OUTm025);
  close(OUTm045);
  close(OUTm065);
  close(OUTm085);

  close(IN);
}