add outputting eigenvectors function

README.md

@@ -59,6 +59,7 @@ The code is written in modern Fortran with modular design, separating physical m
    JOB B
    FPOSCAR cont.vasp
    FEIGEN tb_band.dat
+   EIGNVEC .false.
    FKPOINTS KPOINTS.band
    NCACHE 1
    ```
@@ -105,6 +106,7 @@ mpirun -n 4 ./tbsolver
 ### Output Files
 Depending on the job type, the program produces:
 - **Band structure**: File specified by `FEIGEN` (default `tb_band.dat`) contains k‑points along the path and corresponding eigenvalues.
+- **Wavefunctions / eigenvectors**: If `EIGNVEC .true.` is set, the program also writes the eigenvectors (wavefunctions) to the file specified by `FEIGENVEC` (default `tb_wavef.dat`). The output follows the selected band range defined by `IBAND`.
 - **DOS results**: File specified by `FDOS` (default `tb_dos.dat`) contains energy and DOS columns. With `DOS_NORMALIZE .true.`, the DOS is divided by the number of k points and is reported in states/eV per cell.
 - **EELS results**: Loss function and related quantities are written to files with names derived from the `QTAG` keyword.
 
@@ -117,6 +119,8 @@ All keywords are read by the parser in `src/parser.f90`. Lines starting with `!`
 | **JOB** | `B`, `D`, or `E` | Type of calculation:<br>`B` – Band structure<br>`D` – Density of States (DOS)<br>`E` – Electron Energy Loss Spectroscopy (EELS) |
 | **FPOSCAR** | `<filename>` | Path to the VASP‑format POSCAR file (default: `cont.vasp`) |
 | **FEIGEN** | `<filename>` | Output file for band eigenvalues (default: `tb_band.dat`). Also enables band‑structure writing. |
+| **EIGNVEC** | `<logical>` | Whether to write eigenvectors / wavefunctions for band calculations (default: `.false.`). When enabled, wavefunctions are saved together with the selected band range. |
+| **FEIGENVEC** | `<filename>` | Output file for eigenvector / wavefunction data (default: `tb_wavef.dat`). |
 | **FDOS** | `<filename>` | Output file for DOS data (default: `tb_dos.dat`). |
 | **FKPOINTS** | `<filename>` | VASP‑format file containing k‑point definitions (default: `KPOINTS`). Line mode is used for band paths; Gamma or Monkhorst-Pack mesh mode is used for DOS and EELS meshes. |
 | **FWANNIER** | `<filename>` | Path to the Wannier90 input file (default: `wannier90`). If present, sets model type to Wannier (`wan`). |
@@ -168,10 +172,21 @@ JOB B
 FPOSCAR cont.vasp
 FKPOINTS KPOINTS.band
 FEIGEN tb_band.dat
+EIGNVEC .true.
+FEIGENVEC tb_wavef.dat
 IBAND 1 20
 NCACHE 1
 ```
 
