(git:b279b6b)
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-2024 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 ! **************************************************************************************************
14 MODULE 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_operations, ONLY: dbcsr_deallocate_matrix_set
20  USE cp_log_handling, ONLY: cp_logger_get_default_unit_nr
21  USE dbcsr_api, ONLY: &
22  dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_frobenius_norm, dbcsr_get_block_p, &
23  dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
24  dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
25  dbcsr_scale, dbcsr_set, dbcsr_type
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 
49 CONTAINS
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 :: blk, 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, blk)
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 :: blk, 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, blk)
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 
549 END MODULE admm_dm_methods
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
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_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:857
input_constants::do_admm_basis_projection
integer, parameter, public do_admm_basis_projection
Definition: input_constants.F:857
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:472
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:1710
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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
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