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
10!this routine calculate the terms  \psi_i(r)\_psi_v(sc)(r)
11!and write them on disk on global G grid
12
13  subroutine semicore(n_semicore, num_nbnds,ispin)
14!NOT_TO_BE_INCLUDED_START
15     USE io_global,             ONLY : stdout, ionode,ionode_id
16     USE io_files,              ONLY : diropn,prefix,tmp_dir
17     use pwcom
18     USE wavefunctions,  ONLY : evc
19     USE kinds,                 ONLY : DP
20     USE gvect,        ONLY : ngm_g, ig_l2g, gstart
21     USE mp,           ONLY : mp_sum, mp_barrier, mp_bcast
22     USE mp_wave,      ONLY : mergewf,splitwf
23     USE mp_pools,     ONLY : intra_pool_comm, inter_pool_comm, intra_pool_comm
24     USE mp_world,     ONLY : world_comm, mpime, nproc
25     USE fft_base,             ONLY : dfftp, dffts
26     USE fft_interfaces,       ONLY : fwfft, invfft
27     USE wavefunctions, ONLY : psic
28     USE wvfct,    ONLY : et
29
30    implicit none
31
32    INTEGER, EXTERNAL :: find_free_unit
33    INTEGER, INTENT(in) :: n_semicore!number of semicore states
34    INTEGER, INTENT(in) :: num_nbnds!total KS states considered
35    INTEGER, INTENT(in) :: ispin!spin channel
36
37    REAL(kind=DP), ALLOCATABLE :: psi_sc(:,:)
38    COMPLEX(kind=DP), ALLOCATABLE :: prod(:), prod_g(:,:), prod_g_tot(:)
39    INTEGER :: iv, iun, ii
40    INTEGER :: npwx_g
41
42
43!fft trasform semicore states to R space
44    allocate(psi_sc(dfftp%nnr,n_semicore))
45    allocate(prod(dfftp%nnr), prod_g(npw,2),prod_g_tot(ngm_g))
46
47    do iv=1,n_semicore,2
48       psic(:)=(0.d0,0.d0)
49       if(iv<n_semicore) then
50          psic(dffts%nl(igk_k(1:npw,1)))  = evc(1:npw,iv) + &
51               ( 0.D0, 1.D0 ) * evc(1:npw,iv+1)
52          psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( evc(1:npw,iv) - &
53                    ( 0.D0, 1.D0 ) * evc(1:npw,iv+1) )
54
55       else
56          psic(dffts%nl(igk_k(1:npw,1)))  = evc(1:npw,iv)
57          psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( evc(1:npw,iv) )
58       endif
59       CALL invfft ('Wave', psic, dffts)
60       psi_sc(1:dfftp%nnr,iv)=dble(psic(1:dfftp%nnr))
61       if(iv< n_semicore) psi_sc(1:dfftp%nnr,iv+1)=dimag(psic(1:dfftp%nnr))
62
63    enddo
64
65
66!write header of file with KS energies of all states (just in case..)
67
68    npwx_g=npwx
69    call mp_sum(npwx_g, world_comm)
70
71    if(ionode) then
72       iun =  find_free_unit()
73       if(ispin==1) then
74          open( unit= iun, file=trim(tmp_dir)//trim(prefix)//'.sc_states', status='unknown',form='unformatted')
75       else
76          open( unit= iun, file=trim(tmp_dir)//trim(prefix)//'.sc_states2', status='unknown',form='unformatted')
77       endif
78       write(iun) num_nbnds
79       write(iun) n_semicore
80       write(iun) npwx_g
81       write(iun) et(1:num_nbnds,ispin)
82    endif
83
84
85!loop (double) on KS states
86    do ii=n_semicore+1,num_nbnds
87
88!fft
89
90       psic(:)=(0.d0,0.d0)
91       psic(dffts%nl(igk_k(1:npw,1)))  = evc(1:npw,ii)
92       psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( evc(1:npw,ii) )
93
94       CALL invfft ('Wave', psic, dffts)
95
96!calculate product
97       do iv=1,n_semicore,2
98          if(iv<n_semicore) then
99             prod(1:dfftp%nnr)=dcmplx(dble(psic(1:dfftp%nnr))*psi_sc(1:dfftp%nnr,iv),&
100                  &dble(psic(1:dfftp%nnr))*psi_sc(1:dfftp%nnr,iv+1))
101          else
102             prod(1:dfftp%nnr)=dble(psic(1:dfftp%nnr))
103             prod(1:dfftp%nnr)=prod(1:dfftp%nnr)*psi_sc(1:dfftp%nnr,iv)
104          endif
105!fft back
106
107          CALL fwfft ('Rho', prod, dfftp)
108          if(iv==n_semicore) then
109             prod_g(1:npw,1)=prod(dfftp%nl(1:npw))
110             if(gstart==2) then
111                write(stdout,*) 'Putting to zero:', iv,ii, prod_g(1,1)
112                prod_g(1,1)=(0.d0,0.d0)
113             endif
114          else
115             prod_g(1:npw, 1)= 0.5d0*(prod(dfftp%nl(1:npw))+conjg( prod(dfftp%nlm(1:npw))))
116             prod_g(1:npw, 2)= (0.d0,-0.5d0)*(prod(dfftp%nl(1:npw)) - conjg(prod(dfftp%nlm(1:npw))))
117             if(gstart==2) then
118                write(stdout,*)'Putting to zero:', iv,ii, prod_g(1,1)
119                write(stdout,*)'Putting to zero:', iv+1,ii, prod_g(1,2)
120                prod_g(1,1:2)=(0.d0,0.d0)
121             endif
122          endif
123
124!merge
125
126          prod_g_tot(:)=(0.d0,0.d0)
127          call mergewf(prod_g(:,1),prod_g_tot,npw,ig_l2g,mpime,nproc,ionode_id,intra_pool_comm)
128!put it on disk
129          if(ionode) then
130             write(iun) prod_g_tot(1:npwx_g)
131          endif
132          if(iv<n_semicore) then
133
134!merge
135             prod_g_tot(:)=(0.d0,0.d0)
136             call mergewf(prod_g(:,2),prod_g_tot,npw,ig_l2g,mpime,nproc,ionode_id,intra_pool_comm)
137!put it on disk
138             if(ionode) then
139             write(iun) prod_g_tot(1:npwx_g)
140          endif
141
142
143          endif
144       enddo
145
146
147
148
149    enddo
150
151
152!now put on disk KS states on global G
153    do ii=n_semicore+1,num_nbnds
154       prod_g_tot(:)=(0.d0,0.d0)
155       call mergewf(evc(:,ii),prod_g_tot,npw,ig_l2g,mpime,nproc,ionode_id,intra_pool_comm)
156!put it on disk
157       if(ionode) then
158          write(iun) prod_g_tot(1:npwx_g)
159       endif
160    enddo
161
162
163    deallocate(psi_sc,prod_g_tot)
164    close(iun)
165    return
166!NOT_TO_BE_INCLUDED_END
167  end subroutine semicore
168