[AMBER] problem with GB REMD trajectories and ptraj

From: Mark Abraham <Mark.Abraham.anu.edu.au>
Date: Fri, 21 Jan 2011 13:05:08 +1100

Hi,

I've run non-periodic GB REMD on ACE-ALA9-NHE with Amber 11 and run into
problems trying to use ptraj to get common-temperature ensembles from
the replica trajectories. I've included a pretty exhaustive account of
my workflow, because I think this behaviour is indicative of a mismatch
in the way multisander writes non-periodic REMD trajectories and the
assumptions ptraj uses in reading files.

Having equilibrated at the appropriate temp0, I simulated with

  remd simulation
&cntrl
    imin=0,
    ntc=2, ntf=2,
    cut=999, igb=5, saltcon=0.0, gbsa=1, rgbmax=999
    ntpr=50,
    nstlim = 5000, numexchg=5000, dt=0.002,
    ntt=1, temp0=267.00, tautp=1.0,
    ntx=5, irest=1, ntb=0,
    ntwx=50000,
    nscm = 1,
    ntr=0,
&end

(only the lowest of 8 replicas is shown; the others vary only in temp0)

with group file

-O -A -rem 1 -i ala10_linear0.in -p ../../../leap/ala10_linear.top -c
../equil/equil0.rst -o sim0.out -x sim.mdcrd.000 -r sim0.rst -inf sim0.inf
-O -A -rem 1 -i ala10_linear1.in -p ../../../leap/ala10_linear.top -c
../equil/equil1.rst -o sim1.out -x sim.mdcrd.001 -r sim1.rst -inf sim1.inf
-O -A -rem 1 -i ala10_linear2.in -p ../../../leap/ala10_linear.top -c
../equil/equil2.rst -o sim2.out -x sim.mdcrd.002 -r sim2.rst -inf sim2.inf
-O -A -rem 1 -i ala10_linear3.in -p ../../../leap/ala10_linear.top -c
../equil/equil3.rst -o sim3.out -x sim.mdcrd.003 -r sim3.rst -inf sim3.inf
-O -A -rem 1 -i ala10_linear4.in -p ../../../leap/ala10_linear.top -c
../equil/equil4.rst -o sim4.out -x sim.mdcrd.004 -r sim4.rst -inf sim4.inf
-O -A -rem 1 -i ala10_linear5.in -p ../../../leap/ala10_linear.top -c
../equil/equil5.rst -o sim5.out -x sim.mdcrd.005 -r sim5.rst -inf sim5.inf
-O -A -rem 1 -i ala10_linear6.in -p ../../../leap/ala10_linear.top -c
../equil/equil6.rst -o sim6.out -x sim.mdcrd.006 -r sim6.rst -inf sim6.inf
-O -A -rem 1 -i ala10_linear7.in -p ../../../leap/ala10_linear.top -c
../equil/equil7.rst -o sim7.out -x sim.mdcrd.007 -r sim7.rst -inf sim7.inf

The top of sim.mdcrd.000 is

REMD 1 10 50000 267.000
    6.477 8.997 0.146 6.796 9.341 1.130 7.838 9.116
1.357 6.185
    8.791 1.847 6.438 10.800 1.275 5.280 11.090 0.931
7.351 11.661
    1.737 8.273 11.268 1.860 7.264 13.127 1.662 6.523
13.293 0.880
    6.652 13.691 2.926 7.187 13.338 3.808 6.659 14.779
2.861 5.594
   13.456 3.045 8.635 13.785 1.214 9.644 13.087 1.051
8.654 15.074
    0.935 7.803 15.584 1.123 9.732 15.932 0.396 10.342
15.327 -0.275
    8.917 16.984 -0.419 8.248 17.609 0.171 9.632 17.529
-1.036 8.171
   16.526 -1.068 10.708 16.487 1.470 11.889 16.654 1.178
10.074 16.714
    2.595 9.098 16.518 2.768 10.918 17.063 3.750 11.915
16.701 3.498
   10.880 18.593 3.941 9.817 18.811 4.047 11.345 18.903
4.876 11.296
   19.001 3.020 10.430 16.208 5.012 9.464 15.380 4.982
11.278 16.304
    6.000 11.958 17.045 5.908 11.252 15.771 7.330 10.239
