PLUMED User's Guide
Contents
1. DAFED CV 1 TEMPERATURE 600 MASS 1e5 KAPPA 1e6 TAUTHERMO 0 2 JACOBIAN_FORCE DAFED CV 2 TEMPERATURE 600 MASS 1e3 KAPPA 1e4 TAUTHERMO 0 2 PERIODIC MINUS_PI PLUS_PI DAFED_CONTROL RESTART checkpoint_file WRITE_STATE 1 N_RESPA 1 PRINT W_STRIDE 100 ENDMETA A separate DAFED CONTROL directive contains general controls for the dAFED simulation The dAFED dynamics including all variables described in Eqs 1 3 1 6 can be restarted exactly from a previous run using a check point file Following the WRITE_STATE keyword appears the number of steps after which a checkpoint file is saved A value of 1 implies that a check point file is written only when GROMACS saves its own checkpoint file i e at regular wall clock time intervals The checkpoint file is saved in the current directory with default name DAFED_STATE The optional keyword RESTART is used to specify the path to the checkpoint file from which to restart The integrator for S can be selected on the DAFED_CONTROL line with the keyword INTEGRATOR followed by either GGMT default or LANGEVIN_EM for Langevin evolution with the simple Euler Maruyama integrator or LANGEVIN_CV for Langevin evolution with the integrator of Ciccotti Vanden Eijnden 6 Specifying LANGEVIN defaults to LANGEVIN_EM With a high value of x as required by the dAFED method oscillations of s r can become faster than the fastest mode in the physical system This would in principle require choosing a smaller time2. K is given after keyword TEMPERATURE The thermostat time constant 7 is given in ps after keyword TAUTHERMO The mass mg and harmonic constant K are given after the keywords MASS and KAPPA respectively The units of k and ms depend on the nature of the CV They should always be such that kS and mS are both in units of energy kJ mol amu nm ps see the example below In addition tow optional keywords can be used with the DAFED directive First for periodic CVs such as torsion angles the S variable should also evolve on a periodic interval This is specified by the keyword PERIODIC followed by two numbers for the lower and upper bounds The numbers can be replaced by MINUS_PI PLUS_PI or PLUS_2PI to specify a 7 or 27 respectively The optional keyword JACOBIAN FORCE causes a bias force F 2kpT S to be applied to the dynamics of S This is useful with distance CVs in order to counterbalance the effect of the Jacobian factor and sample a more uniform distribution along the CV Example The following lines couple CV 1 a distance in nm to a meta variable of mass 10 amu with a harmonic constant of 10 kJ mol nm and CV 2 a unitless number to a meta variable with mass 10 amu nm with a harmonic constant of 10 kJ mol For both CV the dAFED temperature is 600 K and the GGMT thermostat time constant is 0 2 ps See text for the optional Keywords JACOBIAN_FORCE and PERIODIC DISTANCE LIST 1 34 TORSION LIST 5 15 29 36
3. can be calculated 2 P Hs S ps n Pn pe nem V S Js kas 1 OLT Wee P s r 2 1 8 ms The effective adiabaticity of the coupling can thus be asserted In addi tion for each extended variable S the quantity Hs Ws should be strictly conserved which provides a quality check for the simulation In addition Let H be the pseudo energy of the physical system including the associ ated thermostats and barostats Then the total energy of the simulation H 4 Hg should be conserved as well As a corollary considering the physical system only the quantity H 2a Ws should be conserved A different way to assess the adiabaticity of the simulation is through the configurational temperature of the collective variable s _ 1 VU s _ s ons kp V2U s kg For a system at equilibrium the configurational temperature should be equiv alent to the kinetic temperature and the heat bath temperature A higher configurational temperature is a signature of nonequilibrium dynamics in which a significant heat flow would take place between extended and physi cal system In practice the choice of ms is subject to some pragmatic considerations The value of mg should be as high as possible to ensure good adiabatic sep aration However given the limited planned simulation time the evolution of S has to be fast enough to correctly sample the CV range of interest By running a short dAFED simulation and plotting the
4. evolution of S one can estimate an average diffusion speed From that a maximum admissible value for mg can be deduced such that S can cross many time the CV range during the simulation The choice of the coupling constant x determines the resolution of the observed free energy surface G S Ideally x should be very large but its value is limited by the requirement of integrating accurately the coupling term Eq The typical period of that harmonic oscillator is given by T 27 u k where u MgmMeg Ms Meg is the reduced mass with Meg Te 8 1 9 the effective mass of the CV s r For a one dimensional CV meg can be expressed as mer fy i l 1 10 Note that if the CV is multidimensional the situation is slightly more complicated and meg is in fact a tensor whose diagonal elements are given by Eq Nevertheless meg can be used to estimate an order of magnitude for the period 7 for each CV and thus estimate which time step is appropriate for a given msg or vice versa It is recommended to check that the time step is appropriate for the chosen by plotting the evolution of s r and S from a short simulation in which the dAFED variables are saved very frequently 1 1 1 Input for d AFED For each CV a DAFED directive is used to define the parameters of the cor responding dynamics On the same line the number of the CV to which the directive applies is specified after the keyword CV The temperature Ts in
5. example shows how to setup a UFED run CV 1 a distance in nm is coupled with a harmonic constant of 10 kJ mol nm to a meta variable of mass 10 a m u at temperature 400 K In this case we have specified a Langevin integrator for this meta variable and we print a compact COLVAR file For UFED we deposit every 1000 steps a hill of height 0 5 kJ mol and width 0 1 nm In addition we have restricted the space of the meta variable with a lower reflection wall at 0 1 nm and an upper reflection wall at 1 5 nm DISTANCE LIST 1 2 SIGMA 0 1 DAFED CV 1 TEMPERATURE 400 MASS 1e6 KAPPA 1e5 TAUTHERMO 0 5 DAFED_CONTROL WRITE_STATE 1 INTEGRATOR LANGEVIN PRINT_NO_DETAILS UFED_HILLS HEIGHT 0 5 W_STRIDE 1000 LREFLECT CV i LIMIT 0 1 UREFLECT CV 1 LIMIT 1 5 PRINT W_STRIDE 100 ENDMETA We note that in addition to G S it is possible to calculate the ensemble 8 average of any observable A r in during a d AFED or UFED simulation even if the distribution of states is different form the canonical ensemble at low physical temperature It can be shown that if the adiabatic separation is effective we have _ Jas Ajs e PAS A fds eE 1 15 where A s is the average value of A accumulated on a grid of fixed S posi tions Here G S is obtained from the same simulation by integrating values of fr s g accumulated on the same grid The integral in Eq 1 15 is per formed numerically a posteriori Using Eq dAFED can for example be com
6. step at the expense of sampling efficiency Instead following a multiple time step approach the dAFED force can be integrated more often than the forces in the physical system This feature is implemented only with GROMACS and the GGMT thermostat The user has to divide the general MD time step by a number Nrespa typically between 2 and 10 The optional keyword N_RESPA followed by the number Nergpa in the DAFED_CONTROL directive instructs GROMACS to evaluate the physical forces only every NRespa steps Together with the md vv integrator of GROMACS this should produce a correct RESPA 7 5 scheme in the NVT ensemble With this nstcalcenergy 1 has to be set in the GROMACS input file Note that this feature is still experimental and energy conservation should be checked 1 1 2 Typical output for dAFED With the dAFED method the COLVAR file will contain the following data if d collective variables are used e time step e value of the collective variable s r sa r Then for each of the S j 1 d appears a set of 5 columns with e the meta variable S e the instantaneous temperature of S in K e the conserved quantity Hs see Eq in kJ mol e the work Ws from S to the physical system see Eq in kJ mol e the effective mass meg according to Eq 1 10 in a m u These fields are labeled sj T sj E sj W sj and Meffj respectively in the COLVAR header line j 1 d Additional collective variables can be monitored du
7. PLUMED User s Guide Complement for the AFED and UFED methods Michel A Cuendet PLUMED 1 3 dAFED 0 7 July 2012 Chapter 1 Running free energy simulations 1 1 Driven Adiabatic Free Energy Dynamics dAFED The driven adiabatic free energy AAFED algorithm I is also called tem perature accelerated molecular dnamics TAMD by other authors 2 In fact AFED TAMD is an improvement over the earlier AFED method 3 4 which required cumbersome coordinates transformations In the AFED TAMD method an extra dynamical variable S is coupled to a collective variable s r where r represents the coordinates of a number N of atoms in the system The coupling is mediated by a potential energy function with harmonic con stant K V S 8 x i S s n 1 1 The dynamics of the S meta variable is adiabatically decoupled from the dynamics of the underlying physical system by choosing a large mass ms gt gt mM where m is a typical mass of the physical system Thanks to the adiabatic separation a temperature Ts gt T can be assigned to the S meta variable With this choice of mg and T the physical system will evolve fast at room temperature T around the instantaneous value of s r S On the other hand S will evolve slowly but have a temperature large enough to drive the system over high free energy barriers In the limit of k oo it can be shown that the free energy surface at temperature T can be recovered from the
8. UFED_HILLS works just as the HILLS directive of metadynamics It must be followed by a key word HEIGHT after which the value of h is specified see Eq 1 14 The hill deposition stride is specified after the keyword W_STRIDE The Gaussian widths g are taken from the keywords SIGMA specified on the line of each CV i 1 d Some collective variables have intrinsic domain limitations such as the number of H bonds that cannot be smaller than zero or the user might want to impose limitations such as the maximum distance to which a ligand can be separated from it s host In these cases it is useful to impose the limitations to the domain of the S variables which are otherwise unbounded This is especially necessary when a bias potential is used One way to do this without perturbing the distribution of S within the range of interest is to use reflective walls at which the momentum Ps is inverted Reflective walls are activated with the directives LREFLECT and UREFLECT corresponding to a lower or upper limit respectively The CV on which the reflective wall acts is specified after the keyword CV and the limit value is given after the keyword LIMIT If UFED_HILLS and LREFLECT or UREFLECT are active extra hills are added at a symmetrical position on the other side of the wall as soon as S is closer than 30 to the wall This prevents the formation of an artificial ditch in the bias potential close to the wall IT Example The following
9. bined with thermodynamic integration A dH d or free en ergy perturbation A exp AH A to calculate alchemical free energy differences in flexible molecules Bibliography 1 5 7 8 J B Abrams M E Tuckerman Efficient and direct generation of mul tidimensional free energy surfaces via adiabatic dynamics without coor dinate transformations J Phys Chem B 112 2008 15742 15757 L Maragliano E Vanden Eijnden A temperature accelerated method for sampling free energy and determining reaction pathways in rare events simulations Chem Phys Lett 426 2006 168 175 L Rosso M E Tuckerman An adiabatic molecular dynamics method for the calculation of free energy profiles Mol Simul 28 1 2002 91 112 L Rosso P Minary Z Zhu M E Tuckerman On the use of adiabatic molecular dynamics to calculate free energy profiles J Chem Phys 116 2002 4389 4402 Y Liu M E Tuckerman Generalized Gaussian moment thermostat ting A new continuous dynamical approach to the canonical ensemble J Chem Phys 112 4 2000 1685 1700 E Vanden Eijnden G Ciccotti Second order integrators for langevin equations with holonomic constraints Chem Phys Lett 429 1 3 2006 310 316 G J Martyna M E Tuckerman D J Tobias M L Klein Explicit reversible integrators for extended systems dynamics Mol Phys 87 5 1996 1117 M Chen M A Cuendet M E Tuckerman Heating and flood
10. density p S sampled at temperature Ts during the adiabatic dAFED simulation using G S kpTs log p S 1 2 This result generalizes well to the case where more than one collective variable is used and G S is a multi dimensional free energy surface The dAFED method requires very efficient thermostatting of the meta variable S In the present implementation S is coupled to a Generalized Gaussian Moment Thermostat GGMT 5 or to a Langevin thermostat 6 If multiple reaction coordinates are used one separate GGMT or Langevin thermostat is associated to each of them In the case of GGMT the meta variable is coupled to two thermostatting variables p and pc with associated masses Q and Qe respectively Given a typical time scale 7 of the thermostated system optimal masses are Q kpT st and Qg kgTs r The order 2 GGMT dynamics for one degree of freedom is 2 Pn Pe 1 De gg L pete See 1 ps V S s 0 Zeps PE ints 528 ps 1 3 p Dn kpTs 1 4 1 p w 3 Ea kpTs 1 5 Ps Pn n De Pn L K 1 6 The implemented integrator for the dynamics above is based on a Trotter decomposition of the corresponding Liouville operator I The quality of the integration can be monitored using the quantity Hs which would be conserved if the dynamics of was decoupled from the physical system Pi p 20 20 The heat transfer Ws from the meta variable to the physical system
11. ing A unified approach for rapid generation of free energy surfaces J Chem Phys 137 2012 024102 10 9 J Kastner W Thiel Bridging the gap between thermodynamic integra tion and umbrella sampling provides a novel analysis method umbrella integration J Chem Phys 123 2005 144104 10 L Maragliano E Vanden Eijnden Single sweep methods for free energy calculations J Chem Phys 128 2008 184110 11 Y Crespo F Marinelli F Pietrucci A Laio Metadynamics conver gence law in a multidimensional system Physical Review E 81 5 2010 55701 12 M A Cuendet M E Tuckerman Alchemical free energy differences in flexible molecules from thermodynamic integration or free energy per turbation combined with driven adiabatic dynamics J Chem Theory Comput 2012 in press doi 10 1021 ct300090z 11
12. ring a dAFED run in which case more columns will appear before the first set of dAFED fields For long production runs the user can choose a more compact form of output in which only the CV s r and the position S of each corresponding extended variable are printed This is achieved by adding the keyword PRINT_NO_DETAILS on the DAFED_CONTROL input line 1 1 3 Unified Free Energy Dynamics UFED UFED is a recent extension 8 of the AFED TAMD method which com bines high temperature extended variables with adaptive bias potentials sim ilar to those used in metadynamics The UFED method rests on the fact that the free energy surface can be reconstructed from the thermodynamic force F S instead of the histogram p gt S Note that the use of the force is also found in some other methods such as the umbrella integration method O or the single sweep method 10 Essentially in the spirit of the well known thermodynamic integration technique we can write FJ VsGp S oo gg fae se 9 Pese 1 12 _ ee 1 13 The last line represents the force exerted by the physical system on the extended variable averaged at a fixed position of S This average can easily be obtained in a post processing phase from the values of S and s r stored in the COLVAR file using a grid in the S space Note that due to the fast oscil lations of s r samples should be collected at high frequency Finally F S is integrated numerically to get the free energ
13. y profile In dimensions greater than one F S will not exactly be a consistent multidimensional gradient due to statistical noise The PMF can however easily be reconstructed as the surface G g S whose discrete derivative best fits F S in the least squares sense This postprocessing step provides the additional benefit of producing a smooth PMF The second fundamental ingredient for the UFED method is that if the adiabatic separation is effective F S does not depend on the actual distri bution of s 8 If this holds we can introduce a bias potential of any kind acting on S We introduce a Gaussian based adaptive potential a S S kr Vil Sty n 5 exp gt 552 kr lt t i 1 1 14 which is similar to the metadynamics bias potential except that it acts on the extended variables S instead of the CV s r UFED has several advantages over its parent methods First the adaptive bias allows using a lower temper ature than dAFED which facilitates obtaining effective adiabatic separation Second the use of the force to construct the free energy surface instead of the sum of hills in metadynamics makes the final accuracy of UFED inde pendent of the hill size and does not require that all basins are filled up with hills In order to activate UFED the user only needs to add a line with the directive UFED_HILLS in addition to the DAFED and DAFED_CONTROL directives described above see example below For most aspects
Download Pdf Manuals
Related Search