1!==============================================================================
2!
3! Modules:
4!
5! chi_convergence_m   (Originally by DVF)            Last Modified: 05/08/2015 (DVF)
6!
7!   This module has functions to create the structures for the convergence tests
8!   with respect to to bands in epsilon, to do the convergence tests, and to
9!   print the results.
10!
11!==============================================================================
12
13#include "f_defs.h"
14
15module chi_convergence_m
16
17  use global_m
18
19  implicit none
20
21  private
22
23  public :: create_chi_converger, free_chi_converger,&
24    chi_convergence_test,chi_convergence_print, &
25    chi_converger_t
26
27  type chi_converger_t
28    SCALAR, allocatable :: head(:)
29    SCALAR, allocatable :: maxG(:)
30    SCALAR, allocatable :: head2(:)
31    SCALAR, allocatable :: maxG2(:)
32    SCALAR, allocatable :: headtotal(:)
33    SCALAR, allocatable :: maxGtotal(:)
34  end type chi_converger_t
35
36contains
37
38  subroutine create_chi_converger(conv,nvalbands,ntotbands)
39    type(chi_converger_t), intent(inout) :: conv !< the chi_converger_t object
40    integer, intent(in) :: nvalbands  !< number of valence bands
41    integer, intent(in) :: ntotbands  !< total number of bands
42
43    PUSH_SUB(create_chi_converger)
44
45    SAFE_ALLOCATE(conv%head,(ntotbands-nvalbands))
46    SAFE_ALLOCATE(conv%maxG,(ntotbands-nvalbands))
47    SAFE_ALLOCATE(conv%head2,(ntotbands-nvalbands))
48    SAFE_ALLOCATE(conv%maxG2,(ntotbands-nvalbands))
49    SAFE_ALLOCATE(conv%headtotal,(ntotbands-nvalbands))
50    SAFE_ALLOCATE(conv%maxGtotal,(ntotbands-nvalbands))
51
52    conv%head = ZERO
53    conv%head2 = ZERO
54    conv%headtotal = ZERO
55    conv%maxG = ZERO
56    conv%maxG2 = ZERO
57    conv%maxGtotal = ZERO
58
59    POP_SUB(create_chi_converger)
60    return
61
62  end subroutine create_chi_converger
63
64  subroutine free_chi_converger(conv)
65    type(chi_converger_t), intent(inout) :: conv !<the chi_converger_t object
66
67    PUSH_SUB(free_chi_converger)
68
69    SAFE_DEALLOCATE(conv%head)
70    SAFE_DEALLOCATE(conv%maxG)
71    SAFE_DEALLOCATE(conv%head2)
72    SAFE_DEALLOCATE(conv%maxG2)
73    SAFE_DEALLOCATE(conv%headtotal)
74    SAFE_DEALLOCATE(conv%maxGtotal)
75
76    POP_SUB(free_chi_converger)
77    return
78
79  end subroutine free_chi_converger
80
81
82  subroutine chi_convergence_test(pol,pht,indt,kp,nrk,nst,nvalbands,ntotbands,fact,conv)
83    type(polarizability), intent(in) :: pol !< polarizability
84    SCALAR, intent(in) :: pht(:, :, :) !< (pol%nmtx,neqmax,nrk) phases
85      !! neqmax = max number of kpoints equivalent by symmetry (across all irr. k-points)
86    integer, intent(in) :: indt(:, :, :) !< (pol%nmtx,neqmax,nrk) indexing array
87    type(kpoints),intent(in) :: kp !<kpoint array
88    integer, intent(in) :: nrk !<number of kpoints
89    integer, intent(in) :: nst(:) !< (nrk) !<number of elements in the star
90    integer, intent(in) :: nvalbands !< number of valence bands
91    integer, intent(in) :: ntotbands !< number of total bands
92    real(DP), intent(in) :: fact !< volume factor used throughout the code
93    type(chi_converger_t), intent(inout) :: conv
94
95    SCALAR :: gmehead,gmemaxG  !<matrix elements for head and tail (maxG)
96    SCALAR :: mod_square_phthead,mod_square_phtmaxG !<mod squared of phase
97    integer :: ispin,irk,iv,ic,istar,isend,iown
98
99    ! JRD/DVF: Test convergence with bands for G=G`=0 and G=G`=pol%nmtx
100    ! The result is chi(0,0), chi(Gmax,Gmax) as a function of conduction
101    ! bands, i.e. for these gvectors the chi summation is done for all
102    ! valence bands and k-points for certain numbers of conduction bands
103    ! and the convergence is tracked as the number of conduction bands
104    ! is increased. These will be printed in chi_convergence_print.
105
106    PUSH_SUB(chi_convergence_test)
107
108    do ispin =1, kp%nspin
109      do irk = 1, nrk
110        do istar = 1, nst(irk)
111
112          mod_square_phthead=pht(1,istar,irk)*MYCONJG(pht(1,istar,irk))
113          mod_square_phtmaxG=pht(pol%nmtx,istar,irk)*MYCONJG(pht(pol%nmtx,istar,irk))
114
115          do iv=1,(nvalbands+pol%ncrit)
116            iown =1
117            do ic=1,ntotbands-nvalbands
118              isend=peinf%global_pairowner(iv,ic)-1
119              if (isend .eq. peinf%inode) then
120                if (iown .gt. peinf%ncownactual) write(6,*) 'iown bigger than ncownactual'
121                gmehead = pol%gme(indt(1,istar,irk),iown,peinf%indexv(iv), &
122                  ispin,irk,pol%nfreq_group)
123                gmemaxG = pol%gme(indt(pol%nmtx,istar,irk),iown,peinf%indexv(iv), &
124                  ispin,irk,pol%nfreq_group)
125                conv%head2(ic) = conv%head2(ic) + gmehead*MYCONJG(gmehead)*mod_square_phthead*fact*(-1D0)
126                conv%maxG2(ic) = conv%maxG2(ic) + gmemaxG*MYCONJG(gmemaxG)*mod_square_phtmaxG*fact*(-1D0)
127                iown=iown+1
128              endif
129            enddo ! ic
130          enddo ! iv
131
132        enddo ! istar
133      enddo ! irk
134    enddo ! ispin
135
136#ifdef MPI
137    call MPI_reduce(conv%head2,conv%head,ntotbands-nvalbands,MPI_SCALAR,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
138    call MPI_reduce(conv%maxG2,conv%maxG,ntotbands-nvalbands,MPI_SCALAR,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
139#else
140    conv%head = conv%head2
141    conv%maxG = conv%maxG2
142#endif
143
144    POP_SUB(chi_convergence_test)
145
146  end subroutine chi_convergence_test
147
148  subroutine chi_convergence_print(pol,iq,nvalbands,ntotbands,conv)
149
150    type(polarizability), intent(in) :: pol !<polarizability
151    integer, intent(in) :: iq !<current qpoint
152    integer, intent(in) :: nvalbands !<number of valence bands
153    integer, intent(in) :: ntotbands !<total number of bands
154    type(chi_converger_t), intent(inout) :: conv
155
156    SCALAR :: extrap_head,extrap_maxG
157    integer :: ic,idis,ncondbands
158    character*24 :: strhead,strtail
159
160    ! JRD: Print out convergence
161
162    PUSH_SUB(chi_convergence_print)
163
164    conv%headtotal(1)=conv%head(1)
165    conv%maxGtotal(1)=conv%maxG(1)
166    ncondbands = ntotbands-nvalbands
167    do ic = 2, ncondbands
168      conv%headtotal(ic)=conv%headtotal(ic-1)+conv%head(ic)
169      conv%maxGtotal(ic)=conv%maxGtotal(ic-1)+conv%maxG(ic)
170    enddo
171
172    ! DVF: Get the value of the matrix elements extrapolated with respect to bands
173    ! idis is some sort of extrapolation "distance"
174
175    idis = (ncondbands)/10
176    if (idis.gt.0) then
177      extrap_head = (ncondbands*conv%headtotal(ncondbands) - (ncondbands-idis)*conv%headtotal(ncondbands-idis)) / idis
178      extrap_maxG = (ncondbands*conv%maxGtotal(ncondbands) - (ncondbands-idis)*conv%maxGtotal(ncondbands-idis)) / idis
179    else
180      extrap_head = ZERO
181      extrap_maxG = ZERO
182    endif
183
184    write(strhead,701)1
185    write(strtail,701)pol%nmtx
186
187    if (pol%fullConvLog .eq. 0) then
188      write(17,'(a,3e16.8,a,i0,a)') '# q=', pol%qpt(:,iq), '  (qpt ', iq, ')'
189      write(17,'(a1,a7,4a20)') '#', 'ncbands', 'Re chi(0,0)', 'extrap', 'Re chi(Gmax,Gmax)', 'extrap'
190
191      do ic = 1, ncondbands
192        write(17,'(i8,4e20.8)') ic, dble(conv%headtotal(ic)), &
193          dble(extrap_head), dble(conv%maxGtotal(ic)), dble(extrap_maxG)
194      enddo
195      write(17,*)
196    elseif (pol%fullConvLog .eq. 1) then
197      write(17,801) pol%qpt(:,iq), iq
198      write(17,802)
199
200      write(17,803)TRUNC(strhead)
201      do ic = 1, ncondbands
202        write(17,805) ic, dble(conv%headtotal(ic))
203      enddo
204
205      write(17,804)TRUNC(strtail)
206      do ic = 1, ncondbands
207        write(17,805) ic, dble(conv%maxGtotal(ic))
208      enddo
209    elseif (pol%fullConvLog .eq. 2) then
210      write(17,801) pol%qpt(:,iq), iq
211      write(17,802)
212
213      write(17,803)TRUNC(strhead)
214      do ic = 1, ncondbands
215        write(17,805) ic, conv%headtotal(ic)
216      enddo
217
218      write(17,804)TRUNC(strtail)
219      do ic = 1, ncondbands
220        write(17,805) ic, conv%maxGtotal(ic)
221      enddo
222    endif
223
224701 format(i16)
225801 format('#',1x,'q =',3f10.6,1x,'iq =',i4)
226802 format(2x,'nbands',1x,'chi0')
227803 format(2x,'head ig = igp =',1x,a)
228804 format(2x,'tail ig = igp =',1x,a)
229805 format(i8,2e16.8)
230
231    SAFE_DEALLOCATE(conv%headtotal)
232    SAFE_DEALLOCATE(conv%maxGtotal)
233
234    POP_SUB(chi_convergence_print)
235
236  end subroutine chi_convergence_print
237
238end module chi_convergence_m
239