Hello,
I have been trying over the last 2 weeks to use the ambermd engine with the charmm forcefield with minimal sucess and I need advice on how to best approach this problem.
My system is Dimeric Hemoglobin (HbI) in the Deoxy state or 4SDH. Before preceding with constructing my system I used modeller to predict the coordinates of the missing N-termini (ie Proline), then I used H++ to predict the protonation states of the titratable residues and changed the 3 letter code to match the correct charge ie HIS -> HSD. Made sure the histidine that links to the Iron center was HSD. I then split all chains to separate pdb files: chain A and B, Heme A and B, and, waters, producing 5 separate pdb files.
I have been successful in producing psf and crd files in charmm and perform minimization of input structure with missing atoms followed by solvation then minimization of the new system. The resulting psf and crd files were then read into the ParmEd program using the chamber command. Noting that I could only get this to work for a cubic box but not a rhombic dodecahedron, I was under the assumption that the rhdo can be defined by x=y=z a=c=60deg b=90deg and that I would have to check the defined box size in both prmtop and inpcrd but received an exeption error regarding mismatched box size.
With the output of ParmEd using parmout I am able to minimize my system, perform NVT heating over 20ps, I however cannot heat over a longer period of time and have my temperature not blow up, 30ps is enough to see the start of the temp blow up from my temp of 300K at 25 ps to 450K at 30ps to Nan at 35ps. I have tried berendsen's thermostat with tautp=1 and 2 and Andersen's thermostat vrand=1000 step and Langevin dynamics with gamma_ln=1 and 2. I have tried heating to 300K with and without posres, I have tried with ntr=1, restraint_wt=99.9, 9 and 2.
If I continue with NPT using the berenden barostat with taup=1 or 2ps and all combinations of the thermostat params described for NVT the system quickly ~50ps behaves erratically with respect to pressure and quickly pressure blows up. Oddly the temperature stays the same, 300K. If I use ntr=1 or position restraints I find myself in the same hole.
For all simulation I have tried dt=2fs, 1fs, 0.5fs; restraining all heavy atoms, restraining heme and protein backbone, restraining backbone and none.
Any and all help is greatly appreciated and my scripts for charmm (free version) and AMBERMD v16 are appended to the end. I did not include multiple versions with differing params usch as showing 2 files for the different gamma_ln used, etc.
Thank you,
Korey
psf and crd generation scripts:
top/par stream:
file = heme_toppar.str
* Load All toppar files for HBI
*
read rtf card name toppar/top_all36_prot.rtf
read para card flex name toppar/par_all36_prot.prm
stream toppar/toppar_water_ions.str
stream toppar/stream/prot/toppar_all36_prot_heme.str
stream toppar/protpatch_protein_toppar36.str
step1 - create crd and psf for prot/heme/wat:
create individual psf/crd files for prot/heme and water
combine psf/crd files
step2 - solvate
use sphere of waters that have been mini/equil
remove waters that overlap prot/heme/wat from step1
remove all waters outside of box
step3 - minimize
minimize with SD 100 and ABNER 1000 with tol <1.0
Parmed
>chamber -top top_all36_prot.rtf -param par_all36_prot.prm -str protpatch_protein_toppar36.str -str toppar_water_ions.str -str toppar_all36_prot_heme.str -psf charmm_system.psf -crd charmm_system.crd -box 80,80,80,90,90,90
AMBER input scripts:
Minimize:
MINIMIZE
&cntrl
imin=1,
ntx=1,
irest=0,
maxcyc=2000,
ncyc=50,
ntpr=50,
cut=12.0,
ntr=1,
restraint_wt=99
restraintmask='.* & !.H='
watnam='TIP3'
owtnm='OH2'
/
HEAT:
Heat
&cntrl
imin=0,
ntx=1,
irest=0,
nstlim=20000,
dt=0.001,
ntf=2,
ntc=2,
tempi=0.0,
temp0=300.0,
ntpr=100,
ntwx=100,
cut=9.0,
ntb=1,
ntp=0,
ntt=1,
ntr=1,
restraint_wt=99
restraintmask='.* & !.H='
gamma_ln=0,
nmropt=0,
ig=-1,
watnam='TIP3'
owtnm='OH2'
/
&wt type='TEMP0', istep1=0, istep2=18000, value1=0.0, value2=300.0 /
&wt type='TEMP0', istep1=18001, istep2=20000, value1=300.0, value2=300.0 /
&wt type='END' /
NPT:
Equilibrate
&cntrl
imin=0,
ntx=5,
irest=1,
nstlim=50000,
dt=0.001,
ntf=2,
ntc=2,
temp0=300.0,
ntpr=2500,
ntwx=2500,
cut=9.0,
ntb=2,
ntp=1,
ntt=1,
ntr=1,
restraint_wt=99
restraintmask='.CA,C,N,O'
gamma_ln=0,
ig=-1,
watnam='TIP3'
owtnm='OH2'
/
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Oct 06 2016 - 11:30:02 PDT
