Line data Source code
1 : MODULE m_olap
2 : USE m_types_hybdat
3 : USE m_types_mat
4 : #ifdef CPP_MPI
5 : use mpi
6 : #endif
7 : private olap_pw_real, olap_pw_cmplx
8 : public olap_pw, olap_pwp, wfolap_inv, wfolap_noinv
9 :
10 : CONTAINS
11 :
12 : ! Calculates plane-wave overlap matrix olap defined by GPT(1:3,1:NGPT).
13 : ! (Muffin-tin spheres are cut out.)
14 : ! olap_pw calculates full overlap matrix
15 :
16 84 : SUBROUTINE olap_pw(olap, gpt, ngpt, atoms, cell, fmpi)
17 : use m_juDFT
18 : USE m_constants
19 : USE m_types_cell
20 : USE m_types_atoms
21 : USE m_types_mpi
22 : USE m_types_mat
23 : IMPLICIT NONE
24 : TYPE(t_cell), INTENT(IN) :: cell
25 : TYPE(t_atoms), INTENT(IN) :: atoms
26 : type(t_mpi), intent(in) :: fmpi
27 :
28 : ! - scalars -
29 : INTEGER, INTENT(IN) :: ngpt
30 : ! - arrays -
31 : INTEGER, INTENT(IN) :: gpt(:, :)
32 : TYPE(t_mat) :: olap
33 :
34 84 : if (olap%l_real) then
35 8 : call olap_pw_real(olap, gpt, ngpt, atoms, cell, fmpi)
36 : else
37 76 : call olap_pw_cmplx(olap, gpt, ngpt, atoms, cell)
38 : endif
39 84 : END SUBROUTINE olap_pw
40 :
41 8 : subroutine olap_pw_real(olap, gpt, ngpt, atoms, cell, fmpi)
42 : use m_juDFT
43 : USE m_constants
44 : USE m_types_cell
45 : USE m_types_atoms
46 : USE m_types_mpi
47 : USE m_types_mat
48 : IMPLICIT NONE
49 : TYPE(t_cell), INTENT(IN) :: cell
50 : TYPE(t_atoms), INTENT(IN) :: atoms
51 : type(t_mpi), intent(in) :: fmpi
52 :
53 : ! - scalars -
54 : INTEGER, INTENT(IN) :: ngpt
55 : ! - arrays -
56 : INTEGER, INTENT(IN) :: gpt(:, :)
57 : TYPE(t_mat) :: olap
58 : ! - local -
59 : INTEGER :: i, j, itype, icent, ierr, root
60 : REAL :: g, r, fgr
61 : COMPLEX, PARAMETER :: img = (0.0, 1.0)
62 : INTEGER :: dg(3)
63 :
64 8 : call timestart("olap_pw_real")
65 : !$OMP PARALLEL default(none) &
66 : !$OMP private(i,j,dg,g,itype, icent, r, fgr)&
67 8 : !$OMP shared(olap, gpt, cell, atoms, ngpt, fmpi)
68 :
69 : !$OMP DO schedule(guided)
70 : DO j = 1+fmpi%n_rank, ngpt, fmpi%n_size
71 : DO i = 1, j
72 : olap%data_r(i,j) = 0.0
73 : dg = gpt(:, j) - gpt(:, i)
74 : g = gptnorm(dg, cell%bmat)
75 : IF (abs(g) < 1e-10) THEN
76 : DO itype = 1, atoms%ntype
77 : r = atoms%rmt(itype)
78 : olap%data_r(i, j) = olap%data_r(i, j) - atoms%neq(itype)*fpi_const*r**3/3/cell%omtil
79 : END DO
80 : ELSE
81 : do icent = 1, atoms%nat
82 : itype = atoms%itype(icent)
83 : r = g*atoms%rmt(itype)
84 : fgr = fpi_const*(sin(r) - r*cos(r))/g**3/cell%omtil
85 : olap%data_r(i, j) = real(olap%data_r(i, j) - fgr*exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent))))
86 : END DO
87 : END IF
88 : END DO
89 : END DO
90 : !$OMP end do
91 :
92 : ! work on diagonal
93 : !$OMP DO
94 : do i = 1+fmpi%n_rank, ngpt, fmpi%n_size
95 : olap%data_r(i,i) = olap%data_r(i,i) + 1
96 : enddo
97 : !$OMP END DO
98 : !$OMP end parallel
99 : #ifdef CPP_MPI
100 952 : do j = 1, ngpt
101 944 : root = mod(j-1, fmpi%n_size)
102 952 : call MPI_Bcast(olap%data_r(1,j), j, MPI_DOUBLE_PRECISION, root, fmpi%sub_comm, ierr)
103 : enddo
104 : #endif
105 8 : call olap%u2l()
106 8 : call timestop("olap_pw_real")
107 8 : END SUBROUTINE olap_pw_real
108 :
109 76 : SUBROUTINE olap_pw_cmplx(olap, gpt, ngpt, atoms, cell)
110 : use m_juDFT
111 : USE m_constants
112 : USE m_types_cell
113 : USE m_types_atoms
114 : IMPLICIT NONE
115 : TYPE(t_cell), INTENT(IN) :: cell
116 : TYPE(t_atoms), INTENT(IN) :: atoms
117 :
118 : ! - scalars -
119 : INTEGER, INTENT(IN) :: ngpt
120 : ! - arrays -
121 : INTEGER, INTENT(IN) :: gpt(:, :)
122 : TYPE(t_mat) :: olap
123 : ! - local -
124 : INTEGER :: i, j, itype, icent
125 : REAL :: g, r, fgr
126 : COMPLEX, PARAMETER :: img = (0.0, 1.0)
127 : INTEGER :: dg(3)
128 :
129 76 : call timestart("olap_pw_cmplx")
130 : !$OMP PARALLEL DO default(none) schedule(guided) &
131 : !$OMP private(i,j,dg,g,itype, icent, r, fgr)&
132 76 : !$OMP shared(olap, gpt, cell, atoms, ngpt)
133 : DO i = 1, ngpt
134 : DO j = 1, i
135 : olap%data_c(i,j) = cmplx_0
136 : dg = gpt(:, j) - gpt(:, i)
137 : g = gptnorm(dg, cell%bmat)
138 : IF (abs(g) < 1e-10) THEN
139 : DO itype = 1, atoms%ntype
140 : r = atoms%rmt(itype)
141 : olap%data_c(i, j) = olap%data_c(i, j) - atoms%neq(itype)*fpi_const*r**3/3/cell%omtil
142 : END DO
143 : ELSE
144 : icent = 0
145 : do icent = 1, atoms%nat
146 : itype = atoms%itype(icent)
147 : r = g*atoms%rmt(itype)
148 : fgr = fpi_const*(sin(r) - r*cos(r))/g**3/cell%omtil
149 : olap%data_c(i, j) = olap%data_c(i, j) - fgr*exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent)))
150 : END DO
151 : END IF
152 : IF (i == j) olap%data_c(i, j) = olap%data_c(i, j) + 1
153 : olap%data_c(j, i) = conjg(olap%data_c(i, j))
154 : END DO
155 : END DO
156 : !$OMP END PARALLEL DO
157 76 : call timestop("olap_pw_cmplx")
158 76 : END SUBROUTINE olap_pw_cmplx
159 :
160 : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161 :
162 : ! olap_pwp calculates upper triangular part of overlap matrix
163 :
164 0 : SUBROUTINE olap_pwp(l_real, olap_r, olap_c, gpt, ngpt, atoms, cell)
165 :
166 : USE m_constants, ONLY: REAL_NOT_INITALIZED, CMPLX_NOT_INITALIZED, &
167 : fpi_const, tpi_const
168 : USE m_types_cell
169 : USE m_types_atoms
170 : IMPLICIT NONE
171 : TYPE(t_cell), INTENT(IN) :: cell
172 : TYPE(t_atoms), INTENT(IN) :: atoms
173 :
174 : ! - scalars -
175 : INTEGER, INTENT(IN) :: ngpt
176 : ! - arrays -
177 : INTEGER, INTENT(IN) :: gpt(:, :)
178 :
179 : LOGICAL, INTENT(IN) :: l_real
180 : REAL, INTENT(INOUT) :: olap_r(ngpt*(ngpt + 1)/2)
181 : COMPLEX, INTENT(INOUT) :: olap_c(ngpt*(ngpt + 1)/2)
182 : ! - local -
183 : INTEGER :: i, j, k, itype, icent, ineq
184 : REAL :: g, r, fgr
185 : COMPLEX, PARAMETER :: img = (0.0, 1.0)
186 : INTEGER :: dg(3)
187 :
188 0 : olap_r = REAL_NOT_INITALIZED; olap_c = CMPLX_NOT_INITALIZED
189 :
190 0 : if (l_real) THEN
191 : k = 0
192 0 : DO i = 1, ngpt
193 0 : DO j = 1, i
194 0 : k = k + 1
195 0 : dg = gpt(:, i) - gpt(:, j)
196 0 : g = gptnorm(dg, cell%bmat)
197 0 : olap_r(k) = 0
198 0 : IF (abs(g) < 1e-10) THEN
199 0 : DO itype = 1, atoms%ntype
200 0 : r = atoms%rmt(itype)
201 0 : olap_r(k) = olap_r(k) - atoms%neq(itype)*fpi_const*r**3/3/cell%omtil
202 : END DO
203 : ELSE
204 0 : icent = 0
205 0 : DO itype = 1, atoms%ntype
206 0 : r = g*atoms%rmt(itype)
207 0 : fgr = fpi_const*(sin(r) - r*cos(r))/g**3/cell%omtil
208 0 : DO ineq = 1, atoms%neq(itype)
209 0 : icent = icent + 1
210 : olap_r(k) = olap_r(k) - real(fgr* &
211 0 : exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent))))
212 : END DO
213 : END DO
214 : END IF
215 0 : IF (i == j) olap_r(k) = olap_r(k) + 1
216 : END DO
217 : END DO
218 : else
219 : k = 0
220 0 : DO i = 1, ngpt
221 0 : DO j = 1, i
222 0 : k = k + 1
223 0 : dg = gpt(:, i) - gpt(:, j)
224 0 : g = gptnorm(dg, cell%bmat)
225 0 : olap_c(k) = 0
226 0 : IF (abs(g) < 1e-10) THEN
227 0 : DO itype = 1, atoms%ntype
228 0 : r = atoms%rmt(itype)
229 0 : olap_c(k) = olap_c(k) - atoms%neq(itype)*fpi_const*r**3/3/cell%omtil
230 : END DO
231 : ELSE
232 0 : icent = 0
233 0 : DO itype = 1, atoms%ntype
234 0 : r = g*atoms%rmt(itype)
235 0 : fgr = fpi_const*(sin(r) - r*cos(r))/g**3/cell%omtil
236 0 : DO ineq = 1, atoms%neq(itype)
237 0 : icent = icent + 1
238 : olap_c(k) = olap_c(k) - fgr* &
239 0 : exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent)))
240 : END DO
241 : END DO
242 : END IF
243 0 : IF (i == j) olap_c(k) = olap_c(k) + 1
244 : END DO
245 : END DO
246 :
247 : endif
248 0 : END SUBROUTINE olap_pwp
249 :
250 : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
251 :
252 : ! SUBROUTINE wfolap_init(olappw, olapmt, gpt, &
253 : ! atoms, mpdata, cell, &
254 : ! bas1, bas2)
255 :
256 : ! USE m_intgrf, ONLY: intgrf, intgrf_init
257 : ! USE m_types
258 : ! IMPLICIT NONE
259 : ! TYPE(t_mpdata), intent(in) :: mpdata
260 : ! TYPE(t_cell), INTENT(IN) :: cell
261 : ! TYPE(t_atoms), INTENT(IN) :: atoms
262 :
263 : ! ! - arrays -
264 : ! INTEGER, INTENT(IN) :: gpt(:, :)!(3,ngpt)
265 : ! REAL, INTENT(IN) :: bas1(atoms%jmtd, maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype), &
266 : ! bas2(atoms%jmtd, maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
267 : ! REAL, INTENT(INOUT) :: olapmt(maxval(mpdata%num_radfun_per_l), &
268 : ! maxval(mpdata%num_radfun_per_l), &
269 : ! 0:atoms%lmaxd, &
270 : ! atoms%ntype)
271 : ! TYPE(t_mat), INTENT(INOUT):: olappw
272 :
273 : ! ! - local -
274 : ! INTEGER :: itype, l, nn, n1, n2
275 :
276 : ! REAL, ALLOCATABLE :: gridf(:, :)
277 :
278 : ! CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, gridf)
279 : ! olapmt = 0
280 : ! DO itype = 1, atoms%ntype
281 : ! DO l = 0, atoms%lmax(itype)
282 : ! nn = mpdata%num_radfun_per_l(l, itype)
283 : ! DO n2 = 1, nn
284 : ! DO n1 = 1, nn!n2
285 : ! !IF( n1 .gt. 2 .or. n2 .gt. 2) CYCLE
286 : ! olapmt(n1, n2, l, itype) = intgrf( &
287 : ! bas1(:, n1, l, itype)*bas1(:, n2, l, itype) &
288 : ! + bas2(:, n1, l, itype)*bas2(:, n2, l, itype), &
289 : ! atoms, itype, gridf)
290 : ! ! olapmt(n2,n1,l,itype) = olapmt(n1,n2,l,itype)
291 : ! END DO
292 : ! END DO
293 : ! END DO
294 : ! END DO
295 :
296 : ! CALL olap_pw(olappw, gpt, size(gpt, 2), atoms, cell)
297 :
298 : ! END SUBROUTINE wfolap_init
299 :
300 0 : FUNCTION wfolap_inv(cmt1, cpw1, cmt2, cpw2, olappw, olapmt, atoms, mpdata)
301 :
302 : USE m_wrapper
303 : USE m_types_mpdata
304 : USE m_types_atoms
305 : IMPLICIT NONE
306 : TYPE(t_mpdata), intent(in) :: mpdata
307 : TYPE(t_atoms), INTENT(IN) :: atoms
308 :
309 : ! - scalars -
310 : COMPLEX :: wfolap_inv
311 : ! - arrays -
312 : COMPLEX, INTENT(IN) :: cmt1(:, :), cmt2(:, :)
313 : REAL, INTENT(IN) :: cpw1(:)
314 : COMPLEX, INTENT(IN) :: cpw2(:)
315 : REAL, INTENT(IN) :: olappw(:, :)
316 : REAL, INTENT(IN) :: olapmt(maxval(mpdata%num_radfun_per_l), maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
317 : ! - local -
318 : INTEGER :: itype, ieq, iatom, l, m, lm, nn
319 :
320 0 : wfolap_inv = 0
321 0 : iatom = 0
322 0 : DO itype = 1, atoms%ntype
323 0 : DO ieq = 1, atoms%neq(itype)
324 0 : iatom = iatom + 1
325 0 : lm = 0
326 0 : DO l = 0, atoms%lmax(itype)
327 0 : DO M = -l, l
328 0 : nn = mpdata%num_radfun_per_l(l, itype)
329 : wfolap_inv = wfolap_inv + &
330 : dot_product(cmt1(lm + 1:lm + nn, iatom), &
331 0 : matmul(olapmt(:nn, :nn, l, itype), &
332 0 : cmt2(lm + 1:lm + nn, iatom)))
333 0 : lm = lm + nn
334 : END DO
335 : END DO
336 : END DO
337 : END DO
338 :
339 0 : wfolap_inv = wfolap_inv + dot_product(cpw1, matmul(olappw, cpw2))
340 :
341 : ! CALL dgemv('N',ngpt1,ngpt2,1.0,olappw,ngpt1,real(cpw2),1,0.0,rarr1,1)
342 : ! CALL dgemv('N',ngpt1,ngpt2,1.0,olappw,ngpt1,aimag(cpw2),1,0.0,rarr2,1)
343 : !
344 : ! rdum1 = dot_product(cpw1,rarr1)
345 : ! rdum2 = dot_product(cpw1,rarr2)
346 : ! cdum = cmplx( rdum1, rdum2 )
347 :
348 : ! wfolap = wfolap + cdum
349 :
350 0 : END FUNCTION wfolap_inv
351 :
352 0 : FUNCTION wfolap_noinv(cmt1, cpw1, cmt2, cpw2, olappw, olapmt, atoms, mpdata)
353 :
354 : USE m_wrapper
355 : USE m_types_mpdata
356 : USE m_types_atoms
357 :
358 : IMPLICIT NONE
359 : TYPE(t_mpdata), intent(in) :: mpdata
360 : TYPE(t_atoms), INTENT(IN) :: atoms
361 :
362 : ! - scalars -
363 : COMPLEX :: wfolap_noinv
364 : ! - arrays -
365 : COMPLEX, INTENT(IN) :: cmt1(:, :), cmt2(:, :)
366 : COMPLEX, INTENT(IN) :: cpw1(:)
367 : COMPLEX, INTENT(IN) :: cpw2(:)
368 : COMPLEX, INTENT(IN) :: olappw(:, :)
369 : REAL, INTENT(IN) :: olapmt(maxval(mpdata%num_radfun_per_l), maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
370 : ! - local -
371 : INTEGER :: itype, ieq, iatom, l, m, lm, nn
372 :
373 0 : wfolap_noinv = 0
374 0 : iatom = 0
375 0 : DO itype = 1, atoms%ntype
376 0 : DO ieq = 1, atoms%neq(itype)
377 0 : iatom = iatom + 1
378 0 : lm = 0
379 0 : DO l = 0, atoms%lmax(itype)
380 0 : DO M = -l, l
381 0 : nn = mpdata%num_radfun_per_l(l, itype)
382 : wfolap_noinv = wfolap_noinv + &
383 : dot_product(cmt1(lm + 1:lm + nn, iatom), &
384 0 : matmul(olapmt(:nn, :nn, l, itype), &
385 0 : cmt2(lm + 1:lm + nn, iatom)))
386 0 : lm = lm + nn
387 : END DO
388 : END DO
389 : END DO
390 : END DO
391 :
392 0 : wfolap_noinv = wfolap_noinv + dot_product(cpw1, matmul(olappw, cpw2))
393 :
394 : ! CALL dgemv('N',ngpt1,ngpt2,1.0,olappw,ngpt1,real(cpw2),1,0.0,rarr1,1)
395 : ! CALL dgemv('N',ngpt1,ngpt2,1.0,olappw,ngpt1,aimag(cpw2),1,0.0,rarr2,1)
396 : !
397 : ! rdum1 = dot_product(cpw1,rarr1)
398 : ! rdum2 = dot_product(cpw1,rarr2)
399 : ! cdum = cmplx( rdum1, rdum2 )
400 :
401 : ! wfolap = wfolap + cdum
402 :
403 0 : END FUNCTION wfolap_noinv
404 0 : END MODULE m_olap
|
