(git:58e3e09)
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 
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
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,&
41  USE kinds, ONLY: dp,&
42  int_4,&
43  int_8
52  lib_pexsi_env,&
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 
69 
70 CONTAINS
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
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
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
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 
627 END MODULE pexsi_methods
arnoldi iteration using dbcsr
Definition: arnoldi_api.F:15
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
DBCSR operations in CP2K.
various routines to log and control the output. The idea is that decisions about where to log should ...
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
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
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
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
Types needed for a linear scaling quickstep SCF run based on the density matrix.
objects that represent the structure of input sections and the data contained in an input section
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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public int_8
Definition: kinds.F:54
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public int_4
Definition: kinds.F:51
Interface to the PEXSI library, providing wrappers for all PEXSI routines that are called inside CP2K...
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.
subroutine, public cp_pexsi_retrieve_real_dft_matrix(plan, DMnzvalLocal, EDMnzvalLocal, FDMnzvalLocal, totalEnergyH, totalEnergyS, totalFreeEnergy)
...
subroutine, public cp_pexsi_set_default_options(pexsi_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.
subroutine, public cp_pexsi_dft_driver(plan, pexsi_options, numElectronExact, muPEXSI, numElectronPEXSI, muMinInertia, muMaxInertia, numTotalInertiaIter, numTotalPEXSIIter)
...
subroutine, public cp_pexsi_load_real_hs_matrix(plan, pexsi_options, nrows, nnz, nnzLocal, numColLocal, colptrLocal, rowindLocal, HnzvalLocal, isSIdentity, SnzvalLocal)
...
Methods using the PEXSI library to calculate the density matrix and related quantities using the Kohn...
Definition: pexsi_methods.F:17
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
subroutine, public pexsi_finalize_scf(pexsi_env, mu_spin)
Deallocations and post-processing after SCF.
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.
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...
subroutine, public pexsi_init_scf(ks_env, pexsi_env, template_matrix)
Initializations needed before SCF.
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 ...
Environment storing all data that is needed in order to call the DFT driver of the PEXSI library with...
Definition: pexsi_types.F:18
integer, parameter, public cp2k_to_pexsi
Definition: pexsi_types.F:48
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
integer, parameter, public pexsi_to_cp2k
Definition: pexsi_types.F:48
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.