diff --git a/src/crystal.f90 b/src/crystal.f90
index 8e820b515..f0dc1afcb 100644
--- a/src/crystal.f90
+++ b/src/crystal.f90
@@ -431,44 +431,34 @@ function crystal_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
   real(pREAL), dimension(sum(Ntwin))               :: characteristicShear
 
   integer :: &
-    a, &                                                                                            !< index of active system
     f, &                                                                                            !< index of my family
-    s                                                                                               !< index of my system in current family
-
-  integer, dimension(size(HP_NTWINSYSTEM)), parameter :: &
-    HP_SHEARTWIN = [&
-      1, &  ! <-10.1>{10.2}
-      2, &  ! <11.6>{-1-1.1}
-      3, &  ! <10.-2>{10.1}
-      4  &  ! <11.-3>{11.2}
-      ]                                                                                             !< indicator to formulas below
+    s, e
 
 
-  a = 0
-  myFamilies: do f = 1,size(Ntwin,1)
-    mySystems: do s = 1,Ntwin(f)
-      a = a + 1
-      select case(lattice)
-        case('cF','cI')
-          characteristicShear(a) = 0.5_pREAL*sqrt(2.0_pREAL)
-        case('hP')
-          if (cOverA < 1.0_pREAL .or. cOverA > 2.0_pREAL) &
-            call IO_error(131,ext_msg='crystal_characteristicShear_Twin')
-          select case(HP_SHEARTWIN(f))                                                              ! from Christian & Mahajan 1995 p.29
-            case (1)                                                                                ! <-10.1>{10.2}
-              characteristicShear(a) = (3.0_pREAL-cOverA**2)/sqrt(3.0_pREAL)/CoverA
-            case (2)                                                                                ! <11.6>{-1-1.1}
-              characteristicShear(a) = 1.0_pREAL/cOverA
-            case (3)                                                                                ! <10.-2>{10.1}
-              characteristicShear(a) = (4.0_pREAL*cOverA**2-9.0_pREAL)/sqrt(48.0_pREAL)/cOverA
-            case (4)                                                                                ! <11.-3>{11.2}
-              characteristicShear(a) = 2.0_pREAL*(cOverA**2-2.0_pREAL)/3.0_pREAL/cOverA
-          end select
-        case default
-          call IO_error(137,ext_msg='crystal_characteristicShear_Twin: '//trim(lattice))
-      end select
-    end do mySystems
-  end do myFamilies
+  select case(lattice)
+    case('cF','cI')                                                                                 ! 10.1016/0079-6425(94)00007-7, Table 1
+      characteristicShear = 0.5_pREAL*sqrt(2.0_pREAL)
+    case('hP')                                                                                      ! 10.1016/0079-6425(94)00007-7, Table 3
+      if (cOverA < 1.0_pREAL .or. cOverA > 2.0_pREAL) &
+        call IO_error(131,ext_msg='crystal_characteristicShear_Twin')
+
+      myFamilies: do f = 1,size(Ntwin,1)
+        s = sum(Ntwin(:f-1)) + 1
+        e = sum(Ntwin(:f))
+        select case(f)
+          case (1)                                                                                  ! <-10.1>{10.2}
+            characteristicShear(s:e) = (3.0_pREAL-cOverA**2)/sqrt(3.0_pREAL)/CoverA
+          case (2)                                                                                  ! <11.6>{-1-1.1}
+            characteristicShear(s:e) = 1.0_pREAL/cOverA
+          case (3)                                                                                  ! <10.-2>{10.1}
+            characteristicShear(s:e) = (4.0_pREAL*cOverA**2-9.0_pREAL)/sqrt(48.0_pREAL)/cOverA
+          case (4)                                                                                  ! <11.-3>{11.2}
+            characteristicShear(s:e) = 2.0_pREAL*(cOverA**2-2.0_pREAL)/3.0_pREAL/cOverA
+        end select
+      end do myFamilies
+    case default
+      call IO_error(137,ext_msg='crystal_characteristicShear_Twin: '//trim(lattice))
+  end select
 
 end function crystal_characteristicShear_Twin