15.549 7.665
   12.181 14.576 7.373 13.076 14.767 6.782 12.429 14.312
8.401 11.772
   13.693 6.882 11.735 16.768 8.390 12.396 17.713 7.993
11.552 16.442
    9.655 11.119 15.563 9.901 11.792 17.402 10.758 11.424
18.336 10.333
   10.890 16.975 11.918 11.190 15.995 12.289 10.960 17.651
12.771 9.823
   16.926 11.701 13.277 17.523 11.180 13.669 18.645 11.548
14.136 16.566
   10.878 13.680 15.717 10.576 15.568 16.484 11.013 15.945
16.777 11.993
   15.916 15.030 10.811 15.555 14.675 9.845 17.003 14.972
10.879 15.547
   14.429 11.642 16.313 17.474 10.012 15.775 17.758 8.959
17.449 18.096
   10.481 17.785 17.805 11.388 18.149 19.154 9.761 17.517
20.037 9.673
   19.325 19.577 10.527 20.140 18.853 10.495 19.721 20.494
10.091 19.076
   19.877 11.545 18.565 18.653 8.412 18.297 19.222 7.387
19.076 17.438
    8.395 19.078 16.990 9.300 19.586 16.693 7.269 20.404
17.233 6.792
   20.149 15.515 7.902 19.379 14.782 8.143 20.936 15.120
7.260 20.644
   15.747 8.845 18.566 16.470 6.096 19.015 15.927 5.081
17.311 16.863
    6.228 17.117 17.380 7.074 16.292 16.777 5.136 16.899
16.805 4.230
   15.670 15.363 5.129 15.642 14.967 6.144 14.725 15.378
4.587 16.302
   14.717 4.519 15.196 17.867 5.065 13.986 17.620 4.934
15.605 19.163
    5.156 16.595 19.363 5.172 14.649 20.258 4.983 13.858
20.118 5.720
   15.406 21.481 5.404 16.195 21.790 4.719 14.708 22.316
5.460 15.771
   21.397 6.428 13.915 20.507 3.600 14.386 19.939 2.617
12.811 21.219
    3.623 12.320 21.467 2.775 12.485 21.617 4.492
REMD 1 20 100000 297.620
   13.662 13.278 -0.053 13.534 14.248 -0.532 13.559 14.268
-1.622 14.392
   14.818 -0.175 12.164 14.755 -0.098 11.481 15.524 -0.780
11.730 14.277

ie. there is no box information, as expected with ntb=0.

Now I followed the example here
http://ambermd.org/tutorials/advanced/tutorial7/ to get the set of
structures at my second replica temperature (297.62K) with

ptraj ../../../leap/ala10_linear.top <<EOF
trajin sim.mdcrd.001 remdtraj remdtrajtemp 297.62
trajout ala10_linear_sim1.pdb pdb nobox parse append
EOF

and got output of

   \-/
   -/- PTRAJ: a utility for processing trajectory files
   /-\
   \-/ Version: "AMBER 11.0 integrated" (4/2010)
   -/- Executable is: "ptraj"
   /-\ Running on 1 processor(s)
   \-/ Residue labels:

  ACE ALA ALA ALA ALA ALA ALA ALA ALA ALA
  ALA NHE


PTRAJ: Processing input from "STDIN" ...

PTRAJ: trajin sim.mdcrd.001 remdtraj remdtrajtemp 297.62
   Checking coordinates: sim.mdcrd.001
checkCoordinates(): Could not predict number of frames for AMBER
trajectory file: sim.mdcrd.001
         If this is not a compressed file then there is a problem
Rank: 0 Atoms: 109 FrameSize: 2650 TitleSize: 42 NumBox: 0 Seekable 0

WARNING: setupREMDTRAJ(): remdtraj requested in trajin but this is not
an amber REMD trajectory!
                    Trajectory will be processed as replica traj.

PTRAJ: trajout ala10_linear_sim1.pdb pdb

PTRAJ: Successfully read the input file.
        Coordinate processing will occur until EOF (unknown number of
frames).
        Summary of I/O and actions follows:

INPUT COORDINATE FILES
   File (sim.mdcrd.001) is an AMBER trajectory (with box info)
