Grand Canonical Monte Carlo (GCMC)¶

General info¶

About Monte Carlo / the Grand Canonical Ensemble

It is best to read a bit about Monte Carlo and ensembles before working with the GCMC code. Almost every book or review text on molecular simulations will do, for example: Frenkel D, Smit B. Understanding molecular simulation: from algorithms to applications. Academic Press; 2002. 672 p.

Wikipedia also has some pages of interest:

• http://en.wikipedia.org/wiki/Monte_Carlo_method
• http://en.wikipedia.org/wiki/Grand_canonical_ensemble

It is important to note that this method heavily relies on random numbers, and simulations are thus non-repeatable in detail, but should converge to the same answer.

About the AMS GCMC code

The GCMC code was originally developed for standalone Reaxff by Thomas Senftle, working as a Graduate Student at Penn State University under the supervision of Dr. Adri van Duin [1] [2]. The original version was a wrapper code that called an external executable to perform the reaxff minimization step and energy calculation, and relied on file modification and parsing to steer the reaxff code and get the results back.

The code was later rewritten by Hans van Schoot (SCM), in close collaboration with Thomas Senftle, to integrate it directly into the ADF-ReaxFF code. The current version is an AMS re-implementation so the method can be used with almost any engine supported by AMS (support for 3D periodic boundary conditions by the engine is currently a requirement).

Method Details¶

Overview

The GCMC method will perform a number of Grand Canonical Monte Carlo trial moves (set by the Iterations option of the GCMC input block), and accept or reject them based on the energy produced by the geometry optimization of the trial geometry for the given engine. The Monte Carlo algorithm will always accept a step if it results in a decrease of the energy, and accept steps that go up in energy with a probability. This section will give some details about how the method works.

MC Moves (Insert/Delete/Displace/ChangeVolume)

The GCMC method currently supports 4 types of MC Moves: Insert, Delete, Displace (sometimes also called Move), and ChangeVolume. The first three MC moves are always available and the ChangeVolume becomes available only ChangeVolume option is set to True. The first three move types change coordinates of atoms in the system, while the ChangeVolume move changes the lattice only.

On every MC iteration, the method first selects one of the molecules defined by the Molecule input blocks at random and then selects an applicable MC move type. If there are no molecules of this type in the system then no Delete or Displace is attempted. If the selected molecule has the NoAddRemove option set then the Insert and Delete moves will not be attempted. If no MC move is possible with the selected molecule type then another one is selected or a VolumeChange is attempted, if allowed. If no moves with any of the provided molecules is possible (i.e. if all molecules have NoAddRemove set to True, there is nothing to displace and the volume is fixed) then the program will stop.

The Insert and Displace MC move will rotate the molecule randomly and put it at a random position, and then check if the minimum interatomic distance between the molecule and the rest of the system is within the MinDistance and MaxDistance boundaries. If the condition is not satisfied, a new set of coordinates is generated and the check is performed again. This is repeated up to NumAttempts times before stopping with an error.

The volume change is controlled by the VolumeChangeMax keyword. This sets the volume change limit, and it should be a value between between 0 and 1. The new volume will be calculated as: V_new = exp(random(-1:1)*VolumeChangeMax)*V_old.

Calculating the chemical potential

The chemical potential of the molecule (or atom) reservoir is used when calculating the Boltzmann accept/reject criteria after a MC move is executed. This value can be derived from first principles using statistical mechanics, or equivalently, it can be determined from thermochemical tables available in literature sources.

For example, the proper chemical potential for a GCMC simulation in which single oxygen atoms are exchanged with a reservoir of O2 gas, should be equal to 1/2 the chemical potential of O2 at the temperature and pressure of the reservoir [1]:

\(\mu^{MC}_{O}(T,P) = \frac {1} {2} \mu^{MC}_{O2}(T,P) = \frac {1} {2} \left [ \mu^{ref}_{O2}(T,P_{ref}) + kT ln \left ( \frac {P} {P_{ref}} \right ) - E^{diss}_{O2} \right ]\)

where the reference chemical potential \(\mu^{ref}_{O2}(T,P_{ref})\) is the experimentally determined chemical potential of O2 at T and P_ref, \(kT ln \left ( \frac {P} {P_{ref}} \right )\) is the pressure correction to the free energy, and \(E^{diss}_{O2}\) is the dissociation energy of the O₂ molecule.

Calculating energies

