1#!/usr/local/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	 "-----------------------------------------------------------------");
20for($c=0; ($c<=$#comm); $c++) {
21    printf("$comm[$c]\n");
22}
23
24# Open input and output files
25open (IN_FILE,"$in") || die ("Cannot open input file $in");
26open (NDX,">$ndx") || die("Opening $ndx for writing");
27open (TEMP,">$temp") || die("Opening $temp for writing");
28
29
30sub pr_order {
31    my $t     = shift;
32    my $nrepl = shift;
33    printf(NDX "%-20.2f",$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
41sub pr_revorder {
42    my $t     = shift;
43    my $nrepl = shift;
44    printf(TEMP "%-20.2f",$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;
57while ($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}
123close IN_FILE;
124close NDX;
125close TEMP;
126
127printf ("Finished writing $ndx and $temp with %d lines\n",$nline);
128