mt_util.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                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
MODULE mt_util
   USE bibliography,                    ONLY: martyna1999,&
                                              cite_reference
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: fourpi,&
                                              oorootpi,&
                                              pi
   USE pw_grid_types,                   ONLY: pw_grid_type
   USE pw_methods,                      ONLY: pw_axpy,&
                                              pw_transfer,&
                                              pw_zero
   USE pw_pool_types,                   ONLY: pw_pool_create,&
                                              pw_pool_release,&
                                              pw_pool_type
   USE pw_types,                        ONLY: pw_c1d_gs_type,&
                                              pw_r3d_rs_type
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mt_util'
 
   INTEGER, PARAMETER, PUBLIC               :: mt2d = 1101, &
                                               mt1d = 1102, &
                                               mt0d = 1103
 
   PUBLIC :: mtin_create_screen_fn, mt2d_kernel_strain
CONTAINS
 
! **************************************************************************************************
!> \brief MT2D Coulomb kernel and its affine strain derivative in the periodic plane.
!> \param g Cartesian reciprocal vector, transformed by the inverse transpose of the strain.
!> \param special_dimension Cartesian direction normal to the slab.
!> \param slab_size Cell thickness in the nonperiodic direction, fixed for in-plane strain.
!> \param kernel Reciprocal Coulomb kernel, including the existing MT2D G=0 convention.
!> \param strain Kernel derivative with respect to in-plane affine strain.
! **************************************************************************************************
   SUBROUTINE mt2d_kernel_strain(g, special_dimension, slab_size, kernel, strain)
      REAL(kind=dp), INTENT(IN)                          :: g(3)
      INTEGER, INTENT(IN)                                :: special_dimension
      REAL(kind=dp), INTENT(IN)                          :: slab_size
      REAL(kind=dp), INTENT(OUT)                         :: kernel, strain(3, 3)
 
      INTEGER                                            :: alpha, beta
      REAL(kind=dp)                                      :: f, g2, k, parallel(3), q, radius
 
      cpassert(special_dimension >= 1 .AND. special_dimension <= 3)
      cpassert(slab_size > 0.0_dp)
      strain = 0.0_dp
      g2 = sum(g**2)
      IF (g2 == 0.0_dp) THEN
         ! This kernel is independent of in-plane strain. Its volume term is
         ! already included in the Hartree energy contribution to the virial.
         kernel = pi*slab_size**2/2.0_dp
         RETURN
      END IF
      parallel = g
      parallel(special_dimension) = 0.0_dp
      k = norm2(parallel)
      radius = slab_size/2.0_dp
      q = cos(g(special_dimension)*radius)*exp(-k*radius)
      f = 1.0_dp - q
      kernel = fourpi/g2*f
      ! For k=0 the reciprocal vector is unchanged by every in-plane strain.
      ! This also avoids dividing by k at the even-normal modes where K=0.
      IF (k == 0.0_dp) RETURN
      DO alpha = 1, 3
         IF (alpha == special_dimension) cycle
         DO beta = alpha, 3
            IF (beta == special_dimension) cycle
            strain(alpha, beta) = fourpi/g2*g(alpha)*g(beta)*(2.0_dp*f/g2 - q*radius/k)
            strain(beta, alpha) = strain(alpha, beta)
         END DO
      END DO
   END SUBROUTINE mt2d_kernel_strain
 
