wannier90_nnkpts.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Read explicit Wannier90 overlap connections, including nonuniform closed loops.
! **************************************************************************************************
MODULE wannier90_nnkpts
   USE cp_files,                        ONLY: close_file,&
                                              open_file
   USE ieee_arithmetic,                 ONLY: ieee_is_finite
   USE kinds,                           ONLY: dp
   USE string_utilities,                ONLY: lowercase
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'wannier90_nnkpts'
   ! Absolute file-format tolerances: Angstrom and reciprocal Angstrom, respectively.
   REAL(KIND=dp), PARAMETER :: real_lattice_tol = 1.e-6_dp, recip_lattice_tol = 1.e-6_dp
   PUBLIC :: read_wannier90_nnkpts
CONTAINS
 
! **************************************************************************************************
!> \brief Read and validate the cell, points, and directed connections of an nnkp file.
!> \param filename input file
!> \param real_lattice direct lattice vectors as columns, in Angstrom
!> \param recip_lattice reciprocal lattice vectors as columns, in inverse Angstrom (including 2*pi)
!> \param kpt fractional k-points, without coordinate wrapping
!> \param nnlist target point for each directed connection
!> \param nncell reciprocal lattice translation for each connection
! **************************************************************************************************
   SUBROUTINE read_wannier90_nnkpts(filename, real_lattice, recip_lattice, kpt, nnlist, nncell)
      CHARACTER(len=*), INTENT(IN)                       :: filename
      REAL(kind=dp), INTENT(IN)                          :: real_lattice(3, 3), recip_lattice(3, 3)
      REAL(kind=dp), ALLOCATABLE, INTENT(OUT)            :: kpt(:, :)
      INTEGER, ALLOCATABLE, INTENT(OUT)                  :: nnlist(:, :), nncell(:, :, :)
 
      CHARACTER(len=1024)                                :: line
      CHARACTER(len=64)                                  :: block, marker
      INTEGER                                            :: i, ios, j, n, nexclude, nneigh, origin, &
                                                            shift(3), TARGET, unit
      INTEGER, ALLOCATABLE                               :: counts(:)
      LOGICAL                                            :: have_real, have_recip
      REAL(kind=dp)                                      :: lattice(3, 3)
 
      have_real = .false.
      have_recip = .false.
      CALL open_file(filename, unit_number=unit, file_status="OLD", file_action="READ")
      DO
         CALL next_line(unit, line, ios)
         IF (ios < 0) EXIT
         IF (ios /= 0) cpabort("NNKP_FILE: read error.")
         marker = ""
         block = ""
         READ (line, *, iostat=ios) marker, block
         IF (ios /= 0) cycle
         CALL lowercase(marker)
         CALL lowercase(block)
         IF (marker /= "begin") cycle
         SELECT CASE (trim(block))
         CASE ("real_lattice", "recip_lattice")
            DO i = 1, 3
               CALL required_line(unit, line)
               READ (line, *, iostat=ios) lattice(:, i)
               IF (ios /= 0) cpabort("NNKP_FILE: invalid lattice vector.")
            END DO
            IF (.NOT. all(ieee_is_finite(lattice))) cpabort("NNKP_FILE: nonfinite lattice.")
            IF (block == "real_lattice") THEN
               IF (have_real) cpabort("NNKP_FILE: duplicate real_lattice block.")
               IF (maxval(abs(lattice - real_lattice)) > real_lattice_tol) THEN
                  cpabort("NNKP_FILE: real lattice differs from the CP2K cell (Angstrom).")
               END IF
               have_real = .true.
            ELSE
               IF (have_recip) cpabort("NNKP_FILE: duplicate recip_lattice block.")
               IF (maxval(abs(lattice - recip_lattice)) > recip_lattice_tol) THEN
                  cpabort("NNKP_FILE: reciprocal lattice differs from the CP2K cell.")
               END IF
               have_recip = .true.
            END IF
            CALL end_block(unit, block)
         CASE ("kpoints")
            IF (ALLOCATED(kpt)) cpabort("NNKP_FILE: duplicate kpoints block.")
            CALL required_line(unit, line)
            READ (line, *, iostat=ios) n
            IF (ios /= 0) cpabort("NNKP_FILE: invalid point count.")
            IF (n < 1) cpabort("NNKP_FILE: point count must be positive.")
            ALLOCATE (kpt(3, n))
            DO i = 1, n
               CALL required_line(unit, line)
               READ (line, *, iostat=ios) kpt(:, i)
               IF (ios /= 0) cpabort("NNKP_FILE: invalid k-point.")
            END DO
            IF (.NOT. all(ieee_is_finite(kpt))) cpabort("NNKP_FILE: nonfinite k-point.")
            CALL end_block(unit, block)
         CASE ("nnkpts")
            IF (.NOT. ALLOCATED(kpt)) cpabort("NNKP_FILE: kpoints must precede nnkpts.")
            IF (ALLOCATED(nnlist)) cpabort("NNKP_FILE: duplicate nnkpts block.")
            CALL required_line(unit, line)
            READ (line, *, iostat=ios) nneigh
            IF (ios /= 0) cpabort("NNKP_FILE: invalid neighbour count.")
            IF (nneigh < 1) cpabort("NNKP_FILE: neighbour count must be positive.")
            n = SIZE(kpt, 2)
            ALLOCATE (nnlist(n, nneigh), nncell(3, n, nneigh), counts(n))
            counts = 0
            DO i = 1, n*nneigh
               CALL required_line(unit, line)
               READ (line, *, iostat=ios) origin, TARGET, shift
               IF (ios /= 0) cpabort("NNKP_FILE: invalid connection.")
               IF (min(origin, TARGET) < 1 .OR. max(origin, TARGET) > n) THEN
                  cpabort("NNKP_FILE: connection point index out of range.")
               END IF
               IF (counts(origin) == nneigh) cpabort("NNKP_FILE: too many neighbours for a point.")
               DO j = 1, counts(origin)
                  IF (nnlist(origin, j) == TARGET .AND. all(nncell(:, origin, j) == shift)) THEN
                     cpabort("NNKP_FILE: duplicate connection.")
                  END IF
               END DO
               counts(origin) = counts(origin) + 1
               nnlist(origin, counts(origin)) = TARGET
               nncell(:, origin, counts(origin)) = shift
            END DO
            IF (any(counts /= nneigh)) cpabort("NNKP_FILE: missing connections.")
            DEALLOCATE (counts)
            CALL end_block(unit, block)
         CASE ("exclude_bands")
            CALL required_line(unit, line)
            READ (line, *, iostat=ios) nexclude
            IF (ios /= 0) cpabort("NNKP_FILE: invalid exclude_bands count.")
            IF (nexclude /= 0) cpabort("NNKP_FILE: specify band exclusions using CP2K EXCLUDE_BANDS.")
            CALL end_block(unit, block)
         CASE DEFAULT
            ! Skip projections and other Wannier90 blocks not required for overlaps.
            DO
               CALL required_line(unit, line)
               marker = ""
               READ (line, *, iostat=ios) marker
               IF (ios /= 0) cpabort("NNKP_FILE: invalid line in unused block.")
               CALL lowercase(marker)
               IF (marker == "end") EXIT
            END DO
         END SELECT
      END DO
      CALL close_file(unit)
      IF (.NOT. have_real .OR. .NOT. have_recip) cpabort("NNKP_FILE: missing lattice blocks.")
      IF (.NOT. ALLOCATED(kpt) .OR. .NOT. ALLOCATED(nnlist)) THEN
         cpabort("NNKP_FILE: missing kpoints or nnkpts block.")
      END IF
   END SUBROUTINE read_wannier90_nnkpts
 
