Branch data 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 : : !!***
|
