Metadynamics sampling¶
CLEASE offers the possibility to calculate free energies as a function of arbitrary collective variables via metadynamics sampling. A common collective variable is the concentration of species, but variable you can think of will work. In short we want to calculate curves as shown below
Let’s move on to the details of how the sampling algorithm works. We start with the definition of the free energy
where k is the Boltzmann constant, T is the temperature, E is the energy of a configuration, Z is the partition function and \(\sigma\) denotes an atomic configuration. Hence, the sum runs over all possible configurations. Furthermore, the probability that the system is in a state \(\sigma\) is given by
The free energy at a given value for a general collective variable q is defined by
the function \(f(\sigma)\) is a mapping from an atomic configuration to the sought collective variable. It might for instance return the concentration of the atomic arrangement. \(\delta\) is a function that is 1 when \(q = f(\sigma)\) and zero otherwise. Thus, the difference is now that contributions from configurations that has a different value of the collective variable is cancelled out. Since, we now summed over all configurations that satisfy \(f(\sigma) = q\), the probability of finding the system in any state that satisfy \(f(\sigma) = q\) can be obtained by dividing by Z
Now, let’s see how the probabilities changes if we subtract an artificial potential V(q) that is only a function of the collective variable. First, we note that this potential can go inside the sum since the sum as only over configurations that has the same value for q. A new free energy F’ can therefore be defined as follows
by comparison it follows that the relation between the two free energies is
Similarly, the probability of occupying any configuration with \(f(\sigma) = q\) in the presence of an artificial potential is
from the above equation, we note that if we are able to select a potential that is such that it is exactly equal to the original free energy, the probability of being in a state satisfying \(f(\sigma) = q\) is
which is constant for all values of q! Hence, if we partition the domain of possible q values into bins, monitor how often the MC sampler visits each bin and adaptively tune the artificial potential V(q) until we visit all bins equally often, we know that we have found the free energy.
Carrying out a metadynamics calculation in practice¶
As before, we first need to define the settings. Let’s once again use our favorite example: AuCu!
>>> from clease.settings import CEBulk, Concentration
>>> conc = Concentration(basis_elements=[['Au', 'Cu']])
>>> settings = CEBulk(crystalstructure='fcc',
... a=3.8,
... supercell_factor=27,
... concentration=conc,
... db_name="aucu_metadyn.db",
... max_cluster_dia=[4.0])
The next thing we need to do is to load the ECIs and attach the calculator
>>> eci = {'c0': -1.0, 'c1_0': 0.1, 'c2_d0000_0_00': -0.2}
>>> atoms = settings.atoms.copy()*(5, 5, 5)
>>> from clease.calculator import attach_calculator
>>> atoms = attach_calculator(settings, atoms=atoms, eci=eci)
In pratice, the collective variables are calculated via one of the observers in CLEASE. If you plan to implement your own observers to use here, please note that there are certain requirements that needs to be satisfied if an observer should be applicable for metadynamics calculations.
The __call__ method needs to support a peak key word. Which is used to check what the collective variable is after a move, without actually performing the move
>>> def __call__(self, system_changes, peak=False):
... pass
It needs to have a method calculate_from_scratch that takes an atoms object as the only argument. This method is used to calculate the collective variable from scratch without making use of fast updates when the system_changes is known
>>> def calculate_from_scratch(self, atoms):
... pass
In this example we are going to use the concentration observer to track the concentration of Au
>>> from clease.montecarlo.observers import ConcentrationObserver
>>> obs = ConcentrationObserver(atoms, element='Au')
Next, we need to define a sampler. Since, the nature of the problem requires that the concentration can change, we will use the Semi-Grand Canonical ensemble
>>> from clease.montecarlo import SGCMonteCarlo
>>> mc = SGCMonteCarlo(atoms, 600, symbols=['Au', 'Cu'])
Then we need to define the artificial bias potential. Here, we are going to use a binned potential, which is a potential that is defined via values on a grid.
>>> from clease.montecarlo import BinnedBiasPotential
>>> bias = BinnedBiasPotential(xmin=0.0, xmax=1.0, nbins=60, getter=obs)
Here, the minimum concentration is set to 0 and the maximum concentration is set to 1, and the domain is partitioned into 60 bins. At last, we pass everything to the metadynamics sampler
>>> from clease.montecarlo import MetaDynamicsSampler
>>> meta_dyn = MetaDynamicsSampler(mc=mc, bias=bias, flat_limit=0.8, mod_factor=0.01,
... fname='aucu_metadyn.json')
>>> meta_dyn.run(max_sweeps=1)
The parameter flat_limit is a threshold used to determine if we have visited all the bins equally likely. In the above example, the algorithm will say that all bins have been visited equally likely if the bins with the fewest visits is visited at least 80% of the average.
The mod_factor tunes how much we should modify the artificial potential when the sampler visits a bin. It is given in units of kT, hence the artifial potential is altered by 0.01*kT everytime the sampler visits a bin. Finally, when we run we set here that the maximum number of sweeps is 1. This is only to avoid that the trial example takes too long running. This number should be much higher. If you set it None, the algorithm will run until it converges.
When you have managed to converge a calculation, you should reload the previous estimate, lower the modification factor and run again. Continue to lower the modification factor until the estimated free energy curve no longer changes.
To load an existing estimate, call this prior to passing the binned potential to the metadynamics sampler
>>> import json
>>> with open('aucu_metadyn.json', 'r') as f:
... data = json.load(f)
>>> bias.from_dict(data['bias_pot'])