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