[gmx-users] AGAIN: demixing REMD trajectories
David van der Spoel
spoel at xray.bmc.uu.se
Mon Feb 25 09:23:38 CET 2008
madeleine.kittner at mpikg.mpg.de wrote:
> Hi gromacs users,
>
> I worked a on my problem (see mail below) and got some new questions:
>
> I did REMD using 50 replicas using gromacs version 3.3, while the data
> were saved every 20 ps and exchanges were attempted every 5 ps. I used
> demux.pl to generate replica_temp.xvg and replica_index.xvg. To get
> continuous trajectories I used trjcat from version 3.3.1 (it's not
> available in older versions)
>
> trjcat -f traj*.xtc -demux replica_index.xvg
>
> I solved the segmentation fault problem by compiling gromacs 3.3.1
> differently. But now I'm stuck with the same problem as other people
> before. trjcat produces only one output file trajout.xtc which contains
> only the box parameters and looks like this:
>
> trajout.xtc frame 0:
> natoms= 0 step= 0 time= 0 prec= 0
> box (3x3):
> box[ 0]={ 5.46920e+00, 0.00000e+00, 0.00000e+00}
> box[ 1]={ 1.82296e+00, 5.15645e+00, 0.00000e+00}
> box[ 2]={-1.82292e+00, 2.57794e+00, 4.46567e+00}
> not available: x
> trajout.xtc frame 1:
> natoms= 0 step= 5000 time= 20 prec= 0
> box (3x3):
> box[ 0]={ 5.46920e+00, 0.00000e+00, 0.00000e+00}
> box[ 1]={ 1.82296e+00, 5.15645e+00, 0.00000e+00}
> box[ 2]={-1.82292e+00, 2.57794e+00, 4.46567e+00}
> not available: x
> trajout.xtc frame 2:
> .... and so on.
>
> Is the problem caused by the different time frames in the traj*.xtc files
> and the replica_index.xvg file? Or is it caused by mixing up the different
> gromacs versions? For other people using the demux.pl script seemed to
> solve all the problems but it does not in my case.
>
> Thanks for your help.
> Cheers, Madeleine
>
>> Hi gromacs users,
>>
>> I ran a REMD simulation of 2 peptides in water using 50 replicas. To check
>> the consistency of the REMD system I want to check the structure
>> distributions of each replica. For that I need to demultiplex the
>> (temperature) trajectories. I attempt exchanges every 5ps and save data
>> every 20 ps.
>>
>> I used the demux.pl script to generate replica_index.xvg and
>> replica_temp.xvg which look fine. But if I try to generate continuous
>> trajectories for each replica using
>>
>> trjcat -f traj*.xtc -demux replica_index.xvg
>>
>> I get...
>>
>> Read 50 sets of 17272 points, dt = 5
>>
>> Reading frame 0 time 0.000 Segmentation fault
>>
>> Is the problem caused by the fact that the number of time frames in the
>> trajectories don't match the number of time frames in the replica index
>> file?
>> I read through the earlier postings on this subject, after using demux.pl
>> everyone seemed satisfied. So, I guess it works somehow and I just do
>> something wrong I can't figure out.
>>
did you read the help text in the perl script?
in the latest version (attached) there is an option for your situation.
>> Thanks for your help.
>> Cheers,
>> Madeleine
>>
>>
>
>
> ---
>
> Madeleine Kittner, PhD student
> Department of Theory and Bio-Systems
> Max Planck Institute of Colloids and Interfaces
> Research Campus Golm
> 14424 Potsdam, Germany
>
>
>
>
> _______________________________________________
> 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 van der Spoel, Ph.D.
Molec. Biophys. group, Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. Fax: +4618511755.
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 2008/01/10 09:47:59 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);
More information about the gmx-users
mailing list