Merge pull request #92 from MineralsCloud/feature/strain

Implement `LinearFitting` module

Project.toml

@@ -7,7 +7,6 @@ version = "3.0.2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
@@ -17,7 +16,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 ConstructionBase = "1.0"
 IterTools = "1.2, 1.3"
 LsqFit = "0.8, 0.9, 0.10"
-Polynomials = "0.5.2, 0.6, 0.7, 0.8"
+Polynomials = "0.5.2, 0.6, 0.7, 0.8, 1.0"
 Roots = "0.8, 1.0"
 Unitful = "0.18, 1.0"
 julia = "1"

docs/src/api/Collections.md

@@ -10,7 +10,7 @@ The current `EquationOfState`s contain
 EquationOfState
 ├─ AntonSchmidt
 ├─ BreenanStacey
-├─ FiniteStrainEquationOfState
+├─ FiniteStrainEOS
 │  ├─ BirchMurnaghan2nd
 │  ├─ BirchMurnaghan3rd
 │  ├─ BirchMurnaghan4th
@@ -264,7 +264,7 @@ Energy
 Pressure
 BulkModulus
 EquationOfState
-FiniteStrainEquationOfState
+FiniteStrainEOS
 Murnaghan
 BirchMurnaghan2nd
 BirchMurnaghan3rd

src/Collections.jl

@@ -23,7 +23,8 @@ export Energy,
     Vinet,
     AntonSchmidt,
     BreenanStacey,
-    Shanker
+    Shanker,
+    PolynomialEOS
 
 # ============================================================================ #
 #                                     Types                                    #
@@ -142,11 +143,11 @@ An abstraction of equations of state, where `T` specifies the elements' common t
 abstract type EquationOfState{T} end
 
 """
-    FiniteStrainEquationOfState{T} <: EquationOfState{T}
+    FiniteStrainEOS{T} <: EquationOfState{T}
 
 An abstraction of finite strain equations of state, where `T` specifies the elements' common type.
 """
-abstract type FiniteStrainEquationOfState{T} <: EquationOfState{T} end
+abstract type FiniteStrainEOS{T} <: EquationOfState{T} end
 
 """
     Murnaghan(v0, b0, b′0, e0)
@@ -212,7 +213,7 @@ julia> BirchMurnaghan2nd(1u"nm^3", 2u"GPa", 3.0u"eV")
 BirchMurnaghan2nd{Quantity{Float64,D,U} where U where D}(1.0 nm³, 2.0 GPa, 3.0 eV)
 ```
 """
-struct BirchMurnaghan2nd{T} <: FiniteStrainEquationOfState{T}
+struct BirchMurnaghan2nd{T} <: FiniteStrainEOS{T}
     v0::T
     b0::T
     e0::T
@@ -250,7 +251,7 @@ julia> BirchMurnaghan3rd(1u"nm^3", 2u"GPa", 4.0, 3u"eV")
 BirchMurnaghan3rd{Quantity{Float64,D,U} where U where D}(1.0 nm³, 2.0 GPa, 4.0, 3.0 eV)
 ```
 """
-struct BirchMurnaghan3rd{T} <: FiniteStrainEquationOfState{T}
+struct BirchMurnaghan3rd{T} <: FiniteStrainEOS{T}
     v0::T
     b0::T
     b′0::T
@@ -290,7 +291,7 @@ julia> BirchMurnaghan4th(1u"nm^3", 2u"GPa", 3.0, 4u"1/GPa", 5u"eV")
 BirchMurnaghan4th{Quantity{Float64,D,U} where U where D}(1.0 nm³, 2.0 GPa, 3.0, 4.0 GPa⁻¹, 5.0 eV)
 ```
 """
-struct BirchMurnaghan4th{T} <: FiniteStrainEquationOfState{T}
+struct BirchMurnaghan4th{T} <: FiniteStrainEOS{T}
     v0::T
     b0::T
     b′0::T
@@ -330,7 +331,7 @@ julia> PoirierTarantola2nd(1u"nm^3", 2u"GPa", 3.0u"eV")
 PoirierTarantola2nd{Quantity{Float64,D,U} where U where D}(1.0 nm³, 2.0 GPa, 3.0 eV)
 ```
 """
-struct PoirierTarantola2nd{T} <: FiniteStrainEquationOfState{T}
+struct PoirierTarantola2nd{T} <: FiniteStrainEOS{T}
     v0::T
     b0::T
     e0::T
@@ -368,7 +369,7 @@ julia> PoirierTarantola3rd(1u"nm^3", 2u"GPa", 3, 4.0u"eV")
 PoirierTarantola3rd{Quantity{Float64,D,U} where U where D}(1.0 nm³, 2.0 GPa, 3.0, 4.0 eV)
 ```
 """
-struct PoirierTarantola3rd{T} <: FiniteStrainEquationOfState{T}
+struct PoirierTarantola3rd{T} <: FiniteStrainEOS{T}
     v0::T
     b0::T
     b′0::T
@@ -408,7 +409,7 @@ julia> PoirierTarantola4th(1u"nm^3", 2u"GPa", 3, 4u"1/GPa", 5.0u"eV")
 PoirierTarantola4th{Quantity{Float64,D,U} where U where D}(1.0 nm³, 2.0 GPa, 3.0, 4.0 GPa⁻¹, 5.0 eV)
 ```
 """
-struct PoirierTarantola4th{T} <: FiniteStrainEquationOfState{T}
+struct PoirierTarantola4th{T} <: FiniteStrainEOS{T}
     v0::T
     b0::T
     b′0::T
@@ -502,6 +503,16 @@ end
 Shanker(v0::Real, b0::Real, b′0::Real) = Shanker(v0, b0, b′0, 0)
 Shanker(v0::AbstractQuantity, b0::AbstractQuantity, b′0) =
     Shanker(v0, b0, b′0, 0 * upreferred(Unitful.J))
+
+struct PolynomialEOS{N,T} <: EquationOfState{T}
+    v0::T
+    b0::NTuple{N,T}
+    e0::T
+end
+function PolynomialEOS(v0, b0, e0)
+    T = Base.promote_typeof(v0, b0..., e0)
+    return PolynomialEOS{length(b0),T}(convert(T, v0), T.(b0), convert(T, e0))
+end
 # =================================== Types ================================== #
 
 # Energy evaluation

src/EquationsOfState.jl

@@ -2,7 +2,6 @@ module EquationsOfState
 
 include("Collections.jl")
 include("NonlinearFitting.jl")
-include("FiniteStrains.jl")
 include("LinearFitting.jl")
 include("Find.jl")
 

src/LinearFitting.jl

@@ -3,89 +3,67 @@ This module provides some linear fitting methods.
 """
 module LinearFitting
 
-using LinearAlgebra: dot
-using Polynomials: polyder, degree, coeffs, Poly
-using Polynomials.PolyCompat: polyfit
+using Polynomials: Polynomial, fit, derivative, roots, coeffs
 
-using ..FiniteStrains
+using ..Collections: PolynomialEOS
 
-export energy_strain_expansion,
-    energy_strain_derivative,
-    strain_volume_derivative,
-    energy_volume_expansion,
-    energy_volume_derivatives,
-    energy_volume_derivative_at_order
+export FiniteStrain, Eulerian, Lagrangian, Natural, Infinitesimal, linearfit
 
-energy_strain_expansion(f::Vector{<:Real}, e::Vector{<:Real}, n::Int)::Poly =
-    polyfit(f, e, n)
+abstract type FiniteStrain end
+struct Eulerian <: FiniteStrain end
+struct Lagrangian <: FiniteStrain end
+struct Natural <: FiniteStrain end
+struct Infinitesimal <: FiniteStrain end
 
-energy_strain_derivative(p::Poly, m::Int)::Poly = polyder(p, m)
-
-function strain_volume_derivative(s::EulerianStrain, v0::Real, v::Real, m::Int)
-    m == 1 && return -1 / (3v) * cbrt(v0 / v)^2
-    -(3m + 2) / (3v) * strain_volume_derivative(s, v0, v, m - 1)
-end
-function strain_volume_derivative(s::LagrangianStrain, v0::Real, v::Real, m::Int)
-    m == 1 && return -1 / (3v) * cbrt(v / v0)^2
-    -(3m - 2) / (3v) * strain_volume_derivative(s, v0, v, m - 1)
-end
-function strain_volume_derivative(s::NaturalStrain, v0::Real, v::Real, m::Int)
-    m == 1 && return 1 / (3v)
-    -m / v * strain_volume_derivative(s, v0, v, m - 1)
-end
-function strain_volume_derivative(s::InfinitesimalStrain, v0::Real, v::Real, m::Int)
-    m == 1 && return (1 - getstrain(s, v0, v))^4 / 3 / v0
-    -(3m + 1) / (3v) * strain_volume_derivative(s, v0, v, m - 1)
+struct StrainFromVolume{T<:FiniteStrain}
+    v0::Any
 end
