GitHub

Example 1: Melting a 13-Atoms Neon Cluster

This is an example calculation for finding the melting temperature of a 13-atoms neon cluster using a parallel-tempering Monte Carlo simulation. First, we load ParallelTemperingMonteCarlo and Plots:

using Plots
using ParallelTemperingMonteCarlo
using LaTeXStrings

Setting up the Model

Firstly, we set the number of atoms:

n_atoms = 13;

Next, we define the potential, here an extended Lennard Jones potential with coefficients for even powers of r, starting from powers of -6 to -16:

c=[-10.5097942564988, 989.725135614556, -101383.865938807, 3918846.12841668, -56234083.4334278, 288738837.441765]
pot = ELJPotentialEven{6}(c)
ELJPotentialEven{6, Float64}([-10.5097942564988, 989.725135614556, -101383.865938807, 3.91884612841668e6, -5.62340834334278e7, 2.88738837441765e8])

We further have to define the starting configuration of the simulation. For Ne13 we choose the icosahedral ground state of Ne13 (from the Cambridge cluster database). The atomic positions are given in Angstrom, which are then converted to Bohr radii (as the program uses atomic units).

pos_ne13 = [[2.825384495892464, 0.928562467914040, 0.505520149314310],
[2.023342172678102,	-2.136126268595355, 0.666071287554958],
[2.033761811732818,	-0.643989413759464, -2.133000349161121],
[0.979777205108572,	2.312002562803556, -1.671909307631893],
[0.962914279874254,	-0.102326586625353, 2.857083360096907],
[0.317957619634043,	2.646768968413408, 1.412132053672896],
[-2.825388342924982, -0.928563755928189, -0.505520471387560],
[-0.317955944853142, -2.646769840660271, -1.412131825293682],
[-0.979776174195320, -2.312003751825495, 1.671909138648006],
[-0.962916072888105, 0.102326392265998,	-2.857083272537599],
[-2.023340541398004, 2.136128558801072,	-0.666071089291685],
[-2.033762834001679, 0.643989905095452, 2.132999911364582],
[0.000002325340981,	0.000000762100600, 0.000000414930733]];
AtoBohr = 1.8897259886;
pos_ne13 = pos_ne13 * AtoBohr
13-element Vector{Vector{Float64}}:
 [5.339202509675499, 1.7547286276557148, 0.9552945639202041]
 [3.8235622875401982, -4.036693324695786, 1.2586922223528678]
 [3.8432525502537267, -1.2169635315645375, -4.030786193502644]
 [1.8515104475315411, 4.369051328639683, -3.1594504692142205]
 [1.8196441394724319, -0.19336921007065874, 5.399104677171738]
 [0.6008527770958447, 5.001668105430829, 2.668542641160762]
 [-5.339209779512827, -1.754731061649526, -0.9552951725503948]
 [-0.6008496122188508, -5.001669753738395, -2.668542209586626]
 [-1.851508499387977, -4.369053575565342, 3.1594501498809775]
 [-1.8196475277773039, 0.19336884278473446, -5.399104511708638]
 [-3.823559204867802, 4.036697652557049, -1.2586918476896083]
 [-3.8432544820617607, 1.2169644600549232, 4.030785366187147]
 [4.394257284152319e-6, 1.4401613097476533e-6, 7.841053896189476e-7]

Finally, we have to choose appropriate boundary conditions, here spherical boundary conditions (solid boundary around the cluster), to suppress atom loss processes. Finding this radius is a non-trivial task, and has to be chosen and tested carefully. A radius chosen too small will exert artificial pressure on the cluster while a too large value leads to atoms being ejected.

bc_ne13 = SphericalBC(radius=5.32*AtoBohr)
SphericalBC{Float64}(101.06969058367278)

We package the initial configuration and boundary conditions into a Config struct:

start_config = Config(pos_ne13, bc_ne13)
Config{13, SphericalBC{Float64}, Float64}(StaticArraysCore.SVector{3, Float64}[[5.339202509675499, 1.7547286276557148, 0.9552945639202041], [3.8235622875401982, -4.036693324695786, 1.2586922223528678], [3.8432525502537267, -1.2169635315645375, -4.030786193502644], [1.8515104475315411, 4.369051328639683, -3.1594504692142205], [1.8196441394724319, -0.19336921007065874, 5.399104677171738], [0.6008527770958447, 5.001668105430829, 2.668542641160762], [-5.339209779512827, -1.754731061649526, -0.9552951725503948], [-0.6008496122188508, -5.001669753738395, -2.668542209586626], [-1.851508499387977, -4.369053575565342, 3.1594501498809775], [-1.8196475277773039, 0.19336884278473446, -5.399104511708638], [-3.823559204867802, 4.036697652557049, -1.2586918476896083], [-3.8432544820617607, 1.2169644600549232, 4.030785366187147], [4.394257284152319e-6, 1.4401613097476533e-6, 7.841053896189476e-7]], SphericalBC{Float64}(101.06969058367278))

Setting up the Simulation Parameters

We first set the temperature grid, which defines the range of temperatures we consider. This is done by defining the upper and lower temperature limits, along with the number of temperatures (also called trajectories) we want to sample. Note, that a geometrical distribution of temperatures is chosen to maximise overlaps in the energy histograms.

