[EMAIL PROTECTED] 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    [email protected]
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 [EMAIL PROTECTED]
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.
[EMAIL PROTECTED]       [EMAIL PROTECTED]   http://folding.bmc.uu.se
#!/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,@order);
                pr_revorder($tinit+$ee,$nrepl,@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,@order);
                pr_revorder($tstep+$ee,$nrepl,@revorder);
                $nline++;
            }
        }
    }
}
close IN_FILE;
close NDX;
close TEMP;

printf ("Finished writing $ndx and $temp with %d lines\n",$nline);
_______________________________________________
gmx-users mailing list    [email protected]
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 [EMAIL PROTECTED]
Can't post? Read http://www.gromacs.org/mailing_lists/users.php

Reply via email to