Tutorials

Harmonic and Anharmonic Thermochemistry

This tutorial explains step by step how to use the panther package to calculate the thermodynamic functions of molecules and solids making use of the standard harmonic vibrational analysis and anharmonic vibrations in the independent mode approximation as explained in detail in references [1], [2], [3].

For each step some of the data will be explicitly printed to illustrate the underlying data structures, however in production runs such level of verbosity is not necessary.

The methanol molecule will be used as and example and VASP code will be used to perform the calculations, however since panther is interfaced with ASE any of the supported calculators can be used instead with appropriate modifications.

Structure relaxation

First of all the molecule needs to be relaxed and it is recommended to converge the forces below 1.0e-5 eV/A. Here the initial structure is read from the methanol.xyz file and the structure relaxation is performed using the LBFGS method implemented in ASE instead of the internal VASP optimizers.

import ase.io
from ase.calculators.vasp import Vasp
from ase.optimizers import LBFGS

meoh = ase.io.read('methanol.xyz')

calc = Vasp(
        prec='Accurate',
        gga='PE',
        lreal=False,
        ediff=1.0e-8,
        encut=600.0,
        nelmin=5,
        nsw=1,
        nelm=100,
        ediffg=-0.001,
        ismear=0,
        ibrion=-1,
        nfree=2,
        isym=0,
        lvdw=True,
        lcharg=False,
        lwave=False,
        istart=0,
        npar=2,
        ialgo=48,
        lplane=True,
        ispin=1,
)

meoh.set_calculator(calc)

optimizer = LBFGS(meoh, trajectory='relaxed.traj',
                  restart='lbfgs.pkl', logfile='optimizer.log')

optimizer.run(fmax=0.00001)

Hessian matrix

Having optimized the structure we will also need the hessian matrix. We will use internal VASP mode (IBRION=5) to generate the hessian using cartesian coordinate displacements, therefore we need to update the calculator’s parameters. After the hessian is calculated it is read from the OUTCAR symmetrized, converted to atomic units and saved as a numpy array for convenience

from panther.io import read_vasp_hessian

# adjust the calcualtor argument for hessian calculation
calc.set(ibrion=5, potim=0.02)

calc.calculate(meoh)

hessian = read_vasp_hessian('OUTCAR', symmetrize=True, convert2au=True, negative=True)
np.save('hessian', hessian)

Harmonic vibrations and thermochemistry

We can now use the hessian to calculate thermochemical functions in the harmonic oscillator approximation, starting by calcualting the frequencies and normal modes

from panther.vibrations import harmonic_vibrational_analysis

frequencies, normal_modes = harmonic_vibrational_analysis(hessian, meoh,
            proj_translations=True, proj_rotations=True, ascomplex=False)

The resulting frequencies are in atomic units and need to be converted to Joules and passed to Thermochemistry to calculate thermochemical functions

from scipy.constants import value, Planck
from panther.thermochemistry import Thermochemistry

vibenergies = Planck * frequencies.real * value('hartree-hertz relationship')
vibenergies = vibenergies[vibenergies > 0.0]

thermo = Thermochemistry(vibenergies, meoh, phase='gas', pointgroup='Cs')
thermo.summary(T=273.15, p=0.1)
================ THERMOCHEMISTRY =================

   @ T = 273.15 K  p =   0.10 MPa

--------------------------------------------------
Partition functions:
ln q                     :          23.802
    ln q_translational   :          15.574
    ln q_rotational      :           7.949
    ln q_vibrational     :           0.280
--------------------------------------------------
Enthalpy (H)             :         140.014  kJ/mol
    H translational      :           3.407  kJ/mol
    H rotational         :           3.407  kJ/mol
    H vibrational        :         130.930  kJ/mol
        @ 0 K (ZPVE)     :         129.733  kJ/mol
        @ 273.15 K       :           1.197  kJ/mol
        pV               :           2.271  kJ/mol
--------------------------------------------------------------------------
                                                               *T
