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