+(x::StrainFromVolume{Eulerian})(v) = (cbrt(x.v0 / v)^2 - 1) / 2
+(x::StrainFromVolume{Lagrangian})(v) = (cbrt(v / x.v0)^2 - 1) / 2
+(x::StrainFromVolume{Natural})(v) = log(v / x.v0) / 3
+(x::StrainFromVolume{Infinitesimal})(v) = 1 - cbrt(x.v0 / v)
 
-function energy_volume_expansion(
-    s::FiniteStrain,
-    v0::Real,
-    v::Real,
-    p::Poly,
-    highest_order::Int = degree(p),
-)
-    # The zeroth order value plus values from the first to the ``highest_order`.
-    p(v) + dot(
-        energy_volume_derivatives(s, v0, v, p, highest_order),
-        getstrain(s, v0, v) .^ collect(1:highest_order),
-    )
+struct VolumeFromStrain{T<:FiniteStrain}
+    v0::Any
 end
+(x::VolumeFromStrain{Eulerian})(f) = x.v0 / (2f + 1)^(3 / 2)
+(x::VolumeFromStrain{Lagrangian})(f) = x.v0 * (2f + 1)^(3 / 2)
+(x::VolumeFromStrain{Natural})(f) = x.v0 * exp(3f)
+(x::VolumeFromStrain{Infinitesimal})(f) = x.v0 / (1 - f)^3
 
-function energy_volume_derivatives(
-    s::FiniteStrain,
-    v0::Real,
-    v::Real,
-    p::Poly,
-    highest_order::Int,
-)
-    0 ≤ highest_order ≤ degree(p) ? (x = 1:highest_order) :
-    throw(DomainError("The `highest_order` must be within 0 to $(degree(p))!"))
-    strain_derivatives = map(m -> strain_volume_derivative(s, v0, v, m), x)
-    energy_derivatives = map(f -> f(v), map(m -> energy_strain_derivative(p, m), x))
-    map(
-        m -> energy_volume_derivative_at_order(m)(strain_derivatives, energy_derivatives),
-        x,
-    )
-end
+_islocalminimum(f, x, δx) = f(x) < f(x - δx) && f(x) < f(x + δx)
 
-function energy_volume_derivative_at_order(m::Int)::Function
-    function (f::Vector{<:Real}, e::Vector{<:Real})
-        if m == 1
-            e[1] * f[1]
-        elseif m == 2
-            e[2] * f[1]^2 + e[1] * f[1]
-        elseif m == 3
-            e[3] * f[1]^3 + 3 * f[1] * f[2] * e[2] + e[1] * f[3]
-        elseif m == 4
-            e[4] * f[1]^4 +
-            6 * f[1]^2 * f[2] * e[3] +
-            (4 * f[1] * f[3] + 3 * f[3]^2) * e[2] +
-            e[1] * f[3]
-        else
-            error("Expansion is not defined at order = $(m)!")
+function _findglobalminimum(f, localminima, δx)
+    if length(localminima) == 0
+        error("no volume minimizes the energy!")
+    else
+        f0, i = findmin(f.(localminima))
+        x0 = localminima[i]
+        return x0, f0
+    end
+end # function _findglobalminimum
+
+function linearfit(volumes, energies, deg = 3)
+    poly = fit(volumes, energies, deg)
+    poly1d = derivative(poly, 1)
+    δx = minimum(diff(volumes)) / 10
+    localminima = eltype(volumes)[]
+    for x in real(filter(isreal, roots(poly1d)))  # Complex volume is meaningless
+        if _islocalminimum(poly, x, δx)
+            push!(localminima, x)
         end
     end
-end
+    v0, e0 = _findglobalminimum(poly, localminima, δx)
+    bs = Tuple(derivative(poly, n)(v0) / factorial(n) for n in 1:deg)
+    return PolynomialEOS(v0, bs, e0)
+end # function linearfit
+
+Base.inv(x::StrainFromVolume{T}) where {T} = VolumeFromStrain{T}(x.v0)
+Base.inv(x::VolumeFromStrain{T}) where {T} = StrainFromVolume{T}(x.v0)
+
+Base.:∘(a::StrainFromVolume{T}, b::VolumeFromStrain{T}) where {T} =
+    a.v0 == b.v0 ? identity : error("operation undefined!")
+Base.:∘(a::VolumeFromStrain{T}, b::StrainFromVolume{T}) where {T} =
+    a.v0 == b.v0 ? identity : error("operation undefined!")
 
 end