Dear all,
I have a question about how the way restraints are handled in TI runs. To be more specific, I am running a simulation on a ligand-protein complex in which the ligand (residue 301) is decoupled in two steps (step 1: charge decoupling; step 2: vdw decoupling). The ligand is restrained with six restraints (1 distance, 2 angles and 3 torsions) and my particular question only pertains to the second step (vdw decoupling).
This is the corresponding pmemd.MPI input file of the second step (vdw decoupling):
===============================================================================
TI - VDW DECOUPLING
&cntrl
! simulation
imin = 0, ! flag to run minization (0=false)
nstlim = 5000000, ! number of MD steps (10.0 ns)
irest = 0, ! flag to restart simulation
t = 0 ! set time to 0
ntx = 1, ! option to read from inpcrd file
ntxo = 2, ! write rst as NetCDF file
dt = 0.002, ! timestep (2 fs)
ntb = 2, ! specify PBC type
nmropt = 1, ! use restraints
! thermostate
ntt = 3, ! switch for time scaling (3 = Langevin)
tempi = 300, ! initial temperature
temp0 = 300, ! reference temperature (if ntt > 0)
gamma_ln = 2.0, ! collision frequency for ntt=3
ig = -1, ! seed for Random Number Generator (RNG) (use -1 with ntt=2,3)
! barostate
ntp = 1, ! constant pressure dynamics (1=do scaling)
barostat = 2, ! barostat to use (2 = Monte Carlo)
mcbarint = 100, ! volume change attempt every n steps
pres0 = 1.00000, ! reference pressure (in bars)
taup = 2.0, ! pressure relation time (in ps)
! output
ioutfm = 1, ! file format (1=binary NetCDF)
iwrap = 0, ! wrap coordination into primary box (1=do)
ntpr = 1000, ! write energies to .mdout/.mdinfo every n steps
ntwr = 0, ! write rstrt file every n steps
! alchemical free energy
cut = 8.0
ntc = 2, ! SHAKE bond length constraints (2=H-atom)
ntf = 1, ! force evaluation (1=evalute all)
icfe = 1, ! turn on free energy calculation
ifsc = 1, ! use softcore parameters (0=undo)
logdvdl = 0, ! print summary of dvdl at end
clambda = 0.40, ! current lambda
timask1 = ':301', ! specify atoms unique to V_0 (charges uncoupled)
scmask1 = ':301', ! softcore atoms in V_0
timask2 = '', ! specify atoms unique to V_1 (NULL)
scmask2 = '', ! softcore atoms in V_1
crgmask = ':301', ! remove charges from ligand
ifmbar = 1, ! FEP
mbar_states = 21, ! Number of lambdas
mbar_lambda = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00,
/
&wt type='REST', istep1=0, istep2=0, value1=1, value2=1, /
&wt type='END' /
DISANG=restraint.rst
===============================================================================
As can be seen from the input file, the charges of the ligand (":301") are set to zero ("crgmask = ':301'"), and the restraints between ligand and protein are contained in a file named "restraint.rst". By setting "timask2 = ''", I am assuming that the VDW and ELECTROSTATIC parameters of the ligand in state 1 (when going from 0 to 1) are set to zero.
My question is now: will the restraints remain included when going from state 0 to state 1? Or are these restraints also removed in state 1? I am asking because from the output it looks that the NMR restraints are zero:
===============================================================================
| TI region 1
NSTEP = 66000 TIME(PS) = 132.000 TEMP(K) = 303.42 PRESS = 0.0
Etot = -73908.8593 EKtot = 18927.5909 EPtot = -92836.4502
BOND = 893.6617 ANGLE = 2485.6343 DIHED = 1581.5798
1-4 NB = 1061.1492 1-4 EEL = 11228.2615 VDWAALS = 9160.1205
EELEC = -119246.8572 EHBOND = 0.0000 RESTRAINT = 0.0000
EKCMT = 0.0000 VIRIAL = 0.0000 VOLUME = 298723.6027
Density = 1.0403
DV/DL = -49.3487
Ewald error estimate: 0.5740E-04
------------------------------------------------------------------------------
Softcore part of the system: 46 atoms, TEMP(K) = 334.22
SC_Etot= 112.6171 SC_EKtot= 39.1857 SC_EPtot = 73.4314
SC_BOND= 9.1324 SC_ANGLE= 42.4318 SC_DIHED = 22.4291
SC_14NB= 5.9907 SC_14EEL= 0.0000 SC_VDW = -6.5526
SC_EEL = 0.0000
SC_RES_DIST= 29.2903 SC_RES_ANG= 0.7853 SC_RES_TORS= 0.2738
SC_EEL_DER= -0.0000 SC_VDW_DER= -92.7074 SC_DERIV = -92.7074
------------------------------------------------------------------------------
| TI region 2
NSTEP = 66000 TIME(PS) = 132.000 TEMP(K) = 303.36 PRESS = 0.0
Etot = -73948.0450 EKtot = 18888.4052 EPtot = -92836.4502
BOND = 893.6617 ANGLE = 2485.6343 DIHED = 1581.5798
1-4 NB = 1061.1492 1-4 EEL = 11228.2615 VDWAALS = 9160.1205
EELEC = -119246.8572 EHBOND = 0.0000 RESTRAINT = 0.0000
EKCMT = 0.0000 VIRIAL = 0.0000 VOLUME = 298723.6027
Density = 1.0381
DV/DL = -49.3487
Ewald error estimate: 0.5740E-04
------------------------------------------------------------------------------
NMR restraints: Bond = 0.000 Angle = 0.000 Torsion = 0.000
===============================================================================
This seems confusing because I would expect that these Bond, Angle and Torsion restraints would be different from 0.
Can anyone provide me with any clues why these restraint values are zero?
Kind regards,
Hans
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu May 20 2021 - 04:00:02 PDT
