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 : !c by use of cartesian x,y,z components of charge density gradients,
10 : !c make the quantities
11 : !cc agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd,gzgr
12 : !cc used to calculate gradient contribution to xc potential and
13 : !cc energy.
14 : !c.....------------------------------------------------------------------
15 : CONTAINS
16 681 : SUBROUTINE mkgxyz3(vl,dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvzx,dvxy,grad)
17 : USE m_types
18 : IMPLICIT NONE
19 : REAL, INTENT (IN) :: vl(:,:)
20 : REAL, INTENT (IN) :: dvx(:,:),dvy(:,:),dvz(:,:)
21 : REAL, INTENT (IN) :: dvxx(:,:),dvyy(:,:),dvzz(:,:)
22 : REAL, INTENT (IN) :: dvyz(:,:),dvzx(:,:),dvxy(:,:)
23 :
24 : TYPE(t_gradients),INTENT(INOUT)::grad
25 :
26 : REAL vlt,dvxt,dvyt,dvzt,dvxxt,dvyyt,dvzzt,dvyzt,dvzxt,dvxyt,&
27 : vlu,dvxu,dvyu,dvzu,dvxxu,dvyyu,dvzzu,dvyzu,dvzxu,dvxyu,&
28 : vld,dvxd,dvyd,dvzd,dvxxd,dvyyd,dvzzd,dvyzd,dvzxd,dvxyd,&
29 : dagrxt,dagrxd,dagrxu,dagryt,dagryd,dagryu,dagrzt,dagrzd,&
30 : dagrzu,dzdx,dzdy,dzdz
31 : REAL, PARAMETER :: sml = 1.e-14
32 : INTEGER i,js,jspins,nsp
33 :
34 681 : nsp=SIZE(dvx,1)
35 681 : jspins=SIZE(dvx,2)
36 :
37 681 : IF (ALLOCATED(grad%gr)) THEN
38 : ! Gradients for libxc
39 240 : DO js=1,jspins
40 1024140 : DO i=1,nsp
41 512100 : grad%gr(:,i,js)=(/dvx(i,js),dvy(i,js),dvz(i,js)/)
42 : ENDDO
43 : END DO
44 40 : IF(ALLOCATED(grad%sigma)) THEN
45 : !Use only contracted gradients for libxc
46 40 : IF (jspins==1) THEN
47 0 : DO i=1,nsp
48 0 : grad%sigma(1,i) = dvx(i,1)*dvx(i,1) + dvy(i,1)*dvy(i,1) + dvz(i,1)*dvz(i,1)
49 : ENDDO
50 : ELSE
51 409640 : DO i=1,nsp
52 204800 : grad%sigma(1,i) = dvx(i,1)*dvx(i,1) + dvy(i,1)*dvy(i,1) + dvz(i,1)*dvz(i,1)
53 204800 : grad%sigma(2,i) = dvx(i,1)*dvx(i,2) + dvy(i,1)*dvy(i,2) + dvz(i,1)*dvz(i,2)
54 204840 : grad%sigma(3,i) = dvx(i,2)*dvx(i,2) + dvy(i,2)*dvy(i,2) + dvz(i,2)*dvz(i,2)
55 : ENDDO
56 : ENDIF
57 : END IF
58 40 : IF(ALLOCATED(grad%laplace)) THEN
59 240 : DO js=1,jspins
60 1024140 : DO i=1,nsp
61 512100 : grad%laplace(i,js)= dvxx(i,js)+dvyy(i,js)+dvzz(i,js)
62 : ENDDO
63 : ENDDO
64 : ENDIF
65 : RETURN
66 : ENDIF
67 :
68 641 : IF (ANY(SHAPE(vl).NE.SHAPE(dvx))) CALL judft_error("Gradients for internal GGA called with inconsistent sizes",hint="This is a bug")
69 :
70 1653082 : DO i = 1,size(grad%agrt)
71 1652441 : grad%agrt(i) = 0.0
72 1652441 : grad%agru(i) = 0.0
73 1652441 : grad%agrd(i) = 0.0
74 1652441 : grad%gggrt(i) = 0.0
75 1652441 : grad%gggru(i) = 0.0
76 1652441 : grad%gggrd(i) = 0.0
77 1652441 : grad%gzgr(i) = 0.0
78 1652441 : grad%g2rt(i) = 0.0
79 1652441 : grad%g2ru(i) = 0.0
80 1653082 : grad%g2rd(i) = 0.0
81 : ENDDO
82 :
83 641 : IF (jspins.eq.1) THEN
84 :
85 1371107 : DO i = 1,nsp
86 :
87 685440 : vlu = max(vl(i,1)/2,sml)
88 685440 : dvxu = dvx(i,1)/2
89 685440 : dvyu = dvy(i,1)/2
90 685440 : dvzu = dvz(i,1)/2
91 685440 : dvxxu = dvxx(i,1)/2
92 685440 : dvyyu = dvyy(i,1)/2
93 685440 : dvzzu = dvzz(i,1)/2
94 685440 : dvyzu = dvyz(i,1)/2
95 685440 : dvzxu = dvzx(i,1)/2
96 685440 : dvxyu = dvxy(i,1)/2
97 :
98 685440 : vld = vlu
99 685440 : dvxd = dvxu
100 685440 : dvyd = dvyu
101 685440 : dvzd = dvzu
102 685440 : dvxxd = dvxxu
103 685440 : dvyyd = dvyyu
104 685440 : dvzzd = dvzzu
105 685440 : dvyzd = dvyzu
106 685440 : dvzxd = dvzxu
107 685440 : dvxyd = dvxyu
108 :
109 :
110 685440 : vlt = vlu + vld
111 :
112 685440 : dvxt = dvxu + dvxd
113 685440 : dvyt = dvyu + dvyd
114 685440 : dvzt = dvzu + dvzd
115 685440 : dvxxt = dvxxu + dvxxd
116 685440 : dvyyt = dvyyu + dvyyd
117 685440 : dvzzt = dvzzu + dvzzd
118 685440 : dvyzt = dvyzu + dvyzd
119 685440 : dvzxt = dvzxu + dvzxd
120 685440 : dvxyt = dvxyu + dvxyd
121 :
122 : ! agr: abs(grad(ro)), t,u,d for total, up and down.
123 :
124 685440 : grad%agrt(i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
125 685440 : grad%agru(i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
126 685440 : grad%agrd(i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
127 :
128 685440 : dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvzxt) / grad%agrt(i)
129 685440 : dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvzxu) / grad%agru(i)
130 685440 : dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvzxd) / grad%agrd(i)
131 :
132 685440 : dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(i)
133 685440 : dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(i)
134 685440 : dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(i)
135 :
136 685440 : dagrzt = (dvxt*dvzxt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(i)
137 685440 : dagrzu = (dvxu*dvzxu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(i)
138 685440 : dagrzd = (dvxd*dvzxd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(i)
139 :
140 685440 : grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
141 685440 : grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
142 685440 : grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
143 :
144 : ! dzdx=d(zeta)/dx,..
145 :
146 685440 : dzdx = (dvxu-dvxd)/vlt - (vlu-vld)*dvxt/vlt**2
147 685440 : dzdy = (dvyu-dvyd)/vlt - (vlu-vld)*dvyt/vlt**2
148 685440 : dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2
149 :
150 : ! gzgr=grad(zeta)*grad(ro).
151 :
152 685440 : grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
153 :
154 : ! g2r: grad(grad(ro))
155 :
156 685440 : grad%g2rt(i) = dvxxt + dvyyt + dvzzt
157 685440 : grad%g2ru(i) = dvxxu + dvyyu + dvzzu
158 685667 : grad%g2rd(i) = dvxxd + dvyyd + dvzzd
159 :
160 :
161 : ENDDO
162 :
163 : ELSE
164 :
165 1934416 : DO i = 1,nsp
166 :
167 967001 : vlu = max(vl(i,1),sml)
168 967001 : dvxu = dvx(i,1)
169 967001 : dvyu = dvy(i,1)
170 967001 : dvzu = dvz(i,1)
171 967001 : dvxxu = dvxx(i,1)
172 967001 : dvyyu = dvyy(i,1)
173 967001 : dvzzu = dvzz(i,1)
174 967001 : dvyzu = dvyz(i,1)
175 967001 : dvzxu = dvzx(i,1)
176 967001 : dvxyu = dvxy(i,1)
177 :
178 967001 : vld = max(vl(i,jspins),sml)
179 967001 : dvxd = dvx(i,jspins)
180 967001 : dvyd = dvy(i,jspins)
181 967001 : dvzd = dvz(i,jspins)
182 967001 : dvxxd = dvxx(i,jspins)
183 967001 : dvyyd = dvyy(i,jspins)
184 967001 : dvzzd = dvzz(i,jspins)
185 967001 : dvyzd = dvyz(i,jspins)
186 967001 : dvzxd = dvzx(i,jspins)
187 967001 : dvxyd = dvxy(i,jspins)
188 :
189 967001 : vlt = vlu + vld
190 :
191 967001 : dvxt = dvxu + dvxd
192 967001 : dvyt = dvyu + dvyd
193 967001 : dvzt = dvzu + dvzd
194 967001 : dvxxt = dvxxu + dvxxd
195 967001 : dvyyt = dvyyu + dvyyd
196 967001 : dvzzt = dvzzu + dvzzd
197 967001 : dvyzt = dvyzu + dvyzd
198 967001 : dvzxt = dvzxu + dvzxd
199 967001 : dvxyt = dvxyu + dvxyd
200 :
201 : !c agr: abs(grad(ro)), t,u,d for total, up and down.
202 :
203 967001 : grad%agrt(i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
204 967001 : grad%agru(i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
205 967001 : grad%agrd(i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
206 :
207 967001 : dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvzxt) / grad%agrt(i)
208 967001 : dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvzxu) / grad%agru(i)
209 967001 : dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvzxd) / grad%agrd(i)
210 :
211 967001 : dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(i)
212 967001 : dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(i)
213 967001 : dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(i)
214 :
215 967001 : dagrzt = (dvxt*dvzxt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(i)
216 967001 : dagrzu = (dvxu*dvzxu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(i)
217 967001 : dagrzd = (dvxd*dvzxd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(i)
218 :
219 967001 : grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
220 967001 : grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
221 967001 : grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
222 :
223 : !c dzdx=d(zeta)/dx,..
224 :
225 967001 : dzdx = (dvxu-dvxd)/vlt - (vlu-vld)*dvxt/vlt**2
226 967001 : dzdy = (dvyu-dvyd)/vlt - (vlu-vld)*dvyt/vlt**2
227 967001 : dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2
228 :
229 : !c gzgr=grad(zeta)*grad(ro).
230 :
231 967001 : grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
232 :
233 : !c g2r: grad(grad(ro))
234 :
235 967001 : grad%g2rt(i) = dvxxt + dvyyt + dvzzt
236 967001 : grad%g2ru(i) = dvxxu + dvyyu + dvzzu
237 967415 : grad%g2rd(i) = dvxxd + dvyyd + dvzzd
238 :
239 : ENDDO
240 :
241 : ENDIF
242 :
243 : END SUBROUTINE mkgxyz3
244 : END MODULE m_mkgxyz3
|
