1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Interface to the PEXSI library, providing wrappers for all PEXSI 8!> routines that are called inside CP2K. Requires PEXSI version 0.10.x. 9!> \par History 10!> 2014.12 created [Patrick Seewald] 11!> \author Patrick Seewald 12! ************************************************************************************************** 13MODULE pexsi_interface 14 15#ifdef __LIBPEXSI 16 USE f_ppexsi_interface, ONLY: f_ppexsi_dft_driver, & 17 f_ppexsi_load_real_hs_matrix, & 18 f_ppexsi_options, & 19 f_ppexsi_plan_finalize, & 20 f_ppexsi_plan_initialize, & 21 f_ppexsi_retrieve_real_dft_matrix, & 22 f_ppexsi_set_default_options 23#endif 24#if defined (__HAS_IEEE_EXCEPTIONS) 25 USE ieee_exceptions, ONLY: ieee_get_halting_mode, & 26 ieee_set_halting_mode, & 27 ieee_all 28#endif 29 USE kinds, ONLY: int_8, & 30 real_8 31#include "./base/base_uses.f90" 32 33 IMPLICIT NONE 34 35 PRIVATE 36 37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pexsi_interface' 38 39 PUBLIC :: cp_pexsi_options, cp_pexsi_plan_initialize, & 40 cp_pexsi_load_real_hs_matrix, cp_pexsi_dft_driver, & 41 cp_pexsi_retrieve_real_dft_matrix, cp_pexsi_plan_finalize, & 42 cp_pexsi_set_options, cp_pexsi_get_options, cp_pexsi_set_default_options 43 44 TYPE cp_pexsi_options 45 PRIVATE 46#ifdef __LIBPEXSI 47 TYPE(f_ppexsi_options) :: options 48#else 49 INTEGER :: unused = -1 50#endif 51 END TYPE cp_pexsi_options 52 53CONTAINS 54 55! ************************************************************************************************** 56!> \brief Set PEXSI internal options 57!> \param pexsi_options ... 58!> \param temperature ... 59!> \param gap ... 60!> \param deltaE ... 61!> \param numPole ... 62!> \param isInertiaCount ... 63!> \param maxPEXSIIter ... 64!> \param muMin0 ... 65!> \param muMax0 ... 66!> \param mu0 ... 67!> \param muInertiaTolerance ... 68!> \param muInertiaExpansion ... 69!> \param muPEXSISafeGuard ... 70!> \param numElectronPEXSITolerance ... 71!> \param matrixType ... 72!> \param isSymbolicFactorize ... 73!> \param ordering ... 74!> \param rowOrdering ... 75!> \param npSymbFact ... 76!> \param verbosity ... 77! ************************************************************************************************** 78 SUBROUTINE cp_pexsi_set_options(pexsi_options, temperature, gap, deltaE, numPole, & 79 isInertiaCount, maxPEXSIIter, muMin0, muMax0, mu0, & 80 muInertiaTolerance, muInertiaExpansion, & 81 muPEXSISafeGuard, numElectronPEXSITolerance, & 82 matrixType, isSymbolicFactorize, ordering, rowOrdering, & 83 npSymbFact, verbosity) 84 85 TYPE(cp_pexsi_options), INTENT(INOUT) :: pexsi_options 86 REAL(KIND=real_8), INTENT(IN), OPTIONAL :: temperature, gap, deltaE 87 INTEGER, INTENT(IN), OPTIONAL :: numPole, isInertiaCount, & 88 maxPEXSIIter 89 REAL(KIND=real_8), INTENT(IN), OPTIONAL :: muMin0, muMax0, mu0, & 90 muInertiaTolerance, muInertiaExpansion, muPEXSISafeGuard, & 91 numElectronPEXSITolerance 92 INTEGER, INTENT(IN), OPTIONAL :: matrixType, & 93 isSymbolicFactorize, & 94 ordering, rowOrdering, npSymbFact, & 95 verbosity 96 97 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_pexsi_set_options', & 98 routineP = moduleN//':'//routineN 99 100#ifdef __LIBPEXSI 101 IF (PRESENT(temperature)) pexsi_options%options%temperature = temperature 102 IF (PRESENT(gap)) pexsi_options%options%gap = gap 103 IF (PRESENT(deltaE)) pexsi_options%options%deltaE = deltaE 104 IF (PRESENT(numPole)) pexsi_options%options%numPole = numPole 105 IF (PRESENT(isInertiaCount)) pexsi_options%options%isInertiaCount = isInertiaCount 106 IF (PRESENT(maxPEXSIIter)) pexsi_options%options%maxPEXSIIter = maxPEXSIIter 107 IF (PRESENT(muMin0)) pexsi_options%options%muMin0 = muMin0 108 IF (PRESENT(muMax0)) pexsi_options%options%muMax0 = muMax0 109 IF (PRESENT(mu0)) pexsi_options%options%mu0 = mu0 110 IF (PRESENT(muInertiaTolerance)) & 111 pexsi_options%options%muInertiaTolerance = muInertiaTolerance 112 IF (PRESENT(muInertiaExpansion)) & 113 pexsi_options%options%muInertiaExpansion = muInertiaExpansion 114 IF (PRESENT(muPEXSISafeGuard)) & 115 pexsi_options%options%muPEXSISafeGuard = muPEXSISafeGuard 116 IF (PRESENT(numElectronPEXSITolerance)) & 117 pexsi_options%options%numElectronPEXSITolerance = numElectronPEXSITolerance 118 IF (PRESENT(matrixType)) pexsi_options%options%matrixType = matrixType 119 IF (PRESENT(isSymbolicFactorize)) & 120 pexsi_options%options%isSymbolicFactorize = isSymbolicFactorize 121 IF (PRESENT(ordering)) pexsi_options%options%ordering = ordering 122 IF (PRESENT(rowOrdering)) pexsi_options%options%rowOrdering = rowOrdering 123 IF (PRESENT(npSymbFact)) pexsi_options%options%npSymbFact = npSymbFact 124 IF (PRESENT(verbosity)) pexsi_options%options%verbosity = verbosity 125#else 126 MARK_USED(pexsi_options) 127 MARK_USED(temperature) 128 MARK_USED(gap) 129 MARK_USED(deltaE) 130 MARK_USED(numPole) 131 MARK_USED(isInertiaCount) 132 MARK_USED(maxPEXSIIter) 133 MARK_USED(muMin0) 134 MARK_USED(muMax0) 135 MARK_USED(mu0) 136 MARK_USED(muInertiaTolerance) 137 MARK_USED(muInertiaExpansion) 138 MARK_USED(muPEXSISafeGuard) 139 MARK_USED(numElectronPEXSITolerance) 140 MARK_USED(matrixType) 141 MARK_USED(isSymbolicFactorize) 142 MARK_USED(ordering) 143 MARK_USED(rowOrdering) 144 MARK_USED(npSymbFact) 145 MARK_USED(verbosity) 146 CPABORT("Requires linking to the PEXSI library.") 147#endif 148 149 ! Additional PEXSI parameters and their defaults not made available here 150 ! because CP2K should always use PEXSI's defaults: 151 ! isConstructCommPattern (=?, pexsi does not even use it) 152 ! symmetric (=1) 153 ! transpose (=0) 154 END SUBROUTINE cp_pexsi_set_options 155 156! ************************************************************************************************** 157!> \brief Access PEXSI internal options 158!> \param pexsi_options ... 159!> \param temperature ... 160!> \param gap ... 161!> \param deltaE ... 162!> \param numPole ... 163!> \param isInertiaCount ... 164!> \param maxPEXSIIter ... 165!> \param muMin0 ... 166!> \param muMax0 ... 167!> \param mu0 ... 168!> \param muInertiaTolerance ... 169!> \param muInertiaExpansion ... 170!> \param muPEXSISafeGuard ... 171!> \param numElectronPEXSITolerance ... 172!> \param matrixType ... 173!> \param isSymbolicFactorize ... 174!> \param ordering ... 175!> \param rowOrdering ... 176!> \param npSymbFact ... 177!> \param verbosity ... 178! ************************************************************************************************** 179 SUBROUTINE cp_pexsi_get_options(pexsi_options, temperature, gap, deltaE, numPole, & 180 isInertiaCount, maxPEXSIIter, muMin0, muMax0, mu0, & 181 muInertiaTolerance, muInertiaExpansion, & 182 muPEXSISafeGuard, numElectronPEXSITolerance, & 183 matrixType, isSymbolicFactorize, ordering, rowOrdering, & 184 npSymbFact, verbosity) 185 TYPE(cp_pexsi_options), INTENT(IN) :: pexsi_options 186 REAL(KIND=real_8), INTENT(OUT), OPTIONAL :: temperature, gap, deltaE 187 INTEGER, INTENT(OUT), OPTIONAL :: numPole, isInertiaCount, & 188 maxPEXSIIter 189 REAL(KIND=real_8), INTENT(OUT), OPTIONAL :: muMin0, muMax0, mu0, & 190 muInertiaTolerance, muInertiaExpansion, muPEXSISafeGuard, & 191 numElectronPEXSITolerance 192 INTEGER, INTENT(OUT), OPTIONAL :: matrixType, & 193 isSymbolicFactorize, & 194 ordering, rowOrdering, npSymbFact, & 195 verbosity 196 197 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_pexsi_get_options', & 198 routineP = moduleN//':'//routineN 199 200#ifdef __LIBPEXSI 201 IF (PRESENT(temperature)) temperature = pexsi_options%options%temperature 202 IF (PRESENT(gap)) gap = pexsi_options%options%gap 203 IF (PRESENT(deltaE)) deltaE = pexsi_options%options%deltaE 204 IF (PRESENT(numPole)) numPole = pexsi_options%options%numPole 205 IF (PRESENT(isInertiaCount)) isInertiaCount = pexsi_options%options%isInertiaCount 206 IF (PRESENT(maxPEXSIIter)) maxPEXSIIter = pexsi_options%options%maxPEXSIIter 207 IF (PRESENT(muMin0)) muMin0 = pexsi_options%options%muMin0 208 IF (PRESENT(muMax0)) muMax0 = pexsi_options%options%muMax0 209 IF (PRESENT(mu0)) mu0 = pexsi_options%options%mu0 210 IF (PRESENT(muInertiaTolerance)) & 211 muInertiaTolerance = pexsi_options%options%muInertiaTolerance 212 IF (PRESENT(muInertiaExpansion)) & 213 muInertiaExpansion = pexsi_options%options%muInertiaExpansion 214 IF (PRESENT(muPEXSISafeGuard)) & 215 muPEXSISafeGuard = pexsi_options%options%muPEXSISafeGuard 216 IF (PRESENT(numElectronPEXSITolerance)) & 217 numElectronPEXSITolerance = pexsi_options%options%numElectronPEXSITolerance 218 IF (PRESENT(matrixType)) matrixType = pexsi_options%options%matrixType 219 IF (PRESENT(isSymbolicFactorize)) & 220 isSymbolicFactorize = pexsi_options%options%isSymbolicFactorize 221 IF (PRESENT(ordering)) ordering = pexsi_options%options%ordering 222 IF (PRESENT(rowOrdering)) rowOrdering = pexsi_options%options%rowOrdering 223 IF (PRESENT(npSymbFact)) npSymbFact = pexsi_options%options%npSymbFact 224 IF (PRESENT(verbosity)) verbosity = pexsi_options%options%verbosity 225#else 226 MARK_USED(pexsi_options) 227 ! assign intent-out arguments to silence compiler warnings 228 IF (PRESENT(temperature)) temperature = 0.0_real_8 229 IF (PRESENT(gap)) gap = 0.0_real_8 230 IF (PRESENT(deltaE)) deltaE = 0.0_real_8 231 IF (PRESENT(numPole)) numPole = -1 232 IF (PRESENT(isInertiaCount)) isInertiaCount = -1 233 IF (PRESENT(maxPEXSIIter)) maxPEXSIIter = -1 234 IF (PRESENT(muMin0)) muMin0 = 0.0_real_8 235 IF (PRESENT(muMax0)) muMax0 = 0.0_real_8 236 IF (PRESENT(mu0)) mu0 = 0.0_real_8 237 IF (PRESENT(muInertiaTolerance)) muInertiaTolerance = 0.0_real_8 238 IF (PRESENT(muInertiaExpansion)) muInertiaExpansion = 0.0_real_8 239 IF (PRESENT(muPEXSISafeGuard)) muPEXSISafeGuard = 0.0_real_8 240 IF (PRESENT(numElectronPEXSITolerance)) numElectronPEXSITolerance = 0.0_real_8 241 IF (PRESENT(matrixType)) matrixType = -1 242 IF (PRESENT(isSymbolicFactorize)) isSymbolicFactorize = -1 243 IF (PRESENT(ordering)) ordering = -1 244 IF (PRESENT(rowOrdering)) rowOrdering = -1 245 IF (PRESENT(npSymbFact)) npSymbFact = -1 246 IF (PRESENT(verbosity)) verbosity = -1 247 CPABORT("Requires linking to the PEXSI library.") 248#endif 249 END SUBROUTINE cp_pexsi_get_options 250 251! ************************************************************************************************** 252!> \brief ... 253!> \param pexsi_options ... 254! ************************************************************************************************** 255 SUBROUTINE cp_pexsi_set_default_options(pexsi_options) 256 TYPE(cp_pexsi_options), INTENT(OUT) :: pexsi_options 257 258 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_pexsi_set_default_options', & 259 routineP = moduleN//':'//routineN 260 261#ifdef __LIBPEXSI 262 CALL f_ppexsi_set_default_options(pexsi_options%options) 263#else 264 CPABORT("Requires linking to the PEXSI library.") 265#endif 266 END SUBROUTINE cp_pexsi_set_default_options 267 268! ************************************************************************************************** 269!> \brief ... 270!> \param comm ... 271!> \param numProcRow ... 272!> \param numProcCol ... 273!> \param outputFileIndex ... 274!> \return ... 275! ************************************************************************************************** 276 FUNCTION cp_pexsi_plan_initialize(comm, numProcRow, numProcCol, outputFileIndex) 277 INTEGER, INTENT(IN) :: comm, numProcRow, numProcCol, & 278 outputFileIndex 279 INTEGER(KIND=int_8) :: cp_pexsi_plan_initialize 280 281 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_pexsi_plan_initialize', & 282 routineP = moduleN//':'//routineN 283 284#ifdef __LIBPEXSI 285 INTEGER :: info, handle 286 CALL timeset(routineN, handle) 287 cp_pexsi_plan_initialize = f_ppexsi_plan_initialize(comm, numProcRow, & 288 numProcCol, outputFileIndex, info) 289 IF (info .NE. 0) & 290 CPABORT("Pexsi returned an error. Consider logPEXSI0 for details.") 291 CALL timestop(handle) 292#else 293 MARK_USED(comm) 294 MARK_USED(numProcRow) 295 MARK_USED(numProcCol) 296 MARK_USED(outputFileIndex) 297 cp_pexsi_plan_initialize = 0 298 CPABORT("Requires linking to the PEXSI library.") 299#endif 300 END FUNCTION cp_pexsi_plan_initialize 301 302! ************************************************************************************************** 303!> \brief ... 304!> \param plan ... 305!> \param pexsi_options ... 306!> \param nrows ... 307!> \param nnz ... 308!> \param nnzLocal ... 309!> \param numColLocal ... 310!> \param colptrLocal ... 311!> \param rowindLocal ... 312!> \param HnzvalLocal ... 313!> \param isSIdentity ... 314!> \param SnzvalLocal ... 315! ************************************************************************************************** 316 SUBROUTINE cp_pexsi_load_real_hs_matrix(plan, pexsi_options, nrows, nnz, & 317 nnzLocal, numColLocal, colptrLocal, & 318 rowindLocal, HnzvalLocal, isSIdentity, & 319 SnzvalLocal) 320 INTEGER(KIND=int_8), INTENT(IN) :: plan 321 TYPE(cp_pexsi_options), INTENT(IN) :: pexsi_options 322 INTEGER, INTENT(IN) :: nrows, nnz, nnzLocal, & 323 numColLocal, colptrLocal(*), & 324 rowindLocal(*) 325 REAL(KIND=real_8), INTENT(IN) :: HnzvalLocal(*) 326 INTEGER, INTENT(IN) :: isSIdentity 327 REAL(KIND=real_8), INTENT(IN) :: SnzvalLocal(*) 328 329 CHARACTER(LEN=*), PARAMETER :: & 330 routineN = 'cp_pexsi_load_real_symmetric_hs_matrix', & 331 routineP = moduleN//':'//routineN 332 333#ifdef __LIBPEXSI 334 INTEGER :: handle, info 335 336 CALL timeset(routineN, handle) 337 CALL f_ppexsi_load_real_hs_matrix(plan, pexsi_options%options, nrows, nnz, nnzLocal, & 338 numColLocal, colptrLocal, rowindLocal, & 339 HnzvalLocal, isSIdentity, SnzvalLocal, info) 340 IF (info .NE. 0) & 341 CPABORT("Pexsi returned an error. Consider logPEXSI0 for details.") 342 CALL timestop(handle) 343#else 344 MARK_USED(plan) 345 MARK_USED(pexsi_options) 346 MARK_USED(nrows) 347 MARK_USED(nnz) 348 MARK_USED(nnzLocal) 349 MARK_USED(numColLocal) 350 MARK_USED(isSIdentity) 351 CPABORT("Requires linking to the PEXSI library.") 352 353 ! MARK_USED macro does not work on assumed shape variables 354 IF (.FALSE.) THEN; DO 355 IF (colptrLocal(1) > rowindLocal(1) .OR. HnzvalLocal(1) > SnzvalLocal(1)) EXIT 356 ENDDO; ENDIF 357#endif 358 END SUBROUTINE cp_pexsi_load_real_hs_matrix 359 360! ************************************************************************************************** 361!> \brief ... 362!> \param plan ... 363!> \param pexsi_options ... 364!> \param numElectronExact ... 365!> \param muPEXSI ... 366!> \param numElectronPEXSI ... 367!> \param muMinInertia ... 368!> \param muMaxInertia ... 369!> \param numTotalInertiaIter ... 370!> \param numTotalPEXSIIter ... 371! ************************************************************************************************** 372 SUBROUTINE cp_pexsi_dft_driver(plan, pexsi_options, numElectronExact, muPEXSI, & 373 numElectronPEXSI, muMinInertia, muMaxInertia, & 374 numTotalInertiaIter, numTotalPEXSIIter) 375 INTEGER(KIND=int_8), INTENT(IN) :: plan 376 TYPE(cp_pexsi_options), INTENT(IN) :: pexsi_options 377 REAL(KIND=real_8), INTENT(IN) :: numElectronExact 378 REAL(KIND=real_8), INTENT(out) :: muPEXSI, numElectronPEXSI, & 379 muMinInertia, muMaxInertia 380 INTEGER, INTENT(out) :: numTotalInertiaIter, & 381 numTotalPEXSIIter 382 383 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_pexsi_dft_driver', & 384 routineP = moduleN//':'//routineN 385 386#ifdef __LIBPEXSI 387 INTEGER :: handle, info 388#if defined (__HAS_IEEE_EXCEPTIONS) 389 LOGICAL, DIMENSION(5) :: halt 390#endif 391 392 CALL timeset(routineN, handle) 393 394 ! Unfortuntatelly, some PEXSI kernels raise IEEE754 exceptions. 395 ! Therefore, we disable floating point traps temporarily. 396#if defined (__HAS_IEEE_EXCEPTIONS) 397 CALL ieee_get_halting_mode(IEEE_ALL, halt) 398 CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.) 399#endif 400 401 CALL f_ppexsi_dft_driver(plan, pexsi_options%options, numElectronExact, muPEXSI, & 402 numElectronPEXSI, muMinInertia, muMaxInertia, & 403 numTotalInertiaIter, numTotalPEXSIIter, info) 404 405#if defined (__HAS_IEEE_EXCEPTIONS) 406 CALL ieee_set_halting_mode(IEEE_ALL, halt) 407#endif 408 409 IF (info .NE. 0) & 410 CPABORT("Pexsi returned an error. Consider logPEXSI0 for details.") 411 CALL timestop(handle) 412#else 413 MARK_USED(plan) 414 MARK_USED(numelectronexact) 415 MARK_USED(pexsi_options) 416 ! assign intent-out arguments to silence compiler warnings 417 muPEXSI = 0.0_real_8 418 numElectronPEXSI = 0.0_real_8 419 muMinInertia = 0.0_real_8 420 muMaxInertia = 0.0_real_8 421 numTotalInertiaIter = -1 422 numTotalPEXSIIter = -1 423 CPABORT("Requires linking to the PEXSI library.") 424#endif 425 END SUBROUTINE cp_pexsi_dft_driver 426 427! ************************************************************************************************** 428!> \brief ... 429!> \param plan ... 430!> \param DMnzvalLocal ... 431!> \param EDMnzvalLocal ... 432!> \param FDMnzvalLocal ... 433!> \param totalEnergyH ... 434!> \param totalEnergyS ... 435!> \param totalFreeEnergy ... 436! ************************************************************************************************** 437 SUBROUTINE cp_pexsi_retrieve_real_dft_matrix(plan, DMnzvalLocal, EDMnzvalLocal, & 438 FDMnzvalLocal, totalEnergyH, & 439 totalEnergyS, totalFreeEnergy) 440 INTEGER(KIND=int_8), INTENT(IN) :: plan 441 REAL(KIND=real_8), INTENT(out) :: DMnzvalLocal(*), EDMnzvalLocal(*), & 442 FDMnzvalLocal(*), totalEnergyH, totalEnergyS, totalFreeEnergy 443 444 CHARACTER(LEN=*), PARAMETER :: & 445 routineN = 'cp_pexsi_retrieve_real_symmetric_dft_matrix', & 446 routineP = moduleN//':'//routineN 447 448#ifdef __LIBPEXSI 449 INTEGER :: handle, info 450 451 CALL timeset(routineN, handle) 452 CALL f_ppexsi_retrieve_real_dft_matrix(plan, DMnzvalLocal, EDMnzvalLocal, & 453 FDMnzvalLocal, totalEnergyH, & 454 totalEnergyS, totalFreeEnergy, info) 455 IF (info .NE. 0) & 456 CPABORT("Pexsi returned an error. Consider logPEXSI0 for details.") 457 CALL timestop(handle) 458#else 459 MARK_USED(plan) 460 ! assign intent-out arguments to silence compiler warnings 461 DMnzvalLocal(1) = 0.0_real_8 462 EDMnzvalLocal(1) = 0.0_real_8 463 FDMnzvalLocal(1) = 0.0_real_8 464 totalEnergyH = 0.0_real_8 465 totalEnergyS = 0.0_real_8 466 totalFreeEnergy = 0.0_real_8 467 468 CPABORT("Requires linking to the PEXSI library.") 469#endif 470 END SUBROUTINE cp_pexsi_retrieve_real_dft_matrix 471 472! ************************************************************************************************** 473!> \brief ... 474!> \param plan ... 475! ************************************************************************************************** 476 SUBROUTINE cp_pexsi_plan_finalize(plan) 477 INTEGER(KIND=int_8), INTENT(IN) :: plan 478 479 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_pexsi_plan_finalize', & 480 routineP = moduleN//':'//routineN 481 482#ifdef __LIBPEXSI 483 INTEGER :: info, handle 484 CALL timeset(routineN, handle) 485 CALL f_ppexsi_plan_finalize(plan, info) 486 IF (info .NE. 0) & 487 CPABORT("Pexsi returned an error. Consider logPEXSI0 for details.") 488 CALL timestop(handle) 489#else 490 MARK_USED(plan) 491 CPABORT("Requires linking to the PEXSI library.") 492#endif 493 END SUBROUTINE 494 495END MODULE pexsi_interface 496