wannier90_projections_unittest.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 Test mesh-wide rank validation and unitary covariance of AO trial selection.
!> \author Thomas D. Kuehne, tkuehne@cp2k.org
! **************************************************************************************************
PROGRAM wannier90_projections_unittest
#if defined (__HAS_IEEE_EXCEPTIONS)
   USE ieee_exceptions, ONLY: ieee_all, ieee_divide_by_zero, ieee_invalid, ieee_overflow, &
                              ieee_get_halting_mode, ieee_set_halting_mode, ieee_support_halting
#endif
   USE kinds, ONLY: dp
   USE message_passing, ONLY: mp_comm_type, &
                              mp_para_env_release, &
                              mp_para_env_type, &
                              mp_world_finalize, &
                              mp_world_init
   USE wannier90_projections, ONLY: wannier90_hybrid_candidates, &
                                    wannier90_window_mask, &
                                    wannier90_projection_rank, &
                                    wannier90_select_projections
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   COMPLEX(KIND=dp), ALLOCATABLE :: candidates(:, :, :)
   COMPLEX(KIND=dp), ALLOCATABLE :: hybrids(:, :, :), rotated_hybrids(:, :, :)
   COMPLEX(KIND=dp) :: quartet(2, 4)
   COMPLEX(KIND=dp) :: amn(2, 2, 2), reference(2, 2, 2), rotation(2, 2, 2)
   COMPLEX(KIND=dp) :: repaired(2, 2, 3)
   COMPLEX(KIND=dp) :: mixed_reference(4, 3, 2)
   LOGICAL :: allowed(4)
   INTEGER :: anchor, ik, nfail, pivots(2), reference_pivots(2)
   INTEGER :: empty_groups(4, 0)
   REAL(kind=dp) :: angle, min_svalue
   TYPE(mp_comm_type) :: comm
   TYPE(mp_para_env_type), POINTER :: para_env
#if defined (__HAS_IEEE_EXCEPTIONS)
   LOGICAL :: saved_halt(size(ieee_all)), expected_halt(size(ieee_all)), actual_halt(size(ieee_all))
#endif
 
   CALL mp_world_init(comm)
   ALLOCATE (para_env)
   CALL para_env%from_dup(comm)
   nfail = 0
#if defined (__HAS_IEEE_EXCEPTIONS)
   CALL ieee_get_halting_mode(ieee_all, saved_halt)
   IF (ieee_support_halting(ieee_invalid)) CALL ieee_set_halting_mode(ieee_invalid, .true.)
   IF (ieee_support_halting(ieee_divide_by_zero)) CALL ieee_set_halting_mode(ieee_divide_by_zero, .true.)
   IF (ieee_support_halting(ieee_overflow)) CALL ieee_set_halting_mode(ieee_overflow, .true.)
   CALL ieee_get_halting_mode(ieee_all, expected_halt)
