Line data Source code
1 : MODULE m_vacden
2 : USE m_juDFT
3 : ! Legacy comments:
4 : ! *************************************************************
5 : ! determines the 2-d star function expansion coefficients of
6 : ! vacuum charge density. speed up by r. wu 1992
7 : ! *************************************************************
8 : !***********************************************************************
9 : ! ****** change vacden(....,q) for vacuum density of states shz Jan.96
10 : ! ****** change vacden(......,vacdos%qstars) for starcoefficients, shz. Jan.99
11 : ! ****** changed for fleur dw
12 : ! In non-collinear calculations the density becomes a hermitian 2x2
13 : ! matrix. This subroutine generates this density matrix in the
14 : ! vacuum region. The diagonal elements of this matrix (n_11 & n_22)
15 : ! are store in den%vacz and den%vacxy, while the real and imaginary part
16 : ! of the off-diagonal element are stored in den%vacz(:,:,3:4) and den%vacxy(:,:,:,3).
17 : !
18 : ! Philipp Kurz 99/07
19 : !***********************************************************************
20 :
21 : !******** ABBREVIATIONS ************************************************
22 : ! qvac : vacuum charge of each eigenstate, needed in in cdnval
23 : ! to determine the vacuum energy parameters
24 : ! vz : non-warping part of the vacuum potential (matrix)
25 : ! collinear : 2. index = ivac (# of vaccum)
26 : ! non-collinear: 2. index = ipot (comp. of pot. matr.)
27 : ! den%vacz : non-warping part of the vacuum density matrix,
28 : ! diagonal elements n_11 and n_22
29 : ! den%vacxy: warping part of the vacuum density matrix,
30 : ! diagonal elements n_11 and n_22
31 : ! den%vacz(:,:,3:4): non-warping part of the vacuum density matrix,
32 : ! off-diagonal elements n_21 (real part in (:,:,3), imaginary part in (:,:,4))
33 : ! den%vacxy(:,:,:,3): warping part of the vacuum density matrix,
34 : ! off-diagonal elements n_21
35 : !***********************************************************************
36 : ! *******************************************************************************
37 : ! layers: no. of layers to be calculated (in vertical direction with z-values as given by izlay)
38 : ! izlay : defines vertical position of layers in delz (=0.1 a.u.) units from begining of vacuum region
39 : ! vacdos: =T: calculate vacuum dos in layers as given by the above
40 : ! integ : =T: integrate in vertical position between izlay(layer,1)..izlay(layer,2)
41 : ! nstm : 0: s-Tip, 1: p_z-Tip, 2: d_z^2-Tip (following Chen's derivative rule) ->rhzgrd.f is used
42 : ! to calculate derivatives
43 : ! tworkf: Workfunction of Tip (in hartree units) is needed for d_z^2-Orbital)
44 : ! starcoeff: =T: star coefficients are calculated at values of izlay for 0th (=q) to nstars-1. star
45 : ! (vacdos%qstars(1..nstars-1))
46 : ! nstars: number of star functions to be used (0th star is given by value of q=charge integrated in 2D)
47 : !
48 : ! further possibility: (readin of locx, locy has to be implemented in flapw7.f or they have to be set explicitly)
49 : !
50 : ! locx and locy can be used to calculate local DOS at a certain vertical position z (or integrated in z)
51 : ! within a restricted area of the 2D unit cell, the corners of this area is given by locx and locy
52 : ! they are defined in internal coordinates, i.e. \vec{r}_1=locx(1)*\vec{a}_1+locy(1)*\vec{a}_2
53 : ! \vec{r}_2=locx(2)*\vec{a}_1+locy(2)*\vec{a}_2
54 : ! \vec{a}_1,2 are the 2D lattice vectors
55 : !
56 : ! **************************************************************************************************
57 : CONTAINS
58 140 : SUBROUTINE vacden(vacuum,stars,input,cell,atoms,noco,nococonv,banddos,&
59 140 : we,ikpt,jspin,vz,ne,ev_list,lapw,evac,den,zMat,vacdos,dos,lapwq,we1,zMat1)
60 : !! Calculates the vacuum part of the density and puts it into den%vac. The variable has
61 : !! four dimensions: The z, star, vacuum and spin index. Recent refactoring combined the
62 : !! real variable vacz and the complex star expansion vacxy into one. vacz is identical
63 : !! to the array for the first star (\(\boldsymbol{G}_{||}=0\)), which was not present in vacxy. The
64 : !! array is bounded by vacuum%nmzd in the first dimension, but only filled up to
65 : !! vacuum%nmzxyd for any star but the first.
66 : !!
67 : !! In practice, the density looks as follows:
68 : !! $$$$
69 : USE m_constants
70 : USE m_grdchlh
71 : USE m_qsf
72 : USE m_cylbes
73 : USE m_dcylbs
74 : USE m_vacuz
75 : USE m_vacudz
76 : USE m_types
77 : USE m_types_vacdos
78 : USE m_types_dos
79 : USE m_npy
80 :
81 : IMPLICIT NONE
82 :
83 : TYPE(t_lapw), INTENT(IN) :: lapw !I just removed the assignment and made lapw INTENT(IN).
84 : TYPE(t_banddos), INTENT(IN) :: banddos
85 : TYPE(t_input), INTENT(IN) :: input
86 : TYPE(t_vacuum), INTENT(IN) :: vacuum
87 : TYPE(t_noco), INTENT(IN) :: noco
88 : TYPE(t_nococonv), INTENT(IN) :: nococonv
89 : TYPE(t_stars), INTENT(IN) :: stars
90 : TYPE(t_cell), INTENT(IN) :: cell
91 : TYPE(t_atoms), INTENT(IN) :: atoms
92 : TYPE(t_mat), INTENT(IN) :: zMat
93 : TYPE(t_potden), INTENT(INOUT) :: den
94 : TYPE(t_vacdos), INTENT(INOUT) :: vacdos
95 : TYPE(t_dos), INTENT(INOUT) :: dos
96 :
97 : INTEGER, INTENT (IN) :: jspin
98 : INTEGER, INTENT (IN) :: ne
99 : INTEGER, INTENT (IN) :: ikpt
100 :
101 : INTEGER, INTENT(IN) :: ev_list(ne)
102 : REAL, INTENT(IN) :: evac(2,input%jspins)
103 : REAL, INTENT(IN) :: we(input%neig)
104 : REAL, INTENT(IN) :: vz(:,:,:) !(vacuum%nmzd,ivac,ispin)
105 :
106 : ! Optional DFPT variables
107 : TYPE(t_lapw), OPTIONAL, INTENT(IN) :: lapwq
108 : TYPE(t_mat), OPTIONAL, INTENT(IN) :: zMat1
109 : REAL, OPTIONAL, INTENT(IN) :: we1(:)
110 :
111 : ! STM/unused arguments
112 : !TYPE(t_gVacMap), INTENT(IN) :: gVacMap ! STM
113 : !REAL, INTENT (IN) :: eig(input%neig)
114 : !TYPE(t_sym), INTENT(IN) :: sym
115 :
116 : COMPLEX :: aa, ab, av, ba, bb, bv, t1, factorx, factory, c_1, tempCmplx
117 : REAL :: arg, const, eps, ev, evacp, phs, phsp, qout, scale, sign, &
118 : uei, uej, ui, uj, wronk, zks, RESULT(1), ui2, uei2, &
119 : k_diff,k_d1,k_d2
120 : INTEGER :: ii, i1, i2, i3, ig3, ig3p, ik, ind2, ind2p, ivac, j, jj, &
121 : jz, k, ikG, ll, ikGPr, n, n2, ispin, kspin, jsp_start, jsp_end, &
122 : ie, isp_start, isp_end, nv2_outer
123 : LOGICAL :: l_dfpt
124 :
125 140 : INTEGER :: nv2(input%jspins)
126 140 : INTEGER :: kvac1(lapw%dim_nv2d(),input%jspins), kvac2(lapw%dim_nv2d(),input%jspins), map2(lapw%dim_nvd(),input%jspins)
127 : REAL :: qssbti(3,2), v(3)
128 :
129 140 : COMPLEX, ALLOCATABLE :: ac(:,:,:), bc(:,:,:), t1jz(:)
130 140 : REAL, ALLOCATABLE :: dt(:), dte(:)
131 140 : REAL, ALLOCATABLE :: t(:), te(:), tei(:,:)
132 140 : REAL, ALLOCATABLE :: u(:,:,:), ue(:,:,:), yy(:)
133 :
134 : ! DFPT stuff:
135 140 : COMPLEX, ALLOCATABLE :: ac1(:,:,:), bc1(:,:,:)
136 140 : REAL, ALLOCATABLE :: dtq(:), dteq(:)
137 140 : REAL, ALLOCATABLE :: tq(:), teq(:), teiq(:,:)
138 140 : REAL, ALLOCATABLE :: uq(:,:,:), ueq(:,:,:)
139 140 : INTEGER :: nv2q(input%jspins)
140 140 : INTEGER :: kvac1q(lapw%dim_nv2d(),input%jspins), kvac2q(lapw%dim_nv2d(),input%jspins), map2q(lapw%dim_nvd(),input%jspins)
141 :
142 : ! local STM/unused variables
143 : !INTEGER :: i,imz,ind1,ind1p,irec2,irec3,m,ipot
144 : !REAL :: ddui,dduj,dduei,dduej,ui_1,uei_1,uj_1,uej_1,wronk_1
145 : !COMPLEX :: aa_1,ab_1,ba_1,bb_1,av_1,bv_1,aae,bbe,abe,bae,aaee,bbee,abee,baee,d
146 : !INTEGER :: mapg2k(lapw%dim_nv2d())
147 : !INTEGER, PARAMETER :: n2max=13
148 : !REAL, PARAMETER :: emax=2.0/hartree_to_ev_const
149 : !REAL, ALLOCATABLE :: du(:),ddu(:,:),due(:),ddue(:,:),dummy(:)
150 :
151 140 : CALL timestart("vacden")
152 :
153 140 : l_dfpt = PRESENT(zMat1)
154 :
155 0 : ALLOCATE(ac(lapw%dim_nv2d(),input%neig,input%jspins), bc(lapw%dim_nv2d(),input%neig,input%jspins), dt(lapw%dim_nv2d()), &
156 : dte(lapw%dim_nv2d()), t(lapw%dim_nv2d()), te(lapw%dim_nv2d()), tei(lapw%dim_nv2d(),input%jspins), &
157 4200 : u(vacuum%nmzd,lapw%dim_nv2d(),input%jspins), ue(vacuum%nmzd,lapw%dim_nv2d(),input%jspins), yy(vacuum%nmzd))
158 140 : IF (l_dfpt) THEN
159 0 : ALLOCATE(ac1(lapw%dim_nv2d(),input%neig,input%jspins), bc1(lapw%dim_nv2d(),input%neig,input%jspins), dtq(lapw%dim_nv2d()), &
160 : dteq(lapw%dim_nv2d()), tq(lapw%dim_nv2d()), teq(lapw%dim_nv2d()), teiq(lapw%dim_nv2d(),input%jspins), &
161 0 : uq(vacuum%nmzd,lapw%dim_nv2d(),input%jspins), ueq(vacuum%nmzd,lapw%dim_nv2d(),input%jspins))
162 : END IF
163 : !ALLOCATE(du(vacuum%nmzd),ddu(vacuum%nmzd,lapw%dim_nv2d()),due(vacuum%nmzd),ddue(vacuum%nmzd,lapw%dim_nv2d()),)
164 :
165 140 : eps=0.01
166 : ! -----> set up mapping arrays
167 140 : IF (noco%l_ss) THEN
168 : jsp_start = 1
169 : jsp_end = 2
170 : ELSE
171 140 : jsp_start = jspin
172 140 : jsp_end = jspin
173 : END IF
174 :
175 280 : DO ispin = jsp_start, jsp_end
176 140 : n2 = 0
177 56166 : k_loop2: DO k = 1, lapw%nv(ispin)
178 1414576 : DO j = 1, n2
179 1414576 : IF (lapw%gvec(1,k,ispin).EQ.kvac1(j,ispin).AND.lapw%gvec(2,k,ispin).EQ.kvac2(j,ispin)) THEN
180 47348 : map2(k,ispin) = j
181 47348 : CYCLE k_loop2
182 : END IF
183 : END DO
184 8678 : n2 = n2 + 1
185 8678 : IF (n2>lapw%dim_nv2d()) CALL juDFT_error("vacden0","vacden")
186 8678 : kvac1(n2,ispin) = lapw%gvec(1,k,ispin)
187 8678 : kvac2(n2,ispin) = lapw%gvec(2,k,ispin)
188 8818 : map2(k,ispin) = n2
189 : END DO k_loop2
190 280 : nv2(ispin) = n2
191 : END DO
192 140 : IF (l_dfpt) THEN
193 0 : DO ispin = jsp_start, jsp_end
194 0 : n2 = 0
195 0 : k_loop2q: DO k = 1, lapwq%nv(ispin)
196 0 : DO j = 1, n2
197 0 : IF (lapwq%gvec(1,k,ispin).EQ.kvac1q(j,ispin).AND.lapwq%gvec(2,k,ispin).EQ.kvac2q(j,ispin)) THEN
198 0 : map2q(k,ispin) = j
199 0 : CYCLE k_loop2q
200 : END IF
201 : END DO
202 0 : n2 = n2 + 1
203 0 : IF (n2>lapw%dim_nv2d()) CALL juDFT_error("vacden0","vacden")
204 0 : kvac1q(n2,ispin) = lapwq%gvec(1,k,ispin)
205 0 : kvac2q(n2,ispin) = lapwq%gvec(2,k,ispin)
206 0 : map2q(k,ispin) = n2
207 : END DO k_loop2q
208 0 : nv2q(ispin) = n2
209 : END DO
210 : END IF
211 :
212 140 : IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
213 : !lapw%nv(2) = lapw%nv(1) ! I just removed this and made lapw INTENT(IN).
214 0 : nv2(2) = nv2(1)
215 0 : DO k = 1,nv2(1)
216 0 : kvac1(k,2) = kvac1(k,1)
217 0 : kvac2(k,2) = kvac2(k,1)
218 : END DO
219 0 : DO k = 1,lapw%nv(1)
220 : !lapw%k3(k,2) = lapw%k3(k,1) ! I just removed this and made lapw INTENT(IN).
221 0 : map2(k,2) = map2(k,1)
222 : END DO
223 : END IF
224 :
225 : ! !+dw
226 : ! ! if tunneling current should be calculated we need to write out
227 : ! ! info on electronic structure: --> mapping from kvac to gvac by mapg2k
228 : ! ! shz, Jan.99
229 : ! IF (.false.) then !vacuum%nstm.EQ.3
230 : ! DO j=1, n2max
231 : ! mapg2k(j)=j
232 : ! DO i=1, nv2(jspin)
233 : ! IF (kvac1(i,jspin).EQ.gVacMap%gvac1d(j).AND.kvac2(i,jspin).EQ.gVacMap%gvac2d(j)) mapg2k(j)=i
234 : ! END DO
235 : ! END DO
236 : ! END IF
237 : ! !
238 : ! !-dw
239 :
240 140 : wronk = 2.0
241 140 : const = 1.0 / ( SQRT(cell%omtil)*wronk )
242 320 : DO ivac = 1,vacuum%nvac
243 460048 : ac(:,:,:) = CMPLX(0.0,0.0)
244 460048 : bc(:,:,:) = CMPLX(0.0,0.0)
245 180 : IF (l_dfpt) THEN
246 0 : ac1(:,:,:) = CMPLX(0.0,0.0)
247 0 : bc1(:,:,:) = CMPLX(0.0,0.0)
248 : END IF
249 180 : sign = 3. - 2.*ivac
250 :
251 180 : IF (noco%l_noco) THEN
252 : !---> In a non-collinear calculation vacden is only called once.
253 : !---> Thus, the vaccum wavefunctions and the A- and B-coeff. (ac bc)
254 : !---> have to be calculated for both spins on that call.
255 : !---> setup the spin-spiral q-vector
256 0 : qssbti(1,1) = - nococonv%qss(1)/2
257 0 : qssbti(2,1) = - nococonv%qss(2)/2
258 0 : qssbti(1,2) = + nococonv%qss(1)/2
259 0 : qssbti(2,2) = + nococonv%qss(2)/2
260 0 : qssbti(3,1) = - nococonv%qss(3)/2
261 0 : qssbti(3,2) = + nococonv%qss(3)/2
262 0 : DO ispin = 1,input%jspins
263 : ! -----> set up vacuum wave functions
264 0 : evacp = evac(ivac,ispin)
265 0 : DO ik = 1,nv2(ispin)
266 0 : v(1) = lapw%bkpt(1) + kvac1(ik,ispin) + qssbti(1,ispin)
267 0 : v(2) = lapw%bkpt(2) + kvac2(ik,ispin) + qssbti(2,ispin)
268 0 : v(3) = 0.
269 :
270 0 : ev = evacp - 0.5*DOT_PRODUCT(v,MATMUL(v,cell%bbmat))
271 :
272 : CALL vacuz(ev,vz(:,ivac,ispin),vz(vacuum%nmz,ivac,ispin),vacuum%nmz,vacuum%delz,t(ik),&
273 0 : dt(ik),u(1,ik,ispin))
274 : CALL vacudz(ev,vz(:,ivac,ispin),vz(vacuum%nmz,ivac,ispin),vacuum%nmz,vacuum%delz,te(ik),&
275 : dte(ik),tei(ik,ispin),ue(1,ik,ispin),dt(ik),&
276 0 : u(1,ik,ispin))
277 :
278 0 : scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
279 0 : te(ik) = scale*te(ik)
280 0 : dte(ik) = scale*dte(ik)
281 0 : tei(ik,ispin) = scale*tei(ik,ispin)
282 0 : DO j = 1,vacuum%nmz
283 0 : ue(j,ik,ispin) = scale*ue(j,ik,ispin)
284 : END DO
285 : END DO
286 : ! -----> construct a and b coefficients
287 0 : DO k = 1,lapw%nv(ispin)
288 : !---> the coefficients of the spin-down basis functions are
289 : !---> stored in the second half of the eigenvector
290 0 : kspin = (lapw%nv(1)+atoms%nlotot)*(ispin-1) + k
291 0 : ikG = map2(k,ispin)
292 0 : zks = lapw%k3(k,ispin)*cell%bmat(3,3)*sign
293 0 : arg = zks*cell%z1
294 0 : c_1 = CMPLX(COS(arg),SIN(arg)) * const
295 0 : av = -c_1 * CMPLX( dte(ikG),zks*te(ikG) )
296 0 : bv = c_1 * CMPLX( dt(ikG),zks* t(ikG) )
297 : ! -----> loop over basis functions
298 0 : IF (zmat%l_real) THEN
299 0 : ac(ikG,:ne,ispin) = ac(ikG,:ne,ispin) + zMat%data_r(kspin,:ne)*av
300 0 : bc(ikG,:ne,ispin) = bc(ikG,:ne,ispin) + zMat%data_r(kspin,:ne)*bv
301 : ELSE
302 0 : ac(ikG,:ne,ispin) = ac(ikG,:ne,ispin) + zMat%data_c(kspin,:ne)*av
303 0 : bc(ikG,:ne,ispin) = bc(ikG,:ne,ispin) + zMat%data_c(kspin,:ne)*bv
304 : END IF
305 : END DO
306 : !---> end of spin loop
307 : END DO
308 : !---> output for testing
309 : ! DO k = 1,10
310 : ! DO n = 1,5
311 : ! DO ispin = 1,jspins
312 : ! write(*,9000)k,n,ispin,ac(k,n,ispin),bc(k,n,ispin)
313 : ! ENDDO
314 : ! ENDDO
315 : ! ENDDO
316 : ! 9000 FORMAT('k=',i3,' ie=',i3,' isp=',i3,
317 : ! + ' ac= (',e12.6,',',e12.6,')',
318 : ! + ' bc= (',e12.6,',',e12.6,')')
319 : ELSE
320 : ! -----> set up vacuum wave functions
321 180 : evacp = evac(ivac,jspin)
322 10938 : DO ik = 1,nv2(jspin)
323 10758 : v(1) = lapw%bkpt(1) + kvac1(ik,jspin)
324 10758 : v(2) = lapw%bkpt(2) + kvac2(ik,jspin)
325 10758 : v(3) = 0.
326 :
327 172128 : ev = evacp - 0.5*DOT_PRODUCT(v,MATMUL(v,cell%bbmat))
328 :
329 10758 : CALL vacuz(ev,vz(:,ivac,jspin),vz(vacuum%nmz,ivac,jspin),vacuum%nmz,vacuum%delz,t(ik),dt(ik),u(1,ik,jspin))
330 : CALL vacudz(ev,vz(:,ivac,jspin),vz(vacuum%nmz,ivac,jspin),vacuum%nmz,vacuum%delz,te(ik),&
331 10758 : dte(ik),tei(ik,jspin),ue(1,ik,jspin),dt(ik),u(1,ik,jspin))
332 :
333 10758 : scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
334 10758 : te(ik) = scale*te(ik)
335 10758 : dte(ik) = scale*dte(ik)
336 10758 : tei(ik,jspin) = scale*tei(ik,jspin)
337 2700438 : DO j = 1,vacuum%nmz
338 2700258 : ue(j,ik,jspin) = scale*ue(j,ik,jspin)
339 : END DO
340 : END DO
341 180 : IF (l_dfpt) THEN
342 0 : DO ik = 1,nv2q(jspin)
343 0 : v(1) = lapwq%bkpt(1) + kvac1q(ik,jspin) + lapwq%qphon(1)
344 0 : v(2) = lapwq%bkpt(2) + kvac2q(ik,jspin) + lapwq%qphon(2)
345 0 : v(3) = 0.
346 :
347 0 : ev = evacp - 0.5*DOT_PRODUCT(v,MATMUL(v,cell%bbmat))
348 :
349 0 : CALL vacuz(ev,vz(:,ivac,jspin),vz(vacuum%nmz,ivac,jspin),vacuum%nmz,vacuum%delz,tq(ik),dtq(ik),uq(1,ik,jspin))
350 : CALL vacudz(ev,vz(:,ivac,jspin),vz(vacuum%nmz,ivac,jspin),vacuum%nmz,vacuum%delz,teq(ik),&
351 0 : dteq(ik),teiq(ik,jspin),ueq(1,ik,jspin),dtq(ik),uq(1,ik,jspin))
352 :
353 0 : scale = wronk/ (teq(ik)*dtq(ik)-dteq(ik)*tq(ik))
354 0 : teq(ik) = scale*teq(ik)
355 0 : dteq(ik) = scale*dteq(ik)
356 0 : teiq(ik,jspin) = scale*teiq(ik,jspin)
357 0 : DO j = 1,vacuum%nmz
358 0 : ueq(j,ik,jspin) = scale*ueq(j,ik,jspin)
359 : END DO
360 : END DO
361 : END IF
362 : ! -----> construct a and b coefficients
363 75006 : DO k = 1,lapw%nv(jspin)
364 74826 : ikG = map2(k,jspin)
365 74826 : zks = lapw%k3(k,jspin)*cell%bmat(3,3)*sign
366 74826 : arg = zks*cell%z1
367 74826 : c_1 = CMPLX(COS(arg),SIN(arg)) * const
368 74826 : av = -c_1 * CMPLX( dte(ikG),zks*te(ikG) )
369 74826 : bv = c_1 * CMPLX( dt(ikG),zks* t(ikG) )
370 : ! -----> loop over basis functions
371 75006 : IF (zmat%l_real) THEN
372 645575 : ac(ikG,:ne,jspin) = ac(ikG,:ne,jspin) + zMat%data_r(k,:ne)*av
373 645575 : bc(ikG,:ne,jspin) = bc(ikG,:ne,jspin) + zMat%data_r(k,:ne)*bv
374 : ELSE
375 600904 : ac(ikG,:ne,jspin) = ac(ikG,:ne,jspin) + zMat%data_c(k,:ne)*av
376 600904 : bc(ikG,:ne,jspin) = bc(ikG,:ne,jspin) + zMat%data_c(k,:ne)*bv
377 : END IF
378 : END DO
379 180 : IF (l_dfpt) THEN
380 0 : DO k = 1,lapwq%nv(jspin)
381 0 : ikG = map2q(k,jspin)
382 0 : zks = lapwq%k3(k,jspin)*cell%bmat(3,3)*sign
383 0 : arg = zks*cell%z1
384 0 : c_1 = CMPLX(COS(arg),SIN(arg)) * const
385 0 : av = -c_1 * CMPLX( dteq(ikG),zks*teq(ikG) )
386 0 : bv = c_1 * CMPLX( dtq(ikG),zks* tq(ikG) )
387 : ! -----> loop over basis functions
388 0 : ac1(ikG,:ne,jspin) = ac1(ikG,:ne,jspin) + 2*zMat1%data_c(k,:ne)*av
389 0 : bc1(ikG,:ne,jspin) = bc1(ikG,:ne,jspin) + 2*zMat1%data_c(k,:ne)*bv
390 : END DO
391 : END IF
392 : END IF
393 : !
394 : ! ----> calculate first and second derivative of u,ue
395 : ! in order to simulate p_z or d_z^2 Tip in Chen's model , shz. 97
396 : !
397 : ! IF (.false.) THEN !vacuum%nstm.GT.0
398 : ! DO ik = 1,nv2(jspin)
399 : ! ! CALL rhzgrd(nmz,delz,u(1,ik,jspin),4,du,ddu(1,ik))
400 : ! ! CALL rhzgrd(nmz,delz,ue(1,ik,jspin),4,due,ddue(1,ik))
401 : !
402 : ! ALLOCATE ( dummy(vacuum%nmz) )
403 : ! CALL grdchlh(&
404 : ! vacuum%delz,u(1:vacuum%nmz,ik,jspin),&
405 : ! du,ddu(:,ik),order=4)
406 : ! CALL grdchlh(&
407 : ! vacuum%delz,ue(1:vacuum%nmz,ik,jspin),&
408 : ! due,ddue(:,ik),order=4)
409 : ! DEALLOCATE ( dummy )
410 : !
411 : ! IF (.FALSE.) THEN !IF (vacuum%nstm.EQ.1) THEN
412 : ! u(:vacuum%nmz,ik,jspin)=du(:vacuum%nmz)
413 : ! ue(:vacuum%nmz,ik,jspin)=due(:vacuum%nmz)
414 : ! END IF
415 : ! ENDDO
416 : ! END IF
417 :
418 : !+dw
419 :
420 : !
421 : ! --> to calculate Tunneling Current between two systems
422 : ! within Bardeens Approach one needs ac(l,n), bc(l,n);
423 : ! they are written to the file vacwave
424 : ! IF nstm=3
425 : ! tworkf is then the fermi energy (in hartree)
426 : !
427 : ! IF (.false.)then !(vacuum%nstm.EQ.3) THEN
428 : !#ifdef CPP_MPI
429 : ! CALL judft_error("nstm==3 does not work in parallel",calledby="vacden")
430 : !#else
431 : ! i=0
432 : ! DO n = 1, ne
433 : ! IF (ABS(eig(n)-banddos%tworkf).LE.emax) i=i+1
434 : ! END DO
435 : ! WRITE (87,FMT=990) lapw%bkpt(1),lapw%bkpt(2), i, n2max
436 : ! DO n = 1, ne
437 : ! IF (ABS(eig(n)-banddos%tworkf).LE.emax) THEN
438 : ! WRITE (87,FMT=1000) eig(n)
439 : ! DO j=1,n2max
440 : ! WRITE (87,FMT=1010) ac(mapg2k(j),n,jspin),&
441 : ! bc(mapg2k(j),n,jspin)
442 : ! END DO
443 : ! END IF
444 : ! END DO
445 : !#endif
446 : ! END IF
447 : !990 FORMAT(2(f8.4,1x),i3,1x,i3)
448 : !1000 FORMAT(e12.4)
449 : !1010 FORMAT(2(2e20.8,1x))
450 : !
451 : ! ------------------------------------------------------------
452 : !-dw
453 : !
454 : !----> non-warping part of the density (g=0 star)
455 : !
456 : ! IF (.false.) then !vacuum%nstm.EQ.2) THEN
457 : ! !
458 : ! ! ----> d_z^2-Tip needs: |d^2(psi)/dz^2 - kappa^2/3 psi|^2
459 : ! !
460 : ! DO l = 1,nv2(jspin)
461 : ! aa = 0.0
462 : ! bb = 0.0
463 : ! ba = 0.0
464 : ! ab = 0.0
465 : ! DO n = 1,ne
466 : ! aa = aa + we(n)*CONJG(ac(l,n,jspin))*ac(l,n,jspin)
467 : ! bb = bb + we(n)*CONJG(bc(l,n,jspin))*bc(l,n,jspin)
468 : ! ab = ab + we(n)*CONJG(ac(l,n,jspin))*bc(l,n,jspin)
469 : ! ba = ba + we(n)*CONJG(bc(l,n,jspin))*ac(l,n,jspin)
470 : ! qout = REAL(CONJG(ac(l,n,jspin))*ac(l,n,jspin)+tei(l,jspin)*CONJG(bc(l,n,jspin))*bc(l,n,jspin))
471 : ! vacdos%qvac(ev_list(n),ivac,ikpt,jspin) = vacdos%qvac(ev_list(n),ivac,ikpt,jspin) + qout*cell%area
472 : ! dos%qTot(ev_list(n),ikpt,jspin) = dos%qTot(ev_list(n),ikpt,jspin) + qout*cell%area
473 : ! END DO
474 : ! aae=-aa*banddos%tworkf*2/3
475 : ! bbe=-bb*banddos%tworkf*2/3
476 : ! abe=-ab*banddos%tworkf*2/3
477 : ! bae=-ba*banddos%tworkf*2/3
478 : ! aaee=aa*banddos%tworkf*banddos%tworkf*4/9
479 : ! bbee=bb*banddos%tworkf*banddos%tworkf*4/9
480 : ! abee=ab*banddos%tworkf*banddos%tworkf*4/9
481 : ! baee=ba*banddos%tworkf*banddos%tworkf*4/9
482 : ! DO jz = 1,vacuum%nmz
483 : ! ui = u(jz,l,jspin)
484 : ! uei = ue(jz,l,jspin)
485 : ! ddui = ddu(jz,l)
486 : ! dduei = ddue(jz,l)
487 : ! den%vacz(jz,ivac,jspin) = den%vacz(jz,ivac,jspin) +&
488 : ! REAL(aaee*ui*ui+bbee*uei*uei+&
489 : ! (abee+baee)*ui*uei+aa*ddui*ddui+&
490 : ! bb*dduei*dduei+(ab+ba)*ddui*dduei+&
491 : ! 2*aae*ui*ddui+2*bbe*uei*dduei+&
492 : ! (abe+bae)*(ui*dduei+uei*ddui))
493 : !
494 : ! ENDDO
495 : ! END DO
496 : ! !
497 : ! ! -----> s-Tip: |psi|^2 and p-Tip: |d(psi)/dz|^2
498 : ! !
499 : ! ELSE
500 180 : IF (noco%l_noco) THEN
501 : !---> diagonal elements of the density matrix, n_11 and n_22
502 : !---> the non-warping part of n_21 is calculated together with
503 : !---> the warping part of n_21
504 0 : DO ispin = 1,input%jspins
505 0 : DO ikG = 1,nv2(ispin)
506 : aa = 0.0
507 : bb = 0.0
508 : ba = 0.0
509 : ab = 0.0
510 0 : DO n = 1,ne
511 0 : aa=aa + we(n)*CONJG(ac(ikG,n,ispin))*ac(ikG,n,ispin)
512 0 : bb=bb + we(n)*CONJG(bc(ikG,n,ispin))*bc(ikG,n,ispin)
513 0 : ab=ab + we(n)*CONJG(ac(ikG,n,ispin))*bc(ikG,n,ispin)
514 0 : ba=ba + we(n)*CONJG(bc(ikG,n,ispin))*ac(ikG,n,ispin)
515 0 : qout = REAL(CONJG(ac(ikG,n,ispin))*ac(ikG,n,ispin)+tei(ikG,ispin)*CONJG(bc(ikG,n,ispin))*bc(ikG,n,ispin))
516 0 : vacdos%qvac(ev_list(n),ivac,ikpt,ispin) = vacdos%qvac(ev_list(n),ivac,ikpt,ispin) + qout*cell%area
517 0 : dos%qTot(ev_list(n),ikpt,ispin) = dos%qTot(ev_list(n),ikpt,ispin) + qout*cell%area
518 : END DO
519 0 : DO jz = 1,vacuum%nmz
520 0 : ui = u(jz,ikG,ispin)
521 0 : uei = ue(jz,ikG,ispin)
522 0 : den%vac(jz,1,ivac,ispin) = den%vac(jz,1,ivac,ispin) + REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei) ! TODO: AN TB; sollte man das REAL killen?
523 : END DO
524 : END DO
525 : END DO
526 720 : ELSE IF ((.NOT.l_dfpt).OR.(norm2(stars%center)<1e-8)) THEN
527 10938 : DO ikG = 1,nv2(jspin)
528 : aa = CMPLX(0.0,0.0)
529 : bb = CMPLX(0.0,0.0)
530 : ba = CMPLX(0.0,0.0)
531 : ab = CMPLX(0.0,0.0)
532 186571 : DO n = 1,ne
533 175813 : IF (.NOT.l_dfpt) THEN
534 175813 : aa = aa + we(n)*CONJG(ac(ikG,n,jspin))*ac(ikG,n,jspin)
535 175813 : bb = bb + we(n)*CONJG(bc(ikG,n,jspin))*bc(ikG,n,jspin)
536 175813 : ab = ab + we(n)*CONJG(ac(ikG,n,jspin))*bc(ikG,n,jspin)
537 175813 : ba = ba + we(n)*CONJG(bc(ikG,n,jspin))*ac(ikG,n,jspin)
538 : ELSE
539 0 : aa = aa + we1(n)*CONJG(ac(ikG,n,jspin))*ac(ikG,n,jspin)
540 0 : bb = bb + we1(n)*CONJG(bc(ikG,n,jspin))*bc(ikG,n,jspin)
541 0 : ab = ab + we1(n)*CONJG(ac(ikG,n,jspin))*bc(ikG,n,jspin)
542 0 : ba = ba + we1(n)*CONJG(bc(ikG,n,jspin))*ac(ikG,n,jspin)
543 0 : aa = aa + we(n)*CONJG(ac(ikG,n,jspin))*ac1(ikG,n,jspin)
544 0 : bb = bb + we(n)*CONJG(bc(ikG,n,jspin))*bc1(ikG,n,jspin)
545 0 : ab = ab + we(n)*CONJG(ac(ikG,n,jspin))*bc1(ikG,n,jspin)
546 0 : ba = ba + we(n)*CONJG(bc(ikG,n,jspin))*ac1(ikG,n,jspin)
547 : END IF
548 175813 : qout = REAL(CONJG(ac(ikG,n,jspin))*ac(ikG,n,jspin)+tei(ikG,jspin)*CONJG(bc(ikG,n,jspin))*bc(ikG,n,jspin))
549 175813 : vacdos%qvac(ev_list(n),ivac,ikpt,jspin) = vacdos%qvac(ev_list(n),ivac,ikpt,jspin) + qout*cell%area
550 186571 : dos%qTot(ev_list(n),ikpt,jspin) = dos%qTot(ev_list(n),ikpt,jspin) + qout*cell%area
551 : END DO
552 2700438 : DO jz = 1,vacuum%nmz
553 2689500 : ui = u(jz,ikG,jspin)
554 2689500 : uei = ue(jz,ikG,jspin)
555 2700258 : IF (.NOT.l_dfpt) THEN
556 2689500 : den%vac(jz,1,ivac,jspin) = den%vac(jz,1,ivac,jspin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei) ! TODO: REAL weg?
557 : ELSE
558 0 : den%vac(jz,1,ivac,jspin) = den%vac(jz,1,ivac,jspin) + aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei
559 : END IF
560 : END DO
561 : END DO
562 : END IF
563 : ! END IF
564 : !
565 : ! ****************** change for vacuum density of states shz Jan.96 ***
566 : !
567 180 : IF (banddos%vacdos) THEN
568 : !
569 : ! ----> d_z^2-Tip needs: |d^2(psi)/dz^2 - kappa^2/3 psi|^2
570 : !
571 : ! IF (.false.) THEN !IF (vacuum%nstm.EQ.2) THEN
572 : ! DO l=1,nv2(jspin)
573 : ! DO n = 1,ne
574 : ! aa = CONJG(ac(l,n,jspin))*ac(l,n,jspin)
575 : ! bb = CONJG(bc(l,n,jspin))*bc(l,n,jspin)
576 : ! ab = CONJG(ac(l,n,jspin))*bc(l,n,jspin)
577 : ! ba = CONJG(bc(l,n,jspin))*ac(l,n,jspin)
578 : ! aae = -banddos%tworkf*aa*2/3
579 : ! bbe = -banddos%tworkf*bb*2/3
580 : ! abe = -banddos%tworkf*ab*2/3
581 : ! bae = -banddos%tworkf*ba*2/3
582 : ! aaee = aa*banddos%tworkf*banddos%tworkf*4/9
583 : ! bbee = bb*banddos%tworkf*banddos%tworkf*4/9
584 : ! abee = ab*banddos%tworkf*banddos%tworkf*4/9
585 : ! baee = ba*banddos%tworkf*banddos%tworkf*4/9
586 : ! DO jj = 1,banddos%layers
587 : ! !
588 : ! ! ----> either integrated LDOS(z1,z2) or LDOS(z1)
589 : ! !
590 : ! IF (input%integ) THEN
591 : ! ll = 1
592 : ! DO ii = banddos%izlay(jj,1),banddos%izlay(jj,2)
593 : ! ui = u(ii,l,jspin)
594 : ! uei = ue(ii,l,jspin)
595 : ! ddui = ddu(ii,l)
596 : ! dduei = ddue(ii,l)
597 : ! yy(ll) = REAL(aaee*ui*ui+bbee*uei*uei+(abee+baee)*ui*uei+aa*ddui*ddui+bb*&
598 : ! dduei*dduei+(ab+ba)*ddui*dduei+2*aae*ui*ddui+2*bbe*uei*dduei+&
599 : ! (abe+bae)*(ui*dduei+uei*ddui))*cell%area
600 : ! ll = ll+1
601 : ! END DO
602 : ! CALL qsf(vacuum%delz,yy,RESULT,ll-1,0)
603 : ! vacdos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) = vacdos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) + RESULT(1)
604 : ! ELSE
605 : ! ui = u(banddos%izlay(jj,1),l,jspin)
606 : ! uei = ue(banddos%izlay(jj,1),l,jspin)
607 : ! ddui = ddu(banddos%izlay(jj,1),l)
608 : ! dduei = ddue(banddos%izlay(jj,1),l)
609 : ! yy(1) = REAL(aaee*ui*ui+bbee*uei*uei+&
610 : ! (abee+baee)*ui*uei+aa*ddui*ddui+&
611 : ! bb*dduei*dduei+(ab+ba)*ddui*dduei+&
612 : ! 2*aae*ui*ddui+2*bbe*uei*dduei+&
613 : ! (abe+bae)*(ui*dduei+uei*ddui))
614 : ! vacdos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) = vacdos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) +yy (1)
615 : ! END IF
616 : ! END DO
617 : ! END DO
618 : ! END DO
619 : ! !
620 : ! ! ----> s-Tip = calculate LDOS and(!) p_z-Tip (since u->du/dz, ue->due/dz)
621 : ! !
622 : ! ELSE
623 0 : IF (ABS(banddos%locx(1)-banddos%locx(2)).LE.eps) THEN
624 : !----> integrated over 2D-unit cell
625 0 : IF (noco%l_noco) THEN
626 : isp_start = 1
627 : isp_end = input%jspins
628 : ELSE
629 0 : isp_start = jspin
630 0 : isp_end = jspin
631 : END IF
632 0 : DO ispin = isp_start, isp_end
633 0 : DO ikG=1,nv2(ispin)
634 0 : DO n = 1,ne
635 0 : aa = CONJG(ac(ikG,n,ispin))*ac(ikG,n,ispin)
636 0 : bb = CONJG(bc(ikG,n,ispin))*bc(ikG,n,ispin)
637 0 : ab = CONJG(ac(ikG,n,ispin))*bc(ikG,n,ispin)
638 0 : ba = CONJG(bc(ikG,n,ispin))*ac(ikG,n,ispin)
639 0 : DO jj = 1,banddos%layers
640 : !---> either integrated (z1,z2) or slice (z1)
641 0 : IF (input%integ) THEN
642 0 : ll = 1
643 0 : DO ii = banddos%izlay(jj,1),banddos%izlay(jj,2)
644 0 : ui = u(ii,ikG,ispin)
645 0 : uei = ue(ii,ikG,ispin)
646 0 : yy(ll) = REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
647 0 : ll = ll+1
648 : END DO
649 0 : CALL qsf(vacuum%delz,yy,RESULT,ll-1,0)
650 0 : vacdos%qvlay(ev_list(n),jj,ivac,ikpt,ispin) = vacdos%qvlay(ev_list(n),jj,ivac,ikpt,ispin) + RESULT(1)
651 : ELSE
652 0 : ui = u(banddos%izlay(jj,1),ikG,ispin)
653 0 : uei = ue(banddos%izlay(jj,1),ikG,ispin)
654 : vacdos%qvlay(ev_list(n),jj,ivac,ikpt,ispin) = vacdos%qvlay(ev_list(n),jj,ivac,ikpt,ispin) &
655 0 : + REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
656 :
657 : END IF
658 : END DO
659 : END DO
660 : END DO
661 : END DO
662 : ELSE
663 : !----> if LDOS should be calculated over restricted area of the 2D-unit cell
664 : ! lower left corner: (locx(1), locy(1)) } in internal
665 : ! upper right corner: (locx(2), locy(2)) } coordinates
666 : !
667 0 : DO ikG=1, nv2(jspin)
668 0 : DO ikGPr=1, nv2(jspin)
669 0 : IF (kvac1(ikG,jspin).EQ.kvac1(ikGPr,jspin)) THEN
670 0 : factorx = CMPLX((banddos%locx(2)-banddos%locx(1)), 0.)
671 : ELSE
672 0 : k_diff=tpi_const*(kvac1(ikG,jspin)-kvac1(ikGPr,jspin))
673 0 : k_d1 = k_diff*banddos%locx(1)
674 0 : k_d2 = k_diff*banddos%locx(2)
675 : factorx=(CMPLX( COS(k_d2), SIN(k_d2)) - &
676 : CMPLX( COS(k_d1), SIN(k_d1)) ) / &
677 0 : CMPLX( 0.,k_diff )
678 : END IF
679 0 : IF (kvac2(ikG,jspin).EQ.kvac2(ikGPr,jspin)) THEN
680 0 : factory = CMPLX((banddos%locy(2)-banddos%locy(1)), 0.)
681 : ELSE
682 0 : k_diff=tpi_const*(kvac2(ikG,jspin)-kvac2(ikGPr,jspin))
683 0 : k_d1 = k_diff*banddos%locy(1)
684 0 : k_d2 = k_diff*banddos%locy(2)
685 : factory=(CMPLX( COS(k_d2), SIN(k_d2)) - &
686 : CMPLX( COS(k_d1), SIN(k_d1)) ) / &
687 0 : CMPLX( 0.,k_diff )
688 : END IF
689 0 : DO n=1, ne
690 0 : aa = CONJG(ac(ikGPr,n,jspin))*ac(ikG,n,jspin)
691 0 : bb = CONJG(bc(ikGPr,n,jspin))*bc(ikG,n,jspin)
692 0 : ab = CONJG(ac(ikGPr,n,jspin))*bc(ikG,n,jspin)
693 0 : ba = CONJG(bc(ikGPr,n,jspin))*ac(ikG,n,jspin)
694 0 : DO jj = 1,banddos%layers
695 : !---> either integrated (z1,z2) or slice (z1)
696 0 : IF (input%integ) THEN
697 0 : ll = 1
698 0 : DO ii = banddos%izlay(jj,1), banddos%izlay(jj,2)
699 0 : ui = u(ii,ikG,jspin)
700 0 : uei = ue(ii,ikG,jspin)
701 0 : uj = u(ii,ikGPr,jspin)
702 0 : uej = ue(ii,ikGPr,jspin)
703 0 : yy(ll) = REAL((aa*ui*uj+bb*uei*uej+ab*uei*uj+ba*ui*uej)*factorx*factory)
704 0 : ll = ll+1
705 : END DO
706 0 : CALL qsf(vacuum%delz,yy,RESULT,ll-1,0)
707 0 : vacdos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) = vacdos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) + RESULT(1)
708 : ELSE
709 0 : ui = u(banddos%izlay(jj,1),ikG,jspin)
710 0 : uei = ue(banddos%izlay(jj,1),ikG,jspin)
711 0 : uj = u(banddos%izlay(jj,1),ikGPr,jspin)
712 0 : uej = ue(banddos%izlay(jj,1),ikGPr,jspin)
713 0 : vacdos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) = REAL((aa*ui*uj + bb*uei*uej+ab*uei*uj+ba*ui**uej)*factorx*factory)
714 : END IF
715 : END DO
716 : END DO
717 : END DO
718 : END DO
719 : END IF
720 : !END IF
721 : END IF
722 : !
723 : ! **********************************************************************
724 : !
725 : !---> warping part of the density (g.ne.0 stars)
726 : !
727 : ! ---> d_z^2-Tip
728 : !
729 : ! if (.false.) then !IF (vacuum%nstm.EQ.2) THEN
730 : ! DO l = 1,nv2(jspin)
731 : ! DO l1 = 1,l - 1
732 : ! i1 = kvac1(l,jspin) - kvac1(l1,jspin)
733 : ! i2 = kvac2(l,jspin) - kvac2(l1,jspin)
734 : ! i3 = 0
735 : ! IF (iabs(i1).GT.stars%mx1) CYCLE
736 : ! IF (iabs(i2).GT.stars%mx2) CYCLE
737 : ! ig3 = stars%ig(i1,i2,i3)
738 : ! IF (ig3.EQ.0) CYCLE
739 : ! phs = stars%rgphs(i1,i2,i3)
740 : ! ig3p = stars%ig(-i1,-i2,i3)
741 : ! phsp = stars%rgphs(-i1,-i2,i3)
742 : ! ind2 = stars%ig2(ig3)
743 : ! ind2p = stars%ig2(ig3p)
744 : ! aa = 0.0
745 : ! bb = 0.0
746 : ! ba = 0.0
747 : ! ab = 0.0
748 : ! DO n = 1,ne
749 : ! aa = aa + we(n)*CONJG(ac(l1,n,jspin))*ac(l,n,jspin)
750 : ! bb = bb + we(n)*CONJG(bc(l1,n,jspin))*bc(l,n,jspin)
751 : ! ab = ab + we(n)*CONJG(ac(l1,n,jspin))*bc(l,n,jspin)
752 : ! ba = ba + we(n)*CONJG(bc(l1,n,jspin))*ac(l,n,jspin)
753 : ! END DO
754 : ! aae=-aa*2/3*banddos%tworkf
755 : ! bbe=-bb*2/3*banddos%tworkf
756 : ! abe=-ab*2/3*banddos%tworkf
757 : ! bae=-ba*2/3*banddos%tworkf
758 : ! aaee=aa*4/9*banddos%tworkf*banddos%tworkf
759 : ! bbee=bb*4/9*banddos%tworkf*banddos%tworkf
760 : ! abee=ab*4/9*banddos%tworkf*banddos%tworkf
761 : ! baee=ba*4/9*banddos%tworkf*banddos%tworkf
762 : ! DO jz = 1,vacuum%nmzxy
763 : ! ui = u(jz,l,jspin)
764 : ! uj = u(jz,l1,jspin)
765 : ! ddui = ddu(jz,l)
766 : ! dduj = ddu(jz,l1)
767 : ! uei = ue(jz,l,jspin)
768 : ! uej = ue(jz,l1,jspin)
769 : ! dduei = ddue(jz,l)
770 : ! dduej = ddue(jz,l1)
771 : ! t1=aaee*ui*uj+bbee*uei*uej+baee*ui*uej+abee*uei*uj&
772 : ! + aae*(ui*dduj+uj*ddui)+bbe*(uei*dduej+uej*dduei)&
773 : ! + abe*(ui*dduej+uj*dduei)+bae*(ddui*uej+dduj*uei)&
774 : ! + aa*ddui*dduj+bb*dduei*dduej+ba*ddui*dduej&
775 : ! + ab*dduei*dduj
776 : ! den%vacxy(jz,ind2-1,ivac,jspin) = den%vacxy(jz,ind2-1,ivac,jspin) + t1*phs/stars%nstr2(ind2)
777 : ! den%vacxy(jz,ind2p-1,ivac,jspin)= den%vacxy(jz,ind2p-1,ivac,jspin) + CONJG(t1)*phsp/stars%nstr2(ind2p)
778 : ! ENDDO
779 : ! ENDDO
780 : ! END DO
781 : ! !
782 : ! ! ---> s-Tip and p_z-Tip
783 : ! !
784 : ! ELSE
785 : !=============================================================
786 : ! continuation of vacden....
787 : !=============================================================
788 180 : IF (noco%l_noco) THEN
789 : !---> diagonal elements of the density matrix, n_11 and n_22
790 0 : CALL timestart("vacden4_noco")
791 0 : DO ispin = 1,input%jspins
792 0 : DO ikG = 1,nv2(ispin)
793 0 : DO ikGPr = 1, ikG - 1
794 0 : i1 = kvac1(ikG,ispin) - kvac1(ikGPr,ispin)
795 0 : i2 = kvac2(ikG,ispin) - kvac2(ikGPr,ispin)
796 0 : i3 = 0
797 0 : IF (iabs(i1).GT.stars%mx1) CYCLE
798 0 : IF (iabs(i2).GT.stars%mx2) CYCLE
799 0 : ig3 = stars%ig(i1,i2,i3)
800 0 : IF (ig3.EQ.0) CYCLE
801 0 : phs = stars%rgphs(i1,i2,i3)
802 0 : ig3p = stars%ig(-i1,-i2,i3)
803 0 : phsp = stars%rgphs(-i1,-i2,i3)
804 0 : ind2 = stars%ig2(ig3)
805 0 : ind2p = stars%ig2(ig3p)
806 0 : aa = 0.0
807 0 : bb = 0.0
808 0 : ba = 0.0
809 0 : ab = 0.0
810 0 : DO n = 1,ne
811 0 : aa=aa+we(n)*CONJG(ac(ikGPr,n,ispin))*ac(ikG,n,ispin)
812 0 : bb=bb+we(n)*CONJG(bc(ikGPr,n,ispin))*bc(ikG,n,ispin)
813 0 : ab=ab+we(n)*CONJG(ac(ikGPr,n,ispin))*bc(ikG,n,ispin)
814 0 : ba=ba+we(n)*CONJG(bc(ikGPr,n,ispin))*ac(ikG,n,ispin)
815 : END DO
816 0 : DO jz = 1,vacuum%nmzxy
817 0 : ui = u(jz,ikG,ispin)
818 0 : uj = u(jz,ikGPr,ispin)
819 0 : uei = ue(jz,ikG,ispin)
820 0 : uej = ue(jz,ikGPr,ispin)
821 0 : t1 = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
822 0 : den%vac(jz,ind2,ivac,ispin) = den%vac(jz,ind2,ivac,ispin) + t1*phs/stars%nstr2(ind2)
823 0 : den%vac(jz,ind2p,ivac,ispin) = den%vac(jz,ind2p,ivac,ispin) + CONJG(t1)*phsp/stars%nstr2(ind2p)
824 : END DO
825 : END DO
826 : END DO
827 : END DO
828 : !---> off-diagonal element of the density matrix, n_21
829 0 : DO ikG = 1,nv2(1)
830 0 : DO ikGPr = 1,nv2(2)
831 0 : i1 = kvac1(ikG,1) - kvac1(ikGPr,2)
832 0 : i2 = kvac2(ikG,1) - kvac2(ikGPr,2)
833 0 : i3 = 0
834 : !---> treat only the warping part
835 0 : IF (iabs(i1).GT.stars%mx1) CYCLE
836 0 : IF (iabs(i2).GT.stars%mx2) CYCLE
837 0 : ig3 = stars%ig(i1,i2,i3)
838 0 : IF (ig3.EQ.0) CYCLE
839 0 : phs = stars%rgphs(i1,i2,i3)
840 0 : ind2 = stars%ig2(ig3)
841 0 : IF ( ind2.EQ.1) THEN
842 : !---> non-warping part (1st star G=0)
843 : aa = 0.0
844 : bb = 0.0
845 : ba = 0.0
846 : ab = 0.0
847 0 : DO ie = 1,ne
848 0 : aa=aa+we(ie)*CONJG(ac(ikGPr,ie,2))*ac(ikG,ie,1)
849 0 : bb=bb+we(ie)*CONJG(bc(ikGPr,ie,2))*bc(ikG,ie,1)
850 0 : ab=ab+we(ie)*CONJG(ac(ikGPr,ie,2))*bc(ikG,ie,1)
851 0 : ba=ba+we(ie)*CONJG(bc(ikGPr,ie,2))*ac(ikG,ie,1)
852 : END DO
853 0 : DO jz = 1,vacuum%nmz
854 0 : ui = u(jz,ikG,1)
855 0 : ui2 = u(jz,ikGPr,2)
856 0 : uei = ue(jz,ikG,1)
857 0 : uei2 = ue(jz,ikGPr,2)
858 0 : tempCmplx = aa*ui2*ui + bb*uei2*uei + ab*ui2*uei + ba*uei2*ui
859 0 : den%vac(jz,1,ivac,3) = den%vac(jz,1,ivac,3) + CONJG(tempCmplx) !!!! MAGIC MINUS
860 : END DO
861 : ELSE
862 : !---> warping part
863 : aa = 0.0
864 : bb = 0.0
865 : ba = 0.0
866 : ab = 0.0
867 0 : DO ie = 1,ne
868 0 : aa=aa + we(ie)*CONJG(ac(ikGPr,ie,2))*ac(ikG,ie,1)
869 0 : bb=bb + we(ie)*CONJG(bc(ikGPr,ie,2))*bc(ikG,ie,1)
870 0 : ab=ab + we(ie)*CONJG(ac(ikGPr,ie,2))*bc(ikG,ie,1)
871 0 : ba=ba + we(ie)*CONJG(bc(ikGPr,ie,2))*ac(ikG,ie,1)
872 : END DO
873 0 : DO jz = 1,vacuum%nmzxy
874 0 : ui = u(jz,ikG,1)
875 0 : uj = u(jz,ikGPr,2)
876 0 : uei = ue(jz,ikG,1)
877 0 : uej = ue(jz,ikGPr,2)
878 0 : t1 = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
879 0 : den%vac(jz,ind2,ivac,3) = den%vac(jz, ind2,ivac,3) + conjg(t1*phs/stars%nstr2(ind2))
880 : END DO
881 : END IF
882 : END DO
883 : END DO
884 0 : CALL timestop("vacden4_noco")
885 : ELSE ! collinear part
886 180 : IF (.NOT.l_dfpt) THEN
887 180 : nv2_outer = nv2(jspin)
888 : ELSE
889 0 : nv2_outer = nv2q(jspin)
890 : END IF
891 : !$OMP PARALLEL DEFAULT(none) &
892 : !$OMP SHARED(nv2,nv2_outer,jspin,kvac1,kvac2,kvac1q,kvac2q,stars,ne,we,we1,vacuum,den,ac,bc,ac1,bc1,u,ue,uq,ueq,ivac,l_dfpt) &
893 : !$OMP PRIVATE(ikGPr,i1,i2,i3,ig3,phs,ig3p,phsp,ind2,ind2p,n,jz,ui,uj,uei,uej)&
894 180 : !$OMP PRIVATE(aa,bb,ab,ba,t1jz,ikG)
895 : ALLOCATE(t1jz(vacuum%nmzxy))
896 : !$OMP DO SCHEDULE(dynamic,5)
897 : DO ikG = 1, nv2_outer
898 : IF (.NOT.l_dfpt) THEN
899 : DO ikGPr = 1, ikG - 1
900 : i1 = kvac1(ikG,jspin) - kvac1(ikGPr,jspin)
901 : i2 = kvac2(ikG,jspin) - kvac2(ikGPr,jspin)
902 : i3 = 0
903 : ig3 = stars%ig(i1,i2,i3)
904 : ind2 = stars%ig2(ig3)
905 : !IF (ind2 .ne.stars%ig2(ig3)) CYCLE
906 : IF (iabs(i1).GT.stars%mx1) CYCLE
907 : IF (iabs(i2).GT.stars%mx2) CYCLE
908 : IF (ig3.EQ.0) CYCLE
909 : phs = stars%rgphs(i1,i2,i3)
910 : ig3p = stars%ig(-i1,-i2,i3)
911 : phsp = stars%rgphs(-i1,-i2,i3)
912 : ind2p = stars%ig2(ig3p)
913 : aa = 0.0
914 : bb = 0.0
915 : ba = 0.0
916 : ab = 0.0
917 : DO n = 1,ne
918 : aa=aa+we(n)*CONJG(ac(ikGPr,n,jspin))*ac(ikG,n,jspin)
919 : bb=bb+we(n)*CONJG(bc(ikGPr,n,jspin))*bc(ikG,n,jspin)
920 : ab=ab+we(n)*CONJG(ac(ikGPr,n,jspin))*bc(ikG,n,jspin)
921 : ba=ba+we(n)*CONJG(bc(ikGPr,n,jspin))*ac(ikG,n,jspin)
922 : END DO
923 : DO jz = 1,vacuum%nmzxy
924 : ui = u(jz,ikG,jspin)
925 : uj = u(jz,ikGPr,jspin)
926 : uei = ue(jz,ikG,jspin)
927 : uej = ue(jz,ikGPr,jspin)
928 : t1jz(jz) = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
929 : END DO
930 : !$OMP CRITICAL ! (denvacxy,denvac)
931 : den%vac(:vacuum%nmzxy,ind2,ivac,jspin) = den%vac(:vacuum%nmzxy,ind2, ivac,jspin) &
932 : + t1jz(:vacuum%nmzxy)*phs/stars%nstr2(ind2)
933 : den%vac(:vacuum%nmzxy,ind2p,ivac,jspin) = den%vac(:vacuum%nmzxy,ind2p,ivac,jspin) &
934 : + CONJG(t1jz(:vacuum%nmzxy))*phsp/stars%nstr2(ind2p)
935 : !$OMP END CRITICAL ! (denvacxy,denvac)
936 : END DO
937 : ELSE
938 : DO ikGPr = 1, nv2(jspin)
939 : i1 = kvac1q(ikG,jspin) - kvac1(ikGPr,jspin)
940 : i2 = kvac2q(ikG,jspin) - kvac2(ikGPr,jspin)
941 : i3 = 0
942 : ig3 = stars%ig(i1,i2,i3)
943 : ind2 = stars%ig2(ig3)
944 : IF ((ind2==1).AND.(norm2(stars%center)<1e-8)) CYCLE
945 : IF (iabs(i1).GT.stars%mx1) CYCLE
946 : IF (iabs(i2).GT.stars%mx2) CYCLE
947 : IF (ig3.EQ.0) CYCLE
948 : phs = stars%rgphs(i1,i2,i3)
949 : aa = 0.0
950 : bb = 0.0
951 : ba = 0.0
952 : ab = 0.0
953 : DO n = 1,ne
954 : aa=aa+we(n)*CONJG(ac(ikGPr,n,jspin))*ac1(ikG,n,jspin)
955 : bb=bb+we(n)*CONJG(bc(ikGPr,n,jspin))*bc1(ikG,n,jspin)
956 : ab=ab+we(n)*CONJG(ac(ikGPr,n,jspin))*bc1(ikG,n,jspin)
957 : ba=ba+we(n)*CONJG(bc(ikGPr,n,jspin))*ac1(ikG,n,jspin)
958 : IF (norm2(stars%center)<1e-8) THEN
959 : aa=aa+we1(n)*CONJG(ac(ikGPr,n,jspin))*ac(ikG,n,jspin)
960 : bb=bb+we1(n)*CONJG(bc(ikGPr,n,jspin))*bc(ikG,n,jspin)
961 : ab=ab+we1(n)*CONJG(ac(ikGPr,n,jspin))*bc(ikG,n,jspin)
962 : ba=ba+we1(n)*CONJG(bc(ikGPr,n,jspin))*ac(ikG,n,jspin)
963 : END IF
964 : END DO
965 : DO jz = 1,vacuum%nmzxy
966 : ui = uq(jz,ikG,jspin)
967 : uj = u(jz,ikGPr,jspin)
968 : uei = ueq(jz,ikG,jspin)
969 : uej = ue(jz,ikGPr,jspin)
970 : t1jz(jz) = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
971 : END DO
972 : !$OMP CRITICAL ! (denvacxy,denvac)
973 : den%vac(:vacuum%nmzxy,ind2,ivac,jspin) = den%vac(:vacuum%nmzxy,ind2, ivac,jspin) &
974 : + t1jz(:vacuum%nmzxy)*phs/stars%nstr2(ind2)
975 : !$OMP END CRITICAL ! (denvacxy,denvac)
976 : END DO
977 : END IF
978 : END DO
979 : !$OMP END DO
980 : DEALLOCATE(t1jz)
981 : !$OMP END PARALLEL
982 : END IF
983 : !END IF
984 : !=============================================================
985 : !
986 : ! calculate 1. to nstars. starcoefficient for each k and energy eigenvalue
987 : ! to vacdos%qstars(ne,layer,ivac,ikpt) if starcoeff=T (the star coefficient values are written to vacdos)
988 : !
989 320 : IF (banddos%starcoeff .AND. banddos%vacdos) THEN
990 0 : DO n=1,ne
991 0 : DO ikG = 1,nv2(jspin)
992 0 : DO ikGPr = 1, ikG - 1
993 0 : i1 = kvac1(ikG,jspin) - kvac1(ikGPr,jspin)
994 0 : i2 = kvac2(ikG,jspin) - kvac2(ikGPr,jspin)
995 0 : i3 = 0
996 0 : IF (iabs(i1).GT.stars%mx1) CYCLE
997 0 : IF (iabs(i2).GT.stars%mx2) CYCLE
998 0 : ig3 = stars%ig(i1,i2,i3)
999 0 : IF (ig3.EQ.0) CYCLE
1000 0 : ind2 = stars%ig2(ig3)
1001 0 : ig3p = stars%ig(-i1,-i2,i3)
1002 0 : ind2p = stars%ig2(ig3p)
1003 0 : IF ((ind2.GE.2.AND.ind2.LE.banddos%nstars).OR.&
1004 0 : (ind2p.GE.2.AND.ind2p.LE.banddos%nstars)) THEN
1005 0 : phs = stars%rgphs(i1,i2,i3)
1006 0 : phsp = stars%rgphs(-i1,-i2,i3)
1007 0 : aa = CONJG(ac(ikGPr,n,jspin))*ac(ikG,n,jspin)
1008 0 : bb = CONJG(bc(ikGPr,n,jspin))*bc(ikG,n,jspin)
1009 0 : ab = CONJG(ac(ikGPr,n,jspin))*bc(ikG,n,jspin)
1010 0 : ba = CONJG(bc(ikGPr,n,jspin))*ac(ikG,n,jspin)
1011 0 : DO jj = 1,banddos%layers
1012 0 : ui = u(banddos%izlay(jj,1),ikG,jspin)
1013 0 : uj = u(banddos%izlay(jj,1),ikGPr,jspin)
1014 0 : uei = ue(banddos%izlay(jj,1),ikG,jspin)
1015 0 : uej = ue(banddos%izlay(jj,1),ikGPr,jspin)
1016 0 : t1 = aa*ui*uj + bb*uei*uej +ba*ui*uej + ab*uei*uj
1017 0 : IF (ind2.GE.2.AND.ind2.LE.banddos%nstars) &
1018 0 : vacdos%qstars(ind2-1,ev_list(n),jj,ivac,ikpt,jspin) = vacdos%qstars(ind2-1,ev_list(n),jj,ivac,ikpt,jspin)+ t1*phs/stars%nstr2(ind2)
1019 0 : IF (ind2p.GE.2.AND.ind2p.LE.banddos%nstars) &
1020 0 : vacdos%qstars(ind2p-1,ev_list(n),jj,ivac,ikpt,jspin) = vacdos%qstars(ind2p-1,ev_list(n),jj,ivac,ikpt,jspin) +CONJG(t1)*phs/stars%nstr2(ind2p)
1021 : END DO
1022 : END IF
1023 : END DO
1024 : END DO
1025 : END DO
1026 : END IF
1027 : END DO ! vacuum loop
1028 140 : DEALLOCATE (ac,bc,dt,dte,t,te,tei,u,ue,yy )
1029 :
1030 : ! DEALLOCATE (du,ddu,due,ddue)
1031 : !IF(vacuum%nvac.EQ.1) THEN
1032 : ! den%vacz(:,2,:) = den%vacz(:,1,:)
1033 : ! IF (sym%invs) THEN
1034 : ! den%vacxy(:,:,2,:) = CONJG(den%vacxy(:,:,1,:))
1035 : ! ELSE
1036 : ! den%vacxy(:,:,2,:) = den%vacxy(:,:,1,:)
1037 : ! END IF
1038 : !END IF
1039 :
1040 140 : CALL timestop("vacden")
1041 :
1042 140 : END SUBROUTINE vacden
1043 : END MODULE m_vacden
|
