diff --git a/docs/gen_api_sources.sh b/docs/gen_api_sources.sh
index 67fd3df..e65108b 100644
--- a/docs/gen_api_sources.sh
+++ b/docs/gen_api_sources.sh
@@ -1,2 +1,2 @@
-# CD must be the reports directory!
-sphinx-apidoc -o source/generated ../ramannoodle
+# CD must be the docs directory
+sphinx-apidoc --module-first --remove-old -f --templatedir=source/_templates -f  --remove-old -o source/generated ../ramannoodle
diff --git a/docs/gen_test_coverage_badges.sh b/docs/gen_test_coverage_badges.sh
index 8dd6e72..a9af6d3 100644
--- a/docs/gen_test_coverage_badges.sh
+++ b/docs/gen_test_coverage_badges.sh
@@ -1,4 +1,4 @@
-# CD must be the reports directory!
+# CD must be the docs directory!
 cd ..
 coverage run -m pytest --junitxml=reports/junit/junit.xml test
 coverage xml -o  reports/coverage/coverage.xml
diff --git a/docs/source/_templates/module.rst.jinja b/docs/source/_templates/module.rst.jinja
new file mode 100644
index 0000000..49cca84
--- /dev/null
+++ b/docs/source/_templates/module.rst.jinja
@@ -0,0 +1,8 @@
+{%- if show_headings %}
+{{- [basename, "module"] | join(' ') | e | heading }}
+
+{% endif -%}
+.. automodule:: {{ qualname }}
+{%- for option in automodule_options %}
+   :{{ option }}:
+{%- endfor %}
diff --git a/docs/source/_templates/package.rst.jinja b/docs/source/_templates/package.rst.jinja
new file mode 100644
index 0000000..63c87c6
--- /dev/null
+++ b/docs/source/_templates/package.rst.jinja
@@ -0,0 +1,57 @@
+{%- macro automodule(modname, options) -%}
+.. automodule:: {{ modname }}
+{%- for option in options %}
+   :{{ option }}:
+{%- endfor %}
+{%- endmacro %}
+
+{%- macro toctree(docnames) -%}
+.. toctree::
+   :maxdepth: {{ maxdepth }}
+{% for docname in docnames %}
+   {{ docname }}
+{%- endfor %}
+{%- endmacro %}
+
+{%- if is_namespace %}
+{{- [pkgname, ""] | join(" ") | e | heading }}
+{% else %}
+{{- [pkgname, ""] | join(" ") | e | heading }}
+{% endif %}
+
+{%- if is_namespace %}
+.. py:module:: {{ pkgname }}
+{% endif %}
+
+{%- if modulefirst and not is_namespace %}
+{{ automodule(pkgname, automodule_options) }}
+{% endif %}
+
+{%- if subpackages %}
+Subpackages
+-----------
+
+{{ toctree(subpackages) }}
+{% endif %}
+
+{%- if submodules %}
+Submodules
+----------
+{% if separatemodules %}
+{{ toctree(submodules) }}
+{% else %}
+{%- for submodule in submodules %}
+{% if show_headings %}
+{{- [submodule, "module"] | join(" ") | e | heading(2) }}
+{% endif %}
+{{ automodule(submodule, automodule_options) }}
+{% endfor %}
+{%- endif %}
+{%- endif %}
+
+{%- if not modulefirst and not is_namespace %}
+Module contents
+---------------
+
+{{ automodule(pkgname, automodule_options) }}
+{% endif %}
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 22ebecc..bdbee6f 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -38,6 +38,7 @@
     'numpy': ('http://docs.scipy.org/doc/numpy', None),
     'scipy': ('http://docs.scipy.org/doc/scipy/reference', None),
     'torch': ('https://pytorch.org/docs/stable/', None),
+    'pymatgen': ('https://pymatgen.org', None),
 }
 
 intersphinx_disabled_domains = ['std']
diff --git a/docs/source/generated/ramannoodle.dataset.rst b/docs/source/generated/ramannoodle.dataset.rst
index 85478dd..da551bc 100644
--- a/docs/source/generated/ramannoodle.dataset.rst
+++ b/docs/source/generated/ramannoodle.dataset.rst
@@ -1,5 +1,10 @@
-ramannoodle.dataset package
-===========================
+ramannoodle.dataset
+====================
+
+.. automodule:: ramannoodle.dataset
+   :members:
+   :undoc-members:
+   :show-inheritance:
 
 Subpackages
 -----------
@@ -8,11 +13,3 @@ Subpackages
    :maxdepth: 4
 
    ramannoodle.dataset.torch
-
-Module contents
----------------
-
-.. automodule:: ramannoodle.dataset
-   :members:
-   :undoc-members:
-   :show-inheritance:
diff --git a/docs/source/generated/ramannoodle.dataset.torch.rst b/docs/source/generated/ramannoodle.dataset.torch.rst
index b3c3d93..ee2c6b7 100644
--- a/docs/source/generated/ramannoodle.dataset.torch.rst
+++ b/docs/source/generated/ramannoodle.dataset.torch.rst
@@ -1,17 +1,14 @@
-ramannoodle.dataset.torch package
-=================================
+ramannoodle.dataset.torch
+==========================
 
-Submodules
-----------
-
-ramannoodle.dataset.torch.dataset module
-----------------------------------------
-
-.. automodule:: ramannoodle.dataset.torch.dataset
+.. automodule:: ramannoodle.dataset.torch
    :members:
    :undoc-members:
    :show-inheritance:
 
