1!
2! Copyright (C) 2001-2013 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!
9!----------------------------------------------------------------------------
10SUBROUTINE wfc_gamma_real(itask,ispin)
11
12!this subroutine writes the wfcs on real space - coarse grid
13!on disk
14!it works only at GAMMA
15!it supports US pseudopotentials
16!it also sets uo the array bec_gw
17
18
19  USE kinds,                ONLY : DP
20  USE gvect,                ONLY : gstart
21  USE gvecs,              ONLY : doublegrid
22   USE io_files,             ONLY : iunwfc, nwordwfc, diropn
23  USE wvfct,                ONLY : nbnd, npwx, npw, wg, et
24  USE mp,                   ONLY : mp_bcast
25  USE io_global,            ONLY : stdout
26  USE klist,                ONLY : lgauss, degauss, ngauss, nks, &
27                                   nkstot, wk, xk, nelec, nelup, neldw, &
28                                   two_fermi_energies, igk_k
29  USE lsda_mod,             ONLY : lsda, nspin, current_spin, isk
30  USE wavefunctions, ONLY : evc, psic
31  USE io_files,             ONLY : diropn
32  USE wannier_gw,           ONLY : becp_gw, becp_gw_c, l_verbose
33  USE uspp
34  USE control_flags,    ONLY : gamma_only
35  USE fft_base,             ONLY : dfftp, dffts
36  USE fft_interfaces,       ONLY : fwfft, invfft
37
38  IMPLICIT NONE
39
40  INTEGER, EXTERNAL :: find_free_unit
41
42  INTEGER, INTENT(in) :: itask !if ==1 consider subspace{c'}
43  INTEGER, INTENT(in) :: ispin!spin variable 1,2
44
45  !
46  INTEGER :: ikb, jkb, ijkb0, ih, jh, ijh, na, np
47  INTEGER :: ir, is, ig, ibnd, ik
48  INTEGER :: iunwfcreal
49  LOGICAL :: exst
50  REAL(kind=DP), ALLOCATABLE :: tmpreal(:)
51
52  if(l_verbose) write(stdout,*) 'FUNCTION WFC_REAL' !ATTENZIONE
53    FLUSH(stdout)
54
55    allocate(tmpreal(dffts%nnr))
56
57  IF(.not.gamma_only) THEN
58       write(stdout,*) ' wfc_gamma_real only for GAMMA'
59       stop
60  ENDIF
61
62  iunwfcreal=find_free_unit()
63  CALL diropn( iunwfcreal, 'real_whole', dffts%nnr, exst )
64
65       !
66
67  !
68
69  if  ( nkb > 0 .and. okvan) then
70     CALL init_us_2( npw, igk_k(1,1), xk(1,1), vkb )
71     if(itask/=1) then
72        !CALL ccalbec( nkb, npwx, npw, nbnd, becp_gw, vkb, evc )
73     else
74        !CALL ccalbec( nkb, npwx, npw, nbnd, becp_gw_c, vkb, evc )
75     endif
76  endif
77     !
78
79          !
80          ! ... here we compute the band energy: the sum of the eigenvalues
81          !
82
83  if(gstart==2) then
84     do ibnd=1,nbnd
85        evc(1,ibnd)=dble(evc(1,ibnd))
86     enddo
87  endif
88
89  DO ibnd = 1, nbnd, 2
90      if(l_verbose) write(stdout,*) 'IBND:',ibnd
91     FLUSH(stdout)
92     !
93     psic(:) = ( 0.D0, 0.D0 )
94        !
95     IF ( ibnd < nbnd ) THEN
96                !
97                ! ... two ffts at the same time
98                !
99        psic(dffts%nl(igk_k(1:npw,1)))  = evc(1:npw,ibnd) + &
100             ( 0.D0, 1.D0 ) * evc(1:npw,ibnd+1)
101        psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( evc(1:npw,ibnd) - &
102             ( 0.D0, 1.D0 ) * evc(1:npw,ibnd+1) )
103                !
104     ELSE
105        !
106        psic(dffts%nl(igk_k(1:npw,1)))  = evc(1:npw,ibnd)
107        psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( evc(1:npw,ibnd) )
108           !
109     END IF
110             !
111      if(l_verbose) write(stdout,*) 'before'
112     FLUSH(stdout)
113
114     CALL invfft ('Wave', psic, dffts)
115                    !
116        !
117        ! ... increment the charge density ...
118        !
119             !
120
121     if(l_verbose) write(stdout,*) 'after'
122     FLUSH(stdout)
123
124     tmpreal(:)= DBLE(psic(:))
125     CALL davcio( tmpreal,dffts%nnr,iunwfcreal,ibnd+(ispin-1)*nbnd,1)
126     if(ibnd+1 <= nbnd) then
127        tmpreal(:)=dimag(psic(:))
128        CALL davcio( tmpreal,dffts%nnr,iunwfcreal,ibnd+1+(ispin-1)*nbnd,1)
129     endif
130
131             !
132  END DO
133          !
134
135  close(iunwfcreal)
136  deallocate(tmpreal)
137       !
138  END SUBROUTINE
139
140
141  SUBROUTINE write_wfc_plot(itask)
142!save wannier functions on disk for plotting
143     USE io_files,             ONLY : nwordwfc, diropn
144     USE wavefunctions, ONLY : evc
145
146    implicit none
147
148    INTEGER, EXTERNAL :: find_free_unit
149
150    INTEGER, INTENT(in) :: itask!0 save MLWF, 1 save ULWF
151
152    INTEGER :: iunplot
153    LOGICAL :: exst
154
155     iunplot=find_free_unit()
156     if(itask==0) then
157        CALL diropn( iunplot, 'wfc_mlwf', nwordwfc, exst )
158     else
159        CALL diropn( iunplot, 'wfc_ulwf', nwordwfc, exst )
160     endif
161     CALL davcio(evc,2*nwordwfc,iunplot,1,1)
162     close(iunplot)
163    return
164  END SUBROUTINE write_wfc_plot
165