Thermodynamic integration between machine-learned force fields

Calculating the free energy difference [math]\displaystyle{ \Delta A }[/math] between two species can be difficult, especially if the structural changes between the species are significant. For this, we can use thermodynamic integration (TI), integrating over a coupling parameter λ, designating one species as the interacting system ([math]\displaystyle{ \lambda = 1 }[/math]) and the other species as the non-interacting system ([math]\displaystyle{ \lambda = 0 }[/math]):

[math]\displaystyle{ \Delta A = \int_0^1 \langle U_1 - U_0 \rangle_{\lambda} d\lambda }[/math]

where U₁ and U₀ are the potential energies of the interacting and non-interacting systems, respectively.

Introduction

Theory

In this how-to, we will go over the example of an Fe²⁺/Fe³⁺ electrochemical half-cell. A molecular dynamics (MD) calculation is performed in parallel using TI. Since this requires many steps, we use two machine-learned force fields (MLFFs) to speed this up add link to training mlffs page. We can check the accuracy of this later by performing another TI from the MLFF to a DFT functional.

During the Fe²⁺/Fe³⁺ half-cell reaction, one electron is transferred from the reservoir. The Helmholtz free energy difference ΔA is:

[math]\displaystyle{ \Delta A = \int_0^1 \langle U_1 - U_0 \rangle_{\lambda} d\lambda - \mu n }[/math]

where [math]\displaystyle{ \langle U_1 - U_0 \rangle }[/math] is the potential energy difference between an interacting system ([math]\displaystyle{ \lambda=1 }[/math]) and non-interacting system ([math]\displaystyle{ \lambda=0 }[/math]) taken from thermodynamic integration with coupling parameter [math]\displaystyle{ \lambda }[/math], [math]\displaystyle{ \mu }[/math] is the reference chemical potential, and [math]\displaystyle{ n }[/math] is the number of electrons involved in the reaction.

There is an accompanying how-to on calculating the reference chemical potential. In this how-to, you will learn how to calculate the remaining [math]\displaystyle{ \int_0^1 \langle U_1 - U_0 \rangle_{\lambda} d\lambda }[/math] term, then add this to [math]\displaystyle{ \mu }[/math] to calculate the free energy difference [math]\displaystyle{ \Delta A }[/math].

Calculation

Two MD calculations need to be run in parallel between two different systems: Fe³⁺ in 64 H₂O (Fe3P_64H2O) and Fe²⁺ in 64 H₂O (Fe2P_64H2O). I.e., Fe³⁺_(aq), Fe²⁺_(aq), and the water slab.

Procedure and input

Procedure

The procedure is split into two parts, the calculation part:

0. Obtain an initial structure
1. Prepare the directory structure for the TI
2. TI calculation

and the post-processing:

3. Extracting the energies
4. Integrate to obtain the free energy

Input files

The POSCAR files can be found in the calculating the electrochemical potential how-to page.

The Gamma-point only is used for the KPOINTS file:

Gamma-point only
 0
Monkhorst Pack
 1 1 1
 0 0 0

The INCAR files are provided in the text and discussed there. Finally, standard POTCAR files are used: PAW_PBE H 15Jun2001, PAW_PBE O 08Apr2002, and PAW_PBE Fe_sv 23Jul2007.

TI from Fe³⁺ to Fe²⁺ (MLFF)

Once the redox level has been calculated, you can perform the thermodynamic integration (TI) calculations. The first step is thermodynamic integration between two species: Fe³⁺ and Fe²⁺ using MLFFs. The procedure is as follows:

Preparation

Step 0: Obtain initial POSCARs for Fe³⁺ and Fe²⁺

A starting structure for the MD simulations in TI should be carefully chosen. In Ref. , a homemade MD simulation program was used to anneal the two systems: Fe³⁺_(aq) and Fe²⁺_(aq) from 1000 K to 400 K in a 1 ns NVT ensemble MD simulation. We begin with the final structure from each of those two simulations.

Tip: You can learn about running an MD simulation in VASP in the MD calculations page or the first part of the MD tutorials.

Step 1: Preparing the directories

To perform TI, you need to use the VCAIMAGES tag. This requires a parent directory from which the TI is run, containing the following INCAR file:

