(git:ba1d7ca)
Loading...
Searching...
No Matches
qs_density_fit.F
Go to the documentation of this file.
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 Constrained fitting of a real-space density by an AO density matrix.
10!> This is a dense reference implementation with cubic spectral projections.
11!> \author CP2K developers group
12! **************************************************************************************************
13MODULE qs_density_fit
14 USE cp_blacs_env, ONLY: cp_blacs_env_type
15 USE cp_control_types, ONLY: dft_control_type
16 USE cp_dbcsr_api, ONLY: dbcsr_copy,&
17 dbcsr_create,&
18 dbcsr_get_info,&
19 dbcsr_p_type,&
20 dbcsr_release,&
21 dbcsr_set,&
22 dbcsr_type,&
23 dbcsr_type_symmetric
24 USE cp_dbcsr_contrib, ONLY: dbcsr_dot
25 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
26 copy_fm_to_dbcsr
27 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
28 cp_fm_scale_and_add,&
29 cp_fm_trace,&
30 cp_fm_uplo_to_full
31 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
32 cp_fm_cholesky_restore
33 USE cp_fm_diag, ONLY: cp_fm_syevd
34 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
35 cp_fm_struct_release,&
36 cp_fm_struct_type
37 USE cp_fm_types, ONLY: cp_fm_create,&
38 cp_fm_release,&
39 cp_fm_to_fm,&
40 cp_fm_type
41 USE cp_log_handling, ONLY: cp_get_default_logger,&
42 cp_logger_get_default_unit_nr,&
43 cp_logger_type
44 USE kinds, ONLY: dp
45 USE message_passing, ONLY: mp_para_env_type
46 USE parallel_gemm_api, ONLY: parallel_gemm
47 USE pw_env_types, ONLY: pw_env_get,&
48 pw_env_type
49 USE pw_methods, ONLY: pw_axpy,&
50 pw_copy,&
51 pw_integral_ab,&
52 pw_integrate_function
53 USE pw_pool_types, ONLY: pw_pool_type
54 USE pw_types, ONLY: pw_c1d_gs_type,&
55 pw_r3d_rs_type
56 USE qs_collocate_density, ONLY: calculate_rho_elec
57 USE qs_environment_types, ONLY: get_qs_env,&
58 qs_environment_type
59 USE qs_integrate_potential, ONLY: integrate_v_rspace
60 USE qs_ks_types, ONLY: qs_ks_env_type
61 USE qs_rho_types, ONLY: qs_rho_get,&
62 qs_rho_type
63#include "./base/base_uses.f90"
64
65 IMPLICIT NONE
66 PRIVATE
67
68 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_density_fit'
69
70 PUBLIC :: fit_constrained_density, fit_relative_entropy_density
71
72CONTAINS
73
74! **************************************************************************************************
75!> \brief Fits the density currently stored in rho_struct with a spin-restricted AO density matrix.
76!> The density matrix is projected in an orthonormal AO representation onto
77!> 0 <= P <= 2 and Tr(P) = Ne after every gradient step.
78!> \param qs_env QS environment
79!> \param rho_struct density structure containing the target density on entry and the fit on exit
80!> \param max_iter maximum number of projected-gradient iterations
81!> \param eps_rms requested RMS density residual
82!> \param step_size initial spectral projected-gradient step
83!> \param max_backtrack maximum number of step halvings per iteration
84! **************************************************************************************************
85 SUBROUTINE fit_constrained_density(qs_env, rho_struct, max_iter, eps_rms, step_size, max_backtrack)
86 TYPE(qs_environment_type), POINTER :: qs_env
87 TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
88 INTEGER, INTENT(IN) :: max_iter
89 REAL(kind=dp), INTENT(IN) :: eps_rms, step_size
90 INTEGER, INTENT(IN) :: max_backtrack
91
92 CALL fit_density_matrix(qs_env, rho_struct, max_iter, eps_rms, step_size, max_backtrack)
93
94 END SUBROUTINE fit_constrained_density
95
96! **************************************************************************************************
97!> \brief Reconstructs an AO density matrix using the Fermi matrix of H[n_cube] as a
98!> fermionic relative-entropy prior.
99!> \param qs_env QS environment
100!> \param rho_struct density structure containing the target density on entry and the fit on exit
101!> \param prior_hamiltonian frozen Kohn-Sham matrix built from the target cube density
102!> \param temperature electronic temperature of the Fermi prior
103!> \param entropy_weight weight of the dimensionless fermionic relative entropy
104!> \param max_iter maximum number of projected-gradient iterations
105!> \param eps_rms requested RMS density residual
106!> \param step_size initial spectral projected-gradient step
107!> \param max_backtrack maximum number of step halvings per iteration
108! **************************************************************************************************
109 SUBROUTINE fit_relative_entropy_density(qs_env, rho_struct, prior_hamiltonian, temperature, &
110 entropy_weight, max_iter, eps_rms, step_size, max_backtrack)
111 TYPE(qs_environment_type), POINTER :: qs_env
112 TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
113 TYPE(dbcsr_type), INTENT(IN) :: prior_hamiltonian
114 REAL(kind=dp), INTENT(IN) :: temperature, entropy_weight
115 INTEGER, INTENT(IN) :: max_iter
116 REAL(kind=dp), INTENT(IN) :: eps_rms, step_size
117 INTEGER, INTENT(IN) :: max_backtrack
118
119 CALL fit_density_matrix(qs_env, rho_struct, max_iter, eps_rms, step_size, max_backtrack, &
120 prior_hamiltonian, temperature, entropy_weight)
121
122 END SUBROUTINE fit_relative_entropy_density
123
124! **************************************************************************************************
125!> \brief Dense implementation shared by the least-squares and relative-entropy fits.
126!> \param qs_env ...
127!> \param rho_struct ...
128!> \param max_iter ...
129!> \param eps_rms ...
130!> \param step_size ...
131!> \param max_backtrack ...
132!> \param prior_hamiltonian ...
133!> \param temperature ...
134!> \param entropy_weight ...
135! **************************************************************************************************
136 SUBROUTINE fit_density_matrix(qs_env, rho_struct, max_iter, eps_rms, step_size, max_backtrack, &
137 prior_hamiltonian, temperature, entropy_weight)
138 TYPE(qs_environment_type), POINTER :: qs_env
139 TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
140 INTEGER, INTENT(IN) :: max_iter
141 REAL(kind=dp), INTENT(IN) :: eps_rms, step_size
142 INTEGER, INTENT(IN) :: max_backtrack
143 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: prior_hamiltonian
144 REAL(kind=dp), INTENT(IN), OPTIONAL :: temperature, entropy_weight
145
146 CHARACTER(LEN=*), PARAMETER :: routinen = 'fit_density_matrix'
147
148 INTEGER :: backtrack, handle, i, info, iteration, &
149 nao, nelectron, unit_nr
150 LOGICAL :: accepted, converged, do_kpoints, &
151 have_previous, use_relative_entropy
152 REAL(kind=dp) :: alpha, chemical_potential, commutator_norm, direction_derivative, &
153 direction_norm2, fitted_grid_trace, idempotency_error, matrix_trace, objective, &
154 objective_trial, prior_fraction, prior_log_one_minus, relative_entropy, &
155 relative_entropy_trial, relative_temperature, relative_weight, rms, rms_trial, ss, &
156 step_length, sy
157 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, prior_eigenvalues, &
158 prior_logit_values, prior_occupations
159 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
160 TYPE(cp_blacs_env_type), POINTER :: blacs_env
161 TYPE(cp_fm_struct_type), POINTER :: fm_struct
162 TYPE(cp_fm_type) :: density_orth, direction_orth, eigenvectors, gradient_orth, &
163 hamiltonian_orth, overlap_chol, previous_density, previous_gradient, prior_logit, &
164 trial_orth, work1, work2
165 TYPE(cp_logger_type), POINTER :: logger
166 TYPE(dbcsr_p_type) :: gradient_ao
167 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
168 TYPE(dft_control_type), POINTER :: dft_control
169 TYPE(mp_para_env_type), POINTER :: para_env
170 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
171 TYPE(pw_env_type), POINTER :: pw_env
172 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
173 TYPE(pw_r3d_rs_type) :: residual_rspace, target_rspace
174 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
175 TYPE(qs_ks_env_type), POINTER :: ks_env
176
177 CALL timeset(routinen, handle)
178
179 use_relative_entropy = PRESENT(prior_hamiltonian)
180 cpassert(use_relative_entropy .EQV. PRESENT(temperature))
181 cpassert(use_relative_entropy .EQV. PRESENT(entropy_weight))
182 IF (use_relative_entropy) THEN
183 relative_temperature = temperature
184 relative_weight = entropy_weight
185 cpassert(relative_temperature > 0.0_dp)
186 cpassert(relative_weight >= 0.0_dp)
187 ELSE
188 relative_temperature = 1.0_dp
189 relative_weight = 0.0_dp
190 END IF
191
192 NULLIFY (auxbas_pw_pool, blacs_env, dft_control, fm_struct, gradient_ao%matrix, &
193 ks_env, matrix_s, para_env, pw_env, rho_ao, rho_g, rho_r, tot_rho_r)
194
195 CALL get_qs_env(qs_env, blacs_env=blacs_env, dft_control=dft_control, &
196 do_kpoints=do_kpoints, ks_env=ks_env, matrix_s=matrix_s, &
197 nelectron_total=nelectron, para_env=para_env, pw_env=pw_env)
198
199 IF (do_kpoints) cpabort("Harris CUBE_FIT is currently available only at the Gamma point")
200 IF (dft_control%nspins /= 1) THEN
201 cpabort("Harris CUBE_FIT currently requires a spin-restricted calculation")
202 END IF
203 IF (dft_control%qs_control%gapw) cpabort("Harris CUBE_FIT currently supports GPW only")
204
205 CALL qs_rho_get(rho_struct, rho_ao=rho_ao, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
206 cpassert(ASSOCIATED(rho_ao) .AND. SIZE(rho_ao) == 1)
207 cpassert(ASSOCIATED(rho_r) .AND. SIZE(rho_r) == 1)
208 cpassert(ASSOCIATED(rho_g) .AND. SIZE(rho_g) == 1)
209
210 CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
211 IF (nelectron < 0 .OR. nelectron > 2*nao) THEN
212 cpabort("Electron count is incompatible with the constrained AO density-matrix bounds")
213 END IF
214
215 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
216 CALL auxbas_pw_pool%create_pw(target_rspace)
217 CALL auxbas_pw_pool%create_pw(residual_rspace)
218 CALL pw_copy(rho_r(1), target_rspace)
219
220 ALLOCATE (gradient_ao%matrix)
221 CALL dbcsr_create(gradient_ao%matrix, template=matrix_s(1)%matrix, &
222 name="Harris density-fit gradient", matrix_type=dbcsr_type_symmetric)
223 CALL dbcsr_copy(gradient_ao%matrix, matrix_s(1)%matrix)
224 CALL dbcsr_set(gradient_ao%matrix, 0.0_dp)
225
226 CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
227 nrow_global=nao, ncol_global=nao)
228 CALL cp_fm_create(overlap_chol, fm_struct, name="density-fit overlap")
229 CALL cp_fm_create(density_orth, fm_struct, name="density fit in orthogonal basis")
230 CALL cp_fm_create(direction_orth, fm_struct, name="density-fit search direction")
231 CALL cp_fm_create(gradient_orth, fm_struct, name="density-fit gradient in orthogonal basis")
232 CALL cp_fm_create(previous_density, fm_struct, name="previous density-fit matrix")
233 CALL cp_fm_create(previous_gradient, fm_struct, name="previous density-fit gradient")
234 CALL cp_fm_create(trial_orth, fm_struct, name="density-fit trial matrix")
235 CALL cp_fm_create(eigenvectors, fm_struct, name="density-fit eigenvectors")
236 IF (use_relative_entropy) THEN
237 CALL cp_fm_create(hamiltonian_orth, fm_struct, name="relative-entropy prior Hamiltonian")
238 CALL cp_fm_create(prior_logit, fm_struct, name="relative-entropy prior logit")
239 END IF
240 CALL cp_fm_create(work1, fm_struct, name="density-fit work matrix 1")
241 CALL cp_fm_create(work2, fm_struct, name="density-fit work matrix 2")
242 CALL cp_fm_struct_release(fm_struct)
243 ALLOCATE (eigenvalues(nao))
244 IF (use_relative_entropy) THEN
245 ALLOCATE (prior_eigenvalues(nao), prior_logit_values(nao), prior_occupations(nao))
246 END IF
247
248 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, overlap_chol)
249 CALL cp_fm_uplo_to_full(overlap_chol, work1)
250 CALL cp_fm_cholesky_decompose(overlap_chol, info_out=info)
251 IF (info /= 0) cpabort("Overlap Cholesky decomposition failed in Harris CUBE_FIT")
252
253 ! X = U P U^T for S = U^T U. Tr(X) is the electron count.
254 IF (use_relative_entropy) THEN
255 ! H_tilde = U^-T H[n_cube] U^-1 in the same orthonormal representation.
256 CALL copy_dbcsr_to_fm(prior_hamiltonian, work1)
257 CALL cp_fm_uplo_to_full(work1, work2)
258 CALL cp_fm_cholesky_restore(work1, nao, overlap_chol, work2, &
259 "SOLVE", pos="LEFT", transa="T")
260 CALL cp_fm_cholesky_restore(work2, nao, overlap_chol, hamiltonian_orth, &
261 "SOLVE", pos="RIGHT")
262 CALL cp_fm_to_fm(hamiltonian_orth, work1)
263 CALL cp_fm_syevd(work1, eigenvectors, prior_eigenvalues)
264 CALL fixed_trace_fermi_occupations(prior_eigenvalues, real(nelectron, dp), &
265 relative_temperature, chemical_potential, prior_occupations)
266 CALL matrix_from_eigensystem(eigenvectors, prior_occupations, density_orth, work1)
267 prior_log_one_minus = 0.0_dp
268 DO i = 1, nao
269 prior_fraction = clipped_occupation_fraction(prior_occupations(i))
270 prior_logit_values(i) = log(prior_fraction/(1.0_dp - prior_fraction))
271 prior_log_one_minus = prior_log_one_minus + log(1.0_dp - prior_fraction)
272 END DO
273 CALL matrix_from_eigensystem(eigenvectors, prior_logit_values, prior_logit, work1)
274 ELSE
275 CALL copy_dbcsr_to_fm(rho_ao(1)%matrix, work1)
276 CALL cp_fm_uplo_to_full(work1, work2)
277 CALL cp_fm_cholesky_restore(work1, nao, overlap_chol, work2, "MULTIPLY", pos="LEFT")
278 CALL cp_fm_cholesky_restore(work2, nao, overlap_chol, density_orth, &
279 "MULTIPLY", pos="RIGHT", transa="T")
280 CALL project_density_matrix(density_orth, eigenvectors, work1, eigenvalues, &
281 REAL(nelectron, dp), 2.0_dp)
282 chemical_potential = 0.0_dp
283 prior_log_one_minus = 0.0_dp
284 END IF
285 CALL orthogonal_to_ao(density_orth, overlap_chol, work1, work2, nao)
286 CALL copy_fm_to_dbcsr(work2, rho_ao(1)%matrix)
287 CALL evaluate_density_fit(qs_env, ks_env, rho_ao, rho_r, rho_g, tot_rho_r, &
288 target_rspace, residual_rspace, objective, rms)
289 relative_entropy = 0.0_dp
290 IF (use_relative_entropy) THEN
291 CALL fermionic_relative_entropy(density_orth, prior_logit, eigenvectors, work1, &
292 eigenvalues, prior_log_one_minus, relative_entropy)
293 objective = objective + relative_weight*relative_entropy
294 END IF
295
296 logger => cp_get_default_logger()
297 IF (logger%para_env%is_source()) THEN
298 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
299 ELSE
300 unit_nr = -1
301 END IF
302 IF (unit_nr > 0) THEN
303 IF (use_relative_entropy) THEN
304 WRITE (unit_nr, "(/,T3,A,I0)") &
305 "HARRIS| Relative-entropy AO density reconstruction; basis functions: ", nao
306 WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Prior electronic temperature [a.u.]: ", &
307 relative_temperature
308 WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Prior chemical potential [a.u.]: ", &
309 chemical_potential
310 WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Relative-entropy weight: ", relative_weight
311 WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Initial fermionic relative entropy: ", &
312 relative_entropy
313 ELSE
314 WRITE (unit_nr, "(/,T3,A,I0)") "HARRIS| Constrained AO density fit; basis functions: ", nao
315 END IF
316 WRITE (unit_nr, "(T3,A,I0)") "HARRIS| Constrained electron count: ", nelectron
317 WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Initial density RMS error: ", rms
318 IF (use_relative_entropy) THEN
319 WRITE (unit_nr, "(T3,A)") &
320 "HARRIS| Iteration RMS error Relative entropy Objective Step"
321 ELSE
322 WRITE (unit_nr, "(T3,A)") "HARRIS| Iteration RMS error Objective Step"
323 END IF
324 END IF
325
326 converged = rms <= eps_rms
327 alpha = step_size
328 have_previous = .false.
329 iteration = 0
330 DO WHILE (iteration < max_iter .AND. .NOT. converged)
331 iteration = iteration + 1
332
333 CALL dbcsr_set(gradient_ao%matrix, 0.0_dp)
334 CALL integrate_v_rspace(v_rspace=residual_rspace, hmat=gradient_ao, &
335 qs_env=qs_env, calculate_forces=.false.)
336 CALL copy_dbcsr_to_fm(gradient_ao%matrix, work1)
337 CALL cp_fm_uplo_to_full(work1, work2)
338
339 ! dF/dX = U^-T (dF/dP) U^-1.
340 CALL cp_fm_cholesky_restore(work1, nao, overlap_chol, work2, &
341 "SOLVE", pos="LEFT", transa="T")
342 CALL cp_fm_cholesky_restore(work2, nao, overlap_chol, gradient_orth, &
343 "SOLVE", pos="RIGHT")
344
345 IF (use_relative_entropy .AND. relative_weight > 0.0_dp) THEN
346 CALL add_relative_entropy_gradient(density_orth, prior_logit, eigenvectors, work1, &
347 work2, eigenvalues, relative_weight, gradient_orth)
348 END IF
349
350 ! Barzilai-Borwein spectral step for the convex projected-gradient iteration.
351 IF (have_previous) THEN
352 CALL cp_fm_to_fm(density_orth, work1)
353 CALL cp_fm_scale_and_add(1.0_dp, work1, beta=-1.0_dp, matrix_b=previous_density)
354 CALL cp_fm_to_fm(gradient_orth, work2)
355 CALL cp_fm_scale_and_add(1.0_dp, work2, beta=-1.0_dp, matrix_b=previous_gradient)
356 CALL cp_fm_trace(work1, work1, ss)
357 CALL cp_fm_trace(work1, work2, sy)
358 IF (sy > 100.0_dp*epsilon(sy)*max(1.0_dp, ss)) THEN
359 alpha = min(100.0_dp*step_size, max(1.0e-6_dp*step_size, ss/sy))
360 END IF
361 END IF
362 CALL cp_fm_to_fm(density_orth, previous_density)
363 CALL cp_fm_to_fm(gradient_orth, previous_gradient)
364 have_previous = .true.
365
366 ! Project once. Backtracking then stays on the feasible line segment between
367 ! the accepted matrix and its projected spectral-gradient point.
368 CALL cp_fm_to_fm(density_orth, trial_orth)
369 CALL cp_fm_scale_and_add(1.0_dp, trial_orth, beta=-alpha, matrix_b=gradient_orth)
370 CALL project_density_matrix(trial_orth, eigenvectors, work1, eigenvalues, &
371 REAL(nelectron, dp), 2.0_dp)
372 CALL cp_fm_to_fm(trial_orth, direction_orth)
373 CALL cp_fm_scale_and_add(1.0_dp, direction_orth, beta=-1.0_dp, matrix_b=density_orth)
374 CALL cp_fm_trace(gradient_orth, direction_orth, direction_derivative)
375 CALL cp_fm_trace(direction_orth, direction_orth, direction_norm2)
376 IF (direction_norm2 <= 100.0_dp*epsilon(direction_norm2)) THEN
377 IF (unit_nr > 0) WRITE (unit_nr, "(T3,A)") &
378 "HARRIS| Density fit stopped: projected gradient reached a stationary point"
379 converged = .true.
380 EXIT
381 END IF
382
383 accepted = .false.
384 step_length = 1.0_dp
385 DO backtrack = 0, max_backtrack
386 CALL cp_fm_to_fm(density_orth, trial_orth)
387 CALL cp_fm_scale_and_add(1.0_dp, trial_orth, beta=step_length, matrix_b=direction_orth)
388 CALL orthogonal_to_ao(trial_orth, overlap_chol, work1, work2, nao)
389 CALL copy_fm_to_dbcsr(work2, rho_ao(1)%matrix)
390 CALL evaluate_density_fit(qs_env, ks_env, rho_ao, rho_r, rho_g, tot_rho_r, &
391 target_rspace, residual_rspace, objective_trial, rms_trial)
392 relative_entropy_trial = 0.0_dp
393 IF (use_relative_entropy) THEN
394 CALL fermionic_relative_entropy(trial_orth, prior_logit, eigenvectors, work1, &
395 eigenvalues, prior_log_one_minus, relative_entropy_trial)
396 objective_trial = objective_trial + relative_weight*relative_entropy_trial
397 END IF
398 IF (objective_trial <= objective + 1.0e-4_dp*step_length*direction_derivative) THEN
399 accepted = .true.
400 EXIT
401 END IF
402 step_length = 0.5_dp*step_length
403 END DO
404
405 IF (.NOT. accepted) THEN
406 ! Leave both the AO matrix and the grids at the last accepted point.
407 CALL orthogonal_to_ao(density_orth, overlap_chol, work1, work2, nao)
408 CALL copy_fm_to_dbcsr(work2, rho_ao(1)%matrix)
409 CALL evaluate_density_fit(qs_env, ks_env, rho_ao, rho_r, rho_g, tot_rho_r, &
410 target_rspace, residual_rspace, objective, rms)
411 IF (use_relative_entropy) THEN
412 CALL fermionic_relative_entropy(density_orth, prior_logit, eigenvectors, work1, &
413 eigenvalues, prior_log_one_minus, relative_entropy)
414 objective = objective + relative_weight*relative_entropy
415 END IF
416 IF (unit_nr > 0) WRITE (unit_nr, "(T3,A)") &
417 "HARRIS| Density fit stopped: projected line search reached a stationary point"
418 converged = .true.
419 EXIT
420 END IF
421
422 CALL cp_fm_to_fm(trial_orth, density_orth)
423 objective = objective_trial
424 rms = rms_trial
425 relative_entropy = relative_entropy_trial
426 IF (unit_nr > 0) THEN
427 IF (use_relative_entropy) THEN
428 WRITE (unit_nr, "(T3,A,I7,4ES19.9)") &
429 "HARRIS| ", iteration, rms, relative_entropy, objective, alpha*step_length
430 ELSE
431 WRITE (unit_nr, "(T3,A,I7,3ES19.9)") &
432 "HARRIS| ", iteration, rms, objective, alpha*step_length
433 END IF
434 END IF
435 converged = rms <= eps_rms
436 END DO
437
438 CALL dbcsr_dot(rho_ao(1)%matrix, matrix_s(1)%matrix, matrix_trace)
439 ! pw_integrate_function performs an MPI reduction and therefore has to be
440 ! called collectively, not only by the rank that owns the output unit.
441 fitted_grid_trace = pw_integrate_function(rho_r(1), isign=1)
442 commutator_norm = 0.0_dp
443 idempotency_error = 0.0_dp
444 IF (use_relative_entropy) THEN
445 CALL matrix_commutator_norm(density_orth, hamiltonian_orth, work1, work2, nao, &
446 commutator_norm)
447 CALL density_matrix_idempotency_error(density_orth, work1, work2, nao, idempotency_error)
448 END IF
449 IF (unit_nr > 0) THEN
450 IF (rms <= eps_rms) THEN
451 WRITE (unit_nr, "(T3,A,I0,A)") "HARRIS| Density fit reached its RMS target after ", iteration, " iterations"
452 ELSE IF (iteration >= max_iter) THEN
453 WRITE (unit_nr, "(T3,A,I0,A)") "HARRIS| Density fit stopped after ", iteration, " iterations"
454 END IF
455 WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Final density RMS error: ", rms
456 WRITE (unit_nr, "(T3,A,F20.10)") "HARRIS| AO density-matrix electron count: ", matrix_trace
457 WRITE (unit_nr, "(T3,A,F20.10)") "HARRIS| Fitted grid electron count: ", fitted_grid_trace
458 IF (use_relative_entropy) THEN
459 WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Final fermionic relative entropy: ", &
460 relative_entropy
461 WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Prior-Hamiltonian commutator norm: ", &
462 commutator_norm
463 WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Occupation idempotency error: ", &
464 idempotency_error
465 END IF
466 END IF
467
468 DEALLOCATE (eigenvalues)
469 IF (use_relative_entropy) THEN
470 DEALLOCATE (prior_eigenvalues, prior_logit_values, prior_occupations)
471 END IF
472 CALL cp_fm_release(overlap_chol)
473 CALL cp_fm_release(density_orth)
474 CALL cp_fm_release(direction_orth)
475 CALL cp_fm_release(gradient_orth)
476 CALL cp_fm_release(previous_density)
477 CALL cp_fm_release(previous_gradient)
478 CALL cp_fm_release(trial_orth)
479 CALL cp_fm_release(eigenvectors)
480 IF (use_relative_entropy) THEN
481 CALL cp_fm_release(hamiltonian_orth)
482 CALL cp_fm_release(prior_logit)
483 END IF
484 CALL cp_fm_release(work1)
485 CALL cp_fm_release(work2)
486 CALL dbcsr_release(gradient_ao%matrix)
487 DEALLOCATE (gradient_ao%matrix)
488 CALL auxbas_pw_pool%give_back_pw(residual_rspace)
489 CALL auxbas_pw_pool%give_back_pw(target_rspace)
490
491 CALL timestop(handle)
492
493 END SUBROUTINE fit_density_matrix
494
495! **************************************************************************************************
496!> \brief Builds fixed-trace, spin-restricted Fermi occupations for a set of eigenvalues.
497!> \param eigenvalues ...
498!> \param trace_target ...
499!> \param temperature ...
500!> \param chemical_potential ...
501!> \param occupations ...
502! **************************************************************************************************
503 SUBROUTINE fixed_trace_fermi_occupations(eigenvalues, trace_target, temperature, &
504 chemical_potential, occupations)
505 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenvalues
506 REAL(kind=dp), INTENT(IN) :: trace_target, temperature
507 REAL(kind=dp), INTENT(OUT) :: chemical_potential
508 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: occupations
509
510 INTEGER :: i, iteration, n
511 REAL(kind=dp) :: mu_high, mu_low, trace_value
512
513 n = SIZE(eigenvalues)
514 cpassert(SIZE(occupations) == n)
515 cpassert(temperature > 0.0_dp)
516 cpassert(trace_target >= 0.0_dp .AND. trace_target <= 2.0_dp*real(n, dp))
517
518 IF (trace_target <= 0.0_dp) THEN
519 occupations = 0.0_dp
520 chemical_potential = minval(eigenvalues) - 100.0_dp*temperature
521 RETURN
522 ELSE IF (trace_target >= 2.0_dp*real(n, dp)) THEN
523 occupations = 2.0_dp
524 chemical_potential = maxval(eigenvalues) + 100.0_dp*temperature
525 RETURN
526 END IF
527
528 mu_low = minval(eigenvalues) - 100.0_dp*temperature - 1.0_dp
529 mu_high = maxval(eigenvalues) + 100.0_dp*temperature + 1.0_dp
530 DO iteration = 1, 200
531 chemical_potential = 0.5_dp*(mu_low + mu_high)
532 DO i = 1, n
533 occupations(i) = fermi_occupation(eigenvalues(i), chemical_potential, temperature)
534 END DO
535 trace_value = sum(occupations)
536 IF (trace_value > trace_target) THEN
537 mu_high = chemical_potential
538 ELSE
539 mu_low = chemical_potential
540 END IF
541 IF (abs(trace_value - trace_target) <= 1.0e-13_dp*max(1.0_dp, trace_target)) EXIT
542 END DO
543 chemical_potential = 0.5_dp*(mu_low + mu_high)
544 DO i = 1, n
545 occupations(i) = fermi_occupation(eigenvalues(i), chemical_potential, temperature)
546 END DO
547
548 END SUBROUTINE fixed_trace_fermi_occupations
549
550! **************************************************************************************************
551!> \brief Numerically stable spin-restricted Fermi occupation.
552!> \param energy ...
553!> \param chemical_potential ...
554!> \param temperature ...
555!> \return ...
556! **************************************************************************************************
557 PURE FUNCTION fermi_occupation(energy, chemical_potential, temperature) RESULT(occupation)
558 REAL(kind=dp), INTENT(IN) :: energy, chemical_potential, temperature
559 REAL(kind=dp) :: occupation
560
561 REAL(kind=dp) :: x
562
563 x = (energy - chemical_potential)/temperature
564 IF (x >= 50.0_dp) THEN
565 occupation = 2.0_dp*exp(-x)
566 ELSE IF (x <= -50.0_dp) THEN
567 occupation = 2.0_dp
568 ELSE
569 occupation = 2.0_dp/(1.0_dp + exp(x))
570 END IF
571
572 END FUNCTION fermi_occupation
573
574! **************************************************************************************************
575!> \brief Converts a spin-restricted occupation to a numerically interior fraction.
576!>
577!> At ordinary electronic temperatures a Fermi occupation can round to exactly zero or two.
578!> Clipping its per-spin fraction keeps the matrix logit finite and, because the same operation
579!> is used for the prior and every trial matrix, preserves a zero relative entropy at the prior.
580!> \param occupation spin-restricted occupation in [0,2]
581!> \return clipped per-spin occupation in (0,1)
582! **************************************************************************************************
583 PURE FUNCTION clipped_occupation_fraction(occupation) RESULT(fraction)
584 REAL(kind=dp), INTENT(IN) :: occupation
585 REAL(kind=dp) :: fraction
586
587 REAL(kind=dp), PARAMETER :: occupation_clip = 1.0e-14_dp
588
589 fraction = min(1.0_dp - occupation_clip, &
590 max(occupation_clip, 0.5_dp*occupation))
591
592 END FUNCTION clipped_occupation_fraction
593
594! **************************************************************************************************
595!> \brief Reconstructs V diag(values) V^T.
596!> \param eigenvectors ...
597!> \param values ...
598!> \param matrix ...
599!> \param work ...
600! **************************************************************************************************
601 SUBROUTINE matrix_from_eigensystem(eigenvectors, values, matrix, work)
602 TYPE(cp_fm_type), INTENT(IN) :: eigenvectors
603 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: values
604 TYPE(cp_fm_type), INTENT(INOUT) :: matrix, work
605
606 INTEGER :: n
607
608 n = SIZE(values)
609 CALL cp_fm_to_fm(eigenvectors, work)
610 CALL cp_fm_column_scale(work, values)
611 CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, work, eigenvectors, 0.0_dp, matrix)
612
613 END SUBROUTINE matrix_from_eigensystem
614
615! **************************************************************************************************
616!> \brief Evaluates the fermionic quantum relative entropy D_F(X/2 || f_0).
617!>
618!> The matrix logit of the representable Fermi prior is precomputed spectrally.
619!> \param density_orth ...
620!> \param prior_logit ...
621!> \param eigenvectors ...
622!> \param work ...
623!> \param eigenvalues ...
624!> \param prior_log_one_minus ...
625!> \param relative_entropy ...
626! **************************************************************************************************
627 SUBROUTINE fermionic_relative_entropy(density_orth, prior_logit, eigenvectors, work, eigenvalues, &
628 prior_log_one_minus, relative_entropy)
629 TYPE(cp_fm_type), INTENT(IN) :: density_orth, prior_logit
630 TYPE(cp_fm_type), INTENT(INOUT) :: eigenvectors, work
631 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: eigenvalues
632 REAL(kind=dp), INTENT(IN) :: prior_log_one_minus
633 REAL(kind=dp), INTENT(OUT) :: relative_entropy
634
635 INTEGER :: i
636 REAL(kind=dp) :: entropy_part, f, trace_xlogit
637
638 CALL cp_fm_to_fm(density_orth, work)
639 CALL cp_fm_syevd(work, eigenvectors, eigenvalues)
640 entropy_part = 0.0_dp
641 DO i = 1, SIZE(eigenvalues)
642 f = clipped_occupation_fraction(eigenvalues(i))
643 entropy_part = entropy_part + f*log(f) + (1.0_dp - f)*log(1.0_dp - f)
644 END DO
645 CALL cp_fm_trace(density_orth, prior_logit, trace_xlogit)
646 relative_entropy = entropy_part - prior_log_one_minus - 0.5_dp*trace_xlogit
647 IF (relative_entropy < 0.0_dp .AND. abs(relative_entropy) < 1.0e-10_dp) relative_entropy = 0.0_dp
648
649 END SUBROUTINE fermionic_relative_entropy
650
651! **************************************************************************************************
652!> \brief Adds weight*dD_F/dX to a density-fit gradient at fixed trace.
653!> \param density_orth ...
654!> \param prior_logit ...
655!> \param eigenvectors ...
656!> \param work1 ...
657!> \param work2 ...
658!> \param eigenvalues ...
659!> \param weight ...
660!> \param gradient ...
661! **************************************************************************************************
662 SUBROUTINE add_relative_entropy_gradient(density_orth, prior_logit, eigenvectors, &
663 work1, work2, eigenvalues, weight, gradient)
664 TYPE(cp_fm_type), INTENT(IN) :: density_orth, prior_logit
665 TYPE(cp_fm_type), INTENT(INOUT) :: eigenvectors, work1, work2
666 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: eigenvalues
667 REAL(kind=dp), INTENT(IN) :: weight
668 TYPE(cp_fm_type), INTENT(INOUT) :: gradient
669
670 INTEGER :: i
671 REAL(kind=dp) :: f
672
673 CALL cp_fm_to_fm(density_orth, work1)
674 CALL cp_fm_syevd(work1, eigenvectors, eigenvalues)
675 DO i = 1, SIZE(eigenvalues)
676 f = clipped_occupation_fraction(eigenvalues(i))
677 eigenvalues(i) = 0.5_dp*log(f/(1.0_dp - f))
678 END DO
679 CALL matrix_from_eigensystem(eigenvectors, eigenvalues, work1, work2)
680 CALL cp_fm_scale_and_add(1.0_dp, gradient, beta=weight, matrix_b=work1)
681 CALL cp_fm_scale_and_add(1.0_dp, gradient, beta=-0.5_dp*weight, matrix_b=prior_logit)
682
683 END SUBROUTINE add_relative_entropy_gradient
684
685! **************************************************************************************************
686!> \brief Frobenius norm of the commutator [X,H] in the orthonormal AO representation.
687!> \param density_orth ...
688!> \param hamiltonian_orth ...
689!> \param work1 ...
690!> \param work2 ...
691!> \param nao ...
692!> \param norm ...
693! **************************************************************************************************
694 SUBROUTINE matrix_commutator_norm(density_orth, hamiltonian_orth, work1, work2, nao, norm)
695 TYPE(cp_fm_type), INTENT(IN) :: density_orth, hamiltonian_orth
696 TYPE(cp_fm_type), INTENT(INOUT) :: work1, work2
697 INTEGER, INTENT(IN) :: nao
698 REAL(kind=dp), INTENT(OUT) :: norm
699
700 REAL(kind=dp) :: norm_squared
701
702 CALL parallel_gemm('N', 'N', nao, nao, nao, 1.0_dp, density_orth, hamiltonian_orth, &
703 0.0_dp, work1)
704 CALL parallel_gemm('N', 'N', nao, nao, nao, 1.0_dp, hamiltonian_orth, density_orth, &
705 0.0_dp, work2)
706 CALL cp_fm_scale_and_add(1.0_dp, work1, beta=-1.0_dp, matrix_b=work2)
707 CALL cp_fm_trace(work1, work1, norm_squared)
708 norm = sqrt(max(0.0_dp, norm_squared))
709
710 END SUBROUTINE matrix_commutator_norm
711
712! **************************************************************************************************
713!> \brief Frobenius norm of f^2-f for f=X/2.
714!> \param density_orth ...
715!> \param work1 ...
716!> \param work2 ...
717!> \param nao ...
718!> \param error ...
719! **************************************************************************************************
720 SUBROUTINE density_matrix_idempotency_error(density_orth, work1, work2, nao, error)
721 TYPE(cp_fm_type), INTENT(IN) :: density_orth
722 TYPE(cp_fm_type), INTENT(INOUT) :: work1, work2
723 INTEGER, INTENT(IN) :: nao
724 REAL(kind=dp), INTENT(OUT) :: error
725
726 REAL(kind=dp) :: error_squared
727
728 CALL parallel_gemm('N', 'N', nao, nao, nao, 1.0_dp, density_orth, density_orth, &
729 0.0_dp, work1)
730 CALL cp_fm_to_fm(work1, work2)
731 CALL cp_fm_scale_and_add(0.25_dp, work2, beta=-0.5_dp, matrix_b=density_orth)
732 CALL cp_fm_trace(work2, work2, error_squared)
733 error = sqrt(max(0.0_dp, error_squared))
734
735 END SUBROUTINE density_matrix_idempotency_error
736
737! **************************************************************************************************
738!> \brief Projects a symmetric matrix onto eigenvalue bounds and a prescribed trace.
739!> \param matrix ...
740!> \param eigenvectors ...
741!> \param work ...
742!> \param eigenvalues ...
743!> \param trace_target ...
744!> \param max_occupation ...
745! **************************************************************************************************
746 SUBROUTINE project_density_matrix(matrix, eigenvectors, work, eigenvalues, trace_target, max_occupation)
747 TYPE(cp_fm_type), INTENT(INOUT) :: matrix, eigenvectors, work
748 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: eigenvalues
749 REAL(kind=dp), INTENT(IN) :: trace_target, max_occupation
750
751 INTEGER :: i, iteration, n
752 REAL(kind=dp) :: tau, tau_high, tau_low, trace_value
753
754 n = SIZE(eigenvalues)
755 cpassert(trace_target >= 0.0_dp .AND. trace_target <= max_occupation*real(n, dp))
756
757 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues)
758
759 IF (trace_target <= 0.0_dp) THEN
760 eigenvalues = 0.0_dp
761 ELSE IF (trace_target >= max_occupation*real(n, dp)) THEN
762 eigenvalues = max_occupation
763 ELSE
764 tau_low = minval(eigenvalues) - max_occupation
765 tau_high = maxval(eigenvalues)
766 DO iteration = 1, 200
767 tau = 0.5_dp*(tau_low + tau_high)
768 trace_value = sum(min(max_occupation, max(0.0_dp, eigenvalues - tau)))
769 IF (trace_value > trace_target) THEN
770 tau_low = tau
771 ELSE
772 tau_high = tau
773 END IF
774 IF (tau_high - tau_low <= 10.0_dp*epsilon(tau)*max(1.0_dp, abs(tau))) EXIT
775 END DO
776 tau = 0.5_dp*(tau_low + tau_high)
777 DO i = 1, n
778 eigenvalues(i) = min(max_occupation, max(0.0_dp, eigenvalues(i) - tau))
779 END DO
780 END IF
781
782 CALL cp_fm_to_fm(eigenvectors, work)
783 CALL cp_fm_column_scale(work, eigenvalues)
784 CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, work, eigenvectors, 0.0_dp, matrix)
785
786 END SUBROUTINE project_density_matrix
787
788! **************************************************************************************************
789!> \brief Transforms X in the orthonormal representation back to P = U^-1 X U^-T.
790!> \param density_orth ...
791!> \param overlap_chol ...
792!> \param work1 ...
793!> \param density_ao ...
794!> \param nao ...
795! **************************************************************************************************
796 SUBROUTINE orthogonal_to_ao(density_orth, overlap_chol, work1, density_ao, nao)
797 TYPE(cp_fm_type), INTENT(IN) :: density_orth, overlap_chol
798 TYPE(cp_fm_type), INTENT(INOUT) :: work1, density_ao
799 INTEGER, INTENT(IN) :: nao
800
801 CALL cp_fm_cholesky_restore(density_orth, nao, overlap_chol, work1, "SOLVE", pos="LEFT")
802 CALL cp_fm_cholesky_restore(work1, nao, overlap_chol, density_ao, &
803 "SOLVE", pos="RIGHT", transa="T")
804
805 END SUBROUTINE orthogonal_to_ao
806
807! **************************************************************************************************
808!> \brief Collocates the current AO density and evaluates its least-squares residual.
809!> \param qs_env ...
810!> \param ks_env ...
811!> \param rho_ao ...
812!> \param rho_r ...
813!> \param rho_g ...
814!> \param tot_rho_r ...
815!> \param target_rspace ...
816!> \param residual_rspace ...
817!> \param objective ...
818!> \param rms ...
819! **************************************************************************************************
820 SUBROUTINE evaluate_density_fit(qs_env, ks_env, rho_ao, rho_r, rho_g, tot_rho_r, &
821 target_rspace, residual_rspace, objective, rms)
822 TYPE(qs_environment_type), POINTER :: qs_env
823 TYPE(qs_ks_env_type), POINTER :: ks_env
824 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
825 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
826 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
827 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
828 TYPE(pw_r3d_rs_type), INTENT(IN) :: target_rspace
829 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: residual_rspace
830 REAL(kind=dp), INTENT(OUT) :: objective, rms
831
832 mark_used(qs_env)
833
834 CALL calculate_rho_elec(matrix_p=rho_ao(1)%matrix, rho=rho_r(1), &
835 rho_gspace=rho_g(1), total_rho=tot_rho_r(1), ks_env=ks_env)
836 CALL pw_copy(rho_r(1), residual_rspace)
837 CALL pw_axpy(target_rspace, residual_rspace, alpha=-1.0_dp, beta=1.0_dp)
838 objective = 0.5_dp*pw_integral_ab(residual_rspace, residual_rspace)
839 rms = sqrt(2.0_dp*objective/rho_r(1)%pw_grid%vol)
840
841 END SUBROUTINE evaluate_density_fit
842
843END MODULE qs_density_fit
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:78
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_integrate_function
Definition pw_methods.F:116
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1690
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:170
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1584
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
apply Cholesky decomposition op can be "SOLVE" (out = U^-1 * in) or "MULTIPLY" (out = U * in) pos can...
Definition cp_fm_cholesky.F:313
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_fm_cholesky.F:61
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_syevd
subroutine, public cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
Computes all eigenvalues and vectors of a real symmetric matrix significantly faster than syevx,...
Definition cp_fm_diag.F:590
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:576
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::calculate_rho_elec
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Definition qs_collocate_density.F:1720
qs_density_fit
Constrained fitting of a real-space density by an AO density matrix. This is a dense reference implem...
Definition qs_density_fit.F:13
qs_density_fit::fit_constrained_density
subroutine, public fit_constrained_density(qs_env, rho_struct, max_iter, eps_rms, step_size, max_backtrack)
Fits the density currently stored in rho_struct with a spin-restricted AO density matrix....
Definition qs_density_fit.F:86
qs_density_fit::fit_relative_entropy_density
subroutine, public fit_relative_entropy_density(qs_env, rho_struct, prior_hamiltonian, temperature, entropy_weight, max_iter, eps_rms, step_size, max_backtrack)
Reconstructs an AO density matrix using the Fermi matrix of H[n_cube] as a fermionic relative-entropy...
Definition qs_density_fit.F:111
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:539
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition qs_integrate_potential.F:22
qs_ks_types
Definition qs_ks_types.F:15
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:758
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:767
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62