FreeBird.jl Examples

Here we provide a few examples of how to use FreeBird.jl. We focus on how to set up a sampling run with various sampling schemes. You may need to fine-tune the parameters for your specific system.

Atomistic system

Before we start, we need to load the FreeBird.jl package.

using FreeBird

Nested Sampling

First, let's generate some initial configurations

configs = generate_initial_configs(120, 562.5, 6)

120-element Vector{AtomsBase.FastSystem{3, AtomsBase.PeriodicCell{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}, AtomsBase.ChemicalSpecies}}:
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 ⋮
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)
 FastSystem(H₆, periodicity = FFF)

Walkers are the objects that will be used in the simulation. They contain the information about the system's configuration. Let's warp the configurations into walkers.

walkers = AtomWalker.(generate_initial_configs(120, 562.5, 6))

Let's define the Lennard-Jones potential

lj = LJParameters(epsilon=0.1, sigma=2.5, cutoff=4.0)

LJParameters(0.1 eV, 2.5 Å, 4.0, -9.763240814208984e-5 eV)

We then create a LJAtomWalkers object that contains the walkers and the Lennard-Jones potential.

ls = LJAtomWalkers(walkers, lj)

We then set up a NestedSamplingParameters object with the desired parameters. See the documentation for more information.

ns_params = NestedSamplingParameters(mc_steps=200, step_size=0.1, random_seed=1234*rand(Int))

NestedSamplingParameters(200, 0.01, 0.1, 1.0e-6, 1.0, (0.25, 0.75), 0, 100, 1169912280906551876)

We then set up how we want to perform the Monte Carlo moves. MCRandomWalkClone is a routine that clones an existing walker and performs random walks.

mc = MCRandomWalkClone()

MCRandomWalkClone([1, 2, 3])

We also set up a SaveEveryN object for saving the output.

save = SaveEveryN(n_traj=10, n_snap=20_000, df_filename="output_df_lj6.csv", n_info=10)

SaveEveryN("output_df_lj6.csv", "output.traj.extxyz", "output.ls.extxyz", 10, 20000, 10)

Finally, we run the nested sampling simulation.

energies, liveset, _ = nested_sampling(ls, ns_params, 1_000, mc, save)

Metropolis Monte Carlo

Similarly to nested sampling, we need to define some MC parameters. We will use the same Lennard-Jones potential and initial configurations.

Let's use a walker from the previous example as the initial configuration.

at = ls.walkers[1]

