Adaptive Biasing Force Algorithm¶

Introduction¶

Relative to most other free energy methods, which are based on adding biases to the Hamiltonian of a system to obtain a flat histogram in converged sampling, the adaptive biasing force method (ABF) is unique, as it seeks a flattening of the generalized force. Like flat histogram methods, this bias force is applied until the landscape seen by a simulation is flat in collective variable (CV) space and freely diffusive sampling is achieved.

In practice, ABF functions by partitioning relevant CV space into histogram- like bins with a grid, and keeping a running tabulation of the number of visits to each bin. Concurrently, a running sum of the instantaneous generalized force experienced by the system. Putting these together gives an estimate for the mean generalized force, and integrating this quantity yields the free energy. Note that this means ABF returns a vector field rather than a free energy surface (which is the output of most methods.) Ref. [4] contains an excellent write-up on the method in general. The specific implementation in SSAGES is discussed in Ref. [6].

Note

An integrator for 1D, 2D and 3D surfaces are provided in SSAGES/Tools/ABF_integrator (requires numpy, scipy and matplotlib). Syntax is below; this is illustrated further in the tutorial found on this page.

./ABF_integrator.py -i <inputfile> -o <outputname> -p <bool> (<bool> <bool>) --interpolate <integer> (<integer> <integer>) --scale <float>

Implementation Notes¶

ABF calculates a generalized force on CVs at each timestep, and biases simulations using the negative of the estimated generalized force. This requires a grid, which requires that you define a CV range. Outside of the CV range, the simulation will continue to run, but the grid does not extend there, so there is no bias applied and no histogram hits are collected. histogram hits will be collected.

ABF can restart from a previous run. Simply include Fworld_cvX and Nworld outputs generated by the previous run in your working directory, and set the "restart" option to true.

Warning

ABF keeps only a single backup of old files, and overwrites older backups with new data if they already exist. If restarting, a backup will not be created, instead ABF will read from and update the newest files. If you want to keep a copy of the output from a previous run after restart, be sure to rename the output or place a copy in a different directory.

If you are using the multiple walkers option, they will read from and write to the same histogram and force estimate during runtime. The resulting histogram and force data are saved in the same way as single-walker simulations.

ABF can optionally define a restraint range, which biases simulations back toward the region of interest using a harmonic restraint with user-chosen spring constant(s). To disable restraints, enter a spring constant k equal to or less than zero.

Warning

The restraint range should be WIDER than the CV range by at least one bin size in each direction.

If restraints are used on a periodic system, one can define the periodic boundaries, so that minimum image convention to CVs can be applied using the commands CV_periodic_boundary_upper_bounds and CV_periodic_boundary_lower_bounds. For example, on a \(-\pi\) to \(\pi\) CV, if the CV is restrained between -3.14 to -2.36 and the system crosses the periodic boundary, setting this will ensure the restraint is applied correctly back towards -3.14 rather than applying an incorrect a large force to push it toward -2.36.

Example Input¶

"methods" : [
    {
        "type" : "ABF",
        "cvs" : [0,1],
        "CV_lower_bounds" : [-3.14, -3.14],
        "CV_upper_bounds" : [3.14, 3.14],
        "CV_bins" : [21,21],
        "CV_restraint_minimums" : [-5,-5],
        "CV_restraint_maximums" : [5,5],
        "CV_restraint_spring_constants" : [0,0],
        "CV_isperiodic" : [false,false],
        "timestep" : 0.002,
        "minimum_count" : 50,
        "output_file" : "F_out",
        "output_frequency" : 1000,
        "unit_conversion" : 1,
        "frequency" : 1
    }
]

Warning

Be sure to follow correct JSON syntax for your input, with a comma after every line except the last within each bracket.

Options & Parameters¶

Define ABF:

In the methods block, define the ABF method through the syntax:

"type" : "ABF"

Define CVs

To define the collective variables:

"cvs" : [0,1]

In the example input, this defines a two-dimensional CV-space to be sampled by ABF, with indices [0,1]. The argument to this must be a list of integers defining the CVs to be operated on by ABF.

Define the grid

To define the bounds:

"CV_lower_bounds" : [-3.14, -3.14]
"CV_upper_bounds" : [3.14, 3.14]

Thee are arrays of doubles whose length is the number of CVs used. This defines the minimum and maximum values for the CVs for the range in which the method will be used in order.

