Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_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 :
76 : CONTAINS
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 114 : 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 114 : REAL(KIND=dp), DIMENSION(:), POINTER :: occ_num_aux, tot_rho_r_aux
94 : TYPE(admm_type), POINTER :: admm_env
95 114 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
96 114 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p_aux, matrix_p_aux_im, &
97 114 : matrix_s_aux_fit, &
98 114 : matrix_s_aux_fit_vs_orb
99 : TYPE(dft_control_type), POINTER :: dft_control
100 114 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
101 : TYPE(mp_para_env_type), POINTER :: para_env
102 114 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
103 114 : 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 114 : CALL timeset(routineN, handle)
110 114 : NULLIFY (admm_env, matrix_p_aux, matrix_p_aux_im, mos, &
111 114 : mos_aux_fit, para_env, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, rho, &
112 114 : 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 114 : 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 114 : rho_aux_fit=rho_aux_fit)
126 114 : gapw = admm_env%do_gapw
127 :
128 114 : nspins = dft_control%nspins
129 :
130 114 : 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 114 : s_mstruct_changed)
135 :
136 266 : 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 152 : tot_rho_r=tot_rho_r_aux)
143 :
144 152 : 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 152 : occ_num_aux, ispin)
150 :
151 : !IF GAPW, only do the soft basis with PW
152 : basis_type = "AUX_FIT"
153 152 : IF (gapw) THEN
154 : basis_type = "AUX_FIT_SOFT"
155 60 : 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 152 : task_list_external=task_list_aux_fit)
165 :
166 : !IF GAPW, also need to atomic densities
167 418 : 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 60 : para_env=para_env)
173 :
174 : CALL prepare_gapw_den(qs_env, local_rho_set=admm_env%admm_gapw_env%local_rho_set, &
175 60 : do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
176 : END IF
177 : END DO
178 114 : CALL set_qs_env(qs_env, admm_env=admm_env)
179 114 : CALL qs_rho_set(rho_aux_fit, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
180 :
181 114 : CALL timestop(handle)
182 :
183 114 : 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 304 : SUBROUTINE rtp_admm_calculate_dm(admm_env, rtp_coeff_aux_fit, density_matrix_aux, &
195 152 : 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 152 : CALL timeset(routineN, handle)
207 :
208 304 : 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 152 : rtp_coeff_aux_fit, occupation, ispin)
212 : CASE DEFAULT
213 152 : CPWARN("only purification NONE possible with RTP/EMD at the moment")
214 : END SELECT
215 :
216 152 : CALL timestop(handle)
217 :
218 152 : 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 228 : SUBROUTINE rtp_admm_fit_mo_coeffs(qs_env, admm_env, admm_control, para_env, matrix_s_aux_fit, matrix_s_mixed, &
235 114 : 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 114 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
252 : TYPE(section_vals_type), POINTER :: input, xc_section
253 :
254 114 : CALL timeset(routineN, handle)
255 :
256 114 : NULLIFY (xc_section, qs_kind_set)
257 :
258 114 : IF (.NOT. (ASSOCIATED(admm_env))) THEN
259 : ! setup admm environment
260 0 : CALL get_qs_env(qs_env, input=input, natom=natoms, qs_kind_set=qs_kind_set)
261 0 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
262 0 : CALL admm_env_create(admm_env, admm_control, mos, para_env, natoms, nao_aux_fit)
263 0 : 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 0 : admm_env=admm_env)
266 :
267 0 : IF (admm_control%method /= do_admm_basis_projection) THEN
268 0 : CPWARN("RTP requires BASIS_PROJECTION.")
269 : END IF
270 : END IF
271 :
272 114 : recalc_S = geometry_did_change .OR. (rtp%iter == 0 .AND. (rtp%istep == rtp%i_start))
273 :
274 228 : 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 114 : mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, recalc_S)
278 : CASE DEFAULT
279 114 : CPWARN("Purification method not implemented in combination with RTP")
280 : END SELECT
281 :
282 114 : CALL timestop(handle)
283 :
284 114 : 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 228 : SUBROUTINE rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, &
305 114 : 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 114 : REAL(KIND=dp), DIMENSION(:), POINTER :: occ_num, occ_num_aux
321 114 : 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 114 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
325 : TYPE(section_vals_type), POINTER :: input, xc_section
326 :
327 114 : CALL timeset(routineN, handle)
328 :
329 114 : NULLIFY (dft_control, qs_kind_set)
330 :
331 114 : IF (.NOT. (ASSOCIATED(admm_env))) THEN
332 0 : CALL get_qs_env(qs_env, input=input, natom=natoms, dft_control=dft_control, qs_kind_set=qs_kind_set)
333 0 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
334 0 : CALL admm_env_create(admm_env, dft_control%admm_control, mos, para_env, natoms, nao_aux_fit)
335 0 : 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 0 : admm_env=admm_env)
338 : END IF
339 :
340 114 : nao_aux_fit = admm_env%nao_aux_fit
341 114 : nao_orb = admm_env%nao_orb
342 114 : nspins = SIZE(mos)
343 :
344 : ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed
345 :
346 114 : IF (geometry_did_change) THEN
347 30 : CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix, admm_env%S_inv)
348 30 : CALL cp_fm_uplo_to_full(admm_env%S_inv, admm_env%work_aux_aux)
349 30 : CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S)
350 :
351 30 : CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix, admm_env%Q)
352 :
353 : !! Calculate S'_inverse
354 30 : CALL cp_fm_cholesky_decompose(admm_env%S_inv)
355 30 : CALL cp_fm_cholesky_invert(admm_env%S_inv)
356 : !! Symmetrize the guy
357 30 : 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 30 : admm_env%A)
362 : END IF
363 :
364 : ! *** Calculate the mo_coeffs for the fitting basis
365 266 : DO ispin = 1, nspins
366 152 : nmo = admm_env%nmo(ispin)
367 152 : IF (nmo == 0) CYCLE
368 : !! Lambda = C^(T)*B*C
369 152 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
370 152 : 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 152 : 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 152 : 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 152 : rtp_coeff_aux_fit(2*ispin))
380 :
381 418 : CALL cp_fm_to_fm(rtp_coeff_aux_fit(2*ispin - 1), mo_coeff_aux_fit)
382 : END DO
383 :
384 114 : CALL timestop(handle)
385 :
386 114 : 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 456 : SUBROUTINE calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, &
397 152 : 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 152 : CALL timeset(routineN, handle)
412 :
413 152 : re = 2*ispin - 1; im = 2*ispin
414 :
415 152 : CALL dbcsr_set(density_matrix_aux, zero)
416 152 : CALL cp_fm_get_info(rtp_coeff_aux_fit(re), ncol_global=ncol)
417 152 : alpha = occupation(1)
418 628 : 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 108 : 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 108 : alpha=alpha)
429 :
430 : ! compute the imaginary part of the dm
431 108 : 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 108 : symmetry_mode=-1)
438 : ELSE
439 44 : CALL cp_fm_create(fm_tmp, rtp_coeff_aux_fit(1)%matrix_struct)
440 44 : CALL cp_fm_to_fm(rtp_coeff_aux_fit(re), fm_tmp)
441 44 : CALL cp_fm_column_scale(fm_tmp, occupation(1:ncol))
442 44 : 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 44 : ncol=ncol, alpha=alpha)
447 : ! It is actually complex conjugate but i*i=-1 therefore it must be added
448 44 : CALL cp_fm_to_fm(rtp_coeff_aux_fit(im), fm_tmp)
449 44 : 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 44 : ncol=ncol, alpha=alpha)
454 : ! compute the imaginary part of the dm
455 44 : 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 44 : symmetry_mode=-1)
461 44 : CALL cp_fm_release(fm_tmp)
462 : END IF
463 :
464 152 : CALL timestop(handle)
465 :
466 152 : END SUBROUTINE calculate_rtp_admm_density
467 :
468 : ! **************************************************************************************************
469 : !> \brief ...
470 : !> \param qs_env ...
471 : ! **************************************************************************************************
472 114 : 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 114 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
480 114 : matrix_ks_aux_fit_im, matrix_ks_im
481 : TYPE(dft_control_type), POINTER :: dft_control
482 :
483 114 : NULLIFY (admm_env, dft_control, matrix_ks, matrix_ks_im, matrix_ks_aux_fit, matrix_ks_aux_fit_im)
484 114 : 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 114 : matrix_ks_im=matrix_ks_im)
491 114 : 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 266 : DO ispin = 1, dft_control%nspins
496 :
497 114 : SELECT CASE (admm_env%purification_method)
498 : CASE (do_admm_purify_none)
499 : CALL rt_merge_ks_matrix_none(ispin, admm_env, &
500 152 : matrix_ks, matrix_ks_aux_fit)
501 : CALL rt_merge_ks_matrix_none(ispin, admm_env, &
502 152 : matrix_ks_im, matrix_ks_aux_fit_im)
503 : CASE DEFAULT
504 152 : CPWARN("only purification NONE possible with RTP/EMD at the moment")
505 : END SELECT
506 :
507 : END DO !spin loop
508 114 : CALL timestop(handle)
509 :
510 114 : 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 304 : 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 304 : CALL timeset(routineN, handle)
534 :
535 304 : counter = counter + 1
536 304 : nao_aux_fit = admm_env%nao_aux_fit
537 304 : nao_orb = admm_env%nao_orb
538 304 : nmo = admm_env%nmo(ispin)
539 : CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks_aux_fit(ispin)%matrix, &
540 304 : matrix_type=dbcsr_type_no_symmetry)
541 304 : CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
542 304 : CALL dbcsr_desymmetrize(matrix_ks_aux_fit(ispin)%matrix, matrix_ks_nosym)
543 :
544 304 : 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 304 : 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 304 : admm_env%work_orb_orb)
554 :
555 304 : CALL dbcsr_get_info(matrix_ks_aux_fit(ispin)%matrix, matrix_type=matrix_type_fit)
556 :
557 : NULLIFY (matrix_k_tilde)
558 304 : ALLOCATE (matrix_k_tilde)
559 : CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
560 304 : name='MATRIX K_tilde', matrix_type=matrix_type_fit)
561 :
562 304 : CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
563 304 : CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
564 304 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE.)
565 :
566 304 : CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
567 :
568 304 : CALL dbcsr_deallocate_matrix(matrix_k_tilde)
569 304 : CALL dbcsr_release(matrix_ks_nosym)
570 :
571 304 : CALL timestop(handle)
572 :
573 304 : END SUBROUTINE rt_merge_ks_matrix_none
574 :
575 : END MODULE rtp_admm_methods
|
