CodePlexProject Hosting for Open Source Software

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();

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.

// 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();

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