From 77dddafc7ac710843a0d1568d2d481485eafe89c Mon Sep 17 00:00:00 2001 From: mfcmquintela Date: Fri, 13 Feb 2026 10:16:57 +0100 Subject: [PATCH 01/10] Energy cutoff for exporting BSE eigensystem Implements a CLI flag "-t" or "--encut" that enhances the previous --states flag for setting a number of states to export by exporting the eigensystem up to an energy set via this flag. Energy cutoff takes priority over number of states cutoff, making it easier to export the states needed to study optical response up to a specific frequency. --- include/xatu/Result.hpp | 6 +++--- include/xatu/Result.tpp | 38 ++++++++++++++++++++++++++++---------- include/xatu/utils.hpp | 12 +++++++++--- main/xatu.cpp | 10 ++++++---- 4 files changed, 46 insertions(+), 20 deletions(-) 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::writeRealspaceAmplitude(int stateindex, int holeIndex, * @return void */ template -void Result::writeEigenvalues(FILE* textfile, int n){ +void Result::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::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::writeEigenvalues(FILE* textfile, int n){ * @return void */ template -void Result::writeStates(FILE* textfile, int n){ - if(n > exciton->excitonbasisdim || n < 0){ +void Result::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::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::writeStates(FILE* textfile, int n){ * @param textfile Textfile where the spins are written. */ template -void Result::writeSpin(int n, FILE* textfile){ - - if(n > exciton->excitonbasisdim || n < 0){ +void Result::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 -void printEnergies(const std::unique_ptr& results, int n = 8, int precision = 6){ +void printEnergies(const std::unique_ptr& 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> 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> 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..1059ddf 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 statesArg("n", "states", "Specify number of exciton states to show.", false, 8, "No. states", cmd); + TCLAP::ValueArg energycutoff("t", "encut", "Specify up to high energy to print excitons (eV).", false, 0.0, "En cutoff", cmd); TCLAP::ValueArg 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 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(); From a0923c56b12b0717979c3d16cef1b0eab12bbdfb Mon Sep 17 00:00:00 2001 From: visitor Date: Tue, 14 Apr 2026 16:43:00 +0200 Subject: [PATCH 02/10] now it always computes and writes the oscillator strength regardless of kubo.in file existing --- src/skubo_w.f90 | 160 +++++++++++++++++++++++++++--------------------- 1 file changed, 90 insertions(+), 70 deletions(-) diff --git a/src/skubo_w.f90 b/src/skubo_w.f90 index c75b1f2..fb4da48 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 @@ -74,96 +75,101 @@ subroutine skubo_w(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,hhop,shop,npointstota !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 @@ -453,8 +459,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 From 17a7c49c2caa9ccb55f0aabe1414d1c62aea5a43 Mon Sep 17 00:00:00 2001 From: visitor Date: Fri, 17 Apr 2026 10:52:11 +0200 Subject: [PATCH 03/10] fixed degeneracy on wannier hamiltonians; now user can choose number of states by energy instead of index --- main/xatu.cpp | 2 +- src/Wannier90Configuration.cpp | 2 +- utility/wannier2xatu/utils.f90 | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/main/xatu.cpp b/main/xatu.cpp index 1059ddf..b4f6e2e 100755 --- a/main/xatu.cpp +++ b/main/xatu.cpp @@ -22,7 +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 statesArg("n", "states", "Specify number of exciton states to show.", false, 8, "No. states", cmd); - TCLAP::ValueArg energycutoff("t", "encut", "Specify up to high energy to print excitons (eV).", false, 0.0, "En cutoff", cmd); + TCLAP::ValueArg energycutoff("t", "ecut", "Specify up to high energy to print excitons (eV).", false, 0.0, "En cutoff", cmd); TCLAP::ValueArg 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 dftArg("d", "dft", "Indicates that the system file is a .outp CRYSTAL file.", false, -1, "No. Fock matrices", cmd); 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/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 ' From 254d11fb939fb5d23ed8b485436ea7190d5aabb3 Mon Sep 17 00:00:00 2001 From: visitor Date: Mon, 11 May 2026 17:11:53 +0200 Subject: [PATCH 04/10] updates tests for new wannier parser with degen --- test/tests.cpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) 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> expectedEnergies = {{2.833347, 1}, - {2.834596, 1}}; + std::vector> 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 +} From 7993a289ab2cd7bb1dd391af1c481390ee5211d7 Mon Sep 17 00:00:00 2001 From: mfcmquintela Date: Tue, 2 Jun 2026 11:31:03 +0200 Subject: [PATCH 05/10] Fixed absorption diverging for 1D systems. Made calculations consistent based on the system's dimensionality --- src/skubo_w.f90 | 242 +++++++++++++++++++++++------------------------- 1 file changed, 115 insertions(+), 127 deletions(-) diff --git a/src/skubo_w.f90 b/src/skubo_w.f90 index fb4da48..3a1f8ed 100644 --- a/src/skubo_w.f90 +++ b/src/skubo_w.f90 @@ -67,9 +67,27 @@ 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, one cannot assume a direction for the periodicity) +! 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) + write(*,*) "1D system" + write(*,*) +else if ( ( NORM2(rkz) == 0 .xor. NORM2(rky) == 0) .or. (NORM2(rkx) == 0 .xor. NORM2(rky) == 0) ) then + 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) + write(*,*) "2D system" + write(*,*) +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=sqrt(R(3,1)*cx+R(3,2)*cy+R(3,3)*cz) + write(*,*) "3D system" + write(*,*) +endif !call fill_nRvec(nR,R,Rvec,nRvec) !get printing parameters @@ -188,7 +206,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)), & @@ -317,9 +335,21 @@ 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, one cannot assume a direction for the periodicity) + ! 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 + 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) + 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=sqrt(R(3,1)*cx+R(3,2)*cy+R(3,3)*cz) + endif !SP arrays allocate (sderhop(3,nR,norb,norb)) @@ -356,7 +386,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 @@ -378,7 +408,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 @@ -386,18 +416,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 @@ -497,7 +527,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) @@ -506,34 +536,43 @@ 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 +! Rx=Rvec(iR,1) +! Ry=Rvec(iR,2) +! if (Rvec(iR,3) /= 0) then +! Rz = Rvec(iR,3) +! else +! Rz = rhop(3, iR, ialpha,ialphap) +! end if - ! 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) and we're not in the first unit cell (there the Rvec is always zero) + ! 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 @@ -580,14 +619,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) @@ -636,7 +679,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) @@ -648,41 +691,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 @@ -701,8 +720,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 @@ -772,7 +791,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 @@ -839,7 +858,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 @@ -885,7 +903,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 @@ -900,40 +918,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 @@ -978,27 +966,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 From 402ea6bc0432d7c24c81638d3ee53dd9221dfdc2 Mon Sep 17 00:00:00 2001 From: mfcmquintela Date: Tue, 2 Jun 2026 11:32:54 +0200 Subject: [PATCH 06/10] More careful cutoffs in Coulomb potential in real space. Added dielectric constant of material parameter --- src/ExcitonTB.cpp | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) 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; } /** From 331666cc353ee26408462f7adf39f013053bc908 Mon Sep 17 00:00:00 2001 From: mfcmquintela Date: Tue, 2 Jun 2026 11:36:50 +0200 Subject: [PATCH 07/10] Removed forgotten debug flag --- src/skubo_w.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/skubo_w.f90 b/src/skubo_w.f90 index 3a1f8ed..fc69d8f 100644 --- a/src/skubo_w.f90 +++ b/src/skubo_w.f90 @@ -75,18 +75,18 @@ subroutine skubo_w(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,hhop,shop,npointstota 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) - write(*,*) "1D system" - write(*,*) +! write(*,*) "1D system" +! write(*,*) else if ( ( NORM2(rkz) == 0 .xor. NORM2(rky) == 0) .or. (NORM2(rkx) == 0 .xor. NORM2(rky) == 0) ) then 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) - write(*,*) "2D system" - write(*,*) +! write(*,*) "2D system" +! write(*,*) 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=sqrt(R(3,1)*cx+R(3,2)*cy+R(3,3)*cz) - write(*,*) "3D system" - write(*,*) +! write(*,*) "3D system" +! write(*,*) endif !call fill_nRvec(nR,R,Rvec,nRvec) From db0c36791f71a284f7a0e00dfddb27f6d872b1e1 Mon Sep 17 00:00:00 2001 From: mfcmquintela Date: Tue, 2 Jun 2026 11:52:49 +0200 Subject: [PATCH 08/10] Fixed (?) failed build in docker due to line length --- src/skubo_w.f90 | 49 +++++++++++++++++++++++++++---------------------- 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/src/skubo_w.f90 b/src/skubo_w.f90 index fc69d8f..4e7d109 100644 --- a/src/skubo_w.f90 +++ b/src/skubo_w.f90 @@ -68,25 +68,28 @@ subroutine skubo_w(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,hhop,shop,npointstota eigval_stack=eigval_stack/27.211385d0 ! 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, one cannot assume a direction for the periodicity) +! 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 +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) -! write(*,*) "1D system" -! write(*,*) -else if ( ( NORM2(rkz) == 0 .xor. NORM2(rky) == 0) .or. (NORM2(rkx) == 0 .xor. NORM2(rky) == 0) ) then - call crossproduct(R(1,1),R(1,2),R(1,3),R(2,1),R(2,2),R(2,3),cx,cy,cz) +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) -! write(*,*) "2D system" -! write(*,*) 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=sqrt(R(3,1)*cx+R(3,2)*cy+R(3,3)*cz) -! write(*,*) "3D system" -! write(*,*) endif !call fill_nRvec(nR,R,Rvec,nRvec) @@ -336,15 +339,24 @@ subroutine exciton_oscillator_strength(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,h !get unit cell volume ! 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, one cannot assume a direction for the periodicity) + ! 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 + 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 - call crossproduct(R(1,1),R(1,2),R(1,3),R(2,1),R(2,2),R(2,3),cx,cy,cz) + 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) @@ -545,15 +557,8 @@ subroutine hoppings_observables_TB(norb,nR,Rvec,shop,hhop, & do iR=1,nR do ialpha=1,norb do ialphap=1,ialpha -! Rx=Rvec(iR,1) -! Ry=Rvec(iR,2) -! if (Rvec(iR,3) /= 0) then -! Rz = Rvec(iR,3) -! else -! Rz = rhop(3, iR, ialpha,ialphap) -! end if - ! Check if each direction is defined (thus periodic) and we're not in the first unit cell (there the Rvec is always zero) + ! 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) From 18e34999bb8392f7a15de24fadae03cde6e24eca Mon Sep 17 00:00:00 2001 From: mfcmquintela Date: Fri, 5 Jun 2026 15:52:48 +0200 Subject: [PATCH 09/10] 3D volume should be magnitude of triple product, not sqrt of said magnitude --- src/skubo_w.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/skubo_w.f90 b/src/skubo_w.f90 index 4e7d109..5d00c5c 100644 --- a/src/skubo_w.f90 +++ b/src/skubo_w.f90 @@ -89,7 +89,7 @@ subroutine skubo_w(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,hhop,shop,npointstota 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=sqrt(R(3,1)*cx+R(3,2)*cy+R(3,3)*cz) + vcell=abs(R(3,1)*cx+R(3,2)*cy+R(3,3)*cz) endif !call fill_nRvec(nR,R,Rvec,nRvec) From be9b8c8d0fd1b089768ca1a0dc0b00697591c209 Mon Sep 17 00:00:00 2001 From: mfcmquintela Date: Sat, 6 Jun 2026 01:09:44 +0200 Subject: [PATCH 10/10] Vcell is defined twice, I only fixed the bug in one of them --- src/skubo_w.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/skubo_w.f90 b/src/skubo_w.f90 index 5d00c5c..b831163 100644 --- a/src/skubo_w.f90 +++ b/src/skubo_w.f90 @@ -360,7 +360,7 @@ subroutine exciton_oscillator_strength(nR,norb,norb_ex,nv_ex,nc_ex,nv,Rvec,R,B,h 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=sqrt(R(3,1)*cx+R(3,2)*cy+R(3,3)*cz) + vcell=abs(R(3,1)*cx+R(3,2)*cy+R(3,3)*cz) endif !SP arrays