To define the number of CV bins used:

"CV_bins" : [21,21]

This array of integers defines the number of histogram bins in each CV dimension in order.

Define the restraints

"CV_restraint_minimums" : [-5,-5],
"CV_restraint_maximums" : [5,5],

These arrays define the minimum and maximum values for the CV restraints in order.

"CV_restraint_spring_constants" : [0,0],

This array defines the spring constants for the CV restraints in order. Enter a value equal to or less than zero to turn restraints off.

"CV_isperiodic" : [false,false],

This array defines whether a given CV is periodic for restraint purposes.
This is only used to apply minimum image convention to CV restraints. The
value can be safely set to ``false`` *even for periodic CVs* if no
restraints are being used.

Warning

If ANY CV is set to periodic, then CV_periodic_boundary_lower_bounds and CV_periodic_boundary_upper_bounds must be provided for ALL CVs. Values entered for non-periodic CVs are not used.

"CV_periodic_boundary_lower_bounds" : [-3.14, -3.14],
"CV_periodic_boundary_upper_bounds" : [3.14, 3.14],

These arrays define the lower and upper end of the period. This only matters if CV_isperiodic is true for the CV.

Define time and unit parameters

"timestep" : 0.002,

The timestep of the simulation. Units depend on the conversion factor that follows. This must be entered correctly, otherwise the generalized force estimate will be incorrect.

"unit_conversion" : 1,

Defines the unit conversion from d(momentum)/d(time) to force for the simulation. For LAMMPS using units real, this is 2390.06 (gram.angstrom/mole.femtosecond^2 -> kcal/mole.angstrom). For GROMACS, this is 1.

"minimum_count" : 50,

This is the number of hits required to a bin in the general histogram before the full bias is active. Below this value, the bias linearly decreases to equal 0 at hits = 0. Default = 200, but user should provide a reasonable value for their system. See [4] and [6] for more details.

Output parameters

"output_frequency" : 1000,

Optional: This defines how many timesteps pass in between output of the generalized force.

"output_file" : "F_out",

This is a string value defining the file name for the adaptive vector force field that is acquired. The default name is “F_out”.

"Fworld_output_file" : "Fworld_cv"

Optional: This is the name of the file to backup raw Fworld force output for use in restarts. There will be separate outputs for each CV. The default filename is Fworld_cv, which saves each CV’s output to Fworld_cvX.

"Nworld_output_file" : "Nworld"

Optional: This is name of the file which backs up the raw histogram data for restart purposes. The default filename is “Nworld”.

Optional Parameters

"mass_weighting" : false,

Turns on/off mass weighing of the adaptive force. The default is false, which turns off the weighting.

Warning

Leave this off if your system has massless sites such as in TIP4P water.

"restart" : false

This boolean determines whether the simulation is a restart. The default value is false. If set to true, ABF will attempt to load a previous state from Nworld and Fworld files.

"frequency" : 1

Leave at 1.

Output¶

The main output of the method is stored in a file specified in ‘filename’. This file will contain the Adaptive Force vector field printed out every ‘backup_frequency’ steps and at the end of a simulation. The method outputs a vector field, with vectors defined on each point on a grid that goes from (CV_lower_bounds) to (CV_upper_bounds) of each CV in its dimension, with (CV_bins) of grid points in each dimension. For example, for 2 CVs defined from (-1,1) and (-1,0) with 3 and 2 bins respectively would be a 3x2 grid (6 grid points). The printout is in the following format: 2*N number of columns, where N is the number of CVs. First N columns are coordinates in CV space, the N+1 to 2N columns are components of the Adaptive Force vectors. An example for N=2 is:

CV1 Coord	CV2 Coord	d(A)/d(CV1)	d(A)/d(CV2)
-1	-1	-1	1
-1	0	2	1
0	-1	1	2
0	0	2	3
1	-1	2	4
1	0	3	5

Tutorial¶

Alanine Dipeptide

For LAMMPS (must be built with RIGID and MOLECULE packages) To build RIGID and MOLECULE:

Go to LAMMPS src folder (/build/hooks/lammps/lammps-download-prefix/src/lammps-download/src/ for -DLAMMPS=YES)
Do:

make yes-RIGID
make yes-MOLECULE

Go to your build folder and make.

