Thermodynamic integration between machine-learned force fields
TI from Fe3+ to Fe2+ (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: Fe3+ and Fe2+ using MLFFs. The procedure is as follows:
Preparation
Step 0: Obtain initial POSCARs for Fe3+ and Fe2+
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: Fe3+(aq) and Fe2+(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:
|
In this case, Fe3+ (the oxidised state Ox) and Fe2+ (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()
You can then plot the free energy against the lambda values and integrate to obtain the thermodynamic integration:
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.