Line data Source code
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 :
17 : MODULE 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 :
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 0 : 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 0 : r_val=temperature)
99 : CALL section_vals_val_get(pexsi_section, "GAP", &
100 0 : r_val=gap)
101 : CALL section_vals_val_get(pexsi_section, "NUM_POLE", &
102 0 : i_val=numPole)
103 : CALL section_vals_val_get(pexsi_section, "IS_INERTIA_COUNT", &
104 0 : l_val=isInertiaCount)
105 : CALL section_vals_val_get(pexsi_section, "MAX_PEXSI_ITER", &
106 0 : i_val=maxPEXSIIter)
107 : CALL section_vals_val_get(pexsi_section, "MU_MIN_0", &
108 0 : r_val=muMin0)
109 : CALL section_vals_val_get(pexsi_section, "MU_MAX_0", &
110 0 : r_val=muMax0)
111 : CALL section_vals_val_get(pexsi_section, "MU_INERTIA_TOLERANCE", &
112 0 : r_val=muInertiaTolerance)
113 : CALL section_vals_val_get(pexsi_section, "MU_INERTIA_EXPANSION", &
114 0 : r_val=muInertiaExpansion)
115 : CALL section_vals_val_get(pexsi_section, "MU_PEXSI_SAFE_GUARD", &
116 0 : r_val=muPEXSISafeGuard)
117 : CALL section_vals_val_get(pexsi_section, "NUM_ELECTRON_INITIAL_TOLERANCE", &
118 0 : r_val=numElectronInitialTolerance)
119 : CALL section_vals_val_get(pexsi_section, "NUM_ELECTRON_PEXSI_TOLERANCE", &
120 0 : r_val=numElectronTargetTolerance)
121 : CALL section_vals_val_get(pexsi_section, "ORDERING", &
122 0 : i_val=ordering)
123 : CALL section_vals_val_get(pexsi_section, "ROW_ORDERING", &
124 0 : i_val=rowOrdering)
125 : CALL section_vals_val_get(pexsi_section, "NP_SYMB_FACT", &
126 0 : i_val=npSymbFact)
127 : CALL section_vals_val_get(pexsi_section, "VERBOSITY", &
128 0 : i_val=verbosity)
129 : CALL section_vals_val_get(pexsi_section, "MIN_RANKS_PER_POLE", &
130 0 : i_val=min_ranks_per_pole)
131 : CALL section_vals_val_get(pexsi_section, "CSR_SCREENING", &
132 0 : l_val=csr_screening)
133 :
134 0 : isInertiaCount_int = MERGE(1, 0, isInertiaCount) ! is integer in PEXSI
135 :
136 : ! Set default options inside PEXSI
137 0 : 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 0 : ordering=ordering, rowOrdering=rowOrdering, npSymbFact=npSymbFact, verbosity=verbosity)
145 :
146 0 : pexsi_env%num_ranks_per_pole = min_ranks_per_pole ! not a PEXSI option
147 0 : pexsi_env%csr_screening = csr_screening
148 :
149 0 : IF (numElectronInitialTolerance < numElectronTargetTolerance) &
150 0 : numElectronInitialTolerance = numElectronTargetTolerance
151 :
152 0 : pexsi_env%tol_nel_initial = numElectronInitialTolerance
153 0 : pexsi_env%tol_nel_target = numElectronTargetTolerance
154 :
155 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
175 :
176 0 : logger => cp_get_default_logger()
177 0 : IF (logger%para_env%is_source()) THEN
178 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
179 : ELSE
180 0 : unit_nr = -1
181 : END IF
182 :
183 : ! Create template matrices fixing sparsity pattern for PEXSI
184 :
185 0 : IF (dbcsr_has_symmetry(template_matrix)) THEN
186 : CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, template_matrix, &
187 0 : "symmetric template matrix for CSR conversion")
188 : CALL dbcsr_desymmetrize(pexsi_env%dbcsr_template_matrix_sym, &
189 0 : pexsi_env%dbcsr_template_matrix_nonsym)
190 : ELSE
191 : CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_nonsym, template_matrix, &
192 0 : "non-symmetric template matrix for CSR conversion")
193 : CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, template_matrix, &
194 0 : "symmetric template matrix for CSR conversion")
195 : END IF
196 :
197 : CALL dbcsr_create(pexsi_env%csr_sparsity, "CSR sparsity", &
198 0 : template=pexsi_env%dbcsr_template_matrix_sym)
199 0 : CALL dbcsr_copy(pexsi_env%csr_sparsity, pexsi_env%dbcsr_template_matrix_sym)
200 :
201 0 : CALL cp_dbcsr_to_csr_screening(ks_env, pexsi_env%csr_sparsity)
202 :
203 0 : 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 0 : numnodes=pexsi_env%num_ranks_per_pole)
209 :
210 0 : IF (unit_nr > 0) WRITE (unit_nr, "(/T2,A)") "SPARSITY OF THE OVERLAP MATRIX IN CSR FORMAT"
211 0 : CALL dbcsr_csr_print_sparsity(pexsi_env%csr_mat_s, unit_nr)
212 :
213 0 : CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_nonsym, pexsi_env%csr_mat_s)
214 :
215 0 : CALL dbcsr_csr_create(pexsi_env%csr_mat_ks, pexsi_env%csr_mat_s)
216 0 : CALL dbcsr_csr_create(pexsi_env%csr_mat_p, pexsi_env%csr_mat_s)
217 0 : CALL dbcsr_csr_create(pexsi_env%csr_mat_E, pexsi_env%csr_mat_s)
218 0 : CALL dbcsr_csr_create(pexsi_env%csr_mat_F, pexsi_env%csr_mat_s)
219 :
220 0 : 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 0 : template=template_matrix, matrix_type=dbcsr_type_no_symmetry)
224 : END DO
225 :
226 0 : CALL cp_pexsi_set_options(pexsi_env%options, numElectronPEXSITolerance=pexsi_env%tol_nel_initial)
227 :
228 0 : CALL timestop(handle)
229 :
230 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
248 :
249 0 : logger => cp_get_default_logger()
250 0 : IF (logger%para_env%is_source()) THEN
251 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
252 : ELSE
253 : unit_nr = -1
254 : END IF
255 :
256 0 : mu_total = SUM(mu_spin(1:pexsi_env%nspin))/REAL(pexsi_env%nspin, KIND=dp)
257 0 : kTS_total = SUM(pexsi_env%kTS)
258 0 : IF (pexsi_env%nspin == 1) kTS_total = kTS_total*2.0_dp
259 :
260 0 : IF (unit_nr > 0) THEN
261 : WRITE (unit_nr, "(/A,T55,F26.15)") &
262 0 : " PEXSI| Electronic entropic energy (a.u.):", kTS_total
263 : WRITE (unit_nr, "(A,T55,F26.15)") &
264 0 : " PEXSI| Chemical potential (a.u.):", mu_total
265 : END IF
266 :
267 0 : CALL dbcsr_release(pexsi_env%dbcsr_template_matrix_sym)
268 0 : CALL dbcsr_release(pexsi_env%dbcsr_template_matrix_nonsym)
269 0 : CALL dbcsr_release(pexsi_env%csr_sparsity)
270 0 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_p)
271 0 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_ks)
272 0 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_s)
273 0 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_E)
274 0 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_F)
275 0 : DO ispin = 1, pexsi_env%nspin
276 0 : CALL dbcsr_release(pexsi_env%max_ev_vector(ispin))
277 0 : CALL dbcsr_release(pexsi_env%matrix_w(ispin)%matrix)
278 : END DO
279 0 : CALL timestop(handle)
280 0 : pexsi_env%tol_nel_initial = pexsi_env%tol_nel_target ! Turn off adaptive threshold for subsequent SCF cycles
281 0 : 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 0 : 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 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: arnoldi_matrices
323 :
324 0 : CALL timeset(routineN, handle)
325 :
326 : ! get a useful output_unit
327 0 : logger => cp_get_default_logger()
328 0 : IF (logger%para_env%is_source()) THEN
329 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
330 : ELSE
331 : unit_nr = -1
332 : END IF
333 :
334 0 : first_call = (iscf == 1) .AND. (ispin == 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 0 : IF (.NOT. dbcsr_has_symmetry(matrix_ks)) &
341 0 : CPABORT("PEXSI interface expects a non-symmetric DBCSR Kohn-Sham matrix")
342 0 : IF (.NOT. dbcsr_has_symmetry(matrix_s)) &
343 0 : 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 /= dbcsr_csr_type_real_8) &
347 0 : .OR. (pexsi_env%csr_mat_ks%nzval_local%data_type /= dbcsr_csr_type_real_8)) &
348 0 : 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 0 : IF (pexsi_env%csr_mat_s%nze_total >= INT(2, kind=int_8)**31) &
353 0 : CPABORT("Total number of non-zero elements of CSR matrix is too large to be handled by PEXSI")
354 : END IF
355 :
356 0 : CALL dbcsr_get_info(matrix_ks, distribution=dist)
357 0 : 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 0 : 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 0 : pexsi_env%csr_mat_s)
365 :
366 0 : 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 0 : pexsi_env%csr_mat_ks)
369 :
370 : ! Get PEXSI input delta_E (upper bound for largest eigenvalue) using Arnoldi
371 0 : NULLIFY (arnoldi_matrices)
372 0 : CALL dbcsr_allocate_matrix_set(arnoldi_matrices, 2)
373 0 : arnoldi_matrices(1)%matrix => matrix_ks
374 0 : 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 0 : generalized_ev=.TRUE., iram=.FALSE.)
378 0 : IF (iscf > 1) CALL set_arnoldi_initial_vector(arnoldi_env, pexsi_env%max_ev_vector(ispin))
379 0 : CALL arnoldi_ev(arnoldi_matrices, arnoldi_env)
380 0 : 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 0 : 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 0 : pexsi_env%max_ev_vector(ispin))
385 0 : CALL deallocate_arnoldi_env(arnoldi_env)
386 0 : DEALLOCATE (arnoldi_matrices)
387 :
388 0 : nelectron_exact_pexsi = nelectron_exact
389 :
390 0 : CALL cp_pexsi_set_options(pexsi_env%options, deltaE=delta_E)
391 :
392 : ! Set PEXSI options appropriately for first SCF iteration
393 0 : IF (iscf == 1) THEN
394 : ! Get option isInertiaCount to reset it later on and set it to 1 for first SCF iteration
395 0 : CALL cp_pexsi_get_options(pexsi_env%options, isInertiaCount=isInertiaCount)
396 : CALL cp_pexsi_set_options(pexsi_env%options, isInertiaCount=1, &
397 0 : 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 0 : 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 0 : IF (unit_nr > 0) WRITE (unit_nr, '(/A,T41,L20)') " PEXSI| Do inertia counting?", &
410 0 : isInertiaCount_out == 1
411 : IF (unit_nr > 0) WRITE (unit_nr, '(A,T50,F5.2,T56,F5.2)') &
412 0 : " 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 0 : " 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==1
419 :
420 : IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,F20.2)') &
421 0 : " 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 0 : 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 0 : 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 0 : n_total_inertia_iter, n_total_pexsi_iter)
445 :
446 : ! Check convergence
447 0 : nelectron_diff = nelectron_out - nelectron_exact_pexsi
448 0 : pexsi_convergence = ABS(nelectron_diff) < nel_tol
449 :
450 0 : IF (unit_nr > 0) THEN
451 0 : IF (pexsi_convergence) THEN
452 0 : WRITE (unit_nr, '(/A)') " PEXSI| Converged"
453 : ELSE
454 0 : 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 0 : WRITE (unit_nr, '(A,T41,F20.6)') " PEXSI| Chemical potential", mu
461 :
462 0 : WRITE (unit_nr, '(A,T41,I20)') " PEXSI| PEXSI iterations", n_total_pexsi_iter
463 0 : WRITE (unit_nr, '(A,T41,I20/)') " PEXSI| Inertia counting iterations", &
464 0 : n_total_inertia_iter
465 : END IF
466 :
467 0 : IF (.NOT. pexsi_convergence) &
468 0 : CPABORT("PEXSI did not converge. Consider logPEXSI0 for more information.")
469 :
470 : ! Retrieve results from PEXSI
471 0 : 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 0 : energy_H, energy_S, free_energy)
478 : ! calculate entropic energy contribution -TS = A - U
479 0 : kTS = (free_energy - energy_H)
480 : END IF
481 :
482 : ! send kTS to all nodes:
483 0 : 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 0 : pexsi_env%csr_mat_p)
488 0 : 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 0 : pexsi_env%csr_mat_E)
491 0 : 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 0 : 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 0 : muMax0=mu_max_out)
500 :
501 : ! Reset isInertiaCount according to user input
502 0 : IF (iscf == 1) THEN
503 : CALL cp_pexsi_set_options(pexsi_env%options, isInertiaCount= &
504 0 : isInertiaCount)
505 : END IF
506 :
507 : ! Turn off symbolic factorization for subsequent calls
508 0 : IF (first_call) THEN
509 0 : CALL cp_pexsi_set_options(pexsi_env%options, isSymbolicFactorize=0)
510 : END IF
511 :
512 0 : CALL timestop(handle)
513 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
546 :
547 0 : tol_nel = pexsi_env%tol_nel_initial
548 :
549 0 : IF (initialize) THEN
550 : pexsi_env%adaptive_nel_alpha = &
551 0 : (pexsi_env%tol_nel_initial - pexsi_env%tol_nel_target)/(ABS(delta_scf) - eps_scf)
552 : pexsi_env%adaptive_nel_beta = &
553 0 : pexsi_env%tol_nel_initial - pexsi_env%adaptive_nel_alpha*ABS(delta_scf)
554 0 : pexsi_env%do_adaptive_tol_nel = .TRUE.
555 : END IF
556 0 : IF (pexsi_env%do_adaptive_tol_nel) THEN
557 0 : tol_nel = pexsi_env%adaptive_nel_alpha*ABS(delta_scf) + pexsi_env%adaptive_nel_beta
558 0 : tol_nel = MAX(tol_nel, pexsi_env%tol_nel_target)
559 0 : tol_nel = MIN(tol_nel, pexsi_env%tol_nel_initial)
560 : END IF
561 :
562 0 : check_convergence = (tol_nel <= pexsi_env%tol_nel_target)
563 :
564 0 : CALL cp_pexsi_set_options(pexsi_env%options, numElectronPEXSITolerance=tol_nel)
565 0 : CALL timestop(handle)
566 0 : END SUBROUTINE pexsi_set_convergence_tolerance
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 0 : 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 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w_qs
592 : TYPE(qs_energy_type), POINTER :: energy
593 :
594 0 : CALL timeset(routineN, handle)
595 :
596 0 : NULLIFY (energy)
597 :
598 : ! get a useful output_unit
599 0 : logger => cp_get_default_logger()
600 0 : IF (logger%para_env%is_source()) THEN
601 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
602 : ELSE
603 : unit_nr = -1
604 : END IF
605 :
606 0 : CALL get_qs_env(qs_env, energy=energy, matrix_w=matrix_w_qs)
607 :
608 0 : IF (PRESENT(matrix_w)) THEN
609 0 : DO ispin = 1, ls_scf_env%nspins
610 : CALL matrix_ls_to_qs(matrix_w_qs(ispin)%matrix, matrix_w(ispin)%matrix, &
611 0 : ls_scf_env%ls_mstruct, covariant=.FALSE.)
612 0 : IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(matrix_w_qs(ispin)%matrix, 2.0_dp)
613 : END DO
614 : END IF
615 :
616 0 : IF (PRESENT(kTS)) THEN
617 0 : kTS_total = SUM(kTS)
618 0 : IF (ls_scf_env%nspins == 1) kTS_total = kTS_total*2.0_dp
619 0 : energy%kTS = kTS_total
620 : END IF
621 :
622 0 : CALL timestop(handle)
623 0 : END SUBROUTINE pexsi_to_qs
624 :
625 : END MODULE pexsi_methods
|
