Nested sampling method

Authors:Noam, Rob, Gabor, Livia
Date:2015

Overview

The nested sampling method was first introduced by John Skilling [1] in the field of applied probability and inference to efficiently sample probability densities in high-dimensional spaces where the regions contributing most of the probability mass are exponentially localized. It is an annealing algorithm that creates a sequence of probability distributions, each more concentrated than the previous one near the high-likelihood region of the parameter space - i.e. in the context of materials, the low-energy region of configuration space.

Since its original inception, nested sampling has also been applied to atomistic systems, and its several advantages mean it became a powerful method to sample atomic configuration spaces.

  • calculate the partition function, hence all thermodynamic properties become accessible
    • calculate the heat capacity and locate the phase transitions
    • with an order parameter calculate the free energy
  • the sampling process itself is independent of temperature, thus the calculation of thermodynamic quantities is a simple post-processing step
  • no prior knowledge is needed of the phases or phase transitions
  • can be used with both constant volume and constant pressure calculations
  • calculate the entire phase diagram in an automated way
  • considerable computational gain over parallel tempering
[1]
  1. Skilling, in Nested Sampling, edited by Rainer Fischer, Roland Preuss, and Udo von Toussaint, AIP Conf. Proc. No. 735 (AIP, New York, 2004), p. 395

Publications

L.B. Partay, A.P. Bartok, G. Csanyi, Efficient Sampling of Atomic Configurational Spaces, J. Phys. Chem. B (2010), 114, 10502–10512, http://pubs.acs.org/doi/abs/10.1021/jp1012973

L.B. Partay, A.P. Bartok, G. Csanyi, Nested sampling for materials: The case of hard spheres, Phys. Rev. E (2014), 89, 022302, http://journals.aps.org/pre/abstract/10.1103/PhysRevE.89.022302

R.J.N. Baldock, L.B. Partay, A.P. Bartok, M.C. Payne, G. Csanyi, Determining pressure-temperature phase diagrams of materials, Phys. Rev. B (2016), 93, 174108, http://journals.aps.org/prb/abstract/10.1103/PhysRevB.93.174108

Algorithm

Nested sampling is an iterative method… Top-down, Size of live set, random walk…etc.

pymatnest

The pymatnest package is a software library written in Fortran 95/python for the purpose of carrying out nested sampling calculations with a variety of options suitable for different systems. It can be used with the supplied fortran potential models and it also has interfaces with the following packages:

  • LAMMPS
  • QUIP

MC and MD step algorithms

rej_free_perturb_velo
  • if atom_velo_rej_free_fully_randomize, pick random velocities consistent with Emax
  • else perturb velocities
    • if current velocity=0, can’t rescale, so pick random velocities consistent with Emax
    • else, pick new random magnitude consistent with Emax, random rotation of current direction with angle uniform in +/- atom_velo_rej_free_perturb_angle

beta\_from\_current\_ns\_pe\_values This routine computes beta from the NS estimate of $`frac{\Delta X}{\Delta E}$. ($E=$potential energy.)

rej\_free\_canonical\_velo
  • pick random velocities from cannonical distribution with inverse T currentbeta
  • If movement_args[‘atom_momentum_retain_fraction’]>0.0
    • Partial momentum update
    • $p’ = alpha p + ( 1- alpha^2 )^{frac{1}{2}}n$ www.mcmchandbook.net eq. (5.34)
    • Else if “rescale_only” rescale velocities to current temperature
    • Else completely refresh velocities
do\_MC\_atom\_velo\_walk
  • If separable_MDNS = T , call rej_free_canonical_velo()
  • Else if MC_atom_velo_walk_rej_free
    • call rej_free_perturb_velo()
  • else do MC pertubation to velocities
do_MD_atom_walk
  • if MD_atom_velo_pre_perturb, call do_MC_atom_velo_walk() for magnitude and rotation
  • propagate in time atom_traj_len time steps of length MD_atom_timestep
  • If MD_atom_reject_energy_violation is set, accept/reject entire move on E deviating by less than MD_atom_energy_fuzz times kinetic energy
  • accept/reject entire move on E < Emax and KE < KEmax
  • if reject
    • set positions, velocities, energy back to value before perturbation (maybe should be after?)
  • else
    • flip velocities if MD_atom_velo_flip_accept
  • if MD_atom_velo_post_perturb, call do_MC_atom_velo_walk() for magnitude and rotation
do_MC_atom_walk
  • if MC_atom_velocities and MC_atom_velocities_pre_perturb, call do_MC_atom_velo_walk() to perturb velocities, magnitude and and rotation
  • if using fortran calculator and not reproducible
    • call fortran MC code f_MC_MD.MC_atom_walk
  • else
    • do python MC
    • if MC_atom_Galilean
      • go Galilean MC in python
    • else
      • loop atom_traj_len times
        • loop over atoms in random order
          • propose single atom move
          • accept/reject on E < Emax
do_MC_swap_step
  • return if all atomic numbers are identical
  • randomly pick a desired cluster size
  • pick two clusters with distinct atomic numbers, backing off to smaller clusters on failure to find appropriate larger clusters, but always pick at least a pair of atoms to swap
  • accept swap if energy < Emax

do_MC_cell_volume_step

do_MC_cell_shear_step

do_MC_cell_stretch_step

do_atom_walk
  • loop n_atom_steps_per_call times, calling do_MC_atom_walk() or do_MD_atom_walk()
walk_single_walker
  • create list
    • do_atom_walk * n_atom_step_n_calls
    • do_cell_volume_step * n_cell_volume_steps
    • do_cell_shear_step * n_cell_shear_steps
    • do_cell_stretch_step * n_cell_stretch_steps
    • do_swap_step * n_swap_steps
  • loop while n_model_calls_used < n_model_calls
    • pick random item from list
    • do move
  • perturb final energy by random_energy_perturbation
full_auto_set_stepsizes
  • Step sizes for each (H)MC move are set via a loop which performs additional exploration moves, calibrating each step size to obtain an acceptance rate inside a specified range.
  • The routine is MPI parallelised, so that the wall time goes as 1/num_of_processes
  • For each (H)MC move type the following is performed
    • Set ‘’movement_args’’ parameters so that only one (H)MC call is made at a time
    • Min and max acceptance rates are copied from parameters MC_adjust_min_rate / MD_adjust_min_rate and MC_adjust_max_rate / MD_adjust_max_rate
    • Step size calibration loop:
      • Repeat the following 200/num_of_MPI_processes times:
        • Copy a configuration from the live set (each MPI process chooses a different configuration)
        • Each MPI processes performs one (H)MC move on its cloned configuration running statistics for the number of accepted/rejected moves on each process are recorded
      • The total number of accepted/rejected moves for this step size (summed across all MPI processes) are estabilshed
        • If the total acceptance rate is within the desired range, return this stepsize
          • If this is NOT the first time the loop has been performed for this (H)MC move AND we previously obtained an acceptance rate on one side of the desired range, and now find an acceptance rate on the other side of the desired range
            • Return the step size that gave an acceptance rate closest to the middle of the desired range.
        • Store step length and acceptance rate
        • update step length, by * or / by MC_adjust_step_factor, to improve acceptance rate
        • Check that step size is not larger than max allowed value (specified by user), and also that step size is not smaller than 10-20(useful for detecting errors).
  • Return step sizes and time taken for routine to run
main: parse arguments
  • process n_atom_steps
    • If break_up_atom_traj
      • n_atom_steps_per_call = 1
      • n_atom_steps_n_calls = n_atom_steps
    • else
      • n_atom_steps_per_call = n_atom_steps
      • n_atom_steps_n_calls = 1