RPQS

Reference-potential with quantum sampling.
Described in JCTC (https://doi.org/10.1021/acs.jctc.6b01217)


Here are two scripts to run RPQS calculations for ligand binding (Meiting 240113):

The first script (pmemd-ti-1step-mm.in) involves the calculation of ΔG_free using TI simulation at the MM level.

This script utilizes 11 windows and is implemented with Amber's pmemd.


The second script(sand_mm_to_DFTB.sh) pertains to the computation of the free energy difference from MM to DFTB using TI simulation. For this script, the sander module in Amber has been employed, and it also encompasses 11 windows.


pmemd-ti-1step-mm.in

#!/bin/bash
wd=`pwd`
pre0='ti'
Y='$SLURM_SUBMIT_DIR'

for X in 0 1 2 3 4 5 6 7 8 9

do

cat << EOF > mdin_mini_l${X}
density minlibration
 &cntrl
  imin = 1,     ntx = 1,
  maxcyc=10000,  ntmin = 2,
  ntpr = 100,
  ntf = 1,      ntc = 2,
  ntb = 1,      cut = 12.0,
  icfe=1,       clambda = 0.${X},
  ifsc=1,
  timask1=':1',timask2=':2',
  scmask1=':1@H6',scmask2=':2@O1,H6',
 /
EOF

cat <<EOF > mdin_heat_l${X}
heat
 &cntrl
  imin=0,irest=0,ntx=1,
  nstlim=50000,dt=0.001,
  ntc=2,ntf=1,ifsc=1,
  cut=12.0,ntb=1,ntp=0,
  ntpr=1000,ntwr=1000,ntwx=0,
  ntt=3,gamma_ln=2.0,
  tempi=0.0,temp0=298.0,
  nmropt=1,icfe=1,clambda=0.${X},
  timask1=':1',timask2=':2',
  scmask1=':1@H6',scmask2=':2@O1,H6',
 /
 &wt TYPE='TEMP0',istep1=0,istep2=50000,value1=0.1,value2=300.0, /
 &wt TYPE='END' /
EOF

cat << EOF > mdin_equi_l${X}
density equilibration
 &cntrl
  imin = 0,     ntx = 5,        irest = 1,
  ntpr = 2500,  ntwr = 10000,   ntwx = 0,
  ntf = 1,      ntc = 2,
  ntb = 1,      cut = 12.0,
  nstlim = 200000,      dt = 0.001,
  temp0 = 298.0,        ntt = 3,        gamma_ln = 5,
  icfe=1,       clambda = 0.${X},
  ifsc=1,
  timask1=':1',timask2=':2',
  scmask1=':1@H6',scmask2=':2@O1,H6',
/
EOF

cat << EOF > mdin_prod_l${X}
NPT production
 &cntrl
  imin = 0,     ntx = 5,        irest = 1,
  ntpr = 1000,  ntwr = 1000,    ntwx = 1000,
  ntf = 1,      ntc = 2,
  ntb = 2,      cut = 12.0,
  nstlim = 5000000,     dt = 0.002,
  temp0 = 298.0,        ntt = 3,        gamma_ln = 2,
  ntp = 1,      pres0 = 1.0,    taup = 2.0,
  icfe=1,       clambda = 0.${X},
  ifsc=1,       ifmbar=1,
  bar_l_min=0.0,bar_l_max=1.0,
  timask1=':1',timask2=':2',
  scmask1=':1@H6',scmask2=':2@O1,H6',
  bar_intervall=1000,ioutfm=1,
 /
EOF

cat << EOF > run.pbs.${X}
#!/bin/bash
#SBATCH -N 1
#SBATCH --job-name=P-L-${X}
#SBATCH -p tetralith
#SBATCH --exclusive
#SBATCH -A snic2022-3-19
#SBATCH -t 168:00:00


module purge
module load buildtool-easybuild/4.2.2-nsc2035f5dd0
module --ignore-cache load "binutils/.2.28"
module  load icc/.2018.1.163-GCC-6.4.0-2.28
module load  impi/.2018.1.163
module load SciPy-bundle/2020.03-Python-3.8.6-nsc1
module load GCC/11.3.0
module load  OpenMPI/4.1.4
module load SciPy-bundle/2022.05-nsc1
module load matplotlib/3.5.2-nsc1
module load Amber/16-nsc1-intel-2017.u7-bare

cd $Y
#mpprun pmemd.MPI  -i mdin_mini_l${X} -O -o mini_l${X}.out  -p ti.parm7 -c ti.rst7          -r mini${X}.rst
#mpprun pmemd.MPI  -i mdin_heat_l${X} -O -o heat_l${X}.out  -p ti.parm7 -c mini${X}.rst     -r heat${X}.rst
#mpprun pmemd.MPI  -i mdin_equi_l${X} -O -o equi_l${X}.out  -p ti.parm7 -c heat${X}.rst     -r equi${X}.rst
mpprun pmemd.MPI  -i mdin_prod_l${X} -O -o prod_l${X}.out  -p ti.parm7 -c equi${X}.rst     -r prod${X}.rst -x prod${X}.crd
EOF
sbatch   run.pbs.${X}
done

sand_mm_to_DFTB.sh

#!/bin/bash
wd=`pwd`
pre0='gal-l03'
pre1='gal-l03'

for X in 0 1 2 3 4  5 6 7 8 9   ##number of windows, could be changed according to Amber mannual
do

cat << EOF > group_min_l${X}
-O -i mdin_min_v0_l${X} -o ${pre0}_min_v0_l${X}.out -p ${pre0}.prm -c ${pre0}.rst -ref ${pre0}.rst -r ${pre0}_min_v0_l${X}.rst
-O -i mdin_min_v1_l${X} -o ${pre1}_min_v1_l${X}.out -p ${pre1}.prm -c ${pre1}.rst -ref ${pre1}.rst -r ${pre1}_min_v1_l${X}.rst
EOF

cat << EOF > group_heat_l${X}
-O -i mdin_heat_v0_l${X} -o ${pre0}_heat_v0_l${X}.out -p ${pre0}.prm -c ${pre0}_min_v0_l${X}.rst -r ${pre0}_heat_v0_l${X}.rst
-O -i mdin_heat_v1_l${X} -o ${pre1}_heat_v1_l${X}.out -p ${pre1}.prm -c ${pre1}_min_v1_l${X}.rst -r ${pre1}_heat_v1_l${X}.rst
EOF

cat << EOF > group_equi_l${X}
-O -i mdin_equi_v0_l${X} -o ${pre0}_equi_v0_l${X}.out -p ${pre0}.prm -c ${pre0}_heat_v0_l${X}.rst -r ${pre0}_equi_v0_l${X}.rst -x ${pre0}_equi_v0_l${X}.crd
-O -i mdin_equi_v1_l${X} -o ${pre1}_equi_v1_l${X}.out -p ${pre1}.prm -c ${pre1}_heat_v1_l${X}.rst -r ${pre1}_equi_v1_l${X}.rst  -x ${pre1}_equi_v1_l${X}.crd
EOF
cat << EOF > group_prod_l${X}
-O -i mdin_prod_v0_l${X} -o ${pre0}_prod_v0_l${X}.out -p ${pre0}.prm -c ${pre0}_equi_v0_l${X}.rst -r ${pre0}_prod_v0_l${X}.rst -x ${pre0}_prod_v0_l${X}.crd -inf mdinf0.${X}
-O -i mdin_prod_v1_l${X} -o ${pre1}_prod_v1_l${X}.out -p ${pre1}.prm -c ${pre1}_equi_v1_l${X}.rst -r ${pre1}_prod_v1_l${X}.rst -x ${pre1}_prod_v1_l${X}.crd -inf mdinf1.${X}
EOF

cat << EOF > mdin_min_v0_l${X}
density minlibration
 &cntrl
  irest=0,ntx=1,
  imin=1,maxcyc=5000,drms=0.0001,ntmin=2,
  ntc=1,ntf=1,
  cut=8.0,
  ntpr=100,ntwx=0,ntwv=0,ntwe=0,
  clambda =0.${X},
  icfe=1,
EOF

cp mdin_min_v0_l${X} mdin_min_v1_l${X}
cat <<EOF >> mdin_min_v0_l${X}
 / 
EOF
cat <<EOF >> mdin_min_v1_l${X}
  ifqnt = 1,
 /
  &qmmm
  qmmask=':L03@2287-2301',
  qmcharge=0,
  qm_theory='DFTB3',
 / 
EOF
cat <<EOF > mdin_heat_v0_l${X}
Heat
 &cntrl
  irest=0,ntx=1,
  nstlim=10000,dt=0.001,
  ntb=1,temp0=300.0,ntt=3,gamma_ln=2.0,
  ntc=1,ntf=1,
  cut=8.0,
  ntpr=5000,ntwx=0,ntwv=0,ntwe=0,
  icfe=1,clambda=0.${X},
EOF
cp  mdin_heat_v0_l${X} mdin_heat_v1_l${X}
cat <<EOF >> mdin_heat_v0_l${X}
 /
EOF
cat <<EOF >> mdin_heat_v1_l${X}
  ifqnt=1,
 /
 &qmmm
  qmmask=':L03@2287-2301',
  qmcharge=0,
  qm_theory='DFTB3',
 / 
EOF

cat << EOF > mdin_equi_v0_l${X}
density equilibration
 &cntrl
  ntx=5,irest=1,
  nstlim=250000,dt=0.001,
  temp0=300.0,ntt=3,gamma_ln=2.0,
  ntc=1,ntf=1,
  cut=8.0,
  ntb=2,ntp=1,pres0=1.0,taup=1.0,
  ntpr=5000,ntwx=0,ntwv=0,ntwe=0,
  icfe=1,clambda=0.${X},
EOF
cp  mdin_equi_v0_l${X} mdin_equi_v1_l${X}
cat <<EOF >> mdin_equi_v0_l${X}
 /
EOF
cat <<EOF >> mdin_equi_v1_l${X}
  ifqnt=1,
 /
 &qmmm
  qmmask=':L03@2287-2301',
  qmcharge=0,
  qm_theory='DFTB3',
 / 
EOF
cat << EOF > mdin_prod_v0_l${X}
NPT productin
 &cntrl
  irest=1,ntx=5,
  nstlim=1000000,dt=0.001,
  temp0=300.0,ntt=3,gamma_ln=2.0,
  ntc=1,ntf=1,
  cut=8.0,
  ntb=2,ntp=1,pres0=1.0,taup=1.0,
  ntpr=1000,ntwx=1000,ntwv=0,ntwe=0,
  icfe=1,clambda=0.${X},ioutfm=0,
  ifmbar=1, bar_intervall=1000,bar_l_min=0.000,bar_l_max=1.000,bar_l_incr=0.1,
EOF
cp mdin_prod_v0_l${X}  mdin_prod_v1_l${X}
cat <<EOF >> mdin_prod_v0_l${X}
 /
EOF
cat <<EOF >> mdin_prod_v1_l${X}
  ifqnt=1,
 /
  &qmmm
  qmmask=':L03@2287-2301',
  qmcharge=0,
  qm_theory='DFTB3',
 / 
EOF

cat << EOF > minrun.sh.${X}
#!/bin/bash
#SBATCH -N 1
#SBATCH --job-name=3SM1-${X}
#SBATCH -p tetralith
#SBATCH --tasks-per-node=2
##SBATCH --exclusive
#SBATCH -A snic2022-3-19
#SBATCH -t 167:30:00


module purge
module load buildtool-easybuild/4.2.2-nsc2035f5dd0
module --ignore-cache load "binutils/.2.28"
module  load icc/.2018.1.163-GCC-6.4.0-2.28
module load  impi/.2018.1.163
module load SciPy-bundle/2020.03-Python-3.8.6-nsc1
module load GCC/11.3.0
module load  OpenMPI/4.1.4
module load SciPy-bundle/2022.05-nsc1
module load matplotlib/3.5.2-nsc1
module load Amber/16-nsc1-intel-2017.u7-bare

#mpirun -bind-to core -np 2  sander.MPI -ng 2 -groupfile group_min_l${X}
#mpirun -bind-to core -np 8 sander.MPI -ng 2 -groupfile group_heat_l${X}
#mpirun -bind-to core -np 8 sander.MPI -ng 2 -groupfile group_equi_l${X}
mpprun  -np 2 sander.MPI -ng 2 -groupfile group_prod_l${X}

EOF
sbatch  minrun.sh.${X}
done