Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief
      10             : !> \author Jan Wilhelm
      11             : !> \date 05.2024
      12             : ! **************************************************************************************************
      13             : MODULE gw_kp_to_real_space_and_back
      14             :    USE cp_cfm_types,                    ONLY: cp_cfm_type
      15             :    USE cp_fm_types,                     ONLY: cp_fm_set_all,&
      16             :                                               cp_fm_type
      17             :    USE kinds,                           ONLY: dp
      18             :    USE kpoint_types,                    ONLY: kpoint_type
      19             :    USE mathconstants,                   ONLY: gaussi,&
      20             :                                               twopi,&
      21             :                                               z_one,&
      22             :                                               z_zero
      23             : #include "./base/base_uses.f90"
      24             : 
      25             :    IMPLICIT NONE
      26             : 
      27             :    PRIVATE
      28             : 
      29             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_kp_to_real_space_and_back'
      30             : 
      31             :    PUBLIC :: fm_trafo_rs_to_ikp, trafo_rs_to_ikp, trafo_ikp_to_rs, fm_add_ikp_to_rs, &
      32             :              add_ikp_to_all_rs
      33             : 
      34             : CONTAINS
      35             : 
      36             : ! **************************************************************************************************
      37             : !> \brief ...
      38             : !> \param cfm_ikp ...
      39             : !> \param fm_rs ...
      40             : !> \param kpoints ...
      41             : !> \param ikp ...
      42             : ! **************************************************************************************************
      43        4364 :    SUBROUTINE fm_trafo_rs_to_ikp(cfm_ikp, fm_rs, kpoints, ikp)
      44             :       TYPE(cp_cfm_type)                                  :: cfm_ikp
      45             :       TYPE(cp_fm_type), DIMENSION(:)                     :: fm_rs
      46             :       TYPE(kpoint_type), POINTER                         :: kpoints
      47             :       INTEGER                                            :: ikp
      48             : 
      49             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_trafo_rs_to_ikp'
      50             : 
      51             :       INTEGER                                            :: handle, img, nimages, nimages_fm_rs
      52             : 
      53        4364 :       CALL timeset(routineN, handle)
      54             : 
      55        4364 :       nimages = SIZE(kpoints%index_to_cell, 1)
      56        4364 :       nimages_fm_rs = SIZE(fm_rs)
      57             : 
      58        4364 :       CPASSERT(nimages == nimages_fm_rs)
      59             : 
      60      246382 :       cfm_ikp%local_data(:, :) = z_zero
      61       43640 :       DO img = 1, nimages
      62             : 
      63             :          CALL add_rs_to_ikp(fm_rs(img)%local_data, cfm_ikp%local_data, kpoints%index_to_cell, &
      64       43640 :                             kpoints%xkp(1:3, ikp), img)
      65             : 
      66             :       END DO
      67             : 
      68        4364 :       CALL timestop(handle)
      69             : 
      70        4364 :    END SUBROUTINE fm_trafo_rs_to_ikp
      71             : 
      72             : ! **************************************************************************************************
      73             : !> \brief ...
      74             : !> \param array_rs ...
      75             : !> \param array_kp ...
      76             : !> \param index_to_cell ...
      77             : !> \param xkp ...
      78             : ! **************************************************************************************************
      79       21408 :    SUBROUTINE trafo_rs_to_ikp(array_rs, array_kp, index_to_cell, xkp)
      80             :       REAL(KIND=dp), DIMENSION(:, :, :)                  :: array_rs
      81             :       COMPLEX(KIND=dp), DIMENSION(:, :)                  :: array_kp
      82             :       INTEGER, DIMENSION(:, :)                           :: index_to_cell
      83             :       REAL(KIND=dp)                                      :: xkp(3)
      84             : 
      85             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'trafo_rs_to_ikp'
      86             : 
      87             :       INTEGER                                            :: handle, i_cell, nimages
      88             : 
      89       21408 :       CALL timeset(routineN, handle)
      90             : 
      91       21408 :       nimages = SIZE(index_to_cell, 1)
      92             : 
      93       21408 :       CPASSERT(nimages == SIZE(array_rs, 3))
      94             : 
      95      727776 :       array_kp(:, :) = 0.0_dp
      96      214080 :       DO i_cell = 1, nimages
      97             : 
      98      214080 :          CALL add_rs_to_ikp(array_rs(:, :, i_cell), array_kp, index_to_cell, xkp, i_cell)
      99             : 
     100             :       END DO
     101             : 
     102       21408 :       CALL timestop(handle)
     103             : 
     104       21408 :    END SUBROUTINE trafo_rs_to_ikp
     105             : 
     106             : ! **************************************************************************************************
     107             : !> \brief ...
     108             : !> \param array_rs ...
     109             : !> \param array_kp ...
     110             : !> \param index_to_cell ...
     111             : !> \param xkp ...
     112             : !> \param i_cell ...
     113             : ! **************************************************************************************************
     114      231948 :    SUBROUTINE add_rs_to_ikp(array_rs, array_kp, index_to_cell, xkp, i_cell)
     115             :       REAL(KIND=dp), DIMENSION(:, :)                     :: array_rs
     116             :       COMPLEX(KIND=dp), DIMENSION(:, :)                  :: array_kp
     117             :       INTEGER, DIMENSION(:, :)                           :: index_to_cell
     118             :       REAL(KIND=dp)                                      :: xkp(3)
     119             :       INTEGER                                            :: i_cell
     120             : 
     121             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'add_rs_to_ikp'
     122             : 
     123             :       COMPLEX(KIND=dp)                                   :: expikR
     124             :       INTEGER                                            :: handle
     125             :       REAL(KIND=dp)                                      :: arg
     126             : 
     127      231948 :       CALL timeset(routineN, handle)
     128             : 
     129             :       arg = REAL(index_to_cell(i_cell, 1), dp)*xkp(1) + &
     130             :             REAL(index_to_cell(i_cell, 2), dp)*xkp(2) + &
     131      231948 :             REAL(index_to_cell(i_cell, 3), dp)*xkp(3)
     132             : 
     133      231948 :       expikR = z_one*COS(twopi*arg) + gaussi*SIN(twopi*arg)
     134             : 
     135     8767422 :       array_kp(:, :) = array_kp(:, :) + expikR*array_rs(:, :)
     136             : 
     137      231948 :       CALL timestop(handle)
     138             : 
     139      231948 :    END SUBROUTINE add_rs_to_ikp
     140             : 
     141             : ! **************************************************************************************************
     142             : !> \brief ...
     143             : !> \param array_kp ...
     144             : !> \param array_rs ...
     145             : !> \param cell ...
     146             : !> \param kpoints ...
     147             : ! **************************************************************************************************
     148           0 :    SUBROUTINE trafo_ikp_to_rs(array_kp, array_rs, cell, kpoints)
     149             :       COMPLEX(KIND=dp), DIMENSION(:, :, :)               :: array_kp
     150             :       REAL(KIND=dp), DIMENSION(:, :)                     :: array_rs
     151             :       INTEGER                                            :: cell(3)
     152             :       TYPE(kpoint_type), POINTER                         :: kpoints
     153             : 
     154             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'trafo_ikp_to_rs'
     155             : 
     156             :       INTEGER                                            :: handle, ikp
     157             : 
     158           0 :       CALL timeset(routineN, handle)
     159             : 
     160           0 :       CPASSERT(kpoints%nkp == SIZE(array_kp, 3))
     161             : 
     162           0 :       array_rs(:, :) = 0.0_dp
     163             : 
     164           0 :       DO ikp = 1, kpoints%nkp
     165             : 
     166           0 :          CALL add_ikp_to_rs(array_kp(:, :, ikp), array_rs, cell, kpoints, ikp)
     167             : 
     168             :       END DO
     169             : 
     170           0 :       CALL timestop(handle)
     171             : 
     172           0 :    END SUBROUTINE trafo_ikp_to_rs
     173             : 
     174             : ! **************************************************************************************************
     175             : !> \brief ...
     176             : !> \param cfm_ikp ...
     177             : !> \param fm_rs ...
     178             : !> \param kpoints ...
     179             : !> \param ikp ...
     180             : ! **************************************************************************************************
     181        3424 :    SUBROUTINE fm_add_ikp_to_rs(cfm_ikp, fm_rs, kpoints, ikp)
     182             :       TYPE(cp_cfm_type)                                  :: cfm_ikp
     183             :       TYPE(cp_fm_type), DIMENSION(:)                     :: fm_rs
     184             :       TYPE(kpoint_type), POINTER                         :: kpoints
     185             :       INTEGER                                            :: ikp
     186             : 
     187             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'fm_add_ikp_to_rs'
     188             : 
     189             :       INTEGER                                            :: handle, img, nimages, nimages_fm_rs
     190        3424 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: index_to_cell
     191             : 
     192        3424 :       CALL timeset(routineN, handle)
     193             : 
     194        3424 :       nimages = SIZE(kpoints%index_to_cell, 1)
     195        3424 :       nimages_fm_rs = SIZE(fm_rs)
     196             : 
     197        3424 :       CPASSERT(nimages == nimages_fm_rs)
     198             : 
     199       10272 :       ALLOCATE (index_to_cell(nimages, 3))
     200      106144 :       index_to_cell(1:nimages, 1:3) = kpoints%index_to_cell(1:nimages, 1:3)
     201             : 
     202       34240 :       DO img = 1, nimages
     203             : 
     204       30816 :          IF (ikp == 1) CALL cp_fm_set_all(fm_rs(img), 0.0_dp)
     205             : 
     206             :          CALL add_ikp_to_rs(cfm_ikp%local_data(:, :), fm_rs(img)%local_data, &
     207      219136 :                             index_to_cell(img, 1:3), kpoints, ikp)
     208             : 
     209             :       END DO
     210             : 
     211        3424 :       CALL timestop(handle)
     212             : 
     213        6848 :    END SUBROUTINE fm_add_ikp_to_rs
     214             : 
     215             : ! **************************************************************************************************
     216             : !> \brief ...
     217             : !> \param array_kp ...
     218             : !> \param array_rs ...
     219             : !> \param kpoints ...
     220             : !> \param ikp ...
     221             : !> \param index_to_cell_ext ...
     222             : ! **************************************************************************************************
     223       17104 :    SUBROUTINE add_ikp_to_all_rs(array_kp, array_rs, kpoints, ikp, index_to_cell_ext)
     224             :       COMPLEX(KIND=dp), DIMENSION(:, :)                  :: array_kp
     225             :       REAL(KIND=dp), DIMENSION(:, :, :)                  :: array_rs
     226             :       TYPE(kpoint_type), POINTER                         :: kpoints
     227             :       INTEGER                                            :: ikp
     228             :       INTEGER, DIMENSION(:, :), OPTIONAL, POINTER        :: index_to_cell_ext
     229             : 
     230             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'add_ikp_to_all_rs'
     231             : 
     232             :       INTEGER                                            :: cell(3), handle, img, nimages
     233       17104 :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell
     234             : 
     235       17104 :       CALL timeset(routineN, handle)
     236             : 
     237       17104 :       IF (PRESENT(index_to_cell_ext)) THEN
     238       16640 :          index_to_cell => index_to_cell_ext
     239             :       ELSE
     240         464 :          index_to_cell => kpoints%index_to_cell
     241             :       END IF
     242             : 
     243       17104 :       nimages = SIZE(index_to_cell, 1)
     244       17104 :       CPASSERT(SIZE(array_rs, 3) == nimages)
     245      171040 :       DO img = 1, nimages
     246             : 
     247      615744 :          cell(1:3) = index_to_cell(img, 1:3)
     248             : 
     249      171040 :          CALL add_ikp_to_rs(array_kp, array_rs(:, :, img), cell, kpoints, ikp)
     250             : 
     251             :       END DO
     252             : 
     253       17104 :       CALL timestop(handle)
     254             : 
     255       17104 :    END SUBROUTINE add_ikp_to_all_rs
     256             : 
     257             : ! **************************************************************************************************
     258             : !> \brief ...
     259             : !> \param array_kp ...
     260             : !> \param array_rs ...
     261             : !> \param cell ...
     262             : !> \param kpoints ...
     263             : !> \param ikp ...
     264             : ! **************************************************************************************************
     265      184752 :    SUBROUTINE add_ikp_to_rs(array_kp, array_rs, cell, kpoints, ikp)
     266             :       COMPLEX(KIND=dp), DIMENSION(:, :)                  :: array_kp
     267             :       REAL(KIND=dp), DIMENSION(:, :)                     :: array_rs
     268             :       INTEGER                                            :: cell(3)
     269             :       TYPE(kpoint_type), POINTER                         :: kpoints
     270             :       INTEGER                                            :: ikp
     271             : 
     272             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'add_ikp_to_rs'
     273             : 
     274             :       INTEGER                                            :: handle
     275             :       REAL(KIND=dp)                                      :: arg, im, re
     276             : 
     277      184752 :       CALL timeset(routineN, handle)
     278             : 
     279             :       arg = REAL(cell(1), dp)*kpoints%xkp(1, ikp) + &
     280             :             REAL(cell(2), dp)*kpoints%xkp(2, ikp) + &
     281      184752 :             REAL(cell(3), dp)*kpoints%xkp(3, ikp)
     282             : 
     283      184752 :       re = COS(twopi*arg)*kpoints%wkp(ikp)
     284      184752 :       im = SIN(twopi*arg)*kpoints%wkp(ikp)
     285             : 
     286     7205040 :       array_rs(:, :) = array_rs(:, :) + re*REAL(array_kp(:, :)) + im*AIMAG(array_kp(:, :))
     287             : 
     288      184752 :       CALL timestop(handle)
     289             : 
     290      184752 :    END SUBROUTINE add_ikp_to_rs
     291             : 
     292             : END MODULE gw_kp_to_real_space_and_back