From c252edab9cc5157cbd75d3701523753e82545f7f Mon Sep 17 00:00:00 2001 From: "tao.zj" Date: Mon, 8 Jun 2026 16:42:02 +0800 Subject: [PATCH 1/3] add outputting eigenvectors function --- README.md | 15 ++++++++++++++ src/constants.f90 | 4 ++++ src/ioutils.f90 | 50 ++++++++++++++++++++++++++++++++++++++++++++++ src/mpi_solver.f90 | 19 ++++++++++++++---- src/parser.f90 | 4 ++++ 5 files changed, 88 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index b6f62d9..d96adfc 100644 --- a/README.md +++ b/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:
`B` – Band structure
`D` – Density of States (DOS)
`E` – Electron Energy Loss Spectroscopy (EELS) | | **FPOSCAR** | `` | Path to the VASP‑format POSCAR file (default: `cont.vasp`) | | **FEIGEN** | `` | Output file for band eigenvalues (default: `tb_band.dat`). Also enables band‑structure writing. | +| **EIGNVEC** | `` | Whether to write eigenvectors / wavefunctions for band calculations (default: `.false.`). When enabled, wavefunctions are saved together with the selected band range. | +| **FEIGENVEC** | `` | Output file for eigenvector / wavefunction data (default: `tb_wavef.dat`). | | **FDOS** | `` | Output file for DOS data (default: `tb_dos.dat`). | | **FKPOINTS** | `` | 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** | `` | 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 ` +- 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 diff --git a/src/constants.f90 b/src/constants.f90 index dd545c7..74a66b3 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -35,11 +35,13 @@ 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_matrix @@ -64,6 +66,7 @@ subroutine init_constants() 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_eels='tb_eels.dat' @@ -74,6 +77,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 diff --git a/src/ioutils.f90 b/src/ioutils.f90 index 56aab4b..87ac099 100644 --- a/src/ioutils.f90 +++ b/src/ioutils.f90 @@ -397,4 +397,54 @@ 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 + end module ioutils diff --git a/src/mpi_solver.f90 b/src/mpi_solver.f90 index 5a8623a..4e981b9 100644 --- a/src/mpi_solver.f90 +++ b/src/mpi_solver.f90 @@ -26,7 +26,7 @@ subroutine init_mpi() end subroutine init_mpi subroutine init_POSCAR_mpi(filename) - use constants, only: rank,nions,nbands,iband,basis,basis_rec,position_frac,position_cart,f_poscar,timer_mpi + use constants, only: rank,nions,nbands,iband,basis,basis_rec,position_frac,position_cart,f_poscar,timer_mpi,eigenvec_out,f_eigvec use ioutils, only: parsePOSCAR integer :: ierr character(len=*), intent(in), optional :: filename @@ -60,6 +60,8 @@ subroutine init_POSCAR_mpi(filename) end if call MPI_Bcast(position_frac, size(position_frac), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) call MPI_Bcast(position_cart, size(position_cart), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + call MPI_Bcast(eigenvec_out, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) + call MPI_Bcast(f_eigvec, len(f_eigvec), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) if (rank == 0 .and. timer_mpi) then call cpu_time(t_end) @@ -274,14 +276,15 @@ subroutine calculate_klist_mpi(klist,eig,wavef,tag) end subroutine calculate_klist_mpi subroutine calculate_band_mpi(kpath, nk) - use constants, only: prec,kpath_default,rank,kpt_select,f_kpoint,nk_path - use ioutils, only: readKPOINTS,writeBand,writeLabels + use constants, only: prec,kpath_default,rank,kpt_select,f_kpoint,nk_path,eigenvec_out,f_eigvec,nions,nbands,iband + use ioutils, only: readKPOINTS,writeBand,writeLabels,writeWavefunc real(prec), allocatable, optional:: kpath(:,:,:) integer, optional :: nk real(prec), allocatable:: paths(:,:,:),klist_frac(:,:),klist_cart(:,:),xlist(:),lable_positions(:),eig(:,:) real(prec), allocatable:: klist_tmp(:,:) character(len=3), allocatable :: labels(:) integer :: ierr,nkpts,nk_local,i,npath + complex(prec), allocatable :: wavef(:,:,:) if (rank == 0) then if (present(kpath)) then @@ -313,10 +316,18 @@ subroutine calculate_band_mpi(kpath, nk) if ( .not. allocated(klist_frac)) allocate(klist_frac(nkpts,3)) call MPI_Bcast(klist_frac, nkpts*3, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - call calculate_klist_mpi(klist_frac,eig) + if (eigenvec_out) then + call calculate_klist_mpi(klist_frac, eig, wavef) + else + call calculate_klist_mpi(klist_frac, eig) + end if if (rank == 0) then call writeBand(xlist, eig) + if (eigenvec_out) then + call writeWavefunc(wavef) + if (allocated(wavef)) deallocate(wavef) + end if call writeLabels(labels, lable_positions) end if end subroutine calculate_band_mpi diff --git a/src/parser.f90 b/src/parser.f90 index 2beb0fb..0462b1e 100644 --- a/src/parser.f90 +++ b/src/parser.f90 @@ -93,6 +93,10 @@ subroutine input_parser(filename) read(line,*,iostat=io) tmp,timer_mpi case('TDEBUG') read(line,*,iostat=io) tmp,timer_debug + case('EIGNVEC') + read(line,*,iostat=io) tmp,eigenvec_out + case('FEIGENVEC') + read(line,*,iostat=io) tmp,f_eigvec case('MODEL') read(line,*,iostat=io) tmp,model_type case('WRITEM') From a6c671faa2acc04df6ecaf200c78baaa1da4571a Mon Sep 17 00:00:00 2001 From: "tao.zj" Date: Wed, 10 Jun 2026 13:36:13 +0800 Subject: [PATCH 2/3] =?UTF-8?q?=E5=A2=9E=E5=BC=BA=20readKPOINTS=20?= =?UTF-8?q?=E5=AD=90=E4=BE=8B=E7=A8=8B=E4=BB=A5=E6=94=AF=E6=8C=81=E8=B7=AF?= =?UTF-8?q?=E5=BE=84=E6=A8=A1=E5=BC=8F=EF=BC=8C=E5=B9=B6=E5=9C=A8=20calcul?= =?UTF-8?q?ate=5Fband=5Fmpi=20=E4=B8=AD=E7=9B=B8=E5=BA=94=E8=B0=83?= =?UTF-8?q?=E6=95=B4=20KPOINTS=20=E5=A4=84=E7=90=86=E9=80=BB=E8=BE=91?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/ioutils.f90 | 54 +++++++++++++++++++++++++++++++++++++++++----- src/mpi_solver.f90 | 45 ++++++++++++++++++++++++++------------ 2 files changed, 80 insertions(+), 19 deletions(-) diff --git a/src/ioutils.f90 b/src/ioutils.f90 index 87ac099..31e09b3 100644 --- a/src/ioutils.f90 +++ b/src/ioutils.f90 @@ -105,7 +105,7 @@ 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) @@ -113,11 +113,13 @@ subroutine readKPOINTS(filename,klist_frac) 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,14 +148,22 @@ 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 @@ -157,6 +171,34 @@ subroutine readKPOINTS(filename,klist_frac) 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) diff --git a/src/mpi_solver.f90 b/src/mpi_solver.f90 index 4e981b9..7a66b1a 100644 --- a/src/mpi_solver.f90 +++ b/src/mpi_solver.f90 @@ -285,30 +285,45 @@ subroutine calculate_band_mpi(kpath, nk) character(len=3), allocatable :: labels(:) integer :: ierr,nkpts,nk_local,i,npath complex(prec), allocatable :: wavef(:,:,:) + logical :: path_mode + path_mode = .true. if (rank == 0) then if (present(kpath)) then paths = kpath else if (kpt_select) then - call readKPOINTS(f_kpoint,klist_tmp) - npath = size(klist_tmp,1)/2 - allocate(paths(npath,2,3)) - do i = 1, npath - paths(i,1,:) = klist_tmp(2*i-1,:) - paths(i,2,:) = klist_tmp(2*i,:) - end do + call readKPOINTS(f_kpoint,klist_tmp,path_mode) + if (path_mode) then + npath = size(klist_tmp,1)/2 + allocate(paths(npath,2,3)) + do i = 1, npath + paths(i,1,:) = klist_tmp(2*i-1,:) + paths(i,2,:) = klist_tmp(2*i,:) + end do + else + allocate(klist_frac(size(klist_tmp,1),3)) + klist_frac = klist_tmp + nkpts = size(klist_frac,1) + allocate(xlist(nkpts)) + do i = 1, nkpts + xlist(i) = real(i-1, prec) + end do + write(*,'(A,I0,A)') '[Main] Using explicit KPOINTS list with ', nkpts, ' points.' + end if else call kpath_default(paths) end if end if - if (present(nk)) then - nk_local = nk - else - nk_local = nk_path + if (path_mode .or. .not. allocated(klist_frac)) then + if (present(nk)) then + nk_local = nk + else + nk_local = nk_path + end if + call generate_kpath(paths, nk_local, klist_frac, klist_cart, xlist, labels, lable_positions) + nkpts = size(klist_frac,1) end if - call generate_kpath(paths, nk_local, klist_frac, klist_cart, xlist, labels, lable_positions) - nkpts = size(klist_frac,1) end if @@ -328,7 +343,9 @@ subroutine calculate_band_mpi(kpath, nk) call writeWavefunc(wavef) if (allocated(wavef)) deallocate(wavef) end if - call writeLabels(labels, lable_positions) + if (allocated(labels) .and. allocated(lable_positions)) then + call writeLabels(labels, lable_positions) + end if end if end subroutine calculate_band_mpi From df02ee0d2208fc19db3d97cfef7c79d33a19b26c Mon Sep 17 00:00:00 2001 From: "tao.zj" Date: Sat, 13 Jun 2026 19:38:58 +0800 Subject: [PATCH 3/3] adding outputing Hamiltonian --- src/constants.f90 | 3 + src/ioutils.f90 | 69 ++++++ src/mpi_solver.f90 | 8 +- src/parser.f90 | 4 + tools/README.md | 53 +++++ tools/__pycache__/compute_g.cpython-312.pyc | Bin 0 -> 12824 bytes tools/compute_g.py | 235 ++++++++++++++++++++ 7 files changed, 370 insertions(+), 2 deletions(-) create mode 100644 tools/README.md create mode 100644 tools/__pycache__/compute_g.cpython-312.pyc create mode 100644 tools/compute_g.py diff --git a/src/constants.f90 b/src/constants.f90 index 74a66b3..c01bafa 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -44,6 +44,7 @@ module constants 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 @@ -62,6 +63,7 @@ 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' @@ -69,6 +71,7 @@ subroutine init_constants() 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 diff --git a/src/ioutils.f90 b/src/ioutils.f90 index 31e09b3..b5f431b 100644 --- a/src/ioutils.f90 +++ b/src/ioutils.f90 @@ -491,4 +491,73 @@ subroutine writeWavefunc(wavef, filename, istart, iend) 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 diff --git a/src/mpi_solver.f90 b/src/mpi_solver.f90 index 7a66b1a..dd18c6d 100644 --- a/src/mpi_solver.f90 +++ b/src/mpi_solver.f90 @@ -26,8 +26,8 @@ subroutine init_mpi() end subroutine init_mpi subroutine init_POSCAR_mpi(filename) - use constants, only: rank,nions,nbands,iband,basis,basis_rec,position_frac,position_cart,f_poscar,timer_mpi,eigenvec_out,f_eigvec - use ioutils, only: parsePOSCAR + use constants, only: rank,nions,nbands,iband,basis,basis_rec,position_frac,position_cart,f_poscar,timer_mpi,eigenvec_out,f_eigvec,write_hr,f_hr + use ioutils, only: parsePOSCAR, writeHamiltonianList integer :: ierr character(len=*), intent(in), optional :: filename real :: t_start,t_end @@ -47,6 +47,10 @@ subroutine init_POSCAR_mpi(filename) write(*, '(4X,A,1X,$)') "[MPI] Updating POSCAR data to all MPI processes ..." end if + if (rank == 0 .and. write_hr) then + call writeHamiltonianList(f_hr) + end if + call MPI_Bcast(nions, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) call MPI_Bcast(nbands, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) call MPI_Bcast(iband, 2, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) diff --git a/src/parser.f90 b/src/parser.f90 index 0462b1e..d7ef43e 100644 --- a/src/parser.f90 +++ b/src/parser.f90 @@ -43,6 +43,10 @@ subroutine input_parser(filename) job_band = .true. write_band = .true. read(line,*,iostat=io) tmp,f_eig + case ('WHR') + read(line,*,iostat=io) tmp,write_hr + case ('FHR') + read(line,*,iostat=io) tmp,f_hr case ('FKPOINTS') kpt_select = .true. read(line,*,iostat=io) tmp,f_kpoint diff --git a/tools/README.md b/tools/README.md new file mode 100644 index 0000000..7a75d1b --- /dev/null +++ b/tools/README.md @@ -0,0 +1,53 @@ +# g-matrix post-processing + +This folder contains a small helper for computing + +$$g_{mn} = \langle \phi^f_{m,k+q} | (H^f - H^r) | \phi^r_{n,k} \rangle$$ + +from real-space hopping lists and wavefunction outputs. + +## 1. Export real-space Hamiltonian lists from TB-tbG + +In the input file, enable: + +```text +WHR T +FHR tb_hr.dat +``` + +Run once for the relaxed structure and once for the frozen-phonon structure. + +## 2. Compute g + +```bash +python tools/compute_g.py \ + --hr relaxed/tb_hr.dat \ + --hf frozen/tb_hr.dat \ + --wf-r relaxed/tb_wavef.dat \ + --wf-f frozen/tb_wavef.dat \ + --kpoints KPOINTS \ + --k-index 1 \ + --out g_matrix.dat +``` + +If you already know the k-vector, you can skip `--kpoints` and use: + +```bash +python tools/compute_g.py \ + --hr relaxed/tb_hr.dat \ + --hf frozen/tb_hr.dat \ + --wf-r relaxed/tb_wavef.dat \ + --wf-f frozen/tb_wavef.dat \ + --kvec 0.0 0.0 0.0 \ + --out g_matrix.dat +``` + +## Output + +`g_matrix.dat` contains one row per band pair: + +- frozen-phonon band index +- relaxed-state band index +- Re(g) +- Im(g) +- |g| diff --git a/tools/__pycache__/compute_g.cpython-312.pyc b/tools/__pycache__/compute_g.cpython-312.pyc new file mode 100644 index 0000000000000000000000000000000000000000..8c6d8276576f49a3ab3a0d1b495242e488af5aba GIT binary patch literal 12824 zcma)idu$s=nrAoJB%4o>67{yNmgI*RTaqQqjv~L-vLrvF#IHP@Nf?UVlqgc9T-}N! zZgZn#=P=a7I8({&iJt6UXvBaeSc}o%ZZR564x8E8!DV(2IBKYDn2rv};&Qmle~KL- zH?#N0eO1jSDOvWU2==4utE#W6tG@dEzN&w=+btA?U;M|XMvpdA)c?dUTF~Vm9!B6{ zp5iH<4pZ0FJAGXT@4B#ljJeLx6tDjwb=|--a5eHa-ta^Ex`{W!8_S#EYUb@c3vm|S z3~yH60#}=0=dGa@ibAZ85^VSLmC@8RJxzNYU*q>d53MlI`=k*+olnL$tg#;c58FV) zd1{!7(^4T_1}%9VuU}&BGgz3S4pVWRtm~s*%~4cbFY8Bf-Yr3^9j7A~$QZT=z$z=~0h$#y?K8ser+W!Xz}@QzJsWu<8`=+PJy7iVR9 zv_dPlc$W>bQ8vkJpwL@1byHy|GsQdcPphA%WLk7cm0C{TB~$NHf1}g-?`>M9uTafY z2PM*zI+}W!qTzS5eoFr`HA#EQV?Vhp_;?NpzHqBJ?(+-WS>ISNEJY(hUxXWpj*kZ; z!(2EhN^KSk2cdszT<}W*&-tU_iLr>-Lt;5D$c=KBJG=qp^K+Lw^Y=sfyTF|rBXSpm z5kcfc$%iCd?qpCJ;hH%g61d=SB#H#Sjfg2U`lDe&3bb}BMm1qkr?5iAALSvNViXXH zB2l56zs7?O;cxLZ9*kd;Cp@6P#-n|aRt&r}JuZBW$L}BDUlKM>er<%;!Pt($XrZA( zL_$H(Fc6)H@LVv$d8`T(9Ty^sL4=fX#X@u!vEmRl9u7*1J`BnFU_?@Mk#WV~L&!I+ z7y{v_Pg0D&amdUomKS{C385dW?J+CNILu*%;eC=%p#usNhF(z&GrF=tgp>3z#B zpPG-Z-1x#>ndr+^RL}R!^}Kg1(f4^pb*k%wlebSU99pYy&pB&y?wa}YbLZc?nCM?G zD^DJNCzd+(d3KKd>pf{h+INS2zkTWOa^>=kpV!75hd%Ab2a$#cm(&E(phSlBOtL5E^ z(>a^-)`i&%DOc)L>PEUb<@nUL@0)6>vg$iMWjXdu9c6Q_J8ClZPp>e`Z)9puWF03n z?8%(9Z05pu4Iq&*AN`|fgPwTz^sYXu{vWM|KDTL;nZ*N&nAyGXFi#atEKmQClBqsQ zGj{j$X7LK;(G6}Y0)N#*@PBv$E*GH7=ON>(Fn50dg~#c*F0PO3iv`QHtkdTCf1u}S zn77poHBJ4kev)3K2Vr#U$9z*BM$xq$L|7}d4=xJr6|t^(z#tEfFdrKhBEr--IuY9y zj06PajRd0+vF%tm>W2|_qOFkL2;yP`T>cOBe`lzSBX);kP6T>sH4qE;{Qtn6)0vFb5w( zH03W{hggUpoqZ6>gkI9kSWEh)xH{;hLatHV=q*8pl`f!@^0u0yM)6qR5|(JKk72cg z&-m&O^527r8{$T=@5LIhR$BsJS~|fFn+33)FDU_!uq^?;vQg6jOKZaiRt7z+TlyH3 zDjU>$KsNy55;p{vI0BMBzXW5%7lxI=9~?(`pn%1XBCx+cNe~afOGE?<0tQpH z;?jsfO!q((NeGsX@e!XWV5v7i1d8Hd0Ja3V&i_8ch5q3emr0c#)!HMxPEW(XH2fLD zYrDrp1`pZ+&x#=?ps1)=uu%&OQyrb%iXKW;bmKuq&qpPX8R6Nb7}W+T1~SYQ-AJfp zb{G+<8UfOR5{Z_ZZLB(H*rADF7>s2IodUt#@E6~K3z$e7#oA}$>(2UAU)I@}Fsz%c zw`yiNcz{E<@3flV-o#!?+5+2`{!bJ zd+!b08CW>7d?dTOXSKQ~>pYq;w3qu=T3X|qF7)GTQ>Tw|O$zHR_eCtX5Lar$IIWApxN2Dk+;gmHz>rFKt!W#94&ytVn*Jo z>BZYH>?(y~aRBTQn!&)?#M|SnRHo&V^;)jDS!QK(2zon7qb``31c;6n0H&4$rf>f7 zu|4F_%7l{i&_hpZp`}C@i?O2hEgh zA1wgfVoyue+OurrHv{kTO=-5$)@#69(+&W%1KJN`rlkG!&9`LJc6G`&2^U7fXr1;{ zyvuf}Mtd8@Mr%u0s?bJupa9W91N;TB2CO#9kS3LW0`Lp8 z006#e6ffCZ0$=q*YN@)Um8L%Gmh+a&uO9}h@hZ&AAHZmiJL4|dIj6g6W~g{Mek=JB zFdJ}V`Az*7uaKFys$@p42t5s>&?Pl$74tPR1Nzm<6+l1=`q5NLk2Pts#k-94a@}%o z(aRNg%WhRKzdr7cyR??XE9G*zGN5~u)FgE`FRQesnUc!28g^62T-w7g!wkOpzbC0} zCBYJqbc>em2g(HXyW-AC=r!^Tx=d1@y4b%y??bS{MI!OGtT#Na28RJ!iArm4r3$HJOF=sRPB5v74gg^wb2U#Ir0o-|9p*eIMi(?B+ zSXqp1=0c=%5GKzEvvh>*D{|utwx2OJiCm0J#zc^^@ck_mU4%#K3t`?A#G?OJ}H*FHQfH zC+jQ}tq%#)G50HzJ`{PC8v!OIAdYn9@mSdbj+D>mBRu7(3QADSm~dfoCn9VW6cf<2 zFmP;yHAnS0%0QZquqi>-*SvMk3N+F{7#)BhyKAQZ3uk450Ty@bqBE>puFK`_!uYRMveu(Vr`;O6^%KYgn(TPxYtw+}nR=|H5##>Bx#} zwW()y$I)!fv83U1n`^^PxvDmtl*2hQkh554POsa_l7@HA&-CXkc6ivID|08W%*H+| zYk<BfocO-i^y=<_QyDC-v z?(y8thU7pFH1ypzB~Rt*>Qll`J-MAt>5khkeel}t*RnfXQ@YO^8gn}ubM-qvIDY&1 zLd9BrD}Lpg?y+~+4=pLfdi}1{#QQx9uAiRV(CK&X$#Hw{+3(nwgysH~#-E?PAN^SO zkJmmqoVj%Qciq47uKwB8)wZi^+%@Rm+Pc)G+r|{VQA^eDND1$IHtML#`uS()o=tb& zMfaxeOf4Hf8v6UOUyiNqyO?QucFjGQW37N**@S-MC}njd4LA3HbAqy!1M^OGr_ZFH zSr}aDTG^YaJ(YFzX4u{@8Ox6@ymeu2aOT1q)BJgTOf{)%>Co?YcC54ZjB|JTTKdew$(5rktr^GJHTK+B z7Tc}%+4jX3*6O+zU$}Dw_^bb?hZbJ=eO*`9a(L!+&gx7yCwekW<##t|p#Q7pqyP5p zL8_+d0i|>8%(<&m9dl=srxo{(4Fg2{+qWFGyX^r*@7$wsyEpW3|9=nNl)d>64=q&v zVeu#o-K6VWhyFL6JI*b(r(mx%(ts54ha@WC>Vx_W&Bu z$8}p7fP@zW65cEzVGOo+r-P?GU|| zhv4i2VlNh4MEqi@Znpd2P{84B8nCYsYOwQmnTbJ69VY)TLl2q&Q9VIRU=*Um<_&L0dJ3Do4zH;#v;C2D!8Jrs0c|H*G)d+a0blPkasz|Vt=P1C!p?O zz#CGji{Zla&z&2*dSzROb0CVwe3IsG#tT~q1}|d(Sn?1D$8~a%iLoJAo1*~@m#ILD zgF++1Ar9Okd51Z062(Ff+rnAlkE@(xKl8qAV30L z$4{5x0#5Fb!i+_EL19AUl86ehVmmiVmAc@o5K#1LZ9&6O(iHr~9+-}>-j{D(??2zW zFY7#V$CuVEbpEV+vHNd(lI8O?b2Wcem-ZzL8Rrp#0=_a?5{LgXmOSnDZe9L zCOhIzxh(FAFtQDv?7&2eR@S3Y(dFuI*%w=%{+4_fgxUL`wai(D1J?peK0fDS&bZq0DR99T6JoXT^< z-f4~(0=|i`#I;NSWdWycOo({HLO?n`0N#4R!wp4Y*K7#+ILrn^g2aIw7}2nAO5j_8 z3`nXcmhjI4OCuSMV{+(4EusB4JO@bDLX4-4BPJ;j4TD(_F{dQ<5z50nk_Uo`M5z@# zz+kY1Fjy%RjkJ!BL?htaRa=0gu{FjG0E3FXzABMC7vNss+Lo(mLU{cEp{}?oU<@T+ zNWiXVa8eUOOB()|_7J+uy*?T6243fU5*NC`wScghsvcG7tDb7bG=j+P7(4{j2L>ua3cd!=C7hED-B3v@I)-1% zwzN($;{6cFqTU0Fnb0OL@`e->`2_G_Da74S*nH#-EjHi{1FAyVybzXr-Vuc!7BQQAl!uWRm5iA&^Cz2wXa;0|0)KH8F2o!=8d8&4`@Wfh zb^EUT)3-*du4CD~JejRMxyn?{oJ&l8Ro}VXwNjIGKo)o^54T6rnk zFtEz*NLZ3*bCnGVORl1Rmi;pSlwhGOyEAE;x6j$ZI6Jv?%e}TcZOg#E4rOg!GyUt< zisV3QU}4Xi^&sp=ZhFf(%9B@TLo;W-WSzI{v-Wq&K4WXa)t)-D$~LEmSJ~E?Q#os8 zays3Zu{LLzX4O*;v#=1bFM?RXi;q}Ad|YBluUQ8;EWyakc6gdE^>7s}Rl!~Ymh|O) zY_N?~V#^ACKad#AOPI)Jk3lvC$mR}so42hK{G5jO8QusuQN5eDwS6WcV;(KF`}wBx zB(20ZY}V{MNn!t{;AI9}`8-42^phh8_>c>wSyMQY@(8}0&im6m=0Sy?Bq#}H9i}p| zgUx`_)T0927#YIe02fgo5CEN)VGo=Rz}HhlBHS#BuA4{wIIFSqU|UD0pwJeNQT3=8 zCQ(ol6te&zPe5Sl6b2t%P?+GDZ&=a$hD47g?*~z5S3oi4rx|(^id4Z2-YZpfV?%lO zMr4TW28nq4t~g`VNn#hl*X14dsum9Ydk|THzxW|sU?)@CDF4*_)PgU^ns1qAO&RCu zHMVcP66VdnRPLMYPn=4i_1dQN)BjStf6kDklbyM$YJe3TDQ0dssmGoT)0`>QnHoy% zOS{sY>7n$#1=qqM411ObGEW^{=~%h)_b>hOrH=!dW6yqa>A#7YD=+->ROY4EGK1cX z*O#do`qb(F#z55__(o4vwR~rwEc-t*KecLpO5M`F)l2u*>Hne5(A#YCgWiRJ8%(A_ z00wX=7#PK7b4%c9STys%AP)YwZJqvlp5YBkMs50l<8^yE6UedOL(U9xmhT~F1v%UI zkhAlSrLy8cevEcbsjyGXyOzrTWNRvTw^W#HeC1M=cGl{#?XNDWq2_zkP|Mfx^-DYc zq$WH0U3|k5w@L0vR2JL^+d&ufcO}~Z+e#Ehf(}Y3oQoz1n>^O>5Wc5a8S;@Xtc$ zaQlnZ9^jyH{%qG>dnGm7WUc5hZBh_OC^PtE7NcwYOH@MD5)EI<2VyT6e{lH z&EQDg2!ekA;oU|Q#TTq9ftQLN+K-;YuR*xLaZ!3)VQ`%&fFm#A3uc`N7mt37FYn+B z#<$z~;27-QM0X+e*I42LU^_srtn)LrVU=yj9GF;R-vGwM7=ASP)?o4oc%s^`&kU|H zuRI~DJI!Z0E@P=HniPq;grhENQAaa-UdIu>Ld4xGnf7OK)SwpCo$gpVoOwEgBg0zc zvAh1v!PjtDi;u1%A)eGi~rSntf5FytJ>jj1sjCj4M6~>r{9`HC&3ZB3YPK=FDtLOJj z_}FteID|0gP<=7@WGL~!px5yeb_)|{1epnB%;KUWzA*ucK7&*7sW+9HvB$s*JrNd8 zphbwoC{e_-26pFY`U|S;3#$BoP&E&X)wJ_L7fp9OWNozL!LDAKu6l4@S4;1E*xg9) mesGpH(meg3d@tR#ahawVTO#&*eeLHA`|k!e^b}J^D*XSdUySVl literal 0 HcmV?d00001 diff --git a/tools/compute_g.py b/tools/compute_g.py new file mode 100644 index 0000000..27c5245 --- /dev/null +++ b/tools/compute_g.py @@ -0,0 +1,235 @@ +import argparse +import numpy as np +from pathlib import Path + + +def read_hopping_list(path): + """Read real-space Hamiltonian hopping list. + + Expected columns: + i j R1_frac R2_frac R3_frac Re Im + Lines starting with # are ignored. + """ + data = [] + with open(path, 'r', encoding='utf-8', errors='ignore') as f: + for line in f: + line = line.strip() + if not line or line.startswith('#'): + continue + parts = line.split() + if len(parts) < 7: + continue + i = int(parts[0]) + j = int(parts[1]) + r = np.array([float(parts[2]), float(parts[3]), float(parts[4])], dtype=float) + val = float(parts[5]) + 1j * float(parts[6]) + data.append((i, j, r, val)) + if not data: + raise ValueError(f'No hopping entries found in {path}') + return data + + +def infer_nions(hops): + return max(max(i, j) for i, j, _, _ in hops) + + +def build_hk(hops, k_frac, nions=None): + """Reconstruct H(k) from hopping list. + + k_frac must be in fractional reciprocal coordinates, consistent with the + Fortran phase convention exp(i 2π k·R_frac). + """ + if nions is None: + nions = infer_nions(hops) + hk = np.zeros((nions, nions), dtype=np.complex128) + for i, j, r_frac, val in hops: + phase = np.exp(2j * np.pi * np.dot(k_frac, r_frac)) + hk[i - 1, j - 1] += val * phase + return hk + + +def parse_wavefunction_file(path, k_index): + """Parse one k-point block from tb_wavef.dat-like file. + + Returns + ------- + bands : list[int] + Global band indices stored in the file. + psi : np.ndarray + Shape (nions, nbands) complex coefficients. + """ + with open(path, 'r', encoding='utf-8', errors='ignore') as f: + lines = f.readlines() + + nions = None + for line in lines: + if line.startswith('# Wavefunctions:'): + # format: # Wavefunctions: nkpts= ... nions= ... bands= istart iend + parts = line.replace('#', '').replace('=', ' ').split() + if 'nions' in parts: + nions = int(parts[parts.index('nions') + 1]) + break + + blocks = [] + idx = 0 + total = len(lines) + while idx < total: + line = lines[idx].strip() + if line.startswith('# kpoint'): + toks = line.split() + # '# kpoint ' + if len(toks) >= 4: + try: + ik = int(toks[2]) + ib = int(toks[3]) + except ValueError: + idx += 1 + continue + if ik == k_index: + coeffs = [] + j = idx + 1 + while j < total: + ln = lines[j].strip() + if not ln: + break + if ln.startswith('#'): + break + parts = ln.split() + if len(parts) >= 2: + try: + coeffs.append(float(parts[0]) + 1j * float(parts[1])) + except ValueError: + pass + j += 1 + blocks.append((ib, np.array(coeffs, dtype=np.complex128))) + idx += 1 + + if not blocks: + raise ValueError(f'k_index={k_index} not found in {path}') + + blocks.sort(key=lambda x: x[0]) + bands = [b for b, _ in blocks] + nions_found = len(blocks[0][1]) + if nions is not None and nions != nions_found: + print(f'[warn] header nions={nions}, parsed nions={nions_found}') + + psi = np.column_stack([vec for _, vec in blocks]) + return bands, psi + + +def parse_kpoints_explicit(path): + """Parse explicit reciprocal KPOINTS list. + + Expected format: + line1 comment + line2 number of k-points + line3 Reciprocal + remaining lines: kx ky kz [weight] + """ + with open(path, 'r', encoding='utf-8', errors='ignore') as f: + lines = [ln.strip() for ln in f if ln.strip()] + if len(lines) < 4: + raise ValueError('KPOINTS file too short') + nk = int(lines[1].split()[0]) + mode = lines[2].lower() + if not mode.startswith('reciprocal'): + raise ValueError('This helper currently expects Reciprocal explicit KPOINTS') + kpts = [] + for line in lines[3:3 + nk]: + parts = line.split() + if len(parts) < 3: + continue + kpts.append([float(parts[0]), float(parts[1]), float(parts[2])]) + if len(kpts) != nk: + raise ValueError(f'Expected {nk} k-points, parsed {len(kpts)}') + return np.array(kpts, dtype=float) + + +def compute_g(hr_list, hf_list, k_index, kq_index, k_frac,wf_r=None, wf_f=None, use_frozen_left=False): + """Compute g matrix. + + By default (use_frozen_left=False) both bra and ket use relaxed-state + wavefunctions: bra = psi_r(k+q), ket = psi_r(k). This follows the + standard first-order electron-phonon matrix element definition. + If `use_frozen_left=True`, the bra is taken from the frozen-state + wavefunction file `wf_f` at kq (legacy option). + """ + hops_r = read_hopping_list(hr_list) + hops_f = read_hopping_list(hf_list) + + # right state: relaxed at k + bands_r_k, psi_r_k = parse_wavefunction_file(wf_r, k_index) + # left state: by default relaxed at k+q + bands_r_kq, psi_r_kq = parse_wavefunction_file(wf_r, kq_index) + + if use_frozen_left: + # legacy: read frozen-state wavefunction for left vector + bands_f_kq, psi_f_kq = parse_wavefunction_file(wf_f, kq_index) + left_bands = bands_f_kq + psi_left = psi_f_kq + else: + left_bands = bands_r_kq + psi_left = psi_r_kq + + # ensure dimensions + nions = max(psi_r_k.shape[0], psi_left.shape[0]) + hr_k = build_hk(hops_r, k_frac, nions=nions) + hf_k = build_hk(hops_f, k_frac, nions=nions) + delta_h = hf_k - hr_k + + g = psi_left.conj().T @ delta_h @ psi_r_k + return left_bands, bands_r_k, g + + +def save_g_matrix(path, bands_f, bands_r, g): + with open(path, 'w', encoding='utf-8') as f: + f.write('# columns: m_band_f n_band_r Re(g) Im(g) Abs(g)\n') + for i, mb in enumerate(bands_f): + for j, nb in enumerate(bands_r): + val = g[i, j] + f.write(f'{mb:8d} {nb:8d} {val.real:20.12e} {val.imag:20.12e} {abs(val):20.12e}\n') + + +def main(): + p = argparse.ArgumentParser(description='Compute g_mn = from hopping lists and wavefunctions') + p.add_argument('--hr', required=True, help='Relaxed-state real-space Hamiltonian list') + p.add_argument('--hf', required=True, help='Frozen-phonon real-space Hamiltonian list') + p.add_argument('--wf-r', required=True, help='Relaxed-state wavefunction file') + p.add_argument('--wf-f', required=True, help='Frozen-state wavefunction file') + p.add_argument('--kpoints', help='Explicit reciprocal KPOINTS file used for the calculation') + p.add_argument('--k-index', type=int, help='1-based k-point index to use') + p.add_argument('--kvec', nargs=3, type=float, help='Fractional k vector if no KPOINTS file is provided') + p.add_argument('--kq-index', type=int, help='1-based k+q point index when using --kpoints') + p.add_argument('--use-frozen-left', action='store_true', help='Use frozen-state wavefunction for left vector (legacy)') + p.add_argument('--out', default='g_matrix.dat', help='Output matrix file') + args = p.parse_args() + + if args.kpoints: + kpts = parse_kpoints_explicit(args.kpoints) + if args.k_index is None: + raise SystemExit('Provide --k-index when using --kpoints') + k_index = args.k_index + if not (1 <= k_index <= len(kpts)): + raise SystemExit(f'k-index out of range: 1..{len(kpts)}') + k_frac = kpts[k_index - 1] + if args.kq_index is None: + kq_index = k_index + else: + kq_index = args.kq_index + if not (1 <= kq_index <= len(kpts)): + raise SystemExit(f'kq-index out of range: 1..{len(kpts)}') + else: + if args.kvec is None: + raise SystemExit('Provide either --kpoints + --k-index, or --kvec') + k_index = args.k_index if args.k_index is not None else 1 + k_frac = np.array(args.kvec, dtype=float) + # when user provides explicit kvec, use same index for k+q unless provided + kq_index = args.kq_index if args.kq_index is not None else k_index + + bands_f, bands_r, g = compute_g(args.hr, args.hf, k_index, kq_index, k_frac, wf_r=args.wf_r, wf_f=args.wf_f, use_frozen_left=args.use_frozen_left) + save_g_matrix(args.out, bands_f, bands_r, g) + print(f'Wrote {args.out} with shape {g.shape}') + + +if __name__ == '__main__': + main()