AMBER Archive (2009)

Subject: [AMBER] REMD using hybrid solvent model

From: Morad Mustafa (alawneh_at_uga.edu)
Date: Wed Aug 12 2009 - 14:05:40 CDT

Hi,

I have carried out REMD using hybrid solvent model based on the tutorial from the web and the
archive notes. In my simulation, I use Amber 10.
I have an issue with reading the REMD trajectories. The ptraj reads only one frame out of 500
frames.
How can I analyze those types of REMD trajectories if ptraj does not read them?

This is the script I used for reading the REMD trajectories:
#!/bin/csh

set Replica_First = 1
set Replica_Last = 1
set PrmTop = "../../Modeling/EGFR-WT_active.prmtop"
set Process = "REMD"

set Replica = $Replica_First

while ($Replica <= $Replica_Last)
        set RepID = `printf "%3.3i" $Replica`
        set TrajIn = "../${Process}-${RepID}.mdcrd"
        set TrajOut = "${Process}-${RepID}-Test.mdcrd"
        ptraj $PrmTop << EOF
                trajin $TrajIn remdtraj
                trajout $TrajOut
EOF
        @ Replica = ($Replica + 1)
end

Here is part of the outputs:
PTRAJ: Processing input from "STDIN" ...

PTRAJ: trajin ../REMD-001.mdcrd remdtraj
  Checking coordinates: ../REMD-001.mdcrd

PTRAJ: trajout REMD-001-Test.mdcrd

PTRAJ: Successfully read the input file.
       Coordinate processing will occur on 500 frames.
       Summary of I/O and actions follows:

INPUT COORDINATE FILES
  File (../REMD-001.mdcrd) is an AMBER trajectory (with box info) with 500 sets

OUTPUT COORDINATE FILE
  File (REMD-001-Test.mdcrd) is an AMBER trajectory (with box info)

NO ACTIONS WERE SPECIFIED

Processing AMBER trajectory file ../REMD-001.mdcrd

Set 1 .
WARNING in readAmberTrajectory(): Set #2 is corrupted (REMD 1 4 )...

PTRAJ: Successfully read in 1 sets and processed 1 sets.
       Dumping accumulated results (if any)

This is the REMD.in
Replica Exchange Molecular Dynamics Run
 &cntrl
        ! Minimization Flags
                imin = 0,
                maxcyc = 1,
                ncyc = 500,
                ntmin = 1,
        ! Input Format
                ntx = 1,
                irest = 0,
        ! Output Format and Frequency
                ntpr = 500,
                ntave = 0,
                ntwr = 5000,
                iwrap = 1,
                ntwx = 1000,
                ntwv = 0,
                ntwe = 0,
                ntwprt = 0,
                idecomp = 0,
        ! Harmonically Restrained Atoms
                ntr = 0,
                !restraint_wt = ,
                !restraintmask= ,
        ! MD Parameters
                nstlim = 500,
                numexchg= 1000,
                nscm = 500,
                t = 0.0,
                dt = 0.002,
                nrespa = 2,
        ! Temperature Regulation
                ntt = 3,
                temp0 = 298.15,
                tempi = 298.15,
                ig = 12993,
                tautp = 1.0,
                gamma_ln= 1.0,
        ! Pressure Regulation
                ntp = 0,
                pres0 = 1.0,
                taup = 2.0,
        ! SHAKE Parameters
                ntc = 2,
        ! Potential Function Parameters
                ntf = 2,
                ntb = 1,
                cut = 12.0,
                scnb = 2.0,
                scee = 1.2,
                nsnb = 25,
                igb = 0,
                hybridgb= 5,
                numwatkeep = 1000,
                ipol = 0,
                ifqnt = 0,
 /
 &ewald
                nfft1 = 128,
                nfft2 = 128,
                nfft3 = 128,
                order = 6,
                verbose = 0,
                ew_type = 0,
                nbflag = 1,
                skinnb = 2.0,
                netfrc = 1,
                vdwmeth = 1,
                column_fft = 1,
 /

This is how the topology file was prepared:
# Load a force field
source "/usr/local/amber10/dat/leap/cmd/leaprc.ff03.r1"

# Use implicit solvent model
set default PBradii bondi

# Load in a pdb file.
Seq_List = {
NGLY GLU ALA PRO ASN GLN ALA LEU LEU ARG ILE LEU LYS
GLU THR GLU PHE LYS LYS ILE LYS VAL LEU GLY SER GLY
ALA PHE GLY THR VAL TYR LYS GLY LEU TRP ILE PRO GLU
GLY GLU LYS VAL LYS ILE PRO VAL ALA ILE LYS GLU LEU
ARG GLU ALA THR SER PRO LYS ALA ASN LYS GLU ILE LEU
ASP GLU ALA TYR VAL MET ALA SER VAL ASP ASN PRO HIS
VAL CYS ARG LEU LEU GLY ILE CYS LEU THR SER THR VAL
GLN LEU ILE THR GLN LEU MET PRO PHE GLY CYS LEU LEU
ASP TYR VAL ARG GLU HIS LYS ASP ASN ILE GLY SER GLN
TYR LEU LEU ASN TRP CYS VAL GLN ILE ALA LYS GLY MET
ASN TYR LEU GLU ASP ARG ARG LEU VAL HIS ARG ASP LEU
ALA ALA ARG ASN VAL LEU VAL LYS THR PRO GLN HIS VAL
LYS ILE THR ASP PHE GLY LEU ALA LYS LEU LEU GLY ALA
GLU GLU LYS GLU TYR HIS ALA GLU GLY GLY LYS VAL PRO
ILE LYS TRP MET ALA LEU GLU SER ILE LEU HIS ARG ILE
TYR THR HIS GLN SER ASP VAL TRP SER TYR GLY VAL THR
VAL TRP GLU LEU MET THR PHE GLY SER LYS PRO TYR ASP
GLY ILE PRO ALA SER GLU ILE SER SER ILE LEU GLU LYS
GLY GLU ARG LEU PRO GLN PRO PRO ILE CYS THR ILE ASP
VAL TYR MET ILE MET VAL LYS CYS TRP MET ILE ASP ALA
ASP SER ARG PRO LYS PHE ARG GLU LEU ILE ILE GLU PHE
SER LYS MET ALA ARG ASP PRO GLN ARG TYR LEU VAL ILE
GLN GLY ASP GLU ARG MET HIS LEU PRO SER PRO THR ASP
SER ASN PHE TYR ARG ALA LEU MET ASP GLU GLU ASP MET
ASP ASP VAL VAL ASP ALA ASP GLU TYR LEU CILE }
# PRO CGLN }
# GLN GLY }

Protein = loadPdbUsingSeq "EGFR-WT_active_equil.pdb" Seq_List

# Neutralize the protein
addIons Protein K+ 0

# Solvate the protein with the right water model
solvateBox Protein TIP3PBOX 5.0

# Save the prmtop and inpcrd files
saveAmberParm Protein "EGFR-WT_active.prmtop" "EGFR-WT_active.mdcrd"
savePdb Protein "EGFR-WT_active.pdb"

quit

Regards,

Dr. Morad Mustafa

_______________________________________________
AMBER mailing list
AMBER_at_ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber