mineralscloud / equationsofstateofsolids.jl Goto Github PK

View Code? Open in Web Editor NEW

13.0 2.0 0.0 4.83 MB

A Julia package for fitting the equation of state of solids, and more

Home Page: https://mineralscloud.github.io/EquationsOfStateOfSolids.jl/

License: MIT License

Julia 100.00%

equation-of-state solid-state-physics fitting-curve thermodynamic-properties solids julia-package solve strains

equationsofstateofsolids.jl's People

Contributors

Stargazers

Watchers

equationsofstateofsolids.jl's Issues

Register version v0.4.1

@JuliaRegistrator register

Fit EOS with `f` rather than `v`

Similar reason as #5, but this may benefit us because logarithmic EOS is very hard to fit if fitting by v not f. We can fit f first (which does not have log involved), and reverse back to find v.

Deprecate `LeastSquaresOptim`, use `LsqFit`?

See #11 (comment), it is hard to use

Deprecate `NTuple{N,T}` as `x0` in `BirchMurnaghan{N,T}`?

As mentioned in #24 and #25, the current implementation has some drawback:

e0 is not properly handled in nonlinfit
Base.getproperty will throw BoundsError
It is ugly to implement model inside nonlinfit.
It is ugly to write e0 = zero(eos.params.v0 * eos.params.b0).

Add `strainform(::EnergyEOS{<:Murnaghan})` to allow calculate from strains?

Interesting exponentiation

For Reals and Complexes, usually ^(2/3) is faster:

