1!============================================================================
2!
3! Routines:
4!
5! (1) mtxel_occ()        Originally By ?         Last Modified 8/15/2015 (FHJ)
6!
7!     Adapted from subroutine mtxel. Subroutine computes required matrix
8!     elements that involve only occupied states, of the form <nk|exp{i g.r}|nk>
9!     They are required for COHSEX and exact CH calculations.
10!
11!     input   n,m                band indices
12!     input   ncoul              number of matrix elements required
13!     input   isrtrq             index array for g-vectors in
14!                                <nk|exp{i g.r}|nk>
15!     output  aqs                matrix elements required
16!
17!============================================================================
18
19#include "f_defs.h"
20
21subroutine mtxel_occ(n,m,gvec,wfnk,ncoul,isrtrq,aqs,ispin,kp)
22
23  use global_m
24  use fftw_m
25  use misc_m
26  implicit none
27
28  integer, intent(in) :: n, m
29  type (gspace), intent(in) :: gvec
30  type (wfnkstates), intent(in) :: wfnk
31  type (kpoints), intent(inout) :: kp
32  integer, intent(in) :: ncoul
33  integer, intent(in) :: isrtrq(gvec%ng)
34  SCALAR, intent(out) :: aqs(ncoul)
35  integer, intent(in) :: ispin
36
37!-------------------------
38! If we are using FFT to calculate matrix elements...
39
40! We use FFT to compute <u_nk|e^(iG.r)|u_nk> elements where
41! u_nk is the periodic part of the wave function.
42! The calculation is done in real space, and integration over
43! the grid is replaced by the sum over the grid points p:
44!
45! <u_nk|e^(iG.r)|u_nk>  =
46!     Volume/Np * sum_p { conj(u_nk(p))*e^(iG.p)*u_nk(p) }
47!
48! Since u_nk(p) = Volume^-0.5 * sum_G { cnk(G)*e^(iG.p) },
49! and FFT is defined as FFT(cnk,+,p) = sum_G { cnk(G)*e^{+iG.p} },
50! we must compute
51!
52! <u_nk|e^(iG.r)|u_nk>
53!   = 1/Np * sum_p { conj(FFT(cnk,+,p))*e^(iG.p)*FFT(cnk,+,p) }
54!   = 1/Np * FFT(conj(FFT(cnk,+,:)).*FFT(cnk,+,:),+,G)
55!
56! where .* is a point by point multiplication on the grid
57
58  complex(DPC), dimension(:,:,:), allocatable :: fftbox1,fftbox2
59  integer :: jsp
60  integer, dimension(3) :: Nfft
61  real(DP) :: scale
62  SCALAR, dimension(:), allocatable :: tmparray
63
64  PUSH_SUB(mtxel_occ)
65
66! Compute size of FFT box we need and scale factor
67
68  call setup_FFT_sizes(gvec%FFTgrid,Nfft,scale)
69
70! Allocate FFT boxes
71
72  SAFE_ALLOCATE(fftbox1, (Nfft(1),Nfft(2),Nfft(3)))
73  SAFE_ALLOCATE(fftbox2, (Nfft(1),Nfft(2),Nfft(3)))
74
75! Put the data for band n into FFT box 1 and do the FFT,zk(:,1)
76
77  SAFE_ALLOCATE(tmparray, (ncoul))
78
79  do jsp = ispin,ispin*kp%nspinor
80
81    call put_into_fftbox(wfnk%nkpt,wfnk%zk((n-1)*wfnk%nkpt+1:,jsp),gvec%components,wfnk%isrtk,fftbox1,Nfft)
82    call do_FFT(fftbox1,Nfft,1)
83
84! We need the complex conjugate of the |nk> band actually
85
86    call conjg_fftbox(fftbox1,Nfft)
87
88! Now we get the matrix elements:
89!  Get n wave function and put it into box 2,
90!  do FFT,
91!  multiply by box1 contents,
92!  do FFT again,
93!  and extract the resulting matrix elements
94
95    call put_into_fftbox(wfnk%nkpt,wfnk%zk((m-1)*wfnk%nkpt+1:,jsp),gvec%components,wfnk%isrtk,fftbox2,Nfft)
96    call do_FFT(fftbox2,Nfft,1)
97    call multiply_fftboxes(fftbox1,fftbox2,Nfft)
98    call do_FFT(fftbox2,Nfft,1)
99    call get_from_fftbox(ncoul,tmparray,gvec%components,isrtrq,fftbox2,Nfft,scale)
100    if (kp%nspinor.eq.1 .or. jsp.eq. 1) then
101      aqs(:) = tmparray(:)
102    else
103      aqs(:) = aqs(:) + tmparray(:)
104    endif
105
106  enddo
107
108  SAFE_DEALLOCATE(tmparray)
109
110! We are done, so deallocate FFT boxes
111
112  SAFE_DEALLOCATE(fftbox1)
113  SAFE_DEALLOCATE(fftbox2)
114
115  POP_SUB(mtxel_occ)
116
117  return
118end subroutine mtxel_occ
119