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?

  • $\begingroup$ 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). $\endgroup$
    – CW Tan
    Aug 6, 2023 at 19:01
  • $\begingroup$ @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? $\endgroup$ Aug 6, 2023 at 20:37

1 Answer 1


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
                    last_but_one_entry = float(line_elements[-2])
                except ValueError:
                    # If conversion to float fails, skip this line
                # 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)
  • 2
    $\begingroup$ This answer would benefit from an overall summary of what the code is doing beyond the few code comments explaining particular parts. $\endgroup$
    – Tyberius
    Aug 9, 2023 at 13:06

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .