diff --git a/include/xatu/Result.hpp b/include/xatu/Result.hpp
index 9161cb7..4a5a11b 100755
--- a/include/xatu/Result.hpp
+++ b/include/xatu/Result.hpp
@@ -53,9 +53,9 @@ class Result {
         void writePhase(int, FILE*);
         void writeExtendedPhase(const arma::cx_vec&, FILE*);
         void writeExtendedPhase(int, FILE*);
-        void writeEigenvalues(FILE*, int n = 0);
-        void writeStates(FILE*, int n = 0);
-        void writeSpin(int, FILE*);
+        void writeEigenvalues(FILE*, int n = 0, double encut = 0.0);
+        void writeStates(FILE*, int n = 0, double encut = 0.0);
+        void writeSpin(int, double, FILE*);
         void writeRealspaceAmplitude(int, int, const arma::rowvec&, FILE*, int);
         virtual void writeRealspaceAmplitude(const arma::cx_vec&, int, const arma::rowvec&, FILE*, int) = 0;
         virtual void writeAbsorptionSpectrum() = 0;
diff --git a/include/xatu/Result.tpp b/include/xatu/Result.tpp
index 09f9939..e4e19e6 100755
--- a/include/xatu/Result.tpp
+++ b/include/xatu/Result.tpp
@@ -401,9 +401,15 @@ void Result<T>::writeRealspaceAmplitude(int stateindex, int holeIndex,
  * @return void 
  */
 template <typename T>
-void Result<T>::writeEigenvalues(FILE* textfile, int n){
+void Result<T>::writeEigenvalues(FILE* textfile, int n, double encut){
 
-    if(n > exciton->excitonbasisdim || n < 0){
+    int newn = n;
+    
+    if (encut != 0.0){
+        newn = (int) arma::abs(eigval - encut).index_min() + 1;
+    }
+    
+    if(newn > exciton->excitonbasisdim || newn < 0){
         throw std::invalid_argument("Optional argument n must be a positive integer equal or below basisdim");
     }
 
@@ -411,7 +417,7 @@ void Result<T>::writeEigenvalues(FILE* textfile, int n){
     fprintf(textfile, "%d\n", exciton->excitonbasisdim);
 
     // second line: how many eigenvalues will follow
-    int maxEigval = (n == 0) ? exciton->excitonbasisdim : n;
+    int maxEigval = (newn == 0) ? exciton->excitonbasisdim : newn;
     fprintf(textfile, "%d\n", maxEigval);
 
     // then one eigenvalue per line
@@ -427,8 +433,14 @@ void Result<T>::writeEigenvalues(FILE* textfile, int n){
  * @return void
  */
 template <typename T>
-void Result<T>::writeStates(FILE* textfile, int n){
-    if(n > exciton->excitonbasisdim || n < 0){
+void Result<T>::writeStates(FILE* textfile, int n, double encut){
+    int newn = n;
+    
+    if (encut != 0.0){
+        newn = (int) arma::abs(eigval - encut).index_min() + 1;
+    }
+    
+    if(newn > exciton->excitonbasisdim || newn < 0){
         throw std::invalid_argument("Optional argument n must be a positive integer equal or below basisdim");
     }
     // First write basis
@@ -442,7 +454,7 @@ void Result<T>::writeStates(FILE* textfile, int n){
                 kpoint(0), kpoint(1), kpoint(2), v, c);
     }
 
-    int nstates = (n == 0) ? exciton->excitonbasisdim : n;  
+    int nstates = (newn == 0) ? exciton->excitonbasisdim : newn;  
     for(unsigned int i = 0; i < nstates; i++){
         for(unsigned int j = 0; j < exciton->excitonbasisdim; j++){
             fprintf(textfile, "%11.7lf\t%11.7lf\t", 
@@ -458,13 +470,19 @@ void Result<T>::writeStates(FILE* textfile, int n){
  * @param textfile Textfile where the spins are written.
  */
 template <typename T>
-void Result<T>::writeSpin(int n, FILE* textfile){
-
-    if(n > exciton->excitonbasisdim || n < 0){
+void Result<T>::writeSpin(int n, double encut, FILE* textfile){
+    
+    int newn = n;
+    
+    if (encut != 0.0){
+        newn = (int) arma::abs(eigval - encut).index_min() + 1;
+    }
+    
+    if(newn > exciton->excitonbasisdim || newn < 0){
         throw std::invalid_argument("Optional argument n must be a positive integer equal or below basisdim");
     }
 
-    int maxState = (n == 0) ? exciton->excitonbasisdim : n;  
+    int maxState = (newn == 0) ? exciton->excitonbasisdim : newn;  
     fprintf(textfile, "n\tSt\tSe\tSh\n");
     for(unsigned int i = 0; i < maxState; i++){
         auto spin = spinX(i);
diff --git a/include/xatu/utils.hpp b/include/xatu/utils.hpp
index 6f69b33..e767146 100755
--- a/include/xatu/utils.hpp
+++ b/include/xatu/utils.hpp
@@ -36,14 +36,20 @@ bool checkIfTriangular(const arma::cx_mat&);
 * @param precision Number of decimals (sets degeneracy threshold)
 **/
 template<typename T>
-void printEnergies(const std::unique_ptr<T>& results, int n = 8, int precision = 6){
+void printEnergies(const std::unique_ptr<T>& results, int n = 8, double encut = 0.0, int precision = 6){
 
     // Print header
     printf("+---------------+-----------------------------+-----------------------------+\n");
     printf("|       N       |          Eigval (eV)        |          Degeneracy         |\n");
     printf("+---------------+-----------------------------+-----------------------------+\n");
-
-    std::vector<std::vector<double>> pairs = detectDegeneracies(results->eigval, n, precision);
+    
+    int newn = n;
+    
+    if (encut != 0.0){
+        newn = (int) arma::abs(results->eigval - encut).index_min() + 1;
+    }
+    
+    std::vector<std::vector<double>> pairs = detectDegeneracies(results->eigval, newn, precision);
     int it = 1;
 
     for(auto pair : pairs){
diff --git a/main/xatu.cpp b/main/xatu.cpp
index f4e43c2..b4f6e2e 100755
--- a/main/xatu.cpp
+++ b/main/xatu.cpp
@@ -22,6 +22,7 @@ int main(int argc, char* argv[]){
     TCLAP::CmdLine cmd("Command line interface options of the Xatu binary. For a more detailed description, refer to the user guide or the API documentation.", ' ', "1.0");
 
     TCLAP::ValueArg<int>    statesArg("n", "states", "Specify number of exciton states to show.", false, 8, "No. states", cmd);
+    TCLAP::ValueArg<double> energycutoff("t", "ecut", "Specify up to high energy to print excitons (eV).", false, 0.0, "En cutoff", cmd);
     TCLAP::ValueArg<int>    precisionArg("p", "precision", "Desired energy precision. Used to compute degeneracies.", false, 6, "No. decimals", cmd);
     TCLAP::SwitchArg        spinArg("s", "spin", "Compute exciton spin and write it to file.", cmd, false);
     TCLAP::ValueArg<int>    dftArg("d", "dft", "Indicates that the system file is a .outp CRYSTAL file.", false, -1, "No. Fock matrices", cmd);
@@ -50,6 +51,7 @@ int main(int argc, char* argv[]){
 
     // Extract information from parsed CLI options
     int nstates        = statesArg.getValue();
+    double encut       = energycutoff.getValue();
     int ncells         = dftArg.getValue();
     int electronNum    = w90Arg.getValue();
     int decimals       = precisionArg.getValue();
@@ -152,7 +154,7 @@ int main(int argc, char* argv[]){
     cout << "|                                    Results                                |" << endl;
     cout << "+---------------------------------------------------------------------------+" << endl;
 
-    xatu::printEnergies(results, nstates, decimals);
+    xatu::printEnergies(results, nstates, encut, decimals);
 
     cout << "+---------------------------------------------------------------------------+" << endl;
     cout << "|                                    Output                                 |" << endl;
@@ -168,7 +170,7 @@ int main(int argc, char* argv[]){
 
         std::cout << "Writing eigvals to file: " << filename_en << std::endl;
         fprintf(textfile_en, "%d\n", excitonConfig->excitonInfo.ncell);
-        results->writeEigenvalues(textfile_en, nstates);
+        results->writeEigenvalues(textfile_en, nstates, encut);
 
         fclose(textfile_en);
     }
@@ -179,7 +181,7 @@ int main(int argc, char* argv[]){
         FILE* textfile_st = fopen(filename_st.c_str(), "w");
 
         std::cout << "Writing states to file: " << filename_st << std::endl;
-        results->writeStates(textfile_st, nstates);
+        results->writeStates(textfile_st, nstates, encut);
 
         fclose(textfile_st);
     }
@@ -231,7 +233,7 @@ int main(int argc, char* argv[]){
         FILE* textfile_spin = fopen(filename_spin.c_str(), "w");
 
         std::cout << "Writing excitons spin fo file: " << filename_spin << std::endl;
-        results->writeSpin(nstates, textfile_spin);
+        results->writeSpin(nstates, encut, textfile_spin);
     }
 
     auto stop = high_resolution_clock::now();
diff --git a/src/ExcitonTB.cpp b/src/ExcitonTB.cpp
index db64500..37b985c 100644
--- a/src/ExcitonTB.cpp
+++ b/src/ExcitonTB.cpp
@@ -420,10 +420,23 @@ double ExcitonTB::keldysh(arma::rowvec r){
 double ExcitonTB::coulomb(arma::rowvec r){
     double cutoff = arma::norm(system->bravaisLattice.row(0)) * cutoff_ + 1E-5;
     double R = abs(arma::norm(r));
-    if (R > cutoff){
-        return 0.0;
+    /*    if (R > cutoff){
+     *       return 0.0;
     }
-    return (R != 0) ? ec/(4E-10*PI*eps0*R) : ec*1E10/(4*PI*eps0*regularization);    
+    return (R != 0) ? ec/(4E-10*PI*eps0*R) : ec*1E10/(4*PI*eps0*regularization);   */ 
+    double eps_bar = (eps_m + eps_s)/2;
+    double potential_value;
+    if(R < 1E-5){
+        potential_value =ec/(4E-10*PI*eps0*eps_bar*regularization);
+    }
+    else if (R > cutoff){
+        potential_value = 0.0;
+    }
+    else{
+        potential_value = ec*1E10/(4*PI*eps0*eps_bar*R);   
+    };
+    
+    return potential_value;
 }
 
 /**
diff --git a/src/Wannier90Configuration.cpp b/src/Wannier90Configuration.cpp
index 648ab4f..f14251a 100644
--- a/src/Wannier90Configuration.cpp
+++ b/src/Wannier90Configuration.cpp
@@ -116,7 +116,7 @@ namespace xatu {
                     iss.clear(); 
                 }
             }
-            // fockMatrices.slice(i) /= Degen(i);  // correcting double-counted neighbors
+            fockMatrices.slice(i) /= Degen(i);  // correcting double-counted neighbors
             std::getline(m_file, line);
         }
 
diff --git a/src/skubo_w.f90 b/src/skubo_w.f90
index c75b1f2..b831163 100644
--- a/src/skubo_w.f90
+++ b/src/skubo_w.f90
@@ -47,6 +47,7 @@ subroutine skubo_w(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,hhop,shop,npointstota
 character(100) file_name_sp_imag
 character(100) file_name_ex_imag
 character(len=:), allocatable :: file_name_strength
+logical :: do_kubo                    ! flag: do Kubo calculation only if input exists
 integer :: p_dot
 
 
@@ -66,104 +67,130 @@ subroutine skubo_w(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,hhop,shop,npointstota
 e_ex=e_ex/27.211385d0
 eigval_stack=eigval_stack/27.211385d0
 
-call crossproduct(R(1,1),R(1,2),R(1,3),R(2,1), &
-R(2,2),R(2,3),cx,cy,cz)
-vcell=sqrt(cx**2+cy**2+cz**2)
+! volume formula depends on the dimensionality of the unit cell 
+! for 1D, we have the magnitude of the vector (we take the magnitude of the three vectors just to be safe)
+! for 2D, we have the crossproduct (as before, but now allowing for yz or xz -oriented lattices)
+! for 3D, we have the triple product dot(z,cross(x,y))
+! this is checked by checking the norm of rkx/rky/rkz to determine how many directions are defined in k
+
+if ( (NORM2(rkx) == 0 .and. NORM2(rky) == 0) .or. &
+      (NORM2(rky) == 0 .and. NORM2(rkz) == 0) .or. &
+      (NORM2(rkx) == 0 .and. NORM2(rkz) == 0) ) then
+  vcell=sqrt(R(1,1)**2+R(1,2)**2+R(1,3)**2+R(2,1)**2+R(2,2)**2+R(2,3)**2+R(3,1)**2+R(3,2)**2+R(3,3)**2)
+else if ( ( NORM2(rkz) == 0 .xor. NORM2(rky) == 0) .or. &
+          (NORM2(rkx) == 0 .xor. NORM2(rky) == 0) ) then
+  if (NORM2(rkx) == 0) then
+    call crossproduct(R(3,1),R(3,2),R(3,3),R(2,1),R(2,2),R(2,3),cx,cy,cz)
+  else if (NORM2(rky) == 0) then
+    call crossproduct(R(1,1),R(1,2),R(1,3),R(3,1),R(3,2),R(3,3),cx,cy,cz)
+  else
+    call crossproduct(R(1,1),R(1,2),R(1,3),R(2,1),R(2,2),R(2,3),cx,cy,cz)
+  endif
+  vcell=sqrt(cx**2+cy**2+cz**2)
+else 
+  call crossproduct(R(1,1),R(1,2),R(1,3),R(2,1),R(2,2),R(2,3),cx,cy,cz)
+  vcell=abs(R(3,1)*cx+R(3,2)*cy+R(3,3)*cz)
+endif
 
 !call fill_nRvec(nR,R,Rvec,nRvec)
 !get printing parameters
 call get_kubo_parameters(w0,wrange,nw,eta,type_broad, &
 file_name_sp,file_name_ex)
-eta=eta/27.211385d0
 
-!SP arrays
-allocate (wp(nw))
-allocate (sigma_w_sp(3,3,nw))
-allocate (hk_ev(norb,nbands))
-allocate (e(nbands))
+!flag indicates whether we actually read parameters
+ do_kubo = (nw .gt. 0)
+ if (.not. do_kubo) then
+   write(*,*) 'Notice: kubo input missing or empty -- conductivity routines will be skipped'
+ endif
+ if (do_kubo) then
+   eta = eta/27.211385d0
+ end if
+
+! allocate the vme array unconditionally – used by the exciton routine
 allocate (vme(npointstotal,3,nbands,nbands))
-wp=0.0d0
-wn_sp=0.0d0
-sigma_w_sp=0.0d0
-do i=1,nw
-  wp(i)=(w0+wrange/dble(nw)*dble(i-1))/27.211385d0
-end do
 
-!exciton arrays
+!exciton arrays (only the part needed for oscillator strengths)
 norb_ex_band=nv_ex*nc_ex !number of electron-hole pairs per k-point
 norb_ex_cut=norb_ex  !total number of optical transitions
 allocate (vme_ex(3,norb_ex_cut,2))
-allocate (sigma_w_ex(3,3,nw))
-allocate (skubo_ex_int(3,3,norb_ex_cut))
 vme_ex=0.0d0
-sigma_w_ex=0.0d0
-skubo_ex_int=0.0d0
 
+!calculate exciton oscillator strengths regardless of Kubo input
 call exciton_oscillator_strength(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,hhop,shop,npointstotal,rkx, &
                                  rky,rkz,fk_ex,e_ex,eigval_stack,eigvec_stack,vme,vme_ex,.false.)
 
+if (do_kubo) then
+  ! -- arrays needed for conductivity evaluation --
+  allocate (wp(nw))
+  allocate (sigma_w_sp(3,3,nw))
+  allocate (hk_ev(norb,nbands))
+  allocate (e(nbands))
+  wp=0.0d0
+  wn_sp=0.0d0
+  sigma_w_sp=0.0d0
+  do i=1,nw
+    wp(i)=(w0+wrange/dble(nw)*dble(i-1))/27.211385d0
+  end do
 
-! Obtain SP Kubo
-do ibz=1,npointstotal
-  e(:) = eigval_stack(:, ibz)
-
-  !get strength for kubo SP
-  call get_kubo_intens(nv_ex,npointstotal,vcell,nbands,e,vme(ibz, :, :, :),nw,wp,sigma_w_sp,eta)
-
-end do
+  allocate (sigma_w_ex(3,3,nw))
+  allocate (skubo_ex_int(3,3,norb_ex_cut))
+  sigma_w_ex=0.0d0
+  skubo_ex_int=0.0d0
 
+  ! Obtain SP Kubo
+  do ibz=1,npointstotal
+    e(:) = eigval_stack(:, ibz)
+    call get_kubo_intens(nv_ex,npointstotal,vcell,nbands,e,vme(ibz, :, :, :),nw,wp,sigma_w_sp,eta)
+  end do
 
-!fill kubo oscillators of EXCITONS
-do nn=1,norb_ex_cut
-  !save oscillator stregths
-  do nj=1,3
-    do njp=1,3
-      skubo_ex_int(nj,njp,nn)=pi/(dble(npointstotal)*vcell) &
-      *conjg(vme_ex(nj,nn,1))*vme_ex(njp,nn,1)/e_ex(nn)   !pick the correct order of operators
+  !fill kubo oscillators of EXCITONS
+  do nn=1,norb_ex_cut
+    do nj=1,3
+      do njp=1,3
+        skubo_ex_int(nj,njp,nn)=pi/(dble(npointstotal)*vcell) &
+        *conjg(vme_ex(nj,nn,1))*vme_ex(njp,nn,1)/e_ex(nn)
+      end do
     end do
   end do
-end do
-
 
-!excitons
-do ialpha=1,3
-  do ialphap=1,3
-    call broad_vector(type_broad,norb_ex_cut,e_ex,skubo_ex_int(ialpha,ialphap,:), &
-    nw,wp,sigma_w_ex(ialpha,ialphap,:),eta)
+  !apply broadening to exciton contributions
+  do ialpha=1,3
+    do ialphap=1,3
+      call broad_vector(type_broad,norb_ex_cut,e_ex,skubo_ex_int(ialpha,ialphap,:), &
+      nw,wp,sigma_w_ex(ialpha,ialphap,:),eta)
+    end do
   end do
-end do
 
-!write frequency dependent conductivity
-open(50,file=file_name_sp)
-open(60,file=file_name_ex)
-do iw=1,nw
-  feps=1.0d0
-  !feps=4.0d0*pi*1.0d0/137.035999084d0*100.0d0   !absorbance units
-  !eps=4.0d0   !\sigma_0 units
-  write(50,*) wp(iw)*27.211385d0, &
-              -realpart(feps*sigma_w_sp(1,1,iw)), &
-              -realpart(feps*sigma_w_sp(1,2,iw)), &
-              -realpart(feps*sigma_w_sp(1,3,iw)), &
-              -realpart(feps*sigma_w_sp(2,1,iw)), &
-              -realpart(feps*sigma_w_sp(2,2,iw)), &
-              -realpart(feps*sigma_w_sp(2,3,iw)), &
-              -realpart(feps*sigma_w_sp(3,1,iw)), &
-              -realpart(feps*sigma_w_sp(3,2,iw)), &
-              -realpart(feps*sigma_w_sp(3,3,iw))
-  write(60,*) wp(iw)*27.211385d0, &
-              realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(1,1,iw)), &
-              realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(1,2,iw)), &
-              realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(1,3,iw)), &
-              realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(2,1,iw)), &
-              realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(2,2,iw)), &
-              realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(2,3,iw)), &
-              realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(3,1,iw)), &
-              realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(3,2,iw)), &
-              realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(3,3,iw))
-end do
+  !write frequency dependent conductivity
+  open(50,file=file_name_sp)
+  open(60,file=file_name_ex)
+  do iw=1,nw
+    feps=1.0d0
+    write(50,*) wp(iw)*27.211385d0, &
+                -realpart(feps*sigma_w_sp(1,1,iw)), &
+                -realpart(feps*sigma_w_sp(1,2,iw)), &
+                -realpart(feps*sigma_w_sp(1,3,iw)), &
+                -realpart(feps*sigma_w_sp(2,1,iw)), &
+                -realpart(feps*sigma_w_sp(2,2,iw)), &
+                -realpart(feps*sigma_w_sp(2,3,iw)), &
+                -realpart(feps*sigma_w_sp(3,1,iw)), &
+                -realpart(feps*sigma_w_sp(3,2,iw)), &
+                -realpart(feps*sigma_w_sp(3,3,iw))
+    write(60,*) wp(iw)*27.211385d0, &
+                realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(1,1,iw)), &
+                realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(1,2,iw)), &
+                realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(1,3,iw)), &
+                realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(2,1,iw)), &
+                realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(2,2,iw)), &
+                realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(2,3,iw)), &
+                realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(3,1,iw)), &
+                realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(3,2,iw)), &
+                realpart(complex(0.0d0,-1.0d0)*feps*sigma_w_ex(3,3,iw))
+  end do
+  close(50)
+  close(60)
+endif
 
-close(50)
-close(60)
 
 p_dot = scan(file_name_ex, '.', .true.)           ! find first “.” from the right
 if (p_dot == 0) then
@@ -182,7 +209,7 @@ subroutine skubo_w(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,hhop,shop,npointstota
 do iw=1,nw
   feps=1.0d0
   !feps=4.0d0*pi*1.0d0/137.035999084d0*100.0d0   !absorbance units
-  !eps=4.0d0   !\sigma_0 units
+  !feps=4.0d0   !\sigma_0 units
   write(50,*) wp(iw)*27.211385d0, &
               -imagpart(feps*sigma_w_sp(1,1,iw)), &
               -imagpart(feps*sigma_w_sp(1,2,iw)), &
@@ -311,9 +338,30 @@ subroutine exciton_oscillator_strength(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,h
   pi=acos(-1.0d0)
 
   !get unit cell volume
-  call crossproduct(R(1,1),R(1,2),R(1,3),R(2,1), &
-  R(2,2),R(2,3),cx,cy,cz)
-  vcell=sqrt(cx**2+cy**2+cz**2)
+  ! volume formula depends on the dimensionality of the unit cell 
+  ! for 1D, we have the magnitude of the vector (we take the magnitude of the three vectors just to be safe)
+  ! for 2D, we have the crossproduct (as before, now allowing the lattice to also be oriented in the yz or xz planes)
+  ! for 3D, we have the triple product dot(z,cross(x,y))
+  ! this is checked by checking the norm of rkx/rky/rkz to determine how many directions are defined in k
+
+  if ( (NORM2(rkx) == 0 .and. NORM2(rky) == 0) .or. &
+       (NORM2(rky) == 0 .and. NORM2(rkz) == 0) .or. &
+       (NORM2(rkx) == 0 .and. NORM2(rkz) == 0) ) then
+    vcell=sqrt(R(1,1)**2+R(1,2)**2+R(1,3)**2+R(2,1)**2+R(2,2)**2+R(2,3)**2+R(3,1)**2+R(3,2)**2+R(3,3)**2)
+  else if ( ( NORM2(rkz) == 0 .xor. NORM2(rky) == 0) .or. &
+            (NORM2(rkx) == 0 .xor. NORM2(rky) == 0) ) then
+    if (NORM2(rkx) == 0) then
+      call crossproduct(R(3,1),R(3,2),R(3,3),R(2,1),R(2,2),R(2,3),cx,cy,cz)
+    else if (NORM2(rky) == 0) then
+      call crossproduct(R(1,1),R(1,2),R(1,3),R(3,1),R(3,2),R(3,3),cx,cy,cz)
+    else
+      call crossproduct(R(1,1),R(1,2),R(1,3),R(2,1),R(2,2),R(2,3),cx,cy,cz)
+    endif
+    vcell=sqrt(cx**2+cy**2+cz**2)
+  else 
+    call crossproduct(R(1,1),R(1,2),R(1,3),R(2,1),R(2,2),R(2,3),cx,cy,cz)
+    vcell=abs(R(3,1)*cx+R(3,2)*cy+R(3,3)*cz)
+  endif
 
   !SP arrays
   allocate (sderhop(3,nR,norb,norb))
@@ -350,7 +398,7 @@ subroutine exciton_oscillator_strength(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,h
     end do
 
   !getting some SP variables
-  call hoppings_observables_TB(norb,nR,Rvec,shop,hhop,rhop,sderhop,hderhop)
+  call hoppings_observables_TB(norb,nR,Rvec,shop,hhop,rhop,sderhop,hderhop,NORM2(rkx),NORM2(rky),NORM2(rkz))
   !11/05/2023 JJEP: fill rhop here. Easier to extend to DFT later
 
   !Brillouin zone sampling
@@ -372,7 +420,7 @@ subroutine exciton_oscillator_strength(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,h
     hhop,rhop,sderhop,hderhop,sderkernel,hderkernel,akernel,pgaugekernel)
     !velocity matrix elements
     call get_eigen_vme(norb,nbands,skernel,hkernel,akernel,hderkernel, &
-    pgaugekernel,hk_ev,e,pgauge,vjseudoa,vjseudob,vme(ibz, :, :, :))
+    pgaugekernel,hk_ev,e,pgauge,vjseudoa,vjseudob,vme(ibz, :, :, :),NORM2(rkx),NORM2(rky),NORM2(rkz))
 
     !fill V_N
     !!$OMP critical
@@ -380,18 +428,18 @@ subroutine exciton_oscillator_strength(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,h
       do iex=1,norb_ex_cut
         !!! Iterate over band index explicitly (AJU 03-06-23)
         do ic=1,nc_ex
-        do iv=1,nv_ex
+          do iv=1,nv_ex
 
-          j = nc_ex * nv_ex * (ibz - 1) + nv_ex * (ic - 1) + iv
-          !get valence/conduction indices
-          !nv_ip=(nv-nv_ex)+ib-int((ib-1)/nv_ex)*nv_ex
-          !nc_ip=(nv+1)+int((ib-1)/nv_ex)
-          !j=(ibz-1)*norb_ex_band+ib
+            j = nc_ex * nv_ex * (ibz - 1) + nv_ex * (ic - 1) + iv
+            !get valence/conduction indices
+            !nv_ip=(nv-nv_ex)+ib-int((ib-1)/nv_ex)*nv_ex
+            !nc_ip=(nv+1)+int((ib-1)/nv_ex)
+            !j=(ibz-1)*norb_ex_band+ib
 
-          vme_ex(nj,iex,1)=vme_ex(nj,iex,1)+fk_ex(j,iex)*vme(ibz, nj,iv,ic + nv_ex)
-          vme_ex(nj,iex,2)=vme_ex(nj,iex,2)+fk_ex(j,iex)*vme(ibz, nj,ic + nv_ex,iv)
+            vme_ex(nj,iex,1)=vme_ex(nj,iex,1)+fk_ex(j,iex)*vme(ibz, nj,iv,ic + nv_ex)
+            vme_ex(nj,iex,2)=vme_ex(nj,iex,2)+fk_ex(j,iex)*vme(ibz, nj,ic + nv_ex,iv)
 
-        end do
+          end do
         end do
       end do
     end do
@@ -453,8 +501,22 @@ subroutine get_kubo_parameters(w0,wrange,nw,eta,type_broad, &
 character(100) type_broad
 character(100) file_name_sp
 character(100) file_name_ex
+logical :: exists_file
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! make sure the input file exists before trying to read it
+inquire(file='kubo_w.in', exist=exists_file)
+if (.not. exists_file) then
+  write(*,*) 'Warning: kubo_w.in not found, skipping Kubo conductivity'
+  nw = 0
+  w0 = 0.0d0
+  wrange = 0.0d0
+  eta = 0.0d0
+  type_broad = 'lorentzian'
+  file_name_sp = 'kubo_sp.dat'
+  file_name_ex = 'kubo_ex.dat'
+  return
+endif
 open(10,file='kubo_w.in')
 read(10,*)
 read(10,*) w0
@@ -477,7 +539,7 @@ subroutine get_kubo_parameters(w0,wrange,nw,eta,type_broad, &
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine hoppings_observables_TB(norb,nR,Rvec,shop,hhop, &
-rhop,sderhop,hderhop)
+rhop,sderhop,hderhop,nrkx,nrky,nrkz)
 
 implicit real*8 (a-h,o-z)
 dimension Rvec(nR,3), rhop(3,nR,norb,norb)
@@ -486,34 +548,36 @@ subroutine hoppings_observables_TB(norb,nR,Rvec,shop,hhop, &
 
 complex*16 hhop
 complex*16 hderhop,sderhop
+real*8 nrkx
+real*8 nrky
+real*8 nrkz
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 hderhop=0.0d0
 sderhop=0.0d0
 do iR=1,nR
   do ialpha=1,norb
-  do ialphap=1,ialpha
-    Rx=Rvec(iR,1)
-    Ry=Rvec(iR,2)
+    do ialphap=1,ialpha
 
-    ! Check if the z direction is defined (thus periodic)
-    if (Rvec(iR,3) /= 0) then
-      Rz = Rvec(iR,3)
-    else
-      Rz = rhop(3, iR, ialpha,ialphap) ! Quintela et. al. (2023) DOI: 10.1103/PhysRevB.107.235416
-    end if
+      ! Check if each direction is defined (thus periodic)
+      ! Quintela et. al. (2023) DOI: 10.1103/PhysRevB.107.235416
 
+      Rx = merge(Rvec(iR,1), rhop(1,iR,ialpha,ialphap), nrkx /= 0)
+      Ry = merge(Rvec(iR,2), rhop(2,iR,ialpha,ialphap), nrky /= 0)
+      Rz = merge(Rvec(iR,3), rhop(3,iR,ialpha,ialphap), nrkz /= 0)
 
-    hderhop(1,iR,ialpha,ialphap)=complex(0.0d0,Rx)*hhop(ialpha,ialphap,iR)
-    hderhop(2,iR,ialpha,ialphap)=complex(0.0d0,Ry)*hhop(ialpha,ialphap,iR)
-    hderhop(3,iR,ialpha,ialphap)=complex(0.0d0,Rz)*hhop(ialpha,ialphap,iR)
 
-    sderhop(1,iR,ialpha,ialphap)=complex(0.0d0,Rx)*shop(ialpha,ialphap,iR)
-    sderhop(2,iR,ialpha,ialphap)=complex(0.0d0,Ry)*shop(ialpha,ialphap,iR)
-    sderhop(3,iR,ialpha,ialphap)=complex(0.0d0,Rz)*shop(ialpha,ialphap,iR)
+      hderhop(1,iR,ialpha,ialphap)=complex(0.0d0,Rx)*hhop(ialpha,ialphap,iR)
+      hderhop(2,iR,ialpha,ialphap)=complex(0.0d0,Ry)*hhop(ialpha,ialphap,iR)
+      hderhop(3,iR,ialpha,ialphap)=complex(0.0d0,Rz)*hhop(ialpha,ialphap,iR)
 
+      sderhop(1,iR,ialpha,ialphap)=complex(0.0d0,Rx)*shop(ialpha,ialphap,iR)
+      sderhop(2,iR,ialpha,ialphap)=complex(0.0d0,Ry)*shop(ialpha,ialphap,iR)
+      sderhop(3,iR,ialpha,ialphap)=complex(0.0d0,Rz)*shop(ialpha,ialphap,iR)
 
+
+    end do
   end do
-  end do
+
 end do
 
 
@@ -560,14 +624,18 @@ subroutine get_vme_kernels(rkx,rky,rkz,nR,Rvec,norb, &
   do ialphap=1,ialpha
 
     do iRp=1,nR
-      Rx=Rvec(iRp,1)
-      Ry=Rvec(iRp,2)
-
-      if (Rvec(iRp, 3) /= 0) then
-        Rz=Rvec(iRp,3)
-      else
-        Rz=0.0d0
-      end if
+!       Rx=Rvec(iRp,1)
+!       Ry=Rvec(iRp,2)
+! 
+!       if (Rvec(iRp, 3) /= 0) then
+!         Rz=Rvec(iRp,3)
+!       else
+!         Rz=0.0d0
+!       end if
+      Rx = merge(Rvec(iRp,1), 0.0d0, Rvec(iRp,1) /= 0.0)
+      Ry = merge(Rvec(iRp,2), 0.0d0, Rvec(iRp,2) /= 0.0)
+      Rz = merge(Rvec(iRp,3), 0.0d0, Rvec(iRp,3) /= 0.0)
+      
 
       phase=complex(0.0d0,rkx*Rx+rky*Ry+rkz*Rz)
       factor=exp(phase)
@@ -616,7 +684,7 @@ subroutine get_vme_kernels(rkx,rky,rkz,nR,Rvec,norb, &
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine get_eigen_vme(norb,nbands,skernel, &
 hkernel,akernel,hderkernel,pgaugekernel, &
-hk_ev,e,pgauge,vjseudoa,vjseudob,vme)
+hk_ev,e,pgauge,vjseudoa,vjseudob,vme,nrkx,nrky,nrkz)
 
 implicit real*8 (a-h,o-z)
 
@@ -628,41 +696,17 @@ subroutine get_eigen_vme(norb,nbands,skernel, &
 dimension ecomplex(norb)
 dimension vjseudoa(3,nbands,nbands),vjseudob(3,nbands,nbands),vme(3,nbands,nbands)
 
+
 complex*16 ecomplex
 complex*16 skernel,hkernel,akernel,hderkernel
 complex*16 hk_ev,vjseudoa,vjseudob,vme,pgaugekernel,pgauge
+real*8 nrkx
+real*8 nrky
+real*8 nrkz
 
 complex*16 amu,amup
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-!!! Remove diagonalization; instead use stored eigvecs and eigvals (AJU 29-05-23)
-! e=0.0d0
-! call diagoz(norb,e,hkernel)
-! !call diagoz_arbg(norb,ecomplex,skernel,hkernel)
-! !call order(norb,ecomplex,e,hkernel)
-! do ii=1,norb
-! do jj=1,norb
-!   hk_ev(ii,jj)=hkernel(ii,jj)
-! end do
-! end do
-
-! !phase election
-! do j=1,norb
-!   aux1=0.0d0
-!   do i=1,norb
-!     aux1=aux1+hk_ev(i,j)
-!   end do
-
-!   !argument of the sym
-!   arg=atan2(aimag(aux1),realpart(aux1))
-!   factor=exp(complex(0.0d0,-arg))
-!   !write(*,*) 'sum is now:',aux1*factor
-!   do ii=1,norb
-!     hk_ev(ii,j)=hk_ev(ii,j)*factor
-! 	!hk_ev(ii,j)=hk_ev(ii,j)*1.0d0
-!   end do
-! end do
-
 !write(*,*) 'computing velocity matrix elements'
 vme=0.0d0
 vjseudoa=0.0d0
@@ -681,8 +725,8 @@ subroutine get_eigen_vme(norb,nbands,skernel, &
     do nj=1,3
       pgauge(nj,nn,nnp)=pgauge(nj,nn,nnp)+ &
       conjg(amu)*amup*pgaugekernel(nj,ialpha,ialphap)
-
-      if (nj == 3) then
+      
+      if ( ( nrkx == 0 .and. nj == 1) .or. ( nrky == 0 .and. nj == 2) .or. ( nrkz == 0 .and. nj == 3) )then
         vjseudoa(nj,nn,nnp)=vjseudoa(nj,nn,nnp)+ &
         conjg(amu)*amup*hderkernel(nj,ialpha,ialphap)*(e(nnp)-e(nn))
       else
@@ -752,7 +796,7 @@ subroutine get_kubo_intens(nv,npointstotal,vcell,nbands,e,vme,nw,wp,sigma_w_sp,e
 	  !lorentzian
       delta_nnp=1.0d0/pi*aimag(1.0d0/(wp(iw)-e(nn)+e(nnp)-complex(0.0d0,eta)))
 	  !exponential
-    ! delta_nnp=1.0d0/eta*1.0d0/sqrt(2.0d0*pi)*exp(-0.5d0/(eta**2)*(wp(iw)-e(nn)+e(nnp))**2)
+     !delta_nnp=1.0d0/eta*1.0d0/sqrt(2.0d0*pi)*exp(-0.5d0/(eta**2)*(wp(iw)-e(nn)+e(nnp))**2)
 
       !save oscillator stregths
       do nj=1,3
@@ -819,7 +863,6 @@ subroutine crossproduct(ax,ay,az,bx,by,bz,cx,cy,cz)
 cx=ay*bz-az*by
 cy=az*bx-ax*bz
 cz=ax*by-ay*bx
-
 return
 end
 
@@ -865,7 +908,7 @@ subroutine diagoz(n,w,h)
 !     Juan Jose Esteve-Paredes                24.04.2017
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-      subroutine diagoz_arbg(n,w,s,h)
+subroutine diagoz_arbg(n,w,s,h)
 !finding the eigenvalues of a complex matrix using LAPACK
       implicit real*8 (a-h,o-z)
 !declarations, notice double precision
@@ -880,40 +923,10 @@ subroutine diagoz_arbg(n,w,s,h)
       saux=s
       haux=h
       !find the solution using the LAPACK routine ZGGEEV
-      !e=0.0d0
-      !do i=1,3
-        !write(*,*) (s(i,j),j=1,3)
-      !end do
-      !pause
-      !call zggev('V','V',n,s,n,h,n,ALPHA,BETA,VL,n,VR,n,WORK,200*n, RWORK,INFO)
       call zhegv( 1, 'V', 'U', n, h, n, s, n, wreal, WORK, 200*n, RWORK, INFO )
-      !do i=1,n
-        !write(*,*) wreal(i),1.0d0/wreal(i)
-      !end do
-      !do i=1,n
-        !w(i)=ALPHA(i)/BETA(i)
-        !write(*,*) ALPHA(i),BETA(i),w(i)
-        !do j=1,n
-          !h(i,j)=VR(i,j)
-        !end do
-      !end do
-      !pause
-
-
-      !ANOTHER SIMILAR TO THE PREVIOUS ONE
-      !call zgegv( 'N', 'V', n, h, n, s, n, ALPHA, BETA, &
-      !QVL, n, VR, n, WORK, 200*n, RWORK, INFO )
-      !write(*,*) 'Diagonalization=',INFO
-      !write(*,*) INFO
+      
       do i=1,n
-        !w(i)=ALPHA(i)/BETA(i)
         w(i)=complex(wreal(i),0.0d0)
-        !write(*,*) ALPHA(i),BETA(i),w(i)
-        !do j=1,n
-          !h(i,j)=VR(i,j)
-        !end do
-        !write(*,*) wreal(i),h(i,1)
-        !write(*,*) ALPHA(i),BETA(i),w(i)
       end do
       !pause
       do k=1,n
@@ -958,27 +971,27 @@ subroutine diagoz_arbg(n,w,s,h)
       !pause
 
       return
-      end
+end
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !subroutine to order the complez eigenvalues and making them real
      !it also orders the eigenvectors
-     subroutine order(n,ecomplex,e,hk)
-     implicit real*8 (a-h,o-z)
-     complex*16 ecomplex,hk,hkaux
-     dimension ecomplex(n),e(n),eaux(n),hk(n,n),hkaux(n,n)
+subroutine order(n,ecomplex,e,hk)
+  implicit real*8 (a-h,o-z)
+  complex*16 ecomplex,hk,hkaux
+  dimension ecomplex(n),e(n),eaux(n),hk(n,n),hkaux(n,n)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-     hkaux=hk
-     do i=1,n
-       eaux(i)=realpart(ecomplex(i))
-     end do
-
-     do i=1,n
-       imin=minloc(eaux,1)
-       e(i)=eaux(imin)
-       hk(:,i)=hkaux(:,imin)
-       eaux(imin)=1.0d8
-     end do
-
-     return
-     end
+  hkaux=hk
+  do i=1,n
+    eaux(i)=realpart(ecomplex(i))
+  end do
+
+  do i=1,n
+    imin=minloc(eaux,1)
+    e(i)=eaux(imin)
+    hk(:,i)=hkaux(:,imin)
+    eaux(imin)=1.0d8
+  end do
+
+  return
+end
diff --git a/test/tests.cpp b/test/tests.cpp
index 934c905..593b4f6 100644
--- a/test/tests.cpp
+++ b/test/tests.cpp
@@ -151,8 +151,8 @@ TEST(FileParsing, W90fileParsing){
     // arma::vec expectedDegenHash = {1, 1, 2, 1, 1};
 
     // Hashes for Fock matrices (real/imag parts handled separately)
-    arma::vec expectedHamiltonianHash = {0.978634, 0.990920, 1.046875, 1.054984, 0.996068, 0.968801, 1.101170, 0.961007,
-        0.968796, 1.043775, 1.032634, 1.017193, 0.940727, 1.172638, 0.890277, 1.110180,
+    arma::vec expectedHamiltonianHash = {0.978634, 0.990920, 1.023349, 1.054984, 0.996068, 0.968801, 1.101170, 0.961007,
+        0.968796, 1.021764, 1.032634, 1.017193, 0.940727, 1.172638, 0.890277, 1.110180,
         1.088959, 1.027729, 0.962491, 0.970735, 1.036763, 1.089453, 1.079106, -0.531977,
         0.946280, 0.965448, 1.031659, 1.011340, 1.016294, 1.009594, 0.999775, 0.951863,
         0.458796, 0.579032, 1.125541, 1.059755, 0.985202, 0.986659, 1.042266, 1.042844,
@@ -160,8 +160,8 @@ TEST(FileParsing, W90fileParsing){
         1.131210, 2.914045, 6.239713, 0.937042, 1.013736, 0.995272, 1.019086, 0.998125,
         1.032483, 0.973986, 1.026194, 1.357262, 1.224249, 1.032389, 1.027929, 0.987408,
         0.982629, 1.003563, 1.071731, 1.044472, 0.880948, 1.118926, 0.963565, 0.998474,
-        1.029212, 1.046875, 0.982783, 0.964047, 1.057069, 0.999934, 0.973145, 1.045742,
-        1.043775, 1.005556, 0.966527
+        1.029212, 1.023349, 0.982783, 0.964047, 1.057069, 0.999934, 0.973145, 1.045742,
+        1.021764, 1.005556, 0.966527
     };
     
     // Hashes for Rhop matrices (if used)
@@ -679,8 +679,8 @@ TEST(hBNTest, WannierCalculation){
 
     auto energies = xatu::detectDegeneracies(results->eigval, nstates, 6);
     
-    std::vector<std::vector<double>> expectedEnergies = {{2.833347, 1}, 
-                                                         {2.834596, 1}};
+    std::vector<std::vector<double>> expectedEnergies = {{2.829980, 1}, 
+                                                         {2.831276, 1}};
                                                          
     for(uint i = 0; i < energies.size(); i++){
         EXPECT_NEAR(energies[i][0], expectedEnergies[i][0], 1E-4);
@@ -692,7 +692,7 @@ TEST(hBNTest, WannierCalculation){
     arma::mat norm_vme_ex = arma::square(arma::abs(vme_ex));
     double cum_norm_vme_ex = arma::accu(norm_vme_ex);
 
-    double expectedTotalOscillator = 103.0325881;
+    double expectedTotalOscillator = 102.7178826790;
     EXPECT_NEAR(cum_norm_vme_ex, expectedTotalOscillator, 1E-6);
 
     // Check reciprocal w.f.
@@ -713,7 +713,7 @@ TEST(hBNTest, WannierCalculation){
 
 
     double kwfHash = xatu::array2hash(kwf);
-    double expectedKwfHash = 1.057666;
+    double expectedKwfHash = 1.0576659386;
     EXPECT_NEAR(kwfHash, expectedKwfHash, 1E-5);
 
     // Check realspace w.f.
@@ -737,7 +737,7 @@ TEST(hBNTest, WannierCalculation){
     }
 
     double rswfHash = xatu::array2hash(rswf);
-    double expectedRSwfHash = 1.003540;
+    double expectedRSwfHash = 1.0035186777;
     EXPECT_NEAR(rswfHash, expectedRSwfHash, 1E-5);
 
 }
@@ -1046,4 +1046,4 @@ TEST(MoS2Test, Conductivity){
 int main(int argc, char **argv) {
     ::testing::InitGoogleTest(&argc, argv);
     return RUN_ALL_TESTS();
-}
\ No newline at end of file
+}
diff --git a/utility/wannier2xatu/utils.f90 b/utility/wannier2xatu/utils.f90
index 617de4a..a2356c0 100644
--- a/utility/wannier2xatu/utils.f90
+++ b/utility/wannier2xatu/utils.f90
@@ -157,7 +157,7 @@ subroutine Export2Xatu
             ! ------------------------------------------------------------------------------------ !
                 write(iunit, '(A)') '# hamiltonian'
                 do i=1, nFock
-                    H(i,:,:) = H(i,:,:)*Degen(i)
+                    H(i,:,:) = H(i,:,:)/Degen(i)
                     do j=1, mSize
                         do k=1,mSize
                             write(iunit, '(F20.15, A, F20.15, A)', advance = 'no') real(H(i, j, k)),' ',aimag(H(i, j, k)), 'j    '