Because the GCMC simulation adds and deletes atoms or molecules during the runtime, it cannot directly compare the AMS energies for the MC acceptance criteria: inserting a molecule will usually lower the total energy of the system, causing the MC to always accept it, and always reject a deletion. To compensate this, the GCMC method calculates a “corrected” MC energy to compare the trial energy with, consisting of the previously accepted AMS energy and a correction depending on the move:

\(E_{old}^{MC} = E_{old}^{AMS} + \mu^{MC}\) for an Insert move;

\(E_{old}^{MC} = E_{old}^{AMS} - \mu^{MC}\) for a Delete move;

\(E_{old}^{MC} = E_{old}^{AMS}\) for a Displace move;

\(E_{old}^{MC} = E_{old}^{AMS} - P (V_{new} - V_{old}) + N_{inserted} ln \left ( \frac {V_{new}^{avail}}{V_{old}^{avail}} \right ) kT\) for a ChangeVolume move.

Here, \(\mu^{MC}\) is the chemical potential of the inserted/deleted molecule, P is the pressure, V is the volume, and N_inserted is the total number of MC molecules. The “new” and “old” subscripts refer to the current and the last accepted values. The V^avail values are calculated from the MC-available volume as described below.

Calculating volumes

The available volume can be calculated in a few different ways, depending on the VolumeOption setting:

• Free: volume = total volume - occupied volume - specified vacuum volume (NonAccessibleVolume)
• Total: volume = total cell volume
• Accessible: volume = specified accessible volume (AccessibleVolume)
• FreeAccessible: volume = specified accessible volume (AccessibleVolume) - occupied volume

Here, the occupied volume is calculated as a sum of volumes of atoms that do not belong to the MC part of the system (i.e. that were not inserted during calculation and are not Removables). The volume of an atom is calculated using the average of the covalent and the Van der Waals radii of the atom defined in the atominfo module used throughout AMS.

The AccessibleVolume and NonAccessibleVolume keywords can be used to get a more accurate available volume.

Acceptance criteria

An MC move is always accepted if the AMS energy is lower than the corrected MC energy of the last accepted MC move, or if the energy increase is small enough. If the new energy is higher, the code generates a random number between 0 and 1, and accepts the move if the random number is greater than:

prob = preFactor * exp(-Beta*deltaE)

The prefactor is calculated (for insert and delete moves) using the deBroglie wavelength of the inserted molecules, the number of inserted molecules and the available MC volume of the system.

Input¶

The GCMC functionality in AMS is triggered using the following Task key:

Task GCMC

GCMC
   AccessibleVolume float
   Box
      Amax float
      Amin float
      Bmax float
      Bmin float
      Cmax float
      Cmin float
   End
   Ensemble [Mu-VT | Mu-PT]
   Iterations integer
   MapAtomsToOriginalCell Yes/No
   MaxDistance float
   MinDistance float
   Molecule
      ChemicalPotential float
      NoAddRemove Yes/No
      SystemName string
   End
   NonAccessibleVolume float
   NumAttempts integer
   Pressure float
   Removables # Non-standard block. See details.
      ...
   End
   Restart string
   SwapAtoms
      Probability float
      Regions string
   End
   Temperature float
   UseGCPreFactor Yes/No
   VolumeChangeMax float
   VolumeOption [Free | Total | Accessible | FreeAccessible]
End

The following keys are common for all GCMC calculations and should always be specified. The ChemicalPotential value should correspond to the \(\mu^{MC}\) expression above, and not to the experimental chemical potential \(\mu^{ref}\), which means it should include the (engine-dependent) free molecule’s energy.

GCMC

Molecule

Type:	Block
Recurring:	True
GUI name:	Molecules
Description:	This block defines the molecule (or atom) that can be inserted/moved/deleted with the MC method. The coordinates should form a reasonable structure. The MC code uses these coordinates during the insertion step by giving them a random rotation, followed by a random translation to generate a random position of the molecule inside the box. Currently, there is no check to make sure all atoms of the molecule stay inside the simulation box. The program does check that the MaxDistance/MinDistance conditions are satisfied.

ChemicalPotential

Type:	Float
Unit:	Hartree
Description:	Chemical potential of the molecule (or atom) reservoir. It is used when calculating the Boltzmann accept/reject criteria after a MC move is executed. This value can be derived from first principles using statistical mechanics, or equivalently, it can be determined from thermochemical tables available in literature sources. For example, the proper chemical potential for a GCMC simulation in which single oxygen atoms are exchanged with a reservoir of O2 gas, should equal 1/2 the chemical potential of O2 at the temperature and pressure of the reservoir: cmpot = Mu_O(T,P) = 1/2*Mu_O2(T,P) = 1/2 * [Mu_ref(T,P_ref) + kT*Log(P/Pref) - E_diss] where the reference chemical potential [Mu_ref(T,P_ref)] is the experimentally determined chemical potential of O2 at T and Pref; kT*Log(P/Pref) is the pressure correction to the free energy, and E_diss is the dissociation energy of the O2 molecule.

