diff --git a/Manifest.toml b/Manifest.toml
index d4b4bec..5b1bf94 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -41,9 +41,9 @@ 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"
+git-tree-sha1 = "ed2c4abadf84c53d9e58510b5fc48912c2336fbb"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "2.1.0"
+version = "2.2.0"
 
 [[ConstructionBase]]
 git-tree-sha1 = "e3efe0a0f49dcd294c8c73e897b4fdf891f0fbd3"
@@ -179,10 +179,10 @@ deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUID
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
 [[Polynomials]]
-deps = ["LinearAlgebra", "SparseArrays", "Test"]
-git-tree-sha1 = "62142bd65d3f8aeb2226ec64dd8493349147df94"
+deps = ["LinearAlgebra", "RecipesBase"]
+git-tree-sha1 = "f7c0c07e82798aef542d60a6e6e85e39f4590750"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
-version = "0.5.2"
+version = "0.5.3"
 
 [[Printf]]
 deps = ["Unicode"]
@@ -202,6 +202,11 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 deps = ["Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
+[[RecipesBase]]
+git-tree-sha1 = "7bdce29bc9b2f5660a6e5e64d64d91ec941f6aa2"
+uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+version = "0.7.0"
+
 [[Reexport]]
 deps = ["Pkg"]
 git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
@@ -302,6 +307,6 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 deps = ["LinearAlgebra", "Random"]
 git-tree-sha1 = "0e9ad8a9621a9e1ec48f7cb48a5d6ded9e2c68f3"
 repo-rev = "master"
-repo-url = "https://github.com/PainterQubits/Unitful.jl.git"
+repo-url = "https://github.com/PainterQubits/Unitful.jl"
 uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
 version = "0.17.0"
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
new file mode 100644
index 0000000..9b6e197
--- /dev/null
+++ b/docs/Manifest.toml
@@ -0,0 +1,104 @@
+# This file is machine-generated - editing it directly is not advised
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[DocStringExtensions]]
+deps = ["LibGit2", "Markdown", "Pkg", "Test"]
+git-tree-sha1 = "88bb0edb352b16608036faadcc071adda068582a"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.8.1"
+
+[[Documenter]]
+deps = ["Base64", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
+git-tree-sha1 = "d45c163c7a3ae293c15361acc52882c0f853f97c"
+uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+version = "0.23.4"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[JSON]]
+deps = ["Dates", "Mmap", "Parsers", "Unicode"]
+git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.21.0"
+
+[[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"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[Parsers]]
+deps = ["Dates", "Test"]
+git-tree-sha1 = "ef0af6c8601db18c282d092ccbd2f01f3f0cd70b"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "0.3.7"
+
+[[Pkg]]
+deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[Test]]
+deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[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/docs/Project.toml b/docs/Project.toml
index dfa65cd..2c28345 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,2 +1,3 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
diff --git a/docs/make.jl b/docs/make.jl
index c925022..2f53032 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,13 +1,16 @@
 using Documenter, EquationsOfState
 
+DocMeta.setdocmeta!(EquationsOfState, :DocTestSetup, :(using EquationsOfState, EquationOfState.Collections, Unitful); recursive=true)
 makedocs(;
     modules=[EquationsOfState],
     format=Documenter.HTML(),
     pages=[
         "Home" => "index.md",
+        "Installation" => "Installation.md",
         "Manual" => [
             "Collections" => "Collections.md",
-            "Nonlinear fitting" => "NonlinearFitting.md"
+            "Nonlinear fitting" => "NonlinearFitting.md",
+            "Find volume" => "Find.md"
         ]
     ],
     repo="https://github.com/MineralsCloud/EquationsOfState.jl/blob/{commit}{path}#L{line}",
diff --git a/docs/src/Collections.md b/docs/src/Collections.md
index 62c2d6e..663045a 100644
--- a/docs/src/Collections.md
+++ b/docs/src/Collections.md
@@ -18,94 +18,67 @@ EquationOfState
 │  ├─ PoirierTarantola3rd
 │  └─ PoirierTarantola4th
 ├─ Murnaghan
-├─ PolynomialEquationOfState
 └─ Vinet
 ```
 
-## Types
-
-```@docs
-EquationOfState
-FiniteStrainEquationOfState
-Murnaghan
-BirchMurnaghan2nd
-BirchMurnaghan3rd
-BirchMurnaghan4th
-PoirierTarantola2nd
-PoirierTarantola3rd
-PoirierTarantola4th
-Vinet
-```
-
 ## Usage
 
 ### Construct an `EquationOfState`
+
 We will use `BirchMurnaghan3rd` as an example.
 
-`BirchMurnaghan3rd` can be constructed from scratch:
+A `BirchMurnaghan3rd` can be constructed from scratch, as shown above. It can
+also be constructed from an existing `BirchMurnaghan3rd`, with
+[`Setfield.jl`](https://github.com/jw3126/Setfield.jl)
+[`@set!`](https://jw3126.github.io/Setfield.jl/stable/#Setfield.@set!-Tuple{Any})
+macro:
 
 ```julia
-julia> BirchMurnaghan3rd(1, 2, 3)
-4-element BirchMurnaghan3rd{Int64}:
- 1
- 2
- 3
- 0
-
-julia> BirchMurnaghan3rd(1, 2, 3, 4)
-4-element BirchMurnaghan3rd{Int64}:
- 1
- 2
- 3
- 4
-
-julia> BirchMurnaghan3rd(1, 2, 3, 4.0)
-4-element BirchMurnaghan3rd{Float64}:
- 1.0
- 2.0
- 3.0
- 4.0
+julia> using Setfield
+
+julia> eos = Murnaghan(1, 2, 3.0)
+Murnaghan{Float64}(1.0, 2.0, 3.0, 0.0)
+
+julia> @set! eos.v0 = 4
+Murnaghan{Float64}(4.0, 2.0, 3.0, 0.0)
+
+julia> eos
+Murnaghan{Float64}(4.0, 2.0, 3.0, 0.0)
 ```
 
-It can also be constructed from an existing `BirchMurnaghan3rd`:
+To modify multiple fields (say, `:v0`, `:bp0`, `:bpp0`, `:e0`) at a time, use
+[`@batchlens`](https://tkf.github.io/Kaleido.jl/stable/#Kaleido.@batchlens) from
+[`Kaleido.jl`](https://github.com/tkf/Kaleido.jl):
 
 ```julia
-julia> BirchMurnaghan3rd(BirchMurnaghan3rd(1, 2, 3, 4.0), b0=10, e0=5)
-4-element BirchMurnaghan3rd{Float64}:
-  1.0
- 10.0
-  3.0
-  5.0
+julia> using Setfield, Kaleido
 
-julia> BirchMurnaghan3rd(BirchMurnaghan3rd(1, 2, 3, 4.0), Dict(:b0=>10, :e0=>5))
-4-element BirchMurnaghan3rd{Float64}:
-  1.0
- 10.0
-  3.0
-  5.0
+julia> lens = @batchlens(begin
+           _.v0
+           _.bp0
+           _.bpp0
+           _.e0
+       end)
+IndexBatchLens(:v0, :bp0, :bpp0, :e0)
 
-julia> BirchMurnaghan3rd(BirchMurnaghan3rd(1, 2, 3, 4.0), (:b0, 10))
-4-element BirchMurnaghan3rd{Float64}:
-  1.0
- 10.0
-  3.0
-  4.0
+julia> eos = BirchMurnaghan4th(1, 2.0, 3, 4)
+BirchMurnaghan4th{Float64}(1.0, 2.0, 3.0, 4.0, 0.0)
+
+julia> set(eos, lens, (5, 6, 7, 8))
+BirchMurnaghan4th{Float64}(5.0, 2.0, 6.0, 7.0, 8.0)
 ```
 
-Users can access `BirchMurnaghan3rd`'s element by either "dot notation" or indexing:
+Users can access `BirchMurnaghan3rd`'s elements by "dot notation":
 
 ```julia
-julia> b = BirchMurnaghan3rd(1, 2, 3, 4.0)
+julia> eos = BirchMurnaghan3rd(1, 2, 3, 4.0)
 4-element BirchMurnaghan3rd{Float64}:
  1.0
  2.0
  3.0
  4.0
 
-julia> b.v0
-1.0
-
-julia> b[1]
+julia> eos.v0
 1.0
 ```
 
@@ -114,6 +87,7 @@ julia> b[1]
 The $E(V)$ relation of equations of state are listed as below:
 
 1. `Murnaghan`:
+
    ```math
    E(V) = E_{0}+K_{0} V_{0}\left[\frac{1}{K_{0}^{\prime}\left(K_{0}^{\prime}-1\right)}\left(\frac{V}{V_{0}}\right)^{1-K_{0}^{\prime}}+\frac{1}{K_{0}^{\prime}} \frac{V}{V_{0}}-\frac{1}{K_{0}^{\prime}-1}\right].
    ```
@@ -130,7 +104,8 @@ The $E(V)$ relation of equations of state are listed as below:
    E(V) = E_{0}+\frac{9}{16} V_{0} B_{0} \frac{\left(x^{2 / 3}-1\right)^{2}}{x^{7 / 3}}\left\{x^{1 / 3}\left(B_{0}^{\prime}-4\right)-x\left(B_{0}^{\prime}-6\right)\right\}.
    ```
 
-   where ``x = V / V_0``, and ``f = \frac{ 1 }{ 2 } \bigg[ \bigg( \frac{ V_0 }{ V } \bigg)^{2/3} - 1 \bigg]``.
+   where ``x = V / V_0``, and
+   ``f = \frac{ 1 }{ 2 } \bigg[ \bigg( \frac{ V_0 }{ V } \bigg)^{2/3} - 1 \bigg]``.
 
 4. `BirchMurnaghan4th`:
 
@@ -157,6 +132,7 @@ The $E(V)$ relation of equations of state are listed as below:
    ```math
    E(V) = E_{0}+\frac{1}{24} B_{0} V_{0} \ln ^{2} x\left\{\left(H+3 B_{0}^{\prime}+3\right) \ln ^{2} x\right. \left.+4\left(B_{0}^{\prime}+2\right) \ln x+12\right\}.
    ```
+
    where ``H = B_0 B_0'' + (B_0')^2``.
 
 8. `Vinet`:
@@ -242,7 +218,7 @@ The $B(V)$ relation of equations of state are listed as below:
 2. `BirchMurnaghan3rd`:
 
    ```math
-   B(V) = \frac{B_{0}}{8 \chi^{10 / 3}}\left\{x^{5 / 3}\left(15 B_{0}^{\prime}-80\right)-x\left(42 B_{0}^{\prime}-196\right)\right.\left.+27 x^{1 / 3}\left(B_{0}^{\prime}-4\right)\right\}.
+   B(V) = \frac{B_{0}}{8 x^{10 / 3}}\left\{x^{5 / 3}\left(15 B_{0}^{\prime}-80\right)-x\left(42 B_{0}^{\prime}-196\right)\right.\left.+27 x^{1 / 3}\left(B_{0}^{\prime}-4\right)\right\}.
    ```
 
 3. `BirchMurnaghan4th`:
@@ -281,38 +257,26 @@ The $B(V)$ relation of equations of state are listed as below:
    B(V) = \beta\left(\frac{V}{V_{0}}\right)^{n}\left[1+n \ln \frac{V}{V_{0}}\right].
    ```
 
+## Types
+
+```@docs
+EquationOfState
+FiniteStrainEquationOfState
+Murnaghan
+BirchMurnaghan2nd
+BirchMurnaghan3rd
+BirchMurnaghan4th
+PoirierTarantola2nd
+PoirierTarantola3rd
+PoirierTarantola4th
+Vinet
+```
 
 ## Public interfaces
 
 ```@docs
 apply(::EnergyForm, eos::EquationOfState)
-apply(::EnergyForm, eos::Murnaghan, v::Real)
-apply(::EnergyForm, eos::BirchMurnaghan2nd, v::Real)
-apply(::EnergyForm, eos::BirchMurnaghan3rd, v::Real)
-apply(::EnergyForm, eos::BirchMurnaghan4th, v::Real)
-apply(::EnergyForm, eos::PoirierTarantola2nd, v::Real)
-apply(::EnergyForm, eos::PoirierTarantola3rd, v::Real)
-apply(::EnergyForm, eos::PoirierTarantola4th, v::Real)
-apply(::EnergyForm, eos::Vinet, v::Real)
-apply(::EnergyForm, eos::AntonSchmidt, v::Real)
-apply(::PressureForm, eos::EquationOfState)
-apply(::PressureForm, eos::Murnaghan, v::Real)
-apply(::PressureForm, eos::BirchMurnaghan2nd, v::Real)
-apply(::PressureForm, eos::BirchMurnaghan3rd, v::Real)
-apply(::PressureForm, eos::BirchMurnaghan4th, v::Real)
-apply(::PressureForm, eos::PoirierTarantola2nd, v::Real)
-apply(::PressureForm, eos::PoirierTarantola3rd, v::Real)
-apply(::PressureForm, eos::PoirierTarantola4th, v::Real)
-apply(::PressureForm, eos::Vinet, v::Real)
-apply(::PressureForm, eos::AntonSchmidt, v::Real)
-apply(::PressureForm, eos::BreenanStacey, v::Real)
-apply(::BulkModulusForm, eos::EquationOfState)
-apply(::BulkModulusForm, eos::BirchMurnaghan2nd, v::Real)
-apply(::BulkModulusForm, eos::BirchMurnaghan3rd, v::Real)
-apply(::BulkModulusForm, eos::BirchMurnaghan4th, v::Real)
-apply(::BulkModulusForm, eos::PoirierTarantola2nd, v::Real)
-apply(::BulkModulusForm, eos::PoirierTarantola3rd, v::Real)
-apply(::BulkModulusForm, eos::PoirierTarantola4th, v::Real)
-apply(::BulkModulusForm, eos::Vinet, v::Real)
-apply(::BulkModulusForm, eos::AntonSchmidt, v::Real)
+apply(::EnergyForm, eos::Murnaghan, v)
+apply(::PressureForm, eos::Murnaghan, v)
+apply(::BulkModulusForm, eos::BirchMurnaghan2nd, v)
 ```
diff --git a/docs/src/Find.md b/docs/src/Find.md
new file mode 100644
index 0000000..e9d479b
--- /dev/null
+++ b/docs/src/Find.md
@@ -0,0 +1,63 @@
+# Find
+
+```@meta
+CurrentModule = EquationsOfState.Find
+```
+
+This module contains a function `findvolume`, which is used to find an
+approximate volume at a given pressure, energy, or bulk modulus based on an
+equation of state. A result is not always guaranteed, especially when the
+equation of state is not a monotonic function of volume. However, according to
+experience, `P(V)` relation is usually a monotonic function. So we suggest using
+`PressureForm` to find the corresponding volume.
+
+## Usage
+
+```julia
+julia> using EquationsOfState, EquationsOfState.Collections, EquationsOfState.Find, Unitful, UnitfulAtomic
+
+julia> pressures = collect(0:20:200) .* u"GPa";
+
+julia> eos = BirchMurnaghan3rd(167u"angstrom^3", 2600u"kbar", 4.0);
+
+julia> volumes = map(
+           p -> findvolume(PressureForm(), eos, p, (eps() * u"bohr^3", eos.v0 * 1.3)),
+           pressures
+       )
+[ Info: Using method "Roots.Bisection"...
+[ Info: Using method "Roots.Bisection"...
+[ Info: Using method "Roots.Bisection"...
+[ Info: Using method "Roots.Bisection"...
+[ Info: Using method "Roots.Bisection"...
+[ Info: Using method "Roots.Bisection"...
+[ Info: Using method "Roots.Bisection"...
+[ Info: Using method "Roots.Bisection"...
+[ Info: Using method "Roots.Bisection"...
+[ Info: Using method "Roots.Bisection"...
+[ Info: Using method "Roots.Bisection"...
+11-element Array{Quantity{Float64,𝐋^3,Unitful.FreeUnits{(Å^3,),𝐋^3,nothing}},1}:
+              167.0 Å^3
+ 156.14036210727835 Å^3
+ 147.99803635986564 Å^3
+ 141.51093713795865 Å^3
+ 136.13864615965332 Å^3
+ 131.56784031939347 Å^3
+ 127.60046278645824 Å^3
+ 124.10332447387113 Å^3
+ 120.98257680606459 Å^3
+ 118.16962836248427 Å^3
+ 115.61284838696814 Å^3
+```
+
+Here we let the algorithm choose the bisection root-finding method to find the
+`volumes` corresponding to `pressures`.
+
+A figure is plotted below to verify our results, and it fits very well.
+
+![findvolume](assets/findvolume.png)
+
+## Public interfaces
+
+```@docs
+findvolume(form::EquationForm, eos::EquationOfState, y, x0, method)
+```
diff --git a/docs/src/Installation.md b/docs/src/Installation.md
new file mode 100644
index 0000000..ecb8edc
--- /dev/null
+++ b/docs/src/Installation.md
@@ -0,0 +1,53 @@
+# Installation
+
+To install this package, first, you need to install a `julia` executable from
+[its official website](https://julialang.org/downloads/). Version `v1.0.0` and
+above is required. This package may not work on `v0.7` and below.
+
+If you are using a Mac, and have [Homebrew](https://brew.sh) installed, open
+`Terminal.app` and type:
+
+```shell
+brew cask install julia
+```
+
+Now I am using [macOS](https://en.wikipedia.org/wiki/MacOS) as a standard
+platform to explain the following steps:
+
+1. Open `Terminal.app`, and type `julia` to start a Julia session.
+
+2. Run
+
+   ```julia
+   julia> using Pkg; Pkg.update()
+
+   julia> Pkg.add(PackageSpec(url="https://github.com/PainterQubits/Unitful.jl"))
+
+   julia> Pkg.add(PackageSpec(url="https://github.com/sostock/UnitfulAtomic.jl"))
+
+   julia> Pkg.add(PackageSpec(url="https://github.com/MineralsCloud/EquationsOfState.jl"))
+   ```
+
+   and wait for its finish.
+
+3. Run
+
+   ```julia
+   julia> using EquationsOfState, EquationsOfState.Collections, EquationsOfState.Find, EquationsOfState.NonlinearFitting, Unitful, UnitfulAtomic
+   ```
+
+   and have fun!
+
+4. While using, please keep this Julia session alive. Restarting may recompile
+   the package and cost some time.
+
+## Reinstall
+
+1. In the same Julia session, run
+
+   ```julia
+   julia> Pkg.rm("EquationsOfState"); Pkg.gc()
+   ```
+
+2. Press `ctrl+d` to quit the current session. Start a new Julia session and
+   repeat the above steps.
diff --git a/docs/src/NonlinearFitting.md b/docs/src/NonlinearFitting.md
index 1a76454..8694d18 100644
--- a/docs/src/NonlinearFitting.md
+++ b/docs/src/NonlinearFitting.md
@@ -6,9 +6,24 @@ CurrentModule = EquationsOfState.NonlinearFitting
 
 From Ref. 1,
 
-> The equations of state depend nonlinearly of a collection of parameters, $E_0$, $V_0$, $B_0$, $B_0'$, ..., that represent physical properties of the solid at equilibrium and can, in principle, be obtained expermentally by independent methods. The use of a given analytical EOS may have significant influence on the results obtained, particularly because the parameters are far from being independent. The number of parameters has to be considered in comparing the goodness of fit of functional forms with different analytical flexibility. The possibility of using too many parameters, beyond what is physically justified by the information contained in the experimental data, is a serious aspect that deserves consideration.
+> The equations of state depend nonlinearly of a collection of parameters,
+> $E_0$, $V_0$, $B_0$, $B_0'$, ..., that represent physical properties of the
+> solid at equilibrium and can, in principle, be obtained expermentally by
+> independent methods. The use of a given analytical EOS may have significant
+> influence on the results obtained, particularly because the parameters are far
+> from being independent. The number of parameters has to be considered in
+> comparing the goodness of fit of functional forms with different analytical
+> flexibility. The possibility of using too many parameters, beyond what is
+> physically justified by the information contained in the experimental data, is
+> a serious aspect that deserves consideration.
 
-In [`EquationsOfState`](https://github.com/MineralsCloud/EquationsOfState.jl), the nonlinear fitting is currently implemented by [`LsqFit`](https://github.com/JuliaNLSolvers/LsqFit.jl), a small library that provides basic least-squares fitting in pure Julia. It only utilizes the *Levenberg-Marquardt algorithm* for non-linear fitting. See its [documentation](https://github.com/JuliaNLSolvers/LsqFit.jl/blob/master/README.md) for more information.
+In [`EquationsOfState`](https://github.com/MineralsCloud/EquationsOfState.jl),
+the nonlinear fitting is currently implemented by
+[`LsqFit`](https://github.com/JuliaNLSolvers/LsqFit.jl), a small library that
+provides basic least-squares fitting in pure Julia. It only utilizes the
+_Levenberg-Marquardt algorithm_ for non-linear fitting. See its
+[documentation](https://github.com/JuliaNLSolvers/LsqFit.jl/blob/master/README.md)
+for more information.
 
 ## Usage
 
@@ -47,8 +62,8 @@ volumes = [
     54.3808371612,
     55.8775030703,
     57.4014349722,
-    58.9526328669
-]
+    58.9526328669,
+];
 energies = [
     -7.63622156576,
     -8.16831294894,
@@ -77,23 +92,30 @@ energies = [
     -10.1197772808,
     -9.99504030111,
     -9.86535084973,
-    -9.73155247952
-]
+    -9.73155247952,
+];
 
-lsqfit(EnergyForm(), BirchMurnaghan3rd(40, 0.5, 4, 0), volumes, energies)
-lsqfit(EnergyForm(), Murnaghan(41, 0.5, 4, 0), volumes, energies)
-lsqfit(EnergyForm(), PoirierTarantola3rd(41, 0.5, 4, 0), volumes, energies)
-lsqfit(EnergyForm(), Vinet(41, 0.5, 4, 0), volumes, energies)
+julia> lsqfit(EnergyForm(), BirchMurnaghan3rd(40, 0.5, 4, 0), volumes, energies)
+BirchMurnaghan3rd{Float64}(40.989265727925826, 0.5369258245608038, 4.1786442319302015, -10.842803908298968)
+
+julia> lsqfit(EnergyForm(), Murnaghan(41, 0.5, 4, 0), volumes, energies)
+Murnaghan{Float64}(41.13757924894751, 0.5144967655882123, 3.912386317519504, -10.836794511015869)
+
+julia> lsqfit(EnergyForm(), PoirierTarantola3rd(41, 0.5, 4, 0), volumes, energies)
+PoirierTarantola3rd{Float64}(40.86770643567383, 0.5667729960008705, 4.331688934942696, -10.851486685029547)
+
+julia> lsqfit(EnergyForm(), Vinet(41, 0.5, 4, 0), volumes, energies)
+Vinet{Float64}(40.91687567368755, 0.5493839427734198, 4.30519294991197, -10.846160810968053)
 ```
+
 Then 4 different equations of state will be fitted.
 
 ## Public interfaces
 
 ```@docs
-lsqfit(::EquationForm, eos::E, xdata::X, ydata::Y; debug = false, kwargs...) where {E<:EquationOfState,X<:AbstractVector,Y<:AbstractVector}
+lsqfit(form::EquationForm, eos::EquationOfState{<:Real}, xdata::AbstractVector{<:Real}, ydata::AbstractVector{<:Real}; debug::Bool, kwargs...)
 ```
 
-
 ## References
 
-1. [A. Otero-De-La-Roza, V. Luaña, *Computer Physics Communications*. **182**, 1708–1720 (2011), doi:10.1016/j.cpc.2011.04.016.](https://www.sciencedirect.com/science/article/pii/S0010465511001470)
+1. [A. Otero-De-La-Roza, V. Luaña, _Computer Physics Communications_. **182**, 1708–1720 (2011), doi:10.1016/j.cpc.2011.04.016.](https://www.sciencedirect.com/science/article/pii/S0010465511001470)
diff --git a/docs/src/assets/findvolume.png b/docs/src/assets/findvolume.png
new file mode 100644
index 0000000..050aea4
Binary files /dev/null and b/docs/src/assets/findvolume.png differ
diff --git a/docs/src/index.md b/docs/src/index.md
index a73596d..45950ed 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -2,25 +2,30 @@
 
 !!! note
 
-    Starting from `v2.0.0`, `EquationsOfState.jl` has been rewritten. 
+    Starting from `v2.0.0`, `EquationsOfState.jl` has been rewritten.
     The former behaviors of `EquationsOfState.jl` (`v1.1.4` and earlier) will be deprecated.
     Please follow the latest documentation.
 
 ## Package Features
 
-- Calculate energy, pressure, and bulk modulus of an `EquationOfState` of a specific volume.
-- Fit an `EquationOfState` on a series of volumes using least-squares fitting method.
+- Calculate energy, pressure, and bulk modulus of an `EquationOfState` on a (an)
+  volume (array of volumes).
+- Fit an `EquationOfState` on a series of volumes using least-squares fitting
+  method.
 - Fit an `EquationOfState` on a series of volumes linearly.
-- Find the corresponding volume on an `EquationOfState` given one of the energy, pressure, and bulk modulus.
+- Find the corresponding volume of an `EquationOfState` given an (a) energy,
+  pressure, and bulk modulus.
 
-See the [Index](@ref main-index) for the complete list of documented functions and types.
+See the [Index](@ref main-index) for the complete list of documented functions
+and types.
 
 ## Manual Outline
 
 ```@contents
 Pages = [
     "Collections.md",
-    "NonlinearFitting.md"
+    "NonlinearFitting.md",
+    "Find.md",
 ]
 Depth = 3
 ```
@@ -28,4 +33,5 @@ Depth = 3
 ### [Index](@id main-index)
 
 ```@index
+
 ```
diff --git a/src/Collections.jl b/src/Collections.jl
index 0dbd21e..af843c1 100644
--- a/src/Collections.jl
+++ b/src/Collections.jl
@@ -1,13 +1,7 @@
 """
-# module Collections
-
-
-
-# Examples
-
-```jldoctest
-julia>
-```
+This module provides `EquationOfState` types and `apply` methods to calculate
+energy, pressure, or bulk modulus of an `EquationOfState` on
+a (an) volume (array of volumes).
 """
 module Collections
 
@@ -36,27 +30,42 @@ export apply,
 """
     EquationOfState{T}
 
-An abstraction of equations of state, where `T` specifies the elements' type.
+An abstraction of equations of state, where `T` specifies the elements' common type.
 """
 abstract type EquationOfState{T} end
 
 """
     FiniteStrainEquationOfState{T} <: EquationOfState{T}
 
-An abstraction of finite strain equations of state.
+An abstraction of finite strain equations of state, where `T` specifies the elements' common type.
 """
 abstract type FiniteStrainEquationOfState{T} <: EquationOfState{T} end
 
 """
-    Murnaghan(v0, b0, bp0, e0=0)
+    Murnaghan(v0, b0, bp0, e0)
 
 Create a Murnaghan equation of state. The elements' type will be handled automatically.
 
+This equation of state can have units. The units are specified in [`Unitful.jl`](https://github.com/PainterQubits/Unitful.jl)'s
+`@u_str` style.
+
 # Arguments
 - `v0`: the volume of solid at zero pressure.
 - `b0`: the bulk modulus of solid at zero pressure.
 - `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`.
+- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.
+
+# Examples
+```jldoctest
+julia> Murnaghan(1, 2, 3.0)
+Murnaghan{Float64}(1.0, 2.0, 3.0, 0.0)
+
+julia> Murnaghan(Int8(1), 2//1, 3.0, 4)
+Murnaghan{Float64}(1.0, 2.0, 3.0, 4.0)
+
+julia> Murnaghan(1u"nm^3", 2u"GPa", 3, 3.0u"eV")
+Murnaghan{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0, 3.0 eV)
+```
 """
 struct Murnaghan{T} <: EquationOfState{T}
     v0::T
@@ -73,14 +82,28 @@ Murnaghan(v0::AbstractQuantity, b0::AbstractQuantity, bp0) =
     Murnaghan(v0, b0, bp0, 0 * upreferred(Unitful.J))
 
 """
-    BirchMurnaghan2nd(v0, b0, e0=0)
+    BirchMurnaghan2nd(v0, b0, e0)
 
 Create a Birch–Murnaghan 2nd order equation of state. The elements' type will be handled automatically.
 
 # Arguments
 - `v0`: the volume of solid at zero pressure.
 - `b0`: the bulk modulus of solid at zero pressure.
-- `e0=0`: the energy of solid at zero pressure. By default is `0`.
+- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.
+
+See also: [`BirchMurnaghan3rd`](@ref), [`BirchMurnaghan4th`](@ref)
+
+# Examples
+```jldoctest
+julia> BirchMurnaghan2nd(1, 2.0)
+BirchMurnaghan2nd{Float64}(1.0, 2.0, 0.0)
+
+julia> BirchMurnaghan2nd(Int8(1), 2//1, 0.0)
+BirchMurnaghan2nd{Float64}(1.0, 2.0, 0.0)
+
+julia> BirchMurnaghan2nd(1u"nm^3", 2u"GPa", 3.0u"eV")
+BirchMurnaghan2nd{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0 eV)
+```
 """
 struct BirchMurnaghan2nd{T} <: FiniteStrainEquationOfState{T}
     v0::T
@@ -96,7 +119,7 @@ BirchMurnaghan2nd(v0::AbstractQuantity, b0::AbstractQuantity) =
     BirchMurnaghan2nd(v0, b0, 0 * upreferred(Unitful.J))
 
 """
-    BirchMurnaghan3rd(v0, b0, bp0, e0=0)
+    BirchMurnaghan3rd(v0, b0, bp0, e0)
 
 Create a Birch–Murnaghan 3rd order equation of state. The elements' type will be handled automatically.
 
@@ -104,7 +127,21 @@ Create a Birch–Murnaghan 3rd order equation of state. The elements' type will
 - `v0`: the volume of solid at zero pressure.
 - `b0`: the bulk modulus of solid at zero pressure.
 - `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`.
+- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.
+
+See also: [`BirchMurnaghan2nd`](@ref), [`BirchMurnaghan4th`](@ref)
+
+# Examples
+```jldoctest
+julia> BirchMurnaghan3rd(1, 2.0, 3)
+BirchMurnaghan3rd{Float64}(1.0, 2.0, 3.0, 0.0)
+
+julia> BirchMurnaghan3rd(Int8(1), 2//1, 4, 0.0)
+BirchMurnaghan3rd{Float64}(1.0, 2.0, 4.0, 0.0)
+
+julia> BirchMurnaghan3rd(1u"nm^3", 2u"GPa", 4.0, 3u"eV")
+BirchMurnaghan3rd{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 4.0, 3.0 eV)
+```
 """
 struct BirchMurnaghan3rd{T} <: FiniteStrainEquationOfState{T}
     v0::T
@@ -121,7 +158,7 @@ BirchMurnaghan3rd(v0::AbstractQuantity, b0::AbstractQuantity, bp0) =
     BirchMurnaghan3rd(v0, b0, bp0, 0 * upreferred(Unitful.J))
 
 """
-    BirchMurnaghan4th(v0, b0, bp0, bpp0, e0=0)
+    BirchMurnaghan4th(v0, b0, bp0, bpp0, e0)
 
 Create a Birch–Murnaghan 4th order equation of state. The elements' type will be handled automatically.
 
@@ -130,7 +167,21 @@ Create a Birch–Murnaghan 4th order equation of state. The elements' type will
 - `b0`: the bulk modulus of solid at zero pressure.
 - `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
 - `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`.
+- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.
+
+See also: [`BirchMurnaghan2nd`](@ref), [`BirchMurnaghan4th`](@ref)
+
+# Examples
+```jldoctest
+julia> BirchMurnaghan4th(1, 2.0, 3, 4)
+BirchMurnaghan4th{Float64}(1.0, 2.0, 3.0, 4.0, 0.0)
+
+julia> BirchMurnaghan4th(Int8(1), 2//1, 4, 5.0, Float16(6))
+BirchMurnaghan4th{Float64}(1.0, 2.0, 4.0, 5.0, 6.0)
+
+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^3, 2.0 GPa, 3.0, 4.0 GPa^-1, 5.0 eV)
+```
 """
 struct BirchMurnaghan4th{T} <: FiniteStrainEquationOfState{T}
     v0::T
@@ -153,14 +204,28 @@ BirchMurnaghan4th(
 ) = BirchMurnaghan4th(v0, b0, bp0, bpp0, 0 * upreferred(Unitful.J))
 
 """
-    PoirierTarantola2nd(v0, b0, e0=0)
+    PoirierTarantola2nd(v0, b0, e0)
 
 Create a Poirier–Tarantola order equation of state. The elements' type will be handled automatically.
 
 # Arguments
 - `v0`: the volume of solid at zero pressure.
 - `b0`: the bulk modulus of solid at zero pressure.
-- `e0=0`: the energy of solid at zero pressure. By default is `0`.
+- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.
+
+See also: [`PoirierTarantola3rd`](@ref), [`PoirierTarantola4th`](@ref)
+
+# Examples
+```jldoctest
+julia> PoirierTarantola2nd(1, 2.0)
+PoirierTarantola2nd{Float64}(1.0, 2.0, 0.0)
+
+julia> PoirierTarantola2nd(Int8(1), 2//1, 3.0)
+PoirierTarantola2nd{Float64}(1.0, 2.0, 3.0)
+
+julia> PoirierTarantola2nd(1u"nm^3", 2u"GPa", 3.0u"eV")
+PoirierTarantola2nd{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0 eV)
+```
 """
 struct PoirierTarantola2nd{T} <: FiniteStrainEquationOfState{T}
     v0::T
@@ -176,7 +241,7 @@ PoirierTarantola2nd(v0::AbstractQuantity, b0::AbstractQuantity) =
     PoirierTarantola2nd(v0, b0, 0 * upreferred(Unitful.J))
 
 """
-    PoirierTarantola3rd(v0, b0, bp0, e0=0)
+    PoirierTarantola3rd(v0, b0, bp0, e0)
 
 Create a Poirier–Tarantola 3rd order equation of state. The elements' type will be handled automatically.
 
@@ -184,7 +249,21 @@ Create a Poirier–Tarantola 3rd order equation of state. The elements' type wil
 - `v0`: the volume of solid at zero pressure.
 - `b0`: the bulk modulus of solid at zero pressure.
 - `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`.
+- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.
+
+See also: [`PoirierTarantola2nd`](@ref), [`PoirierTarantola4th`](@ref)
+
+# Examples
+```jldoctest
+julia> PoirierTarantola3rd(1, 2.0, 3)
+PoirierTarantola3rd{Float64}(1.0, 2.0, 3.0, 0.0)
+
+julia> PoirierTarantola3rd(Int8(1), 2//1, 3.0, Float16(4))
+PoirierTarantola3rd{Float64}(1.0, 2.0, 3.0, 4.0)
+
+julia> PoirierTarantola3rd(1u"nm^3", 2u"GPa", 3, 4.0u"eV")
+PoirierTarantola3rd{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0, 4.0 eV)
+```
 """
 struct PoirierTarantola3rd{T} <: FiniteStrainEquationOfState{T}
     v0::T
@@ -201,7 +280,7 @@ PoirierTarantola3rd(v0::AbstractQuantity, b0::AbstractQuantity, bp0) =
     PoirierTarantola3rd(v0, b0, bp0, 0 * upreferred(Unitful.J))
 
 """
-    PoirierTarantola4th(v0, b0, bp0, bpp0, e0=0)
+    PoirierTarantola4th(v0, b0, bp0, bpp0, e0)
 
 Create a Birch–Murnaghan 4th order equation of state. The elements' type will be handled automatically.
 
@@ -210,7 +289,21 @@ Create a Birch–Murnaghan 4th order equation of state. The elements' type will
 - `b0`: the bulk modulus of solid at zero pressure.
 - `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
 - `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`.
+- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.
+
+See also: [`PoirierTarantola2nd`](@ref), [`PoirierTarantola3rd`](@ref)
+
+# Examples
+```jldoctest
+julia> PoirierTarantola4th(1, 2.0, 3, 4)
+PoirierTarantola4th{Float64}(1.0, 2.0, 3.0, 4.0, 0.0)
+
+julia> PoirierTarantola4th(Int8(1), 2//1, 3.0, Float16(4), 5)
+PoirierTarantola4th{Float64}(1.0, 2.0, 3.0, 4.0, 5.0)
+
+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^3, 2.0 GPa, 3.0, 4.0 GPa^-1, 5.0 eV)
+```
 """
 struct PoirierTarantola4th{T} <: FiniteStrainEquationOfState{T}
     v0::T
@@ -233,7 +326,7 @@ PoirierTarantola4th(
 ) = PoirierTarantola4th(v0, b0, bp0, bpp0, 0 * upreferred(Unitful.J))
 
 """
-    Vinet(v0, b0, bp0, e0=0)
+    Vinet(v0, b0, bp0, e0)
 
 Create a Vinet equation of state. The elements' type will be handled automatically.
 
@@ -241,7 +334,19 @@ Create a Vinet equation of state. The elements' type will be handled automatical
 - `v0`: the volume of solid at zero pressure.
 - `b0`: the bulk modulus of solid at zero pressure.
 - `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`.
+- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.
+
+# Examples
+```jldoctest
+julia> Vinet(1, 2.0, 3)
+Vinet{Float64}(1.0, 2.0, 3.0, 0.0)
+
+julia> Vinet(Int8(1), 2//1, 3.0, Float16(4))
+Vinet{Float64}(1.0, 2.0, 3.0, 4.0)
+
+julia> Vinet(1u"nm^3", 2u"GPa", 3, 4.0u"eV")
+Vinet{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0, 4.0 eV)
+```
 """
 struct Vinet{T} <: EquationOfState{T}
     v0::T
@@ -288,8 +393,10 @@ BreenanStacey(v0::Real, b0::Real, γ0::Real) = BreenanStacey(v0, b0, γ0, 0)
 # ============================================================================ #
 """
     apply(EnergyForm(), eos::EquationOfState)
+    apply(PressureForm(), eos::EquationOfState)
+    apply(BulkModulusForm(), eos::EquationOfState)
 
-Return a function that can take a volume as a parameter, suitable for batch-applying.
+Return a function that takes a volume as a variable, suitable for mapping onto an array.
 
 # Examples
 ```jldoctest
@@ -309,13 +416,44 @@ julia> map(f, 1:1:10)
  1.6017034530570884
  1.6679539823686644
  1.7203642945516917
+
+julia> g = apply(PressureForm(), Vinet(1, 2, 3));
+
+julia> map(g, 1:1:10)
+10-element Array{Float64,1}:
+  0.0
+ -0.45046308428750254
+ -0.3384840350043251
+ -0.24010297221667418
+ -0.17314062272722755
+ -0.12795492664586872
+ -0.09677154467733216
+ -0.07468060255179591
+ -0.05864401631176751
+ -0.04674768462396211
+
+julia> h = apply(BulkModulusForm(), BirchMurnaghan3rd(1, 2, 3));
+
+julia> map(h, 1:1:10)
+10-element Array{Float64,1}:
+ 2.0
+ 0.9216086833346415
+ 0.444903691617472
+ 0.2540009203153288
+ 0.16193296566524193
+ 0.11130584492987289
+ 0.08076305569984538
+ 0.06103515625
+ 0.047609811583958425
+ 0.03808959181078831
 ```
 """
 apply(form::EnergyForm, eos::EquationOfState) = v -> apply(form, eos, v)
 """
-    apply(EnergyForm(), eos::Murnaghan, v)
+    apply(EnergyForm(), eos::EquationOfState, v)
 
-Return the energy of a `Murnaghan` equation of state on volume `v`.
+Return the energy of an `EquationOfState` on volume `v`. If `eos` has units,
+`v` must also has.
 """
 function apply(::EnergyForm, eos::Murnaghan, v)
     v0, b0, bp0, e0 = fieldvalues(eos)
@@ -324,22 +462,12 @@ function apply(::EnergyForm, eos::Murnaghan, v)
     y = (v0 / v)^bp0
     return e0 + b0 / bp0 * v * (y / x + 1) - v0 * b0 / x
 end
-"""
-    apply(EnergyForm(), eos::BirchMurnaghan2nd, v)
-
-Return the energy of a `BirchMurnaghan2nd` equation of state on volume `v`.
-"""
 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)
-
-Return the energy of a `BirchMurnaghan3rd` equation of state on volume `v`.
-"""
 function apply(::EnergyForm, eos::BirchMurnaghan3rd, v)
     v0, b0, bp0, e0 = fieldvalues(eos)
 
@@ -347,34 +475,18 @@ function apply(::EnergyForm, eos::BirchMurnaghan3rd, 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)
-
-Return the energy of a `BirchMurnaghan4th` equation of state on volume `v`.
-"""
 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)
-
-Return the energy of a `PoirierTarantola2nd` equation of state on volume `v`.
-"""
 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)
-
-Return the energy of a `PoirierTarantola3rd` equation of state on volume `v`.
-"""
 function apply(::EnergyForm, eos::PoirierTarantola3rd, v)
     v0, b0, bp0, e0 = fieldvalues(eos)
 
@@ -382,11 +494,6 @@ function apply(::EnergyForm, eos::PoirierTarantola3rd, v)
     xi = -3 * log(x)
     return e0 + b0 / 6 * v0 * xi^2 * ((bp0 - 2) * xi + 3)
 end
-"""
-    apply(EnergyForm(), eos::PoirierTarantola4th, v)
-
-Return the energy of a `PoirierTarantola4th` equation of state on volume `v`.
-"""
 function apply(::EnergyForm, eos::PoirierTarantola4th, v)
     v0, b0, bp0, bpp0, e0 = fieldvalues(eos)
 
@@ -395,11 +502,6 @@ function apply(::EnergyForm, eos::PoirierTarantola4th, v)
     h = b0 * bpp0 + bp0^2
     return e0 + b0 / 24v0 * xi^2 * ((h + 3bp0 + 3) * xi^2 + 4 * (bp0 + 2) * xi + 12)
 end
-"""
-    apply(EnergyForm(), eos::Vinet, v)
-
-Return the energy of a `Vinet` equation of state on volume `v`.
-"""
 function apply(::EnergyForm, eos::Vinet, v)
     v0, b0, bp0, e0 = fieldvalues(eos)
 
@@ -407,11 +509,6 @@ function apply(::EnergyForm, eos::Vinet, v)
     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)
-
-Return the energy of a `AntonSchmidt` equation of state on volume `v`.
-"""
 function apply(::EnergyForm, eos::AntonSchmidt, v)
     v0, β, n, e∞ = fieldvalues(eos)
 
@@ -425,69 +522,30 @@ end
 # ============================================================================ #
 #                              Pressure evaluation                             #
 # ============================================================================ #
-"""
-    apply(PressureForm(), eos::EquationOfState)
-
-Return a function that can take a volume as a parameter, suitable for batch-applying.
-
-# Examples
-```jldoctest
-julia> using EquationsOfState, EquationsOfState.Collections
-
-julia> f = apply(PressureForm(), Vinet(1, 2, 3));
-
-julia> map(f, 1:1:10)
-10-element Array{Float64,1}:
-  0.0                
- -0.45046308428750254
- -0.3384840350043251 
- -0.24010297221667418
- -0.17314062272722755
- -0.12795492664586872
- -0.09677154467733216
- -0.07468060255179591
- -0.05864401631176751
- -0.04674768462396211
-```
-"""
 apply(::PressureForm, eos::EquationOfState) = v -> apply(PressureForm(), eos, v)
 """
-    apply(PressureForm(), eos::Murnaghan, v)
+    apply(PressureForm(), eos::EquationOfState, v)
 
-Return the pressure of a `Murnaghan` equation of state on volume `v`.
+Return the pressure of an `EquationOfState` on volume `v`. If `eos` has units,
+`v` must also has.
 """
 function apply(::PressureForm, eos::Murnaghan, v)
     v0, b0, bp0 = fieldvalues(eos)
 
     return b0 / bp0 * ((v0 / v)^bp0 - 1)
 end
-"""
-    apply(PressureForm(), eos::BirchMurnaghan2nd, v)
-
-Return the pressure of a `BirchMurnaghan2nd` equation of state on volume `v`.
-"""
 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)
-
-Return the pressure of a `BirchMurnaghan3rd` equation of state on volume `v`.
-"""
 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)
-
-Return the pressure of a `BirchMurnaghan4th` equation of state on volume `v`.
-"""
 function apply(::PressureForm, eos::BirchMurnaghan4th, v)
     v0, b0, bp0, bpp0 = fieldvalues(eos)
 
@@ -495,22 +553,12 @@ function apply(::PressureForm, eos::BirchMurnaghan4th, v)
     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)
-
-Return the pressure of a `PoirierTarantola2nd` equation of state on volume `v`.
-"""
 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)
-
-Return the pressure of a `PoirierTarantola3rd` equation of state on volume `v`.
-"""
 function apply(::PressureForm, eos::PoirierTarantola3rd, v)
     v0, b0, bp0 = fieldvalues(eos)
 
@@ -518,11 +566,6 @@ function apply(::PressureForm, eos::PoirierTarantola3rd, v)
     xi = log(x)
     return -b0 * xi / 2x * ((bp0 - 2) * xi - 2)
 end
-"""
-    apply(PressureForm(), eos::PoirierTarantola4th, v)
-
-Return the pressure of a `PoirierTarantola4th` equation of state on volume `v`.
-"""
 function apply(::PressureForm, eos::PoirierTarantola4th, v)
     v0, b0, bp0, bpp0 = fieldvalues(eos)
 
@@ -531,11 +574,6 @@ function apply(::PressureForm, eos::PoirierTarantola4th, v)
     h = b0 * bpp0 + bp0^2
     return -b0 * xi / 6 / x * ((h + 3bp0 + 3) * xi^2 + 3 * (bp0 + 6) * xi + 6)
 end
-"""
-    apply(PressureForm(), eos::Vinet, v)
-
-Return the pressure of a `Vinet` equation of state on volume `v`.
-"""
 function apply(::PressureForm, eos::Vinet, v)
     v0, b0, bp0 = fieldvalues(eos)
 
@@ -543,22 +581,12 @@ function apply(::PressureForm, eos::Vinet, v)
     xi = 3 / 2 * (bp0 - 1)
     return 3b0 / x^2 * (1 - x) * exp(xi * (1 - x))
 end
-"""
-    apply(PressureForm(), eos::AntonSchmidt, v)
-
-Return the pressure of a `AntonSchmidt` equation of state on volume `v`.
-"""
 function apply(::PressureForm, eos::AntonSchmidt, v)
     v0, β, n = fieldvalues(eos)
 
     x = v / v0
     return -β * x^n * log(x)
 end
-"""
-    apply(PressureForm(), eos::BreenanStacey, v)
-
-Return the pressure of a `BreenanStacey` equation of state on volume `v`.
-"""
 function apply(::PressureForm, eos::BreenanStacey, v)
     v0, b0, γ0 = fieldvalues(eos)
 
@@ -571,36 +599,12 @@ end
 # ============================================================================ #
 #                            Bulk modulus evaluation                           #
 # ============================================================================ #
-"""
-    apply(BulkModulusForm(), eos::EquationOfState)
-
-Return a function that can take a volume as a parameter, suitable for batch-applying.
-
-# Examples
-```jldoctest
-julia> using EquationsOfState, EquationsOfState.Collections
-
-julia> f = apply(BulkModulusForm(), BirchMurnaghan3rd(1, 2, 3));
-
-julia> map(f, 1:1:10)
-10-element Array{Float64,1}:
- 2.0                 
- 0.9216086833346415  
- 0.444903691617472   
- 0.2540009203153288  
- 0.16193296566524193 
- 0.11130584492987289 
- 0.08076305569984538 
- 0.06103515625       
- 0.047609811583958425
- 0.03808959181078831 
-```
-"""
 apply(::BulkModulusForm, eos::EquationOfState) = v -> apply(BulkModulusForm(), eos, v)
 """
-    apply(BulkModulusForm(), eos::BirchMurnaghan2nd, v)
+    apply(BulkModulusForm(), eos::EquationOfState, v)
 
-Return the bulk modulus of a `BirchMurnaghan2nd` equation of state on volume `v`.
+Return the bulk modulus of an `EquationOfState` on volume `v`. If `eos` has units,
+`v` must also has.
 """
 function apply(::BulkModulusForm, eos::BirchMurnaghan2nd, v)
     v0, b0 = fieldvalues(eos)
@@ -608,22 +612,12 @@ function apply(::BulkModulusForm, eos::BirchMurnaghan2nd, v)
     f = ((v0 / v)^(2 / 3) - 1) / 2
     return b0 * (7f + 1) * (2f + 1)^(5 / 2)
 end
-"""
-    apply(BulkModulusForm(), eos::BirchMurnaghan3rd, v)
-
-Return the bulk modulus of a `BirchMurnaghan3rd` equation of state on volume `v`.
-"""
 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)
-
-Return the bulk modulus of a `BirchMurnaghan4th` equation of state on volume `v`.
-"""
 function apply(::BulkModulusForm, eos::BirchMurnaghan4th, v)
     v0, b0, bp0, bpp0 = fieldvalues(eos)
 
@@ -632,22 +626,12 @@ function apply(::BulkModulusForm, eos::BirchMurnaghan4th, v)
     return b0 / 6 * (2f + 1)^(5 / 2) *
            ((99h - 693bp0 + 1573) * f^3 + (27h - 108bp0 + 105) * f^2 + 6f * (3bp0 - 5) + 6)
 end
-"""
-    apply(BulkModulusForm(), eos::PoirierTarantola2nd, v)
-
-Return the bulk modulus of a `PoirierTarantola2nd` equation of state on volume `v`.
-"""
 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)
-
-Return the bulk modulus of a `PoirierTarantola3rd` equation of state on volume `v`.
-"""
 function apply(::BulkModulusForm, eos::PoirierTarantola3rd, v)
     v0, b0, bp0 = fieldvalues(eos)
 
@@ -655,11 +639,6 @@ function apply(::BulkModulusForm, eos::PoirierTarantola3rd, v)
     xi = log(x)
     return -b0 / 2x * (((bp0 - 2) * xi + 2 - 2bp0) * xi + 2)
 end
-"""
-    apply(BulkModulusForm(), eos::PoirierTarantola4th, v)
-
-Return the bulk modulus of a `PoirierTarantola4th` equation of state on volume `v`.
-"""
 function apply(::BulkModulusForm, eos::PoirierTarantola4th, v)
     v0, b0, bp0, bpp0 = fieldvalues(eos)
 
@@ -669,11 +648,6 @@ function apply(::BulkModulusForm, eos::PoirierTarantola4th, v)
     return -b0 / (6x) *
            ((h + 3bp0 + 3) * xi^3 - 3 * xi^2 * (h + 2bp0 + 1) - 6xi * (bp0 + 1) - 6)
 end
-"""
-    apply(BulkModulusForm(), eos::Vinet, v)
-
-Return the bulk modulus of a `Vinet` equation of state on volume `v`.
-"""
 function apply(::BulkModulusForm, eos::Vinet, v)
     v0, b0, bp0 = fieldvalues(eos)
 
@@ -681,11 +655,6 @@ function apply(::BulkModulusForm, eos::Vinet, v)
     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)
-
-Return the bulk modulus of a `AntonSchmidt` equation of state on volume `v`.
-"""
 function apply(::BulkModulusForm, eos::AntonSchmidt, v)
     v0, β, n = fieldvalues(eos)
 
diff --git a/src/Find.jl b/src/Find.jl
index b5f8579..c1a7ec2 100644
--- a/src/Find.jl
+++ b/src/Find.jl
@@ -1,13 +1,6 @@
 """
-# module Find
-
-
-
-# Examples
-
-```jldoctest
-julia>
-```
+This module provides `findvolume` methods to find the volume at a given
+pressure, energy, or bulk modulus with(out) units.
 """
 module Find
 
@@ -31,11 +24,34 @@ using ..Collections: EquationOfState, apply
 
 export findvolume
 
+"""
+    findvolume(form, eos, y, x0, method)
+    findvolume(form, eos, y, x0::Union{AbstractVector,Tuple})
+
+Find a volume which leads to the given pressure, energy, or bulk modulus based on an `eos`.
+
+# Arguments
+- `form::EquationForm`: an `EquationForm` instance.
+- `eos::EquationOfState`: an equation of state. If it has units, `y` and `x0` must also have.
+- `y`: a pressure, energy, or bulk modulus.
+- `x0`: can be either a range of volumes (`Vector`, `Tuple`, etc.) or just a single volume.
+    Units can be provided if necessary.
+- `method::Roots.AbstractUnivariateZeroMethod`: a method used to find the root of an equation.
+    If it is omitted, the algorithm will traverse all possible methods of 
+    [Roots.jl](https://github.com/JuliaMath/Roots.jl). And the `x0` parameter must be
+    an array or a tuple, of which only the maximum and minimum values will be used in the
+    root-finding process.
+"""
 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, x0::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
diff --git a/src/FiniteStrains.jl b/src/FiniteStrains.jl
index 5b68efc..b17a4d0 100644
--- a/src/FiniteStrains.jl
+++ b/src/FiniteStrains.jl
@@ -1,13 +1,6 @@
 """
-# module FiniteStrains
-
-
-
-# Examples
-
-```jldoctest
-julia>
-```
+This module provides some `FiniteStrain` types and `getstrain` methods to
+calculate various strains.
 """
 module FiniteStrains
 
diff --git a/src/LinearFitting.jl b/src/LinearFitting.jl
index 491f757..3e3cc91 100644
--- a/src/LinearFitting.jl
+++ b/src/LinearFitting.jl
@@ -1,13 +1,5 @@
 """
-# module LinearFitting
-
-
-
-# Examples
-
-```jldoctest
-julia>
-```
+This module provides some linear fitting methods.
 """
 module LinearFitting
 
@@ -54,8 +46,7 @@ function energy_volume_expansion(
     highest_order::Int = degree(p),
 )
     # The zeroth order value plus values from the first to the ``highest_order`.
-    p(v) +
-    dot(
+    p(v) + dot(
         energy_volume_derivatives(s, v0, v, p, highest_order),
         getstrain(s, v0, v) .^ collect(1:highest_order),
     )
diff --git a/src/NonlinearFitting.jl b/src/NonlinearFitting.jl
index 2acbb24..3106c4f 100644
--- a/src/NonlinearFitting.jl
+++ b/src/NonlinearFitting.jl
@@ -1,13 +1,6 @@
 """
-# module NonlinearFitting
-
-
-
-# Examples
-
-```jldoctest
-julia>
-```
+This module provides `lsqfit` methods for fitting an equation of state
+with(out) units.
 """
 module NonlinearFitting
 
@@ -27,11 +20,12 @@ Fit an equation of state using least-squares fitting method (with the Levenberg-
 
 # Arguments
 - `form::EquationForm`: an `EquationForm` instance. If `EnergyForm`, fit ``E(V)``; if `PressureForm`, fit ``P(V)``; if `BulkModulusForm`, fit ``B(V)``.
-- `eos::EquationOfState`: a trial equation of state.
-- `xdata::AbstractVector`: a vector of volumes.
-- `ydata::AbstractVector`: a vector of energies, pressures, or bulk moduli.
+- `eos::EquationOfState`: a trial equation of state. If it has units, `xdata` and `ydata` must also have.
+- `xdata::AbstractVector`: a vector of volumes (``V``), with(out) units.
+- `ydata::AbstractVector`: a vector of energies (``E``), pressures (``P``), or bulk moduli (``B``), with(out) units. It must be consistent with `form`.
 - `debug::Bool=false`: if `true`, then an `LsqFit.LsqFitResult` is returned, containing estimated Jacobian, residuals, etc.; if `false`, a fitted `EquationOfState` is returned. The default value is `false`.
-- `kwargs`: the rest keyword arguments that will be sent to `LsqFit.curve_fit`. See its [documentation](https://github.com/JuliaNLSolvers/LsqFit.jl/blob/master/README.md).
+- `kwargs`: the rest keyword arguments are the same as that of `LsqFit.curve_fit`. See its [documentation](https://github.com/JuliaNLSolvers/LsqFit.jl/blob/master/README.md)
+    and [tutorial](https://julianlsolvers.github.io/LsqFit.jl/latest/tutorial/).
 """
 function lsqfit(
     form::EquationForm,
@@ -63,13 +57,15 @@ function lsqfit(
     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)]
+    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)]...
-        )
+        return E([data[i] * upreferred(u) |> u for (i, u) in enumerate(original_units)]...)
     end
     return result
 end  # function lsqfit
diff --git a/test/NonlinearFitting.jl b/test/NonlinearFitting.jl
index 3d86801..2eab177 100644
--- a/test/NonlinearFitting.jl
+++ b/test/NonlinearFitting.jl
@@ -1,5 +1,7 @@
 using Test
 
+using Unitful, UnitfulAtomic
+
 using EquationsOfState
 using EquationsOfState.Collections
 using EquationsOfState.NonlinearFitting
@@ -103,7 +105,12 @@ end
 end
 
 @testset "Test fitting bulk modulus with different element types" begin
-    result = BirchMurnaghan3rd(7.218928431312577, 5.007900469653902, 4.06037725509478, 0.0) |> Collections.fieldvalues
+    result = BirchMurnaghan3rd(
+        7.218928431312577,
+        5.007900469653902,
+        4.06037725509478,
+        0.0,
+    ) |> Collections.fieldvalues
     @test isapprox(
         lsqfit(
             BulkModulusForm(),
@@ -597,7 +604,12 @@ end
         BirchMurnaghan3rd(224, 0.0006, 4, -323),
         volumes,
         energies,
-    ) |> Collections.fieldvalues ≈ BirchMurnaghan3rd(224.444565, 0.00062506191050572675, 3.740369, -323.417714) |> Collections.fieldvalues
+    ) |> Collections.fieldvalues ≈ BirchMurnaghan3rd(
+        224.444565,
+        0.00062506191050572675,
+        3.740369,
+        -323.417714,
+    ) |> Collections.fieldvalues
     @test isapprox(
         lsqfit(
             EnergyForm(),
@@ -631,3 +643,92 @@ end
     )
     # @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)
 end
+
+@testset "`Test w2k-lda-na.dat` from `Gibbs2` with units" begin
+    data = [
+        159.9086 -323.4078898
+        162.5738 -323.4089153
+        165.2389 -323.4098546
+        167.9041 -323.410722
+        170.5692 -323.4115195
+        173.2344 -323.4122481
+        175.8995 -323.4129189
+        178.5647 -323.413528
+        181.2298 -323.4140871
+        183.8949 -323.4145889
+        186.5601 -323.4150471
+        189.2252 -323.415459
+        191.8904 -323.4158302
+        194.5555 -323.4161579
+        197.2207 -323.4164498
+        199.8858 -323.4167071
+        202.551 -323.4169305
+        205.2161 -323.4171194
+        207.8812 -323.4172809
+        210.5464 -323.4174144
+        213.2115 -323.4175216
+        215.8767 -323.4176029
+        218.5418 -323.417661
+        221.207 -323.4176975
+        223.8721 -323.41771
+        226.5373 -323.4177051
+        229.2024 -323.417682
+        231.8675 -323.4176375
+        234.5327 -323.417579
+        237.1978 -323.4175048
+        239.863 -323.4174142
+        242.5281 -323.4173101
+        245.1933 -323.4171922
+        247.8584 -323.4170611
+        250.5236 -323.4169184
+        253.1887 -323.4167647
+        255.8538 -323.4166002
+        258.519 -323.4164244
+        261.1841 -323.4162386
+        263.8493 -323.4160446
+        266.5144 -323.4158421
+        269.1796 -323.4156312
+        271.8447 -323.4154125
+        274.5098 -323.4151861
+        277.175 -323.4149528
+        279.8401 -323.4147131
+        282.5053 -323.414467
+        285.1704 -323.414215
+        287.8356 -323.4139583
+        290.5007 -323.4136953
+        293.1659 -323.4134285
+        295.831 -323.4131559
+        298.4961 -323.4128797
+        301.1613 -323.4125984
+        303.8264 -323.4123147
+        306.4916 -323.4120269
+        309.1567 -323.411736
+        311.8219 -323.4114399
+        314.487 -323.4111421
+        317.1522 -323.4108418
+        319.8173 -323.4105393
+    ]
+    volumes = data[:, 1] .* u"bohr^3"
+    energies = data[:, 2] .* u"Ry"
+    fitted_eos = lsqfit(
+        EnergyForm(),
+        BirchMurnaghan3rd(224 * u"bohr^3", 0.0006 * u"Ry/bohr^3", 4, -323 * u"Ry"),
+        volumes,
+        energies,
+    )
+    @test ustrip.(fitted_eos |> Collections.fieldvalues) ≈
+    ustrip.(
+        BirchMurnaghan3rd(
+            224.444565 * u"bohr^3",
+            0.00062506191050572675 * u"Ry/bohr^3",
+            3.740369,
+            -323.417714 * u"Ry",
+        ) |> Collections.fieldvalues
+    )
+    @test ustrip.(
+        lsqfit(EnergyForm(), BirchMurnaghan3rd(224u"bohr^3", 10u"GPa", 3.75, -161u"hartree"), volumes, energies) |> Collections.fieldvalues
+    ) ≈
+    ustrip.(
+        BirchMurnaghan3rd(224.4445656763778u"bohr^3", 9.194980249913018u"GPa", 3.7403684211716297, -161.70885710742223u"hartree") |> Collections.fieldvalues
+    )
+end