pwtools.thermo.Gibbs

class pwtools.thermo.Gibbs(T=None, P=None, etot=None, phdos=None, axes_flat=None, volfunc_ax=None, case=None, **kwds)[source]

Bases: object

Calculate thermodynamic properties on a T-P grid in the quasiharmonic approximation, given some variation grid of unit cell axes (axes_flat) and corresponding phonon DOS data for each grid point.

The axes_flat approach allows us to easily calculate per-axis thermal expansion rates.

We have 3 cases for how unit cell axis can be be varied.

case

axes_flat.shape

cell parameter grid

fitfunc

bulk modulus

1d

(N,) or (N,1)

ax0 = a (-> V) (e.g. cubic or rhombohedral)

G(V), thus any EOS will do

V*d^2G/dV^2

2d

(N,2)

(ax0,ax1) = (a,b), (a,c) or (b,c) (e.g. hexagonal: (a,c))

G(ax0, ax1)

not defined

3d

(N,3)

(ax0,ax1,ax2) (e.g. triclinic)

no default

not defined

Internally, we only deal with ax0 (1d), ax0 + ax1 (2d) and ax0 + ax1 + ax2 (3d), thus for 2d and 3d it doesn’t matter which cell axes is which, you just have to remember :)

Notes

Functions for fitting G(V), V(T) etc are defined in a dictionary attribute fitfunc. This is a dict with functions, which have the signature func(x,y). The functions get x-y type data and return a class instance. The instances must be a Spline-like object with a get_min() method. The __call__() method must accept a keyword der for calculating derivatives in some cases. See self._default_fit_* to get an idea: for 1d Spline or PolyFit1D, for 2d: PolyFit or Interpol2D. See also scipy.interpolate. Use set_fitfunc() to change. This can (and should!) be used to tune fitting methods. See examples below for how to change fit functions.

Here is a list of all valid fitfunc keys and what they do:

key

value

‘1d-G’

fit G(V), __call__(V,der=2) for B(V) = V*d^2G/dV^2

‘2d-G’

fit G(ax0,ax1)

‘3d-G’

fit G(ax0,ax1,ax2), no default, define and test youself, please

‘1d-ax’

fit V(ax0)

‘alpha’

fit x_opt(T), x=ax0,ax1,ax2,V, __call__(T,der=1) for alpha_x = 1/x * dx/dT

‘C’

fit G_opt(T), __call__(T,der=2) for Cp = -T * d^G_opt(T,P)/dT^2

