ensemble (TI_AS)

This keyword is used to set up a nonequilibrium thermodynamic integration integrator with Adiabatic switching (AS) path. It calculates Gibbs free energy along an isothermal path. Please check [Cajahuaringa2022] for more details.


The parameters can be specified as follows:

ensemble ti_rs temp <temperature> tperiod <tau_temperature> <pressure_control> <pmin> <pmax> pperiod <tau_pressure> tswitch <switch_time> tequil <equilibrium_time>
  • <temperature>: The temperature of the AS simulation.

  • <pmin> and <pmax>: The pressure range of the AS simulation.

  • <pressure_control>: Please refer to MTTK ensemble for more information.

  • <tau_temperature> and <tau_pressure>: These parameters are optional. Please refer to MTTK ensemble for more information.

  • <equilibrium_time>: The number timesteps to equilibrate the system.

  • <switch_time>: The number timesteps to vary pressure.

If you do not specify <equilibrium_time> and <switch_time>, they will be automatically set in a 1:4 ratio.


ensemble ti_as temp 300 aniso 10 30 tswitch 10000

This command switch pressure from 10 to 30 GPa in 10000 timesteps.

Output file

This command will produce a csv file. The columns are pressure and volume per atom.

The following code can help you calculate the Gibbs free energy:

from pandas import read_csv
import matplotlib.pyplot as plt
import numpy as np
from ase.units import kB, GPa
import yaml
from scipy.integrate import cumtrapz

with open("ti_spring_300.yaml", "r") as f:
    y =  yaml.safe_load(f)

T0 = y["T"]
G0 = y["G"]

ti = read_csv("ti_as.csv")
n = int(len(ti)/2)
forward = ti[:n]
backward = ti[n:][::-1]
p = forward["p"]
V1 = forward["V"]
V2 = backward["V"]

w = (cumtrapz(V1,p,initial=0) + cumtrapz(V2,p,initial=0))*0.5

G = G0 + w
plt.plot(p/GPa, G, label="AS")
plt.xlabel("P (GPa)")
plt.ylabel("G (eV/atom)")