# Helmholtz vibrational energy from VASP calculations

I would like to calculate the Helmholtz vibrational energy. I have done the phonon calculations. My question is how to extract the vibrational energy from the output of the Vasp calculation. A source code (vaspkit) I used allows me to input the desired temperature and gives the following output:

Temperature (T): 300.0 K
Zero-point energy E_ZPE : 72.232 kcal/mol 3.132266 eV
Thermal correction to U(T): 85.703 kcal/mol 3.716441 eV
Thermal correction to H(T): 85.703 kcal/mol 3.716441 eV
Thermal correction to G(T): 66.293 kcal/mol 2.874750 eV
Entropy S : 270.703 J/(molK) 0.002806 eV/K
Entropy contribution TS : 81210.770 J/(mol) 0.841691 eV


Since F = U - TS, can I calculate the energy from this output?

• I'm guessing that you'd extract the bare "U" as the energy from a typical SCF calculation. Then you'd have to add the "Thermal correction to U(T)" to it. Then subtract the "TS" term (on the last line). Aug 6, 2023 at 19:01
• @CWTan Ok, so my system is a surface with adsorbed oxygen atoms. I froze the substrate and performed the frequency calculation. Should I use the same system (frozen substrate) to perform the SCF or the unfrozen one? Aug 6, 2023 at 20:37

    #!/bin/bash

from ase.thermochemistry import HarmonicThermo

def search_and_process_cm1_lines(file_path):
result_list = []
with open(file_path, 'r') as file:
for line in file:
if "cm-1" in line:
# Split the line by whitespace
line_elements = line.strip().split()
# Convert the last but one entry to a float
try:
last_but_one_entry = float(line_elements[-2])
except ValueError:
# If conversion to float fails, skip this line
continue
# Divide the value by 1000 and append to the result list
result_list.append(last_but_one_entry / 1000)

return result_list

# Example usage
file_path = "OUTCAR"
cm1_values = search_and_process_cm1_lines(file_path)
thermo = HarmonicThermo(vib_energies=cm1_values)

# get thermodynamical parameters at required temperature
s = thermo.get_entropy(temperature=300)
u = thermo.get_internal_energy(temperature=300)
f = thermo.get_helmholtz_energy(temperature=300)

print("Entropy S: ", s)
print("Internal Energy U: ", u)
print("Helmholtz Energy F: ", f)

• This answer would benefit from an overall summary of what the code is doing beyond the few code comments explaining particular parts.
– Tyberius
Aug 9, 2023 at 13:06