NoAddRemove

Type:	Bool
Default value:	No
GUI name:	Fix molecule
Description:	Set to True to tell the GCMC code to keep the number of molecules/atoms of this type fixed. It will thus disable Insert/Delete moves on this type, meaning it can only do a displacement move, or volume change move (for an NPT ensemble).

SystemName

Type:	String
GUI name:	Molecule
Description:	String ID of a named [System] to be inserted. The lattice specified with this System, if any, is ignored and the main system’s lattice is used instead.

Iterations

Type:	Integer
GUI name:	Number of GCMC iterations
Description:	Number of GCMC moves.

Temperature

Type:	Float
Default value:	300.0
Unit:	Kelvin
Description:	Temperature of the simulation. Increase the temperature to improve the chance of accepting steps that result in a higher energy.

The following keys are related to Insert and Displace moves.

GCMC

NumAttempts

Type:	Integer
Default value:	1000
GUI name:	Max tries
Description:	Try inserting/moving the selected molecule up to the specified number of times or until all constraints are satisfied. If all attempts fail a message will be printed and the simulation will stop. If the MaxDistance-MinDistance interval is small this number may have to be large.

MinDistance

Type:	Float
Default value:	0.3
Unit:	Angstrom
GUI name:	Add molecules not closer than
Description:	Keep the minimal distance to other atoms of the system when adding the molecule.

MaxDistance

Type:	Float
Default value:	3.0
Unit:	Angstrom
GUI name:	Add molecules within
Description:	The max distance to other atoms of the system when adding the molecule.

The following keys influence computation of the acceptance probability and of the MC energy correction.

GCMC

UseGCPreFactor

Type:	Bool
Default value:	Yes
GUI name:	Use GC prefactor
Description:	Use the GC pre-exponential factor for probability.

VolumeOption

Type:	Multiple Choice
Default value:	Free
Options:	[Free, Total, Accessible, FreeAccessible]
GUI name:	Volume method
Description:	Specifies the method to calculate the volume used to calculate the GC pre-exponential factor and the energy correction in the Mu-PT ensemble: Free: V = totalVolume - occupiedVolume - NonAccessibleVolume; Total: V = totalVolume; Accessible: V = AccessibleVolume; FreeAccessible: V = AccessibleVolume - occupiedVolume. The AccessibleVolume and NonAccessibleVolume are specified in the input, the occupiedVolume is calculated as a sum of atomic volumes.

AccessibleVolume

Type:	Float
Default value:	0.0
Description:	Volume available to GCMC, in cubic Angstroms. AccessibleVolume should be specified for “Accessible” and “FreeAccessible” [VolumeOption].

NonAccessibleVolume

Type:	Float
Default value:	0.0
GUI name:	Non-accessible volume
Description:	Volume not available to GCMC, in cubic Angstroms. NonAccessibleVolume may be specified for the “Free” [VolumeOption] to reduce the accessible volume.

The following keys apply to the ensemble choice and options for the Mu-PT ensemble.

GCMC

Ensemble

Type:	Multiple Choice
Default value:	Mu-VT
Options:	[Mu-VT, Mu-PT]
Description:	Select the MC ensemble: Mu-VT for fixed volume or Mu-PT for variable volume. When the Mu-PT ensemble is selected the [Pressure] and [VolumeChangeMax] should also be specified.

VolumeChangeMax

Type:	Float
Default value:	0.05
Description:	Fractional value by which logarithm of the volume is allowed to change at each step. The new volume is then calculated as Vnew = exp(random(-1:1)*VolumeChangeMax)*Vold

Pressure

Type:	Float
Default value:	0.0
Unit:	Pascal
Description:	Pressure used to calculate the energy correction in the Mu-PT ensemble. Set it to zero for incompressible solid systems unless at very high pressures.

The GCMC code can insert multiple atom/molecule types in a single simulation, so it needs to keep track of what atom belongs to which insert. This information is automatically stored and updated when insertion/deletion/moving of atoms or molecules during the simulation, but is by default unknown for the atoms of the starting geometry. The GCMC code will therefore by default not modify the atoms in the original input in the MC trial moves. The Restart key and the Removables block are two ways to provide information about Deletable/Movable atoms/molecules in the input structure. If the Restart key is present the Removables block will be ignored.

