(git:0de0cc2)
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-2024 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 ! **************************************************************************************************
16  USE admm_types, ONLY: admm_env_create,&
17  admm_type,&
19  USE cp_control_types, ONLY: admm_control_type,&
20  dft_control_type
27  USE cp_fm_types, ONLY: cp_fm_get_info,&
28  cp_fm_to_fm,&
29  cp_fm_type
30  USE dbcsr_api, ONLY: &
31  dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
32  dbcsr_get_info, dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
37  section_vals_type
38  USE kinds, ONLY: default_string_length,&
39  dp
40  USE mathconstants, ONLY: one,&
41  zero
42  USE message_passing, ONLY: mp_para_env_type
43  USE parallel_gemm_api, ONLY: parallel_gemm
44  USE pw_types, ONLY: pw_c1d_gs_type,&
45  pw_r3d_rs_type
47  USE qs_environment_types, ONLY: get_qs_env,&
48  qs_environment_type,&
51  USE qs_kind_types, ONLY: get_qs_kind_set,&
52  qs_kind_type
53  USE qs_ks_types, ONLY: qs_ks_env_type
54  USE qs_mo_types, ONLY: get_mo_set,&
55  mo_set_type
57  USE qs_rho_types, ONLY: qs_rho_get,&
58  qs_rho_set,&
59  qs_rho_type
60  USE rt_propagation_types, ONLY: get_rtp,&
61  rt_prop_type
62  USE task_list_types, ONLY: task_list_type
63 #include "./base/base_uses.f90"
64 
65  IMPLICIT NONE
66 
67  PRIVATE
68 
69  ! *** Public subroutines ***
71 
72  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rtp_admm_methods'
73 
74 CONTAINS
75 
76 ! **************************************************************************************************
77 !> \brief Compute the ADMM density matrix in case of rtp (complex MO's)
78 !>
79 !> \param qs_env ...
80 !> \par History
81 ! **************************************************************************************************
82  SUBROUTINE rtp_admm_calc_rho_aux(qs_env)
83 
84  TYPE(qs_environment_type), POINTER :: qs_env
85 
86  CHARACTER(LEN=*), PARAMETER :: routinen = 'rtp_admm_calc_rho_aux'
87 
88  CHARACTER(LEN=default_string_length) :: basis_type
89  INTEGER :: handle, ispin, nspins
90  LOGICAL :: gapw, s_mstruct_changed
91  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r_aux
92  TYPE(admm_type), POINTER :: admm_env
93  TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
94  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p_aux, matrix_p_aux_im, &
95  matrix_s_aux_fit, &
96  matrix_s_aux_fit_vs_orb
97  TYPE(dft_control_type), POINTER :: dft_control
98  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
99  TYPE(mp_para_env_type), POINTER :: para_env
100  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
101  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_aux
102  TYPE(qs_ks_env_type), POINTER :: ks_env
103  TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
104  TYPE(rt_prop_type), POINTER :: rtp
105  TYPE(task_list_type), POINTER :: task_list_aux_fit
106 
107  CALL timeset(routinen, handle)
108  NULLIFY (admm_env, matrix_p_aux, matrix_p_aux_im, mos, &
109  mos_aux_fit, para_env, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, rho, &
110  ks_env, dft_control, tot_rho_r_aux, rho_r_aux, rho_g_aux, task_list_aux_fit)
111 
112  CALL get_qs_env(qs_env, &
113  admm_env=admm_env, &
114  ks_env=ks_env, &
115  dft_control=dft_control, &
116  para_env=para_env, &
117  mos=mos, &
118  rtp=rtp, &
119  rho=rho, &
120  s_mstruct_changed=s_mstruct_changed)
121  CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit, task_list_aux_fit=task_list_aux_fit, &
122  matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, mos_aux_fit=mos_aux_fit, &
123  rho_aux_fit=rho_aux_fit)
124  gapw = admm_env%do_gapw
125 
126  nspins = dft_control%nspins
127 
128  CALL get_rtp(rtp=rtp, admm_mos=rtp_coeff_aux_fit)
129  CALL rtp_admm_fit_mo_coeffs(qs_env, admm_env, dft_control%admm_control, para_env, &
130  matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, &
131  mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, &
132  s_mstruct_changed)
133 
134  DO ispin = 1, nspins
135  CALL qs_rho_get(rho_aux_fit, &
136  rho_ao=matrix_p_aux, &
137  rho_ao_im=matrix_p_aux_im, &
138  rho_r=rho_r_aux, &
139  rho_g=rho_g_aux, &
140  tot_rho_r=tot_rho_r_aux)
141 
142  CALL rtp_admm_calculate_dm(admm_env, rtp_coeff_aux_fit, &
143  matrix_p_aux(ispin)%matrix, &
144  matrix_p_aux_im(ispin)%matrix, &
145  ispin)
146 
147  !IF GAPW, only do the soft basis with PW
148  basis_type = "AUX_FIT"
149  IF (gapw) THEN
150  basis_type = "AUX_FIT_SOFT"
151  task_list_aux_fit => admm_env%admm_gapw_env%task_list
152  END IF
153 
154  CALL calculate_rho_elec(matrix_p=matrix_p_aux(ispin)%matrix, &
155  rho=rho_r_aux(ispin), &
156  rho_gspace=rho_g_aux(ispin), &
157  total_rho=tot_rho_r_aux(ispin), &
158  ks_env=ks_env, soft_valid=.false., &
159  basis_type="AUX_FIT", &
160  task_list_external=task_list_aux_fit)
161 
162  !IF GAPW, also need to atomic densities
163  IF (gapw) THEN
164  CALL calculate_rho_atom_coeff(qs_env, matrix_p_aux, &
165  rho_atom_set=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
166  qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
167  oce=admm_env%admm_gapw_env%oce, sab=admm_env%sab_aux_fit, &
168  para_env=para_env)
169 
170  CALL prepare_gapw_den(qs_env, local_rho_set=admm_env%admm_gapw_env%local_rho_set, &
171  do_rho0=.false., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
172  END IF
173  END DO
174  CALL set_qs_env(qs_env, admm_env=admm_env)
175  CALL qs_rho_set(rho_aux_fit, rho_r_valid=.true., rho_g_valid=.true.)
176 
177  CALL timestop(handle)
178 
179  END SUBROUTINE rtp_admm_calc_rho_aux
180 
181 ! **************************************************************************************************
182 !> \brief ...
183 !> \param admm_env ...
184 !> \param rtp_coeff_aux_fit ...
185 !> \param density_matrix_aux ...
186 !> \param density_matrix_aux_im ...
187 !> \param ispin ...
188 ! **************************************************************************************************
189  SUBROUTINE rtp_admm_calculate_dm(admm_env, rtp_coeff_aux_fit, density_matrix_aux, &
190  density_matrix_aux_im, ispin)
191  TYPE(admm_type), POINTER :: admm_env
192  TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
193  TYPE(dbcsr_type), POINTER :: density_matrix_aux, density_matrix_aux_im
194  INTEGER, INTENT(in) :: ispin
195 
196  CHARACTER(len=*), PARAMETER :: routinen = 'rtp_admm_calculate_dm'
197 
198  INTEGER :: handle
199 
200  CALL timeset(routinen, handle)
201 
202  SELECT CASE (admm_env%purification_method)
203  CASE (do_admm_purify_none)
204  CALL calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, &
205  rtp_coeff_aux_fit, ispin)
206  CASE DEFAULT
207  cpwarn("only purification NONE possible with RTP/EMD at the moment")
208  END SELECT
209 
210  CALL timestop(handle)
211 
212  END SUBROUTINE rtp_admm_calculate_dm
213 
214 ! **************************************************************************************************
215 !> \brief ...
216 !> \param qs_env ...
217 !> \param admm_env ...
218 !> \param admm_control ...
219 !> \param para_env ...
220 !> \param matrix_s_aux_fit ...
221 !> \param matrix_s_mixed ...
222 !> \param mos ...
223 !> \param mos_aux_fit ...
224 !> \param rtp ...
225 !> \param rtp_coeff_aux_fit ...
226 !> \param geometry_did_change ...
227 ! **************************************************************************************************
228  SUBROUTINE rtp_admm_fit_mo_coeffs(qs_env, admm_env, admm_control, para_env, matrix_s_aux_fit, matrix_s_mixed, &
229  mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, geometry_did_change)
230 
231  TYPE(qs_environment_type), POINTER :: qs_env
232  TYPE(admm_type), POINTER :: admm_env
233  TYPE(admm_control_type), POINTER :: admm_control
234  TYPE(mp_para_env_type), POINTER :: para_env
235  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_mixed
236  TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mos_aux_fit
237  TYPE(rt_prop_type), POINTER :: rtp
238  TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
239  LOGICAL, INTENT(IN) :: geometry_did_change
240 
241  CHARACTER(LEN=*), PARAMETER :: routinen = 'rtp_admm_fit_mo_coeffs'
242 
243  INTEGER :: handle, nao_aux_fit, natoms
244  LOGICAL :: recalc_s
245  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
246  TYPE(section_vals_type), POINTER :: input, xc_section
247 
248  CALL timeset(routinen, handle)
249 
250  NULLIFY (xc_section, qs_kind_set)
251 
252  IF (.NOT. (ASSOCIATED(admm_env))) THEN
253  ! setup admm environment
254  CALL get_qs_env(qs_env, input=input, natom=natoms, qs_kind_set=qs_kind_set)
255  CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
256  CALL admm_env_create(admm_env, admm_control, mos, para_env, natoms, nao_aux_fit)
257  xc_section => section_vals_get_subs_vals(input, "DFT%XC")
258  CALL create_admm_xc_section(x_data=qs_env%x_data, xc_section=xc_section, &
259  admm_env=admm_env)
260 
261  IF (admm_control%method /= do_admm_basis_projection) THEN
262  cpwarn("RTP requires BASIS_PROJECTION.")
263  END IF
264  END IF
265 
266  recalc_s = geometry_did_change .OR. (rtp%iter == 0 .AND. (rtp%istep == rtp%i_start))
267 
268  SELECT CASE (admm_env%purification_method)
269  CASE (do_admm_purify_none)
270  CALL rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, &
271  mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, recalc_s)
272  CASE DEFAULT
273  cpwarn("Purification method not implemented in combination with RTP")
274  END SELECT
275 
276  CALL timestop(handle)
277 
278  END SUBROUTINE rtp_admm_fit_mo_coeffs
279 ! **************************************************************************************************
280 !> \brief Calculates the MO coefficients for the auxiliary fitting basis set
281 !> by minimizing int (psi_i - psi_aux_i)^2 using Lagrangian Multipliers
282 !>
283 !> \param qs_env ...
284 !> \param admm_env The ADMM env
285 !> \param para_env The parallel env
286 !> \param matrix_s_aux_fit the overlap matrix of the auxiliary fitting basis set
287 !> \param matrix_s_mixed the mixed overlap matrix of the auxiliary fitting basis
288 !> set and the orbital basis set
289 !> \param mos the MO's of the orbital basis set
290 !> \param mos_aux_fit the MO's of the auxiliary fitting basis set
291 !> \param rtp ...
292 !> \param rtp_coeff_aux_fit ...
293 !> \param geometry_did_change flag to indicate if the geomtry changed
294 !> \par History
295 !> 05.2008 created [Manuel Guidon]
296 !> \author Manuel Guidon
297 ! **************************************************************************************************
298  SUBROUTINE rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, &
299  mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, geometry_did_change)
300 
301  TYPE(qs_environment_type), POINTER :: qs_env
302  TYPE(admm_type), POINTER :: admm_env
303  TYPE(mp_para_env_type), POINTER :: para_env
304  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_mixed
305  TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mos_aux_fit
306  TYPE(rt_prop_type), POINTER :: rtp
307  TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
308  LOGICAL, INTENT(IN) :: geometry_did_change
309 
310  CHARACTER(LEN=*), PARAMETER :: routinen = 'rtp_fit_mo_coeffs_none'
311 
312  INTEGER :: handle, ispin, nao_aux_fit, nao_orb, &
313  natoms, nmo, nmo_mos, nspins
314  REAL(kind=dp), DIMENSION(:), POINTER :: occ_num, occ_num_aux
315  TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
316  TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
317  TYPE(dft_control_type), POINTER :: dft_control
318  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
319  TYPE(section_vals_type), POINTER :: input, xc_section
320 
321  CALL timeset(routinen, handle)
322 
323  NULLIFY (dft_control, qs_kind_set)
324 
325  IF (.NOT. (ASSOCIATED(admm_env))) THEN
326  CALL get_qs_env(qs_env, input=input, natom=natoms, dft_control=dft_control, qs_kind_set=qs_kind_set)
327  CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
328  CALL admm_env_create(admm_env, dft_control%admm_control, mos, para_env, natoms, nao_aux_fit)
329  xc_section => section_vals_get_subs_vals(input, "DFT%XC")
330  CALL create_admm_xc_section(x_data=qs_env%x_data, xc_section=xc_section, &
331  admm_env=admm_env)
332  END IF
333 
334  nao_aux_fit = admm_env%nao_aux_fit
335  nao_orb = admm_env%nao_orb
336  nspins = SIZE(mos)
337 
338  ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed
339 
340  IF (geometry_did_change) THEN
341  CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix, admm_env%S_inv)
342  CALL cp_fm_upper_to_full(admm_env%S_inv, admm_env%work_aux_aux)
343  CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S)
344 
345  CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix, admm_env%Q)
346 
347  !! Calculate S'_inverse
348  CALL cp_fm_cholesky_decompose(admm_env%S_inv)
349  CALL cp_fm_cholesky_invert(admm_env%S_inv)
350  !! Symmetrize the guy
351  CALL cp_fm_upper_to_full(admm_env%S_inv, admm_env%work_aux_aux)
352  !! Calculate A=S'^(-1)*P
353  CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
354  1.0_dp, admm_env%S_inv, admm_env%Q, 0.0_dp, &
355  admm_env%A)
356  END IF
357 
358  ! *** Calculate the mo_coeffs for the fitting basis
359  DO ispin = 1, nspins
360  nmo = admm_env%nmo(ispin)
361  IF (nmo == 0) cycle
362  !! Lambda = C^(T)*B*C
363  CALL get_rtp(rtp=rtp, mos_new=mos_new)
364  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_num, nmo=nmo_mos)
365  CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit, &
366  occupation_numbers=occ_num_aux)
367 
368  CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
369  1.0_dp, admm_env%A, mos_new(2*ispin - 1), 0.0_dp, &
370  rtp_coeff_aux_fit(2*ispin - 1))
371  CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
372  1.0_dp, admm_env%A, mos_new(2*ispin), 0.0_dp, &
373  rtp_coeff_aux_fit(2*ispin))
374 
375  CALL cp_fm_to_fm(rtp_coeff_aux_fit(2*ispin - 1), mo_coeff_aux_fit)
376  END DO
377 
378  CALL timestop(handle)
379 
380  END SUBROUTINE rtp_fit_mo_coeffs_none
381 
382 ! **************************************************************************************************
383 !> \brief ...
384 !> \param density_matrix_aux ...
385 !> \param density_matrix_aux_im ...
386 !> \param rtp_coeff_aux_fit ...
387 !> \param ispin ...
388 ! **************************************************************************************************
389  SUBROUTINE calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, &
390  rtp_coeff_aux_fit, ispin)
391 
392  TYPE(dbcsr_type), POINTER :: density_matrix_aux, density_matrix_aux_im
393  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: rtp_coeff_aux_fit
394  INTEGER, INTENT(in) :: ispin
395 
396  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rtp_admm_density'
397  REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
398 
399  INTEGER :: handle, im, ncol, re
400  REAL(kind=dp) :: alpha
401 
402  CALL timeset(routinen, handle)
403 
404  re = 2*ispin - 1; im = 2*ispin
405  alpha = 3*one - real(SIZE(rtp_coeff_aux_fit)/2, dp)
406  CALL dbcsr_set(density_matrix_aux, zero)
407  CALL cp_fm_get_info(rtp_coeff_aux_fit(re), ncol_global=ncol)
408  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
409  matrix_v=rtp_coeff_aux_fit(re), &
410  ncol=ncol, &
411  alpha=alpha)
412 
413  ! It is actually complex conjugate but i*i=-1 therefore it must be added
414  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
415  matrix_v=rtp_coeff_aux_fit(im), &
416  ncol=ncol, &
417  alpha=alpha)
418 
419 ! compute the imaginary part of the dm
420  CALL dbcsr_set(density_matrix_aux_im, zero)
421  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux_im, &
422  matrix_v=rtp_coeff_aux_fit(im), &
423  matrix_g=rtp_coeff_aux_fit(re), &
424  ncol=ncol, &
425  alpha=2.0_dp*alpha, &
426  symmetry_mode=-1)
427 
428  CALL timestop(handle)
429 
430  END SUBROUTINE calculate_rtp_admm_density
431 
432 ! **************************************************************************************************
433 !> \brief ...
434 !> \param qs_env ...
435 ! **************************************************************************************************
436  SUBROUTINE rtp_admm_merge_ks_matrix(qs_env)
437  TYPE(qs_environment_type), POINTER :: qs_env
438 
439  CHARACTER(LEN=*), PARAMETER :: routinen = 'rtp_admm_merge_ks_matrix'
440 
441  INTEGER :: handle, ispin
442  TYPE(admm_type), POINTER :: admm_env
443  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
444  matrix_ks_aux_fit_im, matrix_ks_im
445  TYPE(dft_control_type), POINTER :: dft_control
446 
447  NULLIFY (admm_env, dft_control, matrix_ks, matrix_ks_im, matrix_ks_aux_fit, matrix_ks_aux_fit_im)
448  CALL timeset(routinen, handle)
449 
450  CALL get_qs_env(qs_env, &
451  admm_env=admm_env, &
452  dft_control=dft_control, &
453  matrix_ks=matrix_ks, &
454  matrix_ks_im=matrix_ks_im)
455  CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit, matrix_ks_aux_fit_im=matrix_ks_aux_fit_im)
456 
457  !note: the GAPW contribution to ks_aux_fit taken care of in qs_ks_methods.F/update_admm_ks_atom
458 
459  DO ispin = 1, dft_control%nspins
460 
461  SELECT CASE (admm_env%purification_method)
462  CASE (do_admm_purify_none)
463  CALL rt_merge_ks_matrix_none(ispin, admm_env, &
464  matrix_ks, matrix_ks_aux_fit)
465  CALL rt_merge_ks_matrix_none(ispin, admm_env, &
466  matrix_ks_im, matrix_ks_aux_fit_im)
467  CASE DEFAULT
468  cpwarn("only purification NONE possible with RTP/EMD at the moment")
469  END SELECT
470 
471  END DO !spin loop
472  CALL timestop(handle)
473 
474  END SUBROUTINE rtp_admm_merge_ks_matrix
475 
476 ! **************************************************************************************************
477 !> \brief ...
478 !> \param ispin ...
479 !> \param admm_env ...
480 !> \param matrix_ks ...
481 !> \param matrix_ks_aux_fit ...
482 ! **************************************************************************************************
483  SUBROUTINE rt_merge_ks_matrix_none(ispin, admm_env, &
484  matrix_ks, matrix_ks_aux_fit)
485  INTEGER, INTENT(IN) :: ispin
486  TYPE(admm_type), POINTER :: admm_env
487  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit
488 
489  CHARACTER(LEN=*), PARAMETER :: routinen = 'rt_merge_ks_matrix_none'
490 
491  CHARACTER :: matrix_type_fit
492  INTEGER :: handle, nao_aux_fit, nao_orb, nmo
493  INTEGER, SAVE :: counter = 0
494  TYPE(dbcsr_type) :: matrix_ks_nosym
495  TYPE(dbcsr_type), POINTER :: matrix_k_tilde
496 
497  CALL timeset(routinen, handle)
498 
499  counter = counter + 1
500  nao_aux_fit = admm_env%nao_aux_fit
501  nao_orb = admm_env%nao_orb
502  nmo = admm_env%nmo(ispin)
503  CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks_aux_fit(ispin)%matrix, &
504  matrix_type=dbcsr_type_no_symmetry)
505  CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
506  CALL dbcsr_desymmetrize(matrix_ks_aux_fit(ispin)%matrix, matrix_ks_nosym)
507 
508  CALL copy_dbcsr_to_fm(matrix_ks_nosym, admm_env%K(ispin))
509 
510  !! K*A
511  CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
512  1.0_dp, admm_env%K(ispin), admm_env%A, 0.0_dp, &
513  admm_env%work_aux_orb)
514  !! A^T*K*A
515  CALL parallel_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, &
516  1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
517  admm_env%work_orb_orb)
518 
519  CALL dbcsr_get_info(matrix_ks_aux_fit(ispin)%matrix, matrix_type=matrix_type_fit)
520 
521  NULLIFY (matrix_k_tilde)
522  ALLOCATE (matrix_k_tilde)
523  CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
524  name='MATRIX K_tilde', matrix_type=matrix_type_fit)
525 
526  CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
527  CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
528  CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.true.)
529 
530  CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
531 
532  CALL dbcsr_deallocate_matrix(matrix_k_tilde)
533  CALL dbcsr_release(matrix_ks_nosym)
534 
535  CALL timestop(handle)
536 
537  END SUBROUTINE rt_merge_ks_matrix_none
538 
539 END MODULE rtp_admm_methods
Types and set/get functions for auxiliary density matrix methods.
Definition: admm_types.F:15
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
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
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
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,...
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
basic linear algebra operations for full matrices
subroutine, public cp_fm_upper_to_full(matrix, work)
given an upper triangular matrix computes the corresponding full matrix
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_invert(matrix, n, info_out)
used to replace the cholesky decomposition by the inverse
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,...
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
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:1016
Utilities for hfx and admm methods.
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...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_admm_purify_none
integer, parameter, public do_admm_basis_projection
objects that represent the structure of input sections and the data contained in an input section
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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public one
real(kind=dp), parameter, public zero
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
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
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.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, 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, 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, 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, rhs)
Set the QUICKSTEP environment.
subroutine, public prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external)
...
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
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, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
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)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
subroutine, public calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
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
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
Types and set_get for real time propagation depending on runtype and diagonalization method different...
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)
...
Utilities for rtp in combination with admm methods adapted routines from admm_method (author Manuel G...
subroutine, public rtp_admm_merge_ks_matrix(qs_env)
...
subroutine, public rtp_admm_calc_rho_aux(qs_env)
Compute the ADMM density matrix in case of rtp (complex MO's)
types for task lists