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..c01bafa 100644 --- a/src/constants.f90 +++ b/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 diff --git a/src/ioutils.f90 b/src/ioutils.f90 index 56aab4b..b5f431b 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) @@ -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 diff --git a/src/mpi_solver.f90 b/src/mpi_solver.f90 index 5a8623a..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 - 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) @@ -60,6 +64,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,38 +280,54 @@ 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(:,:,:) + 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 @@ -313,11 +335,21 @@ 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) - call writeLabels(labels, lable_positions) + if (eigenvec_out) then + call writeWavefunc(wavef) + if (allocated(wavef)) deallocate(wavef) + end if + if (allocated(labels) .and. allocated(lable_positions)) then + call writeLabels(labels, lable_positions) + end if end if end subroutine calculate_band_mpi diff --git a/src/parser.f90 b/src/parser.f90 index 2beb0fb..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 @@ -93,6 +97,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') 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 0000000..8c6d827 Binary files /dev/null and b/tools/__pycache__/compute_g.cpython-312.pyc differ 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()