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 pw_types
      10             : !> \author CJM
      11             : ! **************************************************************************************************
      12             : MODULE ewald_pw_types
      13             :    USE ao_util,                         ONLY: exp_radius
      14             :    USE cell_types,                      ONLY: cell_type
      15             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      16             :                                               cp_logger_type
      17             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      18             :                                               cp_print_key_unit_nr
      19             :    USE cp_realspace_grid_init,          ONLY: init_input_type
      20             :    USE dg_types,                        ONLY: dg_create,&
      21             :                                               dg_release,&
      22             :                                               dg_type
      23             :    USE dgs,                             ONLY: dg_pme_grid_setup
      24             :    USE ewald_environment_types,         ONLY: ewald_env_get,&
      25             :                                               ewald_environment_type
      26             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      27             :                                               section_vals_type
      28             :    USE kinds,                           ONLY: dp
      29             :    USE mathconstants,                   ONLY: pi
      30             :    USE message_passing,                 ONLY: mp_comm_self,&
      31             :                                               mp_para_env_type
      32             :    USE pw_grid_types,                   ONLY: HALFSPACE,&
      33             :                                               pw_grid_type
      34             :    USE pw_grids,                        ONLY: pw_grid_create,&
      35             :                                               pw_grid_release
      36             :    USE pw_poisson_methods,              ONLY: pw_poisson_set
      37             :    USE pw_poisson_read_input,           ONLY: pw_poisson_read_parameters
      38             :    USE pw_poisson_types,                ONLY: do_ewald_ewald,&
      39             :                                               do_ewald_none,&
      40             :                                               do_ewald_pme,&
      41             :                                               do_ewald_spme,&
      42             :                                               pw_poisson_parameter_type,&
      43             :                                               pw_poisson_type
      44             :    USE pw_pool_types,                   ONLY: pw_pool_create,&
      45             :                                               pw_pool_p_type,&
      46             :                                               pw_pool_release,&
      47             :                                               pw_pool_type
      48             :    USE realspace_grid_types,            ONLY: &
      49             :         realspace_grid_desc_type, realspace_grid_input_type, realspace_grid_type, rs_grid_create, &
      50             :         rs_grid_create_descriptor, rs_grid_print, rs_grid_release, rs_grid_release_descriptor, &
      51             :         rs_grid_retain_descriptor
      52             : #include "./base/base_uses.f90"
      53             : 
      54             :    IMPLICIT NONE
      55             : 
      56             :    PRIVATE
      57             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewald_pw_types'
      58             :    PUBLIC :: ewald_pw_type, ewald_pw_release, &
      59             :              ewald_pw_create, &
      60             :              ewald_pw_get, ewald_pw_set
      61             : 
      62             : ! **************************************************************************************************
      63             :    TYPE ewald_pw_type
      64             :       PRIVATE
      65             :       TYPE(pw_pool_type), POINTER       :: pw_small_pool => NULL()
      66             :       TYPE(pw_pool_type), POINTER       :: pw_big_pool => NULL()
      67             :       TYPE(realspace_grid_desc_type), POINTER    :: rs_desc => NULL()
      68             :       TYPE(pw_poisson_type), POINTER    :: poisson_env => NULL()
      69             :       TYPE(dg_type), POINTER            :: dg => NULL()
      70             :    END TYPE ewald_pw_type
      71             : 
      72             : CONTAINS
      73             : 
      74             : ! **************************************************************************************************
      75             : !> \brief creates the structure ewald_pw_type
      76             : !> \param ewald_pw ...
      77             : !> \param ewald_env ...
      78             : !> \param cell ...
      79             : !> \param cell_ref ...
      80             : !> \param print_section ...
      81             : ! **************************************************************************************************
      82        4169 :    SUBROUTINE ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section)
      83             :       TYPE(ewald_pw_type), INTENT(OUT)                   :: ewald_pw
      84             :       TYPE(ewald_environment_type), POINTER              :: ewald_env
      85             :       TYPE(cell_type), POINTER                           :: cell, cell_ref
      86             :       TYPE(section_vals_type), POINTER                   :: print_section
      87             : 
      88             :       NULLIFY (ewald_pw%pw_big_pool)
      89             :       NULLIFY (ewald_pw%pw_small_pool)
      90             :       NULLIFY (ewald_pw%rs_desc)
      91             :       NULLIFY (ewald_pw%poisson_env)
      92        4169 :       ALLOCATE (ewald_pw%dg)
      93        4169 :       CALL dg_create(ewald_pw%dg)
      94        4169 :       CALL ewald_pw_init(ewald_pw, ewald_env, cell, cell_ref, print_section)
      95        4169 :    END SUBROUTINE ewald_pw_create
      96             : 
      97             : ! **************************************************************************************************
      98             : !> \brief releases the memory used by the ewald_pw
      99             : !> \param ewald_pw ...
     100             : ! **************************************************************************************************
     101        4169 :    SUBROUTINE ewald_pw_release(ewald_pw)
     102             :       TYPE(ewald_pw_type), INTENT(INOUT)                 :: ewald_pw
     103             : 
     104        4169 :       CALL pw_pool_release(ewald_pw%pw_small_pool)
     105        4169 :       CALL pw_pool_release(ewald_pw%pw_big_pool)
     106        4169 :       CALL rs_grid_release_descriptor(ewald_pw%rs_desc)
     107        4169 :       IF (ASSOCIATED(ewald_pw%poisson_env)) THEN
     108        2302 :          CALL ewald_pw%poisson_env%release()
     109        2302 :          DEALLOCATE (ewald_pw%poisson_env)
     110             :       END IF
     111        4169 :       CALL dg_release(ewald_pw%dg)
     112        4169 :       DEALLOCATE (ewald_pw%dg)
     113             : 
     114        4169 :    END SUBROUTINE ewald_pw_release
     115             : 
     116             : ! **************************************************************************************************
     117             : !> \brief ...
     118             : !> \param ewald_pw ...
     119             : !> \param ewald_env ...
     120             : !> \param cell ...
     121             : !> \param cell_ref ...
     122             : !> \param print_section ...
     123             : !> \par History
     124             : !>      JGH (12-Jan-2001): Added SPME part
     125             : !>      JGH (15-Mar-2001): Work newly distributed between initialize, setup,
     126             : !>                         and force routine
     127             : !> \author CJM
     128             : ! **************************************************************************************************
     129        4169 :    SUBROUTINE ewald_pw_init(ewald_pw, ewald_env, cell, cell_ref, print_section)
     130             :       TYPE(ewald_pw_type), INTENT(INOUT)                 :: ewald_pw
     131             :       TYPE(ewald_environment_type), POINTER              :: ewald_env
     132             :       TYPE(cell_type), POINTER                           :: cell, cell_ref
     133             :       TYPE(section_vals_type), POINTER                   :: print_section
     134             : 
     135             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ewald_pw_init'
     136             : 
     137             :       INTEGER                                            :: bo(2, 3), ewald_type, gmax(3), handle, &
     138             :                                                             npts_s(3), ns_max, o_spline, &
     139             :                                                             output_unit
     140             :       REAL(KIND=dp)                                      :: alpha, alphasq, cutoff_radius, epsilon, &
     141             :                                                             norm
     142             :       TYPE(cp_logger_type), POINTER                      :: logger
     143             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     144             :       TYPE(pw_grid_type), POINTER                        :: pw_big_grid, pw_small_grid
     145             :       TYPE(pw_poisson_parameter_type)                    :: poisson_params
     146        4169 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
     147             :       TYPE(pw_pool_type), POINTER                        :: pw_pool
     148             :       TYPE(realspace_grid_desc_type), POINTER            :: rs_desc
     149             :       TYPE(realspace_grid_input_type)                    :: input_settings
     150             :       TYPE(section_vals_type), POINTER                   :: poisson_section, rs_grid_section
     151             : 
     152        4169 :       CALL timeset(routineN, handle)
     153             : 
     154        4169 :       NULLIFY (pw_big_grid)
     155        4169 :       NULLIFY (pw_small_grid, poisson_section)
     156             : 
     157        4169 :       CPASSERT(ASSOCIATED(ewald_env))
     158        4169 :       CPASSERT(ASSOCIATED(cell))
     159             :       CALL ewald_env_get(ewald_env=ewald_env, &
     160             :                          para_env=para_env, &
     161             :                          gmax=gmax, alpha=alpha, &
     162             :                          ns_max=ns_max, &
     163             :                          ewald_type=ewald_type, &
     164             :                          o_spline=o_spline, &
     165             :                          poisson_section=poisson_section, &
     166        4169 :                          epsilon=epsilon)
     167             : 
     168        4169 :       rs_grid_section => section_vals_get_subs_vals(poisson_section, "EWALD%RS_GRID")
     169             : 
     170         821 :       SELECT CASE (ewald_type)
     171             :       CASE (do_ewald_ewald)
     172             :          ! set up Classic EWALD sum
     173         821 :          logger => cp_get_default_logger()
     174         821 :          output_unit = cp_print_key_unit_nr(logger, print_section, "", extension=".Log")
     175             : 
     176        3284 :          IF (ANY(gmax == 2*(gmax/2))) THEN
     177           0 :             CPABORT("gmax has to be odd.")
     178             :          END IF
     179        3284 :          bo(1, :) = -gmax/2
     180        3284 :          bo(2, :) = +gmax/2
     181             :          CALL pw_grid_create(pw_big_grid, mp_comm_self, cell_ref%hmat, grid_span=HALFSPACE, bounds=bo, spherical=.TRUE., &
     182         821 :                              fft_usage=.FALSE., iounit=output_unit)
     183         821 :          NULLIFY (pw_pool)
     184         821 :          CALL pw_pool_create(pw_pool, pw_grid=pw_big_grid)
     185         821 :          ewald_pw%pw_big_pool => pw_pool
     186         821 :          CALL pw_grid_release(pw_big_grid)
     187         821 :          CALL cp_print_key_finished_output(output_unit, logger, print_section, "")
     188             : 
     189             :       CASE (do_ewald_pme)
     190             :          ! set up Particle-Mesh EWALD sum
     191          86 :          logger => cp_get_default_logger()
     192          86 :          output_unit = cp_print_key_unit_nr(logger, print_section, "", extension=".Log")
     193          86 :          IF (.NOT. ASSOCIATED(ewald_pw%poisson_env)) THEN
     194        1376 :             ALLOCATE (ewald_pw%poisson_env)
     195          86 :             CALL ewald_pw%poisson_env%create()
     196             :          END IF
     197          86 :          IF (ns_max == 2*(ns_max/2)) THEN
     198           0 :             CPABORT("ns_max has to be odd.")
     199             :          END IF
     200         344 :          npts_s(:) = ns_max
     201             :          ! compute cut-off radius
     202          86 :          alphasq = alpha**2
     203          86 :          norm = (2.0_dp*alphasq/pi)**(1.5_dp)
     204          86 :          cutoff_radius = exp_radius(0, 2.0_dp*alphasq, epsilon, norm)
     205             : 
     206             :          CALL dg_pme_grid_setup(cell_ref%hmat, npts_s, cutoff_radius, &
     207             :                                 pw_small_grid, pw_big_grid, para_env, rs_dims=(/para_env%num_pe, 1/), &
     208         258 :                                 iounit=output_unit, fft_usage=.TRUE.)
     209             :          ! Write some useful info
     210          86 :          IF (output_unit > 0) THEN
     211             :             WRITE (output_unit, '( A,T71,E10.4 )') &
     212          43 :                ' EWALD| Gaussian tolerance (effective) ', epsilon
     213             :             WRITE (output_unit, '( A,T63,3I6 )') &
     214          43 :                ' EWALD| Small box grid ', pw_small_grid%npts
     215             :             WRITE (output_unit, '( A,T63,3I6 )') &
     216          43 :                ' EWALD| Full box grid ', pw_big_grid%npts
     217             :          END IF
     218             : 
     219             :          ! pw pools initialized
     220          86 :          NULLIFY (pw_pool)
     221          86 :          CALL pw_pool_create(pw_pool, pw_grid=pw_big_grid)
     222          86 :          ewald_pw%pw_big_pool => pw_pool
     223             : 
     224          86 :          NULLIFY (pw_pool)
     225          86 :          CALL pw_pool_create(pw_pool, pw_grid=pw_small_grid)
     226          86 :          ewald_pw%pw_small_pool => pw_pool
     227             : 
     228          86 :          NULLIFY (rs_desc)
     229             :          CALL init_input_type(input_settings, nsmax=MAXVAL(pw_small_grid%npts(1:3)), &
     230             :                               rs_grid_section=rs_grid_section, ilevel=1, &
     231         344 :                               higher_grid_layout=(/-1, -1, -1/))
     232          86 :          CALL rs_grid_create_descriptor(rs_desc, pw_big_grid, input_settings)
     233             : 
     234         258 :          BLOCK
     235        1376 :             TYPE(realspace_grid_type) :: rs
     236          86 :             CALL rs_grid_create(rs, rs_desc)
     237          86 :             CALL rs_grid_print(rs, output_unit)
     238         258 :             CALL rs_grid_release(rs)
     239             :          END BLOCK
     240             : 
     241          86 :          CALL cp_print_key_finished_output(output_unit, logger, print_section, "")
     242             : 
     243          86 :          ewald_pw%rs_desc => rs_desc
     244             : 
     245          86 :          CALL rs_grid_retain_descriptor(ewald_pw%rs_desc)
     246          86 :          CALL rs_grid_release_descriptor(rs_desc)
     247             : 
     248          86 :          CALL pw_grid_release(pw_small_grid)
     249          86 :          CALL pw_grid_release(pw_big_grid)
     250             : 
     251             :       CASE (do_ewald_spme)
     252             :          ! set up the Smooth-Particle-Mesh EWALD sum
     253        2216 :          logger => cp_get_default_logger()
     254        2216 :          output_unit = cp_print_key_unit_nr(logger, print_section, "", extension=".Log")
     255        2216 :          IF (.NOT. ASSOCIATED(ewald_pw%poisson_env)) THEN
     256       35456 :             ALLOCATE (ewald_pw%poisson_env)
     257        2216 :             CALL ewald_pw%poisson_env%create()
     258             :          END IF
     259        2216 :          npts_s = gmax
     260             :          CALL pw_grid_create(pw_big_grid, para_env, cell_ref%hmat, grid_span=HALFSPACE, npts=npts_s, spherical=.TRUE., &
     261        6648 :                              rs_dims=(/para_env%num_pe, 1/), iounit=output_unit, fft_usage=.TRUE.)
     262             : 
     263             :          ! pw pools initialized
     264        2216 :          NULLIFY (pw_pool)
     265        2216 :          CALL pw_pool_create(pw_pool, pw_grid=pw_big_grid)
     266        2216 :          ewald_pw%pw_big_pool => pw_pool
     267             : 
     268        2216 :          NULLIFY (rs_desc)
     269             :          CALL init_input_type(input_settings, nsmax=o_spline, &
     270             :                               rs_grid_section=rs_grid_section, ilevel=1, &
     271        2216 :                               higher_grid_layout=(/-1, -1, -1/))
     272        2216 :          CALL rs_grid_create_descriptor(rs_desc, pw_big_grid, input_settings)
     273             : 
     274        6648 :          BLOCK
     275       35456 :             TYPE(realspace_grid_type) :: rs
     276             : 
     277        2216 :             CALL rs_grid_create(rs, rs_desc)
     278        2216 :             CALL rs_grid_print(rs, output_unit)
     279        6648 :             CALL rs_grid_release(rs)
     280             :          END BLOCK
     281        2216 :          CALL cp_print_key_finished_output(output_unit, logger, print_section, "")
     282             : 
     283        2216 :          ewald_pw%rs_desc => rs_desc
     284             : 
     285        2216 :          CALL rs_grid_retain_descriptor(ewald_pw%rs_desc)
     286        2216 :          CALL rs_grid_release_descriptor(rs_desc)
     287             : 
     288        2216 :          CALL pw_grid_release(pw_big_grid)
     289             :       CASE (do_ewald_none)
     290             :          ! No EWALD sums..
     291             :       CASE default
     292        4169 :          CPABORT("")
     293             :       END SELECT
     294             :       ! Poisson Environment
     295        4169 :       IF (ASSOCIATED(ewald_pw%poisson_env)) THEN
     296        4604 :          ALLOCATE (pw_pools(1))
     297        2302 :          pw_pools(1)%pool => ewald_pw%pw_big_pool
     298        2302 :          CALL pw_poisson_read_parameters(poisson_section, poisson_params)
     299        2302 :          poisson_params%ewald_type = ewald_type
     300        2302 :          poisson_params%ewald_o_spline = o_spline
     301        2302 :          poisson_params%ewald_alpha = alpha
     302             :          CALL pw_poisson_set(ewald_pw%poisson_env, cell_hmat=cell%hmat, parameters=poisson_params, &
     303        2302 :                              use_level=1, pw_pools=pw_pools)
     304        2302 :          DEALLOCATE (pw_pools)
     305             :       END IF
     306        4169 :       CALL timestop(handle)
     307       29183 :    END SUBROUTINE ewald_pw_init
     308             : 
     309             : ! **************************************************************************************************
     310             : !> \brief get the ewald_pw environment to the correct program.
     311             : !> \param ewald_pw ...
     312             : !> \param pw_big_pool ...
     313             : !> \param pw_small_pool ...
     314             : !> \param rs_desc ...
     315             : !> \param poisson_env ...
     316             : !> \param dg ...
     317             : !> \author CJM
     318             : ! **************************************************************************************************
     319       94075 :    SUBROUTINE ewald_pw_get(ewald_pw, pw_big_pool, pw_small_pool, rs_desc, poisson_env, dg)
     320             : 
     321             :       TYPE(ewald_pw_type), INTENT(IN)                    :: ewald_pw
     322             :       TYPE(pw_pool_type), OPTIONAL, POINTER              :: pw_big_pool, pw_small_pool
     323             :       TYPE(realspace_grid_desc_type), OPTIONAL, POINTER  :: rs_desc
     324             :       TYPE(pw_poisson_type), OPTIONAL, POINTER           :: poisson_env
     325             :       TYPE(dg_type), OPTIONAL, POINTER                   :: dg
     326             : 
     327       94075 :       IF (PRESENT(poisson_env)) poisson_env => ewald_pw%poisson_env
     328       94075 :       IF (PRESENT(pw_big_pool)) pw_big_pool => ewald_pw%pw_big_pool
     329       94075 :       IF (PRESENT(pw_small_pool)) pw_small_pool => ewald_pw%pw_small_pool
     330       94075 :       IF (PRESENT(rs_desc)) rs_desc => ewald_pw%rs_desc
     331       94075 :       IF (PRESENT(dg)) dg => ewald_pw%dg
     332             : 
     333       94075 :    END SUBROUTINE ewald_pw_get
     334             : 
     335             : ! **************************************************************************************************
     336             : !> \brief set the ewald_pw environment to the correct program.
     337             : !> \param ewald_pw ...
     338             : !> \param pw_big_pool ...
     339             : !> \param pw_small_pool ...
     340             : !> \param rs_desc ...
     341             : !> \param dg ...
     342             : !> \param poisson_env ...
     343             : !> \author CJM
     344             : ! **************************************************************************************************
     345       10722 :    SUBROUTINE ewald_pw_set(ewald_pw, pw_big_pool, pw_small_pool, rs_desc, dg, &
     346             :                            poisson_env)
     347             : 
     348             :       TYPE(ewald_pw_type), INTENT(INOUT)                 :: ewald_pw
     349             :       TYPE(pw_pool_type), OPTIONAL, POINTER              :: pw_big_pool, pw_small_pool
     350             :       TYPE(realspace_grid_desc_type), OPTIONAL, POINTER  :: rs_desc
     351             :       TYPE(dg_type), OPTIONAL, POINTER                   :: dg
     352             :       TYPE(pw_poisson_type), OPTIONAL, POINTER           :: poisson_env
     353             : 
     354       10722 :       IF (PRESENT(pw_big_pool)) THEN
     355       10722 :          CALL pw_big_pool%retain()
     356       10722 :          CALL pw_pool_release(ewald_pw%pw_big_pool)
     357       10722 :          ewald_pw%pw_big_pool => pw_big_pool
     358             :       END IF
     359       10722 :       IF (PRESENT(pw_small_pool)) THEN
     360          86 :          CALL pw_small_pool%retain()
     361          86 :          CALL pw_pool_release(ewald_pw%pw_small_pool)
     362          86 :          ewald_pw%pw_small_pool => pw_small_pool
     363             :       END IF
     364       10722 :       IF (PRESENT(rs_desc)) THEN
     365           0 :          CALL rs_grid_retain_descriptor(rs_desc)
     366           0 :          CALL rs_grid_release_descriptor(ewald_pw%rs_desc)
     367           0 :          ewald_pw%rs_desc => rs_desc
     368             :       END IF
     369       10722 :       IF (PRESENT(dg)) THEN
     370           0 :          CALL dg_release(ewald_pw%dg)
     371           0 :          ewald_pw%dg => dg
     372             :       END IF
     373       10722 :       IF (PRESENT(poisson_env)) THEN
     374       10722 :          IF (ASSOCIATED(ewald_pw%poisson_env)) THEN
     375        8872 :          IF (.NOT. ASSOCIATED(ewald_pw%poisson_env, poisson_env)) THEN
     376           0 :             CALL ewald_pw%poisson_env%release()
     377           0 :             DEALLOCATE (ewald_pw%poisson_env)
     378             :          END IF
     379             :          END IF
     380       10722 :          ewald_pw%poisson_env => poisson_env
     381             :       END IF
     382             : 
     383       10722 :    END SUBROUTINE ewald_pw_set
     384             : 
     385           0 : END MODULE ewald_pw_types