1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
|
subroutine plot_excwfn(nstart,nend,vstate_r,fc)
! this subroutine computes and writes on file the nplot-th excitonic wavefunction
! to be read by pp.x
! note that this subroutine is working only for gamma only calculations
USE exciton
USE bse_wannier, ONLY:num_nbndv,&
r_hole,l_plotaverage
use bse_basic_structures
USE pwcom
USE cell_base, ONLY : alat, at, bg
USE fft_custom_gwl
USE io_global, ONLY : stdout,ionode,ionode_id
USE io_files, ONLY : tmp_dir,prefix
USE mp_world, ONLY : mpime, nproc
USE mp, ONLY: mp_sum
USE mp_world, ONLY : world_comm
!USE io_files, ONLY : find_free_unit
USE ions_base, ONLY : nat, tau, atm,ityp
implicit none
INTEGER, EXTERNAL :: find_free_unit
integer :: nplot,nstart,nend
type(v_state_r) :: vstate_r
type(exc_r) :: a_rt
type(fft_cus) :: fc
integer ::nxh,nyh,nzh,nh
INTEGER :: nr3s_start, nr3s_end
real(kind=dp), allocatable :: psi_exc(:)
real(kind=dp), allocatable :: psi_excio(:)
real(kind=dp), allocatable :: psi_excsum(:)
real(kind=dp), allocatable :: v_rh(:)
integer :: iv,ii,iplane,iz,ounit,ix,iy,iip
logical :: debug
CHARACTER(5) :: nfile
call start_clock('plot_excwfn')
debug=.true.
!check if all variables are ok from read_file in main
if (debug) then
if(ionode) then
!bg(:,i) are the reciprocal lattice vectors, b_i,
!in tpiba=2pi/alat units: b_i(:) = bg(:,i)/tpiba
!at(:,i) are the lattice vectors of the simulation cell, a_i,
!in alat units: a_i(:) = at(:,i)/alat
write(stdout,*) 'plotexcwfn bg(:,1)=',bg(1,1),bg(2,1),bg(3,1)
write(stdout,*) 'plotexcwfn bg(:,2)=',bg(1,2),bg(2,2),bg(3,2)
write(stdout,*) 'plotexcwfn bg(:,3)=',bg(1,3),bg(2,3),bg(3,3)
write(stdout,*) 'plotexcwfn alat=',alat
endif
endif
!find FFT grid point (dual grid) closer to r_hole (given in alat units)
nxh = nint ( (r_hole(1)*bg(1,1) + r_hole(2)*bg(2,1) + r_hole(3)*bg(3,1) )*fc%nr1t) + 1
nyh = nint ( (r_hole(1)*bg(1,2) + r_hole(2)*bg(2,2) + r_hole(3)*bg(3,2) )*fc%nr2t) + 1
nzh = nint ( (r_hole(1)*bg(1,3) + r_hole(2)*bg(2,3) + r_hole(3)*bg(3,3) )*fc%nr3t) + 1
allocate(v_rh(num_nbndv(1)))
v_rh(:)=0.d0
!get the valence wavefunctions at the nxh,nyh,nzh (only one processor has it!)
#if !defined(__MPI)
nh=(nzh-1)*fc%nrx1t*fc%nrx2t+(nyh-1)*fc%nrx1t+nxh
v_rh(:)=vstate_r%wfnrt(nh,:,1)
#else
nr3s_start=0
nr3s_end =0
do ii=1,mpime+1
nr3s_start=nr3s_end+1
nr3s_end=nr3s_end+fc%dfftt%nr3p(ii)
enddo
do iplane=1,fc%dfftt%my_nr3p
iz=nr3s_start+iplane-1
if (iz==nzh) then
nh=(iplane-1)*fc%nrx1t*fc%nrx2t+(nyh-1)*fc%nrx1t+nxh
v_rh(:)=vstate_r%wfnrt(nh,:,1)
endif
enddo
call mp_sum(v_rh,world_comm)
#endif
if (debug) then
if(ionode) write(stdout,*) 'plotexcwfn qui'
endif
!stop
!allocate and initialize the excitonic wavefunction
allocate(psi_exc(fc%nrxxt))
psi_exc(1:fc%nrxxt)=0.d0
allocate(psi_excsum(fc%nrx1t*fc%nrx2t*fc%nrx3t))
psi_excsum(1:fc%nrx1t*fc%nrx2t*fc%nrx3t)=0.d0
do nplot=nstart,nend
if (debug) then
if(ionode) write(stdout,*) 'plotexcwfn qui2',nplot
endif
!
!FFT the excitonic wavefunction vector to real space (dual grid)
call initialize_exc_r(a_rt)
call fft_a_exc(bse_spectrum(nplot),fc,a_rt)
!now compute the exitonic wavefunction
do iv=1,num_nbndv(1)
psi_exc(1:fc%nrxxt)=psi_exc(1:fc%nrxxt)+v_rh(iv)*&
a_rt%ar(1:a_rt%nrxxt,iv)
enddo
if (debug) then
if(ionode) write(stdout,*) 'plotexcwfn qui3',nplot
endif
!square modulus
psi_exc(1:fc%nrxxt)=psi_exc(1:fc%nrxxt)**2
if(debug) then
if(ionode) write(stdout,*) 'fc%nr1t, fc%nr2t, fc%nr3t', fc%nr1t, fc%nr2t, fc%nr3t
if(ionode) write(stdout,*) 'fc%nrx1t, fc%nrx2t, fc%nrx3t', fc%nrx1t, fc%nrx2t, fc%nrx3t
endif
!Now gather the excitonic wavefunction from all the processors
!and sum for the l_plotaverage case (when nstart is different from nend)
allocate(psi_excio(fc%nrx1t*fc%nrx2t*fc%nrx3t))
psi_excio(1:fc%nrx1t*fc%nrx3t*fc%nrx3t)=0.d0
do iplane=1,fc%dfftt%my_nr3p
iz=nr3s_start+iplane-1
do iy=1,fc%nr2t
do ix=1,fc%nr1t
ii=(iz-1)*(fc%nrx1t*fc%nrx2t)+(iy-1)*fc%nrx1t+ix
iip=(iplane-1)*fc%nrx1t*fc%nrx2t+(iy-1)*fc%nrx1t+ix
psi_excio(ii)=psi_exc(iip)
enddo
enddo
enddo
call mp_sum(psi_excio,world_comm)
psi_excsum(1:fc%nrx1t*fc%nrx3t*fc%nrx3t)=psi_excsum(1:fc%nrx1t*fc%nrx3t*fc%nrx3t)+&
psi_excio(1:fc%nrx1t*fc%nrx3t*fc%nrx3t)/(real(nend)-real(nstart)+1.d0)
if (debug) then
if(ionode) write(stdout,*) 'plotexcwfn qui3',nplot
endif
call free_memory_exc_a_r(a_rt)
deallocate(psi_excio)
enddo ! nplot
!
! XCRYSDEN FORMAT
!
if(ionode) then
ounit=find_free_unit()
! open(ounit,file='exc_average.xsf',form='formatted')
if(l_plotaverage) open(ounit,file='exc_average.xsf',form='formatted')
if(.not.l_plotaverage) then
write(nfile,'(5i1)') &
& nstart/10000,mod(nstart,10000)/1000,mod(nstart,1000)/100,mod(nstart,100)/10,mod(nstart,10)
open(ounit,file=trim(tmp_dir)//trim(prefix)//'.exc.xsf'//nfile,form='formatted')
endif
CALL xsf_struct (alat, at, nat, tau, atm, ityp, ounit)
CALL xsf_fast_datagrid_3d &
(psi_excsum, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, at, alat, ounit)
! close(ounit)
endif
deallocate(v_rh,psi_exc)
deallocate(psi_excsum)
call stop_clock('plot_excwfn')
end subroutine plot_excwfn
|
