diff --git a/src/tensor_operations/matrix_algebra.jl b/src/tensor_operations/matrix_algebra.jl
index 4099c783c5..eed71ac72a 100644
--- a/src/tensor_operations/matrix_algebra.jl
+++ b/src/tensor_operations/matrix_algebra.jl
@@ -65,8 +65,8 @@ function exp(A::ITensor, Linds, Rinds; ishermitian=false)
   end
 
   # <fermions>
-  auto_fermion_enabled = using_auto_fermion()
-  if auto_fermion_enabled
+  fermionic_itensor = using_auto_fermion() && has_fermionic_subspaces(inds(A))
+  if fermionic_itensor
     # If fermionic, bring indices into i',j',..,dag(j),dag(i)
     # ordering with Out indices coming before In indices
     # Resulting tensor acts like a normal matrix (no extra signs
@@ -95,7 +95,7 @@ function exp(A::ITensor, Linds, Rinds; ishermitian=false)
   expA = (itensor(expAT) * dag(CR)) * dag(CL)
 
   # <fermions>
-  if auto_fermion_enabled
+  if fermionic_itensor
     # Ensure expA indices in "matrix" form before re-enabling fermion system
     expA = permute(expA, ordered_inds)
     enable_auto_fermion()
diff --git a/test/base/test_fermions.jl b/test/base/test_fermions.jl
index b68d273663..3ed5183798 100644
--- a/test/base/test_fermions.jl
+++ b/test/base/test_fermions.jl
@@ -1,3 +1,4 @@
+@eval module $(gensym())
 using ITensors, Test
 import ITensors: Out, In
 using ITensors.SiteTypes: op, siteind, siteinds
@@ -801,7 +802,32 @@ using ITensors.SiteTypes: op, siteind, siteinds
 
     # Check wrong index ordering fails (i.e. we are actually paying attention to it)
     @test norm(id_tensor - exp(0.0 * id_tensor, (dag(s[1]), dag(s[2])), (s[2]', s[1]'))) > 1
+
+    # Test a different, random tensor
+    T = random_itensor(s[1]', dag(s[1]))
+    T = 1/2*(T+swapprime(dag(T), 0=>1))
+    t = 0.01
+    eT = exp(t*T)
+    eT_taylor = (op("I", s[1])+t*T+t^2*apply(T, T)/2)
+    @test norm(eT - eT_taylor) < 1E-5
+
+    #
+    # Test that bosonic tensor exp works with auto fermion enabled
+    #
+    j1 = Index([QN("Nb", 0)=>2, QN("Nb", 1)=>2])
+    j2 = Index([QN("Nb", 0)=>2, QN("Nb", 1)=>2])
+    id_tensor = op("I", j1) * op("I", j2)
+    @test id_tensor ≈ exp(0.0 * id_tensor)
+
+    T = random_itensor(j1', dag(j1))
+    T = 1/2*(T+swapprime(dag(T), 0=>1))
+    t = 0.01
+    eT = exp(t*T)
+    eT_taylor = (op("I", j1)+t*T+t^2*apply(T, T)/2)
+    @test norm(eT - eT_taylor) < 1E-5
   end
 
   ITensors.disable_auto_fermion()
 end
+
+end # module
diff --git a/test/base/test_qnitensor.jl b/test/base/test_qnitensor.jl
index 1a9c5d84f7..820c8ccabe 100644
--- a/test/base/test_qnitensor.jl
+++ b/test/base/test_qnitensor.jl
@@ -1,6 +1,7 @@
+@eval module $(gensym())
 using ITensors
 using ITensors.NDTensors
-using ITensors.SiteTypes: siteind, siteinds
+using ITensors.SiteTypes: op, siteind, siteinds
 using LinearAlgebra
 using Random
 using Test
@@ -1915,3 +1916,5 @@ Random.seed!(1234)
     @test !ITensors.have_same_qns([QN(0) => 1, QN(0) => 2, QN(("Sz", 2)) => 1])
   end
 end
+
+end # module