diff --git a/hexrd/__init__.py b/hexrd/__init__.py
index 69e3be50d..77ca9b19a 100644
--- a/hexrd/__init__.py
+++ b/hexrd/__init__.py
@@ -1 +1,32 @@
-__path__ = __import__('pkgutil').extend_path(__path__, __name__)
+import importlib
+import sys
+
+from .material import crystallography
+from .material import jcpds
+from .material import mksupport
+from .material import spacegroup
+from .material import symbols
+from .material import symmetry
+from .material import unitcell
+
+# These are aliases for import paths, so we don't break old HEXRD scripts.
+# We will verify the alias files *do not* exist, to avoid confusion.
+module_aliases = {
+    'hexrd.crystallography': crystallography,
+    'hexrd.mksupport': mksupport,
+    'hexrd.spacegroup': spacegroup,
+    'hexrd.symbols': symbols,
+    'hexrd.symmetry': symmetry,
+    'hexrd.unitcell': unitcell,
+}
+
+for alias, module in module_aliases.items():
+    try:
+        file_exists = importlib.import_module(alias).__name__ == alias
+    except ImportError:
+        file_exists = False
+
+    if file_exists:
+        raise Exception(f'"{alias}" is an alias path and should not exist')
+
+    sys.modules[alias] = module
diff --git a/hexrd/instrument/detector.py b/hexrd/instrument/detector.py
index fbd3763e1..6cfd8d834 100644
--- a/hexrd/instrument/detector.py
+++ b/hexrd/instrument/detector.py
@@ -9,7 +9,7 @@
 from hexrd import matrixutil as mutil
 from hexrd import xrdutil
 
-from hexrd.crystallography import PlaneData
+from hexrd.material.crystallography import PlaneData
 
 from hexrd.transforms.xfcapi import (
     detectorXYToGvec,
diff --git a/hexrd/instrument/hedm_instrument.py b/hexrd/instrument/hedm_instrument.py
index 1ecbdb5ae..8f5b7e1ca 100644
--- a/hexrd/instrument/hedm_instrument.py
+++ b/hexrd/instrument/hedm_instrument.py
@@ -68,7 +68,7 @@
     unitRowVector,
 )
 from hexrd import xrdutil
