2022 AMBER tutorial 3 with PDBID 1X70
Contents
"introduction"
In this tutorial, we will be doing molecular dynamics (MD) simulations using Amber 16. This will allow us to obtain much more detailed information about the binding affinities of a ligand with a receptor. Like the DOCK tutorial, we should start with setting up the directory. This can be done with the following command.
mkdir 001_structure 002_parameters 003_leap 004_equil 005_production
"Making the structures"
It is recommended to start fresh with the receptor and ligand files. Do the same process as in the DOCK tutorial. Start with importing pdb code 1X70, then isolate one of the proteins and its associated ligand. With an isolated protein, you can do the following command
Select -> residue -> all nonstandard
then do
actions -> atoms/bonds -> delete
This should isolate the protein and remove all of the nonstandard amino acids. With this, we can save this as a pdb and title it 1X70_fresh.pdb
Moving on to the ligand, we can repeat the same process as above but just isolate the ligand instead. So start with the fresh pdb code, then select the ligand and do
Select -> Invert (selected molecules)
and
actions -> atoms/bonds -> delete
This should isolate the ligand. Then we have to make sure it is charged properly and has hydrogens. We can do that with the following two commands.
Tools -> Structure Editing -> Add H Tools -> Structure Editing -> Add Charge -> (have Amber ff14SB and AM1-BCC selected) -> Ok
Once this is done, save the file as 1x70_lig_wH.mol2
After doing this with both files, upload them to seawulf and place them in the 001_structure directory.
"Simulation parameters"
In this section, we are generating the force field parameters that will be used in future steps. The first step is to swap to the parameters directory, 002_parameters, then use antechamber.
antechamber -i ../001_structure/1x70_lig_wH.mol2 -fi mol2 -o 1x70_ligand_antechamber.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc 1
The -nc tag at the end is the total charge of the ligand, so this may be different depending on the system. In this case, the net charge of the ligand is 1. After this step is completed, we need to parmch2 to generate another file.
parmchk2 -i 1x70_ligand_antechamber.mol2 -f mol2 -o 1x70_ligand.am1bcc.frcmod
"TLeap"
The next step is to create the files needed to simulate the system. Copy the following code into a file named leap.in
Last login: Mon Apr 25 06:08:28 on console (base) noahdean@Noahs-MacBook-Pro ~ % ssh nkdean@login.seawulf.stonybrook.edu nkdean@login.seawulf.stonybrook.edu's password: Reading $DUO_PASSCODE...
Pushed a login request to your device... Success. Logging you in... [nkdean@login2 ~]$ AMS_536 [nkdean@login2 noah]$ cd 1x70_AMBER/ [nkdean@login2 1x70_AMBER]$ cd 003_leap/ [nkdean@login2 003_leap]$ ls 1x70.complex.parm7 1x70.gas.ligand.rst7 1x70.wet.complex.pdb 1x70.gas.complex.pdb 1x70.gas.receptor.parm7 1x70.wet.complex.rst7 1x70.gas.complex.rst7 1x70.gas.receptor.rst7 leap.in 1x70.gas.ligand.parm7 1x70.wet.complex.parm7 leap.log [nkdean@login2 003_leap]$ vim leap.in [nkdean@login2 003_leap]$ cd ../001_structure/ [nkdean@login2 001_structure]$ ls 1x70_fresh.pdb 1x70_ligand_wH.mol2 1x70_receptor_fresh.pdb [nkdean@login2 001_structure]$ vim 1x70_fresh.pdb [nkdean@login2 001_structure]$ cd ../003_leap/ [nkdean@login2 003_leap]$ vim leap.in
#!/usr/bin/sh ###load protein force field source leaprc.protein.ff14SB ###load GAFF force field (for our ligand) source leaprc.gaff ###load TIP3P (water) force field source leaprc.water.tip3p ###load ions frcmod for the tip3p model loadamberparams frcmod.ionsjc_tip3p ###needed so we can use igb=8 model set default PBradii mbondi3
###load protein pdb file rec=loadpdb ../001_structure/1x70_fresh.pdb #bond rec.649.SG rec.762.SG #bond rec.454.SG rec.472.SG #bond rec.444.SG rec.447.SG #bond rec.328.SG rec.339.SG #bond rec.385.SG rec.394.SG ###load ligand frcmod/mol2 loadamberparams ../002_parameters/1x70_ligand.am1bcc.frcmod lig=loadmol2 ../002_parameters/1x70_ligand_antechamber.mol2
###create gase-phase complex
gascomplex= combine {rec lig}
###write gas-phase pdb savepdb gascomplex 1x70.gas.complex.pdb ###write gase-phase toplogy and coord files for MMGBSA calc saveamberparm gascomplex 1x70.complex.parm7 1x70.gas.complex.rst7 saveamberparm rec 1x70.gas.receptor.parm7 1x70.gas.receptor.rst7 saveamberparm lig 1x70.gas.ligand.parm7 1x70.gas.ligand.rst7
###create solvated complex (albeit redundant)
solvcomplex= combine {rec lig}
###solvate the system solvateoct solvcomplex TIP3PBOX 12.0
###Neutralize system addions solvcomplex Cl- 0 addions solvcomplex Na+ 0
#write solvated pdb file savepdb solvcomplex 1x70.wet.complex.pdb
###check the system charge solvcomplex check solvcomplex
###write solvated toplogy and coordinate file saveamberparm solvcomplex 1x70.wet.complex.parm7 1x70.wet.complex.rst7 quit
Note the commented out bond section in this code. This is needed for any disulfide bonds that are found in the receptor. For this system, the disulfide bonds were already made in the .pdb file. They have a tag of SSBOND in that file. In the event that your .pdb file does not contain the disulfide bond information, you will have to manually add those bonds in using the same format as above. With the input file copied, we can use the command
tleap -f leap.in
to create the system. Once this is done, there should be 10 parameter/coordinate files in your directory. These can be loaded into chimera with the following command
Tools -> MD/Ensemble Analysis -> MD Movie --> Select the parameter and coordinate file
It is important to check them in Chimera just to make sure that everything is running properly.
"Equilibration"
Before we can move into the primary production MD simulation, we have to equlibrate the system. Our equilibration process is 9 steps. Copy the following files into your 004_equil directory.
vi 01.min.mdin
Minimize all the hydrogens &cntrl imin=1, ! Minimize the initial structure ntmin=2, ! Use steepest descent Ryota Added maxcyc=5000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask="!@H=", ! atoms to be restrained restraint_wt=5.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
vi 02.equil.mdin
MD simulation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntb=2, ! Constant Pressure ntc=1, ! No SHAKE on bonds between hydrogens dt=0.001, ! Timestep (ps) ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":!@H=", ! atoms to be restrained restraint_wt=5.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
vi 03.min.mdin
Minimize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask="!@H=", ! atoms to be restrained restraint_wt=2.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
vi 04.min.mdin
Minimize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask="!@H=", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
vi 05.min.mdin
Minimize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask="!@H=", ! atoms to be restrained restraint_wt=0.05, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
vi 06.equil.mdin
MD simulation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntb=2, ! Constant Pressure ntc=1, ! No SHAKE on bonds between hydrogens dt=0.001, ! Timestep (ps) ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask="!@H=", ! atoms to be restrained restraint_wt=1.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
vi 07.equil.mdin
MD simulation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask="!@H=", ! atoms to be restrained restraint_wt=0.5, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
vi 08.equil.mdin *See below for specific comments on this input
MD simulation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-729@CA,C,N", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
vi 09.equil.mdin *See below for comments on this input
MD simulation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-728@CA,C,N", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
The restraint mask values for steps 8 and 9 are system-specific. The one for step 8 should include both the receptor and the ligand, and the one for step 9 should include only the receptor. The values for this can be found by looking at the wet.complex.pdb file that was made earlier. Search for the residue labelled LIG, and that will be the number for step 8 (in our case it was 729), and for step 9, it will be 1 number less (removes ligand from the mask). To run the equilibration process, create the following input file titled mdequilibration.sh
#!/bin/sh
#SBATCH --job-name=1HW9_equilibration
#SBATCH --ntasks-per-node=40
#SBATCH --nodes=2
#SBATCH --time=8:00:00
#SBATCH -p long-40core
echo "started Equilibration on 'date' "
parm7="../003_leap/1HW9.wet.complex.parm7"
coords="../003_leap/1HW9.wet.complex"
MDINPUTS=( 01.min 02.equil 03.min 04.min 05.min 06.equil 07.equil 08.equil 09.equil)
for input in ${MDINPUTS[@]}; do
pmemd -O -i ${input}.mdin -o ${input}.mdout -p $parm7 -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r ${input}.rst7
coords=$input
done echo " Finished equilibration on 'date' "
And submit it to the queue with the command
sbatch mdequilibration.sh
"Production MD"
With all of the set up and equilibration steps done, we can move onto the production MD for our system. Switch over to the 005_production directory and create the file 10.prod.mdin
MD simulations &cntrl imin=0, ! Perform MD nstlim=5000000, ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=2, ! SHAKE on for bonds with hydrogen dt=0.002, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=2, ! No force evaluation for bonds with hydrogen ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 2500, ! Write to trajectory file every ntwx steps ntpr= 2500, ! Print to mdout every ntpr steps ntwr= 5000000, ! Write a restart file every ntwr steps cut=8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-728@CA,C,N", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on
Note that the restraint mask should be the same as in step 9 of the equilibration. We can now create the slurm file mdproduction.sh with the following script
!/bin/sh
#SBATCH --job-name=1HW9_prod
#SBATCH --ntasks-per-node=24
#SBATCH --nodes=2
#SBATCH --time=7-00:00:00
#SBATCH -p extended-24core
echo "started Equilibration on 'date' "
parm7="../003_leap/1HW9.wet.complex.parm7"
coords="../004_equil/09.equil"
MDINPUTS=( 10.prod)
for input in ${MDINPUTS[@]}; do
pmemd -O -i ${input}.mdin -o ${input}.mdout -p ${parm7} -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r ${input}.rst7
coords=$input
done
echo "Finished Equilibration on `date` "
Submit this to the queue with
sbatch mdproduction.sh
This will take a while to run, so try not to wait until the last minute to get this job in.
