Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Methods using the PEXSI library to calculate the density matrix and
10 : !> related quantities using the Kohn-Sham and overlap matrices from the
11 : !> linear scaling quickstep SCF environment.
12 : !> \par History
13 : !> 2014.11 created [Patrick Seewald]
14 : !> \author Patrick Seewald
15 : ! **************************************************************************************************
16 :
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) THEN
150 0 : numElectronInitialTolerance = numElectronTargetTolerance
151 : END IF
152 :
153 0 : pexsi_env%tol_nel_initial = numElectronInitialTolerance
154 0 : pexsi_env%tol_nel_target = numElectronTargetTolerance
155 :
156 0 : END SUBROUTINE pexsi_init_read_input
157 :
158 : ! **************************************************************************************************
159 : !> \brief Initializations needed before SCF
160 : !> \param ks_env ...
161 : !> \param pexsi_env ...
162 : !> \param template_matrix DBCSR matrix that defines the block structure and
163 : !> sparsity pattern of all matrices that are sent to PEXSI
164 : ! **************************************************************************************************
165 0 : SUBROUTINE pexsi_init_scf(ks_env, pexsi_env, template_matrix)
166 : TYPE(qs_ks_env_type), POINTER :: ks_env
167 : TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
168 : TYPE(dbcsr_type), INTENT(IN) :: template_matrix
169 :
170 : CHARACTER(len=*), PARAMETER :: routineN = 'pexsi_init_scf'
171 :
172 : INTEGER :: handle, ispin, unit_nr
173 : TYPE(cp_logger_type), POINTER :: logger
174 :
175 0 : CALL timeset(routineN, handle)
176 :
177 0 : logger => cp_get_default_logger()
178 0 : IF (logger%para_env%is_source()) THEN
179 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
180 : ELSE
181 0 : unit_nr = -1
182 : END IF
183 :
184 : ! Create template matrices fixing sparsity pattern for PEXSI
185 :
186 0 : IF (dbcsr_has_symmetry(template_matrix)) THEN
187 : CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, template_matrix, &
188 0 : "symmetric template matrix for CSR conversion")
189 : CALL dbcsr_desymmetrize(pexsi_env%dbcsr_template_matrix_sym, &
190 0 : pexsi_env%dbcsr_template_matrix_nonsym)
191 : ELSE
192 : CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_nonsym, template_matrix, &
193 0 : "non-symmetric template matrix for CSR conversion")
194 : CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, template_matrix, &
195 0 : "symmetric template matrix for CSR conversion")
196 : END IF
197 :
198 : CALL dbcsr_create(pexsi_env%csr_sparsity, "CSR sparsity", &
199 0 : template=pexsi_env%dbcsr_template_matrix_sym)
200 0 : CALL dbcsr_copy(pexsi_env%csr_sparsity, pexsi_env%dbcsr_template_matrix_sym)
201 :
202 0 : CALL cp_dbcsr_to_csr_screening(ks_env, pexsi_env%csr_sparsity)
203 :
204 0 : IF (.NOT. pexsi_env%csr_screening) CALL dbcsr_set(pexsi_env%csr_sparsity, 1.0_dp)
205 : CALL dbcsr_csr_create_from_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
206 : pexsi_env%csr_mat_s, &
207 : dbcsr_csr_eqrow_floor_dist, &
208 : csr_sparsity=pexsi_env%csr_sparsity, &
209 0 : numnodes=pexsi_env%num_ranks_per_pole)
210 :
211 0 : IF (unit_nr > 0) WRITE (unit_nr, "(/T2,A)") "SPARSITY OF THE OVERLAP MATRIX IN CSR FORMAT"
212 0 : CALL dbcsr_csr_print_sparsity(pexsi_env%csr_mat_s, unit_nr)
213 :
214 0 : CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_nonsym, pexsi_env%csr_mat_s)
215 :
216 0 : CALL dbcsr_csr_create(pexsi_env%csr_mat_ks, pexsi_env%csr_mat_s)
217 0 : CALL dbcsr_csr_create(pexsi_env%csr_mat_p, pexsi_env%csr_mat_s)
218 0 : CALL dbcsr_csr_create(pexsi_env%csr_mat_E, pexsi_env%csr_mat_s)
219 0 : CALL dbcsr_csr_create(pexsi_env%csr_mat_F, pexsi_env%csr_mat_s)
220 :
221 0 : 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 0 : template=template_matrix, matrix_type=dbcsr_type_no_symmetry)
225 : END DO
226 :
227 0 : CALL cp_pexsi_set_options(pexsi_env%options, numElectronPEXSITolerance=pexsi_env%tol_nel_initial)
228 :
229 0 : CALL timestop(handle)
230 :
231 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
249 :
250 0 : logger => cp_get_default_logger()
251 0 : IF (logger%para_env%is_source()) THEN
252 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
253 : ELSE
254 : unit_nr = -1
255 : END IF
256 :
257 0 : mu_total = SUM(mu_spin(1:pexsi_env%nspin))/REAL(pexsi_env%nspin, KIND=dp)
258 0 : kTS_total = SUM(pexsi_env%kTS)
259 0 : IF (pexsi_env%nspin == 1) kTS_total = kTS_total*2.0_dp
260 :
261 0 : IF (unit_nr > 0) THEN
262 : WRITE (unit_nr, "(/A,T55,F26.15)") &
263 0 : " PEXSI| Electronic entropic energy (a.u.):", kTS_total
264 : WRITE (unit_nr, "(A,T55,F26.15)") &
265 0 : " PEXSI| Chemical potential (a.u.):", mu_total
266 : END IF
267 :
268 0 : CALL dbcsr_release(pexsi_env%dbcsr_template_matrix_sym)
269 0 : CALL dbcsr_release(pexsi_env%dbcsr_template_matrix_nonsym)
270 0 : CALL dbcsr_release(pexsi_env%csr_sparsity)
271 0 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_p)
272 0 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_ks)
273 0 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_s)
274 0 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_E)
275 0 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_F)
276 0 : DO ispin = 1, pexsi_env%nspin
277 0 : CALL dbcsr_release(pexsi_env%max_ev_vector(ispin))
278 0 : CALL dbcsr_release(pexsi_env%matrix_w(ispin)%matrix)
279 : END DO
280 0 : CALL timestop(handle)
281 0 : pexsi_env%tol_nel_initial = pexsi_env%tol_nel_target ! Turn off adaptive threshold for subsequent SCF cycles
282 0 : 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 0 : SUBROUTINE density_matrix_pexsi(pexsi_env, matrix_p, matrix_w, kTS, matrix_ks, matrix_s, &
302 : nelectron_exact, mu, iscf, ispin)
303 : TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
304 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
305 : TYPE(dbcsr_p_type), INTENT(INOUT) :: matrix_w
306 : REAL(KIND=dp), INTENT(OUT) :: kTS
307 : TYPE(dbcsr_type), INTENT(IN), TARGET :: matrix_ks, matrix_s
308 : INTEGER, INTENT(IN) :: nelectron_exact
309 : REAL(KIND=dp), INTENT(OUT) :: mu
310 : INTEGER, INTENT(IN) :: iscf, ispin
311 :
312 : CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_pexsi'
313 : INTEGER, PARAMETER :: S_not_identity = 0
314 :
315 : INTEGER :: handle, is_symbolic_factorize, isInertiaCount, isInertiaCount_out, mynode, &
316 : n_total_inertia_iter, n_total_pexsi_iter, unit_nr
317 : LOGICAL :: first_call, pexsi_convergence
318 : REAL(KIND=dp) :: delta_E, energy_H, energy_S, free_energy, mu_max_in, mu_max_out, mu_min_in, &
319 : mu_min_out, nel_tol, nelectron_diff, nelectron_exact_pexsi, nelectron_out
320 : TYPE(arnoldi_env_type) :: arnoldi_env
321 : TYPE(cp_logger_type), POINTER :: logger
322 : TYPE(dbcsr_distribution_type) :: dist
323 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: arnoldi_matrices
324 :
325 0 : CALL timeset(routineN, handle)
326 :
327 : ! get a useful output_unit
328 0 : logger => cp_get_default_logger()
329 0 : IF (logger%para_env%is_source()) THEN
330 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
331 : ELSE
332 : unit_nr = -1
333 : END IF
334 :
335 0 : first_call = (iscf == 1) .AND. (ispin == 1)
336 :
337 : ! Assert a few things the first time PEXSI is called
338 : IF (first_call) THEN
339 : ! Assertion that matrices have the expected symmetry (both should be symmetric if no
340 : ! S preconditioning and no molecular clustering)
341 0 : IF (.NOT. dbcsr_has_symmetry(matrix_ks)) THEN
342 0 : CPABORT("PEXSI interface expects a non-symmetric DBCSR Kohn-Sham matrix")
343 : END IF
344 0 : IF (.NOT. dbcsr_has_symmetry(matrix_s)) THEN
345 0 : CPABORT("PEXSI interface expects a non-symmetric DBCSR overlap matrix")
346 : END IF
347 :
348 : ! Assertion on datatype
349 : IF ((pexsi_env%csr_mat_s%nzval_local%data_type /= dbcsr_csr_type_real_8) &
350 0 : .OR. (pexsi_env%csr_mat_ks%nzval_local%data_type /= dbcsr_csr_type_real_8)) THEN
351 0 : CPABORT("Complex data type not supported by PEXSI")
352 : END IF
353 :
354 : ! Assertion on number of non-zero elements
355 : !(TODO: update when PEXSI changes to Long Int)
356 0 : IF (pexsi_env%csr_mat_s%nze_total >= INT(2, kind=int_8)**31) THEN
357 0 : CPABORT("Total number of non-zero elements of CSR matrix is too large to be handled by PEXSI")
358 : END IF
359 : END IF
360 :
361 0 : CALL dbcsr_get_info(matrix_ks, distribution=dist)
362 0 : CALL dbcsr_distribution_get(dist, mynode=mynode)
363 :
364 : ! Convert DBCSR matrices to PEXSI CSR format. Intermediate step to template matrix
365 : ! needed in order to retain the initial sparsity pattern that is required for the
366 : ! conversion to CSR format.
367 0 : CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, matrix_s, keep_sparsity=.TRUE.)
368 : CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_sym, &
369 0 : pexsi_env%csr_mat_s)
370 :
371 0 : CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, matrix_ks, keep_sparsity=.TRUE.)
372 : CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_sym, &
373 0 : pexsi_env%csr_mat_ks)
374 :
375 : ! Get PEXSI input delta_E (upper bound for largest eigenvalue) using Arnoldi
376 0 : NULLIFY (arnoldi_matrices)
377 0 : CALL dbcsr_allocate_matrix_set(arnoldi_matrices, 2)
378 0 : arnoldi_matrices(1)%matrix => matrix_ks
379 0 : arnoldi_matrices(2)%matrix => matrix_s
380 : CALL setup_arnoldi_env(arnoldi_env, arnoldi_matrices, max_iter=20, &
381 : threshold=1.0E-2_dp, selection_crit=2, nval_request=1, nrestarts=21, &
382 0 : generalized_ev=.TRUE., iram=.FALSE.)
383 0 : IF (iscf > 1) CALL set_arnoldi_initial_vector(arnoldi_env, pexsi_env%max_ev_vector(ispin))
384 0 : CALL arnoldi_ev(arnoldi_matrices, arnoldi_env)
385 0 : delta_E = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
386 : ! increase delta_E a bit to make sure that it really is an upper bound
387 0 : delta_E = delta_E + 1.0E-2_dp*ABS(delta_E)
388 : CALL get_selected_ritz_vector(arnoldi_env, 1, arnoldi_matrices(1)%matrix, &
389 0 : pexsi_env%max_ev_vector(ispin))
390 0 : CALL deallocate_arnoldi_env(arnoldi_env)
391 0 : DEALLOCATE (arnoldi_matrices)
392 :
393 0 : nelectron_exact_pexsi = nelectron_exact
394 :
395 0 : CALL cp_pexsi_set_options(pexsi_env%options, deltaE=delta_E)
396 :
397 : ! Set PEXSI options appropriately for first SCF iteration
398 0 : IF (iscf == 1) THEN
399 : ! Get option isInertiaCount to reset it later on and set it to 1 for first SCF iteration
400 0 : CALL cp_pexsi_get_options(pexsi_env%options, isInertiaCount=isInertiaCount)
401 : CALL cp_pexsi_set_options(pexsi_env%options, isInertiaCount=1, &
402 0 : isSymbolicFactorize=1)
403 : END IF
404 :
405 : ! Write PEXSI options to output
406 : CALL cp_pexsi_get_options(pexsi_env%options, isInertiaCount=isInertiaCount_out, &
407 : isSymbolicFactorize=is_symbolic_factorize, &
408 : muMin0=mu_min_in, muMax0=mu_max_in, &
409 0 : NumElectronPEXSITolerance=nel_tol)
410 :
411 : ! IF(unit_nr>0) WRITE(unit_nr,'(/A,I4,A,I4)') " PEXSI| SCF", iscf, &
412 : ! ", spin component", ispin
413 :
414 0 : IF (unit_nr > 0) WRITE (unit_nr, '(/A,T41,L20)') " PEXSI| Do inertia counting?", &
415 0 : isInertiaCount_out == 1
416 : IF (unit_nr > 0) WRITE (unit_nr, '(A,T50,F5.2,T56,F5.2)') &
417 0 : " PEXSI| Guess for min mu, max mu", mu_min_in, mu_max_in
418 :
419 : IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,E20.3)') &
420 0 : " PEXSI| Tolerance in number of electrons", nel_tol
421 :
422 : ! IF(unit_nr>0) WRITE(unit_nr,'(A,T41,L20)') &
423 : ! " PEXSI| Do symbolic factorization?", is_symbolic_factorize==1
424 :
425 : IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,F20.2)') &
426 0 : " PEXSI| Arnoldi est. spectral radius", delta_E
427 :
428 : ! Load data into PEXSI
429 : CALL cp_pexsi_load_real_hs_matrix( &
430 : pexsi_env%plan, &
431 : pexsi_env%options, &
432 : pexsi_env%csr_mat_ks%nrows_total, &
433 : INT(pexsi_env%csr_mat_ks%nze_total, kind=int_4), & ! TODO: update when PEXSI changes to Long Int
434 : pexsi_env%csr_mat_ks%nze_local, &
435 : pexsi_env%csr_mat_ks%nrows_local, &
436 : pexsi_env%csr_mat_ks%rowptr_local, &
437 : pexsi_env%csr_mat_ks%colind_local, &
438 : pexsi_env%csr_mat_ks%nzval_local%r_dp, &
439 : S_not_identity, &
440 0 : pexsi_env%csr_mat_s%nzval_local%r_dp)
441 :
442 : ! convert to spin restricted before passing number of electrons to PEXSI
443 : CALL convert_nspin_cp2k_pexsi(cp2k_to_pexsi, &
444 0 : numElectron=nelectron_exact_pexsi)
445 :
446 : ! Call DFT driver of PEXSI doing the actual calculation
447 : CALL cp_pexsi_dft_driver(pexsi_env%plan, pexsi_env%options, &
448 : nelectron_exact_pexsi, mu, nelectron_out, mu_min_out, mu_max_out, &
449 0 : n_total_inertia_iter, n_total_pexsi_iter)
450 :
451 : ! Check convergence
452 0 : nelectron_diff = nelectron_out - nelectron_exact_pexsi
453 0 : pexsi_convergence = ABS(nelectron_diff) < nel_tol
454 :
455 0 : IF (unit_nr > 0) THEN
456 0 : IF (pexsi_convergence) THEN
457 0 : WRITE (unit_nr, '(/A)') " PEXSI| Converged"
458 : ELSE
459 0 : WRITE (unit_nr, '(/A)') " PEXSI| PEXSI did not converge!"
460 : END IF
461 :
462 : ! WRITE(unit_nr,'(A,T41,F20.10,T61,F20.10)') " PEXSI| Number of electrons", &
463 : ! nelectron_out/pexsi_env%nspin, nelectron_diff/pexsi_env%nspin
464 :
465 0 : WRITE (unit_nr, '(A,T41,F20.6)') " PEXSI| Chemical potential", mu
466 :
467 0 : WRITE (unit_nr, '(A,T41,I20)') " PEXSI| PEXSI iterations", n_total_pexsi_iter
468 0 : WRITE (unit_nr, '(A,T41,I20/)') " PEXSI| Inertia counting iterations", &
469 0 : n_total_inertia_iter
470 : END IF
471 :
472 0 : IF (.NOT. pexsi_convergence) THEN
473 0 : CPABORT("PEXSI did not converge. Consider logPEXSI0 for more information.")
474 : END IF
475 :
476 : ! Retrieve results from PEXSI
477 0 : IF (mynode < pexsi_env%mp_dims(1)*pexsi_env%mp_dims(2)) THEN
478 : CALL cp_pexsi_retrieve_real_dft_matrix( &
479 : pexsi_env%plan, &
480 : pexsi_env%csr_mat_p%nzval_local%r_dp, &
481 : pexsi_env%csr_mat_E%nzval_local%r_dp, &
482 : pexsi_env%csr_mat_F%nzval_local%r_dp, &
483 0 : energy_H, energy_S, free_energy)
484 : ! calculate entropic energy contribution -TS = A - U
485 0 : kTS = (free_energy - energy_H)
486 : END IF
487 :
488 : ! send kTS to all nodes:
489 0 : CALL pexsi_env%mp_group%bcast(kTS, 0)
490 :
491 : ! Convert PEXSI CSR matrices to DBCSR matrices
492 : CALL dbcsr_convert_csr_to_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
493 0 : pexsi_env%csr_mat_p)
494 0 : CALL dbcsr_copy(matrix_p, pexsi_env%dbcsr_template_matrix_nonsym)
495 : CALL dbcsr_convert_csr_to_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
496 0 : pexsi_env%csr_mat_E)
497 0 : CALL dbcsr_copy(matrix_w%matrix, pexsi_env%dbcsr_template_matrix_nonsym)
498 :
499 : ! Convert to spin unrestricted
500 : CALL convert_nspin_cp2k_pexsi(pexsi_to_cp2k, matrix_p=matrix_p, &
501 0 : matrix_w=matrix_w, kTS=kTS)
502 :
503 : ! Pass resulting mu as initial guess for next SCF to PEXSI
504 : CALL cp_pexsi_set_options(pexsi_env%options, mu0=mu, muMin0=mu_min_out, &
505 0 : muMax0=mu_max_out)
506 :
507 : ! Reset isInertiaCount according to user input
508 0 : IF (iscf == 1) THEN
509 : CALL cp_pexsi_set_options(pexsi_env%options, isInertiaCount= &
510 0 : isInertiaCount)
511 : END IF
512 :
513 : ! Turn off symbolic factorization for subsequent calls
514 0 : IF (first_call) THEN
515 0 : CALL cp_pexsi_set_options(pexsi_env%options, isSymbolicFactorize=0)
516 : END IF
517 :
518 0 : CALL timestop(handle)
519 0 : END SUBROUTINE density_matrix_pexsi
520 :
521 : ! **************************************************************************************************
522 : !> \brief Set PEXSI convergence tolerance (numElectronPEXSITolerance), adapted
523 : !> to the convergence error of the previous SCF step. The tolerance is
524 : !> calculated using an initial convergence threshold (tol_nel_initial)
525 : !> and a target convergence threshold (tol_nel_target):
526 : !> numElectronPEXSITolerance(delta_scf) = alpha*delta_scf+beta
527 : !> where alpha and beta are chosen such that
528 : !> numElectronPEXSITolerance(delta_scf_0) = tol_nel_initial
529 : !> numElectronPEXSITolerance(eps_scf) = tol_nel_target
530 : !> and delta_scf_0 is the initial SCF convergence error.
531 : !> \param pexsi_env ...
532 : !> \param delta_scf convergence error of previous SCF step
533 : !> \param eps_scf SCF convergence criterion
534 : !> \param initialize whether or not adaptive thresholing should be initialized
535 : !> with delta_scf as initial convergence error
536 : !> \param check_convergence is set to .FALSE. if convergence in number of electrons
537 : !> will not be achieved in next SCF step
538 : ! **************************************************************************************************
539 0 : SUBROUTINE pexsi_set_convergence_tolerance(pexsi_env, delta_scf, eps_scf, initialize, &
540 : check_convergence)
541 : TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
542 : REAL(KIND=dp), INTENT(IN) :: delta_scf, eps_scf
543 : LOGICAL, INTENT(IN) :: initialize
544 : LOGICAL, INTENT(OUT) :: check_convergence
545 :
546 : CHARACTER(len=*), PARAMETER :: routineN = 'pexsi_set_convergence_tolerance'
547 :
548 : INTEGER :: handle
549 : REAL(KIND=dp) :: tol_nel
550 :
551 0 : CALL timeset(routineN, handle)
552 :
553 0 : tol_nel = pexsi_env%tol_nel_initial
554 :
555 0 : IF (initialize) THEN
556 : pexsi_env%adaptive_nel_alpha = &
557 0 : (pexsi_env%tol_nel_initial - pexsi_env%tol_nel_target)/(ABS(delta_scf) - eps_scf)
558 : pexsi_env%adaptive_nel_beta = &
559 0 : pexsi_env%tol_nel_initial - pexsi_env%adaptive_nel_alpha*ABS(delta_scf)
560 0 : pexsi_env%do_adaptive_tol_nel = .TRUE.
561 : END IF
562 0 : IF (pexsi_env%do_adaptive_tol_nel) THEN
563 0 : tol_nel = pexsi_env%adaptive_nel_alpha*ABS(delta_scf) + pexsi_env%adaptive_nel_beta
564 0 : tol_nel = MAX(tol_nel, pexsi_env%tol_nel_target)
565 0 : tol_nel = MIN(tol_nel, pexsi_env%tol_nel_initial)
566 : END IF
567 :
568 0 : check_convergence = (tol_nel <= pexsi_env%tol_nel_target)
569 :
570 0 : CALL cp_pexsi_set_options(pexsi_env%options, numElectronPEXSITolerance=tol_nel)
571 0 : CALL timestop(handle)
572 0 : END SUBROUTINE pexsi_set_convergence_tolerance
573 :
574 : ! **************************************************************************************************
575 : !> \brief Pass energy weighted density matrix and entropic energy contribution
576 : !> to QS environment
577 : !> \param ls_scf_env ...
578 : !> \param qs_env ...
579 : !> \param kTS ...
580 : !> \param matrix_w ...
581 : !> \par History
582 : !> 12.2014 created [Patrick Seewald]
583 : !> \author Patrick Seewald
584 : ! **************************************************************************************************
585 0 : SUBROUTINE pexsi_to_qs(ls_scf_env, qs_env, kTS, matrix_w)
586 : TYPE(ls_scf_env_type) :: ls_scf_env
587 : TYPE(qs_environment_type), INTENT(INOUT), POINTER :: qs_env
588 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: kTS
589 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
590 : OPTIONAL :: matrix_w
591 :
592 : CHARACTER(len=*), PARAMETER :: routineN = 'pexsi_to_qs'
593 :
594 : INTEGER :: handle, ispin, unit_nr
595 : REAL(KIND=dp) :: kTS_total
596 : TYPE(cp_logger_type), POINTER :: logger
597 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w_qs
598 : TYPE(qs_energy_type), POINTER :: energy
599 :
600 0 : CALL timeset(routineN, handle)
601 :
602 0 : NULLIFY (energy)
603 :
604 : ! get a useful output_unit
605 0 : logger => cp_get_default_logger()
606 0 : IF (logger%para_env%is_source()) THEN
607 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
608 : ELSE
609 : unit_nr = -1
610 : END IF
611 :
612 0 : CALL get_qs_env(qs_env, energy=energy, matrix_w=matrix_w_qs)
613 :
614 0 : IF (PRESENT(matrix_w)) THEN
615 0 : DO ispin = 1, ls_scf_env%nspins
616 : CALL matrix_ls_to_qs(matrix_w_qs(ispin)%matrix, matrix_w(ispin)%matrix, &
617 0 : ls_scf_env%ls_mstruct, covariant=.FALSE.)
618 0 : IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(matrix_w_qs(ispin)%matrix, 2.0_dp)
619 : END DO
620 : END IF
621 :
622 0 : IF (PRESENT(kTS)) THEN
623 0 : kTS_total = SUM(kTS)
624 0 : IF (ls_scf_env%nspins == 1) kTS_total = kTS_total*2.0_dp
625 0 : energy%kTS = kTS_total
626 : END IF
627 :
628 0 : CALL timestop(handle)
629 0 : END SUBROUTINE pexsi_to_qs
630 :
631 : END MODULE pexsi_methods
|