GCMC

Restart

Type:

String

Description:

Name of an RKF restart file. Upon restart, the information about the GCMC input parameters, the initial system (atomic coordinates, lattice, charge, etc.) and the MC molecules (both already inserted and to be inserted) are read from the restart file. The global GCMC input parameters and the MC Molecules can be modified from input. Any parameter not specified in the input will use its value from the restart file (i.e. not the default value). Molecules found in the restart file do not have to be present as named Systems in the input, however if there is a System present that matches the name of a molecule from restart then the System’s geometry will replace that found in the restart file. It is also possible to specify new Molecules in the input, which will be added to the pool of the MC molecules from restart.

Removables

Type:

Non-standard block

Description:

The Removables can be used to specify a list of molecules that can be removed or moved during this GCMC calculation. Molecules are specified one per line in the format following format: MoleculeName atom1 atom2 … The MoleculeName must match a name specified in one of the [Molecule] blocks. The atom indices refer to the whole input System and the number of atoms must match that in the specified Molecule. A suitable Removables block is written to the standard output after each accepted MC move. If you do so then you should also replace the initial atomic coordinates with the ones found in the same file. If a [Restart] key is present then the Removables block is ignored.

An example of the Removables block:

Removables
  Oatom 41
  O2  44 45
  Oatom 42
  Oatom 43
End

This example specifies that 5 atoms belong to 4 different GCMC molecules of two different types, Oatom and O2. Thus in addition to the main input System there should be at least two additional Systems defined, one called “Oatom” (containing one atom) and the other “O2” (containing two atoms). The first one was inserted three times (atoms 41, 42, and 43) and the second one was inserted once.

Finally there are more technical keywords:

GCMC

MapAtomsToOriginalCell

Type:	Bool
Default value:	Yes
Description:	Keeps the atom (mostly) in the original cell by mapping them back before the geometry optimizations.

Note that the GeometryOptimization block is also read by the GCMC task, and the settings used for the individual optimizations. The documentation for these keywords can be found in the Geometry Optimization section of this manual.

Output¶

In addition to the standard KF variables in the “History” section on ams.rkf such as “Coords” and “Energy”, the following GCMC-specific variables are also created for each accepted MC step:

• MCMove - integer index of the MC move.
• MCMoveType - string containing the type of the MC move.
• MCMolecule - string containing the name of the inserted/displaced/removed molecule.
• Converged - a Fortran logical value containing the convergence status of the given geometry.

Results of a GCMC calculation are stored in the GCMC section of the RKF file, in a number of variables. The following variables contain a summary of the MC statistics up to and including the latest step:

• NIterMCtried - the latest iteration number.
• NIterMCaccept - the number of accepted MC moves.
• NIterMCreject - the number of rejected MC moves.
• NMCacceptAdd - the number of accepted MC molecule insertions.
• NMCacceptRemove - the number of accepted MC molecule removals.
• NMCacceptMove - the number of accepted MC molecule moves.
• NMCacceptVolume - the number of accepted volume changes.
• NMCrejectAdd - the number of rejected MC molecule insertions.
• NMCrejectRemove - the number of rejected MC molecule removals.
• NMCrejectMove - the number of rejected MC molecule moves.
• NMCrejectVolume - the number of rejected volume changes.

The following variables (actually arrays of the size Iterations) in the GCMC section contain the detailed information about all MC moves in the current simulation. Only the first NIterMCtried elements of each array contain valid data.

• HistoryAccepted - MC move status value (1 - accepted, 0 - rejected, -1 - not done yet).
• HistoryAMSEnergy - the AMS energy (the \(E^{AMS}\) above).
• HistoryMCEnergy - the corrected MC energy (\(E^{MC} = E^{AMS} - \Sigma \mu^{MC}_i\), where \(\Sigma \mu^{MC}_i\) is the total chemical potential of all inserted molecules).
• HistoryVolume - the simulation box volume.
• HistoryMoveType - the MC move type index (0 - insert, 1 - delete, 2 - displace, 3 - change volume). The name of the move type with index i can be found in the MoveType(i) variable.
• HistoryMoleculeType - the inserted/deleted/displaced molecule type index. The name of the molecule type with index i can be found in the MoleculeName(i) variable.
• HistoryMoleculeIndex - the inserted/deleted/displaced molecule index within its type.

References