Self-Guided Langevin Dynamics (SGLD)
and
Self-Guided Molecular Dynamics (SGMD) Simulation Methods
By Xiongwu Wu and Bernard R. Brooks
SGMD/SGLD is a method to enhance conformational searching and sampling
through enhancing the low frequency motion. The main references for the
SGLD simulation method are:
(1) Wu, X., Brooks, B.R., "Self-guided Langevin dynamics simulation method", Chem. Phys. Letter, 381, 512-518(2003).
(2) Wu, X., Brooks, B.R., "Toward Canonical Ensemble Distribution from Self-guided Langevin dynamics simulation", J. Chem. Phys., 134, 134108 (2011).
(3) Wu, X., Brooks, B.R., "Force-momentum based self-guided Langevin dynamics: a rapid sampling method that approaches the canonical ensemble", J. Chem. Phys., 135, 204101 (2011)
(4) Wu, X., Damjanovic, A., Brooks, B.R., "Efficient and unbiased sampling of biomolecular systems in the canonical ensemble: a review of self-guided Langevin dynamics", Adv. Chem. Phys., Vol.150, 255-326 (2012)
(5) Wu, X., Hodoscek, M., Brooks, B.R., "Replica exchange of self-guided Langevin dynamics simulation", J. Chem. Phys., 137, 044106 (2012).
(6) Wu, X., Subramaniam,S., Case, D.A.,Wu,K., Brooks, B.R., "Targeted Conformational Search with Map-Restrained Self-Guided Langevin Dynamics: Application to Flexible Fitting of Electron microscopic density maps", J. Struct. Biol., 183, 429-440 (2013).
(7) Wu, X., Brooks, B.R., Vanden-Ejinden,E. "Self-Guided Langevin Dynamics via Generalized Langevin Equation", J. Comput. Chem., Accepted (2015).
* Menu:
* Syntax:: Syntax of the SGLD dynamics command
* Background:: Description of SGMD/SGLD methods
* Examples:: SGLD usage examples
[Syntax SGLD]
SCALar SGWT SET 0.0 ! clean guiding weight array: SGWT
SCALar SGWT SET 1.0 atom-selection ! select atoms to set guiding weights
! Using leap-frog integrator
DYNAmics {LEAP {[LANG SGLD|SGLDG]} } ! SGLD simulation
{ {[HOOVER|CPT] SGMD|SGMDG]} } ! SGMD for a NVT/NPT ensemble
{ {[ SGMD|SGMDG]} sg-spec} ! SGMD for a NVE ensemble
{ } other-dynamic-spec
! Using velocity-Verlet integrator
DYNAmics { VV2 {[SGLD|SGLDG]} } ! SGLD
{ {[SGMD|SGMDG]} sg-spec } ! SGMD
{ } other-specs
sg-spec::= {[SGFT real][SGBZ]}[TSGAVG real][TREFLF real] other-sg-spec
{[TEMPSG real] }
other-sg-spec::= [SGFF real] [TSGAVP real] [ISGSTA int] [ISGEND int]
[SGSTOPT] [SGSTOPR] [SGCOM]
Keyword Default Purpose
SGLD false Turn on SGLD simulation (for Langevin dynamics, LANG).
Friction constants, FBETA, must be set at first.
It is equivalent to SGMD for non-Langevin dynamics.
SGMD false Turn on SGMD simulation (for MD simulations).
It is equivalent to SGLD for Langevin dynamics.
SGLDG false Turn on SGLDG or SGMDG simulation to maintain ensemble distribution.
SGMDG false Turn on SGLDG or SGMDG simulation to maintain ensemble distribution.
SGBZ false Activate SGMDfp/SGLDfp to maintain a Boltzmann conformational
distribution. When SGBZ is on, SGFF is estimated from
simulation to mantain a Boltzmann distribution.
SGFT 0.2(SGMD) Momentum guiding factor. This is the lamda in the SGLD
1.0(SGLD) equation of motion. When TEMPSG=0, SGFT will be
constant throughout a simulation.
SGFF 0.0 Force guiding factor. Defines how much low frequency forces are included. For SGMDG, SGFT and SGFF are interdependent and only one of them need be set. For SGLDG, SGFF provides force based guiding force in addition to the momentum based guiding force set by SGFT.
TEMPSG 0 K Target self-guiding temperature. Please note that the TEMPSG
definition has been changed since c36. TEMPSG set the
conformatonal search ability comparable to that in a simulation
at temperature TEMPSG. TEMPSG>0 will override SGFT. a TEMPSG
higher than TBATH or TREF will enhance a conformational search
while a TEMPSG lower than TBATH or TREF will slow down a
conformational search. TEMPSG should not be set when SGBZ is set.
When TEMPSG>0, SGFT fluctuates during the simulation to maintain a guiding
temperature of TEMPSG. If a constant SGFT simulation is
desired, one can read the approximate values of SGFT from
the output of a SGMD/SGLD simulation with TEMPSG=2*TREF as a suggested value.
TSGAVG 0.2 ps Local average time. A larger TSGAVG will result in slower
motion to be enhanced. All motions with periods larger than
TSGAVG will be enhanced.
ISGSTA 1 The index of the first atom of the region to apply the
guiding force
ISGEND natom The index of the last atom of the region to apply the
guiding force. If reweighting is to be performed and ISGSTA>1 or
ISGEND<natom, one need turn on analysis with "anal on" before
running dynamics.
SGSTOPT false Stop translation of the center of mass. By default, SGLD
allows the center of mass to have an random motion. By set
SGSTOPT, SGLD stops any translation in the center of mass.
SGSTOPR false Stop rotation around the center of mass. By default, SGLD
allows the center of mass to have an random motion. By set
SGSTOPR, SGLD stops any rotation around the center of mass.
In SGLD simulations, NTRFRQ has no effect if SGSTOPT and
SGSTOPR is set.
SGCOM false Allow net guiding force and torque to the system. By default the net
guiding force and torque on the center of mass is removed. Set SGCOM
allows accelerating the motion of the center of mass.
TSGAVP 10*TSGAVG Convergency time. This time is used to control convergency of
running averages to estimate instaneous parameters like SGFTI,
TEMPSGI, TEMPLFI, etc. A larger TSGAVP will result in a smaller
fluctuation of instaneous parameters.
TREFLF 0 K Reference low frequency temperature in a canonical ensemble.
TREFLF is the low frequency temperature when no guiding force is
applied to a simulation system. TREFLF is required for weighting
factor calculation in SGMD/SGLD simulation or for the guiding
force calculation in SGMDfp/SGLDfp simulations (when SGBZ is
turned on). TREFLF can be determined from a SGMD/SGLD simulation
with SGFT=0 and TEMPSG=0, TREFLF=<TEMPLF>. The average TEMPLF value
can be accessed by querying ?TEMPLF in CHARMM input script.
When TREFLF=0, TREFLF is estimated during a simulation. An actual
value for TREFLF can increase the accuracy in SGMDfp/SGLDfp simulations.
SGWT {1.0} Self-guiding weight array for atoms. Using SCALAR to set SGWT.
SGWT is used to select atoms to apply guiding forces.
Backgrounds
SGMD/SGLD simulation enhances conformational search efficiency through
acceleration of low frequency motions in a molecular system. The low frequency
motion represents the slow conformational change that contributes the most
in conformational search. A guiding force resonant with the low freuqency
motion, which is calculated based on a local average of momentums and/or forces, is introduced
into the equation of motion to enhance the low frequency motion. It has been
demonstrated that within the suggested parameter range, a SGLD simulation
has a dramatically enhanced conformational search efficiency while has little
perturbation on conformational distribution.
Recently, based on the generalized Langevin equation, we developed the SGLD-GLE method that restores the detailed-balance property of the dynamics. This
permits to sample the NVT and NPT ensembles exactly, without the reweighting needed in SGLD. Compared with other sampling enhancement methods, SGLD-GLE has several unique characters. First, SGLD-GLE is size
extensive. The guiding force is calculated from momentum
and random forces, independent of system size. Reweighting is no longer needed which allows the method to be cast in a manner that is size extensive, making it a possible method of choice for million atom systems that may benefit from enhanced sampling. Second, the method can be applied to a part of a simulation system. In other words, the guiding factor can be defined for individual atoms, just like the friction con-
stant can be atom specific. Third, the motion mode to be
enhanced can be controlled by the local averaging time.
Because of these features, we expect that SGLD-GLE, like
SGLD, will be useful in a wide variety of contexts, such as phase transition, ligand docking, and protein folding.
A further development lead to the generalized SGLD (SGLDg), which includes the force based guiding effect and can be applied to MD as SGMDg. SGLDg/SGMDg have enhanced conformational sampling while maintain correct ensemble distribution.
The self-guided dynamics simulation enhances conformational search and sampling
efficiency by accelerating low frequency motion. The conformational distribution in
a self-guided dynamics simulation can be quantitatively described by the so-called
low frequency properties, like low frequency temperature and low frequency energy.
A weighting factor for each frame is printed out during simulation and ensemble average
properties can be calculated as weighted avearges. The self-guided dynamcs simulation
methods are further improved to directly produce a canonical ensemble distribution
by incooperate the low frequency force and low frequency motions in the guiding
force calculation. This new method is referred to as SGMDfp or SGLDfp, depending on
the presents of the friction forces. This extension makes SGMD/SGLD an excellent tool
for both conformational search and conformational sampling.
In CHARMM input script, only a key word, SGMD/SGLD, or SGMDG/SGLDG in a "DYNA ..."
simulation command line is needed to turn on SGMD/SGLD or SGMDg/SGLDg with default values
for all parameters. For LD or SGLD simulation, the friction constant, FBETA,
need be set. Whether to use SGMD or SGLD is automatically decided depending on
whether the friction constant is zero or not.
SGMD/SGLD has two primary input parameters, the guiding factor, SGFT, and
the local average time, TSGAVG, which define the guiding effect in a SGMD/SGLD
simulation. SGFT set the strength of the guiding force and TSGAVG defines the
slow motion mode to be enhanced. A larger TSGAVG will result in a slower
motion (with period longer than TSGAVG) to be enhanced. A larger SGFT will
introduce stronger guiding forces and result in a larger energy barrier
overcoming ability. When SGFT=0, a SGMD/SGLD simulation reduces to a normal
MD/LD simulation. SGFT should be limited to keep low frequency motion from changing
too much. The low frequency motion is measured by the low frequency temperature,
which is calculated from the local average momentum and depends on TSGAVG and SGFT.
A scalar array SGWT can be used to define atomic specific guiding factors. the actual guiding factor of atom i will be SGFT*SGWT(I).
To avoid the difficulty in choosing the value of SGFT, a parameter
called guiding temperature (TEMPSG) is introduced. TEMPSG, as an alternate
to SGFT, specifies the conformational searching ability that is comparable to
a simulation at temperature TEMPSG. When TEMPSG=0, SGFT will be used to calculate
the guiding force. If TEMPSG>0, SGFT will fluctuate during a simulation to
reach a conformational search comparable to a simulation at a temperature of TEMPSG.
TEMPSG>TBATH for a LD simulation or TEMPSG>TREF for a MD simultaion will enhance
conformational search, while TEMPSG<TBATH or TEMPSG<TREF will slow down conformational
search.
The Key word, SGBZ, will turn on the SGMDfp/SGLDfp method to maintain a
canonical ensemble. A input parameter, TREFLF, can be used to input
the reference low frequency temperature, which is the low frequency temperature
when no guiding force is applied. If TREFLF is not input, it will be estimated
during a SGMD/SGLD simulation. When SGBZ is set, TEMPSG should not be used to set
guiding effect, instead, SGFT should be used.
The guiding range, ISGSTA and ISGEND, define the starting and ending atom indexes
to apply the guiding forces. By default, all atoms are included in the guiding range.
Within the range, the guiding weight array, SGWT, can be set with the SCALAR command.
The actual guiding factor is SGFT*SGWT(I) for atom I. By default, SGWT=1.
SGLD can be used for replica exchange, abbreviated as RXSGLD, to obtain
accelerated conformational search while preserves the canonical ensemble. The guiding
effect, SGFT or TEMPSG, is used to define different replicas, with or without
temperature changes. Without temperature changes, RXSGLD can achieve very high
exchange efficiency and can be applied to much larger systems than temperature-based
replica exchange simulations. More details can be found in repdstr.doc.
When SGMD or SGLD is turned on, DYNA simulation output will contain
two lines with the following quantities:
DYNA SGLF> SGFT TEMPSG TEMPLF TREFLF FRCLF EPOTLF SGWT
DYNA SGHF> SGFF SGFD TEMPHF TREFHF FRCHF EPOTHF VIRSG
These quantities are instaneouse values defined as below:
SGFT: Momentum guiding factor
SGFF: Force guiding factor. Adjusted when SGBZ is on
SGFD: Force dumping factor. Adjusted when SGBZ is on.
TEMPSG: Guiding temperature.
SGEXP: Probability exponent. exp(SGEXP) is the weighting factor of current frame.
VIRSG: Virial of the guiding force.
TEMPLF: low frequency temperature
TEMPHF: high frequency temperature
TREFLF: reference low frequency temperature
TREFHF: reference high frequency temperature
FRCLF: low frequency force factor
FRCHF: high frequency force factor
EPOTLF: low frequency potential energy
EPOTHF: high frequency potential energy
The weight of a conformation is calculated by
Weight= exp(SGEXP)
=exp(((FRCLF*TREFLF/TEMPLF-1)*EPOTLF+(FRCHF*TREFHF/TEMPHF-1)*EPOTHF+
VIRSG)/(KBOLTZ*Temp))
The averages of above variables can be queried by the following
parameters:
?SGFT=<SGFTI>, ?SGFF=<SGFFI>, ?TEMPSG=<TEMPSG>, ?SGFD=<SGFD>,?SGEXP=<SGEXP>,
?TEMPLF=<TEMPLF>,?TEMPHF=<TEMPHF>, ?TREFLF=<TREFLF>,?TREFHF=<TREFHF>,
?FRCLF=<FRCLF>,?FRCHF=<FRCHF>,?EPOTLF=<EPOTLF>,?EPOTHF=<EPOTHF>,?VIRSG=<VIRSG>,
In addition to the DYNA output, for each trajectory frame, there is a
SGMD/SGLD output marked by "TRAJ SGLF>" and "TRAJ SGHF>", which can be used for
post-processing ensemble distributions.
SGMD/SGLD simulation results can be weighted to produce canonical ensemble
properties. The weighting information can be extracted from simulation outputs.
A simple command line awk script is provided here to get a weighted average
from a SGMD/SGLD simulation:
%> awk 'BEGIN{n=0;if(N0==0)N0=1000;} \
{if($1=="DYNA>"){ep=$6;n++;} if(n<N0)next; \
if($1=="DYNA"&& $2=="SGLF>"){wt=exp($9);aw+=wt;ae+=ep;aew+=wt*ep;m++;}} \
END{printf("Ep and weighted Ep: %10.4f %10.4f\n",ae/m,aew/aw);}' sgmd-sgld.out
The following awk script will print weighting factors for every frame in
a SGMD/SGLD trajectory:
%> awk '{if($1=="TRAJ"&&$2=="SGLF>")printf(" %10d %10.4f\n",++i,exp($9));}' sgmd-sgld.out
The following awk script will recalculate and print weighting factors for every frame in
a SGMD/SGLD trajectory:
%> cat >sgwt.awk
BEGIN{n=0;if(TEMP==0)TEMP=300;kt=0.001987*TEMP;if(N0==0)N0=10;} \
{if($1=="TRAJ"&&$2=="SGLF>"){templf=$5;treflf=$6;frclf=$7;epotlf=$8;sgwt=$9;} \
if($1=="TRAJ"&&$2=="SGHF>"){temphf=$5;trefhf=$6;frchf=$7;epothf=$8;virsg=$9; \
eplfs[n]=epotlf;ephfs[n]=epothf;virs[n]=virsg;n++;if(n<N0)next; \
avgtlf+=templf;avgrlf+=treflf;avgflf+=frclf;avgelf+=epotlf; \
avgthf+=temphf;avgrhf+=trefhf;avgfhf+=frchf;avgehf+=epothf;avgvir+=virsg;}} \
END{m=n-N0;avgtlf/=m;avgrlf/=m;avgflf/=m;avgelf/=m; \
avgthf/=m;avgrhf/=m;avgfhf/=m;avgehf/=m;avgvir/=m; \
if(TREFLF>0){avgrlf=TREFLF;avgrhf=TEMP-TREFLF;} if(TREFHF==0){avgthf=TEMP-avgtlf;} \
for(i=0;i<n;i++){ \
wt=exp(((avgflf*avgrlf/avgtlf-1)*(eplfs[i]-avgelf)+(avgfhf*avgrhf/avgthf-1)*(ephfs[i]-avgehf)+virs[i]-avgvir)/kt); \
printf(" %6d %10.4f\n", i+1,wt);}}
EOF
Run awk with above script can produce weighting factor for each frame in the trajectory file:
%> awk -f sgwt.awk -v N0=100 -v TREFLF=8.5 -v TEMP=400 sgmd-sgld.out
Here, N0 is equilibrium frame number (default 10), TEMP is simulation temperature (default is 300K),
and TREFLF is the reference low frequency temperature (default is the estimated value from simulation).
Examples
1) SGMD simulation with default setting.
! Perform SGMD simulation
DYNA LEAP STRT CPT NSTE 1000000 TIME 0.002 -
IPRFRQ 10000 ISVFRQ 1000 IHTFRQ 0 IEQFRQ 0 INBFRQ 10 IHBFRQ 0 -
IUNREA -1 IUNWRI 31 IUNCRD 32 IUNVEL 33 KUNIT -1 -
NSAVC 1000 NSAVV 00 NPRINT 1000 ISEED 314159 -
SGMD -
TCON TCOU 0.1 TREF 300 FIRST 260 -
NBXMOD 5 ATOM CDIEL VATOM VDISTANCE EIPS VIPS PXYZ -
CUTNB 12.0 CTOFNB 10.0 EPS 1.0 E14FAC 1.0 WMIN 1.0
2) SGLD simulation with default setting.
! Set friction forces
SCAL FBETA SET 1.0 SELE ALL END
! Perform SGLD simulation
DYNA LANG LEAP STRT NSTE 1000000 TIME 0.002 -
IPRFRQ 10000 ISVFRQ 1000 IHTFRQ 0 IEQFRQ 0 INBFRQ 10 IHBFRQ 0 -
IUNREA -1 IUNWRI 31 IUNCRD 32 IUNVEL 33 KUNIT -1 -
NSAVC 1000 NSAVV 00 NPRINT 1000 ISEED 314159 -
SGLD -
TBATH 300 FIRST 260 -
NBXMOD 5 ATOM CDIEL VATOM VDISTANCE EIPS VIPS PXYZ -
CUTNB 12.0 CTOFNB 10.0 EPS 1.0 E14FAC 1.0 WMIN 1.0
3) SGLD simulation with a constant guiding factor.
! Set friction forces
SCAL FBETA SET 1.0 SELE ALL END
! Perform SGLD simulation
DYNA LANG LEAP STRT NSTE 1000000 TIME 0.002 -
IPRFRQ 10000 ISVFRQ 1000 IHTFRQ 0 IEQFRQ 0 INBFRQ 10 IHBFRQ 0 -
IUNREA -1 IUNWRI 31 IUNCRD 32 IUNVEL 33 KUNIT -1 -
NSAVC 1000 NSAVV 00 NPRINT 1000 ISEED 314159 -
SGLD TSGAVG 0.2 SGFT 1.0 -
TBATH 300 FIRST 260 -
NBXMOD 5 ATOM CDIEL VATOM VDISTANCE EIPS VIPS PXYZ -
CUTNB 12.0 CTOFNB 10.0 EPS 1.0 E14FAC 1.0 WMIN 1.0
4) SGLD simulation with a guiding temperature of 500 K (to reach a searching ability comparable to
a 500 K high temperature simulation).
! Set friction forces
SCAL FBETA SET 1.0 SELE ALL END
! clean self-guiding weights.
SCAL SGWT SET 0.0 SELE ALL END
! apply guiding force to only segment 1
SCAL SGWT SET 1.0 SELE iseg 1 end
! Perform SGLD simulation
DYNA LANG LEAP STRT NSTE 1000000 TIME 0.002 -
IPRFRQ 10000 ISVFRQ 1000 IHTFRQ 0 IEQFRQ 0 INBFRQ 10 IHBFRQ 0 -
IUNREA -1 IUNWRI 31 IUNCRD 32 IUNVEL 33 KUNIT -1 -
NSAVC 1000 NSAVV 00 NPRINT 1000 ISEED 314159 -
SGLD TSGAVG 0.2 TEMPSG 500 -
TBATH 300 FIRST 260 -
NBXMOD 5 ATOM CDIEL VATOM VDISTANCE EIPS VIPS PXYZ -
CUTNB 12.0 CTOFNB 10.0 EPS 1.0 E14FAC 1.0 WMIN 1.0
5) SGLDfp simulation with a constant guiding factor and to maintain ensemble distribution.
! Set friction forces
SCAL FBETA SET 1.0 SELE ALL END
DYNA LANG LEAP STRT NSTE 1000000 TIME 0.002 -
IPRFRQ 10000 ISVFRQ 1000 IHTFRQ 0 IEQFRQ 0 INBFRQ 10 IHBFRQ 0 -
IUNREA -1 IUNWRI 31 IUNCRD 32 IUNVEL 33 KUNIT -1 -
NSAVC 1000 NSAVV 00 NPRINT 1000 ISEED 314159 -
SGLD SGBZ TSGAVG 0.2 SGFT 1.0 -
TBATH 300 FIRST 260 -
NBXMOD 5 ATOM CDIEL VATOM VDISTANCE EIPS VIPS PXYZ -
CUTNB 12.0 CTOFNB 10.0 EPS 1.0 E14FAC 1.0 WMIN 1.0
6> SGLDg simulation with VV2 integrator to maintain ensemble distribution
DYNA VV2 START NSTEP 10000 TIME 0.001 -
SGLDG SGFT 1.0
CHARMM Documentation / Rick_Venable@nih.gov