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_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_desymmetrize, &
29 : dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_info, dbcsr_has_symmetry, &
30 : dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, &
31 : dbcsr_type_real_8
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 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_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 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 470 : 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(pexsi_env%dbcsr_template_matrix_sym, matrix_s, keep_sparsity=.TRUE.)
364 : CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_sym, &
365 94 : pexsi_env%csr_mat_s)
366 :
367 94 : CALL dbcsr_copy(pexsi_env%dbcsr_template_matrix_sym, matrix_ks, keep_sparsity=.TRUE.)
368 : CALL dbcsr_convert_dbcsr_to_csr(pexsi_env%dbcsr_template_matrix_sym, &
369 94 : pexsi_env%csr_mat_ks)
370 :
371 : ! Get PEXSI input delta_E (upper bound for largest eigenvalue) using Arnoldi
372 94 : NULLIFY (arnoldi_matrices)
373 94 : CALL dbcsr_allocate_matrix_set(arnoldi_matrices, 2)
374 94 : arnoldi_matrices(1)%matrix => matrix_ks
375 94 : arnoldi_matrices(2)%matrix => matrix_s
376 : CALL setup_arnoldi_data(my_arnoldi, arnoldi_matrices, max_iter=20, &
377 : threshold=1.0E-2_dp, selection_crit=2, nval_request=1, nrestarts=21, &
378 94 : generalized_ev=.TRUE., iram=.FALSE.)
379 94 : IF (iscf .GT. 1) CALL set_arnoldi_initial_vector(my_arnoldi, pexsi_env%max_ev_vector(ispin))
380 94 : CALL arnoldi_ev(arnoldi_matrices, my_arnoldi)
381 94 : delta_E = REAL(get_selected_ritz_val(my_arnoldi, 1), dp)
382 : ! increase delta_E a bit to make sure that it really is an upper bound
383 94 : delta_E = delta_E + 1.0E-2_dp*ABS(delta_E)
384 : CALL get_selected_ritz_vector(my_arnoldi, 1, arnoldi_matrices(1)%matrix, &
385 94 : pexsi_env%max_ev_vector(ispin))
386 94 : CALL deallocate_arnoldi_data(my_arnoldi)
387 94 : DEALLOCATE (arnoldi_matrices)
388 :
389 94 : nelectron_exact_pexsi = nelectron_exact
390 :
391 94 : CALL cp_pexsi_set_options(pexsi_env%options, deltaE=delta_E)
392 :
393 : ! Set PEXSI options appropriately for first SCF iteration
394 94 : IF (iscf .EQ. 1) THEN
395 : ! Get option isInertiaCount to reset it later on and set it to 1 for first SCF iteration
396 14 : CALL cp_pexsi_get_options(pexsi_env%options, isInertiaCount=isInertiaCount)
397 : CALL cp_pexsi_set_options(pexsi_env%options, isInertiaCount=1, &
398 14 : isSymbolicFactorize=1)
399 : END IF
400 :
401 : ! Write PEXSI options to output
402 : CALL cp_pexsi_get_options(pexsi_env%options, isInertiaCount=isInertiaCount_out, &
403 : isSymbolicFactorize=is_symbolic_factorize, &
404 : muMin0=mu_min_in, muMax0=mu_max_in, &
405 94 : NumElectronPEXSITolerance=nel_tol)
406 :
407 : ! IF(unit_nr>0) WRITE(unit_nr,'(/A,I4,A,I4)') " PEXSI| SCF", iscf, &
408 : ! ", spin component", ispin
409 :
410 141 : IF (unit_nr > 0) WRITE (unit_nr, '(/A,T41,L20)') " PEXSI| Do inertia counting?", &
411 94 : isInertiaCount_out .EQ. 1
412 : IF (unit_nr > 0) WRITE (unit_nr, '(A,T50,F5.2,T56,F5.2)') &
413 47 : " PEXSI| Guess for min mu, max mu", mu_min_in, mu_max_in
414 :
415 : IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,E20.3)') &
416 47 : " PEXSI| Tolerance in number of electrons", nel_tol
417 :
418 : ! IF(unit_nr>0) WRITE(unit_nr,'(A,T41,L20)') &
419 : ! " PEXSI| Do symbolic factorization?", is_symbolic_factorize.EQ.1
420 :
421 : IF (unit_nr > 0) WRITE (unit_nr, '(A,T41,F20.2)') &
422 47 : " PEXSI| Arnoldi est. spectral radius", delta_E
423 :
424 : ! Load data into PEXSI
425 : CALL cp_pexsi_load_real_hs_matrix( &
426 : pexsi_env%plan, &
427 : pexsi_env%options, &
428 : pexsi_env%csr_mat_ks%nrows_total, &
429 : INT(pexsi_env%csr_mat_ks%nze_total, kind=int_4), & ! TODO: update when PEXSI changes to Long Int
430 : pexsi_env%csr_mat_ks%nze_local, &
431 : pexsi_env%csr_mat_ks%nrows_local, &
432 : pexsi_env%csr_mat_ks%rowptr_local, &
433 : pexsi_env%csr_mat_ks%colind_local, &
434 : pexsi_env%csr_mat_ks%nzval_local%r_dp, &
435 : S_not_identity, &
436 94 : pexsi_env%csr_mat_s%nzval_local%r_dp)
437 :
438 : ! convert to spin restricted before passing number of electrons to PEXSI
439 : CALL convert_nspin_cp2k_pexsi(cp2k_to_pexsi, &
440 94 : numElectron=nelectron_exact_pexsi)
441 :
442 : ! Call DFT driver of PEXSI doing the actual calculation
443 : CALL cp_pexsi_dft_driver(pexsi_env%plan, pexsi_env%options, &
444 : nelectron_exact_pexsi, mu, nelectron_out, mu_min_out, mu_max_out, &
445 94 : n_total_inertia_iter, n_total_pexsi_iter)
446 :
447 : ! Check convergence
448 94 : nelectron_diff = nelectron_out - nelectron_exact_pexsi
449 94 : pexsi_convergence = ABS(nelectron_diff) .LT. nel_tol
450 :
451 94 : IF (unit_nr > 0) THEN
452 47 : IF (pexsi_convergence) THEN
453 47 : WRITE (unit_nr, '(/A)') " PEXSI| Converged"
454 : ELSE
455 0 : WRITE (unit_nr, '(/A)') " PEXSI| PEXSI did not converge!"
456 : END IF
457 :
458 : ! WRITE(unit_nr,'(A,T41,F20.10,T61,F20.10)') " PEXSI| Number of electrons", &
459 : ! nelectron_out/pexsi_env%nspin, nelectron_diff/pexsi_env%nspin
460 :
461 47 : WRITE (unit_nr, '(A,T41,F20.6)') " PEXSI| Chemical potential", mu
462 :
463 47 : WRITE (unit_nr, '(A,T41,I20)') " PEXSI| PEXSI iterations", n_total_pexsi_iter
464 47 : WRITE (unit_nr, '(A,T41,I20/)') " PEXSI| Inertia counting iterations", &
465 94 : n_total_inertia_iter
466 : END IF
467 :
468 94 : IF (.NOT. pexsi_convergence) &
469 0 : CPABORT("PEXSI did not converge. Consider logPEXSI0 for more information.")
470 :
471 : ! Retrieve results from PEXSI
472 94 : IF (mynode < pexsi_env%mp_dims(1)*pexsi_env%mp_dims(2)) THEN
473 : CALL cp_pexsi_retrieve_real_dft_matrix( &
474 : pexsi_env%plan, &
475 : pexsi_env%csr_mat_p%nzval_local%r_dp, &
476 : pexsi_env%csr_mat_E%nzval_local%r_dp, &
477 : pexsi_env%csr_mat_F%nzval_local%r_dp, &
478 94 : energy_H, energy_S, free_energy)
479 : ! calculate entropic energy contribution -TS = A - U
480 94 : kTS = (free_energy - energy_H)
481 : END IF
482 :
483 : ! send kTS to all nodes:
484 94 : CALL pexsi_env%mp_group%bcast(kTS, 0)
485 :
486 : ! Convert PEXSI CSR matrices to DBCSR matrices
487 : CALL dbcsr_convert_csr_to_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
488 94 : pexsi_env%csr_mat_p)
489 94 : CALL dbcsr_copy(matrix_p, pexsi_env%dbcsr_template_matrix_nonsym)
490 : CALL dbcsr_convert_csr_to_dbcsr(pexsi_env%dbcsr_template_matrix_nonsym, &
491 94 : pexsi_env%csr_mat_E)
492 94 : CALL dbcsr_copy(matrix_w%matrix, pexsi_env%dbcsr_template_matrix_nonsym)
493 :
494 : ! Convert to spin unrestricted
495 : CALL convert_nspin_cp2k_pexsi(pexsi_to_cp2k, matrix_p=matrix_p, &
496 94 : matrix_w=matrix_w, kTS=kTS)
497 :
498 : ! Pass resulting mu as initial guess for next SCF to PEXSI
499 : CALL cp_pexsi_set_options(pexsi_env%options, mu0=mu, muMin0=mu_min_out, &
500 94 : muMax0=mu_max_out)
501 :
502 : ! Reset isInertiaCount according to user input
503 94 : IF (iscf .EQ. 1) THEN
504 : CALL cp_pexsi_set_options(pexsi_env%options, isInertiaCount= &
505 14 : isInertiaCount)
506 : END IF
507 :
508 : ! Turn off symbolic factorization for subsequent calls
509 94 : IF (first_call) THEN
510 14 : CALL cp_pexsi_set_options(pexsi_env%options, isSymbolicFactorize=0)
511 : END IF
512 :
513 94 : CALL timestop(handle)
514 94 : END SUBROUTINE density_matrix_pexsi
515 :
516 : ! **************************************************************************************************
517 : !> \brief Set PEXSI convergence tolerance (numElectronPEXSITolerance), adapted
518 : !> to the convergence error of the previous SCF step. The tolerance is
519 : !> calculated using an initial convergence threshold (tol_nel_initial)
520 : !> and a target convergence threshold (tol_nel_target):
521 : !> numElectronPEXSITolerance(delta_scf) = alpha*delta_scf+beta
522 : !> where alpha and beta are chosen such that
523 : !> numElectronPEXSITolerance(delta_scf_0) = tol_nel_initial
524 : !> numElectronPEXSITolerance(eps_scf) = tol_nel_target
525 : !> and delta_scf_0 is the initial SCF convergence error.
526 : !> \param pexsi_env ...
527 : !> \param delta_scf convergence error of previous SCF step
528 : !> \param eps_scf SCF convergence criterion
529 : !> \param initialize whether or not adaptive thresholing should be initialized
530 : !> with delta_scf as initial convergence error
531 : !> \param check_convergence is set to .FALSE. if convergence in number of electrons
532 : !> will not be achieved in next SCF step
533 : ! **************************************************************************************************
534 82 : SUBROUTINE pexsi_set_convergence_tolerance(pexsi_env, delta_scf, eps_scf, initialize, &
535 : check_convergence)
536 : TYPE(lib_pexsi_env), INTENT(INOUT) :: pexsi_env
537 : REAL(KIND=dp), INTENT(IN) :: delta_scf, eps_scf
538 : LOGICAL, INTENT(IN) :: initialize
539 : LOGICAL, INTENT(OUT) :: check_convergence
540 :
541 : CHARACTER(len=*), PARAMETER :: routineN = 'pexsi_set_convergence_tolerance'
542 :
543 : INTEGER :: handle
544 : REAL(KIND=dp) :: tol_nel
545 :
546 82 : CALL timeset(routineN, handle)
547 :
548 82 : tol_nel = pexsi_env%tol_nel_initial
549 :
550 82 : IF (initialize) THEN
551 : pexsi_env%adaptive_nel_alpha = &
552 8 : (pexsi_env%tol_nel_initial - pexsi_env%tol_nel_target)/(ABS(delta_scf) - eps_scf)
553 : pexsi_env%adaptive_nel_beta = &
554 8 : pexsi_env%tol_nel_initial - pexsi_env%adaptive_nel_alpha*ABS(delta_scf)
555 8 : pexsi_env%do_adaptive_tol_nel = .TRUE.
556 : END IF
557 82 : IF (pexsi_env%do_adaptive_tol_nel) THEN
558 74 : tol_nel = pexsi_env%adaptive_nel_alpha*ABS(delta_scf) + pexsi_env%adaptive_nel_beta
559 74 : tol_nel = MAX(tol_nel, pexsi_env%tol_nel_target)
560 74 : tol_nel = MIN(tol_nel, pexsi_env%tol_nel_initial)
561 : END IF
562 :
563 82 : check_convergence = (tol_nel .LE. pexsi_env%tol_nel_target)
564 :
565 82 : CALL cp_pexsi_set_options(pexsi_env%options, numElectronPEXSITolerance=tol_nel)
566 82 : CALL timestop(handle)
567 82 : END SUBROUTINE
568 :
569 : ! **************************************************************************************************
570 : !> \brief Pass energy weighted density matrix and entropic energy contribution
571 : !> to QS environment
572 : !> \param ls_scf_env ...
573 : !> \param qs_env ...
574 : !> \param kTS ...
575 : !> \param matrix_w ...
576 : !> \par History
577 : !> 12.2014 created [Patrick Seewald]
578 : !> \author Patrick Seewald
579 : ! **************************************************************************************************
580 102 : SUBROUTINE pexsi_to_qs(ls_scf_env, qs_env, kTS, matrix_w)
581 : TYPE(ls_scf_env_type) :: ls_scf_env
582 : TYPE(qs_environment_type), INTENT(INOUT), POINTER :: qs_env
583 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: kTS
584 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
585 : OPTIONAL :: matrix_w
586 :
587 : CHARACTER(len=*), PARAMETER :: routineN = 'pexsi_to_qs'
588 :
589 : INTEGER :: handle, ispin, unit_nr
590 : REAL(KIND=dp) :: kTS_total
591 : TYPE(cp_logger_type), POINTER :: logger
592 102 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w_qs
593 : TYPE(qs_energy_type), POINTER :: energy
594 :
595 102 : CALL timeset(routineN, handle)
596 :
597 102 : NULLIFY (energy)
598 :
599 : ! get a useful output_unit
600 102 : logger => cp_get_default_logger()
601 102 : IF (logger%para_env%is_source()) THEN
602 51 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
603 : ELSE
604 : unit_nr = -1
605 : END IF
606 :
607 102 : CALL get_qs_env(qs_env, energy=energy, matrix_w=matrix_w_qs)
608 :
609 102 : IF (PRESENT(matrix_w)) THEN
610 16 : DO ispin = 1, ls_scf_env%nspins
611 : CALL matrix_ls_to_qs(matrix_w_qs(ispin)%matrix, matrix_w(ispin)%matrix, &
612 8 : ls_scf_env%ls_mstruct, covariant=.FALSE.)
613 16 : IF (ls_scf_env%nspins .EQ. 1) CALL dbcsr_scale(matrix_w_qs(ispin)%matrix, 2.0_dp)
614 : END DO
615 : END IF
616 :
617 102 : IF (PRESENT(kTS)) THEN
618 188 : kTS_total = SUM(kTS)
619 94 : IF (ls_scf_env%nspins .EQ. 1) kTS_total = kTS_total*2.0_dp
620 94 : energy%kTS = kTS_total
621 : END IF
622 :
623 102 : CALL timestop(handle)
624 102 : END SUBROUTINE pexsi_to_qs
625 :
626 : END MODULE pexsi_methods
|
