sACE2 or decoy ± RBD
For example
CHARMM-GUI: Building of Molecular Dynamics Simulation System
Gromacs 2020.4: Minimization, Equilibration and Production
step4.0_minimization.mdp
define = -DPOSRES -DPOSRES_FC_BB=400.0 -DPOSRES_FC_SC=40.0 -DDIHRES -DDIHRES_FC=4.0
integrator = steep
emtol = 1000.0
nsteps = 5000
nstlist = 10
cutoff-scheme = Verlet
rlist = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
coulombtype = pme
rcoulomb = 1.2
;
constraints = h-bonds
constraint_algorithm = LINCS
step4.1_equilibration.mdp
define = -DPOSRES -DPOSRES_FC_BB=400.0 -DPOSRES_FC_SC=40.0 -DDIHRES -DDIHRES_FC=4.0 Remove this
integrator = md
dt = 0.001 Change this to 0.0001
nsteps = 125000
nstxtcout = 5000
nstvout = 5000
nstfout = 5000
nstcalcenergy = 100
nstenergy = 1000
nstlog = 1000
;
cutoff-scheme = Verlet
nstlist = 20
rlist = 1.2
coulombtype = pme
rcoulomb = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
;
tcoupl = Nose-Hoover Change this to v-rescale
tc_grps = SYSTEM
tau_t = 1.0
ref_t = 310
;
constraints = h-bonds
constraint_algorithm = LINCS
;
nstcomm = 100
comm_mode = linear
comm_grps = SYSTEM
;
gen-vel = yes
gen-temp = 310
gen-seed = -1
;
refcoord_scaling = com
step5_production.mdp
integrator = md
dt = 0.002
nsteps = 500000 Change this to 25000000 (for 50 ns) or to 50000000 (for 100 ns)
nstxtcout = 50000 Change this to 500000
nstvout = 50000 Change this to 500000
nstfout = 50000 Change this to 500000
nstcalcenergy = 100 Change this to 100000
nstenergy = 1000 Change this to 100000
nstlog = 1000 Change this to 100000
;
cutoff-scheme = Verlet
nstlist = 20
rlist = 1.2
coulombtype = pme
rcoulomb = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
;
tcoupl = Nose-Hoover
tc_grps = SYSTEM
tau_t = 1.0
ref_t = 310
;
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 5.0
compressibility = 4.5e-5
ref_p = 1.0
;
constraints = h-bonds
constraint_algorithm = LINCS
continuation = yes
;
nstcomm = 100
comm_mode = linear
comm_grps = SYSTEM
;
refcoord_scaling = com
Clustering and Centering of sACE2 + RBD
3D coordinates of proteins can be expressed in different ways because of periodic boundaries used in the MD simulation. We want to keep sACE2 + RBD together near the center of the periodic MD system. This can be achieved with the follwoing protocol.
module load cuda/11.1.1
source $HOME/Library/gromacs/versions/2020.4/cuda-11.1/gcc-8.4/threads/bin/GMXRC
gmx trjconv -f step5_production.xtc -s step5_production.tpr -o temp.xtc -center -ur compact -pbc cluster
(Select Protein for Clustering and Centering and select System for Output)
gmx trjconv -f temp.xtc -s step5_production.tpr -o sACE2-RBD.xtc -pbc mol -ur compact -center
Optional: gmx trjconv -f sACE2-RBD.xtc -s step5_production.tpr -o aligned.xtc -fit rot+trans
cd ..
(You are in the charmm-gui folder.)
tar -zcvf (output_name).tar.gz ./gromacs
(For example, the output name can be sACE2-RBD-date)
sftp
get sACE2-RBD-date.tar.gz
Calculation of Interaction Energies
Prepartion of Input Files for NAMDEnergy Calculation
First, we want to calculate the interaction energy between the sACE2 and SARS-CoV-2 RBD. Secondly, we want to know the contribution of each residue in sACE2 to this interaction energy.
These energies can be calculated using a script in VMD Tk console. In the Tk console, type "source NAMDEnergy.tcl" (video).