-from hexrd.crystallography import PlaneData
+from hexrd.material.crystallography import PlaneData
 from hexrd import constants as ct
 from hexrd.rotations import (
     angleAxisOfRotMat,
diff --git a/hexrd/material/__init__.py b/hexrd/material/__init__.py
new file mode 100644
index 000000000..8ce0c9625
--- /dev/null
+++ b/hexrd/material/__init__.py
@@ -0,0 +1,7 @@
+from .material import (
+    _angstroms,
+    _kev,
+    load_materials_hdf5,
+    Material,
+    save_materials_hdf5,
+)
diff --git a/hexrd/crystallography.py b/hexrd/material/crystallography.py
similarity index 100%
rename from hexrd/crystallography.py
rename to hexrd/material/crystallography.py
diff --git a/hexrd/material/jcpds.py b/hexrd/material/jcpds.py
new file mode 100644
index 000000000..0affaa8c2
--- /dev/null
+++ b/hexrd/material/jcpds.py
@@ -0,0 +1,327 @@
+import os
+import numpy as np
+
+
+class JCPDS_extend():
+    def __init__(self, filename=None):
+        self.a0 = 0
+        self.b0 = 0
+        self.c0 = 0
+        self.alpha0 = 0
+        self.beta0 = 0
+        self.gamma0 = 0
+        self.v0 = 0
+
+        self.k0 = 100
+        self.k0p = 5  # k0p at 298K
+
+        self.dk0dt = 0
+        self.dk0pdt = 0
+
+        self.symmetry = ''
+        self.alpha_t = 0  # alphat at 298K
+        self.dalpha_t_dt = 0
+
+        self.file = ' '
+        self.name = ' '
+        self.version = 0
+        self.comments = ''
+
+        if filename:
+            self.file = filename
+            self.read_file(self.file)
+            self.update_v0()
+
+    def read_file(self, file):
+        """
+        read a jcpds file
+        """
+        self.file = file
+        # Construct base name = file without path and without extension
+        name = os.path.splitext(os.path.basename(self.file))[0]
+        self.name = name
+#       line = '', nd=0
+        version = 0.
+        self.comments = []
+        self.DiffLines = []
+
+        version_status = ''
+
+        inp = open(file, 'r').readlines()
+#       my_list = [] # get all the text first and throw into my_list
+
+        if inp[0][0] in ('2', '3', '4'):
+            version = int(inp[0])  # JCPDS version number
+            self.version = version
+            header = inp[1]  # header
+            self.comments = header
+
+            item = str.split(inp[2])
+            crystal_system = int(item[0])
+            if crystal_system == 1:
+                self.symmetry = 'cubic'
+            elif crystal_system == 2:
+                self.symmetry = 'hexagonal'
+            elif crystal_system == 3:
+                self.symmetry = 'tetragonal'
+            elif crystal_system == 4:
+                self.symmetry = 'orthorhombic'
+            elif crystal_system == 5:
+                self.symmetry = 'monoclinic'
+            elif crystal_system == 6:
+                self.symmetry = 'triclinic'
+            elif crystal_system == 7:
+                self.symmetry = 'manual'
+            # 1 cubic, 2 hexagonal, 3 tetragonal, 4 orthorhombic
+            # 5 monoclinic, 6 triclinic, 7 manual P, d-sp input
+
+            k0 = float(item[1])
+            k0p = float(item[2])
+            self.k0 = k0
+            self.k0p = k0p
+
+            item = str.split(inp[3])  # line for unit-cell parameters
+
+            if crystal_system == 1:  # cubic
+                a = float(item[0])
+                b = a
+                c = a
+                alpha = 90.
+                beta = 90.
+                gamma = 90.
+            elif crystal_system == 7:  # P, d-sp input
+                a = float(item[0])
+                b = a
+                c = a
+                alpha = 90.
+                beta = 90.
+                gamma = 90.
+            elif crystal_system == 2:  # hexagonal
+                a = float(item[0])
+                c = float(item[1])
+                b = a
+                alpha = 90.
+                beta = 90.
+                gamma = 120.
+            elif crystal_system == 3:  # tetragonal
+                a = float(item[0])
+                c = float(item[1])
+                b = a
+                alpha = 90.
+                beta = 90.
+                gamma = 90.
+            elif crystal_system == 4:  # orthorhombic
+                a = float(item[0])
+                b = float(item[1])
+                c = float(item[2])
+                alpha = 90.
+                beta = 90.
+                gamma = 90.
+            elif crystal_system == 5:  # monoclinic
+                a = float(item[0])
+                b = float(item[1])
+                c = float(item[2])
+                beta = float(item[3])
+                alpha = 90.
+                gamma = 90.
+            elif crystal_system == 6:  # triclinic
+                a = float(item[0])
+                b = float(item[1])
+                c = float(item[2])
+                alpha = float(item[3])
+                beta = float(item[4])
+                gamma = float(item[5])
+
+            self.a0 = a
+            self.b0 = b
+            self.c0 = c
+            self.alpha0 = alpha
+            self.beta0 = beta
+            self.gamma0 = gamma
+
+            item = str.split(inp[4])
+
+            if self.version == 3:
+                alpha_t = 0.
+            else:
+                alpha_t = float(item[0])
+            self.alpha_t = alpha_t
+
+            version_status = 'new'
+
+        elif 'VERSION' in inp[0]:
+            jcpdsfile = open(file, 'r')
+            while True:
+                jcpdsline = jcpdsfile.readline()
+                if jcpdsline == '':
+                    break
+
+                jlinespl = jcpdsline.split()
+
+                if jlinespl[0] == 'VERSION:':
+                    version = int(jlinespl[1])
+                    self.version = version
+
+                if jlinespl[0] == 'COMMENT:':
+                    header = ' '.join(jlinespl[1:])
+                    self.comments = header
+
+                if jlinespl[0] == 'K0:':
+                    k0 = float(jlinespl[1])
+                    self.k0 = k0
+
+                if jlinespl[0] == 'K0P:':
+                    k0p = float(jlinespl[1])
+                    self.k0p = k0p
+
+                if jlinespl[0] == 'DK0DT:':
+                    dk0dt = float(jlinespl[1])
+                    self.dk0dt = dk0dt
+
+                if jlinespl[0] == 'DK0PDT:':
+                    dk0pdt = float(jlinespl[1])
+                    self.dk0pdt = dk0pdt
+
+                if jlinespl[0] == 'SYMMETRY:':
+                    self.symmetry = jlinespl[1].lower()
+
+                if jlinespl[0] == 'A:':
+                    a = float(jlinespl[1])
+                    self.a0 = a
+
+                if jlinespl[0] == 'B:':
+                    b = float(jlinespl[1])
+                    self.b0 = b
+
+                if jlinespl[0] == 'C:':
+                    c = float(jlinespl[1])
+                    self.c0 = c
+
+                if jlinespl[0] == 'ALPHA:':
+                    alpha = float(jlinespl[1])
+                    self.alpha0 = alpha
+
+                if jlinespl[0] == 'BETA:':
+                    beta = float(jlinespl[1])
+                    self.beta0 = beta
+
+                if jlinespl[0] == 'GAMMA:':
+                    gamma = float(jlinespl[1])
+                    self.gamma0 = gamma
+
+                if jlinespl[0] == 'VOLUME:':
+                    v = float(jlinespl[1])
+                    self.v0 = v
+
+                if jlinespl[0] == 'ALPHAT:':
+                    alphat = float(jlinespl[1])
+                    self.alpha_t = alphat
+
+                if jlinespl[0] == 'DALPHATDT:':
+                    dalphatdt = float(jlinespl[1])
+                    self.dalpha_t_dt = dalphatdt
+
+                if jlinespl[0] == 'DIHKL:':
+                    pass
+
+            if self.symmetry == 'cubic':
+                self.b0 = self.a0
+                self.c0 = self.a0
+                self.alpha0 = 90.
+                self.beta0 = 90.
+                self.gamma0 = 90.
+            elif self.symmetry == 'manual':
+                self.b0 = self.a0
+                self.c0 = self.a0
+                self.alpha0 = 90.
+                self.beta0 = 90.
+                self.gamma0 = 90.
+            elif self.symmetry == 'hexagonal' or self.symmetry == 'trigonal':
+                self.b0 = self.a0
+                self.alpha0 = 90.
+                self.beta0 = 90.
+                self.gamma0 = 120.
+            elif self.symmetry == 'tetragonal':
+                self.b0 = self.a0
+                self.alpha0 = 90.
+                self.beta0 = 90.
+                self.gamma0 = 90.
+            elif self.symmetry == 'orthorhombic':
+                self.alpha0 = 90.
+                self.beta0 = 90.
+                self.gamma0 = 90.
+            elif self.symmetry == 'monoclinic':
+                self.alpha0 = 90.
+                self.gamma0 = 90.
+            # elif self.symmetry == 'triclinic':
+
+            jcpdsfile.close()
+
+            version_status = 'new'
+
+        else:
+            version_status = 'old'
+
+        if version_status == 'new':
+            self.v = self.calc_volume_unitcell()
+
+    def verify_symmetry_match(self, mat):
+        if not self.symmetry_matches(mat):
+            msg = (
+                f'The JCPDS symmetry "{self.symmetry}" does not match the '
+                f'symmetry of the material "{mat.latticeType}"!'
+            )
+            raise SymmetryMismatch(msg)
+
+    def symmetry_matches(self, mat):
+        return mat.latticeType == self.symmetry
+
+    @property
+    def _lat_param_names(self):
+        return ['a0', 'b0', 'c0', 'alpha0', 'beta0', 'gamma0']
+
+    @property
+    def lattice_params(self):
+        return [getattr(self, x) for x in self._lat_param_names]
+
+    def matches_material(self, mat):
+        self.verify_symmetry_match(mat)
+        mat_lp = [x.value for x in mat.latticeParameters]
+        for v1, v2 in zip(mat_lp, self.lattice_params):
+            if not np.isclose(v1, v2):
+                return False
+
+    def write_lattice_params_to_material(self, mat):
+        self.verify_symmetry_match(mat)
+        mat.latticeParameters = self.lattice_params
+        mat.reset_v0()
+
+    def write_pt_params_to_material(self, mat):
+        self.verify_symmetry_match(mat)
+        mat.k0 = self.k0
+        mat.k0p = self.k0p
+        mat.dk0dt = self.dk0dt
+        mat.dk0pdt = self.dk0pdt
+        mat.alpha_t = self.alpha_t
+        mat.dalpha_t_dt = self.dalpha_t_dt
+
+    def update_v0(self):
+        self.v0 = self.calc_volume_unitcell()
+
+    def calc_volume_unitcell(self):
+        '''calculate volume of the unitcell
+        Ref: Introduction to conventional TEM
+        Marc De Graef, Appendix 1 pg. 662
+        '''
+        ca = np.cos(np.radians(self.alpha0))
+        cb = np.cos(np.radians(self.beta0))
+        cg = np.cos(np.radians(self.gamma0))
+
+        v0 = self.a0*self.b0*self.c0
+        f = np.sqrt(1 - ca**2 - cb**2 - cg**2
+                    + 2 * ca * cb * cg)
+        return v0*f
+
+
+class SymmetryMismatch(Exception):
+    pass
diff --git a/hexrd/material.py b/hexrd/material/material.py
similarity index 62%
rename from hexrd/material.py
rename to hexrd/material/material.py
index 029cbe52b..abb139f08 100644
--- a/hexrd/material.py
+++ b/hexrd/material/material.py
@@ -32,23 +32,24 @@
 materials are defined by name in materialDict.
 """
 from configparser import SafeConfigParser as Parser
-import numpy
+import numpy as np
 
-from hexrd.crystallography import PlaneData as PData
+from hexrd.material.crystallography import PlaneData as PData
+from hexrd.material import symmetry, unitcell
 from hexrd.valunits import valWUnit
-from hexrd import unitcell
 from hexrd.constants import ptable, ptableinverse, chargestate
-from hexrd import symmetry
-import copy
 
 from os import path
 from pathlib import Path
 from CifFile import ReadCif
 import h5py
 from warnings import warn
-from hexrd.mksupport import Write2H5File
-from hexrd.symbols import xtal_sys_dict
-from hexrd.symbols import Hall_to_sgnum, HM_to_sgnum
+from hexrd.material.mksupport import Write2H5File
+from hexrd.material.symbols import (
+    xtal_sys_dict,
+    Hall_to_sgnum,
+    HM_to_sgnum,
+)
 from hexrd.utils.compatibility import h5py_read_string
 from hexrd.fitting.peakfunctions import _unit_gaussian
 
@@ -61,15 +62,15 @@
 
 
 def _angstroms(x):
-    return valWUnit('lp', 'length',  x, 'angstrom')
+    return valWUnit('lp', 'length', x, 'angstrom')
 
 
 def _degrees(x):
-    return valWUnit('lp', 'angle',   x, 'degrees')
+    return valWUnit('lp', 'angle', x, 'degrees')
 
 
 def _kev(x):
-    return valWUnit('xrayenergy', 'energy',  x, 'keV')
+    return valWUnit('xrayenergy', 'energy', x, 'keV')
 
 
 def _key(x):
@@ -88,17 +89,24 @@ class Material(object):
     space group data.
     default data is for nickel, but name is material
     """
+
     DFLT_NAME = 'material.xtal'
     DFLT_XTAL = 'Ni'
     DFLT_SGNUM = 225
 
-    DFLT_LPARMS = [_angstroms(3.61), _angstroms(3.61), _angstroms(3.61),
-                   _degrees(90.0), _degrees(90.0), _degrees(90.0)]
+    DFLT_LPARMS = [
+        _angstroms(3.61),
+        _angstroms(3.61),
+        _angstroms(3.61),
+        _degrees(90.0),
+        _degrees(90.0),
+        _degrees(90.0),
+    ]
     DFLT_SSMAX = 100
 
     DFLT_KEV = valWUnit('wavelength', 'energy', 80.725e0, 'keV')
     DFLT_STR = 0.0025
-    DFLT_TTH = numpy.radians(0.25)
+    DFLT_TTH = np.radians(0.25)
     DFLT_TTHMAX = None
     """
     ATOMINFO    Fractional Atom Position of an atom in the
@@ -108,10 +116,10 @@ class Material(object):
 
     ATOMTYPE    atomic number of all the different species in the unitcell
     """
-    DFLT_ATOMINFO = numpy.array([[0., 0., 0., 1.]])
-    DFLT_U = numpy.array([6.33E-3])
-    DFLT_ATOMTYPE = numpy.array([28])
-    DFLT_CHARGE = numpy.array(["0"])
+    DFLT_ATOMINFO = np.array([[0.0, 0.0, 0.0, 1.0]])
+    DFLT_U = np.array([6.33e-3])
+    DFLT_ATOMTYPE = np.array([28])
+    DFLT_CHARGE = np.array(["0"])
 
     '''
     the dmin parameter is used to figure out the maximum sampling for g-vectors
@@ -122,7 +130,7 @@ class Material(object):
     '''
     default stiffness tensor in voight notation
     '''
-    DFLT_STIFFNESS = numpy.eye(6)
+    DFLT_STIFFNESS = np.eye(6)
 
     '''
     some materials have more than one space group setting. for ex
@@ -134,11 +142,14 @@ class Material(object):
     '''
     DFLT_SGSETTING = 0
 
-    def __init__(self, name=None,
-                 material_file=None,
-                 dmin=None,
-                 kev=None,
-                 sgsetting=DFLT_SGSETTING):
+    def __init__(
+        self,
+        name=None,
+        material_file=None,
+        dmin=None,
+        kev=None,
+        sgsetting=DFLT_SGSETTING,
+    ):
         """Constructor for Material
 
         name -- (str) name of crystal
@@ -157,7 +168,6 @@ def __init__(self, name=None,
         self.sgsetting = sgsetting
 
         if material_file:
-
             # >> @ date 08/20/2020 SS removing dependence on hklmax
             if isinstance(material_file, (Path, str)):
                 form = Path(material_file).suffix[1:]
@@ -193,6 +203,15 @@ def __init__(self, name=None,
             self._dmin = Material.DFLT_DMIN
             self._beamEnergy = Material.DFLT_KEV
 
+            self.pressure = 0
+            self.temperature = 298
+            self.k0 = 100.0
+            self.k0p = 0.0
+            self.dk0dt = 0.0
+            self.dk0pdt = 0.0
+            self.alpha_t = 0.0
+            self.dalpha_t_dt = 0.0
+
         # If these were specified, they override any other method of
         # obtaining them (including loading them from files).
         if dmin is not None:
@@ -202,6 +221,10 @@ def __init__(self, name=None,
             self._beamEnergy = kev
 
         self._newUnitcell()
+
+        if not hasattr(self, 'v0'):
+            self.reset_v0()
+
         self._newPdata()
         self.update_structure_factor()
 
@@ -227,7 +250,7 @@ def _reset_lparms(self):
         lparms = unitcell._rqpDict[ltype][1](lparms)
         lparms_vu = []
         for i in range(6):
-            if(i < 3):
+            if i < 3:
                 lparms_vu.append(_angstroms(lparms[i]))
             else:
                 lparms_vu.append(_degrees(lparms[i]))
@@ -243,10 +266,16 @@ def _newUnitcell(self):
         """
         self._reset_lparms()
         self._unitcell = unitcell.unitcell(
-            self._lparms, self.sgnum, self._atomtype, self._charge,
-            self._atominfo, self._U,
-            self._dmin.getVal('nm'), self._beamEnergy.value,
-            self._sgsetting)
+            self._lparms,
+            self.sgnum,
+            self._atomtype,
+            self._charge,
+            self._atominfo,
+            self._U,
+            self._dmin.getVal('nm'),
+            self._beamEnergy.value,
+            self._sgsetting,
+        )
 
         if hasattr(self, 'stiffness'):
             self._unitcell.stiffness = self.stiffness
@@ -288,10 +317,10 @@ def _newPdata(self):
             old_indices, new_indices = map_hkls_to_new(old_pdata, self._pData)
 
             # Map the new exclusions to the old ones. New ones default to True.
-            exclusions = numpy.ones(hkls.shape[1], dtype=bool)
+            exclusions = np.ones(hkls.shape[1], dtype=bool)
             exclusions[new_indices] = old_pdata.exclusions[old_indices]
 
-            if numpy.all(exclusions):
+            if np.all(exclusions):
                 # If they are all excluded, just set the default exclusions
                 self.set_default_exclusions()
             else:
@@ -300,13 +329,21 @@ def _newPdata(self):
             # Make the PlaneData object from scratch...
             lprm = self.reduced_lattice_parameters
             laue = self.unitcell._laueGroup
-            self._pData = PData(hkls, lprm, laue,
-                                self._beamEnergy, Material.DFLT_STR,
-                                tThWidth=self._tThWidth,
-                                tThMax=Material.DFLT_TTHMAX)
+            self._pData = PData(
+                hkls,
+                lprm,
+                laue,
+                self._beamEnergy,
+                Material.DFLT_STR,
+                tThWidth=self._tThWidth,
+                tThMax=Material.DFLT_TTHMAX,
+            )
 
             self.set_default_exclusions()
 
+    def reset_v0(self):
+        self.v0 = self.vol
+
     def set_default_exclusions(self):
         if hasattr(self, 'hkl_from_file'):
             # If we loaded hkls from the file, use those
@@ -318,7 +355,7 @@ def set_default_exclusions(self):
     def enable_hkls_from_file(self):
         # Enable hkls from the file
         # 'hkl_from_file' must be an attribute on `self`
-        exclusions = numpy.ones_like(self._pData.exclusions, dtype=bool)
+        exclusions = np.ones_like(self._pData.exclusions, dtype=bool)
         for i, g in enumerate(self._pData.hklDataList):
             if g['hkl'].tolist() in self.hkl_from_file.tolist():
                 exclusions[i] = False
@@ -327,40 +364,43 @@ def enable_hkls_from_file(self):
 
     def enable_hkls_below_index(self, index=5):
         # Enable hkls with indices less than @index
-        exclusions = numpy.ones_like(self._pData.exclusions, dtype=bool)
+        exclusions = np.ones_like(self._pData.exclusions, dtype=bool)
         exclusions[:index] = False
         self._pData.exclusions = exclusions
 
     def enable_hkls_below_tth(self, tth_threshold=90.0):
         '''
-          enable reflections with two-theta less than @tth_threshold degrees
+        enable reflections with two-theta less than @tth_threshold degrees
         '''
-        tth_threshold = numpy.radians(tth_threshold)
+        tth_threshold = np.radians(tth_threshold)
 
-        tth = numpy.array([hkldata['tTheta']
-                           for hkldata in self._pData.hklDataList])
-        dflt_excl = numpy.ones(tth.shape, dtype=numpy.bool)
-        dflt_excl2 = numpy.ones(tth.shape, dtype=numpy.bool)
+        tth = np.array(
+            [hkldata['tTheta'] for hkldata in self._pData.hklDataList]
+        )
+        dflt_excl = np.ones(tth.shape, dtype=np.bool)
+        dflt_excl2 = np.ones(tth.shape, dtype=np.bool)
 
-        if(hasattr(self, 'hkl_from_file')):
+        if hasattr(self, 'hkl_from_file'):
             """
             hkls were read from the file so the exclusions will be set
             based on what hkls are present
             """
             for i, g in enumerate(self._pData.hklDataList):
-                if(g['hkl'].tolist() in self.hkl_from_file.tolist()):
+                if g['hkl'].tolist() in self.hkl_from_file.tolist():
                     dflt_excl[i] = False
 
-            dflt_excl2[~numpy.isnan(tth)] = \
-                ~((tth[~numpy.isnan(tth)] >= 0.0) &
-                  (tth[~numpy.isnan(tth)] <= tth_threshold))
+            dflt_excl2[~np.isnan(tth)] = ~(
+                (tth[~np.isnan(tth)] >= 0.0)
+                & (tth[~np.isnan(tth)] <= tth_threshold)
+            )
 
-            dflt_excl = numpy.logical_or(dflt_excl, dflt_excl2)
+            dflt_excl = np.logical_or(dflt_excl, dflt_excl2)
 
         else:
-            dflt_excl[~numpy.isnan(tth)] = \
-                ~((tth[~numpy.isnan(tth)] >= 0.0) &
-                  (tth[~numpy.isnan(tth)] <= tth_threshold))
+            dflt_excl[~np.isnan(tth)] = ~(
+                (tth[~np.isnan(tth)] >= 0.0)
+                & (tth[~np.isnan(tth)] <= tth_threshold)
+            )
             dflt_excl[0] = False
 
         self._pData.exclusions = dflt_excl
@@ -370,10 +410,9 @@ def update_structure_factor(self):
         sf = self.unitcell.CalcXRSF(hkls)
         self.planeData.set_structFact(sf)
 
-    def compute_powder_overlay(self,
-                               ttharray=numpy.linspace(0, 80, 2000),
-                               fwhm=0.25,
-                               scale=1.0):
+    def compute_powder_overlay(
+        self, ttharray=np.linspace(0, 80, 2000), fwhm=0.25, scale=1.0
+    ):
         """
         this function computes a simulated spectra
         for using in place of lines for the powder
@@ -383,12 +422,12 @@ def compute_powder_overlay(self,
 
         requested feature from Amy Jenei
         """
-        tth = numpy.degrees(self.planeData.getTTh())  # convert to degrees
+        tth = np.degrees(self.planeData.getTTh())  # convert to degrees
         Ip = self.planeData.powder_intensity
-        self.powder_overlay = numpy.zeros_like(ttharray)
+        self.powder_overlay = np.zeros_like(ttharray)
         for t, I in zip(tth, Ip):
             p = [t, fwhm]
-            self.powder_overlay += scale*I*_unit_gaussian(p, ttharray)
+            self.powder_overlay += scale * I * _unit_gaussian(p, ttharray)
 
     def remove_duplicate_atoms(self):
         """
@@ -402,8 +441,130 @@ def remove_duplicate_atoms(self):
         self.charge = self.unitcell.chargestates
         self._hkls_changed()
 
-    def _readCif(self, fcif=DFLT_NAME+'.cif'):
+    def vt(self, temperature=None):
+        '''calculate volume at high
+        temperature
+        '''
+        alpha0 = self.thermal_expansion
+        alpha1 = self.thermal_expansion_dt
+        if temperature is None:
+            vt = self.v0
+        else:
+            delT = temperature - 298
+            delT2 = temperature**2 - 298**2
+            vt = self.v0 * np.exp(alpha0 * delT + 0.5 * alpha1 * delT2)
+        return vt
+
+    def kt(self, temperature=None):
+        '''calculate bulk modulus for
+        high temperature
+        '''
+        k0 = self.k0
+        if temperature is None:
+            return k0
+        else:
+            delT = temperature - 298
+            return k0 + self.dk0dt * delT
+
+    def ktp(self, temperature=None):
+        '''calculate bulk modulus derivative
+        for high temperature
+        '''
+        k0p = self.k0p
+        if temperature is None:
+            return k0p
+        else:
+            delT = temperature - 298
+            return k0p + self.dk0pdt * delT
+
+    @property
+    def pt_lp_factor(self):
+        return (self.unitcell.vol * 1e3 / self.v0) ** (1 / 3)
+
+    @property
+    def lparms0(self):
+        # Get the lattice parameters for 0 pressure and temperature (at v0)
+        lparms = self.lparms
+        return np.array([
+            *(lparms[:3] / self.pt_lp_factor),
+            *lparms[3:],
+        ])
+
+    def calc_pressure(self, volume=None, temperature=None):
+        '''calculate the pressure given the volume
+        and temperature using the third order
+        birch-murnaghan equation of state.
+        '''
+        if volume is None:
+            return 0
+        else:
+            vt = self.vt(temperature=temperature)
+            kt = self.kt(temperature=temperature)
+            ktp = self.ktp(temperature=temperature)
+            f = 0.5 * ((vt / volume) ** (2.0 / 3.0) - 1)
+
+            return (
+                3.0 * kt * f * (1 - 1.5 * (4 - ktp) * f) * (1 + 2 * f) ** 2.5
+            )
+
+    def calc_volume(self, pressure=None, temperature=None):
+        '''solve for volume in the birch-murnaghan EoS to
+        compute the volume. this number will be propagated
+        to the Material object as updated lattice constants.
+        '''
+        vt = self.vt(temperature=temperature)
+        kt = self.kt(temperature=temperature)
+        ktp = self.ktp(temperature=temperature)
+
+        if pressure is None:
+            return vt
+        else:
+            alpha = 0.75 * (ktp - 4)
+            p = np.zeros(
+                [
+                    10,
+                ]
+            )
+            p[0] = 1.0
+            p[2] = (1 - 2 * alpha) / alpha
+            p[4] = (alpha - 1) / alpha
+            p[9] = -2 * pressure / 3 / kt / alpha
+            res = np.roots(p)
+            res = res[np.isreal(res)]
+            res = 1 / np.real(res) ** 3
+
+            mask = np.logical_and(res >= 0.0, res <= 1.0)
+            res = res[mask]
+            if len(res) != 1:
+                msg = 'more than one physically ' 'reasonable solution found!'
+                raise ValueError(msg)
+            return res[0] * vt
 
+    def calc_lp_factor(self, pressure=None, temperature=None):
+        '''calculate the factor to multiply the lattice
+        constants by. only the lengths will be modified, the
+        angles will be kept constant.
+        '''
+        vt = self.vt(temperature=temperature)
+        vpt = self.calc_volume(pressure=pressure, temperature=temperature)
+        return (vpt / vt) ** (1.0 / 3.0)
+
+    def calc_lp_at_PT(self, pressure=None, temperature=None):
+        '''calculate the lattice parameters for a given
+        pressure and temperature using the BM EoS. This
+        is the main function which will be called from
+        the GUI.
+        '''
+        f = self.calc_lp_factor(pressure=pressure, temperature=temperature)
+        lparms0 = self.lparms0
+        return np.array(
+            [
+                *(f * lparms0[:3]),
+                *lparms0[3:],
+            ]
+        )
+
+    def _readCif(self, fcif=DFLT_NAME + '.cif'):
         """
         >> @AUTHOR:     Saransh Singh, Lawrence Livermore National Lab,
                         saransh1@llnl.gov
@@ -422,7 +583,7 @@ def _readCif(self, fcif=DFLT_NAME+'.cif'):
 
         # read the file
         for k in cif.keys():
-            if('_cell_length_a' in cif[k]):
+            if '_cell_length_a' in cif[k]:
                 m = k
                 break
         cifdata = cif[m]
@@ -430,51 +591,58 @@ def _readCif(self, fcif=DFLT_NAME+'.cif'):
 
         # make sure the space group is present in the cif file, either as
         # international table number, hermann-maguain or hall symbol
-        sgkey = ['_space_group_IT_number',
-                 '_symmetry_space_group_name_h-m',
-                 '_symmetry_space_group_name_hall',
-                 '_symmetry_Int_Tables_number']
+        sgkey = [
+            '_space_group_IT_number',
+            '_symmetry_space_group_name_h-m',
+            '_symmetry_space_group_name_hall',
+            '_symmetry_Int_Tables_number',
+        ]
 
         sgdata = False
         for key in sgkey:
             sgdata = sgdata or (key in cifdata)
-            if(sgdata):
+            if sgdata:
                 skey = key
                 break
 
-        if(not(sgdata)):
+        if not (sgdata):
             raise RuntimeError(' No space group information in CIF file! ')
 
         sgnum = 0
         if skey is sgkey[0]:
             sgnum = int(cifdata[sgkey[0]])
-        elif (skey is sgkey[1]):
+        elif skey is sgkey[1]:
             HM = cifdata[sgkey[1]]
             HM = HM.replace(" ", "")
             HM = HM.replace("_", "")
             sgnum = HM_to_sgnum[HM]
-        elif (skey is sgkey[2]):
+        elif skey is sgkey[2]:
             hall = cifdata[sgkey[2]]
             hall = hall.replace(" ", "")
             sgnum = Hall_to_sgnum[HM]
-        elif(skey is sgkey[3]):
+        elif skey is sgkey[3]:
             sgnum = int(cifdata[sgkey[3]])
 
         # lattice parameters
         lparms = []
-        lpkey = ['_cell_length_a', '_cell_length_b',
-                 '_cell_length_c', '_cell_angle_alpha',
-                 '_cell_angle_beta', '_cell_angle_gamma']
+        lpkey = [
+            '_cell_length_a',
+            '_cell_length_b',
+            '_cell_length_c',
+            '_cell_angle_alpha',
+            '_cell_angle_beta',
+            '_cell_angle_gamma',
+        ]
 
         for key in lpkey:
             n = cifdata[key].find('(')
-            if(n != -1):
+            if n != -1:
                 lparms.append(float(cifdata[key][:n]))
             else:
                 lparms.append(float(cifdata[key]))
 
         for i in range(6):
-            if(i < 3):
+            if i < 3:
                 lparms[i] = _angstroms(lparms[i])
             else:
                 lparms[i] = _degrees(lparms[i])
@@ -483,20 +651,27 @@ def _readCif(self, fcif=DFLT_NAME+'.cif'):
         self.sgnum = sgnum
 
         # fractional atomic site, occ and vibration amplitude
-        fracsitekey = ['_atom_site_fract_x', '_atom_site_fract_y',
-                       '_atom_site_fract_z', ]
-
-        occ_U = ['_atom_site_occupancy',
-                 '_atom_site_u_iso_or_equiv', '_atom_site_U_iso_or_equiv']
+        fracsitekey = [
+            '_atom_site_fract_x',
+            '_atom_site_fract_y',
+            '_atom_site_fract_z',
+        ]
+
+        occ_U = [
+            '_atom_site_occupancy',
+            '_atom_site_u_iso_or_equiv',
+            '_atom_site_U_iso_or_equiv',
+        ]
 
         sitedata = True
         for key in fracsitekey:
             sitedata = sitedata and (key in cifdata)
 
-        if(not(sitedata)):
+        if not (sitedata):
             raise RuntimeError(
                 ' fractional site position is not present \
-                or incomplete in the CIF file! ')
+                or incomplete in the CIF file! '
+            )
 
         atompos = []
         for key in fracsitekey:
@@ -506,7 +681,7 @@ def _readCif(self, fcif=DFLT_NAME+'.cif'):
             for p in slist:
                 n = p.find('(')
 
-                if(n != -1):
+                if n != -1:
                     pos.append(p[:n])
                 else:
                     pos.append(p)
@@ -515,8 +690,8 @@ def _readCif(self, fcif=DFLT_NAME+'.cif'):
             sometimes cif files have negative values so need to
             bring them back to fractional coordinates between 0-1
             '''
-            pos = numpy.asarray(pos).astype(numpy.float64)
-            pos, _ = numpy.modf(pos+100.0)
+            pos = np.asarray(pos).astype(np.float64)
+            pos, _ = np.modf(pos + 100.0)
             atompos.append(pos)
 
         """note that the vibration amplitude, U is just the amplitude (in A)
@@ -526,12 +701,12 @@ def _readCif(self, fcif=DFLT_NAME+'.cif'):
             this will be done here so we dont have to worry about it later
         """
 
-        pocc = (occ_U[0] in cifdata.keys())
+        pocc = occ_U[0] in cifdata.keys()
         pU = (occ_U[1] in cifdata.keys()) or (occ_U[2] in cifdata.keys())
 
-        if(not pocc):
+        if not pocc:
             warn('occupation fraction not present. setting it to 1')
-            occ = numpy.ones(atompos[0].shape)
+            occ = np.ones(atompos[0].shape)
             atompos.append(occ)
         else:
             slist = cifdata[occ_U[0]]
@@ -539,25 +714,27 @@ def _readCif(self, fcif=DFLT_NAME+'.cif'):
             for p in slist:
                 n = p.find('(')
 
-                if(n != -1):
+                if n != -1:
                     occ.append(p[:n])
                 else:
                     occ.append(p)
 
-            chkstr = numpy.asarray([isinstance(x, str) for x in occ])
-            occstr = numpy.array(occ)
+            chkstr = np.asarray([isinstance(x, str) for x in occ])
+            occstr = np.array(occ)
             occstr[chkstr] = 1.0
 
-            atompos.append(numpy.asarray(occstr).astype(numpy.float64))
+            atompos.append(np.asarray(occstr).astype(np.float64))
 
-        if(not pU):
-            msg = (f"'Debye-Waller factors not present. "
-                   f"setting to same values for all atoms.'")
+        if not pU:
+            msg = (
+                "'Debye-Waller factors not present. "
+                "setting to same values for all atoms.'"
+            )
             warn(msg)
-            U = self.DFLT_U[0] * numpy.ones(atompos[0].shape)
+            U = self.DFLT_U[0] * np.ones(atompos[0].shape)
             self._U = U
         else:
-            if(occ_U[1] in cifdata.keys()):
+            if occ_U[1] in cifdata.keys():
                 k = occ_U[1]
             else:
                 k = occ_U[2]
@@ -567,12 +744,12 @@ def _readCif(self, fcif=DFLT_NAME+'.cif'):
             for p in slist:
                 n = p.find('(')
 
-                if(n != -1):
+                if n != -1:
                     U.append(p[:n])
                 else:
                     U.append(p)
 
-            chkstr = numpy.asarray([isinstance(x, str) for x in U])
+            chkstr = np.asarray([isinstance(x, str) for x in U])
 
             for ii, x in enumerate(chkstr):
                 if x:
@@ -581,18 +758,18 @@ def _readCif(self, fcif=DFLT_NAME+'.cif'):
                     except ValueError:
                         U[ii] = self.DFLT_U[0]
 
-            self._U = numpy.asarray(U)
+            self._U = np.asarray(U)
         '''
         format everything in the right shape etc.
         '''
-        self._atominfo = numpy.asarray(atompos).T
+        self._atominfo = np.asarray(atompos).T
 
         '''
         get atome types here i.e. the atomic number of atoms at each site
         '''
         atype = '_atom_site_type_symbol'
-        patype = (atype in cifdata)
-        if(not patype):
+        patype = atype in cifdata
+        if not patype:
             raise RuntimeError('atom types not defined in cif file.')
 
         satype = cifdata[atype]
@@ -616,7 +793,7 @@ def _readCif(self, fcif=DFLT_NAME+'.cif'):
             atomtype.append(ptable[ss])
             charge.append(c)
 
-        self._atomtype = numpy.asarray(atomtype).astype(numpy.int32)
+        self._atomtype = np.asarray(atomtype).astype(np.int32)
         self._charge = charge
         self._sgsetting = 0
 
@@ -624,6 +801,19 @@ def _readCif(self, fcif=DFLT_NAME+'.cif'):
         self._beamEnergy = Material.DFLT_KEV
         self._tThWidth = Material.DFLT_TTH
 
+        '''set the Birch-Murnaghan equation of state
+        parameters to default values. These values can
+        be updated by user or by reading a JCPDS file
+        '''
+        self.pressure = 0
+        self.temperature = 298
+        self.k0 = 100.0
+        self.k0p = 0.0
+        self.dk0dt = 0.0
+        self.dk0pdt = 0.0
+        self.alpha_t = 0.0
+        self.dalpha_t_dt = 0.0
+
     def _readHDFxtal(self, fhdf=DFLT_NAME, xtal=DFLT_NAME):
         """
         >> @AUTHOR:     Saransh Singh, Lawrence Livermore National Lab,
@@ -651,12 +841,11 @@ def _readHDFxtal(self, fhdf=DFLT_NAME, xtal=DFLT_NAME):
             raise Exception(f'Unknown type for fhdf: {fhdf}')
 
         if xtal not in root_gid:
-            raise IOError('crystal doesn''t exist in material file.')
+            raise IOError('crystal doesn' 't exist in material file.')
 
         gid = root_gid.get(xtal)
 
-        sgnum = numpy.array(gid.get('SpaceGroupNumber'),
-                            dtype=numpy.int32).item()
+        sgnum = np.array(gid.get('SpaceGroupNumber'), dtype=np.int32).item()
         """
             IMPORTANT NOTE:
             note that the latice parameters is nm by default
@@ -667,8 +856,8 @@ def _readHDFxtal(self, fhdf=DFLT_NAME, xtal=DFLT_NAME):
         lparms = list(gid.get('LatticeParameters'))
 
         for i in range(6):
-            if(i < 3):
-                lparms[i] = _angstroms(lparms[i]*10.0)
+            if i < 3:
+                lparms[i] = _angstroms(lparms[i] * 10.0)
             else:
                 lparms[i] = _degrees(lparms[i])
 
@@ -678,61 +867,123 @@ def _readHDFxtal(self, fhdf=DFLT_NAME, xtal=DFLT_NAME):
         self.sgnum = sgnum
 
         # the U factors are related to B by the relation B = 8pi^2 U
-        self._atominfo = numpy.transpose(numpy.array(
-            gid.get('AtomData'), dtype=numpy.float64))
-        self._U = numpy.transpose(numpy.array(
-            gid.get('U'), dtype=numpy.float64))
+        self._atominfo = np.transpose(
+            np.array(gid.get('AtomData'), dtype=np.float64)
+        )
+        self._U = np.transpose(np.array(gid.get('U'), dtype=np.float64))
 
         # read atom types (by atomic number, Z)
-        self._atomtype = numpy.array(gid.get('Atomtypes'), dtype=numpy.int32)
+        self._atomtype = np.array(gid.get('Atomtypes'), dtype=np.int32)
         if 'ChargeStates' in gid:
             self._charge = h5py_read_string(gid['ChargeStates'])
         else:
-            self._charge = ['0']*self._atomtype.shape[0]
+            self._charge = ['0'] * self._atomtype.shape[0]
         self._atom_ntype = self._atomtype.shape[0]
 
-        self._sgsetting = numpy.array(gid.get('SpaceGroupSetting'),
-                                      dtype=numpy.int32).item()
+        self._sgsetting = np.array(
+            gid.get('SpaceGroupSetting'), dtype=np.int32
+        ).item()
 
-        if('stiffness' in gid):
+        if 'stiffness' in gid:
             # we're assuming the stiffness is in units of GPa
-            self.stiffness = numpy.array(gid.get('stiffness'))
-        elif('compliance' in gid):
+            self.stiffness = np.array(gid.get('stiffness'))
+        elif 'compliance' in gid:
             # we're assuming the compliance is in units of TPa^-1
-            self.stiffness = numpy.linalg.inv(
-                numpy.array(gid.get('compliance'))) * 1.e3
+            self.stiffness = (
+                np.linalg.inv(np.array(gid.get('compliance'))) * 1.0e3
+            )
         else:
-            self.stiffness = numpy.zeros([6, 6])
+            self.stiffness = np.zeros([6, 6])
+
+        '''start reading the Birch-Murnaghan equation of state
+        parameters
+        '''
+        self.pressure = 0
+        if 'pressure' in gid:
+            self.pressure = np.array(gid.get('pressure'),
+                                     dtype=np.float64).item()
+
+        self.temperature = 298
+        if 'temperature' in gid:
+            self.temperature = np.array(gid.get('temperature'),
+                                        dtype=np.float64).item()
+
+        self.k0 = 100.0
+        if 'k0' in gid:
+            # this is the isotropic bulk modulus
+            k0 = np.array(gid.get('k0'), dtype=np.float64).item()
+            self.k0 = k0
+
+        self.k0p = 0.0
+        if 'k0p' in gid:
+            # this is the pressure derivation of
+            # the isotropic bulk modulus
+            k0p = np.array(gid.get('k0p'), dtype=np.float64).item()
+            self.k0p = k0p
+
+        self.dk0dt = 0.0
+        if 'dk0dt' in gid:
+            # this is the temperature derivation of
+            # the isotropic bulk modulus
+            dk0dt = np.array(gid.get('dk0dt'), dtype=np.float64).item()
+            self.dk0dt = dk0dt
+
+        self.dk0pdt = 0.0
+        if 'dk0pdt' in gid:
+            # this is the temperature derivation of
+            # the pressure derivative of isotropic bulk modulus
+            dk0pdt = np.array(gid.get('dk0pdt'), dtype=np.float64).item()
+            self.dk0pdt = dk0pdt
+
+        self.alpha_t = 0.0
+        if 'alpha_t' in gid:
+            # this is the temperature derivation of
+            # the pressure derivative of isotropic bulk modulus
+            alpha_t = np.array(gid.get('alpha_t'), dtype=np.float64).item()
+            self.alpha_t = alpha_t
+
+        self.dalpha_t_dt = 0.0
+        if 'dalpha_t_dt' in gid:
+            # this is the temperature derivation of
+            # the pressure derivative of isotropic bulk modulus
+            dalpha_t_dt = np.array(gid.get('dalpha_t_dt'),
+                                   dtype=np.float64).item()
+            self.dalpha_t_dt = dalpha_t_dt
+
+        '''Finished with the BM EOS
+        '''
+
+        if 'v0' in gid:
+            # this is the isotropic ambient unitcell
+            # volume
+            v0 = np.array(gid.get('v0'), dtype=np.float64).item()
+            self.v0 = v0
 
         # if dmin is present in file:
-        if('dmin' in gid):
+        if 'dmin' in gid:
             # if dmin is present in the HDF5 file, then use that
-            dmin = numpy.array(gid.get('dmin'),
-                               dtype=numpy.float64).item()
-            self._dmin = _angstroms(dmin*10.)
+            dmin = np.array(gid.get('dmin'), dtype=np.float64).item()
+            self._dmin = _angstroms(dmin * 10.0)
         else:
             self._dmin = Material.DFLT_DMIN
 
         # if kev is present in file
-        if('kev' in gid):
-            kev = numpy.array(gid.get('kev'),
-                              dtype=numpy.float64).item()
+        if 'kev' in gid:
+            kev = np.array(gid.get('kev'), dtype=np.float64).item()
             self._beamEnergy = _kev(kev)
         else:
             self._beamEnergy = Material.DFLT_KEV
 
         if 'tThWidth' in gid:
-            tThWidth = numpy.array(gid.get('tThWidth'),
-                                   dtype=numpy.float64).item()
-            tThWidth = numpy.radians(tThWidth)
+            tThWidth = np.array(gid.get('tThWidth'), dtype=np.float64).item()
+            tThWidth = np.radians(tThWidth)
         else:
             tThWidth = Material.DFLT_TTH
 
         self._tThWidth = tThWidth
 
-        if('hkls' in gid):
-            self.hkl_from_file = numpy.array(gid.get('hkls'),
-                                             dtype=numpy.int32)
+        if 'hkls' in gid:
+            self.hkl_from_file = np.array(gid.get('hkls'), dtype=np.int32)
 
         if isinstance(fhdf, (Path, str)):
             # Close the file...
@@ -758,18 +1009,51 @@ def dump_material(self, file, path=None):
         AtomInfo['dmin'] = self.unitcell.dmin
         AtomInfo['kev'] = self.beamEnergy.getVal("keV")
         if self.planeData.tThWidth is None:
-            AtomInfo['tThWidth'] = numpy.degrees(Material.DFLT_TTH)
+            AtomInfo['tThWidth'] = np.degrees(Material.DFLT_TTH)
         else:
-            AtomInfo['tThWidth'] = numpy.degrees(self.planeData.tThWidth)
+            AtomInfo['tThWidth'] = np.degrees(self.planeData.tThWidth)
+
+        AtomInfo['pressure'] = self.pressure
+        AtomInfo['temperature'] = self.temperature
+
+        AtomInfo['k0'] = 100.0
+        if hasattr(self, 'k0'):
+            AtomInfo['k0'] = self.k0
+
+        AtomInfo['k0p'] = 0.0
+        if hasattr(self, 'k0p'):
+            AtomInfo['k0p'] = self.k0p
+
+        AtomInfo['v0'] = self.vol
+        if hasattr(self, 'v0'):
+            AtomInfo['v0'] = self.v0
+
+        AtomInfo['dk0dt'] = 0.0
+        if hasattr(self, 'dk0dt'):
+            AtomInfo['dk0dt'] = self.dk0dt
+
+        AtomInfo['dk0pdt'] = 0.0
+        if hasattr(self, 'dk0pdt'):
+            AtomInfo['dk0pdt'] = self.dk0pdt
+
+        AtomInfo['alpha_t'] = 0.0
+        if hasattr(self, 'alpha_t'):
+            AtomInfo['alpha_t'] = self.alpha_t
+
+        AtomInfo['dalpha_t_dt'] = 0.0
+        if hasattr(self, 'dalpha_t_dt'):
+            AtomInfo['dalpha_t_dt'] = self.dalpha_t_dt
         '''
         lattice parameters
         '''
-        lat_param = {'a': self.unitcell.a,
-                     'b': self.unitcell.b,
-                     'c': self.unitcell.c,
-                     'alpha': self.unitcell.alpha,
-                     'beta': self.unitcell.beta,
-                     'gamma': self.unitcell.gamma}
+        lat_param = {
+            'a': self.unitcell.a,
+            'b': self.unitcell.b,
+            'c': self.unitcell.c,
+            'alpha': self.unitcell.alpha,
+            'beta': self.unitcell.beta,
+            'gamma': self.unitcell.gamma,
+        }
 
         Write2H5File(AtomInfo, lat_param, path)
 
@@ -796,12 +1080,16 @@ def spaceGroup(self):
 
     @property
     def vol(self):
-        return self.unitcell.vol*1e3
+        return self.unitcell.vol * 1e3
 
     @property
     def lparms(self):
-        return numpy.array([x.getVal("nm") if x.isLength() else
-                            x.getVal("degrees") for x in self._lparms])
+        return np.array(
+            [
+                x.getVal("nm") if x.isLength() else x.getVal("degrees")
+                for x in self._lparms
+            ]
+        )
 
     @property
     def latticeType(self):
@@ -829,7 +1117,7 @@ def U(self):
 
     @U.setter
     def U(self, Uarr):
-        Uarr = numpy.array(Uarr)
+        Uarr = np.array(Uarr)
         self._U = Uarr
 
     # property:  sgnum
@@ -850,15 +1138,12 @@ def _set_sgnum(self, v):
             self._newUnitcell()
             self._hkls_changed()
 
-    sgnum = property(_get_sgnum, _set_sgnum, None,
-                     "Space group number")
+    sgnum = property(_get_sgnum, _set_sgnum, None, "Space group number")
 
     @property
     def pgnum(self):
         return self.unitcell.pgnum
 
-    # property:  beamEnergy
-
     def _get_beamEnergy(self):
         """Get method for beamEnergy"""
         return self._beamEnergy
@@ -871,7 +1156,7 @@ def _set_beamEnergy(self, keV):
           float arguments.  Also can take a valWUnit
           instance
         """
-        if(isinstance(keV, valWUnit)):
+        if isinstance(keV, valWUnit):
             self._beamEnergy = keV
         else:
             self._beamEnergy = valWUnit('kev', 'energy', keV, 'keV')
@@ -881,16 +1166,18 @@ def _set_beamEnergy(self, keV):
         @TODO make voltage valWUnit instance so this that we dont have to
         track this manually inn the future
         '''
-        self.unitcell.voltage = self.beamEnergy.value*1e3
+        self.unitcell.voltage = self.beamEnergy.value * 1e3
         self.planeData.wavelength = keV
         self._hkls_changed()
 
-    beamEnergy = property(_get_beamEnergy, _set_beamEnergy, None,
-                          "Beam energy in keV")
+    beamEnergy = property(
+        _get_beamEnergy, _set_beamEnergy, None, "Beam energy in keV"
+    )
 
     """
     03/11/2021 SS 1.0 original
     """
+
     @property
     def unitcell(self):
         return self._unitcell
@@ -908,7 +1195,7 @@ def latticeParameters(self):
     @latticeParameters.setter
     def latticeParameters(self, v):
         """Set method for latticeParameters"""
-        if(len(v) != 6):
+        if len(v) != 6:
             v = unitcell._rqpDict[self.unitcell.latticeType][1](v)
         lp = [_angstroms(v[i]) for i in range(3)]
         for i in range(3, 6):
@@ -932,8 +1219,7 @@ def _set_name(self, v):
 
         return
 
-    name = property(_get_name, _set_name, None,
-                    "Name of material")
+    name = property(_get_name, _set_name, None, "Name of material")
 
     @property
     def dmin(self):
@@ -961,8 +1247,10 @@ def charge(self, vals):
         first make sure the lengths are correct
         """
         if len(vals) != len(self.atomtype):
-            msg = (f"incorrect size of charge. "
-                   f"must be same size as number of aom types.")
+            msg = (
+                "incorrect size of charge. "
+                "must be same size as number of aom types."
+            )
             raise ValueError(msg)
         """
         now check if charge states are actually allowed
@@ -971,9 +1259,11 @@ def charge(self, vals):
             elem = ptableinverse[z]
             cs = chargestate[elem]
             if not vals[ii] in cs:
-                msg = (f"element {elem} does not allow "
-                       f"charge state of {vals[ii]}. "
-                       f"allowed charge states are {cs}.")
+                msg = (
+                    f"element {elem} does not allow "
+                    f"charge state of {vals[ii]}. "
+                    f"allowed charge states are {cs}."
+                )
                 raise ValueError(msg)
 
         self._charge = vals
@@ -1004,8 +1294,11 @@ def _set_atominfo(self, v):
         self._atominfo = v
 
     atominfo = property(
-        _get_atominfo, _set_atominfo, None,
-        "Information about atomic positions and electron number")
+        _get_atominfo,
+        _set_atominfo,
+        None,
+        "Information about atomic positions and electron number",
+    )
 
     # property: "atomtype"
     def _get_atomtype(self):
@@ -1017,8 +1310,24 @@ def _set_atomtype(self, v):
         self._atomtype = v
 
     atomtype = property(
-        _get_atomtype, _set_atomtype, None,
-        "Information about atomic types")
+        _get_atomtype, _set_atomtype, None, "Information about atomic types"
+    )
+
+    @property
+    def thermal_expansion(self):
+        return self.alpha_t
+
+    @thermal_expansion.setter
+    def thermal_expansion(self, val):
+        self.alpha_t = val
+
+    @property
+    def thermal_expansion_dt(self):
+        return self.dalpha_t_dt
+
+    @thermal_expansion_dt.setter
+    def thermal_expansion_dt(self, val):
+        self.dalpha_t_dt = val
 
     def _set_atomdata(self, atomtype, atominfo, U, charge):
         """
@@ -1039,32 +1348,40 @@ def _set_atomdata(self, atomtype, atominfo, U, charge):
         """
 
         # check for consistency of sizes here
-        atomtype = numpy.array(atomtype)
-        atominfo = numpy.array(atominfo)
-        U = numpy.array(U)
+        atomtype = np.array(atomtype)
+        atominfo = np.array(atominfo)
+        U = np.array(U)
 
         if atomtype.shape[0] != atominfo.shape[0]:
-            msg = (f"inconsistent shapes: number of atoms "
-                   f"types passed = {atomtype.shape[0]} \n"
-                   f" number of atom positions passed = {atominfo.shape[0]}")
+            msg = (
+                f"inconsistent shapes: number of atoms "
+                f"types passed = {atomtype.shape[0]} \n"
+                f" number of atom positions passed = {atominfo.shape[0]}"
+            )
             raise ValueError(msg)
 
         if atomtype.shape[0] != U.shape[0]:
-            msg = (f"inconsistent shapes: number of atoms "
-                   f"types passed = {atomtype.shape[0]} \n"
-                   f" U passed for {U.shape[0]} atoms.")
+            msg = (
+                f"inconsistent shapes: number of atoms "
+                f"types passed = {atomtype.shape[0]} \n"
+                f" U passed for {U.shape[0]} atoms."
+            )
             raise ValueError(msg)
 
         if atominfo.shape[0] != U.shape[0]:
-            msg = (f"inconsistent shapes: number of atom "
-                   f"positions passed = {atominfo.shape[0]} \n"
-                   f"U passed for {U.shape[0]} atoms.")
+            msg = (
+                f"inconsistent shapes: number of atom "
+                f"positions passed = {atominfo.shape[0]} \n"
+                f"U passed for {U.shape[0]} atoms."
+            )
             raise ValueError(msg)
 
         if len(charge) != atomtype.shape[0]:
-            msg = (f"inconsistent shapes: number of atoms "
-                   f"types passed = {atomtype.shape[0]} \n"
-                   f"charge value passed for {len(charge)} atoms.")
+            msg = (
+                f"inconsistent shapes: number of atoms "
+                f"types passed = {atomtype.shape[0]} \n"
+                f"charge value passed for {len(charge)} atoms."
+            )
             raise ValueError(msg)
 
         self.atomtype = atomtype
@@ -1074,6 +1391,7 @@ def _set_atomdata(self, atomtype, atominfo, U, charge):
 
         self._newUnitcell()
         self.update_structure_factor()
+
     #
     #  ========== Methods
     #
@@ -1091,8 +1409,7 @@ def _set_atomdata(self, atomtype, atominfo, U, charge):
 def loadMaterialList(cfgFile):
     """Load a list of materials from a file
 
-    The file uses the config file format.  See ConfigParser module.
-"""
+    The file uses the config file format.  See ConfigParser module."""
     p = Parser()
     p.read(cfgFile)
     #
@@ -1106,8 +1423,9 @@ def loadMaterialList(cfgFile):
     return matList
 
 
-def load_materials_hdf5(f, dmin=None, kev=None,
-                        sgsetting=Material.DFLT_SGSETTING):
+def load_materials_hdf5(
+    f, dmin=None, kev=None, sgsetting=Material.DFLT_SGSETTING
+):
     """Load materials from an HDF5 file
 
     The file uses the HDF5 file format.
@@ -1132,10 +1450,7 @@ def load_materials_hdf5(f, dmin=None, kev=None,
         'kev': kev,
         'sgsetting': sgsetting,
     }
-    return {
-        name: Material(name, **kwargs)
-        for name in names
-    }
+    return {name: Material(name, **kwargs) for name in names}
 
 
 def save_materials_hdf5(f, materials, path=None):
@@ -1147,8 +1462,12 @@ def save_materials_hdf5(f, materials, path=None):
 def hkls_match(a, b):
     # Check if hkls match. Expects inputs to have shape (x, 3).
     def sorted_hkls(x):
-        return x[numpy.lexsort((x[:, 2], x[:, 1], x[:, 0]))]
-    return numpy.array_equal(sorted_hkls(a), sorted_hkls(b))
+        return x[np.lexsort((x[:, 2], x[:, 1], x[:, 0]))]
+
+    if a.shape != b.shape:
+        return False
+
+    return np.array_equal(sorted_hkls(a), sorted_hkls(b))
 
 
 def map_hkls_to_new(old_pdata, new_pdata):
@@ -1162,7 +1481,8 @@ def get_hkl_strings(pdata):
         'ar2': get_hkl_strings(new_pdata),
         'return_indices': True,
     }
-    return numpy.intersect1d(**kwargs)[1:]
+    return np.intersect1d(**kwargs)[1:]
+
 
 #
 #  ============================== Executable section for testing
diff --git a/hexrd/mksupport.py b/hexrd/material/mksupport.py
similarity index 92%
rename from hexrd/mksupport.py
rename to hexrd/material/mksupport.py
index 11d97199e..083272624 100644
--- a/hexrd/mksupport.py
+++ b/hexrd/material/mksupport.py
@@ -1,4 +1,4 @@
-from hexrd.symbols import pstr_Elements, sitesym, \
+from hexrd.material.symbols import pstr_Elements, sitesym, \
     tworig, PrintPossibleSG, TRIG, pstr_spacegroup,\
     pstr_mkxtal
 import h5py
@@ -6,7 +6,7 @@
 import numpy as np
 import datetime
 import getpass
-from hexrd.unitcell import _StiffnessDict, _pgDict
+from hexrd.material.unitcell import _StiffnessDict, _pgDict
 
 def mk(filename, xtalname):
 
@@ -515,8 +515,34 @@ def WriteH5Data(fid, AtomInfo, lat_param, path=None):
     did = gid.create_dataset("stiffness", (6, 6), dtype=np.float64)
     did.write_direct(np.array(AtomInfo['stiffness'], dtype=np.float64))
 
+    did = gid.create_dataset("pressure", (1,), dtype=np.float64)
+    did.write_direct(np.array(AtomInfo['pressure'], dtype=np.float64))
 
-    if "tThWidth" in AtomInfo: 
+    did = gid.create_dataset("temperature", (1,), dtype=np.float64)
+    did.write_direct(np.array(AtomInfo['temperature'], dtype=np.float64))
+
+    did = gid.create_dataset("k0", (1,), dtype=np.float64)
+    did.write_direct(np.array([AtomInfo['k0']], dtype=np.float64))
+
+    did = gid.create_dataset("k0p", (1,), dtype=np.float64)
+    did.write_direct(np.array([AtomInfo['k0p']], dtype=np.float64))
+
+    did = gid.create_dataset("dk0dt", (1,), dtype=np.float64)
+    did.write_direct(np.array([AtomInfo['dk0dt']], dtype=np.float64))
+
+    did = gid.create_dataset("dk0pdt", (1,), dtype=np.float64)
+    did.write_direct(np.array([AtomInfo['dk0pdt']], dtype=np.float64))
+
+    did = gid.create_dataset("alpha_t", (1,), dtype=np.float64)
+    did.write_direct(np.array([AtomInfo['alpha_t']], dtype=np.float64))
+
+    did = gid.create_dataset("dalpha_t_dt", (1,), dtype=np.float64)
+    did.write_direct(np.array([AtomInfo['dalpha_t_dt']], dtype=np.float64))
+
+    did = gid.create_dataset("v0", (1,), dtype=np.float64)
+    did.write_direct(np.array([AtomInfo['v0']], dtype=np.float64))
+
+    if "tThWidth" in AtomInfo:
         did = gid.create_dataset("tThWidth", (1,), dtype=np.float64)
         did.write_direct(np.array([AtomInfo['tThWidth']], dtype=np.float64))
 
diff --git a/hexrd/spacegroup.py b/hexrd/material/spacegroup.py
similarity index 99%
rename from hexrd/spacegroup.py
rename to hexrd/material/spacegroup.py
index c1a559392..3efb7212c 100644
--- a/hexrd/spacegroup.py
+++ b/hexrd/material/spacegroup.py
@@ -72,7 +72,8 @@
 from collections import OrderedDict
 from math import sqrt, floor
 
-from hexrd import symbols, constants, symmetry
+from hexrd import constants
+from hexrd.material import symbols, symmetry
 import numpy as np
 #
 __all__ = ['SpaceGroup']
diff --git a/hexrd/symbols.py b/hexrd/material/symbols.py
similarity index 100%
rename from hexrd/symbols.py
rename to hexrd/material/symbols.py
diff --git a/hexrd/symmetry.py b/hexrd/material/symmetry.py
similarity index 99%
rename from hexrd/symmetry.py
rename to hexrd/material/symmetry.py
index 27e420083..a6a566eb4 100644
--- a/hexrd/symmetry.py
+++ b/hexrd/material/symmetry.py
@@ -37,7 +37,7 @@
 # from hexrd.rotations import quatOfAngleAxis, quatProductMatrix, fixQuat
 from hexrd import rotations as rot
 from hexrd import constants
-from hexrd.utils.decorators import numba_njit_if_available
+from hexrd.utils.decorators import memoize, numba_njit_if_available
 
 
 # =============================================================================
@@ -417,6 +417,8 @@ def SYM_fillgen(t):
     mat = np.broadcast_to(mat, [1,4,4])
     return mat
 
+
+@memoize(maxsize=20)
 def GenerateSGSym(sgnum, setting=0):
 
     '''
diff --git a/hexrd/unitcell.py b/hexrd/material/unitcell.py
similarity index 99%
rename from hexrd/unitcell.py
rename to hexrd/material/unitcell.py
index 094981dc2..2acbb22b3 100644
--- a/hexrd/unitcell.py
+++ b/hexrd/material/unitcell.py
@@ -1,8 +1,7 @@
 import importlib.resources
 import numpy as np
 from hexrd import constants
-from hexrd import symmetry, symbols
-from hexrd.spacegroup import Allowed_HKLs
+from hexrd.material import spacegroup, symbols, symmetry
 from hexrd.ipfcolor import sphere_sector, colorspace
 from hexrd.valunits import valWUnit
 import hexrd.resources
@@ -1033,7 +1032,7 @@ def getHKLs(self, dmin):
                             for ik in np.arange(kmax, kmin, -1)
                             for il in np.arange(lmax, lmin, -1)])
 
-        hkl_allowed = Allowed_HKLs(self.sgnum, hkllist)
+        hkl_allowed = spacegroup.Allowed_HKLs(self.sgnum, hkllist)
 
         hkl = []
         dsp = []
diff --git a/hexrd/rotations.py b/hexrd/rotations.py
index c419cc29c..464e9a0af 100644
--- a/hexrd/rotations.py
+++ b/hexrd/rotations.py
@@ -49,7 +49,7 @@
     columnNorm, unitVector, \
     skewMatrixOfVector, findDuplicateVectors, \
     multMatArray, nullSpace
-from hexrd import symmetry
+
 from hexrd.utils.decorators import numba_njit_if_available
 
 # =============================================================================
@@ -445,7 +445,7 @@ def quatAverageCluster(q_in, qsym):
             q_in)
 
         # second, re-cast to FR
-        qrot = symmetry.toFundamentalRegion(qrot.squeeze(), crysSym=qsym)
+        qrot = toFundamentalRegion(qrot.squeeze(), crysSym=qsym)
 
         # compute arithmetic average
         q_bar = unitVector(average(qrot, axis=1).reshape(4, 1))
@@ -456,7 +456,7 @@ def quatAverageCluster(q_in, qsym):
             q_bar)
 
         # re-map
-        q_bar = symmetry.toFundamentalRegion(q_bar, crysSym=qsym)
+        q_bar = toFundamentalRegion(q_bar, crysSym=qsym)
     return q_bar
 
 
diff --git a/hexrd/wppf/phase.py b/hexrd/wppf/phase.py
index 8ef597208..b90ca018c 100644
--- a/hexrd/wppf/phase.py
+++ b/hexrd/wppf/phase.py
@@ -1,9 +1,10 @@
 import numpy as np
 from hexrd.valunits import valWUnit
-from hexrd.spacegroup import Allowed_HKLs, SpaceGroup
-from hexrd import symmetry, symbols, constants
+from hexrd.material.spacegroup import Allowed_HKLs, SpaceGroup
+from hexrd import constants
+from hexrd.material import symmetry, symbols
 from hexrd.material import Material
-from hexrd.unitcell import _rqpDict
+from hexrd.material.unitcell import _rqpDict
 from hexrd.wppf import wppfsupport
 from hexrd.wppf.xtal import _calc_dspacing, _get_tth, _calcxrsf,\
 _calc_extinction_factor, _calc_absorption_factor
@@ -21,7 +22,8 @@ class Material_LeBail:
                     09/14/2020 SS 1.1 class can now be initialized using
                     a material.Material class instance
     >> @DETAILS:    Material_LeBail class is a stripped down version of the
-                    materials.Material class.this is done to keep the class lightweight
+                    material.Material class.this is done to keep the class
+                    lightweight
                     and make sure only the information necessary for the lebail fit is kept
 
     =========================================================================================
@@ -605,9 +607,10 @@ class Material_Rietveld:
     >> @DATE:       05/18/2020 SS 1.0 original
                     02/01/2021 SS 1.1 class can now be initialized using a
                     material.Material class instance
-    >> @DETAILS:    Material_LeBail class is a stripped down version of the materials.Material
-                    class.this is done to keep the class lightweight and make sure only the
-                    information necessary for the Rietveld fit is kept
+    >> @DETAILS:    Material_LeBail class is a stripped down version of the
+                    material.Material class. This is done to keep the class
+                    lightweight and make sure only the information necessary
+                    for the Rietveld fit is kept
     ===========================================================================================
      ==========================================================================================
     """
diff --git a/hexrd/wppf/wppfsupport.py b/hexrd/wppf/wppfsupport.py
index 07e320173..0720e56f7 100644
--- a/hexrd/wppf/wppfsupport.py
+++ b/hexrd/wppf/wppfsupport.py
@@ -32,12 +32,12 @@
 classes are put here to minimize code duplication. Some examples
 include initialize background, generate_default_parameter list etc.
 """
-from hexrd.symbols import pstr_spacegroup
+from hexrd.material.symbols import pstr_spacegroup
 from hexrd.wppf.parameters import Parameters
 from lmfit import Parameters as Parameters_lmfit
 from hexrd.wppf.phase import Phases_LeBail, Phases_Rietveld
 from hexrd.material import Material
-from hexrd.unitcell import _rqpDict
+from hexrd.material.unitcell import _rqpDict
 import hexrd
 import numpy as np
 from hexrd import constants
diff --git a/hexrd/xrdutil/utils.py b/hexrd/xrdutil/utils.py
index 376fbf40d..768e05ddc 100644
--- a/hexrd/xrdutil/utils.py
+++ b/hexrd/xrdutil/utils.py
@@ -32,7 +32,7 @@
 from hexrd import matrixutil as mutil
 from hexrd import gridutil as gutil
 
-from hexrd.crystallography import processWavelength, PlaneData
+from hexrd.material.crystallography import processWavelength, PlaneData
 
 from hexrd.transforms import xf
 from hexrd.transforms import xfcapi
diff --git a/pyproject.toml b/pyproject.toml
index f08833fec..eb28bae5d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,2 +1,6 @@
 [build-system]
 requires = ["setuptools", "wheel", "numpy<1.27", "setuptools_scm[toml]", "pybind11>=2.11.0"]
+
+[tool.black]
+line-length = 79
+skip-string-normalization = true
diff --git a/tests/find_orientations_testing.py b/tests/find_orientations_testing.py
index b9ff117fa..d8a0907f3 100755
--- a/tests/find_orientations_testing.py
+++ b/tests/find_orientations_testing.py
@@ -11,7 +11,7 @@
 
 import numpy as np
 
-from hexrd.crystallography import PlaneData
+from hexrd.material.crystallography import PlaneData
 from hexrd.rotations import misorientation
 
 
diff --git a/tests/test_find_orientations.py b/tests/test_find_orientations.py
index bb4a472fe..c6208f012 100644
--- a/tests/test_find_orientations.py
+++ b/tests/test_find_orientations.py
@@ -10,7 +10,7 @@
 
 from hexrd.findorientations import find_orientations, generate_eta_ome_maps
 from hexrd import config
-from hexrd.crystallography import PlaneData
+from hexrd.material.crystallography import PlaneData
 
 import find_orientations_testing as test_utils