#Thermodynamic integration
VCAIMAGES = 0.0
NCORE_IN_IMAGE1 = 12

and an ML_FF file (since you are using MLFFs). The VCAIMAGES tag runs calculations in two image directories 01 and 02, which contain the two "non-interacting" λ=0 amd "interacting" λ=1 systems, respectively.

Important: The calculation will crash with a warning message if no ML_FF file is present. The file itself can be empty. NCORE_IN_IMAGE1 sets the number of cores used by the first image; the second image uses the remaining cores.

In this case, Fe³⁺ (the oxidised state Ox) and Fe²⁺ (the reduced state Red). Since MLFFs are used, they must contained the ML_FF's trained for the Ox and Red systems, respectively.

Tip: You can link the ML_FF files to the directory where they were refitted using: ln -s ${ML_FFN_directory}/ML_FFN ML_FF

Each 01 and 02 directory must also contain the same POSCAR file, the corresponding POTCAR file, a KPOINTS file, and an INCAR file. The two calculations will be run in these directories.

Calculation

Step 2: Molecular dynamics calculation

The VCAIMAGES runs the 01 and 02 calculations in parallel, with the number of cores for the first image given by NCORE_IN_IMAGE1 in the parent directory. The following INCAR can be used for the MLFFs:

#Molecular dynamics
IBRION = 0
ISYM   = 0
NSW    = 100000
POTIM  = 0.7
TEBEG = 300
TEEND = 500
RANDOM_SEED =         248489752                0                0
MDALGO = 2
ISIF = 2
SMASS = 1.0
POMASS = 8.0 16.0 55.847

#Machine learning
ML_LMLFF = .TRUE.        
ML_MODE = run
ML_MB = 5000
ML_DESC_TYPE = 1

#Thermodynamic integration
VCAIMAGES = 0.0

These will run two parallel MD calculations for the value of λ defined in VCAIMAGES, i.e., 0.0 in this example.

Important: We have used the Nosé-Hoover thermostat (MDALGO = 2). The trajectories in 01 and 02 directories need to be identical, so you cannot use the Langevin thermostat as it introduces random numbers.

You will then need to repeat these calculations for other values of λ, e.g., 0.25, 0.5, 0.75, 1.0. We recommend doing this in 5 separate directories, e.g.:

lambda_0p0  lambda_0p25  lambda_0p5  lambda_0p75  lambda_1p0

for 0.0, 0.25, 0.5, 0.75, and 1.0, respectively. The calculation should be submitted from the parent directory, similar to for NEB

Post-processing

Step 3: Extracting the energies

Each of these calculations will output the energy for each MD step in the following format:

free  energy ML TOTEN  =      -953.27166392 eV

You should take the values for the 01 and 02 directories separately, e.g., by grepping for the energies and temperatures in each of the sub-directories:

grep "free  energy" 01/OUTCAR | awk '{print $5}' > free_E.dat
grep "free  energy ML TOTEN  =" 02/OUTCAR | awk '{print $6}' >> free_E.dat
grep temperature 01/OUTCAR | awk '{print $6}' > T.dat

then calculate the two ensemble averages, before taking the difference between the two. We exclude the first 20000 MD steps to allow time for equilibration. You can do this with the following script, which will plot the probability vs potential energy, the MD step number vs the potential energy, and the MD step number vs the temperature:

from py4vasp import plot
import plotly.graph_objects as go
import numpy as np
import os 
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde

def delta_A(path, directory, lower, upper):
    number_str = directory.split("_")[1] 
    number_str = number_str.replace("p", ".")
    Lambda = float(number_str)
    print(path+directory)
    
    data = np.genfromtxt(path + str(directory) + "/free_E.dat", dtype=None, encoding=None)
    length=len(data)
    
    #U_ox, U_red = data[exclude:int(length/2)], data[int(length/2)+exclude:]
    #if "200000" in directory:
    #    upper = 200000
    U_ox, U_red = data[lower:upper], data[int(length/2)+lower:int(length/2)+upper]
    n_ox, n_red = range(lower, upper+1), range(int(length/2)+lower,int(length/2)+upper+1)
    print(n_ox, n_red)
    print(len(U_ox), len(U_red))
    U_red_av, U_ox_av = np.average(U_red), np.average(U_ox)
    #U_1_0 =  U_red - U_ox
    U_1_0_av = U_red_av - U_ox_av
    #U_1_0_av = np.average(U_1_0)
    return(Lambda, (U_1_0_av), (U_red-U_ox))

