SolidThermodynamicIntegrationHook

class SolidThermodynamicIntegrationHook(configuration, spring_constants, total_steps, temperature=None, equilibration_steps=None, lambda_min=None, measurement_call_interval=None)

Post step hook for calculating the free energy of solid

Parameters:

• configuration (DistributedConfiguration) – The current configuration.

• spring_constants (PhysicalQuantity of type force per distance.) – The atomic spring constants used in the reference system. Can be one spring constant for all atoms, or an array of spring constants, one for each atom.

• total_steps (int) – The total number of simulation steps.

• temperature (float) – The temperature of the system

• equilibration_steps (int) – The number of steps to equilibrate at the reference state.

• lambda_min (float) – Minimum value of the switching parameter used during the simulation.

• measurement_call_interval (int) – The interval between calculating measurement data

callInterval()

Returns:

The call interval of this hook function.

Return type:

int

reference()

Returns the reference system.

Returns:

The reference system

Return type:

BaseReference

referenceFreeEnergy()

Returns the free energy of the reference system.

Returns:

The enegy of the reference system

Return type:

PhysicalQuantity of type energy

uniqueString()

Return a unique string representing the state of the object.

Usage Examples

In the example below, the free energy of two Co crystal phases is calculated at different temperatures. The result can then be used to determine the relative stability of the two phases.

For each phase and temperature, a number of steps are performed:

• Equilibration in the iso-baric ensemble at the given temperature and pressure. During this simulation the volume is sampled, and finally averaged.

• From the iso-baric equilibration, the force constant is calculated and the starting configuration is generated using the sampled volume.

• A production MD simulation in the canonical ensemble is performed. During this simulation a number of properties are sampled, such as the potential energy of both the actual system and the reference system.

The result of each production simulation is the absolute free energy at the given temperature and pressure.

Note that each simulation takes about 15 minutes to complete on a standard workstation. This means that simulations of the two structures at four different temperatures would take about 2 hours to run.

ti_example_workflow.hdf5

General

Thermodynamic integration [1] is a method used to calculate free energy differences between systems. In molecular simulations, it is often used to calculate the free energy difference between two states of a system, such as two different crystal phases of a solid. In that way, one can determine which phase is the most stable under the given conditions.

Thermodynamic integration involves running several simulations where a parameter connecting the systems is varied. This can be computationally expensive and a more efficient way is to utilize non-equilibrium simulations. In non-equilibrium thermodynamic integration [2] [3] the parameter connecting the system is changed in a non-equilibrium manner in a single simulation. From this simulation the irreversible work can be estimated, which in turn enables one to calculate the free energy difference.

By calculating the free energy difference of a system with that of a reference system, with known free energy, one can calculate the absolute free energy of the system. This method is implemented in the two hook function classes SolidThermodynamicIntegrationHook and LiquidThermodynamicIntegrationHook. SolidThermodynamicIntegrationHook is suitable for calculating the free energy of a solid and uses the Einstein crystal as a reference system [3]. When using LiquidThermodynamicIntegrationHook the reference system is an Uhlenbeck-Ford liquid [4].

Notes

• The Uhlenbeck-Ford model depends on two parameters. One that determines the length scale of the system, and one that determines the energy scale. Both of these parameters can usually be taken as their Lennard-Jones counterparts.

• The Einstein crystal depends on a force constant. In the supplied template this force constant is automatically calculated.

[1]

Daan Frenkel and Berend Smit. Understanding molecular simulation: from algorithms to applications. Elsevier, 2023.

[2]

Masakatsu Watanabe and William P. Reinhardt. Direct dynamical calculation of entropy and free energy by adiabatic switching. Phys. Rev. Lett., 65:3301–3304, Dec 1990. URL: https://link.aps.org/doi/10.1103/PhysRevLett.65.3301, doi:10.1103/PhysRevLett.65.3301.

[3] (1,2)

Sarath Menon, Yury Lysogorskiy, Jutta Rogal, and Ralf Drautz. Automated free-energy calculation from atomistic simulations. Phys. Rev. Mater., 5:103801, Oct 2021. URL: https://link.aps.org/doi/10.1103/PhysRevMaterials.5.103801, doi:10.1103/PhysRevMaterials.5.103801.

[4]

Rodolfo Paula Leite, Rodrigo Freitas, Rodolfo Azevedo, and Maurice de Koning. The Uhlenbeck-Ford model: Exact virial coefficients and application as a reference system in fluid-phase free-energy calculations. The Journal of Chemical Physics, 145(19):194101, 11 2016. URL: https://doi.org/10.1063/1.4967775, arXiv:https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4967775/13587547/194101\_1\_online.pdf, doi:10.1063/1.4967775.