Line data Source code
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_api, ONLY: &
20 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_frobenius_norm, dbcsr_get_block_p, &
21 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
22 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
23 : dbcsr_scale, dbcsr_set, dbcsr_type
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 140 : 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 140 : NULLIFY (admm_dm)
65 140 : CALL timeset(routineN, handle)
66 140 : CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
67 :
68 240 : SELECT CASE (admm_dm%method)
69 : CASE (do_admm_basis_projection)
70 100 : CALL map_dm_projection(qs_env)
71 :
72 : CASE (do_admm_blocked_projection)
73 40 : CALL map_dm_blocked(qs_env)
74 :
75 : CASE DEFAULT
76 140 : CPABORT("admm_dm_calc_rho_aux: unknown method")
77 : END SELECT
78 :
79 140 : IF (admm_dm%purify) &
80 20 : CALL purify_mcweeny(qs_env)
81 :
82 140 : CALL update_rho_aux(qs_env)
83 :
84 140 : CALL timestop(handle)
85 140 : 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 140 : 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 140 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge
100 :
101 140 : CALL timeset(routineN, handle)
102 140 : NULLIFY (admm_dm, matrix_ks_merge)
103 :
104 140 : CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
105 :
106 140 : IF (admm_dm%purify) THEN
107 20 : CALL revert_purify_mcweeny(qs_env, matrix_ks_merge)
108 : ELSE
109 120 : CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit=matrix_ks_merge)
110 : END IF
111 :
112 240 : SELECT CASE (admm_dm%method)
113 : CASE (do_admm_basis_projection)
114 100 : CALL merge_dm_projection(qs_env, matrix_ks_merge)
115 :
116 : CASE (do_admm_blocked_projection)
117 40 : CALL merge_dm_blocked(qs_env, matrix_ks_merge)
118 :
119 : CASE DEFAULT
120 140 : CPABORT("admm_dm_merge_ks_matrix: unknown method")
121 : END SELECT
122 :
123 140 : IF (admm_dm%purify) &
124 20 : CALL dbcsr_deallocate_matrix_set(matrix_ks_merge)
125 :
126 140 : CALL timestop(handle)
127 :
128 140 : 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 100 : 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 100 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux, matrix_s_mixed, rho_ao, &
143 100 : 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 100 : NULLIFY (dft_control, admm_dm, matrix_s_aux, matrix_s_mixed, rho, rho_aux)
149 100 : NULLIFY (rho_ao, rho_ao_aux)
150 :
151 100 : 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 100 : matrix_s_aux_fit_vs_orb=matrix_s_mixed, admm_dm=admm_dm)
154 :
155 100 : CALL qs_rho_get(rho, rho_ao=rho_ao)
156 100 : CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux)
157 :
158 100 : IF (s_mstruct_changed) THEN
159 : ! Calculate A = S_aux^(-1) * S_mixed
160 10 : CALL dbcsr_create(matrix_s_aux_inv, template=matrix_s_aux(1)%matrix, matrix_type="N")
161 10 : threshold = MAX(admm_dm%eps_filter, 1.0e-12_dp)
162 10 : CALL invert_Hotelling(matrix_s_aux_inv, matrix_s_aux(1)%matrix, threshold)
163 :
164 10 : IF (.NOT. ASSOCIATED(admm_dm%matrix_A)) THEN
165 10 : ALLOCATE (admm_dm%matrix_A)
166 10 : 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 10 : 0.0_dp, admm_dm%matrix_A)
170 10 : CALL dbcsr_release(matrix_s_aux_inv)
171 : END IF
172 :
173 : ! Calculate P_aux = A * P * A^T
174 100 : CALL dbcsr_create(matrix_tmp, template=admm_dm%matrix_A)
175 260 : DO ispin = 1, dft_control%nspins
176 : CALL dbcsr_multiply("N", "N", 1.0_dp, admm_dm%matrix_A, rho_ao(ispin)%matrix, &
177 160 : 0.0_dp, matrix_tmp)
178 : CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp, admm_dm%matrix_A, &
179 260 : 0.0_dp, rho_ao_aux(ispin)%matrix)
180 : END DO
181 100 : CALL dbcsr_release(matrix_tmp)
182 :
183 100 : 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 40 : 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 40 : REAL(dp), DIMENSION(:, :), POINTER :: sparse_block, sparse_block_aux
196 : TYPE(admm_dm_type), POINTER :: admm_dm
197 : TYPE(dbcsr_iterator_type) :: iter
198 40 : 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 40 : NULLIFY (dft_control, admm_dm, rho, rho_aux, rho_ao, rho_ao_aux)
203 :
204 40 : CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
205 40 : CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux, admm_dm=admm_dm)
206 :
207 40 : CALL qs_rho_get(rho, rho_ao=rho_ao)
208 40 : CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux)
209 :
210 : ! ** set blocked density matrix to 0
211 100 : DO ispin = 1, dft_control%nspins
212 60 : CALL dbcsr_set(rho_ao_aux(ispin)%matrix, 0.0_dp)
213 : ! ** now loop through the list and copy corresponding blocks
214 60 : CALL dbcsr_iterator_start(iter, rho_ao(ispin)%matrix)
215 290 : DO WHILE (dbcsr_iterator_blocks_left(iter))
216 230 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk)
217 290 : IF (admm_dm%block_map(iatom, jatom) == 1) THEN
218 : CALL dbcsr_get_block_p(rho_ao_aux(ispin)%matrix, &
219 160 : row=iatom, col=jatom, BLOCK=sparse_block_aux, found=found)
220 160 : IF (found) &
221 1760 : sparse_block_aux = sparse_block
222 : END IF
223 : END DO
224 160 : CALL dbcsr_iterator_stop(iter)
225 : END DO
226 :
227 40 : END SUBROUTINE map_dm_blocked
228 :
229 : ! **************************************************************************************************
230 : !> \brief Call calculate_rho_elec() for auxiliary density
231 : !> \param qs_env ...
232 : ! **************************************************************************************************
233 140 : SUBROUTINE update_rho_aux(qs_env)
234 : TYPE(qs_environment_type), POINTER :: qs_env
235 :
236 : INTEGER :: ispin
237 140 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r_aux
238 : TYPE(admm_dm_type), POINTER :: admm_dm
239 140 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_aux
240 : TYPE(dft_control_type), POINTER :: dft_control
241 140 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
242 140 : 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 140 : NULLIFY (dft_control, admm_dm, rho_aux, rho_ao_aux, rho_r_aux, rho_g_aux, tot_rho_r_aux, &
248 140 : task_list_aux_fit, ks_env)
249 :
250 140 : 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 140 : 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 140 : tot_rho_r=tot_rho_r_aux)
259 :
260 360 : 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 360 : task_list_external=task_list_aux_fit)
269 : END DO
270 :
271 140 : CALL qs_rho_set(rho_aux, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
272 :
273 140 : 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 100 : 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 100 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
288 : TYPE(dbcsr_type) :: matrix_tmp
289 : TYPE(dft_control_type), POINTER :: dft_control
290 :
291 100 : NULLIFY (admm_dm, dft_control, matrix_ks)
292 :
293 100 : CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks)
294 100 : CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
295 :
296 : ! Calculate K += A^T * K_aux * A
297 100 : CALL dbcsr_create(matrix_tmp, template=admm_dm%matrix_A, matrix_type="N")
298 :
299 260 : 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 160 : 0.0_dp, matrix_tmp)
302 : CALL dbcsr_multiply("T", "N", 1.0_dp, admm_dm%matrix_A, matrix_tmp, &
303 260 : 1.0_dp, matrix_ks(ispin)%matrix)
304 : END DO
305 :
306 100 : CALL dbcsr_release(matrix_tmp)
307 :
308 100 : 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 40 : 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 40 : REAL(dp), DIMENSION(:, :), POINTER :: sparse_block
322 : TYPE(admm_dm_type), POINTER :: admm_dm
323 : TYPE(dbcsr_iterator_type) :: iter
324 40 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
325 : TYPE(dft_control_type), POINTER :: dft_control
326 :
327 40 : NULLIFY (admm_dm, dft_control, matrix_ks)
328 :
329 40 : CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks)
330 40 : CALL get_admm_env(qs_env%admm_env, admm_dm=admm_dm)
331 :
332 100 : DO ispin = 1, dft_control%nspins
333 60 : CALL dbcsr_iterator_start(iter, matrix_ks_merge(ispin)%matrix)
334 290 : DO WHILE (dbcsr_iterator_blocks_left(iter))
335 230 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk)
336 230 : IF (admm_dm%block_map(iatom, jatom) == 0) &
337 600 : sparse_block = 0.0_dp
338 : END DO
339 60 : CALL dbcsr_iterator_stop(iter)
340 160 : CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_ks_merge(ispin)%matrix, 1.0_dp, 1.0_dp)
341 : END DO
342 :
343 40 : 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 20 : 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 20 : 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 20 : CALL timeset(routineN, handle)
366 20 : NULLIFY (dft_control, admm_dm, matrix_s_aux_fit, rho_aux_fit, new_hist_entry, &
367 20 : matrix_p, matrix_s, rho_ao_aux)
368 :
369 20 : unit_nr = cp_logger_get_default_unit_nr()
370 20 : 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 20 : rho_aux_fit=rho_aux_fit, admm_dm=admm_dm)
373 :
374 20 : CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux)
375 :
376 20 : matrix_p => rho_ao_aux(1)%matrix
377 20 : CALL dbcsr_create(matrix_PS, template=matrix_p, matrix_type="N")
378 20 : CALL dbcsr_create(matrix_PSP, template=matrix_p, matrix_type="S")
379 20 : CALL dbcsr_create(matrix_test, template=matrix_p, matrix_type="S")
380 :
381 20 : nspins = dft_control%nspins
382 60 : DO ispin = 1, nspins
383 40 : matrix_p => rho_ao_aux(ispin)%matrix
384 40 : matrix_s => matrix_s_aux_fit(1)%matrix
385 40 : history => admm_dm%mcweeny_history(ispin)%p
386 40 : IF (ASSOCIATED(history)) CPABORT("purify_dm_mcweeny: history already associated")
387 40 : IF (nspins == 1) CALL dbcsr_scale(matrix_p, 0.5_dp)
388 :
389 204 : DO istep = 1, admm_dm%mcweeny_max_steps
390 : ! allocate new element in linked list
391 204 : ALLOCATE (new_hist_entry)
392 204 : new_hist_entry%next => history
393 204 : history => new_hist_entry
394 204 : history%count = istep
395 204 : NULLIFY (new_hist_entry)
396 204 : CALL dbcsr_create(history%m, template=matrix_p, matrix_type="N")
397 204 : 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 204 : 0.0_dp, matrix_ps)
402 :
403 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ps, matrix_p, &
404 204 : 0.0_dp, matrix_psp)
405 :
406 : !test convergence
407 204 : CALL dbcsr_copy(matrix_test, matrix_psp)
408 204 : CALL dbcsr_add(matrix_test, matrix_p, 1.0_dp, -1.0_dp)
409 204 : frob_norm = dbcsr_frobenius_norm(matrix_test)
410 408 : IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,i5,a,f16.8)') "McWeeny-Step", istep, &
411 408 : ": Deviation of idempotency", frob_norm
412 204 : IF (frob_norm < 1000_dp*admm_dm%eps_filter .AND. istep > 1) EXIT
413 :
414 : ! build next P matrix
415 164 : 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 204 : 3.0_dp, matrix_p)
418 : END DO
419 40 : admm_dm%mcweeny_history(ispin)%p => history
420 60 : IF (nspins == 1) CALL dbcsr_scale(matrix_p, 2.0_dp)
421 : END DO
422 :
423 : ! clean up
424 20 : CALL dbcsr_release(matrix_PS)
425 20 : CALL dbcsr_release(matrix_PSP)
426 20 : CALL dbcsr_release(matrix_test)
427 20 : CALL timestop(handle)
428 20 : 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 20 : 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 20 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
445 20 : matrix_s_aux_fit, &
446 20 : 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 20 : CALL timeset(routineN, handle)
452 20 : unit_nr = cp_logger_get_default_unit_nr()
453 20 : NULLIFY (admm_dm, dft_control, matrix_ks, matrix_ks_aux_fit, &
454 20 : matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, &
455 20 : history_next, history_curr, matrix_k)
456 :
457 20 : 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 20 : matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, matrix_ks_aux_fit=matrix_ks_aux_fit)
460 :
461 20 : nspins = dft_control%nspins
462 100 : ALLOCATE (matrix_ks_merge(nspins))
463 :
464 60 : DO ispin = 1, nspins
465 40 : ALLOCATE (matrix_ks_merge(ispin)%matrix)
466 40 : matrix_k => matrix_ks_merge(ispin)%matrix
467 40 : CALL dbcsr_copy(matrix_k, matrix_ks_aux_fit(ispin)%matrix, name="K")
468 40 : history_curr => admm_dm%mcweeny_history(ispin)%p
469 40 : NULLIFY (admm_dm%mcweeny_history(ispin)%p)
470 :
471 : ! reverse McWeeny iteration
472 264 : DO WHILE (ASSOCIATED(history_curr))
473 204 : 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 204 : matrix_p=history_curr%m)
477 204 : CALL dbcsr_release(history_curr%m)
478 204 : history_next => history_curr%next
479 204 : DEALLOCATE (history_curr)
480 204 : history_curr => history_next
481 204 : NULLIFY (history_next)
482 : END DO
483 :
484 : END DO
485 :
486 : ! clean up
487 20 : CALL timestop(handle)
488 :
489 20 : 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 204 : 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 204 : CALL timeset(routineN, handle)
508 204 : CALL dbcsr_create(matrix_ps, template=matrix_p, matrix_type="N")
509 204 : CALL dbcsr_create(matrix_sp, template=matrix_p, matrix_type="N")
510 204 : CALL dbcsr_create(matrix_tmp, template=matrix_p, matrix_type="N")
511 204 : 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 204 : 0.0_dp, matrix_ps)
515 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, matrix_p, &
516 204 : 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 204 : 0.0_dp, matrix_sum)
521 : CALL dbcsr_multiply("N", "N", 3.0_dp, matrix_sp, matrix_k, &
522 204 : 1.0_dp, matrix_sum)
523 :
524 : !matrix_tmp = KPS
525 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_k, matrix_ps, &
526 204 : 0.0_dp, matrix_tmp)
527 : CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_tmp, matrix_ps, &
528 204 : 1.0_dp, matrix_sum)
529 : CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, &
530 204 : 1.0_dp, matrix_sum)
531 :
532 : !matrix_tmp = SPK
533 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sp, matrix_k, &
534 204 : 0.0_dp, matrix_tmp)
535 : CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, &
536 204 : 1.0_dp, matrix_sum)
537 :
538 : ! overwrite matrix_k
539 204 : CALL dbcsr_copy(matrix_k, matrix_sum, name="K from reverse McWeeny")
540 :
541 : ! clean up
542 204 : CALL dbcsr_release(matrix_sum)
543 204 : CALL dbcsr_release(matrix_tmp)
544 204 : CALL dbcsr_release(matrix_ps)
545 204 : CALL dbcsr_release(matrix_sp)
546 204 : CALL timestop(handle)
547 204 : END SUBROUTINE reverse_mcweeny_step
548 :
549 : END MODULE admm_dm_methods
|
