Line data Source code
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 :
17 : MODULE pexsi_methods
18 : USE arnoldi_api, ONLY: arnoldi_data_type,&
19 : arnoldi_ev,&
20 : deallocate_arnoldi_data,&
21 : get_selected_ritz_val,&
22 : get_selected_ritz_vector,&
23 : set_arnoldi_initial_vector,&
24 : setup_arnoldi_data
25 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_to_csr_screening
26 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
27 : USE cp_log_handling, ONLY: cp_get_default_logger,&
28 : cp_logger_get_default_unit_nr,&
29 : cp_logger_type
30 : USE dbcsr_api, ONLY: &
31 : dbcsr_convert_csr_to_dbcsr, dbcsr_convert_dbcsr_to_csr, dbcsr_copy, &
32 : dbcsr_copy_into_existing, dbcsr_create, dbcsr_csr_create, dbcsr_csr_create_from_dbcsr, &
33 : dbcsr_csr_destroy, dbcsr_csr_eqrow_floor_dist, dbcsr_csr_print_sparsity, &
34 : dbcsr_desymmetrize, dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_info, &
35 : dbcsr_has_symmetry, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, &
36 : dbcsr_type_no_symmetry, dbcsr_type_real_8
37 : USE dm_ls_scf_qs, ONLY: matrix_ls_to_qs
38 : USE dm_ls_scf_types, ONLY: ls_scf_env_type
39 : USE input_section_types, ONLY: section_vals_type,&
40 : section_vals_val_get
41 : USE kinds, ONLY: dp,&
42 : int_4,&
43 : int_8
44 : USE pexsi_interface, ONLY: cp_pexsi_dft_driver,&
45 : cp_pexsi_get_options,&
46 : cp_pexsi_load_real_hs_matrix,&
47 : cp_pexsi_retrieve_real_dft_matrix,&
48 : cp_pexsi_set_default_options,&
49 : cp_pexsi_set_options
50 : USE pexsi_types, ONLY: convert_nspin_cp2k_pexsi,&
51 : cp2k_to_pexsi,&
52 : lib_pexsi_env,&
53 : pexsi_to_cp2k
54 : USE qs_energy_types, ONLY: qs_energy_type
55 : USE qs_environment_types, ONLY: get_qs_env,&
56 : qs_environment_type
57 : USE qs_ks_types, ONLY: qs_ks_env_type
58 : #include "./base/base_uses.f90"
59 :
60 : IMPLICIT NONE
61 : PRIVATE
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pexsi_methods'
64 :
65 : LOGICAL, PARAMETER, PRIVATE :: careful_mod = .FALSE.
66 :
67 : PUBLIC :: density_matrix_pexsi, pexsi_init_read_input, pexsi_to_qs, pexsi_init_scf, pexsi_finalize_scf, &
68 : pexsi_set_convergence_tolerance
69 :
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 144 : 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 8 : r_val=temperature)
99 : CALL section_vals_val_get(pexsi_section, "GAP", &
100 8 : r_val=gap)
101 : CALL section_vals_val_get(pexsi_section, "NUM_POLE", &
102 8 : i_val=numPole)
103 : CALL section_vals_val_get(pexsi_section, "IS_INERTIA_COUNT", &
104 8 : l_val=isInertiaCount)
105 : CALL section_vals_val_get(pexsi_section, "MAX_PEXSI_ITER", &
106 8 : i_val=maxPEXSIIter)
107 : CALL section_vals_val_get(pexsi_section, "MU_MIN_0", &
108 8 : r_val=muMin0)
109 : CALL section_vals_val_get(pexsi_section, "MU_MAX_0", &
110 8 : r_val=muMax0)
111 : CALL section_vals_val_get(pexsi_section, "MU_INERTIA_TOLERANCE", &
112 8 : r_val=muInertiaTolerance)
113 : CALL section_vals_val_get(pexsi_section, "MU_INERTIA_EXPANSION", &
114 8 : r_val=muInertiaExpansion)
115 : CALL section_vals_val_get(pexsi_section, "MU_PEXSI_SAFE_GUARD", &
116 8 : r_val=muPEXSISafeGuard)
117 : CALL section_vals_val_get(pexsi_section, "NUM_ELECTRON_INITIAL_TOLERANCE", &
118 8 : r_val=numElectronInitialTolerance)
119 : CALL section_vals_val_get(pexsi_section, "NUM_ELECTRON_PEXSI_TOLERANCE", &
120 8 : r_val=numElectronTargetTolerance)
121 : CALL section_vals_val_get(pexsi_section, "ORDERING", &
122 8 : i_val=ordering)
123 : CALL section_vals_val_get(pexsi_section, "ROW_ORDERING", &
124 8 : i_val=rowOrdering)
125 : CALL section_vals_val_get(pexsi_section, "NP_SYMB_FACT", &
126 8 : i_val=npSymbFact)
127 : CALL section_vals_val_get(pexsi_section, "VERBOSITY", &
128 8 : i_val=verbosity)
129 : CALL section_vals_val_get(pexsi_section, "MIN_RANKS_PER_POLE", &
130 8 : i_val=min_ranks_per_pole)
131 : CALL section_vals_val_get(pexsi_section, "CSR_SCREENING", &
132 8 : l_val=csr_screening)
133 :
134 8 : isInertiaCount_int = MERGE(1, 0, isInertiaCount) ! is integer in PEXSI
135 :
136 : ! Set default options inside PEXSI
137 8 : 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 8 : ordering=ordering, rowOrdering=rowOrdering, npSymbFact=npSymbFact, verbosity=verbosity)
145 :
146 8 : pexsi_env%num_ranks_per_pole = min_ranks_per_pole ! not a PEXSI option
147 8 : pexsi_env%csr_screening = csr_screening
148 :
149 8 : IF (numElectronInitialTolerance .LT. numElectronTargetTolerance) &
150 8 : numElectronInitialTolerance = numElectronTargetTolerance
151 :
152 8 : pexsi_env%tol_nel_initial = numElectronInitialTolerance
153 8 : pexsi_env%tol_nel_target = numElectronTargetTolerance
154 :
155 8 : 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 14 : 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 14 : CALL timeset(routineN, handle)
175 :
176 14 : logger => cp_get_default_logger()
177 14 : IF (logger%para_env%is_source()) THEN
178 7 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
179 : ELSE
180 7 : unit_nr = -1
181 : END IF
182 :
183 : ! Create template matrices fixing sparsity pattern for PEXSI
184 :
185 14 : IF (dbcsr_has_symmetry(template_matrix)) THEN
186 : CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, template_matrix, &
187 14 : "symmetric template matrix for CSR conversion")
188 : CALL dbcsr_desymmetrize(pexsi_env%dbcsr_template_matrix_sym, &
189 14 : 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 : template=pexsi_env%dbcsr_template_matrix_sym, &
199 14 : data_type=dbcsr_type_real_8)
200 14 : CALL dbcsr_copy(pexsi_env%csr_sparsity, pexsi_env%dbcsr_template_matrix_sym)
201 :
202 14 : CALL cp_dbcsr_to_csr_screening(ks_env, pexsi_env%csr_sparsity)
203 :
204 14 : 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 14 : numnodes=pexsi_env%num_ranks_per_pole)
210 :
211 14 : IF (unit_nr > 0) WRITE (unit_nr, "(/T2,A)") "SPARSITY OF THE OVERLAP MATRIX IN CSR FORMAT"
212 14 : CALL dbcsr_csr_print_sparsity(pexsi_env%csr_mat_s, unit_nr)
213 :
214 14 : CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_nonsym, pexsi_env%csr_mat_s)
215 :
216 14 : CALL dbcsr_csr_create(pexsi_env%csr_mat_ks, pexsi_env%csr_mat_s)
217 14 : CALL dbcsr_csr_create(pexsi_env%csr_mat_p, pexsi_env%csr_mat_s)
218 14 : CALL dbcsr_csr_create(pexsi_env%csr_mat_E, pexsi_env%csr_mat_s)
219 14 : CALL dbcsr_csr_create(pexsi_env%csr_mat_F, pexsi_env%csr_mat_s)
220 :
221 28 : 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 28 : template=template_matrix, matrix_type=dbcsr_type_no_symmetry)
225 : END DO
226 :
227 14 : CALL cp_pexsi_set_options(pexsi_env%options, numElectronPEXSITolerance=pexsi_env%tol_nel_initial)
228 :
229 14 : CALL timestop(handle)
230 :
231 14 : 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 14 : 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 14 : CALL timeset(routineN, handle)
249 :
250 14 : logger => cp_get_default_logger()
251 14 : IF (logger%para_env%is_source()) THEN
252 7 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
253 : ELSE
254 : unit_nr = -1
255 : END IF
256 :
257 28 : mu_total = SUM(mu_spin(1:pexsi_env%nspin))/REAL(pexsi_env%nspin, KIND=dp)
258 28 : kTS_total = SUM(pexsi_env%kTS)
259 14 : IF (pexsi_env%nspin .EQ. 1) kTS_total = kTS_total*2.0_dp
260 :
261 14 : IF (unit_nr > 0) THEN
262 : WRITE (unit_nr, "(/A,T55,F26.15)") &
263 7 : " PEXSI| Electronic entropic energy (a.u.):", kTS_total
264 : WRITE (unit_nr, "(A,T55,F26.15)") &
265 7 : " PEXSI| Chemical potential (a.u.):", mu_total
266 : END IF
267 :
268 14 : CALL dbcsr_release(pexsi_env%dbcsr_template_matrix_sym)
269 14 : CALL dbcsr_release(pexsi_env%dbcsr_template_matrix_nonsym)
270 14 : CALL dbcsr_release(pexsi_env%csr_sparsity)
271 14 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_p)
272 14 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_ks)
273 14 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_s)
274 14 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_E)
275 14 : CALL dbcsr_csr_destroy(pexsi_env%csr_mat_F)
276 28 : DO ispin = 1, pexsi_env%nspin
277 14 : CALL dbcsr_release(pexsi_env%max_ev_vector(ispin))
278 28 : CALL dbcsr_release(pexsi_env%matrix_w(ispin)%matrix)
279 : END DO
280 14 : CALL timestop(handle)
281 14 : pexsi_env%tol_nel_initial = pexsi_env%tol_nel_target ! Turn off adaptive threshold for subsequent SCF cycles
282 14 : 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 376 : 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 94 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: arnoldi_matrices
324 :
325 94 : CALL timeset(routineN, handle)
326 :
327 : ! get a useful output_unit
328 94 : logger => cp_get_default_logger()
329 94 : IF (logger%para_env%is_source()) THEN
330 47 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
331 : ELSE
332 : unit_nr = -1
333 : END IF
334 :
335 94 : 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 14 : IF (.NOT. dbcsr_has_symmetry(matrix_ks)) &
342 0 : CPABORT("PEXSI interface expects a non-symmetric DBCSR Kohn-Sham matrix")
343 14 : IF (.NOT. dbcsr_has_symmetry(matrix_s)) &
344 0 : 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 14 : .OR. (pexsi_env%csr_mat_ks%nzval_local%data_type .NE. dbcsr_type_real_8)) &
349 0 : 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 14 : IF (pexsi_env%csr_mat_s%nze_total .GE. INT(2, kind=int_8)**31) &
354 0 : CPABORT("Total number of non-zero elements of CSR matrix is too large to be handled by PEXSI")
355 : END IF
356 :
357 94 : CALL dbcsr_get_info(matrix_ks, distribution=dist)
358 94 : 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 94 : 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 94 : pexsi_env%csr_mat_s)
366 :
367 : CALL dbcsr_copy_into_existing(pexsi_env%dbcsr_template_matrix_sym, &
368 94 : matrix_ks)
369 : CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_sym, &
370 94 : pexsi_env%csr_mat_ks)
371 :
372 : ! Get PEXSI input delta_E (upper bound for largest eigenvalue) using Arnoldi
373 94 : NULLIFY (arnoldi_matrices)
374 94 : CALL dbcsr_allocate_matrix_set(arnoldi_matrices, 2)
375 94 : arnoldi_matrices(1)%matrix => matrix_ks
376 94 : 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 94 : generalized_ev=.TRUE., iram=.FALSE.)
380 94 : IF (iscf .GT. 1) CALL set_arnoldi_initial_vector(my_arnoldi, pexsi_env%max_ev_vector(ispin))
381 94 : CALL arnoldi_ev(arnoldi_matrices, my_arnoldi)
382 94 : 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 94 : 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 94 : pexsi_env%max_ev_vector(ispin))
387 94 : CALL deallocate_arnoldi_data(my_arnoldi)
388 94 : DEALLOCATE (arnoldi_matrices)
389 :
390 94 : nelectron_exact_pexsi = nelectron_exact
391 :
392 94 : CALL cp_pexsi_set_options(pexsi_env%options, deltaE=delta_E)
393 :
394 : ! Set PEXSI options appropriately for first SCF iteration
395 94 : IF (iscf .EQ. 1) THEN
396 : ! Get option isInertiaCount to reset it later on and set it to 1 for first SCF iteration
397 14 : CALL cp_pexsi_get_options(pexsi_env%options, isInertiaCount=isInertiaCount)
398 : CALL cp_pexsi_set_options(pexsi_env%options, isInertiaCount=1, &
399 14 : 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 94 : 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 141 : IF (unit_nr > 0) WRITE (unit_nr, '(/A,T41,L20)') " PEXSI| Do inertia counting?", &
412 94 : isInertiaCount_out .EQ. 1
413 94 : IF (unit_nr > 0) WRITE (unit_nr, '(A,T50,F5.2,T56,F5.2)') &
414 47 : " PEXSI| Guess for min mu, max mu", mu_min_in, mu_max_in
415 :
416 94 : IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,E20.3)') &
417 47 : " 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 94 : IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,F20.2)') &
423 47 : " PEXSI| Arnoldi est. spectral radius", delta_E
424 :
425 : ! Load data into PEXSI
426 : CALL cp_pexsi_load_real_hs_matrix( &
427 : pexsi_env%plan, &
428 : pexsi_env%options, &
429 : pexsi_env%csr_mat_ks%nrows_total, &
430 : INT(pexsi_env%csr_mat_ks%nze_total, kind=int_4), & ! TODO: update when PEXSI changes to Long Int
431 : pexsi_env%csr_mat_ks%nze_local, &
432 : pexsi_env%csr_mat_ks%nrows_local, &
433 : pexsi_env%csr_mat_ks%rowptr_local, &
434 : pexsi_env%csr_mat_ks%colind_local, &
435 : pexsi_env%csr_mat_ks%nzval_local%r_dp, &
436 : S_not_identity, &
437 94 : pexsi_env%csr_mat_s%nzval_local%r_dp)
438 :
439 : ! convert to spin restricted before passing number of electrons to PEXSI
440 : CALL convert_nspin_cp2k_pexsi(cp2k_to_pexsi, &
441 94 : 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 94 : n_total_inertia_iter, n_total_pexsi_iter)
447 :
448 : ! Check convergence
449 94 : nelectron_diff = nelectron_out - nelectron_exact_pexsi
450 94 : pexsi_convergence = ABS(nelectron_diff) .LT. nel_tol
451 :
452 94 : IF (unit_nr > 0) THEN
453 47 : IF (pexsi_convergence) THEN
454 47 : WRITE (unit_nr, '(/A)') " PEXSI| Converged"
455 : ELSE
456 0 : 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 47 : WRITE (unit_nr, '(A,T41,F20.6)') " PEXSI| Chemical potential", mu
463 :
464 47 : WRITE (unit_nr, '(A,T41,I20)') " PEXSI| PEXSI iterations", n_total_pexsi_iter
465 47 : WRITE (unit_nr, '(A,T41,I20/)') " PEXSI| Inertia counting iterations", &
466 94 : n_total_inertia_iter
467 : END IF
468 :
469 94 : IF (.NOT. pexsi_convergence) &
470 0 : CPABORT("PEXSI did not converge. Consider logPEXSI0 for more information.")
471 :
472 : ! Retrieve results from PEXSI
473 94 : IF (mynode < pexsi_env%mp_dims(1)*pexsi_env%mp_dims(2)) THEN
474 : CALL cp_pexsi_retrieve_real_dft_matrix( &
475 : pexsi_env%plan, &
476 : pexsi_env%csr_mat_p%nzval_local%r_dp, &
477 : pexsi_env%csr_mat_E%nzval_local%r_dp, &
478 : pexsi_env%csr_mat_F%nzval_local%r_dp, &
479 94 : energy_H, energy_S, free_energy)
480 : ! calculate entropic energy contribution -TS = A - U
481 94 : kTS = (free_energy - energy_H)
482 : END IF
483 :
484 : ! send kTS to all nodes:
485 94 : 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 94 : pexsi_env%csr_mat_p)
490 94 : 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 94 : pexsi_env%csr_mat_E)
493 94 : 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 94 : 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 94 : muMax0=mu_max_out)
502 :
503 : ! Reset isInertiaCount according to user input
504 94 : IF (iscf .EQ. 1) THEN
505 : CALL cp_pexsi_set_options(pexsi_env%options, isInertiaCount= &
506 14 : isInertiaCount)
507 : END IF
508 :
509 : ! Turn off symbolic factorization for subsequent calls
510 94 : IF (first_call) THEN
511 14 : CALL cp_pexsi_set_options(pexsi_env%options, isSymbolicFactorize=0)
512 : END IF
513 :
514 94 : CALL timestop(handle)
515 94 : 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 82 : 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 82 : CALL timeset(routineN, handle)
548 :
549 82 : tol_nel = pexsi_env%tol_nel_initial
550 :
551 82 : IF (initialize) THEN
552 : pexsi_env%adaptive_nel_alpha = &
553 8 : (pexsi_env%tol_nel_initial - pexsi_env%tol_nel_target)/(ABS(delta_scf) - eps_scf)
554 : pexsi_env%adaptive_nel_beta = &
555 8 : pexsi_env%tol_nel_initial - pexsi_env%adaptive_nel_alpha*ABS(delta_scf)
556 8 : pexsi_env%do_adaptive_tol_nel = .TRUE.
557 : END IF
558 82 : IF (pexsi_env%do_adaptive_tol_nel) THEN
559 74 : tol_nel = pexsi_env%adaptive_nel_alpha*ABS(delta_scf) + pexsi_env%adaptive_nel_beta
560 74 : tol_nel = MAX(tol_nel, pexsi_env%tol_nel_target)
561 74 : tol_nel = MIN(tol_nel, pexsi_env%tol_nel_initial)
562 : END IF
563 :
564 82 : check_convergence = (tol_nel .LE. pexsi_env%tol_nel_target)
565 :
566 82 : CALL cp_pexsi_set_options(pexsi_env%options, numElectronPEXSITolerance=tol_nel)
567 82 : CALL timestop(handle)
568 82 : 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 102 : 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 102 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w_qs
594 : TYPE(qs_energy_type), POINTER :: energy
595 :
596 102 : CALL timeset(routineN, handle)
597 :
598 102 : NULLIFY (energy)
599 :
600 : ! get a useful output_unit
601 102 : logger => cp_get_default_logger()
602 102 : IF (logger%para_env%is_source()) THEN
603 51 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
604 : ELSE
605 : unit_nr = -1
606 : END IF
607 :
608 102 : CALL get_qs_env(qs_env, energy=energy, matrix_w=matrix_w_qs)
609 :
610 102 : IF (PRESENT(matrix_w)) THEN
611 16 : DO ispin = 1, ls_scf_env%nspins
612 : CALL matrix_ls_to_qs(matrix_w_qs(ispin)%matrix, matrix_w(ispin)%matrix, &
613 8 : ls_scf_env%ls_mstruct, covariant=.FALSE.)
614 16 : 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 102 : IF (PRESENT(kTS)) THEN
619 188 : kTS_total = SUM(kTS)
620 94 : IF (ls_scf_env%nspins .EQ. 1) kTS_total = kTS_total*2.0_dp
621 94 : energy%kTS = kTS_total
622 : END IF
623 :
624 102 : CALL timestop(handle)
625 102 : END SUBROUTINE pexsi_to_qs
626 :
627 : END MODULE pexsi_methods
|
