doc/src/compute_fep.txt

   1 "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
   2
   3 :link(lws,http://lammps.sandia.gov)
   4 :link(ld,Manual.html)
   5 :link(lc,Section_commands.html#comm)
   6
   7 :line
   8
   9 compute fep command :h3
  10
  11 [Syntax:]
  12
  13 compute ID group-ID fep temp attribute args ... keyword value ... :pre
  14
  15 ID, group-ID are documented in the "compute"_compute.html command :ulb,l
  16 fep = name of this compute command :l
  17 temp = external temperature (as specified for constant-temperature run) :l
  18 one or more attributes with args may be appended :l
  19 attribute = {pair} or {atom} :l
  20   {pair} args = pstyle pparam I J v_delta
  21     pstyle = pair style name, e.g. lj/cut
  22     pparam = parameter to perturb
  23     I,J = type pair(s) to set parameter for
  24     v_delta = variable with perturbation to apply (in the units of the parameter)
  25   {atom} args = aparam I v_delta
  26     aparam = parameter to perturb
  27     I = type to set parameter for
  28     v_delta = variable with perturbation to apply (in the units of the parameter) :pre
  29 zero or more keyword/value pairs may be appended :l
  30 keyword = {tail} or {volume} :l
  31   {tail} value = {no} or {yes}
  32     {no} = ignore tail correction to pair energies (usually small in fep)
  33     {yes} = include tail correction to pair energies
  34   {volume} value = {no} or {yes}
  35     {no} = ignore volume changes (e.g. in {NVE} or {NVT} trajectories)
  36     {yes} = include volume changes (e.g. in {NpT} trajectories) :pre
  37 :ule
  38
  39 [Examples:]
  40
  41 compute 1 all fep 298 pair lj/cut epsilon 1 * v_delta pair lj/cut sigma 1 * v_delta volume yes
  42 compute 1 all fep 300 atom charge 2 v_delta :pre
  43
  44
  45 [Description:]
  46
  47 Apply a perturbation to parameters of the interaction potential and
  48 recalculate the pair potential energy without changing the atomic
  49 coordinates from those of the reference, unperturbed system. This
  50 compute can be used to calculate free energy differences using several
  51 methods, such as free-energy perturbation (FEP), finite-difference
  52 thermodynamic integration (FDTI) or Bennet's acceptance ratio method
  53 (BAR).
  54
  55 The potential energy of the system is decomposed in three terms: a
  56 background term corresponding to interaction sites whose parameters
  57 remain constant, a reference term \(U_0\) corresponding to the
  58 initial interactions of the atoms that will undergo perturbation, and
  59 a term \(U_1\) corresponding to the final interactions of
  60 these atoms:
  61
  62 :c,image(Eqs/compute_fep_u.jpg)
  63
  64 A coupling parameter \(\lambda\) varying from 0 to 1 connects the
  65 reference and perturbed systems:
  66
  67 :c,image(Eqs/compute_fep_lambda.jpg)
  68
  69 It is possible but not necessary that the coupling parameter (or a
  70 function thereof) appears as a multiplication factor of the potential
  71 energy. Therefore, this compute can apply perturbations to interaction
  72 parameters that are not directly proportional to the potential energy
  73 (e.g. \(\sigma\) in Lennard-Jones potentials).
  74
  75 This command can be combined with "fix adapt"_fix_adapt.html to
  76 perform multistage free-energy perturbation calculations along
  77 stepwise alchemical transformations during a simulation run:
  78
  79 :c,image(Eqs/compute_fep_fep.jpg)
  80
  81 This compute is suitable for the finite-difference thermodynamic
  82 integration (FDTI) method "(Mezei)"_#Mezei, which is based on an
  83 evaluation of the numerical derivative of the free energy by a
  84 perturbation method using a very small \(\delta\):
  85
  86 :c,image(Eqs/compute_fep_fdti.jpg)
  87
  88 where \(w_i\) are weights of a numerical quadrature. The "fix
  89 adapt"_fix_adapt.html command can be used to define the stages of
  90 \(\lambda\) at which the derivative is calculated and averaged.
  91
  92 The compute fep calculates the exponential Boltzmann term and also the
  93 potential energy difference \(U_1 -U_0\). By
  94 choosing a very small perturbation \(\delta\) the thermodynamic
  95 integration method can be implemented using a numerical evaluation of
  96 the derivative of the potential energy with respect to \(\lambda\):
  97
  98 :c,image(Eqs/compute_fep_ti.jpg)
  99
 100 Another technique to calculate free energy differences is the
 101 acceptance ratio method "(Bennet)"_#Bennet, which can be implemented
 102 by calculating the potential energy differences with \(\delta\) = 1.0 on
 103 both the forward and reverse routes:
 104
 105 :c,image(Eqs/compute_fep_bar.jpg)
 106
 107 The value of the free energy difference is determined by numerical
 108 root finding to establish the equality.
 109
 110 Concerning the choice of how the atomic parameters are perturbed in
 111 order to setup an alchemical transformation route, several strategies
 112 are available, such as single-topology or double-topology strategies
 113 "(Pearlman)"_#Pearlman. The latter does not require modification of
 114 bond lengths, angles or other internal coordinates.
 115
 116 NOTES: This compute command does not take kinetic energy into account,
 117 therefore the masses of the particles should not be modified between
 118 the reference and perturbed states, or along the alchemical
 119 transformation route.  This compute command does not change bond
 120 lengths or other internal coordinates "(Boresch,
 121 Karplus)"_#BoreschKarplus.
 122
 123 :line
 124
 125 The {pair} attribute enables various parameters of potentials defined
 126 by the "pair_style"_pair_style.html and "pair_coeff"_pair_coeff.html
 127 commands to be changed, if the pair style supports it.
 128
 129 The {pstyle} argument is the name of the pair style. For example,
 130 {pstyle} could be specified as "lj/cut".  The {pparam} argument is the
 131 name of the parameter to change.  This is a (non-exclusive) list of
 132 pair styles and parameters that can be used with this compute.  See
 133 the doc pages for individual pair styles and their energy formulas for
 134 the meaning of these parameters:
 135
 136 "lj/cut"_pair_lj.html: epsilon,sigma: type pairs:
 137 "lj/cut/coul/cut"_pair_lj.html: epsilon,sigma: type pairs:
 138 "lj/cut/coul/long"_pair_lj.html: epsilon,sigma: type pairs:
 139 "lj/cut/soft"_pair_lj_soft.html: epsilon,sigma,lambda: type pairs:
 140 "coul/cut/soft"_pair_lj_soft.html: lambda: type pairs:
 141 "coul/long/soft"_pair_lj_soft.html: lambda: type pairs:
 142 "lj/cut/coul/cut/soft"_pair_lj_soft.html: epsilon,sigma,lambda: type pairs:
 143 "lj/cut/coul/long/soft"_pair_lj_soft.html: epsilon,sigma,lambda: type pairs:
 144 "lj/cut/tip4p/long/soft"_pair_lj_soft.html: epsilon,sigma,lambda: type pairs:
 145 "tip4p/long/soft"_pair_lj_soft.html: lambda: type pairs:
 146 "lj/charmm/coul/long/soft"_pair_lj_soft.html: epsilon,sigma,lambda: type pairs:
 147 "born"_pair_born.html: a,b,c: type pairs:
 148 "buck"_pair_buck.html: a,c : type pairs :tb(c=3,s=:)
 149
 150 Note that it is easy to add new potentials and their parameters to
 151 this list.  All it typically takes is adding an extract() method to
 152 the pair_*.cpp file associated with the potential.
 153
 154 Similar to the "pair_coeff"_pair_coeff.html command, I and J can be
 155 specified in one of two ways.  Explicit numeric values can be used for
 156 each, as in the 1st example above.  I <= J is required.  LAMMPS sets
 157 the coefficients for the symmetric J,I interaction to the same
 158 values. A wild-card asterisk can be used in place of or in conjunction
 159 with the I,J arguments to set the coefficients for multiple pairs of
 160 atom types.  This takes the form "*" or "*n" or "n*" or "m*n".  If N =
 161 the number of atom types, then an asterisk with no numeric values
 162 means all types from 1 to N.  A leading asterisk means all types from
 163 1 to n (inclusive).  A trailing asterisk means all types from n to N
 164 (inclusive).  A middle asterisk means all types from m to n
 165 (inclusive).  Note that only type pairs with I <= J are considered; if
 166 asterisks imply type pairs where J < I, they are ignored.
 167
 168 If "pair_style hybrid or hybrid/overlay"_pair_hybrid.html is being
 169 used, then the {pstyle} will be a sub-style name.  You must specify
 170 I,J arguments that correspond to type pair values defined (via the
 171 "pair_coeff"_pair_coeff.html command) for that sub-style.
 172
 173 The {v_name} argument for keyword {pair} is the name of an
 174 "equal-style variable"_variable.html which will be evaluated each time
 175 this compute is invoked.  It should be specified as v_name, where name
 176 is the variable name.
 177
 178 :line
 179
 180 The {atom} attribute enables atom properties to be changed.  The
 181 {aparam} argument is the name of the parameter to change.  This is the
 182 current list of atom parameters that can be used with this compute:
 183
 184 charge = charge on particle :ul
 185
 186 The {v_name} argument for keyword {pair} is the name of an
 187 "equal-style variable"_variable.html which will be evaluated each time
 188 this compute is invoked.  It should be specified as v_name, where name
 189 is the variable name.
 190
 191 :line
 192
 193 The {tail} keyword controls the calculation of the tail correction to
 194 "van der Waals" pair energies beyond the cutoff, if this has been
 195 activated via the "pair_modify"_pair_modify.html command. If the
 196 perturbation is small, the tail contribution to the energy difference
 197 between the reference and perturbed systems should be negligible.
 198
 199 If the keyword {volume} = {yes}, then the Boltzmann term is multiplied
 200 by the volume so that correct ensemble averaging can be performed over
 201 trajectories during which the volume fluctuates or changes "(Allen and
 202 Tildesley)"_#AllenTildesley:
 203
 204 :c,image(Eqs/compute_fep_vol.jpg)
 205
 206
 207 :line
 208
 209 [Output info:]
 210
 211 This compute calculates a global vector of length 3 which contains the
 212 energy difference ( \(U_1-U_0\) ) as c_ID\[1\], the
 213 Boltzmann factor \(\exp(-(U_1-U_0)/kT)\), or
 214 \(V \exp(-(U_1-U_0)/kT)\), as c_ID\[2\] and the
 215 volume of the simulation box \(V\) as c_ID\[3\]. \(U_1\) is the
 216 pair potential energy obtained with the perturbed parameters and
 217 \(U_0\) is the pair potential energy obtained with the
 218 unperturbed parameters. The energies include kspace terms if these
 219 are used in the simulation.
 220
 221 These output results can be used by any command that uses a global
 222 scalar or vector from a compute as input.  See "Section
 223 6.15"_Section_howto.html#howto_15 for an overview of LAMMPS output
 224 options. For example, the computed values can be averaged using "fix
 225 ave/time"_fix_ave_time.html.
 226
 227 The values calculated by this compute are "extensive".
 228
 229
 230 [Restrictions:]
 231
 232 This compute is distributed as the USER-FEP package.  It is only
 233 enabled if LAMMPS was built with that package.  See the "Making
 234 LAMMPS"_Section_start.html#start_3 section for more info.
 235
 236 [Related commands:]
 237
 238 "fix adapt/fep"_fix_adapt_fep.html, "fix ave/time"_fix_ave_time.html,
 239 "pair_style lj/soft/coul/soft"_pair_lj_soft.html
 240
 241 [Default:]
 242
 243 The option defaults are {tail} = {no}, {volume} = {no}.
 244
 245 :line
 246
 247 :link(Pearlman)
 248 [(Pearlman)] Pearlman, J Chem Phys, 98, 1487 (1994)
 249
 250 :link(Mezei)
 251 [(Mezei)] Mezei, J Chem Phys, 86, 7084 (1987)
 252
 253 :link(Bennet)
 254 [(Bennet)] Bennet, J Comput Phys, 22, 245 (1976)
 255
 256 :link(BoreschKarplus)
 257 [(BoreschKarplus)] Boresch and Karplus, J Phys Chem A, 103, 103 (1999)
 258
 259 :link(AllenTildesley)
 260 [(AllenTildesley)] Allen and Tildesley, Computer Simulation of
 261 Liquids, Oxford University Press (1987)
 262
 263