#endif
   IF (para_env%is_source()) THEN
      ALLOCATE (candidates(2, 3, 2))
      candidates = cmplx(0.0_dp, 0.0_dp, kind=dp)
      ! The strongest stacked column vanishes at the second k-point.
      candidates(1, 1, 1) = 5.0_dp
      candidates(2, 2, :) = 2.0_dp
      candidates(1, 3, :) = 0.4_dp
   ELSE
      ALLOCATE (candidates(0, 0, 0))
   END IF
   CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue, anchor_kpoint=0)
   IF (min_svalue > 1.0e-14_dp) nfail = nfail + 1
   CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue)
   IF (abs(min_svalue - 0.4_dp) > 1.0e-14_dp) nfail = nfail + 1
   IF (.NOT. all(pivots == [2, 3])) nfail = nfail + 1
   reference = amn
   reference_pivots = pivots
 
   DO ik = 1, 2
      angle = 0.31_dp*real(ik, kind=dp)
      rotation(1, 1, ik) = cos(angle)
      rotation(2, 2, ik) = cos(angle)
      rotation(1, 2, ik) = cmplx(0.0_dp, sin(angle), kind=dp)
      rotation(2, 1, ik) = cmplx(0.0_dp, sin(angle), kind=dp)
      IF (para_env%is_source()) candidates(:, :, ik) = matmul(rotation(:, :, ik), candidates(:, :, ik))
   END DO
   CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue)
   IF (abs(min_svalue - 0.4_dp) > 1.0e-14_dp) nfail = nfail + 1
   IF (.NOT. all(pivots == reference_pivots)) nfail = nfail + 1
   DO ik = 1, 2
      IF (maxval(abs(amn(:, :, ik) - matmul(rotation(:, :, ik), reference(:, :, ik)))) > 1.0e-14_dp) THEN
         nfail = nfail + 1
      END IF
   END DO
   IF (para_env%is_source()) candidates(:, :, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
   CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue)
   IF (min_svalue /= 0.0_dp) nfail = nfail + 1
   DEALLOCATE (candidates)
   IF (para_env%is_source()) THEN
      ALLOCATE (candidates(2, 5, 2))
      DO ik = 1, 2
         DO anchor = 1, 5
            candidates(1, anchor, ik) = cmplx(0.1_dp*anchor, 0.2_dp*ik, kind=dp)
            candidates(2, anchor, ik) = cmplx(-0.3_dp*anchor, 0.05_dp*ik*anchor, kind=dp)
         END DO
      END DO
   ELSE
      ALLOCATE (candidates(0, 0, 0))
   END IF
   CALL wannier90_hybrid_candidates(candidates, reshape([1, 2, 3, 4], [4, 1]), hybrids)
   IF (para_env%is_source()) THEN
      IF (.NOT. all(shape(hybrids) == [2, 5, 2])) nfail = nfail + 1
      IF (maxval(abs(hybrids(:, 1, :) - candidates(:, 5, :))) > 1.0e-14_dp) nfail = nfail + 1
      DO ik = 1, 2
         ! A full tetrahedral quartet is a unitary change of the four AO trials.
         quartet = hybrids(:, 2:5, ik)
         IF (maxval(abs(matmul(quartet, conjg(transpose(quartet))) - &
                        matmul(candidates(:, 1:4, ik), conjg(transpose(candidates(:, 1:4, ik)))))) &
             > 1.0e-14_dp) nfail = nfail + 1
         candidates(:, :, ik) = matmul(rotation(:, :, ik), candidates(:, :, ik))
      END DO
   END IF
   CALL wannier90_hybrid_candidates(candidates, reshape([1, 2, 3, 4], [4, 1]), rotated_hybrids)
   IF (para_env%is_source()) THEN
      DO ik = 1, 2
         IF (maxval(abs(rotated_hybrids(:, :, ik) - matmul(rotation(:, :, ik), hybrids(:, :, ik)))) &
             > 1.0e-14_dp) nfail = nfail + 1
      END DO
   ELSE
      IF (SIZE(hybrids) /= 0 .OR. SIZE(rotated_hybrids) /= 0) nfail = nfail + 1
   END IF
   DEALLOCATE (hybrids, rotated_hybrids)
   CALL wannier90_hybrid_candidates(candidates, reshape([1, 2, 3, 4, 5, 2, 3, 4], [4, 2]), hybrids)
   IF (para_env%is_source()) THEN
      IF (.NOT. all(shape(hybrids) == [2, 8, 2])) nfail = nfail + 1
      DO ik = 1, 2
         IF (maxval(abs(sum(hybrids(:, 1:4, ik), dim=2)*0.5_dp - candidates(:, 1, ik))) &
             > 1.0e-14_dp) nfail = nfail + 1
         IF (maxval(abs(sum(hybrids(:, 5:8, ik), dim=2)*0.5_dp - candidates(:, 5, ik))) &
             > 1.0e-14_dp) nfail = nfail + 1
      END DO
   END IF
   DEALLOCATE (hybrids)
   CALL wannier90_hybrid_candidates(candidates, empty_groups, hybrids)
   IF (para_env%is_source()) THEN
      IF (maxval(abs(hybrids - candidates)) > 1.0e-14_dp) nfail = nfail + 1
   END IF
   DEALLOCATE (hybrids, candidates)
   IF (para_env%is_source()) THEN
      ALLOCATE (candidates(2, 4, 3))
      candidates = cmplx(0.0_dp, 0.0_dp, kind=dp)
      ! Every QR anchor chooses at least one column that vanishes elsewhere.
      ! Exchanging both of those columns recovers the fixed full-rank pair.
      candidates(1, 1, 1) = 10.0_dp
      candidates(1, 1, 3) = 10.0_dp
      candidates(2, 2, 2) = 10.0_dp
      candidates(2, 2, 3) = 10.0_dp
      candidates(1, 3, :) = 1.0_dp
      candidates(2, 4, :) = 1.0_dp
   ELSE
      ALLOCATE (candidates(0, 0, 0))
   END IF
   DO anchor = 0, 3
      CALL wannier90_select_projections(candidates, para_env, repaired, pivots, min_svalue, anchor)
      IF (min_svalue > 1.0e-14_dp) nfail = nfail + 1
   END DO
   CALL wannier90_select_projections(candidates, para_env, repaired, pivots, min_svalue)
   IF (abs(min_svalue - 1.0_dp) > 1.0e-14_dp) nfail = nfail + 1
   IF (.NOT. all(pivots == [3, 4])) nfail = nfail + 1
   DEALLOCATE (candidates)
   IF (para_env%is_source()) THEN
      ALLOCATE (candidates(2, 3, 2))
      candidates = cmplx(0.0_dp, 0.0_dp, kind=dp)
      candidates(1, 1, :) = 1.0_dp
      candidates(1, 2, :) = 1.0_dp + 1.0e-14_dp
      candidates(2, 3, :) = 1.0_dp
   ELSE
      ALLOCATE (candidates(0, 0, 0))
   END IF
   CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue)
   IF (.NOT. all(pivots == [1, 3])) nfail = nfail + 1
   IF (para_env%is_source()) THEN
      DO ik = 1, 2
         candidates(:, :, ik) = matmul(rotation(:, :, ik), candidates(:, :, ik))
      END DO
   END IF
   CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue)
   IF (.NOT. all(pivots == [1, 3])) nfail = nfail + 1
   DEALLOCATE (candidates)
   DO ik = 1, 2
      DO anchor = 1, 3
         mixed_reference(1, anchor, ik) = cmplx(0.1_dp*anchor, 0.2_dp*ik, kind=dp)
         mixed_reference(2, anchor, ik) = cmplx(-0.3_dp*anchor, 0.05_dp*ik*anchor, kind=dp)
         mixed_reference(3, anchor, ik) = cmplx(0.2_dp*anchor**2, 0.3_dp*ik, kind=dp)
         mixed_reference(4, anchor, ik) = cmplx(0.7_dp, -0.1_dp*ik*anchor, kind=dp)
      END DO
   END DO
   min_svalue = wannier90_projection_rank(mixed_reference)
