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 highdimensional 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 highlikelihood region of the parameter space  i.e. in the context of materials, the lowenergy 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 postprocessing 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] 

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 pressuretemperature 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… Topdown, 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