+Submodules
+----------
+
 ramannoodle.dataset.torch.utils module
 --------------------------------------
 
@@ -19,11 +16,3 @@ ramannoodle.dataset.torch.utils module
    :members:
    :undoc-members:
    :show-inheritance:
-
-Module contents
----------------
-
-.. automodule:: ramannoodle.dataset.torch
-   :members:
-   :undoc-members:
-   :show-inheritance:
diff --git a/docs/source/generated/ramannoodle.dynamics.rst b/docs/source/generated/ramannoodle.dynamics.rst
index d1e2699..eef260e 100644
--- a/docs/source/generated/ramannoodle.dynamics.rst
+++ b/docs/source/generated/ramannoodle.dynamics.rst
@@ -1,27 +1,5 @@
-ramannoodle.dynamics package
-============================
-
-Submodules
-----------
-
-ramannoodle.dynamics.phonon module
-----------------------------------
-
-.. automodule:: ramannoodle.dynamics.phonon
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-ramannoodle.dynamics.trajectory module
---------------------------------------
-
-.. automodule:: ramannoodle.dynamics.trajectory
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-Module contents
----------------
+ramannoodle.dynamics
+=====================
 
 .. automodule:: ramannoodle.dynamics
    :members:
diff --git a/docs/source/generated/ramannoodle.io.rst b/docs/source/generated/ramannoodle.io.rst
index cbbc98a..043453f 100644
--- a/docs/source/generated/ramannoodle.io.rst
+++ b/docs/source/generated/ramannoodle.io.rst
@@ -1,5 +1,10 @@
-ramannoodle.io package
-======================
+ramannoodle.io
+===============
+
+.. automodule:: ramannoodle.io
+   :members:
+   :undoc-members:
+   :show-inheritance:
 
 Subpackages
 -----------
@@ -20,18 +25,10 @@ ramannoodle.io.generic module
    :undoc-members:
    :show-inheritance:
 
-ramannoodle.io.utils module
----------------------------
+ramannoodle.io.pymatgen module
+------------------------------
 
-.. automodule:: ramannoodle.io.utils
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-Module contents
----------------
-
-.. automodule:: ramannoodle.io
+.. automodule:: ramannoodle.io.pymatgen
    :members:
    :undoc-members:
    :show-inheritance:
diff --git a/docs/source/generated/ramannoodle.io.vasp.rst b/docs/source/generated/ramannoodle.io.vasp.rst
index 4cf5239..99617d1 100644
--- a/docs/source/generated/ramannoodle.io.vasp.rst
+++ b/docs/source/generated/ramannoodle.io.vasp.rst
@@ -1,5 +1,10 @@
-ramannoodle.io.vasp package
-===========================
+ramannoodle.io.vasp
+====================
+
+.. automodule:: ramannoodle.io.vasp
+   :members:
+   :undoc-members:
+   :show-inheritance:
 
 Submodules
 ----------
@@ -35,11 +40,3 @@ ramannoodle.io.vasp.xdatcar module
    :members:
    :undoc-members:
    :show-inheritance:
-
-Module contents
----------------
-
-.. automodule:: ramannoodle.io.vasp
-   :members:
-   :undoc-members:
-   :show-inheritance:
diff --git a/docs/source/generated/ramannoodle.pmodel.rst b/docs/source/generated/ramannoodle.pmodel.rst
index 5aa303e..0b6f9e9 100644
--- a/docs/source/generated/ramannoodle.pmodel.rst
+++ b/docs/source/generated/ramannoodle.pmodel.rst
@@ -1,5 +1,10 @@
-ramannoodle.pmodel package
-==========================
+ramannoodle.pmodel
+===================
+
+.. automodule:: ramannoodle.pmodel
+   :members:
+   :undoc-members:
+   :show-inheritance:
 
 Subpackages
 -----------
@@ -8,30 +13,3 @@ Subpackages
    :maxdepth: 4
 
    ramannoodle.pmodel.torch
-
-Submodules
-----------
-
-ramannoodle.pmodel.art module
------------------------------
-
-.. automodule:: ramannoodle.pmodel.art
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-ramannoodle.pmodel.interpolation module
----------------------------------------
-
-.. automodule:: ramannoodle.pmodel.interpolation
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-Module contents
----------------
-
-.. automodule:: ramannoodle.pmodel
-   :members:
-   :undoc-members:
-   :show-inheritance:
diff --git a/docs/source/generated/ramannoodle.pmodel.torch.rst b/docs/source/generated/ramannoodle.pmodel.torch.rst
index 2d69113..e7cc105 100644
--- a/docs/source/generated/ramannoodle.pmodel.torch.rst
+++ b/docs/source/generated/ramannoodle.pmodel.torch.rst
@@ -1,35 +1,5 @@
-ramannoodle.pmodel.torch package
-================================
-
-Submodules
-----------
-
-ramannoodle.pmodel.torch.gnn module
------------------------------------
-
-.. automodule:: ramannoodle.pmodel.torch.gnn
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-ramannoodle.pmodel.torch.train module
--------------------------------------
-
-.. automodule:: ramannoodle.pmodel.torch.train
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-ramannoodle.pmodel.torch.utils module
--------------------------------------
-
-.. automodule:: ramannoodle.pmodel.torch.utils
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-Module contents
----------------
+ramannoodle.pmodel.torch
+=========================
 
 .. automodule:: ramannoodle.pmodel.torch
    :members:
