Monte Carlo Sampling

CLEASE currently support two ensembles for Monte Carlo sampling: canonical and semi-grand canonical ensembles. A canonical ensemble has a fixed number of atoms, concentration and temperature while a semi-grand canonical ensemble has a fixed number of atoms, temperature and chemical potential. To use a fitted CE model to run MC sampling we first initialise small cell holding the nessecary information about the lattice and the clusters

1from clease.settings import CEBulk, Concentration
2conc = Concentration(basis_elements=[['Au', 'Cu']])
3settings = CEBulk(crystalstructure='fcc',
4                  a=3.8,
5                  supercell_factor=27,
6                  concentration=conc,
7                  db_name="aucu.db",
8                  max_cluster_dia=[6.0, 5.0])

Next, we need to specify a set if ECIs. These can for instance be loaded from a file, but here we hard code them for simplicity

1eci = {'c0': -1.0, 'c1_0': 0.1, 'c2_d0000_0_00': -0.2}

For efficient initialisation of large cells, CLEASE comes with a convenient helper function called attach_calculator. We create our MC cell by repeating the atoms object of the settings.

1from clease.calculator import attach_calculator
2atoms = settings.atoms.copy()*(5, 5, 5)
3atoms = attach_calculator(settings, atoms=atoms, eci=eci)

Let's insert a few Cu atoms

1atoms[0].symbol = 'Cu'
2atoms[1].symbol = 'Cu'
3atoms[2].symbol = 'Cu'

We are now ready to run a MC calculation

1from clease.montecarlo import Montecarlo
2T = 500
3mc = Montecarlo(atoms, T)
4mc.run(steps=1000)

After a MC run, you can retrieve internal energy, heat capacity etc. by calling

1thermo = mc.get_thermodynamic_quantities()

Monitoring a MC run

In many cases it is useful to be able to monitor the evolution of parameters during a run, and not simply getting the quantities after the run is finished. A good example can be to monitor the evolution of the energy in order to determine whether the system has reached equilibrium. CLEASE comes with a special set of classes called MCObservers for this task. As an example, we can store a value for the energy every 100 iteration by

1from clease.montecarlo.observers import EnergyEvolution
2obs = EnergyEvolution(mc)
3mc.attach(obs, interval=100)
4mc.run(steps=1000)
5energies = obs.energies

Another useful observer is the Snapshot observer. This observers takes snaptshots of the configuration at regular intervals and stores them in a trajectory file.

1from clease.montecarlo.observers import Snapshot
2snap = Snapshot(atoms, fname='snapshot')
3mc.attach(snap, interval=200)
4mc.run(steps=1000)

There are many more observers distributes with CLEASE, for a complete list check the API documentation.

Constraining the MC sampling

In some cases you might want to prevent certain moves to occur. That can for instance be that certain elements should remain fixed. CLEASE offers the possibility to impose arbitrary constraint via its MCConstraint functionality. MCConstraints can be added in a very similar fashion as the observers. To fix one element

1from clease.montecarlo.constraints import FixedElement
2cnst = FixedElement('Cu')
3mc.generator.add_constraint(cnst)

Note, that the usage of a constraint in this system is a bit weird as it has only two elements. Hence, fixing one prevents any move from happening. But the point here is just to illustrate how a constraint can be attached.

Note: If your system has multiple basis, you most likely want to add a [ConstrainSwapByBasis] constraint object, in order to avoid swaps happening across different basis sites. The Montecarlo object will not automatically avoid cross-basis swaps.

Implementing Your Own Observer

You can implement your own observer and monitor whatever quantity you might be interested in. To to so you can create your own class that inherits from the base MCObserver class. To illustrate the usage, let's create an observers that monitor how many Cu atoms there are on average in each (100) layer!

Before we initialise this monitor we need to make sure that the tag of each atom represents the corresponding layer.

1from clease.montecarlo.observers import MCObserver
2from ase.geometry import get_layers
3class LayerMonitor(MCObserver):
4   def __init__(self, atoms):
5       self.layers, _ = get_layers(atoms, [1, 0, 0])
6       self.layer_average = [0 for _ in set(self.layers)]
7       self.num_calls = 1
8       # Initialise the structure
9       for atom in atoms:
10           if atom.symbol == 'Cu':
11               self.layer_average[self.layers[atom.index]] += 1
12   def observe_step(self, step):
13       self.num_calls += 1
14       system_changes = step.last_change
15       for change in system_changes:
16           layer = self.layers[change[0]]
17           if change[2] == 'Cu':
18               self.layer_average[layer] += 1
19           if change[1] == 'Cu':
20               self.layer_average[layer] -= 1
21   def get_averages(self):
22       return {'layer{}'.format(i): x/self.num_calls for i, x in enumerate(self.layer_average)}

When this observer is attached, the `observe_step` method will be executed on every Monte Carlo step. The call signature takes in a [MCStep] instance. The system_changes variable here is a list of the following form `[(10, Au, Cu), (34, Cu, Au)]` which means that the symbol on site 10 changes from Au to Cu and the symbol on site 34 changes from Cu to Au. Hence, in the update algorithm above we check if the last element of a single change is equal to Cu, if so we know that there is one additional Cu atom in the new layer. And if the middle element of a change is equal to Cu, there is one less atom in the corresponding layer. Note that if a MC move is rejected the system_changes will typically be `[(10, Au, Au), (34, Cu, Cu)]`. The get_averages function returns a dictionary. This method is optimal to implement, but if it is implemented the result will automatically be added to the result of get_thermodynamic_quantities.

To use this observer in our calculation

1monitor = LayerMonitor(atoms)
2mc = Montecarlo(atoms, T)
3mc.attach(monitor, interval=1)
4mc.run(steps=1000)

There are a few other methods that can be useful to implement. First, the reset method. This method can be invoked if the reset method of the mc calculation is called.

Implementing Your Own Constraints

If you want to have custom constraints on MC moves, CLEASE lets you implement your own. The idea is to create a class that inherits from the base MCConstraint class and has a function __call__ that returns True if a move is valid and False if a move is not valid. To illustrate this, let's say that we want the atoms on sites less that 25 to remain fixed. The reason for doing so, can be that you have a set of indices that you know constitutes a surface and you want to keep them fixed.

1from clease.montecarlo.constraints import MCConstraint
2class FixedIndices(MCConstraint):
3   def __call__(self, system_changes):
4       for change in system_changes:
5           if change.index <= 25:
6               return False
7       return True

To use this constrain in our calculation

1cnst = FixedIndices()
2mc.generator.add_constraint(cnst)
3mc.run(steps=1000)

Sampling the SGC Ensemble

CLEASE also gives the possibility to perform MC sampling in the semi grand canonical ensemble. Everything that has to do with observers and constraints mentioned above can also be used together with this class. To run a calcualtion in the SGC ensemble

1from clease.montecarlo import SGCMonteCarlo
2sgc_mc = SGCMonteCarlo(atoms, T, symbols=['Au', 'Cu'])
3sgc_mc.run(steps=1000, chem_pot={'c1_0': -0.15})

The chem_pot parameter sets the chemical potentials. It is possible to set one chemical potential for each singlet correlation function (i.e. ECIs that starts with c1).