! **************************************************************************************************
!> \brief Read a nonempty line, removing comments.
!> \param unit input unit
!> \param line resulting line
!> \param ios IO status
! **************************************************************************************************
   SUBROUTINE next_line(unit, line, ios)
      INTEGER, INTENT(IN)                                :: unit
      CHARACTER(len=*), INTENT(OUT)                      :: line
      INTEGER, INTENT(OUT)                               :: ios
 
      INTEGER                                            :: i
 
      DO
         READ (unit, '(A)', iostat=ios) line
         IF (ios /= 0) RETURN
         i = scan(line, "!#")
         IF (i > 0) line(i:) = ""
         line = adjustl(line)
         IF (len_trim(line) > 0) RETURN
      END DO
   END SUBROUTINE next_line
 
! **************************************************************************************************
!> \brief Read a required nonempty line.
!> \param unit input unit
!> \param line resulting line
! **************************************************************************************************
   SUBROUTINE required_line(unit, line)
      INTEGER, INTENT(IN)                                :: unit
      CHARACTER(len=*), INTENT(OUT)                      :: line
 
      INTEGER                                            :: ios
 
      CALL next_line(unit, line, ios)
      IF (ios /= 0) cpabort("NNKP_FILE: unexpected end of file or read error.")
   END SUBROUTINE required_line
 
! **************************************************************************************************
!> \brief Check a block terminator.
!> \param unit input unit
!> \param block expected block name
! **************************************************************************************************
   SUBROUTINE end_block(unit, block)
      INTEGER, INTENT(IN)                                :: unit
      CHARACTER(len=*), INTENT(IN)                       :: block
 
      CHARACTER(len=1024)                                :: line
      CHARACTER(len=64)                                  :: marker, name
      INTEGER                                            :: ios
 
      CALL required_line(unit, line)
      READ (line, *, iostat=ios) marker, name
      IF (ios /= 0) cpabort("NNKP_FILE: missing block terminator.")
      CALL lowercase(marker)
      CALL lowercase(name)
      IF (marker /= "end" .OR. name /= block) cpabort("NNKP_FILE: incorrect block terminator.")
   END SUBROUTINE end_block
 
END MODULE wannier90_nnkpts