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