1!============================================================================== 2! 3! Routines: 4! 5! (1) epsilon Originally by (MSH) Last Edited 5/1/2008 (JRD) 6! 7! Send comments/bugs to jdeslip@berkeley.edu 8! 9! See README file for more details 10! 11!============================================================================== 12 13#include "f_defs.h" 14 15program epsilon 16 17 use global_m 18 use blas_m 19 use chi_summation_m 20 use chi_convergence_m 21 use fftw_m 22 use fullbz_m 23 use genwf_eps_m 24 use gmap_m 25 use input_m 26 use input_q_m 27 use input_utils_m 28 use irrbz_m 29 use mtxel_m 30 use mtxelmultiply_m 31 use read_matrix_m 32 use scalapack_m 33 use sort_m 34 use vcoul_generator_m 35 use write_matrix_m 36 use io_utils_m 37 use epswrite_hdf5_m 38 use tile_m 39#ifdef HDF5 40 use hdf5 41#endif 42 use timing_m, only: common_timing, timing => epsilon_timing 43 44 implicit none 45 46 type (kpoints) :: kp 47 type (kpoints) :: kpq 48 type (symmetry) :: syms 49 type (gspace) :: gvec 50 type (crystal) :: crys 51 type (polarizability) :: pol 52 type (valence_wfns) :: vwfn 53 type (conduction_wfns) :: cwfn 54 type (scalapack) :: scal 55 type (int_wavefunction) :: intwfnv 56 type (int_wavefunction) :: intwfnvq 57 type (int_wavefunction) :: intwfnc 58 type(chi_summator_t) :: chi_summator 59 type(chi_converger_t) :: chi_converger 60 type(wfn_FFT_comm_t) :: wfn_FFT_comm 61 62!----------------------- 63! Arrays for kpoints (fullbz, ...etc) 64 65 integer :: nrk 66 integer :: indst(48) 67 integer, allocatable :: indrk(:),neq(:) 68 type(grid) :: gr 69 70!----------------------- 71! Arrays for polarizability matrix 72 73 integer :: nstar 74 logical :: is_q0, use_WFNq, qpt_done 75 integer :: npcmax,nprmax,ivin,neqmax 76 integer, allocatable :: ind(:),indt(:,:,:) 77 integer, allocatable :: nprdtemp(:),npcdtemp(:) 78 real(DP) :: qq(3) !< The current q-pt under consideration 79 real(DP) :: omega_plasma, kfact 80 real(DP), allocatable :: ekin(:), kweights(:), kvols(:) 81 SCALAR, allocatable :: ph(:) 82 SCALAR, allocatable :: pht(:,:,:) 83 integer, allocatable :: nst(:) 84 85!----------------------- 86! Variables used with HDF5 87 88#ifdef HDF5 89 integer :: hdf5_error 90 integer :: my_iq 91 integer, allocatable :: isorti(:) 92 character(len=13) :: filename_chi_hdf5, filename_eps_hdf5, filename_out_hdf5 93#endif 94 95 character :: tmpstr*100 96 character :: filename*13 97 integer :: initial_access = 0 98 integer :: i,j,k,n,irk,iv,itran,it 99 integer :: ncount,ix,jx,kgq(3) 100 integer :: np,iunit,iq 101 integer :: ig_l, ig_g, ipe 102 integer :: ia, ib, ic, id, ie, if 103 integer :: nmtx_t 104 real(DP) :: tsec(2) 105 real(DP) :: fact,rk(3) 106 integer :: ispin 107 character*24 :: routnam(120) 108 integer, allocatable :: routsrt(:) 109 integer :: Ntime 110 integer :: iunit_v, iunit_c 111 real(DP) :: dnpcr,dnpcrmax 112 real(DP) :: E_rpa 113 114 integer :: jj 115 116 type(progress_info) :: prog_info !< a user-friendly progress report 117 118!--------------- Begin Program --------------------------------------- 119 120 call peinfo_init() 121 122!---------------------- 123! Initialize random numbers 124 125 peinf%jobtypeeval = 0 126 127!-------------------- 128! Initialize timer 129 call timing%init() 130 call common_timing%init() 131 if(peinf%inode .eq. 0) then 132 call timing%start(timing%total) 133 endif 134 135!------------------------ 136! Initialize files 137 138 if(peinf%inode .eq. 0) then 139 call timing%start(timing%job_setup) 140 endif 141 142 call open_file(55,file='epsilon.inp',form='formatted',status='old') 143 if(peinf%inode .eq. 0) then 144 call open_file(7,file='epsilon.log',form='formatted',status='replace') 145 endif 146 147 call write_program_header('Epsilon', .true.) 148 149!----------- Call Input: Read crystal data from unit 25 --------------- 150 151! read parameters and q-points from unit 55 (input file) 152! initialize unit 10 (output for polarizability matrix) 153 154 if(peinf%inode .eq. 0) then 155 call timing%stop(timing%job_setup) 156 endif 157 158 if(peinf%inode .eq. 0) then 159 call timing%start(timing%input) 160 endif 161 162 call input(kp,crys,syms,gvec,pol,cwfn,vwfn,intwfnv,intwfnc,omega_plasma,gr) 163 ! FHJ: consistency check 164 if (pol%os_opt_ffts/=0 .and. (kp%nspin*kp%nspinor/=1 .or. pol%ncrit/=0)) then 165 call die('Cannot use os_opt_fft/=0 for metals or spin-polarized calculations', & 166 only_root_writes=.true.) 167 endif 168 pol%FFTgrid = gvec%FFTgrid 169 if (pol%min_fftgrid) then 170 ! FHJ: Figure our the FFT grid that holds the WFNs 171 call get_wfn_fftgrid(pol, gvec, kp, intwfnv) 172 endif 173 174 if(.not. pol%skip_chi .and. peinf%inode == 0) then 175 call open_file(17,file='chi_converge.dat',form='formatted',status='replace') 176 endif 177 178 if (pol%iwritecoul .eq. 1) then 179 if (peinf%inode .eq. 0) then 180 call open_file(19,file='vcoul',form='formatted',status='replace') 181 endif 182 endif 183 184! CHP: saves the (G=0,G`=0) component of (retarded) epsilon inverse 185 if(peinf%inode==0 .and. pol%freq_dep>0 .and. .not.pol%skip_epsilon) then 186 call open_file(51,file='EpsInvDyn',form='formatted',status='replace') 187 call open_file(52,file='EpsDyn',form='formatted',status='replace') 188 endif 189 190 SAFE_ALLOCATE(vwfn%isort, (gvec%ng)) 191 SAFE_ALLOCATE(cwfn%isort, (gvec%ng)) 192 193 if (pol%nq0>0) then 194 ! FHJ: no q->0 point can have all coordinates set to zero 195 if (pol%icutv==TRUNC_NONE .and. any(all(abs(pol%qpt(:,:pol%nq0))<TOL_ZERO,dim=1))) then 196 call die('No truncation and zero q0', only_root_writes=.true.) 197 endif 198 endif 199 200 if(peinf%inode .eq. 0) then 201 call timing%stop(timing%input) 202 endif 203 204 205!-------------- Read wavefunctions for (k+q) points --------------------- 206 207! SIB: The input_q routine looks the same as the input routine but 208! if reads from a file called WFNq instead of WFN. Presumably 209! these are the shifted (by "q") wave functions. 210 211 if(peinf%inode .eq. 0) then 212 call timing%start(timing%input_q) 213 endif 214 215 if (pol%need_WFNq) then 216 if (peinf%inode .eq. 0) then 217 write(6,*) 'You have a slightly shifted q0 vector and a semiconductor.' 218 write(6,*) 'So, reading from WFNq.' 219 endif 220 call input_q(gvec,kpq,cwfn,vwfn,pol,intwfnvq) 221 elseif (pol%lin_denominator>TOL_Zero) then 222 endif 223 224 if(peinf%inode .eq. 0) then 225 call timing%stop(timing%input_q) 226 endif 227 228 229!-------------- GENERATE FULL BZ ---------------------------------------- 230 231! SIB: fullbz() takes the kpoints in kp%components(1:3,kp%nrk) and applies all 232! the symmetries in syms to them. The resulting set of unique vectors 233! are in gr%f(1:3,gr%nf) (gr%nf of them). 234 235 if (peinf%inode .eq. 0) then 236 write(6,'(1x,a)') 'Summary of the WFN files:' 237 write(6,'(1x,a,i0)') '- Number of k-points in WFN: ', kp%nrk 238 if (pol%need_WFNq) then 239 write(6,'(1x,a,i0)') '- Number of k-points in WFNq: ', kpq%nrk 240 endif 241 write(6,'(1x,a,i0)') '- Number of k-points in the full BZ of WFN: ', gr%nf 242 if (peinf%verb_high) then 243 write(6,'(1x,a,i0,a)') '- Full list of k-points generated with ',syms%ntran,' symmetries:' 244 write(6,'(/7x,a,6x)') 'k-point rk (irr. BZ)' 245 write(6,'(1x,32("-"))') 246 write(6,'(3(1x,f10.6))') (kp%rk(:,iq), iq=1, kp%nrk) 247 write(6,'(/7x,a,6x)') 'k-point fk (full BZ)' 248 write(6,'(1x,32("-"))') 249 write(6,'(3(1x,f10.6))') (gr%f(:,iq), iq=1, gr%nf) 250 endif 251 write(6,'()') 252 endif 253 254#ifdef USEVORO 255 if (pol%non_uniform) then 256 SAFE_ALLOCATE(kvols, (gr%nf)) 257 call get_kpts_volumes(crys%bdot, gr%f, kvols) 258 endif 259#endif 260 SAFE_ALLOCATE(ekin, (gvec%ng)) 261 SAFE_ALLOCATE(scal%nprd, (peinf%npes_freqgrp)) 262 SAFE_ALLOCATE(scal%npcd, (peinf%npes_freqgrp)) 263 264 if (pol%os_opt_ffts==2) then 265 ! FHJ: FFTs opt. level 2 => do all the FFTs using all the processor, int_wfn arrays 266 call genwf_FFT_Isend(wfn_FFT_comm,crys,gvec,syms,kp,kpq,vwfn,pol,cwfn,intwfnv,intwfnvq,intwfnc) 267 !call genwf_FFT(crys,gvec,syms,kp,kpq,vwfn,pol,cwfn,intwfnv,intwfnvq,intwfnc,need_WFNq) 268 endif 269 270!----------- LOOP over q points for which chi and eps are calculated ----- 271 do iq=1,pol%nq 272 273 if (pol%stop_after_qpt>-1 .and. iq>pol%stop_after_qpt) then 274 if (peinf%inode==0) write(6,'(/,1x,a,/)') 'stop_after_qpt: emulating a sudden application stop.' 275 FLUSH(6) 276 FLUSH(0) 277#ifdef MPI 278 call MPI_Barrier(MPI_COMM_WORLD, mpierr) 279 call MPI_Finalize(mpierr) 280#endif 281 MYSLEEP(1) 282 stop 283 endif 284 285 ! SIB: qq(1:3) is the current q vector under consideration 286 qq(:)=pol%qpt(:,iq) 287 if(peinf%inode.eq.0) then 288 call print_dealing_with(iq, pol%nq, qq, 'q') 289 endif 290 is_q0 = iq<=pol%nq0 291 use_WFNq = (is_q0.and.pol%need_WFNq).or.pol%patched_sampling.or.(pol%qflags(iq)==-1) 292 if (peinf%inode==0) then 293 if (is_q0) then 294 write(6,'(1x,a)') 'This is the special q->0 point.' 295 else 296 write(6,'(1x,a)') 'This is a regular non-zero q-point.' 297 endif 298 endif 299 300#ifdef HDF5 301 my_iq = iq 302 if (.not.is_q0) my_iq = iq - pol%nq0 303 if (is_q0) then 304 filename_eps_hdf5 = 'eps0mat.h5' 305 filename_chi_hdf5 = 'chi0mat.h5' 306 else 307 filename_eps_hdf5 = 'epsmat.h5' 308 filename_chi_hdf5 = 'chimat.h5' 309 endif 310 if (pol%skip_epsilon) then 311 filename_out_hdf5 = filename_chi_hdf5 ! Write to/restart chimat file 312 else 313 filename_out_hdf5 = filename_eps_hdf5 ! Write to/restart epsmat file 314 endif 315 if (pol%use_hdf5.and.pol%restart) then 316 if (peinf%inode==0) qpt_done = is_qpt_done(TRUNC(filename_out_hdf5), my_iq) 317#ifdef MPI 318 call MPI_BCAST(qpt_done, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, mpierr) 319#endif 320 if (qpt_done) then 321 if (peinf%inode==0) then 322 write(6,'(/,1x,a,/)') 'This q-point was already calculated: skipping.' 323 endif 324 cycle 325 endif 326 endif 327#endif 328 329 if(peinf%inode .eq. 0) then 330 call timing%start(timing%q_loop_setup) 331 endif 332 333!-------------------- 334! Determine number of matrix elements 335! 336! Calculate kinetic energies |q+G|^2 337 338 if (is_q0) then 339 call kinetic_energies(gvec, crys%bdot, ekin) 340 else 341 call kinetic_energies(gvec, crys%bdot, ekin, qvec = qq) 342 endif 343 344!-------------------- 345! Sort kinetic energies 346! index of ordered kinetic energies in array isrtx 347! 348! SIB: pol%isrtx has the indices for sorted ekin 349 350 SAFE_ALLOCATE(pol%isrtx, (gvec%ng)) 351 352 if(peinf%inode .eq. 0) then 353 call timing%stop(timing%q_loop_setup) 354 endif 355 356 if(peinf%inode .eq. 0) then 357 call timing%start(timing%gvec) 358 endif 359 call logit('sorting gvec') 360 call sortrx(gvec%ng, ekin, pol%isrtx, gvec = gvec%components) 361 if(peinf%inode .eq. 0) then 362 call timing%stop(timing%gvec) 363 endif 364 365!--------------------- 366! Compute inverse array to isrtx 367 368 SAFE_ALLOCATE(pol%isrtxi, (gvec%ng)) 369 do i=1,gvec%ng 370 pol%isrtxi(pol%isrtx(i))=i 371 enddo 372 373! SIB: pol%nmtx becomes the number of matrix elements to be computed; 374! the matrix is computed if its ekin is < ecuts 375 376 if(peinf%inode .eq. 0) then 377 call timing%start(timing%init_cutoff) 378 endif 379 380 pol%nmtx = gcutoff(gvec%ng, ekin, pol%isrtx, pol%ecuts) 381 if(peinf%inode.eq.0) then 382 write(6, '(1x,a,i0)') 'Rank of the polarizability matrix (nmtx): ', pol%nmtx 383 endif 384 ! FHJ: Do we want to use the economical fftgrid/box? If so, we pad the WFN FFTbox 385 ! by the box that holds nmtx gvectors, which is much smaller than the WFN fftbox. 386 if (pol%min_fftgrid.and.pol%os_opt_ffts<2) call get_eps_fftgrid(pol, gvec) 387 388 if(peinf%inode .eq. 0) then 389 call timing%stop(timing%init_cutoff) 390 endif 391 392 if(peinf%inode .eq. 0) then 393 call timing%start(timing%init_scalapack) 394 endif 395 396 SAFE_ALLOCATE(pol%irow, (pol%nmtx)) 397 pol%irow=0 398 399! JRD: Determine size of distributed matrices 400! DVF : when using parallel frequencies and the number of processors is not divisible 401! by the number of frequencies done in parallel, there are excess processors that we 402! don`t want to include in our distributed matrix operations. So, for these processors 403! we call blacs_setup and then reset the values of the ScaLAPACK variables 404! to harmless values that won`t affect any of the global variables obatined via MPI_Allreduce. 405 406 if(pol%nfreq_group .gt. 1) then 407 if (peinf%inode .lt. peinf%npes) then 408 call blacs_setup(scal, pol%nmtx, .false.,peinf%npes_freqgrp,pol%nfreq_group) 409 else 410 call blacs_setup(scal, pol%nmtx, .false.,peinf%npes_freqgrp,pol%nfreq_group,peinf%npes_orig-peinf%npes) 411 scal%npr=0 412 scal%npc=0 413 scal%nbl=1 414 scal%nprow=1 415 scal%npcol=1 416 scal%myprow=1000000 ! DVF: A very large number so that we never take anything from 417 scal%mypcol=1000000 ! these processors in the garbage frequency/mtxel groups 418 ! See where these variables are used in epsinv.f90 to see 419 ! what I mean. 420 endif 421 else 422 call blacs_setup(scal, pol%nmtx, .true.) 423 endif 424#ifdef USESCALAPACK 425 call logit('Initializing ScaLAPACK') 426#endif 427 428#ifdef MPI 429 SAFE_ALLOCATE(nprdtemp, (peinf%npes_freqgrp)) 430 SAFE_ALLOCATE(npcdtemp, (peinf%npes_freqgrp)) 431 432 scal%nprd = 0 433 scal%npcd = 0 434 nprdtemp = 0 435 npcdtemp = 0 436 437! write(*,*) "before nrpdtemp allreduce",peinf%inode 438 if(pol%nfreq_group .eq. 1) then 439 nprdtemp(peinf%inode + 1) = scal%npr 440 npcdtemp(peinf%inode + 1) = scal%npc 441 call MPI_ALLREDUCE(nprdtemp,scal%nprd,peinf%npes_freqgrp,MPI_INTEGER,MPI_SUM, MPI_COMM_WORLD,mpierr) 442 call MPI_ALLREDUCE(npcdtemp,scal%npcd,peinf%npes_freqgrp,MPI_INTEGER,MPI_SUM, MPI_COMM_WORLD,mpierr) 443 elseif(pol%nfreq_group .gt. 1) then 444 if(peinf%inode .lt. peinf%npes) then 445 nprdtemp(peinf%rank_f + 1) = scal%npr 446 npcdtemp(peinf%rank_f + 1) = scal%npc 447 call MPI_ALLREDUCE(nprdtemp,scal%nprd,peinf%npes_freqgrp,MPI_INTEGER,MPI_SUM, peinf%freq_comm,mpierr) 448 call MPI_ALLREDUCE(npcdtemp,scal%npcd,peinf%npes_freqgrp,MPI_INTEGER,MPI_SUM, peinf%freq_comm,mpierr) 449 else 450 endif 451 endif 452 SAFE_DEALLOCATE(nprdtemp) 453 SAFE_DEALLOCATE(npcdtemp) 454#else 455 scal%nprd = scal%npr 456 scal%npcd = scal%npc 457#endif 458 459! DVF: Get the maximum number of columns/rows owned by one of the processors. This is so you can 460! allocate arrays of the right size. For what determines how many rows and columns a procesor 461! has, see blacs_setup routine in Common directory and google the numroc routine of scalapack 462! Numroc is a nifty little routine 463 464 dnpcr = scal%npc*1D0*scal%npr 465 466! write(*,*) "before npr allreduce",peinf%inode 467#ifdef MPI 468 call MPI_ALLREDUCE(scal%npc,npcmax,1,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,mpierr) 469 call MPI_ALLREDUCE(scal%npr,nprmax,1,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,mpierr) 470 call MPI_ALLREDUCE(dnpcr,dnpcrmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,MPI_COMM_WORLD,mpierr) 471#else 472 npcmax = scal%npc 473 nprmax = scal%npr 474#endif 475 476 ! FHJ: I think there`s actually a but in this report.. 477 if (dnpcr>(pol%nmtx*1.5D0*pol%nmtx/(1D0*peinf%npes)) .and. & 478 peinf%inode==0 .and. peinf%verb_high) then 479 write(6,'(/1x,a)') 'NOTE: ScaLAPACK layout is not well load-balanced.' 480 write(6,'(1x,a,f12.0/)') 'Max number of elements owned by one task is:', dnpcr 481 endif 482 483! Allocate scalapack arrays needed later. See scalapack in Common/scalapack.f90 to see what 484! these arrays hold 485 486 SAFE_ALLOCATE(scal%isrtxcol, (scal%npc)) 487 SAFE_ALLOCATE(scal%isrtxrow, (scal%npr)) 488 SAFE_ALLOCATE(scal%imycol, (scal%npc)) 489 SAFE_ALLOCATE(scal%imyrow, (scal%npr)) 490 SAFE_ALLOCATE(scal%imycolinv, (pol%nmtx)) 491 SAFE_ALLOCATE(scal%imyrowinv, (pol%nmtx)) 492 ! FHJ: FIXME - DVF should be sprayed hard for overwriting peinf%npes!! 493 SAFE_ALLOCATE(scal%imycold, (npcmax,peinf%npes_orig)) 494 SAFE_ALLOCATE(scal%imyrowd, (nprmax,peinf%npes_orig)) 495 scal%imycold = 0 496 scal%imyrowd = 0 497 498 ! scal%imyrow(ig_l) = indxl2g(ig_l, scal%nbl, scal%myprow, 0, scal%nprow) 499 ! ipe = indxg2p(ig_g, scal%nbl, 0, 0, scal%npcol) 500 ! scal%imyrowd(ig_l, inode+1) = indxl2g(ig_l, scal%nbl, ipe+1, 0, scal%nprow) 501 ! FHJ: FIXME - DVF should be sprayed hard for overwriting peinf%npes!! 502 if (peinf%inode<peinf%npes_orig) then 503 ! FHJ: FIXME - these indexing arrays are completely useless and should be 504 ! removed. Let`s not reinvent BLACS, plz! 505 do ig_l = 1, scal%npr 506 ig_g = indxl2g(ig_l, scal%nbl, scal%myprow, 0, scal%nprow) 507 scal%isrtxrow(ig_l) = pol%isrtx(ig_g) 508 scal%imyrow(ig_l) = ig_g 509 scal%imyrowd(ig_l, peinf%inode+1) = ig_g 510 scal%imyrowinv(ig_g) = ig_l 511 enddo 512 do ig_l = 1, scal%npc 513 ig_g = indxl2g(ig_l, scal%nbl, scal%mypcol, 0, scal%npcol) 514 scal%isrtxcol(ig_l) = pol%isrtx(ig_g) 515 scal%imycol(ig_l) = ig_g 516 scal%imycold(ig_l, peinf%inode+1) = ig_g 517 scal%imycolinv(ig_g) = ig_l 518 enddo 519 !do ig_g = 1, pol%nmtx 520 ! ig_l = indxl2g(ig_g, scal%nbl, 0, 0, scal%npcol) 521 ! ipcol = indxg2p(ig_g, scal%nbl, 0, 0, scal%npcol) 522 ! scal%imycold(ig_l,ipe+1) = ig_g 523 ! ipe = indxg2p(ig_g, scal%nbl, 0, 0, scal%nprow) 524 ! ig_l = indxl2g(ig_g, scal%nbl, 0, 0, scal%nprow) 525 ! scal%imyrowd(ig_l,ipe+1) = ig_g 526 !enddo 527 endif 528#ifdef MPI 529 call MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, scal%imyrowd, & 530 size(scal%imyrowd,1), MPI_INTEGER, MPI_COMM_WORLD, mpierr) 531 call MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, scal%imycold, & 532 size(scal%imycold,1), MPI_INTEGER, MPI_COMM_WORLD, mpierr) 533#endif 534 535 if(peinf%inode .eq. 0) then 536 call timing%stop(timing%init_scalapack) 537 endif 538 539 if(peinf%inode .eq. 0) then 540 call timing%start(timing%init_arrays) 541 endif 542 543!---------------------- 544! Determine subgroup which leaves qq invariant 545! 546! SIB: figures out which symmetries acting on qq result in qq + integer 547! entries. syms%ntranq is their number, syms%indsub are their indices 548! (pointing to syms%mtrx), and syms%kgzero(1:3,:) are the integer entries. 549 550 if(peinf%inode .eq. 0) then 551 call timing%start(timing%subgrp) 552 endif 553 554 call subgrp(qq,syms) 555 if (pol%patched_sampling) then 556 syms%ntranq = 1 557 endif 558 559 if(peinf%inode .eq. 0) then 560 call timing%stop(timing%subgrp) 561 endif 562 563!----------------------- 564! Determine independent elements of polarizability matrix 565! 566! SIB: independent means due to symmetries. This initializes 567! the polarizability matrix pol%chi to zero (for independent entries) 568! and figure out phases due to symmetries for dependent ones, 569! and points dependent ones to the entries they depend on (pol%kxi indices) 570 571! JRD: we don`t do this anymore 572 573! if(peinf%inode .eq. 0) then 574! call timing%start(timing%wedontdothisanymore) 575! endif 576! call logit('calling indep') 577! call indep(nind,gvec,syms,pol,kp%nspin) 578! 579! JRD: Testing what if we set pol%kxi to zero 580! pol%kxi = 0 581! pol%chi = 0D0 582! nind=pol%nmtx*(pol%nmtx+1)/2 583! 584! if(peinf%inode .eq. 0) then 585! call timing%stop(timing%wedontdothisanymore) 586! endif 587 588!---------------------- 589! Reduce the k-points to the irr. bz with respect to q 590! 591! SIB: figure out k-points in irr. BZ (based on symmetries for current q) 592! nrk is # of irr. points, indrk are their indices in the full zone, 593! and neq is the number of equiv. points for an irr. point. 594! (Full zone vectors come in through gr%f(1:3,1:gr%nf).) 595 596 if(peinf%inode .eq. 0) then 597 call timing%start(timing%irrbz) 598 endif 599 600 SAFE_ALLOCATE(indrk, (gr%nf)) 601 SAFE_ALLOCATE(neq, (gr%nf)) 602 call irrbz(syms,gr%nf,gr%f,nrk,neq,indrk) 603#ifdef USEVORO 604 if (pol%non_uniform) then 605 SAFE_ALLOCATE(kweights, (nrk)) 606 do irk=1,nrk 607 ! FHJ: we divide by gr%nf afterwards... 608 kweights(irk) = kvols(indrk(irk))/sum(kvols) * dble(gr%nf) 609 enddo 610 if (peinf%inode==0) then 611 call open_file(666, file='kweights.dat', form='formatted', status='old') 612 write(666,'(/,1x,a,i0)') 'k-weights for iq=',iq 613 write(666,'(1x,i0,1x,f15.9)') (irk, kweights(irk), irk=1,nrk) 614 write(666,*) 615 call close_file(666) 616 endif 617 endif 618#endif 619 if(peinf%inode .eq. 0) then 620 call timing%stop(timing%irrbz) 621 endif 622 623 neqmax = maxval(neq(1:nrk)) 624 625! write(6,*) peinf%inode, 'neqmax', neq(1), neqmax 626 627!--------------------------- 628! Output points in irr. bz 629 630 if(peinf%inode.eq.0) then 631 write(6,'(1x,a,i0)') 'Number of k-points in the irreducible BZ(q) (nrk): ', nrk 632 if (peinf%verb_medium) then 633 write(6,'(/6x,a,5x,a)') 'k-point rk (irr. BZ)', '#eq/fBZ' 634 write(6,'(1x,29("-"),1x,7("-"))') 635 write(6,'(3(1x,f9.6),1x,i7)') (gr%f(1:3,indrk(j)), neq(j), j=1,nrk) 636 endif 637 write(7,70) (qq(i),i=1,3),pol%nmtx,nrk 63870 format(/ /,5x,'q=',3f7.4,2x,'nmtx=',i8,2x,'nrk=',i3) 639 endif 640 641 if (pol%patched_sampling) then 642 fact=4.0d0/(product(kp%kgrid)*crys%celvol*kp%nspin*kp%nspinor) 643 else 644 fact=4.0d0/(dble(gr%nf)*crys%celvol*kp%nspin*kp%nspinor) 645 endif 646 647 if (pol%freq_dep .eq. 0) then 648 SAFE_ALLOCATE(pol%chi, (scal%npr,scal%npc,kp%nspin)) 649 pol%chi=ZERO 650 endif 651 652 ! some check for subspace truncation method 653 IF(pol%subspace) THEN 654 IF(.NOT.(pol%freq_dep==2 .AND. pol%freq_dep_method==2)) THEN 655 call die('Subspace truncation implemented only for freq_dep=2 and freq_dep_method=2') 656 END IF 657 ! set this flag to false for the meantime, this will be used to decide 658 ! if regenerate the full chi or proceed in the calculation of epsilon 659 ! directly within the subspace 660 pol%need_full_chi = .FALSE. 661 END IF 662 663 ! allocate pol%chiRDyn later 664 IF(.NOT. pol%subspace) THEN 665 if ((pol%freq_dep .eq. 2).and.(pol%freq_dep_method .eq. 0 .or. pol%freq_dep_method .eq. 2)) then 666 SAFE_ALLOCATE(pol%chiRDyn, (scal%npr,scal%npc,pol%nfreq_in_group,kp%nspin)) 667 pol%chiRDyn=(0.0,0.0) 668 endif 669 END IF 670 671 if ((pol%freq_dep .eq. 2).and.(pol%freq_dep_method .eq. 1)) then 672 SAFE_ALLOCATE(pol%chiRDyn, (scal%npr,scal%npc,pol%nfreq_in_group,kp%nspin)) 673 pol%chiRDyn=(0.0,0.0) 674 SAFE_ALLOCATE(pol%chiTDyn, (pol%os_nsfreq_para,scal%npr,scal%npc,kp%nspin)) 675 pol%chiTDyn=(0.0,0.0) 676 endif 677 678 if (pol%freq_dep .eq. 3) then 679 SAFE_ALLOCATE(pol%chiRDyn, (scal%npr,scal%npc,pol%nfreq_in_group,kp%nspin)) 680 pol%chiRDyn=(0.0,0.0) 681 endif 682 683 if (.not. pol%skip_chi) then 684 685!------------------------------------ 686! SIB: allocate space 687 688! write(6,*) peinf%inode, 'allocating pht', neqmax, pol%nmtx, nrk 689 690 SAFE_ALLOCATE(ind, (pol%nmtx)) 691 SAFE_ALLOCATE(ph, (pol%nmtx)) 692 SAFE_ALLOCATE(pht, (pol%nmtx,neqmax,nrk)) 693 SAFE_ALLOCATE(indt, (pol%nmtx,neqmax,nrk)) 694 ind=0 695 696 SAFE_ALLOCATE(nst, (nrk)) 697 698! JRD: Possible Memory Hazard. We can speed this up by possibly 699! only allocating number of bands on current proc and doing send/recvs 700 if(pol%nfreq_group .gt. 1) then 701 if(pol%gcomm .eq. -1) then 702 SAFE_ALLOCATE(pol%gme, (pol%nmtx,peinf%ncownmax,peinf%nvownmax,kp%nspin,nrk,pol%nfreq_group)) 703!$OMP PARALLEL DO collapse(6) 704 do ia = 1, pol%nfreq_group 705 do ib = 1, nrk 706 do ic = 1 , kp%nspin 707 do id = 1, peinf%nvownmax 708 do ie = 1, peinf%ncownmax 709 do if = 1, pol%nmtx 710 pol%gme(if,ie,id,ic,ib,ia)=ZERO 711 enddo 712 enddo 713 enddo 714 enddo 715 enddo 716 enddo 717 else 718 SAFE_ALLOCATE(pol%gme, (pol%nmtx,peinf%ncownmax,peinf%nvownmax,kp%nspin,nrk,1)) 719!$OMP PARALLEL DO collapse(6) 720 do ia = 1, 1 721 do ib = 1, nrk 722 do ic = 1 , kp%nspin 723 do id = 1, peinf%nvownmax 724 do ie = 1, peinf%ncownmax 725 do if = 1, pol%nmtx 726 pol%gme(if,ie,id,ic,ib,ia)=ZERO 727 enddo 728 enddo 729 enddo 730 enddo 731 enddo 732 enddo 733 endif 734 else 735 SAFE_ALLOCATE(pol%gme, (pol%nmtx,peinf%ncownactual,peinf%nvownactual,kp%nspin,nrk,pol%nfreq_group)) 736!$OMP PARALLEL DO collapse(6) 737 do ia = 1, pol%nfreq_group 738 do ib = 1, nrk 739 do ic = 1 , kp%nspin 740 do id = 1, peinf%nvownactual 741 do ie = 1, peinf%ncownactual 742 do if = 1, pol%nmtx 743 pol%gme(if,ie,id,ic,ib,ia)=ZERO 744 enddo 745 enddo 746 enddo 747 enddo 748 enddo 749 enddo 750 endif 751 752 if (pol%freq_dep .eq. 2 .or. pol%freq_dep .eq. 3) then 753 if(pol%nfreq_group .eq. 1) then 754 SAFE_ALLOCATE(pol%edenDyn, (peinf%nvownactual,peinf%ncownactual,kp%nspin,nrk,pol%nfreq_group)) 755 else 756 SAFE_ALLOCATE(pol%edenDyn, (peinf%nvownmax,peinf%ncownmax,kp%nspin,nrk,pol%nfreq_group)) 757 endif 758 endif 759 760 if(peinf%inode .eq. 0) then 761 call timing%stop(timing%init_arrays) 762 endif 763 764!--------- LOOP OVER K-POINTS IN SET RK --------------------------------- 765 ! FHJ: this is to generate nice output / time estimate 766 call progress_init(prog_info, 'calculation of matrix elements', 'transition', nrk*peinf%nvownmax) 767 768! SIB: loop over points in irreducible zone 769 do irk=1,nrk 770 771 rk(:)=gr%f(:,indrk(irk)) ! rk(:) is the current irr. k-point 772! Regenerate star of rk,store the index of the rotation 773! SIB: Star is the set of vectors generated by applying all 774! subgroup operations for the current q-vector to the k-point rk. 775 776 if(peinf%inode .eq. 0) then 777 call timing%start(timing%rqstar) 778 endif 779 780 call rqstar(syms,nstar,indst,rk) 781 if(nstar.ne.neq(irk)) then 782 write(0,*) 'nstar?',irk,nstar,neq(irk) 783 call die('nstar mismatch') 784 endif 785 786 if(peinf%inode .eq. 0) then 787 call timing%stop(timing%rqstar) 788 endif 789! JRD: loop over transfs which generate the star of rk for gmap 790 791 nst(irk) = nstar 792 793 if(peinf%inode .eq. 0) then 794 call timing%start(timing%gmap) 795 endif 796 797 do it=1,nstar 798 799! Map g-vectors in polarizability to r**(-1)(g-gq) 800! note that gmap requires index of transf in full group 801! whereas indst gives index in subgroup 802 803 itran = syms%indsub(indst(it)) 804 kgq(:) = -syms%kgzero(:,indst(it)) ! note minus sign 805 call gmap(gvec,syms,pol%nmtx,itran,kgq,pol%isrtx,pol%isrtxi,ind,ph,.true.) 806 pht(:,it,irk) = ph(:) 807 indt(:,it,irk) = ind(:) 808 809! debug Statement here 810 enddo 811 812 if(peinf%inode .eq. 0) then 813 call timing%stop(timing%gmap) 814 endif 815 816 817!--------- loop over occupied states ------------------------------------- 818 819! SIB: loop over valence states (iv,ispin) where iv is the band index. 820 821 if (pol%os_opt_ffts==2) then 822 if (.not.wfn_FFT_comm%done) call genwf_FFT_Wait(wfn_FFT_comm) 823 !FHJ: TODO free me later! 824 call genwf_lvl2(kp,kpq,vwfn,pol,cwfn) 825 endif 826 do iv=1,peinf%nvownmax 827 ! FHJ : friendly output / running time estimate 828 call progress_step(prog_info) 829 830 if (peinf%verb_debug .and. peinf%inode==0) then 831 write(6,*) 'Doing iv', iv,'of', peinf%nvownmax 832 endif 833 834 if (pol%os_opt_ffts/=2) then 835 call genwf_gen(syms,gvec,crys,kp,kpq,irk,rk,qq,vwfn,pol,cwfn,use_WFNq,intwfnv,intwfnvq,intwfnc,iv) 836 837 if (pol%os_opt_ffts==1) then 838 ! FHJ: FFT opt. level 1: precalculates FFTs of the conduction bands 839 ! each kpt at a time. 840 if (iv==1) then 841 call mtxel_init_FFT_cond(gvec,pol,cwfn,kp) 842 endif 843 endif 844 endif 845 846! SIB: compute matrix elements and energy denominators for (iv,ispin) 847! with all other conduction bands. 848 849 do ispin=1,kp%nspin 850 851 write(tmpstr,'(a,i2,a,i4,a)') "is =", ispin, " iv = ", iv, " calling mtxel" 852 call logit(tmpstr) 853 854 if ( iv .le. peinf%nvownactual) then 855 if(peinf%inode .eq. 0) then 856 call timing%start(timing%mtxel) 857 endif 858 ivin=peinf%invindexv(iv) 859 860 kfact = 1d0 861 call mtxel(ivin,gvec,vwfn,cwfn,pol,ispin,irk,kp,kpq,peinf%rank_mtxel,kfact) 862 if(peinf%inode .eq. 0) then 863 call timing%stop(timing%mtxel) 864 endif 865 endif 866 867 enddo ! ispin 868 869 if (pol%os_opt_ffts<2) then 870 ! FHJ: opt. lvl 2 doesn`t even own the WFNs.. 871 if ( iv .le. peinf%nvownactual) then 872 SAFE_DEALLOCATE_P(vwfn%zv) 873 endif 874 endif 875 876 enddo ! iv 877 878 if (peinf%nvownactual>0) then 879 if (pol%os_opt_ffts<2) then 880 SAFE_DEALLOCATE_P(cwfn%zc) 881 endif 882 if (pol%os_opt_ffts==1) then 883 ! FHJ: destroy FFTs of conduction bands 884 call mtxel_free_FFT_cond(cwfn) 885 endif 886 887 SAFE_DEALLOCATE_P(vwfn%ev) 888 SAFE_DEALLOCATE_P(cwfn%ec) 889 endif 890 891 enddo ! irk 892 call progress_free(prog_info) 893#ifdef USEVORO 894 if (pol%non_uniform) then 895 SAFE_DEALLOCATE(kweights) 896 endif 897#endif 898 899!------------------------------------------------------------------ 900! DVF: if requested, test convergence of chi with conduction bands 901 902 if (peinf%inode .eq. 0 .and. pol%freq_dep .eq. 0 .and. pol%fullConvLog .ne. -1) then 903 write(6,'(1x,a,i0,a)') 'Preparing simple convergence tests with ', & 904 cwfn%nband-vwfn%nband, ' unoccupied bands.' 905 endif 906 907 if(peinf%inode .eq. 0) then 908 call timing%start(timing%converge_tests) 909 endif 910 911 if (pol%freq_dep .eq. 0 .and. pol%fullConvLog .ne. -1) then 912 913 call create_chi_converger(chi_converger,vwfn%nband,cwfn%nband) 914 915 if (peinf%verb_debug .and. peinf%inode==0) then 916 write(6,'(/,1x,"Starting Convergence Tests")') 917 endif 918 call chi_convergence_test(pol,pht,indt,kp,nrk,nst,vwfn%nband,cwfn%nband,fact,chi_converger) 919 920 if(peinf%inode .eq. 0) then 921 call chi_convergence_print(pol,iq,vwfn%nband,cwfn%nband,chi_converger) 922 endif 923 924 call free_chi_converger(chi_converger) 925 926! write(6,*) 'End Convergence Writing' 927 928 endif ! pol%freq_dep .eq. 0 929 930 if(peinf%inode .eq. 0) then 931 call timing%stop(timing%converge_tests) 932 endif 933 934!----------------------------------------------------------------------- 935! Construct part of chi that this proc owns by summing the pol%gme matrices 936 937 if (peinf%verb_debug .and. peinf%inode==0) write(6,'(/,1x,"Doing chi Summation")') 938 if(peinf%inode .eq. 0) then 939 call timing%start(timing%chi_sum_total) 940 endif 941 942 call create_chi_summator(chi_summator, pol, scal, fact, kp%nspin) 943 944#if defined MPI && defined USESCALAPACK 945 IF(pol%subspace) THEN 946 ! here we do subspace truncation 947 IF(pol%gcomm == 0) call die('Subspace truncation implemented only for gcomm=0') 948 949 call chi_summation_sub_trunc(chi_summator,pol,scal,kp,kpq,vwfn,cwfn,& 950 nst,nrk,indt,pht,gvec,crys,& 951 qq,iq) 952 ! copy into right location 953 pol%neigen_of_q(iq) = pol%neig_sub 954 ELSE 955#endif 956 957 if (pol%gcomm .eq. 0) then 958 call chi_summation_comm_elements(chi_summator,& 959 pol,scal,kp,vwfn,cwfn,& 960 nst,nrk,indt,pht) 961 else 962 call chi_summation_comm_matrix(chi_summator,& 963 pol,scal,kp,kpq,vwfn,& 964 nst,nrk,indt,pht) 965 endif 966 967#if defined MPI && defined USESCALAPACK 968 END IF 969#endif 970 971 call free_chi_summator(chi_summator, pol) 972 973 if(peinf%inode .eq. 0) then 974 call timing%stop(timing%chi_sum_total) 975 endif 976 977 if (peinf%verb_debug .and. peinf%inode==0) write(6,'(1x,a)') "Done Polarizability" 978 979! Done ChiSum 980!----------------------------------------------------------------------- 981! Deallocate some arrays no longer needed 982 983 SAFE_DEALLOCATE(pht) 984 SAFE_DEALLOCATE(indt) 985 SAFE_DEALLOCATE(ind) 986 SAFE_DEALLOCATE(ph) 987 SAFE_DEALLOCATE(nst) 988 SAFE_DEALLOCATE_P(pol%gme) 989 990 if (pol%freq_dep .eq. 2 .or. pol%freq_dep .eq. 3) then 991 SAFE_DEALLOCATE_P(pol%edenDyn) 992 endif 993 994 else ! pol%skip_chi 995!DVF: read chi from chi if this is specified 996 if (peinf%inode==0) write(6,'(/1x,a)') 'Reading polarizability matrix from file' 997#ifdef HDF5 998 if (pol%use_hdf5) then 999 if (peinf%inode .eq. 0) then 1000 ! FHJ: FIXME: consistency checks... 1001 endif 1002 ! FHJ: TODO - write diagonal elements (I actually think this is useless) 1003 if (pol%freq_dep .eq. 0) then 1004 do ispin=1,kp%nspin 1005 call read_matrix_d_hdf5(scal, pol%chi(:,:,ispin), pol%nmtx, TRUNC(filename_chi_hdf5), my_iq, ispin) 1006 enddo 1007 else 1008 do ispin=1,kp%nspin 1009 call read_matrix_f_hdf5(scal, pol%nFreq, pol%nfreq_in_group, pol%chiRDyn(:,:,:,ispin), & 1010 pol%nmtx, pol%nfreq_group, TRUNC(filename_chi_hdf5), my_iq, ispin) 1011 enddo 1012 endif 1013 else 1014#endif 1015 1016 if (is_q0) then 1017 1018 iunit=10 1019 if(peinf%inode.eq.0) then 1020 write(6,'(1x,a)') 'Reading from file chi0mat.' 1021 call open_file(unit=10,file='chi0mat',form='unformatted',status='old') 1022 read(10) 1023 read(10) 1024 read(10) 1025 read(10) 1026 read(10) 1027 read(10) 1028 read(10) 1029 read(10) 1030 read(10) 1031 read(10) 1032 read(10) 1033 read(10) nmtx_t 1034 endif 1035#ifdef MPI 1036 call mpi_barrier(MPI_COMM_WORLD,mpierr) 1037#endif 1038 do ispin=1,kp%nspin 1039 if (pol%freq_dep .eq. 0) then 1040 call read_matrix_d(scal,pol%chi(:,:,ispin),pol%nmtx,iunit) 1041 endif ! pol%freq_dep .eq. 0 1042 if (pol%freq_dep .eq. 2) then 1043 call read_matrix_f(scal, pol%nFreq, pol%nfreq_in_group, & 1044 pol%chiRDyn(:,:,:,ispin), pol%nmtx, pol%nfreq_group, iunit) 1045 endif ! pol%freq_dep .eq. 2! 1046 enddo ! ispin 1047 else ! is_q0 1048 1049 iunit=11 1050 if(peinf%inode.eq.0) then 1051 write(6,'(1x,a)') 'Reading from file chimat.' 1052 if (initial_access .eq. 0) then 1053 call open_file(unit=11,file='chimat',form='unformatted',status='old') 1054 read(11) 1055 read(11) 1056 read(11) 1057 read(11) 1058 read(11) 1059 read(11) 1060 read(11) 1061 read(11) 1062 read(11) 1063 read(11) 1064 endif 1065 read(11) 1066 read(11) nmtx_t 1067! write(6,*) 'nmtx_t for chimat', nmtx_t 1068 endif 1069#ifdef MPI 1070 call mpi_barrier(MPI_COMM_WORLD,mpierr) 1071#endif 1072 do ispin=1,kp%nspin 1073 if (pol%freq_dep .eq. 0) then 1074 call read_matrix_d(scal,pol%chi(:,:,ispin),pol%nmtx,iunit) 1075 endif ! pol%freq_dep .eq. 0 1076 if (pol%freq_dep .eq. 2) then 1077 call read_matrix_f(scal, pol%nFreq, pol%nfreq_in_group, & 1078 pol%chiRDyn(:,:,:,ispin), pol%nmtx, pol%nfreq_group, iunit) 1079 endif ! pol%freq_dep .eq. 2 1080 enddo ! ispin 1081 initial_access = 1 1082 endif ! is_q0 1083#ifdef HDF5 1084 endif 1085#endif 1086 if (peinf%inode==0) write(6,'(1x,a/)') 'Ok' 1087 endif ! pol%skip_chi 1088 1089!----------------------------------------------------------------------- 1090! JRD: Now write out elements that Proc 1 owns 1091 1092 do ispin = 1, kp%nspin 1093 1094 if (pol%freq_dep.eq.0 .and. peinf%inode.eq.0) then 1095 write(7,940) ispin,kp%nspin 1096 do i=1,scal%npr 1097 ix=scal%isrtxrow(i) 1098 do j=1,scal%npc 1099 1100! JRD: Diagonal, subdiagonal and wings only 1101 1102 jx=scal%isrtxcol(j) 1103 if(i.eq.j .or. (gvec%components(1,ix) .eq. 0 .and. gvec%components(2,ix) .eq. 0 .and. gvec%components(3,ix) .eq. 0)) & 1104 write(7,950) (gvec%components(k,ix),k=1,3),ekin(ix),(gvec%components(k,jx),k=1,3),ekin(jx), & 1105 pol%chi(i,j,ispin) 1106 enddo 1107 enddo 1108 endif ! pol%freq_dep.eq.0 .and. peinf%inode.eq.0 1109 1110 if ((pol%freq_dep .eq. 2).and.(pol%freq_dep_method .eq. 0 .or. pol%freq_dep_method .eq. 2) .and. peinf%inode.eq.0) then 1111 1112 IF(pol%subspace .AND. (.NOT.pol%need_full_chi) ) THEN 1113 ! write only for omega = 0 1114 ! write(7,940) ispin,kp%nspin 1115 IF(ispin == 1) THEN 1116 write(7,*)'frq=',0 1117 do i=1,scal%npr 1118 ix=scal%isrtxrow(i) 1119 do j=1,scal%npc 1120 jx=scal%isrtxcol(j) 1121 if(i.eq.j) & 1122 write(7,950) (gvec%components(k,ix),k=1,3),ekin(ix),(gvec%components(k,jx),k=1,3),ekin(jx), & 1123 pol%chiRDyn_sym_omega0(i,j) 1124 enddo 1125 enddo 1126 END IF 1127 WRITE(7,*) 1128 WRITE(7,*) 'Eigenvalues symmetrized chi at omega = 0' 1129 DO i = 1, pol%nmtx 1130 WRITE(7,*) i, pol%eigenval_omega0(i) 1131 END DO 1132 WRITE(7,*) 1133 ELSE 1134 1135 write(7,940) ispin,kp%nspin 1136 1137 do jj=1,pol%nfreq_in_group 1138 write(7,*)'frq=',jj 1139! JRD XXX the i and j loops are out of order here.... 1140 do i=1,scal%npr 1141 ix=scal%isrtxrow(i) 1142 do j=1,scal%npc 1143 1144 !!!JRD: Diagonal and subdiagonal only 1145 1146 jx=scal%isrtxcol(j) 1147 if(i.eq.j) & 1148 write(7,950) (gvec%components(k,ix),k=1,3),ekin(ix),(gvec%components(k,jx),k=1,3),ekin(jx), & 1149 pol%chiRDyn(i,j,jj,ispin) 1150 enddo 1151 enddo 1152 enddo 1153 END IF ! IF(pol%subspace 1154 1155 endif ! (pol%freq_dep .eq. 2).and.(pol%freq_dep_method .eq. 0 .or. pol%freq_dep_method .eq. 2) .and. peinf%inode.eq.0 1156 1157 if ((pol%freq_dep .eq. 2).and.(pol%freq_dep_method .eq. 1).and. peinf%inode.eq.0) then 1158 write(7,940) ispin,kp%nspin 1159 1160 do jj=1,pol%nfreq_in_group 1161 write(7,*) 'frq=',jj 1162! JRD XXX the i and j loops are out of order here 1163 do i=1,scal%npr 1164 ix=scal%isrtxrow(i) 1165 do j=1,scal%npc 1166 jx=scal%isrtxcol(j) 1167 if(i.eq.j) & 1168 write(7,950) (gvec%components(k,ix),k=1,3),ekin(ix),(gvec%components(k,jx),k=1,3),ekin(jx), & 1169 pol%chiRDyn(i,j,jj,ispin) 1170 enddo 1171 enddo 1172 enddo 1173 1174 endif ! (pol%freq_dep .eq. 2).and.(pol%freq_dep_method .eq. 1) .and. peinf%inode.eq.0 1175 1176940 format(/,10x,' independent matrix elements of chi', 7x,'spin index= ',1i1,1x,1i1,/,/,& 1177 10x,'g',10x,'g**2',10x,'gp',10x,'gp**2',10x,'chi(g,gp)') 1178! if last value is real, only one of the f13.8 will be used. 1179950 format(3i5,f10.5,5x,3i5,f10.5,5x,2f15.10) 1180 1181 enddo ! ispin (loop over spins) 1182 1183! write(6,*) 'End Element Writing' 1184 1185 1186!--------- write polarizability matrix and crystal info to file --------- 1187 1188 if (pol%skip_epsilon) then 1189 if (peinf%inode==0) write(6,'(/1x,a)') 'Writing polarizability matrix to file' 1190#ifdef HDF5 1191 if (pol%use_hdf5) then 1192 if (peinf%inode .eq. 0) then 1193 SAFE_ALLOCATE(isorti, (gvec%ng)) 1194 do i=1,gvec%ng 1195 isorti(pol%isrtx(i)) = i 1196 enddo 1197 call write_gvec_indices_hdf(gvec%ng,pol%isrtx,isorti,ekin,my_iq,TRUNC(filename_chi_hdf5)) 1198 SAFE_DEALLOCATE(isorti) 1199 endif 1200 ! FHJ: TODO - write diagonal elements (I actually think this is useless) 1201 if (pol%freq_dep .eq. 0) then 1202 do ispin=1,kp%nspin 1203 ! FHJ: FIXME: We should unify these routines! 1204#ifdef USESCALAPACK 1205 call write_matrix_d_par_hdf(scal, pol%chi(:,:,ispin), pol%nmtx, my_iq, ispin, TRUNC(filename_chi_hdf5)) 1206#else 1207 call write_matrix_d_hdf(scal, pol%chi(:,:,ispin), pol%nmtx, my_iq, ispin, TRUNC(filename_chi_hdf5)) 1208#endif 1209 enddo 1210 else 1211 do ispin=1,kp%nspin 1212 ! FHJ: FIXME: We should unify these routines! 1213#ifdef USESCALAPACK 1214 call write_matrix_f_par_hdf(scal, pol%nFreq_in_group, pol%chiRDyn(:,:,:,ispin), & 1215 pol%nmtx, my_iq, ispin, TRUNC(filename_chi_hdf5), pol%nfreq_group) 1216#else 1217 call write_matrix_f_hdf(scal, pol%nFreq, pol%chiRDyn(:,:,:,ispin), & 1218 pol%nmtx, my_iq, ispin, TRUNC(filename_chi_hdf5)) 1219#endif 1220 enddo 1221 endif 1222 else 1223#endif 1224 iunit=11 1225 if (is_q0) iunit=10 1226 if(peinf%inode.eq.0) then 1227 write(iunit) syms%ntranq,(((syms%mtrx(i,j,syms%indsub(n)),i=1,3),j=1,3), & 1228 (syms%tnp(k,syms%indsub(n)),syms%kgzero(k,n),k=1,3),n=1,syms%ntranq) 1229 np=pol%nmtx*(pol%nmtx+1)/2 1230 write(iunit) pol%nmtx,np,(pol%isrtx(i),ekin(i),i=1,gvec%ng),(pol%irow(i),i=1,pol%nmtx) 1231 endif 1232 do ispin=1,kp%nspin 1233 if (pol%freq_dep .eq. 0) then 1234 call write_matrix_d(scal,pol%chi(:,:,ispin),pol%nmtx,iunit) 1235 endif ! pol%freq_dep .eq. 0 1236 if (pol%freq_dep==2 .or. pol%freq_dep==3) then 1237 call write_matrix_f(scal, pol%nFreq, pol%chiRDyn(:,:,:,ispin), & 1238 pol%nmtx, iunit, pol%nfreq_group) 1239 endif 1240 enddo ! ispin 1241#ifdef HDF5 1242 endif 1243#endif 1244 if (peinf%inode==0) write(6,'(1x,a/)') 'Ok' 1245 endif ! pol%skip_epsilon 1246 1247! Use pol%chi(j,1) as sum over spin components 1248! JRD: Why was proc 0 the only one doing this?! 1249 1250 if (pol%freq_dep .eq. 0) then 1251 if(kp%nspin.eq.2) pol%chi(:,:,1)=pol%chi(:,:,1)+pol%chi(:,:,2) 1252 endif ! pol%freq_dep .eq. 0 1253 if (pol%freq_dep .eq. 2) then 1254 if(kp%nspin.eq.2) pol%chiRDyn(:,:,:,1)=pol%chiRDyn(:,:,:,1)+pol%chiRDyn(:,:,:,2) 1255 endif ! pol%freq_dep .eq. 2 1256 1257 if (peinf%inode .eq. 0) then 1258 1259 if(peinf%inode .eq. 0) then 1260 call timing%start(timing%epsinv_total) 1261 endif 1262 1263 call logit('Calling epsinv') 1264 endif ! peinf%inode .eq. 0 1265 1266 if (.not. pol%skip_epsilon) then 1267 if(pol%do_rpa) then 1268 call epsinv(gvec,pol,ekin,qq,is_q0,crys,scal,kp,omega_plasma,iq,E_rpa) 1269 pol%E_rpa_qp(iq) = E_rpa 1270 else 1271 call epsinv(gvec,pol,ekin,qq,is_q0,crys,scal,kp,omega_plasma,iq) 1272 endif 1273 1274 IF(pol%subspace .AND. (.NOT.pol%need_full_chi)) THEN 1275 SAFE_DEALLOCATE(pol%chiRDyn_sym_omega0) 1276 SAFE_DEALLOCATE(pol%eigenvect_omega0) 1277 SAFE_DEALLOCATE(pol%eigenval_omega0) 1278 SAFE_DEALLOCATE(pol%vcoul_sub) 1279 END IF 1280 1281 if(peinf%inode.eq.0) then 1282 call logit('Finished epsinv') 1283 if(peinf%inode .eq. 0) then 1284 call timing%start(timing%epsinv_total) 1285 endif 1286 endif ! peinf%inode.eq.0 1287 1288 endif ! pol%skip_epsilon 1289 1290 if (pol%freq_dep .eq. 0) then 1291 SAFE_DEALLOCATE_P(pol%chi) 1292 endif 1293 if (pol%freq_dep .eq. 2) then 1294 SAFE_DEALLOCATE_P(pol%chiRDyn) 1295 if(pol%freq_dep_method .eq. 1) then 1296 SAFE_DEALLOCATE_P(pol%chiTDyn) 1297 endif 1298 endif 1299 if (pol%freq_dep .eq. 3) then 1300 SAFE_DEALLOCATE_P(pol%chiRDyn) 1301 endif 1302 1303 1304 SAFE_DEALLOCATE(indrk) 1305 SAFE_DEALLOCATE(neq) 1306 1307 SAFE_DEALLOCATE_P(pol%isrtx) 1308 SAFE_DEALLOCATE_P(pol%isrtxi) 1309 SAFE_DEALLOCATE_P(pol%irow) 1310 SAFE_DEALLOCATE_P(scal%isrtxcol) 1311 SAFE_DEALLOCATE_P(scal%isrtxrow) 1312 SAFE_DEALLOCATE_P(scal%imycol) 1313 SAFE_DEALLOCATE_P(scal%imyrow) 1314 SAFE_DEALLOCATE_P(scal%imycolinv) 1315 SAFE_DEALLOCATE_P(scal%imyrowinv) 1316 SAFE_DEALLOCATE_P(scal%imycold) 1317 SAFE_DEALLOCATE_P(scal%imyrowd) 1318 1319#ifdef HDF5 1320 if (pol%use_hdf5 .and. peinf%inode .eq. 0) then 1321 if (.not. pol%skip_epsilon) then 1322 if ( pol%subspace ) then 1323 if ( pol%matrix_in_subspace_basis ) then 1324 call set_subspace_neigen_q(TRUNC(filename_out_hdf5), my_iq, pol%neig_sub ) 1325 end if 1326 end if 1327 end if 1328 end if 1329#endif 1330 1331#ifdef HDF5 1332 if (pol%use_hdf5.and.peinf%inode==0) call set_qpt_done(TRUNC(filename_out_hdf5), my_iq) 1333#endif 1334 1335 enddo ! iq (loop over q points) 1336 1337 SAFE_DEALLOCATE_P(scal%nprd) 1338 SAFE_DEALLOCATE_P(scal%npcd) 1339 call dealloc_grid(gr) 1340 1341! End q point loop! 1342!------------------------------------------------------------------- 1343 1344! XXXXXXXXXXXXXXXX Do BZ sum to get RPA correlation energy 1345 if(pol%do_rpa) then 1346 E_rpa = 0.0D+00 1347 do iq = 1, pol%nq 1348 E_rpa = E_rpa + pol%qw_rpa(iq) * pol%E_rpa_qp(iq) 1349 enddo 1350 if(peinf%inode.eq.0) then 1351 write(*,*) "RPA energy (Ry) = ", E_rpa 1352 write(*,*) "RPA energy (Ha) = ", E_rpa/2.0D+00 1353 write(*,*) "RPA energy (eV) = ", E_rpa*ryd 1354 endif 1355 endif 1356! XXXXXXXXXXXXXXXX 1357 1358!----------- Clean House ------------------------------------------- 1359 1360 call logit('Cleaning up') 1361 if (.not. pol%skip_epsilon) then 1362 call destroy_qran() 1363 endif 1364 1365 if(.not. pol%skip_chi) then 1366 if(peinf%inode == 0) call close_file(17) ! file chi_converge.dat 1367 call destroy_fftw_plans() 1368 endif 1369 if (pol%iwritecoul .eq. 1) then 1370 if (peinf%inode .eq. 0) then 1371 call close_file(19) ! file vcoul 1372 endif 1373 endif 1374 1375 if (peinf%inode==0 .and. pol%freq_dep==2 .and. .not.pol%skip_epsilon) then 1376 call close_file(51) !file EpsInvDyn 1377 call close_file(52) !file EpsDyn 1378 endif 1379 1380#ifdef USEVORO 1381 if (pol%non_uniform) then 1382 SAFE_DEALLOCATE(kvols) 1383 endif 1384#endif 1385 SAFE_DEALLOCATE(ekin) 1386 SAFE_DEALLOCATE_P(kp%w) 1387 SAFE_DEALLOCATE_P(kp%rk) 1388 SAFE_DEALLOCATE_P(kp%el) 1389 1390 if(pol%need_WFNq) then 1391 SAFE_DEALLOCATE_P(kpq%w) 1392 SAFE_DEALLOCATE_P(kpq%rk) 1393 SAFE_DEALLOCATE_P(kpq%el) 1394 endif 1395 SAFE_DEALLOCATE_P(gvec%components) 1396 SAFE_DEALLOCATE_P(gvec%index_vec) 1397 SAFE_DEALLOCATE_P(pol%qpt) 1398 SAFE_DEALLOCATE_P(vwfn%isort) 1399 SAFE_DEALLOCATE_P(cwfn%isort) 1400 SAFE_DEALLOCATE_P(peinf%global_nvown) 1401 SAFE_DEALLOCATE_P(peinf%global_ncown) 1402 SAFE_DEALLOCATE_P(peinf%indexc) 1403 SAFE_DEALLOCATE_P(peinf%indexv) 1404 SAFE_DEALLOCATE_P(peinf%global_indexv) 1405 SAFE_DEALLOCATE_P(peinf%invindexv) 1406 SAFE_DEALLOCATE_P(peinf%invindexc) 1407 SAFE_DEALLOCATE_P(peinf%doiownv) 1408 SAFE_DEALLOCATE_P(peinf%doiownc) 1409 SAFE_DEALLOCATE_P(peinf%does_it_ownv) 1410 SAFE_DEALLOCATE_P(peinf%does_it_ownc) 1411 SAFE_DEALLOCATE_P(peinf%global_pairowner) 1412 SAFE_DEALLOCATE_P(pol%dFreqGrid) 1413 SAFE_DEALLOCATE_P(pol%dFreqBrd) 1414 if(peinf%inode.eq.0) then 1415 SAFE_DEALLOCATE_P(pol%nmtx_of_q) 1416 call close_file(7) ! epsilon.log 1417 1418 if (pol%skip_epsilon.and..not.pol%use_hdf5) then 1419 if (pol%nq0>0) call close_file(10) ! chi0mat 1420 if (pol%nq1>0) call close_file(11) ! chimat 1421 else 1422 if (.not.pol%use_hdf5) then 1423 if (pol%nq0>0) call close_file(12) ! eps0mat 1424 if (pol%nq1>0) call close_file(13) ! epsmat 1425 endif 1426 endif ! pol%skip_epsilon 1427 endif 1428 if(pol%subspace) then 1429 SAFE_DEALLOCATE(pol%neigen_of_q) 1430 end if 1431 1432 call free_wfns(pol, intwfnv, intwfnvq, intwfnc, .true.) 1433 1434!------------- Print Timing Info ----------------------------------------- 1435 1436#ifdef MPI 1437 call MPI_barrier(MPI_COMM_WORLD,mpierr) 1438#endif 1439 call logit('Calculating Timing Info') 1440 1441 call timing%stop(timing%total) 1442 call timing%print(common_timing, .true.) 1443 1444 SAFE_DEALLOCATE(routsrt) 1445 1446 call write_memory_usage() 1447 1448!------------------------------- 1449! JIM: Close HDF interface 1450#ifdef HDF5 1451 if(pol%use_hdf5.or.pol%os_hdf5) call h5close_f(hdf5_error) 1452#endif 1453 1454#ifdef MPI 1455 call MPI_Finalize(mpierr) 1456#endif 1457 1458end program epsilon 1459