??: Re: [gmx-users] demux.pl

[David van der Spoel]

Mark Abraham wrote:
> David van der Spoel wrote:
>> xi zhao wrote:
>>> Sorry, but in  source of v3.3.2 and 3.3.1 of the website , there is 
>>> not demux.pl
>>
>> src/contrib/scripts
You're dead right.

WIll fix this for the next release. It is attached to this mail.

> 
> Actually, Xi Zhao is right here, David. Neither of the 3.3.1 or 3.3.2 
> source distributions has either this subdirectory or demux.pl anywhere 
> else.
> 
> Presumably the CVS distribution will have it, if you want to try 
> downloading that, Xi. See http://wiki.gromacs.org/index.php/CVS_HowTo
> 
> Mark
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the www 
> interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php


-- 
David.
________________________________________________________________________
David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,  	75124 Uppsala, Sweden
phone:	46 18 471 4205		fax: 46 18 511 755
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-------------- next part --------------
#!/usr/bin/perl -w
# $Id: demux.pl,v 1.1.2.4 2007/01/16 12:03:48 spoel Exp $

# in: input filename
$in = shift || die("Please specify input filename");
# If your exchange was every N ps and you saved every M ps you can make for
# the missing frames by setting extra to (N/M - 1). If N/M is not integer,
# you're out of luck and you will not be able to demux your trajectories at all.
$extra = shift || 0;
$ndx  = "replica_index.xvg";
$temp = "replica_temp.xvg";

@comm = ("-----------------------------------------------------------------",
	 "Going to read a file containing the exchange information from",
	 "your mdrun log file ($in).", 
	 "This will produce a file ($ndx) suitable for",
	 "demultiplexing your trajectories using trjcat,",
	 "as well as a replica temperature file ($temp).",
	 "Each entry in the log file will be copied $extra times.",
	 "-----------------------------------------------------------------");
for($c=0; ($c<=$#comm); $c++) {
    printf("$comm[$c]\n");
}

# Open input and output files
open (IN_FILE,"$in") || die ("Cannot open input file $in");
open (NDX,">$ndx") || die("Opening $ndx for writing");
open (TEMP,">$temp") || die("Opening $temp for writing");


sub pr_order {
    my $t     = shift;
    my $nrepl = shift;
    printf(NDX "%-8g",$t);
    for(my $k=0; ($k<$nrepl); $k++) {
	my $oo = shift;
	printf(NDX "  %3d",$oo);
    }
    printf(NDX "\n");
}

sub pr_revorder {
    my $t     = shift;
    my $nrepl = shift;
    printf(TEMP "%-8g",$t);
    for(my $k=0; ($k<$nrepl); $k++) {
	my $oo = shift;
	printf(TEMP "  %3d",$oo);
    }
    printf(TEMP "\n");
}

$nrepl = 0;
$init  = 0;
$tstep = 0;
$nline = 0;
$tinit = 0;
while ($line = <IN_FILE>) {
    chomp($line);
    
    if (index($line,"init_t") >= 0) {
	@log_line = split (' ',$line);
	$tinit = $log_line[2];
    }
    if (index($line,"Repl") == 0) {
	@log_line = split (' ',$line);
	if (index($line,"There") >= 0) {
	    $nrepl = $log_line[3];
	}
	elsif (index($line,"time") >= 0) {
	    $tstep = $log_line[6];
	}
	elsif ((index($line,"Repl ex") == 0) && ($nrepl == 0)) {
            # Determine number of replicas from the exchange information
	    printf("%s\n%s\n",
		   "WARNING: I did not find a statement about number of replicas",
		   "I will try to determine it from the exchange information.");
	    for($k=2; ($k<=$#log_line); $k++) {
		if ($log_line[$k] ne "x") {
		    $nrepl++;
		}
	    }
	}
	if (($init == 0) && ($nrepl > 0)) {
	    printf("There are $nrepl replicas.\n");

	    @order = ();
            @revorder = ();
	    for($k=0; ($k<$nrepl); $k++) {
		$order[$k] = $k;
                $revorder[$k] = $k;
	    }
	    for($ee=0; ($ee<=$extra); $ee++) {
		pr_order($tinit+$ee,$nrepl, at order);
		pr_revorder($tinit+$ee,$nrepl, at revorder);
		$nline++;
	    }
	    $init = 1;
	}

	if (index($line,"Repl ex") == 0) {
	    $k = 0;
	    for($m=3; ($m<$#log_line); $m++) {
		if ($log_line[$m] eq "x") {
		    $revorder[$order[$k]] = $k+1;
		    $revorder[$order[$k+1]] = $k;
		    $tmp = $order[$k];
		    $order[$k] = $order[$k+1];
		    $order[$k+1] = $tmp;
#	    printf ("Swapping %d and %d on line %d\n",$k,$k+1,$line_number); 
		}
		else {
		    $k++;
		}
	    }
	    for($ee=0; ($ee<=$extra); $ee++) {
		pr_order($tstep+$ee,$nrepl, at order);
		pr_revorder($tstep+$ee,$nrepl, at revorder);
		$nline++;
	    }
	}
    }
}
close IN_FILE;
close NDX;
close TEMP;

printf ("Finished writing $ndx and $temp with %d lines\n",$nline);