share/README_FreeEnergyModifications.txt

   1 ======== GROMACS 4.5 Free Energy Modifications ===================
   2
   3 New with this version:
   4
   5 - coulombic, vdw, bonded, and restraint transformation lambdas all decoupled
   6   from each other -- any specified pathway can be taken.
   7 - extended ensemble MC, with a number of different move choices
   8 - Supports either fixed weights, or dynamically adjusted weights
   9 - free energy dependent dihedral, angle, and distance restraints
  10
  11 === BASIC DOCUMENTATION ===
  12
  13 Important options that may have changed behavior in this version: --
  14 defaults listed first, other options in [ ], then parenthetical
  15 comments.
  16
  17 free-energy              = no [yes,mutate,decouple]
  18
  19 lmc-stats            = no [metropolis-transition,barker-transition,wang-landau,gibbs-wang-landau,minvar]  (the
  20 method used to update the weights)
  21
  22 nst-transition-matrix    = 0 [nonnegativeinteger] frequency at which the transition-matrix is output to the log
  23 lmc-mc-move          = gibbs [metropolis,barker,gibbs,metropolized-gibbs] (the method used to for lambda MC moves)
  24
  25 lmc-seed                 = -1 [for lambda mc transformations, -1 means it's taken from the process number]
  26 mc-temperature          = [positive real] [If omitted, set to the same as the ref_t for the 0th system, if there is one]
  27
  28 lmc-gibbsdelta          = -1 [any integer] (the interval [lambda-1,lambda+1] for moves for gibbs sampling on the lambdas.  -1,
  29 the default, means the entire interval [0,1] will be sampled)
  30
  31 initial-wl-delta        = 1.0 [any positive real] (the initial delta factor for Wang-Landau, in kT)
  32 wl-scale                = 0.8 [real between 0 and 1] (the scaling factor for Wang-Landau)
  33 wl-ratio                = 0.8 [real between 0 and 1] ratio of lowest to highest occupancies for Wang-Landau being seen as flat)
  34 symmetrized-transition-matrix=no [yes/no]  Whether to compute a symmetrized version of the transition matrix, by averaging the transpose with the matrix.
  35
  36 lmc-forced-nstart       = 0 [any positive integer] (the number of equilibration steps at each lambda when
  37                             'forcing through', essentially peforming slow growth, to initialize weights
  38
  39 mininum-var-min        = 100 [any positive integer] must have this many samples in each state before reweighting to
  40                          get the number of states per state that minimizes variace.  Preferable to make it
  41                          something decent so that it the initial weights don't get swallowed up in noise,
  42                          as it is using an asymptotic theory.
  43
  44 weight-c-range          = 0 when using minvar, uses the C that is closest to the G*(ln n0/n1).
  45                           0 means it defaults to C=0. Otherwise, it computes the free energies
  46                           along +/- c-range, with kT spacing.  Should eventually be used for
  47                           BAR as well.
  48
  49 lmc-weights-equil       = no [weight-equil-wl-delta,weight-equil-number-all-lambda,weight-equil-number-steps,weight-equil-number-samples,weight-equil-count-ratio]
  50                         The condition to set when we stop updating the weights and start to equilibrate.
  51
  52 lmc-weights-equil       = no [wl-delta,number-all-lambda,number-steps,number-samples,count-ratio]
  53 weight-equil-wl-delta           =  [positive real] stop when wl-delta drops below this value, if lmc-weights-equil  = wl-delta
  54 weight-equil-number-all-lambda  =  [positive integer] stop when we've got this number of samples at all lambda, if lmc-weights-equil  = number-all-lambda
  55 weight-equil-number-steps       =  [positive integer] stop when we've done this number of steps, if lmc-weights-equil  = number-steps
  56 weight-equil-number-samples     =  [positive integer] stop when we've done this number of samples, if lmc-weights-equil  = number-samples
  57 weight-equil-count-ratio        =  [positive real] stop when the ratio of min to max value is greater than this, if lmc-weights-equil  = count-ratio
  58
  59 fep-lambdas         = 0.0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1.0 (array of lambdas)
  60 mass-lambdas        = 0.0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1.0 (array of lambdas)
  61 coul-lambdas        = 0.0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1.0 (array of lambdas)
  62 vdw-lambdas         = 0.0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1.0 (array of lambdas)
  63 bonded-lambdas      = 0.0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1.0 (array of lambdas)
  64 restraint-lambdas   = 0.0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1.0 (array of lambdas)
  65 init-lambda-weights = 0.0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1.0 (initial weights for each lambda)
  66
  67 sc-alpha            = 0.5  [any positive real] (soft core alpha - now, only operates on vdw)
  68 sc-coul             = no [yes/no] -  Controls whether coulomb is also softocore
  69
  70 nstfep              = 20 [any integer multiple of nlist] (the
  71 frequency at which the other energies are determined, and the lambdas are updated.
  72 nstdgdl             = 200 [any integer multiple of nstfep] (the rate at which
  73 the dE terms are output to the dgdl file, per step)
  74 dhdl-print-energy   = Whether to print the energies as well as the energy differences (helps for some analysis codes --
  75                       for example, required when doing umbrella sampling with different temperatures)
  76 Notes on non-self-explanatory terms:
  77
  78 Weight methods:
  79 *barker-transition -- computes <alpha>_forward/<alpha>_reverse, with
  80 alpha being the Fermi function.  Like Bennett, but without self
  81 consistent iteration.  More efficient than metropolis transition
  82 *metropolis-transition, same as above but with the Metropolis function
  83 (min{1,exp(dE}) instead of the Fermi function
  84 *wang-landau- implements standard wang-landau
  85 *gibbs-wang-landau -- uses p(k)delta for nonlocal updating of weights, instead of delta
  86 *mbar.  Highly experimental, recommend against using.
  87 *minvar - optimizes the number of samples at each states with weights that lower the variance, as described
  88 by Escobedo (ref #???)
  89
  90 weight-equilibration criteria options:
  91
  92 *weight-equil-wl-delta - stop when wl-delta drops below this value set in 'weight-equil-wl-delta'.
  93 *weight-equil-number-all-lambda - stop when all the lambdas are equal or greater than 'weight-equil-number-all-lambda'
  94 *weight-equil-number-steps - stop when we've done the number of steps set in 'weight-equil-number-steps'
  95 *weight-equil-number-samples - stop when we've done the number of samples set in 'weight-equil-number-samples'
  96 *weight-equil-count-ratio - stop when the ratio of top to bottom value is greater than the value set in 'weight-equil-count-ratio'
  97
  98 wl-ratio -- For Wang-Landau, scaling is done when everything is
  99 within mc-ratio of the mean.
 100
 101 Right now, nstfep must be an integer multiple of nslist -- this can
 102 possibly be fixed later.  For now, since we update all the other
 103 energies every nstfep steps, it probably would be too slow if done
 104 more frequently.  Previously, I though if it was too fast, you could
 105 get an unlucky run to an unstable configuration, if it does not
 106 equilibrate the momenta at the new level.  I think I fixed a bug that
 107 was causing this problem.  It might be possible to go down to 10,
 108 though at that point, the cost starts to be large (takes about 20%
 109 longer with nstfep = 10 vs. 20 in my sample system).
 110
 111 This does bring up a point that perhaps for strict detailed balance, the
 112 momenta need to be resampled from a Boltzmann distribution at each
 113 transition.  This would slow down dynamics a lot, though -- is it
 114 necessary for complete detailed balance?  (Answer seems to be that if
 115 the velocity distribution function is equal in both states, as it is
 116 if the temperature is the same for both, then we're OK).
 117
 118 fep-lambdas is the default lambda array.  All the other four lambda
 119 arrays (except the serial tempering temperatures) use the values in
 120 this array if they are not specified.  If not specificed, it is
 121 assumed to be zero throughout.
 122
 123 For example, if you wanted to first to change coul, then vdw, while
 124 changing bonded at the same time as vdw, but the restraints throughout
 125 the first two thirds of the simulation, then you'd do something like:
 126
 127 init-lambda-state      = 0
 128 coul-lambdas           = 0.0 0.2 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0
 129 vdw-lambdas            = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
 130 bonded-lambdas         = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
 131 restraint-lambdas      = 0.0 0.0 0.1 0.2 0.3 0.5 0.7 1.0 1.0 1.0
 132
 133 This is equivalent to:
 134
 135 init-lambda-state      = 0
 136 fep-lambdas            = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
 137 coul-lambdas           = 0.0 0.2 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0
 138 restraint-lambdas      = 0.0 0.0 0.1 0.2 0.3 0.5 0.7 1.0 1.0 1.0
 139
 140 The fep-lambda array, in this case, is being used as the default to
 141 fill in the bonded and vdw lambda arrays.  Usually, it's best to fill
 142 in all array explicitly, just to make sure things are properly
 143 assigned.
 144
 145 If you want to turn on only restraints going from A to B, then it would be:
 146
 147 init-lambda-state      = 1
 148 restraint-lambdas      = 0.0 0.1 0.2 0.4 0.6 1.0
 149
 150 The current weights are written out to the logfile every nstlog steps,
 151 in units of kT.
 152
 153 You can load in weights using init-lambda-weights.  Use 0.00 for the first
 154 (zero because it's the reference) entries, then the rest of the weights,
 155 like:
 156
 157 init-lambda-weights =  0.00  0.54  1.45 3.45 6.78
 158
 159 However, most of the update methods don't do a good job with using
 160 given weights if there is no data added in previously.  So, it's best
 161 to use this for entering the weights for a subsequent equilibrium
 162 simulation.  I'll try to figure out if there is a better way.
 163
 164 Make sure you make the spacing in lambda sufficiently close.  It needs
 165 to be close in order to have good overlap.  I'd keep the free energy
 166 between states less that 5 kT or so, perhaps more like 2-3 kT.
 167
 168  ============ Coding ====================
 169
 170 The expanded ensemble changes (other than keeping track of global
 171 variables) are almost mostly contained in src/mdlib/sim_util.c, with
 172 additional changes in readir.c and mdebin.c
 173
 174 Things that would be nice to address:
 175
 176   -nstdgdl not necessarily a multiple of nslist (especially with nstlist = -1)
 177
 178 Test to do:
 179
 180  - throw a number of systems on to find an bugs.
 181
 182 Let me know if there is anything else I should be adding. . .
 183
 184 There's a very large number combinations of different parameters and
 185 approaches that could be adjusted in order to study what works best.
 186
 187 ======== Implementation Notes====
 188
 189 I'm still playing around trying to figure out the best combination of
 190 methods to use.  If the weights are poor, equilibration to a flat
 191 histogram in weights can be very slow.  Right now, I think that the
 192 Gibbs Wang-Landau method might be the fastest to equilibrate for less
 193 overlap.  In the limit of large numbers of samples, it probably
 194 doesn't matter very much.
 195
 196
 197 ======Notes on topologies for restraints=======
 198 ******This section is old, and may have some issues *****
 199 For distance and dihedral restraints, you need entries like this in
 200 the mdp -- for angle restraints, you don't.  This is a strange gromacs
 201 convention.
 202
 203 Dihedral restraint entries look like this -- note the A->B state, the
 204 A state and the B state (last two commented).
 205
 206 [ dihedral_restraints ]
 207 ;  i    j    k    l type label power phiA dphiA  kfacA  phiB  dphiB   kfacB
 208 ;mixed state
 209 1410 1393 1391 2610    1     0     2    38    0    0.00     38     0    41.84
 210 1393 1391 2610 2604    1     1     2   111    0    0.00    111     0    41.84
 211 1391 2610 2604 2606    1     2     2   -39    0    0.00    -39     0    41.84
 212 ; zero state
 213 ;1410 1393 1391 2610    1     0     2   38    0    0.00
 214 ;1393 1391 2610 2604    1     1     2  111    0    0.00
 215 ;1391 2610 2604 2606    1     2     2  -39    0    0.00
 216 ; one state
 217 ;1410 1393 1391 2610    1     0     2   38    0   41.84
 218 ;1393 1391 2610 2604    1     1     2  111    0   41.84
 219 ;1391 2610 2604 2606    1     2     2  -39    0   41.84
 220
 221 Free energy dependent angle restraints look like this:
 222
 223 [ angle_restraints ]
 224 ;  i    j    k    l type theta0A fcA mult  theta0B    fcB
 225 ;mixed state  -- need to have MULT listed twice!!!
 226 1393 1391 2610 1391    1   62      0.0  1     62    41.84    1
 227 1391 2610 2604 2610    1   79      0.0  1     79    41.84    1
 228 ; zero state
 229 ;1393 1391 2610 1391   1   62      0.0  1
 230 ;1391 2610 2604 2610   1   79      0.0  1
 231 ; one state
 232 ;1393 1391 2610 1391   1   62    41.84  1
 233 ;1391 2610 2604 2610   1   79    41.84  1
 234
 235 Note particularly that even though the multiplicity can't actually be
 236 changed in free energy perturbation, you need to have the multiplicity
 237 listed in both the A and B states.  This is a strange effect of using
 238 proper dihedral codes to do the angle restraints -- I tried to find a
 239 way to fix it, but gave up after a few hours -- it works using this
 240 topology framework.
 241
 242 Free energy dependent harmonic distance restraints are implemented as
 243 follows -- R2 is a ridiculously large number in order to obtain purely
 244 harmonic restraints -- otherwise, it becomes linear outside of R2.
 245
 246 [ distance_restraints ]
 247 ;  i    j type label typeprime    r0A   r1A   r2A    fcA  r0B   r1B
 248 r2B      fcB
 249 ;mixed state
 250 1391 2610     1     0         1   0.61  0.61   10.0     0.0  0.61  0.61 10.0  4184.0 ;zero state
 251 ;1391 2610    1     0         1   0.61  0.61   10.0     0.0 ;one state
 252 ;1391 2610    1     0         1   0.61  0.61    0.81 4184.0
 253
 254