Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : ! **************************************************************************************************
15 : MODULE 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_uplo_to_full
28 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
29 : cp_fm_cholesky_invert
30 : USE cp_fm_types, ONLY: cp_fm_get_info,&
31 : cp_fm_to_fm,&
32 : cp_fm_type
33 : USE hfx_admm_utils, ONLY: create_admm_xc_section
34 : USE input_constants, ONLY: do_admm_basis_projection,&
35 : do_admm_purify_none
36 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
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
46 : USE qs_collocate_density, ONLY: calculate_rho_elec
47 : USE qs_environment_types, ONLY: get_qs_env,&
48 : qs_environment_type,&
49 : set_qs_env
50 : USE qs_gapw_densities, ONLY: prepare_gapw_den
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
56 : USE qs_rho_atom_methods, ONLY: calculate_rho_atom_coeff
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 ***
70 : PUBLIC :: rtp_admm_calc_rho_aux, rtp_admm_merge_ks_matrix
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 76 : 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 76 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r_aux
92 : TYPE(admm_type), POINTER :: admm_env
93 76 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
94 76 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p_aux, matrix_p_aux_im, &
95 76 : matrix_s_aux_fit, &
96 76 : matrix_s_aux_fit_vs_orb
97 : TYPE(dft_control_type), POINTER :: dft_control
98 76 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
99 : TYPE(mp_para_env_type), POINTER :: para_env
100 76 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
101 76 : 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 76 : CALL timeset(routineN, handle)
108 76 : NULLIFY (admm_env, matrix_p_aux, matrix_p_aux_im, mos, &
109 76 : mos_aux_fit, para_env, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, rho, &
110 76 : 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 76 : 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 76 : rho_aux_fit=rho_aux_fit)
124 76 : gapw = admm_env%do_gapw
125 :
126 76 : nspins = dft_control%nspins
127 :
128 76 : 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 76 : s_mstruct_changed)
133 :
134 152 : 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 76 : 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 76 : ispin)
146 :
147 : !IF GAPW, only do the soft basis with PW
148 76 : basis_type = "AUX_FIT"
149 76 : IF (gapw) THEN
150 16 : basis_type = "AUX_FIT_SOFT"
151 16 : 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 76 : task_list_external=task_list_aux_fit)
161 :
162 : !IF GAPW, also need to atomic densities
163 152 : 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 16 : para_env=para_env)
169 :
170 : CALL prepare_gapw_den(qs_env, local_rho_set=admm_env%admm_gapw_env%local_rho_set, &
171 16 : do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
172 : END IF
173 : END DO
174 76 : CALL set_qs_env(qs_env, admm_env=admm_env)
175 76 : CALL qs_rho_set(rho_aux_fit, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
176 :
177 76 : CALL timestop(handle)
178 :
179 76 : 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 76 : 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 76 : CALL timeset(routineN, handle)
201 :
202 152 : 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 76 : rtp_coeff_aux_fit, ispin)
206 : CASE DEFAULT
207 76 : CPWARN("only purification NONE possible with RTP/EMD at the moment")
208 : END SELECT
209 :
210 76 : CALL timestop(handle)
211 :
212 76 : 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 152 : SUBROUTINE rtp_admm_fit_mo_coeffs(qs_env, admm_env, admm_control, para_env, matrix_s_aux_fit, matrix_s_mixed, &
229 76 : 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 76 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
246 : TYPE(section_vals_type), POINTER :: input, xc_section
247 :
248 76 : CALL timeset(routineN, handle)
249 :
250 76 : NULLIFY (xc_section, qs_kind_set)
251 :
252 76 : IF (.NOT. (ASSOCIATED(admm_env))) THEN
253 : ! setup admm environment
254 0 : CALL get_qs_env(qs_env, input=input, natom=natoms, qs_kind_set=qs_kind_set)
255 0 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
256 0 : CALL admm_env_create(admm_env, admm_control, mos, para_env, natoms, nao_aux_fit)
257 0 : 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 0 : admm_env=admm_env)
260 :
261 0 : IF (admm_control%method /= do_admm_basis_projection) THEN
262 0 : CPWARN("RTP requires BASIS_PROJECTION.")
263 : END IF
264 : END IF
265 :
266 76 : recalc_S = geometry_did_change .OR. (rtp%iter == 0 .AND. (rtp%istep == rtp%i_start))
267 :
268 152 : 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 76 : mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, recalc_S)
272 : CASE DEFAULT
273 76 : CPWARN("Purification method not implemented in combination with RTP")
274 : END SELECT
275 :
276 76 : CALL timestop(handle)
277 :
278 76 : 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 152 : SUBROUTINE rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, &
299 76 : 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 76 : REAL(KIND=dp), DIMENSION(:), POINTER :: occ_num, occ_num_aux
315 76 : 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 76 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
319 : TYPE(section_vals_type), POINTER :: input, xc_section
320 :
321 76 : CALL timeset(routineN, handle)
322 :
323 76 : NULLIFY (dft_control, qs_kind_set)
324 :
325 76 : IF (.NOT. (ASSOCIATED(admm_env))) THEN
326 0 : CALL get_qs_env(qs_env, input=input, natom=natoms, dft_control=dft_control, qs_kind_set=qs_kind_set)
327 0 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
328 0 : CALL admm_env_create(admm_env, dft_control%admm_control, mos, para_env, natoms, nao_aux_fit)
329 0 : 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 0 : admm_env=admm_env)
332 : END IF
333 :
334 76 : nao_aux_fit = admm_env%nao_aux_fit
335 76 : nao_orb = admm_env%nao_orb
336 76 : nspins = SIZE(mos)
337 :
338 : ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed
339 :
340 76 : IF (geometry_did_change) THEN
341 20 : CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix, admm_env%S_inv)
342 20 : CALL cp_fm_uplo_to_full(admm_env%S_inv, admm_env%work_aux_aux)
343 20 : CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S)
344 :
345 20 : CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix, admm_env%Q)
346 :
347 : !! Calculate S'_inverse
348 20 : CALL cp_fm_cholesky_decompose(admm_env%S_inv)
349 20 : CALL cp_fm_cholesky_invert(admm_env%S_inv)
350 : !! Symmetrize the guy
351 20 : CALL cp_fm_uplo_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 20 : admm_env%A)
356 : END IF
357 :
358 : ! *** Calculate the mo_coeffs for the fitting basis
359 152 : DO ispin = 1, nspins
360 76 : nmo = admm_env%nmo(ispin)
361 76 : IF (nmo == 0) CYCLE
362 : !! Lambda = C^(T)*B*C
363 76 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
364 76 : 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 76 : 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 76 : 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 76 : rtp_coeff_aux_fit(2*ispin))
374 :
375 228 : CALL cp_fm_to_fm(rtp_coeff_aux_fit(2*ispin - 1), mo_coeff_aux_fit)
376 : END DO
377 :
378 76 : CALL timestop(handle)
379 :
380 76 : 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 228 : SUBROUTINE calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, &
390 76 : 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 76 : CALL timeset(routineN, handle)
403 :
404 76 : re = 2*ispin - 1; im = 2*ispin
405 76 : alpha = 3*one - REAL(SIZE(rtp_coeff_aux_fit)/2, dp)
406 76 : CALL dbcsr_set(density_matrix_aux, zero)
407 76 : 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 76 : 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 76 : alpha=alpha)
418 :
419 : ! compute the imaginary part of the dm
420 76 : 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 76 : symmetry_mode=-1)
427 :
428 76 : CALL timestop(handle)
429 :
430 76 : END SUBROUTINE calculate_rtp_admm_density
431 :
432 : ! **************************************************************************************************
433 : !> \brief ...
434 : !> \param qs_env ...
435 : ! **************************************************************************************************
436 76 : 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 76 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
444 76 : matrix_ks_aux_fit_im, matrix_ks_im
445 : TYPE(dft_control_type), POINTER :: dft_control
446 :
447 76 : NULLIFY (admm_env, dft_control, matrix_ks, matrix_ks_im, matrix_ks_aux_fit, matrix_ks_aux_fit_im)
448 76 : 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 76 : matrix_ks_im=matrix_ks_im)
455 76 : 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 152 : DO ispin = 1, dft_control%nspins
460 :
461 76 : SELECT CASE (admm_env%purification_method)
462 : CASE (do_admm_purify_none)
463 : CALL rt_merge_ks_matrix_none(ispin, admm_env, &
464 76 : matrix_ks, matrix_ks_aux_fit)
465 : CALL rt_merge_ks_matrix_none(ispin, admm_env, &
466 76 : matrix_ks_im, matrix_ks_aux_fit_im)
467 : CASE DEFAULT
468 76 : CPWARN("only purification NONE possible with RTP/EMD at the moment")
469 : END SELECT
470 :
471 : END DO !spin loop
472 76 : CALL timestop(handle)
473 :
474 76 : 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 152 : 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 152 : CALL timeset(routineN, handle)
498 :
499 152 : counter = counter + 1
500 152 : nao_aux_fit = admm_env%nao_aux_fit
501 152 : nao_orb = admm_env%nao_orb
502 152 : nmo = admm_env%nmo(ispin)
503 : CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks_aux_fit(ispin)%matrix, &
504 152 : matrix_type=dbcsr_type_no_symmetry)
505 152 : CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
506 152 : CALL dbcsr_desymmetrize(matrix_ks_aux_fit(ispin)%matrix, matrix_ks_nosym)
507 :
508 152 : 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 152 : 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 152 : admm_env%work_orb_orb)
518 :
519 152 : CALL dbcsr_get_info(matrix_ks_aux_fit(ispin)%matrix, matrix_type=matrix_type_fit)
520 :
521 : NULLIFY (matrix_k_tilde)
522 152 : ALLOCATE (matrix_k_tilde)
523 : CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
524 152 : name='MATRIX K_tilde', matrix_type=matrix_type_fit)
525 :
526 152 : CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
527 152 : CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
528 152 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE.)
529 :
530 152 : CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
531 :
532 152 : CALL dbcsr_deallocate_matrix(matrix_k_tilde)
533 152 : CALL dbcsr_release(matrix_ks_nosym)
534 :
535 152 : CALL timestop(handle)
536 :
537 152 : END SUBROUTINE rt_merge_ks_matrix_none
538 :
539 : END MODULE rtp_admm_methods
|
