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