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 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_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 :
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 214 : 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 214 : NULLIFY (admm_dm)
65 214 : CALL timeset(routineN, handle)
66 214 : CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
67 :
68 368 : SELECT CASE (admm_dm%method)
69 : CASE (do_admm_basis_projection)
70 154 : CALL map_dm_projection(qs_env)
71 :
72 : CASE (do_admm_blocked_projection)
73 60 : CALL map_dm_blocked(qs_env)
74 :
75 : CASE DEFAULT
76 214 : CPABORT("admm_dm_calc_rho_aux: unknown method")
77 : END SELECT
78 :
79 214 : IF (admm_dm%purify) THEN
80 38 : CALL purify_mcweeny(qs_env)
81 : END IF
82 :
83 214 : CALL update_rho_aux(qs_env)
84 :
85 214 : CALL timestop(handle)
86 214 : 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 214 : 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 214 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge
101 :
102 214 : CALL timeset(routineN, handle)
103 214 : NULLIFY (admm_dm, matrix_ks_merge)
104 :
105 214 : CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
106 :
107 214 : IF (admm_dm%purify) THEN
108 38 : CALL revert_purify_mcweeny(qs_env, matrix_ks_merge)
109 : ELSE
110 176 : CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit=matrix_ks_merge)
111 : END IF
112 :
113 368 : SELECT CASE (admm_dm%method)
114 : CASE (do_admm_basis_projection)
115 154 : CALL merge_dm_projection(qs_env, matrix_ks_merge)
116 :
117 : CASE (do_admm_blocked_projection)
118 60 : CALL merge_dm_blocked(qs_env, matrix_ks_merge)
119 :
120 : CASE DEFAULT
121 214 : CPABORT("admm_dm_merge_ks_matrix: unknown method")
122 : END SELECT
123 :
124 214 : IF (admm_dm%purify) THEN
125 38 : CALL dbcsr_deallocate_matrix_set(matrix_ks_merge)
126 : END IF
127 :
128 214 : CALL timestop(handle)
129 :
130 214 : 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 154 : 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 154 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux, matrix_s_mixed, rho_ao, &
145 154 : 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 154 : NULLIFY (dft_control, admm_dm, matrix_s_aux, matrix_s_mixed, rho, rho_aux)
151 154 : NULLIFY (rho_ao, rho_ao_aux)
152 :
153 154 : 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 154 : matrix_s_aux_fit_vs_orb=matrix_s_mixed, admm_dm=admm_dm)
156 :
157 154 : CALL qs_rho_get(rho, rho_ao=rho_ao)
158 154 : CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux)
159 :
160 154 : IF (s_mstruct_changed) THEN
161 : ! Calculate A = S_aux^(-1) * S_mixed
162 10 : CALL dbcsr_create(matrix_s_aux_inv, template=matrix_s_aux(1)%matrix, matrix_type="N")
163 10 : threshold = MAX(admm_dm%eps_filter, 1.0e-12_dp)
164 10 : CALL invert_Hotelling(matrix_s_aux_inv, matrix_s_aux(1)%matrix, threshold)
165 :
166 10 : IF (.NOT. ASSOCIATED(admm_dm%matrix_A)) THEN
167 10 : ALLOCATE (admm_dm%matrix_A)
168 10 : 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 10 : 0.0_dp, admm_dm%matrix_A)
172 10 : CALL dbcsr_release(matrix_s_aux_inv)
173 : END IF
174 :
175 : ! Calculate P_aux = A * P * A^T
176 154 : CALL dbcsr_create(matrix_tmp, template=admm_dm%matrix_A)
177 404 : DO ispin = 1, dft_control%nspins
178 : CALL dbcsr_multiply("N", "N", 1.0_dp, admm_dm%matrix_A, rho_ao(ispin)%matrix, &
179 250 : 0.0_dp, matrix_tmp)
180 : CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp, admm_dm%matrix_A, &
181 404 : 0.0_dp, rho_ao_aux(ispin)%matrix)
182 : END DO
183 154 : CALL dbcsr_release(matrix_tmp)
184 :
185 154 : 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 60 : SUBROUTINE map_dm_blocked(qs_env)
193 : TYPE(qs_environment_type), POINTER :: qs_env
194 :
195 : INTEGER :: iatom, ispin, jatom
196 : LOGICAL :: found
197 60 : REAL(dp), DIMENSION(:, :), POINTER :: sparse_block, sparse_block_aux
198 : TYPE(admm_dm_type), POINTER :: admm_dm
199 : TYPE(dbcsr_iterator_type) :: iter
200 60 : 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 60 : NULLIFY (dft_control, admm_dm, rho, rho_aux, rho_ao, rho_ao_aux)
205 :
206 60 : CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
207 60 : CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux, admm_dm=admm_dm)
208 :
209 60 : CALL qs_rho_get(rho, rho_ao=rho_ao)
210 60 : CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux)
211 :
212 : ! ** set blocked density matrix to 0
213 142 : DO ispin = 1, dft_control%nspins
214 82 : CALL dbcsr_set(rho_ao_aux(ispin)%matrix, 0.0_dp)
215 : ! ** now loop through the list and copy corresponding blocks
216 82 : CALL dbcsr_iterator_start(iter, rho_ao(ispin)%matrix)
217 471 : DO WHILE (dbcsr_iterator_blocks_left(iter))
218 389 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block)
219 471 : IF (admm_dm%block_map(iatom, jatom) == 1) THEN
220 : CALL dbcsr_get_block_p(rho_ao_aux(ispin)%matrix, &
221 272 : row=iatom, col=jatom, BLOCK=sparse_block_aux, found=found)
222 272 : IF (found) THEN
223 3456 : sparse_block_aux = sparse_block
224 : END IF
225 : END IF
226 : END DO
227 224 : CALL dbcsr_iterator_stop(iter)
228 : END DO
229 :
230 60 : END SUBROUTINE map_dm_blocked
231 :
232 : ! **************************************************************************************************
233 : !> \brief Call calculate_rho_elec() for auxiliary density
234 : !> \param qs_env ...
235 : ! **************************************************************************************************
236 214 : SUBROUTINE update_rho_aux(qs_env)
237 : TYPE(qs_environment_type), POINTER :: qs_env
238 :
239 : INTEGER :: ispin
240 214 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r_aux
241 : TYPE(admm_dm_type), POINTER :: admm_dm
242 214 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_aux
243 : TYPE(dft_control_type), POINTER :: dft_control
244 214 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
245 214 : 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 214 : NULLIFY (dft_control, admm_dm, rho_aux, rho_ao_aux, rho_r_aux, rho_g_aux, tot_rho_r_aux, &
251 214 : task_list_aux_fit, ks_env)
252 :
253 214 : 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 214 : 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 214 : tot_rho_r=tot_rho_r_aux)
262 :
263 546 : 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 546 : task_list_external=task_list_aux_fit)
272 : END DO
273 :
274 214 : CALL qs_rho_set(rho_aux, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
275 :
276 214 : 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 154 : 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 154 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
291 : TYPE(dbcsr_type) :: matrix_tmp
292 : TYPE(dft_control_type), POINTER :: dft_control
293 :
294 154 : NULLIFY (admm_dm, dft_control, matrix_ks)
295 :
296 154 : CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks)
297 154 : CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
298 :
299 : ! Calculate K += A^T * K_aux * A
300 154 : CALL dbcsr_create(matrix_tmp, template=admm_dm%matrix_A, matrix_type="N")
301 :
302 404 : 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 250 : 0.0_dp, matrix_tmp)
305 : CALL dbcsr_multiply("T", "N", 1.0_dp, admm_dm%matrix_A, matrix_tmp, &
306 404 : 1.0_dp, matrix_ks(ispin)%matrix)
307 : END DO
308 :
309 154 : CALL dbcsr_release(matrix_tmp)
310 :
311 154 : 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 60 : 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 60 : REAL(dp), DIMENSION(:, :), POINTER :: sparse_block
325 : TYPE(admm_dm_type), POINTER :: admm_dm
326 : TYPE(dbcsr_iterator_type) :: iter
327 60 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
328 : TYPE(dft_control_type), POINTER :: dft_control
329 :
330 60 : NULLIFY (admm_dm, dft_control, matrix_ks)
331 :
332 60 : CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks)
333 60 : CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
334 :
335 142 : DO ispin = 1, dft_control%nspins
336 82 : CALL dbcsr_iterator_start(iter, matrix_ks_merge(ispin)%matrix)
337 471 : DO WHILE (dbcsr_iterator_blocks_left(iter))
338 389 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block)
339 471 : IF (admm_dm%block_map(iatom, jatom) == 0) THEN
340 978 : sparse_block = 0.0_dp
341 : END IF
342 : END DO
343 82 : CALL dbcsr_iterator_stop(iter)
344 224 : CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_ks_merge(ispin)%matrix, 1.0_dp, 1.0_dp)
345 : END DO
346 :
347 60 : 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 38 : 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 38 : 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 38 : CALL timeset(routineN, handle)
370 38 : NULLIFY (dft_control, admm_dm, matrix_s_aux_fit, rho_aux_fit, new_hist_entry, &
371 38 : matrix_p, matrix_s, rho_ao_aux)
372 :
373 38 : unit_nr = cp_logger_get_default_unit_nr()
374 38 : 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 38 : rho_aux_fit=rho_aux_fit, admm_dm=admm_dm)
377 :
378 38 : CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux)
379 :
380 38 : matrix_p => rho_ao_aux(1)%matrix
381 38 : CALL dbcsr_create(matrix_PS, template=matrix_p, matrix_type="N")
382 38 : CALL dbcsr_create(matrix_PSP, template=matrix_p, matrix_type="S")
383 38 : CALL dbcsr_create(matrix_test, template=matrix_p, matrix_type="S")
384 :
385 38 : nspins = dft_control%nspins
386 114 : DO ispin = 1, nspins
387 76 : matrix_p => rho_ao_aux(ispin)%matrix
388 76 : matrix_s => matrix_s_aux_fit(1)%matrix
389 76 : history => admm_dm%mcweeny_history(ispin)%p
390 76 : IF (ASSOCIATED(history)) CPABORT("purify_dm_mcweeny: history already associated")
391 76 : IF (nspins == 1) CALL dbcsr_scale(matrix_p, 0.5_dp)
392 :
393 336 : DO istep = 1, admm_dm%mcweeny_max_steps
394 : ! allocate new element in linked list
395 336 : ALLOCATE (new_hist_entry)
396 336 : new_hist_entry%next => history
397 336 : history => new_hist_entry
398 336 : history%count = istep
399 336 : NULLIFY (new_hist_entry)
400 336 : CALL dbcsr_create(history%m, template=matrix_p, matrix_type="N")
401 336 : 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 336 : 0.0_dp, matrix_ps)
406 :
407 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ps, matrix_p, &
408 336 : 0.0_dp, matrix_psp)
409 :
410 : !test convergence
411 336 : CALL dbcsr_copy(matrix_test, matrix_psp)
412 336 : CALL dbcsr_add(matrix_test, matrix_p, 1.0_dp, -1.0_dp)
413 336 : frob_norm = dbcsr_frobenius_norm(matrix_test)
414 672 : IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,i5,a,f16.8)') "McWeeny-Step", istep, &
415 672 : ": Deviation of idempotency", frob_norm
416 336 : IF (frob_norm < 1000_dp*admm_dm%eps_filter .AND. istep > 1) EXIT
417 :
418 : ! build next P matrix
419 260 : 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 336 : 3.0_dp, matrix_p)
422 : END DO
423 76 : admm_dm%mcweeny_history(ispin)%p => history
424 114 : IF (nspins == 1) CALL dbcsr_scale(matrix_p, 2.0_dp)
425 : END DO
426 :
427 : ! clean up
428 38 : CALL dbcsr_release(matrix_PS)
429 38 : CALL dbcsr_release(matrix_PSP)
430 38 : CALL dbcsr_release(matrix_test)
431 38 : CALL timestop(handle)
432 38 : 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 38 : 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 38 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
449 38 : matrix_s_aux_fit, &
450 38 : 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 38 : CALL timeset(routineN, handle)
456 38 : unit_nr = cp_logger_get_default_unit_nr()
457 38 : NULLIFY (admm_dm, dft_control, matrix_ks, matrix_ks_aux_fit, &
458 38 : matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, &
459 38 : history_next, history_curr, matrix_k)
460 :
461 38 : 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 38 : matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, matrix_ks_aux_fit=matrix_ks_aux_fit)
464 :
465 38 : nspins = dft_control%nspins
466 190 : ALLOCATE (matrix_ks_merge(nspins))
467 :
468 114 : DO ispin = 1, nspins
469 76 : ALLOCATE (matrix_ks_merge(ispin)%matrix)
470 76 : matrix_k => matrix_ks_merge(ispin)%matrix
471 76 : CALL dbcsr_copy(matrix_k, matrix_ks_aux_fit(ispin)%matrix, name="K")
472 76 : history_curr => admm_dm%mcweeny_history(ispin)%p
473 76 : NULLIFY (admm_dm%mcweeny_history(ispin)%p)
474 :
475 : ! reverse McWeeny iteration
476 450 : DO WHILE (ASSOCIATED(history_curr))
477 336 : 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 336 : matrix_p=history_curr%m)
481 336 : CALL dbcsr_release(history_curr%m)
482 336 : history_next => history_curr%next
483 336 : DEALLOCATE (history_curr)
484 336 : history_curr => history_next
485 336 : NULLIFY (history_next)
486 : END DO
487 :
488 : END DO
489 :
490 : ! clean up
491 38 : CALL timestop(handle)
492 :
493 38 : 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 336 : 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 336 : CALL timeset(routineN, handle)
512 336 : CALL dbcsr_create(matrix_ps, template=matrix_p, matrix_type="N")
513 336 : CALL dbcsr_create(matrix_sp, template=matrix_p, matrix_type="N")
514 336 : CALL dbcsr_create(matrix_tmp, template=matrix_p, matrix_type="N")
515 336 : 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 336 : 0.0_dp, matrix_ps)
519 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, matrix_p, &
520 336 : 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 336 : 0.0_dp, matrix_sum)
525 : CALL dbcsr_multiply("N", "N", 3.0_dp, matrix_sp, matrix_k, &
526 336 : 1.0_dp, matrix_sum)
527 :
528 : !matrix_tmp = KPS
529 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_k, matrix_ps, &
530 336 : 0.0_dp, matrix_tmp)
531 : CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_tmp, matrix_ps, &
532 336 : 1.0_dp, matrix_sum)
533 : CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, &
534 336 : 1.0_dp, matrix_sum)
535 :
536 : !matrix_tmp = SPK
537 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sp, matrix_k, &
538 336 : 0.0_dp, matrix_tmp)
539 : CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, &
540 336 : 1.0_dp, matrix_sum)
541 :
542 : ! overwrite matrix_k
543 336 : CALL dbcsr_copy(matrix_k, matrix_sum, name="K from reverse McWeeny")
544 :
545 : ! clean up
546 336 : CALL dbcsr_release(matrix_sum)
547 336 : CALL dbcsr_release(matrix_tmp)
548 336 : CALL dbcsr_release(matrix_ps)
549 336 : CALL dbcsr_release(matrix_sp)
550 336 : CALL timestop(handle)
551 336 : END SUBROUTINE reverse_mcweeny_step
552 :
553 : END MODULE admm_dm_methods
|
