1subroutine direct_v_exc(a_in,fc,a_out)
2! computes the v part of the direct term of the exc Hamiltonian
3USE fft_custom_gwl
4use bse_basic_structures
5use exciton
6USE wavefunctions, ONLY :  psic
7USE gvect,                 ONLY : ig_l2g
8USE io_global, ONLY : stdout, ionode, ionode_id
9USE mp_world, ONLY : mpime, nproc
10USE mp_pools, ONLY: intra_pool_comm
11USE mp_wave, ONLY : mergewf,splitwf
12USE io_global, ONLY : stdout,ionode
13USE lsda_mod, ONLY :nspin
14USE gvect, ONLY : gstart
15USE mp, ONLY : mp_sum
16USE mp_world,             ONLY : world_comm
17
18
19
20implicit none
21type(exc), intent(in) :: a_in
22type(exc):: a_out
23type(exc_r):: a_in_rt
24type(exc_r):: a_tmp_rt
25type(ii_mat) :: iimat
26type(vww_prod) :: vww
27type(fft_cus) :: fc
28
29
30COMPLEX(kind=DP), allocatable :: vwwg_t(:,:)
31COMPLEX(kind=DP), ALLOCATABLE :: evc_g(:)
32
33
34!real(kind=dp), allocatable :: phivwwr(:,:)
35COMPLEX(kind=DP) :: csca
36
37integer ii, iv,ispin, iimax
38
39logical debug
40
41call start_clock('direct_v_exc')
42debug=.false.
43
44!if(debug) then
45!   if(ionode) write(stdout,*) 'Direct_v_exc #1'
46!endif
47
48
49! allocate tmp matrix
50call initialize_exc_r(a_tmp_rt)
51a_tmp_rt%nrxxt=fc%nrxxt
52a_tmp_rt%numb_v=a_in%numb_v
53a_tmp_rt%label=12
54allocate(a_tmp_rt%ar(a_tmp_rt%nrxxt,a_tmp_rt%numb_v))
55
56
57! FFT a_in to real space (dual grid)
58call initialize_exc_r(a_in_rt)
59call fft_a_exc(a_in,fc,a_in_rt)
60
61
62! read iimat, that tells us for every w_iv which other valence band ivp is
63! overlapping, and read the corresponding v*w_iv*w_ivp(G) products
64
65call initialize_imat(iimat)
66
67do ispin=1,nspin
68! note that for spin-polarized case, probably this do-loop will have to include
69! also the rest of the subroutine, or something like that
70   call read_iimat(iimat,ispin)
71enddo
72
73
74if(debug) then
75   if(ionode) write(stdout,*) 'Direct_v_exc #5'
76   if(ionode) write(stdout,*) 'a_in%numb_v=',a_in%numb_v
77   if(ionode) write(stdout,*) 'a_in%npw=',a_in%npw
78   if(ionode) write(stdout,*) 'iimat%np_max=',iimat%np_max
79endif
80
81call initialize_vww_prod(vww)
82call read_vww_prod(1,a_in%numb_v,a_in%npw,iimat%np_max,iimat,vww)
83
84!if(debug) then
85!   if(ionode) write(stdout,*) 'Direct_v_exc #6'
86!endif
87
88! for every element iv of the excitonic wavefunction vector, here we FFT all
89! the available v*w_iv*w_ivp(G) products, multiply by a_in_rt%ar(:,ivp)
90! sum over ivp, and FFT back
91
92allocate(vwwg_t(fc%npwt,iimat%np_max))
93!allocate(phivwwr(fc%nrxxt,a_in%numb_v))
94allocate(evc_g(fc%ngmt_g ))
95
96
97
98a_tmp_rt%ar(1:a_tmp_rt%nrxxt,1:a_tmp_rt%numb_v) =0.d0
99do iv=1,a_in%numb_v
100
101   vwwg_t(1:fc%npwt,1:iimat%np_max)=dcmplx(0.d0,0.d0)
102   iimax=0
103   do ii=1,iimat%np_max
104      if (iimat%iimat(ii,iv)>0) then
105         iimax=iimax+1
106      else
107         exit
108      endif
109   enddo
110   if(iimax>0) then
111      call reorderwfp_col(iimax,vww%npw,fc%npwt,vww%vww(1,1,iv),vwwg_t, vww%npw,fc%npwt, &
112           & ig_l2g,fc%ig_l2gt,fc%ngmt_g,mpime, nproc,intra_pool_comm )
113   endif
114!!!!!!!!!!!!!!!!
115!   do ii=1,iimat%np_max
116!      if (iimat%iimat(ii,iv)==0) exit
117!      if(fc%dual_t==4.d0) then
118!         vwwg_t(1:fc%npwt,ii)= vww%vww(1:fc%npwt,ii,iv)
119!      else
120!        call reorderwfp_col(1,vww%npw,fc%npwt,vww%vww(1,ii,iv),vwwg_t(1,ii), vww%npw,fc%npwt, &
121!           & ig_l2g,fc%ig_l2gt,fc%ngmt_g,mpime, nproc,intra_pool_comm )
122!
123!      endif
124!
125!
126!   enddo
127
128
129   do ii=1,iimat%np_max,2
130      if(debug) then
131         if(ionode) write(stdout,*) 'ii,iv,iimat', ii,iv,iimat%iimat(ii,iv)
132      endif
133      if (iimat%iimat(ii,iv)==0) exit
134      psic(:)=(0.d0,0.d0)
135      if ((ii==iimat%np_max).or.(iimat%iimat(ii+1,iv)==0)) then
136          psic(fc%nlt(1:fc%npwt))  = vwwg_t(1:fc%npwt,ii)
137          psic(fc%nltm(1:fc%npwt)) = CONJG( vwwg_t(1:fc%npwt,ii) )
138      else
139          psic(fc%nlt(1:fc%npwt))=vwwg_t(1:fc%npwt,ii)+(0.d0,1.d0)*vwwg_t(1:fc%npwt,ii+1)
140          psic(fc%nltm(1:fc%npwt))=CONJG(vwwg_t(1:fc%npwt,ii))+(0.d0,1.d0)*CONJG(vwwg_t(1:fc%npwt,ii+1))
141      endif
142      CALL cft3t( fc, psic, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, 2 )
143      a_tmp_rt%ar(1:fc%nrxxt,iv)= &
144            &a_tmp_rt%ar(1:fc%nrxxt,iv)+DBLE(psic(1:fc%nrxxt))*a_in_rt%ar(1:a_in_rt%nrxxt,iimat%iimat(ii,iv))
145      if ((ii/=iimat%np_max).and.(iimat%iimat(ii+1,iv)/=0)) then
146         a_tmp_rt%ar(1:fc%nrxxt,iv)=&
147         &a_tmp_rt%ar(1:fc%nrxxt,iv)+DIMAG(psic(1:fc%nrxxt))*a_in_rt%ar(1:a_in_rt%nrxxt,iimat%iimat(ii+1,iv))
148      endif
149   enddo
150enddo
151
152!if(debug) then
153!   if(ionode) write(stdout,*) 'Direct_v_exc #7'
154!endif
155
156
157call free_memory_exc_a_r(a_in_rt)
158
159call fftback_a_exc(a_tmp_rt,fc,a_out)
160
161! free memory
162call free_memory_exc_a_r(a_tmp_rt)
163call free_imat(iimat)
164call free_vww_prod(vww)
165deallocate(vwwg_t)
166!deallocate(vwwr_t)
167deallocate(evc_g)
168
169
170call stop_clock('direct_v_exc')
171end subroutine
172
173
174