Find the following input files in Examples/User/ABF/Example_AlanineDipeptide:

in.ADP_ABF_Example(0-1) (2 files)
example.input
ADP_ABF_1walker.json
ADP_ABF_2walkers.json

Put the contents of ABF_ADP_LAMMPS_Example folder in your ssages build folder
For a single walker example, do:

./ssages ADP_ABF_1walker.json.json

For 2 walkers, do:

mpirun -np 2 ./ssages ADP_ABF_2walkers.json

For GROMACS:

Optional:

adp.gro
topol.top
nvt.mdp

Required:

example_adp(0-1).tpr (2 files)
ADP_ABF_1walker.json
ADP_ABF_2walkers.json

Put the contents of ABF_ADP_Gromacs_Example in your ssages build folder
For a single walker example, do:

./ssages ABF_ADP_1walker.json

For 2 walkers, do:

mpirun -np 2 ./ssages ABF_ADP_2walkers.json

These will run using the pre-prepared input files in .tpr format. If you wish to prepare the input files yourself using GROMACS tools (if compiled with -DGROMACS=YES):

/build/hooks/gromacs/gromacs/bin/gmx_ssages grompp -f nvt.mdp -p topol.top -c adp.gro -o example_adp0.tpr
/build/hooks/gromacs/gromacs/bin/gmx_ssages grompp -f nvt.mdp -p topol.top -c adp.gro -o example_adp1.tpr

Be sure to change the seed in .mdp files for random velocity generation, so walkers can explore different places on the free energy surface.

Multiple walkers initiated from different seeds will explore different regions and will all contribute to the same adaptive force.

After the run is finished, you can check that your output matches the sample outputs given in the examples folders:

Copy ABF_integrator.py (requires numpy, scipy and matplotlib) into your build folder.
Run the integrator:

python ABF_integrator.py --periodic1 True --periodic2 True --interpolate 200

This will output a contour map, a gradient field and a heatmap. Compare these to the sample outputs.

Sodium Chloride

For LAMMPS (must be built with KSPACE and MOLECULE packages) To build RIGID and MOLECULE:

Go to LAMMPS src folder (/build/hooks/lammps/lammps-download-prefix/src/lammps-download/src/ for -DLAMMPS=YES)
Do:

make yes-KSPACE
make yes-MOLECULE

Go to your build folder and make.

Find the following input files in Examples/User/ABF/Example_NaCl/ABF_NaCl_LAMMPS_Example:

in.NaCl_ADP_example(0-1) (2 files)
data.spce
ADP_NaCl_1walker.json
ADP_NaCl_2walkers.json

Put the contents of ABF_NaCl_LAMMPS_Example folder in your ssages build folder
For a single walker example, do:

./ssages ADP_NaCl_1walker.json.json

For 2 walkers, do:

mpirun -np 2 ./ssages ADP_NaCl_2walkers.json

For GROMACS:

Optional:

NaCl.gro
topol.top
npt.mdp

Required:

example_NaCl(0-1).tpr (2 files)
ADP_NaCl_1walker.json
ADP_NaCl_2walkers.json

Put the contents of ABF_NaCl_Gromacs_Example in your ssages build folder
For a single walker example, do:

./ssages ABF_NaCl_1walker.json

For 2 walkers, do:

mpirun -np 2 ./ssages ABF_NaCl_2walkers.json

These will run using the pre-prepared input files in .tpr format. If you wish to prepare the input files yourself using GROMACS tools (if compiled with -DGROMACS=YES):

/build/hooks/gromacs/gromacs/bin/gmx_ssages grompp -f npt.mdp -p topol.top -c NaCl.gro -o example_NaCl0.tpr
/build/hooks/gromacs/gromacs/bin/gmx_ssages grompp -f npt.mdp -p topol.top -c NaCl.gro -o example_NaCl1.tpr

Be sure to change the seed in .mdp files for random velocity generation, so walkers can explore different places on the free energy surface.

Multiple walkers initiated from different seeds will explore different regions and will all contribute to the same adaptive force.

After the run is finished, you can check that your output matches the sample outputs given in the examples folders:

Copy ABF_integrator.py (requires numpy, scipy and matplotlib) into your build folder.
Run the integrator:

python ABF_integrator.py

This will output a Potential of Mean Force graph. Compare this to the sample output.

Developers¶

Emre Sevgen
Hythem Sidky