(git:ed6f26b)
Loading...
Searching...
No Matches
admm_dm_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Contains ADMM methods which only require the density matrix
10!> \par History
11!> 11.2014 created [Ole Schuett]
12!> \author Ole Schuett
13! **************************************************************************************************
14MODULE admm_dm_methods
15 USE admm_dm_types, ONLY: admm_dm_type,&
16 mcweeny_history_type
17 USE admm_types, ONLY: get_admm_env
18 USE cp_control_types, ONLY: dft_control_type
19 USE cp_dbcsr_api, ONLY: &
20 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
21 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
22 dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type
23 USE cp_dbcsr_contrib, ONLY: dbcsr_frobenius_norm
24 USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
25 USE cp_log_handling, ONLY: cp_logger_get_default_unit_nr
26 USE input_constants, ONLY: do_admm_basis_projection,&
27 do_admm_blocked_projection
28 USE iterate_matrix, ONLY: invert_hotelling
29 USE kinds, ONLY: dp
30 USE pw_types, ONLY: pw_c1d_gs_type,&
31 pw_r3d_rs_type
32 USE qs_collocate_density, ONLY: calculate_rho_elec
33 USE qs_environment_types, ONLY: get_qs_env,&
34 qs_environment_type
35 USE qs_ks_types, ONLY: qs_ks_env_type
36 USE qs_rho_types, ONLY: qs_rho_get,&
37 qs_rho_set,&
38 qs_rho_type
39 USE task_list_types, ONLY: task_list_type
40#include "./base/base_uses.f90"
41
42 IMPLICIT NONE
43 PRIVATE
44
45 PUBLIC :: admm_dm_calc_rho_aux, admm_dm_merge_ks_matrix
46
47 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'admm_dm_methods'
48
49CONTAINS
50
51! **************************************************************************************************
52!> \brief Entry methods: Calculates auxiliary density matrix from primary one.
53!> \param qs_env ...
54!> \author Ole Schuett
55! **************************************************************************************************
56 SUBROUTINE admm_dm_calc_rho_aux(qs_env)
57 TYPE(qs_environment_type), POINTER :: qs_env
58
59 CHARACTER(len=*), PARAMETER :: routinen = 'admm_dm_calc_rho_aux'
60
61 INTEGER :: handle
62 TYPE(admm_dm_type), POINTER :: admm_dm
63
64 NULLIFY (admm_dm)
65 CALL timeset(routinen, handle)
66 CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
67
68 SELECT CASE (admm_dm%method)
69 CASE (do_admm_basis_projection)
70 CALL map_dm_projection(qs_env)
71
72 CASE (do_admm_blocked_projection)
73 CALL map_dm_blocked(qs_env)
74
75 CASE DEFAULT
76 cpabort("admm_dm_calc_rho_aux: unknown method")
77 END SELECT
78
79 IF (admm_dm%purify) &
80 CALL purify_mcweeny(qs_env)
81
82 CALL update_rho_aux(qs_env)
83
84 CALL timestop(handle)
85 END SUBROUTINE admm_dm_calc_rho_aux
86
87! **************************************************************************************************
88!> \brief Entry methods: Merges auxiliary Kohn-Sham matrix into primary one.
89!> \param qs_env ...
90!> \author Ole Schuett
91! **************************************************************************************************
92 SUBROUTINE admm_dm_merge_ks_matrix(qs_env)
93 TYPE(qs_environment_type), POINTER :: qs_env
94
95 CHARACTER(LEN=*), PARAMETER :: routinen = 'admm_dm_merge_ks_matrix'
96
97 INTEGER :: handle
98 TYPE(admm_dm_type), POINTER :: admm_dm
99 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge
100
101 CALL timeset(routinen, handle)
102 NULLIFY (admm_dm, matrix_ks_merge)
103
104 CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
105
106 IF (admm_dm%purify) THEN
107 CALL revert_purify_mcweeny(qs_env, matrix_ks_merge)
108 ELSE
109 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit=matrix_ks_merge)
110 END IF
111
112 SELECT CASE (admm_dm%method)
113 CASE (do_admm_basis_projection)
114 CALL merge_dm_projection(qs_env, matrix_ks_merge)
115
116 CASE (do_admm_blocked_projection)
117 CALL merge_dm_blocked(qs_env, matrix_ks_merge)
118
119 CASE DEFAULT
120 cpabort("admm_dm_merge_ks_matrix: unknown method")
121 END SELECT
122
123 IF (admm_dm%purify) &
124 CALL dbcsr_deallocate_matrix_set(matrix_ks_merge)
125
126 CALL timestop(handle)
127
128 END SUBROUTINE admm_dm_merge_ks_matrix
129
130! **************************************************************************************************
131!> \brief Calculates auxiliary density matrix via basis projection.
132!> \param qs_env ...
133!> \author Ole Schuett
134! **************************************************************************************************
135 SUBROUTINE map_dm_projection(qs_env)
136 TYPE(qs_environment_type), POINTER :: qs_env
137
138 INTEGER :: ispin
139 LOGICAL :: s_mstruct_changed
140 REAL(kind=dp) :: threshold
141 TYPE(admm_dm_type), POINTER :: admm_dm
142 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux, matrix_s_mixed, rho_ao, &
143 rho_ao_aux
144 TYPE(dbcsr_type) :: matrix_s_aux_inv, matrix_tmp
145 TYPE(dft_control_type), POINTER :: dft_control
146 TYPE(qs_rho_type), POINTER :: rho, rho_aux
147
148 NULLIFY (dft_control, admm_dm, matrix_s_aux, matrix_s_mixed, rho, rho_aux)
149 NULLIFY (rho_ao, rho_ao_aux)
150
151 CALL get_qs_env(qs_env, dft_control=dft_control, s_mstruct_changed=s_mstruct_changed, rho=rho)
152 CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux, rho_aux_fit=rho_aux, &
153 matrix_s_aux_fit_vs_orb=matrix_s_mixed, admm_dm=admm_dm)
154
155 CALL qs_rho_get(rho, rho_ao=rho_ao)
156 CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux)
157
158 IF (s_mstruct_changed) THEN
159 ! Calculate A = S_aux^(-1) * S_mixed
160 CALL dbcsr_create(matrix_s_aux_inv, template=matrix_s_aux(1)%matrix, matrix_type="N")
161 threshold = max(admm_dm%eps_filter, 1.0e-12_dp)
162 CALL invert_hotelling(matrix_s_aux_inv, matrix_s_aux(1)%matrix, threshold)
163
164 IF (.NOT. ASSOCIATED(admm_dm%matrix_A)) THEN
165 ALLOCATE (admm_dm%matrix_A)
166 CALL dbcsr_create(admm_dm%matrix_A, template=matrix_s_mixed(1)%matrix, matrix_type="N")
167 END IF
168 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_aux_inv, matrix_s_mixed(1)%matrix, &
169 0.0_dp, admm_dm%matrix_A)
170 CALL dbcsr_release(matrix_s_aux_inv)
171 END IF
172
173 ! Calculate P_aux = A * P * A^T
174 CALL dbcsr_create(matrix_tmp, template=admm_dm%matrix_A)
175 DO ispin = 1, dft_control%nspins
176 CALL dbcsr_multiply("N", "N", 1.0_dp, admm_dm%matrix_A, rho_ao(ispin)%matrix, &
177 0.0_dp, matrix_tmp)
178 CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp, admm_dm%matrix_A, &
179 0.0_dp, rho_ao_aux(ispin)%matrix)
180 END DO
181 CALL dbcsr_release(matrix_tmp)
182
183 END SUBROUTINE map_dm_projection
184
185! **************************************************************************************************
186!> \brief Calculates auxiliary density matrix via blocking.
187!> \param qs_env ...
188!> \author Ole Schuett
189! **************************************************************************************************
190 SUBROUTINE map_dm_blocked(qs_env)
191 TYPE(qs_environment_type), POINTER :: qs_env
192
193 INTEGER :: iatom, ispin, jatom
194 LOGICAL :: found
195 REAL(dp), DIMENSION(:, :), POINTER :: sparse_block, sparse_block_aux
196 TYPE(admm_dm_type), POINTER :: admm_dm
197 TYPE(dbcsr_iterator_type) :: iter
198 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_ao_aux
199 TYPE(dft_control_type), POINTER :: dft_control
200 TYPE(qs_rho_type), POINTER :: rho, rho_aux
201
202 NULLIFY (dft_control, admm_dm, rho, rho_aux, rho_ao, rho_ao_aux)
203
204 CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
205 CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux, admm_dm=admm_dm)
206
207 CALL qs_rho_get(rho, rho_ao=rho_ao)
208 CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux)
209
210 ! ** set blocked density matrix to 0
211 DO ispin = 1, dft_control%nspins
212 CALL dbcsr_set(rho_ao_aux(ispin)%matrix, 0.0_dp)
213 ! ** now loop through the list and copy corresponding blocks
214 CALL dbcsr_iterator_start(iter, rho_ao(ispin)%matrix)
215 DO WHILE (dbcsr_iterator_blocks_left(iter))
216 CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block)
217 IF (admm_dm%block_map(iatom, jatom) == 1) THEN
218 CALL dbcsr_get_block_p(rho_ao_aux(ispin)%matrix, &
219 row=iatom, col=jatom, block=sparse_block_aux, found=found)
220 IF (found) &
221 sparse_block_aux = sparse_block
222 END IF
223 END DO
224 CALL dbcsr_iterator_stop(iter)
225 END DO
226
227 END SUBROUTINE map_dm_blocked
228
229! **************************************************************************************************
230!> \brief Call calculate_rho_elec() for auxiliary density
231!> \param qs_env ...
232! **************************************************************************************************
233 SUBROUTINE update_rho_aux(qs_env)
234 TYPE(qs_environment_type), POINTER :: qs_env
235
236 INTEGER :: ispin
237 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r_aux
238 TYPE(admm_dm_type), POINTER :: admm_dm
239 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_aux
240 TYPE(dft_control_type), POINTER :: dft_control
241 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
242 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_aux
243 TYPE(qs_ks_env_type), POINTER :: ks_env
244 TYPE(qs_rho_type), POINTER :: rho_aux
245 TYPE(task_list_type), POINTER :: task_list_aux_fit
246
247 NULLIFY (dft_control, admm_dm, rho_aux, rho_ao_aux, rho_r_aux, rho_g_aux, tot_rho_r_aux, &
248 task_list_aux_fit, ks_env)
249
250 CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control)
251 CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list_aux_fit, rho_aux_fit=rho_aux, &
252 admm_dm=admm_dm)
253
254 CALL qs_rho_get(rho_aux, &
255 rho_ao=rho_ao_aux, &
256 rho_r=rho_r_aux, &
257 rho_g=rho_g_aux, &
258 tot_rho_r=tot_rho_r_aux)
259
260 DO ispin = 1, dft_control%nspins
261 CALL calculate_rho_elec(ks_env=ks_env, &
262 matrix_p=rho_ao_aux(ispin)%matrix, &
263 rho=rho_r_aux(ispin), &
264 rho_gspace=rho_g_aux(ispin), &
265 total_rho=tot_rho_r_aux(ispin), &
266 soft_valid=.false., &
267 basis_type="AUX_FIT", &
268 task_list_external=task_list_aux_fit)
269 END DO
270
271 CALL qs_rho_set(rho_aux, rho_r_valid=.true., rho_g_valid=.true.)
272
273 END SUBROUTINE update_rho_aux
274
275! **************************************************************************************************
276!> \brief Merges auxiliary Kohn-Sham matrix via basis projection.
277!> \param qs_env ...
278!> \param matrix_ks_merge Input: The KS matrix to be merged
279!> \author Ole Schuett
280! **************************************************************************************************
281 SUBROUTINE merge_dm_projection(qs_env, matrix_ks_merge)
282 TYPE(qs_environment_type), POINTER :: qs_env
283 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge
284
285 INTEGER :: ispin
286 TYPE(admm_dm_type), POINTER :: admm_dm
287 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
288 TYPE(dbcsr_type) :: matrix_tmp
289 TYPE(dft_control_type), POINTER :: dft_control
290
291 NULLIFY (admm_dm, dft_control, matrix_ks)
292
293 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks)
294 CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
295
296 ! Calculate K += A^T * K_aux * A
297 CALL dbcsr_create(matrix_tmp, template=admm_dm%matrix_A, matrix_type="N")
298
299 DO ispin = 1, dft_control%nspins
300 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ks_merge(ispin)%matrix, admm_dm%matrix_A, &
301 0.0_dp, matrix_tmp)
302 CALL dbcsr_multiply("T", "N", 1.0_dp, admm_dm%matrix_A, matrix_tmp, &
303 1.0_dp, matrix_ks(ispin)%matrix)
304 END DO
305
306 CALL dbcsr_release(matrix_tmp)
307
308 END SUBROUTINE merge_dm_projection
309
310! **************************************************************************************************
311!> \brief Merges auxiliary Kohn-Sham matrix via blocking.
312!> \param qs_env ...
313!> \param matrix_ks_merge Input: The KS matrix to be merged
314!> \author Ole Schuett
315! **************************************************************************************************
316 SUBROUTINE merge_dm_blocked(qs_env, matrix_ks_merge)
317 TYPE(qs_environment_type), POINTER :: qs_env
318 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge
319
320 INTEGER :: iatom, ispin, jatom
321 REAL(dp), DIMENSION(:, :), POINTER :: sparse_block
322 TYPE(admm_dm_type), POINTER :: admm_dm
323 TYPE(dbcsr_iterator_type) :: iter
324 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
325 TYPE(dft_control_type), POINTER :: dft_control
326
327 NULLIFY (admm_dm, dft_control, matrix_ks)
328
329 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks)
330 CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
331
332 DO ispin = 1, dft_control%nspins
333 CALL dbcsr_iterator_start(iter, matrix_ks_merge(ispin)%matrix)
334 DO WHILE (dbcsr_iterator_blocks_left(iter))
335 CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block)
336 IF (admm_dm%block_map(iatom, jatom) == 0) &
337 sparse_block = 0.0_dp
338 END DO
339 CALL dbcsr_iterator_stop(iter)
340 CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_ks_merge(ispin)%matrix, 1.0_dp, 1.0_dp)
341 END DO
342
343 END SUBROUTINE merge_dm_blocked
344
345! **************************************************************************************************
346!> \brief Apply McWeeny purification to auxiliary density matrix
347!> \param qs_env ...
348!> \author Ole Schuett
349! **************************************************************************************************
350 SUBROUTINE purify_mcweeny(qs_env)
351 TYPE(qs_environment_type), POINTER :: qs_env
352
353 CHARACTER(LEN=*), PARAMETER :: routinen = 'purify_mcweeny'
354
355 INTEGER :: handle, ispin, istep, nspins, unit_nr
356 REAL(kind=dp) :: frob_norm
357 TYPE(admm_dm_type), POINTER :: admm_dm
358 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, rho_ao_aux
359 TYPE(dbcsr_type) :: matrix_ps, matrix_psp, matrix_test
360 TYPE(dbcsr_type), POINTER :: matrix_p, matrix_s
361 TYPE(dft_control_type), POINTER :: dft_control
362 TYPE(mcweeny_history_type), POINTER :: history, new_hist_entry
363 TYPE(qs_rho_type), POINTER :: rho_aux_fit
364
365 CALL timeset(routinen, handle)
366 NULLIFY (dft_control, admm_dm, matrix_s_aux_fit, rho_aux_fit, new_hist_entry, &
367 matrix_p, matrix_s, rho_ao_aux)
368
369 unit_nr = cp_logger_get_default_unit_nr()
370 CALL get_qs_env(qs_env, dft_control=dft_control)
371 CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux_fit, &
372 rho_aux_fit=rho_aux_fit, admm_dm=admm_dm)
373
374 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux)
375
376 matrix_p => rho_ao_aux(1)%matrix
377 CALL dbcsr_create(matrix_ps, template=matrix_p, matrix_type="N")
378 CALL dbcsr_create(matrix_psp, template=matrix_p, matrix_type="S")
379 CALL dbcsr_create(matrix_test, template=matrix_p, matrix_type="S")
380
381 nspins = dft_control%nspins
382 DO ispin = 1, nspins
383 matrix_p => rho_ao_aux(ispin)%matrix
384 matrix_s => matrix_s_aux_fit(1)%matrix
385 history => admm_dm%mcweeny_history(ispin)%p
386 IF (ASSOCIATED(history)) cpabort("purify_dm_mcweeny: history already associated")
387 IF (nspins == 1) CALL dbcsr_scale(matrix_p, 0.5_dp)
388
389 DO istep = 1, admm_dm%mcweeny_max_steps
390 ! allocate new element in linked list
391 ALLOCATE (new_hist_entry)
392 new_hist_entry%next => history
393 history => new_hist_entry
394 history%count = istep
395 NULLIFY (new_hist_entry)
396 CALL dbcsr_create(history%m, template=matrix_p, matrix_type="N")
397 CALL dbcsr_copy(history%m, matrix_p, name="P from McWeeny")
398
399 ! calc PS and PSP
400 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, matrix_s, &
401 0.0_dp, matrix_ps)
402
403 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ps, matrix_p, &
404 0.0_dp, matrix_psp)
405
406 !test convergence
407 CALL dbcsr_copy(matrix_test, matrix_psp)
408 CALL dbcsr_add(matrix_test, matrix_p, 1.0_dp, -1.0_dp)
409 frob_norm = dbcsr_frobenius_norm(matrix_test)
410 IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,i5,a,f16.8)') "McWeeny-Step", istep, &
411 ": Deviation of idempotency", frob_norm
412 IF (frob_norm < 1000_dp*admm_dm%eps_filter .AND. istep > 1) EXIT
413
414 ! build next P matrix
415 CALL dbcsr_copy(matrix_p, matrix_psp, name="P from McWeeny")
416 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_ps, matrix_psp, &
417 3.0_dp, matrix_p)
418 END DO
419 admm_dm%mcweeny_history(ispin)%p => history
420 IF (nspins == 1) CALL dbcsr_scale(matrix_p, 2.0_dp)
421 END DO
422
423 ! clean up
424 CALL dbcsr_release(matrix_ps)
425 CALL dbcsr_release(matrix_psp)
426 CALL dbcsr_release(matrix_test)
427 CALL timestop(handle)
428 END SUBROUTINE purify_mcweeny
429
430! **************************************************************************************************
431!> \brief Prepare auxiliary KS-matrix for merge using reverse McWeeny
432!> \param qs_env ...
433!> \param matrix_ks_merge Output: The KS matrix for the merge
434!> \author Ole Schuett
435! **************************************************************************************************
436 SUBROUTINE revert_purify_mcweeny(qs_env, matrix_ks_merge)
437 TYPE(qs_environment_type), POINTER :: qs_env
438 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge
439
440 CHARACTER(LEN=*), PARAMETER :: routinen = 'revert_purify_mcweeny'
441
442 INTEGER :: handle, ispin, nspins, unit_nr
443 TYPE(admm_dm_type), POINTER :: admm_dm
444 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
445 matrix_s_aux_fit, &
446 matrix_s_aux_fit_vs_orb
447 TYPE(dbcsr_type), POINTER :: matrix_k
448 TYPE(dft_control_type), POINTER :: dft_control
449 TYPE(mcweeny_history_type), POINTER :: history_curr, history_next
450
451 CALL timeset(routinen, handle)
452 unit_nr = cp_logger_get_default_unit_nr()
453 NULLIFY (admm_dm, dft_control, matrix_ks, matrix_ks_aux_fit, &
454 matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, &
455 history_next, history_curr, matrix_k)
456
457 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks)
458 CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux_fit, admm_dm=admm_dm, &
459 matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, matrix_ks_aux_fit=matrix_ks_aux_fit)
460
461 nspins = dft_control%nspins
462 ALLOCATE (matrix_ks_merge(nspins))
463
464 DO ispin = 1, nspins
465 ALLOCATE (matrix_ks_merge(ispin)%matrix)
466 matrix_k => matrix_ks_merge(ispin)%matrix
467 CALL dbcsr_copy(matrix_k, matrix_ks_aux_fit(ispin)%matrix, name="K")
468 history_curr => admm_dm%mcweeny_history(ispin)%p
469 NULLIFY (admm_dm%mcweeny_history(ispin)%p)
470
471 ! reverse McWeeny iteration
472 DO WHILE (ASSOCIATED(history_curr))
473 IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,i5)') "Reverse McWeeny-Step ", history_curr%count
474 CALL reverse_mcweeny_step(matrix_k=matrix_k, &
475 matrix_s=matrix_s_aux_fit(1)%matrix, &
476 matrix_p=history_curr%m)
477 CALL dbcsr_release(history_curr%m)
478 history_next => history_curr%next
479 DEALLOCATE (history_curr)
480 history_curr => history_next
481 NULLIFY (history_next)
482 END DO
483
484 END DO
485
486 ! clean up
487 CALL timestop(handle)
488
489 END SUBROUTINE revert_purify_mcweeny
490
491! **************************************************************************************************
492!> \brief Multiply matrix_k with partial derivative of McWeeny by reversing it.
493!> \param matrix_k ...
494!> \param matrix_s ...
495!> \param matrix_p ...
496!> \author Ole Schuett
497! **************************************************************************************************
498 SUBROUTINE reverse_mcweeny_step(matrix_k, matrix_s, matrix_p)
499 TYPE(dbcsr_type) :: matrix_k, matrix_s, matrix_p
500
501 CHARACTER(LEN=*), PARAMETER :: routinen = 'reverse_mcweeny_step'
502
503 INTEGER :: handle
504 TYPE(dbcsr_type) :: matrix_ps, matrix_sp, matrix_sum, &
505 matrix_tmp
506
507 CALL timeset(routinen, handle)
508 CALL dbcsr_create(matrix_ps, template=matrix_p, matrix_type="N")
509 CALL dbcsr_create(matrix_sp, template=matrix_p, matrix_type="N")
510 CALL dbcsr_create(matrix_tmp, template=matrix_p, matrix_type="N")
511 CALL dbcsr_create(matrix_sum, template=matrix_p, matrix_type="N")
512
513 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, matrix_s, &
514 0.0_dp, matrix_ps)
515 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, matrix_p, &
516 0.0_dp, matrix_sp)
517
518 !TODO: can we exploid more symmetry?
519 CALL dbcsr_multiply("N", "N", 3.0_dp, matrix_k, matrix_ps, &
520 0.0_dp, matrix_sum)
521 CALL dbcsr_multiply("N", "N", 3.0_dp, matrix_sp, matrix_k, &
522 1.0_dp, matrix_sum)
523
524 !matrix_tmp = KPS
525 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_k, matrix_ps, &
526 0.0_dp, matrix_tmp)
527 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_tmp, matrix_ps, &
528 1.0_dp, matrix_sum)
529 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, &
530 1.0_dp, matrix_sum)
531
532 !matrix_tmp = SPK
533 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sp, matrix_k, &
534 0.0_dp, matrix_tmp)
535 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, &
536 1.0_dp, matrix_sum)
537
538 ! overwrite matrix_k
539 CALL dbcsr_copy(matrix_k, matrix_sum, name="K from reverse McWeeny")
540
541 ! clean up
542 CALL dbcsr_release(matrix_sum)
543 CALL dbcsr_release(matrix_tmp)
544 CALL dbcsr_release(matrix_ps)
545 CALL dbcsr_release(matrix_sp)
546 CALL timestop(handle)
547 END SUBROUTINE reverse_mcweeny_step
548
549END MODULE admm_dm_methods
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
iterate_matrix::purify_mcweeny
Definition iterate_matrix.F:58
admm_dm_methods
Contains ADMM methods which only require the density matrix.
Definition admm_dm_methods.F:14
admm_dm_methods::admm_dm_merge_ks_matrix
subroutine, public admm_dm_merge_ks_matrix(qs_env)
Entry methods: Merges auxiliary Kohn-Sham matrix into primary one.
Definition admm_dm_methods.F:93
admm_dm_methods::admm_dm_calc_rho_aux
subroutine, public admm_dm_calc_rho_aux(qs_env)
Entry methods: Calculates auxiliary density matrix from primary one.
Definition admm_dm_methods.F:57
admm_dm_types
Types and set/get functions for auxiliary density matrix methods.
Definition admm_dm_types.F:14
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition admm_types.F:593
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_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1165
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
Definition cp_dbcsr_api.F:956
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:930
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1027
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:679
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1073
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:989
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_frobenius_norm
real(dp) function, public dbcsr_frobenius_norm(matrix)
Compute the frobenius norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:127
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
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:567
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_admm_blocked_projection
integer, parameter, public do_admm_blocked_projection
Definition input_constants.F:863
input_constants::do_admm_basis_projection
integer, parameter, public do_admm_basis_projection
Definition input_constants.F:863
iterate_matrix
Routines useful for iterative matrix calculations.
Definition iterate_matrix.F:13
iterate_matrix::invert_hotelling
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Definition iterate_matrix.F:473
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
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:1711
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, 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, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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_set
subroutine, public qs_rho_set(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)
...
Definition qs_rho_types.F:308
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
task_list_types
types for task lists
Definition task_list_types.F:14
admm_dm_types::admm_dm_type
Definition admm_dm_types.F:41
admm_dm_types::mcweeny_history_type
Definition admm_dm_types.F:31
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:186
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
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:217
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:136
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
task_list_types::task_list_type
Definition task_list_types.F:58