diff --git a/docs/source/generated/ramannoodle.rst b/docs/source/generated/ramannoodle.rst
index 26474fa..d4cace4 100644
--- a/docs/source/generated/ramannoodle.rst
+++ b/docs/source/generated/ramannoodle.rst
@@ -1,5 +1,10 @@
-ramannoodle package
-===================
+ramannoodle
+============
+
+.. automodule:: ramannoodle
+   :members:
+   :undoc-members:
+   :show-inheritance:
 
 Subpackages
 -----------
@@ -40,11 +45,3 @@ ramannoodle.exceptions module
    :members:
    :undoc-members:
    :show-inheritance:
-
-Module contents
----------------
-
-.. automodule:: ramannoodle
-   :members:
-   :undoc-members:
-   :show-inheritance:
diff --git a/docs/source/generated/ramannoodle.spectrum.rst b/docs/source/generated/ramannoodle.spectrum.rst
index 8861e28..09607b0 100644
--- a/docs/source/generated/ramannoodle.spectrum.rst
+++ b/docs/source/generated/ramannoodle.spectrum.rst
@@ -1,17 +1,14 @@
-ramannoodle.spectrum package
-============================
+ramannoodle.spectrum
+=====================
 
-Submodules
-----------
-
-ramannoodle.spectrum.raman module
----------------------------------
-
-.. automodule:: ramannoodle.spectrum.raman
+.. automodule:: ramannoodle.spectrum
    :members:
    :undoc-members:
    :show-inheritance:
 
+Submodules
+----------
+
 ramannoodle.spectrum.utils module
 ---------------------------------
 
@@ -19,11 +16,3 @@ ramannoodle.spectrum.utils module
    :members:
    :undoc-members:
    :show-inheritance:
-
-Module contents
----------------
-
-.. automodule:: ramannoodle.spectrum
-   :members:
-   :undoc-members:
-   :show-inheritance:
diff --git a/docs/source/generated/ramannoodle.structure.rst b/docs/source/generated/ramannoodle.structure.rst
index b131a5e..e76c2e4 100644
--- a/docs/source/generated/ramannoodle.structure.rst
+++ b/docs/source/generated/ramannoodle.structure.rst
@@ -1,32 +1,13 @@
-ramannoodle.structure package
-=============================
+ramannoodle.structure
+======================
 
-Submodules
-----------
-
-ramannoodle.structure.displace module
--------------------------------------
-
-.. automodule:: ramannoodle.structure.displace
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-ramannoodle.structure.reference module
---------------------------------------
-
-.. automodule:: ramannoodle.structure.reference
+.. automodule:: ramannoodle.structure
    :members:
    :undoc-members:
    :show-inheritance:
 
-ramannoodle.structure.symmetry\_utils module
---------------------------------------------
-
-.. automodule:: ramannoodle.structure.symmetry_utils
-   :members:
-   :undoc-members:
-   :show-inheritance:
+Submodules
+----------
 
 ramannoodle.structure.utils module
 ----------------------------------
@@ -35,11 +16,3 @@ ramannoodle.structure.utils module
    :members:
    :undoc-members:
    :show-inheritance:
-
-Module contents
----------------
-
-.. automodule:: ramannoodle.structure
-   :members:
-   :undoc-members:
-   :show-inheritance:
diff --git a/ramannoodle/dataset/torch/__init__.py b/ramannoodle/dataset/torch/__init__.py
index 36c1970..fa15ff2 100644
--- a/ramannoodle/dataset/torch/__init__.py
+++ b/ramannoodle/dataset/torch/__init__.py
@@ -1,7 +1,7 @@
 """Polarizability dataset for pytorch."""
 
 # flake8: noqa: F401
-from ramannoodle.dataset.torch.dataset import PolarizabilityDataset
+from ramannoodle.dataset.torch._dataset import PolarizabilityDataset
 from ramannoodle.dataset.torch import utils
 
 __all__ = ["PolarizabilityDataset", "utils"]
diff --git a/ramannoodle/dataset/torch/dataset.py b/ramannoodle/dataset/torch/_dataset.py
similarity index 100%
rename from ramannoodle/dataset/torch/dataset.py
rename to ramannoodle/dataset/torch/_dataset.py
diff --git a/ramannoodle/dynamics/__init__.py b/ramannoodle/dynamics/__init__.py
index 1923a5c..c3b57d6 100644
--- a/ramannoodle/dynamics/__init__.py
+++ b/ramannoodle/dynamics/__init__.py
@@ -1,7 +1,7 @@
 """Classes representing atomic motions, for example phonons and trajectories."""
 
 # flake8: noqa: F401
-from ramannoodle.dynamics.phonon import Phonons
-from ramannoodle.dynamics.trajectory import Trajectory
+from ramannoodle.dynamics._phonon import Phonons
+from ramannoodle.dynamics._trajectory import Trajectory
 
 __all__ = ["Phonons", "Trajectory"]
diff --git a/ramannoodle/dynamics/phonon.py b/ramannoodle/dynamics/_phonon.py
similarity index 98%
rename from ramannoodle/dynamics/phonon.py
rename to ramannoodle/dynamics/_phonon.py
index 80efb9e..7aa990e 100644
--- a/ramannoodle/dynamics/phonon.py
+++ b/ramannoodle/dynamics/_phonon.py
@@ -6,7 +6,7 @@
 from ramannoodle.abstract import Dynamics
 from ramannoodle.constants import RAMAN_TENSOR_CENTRAL_DIFFERENCE
 from ramannoodle.abstract import PolarizabilityModel
