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