Equation of state

The volumetric dependency of the zero-temperature ground state energies is usually described through an equation of state, represented by the interface IStateEquation in the namespace Atomic.Thermodynamics.StateEquations. The usual C# naming convension is followed here, where specialization modifiers are added in front of the name, IStateEquation is a special kind of equation, IFittedStateEquation is such an equation being fitted to data points. The "I" prefix denotes an interface, i.e. a description of behavior with no actual implementation. This interface has the following methods and properties (a property in C# is just like a variable or parameters):

public interface IStateEquation
{
    double Energy(double volume);
    double Pressure(double volume);
    double BulkModulus(double volume);
    double EquilibriumEnergy { get; }
    double EquilibriumVolume { get; }
}

An equation of state may be fitted from a finite set of data points by the methods BirchMurnaghan4.Fit or BirchMurnaghan5.Fit for 4/5-parameters Birch-Murnaghan models. In this example ground state energies are loaded as VASP results from a directory and an equation of state is fitted:

// Load volume and energy from VASP output files, one from each subdirectory.
List<IVolumeEnergyPoint> points = new List<IVolumeEnergyPoint>();
foreach (DirectoryInfo d in new DirectoryInfo("/home/mbjba/vasp")
    .GetDirectories())
{
    VaspResult r = new VaspResult(d);
    points.Add(new VolumeEnergyPoint(r.FinalStructure.Volume, r.Energy));
}

// Fit a BM4 equation of state (EOS) to the energies.
IStateEquation bm = BirchMurnaghan4.Fit(points);

// Plot fitted EOS using the method Energy defined by IStateEquation.
new Gnuplot()
    .SetXLabel("Volume")
    .SetYLabel("Energy")
    .Plot(
        new DataPlot(points.Select(p => new PlaneVector(p.Volume, p.Energy)),
            LineColor.Black, PointType.Square),
        new FunctionPlot(v => bm.Energy(v),
            points.Min(p => p.Volume), points.Max(p => p.Volume),
            LineColor.Black)
    )
    .Pause()
    .Run();

equation_of_state.png

Potentials

Thermodynamics is all about potentials. The two most important potentials, the Gibbs and the Helmholtz potentials are represented through corresponding C# interfaces, IGibbsPotential and IHelmholtzPotential, for systems in thermal equilibrium with constant pressure and constant volume being imposed, respectively.

potential_cube.png

The interface are idealizations of actual computations. The are defined for every possible continuously defined pressure/volume and temperature. The design philosophy is to provide a general description of what a thermodynamic potential is in general, without specifying the computational method of energies and other derived quantities. This way the same abstract description applies to phonon models as well as simpler Debye models and the given C# object may be passed freely around and consumed by other classes designed to take a potential as an argument. This is the interface definitions:

public interface IHelmholtzPotential
{
    double FreeEnergy(double temperature, double volume);
    double Pressure(double temperature, double volume);
    double BulkModulus(double temperature, double volume);
    double EquilibriumVolume(double temperature, double pressure);
    IGibbsPotential GibbsPotential { get; }
}

public interface IGibbsPotential
{
    double FreeEnergy(double temperature, double pressure);
    double Volume(double temperature, double pressure);
    IHelmholtzPotential HelmholtzPotential { get; }
}

The property GibbsPotential of IHelmholtzPotential is the Gibbs potential obtained by the Legendre transform.

For systems where also the number of particles is allowed to vary, IGibbsGrandPotential or IHelmholtzGrandPotential should be used instead.

Debye model

The Debye model is one instance of a Helmholtz potential. It assumes constant sound velocity along each polarization directions and is usually determined by the stiffness of the material as measured by the bulk modulus. The changing in stiffness by volume is expanded by the Debye-Grüneisen correction, determined by the derivative of the bulk modulus.

// Load again. This time using powerful LINQ expression. Store as VASP objects.
VaspResult[] results = new DirectoryInfo("/home/mbjba/vasp")
    .GetDirectories()
    .Select(d => new VaspResult(d))
    .ToArray();

IVolumeEnergyPoint[] points = results
    .Select(r => new VolumeEnergyPoint(r.FinalStructure.Volume, r.Energy))
    .ToArray();

// Fit a BM4 equation of state (EOS) to the energies.
BirchMurnaghan4 bm = BirchMurnaghan4.Fit(points);

// Extract a structure object to calculate mass and number of atoms. Just take the first.
Structure structure = results.First().FinalStructure;

double scaling = 0.617;
double gamma = (1.0 + bm.EquilibriumBulkModulusPressureDerivative) / 2.0 - 2.0 / 3.0;
double mass = DebyeGruneisenPotential.StructureMassGeometric(structure);
int atoms = structure.Sites.Count;

// The everything is required to fit the Debye model is calculated.
IHelmholtzPotential debye = new DebyeGruneisenPotential(bm.EquilibriumBulkModulus, bm.EquilibriumVolume, mass, atoms, scaling, gamma);

// The total energy is elastic energy and vibrational free energy.
IHelmholtzPotential pot = HelmholtzPotential.Add(HelmholtzPotential.Extend(bm), debye);

// Plot with magnetic moment on the secondary axis.
double t = 600.0;
new Gnuplot()
    .SetXLabel("Volume")
    .SetYLabel("Energy")
    .SetYSecondaryLabel("Magnetic moment")
    .Plot(
        new DataPlot(points.Select(p => new PlaneVector(p.Volume, p.Energy)),
            LineColor.Black, PointType.Square),
        new DataPlot(results.Select(r => new PlaneVector(r.FinalStructure.Volume, r.MagneticMoment)),
            Axes.XYSecondary, LineColor.Blue, PointType.Circle),
        new FunctionPlot(v => bm.Energy(v),
            points.Min(p => p.Volume), points.Max(p => p.Volume),
            LineColor.Black, new PlotLabel("E_0")),
        new FunctionPlot(v => pot.FreeEnergy(t, v),
            points.Min(p => p.Volume), points.Max(p => p.Volume),
            LineColor.Red, new PlotLabel("F_ph"))
    )
    .Pause()
    .Run();

debye_model.png

Temperature-only potential (electronic free energy)

There is also an even simpler interface IPotential for systems without a well-defined volume or where the volume is unimportant. This is the case for electronic potential defined from the electronic density of states as computed by VASP. In this case the potential is only defined for one fixed volume, and the free energy is only a function of temperature. The IPotential interfaces requires the following methods to be defined:

public interface IPotential
{
    FreeEnergy(double temperature);
}

The class ElectronicPotential implements this interface and provides a precise and efficient implementation of the required integration of the density of states.

If multiple electronic potentials have been calculated for a range of volumes,
the free energies may be interpolated similarly to the way an equation of state is constructed. This may be through the static FitPotential method of BirchMurnaghan4 class. The interpolated result is then naturally class implementing IHelmholtzPotential, as for each temperature the same interpolation method of volumes may be performed.

Example of loading electronic potentials and do this interpolation and plotting through Gnuplot for two different temperatures:

xxx

plot.png

Exactly the same applies to phonon computations. Here the PhononPotential is used. However, in this case the phonon density of states must be computed, e.g. using Phonopy with the wrapper class provided in the library.

Last edited Jul 3, 2014 at 2:01 PM by bakkedal, version 13