def diff_cutoff(path, lower, upper):
    files = [d for d in os.listdir(path) if d.startswith("lambda_")]

    lambdas, A = [], []
    for a in ['lambda_0p0', 'lambda_0p25', 'lambda_0p5', 'lambda_0p75', 'lambda_1p0']:
        print(files)
        if "200000" in a:
            upper = 200000
        temp1, temp2, data = delta_A(path, a, lower, upper)
        lambdas.append(temp1)
        A.append(temp2)
        print(lambdas, A)
        #plt.hist(data, bins=30, alpha= 0.5, label=a)
        U = data
        # Plot first graph
        kde = gaussian_kde(data)
        x = np.linspace(min(data), max(data), 20)
        ax1.plot(x, kde(x), '-', linewidth=2, alpha= 0.5, label=a)
        ax1.set_title('P vs. U')
        ax1.set_xlabel(r'$\Delta U_{ML}$')
        ax1.set_ylabel(r'$P(\Delta U_{ML})$')
        ax1.legend()       

        # Plot second graph
        ax2.set_title('MD step vs. U')
        ax2.set_xlabel(r'$\Delta U_{ML}$')
        ax2.set_ylabel('MD step')
        ax2.plot(U, range(lower,upper), '-', linewidth=2, alpha= 0.5, label=a)        

        # Plot T
        data = np.genfromtxt(path + str(a) + "/T.dat", dtype=None, encoding=None)
        #print(len(data))
        ax3.plot(data, range(len(data)), '-', linewidth=2, alpha= 0.5, label=a)
        ax3.set_title('MD step vs. T')
        ax3.set_xlabel('T')
        ax3.set_ylabel('MD step')
        ax3.legend()  
        
    return(lambdas, A)

path = "$PATH_TO_TI_MLFF_MLFF_DIRECTORIES/"

l_cutoff, A_cutoff = [], []

# Create a figure with 1 row and 2 columns
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 4))  # 1 row, 2 columns

lambdas, A = diff_cutoff(path, 20000, 100000)
l_cutoff.append(lambdas)
A_cutoff.append(A)

print((A_cutoff))
print((l_cutoff))

ax1.legend()
ax2.legend()

# Adjust layout so titles/labels don't overlap
plt.tight_layout()
plt.show()

File:mlff mlff test.png

Figure 4. Probability vs potential energy (cf. Supplementary Figure 8), the MD step number vs the potential energy, and the MD step number vs the temperature for between Fe²⁺ (λ = 0) and Fe³⁺ (λ = 1).

You can then plot the free energy against the lambda values and integrate to obtain the thermodynamic integration:

File:mlff mlff int.png

Figure 5. Free energy difference ΔA for thermodynamic integration using a parameter λ between Fe³⁺ (λ = 0) and Fe²⁺ (λ = 1) (cf. Supplementary Figure 8).

Step 4: Integrate to obtain the free energy

With the free energy for each individual λ, you can integrate over them to obtain the free energy of the TI.

[math]\displaystyle{ \Delta A = \int_0^1 \langle U_1 - U_0 \rangle_{\lambda} d\lambda - ne\Delta \bar{\phi} }[/math]

from scipy.integrate import simpson
print('Simpson: ' + str(simpson(A_cutoff, l_cutoff)[0]))

This gave a value of: -1.298 eV. Adding this to the value of [math]\displaystyle{ e \Delta \bar{\phi} }[/math] from earlier gives and considering that the redox potential [math]\displaystyle{ U_{redox} = -\Delta A /e }[/math], [math]\displaystyle{ U_{redox} }[/math] can be calculated as::

[math]\displaystyle{ U_{redox} = -\Delta A/e = -(-1.30 - (3.65))/1 = 4.97 \: \mathrm{ eV} }[/math].

Comparing this to the literature value (cf. Supplementary Table 6) of 4.99 eV, we are in reasonable agreement despite additional approximations.