Line data Source code
1 : MODULE m_checkolap
2 : use m_ylm
3 : CONTAINS
4 :
5 0 : SUBROUTINE checkolap(atoms, hybdat, mpdata, hybinp, nkpti, kpts, fmpi, &
6 : input, sym, noco, nococonv, cell, lapw, jsp)
7 : USE m_util, ONLY: chr, sphbessel, harmonicsr
8 : use m_intgrf, only: intgrf, intgrf_init
9 : use m_calc_cmt
10 : USE m_constants
11 : USE m_types
12 : USE m_io_hybrid
13 : USE m_types_hybdat
14 : use m_calc_l_m_from_lm
15 : #ifdef CPP_MPI
16 : use mpi
17 : #endif
18 :
19 : IMPLICIT NONE
20 :
21 : TYPE(t_hybdat), INTENT(IN) :: hybdat
22 :
23 : TYPE(t_mpi), INTENT(IN) :: fmpi
24 : TYPE(t_mpdata), intent(in) :: mpdata
25 : TYPE(t_hybinp), INTENT(IN) :: hybinp
26 : TYPE(t_input), INTENT(IN) :: input
27 : TYPE(t_noco), INTENT(IN) :: noco
28 : TYPE(t_nococonv), INTENT(IN) :: nococonv
29 : TYPE(t_sym), INTENT(IN) :: sym
30 : TYPE(t_cell), INTENT(IN) :: cell
31 : TYPE(t_kpts), INTENT(IN) :: kpts
32 : TYPE(t_atoms), INTENT(IN) :: atoms
33 :
34 : TYPE(t_lapw), INTENT(INOUT) :: lapw
35 :
36 : ! - scalars -
37 : INTEGER, INTENT(IN) :: jsp
38 : INTEGER, INTENT(IN) :: nkpti
39 :
40 : ! - local scalars -
41 : INTEGER :: i, itype, iatom, ikpt, ineq, igpt, iband
42 : INTEGER :: j, m
43 : INTEGER :: l
44 : INTEGER :: lm, lm1
45 : INTEGER :: n, nbasfcn
46 :
47 : REAL :: rdum, rdum1
48 : REAL :: qnorm
49 :
50 : COMPLEX :: cexp, cdum, pre_fac
51 : COMPLEX, PARAMETER :: img = (0.0, 1.0)
52 :
53 : ! -local arrays -
54 : INTEGER :: iarr(2), gpt(3), ierr
55 0 : INTEGER, ALLOCATABLE :: olapcv_loc(:, :, :, :, :)
56 :
57 0 : REAL :: sphbes(0:atoms%lmaxd)
58 : REAL :: q(3)
59 0 : REAL :: integrand(atoms%jmtd)
60 0 : REAL :: rarr(maxval(hybdat%nbands))
61 0 : REAL, ALLOCATABLE :: olapcb(:)
62 0 : REAL, ALLOCATABLE :: olapcv_avg(:, :, :, :), olapcv_max(:, :, :, :)
63 0 : TYPE(t_mat), ALLOCATABLE :: z(:)
64 :
65 0 : COMPLEX :: cmt(input%neig, hybdat%maxlmindx, atoms%nat, nkpti)
66 0 : COMPLEX :: y((atoms%lmaxd + 1)**2)
67 0 : COMPLEX, ALLOCATABLE :: olapcv(:, :), c_phase(:)
68 0 : COMPLEX, ALLOCATABLE :: carr1(:, :), carr2(:, :), carr3(:, :)
69 :
70 : CHARACTER, PARAMETER :: lchar(0:38) = &
71 : (/'s', 'p', 'd', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', &
72 : 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', &
73 : 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x'/)
74 : LOGICAL, parameter :: l_mism = .false.
75 :
76 0 : call timestart("checkolap")
77 0 : allocate(z(nkpti))
78 0 : DO ikpt = 1, nkpti
79 0 : CALL lapw%init(input, noco, nococonv, kpts, atoms, sym, ikpt, cell)
80 0 : nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco)
81 0 : call z(ikpt)%alloc(sym%invs, nbasfcn, input%neig)
82 : ENDDO
83 :
84 0 : IF(fmpi%irank == 0) WRITE(oUnit, '(//A)') '### checkolap ###'
85 :
86 0 : cmt = 0
87 :
88 : ! initialize gridf -> was done in eigen_HF_init
89 : !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
90 :
91 : ! read in cmt
92 0 : DO ikpt = 1, nkpti
93 0 : if(allocated(c_phase)) deallocate(c_phase)
94 0 : allocate(c_phase(hybdat%nbands(ikpt, jsp)))
95 :
96 0 : if(ikpt /= kpts%bkp(ikpt)) call juDFT_error("We should be reading the parent z-mat here!")
97 : call read_z(atoms, cell, hybdat, kpts, sym, noco, nococonv, input, ikpt, &
98 0 : jsp, z(ikpt), c_phase=c_phase)
99 : #ifdef CPP_MPI
100 : ! call timestart("Post read_z Barrier: checkolap")
101 : ! call MPI_Barrier(MPI_COMM_WORLD, ierr)
102 : ! call timestop("Post read_z Barrier: checkolap")
103 : #endif
104 : call calc_cmt(atoms, cell, input, noco, nococonv, hybinp, hybdat, mpdata, kpts, &
105 : sym, z(kpts%bkp(ikpt)), jsp, ikpt, c_phase, &
106 0 : cmt(:hybdat%nbands(ikpt,jsp), :, :, ikpt))
107 : END DO
108 :
109 0 : IF(fmpi%irank == 0) WRITE(oUnit, '(/A)') ' Overlap <core|core>'
110 0 : DO itype = 1, atoms%ntype
111 0 : IF(atoms%ntype > 1 .AND. fmpi%irank == 0) &
112 0 : WRITE(oUnit, '(A,I3)') ' Atom type', itype
113 0 : DO l = 0, hybdat%lmaxc(itype)
114 0 : DO i = 1, hybdat%nindxc(l, itype)
115 0 : IF(fmpi%irank == 0) &
116 0 : WRITE(oUnit, '(1x,I1,A,2X)', advance='no') i + l, lchar(l)
117 0 : DO j = 1, i
118 : integrand = hybdat%core1(:, i, l, itype)*hybdat%core1(:, j, l, itype) &
119 0 : + hybdat%core2(:, i, l, itype)*hybdat%core2(:, j, l, itype)
120 0 : IF(fmpi%irank == 0) WRITE(oUnit, '(F10.6)', advance='no') &
121 0 : intgrf(integrand, atoms, itype, hybdat%gridf)
122 : END DO
123 0 : IF(fmpi%irank == 0) WRITE(oUnit, *)
124 : END DO
125 : END DO
126 : END DO
127 :
128 0 : IF(fmpi%irank == 0) WRITE(oUnit, '(/A)') ' Overlap <core|basis>'
129 : allocate(olapcb(maxval(mpdata%num_radfun_per_l)), olapcv(maxval(hybdat%nbands), nkpti), &
130 : olapcv_avg(-hybdat%lmaxcd:hybdat%lmaxcd, hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype), &
131 : olapcv_max(-hybdat%lmaxcd:hybdat%lmaxcd, hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype), &
132 0 : olapcv_loc(2, -hybdat%lmaxcd:hybdat%lmaxcd, hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype))
133 :
134 0 : DO itype = 1, atoms%ntype
135 0 : IF(atoms%ntype > 1 .AND. fmpi%irank == 0) &
136 0 : WRITE(oUnit, '(A,I3)') ' Atom type', itype
137 0 : DO l = 0, hybdat%lmaxc(itype)
138 0 : IF(l > atoms%lmax(itype)) EXIT ! very improbable case
139 0 : IF(fmpi%irank == 0) &
140 : WRITE(oUnit, "(9X,'u(',A,')',4X,'udot(',A,')',:,3X,'ulo(',A,"// &
141 0 : "') ...')")(lchar(l), i=1, min(3, mpdata%num_radfun_per_l(l, itype)))
142 0 : DO i = 1, hybdat%nindxc(l, itype)
143 0 : IF(fmpi%irank == 0) &
144 0 : WRITE(oUnit, '(1x,I1,A,2X)', advance='no') i + l, lchar(l)
145 0 : DO j = 1, mpdata%num_radfun_per_l(l, itype)
146 :
147 : integrand = hybdat%core1(:, i, l, itype)*hybdat%bas1(:, j, l, itype) &
148 0 : + hybdat%core2(:, i, l, itype)*hybdat%bas2(:, j, l, itype)
149 :
150 : olapcb(j) = &
151 0 : intgrf(integrand, atoms, itype, hybdat%gridf)
152 :
153 0 : IF(fmpi%irank == 0) &
154 0 : WRITE(oUnit, '(F10.6)', advance='no') olapcb(j)
155 : END DO
156 :
157 0 : lm = sum([(mpdata%num_radfun_per_l(j, itype)*(2*j + 1), j=0, l - 1)])
158 0 : iatom = atoms%firstAtom(itype)
159 0 : DO m = -l, l
160 0 : olapcv = 0
161 0 : DO j = 1, mpdata%num_radfun_per_l(l, itype)
162 0 : lm = lm + 1
163 : olapcv(:, :) = olapcv(:, :) + &
164 0 : olapcb(j)*cmt(:maxval(hybdat%nbands), lm, iatom, :nkpti)
165 : END DO
166 0 : rdum = sum(abs(olapcv(:, :))**2)
167 0 : rdum1 = maxval(abs(olapcv(:, :)))
168 0 : iarr = maxloc(abs(olapcv(:, :)))
169 : olapcv_avg(m, i, l, itype) = &
170 0 : sqrt(rdum/nkpti/sum(hybdat%nbands(:nkpti,jsp))*nkpti)
171 0 : olapcv_max(m, i, l, itype) = rdum1
172 0 : olapcv_loc(:, m, i, l, itype) = iarr
173 : END DO
174 0 : IF(fmpi%irank == 0) WRITE(oUnit, *)
175 :
176 : END DO
177 : END DO
178 : END DO
179 :
180 0 : IF(fmpi%irank == 0) THEN
181 0 : WRITE(oUnit, '(/A)') ' Average overlap <core|val>'
182 0 : DO itype = 1, atoms%ntype
183 0 : IF(atoms%ntype > 1) write(oUnit, '(A,I3)') ' Atom type', itype
184 0 : DO l = 0, hybdat%lmaxc(itype)
185 0 : DO i = 1, hybdat%nindxc(l, itype)
186 0 : WRITE(oUnit, '(1x,I1,A,2X)', advance='no') i + l, lchar(l)
187 : WRITE(oUnit, '('//chr(2*l + 1)//'F10.6)') &
188 0 : olapcv_avg(-l:l, i, l, itype)
189 : END DO
190 : END DO
191 : END DO
192 :
193 0 : WRITE(oUnit, '(/A)') ' Maximum overlap <core|val> at (band/kpoint)'
194 0 : DO itype = 1, atoms%ntype
195 0 : IF(atoms%ntype > 1) write(oUnit, '(A,I3)') ' Atom type', itype
196 0 : DO l = 0, hybdat%lmaxc(itype)
197 0 : DO i = 1, hybdat%nindxc(l, itype)
198 0 : WRITE(oUnit, '(1x,I1,A,2X)', advance='no') i + l, lchar(l)
199 : WRITE(oUnit, '('//chr(2*l + 1)// &
200 : '(F10.6,'' ('',I3.3,''/'',I4.3,'')''))') &
201 0 : (olapcv_max(m, i, l, itype), &
202 0 : olapcv_loc(:, m, i, l, itype), m=-l, l)
203 : END DO
204 : END DO
205 : END DO
206 : END IF ! fmpi%irank == 0
207 :
208 0 : deallocate(olapcb, olapcv, olapcv_avg, olapcv_max, olapcv_loc)
209 :
210 0 : IF(fmpi%irank == 0) WRITE(oUnit, '(/A)') ' Overlap <basis|basis>'
211 :
212 0 : DO itype = 1, atoms%ntype
213 0 : IF(atoms%ntype > 1 .AND. fmpi%irank == 0) &
214 0 : WRITE(oUnit, '(A,I3)') ' Atom type', itype
215 0 : DO l = 0, atoms%lmax(itype)
216 0 : DO i = 1, mpdata%num_radfun_per_l(l, itype)
217 0 : IF(fmpi%irank == 0) THEN
218 0 : SELECT CASE(i)
219 : CASE(1)
220 0 : WRITE(oUnit, '(1x,'' u('',A,'')'')', advance='no') lchar(l)
221 : CASE(2)
222 0 : WRITE(oUnit, '(1x,''udot('',A,'')'')', advance='no') lchar(l)
223 : CASE DEFAULT
224 0 : WRITE(oUnit, '(1x,'' ulo('',A,'')'')', advance='no') lchar(l)
225 : END SELECT
226 : END IF
227 0 : DO j = 1, i
228 : integrand = hybdat%bas1(:, i, l, itype)*hybdat%bas1(:, j, l, itype) &
229 0 : + hybdat%bas2(:, i, l, itype)*hybdat%bas2(:, j, l, itype)
230 :
231 0 : IF(fmpi%irank == 0) WRITE(oUnit, '(F10.6)', advance='no') &
232 0 : intgrf(integrand, atoms, itype, hybdat%gridf)
233 : END DO
234 0 : IF(fmpi%irank == 0) WRITE(oUnit, *)
235 : END DO
236 : END DO
237 : END DO
238 :
239 : IF(l_mism)then
240 :
241 : IF(fmpi%irank == 0) WRITE(oUnit, '(/A)') &
242 : 'Mismatch of wave functions at the MT-sphere boundaries'
243 : allocate(carr1(maxval(hybdat%nbands),(atoms%lmaxd + 1)**2))
244 : allocate(carr2(maxval(hybdat%nbands),(atoms%lmaxd + 1)**2))
245 : allocate(carr3(maxval(hybdat%nbands),(atoms%lmaxd + 1)**2))
246 :
247 : ! create lock for race-condition in coulomb
248 : DO ikpt = 1, nkpti
249 : iatom = 0
250 : DO itype = 1, atoms%ntype
251 : DO ineq = 1, atoms%neq(itype)
252 : iatom = iatom + 1
253 : carr1 = 0; carr2 = 0; carr3 = 0
254 :
255 : ! calculate k1,k2,k3
256 : CALL lapw%init(input, noco, nococonv, kpts, atoms, sym, ikpt, cell)
257 : call timestart("pw part")
258 : ! PW part
259 : !$OMP PARALLEL DO default(none) &
260 : !$OMP private(igpt, gpt, cexp, q, qnorm, sphbes, y, pre_fac, lm, l, m, iband, cdum) &
261 : !$OMP shared(lapw, jsp, atoms, kpts, iatom, ikpt, cell, itype, ineq, z, hybdat) &
262 : !$OMP reduction(+:carr2)
263 : DO igpt = 1, lapw%nv(jsp)
264 : gpt = lapw%gvec(:, igpt, jsp)
265 :
266 : cexp = exp(img*2*pi_const* &
267 : dot_product(kpts%bkf(:, ikpt) + gpt, atoms%taual(:, iatom)))
268 : q = matmul(kpts%bkf(:, ikpt) + gpt, cell%bmat)
269 :
270 : qnorm = norm2(q)
271 : call sphbessel(sphbes, atoms%rmt(itype)*qnorm, atoms%lmax(itype))
272 :
273 : call ylm4(atoms%lmax(itype), q, y)
274 : y = conjg(y)
275 :
276 : pre_fac = fpi_const / sqrt(cell%omtil) * cexp
277 : if(z(1)%l_real) THEN
278 : do lm = 1, (atoms%lmax(itype)+1)**2
279 : call calc_l_m_from_lm(lm, l, m)
280 : DO iband = 1, hybdat%nbands(ikpt,jsp)
281 : cdum = pre_fac * ImagUnit**l * sphbes(l)
282 : carr2(iband, lm) = carr2(iband, lm) + cdum*z(ikpt)%data_r(igpt, iband)*y(lm)
283 : enddo
284 : enddo
285 : else
286 : do lm = 1, (atoms%lmax(itype)+1)**2
287 : call calc_l_m_from_lm(lm, l, m)
288 : DO iband = 1, hybdat%nbands(ikpt,jsp)
289 : cdum = pre_fac * ImagUnit**l * sphbes(l)
290 : carr2(iband, lm) = carr2(iband, lm) + cdum*z(ikpt)%data_c(igpt, iband)*y(lm)
291 : end DO
292 : END DO
293 : endif
294 : enddo
295 : !$OMP END PARALLEL DO
296 : call timestop("pw part")
297 :
298 : call timestart("MT-part")
299 : ! MT
300 : lm1 = 0
301 : do lm = 1,(atoms%lmax(itype)+1)**2
302 : call calc_l_m_from_lm(lm, l, m)
303 : DO n = 1, mpdata%num_radfun_per_l(l, itype)
304 : lm1 = lm1 + 1
305 : rdum = hybdat%bas1(atoms%jri(itype), n, l, itype)/atoms%rmt(itype)
306 : DO iband = 1, hybdat%nbands(ikpt,jsp)
307 : carr3(iband, lm) = carr3(iband, lm) + cmt(iband, lm1, iatom, ikpt)*rdum
308 : END DO
309 : END DO
310 : END DO
311 : call timestop("MT-part")
312 : carr1 = carr2 - carr3
313 :
314 :
315 : rarr = 0
316 : lm = 0
317 : DO l = 0, atoms%lmax(itype)
318 : DO m = -l, l
319 : lm = lm + 1
320 : rarr = rarr + abs(carr1(:, lm))**2
321 : END DO
322 : END DO
323 : rarr = sqrt(rarr/(4*pi_const))
324 :
325 : write (oUnit, '(I6,4X,F14.12,'' ('',F14.12,'')'')') ikpt,sum(rarr(:1)**2/hybdat%nbands(ikpt,jsp)),maxval(rarr(:1))
326 : END DO
327 : END DO
328 : END DO
329 : endif
330 0 : call timestop("checkolap")
331 0 : END SUBROUTINE checkolap
332 :
333 : END MODULE m_checkolap
|
