diff --git a/fpm.toml b/fpm.toml
index bcd024f..4578174 100644
--- a/fpm.toml
+++ b/fpm.toml
@@ -70,6 +70,11 @@ name = "test_periodictable"
 source-dir="test"
 main = "test_periodictable.f90"
 
+[[test]]
+name = "test_vector_allreps"
+source-dir="test"
+main = "test_vector_allreps.f90"
+
 [[test]]
 name = "check_interactions_symmetries"
 source-dir="test"
diff --git a/src/screening.f90 b/src/screening.f90
index acd099d..be0bf0c 100644
--- a/src/screening.f90
+++ b/src/screening.f90
@@ -280,8 +280,8 @@ subroutine calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)
     !Silicon
     !omega_plasma = 1.0e-9*hbar*sqrt(el%conc_el/perm0/crys%epsilon0/(0.267*me)) !eV
     !wGaN
-    !omega_plasma = 1.0e-9*hbar*sqrt(el%conc_el/perm0/crys%epsilon0/(0.2*me)) !eV
-    omega_plasma = 1.0e-9*hbar*sqrt(el%conc_el/perm0/crys%epsiloninf/(0.2*me)) !eV
+    omega_plasma = 1.0e-9*hbar*sqrt(el%conc_el/perm0/crys%epsilon0/(0.2*me)) !eV
+    !omega_plasma = 1.0e-9*hbar*sqrt(el%conc_el/perm0/crys%epsiloninf/(0.2*me)) !eV
     
     if(this_image() == 1) then
        print*, "plasmon energy = ", omega_plasma
@@ -356,7 +356,7 @@ subroutine calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)
           call calculate_Hilbert_weights(&
                w_disc = energylist, &
                w_cont = energylist, &
-               zeroplus = 1.0e-6_r64, &
+               zeroplus = 1.0e-6_r64, & !Can this magic "small" number be removed?
                Hilbert_weights = Hilbert_weights)
 
           !call Re_head_polarizability_3d_T(Reeps, energylist, Imeps, Hilbert_weights)
@@ -376,11 +376,11 @@ subroutine calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)
 !!$               (real(eps) + oneI*pi*spec_eps)/perm0*qe*1.0e9_r64
 
        if(all(qcrys == 0.0_r64)) then
-          diel(iq, 1) = 0.0_r64
-          !DEBUG: This is not the correct limit. Just leaving it here
-          !to get the plasmon peak in the loss function.
-          diel(iq, 2:numomega) = 1.0_r64 - &
-               (omega_plasma/energylist(2:numomega))**2 + oneI*1.0e-9
+          !DEBUG: This is not a generally computable limit since the plasmon
+          !energy expression requires a single effective mass.
+          !Just leaving it here for now to get the plasmon peak in the loss function.
+          diel(iq, 1:numomega) = 1.0_r64 - &
+               (omega_plasma/energylist(1:numomega))**2 + oneI*1.0e-9
        else
           !DBG
           diel(iq, :) = 1.0_r64 - &
@@ -392,8 +392,7 @@ subroutine calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)
     call co_sum(diel)
 
     !Handle Omega = 0 case
-    !diel(2:numq, 1) = crys%epsilon0*&
-    !     (1.0_r64 + (crys%qTF/qmaglist(2:numq))**2)
+    diel(:, 1) = 1.0_r64 + (crys%qTF/qmaglist(:))**2
     
     !Print to file
     call write2file_rank2_real("RPA_dielectric_3D_G0_qpath", qlist)
