(git:ed6f26b)
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-2025 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 .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 CALL dbcsr_copy(pexsi_env%csr_sparsity, pexsi_env%dbcsr_template_matrix_sym)
200
201 CALL cp_dbcsr_to_csr_screening(ks_env, pexsi_env%csr_sparsity)
202
203 IF (.NOT. pexsi_env%csr_screening) CALL dbcsr_set(pexsi_env%csr_sparsity, 1.0_dp)
204 CALL dbcsr_csr_create_from_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
205 pexsi_env%csr_mat_s, &
206 dbcsr_csr_eqrow_floor_dist, &
207 csr_sparsity=pexsi_env%csr_sparsity, &
208 numnodes=pexsi_env%num_ranks_per_pole)
209
210 IF (unit_nr > 0) WRITE (unit_nr, "(/T2,A)") "SPARSITY OF THE OVERLAP MATRIX IN CSR FORMAT"
211 CALL dbcsr_csr_print_sparsity(pexsi_env%csr_mat_s, unit_nr)
212
213 CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_nonsym, pexsi_env%csr_mat_s)
214
215 CALL dbcsr_csr_create(pexsi_env%csr_mat_ks, pexsi_env%csr_mat_s)
216 CALL dbcsr_csr_create(pexsi_env%csr_mat_p, pexsi_env%csr_mat_s)
217 CALL dbcsr_csr_create(pexsi_env%csr_mat_E, pexsi_env%csr_mat_s)
218 CALL dbcsr_csr_create(pexsi_env%csr_mat_F, pexsi_env%csr_mat_s)
219
220 DO ispin = 1, pexsi_env%nspin
221 ! max_ev_vector only initialised to avoid warning in case max_scf==0
222 CALL dbcsr_create(pexsi_env%matrix_w(ispin)%matrix, "W matrix", &
223 template=template_matrix, matrix_type=dbcsr_type_no_symmetry)
224 END DO
225
226 CALL cp_pexsi_set_options(pexsi_env%options, numelectronpexsitolerance=pexsi_env%tol_nel_initial)
227
228 CALL timestop(handle)
229
230 END SUBROUTINE pexsi_init_scf
231
232! **************************************************************************************************
233!> \brief Deallocations and post-processing after SCF
234!> \param pexsi_env ...
235!> \param mu_spin ...
236! **************************************************************************************************
237 SUBROUTINE pexsi_finalize_scf(pexsi_env, mu_spin)
238 TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
239 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: mu_spin
240
241 CHARACTER(len=*), PARAMETER :: routinen = 'pexsi_finalize_scf'
242
243 INTEGER :: handle, ispin, unit_nr
244 REAL(kind=dp) :: kts_total, mu_total
245 TYPE(cp_logger_type), POINTER :: logger
246
247 CALL timeset(routinen, handle)
248
249 logger => cp_get_default_logger()
250 IF (logger%para_env%is_source()) THEN
251 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
252 ELSE
253 unit_nr = -1
254 END IF
255
256 mu_total = sum(mu_spin(1:pexsi_env%nspin))/real(pexsi_env%nspin, kind=dp)
257 kts_total = sum(pexsi_env%kTS)
258 IF (pexsi_env%nspin .EQ. 1) kts_total = kts_total*2.0_dp
259
260 IF (unit_nr > 0) THEN
261 WRITE (unit_nr, "(/A,T55,F26.15)") &
262 " PEXSI| Electronic entropic energy (a.u.):", kts_total
263 WRITE (unit_nr, "(A,T55,F26.15)") &
264 " PEXSI| Chemical potential (a.u.):", mu_total
265 END IF
266
267 CALL dbcsr_release(pexsi_env%dbcsr_template_matrix_sym)
268 CALL dbcsr_release(pexsi_env%dbcsr_template_matrix_nonsym)
269 CALL dbcsr_release(pexsi_env%csr_sparsity)
270 CALL dbcsr_csr_destroy(pexsi_env%csr_mat_p)
271 CALL dbcsr_csr_destroy(pexsi_env%csr_mat_ks)
272 CALL dbcsr_csr_destroy(pexsi_env%csr_mat_s)
273 CALL dbcsr_csr_destroy(pexsi_env%csr_mat_E)
274 CALL dbcsr_csr_destroy(pexsi_env%csr_mat_F)
275 DO ispin = 1, pexsi_env%nspin
276 CALL dbcsr_release(pexsi_env%max_ev_vector(ispin))
277 CALL dbcsr_release(pexsi_env%matrix_w(ispin)%matrix)
278 END DO
279 CALL timestop(handle)
280 pexsi_env%tol_nel_initial = pexsi_env%tol_nel_target ! Turn off adaptive threshold for subsequent SCF cycles
281 END SUBROUTINE pexsi_finalize_scf
282
283! **************************************************************************************************
284!> \brief Calculate density matrix, energy-weighted density matrix and entropic
285!> energy contribution with the DFT driver of the PEXSI library.
286!> \param[in,out] pexsi_env PEXSI environment
287!> \param[in,out] matrix_p density matrix returned by PEXSI
288!> \param[in,out] matrix_w energy-weighted density matrix returned by PEXSI
289!> \param[out] kTS entropic energy contribution returned by PEXSI
290!> \param[in] matrix_ks Kohn-Sham matrix from linear scaling QS environment
291!> \param[in] matrix_s overlap matrix from linear scaling QS environment
292!> \param[in] nelectron_exact exact number of electrons
293!> \param[out] mu chemical potential calculated by PEXSI
294!> \param[in] iscf SCF step
295!> \param[in] ispin Number of spin
296!> \par History
297!> 11.2014 created [Patrick Seewald]
298!> \author Patrick Seewald
299! **************************************************************************************************
300 SUBROUTINE density_matrix_pexsi(pexsi_env, matrix_p, matrix_w, kTS, matrix_ks, matrix_s, &
301 nelectron_exact, mu, iscf, ispin)
302 TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
303 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
304 TYPE(dbcsr_p_type), INTENT(INOUT) :: matrix_w
305 REAL(kind=dp), INTENT(OUT) :: kts
306 TYPE(dbcsr_type), INTENT(IN), TARGET :: matrix_ks, matrix_s
307 INTEGER, INTENT(IN) :: nelectron_exact
308 REAL(kind=dp), INTENT(OUT) :: mu
309 INTEGER, INTENT(IN) :: iscf, ispin
310
311 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_pexsi'
312 INTEGER, PARAMETER :: s_not_identity = 0
313
314 INTEGER :: handle, is_symbolic_factorize, isinertiacount, isinertiacount_out, mynode, &
315 n_total_inertia_iter, n_total_pexsi_iter, unit_nr
316 LOGICAL :: first_call, pexsi_convergence
317 REAL(kind=dp) :: delta_e, energy_h, energy_s, free_energy, mu_max_in, mu_max_out, mu_min_in, &
318 mu_min_out, nel_tol, nelectron_diff, nelectron_exact_pexsi, nelectron_out
319 TYPE(arnoldi_env_type) :: arnoldi_env
320 TYPE(cp_logger_type), POINTER :: logger
321 TYPE(dbcsr_distribution_type) :: dist
322 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: arnoldi_matrices
323
324 CALL timeset(routinen, handle)
325
326 ! get a useful output_unit
327 logger => cp_get_default_logger()
328 IF (logger%para_env%is_source()) THEN
329 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
330 ELSE
331 unit_nr = -1
332 END IF
333
334 first_call = (iscf .EQ. 1) .AND. (ispin .EQ. 1)
335
336 ! Assert a few things the first time PEXSI is called
337 IF (first_call) THEN
338 ! Assertion that matrices have the expected symmetry (both should be symmetric if no
339 ! S preconditioning and no molecular clustering)
340 IF (.NOT. dbcsr_has_symmetry(matrix_ks)) &
341 cpabort("PEXSI interface expects a non-symmetric DBCSR Kohn-Sham matrix")
342 IF (.NOT. dbcsr_has_symmetry(matrix_s)) &
343 cpabort("PEXSI interface expects a non-symmetric DBCSR overlap matrix")
344
345 ! Assertion on datatype
346 IF ((pexsi_env%csr_mat_s%nzval_local%data_type .NE. dbcsr_csr_type_real_8) &
347 .OR. (pexsi_env%csr_mat_ks%nzval_local%data_type .NE. dbcsr_csr_type_real_8)) &
348 cpabort("Complex data type not supported by PEXSI")
349
350 ! Assertion on number of non-zero elements
351 !(TODO: update when PEXSI changes to Long Int)
352 IF (pexsi_env%csr_mat_s%nze_total .GE. int(2, kind=int_8)**31) &
353 cpabort("Total number of non-zero elements of CSR matrix is too large to be handled by PEXSI")
354 END IF
355
356 CALL dbcsr_get_info(matrix_ks, distribution=dist)
357 CALL dbcsr_distribution_get(dist, mynode=mynode)
358
359 ! Convert DBCSR matrices to PEXSI CSR format. Intermediate step to template matrix
360 ! needed in order to retain the initial sparsity pattern that is required for the
361 ! conversion to CSR format.
362 CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, matrix_s, keep_sparsity=.true.)
363 CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_sym, &
364 pexsi_env%csr_mat_s)
365
366 CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, matrix_ks, keep_sparsity=.true.)
367 CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_sym, &
368 pexsi_env%csr_mat_ks)
369
370 ! Get PEXSI input delta_E (upper bound for largest eigenvalue) using Arnoldi
371 NULLIFY (arnoldi_matrices)
372 CALL dbcsr_allocate_matrix_set(arnoldi_matrices, 2)
373 arnoldi_matrices(1)%matrix => matrix_ks
374 arnoldi_matrices(2)%matrix => matrix_s
375 CALL setup_arnoldi_env(arnoldi_env, arnoldi_matrices, max_iter=20, &
376 threshold=1.0e-2_dp, selection_crit=2, nval_request=1, nrestarts=21, &
377 generalized_ev=.true., iram=.false.)
378 IF (iscf .GT. 1) CALL set_arnoldi_initial_vector(arnoldi_env, pexsi_env%max_ev_vector(ispin))
379 CALL arnoldi_ev(arnoldi_matrices, arnoldi_env)
380 delta_e = real(get_selected_ritz_val(arnoldi_env, 1), dp)
381 ! increase delta_E a bit to make sure that it really is an upper bound
382 delta_e = delta_e + 1.0e-2_dp*abs(delta_e)
383 CALL get_selected_ritz_vector(arnoldi_env, 1, arnoldi_matrices(1)%matrix, &
384 pexsi_env%max_ev_vector(ispin))
385 CALL deallocate_arnoldi_env(arnoldi_env)
386 DEALLOCATE (arnoldi_matrices)
387
388 nelectron_exact_pexsi = nelectron_exact
389
390 CALL cp_pexsi_set_options(pexsi_env%options, deltae=delta_e)
391
392 ! Set PEXSI options appropriately for first SCF iteration
393 IF (iscf .EQ. 1) THEN
394 ! Get option isInertiaCount to reset it later on and set it to 1 for first SCF iteration
395 CALL cp_pexsi_get_options(pexsi_env%options, isinertiacount=isinertiacount)
396 CALL cp_pexsi_set_options(pexsi_env%options, isinertiacount=1, &
397 issymbolicfactorize=1)
398 END IF
399
400 ! Write PEXSI options to output
401 CALL cp_pexsi_get_options(pexsi_env%options, isinertiacount=isinertiacount_out, &
402 issymbolicfactorize=is_symbolic_factorize, &
403 mumin0=mu_min_in, mumax0=mu_max_in, &
404 numelectronpexsitolerance=nel_tol)
405
406! IF(unit_nr>0) WRITE(unit_nr,'(/A,I4,A,I4)') " PEXSI| SCF", iscf, &
407! ", spin component", ispin
408
409 IF (unit_nr > 0) WRITE (unit_nr, '(/A,T41,L20)') " PEXSI| Do inertia counting?", &
410 isinertiacount_out .EQ. 1
411 IF (unit_nr > 0) WRITE (unit_nr, '(A,T50,F5.2,T56,F5.2)') &
412 " PEXSI| Guess for min mu, max mu", mu_min_in, mu_max_in
413
414 IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,E20.3)') &
415 " PEXSI| Tolerance in number of electrons", nel_tol
416
417! IF(unit_nr>0) WRITE(unit_nr,'(A,T41,L20)') &
418! " PEXSI| Do symbolic factorization?", is_symbolic_factorize.EQ.1
419
420 IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,F20.2)') &
421 " PEXSI| Arnoldi est. spectral radius", delta_e
422
423 ! Load data into PEXSI
424 CALL cp_pexsi_load_real_hs_matrix( &
425 pexsi_env%plan, &
426 pexsi_env%options, &
427 pexsi_env%csr_mat_ks%nrows_total, &
428 int(pexsi_env%csr_mat_ks%nze_total, kind=int_4), & ! TODO: update when PEXSI changes to Long Int
429 pexsi_env%csr_mat_ks%nze_local, &
430 pexsi_env%csr_mat_ks%nrows_local, &
431 pexsi_env%csr_mat_ks%rowptr_local, &
432 pexsi_env%csr_mat_ks%colind_local, &
433 pexsi_env%csr_mat_ks%nzval_local%r_dp, &
434 s_not_identity, &
435 pexsi_env%csr_mat_s%nzval_local%r_dp)
436
437 ! convert to spin restricted before passing number of electrons to PEXSI
438 CALL convert_nspin_cp2k_pexsi(cp2k_to_pexsi, &
439 numelectron=nelectron_exact_pexsi)
440
441 ! Call DFT driver of PEXSI doing the actual calculation
442 CALL cp_pexsi_dft_driver(pexsi_env%plan, pexsi_env%options, &
443 nelectron_exact_pexsi, mu, nelectron_out, mu_min_out, mu_max_out, &
444 n_total_inertia_iter, n_total_pexsi_iter)
445
446 ! Check convergence
447 nelectron_diff = nelectron_out - nelectron_exact_pexsi
448 pexsi_convergence = abs(nelectron_diff) .LT. nel_tol
449
450 IF (unit_nr > 0) THEN
451 IF (pexsi_convergence) THEN
452 WRITE (unit_nr, '(/A)') " PEXSI| Converged"
453 ELSE
454 WRITE (unit_nr, '(/A)') " PEXSI| PEXSI did not converge!"
455 END IF
456
457! WRITE(unit_nr,'(A,T41,F20.10,T61,F20.10)') " PEXSI| Number of electrons", &
458! nelectron_out/pexsi_env%nspin, nelectron_diff/pexsi_env%nspin
459
460 WRITE (unit_nr, '(A,T41,F20.6)') " PEXSI| Chemical potential", mu
461
462 WRITE (unit_nr, '(A,T41,I20)') " PEXSI| PEXSI iterations", n_total_pexsi_iter
463 WRITE (unit_nr, '(A,T41,I20/)') " PEXSI| Inertia counting iterations", &
464 n_total_inertia_iter
465 END IF
466
467 IF (.NOT. pexsi_convergence) &
468 cpabort("PEXSI did not converge. Consider logPEXSI0 for more information.")
469
470 ! Retrieve results from PEXSI
471 IF (mynode < pexsi_env%mp_dims(1)*pexsi_env%mp_dims(2)) THEN
472 CALL cp_pexsi_retrieve_real_dft_matrix( &
473 pexsi_env%plan, &
474 pexsi_env%csr_mat_p%nzval_local%r_dp, &
475 pexsi_env%csr_mat_E%nzval_local%r_dp, &
476 pexsi_env%csr_mat_F%nzval_local%r_dp, &
477 energy_h, energy_s, free_energy)
478 ! calculate entropic energy contribution -TS = A - U
479 kts = (free_energy - energy_h)
480 END IF
481
482 ! send kTS to all nodes:
483 CALL pexsi_env%mp_group%bcast(kts, 0)
484
485 ! Convert PEXSI CSR matrices to DBCSR matrices
486 CALL dbcsr_convert_csr_to_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
487 pexsi_env%csr_mat_p)
488 CALL dbcsr_copy(matrix_p, pexsi_env%dbcsr_template_matrix_nonsym)
489 CALL dbcsr_convert_csr_to_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
490 pexsi_env%csr_mat_E)
491 CALL dbcsr_copy(matrix_w%matrix, pexsi_env%dbcsr_template_matrix_nonsym)
492
493 ! Convert to spin unrestricted
494 CALL convert_nspin_cp2k_pexsi(pexsi_to_cp2k, matrix_p=matrix_p, &
495 matrix_w=matrix_w, kts=kts)
496
497 ! Pass resulting mu as initial guess for next SCF to PEXSI
498 CALL cp_pexsi_set_options(pexsi_env%options, mu0=mu, mumin0=mu_min_out, &
499 mumax0=mu_max_out)
500
501 ! Reset isInertiaCount according to user input
502 IF (iscf .EQ. 1) THEN
503 CALL cp_pexsi_set_options(pexsi_env%options, isinertiacount= &
504 isinertiacount)
505 END IF
506
507 ! Turn off symbolic factorization for subsequent calls
508 IF (first_call) THEN
509 CALL cp_pexsi_set_options(pexsi_env%options, issymbolicfactorize=0)
510 END IF
511
512 CALL timestop(handle)
513 END SUBROUTINE density_matrix_pexsi
514
515! **************************************************************************************************
516!> \brief Set PEXSI convergence tolerance (numElectronPEXSITolerance), adapted
517!> to the convergence error of the previous SCF step. The tolerance is
518!> calculated using an initial convergence threshold (tol_nel_initial)
519!> and a target convergence threshold (tol_nel_target):
520!> numElectronPEXSITolerance(delta_scf) = alpha*delta_scf+beta
521!> where alpha and beta are chosen such that
522!> numElectronPEXSITolerance(delta_scf_0) = tol_nel_initial
523!> numElectronPEXSITolerance(eps_scf) = tol_nel_target
524!> and delta_scf_0 is the initial SCF convergence error.
525!> \param pexsi_env ...
526!> \param delta_scf convergence error of previous SCF step
527!> \param eps_scf SCF convergence criterion
528!> \param initialize whether or not adaptive thresholing should be initialized
529!> with delta_scf as initial convergence error
530!> \param check_convergence is set to .FALSE. if convergence in number of electrons
531!> will not be achieved in next SCF step
532! **************************************************************************************************
533 SUBROUTINE pexsi_set_convergence_tolerance(pexsi_env, delta_scf, eps_scf, initialize, &
534 check_convergence)
535 TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
536 REAL(kind=dp), INTENT(IN) :: delta_scf, eps_scf
537 LOGICAL, INTENT(IN) :: initialize
538 LOGICAL, INTENT(OUT) :: check_convergence
539
540 CHARACTER(len=*), PARAMETER :: routinen = 'pexsi_set_convergence_tolerance'
541
542 INTEGER :: handle
543 REAL(kind=dp) :: tol_nel
544
545 CALL timeset(routinen, handle)
546
547 tol_nel = pexsi_env%tol_nel_initial
548
549 IF (initialize) THEN
550 pexsi_env%adaptive_nel_alpha = &
551 (pexsi_env%tol_nel_initial - pexsi_env%tol_nel_target)/(abs(delta_scf) - eps_scf)
552 pexsi_env%adaptive_nel_beta = &
553 pexsi_env%tol_nel_initial - pexsi_env%adaptive_nel_alpha*abs(delta_scf)
554 pexsi_env%do_adaptive_tol_nel = .true.
555 END IF
556 IF (pexsi_env%do_adaptive_tol_nel) THEN
557 tol_nel = pexsi_env%adaptive_nel_alpha*abs(delta_scf) + pexsi_env%adaptive_nel_beta
558 tol_nel = max(tol_nel, pexsi_env%tol_nel_target)
559 tol_nel = min(tol_nel, pexsi_env%tol_nel_initial)
560 END IF
561
562 check_convergence = (tol_nel .LE. pexsi_env%tol_nel_target)
563
564 CALL cp_pexsi_set_options(pexsi_env%options, numelectronpexsitolerance=tol_nel)
565 CALL timestop(handle)
566 END SUBROUTINE
567
568! **************************************************************************************************
569!> \brief Pass energy weighted density matrix and entropic energy contribution
570!> to QS environment
571!> \param ls_scf_env ...
572!> \param qs_env ...
573!> \param kTS ...
574!> \param matrix_w ...
575!> \par History
576!> 12.2014 created [Patrick Seewald]
577!> \author Patrick Seewald
578! **************************************************************************************************
579 SUBROUTINE pexsi_to_qs(ls_scf_env, qs_env, kTS, matrix_w)
580 TYPE(ls_scf_env_type) :: ls_scf_env
581 TYPE(qs_environment_type), INTENT(INOUT), POINTER :: qs_env
582 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: kts
583 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
584 OPTIONAL :: matrix_w
585
586 CHARACTER(len=*), PARAMETER :: routinen = 'pexsi_to_qs'
587
588 INTEGER :: handle, ispin, unit_nr
589 REAL(kind=dp) :: kts_total
590 TYPE(cp_logger_type), POINTER :: logger
591 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w_qs
592 TYPE(qs_energy_type), POINTER :: energy
593
594 CALL timeset(routinen, handle)
595
596 NULLIFY (energy)
597
598 ! get a useful output_unit
599 logger => cp_get_default_logger()
600 IF (logger%para_env%is_source()) THEN
601 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
602 ELSE
603 unit_nr = -1
604 END IF
605
606 CALL get_qs_env(qs_env, energy=energy, matrix_w=matrix_w_qs)
607
608 IF (PRESENT(matrix_w)) THEN
609 DO ispin = 1, ls_scf_env%nspins
610 CALL matrix_ls_to_qs(matrix_w_qs(ispin)%matrix, matrix_w(ispin)%matrix, &
611 ls_scf_env%ls_mstruct, covariant=.false.)
612 IF (ls_scf_env%nspins .EQ. 1) CALL dbcsr_scale(matrix_w_qs(ispin)%matrix, 2.0_dp)
613 END DO
614 END IF
615
616 IF (PRESENT(kts)) THEN
617 kts_total = sum(kts)
618 IF (ls_scf_env%nspins .EQ. 1) kts_total = kts_total*2.0_dp
619 energy%kTS = kts_total
620 END IF
621
622 CALL timestop(handle)
623 END SUBROUTINE pexsi_to_qs
624
625END MODULE pexsi_methods
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
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:914
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1165
cp_dbcsr_api::dbcsr_convert_dbcsr_to_csr
subroutine, public dbcsr_convert_dbcsr_to_csr(dbcsr_mat, csr_mat)
...
Definition cp_dbcsr_api.F:349
cp_dbcsr_api::dbcsr_convert_csr_to_dbcsr
subroutine, public dbcsr_convert_csr_to_dbcsr(dbcsr_mat, csr_mat)
...
Definition cp_dbcsr_api.F:333
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:517
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
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:794
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:458
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
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:564
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:306
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:1047
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:580
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:302
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:238
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:535
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, 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, 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, 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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_ks_types
Definition qs_ks_types.F:15
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:180
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
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:217
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:136