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 Contains ADMM methods which require molecular orbitals
10 : !> \par History
11 : !> 04.2008 created [Manuel Guidon]
12 : !> 12.2019 Made GAPW compatible [A. Bussy]
13 : !> \author Manuel Guidon
14 : ! **************************************************************************************************
15 : MODULE admm_methods
16 : USE admm_types, ONLY: admm_gapw_r3d_rs_type,&
17 : admm_type,&
18 : get_admm_env
19 : USE atomic_kind_types, ONLY: atomic_kind_type
20 : USE bibliography, ONLY: Merlot2014,&
21 : cite_reference
22 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale,&
23 : cp_cfm_scale_and_add,&
24 : cp_cfm_uplo_to_full
25 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose,&
26 : cp_cfm_cholesky_invert
27 : USE cp_cfm_types, ONLY: cp_cfm_create,&
28 : cp_cfm_release,&
29 : cp_cfm_to_cfm,&
30 : cp_cfm_to_fm,&
31 : cp_cfm_type,&
32 : cp_fm_to_cfm
33 : USE cp_control_types, ONLY: dft_control_type
34 : USE cp_dbcsr_api, ONLY: &
35 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
36 : dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
37 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
38 : dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
39 : dbcsr_type_no_symmetry, dbcsr_type_symmetric
40 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
41 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
42 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
43 : copy_fm_to_dbcsr,&
44 : cp_dbcsr_plus_fm_fm_t,&
45 : dbcsr_allocate_matrix_set,&
46 : dbcsr_deallocate_matrix_set
47 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
48 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
49 : cp_fm_scale,&
50 : cp_fm_scale_and_add,&
51 : cp_fm_schur_product,&
52 : cp_fm_uplo_to_full
53 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
54 : cp_fm_cholesky_invert,&
55 : cp_fm_cholesky_reduce,&
56 : cp_fm_cholesky_restore
57 : USE cp_fm_diag, ONLY: cp_fm_syevd
58 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
59 : cp_fm_struct_release,&
60 : cp_fm_struct_type
61 : USE cp_fm_types, ONLY: &
62 : copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, &
63 : cp_fm_get_info, cp_fm_release, cp_fm_set_all, cp_fm_set_element, cp_fm_start_copy_general, &
64 : cp_fm_to_fm, cp_fm_type
65 : USE cp_log_handling, ONLY: cp_get_default_logger,&
66 : cp_logger_type,&
67 : cp_to_string
68 : USE cp_output_handling, ONLY: cp_p_file,&
69 : cp_print_key_finished_output,&
70 : cp_print_key_should_output,&
71 : cp_print_key_unit_nr
72 : USE input_constants, ONLY: do_admm_purify_cauchy,&
73 : do_admm_purify_cauchy_subspace,&
74 : do_admm_purify_mo_diag,&
75 : do_admm_purify_mo_no_diag,&
76 : do_admm_purify_none
77 : USE input_section_types, ONLY: section_vals_type,&
78 : section_vals_val_get
79 : USE kinds, ONLY: default_string_length,&
80 : dp
81 : USE kpoint_methods, ONLY: kpoint_density_matrices,&
82 : kpoint_density_transform,&
83 : rskp_transform
84 : USE kpoint_types, ONLY: get_kpoint_env,&
85 : get_kpoint_info,&
86 : kpoint_env_type,&
87 : kpoint_type
88 : USE mathconstants, ONLY: z_one,&
89 : z_zero
90 : USE message_passing, ONLY: mp_para_env_type
91 : USE parallel_gemm_api, ONLY: parallel_gemm
92 : USE pw_types, ONLY: pw_c1d_gs_type,&
93 : pw_r3d_rs_type
94 : USE qs_collocate_density, ONLY: calculate_rho_elec
95 : USE qs_energy_types, ONLY: qs_energy_type
96 : USE qs_environment_types, ONLY: get_qs_env,&
97 : qs_environment_type
98 : USE qs_force_types, ONLY: add_qs_force,&
99 : qs_force_type
100 : USE qs_gapw_densities, ONLY: prepare_gapw_den
101 : USE qs_ks_atom, ONLY: update_ks_atom
102 : USE qs_ks_types, ONLY: qs_ks_env_type
103 : USE qs_local_rho_types, ONLY: local_rho_set_create,&
104 : local_rho_set_release,&
105 : local_rho_type
106 : USE qs_mo_types, ONLY: get_mo_set,&
107 : mo_set_type
108 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
109 : USE qs_overlap, ONLY: build_overlap_force
110 : USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,&
111 : calculate_rho_atom_coeff
112 : USE qs_rho_types, ONLY: qs_rho_get,&
113 : qs_rho_set,&
114 : qs_rho_type
115 : USE qs_scf_types, ONLY: qs_scf_env_type
116 : USE qs_vxc, ONLY: qs_vxc_create
117 : USE qs_vxc_atom, ONLY: calculate_vxc_atom
118 : USE task_list_types, ONLY: task_list_type
119 : #include "./base/base_uses.f90"
120 :
121 : IMPLICIT NONE
122 : PRIVATE
123 :
124 : PUBLIC :: admm_mo_calc_rho_aux, &
125 : admm_mo_calc_rho_aux_kp, &
126 : admm_mo_merge_ks_matrix, &
127 : admm_mo_merge_derivs, &
128 : admm_aux_response_density, &
129 : calc_mixed_overlap_force, &
130 : scale_dm, &
131 : admm_fit_mo_coeffs, &
132 : admm_update_ks_atom, &
133 : calc_admm_mo_derivatives, &
134 : calc_admm_ovlp_forces, &
135 : calc_admm_ovlp_forces_kp, &
136 : admm_projection_derivative, &
137 : kpoint_calc_admm_matrices
138 :
139 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'admm_methods'
140 :
141 : CONTAINS
142 :
143 : ! **************************************************************************************************
144 : !> \brief ...
145 : !> \param qs_env ...
146 : ! **************************************************************************************************
147 13160 : SUBROUTINE admm_mo_calc_rho_aux(qs_env)
148 : TYPE(qs_environment_type), POINTER :: qs_env
149 :
150 : CHARACTER(len=*), PARAMETER :: routineN = 'admm_mo_calc_rho_aux'
151 :
152 : CHARACTER(LEN=default_string_length) :: basis_type
153 : INTEGER :: handle, ispin
154 : LOGICAL :: gapw, s_mstruct_changed
155 13160 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r_aux
156 : TYPE(admm_type), POINTER :: admm_env
157 13160 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_aux_fit, &
158 13160 : matrix_s_aux_fit_vs_orb, rho_ao, &
159 13160 : rho_ao_aux
160 : TYPE(dft_control_type), POINTER :: dft_control
161 13160 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
162 : TYPE(mp_para_env_type), POINTER :: para_env
163 13160 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
164 13160 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_aux
165 : TYPE(qs_ks_env_type), POINTER :: ks_env
166 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
167 : TYPE(task_list_type), POINTER :: task_list
168 :
169 13160 : CALL timeset(routineN, handle)
170 :
171 13160 : NULLIFY (ks_env, admm_env, mos, mos_aux_fit, matrix_s_aux_fit, &
172 13160 : matrix_s_aux_fit_vs_orb, matrix_s, rho, rho_aux_fit, para_env)
173 13160 : NULLIFY (rho_g_aux, rho_r_aux, rho_ao, rho_ao_aux, tot_rho_r_aux, task_list)
174 :
175 : CALL get_qs_env(qs_env, &
176 : ks_env=ks_env, &
177 : admm_env=admm_env, &
178 : dft_control=dft_control, &
179 : mos=mos, &
180 : matrix_s=matrix_s, &
181 : para_env=para_env, &
182 : s_mstruct_changed=s_mstruct_changed, &
183 13160 : rho=rho)
184 : CALL get_admm_env(admm_env, mos_aux_fit=mos_aux_fit, matrix_s_aux_fit=matrix_s_aux_fit, &
185 13160 : matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, rho_aux_fit=rho_aux_fit)
186 :
187 13160 : CALL qs_rho_get(rho, rho_ao=rho_ao)
188 : CALL qs_rho_get(rho_aux_fit, &
189 : rho_ao=rho_ao_aux, &
190 : rho_g=rho_g_aux, &
191 : rho_r=rho_r_aux, &
192 13160 : tot_rho_r=tot_rho_r_aux)
193 :
194 13160 : gapw = admm_env%do_gapw
195 :
196 : ! convert mos from full to dbcsr matrices
197 28754 : DO ispin = 1, dft_control%nspins
198 28754 : IF (mos(ispin)%use_mo_coeff_b) THEN
199 9640 : CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, mos(ispin)%mo_coeff)
200 : END IF
201 : END DO
202 :
203 : ! fit mo coeffcients
204 : CALL admm_fit_mo_coeffs(admm_env, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, &
205 13160 : mos, mos_aux_fit, s_mstruct_changed)
206 :
207 28754 : DO ispin = 1, dft_control%nspins
208 15594 : IF (admm_env%block_dm) THEN
209 : CALL blockify_density_matrix(admm_env, &
210 : density_matrix=rho_ao(ispin)%matrix, &
211 : density_matrix_aux=rho_ao_aux(ispin)%matrix, &
212 : ispin=ispin, &
213 354 : nspins=dft_control%nspins)
214 :
215 : ELSE
216 :
217 : ! Here, the auxiliary DM gets calculated and is written into rho_aux_fit%...
218 : CALL calculate_dm_mo_no_diag(admm_env, &
219 : mo_set=mos(ispin), &
220 : overlap_matrix=matrix_s_aux_fit(1)%matrix, &
221 : density_matrix=rho_ao_aux(ispin)%matrix, &
222 : overlap_matrix_large=matrix_s(1)%matrix, &
223 : density_matrix_large=rho_ao(ispin)%matrix, &
224 15240 : ispin=ispin)
225 :
226 : END IF
227 :
228 15594 : IF (admm_env%purification_method == do_admm_purify_cauchy) THEN
229 : CALL purify_dm_cauchy(admm_env, &
230 : mo_set=mos_aux_fit(ispin), &
231 : density_matrix=rho_ao_aux(ispin)%matrix, &
232 : ispin=ispin, &
233 484 : blocked=admm_env%block_dm)
234 : END IF
235 :
236 : !GPW is the default, PW density is computed using the AUX_FIT basis and task_list
237 : !If GAPW, the we use the AUX_FIT_SOFT basis and task list
238 15594 : basis_type = "AUX_FIT"
239 15594 : task_list => admm_env%task_list_aux_fit
240 15594 : IF (gapw) THEN
241 5594 : basis_type = "AUX_FIT_SOFT"
242 5594 : task_list => admm_env%admm_gapw_env%task_list
243 : END IF
244 :
245 : CALL calculate_rho_elec(ks_env=ks_env, &
246 : matrix_p=rho_ao_aux(ispin)%matrix, &
247 : rho=rho_r_aux(ispin), &
248 : rho_gspace=rho_g_aux(ispin), &
249 : total_rho=tot_rho_r_aux(ispin), &
250 : soft_valid=.FALSE., &
251 : basis_type=basis_type, &
252 28754 : task_list_external=task_list)
253 :
254 : END DO
255 :
256 : !If GAPW, also need to prepare the atomic densities
257 13160 : IF (gapw) THEN
258 :
259 : CALL calculate_rho_atom_coeff(qs_env, rho_ao_aux, &
260 : rho_atom_set=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
261 : qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
262 4848 : oce=admm_env%admm_gapw_env%oce, sab=admm_env%sab_aux_fit, para_env=para_env)
263 :
264 : CALL prepare_gapw_den(qs_env, local_rho_set=admm_env%admm_gapw_env%local_rho_set, &
265 4848 : do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
266 : END IF
267 :
268 13160 : IF (dft_control%nspins == 1) THEN
269 10726 : admm_env%gsi(3) = admm_env%gsi(1)
270 : ELSE
271 2434 : admm_env%gsi(3) = (admm_env%gsi(1) + admm_env%gsi(2))/2.0_dp
272 : END IF
273 :
274 13160 : CALL qs_rho_set(rho_aux_fit, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
275 :
276 13160 : CALL timestop(handle)
277 :
278 13160 : END SUBROUTINE admm_mo_calc_rho_aux
279 : ! **************************************************************************************************
280 : !> \brief ...
281 : !> \param qs_env ...
282 : ! **************************************************************************************************
283 242 : SUBROUTINE admm_mo_calc_rho_aux_kp(qs_env)
284 : TYPE(qs_environment_type), POINTER :: qs_env
285 :
286 : CHARACTER(len=*), PARAMETER :: routineN = 'admm_mo_calc_rho_aux_kp'
287 :
288 : CHARACTER(LEN=default_string_length) :: basis_type
289 : INTEGER :: handle, igroup, ik, ikp, img, indx, &
290 : ispin, kplocal, kpmax, nao_aux_fit, &
291 : nao_orb, natom, nkp, nkp_groups, nmo, &
292 : nspins
293 : INTEGER, DIMENSION(2) :: kp_range
294 242 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
295 242 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
296 : LOGICAL :: gapw, my_kpgrp, pmat_from_rs, &
297 : use_real_wfn
298 : REAL(dp) :: maxval_mos, nelec_aux(2), nelec_orb(2), &
299 : tmp
300 242 : REAL(KIND=dp), DIMENSION(:), POINTER :: occ_num, occ_num_aux, tot_rho_r_aux
301 242 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
302 : TYPE(admm_type), POINTER :: admm_env
303 242 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
304 : TYPE(cp_cfm_type) :: cpmatrix, cwork_aux_aux, cwork_aux_orb
305 : TYPE(cp_fm_struct_type), POINTER :: struct_aux_aux, struct_aux_orb, &
306 : struct_orb_orb
307 : TYPE(cp_fm_type) :: fmdummy, work_aux_orb, work_orb_orb, &
308 : work_orb_orb2
309 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
310 242 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
311 242 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_s_aux_fit, rho_ao_aux, &
312 242 : rho_ao_orb
313 : TYPE(dbcsr_type) :: pmatrix_tmp
314 242 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: pmatrix
315 : TYPE(dft_control_type), POINTER :: dft_control
316 : TYPE(kpoint_env_type), POINTER :: kp
317 : TYPE(kpoint_type), POINTER :: kpoints
318 242 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit, mos_aux_fit_kp, mos_kp
319 : TYPE(mp_para_env_type), POINTER :: para_env
320 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
321 242 : POINTER :: sab_aux_fit, sab_kp
322 242 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
323 242 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_aux
324 : TYPE(qs_ks_env_type), POINTER :: ks_env
325 : TYPE(qs_rho_type), POINTER :: rho_aux_fit, rho_orb
326 : TYPE(qs_scf_env_type), POINTER :: scf_env
327 : TYPE(task_list_type), POINTER :: task_list
328 :
329 242 : CALL timeset(routineN, handle)
330 :
331 242 : NULLIFY (ks_env, admm_env, mos, mos_aux_fit, matrix_s, rho_orb, &
332 242 : matrix_s_aux_fit, rho_aux_fit, rho_ao_orb, &
333 242 : para_env, rho_g_aux, rho_r_aux, rho_ao_aux, tot_rho_r_aux, &
334 242 : kpoints, sab_aux_fit, sab_kp, kp, &
335 242 : struct_orb_orb, struct_aux_orb, struct_aux_aux)
336 :
337 : CALL get_qs_env(qs_env, &
338 : ks_env=ks_env, &
339 : admm_env=admm_env, &
340 : dft_control=dft_control, &
341 : kpoints=kpoints, &
342 : natom=natom, &
343 : scf_env=scf_env, &
344 : matrix_s_kp=matrix_s, &
345 242 : rho=rho_orb)
346 : CALL get_admm_env(admm_env, &
347 : rho_aux_fit=rho_aux_fit, &
348 : matrix_s_aux_fit_kp=matrix_s_aux_fit, &
349 242 : sab_aux_fit=sab_aux_fit)
350 242 : gapw = admm_env%do_gapw
351 :
352 : CALL qs_rho_get(rho_aux_fit, &
353 : rho_ao_kp=rho_ao_aux, &
354 : rho_g=rho_g_aux, &
355 : rho_r=rho_r_aux, &
356 242 : tot_rho_r=tot_rho_r_aux)
357 :
358 242 : CALL qs_rho_get(rho_orb, rho_ao_kp=rho_ao_orb)
359 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
360 : nkp_groups=nkp_groups, kp_dist=kp_dist, &
361 242 : cell_to_index=cell_to_index, sab_nl=sab_kp)
362 :
363 : ! the temporary DBCSR matrices for the rskp_transform we have to manually allocate
364 : ! index 1 => real, index 2 => imaginary
365 726 : ALLOCATE (pmatrix(2))
366 : CALL dbcsr_create(pmatrix(1), template=matrix_s(1, 1)%matrix, &
367 242 : matrix_type=dbcsr_type_symmetric)
368 : CALL dbcsr_create(pmatrix(2), template=matrix_s(1, 1)%matrix, &
369 242 : matrix_type=dbcsr_type_antisymmetric)
370 : CALL dbcsr_create(pmatrix_tmp, template=matrix_s(1, 1)%matrix, &
371 242 : matrix_type=dbcsr_type_no_symmetry)
372 242 : CALL cp_dbcsr_alloc_block_from_nbl(pmatrix(1), sab_kp)
373 242 : CALL cp_dbcsr_alloc_block_from_nbl(pmatrix(2), sab_kp)
374 :
375 242 : nao_aux_fit = admm_env%nao_aux_fit
376 242 : nao_orb = admm_env%nao_orb
377 242 : nspins = dft_control%nspins
378 :
379 : !Create fm and cfm work matrices, for each KP subgroup
380 : CALL cp_fm_struct_create(struct_orb_orb, context=kpoints%blacs_env, para_env=kpoints%para_env_kp, &
381 242 : nrow_global=nao_orb, ncol_global=nao_orb)
382 242 : CALL cp_fm_create(work_orb_orb, struct_orb_orb)
383 242 : CALL cp_fm_create(work_orb_orb2, struct_orb_orb)
384 :
385 : CALL cp_fm_struct_create(struct_aux_aux, context=kpoints%blacs_env, para_env=kpoints%para_env_kp, &
386 242 : nrow_global=nao_aux_fit, ncol_global=nao_aux_fit)
387 :
388 : CALL cp_fm_struct_create(struct_aux_orb, context=kpoints%blacs_env, para_env=kpoints%para_env_kp, &
389 242 : nrow_global=nao_aux_fit, ncol_global=nao_orb)
390 242 : CALL cp_fm_create(work_aux_orb, struct_orb_orb)
391 :
392 242 : IF (.NOT. use_real_wfn) THEN
393 242 : CALL cp_cfm_create(cpmatrix, struct_orb_orb)
394 :
395 242 : CALL cp_cfm_create(cwork_aux_aux, struct_aux_aux)
396 :
397 242 : CALL cp_cfm_create(cwork_aux_orb, struct_aux_orb)
398 :
399 : END IF
400 :
401 242 : CALL cp_fm_struct_release(struct_orb_orb)
402 242 : CALL cp_fm_struct_release(struct_aux_aux)
403 242 : CALL cp_fm_struct_release(struct_aux_orb)
404 :
405 242 : para_env => kpoints%blacs_env_all%para_env
406 242 : kplocal = kp_range(2) - kp_range(1) + 1
407 660 : kpmax = MAXVAL(kp_dist(2, :) - kp_dist(1, :) + 1)
408 :
409 : !We querry the maximum absolute value of the KP MOs to see if they are populated at all. If not, we
410 : !need to get the KP Pmat from the RS ones (happens at first SCF step, for example)
411 242 : maxval_mos = 0.0_dp
412 1554 : DO ikp = 1, kplocal
413 1312 : CALL get_kpoint_env(kpoints%kp_env(ikp)%kpoint_env, mos=mos_kp)
414 2989 : DO ispin = 1, nspins
415 2747 : IF (use_real_wfn) THEN
416 0 : maxval_mos = MAX(maxval_mos, MAXVAL(ABS(mos_kp(ispin)%mo_coeff%local_data)))
417 : ELSE
418 : maxval_mos = MAX(maxval_mos, MAXVAL(ABS(REAL(mos_kp(ispin)%cmo_coeff%local_data, KIND=dp))), &
419 167662 : MAXVAL(ABS(AIMAG(mos_kp(ispin)%cmo_coeff%local_data))))
420 : END IF
421 : END DO
422 : END DO
423 242 : CALL para_env%sum(maxval_mos) !I think para_env is the global one
424 :
425 242 : pmat_from_rs = .FALSE.
426 242 : IF (maxval_mos < EPSILON(0.0_dp)) pmat_from_rs = .TRUE.
427 :
428 : !TODO: issue a warning when doing ADMM with ATOMIC guess. If small number of K-points => leads to bad things
429 :
430 8754 : ALLOCATE (info(nkp*nspins, 2))
431 : !Start communication: only P matrix, and only if required
432 242 : indx = 0
433 242 : IF (pmat_from_rs) THEN
434 216 : DO ikp = 1, kpmax
435 436 : DO ispin = 1, nspins
436 838 : DO igroup = 1, nkp_groups
437 : ! number of current kpoint
438 438 : ik = kp_dist(1, igroup) + ikp - 1
439 438 : IF (ik > kp_dist(2, igroup)) CYCLE
440 424 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
441 424 : indx = indx + 1
442 :
443 : ! FT of matrices P if required, then transfer to FM type
444 424 : IF (use_real_wfn) THEN
445 0 : CALL dbcsr_set(pmatrix(1), 0.0_dp)
446 : CALL rskp_transform(rmatrix=pmatrix(1), rsmat=rho_ao_orb, ispin=ispin, &
447 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_kp)
448 0 : CALL dbcsr_desymmetrize(pmatrix(1), pmatrix_tmp)
449 0 : CALL copy_dbcsr_to_fm(pmatrix_tmp, admm_env%work_orb_orb)
450 : ELSE
451 424 : CALL dbcsr_set(pmatrix(1), 0.0_dp)
452 424 : CALL dbcsr_set(pmatrix(2), 0.0_dp)
453 : CALL rskp_transform(rmatrix=pmatrix(1), cmatrix=pmatrix(2), rsmat=rho_ao_orb, ispin=ispin, &
454 424 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_kp)
455 424 : CALL dbcsr_desymmetrize(pmatrix(1), pmatrix_tmp)
456 424 : CALL copy_dbcsr_to_fm(pmatrix_tmp, admm_env%work_orb_orb)
457 424 : CALL dbcsr_desymmetrize(pmatrix(2), pmatrix_tmp)
458 424 : CALL copy_dbcsr_to_fm(pmatrix_tmp, admm_env%work_orb_orb2)
459 : END IF
460 :
461 644 : IF (my_kpgrp) THEN
462 213 : CALL cp_fm_start_copy_general(admm_env%work_orb_orb, work_orb_orb, para_env, info(indx, 1))
463 213 : IF (.NOT. use_real_wfn) THEN
464 213 : CALL cp_fm_start_copy_general(admm_env%work_orb_orb2, work_orb_orb2, para_env, info(indx, 2))
465 : END IF
466 : ELSE
467 211 : CALL cp_fm_start_copy_general(admm_env%work_orb_orb, fmdummy, para_env, info(indx, 1))
468 211 : IF (.NOT. use_real_wfn) THEN
469 211 : CALL cp_fm_start_copy_general(admm_env%work_orb_orb2, fmdummy, para_env, info(indx, 2))
470 : END IF
471 : END IF !my_kpgrp
472 : END DO
473 : END DO
474 : END DO
475 : END IF !pmat_from_rs
476 :
477 : indx = 0
478 1594 : DO ikp = 1, kpmax
479 3074 : DO ispin = 1, nspins
480 4374 : DO igroup = 1, nkp_groups
481 : ! number of current kpoint
482 2894 : ik = kp_dist(1, igroup) + ikp - 1
483 2894 : IF (ik > kp_dist(2, igroup)) CYCLE
484 2804 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
485 2804 : indx = indx + 1
486 4284 : IF (my_kpgrp .AND. pmat_from_rs) THEN
487 213 : CALL cp_fm_finish_copy_general(work_orb_orb, info(indx, 1))
488 213 : IF (.NOT. use_real_wfn) THEN
489 213 : CALL cp_fm_finish_copy_general(work_orb_orb2, info(indx, 2))
490 213 : CALL cp_fm_to_cfm(work_orb_orb, work_orb_orb2, cpmatrix)
491 : END IF
492 : END IF
493 : END DO
494 :
495 1480 : IF (ikp > kplocal) CYCLE
496 2787 : IF (use_real_wfn) THEN
497 :
498 0 : nmo = admm_env%nmo(ispin)
499 : !! Each kpoint group has now information on a kpoint for which to calculate the MOS_aux
500 0 : CALL get_kpoint_env(kpoints%kp_env(ikp)%kpoint_env, mos=mos_kp)
501 0 : CALL get_kpoint_env(kpoints%kp_aux_env(ikp)%kpoint_env, mos=mos_aux_fit_kp)
502 0 : mos => mos_kp
503 0 : mos_aux_fit => mos_aux_fit_kp
504 :
505 0 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_num)
506 : CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit, &
507 0 : occupation_numbers=occ_num_aux)
508 :
509 0 : kp => kpoints%kp_aux_env(ikp)%kpoint_env
510 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, 1.0_dp, kp%amat, &
511 0 : mo_coeff, 0.0_dp, mo_coeff_aux_fit)
512 :
513 0 : occ_num_aux(1:nmo) = occ_num(1:nmo)
514 :
515 0 : IF (pmat_from_rs) THEN
516 : !We project on the AUX basis: P_aux = A * P *A^T
517 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, 1.0_dp, kp%amat, &
518 0 : work_orb_orb, 0.0_dp, work_aux_orb)
519 : CALL parallel_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nao_orb, 1.0_dp, work_aux_orb, &
520 0 : kp%amat, 0.0_dp, kpoints%kp_aux_env(ikp)%kpoint_env%pmat(1, ispin))
521 : END IF
522 :
523 : ELSE !complex wfn
524 :
525 1435 : nmo = admm_env%nmo(ispin)
526 1435 : CALL get_kpoint_env(kpoints%kp_env(ikp)%kpoint_env, mos=mos_kp)
527 :
528 : !project
529 1435 : kp => kpoints%kp_aux_env(ikp)%kpoint_env
530 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
531 : z_one, kp%camat, mos_kp(ispin)%cmo_coeff, z_zero, &
532 1435 : kp%mos(ispin)%cmo_coeff)
533 :
534 : !write result back to KP MOs
535 1435 : CALL get_kpoint_env(kpoints%kp_aux_env(ikp)%kpoint_env, mos=mos_aux_fit_kp)
536 1435 : CALL get_mo_set(mos_kp(ispin), occupation_numbers=occ_num)
537 1435 : CALL get_mo_set(mos_aux_fit_kp(ispin), occupation_numbers=occ_num_aux)
538 9546 : occ_num_aux(:) = occ_num(:)
539 :
540 1435 : IF (pmat_from_rs) THEN
541 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, z_one, kp%camat, &
542 213 : cpmatrix, z_zero, cwork_aux_orb)
543 : CALL parallel_gemm('N', 'C', nao_aux_fit, nao_aux_fit, nao_orb, z_one, cwork_aux_orb, &
544 213 : kp%camat, z_zero, cwork_aux_aux)
545 :
546 : CALL cp_cfm_to_fm(cwork_aux_aux, mtargetr=kpoints%kp_aux_env(ikp)%kpoint_env%pmat(1, ispin), &
547 213 : mtargeti=kpoints%kp_aux_env(ikp)%kpoint_env%pmat(2, ispin))
548 : END IF
549 : END IF
550 :
551 : END DO
552 : END DO
553 :
554 : ! Symmetry is handled by the reduced-grid transform: one P_aux(k) per
555 : ! irreducible k-point is expanded with the stored operations and weighted
556 : ! by wkp, instead of computing every symmetry-equivalent k-point.
557 242 : IF (.NOT. pmat_from_rs) CALL kpoint_density_matrices(kpoints, for_aux_fit=.TRUE.)
558 : CALL kpoint_density_transform(kpoints, rho_ao_aux, .FALSE., &
559 : matrix_s_aux_fit(1, 1)%matrix, sab_aux_fit, &
560 : admm_env%scf_work_aux_fit, for_aux_fit=.TRUE., &
561 242 : overlap_rs=matrix_s_aux_fit)
562 :
563 : !Clean-up communication
564 242 : IF (pmat_from_rs) THEN
565 460 : DO indx = 1, SIZE(info, 1)
566 424 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
567 460 : IF (.NOT. use_real_wfn) CALL cp_fm_cleanup_copy_general(info(indx, 2))
568 : END DO
569 : END IF
570 :
571 6092 : DEALLOCATE (info)
572 242 : CALL dbcsr_release(pmatrix(1))
573 242 : CALL dbcsr_release(pmatrix(2))
574 242 : CALL dbcsr_release(pmatrix_tmp)
575 :
576 242 : CALL cp_fm_release(work_orb_orb)
577 242 : CALL cp_fm_release(work_orb_orb2)
578 242 : CALL cp_fm_release(work_aux_orb)
579 242 : IF (.NOT. use_real_wfn) THEN
580 242 : CALL cp_cfm_release(cpmatrix)
581 242 : CALL cp_cfm_release(cwork_aux_aux)
582 242 : CALL cp_cfm_release(cwork_aux_orb)
583 : END IF
584 :
585 : !ADMMQ, ADMMP, ADMMS
586 242 : IF (admm_env%do_admmq .OR. admm_env%do_admmp .OR. admm_env%do_admms) THEN
587 :
588 94 : CALL cite_reference(Merlot2014)
589 :
590 94 : nelec_orb = 0.0_dp
591 94 : nelec_aux = 0.0_dp
592 376 : admm_env%n_large_basis = 0.0_dp
593 : !Note: we can take the trace of the symmetric-typed matrices as P_mu^0,nu^b = P_nu^0,mu^-b
594 : ! and because of the sum over all images, all atomic blocks are accounted for
595 3680 : DO img = 1, dft_control%nimages
596 7496 : DO ispin = 1, dft_control%nspins
597 3816 : CALL dbcsr_dot(rho_ao_orb(ispin, img)%matrix, matrix_s(1, img)%matrix, tmp)
598 3816 : nelec_orb(ispin) = nelec_orb(ispin) + tmp
599 3816 : CALL dbcsr_dot(rho_ao_aux(ispin, img)%matrix, matrix_s_aux_fit(1, img)%matrix, tmp)
600 7402 : nelec_aux(ispin) = nelec_aux(ispin) + tmp
601 : END DO
602 : END DO
603 :
604 202 : DO ispin = 1, dft_control%nspins
605 108 : admm_env%n_large_basis(ispin) = nelec_orb(ispin)
606 202 : admm_env%gsi(ispin) = nelec_orb(ispin)/nelec_aux(ispin)
607 : END DO
608 :
609 94 : IF (admm_env%charge_constrain) THEN
610 3102 : DO img = 1, dft_control%nimages
611 6354 : DO ispin = 1, dft_control%nspins
612 6274 : CALL dbcsr_scale(rho_ao_aux(ispin, img)%matrix, admm_env%gsi(ispin))
613 : END DO
614 : END DO
615 : END IF
616 :
617 94 : IF (dft_control%nspins == 1) THEN
618 80 : admm_env%gsi(3) = admm_env%gsi(1)
619 : ELSE
620 14 : admm_env%gsi(3) = (admm_env%gsi(1) + admm_env%gsi(2))/2.0_dp
621 : END IF
622 : END IF
623 :
624 242 : basis_type = "AUX_FIT"
625 242 : task_list => admm_env%task_list_aux_fit
626 242 : IF (gapw) THEN
627 84 : basis_type = "AUX_FIT_SOFT"
628 84 : task_list => admm_env%admm_gapw_env%task_list
629 : END IF
630 :
631 508 : DO ispin = 1, nspins
632 266 : rho_ao => rho_ao_aux(ispin, :)
633 : CALL calculate_rho_elec(ks_env=ks_env, &
634 : matrix_p_kp=rho_ao, &
635 : rho=rho_r_aux(ispin), &
636 : rho_gspace=rho_g_aux(ispin), &
637 : total_rho=tot_rho_r_aux(ispin), &
638 : soft_valid=.FALSE., &
639 : basis_type=basis_type, &
640 508 : task_list_external=task_list)
641 : END DO
642 :
643 242 : IF (gapw) THEN
644 : CALL calculate_rho_atom_coeff(qs_env, rho_ao_aux, &
645 : rho_atom_set=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
646 : qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
647 : oce=admm_env%admm_gapw_env%oce, &
648 84 : sab=admm_env%sab_aux_fit, para_env=para_env)
649 :
650 : CALL prepare_gapw_den(qs_env, local_rho_set=admm_env%admm_gapw_env%local_rho_set, &
651 84 : do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
652 : END IF
653 :
654 242 : CALL qs_rho_set(rho_aux_fit, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
655 :
656 242 : CALL timestop(handle)
657 :
658 968 : END SUBROUTINE admm_mo_calc_rho_aux_kp
659 :
660 : ! **************************************************************************************************
661 : !> \brief Adds the GAPW exchange contribution to the aux_fit ks matrices
662 : !> \param qs_env ...
663 : !> \param calculate_forces ...
664 : ! **************************************************************************************************
665 4970 : SUBROUTINE admm_update_ks_atom(qs_env, calculate_forces)
666 :
667 : TYPE(qs_environment_type), POINTER :: qs_env
668 : LOGICAL, INTENT(IN) :: calculate_forces
669 :
670 : CHARACTER(len=*), PARAMETER :: routineN = 'admm_update_ks_atom'
671 :
672 : INTEGER :: handle, img, ispin
673 : REAL(dp) :: force_fac(2)
674 : TYPE(admm_type), POINTER :: admm_env
675 4970 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit, &
676 4970 : matrix_ks_aux_fit_dft, &
677 4970 : matrix_ks_aux_fit_hfx, rho_ao_aux
678 : TYPE(dft_control_type), POINTER :: dft_control
679 : TYPE(qs_rho_type), POINTER :: rho_aux_fit
680 :
681 4970 : NULLIFY (matrix_ks_aux_fit, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, rho_ao_aux, rho_aux_fit)
682 4970 : NULLIFY (admm_env, dft_control)
683 :
684 4970 : CALL timeset(routineN, handle)
685 :
686 4970 : CALL get_qs_env(qs_env, admm_env=admm_env, dft_control=dft_control)
687 : CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, matrix_ks_aux_fit_kp=matrix_ks_aux_fit, &
688 : matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft, &
689 4970 : matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx)
690 4970 : CALL qs_rho_get(rho_aux_fit, rho_ao_kp=rho_ao_aux)
691 :
692 : !In case of ADMMS or ADMMP, need to scale the forces stemming from DFT exchagne correction
693 14910 : force_fac = 1.0_dp
694 4970 : IF (admm_env%do_admms) THEN
695 298 : DO ispin = 1, dft_control%nspins
696 298 : force_fac(ispin) = admm_env%gsi(ispin)**(2.0_dp/3.0_dp)
697 : END DO
698 4852 : ELSE IF (admm_env%do_admmp) THEN
699 756 : DO ispin = 1, dft_control%nspins
700 756 : force_fac(ispin) = admm_env%gsi(ispin)**2
701 : END DO
702 : END IF
703 :
704 : CALL update_ks_atom(qs_env, matrix_ks_aux_fit, rho_ao_aux, calculate_forces, tddft=.FALSE., &
705 : rho_atom_external=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
706 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
707 : oce_external=admm_env%admm_gapw_env%oce, &
708 4970 : sab_external=admm_env%sab_aux_fit, fscale=force_fac)
709 :
710 : !Following the logic of sum_up_and_integrate to recover the pure DFT exchange contribution
711 13144 : DO img = 1, dft_control%nimages
712 22316 : DO ispin = 1, dft_control%nspins
713 : CALL dbcsr_add(matrix_ks_aux_fit_dft(ispin, img)%matrix, matrix_ks_aux_fit(ispin, img)%matrix, &
714 9172 : 0.0_dp, -1.0_dp)
715 : CALL dbcsr_add(matrix_ks_aux_fit_dft(ispin, img)%matrix, matrix_ks_aux_fit_hfx(ispin, img)%matrix, &
716 17346 : 1.0_dp, 1.0_dp)
717 : END DO
718 : END DO
719 :
720 4970 : CALL timestop(handle)
721 :
722 4970 : END SUBROUTINE admm_update_ks_atom
723 :
724 : ! **************************************************************************************************
725 : !> \brief ...
726 : !> \param qs_env ...
727 : ! **************************************************************************************************
728 13390 : SUBROUTINE admm_mo_merge_ks_matrix(qs_env)
729 : TYPE(qs_environment_type), POINTER :: qs_env
730 :
731 : CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_mo_merge_ks_matrix'
732 :
733 : INTEGER :: handle
734 : TYPE(admm_type), POINTER :: admm_env
735 : TYPE(dft_control_type), POINTER :: dft_control
736 :
737 13390 : CALL timeset(routineN, handle)
738 13390 : NULLIFY (admm_env)
739 :
740 13390 : CALL get_qs_env(qs_env, admm_env=admm_env, dft_control=dft_control)
741 :
742 13680 : SELECT CASE (admm_env%purification_method)
743 : CASE (do_admm_purify_cauchy)
744 290 : CALL merge_ks_matrix_cauchy(qs_env)
745 :
746 : CASE (do_admm_purify_cauchy_subspace)
747 154 : CALL merge_ks_matrix_cauchy_subspace(qs_env)
748 :
749 : CASE (do_admm_purify_none)
750 11204 : IF (dft_control%nimages > 1) THEN
751 242 : CALL merge_ks_matrix_none_kp(qs_env)
752 : ELSE
753 10962 : CALL merge_ks_matrix_none(qs_env)
754 : END IF
755 :
756 : CASE (do_admm_purify_mo_diag, do_admm_purify_mo_no_diag)
757 : !do nothing
758 : CASE DEFAULT
759 13390 : CPABORT("admm_mo_merge_ks_matrix: unknown purification method")
760 : END SELECT
761 :
762 13390 : CALL timestop(handle)
763 :
764 13390 : END SUBROUTINE admm_mo_merge_ks_matrix
765 :
766 : ! **************************************************************************************************
767 : !> \brief ...
768 : !> \param ispin ...
769 : !> \param admm_env ...
770 : !> \param mo_set ...
771 : !> \param mo_coeff ...
772 : !> \param mo_coeff_aux_fit ...
773 : !> \param mo_derivs ...
774 : !> \param mo_derivs_aux_fit ...
775 : !> \param matrix_ks_aux_fit ...
776 : ! **************************************************************************************************
777 8058 : SUBROUTINE admm_mo_merge_derivs(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, mo_derivs, &
778 8058 : mo_derivs_aux_fit, matrix_ks_aux_fit)
779 : INTEGER, INTENT(IN) :: ispin
780 : TYPE(admm_type), POINTER :: admm_env
781 : TYPE(mo_set_type), INTENT(IN) :: mo_set
782 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff, mo_coeff_aux_fit
783 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_derivs, mo_derivs_aux_fit
784 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit
785 :
786 : CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_mo_merge_derivs'
787 :
788 : INTEGER :: handle
789 :
790 8058 : CALL timeset(routineN, handle)
791 :
792 9224 : SELECT CASE (admm_env%purification_method)
793 : CASE (do_admm_purify_mo_diag)
794 : CALL merge_mo_derivs_diag(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, &
795 1166 : mo_derivs, mo_derivs_aux_fit, matrix_ks_aux_fit)
796 :
797 : CASE (do_admm_purify_mo_no_diag)
798 100 : CALL merge_mo_derivs_no_diag(ispin, admm_env, mo_set, mo_derivs, matrix_ks_aux_fit)
799 :
800 : CASE (do_admm_purify_none, do_admm_purify_cauchy, do_admm_purify_cauchy_subspace)
801 : !do nothing
802 : CASE DEFAULT
803 8058 : CPABORT("admm_mo_merge_derivs: unknown purification method")
804 : END SELECT
805 :
806 8058 : CALL timestop(handle)
807 :
808 8058 : END SUBROUTINE admm_mo_merge_derivs
809 :
810 : ! **************************************************************************************************
811 : !> \brief ...
812 : !> \param admm_env ...
813 : !> \param matrix_s_aux_fit ...
814 : !> \param matrix_s_mixed ...
815 : !> \param mos ...
816 : !> \param mos_aux_fit ...
817 : !> \param geometry_did_change ...
818 : ! **************************************************************************************************
819 26328 : SUBROUTINE admm_fit_mo_coeffs(admm_env, matrix_s_aux_fit, matrix_s_mixed, &
820 13164 : mos, mos_aux_fit, geometry_did_change)
821 :
822 : TYPE(admm_type), POINTER :: admm_env
823 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_mixed
824 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mos_aux_fit
825 : LOGICAL, INTENT(IN) :: geometry_did_change
826 :
827 : CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_fit_mo_coeffs'
828 :
829 : INTEGER :: handle
830 :
831 13164 : CALL timeset(routineN, handle)
832 :
833 13164 : IF (geometry_did_change) THEN
834 924 : CALL fit_mo_coeffs(admm_env, matrix_s_aux_fit, matrix_s_mixed)
835 : END IF
836 :
837 13432 : SELECT CASE (admm_env%purification_method)
838 : CASE (do_admm_purify_mo_no_diag, do_admm_purify_cauchy_subspace)
839 268 : CALL purify_mo_cholesky(admm_env, mos, mos_aux_fit)
840 :
841 : CASE (do_admm_purify_mo_diag)
842 1628 : CALL purify_mo_diag(admm_env, mos, mos_aux_fit)
843 :
844 : CASE DEFAULT
845 13164 : CALL purify_mo_none(admm_env, mos, mos_aux_fit)
846 : END SELECT
847 :
848 13164 : CALL timestop(handle)
849 :
850 13164 : END SUBROUTINE admm_fit_mo_coeffs
851 :
852 : ! **************************************************************************************************
853 : !> \brief Calculate S^-1, Q, B full-matrices given sparse S_tilde and Q
854 : !> \param admm_env ...
855 : !> \param matrix_s_aux_fit ...
856 : !> \param matrix_s_mixed ...
857 : ! **************************************************************************************************
858 924 : SUBROUTINE fit_mo_coeffs(admm_env, matrix_s_aux_fit, matrix_s_mixed)
859 : TYPE(admm_type), POINTER :: admm_env
860 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_mixed
861 :
862 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fit_mo_coeffs'
863 :
864 : INTEGER :: handle, iatom, jatom, nao_aux_fit, &
865 : nao_orb
866 924 : REAL(dp), DIMENSION(:, :), POINTER :: sparse_block
867 : TYPE(dbcsr_iterator_type) :: iter
868 : TYPE(dbcsr_type), POINTER :: matrix_s_tilde
869 :
870 924 : CALL timeset(routineN, handle)
871 :
872 924 : nao_aux_fit = admm_env%nao_aux_fit
873 924 : nao_orb = admm_env%nao_orb
874 :
875 : ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed
876 :
877 924 : IF (.NOT. admm_env%block_fit) THEN
878 916 : CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix, admm_env%S_inv)
879 : ELSE
880 : NULLIFY (matrix_s_tilde)
881 8 : ALLOCATE (matrix_s_tilde)
882 : CALL dbcsr_create(matrix_s_tilde, template=matrix_s_aux_fit(1)%matrix, &
883 : name='MATRIX s_tilde', &
884 8 : matrix_type=dbcsr_type_symmetric)
885 :
886 8 : CALL dbcsr_copy(matrix_s_tilde, matrix_s_aux_fit(1)%matrix)
887 :
888 8 : CALL dbcsr_iterator_start(iter, matrix_s_tilde)
889 48 : DO WHILE (dbcsr_iterator_blocks_left(iter))
890 40 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block)
891 48 : IF (admm_env%block_map(iatom, jatom) == 0) THEN
892 102 : sparse_block = 0.0_dp
893 : END IF
894 : END DO
895 8 : CALL dbcsr_iterator_stop(iter)
896 8 : CALL copy_dbcsr_to_fm(matrix_s_tilde, admm_env%S_inv)
897 8 : CALL dbcsr_deallocate_matrix(matrix_s_tilde)
898 : END IF
899 :
900 924 : CALL cp_fm_uplo_to_full(admm_env%S_inv, admm_env%work_aux_aux)
901 924 : CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S)
902 :
903 924 : CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix, admm_env%Q)
904 :
905 : !! Calculate S'_inverse
906 924 : CALL cp_fm_cholesky_decompose(admm_env%S_inv)
907 924 : CALL cp_fm_cholesky_invert(admm_env%S_inv)
908 : !! Symmetrize the guy
909 924 : CALL cp_fm_uplo_to_full(admm_env%S_inv, admm_env%work_aux_aux)
910 :
911 : !! Calculate A=S'^(-1)*Q
912 924 : IF (admm_env%block_fit) THEN
913 8 : CALL cp_fm_set_all(admm_env%A, 0.0_dp, 1.0_dp)
914 : ELSE
915 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
916 : 1.0_dp, admm_env%S_inv, admm_env%Q, 0.0_dp, &
917 916 : admm_env%A)
918 :
919 : ! this multiplication is apparent not need for purify_none
920 : !! B=Q^(T)*A
921 : CALL parallel_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, &
922 : 1.0_dp, admm_env%Q, admm_env%A, 0.0_dp, &
923 916 : admm_env%B)
924 : END IF
925 :
926 924 : CALL timestop(handle)
927 :
928 924 : END SUBROUTINE fit_mo_coeffs
929 :
930 : ! **************************************************************************************************
931 : !> \brief Calculates the MO coefficients for the auxiliary fitting basis set
932 : !> by minimizing int (psi_i - psi_aux_i)^2 using Lagrangian Multipliers
933 : !>
934 : !> \param admm_env The ADMM env
935 : !> \param mos the MO's of the orbital basis set
936 : !> \param mos_aux_fit the MO's of the auxiliary fitting basis set
937 : !> \par History
938 : !> 05.2008 created [Manuel Guidon]
939 : !> \author Manuel Guidon
940 : ! **************************************************************************************************
941 268 : SUBROUTINE purify_mo_cholesky(admm_env, mos, mos_aux_fit)
942 :
943 : TYPE(admm_type), POINTER :: admm_env
944 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mos_aux_fit
945 :
946 : CHARACTER(LEN=*), PARAMETER :: routineN = 'purify_mo_cholesky'
947 :
948 : INTEGER :: handle, ispin, nao_aux_fit, nao_orb, &
949 : nmo, nspins
950 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
951 :
952 268 : CALL timeset(routineN, handle)
953 :
954 268 : nao_aux_fit = admm_env%nao_aux_fit
955 268 : nao_orb = admm_env%nao_orb
956 268 : nspins = SIZE(mos)
957 :
958 : ! *** Calculate the mo_coeffs for the fitting basis
959 670 : DO ispin = 1, nspins
960 402 : nmo = admm_env%nmo(ispin)
961 402 : IF (nmo == 0) CYCLE
962 : !! Lambda = C^(T)*B*C
963 402 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
964 402 : CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
965 : CALL parallel_gemm('N', 'N', nao_orb, nmo, nao_orb, &
966 : 1.0_dp, admm_env%B, mo_coeff, 0.0_dp, &
967 402 : admm_env%work_orb_nmo(ispin))
968 : CALL parallel_gemm('T', 'N', nmo, nmo, nao_orb, &
969 : 1.0_dp, mo_coeff, admm_env%work_orb_nmo(ispin), 0.0_dp, &
970 402 : admm_env%lambda(ispin))
971 402 : CALL cp_fm_to_fm(admm_env%lambda(ispin), admm_env%work_nmo_nmo1(ispin))
972 :
973 402 : CALL cp_fm_cholesky_decompose(admm_env%work_nmo_nmo1(ispin))
974 402 : CALL cp_fm_cholesky_invert(admm_env%work_nmo_nmo1(ispin))
975 : !! Symmetrize the guy
976 402 : CALL cp_fm_uplo_to_full(admm_env%work_nmo_nmo1(ispin), admm_env%lambda_inv(ispin))
977 402 : CALL cp_fm_to_fm(admm_env%work_nmo_nmo1(ispin), admm_env%lambda_inv(ispin))
978 :
979 : !! ** C_hat = AC
980 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
981 : 1.0_dp, admm_env%A, mo_coeff, 0.0_dp, &
982 402 : admm_env%C_hat(ispin))
983 670 : CALL cp_fm_to_fm(admm_env%C_hat(ispin), mo_coeff_aux_fit)
984 :
985 : END DO
986 :
987 268 : CALL timestop(handle)
988 :
989 268 : END SUBROUTINE purify_mo_cholesky
990 :
991 : ! **************************************************************************************************
992 : !> \brief Calculates the MO coefficients for the auxiliary fitting basis set
993 : !> by minimizing int (psi_i - psi_aux_i)^2 using Lagrangian Multipliers
994 : !>
995 : !> \param admm_env The ADMM env
996 : !> \param mos the MO's of the orbital basis set
997 : !> \param mos_aux_fit the MO's of the auxiliary fitting basis set
998 : !> \par History
999 : !> 05.2008 created [Manuel Guidon]
1000 : !> \author Manuel Guidon
1001 : ! **************************************************************************************************
1002 1628 : SUBROUTINE purify_mo_diag(admm_env, mos, mos_aux_fit)
1003 :
1004 : TYPE(admm_type), POINTER :: admm_env
1005 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mos_aux_fit
1006 :
1007 : CHARACTER(LEN=*), PARAMETER :: routineN = 'purify_mo_diag'
1008 :
1009 : INTEGER :: handle, i, ispin, nao_aux_fit, nao_orb, &
1010 : nmo, nspins
1011 1628 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig_work
1012 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
1013 :
1014 1628 : CALL timeset(routineN, handle)
1015 :
1016 1628 : nao_aux_fit = admm_env%nao_aux_fit
1017 1628 : nao_orb = admm_env%nao_orb
1018 1628 : nspins = SIZE(mos)
1019 :
1020 : ! *** Calculate the mo_coeffs for the fitting basis
1021 3628 : DO ispin = 1, nspins
1022 2000 : nmo = admm_env%nmo(ispin)
1023 2000 : IF (nmo == 0) CYCLE
1024 : !! Lambda = C^(T)*B*C
1025 2000 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
1026 2000 : CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
1027 : CALL parallel_gemm('N', 'N', nao_orb, nmo, nao_orb, &
1028 : 1.0_dp, admm_env%B, mo_coeff, 0.0_dp, &
1029 2000 : admm_env%work_orb_nmo(ispin))
1030 : CALL parallel_gemm('T', 'N', nmo, nmo, nao_orb, &
1031 : 1.0_dp, mo_coeff, admm_env%work_orb_nmo(ispin), 0.0_dp, &
1032 2000 : admm_env%lambda(ispin))
1033 2000 : CALL cp_fm_to_fm(admm_env%lambda(ispin), admm_env%work_nmo_nmo1(ispin))
1034 :
1035 : CALL cp_fm_syevd(admm_env%work_nmo_nmo1(ispin), admm_env%R(ispin), &
1036 2000 : admm_env%eigvals_lambda(ispin)%eigvals%data)
1037 6000 : ALLOCATE (eig_work(nmo))
1038 10034 : DO i = 1, nmo
1039 10034 : eig_work(i) = 1.0_dp/SQRT(admm_env%eigvals_lambda(ispin)%eigvals%data(i))
1040 : END DO
1041 2000 : CALL cp_fm_to_fm(admm_env%R(ispin), admm_env%work_nmo_nmo1(ispin))
1042 2000 : CALL cp_fm_column_scale(admm_env%work_nmo_nmo1(ispin), eig_work)
1043 : CALL parallel_gemm('N', 'T', nmo, nmo, nmo, &
1044 : 1.0_dp, admm_env%work_nmo_nmo1(ispin), admm_env%R(ispin), 0.0_dp, &
1045 2000 : admm_env%lambda_inv_sqrt(ispin))
1046 : CALL parallel_gemm('N', 'N', nao_orb, nmo, nmo, &
1047 : 1.0_dp, mo_coeff, admm_env%lambda_inv_sqrt(ispin), 0.0_dp, &
1048 2000 : admm_env%work_orb_nmo(ispin))
1049 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
1050 : 1.0_dp, admm_env%A, admm_env%work_orb_nmo(ispin), 0.0_dp, &
1051 2000 : mo_coeff_aux_fit)
1052 :
1053 2000 : CALL cp_fm_to_fm(mo_coeff_aux_fit, admm_env%C_hat(ispin))
1054 2000 : CALL cp_fm_set_all(admm_env%lambda_inv(ispin), 0.0_dp, 1.0_dp)
1055 3628 : DEALLOCATE (eig_work)
1056 : END DO
1057 :
1058 1628 : CALL timestop(handle)
1059 :
1060 1628 : END SUBROUTINE purify_mo_diag
1061 :
1062 : ! **************************************************************************************************
1063 : !> \brief ...
1064 : !> \param admm_env ...
1065 : !> \param mos ...
1066 : !> \param mos_aux_fit ...
1067 : ! **************************************************************************************************
1068 11268 : SUBROUTINE purify_mo_none(admm_env, mos, mos_aux_fit)
1069 : TYPE(admm_type), POINTER :: admm_env
1070 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos, mos_aux_fit
1071 :
1072 : CHARACTER(LEN=*), PARAMETER :: routineN = 'purify_mo_none'
1073 :
1074 : INTEGER :: handle, ispin, nao_aux_fit, nao_orb, &
1075 : nmo, nmo_mos, nspins
1076 11268 : REAL(KIND=dp), DIMENSION(:), POINTER :: occ_num, occ_num_aux
1077 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
1078 :
1079 11268 : CALL timeset(routineN, handle)
1080 :
1081 11268 : nao_aux_fit = admm_env%nao_aux_fit
1082 11268 : nao_orb = admm_env%nao_orb
1083 11268 : nspins = SIZE(mos)
1084 :
1085 24464 : DO ispin = 1, nspins
1086 13196 : nmo = admm_env%nmo(ispin)
1087 13196 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_num, nmo=nmo_mos)
1088 : CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit, &
1089 13196 : occupation_numbers=occ_num_aux)
1090 :
1091 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
1092 : 1.0_dp, admm_env%A, mo_coeff, 0.0_dp, &
1093 13196 : mo_coeff_aux_fit)
1094 13196 : CALL cp_fm_to_fm(mo_coeff_aux_fit, admm_env%C_hat(ispin))
1095 :
1096 144352 : occ_num_aux(1:nmo) = occ_num(1:nmo)
1097 : ! XXXX should only be done first time XXXX
1098 13196 : CALL cp_fm_set_all(admm_env%lambda(ispin), 0.0_dp, 1.0_dp)
1099 13196 : CALL cp_fm_set_all(admm_env%lambda_inv(ispin), 0.0_dp, 1.0_dp)
1100 37660 : CALL cp_fm_set_all(admm_env%lambda_inv_sqrt(ispin), 0.0_dp, 1.0_dp)
1101 : END DO
1102 :
1103 11268 : CALL timestop(handle)
1104 :
1105 11268 : END SUBROUTINE purify_mo_none
1106 :
1107 : ! **************************************************************************************************
1108 : !> \brief ...
1109 : !> \param admm_env ...
1110 : !> \param mo_set ...
1111 : !> \param density_matrix ...
1112 : !> \param ispin ...
1113 : !> \param blocked ...
1114 : ! **************************************************************************************************
1115 484 : SUBROUTINE purify_dm_cauchy(admm_env, mo_set, density_matrix, ispin, blocked)
1116 :
1117 : TYPE(admm_type), POINTER :: admm_env
1118 : TYPE(mo_set_type), INTENT(IN) :: mo_set
1119 : TYPE(dbcsr_type), POINTER :: density_matrix
1120 : INTEGER :: ispin
1121 : LOGICAL, INTENT(IN) :: blocked
1122 :
1123 : CHARACTER(len=*), PARAMETER :: routineN = 'purify_dm_cauchy'
1124 :
1125 : INTEGER :: handle, i, nao_aux_fit, nao_orb, nmo, &
1126 : nspins
1127 : REAL(KIND=dp) :: pole
1128 : TYPE(cp_fm_type), POINTER :: mo_coeff_aux_fit
1129 :
1130 484 : CALL timeset(routineN, handle)
1131 :
1132 484 : nao_aux_fit = admm_env%nao_aux_fit
1133 484 : nao_orb = admm_env%nao_orb
1134 484 : nmo = admm_env%nmo(ispin)
1135 :
1136 484 : nspins = SIZE(admm_env%P_to_be_purified)
1137 :
1138 484 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff_aux_fit)
1139 :
1140 : !! * For the time beeing, get the P to be purified from the mo_coeffs
1141 : !! * This needs to be replaced with the a block modified P
1142 :
1143 484 : IF (.NOT. blocked) THEN
1144 : CALL parallel_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nmo, &
1145 : 1.0_dp, mo_coeff_aux_fit, mo_coeff_aux_fit, 0.0_dp, &
1146 250 : admm_env%P_to_be_purified(ispin))
1147 : END IF
1148 :
1149 484 : CALL cp_fm_to_fm(admm_env%S, admm_env%work_aux_aux)
1150 484 : CALL cp_fm_to_fm(admm_env%P_to_be_purified(ispin), admm_env%work_aux_aux2)
1151 :
1152 484 : CALL cp_fm_cholesky_decompose(admm_env%work_aux_aux)
1153 :
1154 484 : CALL cp_fm_cholesky_reduce(admm_env%work_aux_aux2, admm_env%work_aux_aux, itype=3)
1155 :
1156 : CALL cp_fm_syevd(admm_env%work_aux_aux2, admm_env%R_purify(ispin), &
1157 484 : admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data)
1158 :
1159 : CALL cp_fm_cholesky_restore(admm_env%R_purify(ispin), nao_aux_fit, admm_env%work_aux_aux, &
1160 484 : admm_env%work_aux_aux3, op="MULTIPLY", pos="LEFT", transa="T")
1161 :
1162 484 : CALL cp_fm_to_fm(admm_env%work_aux_aux3, admm_env%R_purify(ispin))
1163 :
1164 : ! *** Construct Matrix M for Hadamard Product
1165 484 : CALL cp_fm_set_all(admm_env%M_purify(ispin), 0.0_dp)
1166 : pole = 0.0_dp
1167 3140 : DO i = 1, nao_aux_fit
1168 2656 : pole = Heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i) - 0.5_dp)
1169 3140 : CALL cp_fm_set_element(admm_env%M_purify(ispin), i, i, pole)
1170 : END DO
1171 484 : CALL cp_fm_uplo_to_full(admm_env%M_purify(ispin), admm_env%work_aux_aux)
1172 :
1173 484 : CALL copy_dbcsr_to_fm(density_matrix, admm_env%work_aux_aux3)
1174 484 : CALL cp_fm_uplo_to_full(admm_env%work_aux_aux3, admm_env%work_aux_aux)
1175 :
1176 : ! ** S^(-1)*R
1177 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit, &
1178 : 1.0_dp, admm_env%S_inv, admm_env%R_purify(ispin), 0.0_dp, &
1179 484 : admm_env%work_aux_aux)
1180 : ! ** S^(-1)*R*M
1181 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit, &
1182 : 1.0_dp, admm_env%work_aux_aux, admm_env%M_purify(ispin), 0.0_dp, &
1183 484 : admm_env%work_aux_aux2)
1184 : ! ** S^(-1)*R*M*R^T*S^(-1)
1185 : CALL parallel_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nao_aux_fit, &
1186 : 1.0_dp, admm_env%work_aux_aux2, admm_env%work_aux_aux, 0.0_dp, &
1187 484 : admm_env%work_aux_aux3)
1188 :
1189 484 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux3, density_matrix, keep_sparsity=.TRUE.)
1190 :
1191 484 : IF (nspins == 1) THEN
1192 96 : CALL dbcsr_scale(density_matrix, 2.0_dp)
1193 : END IF
1194 :
1195 484 : CALL timestop(handle)
1196 :
1197 484 : END SUBROUTINE purify_dm_cauchy
1198 :
1199 : ! **************************************************************************************************
1200 : !> \brief ...
1201 : !> \param qs_env ...
1202 : ! **************************************************************************************************
1203 290 : SUBROUTINE merge_ks_matrix_cauchy(qs_env)
1204 : TYPE(qs_environment_type), POINTER :: qs_env
1205 :
1206 : CHARACTER(LEN=*), PARAMETER :: routineN = 'merge_ks_matrix_cauchy'
1207 :
1208 : INTEGER :: handle, i, iatom, ispin, j, jatom, &
1209 : nao_aux_fit, nao_orb, nmo
1210 : REAL(dp) :: eig_diff, pole, tmp
1211 290 : REAL(dp), DIMENSION(:, :), POINTER :: sparse_block
1212 : TYPE(admm_type), POINTER :: admm_env
1213 : TYPE(cp_fm_type), POINTER :: mo_coeff
1214 : TYPE(dbcsr_iterator_type) :: iter
1215 290 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit
1216 : TYPE(dbcsr_type), POINTER :: matrix_k_tilde
1217 : TYPE(dft_control_type), POINTER :: dft_control
1218 290 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1219 :
1220 290 : CALL timeset(routineN, handle)
1221 290 : NULLIFY (admm_env, dft_control, matrix_ks, matrix_ks_aux_fit, mos, mo_coeff)
1222 :
1223 : CALL get_qs_env(qs_env, &
1224 : admm_env=admm_env, &
1225 : dft_control=dft_control, &
1226 : matrix_ks=matrix_ks, &
1227 290 : mos=mos)
1228 290 : CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit)
1229 :
1230 774 : DO ispin = 1, dft_control%nspins
1231 484 : nao_aux_fit = admm_env%nao_aux_fit
1232 484 : nao_orb = admm_env%nao_orb
1233 484 : nmo = admm_env%nmo(ispin)
1234 484 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1235 :
1236 484 : IF (.NOT. admm_env%block_dm) THEN
1237 : !** Get P from mo_coeffs, otherwise we have troubles with occupation numbers ...
1238 : CALL parallel_gemm('N', 'T', nao_orb, nao_orb, nmo, &
1239 : 1.0_dp, mo_coeff, mo_coeff, 0.0_dp, &
1240 250 : admm_env%work_orb_orb)
1241 :
1242 : !! A*P
1243 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, &
1244 : 1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
1245 250 : admm_env%work_aux_orb2)
1246 : !! A*P*A^T
1247 : CALL parallel_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nao_orb, &
1248 : 1.0_dp, admm_env%work_aux_orb2, admm_env%A, 0.0_dp, &
1249 250 : admm_env%P_to_be_purified(ispin))
1250 :
1251 : END IF
1252 :
1253 484 : CALL cp_fm_to_fm(admm_env%S, admm_env%work_aux_aux)
1254 484 : CALL cp_fm_to_fm(admm_env%P_to_be_purified(ispin), admm_env%work_aux_aux2)
1255 :
1256 484 : CALL cp_fm_cholesky_decompose(admm_env%work_aux_aux)
1257 :
1258 484 : CALL cp_fm_cholesky_reduce(admm_env%work_aux_aux2, admm_env%work_aux_aux, itype=3)
1259 :
1260 : CALL cp_fm_syevd(admm_env%work_aux_aux2, admm_env%R_purify(ispin), &
1261 484 : admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data)
1262 :
1263 : CALL cp_fm_cholesky_restore(admm_env%R_purify(ispin), nao_aux_fit, admm_env%work_aux_aux, &
1264 484 : admm_env%work_aux_aux3, op="MULTIPLY", pos="LEFT", transa="T")
1265 :
1266 484 : CALL cp_fm_to_fm(admm_env%work_aux_aux3, admm_env%R_purify(ispin))
1267 :
1268 : ! *** Construct Matrix M for Hadamard Product
1269 484 : pole = 0.0_dp
1270 3140 : DO i = 1, nao_aux_fit
1271 14156 : DO j = i, nao_aux_fit
1272 : eig_diff = (admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i) - &
1273 11016 : admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(j))
1274 : ! *** two eigenvalues could be the degenerated. In that case use 2nd order formula for the poles
1275 13672 : IF (ABS(eig_diff) == 0.0_dp) THEN
1276 2754 : pole = delta(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i) - 0.5_dp)
1277 2754 : CALL cp_fm_set_element(admm_env%M_purify(ispin), i, j, pole)
1278 : ELSE
1279 : pole = 1.0_dp/(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i) - &
1280 8262 : admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(j))
1281 8262 : tmp = Heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i) - 0.5_dp)
1282 8262 : tmp = tmp - Heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(j) - 0.5_dp)
1283 8262 : pole = tmp*pole
1284 8262 : CALL cp_fm_set_element(admm_env%M_purify(ispin), i, j, pole)
1285 : END IF
1286 : END DO
1287 : END DO
1288 484 : CALL cp_fm_uplo_to_full(admm_env%M_purify(ispin), admm_env%work_aux_aux)
1289 :
1290 484 : CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
1291 484 : CALL cp_fm_uplo_to_full(admm_env%K(ispin), admm_env%work_aux_aux)
1292 :
1293 : !! S^(-1)*R
1294 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit, &
1295 : 1.0_dp, admm_env%S_inv, admm_env%R_purify(ispin), 0.0_dp, &
1296 484 : admm_env%work_aux_aux)
1297 : !! K*S^(-1)*R
1298 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit, &
1299 : 1.0_dp, admm_env%K(ispin), admm_env%work_aux_aux, 0.0_dp, &
1300 484 : admm_env%work_aux_aux2)
1301 : !! R^T*S^(-1)*K*S^(-1)*R
1302 : CALL parallel_gemm('T', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit, &
1303 : 1.0_dp, admm_env%work_aux_aux, admm_env%work_aux_aux2, 0.0_dp, &
1304 484 : admm_env%work_aux_aux3)
1305 : !! R^T*S^(-1)*K*S^(-1)*R x M
1306 : CALL cp_fm_schur_product(admm_env%work_aux_aux3, admm_env%M_purify(ispin), &
1307 484 : admm_env%work_aux_aux)
1308 :
1309 : !! R^T*A
1310 : CALL parallel_gemm('T', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
1311 : 1.0_dp, admm_env%R_purify(ispin), admm_env%A, 0.0_dp, &
1312 484 : admm_env%work_aux_orb)
1313 :
1314 : !! (R^T*S^(-1)*K*S^(-1)*R x M) * R^T*A
1315 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
1316 : 1.0_dp, admm_env%work_aux_aux, admm_env%work_aux_orb, 0.0_dp, &
1317 484 : admm_env%work_aux_orb2)
1318 : !! A^T*R*(R^T*S^(-1)*K*S^(-1)*R x M) * R^T*A
1319 : CALL parallel_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, &
1320 : 1.0_dp, admm_env%work_aux_orb, admm_env%work_aux_orb2, 0.0_dp, &
1321 484 : admm_env%work_orb_orb)
1322 :
1323 : NULLIFY (matrix_k_tilde)
1324 484 : ALLOCATE (matrix_k_tilde)
1325 : CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
1326 : name='MATRIX K_tilde', &
1327 484 : matrix_type=dbcsr_type_symmetric)
1328 :
1329 484 : CALL cp_fm_to_fm(admm_env%work_orb_orb, admm_env%ks_to_be_merged(ispin))
1330 :
1331 484 : CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
1332 484 : CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
1333 484 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE.)
1334 :
1335 484 : IF (admm_env%block_dm) THEN
1336 : ! ** now loop through the list and nullify blocks
1337 234 : CALL dbcsr_iterator_start(iter, matrix_k_tilde)
1338 851 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1339 617 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block)
1340 851 : IF (admm_env%block_map(iatom, jatom) == 0) THEN
1341 1206 : sparse_block = 0.0_dp
1342 : END IF
1343 : END DO
1344 234 : CALL dbcsr_iterator_stop(iter)
1345 : END IF
1346 :
1347 484 : CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
1348 :
1349 774 : CALL dbcsr_deallocate_matrix(matrix_k_tilde)
1350 :
1351 : END DO !spin-loop
1352 :
1353 290 : CALL timestop(handle)
1354 :
1355 290 : END SUBROUTINE merge_ks_matrix_cauchy
1356 :
1357 : ! **************************************************************************************************
1358 : !> \brief ...
1359 : !> \param qs_env ...
1360 : ! **************************************************************************************************
1361 154 : SUBROUTINE merge_ks_matrix_cauchy_subspace(qs_env)
1362 : TYPE(qs_environment_type), POINTER :: qs_env
1363 :
1364 : CHARACTER(LEN=*), PARAMETER :: routineN = 'merge_ks_matrix_cauchy_subspace'
1365 :
1366 : INTEGER :: handle, ispin, nao_aux_fit, nao_orb, nmo
1367 : TYPE(admm_type), POINTER :: admm_env
1368 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
1369 154 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit
1370 : TYPE(dbcsr_type), POINTER :: matrix_k_tilde
1371 : TYPE(dft_control_type), POINTER :: dft_control
1372 154 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
1373 :
1374 154 : CALL timeset(routineN, handle)
1375 154 : NULLIFY (admm_env, dft_control, matrix_ks, matrix_ks_aux_fit, mos, mos_aux_fit, &
1376 154 : mo_coeff, mo_coeff_aux_fit)
1377 :
1378 : CALL get_qs_env(qs_env, &
1379 : admm_env=admm_env, &
1380 : dft_control=dft_control, &
1381 : matrix_ks=matrix_ks, &
1382 154 : mos=mos)
1383 154 : CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit, mos_aux_fit=mos_aux_fit)
1384 :
1385 366 : DO ispin = 1, dft_control%nspins
1386 212 : nao_aux_fit = admm_env%nao_aux_fit
1387 212 : nao_orb = admm_env%nao_orb
1388 212 : nmo = admm_env%nmo(ispin)
1389 212 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1390 212 : CALL get_mo_set(mo_set=mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
1391 :
1392 : !! Calculate Lambda^{-2}
1393 212 : CALL cp_fm_to_fm(admm_env%lambda(ispin), admm_env%work_nmo_nmo1(ispin))
1394 212 : CALL cp_fm_cholesky_decompose(admm_env%work_nmo_nmo1(ispin))
1395 212 : CALL cp_fm_cholesky_invert(admm_env%work_nmo_nmo1(ispin))
1396 : !! Symmetrize the guy
1397 212 : CALL cp_fm_uplo_to_full(admm_env%work_nmo_nmo1(ispin), admm_env%lambda_inv2(ispin))
1398 : !! Take square
1399 : CALL parallel_gemm('N', 'T', nmo, nmo, nmo, &
1400 : 1.0_dp, admm_env%work_nmo_nmo1(ispin), admm_env%work_nmo_nmo1(ispin), 0.0_dp, &
1401 212 : admm_env%lambda_inv2(ispin))
1402 :
1403 : !! ** C_hat = AC
1404 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
1405 : 1.0_dp, admm_env%A, mo_coeff, 0.0_dp, &
1406 212 : admm_env%C_hat(ispin))
1407 :
1408 : !! calc P_tilde from C_hat
1409 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nmo, &
1410 : 1.0_dp, admm_env%C_hat(ispin), admm_env%lambda_inv(ispin), 0.0_dp, &
1411 212 : admm_env%work_aux_nmo(ispin))
1412 :
1413 : CALL parallel_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nmo, &
1414 : 1.0_dp, admm_env%C_hat(ispin), admm_env%work_aux_nmo(ispin), 0.0_dp, &
1415 212 : admm_env%P_tilde(ispin))
1416 :
1417 : !! ** C_hat*Lambda^{-2}
1418 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nmo, &
1419 : 1.0_dp, admm_env%C_hat(ispin), admm_env%lambda_inv2(ispin), 0.0_dp, &
1420 212 : admm_env%work_aux_nmo(ispin))
1421 :
1422 : !! ** C_hat*Lambda^{-2}*C_hat^T
1423 : CALL parallel_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nmo, &
1424 : 1.0_dp, admm_env%work_aux_nmo(ispin), admm_env%C_hat(ispin), 0.0_dp, &
1425 212 : admm_env%work_aux_aux)
1426 :
1427 : !! ** S*C_hat*Lambda^{-2}*C_hat^T
1428 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit, &
1429 : 1.0_dp, admm_env%S, admm_env%work_aux_aux, 0.0_dp, &
1430 212 : admm_env%work_aux_aux2)
1431 :
1432 212 : CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
1433 212 : CALL cp_fm_uplo_to_full(admm_env%K(ispin), admm_env%work_aux_aux)
1434 :
1435 : !! ** S*C_hat*Lambda^{-2}*C_hat^T*H_tilde
1436 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit, &
1437 : 1.0_dp, admm_env%work_aux_aux2, admm_env%K(ispin), 0.0_dp, &
1438 212 : admm_env%work_aux_aux)
1439 :
1440 : !! ** P_tilde*S
1441 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit, &
1442 : 1.0_dp, admm_env%P_tilde(ispin), admm_env%S, 0.0_dp, &
1443 212 : admm_env%work_aux_aux2)
1444 :
1445 : !! ** -S*C_hat*Lambda^{-2}*C_hat^T*H_tilde*P_tilde*S
1446 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit, &
1447 : -1.0_dp, admm_env%work_aux_aux, admm_env%work_aux_aux2, 0.0_dp, &
1448 212 : admm_env%work_aux_aux3)
1449 :
1450 : !! ** -S*C_hat*Lambda^{-2}*C_hat^T*H_tilde*P_tilde*S+S*C_hat*Lambda^{-2}*C_hat^T*H_tilde
1451 212 : CALL cp_fm_scale_and_add(1.0_dp, admm_env%work_aux_aux3, 1.0_dp, admm_env%work_aux_aux)
1452 :
1453 : !! first_part*A
1454 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
1455 : 1.0_dp, admm_env%work_aux_aux3, admm_env%A, 0.0_dp, &
1456 212 : admm_env%work_aux_orb)
1457 :
1458 : !! + first_part^T*A
1459 : CALL parallel_gemm('T', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
1460 : 1.0_dp, admm_env%work_aux_aux3, admm_env%A, 1.0_dp, &
1461 212 : admm_env%work_aux_orb)
1462 :
1463 : !! A^T*(first+seccond)=H
1464 : CALL parallel_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, &
1465 : 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1466 212 : admm_env%work_orb_orb)
1467 :
1468 : NULLIFY (matrix_k_tilde)
1469 212 : ALLOCATE (matrix_k_tilde)
1470 : CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
1471 : name='MATRIX K_tilde', &
1472 212 : matrix_type=dbcsr_type_symmetric)
1473 :
1474 212 : CALL cp_fm_to_fm(admm_env%work_orb_orb, admm_env%ks_to_be_merged(ispin))
1475 :
1476 212 : CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
1477 212 : CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
1478 212 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE.)
1479 :
1480 : CALL parallel_gemm('N', 'N', nao_orb, nmo, nao_orb, &
1481 : 1.0_dp, admm_env%work_orb_orb, mo_coeff, 0.0_dp, &
1482 212 : admm_env%mo_derivs_tmp(ispin))
1483 :
1484 212 : CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
1485 :
1486 366 : CALL dbcsr_deallocate_matrix(matrix_k_tilde)
1487 :
1488 : END DO !spin loop
1489 154 : CALL timestop(handle)
1490 :
1491 154 : END SUBROUTINE merge_ks_matrix_cauchy_subspace
1492 :
1493 : ! **************************************************************************************************
1494 : !> \brief Calculates the product Kohn-Sham-Matrix x mo_coeff for the auxiliary
1495 : !> basis set and transforms it into the orbital basis. This is needed
1496 : !> in order to use OT
1497 : !>
1498 : !> \param ispin which spin to transform
1499 : !> \param admm_env The ADMM env
1500 : !> \param mo_set ...
1501 : !> \param mo_coeff the MO coefficients from the orbital basis set
1502 : !> \param mo_coeff_aux_fit the MO coefficients from the auxiliary fitting basis set
1503 : !> \param mo_derivs KS x mo_coeff from the orbital basis set to which we add the
1504 : !> auxiliary basis set part
1505 : !> \param mo_derivs_aux_fit ...
1506 : !> \param matrix_ks_aux_fit the Kohn-Sham matrix from the auxiliary fitting basis set
1507 : !> \par History
1508 : !> 05.2008 created [Manuel Guidon]
1509 : !> \author Manuel Guidon
1510 : ! **************************************************************************************************
1511 3498 : SUBROUTINE merge_mo_derivs_diag(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, mo_derivs, &
1512 1166 : mo_derivs_aux_fit, matrix_ks_aux_fit)
1513 : INTEGER, INTENT(IN) :: ispin
1514 : TYPE(admm_type), POINTER :: admm_env
1515 : TYPE(mo_set_type), INTENT(IN) :: mo_set
1516 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff, mo_coeff_aux_fit
1517 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_derivs, mo_derivs_aux_fit
1518 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit
1519 :
1520 : CHARACTER(LEN=*), PARAMETER :: routineN = 'merge_mo_derivs_diag'
1521 :
1522 : INTEGER :: handle, i, j, nao_aux_fit, nao_orb, nmo
1523 : REAL(dp) :: eig_diff, pole, tmp32, tmp52, tmp72, &
1524 : tmp92
1525 1166 : REAL(dp), DIMENSION(:), POINTER :: occupation_numbers, scaling_factor
1526 :
1527 1166 : CALL timeset(routineN, handle)
1528 :
1529 1166 : nao_aux_fit = admm_env%nao_aux_fit
1530 1166 : nao_orb = admm_env%nao_orb
1531 1166 : nmo = admm_env%nmo(ispin)
1532 :
1533 1166 : CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
1534 1166 : CALL cp_fm_uplo_to_full(admm_env%K(ispin), admm_env%work_aux_aux)
1535 :
1536 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_aux_fit, &
1537 : 1.0_dp, admm_env%K(ispin), mo_coeff_aux_fit, 0.0_dp, &
1538 1166 : admm_env%H(ispin))
1539 :
1540 1166 : CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
1541 3498 : ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
1542 10734 : scaling_factor = 2.0_dp*occupation_numbers
1543 :
1544 1166 : CALL cp_fm_column_scale(admm_env%H(ispin), scaling_factor)
1545 :
1546 1166 : CALL cp_fm_to_fm(admm_env%H(ispin), mo_derivs_aux_fit(ispin))
1547 :
1548 : ! *** Add first term
1549 : CALL parallel_gemm('N', 'T', nao_aux_fit, nmo, nmo, &
1550 : 1.0_dp, admm_env%H(ispin), admm_env%lambda_inv_sqrt(ispin), 0.0_dp, &
1551 1166 : admm_env%work_aux_nmo(ispin))
1552 : CALL parallel_gemm('T', 'N', nao_orb, nmo, nao_aux_fit, &
1553 : 1.0_dp, admm_env%A, admm_env%work_aux_nmo(ispin), 0.0_dp, &
1554 1166 : admm_env%mo_derivs_tmp(ispin))
1555 :
1556 : ! *** Construct Matrix M for Hadamard Product
1557 : pole = 0.0_dp
1558 5950 : DO i = 1, nmo
1559 21538 : DO j = i, nmo
1560 : eig_diff = (admm_env%eigvals_lambda(ispin)%eigvals%data(i) - &
1561 15588 : admm_env%eigvals_lambda(ispin)%eigvals%data(j))
1562 : ! *** two eigenvalues could be the degenerated. In that case use 2nd order formula for the poles
1563 20372 : IF (ABS(eig_diff) < 0.0001_dp) THEN
1564 6558 : tmp32 = 1.0_dp/SQRT(admm_env%eigvals_lambda(ispin)%eigvals%data(j))**3
1565 6558 : tmp52 = tmp32/admm_env%eigvals_lambda(ispin)%eigvals%data(j)*eig_diff
1566 6558 : tmp72 = tmp52/admm_env%eigvals_lambda(ispin)%eigvals%data(j)*eig_diff
1567 6558 : tmp92 = tmp72/admm_env%eigvals_lambda(ispin)%eigvals%data(j)*eig_diff
1568 :
1569 6558 : pole = -0.5_dp*tmp32 + 3.0_dp/8.0_dp*tmp52 - 5.0_dp/16.0_dp*tmp72 + 35.0_dp/128.0_dp*tmp92
1570 6558 : CALL cp_fm_set_element(admm_env%M(ispin), i, j, pole)
1571 : ELSE
1572 9030 : pole = 1.0_dp/SQRT(admm_env%eigvals_lambda(ispin)%eigvals%data(i))
1573 9030 : pole = pole - 1.0_dp/SQRT(admm_env%eigvals_lambda(ispin)%eigvals%data(j))
1574 : pole = pole/(admm_env%eigvals_lambda(ispin)%eigvals%data(i) - &
1575 9030 : admm_env%eigvals_lambda(ispin)%eigvals%data(j))
1576 9030 : CALL cp_fm_set_element(admm_env%M(ispin), i, j, pole)
1577 : END IF
1578 : END DO
1579 : END DO
1580 1166 : CALL cp_fm_uplo_to_full(admm_env%M(ispin), admm_env%work_nmo_nmo1(ispin))
1581 :
1582 : ! *** 2nd term to be added to fm_H
1583 :
1584 : !! Part 1: B^(T)*C* R*[R^(T)*c^(T)*A^(T)*H_aux_fit*R x M]*R^(T)
1585 : !! Part 2: B*C*(R*[R^(T)*c^(T)*A^(T)*H_aux_fit*R x M]*R^(T))^(T)
1586 :
1587 : ! *** H'*R
1588 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nmo, &
1589 : 1.0_dp, admm_env%H(ispin), admm_env%R(ispin), 0.0_dp, &
1590 1166 : admm_env%work_aux_nmo(ispin))
1591 : ! *** A^(T)*H'*R
1592 : CALL parallel_gemm('T', 'N', nao_orb, nmo, nao_aux_fit, &
1593 : 1.0_dp, admm_env%A, admm_env%work_aux_nmo(ispin), 0.0_dp, &
1594 1166 : admm_env%work_orb_nmo(ispin))
1595 : ! *** c^(T)*A^(T)*H'*R
1596 : CALL parallel_gemm('T', 'N', nmo, nmo, nao_orb, &
1597 : 1.0_dp, mo_coeff, admm_env%work_orb_nmo(ispin), 0.0_dp, &
1598 1166 : admm_env%work_nmo_nmo1(ispin))
1599 : ! *** R^(T)*c^(T)*A^(T)*H'*R
1600 : CALL parallel_gemm('T', 'N', nmo, nmo, nmo, &
1601 : 1.0_dp, admm_env%R(ispin), admm_env%work_nmo_nmo1(ispin), 0.0_dp, &
1602 1166 : admm_env%work_nmo_nmo2(ispin))
1603 : ! *** R^(T)*c^(T)*A^(T)*H'*R x M
1604 : CALL cp_fm_schur_product(admm_env%work_nmo_nmo2(ispin), &
1605 1166 : admm_env%M(ispin), admm_env%work_nmo_nmo1(ispin))
1606 : ! *** R* (R^(T)*c^(T)*A^(T)*H'*R x M)
1607 : CALL parallel_gemm('N', 'N', nmo, nmo, nmo, &
1608 : 1.0_dp, admm_env%R(ispin), admm_env%work_nmo_nmo1(ispin), 0.0_dp, &
1609 1166 : admm_env%work_nmo_nmo2(ispin))
1610 :
1611 : ! *** R* (R^(T)*c^(T)*A^(T)*H'*R x M) *R^(T)
1612 : CALL parallel_gemm('N', 'T', nmo, nmo, nmo, &
1613 : 1.0_dp, admm_env%work_nmo_nmo2(ispin), admm_env%R(ispin), 0.0_dp, &
1614 1166 : admm_env%R_schur_R_t(ispin))
1615 :
1616 : ! *** B^(T)*c
1617 : CALL parallel_gemm('T', 'N', nao_orb, nmo, nao_orb, &
1618 : 1.0_dp, admm_env%B, mo_coeff, 0.0_dp, &
1619 1166 : admm_env%work_orb_nmo(ispin))
1620 :
1621 : ! *** Add first term to fm_H
1622 : ! *** B^(T)*c* R* (R^(T)*c^(T)*A^(T)*H'*R x M) *R^(T)
1623 : CALL parallel_gemm('N', 'N', nao_orb, nmo, nmo, &
1624 : 1.0_dp, admm_env%work_orb_nmo(ispin), admm_env%R_schur_R_t(ispin), 1.0_dp, &
1625 1166 : admm_env%mo_derivs_tmp(ispin))
1626 :
1627 : ! *** Add second term to fm_H
1628 : ! *** B*C *[ R* (R^(T)*c^(T)*A^(T)*H'*R x M) *R^(T)]^(T)
1629 : CALL parallel_gemm('N', 'T', nao_orb, nmo, nmo, &
1630 : 1.0_dp, admm_env%work_orb_nmo(ispin), admm_env%R_schur_R_t(ispin), 1.0_dp, &
1631 1166 : admm_env%mo_derivs_tmp(ispin))
1632 :
1633 5950 : DO i = 1, SIZE(scaling_factor)
1634 5950 : scaling_factor(i) = 1.0_dp/scaling_factor(i)
1635 : END DO
1636 :
1637 1166 : CALL cp_fm_column_scale(admm_env%mo_derivs_tmp(ispin), scaling_factor)
1638 :
1639 1166 : CALL cp_fm_scale_and_add(1.0_dp, mo_derivs(ispin), 1.0_dp, admm_env%mo_derivs_tmp(ispin))
1640 :
1641 1166 : DEALLOCATE (scaling_factor)
1642 :
1643 1166 : CALL timestop(handle)
1644 :
1645 1166 : END SUBROUTINE merge_mo_derivs_diag
1646 :
1647 : ! **************************************************************************************************
1648 : !> \brief ...
1649 : !> \param qs_env ...
1650 : ! **************************************************************************************************
1651 10962 : SUBROUTINE merge_ks_matrix_none(qs_env)
1652 : TYPE(qs_environment_type), POINTER :: qs_env
1653 :
1654 : CHARACTER(LEN=*), PARAMETER :: routineN = 'merge_ks_matrix_none'
1655 :
1656 : INTEGER :: handle, iatom, ispin, jatom, &
1657 : nao_aux_fit, nao_orb, nmo
1658 10962 : REAL(dp), DIMENSION(:, :), POINTER :: sparse_block
1659 : REAL(KIND=dp) :: ener_k(2), ener_x(2), ener_x1(2), &
1660 : gsi_square, trace_tmp, trace_tmp_two
1661 : TYPE(admm_type), POINTER :: admm_env
1662 : TYPE(dbcsr_iterator_type) :: iter
1663 10962 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
1664 10962 : matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_s, matrix_s_aux_fit, rho_ao, &
1665 10962 : rho_ao_aux
1666 : TYPE(dbcsr_type), POINTER :: matrix_k_tilde, &
1667 : matrix_ks_aux_fit_admms_tmp, &
1668 : matrix_TtsT
1669 : TYPE(dft_control_type), POINTER :: dft_control
1670 : TYPE(mp_para_env_type), POINTER :: para_env
1671 : TYPE(qs_energy_type), POINTER :: energy
1672 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
1673 :
1674 10962 : CALL timeset(routineN, handle)
1675 10962 : NULLIFY (admm_env, dft_control, matrix_ks, matrix_ks_aux_fit, matrix_ks_aux_fit_dft, &
1676 10962 : matrix_ks_aux_fit_hfx, matrix_s, matrix_s_aux_fit, rho_ao, rho_ao_aux, matrix_k_tilde, &
1677 10962 : matrix_TtsT, matrix_ks_aux_fit_admms_tmp, rho, rho_aux_fit, sparse_block, para_env, energy)
1678 :
1679 : CALL get_qs_env(qs_env, &
1680 : admm_env=admm_env, &
1681 : dft_control=dft_control, &
1682 : matrix_ks=matrix_ks, &
1683 : rho=rho, &
1684 : matrix_s=matrix_s, &
1685 : energy=energy, &
1686 10962 : para_env=para_env)
1687 : CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit, matrix_ks_aux_fit_dft=matrix_ks_aux_fit_dft, &
1688 : matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx, rho_aux_fit=rho_aux_fit, &
1689 10962 : matrix_s_aux_fit=matrix_s_aux_fit)
1690 :
1691 10962 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1692 : CALL qs_rho_get(rho_aux_fit, &
1693 10962 : rho_ao=rho_ao_aux)
1694 :
1695 23660 : DO ispin = 1, dft_control%nspins
1696 23660 : IF (admm_env%block_dm) THEN
1697 120 : CALL dbcsr_iterator_start(iter, matrix_ks_aux_fit(ispin)%matrix)
1698 832 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1699 712 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block)
1700 832 : IF (admm_env%block_map(iatom, jatom) == 0) THEN
1701 1890 : sparse_block = 0.0_dp
1702 : END IF
1703 : END DO
1704 120 : CALL dbcsr_iterator_stop(iter)
1705 120 : CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_ks_aux_fit(ispin)%matrix, 1.0_dp, 1.0_dp)
1706 :
1707 : ELSE
1708 :
1709 12578 : nao_aux_fit = admm_env%nao_aux_fit
1710 12578 : nao_orb = admm_env%nao_orb
1711 12578 : nmo = admm_env%nmo(ispin)
1712 :
1713 : ! ADMMS: different matrix for calculating A^(T)*K*A, see Eq. (37) Merlot
1714 12578 : IF (admm_env%do_admms) THEN
1715 : NULLIFY (matrix_ks_aux_fit_admms_tmp)
1716 392 : ALLOCATE (matrix_ks_aux_fit_admms_tmp)
1717 : CALL dbcsr_create(matrix_ks_aux_fit_admms_tmp, template=matrix_ks_aux_fit(ispin)%matrix, &
1718 392 : name='matrix_ks_aux_fit_admms_tmp', matrix_type='s')
1719 : ! matrix_ks_aux_fit_admms_tmp = k(d_Q)
1720 392 : CALL dbcsr_copy(matrix_ks_aux_fit_admms_tmp, matrix_ks_aux_fit_hfx(ispin)%matrix)
1721 :
1722 : ! matrix_ks_aux_fit_admms_tmp = k(d_Q) - gsi^2/3 x(d_Q)
1723 : CALL dbcsr_add(matrix_ks_aux_fit_admms_tmp, matrix_ks_aux_fit_dft(ispin)%matrix, &
1724 392 : 1.0_dp, -(admm_env%gsi(ispin))**(2.0_dp/3.0_dp))
1725 392 : CALL copy_dbcsr_to_fm(matrix_ks_aux_fit_admms_tmp, admm_env%K(ispin))
1726 392 : CALL dbcsr_deallocate_matrix(matrix_ks_aux_fit_admms_tmp)
1727 : ELSE
1728 12186 : CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
1729 : END IF
1730 :
1731 12578 : CALL cp_fm_uplo_to_full(admm_env%K(ispin), admm_env%work_aux_aux)
1732 :
1733 : !! K*A
1734 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
1735 : 1.0_dp, admm_env%K(ispin), admm_env%A, 0.0_dp, &
1736 12578 : admm_env%work_aux_orb)
1737 : !! A^T*K*A
1738 : CALL parallel_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, &
1739 : 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1740 12578 : admm_env%work_orb_orb)
1741 :
1742 : NULLIFY (matrix_k_tilde)
1743 12578 : ALLOCATE (matrix_k_tilde)
1744 : CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
1745 12578 : name='MATRIX K_tilde', matrix_type='S')
1746 12578 : CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
1747 12578 : CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
1748 12578 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE.)
1749 :
1750 : ! Scale matrix_K_tilde here. Then, the scaling has to be done for forces separately
1751 : ! Scale matrix_K_tilde by gsi for ADMMQ and ADMMS (Eqs. (27), (37) in Merlot, 2014)
1752 12578 : IF (admm_env%do_admmq .OR. admm_env%do_admms) THEN
1753 600 : CALL dbcsr_scale(matrix_k_tilde, admm_env%gsi(ispin))
1754 : END IF
1755 :
1756 : ! Scale matrix_K_tilde by gsi^2 for ADMMP (Eq. (35) in Merlot, 2014)
1757 12578 : IF (admm_env%do_admmp) THEN
1758 430 : gsi_square = (admm_env%gsi(ispin))*(admm_env%gsi(ispin))
1759 430 : CALL dbcsr_scale(matrix_k_tilde, gsi_square)
1760 : END IF
1761 :
1762 12578 : admm_env%lambda_merlot(ispin) = 0
1763 :
1764 : ! Calculate LAMBDA according to Merlot, 1. IF: ADMMQ, 2. IF: ADMMP, 3. IF: ADMMS,
1765 12578 : IF (admm_env%do_admmq) THEN
1766 208 : CALL dbcsr_dot(matrix_ks_aux_fit(ispin)%matrix, rho_ao_aux(ispin)%matrix, trace_tmp)
1767 :
1768 : ! Factor of 2 is missing compared to Eq. 28 in Merlot due to
1769 : ! Tr(ds) = N in the code \neq 2N in Merlot
1770 208 : admm_env%lambda_merlot(ispin) = trace_tmp/(admm_env%n_large_basis(ispin))
1771 :
1772 12370 : ELSE IF (admm_env%do_admmp) THEN
1773 430 : IF (dft_control%nspins == 2) THEN
1774 : CALL calc_spin_dep_aux_exch_ener(qs_env=qs_env, admm_env=admm_env, ener_k_ispin=ener_k(ispin), &
1775 : ener_x_ispin=ener_x(ispin), ener_x1_ispin=ener_x1(ispin), &
1776 52 : ispin=ispin)
1777 : admm_env%lambda_merlot(ispin) = 2.0_dp*(admm_env%gsi(ispin))**2* &
1778 : (ener_k(ispin) + ener_x(ispin) + ener_x1(ispin))/ &
1779 52 : (admm_env%n_large_basis(ispin))
1780 :
1781 : ELSE
1782 : admm_env%lambda_merlot(ispin) = 2.0_dp*(admm_env%gsi(ispin))**2* &
1783 : (energy%ex + energy%exc_aux_fit + energy%exc1_aux_fit) &
1784 378 : /(admm_env%n_large_basis(ispin))
1785 : END IF
1786 :
1787 11940 : ELSE IF (admm_env%do_admms) THEN
1788 392 : CALL dbcsr_dot(matrix_ks_aux_fit_hfx(ispin)%matrix, rho_ao_aux(ispin)%matrix, trace_tmp)
1789 392 : CALL dbcsr_dot(matrix_ks_aux_fit_dft(ispin)%matrix, rho_ao_aux(ispin)%matrix, trace_tmp_two)
1790 : ! For ADMMS open-shell case we need k and x (Merlot) separately since gsi(a)\=gsi(b)
1791 392 : IF (dft_control%nspins == 2) THEN
1792 : CALL calc_spin_dep_aux_exch_ener(qs_env=qs_env, admm_env=admm_env, ener_k_ispin=ener_k(ispin), &
1793 : ener_x_ispin=ener_x(ispin), ener_x1_ispin=ener_x1(ispin), &
1794 324 : ispin=ispin)
1795 : admm_env%lambda_merlot(ispin) = &
1796 : (trace_tmp + 2.0_dp/3.0_dp*((admm_env%gsi(ispin))**(2.0_dp/3.0_dp))* &
1797 : (ener_x(ispin) + ener_x1(ispin)) - ((admm_env%gsi(ispin))**(2.0_dp/3.0_dp))* &
1798 324 : trace_tmp_two)/(admm_env%n_large_basis(ispin))
1799 :
1800 : ELSE
1801 : admm_env%lambda_merlot(ispin) = (trace_tmp + (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)* &
1802 : (2.0_dp/3.0_dp*(energy%exc_aux_fit + energy%exc1_aux_fit) - &
1803 68 : trace_tmp_two))/(admm_env%n_large_basis(ispin))
1804 : END IF
1805 : END IF
1806 :
1807 : ! Calculate variational distribution to KS matrix according
1808 : ! to Eqs. (27), (35) and (37) in Merlot, 2014
1809 :
1810 12578 : IF (admm_env%do_admmp .OR. admm_env%do_admmq .OR. admm_env%do_admms) THEN
1811 :
1812 : !! T^T*s_aux*T in (27) Merlot (T=A), as calculating A^T*K*A few lines above
1813 1030 : CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix, admm_env%work_aux_aux4)
1814 1030 : CALL cp_fm_uplo_to_full(admm_env%work_aux_aux4, admm_env%work_aux_aux5)
1815 :
1816 : ! s_aux*T
1817 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
1818 : 1.0_dp, admm_env%work_aux_aux4, admm_env%A, 0.0_dp, &
1819 1030 : admm_env%work_aux_orb3)
1820 : ! T^T*s_aux*T
1821 : CALL parallel_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, &
1822 : 1.0_dp, admm_env%A, admm_env%work_aux_orb3, 0.0_dp, &
1823 1030 : admm_env%work_orb_orb3)
1824 :
1825 : NULLIFY (matrix_TtsT)
1826 1030 : ALLOCATE (matrix_TtsT)
1827 : CALL dbcsr_create(matrix_TtsT, template=matrix_ks(ispin)%matrix, &
1828 1030 : name='MATRIX TtsT', matrix_type='S')
1829 1030 : CALL dbcsr_copy(matrix_TtsT, matrix_ks(ispin)%matrix)
1830 1030 : CALL dbcsr_set(matrix_TtsT, 0.0_dp)
1831 1030 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb3, matrix_TtsT, keep_sparsity=.TRUE.)
1832 :
1833 : !Add -(gsi)*Lambda*TtsT and Lambda*S to the KS matrix according to Merlot2014
1834 :
1835 : CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_TtsT, 1.0_dp, &
1836 1030 : (-admm_env%lambda_merlot(ispin))*admm_env%gsi(ispin))
1837 :
1838 1030 : CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_s(1)%matrix, 1.0_dp, admm_env%lambda_merlot(ispin))
1839 :
1840 1030 : CALL dbcsr_deallocate_matrix(matrix_TtsT)
1841 :
1842 : END IF
1843 :
1844 12578 : CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
1845 :
1846 12578 : CALL dbcsr_deallocate_matrix(matrix_k_tilde)
1847 :
1848 : END IF
1849 : END DO !spin loop
1850 :
1851 : ! Scale energy for ADMMP and ADMMS
1852 10962 : IF (admm_env%do_admmp) THEN
1853 : ! ener_k = ener_k*(admm_env%gsi(1))*(admm_env%gsi(1))
1854 : ! ener_x = ener_x*(admm_env%gsi(1))*(admm_env%gsi(1))
1855 : ! PRINT *, 'energy%ex = ', energy%ex
1856 404 : IF (dft_control%nspins == 2) THEN
1857 26 : energy%exc_aux_fit = 0.0_dp
1858 26 : energy%exc1_aux_fit = 0.0_dp
1859 26 : energy%ex = 0.0_dp
1860 78 : DO ispin = 1, dft_control%nspins
1861 52 : energy%exc_aux_fit = energy%exc_aux_fit + (admm_env%gsi(ispin))**2.0_dp*ener_x(ispin)
1862 52 : energy%exc1_aux_fit = energy%exc1_aux_fit + (admm_env%gsi(ispin))**2.0_dp*ener_x1(ispin)
1863 78 : energy%ex = energy%ex + (admm_env%gsi(ispin))**2.0_dp*ener_k(ispin)
1864 : END DO
1865 : ELSE
1866 378 : energy%exc_aux_fit = (admm_env%gsi(1))**2.0_dp*energy%exc_aux_fit
1867 378 : energy%exc1_aux_fit = (admm_env%gsi(1))**2.0_dp*energy%exc1_aux_fit
1868 378 : energy%ex = (admm_env%gsi(1))**2.0_dp*energy%ex
1869 : END IF
1870 :
1871 10558 : ELSE IF (admm_env%do_admms) THEN
1872 230 : IF (dft_control%nspins == 2) THEN
1873 162 : energy%exc_aux_fit = 0.0_dp
1874 162 : energy%exc1_aux_fit = 0.0_dp
1875 486 : DO ispin = 1, dft_control%nspins
1876 324 : energy%exc_aux_fit = energy%exc_aux_fit + (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)*ener_x(ispin)
1877 486 : energy%exc1_aux_fit = energy%exc1_aux_fit + (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)*ener_x1(ispin)
1878 : END DO
1879 : ELSE
1880 68 : energy%exc_aux_fit = (admm_env%gsi(1))**(2.0_dp/3.0_dp)*energy%exc_aux_fit
1881 68 : energy%exc1_aux_fit = (admm_env%gsi(1))**(2.0_dp/3.0_dp)*energy%exc1_aux_fit
1882 : END IF
1883 : END IF
1884 :
1885 10962 : CALL timestop(handle)
1886 :
1887 10962 : END SUBROUTINE merge_ks_matrix_none
1888 :
1889 : ! **************************************************************************************************
1890 : !> \brief ...
1891 : !> \param qs_env ...
1892 : ! **************************************************************************************************
1893 242 : SUBROUTINE merge_ks_matrix_none_kp(qs_env)
1894 : TYPE(qs_environment_type), POINTER :: qs_env
1895 :
1896 : CHARACTER(LEN=*), PARAMETER :: routineN = 'merge_ks_matrix_none_kp'
1897 :
1898 : COMPLEX(dp) :: fac, fac2
1899 : INTEGER :: handle, i, igroup, ik, ikp, img, indx, &
1900 : ispin, kplocal, kpmax, nao_aux_fit, &
1901 : nao_orb, natom, nkp, nkp_groups, nspins
1902 : INTEGER, DIMENSION(2) :: kp_range
1903 242 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1904 242 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1905 : LOGICAL :: my_kpgrp, use_real_wfn
1906 : REAL(dp) :: ener_k(2), ener_x(2), ener_x1(2), tmp, &
1907 : trace_tmp, trace_tmp_two
1908 242 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1909 : TYPE(admm_type), POINTER :: admm_env
1910 242 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
1911 : TYPE(cp_cfm_type) :: cK, cwork_aux_orb, cwork_orb_orb
1912 : TYPE(cp_fm_struct_type), POINTER :: struct_aux_aux, struct_aux_orb, &
1913 : struct_orb_orb
1914 : TYPE(cp_fm_type) :: fmdummy, work_aux_aux, work_aux_aux2, &
1915 : work_aux_orb
1916 242 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fmwork
1917 242 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_ks
1918 242 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_k_tilde, matrix_ks_aux_fit, &
1919 242 : matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_kp, matrix_s, matrix_s_aux_fit, &
1920 242 : rho_ao_aux
1921 : TYPE(dbcsr_type) :: tmpmatrix_ks
1922 242 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: ksmatrix
1923 : TYPE(dft_control_type), POINTER :: dft_control
1924 : TYPE(kpoint_env_type), POINTER :: kp
1925 : TYPE(kpoint_type), POINTER :: kpoints
1926 : TYPE(mp_para_env_type), POINTER :: para_env
1927 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1928 242 : POINTER :: sab_aux_fit, sab_kp
1929 : TYPE(qs_energy_type), POINTER :: energy
1930 : TYPE(qs_rho_type), POINTER :: rho_aux_fit
1931 : TYPE(qs_scf_env_type), POINTER :: scf_env
1932 :
1933 242 : CALL timeset(routineN, handle)
1934 242 : NULLIFY (admm_env, rho_ao_aux, rho_aux_fit, &
1935 242 : matrix_s_aux_fit, energy, &
1936 242 : para_env, kpoints, sab_aux_fit, &
1937 242 : matrix_k_tilde, matrix_ks_kp, matrix_ks_aux_fit, scf_env, &
1938 242 : struct_orb_orb, struct_aux_orb, struct_aux_aux, kp, &
1939 242 : matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft)
1940 :
1941 : CALL get_qs_env(qs_env, &
1942 : admm_env=admm_env, &
1943 : dft_control=dft_control, &
1944 : matrix_ks_kp=matrix_ks_kp, &
1945 : matrix_s_kp=matrix_s, &
1946 : para_env=para_env, &
1947 : scf_env=scf_env, &
1948 : natom=natom, &
1949 : kpoints=kpoints, &
1950 242 : energy=energy)
1951 :
1952 : CALL get_admm_env(admm_env, &
1953 : matrix_ks_aux_fit_kp=matrix_ks_aux_fit, &
1954 : matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx, &
1955 : matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft, &
1956 : matrix_s_aux_fit_kp=matrix_s_aux_fit, &
1957 : sab_aux_fit=sab_aux_fit, &
1958 242 : rho_aux_fit=rho_aux_fit)
1959 242 : CALL qs_rho_get(rho_aux_fit, rho_ao_kp=rho_ao_aux)
1960 :
1961 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
1962 : nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_kp, &
1963 242 : cell_to_index=cell_to_index)
1964 :
1965 242 : nao_aux_fit = admm_env%nao_aux_fit
1966 242 : nao_orb = admm_env%nao_orb
1967 242 : nspins = dft_control%nspins
1968 :
1969 : !Case study on ADMMQ, ADMMS and ADMMP
1970 :
1971 : !ADMMQ: calculate lamda as in Merlot eq (28)
1972 242 : IF (admm_env%do_admmq) THEN
1973 30 : admm_env%lambda_merlot = 0.0_dp
1974 506 : DO img = 1, dft_control%nimages
1975 1002 : DO ispin = 1, nspins
1976 496 : CALL dbcsr_dot(matrix_ks_aux_fit(ispin, img)%matrix, rho_ao_aux(ispin, img)%matrix, trace_tmp)
1977 992 : admm_env%lambda_merlot(ispin) = admm_env%lambda_merlot(ispin) + trace_tmp/admm_env%n_large_basis(ispin)
1978 : END DO
1979 : END DO
1980 : END IF
1981 :
1982 : !ADMMP: calculate lamda as in Merlot eq (34)
1983 242 : IF (admm_env%do_admmp) THEN
1984 14 : IF (nspins == 1) THEN
1985 : admm_env%lambda_merlot(1) = 2.0_dp*(admm_env%gsi(1))**2* &
1986 : (energy%ex + energy%exc_aux_fit + energy%exc1_aux_fit) &
1987 14 : /(admm_env%n_large_basis(1))
1988 : ELSE
1989 0 : DO ispin = 1, nspins
1990 : CALL calc_spin_dep_aux_exch_ener(qs_env=qs_env, admm_env=admm_env, &
1991 : ener_k_ispin=ener_k(ispin), ener_x_ispin=ener_x(ispin), &
1992 0 : ener_x1_ispin=ener_x1(ispin), ispin=ispin)
1993 : admm_env%lambda_merlot(ispin) = 2.0_dp*(admm_env%gsi(ispin))**2* &
1994 : (ener_k(ispin) + ener_x(ispin) + ener_x1(ispin))/ &
1995 0 : (admm_env%n_large_basis(ispin))
1996 : END DO
1997 : END IF
1998 : END IF
1999 :
2000 : !ADMMS: calculate lambda as in Merlot eq (36)
2001 242 : IF (admm_env%do_admms) THEN
2002 70 : IF (nspins == 1) THEN
2003 56 : trace_tmp = 0.0_dp
2004 56 : trace_tmp_two = 0.0_dp
2005 2352 : DO img = 1, dft_control%nimages
2006 2296 : CALL dbcsr_dot(matrix_ks_aux_fit_hfx(1, img)%matrix, rho_ao_aux(1, img)%matrix, tmp)
2007 2296 : trace_tmp = trace_tmp + tmp
2008 2296 : CALL dbcsr_dot(matrix_ks_aux_fit_dft(1, img)%matrix, rho_ao_aux(1, img)%matrix, tmp)
2009 2352 : trace_tmp_two = trace_tmp_two + tmp
2010 : END DO
2011 : admm_env%lambda_merlot(1) = (trace_tmp + (admm_env%gsi(1))**(2.0_dp/3.0_dp)* &
2012 : (2.0_dp/3.0_dp*(energy%exc_aux_fit + energy%exc1_aux_fit) - &
2013 56 : trace_tmp_two))/(admm_env%n_large_basis(1))
2014 : ELSE
2015 :
2016 42 : DO ispin = 1, nspins
2017 28 : trace_tmp = 0.0_dp
2018 28 : trace_tmp_two = 0.0_dp
2019 488 : DO img = 1, dft_control%nimages
2020 460 : CALL dbcsr_dot(matrix_ks_aux_fit_hfx(ispin, img)%matrix, rho_ao_aux(ispin, img)%matrix, tmp)
2021 460 : trace_tmp = trace_tmp + tmp
2022 460 : CALL dbcsr_dot(matrix_ks_aux_fit_dft(ispin, img)%matrix, rho_ao_aux(ispin, img)%matrix, tmp)
2023 488 : trace_tmp_two = trace_tmp_two + tmp
2024 : END DO
2025 :
2026 : CALL calc_spin_dep_aux_exch_ener(qs_env=qs_env, admm_env=admm_env, &
2027 : ener_k_ispin=ener_k(ispin), ener_x_ispin=ener_x(ispin), &
2028 28 : ener_x1_ispin=ener_x1(ispin), ispin=ispin)
2029 :
2030 : admm_env%lambda_merlot(ispin) = &
2031 : (trace_tmp + 2.0_dp/3.0_dp*((admm_env%gsi(ispin))**(2.0_dp/3.0_dp))* &
2032 : (ener_x(ispin) + ener_x1(ispin)) - ((admm_env%gsi(ispin))**(2.0_dp/3.0_dp))* &
2033 42 : trace_tmp_two)/(admm_env%n_large_basis(ispin))
2034 : END DO
2035 : END IF
2036 :
2037 : !Here we buld the KS matrix: KS_hfx = gsi^2/3*KS_dft, the we then pass as the ususal KS_aux_fit
2038 70 : NULLIFY (matrix_ks_aux_fit)
2039 5562 : ALLOCATE (matrix_ks_aux_fit(nspins, dft_control%nimages))
2040 2596 : DO img = 1, dft_control%nimages
2041 5352 : DO ispin = 1, nspins
2042 2756 : NULLIFY (matrix_ks_aux_fit(ispin, img)%matrix)
2043 2756 : ALLOCATE (matrix_ks_aux_fit(ispin, img)%matrix)
2044 2756 : CALL dbcsr_create(matrix_ks_aux_fit(ispin, img)%matrix, template=matrix_s_aux_fit(1, 1)%matrix)
2045 2756 : CALL dbcsr_copy(matrix_ks_aux_fit(ispin, img)%matrix, matrix_ks_aux_fit_hfx(ispin, img)%matrix)
2046 : CALL dbcsr_add(matrix_ks_aux_fit(ispin, img)%matrix, matrix_ks_aux_fit_dft(ispin, img)%matrix, &
2047 5282 : 1.0_dp, -admm_env%gsi(ispin)**(2.0_dp/3.0_dp))
2048 : END DO
2049 : END DO
2050 : END IF
2051 :
2052 : ! the temporary DBCSR matrices for the rskp_transform we have to manually allocate
2053 726 : ALLOCATE (ksmatrix(2))
2054 : CALL dbcsr_create(ksmatrix(1), template=matrix_ks_aux_fit(1, 1)%matrix, &
2055 242 : matrix_type=dbcsr_type_symmetric)
2056 : CALL dbcsr_create(ksmatrix(2), template=matrix_ks_aux_fit(1, 1)%matrix, &
2057 242 : matrix_type=dbcsr_type_antisymmetric)
2058 : CALL dbcsr_create(tmpmatrix_ks, template=matrix_ks_aux_fit(1, 1)%matrix, &
2059 242 : matrix_type=dbcsr_type_symmetric)
2060 242 : CALL cp_dbcsr_alloc_block_from_nbl(ksmatrix(1), sab_aux_fit)
2061 242 : CALL cp_dbcsr_alloc_block_from_nbl(ksmatrix(2), sab_aux_fit)
2062 :
2063 242 : kplocal = kp_range(2) - kp_range(1) + 1
2064 660 : kpmax = MAXVAL(kp_dist(2, :) - kp_dist(1, :) + 1)
2065 242 : para_env => kpoints%blacs_env_all%para_env
2066 :
2067 : CALL cp_fm_struct_create(struct_aux_aux, context=kpoints%blacs_env, para_env=kpoints%para_env_kp, &
2068 242 : nrow_global=nao_aux_fit, ncol_global=nao_aux_fit)
2069 242 : CALL cp_fm_create(work_aux_aux, struct_aux_aux)
2070 242 : CALL cp_fm_create(work_aux_aux2, struct_aux_aux)
2071 :
2072 : CALL cp_fm_struct_create(struct_aux_orb, context=kpoints%blacs_env, para_env=kpoints%para_env_kp, &
2073 242 : nrow_global=nao_aux_fit, ncol_global=nao_orb)
2074 242 : CALL cp_fm_create(work_aux_orb, struct_aux_orb)
2075 :
2076 : CALL cp_fm_struct_create(struct_orb_orb, context=kpoints%blacs_env, para_env=kpoints%para_env_kp, &
2077 242 : nrow_global=nao_orb, ncol_global=nao_orb)
2078 :
2079 : !Create cfm work matrices
2080 242 : IF (.NOT. use_real_wfn) THEN
2081 242 : CALL cp_cfm_create(cK, struct_aux_aux)
2082 :
2083 242 : CALL cp_cfm_create(cwork_aux_orb, struct_aux_orb)
2084 :
2085 242 : CALL cp_cfm_create(cwork_orb_orb, struct_orb_orb)
2086 : END IF
2087 :
2088 : !We create the fms in which we store the KS ORB matrix at each kp
2089 4636 : ALLOCATE (fm_ks(kplocal, 2, nspins))
2090 508 : DO ispin = 1, nspins
2091 1040 : DO i = 1, 2
2092 3668 : DO ikp = 1, kplocal
2093 3402 : CALL cp_fm_create(fm_ks(ikp, i, ispin), struct_orb_orb)
2094 : END DO
2095 : END DO
2096 : END DO
2097 :
2098 242 : CALL cp_fm_struct_release(struct_aux_aux)
2099 242 : CALL cp_fm_struct_release(struct_aux_orb)
2100 242 : CALL cp_fm_struct_release(struct_orb_orb)
2101 :
2102 8270 : ALLOCATE (info(nkp*nspins, 2))
2103 242 : indx = 0
2104 1594 : DO ikp = 1, kpmax
2105 3074 : DO ispin = 1, nspins
2106 5726 : DO igroup = 1, nkp_groups
2107 : ! number of current kpoint
2108 2894 : ik = kp_dist(1, igroup) + ikp - 1
2109 2894 : IF (ik > kp_dist(2, igroup)) CYCLE
2110 2804 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2111 2804 : indx = indx + 1
2112 :
2113 2804 : IF (use_real_wfn) THEN
2114 0 : CALL dbcsr_set(ksmatrix(1), 0.0_dp)
2115 : CALL rskp_transform(rmatrix=ksmatrix(1), rsmat=matrix_ks_aux_fit, ispin=ispin, &
2116 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
2117 0 : CALL dbcsr_desymmetrize(ksmatrix(1), tmpmatrix_ks)
2118 0 : CALL copy_dbcsr_to_fm(tmpmatrix_ks, admm_env%work_aux_aux)
2119 : ELSE
2120 2804 : CALL dbcsr_set(ksmatrix(1), 0.0_dp)
2121 2804 : CALL dbcsr_set(ksmatrix(2), 0.0_dp)
2122 : CALL rskp_transform(rmatrix=ksmatrix(1), cmatrix=ksmatrix(2), rsmat=matrix_ks_aux_fit, ispin=ispin, &
2123 2804 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
2124 2804 : CALL dbcsr_desymmetrize(ksmatrix(1), tmpmatrix_ks)
2125 2804 : CALL copy_dbcsr_to_fm(tmpmatrix_ks, admm_env%work_aux_aux)
2126 2804 : CALL dbcsr_desymmetrize(ksmatrix(2), tmpmatrix_ks)
2127 2804 : CALL copy_dbcsr_to_fm(tmpmatrix_ks, admm_env%work_aux_aux2)
2128 : END IF
2129 :
2130 4284 : IF (my_kpgrp) THEN
2131 1435 : CALL cp_fm_start_copy_general(admm_env%work_aux_aux, work_aux_aux, para_env, info(indx, 1))
2132 1435 : IF (.NOT. use_real_wfn) THEN
2133 : CALL cp_fm_start_copy_general(admm_env%work_aux_aux2, work_aux_aux2, &
2134 1435 : para_env, info(indx, 2))
2135 : END IF
2136 : ELSE
2137 1369 : CALL cp_fm_start_copy_general(admm_env%work_aux_aux, fmdummy, para_env, info(indx, 1))
2138 1369 : IF (.NOT. use_real_wfn) THEN
2139 1369 : CALL cp_fm_start_copy_general(admm_env%work_aux_aux2, fmdummy, para_env, info(indx, 2))
2140 : END IF
2141 : END IF
2142 : END DO
2143 : END DO
2144 : END DO
2145 :
2146 : indx = 0
2147 1594 : DO ikp = 1, kpmax
2148 3074 : DO ispin = 1, nspins
2149 4374 : DO igroup = 1, nkp_groups
2150 : ! number of current kpoint
2151 2894 : ik = kp_dist(1, igroup) + ikp - 1
2152 2894 : IF (ik > kp_dist(2, igroup)) CYCLE
2153 2804 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2154 1369 : indx = indx + 1
2155 1480 : IF (my_kpgrp) THEN
2156 1435 : CALL cp_fm_finish_copy_general(work_aux_aux, info(indx, 1))
2157 1435 : IF (.NOT. use_real_wfn) THEN
2158 1435 : CALL cp_fm_finish_copy_general(work_aux_aux2, info(indx, 2))
2159 1435 : CALL cp_fm_to_cfm(work_aux_aux, work_aux_aux2, cK)
2160 : END IF
2161 : END IF
2162 : END DO
2163 :
2164 1480 : IF (ikp > kplocal) CYCLE
2165 1435 : kp => kpoints%kp_aux_env(ikp)%kpoint_env
2166 2787 : IF (use_real_wfn) THEN
2167 :
2168 : !! K*A
2169 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
2170 : 1.0_dp, work_aux_aux, kp%amat, 0.0_dp, &
2171 0 : work_aux_orb)
2172 : !! A^T*K*A
2173 : CALL parallel_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, &
2174 : 1.0_dp, kp%amat, work_aux_orb, 0.0_dp, &
2175 0 : fm_ks(ikp, 1, ispin))
2176 : ELSE
2177 :
2178 1435 : IF (admm_env%do_admmq .OR. admm_env%do_admms) THEN
2179 : !Need to subdtract lambda* S_aux to K_aux, and scale the whole thing by gsi
2180 1022 : fac = CMPLX(-admm_env%lambda_merlot(ispin), 0.0_dp, dp)
2181 1022 : CALL cp_cfm_scale_and_add(z_one, cK, fac, kp%csmat)
2182 1022 : CALL cp_cfm_scale(admm_env%gsi(ispin), cK)
2183 : END IF
2184 :
2185 1435 : IF (admm_env%do_admmp) THEN
2186 : !Need to substract labda*gsi*S_aux to gsi**2*K_aux
2187 98 : fac = CMPLX(-admm_env%gsi(ispin)*admm_env%lambda_merlot(ispin), 0.0_dp, dp)
2188 98 : fac2 = CMPLX(admm_env%gsi(ispin)**2, 0.0_dp, dp)
2189 98 : CALL cp_cfm_scale_and_add(fac2, cK, fac, kp%csmat)
2190 : END IF
2191 :
2192 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
2193 1435 : z_one, cK, kp%camat, z_zero, cwork_aux_orb)
2194 :
2195 : CALL parallel_gemm('C', 'N', nao_orb, nao_orb, nao_aux_fit, &
2196 1435 : z_one, kp%camat, cwork_aux_orb, z_zero, cwork_orb_orb)
2197 :
2198 1435 : CALL cp_cfm_to_fm(cwork_orb_orb, mtargetr=fm_ks(ikp, 1, ispin), mtargeti=fm_ks(ikp, 2, ispin))
2199 : END IF
2200 : END DO
2201 : END DO
2202 :
2203 3046 : DO indx = 1, SIZE(info, 1)
2204 2804 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
2205 3046 : IF (.NOT. use_real_wfn) CALL cp_fm_cleanup_copy_general(info(indx, 2))
2206 : END DO
2207 :
2208 5850 : DEALLOCATE (info)
2209 242 : CALL dbcsr_release(ksmatrix(1))
2210 242 : CALL dbcsr_release(ksmatrix(2))
2211 242 : CALL dbcsr_release(tmpmatrix_ks)
2212 :
2213 242 : CALL cp_fm_release(work_aux_aux)
2214 242 : CALL cp_fm_release(work_aux_aux2)
2215 242 : CALL cp_fm_release(work_aux_orb)
2216 242 : IF (.NOT. use_real_wfn) THEN
2217 242 : CALL cp_cfm_release(cK)
2218 242 : CALL cp_cfm_release(cwork_aux_orb)
2219 242 : CALL cp_cfm_release(cwork_orb_orb)
2220 : END IF
2221 :
2222 242 : NULLIFY (matrix_k_tilde)
2223 :
2224 242 : CALL dbcsr_allocate_matrix_set(matrix_k_tilde, dft_control%nspins, dft_control%nimages)
2225 :
2226 508 : DO ispin = 1, nspins
2227 17700 : DO img = 1, dft_control%nimages
2228 17192 : ALLOCATE (matrix_k_tilde(ispin, img)%matrix)
2229 : CALL dbcsr_create(matrix=matrix_k_tilde(ispin, img)%matrix, template=matrix_ks_kp(1, 1)%matrix, &
2230 : name='MATRIX K_tilde '//TRIM(ADJUSTL(cp_to_string(ispin)))//'_'//TRIM(ADJUSTL(cp_to_string(img))), &
2231 17192 : matrix_type=dbcsr_type_symmetric)
2232 17192 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_k_tilde(ispin, img)%matrix, sab_kp)
2233 17458 : CALL dbcsr_set(matrix_k_tilde(ispin, img)%matrix, 0.0_dp)
2234 : END DO
2235 : END DO
2236 :
2237 242 : CALL cp_fm_get_info(admm_env%work_orb_orb, matrix_struct=struct_orb_orb)
2238 726 : ALLOCATE (fmwork(2))
2239 242 : CALL cp_fm_create(fmwork(1), struct_orb_orb)
2240 242 : CALL cp_fm_create(fmwork(2), struct_orb_orb)
2241 :
2242 : ! reuse the density transform to FT the KS matrix
2243 : CALL kpoint_density_transform(kpoints, matrix_k_tilde, .FALSE., &
2244 : matrix_k_tilde(1, 1)%matrix, sab_kp, &
2245 242 : fmwork, for_aux_fit=.FALSE., pmat_ext=fm_ks, overlap_rs=matrix_s)
2246 242 : CALL cp_fm_release(fmwork(1))
2247 242 : CALL cp_fm_release(fmwork(2))
2248 :
2249 508 : DO ispin = 1, nspins
2250 1040 : DO i = 1, 2
2251 3668 : DO ikp = 1, kplocal
2252 3402 : CALL cp_fm_release(fm_ks(ikp, i, ispin))
2253 : END DO
2254 : END DO
2255 : END DO
2256 :
2257 508 : DO ispin = 1, nspins
2258 17700 : DO img = 1, dft_control%nimages
2259 17192 : CALL dbcsr_add(matrix_ks_kp(ispin, img)%matrix, matrix_k_tilde(ispin, img)%matrix, 1.0_dp, 1.0_dp)
2260 17458 : IF (admm_env%do_admmq .OR. admm_env%do_admmp .OR. admm_env%do_admms) THEN
2261 : !In ADMMQ and ADMMP, need to add lambda*S_orb (Merlot eq 27)
2262 : CALL dbcsr_add(matrix_ks_kp(ispin, img)%matrix, matrix_s(1, img)%matrix, &
2263 3816 : 1.0_dp, admm_env%lambda_merlot(ispin))
2264 : END IF
2265 : END DO
2266 : END DO
2267 :
2268 : !Scale the energies
2269 242 : IF (admm_env%do_admmp) THEN
2270 14 : IF (nspins == 1) THEN
2271 14 : energy%exc_aux_fit = (admm_env%gsi(1))**2.0_dp*energy%exc_aux_fit
2272 14 : energy%exc1_aux_fit = (admm_env%gsi(1))**2.0_dp*energy%exc1_aux_fit
2273 14 : energy%ex = (admm_env%gsi(1))**2.0_dp*energy%ex
2274 : ELSE
2275 0 : energy%exc_aux_fit = 0.0_dp
2276 0 : energy%exc1_aux_fit = 0.0_dp
2277 0 : energy%ex = 0.0_dp
2278 0 : DO ispin = 1, dft_control%nspins
2279 0 : energy%exc_aux_fit = energy%exc_aux_fit + (admm_env%gsi(ispin))**2.0_dp*ener_x(ispin)
2280 0 : energy%exc1_aux_fit = energy%exc1_aux_fit + (admm_env%gsi(ispin))**2.0_dp*ener_x1(ispin)
2281 0 : energy%ex = energy%ex + (admm_env%gsi(ispin))**2.0_dp*ener_k(ispin)
2282 : END DO
2283 : END IF
2284 : END IF
2285 :
2286 : !Scale the energies and clean-up
2287 242 : IF (admm_env%do_admms) THEN
2288 70 : IF (nspins == 1) THEN
2289 56 : energy%exc_aux_fit = (admm_env%gsi(1))**(2.0_dp/3.0_dp)*energy%exc_aux_fit
2290 56 : energy%exc1_aux_fit = (admm_env%gsi(1))**(2.0_dp/3.0_dp)*energy%exc1_aux_fit
2291 : ELSE
2292 14 : energy%exc_aux_fit = 0.0_dp
2293 14 : energy%exc1_aux_fit = 0.0_dp
2294 42 : DO ispin = 1, nspins
2295 28 : energy%exc_aux_fit = energy%exc_aux_fit + (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)*ener_x(ispin)
2296 42 : energy%exc1_aux_fit = energy%exc1_aux_fit + (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)*ener_x1(ispin)
2297 : END DO
2298 : END IF
2299 :
2300 70 : CALL dbcsr_deallocate_matrix_set(matrix_ks_aux_fit)
2301 : END IF
2302 :
2303 242 : CALL dbcsr_deallocate_matrix_set(matrix_k_tilde)
2304 :
2305 242 : CALL timestop(handle)
2306 :
2307 968 : END SUBROUTINE merge_ks_matrix_none_kp
2308 :
2309 : ! **************************************************************************************************
2310 : !> \brief Calculate exchange correction energy (Merlot2014 Eqs. 32, 33) for every spin, for KP
2311 : !> \param qs_env ...
2312 : !> \param admm_env ...
2313 : !> \param ener_k_ispin exact ispin (Fock) exchange in auxiliary basis
2314 : !> \param ener_x_ispin ispin DFT exchange in auxiliary basis
2315 : !> \param ener_x1_ispin ispin DFT exchange in auxiliary basis, due to the GAPW atomic contributions
2316 : !> \param ispin ...
2317 : ! **************************************************************************************************
2318 404 : SUBROUTINE calc_spin_dep_aux_exch_ener(qs_env, admm_env, ener_k_ispin, ener_x_ispin, &
2319 : ener_x1_ispin, ispin)
2320 : TYPE(qs_environment_type), POINTER :: qs_env
2321 : TYPE(admm_type), POINTER :: admm_env
2322 : REAL(dp), INTENT(INOUT) :: ener_k_ispin, ener_x_ispin, ener_x1_ispin
2323 : INTEGER, INTENT(IN) :: ispin
2324 :
2325 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_spin_dep_aux_exch_ener'
2326 :
2327 : CHARACTER(LEN=default_string_length) :: basis_type
2328 : INTEGER :: handle, img, myspin, nimg
2329 : LOGICAL :: gapw
2330 : REAL(dp) :: tmp
2331 404 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
2332 : TYPE(admm_gapw_r3d_rs_type), POINTER :: admm_gapw_env
2333 404 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2334 404 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
2335 404 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit_hfx, rho_ao_aux, &
2336 404 : rho_ao_aux_buffer
2337 : TYPE(dft_control_type), POINTER :: dft_control
2338 : TYPE(local_rho_type), POINTER :: local_rho_buffer
2339 : TYPE(mp_para_env_type), POINTER :: para_env
2340 404 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
2341 404 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, v_rspace_dummy, v_tau_rspace_dummy
2342 : TYPE(qs_ks_env_type), POINTER :: ks_env
2343 : TYPE(qs_rho_type), POINTER :: rho_aux_fit, rho_aux_fit_buffer
2344 : TYPE(section_vals_type), POINTER :: xc_section_aux
2345 : TYPE(task_list_type), POINTER :: task_list
2346 :
2347 404 : CALL timeset(routineN, handle)
2348 :
2349 404 : NULLIFY (ks_env, rho_aux_fit, rho_aux_fit_buffer, rho_ao, &
2350 404 : xc_section_aux, v_rspace_dummy, v_tau_rspace_dummy, &
2351 404 : rho_ao_aux, rho_ao_aux_buffer, dft_control, &
2352 404 : matrix_ks_aux_fit_hfx, task_list, local_rho_buffer, admm_gapw_env)
2353 :
2354 404 : NULLIFY (rho_g, rho_r, tot_rho_r)
2355 :
2356 404 : CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control)
2357 : CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, rho_aux_fit_buffer=rho_aux_fit_buffer, &
2358 404 : matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx)
2359 :
2360 : CALL qs_rho_get(rho_aux_fit, &
2361 404 : rho_ao_kp=rho_ao_aux)
2362 :
2363 : CALL qs_rho_get(rho_aux_fit_buffer, &
2364 : rho_ao_kp=rho_ao_aux_buffer, &
2365 : rho_g=rho_g, &
2366 : rho_r=rho_r, &
2367 404 : tot_rho_r=tot_rho_r)
2368 :
2369 404 : gapw = admm_env%do_gapw
2370 404 : nimg = dft_control%nimages
2371 :
2372 : ! Calculate rho_buffer = rho_aux(ispin) to get exchange of ispin electrons
2373 1240 : DO img = 1, nimg
2374 836 : CALL dbcsr_set(rho_ao_aux_buffer(1, img)%matrix, 0.0_dp)
2375 836 : CALL dbcsr_set(rho_ao_aux_buffer(2, img)%matrix, 0.0_dp)
2376 : CALL dbcsr_add(rho_ao_aux_buffer(ispin, img)%matrix, &
2377 1240 : rho_ao_aux(ispin, img)%matrix, 0.0_dp, 1.0_dp)
2378 : END DO
2379 :
2380 : ! By default use standard AUX_FIT basis and task_list. IF GAPW use the soft ones
2381 : basis_type = "AUX_FIT"
2382 404 : task_list => admm_env%task_list_aux_fit
2383 404 : IF (gapw) THEN
2384 : basis_type = "AUX_FIT_SOFT"
2385 124 : task_list => admm_env%admm_gapw_env%task_list
2386 : END IF
2387 :
2388 : ! integration for getting the spin dependent density has to done for both spins!
2389 1212 : DO myspin = 1, dft_control%nspins
2390 :
2391 808 : rho_ao => rho_ao_aux_buffer(myspin, :)
2392 : CALL calculate_rho_elec(ks_env=ks_env, &
2393 : matrix_p_kp=rho_ao, &
2394 : rho=rho_r(myspin), &
2395 : rho_gspace=rho_g(myspin), &
2396 : total_rho=tot_rho_r(myspin), &
2397 : soft_valid=.FALSE., &
2398 : basis_type="AUX_FIT", &
2399 1212 : task_list_external=task_list)
2400 :
2401 : END DO
2402 :
2403 : ! Write changes in buffer density matrix
2404 404 : CALL qs_rho_set(rho_aux_fit_buffer, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
2405 :
2406 404 : xc_section_aux => admm_env%xc_section_aux
2407 :
2408 : ener_x_ispin = 0.0_dp
2409 :
2410 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_aux_fit_buffer, xc_section=xc_section_aux, &
2411 : vxc_rho=v_rspace_dummy, vxc_tau=v_tau_rspace_dummy, exc=ener_x_ispin, &
2412 404 : just_energy=.TRUE.)
2413 :
2414 : !atomic contributions: use the atomic density as stored in admm_env%gapw_env
2415 404 : ener_x1_ispin = 0.0_dp
2416 404 : IF (gapw) THEN
2417 :
2418 124 : admm_gapw_env => admm_env%admm_gapw_env
2419 : CALL get_qs_env(qs_env, &
2420 : atomic_kind_set=atomic_kind_set, &
2421 124 : para_env=para_env)
2422 :
2423 124 : CALL local_rho_set_create(local_rho_buffer)
2424 : CALL allocate_rho_atom_internals(local_rho_buffer%rho_atom_set, atomic_kind_set, &
2425 124 : admm_gapw_env%admm_kind_set, dft_control, para_env)
2426 :
2427 : CALL calculate_rho_atom_coeff(qs_env, rho_ao_aux_buffer, &
2428 : rho_atom_set=local_rho_buffer%rho_atom_set, &
2429 : qs_kind_set=admm_gapw_env%admm_kind_set, &
2430 : oce=admm_gapw_env%oce, sab=admm_env%sab_aux_fit, &
2431 124 : para_env=para_env)
2432 :
2433 : CALL prepare_gapw_den(qs_env, local_rho_set=local_rho_buffer, do_rho0=.FALSE., &
2434 124 : kind_set_external=admm_gapw_env%admm_kind_set)
2435 :
2436 : CALL calculate_vxc_atom(qs_env, energy_only=.TRUE., exc1=ener_x1_ispin, &
2437 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
2438 : xc_section_external=xc_section_aux, &
2439 124 : rho_atom_set_external=local_rho_buffer%rho_atom_set)
2440 :
2441 124 : CALL local_rho_set_release(local_rho_buffer)
2442 : END IF
2443 :
2444 404 : ener_k_ispin = 0.0_dp
2445 :
2446 : !! ** Calculate the exchange energy
2447 1240 : DO img = 1, nimg
2448 836 : CALL dbcsr_dot(matrix_ks_aux_fit_hfx(ispin, img)%matrix, rho_ao_aux_buffer(ispin, img)%matrix, tmp)
2449 1240 : ener_k_ispin = ener_k_ispin + tmp
2450 : END DO
2451 :
2452 : ! Divide exchange for indivivual spin by two, since the ener_k_ispin originally is total
2453 : ! exchange of alpha and beta
2454 404 : ener_k_ispin = ener_k_ispin/2.0_dp
2455 :
2456 404 : CALL timestop(handle)
2457 :
2458 404 : END SUBROUTINE calc_spin_dep_aux_exch_ener
2459 :
2460 : ! **************************************************************************************************
2461 : !> \brief Scale density matrix by gsi(ispin), is needed for force scaling in ADMMP
2462 : !> \param qs_env ...
2463 : !> \param rho_ao_orb ...
2464 : !> \param scale_back ...
2465 : !> \author Jan Wilhelm, 12/2014
2466 : ! **************************************************************************************************
2467 640 : SUBROUTINE scale_dm(qs_env, rho_ao_orb, scale_back)
2468 : TYPE(qs_environment_type), POINTER :: qs_env
2469 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_orb
2470 : LOGICAL, INTENT(IN) :: scale_back
2471 :
2472 : CHARACTER(LEN=*), PARAMETER :: routineN = 'scale_dm'
2473 :
2474 : INTEGER :: handle, img, ispin
2475 : TYPE(admm_type), POINTER :: admm_env
2476 : TYPE(dft_control_type), POINTER :: dft_control
2477 :
2478 640 : CALL timeset(routineN, handle)
2479 :
2480 640 : NULLIFY (admm_env, dft_control)
2481 :
2482 : CALL get_qs_env(qs_env, &
2483 : admm_env=admm_env, &
2484 640 : dft_control=dft_control)
2485 :
2486 : ! only for ADMMP
2487 640 : IF (admm_env%do_admmp) THEN
2488 72 : DO ispin = 1, dft_control%nspins
2489 296 : DO img = 1, dft_control%nimages
2490 264 : IF (scale_back) THEN
2491 112 : CALL dbcsr_scale(rho_ao_orb(ispin, img)%matrix, 1.0_dp/admm_env%gsi(ispin))
2492 : ELSE
2493 112 : CALL dbcsr_scale(rho_ao_orb(ispin, img)%matrix, admm_env%gsi(ispin))
2494 : END IF
2495 : END DO
2496 : END DO
2497 : END IF
2498 :
2499 640 : CALL timestop(handle)
2500 :
2501 640 : END SUBROUTINE scale_dm
2502 :
2503 : ! **************************************************************************************************
2504 : !> \brief ...
2505 : !> \param ispin ...
2506 : !> \param admm_env ...
2507 : !> \param mo_set ...
2508 : !> \param mo_coeff_aux_fit ...
2509 : ! **************************************************************************************************
2510 230 : SUBROUTINE calc_aux_mo_derivs_none(ispin, admm_env, mo_set, mo_coeff_aux_fit)
2511 : INTEGER, INTENT(IN) :: ispin
2512 : TYPE(admm_type), POINTER :: admm_env
2513 : TYPE(mo_set_type), INTENT(IN) :: mo_set
2514 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff_aux_fit
2515 :
2516 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_aux_mo_derivs_none'
2517 :
2518 : INTEGER :: handle, nao_aux_fit, nao_orb, nmo
2519 230 : REAL(dp), DIMENSION(:), POINTER :: occupation_numbers, scaling_factor
2520 230 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit, &
2521 230 : matrix_ks_aux_fit_dft, &
2522 230 : matrix_ks_aux_fit_hfx
2523 : TYPE(dbcsr_type) :: dbcsr_work
2524 :
2525 230 : NULLIFY (matrix_ks_aux_fit, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx)
2526 :
2527 230 : CALL timeset(routineN, handle)
2528 :
2529 230 : nao_aux_fit = admm_env%nao_aux_fit
2530 230 : nao_orb = admm_env%nao_orb
2531 230 : nmo = admm_env%nmo(ispin)
2532 :
2533 : CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit, &
2534 : matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx, &
2535 230 : matrix_ks_aux_fit_dft=matrix_ks_aux_fit_dft)
2536 :
2537 : ! just calculate the mo derivs in the aux basis
2538 : ! only needs to be done on the converged ks matrix for the force calc
2539 : ! Note with OT and purification NONE, the merging of the derivs
2540 : ! happens implicitly because the KS matrices have been already been merged
2541 : ! and adding them here would be double counting.
2542 :
2543 230 : IF (admm_env%do_admms) THEN
2544 : !In ADMMS, we use the K matrix defined as K_hf - gsi^2/3*K_dft
2545 12 : CALL dbcsr_create(dbcsr_work, template=matrix_ks_aux_fit(ispin)%matrix)
2546 12 : CALL dbcsr_copy(dbcsr_work, matrix_ks_aux_fit_hfx(ispin)%matrix)
2547 12 : CALL dbcsr_add(dbcsr_work, matrix_ks_aux_fit_dft(ispin)%matrix, 1.0_dp, -admm_env%gsi(ispin)**(2.0_dp/3.0_dp))
2548 12 : CALL copy_dbcsr_to_fm(dbcsr_work, admm_env%K(ispin))
2549 12 : CALL dbcsr_release(dbcsr_work)
2550 : ELSE
2551 218 : CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
2552 : END IF
2553 230 : CALL cp_fm_uplo_to_full(admm_env%K(ispin), admm_env%work_aux_aux)
2554 :
2555 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_aux_fit, &
2556 : 1.0_dp, admm_env%K(ispin), mo_coeff_aux_fit, 0.0_dp, &
2557 230 : admm_env%H(ispin))
2558 :
2559 230 : CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
2560 690 : ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
2561 :
2562 2194 : scaling_factor = 2.0_dp*occupation_numbers
2563 :
2564 230 : CALL cp_fm_column_scale(admm_env%H(ispin), scaling_factor)
2565 :
2566 230 : DEALLOCATE (scaling_factor)
2567 :
2568 230 : CALL timestop(handle)
2569 :
2570 230 : END SUBROUTINE calc_aux_mo_derivs_none
2571 :
2572 : ! **************************************************************************************************
2573 : !> \brief ...
2574 : !> \param ispin ...
2575 : !> \param admm_env ...
2576 : !> \param mo_set ...
2577 : !> \param mo_derivs ...
2578 : !> \param matrix_ks_aux_fit ...
2579 : ! **************************************************************************************************
2580 100 : SUBROUTINE merge_mo_derivs_no_diag(ispin, admm_env, mo_set, mo_derivs, matrix_ks_aux_fit)
2581 : INTEGER, INTENT(IN) :: ispin
2582 : TYPE(admm_type), POINTER :: admm_env
2583 : TYPE(mo_set_type), INTENT(IN) :: mo_set
2584 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_derivs
2585 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit
2586 :
2587 : CHARACTER(LEN=*), PARAMETER :: routineN = 'merge_mo_derivs_no_diag'
2588 :
2589 : INTEGER :: handle, nao_aux_fit, nao_orb, nmo
2590 100 : REAL(dp), DIMENSION(:), POINTER :: occupation_numbers, scaling_factor
2591 :
2592 100 : CALL timeset(routineN, handle)
2593 :
2594 100 : nao_aux_fit = admm_env%nao_aux_fit
2595 100 : nao_orb = admm_env%nao_orb
2596 100 : nmo = admm_env%nmo(ispin)
2597 :
2598 100 : CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
2599 100 : CALL cp_fm_uplo_to_full(admm_env%K(ispin), admm_env%work_aux_aux)
2600 :
2601 100 : CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
2602 300 : ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
2603 460 : scaling_factor = 0.5_dp
2604 :
2605 : !! ** calculate first part
2606 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nmo, &
2607 : 1.0_dp, admm_env%C_hat(ispin), admm_env%lambda_inv(ispin), 0.0_dp, &
2608 100 : admm_env%work_aux_nmo(ispin))
2609 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_aux_fit, &
2610 : 1.0_dp, admm_env%K(ispin), admm_env%work_aux_nmo(ispin), 0.0_dp, &
2611 100 : admm_env%work_aux_nmo2(ispin))
2612 : CALL parallel_gemm('T', 'N', nao_orb, nmo, nao_aux_fit, &
2613 : 2.0_dp, admm_env%A, admm_env%work_aux_nmo2(ispin), 0.0_dp, &
2614 100 : admm_env%mo_derivs_tmp(ispin))
2615 : !! ** calculate second part
2616 : CALL parallel_gemm('T', 'N', nmo, nmo, nao_aux_fit, &
2617 : 1.0_dp, admm_env%work_aux_nmo(ispin), admm_env%work_aux_nmo2(ispin), 0.0_dp, &
2618 100 : admm_env%work_orb_orb)
2619 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nmo, &
2620 : 1.0_dp, admm_env%C_hat(ispin), admm_env%work_orb_orb, 0.0_dp, &
2621 100 : admm_env%work_aux_orb)
2622 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_aux_fit, &
2623 : 1.0_dp, admm_env%S, admm_env%work_aux_orb, 0.0_dp, &
2624 100 : admm_env%work_aux_nmo(ispin))
2625 : CALL parallel_gemm('T', 'N', nao_orb, nmo, nao_aux_fit, &
2626 : -2.0_dp, admm_env%A, admm_env%work_aux_nmo(ispin), 1.0_dp, &
2627 100 : admm_env%mo_derivs_tmp(ispin))
2628 :
2629 100 : CALL cp_fm_column_scale(admm_env%mo_derivs_tmp(ispin), scaling_factor)
2630 :
2631 100 : CALL cp_fm_scale_and_add(1.0_dp, mo_derivs(ispin), 1.0_dp, admm_env%mo_derivs_tmp(ispin))
2632 :
2633 100 : DEALLOCATE (scaling_factor)
2634 :
2635 100 : CALL timestop(handle)
2636 :
2637 100 : END SUBROUTINE merge_mo_derivs_no_diag
2638 :
2639 : ! **************************************************************************************************
2640 : !> \brief Calculate the derivative of the AUX_FIT mo, based on the ORB mo_derivs
2641 : !> \param qs_env ...
2642 : !> \param mo_derivs the MO derivatives in the orbital basis
2643 : ! **************************************************************************************************
2644 6864 : SUBROUTINE calc_admm_mo_derivatives(qs_env, mo_derivs)
2645 :
2646 : TYPE(qs_environment_type), POINTER :: qs_env
2647 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs
2648 :
2649 : INTEGER :: ispin, nspins
2650 : TYPE(admm_type), POINTER :: admm_env
2651 6864 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_derivs_fm
2652 6864 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_derivs_aux_fit
2653 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
2654 6864 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit
2655 6864 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array, mos_aux_fit
2656 :
2657 6864 : NULLIFY (mo_array, mos_aux_fit, matrix_ks_aux_fit, mo_coeff_aux_fit, &
2658 6864 : mo_derivs_aux_fit, mo_coeff)
2659 :
2660 6864 : CALL get_qs_env(qs_env, admm_env=admm_env, mos=mo_array)
2661 : CALL get_admm_env(admm_env, mos_aux_fit=mos_aux_fit, mo_derivs_aux_fit=mo_derivs_aux_fit, &
2662 6864 : matrix_ks_aux_fit=matrix_ks_aux_fit)
2663 :
2664 6864 : nspins = SIZE(mo_derivs)
2665 28650 : ALLOCATE (mo_derivs_fm(nspins))
2666 14922 : DO ispin = 1, nspins
2667 8058 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff)
2668 14922 : CALL cp_fm_create(mo_derivs_fm(ispin), mo_coeff%matrix_struct)
2669 : END DO
2670 :
2671 14922 : DO ispin = 1, nspins
2672 8058 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff)
2673 8058 : CALL get_mo_set(mo_set=mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
2674 :
2675 8058 : CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, mo_derivs_fm(ispin))
2676 : CALL admm_mo_merge_derivs(ispin, admm_env, mo_array(ispin), mo_coeff, mo_coeff_aux_fit, &
2677 8058 : mo_derivs_fm, mo_derivs_aux_fit, matrix_ks_aux_fit)
2678 14922 : CALL copy_fm_to_dbcsr(mo_derivs_fm(ispin), mo_derivs(ispin)%matrix)
2679 : END DO
2680 :
2681 6864 : CALL cp_fm_release(mo_derivs_fm)
2682 :
2683 13728 : END SUBROUTINE calc_admm_mo_derivatives
2684 :
2685 : ! **************************************************************************************************
2686 : !> \brief Calculate the forces due to the AUX/ORB basis overlap in ADMM
2687 : !> \param qs_env ...
2688 : ! **************************************************************************************************
2689 288 : SUBROUTINE calc_admm_ovlp_forces(qs_env)
2690 : TYPE(qs_environment_type), POINTER :: qs_env
2691 :
2692 : INTEGER :: ispin
2693 : TYPE(admm_type), POINTER :: admm_env
2694 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
2695 288 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_aux_fit_vs_orb
2696 : TYPE(dft_control_type), POINTER :: dft_control
2697 288 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
2698 : TYPE(mo_set_type), POINTER :: mo_set
2699 :
2700 288 : CALL get_qs_env(qs_env, dft_control=dft_control)
2701 :
2702 288 : IF (dft_control%do_admm_dm) THEN
2703 0 : CPABORT("Forces with ADMM DM methods not implemented")
2704 : END IF
2705 288 : IF (dft_control%do_admm_mo .AND. .NOT. qs_env%run_rtp) THEN
2706 258 : NULLIFY (matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, mos_aux_fit, mos, admm_env)
2707 : CALL get_qs_env(qs_env=qs_env, &
2708 : mos=mos, &
2709 258 : admm_env=admm_env)
2710 : CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit, mos_aux_fit=mos_aux_fit, &
2711 258 : matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb)
2712 558 : DO ispin = 1, dft_control%nspins
2713 300 : mo_set => mos(ispin)
2714 300 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff)
2715 : ! if no purification we need to calculate the H matrix for forces
2716 558 : IF (admm_env%purification_method == do_admm_purify_none) THEN
2717 230 : CALL get_mo_set(mo_set=mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
2718 230 : CALL calc_aux_mo_derivs_none(ispin, qs_env%admm_env, mo_set, mo_coeff_aux_fit)
2719 : END IF
2720 : END DO
2721 258 : CALL calc_mixed_overlap_force(qs_env)
2722 : END IF
2723 :
2724 288 : END SUBROUTINE calc_admm_ovlp_forces
2725 :
2726 : ! **************************************************************************************************
2727 : !> \brief Calculate the forces due to the AUX/ORB basis overlap in ADMM, in the KP case
2728 : !> \param qs_env ...
2729 : ! **************************************************************************************************
2730 32 : SUBROUTINE calc_admm_ovlp_forces_kp(qs_env)
2731 : TYPE(qs_environment_type), POINTER :: qs_env
2732 :
2733 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_admm_ovlp_forces_kp'
2734 :
2735 : COMPLEX(dp) :: fac, fac2
2736 : INTEGER :: handle, i, igroup, ik, ikp, img, indx, &
2737 : ispin, kplocal, kpmax, nao_aux_fit, &
2738 : nao_orb, natom, nimg, nkp, nkp_groups, &
2739 : nspins
2740 : INTEGER, DIMENSION(2) :: kp_range
2741 32 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
2742 32 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2743 : LOGICAL :: gapw, my_kpgrp, use_real_wfn
2744 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: admm_force
2745 32 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
2746 : TYPE(admm_type), POINTER :: admm_env
2747 32 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2748 32 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
2749 : TYPE(cp_cfm_type) :: ckmatrix, cpmatrix, cS_inv, &
2750 : cwork_aux_aux, cwork_aux_orb, &
2751 : cwork_aux_orb2
2752 : TYPE(cp_fm_struct_type), POINTER :: struct_aux_aux, struct_aux_orb, &
2753 : struct_orb_orb
2754 : TYPE(cp_fm_type) :: fmdummy, S_inv, work_aux_aux, &
2755 : work_aux_aux2, work_aux_aux3, &
2756 : work_aux_orb
2757 32 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_skap, fm_skapa
2758 32 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
2759 32 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit, matrix_ks_aux_fit_dft, &
2760 32 : matrix_ks_aux_fit_hfx, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, matrix_skap, &
2761 32 : matrix_skapa, rho_ao_orb
2762 : TYPE(dbcsr_type) :: kmatrix_tmp
2763 32 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: kmatrix
2764 : TYPE(dft_control_type), POINTER :: dft_control
2765 : TYPE(kpoint_env_type), POINTER :: kp
2766 : TYPE(kpoint_type), POINTER :: kpoints
2767 : TYPE(mp_para_env_type), POINTER :: para_env
2768 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2769 32 : POINTER :: sab_aux_fit, sab_aux_fit_asymm, &
2770 32 : sab_aux_fit_vs_orb, sab_kp
2771 32 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2772 : TYPE(qs_ks_env_type), POINTER :: ks_env
2773 : TYPE(qs_rho_type), POINTER :: rho
2774 :
2775 32 : CALL timeset(routineN, handle)
2776 :
2777 : !Note: we only treat the case with purification none, there the overlap forces read as:
2778 : !F = 2*Tr[P * A^T * K_aux * S^-1_aux * Q^(x)] - 2*Tr[A * P * A^T * K_aux * S^-1_aux *S_aux^(x)]
2779 : !where P is the density matrix in the ORB basis. As a strategy, we FT all relevant matrices
2780 : !from real space to KP, calculate the matrix products, back FT to real space, and calculate the
2781 : !overlap forces
2782 :
2783 32 : NULLIFY (ks_env, admm_env, matrix_ks_aux_fit, &
2784 32 : matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, rho, force, &
2785 32 : para_env, atomic_kind_set, kpoints, sab_aux_fit, &
2786 32 : sab_aux_fit_vs_orb, sab_aux_fit_asymm, struct_orb_orb, &
2787 32 : struct_aux_orb, struct_aux_aux)
2788 :
2789 : CALL get_qs_env(qs_env, &
2790 : ks_env=ks_env, &
2791 : admm_env=admm_env, &
2792 : dft_control=dft_control, &
2793 : kpoints=kpoints, &
2794 : natom=natom, &
2795 : atomic_kind_set=atomic_kind_set, &
2796 : force=force, &
2797 32 : rho=rho)
2798 32 : nimg = dft_control%nimages
2799 : CALL get_admm_env(admm_env, &
2800 : matrix_s_aux_fit_kp=matrix_s_aux_fit, &
2801 : matrix_s_aux_fit_vs_orb_kp=matrix_s_aux_fit_vs_orb, &
2802 : sab_aux_fit=sab_aux_fit, &
2803 : sab_aux_fit_vs_orb=sab_aux_fit_vs_orb, &
2804 : sab_aux_fit_asymm=sab_aux_fit_asymm, &
2805 : matrix_ks_aux_fit_kp=matrix_ks_aux_fit, &
2806 : matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft, &
2807 32 : matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx)
2808 :
2809 32 : gapw = admm_env%do_gapw
2810 32 : nao_aux_fit = admm_env%nao_aux_fit
2811 32 : nao_orb = admm_env%nao_orb
2812 32 : nspins = dft_control%nspins
2813 :
2814 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
2815 : nkp_groups=nkp_groups, kp_dist=kp_dist, &
2816 32 : cell_to_index=cell_to_index, sab_nl=sab_kp)
2817 :
2818 : !Case study on ADMMQ, ADMMS and ADMMP
2819 32 : IF (admm_env%do_admms) THEN
2820 : !Here we buld the KS matrix: KS_hfx = gsi^2/3*KS_dft, the we then pass as the ususal KS_aux_fit
2821 6 : NULLIFY (matrix_ks_aux_fit)
2822 362 : ALLOCATE (matrix_ks_aux_fit(nspins, dft_control%nimages))
2823 146 : DO img = 1, dft_control%nimages
2824 344 : DO ispin = 1, nspins
2825 198 : NULLIFY (matrix_ks_aux_fit(ispin, img)%matrix)
2826 198 : ALLOCATE (matrix_ks_aux_fit(ispin, img)%matrix)
2827 198 : CALL dbcsr_create(matrix_ks_aux_fit(ispin, img)%matrix, template=matrix_s_aux_fit(1, 1)%matrix)
2828 198 : CALL dbcsr_copy(matrix_ks_aux_fit(ispin, img)%matrix, matrix_ks_aux_fit_hfx(ispin, img)%matrix)
2829 : CALL dbcsr_add(matrix_ks_aux_fit(ispin, img)%matrix, matrix_ks_aux_fit_dft(ispin, img)%matrix, &
2830 338 : 1.0_dp, -admm_env%gsi(ispin)**(2.0_dp/3.0_dp))
2831 : END DO
2832 : END DO
2833 : END IF
2834 :
2835 : ! the temporary DBCSR matrices for the rskp_transform we have to manually allocate
2836 : ! index 1 => real, index 2 => imaginary
2837 96 : ALLOCATE (kmatrix(2))
2838 : CALL dbcsr_create(kmatrix(1), template=matrix_ks_aux_fit(1, 1)%matrix, &
2839 32 : matrix_type=dbcsr_type_symmetric)
2840 : CALL dbcsr_create(kmatrix(2), template=matrix_ks_aux_fit(1, 1)%matrix, &
2841 32 : matrix_type=dbcsr_type_antisymmetric)
2842 : CALL dbcsr_create(kmatrix_tmp, template=matrix_ks_aux_fit(1, 1)%matrix, &
2843 32 : matrix_type=dbcsr_type_no_symmetry)
2844 32 : CALL cp_dbcsr_alloc_block_from_nbl(kmatrix(1), sab_aux_fit)
2845 32 : CALL cp_dbcsr_alloc_block_from_nbl(kmatrix(2), sab_aux_fit)
2846 :
2847 32 : kplocal = kp_range(2) - kp_range(1) + 1
2848 94 : kpmax = MAXVAL(kp_dist(2, :) - kp_dist(1, :) + 1)
2849 32 : para_env => kpoints%blacs_env_all%para_env
2850 1112 : ALLOCATE (info(nkp*nspins, 2))
2851 :
2852 : CALL cp_fm_struct_create(struct_aux_aux, context=kpoints%blacs_env, para_env=kpoints%para_env_kp, &
2853 32 : nrow_global=nao_aux_fit, ncol_global=nao_aux_fit)
2854 32 : CALL cp_fm_create(work_aux_aux, struct_aux_aux)
2855 32 : CALL cp_fm_create(work_aux_aux2, struct_aux_aux)
2856 32 : CALL cp_fm_create(work_aux_aux3, struct_aux_aux)
2857 32 : CALL cp_fm_create(s_inv, struct_aux_aux)
2858 :
2859 : CALL cp_fm_struct_create(struct_aux_orb, context=kpoints%blacs_env, para_env=kpoints%para_env_kp, &
2860 32 : nrow_global=nao_aux_fit, ncol_global=nao_orb)
2861 32 : CALL cp_fm_create(work_aux_orb, struct_aux_orb)
2862 :
2863 : CALL cp_fm_struct_create(struct_orb_orb, context=kpoints%blacs_env, para_env=kpoints%para_env_kp, &
2864 32 : nrow_global=nao_orb, ncol_global=nao_orb)
2865 :
2866 : !Create cfm work matrices
2867 32 : IF (.NOT. use_real_wfn) THEN
2868 32 : CALL cp_cfm_create(cpmatrix, struct_orb_orb)
2869 :
2870 32 : CALL cp_cfm_create(cS_inv, struct_aux_aux)
2871 32 : CALL cp_cfm_create(cwork_aux_aux, struct_aux_aux)
2872 32 : CALL cp_cfm_create(ckmatrix, struct_aux_aux)
2873 :
2874 32 : CALL cp_cfm_create(cwork_aux_orb, struct_aux_orb)
2875 32 : CALL cp_cfm_create(cwork_aux_orb2, struct_aux_orb)
2876 : END IF
2877 :
2878 : !We create the fms in which we store the KP matrix products
2879 1184 : ALLOCATE (fm_skap(kplocal, 2, nspins), fm_skapa(kplocal, 2, nspins))
2880 70 : DO ispin = 1, nspins
2881 146 : DO i = 1, 2
2882 496 : DO ikp = 1, kplocal
2883 382 : CALL cp_fm_create(fm_skap(ikp, i, ispin), struct_aux_orb)
2884 458 : CALL cp_fm_create(fm_skapa(ikp, i, ispin), struct_aux_aux)
2885 : END DO
2886 : END DO
2887 : END DO
2888 :
2889 32 : CALL cp_fm_struct_release(struct_aux_aux)
2890 32 : CALL cp_fm_struct_release(struct_aux_orb)
2891 32 : CALL cp_fm_struct_release(struct_orb_orb)
2892 :
2893 32 : indx = 0
2894 194 : DO ikp = 1, kpmax
2895 390 : DO ispin = 1, nspins
2896 748 : DO igroup = 1, nkp_groups
2897 : ! number of current kpoint
2898 390 : ik = kp_dist(1, igroup) + ikp - 1
2899 390 : IF (ik > kp_dist(2, igroup)) CYCLE
2900 380 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2901 380 : indx = indx + 1
2902 :
2903 : ! FT of matrices KS, then transfer to FM type
2904 380 : IF (use_real_wfn) THEN
2905 0 : CALL dbcsr_set(kmatrix(1), 0.0_dp)
2906 : CALL rskp_transform(rmatrix=kmatrix(1), rsmat=matrix_ks_aux_fit, ispin=ispin, &
2907 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
2908 0 : CALL dbcsr_desymmetrize(kmatrix(1), kmatrix_tmp)
2909 0 : CALL copy_dbcsr_to_fm(kmatrix_tmp, admm_env%work_aux_aux)
2910 : ELSE
2911 380 : CALL dbcsr_set(kmatrix(1), 0.0_dp)
2912 380 : CALL dbcsr_set(kmatrix(2), 0.0_dp)
2913 : CALL rskp_transform(rmatrix=kmatrix(1), cmatrix=kmatrix(2), rsmat=matrix_ks_aux_fit, ispin=ispin, &
2914 380 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
2915 380 : CALL dbcsr_desymmetrize(kmatrix(1), kmatrix_tmp)
2916 380 : CALL copy_dbcsr_to_fm(kmatrix_tmp, admm_env%work_aux_aux)
2917 380 : CALL dbcsr_desymmetrize(kmatrix(2), kmatrix_tmp)
2918 380 : CALL copy_dbcsr_to_fm(kmatrix_tmp, admm_env%work_aux_aux2)
2919 : END IF
2920 :
2921 576 : IF (my_kpgrp) THEN
2922 191 : CALL cp_fm_start_copy_general(admm_env%work_aux_aux, work_aux_aux, para_env, info(indx, 1))
2923 191 : IF (.NOT. use_real_wfn) THEN
2924 191 : CALL cp_fm_start_copy_general(admm_env%work_aux_aux2, work_aux_aux2, para_env, info(indx, 2))
2925 : END IF
2926 : ELSE
2927 189 : CALL cp_fm_start_copy_general(admm_env%work_aux_aux, fmdummy, para_env, info(indx, 1))
2928 189 : IF (.NOT. use_real_wfn) THEN
2929 189 : CALL cp_fm_start_copy_general(admm_env%work_aux_aux2, fmdummy, para_env, info(indx, 2))
2930 : END IF
2931 : END IF
2932 : END DO
2933 : END DO
2934 : END DO
2935 :
2936 : indx = 0
2937 194 : DO ikp = 1, kpmax
2938 390 : DO ispin = 1, nspins
2939 586 : DO igroup = 1, nkp_groups
2940 : ! number of current kpoint
2941 390 : ik = kp_dist(1, igroup) + ikp - 1
2942 390 : IF (ik > kp_dist(2, igroup)) CYCLE
2943 380 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2944 189 : indx = indx + 1
2945 196 : IF (my_kpgrp) THEN
2946 191 : CALL cp_fm_finish_copy_general(work_aux_aux, info(indx, 1))
2947 191 : IF (.NOT. use_real_wfn) THEN
2948 191 : CALL cp_fm_finish_copy_general(work_aux_aux2, info(indx, 2))
2949 191 : CALL cp_fm_to_cfm(work_aux_aux, work_aux_aux2, ckmatrix)
2950 : END IF
2951 : END IF
2952 : END DO
2953 196 : IF (ikp > kplocal) CYCLE
2954 191 : kp => kpoints%kp_aux_env(ikp)%kpoint_env
2955 :
2956 353 : IF (use_real_wfn) THEN
2957 :
2958 : !! Calculate S'_inverse
2959 0 : CALL cp_fm_to_fm(kp%smat, S_inv)
2960 0 : CALL cp_fm_cholesky_decompose(S_inv)
2961 0 : CALL cp_fm_cholesky_invert(S_inv)
2962 : !! Symmetrize the guy
2963 0 : CALL cp_fm_uplo_to_full(S_inv, work_aux_aux3)
2964 :
2965 : !We need to calculate S^-1*K*A*P and S^-1*K*A*P*A^T
2966 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit, 1.0_dp, S_inv, &
2967 0 : work_aux_aux, 0.0_dp, work_aux_aux3) ! S^-1 * K
2968 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, 1.0_dp, work_aux_aux3, &
2969 0 : kp%amat, 0.0_dp, work_aux_orb) ! S^-1 * K * A
2970 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, 1.0_dp, work_aux_orb, &
2971 : kpoints%kp_env(ikp)%kpoint_env%pmat(1, ispin), 0.0_dp, &
2972 0 : fm_skap(ikp, 1, ispin)) ! S^-1 * K * A * P
2973 : CALL parallel_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nao_orb, 1.0_dp, fm_skap(ikp, 1, ispin), &
2974 0 : kp%amat, 0.0_dp, fm_skapa(ikp, 1, ispin))
2975 :
2976 : ELSE !complex wfn
2977 :
2978 191 : IF (admm_env%do_admmq .OR. admm_env%do_admms) THEN
2979 : !Need to subdtract lambda* S_aux to K_aux, and scale the whole thing by gsi
2980 97 : fac = CMPLX(-admm_env%lambda_merlot(ispin), 0.0_dp, dp)
2981 97 : CALL cp_cfm_scale_and_add(z_one, ckmatrix, fac, kp%csmat)
2982 97 : CALL cp_cfm_scale(admm_env%gsi(ispin), ckmatrix)
2983 : END IF
2984 :
2985 191 : IF (admm_env%do_admmp) THEN
2986 : !Need to substract labda*gsi*S_aux to gsi**2*K_aux
2987 28 : fac = CMPLX(-admm_env%gsi(ispin)*admm_env%lambda_merlot(ispin), 0.0_dp, dp)
2988 28 : fac2 = CMPLX(admm_env%gsi(ispin)**2, 0.0_dp, dp)
2989 28 : CALL cp_cfm_scale_and_add(fac2, ckmatrix, fac, kp%csmat)
2990 : END IF
2991 :
2992 191 : CALL cp_cfm_to_cfm(kp%csmat, cS_inv)
2993 191 : CALL cp_cfm_cholesky_decompose(cS_inv)
2994 191 : CALL cp_cfm_cholesky_invert(cS_inv)
2995 191 : CALL cp_cfm_uplo_to_full(cS_inv, cwork_aux_aux)
2996 :
2997 : !Take the ORB density matrix from the kp_env
2998 : CALL cp_fm_to_cfm(kpoints%kp_env(ikp)%kpoint_env%pmat(1, ispin), &
2999 : kpoints%kp_env(ikp)%kpoint_env%pmat(2, ispin), &
3000 191 : cpmatrix)
3001 :
3002 : !Do the same thing as in the real case
3003 : !We need to calculate S^-1*K*A*P and S^-1*K*A*P*A^T
3004 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_aux_fit, nao_aux_fit, z_one, cS_inv, &
3005 191 : ckmatrix, z_zero, cwork_aux_aux) ! S^-1 * K
3006 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, z_one, cwork_aux_aux, &
3007 191 : kp%camat, z_zero, cwork_aux_orb) ! S^-1 * K * A
3008 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, z_one, cwork_aux_orb, &
3009 191 : cpmatrix, z_zero, cwork_aux_orb2) ! S^-1 * K * A * P
3010 : CALL parallel_gemm('N', 'C', nao_aux_fit, nao_aux_fit, nao_orb, z_one, cwork_aux_orb2, &
3011 191 : kp%camat, z_zero, cwork_aux_aux)
3012 :
3013 191 : IF (admm_env%do_admmq .OR. admm_env%do_admmp .OR. admm_env%do_admms) THEN
3014 : !In ADMMQ, ADMMS, and ADMMP, there is an extra lambda*Tq *P* Tq^T matrix to contract with S_aux^(x)
3015 : !we calculate it and add it to fm_skapa (aka cwork_aux_aux)
3016 :
3017 : !factor 0.5 because later multiplied by 2
3018 125 : fac = CMPLX(0.5_dp*admm_env%lambda_merlot(ispin)*admm_env%gsi(ispin), 0.0_dp, dp)
3019 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, z_one, kp%camat, cpmatrix, &
3020 125 : z_zero, cwork_aux_orb)
3021 : CALL parallel_gemm('N', 'C', nao_aux_fit, nao_aux_fit, nao_orb, fac, cwork_aux_orb, &
3022 125 : kp%camat, z_one, cwork_aux_aux)
3023 : END IF
3024 :
3025 191 : CALL cp_cfm_to_fm(cwork_aux_orb2, mtargetr=fm_skap(ikp, 1, ispin), mtargeti=fm_skap(ikp, 2, ispin))
3026 191 : CALL cp_cfm_to_fm(cwork_aux_aux, mtargetr=fm_skapa(ikp, 1, ispin), mtargeti=fm_skapa(ikp, 2, ispin))
3027 :
3028 : END IF
3029 :
3030 : END DO
3031 : END DO
3032 :
3033 412 : DO indx = 1, SIZE(info, 1)
3034 380 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
3035 412 : IF (.NOT. use_real_wfn) CALL cp_fm_cleanup_copy_general(info(indx, 2))
3036 : END DO
3037 :
3038 792 : DEALLOCATE (info)
3039 32 : CALL dbcsr_release(kmatrix(1))
3040 32 : CALL dbcsr_release(kmatrix(2))
3041 32 : CALL dbcsr_release(kmatrix_tmp)
3042 :
3043 32 : CALL cp_fm_release(work_aux_aux)
3044 32 : CALL cp_fm_release(work_aux_aux2)
3045 32 : CALL cp_fm_release(work_aux_aux3)
3046 32 : CALL cp_fm_release(S_inv)
3047 32 : CALL cp_fm_release(work_aux_orb)
3048 32 : IF (.NOT. use_real_wfn) THEN
3049 32 : CALL cp_cfm_release(ckmatrix)
3050 32 : CALL cp_cfm_release(cpmatrix)
3051 32 : CALL cp_cfm_release(cS_inv)
3052 32 : CALL cp_cfm_release(cwork_aux_aux)
3053 32 : CALL cp_cfm_release(cwork_aux_orb)
3054 32 : CALL cp_cfm_release(cwork_aux_orb2)
3055 : END IF
3056 :
3057 : !Back FT to real space
3058 10212 : ALLOCATE (matrix_skap(nspins, nimg), matrix_skapa(nspins, nimg))
3059 2494 : DO img = 1, nimg
3060 5042 : DO ispin = 1, nspins
3061 2548 : ALLOCATE (matrix_skap(ispin, img)%matrix)
3062 : CALL dbcsr_create(matrix_skap(ispin, img)%matrix, template=matrix_s_aux_fit_vs_orb(1, 1)%matrix, &
3063 2548 : matrix_type=dbcsr_type_no_symmetry)
3064 2548 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_skap(ispin, img)%matrix, sab_aux_fit_vs_orb)
3065 :
3066 2548 : ALLOCATE (matrix_skapa(ispin, img)%matrix)
3067 : CALL dbcsr_create(matrix_skapa(ispin, img)%matrix, template=matrix_s_aux_fit(1, 1)%matrix, &
3068 2548 : matrix_type=dbcsr_type_no_symmetry)
3069 5010 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_skapa(ispin, img)%matrix, sab_aux_fit_asymm)
3070 : END DO
3071 : END DO
3072 :
3073 96 : ALLOCATE (fmwork(2))
3074 32 : CALL cp_fm_get_info(admm_env%work_aux_orb, matrix_struct=struct_aux_orb)
3075 32 : CALL cp_fm_create(fmwork(1), struct_aux_orb)
3076 32 : CALL cp_fm_create(fmwork(2), struct_aux_orb)
3077 : CALL kpoint_density_transform(kpoints, matrix_skap, .FALSE., &
3078 : matrix_s_aux_fit_vs_orb(1, 1)%matrix, sab_aux_fit_vs_orb, &
3079 : fmwork, for_aux_fit=.TRUE., pmat_ext=fm_skap, &
3080 32 : overlap_rs=matrix_s_aux_fit, mixed_aux_orb=.TRUE.)
3081 32 : CALL cp_fm_release(fmwork(1))
3082 32 : CALL cp_fm_release(fmwork(2))
3083 :
3084 32 : CALL cp_fm_get_info(admm_env%work_aux_aux, matrix_struct=struct_aux_aux)
3085 32 : CALL cp_fm_create(fmwork(1), struct_aux_aux)
3086 32 : CALL cp_fm_create(fmwork(2), struct_aux_aux)
3087 : CALL kpoint_density_transform(kpoints, matrix_skapa, .FALSE., &
3088 : matrix_s_aux_fit(1, 1)%matrix, sab_aux_fit_asymm, &
3089 32 : fmwork, for_aux_fit=.TRUE., pmat_ext=fm_skapa, overlap_rs=matrix_s_aux_fit)
3090 32 : CALL cp_fm_release(fmwork(1))
3091 32 : CALL cp_fm_release(fmwork(2))
3092 32 : DEALLOCATE (fmwork)
3093 :
3094 2494 : DO img = 1, nimg
3095 5010 : DO ispin = 1, nspins
3096 2548 : CALL dbcsr_scale(matrix_skap(ispin, img)%matrix, -2.0_dp)
3097 5010 : CALL dbcsr_scale(matrix_skapa(ispin, img)%matrix, 2.0_dp)
3098 : END DO
3099 2494 : IF (nspins == 2) THEN
3100 86 : CALL dbcsr_add(matrix_skap(1, img)%matrix, matrix_skap(2, img)%matrix, 1.0_dp, 1.0_dp)
3101 86 : CALL dbcsr_add(matrix_skapa(1, img)%matrix, matrix_skapa(2, img)%matrix, 1.0_dp, 1.0_dp)
3102 : END IF
3103 : END DO
3104 :
3105 96 : ALLOCATE (admm_force(3, natom))
3106 32 : admm_force = 0.0_dp
3107 :
3108 32 : IF (admm_env%do_admmq .OR. admm_env%do_admmp .OR. admm_env%do_admms) THEN
3109 12 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_orb)
3110 310 : DO img = 1, nimg
3111 654 : DO ispin = 1, nspins
3112 654 : CALL dbcsr_scale(rho_ao_orb(ispin, img)%matrix, -admm_env%lambda_merlot(ispin))
3113 : END DO
3114 310 : IF (nspins == 2) CALL dbcsr_add(rho_ao_orb(1, img)%matrix, rho_ao_orb(2, img)%matrix, 1.0_dp, 1.0_dp)
3115 : END DO
3116 :
3117 : !In ADMMQ, ADMMS and ADMMP, there is an extra contribution from lambda*P_orb*S^(x)
3118 : CALL build_overlap_force(qs_env%ks_env, admm_force, basis_type_a="ORB", basis_type_b="ORB", &
3119 12 : sab_nl=sab_kp, matrixkp_p=rho_ao_orb(1, :))
3120 310 : DO img = 1, nimg
3121 298 : IF (nspins == 2) CALL dbcsr_add(rho_ao_orb(1, img)%matrix, rho_ao_orb(2, img)%matrix, 1.0_dp, -1.0_dp)
3122 686 : DO ispin = 1, nspins
3123 654 : CALL dbcsr_scale(rho_ao_orb(ispin, img)%matrix, -1.0_dp/admm_env%lambda_merlot(ispin))
3124 : END DO
3125 : END DO
3126 : END IF
3127 :
3128 : CALL build_overlap_force(qs_env%ks_env, admm_force, basis_type_a="AUX_FIT", basis_type_b="ORB", &
3129 32 : sab_nl=sab_aux_fit_vs_orb, matrixkp_p=matrix_skap(1, :))
3130 : CALL build_overlap_force(qs_env%ks_env, admm_force, basis_type_a="AUX_FIT", basis_type_b="AUX_FIT", &
3131 32 : sab_nl=sab_aux_fit_asymm, matrixkp_p=matrix_skapa(1, :))
3132 :
3133 32 : CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set)
3134 32 : DEALLOCATE (admm_force)
3135 :
3136 70 : DO ispin = 1, nspins
3137 146 : DO i = 1, 2
3138 496 : DO ikp = 1, kplocal
3139 382 : CALL cp_fm_release(fm_skap(ikp, i, ispin))
3140 458 : CALL cp_fm_release(fm_skapa(ikp, i, ispin))
3141 : END DO
3142 : END DO
3143 : END DO
3144 32 : CALL dbcsr_deallocate_matrix_set(matrix_skap)
3145 32 : CALL dbcsr_deallocate_matrix_set(matrix_skapa)
3146 :
3147 32 : IF (admm_env%do_admms) THEN
3148 6 : CALL dbcsr_deallocate_matrix_set(matrix_ks_aux_fit)
3149 : END IF
3150 :
3151 32 : CALL timestop(handle)
3152 :
3153 128 : END SUBROUTINE calc_admm_ovlp_forces_kp
3154 :
3155 : ! **************************************************************************************************
3156 : !> \brief Calculate derivatives terms from overlap matrices
3157 : !> \param qs_env ...
3158 : !> \param matrix_hz Fock matrix part using the response density in admm basis
3159 : !> \param matrix_pz response density in orbital basis
3160 : !> \param fval ...
3161 : ! **************************************************************************************************
3162 880 : SUBROUTINE admm_projection_derivative(qs_env, matrix_hz, matrix_pz, fval)
3163 : TYPE(qs_environment_type), POINTER :: qs_env
3164 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: matrix_hz, matrix_pz
3165 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: fval
3166 :
3167 : CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_projection_derivative'
3168 :
3169 : INTEGER :: handle, ispin, nao, natom, naux, nspins
3170 : REAL(KIND=dp) :: my_fval
3171 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: admm_force
3172 : TYPE(admm_type), POINTER :: admm_env
3173 880 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3174 880 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_aux_fit_vs_orb
3175 : TYPE(dbcsr_type), POINTER :: matrix_w_q, matrix_w_s
3176 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3177 880 : POINTER :: sab_aux_fit_asymm, sab_aux_fit_vs_orb
3178 880 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3179 : TYPE(qs_ks_env_type), POINTER :: ks_env
3180 :
3181 880 : CALL timeset(routineN, handle)
3182 :
3183 880 : CPASSERT(ASSOCIATED(qs_env))
3184 :
3185 880 : CALL get_qs_env(qs_env, ks_env=ks_env, admm_env=admm_env)
3186 : CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit, sab_aux_fit_asymm=sab_aux_fit_asymm, &
3187 880 : matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, sab_aux_fit_vs_orb=sab_aux_fit_vs_orb)
3188 :
3189 880 : my_fval = 2.0_dp
3190 880 : IF (PRESENT(fval)) my_fval = fval
3191 :
3192 880 : ALLOCATE (matrix_w_q)
3193 : CALL dbcsr_copy(matrix_w_q, matrix_s_aux_fit_vs_orb(1)%matrix, &
3194 880 : "W MATRIX AUX Q")
3195 880 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_w_q, sab_aux_fit_vs_orb)
3196 880 : ALLOCATE (matrix_w_s)
3197 : CALL dbcsr_create(matrix_w_s, template=matrix_s_aux_fit(1)%matrix, &
3198 : name='W MATRIX AUX S', &
3199 880 : matrix_type=dbcsr_type_no_symmetry)
3200 880 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_w_s, sab_aux_fit_asymm)
3201 :
3202 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
3203 880 : natom=natom, force=force)
3204 2640 : ALLOCATE (admm_force(3, natom))
3205 880 : admm_force = 0.0_dp
3206 :
3207 880 : nspins = SIZE(matrix_pz)
3208 880 : nao = admm_env%nao_orb
3209 880 : naux = admm_env%nao_aux_fit
3210 :
3211 880 : CALL cp_fm_set_all(admm_env%work_aux_orb2, 0.0_dp)
3212 :
3213 1860 : DO ispin = 1, nspins
3214 980 : CALL copy_dbcsr_to_fm(matrix_hz(ispin)%matrix, admm_env%work_aux_aux)
3215 : CALL parallel_gemm("N", "T", naux, naux, naux, 1.0_dp, admm_env%s_inv, &
3216 980 : admm_env%work_aux_aux, 0.0_dp, admm_env%work_aux_aux2)
3217 : CALL parallel_gemm("N", "N", naux, nao, naux, 1.0_dp, admm_env%work_aux_aux2, &
3218 980 : admm_env%A, 0.0_dp, admm_env%work_aux_orb)
3219 980 : CALL copy_dbcsr_to_fm(matrix_pz(ispin)%matrix, admm_env%work_orb_orb)
3220 : ! admm_env%work_aux_orb2 = S-1*H*A*P
3221 : CALL parallel_gemm("N", "N", naux, nao, nao, 1.0_dp, admm_env%work_aux_orb, &
3222 1860 : admm_env%work_orb_orb, 1.0_dp, admm_env%work_aux_orb2)
3223 : END DO
3224 :
3225 880 : CALL copy_fm_to_dbcsr(admm_env%work_aux_orb2, matrix_w_q, keep_sparsity=.TRUE.)
3226 :
3227 : ! admm_env%work_aux_aux = S-1*H*A*P*A(T)
3228 : CALL parallel_gemm("N", "T", naux, naux, nao, 1.0_dp, admm_env%work_aux_orb2, &
3229 880 : admm_env%A, 0.0_dp, admm_env%work_aux_aux)
3230 880 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_w_s, keep_sparsity=.TRUE.)
3231 :
3232 880 : CALL dbcsr_scale(matrix_w_q, -my_fval)
3233 880 : CALL dbcsr_scale(matrix_w_s, my_fval)
3234 :
3235 : CALL build_overlap_force(ks_env, admm_force, &
3236 : basis_type_a="AUX_FIT", basis_type_b="AUX_FIT", &
3237 880 : sab_nl=sab_aux_fit_asymm, matrix_p=matrix_w_s)
3238 : CALL build_overlap_force(ks_env, admm_force, &
3239 : basis_type_a="AUX_FIT", basis_type_b="ORB", &
3240 880 : sab_nl=sab_aux_fit_vs_orb, matrix_p=matrix_w_q)
3241 :
3242 : ! add forces
3243 880 : CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set)
3244 :
3245 880 : DEALLOCATE (admm_force)
3246 880 : CALL dbcsr_deallocate_matrix(matrix_w_s)
3247 880 : CALL dbcsr_deallocate_matrix(matrix_w_q)
3248 :
3249 880 : CALL timestop(handle)
3250 :
3251 880 : END SUBROUTINE admm_projection_derivative
3252 :
3253 : ! **************************************************************************************************
3254 : !> \brief Calculates contribution of forces due to basis transformation
3255 : !>
3256 : !> dE/dR = dE/dC'*dC'/dR
3257 : !> dE/dC = Ks'*c'*occ = H'
3258 : !>
3259 : !> dC'/dR = - tr(A*lambda^(-1/2)*H'^(T)*S^(-1) * dS'/dR)
3260 : !> - tr(A*C*Y^(T)*C^(T)*Q^(T)*A^(T) * dS'/dR)
3261 : !> + tr(C*lambda^(-1/2)*H'^(T)*S^(-1) * dQ/dR)
3262 : !> + tr(A*C*Y^(T)*c^(T) * dQ/dR)
3263 : !> + tr(C*Y^(T)*C^(T)*A^(T) * dQ/dR)
3264 : !>
3265 : !> where
3266 : !>
3267 : !> A = S'^(-1)*Q
3268 : !> lambda = C^(T)*B*C
3269 : !> B = Q^(T)*A
3270 : !> Y = R*[ (R^(T)*C^(T)*A^(T)*H'*R) xx M ]*R^(T)
3271 : !> lambda = R*D*R^(T)
3272 : !> Mij = Poles-Matrix (see above)
3273 : !> xx = schur product
3274 : !>
3275 : !> \param qs_env the QS environment
3276 : !> \par History
3277 : !> 05.2008 created [Manuel Guidon]
3278 : !> \author Manuel Guidon
3279 : ! **************************************************************************************************
3280 258 : SUBROUTINE calc_mixed_overlap_force(qs_env)
3281 :
3282 : TYPE(qs_environment_type), POINTER :: qs_env
3283 :
3284 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_mixed_overlap_force'
3285 :
3286 : INTEGER :: handle, ispin, iw, nao_aux_fit, nao_orb, &
3287 : natom, neighbor_list_id, nmo
3288 : LOGICAL :: omit_headers
3289 258 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: admm_force
3290 : TYPE(admm_type), POINTER :: admm_env
3291 258 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3292 : TYPE(cp_fm_type), POINTER :: mo_coeff
3293 : TYPE(cp_logger_type), POINTER :: logger
3294 258 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_aux_fit, &
3295 258 : matrix_s_aux_fit_vs_orb, rho_ao, &
3296 258 : rho_ao_aux
3297 : TYPE(dbcsr_type), POINTER :: matrix_rho_aux_desymm_tmp, matrix_w_q, &
3298 : matrix_w_s
3299 : TYPE(dft_control_type), POINTER :: dft_control
3300 258 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3301 : TYPE(mp_para_env_type), POINTER :: para_env
3302 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3303 258 : POINTER :: sab_orb
3304 : TYPE(qs_energy_type), POINTER :: energy
3305 258 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3306 : TYPE(qs_ks_env_type), POINTER :: ks_env
3307 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
3308 :
3309 258 : CALL timeset(routineN, handle)
3310 :
3311 258 : NULLIFY (admm_env, logger, dft_control, para_env, mos, mo_coeff, matrix_w_q, matrix_w_s, &
3312 258 : rho, rho_aux_fit, energy, sab_orb, ks_env, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, matrix_s)
3313 :
3314 : CALL get_qs_env(qs_env, &
3315 : admm_env=admm_env, &
3316 : ks_env=ks_env, &
3317 : dft_control=dft_control, &
3318 : matrix_s=matrix_s, &
3319 : neighbor_list_id=neighbor_list_id, &
3320 : rho=rho, &
3321 : energy=energy, &
3322 : sab_orb=sab_orb, &
3323 : mos=mos, &
3324 258 : para_env=para_env)
3325 : CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit, rho_aux_fit=rho_aux_fit, &
3326 258 : matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb)
3327 :
3328 258 : CALL qs_rho_get(rho, rho_ao=rho_ao)
3329 : CALL qs_rho_get(rho_aux_fit, &
3330 258 : rho_ao=rho_ao_aux)
3331 :
3332 258 : nao_aux_fit = admm_env%nao_aux_fit
3333 258 : nao_orb = admm_env%nao_orb
3334 :
3335 258 : logger => cp_get_default_logger()
3336 :
3337 : ! *** forces are only implemented for mo_diag or none and basis_projection ***
3338 258 : IF (admm_env%block_dm) THEN
3339 0 : CPABORT("ADMM Forces not implemented for blocked projection methods!")
3340 : END IF
3341 :
3342 258 : IF (.NOT. (admm_env%purification_method == do_admm_purify_mo_diag .OR. &
3343 : admm_env%purification_method == do_admm_purify_none)) THEN
3344 0 : CPABORT("ADMM Forces only implemented without purification or for MO_DIAG.")
3345 : END IF
3346 :
3347 : ! *** Create sparse work matrices
3348 :
3349 258 : ALLOCATE (matrix_w_s)
3350 : CALL dbcsr_create(matrix_w_s, template=matrix_s_aux_fit(1)%matrix, &
3351 : name='W MATRIX AUX S', &
3352 258 : matrix_type=dbcsr_type_no_symmetry)
3353 258 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_w_s, admm_env%sab_aux_fit_asymm)
3354 :
3355 258 : ALLOCATE (matrix_w_q)
3356 : CALL dbcsr_copy(matrix_w_q, matrix_s_aux_fit_vs_orb(1)%matrix, &
3357 258 : "W MATRIX AUX Q")
3358 :
3359 558 : DO ispin = 1, dft_control%nspins
3360 300 : nmo = admm_env%nmo(ispin)
3361 300 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
3362 :
3363 : ! *** S'^(-T)*H'
3364 300 : IF (.NOT. admm_env%purification_method == do_admm_purify_none) THEN
3365 : CALL parallel_gemm('T', 'N', nao_aux_fit, nmo, nao_aux_fit, &
3366 : 1.0_dp, admm_env%S_inv, admm_env%mo_derivs_aux_fit(ispin), 0.0_dp, &
3367 70 : admm_env%work_aux_nmo(ispin))
3368 : ELSE
3369 :
3370 : CALL parallel_gemm('T', 'N', nao_aux_fit, nmo, nao_aux_fit, &
3371 : 1.0_dp, admm_env%S_inv, admm_env%H(ispin), 0.0_dp, &
3372 230 : admm_env%work_aux_nmo(ispin))
3373 : END IF
3374 :
3375 : ! *** S'^(-T)*H'*Lambda^(-T/2)
3376 : CALL parallel_gemm('N', 'T', nao_aux_fit, nmo, nmo, &
3377 : 1.0_dp, admm_env%work_aux_nmo(ispin), admm_env%lambda_inv_sqrt(ispin), 0.0_dp, &
3378 300 : admm_env%work_aux_nmo2(ispin))
3379 :
3380 : ! *** C*Lambda^(-1/2)*H'^(T)*S'^(-1) minus sign due to force = -dE/dR
3381 : CALL parallel_gemm('N', 'T', nao_aux_fit, nao_orb, nmo, &
3382 : -1.0_dp, admm_env%work_aux_nmo2(ispin), mo_coeff, 0.0_dp, &
3383 300 : admm_env%work_aux_orb)
3384 :
3385 : ! *** A*C*Lambda^(-1/2)*H'^(T)*S'^(-1), minus sign to recover from above
3386 : CALL parallel_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nao_orb, &
3387 : -1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
3388 300 : admm_env%work_aux_aux)
3389 :
3390 300 : IF (.NOT. (admm_env%purification_method == do_admm_purify_none)) THEN
3391 : ! *** C*Y
3392 : CALL parallel_gemm('N', 'N', nao_orb, nmo, nmo, &
3393 : 1.0_dp, mo_coeff, admm_env%R_schur_R_t(ispin), 0.0_dp, &
3394 70 : admm_env%work_orb_nmo(ispin))
3395 : ! *** C*Y^(T)*C^(T)
3396 : CALL parallel_gemm('N', 'T', nao_orb, nao_orb, nmo, &
3397 : 1.0_dp, mo_coeff, admm_env%work_orb_nmo(ispin), 0.0_dp, &
3398 70 : admm_env%work_orb_orb)
3399 : ! *** A*C*Y^(T)*C^(T) Add to work aux_orb, minus sign due to force = -dE/dR
3400 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, &
3401 : -1.0_dp, admm_env%A, admm_env%work_orb_orb, 1.0_dp, &
3402 70 : admm_env%work_aux_orb)
3403 :
3404 : ! *** C*Y^(T)
3405 : CALL parallel_gemm('N', 'T', nao_orb, nmo, nmo, &
3406 : 1.0_dp, mo_coeff, admm_env%R_schur_R_t(ispin), 0.0_dp, &
3407 70 : admm_env%work_orb_nmo(ispin))
3408 : ! *** C*Y*C^(T)
3409 : CALL parallel_gemm('N', 'T', nao_orb, nao_orb, nmo, &
3410 : 1.0_dp, mo_coeff, admm_env%work_orb_nmo(ispin), 0.0_dp, &
3411 70 : admm_env%work_orb_orb)
3412 : ! *** A*C*Y*C^(T) Add to work aux_orb, minus sign due to -dE/dR
3413 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, &
3414 : -1.0_dp, admm_env%A, admm_env%work_orb_orb, 1.0_dp, &
3415 70 : admm_env%work_aux_orb)
3416 : END IF
3417 :
3418 : ! Add derivative contribution matrix*dQ/dR in additional last term in
3419 : ! Eq. (26,32, 33) in Merlot2014 to the force
3420 : ! ADMMS
3421 300 : IF (admm_env%do_admms) THEN
3422 : ! *** scale admm_env%work_aux_orb by gsi due to inner derivative
3423 12 : CALL cp_fm_scale(admm_env%gsi(ispin), admm_env%work_aux_orb)
3424 : CALL parallel_gemm('N', 'T', nao_orb, nao_orb, nmo, &
3425 : 4.0_dp*(admm_env%gsi(ispin))*admm_env%lambda_merlot(ispin)/dft_control%nspins, &
3426 12 : mo_coeff, mo_coeff, 0.0_dp, admm_env%work_orb_orb2)
3427 :
3428 : ! *** prefactor*A*C*C^(T) Add to work aux_orb
3429 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, &
3430 : 1.0_dp, admm_env%A, admm_env%work_orb_orb2, 1.0_dp, &
3431 12 : admm_env%work_aux_orb)
3432 :
3433 : ! ADMMP
3434 288 : ELSE IF (admm_env%do_admmp) THEN
3435 16 : CALL cp_fm_scale(admm_env%gsi(ispin)**2, admm_env%work_aux_orb)
3436 : ! *** prefactor*C*C^(T), nspins since 2/n_spin*C*C^(T)=P
3437 : CALL parallel_gemm('N', 'T', nao_orb, nao_orb, nmo, &
3438 : 4.0_dp*(admm_env%gsi(ispin))*admm_env%lambda_merlot(ispin)/dft_control%nspins, &
3439 16 : mo_coeff, mo_coeff, 0.0_dp, admm_env%work_orb_orb2)
3440 :
3441 : ! *** prefactor*A*C*C^(T) Add to work aux_orb
3442 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, &
3443 : 1.0_dp, admm_env%A, admm_env%work_orb_orb2, 1.0_dp, &
3444 16 : admm_env%work_aux_orb)
3445 :
3446 : ! ADMMQ
3447 272 : ELSE IF (admm_env%do_admmq) THEN
3448 : ! *** scale admm_env%work_aux_orb by gsi due to inner derivative
3449 12 : CALL cp_fm_scale(admm_env%gsi(ispin), admm_env%work_aux_orb)
3450 : CALL parallel_gemm('N', 'T', nao_orb, nao_orb, nmo, &
3451 : 4.0_dp*(admm_env%gsi(ispin))*admm_env%lambda_merlot(ispin)/dft_control%nspins, &
3452 12 : mo_coeff, mo_coeff, 0.0_dp, admm_env%work_orb_orb2)
3453 :
3454 : ! *** prefactor*A*C*C^(T) Add to work aux_orb
3455 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, &
3456 : 1.0_dp, admm_env%A, admm_env%work_orb_orb2, 1.0_dp, &
3457 12 : admm_env%work_aux_orb)
3458 : END IF
3459 :
3460 : ! *** copy to sparse matrix
3461 300 : CALL copy_fm_to_dbcsr(admm_env%work_aux_orb, matrix_w_q, keep_sparsity=.TRUE.)
3462 :
3463 300 : IF (.NOT. (admm_env%purification_method == do_admm_purify_none)) THEN
3464 : ! *** A*C*Y^(T)*C^(T)
3465 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, &
3466 : 1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
3467 70 : admm_env%work_aux_orb)
3468 : ! *** A*C*Y^(T)*C^(T)*A^(T) add to aux_aux, minus sign cancels
3469 : CALL parallel_gemm('N', 'T', nao_aux_fit, nao_aux_fit, nao_orb, &
3470 : 1.0_dp, admm_env%work_aux_orb, admm_env%A, 1.0_dp, &
3471 70 : admm_env%work_aux_aux)
3472 : END IF
3473 :
3474 : ! *** copy to sparse matrix
3475 300 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_w_s, keep_sparsity=.TRUE.)
3476 :
3477 : ! Add derivative of Eq. (33) with respect to s_aux Merlot2014 to the force
3478 300 : IF (admm_env%do_admmp .OR. admm_env%do_admmq .OR. admm_env%do_admms) THEN
3479 :
3480 : !Create desymmetrized auxiliary density matrix
3481 : NULLIFY (matrix_rho_aux_desymm_tmp)
3482 40 : ALLOCATE (matrix_rho_aux_desymm_tmp)
3483 : CALL dbcsr_create(matrix_rho_aux_desymm_tmp, template=matrix_s_aux_fit(1)%matrix, &
3484 : name='Rho_aux non-symm', &
3485 40 : matrix_type=dbcsr_type_no_symmetry)
3486 :
3487 40 : CALL dbcsr_desymmetrize(rho_ao_aux(ispin)%matrix, matrix_rho_aux_desymm_tmp)
3488 :
3489 : ! ADMMS/Q 1. scale original matrix_w_s by gsi due to inner deriv.
3490 : ! 2. add derivative of variational term with resp. to s
3491 40 : IF (admm_env%do_admms .OR. admm_env%do_admmq) THEN
3492 24 : CALL dbcsr_scale(matrix_w_s, admm_env%gsi(ispin))
3493 : CALL dbcsr_add(matrix_w_s, matrix_rho_aux_desymm_tmp, 1.0_dp, &
3494 24 : -admm_env%lambda_merlot(ispin))
3495 :
3496 : ! ADMMP scale by gsi^2 and add derivative of variational term with resp. to s
3497 16 : ELSE IF (admm_env%do_admmp) THEN
3498 :
3499 16 : CALL dbcsr_scale(matrix_w_s, admm_env%gsi(ispin)**2)
3500 : CALL dbcsr_add(matrix_w_s, matrix_rho_aux_desymm_tmp, 1.0_dp, &
3501 16 : (-admm_env%gsi(ispin))*admm_env%lambda_merlot(ispin))
3502 :
3503 : END IF
3504 :
3505 40 : CALL dbcsr_deallocate_matrix(matrix_rho_aux_desymm_tmp)
3506 :
3507 : END IF
3508 :
3509 : ! allocate force vector
3510 300 : CALL get_qs_env(qs_env=qs_env, natom=natom)
3511 900 : ALLOCATE (admm_force(3, natom))
3512 300 : admm_force = 0.0_dp
3513 : CALL build_overlap_force(ks_env, admm_force, &
3514 : basis_type_a="AUX_FIT", basis_type_b="AUX_FIT", &
3515 300 : sab_nl=admm_env%sab_aux_fit_asymm, matrix_p=matrix_w_s)
3516 : CALL build_overlap_force(ks_env, admm_force, &
3517 : basis_type_a="AUX_FIT", basis_type_b="ORB", &
3518 300 : sab_nl=admm_env%sab_aux_fit_vs_orb, matrix_p=matrix_w_q)
3519 :
3520 : ! Add contribution of original basis set for ADMMQ, P, S
3521 300 : IF (admm_env%do_admmq .OR. admm_env%do_admmp .OR. admm_env%do_admms) THEN
3522 40 : CALL dbcsr_scale(rho_ao(ispin)%matrix, -admm_env%lambda_merlot(ispin))
3523 : CALL build_overlap_force(ks_env, admm_force, &
3524 : basis_type_a="ORB", basis_type_b="ORB", &
3525 40 : sab_nl=sab_orb, matrix_p=rho_ao(ispin)%matrix)
3526 40 : CALL dbcsr_scale(rho_ao(ispin)%matrix, -1.0_dp/admm_env%lambda_merlot(ispin))
3527 : END IF
3528 :
3529 : ! add forces
3530 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
3531 300 : force=force)
3532 300 : CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set)
3533 300 : DEALLOCATE (admm_force)
3534 :
3535 300 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
3536 300 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
3537 : qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT"), cp_p_file)) THEN
3538 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT", &
3539 0 : extension=".Log")
3540 : CALL cp_dbcsr_write_sparse_matrix(matrix_w_s, 4, 6, qs_env, &
3541 0 : para_env, output_unit=iw, omit_headers=omit_headers)
3542 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
3543 0 : "DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT")
3544 : END IF
3545 300 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
3546 558 : qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT"), cp_p_file)) THEN
3547 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT", &
3548 0 : extension=".Log")
3549 : CALL cp_dbcsr_write_sparse_matrix(matrix_w_q, 4, 6, qs_env, &
3550 0 : para_env, output_unit=iw, omit_headers=omit_headers)
3551 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
3552 0 : "DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT")
3553 : END IF
3554 :
3555 : END DO !spin loop
3556 :
3557 : ! *** Deallocated weighted density matrices
3558 258 : CALL dbcsr_deallocate_matrix(matrix_w_s)
3559 258 : CALL dbcsr_deallocate_matrix(matrix_w_q)
3560 :
3561 258 : CALL timestop(handle)
3562 :
3563 516 : END SUBROUTINE calc_mixed_overlap_force
3564 :
3565 : ! **************************************************************************************************
3566 : !> \brief ...
3567 : !> \param admm_env environment of auxiliary DM
3568 : !> \param mo_set ...
3569 : !> \param density_matrix auxiliary DM
3570 : !> \param overlap_matrix auxiliary OM
3571 : !> \param density_matrix_large DM of the original basis
3572 : !> \param overlap_matrix_large overlap matrix of original basis
3573 : !> \param ispin ...
3574 : ! **************************************************************************************************
3575 15240 : SUBROUTINE calculate_dm_mo_no_diag(admm_env, mo_set, density_matrix, overlap_matrix, &
3576 : density_matrix_large, overlap_matrix_large, ispin)
3577 : TYPE(admm_type), POINTER :: admm_env
3578 : TYPE(mo_set_type), INTENT(IN) :: mo_set
3579 : TYPE(dbcsr_type), POINTER :: density_matrix, overlap_matrix, &
3580 : density_matrix_large, &
3581 : overlap_matrix_large
3582 : INTEGER :: ispin
3583 :
3584 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_dm_mo_no_diag'
3585 :
3586 : INTEGER :: handle, nao_aux_fit, nmo
3587 : REAL(KIND=dp) :: alpha, nel_tmp_aux
3588 :
3589 : ! Number of electrons in the aux. DM
3590 :
3591 15240 : CALL timeset(routineN, handle)
3592 :
3593 15240 : CALL dbcsr_set(density_matrix, 0.0_dp)
3594 15240 : nao_aux_fit = admm_env%nao_aux_fit
3595 15240 : nmo = admm_env%nmo(ispin)
3596 15240 : CALL cp_fm_to_fm(admm_env%C_hat(ispin), admm_env%work_aux_nmo(ispin))
3597 15240 : CALL cp_fm_column_scale(admm_env%work_aux_nmo(ispin), mo_set%occupation_numbers(1:mo_set%homo))
3598 :
3599 : CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nmo, &
3600 : 1.0_dp, admm_env%work_aux_nmo(ispin), admm_env%lambda_inv(ispin), 0.0_dp, &
3601 15240 : admm_env%work_aux_nmo2(ispin))
3602 :
3603 : ! The following IF doesn't do anything unless !alpha=mo_set%maxocc is uncommented.
3604 15240 : IF (.NOT. mo_set%uniform_occupation) THEN ! not all orbitals 1..homo are equally occupied
3605 360 : alpha = 1.0_dp
3606 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
3607 : matrix_v=admm_env%C_hat(ispin), &
3608 : matrix_g=admm_env%work_aux_nmo2(ispin), &
3609 : ncol=mo_set%homo, &
3610 360 : alpha=alpha)
3611 : ELSE
3612 14880 : alpha = 1.0_dp
3613 : !alpha=mo_set%maxocc
3614 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
3615 : matrix_v=admm_env%C_hat(ispin), &
3616 : matrix_g=admm_env%work_aux_nmo2(ispin), &
3617 : ncol=mo_set%homo, &
3618 14880 : alpha=alpha)
3619 : END IF
3620 :
3621 : ! The following IF checks whether gsi needs to be calculated. This is the case if
3622 : ! the auxiliary density matrix gets scaled
3623 : ! according to Eq. 22 (Merlot) or a scaling of exchange_correction is employed, Eq. 35 (Merlot).
3624 15240 : IF (admm_env%do_admmp .OR. admm_env%do_admmq .OR. admm_env%do_admms) THEN
3625 :
3626 1030 : CALL cite_reference(Merlot2014)
3627 :
3628 1030 : admm_env%n_large_basis(3) = 0.0_dp
3629 :
3630 : ! Calculate number of electrons in the original density matrix, transposing doesn't matter
3631 : ! since both matrices are symmetric
3632 1030 : CALL dbcsr_dot(density_matrix_large, overlap_matrix_large, admm_env%n_large_basis(ispin))
3633 1030 : admm_env%n_large_basis(3) = admm_env%n_large_basis(3) + admm_env%n_large_basis(ispin)
3634 : ! Calculate number of electrons in the auxiliary density matrix
3635 1030 : CALL dbcsr_dot(density_matrix, overlap_matrix, nel_tmp_aux)
3636 1030 : admm_env%gsi(ispin) = admm_env%n_large_basis(ispin)/nel_tmp_aux
3637 :
3638 1030 : IF (admm_env%do_admmq .OR. admm_env%do_admms) THEN
3639 : ! multiply aux. DM with gsi to get the scaled DM (Merlot, Eq. 21)
3640 600 : CALL dbcsr_scale(density_matrix, admm_env%gsi(ispin))
3641 : END IF
3642 :
3643 : END IF
3644 :
3645 15240 : CALL timestop(handle)
3646 :
3647 15240 : END SUBROUTINE calculate_dm_mo_no_diag
3648 :
3649 : ! **************************************************************************************************
3650 : !> \brief ...
3651 : !> \param admm_env ...
3652 : !> \param density_matrix ...
3653 : !> \param density_matrix_aux ...
3654 : !> \param ispin ...
3655 : !> \param nspins ...
3656 : ! **************************************************************************************************
3657 708 : SUBROUTINE blockify_density_matrix(admm_env, density_matrix, density_matrix_aux, &
3658 : ispin, nspins)
3659 : TYPE(admm_type), POINTER :: admm_env
3660 : TYPE(dbcsr_type), POINTER :: density_matrix, density_matrix_aux
3661 : INTEGER :: ispin, nspins
3662 :
3663 : CHARACTER(len=*), PARAMETER :: routineN = 'blockify_density_matrix'
3664 :
3665 : INTEGER :: handle, iatom, jatom
3666 : LOGICAL :: found
3667 354 : REAL(dp), DIMENSION(:, :), POINTER :: sparse_block, sparse_block_aux
3668 : TYPE(dbcsr_iterator_type) :: iter
3669 :
3670 354 : CALL timeset(routineN, handle)
3671 :
3672 : ! ** set blocked density matrix to 0
3673 354 : CALL dbcsr_set(density_matrix_aux, 0.0_dp)
3674 :
3675 : ! ** now loop through the list and copy corresponding blocks
3676 354 : CALL dbcsr_iterator_start(iter, density_matrix)
3677 1683 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3678 1329 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block)
3679 1683 : IF (admm_env%block_map(iatom, jatom) == 1) THEN
3680 : CALL dbcsr_get_block_p(density_matrix_aux, &
3681 924 : row=iatom, col=jatom, block=sparse_block_aux, found=found)
3682 924 : IF (found) THEN
3683 11016 : sparse_block_aux = sparse_block
3684 : END IF
3685 :
3686 : END IF
3687 : END DO
3688 354 : CALL dbcsr_iterator_stop(iter)
3689 :
3690 354 : CALL copy_dbcsr_to_fm(density_matrix_aux, admm_env%P_to_be_purified(ispin))
3691 354 : CALL cp_fm_uplo_to_full(admm_env%P_to_be_purified(ispin), admm_env%work_orb_orb2)
3692 :
3693 354 : IF (nspins == 1) THEN
3694 114 : CALL cp_fm_scale(0.5_dp, admm_env%P_to_be_purified(ispin))
3695 : END IF
3696 :
3697 354 : CALL timestop(handle)
3698 354 : END SUBROUTINE blockify_density_matrix
3699 :
3700 : ! **************************************************************************************************
3701 : !> \brief ...
3702 : !> \param x ...
3703 : !> \return ...
3704 : ! **************************************************************************************************
3705 2754 : ELEMENTAL FUNCTION delta(x)
3706 : REAL(KIND=dp), INTENT(IN) :: x
3707 : REAL(KIND=dp) :: delta
3708 :
3709 2754 : IF (x == 0.0_dp) THEN !TODO: exact comparison of reals?
3710 : delta = 1.0_dp
3711 : ELSE
3712 2754 : delta = 0.0_dp
3713 : END IF
3714 :
3715 2754 : END FUNCTION delta
3716 :
3717 : ! **************************************************************************************************
3718 : !> \brief ...
3719 : !> \param x ...
3720 : !> \return ...
3721 : ! **************************************************************************************************
3722 19180 : ELEMENTAL FUNCTION Heaviside(x)
3723 : REAL(KIND=dp), INTENT(IN) :: x
3724 : REAL(KIND=dp) :: Heaviside
3725 :
3726 19180 : IF (x < 0.0_dp) THEN
3727 : Heaviside = 0.0_dp
3728 : ELSE
3729 10404 : Heaviside = 1.0_dp
3730 : END IF
3731 19180 : END FUNCTION Heaviside
3732 :
3733 : ! **************************************************************************************************
3734 : !> \brief Calculate ADMM auxiliary response density
3735 : !> \param qs_env ...
3736 : !> \param dm ...
3737 : !> \param dm_admm ...
3738 : ! **************************************************************************************************
3739 2502 : SUBROUTINE admm_aux_response_density(qs_env, dm, dm_admm)
3740 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
3741 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: dm
3742 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: dm_admm
3743 :
3744 : CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_aux_response_density'
3745 :
3746 : INTEGER :: handle, ispin, nao, nao_aux, ncol, nspins
3747 : TYPE(admm_type), POINTER :: admm_env
3748 : TYPE(dft_control_type), POINTER :: dft_control
3749 :
3750 2502 : CALL timeset(routineN, handle)
3751 :
3752 2502 : CALL get_qs_env(qs_env, admm_env=admm_env, dft_control=dft_control)
3753 :
3754 2502 : nspins = dft_control%nspins
3755 :
3756 2502 : CPASSERT(ASSOCIATED(admm_env%A))
3757 2502 : CPASSERT(ASSOCIATED(admm_env%work_orb_orb))
3758 2502 : CPASSERT(ASSOCIATED(admm_env%work_aux_orb))
3759 2502 : CPASSERT(ASSOCIATED(admm_env%work_aux_aux))
3760 2502 : CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux, ncol_global=nao)
3761 :
3762 : ! P1 -> AUX BASIS
3763 2502 : CALL cp_fm_get_info(admm_env%work_orb_orb, nrow_global=nao, ncol_global=ncol)
3764 5298 : DO ispin = 1, nspins
3765 2796 : CALL copy_dbcsr_to_fm(dm(ispin)%matrix, admm_env%work_orb_orb)
3766 : CALL parallel_gemm('N', 'N', nao_aux, ncol, nao, 1.0_dp, admm_env%A, &
3767 2796 : admm_env%work_orb_orb, 0.0_dp, admm_env%work_aux_orb)
3768 : CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%A, &
3769 2796 : admm_env%work_aux_orb, 0.0_dp, admm_env%work_aux_aux)
3770 5298 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, dm_admm(ispin)%matrix, keep_sparsity=.TRUE.)
3771 : END DO
3772 :
3773 2502 : CALL timestop(handle)
3774 :
3775 2502 : END SUBROUTINE admm_aux_response_density
3776 :
3777 : ! **************************************************************************************************
3778 : !> \brief Store the ADMM auxiliary overlap and projection A = S^-1 Q in their native scalar type.
3779 : !> \param qs_env ...
3780 : !> \param calculate_forces ...
3781 : ! **************************************************************************************************
3782 52 : SUBROUTINE kpoint_calc_admm_matrices(qs_env, calculate_forces)
3783 : TYPE(qs_environment_type), POINTER :: qs_env
3784 : LOGICAL :: calculate_forces
3785 :
3786 : INTEGER :: igroup, ik, ikp, indx, kplocal, kpmax, &
3787 : nao_aux_fit, nao_orb, nkp, nkp_groups
3788 : INTEGER, DIMENSION(2) :: kp_range
3789 52 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
3790 52 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3791 : LOGICAL :: my_kpgrp, store_overlap, use_real_wfn
3792 52 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
3793 : TYPE(admm_type), POINTER :: admm_env
3794 52 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
3795 : TYPE(cp_cfm_type) :: cmat_aux_fit, cmat_aux_fit_vs_orb, &
3796 : cwork_aux_fit
3797 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct_aux_fit, &
3798 : matrix_struct_aux_fit_vs_orb
3799 : TYPE(cp_fm_type) :: fmdummy, imat_aux_fit, &
3800 : imat_aux_fit_vs_orb, rmat_aux_fit, &
3801 : rmat_aux_fit_vs_orb, work_aux_fit
3802 52 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fmwork
3803 52 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_aux_fit, matrix_s_aux_fit_vs_orb
3804 52 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: dbcsr_aux_fit, dbcsr_aux_fit_vs_orb
3805 : TYPE(kpoint_env_type), POINTER :: kp
3806 : TYPE(kpoint_type), POINTER :: kpoints
3807 : TYPE(mp_para_env_type), POINTER :: para_env_global, para_env_local
3808 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3809 52 : POINTER :: sab_aux_fit, sab_aux_fit_vs_orb
3810 :
3811 52 : NULLIFY (xkp, kp_dist, para_env_local, cell_to_index, admm_env, kp, &
3812 52 : kpoints, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, sab_aux_fit, sab_aux_fit_vs_orb, &
3813 52 : para_env_global, matrix_struct_aux_fit, matrix_struct_aux_fit_vs_orb)
3814 :
3815 52 : CALL get_qs_env(qs_env, kpoints=kpoints, admm_env=admm_env)
3816 :
3817 : CALL get_admm_env(admm_env, matrix_s_aux_fit_kp=matrix_s_aux_fit, &
3818 : matrix_s_aux_fit_vs_orb_kp=matrix_s_aux_fit_vs_orb, &
3819 : sab_aux_fit=sab_aux_fit, &
3820 52 : sab_aux_fit_vs_orb=sab_aux_fit_vs_orb)
3821 52 : store_overlap = admm_env%do_admmp .OR. admm_env%do_admmq .OR. admm_env%do_admms .OR. calculate_forces
3822 :
3823 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
3824 52 : nkp_groups=nkp_groups, kp_dist=kp_dist, cell_to_index=cell_to_index)
3825 52 : kplocal = kp_range(2) - kp_range(1) + 1
3826 154 : kpmax = MAXVAL(kp_dist(2, :) - kp_dist(1, :) + 1)
3827 :
3828 208 : ALLOCATE (dbcsr_aux_fit(3))
3829 52 : CALL dbcsr_create(dbcsr_aux_fit(1), template=matrix_s_aux_fit(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
3830 52 : CALL dbcsr_create(dbcsr_aux_fit(2), template=matrix_s_aux_fit(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
3831 52 : CALL dbcsr_create(dbcsr_aux_fit(3), template=matrix_s_aux_fit(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
3832 52 : CALL cp_dbcsr_alloc_block_from_nbl(dbcsr_aux_fit(1), sab_aux_fit)
3833 52 : CALL cp_dbcsr_alloc_block_from_nbl(dbcsr_aux_fit(2), sab_aux_fit)
3834 :
3835 156 : ALLOCATE (dbcsr_aux_fit_vs_orb(2))
3836 : CALL dbcsr_create(dbcsr_aux_fit_vs_orb(1), template=matrix_s_aux_fit_vs_orb(1, 1)%matrix, &
3837 52 : matrix_type=dbcsr_type_no_symmetry)
3838 : CALL dbcsr_create(dbcsr_aux_fit_vs_orb(2), template=matrix_s_aux_fit_vs_orb(1, 1)%matrix, &
3839 52 : matrix_type=dbcsr_type_no_symmetry)
3840 52 : CALL cp_dbcsr_alloc_block_from_nbl(dbcsr_aux_fit_vs_orb(1), sab_aux_fit_vs_orb)
3841 52 : CALL cp_dbcsr_alloc_block_from_nbl(dbcsr_aux_fit_vs_orb(2), sab_aux_fit_vs_orb)
3842 :
3843 : !Create global work fm
3844 52 : nao_aux_fit = admm_env%nao_aux_fit
3845 52 : nao_orb = admm_env%nao_orb
3846 52 : para_env_global => kpoints%blacs_env_all%para_env
3847 :
3848 260 : ALLOCATE (fmwork(4))
3849 : CALL cp_fm_struct_create(matrix_struct_aux_fit, context=kpoints%blacs_env_all, para_env=para_env_global, &
3850 52 : nrow_global=nao_aux_fit, ncol_global=nao_aux_fit)
3851 52 : CALL cp_fm_create(fmwork(1), matrix_struct_aux_fit)
3852 52 : IF (.NOT. use_real_wfn) CALL cp_fm_create(fmwork(2), matrix_struct_aux_fit)
3853 52 : CALL cp_fm_struct_release(matrix_struct_aux_fit)
3854 :
3855 : CALL cp_fm_struct_create(matrix_struct_aux_fit_vs_orb, context=kpoints%blacs_env_all, para_env=para_env_global, &
3856 52 : nrow_global=nao_aux_fit, ncol_global=nao_orb)
3857 52 : CALL cp_fm_create(fmwork(3), matrix_struct_aux_fit_vs_orb)
3858 52 : IF (.NOT. use_real_wfn) CALL cp_fm_create(fmwork(4), matrix_struct_aux_fit_vs_orb)
3859 52 : CALL cp_fm_struct_release(matrix_struct_aux_fit_vs_orb)
3860 :
3861 : !Create fm local to the KP groups
3862 52 : nao_aux_fit = admm_env%nao_aux_fit
3863 52 : nao_orb = admm_env%nao_orb
3864 52 : para_env_local => kpoints%blacs_env%para_env
3865 :
3866 : CALL cp_fm_struct_create(matrix_struct_aux_fit, context=kpoints%blacs_env, para_env=para_env_local, &
3867 52 : nrow_global=nao_aux_fit, ncol_global=nao_aux_fit)
3868 52 : CALL cp_fm_create(rmat_aux_fit, matrix_struct_aux_fit)
3869 52 : IF (use_real_wfn) THEN
3870 0 : CALL cp_fm_create(work_aux_fit, matrix_struct_aux_fit)
3871 : ELSE
3872 52 : CALL cp_fm_create(imat_aux_fit, matrix_struct_aux_fit)
3873 52 : CALL cp_cfm_create(cwork_aux_fit, matrix_struct_aux_fit)
3874 52 : CALL cp_cfm_create(cmat_aux_fit, matrix_struct_aux_fit)
3875 : END IF
3876 :
3877 : CALL cp_fm_struct_create(matrix_struct_aux_fit_vs_orb, context=kpoints%blacs_env, para_env=para_env_local, &
3878 52 : nrow_global=nao_aux_fit, ncol_global=nao_orb)
3879 52 : CALL cp_fm_create(rmat_aux_fit_vs_orb, matrix_struct_aux_fit_vs_orb)
3880 52 : IF (.NOT. use_real_wfn) THEN
3881 52 : CALL cp_fm_create(imat_aux_fit_vs_orb, matrix_struct_aux_fit_vs_orb)
3882 52 : CALL cp_cfm_create(cmat_aux_fit_vs_orb, matrix_struct_aux_fit_vs_orb)
3883 : END IF
3884 :
3885 3044 : ALLOCATE (info(nkp, 4))
3886 :
3887 : ! Steup and start all the communication
3888 52 : indx = 0
3889 344 : DO ikp = 1, kpmax
3890 926 : DO igroup = 1, nkp_groups
3891 582 : ik = kp_dist(1, igroup) + ikp - 1
3892 582 : IF (ik > kp_dist(2, igroup)) CYCLE
3893 566 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
3894 566 : indx = indx + 1
3895 :
3896 566 : IF (use_real_wfn) THEN
3897 : !AUX-AUX overlap
3898 0 : CALL dbcsr_set(dbcsr_aux_fit(1), 0.0_dp)
3899 : CALL rskp_transform(rmatrix=dbcsr_aux_fit(1), rsmat=matrix_s_aux_fit, ispin=1, &
3900 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
3901 0 : CALL dbcsr_desymmetrize(dbcsr_aux_fit(1), dbcsr_aux_fit(3))
3902 0 : CALL copy_dbcsr_to_fm(dbcsr_aux_fit(3), fmwork(1))
3903 :
3904 : !AUX-ORB overlap
3905 0 : CALL dbcsr_set(dbcsr_aux_fit_vs_orb(1), 0.0_dp)
3906 : CALL rskp_transform(rmatrix=dbcsr_aux_fit_vs_orb(1), rsmat=matrix_s_aux_fit_vs_orb, ispin=1, &
3907 0 : xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit_vs_orb)
3908 0 : CALL copy_dbcsr_to_fm(dbcsr_aux_fit_vs_orb(1), fmwork(3))
3909 : ELSE
3910 : !AUX-AUX overlap
3911 566 : CALL dbcsr_set(dbcsr_aux_fit(1), 0.0_dp)
3912 566 : CALL dbcsr_set(dbcsr_aux_fit(2), 0.0_dp)
3913 : CALL rskp_transform(rmatrix=dbcsr_aux_fit(1), cmatrix=dbcsr_aux_fit(2), rsmat=matrix_s_aux_fit, &
3914 566 : ispin=1, xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
3915 566 : CALL dbcsr_desymmetrize(dbcsr_aux_fit(1), dbcsr_aux_fit(3))
3916 566 : CALL copy_dbcsr_to_fm(dbcsr_aux_fit(3), fmwork(1))
3917 566 : CALL dbcsr_desymmetrize(dbcsr_aux_fit(2), dbcsr_aux_fit(3))
3918 566 : CALL copy_dbcsr_to_fm(dbcsr_aux_fit(3), fmwork(2))
3919 :
3920 : !AUX-ORB overlap
3921 566 : CALL dbcsr_set(dbcsr_aux_fit_vs_orb(1), 0.0_dp)
3922 566 : CALL dbcsr_set(dbcsr_aux_fit_vs_orb(2), 0.0_dp)
3923 : CALL rskp_transform(rmatrix=dbcsr_aux_fit_vs_orb(1), cmatrix=dbcsr_aux_fit_vs_orb(2), &
3924 : rsmat=matrix_s_aux_fit_vs_orb, ispin=1, xkp=xkp(1:3, ik), &
3925 566 : cell_to_index=cell_to_index, sab_nl=sab_aux_fit_vs_orb)
3926 566 : CALL copy_dbcsr_to_fm(dbcsr_aux_fit_vs_orb(1), fmwork(3))
3927 566 : CALL copy_dbcsr_to_fm(dbcsr_aux_fit_vs_orb(2), fmwork(4))
3928 : END IF
3929 :
3930 858 : IF (my_kpgrp) THEN
3931 284 : CALL cp_fm_start_copy_general(fmwork(1), rmat_aux_fit, para_env_global, info(indx, 1))
3932 284 : CALL cp_fm_start_copy_general(fmwork(3), rmat_aux_fit_vs_orb, para_env_global, info(indx, 3))
3933 284 : IF (.NOT. use_real_wfn) THEN
3934 284 : CALL cp_fm_start_copy_general(fmwork(2), imat_aux_fit, para_env_global, info(indx, 2))
3935 284 : CALL cp_fm_start_copy_general(fmwork(4), imat_aux_fit_vs_orb, para_env_global, info(indx, 4))
3936 : END IF
3937 : ELSE
3938 282 : CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env_global, info(indx, 1))
3939 282 : CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env_global, info(indx, 3))
3940 282 : IF (.NOT. use_real_wfn) THEN
3941 282 : CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env_global, info(indx, 2))
3942 282 : CALL cp_fm_start_copy_general(fmwork(4), fmdummy, para_env_global, info(indx, 4))
3943 : END IF
3944 : END IF
3945 :
3946 : END DO
3947 : END DO
3948 :
3949 : ! Finish communication and store
3950 : indx = 0
3951 344 : DO ikp = 1, kpmax
3952 874 : DO igroup = 1, nkp_groups
3953 582 : ik = kp_dist(1, igroup) + ikp - 1
3954 582 : IF (ik > kp_dist(2, igroup)) CYCLE
3955 566 : my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
3956 282 : indx = indx + 1
3957 :
3958 292 : IF (my_kpgrp) THEN
3959 284 : CALL cp_fm_finish_copy_general(rmat_aux_fit, info(indx, 1))
3960 284 : CALL cp_fm_finish_copy_general(rmat_aux_fit_vs_orb, info(indx, 3))
3961 284 : IF (.NOT. use_real_wfn) THEN
3962 284 : CALL cp_fm_finish_copy_general(imat_aux_fit, info(indx, 2))
3963 284 : CALL cp_fm_finish_copy_general(imat_aux_fit_vs_orb, info(indx, 4))
3964 : END IF
3965 : END IF
3966 : END DO
3967 :
3968 292 : IF (ikp > kplocal) CYCLE
3969 284 : kp => kpoints%kp_aux_env(ikp)%kpoint_env
3970 :
3971 : !Allocate local KP matrices
3972 284 : CALL cp_fm_release(kp%amat)
3973 284 : CALL cp_cfm_release(kp%camat)
3974 :
3975 : !Only need the overlap in case of ADMMP, ADMMQ or ADMMS, or for forces
3976 284 : IF (store_overlap) THEN
3977 284 : CALL cp_fm_release(kp%smat)
3978 284 : CALL cp_cfm_release(kp%csmat)
3979 : END IF
3980 :
3981 336 : IF (use_real_wfn) THEN
3982 0 : CALL cp_fm_create(kp%amat, matrix_struct_aux_fit_vs_orb)
3983 0 : IF (store_overlap) THEN
3984 0 : CALL cp_fm_create(kp%smat, matrix_struct_aux_fit)
3985 0 : CALL cp_fm_to_fm(rmat_aux_fit, kp%smat)
3986 : END IF
3987 :
3988 : !Invert S_aux
3989 0 : CALL cp_fm_cholesky_decompose(rmat_aux_fit)
3990 0 : CALL cp_fm_cholesky_invert(rmat_aux_fit)
3991 0 : CALL cp_fm_uplo_to_full(rmat_aux_fit, work_aux_fit)
3992 :
3993 : !A = S^-1 * Q
3994 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, 1.0_dp, &
3995 0 : rmat_aux_fit, rmat_aux_fit_vs_orb, 0.0_dp, kp%amat)
3996 : ELSE
3997 284 : CALL cp_cfm_create(kp%camat, matrix_struct_aux_fit_vs_orb)
3998 284 : CALL cp_fm_to_cfm(rmat_aux_fit, imat_aux_fit, cmat_aux_fit)
3999 284 : IF (store_overlap) THEN
4000 284 : CALL cp_cfm_create(kp%csmat, matrix_struct_aux_fit)
4001 284 : CALL cp_cfm_to_cfm(cmat_aux_fit, kp%csmat)
4002 : END IF
4003 :
4004 : !Invert S_aux
4005 284 : CALL cp_cfm_cholesky_decompose(cmat_aux_fit)
4006 284 : CALL cp_cfm_cholesky_invert(cmat_aux_fit)
4007 284 : CALL cp_cfm_uplo_to_full(cmat_aux_fit, cwork_aux_fit)
4008 :
4009 : !A = S^-1 * Q
4010 284 : CALL cp_fm_to_cfm(rmat_aux_fit_vs_orb, imat_aux_fit_vs_orb, cmat_aux_fit_vs_orb)
4011 : CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, z_one, &
4012 284 : cmat_aux_fit, cmat_aux_fit_vs_orb, z_zero, kp%camat)
4013 : END IF
4014 : END DO
4015 :
4016 : ! Clean up communication
4017 618 : DO indx = 1, SIZE(info, 1)
4018 566 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
4019 566 : CALL cp_fm_cleanup_copy_general(info(indx, 3))
4020 618 : IF (.NOT. use_real_wfn) THEN
4021 566 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
4022 566 : CALL cp_fm_cleanup_copy_general(info(indx, 4))
4023 : END IF
4024 : END DO
4025 :
4026 52 : CALL cp_fm_release(rmat_aux_fit)
4027 52 : CALL cp_fm_release(imat_aux_fit)
4028 52 : CALL cp_fm_release(work_aux_fit)
4029 52 : CALL cp_cfm_release(cwork_aux_fit)
4030 52 : CALL cp_cfm_release(cmat_aux_fit)
4031 52 : CALL cp_fm_release(rmat_aux_fit_vs_orb)
4032 52 : CALL cp_fm_release(imat_aux_fit_vs_orb)
4033 52 : CALL cp_cfm_release(cmat_aux_fit_vs_orb)
4034 52 : CALL cp_fm_struct_release(matrix_struct_aux_fit)
4035 52 : CALL cp_fm_struct_release(matrix_struct_aux_fit_vs_orb)
4036 :
4037 52 : CALL cp_fm_release(fmwork(1))
4038 52 : CALL cp_fm_release(fmwork(2))
4039 52 : CALL cp_fm_release(fmwork(3))
4040 52 : CALL cp_fm_release(fmwork(4))
4041 :
4042 52 : CALL dbcsr_release(dbcsr_aux_fit(1))
4043 52 : CALL dbcsr_release(dbcsr_aux_fit(2))
4044 52 : CALL dbcsr_release(dbcsr_aux_fit(3))
4045 52 : CALL dbcsr_release(dbcsr_aux_fit_vs_orb(1))
4046 52 : CALL dbcsr_release(dbcsr_aux_fit_vs_orb(2))
4047 :
4048 2524 : END SUBROUTINE kpoint_calc_admm_matrices
4049 :
4050 : END MODULE admm_methods
|