diff --git a/src/vector.f90 b/src/vector.f90
new file mode 100644
index 0000000..fa3dae0
--- /dev/null
+++ b/src/vector.f90
@@ -0,0 +1,90 @@
+module vector_allreps_module
+
+  use precision, only: i64, r64
+  use misc, only: mux_vector, demux_vector
+  
+  implicit none
+
+  private
+  public :: vector_allreps, &
+       vector_allreps_add, vector_allreps_sub, &
+       vector_allreps_print
+  
+  type vector_allreps
+     !! A container for a (3-)vector and relevant arithmetic operations.
+     !! It is convenient to use under certain circumstances, for example,
+     !! when low level vector arithmetic is repetitive and error prone.
+     !! However, the low level method will be more performant.
+     
+     integer(i64) :: muxed_index = -1
+     integer(i64) :: int(3) = 0
+     real(r64) :: frac(3) = 0.0
+     real(r64) :: cart(3) = 0.0
+  end type vector_allreps
+
+  interface vector_allreps
+     module procedure create
+  end interface vector_allreps
+
+contains
+  
+  function create(ind, grid, primitive_vecs) result(vector_obj)
+    integer(i64), intent(in) :: ind, grid(3)
+    real(r64), intent(in) :: primitive_vecs(3, 3)
+    type(vector_allreps) :: vector_obj
+
+    !vector_obj%grid = grid
+
+    !vector_obj%primitive_vecs = primitive_vecs
+    
+    vector_obj%muxed_index = ind
+
+    call demux_vector(ind, vector_obj%int, &
+         grid, 0_i64)
+    
+    vector_obj%frac = dble(vector_obj%int)/grid
+    
+    vector_obj%cart = matmul(primitive_vecs, vector_obj%frac)
+  end function create
+
+  subroutine vector_allreps_print(v)
+    !! Printer
+    
+    type(vector_allreps), intent(in) :: v
+
+    print*, v%muxed_index
+    print*, v%int
+    print*, v%frac
+    print*, v%cart
+  end subroutine vector_allreps_print
+  
+  pure function vector_allreps_add(v1, v2, grid, primitive_vecs) result(v3)
+    !! Adder
+    
+    type(vector_allreps), intent(in) :: v1
+    type(vector_allreps), intent(in) :: v2
+    integer(i64), intent(in) :: grid(3)
+    real(r64), intent(in) :: primitive_vecs(3, 3)
+    type(vector_allreps) :: v3
+    
+    v3%int = modulo(v1%int + v2%int, grid)
+    v3%muxed_index = mux_vector(v3%int, grid, 0_i64)
+    v3%frac = dble(v3%int)/grid
+    v3%cart = matmul(primitive_vecs, v3%frac)
+  end function vector_allreps_add
+
+  pure function vector_allreps_sub(v1, v2, grid, primitive_vecs) result(v3)
+    !! Subtracter
+    
+    type(vector_allreps), intent(in) :: v1
+    type(vector_allreps), intent(in) :: v2
+    integer(i64), intent(in) :: grid(3)
+    real(r64), intent(in) :: primitive_vecs(3, 3)
+    type(vector_allreps) :: v3
+
+    v3%int = modulo(v1%int - v2%int, grid)
+    v3%muxed_index = mux_vector(v3%int, grid, 0_i64)
+    v3%frac = dble(v3%int)/grid
+    v3%cart = matmul(primitive_vecs, v3%frac)
+  end function vector_allreps_sub
+end module vector_allreps_module
diff --git a/test/test_vector_allreps.f90 b/test/test_vector_allreps.f90
new file mode 100644
index 0000000..004c9ce
--- /dev/null
+++ b/test/test_vector_allreps.f90
@@ -0,0 +1,145 @@
+program test_vector_allreps
+
+  use precision, only : r64, i64
+  use testify_m, only : testify
+  use vector_allreps_module, only : &
+       vector_allreps, vector_allreps_add, vector_allreps_sub, &
+       vector_allreps_print
+
+  implicit none
+  
+  integer :: itest
+  integer, parameter :: num_tests = 17
+  type(testify) :: test_array(num_tests), tests_all
+
+  integer(i64) :: imuxed, grid(3)
+  real(r64) :: primitive_vecs(3, 3)
+  !integer(i64), parameter :: N = 3
+  type(vector_allreps) :: v0, v1, v2, v3
+
+  print*, '<<module vector_allreps unit tests>>'
+
+  primitive_vecs(:, 1) = [-0.5_r64, 0.0_r64, 0.5_r64]
+  primitive_vecs(:, 2) = [ 0.0_r64, 0.5_r64, 0.5_r64]
+  primitive_vecs(:, 3) = [-0.5_r64, 0.5_r64, 0.0_r64]
+  
+  grid = [4, 4, 4]
+
+  !Null vector
+  v0 = vector_allreps(1_i64, grid, primitive_vecs)
+  
+  itest = 1
+  test_array(itest) = testify("null vector, integer rep")
+  call test_array(itest)%assert(&
+       v0%int, &
+       [0, 0, 0]*1_i64)
+
+  itest = itest + 1
+  test_array(itest) = testify("null vector, fractional rep")
+  call test_array(itest)%assert(&
+       v0%frac, &
+       [0.0, 0.0, 0.0]*1.0_r64)
+
+  itest = itest + 1
+  test_array(itest) = testify("null vector, Cartesian rep")
+  call test_array(itest)%assert(&
+       v0%cart, &
+       [0.0, 0.0, 0.0]*1.0_r64)
+
+  !A non-null vector
+  v1 = vector_allreps(2_i64, grid, primitive_vecs)
+
+  itest = itest + 1
+  test_array(itest) = testify("a non-null vector, integer rep")
+  call test_array(itest)%assert(&
+       v1%int, &
+       [1, 0, 0]*1_i64)
+
+  itest = itest + 1
+  test_array(itest) = testify("a non-null vector, fractional rep")
+  call test_array(itest)%assert(&
+       v1%frac, &
+       [0.25, 0.0, 0.0]*1.0_r64)
+
+  itest = itest + 1
+  test_array(itest) = testify("a non-null vector, Cartesian rep")
+  call test_array(itest)%assert(&
+       v1%cart, &
+       [-0.125, 0.0, 0.125]*1.0_r64)
+
+  !The last vector
+  v2 = vector_allreps(product(grid), grid, primitive_vecs)
+
+  itest = itest + 1
+  test_array(itest) = testify("another non-null vector, integer rep")
+  call test_array(itest)%assert(&
+       v2%int, &
+       [3, 3, 3]*1_i64)
+
+  itest = itest + 1
+  test_array(itest) = testify("another non-null vector, fractional rep")
+  call test_array(itest)%assert(&
+       v2%frac, &
+       [0.75, 0.75, 0.75]*1.0_r64)
+
+  itest = itest + 1
+  test_array(itest) = testify("another non-null vector, Cartesian rep")
+  call test_array(itest)%assert(&
+       v2%cart, &
+       [-0.75, 0.75, 0.75]*1.0_r64)
+
+  !Addition with Umklapp
+  v3 = vector_allreps_add(v1, v2, grid, primitive_vecs)
+
+  itest = itest + 1
+  test_array(itest) = testify("add vectors, muxed index")
+  call test_array(itest)%assert(v3%muxed_index, 61_i64)
+  
+  itest = itest + 1
+  test_array(itest) = testify("add vectors, integer rep")
+  call test_array(itest)%assert(&
+       v3%int, &
+       [0, 3, 3]*1_i64)
+
+  itest = itest + 1
+  test_array(itest) = testify("add vectors, fractional rep")
+  call test_array(itest)%assert(&
+       v3%frac, &
+       [0.0, 3.0, 3.0]/4.0_r64)
+
+  itest = itest + 1
+  test_array(itest) = testify("add vectors, Cartesian rep")
+  call test_array(itest)%assert(&
+       v3%cart, &
+       [-3.0, 6.0, 3.0]/8.0_r64)
+  
+  !Subtraction with Umklapp
+  v3 = vector_allreps_sub(v1, v2, grid, primitive_vecs)
+
+  itest = itest + 1
+  test_array(itest) = testify("subtract vectors, muxed index")
+  call test_array(itest)%assert(v3%muxed_index, 23_i64)
+  
+  itest = itest + 1
+  test_array(itest) = testify("subtract vectors, integer rep")
+  call test_array(itest)%assert(&
+       v3%int, &
+       [2, 1, 1]*1_i64)
+
+  itest = itest + 1
+  test_array(itest) = testify("subtract vectors, fractional rep")
+  call test_array(itest)%assert(&
+       v3%frac, &
+       [0.5, 0.25, 0.25]*1.0_r64)
+
+  itest = itest + 1
+  test_array(itest) = testify("subtract vectors, Cartesian rep")
+  call test_array(itest)%assert(&
+       v3%cart, &
+       [-3.0, 2.0, 3.0]/8.0_r64)
+      
+  tests_all = testify(test_array)
+  call tests_all%report
+
+  if(tests_all%get_status() .eqv. .false.) error stop -1
+end program test_vector_allreps