ti = 4.;
tf = 16.;
n_traj = 25;
temp = TempGrid{n_traj}(ti,tf)
TempGrid{25, Float64}([4.0, 4.237852377437181, 4.489848193237492, 4.756828460010884, 5.039684199579493, 5.3393594166801375, 5.656854249492381, 5.993228307506726, 6.3496042078727974, 6.727171322029716, 7.127189745122714, 7.550994501453547, 8.0, 8.475704754874362, 8.979696386474984, 9.513656920021768, 10.079368399158986, 10.678718833360275, 11.313708498984761, 11.986456615013452, 12.699208415745597, 13.454342644059432, 14.254379490245428, 15.101989002907096, 16.0], [78943.7463157743, 74512.97429431944, 70330.88239791946, 66383.51328363325, 62657.692973985395, 59140.986890041044, 55821.65835215453, 52688.62940989218, 49731.44387040245, 46940.232401844325, 44305.679595409776, 41818.99287601273, 39471.87315788715, 37256.48714715972, 35165.44119895973, 33191.756641816624, 31328.846486992697, 29570.493445020522, 27910.829176077266, 26344.31470494609, 24865.721935201218, 23470.116200922163, 22152.839797704888, 20909.49643800636, 19735.936578943576])

Now we set the hyperparameters for this simulation:

  • mc_cycles is the number of Monte Carlo cycles we want to run, the longer the more accurate but also more expensive.
  • mc_sample is the number of MC cycles after which the energy of the state is recorded.
  • max_displ_atom determines the maximum displacement of an atom over a cycle. The maximum displacement is automatically adjusted in the program guaranteeing a 40-60% acceptance rate.
  • n_adjust is the number of moves after which the step size of atom moves is adjusted.
mc_cycles = 100000;
mc_sample = 1;
displ_atom = 0.1;
max_displ_atom = [0.1*sqrt(displ_atom*temp.t_grid[i]) for i in 1:n_traj];
n_adjust = 100;

For neatness, all parameters are collected in a MCParams struct:

mc_params = MCParams(mc_cycles, n_traj, n_atoms, mc_sample = mc_sample, n_adjust = n_adjust)
MCParams(100000, 20000, 1, 25, 13, 100, 100, 0.4, 0.6)

We then define the ensemble, here we are using the NVT ensemble (keeping N, the number of atoms, V, the volume, and T, the temperature constant). This allows us to derive a MoveStrategy to feed into the PTMC simulation. Here, we do n_atoms atom displacements of randomy chosen atoms per Monte Carlo cycle.

ensemble = NVT(n_atoms);
move_strat = MoveStrategy(ensemble)
MoveStrategy{13, NVT}(NVT(13, 13, 0), ["atommove", "atommove", "atommove", "atommove", "atommove", "atommove", "atommove", "atommove", "atommove", "atommove", "atommove", "atommove", "atommove"])

Running the Simulation

Finally, we run the simulation. This method returns the current state and results of the simulation. The data is stored in various local files created in the current working directory.

mc_states, results = ptmc_run!(mc_params,temp,start_config,pot,ensemble;save=1000);
params set
equilibration complete
MC loop done.
[19.470371917175708, 20.03812515635061, 20.341734753083536, 20.29236091413404, 20.464365124634217, 20.605044582479856, 21.6760871383111, 21.88246979155003, 22.069535559868307, 23.856478017397755, 24.20040417770317, 23.758640093954295, 24.95990909383904, 26.809221683132506, 28.66498765078339, 34.98868506105601, 44.22043246108181, 54.56024616379064, 64.92080974768346, 79.36199917724007, 78.46125883027908, 68.42015628607878, 53.63974968033182, 43.8236439772934, 40.03111224410605]
done

Post-processing and Analyzing of Results

The raw heat capacity plot is obtained from:

plot(temp.t_grid,results.heat_cap, xlabel="T [K]", ylabel=L"\(C_v\) [A.U]", title="Raw Heat Capacity against Temperature")

and the energy histograms by:

data = [results.en_histogram[i] for i in 1:n_traj]
plot(data)

For post-processing of the data we use the multihistogram method. This method accesses the stored data created from the ptmc_run! method and returns values for the energies, histogram data, temperature, partition function, heat capacity, heat capacity gradient, and entropy, which can be plotted as shown:

energies,histogramdata,T,Z,Cv,dCv,S = postprocess();

Plot of heat capacity against temperature:

plot(T,Cv,label=nothing, xlabel="T [K]", ylabel=L"\(C_v\) [A.U]", title="Heat capacity against Temperature")

Plot of the heat capacity gradient against temperature:

plot(T,dCv,label=nothing, xlabel="T [K]", ylabel=L"\(\frac{dC_v}{dT}\) [A.U]", title="Heat capacity gradient against Temperature")

Plot of the partition function against temperature:

plot(T,Z,label=nothing, xlabel="T [K]", ylabel="Z [1]", title="Partition function against Temperature")

Plot of the entropy of the system against temperature:

plot(T,S,label=nothing, xlabel="T [K]", ylabel="S [A.U]", title="Entropy against Temperature")

This page was generated using Literate.jl.