Line data Source code
1 : !!****m* ABINIT/m_precon
2 : !! NAME
3 : !! m_precon
4 : !!
5 : !! FUNCTION
6 : !! Object used for chi0-based preconditioning of the SCF.
7 : !!
8 : !! SOURCE
9 :
10 : #if defined HAVE_CONFIG_H
11 : #include "config.h"
12 : #endif
13 :
14 : #include "abi_common.h"
15 :
16 : module m_precon
17 :
18 : use iso_c_binding
19 : use defs_abitypes, only : MPI_type
20 : use defs_basis
21 : use m_dtset
22 : use m_dtfil
23 : use m_xmpi
24 :
25 : use defs_datatypes, only : pseudopotential_type
26 : use defs_wvltypes
27 : use m_atomdata, only : atom_length
28 : use m_bandfft_kpt, only : bandfft_kpt, bandfft_kpt_get_ikpt, bandfft_kpt_set_ikpt
29 : use m_cgprj, only : ctocprj
30 : use m_cgtools
31 : use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxc_noncoll
32 : use m_fft, only : fourdp, fourwf, fftpac, zerosym
33 : use m_fftcore, only : sphereboundary
34 : use m_fourier_interpol, only : transgrid
35 : use m_iterative_solvers, only : cg_linear_solver, gmres_linear_solver
36 : use m_kg, only : ph1d3d
37 : use m_mkrho
38 : use m_mpinfo, only : proc_distrb_cycle, proc_distrb_band
39 : use m_occ, only : getnel
40 : use m_paw_dmft
41 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free
42 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_mpi_allgather, pawcprj_free
43 : use m_pawang, only : pawang_type
44 : use m_pawfgr, only : pawfgr_type
45 : use m_pawtab, only : pawtab_type
46 : use m_pawfgrtab, only : pawfgrtab_type
47 : use m_paw_finegrid, only : pawgylmg
48 : use m_paw_occupancies, only : pawmkrhoij
49 : use m_paw_mkrho, only : pawmkrho
50 : use m_paw_nhat, only : pawsushat
51 : use m_prep_kgb, only : prep_getghc, prep_index_wavef_bandpp, prep_fourwf
52 : use m_spacepar, only : symrhg
53 :
54 : implicit none
55 : private
56 :
57 : type, public :: precon_object
58 : integer :: iprcel
59 : real(dp) :: dielng, diemix
60 : !Geometry :
61 : real(dp) :: gprimd(3, 3), rprimd(3, 3), gmet(3, 3), rmet(3, 3)
62 : real(dp) :: ucvol, dvol
63 : !PAW :
64 : type(pseudopotential_type), pointer :: psps
65 : integer :: unpaw
66 : integer, pointer :: dimcprj(:), mcprj, usecprj
67 : type(pawcprj_type), pointer :: cprj(:, :)
68 : type(pawang_type), pointer :: pawang
69 : type(pawfgr_type), pointer :: pawfgr
70 : type(pawfgrtab_type), pointer :: pawfgrtab(:)
71 : type(pawtab_type), pointer :: pawtab(:)
72 : real(dp), pointer :: ylm(:, :)
73 : real(dp), pointer :: ylmgr(:, :, :)
74 : !To compute (weighted) densities and other quantities :
75 : real(dp), pointer :: fermie
76 : real(dp), pointer :: cg(:, :), eigen(:), occ(:), ph1d(:, :), phnons(:, :, :)
77 : integer, pointer :: kg(:, :), npwarr(:), irrzon(:, :, :)
78 : integer, pointer :: atindx(:), atindx1(:), nattyp(:)
79 : integer, pointer :: symrec(:, :, :), indsym(:, :, :)
80 : real(dp), pointer :: xred(:, :)
81 : !For ffts :
82 : integer :: nfftprc ! Number of fft grid points for preconditioned quantities (densities and/or potentials).
83 : integer :: ngfftprc(18) ! All needed information about the 3D FFT for preconditioned quantities.
84 : !For Kxc :
85 : integer :: nkxc
86 : real(dp), pointer :: kxc(:, :)
87 : real(dp), pointer :: rhor(:, :)
88 : real(dp), pointer :: vxc(:, :)
89 :
90 : !Logical variables :
91 : logical :: use_precon
92 : logical :: use_ldos
93 : logical :: use_paw_rhoij
94 : logical :: use_dos
95 : logical :: use_kxc
96 : logical :: use_ridgereg
97 : logical :: use_indices_arrays
98 : logical :: use_precomputed_rhoi
99 : logical :: use_precomputed_psii
100 :
101 : !For LDOS preconditioner :
102 : real(dp) :: tdos
103 : real(dp), allocatable :: ldos(:, :)
104 : real(dp), allocatable :: dos(:)
105 :
106 : !Usefull
107 : integer, allocatable :: cg_indices(:, :, :, :)
108 : integer, allocatable :: kg_indices(:, :)
109 :
110 : !Array to store precomputed ffts of cg :
111 : real(dp), allocatable :: precomputed_rhoi(:, :, :)
112 : integer, allocatable :: precomputed_rhoi_indices(:, :, :)
113 : real(dp), allocatable :: precomputed_psii(:, :, :, :)
114 : integer, allocatable :: precomputed_psii_indices(:, :, :)
115 :
116 : !Preconditioner parameters
117 : integer :: precon_verbose
118 : !Linear solver parameters
119 : integer :: linsolve_maxiter
120 : real(dp) :: linsolve_rtol, ridge_param
121 : !chi0_diag/quasidiag parameters
122 : real(dp) :: deigvals_tol_fp
123 : real(dp) :: precon_tsmear
124 :
125 : contains
126 : procedure :: init => precon_init ! Initialize the precon_object.
127 : procedure :: init_kxc => precon_init_kxc ! Initialize kxc in the precon_object.
128 : procedure :: update => precon_update ! Update the precon_object according to iprcel.
129 : procedure :: free => precon_free ! Dealocate arrays that are allocated in precon_init.
130 :
131 : procedure :: apply_dielmat => apply_dielmat ! Apply the dielectric matrix to an input vector.
132 : procedure :: apply_adjdielmat => apply_adjdielmat ! Apply the adjoint dielectric matrix to an input vector.
133 : procedure :: apply_precon => apply_precon ! Apply the preconditioner to an input vector.
134 :
135 : end type precon_object
136 :
137 : contains
138 :
139 : ! TODO :
140 : ! - Debug apply_chi0_quasidiag
141 : ! - Implement non-collinear magnetism with band paral (Hyprid preconditioner)
142 : ! - Non-collinear magnetism: We assume nspinor=2 => nspden=4 and nsppol=2 => nspden=2, but we could have nspden=1 in both cases.
143 : ! For now, it is enforced in chkinp.
144 : ! - Compute the dos in an efficient way.
145 :
146 : !****f* m_precon/precon_init
147 : !! NAME
148 : !! precon_init
149 : !!
150 : !! FUNCTION
151 : !! Initialize the precon_object.
152 : !!
153 : !! INPUTS
154 : !! dtset = all input variables for this dataset
155 : !! atindx = index table for atoms (see gstate.f)
156 : !! atindx1 = index table for atoms, inverse of atindx (see gstate.f)
157 : !! cg = wf in G space
158 : !! cprj =
159 : !! dimcprj = dimension of the cprj array
160 : !! eigen = array of eigenvalues
161 : !! fermie = fermi energie
162 : !! gprimd = dimensional reciprocal space primitive translations
163 : !! irrzon = irreducible zone data
164 : !! kg = reduced planewave coordinates
165 : !! nattyp = number of atoms of each type.
166 : !! nfftmix = number of planewaves in the mixing/preconditioning grid
167 : !! ngfftmix =
168 : !! npwarr = number of planewaves and boundary planewaves at each k
169 : !! pawang = paw angular mesh and related data
170 : !! pawfgr =
171 : !! pawfgrtab =
172 : !! pawtab =
173 : !! phnons = nonsymmorphic translation phases
174 : !! psps = pseudopotential data
175 : !! rprimd = dimensional real space primitive translations
176 : !! ucvol = unit cell volume
177 : !! xred = reduced dimensionless atomic coordinates
178 : !!
179 : !! SOURCE
180 20754 : subroutine precon_init(this, dtset, atindx, atindx1, cg, cprj, dimcprj, dtfil, eigen, fermie, &
181 20754 : & gmet, gprimd, indsym, irrzon, kg, mcprj, nattyp, nfftmix, ngfftmix, npwarr, occ, pawang, pawfgr, pawfgrtab, &
182 20754 : & pawtab, ph1d, phnons, psps, rhor, rmet, rprimd, symrec, ucvol, usecprj, vxc, xred, ylm)
183 :
184 : !Arguments ------------------------------------
185 : class(precon_object), intent(out) :: this
186 : !scalars
187 : type(dataset_type),intent(in) :: dtset
188 : real(dp), intent(in) :: ucvol
189 : real(dp), intent(in), target :: fermie
190 : integer, intent(in) :: nfftmix
191 : integer, intent(in), target :: mcprj
192 : integer, intent(in), target :: usecprj
193 : !arrays
194 : real(dp), intent(in) :: gprimd(:, :), rprimd(:, :), gmet(:, :), rmet(:, :)
195 : integer, intent(in) :: ngfftmix(:)
196 : integer, intent(in), target :: irrzon(:, :, :), kg(:, :), npwarr(:)
197 : integer, intent(in), target :: atindx(:), atindx1(:), nattyp(:)
198 : integer, intent(in), target :: symrec(:, :, :), indsym(:, :, :)
199 : real(dp), intent(in), target :: cg(:, :), eigen(:), occ(:), phnons(:, :, :), ph1d(:, :)
200 : real(dp), intent(in), target :: rhor(:, :), vxc(:, :)
201 : real(dp), intent(in), target :: xred(:, :)
202 : type(datafiles_type),intent(in) :: dtfil
203 : type(pseudopotential_type), intent(in), target :: psps
204 : integer, intent(in), target :: dimcprj(:)
205 : type(pawcprj_type), intent(in), target :: cprj(:, :)
206 : type(pawang_type), intent(in), target :: pawang
207 : type(pawfgr_type), intent(in), target :: pawfgr
208 : type(pawfgrtab_type), intent(in), target :: pawfgrtab(:)
209 : type(pawtab_type), intent(in), target :: pawtab(:)
210 : real(dp), intent(in), target :: ylm(:, :)
211 :
212 : ! *************************************************************************
213 :
214 6918 : this%iprcel = dtset%iprcel
215 6918 : this%use_precon = .false.
216 6918 : this%use_kxc = .false.
217 :
218 6918 : if (this%iprcel >= 200 .and. this%iprcel < 300) then
219 :
220 2 : this%use_precon = .true.
221 :
222 : !Logical variables that describe the preconditioner :
223 : ! 200 -> LDOS with RPA
224 : ! 201 -> Kerker with DOS+ and DOS-
225 : ! 202 -> Hybrid (LDOS+diag)
226 : ! 203 -> Hybrid for antiferro (LDOS+quasidiag) - WIP - Not documented
227 : ! 210 -> Kerker with dielng - Not documented
228 : ! 211 -> Kerker with DOS - Not documented
229 : ! 212 -> LDOS without RPA - Not documented
230 : ! 299 -> No preconditioning - Not documented
231 :
232 : ! this%use_ldos = .true. activates the computation of the ldos.
233 2 : this%use_ldos = .false.
234 2 : if (this%iprcel == 200) this%use_ldos = .true.
235 2 : if (this%iprcel == 212) this%use_ldos = .true.
236 2 : if (this%iprcel == 202) this%use_ldos = .true.
237 2 : if (this%iprcel == 211) this%use_ldos = .true.
238 2 : if (this%iprcel == 201) this%use_ldos = .true.
239 :
240 2 : this%use_paw_rhoij = .false.
241 :
242 : ! If this%use_dos = .true. we will use this%dos.
243 2 : this%use_dos = .false.
244 2 : if (this%iprcel == 211) this%use_dos = .true.
245 2 : if (this%iprcel == 201) this%use_dos = .true.
246 :
247 : ! this%use_kxc = .true. activates the use of the exchange and correlation kernel.
248 : ! If this%use_kxc = .false. the RPA will be used.
249 : this%use_kxc = .false.
250 2 : if (this%iprcel == 212) this%use_kxc = .true.
251 2 : if (this%iprcel == 202) this%use_kxc = .true.
252 :
253 : ! this%use_ridgereg = .true. activates the use of an adapted linear solver.
254 2 : this%use_ridgereg = .false.
255 :
256 : ! this%use_indices_arrays = .true. indicates that we will use the arrays this%cg_indices and this%kg_indices.
257 2 : this%use_indices_arrays = .false.
258 2 : if (this%iprcel == 202) this%use_indices_arrays = .true.
259 2 : if (this%iprcel == 203) this%use_indices_arrays = .true.
260 :
261 2 : this%use_precomputed_rhoi = .false.
262 2 : if (this%iprcel == 202 .and. dtset%precon_in_memory==1) this%use_precomputed_rhoi = .true.
263 2 : if (this%iprcel == 202 .and. dtset%precon_in_memory==0 .and. dtset%npband>1) then
264 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): With band parallelization precon_in_memory must be 1.")
265 : end if
266 :
267 2 : this%use_precomputed_psii = .false.
268 2 : if (this%iprcel == 203 .and. dtset%precon_in_memory==1) this%use_precomputed_psii = .true.
269 :
270 : ! Other than here, iprcel is only used in apply_chi0, apply_dielmat and apply_adjdielmat.
271 :
272 : !Constant data from dtset
273 2 : this%dielng = dtset%dielng
274 2 : this%diemix = dtset%diemix
275 2 : this%nfftprc = nfftmix ! FFT grid for preconditioned densities and/or potentials :
276 38 : this%ngfftprc = ngfftmix ! same grid as the one used for mixing.
277 : !Other constants
278 2 : this%dvol = ucvol/this%nfftprc ! factor for integrals in real space (on the preconditioning FFT grid) : sum(f) * dvol ~ integral f
279 26 : this%gprimd = gprimd
280 26 : this%rprimd = rprimd
281 26 : this%gmet = gmet
282 26 : this%rmet = rmet
283 2 : this%ucvol = ucvol
284 : !Pointers
285 2 : this%atindx => atindx
286 2 : this%atindx1 => atindx1
287 2 : this%cg => cg
288 2 : this%eigen => eigen
289 2 : this%fermie => fermie
290 2 : this%indsym => indsym
291 2 : this%irrzon => irrzon
292 2 : this%kg => kg
293 2 : this%nattyp => nattyp
294 2 : this%npwarr => npwarr
295 2 : this%occ => occ
296 2 : this%ph1d => ph1d
297 2 : this%phnons => phnons
298 2 : this%psps => psps
299 2 : this%rhor => rhor
300 2 : this%symrec => symrec
301 2 : this%vxc => vxc
302 2 : this%xred => xred
303 :
304 : !PAW :
305 2 : if (psps%usepaw==1) then
306 0 : this%unpaw = dtfil%unpaw
307 0 : this%cprj => cprj
308 0 : this%usecprj => usecprj
309 0 : this%dimcprj => dimcprj
310 0 : this%mcprj => mcprj
311 0 : this%pawang => pawang
312 0 : this%pawfgr => pawfgr
313 0 : this%pawfgrtab => pawfgrtab
314 0 : this%pawtab => pawtab
315 0 : this%ylm => ylm
316 : end if
317 :
318 : !Initializing LDOS specific variables
319 2 : if (this%use_ldos) then
320 : !Allocating the array containing ldos
321 8 : ABI_MALLOC(this%ldos, (this%nfftprc, dtset%nspden))
322 : end if
323 :
324 2 : if (this%use_dos) then
325 : !Allocating the array containing the dos
326 0 : ABI_MALLOC(this%dos, (dtset%nspden))
327 : end if
328 :
329 : !Initializing variables needed for Kxc
330 2 : if (this%use_kxc) then
331 : !Preparing the allocation of Kxc
332 1 : if (dtset%xclevel==1) then !LDA
333 0 : this%nkxc = 2*min(dtset%nspden,2)-1
334 1 : else if (dtset%xclevel==2)then !GGA+...
335 1 : if (dtset%nspden==1) then
336 0 : this%nkxc = 7
337 1 : else if (dtset%nspden==2) then
338 1 : this%nkxc = 19
339 : else
340 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): kxc not implemented for nspden > 2 (non-coll magn) in GGA")
341 : end if
342 : end if
343 : end if
344 :
345 : !Preconditioner parameters
346 2 : this%precon_verbose = dtset%precon_verbose
347 : !Linear solver parameters
348 2 : this%linsolve_maxiter = dtset%precon_ls_maxite
349 2 : this%linsolve_rtol = dtset%precon_ls_rtol
350 : ! For inversion of non positive definite (adjointe) dielectric matrix
351 2 : this%ridge_param = 0.01
352 :
353 : !chi0_diag/quasidiag parameters
354 2 : this%deigvals_tol_fp = tol10
355 2 : this%precon_tsmear = max(dtset%precon_tsmear, dtset%tsmear)
356 :
357 : !Usefull : indices mapping arrays
358 2 : if (this%use_indices_arrays)then
359 6 : ABI_MALLOC(this%cg_indices, (2*dtset%nspinor, dtset%mband, dtset%nkpt, dtset%nsppol))
360 3 : ABI_MALLOC(this%kg_indices, (2, dtset%nkpt))
361 : end if
362 2 : if (this%use_precomputed_rhoi) then
363 5 : ABI_MALLOC(this%precomputed_rhoi_indices, (dtset%mband, dtset%nkpt, dtset%nsppol))
364 : end if
365 2 : if (this%use_precomputed_psii) then
366 0 : ABI_MALLOC(this%precomputed_psii_indices, (dtset%mband, dtset%nkpt, dtset%nsppol))
367 : end if
368 :
369 : end if
370 :
371 6918 : end subroutine precon_init
372 :
373 : !****f* m_precon/precon_init_kxc
374 : !! NAME
375 : !! precon_init_kxc
376 : !!
377 : !! FUNCTION
378 : !! Initialize the exchange and correlation kernel (kxc) in the precon_object.
379 : !! This is needed because kxc needs to be initialized at a specific time.
380 : !!
381 : !! INPUTS
382 : !! kxc = exchange and correlation kernel.
383 : !!
384 : !! SOURCE
385 6918 : subroutine precon_init_kxc(this, kxc)
386 :
387 : !Arguments ------------------------------------
388 : class(precon_object), intent(inout) :: this
389 : real(dp), intent(in), target :: kxc(:, :)
390 :
391 : ! *************************************************************************
392 6918 : if (this%use_precon) then
393 2 : if (this%use_kxc) then
394 1 : this%kxc => kxc
395 : end if
396 : end if
397 :
398 6918 : end subroutine precon_init_kxc
399 :
400 : !****f* m_precon/precon_update
401 : !! NAME
402 : !! precon_update
403 : !!
404 : !! FUNCTION
405 : !! Update the precon_object :
406 : !! For preconditioners using the LDOS (iprcel = 200 or 212) :
407 : !! Compute the new ldos (local density of state) with current wavefunctions
408 : !! and the new tdos (total density of state = integral of ldos).
409 : !!
410 : !! INPUTS
411 : !! dtset = All input variables for this dataset.
412 : !! mpi_enreg = Information about MPI parallelization.
413 : !!
414 : !! SOURCE
415 18 : subroutine precon_update(this, dtset, mpi_enreg)
416 :
417 : !Arguments ------------------------------------
418 : class(precon_object), intent(inout) :: this
419 : type(dataset_type), intent(in) :: dtset
420 : type(MPI_type), intent(in) :: mpi_enreg
421 :
422 : !Local variables-------------------------------
423 : integer :: ispden
424 :
425 : ! *************************************************************************
426 :
427 18 : if (this%use_precon) then
428 :
429 : ! Indices in cg array
430 18 : if (this%use_indices_arrays) then
431 10 : call compute_cg_indices(dtset, mpi_enreg, this%npwarr, this%cg_indices)
432 10 : call compute_kg_indices(dtset, mpi_enreg, this%npwarr, this%kg_indices)
433 : end if
434 :
435 : !LDOS
436 18 : if (this%use_ldos) then
437 : !update ldos
438 18 : call compute_ldos(this, dtset, mpi_enreg, this%ldos)
439 : !update tdos
440 60768 : this%tdos = sum(this%ldos(:, 1)) * this%dvol
441 :
442 : ! TODO : Extremely Inefficient way to compute the DOS ...
443 18 : if (this%use_dos) then
444 0 : do ispden = 1, dtset%nspden
445 0 : this%dos(ispden) = sum(this%ldos(:, ispden)) * this%dvol
446 : end do
447 : end if
448 :
449 : end if
450 :
451 : !Diag/Quasidiag chi0
452 18 : if (this%use_precomputed_psii) then
453 0 : call precompute_psii(this, dtset, mpi_enreg)
454 : ! Will also precompute rhoi.
455 18 : elseif (this%use_precomputed_rhoi) then
456 10 : call precompute_rhoi(this, dtset, mpi_enreg)
457 : end if
458 :
459 : end if
460 18 : end subroutine precon_update
461 :
462 : !****f* m_precon/precon_free_update
463 : !! NAME
464 : !! precon_update
465 : !!
466 : !! FUNCTION
467 : !! Deallocates arrays that might have been allocated in 'precon_update'.
468 : !!
469 : !! INPUTS
470 : !! dtset = All input variables for this dataset.
471 : !! mpi_enreg = Information about MPI parallelization.
472 : !!
473 : !! SOURCE
474 18 : subroutine precon_free_update(this)
475 :
476 : !Arguments ------------------------------------
477 : class(precon_object), intent(inout) :: this
478 :
479 : ! *************************************************************************
480 18 : if (this%use_precon) then
481 :
482 18 : if (this%use_precomputed_rhoi) then
483 10 : ABI_FREE(this%precomputed_rhoi)
484 : end if
485 18 : if (this%use_precomputed_psii) then
486 0 : ABI_FREE(this%precomputed_psii)
487 : end if
488 :
489 : end if
490 18 : end subroutine precon_free_update
491 :
492 : !****f* m_precon/precon_free
493 : !! NAME
494 : !! precon_free
495 : !!
496 : !! FUNCTION
497 : !! Dealocate arrays that are allocated in precon_init.
498 : !!
499 : !! SOURCE
500 6918 : subroutine precon_free(this)
501 :
502 : !Arguments ------------------------------------
503 : class(precon_object), intent(inout) :: this
504 :
505 : ! *************************************************************************
506 6918 : if (this%use_precon) then
507 :
508 2 : if (this%use_indices_arrays) then
509 1 : ABI_FREE(this%cg_indices)
510 1 : ABI_FREE(this%kg_indices)
511 : end if
512 :
513 2 : if (this%use_precomputed_rhoi) then
514 1 : ABI_FREE(this%precomputed_rhoi_indices)
515 : end if
516 2 : if (this%use_precomputed_psii) then
517 0 : ABI_FREE(this%precomputed_psii_indices)
518 : end if
519 :
520 2 : if (this%use_ldos) then
521 : !Deallocating the array containing ldos and tdos
522 2 : ABI_FREE(this%ldos)
523 : end if
524 :
525 2 : if (this%use_dos) then
526 0 : ABI_FREE(this%dos)
527 : end if
528 :
529 : end if
530 6918 : end subroutine precon_free
531 :
532 : !****f* m_precon/compute_r
533 : !! NAME
534 : !! compute_r
535 : !!
536 : !! FUNCTION
537 : !! Computes the array of r-vectors (in REDUCED coordinates).
538 : !!
539 : !! INPUTS
540 : !! ngfft = All needed information about 3D FFT, see ~abinit/doc/variables/gstate/#ngfft.
541 : !!
542 : !! OUTPUTS
543 : !! r_vectors(3, :) = 3 coordinates of the r_vectors.
544 : !!
545 : !! -unused-
546 : !!
547 : !! SOURCE
548 : subroutine compute_r(ngfft, r_vectors)
549 :
550 : !Arguments ------------------------------------
551 : real(dp), intent(out) :: r_vectors(:, :)
552 : integer, intent(in) :: ngfft(:)
553 :
554 : !Local variables-------------------------------
555 : integer :: n1, n2, n3, i1, i2, i3, i_r
556 :
557 : ! *************************************************************************
558 :
559 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
560 : do i3=1,n3
561 : do i2=1,n2
562 : do i1=1,n1
563 : i_r = 1 + (i1-1) + (i2-1)*n1 + (i3-1)*n1*n2
564 : r_vectors(1, i_r) = real(i1-1)/n1
565 : r_vectors(2, i_r) = real(i2-1)/n2
566 : r_vectors(3, i_r) = real(i3-1)/n3
567 : end do
568 : end do
569 : end do
570 :
571 : end subroutine compute_r
572 :
573 : !****f* m_precon/get_r_vector
574 : !! NAME
575 : !! get_r_vector
576 : !!
577 : !! FUNCTION
578 : !! Get the vector r (in REDUCED coordinates) of index ifft.
579 : !!
580 : !! INPUTS
581 : !! ifft = Index of sought vector r.
582 : !! ngfft = All needed information about 3D FFT, see ~abinit/doc/variables/gstate/#ngfft.
583 : !!
584 : !! OUTPUT
585 : !! r(3) = sought vector r
586 : !!
587 : !! SOURCE
588 : function get_r_vector(ifft, ngfft) result(r)
589 :
590 : !Arguments ------------------------------------
591 : integer, intent(in) :: ifft
592 : integer, intent(in) :: ngfft(:)
593 :
594 : !Local variables-------------------------------
595 : integer :: n1, n2, n3, i1, i2, i3
596 :
597 : !Returned variable-------------------------------
598 : real(dp) :: r(3)
599 :
600 : ! *************************************************************************
601 :
602 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
603 : i1 = modulo((ifft-1), n1) + 1
604 : i2 = modulo((ifft-1)/n1, n2) + 1
605 : i3 = ((ifft-1)/n1)/n2 + 1
606 : r(1) = real(i1-1)/n1
607 : r(2) = real(i2-1)/n2
608 : r(3) = real(i3-1)/n3
609 :
610 : end function get_r_vector
611 :
612 : !****f* m_precon/get_g_vector
613 : !! NAME
614 : !! get_g_vector
615 : !!
616 : !! FUNCTION
617 : !! Get the vector g (in REDUCED coordinates) of index ifft
618 : !!
619 : !! INPUTS
620 : !! ifft = Index of sought vector g.
621 : !! ngfft = All needed information about 3D FFT, see ~abinit/doc/variables/gstate/#ngfft.
622 : !!
623 : !! OUTPUT
624 : !! g(3) = sought vector g
625 : !!
626 : !! SOURCE
627 303660 : function get_g_vector(ifft, ngfft) result(g)
628 :
629 : !Arguments ------------------------------------
630 : integer, intent(in) :: ifft
631 : integer, intent(in) :: ngfft(:)
632 :
633 : !Local variables-------------------------------
634 : integer :: n1, n2, n3, i1, i2, i3
635 :
636 : !Returned variable-------------------------------
637 : integer(dp) :: g(3)
638 :
639 : ! *************************************************************************
640 :
641 303660 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
642 : ! (ifft-1) = (i1-1) + (i2-1)*n1 + (i3-1)*n1*n2
643 303660 : i1 = modulo(ifft-1, n1) +1
644 303660 : i2 = modulo((ifft-1)/n1, n2) +1
645 303660 : i3 = ((ifft-1)/n1)/n2 +1
646 303660 : g(1) = i1 - (i1/(n1/2+2))*n1-1
647 303660 : g(2) = i2 - (i2/(n2/2+2))*n2-1
648 303660 : g(3) = i3 - (i3/(n3/2+2))*n3-1
649 :
650 303660 : end function get_g_vector
651 :
652 : !!***
653 : !!****f* ABINIT/to_pauli
654 : !! NAME
655 : !! to_pauli
656 : !!
657 : !! FUNCTION
658 : !! Basis change from the default spin-basis to the Pauli basis for potentials and densities
659 : !! in the direct (r) space.
660 : !!
661 : !! INPUT/OUTPUT
662 : !! opt = 0 : v is a potential
663 : !! 1 : v is a density
664 : !! v(nfft, nspden) = On input : Potential/density in the default spin-basis.
665 : !! On output : Potential/density in the Pauli basis.
666 : !!
667 : !! SOURCE
668 128 : subroutine to_pauli(opt, v)
669 : !Arguments ------------------------------------
670 : real(dp), intent(inout) :: v(:, :)
671 : integer :: opt
672 : !Local variables-------------------------------
673 : integer :: nspden
674 128 : real(dp), allocatable :: temp(:)
675 :
676 : ! *************************************************************************
677 128 : nspden = size(v, 2)
678 :
679 : !sigma_0, ... , sigma_3 are the Pauli matrices.
680 128 : if (opt == 0) then !v is a potential
681 55 : if (nspden == 2) then
682 : !On input v(:, 1) is the spin-up potential and v(:, 2) is the spin-down potential.
683 : !On output the entire potential is v(:, 1)*sigma_0 + v(:, 2)*sigma_3.
684 185680 : v(:, 1) = 0.5_dp*(v(:, 1) + v(:, 2))
685 185680 : v(:, 2) = v(:, 1) - v(:, 2)
686 0 : else if (nspden == 4) then
687 : ! v(:, 1) | v(:, 3) + i*v(:, 4)
688 : !On input the entire potential is --------------------|----------------------
689 : ! v(:, 3) - i*v(:, 4) | v(:, 2)
690 : ! (see dotprod_vn in m_cgtools)
691 : !On output the entire potential is
692 : ! v(:, 1)*sigma_0 + v(:, 2)*sigma_1 + v(:, 3)*sigma_2 + v(:, 4)*sigma_3
693 : ! v(:, 1) + v(:, 4) | v(:, 2) - i*v(:, 3)
694 : ! = --------------------|----------------------
695 : ! v(:, 2) + i*v(:, 3) | v(:, 1) - v(:, 4)
696 :
697 0 : ABI_MALLOC(temp, (size(v, 1)))
698 :
699 0 : v(:, 1) = 0.5_dp*(v(:, 1) + v(:, 2))
700 0 : temp = v(:, 4)
701 0 : v(:, 4) = v(:, 1) - v(:, 2) ! = 0.5_dp*(v(:, 1)_old - v(:, 2))
702 0 : v(:, 2) = v(:, 3)
703 0 : v(:, 3) = -temp
704 :
705 0 : ABI_FREE(temp)
706 : end if
707 :
708 73 : else if (opt == 1) then !v is a density
709 73 : if (nspden == 2) then
710 : !On input v(:, 1) is the total density and v(:, 2) is the spin-up density.
711 : !On output v(:, 1) is the total density and v(:, 2) is the spin density.
712 246448 : v(:, 2) = 2*v(:, 2) - v(:, 1)
713 : end if
714 : !If nspden=4, the density is already given in the Pauli basis.
715 : end if
716 128 : end subroutine to_pauli
717 :
718 : !!***
719 : !!****f* ABINIT/from_pauli
720 : !! NAME
721 : !! from_pauli
722 : !!
723 : !! FUNCTION
724 : !! Basis change from the Pauli basis to the Abinit default spin-basis for potentials and densities
725 : !! in the real space.
726 : !!
727 : !! INPUT/OUTPUT
728 : !! opt = 0 : v is a potential
729 : !! 1 : v is a density
730 : !! v(nfft, nspden) = On input : Potential/density in the Pauli basis
731 : !! On output : Potential/density in the default Abinit spin-basis.
732 : !!
733 : !! SOURCE
734 73 : subroutine from_pauli(opt, v)
735 : !Arguments ------------------------------------
736 : real(dp), intent(inout) :: v(:, :)
737 : integer :: opt
738 : !Local variables-------------------------------
739 : integer :: nspden
740 73 : real(dp), allocatable :: temp(:)
741 :
742 : ! *************************************************************************
743 73 : nspden = size(v, 2)
744 :
745 : !sigma_0, ... , sigma_3 are the Pauli matrices.
746 73 : if (opt == 0) then !v is a potential
747 0 : if (nspden == 2) then
748 : !On input the entire potential is v(:, 1)*sigma_0 + v(:, 2)*sigma_3.
749 : !On output v(:, 1) is the spin-up potential and v(:, 2) is the spin-down potential.
750 0 : v(:, 1) = v(:, 1) + v(:, 2)
751 0 : v(:, 2) = v(:, 1) - 2*v(:, 2)
752 0 : else if (nspden == 4) then
753 :
754 : !On input the entire potential is
755 : ! v(:, 1)*sigma_0 + v(:, 2)*sigma_1 + v(:, 3)*sigma_2 + v(:, 4)*sigma_3
756 : ! v(:, 1) + v(:, 4) | v(:, 2) - i*v(:, 3)
757 : ! = --------------------|----------------------
758 : ! v(:, 2) + i*v(:, 3) | v(:, 1) - v(:, 4)
759 : !
760 : ! v(:, 1) | v(:, 3) + i*v(:, 4)
761 : !On output the entire potential is --------------------|----------------------
762 : ! v(:, 3) - i*v(:, 4) | v(:, 2)
763 : ! (see dotprod_vn in m_cgtools)
764 :
765 0 : ABI_MALLOC(temp, (size(v, 1)))
766 :
767 0 : temp = v(:, 2)
768 0 : v(:, 2) = v(:, 1) - v(:, 4) ! = v(:, 1)_old - v(:, 4)_old
769 0 : v(:, 1) = v(:, 1) + v(:, 4) ! = v(:, 1)_old + v(:, 4)_old
770 0 : v(:, 4) = -v(:, 3) ! = -v(:, 3)_old
771 0 : v(:, 3) = temp ! = v(:, 2)_old
772 :
773 0 : ABI_FREE(temp)
774 :
775 : end if
776 :
777 73 : else if (opt == 1) then !v is a density
778 73 : if (nspden == 2) then
779 : !On input v(:, 1) is the total density and v(:, 2) is the spin density.
780 : !On output v(:, 1) is the total density and v(:, 2) is the spin-up density.
781 246448 : v(:, 2) = 0.5_dp*(v(:, 1) + v(:, 2))
782 : end if
783 : !If nspden=4, the density is already given in the Pauli basis.
784 : end if
785 73 : end subroutine from_pauli
786 :
787 : !****f* m_precon/apply_vc
788 : !! NAME
789 : !! apply_vc
790 : !!
791 : !! FUNCTION
792 : !! Apply the Coulomb kernel vc to a vector (in place) in the Pauli basis.
793 : !!
794 : !! INPUTS
795 : !!
796 : !! SIDE EFFECTS
797 : !! vec_r (nfftprc, nspden) = Vector (in direct space) to which the Coulomb kernel vc is applied (in place).
798 : !! When nspden > 1 vec_r is in the Pauli basis.
799 : !!
800 : !! SOURCE
801 90 : subroutine apply_vc(this, dtset, mpi_enreg, vec_r)
802 : !Arguments ------------------------------------
803 : class(precon_object), intent(in) :: this
804 : type(dataset_type),intent(in) :: dtset
805 : type(MPI_type), intent(in) :: mpi_enreg
806 : !arrays
807 : real(dp), intent(inout) :: vec_r(this%nfftprc, dtset%nspden)
808 :
809 : !Local variables-------------------------------
810 : integer :: ifft, ispden, cplex
811 : real(dp) :: g_cart_2
812 : integer :: n1, n2, n3
813 180 : real(dp) :: vec_g(2, this%nfftprc, 1)
814 :
815 : ! *************************************************************************
816 :
817 : !In the sigma_0, 1, 2, 3 (pauli) basis :
818 : ! The sigma_0 component of the density is multiplied by 4pi/G^2
819 : ! and the rest is 0.
820 :
821 90 : cplex = 1 ! vec is REAL
822 90 : n1=this%ngfftprc(1) ; n2=this%ngfftprc(2) ; n3=this%ngfftprc(3)
823 :
824 : ! FFT
825 90 : call fourdp(cplex, vec_g, vec_r(:, 1), -1, mpi_enreg, this%nfftprc, 1, this%ngfftprc, 0)
826 :
827 90 : ispden = 1
828 303750 : do ifft = 2, this%nfftprc
829 6073200 : g_cart_2 = norm2(two_pi * matmul(this%gprimd, get_g_vector(ifft, this%ngfftprc)))**2
830 911070 : vec_g(:, ifft, ispden) = (2*two_pi/g_cart_2) * vec_g(:, ifft, ispden)
831 : end do
832 :
833 : ! Set contribution of unbalanced components to zero.
834 90 : call zerosym(vec_g(:, :, 1), 2, n1, n2, n3)
835 :
836 : ! iFFT
837 90 : call fourdp(cplex, vec_g, vec_r(:, 1), 1, mpi_enreg, this%nfftprc, 1, this%ngfftprc, 0)
838 :
839 180 : do ispden = 2, dtset%nspden
840 303930 : vec_r(:, ispden) = 0
841 : end do
842 :
843 90 : end subroutine apply_vc
844 :
845 : !****f* m_precon/apply_kxc
846 : !! NAME
847 : !! apply_kxc
848 : !!
849 : !! FUNCTION
850 : !! Apply the exchange and correlation kernel Kxc to a vector (in real space).
851 : !!
852 : !! INPUTS
853 : !! dtset = All input variables for this dataset.
854 : !! mpi_enreg = Informations about MPI parallelization.
855 : !! vec_r (nfftprc, nspden) = Vector (in real space) to which the exchange and correlation kernel Kxc is applied.
856 : !! When nspden > 1 vec_r is in the default Abinit spin-basis.
857 : !!
858 : !! OUTPUTS
859 : !! Kxc_vec_r (nfftprc, nspden) = Resulting vector (in real space) containing the application of Kxc to vec_r.
860 : !! When nspden > 1 Kxc_vec_r is in the Pauli-basis.
861 : !!
862 : !! SOURCE
863 55 : subroutine apply_kxc(this, dtset, mpi_enreg, vec_r, Kxc_vec_r)
864 : !Arguments ------------------------------------
865 : class(precon_object), intent(in) :: this
866 : !scalars
867 : type(dataset_type),intent(in) :: dtset
868 : type(MPI_type),intent(in) :: mpi_enreg
869 : !arrays
870 : real(dp), intent(in) :: vec_r(this%nfftprc, dtset%nspden)
871 : real(dp), intent(inout) :: Kxc_vec_r(this%nfftprc, dtset%nspden)
872 :
873 : !Local variables-------------------------------
874 : !scalars
875 : integer :: cplex, n3xccc, nhatdim, nhat1dim, nhat1grdim, nkxc, option, optnc, usexcnhat
876 : logical :: non_magnetic_xc
877 : !arrays
878 55 : real(dp), allocatable :: vec_r_default(:, :)
879 55 : real(dp), allocatable :: nhat(:, :), nhat1(:, :), nhat1gr(:, :, :)
880 : real(dp) :: dummy_xccc3d1(0), qphon(3)
881 :
882 : ! *************************************************************************
883 :
884 55 : if (size(this%kxc, 1) /= this%nfftprc) then
885 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): size(kxc, 1) /= nfftprc")
886 : end if
887 :
888 : !Applying Kxc :
889 55 : cplex = 1 ! Input vector is real in real (direct) space.
890 55 : non_magnetic_xc = .false.
891 55 : nkxc = size(this%kxc, 2)
892 :
893 55 : usexcnhat = 0 !
894 55 : nhat1dim = 0 !
895 110 : ABI_MALLOC(nhat1, (cplex*this%nfftprc, dtset%nspden*nhat1dim)) ! PAW
896 55 : nhat1grdim = 0 !
897 165 : ABI_MALLOC(nhat1gr, (cplex*this%nfftprc, dtset%nspden, 3*nhat1grdim)) !
898 :
899 55 : option = 2 ! Treats only density change (no core_correction)
900 55 : n3xccc = 0 ! -> Core-correction set to 0.
901 55 : qphon = 0.0_dp ! phonon vector
902 :
903 55 : if (dtset%nspden==1) then
904 : call dfpt_mkvxc(cplex, dtset%ixc ,this%kxc, mpi_enreg, this%nfftprc, this%ngfftprc, nhat1, nhat1dim, &
905 : & nhat1gr, nhat1grdim, nkxc, non_magnetic_xc, dtset%nspden, n3xccc, option, &
906 0 : & qphon, vec_r, this%rprimd, usexcnhat, Kxc_vec_r, dummy_xccc3d1)
907 :
908 55 : else if (dtset%nspden==2) then
909 : ! Basis change to the default Abinit spin-basis for densities in case
910 220 : ABI_MALLOC(vec_r_default, ((this%nfftprc), dtset%nspden))
911 371470 : vec_r_default = vec_r
912 55 : call from_pauli(1, vec_r_default)
913 : call dfpt_mkvxc(cplex, dtset%ixc ,this%kxc, mpi_enreg, this%nfftprc, this%ngfftprc, nhat1, nhat1dim, &
914 : & nhat1gr, nhat1grdim, nkxc, non_magnetic_xc, dtset%nspden, n3xccc, option, &
915 55 : & qphon, vec_r_default, this%rprimd, usexcnhat, Kxc_vec_r, dummy_xccc3d1)
916 55 : ABI_FREE(vec_r_default)
917 :
918 0 : else if (dtset%nspden==4) then
919 : ! In non-collinear spin, the default Abinit spin-basis is the pauli basis. No basis change needed.
920 0 : nhatdim = 0
921 0 : ABI_MALLOC(nhat, (this%nfftprc, dtset%nspden*nhatdim)) !
922 0 : optnc = 1 ! Compute the whole 2x2 Vres matrix
923 : call dfpt_mkvxc_noncoll(cplex, dtset%ixc ,this%kxc, mpi_enreg, this%nfftprc, this%ngfftprc, nhat, nhatdim, &
924 : & nhat1, nhat1dim, nhat1gr, nhat1grdim, nkxc, non_magnetic_xc, dtset%nspden, &
925 : & n3xccc, optnc, option, qphon, this%rhor, vec_r, this%rprimd, usexcnhat, &
926 0 : & this%vxc, Kxc_vec_r, dummy_xccc3d1)
927 0 : ABI_FREE(nhat)
928 : end if
929 :
930 55 : ABI_FREE(nhat1)
931 55 : ABI_FREE(nhat1gr)
932 :
933 55 : call to_pauli(0, Kxc_vec_r)
934 :
935 55 : end subroutine apply_kxc
936 :
937 : !****f* m_precon/apply_kernel
938 : !! NAME
939 : !! apply_kernel
940 : !!
941 : !! FUNCTION
942 : !! Apply the kernel vc or (vc + Kxc) depending in iprcel to a vector (in place).
943 : !!
944 : !! INPUTS
945 : !! dtset = All input variables for this dataset.
946 : !! mpi_enreg = Informations about MPI parallelization.
947 : !!
948 : !! SIDE EFFECTS
949 : !! vec_r (nfftprc, nspden) = Vector (in direct space) to which the kernel is applied (in place).
950 : !! When nspden > 1 vec_r is in the default Abinit spin-basis.
951 : !!
952 : !! SOURCE
953 0 : subroutine apply_kernel(this, dtset, mpi_enreg, vec_r)
954 : !Arguments ------------------------------------
955 : class(precon_object), intent(inout) :: this
956 : !scalars
957 : type(dataset_type),intent(in) :: dtset
958 : type(MPI_type),intent(in) :: mpi_enreg
959 : !arrays
960 : real(dp), intent(inout) :: vec_r(this%nfftprc, dtset%nspden)
961 :
962 : !Local variables-------------------------------
963 0 : real(dp), allocatable :: Kxc_vec_r(:, :)
964 :
965 : ! *************************************************************************
966 :
967 : ! RPA : LDOS/Kerker model - only vc
968 0 : if (.not. this%use_kxc) then
969 :
970 0 : call apply_vc(this, dtset, mpi_enreg, vec_r) ! Apply vc in place
971 :
972 : ! No RPA : vc and Kxc
973 : else
974 :
975 : ! Apply Kxc
976 0 : ABI_MALLOC(Kxc_vec_r, (this%nfftprc, dtset%nspden))
977 0 : call apply_Kxc(this, dtset, mpi_enreg, vec_r, Kxc_vec_r)
978 :
979 : ! Apply vc in place
980 0 : call apply_vc(this, dtset, mpi_enreg, vec_r) ! Apply vc in place
981 :
982 : ! Add Kxc_vec_r to vec_r
983 0 : vec_r = vec_r + Kxc_vec_r
984 0 : ABI_FREE(Kxc_vec_r)
985 :
986 : end if
987 :
988 0 : end subroutine apply_kernel
989 :
990 : !****f* m_precon/derivative_occ
991 : !! NAME
992 : !! derivative_occ
993 : !!
994 : !! FUNCTION
995 : !! Compute the derivative of the occupation (f) of the band corresponding to eigenval
996 : !! with respect to the fermie temperature.
997 : !! f = integral_((eigenval - fermie)/tsmear)^infty delta(t) dt
998 : !! f' = -1/tsmear * delta((eigenval - fermie)/tsmear)
999 : !!
1000 : !! INPUTS
1001 : !! occopt = option for occupancies, determines delta
1002 : !! eigenval = eigenvalue
1003 : !! fermie = fermi energie
1004 : !! tsmear = smearing temperature
1005 : !!
1006 : !! OUTPUT
1007 : !! fprim = occupation derivative
1008 : !!
1009 : !! SOURCE
1010 11898 : function derivative_occ(occopt, eigenval, fermie, tsmear) result(fprim)
1011 :
1012 : !Arguments ------------------------------------
1013 : !scalars
1014 : real(dp), intent(in) :: eigenval, fermie, tsmear
1015 : integer, intent(in) :: occopt
1016 :
1017 : !Local variables-------------------------------
1018 : !scalars
1019 : real(dp) :: x, delta, a
1020 :
1021 : !Returned variable-------------------------------
1022 : real(dp) :: fprim
1023 :
1024 : ! *************************************************************************
1025 :
1026 11898 : x = (eigenval - fermie)/tsmear
1027 :
1028 11898 : if (occopt<=2) then
1029 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): Non-metallic occupation.")
1030 : else if (occopt==3) then
1031 : !Fermi-Dirac smearing
1032 11898 : delta = exp(-abs(x))/(1+exp(-abs(x)))**2 !To avoid overflow of exp.
1033 : else if (occopt==4) then
1034 : !Cold Smearing
1035 0 : a = -0.5634
1036 0 : delta = (1.5+x*(-1.5*a+x*(-1.0+a*x)))*exp(-x**2)/sqrt(pi)
1037 : else if (occopt==5) then
1038 : !Cold Smearing
1039 0 : a = -0.8165
1040 0 : delta = (1.5+x*(-1.5*a+x*(-1.0+a*x)))*exp(-x**2)/sqrt(pi)
1041 : else if (occopt==6) then
1042 : !Smering of Methfessel and Paxton
1043 0 : a = 0.0
1044 0 : delta = (1.5+x*(-1.5*a+x*(-1.0+a*x)))*exp(-x**2)/sqrt(pi)
1045 : else if (occopt==7) then
1046 : !Gaussian smearing
1047 0 : delta = exp(-x**2)/sqrt(pi)
1048 : else if (occopt==8) then
1049 : !Uniform smearing
1050 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): preconditioner needs a smooth smearing function.")
1051 : else if (occopt==9) then
1052 : !Fermi-Dirac occupation is enforced with two distinct quasi-Fermi levels
1053 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): preconditioner not implemented for this smearing function.")
1054 : end if
1055 :
1056 11898 : fprim = -1/tsmear * delta
1057 :
1058 11898 : end function derivative_occ
1059 :
1060 : !****f* m_precon/compute_weighted_density
1061 : !! NAME
1062 : !! compute_weighted_density
1063 : !!
1064 : !! FUNCTION
1065 : !! Wrapper for mkrho, symrhg and PAW :
1066 : !! Compute a density-like quantity where the occupation are replaced by some weights
1067 : !! w_rho = sum_i weight_i |psi_i|^2 .
1068 : !! w_rho that has the same size as the preconditioned density/potential.
1069 : !!
1070 : !! INPUTS
1071 : !! dtset = All input variables for this dataset.
1072 : !! mpi_enreg = Information about MPI parallelization.
1073 : !! weights = Weights that replace the occupations in the computation of density.
1074 : !!
1075 : !! OUTPUT
1076 : !! w_rhor = "Weighted density" in real space, in the Pauli basis.
1077 : !!
1078 : !! SOURCE
1079 18 : subroutine compute_weighted_density(this, dtset, mpi_enreg, weights, w_rhor)
1080 :
1081 : !Arguments ------------------------------------
1082 : !scalars
1083 : class(precon_object), intent(in) :: this
1084 : type(dataset_type), intent(in) :: dtset
1085 : type(MPI_type), intent(in) :: mpi_enreg
1086 : !arrays
1087 : real(dp), intent(in) :: weights(:)
1088 : real(dp), intent(out) :: w_rhor(this%nfftprc, dtset%nspden)
1089 :
1090 : !Local variables-------------------------------
1091 : !scalars
1092 54 : type(pawrhoij_type) :: pawrhoij(mpi_enreg%my_natom*this%psps%usepaw)
1093 : integer :: dummy_int, mband_cprj, my_nspinor, mcprj_tmp
1094 : integer :: mcg, cplex, cplex_rhoij
1095 : real(dp) :: compch_fft
1096 : integer :: optin, optout, optgrid
1097 : !arrays
1098 : real(dp) :: qphon(3)
1099 18 : real(dp), allocatable :: w_rhowfg(:, :), w_rhowfr(:, :)
1100 18 : type(pawcprj_type), allocatable :: cprj_tmp(:,:)
1101 18 : type(paw_dmft_type) :: dummy_paw_dmft
1102 : type(wvl_wf_type) :: dummy_wvl_wfs
1103 : type(wvl_denspot_type) :: dummy_wvl_den
1104 : real(dp) :: dummy_ylmgr(0, 0, 0)
1105 18 : real(dp), allocatable :: dummy_rhog(:, :), dummy_rhogf(:, :)
1106 :
1107 : ! *************************************************************************
1108 :
1109 18 : if (this%psps%usepaw==0) then
1110 72 : ABI_MALLOC(w_rhowfr, (dtset%nfft, dtset%nspden))
1111 54 : ABI_MALLOC(w_rhowfg, (2, dtset%nfft))
1112 : else
1113 0 : ABI_MALLOC(w_rhowfr, (this%pawfgr%nfftc, dtset%nspden))
1114 0 : ABI_MALLOC(w_rhowfg, (2, this%pawfgr%nfftc))
1115 : end if
1116 : !w_rhowfr = zero
1117 : !w_rhowfg = zero
1118 :
1119 : ! Compute the weighted density (w_rhor) using mkrho with weights in place of the occupations.
1120 54 : mcg = size(this%cg)
1121 18 : dummy_paw_dmft%use_dmft = 0
1122 18 : dummy_paw_dmft%use_sc_dmft = 0
1123 : call mkrho(this%cg, dtset, this%gprimd, this%irrzon, this%kg, mcg, mpi_enreg, this%npwarr, weights, &
1124 18 : & dummy_paw_dmft, this%phnons, w_rhowfg, w_rhowfr, this%rprimd, 0, this%ucvol, dummy_wvl_den, dummy_wvl_wfs, option=0, printout=.false.)
1125 : ! symrhg already called in mkrho
1126 :
1127 18 : if (this%psps%usepaw==0) then
1128 : ! In NC : the weighted density is directly w_rhowfr.
1129 18 : if (this%nfftprc == dtset%nfft) then
1130 121554 : w_rhor = w_rhowfr
1131 : else
1132 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): nfftprc /= nfft in Norm-conserving not implemented.")
1133 : end if
1134 : else
1135 : ! In PAW :
1136 : ! First option : Add rhoij terms to w_rhowfr and transfer to fine grid. -UNUSED-
1137 0 : if (this%use_paw_rhoij .and. this%nfftprc == this%pawfgr%nfft) then
1138 :
1139 : !Compute the rhoij equivalent for the weighted density.
1140 : ! Sum_{n,k} {weight(n,k)*<Cnk|p_i><p_j|Cnk>}.
1141 :
1142 0 : my_nspinor = max(1, dtset%nspinor/mpi_enreg%nproc_spinor)
1143 0 : mband_cprj = dtset%mband / mpi_enreg%nproc_band
1144 :
1145 : !Initialize pawrhoij
1146 0 : cplex_rhoij = 1
1147 : call pawrhoij_alloc(pawrhoij, cplex_rhoij, dtset%nspden, dtset%nspinor, &
1148 0 : & dtset%nsppol, dtset%typat, pawtab=this%pawtab)
1149 :
1150 : !Compute pawrhoij
1151 0 : if (this%usecprj == 1) then ! cprj is saved in memory
1152 : call pawmkrhoij(this%atindx, this%atindx1, this%cprj, this%dimcprj, dtset%istwfk, dtset%kptopt, dtset%mband,&
1153 : & mband_cprj, this%mcprj, dtset%mkmem, mpi_enreg, dtset%natom, dtset%nband, dtset%nkpt, dtset%nspden, &
1154 : & dtset%nspinor, dtset%nsppol, weights, dtset%paral_kgb, dummy_paw_dmft, pawrhoij, this%unpaw, &
1155 0 : & dtset%usewvl, dtset%wtk)
1156 : else ! cprj is computed on the fly
1157 0 : mcprj_tmp = my_nspinor * mband_cprj * dtset%mkmem * dtset%nsppol
1158 0 : ABI_MALLOC(cprj_tmp, (dtset%natom, mcprj_tmp))
1159 0 : call pawcprj_alloc(cprj_tmp, 0, this%dimcprj)
1160 : call ctocprj(this%atindx, this%cg, 1, cprj_tmp, this%gmet, this%gprimd, 0, 0, 0, dtset%istwfk, this%kg, &
1161 : & dtset%kptns, mcg, mcprj_tmp, dtset%mgfft, dtset%mkmem, mpi_enreg, this%psps%mpsang, dtset%mpw, &
1162 : & dtset%natom, this%nattyp, dtset%nband, dtset%natom, dtset%ngfft, dtset%nkpt, dtset%nloalg, &
1163 : & this%npwarr, dtset%nspinor, dtset%nsppol, dtset%nsppol, dtset%ntypat, dtset%paral_kgb, this%ph1d, &
1164 0 : & this%psps, this%rmet, dtset%typat, this%ucvol, this%unpaw, this%xred, this%ylm, dummy_ylmgr)
1165 : call pawmkrhoij(this%atindx, this%atindx1, cprj_tmp, this%dimcprj, dtset%istwfk, dtset%kptopt, &
1166 : & dtset%mband, mband_cprj, mcprj_tmp, dtset%mkmem, mpi_enreg, dtset%natom, dtset%nband, dtset%nkpt, &
1167 : & dtset%nspden, dtset%nspinor, dtset%nsppol, weights, dtset%paral_kgb, dummy_paw_dmft, pawrhoij, &
1168 0 : & this%unpaw, dtset%usewvl, dtset%wtk)
1169 0 : call pawcprj_free(cprj_tmp)
1170 0 : ABI_FREE(cprj_tmp)
1171 : end if
1172 :
1173 : !Compute the total weighted density (adding PAW-correction).
1174 0 : cplex = 1
1175 0 : dummy_int=0
1176 0 : qphon = 0
1177 : call pawmkrho(1, compch_fft, cplex, this%gprimd, dummy_int, this%indsym, dummy_int, mpi_enreg, &
1178 : & mpi_enreg%my_natom, dtset%natom, dtset%nspden, dtset%nsym, dtset%ntypat, dtset%paral_kgb, this%pawang, &
1179 : & this%pawfgr, this%pawfgrtab, dtset%pawprtvol, pawrhoij, pawrhoij, this%pawtab, qphon, w_rhowfg, &
1180 0 : & w_rhowfr, w_rhor, this%rprimd, dtset%symafm, this%symrec, dtset%typat, this%ucvol, dtset%usewvl, this%xred)
1181 :
1182 0 : call pawrhoij_free(pawrhoij)
1183 :
1184 : ! Second option : Transfer density from coarse to fine grid without rhoij corrections.
1185 0 : elseif (this%nfftprc == this%pawfgr%nfft) then
1186 : !Transfering the weighted density to the fine (PAW) grid, no rhoij correction added.
1187 0 : cplex = 1
1188 0 : optgrid = 1 ! coarse to fine
1189 0 : optin = 0 ! real space
1190 0 : optout = 0 !
1191 0 : ABI_MALLOC(dummy_rhog, (2, this%pawfgr%nfftc))
1192 0 : ABI_MALLOC(dummy_rhogf, (2, this%pawfgr%nfft))
1193 : call transgrid(cplex, mpi_enreg, dtset%nspden, optgrid, optin, optout, dtset%paral_kgb, this%pawfgr, &
1194 0 : & dummy_rhog, dummy_rhogf, w_rhowfr, w_rhor)
1195 0 : ABI_FREE(dummy_rhog)
1196 0 : ABI_FREE(dummy_rhogf)
1197 : !call symrhg(1, this%gprimd, this%irrzon, mpi_enreg, this%nfftprc, ?nfftot, dtset%ngfft, 1, dtset%nsppol, dtset%nsym, &
1198 : !this%phnons, rhog, rhor, this%rprimd, dtset%symafm, dtset%symrel, dtset%tnons)
1199 : else
1200 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): nfftprc /= pawfgr%nfft in PAW not implemented.")
1201 : end if
1202 :
1203 : end if
1204 :
1205 : !With collinear spins the weighted density is not returned in the Pauli (tot/spin) basis by mkrho.
1206 18 : if (dtset%nspden == 2) then
1207 : !spin = 2up - tot
1208 60768 : w_rhor(:, 2) = 2*w_rhor(:, 2) - w_rhor(:, 1)
1209 : end if
1210 :
1211 18 : ABI_FREE(w_rhowfr)
1212 18 : ABI_FREE(w_rhowfg)
1213 :
1214 18 : end subroutine compute_weighted_density
1215 :
1216 : !****f* m_precon/compute_ldos
1217 : !! NAME
1218 : !! compute_ldos
1219 : !!
1220 : !! FUNCTION
1221 : !! Compute the local density of states defined as
1222 : !! ldos = sum_nk f'_nk |u_nk|^2 .
1223 : !! where f'_nk is the derivative of the occupation (nk) with respect to the fermi energie.
1224 : !! When 'nspden'>1, the ldos is returned in the Pauli-basis.
1225 : !!
1226 : !! INPUTS
1227 : !! dtset = all input variables for this dataset
1228 : !! mpi_enreg = informations about MPI parallelization
1229 : !!
1230 : !! OUTPUT
1231 : !! ldos = local density of state
1232 : !!
1233 : !! SOURCE
1234 18 : subroutine compute_ldos(this, dtset, mpi_enreg, ldos)
1235 :
1236 : !Arguments ------------------------------------
1237 : !scalars
1238 : !scalars
1239 : class(precon_object), intent(in) :: this
1240 : type(dataset_type), intent(in) :: dtset
1241 : type(MPI_type), intent(in) :: mpi_enreg
1242 : !arrays
1243 : real(dp), intent(out) :: ldos(:, :)
1244 :
1245 : !Local variables-------------------------------
1246 : !scalars
1247 : integer :: maxocc, i_eigen
1248 : !integer :: ikpt, iband, isppol, nband_k, i_eigen
1249 : !arrays
1250 18 : real(dp), allocatable :: ldos_weights(:)
1251 :
1252 : ! *************************************************************************
1253 :
1254 : !compute weights
1255 54 : ABI_MALLOC(ldos_weights, (dtset%mband*dtset%nkpt*dtset%nsppol))
1256 18 : maxocc = two / (dtset%nsppol * dtset%nspinor) !Maximum number of occupations (1 or 2)
1257 1170 : ldos_weights = 0
1258 :
1259 1170 : do i_eigen = 1, dtset%mband*dtset%nkpt*dtset%nsppol
1260 1170 : ldos_weights(i_eigen) = -derivative_occ(dtset%occopt, this%eigen(i_eigen), this%fermie, dtset%tsmear) * maxocc
1261 : end do
1262 :
1263 : !Compute ldos using mkrho with ldos_weights in place of the occupations
1264 18 : call compute_weighted_density(this, dtset, mpi_enreg, ldos_weights, ldos)
1265 18 : ABI_FREE(ldos_weights)
1266 :
1267 18 : end subroutine compute_ldos
1268 :
1269 : !****f* m_precon/apply_chi0_dfermie
1270 : !! NAME
1271 : !! apply_chi0_dfermie
1272 : !!
1273 : !! FUNCTION
1274 : !!
1275 : !!
1276 : !! INPUTS
1277 : !! dtset = All input variables for this dataset.
1278 : !! mpi_enreg = Information about MPI parallelization.
1279 : !!
1280 : !! SIDE EFFECTS
1281 : !! vec_r (nfftprc, nspden) = Vector (in real space) to which the model chi0 operator is applied (in place).
1282 : !! When nspden > 1 vec_r is in the Pauli spin-basis.
1283 : !!
1284 : !! SOURCE
1285 90 : subroutine apply_chi0_dfermie(this, dtset, vec_r)
1286 :
1287 : !Arguments ------------------------------------
1288 : class(precon_object), intent(in) :: this
1289 : !scalars
1290 : type(dataset_type),intent(in) :: dtset
1291 : !arrays
1292 : real(dp), intent(inout) :: vec_r(this%nfftprc, dtset%nspden)
1293 :
1294 : !Local variables-------------------------------
1295 : integer :: ispden, jspden
1296 : real(dp) :: delta_fermie
1297 :
1298 : ! *************************************************************************
1299 :
1300 : ! Precompute the dot product between the ldos and vec for each spin coordinate
1301 90 : delta_fermie = zero
1302 270 : do jspden = 1, dtset%nspden
1303 607770 : delta_fermie = delta_fermie + 1/this%tdos * dot_product(this%ldos(:, jspden), vec_r(:, jspden)) * this%dvol
1304 : end do
1305 270 : do ispden = 1, dtset%nspden
1306 607770 : vec_r(:, ispden) = delta_fermie * this%ldos(:, ispden)
1307 : end do
1308 :
1309 90 : end subroutine apply_chi0_dfermie
1310 :
1311 : !****f* m_precon/apply_chi0_ldos
1312 : !! NAME
1313 : !! apply_chi0_ldos
1314 : !!
1315 : !! FUNCTION
1316 : !! Apply the ldos model chi0 operator to the vector vec_r (in place) in the Pauli basis.
1317 : !!
1318 : !! INPUTS
1319 : !! dtset =
1320 : !!
1321 : !! SIDE EFFECTS
1322 : !! vec_r (nfftprc, nspden) = Vector (in direct space) to which the model chi0 operator is applied (in place).
1323 : !! When nspden > 1 vec_r is in the Pauli basis.
1324 : !!
1325 : !! SOURCE
1326 90 : subroutine apply_chi0_ldos(this, dtset, vec_r)
1327 :
1328 : !Arguments ------------------------------------
1329 : class(precon_object), intent(in) :: this
1330 : type(dataset_type),intent(in) :: dtset
1331 : !arrays
1332 : real(dp), intent(inout) :: vec_r(this%nfftprc, dtset%nspden)
1333 :
1334 : !Local variables-------------------------------
1335 : !scalars
1336 : integer :: ispden
1337 : !arrays
1338 90 : real(dp), allocatable :: work_r(:, :)
1339 :
1340 : ! *************************************************************************
1341 :
1342 90 : if (abs(this%tdos) > epsilon(this%tdos)) then !Checking that tdos is not 0.
1343 360 : ABI_MALLOC(work_r, (this%nfftprc, dtset%nspden))
1344 :
1345 : !1) chi0(v)(r)_1 = -sum_ispden ldos_ispden(r)*v_ispden(r).
1346 303840 : work_r(:, 1) = 0
1347 270 : do ispden = 1, dtset%nspden
1348 607770 : work_r(:, 1) = work_r(:, 1) - this%ldos(:, ispden)*vec_r(:, ispden)
1349 : end do
1350 : !2) chi0(v)(r)_ispden = -ldos_ispden(r)*v_1(r) for ispden>1.
1351 180 : do ispden = 2, dtset%nspden
1352 303930 : work_r(:, ispden) = -this%ldos(:, ispden)*vec_r(:, 1)
1353 : end do
1354 : !3) Apply the part comming from the variations of the Fermi-level.
1355 90 : call apply_chi0_dfermie(this, dtset, vec_r)
1356 607770 : vec_r = work_r + vec_r
1357 :
1358 90 : ABI_FREE(work_r)
1359 : else
1360 0 : vec_r = 0
1361 : end if
1362 :
1363 90 : end subroutine apply_chi0_ldos
1364 :
1365 : !!***
1366 : !!****f* ABINIT/apply_adjdielmat_ldos
1367 : !! NAME
1368 : !! apply_adjdielmat
1369 : !!
1370 : !! FUNCTION
1371 : !! Apply the ldos- adjoint dielectric matrix I-chi0_ldos*vc to the density rho_r (given in the direct space).
1372 : !!
1373 : !! INPUTS
1374 : !! dtset = All input variables for this dataset.
1375 : !! mpi_enreg = Information about MPI parallelization.
1376 : !! rho_r = Density vector (in direct space).
1377 : !!
1378 : !! OUTPUT
1379 : !! adjdielmat_rho_r = adjdielmat * rho_r
1380 : !!
1381 : !! NOTES
1382 : !!
1383 : !! SOURCE
1384 90 : subroutine apply_adjdielmat_ldos(this, dtset, mpi_enreg, rho_r, adjdielmat_rho_r)
1385 :
1386 : !Arguments ------------------------------------
1387 : class(precon_object) :: this
1388 : !scalars
1389 : type(dataset_type),intent(in) :: dtset
1390 : type(MPI_type),intent(in) :: mpi_enreg
1391 : !arrays
1392 : real(dp), intent(in) :: rho_r(this%nfftprc, dtset%nspden)
1393 : real(dp), intent(inout) :: adjdielmat_rho_r(this%nfftprc, dtset%nspden)
1394 :
1395 : ! *************************************************************************
1396 :
1397 607770 : adjdielmat_rho_r = rho_r
1398 : !1) Apply vc (in the Pauli basis)
1399 90 : call apply_vc(this, dtset, mpi_enreg, adjdielmat_rho_r)
1400 : !2) Apply chi0_ldos (in the Pauli basis)
1401 90 : call apply_chi0_ldos(this, dtset, adjdielmat_rho_r)
1402 : !3) adjdielmat_rho_r = rho_r - vc * chi0 * rho_r = adjdielmat * rho_r
1403 607770 : adjdielmat_rho_r = rho_r - adjdielmat_rho_r
1404 :
1405 90 : end subroutine apply_adjdielmat_ldos
1406 :
1407 : !!***
1408 : !!****f* ABINIT/apply_dielmat_ldos
1409 : !! NAME
1410 : !! apply_adjdielmat
1411 : !!
1412 : !! FUNCTION
1413 : !! Apply the ldos- dielectric matrix I-vc*chi0_ldos to the potential v_r (given in the direct space).
1414 : !!
1415 : !! INPUTS
1416 : !! dtset = All input variables for this dataset.
1417 : !! mpi_enreg = Information about MPI parallelization.
1418 : !! v_r = Density vector (in direct space).
1419 : !!
1420 : !! OUTPUT
1421 : !! dielmat_v_r = dielmat * v_r
1422 : !!
1423 : !! NOTES
1424 : !!
1425 : !! SOURCE
1426 0 : subroutine apply_dielmat_ldos(this, dtset, mpi_enreg, v_r, dielmat_v_r)
1427 :
1428 : !Arguments ------------------------------------
1429 : class(precon_object) :: this
1430 : !scalars
1431 : type(dataset_type),intent(in) :: dtset
1432 : type(MPI_type),intent(in) :: mpi_enreg
1433 : !arrays
1434 : real(dp), intent(in) :: v_r(this%nfftprc, dtset%nspden)
1435 : real(dp), intent(inout) :: dielmat_v_r(this%nfftprc, dtset%nspden)
1436 :
1437 : ! *************************************************************************
1438 :
1439 0 : dielmat_v_r = v_r
1440 : !1) Apply chi0_ldos (in the Pauli basis)
1441 0 : call apply_chi0_ldos(this, dtset, dielmat_v_r)
1442 : !2) Apply vc (in the Pauli basis)
1443 0 : call apply_vc(this, dtset, mpi_enreg, dielmat_v_r)
1444 : !3) dielmat_v_r = v_r - vc * chi0 * v_r = dielmat * v_r
1445 0 : dielmat_v_r = v_r - dielmat_v_r
1446 :
1447 0 : end subroutine apply_dielmat_ldos
1448 :
1449 : !****f* m_precon/get_eigen_index
1450 : !! NAME
1451 : !! get_eigen_index
1452 : !!
1453 : !! FUNCTION
1454 : !! Return the index of the eigenvalue corresponding to (iband, ikpt, isppol)
1455 : !! in the flat 'eigen' array.
1456 : !!
1457 : !! INPUTS
1458 : !! dtset = All input variables for this dataset.
1459 : !! iband = Band index.
1460 : !! ikpt = K-point index.
1461 : !! isppol = Spin-polarization index.
1462 : !!
1463 : !! OUTPUT
1464 : !! i_eigen = Index of (iband, ikpt, isppol) in the eigen array.
1465 : !!
1466 : !! SOURCE
1467 10746 : function get_eigen_index(dtset, iband, ikpt, isppol) result(i_eigen)
1468 :
1469 : !Arguments ------------------------------------
1470 : !scalars
1471 : type(dataset_type),intent(in) :: dtset
1472 : integer, intent(in) :: iband, ikpt, isppol
1473 :
1474 : !Returned variable-------------------------------
1475 : integer :: i_eigen
1476 :
1477 : ! *************************************************************************
1478 :
1479 10746 : i_eigen = iband + (ikpt-1)*dtset%mband + (isppol-1)*dtset%mband*dtset%nkpt
1480 :
1481 : end function get_eigen_index
1482 :
1483 : !****f* m_precon/compute_cg_indices
1484 : !! NAME
1485 : !! compute_cg_indices
1486 : !!
1487 : !! FUNCTION
1488 : !! Compute the index ranges of each wavefunction (iband, ikpt, isppol) in the
1489 : !! flat 'cg' array. For nspinor=2, separate ranges are stored for the
1490 : !! spin-up and spin-down spinor components.
1491 : !!
1492 : !! INPUTS
1493 : !! dtset = All input variables for this dataset.
1494 : !! mpi_enreg = Information about MPI parallelization.
1495 : !! npwarr = Number of plane-waves at each k-point.
1496 : !!
1497 : !! OUTPUTS
1498 : !! cg_indices(2*nspinor, mband, nkpt, nsppol) = Index ranges in the cg array.
1499 : !! cg_indices(1, iband, ikpt, isppol) : start index of wavefunction (iband, ikpt, isppol)
1500 : !! (spin-up spinor component if nspinor=2).
1501 : !! cg_indices(2, iband, ikpt, isppol) : end index.
1502 : !! cg_indices(3, iband, ikpt, isppol) : start index of the spin-down spinor component (nspinor=2 only).
1503 : !! cg_indices(4, iband, ikpt, isppol) : end index of the spin-down spinor component (nspinor=2 only).
1504 : !!
1505 : !! SOURCE
1506 10 : subroutine compute_cg_indices(dtset, mpi_enreg, npwarr, cg_indices)
1507 :
1508 : !Arguments ------------------------------------
1509 : type(dataset_type),intent(in) :: dtset
1510 : type(MPI_type), intent(in) :: mpi_enreg
1511 : integer, intent(in) :: npwarr(:)
1512 : integer :: cg_indices(2*dtset%nspinor, dtset%mband, dtset%nkpt, dtset%nsppol)
1513 : ! If nspinor=2 : cg_indices(1, iband, ikpt, isppol):cg_indices(2, iband, ikpt, isppol) is the range of the spin up
1514 : ! cg_indices(3, iband, ikpt, isppol):cg_indices(4, iband, ikpt, isppol) is the range of the spin down
1515 : ! for the band of indices (iband, ikpt, isppol).
1516 :
1517 : !Local variables ------------------------------
1518 : integer :: iband, ikpt, isppol, i_cg
1519 :
1520 : ! *************************************************************************
1521 :
1522 1990 : cg_indices = zero
1523 10 : i_cg = 1
1524 30 : do isppol =1, dtset%nsppol
1525 70 : do ikpt = 1, dtset%nkpt
1526 700 : do iband = 1, dtset%nband(ikpt)
1527 640 : if (proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, iband, iband, isppol, mpi_enreg%me_kpt)) then
1528 : cycle
1529 : end if
1530 640 : cg_indices(1, iband, ikpt, isppol) = i_cg
1531 640 : cg_indices(2, iband, ikpt, isppol) = cg_indices(1, iband, ikpt, isppol) + npwarr(ikpt) - 1
1532 640 : i_cg = cg_indices(2, iband, ikpt, isppol) + 1
1533 680 : if (dtset%nspinor==2) then
1534 0 : cg_indices(3, iband, ikpt, isppol) = i_cg
1535 0 : cg_indices(4, iband, ikpt, isppol) = cg_indices(3, iband, ikpt, isppol) + npwarr(ikpt) - 1
1536 0 : i_cg = cg_indices(4, iband, ikpt, isppol) + 1
1537 : end if
1538 : end do
1539 : end do
1540 : end do
1541 :
1542 10 : end subroutine compute_cg_indices
1543 :
1544 : !****f* m_precon/compute_kg_indices
1545 : !! NAME
1546 : !! compute_kg_indices
1547 : !!
1548 : !! FUNCTION
1549 : !! Compute the index ranges of the plane-wave coordinates of each k-point
1550 : !! in the flat 'kg' array.
1551 : !!
1552 : !! INPUTS
1553 : !! dtset = All input variables for this dataset.
1554 : !! mpi_enreg = Information about MPI parallelization.
1555 : !! npwarr = Number of plane-waves at each k-point.
1556 : !!
1557 : !! OUTPUTS
1558 : !! kg_indices(2, nkpt) = Index ranges in the kg array.
1559 : !! kg_indices(1, ikpt) : start index of k-point ikpt in the kg array.
1560 : !! kg_indices(2, ikpt) : end index.
1561 : !!
1562 : !! SOURCE
1563 10 : subroutine compute_kg_indices(dtset, mpi_enreg, npwarr, kg_indices)
1564 :
1565 : !Arguments ------------------------------------
1566 : type(dataset_type),intent(in) :: dtset
1567 : type(MPI_type), intent(in) :: mpi_enreg
1568 : integer, intent(in) :: npwarr(:)
1569 : integer, intent(inout) :: kg_indices(2, dtset%nkpt)
1570 : !Returned variable ----------------------------
1571 : integer :: i_kg, ikpt, isppol
1572 :
1573 : ! *************************************************************************
1574 :
1575 70 : kg_indices = zero
1576 30 : do isppol =1, dtset%nsppol
1577 20 : i_kg = 1
1578 70 : do ikpt = 1, dtset%nkpt
1579 40 : if (proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, 1, dtset%nband(ikpt+(isppol-1)*dtset%nkpt), isppol, mpi_enreg%me_kpt)) then
1580 : cycle
1581 : end if
1582 40 : kg_indices(1, ikpt) = i_kg
1583 40 : kg_indices(2, ikpt) = kg_indices(1, ikpt) + npwarr(ikpt) - 1
1584 60 : i_kg = kg_indices(2, ikpt) + 1
1585 : end do
1586 : end do
1587 :
1588 10 : end subroutine compute_kg_indices
1589 :
1590 : !****f* m_precon/transfer_grid
1591 : !! NAME
1592 : !! transfer_grid
1593 : !!
1594 : !! FUNCTION
1595 : !! Wrapper for fftpac and transgrid (PAW).
1596 : !! Transfer one spin component of a real-space density from the augmented
1597 : !! wavefunction FFT grid to the preconditioning FFT grid.
1598 : !! In norm-conserving calculations, the two grids are identical and fftpac is used.
1599 : !! In PAW calculations, the density is first packed to the coarse grid with fftpac
1600 : !! and then interpolated to the fine grid with transgrid.
1601 : !!
1602 : !! INPUTS
1603 : !! dtset = All input variables for this dataset.
1604 : !! mpi_enreg = Information about MPI parallelization.
1605 : !! ispden = Spin component index to write into rho_r.
1606 : !! rho_aug_r(n4, n5, n6) = Density on the augmented wavefunction FFT grid.
1607 : !!
1608 : !! SIDE EFFECTS
1609 : !! rho_r(nfftprc, nspden) = On output, column ispden is filled with the
1610 : !! density transferred to the preconditioning grid.
1611 : !!
1612 : !! SOURCE
1613 216 : subroutine transfer_grid(this, dtset, mpi_enreg, ispden, rho_aug_r, rho_r)
1614 :
1615 : !Arguments ------------------------------------
1616 : class(precon_object), intent(inout) :: this
1617 : !scalars
1618 : type(dataset_type),intent(in) :: dtset
1619 : type(MPI_type), intent(in) :: mpi_enreg
1620 : integer :: ispden
1621 : !arrays
1622 : real(dp), intent(inout) :: rho_r(:, :)
1623 : real(dp), intent(inout) :: rho_aug_r(:, :, :)
1624 :
1625 : !Local variables-------------------------------
1626 : integer :: n1, n2, n3, n4, n5, n6
1627 : integer :: cplex, optgrid, optin, optout
1628 216 : real(dp), allocatable :: rho_coarse_r(:, :)
1629 216 : real(dp), allocatable :: dummy_rhog(:, :), dummy_rhogf(:, :)
1630 :
1631 : ! *************************************************************************
1632 :
1633 216 : n1 = dtset%ngfft(1)
1634 216 : n2 = dtset%ngfft(2)
1635 216 : n3 = dtset%ngfft(3)
1636 216 : n4 = dtset%ngfft(4)
1637 216 : n5 = dtset%ngfft(5)
1638 216 : n6 = dtset%ngfft(6)
1639 :
1640 216 : if (this%psps%usepaw==0) then
1641 : ! In NC, the preconditioning grid should be the density/potential grid.
1642 216 : if (this%nfftprc == n1*n2*n3) then
1643 216 : call fftpac(ispden, mpi_enreg, 1, n1, n2, n3, n4, n5, n6, dtset%ngfft, rho_r, rho_aug_r, 1)
1644 : else
1645 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): nfftprc /= nfft in norm-conserving not implemented.")
1646 : end if
1647 : else
1648 : ! In PAW, the preconditioning grid should be the fine grid.
1649 0 : if (this%nfftprc == this%pawfgr%nfft) then
1650 0 : ABI_MALLOC(rho_coarse_r, (dtset%nfft, 1))
1651 : ! Augmented grid to coarse grid :
1652 0 : call fftpac(1, mpi_enreg, 1, n1, n2, n3, n4, n5, n6, dtset%ngfft, rho_coarse_r, rho_aug_r, 1)
1653 : ! Coarse grid to fine grid :
1654 0 : cplex = 1
1655 0 : optgrid = 1 ! coarse to fine
1656 0 : optin = 0 ! real space
1657 0 : optout = 0 !
1658 0 : ABI_MALLOC(dummy_rhog, (2, this%pawfgr%nfftc))
1659 0 : ABI_MALLOC(dummy_rhogf, (2, this%pawfgr%nfft))
1660 0 : call transgrid(cplex, mpi_enreg, 1, optgrid, optin, optout, dtset%paral_kgb, this%pawfgr, dummy_rhog, dummy_rhogf, rho_coarse_r, rho_r(:, ispden))
1661 0 : ABI_FREE(dummy_rhog)
1662 0 : ABI_FREE(dummy_rhogf)
1663 0 : ABI_FREE(rho_coarse_r)
1664 : else
1665 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): nfftprc /= pawfgr%nfft in PAW not implemented.")
1666 : end if
1667 :
1668 : end if
1669 :
1670 216 : end subroutine transfer_grid
1671 :
1672 : ! Compute rho_i in collinear case without PAW corrections
1673 : !****f* m_precon/compute_rhoi_coll
1674 : !! NAME
1675 : !! compute_rhoi_coll
1676 : !!
1677 : !! FUNCTION
1678 : !! Compute the normalized orbital density
1679 : !! rho_i(r) = |psi_i(r)|^2
1680 : !! for a given band, k-point and spin channel in the collinear case.
1681 : !! The density is transferred to the preconditioning grid and normalized.
1682 : !!
1683 : !! INPUTS
1684 : !! dtset = All input variables for this dataset.
1685 : !! mpi_enreg = Information about MPI parallelization.
1686 : !! iband = Band index.
1687 : !! ikpt = K-point index.
1688 : !! isppol = Spin-polarization index.
1689 : !!
1690 : !! OUTPUTS
1691 : !! rhoi_r(nfftprc,1) = Normalized orbital density on the preconditioning grid.
1692 : !!
1693 : !! SOURCE
1694 216 : subroutine compute_rhoi_coll(this, dtset, mpi_enreg, iband, ikpt, isppol, rhoi_r)
1695 :
1696 : !Arguments ------------------------------------
1697 : class(precon_object), intent(inout) :: this
1698 : !scalars
1699 : type(dataset_type),intent(in) :: dtset
1700 : type(MPI_type), intent(in) :: mpi_enreg
1701 : integer, intent(in) :: isppol, ikpt, iband
1702 : real(dp), intent(inout) :: rhoi_r(this%nfftprc, 1)
1703 :
1704 : !Local variables-------------------------------
1705 : !scalars
1706 : integer :: ndat, option, tim_fourwf
1707 : integer :: i_cg(2), i_kg(2)
1708 : integer :: n1, n2, n3, n4, n5, n6
1709 : integer :: istwf_k, npw_k
1710 : !arrays
1711 432 : integer :: gbound(2*dtset%mgfft+8,2)
1712 216 : integer, allocatable :: kg_k(:, :)
1713 216 : real(dp), allocatable :: rhoi_aug_r(:, :, :)
1714 : !dummy arguments
1715 : integer :: dummy_int
1716 432 : real(dp) :: dummy_fofgout(2, 0), dummy_fofrout(2, dtset%ngfft(4), dtset%ngfft(5), dtset%ngfft(6))
1717 :
1718 : ! *************************************************************************
1719 :
1720 216 : n1 = dtset%ngfft(1)
1721 216 : n2 = dtset%ngfft(2)
1722 216 : n3 = dtset%ngfft(3)
1723 216 : n4 = dtset%ngfft(4)
1724 216 : n5 = dtset%ngfft(5)
1725 216 : n6 = dtset%ngfft(6)
1726 :
1727 : ! No spin or collinear spins - Wafefunctions have one spin component.
1728 216 : if (dtset%nspinor == 1) then
1729 :
1730 : ! 1) Compute orbital density in real space (augmented basis) :
1731 :
1732 : ! Input parameters for fourwf :
1733 216 : option = 1 ! Computes the density.
1734 216 : ndat = 1 ! Only one FFT.
1735 216 : tim_fourwf = 0
1736 1080 : ABI_MALLOC(rhoi_aug_r, (n4, n5, n6))
1737 781056 : rhoi_aug_r = zero ! Initialization for fourwf (accumulation).
1738 216 : istwf_k = dtset%istwfk(ikpt) ! Option parameter that describes the storage of wfs at this kpt.
1739 :
1740 216 : if (.not. (dtset%paral_kgb == 1 .and. dtset%npband > 1)) then
1741 216 : npw_k = this%npwarr(ikpt) ! Number of plane-wave at this kpt.
1742 648 : ABI_MALLOC(kg_k, (3, npw_k))
1743 648 : i_kg = this%kg_indices(:, ikpt)
1744 93344 : kg_k = this%kg(:, i_kg(1):i_kg(2)) ! Reduced plane-wave coordinate (k+G) of this kpt.
1745 216 : call sphereboundary(gbound, istwf_k, kg_k, dtset%mgfft, npw_k)
1746 648 : i_cg = this%cg_indices(:, iband, ikpt, isppol)
1747 : call fourwf(1, rhoi_aug_r, this%cg(:, i_cg(1):i_cg(2)), dummy_fofgout, dummy_fofrout, &
1748 : & gbound, gbound, istwf_k, kg_k, kg_k, dtset%mgfft, mpi_enreg, ndat, dtset%ngfft, npw_k, &
1749 216 : & dummy_int, n4, n5, n6, option, tim_fourwf, one, one)
1750 : else
1751 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): 'compute_rhoi_coll' should not be called with band parallelization.")
1752 : end if
1753 216 : ABI_FREE(kg_k)
1754 :
1755 : ! 2) Transfer rhoi_aug_r defined on the augmented (wavefunction) fft-grid to the preconditioning fft-grid.
1756 216 : call transfer_grid(this, dtset, mpi_enreg, 1, rhoi_aug_r, rhoi_r)
1757 216 : ABI_FREE(rhoi_aug_r)
1758 :
1759 : !3) Normalize rhoi_r.
1760 1458216 : rhoi_r(:, 1) = rhoi_r(:, 1) / (sum(rhoi_r(:, 1)) * this%dvol) !Normalizing rho_ii_r.
1761 :
1762 : else
1763 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): compute_rhoi_coll called with non-collinear magnetism.")
1764 : end if
1765 :
1766 216 : end subroutine compute_rhoi_coll
1767 :
1768 : !****f* m_precon/build_non_coll_density
1769 : !! NAME
1770 : !! build_non_coll_density
1771 : !!
1772 : !! FUNCTION
1773 : !! Construct the four Pauli-basis components of the orbital density
1774 : !! associated with a non-collinear spinor wavefunction and normalize
1775 : !! the resulting density.
1776 : !!
1777 : !! INPUTS
1778 : !! dtset = All input variables for this dataset.
1779 : !! mpi_enreg = Information about MPI parallelization.
1780 : !! psi_r_up = Spin-up component of the wavefunction in real space.
1781 : !! psi_r_down = Spin-down component of the wavefunction in real space.
1782 : !!
1783 : !! OUTPUTS
1784 : !! rhoi_r(nfftprc,4) = Normalized orbital density in the Pauli basis.
1785 : !!
1786 : !! SOURCE
1787 0 : subroutine build_non_coll_density(this, dtset, mpi_enreg, psi_r_up, psi_r_down, rhoi_r)
1788 : !Arguments ------------------------------------
1789 : class(precon_object), intent(inout) :: this
1790 : !scalars
1791 : type(dataset_type),intent(in) :: dtset
1792 : type(MPI_type), intent(in) :: mpi_enreg
1793 : real(dp), intent(in) :: psi_r_up(:, :, :, :), psi_r_down(:, :, :, :)
1794 : real(dp), intent(inout) :: rhoi_r(this%nfftprc, 4)
1795 :
1796 : !Local variables-------------------------------
1797 : !scalars
1798 : integer :: n4, n5, n6
1799 : integer :: ispden
1800 : real(dp) :: norm_tot
1801 : !arrays
1802 0 : real(dp), allocatable :: rhoi_aug_r(:, :, :, :)
1803 :
1804 : ! *************************************************************************
1805 :
1806 0 : n4 = dtset%ngfft(4)
1807 0 : n5 = dtset%ngfft(5)
1808 0 : n6 = dtset%ngfft(6)
1809 :
1810 0 : ABI_MALLOC(rhoi_aug_r, (n4, n5, n6, 4))
1811 0 : ispden = 1 ! rho_sigma0(r) = |psi_up(r)|^2 + |psi_down(r)|^2
1812 0 : rhoi_aug_r(:, :, :, ispden) = psi_r_up(1, :, :, :)**2 + psi_r_up(2, :, :, :)**2 + psi_r_down(1, :, :, :)**2 + psi_r_down(2, :, :, :)**2
1813 : ispden = 2 ! rho_sigma1(r) = psi_up(r)* . psi_down(r) + psi_down(r)* . psi_up(r) = 2 Re(psi_up(r)*psi_down(r))
1814 : ! (* = conjugate)
1815 0 : rhoi_aug_r(:, :, :, ispden) = 2*( psi_r_up(1, :, :, :)*psi_r_down(1, :, :, :) + psi_r_up(2, :, :, :)*psi_r_down(2, :, :, :) )
1816 0 : ispden = 3 ! rho_sigma2(r) = i*(psi_down(r)* psi_up(r) - psi_up(r)* psi_down(r)) = 2 Im(psi_up(r)*psi_down(r))
1817 0 : rhoi_aug_r(:, :, :, ispden) = 2*( psi_r_up(2, :, :, :)*psi_r_down(1, :, :, :) - psi_r_up(1, :, :, :)*psi_r_down(2, :, :, :) )
1818 : ispden = 4 ! rho_sigma3(r) = |psi_up|^2 - |psi_down|^2
1819 0 : rhoi_aug_r(:, :, :, ispden) = psi_r_up(1, :, :, :)**2 + psi_r_up(2, :, :, :)**2 - psi_r_down(1, :, :, :)**2 - psi_r_down(2, :, :, :)**2
1820 :
1821 : ! Change grid
1822 0 : do ispden = 1, 4
1823 0 : call transfer_grid(this, dtset, mpi_enreg, ispden, rhoi_aug_r(:, :, :, ispden), rhoi_r)
1824 : end do
1825 : ! TODO : here deal with cases where nspinor = 2 but nspden != 4 ...
1826 :
1827 0 : ABI_FREE(rhoi_aug_r)
1828 :
1829 : !3) Normalize rhoi_r.
1830 0 : norm_tot = 0.5_dp * sum(rhoi_r(:, 1)) * this%dvol
1831 0 : do ispden = 1, 4
1832 0 : rhoi_r(:, ispden) = rhoi_r(:, ispden) / norm_tot !Normalizing rho_ii_r.
1833 : end do
1834 :
1835 0 : end subroutine build_non_coll_density
1836 :
1837 : !****f* m_precon/compute_rhoi_noncoll
1838 : !! NAME
1839 : !! compute_rhoi_noncoll
1840 : !!
1841 : !! FUNCTION
1842 : !! Compute the normalized orbital density associated with a given
1843 : !! band, k-point and spin channel in the non-collinear case.
1844 : !! The density is returned in the Pauli basis on the preconditioning grid.
1845 : !!
1846 : !! INPUTS
1847 : !! dtset = All input variables for this dataset.
1848 : !! mpi_enreg = Information about MPI parallelization.
1849 : !! iband = Band index.
1850 : !! ikpt = K-point index.
1851 : !! isppol = Spin-polarization index.
1852 : !!
1853 : !! OUTPUTS
1854 : !! rhoi_r(nfftprc,4) = Normalized orbital density in the Pauli basis.
1855 : !!
1856 : !! SOURCE
1857 0 : subroutine compute_rhoi_noncoll(this, dtset, mpi_enreg, iband, ikpt, isppol, rhoi_r)
1858 : !Arguments ------------------------------------
1859 : class(precon_object), intent(inout) :: this
1860 : !scalars
1861 : type(dataset_type),intent(in) :: dtset
1862 : type(MPI_type), intent(in) :: mpi_enreg
1863 : integer, intent(in) :: isppol, ikpt, iband
1864 : real(dp), intent(inout) :: rhoi_r(this%nfftprc, 4)
1865 :
1866 : !Local variables-------------------------------
1867 : !scalars
1868 : integer :: ndat, option, tim_fourwf
1869 : integer :: i_cg(4), i_kg(2)
1870 : integer :: n1, n2, n3, n4, n5, n6
1871 : integer :: istwf_k, npw_k
1872 : !arrays
1873 0 : integer :: gbound(2*dtset%mgfft+8,2)
1874 0 : integer, allocatable :: kg_k(:, :)
1875 0 : real(dp), allocatable :: psi_r_up(:, :, :, :), psi_r_down(:, :, :, :)
1876 : !dummy arguments
1877 : integer :: dummy_int
1878 : real(dp) :: dummy_fofgout(2, 0)
1879 0 : real(dp) :: dummy_denpot(dtset%ngfft(4), dtset%ngfft(5), dtset%ngfft(6))
1880 :
1881 : ! *************************************************************************
1882 :
1883 0 : n1 = dtset%ngfft(1)
1884 0 : n2 = dtset%ngfft(2)
1885 0 : n3 = dtset%ngfft(3)
1886 0 : n4 = dtset%ngfft(4)
1887 0 : n5 = dtset%ngfft(5)
1888 0 : n6 = dtset%ngfft(6)
1889 :
1890 : ! Non collinear spins - Wafefunctions have two spin component.
1891 0 : if (dtset%nspinor == 2) then
1892 :
1893 : !1) Compute psi_up and psi_down in real space :
1894 0 : ABI_MALLOC(psi_r_up, (2, n4, n5, n6))
1895 0 : ABI_MALLOC(psi_r_down, (2, n4, n5, n6))
1896 : ! Input parameters for fourwf :
1897 0 : option = 0 ! Only do the FFT.
1898 0 : ndat = 1
1899 0 : tim_fourwf = 0
1900 0 : i_cg = this%cg_indices(:, iband, ikpt, isppol)
1901 0 : istwf_k = dtset%istwfk(ikpt) ! Option parameter that describes the storage of wfs at this kpt.
1902 0 : npw_k = this%npwarr(ikpt) ! Number of plane-wave at this kpt.
1903 0 : ABI_MALLOC(kg_k, (3, npw_k))
1904 0 : i_kg = this%kg_indices(:, ikpt)
1905 0 : kg_k = this%kg(:, i_kg(1):i_kg(2)) ! Reduced plane-wave coordinate (k+G) of this kpt.
1906 0 : call sphereboundary(gbound, istwf_k, kg_k, dtset%mgfft, npw_k)
1907 : !FFT for psi_up
1908 : call fourwf(dummy_int, dummy_denpot, this%cg(:, i_cg(1):i_cg(2)), dummy_fofgout, psi_r_up, &
1909 : & gbound, gbound, istwf_k, kg_k, kg_k, dtset%mgfft, mpi_enreg, ndat, dtset%ngfft, npw_k, &
1910 0 : & dummy_int, n4, n5, n6, option, tim_fourwf, one, one)
1911 : !FFT for psi_down
1912 : call fourwf(dummy_int, dummy_denpot, this%cg(:, i_cg(3):i_cg(4)), dummy_fofgout, psi_r_down, &
1913 : & gbound, gbound, istwf_k, kg_k, kg_k, dtset%mgfft, mpi_enreg, ndat, dtset%ngfft, npw_k, &
1914 0 : & dummy_int, n4, n5, n6, option, tim_fourwf, one, one)
1915 : ! (done separately for convenience & readability)
1916 0 : ABI_FREE(kg_k)
1917 :
1918 : !2) Build the 4 components of the orbital density rhoi_r:
1919 0 : call build_non_coll_density(this, dtset, mpi_enreg, psi_r_up, psi_r_down, rhoi_r)
1920 0 : ABI_FREE(psi_r_up)
1921 0 : ABI_FREE(psi_r_down)
1922 :
1923 : else
1924 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): compute_rhoi_noncoll called with collinear magnetism.")
1925 : end if
1926 :
1927 0 : end subroutine compute_rhoi_noncoll
1928 :
1929 : !****f* m_precon/cycle_band
1930 : !! NAME
1931 : !! cycle_band
1932 : !!
1933 : !! FUNCTION
1934 : !! Determine whether a band should be skipped on the current MPI process
1935 : !! according to the band-distribution scheme.
1936 : !!
1937 : !! INPUTS
1938 : !! mpi_enreg = Information about MPI parallelization.
1939 : !! nband_k = Number of bands at the current k-point.
1940 : !! iband = Band index.
1941 : !!
1942 : !! OUTPUTS
1943 : !! do_not_belong_to_proc = .true. if the band is not assigned to the
1944 : !! current process, .false. otherwise.
1945 : !!
1946 : !! SOURCE
1947 10044 : function cycle_band(mpi_enreg, nband_k, iband) result(do_not_belong_to_proc)
1948 : !Arguments ------------------------------------
1949 : !type(dataset_type),intent(in) :: dtset
1950 : type(MPI_type), intent(in) :: mpi_enreg
1951 : integer, intent(in) :: iband, nband_k
1952 :
1953 : !Local variables-------------------------------
1954 : !scalars
1955 : integer :: rank, nbdblock, blocksize
1956 : !arrays
1957 :
1958 : !Returned variable-------------------------------
1959 : logical :: do_not_belong_to_proc
1960 :
1961 : ! *************************************************************************
1962 :
1963 10044 : rank = xmpi_comm_rank(mpi_enreg%comm_bandfft)
1964 10044 : nbdblock = nband_k / (mpi_enreg%nproc_band * mpi_enreg%bandpp)
1965 10044 : blocksize = nband_k / nbdblock
1966 : ! Check if this band belong to the current processor (assuming bands are distributed in order).
1967 10044 : do_not_belong_to_proc = .not.(1 + mpi_enreg%bandpp*rank <= mod(iband, blocksize) .and. iband <=mpi_enreg%bandpp*(rank+1))
1968 :
1969 10044 : end function cycle_band
1970 :
1971 : !****f* m_precon/precompute_rhoi
1972 : !! NAME
1973 : !! precompute_rhoi
1974 : !!
1975 : !! FUNCTION
1976 : !! Precompute and store the normalized orbital densities rho_i(r) = |psi_i(r)|^2
1977 : !! for all bands with a non-negligible occupation derivative, to avoid
1978 : !! redundant FFTs during the iterative application of chi0_diag.
1979 : !! Results are stored in this%precomputed_rhoi and indexed by
1980 : !! this%precomputed_rhoi_indices.
1981 : !!
1982 : !! INPUTS
1983 : !! dtset = All input variables for this dataset.
1984 : !! mpi_enreg = Information about MPI parallelization.
1985 : !!
1986 : !! SOURCE
1987 10 : subroutine precompute_rhoi(this, dtset, mpi_enreg)
1988 :
1989 : !Arguments ------------------------------------
1990 : class(precon_object), intent(inout) :: this
1991 : !scalars
1992 : type(dataset_type),intent(in) :: dtset
1993 : type(MPI_type), intent(in) :: mpi_enreg
1994 :
1995 : !Local variables-------------------------------
1996 : !scalars
1997 : integer :: nspin, i_rhoi, isppol, ikpt, i_kpt_sppol, nband_k, iband, iband1, iband2
1998 : !arrays
1999 10 : integer, allocatable :: needed_bands_bounds(:, :)
2000 10 : integer, allocatable :: needed_bands_number(:)
2001 : !for band parall
2002 : integer :: option_fourwf, ndat, blocksize, iblock, ibandblock1, ibandblock2, nbdblock, nfft_blocks
2003 : integer :: n1, n2, n3, n4, n5, n6
2004 : integer :: idat, idat_down
2005 20 : integer :: i_cg_ibandblock1(2*dtset%nspinor), i_cg_ibandblock2(2*dtset%nspinor)
2006 10 : real(dp), allocatable :: dummy_occ_k(:)
2007 10 : real(dp), allocatable :: dummy_denpot(:, :, :)
2008 10 : real(dp), allocatable :: rhoi_aug(:, :, :, :)
2009 10 : real(dp), allocatable :: psii_aug(:, :, :, :)
2010 :
2011 : ! *************************************************************************
2012 :
2013 30 : ABI_MALLOC(needed_bands_bounds, (2, dtset%nkpt*dtset%nsppol))
2014 30 : ABI_MALLOC(needed_bands_number, (dtset%nkpt*dtset%nsppol))
2015 10 : call get_needed_bands_chi0diag(this, dtset, mpi_enreg, needed_bands_bounds, needed_bands_number)
2016 :
2017 : ! Allocate the array containing the precomputed rhoi
2018 10 : if (dtset%nspinor==1) then
2019 : nspin = 1 ! Number of spin components in the orbital densities (rhoi).
2020 0 : else if (dtset%nspinor==2) then
2021 : nspin = 4
2022 : else
2023 0 : ABI_BUG("nspinor /= 1 or 2")
2024 : end if
2025 90 : ABI_MALLOC(this%precomputed_rhoi, (this%nfftprc, nspin, sum(needed_bands_number)))
2026 710 : this%precomputed_rhoi_indices = zero
2027 10 : i_rhoi = 1
2028 :
2029 : !Loop over spins and kpoints
2030 30 : do isppol =1, dtset%nsppol
2031 70 : do ikpt = 1, dtset%nkpt
2032 40 : i_kpt_sppol = ikpt+(isppol-1)*dtset%nkpt
2033 :
2034 : ! MPI parallelization over kpoints : cycle if kpt does not belong to current processor.
2035 40 : nband_k = dtset%nband(i_kpt_sppol)
2036 40 : if (proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, 1, nband_k, isppol, mpi_enreg%me_kpt)) then
2037 : cycle
2038 : end if
2039 :
2040 40 : iband1 = needed_bands_bounds(1, i_kpt_sppol)
2041 40 : iband2 = needed_bands_bounds(2, i_kpt_sppol)
2042 :
2043 60 : if (.not. (dtset%paral_kgb == 1 .and. dtset%npband > 1)) then
2044 : ! No parallelization over band : We loop over needed bands to compute and save the orbital densities.
2045 :
2046 256 : do iband = iband1, iband2
2047 :
2048 : ! No spin or collinear spins - Wafefunctions have one spin component.
2049 216 : if (dtset%nspinor == 1) then
2050 216 : call compute_rhoi_coll(this, dtset, mpi_enreg, iband, ikpt, isppol, this%precomputed_rhoi(:, :, i_rhoi))
2051 216 : this%precomputed_rhoi_indices(iband, ikpt, isppol) = i_rhoi
2052 216 : i_rhoi = i_rhoi + 1
2053 : end if
2054 :
2055 : ! Non collinear spins - Wavefunctions have two spins components.
2056 256 : if (dtset%nspinor == 2) then
2057 0 : call compute_rhoi_noncoll(this, dtset, mpi_enreg, iband, ikpt, isppol, this%precomputed_rhoi(:, :, i_rhoi))
2058 0 : this%precomputed_rhoi_indices(iband, ikpt, isppol) = i_rhoi
2059 0 : i_rhoi = i_rhoi + 1
2060 : end if
2061 :
2062 : end do !iband
2063 :
2064 : else
2065 : ! Parallelization over band : We first do all the needed fft, with 'prep_fourwf' that will take care of
2066 : ! the transpose from the linalg representation to the fft representation.
2067 : ! Then we fill the 'precomputed_rhoi' array band per band.
2068 :
2069 0 : n1 = dtset%ngfft(1)
2070 0 : n2 = dtset%ngfft(2)
2071 0 : n3 = dtset%ngfft(3)
2072 0 : n4 = dtset%ngfft(4)
2073 0 : n5 = dtset%ngfft(5)
2074 0 : n6 = dtset%ngfft(6)
2075 :
2076 0 : nbdblock=nband_k / (mpi_enreg%nproc_band * mpi_enreg%bandpp)
2077 0 : blocksize=nband_k / nbdblock
2078 :
2079 0 : ndat = mpi_enreg%bandpp
2080 0 : ABI_MALLOC(psii_aug, (2, n4, n5, n6*ndat*dtset%nspinor))
2081 0 : ABI_MALLOC(rhoi_aug, (n4, n5, n6, nspin))
2082 :
2083 0 : option_fourwf = 0
2084 :
2085 0 : do iblock = 1, nbdblock ! Loop over (LOBPCG) blocks
2086 :
2087 : ! 1) FFTs
2088 0 : if (dtset%nspinor==1) then
2089 :
2090 0 : ibandblock1 = blocksize*(iblock-1) + 1
2091 0 : ibandblock2 = blocksize*(iblock)
2092 0 : i_cg_ibandblock1 = this%cg_indices(:, ibandblock1, ikpt, isppol)
2093 0 : i_cg_ibandblock2 = this%cg_indices(:, ibandblock2, ikpt, isppol)
2094 :
2095 0 : ABI_MALLOC(dummy_occ_k, (nband_k))
2096 0 : ABI_MALLOC(dummy_denpot, (n4, n5, n6))
2097 :
2098 0 : call bandfft_kpt_set_ikpt(ikpt, mpi_enreg)
2099 0 : nfft_blocks = 1 ! TODO : what is that ??
2100 : call prep_fourwf(dummy_denpot, blocksize, this%cg(:, i_cg_ibandblock1(1):i_cg_ibandblock2(2)), &
2101 : & psii_aug, iblock, dtset%istwfk(ikpt), dtset%mgfft, mpi_enreg, nband_k, &
2102 : & ndat, dtset%ngfft, this%npwarr(ikpt), &
2103 0 : & n4, n5, n6, dummy_occ_k, option_fourwf, this%ucvol, dtset%wtk(ikpt), nfft_blocks)
2104 :
2105 0 : ABI_FREE(dummy_occ_k)
2106 0 : ABI_FREE(dummy_denpot)
2107 : else
2108 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): non-collinear magnetisme with band parall - TODO")
2109 : ! TODO : How is prep_fourwf suppsed to be called with nspinor= = 2 ????
2110 : end if
2111 :
2112 : ! 2) Fill precomputed_rhoi
2113 0 : idat = 0
2114 :
2115 : ! Loop over the bands of this proc (assuming bands are distributed in order).
2116 0 : do iband = ibandblock1, ibandblock2
2117 :
2118 0 : if (cycle_band(mpi_enreg, nband_k, iband)) then
2119 : cycle
2120 : end if ! Checks if the band belongs to the current processor.
2121 :
2122 0 : idat = idat + 1
2123 :
2124 : ! Check if this band is needed
2125 0 : if (.not. (iband1 <= iband .and. iband <= iband2)) then
2126 : cycle
2127 : end if
2128 :
2129 0 : if (dtset%nspinor==1) then
2130 :
2131 : !Compute rhoi_aug from psii_aug
2132 0 : rhoi_aug = zero
2133 0 : call cg_addtorho(n1, n2, n3, n4, n5, n6, 1, one, one, psii_aug(:, :, :, (idat-1)*n6+1:idat*n6), rhoi_aug(:, :, :, 1))
2134 :
2135 : ! Grid transfer
2136 0 : call transfer_grid(this, dtset, mpi_enreg, 1, rhoi_aug(:, :, :, 1), this%precomputed_rhoi(:, :, i_rhoi))
2137 :
2138 : ! Normalize
2139 : this%precomputed_rhoi(:, 1, i_rhoi) = this%precomputed_rhoi(:, 1, i_rhoi) / &
2140 0 : & (sum(this%precomputed_rhoi(:, 1, i_rhoi)) * this%dvol)
2141 :
2142 : else ! Non-collinear case.
2143 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): non-collinear magnetisme with band parall - TODO")
2144 :
2145 : ! This assumes that in the FFT, the two spinorial components consecutively stored.
2146 0 : idat_down = idat +1
2147 : call build_non_coll_density(this, dtset, mpi_enreg, psii_aug(:, :, :, (idat-1)*n6+1:idat*n6), &
2148 0 : & psii_aug(:, :, :, (idat_down-1)*n6+1:idat_down*n6), this%precomputed_rhoi(:, :, i_rhoi))
2149 0 : idat = idat_down ! Because we used two consecutive blocks of psii_aug for the two spin components of the same band.
2150 :
2151 : end if
2152 : ! Save the index
2153 0 : this%precomputed_rhoi_indices(iband, ikpt, isppol) = i_rhoi
2154 0 : i_rhoi = i_rhoi + 1
2155 :
2156 : end do
2157 :
2158 : end do
2159 :
2160 0 : ABI_FREE(rhoi_aug)
2161 0 : ABI_FREE(psii_aug)
2162 :
2163 : end if
2164 :
2165 : end do !ikpt
2166 : end do !isppol
2167 :
2168 10 : ABI_FREE(needed_bands_number)
2169 10 : ABI_FREE(needed_bands_bounds)
2170 :
2171 10 : end subroutine precompute_rhoi
2172 :
2173 : !****f* m_precon/get_needed_bands_chi0diag
2174 : !! NAME
2175 : !! get_needed_bands_chi0diag
2176 : !!
2177 : !! FUNCTION
2178 : !! Determine the range of bands that have a non-negligible occupation derivative
2179 : !! f'(e_nk - e_F) at each (ikpt, isppol), i.e. the bands that contribute to
2180 : !! chi0_diag. Bands outside this range are skipped in the loops of
2181 : !! compute_delta_occ.
2182 : !!
2183 : !! INPUTS
2184 : !! dtset = All input variables for this dataset.
2185 : !! mpi_enreg = Information about MPI parallelization.
2186 : !!
2187 : !! OUTPUTS
2188 : !! needed_bands_bounds(2, nkpt*nsppol) = For each (ikpt, isppol), the minimum
2189 : !! (index 1) and maximum (index 2) band
2190 : !! indices with |f'| > deigvals_tol_fp.
2191 : !! needed_bands_number(nkpt*nsppol) = Number of needed bands at each (ikpt, isppol).
2192 : !!
2193 : !! SOURCE
2194 120 : subroutine get_needed_bands_chi0diag(this, dtset, mpi_enreg, needed_bands_bounds, needed_bands_number)
2195 : !Arguments ------------------------------------
2196 : class(precon_object), intent(in) :: this
2197 : !scalars
2198 : type(dataset_type),intent(in) :: dtset
2199 : type(MPI_type), intent(in) :: mpi_enreg
2200 : !arrays
2201 : integer :: needed_bands_number(dtset%nsppol*dtset%nkpt)
2202 : integer :: needed_bands_bounds(2, dtset%nsppol*dtset%nkpt)
2203 :
2204 : !Local variables-------------------------------
2205 : !scalars
2206 : integer :: i_eigen, i_kpt_sppol, ikpt, isppol, iband, nband_k
2207 : real(dp) :: fp, maxocc
2208 :
2209 : ! *************************************************************************
2210 120 : maxocc = two / (dtset%nsppol * dtset%nspinor) !Maximum number of occupations (1 or 2)
2211 600 : needed_bands_number = zero
2212 :
2213 360 : do isppol =1, dtset%nsppol
2214 840 : do ikpt = 1, dtset%nkpt
2215 :
2216 480 : i_kpt_sppol = ikpt+(isppol-1)*dtset%nkpt
2217 480 : nband_k = dtset%nband(i_kpt_sppol)
2218 480 : needed_bands_bounds(1, i_kpt_sppol) = nband_k + 1
2219 480 : needed_bands_bounds(2, i_kpt_sppol) = 0
2220 :
2221 480 : if (proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, 1, nband_k, isppol, mpi_enreg%me_kpt)) then
2222 : cycle
2223 : end if
2224 :
2225 8160 : do iband = 1, nband_k
2226 :
2227 7680 : if (cycle_band(mpi_enreg, nband_k, iband)) then
2228 : cycle
2229 : end if ! Check if the band belongs to the current processor.
2230 :
2231 7200 : i_eigen = get_eigen_index(dtset, iband, ikpt, isppol) ! Index of (iband, ikpt, isppol) in eigen array.
2232 7200 : fp = derivative_occ(dtset%occopt, this%eigen(i_eigen), this%fermie, this%precon_tsmear) * maxocc
2233 :
2234 7680 : if (abs(fp) > this%deigvals_tol_fp) then
2235 2580 : needed_bands_bounds(1, i_kpt_sppol) = min(needed_bands_bounds(1, i_kpt_sppol), iband) ! iband_min
2236 2580 : needed_bands_bounds(2, i_kpt_sppol) = max(needed_bands_bounds(2, i_kpt_sppol), iband) ! iband_max
2237 : end if
2238 :
2239 : end do
2240 :
2241 720 : needed_bands_number(i_kpt_sppol) = max(needed_bands_bounds(2, i_kpt_sppol) - needed_bands_bounds(1, i_kpt_sppol) + 1, 0)
2242 :
2243 : end do !ikpt
2244 : end do !isppol
2245 :
2246 120 : end subroutine get_needed_bands_chi0diag
2247 :
2248 : !****f* m_precon/compute_delta_occ
2249 : !! NAME
2250 : !! compute_delta_occ
2251 : !!
2252 : !! FUNCTION
2253 : !! Compute the first-order variation of the band occupations induced by a
2254 : !! potential perturbation delta_V, within the diagonal approximation of chi0:
2255 : !! delta_occ(nk) = f'(e_nk - e_F) * <psi_nk| delta_V |psi_nk>
2256 : !! A Fermi-level correction is then applied to enforce electron-number conservation:
2257 : !! delta_occ(nk) -= f'(e_nk - e_F) * delta_e_F
2258 : !! where delta_e_F = sum_nk delta_occ(nk) / DOS(e_F).
2259 : !!
2260 : !! INPUTS
2261 : !! dtset = All input variables for this dataset.
2262 : !! mpi_enreg = Information about MPI parallelization.
2263 : !! delta_V(nfftprc, nspden) = Potential perturbation in the Pauli basis.
2264 : !!
2265 : !! OUTPUTS
2266 : !! delta_occ(mband*nkpt*nsppol) = First-order variation of the occupations.
2267 : !!
2268 : !! SOURCE
2269 55 : subroutine compute_delta_occ(this, dtset, mpi_enreg, delta_V, delta_occ)
2270 : !Arguments ------------------------------------
2271 : class(precon_object), intent(inout) :: this
2272 : !scalars
2273 : type(dataset_type),intent(in) :: dtset
2274 : type(MPI_type), intent(in) :: mpi_enreg
2275 : !arrays
2276 : real(dp), intent(in) :: delta_V(this%nfftprc, dtset%nspden) ! In Pauli basis
2277 : real(dp), intent(inout) :: delta_occ(:)
2278 :
2279 : !Local variables-------------------------------
2280 : !scalars
2281 : integer :: nband_k, nspin
2282 : integer :: i_eigen, ikpt, iband, isppol, ier, ispden
2283 : integer :: iband1, iband2, i_kpt_sppol
2284 : real(dp) :: fp, eigenval, maxocc
2285 : real(dp) :: dos_fermie, delta_occ_tot, delta_fermie
2286 : !arrays
2287 55 : integer, allocatable :: needed_bands_bounds(:, :)
2288 55 : integer, allocatable :: needed_bands_number(:)
2289 55 : real(dp), allocatable :: rhoi_r(:, :)
2290 : !real(dp), allocatable :: doccde(:), occ(:)
2291 :
2292 : ! *************************************************************************
2293 :
2294 165 : ABI_MALLOC(needed_bands_bounds, (2, dtset%nkpt*dtset%nsppol))
2295 165 : ABI_MALLOC(needed_bands_number, (dtset%nkpt*dtset%nsppol))
2296 55 : call get_needed_bands_chi0diag(this, dtset, mpi_enreg, needed_bands_bounds, needed_bands_number)
2297 :
2298 : ! Compute the delta_occ = fi' * <rhoii, vec>
2299 3575 : delta_occ = zero
2300 55 : maxocc = two / (dtset%nsppol * dtset%nspinor) !Maximum number of occupations (1 or 2)
2301 :
2302 : ! Allocate the arrays that will contain rhoi
2303 55 : if (dtset%nspinor==1) then
2304 : nspin = 1 ! Number of spin components in the orbital densities (rhoi).
2305 0 : else if (dtset%nspinor==2) then
2306 : nspin = 4
2307 : else
2308 0 : ABI_BUG("nspinor /= 1 or 2")
2309 : end if
2310 220 : ABI_MALLOC(rhoi_r, (this%nfftprc, nspin))
2311 55 : dos_fermie = zero
2312 55 : delta_occ_tot = zero
2313 :
2314 : ! Occupation derivatives could also be computetd with 'getnel' ... What is best ?
2315 : !ABI_MALLOC(doccde, (size(this%occ)))
2316 : !ABI_MALLOC(occ, (size(this%occ)))
2317 : !option=1
2318 : !call getnel(doccde, dummy_real, this%eigen, entropy, this%fermie, this%fermie, maxocc, &
2319 : !& dtset%mband, dtset%nband, nelect, dtset%nkpt, dtset%nsppol, occ, dtset%occopt, &
2320 : !& option, dtset%tphysel, dtset%tsmear, dummy_int, dtset%wtk)
2321 : !ABI_FREE(occ)
2322 :
2323 : ! 1) Eigenvalue variations
2324 : !Loop over spins and kpoints
2325 165 : do isppol =1, dtset%nsppol
2326 :
2327 385 : do ikpt = 1, dtset%nkpt
2328 :
2329 220 : i_kpt_sppol = ikpt+(isppol-1)*dtset%nkpt
2330 220 : nband_k = dtset%nband(i_kpt_sppol)
2331 :
2332 : ! MPI parallelization over kpoints : cycle if kpt does not belong to current processor.
2333 220 : if (proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, 1, nband_k, isppol, mpi_enreg%me_kpt)) then
2334 : cycle
2335 : end if
2336 :
2337 220 : iband1 = needed_bands_bounds(1, i_kpt_sppol)
2338 220 : iband2 = needed_bands_bounds(2, i_kpt_sppol)
2339 :
2340 1512 : do iband = iband1, iband2 ! Loop over needed bands
2341 :
2342 1182 : if (cycle_band(mpi_enreg, nband_k, iband)) then
2343 : cycle
2344 : end if ! Check if this band belong to current processor.
2345 :
2346 : !Indices
2347 1182 : i_eigen = get_eigen_index(dtset, iband, ikpt, isppol) ! Index of (iband, ikpt, isppol) in eigen array.
2348 :
2349 : !2.1) Computing f'(eig_i - fermie).
2350 1182 : eigenval = this%eigen(i_eigen)
2351 1182 : fp = derivative_occ(dtset%occopt, eigenval, this%fermie, this%precon_tsmear) * maxocc
2352 : !fp = doccde(i_eigen) ! Same as derivative_occ (probably more robust)
2353 1182 : dos_fermie = dos_fermie + fp * dtset%wtk(ikpt)
2354 :
2355 : ! No spin or collinear spins - Wafefunctions have one spin component.
2356 1182 : if (dtset%nspinor == 1) then
2357 :
2358 : !2.2) Computing rho_i = |psi_i|^2 using fourwf (if fp is not 0).
2359 1182 : if (this%use_precomputed_rhoi) then
2360 3992796 : rhoi_r = this%precomputed_rhoi(:, :, this%precomputed_rhoi_indices(iband, ikpt, isppol)) ! TODO : useless copy here, use a pointer ?
2361 : else
2362 0 : call compute_rhoi_coll(this, dtset, mpi_enreg, iband, ikpt, isppol, rhoi_r)
2363 : end if
2364 :
2365 : !2.3) delta_occ(i) = fp_i * dot(rho_i, delta_V)
2366 : ! dot-product in the up/down basis (equal to the dot product in the Pauli basis) :
2367 : ! rhoi_r has only one spin-component corresponding to isppol (up=down if nsppol=1, up or down if nsppol=2).
2368 1182 : if (dtset%nspden == 1) then
2369 0 : delta_occ(i_eigen) = fp * dot_product(rhoi_r(:, 1), delta_V(:, isppol)) * this%dvol
2370 1182 : elseif (dtset%nspden == 2) then
2371 3990432 : delta_occ(i_eigen) = fp * dot_product(rhoi_r(:, 1), delta_V(:, 1) + (1-2*(isppol-1)) * delta_V(:, 2)) * this%dvol
2372 : ! delta_V(:, 1) + (1-2*(isppol-1)) * delta_V(:, 2) is delta_V in the (up/down) coordinate 'isppol'.
2373 : end if
2374 :
2375 : ! Non collinear spins - Wavefunctions have two spins components.
2376 0 : elseif (dtset%nspinor == 2 .and. dtset%nspden==4) then
2377 :
2378 : !2.2) Computing rho_i (4-dim, in pauli basis).
2379 0 : if (this%use_precomputed_rhoi) then
2380 0 : rhoi_r = this%precomputed_rhoi(:, :, this%precomputed_rhoi_indices(iband, ikpt, isppol)) ! TODO : useless copy here, use a pointer ?
2381 : else
2382 0 : call compute_rhoi_noncoll(this, dtset, mpi_enreg, iband, ikpt, isppol, rhoi_r)
2383 : end if
2384 :
2385 : ! dot product in Pauli basis :
2386 0 : delta_occ(i_eigen) = zero
2387 0 : do ispden = 1, 4
2388 : ! rhoi_r has 4 spin-components in the pauli basis that all needs to be multiplied to the corresponding component in delta_V
2389 0 : delta_occ(i_eigen) = delta_occ(i_eigen) + fp * dot_product(rhoi_r(:, ispden), delta_V(:, ispden)) * this%dvol
2390 : end do
2391 :
2392 : else
2393 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): TODO non-collinear magnetism with nspden=/4.")
2394 : ! The calculation above is probably true even when nspden=1, TODO: check
2395 : end if
2396 :
2397 1402 : delta_occ_tot = delta_occ_tot + delta_occ(i_eigen) * dtset%wtk(ikpt)
2398 :
2399 : end do
2400 :
2401 : end do !ikpt
2402 : end do !isppol
2403 :
2404 55 : ABI_FREE(rhoi_r)
2405 :
2406 : !MPI parallelization over kpoints : sum delta_occ on all processors.
2407 : ier = 0
2408 55 : call xmpi_sum(delta_occ, mpi_enreg%comm_kptband, ier)
2409 55 : call xmpi_sum(delta_occ_tot, mpi_enreg%comm_kptband, ier)
2410 55 : call xmpi_sum(dos_fermie, mpi_enreg%comm_kptband, ier)
2411 :
2412 : ! 2) Fermi-level variation
2413 :
2414 55 : delta_fermie = delta_occ_tot / dos_fermie
2415 :
2416 : !Loop over spins and kpoints to do delta_occ(eigenvalue) -= f'(eigenvalue) * delta_fermie
2417 165 : do isppol =1, dtset%nsppol
2418 385 : do ikpt = 1, dtset%nkpt
2419 220 : i_kpt_sppol = ikpt+(isppol-1)*dtset%nkpt
2420 220 : nband_k = dtset%nband(i_kpt_sppol)
2421 :
2422 220 : iband1 = needed_bands_bounds(1, i_kpt_sppol)
2423 220 : iband2 = needed_bands_bounds(2, i_kpt_sppol)
2424 1512 : do iband = iband1, iband2 ! Loop over needed bands
2425 :
2426 1182 : i_eigen = get_eigen_index(dtset, iband, ikpt, isppol)
2427 1182 : eigenval = this%eigen(i_eigen)
2428 1182 : fp = derivative_occ(dtset%occopt, eigenval, this%fermie, this%precon_tsmear) * maxocc
2429 : !fp = doccde(i_eigen)
2430 1402 : delta_occ(i_eigen) = delta_occ(i_eigen) - fp * delta_fermie
2431 :
2432 : end do
2433 : end do !ikpt
2434 : end do !isppol
2435 : !ABI_FREE(doccde)
2436 :
2437 55 : ABI_FREE(needed_bands_bounds)
2438 55 : ABI_FREE(needed_bands_number)
2439 :
2440 55 : end subroutine compute_delta_occ
2441 :
2442 : !****f* m_precon/compute_delta_rho_from_delta_occ_only
2443 : !! NAME
2444 : !! compute_delta_rho_from_delta_occ_only
2445 : !!
2446 : !! FUNCTION
2447 : !! Compute the first-order density variation induced by a set of occupation
2448 : !! variations delta_occ, neglecting wavefunction variations:
2449 : !! delta_rho(r) = sum_nk wtk(k) * delta_occ(nk) * |psi_nk(r)|^2
2450 : !! The result is symmetrized with symrhg and returned in the Pauli basis.
2451 : !!
2452 : !! INPUTS
2453 : !! dtset = All input variables for this dataset.
2454 : !! mpi_enreg = Information about MPI parallelization.
2455 : !! delta_occ(mband*nkpt*nsppol) = First-order variation of the occupations.
2456 : !!
2457 : !! OUTPUTS
2458 : !! delta_rho(nfftprc, nspden) = First-order density variation in the Pauli basis.
2459 : !!
2460 : !! SOURCE
2461 55 : subroutine compute_delta_rho_from_delta_occ_only(this, dtset, mpi_enreg, delta_occ, delta_rho)
2462 : !Arguments ------------------------------------
2463 : class(precon_object) :: this
2464 : !scalars
2465 : type(dataset_type),intent(in) :: dtset
2466 : type(MPI_type), intent(in) :: mpi_enreg
2467 : !arrays
2468 : real(dp), intent(in) :: delta_occ(size(this%eigen))
2469 : real(dp), intent(inout) :: delta_rho(this%nfftprc, dtset%nspden)
2470 :
2471 : !Local variables-------------------------------
2472 : integer :: iband, isppol, ispden, ikpt, i_eigen, nband_k
2473 : integer :: ier
2474 : integer :: maxocc
2475 : integer :: iband1, iband2, i_kpt_sppol
2476 : real(dp) :: fp
2477 : !arrays
2478 55 : integer, allocatable :: needed_bands_bounds(:, :)
2479 55 : integer, allocatable :: needed_bands_number(:)
2480 55 : real(dp), allocatable :: rhoi_r(:, :), delta_rho_g(:, :)
2481 :
2482 : ! *************************************************************************
2483 :
2484 165 : ABI_MALLOC(needed_bands_bounds, (2, dtset%nkpt*dtset%nsppol))
2485 165 : ABI_MALLOC(needed_bands_number, (dtset%nkpt*dtset%nsppol))
2486 55 : call get_needed_bands_chi0diag(this, dtset, mpi_enreg, needed_bands_bounds, needed_bands_number)
2487 :
2488 55 : maxocc = two / (dtset%nsppol * dtset%nspinor) !Maximum number of occupations (1 or 2)
2489 371415 : delta_rho = zero
2490 :
2491 : !Loop over spins and kpoints
2492 165 : do isppol =1, dtset%nsppol
2493 385 : do ikpt = 1, dtset%nkpt
2494 :
2495 220 : i_kpt_sppol = ikpt+(isppol-1)*dtset%nkpt
2496 220 : nband_k = dtset%nband(i_kpt_sppol)
2497 :
2498 : ! MPI parallelization over kpoints : cycle if kpt does not belong to current processor.
2499 220 : if (proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, 1, nband_k, isppol, mpi_enreg%me_kpt)) then
2500 : cycle
2501 : end if
2502 :
2503 220 : iband1 = needed_bands_bounds(1, i_kpt_sppol)
2504 220 : iband2 = needed_bands_bounds(2, i_kpt_sppol)
2505 :
2506 1512 : do iband = iband1, iband2 ! Loop over needed bands
2507 :
2508 1182 : if (cycle_band(mpi_enreg, nband_k, iband)) then
2509 : cycle
2510 : end if ! Check if the band belongs to the current processor.
2511 :
2512 :
2513 1182 : i_eigen = get_eigen_index(dtset, iband, ikpt, isppol) ! Index of (iband, ikpt, isppol) in eigen array.
2514 1182 : fp = derivative_occ(dtset%occopt, this%eigen(i_eigen), this%fermie, this%precon_tsmear) * maxocc
2515 :
2516 : ! Compute and add the contribution to delta_rhol.
2517 1402 : if (dtset%nspinor == 1) then
2518 : ! Collinear magnetism or no magnetism
2519 :
2520 1182 : ispden = isppol
2521 1182 : if (this%use_precomputed_rhoi) then
2522 : ! We use the precomputed orbital density.
2523 : delta_rho(:, ispden) = delta_rho(:, ispden) + &
2524 : & dtset%wtk(ikpt) * delta_occ(i_eigen) * &
2525 3990432 : & this%precomputed_rhoi(:, 1, this%precomputed_rhoi_indices(iband, ikpt, isppol))
2526 : else
2527 : ! We need to recompute the orbital density.
2528 0 : ABI_MALLOC(rhoi_r, (this%nfftprc, 1))
2529 0 : call compute_rhoi_coll(this, dtset, mpi_enreg, iband, ikpt, isppol, rhoi_r)
2530 : delta_rho(:, ispden) = delta_rho(:, ispden) + &
2531 0 : & dtset%wtk(ikpt) * delta_occ(i_eigen) * rhoi_r(:, 1)
2532 0 : ABI_FREE(rhoi_r)
2533 : end if
2534 : ! delta_rho in up/down representation = expected representation by symrhg
2535 :
2536 : ! TODO : here it is assumed that nsppol=2 => nspden=2 ...
2537 :
2538 0 : elseif (dtset%nspden == 4) then
2539 : ! Non collinear magnetism
2540 :
2541 0 : if (this%use_precomputed_rhoi) then
2542 : ! We use the precomputed orbital density.
2543 0 : do ispden = 1, 4
2544 : delta_rho(:, ispden) = delta_rho(:, ispden) + &
2545 : & dtset%wtk(ikpt) * delta_occ(i_eigen) * &
2546 0 : & this%precomputed_rhoi(:, ispden, this%precomputed_rhoi_indices(iband, ikpt, isppol))
2547 : end do
2548 : else
2549 : ! We need to recompute the orbital density.
2550 0 : ABI_MALLOC(rhoi_r, (this%nfftprc, 4))
2551 0 : call compute_rhoi_noncoll(this, dtset, mpi_enreg, iband, ikpt, isppol, rhoi_r)
2552 0 : do ispden = 1, 4
2553 : delta_rho(:, ispden) = delta_rho(:, ispden) + &
2554 0 : & dtset%wtk(ikpt) * delta_occ(i_eigen) * rhoi_r(:, ispden)
2555 : end do
2556 0 : ABI_FREE(rhoi_r)
2557 : end if
2558 : ! delta_rho in Pauli representation = NOT the expected representation by symrhg
2559 :
2560 : else
2561 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): TODO non-collinear magnetism with nspden=/4.")
2562 : end if
2563 :
2564 : end do ! iband
2565 :
2566 : end do ! ikpt
2567 : end do ! isppol
2568 :
2569 55 : ABI_FREE(needed_bands_bounds)
2570 55 : ABI_FREE(needed_bands_number)
2571 :
2572 : ! MPI parallelization over kpoints and bands : sum delta_rho on all processors.
2573 : ier = 0
2574 55 : call xmpi_sum(delta_rho, mpi_enreg%comm_kptband, ier)
2575 :
2576 : ! Symmetrization
2577 165 : ABI_MALLOC(delta_rho_g, (2, this%nfftprc))
2578 55 : if (dtset%nspinor == 1) then
2579 : ! In collinear magnetism, we can readily apply symrhg to delta_rho (correct spin representation).
2580 : call symrhg(1, this%gprimd, this%irrzon, mpi_enreg, dtset%nfft, dtset%nfft, dtset%ngfft, dtset%nspden, dtset%nsppol, &
2581 55 : & dtset%nsym, this%phnons, delta_rho_g, delta_rho, this%rprimd, dtset%symafm, dtset%symrel, dtset%tnons)
2582 : else
2583 : ! In non-collinear magnetism, we apply symrhg independantly to each spin component,
2584 : ! to avoid having to change the spin representation of delta_rho.
2585 0 : do ispden = 1, dtset%nspden
2586 : call symrhg(1, this%gprimd, this%irrzon, mpi_enreg, dtset%nfft, dtset%nfft, dtset%ngfft, 1, 1, &
2587 0 : & dtset%nsym, this%phnons, delta_rho_g, delta_rho(:, ispden:ispden), this%rprimd, dtset%symafm, dtset%symrel, dtset%tnons)
2588 : end do
2589 : end if
2590 55 : ABI_FREE(delta_rho_g)
2591 :
2592 55 : call to_pauli(1, delta_rho)
2593 :
2594 55 : end subroutine compute_delta_rho_from_delta_occ_only
2595 :
2596 : !****f* m_precon/precompute_psii
2597 : !! NAME
2598 : !! precompute_psii
2599 : !!
2600 : !! FUNCTION
2601 : !! NOT WORKING - WIP
2602 : !! Precompute and store the real-space wavefunctions psi_i(r) as well as the
2603 : !! corresponding orbital densities rho_i(r) = |psi_i(r)|^2 for all bands with
2604 : !! a non-negligible occupation derivative, to avoid redundant FFTs during the
2605 : !! iterative application of chi0_quasidiag.
2606 : !! Results are stored in this%precomputed_psii / this%precomputed_rhoi and
2607 : !! indexed by this%precomputed_psii_indices / this%precomputed_rhoi_indices.
2608 : !!
2609 : !! INPUTS
2610 : !! dtset = All input variables for this dataset.
2611 : !! mpi_enreg = Information about MPI parallelization.
2612 : !!
2613 : !! SOURCE
2614 0 : subroutine precompute_psii(this, dtset, mpi_enreg)
2615 :
2616 : !Arguments ------------------------------------
2617 : class(precon_object), intent(inout) :: this
2618 : !scalars
2619 : type(dataset_type),intent(in) :: dtset
2620 : type(MPI_type), intent(in) :: mpi_enreg
2621 :
2622 : !Local variables-------------------------------
2623 : !scalars
2624 : integer :: n1, n2, n3, n4, n5, n6
2625 : integer :: isppol, ikpt, i_kpt_sppol, nband_k, iband, iband1, iband2
2626 : integer :: i_psii, option, ndat, istwf_k, npw_k, idat, ispinor, icplex, tim_fourwf
2627 0 : integer :: i_kg(2), i_cg_iband1(2*dtset%nspinor), i_cg_iband2(2*dtset%nspinor)
2628 : real(dp) :: sum_rhoi_r
2629 : integer :: ifft
2630 : logical :: band_paral
2631 : !arrays
2632 0 : integer, allocatable :: needed_bands_bounds(:, :)
2633 0 : integer, allocatable :: needed_bands_number(:)
2634 0 : integer, allocatable :: kg_k(:, :)
2635 0 : integer :: gbound_k(2*dtset%mgfft+8,2)
2636 0 : real(dp), allocatable :: psii_aug(:, :, :, :)
2637 0 : real(dp), allocatable :: rhoi_aug(:, :, :)
2638 : !for band parall
2639 : integer :: option_fourwf, blocksize, iblock, ibandblock1, ibandblock2, nbdblock, nfft_blocks
2640 0 : integer :: i_cg_ibandblock1(2*dtset%nspinor), i_cg_ibandblock2(2*dtset%nspinor)
2641 0 : real(dp), allocatable :: dummy_occ_k(:)
2642 : !dummy arguments
2643 0 : real(dp) :: dummy_denpot(0, dtset%ngfft(5), dtset%ngfft(6)), dummy_fofgout(2, 0)
2644 :
2645 : ! *************************************************************************
2646 0 : ABI_BUG("WIP - precompute_psii")
2647 :
2648 0 : band_paral = (dtset%paral_kgb == 1 .and. dtset%npband > 1)
2649 :
2650 0 : n1 = dtset%ngfft(1)
2651 0 : n2 = dtset%ngfft(2)
2652 0 : n3 = dtset%ngfft(3)
2653 0 : n4 = dtset%ngfft(4)
2654 0 : n5 = dtset%ngfft(5)
2655 0 : n6 = dtset%ngfft(6)
2656 :
2657 0 : ABI_MALLOC(needed_bands_bounds, (2, dtset%nkpt*dtset%nsppol))
2658 0 : ABI_MALLOC(needed_bands_number, (dtset%nkpt*dtset%nsppol))
2659 0 : call get_needed_bands_chi0diag(this, dtset, mpi_enreg, needed_bands_bounds, needed_bands_number)
2660 :
2661 : ! Allocate the array containing the precomputed psii
2662 0 : ABI_MALLOC(this%precomputed_psii, (2, this%nfftprc, dtset%nspinor, sum(needed_bands_number)))
2663 0 : this%precomputed_psii_indices = zero
2664 0 : ABI_MALLOC(this%precomputed_rhoi, (this%nfftprc, dtset%nspinor, sum(needed_bands_number)))
2665 0 : this%precomputed_rhoi_indices = zero
2666 0 : i_psii = 1
2667 :
2668 0 : ABI_MALLOC(psii_aug, (2, n4, n5, n6*dtset%mband))
2669 0 : ABI_MALLOC(rhoi_aug, (n4, n5, n6))
2670 :
2671 : !Loop over spins and kpoints
2672 0 : do isppol = 1, dtset%nsppol
2673 0 : do ikpt = 1, dtset%nkpt
2674 0 : i_kpt_sppol = ikpt+(isppol-1)*dtset%nkpt
2675 :
2676 : ! MPI parallelization over kpoints : cycle if kpt does not belong to current processor.
2677 0 : nband_k = dtset%nband(i_kpt_sppol)
2678 0 : if (proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, 1, nband_k, isppol, mpi_enreg%me_kpt)) then
2679 : cycle
2680 : end if
2681 :
2682 0 : iband1 = needed_bands_bounds(1, i_kpt_sppol)
2683 0 : iband2 = needed_bands_bounds(2, i_kpt_sppol)
2684 :
2685 0 : if (band_paral) then
2686 :
2687 0 : if (dtset%nspinor==1) then
2688 0 : option_fourwf = 0
2689 :
2690 0 : nbdblock = nband_k / (mpi_enreg%nproc_band * mpi_enreg%bandpp)
2691 0 : blocksize = nband_k / nbdblock
2692 0 : ndat = mpi_enreg%bandpp
2693 :
2694 0 : do iblock = 1, nbdblock ! Loop over (LOBPCG) blocks
2695 :
2696 0 : ibandblock1 = blocksize*(iblock-1) + 1
2697 0 : ibandblock2 = blocksize*(iblock)
2698 0 : i_cg_ibandblock1 = this%cg_indices(:, ibandblock1, ikpt, isppol)
2699 0 : i_cg_ibandblock2 = this%cg_indices(:, ibandblock2, ikpt, isppol) ! Changer
2700 :
2701 0 : ABI_MALLOC(psii_aug, (2, n4, n5, n6*ndat))
2702 0 : ABI_MALLOC(dummy_occ_k, (nband_k))
2703 : !ABI_MALLOC(dummy_denpot, (n4, n5, n6))
2704 :
2705 0 : call bandfft_kpt_set_ikpt(ikpt, mpi_enreg)
2706 0 : nfft_blocks = 1 ! TODO : what is that ??
2707 : call prep_fourwf(dummy_denpot, blocksize, this%cg(:, i_cg_ibandblock1(1):i_cg_ibandblock2(2)), &
2708 : & psii_aug(:, :, :, 1:n6*ndat), iblock, dtset%istwfk(ikpt), dtset%mgfft, mpi_enreg, nband_k, &
2709 : & ndat, dtset%ngfft, this%npwarr(ikpt), &
2710 0 : & n4, n5, n6, dummy_occ_k, option_fourwf, this%ucvol, dtset%wtk(ikpt), nfft_blocks)
2711 :
2712 0 : ABI_FREE(dummy_occ_k)
2713 : !ABI_FREE(dummy_denpot)
2714 :
2715 : end do
2716 : else
2717 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): TODO non-collinear magnetisme in precon with band paral")
2718 : end if
2719 :
2720 : else
2721 :
2722 0 : if (dtset%nspinor==1) then
2723 :
2724 0 : option = 0
2725 0 : ndat = needed_bands_number(i_kpt_sppol)
2726 0 : istwf_k = dtset%istwfk(ikpt) ! Option parameter that describes the storage of wfs at this kpt.
2727 0 : npw_k = this%npwarr(ikpt) ! Number of plane-wave at this kpt.
2728 0 : ABI_MALLOC(kg_k, (3, npw_k))
2729 0 : i_kg = this%kg_indices(:, ikpt)
2730 0 : kg_k = this%kg(:, i_kg(1):i_kg(2))
2731 0 : call sphereboundary(gbound_k, istwf_k, kg_k, dtset%mgfft, npw_k) ! Computes gbound.
2732 0 : tim_fourwf = 0
2733 0 : i_cg_iband1 = this%cg_indices(:, iband1, ikpt, isppol)
2734 0 : i_cg_iband2 = this%cg_indices(:, iband2, ikpt, isppol)
2735 :
2736 : call fourwf(0, dummy_denpot, this%cg(:, i_cg_iband1(1):i_cg_iband2(2)), dummy_fofgout, &
2737 : & psii_aug(:, :, :, 1:n6*ndat), gbound_k, gbound_k, istwf_k, kg_k, kg_k, &
2738 : & dtset%mgfft, mpi_enreg, ndat, dtset%ngfft, npw_k, npw_k, &
2739 0 : & n4, n5, n6, option, tim_fourwf, one, one)
2740 :
2741 0 : ABI_FREE(kg_k)
2742 :
2743 : else
2744 0 : ABI_BUG("chi0-based preconditioner (iprcel=2**): TODO non-collinear magnetisme in precon")
2745 : end if
2746 : end if
2747 :
2748 0 : do iband = 1, nband_k
2749 :
2750 0 : if (cycle_band(mpi_enreg, nband_k, iband)) then
2751 : cycle
2752 : end if ! Checks if the band belongs to the current processor.
2753 :
2754 0 : idat = idat + 1
2755 :
2756 : ! Check if this band is needed
2757 0 : if (.not. (iband1 <= iband .and. iband <= iband2)) then
2758 : cycle
2759 : end if
2760 :
2761 0 : rhoi_aug = zero
2762 0 : call cg_addtorho(n1, n2, n3, n4, n5, n6, 1, one, one, psii_aug(:, :, :, (idat-1)*n6+1:idat*n6), rhoi_aug)
2763 :
2764 : ! Grid transfers:
2765 0 : ispinor = 1
2766 0 : do icplex = 1, 2
2767 0 : call transfer_grid(this, dtset, mpi_enreg, 1, psii_aug(icplex, :, :, idat:idat+n6-1), this%precomputed_psii(icplex, :, ispinor:ispinor, i_psii))
2768 : ! TODO : this will create a copy ... Maybe change dim order ?
2769 : end do
2770 0 : call transfer_grid(this, dtset, mpi_enreg, 1, rhoi_aug, this%precomputed_rhoi(:, :, i_psii))
2771 :
2772 : ! Normalize psii:
2773 0 : sum_rhoi_r = 0
2774 0 : do ispinor = 1, dtset%nspinor
2775 0 : do ifft=1, this%nfftprc
2776 0 : sum_rhoi_r = sum_rhoi_r + this%precomputed_rhoi(ifft, ispinor, i_psii)
2777 : end do
2778 0 : this%precomputed_psii(:, :, ispinor, i_psii) = this%precomputed_psii(:, :, ispinor, i_psii) / sqrt(sum_rhoi_r * this%dvol)
2779 0 : this%precomputed_rhoi(:, ispinor, i_psii) = this%precomputed_rhoi(:, ispinor, i_psii) / (sum_rhoi_r * this%dvol)
2780 : end do
2781 :
2782 : ! Save the index for iband, ikpt, isppol in precomputed_psii and precomputed_rhoi
2783 0 : this%precomputed_psii_indices(iband, ikpt, isppol) = i_psii
2784 0 : this%precomputed_rhoi_indices(iband, ikpt, isppol) = i_psii
2785 0 : i_psii = i_psii + 1
2786 :
2787 : end do !iband
2788 :
2789 :
2790 : end do !ikpt
2791 : end do !isppol
2792 :
2793 0 : ABI_FREE(rhoi_aug)
2794 0 : ABI_FREE(psii_aug)
2795 0 : ABI_FREE(needed_bands_number)
2796 0 : ABI_FREE(needed_bands_bounds)
2797 :
2798 0 : end subroutine precompute_psii
2799 :
2800 : !****f* m_precon/compute_delta_wf
2801 : !! NAME
2802 : !! compute_delta_wf
2803 : !!
2804 : !! FUNCTION
2805 : !! NOT WORKING - WIP
2806 : !!
2807 : !! INPUTS
2808 : !! dtset = All input variables for this dataset.
2809 : !! mpi_enreg = Information about MPI parallelization.
2810 : !! delta_V =
2811 : !!
2812 : !! OUTPUTS
2813 : !! delta_wf =
2814 : !!
2815 : !! SOURCE
2816 : subroutine compute_delta_wf(this, dtset, mpi_enreg, delta_V, delta_wf)
2817 : !Arguments ------------------------------------
2818 : class(precon_object), intent(inout) :: this
2819 : type(dataset_type), intent(in) :: dtset
2820 : type(MPI_type), intent(in) :: mpi_enreg
2821 : real(dp), intent(inout) :: delta_V(this%nfftprc, dtset%nspden)
2822 : real(dp), intent(inout) :: delta_wf(:, :)
2823 :
2824 : !Local variables-------------------------------
2825 : !scalars
2826 : integer :: n1, n2, n3, n4, n5, n6
2827 : integer :: isppol, ispden, ikpt, nband_k, npw_k, istwf_k, iband, jband
2828 : integer :: i_eigen, j_eigen, ier, ndat, option
2829 : integer :: i_kg(2), j_cg(2), i_cg(2)
2830 : real(dp) :: fi, fj, ddiff, coeff
2831 : integer :: tim_fourwf
2832 : integer :: gbound(2*dtset%mgfft+8,2)
2833 : real(dp) :: dotr, doti
2834 : !arrays
2835 : real(dp), allocatable :: delta_V_wf_i_r(:, :,:,:)
2836 : real(dp), allocatable :: delta_V_wf_i(:,:)
2837 : real(dp), allocatable :: delta_V_aug(:,:,:)
2838 : !dummy
2839 : integer :: dummy_int
2840 : real(dp) :: dummy_real
2841 : real(dp) :: dummy_denpot(0, dtset%ngfft(5), dtset%ngfft(6)), dummy_fofg(2, 0)
2842 :
2843 : ! ***************************************************************************
2844 : ABI_BUG("WIP - compute_delta_wf")
2845 :
2846 : n1 = dtset%ngfft(1)
2847 : n2 = dtset%ngfft(2)
2848 : n3 = dtset%ngfft(3)
2849 : n4 = dtset%ngfft(4)
2850 : n5 = dtset%ngfft(5)
2851 : n6 = dtset%ngfft(6)
2852 :
2853 : delta_wf = zero
2854 :
2855 : ABI_MALLOC(delta_V_wf_i_r, (2, n4, n5, n6))
2856 : ABI_MALLOC(delta_V_wf_i, (2, dtset%mpw))
2857 : ABI_MALLOC(delta_V_aug, (n4, n5, n6))
2858 :
2859 : !Loop over spins and kpoints
2860 : do isppol = 1, dtset%nsppol
2861 :
2862 : ispden = isppol !TODO noncoll
2863 : ! Transfer delta_V to the augmented (wavefunction) fft-grid
2864 : call fftpac(ispden, mpi_enreg, dtset%nspden, n1, n2, n3, n4, n5, n6, dtset%ngfft, delta_V, delta_V_aug, 2)
2865 :
2866 : do ikpt = 1, dtset%nkpt
2867 : nband_k = dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
2868 : !MPI parallelization over kpoints : cycle if kpt does not belong to current processor.
2869 : if (proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, 1, nband_k, isppol, mpi_enreg%me_kpt)) then
2870 : cycle
2871 : end if
2872 : npw_k = this%npwarr(ikpt) ! Number of plane-wave at this kpt.
2873 : istwf_k = dtset%istwfk(ikpt) ! Option parameter that describes the storage of wfs at this kpt.
2874 : i_kg = this%kg_indices(:, ikpt)
2875 : call sphereboundary(gbound, istwf_k, this%kg(:, i_kg(1):i_kg(2)), dtset%mgfft, npw_k) ! Computes gbound.
2876 :
2877 : do iband = 1, nband_k
2878 :
2879 : i_eigen = get_eigen_index(dtset, iband, ikpt, isppol) ! Index of (iband, ikpt, isppol) in eigen array.
2880 : if (abs(this%eigen(i_eigen) - this%fermie) > tol1) then ! TODO : condition as input
2881 : cycle
2882 : end if
2883 : fi = this%occ(i_eigen)
2884 :
2885 : i_cg = this%cg_indices(:, iband, ikpt, isppol)
2886 :
2887 : !Input parameters for fourwf :
2888 : ndat = 1 ! Only one FFT.
2889 : tim_fourwf = 0
2890 : ! Multiply the wave function by the potential variation delta_V (ifft -> real space multiplication -> fft)
2891 : option = 0 ! ifft
2892 : call fourwf(dummy_int, dummy_denpot, this%cg(:, i_cg(1):i_cg(2)), dummy_fofg, delta_V_wf_i_r, &
2893 : & gbound, gbound, istwf_k, this%kg(:, i_kg(1):i_kg(2)), this%kg(:, i_kg(1):i_kg(2)), &
2894 : & dtset%mgfft, mpi_enreg, ndat, dtset%ngfft, npw_k, &
2895 : & dummy_int, n4, n5, n6, option, tim_fourwf, dummy_real, dummy_real)
2896 : call cg_vlocpsi(n4, n5, n6, n4, n5, n6, 1, 1, delta_V_aug, delta_V_wf_i_r)
2897 : option = 3 ! fft
2898 : call fourwf(dummy_int, dummy_denpot, dummy_fofg, delta_V_wf_i(:, 1:npw_k), delta_V_wf_i_r, &
2899 : & gbound, gbound, istwf_k, this%kg(:, i_kg(1):i_kg(2)), this%kg(:, i_kg(1):i_kg(2)), &
2900 : & dtset%mgfft, mpi_enreg, ndat, dtset%ngfft, npw_k, &
2901 : & dummy_int, n4, n5, n6, option, tim_fourwf, dummy_real, dummy_real)
2902 : !write(6,*)'chi0diel compute_delta_wf - ok2 '; flush(6) !DEBUG
2903 : ! result stored in delta_V_wf_i(:, i_cg(1):i_cg(2))
2904 :
2905 : do jband = 1, nband_k
2906 : !write(6,*)'chi0diel compute_delta_wf - if yes '; flush(6) !DEBUG
2907 : j_eigen = get_eigen_index(dtset, jband, ikpt, isppol) ! Index of (jband, ikpt, isppol) in eigen array.
2908 :
2909 : if (abs(this%eigen(j_eigen) - this%fermie) > tol1) then ! TODO : condition as input
2910 : cycle
2911 : end if
2912 :
2913 : if (iband == jband) then
2914 : cycle
2915 : end if ! i=j contribution computed in compute_delta_occ
2916 :
2917 : j_cg = this%cg_indices(:, jband, ikpt, isppol)
2918 : fj = this%occ(j_eigen)
2919 :
2920 : ! Compute an equivalent of coeff=1/(eigen_i - eigen_j) that ensure correct compensation of the terms in compute_delta_occ
2921 : if (abs(this%eigen(i_eigen) - this%eigen(j_eigen)) < dtset%tsmear * tol10) then
2922 : ddiff = derivative_occ(dtset%occopt, (this%eigen(i_eigen)+this%eigen(j_eigen))/2, this%fermie, dtset%tsmear)
2923 : else
2924 : ddiff = (fi - fj)/(this%eigen(i_eigen) - this%eigen(j_eigen))
2925 : end if
2926 : coeff = ddiff * fi/(fi**2+fj**2) ! From DFTK
2927 : !write(6,*)'chi0diel compute_delta_wf - ok3 '; flush(6) !DEBUG
2928 :
2929 : ! Compute dot product between wavefunction (j) and delta_V
2930 : call dotprod_g(dotr, doti, istwf_k, npw_k, 2, this%cg(:, j_cg(1):j_cg(2)), delta_V_wf_i(:, 1:npw_k), 0, mpi_enreg%comm_spinorfft) ! TODO : check dotprof(psi_i, psi_i) = 1
2931 : !write(6,*)'chi0diel compute_delta_wf - ok4 '; flush(6) !DEBUG
2932 :
2933 : delta_wf(1, i_cg(1):i_cg(2)) = delta_wf(1, i_cg(1):i_cg(2)) + &
2934 : & coeff * ( dotr * this%cg(1, j_cg(1):j_cg(2)) - doti * this%cg(2, j_cg(1):j_cg(2)) )
2935 : delta_wf(2, i_cg(1):i_cg(2)) = delta_wf(2, i_cg(1):i_cg(2)) + &
2936 : & coeff * ( dotr * this%cg(2, j_cg(1):j_cg(2)) + doti * this%cg(1, j_cg(1):j_cg(2)) )
2937 : !write(6,*)'chi0diel compute_delta_wf - ok5 '; flush(6) !DEBUG
2938 :
2939 : end do
2940 : end do
2941 : end do
2942 : end do
2943 :
2944 : ABI_FREE(delta_V_wf_i_r)
2945 : ABI_FREE(delta_V_wf_i)
2946 : ABI_FREE(delta_V_aug)
2947 :
2948 : !MPI parallelization over kpoints : sum delta_occ on all processors.
2949 : ier = 0
2950 : call xmpi_sum(delta_wf, mpi_enreg%comm_kpt, ier)
2951 :
2952 : end subroutine compute_delta_wf
2953 :
2954 : !****f* m_precon/compute_delta_rho
2955 : !! NAME
2956 : !! compute_delta_rho
2957 : !!
2958 : !! FUNCTION
2959 : !! NOT WORKING - WIP
2960 : !! Compute the first-order density variation induced by both occupation
2961 : !! variations (delta_occ) and wavefunction variations (delta_wf):
2962 : !! delta_rho(r) = sum_nk wtk(k) * [ delta_occ(nk) * |psi_nk(r)|^2
2963 : !! + 2 * occ(nk) * Re(delta_psi_nk(r)* . psi_nk(r)) ]
2964 : !! The result is symmetrized with symrhg and returned in the Pauli basis.
2965 : !!
2966 : !! INPUTS
2967 : !! dtset = All input variables for this dataset.
2968 : !! mpi_enreg = Information about MPI parallelization.
2969 : !! delta_occ(mband*nkpt*nsppol) = First-order variation of the occupations.
2970 : !! delta_wf(2, mcg) = First-order variation of the wavefunctions in G-space.
2971 : !!
2972 : !! OUTPUTS
2973 : !! delta_rho(nfftprc, nspden) = First-order density variation in the Pauli basis.
2974 : !!
2975 : !! SOURCE
2976 0 : subroutine compute_delta_rho(this, dtset, mpi_enreg, delta_occ, delta_wf, delta_rho)
2977 : !Arguments ------------------------------------
2978 : class(precon_object) :: this
2979 : !scalars
2980 : type(dataset_type),intent(in) :: dtset
2981 : type(MPI_type), intent(in) :: mpi_enreg
2982 : !arrays
2983 : real(dp), intent(in) :: delta_occ(size(this%eigen))
2984 : real(dp), intent(inout) :: delta_wf(2, size(this%cg, 2))
2985 : real(dp), intent(inout) :: delta_rho(this%nfftprc, dtset%nspden)
2986 :
2987 : !Local variables-------------------------------
2988 : integer :: iband, isppol, ispden, ikpt, istwf_k, i_eigen, nband_k, npw_k
2989 : integer :: ier
2990 : integer :: n1, n2, n3, n4, n5, n6
2991 : integer :: tim_fourwf, ndat, option
2992 : integer :: maxocc
2993 : real(dp) :: fp
2994 : !arrays
2995 : integer :: i_cg(2), i_kg(2)
2996 0 : integer :: gbound(2*dtset%mgfft+8,2)
2997 0 : integer, allocatable :: kg_k(:, :)
2998 0 : real(dp), allocatable :: rho_aug_r_i(:, :, :), wf_aug_r_i(:, :, :, :), delta_wf_aug_r_i(:, :, :, :), delta_rho_aug_r(:, :, :, :)
2999 0 : real(dp), allocatable :: delta_rho_g(:, :)
3000 : !dummy arguments
3001 : integer :: dummy_int
3002 : real(dp) :: dummy_real
3003 0 : real(dp) :: dummy_denpot(0, dtset%ngfft(5), dtset%ngfft(6)), dummy_fofgout(2, 0)
3004 :
3005 : ! *************************************************************************
3006 0 : ABI_BUG("WIP - compute_delta_rho")
3007 :
3008 0 : n1 = dtset%ngfft(1)
3009 0 : n2 = dtset%ngfft(2)
3010 0 : n3 = dtset%ngfft(3)
3011 0 : n4 = dtset%ngfft(4)
3012 0 : n5 = dtset%ngfft(5)
3013 0 : n6 = dtset%ngfft(6)
3014 0 : maxocc = two / (dtset%nsppol * dtset%nspinor) !Maximum number of occupations (1 or 2)
3015 :
3016 0 : ABI_MALLOC(delta_rho_aug_r, (n4, n5, n6, dtset%nspden))
3017 0 : delta_rho_aug_r = zero
3018 0 : ABI_MALLOC(rho_aug_r_i, (n4, n5, n6))
3019 0 : ABI_MALLOC(wf_aug_r_i, (2, n4, n5, n6))
3020 0 : ABI_MALLOC(delta_wf_aug_r_i, (2, n4, n5, n6))
3021 :
3022 : !Loop over spins and kpoints
3023 0 : do isppol =1, dtset%nsppol
3024 0 : do ikpt = 1, dtset%nkpt
3025 :
3026 0 : nband_k = dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
3027 :
3028 : !MPI parallelization over kpoints : cycle if kpt does not belong to current processor.
3029 0 : if (proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, 1, nband_k, isppol, mpi_enreg%me_kpt)) then
3030 : cycle
3031 : end if
3032 :
3033 0 : istwf_k = dtset%istwfk(ikpt) ! Option parameter that describes the storage of wfs at this kpt.
3034 0 : if (.not. (dtset%paral_kgb == 1 .and. dtset%npband > 1)) then
3035 0 : npw_k = this%npwarr(ikpt) ! Number of plane-wave at this kpt.
3036 0 : ABI_MALLOC(kg_k, (3, npw_k))
3037 0 : i_kg = this%kg_indices(:, ikpt)
3038 0 : kg_k = this%kg(:, i_kg(1):i_kg(2))
3039 0 : call sphereboundary(gbound, istwf_k, kg_k, dtset%mgfft, npw_k) ! Computes gbound.
3040 : else
3041 0 : npw_k = bandfft_kpt(ikpt)%npw_tot ! Number of plane-wave (after transpose) at this kpt.
3042 0 : ABI_MALLOC(kg_k, (3, npw_k))
3043 0 : kg_k = bandfft_kpt(ikpt)%kg_k_gather ! Reduced plane-wave coordinate (k+G) of this kpt (after transpose).
3044 0 : call sphereboundary(gbound, istwf_k, kg_k, dtset%mgfft, npw_k)
3045 : end if ! TODO : Not very nice to have a if here ....
3046 :
3047 0 : do iband = 1, nband_k
3048 :
3049 0 : if (cycle_band(mpi_enreg, nband_k, iband)) then
3050 : cycle
3051 : end if ! Check if the band belongs to the current processor.
3052 :
3053 : !Indices
3054 0 : i_eigen = get_eigen_index(dtset, iband, ikpt, isppol) ! Index of (iband, ikpt, isppol) in eigen array.
3055 0 : fp = derivative_occ(dtset%occopt, this%eigen(i_eigen), this%fermie, this%precon_tsmear) * maxocc
3056 :
3057 0 : if (abs(fp) > this%deigvals_tol_fp) then
3058 :
3059 : !Input parameters for fourwf :
3060 0 : ndat = 1 ! Only one FFT.
3061 0 : tim_fourwf = 0
3062 0 : rho_aug_r_i = zero
3063 0 : option = 1
3064 :
3065 : !1) Contribution of the occupation variation
3066 : ! IFFT for wavefunction
3067 0 : if (.not. (dtset%paral_kgb == 1 .and. dtset%npband > 1)) then
3068 0 : i_cg = this%cg_indices(:, iband, ikpt, isppol) ! Indices range of (iband, ikpt, isppol) in cg array
3069 : call fourwf(1, rho_aug_r_i, this%cg(:, i_cg(1):i_cg(2)), dummy_fofgout, wf_aug_r_i, &
3070 : & gbound, gbound, istwf_k, kg_k, kg_k, &
3071 : & dtset%mgfft, mpi_enreg, ndat, dtset%ngfft, npw_k, &
3072 0 : & dummy_int, n4, n5, n6, option, tim_fourwf, one, one)
3073 : else
3074 : !i_cg = this%cg_fft_indices(:, iband, ikpt, isppol) ! Indices range of (iband, ikpt, isppol) in cg array
3075 : !call fourwf(1, rho_aug_r_i, this%cg_fft(:, i_cg(1):i_cg(2)), dummy_fofgout, wf_aug_r_i, &
3076 : !& gbound, gbound, istwf_k, kg_k, kg_k, &
3077 : !& dtset%mgfft, mpi_enreg, ndat, dtset%ngfft, npw_k, &
3078 : !& dummy_int, n4, n5, n6, option, tim_fourwf, one, one)
3079 : end if ! TODO : change this to remove if ! (use pointers ?)
3080 :
3081 : ! Compute and add the contribution to delta_rho
3082 0 : if (dtset%nspinor == 1) then
3083 0 : ispden = isppol
3084 : delta_rho_aug_r(:, :, :, ispden) = delta_rho_aug_r(:, :, :, ispden) + &
3085 0 : & dtset%wtk(ikpt)/this%ucvol * delta_occ(i_eigen) * rho_aug_r_i
3086 : else
3087 0 : ABI_BUG("Non collinear not implemented")
3088 : end if
3089 :
3090 : !2) Contribution of wavefunction variation (if applicable)
3091 : if (.false.) then
3092 : !if (norm2(delta_wf(:, i_cg(1):i_cg(2))) > tol14) then ! TODO : what tol ?
3093 : ! TODO : his part is probably incorrect !
3094 : ! IFFT for delta_wf (wavefunction variation)
3095 : call fourwf(1, dummy_denpot, delta_wf(:, i_cg(1):i_cg(2)), dummy_fofgout, delta_wf_aug_r_i, &
3096 : & gbound, gbound, istwf_k, this%kg(:, i_kg(1):i_kg(2)), this%kg(:, i_kg(1):i_kg(2)), &
3097 : & dtset%mgfft, mpi_enreg, ndat, dtset%ngfft, npw_k, &
3098 : & dummy_int, n4, n5, n6, 0, tim_fourwf, dummy_real, dummy_real)
3099 :
3100 : ! Sum over MPI processes (for parallelization over bands).
3101 : ier = 0
3102 : call xmpi_sum(delta_wf_aug_r_i, mpi_enreg%comm_bandfft, ier)
3103 :
3104 : ! Compute and add the contribution to delta_rho
3105 : if (dtset%nspinor == 1) then
3106 : ispden = isppol
3107 : delta_rho_aug_r(:, :, :, ispden) = delta_rho_aug_r(:, :, :, ispden) + &
3108 : & dtset%wtk(ikpt) * this%occ(i_eigen) * &
3109 : & ((delta_wf_aug_r_i(1, :, : , :)*wf_aug_r_i(1, :, :, :)) + &
3110 : & (delta_wf_aug_r_i(2, :, : , :)*wf_aug_r_i(2, :, :, :))) ! TODO : false ?
3111 : else
3112 : ABI_BUG("Non collinear not implemented")
3113 : end if
3114 : end if
3115 : end if
3116 :
3117 : end do ! iband
3118 :
3119 0 : ABI_FREE(kg_k)
3120 : end do ! ikpt
3121 : end do ! isppol
3122 :
3123 0 : ABI_FREE(rho_aug_r_i)
3124 0 : ABI_FREE(wf_aug_r_i)
3125 0 : ABI_FREE(delta_wf_aug_r_i)
3126 :
3127 0 : do ispden=1, dtset%nspden
3128 0 : call transfer_grid(this, dtset, mpi_enreg, ispden, delta_rho_aug_r(:, :, :, ispden), delta_rho)
3129 : end do
3130 :
3131 0 : ABI_FREE(delta_rho_aug_r)
3132 :
3133 : ! MPI parallelization over kpoints and bands : sum delta_rho on all processors.
3134 : ier = 0
3135 0 : call xmpi_sum(delta_rho, mpi_enreg%comm_kptband, ier)
3136 :
3137 0 : ABI_MALLOC(delta_rho_g, (2, this%nfftprc))
3138 : call symrhg(1, this%gprimd, this%irrzon, mpi_enreg, dtset%nfft, dtset%nfft, dtset%ngfft, dtset%nspden, dtset%nsppol, &
3139 0 : & dtset%nsym, this%phnons, delta_rho_g, delta_rho, this%rprimd, dtset%symafm, dtset%symrel, dtset%tnons) ! TODO : carefull with tnons and symrel (do they get updated during structure optimization ?)
3140 0 : ABI_FREE(delta_rho_g)
3141 : ! TODO : deal with symmetries when spin (nsppol in non coll)
3142 :
3143 0 : call to_pauli(1, delta_rho)
3144 :
3145 0 : end subroutine compute_delta_rho
3146 :
3147 : !****f* m_precon/apply_chi0_diag
3148 : !! NAME
3149 : !! apply_chi0_diag
3150 : !!
3151 : !! FUNCTION
3152 : !! Apply the diagonal approximation of the independent-particle susceptibility
3153 : !! chi0 to a vector vec_r (in place). The diagonal approximation retains only
3154 : !! the intraband (n=n') terms:
3155 : !! (chi0_diag * v)(r) = sum_nk wtk(k) * f'(e_nk - e_F) * <psi_nk|v|psi_nk> * |psi_nk(r)|^2
3156 : !! with a Fermi-level correction to conserve the electron number.
3157 : !!
3158 : !! INPUTS
3159 : !! dtset = All input variables for this dataset.
3160 : !! mpi_enreg = Information about MPI parallelization.
3161 : !!
3162 : !! SIDE EFFECTS
3163 : !! vec_r(nfftprc, nspden) = On input: potential vector in the Pauli basis.
3164 : !! On output: chi0_diag * vec_r in the Pauli basis (density).
3165 : !!
3166 : !! SOURCE
3167 55 : subroutine apply_chi0_diag(this, dtset, mpi_enreg, vec_r)
3168 :
3169 : !Arguments ------------------------------------
3170 : class(precon_object) :: this
3171 : !scalars
3172 : type(dataset_type),intent(in) :: dtset
3173 : type(MPI_type), intent(in) :: mpi_enreg
3174 : !arrays
3175 : real(dp), intent(inout) :: vec_r(this%nfftprc, dtset%nspden)
3176 :
3177 : !Local variables-------------------------------
3178 : !arrays
3179 55 : real(dp), allocatable :: delta_occ(:)
3180 :
3181 : ! *************************************************************************
3182 :
3183 165 : ABI_MALLOC(delta_occ, (size(this%eigen)))
3184 55 : call compute_delta_occ(this, dtset, mpi_enreg, vec_r, delta_occ)
3185 55 : call compute_delta_rho_from_delta_occ_only(this, dtset, mpi_enreg, delta_occ, vec_r)
3186 :
3187 55 : ABI_FREE(delta_occ)
3188 :
3189 55 : end subroutine apply_chi0_diag
3190 :
3191 : !****f* m_precon/apply_chi0_quasidiag
3192 : !! NAME
3193 : !! apply_chi0_quasidiag
3194 : !!
3195 : !! FUNCTION
3196 : !! NOT WORKING - WIP
3197 : !! Apply the quasidiagonal approximation of chi0 to a vector vec_r (in place).
3198 : !! Like chi0_diag but also include some wavefunction variation contributions.
3199 : !!
3200 : !! INPUTS
3201 : !! dtset = All input variables for this dataset.
3202 : !! mpi_enreg = Information about MPI parallelization.
3203 : !!
3204 : !! SIDE EFFECTS
3205 : !! vec_r(nfftprc, nspden) = On input: potential vector in the Pauli basis.
3206 : !! On output: chi0_quasidiag * vec_r in the Pauli basis.
3207 : !!
3208 : !! SOURCE
3209 0 : subroutine apply_chi0_quasidiag(this, dtset, mpi_enreg, vec_r)
3210 :
3211 : !Arguments ------------------------------------
3212 : class(precon_object) :: this
3213 : !scalars
3214 : type(dataset_type),intent(in) :: dtset
3215 : type(MPI_type), intent(in) :: mpi_enreg
3216 : !arrays
3217 : real(dp), intent(inout) :: vec_r(this%nfftprc, dtset%nspden)
3218 :
3219 : !Local variables-------------------------------
3220 : !arrays
3221 0 : real(dp), allocatable :: delta_occ(:), delta_wf(:, :)
3222 :
3223 : ! *************************************************************************
3224 :
3225 0 : ABI_MALLOC(delta_occ, (size(this%eigen)))
3226 0 : ABI_MALLOC(delta_wf, (2, size(this%cg, 2)))
3227 0 : call compute_delta_occ(this, dtset, mpi_enreg, vec_r, delta_occ)
3228 0 : delta_wf = zero
3229 : !call compute_delta_wf(this, dtset, mpi_enreg, vec_r, delta_wf)
3230 0 : vec_r = zero
3231 0 : call compute_delta_rho(this, dtset, mpi_enreg, delta_occ, delta_wf, vec_r)
3232 :
3233 0 : ABI_FREE(delta_occ)
3234 0 : ABI_FREE(delta_wf)
3235 :
3236 0 : end subroutine apply_chi0_quasidiag
3237 :
3238 : !****f* m_precon/apply_chi0
3239 : !! NAME
3240 : !! apply_chi0
3241 : !!
3242 : !! FUNCTION
3243 : !! NOT WORKING -WIP
3244 : !! Apply the model (determined by iprcel) chi0 operator to the vector vec_r (in place).
3245 : !!
3246 : !! INPUTS
3247 : !! dtset = All input variables for this dataset.
3248 : !! mpi_enreg = Information about MPI parallelization.
3249 : !! ispden = Index of spin-density component.
3250 : !!
3251 : !! SIDE EFFECTS
3252 : !! vec_r (nfftprc, nspden) = Vector (in direct space) to which the model chi0 operator is applied (in place).
3253 : !! When nspden > 1 vec_r is in the Pauli basis.
3254 : !!
3255 : !! SOURCE
3256 0 : subroutine apply_chi0(this, dtset, vec_r)
3257 :
3258 : !Arguments ------------------------------------
3259 : class(precon_object) :: this
3260 : !scalars
3261 : type(dataset_type),intent(in) :: dtset
3262 : !arrays
3263 : real(dp), intent(inout) :: vec_r(this%nfftprc, dtset%nspden)
3264 :
3265 : !Local variables-------------------------------
3266 : !scalars
3267 : integer :: ispden
3268 :
3269 : ! *************************************************************************
3270 :
3271 : !Kerker with user_defined parameter dielng
3272 0 : if (this%iprcel == 210) then
3273 0 : ispden = 1
3274 0 : vec_r(:, ispden) = (-1/(4*pi*(this%dielng)**2)) * vec_r(:, ispden)
3275 0 : do ispden = 2, dtset%nspden
3276 0 : vec_r(:, ispden) = 0
3277 : end do
3278 : end if
3279 :
3280 : !Kerker based on the DOS.
3281 0 : if (this%iprcel == 211) then
3282 0 : ispden = 1
3283 0 : vec_r(:, ispden) = -this%dos(1) * vec_r(:, ispden)
3284 0 : do ispden = 2, dtset%nspden
3285 0 : vec_r(:, ispden) = 0
3286 : end do
3287 : end if
3288 :
3289 : !Kerker based on the DOS, that can be different between the up and down chanels in spin-polarized cases.
3290 0 : if (this%iprcel == 201) then
3291 0 : vec_r(:, 1) = 0
3292 0 : do ispden = 1, dtset%nspden
3293 0 : vec_r(:, 1) = vec_r(:, 1) - this%dos(ispden) * vec_r(:, ispden)
3294 : end do
3295 0 : do ispden = 2, dtset%nspden
3296 0 : vec_r(:, ispden) = - this%dos(ispden) * vec_r(:, 1)
3297 : end do
3298 : end if
3299 :
3300 : !LDOS model
3301 0 : if (this%iprcel == 200 .or. this%iprcel == 212) then
3302 0 : call apply_chi0_ldos(this, dtset, vec_r)
3303 : end if
3304 :
3305 0 : end subroutine apply_chi0
3306 :
3307 : !!***
3308 : !!****f* ABINIT/apply_adjdielmat
3309 : !! NAME
3310 : !! apply_adjdielmat
3311 : !!
3312 : !! FUNCTION
3313 : !! Apply the adjoint dielectric matrix I-chi0*vc to the density rho_r (given in the Fourier space)
3314 : !! with a model chi0 operator.
3315 : !!
3316 : !! INPUTS
3317 : !! dtset = All input variables for this dataset.
3318 : !! mpi_enreg = Information about MPI parallelization.
3319 : !! rho_r = Density vector (in direct space, in Pauli basis).
3320 : !!
3321 : !! OUTPUT
3322 : !! adjdielmat_rho_r = adjdielmat * rho_r
3323 : !!
3324 : !! NOTES
3325 : !!
3326 : !! SOURCE
3327 90 : subroutine apply_adjdielmat(this, dtset, mpi_enreg, rho_r, adjdielmat_rho_r)
3328 :
3329 : !Arguments ------------------------------------
3330 : class(precon_object) :: this
3331 : !scalars
3332 : type(dataset_type),intent(in) :: dtset
3333 : type(MPI_type),intent(in) :: mpi_enreg
3334 : !arrays
3335 : real(dp), intent(in) :: rho_r(this%nfftprc, dtset%nspden)
3336 : real(dp), intent(inout) :: adjdielmat_rho_r(this%nfftprc, dtset%nspden)
3337 : !Local variables-------------------------------
3338 90 : real(dp), allocatable :: chi0_kxc_rho_r(:, :)
3339 :
3340 : ! *************************************************************************
3341 90 : if (this%use_precon) then
3342 :
3343 90 : if (this%iprcel == 299) then
3344 : ! P=I : No preconditioning
3345 0 : adjdielmat_rho_r = rho_r
3346 :
3347 90 : elseif (this%iprcel == 200) then
3348 : ! More efficient implementation for the LDOS preconditioner.
3349 35 : call apply_adjdielmat_ldos(this, dtset, mpi_enreg, rho_r, adjdielmat_rho_r)
3350 :
3351 55 : elseif (this%iprcel == 202) then
3352 : ! When iprcel = 202, P = (I - chi0_ldos*vc - chi0_diag*Kxc)
3353 :
3354 : !1) Compute adjdielmat_rho_r = rho_r - chi0_ldos * vc *rho_r
3355 55 : call apply_adjdielmat_ldos(this, dtset, mpi_enreg, rho_r, adjdielmat_rho_r)
3356 :
3357 : !2) Add -(chi0_diag * Kxc * rho_r) to adjdielmat_rho_r
3358 220 : ABI_MALLOC(chi0_kxc_rho_r, (this%nfftprc, dtset%nspden))
3359 :
3360 : !2.2) Apply Kxc to vec_r
3361 55 : call apply_kxc(this, dtset, mpi_enreg, rho_r, chi0_kxc_rho_r)
3362 : !2.3) Apply chi0_diag to Kxc*vec_r (in place)
3363 55 : call apply_chi0_diag(this, dtset, mpi_enreg, chi0_kxc_rho_r)
3364 : !2.5) Add this contribution to adjdielmat_rho_r
3365 371415 : adjdielmat_rho_r = adjdielmat_rho_r - chi0_kxc_rho_r
3366 :
3367 55 : ABI_FREE(chi0_kxc_rho_r)
3368 :
3369 0 : elseif (this%iprcel == 203) then
3370 : ! When iprcel = 203, P = (I - chi0_ldos*vc - chi0_quasidiag*Kxc)
3371 :
3372 : !1) Compute adjdielmat_rho_r = rho_r - chi0_ldos * vc *rho_r
3373 0 : call apply_adjdielmat_ldos(this, dtset, mpi_enreg, rho_r, adjdielmat_rho_r)
3374 :
3375 : !2) Add -(chi0_diag * Kxc * rho_r) to adjdielmat_rho_r
3376 0 : ABI_MALLOC(chi0_kxc_rho_r, (this%nfftprc, dtset%nspden))
3377 :
3378 : !2.2) Apply Kxc to vec_r
3379 0 : call apply_kxc(this, dtset, mpi_enreg, rho_r, chi0_kxc_rho_r)
3380 : !2.3) Apply chi0_quasidiag to Kxc*vec_r (in place)
3381 0 : call apply_chi0_quasidiag(this, dtset, mpi_enreg, chi0_kxc_rho_r)
3382 : !2.5) Add this contribution to adjdielmat_rho_r
3383 0 : adjdielmat_rho_r = adjdielmat_rho_r - chi0_kxc_rho_r
3384 :
3385 0 : ABI_FREE(chi0_kxc_rho_r)
3386 :
3387 : else
3388 : ! In the general case, P = (I - K*chi0_model) where K and chi0_model are defined
3389 : ! in the subroutine apply_kernel and apply_chi0 (depending on iprcel).
3390 :
3391 0 : adjdielmat_rho_r = rho_r
3392 : !1) Apply the Kernel (vc or vc + Kxc depending on iprcel)
3393 0 : call apply_kernel(this, dtset, mpi_enreg, adjdielmat_rho_r)
3394 : !2) Applythe model chi0 operator
3395 0 : call apply_chi0(this, dtset, adjdielmat_rho_r)
3396 : !3) adjdielmat_rho_r = rho_r - K * chi0 * rho_r = adjdielmat * rho_r
3397 0 : adjdielmat_rho_r = rho_r - adjdielmat_rho_r
3398 :
3399 : end if
3400 :
3401 : end if
3402 90 : end subroutine apply_adjdielmat
3403 : !!***
3404 :
3405 : !!****f* ABINIT/apply_dielmat
3406 : !! NAME
3407 : !! apply_dielmat
3408 : !!
3409 : !! FUNCTION
3410 : !! Apply the model dielectric matrix (I - K*chi0) to the potential v_r (given in the direct/real space),
3411 : !! with a model chi0 operator (contained in precon).
3412 : !!
3413 : !! INPUTS
3414 : !! dtset = All input variables for this dataset.
3415 : !! mpi_enreg = Information about MPI parallelization.
3416 : !! v_r = Potential vector (in real space, in Pauli basis)
3417 : !!
3418 : !! OUTPUT
3419 : !! dielmat_v_r = dielmat * v_r
3420 : !!
3421 : !! NOTES
3422 : !!
3423 : !! SOURCE
3424 0 : subroutine apply_dielmat(this, dtset, mpi_enreg, v_r, dielmat_v_r)
3425 :
3426 : !Arguments ------------------------------------
3427 : class(precon_object) :: this
3428 : !scalars
3429 : type(dataset_type),intent(in) :: dtset
3430 : type(MPI_type),intent(in) :: mpi_enreg
3431 : !arrays
3432 : real(dp), intent(in) :: v_r(this%nfftprc, dtset%nspden)
3433 : real(dp), intent(inout) :: dielmat_v_r(this%nfftprc, dtset%nspden)
3434 : !Local variables-------------------------------
3435 0 : real(dp), allocatable :: kxc_chi0_v_r(:, :), chi0_v_r(:, :)
3436 :
3437 : ! *************************************************************************
3438 0 : if (this%use_precon) then
3439 :
3440 0 : if (this%iprcel == 299) then
3441 : ! P=I : No preconditioning
3442 0 : dielmat_v_r = v_r
3443 :
3444 0 : elseif (this%iprcel == 200) then
3445 : ! More efficient implementation for the LDOS preconditioner.
3446 0 : call apply_dielmat_ldos(this, dtset, mpi_enreg, v_r, dielmat_v_r)
3447 :
3448 0 : elseif (this%iprcel == 202) then
3449 : ! When iprcel = 202, P = (I - vc*chi0_ldos - Kxc*chi0_diag)
3450 :
3451 : !1) Compute dielmat_v_r = v_r - vc * chi0_ldos *v_r
3452 0 : call apply_dielmat_ldos(this, dtset, mpi_enreg, v_r, dielmat_v_r)
3453 :
3454 : !2) Add -(Kxc * chi0_quasidiag * v_r) to dielmat_v_r
3455 :
3456 : !2.1) Apply chi0_quasidiag to v_r
3457 0 : ABI_MALLOC(chi0_v_r, (this%nfftprc, dtset%nspden))
3458 0 : chi0_v_r = v_r
3459 0 : call apply_chi0_diag(this, dtset, mpi_enreg, chi0_v_r)
3460 : !2.2) Apply Kxc to chi0*v_r
3461 0 : ABI_MALLOC(kxc_chi0_v_r, (this%nfftprc, dtset%nspden))
3462 0 : call apply_kxc(this, dtset, mpi_enreg, chi0_v_r, kxc_chi0_v_r)
3463 : !2.3) Add this contribution to dielmat_v_r
3464 0 : dielmat_v_r = dielmat_v_r - kxc_chi0_v_r
3465 0 : ABI_FREE(kxc_chi0_v_r)
3466 0 : ABI_FREE(chi0_v_r)
3467 :
3468 0 : elseif (this%iprcel == 203) then
3469 : ! When iprcel = 203 , P = (I - vc*chi0_ldos - Kxc*chi0_quasidiag)
3470 :
3471 : !1) Compute dielmat_v_r = v_r - vc * chi0_ldos *v_r
3472 0 : call apply_dielmat_ldos(this, dtset, mpi_enreg, v_r, dielmat_v_r)
3473 :
3474 : !2) Add -(Kxc * chi0_quasidiag * v_r) to dielmat_v_r
3475 :
3476 : !2.1) Apply chi0_quasidiag to v_r
3477 0 : ABI_MALLOC(chi0_v_r, (this%nfftprc, dtset%nspden))
3478 0 : chi0_v_r = v_r
3479 0 : call apply_chi0_quasidiag(this, dtset, mpi_enreg, chi0_v_r)
3480 : !2.2) Apply Kxc to chi0*v_r
3481 0 : ABI_MALLOC(kxc_chi0_v_r, (this%nfftprc, dtset%nspden))
3482 0 : call apply_kxc(this, dtset, mpi_enreg, chi0_v_r, kxc_chi0_v_r)
3483 : !2.3) Add this contribution to dielmat_v_r
3484 0 : dielmat_v_r = dielmat_v_r - kxc_chi0_v_r
3485 0 : ABI_FREE(kxc_chi0_v_r)
3486 0 : ABI_FREE(chi0_v_r)
3487 :
3488 : else
3489 : ! In the general case, P = (I - K*chi0_model) where K and chi0_model are defined
3490 : ! in the subroutine apply_kernel and apply_chi0 (depending on iprcel).
3491 :
3492 0 : dielmat_v_r = v_r
3493 : !1) Apply the model chi0 operator
3494 0 : call apply_chi0(this, dtset, dielmat_v_r)
3495 : !2) Apply the Kernel (vc or vc + Kxc depending on iprcel)
3496 0 : call apply_kernel(this, dtset, mpi_enreg, dielmat_v_r)
3497 : !3) dielmat_v_r = v_r - K * chi0 * v_r = dielmat * v_r
3498 0 : dielmat_v_r = v_r - dielmat_v_r
3499 :
3500 : end if
3501 :
3502 : end if
3503 0 : end subroutine apply_dielmat
3504 : !!***
3505 :
3506 : !****f* m_precon/apply_precon
3507 : !! NAME
3508 : !! apply_precon
3509 : !!
3510 : !! FUNCTION
3511 : !! Apply a the preconditioner P^-1 to a given input vector 'vresid' where P is a model for the
3512 : !! dielectric matrix or its adjoint based of a model of the non interacting susceptibility chi0.
3513 : !! More precisely,
3514 : !! - if we are preconditioning potentials ('optres'=0), P is a model of the dielectric matrix (I-K*chi0)
3515 : !! - if we are preconditioning densities ('optres'=1), P is a model of the adjoint dielectric matrix (I-chi0*K).
3516 : !! The approximation are defined in 'apply_dielmat' and 'apply_adjdielmat' by Abinit input 'iprcel'.
3517 : !!
3518 : !! The preconditioner is applied by solving the linear equation P * 'vrespc' = 'vresid' iteratively,
3519 : !! either using the GMRES method if P is well conditioned (positive definite) or using Ridge/Tikhonov regularization
3520 : !! with the conjugate gradient if P might be ill-conditioned.
3521 : !!
3522 : !! INPUTS
3523 : !! dtset = All input variables for this dataset.
3524 : !! mpi_enreg = Information about MPI parallelization.
3525 : !! optreal = Integer flag indicating whether the input is in real space (1) or reciprocal space (2).
3526 : !! optres = Integer flag indicating whether we are preconditioning densities (1) or potentials (0).
3527 : !! vresid = Residual vector to which the preconditioner is applied.
3528 : !!
3529 : !! OUTPUTS
3530 : !! vrespc = Preconditioned residual vector.
3531 : !!
3532 : !! NOTES
3533 : !! 'vresid' and 'vrespc' have a different shape than the typical density/potential vectors
3534 : !! in the rest of this file, to match the shape needed in 'm_prcref'.
3535 : !!
3536 : !! SOURCE
3537 18 : subroutine apply_precon(this, dtset, mpi_enreg, optreal, optres, vresid, vrespc)
3538 : !Arguments ------------------------------------
3539 : class(precon_object), intent(inout) :: this
3540 : type(dataset_type),intent(in) :: dtset
3541 : type(MPI_type),intent(in) :: mpi_enreg
3542 : integer :: optreal, optres
3543 : !arrays
3544 : real(dp), intent(in) :: vresid(optreal*this%nfftprc, dtset%nspden)
3545 : real(dp), intent(inout) :: vrespc(optreal*this%nfftprc, dtset%nspden)
3546 : !Local variables-------------------------------
3547 : !scalars
3548 : integer :: ispden, start_ispden, end_ispden, n
3549 : !arrays
3550 18 : real(dp), allocatable :: rhs(:), est(:), P_rhs(:)
3551 18 : real(dp), allocatable :: work_g(:, :, :)
3552 :
3553 : ! *************************************************************************
3554 36 : if (this%use_precon) then
3555 :
3556 : !0.0) Update precon object
3557 18 : call precon_update(this, dtset, mpi_enreg)
3558 :
3559 : ! The preconditioned density/potential residual vrespc = P^-1 * vresid is computed
3560 : ! by sovling the linear equation P * vrespc = vresid approximately with GMRES.
3561 :
3562 18 : n = dtset%nspden*1*this%nfftprc
3563 :
3564 : !0.1) Convert the input to direct/real space if needed.
3565 18 : if (optreal==0) then
3566 : ! vresid is given in the Fourier space : We need to do a ifft.
3567 0 : ABI_MALLOC(work_g, (2, this%nfftprc, dtset%nspden))
3568 0 : work_g = reshape(vresid, (/2,this%nfftprc, dtset%nspden/))
3569 0 : call fourdp(1, work_g, vrespc(1:this%nfftprc, :), 1, mpi_enreg, this%nfftprc, dtset%nspden, this%ngfftprc, 0)
3570 : else
3571 121554 : vrespc = vresid
3572 : end if
3573 :
3574 : !0.2) Convert the input to the Pauli basis
3575 18 : call to_pauli(optres, vrespc)
3576 :
3577 : !1) Right-hand side : rhs is vresid (flattened) in the direct/real space.
3578 54 : ABI_MALLOC(rhs, (n))
3579 54 : do ispden = 1, dtset%nspden
3580 : ! Indices of the ispden component in the flattened (this%nfftprc, dtset%nspden)-array 'rhs'.
3581 36 : start_ispden = 1+(ispden-1)*this%nfftprc
3582 36 : end_ispden = ispden*this%nfftprc
3583 121554 : rhs(start_ispden:end_ispden) = vrespc(1:this%nfftprc, ispden)
3584 : end do
3585 :
3586 : !2) Initial guess :
3587 36 : ABI_MALLOC(est, (n))
3588 121518 : est = 0
3589 : ! Is est = rhs a better starting point ?
3590 : !est = rhs
3591 :
3592 : !3) Resolution of the linear system :
3593 18 : if (this%precon_verbose>0) call wrtout(std_out, '-- chi0-based preconditioning - Linear Solver --')
3594 18 : if (this%use_ridgereg) then
3595 : ! P is ill-conditionned :
3596 : ! Ridge/Tikhonov regularization and CG :
3597 : ! We solve (P^*P + ridge_param*I) * est = P * rhs
3598 : ! (P^*P + ridge_param*I) is self-adjoint and can be solved with CG.
3599 0 : ABI_MALLOC(P_rhs, (n))
3600 0 : call ridge_matvec(n, rhs, P_rhs)
3601 0 : call cg_linear_solver(n, ridge_matvec, P_rhs, est, (this%linsolve_maxiter-1)/2+1, this%linsolve_rtol, this%precon_verbose>0)
3602 0 : ABI_FREE(P_rhs)
3603 : else
3604 : ! P is well conditionned :
3605 : ! GMRES (P is not self-adjoint)
3606 18 : call gmres_linear_solver(n, matvec, rhs, est, this%linsolve_maxiter, this%linsolve_rtol, this%precon_verbose>0)
3607 :
3608 : end if
3609 :
3610 : !4.0) Free arrays that might have been allocated by precon_update
3611 18 : call precon_free_update(this)
3612 :
3613 : !4) Reshaping the final result :
3614 54 : do ispden = 1, dtset%nspden
3615 : ! Indices of the ispden component in the flattened (this%nfftprc, dtset%nspden)-array 'est'.
3616 36 : start_ispden = 1+(ispden-1)*this%nfftprc
3617 36 : end_ispden = ispden*this%nfftprc
3618 54 : if (optreal==1) then
3619 : ! vrespc must be returned in the direct/real space.
3620 121536 : vrespc(:, ispden) = est(start_ispden:end_ispden)
3621 : else
3622 : ! vrespc must be returned in the fourier space, we need to do a fft.
3623 0 : call fourdp(1, work_g, est(start_ispden:end_ispden), -1, mpi_enreg, this%nfftprc, 1, this%ngfftprc, 0)
3624 0 : vrespc(:, ispden) = reshape(work_g, (/2*this%nfftprc/))
3625 : end if
3626 : end do
3627 : ! vrespc must be returned in the default Abinit spin-basis.
3628 18 : call from_pauli(optres, vrespc)
3629 :
3630 18 : ABI_FREE(rhs)
3631 18 : ABI_FREE(est)
3632 :
3633 : end if
3634 :
3635 : contains
3636 :
3637 : ! Subroutine matvec that applies the model (adjoint-) dielectric matrix. -----------
3638 90 : subroutine matvec(n_, x, y)
3639 : integer, intent(in) :: n_
3640 : real(dp), intent(inout), target :: x(n_), y(n_)
3641 : type(c_ptr) :: x_c, y_c
3642 90 : real(dp), pointer :: x_2d(:, :), y_2d(:, :)
3643 :
3644 : ! **********************************************************************************
3645 :
3646 : ! C-pointers to match the flattened arrays x and y to their 3D versions needed by
3647 : ! 'apply_adjdielmat' and 'apply_dielmat'.
3648 90 : x_c = c_loc(x)
3649 270 : call c_f_pointer(x_c, x_2d, shape=[this%nfftprc, dtset%nspden])
3650 90 : y_c = c_loc(y)
3651 270 : call c_f_pointer(y_c, y_2d, shape=[this%nfftprc, dtset%nspden])
3652 :
3653 90 : if (optres==1) then
3654 : ! We are preconditioning density residual so P models the adjoint dielectric matrix.
3655 90 : call this%apply_adjdielmat(dtset, mpi_enreg, x_2d, y_2d)
3656 0 : else if (optres==0) then
3657 : ! We are preconditioning potential residual so P models the dielectric matrix.
3658 0 : call this%apply_dielmat(dtset, mpi_enreg, x_2d, y_2d)
3659 : end if
3660 :
3661 90 : end subroutine matvec ! ------------------------------------------------------------
3662 :
3663 : ! Subroutine ridge_matvec that applies the operator (P^*P + ridge_param*I) needed for ridge regularization.
3664 0 : subroutine ridge_matvec(n_, x, y)
3665 : integer, intent(in) :: n_
3666 : real(dp), intent(inout), target :: x(n_), y(n_)
3667 : type(c_ptr) :: x_c, y_c
3668 0 : real(dp), pointer :: x_2d(:, :), y_2d(:, :), temp_2d(:, :)
3669 :
3670 : ! **********************************************************************************
3671 :
3672 : ! C pointers to match the flattened arrays x and y to their 3D versions needed by
3673 : ! 'apply_adjdielmat' and 'apply_dielmat'.
3674 0 : x_c = c_loc(x)
3675 0 : call c_f_pointer(x_c, x_2d, shape=[this%nfftprc, dtset%nspden])
3676 0 : y_c = c_loc(y)
3677 0 : call c_f_pointer(y_c, y_2d, shape=[this%nfftprc, dtset%nspden])
3678 0 : ABI_MALLOC(temp_2d, (this%nfftprc, dtset%nspden)) ! Temporary array for intermediate result
3679 :
3680 0 : if (optres==1) then
3681 : ! We are preconditioning density residual so P models the adjoint dielectric matrix.
3682 0 : call this%apply_adjdielmat(dtset, mpi_enreg, x_2d, temp_2d)
3683 0 : call this%apply_dielmat(dtset, mpi_enreg, temp_2d, y_2d)
3684 0 : else if (optres==0) then
3685 : ! We are preconditioning potential residual so P models the dielectric matrix.
3686 0 : call this%apply_dielmat(dtset, mpi_enreg, x_2d, temp_2d)
3687 0 : call this%apply_adjdielmat(dtset, mpi_enreg, temp_2d, y_2d)
3688 : end if
3689 0 : y_2d = y_2d + this%ridge_param*x_2d
3690 0 : ABI_FREE(temp_2d)
3691 :
3692 0 : end subroutine ridge_matvec ! ---------------------------------------------------------
3693 :
3694 : end subroutine apply_precon
3695 :
3696 324414 : end module m_precon
|
