Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : MODULE m_mt_tofrom_grid
7 : USE m_types
8 : PRIVATE
9 : REAL, PARAMETER :: d_15 = 1.e-15
10 : INTEGER, PARAMETER :: ndvgrd = 6 ! this should be consistent across GGA derivative routines
11 : REAL, ALLOCATABLE :: ylh(:, :, :), ylht(:, :, :), ylhtt(:, :, :)
12 : REAL, ALLOCATABLE :: ylhf(:, :, :), ylhff(:, :, :), ylhtf(:, :, :)
13 : REAL, ALLOCATABLE :: wt(:), rx(:, :), thet(:)
14 : PUBLIC :: init_mt_grid, mt_to_grid, mt_from_grid, finish_mt_grid
15 : CONTAINS
16 340 : SUBROUTINE init_mt_grid(jspins, atoms, sphhar, xcpot, sym)
17 : USE m_gaussp
18 : USE m_lhglptg
19 : USE m_lhglpts
20 : IMPLICIT NONE
21 : INTEGER, INTENT(IN) :: jspins
22 : TYPE(t_atoms), INTENT(IN) :: atoms
23 : TYPE(t_sphhar), INTENT(IN) :: sphhar
24 : CLASS(t_xcpot), INTENT(IN) :: xcpot
25 : TYPE(t_sym), INTENT(IN) :: sym
26 :
27 : ! generate nspd points on a sherical shell with radius 1.0
28 : ! angular mesh equidistant in phi,
29 : ! theta are zeros of the legendre polynomials
30 340 : ALLOCATE (wt(atoms%nsp()), rx(3, atoms%nsp()), thet(atoms%nsp()))
31 340 : CALL gaussp(atoms%lmaxd, rx, wt)
32 : ! generate the lattice harmonics on the angular mesh
33 340 : ALLOCATE (ylh(atoms%nsp(), 0:sphhar%nlhd, sphhar%ntypsd))
34 340 : IF (xcpot%needs_grad()) THEN
35 322 : ALLOCATE (ylht, MOLD=ylh)
36 322 : ALLOCATE (ylhtt, MOLD=ylh)
37 322 : ALLOCATE (ylhf, MOLD=ylh)
38 322 : ALLOCATE (ylhff, MOLD=ylh)
39 322 : ALLOCATE (ylhtf, MOLD=ylh)
40 :
41 : CALL lhglptg(sphhar, atoms, rx, atoms%nsp(), xcpot, sym, &
42 322 : ylh, thet, ylht, ylhtt, ylhf, ylhff, ylhtf)
43 : ELSE
44 18 : CALL lhglpts(sphhar, atoms, rx, atoms%nsp(), sym, ylh)
45 : END IF
46 : !ENDIF
47 340 : END SUBROUTINE init_mt_grid
48 :
49 369 : SUBROUTINE mt_to_grid(xcpot, jspins, atoms, sphhar, den_mt, n, grad, ch)
50 : USE m_grdchlh
51 : USE m_mkgylm
52 : IMPLICIT NONE
53 : CLASS(t_xcpot), INTENT(IN) :: xcpot
54 : TYPE(t_atoms), INTENT(IN) :: atoms
55 : TYPE(t_sphhar), INTENT(IN) :: sphhar
56 : REAL, INTENT(IN) :: den_mt(:, 0:, :)
57 : INTEGER, INTENT(IN) :: n, jspins
58 : REAL, INTENT(OUT), OPTIONAL :: ch(:, :)
59 : TYPE(t_gradients), INTENT(INOUT):: grad
60 :
61 369 : REAL, ALLOCATABLE :: chlh(:, :, :), chlhdr(:, :, :), chlhdrr(:, :, :)
62 1107 : REAL, ALLOCATABLE :: chdr(:, :), chdt(:, :), chdf(:, :), ch_tmp(:, :)
63 1107 : REAL, ALLOCATABLE :: chdrr(:, :), chdtt(:, :), chdff(:, :), chdtf(:, :)
64 738 : REAL, ALLOCATABLE :: chdrt(:, :), chdrf(:, :)
65 : INTEGER:: nd, lh, js, jr, kt, k, nsp
66 :
67 369 : nd = atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1)
68 369 : nsp = atoms%nsp()
69 :
70 369 : ALLOCATE (chlh(atoms%jmtd, 0:sphhar%nlhd, jspins))
71 369 : ALLOCATE (ch_tmp(nsp, jspins))
72 369 : IF (xcpot%needs_grad()) THEN
73 : ALLOCATE (chdr(nsp, jspins), chdt(nsp, jspins), chdf(nsp, jspins), chdrr(nsp, jspins), &
74 : chdtt(nsp, jspins), chdff(nsp, jspins), chdtf(nsp, jspins), chdrt(nsp, jspins), &
75 351 : chdrf(nsp, jspins))
76 351 : ALLOCATE (chlhdr(atoms%jmtd, 0:sphhar%nlhd, jspins))
77 351 : ALLOCATE (chlhdrr(atoms%jmtd, 0:sphhar%nlhd, jspins))
78 : ENDIF
79 :
80 6045 : DO lh = 0, sphhar%nlh(nd)
81 : ! calculates gradients of radial charge densities of l=> 0.
82 : ! rho*ylh/r**2 is charge density. chlh=rho/r**2.
83 : ! charge density=sum(chlh*ylh).
84 : ! chlhdr=d(chlh)/dr, chlhdrr=dd(chlh)/drr.
85 :
86 26973 : DO js = 1, jspins
87 6257288 : DO jr = 1, atoms%jri(n)
88 6257288 : chlh(jr, lh, js) = den_mt(jr, lh, js)/(atoms%rmsh(jr, n)*atoms%rmsh(jr, n))
89 : ENDDO
90 10464 : IF (xcpot%needs_grad()) CALL grdchlh(1, 1, atoms%jri(n), atoms%dx(n), atoms%rmsh(1, n), &
91 15852 : chlh(1, lh, js), ndvgrd, chlhdr(1, lh, js), chlhdrr(1, lh, js))
92 :
93 : ENDDO ! js
94 : ENDDO ! lh
95 :
96 369 : kt = 0
97 224120 : DO jr = 1, atoms%jri(n)
98 : ! charge density (on extended grid for all jr)
99 : ! following are at points on jr-th sphere.
100 639841 : ch_tmp(:, :) = 0.0
101 : ! generate the densities on an angular mesh
102 1055931 : DO js = 1, jspins
103 6886665 : DO lh = 0, sphhar%nlh(nd)
104 2324514066 : DO k = 1, nsp
105 1165172400 : ch_tmp(k, js) = ch_tmp(k, js) + ylh(k, lh, nd)*chlh(jr, lh, js)
106 : ENDDO
107 : ENDDO
108 : ENDDO
109 223751 : IF (xcpot%needs_grad()) THEN
110 206921 : chdr(:, :) = 0.0 ! d(ch)/dr
111 206921 : chdt(:, :) = 0.0 ! d(ch)/dtheta
112 206921 : chdf(:, :) = 0.0 ! d(ch)/dfai
113 206921 : chdrr(:, :) = 0.0 ! dd(ch)/drr
114 206921 : chdtt(:, :) = 0.0 ! dd(ch)/dtt
115 206921 : chdff(:, :) = 0.0 ! dd(ch)/dff
116 206921 : chdtf(:, :) = 0.0 ! dd(ch)/dtf
117 206921 : chdrt(:, :) = 0.0 ! d(d(ch)/dr)dt
118 206921 : chdrf(:, :) = 0.0 ! d(d(ch)/dr)df
119 : ! generate the derivatives on an angular mesh
120 1005441 : DO js = 1, jspins
121 6583725 : DO lh = 0, sphhar%nlh(nd)
122 : !
123 2135332696 : DO k = 1, nsp
124 1064677576 : chdr(k, js) = chdr(k, js) + ylh(k, lh, nd)*chlhdr(jr, lh, js)
125 1070655120 : chdrr(k, js) = chdrr(k, js) + ylh(k, lh, nd)*chlhdrr(jr, lh, js)
126 : ENDDO
127 :
128 2135731956 : DO k = 1, nsp
129 1064677576 : chdrt(k, js) = chdrt(k, js) + ylht(k, lh, nd)*chlhdr(jr, lh, js)
130 1064677576 : chdrf(k, js) = chdrf(k, js) + ylhf(k, lh, nd)*chlhdr(jr, lh, js)
131 1064677576 : chdt(k, js) = chdt(k, js) + ylht(k, lh, nd)*chlh(jr, lh, js)
132 1064677576 : chdf(k, js) = chdf(k, js) + ylhf(k, lh, nd)*chlh(jr, lh, js)
133 1064677576 : chdtt(k, js) = chdtt(k, js) + ylhtt(k, lh, nd)*chlh(jr, lh, js)
134 1064677576 : chdff(k, js) = chdff(k, js) + ylhff(k, lh, nd)*chlh(jr, lh, js)
135 1070655120 : chdtf(k, js) = chdtf(k, js) + ylhtf(k, lh, nd)*chlh(jr, lh, js)
136 : ENDDO
137 : ENDDO ! lh
138 : ENDDO ! js
139 :
140 : CALL mkgylm(jspins, atoms%rmsh(jr, n), thet, nsp, &
141 206921 : ch_tmp, chdr, chdt, chdf, chdrr, chdtt, chdff, chdtf, chdrt, chdrf, grad, kt)
142 : ENDIF
143 : !Set charge to minimum value
144 223751 : IF (PRESENT(ch)) THEN
145 549441 : WHERE (ABS(ch_tmp) < d_15) ch_tmp = d_15
146 201151 : ch(kt + 1:kt + nsp, :) = ch_tmp(:nsp, :)
147 : ENDIF
148 224120 : kt = kt + nsp
149 : END DO
150 :
151 369 : END SUBROUTINE mt_to_grid
152 :
153 1027 : SUBROUTINE mt_from_grid(atoms, sphhar, n, jspins, v_in, vr)
154 : IMPLICIT NONE
155 : TYPE(t_atoms), INTENT(IN) :: atoms
156 : TYPE(t_sphhar), INTENT(IN):: sphhar
157 : INTEGER, INTENT(IN) :: jspins, n
158 : REAL, INTENT(IN) :: v_in(:, :)
159 : REAL, INTENT(INOUT) :: vr(:, 0:, :)
160 :
161 2054 : REAL :: vpot(atoms%nsp()), vlh
162 : INTEGER :: js, kt, lh, jr, nd, nsp
163 :
164 1027 : nsp = atoms%nsp()
165 1027 : nd = atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1)
166 :
167 2640 : DO js = 1, jspins
168 : !
169 1613 : kt = 0
170 968171 : DO jr = 1, atoms%jri(n)
171 965531 : vpot = v_in(kt + 1:kt + nsp, js)*wt(:)! multiplicate v_in with the weights of the k-points
172 :
173 16118855 : DO lh = 0, sphhar%nlh(nd)
174 : !
175 : ! ---> determine the corresponding potential number
176 : !c through gauss integration
177 : !
178 15153324 : vlh = dot_PRODUCT(vpot(:), ylh(:nsp, lh, nd))
179 16118855 : vr(jr, lh, js) = vr(jr, lh, js) + vlh
180 : ENDDO ! lh
181 967144 : kt = kt + nsp
182 : ENDDO ! jr
183 : ENDDO
184 :
185 1027 : END SUBROUTINE mt_from_grid
186 :
187 340 : SUBROUTINE finish_mt_grid()
188 340 : DEALLOCATE (ylh, wt, rx, thet)
189 340 : IF (ALLOCATED(ylht)) DEALLOCATE (ylht, ylhtt, ylhf, ylhff, ylhtf)
190 340 : END SUBROUTINE finish_mt_grid
191 :
192 : END MODULE m_mt_tofrom_grid
|
