(git:374b731)
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-2024 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_data_type,&
19 arnoldi_ev,&
20 deallocate_arnoldi_data,&
21 get_selected_ritz_val,&
22 get_selected_ritz_vector,&
23 set_arnoldi_initial_vector,&
24 setup_arnoldi_data
25 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_to_csr_screening
26 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
27 USE cp_log_handling, ONLY: cp_get_default_logger,&
28 cp_logger_get_default_unit_nr,&
29 cp_logger_type
30 USE dbcsr_api, ONLY: &
31 dbcsr_convert_csr_to_dbcsr, dbcsr_convert_dbcsr_to_csr, dbcsr_copy, &
32 dbcsr_copy_into_existing, dbcsr_create, dbcsr_csr_create, dbcsr_csr_create_from_dbcsr, &
33 dbcsr_csr_destroy, dbcsr_csr_eqrow_floor_dist, dbcsr_csr_print_sparsity, &
34 dbcsr_desymmetrize, dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_info, &
35 dbcsr_has_symmetry, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, &
36 dbcsr_type_no_symmetry, dbcsr_type_real_8
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 .LT. numelectrontargettolerance) &
150 numelectroninitialtolerance = numelectrontargettolerance
151
152 pexsi_env%tol_nel_initial = numelectroninitialtolerance
153 pexsi_env%tol_nel_target = numelectrontargettolerance
154
155 END SUBROUTINE pexsi_init_read_input
156
157! **************************************************************************************************
158!> \brief Initializations needed before SCF
159!> \param ks_env ...
160!> \param pexsi_env ...
161!> \param template_matrix DBCSR matrix that defines the block structure and
162!> sparsity pattern of all matrices that are sent to PEXSI
163! **************************************************************************************************
164 SUBROUTINE pexsi_init_scf(ks_env, pexsi_env, template_matrix)
165 TYPE(qs_ks_env_type), POINTER :: ks_env
166 TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
167 TYPE(dbcsr_type), INTENT(IN) :: template_matrix
168
169 CHARACTER(len=*), PARAMETER :: routinen = 'pexsi_init_scf'
170
171 INTEGER :: handle, ispin, unit_nr
172 TYPE(cp_logger_type), POINTER :: logger
173
174 CALL timeset(routinen, handle)
175
176 logger => cp_get_default_logger()
177 IF (logger%para_env%is_source()) THEN
178 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
179 ELSE
180 unit_nr = -1
181 END IF
182
183 ! Create template matrices fixing sparsity pattern for PEXSI
184
185 IF (dbcsr_has_symmetry(template_matrix)) THEN
186 CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, template_matrix, &
187 "symmetric template matrix for CSR conversion")
188 CALL dbcsr_desymmetrize(pexsi_env%dbcsr_template_matrix_sym, &
189 pexsi_env%dbcsr_template_matrix_nonsym)
190 ELSE
191 CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_nonsym, template_matrix, &
192 "non-symmetric template matrix for CSR conversion")
193 CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, template_matrix, &
194 "symmetric template matrix for CSR conversion")
195 END IF
196
197 CALL dbcsr_create(pexsi_env%csr_sparsity, "CSR sparsity", &
198 template=pexsi_env%dbcsr_template_matrix_sym, &
199 data_type=dbcsr_type_real_8)
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)
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 .EQ. 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_data_type) :: my_arnoldi
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 .EQ. 1) .AND. (ispin .EQ. 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)) &
342 cpabort("PEXSI interface expects a non-symmetric DBCSR Kohn-Sham matrix")
343 IF (.NOT. dbcsr_has_symmetry(matrix_s)) &
344 cpabort("PEXSI interface expects a non-symmetric DBCSR overlap matrix")
345
346 ! Assertion on datatype
347 IF ((pexsi_env%csr_mat_s%nzval_local%data_type .NE. dbcsr_type_real_8) &
348 .OR. (pexsi_env%csr_mat_ks%nzval_local%data_type .NE. dbcsr_type_real_8)) &
349 cpabort("Complex data type not supported by PEXSI")
350
351 ! Assertion on number of non-zero elements
352 !(TODO: update when PEXSI changes to Long Int)
353 IF (pexsi_env%csr_mat_s%nze_total .GE. int(2, kind=int_8)**31) &
354 cpabort("Total number of non-zero elements of CSR matrix is too large to be handled by PEXSI")
355 END IF
356
357 CALL dbcsr_get_info(matrix_ks, distribution=dist)
358 CALL dbcsr_distribution_get(dist, mynode=mynode)
359
360 ! Convert DBCSR matrices to PEXSI CSR format. Intermediate step to template matrix
361 ! needed in order to retain the initial sparsity pattern that is required for the
362 ! conversion to CSR format.
363 CALL dbcsr_copy_into_existing(pexsi_env%dbcsr_template_matrix_sym, matrix_s)
364 CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_sym, &
365 pexsi_env%csr_mat_s)
366
367 CALL dbcsr_copy_into_existing(pexsi_env%dbcsr_template_matrix_sym, &
368 matrix_ks)
369 CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_sym, &
370 pexsi_env%csr_mat_ks)
371
372 ! Get PEXSI input delta_E (upper bound for largest eigenvalue) using Arnoldi
373 NULLIFY (arnoldi_matrices)
374 CALL dbcsr_allocate_matrix_set(arnoldi_matrices, 2)
375 arnoldi_matrices(1)%matrix => matrix_ks
376 arnoldi_matrices(2)%matrix => matrix_s
377 CALL setup_arnoldi_data(my_arnoldi, arnoldi_matrices, max_iter=20, &
378 threshold=1.0e-2_dp, selection_crit=2, nval_request=1, nrestarts=21, &
379 generalized_ev=.true., iram=.false.)
380 IF (iscf .GT. 1) CALL set_arnoldi_initial_vector(my_arnoldi, pexsi_env%max_ev_vector(ispin))
381 CALL arnoldi_ev(arnoldi_matrices, my_arnoldi)
382 delta_e = real(get_selected_ritz_val(my_arnoldi, 1), dp)
383 ! increase delta_E a bit to make sure that it really is an upper bound
384 delta_e = delta_e + 1.0e-2_dp*abs(delta_e)
385 CALL get_selected_ritz_vector(my_arnoldi, 1, arnoldi_matrices(1)%matrix, &
386 pexsi_env%max_ev_vector(ispin))
387 CALL deallocate_arnoldi_data(my_arnoldi)
388 DEALLOCATE (arnoldi_matrices)
389
390 nelectron_exact_pexsi = nelectron_exact
391
392 CALL cp_pexsi_set_options(pexsi_env%options, deltae=delta_e)
393
394 ! Set PEXSI options appropriately for first SCF iteration
395 IF (iscf .EQ. 1) THEN
396 ! Get option isInertiaCount to reset it later on and set it to 1 for first SCF iteration
397 CALL cp_pexsi_get_options(pexsi_env%options, isinertiacount=isinertiacount)
398 CALL cp_pexsi_set_options(pexsi_env%options, isinertiacount=1, &
399 issymbolicfactorize=1)
400 END IF
401
402 ! Write PEXSI options to output
403 CALL cp_pexsi_get_options(pexsi_env%options, isinertiacount=isinertiacount_out, &
404 issymbolicfactorize=is_symbolic_factorize, &
405 mumin0=mu_min_in, mumax0=mu_max_in, &
406 numelectronpexsitolerance=nel_tol)
407
408! IF(unit_nr>0) WRITE(unit_nr,'(/A,I4,A,I4)') " PEXSI| SCF", iscf, &
409! ", spin component", ispin
410
411 IF (unit_nr > 0) WRITE (unit_nr, '(/A,T41,L20)') " PEXSI| Do inertia counting?", &
412 isinertiacount_out .EQ. 1
413 IF (unit_nr > 0) WRITE (unit_nr, '(A,T50,F5.2,T56,F5.2)') &
414 " PEXSI| Guess for min mu, max mu", mu_min_in, mu_max_in
415
416 IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,E20.3)') &
417 " PEXSI| Tolerance in number of electrons", nel_tol
418
419! IF(unit_nr>0) WRITE(unit_nr,'(A,T41,L20)') &
420! " PEXSI| Do symbolic factorization?", is_symbolic_factorize.EQ.1
421
422 IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,F20.2)') &
423 " PEXSI| Arnoldi est. spectral radius", delta_e
424
425 ! Load data into PEXSI
426 CALL cp_pexsi_load_real_hs_matrix( &
427 pexsi_env%plan, &
428 pexsi_env%options, &
429 pexsi_env%csr_mat_ks%nrows_total, &
430 int(pexsi_env%csr_mat_ks%nze_total, kind=int_4), & ! TODO: update when PEXSI changes to Long Int
431 pexsi_env%csr_mat_ks%nze_local, &
432 pexsi_env%csr_mat_ks%nrows_local, &
433 pexsi_env%csr_mat_ks%rowptr_local, &
434 pexsi_env%csr_mat_ks%colind_local, &
435 pexsi_env%csr_mat_ks%nzval_local%r_dp, &
436 s_not_identity, &
437 pexsi_env%csr_mat_s%nzval_local%r_dp)
438
439 ! convert to spin restricted before passing number of electrons to PEXSI
440 CALL convert_nspin_cp2k_pexsi(cp2k_to_pexsi, &
441 numelectron=nelectron_exact_pexsi)
442
443 ! Call DFT driver of PEXSI doing the actual calculation
444 CALL cp_pexsi_dft_driver(pexsi_env%plan, pexsi_env%options, &
445 nelectron_exact_pexsi, mu, nelectron_out, mu_min_out, mu_max_out, &
446 n_total_inertia_iter, n_total_pexsi_iter)
447
448 ! Check convergence
449 nelectron_diff = nelectron_out - nelectron_exact_pexsi
450 pexsi_convergence = abs(nelectron_diff) .LT. nel_tol
451
452 IF (unit_nr > 0) THEN
453 IF (pexsi_convergence) THEN
454 WRITE (unit_nr, '(/A)') " PEXSI| Converged"
455 ELSE
456 WRITE (unit_nr, '(/A)') " PEXSI| PEXSI did not converge!"
457 END IF
458
459! WRITE(unit_nr,'(A,T41,F20.10,T61,F20.10)') " PEXSI| Number of electrons", &
460! nelectron_out/pexsi_env%nspin, nelectron_diff/pexsi_env%nspin
461
462 WRITE (unit_nr, '(A,T41,F20.6)') " PEXSI| Chemical potential", mu
463
464 WRITE (unit_nr, '(A,T41,I20)') " PEXSI| PEXSI iterations", n_total_pexsi_iter
465 WRITE (unit_nr, '(A,T41,I20/)') " PEXSI| Inertia counting iterations", &
466 n_total_inertia_iter
467 END IF
468
469 IF (.NOT. pexsi_convergence) &
470 cpabort("PEXSI did not converge. Consider logPEXSI0 for more information.")
471
472 ! Retrieve results from PEXSI
473 IF (mynode < pexsi_env%mp_dims(1)*pexsi_env%mp_dims(2)) THEN
474 CALL cp_pexsi_retrieve_real_dft_matrix( &
475 pexsi_env%plan, &
476 pexsi_env%csr_mat_p%nzval_local%r_dp, &
477 pexsi_env%csr_mat_E%nzval_local%r_dp, &
478 pexsi_env%csr_mat_F%nzval_local%r_dp, &
479 energy_h, energy_s, free_energy)
480 ! calculate entropic energy contribution -TS = A - U
481 kts = (free_energy - energy_h)
482 END IF
483
484 ! send kTS to all nodes:
485 CALL pexsi_env%mp_group%bcast(kts, 0)
486
487 ! Convert PEXSI CSR matrices to DBCSR matrices
488 CALL dbcsr_convert_csr_to_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
489 pexsi_env%csr_mat_p)
490 CALL dbcsr_copy(matrix_p, pexsi_env%dbcsr_template_matrix_nonsym)
491 CALL dbcsr_convert_csr_to_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
492 pexsi_env%csr_mat_E)
493 CALL dbcsr_copy(matrix_w%matrix, pexsi_env%dbcsr_template_matrix_nonsym)
494
495 ! Convert to spin unrestricted
496 CALL convert_nspin_cp2k_pexsi(pexsi_to_cp2k, matrix_p=matrix_p, &
497 matrix_w=matrix_w, kts=kts)
498
499 ! Pass resulting mu as initial guess for next SCF to PEXSI
500 CALL cp_pexsi_set_options(pexsi_env%options, mu0=mu, mumin0=mu_min_out, &
501 mumax0=mu_max_out)
502
503 ! Reset isInertiaCount according to user input
504 IF (iscf .EQ. 1) THEN
505 CALL cp_pexsi_set_options(pexsi_env%options, isinertiacount= &
506 isinertiacount)
507 END IF
508
509 ! Turn off symbolic factorization for subsequent calls
510 IF (first_call) THEN
511 CALL cp_pexsi_set_options(pexsi_env%options, issymbolicfactorize=0)
512 END IF
513
514 CALL timestop(handle)
515 END SUBROUTINE density_matrix_pexsi
516
517! **************************************************************************************************
518!> \brief Set PEXSI convergence tolerance (numElectronPEXSITolerance), adapted
519!> to the convergence error of the previous SCF step. The tolerance is
520!> calculated using an initial convergence threshold (tol_nel_initial)
521!> and a target convergence threshold (tol_nel_target):
522!> numElectronPEXSITolerance(delta_scf) = alpha*delta_scf+beta
523!> where alpha and beta are chosen such that
524!> numElectronPEXSITolerance(delta_scf_0) = tol_nel_initial
525!> numElectronPEXSITolerance(eps_scf) = tol_nel_target
526!> and delta_scf_0 is the initial SCF convergence error.
527!> \param pexsi_env ...
528!> \param delta_scf convergence error of previous SCF step
529!> \param eps_scf SCF convergence criterion
530!> \param initialize whether or not adaptive thresholing should be initialized
531!> with delta_scf as initial convergence error
532!> \param check_convergence is set to .FALSE. if convergence in number of electrons
533!> will not be achieved in next SCF step
534! **************************************************************************************************
535 SUBROUTINE pexsi_set_convergence_tolerance(pexsi_env, delta_scf, eps_scf, initialize, &
536 check_convergence)
537 TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
538 REAL(kind=dp), INTENT(IN) :: delta_scf, eps_scf
539 LOGICAL, INTENT(IN) :: initialize
540 LOGICAL, INTENT(OUT) :: check_convergence
541
542 CHARACTER(len=*), PARAMETER :: routinen = 'pexsi_set_convergence_tolerance'
543
544 INTEGER :: handle
545 REAL(kind=dp) :: tol_nel
546
547 CALL timeset(routinen, handle)
548
549 tol_nel = pexsi_env%tol_nel_initial
550
551 IF (initialize) THEN
552 pexsi_env%adaptive_nel_alpha = &
553 (pexsi_env%tol_nel_initial - pexsi_env%tol_nel_target)/(abs(delta_scf) - eps_scf)
554 pexsi_env%adaptive_nel_beta = &
555 pexsi_env%tol_nel_initial - pexsi_env%adaptive_nel_alpha*abs(delta_scf)
556 pexsi_env%do_adaptive_tol_nel = .true.
557 END IF
558 IF (pexsi_env%do_adaptive_tol_nel) THEN
559 tol_nel = pexsi_env%adaptive_nel_alpha*abs(delta_scf) + pexsi_env%adaptive_nel_beta
560 tol_nel = max(tol_nel, pexsi_env%tol_nel_target)
561 tol_nel = min(tol_nel, pexsi_env%tol_nel_initial)
562 END IF
563
564 check_convergence = (tol_nel .LE. pexsi_env%tol_nel_target)
565
566 CALL cp_pexsi_set_options(pexsi_env%options, numelectronpexsitolerance=tol_nel)
567 CALL timestop(handle)
568 END SUBROUTINE
569
570! **************************************************************************************************
571!> \brief Pass energy weighted density matrix and entropic energy contribution
572!> to QS environment
573!> \param ls_scf_env ...
574!> \param qs_env ...
575!> \param kTS ...
576!> \param matrix_w ...
577!> \par History
578!> 12.2014 created [Patrick Seewald]
579!> \author Patrick Seewald
580! **************************************************************************************************
581 SUBROUTINE pexsi_to_qs(ls_scf_env, qs_env, kTS, matrix_w)
582 TYPE(ls_scf_env_type) :: ls_scf_env
583 TYPE(qs_environment_type), INTENT(INOUT), POINTER :: qs_env
584 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: kts
585 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
586 OPTIONAL :: matrix_w
587
588 CHARACTER(len=*), PARAMETER :: routinen = 'pexsi_to_qs'
589
590 INTEGER :: handle, ispin, unit_nr
591 REAL(kind=dp) :: kts_total
592 TYPE(cp_logger_type), POINTER :: logger
593 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w_qs
594 TYPE(qs_energy_type), POINTER :: energy
595
596 CALL timeset(routinen, handle)
597
598 NULLIFY (energy)
599
600 ! get a useful output_unit
601 logger => cp_get_default_logger()
602 IF (logger%para_env%is_source()) THEN
603 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
604 ELSE
605 unit_nr = -1
606 END IF
607
608 CALL get_qs_env(qs_env, energy=energy, matrix_w=matrix_w_qs)
609
610 IF (PRESENT(matrix_w)) THEN
611 DO ispin = 1, ls_scf_env%nspins
612 CALL matrix_ls_to_qs(matrix_w_qs(ispin)%matrix, matrix_w(ispin)%matrix, &
613 ls_scf_env%ls_mstruct, covariant=.false.)
614 IF (ls_scf_env%nspins .EQ. 1) CALL dbcsr_scale(matrix_w_qs(ispin)%matrix, 2.0_dp)
615 END DO
616 END IF
617
618 IF (PRESENT(kts)) THEN
619 kts_total = sum(kts)
620 IF (ls_scf_env%nspins .EQ. 1) kts_total = kts_total*2.0_dp
621 energy%kTS = kts_total
622 END IF
623
624 CALL timestop(handle)
625 END SUBROUTINE pexsi_to_qs
626
627END MODULE pexsi_methods
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:81
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:15
arnoldi_api::arnoldi_ev
subroutine, public arnoldi_ev(matrix, arnoldi_data)
Driver routine for different arnoldi eigenvalue methods the selection which one is to be taken is mad...
Definition arnoldi_api.F:69
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:567
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:307
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:1040
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:430
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:186
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:315
pexsi_interface::cp_pexsi_set_default_options
subroutine, public cp_pexsi_set_default_options(pexsi_options)
...
Definition pexsi_interface.F:254
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:367
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:582
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:165
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:537
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:233
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_ks_types
Definition qs_ks_types.F:15
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:126
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:215
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:135