Branch data 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 [ - + - + : 20859 : subroutine precon_init(this, dtset, atindx, atindx1, cg, cprj, dimcprj, dtfil, eigen, fermie, &
- + - + ]
181 [ - + - + : 20859 : & gmet, gprimd, indsym, irrzon, kg, mcprj, nattyp, nfftmix, ngfftmix, npwarr, occ, pawang, pawfgr, pawfgrtab, &
- + - + -
+ - + - +
- + - + ]
182 [ - + - + : 20859 : & 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 : 6953 : this%iprcel = dtset%iprcel
215 : 6953 : this%use_precon = .false.
216 : 6953 : this%use_kxc = .false.
217 : :
218 [ + + ]: 6953 : 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 : 6953 : 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 [ - + ]: 6953 : 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 [ + + ]: 6953 : 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 : 6953 : 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 : 6953 : subroutine precon_free(this)
501 : :
502 : : !Arguments ------------------------------------
503 : : class(precon_object), intent(inout) :: this
504 : :
505 : : ! *************************************************************************
506 [ + + ]: 6953 : 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 : 6953 : 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 [ - + - - : 324519 : end module m_precon
+ + - + +
- - + - +
+ - - + -
+ - + - +
- + - + ]
|