! **************************************************************************************************
!> \brief Initialize the Martyna && Tuckerman Poisson Solver
!> \param screen_function ...
!> \param pw_pool ...
!> \param method ...
!> \param alpha ...
!> \param special_dimension ...
!> \param slab_size ...
!> \param super_ref_pw_grid ...
!> \author Teodoro Laino (16.06.2004)
! **************************************************************************************************
   SUBROUTINE mtin_create_screen_fn(screen_function, pw_pool, method, alpha, &
                                    special_dimension, slab_size, super_ref_pw_grid)
      TYPE(pw_c1d_gs_type), POINTER                      :: screen_function
      TYPE(pw_pool_type), POINTER                        :: pw_pool
      INTEGER, INTENT(IN)                                :: method
      REAL(kind=dp), INTENT(in)                          :: alpha
      INTEGER, INTENT(IN)                                :: special_dimension
      REAL(kind=dp), INTENT(in)                          :: slab_size
      TYPE(pw_grid_type), POINTER                        :: super_ref_pw_grid
 
      CHARACTER(len=*), PARAMETER :: routinen = 'MTin_create_screen_fn'
 
      INTEGER                                            :: handle, ig, iz
      REAL(kind=dp)                                      :: alpha2, g2, g3d, gxy, gz, zlength
      TYPE(pw_c1d_gs_type), POINTER                      :: vlocg
      TYPE(pw_pool_type), POINTER                        :: pw_pool_aux
      TYPE(pw_r3d_rs_type), POINTER                      :: vloc
 
      CALL timeset(routinen, handle)
      NULLIFY (vloc, vlocg, pw_pool_aux)
      !
      ! For Martyna-Tuckerman we set up an auxiliary pw_pool at an higher cutoff
      !
      CALL cite_reference(martyna1999)
      IF (ASSOCIATED(super_ref_pw_grid)) THEN
         CALL pw_pool_create(pw_pool_aux, pw_grid=super_ref_pw_grid)
      END IF
      NULLIFY (screen_function)
      ALLOCATE (screen_function)
      CALL pw_pool%create_pw(screen_function)
      CALL pw_zero(screen_function)
      SELECT CASE (method)
      CASE (mt0d)
         NULLIFY (vloc, vlocg)
         ALLOCATE (vloc, vlocg)
         IF (ASSOCIATED(pw_pool_aux)) THEN
            CALL pw_pool_aux%create_pw(vloc)
            CALL pw_pool_aux%create_pw(vlocg)
         ELSE
            CALL pw_pool%create_pw(vloc)
            CALL pw_pool%create_pw(vlocg)
         END IF
         CALL mt0din(vloc, alpha)
         CALL pw_transfer(vloc, vlocg)
         CALL pw_axpy(vlocg, screen_function)
         IF (ASSOCIATED(pw_pool_aux)) THEN
            CALL pw_pool_aux%give_back_pw(vloc)
            CALL pw_pool_aux%give_back_pw(vlocg)
         ELSE
            CALL pw_pool%give_back_pw(vloc)
            CALL pw_pool%give_back_pw(vlocg)
         END IF
         DEALLOCATE (vloc, vlocg)
         !
         ! Get rid of the analytical FT of the erf(a*r)/r
         !
         alpha2 = alpha*alpha
         DO ig = screen_function%pw_grid%first_gne0, screen_function%pw_grid%ngpts_cut_local
            g2 = screen_function%pw_grid%gsq(ig)
            g3d = fourpi/g2
            screen_function%array(ig) = screen_function%array(ig) - g3d*exp(-g2/(4.0e0_dp*alpha2))
         END DO
         IF (screen_function%pw_grid%have_g0) THEN
            screen_function%array(1) = screen_function%array(1) + fourpi/(4.0e0_dp*alpha2)
         END IF
      CASE (mt2d)
         iz = special_dimension ! iz is the direction with NO PBC
         zlength = slab_size ! zlength is the thickness of the cell
         DO ig = screen_function%pw_grid%first_gne0, screen_function%pw_grid%ngpts_cut_local
            gz = screen_function%pw_grid%g(iz, ig)
            g2 = screen_function%pw_grid%gsq(ig)
            gxy = sqrt(abs(g2 - gz*gz))
            g3d = fourpi/g2
            screen_function%array(ig) = -g3d*cos(gz*zlength/2.0_dp)*exp(-gxy*zlength/2.0_dp)
         END DO
         IF (screen_function%pw_grid%have_g0) screen_function%array(1) = pi*zlength*zlength/2.0_dp
      CASE (mt1d)
         iz = special_dimension ! iz is the direction with PBC
         CALL mt1din(screen_function)
         cpabort("MT1D unimplemented")
      END SELECT
      CALL pw_pool_release(pw_pool_aux)
      CALL timestop(handle)
   END SUBROUTINE mtin_create_screen_fn
 
