1MODULE contract_w
2
3USE kinds, ONLY : DP,sgl
4use bse_basic_structures, ONLY : ii_mat
5
6SAVE
7
8type(ii_mat) :: iimat_contract
9REAL(kind=sgl), POINTER  :: vphipizeta_save(:,:)
10REAL(kind=sgl), POINTER  :: vww_save(:,:)
11COMPLEX(kind=sgl), POINTER  :: vphipizeta_save_g(:,:)
12COMPLEX(kind=sgl), POINTER  :: vww_save_g(:,:)
13INTEGER, POINTER ::  vpmax_ii(:),vpmax_ii_start(:),vpmax_ii_end(:)
14INTEGER :: vpmax_tot
15
16CONTAINS
17
18subroutine free_memory_contrac_w
19  use bse_wannier, ONLY : l_gtrick
20  implicit none
21  if(.not.l_gtrick) then
22     deallocate(vphipizeta_save)
23  else
24     deallocate(vphipizeta_save_g)
25  endif
26  deallocate(vpmax_ii)
27  deallocate(vpmax_ii_start,vpmax_ii_end)
28  if(.not.l_gtrick) then
29     deallocate(vww_save)
30  else
31     deallocate(vww_save_g)
32  endif
33end subroutine free_memory_contrac_w
34
35subroutine contract_w_build(fc)
36! this subroutine computes the w part of the direct term of the exc Hamiltonian
37
38USE wvfct,                 ONLY :  npw
39USE fft_custom_gwl
40use bse_basic_structures
41use exciton
42USE gvect
43use bse_wannier, ONLY: l_truncated_coulomb, &
44           truncation_radius, l_gtrick
45USE constants,        ONLY : e2, fpi
46USE cell_base,        ONLY : tpiba,omega,tpiba2
47!USE io_files,             ONLY : find_free_unit, prefix, diropn
48USE io_files,             ONLY :  prefix, diropn
49USE wavefunctions, ONLY :  psic
50USE io_global, ONLY : stdout, ionode, ionode_id
51USE mp_world, ONLY : mpime, nproc
52USE mp_pools, ONLY: intra_pool_comm
53USE mp_wave, ONLY : mergewf,splitwf
54USE polarization
55USE lsda_mod, ONLY :nspin
56USE io_global, ONLY : stdout,ionode
57USE mp,          ONLY :mp_barrier
58USE mp_world,             ONLY : world_comm
59
60
61
62
63
64implicit none
65INTEGER, EXTERNAL :: find_free_unit
66logical :: debug=.false.
67
68type(bse_z) :: z
69type(polaw) :: pw
70
71type(fft_cus) :: fc
72type(ii_mat) :: iimat
73
74
75
76REAL(kind=DP), ALLOCATABLE :: fac(:)
77COMPLEX(kind=DP), ALLOCATABLE :: p_basis(:,:)
78COMPLEX(kind=DP), ALLOCATABLE :: p_basis_t(:,:)
79REAL(kind=DP), ALLOCATABLE :: p_basis_r(:,:)
80REAL(kind=DP), ALLOCATABLE :: zvphi(:)
81REAL(kind=DP), ALLOCATABLE :: zvv(:)
82COMPLEX(kind=DP), ALLOCATABLE :: evc_g(:)
83
84INTEGER ::iungprod
85INTEGER :: ig,ii,iv,ispin
86REAL(kind=DP) :: qq
87LOGICAL       :: exst
88
89
90INTEGER :: vpmax,k
91REAL(kind=DP), allocatable :: zp(:,:)
92REAL(kind=DP), allocatable :: pizeta(:,:)
93REAL(kind=DP), allocatable :: vphipizeta(:,:)
94
95INTEGER ::  kilobytes
96
97
98call start_clock('direct_w_exc')
99!CALL memstat( kilobytes )
100!write(stdout,*) 'memory0', kilobytes
101!FLUSH(stdout)
102
103
104
105
106! read  iimat
107call initialize_imat(iimat)
108
109do ispin=1,nspin
110   call read_iimat(iimat,ispin)
111enddo
112
113! read z terms
114call initialize_bse_z(z)
115call read_z(1,iimat,z)
116
117FLUSH( stdout )
118
119! get Coulomb potential
120allocate(fac(npw))
121if(l_truncated_coulomb) then
122   do ig=1,npw
123      qq = g(1,ig)**2.d0 + g(2,ig)**2.d0 + g(3,ig)**2.d0
124      if (qq > 1.d-8) then
125         fac(ig)=(e2*fpi/(tpiba2*qq))*(1.d0-dcos(dsqrt(qq)*truncation_radius*tpiba))
126      else
127         fac(ig)=e2*fpi*(truncation_radius**2.d0/2.d0)
128     endif
129   enddo
130    fac(:)=fac(:)/omega
131else
132
133   fac(:)=0.d0
134   fac(1:npw)=vg_q(1:npw)
135endif
136
137
138! read polarization basis and multiply per V
139
140iungprod = find_free_unit()
141allocate(p_basis(npw,z%numw_prod))
142CALL diropn( iungprod, 'wiwjwfc_red', npw*2, exst )
143
144do ii=1,z%numw_prod
145   call davcio(p_basis(:,ii),npw*2,iungprod,ii,-1)
146   p_basis(1:npw,ii)=p_basis(1:npw,ii)*dcmplx(fac(1:npw))
147enddo
148
149call mp_barrier(world_comm)
150
151close(iungprod)
152!CALL memstat( kilobytes )
153!write(stdout,*) 'memory1', kilobytes
154!FLUSH(stdout)
155
156! FFT to real space (dual grid)
157allocate(p_basis_t(fc%npwt,z%numw_prod))
158allocate(p_basis_r(fc%nrxxt,z%numw_prod))
159allocate(evc_g(fc%ngmt_g ))
160
161
162if(fc%dual_t==4.d0) then
163   p_basis_t(1:fc%npwt,1:z%numw_prod)=p_basis(1:npw,1:z%numw_prod)
164else
165    call reorderwfp_col(z%numw_prod,npw,fc%npwt,p_basis(1,1),p_basis_t(1,1),npw,fc%npwt, &
166           & ig_l2g,fc%ig_l2gt,fc%ngmt_g,mpime, nproc,intra_pool_comm )
167
168!   do ii=1,z%numw_prod
169!      call mergewf(p_basis(:,ii),evc_g,a_in%npw,ig_l2g,mpime,nproc,ionode_id,intra_pool_comm)
170!      call splitwf(p_basis_t(:,ii),evc_g,fc%npwt,fc%ig_l2gt,mpime,nproc,ionode_id,intra_pool_comm)
171!   enddo
172endif
173
174deallocate(evc_g)
175deallocate(p_basis)
176
177call start_clock('direct_w_cft3t')
178do ii=1,z%numw_prod,2
179   psic(1:fc%nrxxt)=(0.d0,0.d0)
180   if (ii==z%numw_prod) then
181      psic(fc%nlt(1:fc%npwt))  = p_basis_t(1:fc%npwt,ii)
182      psic(fc%nltm(1:fc%npwt)) = CONJG( p_basis_t(1:fc%npwt,ii) )
183   else
184      psic(fc%nlt(1:fc%npwt))=p_basis_t(1:fc%npwt,ii)+(0.d0,1.d0)*p_basis_t(1:fc%npwt,ii+1)
185      psic(fc%nltm(1:fc%npwt))=CONJG(p_basis_t(1:fc%npwt,ii))+(0.d0,1.d0)*CONJG(p_basis_t(1:fc%npwt,ii+1))
186   endif
187      CALL cft3t( fc, psic, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, 2 )
188       p_basis_r(1:fc%nrxxt,ii)= DBLE(psic(1:fc%nrxxt))
189      if(ii/=z%numw_prod) p_basis_r(1:fc%nrxxt,ii+1)= DIMAG(psic(1:fc%nrxxt))
190enddo
191call stop_clock('direct_w_cft3t')
192
193deallocate(p_basis_t)
194
195!read P
196call initialize_polaw(pw)
197call read_polaw_global(0, pw)
198
199
200call mp_barrier(world_comm)
201
202!CALL memstat( kilobytes )
203!write(stdout,*) 'memory2', kilobytes
204!FLUSH(stdout)
205
206
207! allocate tmp matrix
208
209
210!compute line by line the output excitonic vector
211
212!!!!!!!!!!!!!!!!!dgemm subroutine!!!!!!!!!!!!!!!!!!!!!
213call start_clock('direct_w_dgemv')
214!write(stdout,*) 'memory2',z%numw_prod,iimat%np_max
215allocate(zp(z%numw_prod,iimat%np_max))
216!!allocate(pizeta(z%numw_prod,iimat%np_max))
217!!allocate(vphipizeta(fc%nrxxt,iimat%np_max))
218!calculate vpmax_tot
219vpmax_tot=0
220allocate(vpmax_ii( iimat%numb_v))
221allocate(vpmax_ii_start( iimat%numb_v))
222allocate(vpmax_ii_end( iimat%numb_v))
223do iv=1, iimat%numb_v
224
225   vpmax=0
226   vpmax_ii_start(iv)=vpmax_tot+1
227   do ii=1, iimat%np_max
228      if (iimat%iimat(ii,iv)==0) cycle
229      vpmax=vpmax+1
230   enddo
231   vpmax_ii(iv)=vpmax
232   vpmax_tot=vpmax_tot+vpmax
233   vpmax_ii_end(iv)=vpmax_tot
234enddo
235write(stdout,*)  'VPHIZETA_SAVE :', fc%nrxxt,vpmax_tot
236FLUSH(stdout)
237if(.not.l_gtrick) then
238   allocate(vphipizeta_save(fc%nrxxt,vpmax_tot))
239else
240   allocate(vphipizeta_save_g(fc%npwt,vpmax_tot))
241endif
242write(stdout,*)  'VPHIZETA_SAVE :', fc%nrxxt,vpmax_tot
243!
244!CALL memstat( kilobytes )
245!write(stdout,*) 'memory3', kilobytes
246!FLUSH(stdout)
247do iv=1, iimat%numb_v
248   zp(1:z%numw_prod,1:iimat%np_max)=0.d0
249   vpmax=0
250
251   call start_clock('dgemv1')
252   do ii=1, iimat%np_max
253      if (iimat%iimat(ii,iv)==0) cycle
254      vpmax=vpmax+1
255      do k=1, z%numw_prod
256         zp(k,vpmax)=z%z(k,ii,iv)!ATTENZIONE era zp(ii)
257      enddo
258   enddo
259
260  if(vpmax>0) then
261     allocate(pizeta(z%numw_prod,vpmax))
262     allocate(vphipizeta(fc%nrxxt,vpmax))
263  else
264     allocate(pizeta(z%numw_prod,1))
265     allocate(vphipizeta(fc%nrxxt,1))
266  endif
267   call stop_clock('dgemv1')
268
269
270   call start_clock('dgemv2')
271
272
273   call dgemm('N','N', z%numw_prod,vpmax, z%numw_prod,1.d0,pw%pw,z%numw_prod,zp,&
274              z%numw_prod,0.d0,pizeta,z%numw_prod)
275   call stop_clock('dgemv2')
276
277
278   call start_clock('dgemv3')
279   call dgemm('N','N', fc%nrxxt, vpmax, z%numw_prod,1.d0,p_basis_r(1,1),fc%nrxxt,&
280              pizeta(1,1), z%numw_prod, 0.d0, vphipizeta(1,1),fc%nrxxt)
281   call stop_clock('dgemv3')
282   if(.not.l_gtrick) then
283      vphipizeta_save(1:fc%nrxxt,vpmax_ii_start(iv):vpmax_ii_end(iv))= real(vphipizeta(1:fc%nrxxt,1:vpmax))
284   else
285      do ii=1,vpmax,2
286         if(ii==vpmax) then
287            psic(1:fc%nrxxt)=dcmplx(vphipizeta(1:fc%nrxxt,ii),0.d0)
288         else
289            psic(1:fc%nrxxt)=dcmplx(vphipizeta(1:fc%nrxxt,ii),vphipizeta(1:fc%nrxxt,ii+1))
290         endif
291         CALL cft3t(fc, psic, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, -2 )
292         if(ii==vpmax) then
293            vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1)=cmplx(psic(fc%nlt(1:fc%npwt)))
294         else
295            vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1)=&
296                 &cmplx(0.5d0*(psic(fc%nlt(1:fc%npwt))+conjg( psic(fc%nltm(1:fc%npwt)))))
297            vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1+1)=&
298                 &cmplx((0.d0,-0.5d0)*(psic(fc%nlt(1:fc%npwt)) - conjg(psic(fc%nltm(1:fc%npwt)))))
299         endif
300      enddo
301   endif
302
303   !
304!
305   deallocate(pizeta)
306   deallocate(vphipizeta)
307call mp_barrier(world_comm)
308enddo
309call stop_clock('direct_w_dgemv')
310
311deallocate(zp)
312deallocate(p_basis_r)
313
314call free_bse_z(z)
315call free_memory_polaw(pw)
316call free_imat(iimat)
317
318
319!CALL memstat( kilobytes )
320!write(stdout,*) 'memory4', kilobytes
321!FLUSH(stdout)
322call stop_clock('direct_w_exc')
323
324return
325end subroutine contract_w_build
326
327
328
329subroutine contract_w_apply(a_in,fc,a_out)
330! this subroutine computes the w part of the direct term of the exc Hamiltonian
331
332USE fft_custom_gwl
333use bse_basic_structures
334use exciton
335USE gvect
336use bse_wannier, ONLY: l_truncated_coulomb, &
337           truncation_radius,l_gtrick
338USE constants,        ONLY : e2, fpi
339USE cell_base,        ONLY : tpiba,omega,tpiba2
340!USE io_files,             ONLY : find_free_unit, prefix, diropn
341USE io_files,             ONLY :  prefix, diropn
342USE wavefunctions, ONLY :  psic
343USE io_global, ONLY : stdout, ionode, ionode_id
344USE mp_world, ONLY : mpime, nproc
345USE mp_pools, ONLY: intra_pool_comm
346USE mp_wave, ONLY : mergewf,splitwf
347USE polarization
348USE lsda_mod, ONLY :nspin
349USE io_global, ONLY : stdout,ionode
350USE mp,          ONLY :mp_barrier
351USE mp_world,             ONLY : world_comm
352
353
354
355
356
357implicit none
358INTEGER, EXTERNAL :: find_free_unit
359
360type(polaw) :: pw
361type(exc):: a_in
362type(exc):: a_out
363type(exc_r):: a_in_rt
364type(exc_r):: a_tmp_rt
365type(fft_cus) :: fc
366
367
368
369
370
371INTEGER ::iungprod
372INTEGER :: ig,ii,iv,ispin
373REAL(kind=DP) :: qq
374LOGICAL       :: exst
375
376
377INTEGER ::k
378
379
380logical debug
381
382call start_clock('direct_w_contract')
383
384debug=.false.
385! read  iimat
386
387!FFT the input excitonic vector to real space (dual grid)
388call initialize_exc_r(a_in_rt)
389call fft_a_exc(a_in,fc,a_in_rt)
390
391call mp_barrier(world_comm)
392
393! allocate tmp matrix
394call initialize_exc_r(a_tmp_rt)
395a_tmp_rt%nrxxt=fc%nrxxt
396a_tmp_rt%numb_v=a_in%numb_v
397a_tmp_rt%label=12
398allocate(a_tmp_rt%ar(a_tmp_rt%nrxxt,a_tmp_rt%numb_v))
399
400!compute line by line the output excitonic vector
401
402!!!!!!!!!!!!!!!!!dgemm subroutine!!!!!!!!!!!!!!!!!!!!!
403call start_clock('contract_w_dgemv')
404
405
406!
407a_tmp_rt%ar(1:a_tmp_rt%nrxxt,1:a_tmp_rt%numb_v) =0.d0
408!
409do iv=1, a_in%numb_v
410
411   if(.not.l_gtrick) then
412      do ii=1,vpmax_ii(iv)
413         a_tmp_rt%ar(1:fc%nrxxt,iv)= &
414              &a_tmp_rt%ar(1:fc%nrxxt,iv)+a_in_rt%ar(1:fc%nrxxt,iimat_contract%iimat(ii,iv))*&
415              &dble(vphipizeta_save(1:fc%nrxxt,vpmax_ii_start(iv)-1+ii))
416      enddo
417   else
418      do ii=1,vpmax_ii(iv),2
419         psic(1:fc%nrxxt)=(0.d0,0.d0)
420         if (ii==vpmax_ii(iv)) then
421            psic(fc%nlt(1:fc%npwt))  = dcmplx(vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii))
422            psic(fc%nltm(1:fc%npwt)) = dcmplx(CONJG( vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii) ))
423         else
424            psic(fc%nlt(1:fc%npwt))= dcmplx(vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii))+&
425                 &(0.0,1.0)* dcmplx(vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii+1))
426            psic(fc%nltm(1:fc%npwt))=DCMPLX(CONJG( vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii))+&
427                 &(0.0,1.0)*CONJG( vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii+1)))
428         endif
429         CALL cft3t( fc, psic, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, 2 )
430         a_tmp_rt%ar(1:fc%nrxxt,iv)= &
431              &a_tmp_rt%ar(1:fc%nrxxt,iv)+a_in_rt%ar(1:fc%nrxxt,iimat_contract%iimat(ii,iv))*&
432              &dble(psic(1:fc%nrxxt))
433         if (ii/=vpmax_ii(iv)) then
434             a_tmp_rt%ar(1:fc%nrxxt,iv)= &
435                  &a_tmp_rt%ar(1:fc%nrxxt,iv)+a_in_rt%ar(1:fc%nrxxt,iimat_contract%iimat(ii+1,iv))*&
436                  &dimag(psic(1:fc%nrxxt))
437
438         endif
439      enddo
440   endif
441
442
443
444enddo
445call stop_clock('contract_w_dgemv')
446
447
448
449
450call free_memory_exc_a_r(a_in_rt)
451
452
453
454call start_clock('wdirect_fftback')
455!FFT back to provide the output excitonic wave vector in G-space
456call fftback_a_exc(a_tmp_rt,fc,a_out)
457call stop_clock('wdirect_fftback')
458
459call free_memory_exc_a_r(a_tmp_rt)
460
461FLUSH( stdout )
462call stop_clock('direct_w_contract')
463
464return
465end subroutine contract_w_apply
466
467
468subroutine contract_v_build(fc)
469!TO BE CALLED AFTER CONTRACT_W_BUILD
470!IMPORTANT: REQUIRES IIMAT_CONTRACT
471! computes the v part of the direct term of the exc Hamiltonian
472USE fft_custom_gwl
473use bse_basic_structures
474use exciton
475USE wavefunctions, ONLY :  psic
476USE gvect,                 ONLY : ig_l2g
477USE io_global, ONLY : stdout, ionode, ionode_id
478USE mp_world, ONLY : mpime, nproc
479USE mp_pools, ONLY: intra_pool_comm
480USE mp_wave, ONLY : mergewf,splitwf
481USE io_global, ONLY : stdout,ionode
482USE lsda_mod, ONLY :nspin
483USE gvect, ONLY : gstart
484USE mp, ONLY : mp_sum
485USE mp_world,             ONLY : world_comm
486USE wvfct,                 ONLY :  npw
487USE bse_wannier, ONLY : l_gtrick
488
489implicit none
490
491
492type(vww_prod) :: vww
493type(fft_cus) :: fc
494
495
496COMPLEX(kind=DP), allocatable :: vwwg_t(:,:)
497COMPLEX(kind=DP), ALLOCATABLE :: evc_g(:)
498
499
500!real(kind=dp), allocatable :: phivwwr(:,:)
501COMPLEX(kind=DP) :: csca
502
503integer ii, iv,ispin, iimax
504
505logical debug
506
507write(stdout,*) 'Routine contract_v_build'
508
509debug=.false.
510
511!if(debug) then
512!   if(ionode) write(stdout,*) 'Direct_v_exc #1'
513!endif
514write(stdout,*) 'VWW_SAVE :', fc%nrxxt,vpmax_tot
515if(.not.l_gtrick) then
516   allocate(vww_save(fc%nrxxt,vpmax_tot))
517else
518   allocate(vww_save_g(fc%npwt,vpmax_tot))
519endif
520
521
522call initialize_vww_prod(vww)
523call read_vww_prod(1,iimat_contract%numb_v,npw,iimat_contract%np_max,iimat_contract,vww)
524
525!if(debug) then
526!   if(ionode) write(stdout,*) 'Direct_v_exc #6'
527!endif
528
529! for every element iv of the excitonic wavefunction vector, here we FFT all
530! the available v*w_iv*w_ivp(G) products, multiply by a_in_rt%ar(:,ivp)
531! sum over ivp, and FFT back
532
533allocate(vwwg_t(fc%npwt,iimat_contract%np_max))
534!allocate(phivwwr(fc%nrxxt,a_in%numb_v))
535allocate(evc_g(fc%ngmt_g ))
536
537
538write(stdout,*) 'ATT1'
539do iv=1,iimat_contract%numb_v
540
541   vwwg_t(1:fc%npwt,1:iimat_contract%np_max)=dcmplx(0.d0,0.d0)
542   iimax=0
543   do ii=1,iimat_contract%np_max
544      if (iimat_contract%iimat(ii,iv)>0) then
545         iimax=iimax+1
546      else
547         exit
548      endif
549   enddo
550   if(iimax>0) then
551      call reorderwfp_col(iimax,vww%npw,fc%npwt,vww%vww(1,1,iv),vwwg_t, vww%npw,fc%npwt, &
552           & ig_l2g,fc%ig_l2gt,fc%ngmt_g,mpime, nproc,intra_pool_comm )
553   endif
554
555   if(.not.l_gtrick) then
556      do ii=1,vpmax_ii(iv),2
557
558
559         psic(1:fc%nrxxt)=(0.d0,0.d0)
560         if (ii==vpmax_ii(iv)) then
561            psic(fc%nlt(1:fc%npwt))  = vwwg_t(1:fc%npwt,ii)
562            psic(fc%nltm(1:fc%npwt)) = CONJG( vwwg_t(1:fc%npwt,ii) )
563         else
564            psic(fc%nlt(1:fc%npwt))=vwwg_t(1:fc%npwt,ii)+(0.d0,1.d0)*vwwg_t(1:fc%npwt,ii+1)
565            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))
566         endif
567         CALL cft3t( fc, psic, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, 2 )
568         vww_save(1:fc%nrxxt,vpmax_ii_start(iv)+ii-1)=real(psic(1:fc%nrxxt))
569         if (ii/=vpmax_ii(iv)) then
570            vww_save(1:fc%nrxxt,vpmax_ii_start(iv)+ii)=aimag(psic(1:fc%nrxxt))
571         endif
572
573      enddo
574   else
575      do ii=1,vpmax_ii(iv)
576         vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1)=cmplx(vwwg_t(1:fc%npwt,ii))
577      end do
578   endif
579enddo
580
581!if(debug) then
582!   if(ionode) write(stdout,*) 'Direct_v_exc #7'
583!endif
584
585FLUSH( stdout )
586
587
588call free_vww_prod(vww)
589deallocate(vwwg_t)
590!deallocate(vwwr_t)
591deallocate(evc_g)
592
593
594end subroutine
595
596
597subroutine contract_v_apply(a_in,fc,a_out)
598! computes the v part of the direct term of the exc Hamiltonian
599USE fft_custom_gwl
600use bse_basic_structures
601use exciton
602USE wavefunctions, ONLY :  psic
603USE gvect,                 ONLY : ig_l2g
604USE io_global, ONLY : stdout, ionode, ionode_id
605USE mp_world, ONLY : mpime, nproc
606USE mp_pools, ONLY: intra_pool_comm
607USE mp_wave, ONLY : mergewf,splitwf
608USE io_global, ONLY : stdout,ionode
609USE lsda_mod, ONLY :nspin
610USE gvect, ONLY : gstart
611USE mp, ONLY : mp_sum
612USE mp_world,             ONLY : world_comm
613USE bse_wannier, ONLY : l_gtrick
614
615
616implicit none
617type(exc), intent(in) :: a_in
618type(exc):: a_out
619type(exc_r):: a_in_rt
620type(exc_r):: a_tmp_rt
621
622
623type(fft_cus) :: fc
624
625
626COMPLEX(kind=DP), allocatable :: vwwg_t(:,:)
627COMPLEX(kind=DP), ALLOCATABLE :: evc_g(:)
628
629COMPLEX(kind=DP) :: csca
630
631integer ii, iv,ispin, iimax
632
633logical debug
634
635call start_clock('direct_v_contract')
636debug=.false.
637
638call initialize_exc_r(a_tmp_rt)
639a_tmp_rt%nrxxt=fc%nrxxt
640a_tmp_rt%numb_v=a_in%numb_v
641a_tmp_rt%label=12
642allocate(a_tmp_rt%ar(a_tmp_rt%nrxxt,a_tmp_rt%numb_v))
643
644
645! FFT a_in to real space (dual grid)
646call initialize_exc_r(a_in_rt)
647call fft_a_exc(a_in,fc,a_in_rt)
648
649
650
651
652
653a_tmp_rt%ar(1:a_tmp_rt%nrxxt,1:a_tmp_rt%numb_v) =0.d0
654do iv=1,a_in%numb_v
655
656
657   if(.not.l_gtrick) then
658      do ii=1,vpmax_ii(iv)
659         a_tmp_rt%ar(1:fc%nrxxt,iv)= &
660              &a_tmp_rt%ar(1:fc%nrxxt,iv)+DBLE(vww_save(1:fc%nrxxt,vpmax_ii_start(iv)+ii-1))*&
661              &a_in_rt%ar(1:a_in_rt%nrxxt,iimat_contract%iimat(ii,iv))
662      enddo
663   else
664      do ii=1,vpmax_ii(iv),2
665         psic(1:fc%nrxxt)=(0.d0,0.d0)
666         if (ii==vpmax_ii(iv)) then
667            psic(fc%nlt(1:fc%npwt))  = dcmplx(vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1))
668            psic(fc%nltm(1:fc%npwt)) = dcmplx(CONJG(  vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1) ))
669         else
670            psic(fc%nlt(1:fc%npwt))=dcmplx(vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1)+&
671                 &(0.0,1.0)*vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1+1))
672            psic(fc%nltm(1:fc%npwt))=dcmplx(CONJG(vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1))&
673                 &+(0.0,1.0)*CONJG(vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1+1)))
674         endif
675         call start_clock('d_v_fft')
676         CALL cft3t( fc, psic, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, 2 )
677         call stop_clock('d_v_fft')
678         a_tmp_rt%ar(1:fc%nrxxt,iv)= &
679              &a_tmp_rt%ar(1:fc%nrxxt,iv)+DBLE(psic(1:fc%nrxxt))*&
680              &a_in_rt%ar(1:a_in_rt%nrxxt,iimat_contract%iimat(ii,iv))
681         if (ii/=vpmax_ii(iv)) then
682            a_tmp_rt%ar(1:fc%nrxxt,iv)= &
683                 &a_tmp_rt%ar(1:fc%nrxxt,iv)+DIMAG(psic(1:fc%nrxxt))*&
684                 &a_in_rt%ar(1:a_in_rt%nrxxt,iimat_contract%iimat(ii+1,iv))
685         endif
686      enddo
687   endif
688enddo
689
690
691
692call free_memory_exc_a_r(a_in_rt)
693
694call fftback_a_exc(a_tmp_rt,fc,a_out)
695
696! free memory
697call free_memory_exc_a_r(a_tmp_rt)
698
699
700
701call stop_clock('direct_v_contract')
702end subroutine
703
704
705
706
707
708
709
710END MODULE contract_w
711