-from ramannoodle.spectrum.raman import PhononRamanSpectrum
+from ramannoodle.spectrum._raman import PhononRamanSpectrum
 from ramannoodle.exceptions import verify_ndarray_shape
 
 
diff --git a/ramannoodle/dynamics/trajectory.py b/ramannoodle/dynamics/_trajectory.py
similarity index 98%
rename from ramannoodle/dynamics/trajectory.py
rename to ramannoodle/dynamics/_trajectory.py
index 98819fe..9c88c9e 100644
--- a/ramannoodle/dynamics/trajectory.py
+++ b/ramannoodle/dynamics/_trajectory.py
@@ -9,7 +9,7 @@
 from ramannoodle.abstract import Dynamics
 from ramannoodle.abstract import PolarizabilityModel
 from ramannoodle.exceptions import verify_ndarray_shape, get_type_error
-from ramannoodle.spectrum.raman import MDRamanSpectrum
+from ramannoodle.spectrum._raman import MDRamanSpectrum
 from ramannoodle.structure.utils import apply_pbc
 
 
diff --git a/ramannoodle/io/__init__.py b/ramannoodle/io/__init__.py
index 049b897..b474940 100644
--- a/ramannoodle/io/__init__.py
+++ b/ramannoodle/io/__init__.py
@@ -3,6 +3,6 @@
 # flake8: noqa: F401
 from ramannoodle.io import generic
 from ramannoodle.io import vasp
-from ramannoodle.io import utils
+from ramannoodle.io import _utils
 
-__all__ = ["generic", "vasp", "utils"]
+__all__ = ["generic", "vasp", "_utils"]
diff --git a/ramannoodle/io/utils.py b/ramannoodle/io/_utils.py
similarity index 98%
rename from ramannoodle/io/utils.py
rename to ramannoodle/io/_utils.py
index 869ee78..5d39667 100644
--- a/ramannoodle/io/utils.py
+++ b/ramannoodle/io/_utils.py
@@ -20,7 +20,7 @@
 
 TORCH_PRESENT = True
 try:
-    from ramannoodle.dataset.torch.dataset import PolarizabilityDataset
+    from ramannoodle.dataset.torch._dataset import PolarizabilityDataset
 except UserError:
     TORCH_PRESENT = False
 
diff --git a/ramannoodle/io/generic.py b/ramannoodle/io/generic.py
index e46930e..edc1ec5 100644
--- a/ramannoodle/io/generic.py
+++ b/ramannoodle/io/generic.py
@@ -13,15 +13,15 @@
 import numpy as np
 from numpy.typing import NDArray
 
-from ramannoodle.dynamics.phonon import Phonons
-from ramannoodle.dynamics.trajectory import Trajectory
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.dynamics._phonon import Phonons
+from ramannoodle.dynamics._trajectory import Trajectory
+from ramannoodle.structure._reference import ReferenceStructure
 from ramannoodle.exceptions import UserError, get_torch_missing_error
 import ramannoodle.io.vasp as vasp_io
 
 TORCH_PRESENT = True
 try:
-    from ramannoodle.dataset.torch.dataset import PolarizabilityDataset
+    from ramannoodle.dataset.torch._dataset import PolarizabilityDataset
 except UserError:
     TORCH_PRESENT = False
 
diff --git a/ramannoodle/io/pymatgen/pymatgen.py b/ramannoodle/io/pymatgen.py
similarity index 96%
rename from ramannoodle/io/pymatgen/pymatgen.py
rename to ramannoodle/io/pymatgen.py
index 6f1fa79..81d0354 100644
--- a/ramannoodle/io/pymatgen/pymatgen.py
+++ b/ramannoodle/io/pymatgen.py
@@ -11,11 +11,11 @@
     get_type_error,
     IncompatibleStructureException,
 )
-from ramannoodle.dynamics.trajectory import Trajectory
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.dynamics._trajectory import Trajectory
+from ramannoodle.structure._reference import ReferenceStructure
 
 try:
-    from ramannoodle.dataset.torch.dataset import PolarizabilityDataset
+    from ramannoodle.dataset.torch._dataset import PolarizabilityDataset
 except UserError:
     pass
 
diff --git a/ramannoodle/io/pymatgen/__init__.py b/ramannoodle/io/pymatgen/__init__.py
deleted file mode 100644
index 5d078bd..0000000
--- a/ramannoodle/io/pymatgen/__init__.py
+++ /dev/null
@@ -1,18 +0,0 @@
-"""Functions for interacting with pymatgen."""
-
-# flake8: noqa: F401
-from ramannoodle.io.pymatgen.pymatgen import (
-    get_positions,
-    get_structure,
-    construct_polarizability_dataset,
-    construct_ref_structure,
-    construct_trajectory,
-)
-
-__all__ = [
-    "get_positions",
-    "get_structure",
-    "construct_polarizability_dataset",
-    "construct_ref_structure",
-    "construct_trajectory",
-]
diff --git a/ramannoodle/io/vasp/outcar.py b/ramannoodle/io/vasp/outcar.py
index 99c721f..b552dea 100644
--- a/ramannoodle/io/vasp/outcar.py
+++ b/ramannoodle/io/vasp/outcar.py
@@ -5,7 +5,7 @@
 import numpy as np
 from numpy.typing import NDArray
 
