1!--------------------------------------------------------------------------------------------------! 2! DFTB+: general package for performing fast atomistic simulations ! 3! Copyright (C) 2006 - 2019 DFTB+ developers group ! 4! ! 5! See the LICENSE file for terms of usage and distribution. ! 6!--------------------------------------------------------------------------------------------------! 7 8#:include 'common.fypp' 9 10!> Contains the interface to the ELSI solvers 11module dftbp_elsisolver 12#:if WITH_MPI 13 use dftbp_mpifx 14#:endif 15 use dftbp_accuracy, only : dp, lc 16 use dftbp_environment, only : TEnvironment, globalTimers 17 use dftbp_elecsolvertypes, only : electronicSolverTypes 18 use dftbp_elsiiface 19 use dftbp_elsicsc 20 use dftbp_densedescr 21 use dftbp_periodic 22 use dftbp_orbitals 23 use dftbp_message, only : error, cleanshutdown 24 use dftbp_commontypes, only : TParallelKS, TOrbitals 25 use dftbp_energies, only : TEnergies 26 use dftbp_sparse2dense 27 use dftbp_assert 28 use dftbp_spin, only : ud2qm 29 use dftbp_angmomentum, only : getLOnsite 30 use dftbp_spinorbit, only : addOnsiteSpinOrbitHam, getOnsiteSpinOrbitEnergy 31 use dftbp_potentials, only : TPotentials 32 implicit none 33 private 34 35 public :: TElsiSolverInp 36 public :: TElsiSolver, TElsiSolver_init, TElsiSolver_final 37 38 39 !> Input data for the ELSI solvers 40 type :: TElsiSolverInp 41 42 !> Choice of the solver 43 integer :: iSolver 44 45 !> Choice of ELPA solver 46 integer :: elpaSolver = 2 47 48 !> Iterations of ELPA solver before OMM minimization 49 integer :: ommIterationsElpa = 5 50 51 !> Halting tolerance for OMM iterations 52 real(dp) :: ommTolerance = 1.0E-10_dp 53 54 !> Should the overlap be factorized before minimization 55 logical :: ommCholesky = .true. 56 57 !> number of poles for PEXSI expansion 58 integer :: pexsiNPole = 20 59 60 !> number of processors per pole for PEXSI 61 integer :: pexsiNpPerPole = 1 62 63 !> number of interpolation points for mu (Fermi) search 64 integer :: pexsiNMu = 2 65 66 !> number of processors for symbolic factorisation 67 integer :: pexsiNpSymbo = 1 68 69 !> spectral radius (range of eigenvalues) if available 70 real(dp) :: pexsiDeltaE = 10.0_dp 71 72 !> density matrix purification algorithm 73 integer :: ntpolyMethod = 2 74 75 !> truncation threshold for sparse matrix multiplication 76 real(dp) :: ntpolyTruncation = 1.0E-10_dp 77 78 !> convergence tolerance for density matrix purification 79 real(dp) :: ntpolyTolerance = 1.0E-5_dp 80 81 !> Use sparse CSR format 82 logical :: elsiCsr = .false. 83 84 !> Tolerance for converting from dense matrices to internal sparse storage for libOMM, PEXSI and 85 !> NTPoly. 86 real(dp) :: elsi_zero_def 87 88 end type TElsiSolverInp 89 90 91 !> Contains settings for the solvers of the ELSI library. See ELSI manual for detailed meaning of 92 !> variables 93 type :: TElsiSolver 94 95 private 96 97 !> should the code write matrices and stop 98 logical, public :: tWriteHS 99 100 !> handle for matrix IO 101 type(elsi_rw_handle), public :: rwHandle 102 103 !> Solver id (using the central DFTB+ solver type) 104 integer :: iSolver 105 106 !> Handle for the ELSI library 107 type(elsi_handle), public :: handle 108 109 !> Use sparse CSR format 110 logical, public :: isSparse = .false. 111 112 !> Tolerance for converting from dense matrices to internal sparse storage for libOMM, PEXSI and 113 !> NTPoly. 114 real(dp) :: elsi_zero_def 115 116 !> ELSI Solver choice 117 integer :: solver 118 119 !> Level of solver information output 120 integer :: outputLevel 121 122 !> See ELSI manual - parallelisation strategy 123 integer :: parallel 124 125 !> See ELSI manual - type of parallel data decomposition 126 integer :: denseBlacs 127 128 !> Number of basis functions 129 integer :: nBasis 130 131 !> Number of electrons in the system 132 real(dp) :: nElectron 133 134 !> Maximum spin occupation for levels 135 real(dp) :: spinDegeneracy 136 137 !> Number of states to solve when relevant 138 integer :: nState 139 140 !> Global comm 141 integer :: mpiCommWorld 142 143 !> Group comm 144 integer :: myCommWorld 145 146 !> BLACS matrix context in use 147 integer :: myBlacsCtx 148 149 !> BLACS block sizes 150 integer :: BlacsBlockSize 151 152 !> CSC blocksize 153 integer :: csrBlockSize 154 155 !> Number of independent spins 156 integer :: nSpin 157 158 !> Number of independent k-points 159 integer :: nKPoint 160 161 !> Index for current spin processed here 162 integer :: iSpin 163 164 !> Index for current k-point processed here 165 integer :: iKPoint 166 167 !> Weighting for current k-point 168 real(dp) :: kWeight 169 170 !> Choice of broadening function 171 integer :: muBroadenScheme 172 173 !> If Meth-Paxton, order of scheme 174 integer :: muMpOrder 175 176 !> Whether solver had been already initialised 177 logical :: tSolverInitialised = .false. 178 179 !> Has overlap been factorized 180 logical :: tCholeskyDecomposed 181 182 !> Is this the first calls for this geometry 183 logical :: tFirstCalc = .true. 184 185 !> Sparse format data structure 186 type(TElsiCsc), allocatable :: elsiCsc 187 188 189 !! ELPA settings 190 191 !> ELPA solver choice 192 integer :: elpaSolverOption 193 194 !! OMM settings 195 196 !> Starting iterations with ELPA 197 integer :: ommIter 198 199 !> Tolerance for minimisation 200 real(dp) :: ommTolerance 201 202 !> Whether to Cholesky factorize and transform or work with general 203 logical :: ommCholesky 204 205 !! PEXSI settings 206 207 !> Minimum contour range 208 real(dp) :: pexsiMuMin 209 210 !> Maximum contour range 211 real(dp) :: pexsiMuMax 212 213 !> Most negative potential 214 real(dp) :: pexsiDeltaVMin 215 216 !> Most positive potential 217 real(dp) :: pexsiDeltaVMax 218 219 !> Previous potentials 220 real(dp), allocatable :: pexsiVOld(:) 221 222 !> Number of poles for expansion 223 integer :: pexsiNPole 224 225 !> Processors per pole 226 integer :: pexsiNpPerPole 227 228 !> Processes for chemical potential search 229 integer :: pexsiNMu 230 231 !> Processors used for symbolic factorization stage 232 integer :: pexsiNpSymbo 233 234 !> Spectral radius 235 real(dp) :: pexsiDeltaE 236 237 !! NTPoly settings 238 239 !> Choice of minimisation method 240 integer :: ntpolyMethod 241 242 !> Truncation threshold for elements 243 real(dp) :: ntpolyTruncation 244 245 !> Convergence tolerance 246 real(dp) :: ntpolyTolerance 247 248 contains 249 250 procedure :: reset => TElsiSolver_reset 251 procedure :: updateGeometry => TElsiSolver_updateGeometry 252 procedure :: updateElectronicTemp => TElsiSolver_updateElectronicTemp 253 procedure :: getDensity => TElsiSolver_getDensity 254 procedure :: getEDensity => TElsiSolver_getEDensity 255 procedure :: initPexsiDeltaVRanges => TElsiSolver_initPexsiDeltaVRanges 256 procedure :: updatePexsiDeltaVRanges => TElsiSolver_updatePexsiDeltaVRanges 257 procedure :: getSolverName => TElsiSolver_getSolverName 258 259 end type TElsiSolver 260 261 262contains 263 264 265 !> Initialise ELSI solver 266 subroutine TElsiSolver_init(this, inp, env, nBasisFn, nEl, iDistribFn, nSpin, iSpin, nKPoint,& 267 & iKPoint, kWeight, tWriteHS) 268 269 !> control structure for solvers, including ELSI data 270 class(TElsiSolver), intent(out) :: this 271 272 !> input structure for ELSI 273 type(TElsiSolverInp), intent(in) :: inp 274 275 !> Environment settings 276 type(TEnvironment), intent(in) :: env 277 278 !> number of orbitals in the system 279 integer, intent(in) :: nBasisFn 280 281 !> number of electrons 282 real(dp), intent(in) :: nEl(:) 283 284 !> filling function 285 integer, intent(in) :: iDistribFn 286 287 !> total number of spin channels. 288 integer, intent(in) :: nSpin 289 290 !> spin channel processed by current process 291 integer, intent(in) :: iSpin 292 293 !> total number of k-points 294 integer, intent(in) :: nKPoint 295 296 !> K-point processed by current process. 297 integer, intent(in) :: iKPoint 298 299 !> Weight of current k-point 300 real(dp), intent(in) :: kWeight 301 302 !> Should the matrices be written out 303 logical, intent(in) :: tWriteHS 304 305 #:if WITH_ELSI 306 307 this%iSolver = inp%iSolver 308 309 select case(this%iSolver) 310 311 case (electronicSolverTypes%elpa) 312 this%solver = 1 313 ! ELPA is asked for all states 314 this%nState = nBasisFn 315 316 case (electronicSolverTypes%omm) 317 this%solver = 2 318 ! OMM solves only over occupied space 319 ! spin degeneracies for closed shell 320 if (nSpin == 1) then 321 this%nState = nint(sum(nEl)*0.5_dp) 322 else 323 this%nState = nint(sum(nEl)) 324 end if 325 326 case (electronicSolverTypes%pexsi) 327 this%solver = 3 328 ! ignored by PEXSI: 329 this%nState = nBasisFn 330 331 case (electronicSolverTypes%ntpoly) 332 this%solver = 6 333 ! ignored by NTPoly, but set anyway: 334 this%nState = nBasisFn 335 336 end select 337 338 ! parallelism with multiple processes 339 this%parallel = 1 340 341 ! number of basis functions 342 this%nBasis = nBasisFn 343 344 ! number of electrons in the problem 345 this%nElectron = sum(nEl) 346 347 ! should the code write the hamiltonian and overlap then stop 348 this%tWriteHS = tWriteHS 349 350 if (nSpin == 2) then 351 this%spinDegeneracy = nEl(1) - nEl(2) 352 else 353 this%spinDegeneracy = 0.0_dp 354 end if 355 356 ! Number of spin channels passed to the ELSI library 357 this%nSpin = min(nSpin, 2) 358 this%iSpin = iSpin 359 this%nKPoint = nKPoint 360 this%iKPoint = iKPoint 361 this%kWeight = kWeight 362 363 this%mpiCommWorld = env%mpi%globalComm%id 364 this%myCommWorld = env%mpi%groupComm%id 365 this%myBlacsCtx = env%blacs%orbitalGrid%ctxt 366 367 ! assumes row and column sizes the same 368 this%BlacsBlockSize = env%blacs%rowBlockSize 369 370 this%csrBlockSize = this%nBasis / env%mpi%groupComm%size 371 372 if (this%csrBlockSize * env%mpi%groupComm%size < this%nBasis) then 373 this%csrBlockSize = this%csrBlockSize + 1 374 end if 375 376 this%muBroadenScheme = min(iDistribFn,2) 377 if (iDistribFn > 1) then 378 ! set Meth-Pax order 379 this%muMpOrder = iDistribFn - 2 380 else 381 this%muMpOrder = 0 382 end if 383 384 ! ELPA settings 385 this%elpaSolverOption = inp%elpaSolver 386 387 ! OMM settings 388 this%ommIter = inp%ommIterationsElpa 389 this%ommTolerance = inp%ommTolerance 390 this%ommCholesky = inp%ommCholesky 391 392 ! PEXSI settings 393 this%pexsiNPole = inp%pexsiNPole 394 this%pexsiNpPerPole = inp%pexsiNpPerPole 395 this%pexsiNMu = inp%pexsiNMu 396 this%pexsiNpSymbo = inp%pexsiNpSymbo 397 this%pexsiDeltaE = inp%pexsiDeltaE 398 399 ! NTPoly settings 400 this%ntpolyMethod = inp%ntpolyMethod 401 this%ntpolyTruncation = inp%ntpolyTruncation 402 this%ntpolyTolerance = inp%ntpolyTolerance 403 404 this%isSparse = inp%elsiCsr 405 406 this%elsi_zero_def = inp%elsi_zero_def 407 408 ! data as dense BLACS blocks 409 if (this%isSparse) then 410 ! CSR format 411 this%denseBlacs = 2 412 else 413 this%denseBlacs = 0 414 end if 415 416 if (this%isSparse) then 417 allocate(this%elsiCsc) 418 end if 419 420 ! customize output level, note there are levels 0..3 not DFTB+ 0..2 421 this%OutputLevel = 0 422 #:call DEBUG_CODE 423 this%OutputLevel = 3 424 #:endcall DEBUG_CODE 425 426 if (nSpin == 4 .and. (this%isSparse .and. this%solver /= 1)) then 427 call error("Current solver configuration not avaible for two component complex& 428 & hamiltonians") 429 end if 430 431 this%tCholeskyDecomposed = .false. 432 433 #:else 434 435 call error("Internal error: TElsiSolver_init() called despite missing ELSI support") 436 437 #:endif 438 439 end subroutine TElsiSolver_init 440 441 442 !> Finalizes the ELSI solver. 443 subroutine TELsiSolver_final(this) 444 445 !> Instance 446 type(TElsiSolver), intent(inout) :: this 447 448 #:if WITH_ELSI 449 450 call elsi_finalize(this%handle) 451 452 #:else 453 454 call error("Internal error: TELsiSolver_final() called despite missing ELSI& 455 & support") 456 457 #:endif 458 459 end subroutine TELsiSolver_final 460 461 462 !> reset the ELSI solver - safer to do this on geometry change, due to the lack of a Choleskii 463 !> refactorization option 464 subroutine TElsiSolver_reset(this) 465 466 !> Instance 467 class(TElsiSolver), intent(inout) :: this 468 469 470 #:if WITH_ELSI 471 472 if (this%tSolverInitialised) then 473 474 ! reset previous instance of solver 475 call elsi_reinit(this%handle) 476 477 if (this%iSolver == electronicSolverTypes%pexsi) then 478 ! reset PEXSI chemical potential search 479 this%pexsiMuMin = -10.0_dp 480 this%pexsiMuMax = 10.0_dp 481 this%pexsiDeltaVMin = 0.0_dp 482 this%pexsiDeltaVMax = 0.0_dp 483 end if 484 485 else 486 487 ! initialise solver 488 489 call elsi_init(this%handle, this%solver, this%parallel, this%denseBlacs, this%nBasis,& 490 & this%nElectron, this%nState) 491 492 call elsi_set_mpi_global(this%handle, this%mpiCommWorld) 493 call elsi_set_sing_check(this%handle, 0) ! disable singularity check 494 call elsi_set_mpi(this%handle, this%myCommWorld) 495 496 if (this%isSparse) then 497 call elsi_set_csc_blk(this%handle, this%csrBlockSize) 498 else 499 call elsi_set_zero_def(this%handle, this%elsi_zero_def) 500 end if 501 call elsi_set_blacs(this%handle, this%myBlacsCtx, this%BlacsBlockSize) 502 503 if (this%tWriteHS) then 504 ! setup to write a matrix 505 call elsi_init_rw(this%rwHandle, 1, this%parallel, this%nBasis, this%nElectron) 506 ! MPI comm 507 call elsi_set_rw_mpi(this%rwHandle, this%mpiCommWorld) 508 if (.not.this%isSparse) then 509 ! dense matrices 510 call elsi_set_rw_blacs(this%rwHandle, this%myBlacsCtx, this%BlacsBlockSize) 511 ! use default of 1E-15 for the moment, but could be over-ridden if needed, probably with a 512 ! different variable name: 513 !call elsi_set_rw_zero_def(this%rwHandle, this%elsi_zero_def) 514 end if 515 end if 516 517 518 select case(this%iSolver) 519 case(electronicSolverTypes%elpa) 520 521 select case(this%elpaSolverOption) 522 case(1) 523 ! single stage 524 call elsi_set_elpa_solver(this%handle, 1) 525 case(2) 526 ! two stage 527 call elsi_set_elpa_solver(this%handle, 2) 528 case default 529 call error("Unknown ELPA solver modes") 530 end select 531 532 case(electronicSolverTypes%omm) 533 ! libOMM 534 if (this%ommCholesky) then 535 call elsi_set_omm_flavor(this%handle, 2) 536 else 537 call elsi_set_omm_flavor(this%handle, 0) 538 end if 539 call elsi_set_omm_n_elpa(this%handle, this%ommIter) 540 call elsi_set_omm_tol(this%handle, this%ommTolerance) 541 542 case(electronicSolverTypes%pexsi) 543 544 this%pexsiMuMin = -10.0_dp 545 this%pexsiMuMax = 10.0_dp 546 this%pexsiDeltaVMin = 0.0_dp 547 this%pexsiDeltaVMax = 0.0_dp 548 549 ! processors per pole to invert for 550 call elsi_set_pexsi_np_per_pole(this%handle, this%pexsiNpPerPole) 551 552 call elsi_set_pexsi_mu_min(this%handle, this%pexsiMuMin +this%pexsiDeltaVMin) 553 call elsi_set_pexsi_mu_max(this%handle, this%pexsiMuMax +this%pexsiDeltaVMax) 554 555 ! number of poles for the expansion 556 call elsi_set_pexsi_n_pole(this%handle, this%pexsiNPole) 557 558 ! number of interpolation points for mu 559 call elsi_set_pexsi_n_mu(this%handle, this%pexsiNMu) 560 561 ! number of processors for symbolic factorisation task 562 call elsi_set_pexsi_np_symbo(this%handle, this%pexsiNpSymbo) 563 564 ! spectral radius (range of eigenspectrum, if known, otherwise default usually fine) 565 call elsi_set_pexsi_delta_e(this%handle, this%pexsiDeltaE) 566 567 case(electronicSolverTypes%ntpoly) 568 569 ! NTPoly 570 ! set purification method 571 call elsi_set_ntpoly_method(this%handle, this%ntpolyMethod) 572 573 ! set truncation tolerance for sparse matrix multiplications 574 call elsi_set_ntpoly_filter(this%handle, this%ntpolyTruncation) 575 576 ! set purification convergence threshold 577 call elsi_set_ntpoly_tol(this%handle, this%ntpolyTolerance) 578 579 end select 580 581 if (any(this%iSolver == [electronicSolverTypes%omm,& 582 & electronicSolverTypes%pexsi, electronicSolverTypes%ntpoly])) then 583 ! density matrix build needs to know the number of spin channels to normalize against 584 select case(this%nSpin) 585 case(1) 586 call elsi_set_spin(this%handle, 1, this%iSpin) 587 case(2) 588 call elsi_set_spin(this%handle, 2, this%iSpin) 589 case(4) 590 call elsi_set_spin(this%handle, 2, this%iSpin) 591 end select 592 if (this%nKPoint > 1) then 593 call elsi_set_kpoint(this%handle, this%nKPoint, this%iKPoint, this%kWeight) 594 end if 595 end if 596 597 call elsi_set_output(this%handle, this%OutputLevel) 598 if (this%OutputLevel == 3) then 599 call elsi_set_output_log(this%handle, 1) 600 end if 601 this%tSolverInitialised = .true. 602 end if 603 604 this%tCholeskyDecomposed = .false. 605 this%tFirstCalc = .true. 606 607 #:else 608 609 call error("Internal error: TElsiSolver_reset() called despite missing ELSI support") 610 611 #:endif 612 613 end subroutine TElsiSolver_reset 614 615 616 !> Resets solver due to geometry changes 617 subroutine TElsiSolver_updateGeometry(this, env, neighList, nNeighbourSK, iAtomStart,& 618 & iSparseStart, img2CentCell) 619 620 !> Instance 621 class(TElsiSolver), intent(inout) :: this 622 623 !> Environment settings 624 type(TEnvironment), intent(inout) :: env 625 626 !> Neighbour list 627 type(TNeighbourList), intent(in) :: neighList 628 629 !> Number of neighbours for each of the atoms 630 integer, intent(in) :: nNeighbourSK(:) 631 632 !> Atom offset for the squared matrix 633 integer, intent(in) :: iAtomStart(:) 634 635 !> indexing array for the sparse Hamiltonian 636 integer, intent(in) :: iSparseStart(0:,:) 637 638 !> Mapping between image atoms and corresponding atom in the central cell. 639 integer, intent(in) :: img2CentCell(:) 640 641 #:if WITH_ELSI 642 643 call TElsiCsc_init(this%elsiCsc, env, this%nBasis, this%csrBlockSize, neighList,& 644 & nNeighbourSK, iAtomStart, iSparseStart, img2CentCell) 645 646 #:else 647 648 call error("Internal error: TELsiSolver_updateGeometry() called despite missing ELSI& 649 & support") 650 651 #:endif 652 653 end subroutine TElsiSolver_updateGeometry 654 655 656 !> Updated the electronic temperature 657 subroutine TElsiSolver_updateElectronicTemp(this, tempElec) 658 659 !> Instance 660 class(TElsiSolver), intent(inout) :: this 661 662 !> Electronic temperature 663 real(dp), intent(in) :: tempElec 664 665 #:if WITH_ELSI 666 667 call elsi_set_mu_broaden_scheme(this%handle, this%muBroadenScheme) 668 if (this%muBroadenScheme == 2) then 669 call elsi_set_mu_mp_order(this%handle, this%muMpOrder) 670 end if 671 call elsi_set_mu_broaden_width(this%handle, tempElec) 672 if (this%iSolver == electronicSolverTypes%pexsi) then 673 call elsi_set_pexsi_temp(this%handle, tempElec) 674 end if 675 676 #:else 677 678 call error("Internal error: TELsiSolver_updateElectronicTemp() called despite missing ELSI& 679 & support") 680 681 #:endif 682 683 end subroutine TElsiSolver_updateElectronicTemp 684 685 686 !> Returns the density matrix using ELSI non-diagonalisation routines. 687 subroutine TElsiSolver_getDensity(this, env, denseDesc, ham, over, neighbourList, nNeighbourSK,& 688 & iSparseStart, img2CentCell, iCellVec, cellVec, kPoint, kWeight, orb, species, tRealHS,& 689 & tSpinSharedEf, tSpinOrbit, tDualSpinOrbit, tMulliken, parallelKS, Ef, energy, rhoPrim,& 690 & Eband, TS, iHam, xi, orbitalL, HSqrReal, SSqrReal, iRhoPrim, HSqrCplx, SSqrCplx) 691 692 !> Electronic solver information 693 class(TElsiSolver), intent(inout) :: this 694 695 !> Environment settings 696 type(TEnvironment), intent(inout) :: env 697 698 !> Dense matrix descriptor 699 type(TDenseDescr), intent(in) :: denseDesc 700 701 !> hamiltonian in sparse storage 702 real(dp), intent(in) :: ham(:,:) 703 704 !> sparse overlap matrix 705 real(dp), intent(in) :: over(:) 706 707 !> list of neighbours for each atom 708 type(TNeighbourList), intent(in) :: neighbourList 709 710 !> Number of neighbours for each of the atoms 711 integer, intent(in) :: nNeighbourSK(:) 712 713 !> Index array for the start of atomic blocks in sparse arrays 714 integer, intent(in) :: iSparseStart(:,:) 715 716 !> map from image atoms to the original unique atom 717 integer, intent(in) :: img2CentCell(:) 718 719 !> Index for which unit cell atoms are associated with 720 integer, intent(in) :: iCellVec(:) 721 722 !> Vectors (in units of the lattice constants) to cells of the lattice 723 real(dp), intent(in) :: cellVec(:,:) 724 725 !> k-points 726 real(dp), intent(in) :: kPoint(:,:) 727 728 !> Weights for k-points 729 real(dp), intent(in) :: kWeight(:) 730 731 !> Atomic orbital information 732 type(TOrbitals), intent(in) :: orb 733 734 !> species of all atoms in the system 735 integer, intent(in) :: species(:) 736 737 !> Is the hamitonian real (no k-points/molecule/gamma point)? 738 logical, intent(in) :: tRealHS 739 740 !> Is the Fermi level common accross spin channels? 741 logical, intent(in) :: tSpinSharedEf 742 743 !> Are spin orbit interactions present 744 logical, intent(in) :: tSpinOrbit 745 746 !> Are block population spin orbit interactions present 747 logical, intent(in) :: tDualSpinOrbit 748 749 !> Should Mulliken populations be generated/output 750 logical, intent(in) :: tMulliken 751 752 !> K-points and spins to process 753 type(TParallelKS), intent(in) :: parallelKS 754 755 !> Fermi level(s) 756 real(dp), intent(inout) :: Ef(:) 757 758 !> Energy contributions and total 759 type(TEnergies), intent(inout) :: energy 760 761 !> sparse density matrix 762 real(dp), intent(out) :: rhoPrim(:,:) 763 764 !> band structure energy 765 real(dp), intent(out) :: Eband(:) 766 767 !> electronic entropy times temperature 768 real(dp), intent(out) :: TS(:) 769 770 !> imaginary part of hamitonian 771 real(dp), intent(in), allocatable :: iHam(:,:) 772 773 !> spin orbit constants 774 real(dp), intent(in), allocatable :: xi(:,:) 775 776 !> orbital moments of atomic shells 777 real(dp), intent(inout), allocatable :: orbitalL(:,:,:) 778 779 !> imaginary part of density matrix 780 real(dp), intent(inout), allocatable :: iRhoPrim(:,:) 781 782 !> dense real hamiltonian storage 783 real(dp), intent(inout), allocatable :: HSqrReal(:,:) 784 785 !> dense real overlap storage 786 real(dp), intent(inout), allocatable :: SSqrReal(:,:) 787 788 !> dense complex (k-points) hamiltonian storage 789 complex(dp), intent(inout), allocatable :: HSqrCplx(:,:) 790 791 !> dense complex (k-points) overlap storage 792 complex(dp), intent(inout), allocatable :: SSqrCplx(:,:) 793 794 #:if WITH_ELSI 795 796 integer :: nSpin 797 integer :: iKS, iK, iS 798 799 ! as each spin and k-point combination forms a separate group for this solver, then iKS = 1 800 iKS = 1 801 iK = parallelKS%localKS(1, iKS) 802 iS = parallelKS%localKS(2, iKS) 803 804 nSpin = size(ham, dim=2) 805 if (nSpin == 2 .and. .not. tSpinSharedEf) then 806 call error("ELSI currently requires shared Fermi levels over spin") 807 end if 808 809 rhoPrim(:,:) = 0.0_dp 810 if (allocated(iRhoPrim)) then 811 iRhoPrim(:,:) = 0.0_dp 812 end if 813 814 if (this%iSolver == electronicSolverTypes%pexsi) then 815 call elsi_set_pexsi_mu_min(this%handle, this%pexsiMuMin + this%pexsiDeltaVMin) 816 call elsi_set_pexsi_mu_max(this%handle, this%pexsiMuMax + this%pexsiDeltaVMax) 817 end if 818 819 if (nSpin /= 4) then 820 if (tRealHS) then 821 if (this%isSparse) then 822 call getDensityRealSparse(this, parallelKS, ham, over, neighbourList%iNeighbour,& 823 & nNeighbourSK, denseDesc%iAtomStart, iSparseStart, img2CentCell, orb, rhoPrim, Eband) 824 else 825 call getDensityRealDense(this, env, denseDesc, ham, over, neighbourList,& 826 & nNeighbourSK, iSparseStart, img2CentCell, orb, parallelKS, rhoPrim, Eband,& 827 & HSqrReal, SSqrReal) 828 end if 829 else 830 if (this%isSparse) then 831 call getDensityCmplxSparse(this, parallelKS, kPoint(:,iK), kWeight(iK), iCellVec,& 832 & cellVec, ham, over, neighbourList%iNeighbour, nNeighbourSK, denseDesc%iAtomStart,& 833 & iSparseStart, img2CentCell, orb, rhoPrim, Eband) 834 else 835 call getDensityCmplxDense(this, env, denseDesc, ham, over, neighbourList,& 836 & nNeighbourSK, iSparseStart, img2CentCell, iCellVec, cellVec, kPoint, kWeight, orb,& 837 & parallelKS, rhoPrim, Eband, HSqrCplx, SSqrCplx) 838 end if 839 end if 840 call ud2qm(rhoPrim) 841 else 842 if (this%isSparse) then 843 call error("Internal error: getDensityFromElsiSolver : sparse pauli not yet implemented") 844 else 845 call getDensityPauliDense(this, env, denseDesc, ham, over, neighbourList,& 846 & nNeighbourSK, iSparseStart, img2CentCell, iCellVec, cellVec, kPoint, kWeight, orb,& 847 & species, tSpinOrbit, tDualSpinOrbit, tMulliken, parallelKS, energy, rhoPrim, Eband,& 848 & iHam, xi, orbitalL, iRhoPrim, HSqrCplx, SSqrCplx) 849 end if 850 end if 851 852 Ef(:) = 0.0_dp 853 TS(:) = 0.0_dp 854 if (env%mpi%tGroupMaster) then 855 call elsi_get_mu(this%handle, Ef(iS)) 856 call elsi_get_entropy(this%handle, TS(iS)) 857 end if 858 call mpifx_allreduceip(env%mpi%globalComm, Ef, MPI_SUM) 859 call mpifx_allreduceip(env%mpi%globalComm, TS, MPI_SUM) 860 861 if (this%iSolver == electronicSolverTypes%pexsi) then 862 call elsi_get_pexsi_mu_min(this%handle, this%pexsiMuMin) 863 call elsi_get_pexsi_mu_max(this%handle, this%pexsiMuMax) 864 end if 865 866 ! Add up and distribute density matrix contribution from each group 867 call mpifx_allreduceip(env%mpi%globalComm, rhoPrim, MPI_SUM) 868 if (allocated(iRhoPrim)) then 869 call mpifx_allreduceip(env%mpi%globalComm, iRhoPrim, MPI_SUM) 870 end if 871 872 if (tSpinOrbit) then 873 call mpifx_allreduceip(env%mpi%globalComm, energy%atomLS, MPI_SUM) 874 energy%atomLS(:) = 2.0_dp * energy%atomLS 875 end if 876 if (tMulliken .and. tSpinOrbit .and. .not. tDualSpinOrbit) then 877 call mpifx_allreduceip(env%mpi%globalComm, orbitalL, MPI_SUM) 878 orbitalL(:,:,:) = 2.0_dp * orbitalL 879 end if 880 if (tSpinOrbit .and. .not. tDualSpinOrbit) then 881 energy%ELS = sum(energy%atomLS) 882 end if 883 884 call env%globalTimer%stopTimer(globalTimers%densityMatrix) 885 886 #:else 887 888 call error("Internal error: TELsiSolver_getDensity() called despite missing ELSI support") 889 890 #:endif 891 892 end subroutine TElsiSolver_getDensity 893 894 895 ! Returns the energy weighted density matrix using ELSI non-diagonalisation routines. 896 subroutine TElsiSolver_getEDensity(this, env, denseDesc, nSpin, kPoint, kWeight, neighbourList,& 897 & nNeighbourSK, orb, iSparseStart, img2CentCell, iCellVec, cellVec, tRealHS, parallelKS,& 898 & ERhoPrim, SSqrReal, SSqrCplx) 899 900 !> Electronic solver information 901 class(TElsiSolver), intent(inout) :: this 902 903 !> Environment settings 904 type(TEnvironment), intent(inout) :: env 905 906 !> Dense matrix descriptor 907 type(TDenseDescr), intent(in) :: denseDesc 908 909 !> Number of spin channels 910 integer, intent(in) :: nSpin 911 912 !> K-points 913 real(dp), intent(in) :: kPoint(:,:) 914 915 !> Weights for k-points 916 real(dp), intent(in) :: kWeight(:) 917 918 !> list of neighbours for each atom 919 type(TNeighbourList), intent(in) :: neighbourList 920 921 !> Number of neighbours for each of the atoms 922 integer, intent(in) :: nNeighbourSK(:) 923 924 !> Atomic orbital information 925 type(TOrbitals), intent(in) :: orb 926 927 !> Index array for the start of atomic blocks in sparse arrays 928 integer, intent(in) :: iSparseStart(:,:) 929 930 !> map from image atoms to the original unique atom 931 integer, intent(in) :: img2CentCell(:) 932 933 !> Index for which unit cell atoms are associated with 934 integer, intent(in) :: iCellVec(:) 935 936 !> Vectors (in units of the lattice constants) to cells of the lattice 937 real(dp), intent(in) :: cellVec(:,:) 938 939 !> Is the hamitonian real (no k-points/molecule/gamma point)? 940 logical, intent(in) :: tRealHS 941 942 !> K-points and spins to process 943 type(TParallelKS), intent(in) :: parallelKS 944 945 !> Energy weighted sparse matrix 946 real(dp), intent(out) :: ERhoPrim(:) 947 948 !> Storage for dense overlap matrix 949 real(dp), intent(inout), allocatable :: SSqrReal(:,:) 950 951 !> Storage for dense overlap matrix (complex case) 952 complex(dp), intent(inout), allocatable :: SSqrCplx(:,:) 953 954 #:if WITH_ELSI 955 956 if (nSpin == 4) then 957 call getEDensityMtxPauli(this, env, denseDesc, kPoint, kWeight,& 958 & neighbourList, nNeighbourSK, orb, iSparseStart, img2CentCell, iCellVec, cellVec,& 959 & parallelKS, ERhoPrim, SSqrCplx) 960 else 961 if (tRealHS) then 962 call getEDensityMtxReal(this, env, denseDesc, neighbourList, nNeighbourSK, orb,& 963 & iSparseStart, img2CentCell, ERhoPrim, SSqrReal) 964 else 965 call getEDensityMtxCmplx(this, env, denseDesc, kPoint, kWeight, neighbourList,& 966 & nNeighbourSK, orb, iSparseStart, img2CentCell, iCellVec, cellVec, parallelKS,& 967 & ERhoPrim, SSqrCplx) 968 end if 969 end if 970 971 #:else 972 973 call error("Internal error: TELsiSolver_getDensity() called despite missing ELSI support") 974 975 #:endif 976 977 end subroutine TElsiSolver_getEDensity 978 979 980 !> Initializes the PEXSI potentials. 981 subroutine TElsiSolver_initPexsiDeltaVRanges(this, tSccCalc, potential) 982 983 !> Electronic solver information 984 class(TElsiSolver), intent(inout) :: this 985 986 !> Whether we have an SCC calculation 987 logical, intent(in) :: tSccCalc 988 989 !> Potentials acting 990 type(TPotentials), intent(in) :: potential 991 992 #:if WITH_ELSI 993 994 if (tSccCalc) then 995 if (.not.allocated(this%pexsiVOld)) then 996 allocate(this%pexsiVOld(size(potential%intBlock(:,:,:,1)))) 997 end if 998 this%pexsiVOld = reshape(potential%intBlock(:,:,:,1) + potential%extBlock(:,:,:,1),& 999 & [size(potential%extBlock(:,:,:,1))]) 1000 end if 1001 this%pexsiDeltaVMin = 0.0_dp 1002 this%pexsiDeltaVMax = 0.0_dp 1003 1004 #:else 1005 1006 call error("Internal error: TELsiSolver_initPexsiDeltaVRanges() called despite missing ELSI& 1007 & support") 1008 1009 #:endif 1010 1011 end subroutine TElsiSolver_initPexsiDeltaVRanges 1012 1013 1014 !> Updates the PEXSI potentials. 1015 subroutine TElsiSolver_updatePexsiDeltaVRanges(this, potential) 1016 1017 !> Electronic solver information 1018 class(TElsiSolver), intent(inout) :: this 1019 1020 !> Potentials acting 1021 type(TPotentials), intent(in) :: potential 1022 1023 #:if WITH_ELSI 1024 1025 this%pexsiDeltaVMin = minval(this%pexsiVOld& 1026 & - reshape(potential%intBlock(:,:,:,1) + potential%extBlock(:,:,:,1),& 1027 & [size(potential%extBlock(:,:,:,1))])) 1028 this%pexsiDeltaVMax = maxval(this%pexsiVOld& 1029 & - reshape(potential%intBlock(:,:,:,1) + potential%extBlock(:,:,:,1),& 1030 & [size(potential%extBlock(:,:,:,1))])) 1031 this%pexsiVOld = reshape(potential%intBlock(:,:,:,1) + potential%extBlock(:,:,:,1),& 1032 & [size(potential%extBlock(:,:,:,1))]) 1033 1034 #:else 1035 1036 call error("Internal error: TElsiSolver_updatePexsiDeltaVRanges() called despite missing ELSI& 1037 & support") 1038 1039 #:endif 1040 1041 end subroutine TElsiSolver_updatePexsiDeltaVRanges 1042 1043 1044 !> Returns the name of the solver 1045 function TElsiSolver_getSolverName(this) result(solverName) 1046 1047 !> Instance. 1048 class(TElsiSolver), intent(in) :: this 1049 1050 !> Name of the solver. 1051 character(:), allocatable :: solverName 1052 1053 #:if WITH_ELSI 1054 1055 character(lc) :: buffer 1056 1057 select case (this%iSolver) 1058 1059 case(electronicSolverTypes%elpa) 1060 if (this%elpaSolverOption == 1) then 1061 write(buffer, "(A)") "ELSI interface to the 1 stage ELPA solver" 1062 else 1063 write(buffer, "(A)") "ELSI interface to the 2 stage ELPA solver" 1064 end if 1065 1066 case(electronicSolverTypes%omm) 1067 write(buffer, "(A,I0,A,E9.2)") "ELSI solver libOMM with ",& 1068 & this%ommIter, " ELPA iterations",this%ommTolerance 1069 if (this%isSparse) then 1070 write(buffer, "(A)") "ELSI solver libOMM Sparse" 1071 else 1072 write(buffer, "(A)") "ELSI solver libOMM Dense" 1073 end if 1074 1075 case(electronicSolverTypes%pexsi) 1076 if (this%isSparse) then 1077 write(buffer, "(A)") "ELSI solver PEXSI Sparse" 1078 else 1079 write(buffer, "(A)") "ELSI solver PEXSI Dense" 1080 end if 1081 1082 case(electronicSolverTypes%ntpoly) 1083 if (this%isSparse) then 1084 write(buffer, "(A)") "ELSI solver NTPoly Sparse" 1085 else 1086 write(buffer, "(A)") "ELSI solver NTPoly Dense" 1087 end if 1088 1089 case default 1090 write(buffer, "(A,I0,A)") "Invalid electronic solver! (iSolver = ", this%iSolver, ")" 1091 1092 end select 1093 1094 solverName = trim(buffer) 1095 1096 #:else 1097 1098 solverName = "" 1099 1100 call error("Internal error: TElsiSolver_getSolverName() called despite missing ELSI support") 1101 1102 #:endif 1103 1104 end function TElsiSolver_getSolverName 1105 1106 1107!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1108!!! Private routines 1109!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1110 1111#:if WITH_ELSI 1112 1113 !> Returns the density matrix using ELSI non-diagonalisation routines (real dense case). 1114 subroutine getDensityRealDense(this, env, denseDesc, ham, over, neighbourList,& 1115 & nNeighbourSK, iSparseStart, img2CentCell, orb, parallelKS, rhoPrim, Eband, HSqrReal,& 1116 & SSqrReal) 1117 1118 !> Electronic solver information 1119 type(TElsiSolver), intent(inout) :: this 1120 1121 !> Environment settings 1122 type(TEnvironment), intent(inout) :: env 1123 1124 !> Dense matrix descriptor 1125 type(TDenseDescr), intent(in) :: denseDesc 1126 1127 !> hamiltonian in sparse storage 1128 real(dp), intent(in) :: ham(:,:) 1129 1130 !> sparse overlap matrix 1131 real(dp), intent(in) :: over(:) 1132 1133 !> list of neighbours for each atom 1134 type(TNeighbourList), intent(in) :: neighbourList 1135 1136 !> Number of neighbours for each of the atoms 1137 integer, intent(in) :: nNeighbourSK(:) 1138 1139 !> Index array for the start of atomic blocks in sparse arrays 1140 integer, intent(in) :: iSparseStart(:,:) 1141 1142 !> map from image atoms to the original unique atom 1143 integer, intent(in) :: img2CentCell(:) 1144 1145 !> Atomic orbital information 1146 type(TOrbitals), intent(in) :: orb 1147 1148 !> K-points and spins to process 1149 type(TParallelKS), intent(in) :: parallelKS 1150 1151 !> sparse density matrix 1152 real(dp), intent(inout) :: rhoPrim(:,:) 1153 1154 !> band structure energy 1155 real(dp), intent(out) :: Eband(:) 1156 1157 !> dense real hamiltonian storage 1158 real(dp), intent(inout), allocatable :: HSqrReal(:,:) 1159 1160 !> dense real overlap storage 1161 real(dp), intent(inout), allocatable :: SSqrReal(:,:) 1162 1163 real(dp), allocatable :: rhoSqrReal(:,:) 1164 integer :: iKS, iS 1165 1166 ! as each spin and k-point combination forms a separate group for this solver, then iKS = 1 1167 iKS = 1 1168 iS = parallelKS%localKS(2, iKS) 1169 1170 call unpackHSRealBlacs(env%blacs, ham(:,iS), neighbourList%iNeighbour, nNeighbourSK,& 1171 & iSparseStart, img2CentCell, denseDesc, HSqrReal) 1172 if (.not. this%tCholeskyDecomposed) then 1173 call unpackHSRealBlacs(env%blacs, over, neighbourList%iNeighbour, nNeighbourSK,& 1174 & iSparseStart, img2CentCell, denseDesc, SSqrReal) 1175 if (this%ommCholesky .or. this%iSolver == electronicSolverTypes%pexsi) then 1176 this%tCholeskyDecomposed = .true. 1177 end if 1178 end if 1179 1180 allocate(rhoSqrReal(size(HSqrReal, dim=1), size(HSqrReal, dim=2))) 1181 rhoSqrReal(:,:) = 0.0_dp 1182 Eband(iS) = 0.0_dp 1183 if (this%tWriteHS) then 1184 call elsi_write_mat_real(this%rwHandle, "ELSI_Hreal.bin", HSqrReal) 1185 call elsi_write_mat_real(this%rwHandle, "ELSI_Sreal.bin", SSqrReal) 1186 call elsi_finalize_rw(this%rwHandle) 1187 call cleanShutdown("Finished dense matrix write") 1188 end if 1189 call elsi_dm_real(this%handle, HSqrReal, SSqrReal, rhoSqrReal, Eband(iS)) 1190 1191 call packRhoRealBlacs(env%blacs, denseDesc, rhoSqrReal, neighbourList%iNeighbour,& 1192 & nNeighbourSK, orb%mOrb, iSparseStart, img2CentCell, rhoPrim(:,iS)) 1193 1194 end subroutine getDensityRealDense 1195 1196 1197 !> Returns the density matrix using ELSI non-diagonalisation routines (complex dense case). 1198 subroutine getDensityCmplxDense(this, env, denseDesc, ham, over, neighbourList,& 1199 & nNeighbourSK, iSparseStart, img2CentCell, iCellVec, cellVec, kPoint, kWeight, orb,& 1200 & parallelKS, rhoPrim, Eband, HSqrCplx, SSqrCplx) 1201 1202 !> Electronic solver information 1203 type(TElsiSolver), intent(inout) :: this 1204 1205 !> Environment settings 1206 type(TEnvironment), intent(inout) :: env 1207 1208 !> Dense matrix descriptor 1209 type(TDenseDescr), intent(in) :: denseDesc 1210 1211 !> hamiltonian in sparse storage 1212 real(dp), intent(in) :: ham(:,:) 1213 1214 !> sparse overlap matrix 1215 real(dp), intent(in) :: over(:) 1216 1217 !> list of neighbours for each atom 1218 type(TNeighbourList), intent(in) :: neighbourList 1219 1220 !> Number of neighbours for each of the atoms 1221 integer, intent(in) :: nNeighbourSK(:) 1222 1223 !> Index array for the start of atomic blocks in sparse arrays 1224 integer, intent(in) :: iSparseStart(:,:) 1225 1226 !> map from image atoms to the original unique atom 1227 integer, intent(in) :: img2CentCell(:) 1228 1229 !> Index for which unit cell atoms are associated with 1230 integer, intent(in) :: iCellVec(:) 1231 1232 !> Vectors (in units of the lattice constants) to cells of the lattice 1233 real(dp), intent(in) :: cellVec(:,:) 1234 1235 !> k-points 1236 real(dp), intent(in) :: kPoint(:,:) 1237 1238 !> Weights for k-points 1239 real(dp), intent(in) :: kWeight(:) 1240 1241 !> Atomic orbital information 1242 type(TOrbitals), intent(in) :: orb 1243 1244 !> K-points and spins to process 1245 type(TParallelKS), intent(in) :: parallelKS 1246 1247 !> sparse density matrix 1248 real(dp), intent(inout) :: rhoPrim(:,:) 1249 1250 !> band structure energy 1251 real(dp), intent(out) :: Eband(:) 1252 1253 !> electronic entropy times temperature 1254 !> dense complex (k-points) hamiltonian storage 1255 complex(dp), intent(inout), allocatable :: HSqrCplx(:,:) 1256 1257 !> dense complex (k-points) overlap storage 1258 complex(dp), intent(inout), allocatable :: SSqrCplx(:,:) 1259 1260 complex(dp), allocatable :: rhoSqrCplx(:,:) 1261 integer :: iKS, iK, iS 1262 1263 ! as each spin and k-point combination forms a separate group for this solver, then iKS = 1 1264 iKS = 1 1265 iK = parallelKS%localKS(1, iKS) 1266 iS = parallelKS%localKS(2, iKS) 1267 1268 HSqrCplx(:,:) = 0.0_dp 1269 call unpackHSCplxBlacs(env%blacs, ham(:,iS), kPoint(:,iK), neighbourList%iNeighbour,& 1270 & nNeighbourSK, iCellVec, cellVec, iSparseStart, img2CentCell, denseDesc, HSqrCplx) 1271 if (.not. this%tCholeskyDecomposed) then 1272 SSqrCplx(:,:) = 0.0_dp 1273 call unpackHSCplxBlacs(env%blacs, over, kPoint(:,iK), neighbourList%iNeighbour,& 1274 & nNeighbourSK, iCellVec, cellVec, iSparseStart, img2CentCell, denseDesc,& 1275 & SSqrCplx) 1276 if (this%ommCholesky .or. this%iSolver == electronicSolverTypes%pexsi) then 1277 this%tCholeskyDecomposed = .true. 1278 end if 1279 end if 1280 1281 allocate(rhoSqrCplx(size(HSqrCplx,dim=1),size(HSqrCplx,dim=2))) 1282 rhoSqrCplx(:,:) = 0.0_dp 1283 if (this%tWriteHS) then 1284 call elsi_write_mat_complex(this%rwHandle, "ELSI_Hcmplex.bin",& 1285 & HSqrCplx) 1286 call elsi_write_mat_complex(this%rwHandle, "ELSI_Scmplx.bin",& 1287 & SSqrCplx) 1288 call elsi_finalize_rw(this%rwHandle) 1289 call cleanShutdown("Finished dense matrix write") 1290 end if 1291 call elsi_dm_complex(this%handle, HSqrCplx, SSqrCplx, rhoSqrCplx,& 1292 & Eband(iS)) 1293 call packRhoCplxBlacs(env%blacs, denseDesc, rhoSqrCplx, kPoint(:,iK), kWeight(iK),& 1294 & neighbourList%iNeighbour, nNeighbourSK, orb%mOrb, iCellVec, cellVec,& 1295 & iSparseStart, img2CentCell, rhoPrim(:,iS)) 1296 1297 end subroutine getDensityCmplxDense 1298 1299 1300 !> Returns the density matrix using ELSI non-diagonalisation routines (Pauli dense case). 1301 subroutine getDensityPauliDense(this, env, denseDesc, ham, over, neighbourList,& 1302 & nNeighbourSK, iSparseStart, img2CentCell, iCellVec, cellVec, kPoint, kWeight, orb, species,& 1303 & tSpinOrbit, tDualSpinOrbit, tMulliken, parallelKS, energy, rhoPrim, Eband, iHam, xi,& 1304 & orbitalL, iRhoPrim, HSqrCplx, SSqrCplx) 1305 1306 !> Electronic solver information 1307 type(TElsiSolver), intent(inout) :: this 1308 1309 !> Environment settings 1310 type(TEnvironment), intent(inout) :: env 1311 1312 !> Dense matrix descriptor 1313 type(TDenseDescr), intent(in) :: denseDesc 1314 1315 !> hamiltonian in sparse storage 1316 real(dp), intent(in) :: ham(:,:) 1317 1318 !> sparse overlap matrix 1319 real(dp), intent(in) :: over(:) 1320 1321 !> list of neighbours for each atom 1322 type(TNeighbourList), intent(in) :: neighbourList 1323 1324 !> Number of neighbours for each of the atoms 1325 integer, intent(in) :: nNeighbourSK(:) 1326 1327 !> Index array for the start of atomic blocks in sparse arrays 1328 integer, intent(in) :: iSparseStart(:,:) 1329 1330 !> map from image atoms to the original unique atom 1331 integer, intent(in) :: img2CentCell(:) 1332 1333 !> Index for which unit cell atoms are associated with 1334 integer, intent(in) :: iCellVec(:) 1335 1336 !> Vectors (in units of the lattice constants) to cells of the lattice 1337 real(dp), intent(in) :: cellVec(:,:) 1338 1339 !> k-points 1340 real(dp), intent(in) :: kPoint(:,:) 1341 1342 !> Weights for k-points 1343 real(dp), intent(in) :: kWeight(:) 1344 1345 !> Atomic orbital information 1346 type(TOrbitals), intent(in) :: orb 1347 1348 !> species of all atoms in the system 1349 integer, intent(in) :: species(:) 1350 1351 !> Are spin orbit interactions present 1352 logical, intent(in) :: tSpinOrbit 1353 1354 !> Are block population spin orbit interactions present 1355 logical, intent(in) :: tDualSpinOrbit 1356 1357 !> Should Mulliken populations be generated/output 1358 logical, intent(in) :: tMulliken 1359 1360 !> K-points and spins to process 1361 type(TParallelKS), intent(in) :: parallelKS 1362 1363 !> Energy contributions and total 1364 type(TEnergies), intent(inout) :: energy 1365 1366 !> sparse density matrix 1367 real(dp), intent(inout) :: rhoPrim(:,:) 1368 1369 !> band structure energy 1370 real(dp), intent(out) :: Eband(:) 1371 1372 !> imaginary part of hamitonian 1373 real(dp), intent(in), allocatable :: iHam(:,:) 1374 1375 !> spin orbit constants 1376 real(dp), intent(in), allocatable :: xi(:,:) 1377 1378 !> orbital moments of atomic shells 1379 real(dp), intent(inout), allocatable :: orbitalL(:,:,:) 1380 1381 !> imaginary part of density matrix 1382 real(dp), intent(inout), allocatable :: iRhoPrim(:,:) 1383 1384 !> dense complex (k-points) hamiltonian storage 1385 complex(dp), intent(inout), allocatable :: HSqrCplx(:,:) 1386 1387 !> dense complex (k-points) overlap storage 1388 complex(dp), intent(inout), allocatable :: SSqrCplx(:,:) 1389 1390 complex(dp), allocatable :: rhoSqrCplx(:,:) 1391 real(dp), allocatable :: rVecTemp(:), orbitalLPart(:,:,:) 1392 integer :: iKS, iK, iS 1393 integer :: nAtom 1394 1395 ! as each spin and k-point combination forms a separate group for this solver, then iKS = 1 1396 iKS = 1 1397 iK = parallelKS%localKS(1, iKS) 1398 iS = parallelKS%localKS(2, iKS) 1399 1400 nAtom = size(nNeighbourSK) 1401 1402 if (tSpinOrbit .and. .not. tDualSpinOrbit) then 1403 energy%atomLS(:) = 0.0_dp 1404 allocate(rVecTemp(nAtom)) 1405 end if 1406 1407 if (tSpinOrbit .and. .not. tDualSpinOrbit) then 1408 allocate(orbitalLPart(3, orb%mShell, nAtom)) 1409 orbitalL(:,:,:) = 0.0_dp 1410 end if 1411 1412 if (allocated(iHam)) then 1413 call unpackHPauliBlacs(env%blacs, ham, kPoint(:,iK), neighbourList%iNeighbour,& 1414 & nNeighbourSK, iCellVec, cellVec, iSparseStart, img2CentCell, orb%mOrb, denseDesc,& 1415 & HSqrCplx, iorig=iHam) 1416 else 1417 call unpackHPauliBlacs(env%blacs, ham, kPoint(:,iK), neighbourList%iNeighbour,& 1418 & nNeighbourSK, iCellVec, cellVec, iSparseStart, img2CentCell, orb%mOrb, denseDesc,& 1419 & HSqrCplx) 1420 end if 1421 if (.not. this%tCholeskyDecomposed) then 1422 SSqrCplx(:,:) = 0.0_dp 1423 call unpackSPauliBlacs(env%blacs, over, kPoint(:,iK), neighbourList%iNeighbour,& 1424 & nNeighbourSK, iCellVec, cellVec, iSparseStart, img2CentCell, orb%mOrb, denseDesc,& 1425 & SSqrCplx) 1426 if (this%ommCholesky .or. this%iSolver == electronicSolverTypes%pexsi) then 1427 this%tCholeskyDecomposed = .true. 1428 end if 1429 end if 1430 if (allocated(xi) .and. .not. allocated(iHam)) then 1431 call addOnsiteSpinOrbitHam(env, xi, species, orb, denseDesc, HSqrCplx) 1432 end if 1433 allocate(rhoSqrCplx(size(HSqrCplx,dim=1),size(HSqrCplx,dim=2))) 1434 rhoSqrCplx(:,:) = 0.0_dp 1435 if (this%tWriteHS) then 1436 call elsi_write_mat_complex(this%rwHandle, "ELSI_Hcmplex.bin",& 1437 & HSqrCplx) 1438 call elsi_write_mat_complex(this%rwHandle, "ELSI_Scmplx.bin",& 1439 & SSqrCplx) 1440 call elsi_finalize_rw(this%rwHandle) 1441 call cleanShutdown("Finished dense matrix write") 1442 end if 1443 call elsi_dm_complex(this%handle, HSqrCplx, SSqrCplx, rhoSqrCplx,& 1444 & Eband(iS)) 1445 if (tSpinOrbit .and. .not. tDualSpinOrbit) then 1446 call getOnsiteSpinOrbitEnergy(env, rVecTemp, rhoSqrCplx, denseDesc, xi, orb, species) 1447 energy%atomLS = energy%atomLS + kWeight(iK) * rVecTemp 1448 if (tMulliken) then 1449 orbitalLPart(:,:,:) = 0.0_dp 1450 call getLOnsite(env, orbitalLPart, rhoSqrCplx, denseDesc, orb, species) 1451 orbitalL(:,:,:) = orbitalL + kWeight(iK) * orbitalLPart 1452 end if 1453 end if 1454 if (allocated(iHam)) then 1455 call packRhoPauliBlacs(env%blacs, denseDesc, rhoSqrCplx, kPoint(:,iK), kWeight(iK),& 1456 & neighbourList%iNeighbour, nNeighbourSK, orb%mOrb, iCellVec, cellVec,& 1457 & iSparseStart, img2CentCell, rhoPrim, iRhoPrim) 1458 iRhoPrim(:,:) = 2.0_dp * iRhoPrim 1459 else 1460 call packRhoPauliBlacs(env%blacs, denseDesc, rhoSqrCplx, kPoint(:,iK), kWeight(iK),& 1461 & neighbourList%iNeighbour, nNeighbourSK, orb%mOrb, iCellVec, cellVec,& 1462 & iSparseStart, img2CentCell, rhoPrim) 1463 end if 1464 RhoPrim(:,:) = 2.0_dp * RhoPrim 1465 1466 end subroutine getDensityPauliDense 1467 1468 1469 !> Calculates density matrix using the elsi routine. 1470 subroutine getDensityRealSparse(this, parallelKS, ham, over, iNeighbour, nNeighbourSK,& 1471 & iAtomStart, iSparseStart, img2CentCell, orb, rho, Eband) 1472 1473 !> Electronic solver information 1474 type(TElsiSolver), intent(inout) :: this 1475 1476 !> Contains (iK, iS) tuples to be processed in parallel by various processor groups 1477 type(TParallelKS), intent(in) :: parallelKS 1478 1479 !> hamiltonian in sparse storage 1480 real(dp), intent(in) :: ham(:,:) 1481 1482 !> overlap matrix in sparse storage 1483 real(dp), intent(in) :: over(:) 1484 1485 !> Neighbour list for the atoms (First index from 0!) 1486 integer, intent(in) :: iNeighbour(0:,:) 1487 1488 !> Nr. of neighbours for the atoms. 1489 integer, intent(in) :: nNeighbourSK(:) 1490 1491 !> Atom offset for the squared matrix 1492 integer, intent(in) :: iAtomStart(:) 1493 1494 !> indexing array for the sparse Hamiltonian 1495 integer, intent(in) :: iSparseStart(0:,:) 1496 1497 !> Mapping between image atoms and corresponding atom in the central cell. 1498 integer, intent(in) :: img2CentCell(:) 1499 1500 !> data structure with atomic orbital information 1501 type(TOrbitals), intent(in) :: orb 1502 1503 !> Density matrix in DFTB+ sparse format 1504 real(dp), intent(out) :: rho(:,:) 1505 1506 !> Band energy 1507 real(dp), intent(out) :: Eband(:) 1508 1509 integer :: iS 1510 1511 real(dp), allocatable :: HnzValLocal(:), SnzValLocal(:) 1512 real(dp), allocatable :: DMnzValLocal(:) 1513 1514 allocate(HnzValLocal(this%elsiCsc%nnzLocal)) 1515 allocate(SnzValLocal(this%elsiCsc%nnzLocal)) 1516 allocate(DMnzValLocal(this%elsiCsc%nnzLocal)) 1517 1518 call this%elsiCsc%convertPackedToElsiReal(over, iNeighbour, nNeighbourSK, iAtomStart,& 1519 & iSparseStart, img2CentCell, SnzValLocal) 1520 if (this%tFirstCalc) then 1521 call elsi_set_csc(this%handle, this%elsiCsc%nnzGlobal, this%elsiCsc%nnzLocal,& 1522 & this%elsiCsc%numColLocal, this%elsiCsc%rowIndLocal, this%elsiCsc%colPtrLocal) 1523 this%tFirstCalc = .false. 1524 end if 1525 1526 iS = parallelKS%localKS(2, 1) 1527 call this%elsiCsc%convertPackedToElsiReal(ham(:,iS), iNeighbour, nNeighbourSK, iAtomStart,& 1528 & iSparseStart, img2CentCell, HnzValLocal) 1529 1530 if (this%tWriteHS) then 1531 call elsi_set_rw_csc(this%rwHandle, this%elsiCsc%nnzGlobal,& 1532 & this%elsiCsc%nnzLocal, this%elsiCsc%numColLocal) 1533 call elsi_write_mat_real_sparse(this%rwHandle, "ELSI_HrealSparse.bin",& 1534 & this%elsiCsc%rowIndLocal, this%elsiCsc%colPtrLocal, HnzValLocal) 1535 call elsi_write_mat_real_sparse(this%rwHandle, "ELSI_SrealSparse.bin",& 1536 & this%elsiCsc%rowIndLocal, this%elsiCsc%colPtrLocal, SnzValLocal) 1537 call elsi_finalize_rw(this%rwHandle) 1538 call cleanShutdown("Finished matrix write") 1539 end if 1540 1541 call elsi_dm_real_sparse(this%handle, HnzValLocal, SnzValLocal, DMnzValLocal,& 1542 & Eband(iS)) 1543 1544 rho(:,:) = 0.0_dp 1545 call this%elsiCsc%convertElsiToPackedReal(iNeighbour, nNeighbourSK, orb%mOrb,& 1546 & iAtomStart, iSparseStart, img2CentCell, DMnzValLocal, rho(:,iS)) 1547 1548 end subroutine getDensityRealSparse 1549 1550 1551 !> Calculates density matrix using the elsi routine. 1552 subroutine getDensityCmplxSparse(this, parallelKS, kPoint, kWeight, iCellVec, cellVec, ham,& 1553 & over, iNeighbour, nNeighbourSK, iAtomStart, iSparseStart, img2CentCell, orb, rho, Eband) 1554 1555 !> Electronic solver information 1556 type(TElsiSolver), intent(inout) :: this 1557 1558 !> Contains (iK, iS) tuples to be processed in parallel by various processor groups 1559 type(TParallelKS), intent(in) :: parallelKS 1560 1561 !> Current k-point 1562 real(dp), intent(in) :: kPoint(:) 1563 1564 !> Weight for current k-points 1565 real(dp), intent(in) :: kWeight 1566 1567 !> Index for which unit cell atoms are associated with 1568 integer, intent(in) :: iCellVec(:) 1569 1570 !> Vectors (in units of the lattice constants) to cells of the lattice 1571 real(dp), intent(in) :: cellVec(:,:) 1572 1573 !> hamiltonian in sparse storage 1574 real(dp), intent(in) :: ham(:,:) 1575 1576 !> overlap matrix in sparse storage 1577 real(dp), intent(in) :: over(:) 1578 1579 !> Neighbour list for the atoms (First index from 0!) 1580 integer, intent(in) :: iNeighbour(0:,:) 1581 1582 !> Nr. of neighbours for the atoms. 1583 integer, intent(in) :: nNeighbourSK(:) 1584 1585 !> Atom offset for the squared matrix 1586 integer, intent(in) :: iAtomStart(:) 1587 1588 !> indexing array for the sparse Hamiltonian 1589 integer, intent(in) :: iSparseStart(0:,:) 1590 1591 !> Mapping between image atoms and corresponding atom in the central cell. 1592 integer, intent(in) :: img2CentCell(:) 1593 1594 !> data structure with atomic orbital information 1595 type(TOrbitals), intent(in) :: orb 1596 1597 !> Density matrix in DFTB+ sparse format 1598 real(dp), intent(out) :: rho(:,:) 1599 1600 !> Band energy 1601 real(dp), intent(out) :: Eband(:) 1602 1603 integer :: iS 1604 1605 complex(dp), allocatable :: HnzValLocal(:), SnzValLocal(:) 1606 complex(dp), allocatable :: DMnzValLocal(:) 1607 1608 allocate(HnzValLocal(this%elsiCsc%nnzLocal)) 1609 allocate(SnzValLocal(this%elsiCsc%nnzLocal)) 1610 allocate(DMnzValLocal(this%elsiCsc%nnzLocal)) 1611 1612 call this%elsiCsc%convertPackedToElsiCmplx(over, iNeighbour, nNeighbourSK, iAtomStart,& 1613 & iSparseStart, img2CentCell, kPoint, kWeight, iCellVec, cellVec, SnzValLocal) 1614 if (this%tFirstCalc) then 1615 call elsi_set_csc(this%handle, this%elsiCsc%nnzGlobal, this%elsiCsc%nnzLocal,& 1616 & this%elsiCsc%numColLocal, this%elsiCsc%rowIndLocal, this%elsiCsc%colPtrLocal) 1617 this%tFirstCalc = .false. 1618 end if 1619 1620 iS = parallelKS%localKS(2, 1) 1621 call this%elsiCsc%convertPackedToElsiCmplx(ham(:,iS), iNeighbour, nNeighbourSK, iAtomStart,& 1622 & iSparseStart, img2CentCell, kPoint, kWeight, iCellVec, cellVec, HnzValLocal) 1623 1624 if (this%tWriteHS) then 1625 call elsi_set_rw_csc(this%rwHandle, this%elsiCsc%nnzGlobal, this%elsiCsc%nnzLocal,& 1626 & this%elsiCsc%numColLocal) 1627 call elsi_write_mat_complex_sparse(this%rwHandle, "ELSI_HcmplxSparse.bin",& 1628 & this%elsiCsc%rowIndLocal, this%elsiCsc%colPtrLocal, HnzValLocal) 1629 call elsi_write_mat_complex_sparse(this%rwHandle, "ELSI_ScmplxSparse.bin",& 1630 & this%elsiCsc%rowIndLocal, this%elsiCsc%colPtrLocal, SnzValLocal) 1631 call elsi_finalize_rw(this%rwHandle) 1632 call cleanShutdown("Finished matrix write") 1633 end if 1634 1635 call elsi_dm_complex_sparse(this%handle, HnzValLocal, SnzValLocal, DMnzValLocal, Eband(iS)) 1636 1637 rho(:,:) = 0.0_dp 1638 call this%elsiCsc%convertElsiToPackedCmplx(iNeighbour, nNeighbourSK, orb%mOrb, iAtomStart,& 1639 & iSparseStart, img2CentCell, kPoint, kWeight, iCellVec, cellVec, DMnzValLocal, rho(:,iS)) 1640 1641 end subroutine getDensityCmplxSparse 1642 1643 1644 !> Returns the energy weighted density matrix using ELSI non-diagonalisation routines. 1645 subroutine getEDensityMtxReal(this, env, denseDesc, neighbourList, nNeighbourSK, orb,& 1646 & iSparseStart, img2CentCell, ERhoPrim, SSqrReal) 1647 1648 !> Electronic solver information 1649 type(TElsiSolver), intent(inout) :: this 1650 1651 !> Environment settings 1652 type(TEnvironment), intent(inout) :: env 1653 1654 !> Dense matrix descriptor 1655 type(TDenseDescr), intent(in) :: denseDesc 1656 1657 !> list of neighbours for each atom 1658 type(TNeighbourList), intent(in) :: neighbourList 1659 1660 !> Number of neighbours for each of the atoms 1661 integer, intent(in) :: nNeighbourSK(:) 1662 1663 !> Atomic orbital information 1664 type(TOrbitals), intent(in) :: orb 1665 1666 !> Index array for the start of atomic blocks in sparse arrays 1667 integer, intent(in) :: iSparseStart(:,:) 1668 1669 !> map from image atoms to the original unique atom 1670 integer, intent(in) :: img2CentCell(:) 1671 1672 !> Energy weighted sparse matrix 1673 real(dp), intent(out) :: ERhoPrim(:) 1674 1675 !> Storage for dense overlap matrix 1676 real(dp), intent(inout), allocatable :: SSqrReal(:,:) 1677 1678 real(dp), allocatable :: EDMnzValLocal(:) 1679 1680 ERhoPrim(:) = 0.0_dp 1681 if (this%isSparse) then 1682 allocate(EDMnzValLocal(this%elsiCsc%nnzLocal)) 1683 call elsi_get_edm_real_sparse(this%handle, EDMnzValLocal) 1684 call this%elsiCsc%convertElsiToPackedReal(neighbourList%iNeighbour, nNeighbourSK, orb%mOrb,& 1685 & denseDesc%iAtomStart, iSparseStart, img2CentCell, EDMnzValLocal, ErhoPrim) 1686 else 1687 call elsi_get_edm_real(this%handle, SSqrReal) 1688 call packRhoRealBlacs(env%blacs, denseDesc, SSqrReal, neighbourList%iNeighbour, nNeighbourSK,& 1689 & orb%mOrb, iSparseStart, img2CentCell, ERhoPrim) 1690 end if 1691 1692 ! add contributions from different spin channels together if necessary 1693 call mpifx_allreduceip(env%mpi%globalComm, ERhoPrim, MPI_SUM) 1694 1695 end subroutine getEDensityMtxReal 1696 1697 1698 !> Returns the energy weighted density matrix using ELSI non-diagonalisation routines. 1699 subroutine getEDensityMtxCmplx(this, env, denseDesc, kPoint, kWeight, neighbourList,& 1700 & nNeighbourSK, orb, iSparseStart, img2CentCell, iCellVec, cellVec, parallelKS, ERhoPrim,& 1701 & SSqrCplx) 1702 1703 !> Electronic solver information 1704 type(TElsiSolver), intent(inout) :: this 1705 1706 !> Environment settings 1707 type(TEnvironment), intent(inout) :: env 1708 1709 !> Dense matrix descriptor 1710 type(TDenseDescr), intent(in) :: denseDesc 1711 1712 !> K-points 1713 real(dp), intent(in) :: kPoint(:,:) 1714 1715 !> Weights for k-points 1716 real(dp), intent(in) :: kWeight(:) 1717 1718 !> list of neighbours for each atom 1719 type(TNeighbourList), intent(in) :: neighbourList 1720 1721 !> Number of neighbours for each of the atoms 1722 integer, intent(in) :: nNeighbourSK(:) 1723 1724 !> Atomic orbital information 1725 type(TOrbitals), intent(in) :: orb 1726 1727 !> Index array for the start of atomic blocks in sparse arrays 1728 integer, intent(in) :: iSparseStart(:,:) 1729 1730 !> map from image atoms to the original unique atom 1731 integer, intent(in) :: img2CentCell(:) 1732 1733 !> Index for which unit cell atoms are associated with 1734 integer, intent(in) :: iCellVec(:) 1735 1736 !> Vectors (in units of the lattice constants) to cells of the lattice 1737 real(dp), intent(in) :: cellVec(:,:) 1738 1739 !> K-points and spins to process 1740 type(TParallelKS), intent(in) :: parallelKS 1741 1742 !> Energy weighted sparse matrix 1743 real(dp), intent(out) :: ERhoPrim(:) 1744 1745 !> Storage for dense overlap matrix (complex case) 1746 complex(dp), intent(inout), allocatable :: SSqrCplx(:,:) 1747 1748 complex(dp), allocatable :: EDMnzValLocal(:) 1749 integer :: iK 1750 1751 ERhoPrim(:) = 0.0_dp 1752 ! iKS always 1 1753 iK = parallelKS%localKS(1, 1) 1754 if (this%isSparse) then 1755 allocate(EDMnzValLocal(this%elsiCsc%nnzLocal)) 1756 call elsi_get_edm_complex_sparse(this%handle, EDMnzValLocal) 1757 call this%elsiCsc%convertElsiToPackedCmplx(neighbourList%iNeighbour, nNeighbourSK, orb%mOrb,& 1758 & denseDesc%iAtomStart, iSparseStart, img2CentCell, kPoint(:,iK), kWeight(iK), iCellVec,& 1759 & cellVec, EDMnzValLocal, ERhoPrim) 1760 else 1761 call elsi_get_edm_complex(this%handle, SSqrCplx) 1762 call packRhoCplxBlacs(env%blacs, denseDesc, SSqrCplx, kPoint(:,iK), kWeight(iK),& 1763 & neighbourList%iNeighbour, nNeighbourSK, orb%mOrb, iCellVec, cellVec, iSparseStart,& 1764 & img2CentCell, ERhoPrim) 1765 end if 1766 call mpifx_allreduceip(env%mpi%globalComm, ERhoPrim, MPI_SUM) 1767 1768 end subroutine getEDensityMtxCmplx 1769 1770 1771 !> Returns the energy weighted density matrix using ELSI non-diagonalisation routines. 1772 subroutine getEDensityMtxPauli(this, env, denseDesc, kPoint, kWeight, neighbourList,& 1773 & nNeighbourSK, orb, iSparseStart, img2CentCell, iCellVec, cellVec, parallelKS, ERhoPrim,& 1774 & SSqrCplx) 1775 1776 !> Electronic solver information 1777 type(TElsiSolver), intent(inout) :: this 1778 1779 !> Environment settings 1780 type(TEnvironment), intent(inout) :: env 1781 1782 !> Dense matrix descriptor 1783 type(TDenseDescr), intent(in) :: denseDesc 1784 1785 !> K-points 1786 real(dp), intent(in) :: kPoint(:,:) 1787 1788 !> Weights for k-points 1789 real(dp), intent(in) :: kWeight(:) 1790 1791 !> list of neighbours for each atom 1792 type(TNeighbourList), intent(in) :: neighbourList 1793 1794 !> Number of neighbours for each of the atoms 1795 integer, intent(in) :: nNeighbourSK(:) 1796 1797 !> Atomic orbital information 1798 type(TOrbitals), intent(in) :: orb 1799 1800 !> Index array for the start of atomic blocks in sparse arrays 1801 integer, intent(in) :: iSparseStart(:,:) 1802 1803 !> map from image atoms to the original unique atom 1804 integer, intent(in) :: img2CentCell(:) 1805 1806 !> Index for which unit cell atoms are associated with 1807 integer, intent(in) :: iCellVec(:) 1808 1809 !> Vectors (in units of the lattice constants) to cells of the lattice 1810 real(dp), intent(in) :: cellVec(:,:) 1811 1812 !> K-points and spins to process 1813 type(TParallelKS), intent(in) :: parallelKS 1814 1815 !> Energy weighted sparse matrix 1816 real(dp), intent(out) :: ERhoPrim(:) 1817 1818 !> Storage for dense overlap matrix 1819 complex(dp), intent(inout), allocatable :: SSqrCplx(:,:) 1820 1821 integer :: iK 1822 1823 if (this%isSparse) then 1824 call error("EDensity via sparse ELSI solver not supported for Pauli-matrices") 1825 else 1826 ! iKS always 1, as number of groups matches the number of k-points 1827 iK = parallelKS%localKS(1, 1) 1828 call elsi_get_edm_complex(this%handle, SSqrCplx) 1829 call packERhoPauliBlacs(env%blacs, denseDesc, SSqrCplx, kPoint(:,iK), kWeight(iK),& 1830 & neighbourList%iNeighbour, nNeighbourSK, orb%mOrb, iCellVec, cellVec, iSparseStart,& 1831 & img2CentCell, ERhoPrim) 1832 call mpifx_allreduceip(env%mpi%globalComm, ERhoPrim, MPI_SUM) 1833 end if 1834 1835 end subroutine getEDensityMtxPauli 1836 1837#:endif 1838 1839end module dftbp_elsisolver 1840