Temperature estimation from the measurement of the frequencies of specific vibrational modes

Ab initio computed vibrational frequencies of a mineral, at different cell volumes, are used for temperature estimation of a body whose IR spectrum has been recorded, within the framework of the Quasi-Harmonic Approximation (QHA)

It is advised to follow the tutorials Basic Eos and Dealing with the frequencies to learn the basics of the package usage.

The first command assures that graphical output is plot directly in the notebook:

In [1]:
%matplotlib inline

Start the program in the notebook. The example refers to the case of pyrope, for which a folder containing all the required ab initio data (static energies and frequencies as function of the unit cell volume) must be present, together with an input.txt file:

In [6]:
%run bm3_thermal_2.py



File quick_start.txt found in the master folder
Input files in 'pyrope' folder

Input file

STATIC
pyrope_static.dat
VOLUME
volume.dat
FREQ
pyrope_freq_2.dat
EXP
experimental.txt
FITVOL
POLY
725. 773. 16  3
FU
4 20
SET
0 1 4 5 8 10 15 21
CP
0. 2 -2. -1 -0.5
ALPHA
0. 2 -2  -1 -0.5
END

-------- End of input file -------


Static BM3 EoS

Bulk Modulus: 173.95 (0.06) GPa
Kp:             4.30 (0.03)
V0:           756.1932 (0.00) A^3

Suppose that the chosen mode for the investigation is the n. 85 (this is one of the most intense mode in the IR spectrum). Check how its frequency changes with the unit cell volume. Since a polynomial fit has been chosen (from the input.txt file, with the keyword POLY under FITVOL), use the function check_poly:

In [26]:
check_poly(81)

To check for consistency of the ab initio data, compute the equation of state at 300K (3^rd order Birch-Murnaghan EoS), and compare the result with the experimental bulk modulus:

In [9]:
eos_temp(300, prt=False)

 ** BM3 fit **

EoS at the temperature of 300.00 K

Bulk Modulus: 159.97   (0.08) GPa
Kp:             4.67   (0.02)
V0:           767.6963 (0.004) A^3

Fix status: False
Fit of frequencies is active
Polynomial fit: degree  3
Volume range: 725.00 - 773.00, points=16

For subsequent computations it is advised to keep Kp fixed to the values optimized at 300K (Kp = 4.67)

In [20]:
set_fix(4.67)

The current settings can be printed at any time by using the info.show command:

In [21]:
info.show()

Current settings and results

Static data            ** min, max volumes: 725.7483, 772.6196; points: 8
Frequency volume range ** min, max volumes: 725.7483, 772.6196; points: 22
Selected freq. sets    ** min, max volumes: 725.7483, 772.6196; points: 8
Frequency sets:           [0, 1, 4, 5, 8, 10, 15, 21]

Fit of frequencies     ** type: poly, degree: 3
                          min, max volumes: 725.0000, 773.0000; points 16

** Static EoS (BM3) **
K0: 173.95 GPa, Kp: 4.30, V0: 756.1932 A^3

** BM3 EoS from the last computation, at the temperature of 300.00 K **
K0: 159.97 GPa, Kp: 4.67, V0: 767.6963 A^3
Kp fixed

No excluded modes

Kieffer model off

Within the limits of QHA, frequencies do change with T due to thermal expansion: as T increases, the cell volume also increases (generally) and frequencies change accordingly as an effect of the volume variation. Check the thermal expansion at 300K, by using the command thermal_exp_p; to have a description of this command, just type:

In [23]:
help(thermal_exp_p)

Help on function thermal_exp_p in module __main__:

thermal_exp_p(tt, pp, plot=False, exit=False, **kwargs)
    Thermal expansion at given temperature and pressure
    
    Args:
        tt:               temperature
        pp:               pressure
        plot (optional):  plots pressure vs T values (see help to
                          the thermal_exp_v function)
    **kwargs: 
        fix: if fix is provided, it controls (and overrides the setting 
             possibly chosen by set_fix) the optimization of kp in BM3; 
             if fix > 0.1, kp = fix and it is not optimized.
             
    Note: see help for the thermal_exp_v function
In [24]:
thermal_exp_p(300,0)

Thermal expansion: 3.24e-05 K^-1
Bulk modulus:        159.96 GPa
Pressure:              0.00 GPa
Volume:              767.6964 A^3

Compute the frequency of the mode n. 85 at 300K:

In [25]:
help(freq_poly_p)

Help on function freq_poly_p in module __main__:

freq_poly_p(ifr, tt=300.0, p0=0.0, plot=True, prt=True, **kwargs)
    Prints the frequency of a given mode at some temperature
    and pressure if a spline fitting method has been chosen
    
    Args:
        ifr: mode number (starting from 0)
        tt: temperature (K)
        pp: pressure (GPa)
        
    kwargs:
        fix (optional): Kp value fixed to *fix* if *fix* > 0.1
        
    A polynomial fitting must be active
In [27]:
freq_poly_p(81, 300, 0)

Frequency: 858.86 cm^-1
Pressure 0.00 GPa, temperature 300.00 K, volume 767.6964 A^3

So, the frequency of that mode, at T=300K (and P=0 GPa) is 858.86 cm^-1. To check the reliability of the code and the possible error on the temperature estimation due to the fits, use the function temperature_freq, whose description is:

In [30]:
help(temperature_freq)

Help on function temperature_freq in module __main__:

temperature_freq(ifr, freq, tmin, tmax, npt, pp, degree=2, **kwargs)
    Computes the temperature given the frequency of a normal mode, at a fixed
    pressure. A T range must be specified
    
    Args:
        ifr:               normal mode number
        freq:              value of the frequency
        tmin:              minimum value of T
        tmax:              maximum value of T
        npt:               number of T points in the range
        pp:                pressure
        degree (optional): degree of the polynomial fitting the P/freq
                           values from the pressure_freq_list function
        
    **kwargs: 
        fix: Kp fixed, if fix=Kp > 0.1
        
    Note: it is advised to keep Kp fixed by either specifying fix, or
          by using set_fix.
          For "noisy" modes, use polynomial fits (set_poly), or 
          a spline fit (set_spline) with a large smooth parameter.
In [29]:
temperature_freq(81, 858.86, 200, 400, 12, 0, degree=3)

Temperature: 299.93 K
Mode: 81, frequency: 858.86 cm^-1, Pressure:   0.00 GPa

By providing exactly the value of 858.96 cm^-1 as input to the function temperature_freq, the returned value for T is 299.93 which is very close to the actual value of 300K. We can conclude that the error on the temperature estimation (due to the fitting procedure only) is in this case less than 0.1K.

Suppose that the value of the measured frequency is 855.5 cm^-1; the temperature estimation is:

In [34]:
temperature_freq(81, 855.5, 200, 600, 12, 0, degree=3)

Temperature: 429.36 K
Mode: 81, frequency: 855.50 cm^-1, Pressure:   0.00 GPa