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 Environment storing all data that is needed in order to call the DFT
      10             : !>        driver of the PEXSI library with data from the linear scaling quickstep
      11             : !>        SCF environment, mainly parameters and intermediate data for the matrix
      12             : !>        conversion between DBCSR and CSR format.
      13             : !> \par History
      14             : !>       2014.11 created [Patrick Seewald]
      15             : !> \author Patrick Seewald
      16             : ! **************************************************************************************************
      17             : 
      18             : MODULE pexsi_types
      19             :    USE ISO_C_BINDING,                   ONLY: C_INTPTR_T
      20             :    USE bibliography,                    ONLY: Lin2009,&
      21             :                                               Lin2013,&
      22             :                                               cite_reference
      23             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      24             :                                               cp_logger_get_default_unit_nr,&
      25             :                                               cp_logger_type
      26             :    USE dbcsr_api,                       ONLY: dbcsr_csr_type,&
      27             :                                               dbcsr_p_type,&
      28             :                                               dbcsr_scale,&
      29             :                                               dbcsr_type
      30             :    USE kinds,                           ONLY: dp
      31             :    USE message_passing,                 ONLY: mp_comm_type,&
      32             :                                               mp_dims_create
      33             :    USE pexsi_interface,                 ONLY: cp_pexsi_get_options,&
      34             :                                               cp_pexsi_options,&
      35             :                                               cp_pexsi_plan_finalize,&
      36             :                                               cp_pexsi_plan_initialize,&
      37             :                                               cp_pexsi_set_options
      38             : #include "./base/base_uses.f90"
      39             : 
      40             :    IMPLICIT NONE
      41             : 
      42             :    PRIVATE
      43             : 
      44             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pexsi_types'
      45             : 
      46             :    LOGICAL, PARAMETER, PRIVATE          :: careful_mod = .FALSE.
      47             : 
      48             :    INTEGER, PARAMETER, PUBLIC           :: cp2k_to_pexsi = 1, pexsi_to_cp2k = 2
      49             : 
      50             :    PUBLIC :: lib_pexsi_env, lib_pexsi_init, lib_pexsi_finalize, &
      51             :              convert_nspin_cp2k_pexsi
      52             : 
      53             : ! **************************************************************************************************
      54             : !> \brief All PEXSI related data
      55             : !> \param options                       PEXSI options
      56             : !> \param plan                          PEXSI plan
      57             : !> \param mp_group                      message-passing group ID
      58             : !> \param mp_dims                       dimensions of the MPI cartesian grid used
      59             : !>                                      for PEXSI
      60             : !> \param num_ranks_per_pole            number of MPI ranks per pole in PEXSI
      61             : !> \param kTS                           entropic energy contribution
      62             : !> \param matrix_w                      energy-weighted density matrix as needed
      63             : !>                                      for the forces
      64             : !> \param csr_mat                       intermediate matrices in CSR format
      65             : !> \param dbcsr_template_matrix_sym     Symmetric template matrix fixing DBCSR
      66             : !>                                      sparsity pattern
      67             : !> \param dbcsr_template_matrix_nonsym  Nonsymmetric template matrix fixing
      68             : !>                                      DBCSR sparsity pattern
      69             : !> \param csr_sparsity                  DBCSR matrix defining CSR sparsity
      70             : !> \param csr_screening                 whether distance screening should be
      71             : !>                                      applied to CSR matrices
      72             : !> \param max_ev_vector                 eigenvector corresponding to the largest
      73             : !>                                      energy eigenvalue,
      74             : !>                                      returned by the Arnoldi method used to
      75             : !>                                      determine the spectral radius deltaE
      76             : !> \param nspin                         number of spins
      77             : !> \param do_adaptive_tol_nel           Whether or not to use adaptive threshold
      78             : !>                                      for PEXSI convergence
      79             : !> \param adaptive_nel_alpha            constants for adaptive thresholding
      80             : !> \param adaptive_nel_beta             ...
      81             : !> \param tol_nel_initial               Initial convergence threshold (in number of electrons)
      82             : !> \param tol_nel_target                Target convergence threshold (in number of electrons)
      83             : !> \par History
      84             : !>      11.2014 created [Patrick Seewald]
      85             : !> \author Patrick Seewald
      86             : ! **************************************************************************************************
      87             :    TYPE lib_pexsi_env
      88             :       TYPE(dbcsr_type)                         :: dbcsr_template_matrix_sym, &
      89             :                                                   dbcsr_template_matrix_nonsym
      90             :       TYPE(dbcsr_csr_type)                     :: csr_mat_p, csr_mat_ks, csr_mat_s, &
      91             :                                                   csr_mat_E, csr_mat_F
      92             :       TYPE(cp_pexsi_options)                   :: options
      93             :       REAL(KIND=dp), DIMENSION(:), POINTER     :: kTS => NULL()
      94             :       TYPE(dbcsr_p_type), DIMENSION(:), &
      95             :          POINTER                               :: matrix_w => NULL()
      96             :       INTEGER(KIND=C_INTPTR_T)                 :: plan
      97             :       INTEGER                                  :: nspin, num_ranks_per_pole
      98             :       TYPE(mp_comm_type) :: mp_group
      99             :       TYPE(dbcsr_type), DIMENSION(:), &
     100             :          POINTER                               :: max_ev_vector
     101             :       TYPE(dbcsr_type)                         :: csr_sparsity
     102             :       INTEGER, DIMENSION(2)                    :: mp_dims
     103             : 
     104             :       LOGICAL                                  :: csr_screening, do_adaptive_tol_nel
     105             :       REAL(KIND=dp)                            :: adaptive_nel_alpha, adaptive_nel_beta, &
     106             :                                                   tol_nel_initial, tol_nel_target
     107             :    END TYPE lib_pexsi_env
     108             : 
     109             : CONTAINS
     110             : 
     111             : ! **************************************************************************************************
     112             : !> \brief Initialize PEXSI
     113             : !> \param pexsi_env All data needed by PEXSI
     114             : !> \param mp_group message-passing group ID
     115             : !> \param nspin number of spins
     116             : !> \par History
     117             : !>      11.2014 created [Patrick Seewald]
     118             : !> \author Patrick Seewald
     119             : ! **************************************************************************************************
     120           8 :    SUBROUTINE lib_pexsi_init(pexsi_env, mp_group, nspin)
     121             :       TYPE(lib_pexsi_env), INTENT(INOUT)                 :: pexsi_env
     122             : 
     123             :       CLASS(mp_comm_type), INTENT(IN)                     :: mp_group
     124             :       INTEGER, INTENT(IN)                                :: nspin
     125             : 
     126             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'lib_pexsi_init'
     127             : 
     128             :       INTEGER                                            :: handle, ispin, npSymbFact, &
     129             :                                                             unit_nr
     130             :       TYPE(cp_logger_type), POINTER                      :: logger
     131             : 
     132           8 :       logger => cp_get_default_logger()
     133           8 :       IF (logger%para_env%is_source()) THEN
     134           4 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     135             :       ELSE
     136           4 :          unit_nr = -1
     137             :       END IF
     138             : 
     139           8 :       CALL timeset(routineN, handle)
     140             : 
     141           8 :       CALL cite_reference(Lin2009)
     142           8 :       CALL cite_reference(Lin2013)
     143             : 
     144           8 :       pexsi_env%mp_group = mp_group
     145             :       ASSOCIATE (numnodes => pexsi_env%mp_group%num_pe, mynode => pexsi_env%mp_group%mepos)
     146             : 
     147             :          ! Use all nodes available per pole by default or if the user tries to use
     148             :          ! more nodes than MPI size
     149             :          IF ((pexsi_env%num_ranks_per_pole .GT. numnodes) &
     150           2 :              .OR. (pexsi_env%num_ranks_per_pole .EQ. 0)) THEN
     151           8 :             pexsi_env%num_ranks_per_pole = numnodes
     152             :          END IF
     153             : 
     154             :          ! Find num_ranks_per_pole from user input MIN_RANKS_PER_POLE s.t. num_ranks_per_pole
     155             :          ! is the smallest number greater or equal to min_ranks_per_pole that divides
     156             :          ! MPI size without remainder.
     157           8 :          DO WHILE (MOD(numnodes, pexsi_env%num_ranks_per_pole) .NE. 0)
     158           0 :             pexsi_env%num_ranks_per_pole = pexsi_env%num_ranks_per_pole + 1
     159             :          END DO
     160             : 
     161           8 :          CALL cp_pexsi_get_options(pexsi_env%options, npSymbFact=npSymbFact)
     162           8 :          IF ((npSymbFact .GT. pexsi_env%num_ranks_per_pole) .OR. (npSymbFact .EQ. 0)) THEN
     163             :             ! Use maximum possible number of ranks for symbolic factorization
     164           0 :             CALL cp_pexsi_set_options(pexsi_env%options, npSymbFact=pexsi_env%num_ranks_per_pole)
     165             :          END IF
     166             : 
     167             :          ! Create dimensions for MPI cartesian grid for PEXSI
     168          24 :          pexsi_env%mp_dims = 0
     169           8 :          CALL mp_dims_create(pexsi_env%num_ranks_per_pole, pexsi_env%mp_dims)
     170             : 
     171             :          ! allocations with nspin
     172           8 :          pexsi_env%nspin = nspin
     173          24 :          ALLOCATE (pexsi_env%kTS(nspin))
     174          16 :          pexsi_env%kTS(:) = 0.0_dp
     175          32 :          ALLOCATE (pexsi_env%max_ev_vector(nspin))
     176          32 :          ALLOCATE (pexsi_env%matrix_w(nspin))
     177          16 :          DO ispin = 1, pexsi_env%nspin
     178          16 :             ALLOCATE (pexsi_env%matrix_w(ispin)%matrix)
     179             :          END DO
     180             : 
     181             :          ! Initialize PEXSI
     182             :          pexsi_env%plan = cp_pexsi_plan_initialize(pexsi_env%mp_group, pexsi_env%mp_dims(1), &
     183          24 :                                                    pexsi_env%mp_dims(2), mynode)
     184             :       END ASSOCIATE
     185             : 
     186           8 :       pexsi_env%do_adaptive_tol_nel = .FALSE.
     187             : 
     188             :       ! Print PEXSI infos
     189           8 :       IF (unit_nr > 0) CALL print_pexsi_info(pexsi_env, unit_nr)
     190             : 
     191           8 :       CALL timestop(handle)
     192           8 :    END SUBROUTINE lib_pexsi_init
     193             : 
     194             : ! **************************************************************************************************
     195             : !> \brief Release all PEXSI data
     196             : !> \param pexsi_env ...
     197             : !> \par History
     198             : !>      11.2014 created [Patrick Seewald]
     199             : !> \author Patrick Seewald
     200             : ! **************************************************************************************************
     201           8 :    SUBROUTINE lib_pexsi_finalize(pexsi_env)
     202             :       TYPE(lib_pexsi_env), INTENT(INOUT)                 :: pexsi_env
     203             : 
     204             :       CHARACTER(len=*), PARAMETER :: routineN = 'lib_pexsi_finalize'
     205             : 
     206             :       INTEGER                                            :: handle, ispin
     207             : 
     208           8 :       CALL timeset(routineN, handle)
     209           8 :       CALL cp_pexsi_plan_finalize(pexsi_env%plan)
     210           8 :       DEALLOCATE (pexsi_env%kTS)
     211           8 :       DEALLOCATE (pexsi_env%max_ev_vector)
     212          16 :       DO ispin = 1, pexsi_env%nspin
     213          16 :          DEALLOCATE (pexsi_env%matrix_w(ispin)%matrix)
     214             :       END DO
     215           8 :       DEALLOCATE (pexsi_env%matrix_w)
     216           8 :       CALL timestop(handle)
     217           8 :    END SUBROUTINE lib_pexsi_finalize
     218             : 
     219             : ! **************************************************************************************************
     220             : !> \brief Scale various quantities with factors of 2. This converts spin restricted
     221             : !>        DFT calculations (PEXSI) to the unrestricted case (as is the case where
     222             : !>        the density matrix method is called in the linear scaling code).
     223             : !> \param[in] direction ...
     224             : !> \param[in,out] numElectron ...
     225             : !> \param matrix_p ...
     226             : !> \param matrix_w ...
     227             : !> \param kTS ...
     228             : !> \par History
     229             : !>      01.2015 created [Patrick Seewald]
     230             : !> \author Patrick Seewald
     231             : ! **************************************************************************************************
     232         188 :    SUBROUTINE convert_nspin_cp2k_pexsi(direction, numElectron, matrix_p, matrix_w, kTS)
     233             :       INTEGER, INTENT(IN)                                :: direction
     234             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: numElectron
     235             :       TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL          :: matrix_p
     236             :       TYPE(dbcsr_p_type), INTENT(INOUT), OPTIONAL        :: matrix_w
     237             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: kTS
     238             : 
     239             :       CHARACTER(len=*), PARAMETER :: routineN = 'convert_nspin_cp2k_pexsi'
     240             : 
     241             :       INTEGER                                            :: handle
     242             :       REAL(KIND=dp)                                      :: scaling
     243             : 
     244         188 :       CALL timeset(routineN, handle)
     245             : 
     246         282 :       SELECT CASE (direction)
     247             :       CASE (cp2k_to_pexsi)
     248          94 :          scaling = 2.0_dp
     249             :       CASE (pexsi_to_cp2k)
     250         188 :          scaling = 0.5_dp
     251             :       END SELECT
     252             : 
     253         188 :       IF (PRESENT(numElectron)) numElectron = scaling*numElectron
     254         188 :       IF (PRESENT(matrix_p)) CALL dbcsr_scale(matrix_p, scaling)
     255         188 :       IF (PRESENT(matrix_w)) CALL dbcsr_scale(matrix_w%matrix, scaling)
     256         188 :       IF (PRESENT(kTS)) kTS = scaling*kTS
     257             : 
     258         188 :       CALL timestop(handle)
     259         188 :    END SUBROUTINE convert_nspin_cp2k_pexsi
     260             : 
     261             : ! **************************************************************************************************
     262             : !> \brief Print relevant options of PEXSI
     263             : !> \param pexsi_env ...
     264             : !> \param unit_nr ...
     265             : ! **************************************************************************************************
     266           4 :    SUBROUTINE print_pexsi_info(pexsi_env, unit_nr)
     267             :       TYPE(lib_pexsi_env), INTENT(IN)                    :: pexsi_env
     268             :       INTEGER, INTENT(IN)                                :: unit_nr
     269             : 
     270             :       INTEGER                                            :: npSymbFact, numnodes, numPole, ordering, &
     271             :                                                             rowOrdering
     272             :       REAL(KIND=dp)                                      :: gap, muInertiaExpansion, &
     273             :                                                             muInertiaTolerance, muMax0, muMin0, &
     274             :                                                             muPEXSISafeGuard, temperature
     275             : 
     276           4 :       numnodes = pexsi_env%mp_group%num_pe
     277             : 
     278             :       CALL cp_pexsi_get_options(pexsi_env%options, temperature=temperature, gap=gap, &
     279             :                                 numPole=numPole, muMin0=muMin0, muMax0=muMax0, muInertiaTolerance= &
     280             :                                 muInertiaTolerance, muInertiaExpansion=muInertiaExpansion, &
     281             :                                 muPEXSISafeGuard=muPEXSISafeGuard, ordering=ordering, rowOrdering=rowOrdering, &
     282           4 :                                 npSymbFact=npSymbFact)
     283             : 
     284           4 :       WRITE (unit_nr, '(/A)') " PEXSI| Initialized with parameters"
     285           4 :       WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Electronic temperature", temperature
     286           4 :       WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Spectral gap", gap
     287           4 :       WRITE (unit_nr, '(A,T61,I20)') " PEXSI|   Number of poles", numPole
     288           4 :       WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Target tolerance in number of electrons", &
     289           8 :          pexsi_env%tol_nel_target
     290           4 :       WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Convergence tolerance for inertia counting in mu", &
     291           8 :          muInertiaTolerance
     292           4 :       WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Restart tolerance for inertia counting in mu", &
     293           8 :          muPEXSISafeGuard
     294           4 :       WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Expansion of mu interval for inertia counting", &
     295           8 :          muInertiaExpansion
     296             : 
     297           4 :       WRITE (unit_nr, '(/A)') " PEXSI| Parallelization scheme"
     298           4 :       WRITE (unit_nr, '(A,T61,I20)') " PEXSI|   Number of poles processed in parallel", &
     299           8 :          numnodes/pexsi_env%num_ranks_per_pole
     300           4 :       WRITE (unit_nr, '(A,T61,I20)') " PEXSI|   Number of processors per pole", &
     301           8 :          pexsi_env%num_ranks_per_pole
     302           4 :       WRITE (unit_nr, '(A,T71,I5,T76,I5)') " PEXSI|   Process grid dimensions", &
     303           8 :          pexsi_env%mp_dims(1), pexsi_env%mp_dims(2)
     304           4 :       SELECT CASE (ordering)
     305             :       CASE (0)
     306           4 :          WRITE (unit_nr, '(A,T61,A20)') " PEXSI|   Ordering strategy", "parallel"
     307             :       CASE (1)
     308           0 :          WRITE (unit_nr, '(A,T61,A20)') " PEXSI|   Ordering strategy", "sequential"
     309             :       CASE (2)
     310           4 :          WRITE (unit_nr, '(A,T61,A20)') " PEXSI|   Ordering strategy", "multiple minimum degree"
     311             :       END SELECT
     312           4 :       SELECT CASE (rowOrdering)
     313             :       CASE (0)
     314           4 :          WRITE (unit_nr, '(A,T61,A20)') " PEXSI|   Row permutation strategy", "no row permutation"
     315             :       CASE (1)
     316           4 :          WRITE (unit_nr, '(A,T61,A20)') " PEXSI|   Row permutation strategy", "make diagonal entry larger than off diagonal"
     317             :       END SELECT
     318           4 :       IF (ordering .EQ. 0) WRITE (unit_nr, '(A,T61,I20/)') &
     319           4 :          " PEXSI|   Number of processors for symbolic factorization", npSymbFact
     320             : 
     321           4 :    END SUBROUTINE print_pexsi_info
     322           0 : END MODULE pexsi_types