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.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Aug 19 2009 - 22:17:45 PDT