-from ramannoodle.io.utils import (
+from ramannoodle.io._utils import (
     _skip_file_until_line_contains,
     pathify,
     _read_polarizability_dataset,
@@ -13,12 +13,12 @@
 from ramannoodle.exceptions import InvalidFileException, NoMatchingLineFoundException
 from ramannoodle.constants import ATOMIC_WEIGHTS, ATOMIC_NUMBERS
 from ramannoodle.exceptions import get_type_error, UserError
-from ramannoodle.dynamics.phonon import Phonons
-from ramannoodle.dynamics.trajectory import Trajectory
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.dynamics._phonon import Phonons
+from ramannoodle.dynamics._trajectory import Trajectory
+from ramannoodle.structure._reference import ReferenceStructure
 
 try:
-    from ramannoodle.dataset.torch.dataset import PolarizabilityDataset
+    from ramannoodle.dataset.torch._dataset import PolarizabilityDataset
 except UserError:
     pass
 
diff --git a/ramannoodle/io/vasp/poscar.py b/ramannoodle/io/vasp/poscar.py
index 607eb58..24b1ff0 100644
--- a/ramannoodle/io/vasp/poscar.py
+++ b/ramannoodle/io/vasp/poscar.py
@@ -6,10 +6,10 @@
 import numpy as np
 from numpy.typing import NDArray
 
-from ramannoodle.io.utils import verify_structure, pathify
+from ramannoodle.io._utils import verify_structure, pathify
 from ramannoodle.exceptions import InvalidFileException
 from ramannoodle.constants import ATOM_SYMBOLS, ATOMIC_NUMBERS
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.structure._reference import ReferenceStructure
 from ramannoodle.io.vasp.outcar import _get_lattice_vector_from_outcar_line
 
 # pylint: disable=R0801
diff --git a/ramannoodle/io/vasp/vasprun.py b/ramannoodle/io/vasp/vasprun.py
index 5cbcee1..5559f8c 100644
--- a/ramannoodle/io/vasp/vasprun.py
+++ b/ramannoodle/io/vasp/vasprun.py
@@ -8,15 +8,15 @@
 import numpy as np
 from numpy.typing import NDArray
 
-from ramannoodle.io.utils import pathify, _read_polarizability_dataset
+from ramannoodle.io._utils import pathify, _read_polarizability_dataset
 from ramannoodle.exceptions import InvalidFileException, UserError
 from ramannoodle.constants import ATOMIC_WEIGHTS, ATOMIC_NUMBERS
-from ramannoodle.dynamics.phonon import Phonons
-from ramannoodle.dynamics.trajectory import Trajectory
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.dynamics._phonon import Phonons
+from ramannoodle.dynamics._trajectory import Trajectory
+from ramannoodle.structure._reference import ReferenceStructure
 
 try:
-    from ramannoodle.dataset.torch.dataset import PolarizabilityDataset
+    from ramannoodle.dataset.torch._dataset import PolarizabilityDataset
 except UserError:
     pass
 
diff --git a/ramannoodle/io/vasp/xdatcar.py b/ramannoodle/io/vasp/xdatcar.py
index a72d379..dd2e6fc 100644
--- a/ramannoodle/io/vasp/xdatcar.py
+++ b/ramannoodle/io/vasp/xdatcar.py
@@ -5,9 +5,9 @@
 import numpy as np
 from numpy.typing import NDArray
 
-from ramannoodle.io.utils import verify_trajectory, pathify
+from ramannoodle.io._utils import verify_trajectory, pathify
 from ramannoodle.exceptions import InvalidFileException
-from ramannoodle.dynamics.trajectory import Trajectory
+from ramannoodle.dynamics._trajectory import Trajectory
 from ramannoodle.io.vasp.poscar import (
     _read_lattice,
     _read_atomic_symbols,
diff --git a/ramannoodle/pmodel/__init__.py b/ramannoodle/pmodel/__init__.py
index b3a033f..2c2c895 100644
--- a/ramannoodle/pmodel/__init__.py
+++ b/ramannoodle/pmodel/__init__.py
@@ -1,7 +1,7 @@
 """Classes for various polarizability models."""
 
 # flake8: noqa: F401
-from ramannoodle.pmodel.art import ARTModel
-from ramannoodle.pmodel.interpolation import InterpolationModel
+from ramannoodle.pmodel._art import ARTModel
+from ramannoodle.pmodel._interpolation import InterpolationModel
 
 __all__ = ["ARTModel", "InterpolationModel"]
diff --git a/ramannoodle/pmodel/art.py b/ramannoodle/pmodel/_art.py
similarity index 99%
rename from ramannoodle/pmodel/art.py
rename to ramannoodle/pmodel/_art.py
index b3f0f74..8e5752d 100644
--- a/ramannoodle/pmodel/art.py
+++ b/ramannoodle/pmodel/_art.py
@@ -9,7 +9,7 @@
 from tabulate import tabulate
 
 from ramannoodle.constants import ANSICOLORS
-from ramannoodle.pmodel.interpolation import InterpolationModel
+from ramannoodle.pmodel._interpolation import InterpolationModel
 from ramannoodle.exceptions import (
     get_type_error,
     get_shape_error,
diff --git a/ramannoodle/pmodel/interpolation.py b/ramannoodle/pmodel/_interpolation.py
similarity index 99%
rename from ramannoodle/pmodel/interpolation.py
rename to ramannoodle/pmodel/_interpolation.py
index 2a18fec..3714a8a 100644
--- a/ramannoodle/pmodel/interpolation.py
+++ b/ramannoodle/pmodel/_interpolation.py
@@ -17,11 +17,11 @@
 from ramannoodle.constants import ANSICOLORS
 from ramannoodle.abstract import PolarizabilityModel
 from ramannoodle.structure.utils import calc_displacement
-from ramannoodle.structure.symmetry_utils import (
+from ramannoodle.structure._symmetry_utils import (
     is_orthogonal_to_all,
     is_collinear_with_all,
 )
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.structure._reference import ReferenceStructure
 from ramannoodle.exceptions import (
     InvalidDOFException,
     get_type_error,
@@ -33,7 +33,7 @@
     get_pymatgen_missing_error,
 )
 import ramannoodle.io.generic as generic_io
-from ramannoodle.io.utils import pathify_as_list
+from ramannoodle.io._utils import pathify_as_list
 
 PYMATGEN_PRESENT = True
 try:
diff --git a/ramannoodle/pmodel/torch/__init__.py b/ramannoodle/pmodel/torch/__init__.py
index e8550f6..f6183e7 100644
--- a/ramannoodle/pmodel/torch/__init__.py
+++ b/ramannoodle/pmodel/torch/__init__.py
@@ -1,8 +1,8 @@
 """Modules for polarizability models implemented with PyTorch."""
 
 # flake8: noqa: F401
-from ramannoodle.pmodel.torch.gnn import PotGNN
-from ramannoodle.pmodel.torch.train import train_single_epoch
-from ramannoodle.pmodel.torch import utils
+from ramannoodle.pmodel.torch._gnn import PotGNN
+from ramannoodle.pmodel.torch._train import train_single_epoch
+from ramannoodle.pmodel.torch import _utils
 
-__all__ = ["PotGNN", "train_single_epoch", "utils"]
+__all__ = ["PotGNN", "train_single_epoch", "_utils"]
diff --git a/ramannoodle/pmodel/torch/gnn.py b/ramannoodle/pmodel/torch/_gnn.py
similarity index 99%
rename from ramannoodle/pmodel/torch/gnn.py
rename to ramannoodle/pmodel/torch/_gnn.py
index 82d3e5e..802d11b 100644
--- a/ramannoodle/pmodel/torch/gnn.py
+++ b/ramannoodle/pmodel/torch/_gnn.py
@@ -9,7 +9,7 @@
 
 from tqdm import tqdm
 
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.structure._reference import ReferenceStructure
 from ramannoodle.exceptions import (
     get_torch_missing_error,
     UserError,
@@ -32,7 +32,7 @@
     from torch_geometric.nn.inits import reset
     from torch_geometric.nn.models.schnet import ShiftedSoftplus
     from torch_geometric.utils import scatter
-    import ramannoodle.pmodel.torch.utils as rn_torch_utils
+    import ramannoodle.pmodel.torch._utils as rn_torch_utils
     from ramannoodle.dataset.torch.utils import (
         polarizability_vectors_to_tensors,
         polarizability_tensors_to_vectors,
diff --git a/ramannoodle/pmodel/torch/train.py b/ramannoodle/pmodel/torch/_train.py
similarity index 95%
rename from ramannoodle/pmodel/torch/train.py
rename to ramannoodle/pmodel/torch/_train.py
index 309ea5e..5af8808 100644
--- a/ramannoodle/pmodel/torch/train.py
+++ b/ramannoodle/pmodel/torch/_train.py
@@ -14,8 +14,8 @@
     from torch.nn.modules.loss import _Loss
     from torch.optim.optimizer import Optimizer
 
-    from ramannoodle.pmodel.torch.gnn import PotGNN
-    from ramannoodle.dataset.torch.dataset import PolarizabilityDataset
+    from ramannoodle.pmodel.torch._gnn import PotGNN
+    from ramannoodle.dataset.torch._dataset import PolarizabilityDataset
 except (ModuleNotFoundError, UserError) as exc:
     raise get_torch_missing_error() from exc
 
diff --git a/ramannoodle/pmodel/torch/utils.py b/ramannoodle/pmodel/torch/_utils.py
similarity index 100%
rename from ramannoodle/pmodel/torch/utils.py
rename to ramannoodle/pmodel/torch/_utils.py
diff --git a/ramannoodle/spectrum/__init__.py b/ramannoodle/spectrum/__init__.py
index b8adbb4..1f3862f 100644
--- a/ramannoodle/spectrum/__init__.py
+++ b/ramannoodle/spectrum/__init__.py
@@ -1,7 +1,7 @@
 """Classes and functions for calculating and manipulating spectra."""
 
 # flake8: noqa: F401
-from ramannoodle.spectrum.raman import (
+from ramannoodle.spectrum._raman import (
     PhononRamanSpectrum,
     MDRamanSpectrum,
 )
diff --git a/ramannoodle/spectrum/raman.py b/ramannoodle/spectrum/_raman.py
similarity index 100%
rename from ramannoodle/spectrum/raman.py
rename to ramannoodle/spectrum/_raman.py
diff --git a/ramannoodle/structure/__init__.py b/ramannoodle/structure/__init__.py
index 78da70c..9c23428 100644
--- a/ramannoodle/structure/__init__.py
+++ b/ramannoodle/structure/__init__.py
@@ -1,14 +1,14 @@
 """Classes and functions for atomic structures and structural symmetries."""
 
 # flake8: noqa: F401
-from ramannoodle.structure.reference import ReferenceStructure
-from ramannoodle.structure.displace import (
+from ramannoodle.structure._reference import ReferenceStructure
+from ramannoodle.structure._displace import (
     write_ast_displaced_structures,
     get_ast_displaced_positions,
     write_displaced_structures,
     get_displaced_positions,
 )
-from ramannoodle.structure import symmetry_utils
+from ramannoodle.structure import _symmetry_utils
 from ramannoodle.structure import utils
 
 
@@ -18,6 +18,6 @@
     "get_ast_displaced_positions",
     "write_displaced_structures",
     "get_displaced_positions",
-    "symmetry_utils",
+    "_symmetry_utils",
     "utils",
 ]
diff --git a/ramannoodle/structure/displace.py b/ramannoodle/structure/_displace.py
similarity index 98%
rename from ramannoodle/structure/displace.py
rename to ramannoodle/structure/_displace.py
index 6d10f69..184450c 100644
--- a/ramannoodle/structure/displace.py
+++ b/ramannoodle/structure/_displace.py
@@ -16,7 +16,7 @@
 import numpy as np
 from numpy.typing import NDArray
 
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.structure._reference import ReferenceStructure
 from ramannoodle.structure.utils import displace_positions
 from ramannoodle.exceptions import (
     get_type_error,
@@ -24,7 +24,7 @@
     verify_ndarray_shape,
     verify_list_len,
 )
-from ramannoodle.io.utils import pathify_as_list
+from ramannoodle.io._utils import pathify_as_list
 import ramannoodle.io.generic as generic_io
 
 
diff --git a/ramannoodle/structure/reference.py b/ramannoodle/structure/_reference.py
similarity index 98%
rename from ramannoodle/structure/reference.py
rename to ramannoodle/structure/_reference.py
index 36bac1f..98cbd7e 100644
--- a/ramannoodle/structure/reference.py
+++ b/ramannoodle/structure/_reference.py
@@ -20,7 +20,7 @@
     calc_displacement,
 )
 from ramannoodle.constants import ATOM_SYMBOLS
-from ramannoodle.structure import symmetry_utils
+from ramannoodle.structure import _symmetry_utils
 
 
 def _compute_permutation_matrices(
@@ -223,7 +223,7 @@ def get_equivalent_displacements(
                 candidate_positions,
             )
 
-            orthogonal_result = symmetry_utils.is_orthogonal_to_all(
+            orthogonal_result = _symmetry_utils.is_orthogonal_to_all(
                 candidate_displacement.flatten(),
                 [item.flatten() for item in orthogonal_displacements],
             )
@@ -239,7 +239,7 @@ def get_equivalent_displacements(
                 orthogonal_displacements.append(candidate_displacement)
             else:
                 # Candidate can be collinear to maximum of one orthogonal vector
-                collinear_result = symmetry_utils.is_non_collinear_with_all(
+                collinear_result = _symmetry_utils.is_non_collinear_with_all(
                     candidate_displacement,
                     orthogonal_displacements,
                 )
diff --git a/ramannoodle/structure/symmetry_utils.py b/ramannoodle/structure/_symmetry_utils.py
similarity index 100%
rename from ramannoodle/structure/symmetry_utils.py
rename to ramannoodle/structure/_symmetry_utils.py
diff --git a/test/conftest.py b/test/conftest.py
index e24f9f9..6939b4d 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -7,7 +7,7 @@
 import pytest
 from pytest import FixtureRequest
 
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.structure._reference import ReferenceStructure
 import ramannoodle.io.generic as generic_io
 
 
diff --git a/test/tests/test_art.py b/test/tests/test_art.py
index 3f6a9ef..3c665bf 100644
--- a/test/tests/test_art.py
+++ b/test/tests/test_art.py
@@ -7,9 +7,9 @@
 from numpy.typing import NDArray
 import pytest
 
-from ramannoodle.pmodel.art import ARTModel
+from ramannoodle.pmodel._art import ARTModel
 from ramannoodle.exceptions import InvalidDOFException, UserError
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.structure._reference import ReferenceStructure
 
 # pylint: disable=protected-access
 # pylint: disable=too-many-arguments
diff --git a/test/tests/test_displace.py b/test/tests/test_displace.py
index 02f8949..b4cf098 100644
--- a/test/tests/test_displace.py
+++ b/test/tests/test_displace.py
@@ -9,8 +9,8 @@
 import pytest
 
 import ramannoodle.io.vasp as vasp_io
-from ramannoodle.structure.reference import ReferenceStructure
-from ramannoodle.structure.displace import (
+from ramannoodle.structure._reference import ReferenceStructure
+from ramannoodle.structure._displace import (
     write_ast_displaced_structures,
     get_ast_displaced_positions,
     write_displaced_structures,
diff --git a/test/tests/test_interpolation.py b/test/tests/test_interpolation.py
index 8679214..8d7800a 100644
--- a/test/tests/test_interpolation.py
+++ b/test/tests/test_interpolation.py
@@ -7,10 +7,10 @@
 from numpy.typing import NDArray
 import pytest
 
-from ramannoodle.pmodel.interpolation import find_duplicates
-from ramannoodle.pmodel.interpolation import InterpolationModel
+from ramannoodle.pmodel._interpolation import find_duplicates
+from ramannoodle.pmodel._interpolation import InterpolationModel
 from ramannoodle.exceptions import InvalidDOFException
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.structure._reference import ReferenceStructure
 
 # pylint: disable=protected-access
 # pylint: disable=too-many-arguments
diff --git a/test/tests/test_io.py b/test/tests/test_io.py
index ad581a0..f2c9c56 100644
--- a/test/tests/test_io.py
+++ b/test/tests/test_io.py
@@ -9,7 +9,7 @@
 from numpy.typing import NDArray
 
 import ramannoodle.io.generic as generic_io
-from ramannoodle.io.utils import pathify_as_list
+from ramannoodle.io._utils import pathify_as_list
 
 
 @pytest.mark.parametrize(
diff --git a/test/tests/test_phonon.py b/test/tests/test_phonon.py
index 3fd6de0..71bd237 100644
--- a/test/tests/test_phonon.py
+++ b/test/tests/test_phonon.py
@@ -9,9 +9,9 @@
 import pytest
 
 import ramannoodle.io.generic
-from ramannoodle.dynamics.phonon import Phonons
-from ramannoodle.structure.reference import ReferenceStructure
-from ramannoodle.pmodel.art import ARTModel
+from ramannoodle.dynamics._phonon import Phonons
+from ramannoodle.structure._reference import ReferenceStructure
+from ramannoodle.pmodel._art import ARTModel
 
 
 @pytest.mark.parametrize(
diff --git a/test/tests/test_phonon_spectrum.py b/test/tests/test_phonon_spectrum.py
index c43d5da..8b337a1 100644
--- a/test/tests/test_phonon_spectrum.py
+++ b/test/tests/test_phonon_spectrum.py
@@ -12,13 +12,13 @@
 import pymatgen.core
 
 import ramannoodle.io.generic
-from ramannoodle.pmodel.interpolation import InterpolationModel
-from ramannoodle.pmodel.art import ARTModel
-from ramannoodle.spectrum.raman import (
+from ramannoodle.pmodel._interpolation import InterpolationModel
+from ramannoodle.pmodel._art import ARTModel
+from ramannoodle.spectrum._raman import (
     get_bose_einstein_correction,
     get_laser_correction,
 )
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.structure._reference import ReferenceStructure
 from ramannoodle.spectrum.utils import convolve_spectrum
 
 # pylint: disable=protected-access,too-many-locals
diff --git a/test/tests/test_structure.py b/test/tests/test_structure.py
index 6184038..422ab3a 100644
--- a/test/tests/test_structure.py
+++ b/test/tests/test_structure.py
@@ -8,13 +8,13 @@
 import pytest
 
 from ramannoodle.structure.utils import apply_pbc, apply_pbc_displacement
-from ramannoodle.structure.symmetry_utils import (
+from ramannoodle.structure._symmetry_utils import (
     is_collinear_with_all,
     is_non_collinear_with_all,
     are_collinear,
     is_orthogonal_to_all,
 )
-from ramannoodle.structure.reference import (
+from ramannoodle.structure._reference import (
     ReferenceStructure,
     _get_positions_permutation_matrix,
 )
diff --git a/test/tests/test_trajectory_spectrum.py b/test/tests/test_trajectory_spectrum.py
index b95c948..75428dc 100644
--- a/test/tests/test_trajectory_spectrum.py
+++ b/test/tests/test_trajectory_spectrum.py
@@ -9,8 +9,8 @@
 
 from ramannoodle.spectrum.utils import calc_signal_spectrum
 import ramannoodle.io.vasp as vasp_io
-from ramannoodle.structure.reference import ReferenceStructure
-from ramannoodle.pmodel.interpolation import InterpolationModel
+from ramannoodle.structure._reference import ReferenceStructure
+from ramannoodle.pmodel._interpolation import InterpolationModel
 
 # pylint: disable=R0801
 
diff --git a/test/tests/torch/test_dataset.py b/test/tests/torch/test_dataset.py
index 1e6bffd..42961a8 100644
--- a/test/tests/torch/test_dataset.py
+++ b/test/tests/torch/test_dataset.py
@@ -9,7 +9,7 @@
 import pytest
 
 import ramannoodle.io.generic as generic_io
-from ramannoodle.dataset.torch.dataset import PolarizabilityDataset
+from ramannoodle.dataset.torch._dataset import PolarizabilityDataset
 
 
 @pytest.mark.parametrize(
diff --git a/test/tests/torch/test_gnn.py b/test/tests/torch/test_gnn.py
index 6911370..801a6fa 100644
--- a/test/tests/torch/test_gnn.py
+++ b/test/tests/torch/test_gnn.py
@@ -12,18 +12,18 @@
 
 import ramannoodle as rn
 
-from ramannoodle.pmodel.torch.gnn import PotGNN
-from ramannoodle.pmodel.torch.utils import (
+from ramannoodle.pmodel.torch._gnn import PotGNN
+from ramannoodle.pmodel.torch._utils import (
     _radius_graph_pbc,
     get_rotations,
 )
-from ramannoodle.pmodel.torch.train import train_single_epoch
+from ramannoodle.pmodel.torch._train import train_single_epoch
 from ramannoodle.dataset.torch.utils import polarizability_vectors_to_tensors
 
 
 # import ramannoodle.io.vasp as vasp_io
 # from ramannoodle.structure.structure_utils import apply_pbc
-from ramannoodle.structure.reference import ReferenceStructure
+from ramannoodle.structure._reference import ReferenceStructure
 
 # pylint: disable=protected-access, too-many-arguments, not-callable