LCOV - code coverage report
Current view: top level - src - ewald_environment_types.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 98.1 % 207 203
Test Date: 2025-07-25 12:55:17 Functions: 87.5 % 8 7

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \par History
      10              : !>      JGH FEB-13-2007 : Distributed/replicated realspace grids
      11              : !>      Teodoro Laino [tlaino] - University of Zurich - 12.2007
      12              : !> \author CJM NOV-30-2003
      13              : ! **************************************************************************************************
      14              : MODULE ewald_environment_types
      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_units,                        ONLY: cp_unit_from_cp2k
      20              :    USE input_cp2k_poisson,              ONLY: create_ewald_section
      21              :    USE input_enumeration_types,         ONLY: enum_i2c,&
      22              :                                               enumeration_type
      23              :    USE input_keyword_types,             ONLY: keyword_get,&
      24              :                                               keyword_type
      25              :    USE input_section_types,             ONLY: section_get_keyword,&
      26              :                                               section_release,&
      27              :                                               section_type,&
      28              :                                               section_vals_get_subs_vals,&
      29              :                                               section_vals_release,&
      30              :                                               section_vals_retain,&
      31              :                                               section_vals_type,&
      32              :                                               section_vals_val_get
      33              :    USE kinds,                           ONLY: dp
      34              :    USE mathconstants,                   ONLY: twopi
      35              :    USE mathlib,                         ONLY: det_3x3
      36              :    USE message_passing,                 ONLY: mp_comm_type,&
      37              :                                               mp_para_env_release,&
      38              :                                               mp_para_env_type
      39              :    USE pw_grid_info,                    ONLY: pw_grid_n_for_fft
      40              :    USE pw_poisson_types,                ONLY: do_ewald_ewald,&
      41              :                                               do_ewald_none,&
      42              :                                               do_ewald_pme,&
      43              :                                               do_ewald_spme
      44              : #include "./base/base_uses.f90"
      45              : 
      46              :    IMPLICIT NONE
      47              :    PRIVATE
      48              : 
      49              : ! **************************************************************************************************
      50              : !> \brief to build arrays of pointers
      51              : !> \param ewald_env the pointer to the ewald_env
      52              : !> \par History
      53              : !>      11/03
      54              : !> \author CJM
      55              : ! **************************************************************************************************
      56              :    TYPE ewald_environment_type
      57              :       PRIVATE
      58              :       LOGICAL   :: do_multipoles = .FALSE. ! Flag for using the multipole code
      59              :       INTEGER   :: do_ipol = -1 ! Solver for induced dipoles
      60              :       INTEGER   :: max_multipole = -1 ! max expansion in the multipoles
      61              :       INTEGER   :: max_ipol_iter = -1 ! max number of interaction for induced dipoles
      62              :       INTEGER   :: ewald_type = -1 ! type of ewald
      63              :       INTEGER   :: gmax(3) = -1 ! max Miller index
      64              :       INTEGER   :: ns_max = -1 ! # grid points for small grid (PME)
      65              :       INTEGER   :: o_spline = -1 ! order of spline (SPME)
      66              :       REAL(KIND=dp) :: precs = 0.0_dp ! precision achieved when evaluating the real-space part
      67              :       REAL(KIND=dp) :: alpha = 0.0_dp, rcut = 0.0_dp ! ewald alpha and real-space cutoff
      68              :       REAL(KIND=dp) :: epsilon = 0.0_dp ! tolerance for small grid (PME)
      69              :       REAL(KIND=dp) :: eps_pol = 0.0_dp ! tolerance for convergence of induced dipoles
      70              :       TYPE(mp_para_env_type), POINTER          :: para_env => NULL()
      71              :       TYPE(section_vals_type), POINTER         :: poisson_section => NULL()
      72              :       ! interaction cutoff is required to make the electrostatic interaction
      73              :       ! continuous at a pair distance equal to rcut. this is ignored by the
      74              :       ! multipole code and is otherwise only active when SHIFT_CUTOFF is used.
      75              :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: interaction_cutoffs => NULL()
      76              :       ! store current cell, used to rebuild lazily.
      77              :       REAL(KIND=dp), DIMENSION(3, 3)          :: cell_hmat = -1.0_dp
      78              :    END TYPE ewald_environment_type
      79              : 
      80              : ! *** Public data types ***
      81              :    PUBLIC :: ewald_environment_type
      82              : 
      83              : ! *** Public subroutines ***
      84              :    PUBLIC :: ewald_env_get, &
      85              :              ewald_env_set, &
      86              :              ewald_env_create, &
      87              :              ewald_env_release, &
      88              :              read_ewald_section, &
      89              :              read_ewald_section_tb
      90              : 
      91              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewald_environment_types'
      92              : 
      93              : CONTAINS
      94              : 
      95              : ! **************************************************************************************************
      96              : !> \brief Purpose: Get the EWALD environment.
      97              : !> \param ewald_env the pointer to the ewald_env
      98              : !> \param ewald_type ...
      99              : !> \param alpha ...
     100              : !> \param eps_pol ...
     101              : !> \param epsilon ...
     102              : !> \param gmax ...
     103              : !> \param ns_max ...
     104              : !> \param o_spline ...
     105              : !> \param group ...
     106              : !> \param para_env ...
     107              : !> \param poisson_section ...
     108              : !> \param precs ...
     109              : !> \param rcut ...
     110              : !> \param do_multipoles ...
     111              : !> \param max_multipole ...
     112              : !> \param do_ipol ...
     113              : !> \param max_ipol_iter ...
     114              : !> \param interaction_cutoffs ...
     115              : !> \param cell_hmat ...
     116              : !> \par History
     117              : !>      11/03
     118              : !> \author CJM
     119              : ! **************************************************************************************************
     120       535506 :    SUBROUTINE ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, &
     121              :                             gmax, ns_max, o_spline, group, para_env, poisson_section, precs, &
     122              :                             rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, &
     123              :                             interaction_cutoffs, cell_hmat)
     124              :       TYPE(ewald_environment_type), INTENT(IN)           :: ewald_env
     125              :       INTEGER, OPTIONAL                                  :: ewald_type
     126              :       REAL(KIND=dp), OPTIONAL                            :: alpha, eps_pol, epsilon
     127              :       INTEGER, OPTIONAL                                  :: gmax(3), ns_max, o_spline
     128              :       TYPE(mp_comm_type), INTENT(OUT), OPTIONAL          :: group
     129              :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env
     130              :       TYPE(section_vals_type), OPTIONAL, POINTER         :: poisson_section
     131              :       REAL(KIND=dp), OPTIONAL                            :: precs, rcut
     132              :       LOGICAL, INTENT(OUT), OPTIONAL                     :: do_multipoles
     133              :       INTEGER, INTENT(OUT), OPTIONAL                     :: max_multipole, do_ipol, max_ipol_iter
     134              :       REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
     135              :          POINTER                                         :: interaction_cutoffs
     136              :       REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL           :: cell_hmat
     137              : 
     138       535506 :       IF (PRESENT(ewald_type)) ewald_type = ewald_env%ewald_type
     139       535506 :       IF (PRESENT(do_multipoles)) do_multipoles = ewald_env%do_multipoles
     140       535506 :       IF (PRESENT(do_ipol)) do_ipol = ewald_env%do_ipol
     141       535506 :       IF (PRESENT(max_multipole)) max_multipole = ewald_env%max_multipole
     142       535506 :       IF (PRESENT(max_ipol_iter)) max_ipol_iter = ewald_env%max_ipol_iter
     143       535506 :       IF (PRESENT(alpha)) alpha = ewald_env%alpha
     144       535506 :       IF (PRESENT(precs)) precs = ewald_env%precs
     145       535506 :       IF (PRESENT(rcut)) rcut = ewald_env%rcut
     146       535506 :       IF (PRESENT(epsilon)) epsilon = ewald_env%epsilon
     147       535506 :       IF (PRESENT(eps_pol)) eps_pol = ewald_env%eps_pol
     148       552622 :       IF (PRESENT(gmax)) gmax = ewald_env%gmax
     149       535506 :       IF (PRESENT(ns_max)) ns_max = ewald_env%ns_max
     150       535506 :       IF (PRESENT(o_spline)) o_spline = ewald_env%o_spline
     151       535506 :       IF (PRESENT(group)) group = ewald_env%para_env
     152       535506 :       IF (PRESENT(para_env)) para_env => ewald_env%para_env
     153       535506 :       IF (PRESENT(poisson_section)) poisson_section => ewald_env%poisson_section
     154       535506 :       IF (PRESENT(interaction_cutoffs)) interaction_cutoffs => &
     155        79000 :          ewald_env%interaction_cutoffs
     156      1613765 :       IF (PRESENT(cell_hmat)) cell_hmat = ewald_env%cell_hmat
     157       535506 :    END SUBROUTINE ewald_env_get
     158              : 
     159              : ! **************************************************************************************************
     160              : !> \brief Purpose: Set the EWALD environment.
     161              : !> \param ewald_env the pointer to the ewald_env
     162              : !> \param ewald_type ...
     163              : !> \param alpha ...
     164              : !> \param epsilon ...
     165              : !> \param eps_pol ...
     166              : !> \param gmax ...
     167              : !> \param ns_max ...
     168              : !> \param precs ...
     169              : !> \param o_spline ...
     170              : !> \param para_env ...
     171              : !> \param poisson_section ...
     172              : !> \param interaction_cutoffs ...
     173              : !> \param cell_hmat ...
     174              : !> \par History
     175              : !>      11/03
     176              : !> \author CJM
     177              : ! **************************************************************************************************
     178        23020 :    SUBROUTINE ewald_env_set(ewald_env, ewald_type, alpha, epsilon, eps_pol, &
     179              :                             gmax, ns_max, precs, o_spline, para_env, poisson_section, &
     180              :                             interaction_cutoffs, cell_hmat)
     181              : 
     182              :       TYPE(ewald_environment_type), INTENT(INOUT)        :: ewald_env
     183              :       INTEGER, OPTIONAL                                  :: ewald_type
     184              :       REAL(KIND=dp), OPTIONAL                            :: alpha, epsilon, eps_pol
     185              :       INTEGER, OPTIONAL                                  :: gmax(3), ns_max
     186              :       REAL(KIND=dp), OPTIONAL                            :: precs
     187              :       INTEGER, OPTIONAL                                  :: o_spline
     188              :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env
     189              :       TYPE(section_vals_type), OPTIONAL, POINTER         :: poisson_section
     190              :       REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
     191              :          POINTER                                         :: interaction_cutoffs
     192              :       REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL           :: cell_hmat
     193              : 
     194        23020 :       IF (PRESENT(ewald_type)) ewald_env%ewald_type = ewald_type
     195        23020 :       IF (PRESENT(alpha)) ewald_env%alpha = alpha
     196        23020 :       IF (PRESENT(precs)) ewald_env%precs = precs
     197        23020 :       IF (PRESENT(epsilon)) ewald_env%epsilon = epsilon
     198        23020 :       IF (PRESENT(eps_pol)) ewald_env%eps_pol = eps_pol
     199        23020 :       IF (PRESENT(gmax)) ewald_env%gmax = gmax
     200        23020 :       IF (PRESENT(ns_max)) ewald_env%ns_max = ns_max
     201        23020 :       IF (PRESENT(o_spline)) ewald_env%o_spline = o_spline
     202        23020 :       IF (PRESENT(para_env)) ewald_env%para_env => para_env
     203        23020 :       IF (PRESENT(poisson_section)) THEN
     204         4271 :          CALL section_vals_retain(poisson_section)
     205         4271 :          CALL section_vals_release(ewald_env%poisson_section)
     206         4271 :          ewald_env%poisson_section => poisson_section
     207              :       END IF
     208        23020 :       IF (PRESENT(interaction_cutoffs)) ewald_env%interaction_cutoffs => &
     209         2643 :          interaction_cutoffs
     210       232398 :       IF (PRESENT(cell_hmat)) ewald_env%cell_hmat = cell_hmat
     211        23020 :    END SUBROUTINE ewald_env_set
     212              : 
     213              : ! **************************************************************************************************
     214              : !> \brief allocates and intitializes a ewald_env
     215              : !> \param ewald_env the object to create
     216              : !> \param para_env ...
     217              : !> \par History
     218              : !>      12.2002 created [fawzi]
     219              : !> \author Fawzi Mohamed
     220              : ! **************************************************************************************************
     221        72607 :    SUBROUTINE ewald_env_create(ewald_env, para_env)
     222              :       TYPE(ewald_environment_type), INTENT(OUT)          :: ewald_env
     223              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     224              : 
     225              :       NULLIFY (ewald_env%poisson_section)
     226         4271 :       CALL para_env%retain()
     227         4271 :       ewald_env%para_env => para_env
     228         4271 :       NULLIFY (ewald_env%interaction_cutoffs) ! allocated and initialized later
     229         4271 :    END SUBROUTINE ewald_env_create
     230              : 
     231              : ! **************************************************************************************************
     232              : !> \brief releases the given ewald_env (see doc/ReferenceCounting.html)
     233              : !> \param ewald_env the object to release
     234              : !> \par History
     235              : !>      12.2002 created [fawzi]
     236              : !> \author Fawzi Mohamed
     237              : ! **************************************************************************************************
     238         4271 :    SUBROUTINE ewald_env_release(ewald_env)
     239              :       TYPE(ewald_environment_type), INTENT(INOUT)        :: ewald_env
     240              : 
     241         4271 :       CALL mp_para_env_release(ewald_env%para_env)
     242         4271 :       CALL section_vals_release(ewald_env%poisson_section)
     243         4271 :       IF (ASSOCIATED(ewald_env%interaction_cutoffs)) THEN
     244         2643 :          DEALLOCATE (ewald_env%interaction_cutoffs)
     245              :       END IF
     246              : 
     247         4271 :    END SUBROUTINE ewald_env_release
     248              : 
     249              : ! **************************************************************************************************
     250              : !> \brief Purpose: read the EWALD section
     251              : !> \param ewald_env the pointer to the ewald_env
     252              : !> \param ewald_section ...
     253              : !> \author Teodoro Laino [tlaino] -University of Zurich - 2005
     254              : ! **************************************************************************************************
     255         3835 :    SUBROUTINE read_ewald_section(ewald_env, ewald_section)
     256              :       TYPE(ewald_environment_type), INTENT(INOUT)        :: ewald_env
     257              :       TYPE(section_vals_type), POINTER                   :: ewald_section
     258              : 
     259              :       INTEGER                                            :: iw
     260         3835 :       INTEGER, DIMENSION(:), POINTER                     :: gmax_read
     261              :       LOGICAL                                            :: explicit
     262              :       REAL(KIND=dp)                                      :: dummy
     263              :       TYPE(cp_logger_type), POINTER                      :: logger
     264              :       TYPE(enumeration_type), POINTER                    :: enum
     265              :       TYPE(keyword_type), POINTER                        :: keyword
     266              :       TYPE(section_type), POINTER                        :: section
     267              :       TYPE(section_vals_type), POINTER                   :: multipole_section
     268              : 
     269         3835 :       NULLIFY (enum, keyword, section, multipole_section)
     270         7670 :       logger => cp_get_default_logger()
     271         3835 :       CALL section_vals_val_get(ewald_section, "EWALD_TYPE", i_val=ewald_env%ewald_type)
     272         3835 :       CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha)
     273         3835 :       CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
     274              : 
     275         3835 :       IF (ewald_env%ewald_type == do_ewald_none) THEN
     276         1050 :          ewald_env%rcut = 0.0_dp
     277              :       ELSE
     278         2785 :          CALL section_vals_val_get(ewald_section, "RCUT", explicit=explicit)
     279         2785 :          IF (explicit) THEN
     280            8 :             CALL section_vals_val_get(ewald_section, "RCUT", r_val=ewald_env%rcut)
     281              :          ELSE
     282         2777 :             ewald_env%rcut = find_ewald_optimal_value(ewald_env%precs)/ewald_env%alpha
     283              :          END IF
     284              :       END IF
     285              :       ! we have no defaults for gmax, gmax is only needed for ewald and spme
     286         6534 :       SELECT CASE (ewald_env%ewald_type)
     287              :       CASE (do_ewald_ewald, do_ewald_spme)
     288         2699 :          CALL section_vals_val_get(ewald_section, "GMAX", i_vals=gmax_read)
     289         5088 :          SELECT CASE (SIZE(gmax_read, 1))
     290              :          CASE (1)
     291         9556 :             ewald_env%gmax = gmax_read(1)
     292              :          CASE (3)
     293         1240 :             ewald_env%gmax = gmax_read
     294              :          CASE DEFAULT
     295         2699 :             CPABORT("")
     296              :          END SELECT
     297         2699 :          IF (ewald_env%ewald_type == do_ewald_spme) THEN
     298         1878 :             CALL section_vals_val_get(ewald_section, "O_SPLINE", i_val=ewald_env%o_spline)
     299              :          END IF
     300              :       CASE (do_ewald_pme)
     301           86 :          CALL section_vals_val_get(ewald_section, "NS_MAX", i_val=ewald_env%ns_max)
     302           86 :          CALL section_vals_val_get(ewald_section, "EPSILON", r_val=ewald_env%epsilon)
     303              :       CASE DEFAULT
     304              :          ! this should not be used for do_ewald_none
     305         4200 :          ewald_env%gmax = HUGE(0)
     306         4885 :          ewald_env%ns_max = HUGE(0)
     307              :       END SELECT
     308              : 
     309              :       ! Multipoles
     310         3835 :       multipole_section => section_vals_get_subs_vals(ewald_section, "MULTIPOLES")
     311         3835 :       CALL section_vals_val_get(multipole_section, "_SECTION_PARAMETERS_", l_val=ewald_env%do_multipoles)
     312         3835 :       CALL section_vals_val_get(multipole_section, "POL_SCF", i_val=ewald_env%do_ipol)
     313         3835 :       CALL section_vals_val_get(multipole_section, "EPS_POL", r_val=ewald_env%eps_pol)
     314         3835 :       IF (ewald_env%do_multipoles) THEN
     315          336 :          SELECT CASE (ewald_env%ewald_type)
     316              :          CASE (do_ewald_ewald)
     317              :             CALL section_vals_val_get(multipole_section, "MAX_MULTIPOLE_EXPANSION", &
     318          168 :                                       i_val=ewald_env%max_multipole)
     319          168 :             CALL section_vals_val_get(multipole_section, "MAX_IPOL_ITER", i_val=ewald_env%max_ipol_iter)
     320              :          CASE DEFAULT
     321          168 :             CPABORT("Multipole code works at the moment only with standard EWALD sums.")
     322              :          END SELECT
     323              :       END IF
     324              : 
     325              :       iw = cp_print_key_unit_nr(logger, ewald_section, "PRINT%PROGRAM_RUN_INFO", &
     326         3835 :                                 extension=".log")
     327         3835 :       IF (iw > 0) THEN
     328         1906 :          NULLIFY (keyword, enum)
     329         1906 :          CALL create_ewald_section(section)
     330         1906 :          IF (ewald_env%ewald_type /= do_ewald_none) THEN
     331         1419 :             keyword => section_get_keyword(section, "EWALD_TYPE")
     332         1419 :             CALL keyword_get(keyword, enum=enum)
     333         1419 :             WRITE (iw, '(/,T2,"EWALD| ",A,T67,A14 )') 'Summation is done by:', &
     334         2838 :                ADJUSTR(TRIM(enum_i2c(enum, ewald_env%ewald_type)))
     335         1419 :             IF (ewald_env%do_multipoles) THEN
     336           84 :                NULLIFY (keyword, enum)
     337           84 :                keyword => section_get_keyword(section, "MULTIPOLES%MAX_MULTIPOLE_EXPANSION")
     338           84 :                CALL keyword_get(keyword, enum=enum)
     339           84 :                WRITE (iw, '( T2,"EWALD| ",A )') 'Enabled Multipole Method'
     340           84 :                WRITE (iw, '( T2,"EWALD| ",A,T67,A14 )') 'Max Term in Multipole Expansion :', &
     341          168 :                   ADJUSTR(TRIM(enum_i2c(enum, ewald_env%max_multipole)))
     342           84 :                WRITE (iw, '( T2,"EWALD| ",A,T67,3I10 )') 'Max number Iterations for IPOL :', &
     343          168 :                   ewald_env%max_ipol_iter
     344              :             END IF
     345         1419 :             dummy = cp_unit_from_cp2k(ewald_env%alpha, "angstrom^-1")
     346              :             WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
     347         1419 :                'Alpha parameter [', 'ANGSTROM^-1', ']', dummy
     348         1419 :             dummy = cp_unit_from_cp2k(ewald_env%rcut, "angstrom")
     349              :             WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
     350         1419 :                'Real Space Cutoff [', 'ANGSTROM', ']', dummy
     351              : 
     352         1869 :             SELECT CASE (ewald_env%ewald_type)
     353              :             CASE (do_ewald_ewald)
     354              :                WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
     355          450 :                   'G-space max. Miller index', ewald_env%gmax
     356              :             CASE (do_ewald_pme)
     357              :                WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
     358           43 :                   'Max small-grid points (input) ', ewald_env%ns_max
     359              :                WRITE (iw, '( T2,"EWALD| ",A,T71,E10.4 )') &
     360           43 :                   'Gaussian tolerance (input) ', ewald_env%epsilon
     361              :             CASE (do_ewald_spme)
     362              :                WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
     363          926 :                   'G-space max. Miller index', ewald_env%gmax
     364              :                WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
     365          926 :                   'Spline interpolation order ', ewald_env%o_spline
     366              :             CASE DEFAULT
     367         1419 :                CPABORT("")
     368              :             END SELECT
     369              :          ELSE
     370          487 :             WRITE (iw, '( T2,"EWALD| ",T73, A )') 'not used'
     371              :          END IF
     372         1906 :          CALL section_release(section)
     373              :       END IF
     374              :       CALL cp_print_key_finished_output(iw, logger, ewald_section, &
     375         3835 :                                         "PRINT%PROGRAM_RUN_INFO")
     376              : 
     377         3835 :    END SUBROUTINE read_ewald_section
     378              : 
     379              : ! **************************************************************************************************
     380              : !> \brief Purpose: read the EWALD section for TB methods
     381              : !> \param ewald_env the pointer to the ewald_env
     382              : !> \param ewald_section ...
     383              : !> \param hmat ...
     384              : !> \param silent ...
     385              : !> \param pset ...
     386              : !> \author JGH
     387              : ! **************************************************************************************************
     388         3052 :    SUBROUTINE read_ewald_section_tb(ewald_env, ewald_section, hmat, silent, pset)
     389              :       TYPE(ewald_environment_type), INTENT(INOUT)        :: ewald_env
     390              :       TYPE(section_vals_type), POINTER                   :: ewald_section
     391              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat
     392              :       LOGICAL, INTENT(IN), OPTIONAL                      :: silent
     393              :       CHARACTER(LEN=*), OPTIONAL                         :: pset
     394              : 
     395              :       CHARACTER(LEN=5)                                   :: param
     396              :       INTEGER                                            :: i, iw, n(3)
     397          436 :       INTEGER, DIMENSION(:), POINTER                     :: gmax_read
     398              :       LOGICAL                                            :: do_print, explicit
     399              :       REAL(KIND=dp)                                      :: alat, cutoff, dummy, omega
     400              :       TYPE(cp_logger_type), POINTER                      :: logger
     401              : 
     402          872 :       logger => cp_get_default_logger()
     403          436 :       do_print = .TRUE.
     404          436 :       IF (PRESENT(silent)) do_print = .NOT. silent
     405          436 :       param = "none"
     406          436 :       IF (PRESENT(pset)) param = pset
     407              : 
     408          436 :       ewald_env%do_multipoles = .FALSE.
     409          436 :       ewald_env%do_ipol = 0
     410          436 :       ewald_env%eps_pol = 1.e-12_dp
     411          436 :       ewald_env%max_multipole = 0
     412          436 :       ewald_env%max_ipol_iter = 0
     413          436 :       ewald_env%epsilon = 1.e-12_dp
     414          436 :       ewald_env%ns_max = HUGE(0)
     415              : 
     416          436 :       CALL section_vals_val_get(ewald_section, "EWALD_TYPE", explicit=explicit)
     417          436 :       IF (explicit) THEN
     418          160 :          CALL section_vals_val_get(ewald_section, "EWALD_TYPE", i_val=ewald_env%ewald_type)
     419          160 :          IF (ewald_env%ewald_type /= do_ewald_spme) THEN
     420            0 :             CPABORT("TB needs EWALD_TYPE SPME")
     421              :          END IF
     422              :       ELSE
     423          276 :          ewald_env%ewald_type = do_ewald_spme
     424              :       END IF
     425              : 
     426          436 :       CALL section_vals_val_get(ewald_section, "ALPHA", explicit=explicit)
     427          436 :       IF (explicit) THEN
     428          158 :          CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha)
     429              :       ELSE
     430          198 :          SELECT CASE (param)
     431              :          CASE DEFAULT
     432          198 :             ewald_env%alpha = 1.0_dp
     433              :          CASE ("EEQ")
     434           80 :             omega = ABS(det_3x3(hmat))
     435          278 :             ewald_env%alpha = SQRT(twopi)/omega**(1./3.)
     436              :          END SELECT
     437              :       END IF
     438              : 
     439          436 :       CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", explicit=explicit)
     440          436 :       IF (explicit) THEN
     441            0 :          CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
     442              :       ELSE
     443          436 :          CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs)
     444              :       END IF
     445              : 
     446          436 :       CALL section_vals_val_get(ewald_section, "O_SPLINE", explicit=explicit)
     447          436 :       IF (explicit) THEN
     448          108 :          CALL section_vals_val_get(ewald_section, "O_SPLINE", i_val=ewald_env%o_spline)
     449              :       ELSE
     450          248 :          SELECT CASE (param)
     451              :          CASE DEFAULT
     452          248 :             ewald_env%o_spline = 6
     453              :          CASE ("EEQ")
     454          328 :             ewald_env%o_spline = 4
     455              :          END SELECT
     456              :       END IF
     457              : 
     458          436 :       CALL section_vals_val_get(ewald_section, "RCUT", explicit=explicit)
     459          436 :       IF (explicit) THEN
     460            0 :          CALL section_vals_val_get(ewald_section, "RCUT", r_val=ewald_env%rcut)
     461              :       ELSE
     462          436 :          ewald_env%rcut = find_ewald_optimal_value(ewald_env%precs)/ewald_env%alpha
     463              :       END IF
     464              : 
     465          436 :       CALL section_vals_val_get(ewald_section, "GMAX", explicit=explicit)
     466          436 :       IF (explicit) THEN
     467          156 :          CALL section_vals_val_get(ewald_section, "GMAX", i_vals=gmax_read)
     468          254 :          SELECT CASE (SIZE(gmax_read, 1))
     469              :          CASE (1)
     470          392 :             ewald_env%gmax = gmax_read(1)
     471              :          CASE (3)
     472          232 :             ewald_env%gmax = gmax_read
     473              :          CASE DEFAULT
     474          156 :             CPABORT("")
     475              :          END SELECT
     476              :       ELSE
     477          200 :          SELECT CASE (param)
     478              :          CASE DEFAULT
     479              :             ! set GMAX using ECUT=alpha*45 Ry
     480          200 :             cutoff = 45._dp*ewald_env%alpha
     481              :          CASE ("EEQ")
     482              :             ! set GMAX using ECUT=alpha*45 Ry
     483          280 :             cutoff = 30._dp*ewald_env%alpha
     484              :          END SELECT
     485         1120 :          DO i = 1, 3
     486         3360 :             alat = SUM(hmat(:, i)**2)
     487          840 :             CPASSERT(alat /= 0._dp)
     488         1120 :             ewald_env%gmax(i) = 2*FLOOR(SQRT(2.0_dp*cutoff*alat)/twopi) + 1
     489              :          END DO
     490              :       END IF
     491         1744 :       n = ewald_env%gmax
     492         1744 :       ewald_env%gmax = pw_grid_n_for_fft(n, odd=.TRUE.)
     493              : 
     494              :       iw = cp_print_key_unit_nr(logger, ewald_section, "PRINT%PROGRAM_RUN_INFO", &
     495          436 :                                 extension=".log")
     496          436 :       IF (iw > 0 .AND. do_print) THEN
     497          193 :          WRITE (iw, '(/,T2,"EWALD| ",A,T67,A14 )') 'Summation is done by:', ADJUSTR("SPME")
     498          193 :          dummy = cp_unit_from_cp2k(ewald_env%alpha, "angstrom^-1")
     499              :          WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
     500          193 :             'Alpha parameter [', 'ANGSTROM^-1', ']', dummy
     501          193 :          dummy = cp_unit_from_cp2k(ewald_env%rcut, "angstrom")
     502              :          WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') &
     503          193 :             'Real Space Cutoff [', 'ANGSTROM', ']', dummy
     504              :          WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') &
     505          193 :             'G-space max. Miller index', ewald_env%gmax
     506              :          WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') &
     507          193 :             'Spline interpolation order ', ewald_env%o_spline
     508              :       END IF
     509              :       CALL cp_print_key_finished_output(iw, logger, ewald_section, &
     510          436 :                                         "PRINT%PROGRAM_RUN_INFO")
     511              : 
     512          436 :    END SUBROUTINE read_ewald_section_tb
     513              : 
     514              : ! **************************************************************************************************
     515              : !> \brief triggers (by bisection) the optimal value for EWALD parameter x
     516              : !>      EXP(-x^2)/x^2 = EWALD_ACCURACY
     517              : !> \param precs ...
     518              : !> \return ...
     519              : !> \author Teodoro Laino [tlaino] - University of Zurich - 12.2007
     520              : ! **************************************************************************************************
     521         3213 :    FUNCTION find_ewald_optimal_value(precs) RESULT(value)
     522              :       REAL(KIND=dp)                                      :: precs, value
     523              : 
     524              :       REAL(KIND=dp)                                      :: func, func1, func2, s, s1, s2
     525              : 
     526         3213 :       s = 0.1_dp
     527         3213 :       func = EXP(-s**2)/s**2 - precs
     528         3213 :       CPASSERT(func > 0.0_dp)
     529       105487 :       DO WHILE (func > 0.0_dp)
     530       102274 :          s = s + 0.1_dp
     531       105487 :          func = EXP(-s**2)/s**2 - precs
     532              :       END DO
     533         3213 :       s2 = s
     534         3213 :       s1 = s - 0.1_dp
     535              :       ! Start bisection
     536              :       DO WHILE (.TRUE.)
     537        80353 :          func2 = EXP(-s2**2)/s2**2 - precs
     538        80353 :          func1 = EXP(-s1**2)/s1**2 - precs
     539        80353 :          CPASSERT(func1 >= 0)
     540        80353 :          CPASSERT(func2 <= 0)
     541        80353 :          s = 0.5_dp*(s1 + s2)
     542        80353 :          func = EXP(-s**2)/s**2 - precs
     543        80353 :          IF (func > 0.0_dp) THEN
     544              :             s1 = s
     545        29166 :          ELSE IF (func < 0.0_dp) THEN
     546        29166 :             s2 = s
     547              :          END IF
     548        80353 :          IF (ABS(func) < 100.0_dp*EPSILON(0.0_dp)) EXIT
     549              :       END DO
     550         3213 :       value = s
     551         3213 :    END FUNCTION find_ewald_optimal_value
     552              : 
     553            0 : END MODULE ewald_environment_types
     554              : 
        

Generated by: LCOV version 2.0-1