Line data Source code
1 : module m_structureconstant
2 : USE m_types
3 : USE m_juDFT
4 : USE m_constants
5 : use m_ylm
6 : use m_sort
7 : #ifdef CPP_MPI
8 : use mpi
9 : #endif
10 : contains
11 : ! -----------------------------------------------------------------------------------------------
12 :
13 : ! Calculates the structure constant
14 : ! 1 * ^
15 : ! structconst(lm,ic1,ic2,k) = SUM exp(ikT) ----------------------- Y ( T + R(ic) )
16 : ! T | T + R(ic1) - R(ic2) | lm
17 : !
18 : ! with T = lattice vectors
19 : !
20 : ! An Ewald summation method devised by O.K. Andersen is used for l<5
21 : ! (see e.g. H.L. Skriver, "The LMTO method", Springer 1984).
22 : ! (The real-space function G can be calculated with gfunction.f)
23 : !
24 :
25 12 : SUBROUTINE structureconstant(structconst, cell, hybinp, atoms, kpts, fmpi)
26 : IMPLICIT NONE
27 : TYPE(t_mpi), INTENT(IN) :: fmpi
28 : TYPE(t_hybinp), INTENT(IN) :: hybinp
29 : TYPE(t_cell), INTENT(IN) :: cell
30 : TYPE(t_atoms), INTENT(IN) :: atoms
31 : TYPE(t_kpts), INTENT(IN) :: kpts
32 : ! - scalars -
33 :
34 : ! - arrays -
35 : COMPLEX, INTENT(INOUT) :: structconst(:, :, :, :)
36 :
37 : ! - local scalars -
38 : INTEGER :: i, ic1, ic2, lm, ikpt, l, ishell, nshell
39 : INTEGER :: m
40 : INTEGER :: nptsh, maxl
41 :
42 : REAL :: rad, rrad, rdum
43 : REAL :: a, a1, aa
44 : REAL :: pref, rexp
45 : REAL :: scale
46 :
47 : COMPLEX :: cexp
48 :
49 : LOGICAL, SAVE :: first = .TRUE.
50 : logical :: run_loop
51 : ! - local arrays -
52 12 : INTEGER :: conv(0:2*hybinp%lexp), ierr, buf_sz, root
53 12 : INTEGER, ALLOCATABLE :: ptsh(:, :)
54 :
55 : REAL :: k(3), ki(3), ka(3)
56 12 : REAL :: convpar(0:2*hybinp%lexp), g(0:2*hybinp%lexp)
57 12 : REAL, ALLOCATABLE :: radsh(:)
58 :
59 12 : COMPLEX :: y((2*hybinp%lexp + 1)**2)
60 : REAL, PARAMETER :: CONVPARAM = 1e-18
61 : ! Do some additional shells ( real-space and Fourier-space sum )
62 : INTEGER, PARAMETER :: ADDSHELL2 = 0
63 :
64 12 : call timestart("calc struc_const.")
65 :
66 12 : IF (fmpi%irank /= 0) first = .FALSE.
67 :
68 12 : rdum = cell%vol**(1.0/3) ! define "average lattice parameter"
69 :
70 : ! ewaldlambda = ewaldscale
71 12 : scale = hybinp%ewaldlambda/rdum
72 :
73 : ! lambda = ewaldlambda / rdum
74 :
75 12 : pref = fpi_const/(scale**3*cell%vol)
76 :
77 408 : DO l = 0, 2*hybinp%lexp
78 408 : convpar(l) = CONVPARAM/scale**(l + 1)
79 : END DO
80 :
81 12 : IF (first) THEN
82 3 : WRITE (oUnit, '(//A)') '### subroutine: structureconstant ###'
83 3 : WRITE (oUnit, '(/A)') 'Real-space sum:'
84 : END IF
85 :
86 : !
87 : ! Determine cutoff radii for real-space and Fourier-space summation
88 : ! (1) real space
89 12 : call timestart("determine cutoff radii")
90 :
91 12 : a = 0
92 12 : run_loop = .True.
93 820 : do while(run_loop)
94 808 : a = a + 1
95 808 : rexp = EXP(-a)
96 808 : g(0) = rexp/a*(1 + a*11/16*(1 + a*3/11*(1 + a/9)))
97 808 : g(1) = rexp/a**2*(1 + a*(1 + a/2*(1 + a*7/24*(1 + a/7))))
98 : g(2) = rexp/a**3*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a*3/16 &
99 808 : *(1 + a/9))))))
100 : g(3) = rexp/a**4*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6 &
101 808 : *(1 + a/8)))))))
102 : g(4) = rexp/a**5*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6 &
103 808 : *(1 + a/7*(1 + a/8*(1 + a/10)))))))))
104 : g(5) = rexp/a**6*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6 &
105 808 : *(1 + a/7*(1 + a/8*(1 + a/9*(1 + a/10))))))))))
106 : g(6) = rexp/a**7*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6 &
107 808 : *(1 + a/7*(1 + a/8*(1 + a/9*(1 + a/10*(1 + a/11*(1 + a/12))))))))))))
108 : g(7) = rexp/a**8*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6 &
109 808 : *(1 + a/7*(1 + a/8*(1 + a/9*(1 + a/10*(1 + a/11*(1 + a/12*(1 + a/13)))))))))))))
110 21008 : DO l = 8, 2*hybinp%lexp
111 21008 : g(l) = a**(-l - 1)
112 : END DO
113 3136 : run_loop = ANY(g > convpar/10)
114 : enddo
115 12 : rad = a/scale
116 12 : call timestop("determine cutoff radii")
117 :
118 : ! (2) Fourier space
119 12 : call timestart("fourier space")
120 12 : a = 0
121 12 : run_loop = .True.
122 964 : do while(run_loop)
123 952 : a = a + 1
124 952 : aa = (1 + a**2)**(-1)
125 952 : g(0) = pref*aa**4/a**2
126 952 : g(1) = pref*aa**4/a
127 952 : g(2) = pref*aa**5/3
128 952 : g(3) = pref*aa**5*a/15
129 952 : g(4) = pref*aa**6*a**2/105
130 952 : g(5) = pref*aa**6*a**3/945
131 952 : g(6) = pref*aa**7*a**4/10395
132 952 : g(7) = pref*aa**7*a**5/135135
133 1652 : run_loop = ANY(g > convpar)
134 : enddo
135 : ! IF (ANY(g > convpar)) THEN
136 : ! a = a + 1
137 : ! GOTO 2
138 : ! END IF
139 12 : rrad = a*scale
140 12 : call timestop("fourier space")
141 :
142 12 : IF (first) THEN
143 3 : WRITE (oUnit, '(/A,2F10.5)') 'Cutoff radii: ', rad, rrad
144 3 : WRITE (oUnit, '(/A)') 'Real-space sum'
145 : END IF
146 :
147 12 : call realspace_sum(atoms, cell, hybinp, fmpi, kpts, first, scale, convpar, g, a, a1, rad, structconst)
148 :
149 12 : IF (first) WRITE (oUnit, '(/A)') 'Fourier-space sum'
150 :
151 : !
152 : ! Determine reciprocal shells
153 : !
154 12 : call timestart("determince reciproc. shell")
155 12 : CALL getshells(ptsh, nptsh, radsh, nshell, rrad, cell%bmat, first)
156 12 : call timestop("determince reciproc. shell")
157 : ! minimum nonzero reciprocal-shell radius (needed in routines concerning the non-local hartree-fock exchange)
158 : !hybinp%radshmin = radsh(2)
159 : !
160 : ! Fourier-space sum
161 : !
162 12 : call timestart("fourierspace sum")
163 48 : DO ikpt = 1, kpts%nkpt
164 144 : k = kpts%bk(:, ikpt)
165 36 : maxl = MIN(7, hybinp%lexp*2)
166 36 : ishell = 1
167 1224 : conv = HUGE(i)
168 8286936 : DO i = 1, nptsh
169 8286900 : IF (i > 1) THEN
170 8286864 : IF (ABS(radsh(i) - radsh(i - 1)) > 1e-10) ishell = ishell + 1
171 : ENDIF
172 33147600 : ki = ptsh(:, i) + k - NINT(k) ! -nint(...) transforms to Wigner-Seitz cell ( i.e. -0.5 <= x,y,z < 0.5 )
173 107729700 : ka = MATMUL(ki, cell%bmat)
174 33147600 : a = norm2(ka)/scale
175 8286900 : aa = (1 + a**2)**(-1)
176 8286900 : IF (ABS(a - a1) > 1e-10) THEN
177 4230044 : a1 = a
178 4230044 : IF (abs(a) < 1e-12) THEN
179 12 : g(0) = pref*(-4)
180 12 : g(1) = 0
181 : ELSE
182 4230032 : IF (ishell <= conv(0)) g(0) = pref*aa**4/a**2
183 4230032 : IF (ishell <= conv(1)) g(1) = pref*aa**4/a
184 : END IF
185 4230044 : IF (ishell <= conv(2)) g(2) = pref*aa**5/3
186 4230044 : IF (ishell <= conv(3)) g(3) = pref*aa**5*a/15
187 4230044 : IF (ishell <= conv(4)) g(4) = pref*aa**6*a**2/105
188 4230044 : IF (ishell <= conv(5)) g(5) = pref*aa**6*a**3/945
189 4230044 : IF (ishell <= conv(6)) g(6) = pref*aa**7*a**4/10395
190 4230044 : IF (ishell <= conv(7)) g(7) = pref*aa**7*a**5/135135
191 4230044 : IF (ishell > 1) THEN
192 38070072 : DO l = 0, 7
193 38070072 : IF (conv(l) == HUGE(i) .AND. g(l) < convpar(l)) conv(l) = ishell + ADDSHELL2
194 : END DO
195 : END IF
196 : END IF
197 :
198 8286900 : IF (ishell > conv(maxl) .AND. maxl /= 0) maxl = maxl - 1
199 8286900 : call ylm4(maxl, ka, y)
200 33148356 : IF (norm2(ka(:)) .LT. 1.0e-16) y(2:(maxl + 1)**2) = CMPLX(0.0, 0.0)
201 8286900 : lm = 0
202 : !$OMP PARALLEL default(none) &
203 : !$OMP private(l, M, lm, ic1, ic2, cexp) &
204 8286936 : !$OMP shared(ishell, conv, g, y, maxl, structconst, atoms, ikpt, ki, fmpi)
205 :
206 : !$OMP DO schedule(dynamic)
207 : DO l = 0, maxl
208 : lm = l**2
209 : IF (ishell <= conv(l)) THEN
210 : DO M = -l, l
211 : lm = lm + 1
212 : y(lm) = CONJG(y(lm))*g(l)*ImagUnit**l
213 : END DO
214 : ELSE
215 : y(lm + 1:lm + 2*l + 1) = 0
216 : END IF
217 : END DO
218 : !$OMP END DO
219 :
220 : !$OMP DO schedule(dynamic) collapse(2)
221 : DO ic2 = 1+fmpi%irank, atoms%nat, fmpi%isize
222 : DO ic1 = 1, atoms%nat
223 : IF (ic2 /= 1 .AND. ic1 == ic2) CYCLE
224 : cexp = EXP(ImagUnit*tpi_const*dot_PRODUCT(ki, atoms%taual(:, ic1) - atoms%taual(:, ic2)))
225 : DO lm = 1, (maxl + 1)**2
226 : structconst(lm, ic1, ic2, ikpt) = structconst(lm, ic1, ic2, ikpt) + cexp*y(lm)
227 : END DO
228 : END DO
229 : END DO
230 : !$OMP END DO
231 : !$OMP END PARALLEL
232 :
233 : END DO
234 : #ifdef CPP_MPI
235 36 : call timestart("bcast fouriersum")
236 36 : buf_sz = size(structconst,1) * size(structconst,2)
237 96 : DO ic2 = 1, atoms%nat
238 60 : root = mod(ic2-1, fmpi%isize)
239 96 : call MPI_Bcast(structconst(1,1,ic2,ikpt), buf_sz, MPI_DOUBLE_COMPLEX, root, fmpi%mpi_comm, ierr)
240 : enddo
241 48 : call timestop("bcast fouriersum")
242 : #endif
243 : END DO
244 12 : call timestop("fourierspace sum")
245 : !
246 : ! Add contribution for l=0 to diagonal elements and rescale structure constants
247 : !
248 48 : structconst(1, 1, 1, :) = structconst(1, 1, 1, :) - 5.0/16/SQRT(fpi_const)
249 20 : DO i = 2, atoms%nat
250 26180 : structconst(:, i, i, :) = structconst(:, 1, 1, :)
251 : END DO
252 408 : DO l = 0, 2*hybinp%lexp
253 124752 : structconst(l**2 + 1:(l + 1)**2, :, :, :) = structconst(l**2 + 1:(l + 1)**2, :, :, :)*scale**(l + 1)
254 : END DO
255 :
256 12 : rad = (cell%vol*3/4/pi_const)**(1.0/3) ! Wigner-Seitz radius (rad is recycled)
257 :
258 : ! Calculate accuracy of Gamma-decomposition
259 12 : IF (ALL(abs(kpts%bk) > 1e-12)) THEN
260 0 : a = 1e30 ! ikpt = index of shortest non-zero k-point
261 0 : DO i = 2, kpts%nkpt
262 0 : rdum = SUM(MATMUL(kpts%bk(:, i), cell%bmat)**2)
263 0 : IF (rdum < a) THEN
264 0 : ikpt = i
265 0 : a = rdum
266 : END IF
267 : END DO
268 0 : rdum = norm2(MATMUL(kpts%bk(:, ikpt), cell%bmat))
269 0 : a = 0
270 0 : DO ic2 = 1, atoms%nat
271 0 : DO ic1 = 1, MAX(1, ic2 - 1)
272 : a = a + ABS(structconst(1, ic1, ic2, ikpt) - &
273 : (structconst(1, ic1, ic2, 1) + SQRT(fpi_const)/cell%vol/rdum**2* &
274 : EXP(-CMPLX(0.0, 1.0)*tpi_const*dot_PRODUCT( &
275 0 : kpts%bk(:, ikpt), atoms%taual(:, ic2) - atoms%taual(:, ic1)))))**2
276 : END DO
277 : END DO
278 0 : a = SQRT(a/atoms%nat**2)
279 0 : aa = SQRT(SUM(ABS(structconst(1, :, :, ikpt))**2)/atoms%nat**2)
280 0 : IF (first) WRITE (oUnit, '(/A,F8.5,A,F8.5,A)') 'Accuracy of Gamma-decomposition (structureconstant):', a, ' (abs)', a/aa, ' (rel)'
281 : ENDIF
282 :
283 12 : deallocate (ptsh, radsh)
284 :
285 12 : first = .FALSE.
286 :
287 12 : call timestop("calc struc_const.")
288 12 : END SUBROUTINE structureconstant
289 :
290 : ! -----------------
291 :
292 : ! Determines all shells of the crystal defined by lat and vol with radii smaller than rad.
293 : ! The lattice points (number = nptsh) are stored in ptsh, their corresponding lengths (shell radii) in radsh.
294 :
295 24 : SUBROUTINE getshells(ptsh, nptsh, radsh, nshell, rad, lat, lwrite)
296 :
297 : USE m_juDFT
298 : USE m_constants
299 :
300 : IMPLICIT NONE
301 :
302 : LOGICAL, INTENT(IN) :: lwrite
303 : INTEGER, INTENT(INOUT) :: nptsh, nshell
304 : INTEGER, ALLOCATABLE :: ptsh(:, :)
305 : REAL, ALLOCATABLE :: radsh(:)
306 : REAL, INTENT(IN) :: rad, lat(:, :)
307 : REAL :: r(3), rdum
308 24 : INTEGER, ALLOCATABLE :: pnt(:)
309 : INTEGER :: n, i, ix, iy, iz, ok
310 : LOGICAL :: found
311 24 : INTEGER, ALLOCATABLE :: ihelp(:, :)
312 24 : REAL, ALLOCATABLE :: rhelp(:)
313 :
314 24 : allocate (ptsh(3, 100000), radsh(100000), stat=ok)
315 24 : IF (ok /= 0) call judft_error('getshells: failure allocation ptsh/radsh')
316 :
317 9600024 : ptsh = 0
318 2400024 : radsh = 0
319 :
320 : n = 0
321 : i = 0
322 1800 : DO
323 1824 : found = .FALSE.
324 156952 : DO ix = -n, n
325 9534216 : DO iy = -(n - ABS(ix)), n - ABS(ix)
326 9377264 : iz = n - ABS(ix) - ABS(iy)
327 73791584 : 1 r = ix*lat(:, 1) + iy*lat(:, 2) + iz*lat(:, 3)
328 73791584 : rdum = SUM(r**2)
329 18447896 : IF (rdum < rad**2) THEN
330 3451944 : found = .TRUE.
331 3451944 : i = i + 1
332 3451944 : IF (i > SIZE(radsh)) THEN
333 120 : allocate (rhelp(SIZE(radsh)), ihelp(3, SIZE(ptsh, 2)), stat=ok)
334 24 : IF (ok /= 0) call judft_error('getshells: failure allocation rhelp/ihelp')
335 3200048 : rhelp = radsh
336 12800048 : ihelp = ptsh
337 24 : deallocate (radsh, ptsh)
338 120 : allocate (radsh(SIZE(rhelp) + 100000), ptsh(3, SIZE(ihelp, 2) + 100000), stat=ok)
339 24 : IF (ok /= 0) call judft_error('getshells: failure re-allocation ptsh/radsh')
340 3200024 : radsh(1:SIZE(rhelp)) = rhelp
341 12800024 : ptsh(:, 1:SIZE(ihelp, 2)) = ihelp
342 24 : deallocate (rhelp, ihelp)
343 : END IF
344 13807776 : ptsh(:, i) = [ix, iy, iz]
345 3451944 : radsh(i) = SQRT(rdum)
346 : END IF
347 18603024 : IF (iz > 0) THEN
348 9070632 : iz = -iz
349 9070632 : GOTO 1
350 : END IF
351 : END DO
352 : END DO
353 1824 : IF (.NOT. found) EXIT
354 1800 : n = n + 1
355 : END DO
356 24 : nptsh = i
357 :
358 72 : allocate (pnt(nptsh))
359 :
360 : !reallocate radsh ptsh
361 120 : allocate (rhelp(nptsh), ihelp(3, nptsh))
362 3451992 : rhelp = radsh(1:nptsh)
363 13807824 : ihelp = ptsh(:, 1:nptsh)
364 24 : deallocate (radsh, ptsh)
365 72 : allocate (radsh(nptsh), ptsh(3, nptsh))
366 3451992 : radsh = rhelp
367 13807824 : ptsh = ihelp
368 24 : deallocate (rhelp, ihelp)
369 :
370 6903936 : call sort(pnt, radsh, [(1.0*i,i=1,nptsh)])
371 :
372 6903960 : radsh = radsh(pnt)
373 27615600 : ptsh = ptsh(:, pnt)
374 24 : nshell = 1
375 3451944 : DO i = 2, nptsh
376 3451944 : IF (radsh(i) - radsh(i - 1) > 1e-10) nshell = nshell + 1
377 : END DO
378 :
379 24 : IF (lwrite) &
380 : WRITE (oUnit, '(A,F10.5,A,I7,A,I5,A)') &
381 6 : ' Sphere of radius', rad, ' contains', &
382 12 : nptsh, ' lattice points and', nshell, ' shells.'
383 :
384 24 : END SUBROUTINE getshells
385 :
386 12 : subroutine realspace_sum(atoms, cell, hybinp, fmpi, kpts, first, scale, convpar, g, a, a1, rad, structconst)
387 : use ieee_arithmetic
388 : implicit none
389 : type(t_atoms), intent(in) :: atoms
390 : type(t_cell), intent(in) :: cell
391 : type(t_hybinp), intent(in):: hybinp
392 : TYPE(t_mpi), INTENT(IN) :: fmpi
393 : type(t_kpts), intent(in) :: kpts
394 : logical, intent(in) :: first
395 : real, intent(in) :: rad, scale, convpar(0:2*hybinp%lexp)
396 : real, intent(inout) :: g(0:2*hybinp%lexp), a, a1
397 : complex, intent(inout) :: structconst(:,:,:,:)
398 :
399 : integer :: ic2, ic1, i, ishell, l, m, maxl, lm, ikpt, nptsh, nshell, ierr, s
400 12 : integer :: conv(0:2*hybinp%lexp)
401 12 : integer, allocatable :: pnt(:), ptsh(:,:)
402 : INTEGER, PARAMETER :: ADDSHELL1 = 40
403 :
404 : real :: rdum, rexp, ra(3), rc(3), tmp_vec(3)
405 12 : real, allocatable :: radsh(:)
406 :
407 12 : complex :: shlp((2*hybinp%lexp + 1)**2, kpts%nkpt)
408 12 : COMPLEX :: cdum, cexp, y((2*hybinp%lexp + 1)**2)
409 :
410 12 : rdum = cell%vol**(1.0/3)
411 :
412 : !
413 : ! Determine atomic shells
414 12 : call timestart("determine atomic shell")
415 12 : CALL getshells(ptsh, nptsh, radsh, nshell, rad, cell%amat, first)
416 12 : call timestop("determine atomic shell")
417 :
418 36 : allocate (pnt(nptsh))
419 117828 : structconst = 0
420 12 : a1 = 0
421 :
422 : !
423 : ! Real-space sum
424 : !
425 12 : call timestart("realspace sum")
426 22 : DO ic2 = 1+fmpi%irank, atoms%nat, fmpi%isize
427 : ! !$OMP PARALLEL DO default(none) &
428 : ! !$OMP shared(ic2, atoms, cell, nptsh, structconst, hybinp, kpts, scale, convpar) &
429 : ! !$OMP private(ic1, tmp_vec, i, ra, rc, a, pnt, maxl, l, conv, shlp, ishell, rexp, g, y) &
430 : ! !$OMP private(rdum, cexp, lm, cdum)&
431 : ! !$OMP firstprivate(ptsh, radsh) &
432 : ! !$OMP reduction(max:a1)
433 40 : DO ic1 = 1, atoms%nat
434 18 : IF (ic2 /= 1 .AND. ic1 == ic2) CYCLE
435 : !MATMUL(cell%amat, (atoms%taual(:, ic2) - atoms%taual(:, ic1)))
436 56 : tmp_vec = atoms%taual(:, ic2) - atoms%taual(:, ic1)
437 14 : call dgemv("N", 3, 3, 1.0, cell%amat, 3, tmp_vec, 1, 0.0, rc, 1)
438 552380 : DO i = 1, nptsh
439 : !ra = MATMUL(cell%amat, ptsh(:, i)) + rc
440 2209464 : tmp_vec = real(ptsh(:, i))
441 552366 : call dgemv("N", 3, 3, 1.0, cell%amat, 3, tmp_vec, 1, 0.0, ra, 1)
442 2209464 : ra = ra + rc
443 2209478 : radsh(i) = norm2(ra)
444 : END DO
445 1104760 : call sort(pnt, radsh, [(1.0*s, s=1,nptsh)])
446 4418956 : ptsh = ptsh(:, pnt)
447 1104774 : radsh = radsh(pnt)
448 14 : maxl = 2*hybinp%lexp
449 14 : a1 = 1e30 ! stupid initial value
450 14 : ishell = 1
451 476 : conv = HUGE(i)
452 45794 : shlp = 0
453 536458 : DO i = 1, nptsh
454 3026520 : IF (ALL(conv /= HUGE(i))) EXIT
455 536444 : IF (i /= 1) THEN
456 536430 : IF (ABS(radsh(i) - radsh(i - 1)) > 1e-10) ishell = ishell + 1
457 : ENDIF
458 : !ra = MATMUL(cell%amat, ptsh(:, i)) + rc
459 2145776 : tmp_vec = real(ptsh(:, i))
460 536444 : call dgemv("N", 3, 3, 1.0, cell%amat, 3, tmp_vec, 1, 0.0, ra, 1)
461 2145776 : ra = ra + rc
462 :
463 2145776 : a = scale*norm2(ra)
464 536444 : IF (abs(a) < 1e-12) THEN
465 : CYCLE
466 536438 : ELSE IF (ABS(a - a1) > 1e-10) THEN
467 4438 : a1 = a
468 4438 : rexp = EXP(-a)
469 4438 : IF (ishell <= conv(0)) g(0) = rexp/a &
470 3036 : *(1 + a*11/16*(1 + a*3/11*(1 + a/9)))
471 4438 : IF (ishell <= conv(1)) g(1) = rexp/a**2 &
472 2974 : *(1 + a*(1 + a/2*(1 + a*7/24*(1 + a/7))))
473 4438 : IF (ishell <= conv(2)) g(2) = rexp/a**3 &
474 2976 : *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a*3/16*(1 + a/9))))))
475 4438 : IF (ishell <= conv(3)) g(3) = rexp/a**4 &
476 2718 : *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6*(1 + a/8)))))))
477 4438 : IF (ishell <= conv(4)) g(4) = rexp/a**5 &
478 : *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6*(1 + a/7*(1 + a/8 &
479 2594 : *(1 + a/10)))))))))
480 4438 : IF (ishell <= conv(5)) g(5) = rexp/a**6 &
481 : *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6*(1 + a/7*(1 + a/8*(1 + a/9 &
482 2294 : *(1 + a/10))))))))))
483 4438 : IF (ishell <= conv(6)) g(6) = rexp/a**7 &
484 : *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6*(1 + a/7*(1 + a/8*(1 + a/9 &
485 2106 : *(1 + a/10*(1 + a/11*(1 + a/12))))))))))))
486 4438 : IF (ishell <= conv(7)) g(7) = rexp/a**8 &
487 : *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6*(1 + a/7*(1 + a/8*(1 + a/9 &
488 1808 : *(1 + a/10*(1 + a/11*(1 + a/12*(1 + a/13)))))))))))))
489 26556 : DO l = 8, maxl
490 26556 : IF (ishell <= conv(l)) g(l) = a**(-l - 1)
491 : END DO
492 62060 : DO l = 0, maxl
493 62060 : IF (conv(l) == HUGE(i) .AND. g(l) < convpar(l)/10) conv(l) = ishell + ADDSHELL1
494 : END DO
495 : END IF
496 536438 : IF (ishell > conv(maxl) .AND. maxl /= 0) maxl = maxl - 1
497 536438 : call ylm4(maxl, ra, y)
498 584717420 : y = CONJG(y)
499 :
500 2145766 : DO ikpt = 1, kpts%nkpt
501 19378538 : DO l = 0, maxl
502 68931120 : rdum = dot_product(kpts%bk(:, ikpt), ptsh(:, i))
503 17232780 : cexp = EXP(ImagUnit*tpi_const*rdum)
504 17232780 : lm = l**2
505 18842094 : IF (ishell <= conv(l)) THEN
506 9994440 : cdum = cexp*g(l)
507 173057448 : DO M = -l, l
508 163063008 : lm = lm + 1
509 173057448 : shlp(lm, ikpt) = shlp(lm, ikpt) + cdum*y(lm)
510 : END DO
511 : END IF
512 : END DO
513 : END DO
514 : END DO
515 45804 : structconst(:, ic1, ic2, :) = shlp
516 : END DO
517 : ! !$OMP END PARALLEL DO
518 : END DO
519 : #ifdef CPP_MPI
520 12 : call MPI_ALLREDUCE(MPI_IN_PLACE, a1, 1, MPI_DOUBLE_PRECISION, MPI_MAX, fmpi%mpi_comm, ierr)
521 60 : CALL MPI_ALLREDUCE(MPI_IN_PLACE, structconst, size(structconst), MPI_DOUBLE_COMPLEX,MPI_SUM,fmpi%mpi_comm,ierr)
522 : #endif
523 12 : call timestop("realspace sum")
524 12 : deallocate (ptsh, radsh)
525 12 : end subroutine realspace_sum
526 : end module m_structureconstant
|