Entropy (S)              :           0.2355 kJ/mol*K        64.3395 kJ/mol
    S translational      :           0.1503 kJ/mol*K        41.0476 kJ/mol
    S rotational         :           0.0786 kJ/mol*K        21.4591 kJ/mol
    S vibrational        :           0.0067 kJ/mol*K         1.8328 kJ/mol
--------------------------------------------------------------------------
U - T*S                  :          75.6749 kJ/mol
--------------------------------------------------
Electronic energy        :       -2918.9516 kJ/mol

Normal Mode Relaxation

In some cases it is advantegeous to refine the structure using displacements along normal modes of vibrations, such a functionality is provided through the NormalModeBFGS which is based on the Optimizer class from the ASE pakckage. The method requires an initial guess for the hessian matrix for which we’ll reuse the hessian calculated in one of the previous steps, however in this case the hessian should be in eV/Angstrom^2 units, therefore it needs to be converted.

from panther.nmrelaxation import NormalModeBFGS

from scipy.constants import angstrom, value

ang2bohr = angstrom / value('atomic unit of length')
ev2hartree = value('electron volt-hartree relationship')

hessian = hessian * (ang2bohr**2) /ev2hartree

# create the optimizer
optimizer = NormalModeBFGS(meoh, 'gas', hessian, logfile='optimizer.log',
                           trajectory='relaxed.traj', proj_translations=True,
                           proj_rotations=True)

# start the relaxation
optimizer.run(fmax=0.001)

Anharmonic Thermochemistry

Internal coordinate displacements

With frequencies and normal modes we can further generate a grid of displacements along each normal mode using internal coordinates to improve the sampling of the potential energy surface. This is done using the calculate_displacements function. The function returns a nested OrderedDict of structures as ase.Atoms objects with mode number and displacement sample number as keys. For example if npoints=4 is given as an argument there will be 8 structures per mode labeled with numbers 1, 2, 3, 4, -1, -2, -3, -4 signifying the direction and the magnitude of the displacement.

from panther.displacements import calculate_displacements

images, modeinfo = calculate_displacements(meoh, hessian, frequencies, normal_modes, npoints=4)

print(modeinfo.to_string())

           HOfreq  effective_mass displacement is_stretch vibration     P_stretch        P_bend     P_torsion  P_longrange
mode
0     3748.362703     1944.298846      3.05268       True      True  1.000538e+00  9.633947e-07  3.581071e-08          0.0
1     3033.988514     2003.111476      3.39309       True      True  1.000487e+00  1.444002e-03  1.358614e-08          0.0
2     2956.839029     2015.939983      3.43707       True      True  1.000163e+00  2.112104e-03  0.000000e+00          0.0
3     2897.901987     1886.235715      3.47185       True      True  1.002597e+00  4.553844e-05  0.000000e+00          0.0
4     1445.646111     1896.588719      2.45777      False      True  2.423247e-04  8.382427e-01  1.620089e-01          0.0
5     1430.743791     1910.050531      2.47054      False      True  5.828069e-07  7.696362e-01  2.314594e-01          0.0
6     1413.372913     2064.662087      2.48567      False      True  1.101886e-05  1.013066e+00  1.246737e-02          0.0
7     1320.779390     2344.971044      2.57133      False      True  2.648486e-03  8.819130e-01  1.163030e-01          0.0
8     1122.078049     2310.188376      2.78972      False      True  8.856000e-04  8.562025e-01  1.423684e-01          0.0
9     1043.856152     2998.235955      2.89236      False      True  4.590826e-01  4.698373e-01  7.957428e-02          0.0
10     999.428001     4154.928137      2.95595      False      True  5.647864e-01  3.944914e-01  4.924659e-02          0.0
11     276.606858     1950.430919      5.61876      False      True  5.101834e-06  1.465507e-04  1.002540e+00          0.0
12       0.000000     8792.819312          inf      False     False           NaN           inf           NaN          0.0
13       0.000000     4750.377947          inf       True     False           inf           NaN           NaN          0.0
14       0.000000     6312.514911          inf       True     False           inf           NaN           NaN          0.0
15       0.000000     5243.620927          inf       True     False           inf           NaN           NaN          0.0
16       0.000000     2177.259022          inf      False     False           NaN           inf           NaN          0.0
17       0.000000     8532.108968          inf       True     False           inf           NaN           NaN          0.0