OUTPUT COORDINATE FILE
   File (ala10_linear_sim1.pdb) is a PDB file

NO ACTIONS WERE SPECIFIED

Processing AMBER trajectory file sim.mdcrd.001


Set 1 .................................................
Set 50 .................................................
Set 100 .................................................
Set 150 .................................................
Set 200 .......................................... 100%


PTRAJ: Successfully read in 242 sets and processed 242 sets.

It was clear from the above output and reading the AmberTools14 code in
ptraj.c that an REMD trajectory can only be recognized if "REMD" occurs
in the second line of the mdcrd file. In my case, this fails because
"REMD" is in the first line. If I manually edit in a fake box line with
"50.0 50.0 50.0" before every REMD line in all the mdcrd files, then the
file is correctly detected as an REMD trajectory:

<snip>
INPUT COORDINATE FILES
   File (sim.mdcrd.000) is an AMBER REMD (new format) trajectory (with
box info) with 238 sets
   Replica processing by temperature will occur.
   8 files total (First index is 000), frames at 297.620000 K will be used.
<snip>

Now things look like they might be good. However, when I go to look at
the temperature ensemble in the PDB file, the tops of the first few
frames are:

MODEL 1
REMARK 1 PDB file generated by ptraj (set 1)
ATOM 1 1HH3 ACE 1 6.477 8.997 0.146 0.1123 1.0000
ATOM 2 CH3 ACE 1 6.796 9.341 1.130 -0.3662 1.7000
ATOM 3 2HH3 ACE 1 7.838 9.116 1.357 0.1123 1.0000
ATOM 4 3HH3 ACE 1 6.185 8.791 1.847 0.1123 1.0000
ATOM 5 C ACE 1 6.438 10.800 1.275 0.5972 1.7000
ATOM 6 O ACE 1 5.280 11.090 0.931 -0.5679 1.4000
<snip>
ATOM 108 HN1 NHE 12 12.320 21.467 2.775 0.2315 1.0000
ATOM 109 HN2 NHE 12 12.485 21.617 4.492 0.2315 1.0000
TER
ENDMDL
MODEL 2
REMARK 1 PDB file generated by ptraj (set 2)
ATOM 1 1HH3 ACE 1 0.000 1.000 20.000 0.1123 1.0000
ATOM 2 CH3 ACE 1 100000.000 297.620 13.662 -0.3662 1.7000
ATOM 3 2HH3 ACE 1 13.278 -0.053 13.534 0.1123 1.0000
ATOM 4 3HH3 ACE 1 14.248 -0.532 13.559 0.1123 1.0000
ATOM 5 C ACE 1 14.268 -1.622 14.392 0.5972 1.7000
ATOM 6 O ACE 1 14.818 -0.175 12.164 -0.5679 1.4000
<snip>
ATOM 108 HN1 NHE 12 13.555 5.019 17.455 0.2315 1.0000
ATOM 109 HN2 NHE 12 11.596 5.165 16.791 0.2315 1.0000
TER
ENDMDL
MODEL 3
REMARK 1 PDB file generated by ptraj (set 3)
ATOM 1 1HH3 ACE 1 11.367 4.439 17.796 0.1123 1.0000
ATOM 2 CH3 ACE 1 10.903 5.816 0.000 -0.3662 1.7000
ATOM 3 2HH3 ACE 1 1.000 30.000150000.000 0.1123 1.0000
ATOM 4 3HH3 ACE 1 267.000 8.315 8.880 0.1123 1.0000
ATOM 5 C ACE 1 1.081 8.228 9.808 0.5972 1.7000
ATOM 6 O ACE 1 1.646 9.128 9.867 -0.5679 1.4000

Clearly there's a frameshift error occuring when actually reading the
data, such that the REMD lines are being read in as coordinates.

I tried editing in my fake box info *after* the REMD line, but this
failed because the trajectory is no longer detected as REMD because REMD
has to be in the second line of the whole file. Subsequently prepending
a dummy line to the whole trajectory succeeded - it was recognized as
REMD, and with box information, but the same frameshift error occurs.

1) Why can't ptraj read the format produced by multisander?
2) How can I fix my mdcrd files so that ptraj can read them?

Thanks,

Mark

_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jan 20 2011 - 18:30:02 PST
Custom Search