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