+The wavefunction output file is plain text. It is grouped by k point and band:
+- Header line: `# Wavefunctions: nkpts=... nions=... bands=...`
+- For each block: a comment line `# kpoint <ik> <band>`
+- Then `nions` lines, each containing two columns:
+   - real part of the coefficient
+   - imaginary part of the coefficient
+
+This format is convenient for post-processing in Python, MATLAB, or Fortran.
+
 Example `KPOINTS.band`:
 ```plaintext
 High-symmetry path

src/constants.f90

@@ -35,13 +35,16 @@ module constants
     logical :: kpt_select
     logical :: job_band
     logical :: write_band
+    logical :: eigenvec_out
     logical :: write_matrix
     character(len=3) :: model_type
     character(len=64) :: f_input
     character(len=64) :: f_poscar
     character(len=64) :: f_eig
+    character(len=64) :: f_eigvec
     character(len=64) :: f_kpoint
     character(len=64) :: f_wannier
+    character(len=64) :: f_hr
     character(len=64) :: f_matrix
     character(len=64) :: f_eels
     character(len=64) :: f_dos
@@ -60,12 +63,15 @@ subroutine init_constants()
         kpt_select = .false.
         job_band = .false.
         write_band = .false.
+        write_hr = .false.
         write_matrix = .false.
         f_input='input.in'
         f_poscar='cont.vasp'
         f_eig='tb_band.dat'
+        f_eigvec='tb_wavef.dat'
         f_kpoint='KPOINTS'
         f_matrix='tb_eels.mat'
+        f_hr='tb_hr.dat'
         f_eels='tb_eels.dat'
         f_dos='tb_dos.dat'
         nk_path = 10
@@ -74,6 +80,7 @@ subroutine init_constants()
         q_tag = "D"
         delta = 1e-3
         calc_iqr = .false.
+        eigenvec_out = .false.
         model_type = "tbg"
         f_wannier = 'wannier90'
         eels_mode = 0

src/ioutils.f90

@@ -105,19 +105,21 @@ subroutine parsePOSCAR(filename)
 
     end subroutine parsePOSCAR
 
-    subroutine readKPOINTS(filename,klist_frac)
+    subroutine readKPOINTS(filename,klist_frac,is_path_mode)
         use constants, only:prec,f_kpoint
         type :: node
             real(prec) :: data(3)
             type(node), pointer :: next => null()
         end type node
         character(len=64), optional, intent(in):: filename
         real(prec), allocatable,optional, intent(out) :: klist_frac(:,:)
+        logical, optional, intent(out) :: is_path_mode
         character(len=64) :: fin,mode,tag
         character(len=128) :: line
-        integer :: i,unit,nk(3),nkpts
+        integer :: i,unit,nk(3),nkpts,declared_nkpts,io
         type(node), pointer :: head, current
-        real(prec) :: values(3)
+        real(prec) :: values(3), weight
+        logical :: path_mode, cartesian_input
 
         unit = 102
 
@@ -128,12 +130,16 @@ subroutine readKPOINTS(filename,klist_frac)
         end if
 
         nkpts = 0
+        declared_nkpts = 0
+        path_mode = .false.
+        cartesian_input = .false.
         open(unit, file=fin, status='old', action='read')
             read(unit, '(A)') line
-            read(unit, '(I)') nk(1)
+            read(unit, *) nk(1)
             if (nk(1) == 0) then
                 ! determine k-mesh automatically
                 read(unit, '(A)') mode
+                mode = adjustl(mode)
                 read(unit, *) nk
                 if (mode(1:1)=='G' .or. mode(1:1)=='g') then
                     call generate_mesh_gamma(nk,klist_frac)
@@ -142,21 +148,57 @@ subroutine readKPOINTS(filename,klist_frac)
                 end if
             else 
                 read(unit, '(A)') mode
+                mode = adjustl(mode)
                 if (mode(1:1)=='L' .or. mode(1:1)=='l') then
+                    path_mode = .true.
                     nk_path = nk(1)
                     read(unit, '(A)') tag
+                    tag = adjustl(tag)
                     nullify(head)
                     do while (.not. eof(unit))
                         read(unit, '(A)') line
                         if (len_trim(line)==0) cycle ! skip empty lines
-                        read(line,*) values
+                        if (line(1:1)=='!' .or. line(1:1)=='#' .or. line(1:1)=='/') cycle
+                        read(line,*,iostat=io) values
+                        if (io /= 0) cycle
+                        if (tag(1:1)=='C' .or. tag(1:1)=='c') then
+                            values = k2frac(values)
+                        end if
                         allocate(current)
                         current%data = values
                         current%next => head
                         head => current
                         nkpts = nkpts + 1
 
                     end do
+                else
+                    declared_nkpts = nk(1)
+                    cartesian_input = (mode(1:1)=='C' .or. mode(1:1)=='c')
+                    nullify(head)
+                    do while (.not. eof(unit))
+                        read(unit, '(A)', iostat=io) line
+                        if (io /= 0) exit
+                        if (len_trim(line)==0) cycle
+                        if (line(1:1)=='!' .or. line(1:1)=='#' .or. line(1:1)=='/') cycle
+                        read(line,*,iostat=io) values(1), values(2), values(3), weight
+                        if (io /= 0) then
+                            read(line,*,iostat=io) values
+                            if (io /= 0) cycle
+                            weight = 1.0_prec
+                        end if
+                        if (cartesian_input) then
+                            values = k2frac(values)
+                        end if
+                        allocate(current)
+                        current%data = values
+                        current%next => head
+                        head => current
+                        nkpts = nkpts + 1
+                        if (declared_nkpts > 0 .and. nkpts >= declared_nkpts) exit
+                    end do
+                    if (declared_nkpts > 0 .and. nkpts /= declared_nkpts) then
+                        write(*,'(A,1X,I0,1X,A,1X,I0)') '[IO] Warning: declared k-point count is', declared_nkpts, 'but read', nkpts
+                    end if
                 end if
             end if 
         close(unit)
@@ -176,6 +218,8 @@ subroutine readKPOINTS(filename,klist_frac)
             end do
         end if
 
+        if (present(is_path_mode)) is_path_mode = path_mode
+
     end subroutine readKPOINTS
 
     subroutine readWannier(seedname)
@@ -397,4 +441,123 @@ subroutine writeDOS(energy_grid, dos, filename)
         write(*, '(A,1X,F8.3,A)') "[IO] Done! Time elapsed (s): ", t_end - t_start
     end subroutine writeDOS
 
+    subroutine writeWavefunc(wavef, filename, istart, iend)
+        use constants, only: prec, f_eigvec, iband, nbands, nions
+        complex(prec), intent(in) :: wavef(:,:,:)
+        character(len=64), intent(in), optional :: filename
+        integer, intent(in), optional :: istart, iend
+        character(len=64) :: fout
+        integer :: unit, nkpts, ib_start, ib_end
+        integer :: ik, ib, i
+        real :: t_start, t_end
+
+        call cpu_time(t_start)
+
+        if (present(filename)) then
+            fout = filename
+        else
+            fout = f_eigvec
+        end if
+
+        nkpts = size(wavef, 3)
+        if (present(istart)) then
+            ib_start = istart
+        else
+            ib_start = iband(1)
+        end if
+        if (present(iend)) then
+            ib_end = iend
+        else
+            ib_end = iband(2)
+        end if
+
+        write(*,'(A,1X,A)') "[IO] Writing wavefunctions to file:", trim(fout)
+        unit = 600
+        open(unit, file=fout, status='replace', action='write')
+            write(unit,'(A,I8,A,I8,A,I8)') '# Wavefunctions: nkpts=', nkpts, ' nions=', nions, ' bands=', ib_start, ib_end
+            do ik = 1, nkpts
+                do ib = ib_start, ib_end
+                    write(unit,'(A,I6,1X,I6)') '# kpoint', ik, ib
+                    do i = 1, nions
+                        write(unit,'(2F18.12)') real(wavef(i, ib-ib_start+1, ik)), aimag(wavef(i, ib-ib_start+1, ik))
+                    end do
+                    write(unit,*)
+                end do
+            end do
+        close(unit)
+
+        call cpu_time(t_end)
+        write(*, '(A,1X,F8.3)') "[IO] Done! Time elapsed (s): ", t_end - t_start
+
+    end subroutine writeWavefunc
+
+    subroutine writeHamiltonianList(filename)
+        use constants, only: prec, nions, position_frac, basis, onsite, r_c, f_hr, timer_cpu
+        use tbmodel, only: tbg_hopping
+        character(len=64), intent(in), optional :: filename
+        character(len=64) :: fout
+        integer :: unit, i, j, nx, ny, nz
+        integer :: nxmax, nymax, nzmax
+        real(prec) :: a_len(3), R_frac(3), R_cart(3)
+        real(prec) :: tol
+        real :: t_start, t_end
+        complex(prec) :: hij
+
+        tol = 1.0e-12_prec
+        if (present(filename)) then
+            fout = filename
+        else
+            fout = f_hr
+        end if
+
+        if (timer_cpu) then
+            call cpu_time(t_start)
+            write(*,'(A,1X,A,1X,$)') '[IO] Writing real-space Hamiltonian list to', trim(fout)
+        end if
+
+        a_len(1) = norm2(basis(1,:))
+        a_len(2) = norm2(basis(2,:))
+        a_len(3) = norm2(basis(3,:))
+        nxmax = ceiling(r_c / max(a_len(1), tol)) + 1
+        nymax = ceiling(r_c / max(a_len(2), tol)) + 1
+        nzmax = ceiling(r_c / max(a_len(3), tol)) + 1
+
+        unit = 700
+        open(unit, file=fout, status='replace', action='write')
+            write(unit,'(A)') '# TB-tbG real-space Hamiltonian list'
+            write(unit,'(A,I0)') '# nions = ', nions
+            write(unit,'(A,F20.12)') '# cutoff_r_c = ', r_c
+            write(unit,'(A)') '# columns: i j R1_frac R2_frac R3_frac Re(H_ij) Im(H_ij)'
+
+            do i = 1, nions
+                do j = 1, nions
+                    do nx = -nxmax, nxmax
+                        do ny = -nymax, nymax
+                            do nz = -nzmax, nzmax
+                                R_frac = position_frac(j,:) - position_frac(i,:) + &
+                                         [real(nx,prec), real(ny,prec), real(nz,prec)]
+                                R_cart = matmul(transpose(basis), R_frac)
+                                if (norm2(R_cart) <= r_c + tol) then
+                                    if (i == j .and. nx == 0 .and. ny == 0 .and. nz == 0) then
+                                        hij = cmplx(onsite, 0.0_prec, kind=prec)
+                                    else if (norm2(R_cart) > tol) then
+                                        hij = cmplx(tbg_hopping(R_cart), 0.0_prec, kind=prec)
+                                    else
+                                        cycle
+                                    end if
+                                    write(unit,'(2I8,3F20.12,2F20.12)') i, j, R_frac, real(hij), aimag(hij)
+                                end if
+                            end do
+                        end do
+                    end do
+                end do
+            end do
+        close(unit)
+
+        if (timer_cpu) then
+            call cpu_time(t_end)
+            write(*,'(A,1X,F8.3)') 'Done! Time elapsed (s):', t_end - t_start
+        end if
+    end subroutine writeHamiltonianList
+
 end module ioutils