diff --git a/VERSION b/VERSION
index ee90284c..90a27f9c 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-1.0.4
+1.0.5
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 5fc4a885..3010e5fc 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -23,7 +23,7 @@
 author = 'Bernardo Ferreira'
 copyright = '2020, Bernardo Ferreira'
 version = '1.0'
-release = '1.0.4'
+release = '1.0.5'
 
 
 # -- General configuration ----------------------------------------------------
diff --git a/src/cratepy/clustering/clusteringdata.py b/src/cratepy/clustering/clusteringdata.py
index d7504ab4..6b1a1d05 100644
--- a/src/cratepy/clustering/clusteringdata.py
+++ b/src/cratepy/clustering/clusteringdata.py
@@ -367,7 +367,7 @@ def def_gradient_from_log_strain(log_strain):
     """
     # Perform spectral decomposition of material logarithmic strain tensor
     log_eigenvalues, log_eigenvectors, _, _ = \
-        top.spectral_decomposition(log_strain)
+        top.spectral_decomposition(log_strain, is_real_if_close=True)
     # Compute deformation gradient eigenvalues
     dg_eigenvalues = np.exp(log_eigenvalues)
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/src/cratepy/ioput/info.py b/src/cratepy/ioput/info.py
index e4fd97c0..5fd2968c 100644
--- a/src/cratepy/ioput/info.py
+++ b/src/cratepy/ioput/info.py
@@ -93,7 +93,7 @@ def displayinfo(code, *args, **kwargs):
     elif code == '0':
         arguments = \
             ['CRATE - Clustering-based Nonlinear Analysis of Materials',
-             'Created by Bernardo P. Ferreira', 'Release 1.0.4 (Oct 2023)'] \
+             'Created by Bernardo P. Ferreira', 'Release 1.0.5 (Oct 2023)'] \
             + 2*[args[0], ] + list(args[1:3])
         info = tuple(arguments)
         template = '\n' + colorama.Fore.WHITE + tilde_line \
diff --git a/src/cratepy/tensor/tensoroperations.py b/src/cratepy/tensor/tensoroperations.py
index 46fc8a63..e270cb2a 100644
--- a/src/cratepy/tensor/tensoroperations.py
+++ b/src/cratepy/tensor/tensoroperations.py
@@ -211,7 +211,7 @@ def get_id_operators(n_dim):
 #
 #                                                        Spectral decomposition
 # =============================================================================
-def spectral_decomposition(x):
+def spectral_decomposition(x, is_real_if_close=False):
     """Perform spectral decomposition of symmetric second-order tensor.
 
     The computational implementation of the spectral decomposition follows the
@@ -229,6 +229,10 @@ def spectral_decomposition(x):
     x : numpy.ndarray (2d)
         Second-order tensor (square array) whose eigenvalues and eigenvectors
         are computed.
+    is_real_if_close : bool, default=False
+        If True, then drop imaginary parts of eigenvalues and eigenvectors if
+        these are close to zero (tolerance with respect to machine epsilon for
+        input type) and convert to real type.
 
     Returns
     -------
@@ -254,6 +258,12 @@ def spectral_decomposition(x):
     # Perform spectral decomposition
     eigenvalues, eigenvectors = np.linalg.eig(x)
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    # If imaginary parts are close to zero (tolerance with respect to machine
+    # epsilon for input type), then drop imaginary part and convert to real
+    if is_real_if_close:
+        eigenvalues = np.real_if_close(eigenvalues)
+        eigenvectors = np.real_if_close(eigenvectors)
+    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     # Get eigenvalues sorted in descending order
     sort_idxs = np.argsort(eigenvalues)[::-1]
     # Sort eigenvalues in descending order and eigenvectors accordingly
@@ -408,8 +418,8 @@ def isotropic_tensor(mode, x):
         raise RuntimeError('Unknown scalar function.')
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     # Perform spectral decomposition
-    eigenvalues, eigenvectors, eig_multiplicity, eigenprojections = \
-        spectral_decomposition(x)
+    eigenvalues, _, _, eigenprojections = \
+        spectral_decomposition(x, is_real_if_close=True)
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     # Initialize isotropic symmetric tensor-valued function
     y = np.zeros(x.shape)
@@ -467,8 +477,8 @@ def derivative_isotropic_tensor(mode, x):
         raise RuntimeError('Unknown scalar function.')
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     # Perform spectral decomposition
-    eigenvalues, eigenvectors, eig_multiplicity, eigenprojections = \
-        spectral_decomposition(x)
+    eigenvalues, _, eig_multiplicity, eigenprojections = \
+        spectral_decomposition(x, is_real_if_close=True)
     # Compute number of distinct eigenvalues
     n_eig_distinct = n_dim - \
         np.sum([1 for key, val in eig_multiplicity.items() if val == 0])