Line data    Source code

       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Inversion action on CP2K Gaussian Bloch coefficients, including lattice phases.
      10              : ! **************************************************************************************************
      11              : MODULE topology_inversion
      12              :    USE atomic_kind_types,               ONLY: get_atomic_kind
      13              :    USE basis_set_types,                 ONLY: gto_basis_set_type
      14              :    USE cell_types,                      ONLY: cell_type,&
      15              :                                               pbc
      16              :    USE ieee_arithmetic,                 ONLY: ieee_is_finite
      17              :    USE kinds,                           ONLY: dp
      18              :    USE mathconstants,                   ONLY: twopi
      19              :    USE particle_types,                  ONLY: particle_type
      20              :    USE qs_environment_types,            ONLY: get_qs_env,&
      21              :                                               qs_environment_type
      22              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      23              :                                               qs_kind_type
      24              : 
      25              :    IMPLICIT NONE
      26              :    PRIVATE
      27              :    PUBLIC :: gaussian_inversion_action
      28              : CONTAINS
      29              : 
      30              : ! **************************************************************************************************
      31              : !> \brief Construct sparse AO inversion at a TRIM for the operation r -> 2*center-r.
      32              : !> \param qs_env Electronic environment
      33              : !> \param k Fractional reciprocal coordinate at a TRIM
      34              : !> \param center Fractional inversion center in the user's cell setting
      35              : !> \param mapping Target AO row for every source row
      36              : !> \param phase AO parity times exp(+2*pi*i*k.L), with 2c-r_a=r_b+L
      37              : !> \param tolerance Geometry tolerance in Bohr
      38              : !> \param status Zero on success, negative if geometry/basis/periodicity is incompatible
      39              : ! **************************************************************************************************
      40           32 :    SUBROUTINE gaussian_inversion_action(qs_env, k, center, mapping, phase, tolerance, status)
      41              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      42              :       REAL(KIND=dp), INTENT(IN)                          :: k(3), center(3)
      43              :       INTEGER, INTENT(OUT)                               :: mapping(:)
      44              :       COMPLEX(KIND=dp), INTENT(OUT)                      :: phase(:)
      45              :       REAL(KIND=dp), INTENT(IN)                          :: tolerance
      46              :       INTEGER, INTENT(OUT)                               :: status
      47              : 
      48              :       INTEGER                                            :: a, b, found, i, iset, l, n, shell, total
      49           32 :       INTEGER, ALLOCATABLE                               :: kind(:), offset(:), TARGET(:)
      50              :       REAL(KIND=dp)                                      :: angle, delta(3), image(3)
      51           32 :       REAL(KIND=dp), ALLOCATABLE                         :: position(:, :)
      52              :       TYPE(cell_type), POINTER                           :: cell
      53              :       TYPE(gto_basis_set_type), POINTER                  :: basis
      54           32 :       TYPE(particle_type), POINTER                       :: particles(:)
      55           32 :       TYPE(qs_kind_type), POINTER                        :: kinds_set(:)
      56              : 
      57           32 :       status = -1
      58          552 :       mapping = 0
      59          552 :       phase = CMPLX(0, 0, dp)
      60           32 :       IF (tolerance <= 0.0_dp .OR. SIZE(phase) /= SIZE(mapping)) RETURN
      61          128 :       IF (.NOT. ieee_is_finite(tolerance) .OR. .NOT. ALL(ieee_is_finite(center))) RETURN
      62          128 :       IF (.NOT. ALL(ieee_is_finite(k))) RETURN
      63          128 :       IF (MAXVAL(ABS(2*k - NINT(2*k))) > 1.e-8_dp) RETURN
      64           32 :       CALL get_qs_env(qs_env, cell=cell, particle_set=particles, qs_kind_set=kinds_set)
      65          128 :       IF (ANY(cell%perd /= 1)) RETURN
      66           32 :       n = SIZE(particles)
      67          224 :       ALLOCATE (KIND(n), offset(n), TARGET(n), position(3, n))
      68           72 :       total = 0
      69           72 :       DO a = 1, n
      70           40 :          CALL get_atomic_kind(particles(a)%atomic_kind, kind_number=KIND(a))
      71           40 :          CALL get_qs_kind(kinds_set(KIND(a)), basis_set=basis)
      72           40 :          IF (.NOT. ASSOCIATED(basis)) RETURN
      73           40 :          offset(a) = total
      74           40 :          total = total + basis%nsgf
      75          672 :          position(:, a) = MATMUL(cell%h_inv, pbc(particles(a)%r, cell))
      76              :       END DO
      77           32 :       IF (total /= SIZE(mapping)) RETURN
      78           72 :       DO a = 1, n
      79              :          found = 0
      80           96 :          DO b = 1, n
      81           56 :             IF (KIND(a) /= KIND(b)) CYCLE
      82          224 :             delta = 2*center - position(:, a) - position(:, b)
      83         1064 :             IF (SQRT(SUM(MATMUL(cell%hmat, delta - NINT(delta))**2)) > tolerance) CYCLE
      84           40 :             found = found + 1
      85           80 :             TARGET(a) = b
      86              :          END DO
      87           72 :          IF (found /= 1) RETURN
      88              :       END DO
      89           72 :       DO a = 1, n
      90           40 :          b = TARGET(a)
      91           40 :          IF (TARGET(b) /= a) RETURN
      92          160 :          image = REAL(NINT(2*center - position(:, a) - position(:, b)), dp)
      93          160 :          angle = twopi*DOT_PRODUCT(k, image)
      94           40 :          CALL get_qs_kind(kinds_set(KIND(a)), basis_set=basis)
      95          136 :          DO iset = 1, basis%nset
      96          304 :             DO shell = 1, basis%nshell(iset)
      97          200 :                l = basis%l(shell, iset)
      98          784 :                DO i = basis%first_sgf(shell, iset), basis%last_sgf(shell, iset)
      99          520 :                   mapping(offset(a) + i) = offset(b) + i
     100          720 :                   phase(offset(a) + i) = (-1)**l*CMPLX(COS(angle), SIN(angle), dp)
     101              :                END DO
     102              :             END DO
     103              :          END DO
     104              :       END DO
     105          552 :       IF (ANY(mapping == 0)) RETURN
     106           32 :       status = 0
     107           64 :    END SUBROUTINE gaussian_inversion_action
     108              : END MODULE topology_inversion