At brief glance things look fine. In my experience, the only time I've
made a mistake during system construction that *wasn't* immediately obvious
when I started running simulations is when I forgot to create disulfide
bonds.
There is a "checkValidity" command in ParmEd that will look over a topology
file and make sure there are no issues with it. This catches things like
missing disulfide bonds and some rather nefarious LEaP bugs that linger in
certain applications. (In particular, I would encourage using this
functionality whenever you use the "bond" command in tleap, since there is
a known bug involving the cross-linking of non-adjacent molecules.)
So something like this:
parmed.py -p CYP2J2_FE_box.prmtop -c CYP2J2_FE_box.inpcrd -i parmed.in
Where the contents of parmed.in are:
checkValidity
You can also be clever and pipe the command through stdin:
echo "checkValidity" | parmed.py -p CYP2J2_FE_box.prmtop -c
CYP2J2_FE_box.inpcrd
Nothing in your input files jump out as incorrect. And since you're
already using ioutfm=1 and ntxo=2, I can't even deliver my canned speech
about using NetCDF instead of ASCII.
HTH,
Jason
On Fri, Oct 9, 2015 at 6:47 AM, Abelak, Kavin <kabelak.rvc.ac.uk> wrote:
> Dear AMBER users,
>
> I was hoping for another set of eyes to gaze upon my input files. My
> simulation works and I just want to check whether I am missing something
> and/or over/under-complicating things as I am the only AMBER user in my
> department and thus no-one can point out any flaws/mistakes in what I have
> done.
>
> I am using AMBER to minimise and equilibrate a homology model of a
> Cytochrome P450 enzyme in explicit solvent. I used parameters in Shahrokh
> et al. 2012 for the heme and Cys-Fe bond. The overall aim is to carry out
> docking studies to identify residues important for docking and/or catalysis.
>
> My leap file input is:
> source leaprc.ff14SB
> source leaprc.gaff
> loadamberparams frcmod.ionsjc_tip3p
> loadamberparams frcmod.hem
> HEM = loadmol2 hem.mol2
> CYP = loadmol2 cyp.mol2
> 2j2 = loadpdb CYP2J2.pdb
> set 2j2.460.28 element "Fe"
> bond 2j2.460.28 2j2.405.8 #BOND CYS-S to HEM-FE
> set 2j2.405 connect0 2j2.405.1
> set 2j2.405 connect1 2j2.405.9
> bond 2j2.404.9 2j2.405.1 #BOND CYS to the protein (new CYP residue)
> bond 2j2.405.9 2j2.406.1 #As above
> savemol2 2j2 CYP2J2_FE.mol2 0
> saveamberparm 2j2 CYP2J2_FE_gas.prmtop CYP2J2_FE_gas.rst7
> addions 2j2 Na+ 0
> solvateoct 2j2 TIP3PBOX 10
> saveamberparm 2j2 CYP2J2_FE_box.prmtop CYP2J2_FE_box.inpcrd
>
>
> The system then has 460 residues and 13637 waters. I then minimise while
> restraining various parts in successive minimisations:
> Initial Minimisation with solute fixed
> &cntrl
> imin=1, irest=0, ntx=1, !minimise, no restart, no vels
> ntmin=1, maxcyc=2000, ncyc=500, !steepest followed by conjugated gradient
> minimisation
> ntb=1, !PBCs
> ntr=1, restraintmask=':1-460', restraint_wt=100.0, !apply restraint mask
> cut=8.0, !set nonbond cutoff at 8.0A
> ntxo=1 !set output to ASCII
> /
>
> CYP2J2: Minimisation 2 with heavy atoms fixed
> &cntrl
> imin=1, irest=0, ntx=1, !minimise, no restart, no vels
> ntmin=1, maxcyc=2000, ncyc=500,!steepest followed by conjugated gradient
> minimisation
> ntb=1, !PBCs
> ntr=1, restraintmask=':1-460 & !.H=', restraint_wt=10.0 !restrain heavy
> atoms
> cut=8.0, !set nonbond cutoff at 8.0A
> ntxo=1 !set output to ASCII
> /
>
> CYP2J2: Minimisation 3 with backbone fixed
> &cntrl
> imin=1, irest=0, ntx=1, !minimise, no restart, no vels
> ntmin=1, maxcyc=2000, ncyc=500,!steepest followed by conjugated gradient
> minimisation
> ntb=1, !PBCs
> ntr=1, restraintmask=':1-459 & .CA,C,O,N,H', restraint_wt=10.0 !restrain
> backbone
> cut=8.0, !set nonbond cutoff at 8.0A
> ntxo=1 !set output to ASCII
> /
>
> CYP2J2: Minimisation with no restraints
> &cntrl
> imin=1, irest=0, ntx=1, !minimise, no restart, no vels
> ntmin=1, maxcyc=10000, ncyc=3000,!steepest followed by conjugated
> gradient minimisation
> ntb=1, !PBCs
> ntr=0, !no restraints
> cut=8.0, !set nonbond cutoff at 8.0A
> ntxo=1 !set output to ASCII
> /
>
> I then heat:
> CYP2J2: Heat while restraining solute
> &cntrl
> imin=0, irest=0, ntx=1, !No minimise, no restart, no vels
> nstlim=100000, dt=0.002, !100,000 steps of 0.002ps
> ntf=2, ntc=2, !Set SHAKE to constrain H
> tempi=0.0, temp0=300.0, !Set initial and final temps
> ntpr=500, ntwx=500, !Write energy and coordinates every 500 steps
> ntb=1, ntp=0, !PBCs, no pressure scaling,
> ntt=3, gamma_ln=2.0, ig=-1, !Langevin dynamics, collision freq, random
> seed
> nmropt=1, !Use varying conditions
> cut=8.0, !Set nonbond cutoff to 8.0A
> ntr=1, restraintmask=':1-460', restraint_wt=10.0, !Apply restraint to
> solute
> ntxo=2 !Output in netCDF format
> /
> &wt type='TEMP0', istep1=0, istep2=90000, value1=0.0, value2=300.0 /
> &wt type='TEMP0', istep1=90001, istep2=100000, value1=300.0, value2=300.0 /
> &wt type='END' /
>
> and equilibrate, again by restraining various bits in successive runs:
> CYP2J2: Equilibration phase
> &cntrl
> imin=0, irest=1, ntx=5, !No minimise, no restart, no vels
> nstlim=10000, dt=0.002, !10,000 steps of 0.002ps
> ntf=2, ntc=2, !Set SHAKE to constrain H
> temp0=300.0, !Set temp
> ntpr=500, ntwx=500, !Write energy and coordinates every 500 steps
> ntwr=10000, !Write restart every 10000 steps
> ioutfm=1, !Output mdcrd and mdvel in netCDF format
> ntb=2, ntp=1, pres0=1.0, taup=2.0, !PBCs auto ntb, pressure scaling,
> ntt=3, gamma_ln=2.0, ig=-1, !Langevin dynamics, collision freq, random
> seed
> cut=8.0, !Set nonbond cutoff t 8.0A
> ntr=1, restraintmask=':1-460', restraint_wt=10.0, !Apply restraint
> ntxo=2 !Output in netCDF format
> /
>
> CYP2J2: Equilibration phase
> &cntrl
> imin=0, irest=1, ntx=5, !No minimise, no restart, no vels
> nstlim=10000, dt=0.002, !10,000 steps of 0.002ps
> ntf=2, ntc=2, !Set SHAKE to constrain H
> temp0=300.0, !Set temp
> ntpr=500, ntwx=500, !Write energy and coordinates every 500 steps
> ntwr=10000, !Write restart every 10000 steps
> ioutfm=1, !Output mdcrd and mdvel in netCDF format
> ntb=2, ntp=1, pres0=1.0, taup=2.0, !PBCs auto ntb, pressure scaling,
> ntt=3, gamma_ln=2.0, ig=-1, !Langevin dynamics, collision freq, random
> seed
> cut=8.0, !Set nonbond cutoff t 8.0A
> ntr=1, restraintmask=':1-459 & .CA,C,O,N,H', restraint_wt=10.0, !Apply
> restraint
> ntxo=2 !Output in netCDF format
> /
>
> CYP2J2: Equilibration phase
> &cntrl
> imin=0, irest=1, ntx=5, !No minimise, no restart, no vels
> nstlim=50000, dt=0.002, !50,000 steps of 0.002ps
> ntf=2, ntc=2, !Set SHAKE to constrain H
> temp0=300.0, !Set temp
> ntpr=500, ntwx=500, !Write energy and coordinates every 500 steps
> ntwr=10000, !Write restart every 10000 steps
> ioutfm=1, !Output mdcrd and mdvel in netCDF format
> ntb=2, ntp=1, pres0=1.0, taup=2.0, !PBCs auto ntb, pressure scaling,
> ntt=3, gamma_ln=2.0, ig=-1, !Langevin dynamics, collision freq, random
> seed
> cut=8.0, !Set nonbond cutoff t 8.0A
> ntxo=2 !Output in netCDF format
> /
>
> And finally the production run:
> CYP2J2: Production Phase
> &cntrl
> imin=0, irest=1, ntx=5, !No minimise, restart, vels
> nstlim=10000000, dt=0.002, !10,000,000 steps of 0.002ps
> ntf=2, ntc=2, !Set SHAKE to constrain H
> temp0=300.0, !Set temp
> ntpr=1000, ntwx=1000, !Write energy and coordinates every 1000 steps
> ntwr=100000, !Write restart every 100000 steps
> ioutfm=1, !Output mdcrd and mdvel in netCDF format
> ntb=2, ntp=1, pres0=1.0, taup=1.0, !PBCs auto ntb, pressure scaling
> ntt=3, gamma_ln=2.0, ig=-1, !Langevin dynamics, collision freq, random
> seed
> cut=8.0, !Set nonbond cutoff to 8.0A
> ntxo=2 !Output in netCDF format
> /
>
> My submission script then looks something like this:
> cd CYP2J2_Sim1
> RunName=CYP2J2_FE_box
>
> cd 001.leap
> cp -t ../002.min $RunName.prmtop $RunName.rst7
>
> cd ../002.min
>
> $AMBERHOME/bin/pmemd.cuda -O -i min1.in -o min1.out -p $RunName.prmtop -c
> $RunName.rst7 -r min1.rst7 -ref $RunName.rst7 -inf ../mdinfo
> $AMBERHOME/bin/pmemd.cuda -O -i min2.in -o min2.out -p $RunName.prmtop -c
> min1.rst7 -r min2.rst7 -ref min1.rst7 -inf ../mdinfo
> $AMBERHOME/bin/pmemd.cuda -O -i min3.in -o min3.out -p $RunName.prmtop -c
> min2.rst7 -r min3.rst7 -ref min2.rst7 -inf ../mdinfo
> $AMBERHOME/bin/pmemd.cuda -O -i min4.in -o min4.out -p $RunName.prmtop -c
> min3.rst7 -r min4.rst7 -inf ../mdinfo
>
> cp -t ../003.heat $RunName.prmtop min4.rst7
> cd ../003.heat
>
> $AMBERHOME/bin/pmemd.cuda -O -i heat.in -o heat_res.out -p
> $RunName.prmtop -c min4.rst7 -r heat_res.rst7 -x heat_res.mdcrd -ref
> min4.rst7 -inf ../mdinfo
>
> cp -t ../004.equil $RunName.prmtop heat_res.rst7
> cd ../004.equil
>
> $AMBERHOME/bin/pmemd.cuda -O -i equil1.in -o equil1.out -p
> $RunName.prmtop -c heat_res.rst7 -r equil1.rst7 -x equil1.mdcrd -ref
> heat_res.rst7 -inf ../mdinfo
> $AMBERHOME/bin/pmemd.cuda -O -i equil2.in -o equil2.out -p
> $RunName.prmtop -c equil1.rst7 -r equil2.rst7 -x equil2.mdcrd -ref
> equil1.rst7 -inf ../mdinfo
> $AMBERHOME/bin/pmemd.cuda -O -i equil3.in -o equil3.out -p
> $RunName.prmtop -c equil2.rst7 -r equil3.rst7 -x equil3.mdcrd -inf ../mdinfo
>
> cp -t ../005.md $RunName.prmtop equil3.rst7
> cd ../005.md
>
> $AMBERHOME/bin/pmemd.cuda -O -i md.in -o md.out -p $RunName.prmtop -c
> equil3.rst7 -r md.rst7 -x md.mdcrd -inf ../mdinfo
>
>
> I just want to point out again that it runs without any discernible
> errors, and the RMSD plots and energies look reasonable. From my limited
> experience in simulations though, I found that this does not always imply a
> healthy simulation. I will appreciate any input and criticisms.
>
> Many thanks,
> Kavin
>
> -----------------------------
> Kavin Abelak
> PhD Student
> Comparative Biomedical Sciences
> Royal Veterinary College
> 4 Royal College St
> London NW1 0TU
>
>
>
> [http://www.rvc.ac.uk/media/default/images/social-media/change.gif]
>
> [RVC Logo - link to RVC Website]<http://www.rvc.ac.uk> [Twitter icon -
> link to RVC (Official) Twitter] <http://twitter.com/RoyalVetCollege>
> [Facebook icon - link to RVC (Official) Facebook] <
> http://www.facebook.com/theRVC> [YouTube icon - link to RVC YouTube] <
> http://www.youtube.com/user/RoyalVetsLondon?feature=mhee> [Pinterest
> icon - link to RVC Pinterest] <http://pinterest.com/royalvetcollege/>
> [Instagram icon - link to RVC Instagram] <
> http://instagram.com/royalvetcollege>
>
> This message, together with any attachments, is intended for the stated
> addressee(s) only and may contain privileged or confidential information.
> Any views or opinions presented are solely those of the author and do not
> necessarily represent those of the Royal Veterinary College (RVC). If you
> are not the intended recipient, please notify the sender and be advised
> that you have received this message in error and that any use,
> dissemination, forwarding, printing, or copying is strictly prohibited.
> Unless stated expressly in this email, this email does not create, form
> part of, or vary any contractual or unilateral obligation. Email
> communication cannot be guaranteed to be secure or error free as
> information could be intercepted, corrupted, amended, lost, destroyed,
> incomplete or contain viruses. Therefore, we do not accept liability for
> any such matters or their consequences. Communication with us by email will
> be taken as acceptance of the risks inherent in doing so.
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Oct 09 2015 - 06:30:06 PDT
