(git:f2099e5)
Loading...
Searching...
No Matches
rtp_admm_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 Utilities for rtp in combination with admm methods
10!> adapted routines from admm_method (author Manuel Guidon)
11!>
12!> \par History Use new "force only" overlap routine [07.2014,JGH]
13!> \author Florian Schiffmann
14! **************************************************************************************************
15MODULE rtp_admm_methods
16 USE admm_types, ONLY: admm_env_create,&
17 admm_type,&
18 get_admm_env
19 USE cp_control_types, ONLY: admm_control_type,&
20 dft_control_type
21 USE cp_dbcsr_api, ONLY: &
22 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
23 dbcsr_get_info, dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
24 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
25 copy_fm_to_dbcsr,&
26 cp_dbcsr_plus_fm_fm_t
27 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
28 cp_fm_uplo_to_full
29 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
30 cp_fm_cholesky_invert
31 USE cp_fm_types, ONLY: cp_fm_create,&
32 cp_fm_get_info,&
33 cp_fm_release,&
34 cp_fm_to_fm,&
35 cp_fm_type
36 USE hfx_admm_utils, ONLY: create_admm_xc_section
37 USE input_constants, ONLY: do_admm_basis_projection,&
38 do_admm_purify_none
39 USE input_section_types, ONLY: section_vals_get_subs_vals,&
40 section_vals_type
41 USE kinds, ONLY: default_string_length,&
42 dp
43 USE mathconstants, ONLY: zero
44 USE message_passing, ONLY: mp_para_env_type
45 USE parallel_gemm_api, ONLY: parallel_gemm
46 USE pw_types, ONLY: pw_c1d_gs_type,&
47 pw_r3d_rs_type
48 USE qs_collocate_density, ONLY: calculate_rho_elec
49 USE qs_environment_types, ONLY: get_qs_env,&
50 qs_environment_type,&
51 set_qs_env
52 USE qs_gapw_densities, ONLY: prepare_gapw_den
53 USE qs_kind_types, ONLY: get_qs_kind_set,&
54 qs_kind_type
55 USE qs_ks_types, ONLY: qs_ks_env_type
56 USE qs_mo_types, ONLY: get_mo_set,&
57 mo_set_type
58 USE qs_rho_atom_methods, ONLY: calculate_rho_atom_coeff
59 USE qs_rho_types, ONLY: qs_rho_get,&
60 qs_rho_set,&
61 qs_rho_type
62 USE rt_propagation_types, ONLY: get_rtp,&
63 rt_prop_type
64 USE task_list_types, ONLY: task_list_type
65#include "./base/base_uses.f90"
66
67 IMPLICIT NONE
68
69 PRIVATE
70
71 ! *** Public subroutines ***
72 PUBLIC :: rtp_admm_calc_rho_aux, rtp_admm_merge_ks_matrix
73
74 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rtp_admm_methods'
75
76CONTAINS
77
78! **************************************************************************************************
79!> \brief Compute the ADMM density matrix in case of rtp (complex MO's)
80!>
81!> \param qs_env ...
82!> \par History
83! **************************************************************************************************
84 SUBROUTINE rtp_admm_calc_rho_aux(qs_env)
85
86 TYPE(qs_environment_type), POINTER :: qs_env
87
88 CHARACTER(LEN=*), PARAMETER :: routinen = 'rtp_admm_calc_rho_aux'
89
90 CHARACTER(LEN=default_string_length) :: basis_type
91 INTEGER :: handle, ispin, nmo_aux, nspins
92 LOGICAL :: gapw, s_mstruct_changed
93 REAL(kind=dp), DIMENSION(:), POINTER :: occ_num_aux, tot_rho_r_aux
94 TYPE(admm_type), POINTER :: admm_env
95 TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
96 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p_aux, matrix_p_aux_im, &
97 matrix_s_aux_fit, &
98 matrix_s_aux_fit_vs_orb
99 TYPE(dft_control_type), POINTER :: dft_control
100 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
101 TYPE(mp_para_env_type), POINTER :: para_env
102 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
103 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_aux
104 TYPE(qs_ks_env_type), POINTER :: ks_env
105 TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
106 TYPE(rt_prop_type), POINTER :: rtp
107 TYPE(task_list_type), POINTER :: task_list_aux_fit
108
109 CALL timeset(routinen, handle)
110 NULLIFY (admm_env, matrix_p_aux, matrix_p_aux_im, mos, &
111 mos_aux_fit, para_env, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, rho, &
112 ks_env, dft_control, tot_rho_r_aux, rho_r_aux, rho_g_aux, task_list_aux_fit)
113
114 CALL get_qs_env(qs_env, &
115 admm_env=admm_env, &
116 ks_env=ks_env, &
117 dft_control=dft_control, &
118 para_env=para_env, &
119 mos=mos, &
120 rtp=rtp, &
121 rho=rho, &
122 s_mstruct_changed=s_mstruct_changed)
123 CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit, task_list_aux_fit=task_list_aux_fit, &
124 matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, mos_aux_fit=mos_aux_fit, &
125 rho_aux_fit=rho_aux_fit)
126 gapw = admm_env%do_gapw
127
128 nspins = dft_control%nspins
129
130 CALL get_rtp(rtp=rtp, admm_mos=rtp_coeff_aux_fit)
131 CALL rtp_admm_fit_mo_coeffs(qs_env, admm_env, dft_control%admm_control, para_env, &
132 matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, &
133 mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, &
134 s_mstruct_changed)
135
136 DO ispin = 1, nspins
137 CALL qs_rho_get(rho_aux_fit, &
138 rho_ao=matrix_p_aux, &
139 rho_ao_im=matrix_p_aux_im, &
140 rho_r=rho_r_aux, &
141 rho_g=rho_g_aux, &
142 tot_rho_r=tot_rho_r_aux)
143
144 CALL get_mo_set(mos_aux_fit(ispin), occupation_numbers=occ_num_aux, nmo=nmo_aux)
145
146 CALL rtp_admm_calculate_dm(admm_env, rtp_coeff_aux_fit, &
147 matrix_p_aux(ispin)%matrix, &
148 matrix_p_aux_im(ispin)%matrix, &
149 occ_num_aux, ispin)
150
151 !IF GAPW, only do the soft basis with PW
152 basis_type = "AUX_FIT"
153 IF (gapw) THEN
154 basis_type = "AUX_FIT_SOFT"
155 task_list_aux_fit => admm_env%admm_gapw_env%task_list
156 END IF
157
158 CALL calculate_rho_elec(matrix_p=matrix_p_aux(ispin)%matrix, &
159 rho=rho_r_aux(ispin), &
160 rho_gspace=rho_g_aux(ispin), &
161 total_rho=tot_rho_r_aux(ispin), &
162 ks_env=ks_env, soft_valid=.false., &
163 basis_type="AUX_FIT", &
164 task_list_external=task_list_aux_fit)
165
166 !IF GAPW, also need to atomic densities
167 IF (gapw) THEN
168 CALL calculate_rho_atom_coeff(qs_env, matrix_p_aux, &
169 rho_atom_set=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
170 qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
171 oce=admm_env%admm_gapw_env%oce, sab=admm_env%sab_aux_fit, &
172 para_env=para_env)
173
174 CALL prepare_gapw_den(qs_env, local_rho_set=admm_env%admm_gapw_env%local_rho_set, &
175 do_rho0=.false., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
176 END IF
177 END DO
178 CALL set_qs_env(qs_env, admm_env=admm_env)
179 CALL qs_rho_set(rho_aux_fit, rho_r_valid=.true., rho_g_valid=.true.)
180
181 CALL timestop(handle)
182
183 END SUBROUTINE rtp_admm_calc_rho_aux
184
185! **************************************************************************************************
186!> \brief ...
187!> \param admm_env ...
188!> \param rtp_coeff_aux_fit ...
189!> \param density_matrix_aux ...
190!> \param density_matrix_aux_im ...
191!> \param occupation ...
192!> \param ispin ...
193! **************************************************************************************************
194 SUBROUTINE rtp_admm_calculate_dm(admm_env, rtp_coeff_aux_fit, density_matrix_aux, &
195 density_matrix_aux_im, occupation, ispin)
196 TYPE(admm_type), POINTER :: admm_env
197 TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
198 TYPE(dbcsr_type), POINTER :: density_matrix_aux, density_matrix_aux_im
199 REAL(kind=dp), DIMENSION(:), INTENT(in) :: occupation
200 INTEGER, INTENT(in) :: ispin
201
202 CHARACTER(len=*), PARAMETER :: routinen = 'rtp_admm_calculate_dm'
203
204 INTEGER :: handle
205
206 CALL timeset(routinen, handle)
207
208 SELECT CASE (admm_env%purification_method)
209 CASE (do_admm_purify_none)
210 CALL calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, &
211 rtp_coeff_aux_fit, occupation, ispin)
212 CASE DEFAULT
213 cpwarn("only purification NONE possible with RTP/EMD at the moment")
214 END SELECT
215
216 CALL timestop(handle)
217
218 END SUBROUTINE rtp_admm_calculate_dm
219
220! **************************************************************************************************
221!> \brief ...
222!> \param qs_env ...
223!> \param admm_env ...
224!> \param admm_control ...
225!> \param para_env ...
226!> \param matrix_s_aux_fit ...
227!> \param matrix_s_mixed ...
228!> \param mos ...
229!> \param mos_aux_fit ...
230!> \param rtp ...
231!> \param rtp_coeff_aux_fit ...
232!> \param geometry_did_change ...
233! **************************************************************************************************
234 SUBROUTINE rtp_admm_fit_mo_coeffs(qs_env, admm_env, admm_control, para_env, matrix_s_aux_fit, matrix_s_mixed, &
235 mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, geometry_did_change)
236
237 TYPE(qs_environment_type), POINTER :: qs_env
238 TYPE(admm_type), POINTER :: admm_env
239 TYPE(admm_control_type), POINTER :: admm_control
240 TYPE(mp_para_env_type), POINTER :: para_env
241 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_mixed
242 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mos_aux_fit
243 TYPE(rt_prop_type), POINTER :: rtp
244 TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
245 LOGICAL, INTENT(IN) :: geometry_did_change
246
247 CHARACTER(LEN=*), PARAMETER :: routinen = 'rtp_admm_fit_mo_coeffs'
248
249 INTEGER :: handle, nao_aux_fit, natoms
250 LOGICAL :: recalc_s
251 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
252 TYPE(section_vals_type), POINTER :: input, xc_section
253
254 CALL timeset(routinen, handle)
255
256 NULLIFY (xc_section, qs_kind_set)
257
258 IF (.NOT. (ASSOCIATED(admm_env))) THEN
259 ! setup admm environment
260 CALL get_qs_env(qs_env, input=input, natom=natoms, qs_kind_set=qs_kind_set)
261 CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
262 CALL admm_env_create(admm_env, admm_control, mos, para_env, natoms, nao_aux_fit)
263 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
264 CALL create_admm_xc_section(x_data=qs_env%x_data, xc_section=xc_section, &
265 admm_env=admm_env)
266
267 IF (admm_control%method /= do_admm_basis_projection) THEN
268 cpwarn("RTP requires BASIS_PROJECTION.")
269 END IF
270 END IF
271
272 recalc_s = geometry_did_change .OR. (rtp%iter == 0 .AND. (rtp%istep == rtp%i_start))
273
274 SELECT CASE (admm_env%purification_method)
275 CASE (do_admm_purify_none)
276 CALL rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, &
277 mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, recalc_s)
278 CASE DEFAULT
279 cpwarn("Purification method not implemented in combination with RTP")
280 END SELECT
281
282 CALL timestop(handle)
283
284 END SUBROUTINE rtp_admm_fit_mo_coeffs
285! **************************************************************************************************
286!> \brief Calculates the MO coefficients for the auxiliary fitting basis set
287!> by minimizing int (psi_i - psi_aux_i)^2 using Lagrangian Multipliers
288!>
289!> \param qs_env ...
290!> \param admm_env The ADMM env
291!> \param para_env The parallel env
292!> \param matrix_s_aux_fit the overlap matrix of the auxiliary fitting basis set
293!> \param matrix_s_mixed the mixed overlap matrix of the auxiliary fitting basis
294!> set and the orbital basis set
295!> \param mos the MO's of the orbital basis set
296!> \param mos_aux_fit the MO's of the auxiliary fitting basis set
297!> \param rtp ...
298!> \param rtp_coeff_aux_fit ...
299!> \param geometry_did_change flag to indicate if the geomtry changed
300!> \par History
301!> 05.2008 created [Manuel Guidon]
302!> \author Manuel Guidon
303! **************************************************************************************************
304 SUBROUTINE rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, &
305 mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, geometry_did_change)
306
307 TYPE(qs_environment_type), POINTER :: qs_env
308 TYPE(admm_type), POINTER :: admm_env
309 TYPE(mp_para_env_type), POINTER :: para_env
310 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_mixed
311 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mos_aux_fit
312 TYPE(rt_prop_type), POINTER :: rtp
313 TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
314 LOGICAL, INTENT(IN) :: geometry_did_change
315
316 CHARACTER(LEN=*), PARAMETER :: routinen = 'rtp_fit_mo_coeffs_none'
317
318 INTEGER :: handle, ispin, nao_aux_fit, nao_orb, &
319 natoms, nmo, nmo_mos, nspins
320 REAL(kind=dp), DIMENSION(:), POINTER :: occ_num, occ_num_aux
321 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
322 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
323 TYPE(dft_control_type), POINTER :: dft_control
324 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
325 TYPE(section_vals_type), POINTER :: input, xc_section
326
327 CALL timeset(routinen, handle)
328
329 NULLIFY (dft_control, qs_kind_set)
330
331 IF (.NOT. (ASSOCIATED(admm_env))) THEN
332 CALL get_qs_env(qs_env, input=input, natom=natoms, dft_control=dft_control, qs_kind_set=qs_kind_set)
333 CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
334 CALL admm_env_create(admm_env, dft_control%admm_control, mos, para_env, natoms, nao_aux_fit)
335 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
336 CALL create_admm_xc_section(x_data=qs_env%x_data, xc_section=xc_section, &
337 admm_env=admm_env)
338 END IF
339
340 nao_aux_fit = admm_env%nao_aux_fit
341 nao_orb = admm_env%nao_orb
342 nspins = SIZE(mos)
343
344 ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed
345
346 IF (geometry_did_change) THEN
347 CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix, admm_env%S_inv)
348 CALL cp_fm_uplo_to_full(admm_env%S_inv, admm_env%work_aux_aux)
349 CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S)
350
351 CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix, admm_env%Q)
352
353 !! Calculate S'_inverse
354 CALL cp_fm_cholesky_decompose(admm_env%S_inv)
355 CALL cp_fm_cholesky_invert(admm_env%S_inv)
356 !! Symmetrize the guy
357 CALL cp_fm_uplo_to_full(admm_env%S_inv, admm_env%work_aux_aux)
358 !! Calculate A=S'^(-1)*P
359 CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
360 1.0_dp, admm_env%S_inv, admm_env%Q, 0.0_dp, &
361 admm_env%A)
362 END IF
363
364 ! *** Calculate the mo_coeffs for the fitting basis
365 DO ispin = 1, nspins
366 nmo = admm_env%nmo(ispin)
367 IF (nmo == 0) cycle
368 !! Lambda = C^(T)*B*C
369 CALL get_rtp(rtp=rtp, mos_new=mos_new)
370 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_num, nmo=nmo_mos)
371 CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit, &
372 occupation_numbers=occ_num_aux)
373
374 CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
375 1.0_dp, admm_env%A, mos_new(2*ispin - 1), 0.0_dp, &
376 rtp_coeff_aux_fit(2*ispin - 1))
377 CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
378 1.0_dp, admm_env%A, mos_new(2*ispin), 0.0_dp, &
379 rtp_coeff_aux_fit(2*ispin))
380
381 CALL cp_fm_to_fm(rtp_coeff_aux_fit(2*ispin - 1), mo_coeff_aux_fit)
382 END DO
383
384 CALL timestop(handle)
385
386 END SUBROUTINE rtp_fit_mo_coeffs_none
387
388! **************************************************************************************************
389!> \brief ...
390!> \param density_matrix_aux ...
391!> \param density_matrix_aux_im ...
392!> \param rtp_coeff_aux_fit ...
393!> \param occupation ...
394!> \param ispin ...
395! **************************************************************************************************
396 SUBROUTINE calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, &
397 rtp_coeff_aux_fit, occupation, ispin)
398
399 TYPE(dbcsr_type), POINTER :: density_matrix_aux, density_matrix_aux_im
400 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: rtp_coeff_aux_fit
401 REAL(kind=dp), DIMENSION(:), INTENT(in) :: occupation
402 INTEGER, INTENT(in) :: ispin
403
404 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rtp_admm_density'
405 REAL(kind=dp), PARAMETER :: zero = 0.0_dp
406
407 INTEGER :: handle, im, ncol, re
408 REAL(kind=dp) :: alpha
409 TYPE(cp_fm_type) :: fm_tmp
410
411 CALL timeset(routinen, handle)
412
413 re = 2*ispin - 1; im = 2*ispin
414
415 CALL dbcsr_set(density_matrix_aux, zero)
416 CALL cp_fm_get_info(rtp_coeff_aux_fit(re), ncol_global=ncol)
417 alpha = occupation(1)
418 IF (all(occupation == alpha)) THEN
419 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
420 matrix_v=rtp_coeff_aux_fit(re), &
421 ncol=ncol, &
422 alpha=alpha)
423
424 ! It is actually complex conjugate but i*i=-1 therefore it must be added
425 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
426 matrix_v=rtp_coeff_aux_fit(im), &
427 ncol=ncol, &
428 alpha=alpha)
429
430 ! compute the imaginary part of the dm
431 CALL dbcsr_set(density_matrix_aux_im, zero)
432 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux_im, &
433 matrix_v=rtp_coeff_aux_fit(im), &
434 matrix_g=rtp_coeff_aux_fit(re), &
435 ncol=ncol, &
436 alpha=2.0_dp*alpha, &
437 symmetry_mode=-1)
438 ELSE
439 CALL cp_fm_create(fm_tmp, rtp_coeff_aux_fit(1)%matrix_struct)
440 CALL cp_fm_to_fm(rtp_coeff_aux_fit(re), fm_tmp)
441 CALL cp_fm_column_scale(fm_tmp, occupation(1:ncol))
442 alpha = 1.0_dp
443 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
444 matrix_v=rtp_coeff_aux_fit(re), &
445 matrix_g=fm_tmp, &
446 ncol=ncol, alpha=alpha)
447 ! It is actually complex conjugate but i*i=-1 therefore it must be added
448 CALL cp_fm_to_fm(rtp_coeff_aux_fit(im), fm_tmp)
449 CALL cp_fm_column_scale(fm_tmp, occupation(1:ncol))
450 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
451 matrix_v=rtp_coeff_aux_fit(im), &
452 matrix_g=fm_tmp, &
453 ncol=ncol, alpha=alpha)
454 ! compute the imaginary part of the dm
455 CALL dbcsr_set(density_matrix_aux_im, zero)
456 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux_im, &
457 matrix_v=fm_tmp, &
458 matrix_g=rtp_coeff_aux_fit(re), &
459 ncol=ncol, alpha=2.0_dp*alpha, &
460 symmetry_mode=-1)
461 CALL cp_fm_release(fm_tmp)
462 END IF
463
464 CALL timestop(handle)
465
466 END SUBROUTINE calculate_rtp_admm_density
467
468! **************************************************************************************************
469!> \brief ...
470!> \param qs_env ...
471! **************************************************************************************************
472 SUBROUTINE rtp_admm_merge_ks_matrix(qs_env)
473 TYPE(qs_environment_type), POINTER :: qs_env
474
475 CHARACTER(LEN=*), PARAMETER :: routinen = 'rtp_admm_merge_ks_matrix'
476
477 INTEGER :: handle, ispin
478 TYPE(admm_type), POINTER :: admm_env
479 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
480 matrix_ks_aux_fit_im, matrix_ks_im
481 TYPE(dft_control_type), POINTER :: dft_control
482
483 NULLIFY (admm_env, dft_control, matrix_ks, matrix_ks_im, matrix_ks_aux_fit, matrix_ks_aux_fit_im)
484 CALL timeset(routinen, handle)
485
486 CALL get_qs_env(qs_env, &
487 admm_env=admm_env, &
488 dft_control=dft_control, &
489 matrix_ks=matrix_ks, &
490 matrix_ks_im=matrix_ks_im)
491 CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit, matrix_ks_aux_fit_im=matrix_ks_aux_fit_im)
492
493 !note: the GAPW contribution to ks_aux_fit taken care of in qs_ks_methods.F/update_admm_ks_atom
494
495 DO ispin = 1, dft_control%nspins
496
497 SELECT CASE (admm_env%purification_method)
498 CASE (do_admm_purify_none)
499 CALL rt_merge_ks_matrix_none(ispin, admm_env, &
500 matrix_ks, matrix_ks_aux_fit)
501 CALL rt_merge_ks_matrix_none(ispin, admm_env, &
502 matrix_ks_im, matrix_ks_aux_fit_im)
503 CASE DEFAULT
504 cpwarn("only purification NONE possible with RTP/EMD at the moment")
505 END SELECT
506
507 END DO !spin loop
508 CALL timestop(handle)
509
510 END SUBROUTINE rtp_admm_merge_ks_matrix
511
512! **************************************************************************************************
513!> \brief ...
514!> \param ispin ...
515!> \param admm_env ...
516!> \param matrix_ks ...
517!> \param matrix_ks_aux_fit ...
518! **************************************************************************************************
519 SUBROUTINE rt_merge_ks_matrix_none(ispin, admm_env, &
520 matrix_ks, matrix_ks_aux_fit)
521 INTEGER, INTENT(IN) :: ispin
522 TYPE(admm_type), POINTER :: admm_env
523 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit
524
525 CHARACTER(LEN=*), PARAMETER :: routinen = 'rt_merge_ks_matrix_none'
526
527 CHARACTER :: matrix_type_fit
528 INTEGER :: handle, nao_aux_fit, nao_orb, nmo
529 INTEGER, SAVE :: counter = 0
530 TYPE(dbcsr_type) :: matrix_ks_nosym
531 TYPE(dbcsr_type), POINTER :: matrix_k_tilde
532
533 CALL timeset(routinen, handle)
534
535 counter = counter + 1
536 nao_aux_fit = admm_env%nao_aux_fit
537 nao_orb = admm_env%nao_orb
538 nmo = admm_env%nmo(ispin)
539 CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks_aux_fit(ispin)%matrix, &
540 matrix_type=dbcsr_type_no_symmetry)
541 CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
542 CALL dbcsr_desymmetrize(matrix_ks_aux_fit(ispin)%matrix, matrix_ks_nosym)
543
544 CALL copy_dbcsr_to_fm(matrix_ks_nosym, admm_env%K(ispin))
545
546 !! K*A
547 CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
548 1.0_dp, admm_env%K(ispin), admm_env%A, 0.0_dp, &
549 admm_env%work_aux_orb)
550 !! A^T*K*A
551 CALL parallel_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, &
552 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
553 admm_env%work_orb_orb)
554
555 CALL dbcsr_get_info(matrix_ks_aux_fit(ispin)%matrix, matrix_type=matrix_type_fit)
556
557 NULLIFY (matrix_k_tilde)
558 ALLOCATE (matrix_k_tilde)
559 CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
560 name='MATRIX K_tilde', matrix_type=matrix_type_fit)
561
562 CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
563 CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
564 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.true.)
565
566 CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
567
568 CALL dbcsr_deallocate_matrix(matrix_k_tilde)
569 CALL dbcsr_release(matrix_ks_nosym)
570
571 CALL timestop(handle)
572
573 END SUBROUTINE rt_merge_ks_matrix_none
574
575END MODULE rtp_admm_methods
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
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
admm_types::admm_env_create
subroutine, public admm_env_create(admm_env, admm_control, mos, para_env, natoms, nao_aux_fit, blacs_env_ext)
creates ADMM environment, initializes the basic types
Definition admm_types.F:220
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_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
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_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
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_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:232
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Definition cp_dbcsr_operations.F:944
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:129
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1690
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1584
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_invert
subroutine, public cp_fm_cholesky_invert(matrix, n, info_out)
used to replace the cholesky decomposition by the inverse
Definition cp_fm_cholesky.F:182
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_fm_cholesky.F:61
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
hfx_admm_utils
Utilities for hfx and admm methods.
Definition hfx_admm_utils.F:17
hfx_admm_utils::create_admm_xc_section
subroutine, public create_admm_xc_section(x_data, xc_section, admm_env)
This routine modifies the xc section depending on the potential type used for the HF exchange and the...
Definition hfx_admm_utils.F:1471
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_purify_none
integer, parameter, public do_admm_purify_none
Definition input_constants.F:906
input_constants::do_admm_basis_projection
integer, parameter, public do_admm_basis_projection
Definition input_constants.F:915
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::zero
real(kind=dp), parameter, public zero
Definition mathconstants.F:123
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
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:545
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, mimic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1135
qs_gapw_densities
Definition qs_gapw_densities.F:9
qs_gapw_densities::prepare_gapw_den
subroutine, public prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external, pw_env_sub)
...
Definition qs_gapw_densities.F:57
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:984
qs_ks_types
Definition qs_ks_types.F:15
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:445
qs_rho_atom_methods
Definition qs_rho_atom_methods.F:7
qs_rho_atom_methods::calculate_rho_atom_coeff
subroutine, public calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
...
Definition qs_rho_atom_methods.F:561
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
rt_propagation_types
Types and set_get for real time propagation depending on runtype and diagonalization method different...
Definition rt_propagation_types.F:21
rt_propagation_types::get_rtp
subroutine, public get_rtp(rtp, exp_h_old, exp_h_new, h_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, s_inv, s_half, s_minus_half, b_mat, c_mat, propagator_matrix, mixing, mixing_factor, s_der, dt, nsteps, sinvh, sinvh_imag, sinvb, admm_mos)
...
Definition rt_propagation_types.F:416
rtp_admm_methods
Utilities for rtp in combination with admm methods adapted routines from admm_method (author Manuel G...
Definition rtp_admm_methods.F:15
rtp_admm_methods::rtp_admm_merge_ks_matrix
subroutine, public rtp_admm_merge_ks_matrix(qs_env)
...
Definition rtp_admm_methods.F:473
rtp_admm_methods::rtp_admm_calc_rho_aux
subroutine, public rtp_admm_calc_rho_aux(qs_env)
Compute the ADMM density matrix in case of rtp (complex MO's).
Definition rtp_admm_methods.F:85
task_list_types
types for task lists
Definition task_list_types.F:14
admm_types::admm_type
stores some data used in wavefunction fitting
Definition admm_types.F:120
cp_control_types::admm_control_type
Definition cp_control_types.F:464
cp_control_types::dft_control_type
Definition cp_control_types.F:780
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
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:224
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
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_mo_types::mo_set_type
Definition qs_mo_types.F:53
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
rt_propagation_types::rt_prop_type
Definition rt_propagation_types.F:76
task_list_types::task_list_type
Definition task_list_types.F:58