diff --git a/docs/src/features.md b/docs/src/features.md
index a779cd57c..340bb056e 100644
--- a/docs/src/features.md
+++ b/docs/src/features.md
@@ -39,12 +39,11 @@
 
 * Third-party integrations:
     - Seamless integration with many standard [Input and output formats](@ref).
-    - Integration with [ASE](https://wiki.fysik.dtu.dk/ase/) and
-      [AtomsBase](https://github.com/JuliaMolSim/AtomsBase.jl) for passing
-      atomic structures (see [AtomsBase integration](@ref)).
+    - Integration with [AtomsBase](https://github.com/JuliaMolSim/AtomsBase.jl) for passing
+      atomic structures (see [AtomsBase integration](@ref)) as well as
+      the [Atomistic simulation environment](@ref).
     - [Wannierization using Wannier.jl or Wannier90](@ref)
 
-
 * Runs out of the box on Linux, macOS and Windows
 
 Missing a feature? Look for an open issue or [create a new one](https://github.com/JuliaMolSim/DFTK.jl/issues).
diff --git a/examples/geometry_optimization.jl b/examples/geometry_optimization.jl
index bc6db5775..ee0093696 100644
--- a/examples/geometry_optimization.jl
+++ b/examples/geometry_optimization.jl
@@ -3,32 +3,35 @@
 # We use DFTK and the [GeometryOptimization](https://github.com/JuliaMolSim/GeometryOptimization.jl/)
 # package to find the minimal-energy bond length of the ``H_2`` molecule.
 # First we set up an appropriate `DFTKCalculator` (see [AtomsCalculators integration](@ref)),
-# for which we use the LDA model just like in the [Tutorial](@ref) for silicon.
+# for which we use the LDA model just like in the [Tutorial](@ref) for silicon
+# in combination with a pseudodojo pseudopotential (see [Pseudopotentials](@ref)).
 
 using DFTK
+using PseudoPotentialData
 
+pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.oncvpsp3.standard.upf")
 calc = DFTKCalculator(;
-    model_kwargs = (; functionals=LDA()),        # model_DFT keyword arguments
+    model_kwargs = (; functionals=LDA(), pseudopotentials),  # model_DFT keyword arguments
     basis_kwargs = (; kgrid=[1, 1, 1], Ecut=10)  # PlaneWaveBasis keyword arguments
 )
 
 # Next we set up an initial hydrogen molecule within a box of vacuum.
 # We use the parameters of the
-# [equivalent tutorial from ABINIT](https://docs.abinit.org/tutorial/base1/),
-# that is a simulation box of 10 bohr times 10 bohr times 10 bohr and a
+# [equivalent tutorial from ABINIT](https://docs.abinit.org/tutorial/base1/)
+# and DFTK's [AtomsBase integration](@ref) to setup the hydrogen molecule.
+# We employ a simulation box of 10 bohr times 10 bohr times 10 bohr and a
 # pseudodojo pseudopotential.
 using LinearAlgebra
-using PseudoPotentialData
+using Unitful
+using UnitfulAtomic
 
 r0 = 1.4   # Initial bond length in Bohr
 a  = 10.0  # Box size in Bohr
 
-lattice = a * I(3)
-H = ElementPsp(:H, PseudoFamily("dojo.nc.sr.pbe.v0_4_1.oncvpsp3.standard.upf"))
-atoms = [H, H]
-positions = [zeros(3), lattice \ [r0, 0., 0.]]
-
-h2_crude = periodic_system(lattice, atoms, positions)
+cell_vectors = [[a, 0, 0]u"bohr", [0, a, 0]u"bohr", [0, 0, a]u"bohr"]
+h2_crude = periodic_system([:H => [0, 0, 0.]u"bohr",
+                            :H => [r0, 0, 0]u"bohr"],
+                           cell_vectors)
 
 # Finally we call `minimize_energy!` to start the geometry optimisation.
 # We use `verbosity=2` to get some insight into the minimisation.
diff --git a/src/standard_models.jl b/src/standard_models.jl
index cdb9a4f7d..a01d843d2 100644
--- a/src/standard_models.jl
+++ b/src/standard_models.jl
@@ -32,8 +32,12 @@ julia> model_atomic(system; pseudopotentials=Dict(:Si => "hgh/lda/si-q4"))
 ```
 same thing, but specify the pseudopotentials explicitly in a dictionary.
 """
-function model_atomic(system::AbstractSystem;
-                      pseudopotentials=fill(nothing, length(system)), kwargs...)
+function model_atomic(system::AbstractSystem; pseudopotentials, kwargs...)
+    # Note: We are enforcing to specify pseudopotentials at this interface
+    # (unlike the lower-level Model interface) because the argument is that
+    # automatically defaulting to the Coulomb potential will generally trip
+    # people over and could too easily lead to garbage results
+    #
     parsed = parse_system(system, pseudopotentials)
     model_atomic(parsed.lattice, parsed.atoms, parsed.positions;
                  parsed.magnetic_moments, kwargs...)
@@ -72,22 +76,31 @@ All other keyword arguments
 but `functional` are passed to [`model_atomic`](@ref) and from
 there to [`Model`](@ref).
 
+Note in particular that the `pseudopotential` keyword
+argument is mandatory to specify pseudopotential information. This can be easily
+achieved for example using the `PseudoFamily` struct from the `PseudoPotentialData`
+package as shown below:
+
 # Examples
 ```julia-repl
-julia> model_DFT(system; functionals=LDA(), temperature=0.01)
+julia> model_DFT(system; functionals=LDA(), temperature=0.01,
+                 pseudopotentials=PseudoFamily("dojo.nc.sr.lda.v0_4_1.oncvpsp3.standard.upf"))
+
 ```
 builds an [`LDA`](@ref) model for a passed system
 with specified smearing temperature.
 
 ```julia-repl
-julia> model_DFT(system; functionals=[:lda_x, :lda_c_pw], temperature=0.01)
+julia> model_DFT(system; functionals=[:lda_x, :lda_c_pw], temperature=0.01,
+                 pseudopotentials=PseudoFamily("dojo.nc.sr.lda.v0_4_1.oncvpsp3.standard.upf"))
 ```
 Alternative syntax specifying the functionals directly
 via their libxc codes.
 """
-function model_DFT(system::AbstractSystem;
-                   pseudopotentials=fill(nothing, length(system)),
-                   functionals, kwargs...)
+function model_DFT(system::AbstractSystem; pseudopotentials, functionals, kwargs...)
+    # Note: We are deliberately enforcing the user to specify pseudopotentials here.
+    # See the implementation of model_atomic for a rationale why
+    #
     # TODO Could check consistency between pseudos and passed functionals
     parsed = parse_system(system, pseudopotentials)
     _model_DFT(functionals, parsed.lattice, parsed.atoms, parsed.positions;