diff --git a/app/runtime_error_detection/Randomized.F90 b/app/runtime_error_detection/Randomized.F90
index 6e4fa84..32e022c 100644
--- a/app/runtime_error_detection/Randomized.F90
+++ b/app/runtime_error_detection/Randomized.F90
@@ -2,13 +2,16 @@ program Randomized
   use iso_fortran_env, only: error_unit
   use MAC, only: container
   use WannInt_kinds, only: wp => dp
+  use WannInt_definitions, only: cmplx_1
   use WannInt, only: crystal
 
   implicit none
 
-  integer, parameter :: minbnd = 20, maxbnd = 300, &
-                        minnr = 500, maxnr = 3000, &
-                        minder = 1, maxder = 10
+  integer, parameter :: minbnd = 20, maxbnd = 200, &
+                        minnr = 500, maxnr = 2000, &
+                        minder = 1, maxder = 6
+
+  real(wp), parameter :: tol = 1.0E2_wp
 
   real(wp) :: rnd, &
               rnbnd, rnr, rnder
@@ -108,6 +111,9 @@ subroutine random_exec_h_b1(nbnd, nr, nder)
     dlb(2, :) = [0.0_wp, 1.0_wp, 0.0_wp]
     dlb(3, :) = [0.0_wp, 0.0_wp, 1.0_wp]
 