AtomWalker{1}(
    configuration      : FastSystem(H₆, periodicity = FFF)
    energy             : -0.0451529120340907 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

Define the temperature grid, the number of equilibration steps, the number of sampling steps, and the step size.

temperatures = collect(1000.0:-50:0)
num_equilibration_steps = 100_000
num_sampling_steps = 100_000
step_size = 1.0

1.0

Pass the parameters to the MetropolisMCParameters object.

mc_params = MetropolisMCParameters(
    temperatures,
    equilibrium_steps=num_equilibration_steps,
    sampling_steps=num_sampling_steps,
    step_size=step_size,
    step_size_up=1.0,
    accept_range=(0.5,0.5)
)

MetropolisMCParameters([1000.0, 950.0, 900.0, 850.0, 800.0, 750.0, 700.0, 650.0, 600.0, 550.0  …  450.0, 400.0, 350.0, 300.0, 250.0, 200.0, 150.0, 100.0, 50.0, 0.0], 100000, 100000, 1.0, 0.001, 1.0, (0.5, 0.5), 1234)

Run the Monte Carlo simulation.

mc_energies, mc_ls, mc_cvs, acceptance_rates = monte_carlo_sampling(at, lj, mc_params)

Wang-Landau Sampling

It's also very easy to set up a Wang-Landau simulation. We will use the same Lennard-Jones potential and initial configuration.

at = ls.walkers[1]

AtomWalker{1}(
    configuration      : FastSystem(H₆, periodicity = FFF)
    energy             : -0.0451529120340907 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

Define the step size and the energy range.

step_size = 1.0
energy_min = -1.26
energy_max = 0.0
num_energy_bins = 1000

1000

Pass the parameters to the WangLandauParameters object.

wl_params = WangLandauParameters(num_steps=10_000,
                            energy_min=energy_min,
                            energy_max=energy_max,
                            num_energy_bins=num_energy_bins,
                            step_size=step_size,
                            max_iter=10000,
                            f_min=1.00001)

WangLandauParameters(10000, 0.8, 2.718281828459045, 1.00001, [-1.26, -1.2587387387387388, -1.2574774774774775, -1.2562162162162163, -1.254954954954955, -1.2536936936936938, -1.2524324324324325, -1.2511711711711713, -1.2499099099099098, -1.2486486486486486  …  -0.011351351351351352, -0.01009009009009009, -0.008828828828828829, -0.0075675675675675675, -0.006306306306306306, -0.005045045045045045, -0.0037837837837837837, -0.0025225225225225223, -0.0012612612612612612, 0.0], 10000, 1.0, 1234)

Viola! We can now run the Wang-Landau simulation.

energies_wl, configs, wl_params, S, H = wang_landau(at, lj, wl_params)

Lattice system

We can also perform simulations on lattice systems.

Let's set up a 2D square lattice with a supercell of dimensions 4x4x1. Make the first four sites occupied.

initial_lattice = SLattice{SquareLattice}(components=[[1,2,3,4]])

SLattice{SquareLattice}
    lattice_vectors      : [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    positions            : 16 grid points
    supercell_dimensions : (4, 4, 1)
    basis                : [(0.0, 0.0, 0.0)]
    periodicity          : (true, true, false)
    cutoff radii         : 2 nearest neighbors cutoffs [1.1, 1.5]
    occupations          : 
      ● ● ● ● 
      ○ ○ ○ ○ 
      ○ ○ ○ ○ 
      ○ ○ ○ ○ 
    adsorptions          : full adsorption

Define the Hamiltonian with the adsorption energy and the nearest-neighbor and next-nearest-neighbor energies.

adsorption_energy = -0.04
nn_energy = -0.01
nnn_energy = -0.0025
h = GenericLatticeHamiltonian(adsorption_energy, [nn_energy, nnn_energy], u"eV")

GenericLatticeHamiltonian{2,Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}}:
    on_site_interaction:      -0.04 eV
    nth_neighbor_interactions: [-0.01, -0.0025] eV

Wang-Landau Sampling

This time, let's do a Wang-Landau simulation first, as it's the most straightforward to set up.

energy_min = 20.5 * nn_energy + nn_energy / 8
energy_max = 16 * nn_energy - nn_energy / 8
num_energy_bins = 100

100

Pass the parameters to the WangLandauParameters object.

wl_params = WangLandauParameters(
    energy_min=energy_min,
    energy_max=energy_max,
    random_seed=Int(round(time() * 1000)),)

WangLandauParameters(100, 0.8, 2.718281828459045, 1.00000001, [-0.20625000000000002, -0.20577020202020205, -0.20529040404040405, -0.20481060606060608, -0.20433080808080809, -0.20385101010101012, -0.20337121212121215, -0.20289141414141415, -0.20241161616161618, -0.20193181818181818  …  -0.16306818181818183, -0.16258838383838384, -0.16210858585858587, -0.16162878787878787, -0.1611489898989899, -0.16066919191919193, -0.16018939393939394, -0.15970959595959597, -0.15922979797979797, -0.15875], 1000, 0.01, 1742236242747)

Run the Wang-Landau simulation.

energies_wl, configs, wl_params, S, H = wang_landau(initial_lattice, h, wl_params)

Metroplis Monte Carlo

Let's make a Metropolis Monte Carlo simulation on the same lattice system.

mc_lattice = deepcopy(initial_lattice)

SLattice{SquareLattice}
    lattice_vectors      : [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    positions            : 16 grid points
    supercell_dimensions : (4, 4, 1)
    basis                : [(0.0, 0.0, 0.0)]
    periodicity          : (true, true, false)
    cutoff radii         : 2 nearest neighbors cutoffs [1.1, 1.5]
    occupations          : 
      ● ● ● ● 
      ○ ○ ○ ○ 
      ○ ○ ○ ○ 
      ○ ○ ○ ○ 
    adsorptions          : full adsorption

Define the temperature grid, the number of equilibration steps, and the number of sampling steps.

temperatures = collect(200.0:-10:10)  # 200:-1:1
num_equilibration_steps = 25_000
num_sampling_steps = 25_000

25000

Pass the parameters to the MetropolisMCParameters object.

mc_params = MetropolisMCParameters(
    temperatures,
    equilibrium_steps=num_equilibration_steps,
    sampling_steps=num_sampling_steps,
    random_seed=Int(round(time()))
)

MetropolisMCParameters([200.0, 190.0, 180.0, 170.0, 160.0, 150.0, 140.0, 130.0, 120.0, 110.0, 100.0, 90.0, 80.0, 70.0, 60.0, 50.0, 40.0, 30.0, 20.0, 10.0], 25000, 25000, 0.01, 0.001, 1.0, (0.25, 0.75), 1742236243)

Run the Monte Carlo simulation.

mc_energies, mc_configs, mc_cvs, acceptance_rates = monte_carlo_sampling(mc_lattice, h, mc_params)

Nested Sampling

Finally, let's set up a nested sampling simulation. It's a bit more complicated than the previous two, and nested sampling is not the optimal method for lattice systems with many degenerate states.

Let's use the Distributions.jl package to generate the random initial configurations.

using Distributions

Make 1000 copies of the initial lattice.

walkers = [deepcopy(initial_lattice) for i in 1:1000]

Now using the sample function from the Distributions.jl package, we can randomly select 4 sites to be occupied.

for walker in walkers
    walker.components[1] = [false for i in 1:length(walker.components[1])]
    for i in sample(1:length(walker.components[1]), 4, replace=false)
        walker.components[1][i] = true
    end
end

Warp the walkers into LatticeWalker objects.

walkers = [LatticeWalker(walker) for walker in walkers]

Again, we use the same Hamiltonian as before. We then create a LatticeGasWalkers object that contains the walkers and the Hamiltonian.

ls = LatticeGasWalkers(walkers, h)

The rejection routine for nested sampling is the best choice for lattice systems.

mc = MCRejectionSampling()

MCRejectionSampling()

We then set up a LatticeNestedSamplingParameters object with the desired parameters. We will use the same parameters as before, but we will set the number of Monte Carlo steps to 100. We will also set the allowed fail count to 1000000, which is a bit high to make sure we explore the configuration space. We will also set the random seed to a random number.

ns_params = LatticeNestedSamplingParameters(mc_steps=100,allowed_fail_count=1_000_000,random_seed=Int(floor(1234*rand())))

LatticeNestedSamplingParameters(100, 1.0e-12, 0, 1000000, 98)

Redefine the output file name.

save = SaveEveryN(df_filename="output_ns_lattice2d.csv")

SaveEveryN("output_ns_lattice2d.csv", "output.traj.extxyz", "output.ls.extxyz", 100, 1000, 1)

And finally, we can run the nested sampling simulation

ns_energies, ls, _ = nested_sampling(ls, ns_params, 10_000, mc, save) # src

This page was generated using Literate.jl.