(git:98357aa)
Loading...
Searching...
No Matches
pexsi_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Methods using the PEXSI library to calculate the density matrix and
10!> related quantities using the Kohn-Sham and overlap matrices from the
11!> linear scaling quickstep SCF environment.
12!> \par History
13!> 2014.11 created [Patrick Seewald]
14!> \author Patrick Seewald
15! **************************************************************************************************
16
17MODULE pexsi_methods
18 USE arnoldi_api, ONLY: arnoldi_env_type,&
19 arnoldi_ev,&
20 deallocate_arnoldi_env,&
21 get_selected_ritz_val,&
22 get_selected_ritz_vector,&
23 set_arnoldi_initial_vector,&
24 setup_arnoldi_env
25 USE cp_dbcsr_api, ONLY: &
26 dbcsr_convert_csr_to_dbcsr, dbcsr_convert_dbcsr_to_csr, dbcsr_copy, dbcsr_create, &
27 dbcsr_csr_create, dbcsr_csr_create_from_dbcsr, dbcsr_csr_destroy, &
28 dbcsr_csr_eqrow_floor_dist, dbcsr_csr_print_sparsity, dbcsr_csr_type_real_8, &
29 dbcsr_desymmetrize, dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_info, &
30 dbcsr_has_symmetry, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, &
31 dbcsr_type_no_symmetry
32 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_to_csr_screening
33 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
34 USE cp_log_handling, ONLY: cp_get_default_logger,&
35 cp_logger_get_default_unit_nr,&
36 cp_logger_type
37 USE dm_ls_scf_qs, ONLY: matrix_ls_to_qs
38 USE dm_ls_scf_types, ONLY: ls_scf_env_type
39 USE input_section_types, ONLY: section_vals_type,&
40 section_vals_val_get
41 USE kinds, ONLY: dp,&
42 int_4,&
43 int_8
44 USE pexsi_interface, ONLY: cp_pexsi_dft_driver,&
45 cp_pexsi_get_options,&
46 cp_pexsi_load_real_hs_matrix,&
47 cp_pexsi_retrieve_real_dft_matrix,&
48 cp_pexsi_set_default_options,&
49 cp_pexsi_set_options
50 USE pexsi_types, ONLY: convert_nspin_cp2k_pexsi,&
51 cp2k_to_pexsi,&
52 lib_pexsi_env,&
53 pexsi_to_cp2k
54 USE qs_energy_types, ONLY: qs_energy_type
55 USE qs_environment_types, ONLY: get_qs_env,&
56 qs_environment_type
57 USE qs_ks_types, ONLY: qs_ks_env_type
58#include "./base/base_uses.f90"
59
60 IMPLICIT NONE
61 PRIVATE
62
63 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pexsi_methods'
64
65 LOGICAL, PARAMETER, PRIVATE :: careful_mod = .false.
66
67 PUBLIC :: density_matrix_pexsi, pexsi_init_read_input, pexsi_to_qs, pexsi_init_scf, pexsi_finalize_scf, &
68 pexsi_set_convergence_tolerance
69
70CONTAINS
71
72! **************************************************************************************************
73!> \brief Read CP2K input section PEXSI and pass it to the PEXSI environment
74!> \param pexsi_section ...
75!> \param pexsi_env ...
76!> \par History
77!> 11.2014 created [Patrick Seewald]
78!> \author Patrick Seewald
79! **************************************************************************************************
80 SUBROUTINE pexsi_init_read_input(pexsi_section, pexsi_env)
81 TYPE(section_vals_type), INTENT(IN), POINTER :: pexsi_section
82 TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
83
84 INTEGER :: isinertiacount_int, maxpexsiiter, &
85 min_ranks_per_pole, npsymbfact, &
86 numpole, ordering, rowordering, &
87 verbosity
88 LOGICAL :: csr_screening, isinertiacount
89 REAL(kind=dp) :: gap, muinertiaexpansion, muinertiatolerance, mumax0, mumin0, &
90 mupexsisafeguard, numelectroninitialtolerance, numelectrontargettolerance, temperature
91
92! Note: omitting the following PEXSI options: deltaE (estimated by Arnoldi
93! before invoking PEXSI), mu0 (taken from previous SCF step), matrixType
94! (not implemented in PEXSI yet), isSymbolicFactorize (not needed because
95! of fixed sparsity pattern)
96
97 CALL section_vals_val_get(pexsi_section, "TEMPERATURE", &
98 r_val=temperature)
99 CALL section_vals_val_get(pexsi_section, "GAP", &
100 r_val=gap)
101 CALL section_vals_val_get(pexsi_section, "NUM_POLE", &
102 i_val=numpole)
103 CALL section_vals_val_get(pexsi_section, "IS_INERTIA_COUNT", &
104 l_val=isinertiacount)
105 CALL section_vals_val_get(pexsi_section, "MAX_PEXSI_ITER", &
106 i_val=maxpexsiiter)
107 CALL section_vals_val_get(pexsi_section, "MU_MIN_0", &
108 r_val=mumin0)
109 CALL section_vals_val_get(pexsi_section, "MU_MAX_0", &
110 r_val=mumax0)
111 CALL section_vals_val_get(pexsi_section, "MU_INERTIA_TOLERANCE", &
112 r_val=muinertiatolerance)
113 CALL section_vals_val_get(pexsi_section, "MU_INERTIA_EXPANSION", &
114 r_val=muinertiaexpansion)
115 CALL section_vals_val_get(pexsi_section, "MU_PEXSI_SAFE_GUARD", &
116 r_val=mupexsisafeguard)
117 CALL section_vals_val_get(pexsi_section, "NUM_ELECTRON_INITIAL_TOLERANCE", &
118 r_val=numelectroninitialtolerance)
119 CALL section_vals_val_get(pexsi_section, "NUM_ELECTRON_PEXSI_TOLERANCE", &
120 r_val=numelectrontargettolerance)
121 CALL section_vals_val_get(pexsi_section, "ORDERING", &
122 i_val=ordering)
123 CALL section_vals_val_get(pexsi_section, "ROW_ORDERING", &
124 i_val=rowordering)
125 CALL section_vals_val_get(pexsi_section, "NP_SYMB_FACT", &
126 i_val=npsymbfact)
127 CALL section_vals_val_get(pexsi_section, "VERBOSITY", &
128 i_val=verbosity)
129 CALL section_vals_val_get(pexsi_section, "MIN_RANKS_PER_POLE", &
130 i_val=min_ranks_per_pole)
131 CALL section_vals_val_get(pexsi_section, "CSR_SCREENING", &
132 l_val=csr_screening)
133
134 isinertiacount_int = merge(1, 0, isinertiacount) ! is integer in PEXSI
135
136 ! Set default options inside PEXSI
137 CALL cp_pexsi_set_default_options(pexsi_env%options)
138
139 ! Pass CP2K input to PEXSI options
140 CALL cp_pexsi_set_options(pexsi_env%options, temperature=temperature, gap=gap, &
141 numpole=numpole, isinertiacount=isinertiacount_int, maxpexsiiter=maxpexsiiter, &
142 mumin0=mumin0, mumax0=mumax0, muinertiatolerance=muinertiatolerance, &
143 muinertiaexpansion=muinertiaexpansion, mupexsisafeguard=mupexsisafeguard, &
144 ordering=ordering, rowordering=rowordering, npsymbfact=npsymbfact, verbosity=verbosity)
145
146 pexsi_env%num_ranks_per_pole = min_ranks_per_pole ! not a PEXSI option
147 pexsi_env%csr_screening = csr_screening
148
149 IF (numelectroninitialtolerance < numelectrontargettolerance) THEN
150 numelectroninitialtolerance = numelectrontargettolerance
151 END IF
152
153 pexsi_env%tol_nel_initial = numelectroninitialtolerance
154 pexsi_env%tol_nel_target = numelectrontargettolerance
155
156 END SUBROUTINE pexsi_init_read_input
157
158! **************************************************************************************************
159!> \brief Initializations needed before SCF
160!> \param ks_env ...
161!> \param pexsi_env ...
162!> \param template_matrix DBCSR matrix that defines the block structure and
163!> sparsity pattern of all matrices that are sent to PEXSI
164! **************************************************************************************************
165 SUBROUTINE pexsi_init_scf(ks_env, pexsi_env, template_matrix)
166 TYPE(qs_ks_env_type), POINTER :: ks_env
167 TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
168 TYPE(dbcsr_type), INTENT(IN) :: template_matrix
169
170 CHARACTER(len=*), PARAMETER :: routinen = 'pexsi_init_scf'
171
172 INTEGER :: handle, ispin, unit_nr
173 TYPE(cp_logger_type), POINTER :: logger
174
175 CALL timeset(routinen, handle)
176
177 logger => cp_get_default_logger()
178 IF (logger%para_env%is_source()) THEN
179 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
180 ELSE
181 unit_nr = -1
182 END IF
183
184 ! Create template matrices fixing sparsity pattern for PEXSI
185
186 IF (dbcsr_has_symmetry(template_matrix)) THEN
187 CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, template_matrix, &
188 "symmetric template matrix for CSR conversion")
189 CALL dbcsr_desymmetrize(pexsi_env%dbcsr_template_matrix_sym, &
190 pexsi_env%dbcsr_template_matrix_nonsym)
191 ELSE
192 CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_nonsym, template_matrix, &
193 "non-symmetric template matrix for CSR conversion")
194 CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, template_matrix, &
195 "symmetric template matrix for CSR conversion")
196 END IF
197
198 CALL dbcsr_create(pexsi_env%csr_sparsity, "CSR sparsity", &
199 template=pexsi_env%dbcsr_template_matrix_sym)
200 CALL dbcsr_copy(pexsi_env%csr_sparsity, pexsi_env%dbcsr_template_matrix_sym)
201
202 CALL cp_dbcsr_to_csr_screening(ks_env, pexsi_env%csr_sparsity)
203
204 IF (.NOT. pexsi_env%csr_screening) CALL dbcsr_set(pexsi_env%csr_sparsity, 1.0_dp)
205 CALL dbcsr_csr_create_from_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
206 pexsi_env%csr_mat_s, &
207 dbcsr_csr_eqrow_floor_dist, &
208 csr_sparsity=pexsi_env%csr_sparsity, &
209 numnodes=pexsi_env%num_ranks_per_pole)
210
211 IF (unit_nr > 0) WRITE (unit_nr, "(/T2,A)") "SPARSITY OF THE OVERLAP MATRIX IN CSR FORMAT"
212 CALL dbcsr_csr_print_sparsity(pexsi_env%csr_mat_s, unit_nr)
213
214 CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_nonsym, pexsi_env%csr_mat_s)
215
216 CALL dbcsr_csr_create(pexsi_env%csr_mat_ks, pexsi_env%csr_mat_s)
217 CALL dbcsr_csr_create(pexsi_env%csr_mat_p, pexsi_env%csr_mat_s)
218 CALL dbcsr_csr_create(pexsi_env%csr_mat_E, pexsi_env%csr_mat_s)
219 CALL dbcsr_csr_create(pexsi_env%csr_mat_F, pexsi_env%csr_mat_s)
220
221 DO ispin = 1, pexsi_env%nspin
222 ! max_ev_vector only initialised to avoid warning in case max_scf==0
223 CALL dbcsr_create(pexsi_env%matrix_w(ispin)%matrix, "W matrix", &
224 template=template_matrix, matrix_type=dbcsr_type_no_symmetry)
225 END DO
226
227 CALL cp_pexsi_set_options(pexsi_env%options, numelectronpexsitolerance=pexsi_env%tol_nel_initial)
228
229 CALL timestop(handle)
230
231 END SUBROUTINE pexsi_init_scf
232
233! **************************************************************************************************
234!> \brief Deallocations and post-processing after SCF
235!> \param pexsi_env ...
236!> \param mu_spin ...
237! **************************************************************************************************
238 SUBROUTINE pexsi_finalize_scf(pexsi_env, mu_spin)
239 TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
240 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: mu_spin
241
242 CHARACTER(len=*), PARAMETER :: routinen = 'pexsi_finalize_scf'
243
244 INTEGER :: handle, ispin, unit_nr
245 REAL(kind=dp) :: kts_total, mu_total
246 TYPE(cp_logger_type), POINTER :: logger
247
248 CALL timeset(routinen, handle)
249
250 logger => cp_get_default_logger()
251 IF (logger%para_env%is_source()) THEN
252 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
253 ELSE
254 unit_nr = -1
255 END IF
256
257 mu_total = sum(mu_spin(1:pexsi_env%nspin))/real(pexsi_env%nspin, kind=dp)
258 kts_total = sum(pexsi_env%kTS)
259 IF (pexsi_env%nspin == 1) kts_total = kts_total*2.0_dp
260
261 IF (unit_nr > 0) THEN
262 WRITE (unit_nr, "(/A,T55,F26.15)") &
263 " PEXSI| Electronic entropic energy (a.u.):", kts_total
264 WRITE (unit_nr, "(A,T55,F26.15)") &
265 " PEXSI| Chemical potential (a.u.):", mu_total
266 END IF
267
268 CALL dbcsr_release(pexsi_env%dbcsr_template_matrix_sym)
269 CALL dbcsr_release(pexsi_env%dbcsr_template_matrix_nonsym)
270 CALL dbcsr_release(pexsi_env%csr_sparsity)
271 CALL dbcsr_csr_destroy(pexsi_env%csr_mat_p)
272 CALL dbcsr_csr_destroy(pexsi_env%csr_mat_ks)
273 CALL dbcsr_csr_destroy(pexsi_env%csr_mat_s)
274 CALL dbcsr_csr_destroy(pexsi_env%csr_mat_E)
275 CALL dbcsr_csr_destroy(pexsi_env%csr_mat_F)
276 DO ispin = 1, pexsi_env%nspin
277 CALL dbcsr_release(pexsi_env%max_ev_vector(ispin))
278 CALL dbcsr_release(pexsi_env%matrix_w(ispin)%matrix)
279 END DO
280 CALL timestop(handle)
281 pexsi_env%tol_nel_initial = pexsi_env%tol_nel_target ! Turn off adaptive threshold for subsequent SCF cycles
282 END SUBROUTINE pexsi_finalize_scf
283
284! **************************************************************************************************
285!> \brief Calculate density matrix, energy-weighted density matrix and entropic
286!> energy contribution with the DFT driver of the PEXSI library.
287!> \param[in,out] pexsi_env PEXSI environment
288!> \param[in,out] matrix_p density matrix returned by PEXSI
289!> \param[in,out] matrix_w energy-weighted density matrix returned by PEXSI
290!> \param[out] kTS entropic energy contribution returned by PEXSI
291!> \param[in] matrix_ks Kohn-Sham matrix from linear scaling QS environment
292!> \param[in] matrix_s overlap matrix from linear scaling QS environment
293!> \param[in] nelectron_exact exact number of electrons
294!> \param[out] mu chemical potential calculated by PEXSI
295!> \param[in] iscf SCF step
296!> \param[in] ispin Number of spin
297!> \par History
298!> 11.2014 created [Patrick Seewald]
299!> \author Patrick Seewald
300! **************************************************************************************************
301 SUBROUTINE density_matrix_pexsi(pexsi_env, matrix_p, matrix_w, kTS, matrix_ks, matrix_s, &
302 nelectron_exact, mu, iscf, ispin)
303 TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
304 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
305 TYPE(dbcsr_p_type), INTENT(INOUT) :: matrix_w
306 REAL(kind=dp), INTENT(OUT) :: kts
307 TYPE(dbcsr_type), INTENT(IN), TARGET :: matrix_ks, matrix_s
308 INTEGER, INTENT(IN) :: nelectron_exact
309 REAL(kind=dp), INTENT(OUT) :: mu
310 INTEGER, INTENT(IN) :: iscf, ispin
311
312 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_pexsi'
313 INTEGER, PARAMETER :: s_not_identity = 0
314
315 INTEGER :: handle, is_symbolic_factorize, isinertiacount, isinertiacount_out, mynode, &
316 n_total_inertia_iter, n_total_pexsi_iter, unit_nr
317 LOGICAL :: first_call, pexsi_convergence
318 REAL(kind=dp) :: delta_e, energy_h, energy_s, free_energy, mu_max_in, mu_max_out, mu_min_in, &
319 mu_min_out, nel_tol, nelectron_diff, nelectron_exact_pexsi, nelectron_out
320 TYPE(arnoldi_env_type) :: arnoldi_env
321 TYPE(cp_logger_type), POINTER :: logger
322 TYPE(dbcsr_distribution_type) :: dist
323 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: arnoldi_matrices
324
325 CALL timeset(routinen, handle)
326
327 ! get a useful output_unit
328 logger => cp_get_default_logger()
329 IF (logger%para_env%is_source()) THEN
330 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
331 ELSE
332 unit_nr = -1
333 END IF
334
335 first_call = (iscf == 1) .AND. (ispin == 1)
336
337 ! Assert a few things the first time PEXSI is called
338 IF (first_call) THEN
339 ! Assertion that matrices have the expected symmetry (both should be symmetric if no
340 ! S preconditioning and no molecular clustering)
341 IF (.NOT. dbcsr_has_symmetry(matrix_ks)) THEN
342 cpabort("PEXSI interface expects a non-symmetric DBCSR Kohn-Sham matrix")
343 END IF
344 IF (.NOT. dbcsr_has_symmetry(matrix_s)) THEN
345 cpabort("PEXSI interface expects a non-symmetric DBCSR overlap matrix")
346 END IF
347
348 ! Assertion on datatype
349 IF ((pexsi_env%csr_mat_s%nzval_local%data_type /= dbcsr_csr_type_real_8) &
350 .OR. (pexsi_env%csr_mat_ks%nzval_local%data_type /= dbcsr_csr_type_real_8)) THEN
351 cpabort("Complex data type not supported by PEXSI")
352 END IF
353
354 ! Assertion on number of non-zero elements
355 !(TODO: update when PEXSI changes to Long Int)
356 IF (pexsi_env%csr_mat_s%nze_total >= int(2, kind=int_8)**31) THEN
357 cpabort("Total number of non-zero elements of CSR matrix is too large to be handled by PEXSI")
358 END IF
359 END IF
360
361 CALL dbcsr_get_info(matrix_ks, distribution=dist)
362 CALL dbcsr_distribution_get(dist, mynode=mynode)
363
364 ! Convert DBCSR matrices to PEXSI CSR format. Intermediate step to template matrix
365 ! needed in order to retain the initial sparsity pattern that is required for the
366 ! conversion to CSR format.
367 CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, matrix_s, keep_sparsity=.true.)
368 CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_sym, &
369 pexsi_env%csr_mat_s)
370
371 CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, matrix_ks, keep_sparsity=.true.)
372 CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_sym, &
373 pexsi_env%csr_mat_ks)
374
375 ! Get PEXSI input delta_E (upper bound for largest eigenvalue) using Arnoldi
376 NULLIFY (arnoldi_matrices)
377 CALL dbcsr_allocate_matrix_set(arnoldi_matrices, 2)
378 arnoldi_matrices(1)%matrix => matrix_ks
379 arnoldi_matrices(2)%matrix => matrix_s
380 CALL setup_arnoldi_env(arnoldi_env, arnoldi_matrices, max_iter=20, &
381 threshold=1.0e-2_dp, selection_crit=2, nval_request=1, nrestarts=21, &
382 generalized_ev=.true., iram=.false.)
383 IF (iscf > 1) CALL set_arnoldi_initial_vector(arnoldi_env, pexsi_env%max_ev_vector(ispin))
384 CALL arnoldi_ev(arnoldi_matrices, arnoldi_env)
385 delta_e = real(get_selected_ritz_val(arnoldi_env, 1), dp)
386 ! increase delta_E a bit to make sure that it really is an upper bound
387 delta_e = delta_e + 1.0e-2_dp*abs(delta_e)
388 CALL get_selected_ritz_vector(arnoldi_env, 1, arnoldi_matrices(1)%matrix, &
389 pexsi_env%max_ev_vector(ispin))
390 CALL deallocate_arnoldi_env(arnoldi_env)
391 DEALLOCATE (arnoldi_matrices)
392
393 nelectron_exact_pexsi = nelectron_exact
394
395 CALL cp_pexsi_set_options(pexsi_env%options, deltae=delta_e)
396
397 ! Set PEXSI options appropriately for first SCF iteration
398 IF (iscf == 1) THEN
399 ! Get option isInertiaCount to reset it later on and set it to 1 for first SCF iteration
400 CALL cp_pexsi_get_options(pexsi_env%options, isinertiacount=isinertiacount)
401 CALL cp_pexsi_set_options(pexsi_env%options, isinertiacount=1, &
402 issymbolicfactorize=1)
403 END IF
404
405 ! Write PEXSI options to output
406 CALL cp_pexsi_get_options(pexsi_env%options, isinertiacount=isinertiacount_out, &
407 issymbolicfactorize=is_symbolic_factorize, &
408 mumin0=mu_min_in, mumax0=mu_max_in, &
409 numelectronpexsitolerance=nel_tol)
410
411! IF(unit_nr>0) WRITE(unit_nr,'(/A,I4,A,I4)') " PEXSI| SCF", iscf, &
412! ", spin component", ispin
413
414 IF (unit_nr > 0) WRITE (unit_nr, '(/A,T41,L20)') " PEXSI| Do inertia counting?", &
415 isinertiacount_out == 1
416 IF (unit_nr > 0) WRITE (unit_nr, '(A,T50,F5.2,T56,F5.2)') &
417 " PEXSI| Guess for min mu, max mu", mu_min_in, mu_max_in
418
419 IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,E20.3)') &
420 " PEXSI| Tolerance in number of electrons", nel_tol
421
422! IF(unit_nr>0) WRITE(unit_nr,'(A,T41,L20)') &
423! " PEXSI| Do symbolic factorization?", is_symbolic_factorize==1
424
425 IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,F20.2)') &
426 " PEXSI| Arnoldi est. spectral radius", delta_e
427
428 ! Load data into PEXSI
429 CALL cp_pexsi_load_real_hs_matrix( &
430 pexsi_env%plan, &
431 pexsi_env%options, &
432 pexsi_env%csr_mat_ks%nrows_total, &
433 int(pexsi_env%csr_mat_ks%nze_total, kind=int_4), & ! TODO: update when PEXSI changes to Long Int
434 pexsi_env%csr_mat_ks%nze_local, &
435 pexsi_env%csr_mat_ks%nrows_local, &
436 pexsi_env%csr_mat_ks%rowptr_local, &
437 pexsi_env%csr_mat_ks%colind_local, &
438 pexsi_env%csr_mat_ks%nzval_local%r_dp, &
439 s_not_identity, &
440 pexsi_env%csr_mat_s%nzval_local%r_dp)
441
442 ! convert to spin restricted before passing number of electrons to PEXSI
443 CALL convert_nspin_cp2k_pexsi(cp2k_to_pexsi, &
444 numelectron=nelectron_exact_pexsi)
445
446 ! Call DFT driver of PEXSI doing the actual calculation
447 CALL cp_pexsi_dft_driver(pexsi_env%plan, pexsi_env%options, &
448 nelectron_exact_pexsi, mu, nelectron_out, mu_min_out, mu_max_out, &
449 n_total_inertia_iter, n_total_pexsi_iter)
450
451 ! Check convergence
452 nelectron_diff = nelectron_out - nelectron_exact_pexsi
453 pexsi_convergence = abs(nelectron_diff) < nel_tol
454
455 IF (unit_nr > 0) THEN
456 IF (pexsi_convergence) THEN
457 WRITE (unit_nr, '(/A)') " PEXSI| Converged"
458 ELSE
459 WRITE (unit_nr, '(/A)') " PEXSI| PEXSI did not converge!"
460 END IF
461
462! WRITE(unit_nr,'(A,T41,F20.10,T61,F20.10)') " PEXSI| Number of electrons", &
463! nelectron_out/pexsi_env%nspin, nelectron_diff/pexsi_env%nspin
464
465 WRITE (unit_nr, '(A,T41,F20.6)') " PEXSI| Chemical potential", mu
466
467 WRITE (unit_nr, '(A,T41,I20)') " PEXSI| PEXSI iterations", n_total_pexsi_iter
468 WRITE (unit_nr, '(A,T41,I20/)') " PEXSI| Inertia counting iterations", &
469 n_total_inertia_iter
470 END IF
471
472 IF (.NOT. pexsi_convergence) THEN
473 cpabort("PEXSI did not converge. Consider logPEXSI0 for more information.")
474 END IF
475
476 ! Retrieve results from PEXSI
477 IF (mynode < pexsi_env%mp_dims(1)*pexsi_env%mp_dims(2)) THEN
478 CALL cp_pexsi_retrieve_real_dft_matrix( &
479 pexsi_env%plan, &
480 pexsi_env%csr_mat_p%nzval_local%r_dp, &
481 pexsi_env%csr_mat_E%nzval_local%r_dp, &
482 pexsi_env%csr_mat_F%nzval_local%r_dp, &
483 energy_h, energy_s, free_energy)
484 ! calculate entropic energy contribution -TS = A - U
485 kts = (free_energy - energy_h)
486 END IF
487
488 ! send kTS to all nodes:
489 CALL pexsi_env%mp_group%bcast(kts, 0)
490
491 ! Convert PEXSI CSR matrices to DBCSR matrices
492 CALL dbcsr_convert_csr_to_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
493 pexsi_env%csr_mat_p)
494 CALL dbcsr_copy(matrix_p, pexsi_env%dbcsr_template_matrix_nonsym)
495 CALL dbcsr_convert_csr_to_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
496 pexsi_env%csr_mat_E)
497 CALL dbcsr_copy(matrix_w%matrix, pexsi_env%dbcsr_template_matrix_nonsym)
498
499 ! Convert to spin unrestricted
500 CALL convert_nspin_cp2k_pexsi(pexsi_to_cp2k, matrix_p=matrix_p, &
501 matrix_w=matrix_w, kts=kts)
502
503 ! Pass resulting mu as initial guess for next SCF to PEXSI
504 CALL cp_pexsi_set_options(pexsi_env%options, mu0=mu, mumin0=mu_min_out, &
505 mumax0=mu_max_out)
506
507 ! Reset isInertiaCount according to user input
508 IF (iscf == 1) THEN
509 CALL cp_pexsi_set_options(pexsi_env%options, isinertiacount= &
510 isinertiacount)
511 END IF
512
513 ! Turn off symbolic factorization for subsequent calls
514 IF (first_call) THEN
515 CALL cp_pexsi_set_options(pexsi_env%options, issymbolicfactorize=0)
516 END IF
517
518 CALL timestop(handle)
519 END SUBROUTINE density_matrix_pexsi
520
521! **************************************************************************************************
522!> \brief Set PEXSI convergence tolerance (numElectronPEXSITolerance), adapted
523!> to the convergence error of the previous SCF step. The tolerance is
524!> calculated using an initial convergence threshold (tol_nel_initial)
525!> and a target convergence threshold (tol_nel_target):
526!> numElectronPEXSITolerance(delta_scf) = alpha*delta_scf+beta
527!> where alpha and beta are chosen such that
528!> numElectronPEXSITolerance(delta_scf_0) = tol_nel_initial
529!> numElectronPEXSITolerance(eps_scf) = tol_nel_target
530!> and delta_scf_0 is the initial SCF convergence error.
531!> \param pexsi_env ...
532!> \param delta_scf convergence error of previous SCF step
533!> \param eps_scf SCF convergence criterion
534!> \param initialize whether or not adaptive thresholing should be initialized
535!> with delta_scf as initial convergence error
536!> \param check_convergence is set to .FALSE. if convergence in number of electrons
537!> will not be achieved in next SCF step
538! **************************************************************************************************
539 SUBROUTINE pexsi_set_convergence_tolerance(pexsi_env, delta_scf, eps_scf, initialize, &
540 check_convergence)
541 TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
542 REAL(kind=dp), INTENT(IN) :: delta_scf, eps_scf
543 LOGICAL, INTENT(IN) :: initialize
544 LOGICAL, INTENT(OUT) :: check_convergence
545
546 CHARACTER(len=*), PARAMETER :: routinen = 'pexsi_set_convergence_tolerance'
547
548 INTEGER :: handle
549 REAL(kind=dp) :: tol_nel
550
551 CALL timeset(routinen, handle)
552
553 tol_nel = pexsi_env%tol_nel_initial
554
555 IF (initialize) THEN
556 pexsi_env%adaptive_nel_alpha = &
557 (pexsi_env%tol_nel_initial - pexsi_env%tol_nel_target)/(abs(delta_scf) - eps_scf)
558 pexsi_env%adaptive_nel_beta = &
559 pexsi_env%tol_nel_initial - pexsi_env%adaptive_nel_alpha*abs(delta_scf)
560 pexsi_env%do_adaptive_tol_nel = .true.
561 END IF
562 IF (pexsi_env%do_adaptive_tol_nel) THEN
563 tol_nel = pexsi_env%adaptive_nel_alpha*abs(delta_scf) + pexsi_env%adaptive_nel_beta
564 tol_nel = max(tol_nel, pexsi_env%tol_nel_target)
565 tol_nel = min(tol_nel, pexsi_env%tol_nel_initial)
566 END IF
567
568 check_convergence = (tol_nel <= pexsi_env%tol_nel_target)
569
570 CALL cp_pexsi_set_options(pexsi_env%options, numelectronpexsitolerance=tol_nel)
571 CALL timestop(handle)
572 END SUBROUTINE pexsi_set_convergence_tolerance
573
574! **************************************************************************************************
575!> \brief Pass energy weighted density matrix and entropic energy contribution
576!> to QS environment
577!> \param ls_scf_env ...
578!> \param qs_env ...
579!> \param kTS ...
580!> \param matrix_w ...
581!> \par History
582!> 12.2014 created [Patrick Seewald]
583!> \author Patrick Seewald
584! **************************************************************************************************
585 SUBROUTINE pexsi_to_qs(ls_scf_env, qs_env, kTS, matrix_w)
586 TYPE(ls_scf_env_type) :: ls_scf_env
587 TYPE(qs_environment_type), INTENT(INOUT), POINTER :: qs_env
588 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: kts
589 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
590 OPTIONAL :: matrix_w
591
592 CHARACTER(len=*), PARAMETER :: routinen = 'pexsi_to_qs'
593
594 INTEGER :: handle, ispin, unit_nr
595 REAL(kind=dp) :: kts_total
596 TYPE(cp_logger_type), POINTER :: logger
597 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w_qs
598 TYPE(qs_energy_type), POINTER :: energy
599
600 CALL timeset(routinen, handle)
601
602 NULLIFY (energy)
603
604 ! get a useful output_unit
605 logger => cp_get_default_logger()
606 IF (logger%para_env%is_source()) THEN
607 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
608 ELSE
609 unit_nr = -1
610 END IF
611
612 CALL get_qs_env(qs_env, energy=energy, matrix_w=matrix_w_qs)
613
614 IF (PRESENT(matrix_w)) THEN
615 DO ispin = 1, ls_scf_env%nspins
616 CALL matrix_ls_to_qs(matrix_w_qs(ispin)%matrix, matrix_w(ispin)%matrix, &
617 ls_scf_env%ls_mstruct, covariant=.false.)
618 IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(matrix_w_qs(ispin)%matrix, 2.0_dp)
619 END DO
620 END IF
621
622 IF (PRESENT(kts)) THEN
623 kts_total = sum(kts)
624 IF (ls_scf_env%nspins == 1) kts_total = kts_total*2.0_dp
625 energy%kTS = kts_total
626 END IF
627
628 CALL timestop(handle)
629 END SUBROUTINE pexsi_to_qs
630
631END MODULE pexsi_methods
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:16
arnoldi_api::arnoldi_ev
subroutine, public arnoldi_ev(matrix, arnoldi_env)
Driver routine for different arnoldi eigenvalue methods the selection which one is to be taken is mad...
Definition arnoldi_api.F:69
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_has_symmetry
logical function, public dbcsr_has_symmetry(matrix)
...
Definition cp_dbcsr_api.F:928
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1188
cp_dbcsr_api::dbcsr_convert_dbcsr_to_csr
subroutine, public dbcsr_convert_dbcsr_to_csr(dbcsr_mat, csr_mat)
...
Definition cp_dbcsr_api.F:352
cp_dbcsr_api::dbcsr_convert_csr_to_dbcsr
subroutine, public dbcsr_convert_csr_to_dbcsr(dbcsr_mat, csr_mat)
...
Definition cp_dbcsr_api.F:336
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_csr_create_from_dbcsr
subroutine, public dbcsr_csr_create_from_dbcsr(dbcsr_mat, csr_mat, dist_format, csr_sparsity, numnodes)
...
Definition cp_dbcsr_api.F:469
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_distribution_get
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Definition cp_dbcsr_api.F:572
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_to_csr_screening
subroutine, public cp_dbcsr_to_csr_screening(ks_env, csr_sparsity)
Apply distance screening to refine sparsity pattern of matrices in CSR format (using eps_pgf_orb)....
Definition cp_dbcsr_cp2k_link.F:560
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:576
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
dm_ls_scf_qs
Routines for a linear scaling quickstep SCF run based on the density matrix, with a focus on the inte...
Definition dm_ls_scf_qs.F:15
dm_ls_scf_qs::matrix_ls_to_qs
subroutine, public matrix_ls_to_qs(matrix_qs, matrix_ls, ls_mstruct, covariant, keep_sparsity)
second link to QS, copy a LS matrix to QS matrix used to isolate QS style matrices from LS style will...
Definition dm_ls_scf_qs.F:308
dm_ls_scf_types
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Definition dm_ls_scf_types.F:15
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::int_4
integer, parameter, public int_4
Definition kinds.F:51
pexsi_interface
Interface to the PEXSI library, providing wrappers for all PEXSI routines that are called inside CP2K...
Definition pexsi_interface.F:15
pexsi_interface::cp_pexsi_set_options
subroutine, public cp_pexsi_set_options(pexsi_options, temperature, gap, deltae, numpole, isinertiacount, maxpexsiiter, mumin0, mumax0, mu0, muinertiatolerance, muinertiaexpansion, mupexsisafeguard, numelectronpexsitolerance, matrixtype, issymbolicfactorize, ordering, rowordering, npsymbfact, verbosity)
Set PEXSI internal options.
Definition pexsi_interface.F:88
pexsi_interface::cp_pexsi_retrieve_real_dft_matrix
subroutine, public cp_pexsi_retrieve_real_dft_matrix(plan, dmnzvallocal, edmnzvallocal, fdmnzvallocal, totalenergyh, totalenergys, totalfreeenergy)
...
Definition pexsi_interface.F:445
pexsi_interface::cp_pexsi_get_options
subroutine, public cp_pexsi_get_options(pexsi_options, temperature, gap, deltae, numpole, isinertiacount, maxpexsiiter, mumin0, mumax0, mu0, muinertiatolerance, muinertiaexpansion, mupexsisafeguard, numelectronpexsitolerance, matrixtype, issymbolicfactorize, ordering, rowordering, npsymbfact, verbosity)
Access PEXSI internal options.
Definition pexsi_interface.F:191
pexsi_interface::cp_pexsi_load_real_hs_matrix
subroutine, public cp_pexsi_load_real_hs_matrix(plan, pexsi_options, nrows, nnz, nnzlocal, numcollocal, colptrlocal, rowindlocal, hnzvallocal, issidentity, snzvallocal)
...
Definition pexsi_interface.F:326
pexsi_interface::cp_pexsi_set_default_options
subroutine, public cp_pexsi_set_default_options(pexsi_options)
...
Definition pexsi_interface.F:264
pexsi_interface::cp_pexsi_dft_driver
subroutine, public cp_pexsi_dft_driver(plan, pexsi_options, numelectronexact, mupexsi, numelectronpexsi, mumininertia, mumaxinertia, numtotalinertiaiter, numtotalpexsiiter)
...
Definition pexsi_interface.F:381
pexsi_methods
Methods using the PEXSI library to calculate the density matrix and related quantities using the Kohn...
Definition pexsi_methods.F:17
pexsi_methods::pexsi_to_qs
subroutine, public pexsi_to_qs(ls_scf_env, qs_env, kts, matrix_w)
Pass energy weighted density matrix and entropic energy contribution to QS environment.
Definition pexsi_methods.F:586
pexsi_methods::density_matrix_pexsi
subroutine, public density_matrix_pexsi(pexsi_env, matrix_p, matrix_w, kts, matrix_ks, matrix_s, nelectron_exact, mu, iscf, ispin)
Calculate density matrix, energy-weighted density matrix and entropic energy contribution with the DF...
Definition pexsi_methods.F:303
pexsi_methods::pexsi_init_read_input
subroutine, public pexsi_init_read_input(pexsi_section, pexsi_env)
Read CP2K input section PEXSI and pass it to the PEXSI environment.
Definition pexsi_methods.F:81
pexsi_methods::pexsi_finalize_scf
subroutine, public pexsi_finalize_scf(pexsi_env, mu_spin)
Deallocations and post-processing after SCF.
Definition pexsi_methods.F:239
pexsi_methods::pexsi_init_scf
subroutine, public pexsi_init_scf(ks_env, pexsi_env, template_matrix)
Initializations needed before SCF.
Definition pexsi_methods.F:166
pexsi_methods::pexsi_set_convergence_tolerance
subroutine, public pexsi_set_convergence_tolerance(pexsi_env, delta_scf, eps_scf, initialize, check_convergence)
Set PEXSI convergence tolerance (numElectronPEXSITolerance), adapted to the convergence error of the ...
Definition pexsi_methods.F:541
pexsi_types
Environment storing all data that is needed in order to call the DFT driver of the PEXSI library with...
Definition pexsi_types.F:18
pexsi_types::convert_nspin_cp2k_pexsi
subroutine, public convert_nspin_cp2k_pexsi(direction, numelectron, matrix_p, matrix_w, kts)
Scale various quantities with factors of 2. This converts spin restricted DFT calculations (PEXSI) to...
Definition pexsi_types.F:238
pexsi_types::cp2k_to_pexsi
integer, parameter, public cp2k_to_pexsi
Definition pexsi_types.F:48
pexsi_types::pexsi_to_cp2k
integer, parameter, public pexsi_to_cp2k
Definition pexsi_types.F:48
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_ks_types
Definition qs_ks_types.F:15
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
dm_ls_scf_types::ls_scf_env_type
Definition dm_ls_scf_types.F:88
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
pexsi_types::lib_pexsi_env
All PEXSI related data.
Definition pexsi_types.F:87
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137