Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2018 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_metagga
7 : PUBLIC :: calc_EnergyDen
8 : PRIVATE :: calc_EnergyDen_auxillary_weights, &
9 : calc_kinEnergyDen_pw, &
10 : calc_kinEnergyDen_mt
11 :
12 : type t_RS_potden
13 : REAL, ALLOCATABLE :: is(:,:), mt(:,:)
14 : end type t_RS_potden
15 :
16 : TYPE t_kinED
17 : logical :: set=.false.
18 : real, allocatable :: is(:,:) ! (nsp*jmtd, jspins)
19 : real, allocatable :: mt(:,:,:) ! (nsp*jmtd, jspins, local num of types)
20 : contains
21 : procedure :: alloc_mt => kED_alloc_mt
22 : END TYPE t_kinED
23 :
24 : CONTAINS
25 :
26 0 : SUBROUTINE kED_alloc_mt(kED,nsp_x_jmtd, jspins, n_start, n_types, n_stride)
27 : IMPLICIT NONE
28 : class(t_kinED), intent(inout) :: kED
29 : integer, intent(in) :: nsp_x_jmtd, jspins, n_start, n_types, n_stride
30 : integer :: cnt, n
31 :
32 0 : if(.not. allocated(kED%mt)) then
33 0 : cnt = 0
34 0 : do n = n_start,n_types,n_stride
35 0 : cnt = cnt + 1
36 : enddo
37 0 : allocate(kED%mt(nsp_x_jmtd, jspins, cnt), source=0.0)
38 : endif
39 0 : END SUBROUTINE kED_alloc_mt
40 :
41 :
42 0 : SUBROUTINE calc_kinEnergyDen_pw(EnergyDen_rs, vTot_rs, den_rs, kinEnergyDen_RS)
43 : USE m_juDFT_stop
44 : !use m_cdngen
45 : IMPLICIT NONE
46 : REAL, INTENT(in) :: den_RS(:,:), EnergyDen_RS(:,:), vTot_RS(:,:)
47 : REAL, INTENT(inout), allocatable :: kinEnergyDen_RS(:,:)
48 : #ifdef CPP_LIBXC
49 : REAL, PARAMETER :: eps = 1e-15
50 :
51 0 : kinEnergyDen_RS = EnergyDen_RS - vTot_RS * den_RS
52 : #else
53 : CALL juDFT_error("MetaGGA require LibXC",hint="compile Fleur with LibXC (e.g. by giving '-external libxc' to ./configure")
54 : #endif
55 0 : END SUBROUTINE calc_kinEnergyDen_pw
56 :
57 0 : SUBROUTINE calc_kinEnergyDen_mt(EnergyDen_RS, vTot_rs, vTot0_rs, core_den_rs, val_den_rs, &
58 0 : kinEnergyDen_RS)
59 : USE m_juDFT_stop
60 : implicit none
61 : REAL, INTENT(in) :: EnergyDen_RS(:,:), vTot_rs(:,:), vTot0_rs(:,:), core_den_rs(:,:), val_den_rs(:,:)
62 : REAL, INTENT(inout) :: kinEnergyDen_RS(:,:)
63 :
64 : #ifdef CPP_LIBXC
65 0 : kinEnergyDen_RS = EnergyDen_RS - (vTot0_rs * core_den_rs + vTot_rs * val_den_rs)
66 : #else
67 : CALL juDFT_error("MetaGGA require LibXC",hint="compile Fleur with LibXC (e.g. by giving '-external libxc' to ./configure")
68 : #endif
69 0 : END SUBROUTINE calc_kinEnergyDen_mt
70 :
71 :
72 0 : SUBROUTINE calc_EnergyDen(eig_id, fmpi, kpts, noco,nococonv, input, banddos, cell, atoms, enpara, stars, &
73 : vacuum, sphhar, sym, gfinp, hub1inp, vTot, results, EnergyDen)
74 : ! calculates the energy density
75 : ! EnergyDen = \sum_i n_i(r) \varepsilon_i
76 : ! where n_i(r) is the one-particle density
77 : ! and \varepsilon_i are the eigenenergies
78 :
79 : USE m_types_setup
80 : USE m_types_potden
81 : USE m_types_kpts
82 : USE m_types_mpi
83 : USE m_types_enpara
84 : USE m_types_misc
85 : USE m_types_regionCharges
86 : USE m_types_dos
87 : USE m_types_vacdos
88 : USE m_types_cdnval
89 : USE m_cdnval
90 : use m_types_nococonv
91 : IMPLICIT NONE
92 :
93 : INTEGER, INTENT(in) :: eig_id
94 : TYPE(t_mpi), INTENT(in) :: fmpi
95 : TYPE(t_kpts), INTENT(in) :: kpts
96 : TYPE(t_noco), INTENT(in) :: noco
97 : TYPE(t_nococonv), INTENT(in) :: nococonv
98 : TYPE(t_input), INTENT(in) :: input
99 : TYPE(t_banddos), INTENT(in) :: banddos
100 : TYPE(t_cell), INTENT(in) :: cell
101 : TYPE(t_atoms), INTENT(in) :: atoms
102 : TYPE(t_enpara), INTENT(in) :: enpara
103 : TYPE(t_stars), INTENT(in) :: stars
104 : TYPE(t_vacuum), INTENT(in) :: vacuum
105 :
106 : TYPE(t_sphhar), INTENT(in) :: sphhar
107 : TYPE(t_sym), INTENT(in) :: sym
108 : TYPE(t_gfinp), INTENT(in) :: gfinp
109 : TYPE(t_hub1inp), INTENT(in) :: hub1inp
110 : TYPE(t_potden), INTENT(in) :: vTot
111 :
112 : TYPE(t_results), INTENT(in) :: results
113 : TYPE(t_potden), INTENT(inout) :: EnergyDen
114 :
115 : ! local
116 : INTEGER :: jspin
117 :
118 0 : TYPE(t_regionCharges) :: regCharges
119 0 : TYPE(t_dos) :: dos
120 0 : TYPE(t_vacdos) :: vacdos
121 :
122 0 : TYPE(t_moments) :: moments
123 0 : TYPE(t_results) :: tmp_results
124 0 : TYPE(t_cdnvalJob) :: cdnvalJob
125 : TYPE(t_potden) :: aux_den, real_den
126 :
127 :
128 0 : CALL regCharges%init(input, atoms)
129 0 : CALL dos%init(input, atoms, kpts, banddos,results%eig)
130 0 : CALL vacdos%init(input, atoms, kpts, banddos,results%eig)
131 : ! CALL moments%init(input, atoms)
132 0 : CALL moments%init(fmpi,input,sphhar,atoms)
133 0 : tmp_results = results
134 :
135 0 : DO jspin = 1,input%jspins
136 0 : CALL cdnvalJob%init(fmpi,input,kpts,noco,results,jspin)
137 :
138 :
139 : ! replace brillouin weights with auxillary weights
140 0 : CALL calc_EnergyDen_auxillary_weights(eig_id, kpts, jspin, cdnvalJob%weights)
141 :
142 : CALL cdnval(eig_id, fmpi, kpts, jspin, noco,nococonv, input, banddos, cell, atoms, &
143 : enpara, stars, vacuum, sphhar, sym, vTot, cdnvalJob, &
144 0 : EnergyDen, regCharges, dos, vacdos,tmp_results, moments, gfinp, hub1inp)
145 : ENDDO
146 :
147 0 : END SUBROUTINE calc_EnergyDen
148 :
149 0 : SUBROUTINE calc_EnergyDen_auxillary_weights(eig_id, kpts, jspin, f_ik)
150 : USE m_types_kpts
151 : USE m_eig66_io
152 : IMPLICIT NONE
153 : ! calculates new (auxillary-)weights as
154 : ! f_iks = w_iks * E_iks
155 : !, where f_iks are the new (auxillary-)weights
156 : ! w_iks are the weights used in brillouin zone integration
157 : ! E_iks are the eigen energies
158 :
159 : INTEGER, INTENT(in) :: eig_id
160 : INTEGER, INTENT(in) :: jspin
161 : TYPE(t_kpts), INTENT(in) :: kpts
162 : REAL, INTENT(inout) :: f_ik(:,:) ! f_ik(band_idx, kpt_idx)
163 :
164 : ! local vars
165 0 : REAL :: e_i(SIZE(f_ik,dim=1))
166 : INTEGER :: ikpt
167 :
168 0 : DO ikpt = 1,kpts%nkpt
169 0 : CALL read_eig(eig_id,ikpt,jspin, eig=e_i)
170 0 : f_ik(:,ikpt) = f_ik(:,ikpt) * e_i
171 : ENDDO
172 0 : END SUBROUTINE calc_EnergyDen_auxillary_weights
173 :
174 0 : subroutine set_zPrime(dim_idx, zMat, kpt, lapw, cell, zPrime)
175 : USE m_types
176 : USE m_constants
177 : implicit none
178 : INTEGER, intent(in) :: dim_idx
179 : TYPE (t_mat), intent(in) :: zMat
180 : REAL, intent(in) :: kpt(3)
181 : TYPE(t_lapw), intent(in) :: lapw
182 : TYPE(t_cell), intent(in) :: cell
183 : TYPE (t_mat) :: zPrime
184 :
185 : REAL :: k_plus_g(3), fac
186 : INTEGER :: basis_idx
187 :
188 0 : call zPrime%free()
189 0 : call zPrime%init(zMat)
190 :
191 0 : do basis_idx = 1,size(lapw%gvec,dim=2)
192 0 : k_plus_g = kpt + lapw%gvec(:,basis_idx,1)
193 0 : k_plus_g = internal_to_rez(cell, k_plus_g)
194 :
195 0 : fac = k_plus_g(dim_idx)
196 0 : if(zPrime%l_real) then
197 0 : zPrime%data_r(basis_idx,:) = fac * zMat%data_r(basis_idx,:)
198 : else
199 0 : zPrime%data_c(basis_idx,:) = ImagUnit * fac * zMat%data_c(basis_idx,:)
200 : endif
201 : enddo
202 0 : end subroutine set_zPrime
203 :
204 0 : function internal_to_rez(cell, vec) result(res)
205 : use m_types
206 : implicit none
207 : type(t_cell), intent(in) :: cell
208 : real, intent(in) :: vec(3)
209 : real :: res(3)
210 :
211 0 : res = matmul(transpose(cell%bmat), vec)
212 0 : end function internal_to_rez
213 :
214 0 : subroutine undo_vgen_finalize(vtot, atoms, noco, stars)
215 : use m_types
216 : use m_constants
217 : use m_judft
218 : implicit none
219 : TYPE(t_potden), intent(inout) :: vtot
220 : type(t_atoms), intent(in) :: atoms
221 : type(t_noco), intent(in) :: noco
222 : type(t_stars), intent(in) :: stars
223 :
224 : integer :: js, n, st
225 :
226 0 : do js = 1,size(vtot%mt,4)
227 0 : do n = 1,atoms%ntype
228 : vTot%mt(:atoms%jri(n),0,n,js) = vtot%mt(:atoms%jri(n),0,n,js) &
229 0 : / (atoms%rmsh(:atoms%jri(n),n) / sfp_const )
230 : enddo
231 : enddo
232 :
233 0 : if(.not. noco%l_noco) then
234 0 : do js=1,size(vtot%pw_w,2)
235 0 : do st=1,stars%ng3
236 0 : vTot%pw_w(st,js) = vTot%pw_w(st,js) * stars%nstr(st)
237 : enddo
238 : enddo
239 : else
240 0 : call juDFT_error("undo vgen_finalize not implemented for noco")
241 : endif
242 0 : end subroutine undo_vgen_finalize
243 :
244 674 : subroutine set_kinED(fmpi, sphhar, atoms, sym, xcpot, &
245 : input, noco, stars, vacuum ,cell, den, EnergyDen, vTot,kinED)
246 : use m_types
247 : use m_cdn_io
248 : implicit none
249 : TYPE(t_mpi),INTENT(IN) :: fmpi
250 : TYPE(t_sphhar),INTENT(IN) :: sphhar
251 : TYPE(t_atoms),INTENT(IN) :: atoms
252 : TYPE(t_sym), INTENT(IN) :: sym
253 : CLASS(t_xcpot),INTENT(IN) :: xcpot
254 : TYPE(t_input),INTENT(IN) :: input
255 : TYPE(t_noco),INTENT(IN) :: noco
256 : TYPE(t_stars),INTENT(IN) :: stars
257 : TYPE(t_vacuum),INTENT(IN) :: vacuum
258 :
259 : TYPE(t_cell),INTENT(IN) :: cell
260 : TYPE(t_potden),INTENT(IN) :: den, EnergyDen, vTot
261 : TYPE(t_kinED),INTENT(OUT) :: kinED
262 :
263 674 : TYPE(t_potden) :: vTot_corrected
264 :
265 : LOGICAL :: perform_MetaGGA
266 674 : TYPE(t_potden) :: core_den, val_den
267 : real :: rdum, tempDistance
268 : logical :: ldum
269 : perform_MetaGGA = ALLOCATED(EnergyDen%mt) &
270 674 : .AND. (xcpot%exc_is_MetaGGA() .or. xcpot%vx_is_MetaGGA())
271 674 : if(.not.perform_MetaGGA) return
272 : #ifdef CPP_LIBXC
273 0 : core_den=den
274 0 : CALL core_den%resetPotDen
275 : call readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,&
276 : CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
277 0 : 0,rdum,tempDistance,ldum,core_den,inFilename='cdnc')
278 0 : call val_den%subPotDen(den,core_den)
279 :
280 0 : call vTot_corrected%copyPotDen(vTot)
281 0 : call undo_vgen_finalize(vTot_corrected, atoms, noco, stars)
282 :
283 0 : call set_kinED_is(xcpot, input, noco, stars, sym, cell, den, EnergyDen, vTot_corrected,kinED)
284 : call set_kinED_mt(fmpi, sphhar, atoms, sym, noco,core_den, val_den, &
285 0 : xcpot, EnergyDen, input, vTot_corrected,kinED)
286 : #endif
287 674 : end subroutine set_kinED
288 : #ifdef CPP_LIBXC
289 0 : subroutine set_kinED_is(xcpot, input, noco, stars, sym, cell, den, EnergyDen, vTot,kinED)
290 : use m_types
291 : use m_pw_tofrom_grid
292 : implicit none
293 : CLASS(t_xcpot),INTENT(IN) :: xcpot
294 : TYPE(t_input),INTENT(IN) :: input
295 : TYPE(t_noco),INTENT(IN) :: noco
296 : TYPE(t_stars),INTENT(IN) :: stars
297 : TYPE(t_sym), INTENT(IN) :: sym
298 : TYPE(t_cell),INTENT(IN) :: cell
299 : TYPE(t_potden),INTENT(IN) :: den, EnergyDen, vTot
300 : TYPE(t_kinED),INTENT(INOUT) :: kinED
301 :
302 : !local arrays
303 0 : REAL, ALLOCATABLE :: den_rs(:,:), ED_rs(:,:), vTot_rs(:,:)
304 0 : TYPE(t_gradients) :: tmp_grad
305 :
306 0 : CALL init_pw_grid(stars,sym,cell,xcpot)
307 :
308 : CALL pw_to_grid(xcpot%needs_grad(), input%jspins, noco%l_noco, stars, &
309 0 : cell, EnergyDen%pw, tmp_grad, xcpot, ED_rs)
310 : CALL pw_to_grid(xcpot%needs_grad(), input%jspins, noco%l_noco, stars, &
311 0 : cell, vTot%pw, tmp_grad, xcpot, vTot_rs)
312 : CALL pw_to_grid(xcpot%needs_grad(), input%jspins, noco%l_noco, stars, &
313 0 : cell, den%pw, tmp_grad, xcpot, den_rs)
314 :
315 0 : CALL finish_pw_grid()
316 :
317 0 : call calc_kinEnergyDen_pw(ED_rs, vTot_rs, den_rs, kinED%is)
318 : !xcpot%kinED%is = ED_RS - vTot_RS * den_RS
319 0 : kinED%set = .True.
320 0 : end subroutine set_kinED_is
321 :
322 0 : subroutine set_kinED_mt(fmpi, sphhar, atoms, sym, noco,core_den, val_den, &
323 : xcpot, EnergyDen, input, vTot,kinED)
324 : use m_types
325 : use m_mt_tofrom_grid
326 : implicit none
327 : TYPE(t_mpi),INTENT(IN) :: fmpi
328 : TYPE(t_sphhar),INTENT(IN) :: sphhar
329 : TYPE(t_atoms),INTENT(IN) :: atoms
330 : TYPE(t_sym), INTENT(IN) :: sym
331 : TYPE(t_noco), INTENT(IN) :: noco
332 : TYPE(t_potden),INTENT(IN) :: core_den, val_den, EnergyDen, vTot
333 : CLASS(t_xcpot),INTENT(IN) :: xcpot
334 : TYPE(t_input),INTENT(IN) :: input
335 : TYPE(t_kinED),INTENT(INOUT) :: kinED
336 : INTEGER :: jr, loc_n, n, n_start, n_stride, cnt
337 0 : REAL,ALLOCATABLE :: vTot_mt(:,:,:), ED_rs(:,:), vTot_rs(:,:), vTot0_rs(:,:),&
338 : core_den_rs(:,:), val_den_rs(:,:)
339 0 : TYPE(t_gradients) :: tmp_grad
340 0 : TYPE(t_sphhar) :: tmp_sphhar
341 :
342 : #ifdef CPP_MPI
343 0 : n_start=fmpi%irank+1
344 0 : n_stride=fmpi%isize
345 : #else
346 : n_start=1
347 : n_stride=1
348 : #endif
349 0 : CALL init_mt_grid(input%jspins,atoms,sphhar,xcpot%needs_grad(),sym)
350 0 : loc_n = 0
351 0 : allocate(ED_rs(atoms%nsp()*atoms%jmtd, input%jspins))
352 0 : allocate(vTot_rs, mold=ED_rs)
353 0 : allocate(vTot0_rs, mold=ED_rs)
354 0 : allocate(core_den_rs, mold=ED_rs)
355 0 : allocate(val_den_rs, mold=ED_rs)
356 :
357 : call kinED%alloc_mt(atoms%nsp()*atoms%jmtd, input%jspins, &
358 0 : n_start, atoms%ntype, n_stride)
359 0 : loc_n = 0
360 0 : do n = n_start,atoms%ntype,n_stride
361 0 : loc_n = loc_n + 1
362 :
363 0 : if(.not. allocated(vTot_mt)) then
364 : allocate(vTot_mt(lbound(vTot%mt, dim=1):ubound(vTot%mt, dim=1),&
365 : lbound(vTot%mt, dim=2):ubound(vTot%mt, dim=2),&
366 0 : lbound(vTot%mt, dim=4):ubound(vTot%mt, dim=4)))
367 : endif
368 :
369 0 : do jr=1,atoms%jri(n)
370 0 : vTot_mt(jr,0:,:) = vTot%mt(jr,0:,n,:) * atoms%rmsh(jr,n)**2
371 : enddo
372 : CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms,sym, sphhar,.TRUE., EnergyDen%mt(:, 0:, n, :), &
373 0 : n, noco, tmp_grad, ED_rs)
374 : CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,sphhar,.TRUE., vTot_mt(:,0:,:), &
375 0 : n, noco,tmp_grad, vTot_rs)
376 :
377 0 : tmp_sphhar%nlhd = sphhar%nlhd
378 0 : tmp_sphhar%nlh = [(0, cnt=1,size(sphhar%nlh))]
379 :
380 : CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,tmp_sphhar,.TRUE., vTot_mt(:,0:0,:), &
381 0 : n, noco, tmp_grad, vTot0_rs)
382 : CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,sphhar,.TRUE., &
383 0 : core_den%mt(:,0:,n,:), n,noco, tmp_grad, core_den_rs)
384 : CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,sphhar,.TRUE., &
385 0 : val_den%mt(:,0:,n,:), n,noco, tmp_grad, val_den_rs)
386 :
387 : call calc_kinEnergyDen_mt(ED_RS, vTot_rs, vTot0_rs, core_den_rs, val_den_rs, &
388 0 : kinED%mt(:,:,loc_n))
389 : !xcpot%kinED%mt(:,:,loc_n) = ED_RS - (vTot0_rs * core_den_rs + vTot_rs * val_den_rs)
390 : enddo
391 0 : kinED%set = .True.
392 0 : CALL finish_mt_grid()
393 0 : end subroutine set_kinED_mt
394 : #endif
395 0 : END MODULE m_metagga
|