The function also returns modeinfo DataFrame with additional characteristics of the mode such as displacement, is_stretch and effective_mass and components of the vibrational population analysis.

Calculating energies for the displaced structures

Per each displaced structure we can calculate the energy, in this example using the VASP calculator again in the single point calculation mode

from panther.pes import calculate_energies

# set the calculator in single point mode
calc.set(ibrion=-1)

energies = calculate_energies(images, calc, modes='all')

This will return a DataFrame with npoints * 2 energies per mode.

The energies are missing the equilibrium structure energy which can be easily set through

energies['E_0'] = meoh.get_potential_energy()

print(energies.to_string())

         E_-4       E_-3       E_-2       E_-1        E_0        E_1        E_2        E_3        E_4
0  -29.838210 -29.999174 -30.129845 -30.219158 -30.252801 -30.212111 -30.072575 -29.801815 -29.357050
1  -29.887274 -30.033740 -30.148793 -30.224943 -30.252801 -30.220603 -30.113614 -29.913317 -29.596341
2  -29.765209 -29.983697 -30.134831 -30.223555 -30.252801 -30.223573 -30.135003 -29.984312 -29.766713
3  -29.880739 -30.032824 -30.149891 -30.225671 -30.252801 -30.222657 -30.125177 -29.948635 -29.679366
4  -30.194580 -30.220226 -30.238397 -30.249217 -30.252801 -30.249248 -30.238656 -30.221115 -30.196709
5  -30.196107 -30.220941 -30.238652 -30.249266 -30.252801 -30.249268 -30.238667 -30.220985 -30.196211
6  -30.199391 -30.222460 -30.239171 -30.249354 -30.252801 -30.249264 -30.238446 -30.219995 -30.193484
7  -30.202508 -30.224242 -30.239987 -30.249567 -30.252801 -30.249516 -30.239548 -30.222731 -30.198903
8  -30.208386 -30.227840 -30.241715 -30.250030 -30.252801 -30.250030 -30.241714 -30.227847 -30.208412
9  -30.209316 -30.228681 -30.242227 -30.250194 -30.252801 -30.250259 -30.242766 -30.230511 -30.213676
10 -30.210619 -30.229477 -30.242607 -30.250294 -30.252801 -30.250378 -30.243261 -30.231673 -30.215822
11 -30.241961 -30.246534 -30.249960 -30.252081 -30.252801 -30.252072 -30.249935 -30.246486 -30.241889

Calculating the frequencies

Frequencies can now be calculated using the finite difference method implemented in panther.pes.differentiate() function and appended as a frequency column to the modeinfo. The returned vibs matrix contains four columns corresponding to derivatives calculated with the central formula using 2, 4, 6 and 8 points

from panther.pes import differentiate
from scipy.constants import value

dsp = modeinfo.loc[modeinfo['vibration'], 'displacement'].astype(float).values
vibs = differentiate(dsp, energies, order=2)

au2invcm = 0.01 * value('hartree-inverse meter relationship')
np.sqrt(vibs) * au2invcm

array([[ 3757.6949986 ,  3745.36173164,  3745.5117494 ,  3745.51978786],
       [ 3038.75202112,  3032.49074659,  3032.55620542,  3032.56159197],
       [ 2960.0679129 ,  2956.11388526,  2956.13869496,  2956.14205087],
       [ 2900.20373811,  2897.17093192,  2897.18620368,  2897.18888617],
       [ 1446.18374932,  1446.16517443,  1446.1824322 ,  1446.19041632],
       [ 1431.77027217,  1431.68206976,  1431.68116942,  1431.68294377],
       [ 1414.58204596,  1414.17987052,  1414.182009  ,  1414.18039006],
       [ 1321.14267911,  1321.22424463,  1321.24026442,  1321.2470148 ],
       [ 1122.70558461,  1122.64210276,  1122.63752157,  1122.63384929],
       [ 1043.87393741,  1043.78448965,  1043.78296923,  1043.78025961],
       [  999.41759187,   999.30803727,   999.31024481,   999.30971807],
       [  285.06040099,   285.79248818,   285.87505212,   285.9193547 ]])