! **************************************************************************************************
!> \brief Calculates the Tuckerman Green's function in reciprocal space
!>      according the scheme published on:
!>      Martyna and Tuckerman, J. Chem. Phys. Vol. 110, No. 6, 2810-2821
!> \param Vloc ...
!> \param alpha ...
!> \author Teodoro Laino (09.03.2005)
! **************************************************************************************************
   SUBROUTINE mt0din(Vloc, alpha)
      TYPE(pw_r3d_rs_type), POINTER                      :: vloc
      REAL(kind=dp), INTENT(in)                          :: alpha
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'mt0din'
 
      INTEGER                                            :: handle, i, ii, j, jj, k, kk
      INTEGER, DIMENSION(:), POINTER                     :: glb
      INTEGER, DIMENSION(:, :), POINTER                  :: bo
      REAL(kind=dp)                                      :: dx, dy, dz, fact, omega, r, r2, x, y, &
                                                            y2, z, z2
      REAL(kind=dp), DIMENSION(3)                        :: box, box2
      TYPE(pw_grid_type), POINTER                        :: grid
 
      CALL timeset(routinen, handle)
 
      grid => vloc%pw_grid
      bo => grid%bounds_local
      glb => grid%bounds(1, :)
      vloc%array = 0.0_dp
      box = real(grid%npts, kind=dp)*grid%dr
      box2 = box/2.0_dp
      omega = product(box)
      fact = omega
      dx = grid%dr(1)
      dy = grid%dr(2)
      dz = grid%dr(3)
      kk = bo(1, 3)
      DO k = bo(1, 3), bo(2, 3)
         z = real(k - glb(3), dp)*dz; IF (z > box2(3)) z = box(3) - z
         z2 = z*z
         jj = bo(1, 2)
         DO j = bo(1, 2), bo(2, 2)
            y = real(j - glb(2), dp)*dy; IF (y > box2(2)) y = box(2) - y
            y2 = y*y
            ii = bo(1, 1)
            DO i = bo(1, 1), bo(2, 1)
               x = real(i - glb(1), dp)*dx; IF (x > box2(1)) x = box(1) - x
               r2 = x*x + y2 + z2
               r = sqrt(r2)
               IF (r > 1.0e-10_dp) THEN
                  vloc%array(ii, jj, kk) = erf(alpha*r)/r*fact
               ELSE
                  vloc%array(ii, jj, kk) = 2.0_dp*alpha*oorootpi*fact
               END IF
               ii = ii + 1
            END DO
            jj = jj + 1
         END DO
         kk = kk + 1
      END DO
      CALL timestop(handle)
   END SUBROUTINE mt0din
 
! **************************************************************************************************
!> \brief Calculates the Tuckerman Green's function in reciprocal space
!>      according the scheme published on:
!>      Martyna and Tuckerman, J. Chem. Phys. Vol. 121, No. 23, 11949
!> \param screen_function ...
!> \author Teodoro Laino (11.2005)
! **************************************************************************************************
   SUBROUTINE mt1din(screen_function)
      TYPE(pw_c1d_gs_type), POINTER                      :: screen_function
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'mt1din'
 
      INTEGER                                            :: handle
      REAL(kind=dp)                                      :: dx, dy, dz, omega
      REAL(kind=dp), DIMENSION(3)                        :: box, box2
      TYPE(pw_grid_type), POINTER                        :: grid
 
      CALL timeset(routinen, handle)
      grid => screen_function%pw_grid
      box = real(grid%npts, kind=dp)*grid%dr
      box2 = box/2.0_dp
      omega = product(box)
      dx = grid%dr(1)
      dy = grid%dr(2)
      dz = grid%dr(3)
 
      CALL timestop(handle)
   END SUBROUTINE mt1din
 
END MODULE mt_util