diff --git a/.github/PULL_REQUEST_TEMPLATE/conversational.md b/.github/PULL_REQUEST_TEMPLATE.md
old mode 100755
new mode 100644
similarity index 100%
rename from .github/PULL_REQUEST_TEMPLATE/conversational.md
rename to .github/PULL_REQUEST_TEMPLATE.md
diff --git a/.github/PULL_REQUEST_TEMPLATE/bugs-only.md b/.github/PULL_REQUEST_TEMPLATE/bugs-only.md
deleted file mode 100755
index 61262fa..0000000
--- a/.github/PULL_REQUEST_TEMPLATE/bugs-only.md
+++ /dev/null
@@ -1,23 +0,0 @@
-Bug Fixes **ONLY**. NO NEW FEATURE ACCEPTED!
-
-<!--- Provide a general summary of your changes in the Title above -->
-
-## Description
-<!--- Describe your changes in detail -->
-
-## Related Issue
-<!--- This project only accepts pull requests related to open issues -->
-<!--- If suggesting a new feature or change, please discuss it in an issue first -->
-<!--- If fixing a bug, there should be an issue describing it with steps to reproduce -->
-<!--- Please link to the issue here: -->
-
-## Motivation and Context
-<!--- Why is this change required? What problem does it solve? -->
-<!--- If it fixes an open issue, please link to the issue here. -->
-
-## How Has This Been Tested?
-<!--- Please describe in detail how you tested your changes. -->
-<!--- Include details of your testing environment, and the tests you ran to -->
-<!--- see how your change affects other areas of the code, etc. -->
-
-## Screenshots (if appropriate):
diff --git a/Manifest.toml b/Manifest.toml
new file mode 100644
index 0000000..d4b4bec
--- /dev/null
+++ b/Manifest.toml
@@ -0,0 +1,307 @@
+# This file is machine-generated - editing it directly is not advised
+
+[[Arpack]]
+deps = ["BinaryProvider", "Libdl", "LinearAlgebra"]
+git-tree-sha1 = "07a2c077bdd4b6d23a40342a8a108e2ee5e58ab6"
+uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+version = "0.3.1"
+
+[[ArrayInterface]]
+deps = ["LinearAlgebra", "Requires", "SparseArrays"]
+git-tree-sha1 = "981354dab938901c2b607a213e62d9defa50b698"
+uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+version = "1.2.1"
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[BinDeps]]
+deps = ["Compat", "Libdl", "SHA", "URIParser"]
+git-tree-sha1 = "12093ca6cdd0ee547c39b1870e0c9c3f154d9ca9"
+uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee"
+version = "0.8.10"
+
+[[BinaryProvider]]
+deps = ["Libdl", "Logging", "SHA"]
+git-tree-sha1 = "c7361ce8a2129f20b0e05a89f7070820cfed6648"
+uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
+version = "0.5.6"
+
+[[Calculus]]
+deps = ["Compat"]
+git-tree-sha1 = "bd8bbd105ba583a42385bd6dc4a20dad8ab3dc11"
+uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
+version = "0.5.0"
+
+[[CommonSubexpressions]]
+deps = ["Test"]
+git-tree-sha1 = "efdaf19ab11c7889334ca247ff4c9f7c322817b0"
+uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
+version = "0.2.0"
+
+[[Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "84aa74986c5b9b898b0d1acaf3258741ee64754f"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "2.1.0"
+
+[[ConstructionBase]]
+git-tree-sha1 = "e3efe0a0f49dcd294c8c73e897b4fdf891f0fbd3"
+uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+version = "0.1.0"
+
+[[DataAPI]]
+git-tree-sha1 = "674b67f344687a88310213ddfa8a2b3c76cc4252"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.1.0"
+
+[[DataStructures]]
+deps = ["InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "f94423c68f2e47db0d6f626a26d4872266e0ec3d"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.17.2"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[DiffEqDiffTools]]
+deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
+git-tree-sha1 = "21b855cb29ec4594f9651e0e9bdc0cdcfdcd52c1"
+uuid = "01453d9d-ee7c-5054-8395-0335cb756afa"
+version = "1.3.0"
+
+[[DiffResults]]
+deps = ["Compat", "StaticArrays"]
+git-tree-sha1 = "34a4a1e8be7bc99bc9c611b895b5baf37a80584c"
+uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
+version = "0.0.4"
+
+[[DiffRules]]
+deps = ["Random", "Test"]
+git-tree-sha1 = "dc0869fb2f5b23466b32ea799bd82c76480167f7"
+uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
+version = "0.0.10"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[Distributions]]
+deps = ["LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
+git-tree-sha1 = "b419fcf95ef9c8cf4d6610cd323890ad66d64240"
+uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
+version = "0.21.3"
+
+[[ForwardDiff]]
+deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "InteractiveUtils", "LinearAlgebra", "NaNMath", "Random", "SparseArrays", "SpecialFunctions", "StaticArrays", "Test"]
+git-tree-sha1 = "4c4d727f1b7e0092134fabfab6396b8945c1ea5b"
+uuid = "f6369f11-7733-5829-9624-2563aa707210"
+version = "0.10.3"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[LibGit2]]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[LsqFit]]
+deps = ["Distributions", "LinearAlgebra", "NLSolversBase", "OptimBase", "Random", "StatsBase", "Test"]
+git-tree-sha1 = "186c2afbdb3cd52191078cfc6176f7084ed9dfb7"
+uuid = "2fda8390-95c7-5789-9bda-21331edee243"
+version = "0.8.1"
+
+[[MLStyle]]
+git-tree-sha1 = "67f9a88611bc79f992aa705d9bbc833a2547dec7"
+uuid = "d8e11817-5142-5d16-987a-aa16d5891078"
+version = "0.3.1"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[Missings]]
+deps = ["DataAPI"]
+git-tree-sha1 = "de0a5ce9e5289f27df672ffabef4d1e5861247d5"
+uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
+version = "0.4.3"
+
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[NLSolversBase]]
+deps = ["Calculus", "DiffEqDiffTools", "DiffResults", "Distributed", "ForwardDiff"]
+git-tree-sha1 = "f1b8ed89fa332f410cfc7c937682eb4d0b361521"
+uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
+version = "7.5.0"
+
+[[NaNMath]]
+deps = ["Compat"]
+git-tree-sha1 = "ce3b85e484a5d4c71dd5316215069311135fa9f2"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "0.3.2"
+
+[[OptimBase]]
+deps = ["Compat", "NLSolversBase", "Printf", "Reexport", "Test"]
+git-tree-sha1 = "92667ab46a66ad502ec3044f65c41ea68b2e0e9c"
+uuid = "87e2bd06-a317-5318-96d9-3ecbac512eee"
+version = "2.0.0"
+
+[[OrderedCollections]]
+deps = ["Random", "Serialization", "Test"]
+git-tree-sha1 = "c4c13474d23c60d20a67b217f1d7f22a40edf8f1"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.1.0"
+
+[[PDMats]]
+deps = ["Arpack", "LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"]
+git-tree-sha1 = "035f8d60ba2a22cb1d2580b1e0e5ce0cb05e4563"
+uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
+version = "0.9.10"
+
+[[Pkg]]
+deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[Polynomials]]
+deps = ["LinearAlgebra", "SparseArrays", "Test"]
+git-tree-sha1 = "62142bd65d3f8aeb2226ec64dd8493349147df94"
+uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+version = "0.5.2"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[QuadGK]]
+deps = ["DataStructures", "LinearAlgebra", "Test"]
+git-tree-sha1 = "3ce467a8e76c6030d4c3786e7d3a73442017cdc0"
+uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+version = "2.0.3"
+
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[Reexport]]
+deps = ["Pkg"]
+git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "0.2.0"
+
+[[Requires]]
+deps = ["Test"]
+git-tree-sha1 = "f6fbf4ba64d295e146e49e021207993b6b48c7d1"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "0.5.2"
+
+[[Rmath]]
+deps = ["BinaryProvider", "Libdl", "Random", "Statistics", "Test"]
+git-tree-sha1 = "9a6c758cdf73036c3239b0afbea790def1dabff9"
+uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
+version = "0.5.0"
+
+[[Roots]]
+deps = ["Printf"]
+git-tree-sha1 = "9cc4b586c71f9aea25312b94be8c195f119b0ec3"
+uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+version = "0.8.3"
+
+[[SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[SortingAlgorithms]]
+deps = ["DataStructures", "Random", "Test"]
+git-tree-sha1 = "03f5898c9959f8115e30bc7226ada7d0df554ddd"
+uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
+version = "0.3.1"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[SpecialFunctions]]
+deps = ["BinDeps", "BinaryProvider", "Libdl", "Test"]
+git-tree-sha1 = "0b45dc2e45ed77f445617b99ff2adf0f5b0f23ea"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "0.7.2"
+
+[[StaticArrays]]
+deps = ["LinearAlgebra", "Random", "Statistics"]
+git-tree-sha1 = "db23bbf50064c582b6f2b9b043c8e7e98ea8c0c6"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "0.11.0"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[StatsBase]]
+deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"]
+git-tree-sha1 = "c53e809e63fe5cf5de13632090bc3520649c9950"
+uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+version = "0.32.0"
+
+[[StatsFuns]]
+deps = ["Rmath", "SpecialFunctions", "Test"]
+git-tree-sha1 = "b3a4e86aa13c732b8a8c0ba0c3d3264f55e6bb3e"
+uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+version = "0.8.0"
+
+[[SuiteSparse]]
+deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
+uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+
+[[Test]]
+deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[URIParser]]
+deps = ["Test", "Unicode"]
+git-tree-sha1 = "6ddf8244220dfda2f17539fa8c9de20d6c575b69"
+uuid = "30578b45-9adc-5946-b283-645ec420af67"
+version = "0.4.0"
+
+[[UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[Unitful]]
+deps = ["LinearAlgebra", "Random"]
+git-tree-sha1 = "0e9ad8a9621a9e1ec48f7cb48a5d6ded9e2c68f3"
+repo-rev = "master"
+repo-url = "https://github.com/PainterQubits/Unitful.jl.git"
+uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
+version = "0.17.0"
diff --git a/Project.toml b/Project.toml
index 6b9b24e..6af4428 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,19 +4,21 @@ authors = ["Qi Zhang <singularitti@outlook.com>"]
 version = "2.0.0"
 
 [deps]
+ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
 [targets]
-test = ["Test"]
+test = ["Test", "UnitfulAtomic"]
diff --git a/README.md b/README.md
index 1051c6f..1a9b709 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
   <img src="./docs/src/assets/logo.png" height="200"><br>
 </div>
 
-------
+---
 
 # EquationsOfState.jl
 
@@ -17,30 +17,33 @@
 ![GitHub release (latest by date including pre-releases)](https://img.shields.io/github/v/release/MineralsCloud/EquationsOfState.jl?include_prereleases)
 [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle)
 
-This package implements some _equations of state_ (EOS) of solids which are useful in research. It currently includes:
+This package implements some _equations of state_ (EOS) of solids which are
+useful in research. It currently includes:
 
 1. `Murnaghan` EOS
 2. Birch–Murnaghan EOS family:
-    1. `BirchMurnaghan2nd`
-    2. `BirchMurnaghan3rd`
-    3. `BirchMurnaghan4th`
+   1. `BirchMurnaghan2nd`
+   2. `BirchMurnaghan3rd`
+   3. `BirchMurnaghan4th`
 3. `Vinet` EOS
 4. Poirier–Tarantola EOS family:
-    1. `PoirierTarantola2nd`
-    2. `PoirierTarantola3rd`
-    3. `PoirierTarantola4th`
+   1. `PoirierTarantola2nd`
+   2. `PoirierTarantola3rd`
+   3. `PoirierTarantola4th`
 5. `AntonSchmidt` EOS (experimental)
 6. `BreenanStacey` EOS (experimental)
 
 The formula are referenced from Ref. 1.
 
-This package also includes linear and nonlinear fitting methods, also referenced from Ref. 1.
+This package also includes linear and nonlinear fitting methods, also referenced
+from Ref. 1.
 
 ## Compatibility
 
 - Julia version: `v1.0.0` and above
-- Dependencies: see [`Project.toml`]([Project.toml](https://github.com/MineralsCloud/EquationsOfState.jl/blob/master/Project.toml))
-- OS versions: `macOS`, `Linux`, and `Windows`
+- Dependencies: see
+  [`Project.toml`](https://github.com/MineralsCloud/EquationsOfState.jl/blob/master/Project.toml)
+- OS versions: macOS, Linux, and Windows
 
 ## Installation
 
@@ -52,7 +55,7 @@ This package also includes linear and nonlinear fitting methods, also referenced
 
    ```julia
    julia> using Pkg
-   
+
    julia> Pkg.add("https://github.com/MineralsCloud/EquationsOfState.jl")
    ```
 
@@ -60,7 +63,7 @@ This package also includes linear and nonlinear fitting methods, also referenced
 
    ```julia
    julia> using Pkg
-   
+
    julia> Pkg.add("EquationsOfState")
    ```
 
@@ -68,14 +71,15 @@ This package also includes linear and nonlinear fitting methods, also referenced
 
 ## TODOs
 
-- [ ] Implement nonlinear fitting using [CMPFit.jl](https://github.com/gcalderone/CMPFit.jl).
+- [ ] Implement nonlinear fitting using
+      [CMPFit.jl](https://github.com/gcalderone/CMPFit.jl).
 - [ ] Finish [docs](https://mineralscloud.github.io/EquationsOfState.jl/)
 
 ## Related packages
 
 1. [CommandLineEquationsOfState.jl](https://github.com/MineralsCloud/CommandLineEquationsOfState.jl)
-2. [ExtendedEquationsOfState.jl](https://github.com/MineralsCloud/ExtendedEquationsOfState.jl)
 
 ## References
 
-1. A. Otero-De-La-Roza, V. Luaña, *Comput. Phys. Commun.* **182**, 1708–1720 (2011).
+1. A. Otero-De-La-Roza, V. Luaña, _Comput. Phys. Commun._ **182**, 1708–1720
+   (2011).
diff --git a/src/Collections.jl b/src/Collections.jl
index 062fc38..6af09b2 100644
--- a/src/Collections.jl
+++ b/src/Collections.jl
@@ -12,13 +12,14 @@ julia>
 module Collections
 
 using InteractiveUtils
+using Unitful: AbstractQuantity, @u_str, Dimension, Dimensions, upreferred
+import Unitful
 
-using EquationsOfState
+using EquationsOfState: EnergyForm, PressureForm, BulkModulusForm
 
 export apply,
        EquationOfState,
        FiniteStrainEquationOfState,
-    #    PolynomialEquationOfState,
        Murnaghan,
        BirchMurnaghan2nd,
        BirchMurnaghan3rd,
@@ -47,18 +48,6 @@ An abstraction of finite strain equations of state.
 """
 abstract type FiniteStrainEquationOfState{T} <: EquationOfState{T} end
 
-# struct PolynomialEquationOfState{T<:Real,N} <: EquationOfState{T,N}
-#     data::NTuple{N,T}
-#     function PolynomialEquationOfState{T,N}(args::NTuple{N,T}) where {T,N}
-#         @assert N ≤ 10
-#         new(args)
-#     end
-# end
-# function PolynomialEquationOfState(args...)
-#     T = Base.promote_typeof(args...)
-#     PolynomialEquationOfState{T,length(args)}(args)
-# end
-
 """
     Murnaghan(v0, b0, bp0, e0=0)
 
@@ -70,17 +59,19 @@ Create a Murnaghan equation of state. The elements' type will be handled automat
 - `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
 - `e0=0`: the energy of solid at zero pressure. By default is `0`.
 """
-struct Murnaghan{T<:Real} <: EquationOfState{T}
+struct Murnaghan{T} <: EquationOfState{T}
     v0::T
     b0::T
     bp0::T
     e0::T
 end
-function Murnaghan(v0::Real, b0::Real, bp0::Real, e0::Real)
+function Murnaghan(v0, b0, bp0, e0)
     T = Base.promote_typeof(v0, b0, bp0, e0)
-    Murnaghan{T}(v0, b0, bp0, e0)
+    return Murnaghan{T}(convert.(T, [v0, b0, bp0, e0])...)
 end
-Murnaghan(v0, b0, bp0) = Murnaghan(v0, b0, bp0, 0)
+Murnaghan(v0::Real, b0::Real, bp0::Real) = Murnaghan(v0, b0, bp0, 0)
+Murnaghan(v0::AbstractQuantity, b0::AbstractQuantity, bp0) =
+    Murnaghan(v0, b0, bp0, 0 * upreferred(Unitful.J))
 
 """
     BirchMurnaghan2nd(v0, b0, e0=0)
@@ -92,16 +83,18 @@ Create a Birch–Murnaghan 2nd order equation of state. The elements' type will
 - `b0`: the bulk modulus of solid at zero pressure.
 - `e0=0`: the energy of solid at zero pressure. By default is `0`.
 """
-struct BirchMurnaghan2nd{T<:Real} <: FiniteStrainEquationOfState{T}
+struct BirchMurnaghan2nd{T} <: FiniteStrainEquationOfState{T}
     v0::T
     b0::T
     e0::T
 end
-function BirchMurnaghan2nd(v0::Real, b0::Real, e0::Real)
+function BirchMurnaghan2nd(v0, b0, e0)
     T = Base.promote_typeof(v0, b0, e0)
-    BirchMurnaghan2nd{T}(v0, b0, e0)
+    return BirchMurnaghan2nd{T}(convert.(T, [v0, b0, e0])...)
 end
-BirchMurnaghan2nd(v0, b0) = BirchMurnaghan2nd(v0, b0, 0)
+BirchMurnaghan2nd(v0::Real, b0::Real) = BirchMurnaghan2nd(v0, b0, 0)
+BirchMurnaghan2nd(v0::AbstractQuantity, b0::AbstractQuantity) =
+    BirchMurnaghan2nd(v0, b0, 0 * upreferred(Unitful.J))
 
 """
     BirchMurnaghan3rd(v0, b0, bp0, e0=0)
@@ -114,17 +107,19 @@ Create a Birch–Murnaghan 3rd order equation of state. The elements' type will
 - `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
 - `e0=0`: the energy of solid at zero pressure. By default is `0`.
 """
-struct BirchMurnaghan3rd{T<:Real} <: FiniteStrainEquationOfState{T}
+struct BirchMurnaghan3rd{T} <: FiniteStrainEquationOfState{T}
     v0::T
     b0::T
     bp0::T
     e0::T
 end
-function BirchMurnaghan3rd(v0::Real, b0::Real, bp0::Real, e0::Real)
+function BirchMurnaghan3rd(v0, b0, bp0, e0)
     T = Base.promote_typeof(v0, b0, bp0, e0)
-    BirchMurnaghan3rd{T}(v0, b0, bp0, e0)
+    return BirchMurnaghan3rd{T}(convert.(T, [v0, b0, bp0, e0])...)
 end
-BirchMurnaghan3rd(v0, b0, bp0) = BirchMurnaghan3rd(v0, b0, bp0, 0)
+BirchMurnaghan3rd(v0::Real, b0::Real, bp0::Real) = BirchMurnaghan3rd(v0, b0, bp0, 0)
+BirchMurnaghan3rd(v0::AbstractQuantity, b0::AbstractQuantity, bp0) =
+    BirchMurnaghan3rd(v0, b0, bp0, 0 * upreferred(Unitful.J))
 
 """
     BirchMurnaghan4th(v0, b0, bp0, bpp0, e0=0)
@@ -138,18 +133,25 @@ Create a Birch–Murnaghan 4th order equation of state. The elements' type will
 - `bpp0`: the second-order pressure-derivative bulk modulus of solid at zero pressure.
 - `e0=0`: the energy of solid at zero pressure. By default is `0`.
 """
-struct BirchMurnaghan4th{T<:Real} <: FiniteStrainEquationOfState{T}
+struct BirchMurnaghan4th{T} <: FiniteStrainEquationOfState{T}
     v0::T
     b0::T
     bp0::T
     bpp0::T
     e0::T
 end
-function BirchMurnaghan4th(v0::Real, b0::Real, bp0::Real, bpp0::Real, e0::Real)
+function BirchMurnaghan4th(v0, b0, bp0, bpp0, e0)
     T = Base.promote_typeof(v0, b0, bp0, bpp0, e0)
-    BirchMurnaghan4th{T}(v0, b0, bp0, bpp0, e0)
+    return BirchMurnaghan4th{T}(convert.(T, [v0, b0, bp0, bpp0, e0])...)
 end
-BirchMurnaghan4th(v0, b0, bp0, bpp0) = BirchMurnaghan4th(v0, b0, bp0, bpp0, 0)
+BirchMurnaghan4th(v0::Real, b0::Real, bp0::Real, bpp0::Real) =
+    BirchMurnaghan4th(v0, b0, bp0, bpp0, 0)
+BirchMurnaghan4th(
+    v0::AbstractQuantity,
+    b0::AbstractQuantity,
+    bp0,
+    bpp0::AbstractQuantity,
+) = BirchMurnaghan4th(v0, b0, bp0, bpp0, 0 * upreferred(Unitful.J))
 
 """
     PoirierTarantola2nd(v0, b0, e0=0)
@@ -161,16 +163,18 @@ Create a Poirier–Tarantola order equation of state. The elements' type will be
 - `b0`: the bulk modulus of solid at zero pressure.
 - `e0=0`: the energy of solid at zero pressure. By default is `0`.
 """
-struct PoirierTarantola2nd{T<:Real} <: FiniteStrainEquationOfState{T}
+struct PoirierTarantola2nd{T} <: FiniteStrainEquationOfState{T}
     v0::T
     b0::T
     e0::T
 end
-function PoirierTarantola2nd(v0::Real, b0::Real, e0::Real)
+function PoirierTarantola2nd(v0, b0, e0)
     T = Base.promote_typeof(v0, b0, e0)
-    PoirierTarantola2nd{T}(v0, b0, e0)
+    return PoirierTarantola2nd{T}(convert.(T, [v0, b0, e0])...)
 end
-PoirierTarantola2nd(v0, b0) = PoirierTarantola2nd(v0, b0, 0)
+PoirierTarantola2nd(v0::Real, b0::Real) = PoirierTarantola2nd(v0, b0, 0)
+PoirierTarantola2nd(v0::AbstractQuantity, b0::AbstractQuantity) =
+    PoirierTarantola2nd(v0, b0, 0 * upreferred(Unitful.J))
 
 """
     PoirierTarantola3rd(v0, b0, bp0, e0=0)
@@ -183,17 +187,19 @@ Create a Poirier–Tarantola 3rd order equation of state. The elements' type wil
 - `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
 - `e0=0`: the energy of solid at zero pressure. By default is `0`.
 """
-struct PoirierTarantola3rd{T<:Real} <: FiniteStrainEquationOfState{T}
+struct PoirierTarantola3rd{T} <: FiniteStrainEquationOfState{T}
     v0::T
     b0::T
     bp0::T
     e0::T
 end
-function PoirierTarantola3rd(v0::Real, b0::Real, bp0::Real, e0::Real)
+function PoirierTarantola3rd(v0, b0, bp0, e0)
     T = Base.promote_typeof(v0, b0, bp0, e0)
-    PoirierTarantola3rd{T}(v0, b0, bp0, e0)
+    return PoirierTarantola3rd{T}(convert.(T, [v0, b0, bp0, e0])...)
 end
-PoirierTarantola3rd(v0, b0, bp0) = PoirierTarantola3rd(v0, b0, bp0, 0)
+PoirierTarantola3rd(v0::Real, b0::Real, bp0::Real) = PoirierTarantola3rd(v0, b0, bp0, 0)
+PoirierTarantola3rd(v0::AbstractQuantity, b0::AbstractQuantity, bp0) =
+    PoirierTarantola3rd(v0, b0, bp0, 0 * upreferred(Unitful.J))
 
 """
     PoirierTarantola4th(v0, b0, bp0, bpp0, e0=0)
@@ -207,18 +213,25 @@ Create a Birch–Murnaghan 4th order equation of state. The elements' type will
 - `bpp0`: the second-order pressure-derivative bulk modulus of solid at zero pressure.
 - `e0=0`: the energy of solid at zero pressure. By default is `0`.
 """
-struct PoirierTarantola4th{T<:Real} <: FiniteStrainEquationOfState{T}
+struct PoirierTarantola4th{T} <: FiniteStrainEquationOfState{T}
     v0::T
     b0::T
     bp0::T
     bpp0::T
     e0::T
 end
-function PoirierTarantola4th(v0::Real, b0::Real, bp0::Real, bpp0::Real, e0::Real)
+function PoirierTarantola4th(v0, b0, bp0, bpp0, e0)
     T = Base.promote_typeof(v0, b0, bp0, bpp0, e0)
-    PoirierTarantola4th{T}(v0, b0, bp0, bpp0, e0)
+    return PoirierTarantola4th{T}(convert.(T, [v0, b0, bp0, bpp0, e0])...)
 end
-PoirierTarantola4th(v0, b0, bp0, bpp0) = PoirierTarantola4th(v0, b0, bp0, bpp0, 0)
+PoirierTarantola4th(v0::Real, b0::Real, bp0::Real, bpp0::Real) =
+    PoirierTarantola4th(v0, b0, bp0, bpp0, 0)
+PoirierTarantola4th(
+    v0::AbstractQuantity,
+    b0::AbstractQuantity,
+    bp0,
+    bpp0::AbstractQuantity,
+) = PoirierTarantola4th(v0, b0, bp0, bpp0, 0 * upreferred(Unitful.J))
 
 """
     Vinet(v0, b0, bp0, e0=0)
@@ -231,41 +244,43 @@ Create a Vinet equation of state. The elements' type will be handled automatical
 - `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
 - `e0=0`: the energy of solid at zero pressure. By default is `0`.
 """
-struct Vinet{T<:Real} <: EquationOfState{T}
+struct Vinet{T} <: EquationOfState{T}
     v0::T
     b0::T
     bp0::T
     e0::T
 end
-function Vinet(v0::Real, b0::Real, bp0::Real, e0::Real)
+function Vinet(v0, b0, bp0, e0)
     T = Base.promote_typeof(v0, b0, bp0, e0)
-    Vinet{T}(v0, b0, bp0, e0)
+    return Vinet{T}(convert.(T, [v0, b0, bp0, e0])...)
 end
-Vinet(v0, b0, bp0) = Vinet(v0, b0, bp0, 0)
+Vinet(v0::Real, b0::Real, bp0::Real) = Vinet(v0, b0, bp0, 0)
+Vinet(v0::AbstractQuantity, b0::AbstractQuantity, bp0) =
+    Vinet(v0, b0, bp0, 0 * upreferred(Unitful.J))
 
-struct AntonSchmidt{T<:Real} <: EquationOfState{T}
+struct AntonSchmidt{T} <: EquationOfState{T}
     v0::T
     β::T
     n::T
     e∞::T
 end
-function AntonSchmidt(v0::Real, β::Real, n::Real, e∞::Real)
+function AntonSchmidt(v0, β, n, e∞)
     T = Base.promote_typeof(v0, β, n, e∞)
-    AntonSchmidt{T}(v0, β, n, e∞)
+    return AntonSchmidt{T}(convert.(T, [v0, β, n, e∞])...)
 end
-AntonSchmidt(v0, β, n) = AntonSchmidt(v0, β, n, 0)
+AntonSchmidt(v0::Real, β::Real, n::Real) = AntonSchmidt(v0, β, n, 0)
 
-struct BreenanStacey{T<:Real} <: EquationOfState{T}
+struct BreenanStacey{T} <: EquationOfState{T}
     v0::T
     b0::T
     γ0::T
     e0::T
 end
-function BreenanStacey(v0::Real, b0::Real, γ0::Real, e0::Real)
+function BreenanStacey(v0, b0, γ0, e0)
     T = Base.promote_typeof(v0, b0, γ0, e0)
-    BreenanStacey{T}(v0, b0, γ0, e0)
+    return BreenanStacey{T}(convert.(T, [v0, b0, γ0, e0])...)
 end
-BreenanStacey(v0, b0, γ0) = BreenanStacey(v0, b0, γ0, 0)
+BreenanStacey(v0::Real, b0::Real, γ0::Real) = BreenanStacey(v0, b0, γ0, 0)
 # =================================== Types ================================== #
 
 
@@ -299,81 +314,82 @@ julia> map(f, 1:1:10)
 """
 apply(form::EnergyForm, eos::EquationOfState) = v -> apply(form, eos, v)
 """
-    apply(EnergyForm(), eos::Murnaghan, v::Real)
+    apply(EnergyForm(), eos::Murnaghan, v)
 
 Return the energy of a `Murnaghan` equation of state on volume `v`.
 """
-function apply(::EnergyForm, eos::Murnaghan, v::Real)
-    v0, b0, bp0, e0 = collect(eos)
+function apply(::EnergyForm, eos::Murnaghan, v)
+    v0, b0, bp0, e0 = fieldvalues(eos)
 
     x = bp0 - 1
     y = (v0 / v)^bp0
     return e0 + b0 / bp0 * v * (y / x + 1) - v0 * b0 / x
 end
 """
-    apply(EnergyForm(), eos::BirchMurnaghan2nd, v::Real)
+    apply(EnergyForm(), eos::BirchMurnaghan2nd, v)
 
 Return the energy of a `BirchMurnaghan2nd` equation of state on volume `v`.
 """
-function apply(::EnergyForm, eos::BirchMurnaghan2nd, v::Real)
-    v0, b0, e0 = collect(eos)
+function apply(::EnergyForm, eos::BirchMurnaghan2nd, v)
+    v0, b0, e0 = fieldvalues(eos)
 
     f = (cbrt(v0 / v)^2 - 1) / 2
     return e0 + 9 / 2 * b0 * v0 * f^2
 end
 """
-    apply(EnergyForm(), eos::BirchMurnaghan3rd, v::Real)
+    apply(EnergyForm(), eos::BirchMurnaghan3rd, v)
 
 Return the energy of a `BirchMurnaghan3rd` equation of state on volume `v`.
 """
-function apply(::EnergyForm, eos::BirchMurnaghan3rd, v::Real)
-    v0, b0, bp0, e0 = collect(eos)
+function apply(::EnergyForm, eos::BirchMurnaghan3rd, v)
+    v0, b0, bp0, e0 = fieldvalues(eos)
 
     eta = cbrt(v0 / v)
     xi = eta^2 - 1
     return e0 + 9 / 16 * b0 * v0 * xi^2 * (6 + bp0 * xi - 4 * eta^2)
 end
 """
-    apply(EnergyForm(), eos::BirchMurnaghan4th, v::Real)
+    apply(EnergyForm(), eos::BirchMurnaghan4th, v)
 
 Return the energy of a `BirchMurnaghan4th` equation of state on volume `v`.
 """
-function apply(::EnergyForm, eos::BirchMurnaghan4th, v::Real)
-    v0, b0, bp0, bpp0, e0 = collect(eos)
+function apply(::EnergyForm, eos::BirchMurnaghan4th, v)
+    v0, b0, bp0, bpp0, e0 = fieldvalues(eos)
 
     f = (cbrt(v0 / v)^2 - 1) / 2
     h = b0 * bpp0 + bp0^2
-    return e0 + 3 / 8 * v0 * b0 * f^2 * ((9h - 63bp0 + 143) * f^2 + 12 * (bp0 - 4) * f + 12)
+    return e0 +
+           3 / 8 * v0 * b0 * f^2 * ((9h - 63bp0 + 143) * f^2 + 12 * (bp0 - 4) * f + 12)
 end
 """
-    apply(EnergyForm(), eos::PoirierTarantola2nd, v::Real)
+    apply(EnergyForm(), eos::PoirierTarantola2nd, v)
 
 Return the energy of a `PoirierTarantola2nd` equation of state on volume `v`.
 """
-function apply(::EnergyForm, eos::PoirierTarantola2nd, v::Real)
-    v0, b0, e0 = collect(eos)
+function apply(::EnergyForm, eos::PoirierTarantola2nd, v)
+    v0, b0, e0 = fieldvalues(eos)
 
     return e0 + b0 / 2 * v0 * log(v / v0)^(2 / 3)
 end
 """
-    apply(EnergyForm(), eos::PoirierTarantola3rd, v::Real)
+    apply(EnergyForm(), eos::PoirierTarantola3rd, v)
 
 Return the energy of a `PoirierTarantola3rd` equation of state on volume `v`.
 """
-function apply(::EnergyForm, eos::PoirierTarantola3rd, v::Real)
-    v0, b0, bp0, e0 = collect(eos)
+function apply(::EnergyForm, eos::PoirierTarantola3rd, v)
+    v0, b0, bp0, e0 = fieldvalues(eos)
 
     x = cbrt(v / v0)
     xi = -3 * log(x)
     return e0 + b0 / 6 * v0 * xi^2 * ((bp0 - 2) * xi + 3)
 end
 """
-    apply(EnergyForm(), eos::PoirierTarantola4th, v::Real)
+    apply(EnergyForm(), eos::PoirierTarantola4th, v)
 
 Return the energy of a `PoirierTarantola4th` equation of state on volume `v`.
 """
-function apply(::EnergyForm, eos::PoirierTarantola4th, v::Real)
-    v0, b0, bp0, bpp0, e0 = collect(eos)
+function apply(::EnergyForm, eos::PoirierTarantola4th, v)
+    v0, b0, bp0, bpp0, e0 = fieldvalues(eos)
 
     x = cbrt(v / v0)
     xi = log(x)
@@ -381,24 +397,24 @@ function apply(::EnergyForm, eos::PoirierTarantola4th, v::Real)
     return e0 + b0 / 24v0 * xi^2 * ((h + 3bp0 + 3) * xi^2 + 4 * (bp0 + 2) * xi + 12)
 end
 """
-    apply(EnergyForm(), eos::Vinet, v::Real)
+    apply(EnergyForm(), eos::Vinet, v)
 
 Return the energy of a `Vinet` equation of state on volume `v`.
 """
-function apply(::EnergyForm, eos::Vinet, v::Real)
-    v0, b0, bp0, e0 = collect(eos)
+function apply(::EnergyForm, eos::Vinet, v)
+    v0, b0, bp0, e0 = fieldvalues(eos)
 
     x = cbrt(v / v0)
     xi = 3 / 2 * (bp0 - 1)
     return e0 + 9b0 * v0 / xi^2 * (1 + (xi * (1 - x) - 1) * exp(xi * (1 - x)))
 end
 """
-    apply(EnergyForm(), eos::AntonSchmidt, v::Real)
+    apply(EnergyForm(), eos::AntonSchmidt, v)
 
 Return the energy of a `AntonSchmidt` equation of state on volume `v`.
 """
-function apply(::EnergyForm, eos::AntonSchmidt, v::Real)
-    v0, β, n, e∞ = collect(eos)
+function apply(::EnergyForm, eos::AntonSchmidt, v)
+    v0, β, n, e∞ = fieldvalues(eos)
 
     x = v / v0
     η = n + 1
@@ -437,79 +453,79 @@ julia> map(f, 1:1:10)
 """
 apply(::PressureForm, eos::EquationOfState) = v -> apply(PressureForm(), eos, v)
 """
-    apply(PressureForm(), eos::Murnaghan, v::Real)
+    apply(PressureForm(), eos::Murnaghan, v)
 
 Return the pressure of a `Murnaghan` equation of state on volume `v`.
 """
-function apply(::PressureForm, eos::Murnaghan, v::Real)
-    v0, b0, bp0 = collect(eos)
+function apply(::PressureForm, eos::Murnaghan, v)
+    v0, b0, bp0 = fieldvalues(eos)
 
     return b0 / bp0 * ((v0 / v)^bp0 - 1)
 end
 """
-    apply(PressureForm(), eos::BirchMurnaghan2nd, v::Real)
+    apply(PressureForm(), eos::BirchMurnaghan2nd, v)
 
 Return the pressure of a `BirchMurnaghan2nd` equation of state on volume `v`.
 """
-function apply(::PressureForm, eos::BirchMurnaghan2nd, v::Real)
-    v0, b0 = collect(eos)
+function apply(::PressureForm, eos::BirchMurnaghan2nd, v)
+    v0, b0 = fieldvalues(eos)
 
     f = ((v0 / v)^(2 / 3) - 1) / 2
     return 3b0 * f * (1 + 2f)^(5 / 2)
 end
 """
-    apply(PressureForm(), eos::BirchMurnaghan3rd, v::Real)
+    apply(PressureForm(), eos::BirchMurnaghan3rd, v)
 
 Return the pressure of a `BirchMurnaghan3rd` equation of state on volume `v`.
 """
-function apply(::PressureForm, eos::BirchMurnaghan3rd, v::Real)
-    v0, b0, bp0 = collect(eos)
+function apply(::PressureForm, eos::BirchMurnaghan3rd, v)
+    v0, b0, bp0 = fieldvalues(eos)
 
     eta = (v0 / v)^(1 / 3)
     return 3 / 2 * b0 * (eta^7 - eta^5) * (1 + 3 / 4 * (bp0 - 4) * (eta^2 - 1))
 end
 """
-    apply(PressureForm(), eos::BirchMurnaghan4th, v::Real)
+    apply(PressureForm(), eos::BirchMurnaghan4th, v)
 
 Return the pressure of a `BirchMurnaghan4th` equation of state on volume `v`.
 """
-function apply(::PressureForm, eos::BirchMurnaghan4th, v::Real)
-    v0, b0, bp0, bpp0 = collect(eos)
+function apply(::PressureForm, eos::BirchMurnaghan4th, v)
+    v0, b0, bp0, bpp0 = fieldvalues(eos)
 
     f = ((v0 / v)^(2 / 3) - 1) / 2
     h = b0 * bpp0 + bp0^2
     return b0 / 2 * (2f + 1)^(5 / 2) * ((9h - 63bp0 + 143) * f^2 + 9 * (bp0 - 4) * f + 6)
 end
 """
-    apply(PressureForm(), eos::PoirierTarantola2nd, v::Real)
+    apply(PressureForm(), eos::PoirierTarantola2nd, v)
 
 Return the pressure of a `PoirierTarantola2nd` equation of state on volume `v`.
 """
-function apply(::PressureForm, eos::PoirierTarantola2nd, v::Real)
-    v0, b0 = collect(eos)
+function apply(::PressureForm, eos::PoirierTarantola2nd, v)
+    v0, b0 = fieldvalues(eos)
 
     x = (v / v0)^(1 / 3)
     return -b0 / x * log(x)
 end
 """
-    apply(PressureForm(), eos::PoirierTarantola3rd, v::Real)
+    apply(PressureForm(), eos::PoirierTarantola3rd, v)
 
 Return the pressure of a `PoirierTarantola3rd` equation of state on volume `v`.
 """
-function apply(::PressureForm, eos::PoirierTarantola3rd, v::Real)
-    v0, b0, bp0 = collect(eos)
+function apply(::PressureForm, eos::PoirierTarantola3rd, v)
+    v0, b0, bp0 = fieldvalues(eos)
 
     x = v / v0
     xi = log(x)
     return -b0 * xi / 2x * ((bp0 - 2) * xi - 2)
 end
 """
-    apply(PressureForm(), eos::PoirierTarantola4th, v::Real)
+    apply(PressureForm(), eos::PoirierTarantola4th, v)
 
 Return the pressure of a `PoirierTarantola4th` equation of state on volume `v`.
 """
-function apply(::PressureForm, eos::PoirierTarantola4th, v::Real)
-    v0, b0, bp0, bpp0 = collect(eos)
+function apply(::PressureForm, eos::PoirierTarantola4th, v)
+    v0, b0, bp0, bpp0 = fieldvalues(eos)
 
     x = (v / v0)^(1 / 3)
     xi = log(x)
@@ -517,35 +533,35 @@ function apply(::PressureForm, eos::PoirierTarantola4th, v::Real)
     return -b0 * xi / 6 / x * ((h + 3bp0 + 3) * xi^2 + 3 * (bp0 + 6) * xi + 6)
 end
 """
-    apply(PressureForm(), eos::Vinet, v::Real)
+    apply(PressureForm(), eos::Vinet, v)
 
 Return the pressure of a `Vinet` equation of state on volume `v`.
 """
-function apply(::PressureForm, eos::Vinet, v::Real)
-    v0, b0, bp0 = collect(eos)
+function apply(::PressureForm, eos::Vinet, v)
+    v0, b0, bp0 = fieldvalues(eos)
 
     x = (v / v0)^(1 / 3)
     xi = 3 / 2 * (bp0 - 1)
     return 3b0 / x^2 * (1 - x) * exp(xi * (1 - x))
 end
 """
-    apply(PressureForm(), eos::AntonSchmidt, v::Real)
+    apply(PressureForm(), eos::AntonSchmidt, v)
 
 Return the pressure of a `AntonSchmidt` equation of state on volume `v`.
 """
-function apply(::PressureForm, eos::AntonSchmidt, v::Real)
-    v0, β, n = collect(eos)
+function apply(::PressureForm, eos::AntonSchmidt, v)
+    v0, β, n = fieldvalues(eos)
 
     x = v / v0
     return -β * x^n * log(x)
 end
 """
-    apply(PressureForm(), eos::BreenanStacey, v::Real)
+    apply(PressureForm(), eos::BreenanStacey, v)
 
 Return the pressure of a `BreenanStacey` equation of state on volume `v`.
 """
-function apply(::PressureForm, eos::BreenanStacey, v::Real)
-    v0, b0, γ0 = collect(eos)
+function apply(::PressureForm, eos::BreenanStacey, v)
+    v0, b0, γ0 = fieldvalues(eos)
 
     x = v0 / v
     return b0 / 2 / γ0 * x^(4 / 3) * (exp(2γ0 * (1 - x)) - 1)
@@ -583,34 +599,34 @@ julia> map(f, 1:1:10)
 """
 apply(::BulkModulusForm, eos::EquationOfState) = v -> apply(BulkModulusForm(), eos, v)
 """
-    apply(BulkModulusForm(), eos::BirchMurnaghan2nd, v::Real)
+    apply(BulkModulusForm(), eos::BirchMurnaghan2nd, v)
 
 Return the bulk modulus of a `BirchMurnaghan2nd` equation of state on volume `v`.
 """
-function apply(::BulkModulusForm, eos::BirchMurnaghan2nd, v::Real)
-    v0, b0 = collect(eos)
+function apply(::BulkModulusForm, eos::BirchMurnaghan2nd, v)
+    v0, b0 = fieldvalues(eos)
 
     f = ((v0 / v)^(2 / 3) - 1) / 2
     return b0 * (7f + 1) * (2f + 1)^(5 / 2)
 end
 """
-    apply(BulkModulusForm(), eos::BirchMurnaghan3rd, v::Real)
+    apply(BulkModulusForm(), eos::BirchMurnaghan3rd, v)
 
 Return the bulk modulus of a `BirchMurnaghan3rd` equation of state on volume `v`.
 """
-function apply(::BulkModulusForm, eos::BirchMurnaghan3rd, v::Real)
-    v0, b0, bp0 = collect(eos)
+function apply(::BulkModulusForm, eos::BirchMurnaghan3rd, v)
+    v0, b0, bp0 = fieldvalues(eos)
 
     f = ((v0 / v)^(2 / 3) - 1) / 2
     return b0 / 2 * (2f + 1)^(5 / 2) * ((27 * f^2 + 6f) * (bp0 - 4) - 4f + 2)
 end
 """
-    apply(BulkModulusForm(), eos::BirchMurnaghan4th, v::Real)
+    apply(BulkModulusForm(), eos::BirchMurnaghan4th, v)
 
 Return the bulk modulus of a `BirchMurnaghan4th` equation of state on volume `v`.
 """
-function apply(::BulkModulusForm, eos::BirchMurnaghan4th, v::Real)
-    v0, b0, bp0, bpp0 = collect(eos)
+function apply(::BulkModulusForm, eos::BirchMurnaghan4th, v)
+    v0, b0, bp0, bpp0 = fieldvalues(eos)
 
     f = ((v0 / v)^(2 / 3) - 1) / 2
     h = b0 * bpp0 + bp0^2
@@ -618,35 +634,35 @@ function apply(::BulkModulusForm, eos::BirchMurnaghan4th, v::Real)
            ((99h - 693bp0 + 1573) * f^3 + (27h - 108bp0 + 105) * f^2 + 6f * (3bp0 - 5) + 6)
 end
 """
-    apply(BulkModulusForm(), eos::PoirierTarantola2nd, v::Real)
+    apply(BulkModulusForm(), eos::PoirierTarantola2nd, v)
 
 Return the bulk modulus of a `PoirierTarantola2nd` equation of state on volume `v`.
 """
-function apply(::BulkModulusForm, eos::PoirierTarantola2nd, v::Real)
-    v0, b0 = collect(eos)
+function apply(::BulkModulusForm, eos::PoirierTarantola2nd, v)
+    v0, b0 = fieldvalues(eos)
 
     x = (v / v0)^(1 / 3)
     return b0 / x * (1 - log(x))
 end
 """
-    apply(BulkModulusForm(), eos::PoirierTarantola3rd, v::Real)
+    apply(BulkModulusForm(), eos::PoirierTarantola3rd, v)
 
 Return the bulk modulus of a `PoirierTarantola3rd` equation of state on volume `v`.
 """
-function apply(::BulkModulusForm, eos::PoirierTarantola3rd, v::Real)
-    v0, b0, bp0 = collect(eos)
+function apply(::BulkModulusForm, eos::PoirierTarantola3rd, v)
+    v0, b0, bp0 = fieldvalues(eos)
 
     x = v / v0
     xi = log(x)
     return -b0 / 2x * (((bp0 - 2) * xi + 2 - 2bp0) * xi + 2)
 end
 """
-    apply(BulkModulusForm(), eos::PoirierTarantola4th, v::Real)
+    apply(BulkModulusForm(), eos::PoirierTarantola4th, v)
 
 Return the bulk modulus of a `PoirierTarantola4th` equation of state on volume `v`.
 """
-function apply(::BulkModulusForm, eos::PoirierTarantola4th, v::Real)
-    v0, b0, bp0, bpp0 = collect(eos)
+function apply(::BulkModulusForm, eos::PoirierTarantola4th, v)
+    v0, b0, bp0, bpp0 = fieldvalues(eos)
 
     x = (v / v0)^(1 / 3)
     xi = log(x)
@@ -655,24 +671,24 @@ function apply(::BulkModulusForm, eos::PoirierTarantola4th, v::Real)
            ((h + 3bp0 + 3) * xi^3 - 3 * xi^2 * (h + 2bp0 + 1) - 6xi * (bp0 + 1) - 6)
 end
 """
-    apply(BulkModulusForm(), eos::Vinet, v::Real)
+    apply(BulkModulusForm(), eos::Vinet, v)
 
 Return the bulk modulus of a `Vinet` equation of state on volume `v`.
 """
-function apply(::BulkModulusForm, eos::Vinet, v::Real)
-    v0, b0, bp0 = collect(eos)
+function apply(::BulkModulusForm, eos::Vinet, v)
+    v0, b0, bp0 = fieldvalues(eos)
 
     x = (v / v0)^(1 / 3)
     xi = 3 / 2 * (bp0 - 1)
     return -b0 / (2 * x^2) * (3x * (x - 1) * (bp0 - 1) + 2 * (x - 2)) * exp(-xi * (x - 1))
 end
 """
-    apply(BulkModulusForm(), eos::AntonSchmidt, v::Real)
+    apply(BulkModulusForm(), eos::AntonSchmidt, v)
 
 Return the bulk modulus of a `AntonSchmidt` equation of state on volume `v`.
 """
-function apply(::BulkModulusForm, eos::AntonSchmidt, v::Real)
-    v0, β, n = collect(eos)
+function apply(::BulkModulusForm, eos::AntonSchmidt, v)
+    v0, β, n = fieldvalues(eos)
 
     x = v / v0
     return β * x^n * (1 + n * log(x))
@@ -683,32 +699,20 @@ end
 # ============================================================================ #
 #                                 Miscellaneous                                #
 # ============================================================================ #
-function allsubtypes(T::Type, types = Type[])::Vector{Type}
-    for S in subtypes(T)
-        types = allsubtypes(S, push!(types, S))
-    end
-    types
-end
-
-nonabstract(T::Type)::Vector{Type} = filter(!isabstracttype, allsubtypes(T))
-
-for E in nonabstract(EquationOfState)
-    eval(quote
-        similar_type(::Type{A}, ::Type{T}) where {A<:$E,T} = $E{T}
-    end)
-end
-
-Base.iterate(eos::EquationOfState) = getfield(eos, 1), 2
-function Base.iterate(eos::EquationOfState, state::Integer)
-    state > length(eos) && return nothing
-    return getfield(eos, state), state + 1
-end # function Base.iterate
-Base.length(eos::EquationOfState) = nfields(eos)
-Base.size(eos::EquationOfState) = (length(eos),)
-Base.eltype(::Type{<:EquationOfState{T}}) where {T} = T
-
-Base.isapprox(a::T, b::T; kwargs...) where {T<:EquationOfState} = isapprox(collect(a), collect(b); kwargs...)
-# Base.getindex(eos::PolynomialEquationOfState{T,N}, index::Int64) where {T,N} = getindex(eos.data, index)
+# This is a helper function and should not be exported.
+fieldvalues(eos::EquationOfState) = [getfield(eos, i) for i in 1:nfields(eos)]
+
+Unitful.upreferred(::Dimensions{(
+    Dimension{:Length}(2 // 1),
+    Dimension{:Mass}(1 // 1),
+    Dimension{:Time}(-2 // 1),
+)}) = u"eV"
+Unitful.upreferred(::Dimensions{(Dimension{:Length}(3 // 1),)}) = u"angstrom^3"
+Unitful.upreferred(::Dimensions{(
+    Dimension{:Length}(-1 // 1),
+    Dimension{:Mass}(1 // 1),
+    Dimension{:Time}(-2 // 1),
+)}) = u"eV/angstrom^3"
 # =============================== Miscellaneous ============================== #
 
 end
diff --git a/src/Find.jl b/src/Find.jl
index dfd2317..b5f8579 100644
--- a/src/Find.jl
+++ b/src/Find.jl
@@ -12,7 +12,7 @@ julia>
 module Find
 
 using InteractiveUtils: subtypes
-using Statistics: median
+using Unitful: AbstractQuantity, ustrip
 
 using Roots: find_zero,
              AbstractBracketing,
@@ -31,35 +31,22 @@ using ..Collections: EquationOfState, apply
 
 export findvolume
 
-function findvolume(
-    form::EquationForm,
-    eos::EquationOfState,
-    y::Real,
-    domain::Union{AbstractVector,Tuple},
-    method::AbstractBracketing,
-)
-    f(v) = apply(form, eos, v) - y
-    return find_zero(f, (minimum(domain), maximum(domain)), method)
+function findvolume(form::EquationForm, eos::EquationOfState, y, x0, method)
+    f = v -> apply(form, eos, v) - y
+    return find_zero(f, x0, method)
 end # function findvolume
-function findvolume(
-    form::EquationForm,
-    eos::EquationOfState,
-    y::Real,
-    domain::Union{AbstractVector,Tuple},
-    method::Union{AbstractNonBracketing,AbstractHalleyLikeMethod,AbstractNewtonLikeMethod},
-)
-    f(v) = apply(form, eos, v) - y
-    return find_zero(f, median(domain), method)
-end # function findvolume
-function findvolume(
-    form::EquationForm,
-    eos::EquationOfState,
-    y::Real,
-    domain::Union{AbstractVector,Tuple},
-)
+function findvolume(form::EquationForm, eos::EquationOfState, y, x0::Union{AbstractVector,Tuple})
+    for T in [subtypes(AbstractBisection); subtypes(AbstractAlefeldPotraShi)]
+        @info("Using method \"$T\"...")
+        try
+            # `maximum` and `minimum` also works with `AbstractQuantity`s.
+            return findvolume(form, eos, y, (minimum(x0), maximum(x0)), T())
+        catch e
+            @info("Method \"$T\" failed because of $e.")
+            continue
+        end
+    end
     for T in [
-        subtypes(AbstractAlefeldPotraShi)
-        subtypes(AbstractBisection)
         Brent
         subtypes(AbstractHalleyLikeMethod)
         Newton
@@ -67,7 +54,7 @@ function findvolume(
     ]
         @info("Using method \"$T\"...")
         try
-            return findvolume(form, eos, y, domain, T())
+            return findvolume(form, eos, y, (minimum(x0) + maximum(x0)) / 2, T())
         catch e
             @info("Method \"$T\" failed because of $e.")
             continue
diff --git a/src/NonlinearFitting.jl b/src/NonlinearFitting.jl
index 260a434..2acbb24 100644
--- a/src/NonlinearFitting.jl
+++ b/src/NonlinearFitting.jl
@@ -11,7 +11,9 @@ julia>
 """
 module NonlinearFitting
 
+using ConstructionBase: constructorof
 using LsqFit: curve_fit
+using Unitful: AbstractQuantity, upreferred, ustrip, unit
 
 import ..EquationForm
 using ..Collections
@@ -33,17 +35,43 @@ Fit an equation of state using least-squares fitting method (with the Levenberg-
 """
 function lsqfit(
     form::EquationForm,
-    eos::E,
-    xdata::AbstractVector,
-    ydata::AbstractVector;
+    eos::EquationOfState{<:Real},
+    xdata::AbstractVector{<:Real},
+    ydata::AbstractVector{<:Real};
     debug = false,
-    kwargs...
-) where {E<:EquationOfState}
-    T = promote_type(eltype(eos), eltype(xdata), eltype(ydata), Float64)
-    P = Collections.similar_type(E, T)
-    model(x, p) = map(apply(form, P(p...)), x)
-    fitted = curve_fit(model, T.(xdata), T.(ydata), T.(collect(eos)); kwargs...)
-    return debug ? fitted : P(fitted.param...)
+    kwargs...,
+)
+    T = promote_type(eltype(xdata), eltype(ydata), Float64)
+    E = constructorof(typeof(eos))
+    model = (x, p) -> map(apply(form, E(p...)), x)
+    fitted = curve_fit(
+        model,
+        T.(xdata),
+        T.(ydata),
+        T.(Collections.fieldvalues(eos)),
+        kwargs...,
+    )
+    return debug ? fitted : E(fitted.param...)
+end  # function lsqfit
+function lsqfit(
+    form::EquationForm,
+    eos::EquationOfState{<:AbstractQuantity},
+    xdata::AbstractVector{<:AbstractQuantity},
+    ydata::AbstractVector{<:AbstractQuantity};
+    kwargs...,
+)
+    E = constructorof(typeof(eos))
+    values = Collections.fieldvalues(eos)
+    original_units = map(unit, values)
+    trial_params, xdata, ydata = [map(ustrip ∘ upreferred, x) for x in (values, xdata, ydata)]
+    result = lsqfit(form, E(trial_params...), xdata, ydata; kwargs...)
+    if result isa EquationOfState
+        data = Collections.fieldvalues(result)
+        return E(
+            [data[i] * upreferred(u) |> u for (i, u) in enumerate(original_units)]...
+        )
+    end
+    return result
 end  # function lsqfit
 
 end
diff --git a/test/Collections.jl b/test/Collections.jl
index 51eabb9..9d09396 100644
--- a/test/Collections.jl
+++ b/test/Collections.jl
@@ -1,5 +1,7 @@
 using Test
 
+using Unitful
+
 using EquationsOfState.Collections
 
 @testset "Test EOS promotion" begin
@@ -13,6 +15,10 @@ using EquationsOfState.Collections
     @test typeof(PoirierTarantola4th(1, 2, 3, 4, 0)) == PoirierTarantola4th{Int}
     @test typeof(AntonSchmidt(1, 2, 3.0, 0)) == AntonSchmidt{Float64}
     @test typeof(BreenanStacey(1, 2, 3.0, 0)) == BreenanStacey{Float64}
+    @test Murnaghan(1, Int32(2), Int8(3), 0) == Murnaghan{Int64}(1, 2, 3, 0)
+    @test Murnaghan(1, 2//1, Int8(3), 0) == Murnaghan{Rational{Int64}}(1//1, 2//1, 3//1, 0//1)
+    @test typeof(Murnaghan(1u"angstrom^3", 2u"eV/angstrom^3", 3.0, 4u"eV")) == Murnaghan{Quantity{Float64}}
+    @test typeof(Murnaghan(1u"angstrom^3", 2u"eV/angstrom^3", 3//2, 4u"eV")) == Murnaghan{Quantity{Rational{Int64}}}
 end
 
 @testset "Test default EOS parameter `e0` and promotion" begin
@@ -26,4 +32,7 @@ end
     @test PoirierTarantola4th(1, 2, 3, 4) == PoirierTarantola4th(1, 2, 3, 4, 0)
     @test AntonSchmidt(1, 2, 3.0) == AntonSchmidt(1.0, 2.0, 3.0, 0.0)
     @test BreenanStacey(1, 2, 3.0) == BreenanStacey(1.0, 2.0, 3.0, 0.0)
+    @test typeof(Murnaghan(1u"angstrom^3", 2u"eV/angstrom^3", 3)) == Murnaghan{Quantity{Int64}}
+    @test typeof(Murnaghan(1u"angstrom^3", 2u"eV/angstrom^3", 3.0)) == Murnaghan{Quantity{Float64}}
+    @test typeof(Murnaghan(1.0u"angstrom^3", 2u"eV/angstrom^3", 3.0)) == Murnaghan{Quantity{Float64}}
 end
diff --git a/test/Find.jl b/test/Find.jl
index 567109e..ac99aca 100644
--- a/test/Find.jl
+++ b/test/Find.jl
@@ -1,6 +1,8 @@
 using Test
 
 using Roots
+using Unitful
+using UnitfulAtomic
 
 using EquationsOfState
 using EquationsOfState.Collections
@@ -38,7 +40,7 @@ using EquationsOfState.Find
         55.8775030703,
         57.4014349722,
         58.9526328669,
-    ]
+    ] .* u"angstrom^3"
     energies = [
         -7.63622156576,
         -8.16831294894,
@@ -68,23 +70,36 @@ using EquationsOfState.Find
         -9.99504030111,
         -9.86535084973,
         -9.73155247952,
-    ]
-    @test isapprox(
+    ] .* u"eV"
+    isapprox(
         map(
             e -> findvolume(
                 EnergyForm(),
                 BirchMurnaghan3rd(
-                    40.98926572528106,
-                    0.5369258245417454,
-                    4.178644235500821,
-                    -10.842803908240892,
+                    40.98926572528106u"angstrom^3",
+                    0.5369258245417454u"eV/angstrom^3",
+                    4.178644235500821u"1000mm/m",
+                    -10.842803908240892u"eV",
                 ),
                 e,
-                (25, 60),
-                Order16(),
+                (eps(), 100) .* u"angstrom^3",
             ),
             energies,
         ),
         results,
     )
 end
+
+@testset "Test `findvolume` with random unit" begin
+    pressures = collect(0:20:200) .* u"GPa"
+    eos = BirchMurnaghan3rd(167u"angstrom^3", 2600u"kbar", 4.0u"1000mm/m")
+    volumes = map(
+        p -> findvolume(PressureForm(), eos, p, (eps() * u"bohr^3", eos.v0 * 1.3)),
+        pressures,
+    )
+    @test isapprox(
+        ustrip.(map(apply(PressureForm(), eos), volumes) - pressures),
+        zeros(11),
+        atol = 1e-5,
+    )
+end # testset
diff --git a/test/NonlinearFitting.jl b/test/NonlinearFitting.jl
index 9b4861f..3d86801 100644
--- a/test/NonlinearFitting.jl
+++ b/test/NonlinearFitting.jl
@@ -10,14 +10,14 @@ using EquationsOfState.NonlinearFitting
         103.58772269057364,
         -144.45152457521132,
         -40.31992619868024,
-    )
+    ) |> Collections.fieldvalues
     @test isapprox(
         lsqfit(
             EnergyForm(),
             BirchMurnaghan3rd(1, 2, 3.0, 0),
             [1, 2, 3, 4, 5],
             [5, 6, 9, 8, 7],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-5,
     )
@@ -27,7 +27,7 @@ using EquationsOfState.NonlinearFitting
             BirchMurnaghan3rd(1, 2, 3, 0),
             [1, 2, 3, 4, 5.0],
             [5, 6, 9, 8, 7],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-5,
     )
@@ -37,7 +37,7 @@ using EquationsOfState.NonlinearFitting
             BirchMurnaghan3rd(1, 2, 3.0, 0),
             [1, 2, 3, 4, 5],
             [5, 6, 9, 8, 7.0],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-5,
     )
@@ -47,7 +47,7 @@ using EquationsOfState.NonlinearFitting
             BirchMurnaghan3rd(1, 2, 3, 0),
             [1, 2, 3, 4, 5],
             [5, 6, 9, 8, 7],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-5,
     )
@@ -59,14 +59,14 @@ end
         29.30861698140365,
         12.689089871112746,
         0.0,
-    )
+    ) |> Collections.fieldvalues
     @test isapprox(
         lsqfit(
             PressureForm(),
             BirchMurnaghan3rd(1, 2, 3.0, 0),
             [1, 2, 3, 4, 5],
             [5, 6, 9, 8, 7],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-6,
     )
@@ -76,7 +76,7 @@ end
             BirchMurnaghan3rd(1, 2, 3, 0),
             [1, 2, 3, 4, 5.0],
             [5, 6, 9, 8, 7],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-6,
     )
@@ -86,7 +86,7 @@ end
             BirchMurnaghan3rd(1, 2, 3.0, 0),
             [1, 2, 3, 4, 5],
             [5, 6, 9, 8, 7.0],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-6,
     )
@@ -96,21 +96,21 @@ end
             BirchMurnaghan3rd(1, 2, 3, 0),
             [1, 2, 3, 4, 5],
             [5, 6, 9, 8, 7],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-6,
     )
 end
 
 @testset "Test fitting bulk modulus with different element types" begin
-    result = BirchMurnaghan3rd(7.218928431312577, 5.007900469653902, 4.06037725509478, 0.0)
+    result = BirchMurnaghan3rd(7.218928431312577, 5.007900469653902, 4.06037725509478, 0.0) |> Collections.fieldvalues
     @test isapprox(
         lsqfit(
             BulkModulusForm(),
             BirchMurnaghan3rd(1, 2, 3.0, 0),
             [1, 2, 3, 4, 5],
             [5, 6, 9, 8, 7],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-5,
     )
@@ -120,7 +120,7 @@ end
             BirchMurnaghan3rd(1, 2, 3, 0),
             [1, 2, 3, 4, 5.0],
             [5, 6, 9, 8, 7],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-5,
     )
@@ -130,7 +130,7 @@ end
             BirchMurnaghan3rd(1, 2, 3.0, 0),
             [1, 2, 3, 4, 5],
             [5, 6, 9, 8, 7.0],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-5,
     )
@@ -140,7 +140,7 @@ end
             BirchMurnaghan3rd(1, 2, 3, 0),
             [1, 2, 3, 4, 5],
             [5, 6, 9, 8, 7],
-        ),
+        ) |> Collections.fieldvalues,
         result;
         atol = 1e-5,
     )
@@ -210,40 +210,40 @@ end
         -9.73155247952,
     ]
     @test isapprox(
-        lsqfit(EnergyForm(), BirchMurnaghan3rd(40, 0.5, 4, 0), volumes, energies),
+        lsqfit(EnergyForm(), BirchMurnaghan3rd(40, 0.5, 4, 0), volumes, energies) |> Collections.fieldvalues,
         BirchMurnaghan3rd(
             40.98926572528106,
             0.5369258245417454,
             4.178644235500821,
             -10.842803908240892,
-        ),
+        ) |> Collections.fieldvalues,
     )
     @test isapprox(
-        lsqfit(EnergyForm(), Murnaghan(41, 0.5, 4, 0), volumes, energies),
+        lsqfit(EnergyForm(), Murnaghan(41, 0.5, 4, 0), volumes, energies) |> Collections.fieldvalues,
         Murnaghan(
             41.13757930387086,
             0.5144967693786603,
             3.9123862262572264,
             -10.836794514626673,
-        ),
+        ) |> Collections.fieldvalues,
     )
     @test isapprox(
-        lsqfit(EnergyForm(), PoirierTarantola3rd(41, 0.5, 4, 0), volumes, energies),
+        lsqfit(EnergyForm(), PoirierTarantola3rd(41, 0.5, 4, 0), volumes, energies) |> Collections.fieldvalues,
         PoirierTarantola3rd(
             40.86770643373908,
             0.5667729960804602,
             4.331688936974368,
             -10.851486685041658,
-        ),
+        ) |> Collections.fieldvalues,
     )
     @test isapprox(
-        lsqfit(EnergyForm(), Vinet(41, 0.5, 4, 0), volumes, energies),
+        lsqfit(EnergyForm(), Vinet(41, 0.5, 4, 0), volumes, energies) |> Collections.fieldvalues,
         Vinet(
             40.916875663779784,
             0.5493839425156859,
             4.3051929654936885,
             -10.846160810560756,
-        ),
+        ) |> Collections.fieldvalues,
     )
     # 'deltafactor': {'b0': 0.5369258245611414,
 #             'b1': 4.178644231924639,
@@ -331,13 +331,13 @@ end
 
     fitted_eos = lsqfit(EnergyForm(), Vinet(23, 0.5, 4, -2), mp153_volumes, mp153_energies)
     @test isapprox(
-        fitted_eos,
+        fitted_eos |> Collections.fieldvalues,
         Vinet(
             22.95764559358769,
             0.2257091141420788,
             4.060543387224629,
             -1.5944292606251582,
-        ),
+        ) |> Collections.fieldvalues,
     )
     @test isapprox(
         map(apply(EnergyForm(), fitted_eos), mp153_volumes),
@@ -421,13 +421,13 @@ end
 
     fitted_eos = lsqfit(EnergyForm(), Vinet(20, 0.5, 4, -5), mp149_volumes, mp149_energies)
     @test isapprox(
-        fitted_eos,
+        fitted_eos |> Collections.fieldvalues,
         Vinet(
             20.446696754873944,
             0.5516638521306302,
             4.324373909783161,
             -5.424963389876503,
-        ),
+        ) |> Collections.fieldvalues,
     )
     @test isapprox(
         map(apply(EnergyForm(), fitted_eos), mp149_volumes),
@@ -511,13 +511,13 @@ end
 
     fitted_eos = lsqfit(EnergyForm(), Vinet(17, 0.5, 4, -7), mp72_volumes, mp72_energies)
     @test isapprox(
-        fitted_eos,
+        fitted_eos |> Collections.fieldvalues,
         Vinet(
             17.13223026131245,
             0.7029766224730147,
             3.6388077563621812,
             -7.897414959124461,
-        ),
+        ) |> Collections.fieldvalues,
     )
     @test isapprox(
         map(apply(EnergyForm(), fitted_eos), mp72_volumes),
@@ -597,26 +597,26 @@ end
         BirchMurnaghan3rd(224, 0.0006, 4, -323),
         volumes,
         energies,
-    ) ≈ BirchMurnaghan3rd(224.444565, 0.00062506191050572675, 3.740369, -323.417714)
+    ) |> Collections.fieldvalues ≈ BirchMurnaghan3rd(224.444565, 0.00062506191050572675, 3.740369, -323.417714) |> Collections.fieldvalues
     @test isapprox(
         lsqfit(
             EnergyForm(),
             BirchMurnaghan4th(224, 0.0006, 4, -5460, -323),
             volumes,
             energies,
-        ),
+        ) |> Collections.fieldvalues,
         BirchMurnaghan4th(
             224.45756247137314,
             0.0006229382259822287,
             3.730991473426449,
             -5322.693307704408,
             -323.4177113158418,
-        );
-        atol = 1e-3,
+        ) |> Collections.fieldvalues;
+        rtol = 1e-3,
     )
     @test isapprox(
-        lsqfit(EnergyForm(), Murnaghan(224, 0.006, 4, -323), volumes, energies),
-        Murnaghan(224.501825, 0.00060479524074699499, 3.723835, -323.417686);
+        lsqfit(EnergyForm(), Murnaghan(224, 0.006, 4, -323), volumes, energies) |> Collections.fieldvalues,
+        Murnaghan(224.501825, 0.00060479524074699499, 3.723835, -323.417686) |> Collections.fieldvalues;
         atol = 1e-5,
     )
     @test isapprox(
@@ -625,8 +625,8 @@ end
             PoirierTarantola3rd(100, 0.0006, 3.7, -323),
             volumes,
             energies,
-        ),
-        PoirierTarantola3rd(224.509208, 0.000635892264159838, 3.690448, -323.41773);
+        ) |> Collections.fieldvalues,
+        PoirierTarantola3rd(224.509208, 0.000635892264159838, 3.690448, -323.41773) |> Collections.fieldvalues;
         atol = 1e-5,
     )
     # @test lsqfit(EnergyForm(), PoirierTarantola4th(220, 0.0006, 3.7, -5500, -323), volumes, energies; lower = Float64[220, 0, 3, -6000, -400], upper = Float64[300, 0.01, 5, -5000, -300]) ≈ PoirierTarantola4th(224.430182, 0.0006232241765069493, 3.758360, -5493.859729817176, -323.417712)