1!=============================================================================== 2! 3! Routines: 4! 5! (1) input() Originally By ? Last Modified 10/5/2009 (gsm) 6! 7! Sets up various data structures. Reads in and distributes wavefunctions. 8! Stores wavefunctions in structures distributed in memory. 9! 10!=============================================================================== 11 12#include "f_defs.h" 13 14module input_m 15 16 use global_m 17 use checkbz_m 18 use createpools_m 19 use eqpcor_m 20 use fftw_m 21 use fullbz_m 22 use input_utils_m 23 use misc_m 24 use read_rho_vxc_m 25 use scissors_m 26 use sort_m 27 use wfn_rho_vxc_io_m 28 use io_utils_m 29 use epsread_hdf5_m 30#ifdef HDF5 31 use hdf5 32#endif 33 use timing_m, only: timing => sigma_timing 34 35 implicit none 36 37 private 38 public :: input 39 40contains 41 42subroutine input(crys,gvec,syms,kp,wpg,sig,wfnk,iunit_c,iunit_k,fnc,fnk,wfnkqmpi,wfnkmpi, wfnk_phonq, wfnk_phonq_mpi, ep_read_in) 43 type (crystal), intent(out) :: crys 44 type (gspace), intent(out) :: gvec 45 type (symmetry), intent(out) :: syms 46 type (kpoints), intent(out) :: kp 47 type (wpgen), intent(out) :: wpg 48 type (siginfo), intent(inout) :: sig ! ZL modifies 49 type (wfnkstates), intent(out) :: wfnk, wfnk_phonq ! ZL adds wfnk_phonq 50 integer, intent(in) :: iunit_c, iunit_k 51 character*20, intent(in) :: fnc, fnk 52 type (wfnkqmpiinfo), intent(out) :: wfnkqmpi 53 type (wfnkmpiinfo), intent(out) :: wfnkmpi, wfnk_phonq_mpi ! ZL adds wfnk_phonq_mpi 54 ! ZL: EP temporary variables 55 type (crystal) :: crys_tmp 56 type (gspace) :: gvec_tmp 57 type (symmetry) :: syms_tmp 58 type (kpoints) :: kp_tmp 59 60 61 character :: fncor*32, tmpfn*16 62 character :: tmpstr*100,tmpstr1*16,tmpstr2*16 63 character :: filenameh5*80 64 integer :: i,ierr,ii,j,jj,k,kk 65 integer :: ikn, irk, nq 66 integer :: ncoul,ncoulb,ncouls,nmtx,nmtx_l,nmtx_col 67 integer, allocatable :: isrt(:), isrti(:) 68 integer :: ipe 69 integer :: Nrod,Nplane,Nfft(3),dNfft(3),dkmax(3),nmpinode 70 integer :: ntband,npools,ndiag_max,noffdiag_max,ntband_max 71 real(DP) :: scalarsize,qk(3),vcell,qtot,omega_plasma,rhog0 72 real(DP) :: mem,fmem,rmem,rmem2,scale,dscale 73 integer, allocatable :: gvecktmp(:,:) 74 logical :: do_ch_sum, skip_checkbz 75 type(grid) :: gr 76 real(DP) :: del 77 integer :: g0(3), i_looper ! ZL: to map EP sig%indk_phonq_g0 78 79 character(len=3) :: sheader 80 integer :: iflavor, freq_dep, error 81 82 integer :: ng_eps, nq_eps, nmtxmax_eps, qgrid_eps(3), nfreq_eps, nfreq_imag_eps 83 real(DP) :: ecuts 84 85 logical :: file_exists 86 87 ! ZL: electron phonon (EP) control 88 logical, intent(in), optional :: ep_read_in 89 logical :: ep_read 90 logical :: do_phonq ! ZL: if do_phonq = .false., wfnk_phonq and wfnk_phonq_mpi will not 91 ! be used, nor allocated or evaluated 92 93 PUSH_SUB(input) 94 ! ZL: set up EP 95 do_phonq = .false. 96 ep_read = .false. 97 98!--------------- 99! Read sigma.inp 100! write(6,*) 'before inread' 101 if (.not.ep_read) then 102 call inread(sig) ! ZL: in inread(), if nphonq .eq. 1, program exits. 103 ! ZL: so after inread(), simply use sig%elph 104 else 105 if (peinf%inode.eq.0) write(6,*) 'Skip inread sigma for dWFN.' 106 endif 107! write(6,*) 'after inread' 108 ! ZL: conditions for generating wfnk_phonq, wfnk_phonq_mpi 109 110!------------------------ 111! Initialize HDF5 112#ifdef HDF5 113 if(sig%use_hdf5 .or. sig%wfn_hdf5) call h5open_f(error) 114#endif 115 116 if(sig%nkn .lt. 1) then 117 if(peinf%inode.eq.0) write(0,*) 'sig%nkn < 1' 118 call die('sig%nkn') 119 endif 120 121!--------------------------------- 122! Determine the available memory 123 124 call procmem(mem,nmpinode) 125 126 if(peinf%inode.eq.0) then 127 write(6,998) mem/1024.0d0**2 128 endif 129998 format(1x,'Memory available: ',f0.1,' MB per PE') 130 131 fmem=mem/8.0d0 132 133!------------------------------- 134! (gsm) Estimate the required memory 135 136 if(peinf%inode.eq.0) then 137 138! determine number of frequencies and number of q-points 139 if (sig%freq_dep.eq.-1) then 140 nq=sig%nq 141 sig%nFreq=1 ! ZL: for EP, this is fine reevaluate 142 endif 143 144 if (sig%freq_dep.eq.0.or.sig%freq_dep.eq.1.or.sig%freq_dep.eq.2.or.sig%freq_dep.eq.3) then 145 146#ifdef HDF5 147 if(sig%use_hdf5) then 148 filenameh5 = 'eps0mat.h5' 149 call read_eps_grid_sizes_hdf5(ng_eps, nq_eps, ecuts, nfreq_eps, & 150 nfreq_imag_eps, nmtxmax_eps, qgrid_eps, freq_dep, filenameh5) 151 write(6,*) 152 write(6,'(1x,a)') 'Read from eps0mat.h5:' 153 write(6,'(1x,a,i0)') '- Number of G-vectors in the charge-density grid: ', ng_eps 154 write(6,'(1x,a,i0)') '- Max. number of G-vectors for dielectric matrices: ', nmtxmax_eps 155 write(6,'(1x,a,f0.3)') '- Cutoff of the dielectric matrix (Ry): ', ecuts 156 write(6,'(1x,a,i0)') '- Number of imaginary frequencies: ', nfreq_imag_eps 157 write(6,'(1x,a,i0)') '- Total number of frequencies: ', nfreq_eps 158 write(6,'(1x,a,i0)') '- Frequency dependence: ', freq_dep 159 write(6,'(1x,a,i0)') '- Number of q-points: ', nq_eps 160 write(6,'(1x,a,3(i0,1x))') '- Q grid: ', qgrid_eps 161 write(6,*) 162 sig%nFreq = nfreq_eps 163 sig%nfreq_imag = nfreq_imag_eps 164 nq = nq_eps 165 166 INQUIRE(FILE="epsmat.h5", EXIST=file_exists) 167 if (file_exists) then 168 filenameh5 = 'epsmat.h5' 169 call read_eps_grid_sizes_hdf5(ng_eps, nq_eps, ecuts, nfreq_eps, nfreq_imag_eps, & 170 nmtxmax_eps, qgrid_eps, freq_dep, filenameh5) 171 write(6,'(1x,a)') 'Read from epsmat.h5:' 172 write(6,'(1x,a,i0)') '- Number of G-vectors in the charge-density grid: ', ng_eps 173 write(6,'(1x,a,i0)') '- Max. number of G-vectors for dielectric matrices: ', nmtxmax_eps 174 write(6,'(1x,a,f0.3)') '- Cutoff of the dielectric matrix: ', ecuts 175 write(6,'(1x,a,i0)') '- Number of imaginary frequencies: ', nfreq_imag_eps 176 write(6,'(1x,a,i0)') '- Total number of frequencies: ', nfreq_eps 177 write(6,'(1x,a,i0)') '- Frequency dependence: ', freq_dep 178 write(6,'(1x,a,i0)') '- Number of q-points: ', nq_eps 179 write(6,'(1x,a,3(i0,1x))') '- Q grid: ', qgrid_eps 180 write(6,*) 181 nq = nq + nq_eps 182 endif 183 else 184#endif 185 186 call open_file(10,file='eps0mat',form='unformatted',status='old') 187 read(10) 188 read(10) i, sig%nFreq 189 read(10) 190 read(10) 191 read(10) 192 read(10) 193 read(10) ecuts 194 read(10) nq 195 call close_file(10) 196 call open_file(11,file='epsmat',form='unformatted',status='old',iostat=ierr) 197 if (ierr.eq.0) then 198 read(11) 199 read(11) 200 read(11) 201 read(11) 202 read(11) 203 read(11) 204 read(11) 205 read(11) nq_eps 206 nq = nq + nq_eps 207 call close_file(11) 208 endif 209 210#ifdef HDF5 211 endif 212#endif 213 214 endif 215 216 endif 217 218! FHJ: Read header of WFN file (and WFNq, if requred), perform consistency 219! checks, and figure out the number of bands and ncrit. 220 sheader = 'WFN' 221 iflavor = 0 222 if(sig%wfn_hdf5) then 223 else 224 if(peinf%inode == 0) then 225 if(.not.ep_read) then 226 write(6,'(a)') ' Reading header of WFN_inner' 227 call open_file(25,file='WFN_inner',form='unformatted',status='old') 228 else 229 ! ZL: open dWFN 230 write(6,'(a)') ' Reading header of dWFN' 231 call open_file(25,file='dWFN',form='unformatted',status='old') 232 endif 233 endif 234 call read_binary_header_type(25, sheader, iflavor, kp, gvec, syms, crys) 235 endif 236 call check_trunc_kpts(sig%icutv, kp) ! This does not affect EP 237 238 if (sig%ecutb<TOL_ZERO) then 239 sig%ecutb = kp%ecutwfc 240 endif 241 if (sig%freq_dep/=-1 .and. sig%ecuts<TOL_ZERO) then 242#ifdef MPI 243 call MPI_Bcast(ecuts, 1, MPI_REAL_DP, 0, MPI_COMM_WORLD, mpierr) 244#endif 245 sig%ecuts = ecuts 246 endif 247 248 if (sig%ntband==0 .and. (.not.ep_read)) sig%ntband = kp%mnband - 1 249 250! DVF: Here we re-define the number of bands in the dynamical self energy 251! summation to exclude the deep core states. This ensures that the pools will be 252! set up properly. This achieves most of what is needed to exclude core states. 253! There is a little more work in the i/o routines to make sure you`re getting the 254! higher valence states. 255! ZL: we disable this operation if doing elph to avoid possible mismaps of band indices 256 if ((sig%ncore_excl.ne.0) .and. sig%elph) then 257 call die('Electron-Phonon calculation is not consistent with sig%ncore_excl.ne.0 .') 258 endif 259 sig%ntband=sig%ntband-sig%ncore_excl 260 261 SAFE_ALLOCATE(kp%elda, (sig%ntband,kp%nrk,kp%nspin)) 262 kp%el(:,:,:) = kp%el(:,:,:) - sig%avgpot / ryd 263! DVF : add sig%ncore_excl to make sure we are getting the right energies. kp%el is setup in 264! wfn_rho_vxc_io.f90 routine (see .p.f for something comprehensible) and includes core states. 265 kp%elda(1:sig%ntband, 1:kp%nrk, 1:kp%nspin)=ryd*kp%el(1+sig%ncore_excl:sig%ntband+sig%ncore_excl, 1:kp%nrk, 1:kp%nspin) 266 call scissors_shift(kp, sig%scis, sig%spl_tck) 267 ! If QP corrections requested, read the corrected QP energies from file (in eV) 268 if(sig%eqp_corrections .and. (.not.ep_read)) then 269 fncor='eqp.dat' 270! DVF : add sig%ncore_excl to make sure we are getting the right energies. The eqp.dat file 271! will have band indices referenced to the total number of states, including the core states. 272 call eqpcor(fncor,peinf%inode,peinf%npes,kp,1+sig%ncore_excl,sig%ntband+sig%ncore_excl,0,0,kp%el,kp%el,kp%el,1,0) 273 ! note: want in Ry since conversion occurs below 274 endif 275 if (peinf%inode==0 .and. (.not.ep_read)) then 276 call find_efermi(sig%rfermi, sig%efermi, sig%efermi_input, kp, sig%ntband+sig%ncore_excl, 1+sig%ncore_excl, & 277 "unshifted grid", should_search=.true., should_update=.true., write7=.false.) 278 endif 279 280 if (sig%nvband==0 .and. (.not.ep_read)) then 281 sig%nvband = minval(kp%ifmax) 282 sig%ncrit = maxval(kp%ifmax) - minval(kp%ifmax) 283 endif 284 285! DVF: Here we re-define the number of bands in the bare exchange 286! summation to exclude the deep core states. 287! ZL: this has been disabled already for EP 288 sig%nvband=sig%nvband-sig%ncore_excl 289 290#ifdef MPI 291 call MPI_Bcast(sig%efermi, 1, MPI_REAL_DP, 0, MPI_COMM_WORLD, mpierr) 292 call MPI_Bcast(sig%nvband, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr) 293 call MPI_Bcast(sig%ncrit, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr) 294#endif 295 296! DVF: again, this is to exclude the core states. 297 if ((sig%freq_dep==0 .and. sig%exact_ch==1) .or. sig%freq_dep==-1) then 298 if (sig%ntband + sig%ncore_excl > max(maxval(sig%diag(1:sig%ndiag)), sig%nvband + sig%ncore_excl + sig%ncrit)) then 299 ! Do not waste time and memory with the unoccupied bands which will not be used 300 sig%ntband = max(maxval(sig%diag(1:sig%ndiag))-sig%ncore_excl, sig%nvband + sig%ncrit) 301 if (peinf%inode==0 .and. sig%freq_dep==0) then 302 write(0,'(a)') "WARNING: Resetting number_bands to number of bands actually needed for exact static CH." 303 else if (peinf%inode==0) then 304 write(0,'(a)') "WARNING: Resetting number_bands to number of bands actually needed for exchange." 305 endif 306 endif 307 endif 308 309 if(peinf%inode == 0) then 310 if(.not.ep_read) then 311 call logit('Reading WFN_inner -- gvec info') 312 else 313 call logit('Reading dWFN -- gvec info') 314 endif 315 endif 316 SAFE_ALLOCATE(gvec%components, (3, gvec%ng)) 317 if(sig%wfn_hdf5) then 318 else 319 call read_binary_gvectors(25, gvec%ng, gvec%ng, gvec%components) 320 endif 321 322!----------------------- 323! sort G-vectors with respect to their kinetic energy 324 SAFE_ALLOCATE(isrt, (gvec%ng)) 325 SAFE_ALLOCATE(isrti, (gvec%ng)) 326 SAFE_ALLOCATE(gvec%ekin, (gvec%ng)) 327 call kinetic_energies(gvec, crys%bdot, gvec%ekin) 328 call sortrx(gvec%ng, gvec%ekin, isrt, gvec = gvec%components) 329 ncouls = gcutoff(gvec%ng, gvec%ekin, isrt, sig%ecuts) 330 ncoulb = gcutoff(gvec%ng, gvec%ekin, isrt, sig%ecutb) 331 SAFE_DEALLOCATE_P(gvec%ekin) 332 333 if (peinf%inode==0) then 334 write(6,*) 335 write(6,'(a)') ' Calculation parameters:' 336 write(6,'(1x,a,f0.2)') '- Cutoff of the bare Coulomb interaction (Ry): ', sig%ecutb 337 write(6,'(1x,a,f0.2)') '- Cutoff of the screened Coulomb interaction (Ry): ', sig%ecuts 338 write(6,'(1x,a,i0)') '- Number of G-vectors up to the bare int. cutoff: ', ncoulb 339 write(6,'(1x,a,i0)') '- Number of G-vectors up to the screened int. cutoff: ', ncouls 340 write(6,'(1x,a,i0)') '- Total number of bands in the calculation: ', sig%ntband 341 write(6,'(1x,a,i0)') '- Number of fully occupied valence bands: ', sig%nvband 342 write(6,'(1x,a,i0)') '- Number of partially occ. conduction bands: ', sig%ncrit 343 write(6,*) 344 endif 345 346 if(sig%ecuts > sig%ecutb) then 347 call die("The screened_coulomb_cutoff cannot be greater than the bare_coulomb_cutoff.", & 348 only_root_writes = .true.) 349 endif 350 if(sig%ntband > kp%mnband) then 351 if(.not.ep_read) then 352 call die("number_bands is larger than are available in WFN_inner", only_root_writes = .true.) 353 else 354 call die("number_bands is larger than are available in dWFN", only_root_writes = .true.) 355 endif 356 endif 357 do ii = 1, sig%ndiag 358 if (sig%diag(ii)>sig%ntband) then 359 write(0,'(a,3(i0,1x))') 'ERROR: the band_index_max cannot be greater than number_bands: ', & 360 ii, sig%diag(ii), sig%ntband 361 call die('band_index_max cannot be greater than number_bands', only_root_writes = .true.) 362 endif 363 enddo 364 365! FHJ: do vanilla WFN reading stuff 366 367 SAFE_ALLOCATE(gvecktmp, (3,gvec%ng)) 368 gvecktmp(:,:)=gvec%components(:,:) 369 do i=1,gvec%ng 370 gvec%components(:,i)=gvecktmp(:,isrt(i)) 371 isrti(isrt(i)) = i 372 enddo 373 SAFE_DEALLOCATE(gvecktmp) 374 375 call gvec_index(gvec) 376 377!------------------------------------------------------------------------------ 378! Pools, distribution and memory estimation 379!------------------------------------------------------------------------------ 380 381 if(peinf%inode == 0) then 382 ncoul=max(ncouls,ncoulb) 383 nmtx=ncouls 384 385! divide bands over processors (this is repeated below) 386 ntband=min(sig%ntband,kp%mnband) 387 if (peinf%npools .le. 0 .or. peinf%npools .gt. peinf%npes) then 388 call createpools(sig%ndiag,ntband,peinf%npes,npools,ndiag_max,ntband_max) 389 else 390 npools = peinf%npools 391 if (mod(sig%ndiag,npools).eq.0) then 392 ndiag_max=sig%ndiag/npools 393 else 394 ndiag_max=sig%ndiag/npools+1 395 endif 396 if (mod(ntband,peinf%npes/npools).eq.0) then 397 ntband_max=ntband/(peinf%npes/npools) 398 else 399 ntband_max=ntband/(peinf%npes/npools)+1 400 endif 401 endif 402 if (sig%noffdiag.gt.0) then 403 if (mod(sig%noffdiag,npools).eq.0) then 404 noffdiag_max=sig%noffdiag/npools 405 else 406 noffdiag_max=sig%noffdiag/npools+1 407 endif 408 endif 409 410 nmtx_l=int(sqrt(dble(nmtx)**2/dble(peinf%npes/npools))) 411 nmtx_col=int(dble(nmtx)/dble(peinf%npes/npools))+1 412 413 scalarsize = sizeof_scalar() 414 415! required memory 416 rmem=0.0d0 417! arrays eps and epstemp in program main 418 if (sig%freq_dep.eq.0.or.sig%freq_dep.eq.1) then 419 rmem=rmem+(dble(nmtx_l)**2+dble(nmtx))*scalarsize 420 rmem=rmem+(dble(nmtx_col)*nmtx)*scalarsize 421 endif 422! wtilde_array, imagloop_ig, imaginary_flag in mtxel_cor 423 if (sig%freq_dep.eq.1) then 424 rmem=rmem+(dble(nmtx_col)*nmtx)*16 425 rmem=rmem+(dble(nmtx_col)*nmtx)*4 426 rmem=rmem+(dble(nmtx_col)*nmtx)*1 427 endif 428! arrays epsR, epsRtemp, epsA, and epsAtemp in program main 429 if (sig%freq_dep.eq.2 .or.sig%freq_dep.eq.3) then 430 rmem=rmem+(dble(nmtx_l)**2+dble(nmtx)) & 431 *dble(sig%nFreq)*2.0d0*scalarsize 432 endif 433! array epsmpi%eps in subroutine epscopy 434 if (sig%freq_dep.eq.0.or.sig%freq_dep.eq.1) then 435 rmem=rmem+(dble(nmtx_col)*nmtx)*dble(nq)*scalarsize 436 endif 437! arrays epsmpi%epsR and epsmpi%epsA in subroutine epscopy 438 if (sig%freq_dep.eq.2 .or.sig%freq_dep.eq.3) then 439 rmem=rmem+(nmtx_col*nmtx)*dble(nq) & 440 *dble(sig%nFreq)*2.0d0*scalarsize 441 endif 442! array aqs in program main 443 rmem=rmem+dble(ntband_max)*dble(ncoul)*scalarsize 444! array aqsaug in program main 445 if (sig%noffdiag.gt.0) then 446 rmem=rmem+dble(ntband_max)*dble(ncoul) & 447 *dble(sig%ndiag)*dble(sig%nspin)*scalarsize 448 endif 449! array aqsch in program main 450 if (sig%exact_ch.eq.1) then 451 rmem=rmem+dble(ncoul)*scalarsize 452 endif 453! arrays aqsaugchd and aqsaugcho in program main 454 if (sig%exact_ch.eq.1.and.nq.gt.1) then 455 rmem=rmem+dble(ncoul)*dble(ndiag_max) & 456 *dble(sig%nspin)*scalarsize 457 if (sig%noffdiag.gt.0) then 458 rmem=rmem+dble(ncoul)*dble(noffdiag_max)* & 459 dble(sig%nspin)*scalarsize 460 endif 461 endif 462! array wfnk%zk in subroutine input 463 rmem=rmem+dble(ndiag_max)*dble(kp%ngkmax)* & 464 dble(sig%nspin*kp%nspinor)*scalarsize 465! array wfnkq%zkq in subroutine genwf 466 rmem=rmem+dble(ntband_max)*dble(kp%ngkmax) & 467 *dble(sig%nspin*kp%nspinor)*scalarsize 468! arrays zc in subroutine input and zin in subroutine genwf 469 rmem=rmem+dble(kp%ngkmax)*dble(kp%nspin*kp%nspinor)*scalarsize 470! array wfnkoff%zk in program main and subroutines mtxel_vxc and mtxel_cor 471 if ((sig%exact_ch.eq.1.and.sig%noffdiag.gt.0).or. & 472 (.not.sig%use_vxcdat .and. .not.sig%sigma_correction .and. .not. sig%is_EXX)) then 473 rmem=rmem+2.0d0*dble(kp%ngkmax)*scalarsize 474 endif 475! array gvec%index_vec in input 476 rmem=rmem+dble(gvec%nFFTgridpts)*4.0d0 477! arrays fftbox1 and fftbox2 in subroutines mtxel and mtxel_occ 478 call setup_FFT_sizes(gvec%FFTgrid,Nfft,scale) 479 rmem=rmem+dble(Nfft(1)*Nfft(2)*Nfft(3))*32.0d0 480! arrays wfnkqmpi%isort and wfnkqmpi%cg in subroutine input 481 rmem=rmem+dble(kp%ngkmax)*dble(kp%nrk)*4.0d0+ & 482 dble(kp%ngkmax)*dble(ntband_max)*dble(sig%nspin*kp%nspinor)* & 483 dble(kp%nrk)*dble(scalarsize) 484! arrays wfnkmpi%isort and wfnkmpi%cg in subroutine input 485 rmem=rmem+dble(kp%ngkmax)*dble(sig%nkn)*4.0d0+ & 486 dble((ndiag_max*kp%ngkmax)/(peinf%npes/npools))* & 487 dble(sig%nspin*kp%nspinor)*dble(sig%nkn)* & 488 dble(scalarsize) 489 490989 format(1x,'Memory required for execution: ',f0.1,' MB per PE') 491 write(6,989) rmem/1024.0d0**2 492 493! random numbers 494 rmem=0.0D0 495 if (sig%icutv/=TRUNC_BOX) then 496! arrays ran, qran, and qran2 497! (ran is deallocated before qran2 is allocated) 498 rmem=rmem+6.0D0*dble(nmc)*8.0D0 499 endif 500! various truncation schemes 501 rmem2=0.0d0 502! cell wire truncation 503 if (sig%icutv==TRUNC_WIRE) then 504 dkmax(1) = gvec%FFTgrid(1) * n_in_wire 505 dkmax(2) = gvec%FFTgrid(2) * n_in_wire 506 dkmax(3) = 1 507 call setup_FFT_sizes(dkmax,dNfft,dscale) 508! array fftbox_2D 509 rmem2=rmem2+dble(dNfft(1))*dble(dNfft(2))*16.0d0 510! array inv_indx 511 rmem2=rmem2+dble(Nfft(1))*dble(Nfft(2))*dble(Nfft(3))* & 512 4.0d0 513! array qran 514 rmem2=rmem2+3.0D0*dble(nmc)*8.0D0 515 endif 516! cell box truncation (parallel version only) 517 if (sig%icutv==TRUNC_BOX) then 518 dkmax(1) = gvec%FFTgrid(1) * n_in_box 519 dkmax(2) = gvec%FFTgrid(2) * n_in_box 520 dkmax(3) = gvec%FFTgrid(3) * n_in_box 521 call setup_FFT_sizes(dkmax,dNfft,dscale) 522 if (mod(dNfft(3),peinf%npes) == 0) then 523 Nplane = dNfft(3)/peinf%npes 524 else 525 Nplane = dNfft(3)/peinf%npes+1 526 endif 527 if (mod(dNfft(1)*dNfft(2),peinf%npes) == 0) then 528 Nrod = (dNfft(1)*dNfft(2))/peinf%npes 529 else 530 Nrod = (dNfft(1)*dNfft(2))/peinf%npes+1 531 endif 532! array fftbox_2D 533 rmem2=rmem2+dble(dNfft(1))*dble(dNfft(2))*dble(Nplane)* & 534 16.0d0 535! array fftbox_1D 536 rmem2=rmem2+dble(dNfft(3))*dble(Nrod)*16.0d0 537! array dummy 538! rmem2=rmem2+dble(dNfft(1))*dble(dNfft(2))*16.0d0 539! arrays dummy1 and dummy2 540 rmem2=rmem2+dble(Nrod)*dble(peinf%npes+1)*16.0d0 541! array inv_indx 542 rmem2=rmem2+dble(Nfft(1))*dble(Nfft(2))*dble(Nfft(3))* & 543 4.0d0 544 endif 545 if (rmem2 .gt. rmem) rmem = rmem2 546988 format(1x,'Memory required for vcoul: ',f0.1,' MB per PE') 547 write(6,988) rmem/1024.0d0**2 548 write(6,*) 549 endif 550 551 if(peinf%inode == 0) then 552!---------------------------------------------- 553! Compute cell volume from wave number metric 554 call get_volume(vcell,crys%bdot) 555 if (abs(crys%celvol-vcell).gt.TOL_Small) then 556 call die('volume mismatch') 557 endif 558 559! (gsm) check consistency of spin indices 560 561 do k=1,sig%nspin 562 if (sig%spin_index(k).lt.1.or.sig%spin_index(k).gt.kp%nspin) & 563 call die('inconsistent spin indices') 564 enddo 565 566 if(sig%ntband.gt.kp%mnband) then 567 write(tmpstr1,660) sig%ntband 568 write(tmpstr2,660) kp%mnband 569 write(0,666) TRUNC(tmpstr1), TRUNC(tmpstr2) 570 if(.not.ep_read) then 571 call die('More bands specified in sigma.inp than available in WFN_inner.') 572 else 573 call die('More bands specified in sigma.inp than available in dWFN.') 574 endif 575660 format(i16) 576666 format(1x,'The total number of bands (',a,') specified in sigma.inp',/, & 577 3x,'is larger than the number of bands (',a,') available in WFN_inner.',/) 578 endif 579 do_ch_sum = .not.((sig%freq_dep == -1) .or. (sig%freq_dep == 0 .and. sig%exact_ch == 1)) 580 if(sig%ntband.eq.kp%mnband) then 581 if(.not.ep_read) then 582 call die("You must provide one more band in WFN_inner than used in sigma.inp number_bands in order to assess degeneracy.") 583 else 584 call die("You must provide one more band in dWFN than used in sigma.inp number_bands in order to assess degeneracy.") 585 endif 586 endif 587 588!-------------------------------------------------------------------- 589! SIB: Find the k-points in sig%kpt in the list kp%rk (die if not found). 590! sig%indkn holds the index of the k-points in sig%kpt in kp%rk, i.e. 591! kp%rk(sig%indkn(ikn))=sig%kpt(ikn) 592 593 SAFE_ALLOCATE(sig%indkn, (sig%nkn)) 594 do ikn=1,sig%nkn 595 sig%indkn(ikn)=0 596 qk(:)=sig%kpt(:,ikn) 597 do irk=1,kp%nrk 598 if(all(abs(kp%rk(1:3,irk)-qk(1:3)) .lt. TOL_Small)) sig%indkn(ikn)=irk 599 enddo 600 601 ! ZL: print information 602 if(do_phonq .and. peinf%inode.eq.0) write(6,*) "sig%indkn(",ikn,") =", sig%indkn(ikn) 603 604 if(sig%indkn(ikn) .eq. 0) then 605 if(.not.ep_read) then 606 write(0,'(a,3f10.5,a)') 'Could not find the k-point ', (qk(i),i=1,3), ' among those read from WFN_inner :' 607 else 608 write(0,'(a,3f10.5,a)') 'Could not find the k-point ', (qk(i),i=1,3), ' among those read from dWFN :' 609 endif 610 write(0,'(3f10.5)') ((kp%rk(i,irk),i=1,3),irk=1,kp%nrk) 611 call die('k-point in sigma.inp k_points block not available.') 612 endif 613 enddo 614 615 ! ZL: EP, when reading regular WFN_inner, NOT dWFN, we map k+phonq points 616 if(do_phonq .and. (.not. ep_read) .and. (sig%nphonq .eq. 1)) then ! ZL: only nphonq=1 supported for now 617 SAFE_ALLOCATE(sig%indk_phonq, (sig%nkn * sig%nphonq)) 618 SAFE_ALLOCATE(sig%indk_phonq_g0, (3, sig%nkn * sig%nphonq)) 619 do ikn=1,sig%nkn*sig%nphonq 620 sig%indk_phonq(ikn)=0 621 qk(:)=sig%k_phonq(:,ikn) 622 irk_loop_g0: do irk=1,kp%nrk 623 do i_looper=1,3 624 ! ZL: this idea is taken from genwf_mpi.f90 625 ! the most interesting part is that 626 ! it uses the fact that g0(3) is INTEGER 627 ! therefore it will automatically drop the fractional part 628 ! and keep only the integer part 629 ! This is also guaranteed by adding TOL_small so that 0.999999 630 ! is indeed 1 631 ! Then we just need to compare if del is close to integer 632 ! NOTE: gmap uses a different idea, it maps between old and new wfn 633 ! but here we map to master G-list 634 del = qk(i_looper) - kp%rk(i_looper,irk) 635 if (del .ge. 0.0d0) g0(i_looper) = del + TOL_Small 636 if (del .lt. 0.0d0) g0(i_looper) = del - TOL_Small 637 if (abs(del-g0(i_looper)) .gt. TOL_Small) cycle irk_loop_g0 638 enddo 639 sig%indk_phonq(ikn)=irk 640 sig%indk_phonq_g0(:,ikn)=g0(:) 641 exit irk_loop_g0 642 enddo irk_loop_g0 643 644 ! ZL: print information 645 if(peinf%inode .eq. 0) then 646 write(6,*) "sig%indk_phonq(",ikn,") =", sig%indk_phonq(ikn) 647 write(6,*) "sig%indk_phonq_g0(",ikn,") =", sig%indk_phonq_g0(:,ikn) 648 endif 649 650 if(sig%indk_phonq(ikn) .eq. 0) then 651 write(0,'(a,3f10.5,a)') 'Could not find the k_phonq-point ', (qk(i),i=1,3), ' among those in WFN_inner :' 652 write(0,'(3f10.5)') ((kp%rk(i,irk),i=1,3),irk=1,kp%nrk) 653 call die('k_phonq-point from sigma.inp k_points+phonq not available.') 654 endif 655 enddo 656 endif ! EP k+phonq mapping 657 658 659 if(do_ch_sum .and. (.not.ep_read)) then 660 if(any(abs(kp%el(sig%ntband+sig%ncore_excl, 1:kp%nrk, 1:kp%nspin) - & 661 kp%el(sig%ntband+sig%ncore_excl+1, 1:kp%nrk, 1:kp%nspin)) .lt. sig%tol/RYD)) then 662 if(sig%degeneracy_check_override) then 663 write(0,'(a)') & 664 "WARNING: Selected number of bands for CH sum (number_bands) breaks degenerate subspace. " // & 665 "Run degeneracy_check.x for allowable numbers." 666 write(0,*) 667 else 668 write(0,'(a)') & 669 "Run degeneracy_check.x for allowable numbers, or use keyword " // & 670 "degeneracy_check_override to run anyway (at your peril!)." 671 call die("Selected number of bands for CH sum (number_bands) breaks degenerate subspace.") 672 endif 673 endif 674 endif 675 endif 676 677!---------------------------------------------------------------- 678! Check consistency 679 680 if (kp%mnband < maxval(sig%diag) .and. peinf%inode == 0) then 681 write(0,*) 'The highest requested band is ', maxval(sig%diag),' but WFN_inner contains only ', kp%mnband,' bands.' 682 call die('too few bands') 683 endif 684 685 ! qgridsym has no meaning when there is only one q-point (Gamma) 686 ! setting this to false means warnings about its applicability will not be triggered 687 if(all(abs(kp%rk(1:3,1:kp%nrk)) < TOL_Zero)) sig%qgridsym = .false. 688 689 if(peinf%inode == 0) then 690 ! ZL: assess_degeneracies will allocate kp%degeneracy 691 ! kp%el is still in Ry here, but sig%tol is in eV. 692 call assess_degeneracies(kp, kp%el(sig%ntband+sig%ncore_excl+1, :, :), sig%ntband, sig%efermi, sig%tol/RYD, sig = sig, & 693 ncore_excl=sig%ncore_excl) 694 if(.not.ep_read) then 695 call calc_qtot(kp, crys%celvol, sig%efermi, qtot, omega_plasma, write7=.false.) 696 endif 697 endif 698 699 if ( peinf%inode == 0 ) call timing%start(timing%fullbz) 700 gr%nr = kp%nrk 701 SAFE_ALLOCATE(gr%r, (3, gr%nr)) 702 gr%r = kp%rk 703 call fullbz(crys,syms,gr,syms%ntran,skip_checkbz,wigner_seitz=.false.,paranoid=.true.) 704 if(.not.ep_read) then 705 tmpfn='WFN_inner' 706 else 707 tmpfn='dWFN' 708 endif 709 if (.not. skip_checkbz) then 710 call checkbz(gr%nf,gr%f,kp%kgrid,kp%shift,crys%bdot, & 711 tmpfn,'k',.false.,sig%freplacebz,sig%fwritebz) 712 endif 713 call dealloc_grid(gr) 714 if ( peinf%inode == 0 ) call timing%stop(timing%fullbz) 715 716 ! For discussion of how q-symmetry may (and may not) be used with degenerate states, 717 ! see Hybertsen and Louie, Phys. Rev. B 35, 5585 (1987) Appendix A 718 if(sig%qgridsym .and. sig%noffdiag > 0) then 719 if(peinf%inode == 0) then 720 write(0,'(a)') "WARNING: Cannot calculate offdiagonal elements unless no_symmetries_q_grid is set." 721 write(0,'(a)') "This flag is being reset to enable the calculation." 722 endif 723 sig%qgridsym = .false. 724 endif 725 726 ! ZL: for EP, turn off q-symmetry. Interestingly, this coincides with off-diag 727 if(sig%qgridsym .and. ep_read) then 728 if(peinf%inode == 0) then 729 write(0,'(a)') "WARNING: Cannot calculate electron phonon unless no_symmetries_q_grid is set." 730 write(0,'(a)') "This flag is being reset to enable the calculation." 731 endif 732 sig%qgridsym = .false. 733 endif 734 735 if(peinf%inode == 0) then 736 if(sig%qgridsym) then 737 write(6,'(1x,a)') 'Q-grid symmetries are being used.' 738 else 739 write(6,'(1x,a)') 'Q-grid symmetries are not being used.' 740 endif 741 endif 742 743 kp%el(:,:,:)=kp%el(:,:,:)*ryd 744 745 !---------------------------------------------------------------- 746 ! Distribute data 747 748#ifdef MPI 749 if(sig%qgridsym .or. sig%noffdiag > 0 .or. sig%ncrit > 0) then 750 if(peinf%inode.ne.0) then 751 SAFE_ALLOCATE(kp%degeneracy, (sig%ntband, kp%nrk, kp%nspin)) 752 endif 753 call MPI_Bcast(kp%degeneracy(1,1,1), sig%ntband * kp%nrk * kp%nspin, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr) 754 endif 755 756 if(peinf%inode.ne.0) then 757 SAFE_ALLOCATE(sig%indkn, (sig%nkn)) 758 if(do_phonq) then 759 SAFE_ALLOCATE(sig%indk_phonq, (sig%nkn*sig%nphonq)) 760 SAFE_ALLOCATE(sig%indk_phonq_g0, (3, sig%nkn*sig%nphonq)) 761 endif 762 endif 763 call MPI_Bcast(sig%indkn, sig%nkn, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr) 764 if(do_phonq) then 765 call MPI_Bcast(sig%indk_phonq, sig%nkn*sig%nphonq, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr) 766 call MPI_Bcast(sig%indk_phonq_g0, 3*sig%nkn*sig%nphonq, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr) 767 endif 768#endif 769 770!---------------------------------------------------------------- 771! CHP: Set the init_freq_eval relative to the Fermi energy in case 772! of full frequency. This part should NOT be put before the 773! Fermi energy is set. 774 775 if ((sig%freq_dep==2 .or. (sig%freq_dep==1.and.sig%fdf==-3)) .and. sig%freq_grid_shift==1) then 776 ! FHJ: only do this for FE-based freq. grids 777 sig%freqevalmin = sig%freqevalmin + sig%efermi 778 endif 779 780!-------------------------------------------------------------------- 781! Read in the exchange-correlation potential and store in array vxc 782 783 if(.not.sig%use_vxcdat .and. .not.sig%sigma_correction .and. .not. sig%is_EXX) then 784 785 call read_vxc(sig, gvec, kp, syms, crys, isrti, isrt, 1) 786 787 endif ! not using vxc.dat 788 789 if(.not. sig%use_vxc2dat) then 790 ! This is for hybrid functional like calculation (one shot) 791 if (sig%freq_dep.eq.-1 .and. ((1.0d0 - sig%coulomb_mod%long_range_frac_fock > TOL_SMALL) .or. & 792 (1.0d0 - sig%coulomb_mod%short_range_frac_fock > TOL_SMALL))) & 793 call read_vxc(sig, gvec, kp, syms, crys, isrti, isrt, 2) 794 795 endif ! not using vxc2.dat 796 797 ! ZL: read dVXC when reading dWFN 798 if (sig%elph .and. ep_read) then 799 call read_vxc(sig, gvec, kp, syms, crys, isrti, isrt, 3) 800 endif 801 802!-------------------------------------------------------------------- 803! Read in the charge density and store in array rho (formerly known as CD95) 804! CD95 Ref: http://www.nature.com/nature/journal/v471/n7337/full/nature09897.html 805 806 if(sig%freq_dep.eq.1) then 807 call timing%start(timing%input_read) 808 809 call read_rho(wpg, gvec, kp, syms, crys, isrti, isrt, 'WFN_inner') ! ZL: only read from WFN_inner, NOT dWFN, for RHO 810 811 rhog0 = sum(wpg%nelec(1:kp%nspin)) 812 if(peinf%inode == 0 .and. abs(rhog0 - qtot) > TOL_Small) then 813 write(0,'(a,g12.6,a,g12.6)') 'WARNING: RHO has total charge ', rhog0, ' but WFN_inner occupations give charge ', qtot 814 endif 815 816 call timing%stop(timing%input_read) 817 endif ! sig%freq_dep 818 819!------------------------------------------------------------------------ 820! Divide bands over processors 821! 822! sig%ntband number of total (valence and conduction) bands 823! peinf%npools number of parallel sigma calculations 824! peinf%ndiag_max maximum number of diag sigma calculations 825! peinf%noffdiag_max maximum number of offdiag sigma calculations 826! peinf%ntband_max maximum number of total bands per node 827! 828! ntband_node number of total bands per current node 829! nvband_node number of valence bands per current node 830! peinf%indext(itb) indices of total bands belonging to current node 831! 832! peinf%index_diag band index for diag sigma calculation 833! peinf%flag_diag flag for storing diag sigma calculation 834! peinf%index_offdiag band index for offdiag sigma calculation 835! peinf%flag_offdiag flag for storing offdiag sigma calculation 836! 837! flags are needed in case of uneven distribution, nodes still do the 838! calculation because they share epsilon and wavefunctions and need to 839! participate in global communications, but the result is not stored 840 841 if (peinf%npools .le. 0 .or. peinf%npools .gt. peinf%npes) then 842 call createpools(sig%ndiag,sig%ntband,peinf%npes,npools,ndiag_max,ntband_max) 843 peinf%npools = npools 844 peinf%ndiag_max = ndiag_max 845 peinf%ntband_max = ntband_max 846 else 847 if (mod(sig%ndiag,peinf%npools).eq.0) then 848 peinf%ndiag_max=sig%ndiag/peinf%npools 849 else 850 peinf%ndiag_max=sig%ndiag/peinf%npools+1 851 endif 852 853 if (mod(sig%ntband,peinf%npes/peinf%npools).eq.0) then 854 peinf%ntband_max=sig%ntband/(peinf%npes/peinf%npools) 855 else 856 peinf%ntband_max=sig%ntband/(peinf%npes/peinf%npools)+1 857 endif 858 endif 859 860! JRD Create Separate Comms for Each Pool 861 862 peinf%npes_pool = peinf%npes/peinf%npools 863 peinf%my_pool=peinf%inode/peinf%npes_pool 864 peinf%pool_rank = mod(peinf%inode,peinf%npes_pool) 865 866#ifdef MPI 867! call MPI_Comm_Group(MPI_COMM_WORLD,orig_group,mpierr) 868! 869! SAFE_ALLOCATE(ranks,(peinf%npes)) 870! do ii=1,peinf%npes 871! ranks(ii)=ii-1 872! enddo 873! 874! if (peinf%my_pool .ne. -1) then 875! call MPI_Group_incl(orig_group, peinf%npes_pool, & 876! ranks(peinf%my_pool*peinf%npes_pool+1:(peinf%my_pool+1)*peinf%npes_pool), peinf%pool_group,mpierr) 877! else 878!! We need a group to send comm_create below. If you are not in a legit pool, we just use the first group 879! call MPI_Group_incl(orig_group, peinf%npes_pool, & 880! ranks(0*peinf%npes_pool+1:(0+1)*peinf%npes_pool), peinf%pool_group,mpierr) 881! endif 882! 883!! All procs in original comm need to do this 884! 885! call MPI_Comm_create(MPI_COMM_WORLD,peinf%pool_group,peinf%pool_comm,mpierr) 886 887 call MPI_Comm_split(MPI_COMM_WORLD, peinf%my_pool, peinf%pool_rank, peinf%pool_comm, mpierr) 888 889 if (peinf%my_pool .gt. (peinf%npools-1)) then 890 peinf%my_pool = -1 891 endif 892 893! !write(6,*) peinf%inode,"Created groups comms" 894 895! SAFE_DEALLOCATE(ranks) 896#endif 897 898 if (sig%noffdiag.gt.0) then 899 if (mod(sig%noffdiag,peinf%npools).eq.0) then 900 peinf%noffdiag_max=sig%noffdiag/peinf%npools 901 else 902 peinf%noffdiag_max=sig%noffdiag/peinf%npools+1 903 endif 904 endif 905 906 SAFE_ALLOCATE(peinf%index_diag, (peinf%ndiag_max)) 907 SAFE_ALLOCATE(peinf%flag_diag, (peinf%ndiag_max)) 908 peinf%index_diag=1 909 peinf%flag_diag=.false. 910 do ii=1,peinf%ndiag_max*peinf%npools 911 jj=(ii-1)/peinf%npools+1 ! which number in the pool it is 912 kk=mod(ii-1,peinf%npools) ! which pool it is in 913 if (peinf%inode/(peinf%npes/peinf%npools).eq.kk) then 914 if (ii.le.sig%ndiag) then 915 peinf%index_diag(jj)=ii 916 peinf%flag_diag(jj)=.true. 917 else 918 peinf%index_diag(jj)=1 919 peinf%flag_diag(jj)=.false. 920 endif 921 endif 922 enddo 923 924 if (sig%noffdiag.gt.0) then 925 SAFE_ALLOCATE(peinf%index_offdiag, (peinf%noffdiag_max)) 926 SAFE_ALLOCATE(peinf%flag_offdiag, (peinf%noffdiag_max)) 927 peinf%index_offdiag=1 928 peinf%flag_offdiag=.false. 929 do ii=1,peinf%noffdiag_max*peinf%npools 930 jj=(ii-1)/peinf%npools+1 931 kk=mod(ii-1,peinf%npools) 932 if (peinf%inode/(peinf%npes/peinf%npools).eq.kk) then 933 if (ii.le.sig%noffdiag) then 934 peinf%index_offdiag(jj)=ii 935 peinf%flag_offdiag(jj)=.true. 936 else 937 peinf%index_offdiag(jj)=1 938 peinf%flag_offdiag(jj)=.false. 939 endif 940 endif 941 enddo 942 endif 943 944 SAFE_ALLOCATE(peinf%indext, (peinf%ntband_max)) 945 946!JRD Now is over procs per pool 947 948 SAFE_ALLOCATE(peinf%indext_dist, (peinf%ntband_max,peinf%npes_pool)) 949 SAFE_ALLOCATE(peinf%ntband_dist, (peinf%npes_pool)) 950 951 peinf%ntband_node=0 952 peinf%ntband_dist=0 953 peinf%nvband_node=0 954 peinf%indext=0 955 peinf%indext_dist=0 956 957 do ii=1,sig%ntband 958 ! FHJ: The bands are assigned to the processors in a sequential fashion. 959 ! This is more efficient for the parallel IO. Eg: for 2 bands/processor, 960 ! v1->PE0, v2->PE0, v3->PE1, etc. 961 ipe = (ii-1)/peinf%ntband_max 962 if (peinf%pool_rank.eq.ipe) then 963 peinf%ntband_node=peinf%ntband_node+1 964 if(ii.le.(sig%nvband+sig%ncrit)) then 965 peinf%nvband_node=peinf%nvband_node+1 966 endif 967 peinf%indext(peinf%ntband_node)=ii 968 endif 969 peinf%ntband_dist(ipe+1)=peinf%ntband_dist(ipe+1)+1 970 peinf%indext_dist(peinf%ntband_dist(ipe+1),ipe+1)=ii 971 enddo 972#ifdef MPI 973 if (peinf%verb_debug) then 974 call logit('Begin distrib report') 975 FLUSH(0) 976 FLUSH(6) 977 call MPI_Barrier(MPI_COMM_WORLD, mpierr) 978 if (peinf%inode==0) then 979 do ipe=1,peinf%npes_pool 980 write(6,'(/,a,i0,2x,a,i0,a)') '### pool_ipe=',ipe,' ; num bands=',peinf%ntband_dist(ipe),' ###' 981 do ii=1,peinf%ntband_dist(ipe) 982 write(6,'(a,i0,a,i0)') 'local:',ii,' -> global:', peinf%indext_dist(ii,ipe) 983 enddo 984 enddo 985 write(6,*) 986 endif 987 FLUSH(6) 988 call MPI_Barrier(MPI_COMM_WORLD, mpierr) 989 call logit('End distrib report') 990 endif 991#endif 992 993!------------------------------------------------------------------------ 994! Report distribution of bands over processors 995 996 if (peinf%inode==0) then 997 write(6,*) 998 write(6,'(1x,a)') 'Parallelization report:' 999701 format(1x,'- Using ',i0,' processor(s), ' ,i0,' pool(s), ',i0,' processor(s) per pool.') 1000 write(6,701) peinf%npes, peinf%npools, & 1001 peinf%npes/peinf%npools 1002702 format(1x,' - Note: distribution is not ideal; ',i0,' processor(s) is/are idle.') 1003 if (mod(peinf%npes,peinf%npools).ne.0) & 1004 write(6,702) peinf%npes-(peinf%npes/peinf%npools)* & 1005 peinf%npools 1006 if (mod(sig%ndiag,peinf%npools).eq.0) then 1007703 format(1x,'- Each pool is computing ',i0,' diagonal sigma matrix element(s).') 1008 write(6,703) peinf%ndiag_max 1009 else 1010704 format(1x,'- Each pool is computing ',i0,' to ',i0,' diagonal sigma matrix element(s).') 1011 write(6,704) peinf%ndiag_max-1, peinf%ndiag_max 1012705 format(1x,'- Note: distribution is not ideal because the number of diagonal sigma',/,& 1013 3x,'matrix elements (',i0,') is not divisible by the number of pools (',i0,').') 1014 write(6,705) sig%ndiag, peinf%npools 1015 endif 1016 if (mod(sig%noffdiag,peinf%npools).eq.0) then 1017706 format(1x,'- Each pool is computing ',i0,' off-diagonal sigma matrix element(s).') 1018 write(6,706) peinf%noffdiag_max 1019 else 1020707 format(1x,'- Each pool is computing ',i0,' to ',i0,' off-diagonal sigma matrix element(s).') 1021 write(6,707) peinf%noffdiag_max-1, peinf%noffdiag_max 1022708 format(1x,'- Note: distribution is not ideal because the number of off-diagonal sigma',/,& 1023 3x,'matrix elements (',i0,') is not divisible by the number of pools (',i0,').') 1024 write(6,708) sig%noffdiag, peinf%npools 1025 endif 1026 if (mod(sig%ntband,peinf%npes/peinf%npools).eq.0) then 1027709 format(1x,'- Each processor is holding ',i0, ' band(s).') 1028 write(6,709) peinf%ntband_max 1029 else 1030710 format(1x,'- Each pool is holding ',i0,' to ',i0,' band(s).') 1031 write(6,710) minval(peinf%ntband_dist), maxval(peinf%ntband_dist) 1032711 format(1x,'- Note: distribution is not ideal because the total number of bands',/,& 1033 3x,'(',i0,') is not divisible by the number of processors per pool (',i0,').') 1034 write(6,711) sig%ntband, peinf%npes/peinf%npools 1035 endif 1036 write(6,*) 1037 endif 1038 1039!----------------------------------------------------------- 1040! 1041! LOOP OVER K-POINT GRID AND READ IN WAVEFUNCTIONS 1042! 1043!----------------------------------------------------------- 1044 1045 if(sig%wfn_hdf5) then 1046 else 1047 ! ZL: we call read_wavefunctions() only once in input() 1048 if(.not.ep_read) then 1049 if(.not.do_phonq) then 1050 call read_wavefunctions(kp, gvec, sig, wfnk, iunit_c, iunit_k, fnc, fnk, wfnkqmpi, wfnkmpi) 1051 else 1052 call read_wavefunctions(kp, gvec, sig, wfnk, iunit_c, iunit_k, fnc, fnk, wfnkqmpi, wfnkmpi,& 1053 ep_read_in=.false., do_phonq_in=do_phonq, wfnk_phonq=wfnk_phonq, wfnk_phonq_mpi=wfnk_phonq_mpi) 1054 endif 1055 else 1056 call read_wavefunctions(kp, gvec, sig, wfnk, iunit_c, iunit_k, fnc, fnk, wfnkqmpi, wfnkmpi, ep_read_in = ep_read) 1057 endif 1058 if(peinf%inode.eq.0) then 1059 call close_file(25) 1060 endif 1061 endif 1062 1063 if(peinf%inode.eq.0) then 1064! Write out information about self-energy calculation 1065 call scissors_write(6, sig%scis) 1066 call scissors_write(6, sig%scis_outer, "outer") 1067 write(6,'(a)') 1068 1069 if(sig%freq_dep.eq.1) then 1070 do k=1,sig%nspin 1071 write(6,160) sig%spin_index(k),wpg%nelec & 1072 (sig%spin_index(k)),sqrt(wpg%wpsq(sig%spin_index(k))) 1073160 format(/,1x,'Data for sum rule:',/,& 1074 1x,'- rho(0,',i0,') = ',f0.6,' electrons',/,& 1075 1x,'- wp = ',f0.6,' eV',/) 1076 enddo 1077 endif ! sig%freq_dep 1078 write(6,'(1x,a,i0)') 'Number of bands to compute diagonal self-energy matrix elements: ', sig%ndiag 1079 write(6,'(1x,a)') 'Bands:' 1080 write(6,'(1x,"- ",i0)') sig%diag(:) 1081 write(6,'(1x,a,i0)') 'Number of off-diagonal bands to compute self-energy matrix elements: ', sig%noffdiag 1082 if (sig%noffdiag>0) then 1083 write(6,'(1x,a)') 'Off-diagonal bands:' 1084 do k=1,sig%noffdiag 1085174 format(1x,'- offmap(:, ',i6,') =',3i6) 1086 write(6,174) k, (sig%offmap(k,ii), ii = 1, 3) 1087 enddo 1088 endif 1089 write(6,*) 1090 1091 endif ! node 0 1092 1093 SAFE_DEALLOCATE(isrt) 1094 SAFE_DEALLOCATE(isrti) 1095 SAFE_DEALLOCATE_P(kp%w) 1096 SAFE_DEALLOCATE_P(kp%el) 1097 SAFE_DEALLOCATE_P(kp%elda) 1098 1099 POP_SUB(input) 1100 1101 return 1102 1103end subroutine input 1104 1105 1106subroutine read_wavefunctions(kp,gvec,sig,wfnk,iunit_c,iunit_k,fnc,fnk,wfnkqmpi,wfnkmpi, & 1107 ep_read_in,do_phonq_in,wfnk_phonq,wfnk_phonq_mpi) 1108 type (kpoints), intent(in) :: kp 1109 type (gspace), intent(in) :: gvec 1110 type (siginfo), intent(in) :: sig 1111 type (wfnkstates), intent(out) :: wfnk 1112 integer, intent(in) :: iunit_c, iunit_k 1113 character*20, intent(in) :: fnc, fnk 1114 type (wfnkqmpiinfo), intent(out) :: wfnkqmpi 1115 type (wfnkmpiinfo), intent(out) :: wfnkmpi 1116 ! ZL: for EP 1117 type (wfnkstates), optional, intent(out) :: wfnk_phonq ! similar to wfnk 1118 type (wfnkmpiinfo), optional, intent(out) :: wfnk_phonq_mpi ! similar to wfnkmpi 1119 1120 character :: tmpstr*100,tmpstr1*16,tmpstr2*16 1121 integer :: i,ii,j,k 1122 integer :: ikn, irk 1123 integer :: istore,nstore,inum,g1,g2,iknstore, istore_ep, iknstore_ep ! ZL adds istore_ep, iknstore_ep 1124 integer :: ndvmax, ndvmax_l 1125 integer, allocatable :: isort(:) 1126 real(DP) :: qk(3) 1127 SCALAR, allocatable :: zc(:,:) 1128 logical :: dont_read 1129 type(gspace) :: gvec_kpt 1130 1131 type(progress_info) :: prog_info !< a user-friendly progress report 1132 1133 ! ZL: electron phonon (EP) 1134 logical, intent(in), optional :: ep_read_in, do_phonq_in 1135 logical :: ep_read, do_phonq 1136 integer :: induse ! index in use, for indkn or indk_phonq 1137 integer :: tmpngk 1138 1139 PUSH_SUB(read_wavefunctions) 1140 1141 ! ZL: set up EP 1142 ep_read = .false. 1143 if(present(ep_read_in)) then 1144 ep_read = ep_read_in 1145 endif 1146 1147 do_phonq = .false. 1148 if(present(do_phonq_in)) then 1149 do_phonq = do_phonq_in 1150 endif 1151 1152 ! ZL: although it should never happen, in case someone mis-uses, double check here 1153 if(do_phonq .and. ep_read) then 1154 call die("Doing k+phonq points (designed for outer) is NOT allowed for dWFN (only used as inner).", & 1155 only_root_writes = .true.) 1156 endif 1157 1158 if (do_phonq) then 1159 if ((.not.present(wfnk_phonq)) .or. (.not.present(wfnk_phonq_mpi))) then 1160 call die("Doing k_phonq wfn needs wfklkrq and wfnk_phonq_mpi as arguments!", only_root_writes = .true.) 1161 endif 1162 else 1163 if (present(wfnk_phonq) .or. present(wfnk_phonq_mpi)) then 1164 call die("You are not doing k_phonq wfn, then do not pass wfnk_phonq or wfnk_phonq_mpi.", only_root_writes=.true.) 1165 endif 1166 endif 1167 ! ZL: done EP setup 1168 1169 SAFE_ALLOCATE(isort, (gvec%ng)) 1170 1171 SAFE_ALLOCATE(wfnkqmpi%nkptotal, (kp%nrk)) 1172 SAFE_ALLOCATE(wfnkqmpi%isort, (kp%ngkmax,kp%nrk)) 1173 SAFE_ALLOCATE(wfnkqmpi%band_index, (peinf%ntband_max,kp%nrk)) 1174 SAFE_ALLOCATE(wfnkqmpi%qk, (3,kp%nrk)) 1175 SAFE_ALLOCATE(wfnkqmpi%el, (sig%ntband,sig%nspin,kp%nrk)) 1176 ! ZL: ntband_max is total number of bands distributed to each processor 1177 ! wfnkqmpi%cg is distributed over bands 1178 SAFE_ALLOCATE(wfnkqmpi%cg, (kp%ngkmax,peinf%ntband_max,sig%nspin*kp%nspinor,kp%nrk)) 1179 1180 if (sig%nkn.gt.1) then ! ZL: same logic for EP 1181 ! it is using ngkmax, so same for EP, wfnkmpi%cg is distributed over ndv 1182 ndvmax=peinf%ndiag_max*kp%ngkmax 1183 if (mod(ndvmax,peinf%npes/peinf%npools).eq.0) then 1184 ndvmax_l=ndvmax/(peinf%npes/peinf%npools) 1185 else 1186 ndvmax_l=ndvmax/(peinf%npes/peinf%npools)+1 1187 endif 1188 SAFE_ALLOCATE(wfnkmpi%nkptotal, (sig%nkn)) 1189 SAFE_ALLOCATE(wfnkmpi%isort, (kp%ngkmax,sig%nkn)) 1190 SAFE_ALLOCATE(wfnkmpi%qk, (3,sig%nkn)) 1191 SAFE_ALLOCATE(wfnkmpi%el, (sig%ntband,sig%nspin,sig%nkn)) 1192 SAFE_ALLOCATE(wfnkmpi%elda, (sig%ntband,sig%nspin,sig%nkn)) 1193 SAFE_ALLOCATE(wfnkmpi%cg, (ndvmax_l,sig%nspin*kp%nspinor,sig%nkn)) 1194 1195 ! ZL: allocate wfnk_phonq_mpi 1196 ! if sig%nkn=1, then no mpi, only use wfnk_phonq which stores k+phonq 1197 if (do_phonq) then 1198 ! ZL: the size should in principle be sig%nkn*sig%nphonq, 1199 ! since we consider only nphonq=1, we omit it here 1200 SAFE_ALLOCATE(wfnk_phonq_mpi%nkptotal, (sig%nkn)) 1201 SAFE_ALLOCATE(wfnk_phonq_mpi%isort, (kp%ngkmax,sig%nkn)) 1202 SAFE_ALLOCATE(wfnk_phonq_mpi%qk, (3,sig%nkn)) 1203 SAFE_ALLOCATE(wfnk_phonq_mpi%el, (sig%ntband,sig%nspin,sig%nkn)) 1204 SAFE_ALLOCATE(wfnk_phonq_mpi%elda, (sig%ntband,sig%nspin,sig%nkn)) 1205 SAFE_ALLOCATE(wfnk_phonq_mpi%cg, (ndvmax_l,sig%nspin*kp%nspinor,sig%nkn)) 1206 endif 1207 endif 1208 1209!----------------------------------- 1210! Read in wavefunction information 1211 1212 SAFE_ALLOCATE(wfnk%isrtk, (gvec%ng)) 1213 SAFE_ALLOCATE(wfnk%ek, (sig%ntband,sig%nspin)) 1214 SAFE_ALLOCATE(wfnk%elda, (sig%ntband,sig%nspin)) 1215 1216 ! ZL: allocate wfnk_phonq 1217 if (do_phonq) then 1218 SAFE_ALLOCATE(wfnk_phonq%isrtk, (gvec%ng)) 1219 SAFE_ALLOCATE(wfnk_phonq%ek, (sig%ntband,sig%nspin)) 1220 SAFE_ALLOCATE(wfnk_phonq%elda, (sig%ntband,sig%nspin)) 1221 endif 1222 1223 if(.not.ep_read) then 1224 call progress_init(prog_info, 'reading wavefunctions (WFN_inner)', 'state', kp%nrk*sig%ntband) 1225 else 1226 call progress_init(prog_info, 'reading wavefunctions (dWFN)', 'state', kp%nrk*sig%ntband) 1227 endif 1228 1229 ! ZL: outmost loop, irk over kp%nrk 1230 do irk=1,kp%nrk 1231 if(.not.ep_read) then 1232 write(tmpstr,*) 'Reading WFN_inner -> cond/val wfns irk=',irk 1233 else 1234 write(tmpstr,*) 'Reading dWFN -> cond/val wfns irk=',irk 1235 endif 1236 call logit(tmpstr) 1237 qk(:)=kp%rk(:,irk) 1238 1239!---------------------------- 1240! Read in and sort gvectors 1241 SAFE_ALLOCATE(gvec_kpt%components, (3, kp%ngk(irk))) 1242 if ( peinf%inode == 0 ) call timing%start(timing%input_read) 1243 call read_binary_gvectors(25, kp%ngk(irk), kp%ngk(irk), gvec_kpt%components) 1244 if ( peinf%inode == 0 ) call timing%stop(timing%input_read) 1245 1246 do i = 1, kp%ngk(irk) 1247 call findvector(isort(i), gvec_kpt%components(:, i), gvec) 1248 if (isort(i) .eq. 0) then 1249 write(0,*) 'could not find gvec', kp%ngk(irk), i, gvec_kpt%components(1:3, i) 1250 call die('findvector') 1251 endif 1252 enddo 1253 SAFE_DEALLOCATE_P(gvec_kpt%components) 1254 1255!-------------------------------------------------------- 1256! Determine if Sigma must be computed for this k-point. 1257! If so, store the bands and wavefunctions on file iunit_k. 1258! If there is only one k-point, store directly in wfnk. 1259 1260 istore=0 1261 if(do_phonq) istore_ep=0 ! ZL: initialize 1262 do ikn=1,sig%nkn ! ZL: here we take advantage that we implement for one phonq point for now 1263 ! otherwize, need another outer loop: do ikn=1,sig%nkn*sig%nphonq 1264 1265 if(sig%indkn(ikn).eq.irk) then ! ZL: KEY if-statement: here it judges whether store wfnkmpi 1266 istore=1 ! ZL: for a given irk point, only one ind(ikn) index can match 1267 iknstore=ikn ! ZL: and this ind(ikn) is recorded as iknstore 1268 1269 wfnk%nkpt=kp%ngk(irk) 1270 wfnk%ndv=peinf%ndiag_max*kp%ngk(irk) 1271 wfnk%k(:)=qk(:) 1272 wfnk%isrtk(:)=isort(:) 1273 do k=1,sig%nspin 1274 wfnk%ek(1:sig%ntband,k)= & 1275 kp%el(1+ sig%ncore_excl:sig%ntband+ sig%ncore_excl,irk,sig%spin_index(k)) 1276 wfnk%elda(1:sig%ntband,k)= & 1277 kp%elda(1:sig%ntband,irk,sig%spin_index(k)) 1278 enddo 1279 SAFE_ALLOCATE(wfnk%zk, (wfnk%ndv,sig%nspin*kp%nspinor)) 1280 wfnk%zk=ZERO 1281 endif 1282 1283 ! ZL: for do_phonq 1284 if(do_phonq) then 1285 if (sig%indk_phonq(ikn).eq.irk) then 1286 istore_ep=1 1287 iknstore_ep=ikn 1288 1289 wfnk_phonq%nkpt=kp%ngk(irk) 1290 wfnk_phonq%ndv=peinf%ndiag_max*kp%ngk(irk) 1291 wfnk_phonq%k(:)=qk(:) + sig%indk_phonq_g0(:,ikn) 1292 wfnk_phonq%isrtk(:)=isort(:) 1293 1294 call shift_phonq_g0(sig%indk_phonq_g0(:,ikn), wfnk_phonq, gvec) 1295 1296 do k=1,sig%nspin 1297 wfnk_phonq%ek(1:sig%ntband,k)= & 1298 kp%el(1+ sig%ncore_excl:sig%ntband+ sig%ncore_excl,irk,sig%spin_index(k)) 1299 wfnk_phonq%elda(1:sig%ntband,k)= & 1300 kp%elda(1:sig%ntband,irk,sig%spin_index(k)) 1301 enddo 1302 SAFE_ALLOCATE(wfnk_phonq%zk, (wfnk_phonq%ndv,sig%nspin*kp%nspinor)) 1303 wfnk_phonq%zk=ZERO 1304 endif 1305 endif 1306 1307 enddo ! ikn 1308 1309 ! ZL: nkptotal is ngk 1310 wfnkqmpi%nkptotal(irk) = kp%ngk(irk) 1311 wfnkqmpi%isort(1:kp%ngk(irk),irk) = isort(1:kp%ngk(irk)) 1312 do k = 1, sig%nspin 1313 wfnkqmpi%el(1:sig%ntband,k,irk) = & 1314 kp%el(1+ sig%ncore_excl:sig%ntband+ sig%ncore_excl,irk,sig%spin_index(k)) 1315 enddo 1316 wfnkqmpi%qk(1:3,irk) = qk(1:3) 1317 1318!------------------------------------------------------------------------- 1319! SIB: Read wave functions from file WFN_inner (ZL: or dWFN) (unit=25) and have 1320! the appropriate processor write it to iunit_c after checking norm. 1321! The wavefunctions for bands where Sigma matrix elements are 1322! requested are stored in wfnk%zk for later writing to unit iunit_k. 1323! If band index is greater than sig%ntband, we will actually 1324! not do anything with this band (see code below). 1325! We still have to read it though, in order 1326! to advance the file to get to the data for the next k-point. 1327 1328 SAFE_ALLOCATE(zc, (kp%ngk(irk), kp%nspin*kp%nspinor)) 1329 1330 inum=0 1331 do i=1,kp%mnband 1332 1333 if ( peinf%inode == 0 ) call timing%start(timing%input_read) 1334 1335! DVF: don`t read deep core states. 1336 dont_read = (i > sig%ntband+sig%ncore_excl .or. i <= sig%ncore_excl) 1337 1338 call read_binary_data(25, kp%ngk(irk), kp%ngk(irk), kp%nspin*kp%nspinor, zc, dont_read = dont_read) 1339 1340 if ( peinf%inode == 0 ) call timing%stop(timing%input_read) 1341 1342 if (.not.dont_read) then 1343 call progress_step(prog_info, sig%ntband*(irk-1) + i) 1344 if(peinf%inode.eq.0) then 1345 do k=1,sig%nspin 1346 if(.not.ep_read) then 1347 call checknorm('WFN_inner',i,irk,kp%ngk(irk),k,kp%nspinor,zc(:,:)) 1348 endif 1349 enddo 1350 endif 1351 1352 nstore=0 ! ZL: for bands 1353 do j=1,peinf%ndiag_max 1354 if (.not.peinf%flag_diag(j)) cycle 1355 if (i==sig%diag(peinf%index_diag(j))) nstore=j 1356 enddo 1357 1358 ! ZL: if this is the kpoint in outer, save in wfnk 1359 if((istore.eq.1).and.(nstore.ne.0)) then 1360 do k=1,sig%nspin*kp%nspinor 1361 do j=1,kp%ngk(irk) 1362 wfnk%zk((nstore-1)*kp%ngk(irk)+j,k) = zc(j,sig%spin_index(k)) 1363 enddo 1364 enddo 1365 endif 1366 1367 ! ZL: for do_phonq 1368 if(do_phonq.and.(istore_ep.eq.1).and.(nstore.ne.0)) then 1369 do k=1,sig%nspin*kp%nspinor 1370 do j=1,kp%ngk(irk) 1371 wfnk_phonq%zk((nstore-1)*kp%ngk(irk)+j,k) = zc(j,sig%spin_index(k)) 1372 enddo 1373 enddo 1374 endif 1375 1376! DVF : i is indexed including the deep core states (since it runs over all the bands 1377! in the wavefunction), while the indexing arrays 1378! are not. We have to subtract sig%ncore_excl from i to get the right states. 1379 j=0 ! ZL: index to decide if save this band 1380 do ii=1,peinf%ntband_node 1381 if (peinf%indext(ii)==i-sig%ncore_excl) j=1 1382 enddo 1383 1384 if ( peinf%inode == 0 ) call timing%start(timing%input_write) 1385 1386 if (j.eq.1) then 1387 inum=inum+1 1388 wfnkqmpi%band_index(inum,irk)=i-sig%ncore_excl 1389 do k=1,sig%nspin*kp%nspinor 1390 wfnkqmpi%cg(1:kp%ngk(irk),inum,k,irk)= & 1391 zc(1:kp%ngk(irk),sig%spin_index(k)) 1392 enddo 1393 endif 1394 1395 if ( peinf%inode == 0 ) call timing%stop(timing%input_write) 1396 else 1397 ! FHJ: the following lines were introduced in r6294 and are supposed to 1398 ! be a shortcut if we are past the last band of the last k-point. However, 1399 ! in light of a previous bug (#223), this feature is commented out for now. 1400 !! FHJ: shortcut if this is past the last band of the last k-point 1401 !if (irk==kp%nrk) exit 1402 endif 1403 1404 enddo ! i (loop over bands) ZL: kp%mnband 1405 SAFE_DEALLOCATE(zc) 1406 1407 ! ZL: this is the needed kpoint for wfnkmpi, copy from wfnk to wfnkmpi 1408 if((istore.eq.1).and.(sig%nkn.gt.1)) then ! ZL: here it stores wfnkmpi 1409 1410 if ( peinf%inode == 0 ) call timing%start(timing%input_write) 1411 1412 ikn=iknstore 1413 wfnkmpi%nkptotal(ikn)=kp%ngk(irk) 1414 wfnkmpi%isort(1:kp%ngk(irk),ikn)=wfnk%isrtk(1:kp%ngk(irk)) 1415 wfnkmpi%qk(1:3,ikn)=qk(1:3) 1416 wfnkmpi%el(1:sig%ntband,1:sig%nspin,ikn)= & 1417 wfnk%ek(1:sig%ntband,1:sig%nspin) 1418 wfnkmpi%elda(1:sig%ntband,1:sig%nspin,ikn)= & 1419 wfnk%elda(1:sig%ntband,1:sig%nspin) 1420 do k=1,sig%nspin*kp%nspinor 1421#ifdef MPI 1422 i=mod(peinf%inode,peinf%npes/peinf%npools) 1423 if (mod(wfnk%ndv,peinf%npes/peinf%npools).eq.0) then 1424 j=wfnk%ndv/(peinf%npes/peinf%npools) 1425 else 1426 j=wfnk%ndv/(peinf%npes/peinf%npools)+1 1427 endif 1428 g1=1+i*j 1429 g2=min(j+i*j,wfnk%ndv) 1430 if (g2.ge.g1) then 1431 wfnkmpi%cg(1:g2-g1+1,k,ikn)=wfnk%zk(g1:g2,k) 1432 endif ! g2.ge.g1 1433#else 1434 wfnkmpi%cg(1:wfnk%ndv,k,ikn)=wfnk%zk(1:wfnk%ndv,k) 1435#endif 1436 enddo 1437 1438 if ( peinf%inode == 0 ) call timing%stop(timing%input_write) 1439 1440 ! ZL: if sig%nkn.gt.1, clean wfnk%zk, otherwise (only one kpoint in outer) we keep it 1441 SAFE_DEALLOCATE_P(wfnk%zk) 1442 endif 1443 1444 ! ZL: for do_phonq 1445 if(do_phonq.and.(istore_ep.eq.1).and.(sig%nkn.gt.1)) then 1446 1447 if ( peinf%inode == 0 ) call timing%start(timing%input_write) 1448 1449 ikn=iknstore_ep 1450 wfnk_phonq_mpi%nkptotal(ikn)=kp%ngk(irk) 1451 wfnk_phonq_mpi%isort(1:kp%ngk(irk),ikn)=wfnk_phonq%isrtk(1:kp%ngk(irk)) 1452 wfnk_phonq_mpi%qk(1:3,ikn)=wfnk_phonq%k(:) 1453 wfnk_phonq_mpi%el(1:sig%ntband,1:sig%nspin,ikn)= & 1454 wfnk_phonq%ek(1:sig%ntband,1:sig%nspin) 1455 wfnk_phonq_mpi%elda(1:sig%ntband,1:sig%nspin,ikn)= & 1456 wfnk_phonq%elda(1:sig%ntband,1:sig%nspin) 1457 do k=1,sig%nspin*kp%nspinor 1458#ifdef MPI 1459 i=mod(peinf%inode,peinf%npes/peinf%npools) 1460 if (mod(wfnk_phonq%ndv,peinf%npes/peinf%npools).eq.0) then 1461 j=wfnk_phonq%ndv/(peinf%npes/peinf%npools) 1462 else 1463 j=wfnk_phonq%ndv/(peinf%npes/peinf%npools)+1 1464 endif 1465 g1=1+i*j 1466 g2=min(j+i*j,wfnk_phonq%ndv) 1467 if (g2.ge.g1) then 1468 wfnk_phonq_mpi%cg(1:g2-g1+1,k,ikn)=wfnk_phonq%zk(g1:g2,k) 1469 endif ! g2.ge.g1 1470#else 1471 wfnk_phonq_mpi%cg(1:wfnk_phonq%ndv,k,ikn)=wfnk_phonq%zk(1:wfnk_phonq%ndv,k) 1472#endif 1473 enddo 1474 1475 if ( peinf%inode == 0 ) call timing%stop(timing%input_write) 1476 1477 ! ZL: if sig%nkn.gt.1, remove wfnk%zk, otherwise (only one kpoint in outer) we keep it 1478 SAFE_DEALLOCATE_P(wfnk_phonq%zk) 1479 endif 1480 1481 enddo ! irk (loop over k-points) 1482 call progress_free(prog_info) 1483 1484 SAFE_DEALLOCATE(isort) 1485 1486 PUSH_SUB(read_wavefunctions) 1487 1488end subroutine read_wavefunctions 1489 1490 1491end module input_m 1492 1493 1494