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_purify_none
38 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
39 : section_vals_type
40 : USE kinds, ONLY: default_string_length,&
41 : dp
42 : USE mathconstants, ONLY: zero
43 : USE message_passing, ONLY: mp_para_env_type
44 : USE parallel_gemm_api, ONLY: parallel_gemm
45 : USE pw_types, ONLY: pw_c1d_gs_type,&
46 : pw_r3d_rs_type
47 : USE qs_collocate_density, ONLY: calculate_rho_elec
48 : USE qs_environment_types, ONLY: get_qs_env,&
49 : qs_environment_type,&
50 : set_qs_env
51 : USE qs_gapw_densities, ONLY: prepare_gapw_den
52 : USE qs_kind_types, ONLY: get_qs_kind_set,&
53 : qs_kind_type
54 : USE qs_ks_types, ONLY: qs_ks_env_type
55 : USE qs_mo_types, ONLY: get_mo_set,&
56 : mo_set_type
57 : USE qs_rho_atom_methods, ONLY: calculate_rho_atom_coeff
58 : USE qs_rho_types, ONLY: qs_rho_get,&
59 : qs_rho_set,&
60 : qs_rho_type
61 : USE rt_propagation_types, ONLY: get_rtp,&
62 : rt_prop_type
63 : USE task_list_types, ONLY: task_list_type
64 : #include "./base/base_uses.f90"
65 :
66 : IMPLICIT NONE
67 :
68 : PRIVATE
69 :
70 : ! *** Public subroutines ***
71 : PUBLIC :: rtp_admm_calc_rho_aux, rtp_admm_merge_ks_matrix
72 :
73 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rtp_admm_methods'
74 :
75 : CONTAINS
76 :
77 : ! **************************************************************************************************
78 : !> \brief Compute the ADMM density matrix in case of rtp (complex MO's)
79 : !>
80 : !> \param qs_env ...
81 : !> \par History
82 : ! **************************************************************************************************
83 114 : SUBROUTINE rtp_admm_calc_rho_aux(qs_env)
84 :
85 : TYPE(qs_environment_type), POINTER :: qs_env
86 :
87 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_admm_calc_rho_aux'
88 :
89 : CHARACTER(LEN=default_string_length) :: basis_type
90 : INTEGER :: handle, ispin, nmo_aux, nspins
91 : LOGICAL :: gapw, s_mstruct_changed
92 114 : REAL(KIND=dp), DIMENSION(:), POINTER :: occ_num_aux, tot_rho_r_aux
93 : TYPE(admm_type), POINTER :: admm_env
94 114 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
95 114 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p_aux, matrix_p_aux_im, &
96 114 : matrix_s_aux_fit, &
97 114 : matrix_s_aux_fit_vs_orb
98 : TYPE(dft_control_type), POINTER :: dft_control
99 114 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
100 : TYPE(mp_para_env_type), POINTER :: para_env
101 114 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
102 114 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_aux
103 : TYPE(qs_ks_env_type), POINTER :: ks_env
104 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
105 : TYPE(rt_prop_type), POINTER :: rtp
106 : TYPE(task_list_type), POINTER :: task_list_aux_fit
107 :
108 114 : CALL timeset(routineN, handle)
109 114 : NULLIFY (admm_env, matrix_p_aux, matrix_p_aux_im, mos, &
110 114 : mos_aux_fit, para_env, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, rho, &
111 114 : ks_env, dft_control, tot_rho_r_aux, rho_r_aux, rho_g_aux, task_list_aux_fit)
112 :
113 : CALL get_qs_env(qs_env, &
114 : admm_env=admm_env, &
115 : ks_env=ks_env, &
116 : dft_control=dft_control, &
117 : para_env=para_env, &
118 : mos=mos, &
119 : rtp=rtp, &
120 : rho=rho, &
121 114 : s_mstruct_changed=s_mstruct_changed)
122 : CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit, task_list_aux_fit=task_list_aux_fit, &
123 : matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, mos_aux_fit=mos_aux_fit, &
124 114 : rho_aux_fit=rho_aux_fit)
125 114 : gapw = admm_env%do_gapw
126 :
127 114 : nspins = dft_control%nspins
128 :
129 114 : CALL get_rtp(rtp=rtp, admm_mos=rtp_coeff_aux_fit)
130 : CALL rtp_admm_fit_mo_coeffs(qs_env, admm_env, dft_control%admm_control, para_env, &
131 : matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, &
132 : mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, &
133 114 : s_mstruct_changed)
134 :
135 266 : DO ispin = 1, nspins
136 : CALL qs_rho_get(rho_aux_fit, &
137 : rho_ao=matrix_p_aux, &
138 : rho_ao_im=matrix_p_aux_im, &
139 : rho_r=rho_r_aux, &
140 : rho_g=rho_g_aux, &
141 152 : tot_rho_r=tot_rho_r_aux)
142 :
143 152 : CALL get_mo_set(mos_aux_fit(ispin), occupation_numbers=occ_num_aux, nmo=nmo_aux)
144 :
145 : CALL rtp_admm_calculate_dm(admm_env, rtp_coeff_aux_fit, &
146 : matrix_p_aux(ispin)%matrix, &
147 : matrix_p_aux_im(ispin)%matrix, &
148 152 : occ_num_aux, ispin)
149 :
150 : !IF GAPW, only do the soft basis with PW
151 : basis_type = "AUX_FIT"
152 152 : IF (gapw) THEN
153 : basis_type = "AUX_FIT_SOFT"
154 60 : task_list_aux_fit => admm_env%admm_gapw_env%task_list
155 : END IF
156 :
157 : CALL calculate_rho_elec(matrix_p=matrix_p_aux(ispin)%matrix, &
158 : rho=rho_r_aux(ispin), &
159 : rho_gspace=rho_g_aux(ispin), &
160 : total_rho=tot_rho_r_aux(ispin), &
161 : ks_env=ks_env, soft_valid=.FALSE., &
162 : basis_type="AUX_FIT", &
163 152 : task_list_external=task_list_aux_fit)
164 :
165 : !IF GAPW, also need to atomic densities
166 418 : IF (gapw) THEN
167 : CALL calculate_rho_atom_coeff(qs_env, matrix_p_aux, &
168 : rho_atom_set=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
169 : qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
170 : oce=admm_env%admm_gapw_env%oce, sab=admm_env%sab_aux_fit, &
171 60 : para_env=para_env)
172 :
173 : CALL prepare_gapw_den(qs_env, local_rho_set=admm_env%admm_gapw_env%local_rho_set, &
174 60 : do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
175 : END IF
176 : END DO
177 114 : CALL set_qs_env(qs_env, admm_env=admm_env)
178 114 : CALL qs_rho_set(rho_aux_fit, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
179 :
180 114 : CALL timestop(handle)
181 :
182 114 : END SUBROUTINE rtp_admm_calc_rho_aux
183 :
184 : ! **************************************************************************************************
185 : !> \brief ...
186 : !> \param admm_env ...
187 : !> \param rtp_coeff_aux_fit ...
188 : !> \param density_matrix_aux ...
189 : !> \param density_matrix_aux_im ...
190 : !> \param occupation ...
191 : !> \param ispin ...
192 : ! **************************************************************************************************
193 304 : SUBROUTINE rtp_admm_calculate_dm(admm_env, rtp_coeff_aux_fit, density_matrix_aux, &
194 152 : density_matrix_aux_im, occupation, ispin)
195 : TYPE(admm_type), POINTER :: admm_env
196 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
197 : TYPE(dbcsr_type), POINTER :: density_matrix_aux, density_matrix_aux_im
198 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: occupation
199 : INTEGER, INTENT(in) :: ispin
200 :
201 : CHARACTER(len=*), PARAMETER :: routineN = 'rtp_admm_calculate_dm'
202 :
203 : INTEGER :: handle
204 :
205 152 : CALL timeset(routineN, handle)
206 :
207 304 : SELECT CASE (admm_env%purification_method)
208 : CASE (do_admm_purify_none)
209 : CALL calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, &
210 152 : rtp_coeff_aux_fit, occupation, ispin)
211 : END SELECT
212 :
213 152 : CALL timestop(handle)
214 :
215 152 : END SUBROUTINE rtp_admm_calculate_dm
216 :
217 : ! **************************************************************************************************
218 : !> \brief ...
219 : !> \param qs_env ...
220 : !> \param admm_env ...
221 : !> \param admm_control ...
222 : !> \param para_env ...
223 : !> \param matrix_s_aux_fit ...
224 : !> \param matrix_s_mixed ...
225 : !> \param mos ...
226 : !> \param mos_aux_fit ...
227 : !> \param rtp ...
228 : !> \param rtp_coeff_aux_fit ...
229 : !> \param geometry_did_change ...
230 : ! **************************************************************************************************
231 228 : SUBROUTINE rtp_admm_fit_mo_coeffs(qs_env, admm_env, admm_control, para_env, matrix_s_aux_fit, matrix_s_mixed, &
232 114 : mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, geometry_did_change)
233 :
234 : TYPE(qs_environment_type), POINTER :: qs_env
235 : TYPE(admm_type), POINTER :: admm_env
236 : TYPE(admm_control_type), POINTER :: admm_control
237 : TYPE(mp_para_env_type), POINTER :: para_env
238 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_mixed
239 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mos_aux_fit
240 : TYPE(rt_prop_type), POINTER :: rtp
241 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
242 : LOGICAL, INTENT(IN) :: geometry_did_change
243 :
244 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_admm_fit_mo_coeffs'
245 :
246 : INTEGER :: handle, nao_aux_fit, natoms
247 : LOGICAL :: recalc_S
248 114 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
249 : TYPE(section_vals_type), POINTER :: input, xc_section
250 :
251 114 : CALL timeset(routineN, handle)
252 :
253 114 : NULLIFY (xc_section, qs_kind_set)
254 :
255 114 : IF (.NOT. (ASSOCIATED(admm_env))) THEN
256 : ! setup admm environment
257 0 : CALL get_qs_env(qs_env, input=input, natom=natoms, qs_kind_set=qs_kind_set)
258 0 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
259 0 : CALL admm_env_create(admm_env, admm_control, mos, para_env, natoms, nao_aux_fit)
260 0 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
261 : CALL create_admm_xc_section(x_data=qs_env%x_data, xc_section=xc_section, &
262 0 : admm_env=admm_env)
263 : END IF
264 :
265 114 : recalc_S = geometry_did_change .OR. (rtp%iter == 0 .AND. (rtp%istep == rtp%i_start))
266 :
267 228 : SELECT CASE (admm_env%purification_method)
268 : CASE (do_admm_purify_none)
269 : CALL rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, &
270 114 : mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, recalc_S)
271 : END SELECT
272 :
273 114 : CALL timestop(handle)
274 :
275 114 : END SUBROUTINE rtp_admm_fit_mo_coeffs
276 : ! **************************************************************************************************
277 : !> \brief Calculates the MO coefficients for the auxiliary fitting basis set
278 : !> by minimizing int (psi_i - psi_aux_i)^2 using Lagrangian Multipliers
279 : !>
280 : !> \param qs_env ...
281 : !> \param admm_env The ADMM env
282 : !> \param para_env The parallel env
283 : !> \param matrix_s_aux_fit the overlap matrix of the auxiliary fitting basis set
284 : !> \param matrix_s_mixed the mixed overlap matrix of the auxiliary fitting basis
285 : !> set and the orbital basis set
286 : !> \param mos the MO's of the orbital basis set
287 : !> \param mos_aux_fit the MO's of the auxiliary fitting basis set
288 : !> \param rtp ...
289 : !> \param rtp_coeff_aux_fit ...
290 : !> \param geometry_did_change flag to indicate if the geomtry changed
291 : !> \par History
292 : !> 05.2008 created [Manuel Guidon]
293 : !> \author Manuel Guidon
294 : ! **************************************************************************************************
295 228 : SUBROUTINE rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, &
296 114 : mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, geometry_did_change)
297 :
298 : TYPE(qs_environment_type), POINTER :: qs_env
299 : TYPE(admm_type), POINTER :: admm_env
300 : TYPE(mp_para_env_type), POINTER :: para_env
301 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_mixed
302 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mos_aux_fit
303 : TYPE(rt_prop_type), POINTER :: rtp
304 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: rtp_coeff_aux_fit
305 : LOGICAL, INTENT(IN) :: geometry_did_change
306 :
307 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_fit_mo_coeffs_none'
308 :
309 : INTEGER :: handle, ispin, nao_aux_fit, nao_orb, &
310 : natoms, nmo, nmo_mos, nspins
311 114 : REAL(KIND=dp), DIMENSION(:), POINTER :: occ_num, occ_num_aux
312 114 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
313 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
314 : TYPE(dft_control_type), POINTER :: dft_control
315 114 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
316 : TYPE(section_vals_type), POINTER :: input, xc_section
317 :
318 114 : CALL timeset(routineN, handle)
319 :
320 114 : NULLIFY (dft_control, qs_kind_set)
321 :
322 114 : IF (.NOT. (ASSOCIATED(admm_env))) THEN
323 0 : CALL get_qs_env(qs_env, input=input, natom=natoms, dft_control=dft_control, qs_kind_set=qs_kind_set)
324 0 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
325 0 : CALL admm_env_create(admm_env, dft_control%admm_control, mos, para_env, natoms, nao_aux_fit)
326 0 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
327 : CALL create_admm_xc_section(x_data=qs_env%x_data, xc_section=xc_section, &
328 0 : admm_env=admm_env)
329 : END IF
330 :
331 114 : nao_aux_fit = admm_env%nao_aux_fit
332 114 : nao_orb = admm_env%nao_orb
333 114 : nspins = SIZE(mos)
334 :
335 : ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed
336 :
337 114 : IF (geometry_did_change) THEN
338 30 : CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix, admm_env%S_inv)
339 30 : CALL cp_fm_uplo_to_full(admm_env%S_inv, admm_env%work_aux_aux)
340 30 : CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S)
341 :
342 30 : CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix, admm_env%Q)
343 :
344 : !! Calculate S'_inverse
345 30 : CALL cp_fm_cholesky_decompose(admm_env%S_inv)
346 30 : CALL cp_fm_cholesky_invert(admm_env%S_inv)
347 : !! Symmetrize the guy
348 30 : CALL cp_fm_uplo_to_full(admm_env%S_inv, admm_env%work_aux_aux)
349 : !! Calculate A=S'^(-1)*P
350 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
351 : 1.0_dp, admm_env%S_inv, admm_env%Q, 0.0_dp, &
352 30 : admm_env%A)
353 : END IF
354 :
355 : ! *** Calculate the mo_coeffs for the fitting basis
356 266 : DO ispin = 1, nspins
357 152 : nmo = admm_env%nmo(ispin)
358 152 : IF (nmo == 0) CYCLE
359 : !! Lambda = C^(T)*B*C
360 152 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
361 152 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_num, nmo=nmo_mos)
362 : CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit, &
363 152 : occupation_numbers=occ_num_aux)
364 :
365 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
366 : 1.0_dp, admm_env%A, mos_new(2*ispin - 1), 0.0_dp, &
367 152 : rtp_coeff_aux_fit(2*ispin - 1))
368 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
369 : 1.0_dp, admm_env%A, mos_new(2*ispin), 0.0_dp, &
370 152 : rtp_coeff_aux_fit(2*ispin))
371 :
372 418 : CALL cp_fm_to_fm(rtp_coeff_aux_fit(2*ispin - 1), mo_coeff_aux_fit)
373 : END DO
374 :
375 114 : CALL timestop(handle)
376 :
377 114 : END SUBROUTINE rtp_fit_mo_coeffs_none
378 :
379 : ! **************************************************************************************************
380 : !> \brief ...
381 : !> \param density_matrix_aux ...
382 : !> \param density_matrix_aux_im ...
383 : !> \param rtp_coeff_aux_fit ...
384 : !> \param occupation ...
385 : !> \param ispin ...
386 : ! **************************************************************************************************
387 456 : SUBROUTINE calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, &
388 152 : rtp_coeff_aux_fit, occupation, ispin)
389 :
390 : TYPE(dbcsr_type), POINTER :: density_matrix_aux, density_matrix_aux_im
391 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: rtp_coeff_aux_fit
392 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: occupation
393 : INTEGER, INTENT(in) :: ispin
394 :
395 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rtp_admm_density'
396 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
397 :
398 : INTEGER :: handle, im, ncol, re
399 : REAL(KIND=dp) :: alpha
400 : TYPE(cp_fm_type) :: fm_tmp
401 :
402 152 : CALL timeset(routineN, handle)
403 :
404 152 : re = 2*ispin - 1; im = 2*ispin
405 :
406 152 : CALL dbcsr_set(density_matrix_aux, zero)
407 152 : CALL cp_fm_get_info(rtp_coeff_aux_fit(re), ncol_global=ncol)
408 152 : alpha = occupation(1)
409 628 : IF (ALL(occupation == alpha)) THEN
410 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
411 : matrix_v=rtp_coeff_aux_fit(re), &
412 : ncol=ncol, &
413 108 : alpha=alpha)
414 :
415 : ! It is actually complex conjugate but i*i=-1 therefore it must be added
416 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
417 : matrix_v=rtp_coeff_aux_fit(im), &
418 : ncol=ncol, &
419 108 : alpha=alpha)
420 :
421 : ! compute the imaginary part of the dm
422 108 : CALL dbcsr_set(density_matrix_aux_im, zero)
423 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux_im, &
424 : matrix_v=rtp_coeff_aux_fit(im), &
425 : matrix_g=rtp_coeff_aux_fit(re), &
426 : ncol=ncol, &
427 : alpha=2.0_dp*alpha, &
428 108 : symmetry_mode=-1)
429 : ELSE
430 44 : CALL cp_fm_create(fm_tmp, rtp_coeff_aux_fit(1)%matrix_struct)
431 44 : CALL cp_fm_to_fm(rtp_coeff_aux_fit(re), fm_tmp)
432 44 : CALL cp_fm_column_scale(fm_tmp, occupation(1:ncol))
433 44 : alpha = 1.0_dp
434 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
435 : matrix_v=rtp_coeff_aux_fit(re), &
436 : matrix_g=fm_tmp, &
437 44 : ncol=ncol, alpha=alpha)
438 : ! It is actually complex conjugate but i*i=-1 therefore it must be added
439 44 : CALL cp_fm_to_fm(rtp_coeff_aux_fit(im), fm_tmp)
440 44 : CALL cp_fm_column_scale(fm_tmp, occupation(1:ncol))
441 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
442 : matrix_v=rtp_coeff_aux_fit(im), &
443 : matrix_g=fm_tmp, &
444 44 : ncol=ncol, alpha=alpha)
445 : ! compute the imaginary part of the dm
446 44 : CALL dbcsr_set(density_matrix_aux_im, zero)
447 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux_im, &
448 : matrix_v=fm_tmp, &
449 : matrix_g=rtp_coeff_aux_fit(re), &
450 : ncol=ncol, alpha=2.0_dp*alpha, &
451 44 : symmetry_mode=-1)
452 44 : CALL cp_fm_release(fm_tmp)
453 : END IF
454 :
455 152 : CALL timestop(handle)
456 :
457 152 : END SUBROUTINE calculate_rtp_admm_density
458 :
459 : ! **************************************************************************************************
460 : !> \brief ...
461 : !> \param qs_env ...
462 : ! **************************************************************************************************
463 114 : SUBROUTINE rtp_admm_merge_ks_matrix(qs_env)
464 : TYPE(qs_environment_type), POINTER :: qs_env
465 :
466 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_admm_merge_ks_matrix'
467 :
468 : INTEGER :: handle, ispin
469 : TYPE(admm_type), POINTER :: admm_env
470 114 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
471 114 : matrix_ks_aux_fit_im, matrix_ks_im
472 : TYPE(dft_control_type), POINTER :: dft_control
473 :
474 114 : NULLIFY (admm_env, dft_control, matrix_ks, matrix_ks_im, matrix_ks_aux_fit, matrix_ks_aux_fit_im)
475 114 : CALL timeset(routineN, handle)
476 :
477 : CALL get_qs_env(qs_env, &
478 : admm_env=admm_env, &
479 : dft_control=dft_control, &
480 : matrix_ks=matrix_ks, &
481 114 : matrix_ks_im=matrix_ks_im)
482 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)
483 :
484 : !note: the GAPW contribution to ks_aux_fit taken care of in qs_ks_methods.F/update_admm_ks_atom
485 :
486 266 : DO ispin = 1, dft_control%nspins
487 :
488 114 : SELECT CASE (admm_env%purification_method)
489 : CASE (do_admm_purify_none)
490 : CALL rt_merge_ks_matrix_none(ispin, admm_env, &
491 152 : matrix_ks, matrix_ks_aux_fit)
492 : CALL rt_merge_ks_matrix_none(ispin, admm_env, &
493 304 : matrix_ks_im, matrix_ks_aux_fit_im)
494 : END SELECT
495 :
496 : END DO !spin loop
497 114 : CALL timestop(handle)
498 :
499 114 : END SUBROUTINE rtp_admm_merge_ks_matrix
500 :
501 : ! **************************************************************************************************
502 : !> \brief ...
503 : !> \param ispin ...
504 : !> \param admm_env ...
505 : !> \param matrix_ks ...
506 : !> \param matrix_ks_aux_fit ...
507 : ! **************************************************************************************************
508 304 : SUBROUTINE rt_merge_ks_matrix_none(ispin, admm_env, &
509 : matrix_ks, matrix_ks_aux_fit)
510 : INTEGER, INTENT(IN) :: ispin
511 : TYPE(admm_type), POINTER :: admm_env
512 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit
513 :
514 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rt_merge_ks_matrix_none'
515 :
516 : CHARACTER :: matrix_type_fit
517 : INTEGER :: handle, nao_aux_fit, nao_orb, nmo
518 : INTEGER, SAVE :: counter = 0
519 : TYPE(dbcsr_type) :: matrix_ks_nosym
520 : TYPE(dbcsr_type), POINTER :: matrix_k_tilde
521 :
522 304 : CALL timeset(routineN, handle)
523 :
524 304 : counter = counter + 1
525 304 : nao_aux_fit = admm_env%nao_aux_fit
526 304 : nao_orb = admm_env%nao_orb
527 304 : nmo = admm_env%nmo(ispin)
528 : CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks_aux_fit(ispin)%matrix, &
529 304 : matrix_type=dbcsr_type_no_symmetry)
530 304 : CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
531 304 : CALL dbcsr_desymmetrize(matrix_ks_aux_fit(ispin)%matrix, matrix_ks_nosym)
532 :
533 304 : CALL copy_dbcsr_to_fm(matrix_ks_nosym, admm_env%K(ispin))
534 :
535 : !! K*A
536 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
537 : 1.0_dp, admm_env%K(ispin), admm_env%A, 0.0_dp, &
538 304 : admm_env%work_aux_orb)
539 : !! A^T*K*A
540 : CALL parallel_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, &
541 : 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
542 304 : admm_env%work_orb_orb)
543 :
544 304 : CALL dbcsr_get_info(matrix_ks_aux_fit(ispin)%matrix, matrix_type=matrix_type_fit)
545 :
546 : NULLIFY (matrix_k_tilde)
547 304 : ALLOCATE (matrix_k_tilde)
548 : CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
549 304 : name='MATRIX K_tilde', matrix_type=matrix_type_fit)
550 :
551 304 : CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
552 304 : CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
553 304 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE.)
554 :
555 304 : CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
556 :
557 304 : CALL dbcsr_deallocate_matrix(matrix_k_tilde)
558 304 : CALL dbcsr_release(matrix_ks_nosym)
559 :
560 304 : CALL timestop(handle)
561 :
562 304 : END SUBROUTINE rt_merge_ks_matrix_none
563 :
564 : END MODULE rtp_admm_methods
|
