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