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             : MODULE mt_util
      10             :    USE bibliography,                    ONLY: Martyna1999,&
      11             :                                               cite_reference
      12             :    USE kinds,                           ONLY: dp
      13             :    USE mathconstants,                   ONLY: fourpi,&
      14             :                                               oorootpi,&
      15             :                                               pi
      16             :    USE pw_grid_types,                   ONLY: pw_grid_type
      17             :    USE pw_methods,                      ONLY: pw_axpy,&
      18             :                                               pw_transfer,&
      19             :                                               pw_zero
      20             :    USE pw_pool_types,                   ONLY: pw_pool_create,&
      21             :                                               pw_pool_release,&
      22             :                                               pw_pool_type
      23             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      24             :                                               pw_r3d_rs_type
      25             : #include "../base/base_uses.f90"
      26             : 
      27             :    IMPLICIT NONE
      28             : 
      29             :    PRIVATE
      30             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      31             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mt_util'
      32             : 
      33             :    INTEGER, PARAMETER, PUBLIC               :: MT2D = 1101, &
      34             :                                                MT1D = 1102, &
      35             :                                                MT0D = 1103
      36             : 
      37             :    PUBLIC :: MTin_create_screen_fn
      38             : CONTAINS
      39             : 
      40             : ! **************************************************************************************************
      41             : !> \brief Initialize the Martyna && Tuckerman Poisson Solver
      42             : !> \param screen_function ...
      43             : !> \param pw_pool ...
      44             : !> \param method ...
      45             : !> \param alpha ...
      46             : !> \param special_dimension ...
      47             : !> \param slab_size ...
      48             : !> \param super_ref_pw_grid ...
      49             : !> \author Teodoro Laino (16.06.2004)
      50             : ! **************************************************************************************************
      51        1106 :    SUBROUTINE MTin_create_screen_fn(screen_function, pw_pool, method, alpha, &
      52             :                                     special_dimension, slab_size, super_ref_pw_grid)
      53             :       TYPE(pw_c1d_gs_type), POINTER                      :: screen_function
      54             :       TYPE(pw_pool_type), POINTER                        :: pw_pool
      55             :       INTEGER, INTENT(IN)                                :: method
      56             :       REAL(KIND=dp), INTENT(in)                          :: alpha
      57             :       INTEGER, INTENT(IN)                                :: special_dimension
      58             :       REAL(KIND=dp), INTENT(in)                          :: slab_size
      59             :       TYPE(pw_grid_type), POINTER                        :: super_ref_pw_grid
      60             : 
      61             :       CHARACTER(len=*), PARAMETER :: routineN = 'MTin_create_screen_fn'
      62             : 
      63             :       INTEGER                                            :: handle, ig, iz
      64             :       REAL(KIND=dp)                                      :: alpha2, g2, g3d, gxy, gz, zlength
      65             :       TYPE(pw_c1d_gs_type), POINTER                      :: Vlocg
      66             :       TYPE(pw_pool_type), POINTER                        :: pw_pool_aux
      67             :       TYPE(pw_r3d_rs_type), POINTER                      :: Vloc
      68             : 
      69        1106 :       CALL timeset(routineN, handle)
      70        1106 :       NULLIFY (Vloc, Vlocg, pw_pool_aux)
      71             :       !
      72             :       ! For Martyna-Tuckerman we set up an auxiliary pw_pool at an higher cutoff
      73             :       !
      74        1106 :       CALL cite_reference(Martyna1999)
      75        1106 :       IF (ASSOCIATED(super_ref_pw_grid)) THEN
      76        1100 :          CALL pw_pool_create(pw_pool_aux, pw_grid=super_ref_pw_grid)
      77             :       END IF
      78             :       NULLIFY (screen_function)
      79        1106 :       ALLOCATE (screen_function)
      80        1106 :       CALL pw_pool%create_pw(screen_function)
      81        1106 :       CALL pw_zero(screen_function)
      82        2210 :       SELECT CASE (method)
      83             :       CASE (MT0D)
      84        1104 :          NULLIFY (Vloc, Vlocg)
      85        1104 :          ALLOCATE (Vloc, Vlocg)
      86        1104 :          IF (ASSOCIATED(pw_pool_aux)) THEN
      87        1098 :             CALL pw_pool_aux%create_pw(Vloc)
      88        1098 :             CALL pw_pool_aux%create_pw(Vlocg)
      89             :          ELSE
      90           6 :             CALL pw_pool%create_pw(Vloc)
      91           6 :             CALL pw_pool%create_pw(Vlocg)
      92             :          END IF
      93        1104 :          CALL mt0din(Vloc, alpha)
      94        1104 :          CALL pw_transfer(Vloc, Vlocg)
      95        1104 :          CALL pw_axpy(Vlocg, screen_function)
      96        1104 :          IF (ASSOCIATED(pw_pool_aux)) THEN
      97        1098 :             CALL pw_pool_aux%give_back_pw(Vloc)
      98        1098 :             CALL pw_pool_aux%give_back_pw(Vlocg)
      99             :          ELSE
     100           6 :             CALL pw_pool%give_back_pw(Vloc)
     101           6 :             CALL pw_pool%give_back_pw(Vlocg)
     102             :          END IF
     103        1104 :          DEALLOCATE (Vloc, Vlocg)
     104             :          !
     105             :          ! Get rid of the analytical FT of the erf(a*r)/r
     106             :          !
     107        1104 :          alpha2 = alpha*alpha
     108    73270559 :          DO ig = screen_function%pw_grid%first_gne0, screen_function%pw_grid%ngpts_cut_local
     109    73269455 :             g2 = screen_function%pw_grid%gsq(ig)
     110    73269455 :             g3d = fourpi/g2
     111    73270559 :             screen_function%array(ig) = screen_function%array(ig) - g3d*EXP(-g2/(4.0E0_dp*alpha2))
     112             :          END DO
     113        1104 :          IF (screen_function%pw_grid%have_g0) &
     114         558 :             screen_function%array(1) = screen_function%array(1) + fourpi/(4.0E0_dp*alpha2)
     115             :       CASE (MT2D)
     116           2 :          iz = special_dimension ! iz is the direction with NO PBC
     117           2 :          zlength = slab_size ! zlength is the thickness of the cell
     118      194401 :          DO ig = screen_function%pw_grid%first_gne0, screen_function%pw_grid%ngpts_cut_local
     119      194399 :             gz = screen_function%pw_grid%g(iz, ig)
     120      194399 :             g2 = screen_function%pw_grid%gsq(ig)
     121      194399 :             gxy = SQRT(ABS(g2 - gz*gz))
     122      194399 :             g3d = fourpi/g2
     123      194401 :             screen_function%array(ig) = -g3d*COS(gz*zlength/2.0_dp)*EXP(-gxy*zlength/2.0_dp)
     124             :          END DO
     125           2 :          IF (screen_function%pw_grid%have_g0) screen_function%array(1) = pi*zlength*zlength/2.0_dp
     126             :       CASE (MT1D)
     127           0 :          iz = special_dimension ! iz is the direction with PBC
     128           0 :          CALL mt1din(screen_function)
     129        1106 :          CPABORT("MT1D unimplemented")
     130             :       END SELECT
     131        1106 :       CALL pw_pool_release(pw_pool_aux)
     132        1106 :       CALL timestop(handle)
     133        1106 :    END SUBROUTINE MTin_create_screen_fn
     134             : 
     135             : ! **************************************************************************************************
     136             : !> \brief Calculates the Tuckerman Green's function in reciprocal space
     137             : !>      according the scheme published on:
     138             : !>      Martyna and Tuckerman, J. Chem. Phys. Vol. 110, No. 6, 2810-2821
     139             : !> \param Vloc ...
     140             : !> \param alpha ...
     141             : !> \author Teodoro Laino (09.03.2005)
     142             : ! **************************************************************************************************
     143        1104 :    SUBROUTINE mt0din(Vloc, alpha)
     144             :       TYPE(pw_r3d_rs_type), POINTER                      :: Vloc
     145             :       REAL(KIND=dp), INTENT(in)                          :: alpha
     146             : 
     147             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'mt0din'
     148             : 
     149             :       INTEGER                                            :: handle, i, ii, j, jj, k, kk
     150        1104 :       INTEGER, DIMENSION(:), POINTER                     :: glb
     151        1104 :       INTEGER, DIMENSION(:, :), POINTER                  :: bo
     152             :       REAL(KIND=dp)                                      :: dx, dy, dz, fact, omega, r, r2, x, y, &
     153             :                                                             y2, z, z2
     154             :       REAL(KIND=dp), DIMENSION(3)                        :: box, box2
     155             :       TYPE(pw_grid_type), POINTER                        :: grid
     156             : 
     157        1104 :       CALL timeset(routineN, handle)
     158             : 
     159        1104 :       grid => Vloc%pw_grid
     160        1104 :       bo => grid%bounds_local
     161        1104 :       glb => grid%bounds(1, :)
     162   197052523 :       Vloc%array = 0.0_dp
     163        4416 :       box = REAL(grid%npts, kind=dp)*grid%dr
     164        4416 :       box2 = box/2.0_dp
     165        4416 :       omega = PRODUCT(box)
     166        1104 :       fact = omega
     167        1104 :       dx = grid%dr(1)
     168        1104 :       dy = grid%dr(2)
     169        1104 :       dz = grid%dr(3)
     170        1104 :       kk = bo(1, 3)
     171       72192 :       DO k = bo(1, 3), bo(2, 3)
     172       71088 :          z = REAL(k - glb(3), dp)*dz; IF (z .GT. box2(3)) z = box(3) - z
     173       71088 :          z2 = z*z
     174       71088 :          jj = bo(1, 2)
     175     5021000 :          DO j = bo(1, 2), bo(2, 2)
     176     4949912 :             y = REAL(j - glb(2), dp)*dy; IF (y .GT. box2(2)) y = box(2) - y
     177     4949912 :             y2 = y*y
     178     4949912 :             ii = bo(1, 1)
     179   196980331 :             DO i = bo(1, 1), bo(2, 1)
     180   192030419 :                x = REAL(i - glb(1), dp)*dx; IF (x .GT. box2(1)) x = box(1) - x
     181   192030419 :                r2 = x*x + y2 + z2
     182   192030419 :                r = SQRT(r2)
     183   192030419 :                IF (r .GT. 1.0E-10_dp) THEN
     184   192029861 :                   Vloc%array(ii, jj, kk) = erf(alpha*r)/r*fact
     185             :                ELSE
     186         558 :                   Vloc%array(ii, jj, kk) = 2.0_dp*alpha*oorootpi*fact
     187             :                END IF
     188   196980331 :                ii = ii + 1
     189             :             END DO
     190     5021000 :             jj = jj + 1
     191             :          END DO
     192       72192 :          kk = kk + 1
     193             :       END DO
     194        1104 :       CALL timestop(handle)
     195        1104 :    END SUBROUTINE Mt0din
     196             : 
     197             : ! **************************************************************************************************
     198             : !> \brief Calculates the Tuckerman Green's function in reciprocal space
     199             : !>      according the scheme published on:
     200             : !>      Martyna and Tuckerman, J. Chem. Phys. Vol. 121, No. 23, 11949
     201             : !> \param screen_function ...
     202             : !> \author Teodoro Laino (11.2005)
     203             : ! **************************************************************************************************
     204           0 :    SUBROUTINE mt1din(screen_function)
     205             :       TYPE(pw_c1d_gs_type), POINTER                      :: screen_function
     206             : 
     207             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'mt1din'
     208             : 
     209             :       INTEGER                                            :: handle
     210             :       REAL(KIND=dp)                                      :: dx, dy, dz, omega
     211             :       REAL(KIND=dp), DIMENSION(3)                        :: box, box2
     212             :       TYPE(pw_grid_type), POINTER                        :: grid
     213             : 
     214           0 :       CALL timeset(routineN, handle)
     215           0 :       grid => screen_function%pw_grid
     216           0 :       box = REAL(grid%npts, kind=dp)*grid%dr
     217           0 :       box2 = box/2.0_dp
     218           0 :       omega = PRODUCT(box)
     219           0 :       dx = grid%dr(1)
     220           0 :       dy = grid%dr(2)
     221           0 :       dz = grid%dr(3)
     222             : 
     223           0 :       CALL timestop(handle)
     224           0 :    END SUBROUTINE mt1din
     225             : 
     226             : END MODULE mt_util