Line data Source code
1 : !!****m* ABINIT/m_rmm_diis
2 : !! NAME
3 : !! m_rmm_diis
4 : !!
5 : !! FUNCTION
6 : !! This module contains routines for the RMM-DIIS eigenvalue solver.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2020-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_rmm_diis
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_xmpi
27 : use m_abicore
28 : use m_dtset
29 : use m_cgtools
30 : use m_hide_blas
31 : use m_yaml
32 : use m_linalg_interfaces
33 : use m_prep_kgb
34 :
35 : use defs_abitypes, only : mpi_type
36 : use m_fstrings, only : sjoin, itoa, ftoa
37 : use m_time, only : timab, cwtime, cwtime_report
38 : use m_numeric_tools, only : pack_matrix, imin_loc, stats_t, stats_eval
39 : use m_hide_lapack, only : xhegv_cplex, xhesv_cplex
40 : use m_pair_list, only : pair_list
41 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free
42 : use m_fftcore, only : fftcore_set_mixprec
43 : use m_hamiltonian, only : gs_hamiltonian_type
44 : use m_getghc, only : getghc
45 : use m_nonlop, only : nonlop
46 : use m_cgtk, only : cgtk_fixphase
47 : use m_abi_linalg, only : abi_zgemm_2r
48 : !use m_fock, only : fock_set_ieigen, fock_set_getghc_call
49 :
50 : implicit none
51 :
52 : private
53 : !!***
54 :
55 : public :: rmm_diis
56 : public :: subspace_rotation ! rayleigh-ritz procedure from gs_hamk
57 : !!***
58 :
59 : type,private :: rmm_diis_t
60 :
61 : integer :: accuracy_level
62 : ! Defines tolerances, activates/deactivates tricks.
63 :
64 : integer :: usepaw
65 : ! 1 if we are running PAW.
66 :
67 : integer :: istwf_k
68 : ! wavefunction storage mode.
69 :
70 : integer :: cplex
71 : ! 1 if matrices are real (e.g. Gamma-point), 2 for complex
72 :
73 : integer :: bsize
74 : ! (Max) block size for bands
75 :
76 : integer :: max_niter
77 : ! Maximum number of iterations
78 :
79 : integer :: npwsp
80 : ! Total number of planewaves treated by this proc
81 : ! npw * my_nspinor
82 :
83 : integer :: prtvol
84 : ! vervosity level
85 :
86 : integer :: last_iter
87 : ! Last RMM-DIIS iteration performed.
88 :
89 : real(dp) :: tol_occupied
90 : ! Tolerance for partial occupied states
91 :
92 : type(pair_list) :: stats
93 :
94 : real(dp),allocatable :: hist_ene(:,:)
95 : real(dp),allocatable :: hist_resid(:,:)
96 : real(dp),allocatable :: hist_enlx(:,:)
97 : character(len=7),allocatable :: step_type(:,:)
98 : ! (0:max_niter+2, bsize)
99 : ! 0 is the initial step, then DIIS iterations whose number may depend on the block
100 : ! followed by the computation of eigens after ortho.
101 :
102 : real(dp),allocatable :: resmat(:,:,:,:)
103 : ! (2, 0:max_niter, 0:max_niter, bsize))
104 :
105 : real(dp),allocatable :: chain_phi(:,:,:,:)
106 : real(dp),allocatable :: chain_sphi(:,:,:,:)
107 : real(dp),allocatable :: chain_resv(:,:,:,:)
108 : ! (2, npwsp, 0:max_niter, bsize))
109 :
110 : contains
111 : procedure :: free => rmm_diis_free ! Free dynamic memory
112 : procedure :: update_block => rmm_diis_update_block ! DIIS uppdate of wavefuntions and residuals.
113 : procedure :: eval_mats => rmm_diis_eval_mats ! Compute DIIS matrices
114 : procedure :: exit_iter => rmm_diis_exit_iter ! Return True if can exit the DIIS iteration.
115 : procedure :: print_block => rmm_diis_print_block ! Print energies, residuals and diffs for a given block.
116 : ! TODO: Fix problem with last_iter and hist
117 : procedure :: push_iter => rmm_diis_push_iter ! Save results required by the DIIS algorithm
118 :
119 : end type rmm_diis_t
120 :
121 : integer,parameter, private :: level = 432
122 : logical,parameter, private :: timeit = .False.
123 : !logical,parameter, private :: timeit = .True.
124 :
125 : contains
126 : !!***
127 :
128 : !!****f* ABINIT/rmm_diis
129 : !! NAME
130 : !! rmm_diis
131 : !!
132 : !! FUNCTION
133 : !! This routine updates the wave functions at a given (k-point, spin), using the RMM-DIIS method.
134 : !!
135 : !! INPUTS
136 : !! istep,ikpt,isppol=Iteration step, k-point index, spin index (mainly for printing purposes).
137 : !! dtset <type(dataset_type)>=all input variables for this dataset
138 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the hamiltonian at k
139 : !! kinpw(npw)=(modified) kinetic energy for each plane wave (hartree)
140 : !! mpi_enreg=information about MPI parallelization
141 : !! nband=number of bands at this k point for that spin polarization
142 : !! npw=number of plane waves at this k point
143 : !! my_nspinor=number of spinors treated by this MPI proc
144 : !!
145 : !! OUTPUT
146 : !! eig(nband)=array for holding eigenvalues (hartree)
147 : !! If usepaw==1:
148 : !! gsc(2,*)=<g|s|c> matrix elements (s=overlap)
149 : !! If usepaw==0
150 : !! enlx(nband)=contribution from each band to nonlocal psp + potential Fock ACE part
151 : !! of total energy, at this k-point
152 : !!
153 : !! SIDE EFFECTS
154 : !! cg(2,*)=updated wavefunctions
155 : !! resid(nband)=residuals for each states. In input: previous residuals for this k-point, spin.
156 : !! In output: new residuals.
157 : !! rmm_diis_status(2): Status of the eigensolver.
158 : !! The first entry gives the previous accuracy.
159 : !! The second entry gives the number of iterations already performed with this level.
160 : !!
161 : !! SOURCE
162 :
163 653 : subroutine rmm_diis(istep, ikpt, isppol, cg, dtset, eig, occ, enlx, gs_hamk, kinpw, gsc, &
164 653 : mpi_enreg, nband, npw, my_nspinor, resid, rmm_diis_status)
165 :
166 : !Arguments ------------------------------------
167 : integer,intent(in) :: istep, ikpt, isppol, nband, npw, my_nspinor
168 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
169 : type(dataset_type),intent(in) :: dtset
170 : type(mpi_type),intent(in) :: mpi_enreg
171 : real(dp),target,intent(inout) :: cg(2,npw*my_nspinor*nband)
172 : real(dp),target,intent(inout) :: gsc(2,npw*my_nspinor*nband*dtset%usepaw)
173 : real(dp),intent(inout) :: enlx(nband), resid(nband)
174 : real(dp),intent(in) :: occ(nband), kinpw(npw)
175 : real(dp),intent(out) :: eig(nband)
176 : integer,intent(inout) :: rmm_diis_status(2)
177 :
178 : !Local variables-------------------------------
179 : integer,parameter :: type_calc0 = 0, option1 = 1, option2 = 2, tim_getghc = 0
180 : integer,parameter :: choice1 = 1, signs1 = 1, signs2 = 2, tim_nonlop = 0, paw_opt0 = 0, paw_opt3 = 3
181 : integer :: ierr, prtvol, bsize, nblocks, iblock, npwsp, ndat, ib_start, ib_stop, idat, paral_kgb !, ortalgo
182 : integer :: cpopt, sij_opt, igs, ige, mcg, mgsc, istwf_k, optekin, usepaw, iter, max_niter, max_niter_block
183 : integer :: me_g0, nb_pocc, jj, kk, accuracy_level, raise_acc, prev_mixprec, after_ortho, me_cell
184 : integer :: comm_bsf, prev_accuracy_level, ncalls_with_prev_accuracy, signs, paw_opt, savemem
185 : logical :: first_call, use_fft_mixprec, has_fock
186 : real(dp),parameter :: rdummy = zero
187 : real(dp) :: accuracy_ene, max_res_pocc, tol_occupied, lock_tolwfr !, max_absimag
188 : real(dp) :: cpu, wall, gflops, cpu_all, wall_all, gflops_all
189 : character(len=500) :: msg
190 : type(yamldoc_t) :: rmm_ydoc
191 653 : type(rmm_diis_t) :: diis
192 : type(stats_t) :: res_stats
193 : !arrays
194 : real(dp) :: tsec(2)
195 : real(dp),target :: fake_gsc_bk(0,0)
196 653 : real(dp),allocatable :: lambda_bk(:), kres_bk(:,:), dots_bk(:,:), residv_bk(:,:)
197 653 : real(dp),allocatable :: umat(:,:,:), gwork(:,:), dots(:, :)
198 653 : real(dp),target,allocatable :: ghc(:,:), gvnlxc(:,:)
199 1306 : real(dp),contiguous, pointer :: gsc_bk(:,:), cg_bk(:,:), ghc_bk(:,:), gvnlxc_bk(:,:)
200 4571 : type(pawcprj_type) :: cprj_dum(1,1)
201 : ! *************************************************************************
202 :
203 : ! Define useful vars.
204 653 : usepaw = dtset%usepaw; istwf_k = gs_hamk%istwf_k; paral_kgb = mpi_enreg%paral_kgb
205 653 : me_g0 = mpi_enreg%me_g0; comm_bsf = mpi_enreg%comm_bandspinorfft
206 653 : npwsp = npw * my_nspinor; mcg = npwsp * nband; mgsc = npwsp * nband * usepaw
207 653 : me_cell = mpi_enreg%me_cell; prtvol = dtset%prtvol !; prtvol = -level
208 653 : has_fock = associated(gs_hamk%fockcommon)
209 :
210 : if (timeit) then
211 : call cwtime(cpu_all, wall_all, gflops_all, "start")
212 : call cwtime(cpu, wall, gflops, "start")
213 : end if
214 :
215 : ! =================
216 : ! Prepare DIIS loop
217 : ! =================
218 : ! accuracy_level is computed from the maxval of the previous residuals received in input.
219 : ! The different levels are:
220 : !
221 : ! 1: Used at the beginning of the SCF cycle. Use loosy convergence criteria in order
222 : ! to reduce the number of H|psi> applications as much as possible so that
223 : ! we can start to mix densities/potentials.
224 : ! Allow for incosistent data in output rediduals and Vnl matrix elements.
225 : ! Move to the next level after Max 15 iterations.
226 : !
227 : ! 2: Intermediate step. Decrease convergence criteria in order to perform more wavefuction iterations.
228 : ! Allow for incosistent data in output rediduals and Vnl matrix elements.
229 : ! Move to the next level after Max 25 iterations.
230 : !
231 : ! 3: Approaching convergence. Use stricter convergence criteria.
232 : ! Move to the next level after Max 25 iterations.
233 : !
234 : ! 4: Ultimate precision. Try to reach the same accuracy as the other eigenvalue solvers.
235 : ! This implies: using similar convergence criteria as in the other solvers.
236 :
237 : ! Note:
238 : !
239 : ! * Accuracy_level is not allowed to increase during the SCF cycle.
240 : !
241 : ! * Since we operate on blocks of bands, all the states in the block will receive the same treatment.
242 : ! This means that one can observe different convergence behaviour depending on bsize.
243 : !
244 917 : if (all(rmm_diis_status == 0)) then
245 : ! This is the first time we call rmm_diis for this (k-point, spin)
246 : prev_accuracy_level = 1; ncalls_with_prev_accuracy = 0
247 : first_call = .True.
248 : else
249 521 : prev_accuracy_level = rmm_diis_status(1); ncalls_with_prev_accuracy = rmm_diis_status(2)
250 521 : first_call = .False.
251 : end if
252 :
253 : ! Decide whether we should move to the next level.
254 653 : raise_acc = 0
255 653 : if (prev_accuracy_level == 1 .and. ncalls_with_prev_accuracy >= 15) raise_acc = 2
256 653 : if (prev_accuracy_level == 2 .and. ncalls_with_prev_accuracy >= 25) raise_acc = 3
257 653 : if (prev_accuracy_level == 3 .and. ncalls_with_prev_accuracy >= 25) raise_acc = 4
258 653 : if (raise_acc > 0) then
259 0 : ABI_COMMENT("Accuracy_level is automatically increased as we reached the max number of NSCF iterations.")
260 : end if
261 653 : raise_acc = max(raise_acc, prev_accuracy_level)
262 :
263 : ! Define tolerance for occupied states on the basis of prev_accuracy_level
264 : ! and compute max of residuals for these bands.
265 653 : tol_occupied = zero
266 653 : if (dtset%iscf > 0) then
267 588 : tol_occupied = tol3; if (any(prev_accuracy_level == [1])) tol_occupied = tol2
268 : end if
269 13141 : nb_pocc = count(occ >= tol_occupied)
270 10498 : max_res_pocc = maxval(resid(1:nb_pocc))
271 :
272 : ! Define accuracy_level for this run.
273 653 : accuracy_level = 1
274 653 : if (max_res_pocc < tol8) accuracy_level = 2
275 522 : if (max_res_pocc < tol12) accuracy_level = 3
276 195 : if (max_res_pocc < tol16) accuracy_level = 4
277 : !if (max_res_pocc < tol18) accuracy_level = 4
278 653 : accuracy_level = max(prev_accuracy_level, accuracy_level, raise_acc)
279 653 : if (istep == 1) accuracy_level = 2 ! FIXME: Differenciate between restart or rmm_diis - 3.
280 653 : if (first_call .and. max_res_pocc == zero) accuracy_level = 1
281 : !print *, "rmm_diis_status:", rmm_diis_status
282 : !print *, "rmm_prev_acc:", prev_accuracy_level, "rmm_raise_acc:", raise_acc
283 : !print *, "accuracy_level:", accuracy_level, "rmm_raise_acc:", raise_acc
284 :
285 : ! Update rmm_diis_status. Reset number of calls if we've just moved to a new accuracy_level.
286 653 : rmm_diis_status(1) = accuracy_level
287 653 : if (accuracy_level /= prev_accuracy_level) rmm_diis_status(2) = 0
288 653 : rmm_diis_status(2) = rmm_diis_status(2) + 1
289 :
290 : ! Will perform max_niter DIIS steps. Usually 3 as nline by default is 4.
291 : ! Note that, unlike in Vasp's recipe, here we don't end with a trial step after DIIS.
292 653 : max_niter = max(dtset%nline - 1, 1)
293 653 : if (accuracy_level >= 4) max_niter = dtset%nline
294 653 : if (dtset%iscf < 0) max_niter = dtset%nline + 1
295 :
296 : ! Define accuracy_ene for SCF.
297 653 : accuracy_ene = zero
298 653 : if (dtset%iscf > 0) then
299 588 : if (dtset%toldfe /= zero) then
300 200 : accuracy_ene = dtset%toldfe * ten**(-accuracy_level + 2) / nb_pocc
301 : else
302 : ! We are not using toldfe to stop the SCF cycle
303 : ! so we are forced to hardcode a tolerance for the absolute diff in the KS eigenvalue.
304 388 : accuracy_ene = tol8 * ten**(-accuracy_level + 2) / nb_pocc
305 : end if
306 : end if
307 :
308 : ! Tolerance on residuals used for band locking after subdiago.
309 653 : if (dtset%tolwfr > zero) then
310 117 : lock_tolwfr = tol2 * dtset%tolwfr
311 : else
312 536 : lock_tolwfr = tol14
313 536 : if (accuracy_level >= 2) lock_tolwfr = tol16
314 536 : if (accuracy_level >= 3) lock_tolwfr = tol18
315 536 : if (accuracy_level >= 4) lock_tolwfr = tol20 * tol2
316 : end if
317 :
318 : ! Use mixed precisions if requested by the user but only for low accuracy_level
319 653 : use_fft_mixprec = dtset%mixprec == 1 .and. accuracy_level < 2
320 0 : if (use_fft_mixprec) prev_mixprec = fftcore_set_mixprec(1)
321 :
322 : ! Select preconditioning.
323 : optekin = 0; if (dtset%wfoptalg >= 10) optekin = 1
324 653 : optekin = 1 ! optekin = 0
325 :
326 : ! Will treat states in groups of bsize bands even when paral_kgb = 0
327 653 : bsize = 8; if (paral_kgb == 1) bsize = mpi_enreg%nproc_band * mpi_enreg%bandpp
328 653 : nblocks = nband / bsize; if (mod(nband, bsize) /= 0) nblocks = nblocks + 1
329 :
330 : ! Build DIIS object.
331 653 : diis = rmm_diis_new(accuracy_level, usepaw, istwf_k, npwsp, max_niter, bsize, prtvol)
332 653 : diis%tol_occupied = tol_occupied
333 : !call wrtout(std_out, sjoin(" Using Max", itoa(max_niter), "RMM-DIIS iterations"))
334 : !call wrtout(std_out, sjoin( &
335 : ! " Max_input_resid_pocc", ftoa(max_res_pocc), "accuracy_level:", itoa(accuracy_level), &
336 : ! ", accuracy_ene: ", ftoa(accuracy_ene)))
337 653 : call timab(1634, 1, tsec) ! "rmm_diis:band_opt"
338 :
339 653 : rmm_ydoc = yamldoc_open("RMM-DIIS", with_iter_state=.False.)
340 3265 : call rmm_ydoc%add_ints("ikpt, isppol, istep, accuracy_level", [ikpt, isppol, istep, accuracy_level])
341 653 : call rmm_ydoc%open_tabular("RESIDS_POCC") !, tag, indent, newline, comment)
342 653 : write(msg, "(1x, a12, 4(a10))")"level", "mean", "min", "max", "stdev"
343 653 : call rmm_ydoc%add_tabular_line(msg, indent=0)
344 653 : call rmm_ydoc%add_tabular_line(resids2str("input"), indent=0)
345 :
346 : ! =========================
347 : ! === Subspace rotation ===
348 : ! =========================
349 : ! Allocate big (scalable) array with <G|H|C> for all nband so that we can recompute the residuals after the rotation.
350 : ! This approach requires more memory but we avoid one extra call to H|Psi> per band.
351 : ! Alternatively, one can compute ghc and the residuals by applying H|psi>
352 : ! inside the loop over blocks (less memory but slower).
353 653 : savemem = dtset%rmm_diis_savemem
354 : !savemem = 1
355 : !if (savemem == 0) then
356 : ! ABI_MALLOC_OR_DIE(ghc, (2, npwsp*nband), ierr)
357 : ! ABI_MALLOC_OR_DIE(gvnlxc, (2, npwsp*nband), ierr)
358 : !end if
359 :
360 : call subspace_rotation(gs_hamk, dtset%prtvol, mpi_enreg, nband, npw, my_nspinor, savemem, &
361 653 : enlx, eig, cg, gsc, ghc, gvnlxc)
362 :
363 653 : gsc_bk => fake_gsc_bk
364 653 : cpopt = -1; sij_opt = 0
365 653 : if (usepaw == 1) then
366 160 : sij_opt = 1 ! matrix elements <G|S|C> have to be computed in gsc in addition to ghc
367 : cpopt = -1 ! <p_lmn|in> (and derivatives) are computed here (and not saved)
368 : end if
369 :
370 1959 : ABI_MALLOC(lambda_bk, (bsize))
371 1959 : ABI_MALLOC(dots_bk, (2, bsize))
372 1959 : ABI_MALLOC(residv_bk, (2, npwsp*bsize))
373 1306 : ABI_MALLOC(kres_bk, (2, npwsp*bsize))
374 :
375 653 : if (savemem == 1) then
376 96 : ABI_MALLOC(ghc_bk, (2, npwsp*bsize))
377 96 : ABI_MALLOC(gvnlxc_bk, (2, npwsp*bsize))
378 : end if
379 : !write(msg, "(a,f8.1,a)")" Memory required: ", 2 * natom3**2 * (my_q2 - my_q1 + 1) * dp * b2Mb, " [Mb] <<< MEM"
380 : !call wrtout(std_out, msg)
381 :
382 : ! We loop over nblocks, each block contains ndat states.
383 : !
384 : ! - Convergence behaviour may depend on bsize as branches are taken according to
385 : ! the status of all bands in the block.
386 : ! TODO: Transpose only once per block and then work with already_transposed = .True.
387 : if (timeit) call cwtime(cpu, wall, gflops, "start")
388 :
389 1546 : do iblock=1,nblocks
390 893 : igs = 1 + (iblock - 1) * npwsp * bsize; ige = min(iblock * npwsp * bsize, npwsp * nband)
391 893 : ndat = (ige - igs + 1) / npwsp
392 893 : ib_start = 1 + (iblock - 1) * bsize; ib_stop = min(iblock * bsize, nband)
393 :
394 : ! Reduce number of niter iterations if block contains "empty" states.
395 : ! This should happen only if npband is small wrt nband and nband >> nbocc.
396 : ! TODO: Don't reduce niter if MD
397 893 : max_niter_block = max_niter
398 893 : if (dtset%iscf > 0) then
399 1884 : if (all(occ(ib_start:ib_stop) < diis%tol_occupied)) max_niter_block = max(1 + max_niter / 2, 2)
400 : end if
401 :
402 : ! Compute H |phi_0> with cg block after subdiago.
403 893 : cg_bk => cg(:,igs:ige); if (usepaw == 1) gsc_bk => gsc(1:2,igs:ige)
404 :
405 : ! Compute residual vectors after subspace_rotation.
406 893 : if (savemem == 0) then
407 845 : ghc_bk => ghc(:,igs:ige); gvnlxc_bk => gvnlxc(:,igs:ige)
408 845 : call cg_get_residvecs(usepaw, npwsp, ndat, eig(ib_start), cg_bk, ghc_bk, gsc_bk, residv_bk)
409 845 : call cg_norm2g(istwf_k, npwsp, ndat, residv_bk, resid(ib_start), me_g0, comm_bsf)
410 : else
411 : call getghc_eigresid(gs_hamk, npw, my_nspinor, ndat, cg_bk, ghc_bk, gsc_bk, mpi_enreg, prtvol, &
412 48 : eig(ib_start), resid(ib_start), enlx(ib_start), residv_bk, gvnlxc_bk, normalize=.False.)
413 : end if
414 :
415 : ! Band locking.
416 1073 : if (all(resid(ib_start:ib_stop) < lock_tolwfr)) then
417 0 : call diis%stats%increment("locked", ndat)
418 0 : cycle ! iblock
419 : end if
420 :
421 : ! Save <R0|R0> and <phi_0|S|phi_0>, |phi_0>, |S phi_0>. Assume input cg_bk is already S-normalized.
422 893 : call diis%push_iter(0, ndat, eig(ib_start), resid(ib_start), enlx(ib_start), cg_bk, residv_bk, gsc_bk, "SDIAG")
423 :
424 : ! Line minimization with preconditioned steepest descent:
425 : !
426 : ! |phi_1> = |phi_0> + lambda |K R_0>
427 : !
428 : ! where lambda minimizes the residual (we don't try to find the stationary
429 : ! point of the Rayleigh quotient as in Kresse's paper).
430 : !
431 : ! lambda = - Re{<R_0|(H - e_0 S)} |K R_0>} / |(H - e_0 S) |K R_0>|**2
432 : !
433 : ! more expensive than finding the stationary point of the Rayleigh quotient as it requires
434 : ! an extra H application but it should be more stable and more consistent with the RMM approach.
435 : !
436 : ! Precondition |R_0>, output in kres_bk = |K R_0>
437 893 : call cg_zcopy(npwsp * ndat, residv_bk, kres_bk)
438 893 : call cg_precon_many(istwf_k, npw, my_nspinor, ndat, cg_bk, optekin, kinpw, kres_bk, me_g0, comm_bsf)
439 :
440 : ! Compute H |K R_0>
441 893 : if (paral_kgb == 0) then
442 : call getghc(cpopt, kres_bk, cprj_dum, ghc_bk, gsc_bk, gs_hamk, gvnlxc_bk, &
443 393 : rdummy, mpi_enreg, ndat, prtvol, sij_opt, tim_getghc, type_calc0)
444 : else
445 : call prep_getghc(kres_bk, gs_hamk, gvnlxc_bk, ghc_bk, gsc_bk, rdummy, ndat, &
446 500 : mpi_enreg, prtvol, sij_opt, cpopt, cprj_dum, already_transposed=.False.)
447 : end if
448 :
449 : ! Compute (H - e_0 S) |K R_0>
450 893 : call cg_get_residvecs(usepaw, npwsp, ndat, eig(ib_start), kres_bk, ghc_bk, gsc_bk, residv_bk)
451 893 : call cg_norm2g(istwf_k, npwsp, ndat, residv_bk, lambda_bk, me_g0, comm_bsf)
452 :
453 : ! Compute lambda
454 38645 : dots_bk = zero
455 : !$OMP PARALLEL DO PRIVATE(jj, kk)
456 13381 : do idat=1,ndat
457 12488 : jj = 1 + (idat - 1) * npwsp; kk = idat * npwsp
458 : call dotprod_g(dots_bk(1,idat), dots_bk(2,idat), istwf_k, npwsp, option1, &
459 13381 : diis%chain_resv(:,:,0,idat), residv_bk(:,jj), me_g0, xmpi_comm_self)
460 : end do
461 :
462 893 : call xmpi_sum(dots_bk, comm_bsf, ierr)
463 :
464 : ! Build |Psi_1> = |Phi_0> + lambda |K R_0>
465 : !$OMP PARALLEL DO PRIVATE(jj, kk)
466 13381 : do idat=1,ndat
467 12488 : lambda_bk(idat) = -dots_bk(1,idat) / lambda_bk(idat)
468 12488 : jj = 1 + (idat - 1) * npwsp; kk = idat * npwsp
469 41321029 : cg_bk(:,jj:kk) = diis%chain_phi(:,:,0,idat) + lambda_bk(idat) * kres_bk(:,jj:kk)
470 : end do
471 :
472 : ! ===============
473 : ! DIIS iterations
474 : ! ===============
475 3215 : iter_loop: do iter=1,max_niter_block
476 :
477 2454 : if (iter > 1) then
478 : ! Solve DIIS equations and update cg_bk and residv_bk for iter > 1
479 1561 : call diis%update_block(iter, npwsp, ndat, cg_bk, residv_bk, comm_bsf)
480 :
481 : ! Precondition residual, output in kres_bk.
482 1561 : call cg_zcopy(npwsp * ndat, residv_bk, kres_bk)
483 1561 : call cg_precon_many(istwf_k, npw, my_nspinor, ndat, cg_bk, optekin, kinpw, kres_bk, me_g0, comm_bsf)
484 :
485 : ! Compute cg_bk with the same lambda(ndat) obtained at iteration #0
486 1561 : call cg_zaxpy_many_areal(npwsp, ndat, lambda_bk, kres_bk, cg_bk)
487 : end if
488 :
489 : ! Compute H |phi_now> and evaluate new enlx for NC.
490 : call getghc_eigresid(gs_hamk, npw, my_nspinor, ndat, cg_bk, ghc_bk, gsc_bk, mpi_enreg, prtvol, &
491 2454 : eig(ib_start), resid(ib_start), enlx(ib_start), residv_bk, gvnlxc_bk, normalize=.True.)
492 :
493 : ! Store new wavevefunctions and residuals.
494 2454 : call diis%push_iter(iter, ndat, eig(ib_start), resid(ib_start), enlx(ib_start), cg_bk, residv_bk, gsc_bk, "DIIS")
495 :
496 : ! CHECK FOR CONVERGENCE
497 2454 : if (diis%exit_iter(iter, ndat, max_niter_block, occ(ib_start), accuracy_ene, dtset, comm_bsf)) exit iter_loop
498 :
499 : ! Compute <R_i|R_j> and <i|S|j> for j=iter
500 3083 : if (iter /= max_niter_block) call diis%eval_mats(iter, ndat, me_g0, comm_bsf)
501 : end do iter_loop
502 :
503 2439 : if (prtvol == -level) call diis%print_block(ib_start, ndat, istep, ikpt, isppol)
504 : end do ! iblock
505 :
506 653 : call timab(1634, 2, tsec) !"rmm_diis:band_opt"
507 : if (timeit) call cwtime_report(" rmm_diis:band_opt", cpu, wall, gflops)
508 653 : call rmm_ydoc%add_tabular_line(resids2str("rmm-diis"), indent=0)
509 :
510 : ! ===============================
511 : ! Orthogonalize states after DIIS
512 : ! ===============================
513 653 : call timab(583,1,tsec) ! "vtowfk(pw_orthon)"
514 :
515 : !ortalgo = 3 !; ortalgo = mpi_enreg%paral_kgb
516 : !call pw_orthon(0, 0, istwf_k, mcg, mgsc, npwsp, nband, ortalgo, gsc, usepaw, cg, me_g0, comm_bsf)
517 :
518 : ! TODO: Merge the two routines.
519 653 : if (usepaw == 1) then
520 : !call cgtk_fixphase(cg, gsc, 0, 0, istwf_k, mcg, mgsc, mpi_enreg, nband, npwsp, usepaw)
521 : !call cgpaw_normalize(npwsp, nband, cg, gsc, istwf_k, me_g0, comm_bsf)
522 :
523 160 : call cgpaw_cholesky(npwsp, nband, cg, gsc, istwf_k, me_g0, comm_bsf, umat=umat)
524 :
525 : !call cgtk_fixphase(cg, gsc, 0, 0, istwf_k, mcg, mgsc, mpi_enreg, nband, npwsp, usepaw)
526 : !call cgpaw_normalize(npwsp, nband, cg, gsc, istwf_k, me_g0, comm_bsf)
527 : else
528 493 : call cgnc_cholesky(npwsp, nband, cg, istwf_k, me_g0, comm_bsf, use_gemm=.False., umat=umat)
529 : end if
530 :
531 653 : call timab(583,2,tsec)
532 : if (timeit) call cwtime_report(" pw_orthon ", cpu, wall, gflops)
533 :
534 : ! Recompute eigenvalues, residuals, and enlx after orthogonalization.
535 : ! This step is important to improve the convergence of the NC total energy
536 : ! and it guarantees that eigenvalues and residuals are consistent with the output wavefunctions.
537 : ! but we try to avoid it at the beginning of the SCF cycle.
538 : ! NB: In principle, one can rotate Vnl(b,b') using U^-1 from the Cholesky decomposition
539 : ! but the full Vnl matrix should be computed before the ortho step.
540 :
541 : ! Select value of after_ortho:
542 : !
543 : ! 0: return with inconsistent eigenvalues, residuals and enlx_bk to avoid final H |Psi>.
544 : ! 1: recompute enlx_bx after ortho. Return inconsistent eigens and residuals (last DIIS iteration).
545 : ! 2: fully consistent mode: execute final H|Psi> after ortho step to update enlx_bx, eigens, residuals
546 : !
547 : ! Total number of H |Psi> applications:
548 : !
549 : ! 1 for subdiago.
550 : ! 1 for preconditioned steepest descent.
551 : ! (nline - 1) for DIIS or nline if ultimate accuracy is reached.
552 : ! 1 if after_ortho > 0
553 : !
554 653 : after_ortho = 0
555 653 : if (accuracy_level >= 2) after_ortho = 1
556 653 : if (accuracy_level >= 4) after_ortho = 2
557 653 : if (after_ortho >= 1 .and. savemem == 0) after_ortho = 1
558 : ! It seems that PAW is more sensitive to after_ortho. Perhaps I can avoid the final H|phi> if accuracy_level == 1
559 : !if (usepaw == 1) after_ortho = 1
560 : !if (usepaw == 1) after_ortho = 0
561 653 : if (usepaw == 1) after_ortho = 2 ! FIXME ??
562 :
563 653 : if (after_ortho == 0) then
564 90 : call wrtout(std_out, " VERY-FAST: Won't recompute data after orthogonalization.")
565 :
566 : !else if (after_ortho == 1 .and. savemem == 0) then
567 563 : else if (after_ortho == 1 .and. savemem == 0 .and. usepaw == 0) then
568 : !if (prtvol == -level)
569 363 : call wrtout(std_out, " FAST: Recomputing data by rotating matrix elements.")
570 :
571 : if (usepaw == 0 .or. has_fock) then
572 : ! Rotate gvnlxc by solving X_new U = Y_old for X with U upper triangle.
573 : ! Compute enlx with rotated cg and gvnlxc.
574 363 : if (istwf_k == 1) then
575 267 : call ZTRSM('R', 'U', 'N', 'N', npwsp, nband, cone, umat, nband, gvnlxc, npwsp)
576 : else
577 96 : call DTRSM('R', 'U', 'N', 'N', 2*npwsp, nband, one, umat, nband, gvnlxc, 2*npwsp)
578 : end if
579 1089 : ABI_MALLOC(dots, (2, nband))
580 363 : call cg_zdotg_zip(istwf_k, npwsp, nband, option1, cg, gvnlxc, dots, me_g0, comm_bsf)
581 8387 : enlx = dots(1,:)
582 363 : ABI_FREE(dots)
583 : end if
584 :
585 : if (.False.) then
586 : !if (usepaw == 1) then
587 : ! Compute new eigenvalues, residual vectors and norms.
588 : ! Rotate ghc by solving X_new U = Y_old for X with U upper triangle.
589 : if (istwf_k == 1) then
590 : call ZTRSM('R', 'U', 'N', 'N', npwsp, nband, cone, umat, nband, ghc, npwsp)
591 : else
592 : call DTRSM('R', 'U', 'N', 'N', 2*npwsp, nband, one, umat, nband, ghc, 2*npwsp)
593 : end if
594 : ABI_MALLOC_OR_DIE(gwork, (2, npwsp*nband), ierr)
595 : call cg_get_eigens(usepaw, istwf_k, npwsp, nband, cg, ghc, gsc, eig, me_g0, comm_bsf)
596 : call cg_get_residvecs(usepaw, npwsp, nband, eig, cg, ghc, gsc, gwork)
597 : call cg_norm2g(istwf_k, npwsp, nband, gwork, resid, me_g0, comm_bsf)
598 : call rmm_ydoc%add_tabular_line(resids2str("ortho_rot"), indent=0)
599 : ABI_FREE(gwork)
600 : end if
601 :
602 : else
603 200 : if (after_ortho == 1) call wrtout(std_out, " SLOW: Recomputing enlx gvnlx by calling nonlop.")
604 200 : if (after_ortho == 2) call wrtout(std_out, " VERY-SLOW: Recomputing eigens and residues by calling getghc.")
605 :
606 496 : do iblock=1,nblocks
607 296 : igs = 1 + (iblock - 1) * npwsp * bsize; ige = min(iblock * npwsp * bsize, npwsp * nband)
608 296 : ndat = (ige - igs + 1) / npwsp
609 296 : ib_start = 1 + (iblock - 1) * bsize; ib_stop = min(iblock * bsize, nband)
610 296 : cg_bk => cg(:,igs:ige); if (usepaw == 1) gsc_bk => gsc(1:2,igs:ige)
611 296 : if (savemem == 0) then
612 256 : ghc_bk => ghc(:,igs:ige); gvnlxc_bk => gvnlxc(:,igs:ige)
613 : end if
614 :
615 200 : select case (after_ortho)
616 : case (1)
617 : ! recompute NC enlx_bx after ortho.
618 : ! eigens and residuals are inconsistent as they have been computed before pw_orthon.
619 40 : signs = 1; paw_opt = 0
620 40 : if (usepaw == 1) then
621 0 : signs = 2; paw_opt = 3
622 : end if
623 40 : if (paral_kgb == 0) then
624 : call nonlop(choice1, cpopt, cprj_dum, enlx(ib_start:), gs_hamk, 0, eig(ib_start), &
625 40 : mpi_enreg, ndat, 1, paw_opt, signs, gsc_bk, tim_nonlop, cg_bk, gvnlxc_bk)
626 : else
627 : call prep_nonlop(choice1, cpopt, cprj_dum, enlx(ib_start), gs_hamk, 0, eig(ib_start), &
628 : ndat, mpi_enreg, 1, paw_opt, signs, gsc_bk, tim_nonlop, &
629 0 : cg_bk, gvnlxc_bk, already_transposed=.False.)
630 : end if
631 :
632 : case (2)
633 : ! Consistent mode: update enlx_bx, eigens, residuals after orthogonalizalization.
634 : call getghc_eigresid(gs_hamk, npw, my_nspinor, ndat, cg_bk, ghc_bk, gsc_bk, mpi_enreg, prtvol, &
635 : eig(ib_start), resid(ib_start), enlx(ib_start), residv_bk, gvnlxc_bk, &
636 256 : normalize=usepaw == 1)
637 : case default
638 296 : ABI_BUG(sjoin("Wrong after_ortho:", itoa(after_ortho)))
639 : end select
640 : end do ! iblock
641 :
642 200 : call rmm_ydoc%add_tabular_line(resids2str("after_ortho"), indent=0)
643 : end if ! after_ortho > 0
644 :
645 : !if (usepaw == 1) then
646 : ! !call cgtk_fixphase(cg, gsc, 0, 0, istwf_k, mcg, mgsc, mpi_enreg, nband, npwsp, usepaw)
647 : ! call cg_set_imag0_to_zero(istwf_k, me_g0, npwsp, nband, cg, max_absimag)
648 : ! call cg_set_imag0_to_zero(istwf_k, me_g0, npwsp, nband, gsc, max_absimag)
649 : ! call cgpaw_normalize(npwsp, nband, cg, gsc, istwf_k, me_g0, comm_bsf)
650 : !end if
651 :
652 : if (timeit) call cwtime_report(" after_ortho ", cpu, wall, gflops)
653 :
654 : ! Revert mixprec to previous status before returning.
655 653 : if (use_fft_mixprec) prev_mixprec = fftcore_set_mixprec(prev_mixprec)
656 :
657 : !if (dtset%prtvol > 0) then
658 653 : if (diis%stats%length() > 0) call rmm_ydoc%add_dict("skip_stats", diis%stats)
659 653 : call rmm_ydoc%write_and_free(std_out)
660 :
661 : if (timeit) call cwtime_report(" rmm_diis total: ", cpu_all, wall_all, gflops_all)
662 :
663 : ! Final cleanup.
664 653 : ABI_FREE(lambda_bk)
665 653 : ABI_FREE(dots_bk)
666 653 : ABI_FREE(residv_bk)
667 653 : ABI_FREE(kres_bk)
668 653 : ABI_FREE(umat)
669 653 : if (savemem == 0) then
670 605 : ABI_FREE(ghc)
671 605 : ABI_FREE(gvnlxc)
672 : else
673 48 : ABI_FREE(ghc_bk)
674 48 : ABI_FREE(gvnlxc_bk)
675 : end if
676 1959 : call diis%free()
677 :
678 : contains
679 :
680 1506 : function resids2str(level) result(str)
681 : character(len=*),intent(in) :: level
682 : character(len=500) :: str
683 1506 : res_stats = stats_eval(resid(1:nb_pocc))
684 : !res_stats = stats_eval(resid(1:nband))
685 1506 : write(str, "(1x, a12, 4(es10.3))") trim(level), res_stats%mean, res_stats%min, res_stats%max, res_stats%stdev
686 1506 : end function resids2str
687 :
688 : end subroutine rmm_diis
689 : !!***
690 :
691 : !!****f* m_rmm_diis/rmm_diis_push_iter
692 : !! NAME
693 : !! rmm_diis_push_iter
694 : !!
695 : !! FUNCTION
696 : !! Save one iteration of the DIIS algorithm.
697 : !!
698 : !! INPUTS
699 : !!
700 : !! OUTPUT
701 : !!
702 : !! SOURCE
703 :
704 3347 : subroutine rmm_diis_push_iter(diis, iter, ndat, eig_bk, resid_bk, enlx_bk, cg_bk, residv_bk, gsc_bk, tag)
705 :
706 : class(rmm_diis_t),intent(inout) :: diis
707 : integer,intent(in) :: iter, ndat
708 : real(dp),intent(in) :: eig_bk(ndat), resid_bk(ndat), enlx_bk(ndat)
709 : real(dp),intent(in) :: cg_bk(2, diis%npwsp*ndat), residv_bk(2, diis%npwsp*ndat), gsc_bk(2, diis%npwsp*ndat*diis%usepaw)
710 : character(len=*),intent(in) :: tag
711 :
712 : !Local variables-------------------------------
713 : integer :: idat, ibk
714 : ! *************************************************************************
715 :
716 3347 : diis%last_iter = iter
717 51531 : diis%hist_ene(iter, 1:ndat) = eig_bk
718 51531 : diis%hist_resid(iter, 1:ndat) = resid_bk
719 51531 : diis%hist_enlx(iter, 1:ndat) = enlx_bk
720 51531 : diis%step_type(iter, 1:ndat) = tag
721 :
722 51531 : do idat=1,ndat
723 48184 : if (iter == 0) then
724 12488 : if (diis%cplex == 2) then
725 27804 : diis%resmat(:, 0, 0, idat) = [resid_bk(idat), zero]
726 : else
727 6440 : diis%resmat(:, 0, 0, idat) = resid_bk(idat)
728 : end if
729 : end if
730 : !write(std_out, *)"res0", diis%resmat(:, 0, 0, idat)
731 48184 : diis%step_type(iter, idat) = tag
732 48184 : ibk = 1 + (idat - 1) * diis%npwsp
733 48184 : call cg_zcopy(diis%npwsp, cg_bk(:,ibk), diis%chain_phi(:,:,iter,idat))
734 48184 : call cg_zcopy(diis%npwsp, residv_bk(:,ibk), diis%chain_resv(:,:,iter,idat))
735 51531 : if (diis%usepaw == 1) call cg_zcopy(diis%npwsp, gsc_bk(:,ibk), diis%chain_sphi(:,:,iter,idat))
736 : end do
737 :
738 3347 : end subroutine rmm_diis_push_iter
739 : !!***
740 :
741 : !!****f* m_rmm_diis/rmm_diis_exit_iter
742 : !! NAME
743 : !! rmm_diis_exit_iter
744 : !!
745 : !! FUNCTION
746 : !! Return true if we can exit the DIIS iteration
747 : !!
748 : !! INPUTS
749 : !!
750 : !! OUTPUT
751 : !!
752 : !! SOURCE
753 :
754 2454 : logical function rmm_diis_exit_iter(diis, iter, ndat, niter_block, occ_bk, accuracy_ene, dtset, comm) result(ans)
755 :
756 : class(rmm_diis_t),intent(inout) :: diis
757 : integer,intent(in) :: iter, ndat, niter_block, comm
758 : real(dp),intent(in) :: occ_bk(ndat)
759 : real(dp),intent(in) :: accuracy_ene
760 : type(dataset_type),intent(in) :: dtset
761 :
762 : !Local variables-------------------------------
763 : integer,parameter :: master = 0
764 4908 : integer :: idat, ierr, nok, checks(ndat) !nbocc,
765 : real(dp) :: resid, deltae, deold , fact
766 4908 : character(len=50) :: msg_list(ndat)
767 : ! *************************************************************************
768 :
769 2454 : diis%last_iter = iter !; ans = .False.; return
770 2454 : if (xmpi_comm_rank(comm) /= master) goto 10
771 :
772 : ! Tolerances depend on accuracy_level and occupation of the state.
773 16409 : checks = 0
774 :
775 16409 : do idat=1,ndat
776 15016 : resid = diis%hist_resid(iter, idat)
777 15016 : deold = diis%hist_ene(1, idat) - diis%hist_ene(0, idat)
778 15016 : deltae = diis%hist_ene(iter, idat) - diis%hist_ene(iter-1, idat)
779 :
780 : ! Relative criterion on eigenvalue differerence.
781 : ! Abinit default in the CG part is 0.005 that is really low (0.3 in V).
782 : ! Here we increase it depending whether the state is occupied or empty
783 15016 : fact = one !; if (dtset%iscf > 0 .and. abs(occ_bk(idat)) < diis%tol_occupied) fact = three
784 15016 : if (diis%accuracy_level == 1) fact = fact * 18
785 15016 : if (diis%accuracy_level == 2) fact = fact * 12
786 15016 : if (diis%accuracy_level == 3) fact = fact * 6
787 15016 : if (abs(deltae) < fact * dtset%tolrde * abs(deold)) then
788 1528 : checks(idat) = 1; msg_list(idat) = "deltae < fact * tolrde * deold"; cycle
789 : end if
790 :
791 14881 : if (dtset%iscf < 0) then
792 : ! This is the only condition available for NSCF run.
793 869 : if (resid < dtset%tolwfr) then
794 419 : checks(idat) = 1; msg_list(idat) = 'resid < tolwfr'; cycle
795 : end if
796 :
797 : else
798 : ! Conditions available for SCF run.
799 12619 : if (resid < dtset%tolwfr) then
800 102 : checks(idat) = 1; msg_list(idat) = 'resid < tolwfr'; cycle
801 : end if
802 :
803 : ! Absolute criterion on eigenvalue difference. Assuming error on Etot ~ band_energy.
804 12517 : fact = one; if (dtset%iscf > 0 .and. abs(occ_bk(idat)) < diis%tol_occupied) fact = ten
805 12517 : if (sqrt(abs(resid)) < fact * accuracy_ene) then
806 0 : checks(idat) = 1; msg_list(idat) = 'resid < accuracy_ene'; cycle
807 : end if
808 : end if
809 : end do ! idat
810 :
811 : ! Depending on the accuracy_level either full block or a fraction of it must pass the test in order to exit.
812 16409 : nok = count(checks /= 0)
813 1393 : if (diis%accuracy_level == 1) ans = nok >= 0.65_dp * ndat
814 1393 : if (diis%accuracy_level == 2) ans = nok >= 0.75_dp * ndat
815 1393 : if (diis%accuracy_level == 3) ans = nok >= 0.90_dp * ndat
816 1393 : if (diis%accuracy_level == 4) ans = nok == ndat
817 :
818 1393 : if (ans) then
819 : ! Log exit only if this is not the last iteration.
820 108 : if (iter /= niter_block) then
821 646 : do idat=1,ndat
822 646 : if (checks(idat) /= 0) call diis%stats%increment(msg_list(idat), 1)
823 : end do
824 : end if
825 : end if
826 :
827 : ! Broadcast final decision to all ranks.
828 2454 : 10 call xmpi_bcast(ans, master, comm, ierr)
829 :
830 2454 : end function rmm_diis_exit_iter
831 : !!***
832 :
833 : !!****f* m_rmm_diis/rmm_diis_print_block
834 : !! NAME
835 : !! rmm_diis_print_block
836 : !!
837 : !! FUNCTION
838 : !! Print energies, residuals for a block of states.
839 : !!
840 : !! INPUTS
841 : !!
842 : !! OUTPUT
843 : !!
844 : !! SOURCE
845 :
846 0 : subroutine rmm_diis_print_block(diis, ib_start, ndat, istep, ikpt, isppol)
847 :
848 : class(rmm_diis_t),intent(in) :: diis
849 : integer,intent(in) :: ib_start, ndat, istep, ikpt, isppol
850 :
851 : !Local variables-------------------------------
852 : integer :: iter, idat, iband
853 : real(dp) :: deltae, deold, dedold, absdiff
854 : character(len=500) :: msg
855 : ! *************************************************************************
856 :
857 : call wrtout(std_out, &
858 : sjoin("<BEGIN RMM-DIIS-BLOCK, istep:", itoa(istep), ", ikpt:", itoa(ikpt), ", spin: ", itoa(isppol), ">"), &
859 0 : pre_newlines=1)
860 :
861 0 : do idat=1,ndat
862 : write(msg,'(1a, 2(a5), 4(a14), 1x, a6)') &
863 0 : "#", 'iter', "band", "eigen_eV", "eigde_meV", "de/dold", "resid", "type"; call wrtout(std_out, msg)
864 :
865 0 : iband = ib_start + idat - 1
866 0 : deold = diis%hist_ene(1, idat) - diis%hist_ene(0, idat)
867 0 : do iter=0,diis%last_iter
868 0 : dedold = zero; absdiff = zero
869 0 : if (iter > 0) then
870 0 : deltae = diis%hist_ene(iter, idat) - diis%hist_ene(iter-1, idat)
871 0 : dedold = deltae / deold
872 0 : absdiff = (diis%hist_ene(iter, idat) - diis%hist_ene(iter-1, idat))
873 : end if
874 :
875 : write(msg,"(1x, 2(i5), 4(es14.6), 1x, a6)") &
876 0 : iter, iband, diis%hist_ene(iter, idat) * Ha_eV, absdiff * Ha_meV, dedold, &
877 0 : diis%hist_resid(iter, idat), diis%step_type(iter, idat); call wrtout(std_out, msg)
878 : end do
879 : end do
880 :
881 0 : call wrtout(std_out, "<END RMM-DIIS-BLOCK>", newlines=1)
882 :
883 0 : end subroutine rmm_diis_print_block
884 : !!***
885 :
886 : !!****f* m_rmm_diis/getghc_eigresid
887 : !! NAME
888 : !! getghc_eigresid
889 : !!
890 : !! FUNCTION
891 : !! Compute new eigenvalues, residuals, H |psi> and enlx from cg and gsc.
892 : !!
893 : !! INPUTS
894 : !!
895 : !! OUTPUT
896 : !!
897 : !! SOURCE
898 :
899 2758 : subroutine getghc_eigresid(gs_hamk, npw, my_nspinor, ndat, cg, ghc, gsc, mpi_enreg, prtvol, &
900 2758 : eig, resid, enlx, residvecs, gvnlxc, normalize)
901 :
902 : !Arguments ------------------------------------
903 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
904 : integer,intent(in) :: npw, my_nspinor, ndat, prtvol
905 : real(dp),intent(inout) :: cg(2, npw*my_nspinor*ndat)
906 : real(dp),intent(out) :: ghc(2,npw*my_nspinor*ndat), gsc(2,npw*my_nspinor*ndat*gs_hamk%usepaw)
907 : type(mpi_type),intent(in) :: mpi_enreg
908 : real(dp),intent(out) :: eig(ndat), resid(ndat), enlx(ndat)
909 : real(dp),intent(out) :: residvecs(2, npw*my_nspinor*ndat), gvnlxc(2, npw*my_nspinor*ndat)
910 : logical,optional,intent(in) :: normalize
911 :
912 : !Local variables-------------------------------
913 : integer,parameter :: type_calc0 = 0, option1 = 1, option2 = 2, tim_getghc = 0
914 : integer :: istwf_k, usepaw, cpopt, sij_opt, npwsp, me_g0, comm_bsf
915 : real(dp),parameter :: rdummy = zero
916 : !real(dp) :: cpu, wall, gflops
917 : logical :: normalize_, has_fock
918 : !arrays
919 2758 : real(dp) :: dots(2, ndat)
920 22064 : type(pawcprj_type) :: cprj_dum(1,1)
921 : ! *************************************************************************
922 :
923 : !if (timeit) call cwtime(cpu, wall, gflops, "start")
924 2758 : normalize_ = .True.; if (present(normalize)) normalize_ = normalize
925 2758 : npwsp = npw * my_nspinor
926 2758 : usepaw = gs_hamk%usepaw; istwf_k = gs_hamk%istwf_k; me_g0 = mpi_enreg%me_g0
927 2758 : comm_bsf = mpi_enreg%comm_spinorfft; if (mpi_enreg%paral_kgb == 1) comm_bsf = mpi_enreg%comm_bandspinorfft
928 2758 : has_fock = associated(gs_hamk%fockcommon)
929 :
930 2758 : cpopt = -1; sij_opt = 0
931 2758 : if (usepaw == 1) then
932 991 : sij_opt = 1 ! matrix elements <G|S|C> have to be computed in gsc in addition to ghc
933 : cpopt = -1 ! <p_lmn|in> (and derivatives) are computed here (and not saved)
934 : end if
935 :
936 : ! NC normalization.
937 2758 : if (usepaw == 0 .and. normalize_) call cgnc_normalize(npwsp, ndat, cg, istwf_k, me_g0, comm_bsf)
938 :
939 : ! Compute H |cg>
940 : !call fock_set_ieigen(gs_hamk%fockcommon, iband)
941 2758 : if (mpi_enreg%paral_kgb == 0) then
942 : call getghc(cpopt, cg, cprj_dum, ghc, gsc, gs_hamk, gvnlxc, &
943 1210 : rdummy, mpi_enreg, ndat, prtvol, sij_opt, tim_getghc, type_calc0)
944 : else
945 : call prep_getghc(cg, gs_hamk, gvnlxc, ghc, gsc, rdummy, ndat, &
946 1548 : mpi_enreg, prtvol, sij_opt, cpopt, cprj_dum, already_transposed=.False.)
947 : end if
948 :
949 : ! PAW normalization must be done here once gsc is known.
950 2758 : if (usepaw == 1 .and. normalize_) call cgpaw_normalize(npwsp, ndat, cg, gsc, istwf_k, me_g0, comm_bsf)
951 :
952 : ! Compute new eigenvalues, residual vectors and norms.
953 2758 : call cg_get_eigens(usepaw, istwf_k, npwsp, ndat, cg, ghc, gsc, eig, me_g0, comm_bsf)
954 2758 : call cg_get_residvecs(usepaw, npwsp, ndat, eig, cg, ghc, gsc, residvecs)
955 2758 : call cg_norm2g(istwf_k, npwsp, ndat, residvecs, resid, me_g0, comm_bsf)
956 :
957 2758 : if (usepaw == 0 .or. has_fock) then
958 : ! Evaluate new enlx from gvnlxc.
959 1767 : call cg_zdotg_zip(istwf_k, npwsp, ndat, option1, cg, gvnlxc, dots, me_g0, comm_bsf)
960 33150 : enlx = dots(1,:)
961 : end if
962 :
963 : !if (timeit) call cwtime_report(" getghc_eigresid", cpu, wall, gflops)
964 :
965 8274 : end subroutine getghc_eigresid
966 : !!***
967 :
968 : !!****f* m_rmm_diis/rmm_diis_new
969 : !! NAME
970 : !! rmm_diis_new
971 : !!
972 : !! FUNCTION
973 : !! Build new rmm_diis_t instance.
974 : !!
975 : !! INPUTS
976 : !!
977 : !! OUTPUT
978 : !!
979 : !! SOURCE
980 :
981 653 : type(rmm_diis_t) function rmm_diis_new(accuracy_level, usepaw, istwf_k, npwsp, max_niter, bsize, prtvol) result(diis)
982 :
983 : !Arguments ------------------------------------
984 : integer,intent(in) :: accuracy_level, usepaw, istwf_k, npwsp, max_niter, bsize, prtvol
985 : ! *************************************************************************
986 :
987 653 : diis%accuracy_level = accuracy_level
988 653 : diis%usepaw = usepaw
989 653 : diis%istwf_k = istwf_k
990 126 : diis%cplex = 2; if (istwf_k == 2) diis%cplex = 1
991 653 : diis%npwsp = npwsp
992 653 : diis%max_niter = max_niter
993 653 : diis%bsize = bsize
994 653 : diis%prtvol = prtvol
995 :
996 2612 : ABI_MALLOC(diis%hist_ene, (0:max_niter+2, bsize))
997 1959 : ABI_MALLOC(diis%hist_resid, (0:max_niter+2, bsize))
998 78620 : ABI_CALLOC(diis%hist_enlx, (0:max_niter+2, bsize))
999 2612 : ABI_MALLOC(diis%step_type, (0:max_niter+2, bsize))
1000 3265 : ABI_MALLOC(diis%chain_phi, (2, npwsp, 0:max_niter, bsize))
1001 3265 : ABI_MALLOC(diis%chain_sphi, (2, npwsp*usepaw, 0:max_niter, bsize))
1002 2612 : ABI_MALLOC(diis%chain_resv, (2, npwsp, 0:max_niter, bsize))
1003 560566 : ABI_CALLOC(diis%resmat, (diis%cplex, 0:max_niter, 0:max_niter, bsize)) ! <R_i|R_j>
1004 :
1005 653 : end function rmm_diis_new
1006 : !!***
1007 :
1008 : !!****f* m_rmm_diis/rmm_diis_free
1009 : !! NAME
1010 : !! rmm_diis_free
1011 : !!
1012 : !! FUNCTION
1013 : !! Free dynamic memory.
1014 : !!
1015 : !! INPUTS
1016 : !!
1017 : !! OUTPUT
1018 : !!
1019 : !! SOURCE
1020 :
1021 653 : subroutine rmm_diis_free(diis)
1022 :
1023 : !Arguments ------------------------------------
1024 : class(rmm_diis_t),intent(inout) :: diis
1025 : ! *************************************************************************
1026 :
1027 653 : ABI_SFREE(diis%hist_ene)
1028 653 : ABI_SFREE(diis%hist_resid)
1029 653 : ABI_SFREE(diis%hist_enlx)
1030 653 : ABI_SFREE(diis%step_type)
1031 653 : ABI_SFREE(diis%chain_phi)
1032 653 : ABI_SFREE(diis%chain_sphi)
1033 653 : ABI_SFREE(diis%chain_resv)
1034 653 : ABI_SFREE(diis%resmat)
1035 :
1036 653 : call diis%stats%free()
1037 :
1038 653 : end subroutine rmm_diis_free
1039 : !!***
1040 :
1041 : !!****f* m_rmm_diis/rmm_diis_update_block
1042 : !! NAME
1043 : !! rmm_diis_update_block
1044 : !!
1045 : !! FUNCTION
1046 : !! Compute new trial wavefunctions and residuals from the DIIS chain.
1047 : !!
1048 : !! INPUTS
1049 : !!
1050 : !! OUTPUT
1051 : !!
1052 : !! SOURCE
1053 :
1054 1561 : subroutine rmm_diis_update_block(diis, iter, npwsp, ndat, cg_bk, residv_bk, comm)
1055 :
1056 : !Arguments ------------------------------------
1057 : class(rmm_diis_t),intent(in) :: diis
1058 : integer,intent(in) :: iter, npwsp, comm, ndat
1059 : real(dp),intent(inout) :: cg_bk(2, npwsp, ndat), residv_bk(2, npwsp, ndat)
1060 :
1061 : !local variables
1062 : integer,parameter :: master = 0
1063 : integer :: cplex, ierr, nprocs, my_rank, idat
1064 : real(dp) :: noise !, cpu, wall, gflops
1065 1561 : real(dp),allocatable :: wmat1(:,:,:), wvec(:,:,:), alphas(:,:)
1066 : character(len=500) :: msg
1067 : ! *************************************************************************
1068 :
1069 : !if (timeit) call cwtime(cpu, wall, gflops, "start")
1070 1561 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1071 1561 : cplex = diis%cplex
1072 :
1073 : ! Solve system of linear equations.
1074 : ! Only master works so that we are sure we have the same solution.
1075 252929 : ABI_CALLOC(wvec, (cplex, 0:iter, ndat))
1076 :
1077 1561 : if (my_rank == master) then
1078 4355 : ABI_MALLOC(wmat1, (cplex, 0:iter, 0:iter))
1079 :
1080 10447 : do idat=1,ndat
1081 : !if (mod(idat, nprocs) /= my_rank) cycle ! MPI parallelism
1082 9576 : wvec(1, iter, idat) = -one
1083 382676 : wmat1 = zero
1084 43224 : wmat1(1,:,iter) = -one
1085 43224 : wmat1(1,iter,:) = -one
1086 9576 : wmat1(1,iter,iter) = zero
1087 211976 : wmat1(:,0:iter-1, 0:iter-1) = diis%resmat(:, 0:iter-1, 0:iter-1, idat)
1088 :
1089 9576 : call xhesv_cplex("U", cplex, iter+1, 1, wmat1, wvec(:,:,idat), msg, ierr)
1090 9576 : ABI_CHECK(ierr == 0, msg)
1091 :
1092 : !if (diis%prtvol == -level) then
1093 : ! write(std_out,*)"wvec:", wvec(:,:,idat)
1094 : ! write(std_out,*)"sum(wvec):", sum(wvec(:, 0:iter-1, idat), dim=2)
1095 : !end if
1096 10447 : if (cplex == 2) then
1097 : ! coefficients should sum up to 1 but sometimes we get a small imaginary part. here we remove it
1098 26628 : noise = sum(wvec(2, 0:iter-1, idat))
1099 26628 : wvec(2, 0:iter-1, idat) = wvec(2, 0:iter-1, idat) - noise * iter
1100 : end if
1101 :
1102 : end do
1103 871 : ABI_FREE(wmat1)
1104 : end if
1105 :
1106 : ! Master broadcasts data.
1107 1561 : if (nprocs > 1) call xmpi_bcast(wvec, master, comm, ierr)
1108 : !if (nprocs > 1) call xmpi_sum(wvec, comm, ierr)
1109 :
1110 : ! Take linear combination of chain_phi and chain_resv.
1111 : !$omp parallel private(idat, alphas)
1112 1561 : if (cplex /= 2) then
1113 894 : ABI_MALLOC(alphas, (1, 0:iter))
1114 : end if
1115 :
1116 : !$omp do
1117 24769 : do idat=1,ndat
1118 24769 : if (cplex == 2) then
1119 16932 : call cg_zgemv("N", npwsp, iter, diis%chain_phi(:,:,:,idat), wvec(:,:,idat), cg_bk(:,:,idat))
1120 16932 : call cg_zgemv("N", npwsp, iter, diis%chain_resv(:,:,:,idat), wvec(:,:,idat), residv_bk(:,:,idat))
1121 : else
1122 : ! coefficients are real --> use DGEMV.
1123 28368 : alphas(1,:) = wvec(1,:,idat)
1124 6276 : call dgemv("N", 2*npwsp, iter, one, diis%chain_phi(:,:,:,idat), 2*npwsp, alphas, 1, zero, cg_bk(:,:,idat), 1)
1125 6276 : call dgemv("N", 2*npwsp, iter, one, diis%chain_resv(:,:,:,idat), 2*npwsp, alphas, 1, zero, residv_bk(:,:,idat), 1)
1126 : end if
1127 : end do ! idat
1128 : !$omp end do
1129 :
1130 1561 : if (cplex /= 2) then
1131 298 : ABI_FREE(alphas)
1132 : end if
1133 : !$omp end parallel
1134 :
1135 1561 : ABI_FREE(wvec)
1136 : !if (timeit) call cwtime_report(" update_block", cpu, wall, gflops)
1137 :
1138 1561 : end subroutine rmm_diis_update_block
1139 : !!***
1140 :
1141 : !!****f* m_rmm_diis/rmm_diis_eval_mats
1142 : !! NAME
1143 : !! rmm_diis_eval_mats
1144 : !!
1145 : !! FUNCTION
1146 : !! Compute matrix elements required by the RMM-DIIS method.
1147 : !!
1148 : !! INPUTS
1149 : !!
1150 : !! OUTPUT
1151 : !!
1152 : !! SOURCE
1153 :
1154 1561 : subroutine rmm_diis_eval_mats(diis, iter, ndat, me_g0, comm)
1155 :
1156 : !Arguments ------------------------------------
1157 : class(rmm_diis_t),intent(inout) :: diis
1158 : integer,intent(in) :: iter, ndat, me_g0, comm
1159 :
1160 : !local variables
1161 : integer :: ii, ierr, idat, nprocs, option
1162 : real(dp) :: dotr, doti !, cpu, wall, gflops
1163 : !integer :: requests(ndat)
1164 : ! *************************************************************************
1165 :
1166 : !if (timeit) call cwtime(cpu, wall, gflops, "start")
1167 1561 : nprocs = xmpi_comm_size(comm)
1168 1561 : option = 2; if (diis%cplex == 1) option = 1
1169 :
1170 24769 : do idat=1,ndat
1171 :
1172 82897 : do ii=0,iter
1173 : ! <R_i|R_j>
1174 : call dotprod_g(dotr, doti, diis%istwf_k, diis%npwsp, option, &
1175 58128 : diis%chain_resv(:,:,ii,idat), diis%chain_resv(:,:,iter,idat), me_g0, xmpi_comm_self)
1176 58128 : if (ii == iter) doti = zero
1177 81336 : if (diis%cplex == 2) then
1178 126936 : diis%resmat(:, ii, iter, idat) = [dotr, doti]
1179 : else
1180 15816 : diis%resmat(1, ii, iter, idat) = dotr
1181 : end if
1182 : end do ! ii
1183 :
1184 : !if (nprocs > 1) then
1185 : ! call xmpi_sum(diis%resmat(:,0:iter,iter,idat), comm, ierr)
1186 : ! !call xmpi_isum_ip(diis%resmat(:,0:iter,iter,idat), comm, requests(idat), ierr)
1187 : !endif
1188 : !if (diis%prtvol == -level) write(std_out,*)"iter, idat, resmat:", iter, idat, diis%resmat(:,0:iter,iter,idat)
1189 : end do ! idat
1190 :
1191 1561 : if (nprocs > 1) call xmpi_sum(diis%resmat(:,0:iter,iter,1:ndat), comm, ierr)
1192 : !if (nprocs > 1) call xmpi_waitall(requests, ierr)
1193 : !if (timeit) call cwtime_report(" eval_mats", cpu, wall, gflops)
1194 :
1195 1561 : end subroutine rmm_diis_eval_mats
1196 : !!***
1197 :
1198 : !!****f* m_numeric_tools/my_pack_matrix
1199 : !! NAME
1200 : !! my_pack_matrix
1201 : !!
1202 : !! FUNCTION
1203 : !! Packs a matrix into hermitian format
1204 : !!
1205 : !! INPUTS
1206 : !! N: size of matrix
1207 : !! cplx: 2 if matrix is complex, 1 for real matrix.
1208 : !! mat_in(cplx, N*N)= matrix to be packed
1209 : !!
1210 : !! OUTPUT
1211 : !! mat_out(cplx*N*N+1/2)= packed matrix (upper triangle)
1212 : !!
1213 : !! SOURCE
1214 :
1215 160 : subroutine my_pack_matrix(n, mat_in, mat_out)
1216 :
1217 : !Arguments ------------------------------------
1218 : integer, intent(in) :: N
1219 : real(dp), intent(in) :: mat_in(N, N)
1220 : real(dp), intent(out) :: mat_out(2, N*(N+1)/2)
1221 :
1222 : !local variables
1223 : integer :: isubh, i, j
1224 : ! *************************************************************************
1225 :
1226 160 : isubh = 1
1227 4348 : do j=1,N
1228 65134 : do i=1,j
1229 60786 : mat_out(1,isubh) = mat_in(i, j)
1230 60786 : mat_out(2,isubh) = zero
1231 64974 : isubh = isubh + 1
1232 : end do
1233 : end do
1234 :
1235 160 : end subroutine my_pack_matrix
1236 : !!***
1237 :
1238 : !!****f* ABINIT/subspace_rotation
1239 : !! NAME
1240 : !! subspace_rotation
1241 : !!
1242 : !! FUNCTION
1243 : !! This routine computes the <i|H|j> matrix elements and then performs the subspace rotation
1244 : !! of the orbitals (Rayleigh-Ritz procedure)
1245 : !! The main difference with respect to other similar routines is that this implementation does not require
1246 : !! the <i|H|j> matrix elements as input so it can be used before starting the wavefunction optimation
1247 : !! as required e.g. by the RMM-DIIS method.
1248 : !! Moreover, the routine computes the new residuals after the subspace rotation by rotating the
1249 : !! matrix elements of the Hamiltonian in the new basis (requires more memory but client code
1250 : !! can avoid calling getghc after subspace_rotation.
1251 : !!
1252 : !! INPUTS
1253 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the hamiltonian at k
1254 : !! ptrvol
1255 : !! mpi_enreg=information about MPI parallelization
1256 : !! nband=number of bands at this k point and spin
1257 : !! npw=number of plane waves at this k point
1258 : !! my_nspinor=number of spinor components treated by this MPI proc.
1259 : !!
1260 : !! OUTPUT
1261 : !! eig(nband): New eigvalues from subspace rotation.
1262 : !! If usepaw==1:
1263 : !! gsc(2,*)=<g|S|c> matrix elements (S=overlap)
1264 : !! enlx(nband)=contribution from each band to nonlocal psp + potential Fock ACE part of total energy, at this k-point
1265 : !! ghc
1266 : !!
1267 : !! SIDE EFFECTS
1268 : !! cg(2,*)=updated wavefunctions
1269 : !! gsc(2,*)=update <G|S|C>
1270 : !!
1271 : !! SOURCE
1272 :
1273 657 : subroutine subspace_rotation(gs_hamk, prtvol, mpi_enreg, nband, npw, my_nspinor, savemem, enlx, eig, cg, gsc, ghc, gvnlxc)
1274 :
1275 : !Arguments ------------------------------------
1276 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
1277 : integer,intent(in) :: prtvol, nband, npw, my_nspinor, savemem
1278 : type(mpi_type),intent(in) :: mpi_enreg
1279 : real(dp),target,intent(inout) :: cg(2,npw*my_nspinor*nband)
1280 : real(dp),target,intent(inout) :: gsc(2,npw*my_nspinor*nband*gs_hamk%usepaw)
1281 : !real(dp),target,intent(out) :: ghc(2,npw*my_nspinor*nband), gvnlxc(2,npw*my_nspinor*nband)
1282 : real(dp),target,allocatable,intent(out) :: ghc(:,:), gvnlxc(:,:)
1283 : real(dp),intent(out) :: eig(nband), enlx(nband)
1284 :
1285 : !Local variables-------------------------------
1286 : integer,parameter :: type_calc0 = 0, tim_getghc = 0, use_subovl0 = 0, option1 = 1
1287 : integer :: ig, ig0, ib, ierr, bsize, nblocks, iblock, npwsp, ndat, ib_start, ib_stop, paral_kgb
1288 : integer :: iband, cpopt, sij_opt, igs, ige, mcg, mgsc, istwf_k, usepaw, me_g0, cplex, comm_bsf
1289 : logical :: has_fock
1290 : real(dp),parameter :: rdummy = zero
1291 : real(dp) :: cpu, wall, gflops
1292 : !arrays
1293 : real(dp),target :: fake_gsc_bk(0,0)
1294 657 : real(dp) :: subovl(use_subovl0), dots(2, nband)
1295 657 : real(dp),allocatable :: subham(:), h_ij(:,:,:), evec(:,:,:), evec_re(:,:), gwork(:,:)
1296 657 : real(dp),contiguous, pointer :: ghc_bk(:,:), gvnlxc_bk(:,:), gsc_bk(:,:)
1297 4599 : type(pawcprj_type) :: cprj_dum(1,1)
1298 : ! *************************************************************************
1299 :
1300 : if (timeit) call cwtime(cpu, wall, gflops, "start")
1301 :
1302 657 : usepaw = gs_hamk%usepaw; istwf_k = gs_hamk%istwf_k
1303 657 : paral_kgb = mpi_enreg%paral_kgb; me_g0 = mpi_enreg%me_g0
1304 428 : comm_bsf = mpi_enreg%comm_spinorfft; if (mpi_enreg%paral_kgb == 1) comm_bsf = mpi_enreg%comm_bandspinorfft
1305 657 : npwsp = npw * my_nspinor
1306 657 : has_fock = associated(gs_hamk%fockcommon)
1307 :
1308 : ! =======================================
1309 : ! Apply H to input cg to compute <i|H|j>
1310 : ! =======================================
1311 657 : gsc_bk => fake_gsc_bk
1312 657 : cpopt = -1; sij_opt = 0
1313 657 : if (usepaw == 1) then
1314 160 : sij_opt = 1 ! matrix elements <G|S|C> have to be computed in gsc in addition to ghc
1315 : cpopt = -1 ! <p_lmn|in> (and derivatives) are computed here (and not saved)
1316 : end if
1317 :
1318 : ! Treat states in groups of bsize bands even when paral_kgb = 0
1319 657 : bsize = 8; if (paral_kgb == 1) bsize = mpi_enreg%nproc_band * mpi_enreg%bandpp
1320 657 : nblocks = nband / bsize; if (mod(nband, bsize) /= 0) nblocks = nblocks + 1
1321 :
1322 657 : cplex = 2; if (istwf_k /= 1) cplex = 1
1323 : !cplex = 2; if (istwf_k == 2) cplex = 1
1324 :
1325 840085 : ABI_CALLOC(h_ij, (cplex, nband, nband))
1326 :
1327 : ! Allocate full ghc and gvnlxc to be able to rotate residuals and Vnlx matrix elements
1328 : ! after subdiago. More memory but we can save a call to H|psi>.
1329 657 : if (savemem == 0) then
1330 1815 : ABI_MALLOC_OR_DIE(ghc, (2, npwsp*nband), ierr)
1331 1210 : ABI_MALLOC_OR_DIE(gvnlxc, (2, npwsp*nband), ierr)
1332 52 : else if (savemem == 1) then
1333 156 : ABI_MALLOC(ghc_bk, (2, npwsp*bsize))
1334 104 : ABI_MALLOC(gvnlxc_bk, (2, npwsp*bsize))
1335 : else
1336 0 : ABI_ERROR(sjoin("Invalid savemem:", itoa(savemem)))
1337 : end if
1338 :
1339 1566 : do iblock=1,nblocks
1340 909 : igs = 1 + (iblock - 1) * npwsp * bsize; ige = min(iblock * npwsp * bsize, npwsp * nband)
1341 909 : ndat = (ige - igs + 1) / npwsp
1342 909 : ib_start = 1 + (iblock - 1) * bsize; ib_stop = min(iblock * bsize, nband)
1343 :
1344 909 : if (usepaw == 1) gsc_bk => gsc(1:2,igs:ige)
1345 909 : if (savemem == 0) then
1346 845 : ghc_bk => ghc(:, igs:ige); gvnlxc_bk => gvnlxc(:, igs:ige)
1347 : end if
1348 :
1349 909 : if (paral_kgb == 0) then
1350 : call getghc(cpopt, cg(:,igs:ige), cprj_dum, ghc_bk, gsc_bk, gs_hamk, gvnlxc_bk, &
1351 409 : rdummy, mpi_enreg, ndat, prtvol, sij_opt, tim_getghc, type_calc0)
1352 : else
1353 : call prep_getghc(cg(:,igs:ige), gs_hamk, gvnlxc_bk, ghc_bk, gsc_bk, rdummy, ndat, &
1354 500 : mpi_enreg, prtvol, sij_opt, cpopt, cprj_dum, already_transposed=.False.)
1355 : end if
1356 :
1357 : ! Compute <i|H|j> for i=1,nband and all j in block
1358 909 : if (cplex == 2) then
1359 629 : call cg_zgemm("C", "N", npwsp, nband, ndat, cg, ghc_bk, h_ij(:,:,ib_start))
1360 : else
1361 280 : call dgemm("T", "N", nband, ndat, 2*npwsp, one, cg, 2*npwsp, ghc_bk, 2*npwsp, zero, h_ij(:,:,ib_start), nband)
1362 : end if
1363 :
1364 1566 : if (istwf_k /= 1) then
1365 4468 : do iband=ib_start, ib_stop
1366 238956 : h_ij(:,:,iband) = two * h_ij(:,:,iband)
1367 :
1368 4188 : if (istwf_k == 2 .and. me_g0 == 1) then
1369 : ! Gamma k-point and I have G=0. Remove double counting term.
1370 1084 : ig = 1 + (iband - ib_start) * npwsp
1371 31188 : do ib=1,nband
1372 30104 : ig0 = 1 + npwsp * (ib - 1)
1373 : #if defined FC_NVHPC
1374 : if (ig<0) write(100,*) ig,ig0,h_ij(1,ib,iband),cg(1,ig0),ghc_bk(1,ig)
1375 : #endif
1376 31188 : h_ij(1,ib,iband) = h_ij(1,ib,iband) - cg(1,ig0) * ghc_bk(1,ig)
1377 : end do
1378 : end if
1379 :
1380 : ! Force real matrix.
1381 4468 : if (cplex == 2) h_ij(2,:,iband) = zero
1382 : end do
1383 : end if
1384 :
1385 : end do ! iblock
1386 :
1387 : ! Pack <i|H|j> to prepare call to subdiago.
1388 1971 : ABI_MALLOC(subham, (nband*(nband+1)))
1389 657 : if (cplex == 2) then
1390 8925 : do iband=1,nband
1391 8925 : h_ij(2,iband,iband) = zero ! Force diagonal elements to be real
1392 : end do
1393 : end if
1394 :
1395 : if (cplex == 2) then
1396 497 : call pack_matrix(h_ij, subham, nband, 2)
1397 : else
1398 160 : call my_pack_matrix(nband, h_ij, subham)
1399 : end if
1400 :
1401 657 : ABI_FREE(h_ij)
1402 657 : call xmpi_sum(subham, comm_bsf, ierr)
1403 : if (timeit) call cwtime_report(" subspace build Hij", cpu, wall, gflops)
1404 :
1405 : ! ========================
1406 : ! Subspace diagonalization
1407 : ! =======================
1408 : ! Rotate cg, gsc and compute new eigenvalues.
1409 2628 : ABI_MALLOC(evec, (2, nband, nband))
1410 657 : mcg = npwsp * nband; mgsc = npwsp * nband * usepaw
1411 : call subdiago(cg, eig, evec, gsc, 0, 0, istwf_k, mcg, mgsc, nband, npw, my_nspinor, paral_kgb, &
1412 657 : subham, subovl, use_subovl0, usepaw, me_g0)
1413 :
1414 657 : ABI_FREE(subham)
1415 : if (timeit) call cwtime_report(" subspace subdiago", cpu, wall, gflops)
1416 :
1417 657 : if (savemem == 0) then
1418 : ! Rotate ghc matrix in the new subspace:
1419 : !
1420 : ! new_{g,b} = old_{g,i} evec_{i,b}
1421 : !
1422 : ! cg and PAW gsc have been already rotated in subdiago
1423 : !
1424 605 : if (cplex == 1) then
1425 : ! Eigenvectors are real.
1426 624 : ABI_MALLOC(evec_re, (nband, nband))
1427 117660 : evec_re = evec(1,:,:)
1428 : end if
1429 :
1430 1815 : ABI_MALLOC_OR_DIE(gwork, (2, npwsp*nband), ierr)
1431 605 : if (cplex == 1) then
1432 156 : call DGEMM("N", "N", 2*npwsp, nband, nband, one, ghc, 2*npwsp, evec_re, nband, zero, gwork, 2*npwsp)
1433 : else
1434 449 : call abi_zgemm_2r("N", "N", npwsp, nband, nband, cone, ghc, npwsp, evec, nband, czero, gwork, npwsp)
1435 : end if
1436 605 : call cg_zcopy(npwsp * nband, gwork, ghc)
1437 :
1438 : ! Rotate <G|Vnlx|Psi_n> and evaluate new enlx for NC.
1439 605 : if (usepaw == 0 .or. has_fock) then
1440 445 : if (cplex == 1) then
1441 124 : call DGEMM("N", "N", 2*npwsp, nband, nband, one, gvnlxc, 2*npwsp, evec_re, nband, zero, gwork, 2*npwsp)
1442 : else
1443 321 : call abi_zgemm_2r("N", "N", npwsp, nband, nband, cone, gvnlxc, npwsp, evec, nband, czero, gwork, npwsp)
1444 : end if
1445 : !call abi_xgemm('N','N', vectsize, nband, nband, cone, gvnlxc, vectsize, evec, nband, czero, gwork, vectsize, x_cplx=cplx)
1446 445 : call cg_zcopy(npwsp * nband, gwork, gvnlxc)
1447 445 : call cg_zdotg_zip(istwf_k, npwsp, nband, option1, cg, gvnlxc, dots, me_g0, comm_bsf)
1448 10757 : enlx = dots(1,:)
1449 : end if
1450 605 : ABI_FREE(gwork)
1451 605 : ABI_SFREE(evec_re)
1452 : end if
1453 :
1454 657 : ABI_FREE(evec)
1455 657 : if (savemem == 1) then
1456 52 : ABI_FREE(ghc_bk)
1457 52 : ABI_FREE(gvnlxc_bk)
1458 : end if
1459 :
1460 : if (timeit) call cwtime_report(" subspace final rotation", cpu, wall, gflops)
1461 :
1462 1314 : end subroutine subspace_rotation
1463 : !!***
1464 :
1465 0 : end module m_rmm_diis
1466 : !!***
|