julia> @btime 1901^(2//3);
  29.008 ns (0 allocations: 0 bytes)

julia> @btime 1901^(2/3);
  1.225 ns (0 allocations: 0 bytes)

julia> @btime 1901im^(2/3);
  24.609 ns (0 allocations: 0 bytes)

julia> @btime 1901im^(2//3);
  36.309 ns (0 allocations: 0 bytes)

But for SymPy.Syms, ^(2//3) is faster:

julia> v, v0 = SymPy.symbols("v, v0")
(v, v0)

julia> v = SymPy.symbols("v")
v

julia> @btime v^(3/2)
  17.177 μs (15 allocations: 432 bytes)
 1.5
v

julia> @btime v^(3//2)
  7.786 μs (21 allocations: 576 bytes)
 3/2
v

How to get fields from `BirchMurnaghan`?

The current implementation

struct BirchMurnaghan{N,T} <: FiniteStrainEossParameters{N,T}
    x0::NTuple{N,T}
end

function Base.getproperty(value::FiniteStrainEossParameters, name::Symbol)
    if name == :v0
        return value.x0[1]
    elseif name == :b0
        return value.x0[2]
    elseif name == :b′0
        return value.x0[3]
    else
        return getfield(value, name)
    end
end

will throw an BoundsError in the following case

julia> BirchMurnaghan([1,2]).b′0
ERROR: BoundsError: attempt to access (1, 2)
  at index [3]

Register version v0.2.2

@JuliaRegistrator register

Deprecate #27, use the old way like `BirchMurnaghan3rd`

struct BirchMurnaghan3rd{T} <: FiniteStrainEOS{T}
    v0::T
    b0::T
    b′0::T
    e0::T
end

Use `AutoHashEquals.@auto_hash_equals` to check equality

If not applied on Parameters, we will have the following results:

julia> BirchMurnaghan3rd([1u"m^3",2u"GPa",3.0,4u"J"]) == BirchMurnaghan3rd([1u"m^3",2u"GPa",3.0,4u"J"])
true

julia> BirchMurnaghan3rd([1u"m^3",2u"GPa",big(3.0),4u"J"]) == BirchMurnaghan3rd([1u"m^3",2u"GPa",3.0,4u"J"])
false

julia> BirchMurnaghan3rd([1u"m^3",2u"GPa",big(3.0),4u"J"]) == BirchMurnaghan3rd([1u"m^3",2u"GPa",big(3.0),4u"J"])
false

Which is kind of weird, i.e., they cannot be equal if one element is a Big one. With @auto_hash_equals, this problem is solved:

julia> BirchMurnaghan3rd([1u"m^3",2u"GPa",big(3.0),4u"J"]) == BirchMurnaghan3rd([1u"m^3",2u"GPa",big(3.0),4u"J"])
true

Which is faster? Closure of type?

We have code

abstract type Parameters{T} end

abstract type ParametersFiniteStrain{N,T} <: Parameters{T} end

struct BirchMurnaghan{N,T} <: ParametersFiniteStrain{N,T}
    x0::NTuple{N,T}
end
function BirchMurnaghan{3}(v0, b0, b′0)
    T = Base.promote_typeof(v0, b0, b′0)
    return BirchMurnaghan{3,T}(Tuple(convert(T, x) for x in (v0, b0, b′0)))
end

Consider which will be faster?

The type approach

First, using a type for the real EquationOfStateOfSolids:

abstract type EquationOfStateOfSolids{T<:Parameters} end
struct EnergyEOSS{T} <: EquationOfStateOfSolids{T}
    params::T
end
struct PressureEOSS{T} <: EquationOfStateOfSolids{T}
    params::T
end
struct BulkModulusEOSS{T} <: EquationOfStateOfSolids{T}
    params::T
end

function (eos::EnergyEOSS{<:BirchMurnaghan3rd})(v)
    v0, b0, b′0 = eos.params.x0
    x = cbrt(v0 / v)
    y = x^2 - 1
    return 9 / 16 * b0 * v0 * y^2 * (6 - 4 * x^2 + b′0 * y)
end

To make users feel a clean code, we also want to define

energyeq(p::Parameters) = EnergyEOSS(p)
pressureeq(p::Parameters) = PressureEOSS(p)
bulkmoduluseq(p::Parameters) = BulkModulusEOSS(p)

Now

julia> x = BirchMurnaghan{3}(1, 20, 300)
BirchMurnaghan{3,Int64}((1, 20, 300))

julia> eos = energyeq(x)
EquationsOfStateOfSolids.Collections.EnergyEOSS{BirchMurnaghan{3,Int64}}(BirchMurnaghan{3,Int64}((1, 20, 300)))

julia> eos(3)
-460.13550846092517

julia> using BenchmarkTools

julia> @btime eos(3)
  42.890 ns (1 allocation: 16 bytes)
-460.13550846092517

julia> @btime $eos(3)
  34.429 ns (0 allocations: 0 bytes)
-460.13550846092517

julia> @code_lowered eos(3)
CodeInfo(
1 ─ %1  = Base.getproperty(eos, :params)
│   %2  = Base.getproperty(%1, :x0)
│   %3  = Base.indexed_iterate(%2, 1)
│         v0 = Core.getfield(%3, 1)
│         @_5 = Core.getfield(%3, 2)
│   %6  = Base.indexed_iterate(%2, 2, @_5)
│         b0 = Core.getfield(%6, 1)
│         @_5 = Core.getfield(%6, 2)
│   %9  = Base.indexed_iterate(%2, 3, @_5)
│         b′0 = Core.getfield(%9, 1)
│   %11 = v0 / v
│         x = EquationsOfStateOfSolids.Collections.cbrt(%11)
│   %13 = x
│   %14 = Core.apply_type(Base.Val, 2)
│   %15 = (%14)()
│   %16 = Base.literal_pow(EquationsOfStateOfSolids.Collections.:^, %13, %15)
│         y = %16 - 1
│   %18 = 9 / 16
│   %19 = b0
│   %20 = v0
│   %21 = y
│   %22 = Core.apply_type(Base.Val, 2)
│   %23 = (%22)()
│   %24 = Base.literal_pow(EquationsOfStateOfSolids.Collections.:^, %21, %23)
│   %25 = x
│   %26 = Core.apply_type(Base.Val, 2)
│   %27 = (%26)()
│   %28 = Base.literal_pow(EquationsOfStateOfSolids.Collections.:^, %25, %27)
│   %29 = 4 * %28
│   %30 = 6 - %29
│   %31 = b′0 * y
│   %32 = %30 + %31
│   %33 = %18 * %19 * %20 * %24 * %32
└──       return %33
)

The closure approach

This approach is much simpler, we just define

function energyeq(p)
    v0, b0, b′0 = p.x0
    function (v)
        x = cbrt(v0 / v)
        y = x^2 - 1
        return 9 / 16 * b0 * v0 * y^2 * (6 - 4 * x^2 + b′0 * y)
    end
end

and no need for EquationOfStateOfSolids to be defined. But will closures affect speed?

julia> x = BirchMurnaghan{3}(1, 20, 300)
BirchMurnaghan{3,Int64}((1, 20, 300))

julia> eos = energyeq(x)
#3 (generic function with 1 method)

julia> typeof(eos)
EquationsOfStateOfSolids.Collections.var"#3#4"{Int64,Int64,Int64}

It seems from here that the answer is clear: the only difference between the "type" approach and the "closure" approach is that the "type" approach officially defines several types (EquationOfStateOfSolids, etc.) but the other one just uses an internal name var"#3#4"{Int64,Int64,Int64}. I have closed and reopened the REPL multiple times and the results are the same. But we still want to perform a benchmark:

julia> @btime eos(3)
  43.101 ns (1 allocation: 16 bytes)
-460.13550846092517

julia> @btime $eos(3)
  34.427 ns (0 allocations: 0 bytes)
-460.13550846092517

julia> @code_lowered eos(3)
CodeInfo(
1 ─ %1  = Core.getfield(#self#, :v0)
│   %2  = %1 / v
│         x = EquationsOfStateOfSolids.Collections.cbrt(%2)
│   %4  = x
│   %5  = Core.apply_type(Base.Val, 2)
│   %6  = (%5)()
│   %7  = Base.literal_pow(EquationsOfStateOfSolids.Collections.:^, %4, %6)
│         y = %7 - 1
│   %9  = 9 / 16
│   %10 = Core.getfield(#self#, :b0)
│   %11 = Core.getfield(#self#, :v0)
│   %12 = y
│   %13 = Core.apply_type(Base.Val, 2)
│   %14 = (%13)()
│   %15 = Base.literal_pow(EquationsOfStateOfSolids.Collections.:^, %12, %14)
│   %16 = x
│   %17 = Core.apply_type(Base.Val, 2)
│   %18 = (%17)()
│   %19 = Base.literal_pow(EquationsOfStateOfSolids.Collections.:^, %16, %18)
│   %20 = 4 * %19
│   %21 = 6 - %20
│   %22 = Core.getfield(#self#, :b′0)
│   %23 = %22 * y
│   %24 = %21 + %23
│   %25 = %9 * %10 * %11 * %15 * %24
└──       return %25
)

Summary

By comparing their @code_lowered, the closure approach maybe even faster, because the type approach need to unwrap
everytime:

v0, b0, b′0 = eos.params.x0

But this is done only once for the closure. That's why it has a shorter @code_lowered.

Use `AbstractArray{T}` to replace manually `Base.promote_typeof`

The current way of making a Parameters of different types is

function BirchMurnaghan{3}(v0, b0, b′0)
    T = Base.promote_typeof(v0, b0, b′0)
    return BirchMurnaghan{3,T}(Tuple(convert(T, x) for x in (v0, b0, b′0)))
end
BirchMurnaghan{3}(arr::AbstractArray) = BirchMurnaghan{3}(arr...)

That is, manually Base.promote_typeof and convert. But an AbstractArray{T} automatically promotes them when constructing

julia> [1u"m^3", 3u"GPa", 4]
3-element Array{Quantity{Int64,D,U} where U where D,1}:
...

julia> Base.promote_typeof(1u"m^3", 3u"GPa", 4)
Quantity{Int64,D,U} where U where D

julia> [1..10, 2..5, 3]
3-element Array{Any,1}:
  Interval{Int64,Closed,Closed}(1, 10)
  Interval{Int64,Closed,Closed}(2, 5)
 3

julia> Base.promote_typeof(1..10, 2..5, 3)
Any

julia> [measurement("1 +- 0.1"), 3.0, 2, measurement("-123.4(56)")]
4-element Array{Measurement{Float64},1}:
    1.0 ± 0.1
    3.0 ± 0.0
    2.0 ± 0.0
 -123.4 ± 5.6

julia> Base.promote_typeof(measurement("1 +- 0.1"), 3.0, 2, measurement("-123.4(56)"))
Measurement{Float64}

So instead of the original implementation, we could write

BirchMurnaghan{3}(arr::AbstractArray{T}) where {T} = BirchMurnaghan{3,T}(NTuple{3,T}(arr))
BirchMurnaghan{3}(v0, b0, b′0) = BirchMurnaghan{3}([v0, b0, b′0])

and bingo!

Register v0.1.0

Register version v0.1.2

@JuliaRegistrator register

Check domains of `StrainToVolume`

For the 4 types of strains, usually they are reals, but sometimes, even if they are complex numbers, they still can produce real volumes.

`nonlinfit` return `e0` is not properly handled

In #23 I implemented nonlinfit on EnergyEoss. But e0 as a parameter of fit.param is directly fed into constructorof(T) (the BirchMurnaghan, etc.), which will result in an extra parameter for BirchMurnaghan, and cause it upgraded by order 1 (i.e., BirchMurnaghan{3} to BirchMurnaghan{4}). This is wrong.

However, if I write return (constructorof(T)(fit.param), fit.param[end]), fit, this will return a Tuple{BirchMurnaghan,Number} as the 1st element, which has different type of the result of nonlinfit for PressureEoss, which is just a BirchMurnaghan.

How to generate `e0`

Use zero(eltype(b0 * v0)).

Then methods like

BirchMurnaghan{3}(v0::Real, b0::Real, b′0::Real) = BirchMurnaghan{3}(v0, b0, b′0, 0)
BirchMurnaghan{3}(v0::AbstractQuantity, b0::AbstractQuantity, b′0) = BirchMurnaghan{3}(v0, b0, b′0, 0 * u"eV")

is not needed.

But we still need to make sure b0 has the unit of pressure when v0 is an AbstractQuantity.

Small change in floating point exponentiation causes a test failure on Julia 1.8.

From the PkgEval log at https://s3.amazonaws.com/julialang-reports/nanosoldier/pkgeval/by_hash/9b7990d_vs_8611a64/EquationsOfStateOfSolids.primary.log it shows this package has a test failure on the upcoming Julia 1.8 release.

I looked at it a little bit and it seems that a small precision change to floating point power is causing the problem. The test that errors is at

EquationsOfStateOfSolids.jl/test/Fitting.jl

Lines 650 to 654 in 3d9cc31

 nonlinfit( 

 EnergyEquation(Murnaghan1st(19u"angstrom^3", 100u"GPa", 4)), 

 volumes, 

 energies, 

 ),

If we print out v0_prev and v0 we can see that it just oscillates back and forth:

(v0_prev, v0) = (19.361802843683524, 19.36180284371997)
(v0_prev, v0) = (19.36180284371997, 19.361802843683524)
(v0_prev, v0) = (19.361802843683524, 19.36180284371997)
(v0_prev, v0) = (19.36180284371997, 19.361802843683524)
(v0_prev, v0) = (19.361802843683524, 19.36180284371997)
(v0_prev, v0) = (19.36180284371997, 19.361802843683524)
(v0_prev, v0) = (19.361802843683524, 19.36180284371997)
(v0_prev, v0) = (19.36180284371997, 19.361802843683524)
(v0_prev, v0) = (19.361802843683524, 19.36180284371997)
(v0_prev, v0) = (19.36180284371997, 19.361802843683524)

and the convergence criteria is never satisfied. So it seems this algorithm used here need to be a bit more robust to handle cases like this.

The difference in the result comes from:

S = EquationsOfStateOfSolids.FiniteStrains.EulerianStrain
volumes = ([17.789658, 18.382125, 18.987603, 19.336585, 19.606232, 20.238155, 20.883512])
v0 = 19.361802843683524
strains = map(EquationsOfStateOfSolids.FiniteStrains.ToStrain{S}(v0), volumes)

On 1.7:

julia> strains[1]
0.029040354449144212

On 1.8:

 julia> strains[1]
0.029040354449144323

Use simplified EOS with variable `f` (strain) instead of `x = v0 / v` (volume ratio)

The first one usually is simpler and maybe faster.

Register version v0.2.1

@JuliaRegistrator register

Use `UnPack.@unpack` to unpack `EossParameters`?

I used to use IterTools.fieldvalues to unpack EOS parameters, but it seems nice to use UnPack.@unpack? It looks nice and just does one thing and does it well.

Add `v2p`

will be useful in QuasiHarmonicApprox.jl.

Register version v0.2.0

@JuliaRegistrator register

Fitting not consider units when meeting `Murnaghan1st{Number}`

Sometimes the EOSs we read in have Number as eltype,

Murnaghan1st{Number}
 v0 = 300.44 a₀³
 b0 = 74.88 GPa
 b′0 = 4.82
 e0 = -612.43149513 Ry

When fitting with this kind of EOSs, method _unormal(p::Parameters{<:AbstractQuantity}) will not be called, thus without unit conversions, could result in wrong results. So I have to deprecate the methods for unit detection, just do it for all EOSs.

How to fit with `nonlinfit` without type `EquationOfStateOfSolids`?

Use the length of parameters to determine which order of the equation is

For example, BirchMurnaghan(1, 20, 300) would be a third-order EOS.

Set `e0` to be the minimum energy when fitting

Set e0 to be the minimum energy when fitting if eos.param.e0 is not provided (by default, zero).

Add boilerplate outer constrcutors to `EossParam`s

In the past I used to add outer constructors using eval-quote to types that are concrete subtypes of a type, that was slow to load and a bit o f ugly. Or I can list all concrete subtypes by hand, it's maybe a bit faster but it requires manual updating.

Make `FromStrain` do not throw errors

I often come across DomainErrors when find_zero of equations of state using Roots.jl, in the current code. I guess it's just that their root-finding algorithms has to search around and it goes beyond f < - 1 / 2. So I am going to remove this restriction.

Register version 0.1.1

@JuliaRegistrator register

Use `e0` as a separate field together with `x0::NTuple` in `BirchMurnaghan`

In #26 I said I wanted to deprecate NTuple as x0 in BirchMurnaghan{N,T}, but there
is still 1 possibility left:

struct BirchMurnaghan{N,T} <: FiniteStrainEossParameters{N,T}
    x0::NTuple{N,T}
    e0::T
end

Let's see how far can we go.

Construct `BirchMurnaghan{N}` by `BirchMurnaghan(N)`

For example, BirchMurnaghan(3) is a function that would produce BirchMurnaghan{3} and BirchMurnaghan{3} is a constructor now.

BirchMurnaghan(N::Integer) = BirchMurnaghan{N}

Just like how I used to construct Step.

Add `e0 = zero(eltype(v0 * b0))` in `energyeos` closure

Since we deprecate EquationOfStateOfSolids in #3, we cannot just add e0 = 0 to the closure. Since what if the EOS has units? But this can be solved by using e0 = zero(eltype(v0 * b0)), as mentioned in #9.

function energyeos(p::BirchMurnaghan{3})
    v0, b0, b′0 = p.x0
    function (v, e0 = zero(eltype(v0 * b0)))
        x = cbrt(v0 / v)
        y = x^2 - 1
        return 9 / 16 * b0 * v0 * y^2 * (6 - 4 * x^2 + b′0 * y)
    end
end

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Should I restrict the element type of `BirchMurnaghan` to `Number`?

No. Then we cannot use Symbols and Intervals in it.

Will it be faster? If I take `_model` out of `nonlinfit` (separate kernel functions)

I.e., write

function nonlinfit(eos::EquationOfStateOfSolids{T}, xs, ys; kwargs...) where {T}
    model = _model(eos)
    fit =
        curve_fit(model, float.(xs), float.(ys), float.(fieldvalues(eos.param)); kwargs...)
    if fit.converged
        return constructorof(T)(fit.param), fit
    else
        @error "fitting is not converged! Check your initial parameters!"
        return nothing, fit
    end
end # function nonlinfit

function _model(::S) where {T,S<:EquationOfStateOfSolids{T}}
    constructor = constructorof(S) ∘ constructorof(T)
    return (x, p) -> map(constructor(p), x)
end

instead of

function nonlinfit(eos::S, xs, ys; kwargs...) where {T,S<:EquationOfStateOfSolids{T}}
    constructor = constructorof(S) ∘ constructorof(T)
    model(x, p) = map(constructor(p), x)
    fit =
        curve_fit(model, float.(xs), float.(ys), float.(fieldvalues(eos.param)); kwargs...)
    if fit.converged
        return constructorof(T)(fit.param), fit
    else
        @error "fitting is not converged! Check your initial parameters!"
        return nothing, fit
    end
end # function nonlinfit

Will `strain_from_volume` speedup calculations? If so, where to put it?

Which of the following will run faster (when applied to a volume)?

const BirchMurnaghan3rdStrain = strain_from_volume(Eulerian())

function pressureeos(p::BirchMurnaghan{3})
    v0, b0, b′0 = p.x0
    function (v)
        f = (cbrt(v0 / v)^2 - 1) / 2
        return 3f / 2 * b0 * sqrt(2f + 1)^5 * (2 + 3f * (b′0 - 4))
    end
end
function pressureeos2(p::BirchMurnaghan{3})
    v0, b0, b′0 = p.x0
    function (v)
        f = strain_from_volume(p)(v0, v)
        return 3f / 2 * b0 * sqrt(2f + 1)^5 * (2 + 3f * (b′0 - 4))
    end
end
function pressureeos3(p::BirchMurnaghan{3})
    v0, b0, b′0 = p.x0
    func = strain_from_volume(p)
    function (v)
        f = func(v0, v)
        return 3f / 2 * b0 * sqrt(2f + 1)^5 * (2 + 3f * (b′0 - 4))
    end
end
function pressureeos4(p::BirchMurnaghan{3})
    v0, b0, b′0 = p.x0
    function (v)
        f = BirchMurnaghan3rdStrain(v0, v)
        return 3f / 2 * b0 * sqrt(2f + 1)^5 * (2 + 3f * (b′0 - 4))
    end
end

I thought the 4th would be the fastest but they are too close so that the result is undetermined:

julia> eos = pressureeos(x);

julia> eos2 = pressureeos2(x);

julia> eos3 = pressureeos3(x);

julia> eos4 = pressureeos4(x);

julia> @btime $eos(4)
  71.107 ns (0 allocations: 0 bytes)
238.5808498748582

julia> @btime $eos2(4)
  71.128 ns (0 allocations: 0 bytes)
238.5808498748582

julia> @btime $eos3(4)
  71.005 ns (0 allocations: 0 bytes)
238.5808498748582

julia> @btime $eos4(4)
  71.027 ns (0 allocations: 0 bytes)
238.5808498748582

julia> @btime $eos(4)
  70.990 ns (0 allocations: 0 bytes)
238.5808498748582

julia> @btime $eos2(4)
  70.978 ns (0 allocations: 0 bytes)
238.5808498748582

julia> @btime $eos3(4)
  71.102 ns (0 allocations: 0 bytes)
238.5808498748582

julia> @btime $eos4(4)
  70.987 ns (0 allocations: 0 bytes)
238.5808498748582

Find volume from pressures or energies

Volume can be find from pressures or energies for some EOSs like Murnaghan, etc.

	nonlinfit(
	EnergyEquation(Murnaghan1st(19u"angstrom^3", 100u"GPa", 4)),
	volumes,
	energies,
	),

mineralscloud / equationsofstateofsolids.jl Goto Github PK

equationsofstateofsolids.jl's People

Contributors

Stargazers

Watchers

equationsofstateofsolids.jl's Issues

The type approach

The closure approach

Summary

Recommend Projects

Recommend Topics

Recommend Org