Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 : ! **************************************************************************************************
13 : MODULE 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 :
72 : CONTAINS
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 8 : 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 8 : CALL fit_density_matrix(qs_env, rho_struct, max_iter, eps_rms, step_size, max_backtrack)
93 :
94 8 : 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 8 : 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 8 : prior_hamiltonian, temperature, entropy_weight)
121 :
122 8 : 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 16 : 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 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, prior_eigenvalues, &
158 16 : prior_logit_values, prior_occupations
159 16 : 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 16 : 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 16 : 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 16 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
175 : TYPE(qs_ks_env_type), POINTER :: ks_env
176 :
177 16 : CALL timeset(routineN, handle)
178 :
179 16 : use_relative_entropy = PRESENT(prior_hamiltonian)
180 16 : CPASSERT(use_relative_entropy .EQV. PRESENT(temperature))
181 16 : CPASSERT(use_relative_entropy .EQV. PRESENT(entropy_weight))
182 16 : IF (use_relative_entropy) THEN
183 8 : relative_temperature = temperature
184 8 : relative_weight = entropy_weight
185 8 : CPASSERT(relative_temperature > 0.0_dp)
186 8 : CPASSERT(relative_weight >= 0.0_dp)
187 : ELSE
188 8 : relative_temperature = 1.0_dp
189 8 : relative_weight = 0.0_dp
190 : END IF
191 :
192 16 : NULLIFY (auxbas_pw_pool, blacs_env, dft_control, fm_struct, gradient_ao%matrix, &
193 16 : 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 16 : nelectron_total=nelectron, para_env=para_env, pw_env=pw_env)
198 :
199 16 : IF (do_kpoints) CPABORT("Harris CUBE_FIT is currently available only at the Gamma point")
200 16 : IF (dft_control%nspins /= 1) THEN
201 0 : CPABORT("Harris CUBE_FIT currently requires a spin-restricted calculation")
202 : END IF
203 16 : IF (dft_control%qs_control%gapw) CPABORT("Harris CUBE_FIT currently supports GPW only")
204 :
205 16 : 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 16 : CPASSERT(ASSOCIATED(rho_ao) .AND. SIZE(rho_ao) == 1)
207 16 : CPASSERT(ASSOCIATED(rho_r) .AND. SIZE(rho_r) == 1)
208 16 : CPASSERT(ASSOCIATED(rho_g) .AND. SIZE(rho_g) == 1)
209 :
210 16 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
211 16 : IF (nelectron < 0 .OR. nelectron > 2*nao) THEN
212 0 : CPABORT("Electron count is incompatible with the constrained AO density-matrix bounds")
213 : END IF
214 :
215 16 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
216 16 : CALL auxbas_pw_pool%create_pw(target_rspace)
217 16 : CALL auxbas_pw_pool%create_pw(residual_rspace)
218 16 : CALL pw_copy(rho_r(1), target_rspace)
219 :
220 16 : ALLOCATE (gradient_ao%matrix)
221 : CALL dbcsr_create(gradient_ao%matrix, template=matrix_s(1)%matrix, &
222 16 : name="Harris density-fit gradient", matrix_type=dbcsr_type_symmetric)
223 16 : CALL dbcsr_copy(gradient_ao%matrix, matrix_s(1)%matrix)
224 16 : 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 16 : nrow_global=nao, ncol_global=nao)
228 16 : CALL cp_fm_create(overlap_chol, fm_struct, name="density-fit overlap")
229 16 : CALL cp_fm_create(density_orth, fm_struct, name="density fit in orthogonal basis")
230 16 : CALL cp_fm_create(direction_orth, fm_struct, name="density-fit search direction")
231 16 : CALL cp_fm_create(gradient_orth, fm_struct, name="density-fit gradient in orthogonal basis")
232 16 : CALL cp_fm_create(previous_density, fm_struct, name="previous density-fit matrix")
233 16 : CALL cp_fm_create(previous_gradient, fm_struct, name="previous density-fit gradient")
234 16 : CALL cp_fm_create(trial_orth, fm_struct, name="density-fit trial matrix")
235 16 : CALL cp_fm_create(eigenvectors, fm_struct, name="density-fit eigenvectors")
236 16 : IF (use_relative_entropy) THEN
237 8 : CALL cp_fm_create(hamiltonian_orth, fm_struct, name="relative-entropy prior Hamiltonian")
238 8 : CALL cp_fm_create(prior_logit, fm_struct, name="relative-entropy prior logit")
239 : END IF
240 16 : CALL cp_fm_create(work1, fm_struct, name="density-fit work matrix 1")
241 16 : CALL cp_fm_create(work2, fm_struct, name="density-fit work matrix 2")
242 16 : CALL cp_fm_struct_release(fm_struct)
243 48 : ALLOCATE (eigenvalues(nao))
244 16 : IF (use_relative_entropy) THEN
245 32 : ALLOCATE (prior_eigenvalues(nao), prior_logit_values(nao), prior_occupations(nao))
246 : END IF
247 :
248 16 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, overlap_chol)
249 16 : CALL cp_fm_uplo_to_full(overlap_chol, work1)
250 16 : CALL cp_fm_cholesky_decompose(overlap_chol, info_out=info)
251 16 : 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 16 : IF (use_relative_entropy) THEN
255 : ! H_tilde = U^-T H[n_cube] U^-1 in the same orthonormal representation.
256 8 : CALL copy_dbcsr_to_fm(prior_hamiltonian, work1)
257 8 : CALL cp_fm_uplo_to_full(work1, work2)
258 : CALL cp_fm_cholesky_restore(work1, nao, overlap_chol, work2, &
259 8 : "SOLVE", pos="LEFT", transa="T")
260 : CALL cp_fm_cholesky_restore(work2, nao, overlap_chol, hamiltonian_orth, &
261 8 : "SOLVE", pos="RIGHT")
262 8 : CALL cp_fm_to_fm(hamiltonian_orth, work1)
263 8 : CALL cp_fm_syevd(work1, eigenvectors, prior_eigenvalues)
264 : CALL fixed_trace_fermi_occupations(prior_eigenvalues, REAL(nelectron, dp), &
265 8 : relative_temperature, chemical_potential, prior_occupations)
266 8 : CALL matrix_from_eigensystem(eigenvectors, prior_occupations, density_orth, work1)
267 8 : prior_log_one_minus = 0.0_dp
268 48 : DO i = 1, nao
269 40 : prior_fraction = clipped_occupation_fraction(prior_occupations(i))
270 40 : prior_logit_values(i) = LOG(prior_fraction/(1.0_dp - prior_fraction))
271 48 : prior_log_one_minus = prior_log_one_minus + LOG(1.0_dp - prior_fraction)
272 : END DO
273 8 : CALL matrix_from_eigensystem(eigenvectors, prior_logit_values, prior_logit, work1)
274 : ELSE
275 8 : CALL copy_dbcsr_to_fm(rho_ao(1)%matrix, work1)
276 8 : CALL cp_fm_uplo_to_full(work1, work2)
277 8 : 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 8 : "MULTIPLY", pos="RIGHT", transa="T")
280 : CALL project_density_matrix(density_orth, eigenvectors, work1, eigenvalues, &
281 8 : REAL(nelectron, dp), 2.0_dp)
282 8 : chemical_potential = 0.0_dp
283 8 : prior_log_one_minus = 0.0_dp
284 : END IF
285 16 : CALL orthogonal_to_ao(density_orth, overlap_chol, work1, work2, nao)
286 16 : 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 16 : target_rspace, residual_rspace, objective, rms)
289 16 : relative_entropy = 0.0_dp
290 16 : IF (use_relative_entropy) THEN
291 : CALL fermionic_relative_entropy(density_orth, prior_logit, eigenvectors, work1, &
292 8 : eigenvalues, prior_log_one_minus, relative_entropy)
293 8 : objective = objective + relative_weight*relative_entropy
294 : END IF
295 :
296 16 : logger => cp_get_default_logger()
297 16 : IF (logger%para_env%is_source()) THEN
298 8 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
299 : ELSE
300 : unit_nr = -1
301 : END IF
302 8 : IF (unit_nr > 0) THEN
303 8 : IF (use_relative_entropy) THEN
304 : WRITE (unit_nr, "(/,T3,A,I0)") &
305 4 : "HARRIS| Relative-entropy AO density reconstruction; basis functions: ", nao
306 4 : WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Prior electronic temperature [a.u.]: ", &
307 8 : relative_temperature
308 4 : WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Prior chemical potential [a.u.]: ", &
309 8 : chemical_potential
310 4 : WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Relative-entropy weight: ", relative_weight
311 4 : WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Initial fermionic relative entropy: ", &
312 8 : relative_entropy
313 : ELSE
314 4 : WRITE (unit_nr, "(/,T3,A,I0)") "HARRIS| Constrained AO density fit; basis functions: ", nao
315 : END IF
316 8 : WRITE (unit_nr, "(T3,A,I0)") "HARRIS| Constrained electron count: ", nelectron
317 8 : WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Initial density RMS error: ", rms
318 8 : IF (use_relative_entropy) THEN
319 : WRITE (unit_nr, "(T3,A)") &
320 4 : "HARRIS| Iteration RMS error Relative entropy Objective Step"
321 : ELSE
322 4 : WRITE (unit_nr, "(T3,A)") "HARRIS| Iteration RMS error Objective Step"
323 : END IF
324 : END IF
325 :
326 16 : converged = rms <= eps_rms
327 16 : alpha = step_size
328 16 : have_previous = .FALSE.
329 16 : iteration = 0
330 48 : DO WHILE (iteration < max_iter .AND. .NOT. converged)
331 32 : iteration = iteration + 1
332 :
333 32 : CALL dbcsr_set(gradient_ao%matrix, 0.0_dp)
334 : CALL integrate_v_rspace(v_rspace=residual_rspace, hmat=gradient_ao, &
335 32 : qs_env=qs_env, calculate_forces=.FALSE.)
336 32 : CALL copy_dbcsr_to_fm(gradient_ao%matrix, work1)
337 32 : 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 32 : "SOLVE", pos="LEFT", transa="T")
342 : CALL cp_fm_cholesky_restore(work2, nao, overlap_chol, gradient_orth, &
343 32 : "SOLVE", pos="RIGHT")
344 :
345 32 : IF (use_relative_entropy .AND. relative_weight > 0.0_dp) THEN
346 : CALL add_relative_entropy_gradient(density_orth, prior_logit, eigenvectors, work1, &
347 16 : work2, eigenvalues, relative_weight, gradient_orth)
348 : END IF
349 :
350 : ! Barzilai-Borwein spectral step for the convex projected-gradient iteration.
351 32 : IF (have_previous) THEN
352 16 : CALL cp_fm_to_fm(density_orth, work1)
353 16 : CALL cp_fm_scale_and_add(1.0_dp, work1, beta=-1.0_dp, matrix_b=previous_density)
354 16 : CALL cp_fm_to_fm(gradient_orth, work2)
355 16 : CALL cp_fm_scale_and_add(1.0_dp, work2, beta=-1.0_dp, matrix_b=previous_gradient)
356 16 : CALL cp_fm_trace(work1, work1, ss)
357 16 : CALL cp_fm_trace(work1, work2, sy)
358 16 : IF (sy > 100.0_dp*EPSILON(sy)*MAX(1.0_dp, ss)) THEN
359 16 : alpha = MIN(100.0_dp*step_size, MAX(1.0E-6_dp*step_size, ss/sy))
360 : END IF
361 : END IF
362 32 : CALL cp_fm_to_fm(density_orth, previous_density)
363 32 : CALL cp_fm_to_fm(gradient_orth, previous_gradient)
364 32 : 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 32 : CALL cp_fm_to_fm(density_orth, trial_orth)
369 32 : 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 32 : REAL(nelectron, dp), 2.0_dp)
372 32 : CALL cp_fm_to_fm(trial_orth, direction_orth)
373 32 : CALL cp_fm_scale_and_add(1.0_dp, direction_orth, beta=-1.0_dp, matrix_b=density_orth)
374 32 : CALL cp_fm_trace(gradient_orth, direction_orth, direction_derivative)
375 32 : CALL cp_fm_trace(direction_orth, direction_orth, direction_norm2)
376 32 : IF (direction_norm2 <= 100.0_dp*EPSILON(direction_norm2)) THEN
377 0 : IF (unit_nr > 0) WRITE (unit_nr, "(T3,A)") &
378 0 : "HARRIS| Density fit stopped: projected gradient reached a stationary point"
379 : converged = .TRUE.
380 : EXIT
381 : END IF
382 :
383 32 : accepted = .FALSE.
384 32 : step_length = 1.0_dp
385 32 : DO backtrack = 0, max_backtrack
386 32 : CALL cp_fm_to_fm(density_orth, trial_orth)
387 32 : CALL cp_fm_scale_and_add(1.0_dp, trial_orth, beta=step_length, matrix_b=direction_orth)
388 32 : CALL orthogonal_to_ao(trial_orth, overlap_chol, work1, work2, nao)
389 32 : 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 32 : target_rspace, residual_rspace, objective_trial, rms_trial)
392 32 : relative_entropy_trial = 0.0_dp
393 32 : IF (use_relative_entropy) THEN
394 : CALL fermionic_relative_entropy(trial_orth, prior_logit, eigenvectors, work1, &
395 16 : eigenvalues, prior_log_one_minus, relative_entropy_trial)
396 16 : objective_trial = objective_trial + relative_weight*relative_entropy_trial
397 : END IF
398 32 : IF (objective_trial <= objective + 1.0E-4_dp*step_length*direction_derivative) THEN
399 : accepted = .TRUE.
400 : EXIT
401 : END IF
402 32 : step_length = 0.5_dp*step_length
403 : END DO
404 :
405 32 : IF (.NOT. accepted) THEN
406 : ! Leave both the AO matrix and the grids at the last accepted point.
407 0 : CALL orthogonal_to_ao(density_orth, overlap_chol, work1, work2, nao)
408 0 : 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 0 : target_rspace, residual_rspace, objective, rms)
411 0 : IF (use_relative_entropy) THEN
412 : CALL fermionic_relative_entropy(density_orth, prior_logit, eigenvectors, work1, &
413 0 : eigenvalues, prior_log_one_minus, relative_entropy)
414 : objective = objective + relative_weight*relative_entropy
415 : END IF
416 0 : IF (unit_nr > 0) WRITE (unit_nr, "(T3,A)") &
417 0 : "HARRIS| Density fit stopped: projected line search reached a stationary point"
418 16 : converged = .TRUE.
419 : EXIT
420 : END IF
421 :
422 32 : CALL cp_fm_to_fm(trial_orth, density_orth)
423 32 : objective = objective_trial
424 32 : rms = rms_trial
425 32 : relative_entropy = relative_entropy_trial
426 32 : IF (unit_nr > 0) THEN
427 16 : IF (use_relative_entropy) THEN
428 : WRITE (unit_nr, "(T3,A,I7,4ES19.9)") &
429 8 : "HARRIS| ", iteration, rms, relative_entropy, objective, alpha*step_length
430 : ELSE
431 : WRITE (unit_nr, "(T3,A,I7,3ES19.9)") &
432 8 : "HARRIS| ", iteration, rms, objective, alpha*step_length
433 : END IF
434 : END IF
435 48 : converged = rms <= eps_rms
436 : END DO
437 :
438 16 : 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 16 : fitted_grid_trace = pw_integrate_function(rho_r(1), isign=1)
442 16 : commutator_norm = 0.0_dp
443 16 : idempotency_error = 0.0_dp
444 16 : IF (use_relative_entropy) THEN
445 : CALL matrix_commutator_norm(density_orth, hamiltonian_orth, work1, work2, nao, &
446 8 : commutator_norm)
447 8 : CALL density_matrix_idempotency_error(density_orth, work1, work2, nao, idempotency_error)
448 : END IF
449 16 : IF (unit_nr > 0) THEN
450 8 : IF (rms <= eps_rms) THEN
451 8 : WRITE (unit_nr, "(T3,A,I0,A)") "HARRIS| Density fit reached its RMS target after ", iteration, " iterations"
452 0 : ELSE IF (iteration >= max_iter) THEN
453 0 : WRITE (unit_nr, "(T3,A,I0,A)") "HARRIS| Density fit stopped after ", iteration, " iterations"
454 : END IF
455 8 : WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Final density RMS error: ", rms
456 8 : WRITE (unit_nr, "(T3,A,F20.10)") "HARRIS| AO density-matrix electron count: ", matrix_trace
457 8 : WRITE (unit_nr, "(T3,A,F20.10)") "HARRIS| Fitted grid electron count: ", fitted_grid_trace
458 8 : IF (use_relative_entropy) THEN
459 4 : WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Final fermionic relative entropy: ", &
460 8 : relative_entropy
461 4 : WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Prior-Hamiltonian commutator norm: ", &
462 8 : commutator_norm
463 4 : WRITE (unit_nr, "(T3,A,ES16.8)") "HARRIS| Occupation idempotency error: ", &
464 8 : idempotency_error
465 : END IF
466 : END IF
467 :
468 16 : DEALLOCATE (eigenvalues)
469 16 : IF (use_relative_entropy) THEN
470 8 : DEALLOCATE (prior_eigenvalues, prior_logit_values, prior_occupations)
471 : END IF
472 16 : CALL cp_fm_release(overlap_chol)
473 16 : CALL cp_fm_release(density_orth)
474 16 : CALL cp_fm_release(direction_orth)
475 16 : CALL cp_fm_release(gradient_orth)
476 16 : CALL cp_fm_release(previous_density)
477 16 : CALL cp_fm_release(previous_gradient)
478 16 : CALL cp_fm_release(trial_orth)
479 16 : CALL cp_fm_release(eigenvectors)
480 16 : IF (use_relative_entropy) THEN
481 8 : CALL cp_fm_release(hamiltonian_orth)
482 8 : CALL cp_fm_release(prior_logit)
483 : END IF
484 16 : CALL cp_fm_release(work1)
485 16 : CALL cp_fm_release(work2)
486 16 : CALL dbcsr_release(gradient_ao%matrix)
487 16 : DEALLOCATE (gradient_ao%matrix)
488 16 : CALL auxbas_pw_pool%give_back_pw(residual_rspace)
489 16 : CALL auxbas_pw_pool%give_back_pw(target_rspace)
490 :
491 16 : CALL timestop(handle)
492 :
493 128 : 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 8 : SUBROUTINE fixed_trace_fermi_occupations(eigenvalues, trace_target, temperature, &
504 8 : 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 8 : n = SIZE(eigenvalues)
514 8 : CPASSERT(SIZE(occupations) == n)
515 8 : CPASSERT(temperature > 0.0_dp)
516 8 : CPASSERT(trace_target >= 0.0_dp .AND. trace_target <= 2.0_dp*REAL(n, dp))
517 :
518 8 : IF (trace_target <= 0.0_dp) THEN
519 0 : occupations = 0.0_dp
520 0 : chemical_potential = MINVAL(eigenvalues) - 100.0_dp*temperature
521 0 : RETURN
522 8 : ELSE IF (trace_target >= 2.0_dp*REAL(n, dp)) THEN
523 0 : occupations = 2.0_dp
524 0 : chemical_potential = MAXVAL(eigenvalues) + 100.0_dp*temperature
525 0 : RETURN
526 : END IF
527 :
528 48 : mu_low = MINVAL(eigenvalues) - 100.0_dp*temperature - 1.0_dp
529 48 : mu_high = MAXVAL(eigenvalues) + 100.0_dp*temperature + 1.0_dp
530 424 : DO iteration = 1, 200
531 424 : chemical_potential = 0.5_dp*(mu_low + mu_high)
532 2544 : DO i = 1, n
533 2544 : occupations(i) = fermi_occupation(eigenvalues(i), chemical_potential, temperature)
534 : END DO
535 2544 : trace_value = SUM(occupations)
536 424 : IF (trace_value > trace_target) THEN
537 : mu_high = chemical_potential
538 : ELSE
539 224 : mu_low = chemical_potential
540 : END IF
541 424 : IF (ABS(trace_value - trace_target) <= 1.0E-13_dp*MAX(1.0_dp, trace_target)) EXIT
542 : END DO
543 8 : chemical_potential = 0.5_dp*(mu_low + mu_high)
544 48 : DO i = 1, n
545 48 : 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 2160 : 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 2160 : x = (energy - chemical_potential)/temperature
564 2160 : IF (x >= 50.0_dp) THEN
565 160 : occupation = 2.0_dp*EXP(-x)
566 2000 : ELSE IF (x <= -50.0_dp) THEN
567 : occupation = 2.0_dp
568 : ELSE
569 2000 : occupation = 2.0_dp/(1.0_dp + EXP(x))
570 : END IF
571 :
572 2160 : 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 240 : 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 240 : MAX(occupation_clip, 0.5_dp*occupation))
591 :
592 240 : 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 32 : 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 32 : n = SIZE(values)
609 32 : CALL cp_fm_to_fm(eigenvectors, work)
610 32 : CALL cp_fm_column_scale(work, values)
611 32 : CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, work, eigenvectors, 0.0_dp, matrix)
612 :
613 32 : 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 24 : 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 24 : CALL cp_fm_to_fm(density_orth, work)
639 24 : CALL cp_fm_syevd(work, eigenvectors, eigenvalues)
640 24 : entropy_part = 0.0_dp
641 144 : DO i = 1, SIZE(eigenvalues)
642 120 : f = clipped_occupation_fraction(eigenvalues(i))
643 144 : entropy_part = entropy_part + f*LOG(f) + (1.0_dp - f)*LOG(1.0_dp - f)
644 : END DO
645 24 : CALL cp_fm_trace(density_orth, prior_logit, trace_xlogit)
646 24 : relative_entropy = entropy_part - prior_log_one_minus - 0.5_dp*trace_xlogit
647 24 : IF (relative_entropy < 0.0_dp .AND. ABS(relative_entropy) < 1.0E-10_dp) relative_entropy = 0.0_dp
648 :
649 24 : 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 16 : SUBROUTINE add_relative_entropy_gradient(density_orth, prior_logit, eigenvectors, &
663 16 : 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 16 : CALL cp_fm_to_fm(density_orth, work1)
674 16 : CALL cp_fm_syevd(work1, eigenvectors, eigenvalues)
675 96 : DO i = 1, SIZE(eigenvalues)
676 80 : f = clipped_occupation_fraction(eigenvalues(i))
677 96 : eigenvalues(i) = 0.5_dp*LOG(f/(1.0_dp - f))
678 : END DO
679 16 : CALL matrix_from_eigensystem(eigenvectors, eigenvalues, work1, work2)
680 16 : CALL cp_fm_scale_and_add(1.0_dp, gradient, beta=weight, matrix_b=work1)
681 16 : CALL cp_fm_scale_and_add(1.0_dp, gradient, beta=-0.5_dp*weight, matrix_b=prior_logit)
682 :
683 16 : 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 16 : 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 8 : 0.0_dp, work1)
704 : CALL parallel_gemm('N', 'N', nao, nao, nao, 1.0_dp, hamiltonian_orth, density_orth, &
705 8 : 0.0_dp, work2)
706 8 : CALL cp_fm_scale_and_add(1.0_dp, work1, beta=-1.0_dp, matrix_b=work2)
707 8 : CALL cp_fm_trace(work1, work1, norm_squared)
708 8 : norm = SQRT(MAX(0.0_dp, norm_squared))
709 :
710 8 : 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 16 : 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 8 : 0.0_dp, work1)
730 8 : CALL cp_fm_to_fm(work1, work2)
731 8 : CALL cp_fm_scale_and_add(0.25_dp, work2, beta=-0.5_dp, matrix_b=density_orth)
732 8 : CALL cp_fm_trace(work2, work2, error_squared)
733 8 : error = SQRT(MAX(0.0_dp, error_squared))
734 :
735 8 : 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 40 : 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 40 : n = SIZE(eigenvalues)
755 40 : CPASSERT(trace_target >= 0.0_dp .AND. trace_target <= max_occupation*REAL(n, dp))
756 :
757 40 : CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues)
758 :
759 40 : IF (trace_target <= 0.0_dp) THEN
760 0 : eigenvalues = 0.0_dp
761 40 : ELSE IF (trace_target >= max_occupation*REAL(n, dp)) THEN
762 0 : eigenvalues = max_occupation
763 : ELSE
764 240 : tau_low = MINVAL(eigenvalues) - max_occupation
765 240 : tau_high = MAXVAL(eigenvalues)
766 2048 : DO iteration = 1, 200
767 2048 : tau = 0.5_dp*(tau_low + tau_high)
768 12288 : trace_value = SUM(MIN(max_occupation, MAX(0.0_dp, eigenvalues - tau)))
769 2048 : IF (trace_value > trace_target) THEN
770 : tau_low = tau
771 : ELSE
772 880 : tau_high = tau
773 : END IF
774 2048 : IF (tau_high - tau_low <= 10.0_dp*EPSILON(tau)*MAX(1.0_dp, ABS(tau))) EXIT
775 : END DO
776 40 : tau = 0.5_dp*(tau_low + tau_high)
777 240 : DO i = 1, n
778 240 : eigenvalues(i) = MIN(max_occupation, MAX(0.0_dp, eigenvalues(i) - tau))
779 : END DO
780 : END IF
781 :
782 40 : CALL cp_fm_to_fm(eigenvectors, work)
783 40 : CALL cp_fm_column_scale(work, eigenvalues)
784 40 : CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, work, eigenvectors, 0.0_dp, matrix)
785 :
786 40 : 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 48 : 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 48 : 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 48 : "SOLVE", pos="RIGHT", transa="T")
804 :
805 48 : 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 48 : 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 48 : rho_gspace=rho_g(1), total_rho=tot_rho_r(1), ks_env=ks_env)
836 48 : CALL pw_copy(rho_r(1), residual_rspace)
837 48 : CALL pw_axpy(target_rspace, residual_rspace, alpha=-1.0_dp, beta=1.0_dp)
838 48 : objective = 0.5_dp*pw_integral_ab(residual_rspace, residual_rspace)
839 48 : rms = SQRT(2.0_dp*objective/rho_r(1)%pw_grid%vol)
840 :
841 48 : END SUBROUTINE evaluate_density_fit
842 :
843 : END MODULE qs_density_fit
|