# assign the frequencies fitted with 8 points to a frequency column
# in the modeinfo
modeinfo.loc[modeinfo['vibration'], 'frequency'] = (np.sqrt(vibs)*au2invcm)[:, 3]

Fitting the potentials

The last this is to fit the potential energy surfaces as 6th and 4th order polynomials

from panther.pes import fit_potentials

# fit the potentials on 6th and 4th order polynomials
c6o, c4o = fit_potentials(modeinfo, energies)

The two DataFrame objects c6o and c4o contain fitted polynomial coefficients for each mode. We can use the energies and the polynomial coefficients to plot the PES and the fitted potentials, here as an example, second mode (mode=1 since the modes are indexed from 0) is plotted

from panther.plotting import plotmode

plotmode(1, energies, modeinfo, c6o, c4o)
Plot of the mode potential

Anharmonic frequencies from 1-D Schrodinger Equation

Anharmonic frequencies are calculated first by solving the 1-D Schrodinger equation per mode as exaplained in reference [4] and then those frequencies are used to calculate the thermodynamic functions

from panther.anharmonicity import anharmonic_frequencies, harmonic_df, merge_vibs
from panther.thermochemistry import AnharmonicThermo

anh6o = anharmonic_frequencies(meoh, 273.15, c6o, modeinfo)
anh4o = anharmonic_frequencies(meoh, 273.15, c4o, modeinfo)

harmonicdf = harmonic_df(modeinfo, 273.15)
finaldf = merge_vibs(df6, df4, hdf, verbose=False)

at = AnharmonicThermo(fdf, meoh, phase='gas', pointgroup='Cs')
at.summary(T=273.15, p=0.1)
================ THERMOCHEMISTRY =================

  @ T = 273.15 K  p =   0.10 MPa

--------------------------------------------------
Partition functions:
ln q                     :          23.667
    ln qtranslational    :          15.574
    ln qrotational       :           7.949
    ln qvibrational      :           0.145
--------------------------------------------------
Enthalpy (H)             :         138.890  kJ/mol
    H translational      :           3.407  kJ/mol
    H rotational         :           3.407  kJ/mol
    H vibrational        :         129.806  kJ/mol
        @ 0 K (ZPVE)     :         129.675  kJ/mol
        @ 273.15 K       :           0.131  kJ/mol
        pV               :           2.271  kJ/mol
--------------------------------------------------------------------------
                                                               *T
Entropy (S)              :           0.2375 kJ/mol*K        64.8694 kJ/mol
    S translational      :           0.1503 kJ/mol*K        41.0476 kJ/mol
    S rotational         :           0.0786 kJ/mol*K        21.4591 kJ/mol
    S vibrational        :           0.0086 kJ/mol*K         2.3627 kJ/mol
--------------------------------------------------------------------------
H - T*S                  :          74.0209 kJ/mol
--------------------------------------------------
Electronic energy        :       -2918.9516 kJ/mol
[1]Piccini, G., Alessio, M., Sauer, J., Zhi, Y., Liu, Y., Kolvenbach, R., Jentys, A., Lercher, J. A. (2015). Accurate Adsorption Thermodynamics of Small Alkanes in Zeolites. Ab initio Theory and Experiment for H-Chabazite. The Journal of Physical Chemistry C, 119(11), 6128–6137. doi:10.1021/acs.jpcc.5b01739
[2]Piccini, G., & Sauer, J. (2014). Effect of anharmonicity on adsorption thermodynamics. Journal of Chemical Theory and Computation, 10, 2479–2487. doi:10.1021/ct500291x
[3]Piccini, G., & Sauer, J. (2013). Quantum Chemical Free Energies: Structure Optimization and Vibrational Frequencies in Normal Modes. Journal of Chemical Theory and Computation, 9(11), 5038–5045. doi:10.1021/ct4005504
[4]Beste, A. (2010). One-dimensional anharmonic oscillator: Quantum versus classical vibrational partition functions. Chemical Physics Letters, 493(1-3), 200–205. doi:10.1016/j.cplett.2010.05.036