Hints for which fitfuncs to choose: The defaults are pretty good, but nevertheless testing is mandatory! Good choices are PolyFit / PolyFit1D for all axes-related grids (1d-G, 2d-G, 1d-ax) since you won’t have too many points here. Then a fit is better than a Spline. Always use PolyFit(..., scale=True) or scale the x-y data manually before fitting to the same order of magnitude! For T grids, choose a very fine T-axis (e.g. T=linspace(.., num=300) and use a Spline. Needed to resolve details of alpha(T) and C(T).

The methods calc_F(), calc_H() and calc_G() return dicts with nd-arrays holding calculated thermodynamic properties. Naming convention for dict keys returned by methods: The keys (strings) mimic HDF5 path names, e.g. /group0/group1/array, thus the last name in the path is the name of the array (z in the examples below). All previous names define the grid and thus the dimension of the array. A group name starting with “#” is just a prefix. Below, na=len(a), etc.

key

ndim

shape

description

/a/z

1d

(na,)

z on 1d a-grid

/a/b/z

2d

(na,nb)

z on 2d a-b grid

/a/b/c/z

3d

(na,nb,nc)

z on 3d a-b-c grid

/a-b/z

1d

(na*nb,)

z on flattened a-b grid

/a-b-c/z

1d

(na*nb*nc,)

z on flattened a-b-c grid

/#opt/a/z

1d

(na,)

optimized z on a-grid

axes_flat : Usually, flat grids like “ax0-ax1” or “ax0-ax1-ax2” are created by nested loops [(ax0_i,ax1_i) for ax0_i in ax0 for ax1_i in ax1] and therefore have shape (nax0*nax1,2) or (nax0*nax1*nax2,3) . But that is not required. They can be completely unstructured (e.g. if points have been added later to the grid manually) – only fitfunc must be able to handle that.

Units

B,P

GPa

T

K

F,G

eV

ax{0,1,2}

Ang

Cv,Cp

R (8.4314 J/(mol*K))

alpha_{V,ax{0,1,2}}

1/K

Examples

See test/test_gibbs.py for worked examples using fake data. Really, go there and have a look. Now!

>>> from pwtools import mpl, crys, eos
>>> # cubic cell, 1d case
>>> volfunc_ax = lambda x: crys.volume_cc(np.array([[x[0]]*3 + [90]*3]))
>>> gibbs=Gibbs(axes_flat=..., etot=..., phdos=...,
...             T=linspace(5,2500,100), P=linspace(0,20,5),
...             volfunc_ax=volfunc_ax)
>>> # EOS fit for G(V), polynomial for heat capacity (default is Spline)
>>> gibbs.set_fitfunc('C', lambda x,y: num.PolyFit1D(x,y,deg=5,scale=True))
>>> gibbs.set_fitfunc('1d-G', lambda x,y: eos.EosFit(x,y))
>>> g = gibbs.calc_G(calc_all=True)
>>> # 1d case
>>> V = g['/ax0/V']; G=g['/T/P/ax0/G']; T=g['/T/T']
>>> # plot G(V) for all T and P=20 GPa
>>> plot(V, G[:,-1,:].T)
>>> # Cp(T) for all P
>>> plot(T,g['/#opt/T/P/Cp'])
>>> # 2d case plot G(ax0,ax1) for T=2500 K, P=0 GPa
>>> G=g['/T/P/ax0-ax1/G']
>>> d=mpl.Data2D(XY=axes_flat, zz=G[-1,0,:])
>>> fig,ax=mpl.fig_ax3d(); ax.scatter(d.xx,d.yy,d.zz); show()
__init__(T=None, P=None, etot=None, phdos=None, axes_flat=None, volfunc_ax=None, case=None, **kwds)[source]
Parameters:

  • T (1d array) – temperature [K]

  • P (1d array) – pressure [GPa]

  • etot (1d array, (axes_flat.shape[0],)) – Total energy [eV] for each axes_flat[i,…]

  • phdos (sequence of 2d arrays (axes_flat.shape[0],)) – Phonon dos arrays (nfreq, 2) for each axes grid point. axes_flat[i,…] -> phdos[i] = <2d array (nfreq,2)>, units see HarmonicThermo. For each array phdos[i][:,0] = freq, phdos[i][:,1] = dos. Those are passed to HarmonicThermo.

  • axes_flat (1d or 2d array) –

    Flattened cell axes variation grid (for example result of nested loop over axes to vary). Will be cast to shape (N,1) if 1d with shape (N,) .

    1d: (N,) or (N,1) -> vary one cell axis, i.e. cubic cell
    2d: (N,2) -> vary 2 (e.g. a and c for hexagonal)
    example: itertools.product(ax0, ax1)
    3d: (N,3) -> vary a,b,c (general triclinic)
    example: itertools.product(ax0, ax1, ax2)

  • volfunc_ax (callable) –

    calculate cell volume based on cell axes, V[i] = volfunc_ax(axes_flat[i,…]) where axes_flat[i,…]:

    1d: [a0]
    2d: [a0, a1]
    3d: [a0, a1, a2]

    with a0,a1,a2 the length of the unit cell axes.

  • case (str, optional) – ‘1d’, ‘2d’, ‘3d’ or None. If None then it will be determined from axes_flat.shape[1]. Can be used to evaluate “fake” 1d data: set case=’1d’ but let axes_flat be (N,2) or (N,3)

  • **kwds (keywords) – passed to HarmonicThermo and added here as self.<key>=<value>

Methods

calc_F([calc_all])

Free energy properties along T axis for each axes grid point (ax0,ax1,ax2) in self.axes_flat.

calc_G([ret, calc_all])

Gibbs free energy and related properties on T-P grid.

calc_H([calc_all])

Enthalpy H=E+P*V, B and related properties on P grid without zero-point contributions.

set_fitfunc(what, func)

Update dict with fitting functions: self.fitfunc[what] = func.