(git:98357aa)
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-2026 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) THEN
80 CALL purify_mcweeny(qs_env)
81 END IF
82
83 CALL update_rho_aux(qs_env)
84
85 CALL timestop(handle)
86 END SUBROUTINE admm_dm_calc_rho_aux
87
88! **************************************************************************************************
89!> \brief Entry methods: Merges auxiliary Kohn-Sham matrix into primary one.
90!> \param qs_env ...
91!> \author Ole Schuett
92! **************************************************************************************************
93 SUBROUTINE admm_dm_merge_ks_matrix(qs_env)
94 TYPE(qs_environment_type), POINTER :: qs_env
95
96 CHARACTER(LEN=*), PARAMETER :: routinen = 'admm_dm_merge_ks_matrix'
97
98 INTEGER :: handle
99 TYPE(admm_dm_type), POINTER :: admm_dm
100 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge
101
102 CALL timeset(routinen, handle)
103 NULLIFY (admm_dm, matrix_ks_merge)
104
105 CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
106
107 IF (admm_dm%purify) THEN
108 CALL revert_purify_mcweeny(qs_env, matrix_ks_merge)
109 ELSE
110 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit=matrix_ks_merge)
111 END IF
112
113 SELECT CASE (admm_dm%method)
114 CASE (do_admm_basis_projection)
115 CALL merge_dm_projection(qs_env, matrix_ks_merge)
116
117 CASE (do_admm_blocked_projection)
118 CALL merge_dm_blocked(qs_env, matrix_ks_merge)
119
120 CASE DEFAULT
121 cpabort("admm_dm_merge_ks_matrix: unknown method")
122 END SELECT
123
124 IF (admm_dm%purify) THEN
125 CALL dbcsr_deallocate_matrix_set(matrix_ks_merge)
126 END IF
127
128 CALL timestop(handle)
129
130 END SUBROUTINE admm_dm_merge_ks_matrix
131
132! **************************************************************************************************
133!> \brief Calculates auxiliary density matrix via basis projection.
134!> \param qs_env ...
135!> \author Ole Schuett
136! **************************************************************************************************
137 SUBROUTINE map_dm_projection(qs_env)
138 TYPE(qs_environment_type), POINTER :: qs_env
139
140 INTEGER :: ispin
141 LOGICAL :: s_mstruct_changed
142 REAL(kind=dp) :: threshold
143 TYPE(admm_dm_type), POINTER :: admm_dm
144 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux, matrix_s_mixed, rho_ao, &
145 rho_ao_aux
146 TYPE(dbcsr_type) :: matrix_s_aux_inv, matrix_tmp
147 TYPE(dft_control_type), POINTER :: dft_control
148 TYPE(qs_rho_type), POINTER :: rho, rho_aux
149
150 NULLIFY (dft_control, admm_dm, matrix_s_aux, matrix_s_mixed, rho, rho_aux)
151 NULLIFY (rho_ao, rho_ao_aux)
152
153 CALL get_qs_env(qs_env, dft_control=dft_control, s_mstruct_changed=s_mstruct_changed, rho=rho)
154 CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux, rho_aux_fit=rho_aux, &
155 matrix_s_aux_fit_vs_orb=matrix_s_mixed, admm_dm=admm_dm)
156
157 CALL qs_rho_get(rho, rho_ao=rho_ao)
158 CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux)
159
160 IF (s_mstruct_changed) THEN
161 ! Calculate A = S_aux^(-1) * S_mixed
162 CALL dbcsr_create(matrix_s_aux_inv, template=matrix_s_aux(1)%matrix, matrix_type="N")
163 threshold = max(admm_dm%eps_filter, 1.0e-12_dp)
164 CALL invert_hotelling(matrix_s_aux_inv, matrix_s_aux(1)%matrix, threshold)
165
166 IF (.NOT. ASSOCIATED(admm_dm%matrix_A)) THEN
167 ALLOCATE (admm_dm%matrix_A)
168 CALL dbcsr_create(admm_dm%matrix_A, template=matrix_s_mixed(1)%matrix, matrix_type="N")
169 END IF
170 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_aux_inv, matrix_s_mixed(1)%matrix, &
171 0.0_dp, admm_dm%matrix_A)
172 CALL dbcsr_release(matrix_s_aux_inv)
173 END IF
174
175 ! Calculate P_aux = A * P * A^T
176 CALL dbcsr_create(matrix_tmp, template=admm_dm%matrix_A)
177 DO ispin = 1, dft_control%nspins
178 CALL dbcsr_multiply("N", "N", 1.0_dp, admm_dm%matrix_A, rho_ao(ispin)%matrix, &
179 0.0_dp, matrix_tmp)
180 CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp, admm_dm%matrix_A, &
181 0.0_dp, rho_ao_aux(ispin)%matrix)
182 END DO
183 CALL dbcsr_release(matrix_tmp)
184
185 END SUBROUTINE map_dm_projection
186
187! **************************************************************************************************
188!> \brief Calculates auxiliary density matrix via blocking.
189!> \param qs_env ...
190!> \author Ole Schuett
191! **************************************************************************************************
192 SUBROUTINE map_dm_blocked(qs_env)
193 TYPE(qs_environment_type), POINTER :: qs_env
194
195 INTEGER :: iatom, ispin, jatom
196 LOGICAL :: found
197 REAL(dp), DIMENSION(:, :), POINTER :: sparse_block, sparse_block_aux
198 TYPE(admm_dm_type), POINTER :: admm_dm
199 TYPE(dbcsr_iterator_type) :: iter
200 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_ao_aux
201 TYPE(dft_control_type), POINTER :: dft_control
202 TYPE(qs_rho_type), POINTER :: rho, rho_aux
203
204 NULLIFY (dft_control, admm_dm, rho, rho_aux, rho_ao, rho_ao_aux)
205
206 CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
207 CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux, admm_dm=admm_dm)
208
209 CALL qs_rho_get(rho, rho_ao=rho_ao)
210 CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux)
211
212 ! ** set blocked density matrix to 0
213 DO ispin = 1, dft_control%nspins
214 CALL dbcsr_set(rho_ao_aux(ispin)%matrix, 0.0_dp)
215 ! ** now loop through the list and copy corresponding blocks
216 CALL dbcsr_iterator_start(iter, rho_ao(ispin)%matrix)
217 DO WHILE (dbcsr_iterator_blocks_left(iter))
218 CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block)
219 IF (admm_dm%block_map(iatom, jatom) == 1) THEN
220 CALL dbcsr_get_block_p(rho_ao_aux(ispin)%matrix, &
221 row=iatom, col=jatom, block=sparse_block_aux, found=found)
222 IF (found) THEN
223 sparse_block_aux = sparse_block
224 END IF
225 END IF
226 END DO
227 CALL dbcsr_iterator_stop(iter)
228 END DO
229
230 END SUBROUTINE map_dm_blocked
231
232! **************************************************************************************************
233!> \brief Call calculate_rho_elec() for auxiliary density
234!> \param qs_env ...
235! **************************************************************************************************
236 SUBROUTINE update_rho_aux(qs_env)
237 TYPE(qs_environment_type), POINTER :: qs_env
238
239 INTEGER :: ispin
240 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r_aux
241 TYPE(admm_dm_type), POINTER :: admm_dm
242 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_aux
243 TYPE(dft_control_type), POINTER :: dft_control
244 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
245 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_aux
246 TYPE(qs_ks_env_type), POINTER :: ks_env
247 TYPE(qs_rho_type), POINTER :: rho_aux
248 TYPE(task_list_type), POINTER :: task_list_aux_fit
249
250 NULLIFY (dft_control, admm_dm, rho_aux, rho_ao_aux, rho_r_aux, rho_g_aux, tot_rho_r_aux, &
251 task_list_aux_fit, ks_env)
252
253 CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control)
254 CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list_aux_fit, rho_aux_fit=rho_aux, &
255 admm_dm=admm_dm)
256
257 CALL qs_rho_get(rho_aux, &
258 rho_ao=rho_ao_aux, &
259 rho_r=rho_r_aux, &
260 rho_g=rho_g_aux, &
261 tot_rho_r=tot_rho_r_aux)
262
263 DO ispin = 1, dft_control%nspins
264 CALL calculate_rho_elec(ks_env=ks_env, &
265 matrix_p=rho_ao_aux(ispin)%matrix, &
266 rho=rho_r_aux(ispin), &
267 rho_gspace=rho_g_aux(ispin), &
268 total_rho=tot_rho_r_aux(ispin), &
269 soft_valid=.false., &
270 basis_type="AUX_FIT", &
271 task_list_external=task_list_aux_fit)
272 END DO
273
274 CALL qs_rho_set(rho_aux, rho_r_valid=.true., rho_g_valid=.true.)
275
276 END SUBROUTINE update_rho_aux
277
278! **************************************************************************************************
279!> \brief Merges auxiliary Kohn-Sham matrix via basis projection.
280!> \param qs_env ...
281!> \param matrix_ks_merge Input: The KS matrix to be merged
282!> \author Ole Schuett
283! **************************************************************************************************
284 SUBROUTINE merge_dm_projection(qs_env, matrix_ks_merge)
285 TYPE(qs_environment_type), POINTER :: qs_env
286 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge
287
288 INTEGER :: ispin
289 TYPE(admm_dm_type), POINTER :: admm_dm
290 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
291 TYPE(dbcsr_type) :: matrix_tmp
292 TYPE(dft_control_type), POINTER :: dft_control
293
294 NULLIFY (admm_dm, dft_control, matrix_ks)
295
296 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks)
297 CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
298
299 ! Calculate K += A^T * K_aux * A
300 CALL dbcsr_create(matrix_tmp, template=admm_dm%matrix_A, matrix_type="N")
301
302 DO ispin = 1, dft_control%nspins
303 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ks_merge(ispin)%matrix, admm_dm%matrix_A, &
304 0.0_dp, matrix_tmp)
305 CALL dbcsr_multiply("T", "N", 1.0_dp, admm_dm%matrix_A, matrix_tmp, &
306 1.0_dp, matrix_ks(ispin)%matrix)
307 END DO
308
309 CALL dbcsr_release(matrix_tmp)
310
311 END SUBROUTINE merge_dm_projection
312
313! **************************************************************************************************
314!> \brief Merges auxiliary Kohn-Sham matrix via blocking.
315!> \param qs_env ...
316!> \param matrix_ks_merge Input: The KS matrix to be merged
317!> \author Ole Schuett
318! **************************************************************************************************
319 SUBROUTINE merge_dm_blocked(qs_env, matrix_ks_merge)
320 TYPE(qs_environment_type), POINTER :: qs_env
321 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge
322
323 INTEGER :: iatom, ispin, jatom
324 REAL(dp), DIMENSION(:, :), POINTER :: sparse_block
325 TYPE(admm_dm_type), POINTER :: admm_dm
326 TYPE(dbcsr_iterator_type) :: iter
327 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
328 TYPE(dft_control_type), POINTER :: dft_control
329
330 NULLIFY (admm_dm, dft_control, matrix_ks)
331
332 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks)
333 CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
334
335 DO ispin = 1, dft_control%nspins
336 CALL dbcsr_iterator_start(iter, matrix_ks_merge(ispin)%matrix)
337 DO WHILE (dbcsr_iterator_blocks_left(iter))
338 CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block)
339 IF (admm_dm%block_map(iatom, jatom) == 0) THEN
340 sparse_block = 0.0_dp
341 END IF
342 END DO
343 CALL dbcsr_iterator_stop(iter)
344 CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_ks_merge(ispin)%matrix, 1.0_dp, 1.0_dp)
345 END DO
346
347 END SUBROUTINE merge_dm_blocked
348
349! **************************************************************************************************
350!> \brief Apply McWeeny purification to auxiliary density matrix
351!> \param qs_env ...
352!> \author Ole Schuett
353! **************************************************************************************************
354 SUBROUTINE purify_mcweeny(qs_env)
355 TYPE(qs_environment_type), POINTER :: qs_env
356
357 CHARACTER(LEN=*), PARAMETER :: routinen = 'purify_mcweeny'
358
359 INTEGER :: handle, ispin, istep, nspins, unit_nr
360 REAL(kind=dp) :: frob_norm
361 TYPE(admm_dm_type), POINTER :: admm_dm
362 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, rho_ao_aux
363 TYPE(dbcsr_type) :: matrix_ps, matrix_psp, matrix_test
364 TYPE(dbcsr_type), POINTER :: matrix_p, matrix_s
365 TYPE(dft_control_type), POINTER :: dft_control
366 TYPE(mcweeny_history_type), POINTER :: history, new_hist_entry
367 TYPE(qs_rho_type), POINTER :: rho_aux_fit
368
369 CALL timeset(routinen, handle)
370 NULLIFY (dft_control, admm_dm, matrix_s_aux_fit, rho_aux_fit, new_hist_entry, &
371 matrix_p, matrix_s, rho_ao_aux)
372
373 unit_nr = cp_logger_get_default_unit_nr()
374 CALL get_qs_env(qs_env, dft_control=dft_control)
375 CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux_fit, &
376 rho_aux_fit=rho_aux_fit, admm_dm=admm_dm)
377
378 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux)
379
380 matrix_p => rho_ao_aux(1)%matrix
381 CALL dbcsr_create(matrix_ps, template=matrix_p, matrix_type="N")
382 CALL dbcsr_create(matrix_psp, template=matrix_p, matrix_type="S")
383 CALL dbcsr_create(matrix_test, template=matrix_p, matrix_type="S")
384
385 nspins = dft_control%nspins
386 DO ispin = 1, nspins
387 matrix_p => rho_ao_aux(ispin)%matrix
388 matrix_s => matrix_s_aux_fit(1)%matrix
389 history => admm_dm%mcweeny_history(ispin)%p
390 IF (ASSOCIATED(history)) cpabort("purify_dm_mcweeny: history already associated")
391 IF (nspins == 1) CALL dbcsr_scale(matrix_p, 0.5_dp)
392
393 DO istep = 1, admm_dm%mcweeny_max_steps
394 ! allocate new element in linked list
395 ALLOCATE (new_hist_entry)
396 new_hist_entry%next => history
397 history => new_hist_entry
398 history%count = istep
399 NULLIFY (new_hist_entry)
400 CALL dbcsr_create(history%m, template=matrix_p, matrix_type="N")
401 CALL dbcsr_copy(history%m, matrix_p, name="P from McWeeny")
402
403 ! calc PS and PSP
404 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, matrix_s, &
405 0.0_dp, matrix_ps)
406
407 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ps, matrix_p, &
408 0.0_dp, matrix_psp)
409
410 !test convergence
411 CALL dbcsr_copy(matrix_test, matrix_psp)
412 CALL dbcsr_add(matrix_test, matrix_p, 1.0_dp, -1.0_dp)
413 frob_norm = dbcsr_frobenius_norm(matrix_test)
414 IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,i5,a,f16.8)') "McWeeny-Step", istep, &
415 ": Deviation of idempotency", frob_norm
416 IF (frob_norm < 1000_dp*admm_dm%eps_filter .AND. istep > 1) EXIT
417
418 ! build next P matrix
419 CALL dbcsr_copy(matrix_p, matrix_psp, name="P from McWeeny")
420 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_ps, matrix_psp, &
421 3.0_dp, matrix_p)
422 END DO
423 admm_dm%mcweeny_history(ispin)%p => history
424 IF (nspins == 1) CALL dbcsr_scale(matrix_p, 2.0_dp)
425 END DO
426
427 ! clean up
428 CALL dbcsr_release(matrix_ps)
429 CALL dbcsr_release(matrix_psp)
430 CALL dbcsr_release(matrix_test)
431 CALL timestop(handle)
432 END SUBROUTINE purify_mcweeny
433
434! **************************************************************************************************
435!> \brief Prepare auxiliary KS-matrix for merge using reverse McWeeny
436!> \param qs_env ...
437!> \param matrix_ks_merge Output: The KS matrix for the merge
438!> \author Ole Schuett
439! **************************************************************************************************
440 SUBROUTINE revert_purify_mcweeny(qs_env, matrix_ks_merge)
441 TYPE(qs_environment_type), POINTER :: qs_env
442 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge
443
444 CHARACTER(LEN=*), PARAMETER :: routinen = 'revert_purify_mcweeny'
445
446 INTEGER :: handle, ispin, nspins, unit_nr
447 TYPE(admm_dm_type), POINTER :: admm_dm
448 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
449 matrix_s_aux_fit, &
450 matrix_s_aux_fit_vs_orb
451 TYPE(dbcsr_type), POINTER :: matrix_k
452 TYPE(dft_control_type), POINTER :: dft_control
453 TYPE(mcweeny_history_type), POINTER :: history_curr, history_next
454
455 CALL timeset(routinen, handle)
456 unit_nr = cp_logger_get_default_unit_nr()
457 NULLIFY (admm_dm, dft_control, matrix_ks, matrix_ks_aux_fit, &
458 matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, &
459 history_next, history_curr, matrix_k)
460
461 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks)
462 CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux_fit, admm_dm=admm_dm, &
463 matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, matrix_ks_aux_fit=matrix_ks_aux_fit)
464
465 nspins = dft_control%nspins
466 ALLOCATE (matrix_ks_merge(nspins))
467
468 DO ispin = 1, nspins
469 ALLOCATE (matrix_ks_merge(ispin)%matrix)
470 matrix_k => matrix_ks_merge(ispin)%matrix
471 CALL dbcsr_copy(matrix_k, matrix_ks_aux_fit(ispin)%matrix, name="K")
472 history_curr => admm_dm%mcweeny_history(ispin)%p
473 NULLIFY (admm_dm%mcweeny_history(ispin)%p)
474
475 ! reverse McWeeny iteration
476 DO WHILE (ASSOCIATED(history_curr))
477 IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,i5)') "Reverse McWeeny-Step ", history_curr%count
478 CALL reverse_mcweeny_step(matrix_k=matrix_k, &
479 matrix_s=matrix_s_aux_fit(1)%matrix, &
480 matrix_p=history_curr%m)
481 CALL dbcsr_release(history_curr%m)
482 history_next => history_curr%next
483 DEALLOCATE (history_curr)
484 history_curr => history_next
485 NULLIFY (history_next)
486 END DO
487
488 END DO
489
490 ! clean up
491 CALL timestop(handle)
492
493 END SUBROUTINE revert_purify_mcweeny
494
495! **************************************************************************************************
496!> \brief Multiply matrix_k with partial derivative of McWeeny by reversing it.
497!> \param matrix_k ...
498!> \param matrix_s ...
499!> \param matrix_p ...
500!> \author Ole Schuett
501! **************************************************************************************************
502 SUBROUTINE reverse_mcweeny_step(matrix_k, matrix_s, matrix_p)
503 TYPE(dbcsr_type) :: matrix_k, matrix_s, matrix_p
504
505 CHARACTER(LEN=*), PARAMETER :: routinen = 'reverse_mcweeny_step'
506
507 INTEGER :: handle
508 TYPE(dbcsr_type) :: matrix_ps, matrix_sp, matrix_sum, &
509 matrix_tmp
510
511 CALL timeset(routinen, handle)
512 CALL dbcsr_create(matrix_ps, template=matrix_p, matrix_type="N")
513 CALL dbcsr_create(matrix_sp, template=matrix_p, matrix_type="N")
514 CALL dbcsr_create(matrix_tmp, template=matrix_p, matrix_type="N")
515 CALL dbcsr_create(matrix_sum, template=matrix_p, matrix_type="N")
516
517 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, matrix_s, &
518 0.0_dp, matrix_ps)
519 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, matrix_p, &
520 0.0_dp, matrix_sp)
521
522 !TODO: can we exploid more symmetry?
523 CALL dbcsr_multiply("N", "N", 3.0_dp, matrix_k, matrix_ps, &
524 0.0_dp, matrix_sum)
525 CALL dbcsr_multiply("N", "N", 3.0_dp, matrix_sp, matrix_k, &
526 1.0_dp, matrix_sum)
527
528 !matrix_tmp = KPS
529 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_k, matrix_ps, &
530 0.0_dp, matrix_tmp)
531 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_tmp, matrix_ps, &
532 1.0_dp, matrix_sum)
533 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, &
534 1.0_dp, matrix_sum)
535
536 !matrix_tmp = SPK
537 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sp, matrix_k, &
538 0.0_dp, matrix_tmp)
539 CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, &
540 1.0_dp, matrix_sum)
541
542 ! overwrite matrix_k
543 CALL dbcsr_copy(matrix_k, matrix_sum, name="K from reverse McWeeny")
544
545 ! clean up
546 CALL dbcsr_release(matrix_sum)
547 CALL dbcsr_release(matrix_tmp)
548 CALL dbcsr_release(matrix_ps)
549 CALL dbcsr_release(matrix_sp)
550 CALL timestop(handle)
551 END SUBROUTINE reverse_mcweeny_step
552
553END MODULE admm_dm_methods
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
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:94
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:599
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:1188
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
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_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
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:1096
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, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
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_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
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:576
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::do_admm_blocked_projection
integer, parameter, public do_admm_blocked_projection
Definition input_constants.F:897
input_constants::do_admm_basis_projection
integer, parameter, public do_admm_basis_projection
Definition input_constants.F:897
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:1720
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:534
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:757
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
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
task_list_types::task_list_type
Definition task_list_types.F:58