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_mkgxyz3
7 : USE m_judft
8 : !-----------------------------------------------------------------------------
9 : ! Using the cartesian components and derivatives of a charge density rho on
10 : ! the real space grid, make the following quantaties:
11 : !
12 : ! gr(js): grad(rho_js)
13 : ! sigma: |grad(rho)|^2 for jspins==1, otherwise three components, namely
14 : ! |grad(rho_up)|^2,grad(rho_up)*grad(rho_down),|grad(rho_down)|^2
15 : ! laplace(js): laplace(rho_js)
16 : !
17 : ! and these older components:
18 : ! agrt/u/d: |grad(rho)| for total density/spin-up/spin-down
19 : ! g2rt/u/d: laplace(rho)
20 : ! gggrt/u/d: (grad(rho))*(grad(|grad(rho)|)) [scalar product]
21 : ! gzgr: (grad(zeta))*(grad(rho)) for zeta=(rho_up-rho_down)/rho
22 : !
23 : ! which are used to calculate gradient contributions to the xc potential and
24 : ! energy.
25 : !
26 : ! vl is rho, dv[i][j] are the partial derivatives along one/two directions.
27 : !
28 : ! Modified so only allocated old quantities require calculations. A.N. 2019
29 : !
30 : ! Quantities fo libxc are calculated as well. D.W./M.R. 2018
31 : !
32 : ! Original script by T.A. 1996
33 : !-----------------------------------------------------------------------------
34 : CONTAINS
35 2984 : SUBROUTINE mkgxyz3(vl,dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvxz,dvxy,idx,grad)
36 : USE m_types
37 : IMPLICIT NONE
38 : REAL, INTENT (IN) :: vl(:,:)
39 : REAL, INTENT (IN) :: dvx(:,:),dvy(:,:),dvz(:,:)
40 : REAL, INTENT (IN) :: dvxx(:,:),dvyy(:,:),dvzz(:,:)
41 : REAL, INTENT (IN) :: dvyz(:,:),dvxz(:,:),dvxy(:,:)
42 : INTEGER ,intent(in) :: idx
43 : TYPE(t_gradients), INTENT(INOUT) :: grad
44 :
45 : REAL vlt,dvxt,dvyt,dvzt,dvxxt,dvyyt,dvzzt,dvyzt,dvxzt,dvxyt, &
46 : vlu,dvxu,dvyu,dvzu,dvxxu,dvyyu,dvzzu,dvyzu,dvxzu,dvxyu, &
47 : vld,dvxd,dvyd,dvzd,dvxxd,dvyyd,dvzzd,dvyzd,dvxzd,dvxyd, &
48 : dagrxt,dagrxd,dagrxu,dagryt,dagryd,dagryu,dagrzt,dagrzd,&
49 : dagrzu,dzdx,dzdy,dzdz
50 : REAL, PARAMETER :: sml = 1.e-14
51 : INTEGER i,js,jspins,nsp
52 :
53 2984 : nsp=SIZE(dvx,1)
54 2984 : jspins=SIZE(dvx,2)
55 :
56 : ! Gradients for libxc and sourcefree calculations.
57 2984 : IF (ALLOCATED(grad%gr)) THEN
58 39 : DO js=1,jspins
59 135591 : DO i=1,nsp
60 542235 : grad%gr(:,i+idx,js)=(/dvx(i,js),dvy(i,js),dvz(i,js)/)
61 : END DO
62 : END DO
63 : ! Use contracted gradients only for libxc.
64 12 : IF(ALLOCATED(grad%sigma)) THEN
65 12 : IF (jspins==1) THEN
66 41475 : DO i=1,nsp
67 41475 : grad%sigma(1,i+idx) = dvx(i,1)*dvx(i,1) + dvy(i,1)*dvy(i,1) + dvz(i,1)*dvz(i,1)
68 : END DO
69 : ELSE
70 35289 : DO i=1,nsp
71 35280 : grad%sigma(1,i+idx) = dvx(i,1)*dvx(i,1) + dvy(i,1)*dvy(i,1) + dvz(i,1)*dvz(i,1)
72 35280 : grad%sigma(2,i+idx) = dvx(i,1)*dvx(i,2) + dvy(i,1)*dvy(i,2) + dvz(i,1)*dvz(i,2)
73 35289 : grad%sigma(3,i+idx) = dvx(i,2)*dvx(i,2) + dvy(i,2)*dvy(i,2) + dvz(i,2)*dvz(i,2)
74 : END DO
75 : END IF
76 : END IF
77 12 : IF(ALLOCATED(grad%laplace)) THEN
78 39 : DO js=1,jspins
79 135591 : DO i=1,nsp
80 135579 : grad%laplace(i+idx,js)= dvxx(i,js)+dvyy(i,js)+dvzz(i,js)
81 : END DO
82 : END DO
83 : END IF
84 12 : RETURN ! Do not calculate arrays for in-build GGA.
85 : END IF
86 :
87 8916 : IF (ANY(SHAPE(vl).NE.SHAPE(dvx))) CALL judft_error("Gradients for internal GGA called with inconsistent sizes",hint="This is a bug")
88 :
89 2972 : IF(ALLOCATED(grad%agrt)) THEN
90 9964713 : DO i = 1,nsp
91 9961741 : grad%agrt(idx+i) = 0.0
92 9961741 : grad%agru(idx+i) = 0.0
93 9961741 : grad%agrd(idx+i) = 0.0
94 9961741 : grad%gggrt(idx+i) = 0.0
95 9961741 : grad%gggru(idx+i) = 0.0
96 9961741 : grad%gggrd(idx+i) = 0.0
97 9961741 : grad%gzgr(idx+i) = 0.0
98 9961741 : grad%g2rt(idx+i) = 0.0
99 9961741 : grad%g2ru(idx+i) = 0.0
100 9964713 : grad%g2rd(idx+i) = 0.0
101 : END DO
102 :
103 2972 : IF (jspins.eq.1) THEN
104 :
105 8641275 : DO i = 1,nsp
106 :
107 8638845 : vlu = max(vl(i,1)/2,sml)
108 8638845 : dvxu = dvx(i,1)/2
109 8638845 : dvyu = dvy(i,1)/2
110 8638845 : dvzu = dvz(i,1)/2
111 8638845 : dvxxu = dvxx(i,1)/2
112 8638845 : dvyyu = dvyy(i,1)/2
113 8638845 : dvzzu = dvzz(i,1)/2
114 8638845 : dvyzu = dvyz(i,1)/2
115 8638845 : dvxzu = dvxz(i,1)/2
116 8638845 : dvxyu = dvxy(i,1)/2
117 :
118 8638845 : vld = vlu
119 8638845 : dvxd = dvxu
120 8638845 : dvyd = dvyu
121 8638845 : dvzd = dvzu
122 8638845 : dvxxd = dvxxu
123 8638845 : dvyyd = dvyyu
124 8638845 : dvzzd = dvzzu
125 8638845 : dvyzd = dvyzu
126 8638845 : dvxzd = dvxzu
127 8638845 : dvxyd = dvxyu
128 :
129 8638845 : vlt = vlu + vld
130 8638845 : dvxt = dvxu + dvxd
131 8638845 : dvyt = dvyu + dvyd
132 8638845 : dvzt = dvzu + dvzd
133 8638845 : dvxxt = dvxxu + dvxxd
134 8638845 : dvyyt = dvyyu + dvyyd
135 8638845 : dvzzt = dvzzu + dvzzd
136 8638845 : dvyzt = dvyzu + dvyzd
137 8638845 : dvxzt = dvxzu + dvxzd
138 8638845 : dvxyt = dvxyu + dvxyd
139 :
140 8638845 : grad%agrt(idx+i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
141 8638845 : grad%agru(idx+i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
142 8638845 : grad%agrd(idx+i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
143 :
144 8638845 : dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvxzt) / grad%agrt(idx+i)
145 8638845 : dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvxzu) / grad%agru(idx+i)
146 8638845 : dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvxzd) / grad%agrd(idx+i)
147 :
148 8638845 : dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(idx+i)
149 8638845 : dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(idx+i)
150 8638845 : dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(idx+i)
151 :
152 8638845 : dagrzt = (dvxt*dvxzt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(idx+i)
153 8638845 : dagrzu = (dvxu*dvxzu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(idx+i)
154 8638845 : dagrzd = (dvxd*dvxzd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(idx+i)
155 :
156 8638845 : grad%gggrt(idx+i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
157 8638845 : grad%gggru(idx+i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
158 8638845 : grad%gggrd(idx+i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
159 :
160 8638845 : dzdx = (dvxu-dvxd)/vlt - (vlu-vld)*dvxt/vlt**2
161 8638845 : dzdy = (dvyu-dvyd)/vlt - (vlu-vld)*dvyt/vlt**2
162 8638845 : dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2
163 :
164 8638845 : grad%gzgr(idx+i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
165 :
166 8638845 : grad%g2rt(idx+i) = dvxxt + dvyyt + dvzzt
167 8638845 : grad%g2ru(idx+i) = dvxxu + dvyyu + dvzzu
168 8641275 : grad%g2rd(idx+i) = dvxxd + dvyyd + dvzzd
169 : ENDDO
170 : ELSE
171 1323438 : DO i = 1,nsp
172 :
173 1322896 : vlu = max(vl(i,1),sml)
174 1322896 : dvxu = dvx(i,1)
175 1322896 : dvyu = dvy(i,1)
176 1322896 : dvzu = dvz(i,1)
177 1322896 : dvxxu = dvxx(i,1)
178 1322896 : dvyyu = dvyy(i,1)
179 1322896 : dvzzu = dvzz(i,1)
180 1322896 : dvyzu = dvyz(i,1)
181 1322896 : dvxzu = dvxz(i,1)
182 1322896 : dvxyu = dvxy(i,1)
183 :
184 1322896 : vld = max(vl(i,jspins),sml)
185 1322896 : dvxd = dvx(i,jspins)
186 1322896 : dvyd = dvy(i,jspins)
187 1322896 : dvzd = dvz(i,jspins)
188 1322896 : dvxxd = dvxx(i,jspins)
189 1322896 : dvyyd = dvyy(i,jspins)
190 1322896 : dvzzd = dvzz(i,jspins)
191 1322896 : dvyzd = dvyz(i,jspins)
192 1322896 : dvxzd = dvxz(i,jspins)
193 1322896 : dvxyd = dvxy(i,jspins)
194 :
195 1322896 : vlt = vlu + vld
196 :
197 1322896 : dvxt = dvxu + dvxd
198 1322896 : dvyt = dvyu + dvyd
199 1322896 : dvzt = dvzu + dvzd
200 1322896 : dvxxt = dvxxu + dvxxd
201 1322896 : dvyyt = dvyyu + dvyyd
202 1322896 : dvzzt = dvzzu + dvzzd
203 1322896 : dvyzt = dvyzu + dvyzd
204 1322896 : dvxzt = dvxzu + dvxzd
205 1322896 : dvxyt = dvxyu + dvxyd
206 :
207 1322896 : grad%agrt(idx+i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
208 1322896 : grad%agru(idx+i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
209 1322896 : grad%agrd(idx+i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
210 :
211 1322896 : dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvxzt) / grad%agrt(idx+i)
212 1322896 : dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvxzu) / grad%agru(idx+i)
213 1322896 : dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvxzd) / grad%agrd(idx+i)
214 :
215 1322896 : dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(idx+i)
216 1322896 : dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(idx+i)
217 1322896 : dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(idx+i)
218 :
219 1322896 : dagrzt = (dvxt*dvxzt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(idx+i)
220 1322896 : dagrzu = (dvxu*dvxzu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(idx+i)
221 1322896 : dagrzd = (dvxd*dvxzd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(idx+i)
222 :
223 1322896 : grad%gggrt(idx+i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
224 1322896 : grad%gggru(idx+i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
225 1322896 : grad%gggrd(idx+i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
226 :
227 1322896 : dzdx = (dvxu-dvxd)/vlt - (vlu-vld)*dvxt/vlt**2
228 1322896 : dzdy = (dvyu-dvyd)/vlt - (vlu-vld)*dvyt/vlt**2
229 1322896 : dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2
230 :
231 1322896 : grad%gzgr(idx+i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
232 :
233 1322896 : grad%g2rt(idx+i) = dvxxt + dvyyt + dvzzt
234 1322896 : grad%g2ru(idx+i) = dvxxu + dvyyu + dvzzu
235 1323438 : grad%g2rd(idx+i) = dvxxd + dvyyd + dvzzd
236 :
237 : END DO
238 : END IF
239 : END IF
240 :
241 : END SUBROUTINE mkgxyz3
242 : END MODULE m_mkgxyz3
|
