Expanded Ensemble mode

Expanded ensemble (EE) mode can be used for calculations of solvation free energies by gradual removal or insertion of molecules in the system. The theory and computational details behind the expanded ensemble method are described in papers (J.Chem.Phys.96, 1776 (1992); Molec.Simulations , 18, 43 (1996); J.Chem.Phys.,108,227 (1998)). Implementation of the EE approach here differs from the cited above articles, in a way that interaction parameters of the choosen solute molecule with the rest of the system are scaled with the coupling parameter $\alpha$, $0 \le \alpha \le 1$ according the folowing rules:

$$\sigma = const$$ (12)

$$\epsilon \to \alpha^{n1} \epsilon$$

$$q \to \alpha^{n2} q$$

The coupling parameter $\alpha$ assumes a number of fixed values between 0 and 1, where ``1'' corresponds to the solute molecule properly inserted in the solvent, while value ``0'' correspons to the fully eleminated solute (which then can be considered as an a gas phase). The program computes probability distribution over subensembles with different values of $\alpha$, from which the free energy deffirence can be obtained. Additional bias for transition probabilities between subensembles is given by the so-called balancing factors, which need to be tuned in order to provide reasonable homogenious distribution of probabilities over subensembles. The Wang-Landau procedure is by default used for tuning of the balancing factors, but there is also possibility for manual tuning.

Note 1: Special care should be taken in MDEE simulations of molecules having zero Lennard-Jones potential for hydrogen atoms (like SPC or TIP3P water). Simulation may become unstable when repulsive potential between oxygens becoms weeker while electrostatic attraction between hydrogen and oxygen still present. It is advisable to prescribe some small Lennard-Jones potential to such hydrogens.

Note 2: It is always recommended to use option ``Cut_forces'' with threshold level 0.05 - 0.1 in the case of EE simulation to avoid ``collisions'' with ``almost deleted'' particles

The expanded ensemble is specified by the following keywords:

• MDEE <num_ensembles> <mode>

  Run expanded ensemble simulations with <num_ensembles> subensembles, Parameter <mode> can be ``auto'' ($\alpha$ poins are set with equal separation between them within the [0:1] interval, and balancing factors are tuned automatically by the Wang-Landau algorithm, Wang F., Landau D. P., Phys. Rev. Lett. 86, P. 2050-2053 (2001)), ``WL'' ($\alpha$ points are specified by the user while balancing factors are optimized by the Wang-Landau procedure, or ``manual'' (manual setup of EE parameters).

  In the case of manual or WL mode the MDEE command should be followed by exactly num_ensembles lines, with two number on each: the value of coupling parameter $\alpha$ and balancing factor, corresponding to this subensemble. Values of balancing factors are given per inserted molecule in units of ``$kT$''.

  In the case of automatic setup of the expanded ensemble parameters, the set of coupling parameters $\alpha$ is choosen uniformly in the range [0:1] (that is, $\alpha_m = (M-m)/(M-1)$ where $M$ is the total number of subensembles <num_ensembles>.

  Default mode: auto

  In the case of automatic optimization of the balancing factors (mode auto or WL), additional keywords can be used:

• Num_pass <value>

  Defines how many times the system should walk between the extreme subemsembles and back in order to complete one WL iteration. Default is 2.

• Num_WL <value>

  Defines the number of WL iterations. After this number of WL iterations is carried out, the optimization of balancing factors is completed, and the program proceeds with EE procedure with fixed balancing factors. Default is 10 iteration.

• WL_parameters <increment> <factor>

  Sets up initial increment in the Wang-Landau optimization of balancing factors, and scaling factor by which the increment is scaled after each WL iteration.

  Default values: increment = 0.1, factor = 0.5.

  In case of automatic WL tuning of the balancing factors (mode auto) the program starts expended ensemble simulation with zero balancing factors. in case of mode WL the initial values of balancing factors are introduced from the input. After each attempt of visiting a subensemble, an increment is added to the corresponding balancing factor. After walking between extreme subensembles at least Num_pass times, a WL series (sweep) finishes and the value of the increment is decreased by a factor (default value 0.5, see keyword WL_parameters).

  Notes:

  1. If WL procedure did not perform necessary number of iterations (keyword Num_WL), the program given warning in the output ``Expanded ensemble simulation not converged yet''. Then the computed free energy values may be not accurate or not relevant et al. It is always recommended to continue simulation long enough after the balancing factors are optimized.

  2. While evaluating results, ensure that the system has walked at least several times between the extreeme subensembles (the first and the last). Look at the ``table of transition'' in the end of the output.

  3. Pay attention at the acceptance ratios for transitions between subensembles, they can be found in the last two columns of the section ``Distribution over subensembles'' in the output. The acceptance ratios should not be too low. If they are two low (a few percent or less), the number of subensembles should be increased.

• EE_typ <i_typ>

  Specifies that the molecules of type i_typ are those which are inserted/deleted in the EE procedure. This keyword is required!

• EE_start-in <start_in>

  Start EE simulation in subensemble number <start_in>. Default: in the last subensemble (that is, with completely deleted EE molecule(s))

• EE_freq <steps> Try EE transition after steps MD steps.

  Default: 10

• EE_new

  Change balancing factors after restart (manual mode only). Otherwize old values written in the restart file are used.

• EE_scaling <n1> <n2>

  Powers of scaling of Lennerd-Jones and electrostatic interactions, see (12) Default: n1=4, n2=2

• Probability_threshold <P_low> Assing value P_low to computed subensembles probabilities which are lower than this number

  Default: 0.001