+    tunnellings = cmplx_1
+    dipoles = cmplx_1
+
     call dummy%construct(name="dummy", &
                          direct_lattice_basis=dlb, &
                          num_bands=nbnd, &
@@ -118,6 +124,8 @@ subroutine random_exec_h_b1(nbnd, nr, nder)
 
     write (error_unit, fmt="(A)") "All variables allocated."
     alloc_container_res = dummy%hamiltonian(kpt=[0.0_wp, 0.0_wp, 0.0_wp], derivative=nder, all=.true.)
+    if (abs(abs(alloc_container_res(1)%cdp_storage(1))/real(dummy%nrpts()) - 1.0_wp) > tol*epsilon(1.0_wp)) &
+    error stop "Mismatch with reference."
 
   end subroutine random_exec_h_b1
 
@@ -139,6 +147,9 @@ subroutine random_exec_h_b2(nbnd, nr, nder)
     dlb(2, :) = [0.0_wp, 1.0_wp, 0.0_wp]
     dlb(3, :) = [0.0_wp, 0.0_wp, 1.0_wp]
 
+    tunnellings = cmplx_1
+    dipoles = cmplx_1
+
     call dummy%construct(name="dummy", &
                          direct_lattice_basis=dlb, &
                          num_bands=nbnd, &
@@ -149,6 +160,8 @@ subroutine random_exec_h_b2(nbnd, nr, nder)
 
     write (error_unit, fmt="(A)") "All variables allocated."
     contained_res = dummy%hamiltonian(kpt=[0.0_wp, 0.0_wp, 0.0_wp], derivative=nder)
+    if (abs(abs(contained_res%cdp_storage(1))/real(dummy%nrpts())) > tol*epsilon(1.0_wp)) &
+    error stop "Mismatch with reference."
 
   end subroutine random_exec_h_b2
 
@@ -170,6 +183,9 @@ subroutine random_exec_h_b3(nbnd, nr)
     dlb(2, :) = [0.0_wp, 1.0_wp, 0.0_wp]
     dlb(3, :) = [0.0_wp, 0.0_wp, 1.0_wp]
 
+    tunnellings = cmplx_1
+    dipoles = cmplx_1
+
     call dummy%construct(name="dummy", &
                          direct_lattice_basis=dlb, &
                          num_bands=nbnd, &
@@ -180,6 +196,8 @@ subroutine random_exec_h_b3(nbnd, nr)
 
     write (error_unit, fmt="(A)") "All variables allocated."
     resH = dummy%hamiltonian(kpt=[0.0_wp, 0.0_wp, 0.0_wp])
+    if(abs(abs(resH(1, 1))/real(dummy%nrpts()) - 1.0_wp) > tol*epsilon(1.0_wp)) &
+    error stop "Mismatch with reference."
 
   end subroutine random_exec_h_b3
 
@@ -201,6 +219,9 @@ subroutine random_exec_a_b1(nbnd, nr, nder)
     dlb(2, :) = [0.0_wp, 1.0_wp, 0.0_wp]
     dlb(3, :) = [0.0_wp, 0.0_wp, 1.0_wp]
 
+    tunnellings = cmplx_1
+    dipoles = cmplx_1
+
     call dummy%construct(name="dummy", &
                          direct_lattice_basis=dlb, &
                          num_bands=nbnd, &
@@ -211,6 +232,8 @@ subroutine random_exec_a_b1(nbnd, nr, nder)
 
     write (error_unit, fmt="(A)") "All variables allocated."
     alloc_container_res = dummy%berry_connection(kpt=[0.0_wp, 0.0_wp, 0.0_wp], derivative=nder, all=.true.)
+    if(abs(abs(alloc_container_res(1)%cdp_storage(1))/real(dummy%nrpts()) - 1.0_wp) > tol*epsilon(1.0_wp)) &
+    error stop "Mismatch with reference."
 
   end subroutine random_exec_a_b1
 
@@ -232,6 +255,9 @@ subroutine random_exec_a_b2(nbnd, nr, nder)
     dlb(2, :) = [0.0_wp, 1.0_wp, 0.0_wp]
     dlb(3, :) = [0.0_wp, 0.0_wp, 1.0_wp]
 
+    tunnellings = cmplx_1
+    dipoles = cmplx_1
+
     call dummy%construct(name="dummy", &
                          direct_lattice_basis=dlb, &
                          num_bands=nbnd, &
@@ -242,6 +268,8 @@ subroutine random_exec_a_b2(nbnd, nr, nder)
 
     write (error_unit, fmt="(A)") "All variables allocated."
     contained_res = dummy%berry_connection(kpt=[0.0_wp, 0.0_wp, 0.0_wp], derivative=nder)
+    if(abs(abs(contained_res%cdp_storage(1))/real(dummy%nrpts())) > tol*epsilon(1.0_wp)) &
+    error stop "Mismatch with reference."
 
   end subroutine random_exec_a_b2
 
@@ -263,6 +291,9 @@ subroutine random_exec_a_b3(nbnd, nr)
     dlb(2, :) = [0.0_wp, 1.0_wp, 0.0_wp]
     dlb(3, :) = [0.0_wp, 0.0_wp, 1.0_wp]
 
+    tunnellings = cmplx_1
+    dipoles = cmplx_1
+
     call dummy%construct(name="dummy", &
                          direct_lattice_basis=dlb, &
                          num_bands=nbnd, &
@@ -273,6 +304,8 @@ subroutine random_exec_a_b3(nbnd, nr)
 
     write (error_unit, fmt="(A)") "All variables allocated."
     resA = dummy%berry_connection(kpt=[0.0_wp, 0.0_wp, 0.0_wp])
+    if(abs(abs(resA(1, 1, 1))/real(dummy%nrpts()) - 1.0_wp) > tol*epsilon(1.0_wp)) &
+    error stop "Mismatch with reference."
 
   end subroutine random_exec_a_b3
 
diff --git a/src/WannInt.F90 b/src/WannInt.F90
index 69b3711..8ad3620 100644
--- a/src/WannInt.F90
+++ b/src/WannInt.F90
@@ -143,14 +143,14 @@ subroutine construct_from_input(self, name, &
 
     self%num_R_points = size(R_points(:, 1))
 
-    allocate (self%R_points(self%num_R_points, 3), &
+    allocate (self%R_points(3, self%num_R_points), &
               self%deg_R_points(self%num_R_points), stat=istat)
     if (istat /= 0) then
       write (errormsg, "(i20)") istat
       errormsg = "WannInt: Error #5: failure allocating R_points and deg_R_points. stat = "//trim(adjustl(errormsg))//"."
       error stop trim(errormsg)
     endif
-    self%R_points = R_points
+    self%R_points = transpose(R_points)
     self%deg_R_points = 1
 
     allocate (self%real_space_hamiltonian_elements(self%num_R_points, self%bands, self%bands), &
@@ -237,7 +237,7 @@ subroutine construct_from_file(self, name, &
 
     read (unit=stdin, fmt=*) self%num_R_points
 
-    allocate (self%R_points(self%num_R_points, 3), &
+    allocate (self%R_points(3, self%num_R_points), &
               self%deg_R_points(self%num_R_points), stat=istat)
     if (istat /= 0) then
       write (errormsg, "(i20)") istat
@@ -279,7 +279,7 @@ subroutine construct_from_file(self, name, &
     endif
     self%real_space_hamiltonian_elements = cmplx_0
     do irpts = 1, self%num_R_points
-      read (unit=stdin, fmt=*) (self%R_points(irpts, i), i=1, 3)
+      read (unit=stdin, fmt=*) (self%R_points(i, irpts), i=1, 3)
       do i = 1, self%bands
         do j = 1, self%bands
           read (unit=stdin, fmt=*) dummy1, dummy2, dummyR(1), dummyR(2)
@@ -386,7 +386,7 @@ pure function rpts(self)
     integer :: rpts(self%num_R_points, 3)
     if (.not. (self%is_initialized)) error stop &
       "WannInt: Error #6: crystal is not initialized."
-    rpts = self%R_points
+    rpts = transpose(self%R_points)
   end function rpts
 
   pure function deg_rpts(self)
@@ -546,12 +546,12 @@ function hamiltonian_der_driver(self, k, der)
 
         !Compute factor appearing in the exponential.
         !(k is in coords relative to recip. lattice vectors).
-        kdotr = 2.0_wp*pi*dot_product(self%R_points(irpts, :), k)
+        kdotr = 2.0_wp*pi*dot_product(self%R_points(:, irpts), k)
 
         !Compute Bravais lattice vector for label irpts.
         vec = 0.0_wp
         do i = 1, 3
-          vec = vec + self%R_points(irpts, i)*self%direct_l_b(i, :)
+          vec = vec + self%R_points(i, irpts)*self%direct_l_b(i, :)
         enddo
 
         prod_r = 1.0_wp
@@ -694,12 +694,12 @@ function berry_der_driver(self, k, der)
 
         !Compute factor appearing in the exponential.
         !(k is in coords relative to recip. lattice vectors).
-        kdotr = 2.0_wp*pi*dot_product(self%R_points(irpts, :), k)
+        kdotr = 2.0_wp*pi*dot_product(self%R_points(:, irpts), k)
 
         !Compute Bravais lattice vector for label irpts.
         vec = 0.0_wp
         do i = 1, 3
-          vec = vec + self%R_points(irpts, i)*self%direct_l_b(i, :)
+          vec = vec + self%R_points(i, irpts)*self%direct_l_b(i, :)
         enddo
 
         prod_r = 1.0_wp