scripts/demux.pl

   1 #!/usr/bin/perl -w
   2
   3 # in: input filename
   4 $in = shift || die("Please specify input filename");
   5 # If your exchange was every N ps and you saved every M ps you can make for
   6 # the missing frames by setting extra to (N/M - 1). If N/M is not integer,
   7 # you're out of luck and you will not be able to demux your trajectories at all.
   8 $extra = shift || 0;
   9 $ndx  = "replica_index.xvg";
  10 $temp = "replica_temp.xvg";
  11
  12 @comm = ("-----------------------------------------------------------------",
  13          "Going to read a file containing the exchange information from",
  14          "your mdrun log file ($in).",
  15          "This will produce a file ($ndx) suitable for",
  16          "demultiplexing your trajectories using trjcat,",
  17          "as well as a replica temperature file ($temp).",
  18          "Each entry in the log file will be copied $extra times.",
  19          "-----------------------------------------------------------------");
  20 for($c=0; ($c<=$#comm); $c++) {
  21     printf("$comm[$c]\n");
  22 }
  23
  24 # Open input and output files
  25 open (IN_FILE,"$in") || die ("Cannot open input file $in");
  26 open (NDX,">$ndx") || die("Opening $ndx for writing");
  27 open (TEMP,">$temp") || die("Opening $temp for writing");
  28
  29
  30 sub pr_order {
  31     my $t     = shift;
  32     my $nrepl = shift;
  33     printf(NDX "%-20g",$t);
  34     for(my $k=0; ($k<$nrepl); $k++) {
  35         my $oo = shift;
  36         printf(NDX "  %3d",$oo);
  37     }
  38     printf(NDX "\n");
  39 }
  40
  41 sub pr_revorder {
  42     my $t     = shift;
  43     my $nrepl = shift;
  44     printf(TEMP "%-20g",$t);
  45     for(my $k=0; ($k<$nrepl); $k++) {
  46         my $oo = shift;
  47         printf(TEMP "  %3d",$oo);
  48     }
  49     printf(TEMP "\n");
  50 }
  51
  52 $nrepl = 0;
  53 $init  = 0;
  54 $tstep = 0;
  55 $nline = 0;
  56 $tinit = 0;
  57 while ($line = <IN_FILE>) {
  58     chomp($line);
  59
  60     if (index($line,"init_t") >= 0) {
  61         @log_line = split (' ',$line);
  62         $tinit = $log_line[2];
  63     }
  64     if (index($line,"Repl") == 0) {
  65         @log_line = split (' ',$line);
  66         if (index($line,"There") >= 0) {
  67             $nrepl = $log_line[3];
  68         }
  69         elsif (index($line,"time") >= 0) {
  70             $tstep = $log_line[6];
  71         }
  72         elsif ((index($line,"Repl ex") == 0) && ($nrepl == 0)) {
  73             # Determine number of replicas from the exchange information
  74             printf("%s\n%s\n",
  75                    "WARNING: I did not find a statement about number of replicas",
  76                    "I will try to determine it from the exchange information.");
  77             for($k=2; ($k<=$#log_line); $k++) {
  78                 if ($log_line[$k] ne "x") {
  79                     $nrepl++;
  80                 }
  81             }
  82         }
  83         if (($init == 0) && ($nrepl > 0)) {
  84             printf("There are $nrepl replicas.\n");
  85
  86             @order = ();
  87             @revorder = ();
  88             for($k=0; ($k<$nrepl); $k++) {
  89                 $order[$k] = $k;
  90                 $revorder[$k] = $k;
  91             }
  92             for($ee=0; ($ee<=$extra); $ee++) {
  93                 pr_order($tinit+$ee,$nrepl,@order);
  94                 pr_revorder($tinit+$ee,$nrepl,@revorder);
  95                 $nline++;
  96             }
  97             $init = 1;
  98         }
  99
 100         if (index($line,"Repl ex") == 0) {
 101             $k = 0;
 102             for($m=3; ($m<$#log_line); $m++) {
 103                 if ($log_line[$m] eq "x") {
 104                     $revorder[$order[$k]] = $k+1;
 105                     $revorder[$order[$k+1]] = $k;
 106                     $tmp = $order[$k];
 107                     $order[$k] = $order[$k+1];
 108                     $order[$k+1] = $tmp;
 109 #           printf ("Swapping %d and %d on line %d\n",$k,$k+1,$line_number);
 110                 }
 111                 else {
 112                     $k++;
 113                 }
 114             }
 115             for($ee=0; ($ee<=$extra); $ee++) {
 116                 pr_order($tstep+$ee,$nrepl,@order);
 117                 pr_revorder($tstep+$ee,$nrepl,@revorder);
 118                 $nline++;
 119             }
 120         }
 121     }
 122 }
 123 close IN_FILE;
 124 close NDX;
 125 close TEMP;
 126
 127 printf ("Finished writing $ndx and $temp with %d lines\n",$nline);