#if defined (__HAS_IEEE_EXCEPTIONS)
   CALL ieee_get_halting_mode(ieee_all, actual_halt)
   IF (any(actual_halt .NEQV. expected_halt)) nfail = nfail + 1
   ! Also check that a caller with traps disabled is left unchanged.
   CALL ieee_set_halting_mode(ieee_all, .false.)
   IF (abs(wannier90_projection_rank(mixed_reference) - min_svalue) > 1.0e-14_dp) nfail = nfail + 1
   CALL ieee_get_halting_mode(ieee_all, actual_halt)
   IF (any(actual_halt)) nfail = nfail + 1
   CALL ieee_set_halting_mode(ieee_all, expected_halt)
#endif
   CALL wannier90_window_mask([-2.0_dp, -1.0_dp, 1.0_dp, 2.0_dp], 2, &
                              [-1.5_dp, 1.5_dp, 0.0_dp, 0.0_dp], &
                              [.true., .true., .false., .false.], allowed)
   IF (any(allowed .NEQV. [.false., .true., .true., .false.])) nfail = nfail + 1
   CALL wannier90_window_mask([-2.0_dp, -1.0_dp, 1.0_dp, 2.0_dp], 2, &
                              [-3.0_dp, 3.0_dp, 0.0_dp, 0.0_dp], &
                              [.true., .true., .false., .true.], allowed)
   IF (any(allowed .NEQV. [.true., .true., .false., .false.])) nfail = nfail + 1
#if defined (__HAS_IEEE_EXCEPTIONS)
   CALL ieee_get_halting_mode(ieee_all, actual_halt)
   IF (any(actual_halt .NEQV. expected_halt)) nfail = nfail + 1
   CALL ieee_set_halting_mode(ieee_all, saved_halt)
#endif
   CALL comm%sum(nfail)
   cpassert(nfail == 0)
   CALL mp_para_env_release(para_env)
   CALL mp_world_finalize()
END PROGRAM wannier90_projections_unittest