1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 module writewave 9 10 USE precision, only : dp 11 use alloc, only : re_alloc, de_alloc 12 use sys, only: die, message 13 implicit none 14 15 private 16 character(len=192), save, public :: wfs_filename 17 18 integer, save :: nwf = 1 19 logical, save :: wwf = .false. 20 logical, save :: scf_set = .false. 21 logical, save :: debug = .false. 22 23 24 integer, save, dimension(:), pointer :: nwflist => null() 25 integer, save, dimension(:,:), pointer :: iwf => null() 26 27 integer, public, save :: nwk 28 real(dp), public, pointer, save :: wfk(:,:) => null() 29 logical, public, save :: gamma_wavefunctions 30 31 real(dp), save :: wfs_energy_min = -huge(1.0_dp) 32 real(dp), save :: wfs_energy_max = huge(1.0_dp) 33 logical, save :: wfs_energy_window = .true. 34 35 public :: setup_wfs_list, setup_wf_kpoints, wwave, writew 36 37 CONTAINS 38 39 subroutine setup_wf_kpoints() 40 USE alloc, only: re_alloc 41 USE sys, only: die 42 USE atomlist, only: no_u 43 USE m_spin, only: spin 44 45 implicit none 46 47 integer :: max_n_states 48 49 ! Find number of k-points for wavefunction printout 50 51 nullify(wfk) 52 53 nwk = 0 54 if ( spin%Grid == 4 ) then 55 max_n_states = 2 * no_u 56 else 57 max_n_states = no_u 58 end if 59 call initwave( max_n_states, nwk, wfk ) 60 61 if (nwk .eq. 0) then 62 gamma_wavefunctions = .true. 63 else if (nwk .eq. 1) then 64 if (dot_product(wfk(:,1),wfk(:,1)) < 1.0e-20_dp) then 65 gamma_wavefunctions = .true. 66 else 67 gamma_wavefunctions = .false. 68 endif 69 else 70 gamma_wavefunctions = .false. 71 endif 72 73 end subroutine setup_wf_kpoints 74 75 subroutine setup_wfs_list(nk,norb,wfs_band_min,wfs_band_max, 76 $ use_scf_weights,use_energy_window) 77 use fdf 78 ! 79 ! Initialize structures if the number of bands 80 ! to output is the same for all k-points 81 ! (e.g., for the SCF k-point set, or for the new 'bands' option 82 83 integer, intent(in) :: nk, norb 84 integer, intent(in) :: wfs_band_min, wfs_band_max ! Band range 85 logical, intent(in) :: use_scf_weights 86 logical, intent(in) :: use_energy_window 87 88 integer :: j, nwfs 89 90 call re_alloc(nwflist, 1, nk, name="nwflist", 91 $ routine="setup_wfs_list") 92 call re_alloc(iwf, 1, nk, 1, norb, name="iwf", 93 $ routine="setup_wfs_list") 94 95 scf_set = use_scf_weights 96 wfs_energy_window = use_energy_window 97! 98 debug = fdf_get("WriteWaveDebug", .false.) 99 100! Here we request output of all the eigenstates in a k-point 101! except if the user selects a band or energy range. 102! 103 if (wfs_band_min < 1) then 104 call die("WFS.BandMin < 1") 105 endif 106 if (wfs_band_max > norb) then 107 call die("WFS.BandMax > maximum number of states") 108 endif 109 if (wfs_band_max < wfs_band_min) then 110 call message("WARNING","WFS.BandMax < WFS.BandMin." // 111 $ " Nothing written to file") 112 endif 113 114 nwfs = max(0,wfs_band_max - wfs_band_min + 1) 115 nwflist(1:nk) = nwfs 116 do j=1, nwfs 117 iwf(1:nk,j) = wfs_band_min - 1 + j 118 enddo 119 120! The energy window, if allowed, will take effect inside the routine 121! that writes to file, as the energies are not yet known 122 123 wfs_energy_min = fdf_get("WFS.EnergyMin", 124 $ -huge(1.0_dp),"Ry") 125 wfs_energy_max = fdf_get("WFS.EnergyMax", 126 $ huge(1.0_dp),"Ry") 127 128 end subroutine setup_wfs_list 129 130 131 subroutine initwave( norb, nk, kpoint ) 132C ********************************************************************* 133C Finds k-points for wavefunction printout 134C Based on initband routine by J.Soler 135C Written by P. Ordejon, June 2003 136C **************************** INPUT ********************************** 137C integer norb : Number of orbitals (actually max number of bands) 138C *************************** OUTPUT ********************************** 139C integer nk : Number k points to compute wavefunctions 140C real*8 kpoint(3,maxk) : k point vectors 141C *************************** UNITS *********************************** 142C Lengths in atomic units (Bohr). 143C k vectors in reciprocal atomic units. 144C ***************** BEHAVIOUR ***************************************** 145C - If nk=0 on input, k-points are read from labels WaveFuncKPoints and 146C WaveFuncKPointsScale from the input fdf data file. If these labels 147C are not present, it returns with nk=0. 148C - Allowed values for WaveFuncKPointsScale are ReciprocalLatticeVectors 149C and pi/a (default). If another value is given, it returns with nk=0 150C after printing a warning. 151C - If nk>maxk, k-points and wavefunctions are not calculated and no 152C warning is printed before return 153C ***************** USAGE ********************************************* 154C Example of fdf wavefunction k-points specification for an FCC lattice. 155C 156C WaveFuncKPointsScale pi/a 157C %block WaveFuncKPoints # These are comments 158C 0.000 0.000 0.000 from 1 to 10 # eigenstates 1-10 of Gamma 159C 2.000 0.000 0.000 1 3 5 # eigenstates 1,3,5 of X 160C 1.500 1.500 1.500 # all eigenstates of K 161C %endblock WaveFuncKPoints 162C 163C If only given points (not lines) are desired, simply specify 1 as 164C the number of points along the line. 165C ********************************************************************* 166C 167C Modules 168C 169 use precision 170 use parallel, only : Node 171 use fdf 172 use sys, only : die 173 use alloc, only : re_alloc 174 use m_get_kpoints_scale, only: get_kpoints_scale 175 176 implicit none 177 178 integer :: norb, nk 179 real(dp), pointer :: kpoint(:,:) 180 external :: memory 181C ********************************************************************* 182 183 character 184 . name1*10, name2*10, name3*10, 185 . msg*120 186 187 logical :: outlng 188 logical :: WaveFuncPresent 189 190 integer 191 . i, ix, iw, iw1, iw2, iw3, ierr, 192 . ni, nn, nr, nv 193 194 real(dp) 195 . caux(3,3), 196 . rcell(3,3) 197 198 type(block_fdf) :: bfdf 199 type(parsed_line), pointer :: pline 200 201 wfs_energy_window = .true. 202C Find if there are k-points data 203 outlng = fdf_boolean('LongOutput', .false.) 204 wwf = fdf_boolean('WriteWaveFunctions',outlng) 205 206C Find if there are k-points data 207 WaveFuncPresent = fdf_block('WaveFuncKPoints',bfdf) 208 if ( WaveFuncPresent ) then 209 210 call get_kpoints_scale('WaveFuncKPointsScale',rcell,ierr) 211 212 if (ierr /= 0) then 213 nk = 0 214 RETURN 215 endif 216 217C Find max number of k points and max number of bands per k-point 218C Loop on data lines 219 nk = 0 220 nwf = 0 221 do while(fdf_bline(bfdf,pline)) 222C Read and parse data line 223 nn = fdf_bnnames(pline) 224 nv = fdf_bnvalues(pline) 225 ni = fdf_bnintegers(pline) 226 nr = fdf_bnreals(pline) 227 228C Check if syntax line is correct 229 if (nv .ge. 3) then 230 231C Check syntax 232 if (nr .ne. 3) then 233 write(msg,'(a,/,a)') 234 . 'writewave: syntax ERROR in %block WaveFuncKPoints:', 235 . trim(fdf_getline(bfdf%mark)) 236 call die(msg) 237 endif 238 239C Add this point to total number of k points 240 nk = nk + 1 241 242C Find which eigenvectors should be printed 243C Store indexes of wave functions to printout 244 245C #X.XXX #Y.YYY #Z.ZZZ from #start to #end [step] #step 246 if (nn .ge. 1) then 247 248C Check that line contains 'from', 'to' and maybe 'step' 249 if ((nn .ne. 2) .and. (nn .ne. 3)) then 250 write(msg,'(a,/,a)') 251 . 'writewave: syntax ERROR in %block WaveFuncKPoints:', 252 . trim(fdf_getline(bfdf%mark)) 253 call die(msg) 254 endif 255 name1 = fdf_bnames(pline,1) 256 name2 = fdf_bnames(pline,2) 257 if (nn .eq. 3) name3 = fdf_bnames(pline,3) 258 if (name1 .ne. 'from' .or. name2 .ne. 'to') then 259 write(msg,'(a,/,a)') 260 . 'writewave: syntax ERROR in %block WaveFuncKPoints:', 261 . trim(fdf_getline(bfdf%mark)) 262 call die(msg) 263 endif 264 if (nn .eq. 3) then 265 if (name3 .ne. 'step') then 266 write(msg,'(a,/,a)') 267 . 'writewave: syntax ERROR in %block WaveFuncKPoints:', 268 . trim(fdf_getline(bfdf%mark)) 269 call die(msg) 270 endif 271 endif 272 273 iw1 = fdf_bintegers(pline,1) 274 iw2 = fdf_bintegers(pline,2) 275 if (iw1 .lt. 0) iw1 = norb + iw1 + 1 276 if (iw2 .lt. 0) iw2 = norb + iw2 + 1 277 if (nn .eq. 3) then 278 iw3 = abs(fdf_bintegers(pline,3)) 279 else 280 iw3 = 1 281 endif 282 ni = 0 283 do iw = min(iw1,iw2), max(iw1,iw2), iw3 284 ni = ni + 1 285 enddo 286 nwf = max(nwf,ni) 287 288C #X.XXX #Y.YYY #Z.ZZZ #eigen1 [#eigen2 ...] 289 elseif (ni .ne. 0) then 290 nwf = max(nwf,ni) 291 292C #X.XXX #Y.YYY #Z.ZZZ 293 elseif (ni .eq. 0) then 294 nwf = max(nwf,norb) 295 endif 296 else 297C Bad syntax 298 write(msg,'(a,/,a)') 299 . 'writewave: syntax ERROR in %block WaveFuncKPoints:', 300 . trim(fdf_getline(bfdf%mark)) 301 call die(msg) 302 endif 303 enddo 304 305 ! Minimize to the maximum number of eigenvalues 306 nwf = min(nwf, norb) 307 308C Allocate nwflist, iwf and kpoints structures according to nk and nwf 309 call re_alloc( nwflist, 1, nk, 'nwflist', 'initwave' ) 310 call re_alloc( iwf, 1, nk, 1, nwf, 'iwf', 'initwave' ) 311 312 call re_alloc( kpoint, 1, 3, 1, nk, 'kpoint', 'initwave' ) 313 314C Fill nwflist, iwf and kpoints structures 315C Loop on data lines 316 nk = 0 317 call fdf_brewind(bfdf) 318 do while(fdf_bline(bfdf,pline)) 319C Read and parse data line 320 nn = fdf_bnnames(pline) 321 nv = fdf_bnvalues(pline) 322 ni = fdf_bnintegers(pline) 323 nr = fdf_bnreals(pline) 324 325C Add this point to total number of k points 326 nk = nk + 1 327 328C Find coordinates of k point 329 330 do ix = 1,3 331 kpoint(ix,nk) = rcell(ix,1) * fdf_bvalues(pline,1) + 332 . rcell(ix,2) * fdf_bvalues(pline,2) + 333 . rcell(ix,3) * fdf_bvalues(pline,3) 334 enddo 335 336C Find which eigenvectors should be printed 337C Store indexes of wave functions to printout 338 339C #X.XXX #Y.YYY #Z.ZZZ from #start to #end [step] #step 340 if (nn .ge. 2) then 341 342 name1 = fdf_bnames(pline,1) 343 name2 = fdf_bnames(pline,2) 344 if (nn .eq. 3) name3 = fdf_bnames(pline,3) 345 346 iw1 = fdf_bintegers(pline,1) 347 iw2 = fdf_bintegers(pline,2) 348 if (iw1 .lt. 0) iw1 = norb + iw1 + 1 349 if (iw2 .lt. 0) iw2 = norb + iw2 + 1 350 ! Reduce to maximally the number of orbitals 351 iw2 = min(iw2, norb) 352 if (nn .eq. 3) then 353 iw3 = abs(fdf_bintegers(pline,3)) 354 else 355 iw3 = 1 356 endif 357 ni = 0 358 do iw = min(iw1,iw2), max(iw1,iw2), iw3 359 ni = ni + 1 360 if (iw .lt. 0) then 361 iwf(nk,ni) = norb + iw + 1 362 else 363 iwf(nk,ni) = iw 364 endif 365 enddo 366 367C #X.XXX #Y.YYY #Z.ZZZ #eigen1 [#eigen2 ...] 368 elseif (ni .ne. 0) then 369 do i= 1, ni 370 iw = fdf_bintegers(pline,i) 371 if (iw .lt. 0) then 372 iw = norb + iw + 1 373 endif 374 if ( iw <= norb ) then 375 iwf(nk,i) = iw 376 end if 377 enddo 378 379C #X.XXX #Y.YYY #Z.ZZZ 380 elseif (ni .eq. 0) then 381 ni = norb 382 do i= 1, ni 383 iwf(nk,i) = i 384 enddo 385 endif 386 387 nwflist(nk) = ni 388 enddo 389 390 endif 391 392 end subroutine initwave 393 394 subroutine wwave( no, spin, maxo, maxuo, maxnh, 395 . maxk, 396 . numh, listhptr, listh, H, S, ef, xij, indxuo, 397 . gamma, nk, kpoint, nuotot, occtol) 398C ********************************************************************* 399C Finds wavefunctions at selected k-points. 400C Written by P. Ordejon, June 2003 401C from routine 'bands' written by J.M.Soler 402C Changed to use tSpin, N. Papior, August 2019 403C **************************** INPUT ********************************** 404C integer no : Number of basis orbitals 405C tSpin spin : Contains spin information 406C integer maxo : First dimension of ek 407C integer maxuo : Second dimension of H and S 408C integer maxnh : Maximum number of orbitals interacting 409C with any orbital 410C integer maxk : Last dimension of kpoint and ek 411C integer numh(nuo) : Number of nonzero elements of each row 412C of hamiltonian matrix 413C integer listhptr(nuo) : Pointer to start of each row of the 414C hamiltonian matrix 415C integer listh(maxlh) : Nonzero hamiltonian-matrix element 416C column indexes for each matrix row 417C real*8 H(maxnh,spin%H) : Hamiltonian in sparse form 418C real*8 S(maxnh) : Overlap in sparse form 419C real*8 ef : Fermi energy 420C real*8 xij(3,maxnh) : Vectors between orbital centers (sparse) 421C (not used if only gamma point) 422C integer indxuo(no) : Index of equivalent orbital in unit cell 423C Unit cell orbitals must be the first in 424C orbital lists, i.e. indxuo.le.nuo, with 425C nuo the number of orbitals in unit cell 426C logical Gamma : whether only the Gamma-point is sampled 427C integer nk : Number of band k points 428C real*8 kpoint(3,maxk) : k point vectors 429C integer nuotot : Total number of orbitals in unit cell 430C real*8 occtol : Occupancy threshold for DM build 431C *************************** OUTPUT ********************************** 432C None; output is dumped to wave functions file SystemLabel.WFSX 433C *************************** UNITS *********************************** 434C Lengths in atomic units (Bohr). 435C k vectors in reciprocal atomic units. 436C Energies in Rydbergs. 437C 438C Modules 439C 440 use precision 441 use parallel, only : Node, Nodes 442 use fdf 443 use densematrix, only : allocDenseMatrix, resetDenseMatrix 444 use densematrix, only : Haux, Saux, psi 445 use alloc 446 use atmfuncs, only : symfio, cnfigfio, labelfis, nofis 447 use atomlist, only : iaorb, iphorb 448 use siesta_geom, only : isa 449 use t_spin, only: tSpin 450 use units, only : eV 451 452#ifdef MPI 453 use parallelsubs, only : GetNodeOrbs 454#endif 455 use m_diag_option, only: Serial, ParallelOverK 456 use m_diag, only: diag_init 457 458 implicit none 459 460 type(tSpin), intent(in) :: spin 461 logical, intent(in) :: Gamma 462 integer maxk, maxnh, maxo, maxuo, nk, no, 463 . nuotot, indxuo(no), listh(maxnh), numh(*), 464 . listhptr(*) 465 real(dp) ef, H(maxnh,spin%H), kpoint(3,maxk), 466 . S(maxnh), xij(3,maxnh), occtol 467 468 external io_assign, io_close, memory 469C ********************************************************************* 470 471 logical 472 . fixspin 473 474 integer 475 . io, iu, iuo, naux, nuo, j, nhs 476 477 real(dp) 478 . Dnew, qs(2), e1, e2, efs(2), Enew, qk, qtot, 479 . temp, wk, Entropy 480 481C Dynamic arrays 482 logical, parameter :: getD = .false. 483 logical, parameter :: getPSI = .true. 484 real(dp), allocatable :: aux(:) 485 real(dp), allocatable :: ek(:,:,:) 486 487 save Dnew, Enew, e1, e2, qk, qtot, temp, wk 488 data Dnew, Enew, e1, e2, qk, qtot, temp, wk /8*0.d0/ 489 490 logical :: SaveParallelOverK 491 492 493C Get local number of orbitals 494#ifdef MPI 495 call GetNodeOrbs(nuotot,Node,Nodes,nuo) 496#else 497 nuo = nuotot 498#endif 499 500C Start time counter 501 call timer( 'writewave', 1 ) 502 503C Check parameter maxk 504 if (nk .gt. maxk) then 505 if (Node.eq.0) then 506 write(6,'(/,a,/,a)') 507 . 'writewave: WARNING: parameter maxk too small', 508 . 'writewave: No wavefunction calculation performed' 509 endif 510 goto 999 511 endif 512 513C Check spin 514 515C Allocate local arrays - only aux is relevant here 516 if ( spin%Grid == 4 ) then 517 nhs = 2 * (2*nuotot) * (2*nuo) 518 else 519 nhs = 2 * nuotot * nuo 520 endif 521 call allocDenseMatrix(nhs, nhs, nhs) 522 523 allocate(ek(nuotot,spin%spinor,nk)) 524 call memory('A','D',nuotot*spin%spinor*nk,'writewave') 525 naux = 2*nuotot*5 526 allocate(aux(naux)) 527 call memory('A','D',naux,'writewave') 528 529C Open file 530 531 if (Node.eq.0) then 532 533 call io_assign( iu ) 534 open(iu, file=wfs_filename,form="unformatted",status='unknown') 535 536 rewind (iu) 537 538 if (wwf) then 539 write(6,*) 540 write(6,'(a)') 'writewave: Wave Functions Coefficients' 541 write(6,*) 542 write(6,'(a26,2x,i6)') 'Number of k-points = ', nk 543 write(6,'(a26,2x,i6)') 'Number of Spins = ', spin%H 544 write(6,'(a26,2x,i6)') 'Number of basis orbs = ',nuotot 545 write(6,*) 546 endif 547 write(iu) nk, gamma 548 write(iu) spin%H ! This will be 1, 2, 4 or 8 (to tell NC from SOC) 549 write(iu) nuotot 550 write(iu) (iaorb(j),labelfis(isa(iaorb(j))), 551 . iphorb(j), cnfigfio(isa(iaorb(j)),iphorb(j)), 552 . symfio(isa(iaorb(j)),iphorb(j)), j=1,nuotot) 553 554 call io_close(iu) 555 556 endif 557 558C Find the eigenvectors 559c fixspin and qs are not used in diagk, since getD=.false. ... 560 fixspin = .false. 561 qs(1) = 0.0d0 562 qs(2) = 0.0d0 563 564C Handle parallel over K points option which is not allowed for here 565 SaveParallelOverK = ParallelOverK 566 if ( ParallelOverK ) then 567 if ( Node .eq. 0 ) then 568 write(*,"(a)") 569 $ "*** Note: ParallelOverK option not used for "// 570 & "WF info generation" 571 end if 572 ParallelOverK = .false. 573 Serial = (Nodes == 1) 574 call diag_init() 575 end if 576 577C Call appropriate diagonalization routine 578 579 if ( spin%none .or. spin%Col ) then 580 if (gamma) then 581 call diagg( spin%H, nuo, maxuo, maxnh, maxnh, 582 . maxo, numh, listhptr, listh, numh, listhptr, 583 . listh, H, S, getD, getPSI, 584 . fixspin, qtot, qs, temp, 585 . e1, e2, ek, qk, Dnew, Enew, ef, efs, Entropy, 586 . Haux, Saux, psi, nuotot, occtol, 1, nuotot ) 587 else 588 call diagk( spin%H, nuo, no, spin%spinor, maxnh, maxnh, 589 . maxo, numh, listhptr, listh, numh, listhptr, 590 . listh, H, S, getD, getPSI, 591 . fixspin, qtot, qs, temp, 592 . e1, e2, xij, indxuo, nk, kpoint, wk, 593 . ek, qk, Dnew, Enew, ef, efs, Entropy, 594 . Haux, Saux, psi, Haux, Saux, aux, nuotot, 595 . occtol, 1, nuotot ) 596 endif 597 598 elseif ( spin%NCol) then 599 if (gamma) then 600 call diag2g( nuo, no, maxnh, maxnh, maxo, numh, 601 . listhptr, listh, numh, listhptr, listh, 602 . H, S, getD, getPSI, qtot, temp, e1, e2, ek, qk, 603 . Dnew, Enew, ef, Entropy, 604 . Haux, Saux, psi, aux, 605 . nuotot, occtol, 1, nuotot ) 606 else 607 call diag2k( nuo, no, maxnh, maxnh, maxo, numh, 608 . listhptr, listh, numh, listhptr, listh, 609 . H, S, getD, getPSI, qtot, temp, e1, e2, xij, 610 . indxuo, nk, kpoint, wk, ek, qk, Dnew, 611 . Enew, ef, Entropy, 612 . Haux, Saux, psi, Haux, Saux, aux, 613 . nuotot, occtol, 1, nuotot ) 614 endif 615 elseif ( spin%SO ) then 616 if (gamma) then 617 call diag3g( nuo, no, maxnh, maxnh, maxo, numh, 618 . listhptr, listh, numh, listhptr, listh, 619 . H, S, getD, getPSI, qtot, temp, e1, e2, ek, qk, 620 . Dnew, Enew, ef, Entropy, 621 . Haux, Saux, psi, aux, 622 . nuotot, occtol, 1, nuotot ) 623 else 624 call diag3k( nuo, no, maxnh, maxnh, maxo, numh, 625 . listhptr, listh, numh, listhptr, listh, 626 . H, S, getD, getPsi, qtot, temp, e1, e2, xij, 627 . indxuo, nk, kpoint, wk, ek, qk, Dnew, 628 . Enew, ef, Entropy, 629 . Haux, Saux, psi, Haux, Saux, aux, 630 . nuotot, occtol, 1, nuotot ) 631 endif 632 endif 633 634 ! Restore ParallelOverK 635 ParallelOverK = SaveParallelOverK 636 if ( ParallelOverK ) Serial = .true. 637 638 call resetDenseMatrix() 639C Free local arrays 640 call memory('D','D',size(aux),'writewave') 641 deallocate(aux) 642 call memory('D','D',size(ek),'writewave') 643 deallocate(ek) 644 645C This is the only exit point 646 999 continue 647 call timer( 'writewave', 2 ) 648 649 end subroutine wwave 650 651 652 subroutine writew(nuotot,nuo,ik,k,ispin,eo,psi, 653 $ gamma,non_coll,blocksize) 654 655 use precision 656 use sys 657 use parallel, only : Node, Nodes 658 use parallel, only : Module_Blocksize=>BlockSize 659 use parallelsubs, only : GlobalToLocalOrb, WhichNodeOrb 660 use units, only : eV 661 use kpoint_grid, only : kweight 662#ifdef MPI 663 use mpi_siesta 664#endif 665 666 implicit none 667 668#ifdef MPI 669 integer MPIerror, mpistatus(MPI_STATUS_SIZE), tag 670#endif 671 672 integer, intent(in) :: nuotot, nuo, ispin, ik 673 logical, intent(in) :: non_coll, gamma 674 real(dp), intent(in) :: eo(*), k(3) 675 real(dp), intent(in) :: psi(:) 676 integer, intent(in) :: blocksize 677 678 679C Internal variables ............................................. 680 integer BNode, iie, iw, indwf, j, ind, iu 681 integer :: blocksize_saved 682 integer number_of_wfns 683 integer, allocatable :: ind_wfn(:) 684 real(dp) :: kpoint_weight 685 real(SP), dimension(:,:), allocatable :: aux !! NOTE SP 686 687 external io_assign, io_close 688 689C ................... 690 691C Allocate auxiliary arrays 692 693 if (non_coll) then 694 allocate(aux(4,nuotot)) 695 call memory('A','D',4*nuotot,'writewave') 696 else 697 if (gamma) then 698 allocate(aux(1,nuotot)) 699 call memory('A','D',nuotot,'writewave') 700 else 701 allocate(aux(2,nuotot)) 702 call memory('A','D',2*nuotot,'writewave') 703 endif 704 endif 705 706C Find file name, and open for Node 0 707 708 if (Node .eq. 0) then 709 call io_assign( iu ) 710 open(iu,file=wfs_filename,form="unformatted",position='append', 711 . status='old') 712 endif 713 714 if (wfs_energy_window) then 715 allocate(ind_wfn(nwflist(ik))) 716 number_of_wfns = 0 717 do iw = 1,nwflist(ik) 718 indwf = iwf(ik,iw) 719 if (eo(indwf) >= wfs_energy_min .AND. 720 $ eo(indwf) <= wfs_energy_max ) then 721 number_of_wfns = number_of_wfns + 1 722 ind_wfn(number_of_wfns) = indwf 723 endif 724 enddo 725 ! Replace indexes with new values 726 nwflist(ik) = number_of_wfns 727 do iw = 1, nwflist(ik) 728 iwf(ik,iw) = ind_wfn(iw) 729 enddo 730 deallocate(ind_wfn) 731 endif 732 733C First print the index and value of k-point 734 735 if (Node .eq. 0) then 736 if (wwf) then 737 write(6,*) 738 write(6,'(a22,2x,i6,2x,3f10.6)') 'k-point = ',ik, k(1:3) 739 if (non_coll) then 740 write(6,'(a22,2x)') 'Spinors have mixed spins' 741 else 742 write(6,'(a22,2x,i6)') 'Spin component = ',ispin 743 endif 744 write(6,'(a22,2x,i6)') 'Num. wavefunctions = ',nwflist(ik) 745 write(6,'(a)') 'Use readwfsx utility to print ' 746 $ // "wavefunction coefficients from WFSX file" 747 endif 748 if (scf_set) then 749 kpoint_weight = kweight(ik) 750 else 751 kpoint_weight = 1.0_dp 752 endif 753 write(iu) ik, k(1:3), kpoint_weight 754 write(iu) ispin 755 write(iu) nwflist(ik) 756 endif 757 758 ! Loop over wavefunctions that should be stored 759 760 ! The effective blocksize to use is passed as an argument. 761 ! Module_blocksize is the variable in the 'parallel' module that 762 ! holds the value that all the helper functions (WhichNodeOrb, 763 ! GetNodeOrb, etc) see. Hence, we must temporarily set 764 ! Module_blocksize to the passed blocksize. 765 766 ! Examples: non-collinear calculations need a wf blocksize double 767 ! the 'orbital' one 768 769 blocksize_saved = Module_blocksize 770 Module_blocksize = blocksize 771 772 if (Node .eq. 0 .and. debug) then 773 print *, " -*- Will print ", nwflist(ik), " wfns." 774 print *, " -*- ", iwf(ik,1:nwflist(ik)) 775 endif 776 do iw = 1,nwflist(ik) 777 indwf = iwf(ik,iw) 778 779 ! Determine which node handles this wavefunction 780 ! Note that the distribution is block cyclic, 781 ! so the same 'orbital' helper functions apply 782 783 call WhichNodeOrb(indwf,Nodes,BNode) 784 785 if (Node .eq. 0 .and. debug) then 786 print *, " -*- About to print wfn ", indwf, 787 $ "from node ", Bnode 788 endif 789 790 if (Node.eq.BNode) then 791 792 ! Determine the index of the wavefunction in the local node 793 call GlobalToLocalOrb( indwf, BNode, Nodes, iie) 794 795 ! Save wavefunction in aux array 796 797 ! Despite being passed as a one-dimensional array, psi has 798 ! different underlying structures for different cases, so the 799 ! indexing must be handled differently 800 801 if (.not. non_coll) then 802 if (gamma) then 803 ! Real functions 804 do j = 1,nuotot 805 ind = j + (iie-1)*nuotot 806 aux(1,j) = real(psi(ind),kind=sp) 807 enddo 808 else 809 ! Complex functions: two reals per orbital coefficient 810 do j = 1,nuotot 811 ind = 1+(j-1)*2+(iie-1)*2*nuotot 812 aux(1,j) = real(psi(ind),kind=sp) 813 aux(2,j) = real(psi(ind+1),kind=sp) 814 enddo 815 endif 816 else ! non-collinear case, including SOC 817 ! Each 'orbital' coefficient is actually four reals: complex up/down 818 ! (spinor) 819 do j = 1,nuotot 820 ind = 1 + (j-1)*4 + (iie-1)*4*nuotot 821 aux(1,j) = real(psi(ind),kind=sp) 822 aux(2,j) = real(psi(ind+1),kind=sp) 823 aux(3,j) = real(psi(ind+2),kind=sp) 824 aux(4,j) = real(psi(ind+3),kind=sp) 825 enddo 826 endif ! non-coll or not 827 828 endif ! Node == BNode 829 830#ifdef MPI 831C Pass the wf to the master node 832 if (BNode.ne.0) then 833 tag = indwf 834 if (Node .eq. BNode) then 835 call MPI_Send(aux(1,1),size(aux),MPI_real, 836 . 0,tag,MPI_Comm_World,MPIerror) 837 if (debug) print *, Node, " sent msg to 0 with tag ",tag 838 elseif (Node .eq. 0) then 839 call MPI_Recv(aux(1,1),size(aux),MPI_real, 840 . BNode,tag,MPI_Comm_World,MPIStatus,MPIerror) 841 if (debug) print *, "Got msg with tag ",tag, 842 $ " from node ", BNode, ". status: ", mpistatus 843 endif 844 call MPI_Barrier(MPI_Comm_World, MPIError) ! Just to be safe 845 endif 846#endif 847 848C eigenvector is now stored in aux in Node 0, and can be printed 849 850 if (Node .eq. 0) then 851 write(iu) indwf 852 write(iu) eo(indwf)/eV 853 write(iu) (aux(1:,j), j=1,nuotot) 854 endif 855 856 enddo 857 858C Close output file 859 860 if (Node .eq. 0) then 861 call io_close(iu) 862 endif 863 864 ! Restore blocksize value in module 865 Module_blocksize = blocksize_saved 866 867C Deallocate auxiliary arrays 868 869 call memory('D','D',size(aux),'writewave') 870 deallocate(aux) 871 872 end subroutine writew 873 874 end module writewave 875 876