(git:34ef472)
admm_methods.F
Go to the documentation of this file.
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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
36  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
37  copy_fm_to_dbcsr,&
38  cp_dbcsr_plus_fm_fm_t,&
39  dbcsr_allocate_matrix_set,&
40  dbcsr_deallocate_matrix_set
41  USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
42  USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
43  cp_fm_scale,&
44  cp_fm_scale_and_add,&
45  cp_fm_schur_product,&
46  cp_fm_upper_to_full
47  USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
48  cp_fm_cholesky_invert,&
49  cp_fm_cholesky_reduce,&
50  cp_fm_cholesky_restore
51  USE cp_fm_diag, ONLY: cp_fm_syevd
52  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
53  cp_fm_struct_release,&
54  cp_fm_struct_type
55  USE cp_fm_types, ONLY: &
56  copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, &
57  cp_fm_get_info, cp_fm_release, cp_fm_set_all, cp_fm_set_element, cp_fm_start_copy_general, &
58  cp_fm_to_fm, cp_fm_type
59  USE cp_log_handling, ONLY: cp_get_default_logger,&
60  cp_logger_type,&
61  cp_to_string
62  USE cp_output_handling, ONLY: cp_p_file,&
63  cp_print_key_finished_output,&
64  cp_print_key_should_output,&
65  cp_print_key_unit_nr
66  USE dbcsr_api, ONLY: &
67  dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
68  dbcsr_dot, dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
69  dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
70  dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
71  dbcsr_type_no_symmetry, dbcsr_type_symmetric
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  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  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r_aux
157  TYPE(admm_type), POINTER :: admm_env
158  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_aux_fit, &
159  matrix_s_aux_fit_vs_orb, rho_ao, &
160  rho_ao_aux
161  TYPE(dft_control_type), POINTER :: dft_control
162  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
163  TYPE(mp_para_env_type), POINTER :: para_env
164  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
165  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  CALL timeset(routinen, handle)
171 
172  NULLIFY (ks_env, admm_env, mos, mos_aux_fit, matrix_s_aux_fit, &
173  matrix_s_aux_fit_vs_orb, matrix_s, rho, rho_aux_fit, para_env)
174  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  rho=rho)
185  CALL get_admm_env(admm_env, mos_aux_fit=mos_aux_fit, matrix_s_aux_fit=matrix_s_aux_fit, &
186  matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, rho_aux_fit=rho_aux_fit)
187 
188  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  tot_rho_r=tot_rho_r_aux)
194 
195  gapw = admm_env%do_gapw
196 
197  ! convert mos from full to dbcsr matrices
198  DO ispin = 1, dft_control%nspins
199  IF (mos(ispin)%use_mo_coeff_b) THEN
200  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  mos, mos_aux_fit, s_mstruct_changed)
207 
208  DO ispin = 1, dft_control%nspins
209  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  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  ispin=ispin)
226 
227  END IF
228 
229  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  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  basis_type = "AUX_FIT"
239  task_list => admm_env%task_list_aux_fit
240  IF (gapw) THEN
241  basis_type = "AUX_FIT_SOFT"
242  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  task_list_external=task_list)
253 
254  END DO
255 
256  !If GAPW, also need to prepare the atomic densities
257  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  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  do_rho0=.false., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
266  END IF
267 
268  IF (dft_control%nspins == 1) THEN
269  admm_env%gsi(3) = admm_env%gsi(1)
270  ELSE
271  admm_env%gsi(3) = (admm_env%gsi(1) + admm_env%gsi(2))/2.0_dp
272  END IF
273 
274  CALL qs_rho_set(rho_aux_fit, rho_r_valid=.true., rho_g_valid=.true.)
275 
276  CALL timestop(handle)
277 
278  END SUBROUTINE admm_mo_calc_rho_aux
279 
280 ! **************************************************************************************************
281 !> \brief ...
282 !> \param qs_env ...
283 ! **************************************************************************************************
284  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  INTEGER, DIMENSION(:, :), POINTER :: kp_dist
295  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  REAL(kind=dp), DIMENSION(:), POINTER :: occ_num, occ_num_aux, tot_rho_r_aux
301  REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
302  TYPE(admm_type), POINTER :: admm_env
303  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
313  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_s_aux_fit, rho_ao_aux, &
314  rho_ao_orb
315  TYPE(dbcsr_type) :: pmatrix_tmp
316  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  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
321  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  POINTER :: sab_aux_fit, sab_kp
325  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
326  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  CALL timeset(routinen, handle)
333 
334  NULLIFY (ks_env, admm_env, mos, mos_aux_fit, matrix_s, rho_orb, &
335  matrix_s_aux_fit, rho_aux_fit, rho_ao_orb, &
336  para_env, rho_g_aux, rho_r_aux, rho_ao_aux, tot_rho_r_aux, &
337  kpoints, sab_aux_fit, sab_kp, kp, &
338  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  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  sab_aux_fit=sab_aux_fit)
353  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  tot_rho_r=tot_rho_r_aux)
360 
361  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  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  ALLOCATE (pmatrix(2))
369  CALL dbcsr_create(pmatrix(1), template=matrix_s(1, 1)%matrix, &
370  matrix_type=dbcsr_type_symmetric)
371  CALL dbcsr_create(pmatrix(2), template=matrix_s(1, 1)%matrix, &
372  matrix_type=dbcsr_type_antisymmetric)
373  CALL dbcsr_create(pmatrix_tmp, template=matrix_s(1, 1)%matrix, &
374  matrix_type=dbcsr_type_no_symmetry)
375  CALL cp_dbcsr_alloc_block_from_nbl(pmatrix(1), sab_kp)
376  CALL cp_dbcsr_alloc_block_from_nbl(pmatrix(2), sab_kp)
377 
378  nao_aux_fit = admm_env%nao_aux_fit
379  nao_orb = admm_env%nao_orb
380  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  nrow_global=nao_orb, ncol_global=nao_orb)
385  CALL cp_fm_create(work_orb_orb, struct_orb_orb)
386  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  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  nrow_global=nao_aux_fit, ncol_global=nao_orb)
393  CALL cp_fm_create(work_aux_orb, struct_orb_orb)
394 
395  IF (.NOT. use_real_wfn) THEN
396  CALL cp_cfm_create(cpmatrix, struct_orb_orb)
397 
398  CALL cp_cfm_create(cwork_aux_aux, struct_aux_aux)
399 
400  CALL cp_cfm_create(ca, struct_aux_orb)
401  CALL cp_cfm_create(cwork_aux_orb, struct_aux_orb)
402 
403  CALL get_kpoint_env(kpoints%kp_env(1)%kpoint_env, mos=mos_kp)
404  mos => mos_kp(1, :)
405  CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
406  CALL cp_fm_get_info(mo_coeff, matrix_struct=mo_struct)
407  CALL cp_cfm_create(cmo_coeff, mo_struct)
408 
409  CALL get_kpoint_env(kpoints%kp_aux_env(1)%kpoint_env, mos=mos_aux_fit_kp)
410  mos => mos_aux_fit_kp(1, :)
411  CALL get_mo_set(mos(1), mo_coeff=mo_coeff_aux_fit)
412  CALL cp_fm_get_info(mo_coeff_aux_fit, matrix_struct=mo_struct_aux_fit)
413  CALL cp_cfm_create(cmo_coeff_aux_fit, mo_struct_aux_fit)
414  END IF
415 
416  CALL cp_fm_struct_release(struct_orb_orb)
417  CALL cp_fm_struct_release(struct_aux_aux)
418  CALL cp_fm_struct_release(struct_aux_orb)
419 
420  para_env => kpoints%blacs_env_all%para_env
421  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  maxval_mos = 0.0_dp
426  indx = 0
427  DO ikp = 1, kplocal
428  DO ispin = 1, nspins
429  DO igroup = 1, nkp_groups
430  ! number of current kpoint
431  ik = kp_dist(1, igroup) + ikp - 1
432  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
433  indx = indx + 1
434 
435  CALL get_kpoint_env(kpoints%kp_env(ikp)%kpoint_env, mos=mos_kp)
436  mos => mos_kp(1, :)
437  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
438  maxval_mos = max(maxval_mos, maxval(abs(mo_coeff%local_data)))
439 
440  IF (.NOT. use_real_wfn) THEN
441  mos => mos_kp(2, :)
442  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
443  maxval_mos = max(maxval_mos, maxval(abs(mo_coeff%local_data)))
444  END IF
445  END DO
446  END DO
447  END DO
448  CALL para_env%sum(maxval_mos) !I think para_env is the global one
449 
450  pmat_from_rs = .false.
451  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  ALLOCATE (info(kplocal*nspins*nkp_groups, 2))
456  !Start communication: only P matrix, and only if required
457  indx = 0
458  IF (pmat_from_rs) THEN
459  DO ikp = 1, kplocal
460  DO ispin = 1, nspins
461  DO igroup = 1, nkp_groups
462  ! number of current kpoint
463  ik = kp_dist(1, igroup) + ikp - 1
464  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
465  indx = indx + 1
466 
467  ! FT of matrices P if required, then transfer to FM type
468  IF (use_real_wfn) THEN
469  CALL dbcsr_set(pmatrix(1), 0.0_dp)
470  CALL rskp_transform(rmatrix=pmatrix(1), rsmat=rho_ao_orb, ispin=ispin, &
471  xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_kp)
472  CALL dbcsr_desymmetrize(pmatrix(1), pmatrix_tmp)
473  CALL copy_dbcsr_to_fm(pmatrix_tmp, admm_env%work_orb_orb)
474  ELSE
475  CALL dbcsr_set(pmatrix(1), 0.0_dp)
476  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  xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_kp)
479  CALL dbcsr_desymmetrize(pmatrix(1), pmatrix_tmp)
480  CALL copy_dbcsr_to_fm(pmatrix_tmp, admm_env%work_orb_orb)
481  CALL dbcsr_desymmetrize(pmatrix(2), pmatrix_tmp)
482  CALL copy_dbcsr_to_fm(pmatrix_tmp, admm_env%work_orb_orb2)
483  END IF
484 
485  IF (my_kpgrp) THEN
486  CALL cp_fm_start_copy_general(admm_env%work_orb_orb, work_orb_orb, para_env, info(indx, 1))
487  IF (.NOT. use_real_wfn) THEN
488  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  CALL cp_fm_start_copy_general(admm_env%work_orb_orb, fmdummy, para_env, info(indx, 1))
492  IF (.NOT. use_real_wfn) THEN
493  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  DO ikp = 1, kplocal
503  DO ispin = 1, nspins
504  DO igroup = 1, nkp_groups
505  ! number of current kpoint
506  ik = kp_dist(1, igroup) + ikp - 1
507  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
508  indx = indx + 1
509  IF (my_kpgrp .AND. pmat_from_rs) THEN
510  CALL cp_fm_finish_copy_general(work_orb_orb, info(indx, 1))
511  IF (.NOT. use_real_wfn) THEN
512  CALL cp_fm_finish_copy_general(work_orb_orb2, info(indx, 2))
513  CALL cp_fm_to_cfm(work_orb_orb, work_orb_orb2, cpmatrix)
514  END IF
515  END IF
516  END DO
517 
518  IF (use_real_wfn) THEN
519 
520  nmo = admm_env%nmo(ispin)
521  !! Each kpoint group has now information on a kpoint for which to calculate the MOS_aux
522  CALL get_kpoint_env(kpoints%kp_env(ikp)%kpoint_env, mos=mos_kp)
523  CALL get_kpoint_env(kpoints%kp_aux_env(ikp)%kpoint_env, mos=mos_aux_fit_kp)
524  mos => mos_kp(1, :)
525  mos_aux_fit => mos_aux_fit_kp(1, :)
526 
527  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  occupation_numbers=occ_num_aux)
530 
531  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  mo_coeff, 0.0_dp, mo_coeff_aux_fit)
534 
535  occ_num_aux(1:nmo) = occ_num(1:nmo)
536 
537  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  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  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  nmo = admm_env%nmo(ispin)
549  CALL get_kpoint_env(kpoints%kp_env(ikp)%kpoint_env, mos=mos_kp)
550  mos => mos_kp(1, :) !real
551  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
552  CALL cp_cfm_scale_and_add_fm(z_zero, cmo_coeff, z_one, mo_coeff)
553  mos => mos_kp(2, :) !complex
554  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
555  CALL cp_cfm_scale_and_add_fm(z_one, cmo_coeff, gaussi, mo_coeff)
556 
557  !project
558  kp => kpoints%kp_aux_env(ikp)%kpoint_env
559  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  z_one, ca, cmo_coeff, z_zero, cmo_coeff_aux_fit)
562 
563  !write result back to KP MOs
564  CALL get_kpoint_env(kpoints%kp_aux_env(ikp)%kpoint_env, mos=mos_aux_fit_kp)
565  mos_aux_fit => mos_aux_fit_kp(1, :)
566  CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
567  CALL cp_cfm_to_fm(cmo_coeff_aux_fit, mtargetr=mo_coeff_aux_fit)
568  mos_aux_fit => mos_aux_fit_kp(2, :)
569  CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
570  CALL cp_cfm_to_fm(cmo_coeff_aux_fit, mtargeti=mo_coeff_aux_fit)
571 
572  DO i = 1, 2
573  mos => mos_kp(i, :)
574  CALL get_mo_set(mos(ispin), occupation_numbers=occ_num)
575  mos_aux_fit => mos_aux_fit_kp(i, :)
576  CALL get_mo_set(mos_aux_fit(ispin), occupation_numbers=occ_num_aux)
577  occ_num_aux(:) = occ_num(:)
578  END DO
579 
580  IF (pmat_from_rs) THEN
581  CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_orb, z_one, ca, &
582  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  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  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  IF (pmat_from_rs) THEN
596  indx = 0
597  DO ikp = 1, kplocal
598  DO ispin = 1, nspins
599  DO igroup = 1, nkp_groups
600  ! number of current kpoint
601  ik = kp_dist(1, igroup) + ikp - 1
602  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
603  indx = indx + 1
604 
605  CALL cp_fm_cleanup_copy_general(info(indx, 1))
606  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  DEALLOCATE (info)
613  CALL dbcsr_release(pmatrix(1))
614  CALL dbcsr_release(pmatrix(2))
615  CALL dbcsr_release(pmatrix_tmp)
616 
617  CALL cp_fm_release(work_orb_orb)
618  CALL cp_fm_release(work_orb_orb2)
619  CALL cp_fm_release(work_aux_orb)
620  IF (.NOT. use_real_wfn) THEN
621  CALL cp_cfm_release(cpmatrix)
622  CALL cp_cfm_release(cwork_aux_aux)
623  CALL cp_cfm_release(cwork_aux_orb)
624  CALL cp_cfm_release(ca)
625  CALL cp_cfm_release(cmo_coeff)
626  CALL cp_cfm_release(cmo_coeff_aux_fit)
627  END IF
628 
629  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  admm_env%scf_work_aux_fit, for_aux_fit=.true.)
633 
634  !ADMMQ, ADMMP, ADMMS
635  IF (admm_env%do_admmq .OR. admm_env%do_admmp .OR. admm_env%do_admms) THEN
636 
637  CALL cite_reference(merlot2014)
638 
639  nelec_orb = 0.0_dp
640  nelec_aux = 0.0_dp
641  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  DO img = 1, dft_control%nimages
645  DO ispin = 1, dft_control%nspins
646  CALL dbcsr_dot(rho_ao_orb(ispin, img)%matrix, matrix_s(1, img)%matrix, tmp)
647  nelec_orb(ispin) = nelec_orb(ispin) + tmp
648  CALL dbcsr_dot(rho_ao_aux(ispin, img)%matrix, matrix_s_aux_fit(1, img)%matrix, tmp)
649  nelec_aux(ispin) = nelec_aux(ispin) + tmp
650  END DO
651  END DO
652 
653  DO ispin = 1, dft_control%nspins
654  admm_env%n_large_basis(ispin) = nelec_orb(ispin)
655  admm_env%gsi(ispin) = nelec_orb(ispin)/nelec_aux(ispin)
656  END DO
657 
658  IF (admm_env%charge_constrain) THEN
659  DO img = 1, dft_control%nimages
660  DO ispin = 1, dft_control%nspins
661  CALL dbcsr_scale(rho_ao_aux(ispin, img)%matrix, admm_env%gsi(ispin))
662  END DO
663  END DO
664  END IF
665 
666  IF (dft_control%nspins == 1) THEN
667  admm_env%gsi(3) = admm_env%gsi(1)
668  ELSE
669  admm_env%gsi(3) = (admm_env%gsi(1) + admm_env%gsi(2))/2.0_dp
670  END IF
671  END IF
672 
673  basis_type = "AUX_FIT"
674  task_list => admm_env%task_list_aux_fit
675  IF (gapw) THEN
676  basis_type = "AUX_FIT_SOFT"
677  task_list => admm_env%admm_gapw_env%task_list
678  END IF
679 
680  DO ispin = 1, nspins
681  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  task_list_external=task_list)
690  END DO
691 
692  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  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  do_rho0=.false., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
701  END IF
702 
703  CALL qs_rho_set(rho_aux_fit, rho_r_valid=.true., rho_g_valid=.true.)
704 
705  CALL timestop(handle)
706 
707  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  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  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit, &
725  matrix_ks_aux_fit_dft, &
726  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  NULLIFY (matrix_ks_aux_fit, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, rho_ao_aux, rho_aux_fit)
731  NULLIFY (admm_env, dft_control)
732 
733  CALL timeset(routinen, handle)
734 
735  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  matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx)
739  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  force_fac = 1.0_dp
743  IF (admm_env%do_admms) THEN
744  DO ispin = 1, dft_control%nspins
745  force_fac(ispin) = admm_env%gsi(ispin)**(2.0_dp/3.0_dp)
746  END DO
747  ELSE IF (admm_env%do_admmp) THEN
748  DO ispin = 1, dft_control%nspins
749  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  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  DO img = 1, dft_control%nimages
761  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  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  1.0_dp, 1.0_dp)
766  END DO
767  END DO
768 
769  CALL timestop(handle)
770 
771  END SUBROUTINE admm_update_ks_atom
772 
773 ! **************************************************************************************************
774 !> \brief ...
775 !> \param qs_env ...
776 ! **************************************************************************************************
777  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  CALL timeset(routinen, handle)
787  NULLIFY (admm_env)
788 
789  CALL get_qs_env(qs_env, admm_env=admm_env, dft_control=dft_control)
790 
791  SELECT CASE (admm_env%purification_method)
792  CASE (do_admm_purify_cauchy)
793  CALL merge_ks_matrix_cauchy(qs_env)
794 
795  CASE (do_admm_purify_cauchy_subspace)
796  CALL merge_ks_matrix_cauchy_subspace(qs_env)
797 
798  CASE (do_admm_purify_none)
799  IF (dft_control%nimages > 1) THEN
800  CALL merge_ks_matrix_none_kp(qs_env)
801  ELSE
802  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  cpabort("admm_mo_merge_ks_matrix: unknown purification method")
809  END SELECT
810 
811  CALL timestop(handle)
812 
813  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  SUBROUTINE admm_mo_merge_derivs(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, mo_derivs, &
827  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  CALL timeset(routinen, handle)
840 
841  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  mo_derivs, mo_derivs_aux_fit, matrix_ks_aux_fit)
845 
846  CASE (do_admm_purify_mo_no_diag)
847  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  cpabort("admm_mo_merge_derivs: unknown purification method")
853  END SELECT
854 
855  CALL timestop(handle)
856 
857  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  SUBROUTINE admm_fit_mo_coeffs(admm_env, matrix_s_aux_fit, matrix_s_mixed, &
869  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  CALL timeset(routinen, handle)
881 
882  IF (geometry_did_change) THEN
883  CALL fit_mo_coeffs(admm_env, matrix_s_aux_fit, matrix_s_mixed)
884  END IF
885 
886  SELECT CASE (admm_env%purification_method)
887  CASE (do_admm_purify_mo_no_diag, do_admm_purify_cauchy_subspace)
888  CALL purify_mo_cholesky(admm_env, mos, mos_aux_fit)
889 
890  CASE (do_admm_purify_mo_diag)
891  CALL purify_mo_diag(admm_env, mos, mos_aux_fit)
892 
893  CASE DEFAULT
894  CALL purify_mo_none(admm_env, mos, mos_aux_fit)
895  END SELECT
896 
897  CALL timestop(handle)
898 
899  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  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  REAL(dp), DIMENSION(:, :), POINTER :: sparse_block
916  TYPE(dbcsr_iterator_type) :: iter
917  TYPE(dbcsr_type), POINTER :: matrix_s_tilde
918 
919  CALL timeset(routinen, handle)
920 
921  nao_aux_fit = admm_env%nao_aux_fit
922  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  IF (.NOT. admm_env%block_fit) THEN
927  CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix, admm_env%S_inv)
928  ELSE
929  NULLIFY (matrix_s_tilde)
930  ALLOCATE (matrix_s_tilde)
931  CALL dbcsr_create(matrix_s_tilde, template=matrix_s_aux_fit(1)%matrix, &
932  name='MATRIX s_tilde', &
933  matrix_type=dbcsr_type_symmetric)
934 
935  CALL dbcsr_copy(matrix_s_tilde, matrix_s_aux_fit(1)%matrix)
936 
937  CALL dbcsr_iterator_start(iter, matrix_s_tilde)
938  DO WHILE (dbcsr_iterator_blocks_left(iter))
939  CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk)
940  IF (admm_env%block_map(iatom, jatom) == 0) THEN
941  sparse_block = 0.0_dp
942  END IF
943  END DO
944  CALL dbcsr_iterator_stop(iter)
945  CALL copy_dbcsr_to_fm(matrix_s_tilde, admm_env%S_inv)
946  CALL dbcsr_deallocate_matrix(matrix_s_tilde)
947  END IF
948 
949  CALL cp_fm_upper_to_full(admm_env%S_inv, admm_env%work_aux_aux)
950  CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S)
951 
952  CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix, admm_env%Q)
953 
954  !! Calculate S'_inverse
955  CALL cp_fm_cholesky_decompose(admm_env%S_inv)
956  CALL cp_fm_cholesky_invert(admm_env%S_inv)
957  !! Symmetrize the guy
958  CALL cp_fm_upper_to_full(admm_env%S_inv, admm_env%work_aux_aux)
959 
960  !! Calculate A=S'^(-1)*Q
961  IF (admm_env%block_fit) THEN
962  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  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  admm_env%B)
973  END IF
974 
975  CALL timestop(handle)
976 
977  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  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  CALL timeset(routinen, handle)
1002 
1003  nao_aux_fit = admm_env%nao_aux_fit
1004  nao_orb = admm_env%nao_orb
1005  nspins = SIZE(mos)
1006 
1007  ! *** Calculate the mo_coeffs for the fitting basis
1008  DO ispin = 1, nspins
1009  nmo = admm_env%nmo(ispin)
1010  IF (nmo == 0) cycle
1011  !! Lambda = C^(T)*B*C
1012  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
1013  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  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  admm_env%lambda(ispin))
1020  CALL cp_fm_to_fm(admm_env%lambda(ispin), admm_env%work_nmo_nmo1(ispin))
1021 
1022  CALL cp_fm_cholesky_decompose(admm_env%work_nmo_nmo1(ispin))
1023  CALL cp_fm_cholesky_invert(admm_env%work_nmo_nmo1(ispin))
1024  !! Symmetrize the guy
1025  CALL cp_fm_upper_to_full(admm_env%work_nmo_nmo1(ispin), admm_env%lambda_inv(ispin))
1026  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  admm_env%C_hat(ispin))
1032  CALL cp_fm_to_fm(admm_env%C_hat(ispin), mo_coeff_aux_fit)
1033 
1034  END DO
1035 
1036  CALL timestop(handle)
1037 
1038  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  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  REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig_work
1061  TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
1062 
1063  CALL timeset(routinen, handle)
1064 
1065  nao_aux_fit = admm_env%nao_aux_fit
1066  nao_orb = admm_env%nao_orb
1067  nspins = SIZE(mos)
1068 
1069  ! *** Calculate the mo_coeffs for the fitting basis
1070  DO ispin = 1, nspins
1071  nmo = admm_env%nmo(ispin)
1072  IF (nmo == 0) cycle
1073  !! Lambda = C^(T)*B*C
1074  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
1075  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  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  admm_env%lambda(ispin))
1082  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  admm_env%eigvals_lambda(ispin)%eigvals%data)
1086  ALLOCATE (eig_work(nmo))
1087  DO i = 1, nmo
1088  eig_work(i) = 1.0_dp/sqrt(admm_env%eigvals_lambda(ispin)%eigvals%data(i))
1089  END DO
1090  CALL cp_fm_to_fm(admm_env%R(ispin), admm_env%work_nmo_nmo1(ispin))
1091  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  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  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  mo_coeff_aux_fit)
1101 
1102  CALL cp_fm_to_fm(mo_coeff_aux_fit, admm_env%C_hat(ispin))
1103  CALL cp_fm_set_all(admm_env%lambda_inv(ispin), 0.0_dp, 1.0_dp)
1104  DEALLOCATE (eig_work)
1105  END DO
1106 
1107  CALL timestop(handle)
1108 
1109  END SUBROUTINE purify_mo_diag
1110 
1111 ! **************************************************************************************************
1112 !> \brief ...
1113 !> \param admm_env ...
1114 !> \param mos ...
1115 !> \param mos_aux_fit ...
1116 ! **************************************************************************************************
1117  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  REAL(kind=dp), DIMENSION(:), POINTER :: occ_num, occ_num_aux
1126  TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
1127 
1128  CALL timeset(routinen, handle)
1129 
1130  nao_aux_fit = admm_env%nao_aux_fit
1131  nao_orb = admm_env%nao_orb
1132  nspins = SIZE(mos)
1133 
1134  DO ispin = 1, nspins
1135  nmo = admm_env%nmo(ispin)
1136  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  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  mo_coeff_aux_fit)
1143  CALL cp_fm_to_fm(mo_coeff_aux_fit, admm_env%C_hat(ispin))
1144 
1145  occ_num_aux(1:nmo) = occ_num(1:nmo)
1146  ! XXXX should only be done first time XXXX
1147  CALL cp_fm_set_all(admm_env%lambda(ispin), 0.0_dp, 1.0_dp)
1148  CALL cp_fm_set_all(admm_env%lambda_inv(ispin), 0.0_dp, 1.0_dp)
1149  CALL cp_fm_set_all(admm_env%lambda_inv_sqrt(ispin), 0.0_dp, 1.0_dp)
1150  END DO
1151 
1152  CALL timestop(handle)
1153 
1154  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  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  CALL timeset(routinen, handle)
1180 
1181  nao_aux_fit = admm_env%nao_aux_fit
1182  nao_orb = admm_env%nao_orb
1183  nmo = admm_env%nmo(ispin)
1184 
1185  nspins = SIZE(admm_env%P_to_be_purified)
1186 
1187  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  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  admm_env%P_to_be_purified(ispin))
1196  END IF
1197 
1198  CALL cp_fm_to_fm(admm_env%S, admm_env%work_aux_aux)
1199  CALL cp_fm_to_fm(admm_env%P_to_be_purified(ispin), admm_env%work_aux_aux2)
1200 
1201  CALL cp_fm_cholesky_decompose(admm_env%work_aux_aux)
1202 
1203  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  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  admm_env%work_aux_aux3, op="MULTIPLY", pos="LEFT", transa="T")
1210 
1211  CALL cp_fm_to_fm(admm_env%work_aux_aux3, admm_env%R_purify(ispin))
1212 
1213  ! *** Construct Matrix M for Hadamard Product
1214  CALL cp_fm_set_all(admm_env%M_purify(ispin), 0.0_dp)
1215  pole = 0.0_dp
1216  DO i = 1, nao_aux_fit
1217  pole = heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i) - 0.5_dp)
1218  CALL cp_fm_set_element(admm_env%M_purify(ispin), i, i, pole)
1219  END DO
1220  CALL cp_fm_upper_to_full(admm_env%M_purify(ispin), admm_env%work_aux_aux)
1221 
1222  CALL copy_dbcsr_to_fm(density_matrix, admm_env%work_aux_aux3)
1223  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  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  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  admm_env%work_aux_aux3)
1237 
1238  CALL copy_fm_to_dbcsr(admm_env%work_aux_aux3, density_matrix, keep_sparsity=.true.)
1239 
1240  IF (nspins == 1) THEN
1241  CALL dbcsr_scale(density_matrix, 2.0_dp)
1242  END IF
1243 
1244  CALL timestop(handle)
1245 
1246  END SUBROUTINE purify_dm_cauchy
1247 
1248 ! **************************************************************************************************
1249 !> \brief ...
1250 !> \param qs_env ...
1251 ! **************************************************************************************************
1252  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  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  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  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1268 
1269  CALL timeset(routinen, handle)
1270  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  mos=mos)
1277  CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit)
1278 
1279  DO ispin = 1, dft_control%nspins
1280  nao_aux_fit = admm_env%nao_aux_fit
1281  nao_orb = admm_env%nao_orb
1282  nmo = admm_env%nmo(ispin)
1283  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1284 
1285  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  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  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  admm_env%P_to_be_purified(ispin))
1299 
1300  END IF
1301 
1302  CALL cp_fm_to_fm(admm_env%S, admm_env%work_aux_aux)
1303  CALL cp_fm_to_fm(admm_env%P_to_be_purified(ispin), admm_env%work_aux_aux2)
1304 
1305  CALL cp_fm_cholesky_decompose(admm_env%work_aux_aux)
1306 
1307  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  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  admm_env%work_aux_aux3, op="MULTIPLY", pos="LEFT", transa="T")
1314 
1315  CALL cp_fm_to_fm(admm_env%work_aux_aux3, admm_env%R_purify(ispin))
1316 
1317  ! *** Construct Matrix M for Hadamard Product
1318  pole = 0.0_dp
1319  DO i = 1, nao_aux_fit
1320  DO j = i, nao_aux_fit
1321  eig_diff = (admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i) - &
1322  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  IF (abs(eig_diff) == 0.0_dp) THEN
1325  pole = delta(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i) - 0.5_dp)
1326  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  admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(j))
1330  tmp = heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(i) - 0.5_dp)
1331  tmp = tmp - heaviside(admm_env%eigvals_P_to_be_purified(ispin)%eigvals%data(j) - 0.5_dp)
1332  pole = tmp*pole
1333  CALL cp_fm_set_element(admm_env%M_purify(ispin), i, j, pole)
1334  END IF
1335  END DO
1336  END DO
1337  CALL cp_fm_upper_to_full(admm_env%M_purify(ispin), admm_env%work_aux_aux)
1338 
1339  CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
1340  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  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  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  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  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  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  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  admm_env%work_orb_orb)
1371 
1372  NULLIFY (matrix_k_tilde)
1373  ALLOCATE (matrix_k_tilde)
1374  CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
1375  name='MATRIX K_tilde', &
1376  matrix_type=dbcsr_type_symmetric)
1377 
1378  CALL cp_fm_to_fm(admm_env%work_orb_orb, admm_env%ks_to_be_merged(ispin))
1379 
1380  CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
1381  CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
1382  CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.true.)
1383 
1384  IF (admm_env%block_dm) THEN
1385  ! ** now loop through the list and nullify blocks
1386  CALL dbcsr_iterator_start(iter, matrix_k_tilde)
1387  DO WHILE (dbcsr_iterator_blocks_left(iter))
1388  CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk)
1389  IF (admm_env%block_map(iatom, jatom) == 0) THEN
1390  sparse_block = 0.0_dp
1391  END IF
1392  END DO
1393  CALL dbcsr_iterator_stop(iter)
1394  END IF
1395 
1396  CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
1397 
1398  CALL dbcsr_deallocate_matrix(matrix_k_tilde)
1399 
1400  END DO !spin-loop
1401 
1402  CALL timestop(handle)
1403 
1404  END SUBROUTINE merge_ks_matrix_cauchy
1405 
1406 ! **************************************************************************************************
1407 !> \brief ...
1408 !> \param qs_env ...
1409 ! **************************************************************************************************
1410  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  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  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
1422 
1423  CALL timeset(routinen, handle)
1424  NULLIFY (admm_env, dft_control, matrix_ks, matrix_ks_aux_fit, mos, mos_aux_fit, &
1425  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  mos=mos)
1432  CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit, mos_aux_fit=mos_aux_fit)
1433 
1434  DO ispin = 1, dft_control%nspins
1435  nao_aux_fit = admm_env%nao_aux_fit
1436  nao_orb = admm_env%nao_orb
1437  nmo = admm_env%nmo(ispin)
1438  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1439  CALL get_mo_set(mo_set=mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
1440 
1441  !! Calculate Lambda^{-2}
1442  CALL cp_fm_to_fm(admm_env%lambda(ispin), admm_env%work_nmo_nmo1(ispin))
1443  CALL cp_fm_cholesky_decompose(admm_env%work_nmo_nmo1(ispin))
1444  CALL cp_fm_cholesky_invert(admm_env%work_nmo_nmo1(ispin))
1445  !! Symmetrize the guy
1446  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  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  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  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  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  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  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  admm_env%work_aux_aux2)
1480 
1481  CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
1482  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  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  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  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  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  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  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  admm_env%work_orb_orb)
1516 
1517  NULLIFY (matrix_k_tilde)
1518  ALLOCATE (matrix_k_tilde)
1519  CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
1520  name='MATRIX K_tilde', &
1521  matrix_type=dbcsr_type_symmetric)
1522 
1523  CALL cp_fm_to_fm(admm_env%work_orb_orb, admm_env%ks_to_be_merged(ispin))
1524 
1525  CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
1526  CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
1527  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  admm_env%mo_derivs_tmp(ispin))
1532 
1533  CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
1534 
1535  CALL dbcsr_deallocate_matrix(matrix_k_tilde)
1536 
1537  END DO !spin loop
1538  CALL timestop(handle)
1539 
1540  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  SUBROUTINE merge_mo_derivs_diag(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, mo_derivs, &
1561  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  REAL(dp), DIMENSION(:), POINTER :: occupation_numbers, scaling_factor
1575 
1576  CALL timeset(routinen, handle)
1577 
1578  nao_aux_fit = admm_env%nao_aux_fit
1579  nao_orb = admm_env%nao_orb
1580  nmo = admm_env%nmo(ispin)
1581 
1582  CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
1583  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  admm_env%H(ispin))
1588 
1589  CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
1590  ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
1591  scaling_factor = 2.0_dp*occupation_numbers
1592 
1593  CALL cp_fm_column_scale(admm_env%H(ispin), scaling_factor)
1594 
1595  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  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  admm_env%mo_derivs_tmp(ispin))
1604 
1605  ! *** Construct Matrix M for Hadamard Product
1606  pole = 0.0_dp
1607  DO i = 1, nmo
1608  DO j = i, nmo
1609  eig_diff = (admm_env%eigvals_lambda(ispin)%eigvals%data(i) - &
1610  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  IF (abs(eig_diff) < 0.0001_dp) THEN
1613  tmp32 = 1.0_dp/sqrt(admm_env%eigvals_lambda(ispin)%eigvals%data(j))**3
1614  tmp52 = tmp32/admm_env%eigvals_lambda(ispin)%eigvals%data(j)*eig_diff
1615  tmp72 = tmp52/admm_env%eigvals_lambda(ispin)%eigvals%data(j)*eig_diff
1616  tmp92 = tmp72/admm_env%eigvals_lambda(ispin)%eigvals%data(j)*eig_diff
1617 
1618  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  CALL cp_fm_set_element(admm_env%M(ispin), i, j, pole)
1620  ELSE
1621  pole = 1.0_dp/sqrt(admm_env%eigvals_lambda(ispin)%eigvals%data(i))
1622  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  admm_env%eigvals_lambda(ispin)%eigvals%data(j))
1625  CALL cp_fm_set_element(admm_env%M(ispin), i, j, pole)
1626  END IF
1627  END DO
1628  END DO
1629  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  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  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  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  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  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  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  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  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  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  admm_env%mo_derivs_tmp(ispin))
1681 
1682  DO i = 1, SIZE(scaling_factor)
1683  scaling_factor(i) = 1.0_dp/scaling_factor(i)
1684  END DO
1685 
1686  CALL cp_fm_column_scale(admm_env%mo_derivs_tmp(ispin), scaling_factor)
1687 
1688  CALL cp_fm_scale_and_add(1.0_dp, mo_derivs(ispin), 1.0_dp, admm_env%mo_derivs_tmp(ispin))
1689 
1690  DEALLOCATE (scaling_factor)
1691 
1692  CALL timestop(handle)
1693 
1694  END SUBROUTINE merge_mo_derivs_diag
1695 
1696 ! **************************************************************************************************
1697 !> \brief ...
1698 !> \param qs_env ...
1699 ! **************************************************************************************************
1700  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  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
1713  matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_s, matrix_s_aux_fit, rho_ao, &
1714  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  CALL timeset(routinen, handle)
1724  NULLIFY (admm_env, dft_control, matrix_ks, matrix_ks_aux_fit, matrix_ks_aux_fit_dft, &
1725  matrix_ks_aux_fit_hfx, matrix_s, matrix_s_aux_fit, rho_ao, rho_ao_aux, matrix_k_tilde, &
1726  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  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  matrix_s_aux_fit=matrix_s_aux_fit)
1739 
1740  CALL qs_rho_get(rho, rho_ao=rho_ao)
1741  CALL qs_rho_get(rho_aux_fit, &
1742  rho_ao=rho_ao_aux)
1743 
1744  DO ispin = 1, dft_control%nspins
1745  IF (admm_env%block_dm) THEN
1746  CALL dbcsr_iterator_start(iter, matrix_ks_aux_fit(ispin)%matrix)
1747  DO WHILE (dbcsr_iterator_blocks_left(iter))
1748  CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk)
1749  IF (admm_env%block_map(iatom, jatom) == 0) THEN
1750  sparse_block = 0.0_dp
1751  END IF
1752  END DO
1753  CALL dbcsr_iterator_stop(iter)
1754  CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_ks_aux_fit(ispin)%matrix, 1.0_dp, 1.0_dp)
1755 
1756  ELSE
1757 
1758  nao_aux_fit = admm_env%nao_aux_fit
1759  nao_orb = admm_env%nao_orb
1760  nmo = admm_env%nmo(ispin)
1761 
1762  ! ADMMS: different matrix for calculating A^(T)*K*A, see Eq. (37) Merlot
1763  IF (admm_env%do_admms) THEN
1764  NULLIFY (matrix_ks_aux_fit_admms_tmp)
1765  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  name='matrix_ks_aux_fit_admms_tmp', matrix_type='s')
1768  ! matrix_ks_aux_fit_admms_tmp = k(d_Q)
1769  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  1.0_dp, -(admm_env%gsi(ispin))**(2.0_dp/3.0_dp))
1774  CALL copy_dbcsr_to_fm(matrix_ks_aux_fit_admms_tmp, admm_env%K(ispin))
1775  CALL dbcsr_deallocate_matrix(matrix_ks_aux_fit_admms_tmp)
1776  ELSE
1777  CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
1778  END IF
1779 
1780  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  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  admm_env%work_orb_orb)
1790 
1791  NULLIFY (matrix_k_tilde)
1792  ALLOCATE (matrix_k_tilde)
1793  CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
1794  name='MATRIX K_tilde', matrix_type='S')
1795  CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
1796  CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
1797  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  IF (admm_env%do_admmq .OR. admm_env%do_admms) THEN
1802  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  IF (admm_env%do_admmp) THEN
1807  gsi_square = (admm_env%gsi(ispin))*(admm_env%gsi(ispin))
1808  CALL dbcsr_scale(matrix_k_tilde, gsi_square)
1809  END IF
1810 
1811  admm_env%lambda_merlot(ispin) = 0
1812 
1813  ! Calculate LAMBDA according to Merlot, 1. IF: ADMMQ, 2. IF: ADMMP, 3. IF: ADMMS,
1814  IF (admm_env%do_admmq) THEN
1815  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  admm_env%lambda_merlot(ispin) = trace_tmp/(admm_env%n_large_basis(ispin))
1820 
1821  ELSE IF (admm_env%do_admmp) THEN
1822  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  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  (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  /(admm_env%n_large_basis(ispin))
1834  END IF
1835 
1836  ELSE IF (admm_env%do_admms) THEN
1837  CALL dbcsr_dot(matrix_ks_aux_fit_hfx(ispin)%matrix, rho_ao_aux(ispin)%matrix, trace_tmp)
1838  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  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  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  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  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  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  CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix, admm_env%work_aux_aux4)
1863  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  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  admm_env%work_orb_orb3)
1873 
1874  NULLIFY (matrix_ttst)
1875  ALLOCATE (matrix_ttst)
1876  CALL dbcsr_create(matrix_ttst, template=matrix_ks(ispin)%matrix, &
1877  name='MATRIX TtsT', matrix_type='S')
1878  CALL dbcsr_copy(matrix_ttst, matrix_ks(ispin)%matrix)
1879  CALL dbcsr_set(matrix_ttst, 0.0_dp)
1880  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  (-admm_env%lambda_merlot(ispin))*admm_env%gsi(ispin))
1886 
1887  CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_s(1)%matrix, 1.0_dp, admm_env%lambda_merlot(ispin))
1888 
1889  CALL dbcsr_deallocate_matrix(matrix_ttst)
1890 
1891  END IF
1892 
1893  CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
1894 
1895  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  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  IF (dft_control%nspins == 2) THEN
1906  energy%exc_aux_fit = 0.0_dp
1907  energy%exc1_aux_fit = 0.0_dp
1908  energy%ex = 0.0_dp
1909  DO ispin = 1, dft_control%nspins
1910  energy%exc_aux_fit = energy%exc_aux_fit + (admm_env%gsi(ispin))**2.0_dp*ener_x(ispin)
1911  energy%exc1_aux_fit = energy%exc1_aux_fit + (admm_env%gsi(ispin))**2.0_dp*ener_x1(ispin)
1912  energy%ex = energy%ex + (admm_env%gsi(ispin))**2.0_dp*ener_k(ispin)
1913  END DO
1914  ELSE
1915  energy%exc_aux_fit = (admm_env%gsi(1))**2.0_dp*energy%exc_aux_fit
1916  energy%exc1_aux_fit = (admm_env%gsi(1))**2.0_dp*energy%exc1_aux_fit
1917  energy%ex = (admm_env%gsi(1))**2.0_dp*energy%ex
1918  END IF
1919 
1920  ELSE IF (admm_env%do_admms) THEN
1921  IF (dft_control%nspins == 2) THEN
1922  energy%exc_aux_fit = 0.0_dp
1923  energy%exc1_aux_fit = 0.0_dp
1924  DO ispin = 1, dft_control%nspins
1925  energy%exc_aux_fit = energy%exc_aux_fit + (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)*ener_x(ispin)
1926  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  energy%exc_aux_fit = (admm_env%gsi(1))**(2.0_dp/3.0_dp)*energy%exc_aux_fit
1930  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  CALL timestop(handle)
1935 
1936  END SUBROUTINE merge_ks_matrix_none
1937 
1938 ! **************************************************************************************************
1939 !> \brief ...
1940 !> \param qs_env ...
1941 ! **************************************************************************************************
1942  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  INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1953  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  REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
1958  TYPE(admm_type), POINTER :: admm_env
1959  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  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fmwork
1967  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_ks
1968  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_k_tilde, matrix_ks_aux_fit, &
1969  matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_kp, matrix_s, matrix_s_aux_fit, &
1970  rho_ao_aux
1971  TYPE(dbcsr_type) :: tmpmatrix_ks
1972  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  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  CALL timeset(routinen, handle)
1984  NULLIFY (admm_env, rho_ao_aux, rho_aux_fit, &
1985  matrix_s_aux_fit, energy, &
1986  para_env, kpoints, sab_aux_fit, &
1987  matrix_k_tilde, matrix_ks_kp, matrix_ks_aux_fit, scf_env, &
1988  struct_orb_orb, struct_aux_orb, struct_aux_aux, kp, &
1989  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  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  rho_aux_fit=rho_aux_fit)
2009  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  cell_to_index=cell_to_index)
2014 
2015  nao_aux_fit = admm_env%nao_aux_fit
2016  nao_orb = admm_env%nao_orb
2017  nspins = dft_control%nspins
2018 
2019  !Case study on ADMMQ, ADMMS and ADMMP
2020 
2021  !ADMMQ: calculate lamda as in Merlot eq (28)
2022  IF (admm_env%do_admmq) THEN
2023  admm_env%lambda_merlot = 0.0_dp
2024  DO img = 1, dft_control%nimages
2025  DO ispin = 1, nspins
2026  CALL dbcsr_dot(matrix_ks_aux_fit(ispin, img)%matrix, rho_ao_aux(ispin, img)%matrix, trace_tmp)
2027  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  IF (admm_env%do_admmp) THEN
2034  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  /(admm_env%n_large_basis(1))
2038  ELSE
2039  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  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  (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  IF (admm_env%do_admms) THEN
2052  IF (nspins == 1) THEN
2053  trace_tmp = 0.0_dp
2054  trace_tmp_two = 0.0_dp
2055  DO img = 1, dft_control%nimages
2056  CALL dbcsr_dot(matrix_ks_aux_fit_hfx(1, img)%matrix, rho_ao_aux(1, img)%matrix, tmp)
2057  trace_tmp = trace_tmp + tmp
2058  CALL dbcsr_dot(matrix_ks_aux_fit_dft(1, img)%matrix, rho_ao_aux(1, img)%matrix, tmp)
2059  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  trace_tmp_two))/(admm_env%n_large_basis(1))
2064  ELSE
2065 
2066  DO ispin = 1, nspins
2067  trace_tmp = 0.0_dp
2068  trace_tmp_two = 0.0_dp
2069  DO img = 1, dft_control%nimages
2070  CALL dbcsr_dot(matrix_ks_aux_fit_hfx(ispin, img)%matrix, rho_ao_aux(ispin, img)%matrix, tmp)
2071  trace_tmp = trace_tmp + tmp
2072  CALL dbcsr_dot(matrix_ks_aux_fit_dft(ispin, img)%matrix, rho_ao_aux(ispin, img)%matrix, tmp)
2073  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  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  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  NULLIFY (matrix_ks_aux_fit)
2089  ALLOCATE (matrix_ks_aux_fit(nspins, dft_control%nimages))
2090  DO img = 1, dft_control%nimages
2091  DO ispin = 1, nspins
2092  NULLIFY (matrix_ks_aux_fit(ispin, img)%matrix)
2093  ALLOCATE (matrix_ks_aux_fit(ispin, img)%matrix)
2094  CALL dbcsr_create(matrix_ks_aux_fit(ispin, img)%matrix, template=matrix_s_aux_fit(1, 1)%matrix)
2095  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  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  ALLOCATE (ksmatrix(2))
2104  CALL dbcsr_create(ksmatrix(1), template=matrix_ks_aux_fit(1, 1)%matrix, &
2105  matrix_type=dbcsr_type_symmetric)
2106  CALL dbcsr_create(ksmatrix(2), template=matrix_ks_aux_fit(1, 1)%matrix, &
2107  matrix_type=dbcsr_type_antisymmetric)
2108  CALL dbcsr_create(tmpmatrix_ks, template=matrix_ks_aux_fit(1, 1)%matrix, &
2109  matrix_type=dbcsr_type_symmetric)
2110  CALL cp_dbcsr_alloc_block_from_nbl(ksmatrix(1), sab_aux_fit)
2111  CALL cp_dbcsr_alloc_block_from_nbl(ksmatrix(2), sab_aux_fit)
2112 
2113  kplocal = kp_range(2) - kp_range(1) + 1
2114  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  nrow_global=nao_aux_fit, ncol_global=nao_aux_fit)
2118  CALL cp_fm_create(work_aux_aux, struct_aux_aux)
2119  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  nrow_global=nao_aux_fit, ncol_global=nao_orb)
2123  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  nrow_global=nao_orb, ncol_global=nao_orb)
2127 
2128  !Create cfm work matrices
2129  IF (.NOT. use_real_wfn) THEN
2130  CALL cp_cfm_create(cs, struct_aux_aux)
2131  CALL cp_cfm_create(ck, struct_aux_aux)
2132  CALL cp_cfm_create(cwork_aux_aux, struct_aux_aux)
2133 
2134  CALL cp_cfm_create(ca, struct_aux_orb)
2135  CALL cp_cfm_create(cwork_aux_orb, struct_aux_orb)
2136 
2137  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  ALLOCATE (fm_ks(kplocal, 2, nspins))
2142  DO ispin = 1, nspins
2143  DO i = 1, 2
2144  DO ikp = 1, kplocal
2145  CALL cp_fm_create(fm_ks(ikp, i, ispin), struct_orb_orb)
2146  END DO
2147  END DO
2148  END DO
2149 
2150  CALL cp_fm_struct_release(struct_aux_aux)
2151  CALL cp_fm_struct_release(struct_aux_orb)
2152  CALL cp_fm_struct_release(struct_orb_orb)
2153 
2154  ALLOCATE (info(kplocal*nspins*nkp_groups, 2))
2155  indx = 0
2156  DO ikp = 1, kplocal
2157  DO ispin = 1, nspins
2158  DO igroup = 1, nkp_groups
2159  ! number of current kpoint
2160  ik = kp_dist(1, igroup) + ikp - 1
2161  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2162  indx = indx + 1
2163 
2164  IF (use_real_wfn) THEN
2165  CALL dbcsr_set(ksmatrix(1), 0.0_dp)
2166  CALL rskp_transform(rmatrix=ksmatrix(1), rsmat=matrix_ks_aux_fit, ispin=ispin, &
2167  xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
2168  CALL dbcsr_desymmetrize(ksmatrix(1), tmpmatrix_ks)
2169  CALL copy_dbcsr_to_fm(tmpmatrix_ks, admm_env%work_aux_aux)
2170  ELSE
2171  CALL dbcsr_set(ksmatrix(1), 0.0_dp)
2172  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  xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
2175  CALL dbcsr_desymmetrize(ksmatrix(1), tmpmatrix_ks)
2176  CALL copy_dbcsr_to_fm(tmpmatrix_ks, admm_env%work_aux_aux)
2177  CALL dbcsr_desymmetrize(ksmatrix(2), tmpmatrix_ks)
2178  CALL copy_dbcsr_to_fm(tmpmatrix_ks, admm_env%work_aux_aux2)
2179  END IF
2180 
2181  IF (my_kpgrp) THEN
2182  CALL cp_fm_start_copy_general(admm_env%work_aux_aux, work_aux_aux, para_env, info(indx, 1))
2183  IF (.NOT. use_real_wfn) &
2184  CALL cp_fm_start_copy_general(admm_env%work_aux_aux2, work_aux_aux2, &
2185  para_env, info(indx, 2))
2186  ELSE
2187  CALL cp_fm_start_copy_general(admm_env%work_aux_aux, fmdummy, para_env, info(indx, 1))
2188  IF (.NOT. use_real_wfn) &
2189  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  DO ikp = 1, kplocal
2197  DO ispin = 1, nspins
2198  DO igroup = 1, nkp_groups
2199  ! number of current kpoint
2200  ik = kp_dist(1, igroup) + ikp - 1
2201  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2202  indx = indx + 1
2203  IF (my_kpgrp) THEN
2204  CALL cp_fm_finish_copy_general(work_aux_aux, info(indx, 1))
2205  IF (.NOT. use_real_wfn) THEN
2206  CALL cp_fm_finish_copy_general(work_aux_aux2, info(indx, 2))
2207  CALL cp_fm_to_cfm(work_aux_aux, work_aux_aux2, ck)
2208  END IF
2209  END IF
2210  END DO
2211 
2212  kp => kpoints%kp_aux_env(ikp)%kpoint_env
2213  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  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  fm_ks(ikp, 1, ispin))
2223  ELSE
2224 
2225  IF (admm_env%do_admmq .OR. admm_env%do_admms) THEN
2226  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  fac = cmplx(-admm_env%lambda_merlot(ispin), 0.0_dp, dp)
2230  CALL cp_cfm_scale_and_add(z_one, ck, fac, cs)
2231  CALL cp_cfm_scale(admm_env%gsi(ispin), ck)
2232  END IF
2233 
2234  IF (admm_env%do_admmp) THEN
2235  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  fac = cmplx(-admm_env%gsi(ispin)*admm_env%lambda_merlot(ispin), 0.0_dp, dp)
2239  fac2 = cmplx(admm_env%gsi(ispin)**2, 0.0_dp, dp)
2240  CALL cp_cfm_scale_and_add(fac2, ck, fac, cs)
2241  END IF
2242 
2243  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  z_one, ck, ca, z_zero, cwork_aux_orb)
2246 
2247  CALL parallel_gemm('C', 'N', nao_orb, nao_orb, nao_aux_fit, &
2248  z_one, ca, cwork_aux_orb, z_zero, cwork_orb_orb)
2249 
2250  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  DO ikp = 1, kplocal
2257  DO ispin = 1, nspins
2258  DO igroup = 1, nkp_groups
2259  ! number of current kpoint
2260  ik = kp_dist(1, igroup) + ikp - 1
2261  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2262  indx = indx + 1
2263  CALL cp_fm_cleanup_copy_general(info(indx, 1))
2264  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  DEALLOCATE (info)
2270  CALL dbcsr_release(ksmatrix(1))
2271  CALL dbcsr_release(ksmatrix(2))
2272  CALL dbcsr_release(tmpmatrix_ks)
2273 
2274  CALL cp_fm_release(work_aux_aux)
2275  CALL cp_fm_release(work_aux_aux2)
2276  CALL cp_fm_release(work_aux_orb)
2277  IF (.NOT. use_real_wfn) THEN
2278  CALL cp_cfm_release(cs)
2279  CALL cp_cfm_release(ck)
2280  CALL cp_cfm_release(cwork_aux_aux)
2281  CALL cp_cfm_release(ca)
2282  CALL cp_cfm_release(cwork_aux_orb)
2283  CALL cp_cfm_release(cwork_orb_orb)
2284  END IF
2285 
2286  NULLIFY (matrix_k_tilde)
2287 
2288  CALL dbcsr_allocate_matrix_set(matrix_k_tilde, dft_control%nspins, dft_control%nimages)
2289 
2290  DO ispin = 1, nspins
2291  DO img = 1, dft_control%nimages
2292  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  matrix_type=dbcsr_type_symmetric, nze=0)
2296  CALL cp_dbcsr_alloc_block_from_nbl(matrix_k_tilde(ispin, img)%matrix, sab_kp)
2297  CALL dbcsr_set(matrix_k_tilde(ispin, img)%matrix, 0.0_dp)
2298  END DO
2299  END DO
2300 
2301  CALL cp_fm_get_info(admm_env%work_orb_orb, matrix_struct=struct_orb_orb)
2302  ALLOCATE (fmwork(2))
2303  CALL cp_fm_create(fmwork(1), struct_orb_orb)
2304  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  fmwork, for_aux_fit=.false., pmat_ext=fm_ks)
2310  CALL cp_fm_release(fmwork(1))
2311  CALL cp_fm_release(fmwork(2))
2312 
2313  DO ispin = 1, nspins
2314  DO i = 1, 2
2315  DO ikp = 1, kplocal
2316  CALL cp_fm_release(fm_ks(ikp, i, ispin))
2317  END DO
2318  END DO
2319  END DO
2320 
2321  DO ispin = 1, nspins
2322  DO img = 1, dft_control%nimages
2323  CALL dbcsr_add(matrix_ks_kp(ispin, img)%matrix, matrix_k_tilde(ispin, img)%matrix, 1.0_dp, 1.0_dp)
2324  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  1.0_dp, admm_env%lambda_merlot(ispin))
2328  END IF
2329  END DO
2330  END DO
2331 
2332  !Scale the energies
2333  IF (admm_env%do_admmp) THEN
2334  IF (nspins == 1) THEN
2335  energy%exc_aux_fit = (admm_env%gsi(1))**2.0_dp*energy%exc_aux_fit
2336  energy%exc1_aux_fit = (admm_env%gsi(1))**2.0_dp*energy%exc1_aux_fit
2337  energy%ex = (admm_env%gsi(1))**2.0_dp*energy%ex
2338  ELSE
2339  energy%exc_aux_fit = 0.0_dp
2340  energy%exc1_aux_fit = 0.0_dp
2341  energy%ex = 0.0_dp
2342  DO ispin = 1, dft_control%nspins
2343  energy%exc_aux_fit = energy%exc_aux_fit + (admm_env%gsi(ispin))**2.0_dp*ener_x(ispin)
2344  energy%exc1_aux_fit = energy%exc1_aux_fit + (admm_env%gsi(ispin))**2.0_dp*ener_x1(ispin)
2345  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  IF (admm_env%do_admms) THEN
2352  IF (nspins == 1) THEN
2353  energy%exc_aux_fit = (admm_env%gsi(1))**(2.0_dp/3.0_dp)*energy%exc_aux_fit
2354  energy%exc1_aux_fit = (admm_env%gsi(1))**(2.0_dp/3.0_dp)*energy%exc1_aux_fit
2355  ELSE
2356  energy%exc_aux_fit = 0.0_dp
2357  energy%exc1_aux_fit = 0.0_dp
2358  DO ispin = 1, nspins
2359  energy%exc_aux_fit = energy%exc_aux_fit + (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)*ener_x(ispin)
2360  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  CALL dbcsr_deallocate_matrix_set(matrix_ks_aux_fit)
2365  END IF
2366 
2367  CALL dbcsr_deallocate_matrix_set(matrix_k_tilde)
2368 
2369  CALL timestop(handle)
2370 
2371  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
2396  TYPE(admm_gapw_r3d_rs_type), POINTER :: admm_gapw_env
2397  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2398  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
2399  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit_hfx, rho_ao_aux, &
2400  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  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
2405  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  CALL timeset(routinen, handle)
2412 
2413  NULLIFY (ks_env, rho_aux_fit, rho_aux_fit_buffer, rho_ao, &
2414  xc_section_aux, v_rspace_dummy, v_tau_rspace_dummy, &
2415  rho_ao_aux, rho_ao_aux_buffer, dft_control, &
2416  matrix_ks_aux_fit_hfx, task_list, local_rho_buffer, admm_gapw_env)
2417 
2418  NULLIFY (rho_g, rho_r, tot_rho_r)
2419 
2420  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  matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx)
2423 
2424  CALL qs_rho_get(rho_aux_fit, &
2425  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  tot_rho_r=tot_rho_r)
2432 
2433  gapw = admm_env%do_gapw
2434  nimg = dft_control%nimages
2435 
2436 ! Calculate rho_buffer = rho_aux(ispin) to get exchange of ispin electrons
2437  DO img = 1, nimg
2438  CALL dbcsr_set(rho_ao_aux_buffer(1, img)%matrix, 0.0_dp)
2439  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  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  basis_type = "AUX_FIT"
2446  task_list => admm_env%task_list_aux_fit
2447  IF (gapw) THEN
2448  basis_type = "AUX_FIT_SOFT"
2449  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  DO myspin = 1, dft_control%nspins
2454 
2455  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  task_list_external=task_list)
2464 
2465  END DO
2466 
2467  ! Write changes in buffer density matrix
2468  CALL qs_rho_set(rho_aux_fit_buffer, rho_r_valid=.true., rho_g_valid=.true.)
2469 
2470  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  just_energy=.true.)
2477 
2478  !atomic contributions: use the atomic density as stored in admm_env%gapw_env
2479  ener_x1_ispin = 0.0_dp
2480  IF (gapw) THEN
2481 
2482  admm_gapw_env => admm_env%admm_gapw_env
2483  CALL get_qs_env(qs_env, &
2484  atomic_kind_set=atomic_kind_set, &
2485  para_env=para_env)
2486 
2487  CALL local_rho_set_create(local_rho_buffer)
2488  CALL allocate_rho_atom_internals(local_rho_buffer%rho_atom_set, atomic_kind_set, &
2489  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  para_env=para_env)
2496 
2497  CALL prepare_gapw_den(qs_env, local_rho_set=local_rho_buffer, do_rho0=.false., &
2498  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  rho_atom_set_external=local_rho_buffer%rho_atom_set)
2504 
2505  CALL local_rho_set_release(local_rho_buffer)
2506  END IF
2507 
2508  ener_k_ispin = 0.0_dp
2509 
2510  !! ** Calculate the exchange energy
2511  DO img = 1, nimg
2512  CALL dbcsr_dot(matrix_ks_aux_fit_hfx(ispin, img)%matrix, rho_ao_aux_buffer(ispin, img)%matrix, tmp)
2513  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  ener_k_ispin = ener_k_ispin/2.0_dp
2519 
2520  CALL timestop(handle)
2521 
2522  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  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  CALL timeset(routinen, handle)
2543 
2544  NULLIFY (admm_env, dft_control)
2545 
2546  CALL get_qs_env(qs_env, &
2547  admm_env=admm_env, &
2548  dft_control=dft_control)
2549 
2550  ! only for ADMMP
2551  IF (admm_env%do_admmp) THEN
2552  DO ispin = 1, dft_control%nspins
2553  DO img = 1, dft_control%nimages
2554  IF (scale_back) THEN
2555  CALL dbcsr_scale(rho_ao_orb(ispin, img)%matrix, 1.0_dp/admm_env%gsi(ispin))
2556  ELSE
2557  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  CALL timestop(handle)
2564 
2565  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  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  REAL(dp), DIMENSION(:), POINTER :: occupation_numbers, scaling_factor
2584  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit, &
2585  matrix_ks_aux_fit_dft, &
2586  matrix_ks_aux_fit_hfx
2587  TYPE(dbcsr_type) :: dbcsr_work
2588 
2589  NULLIFY (matrix_ks_aux_fit, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx)
2590 
2591  CALL timeset(routinen, handle)
2592 
2593  nao_aux_fit = admm_env%nao_aux_fit
2594  nao_orb = admm_env%nao_orb
2595  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  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  IF (admm_env%do_admms) THEN
2608  !In ADMMS, we use the K matrix defined as K_hf - gsi^2/3*K_dft
2609  CALL dbcsr_create(dbcsr_work, template=matrix_ks_aux_fit(ispin)%matrix)
2610  CALL dbcsr_copy(dbcsr_work, matrix_ks_aux_fit_hfx(ispin)%matrix)
2611  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  CALL copy_dbcsr_to_fm(dbcsr_work, admm_env%K(ispin))
2613  CALL dbcsr_release(dbcsr_work)
2614  ELSE
2615  CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
2616  END IF
2617  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  admm_env%H(ispin))
2622 
2623  CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
2624  ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
2625 
2626  scaling_factor = 2.0_dp*occupation_numbers
2627 
2628  CALL cp_fm_column_scale(admm_env%H(ispin), scaling_factor)
2629 
2630  DEALLOCATE (scaling_factor)
2631 
2632  CALL timestop(handle)
2633 
2634  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  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  REAL(dp), DIMENSION(:), POINTER :: occupation_numbers, scaling_factor
2655 
2656  CALL timeset(routinen, handle)
2657 
2658  nao_aux_fit = admm_env%nao_aux_fit
2659  nao_orb = admm_env%nao_orb
2660  nmo = admm_env%nmo(ispin)
2661 
2662  CALL copy_dbcsr_to_fm(matrix_ks_aux_fit(ispin)%matrix, admm_env%K(ispin))
2663  CALL cp_fm_upper_to_full(admm_env%K(ispin), admm_env%work_aux_aux)
2664 
2665  CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
2666  ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
2667  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  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  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  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  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  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  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  admm_env%mo_derivs_tmp(ispin))
2692 
2693  CALL cp_fm_column_scale(admm_env%mo_derivs_tmp(ispin), scaling_factor)
2694 
2695  CALL cp_fm_scale_and_add(1.0_dp, mo_derivs(ispin), 1.0_dp, admm_env%mo_derivs_tmp(ispin))
2696 
2697  DEALLOCATE (scaling_factor)
2698 
2699  CALL timestop(handle)
2700 
2701  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  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  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_derivs_fm
2716  TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_derivs_aux_fit
2717  TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit
2718  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit
2719  TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array, mos_aux_fit
2720 
2721  NULLIFY (mo_array, mos_aux_fit, matrix_ks_aux_fit, mo_coeff_aux_fit, &
2722  mo_derivs_aux_fit, mo_coeff)
2723 
2724  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  matrix_ks_aux_fit=matrix_ks_aux_fit)
2727 
2728  nspins = SIZE(mo_derivs)
2729  ALLOCATE (mo_derivs_fm(nspins))
2730  DO ispin = 1, nspins
2731  CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff)
2732  CALL cp_fm_create(mo_derivs_fm(ispin), mo_coeff%matrix_struct)
2733  END DO
2734 
2735  DO ispin = 1, nspins
2736  CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff)
2737  CALL get_mo_set(mo_set=mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
2738 
2739  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  mo_derivs_fm, mo_derivs_aux_fit, matrix_ks_aux_fit)
2742  CALL copy_fm_to_dbcsr(mo_derivs_fm(ispin), mo_derivs(ispin)%matrix)
2743  END DO
2744 
2745  CALL cp_fm_release(mo_derivs_fm)
2746 
2747  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  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  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  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
2762  TYPE(mo_set_type), POINTER :: mo_set
2763 
2764  CALL get_qs_env(qs_env, dft_control=dft_control)
2765 
2766  IF (dft_control%do_admm_dm) THEN
2767  cpabort("Forces with ADMM DM methods not implemented")
2768  END IF
2769  IF (dft_control%do_admm_mo .AND. .NOT. qs_env%run_rtp) THEN
2770  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  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  matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb)
2776  DO ispin = 1, dft_control%nspins
2777  mo_set => mos(ispin)
2778  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  IF (admm_env%purification_method == do_admm_purify_none) THEN
2781  CALL get_mo_set(mo_set=mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
2782  CALL calc_aux_mo_derivs_none(ispin, qs_env%admm_env, mo_set, mo_coeff_aux_fit)
2783  END IF
2784  END DO
2785  CALL calc_mixed_overlap_force(qs_env)
2786  END IF
2787 
2788  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  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  INTEGER, DIMENSION(:, :), POINTER :: kp_dist
2805  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2806  LOGICAL :: gapw, my_kpgrp, use_real_wfn
2807  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: admm_force
2808  REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
2809  TYPE(admm_type), POINTER :: admm_env
2810  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2811  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  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_skap, fm_skapa
2821  TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
2822  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit, matrix_ks_aux_fit_dft, &
2823  matrix_ks_aux_fit_hfx, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, matrix_skap, &
2824  matrix_skapa, rho_ao_orb
2825  TYPE(dbcsr_type) :: kmatrix_tmp
2826  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  POINTER :: sab_aux_fit, sab_aux_fit_asymm, &
2833  sab_aux_fit_vs_orb, sab_kp
2834  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  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  NULLIFY (ks_env, admm_env, matrix_ks_aux_fit, &
2847  matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, rho, force, &
2848  para_env, atomic_kind_set, kpoints, sab_aux_fit, &
2849  sab_aux_fit_vs_orb, sab_aux_fit_asymm, struct_orb_orb, &
2850  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  rho=rho)
2861  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  matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx)
2871 
2872  gapw = admm_env%do_gapw
2873  nao_aux_fit = admm_env%nao_aux_fit
2874  nao_orb = admm_env%nao_orb
2875  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  cell_to_index=cell_to_index, sab_nl=sab_kp)
2880 
2881  !Case study on ADMMQ, ADMMS and ADMMP
2882  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  NULLIFY (matrix_ks_aux_fit)
2885  ALLOCATE (matrix_ks_aux_fit(nspins, dft_control%nimages))
2886  DO img = 1, dft_control%nimages
2887  DO ispin = 1, nspins
2888  NULLIFY (matrix_ks_aux_fit(ispin, img)%matrix)
2889  ALLOCATE (matrix_ks_aux_fit(ispin, img)%matrix)
2890  CALL dbcsr_create(matrix_ks_aux_fit(ispin, img)%matrix, template=matrix_s_aux_fit(1, 1)%matrix)
2891  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  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  ALLOCATE (kmatrix(2))
2901  CALL dbcsr_create(kmatrix(1), template=matrix_ks_aux_fit(1, 1)%matrix, &
2902  matrix_type=dbcsr_type_symmetric)
2903  CALL dbcsr_create(kmatrix(2), template=matrix_ks_aux_fit(1, 1)%matrix, &
2904  matrix_type=dbcsr_type_antisymmetric)
2905  CALL dbcsr_create(kmatrix_tmp, template=matrix_ks_aux_fit(1, 1)%matrix, &
2906  matrix_type=dbcsr_type_no_symmetry)
2907  CALL cp_dbcsr_alloc_block_from_nbl(kmatrix(1), sab_aux_fit)
2908  CALL cp_dbcsr_alloc_block_from_nbl(kmatrix(2), sab_aux_fit)
2909 
2910  kplocal = kp_range(2) - kp_range(1) + 1
2911  para_env => kpoints%blacs_env_all%para_env
2912  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  nrow_global=nao_aux_fit, ncol_global=nao_aux_fit)
2916  CALL cp_fm_create(work_aux_aux, struct_aux_aux)
2917  CALL cp_fm_create(work_aux_aux2, struct_aux_aux)
2918  CALL cp_fm_create(work_aux_aux3, struct_aux_aux)
2919  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  nrow_global=nao_aux_fit, ncol_global=nao_orb)
2923  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  nrow_global=nao_orb, ncol_global=nao_orb)
2927 
2928  !Create cfm work matrices
2929  IF (.NOT. use_real_wfn) THEN
2930  CALL cp_cfm_create(cpmatrix, struct_orb_orb)
2931 
2932  CALL cp_cfm_create(cs_inv, struct_aux_aux)
2933  CALL cp_cfm_create(cs, struct_aux_aux)
2934  CALL cp_cfm_create(cwork_aux_aux, struct_aux_aux)
2935  CALL cp_cfm_create(ckmatrix, struct_aux_aux)
2936 
2937  CALL cp_cfm_create(ca, struct_aux_orb)
2938  CALL cp_cfm_create(cq, struct_aux_orb)
2939  CALL cp_cfm_create(cwork_aux_orb, struct_aux_orb)
2940  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  ALLOCATE (fm_skap(kplocal, 2, nspins), fm_skapa(kplocal, 2, nspins))
2945  DO ispin = 1, nspins
2946  DO i = 1, 2
2947  DO ikp = 1, kplocal
2948  CALL cp_fm_create(fm_skap(ikp, i, ispin), struct_aux_orb)
2949  CALL cp_fm_create(fm_skapa(ikp, i, ispin), struct_aux_aux)
2950  END DO
2951  END DO
2952  END DO
2953 
2954  CALL cp_fm_struct_release(struct_aux_aux)
2955  CALL cp_fm_struct_release(struct_aux_orb)
2956  CALL cp_fm_struct_release(struct_orb_orb)
2957 
2958  indx = 0
2959  DO ikp = 1, kplocal
2960  DO ispin = 1, nspins
2961  DO igroup = 1, nkp_groups
2962  ! number of current kpoint
2963  ik = kp_dist(1, igroup) + ikp - 1
2964  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
2965  indx = indx + 1
2966 
2967  ! FT of matrices KS, then transfer to FM type
2968  IF (use_real_wfn) THEN
2969  CALL dbcsr_set(kmatrix(1), 0.0_dp)
2970  CALL rskp_transform(rmatrix=kmatrix(1), rsmat=matrix_ks_aux_fit, ispin=ispin, &
2971  xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
2972  CALL dbcsr_desymmetrize(kmatrix(1), kmatrix_tmp)
2973  CALL copy_dbcsr_to_fm(kmatrix_tmp, admm_env%work_aux_aux)
2974  ELSE
2975  CALL dbcsr_set(kmatrix(1), 0.0_dp)
2976  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  xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
2979  CALL dbcsr_desymmetrize(kmatrix(1), kmatrix_tmp)
2980  CALL copy_dbcsr_to_fm(kmatrix_tmp, admm_env%work_aux_aux)
2981  CALL dbcsr_desymmetrize(kmatrix(2), kmatrix_tmp)
2982  CALL copy_dbcsr_to_fm(kmatrix_tmp, admm_env%work_aux_aux2)
2983  END IF
2984 
2985  IF (my_kpgrp) THEN
2986  CALL cp_fm_start_copy_general(admm_env%work_aux_aux, work_aux_aux, para_env, info(indx, 1))
2987  IF (.NOT. use_real_wfn) &
2988  CALL cp_fm_start_copy_general(admm_env%work_aux_aux2, work_aux_aux2, para_env, info(indx, 2))
2989  ELSE
2990  CALL cp_fm_start_copy_general(admm_env%work_aux_aux, fmdummy, para_env, info(indx, 1))
2991  IF (.NOT. use_real_wfn) &
2992  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  DO ikp = 1, kplocal
3000  DO ispin = 1, nspins
3001  DO igroup = 1, nkp_groups
3002  ! number of current kpoint
3003  ik = kp_dist(1, igroup) + ikp - 1
3004  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
3005  indx = indx + 1
3006  IF (my_kpgrp) THEN
3007  CALL cp_fm_finish_copy_general(work_aux_aux, info(indx, 1))
3008  IF (.NOT. use_real_wfn) THEN
3009  CALL cp_fm_finish_copy_general(work_aux_aux2, info(indx, 2))
3010  CALL cp_fm_to_cfm(work_aux_aux, work_aux_aux2, ckmatrix)
3011  END IF
3012  END IF
3013  END DO
3014  kp => kpoints%kp_aux_env(ikp)%kpoint_env
3015 
3016  IF (use_real_wfn) THEN
3017 
3018  !! Calculate S'_inverse
3019  CALL cp_fm_to_fm(kp%smat(1, 1), s_inv)
3020  CALL cp_fm_cholesky_decompose(s_inv)
3021  CALL cp_fm_cholesky_invert(s_inv)
3022  !! Symmetrize the guy
3023  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  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  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  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  kp%amat(1, 1), 0.0_dp, fm_skapa(ikp, 1, ispin))
3035 
3036  ELSE !complex wfn
3037 
3038  IF (admm_env%do_admmq .OR. admm_env%do_admms) THEN
3039  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  fac = cmplx(-admm_env%lambda_merlot(ispin), 0.0_dp, dp)
3043  CALL cp_cfm_scale_and_add(z_one, ckmatrix, fac, cs)
3044  CALL cp_cfm_scale(admm_env%gsi(ispin), ckmatrix)
3045  END IF
3046 
3047  IF (admm_env%do_admmp) THEN
3048  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  fac = cmplx(-admm_env%gsi(ispin)*admm_env%lambda_merlot(ispin), 0.0_dp, dp)
3052  fac2 = cmplx(admm_env%gsi(ispin)**2, 0.0_dp, dp)
3053  CALL cp_cfm_scale_and_add(fac2, ckmatrix, fac, cs)
3054  END IF
3055 
3056  CALL cp_fm_to_cfm(kp%smat(1, 1), kp%smat(2, 1), cs_inv)
3057  CALL cp_cfm_cholesky_decompose(cs_inv)
3058  CALL cp_cfm_cholesky_invert(cs_inv)
3059  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  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  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  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  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  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  ca, z_zero, cwork_aux_aux)
3077 
3078  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  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  z_zero, cwork_aux_orb)
3086  CALL parallel_gemm('N', 'C', nao_aux_fit, nao_aux_fit, nao_orb, fac, cwork_aux_orb, &
3087  ca, z_one, cwork_aux_aux)
3088  END IF
3089 
3090  CALL cp_cfm_to_fm(cwork_aux_orb2, mtargetr=fm_skap(ikp, 1, ispin), mtargeti=fm_skap(ikp, 2, ispin))
3091  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  DO ikp = 1, kplocal
3100  DO ispin = 1, nspins
3101  DO igroup = 1, nkp_groups
3102  ! number of current kpoint
3103  ik = kp_dist(1, igroup) + ikp - 1
3104  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
3105  indx = indx + 1
3106  CALL cp_fm_cleanup_copy_general(info(indx, 1))
3107  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  DEALLOCATE (info)
3113  CALL dbcsr_release(kmatrix(1))
3114  CALL dbcsr_release(kmatrix(2))
3115  CALL dbcsr_release(kmatrix_tmp)
3116 
3117  CALL cp_fm_release(work_aux_aux)
3118  CALL cp_fm_release(work_aux_aux2)
3119  CALL cp_fm_release(work_aux_aux3)
3120  CALL cp_fm_release(s_inv)
3121  CALL cp_fm_release(work_aux_orb)
3122  IF (.NOT. use_real_wfn) THEN
3123  CALL cp_cfm_release(ckmatrix)
3124  CALL cp_cfm_release(cpmatrix)
3125  CALL cp_cfm_release(cs_inv)
3126  CALL cp_cfm_release(cs)
3127  CALL cp_cfm_release(cwork_aux_aux)
3128  CALL cp_cfm_release(cwork_aux_orb)
3129  CALL cp_cfm_release(cwork_aux_orb2)
3130  CALL cp_cfm_release(ca)
3131  CALL cp_cfm_release(cq)
3132  END IF
3133 
3134  !Back FT to real space
3135  ALLOCATE (matrix_skap(nspins, nimg), matrix_skapa(nspins, nimg))
3136  DO img = 1, nimg
3137  DO ispin = 1, nspins
3138  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  matrix_type=dbcsr_type_no_symmetry)
3141  CALL cp_dbcsr_alloc_block_from_nbl(matrix_skap(ispin, img)%matrix, sab_aux_fit_vs_orb)
3142 
3143  ALLOCATE (matrix_skapa(ispin, img)%matrix)
3144  CALL dbcsr_create(matrix_skapa(ispin, img)%matrix, template=matrix_s_aux_fit(1, 1)%matrix, &
3145  matrix_type=dbcsr_type_no_symmetry)
3146  CALL cp_dbcsr_alloc_block_from_nbl(matrix_skapa(ispin, img)%matrix, sab_aux_fit_asymm)
3147  END DO
3148  END DO
3149 
3150  ALLOCATE (fmwork(2))
3151  CALL cp_fm_get_info(admm_env%work_aux_orb, matrix_struct=struct_aux_orb)
3152  CALL cp_fm_create(fmwork(1), struct_aux_orb)
3153  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  fmwork, for_aux_fit=.true., pmat_ext=fm_skap)
3157  CALL cp_fm_release(fmwork(1))
3158  CALL cp_fm_release(fmwork(2))
3159 
3160  CALL cp_fm_get_info(admm_env%work_aux_aux, matrix_struct=struct_aux_aux)
3161  CALL cp_fm_create(fmwork(1), struct_aux_aux)
3162  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  fmwork, for_aux_fit=.true., pmat_ext=fm_skapa)
3166  CALL cp_fm_release(fmwork(1))
3167  CALL cp_fm_release(fmwork(2))
3168  DEALLOCATE (fmwork)
3169 
3170  DO img = 1, nimg
3171  DO ispin = 1, nspins
3172  CALL dbcsr_scale(matrix_skap(ispin, img)%matrix, -2.0_dp)
3173  CALL dbcsr_scale(matrix_skapa(ispin, img)%matrix, 2.0_dp)
3174  END DO
3175  IF (nspins == 2) THEN
3176  CALL dbcsr_add(matrix_skap(1, img)%matrix, matrix_skap(2, img)%matrix, 1.0_dp, 1.0_dp)
3177  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  ALLOCATE (admm_force(3, natom))
3182  admm_force = 0.0_dp
3183 
3184  IF (admm_env%do_admmq .OR. admm_env%do_admmp .OR. admm_env%do_admms) THEN
3185  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_orb)
3186  DO img = 1, nimg
3187  DO ispin = 1, nspins
3188  CALL dbcsr_scale(rho_ao_orb(ispin, img)%matrix, -admm_env%lambda_merlot(ispin))
3189  END DO
3190  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  sab_nl=sab_kp, matrixkp_p=rho_ao_orb(1, :))
3196  DO img = 1, nimg
3197  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  DO ispin = 1, nspins
3199  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  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  sab_nl=sab_aux_fit_asymm, matrixkp_p=matrix_skapa(1, :))
3208 
3209  CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set)
3210  DEALLOCATE (admm_force)
3211 
3212  DO ispin = 1, nspins
3213  DO i = 1, 2
3214  DO ikp = 1, kplocal
3215  CALL cp_fm_release(fm_skap(ikp, i, ispin))
3216  CALL cp_fm_release(fm_skapa(ikp, i, ispin))
3217  END DO
3218  END DO
3219  END DO
3220  CALL dbcsr_deallocate_matrix_set(matrix_skap)
3221  CALL dbcsr_deallocate_matrix_set(matrix_skapa)
3222 
3223  IF (admm_env%do_admms) THEN
3224  CALL dbcsr_deallocate_matrix_set(matrix_ks_aux_fit)
3225  END IF
3226 
3227  CALL timestop(handle)
3228 
3229  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  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  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3250  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  POINTER :: sab_aux_fit_asymm, sab_aux_fit_vs_orb
3254  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3255  TYPE(qs_ks_env_type), POINTER :: ks_env
3256 
3257  CALL timeset(routinen, handle)
3258 
3259  cpassert(ASSOCIATED(qs_env))
3260 
3261  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  matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, sab_aux_fit_vs_orb=sab_aux_fit_vs_orb)
3264 
3265  my_fval = 2.0_dp
3266  IF (PRESENT(fval)) my_fval = fval
3267 
3268  ALLOCATE (matrix_w_q)
3269  CALL dbcsr_copy(matrix_w_q, matrix_s_aux_fit_vs_orb(1)%matrix, &
3270  "W MATRIX AUX Q")
3271  CALL cp_dbcsr_alloc_block_from_nbl(matrix_w_q, sab_aux_fit_vs_orb)
3272  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  matrix_type=dbcsr_type_no_symmetry)
3276  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  natom=natom, force=force)
3280  ALLOCATE (admm_force(3, natom))
3281  admm_force = 0.0_dp
3282 
3283  nspins = SIZE(matrix_pz)
3284  nao = admm_env%nao_orb
3285  naux = admm_env%nao_aux_fit
3286 
3287  CALL cp_fm_set_all(admm_env%work_aux_orb2, 0.0_dp)
3288 
3289  DO ispin = 1, nspins
3290  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  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  admm_env%A, 0.0_dp, admm_env%work_aux_orb)
3295  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  admm_env%work_orb_orb, 1.0_dp, admm_env%work_aux_orb2)
3299  END DO
3300 
3301  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  admm_env%A, 0.0_dp, admm_env%work_aux_aux)
3306  CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_w_s, keep_sparsity=.true.)
3307 
3308  CALL dbcsr_scale(matrix_w_q, -my_fval)
3309  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  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  sab_nl=sab_aux_fit_vs_orb, matrix_p=matrix_w_q)
3317 
3318  ! add forces
3319  CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set)
3320 
3321  DEALLOCATE (admm_force)
3322  CALL dbcsr_deallocate_matrix(matrix_w_s)
3323  CALL dbcsr_deallocate_matrix(matrix_w_q)
3324 
3325  CALL timestop(handle)
3326 
3327  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: admm_force
3366  TYPE(admm_type), POINTER :: admm_env
3367  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_aux_fit, &
3371  matrix_s_aux_fit_vs_orb, rho_ao, &
3372  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  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  POINTER :: sab_orb
3380  TYPE(qs_energy_type), POINTER :: energy
3381  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  CALL timeset(routinen, handle)
3386 
3387  NULLIFY (admm_env, logger, dft_control, para_env, mos, mo_coeff, matrix_w_q, matrix_w_s, &
3388  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  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  matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb)
3403 
3404  CALL qs_rho_get(rho, rho_ao=rho_ao)
3405  CALL qs_rho_get(rho_aux_fit, &
3406  rho_ao=rho_ao_aux)
3407 
3408  nao_aux_fit = admm_env%nao_aux_fit
3409  nao_orb = admm_env%nao_orb
3410 
3411  logger => cp_get_default_logger()
3412 
3413  ! *** forces are only implemented for mo_diag or none and basis_projection ***
3414  IF (admm_env%block_dm) THEN
3415  cpabort("ADMM Forces not implemented for blocked projection methods!")
3416  END IF
3417 
3418  IF (.NOT. (admm_env%purification_method == do_admm_purify_mo_diag .OR. &
3419  admm_env%purification_method == do_admm_purify_none)) THEN
3420  cpabort("ADMM Forces only implemented without purification or for MO_DIAG.")
3421  END IF
3422 
3423  ! *** Create sparse work matrices
3424 
3425  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  matrix_type=dbcsr_type_no_symmetry)
3429  CALL cp_dbcsr_alloc_block_from_nbl(matrix_w_s, admm_env%sab_aux_fit_asymm)
3430 
3431  ALLOCATE (matrix_w_q)
3432  CALL dbcsr_copy(matrix_w_q, matrix_s_aux_fit_vs_orb(1)%matrix, &
3433  "W MATRIX AUX Q")
3434 
3435  DO ispin = 1, dft_control%nspins
3436  nmo = admm_env%nmo(ispin)
3437  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
3438 
3439  ! *** S'^(-T)*H'
3440  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  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  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  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  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  admm_env%work_aux_aux)
3465 
3466  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  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  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  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  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  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  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  IF (admm_env%do_admms) THEN
3498  ! *** scale admm_env%work_aux_orb by gsi due to inner derivative
3499  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  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  admm_env%work_aux_orb)
3508 
3509  ! ADMMP
3510  ELSE IF (admm_env%do_admmp) THEN
3511  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  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  admm_env%work_aux_orb)
3521 
3522  ! ADMMQ
3523  ELSE IF (admm_env%do_admmq) THEN
3524  ! *** scale admm_env%work_aux_orb by gsi due to inner derivative
3525  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  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  admm_env%work_aux_orb)
3534  END IF
3535 
3536  ! *** copy to sparse matrix
3537  CALL copy_fm_to_dbcsr(admm_env%work_aux_orb, matrix_w_q, keep_sparsity=.true.)
3538 
3539  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  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  admm_env%work_aux_aux)
3548  END IF
3549 
3550  ! *** copy to sparse matrix
3551  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  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  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  matrix_type=dbcsr_type_no_symmetry)
3562 
3563  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  IF (admm_env%do_admms .OR. admm_env%do_admmq) THEN
3568  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  -admm_env%lambda_merlot(ispin))
3571 
3572  ! ADMMP scale by gsi^2 and add derivative of variational term with resp. to s
3573  ELSE IF (admm_env%do_admmp) THEN
3574 
3575  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  (-admm_env%gsi(ispin))*admm_env%lambda_merlot(ispin))
3578 
3579  END IF
3580 
3581  CALL dbcsr_deallocate_matrix(matrix_rho_aux_desymm_tmp)
3582 
3583  END IF
3584 
3585  ! allocate force vector
3586  CALL get_qs_env(qs_env=qs_env, natom=natom)
3587  ALLOCATE (admm_force(3, natom))
3588  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  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  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  IF (admm_env%do_admmq .OR. admm_env%do_admmp .OR. admm_env%do_admms) THEN
3598  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  sab_nl=sab_orb, matrix_p=rho_ao(ispin)%matrix)
3602  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  force=force)
3608  CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set)
3609  DEALLOCATE (admm_force)
3610 
3611  CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
3612  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  extension=".Log")
3616  CALL cp_dbcsr_write_sparse_matrix(matrix_w_s, 4, 6, qs_env, &
3617  para_env, output_unit=iw, omit_headers=omit_headers)
3618  CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
3619  "DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT")
3620  END IF
3621  IF (btest(cp_print_key_should_output(logger%iter_info, &
3622  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  extension=".Log")
3625  CALL cp_dbcsr_write_sparse_matrix(matrix_w_q, 4, 6, qs_env, &
3626  para_env, output_unit=iw, omit_headers=omit_headers)
3627  CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
3628  "DFT%PRINT%AO_MATRICES/W_MATRIX_AUX_FIT")
3629  END IF
3630 
3631  END DO !spin loop
3632 
3633  ! *** Deallocated weighted density matrices
3634  CALL dbcsr_deallocate_matrix(matrix_w_s)
3635  CALL dbcsr_deallocate_matrix(matrix_w_q)
3636 
3637  CALL timestop(handle)
3638 
3639  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  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  CALL timeset(routinen, handle)
3668 
3669  CALL dbcsr_set(density_matrix, 0.0_dp)
3670  nao_aux_fit = admm_env%nao_aux_fit
3671  nmo = admm_env%nmo(ispin)
3672  CALL cp_fm_to_fm(admm_env%C_hat(ispin), admm_env%work_aux_nmo(ispin))
3673  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  admm_env%work_aux_nmo2(ispin))
3678 
3679  ! The following IF doesn't do anything unless !alpha=mo_set%maxocc is uncommented.
3680  IF (.NOT. mo_set%uniform_occupation) THEN ! not all orbitals 1..homo are equally occupied
3681  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  alpha=alpha)
3687  ELSE
3688  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  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  IF (admm_env%do_admmp .OR. admm_env%do_admmq .OR. admm_env%do_admms) THEN
3701 
3702  CALL cite_reference(merlot2014)
3703 
3704  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  CALL dbcsr_dot(density_matrix_large, overlap_matrix_large, admm_env%n_large_basis(ispin))
3709  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  CALL dbcsr_dot(density_matrix, overlap_matrix, nel_tmp_aux)
3712  admm_env%gsi(ispin) = admm_env%n_large_basis(ispin)/nel_tmp_aux
3713 
3714  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  CALL dbcsr_scale(density_matrix, admm_env%gsi(ispin))
3717  END IF
3718 
3719  END IF
3720 
3721  CALL timestop(handle)
3722 
3723  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  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  REAL(dp), DIMENSION(:, :), POINTER :: sparse_block, sparse_block_aux
3744  TYPE(dbcsr_iterator_type) :: iter
3745 
3746  CALL timeset(routinen, handle)
3747 
3748  ! ** set blocked density matrix to 0
3749  CALL dbcsr_set(density_matrix_aux, 0.0_dp)
3750 
3751  ! ** now loop through the list and copy corresponding blocks
3752  CALL dbcsr_iterator_start(iter, density_matrix)
3753  DO WHILE (dbcsr_iterator_blocks_left(iter))
3754  CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk)
3755  IF (admm_env%block_map(iatom, jatom) == 1) THEN
3756  CALL dbcsr_get_block_p(density_matrix_aux, &
3757  row=iatom, col=jatom, block=sparse_block_aux, found=found)
3758  IF (found) THEN
3759  sparse_block_aux = sparse_block
3760  END IF
3761 
3762  END IF
3763  END DO
3764  CALL dbcsr_iterator_stop(iter)
3765 
3766  CALL copy_dbcsr_to_fm(density_matrix_aux, admm_env%P_to_be_purified(ispin))
3767  CALL cp_fm_upper_to_full(admm_env%P_to_be_purified(ispin), admm_env%work_orb_orb2)
3768 
3769  IF (nspins == 1) THEN
3770  CALL cp_fm_scale(0.5_dp, admm_env%P_to_be_purified(ispin))
3771  END IF
3772 
3773  CALL timestop(handle)
3774  END SUBROUTINE blockify_density_matrix
3775 
3776 ! **************************************************************************************************
3777 !> \brief ...
3778 !> \param x ...
3779 !> \return ...
3780 ! **************************************************************************************************
3781  ELEMENTAL FUNCTION delta(x)
3782  REAL(kind=dp), INTENT(IN) :: x
3783  REAL(kind=dp) :: delta
3784 
3785  IF (x == 0.0_dp) THEN !TODO: exact comparison of reals?
3786  delta = 1.0_dp
3787  ELSE
3788  delta = 0.0_dp
3789  END IF
3790 
3791  END FUNCTION delta
3792 
3793 ! **************************************************************************************************
3794 !> \brief ...
3795 !> \param x ...
3796 !> \return ...
3797 ! **************************************************************************************************
3798  ELEMENTAL FUNCTION heaviside(x)
3799  REAL(kind=dp), INTENT(IN) :: x
3800  REAL(kind=dp) :: heaviside
3801 
3802  IF (x < 0.0_dp) THEN
3803  heaviside = 0.0_dp
3804  ELSE
3805  heaviside = 1.0_dp
3806  END IF
3807  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  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  CALL timeset(routinen, handle)
3827 
3828  CALL get_qs_env(qs_env, admm_env=admm_env, dft_control=dft_control)
3829 
3830  nspins = dft_control%nspins
3831 
3832  cpassert(ASSOCIATED(admm_env%A))
3833  cpassert(ASSOCIATED(admm_env%work_orb_orb))
3834  cpassert(ASSOCIATED(admm_env%work_aux_orb))
3835  cpassert(ASSOCIATED(admm_env%work_aux_aux))
3836  CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux, ncol_global=nao)
3837 
3838  ! P1 -> AUX BASIS
3839  CALL cp_fm_get_info(admm_env%work_orb_orb, nrow_global=nao, ncol_global=ncol)
3840  DO ispin = 1, nspins
3841  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  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  admm_env%work_aux_orb, 0.0_dp, admm_env%work_aux_aux)
3846  CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, dm_admm(ispin)%matrix, keep_sparsity=.true.)
3847  END DO
3848 
3849  CALL timestop(handle)
3850 
3851  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  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  INTEGER, DIMENSION(:, :), POINTER :: kp_dist
3867  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3868  LOGICAL :: my_kpgrp, use_real_wfn
3869  REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
3870  TYPE(admm_type), POINTER :: admm_env
3871  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  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fmwork
3880  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_aux_fit, matrix_s_aux_fit_vs_orb
3881  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  POINTER :: sab_aux_fit, sab_aux_fit_vs_orb
3887 
3888  NULLIFY (xkp, kp_dist, para_env_local, cell_to_index, admm_env, kp, &
3889  kpoints, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, sab_aux_fit, sab_aux_fit_vs_orb, &
3890  para_env_global, matrix_struct_aux_fit, matrix_struct_aux_fit_vs_orb)
3891 
3892  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  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  nkp_groups=nkp_groups, kp_dist=kp_dist, cell_to_index=cell_to_index)
3901  kplocal = kp_range(2) - kp_range(1) + 1
3902  nc = 1
3903  IF (.NOT. use_real_wfn) nc = 2
3904 
3905  ALLOCATE (dbcsr_aux_fit(3))
3906  CALL dbcsr_create(dbcsr_aux_fit(1), template=matrix_s_aux_fit(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
3907  CALL dbcsr_create(dbcsr_aux_fit(2), template=matrix_s_aux_fit(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
3908  CALL dbcsr_create(dbcsr_aux_fit(3), template=matrix_s_aux_fit(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
3909  CALL cp_dbcsr_alloc_block_from_nbl(dbcsr_aux_fit(1), sab_aux_fit)
3910  CALL cp_dbcsr_alloc_block_from_nbl(dbcsr_aux_fit(2), sab_aux_fit)
3911 
3912  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  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  matrix_type=dbcsr_type_no_symmetry)
3917  CALL cp_dbcsr_alloc_block_from_nbl(dbcsr_aux_fit_vs_orb(1), sab_aux_fit_vs_orb)
3918  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  nao_aux_fit = admm_env%nao_aux_fit
3922  nao_orb = admm_env%nao_orb
3923  para_env_global => kpoints%blacs_env_all%para_env
3924 
3925  ALLOCATE (fmwork(4))
3926  CALL cp_fm_struct_create(matrix_struct_aux_fit, context=kpoints%blacs_env_all, para_env=para_env_global, &
3927  nrow_global=nao_aux_fit, ncol_global=nao_aux_fit)
3928  CALL cp_fm_create(fmwork(1), matrix_struct_aux_fit)
3929  CALL cp_fm_create(fmwork(2), matrix_struct_aux_fit)
3930  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  nrow_global=nao_aux_fit, ncol_global=nao_orb)
3934  CALL cp_fm_create(fmwork(3), matrix_struct_aux_fit_vs_orb)
3935  CALL cp_fm_create(fmwork(4), matrix_struct_aux_fit_vs_orb)
3936  CALL cp_fm_struct_release(matrix_struct_aux_fit_vs_orb)
3937 
3938  !Create fm local to the KP groups
3939  nao_aux_fit = admm_env%nao_aux_fit
3940  nao_orb = admm_env%nao_orb
3941  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  nrow_global=nao_aux_fit, ncol_global=nao_aux_fit)
3945  CALL cp_fm_create(rmat_aux_fit, matrix_struct_aux_fit)
3946  CALL cp_fm_create(imat_aux_fit, matrix_struct_aux_fit)
3947  CALL cp_fm_create(work_aux_fit, matrix_struct_aux_fit)
3948  CALL cp_cfm_create(cwork_aux_fit, matrix_struct_aux_fit)
3949  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  nrow_global=nao_aux_fit, ncol_global=nao_orb)
3953  CALL cp_fm_create(rmat_aux_fit_vs_orb, matrix_struct_aux_fit_vs_orb)
3954  CALL cp_fm_create(imat_aux_fit_vs_orb, matrix_struct_aux_fit_vs_orb)
3955  CALL cp_cfm_create(cwork_aux_fit_vs_orb, matrix_struct_aux_fit_vs_orb)
3956  CALL cp_cfm_create(cmat_aux_fit_vs_orb, matrix_struct_aux_fit_vs_orb)
3957 
3958  ALLOCATE (info(kplocal*nkp_groups, 4))
3959 
3960  ! Steup and start all the communication
3961  indx = 0
3962  DO ikp = 1, kplocal
3963  DO igroup = 1, nkp_groups
3964  ik = kp_dist(1, igroup) + ikp - 1
3965  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
3966  indx = indx + 1
3967 
3968  IF (use_real_wfn) THEN
3969  !AUX-AUX overlap
3970  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  xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
3973  CALL dbcsr_desymmetrize(dbcsr_aux_fit(1), dbcsr_aux_fit(3))
3974  CALL copy_dbcsr_to_fm(dbcsr_aux_fit(3), fmwork(1))
3975 
3976  !AUX-ORB overlap
3977  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  xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit_vs_orb)
3980  CALL copy_dbcsr_to_fm(dbcsr_aux_fit_vs_orb(1), fmwork(3))
3981  ELSE
3982  !AUX-AUX overlap
3983  CALL dbcsr_set(dbcsr_aux_fit(1), 0.0_dp)
3984  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  ispin=1, xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_aux_fit)
3987  CALL dbcsr_desymmetrize(dbcsr_aux_fit(1), dbcsr_aux_fit(3))
3988  CALL copy_dbcsr_to_fm(dbcsr_aux_fit(3), fmwork(1))
3989  CALL dbcsr_desymmetrize(dbcsr_aux_fit(2), dbcsr_aux_fit(3))
3990  CALL copy_dbcsr_to_fm(dbcsr_aux_fit(3), fmwork(2))
3991 
3992  !AUX-ORB overlap
3993  CALL dbcsr_set(dbcsr_aux_fit_vs_orb(1), 0.0_dp)
3994  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  cell_to_index=cell_to_index, sab_nl=sab_aux_fit_vs_orb)
3998  CALL copy_dbcsr_to_fm(dbcsr_aux_fit_vs_orb(1), fmwork(3))
3999  CALL copy_dbcsr_to_fm(dbcsr_aux_fit_vs_orb(2), fmwork(4))
4000  END IF
4001 
4002  IF (my_kpgrp) THEN
4003  CALL cp_fm_start_copy_general(fmwork(1), rmat_aux_fit, para_env_global, info(indx, 1))
4004  CALL cp_fm_start_copy_general(fmwork(3), rmat_aux_fit_vs_orb, para_env_global, info(indx, 3))
4005  IF (.NOT. use_real_wfn) THEN
4006  CALL cp_fm_start_copy_general(fmwork(2), imat_aux_fit, para_env_global, info(indx, 2))
4007  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  CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env_global, info(indx, 1))
4011  CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env_global, info(indx, 3))
4012  IF (.NOT. use_real_wfn) THEN
4013  CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env_global, info(indx, 2))
4014  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  DO ikp = 1, kplocal
4024  DO igroup = 1, nkp_groups
4025  ik = kp_dist(1, igroup) + ikp - 1
4026  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
4027  indx = indx + 1
4028 
4029  IF (my_kpgrp) THEN
4030  CALL cp_fm_finish_copy_general(rmat_aux_fit, info(indx, 1))
4031  CALL cp_fm_finish_copy_general(rmat_aux_fit_vs_orb, info(indx, 3))
4032  IF (.NOT. use_real_wfn) THEN
4033  CALL cp_fm_finish_copy_general(imat_aux_fit, info(indx, 2))
4034  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  kp => kpoints%kp_aux_env(ikp)%kpoint_env
4040 
4041  !Allocate local KP matrices
4042  CALL cp_fm_release(kp%amat)
4043  ALLOCATE (kp%amat(nc, 1))
4044  DO ic = 1, nc
4045  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  IF (admm_env%do_admmp .OR. admm_env%do_admmq .OR. admm_env%do_admms .OR. calculate_forces) THEN
4050  CALL cp_fm_release(kp%smat)
4051  ALLOCATE (kp%smat(nc, 1))
4052  DO ic = 1, nc
4053  CALL cp_fm_create(kp%smat(ic, 1), matrix_struct_aux_fit)
4054  END DO
4055  CALL cp_fm_to_fm(rmat_aux_fit, kp%smat(1, 1))
4056  IF (.NOT. use_real_wfn) CALL cp_fm_to_fm(imat_aux_fit, kp%smat(2, 1))
4057  END IF
4058 
4059  IF (use_real_wfn) THEN
4060  !Invert S_aux
4061  CALL cp_fm_cholesky_decompose(rmat_aux_fit)
4062  CALL cp_fm_cholesky_invert(rmat_aux_fit)
4063  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  rmat_aux_fit, rmat_aux_fit_vs_orb, 0.0_dp, kp%amat(1, 1))
4068  ELSE
4069 
4070  !Invert S_aux
4071  CALL cp_fm_to_cfm(rmat_aux_fit, imat_aux_fit, cmat_aux_fit)
4072  CALL cp_cfm_cholesky_decompose(cmat_aux_fit)
4073  CALL cp_cfm_cholesky_invert(cmat_aux_fit)
4074  CALL own_cfm_upper_to_full(cmat_aux_fit, cwork_aux_fit)
4075 
4076  !A = S^-1 * Q
4077  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  cmat_aux_fit, cmat_aux_fit_vs_orb, z_zero, cwork_aux_fit_vs_orb)
4080  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  DO ikp = 1, kplocal
4087  DO igroup = 1, nkp_groups
4088  ik = kp_dist(1, igroup) + ikp - 1
4089  my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
4090  indx = indx + 1
4091 
4092  IF (my_kpgrp) THEN
4093  CALL cp_fm_cleanup_copy_general(info(indx, 1))
4094  CALL cp_fm_cleanup_copy_general(info(indx, 3))
4095  IF (.NOT. use_real_wfn) THEN
4096  CALL cp_fm_cleanup_copy_general(info(indx, 2))
4097  CALL cp_fm_cleanup_copy_general(info(indx, 4))
4098  END IF
4099  END IF
4100 
4101  END DO
4102  END DO
4103 
4104  CALL cp_fm_release(rmat_aux_fit)
4105  CALL cp_fm_release(imat_aux_fit)
4106  CALL cp_fm_release(work_aux_fit)
4107  CALL cp_cfm_release(cwork_aux_fit)
4108  CALL cp_cfm_release(cmat_aux_fit)
4109  CALL cp_fm_release(rmat_aux_fit_vs_orb)
4110  CALL cp_fm_release(imat_aux_fit_vs_orb)
4111  CALL cp_cfm_release(cwork_aux_fit_vs_orb)
4112  CALL cp_cfm_release(cmat_aux_fit_vs_orb)
4113  CALL cp_fm_struct_release(matrix_struct_aux_fit)
4114  CALL cp_fm_struct_release(matrix_struct_aux_fit_vs_orb)
4115 
4116  CALL cp_fm_release(fmwork(1))
4117  CALL cp_fm_release(fmwork(2))
4118  CALL cp_fm_release(fmwork(3))
4119  CALL cp_fm_release(fmwork(4))
4120 
4121  CALL dbcsr_release(dbcsr_aux_fit(1))
4122  CALL dbcsr_release(dbcsr_aux_fit(2))
4123  CALL dbcsr_release(dbcsr_aux_fit(3))
4124  CALL dbcsr_release(dbcsr_aux_fit_vs_orb(1))
4125  CALL dbcsr_release(dbcsr_aux_fit_vs_orb(2))
4126 
4127  END SUBROUTINE kpoint_calc_admm_matrices
4128 
4129 ! **************************************************************************************************
4130 !> \brief ...
4131 !> \param cfm_mat_Q ...
4132 !> \param cfm_mat_work ...
4133 ! **************************************************************************************************
4134  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  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
4143 
4144  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  col_indices=col_indices)
4154 
4155  DO jjb = 1, ncol_local
4156  j_global = col_indices(jjb)
4157  DO iib = 1, nrow_local
4158  i_global = row_indices(iib)
4159  IF (j_global < i_global) THEN
4160  cfm_mat_q%local_data(iib, jjb) = z_zero
4161  END IF
4162  IF (j_global == i_global) THEN
4163  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  CALL cp_cfm_transpose(cfm_mat_q, 'C', cfm_mat_work)
4169 
4170  CALL cp_cfm_scale_and_add(z_one, cfm_mat_q, z_one, cfm_mat_work)
4171 
4172  CALL timestop(handle)
4173 
4174  END SUBROUTINE
4175 
4176 END MODULE admm_methods
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
admm_methods
Contains ADMM methods which require molecular orbitals.
Definition: admm_methods.F:15
admm_methods::admm_mo_calc_rho_aux_kp
subroutine, public admm_mo_calc_rho_aux_kp(qs_env)
...
Definition: admm_methods.F:285
admm_methods::admm_mo_merge_derivs
subroutine, public admm_mo_merge_derivs(ispin, admm_env, mo_set, mo_coeff, mo_coeff_aux_fit, mo_derivs, mo_derivs_aux_fit, matrix_ks_aux_fit)
...
Definition: admm_methods.F:828
admm_methods::admm_mo_merge_ks_matrix
subroutine, public admm_mo_merge_ks_matrix(qs_env)
...
Definition: admm_methods.F:778
admm_methods::admm_update_ks_atom
subroutine, public admm_update_ks_atom(qs_env, calculate_forces)
Adds the GAPW exchange contribution to the aux_fit ks matrices.
Definition: admm_methods.F:715
admm_methods::calc_admm_ovlp_forces_kp
subroutine, public calc_admm_ovlp_forces_kp(qs_env)
Calculate the forces due to the AUX/ORB basis overlap in ADMM, in the KP case.
Definition: admm_methods.F:2795
admm_methods::admm_fit_mo_coeffs
subroutine, public admm_fit_mo_coeffs(admm_env, matrix_s_aux_fit, matrix_s_mixed, mos, mos_aux_fit, geometry_did_change)
...
Definition: admm_methods.F:870
admm_methods::admm_mo_calc_rho_aux
subroutine, public admm_mo_calc_rho_aux(qs_env)
...
Definition: admm_methods.F:149
admm_methods::calc_admm_ovlp_forces
subroutine, public calc_admm_ovlp_forces(qs_env)
Calculate the forces due to the AUX/ORB basis overlap in ADMM.
Definition: admm_methods.F:2754
admm_methods::admm_projection_derivative
subroutine, public admm_projection_derivative(qs_env, matrix_hz, matrix_pz, fval)
Calculate derivatives terms from overlap matrices.
Definition: admm_methods.F:3239
admm_methods::admm_aux_response_density
subroutine, public admm_aux_response_density(qs_env, dm, dm_admm)
Calculate ADMM auxiliary response density.
Definition: admm_methods.F:3816
admm_methods::scale_dm
subroutine, public scale_dm(qs_env, rho_ao_orb, scale_back)
Scale density matrix by gsi(ispin), is needed for force scaling in ADMMP.
Definition: admm_methods.F:2532
admm_methods::kpoint_calc_admm_matrices
subroutine, public kpoint_calc_admm_matrices(qs_env, calculate_forces)
Fill the ADMM overlp and basis change matrices in the KP env based on the real-space array.
Definition: admm_methods.F:3859
admm_methods::calc_admm_mo_derivatives
subroutine, public calc_admm_mo_derivatives(qs_env, mo_derivs)
Calculate the derivative of the AUX_FIT mo, based on the ORB mo_derivs.
Definition: admm_methods.F:2709
admm_methods::calc_mixed_overlap_force
subroutine, public calc_mixed_overlap_force(qs_env)
Calculates contribution of forces due to basis transformation.
Definition: admm_methods.F:3357
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition: admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition: admm_types.F:593
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::merlot2014
integer, save, public merlot2014
Definition: bibliography.F:43
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition: cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition: cp_cfm_basic_linalg.F:244
cp_cfm_basic_linalg::cp_cfm_transpose
subroutine, public cp_cfm_transpose(matrix, trans, matrixt)
Transposes a BLACS distributed complex matrix.
Definition: cp_cfm_basic_linalg.F:1239
cp_cfm_basic_linalg::cp_cfm_scale_and_add_fm
subroutine, public cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b). where b is a real matrix (adapted from cp_cf...
Definition: cp_cfm_basic_linalg.F:355
cp_cfm_basic_linalg::cp_cfm_cholesky_invert
subroutine, public cp_cfm_cholesky_invert(matrix, n, info_out)
Used to replace Cholesky decomposition by the inverse.
Definition: cp_cfm_basic_linalg.F:982
cp_cfm_basic_linalg::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition: cp_cfm_basic_linalg.F:926
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition: cp_cfm_types.F:12
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition: cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition: cp_cfm_types.F:159
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition: cp_cfm_types.F:817
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition: cp_cfm_types.F:607
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition: cp_cfm_types.F:765
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Definition: cp_dbcsr_operations.F:887
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition: cp_dbcsr_operations.F:118
cp_dbcsr_output
DBCSR output in CP2K.
Definition: cp_dbcsr_output.F:17
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers)
...
Definition: cp_dbcsr_output.F:163
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition: cp_fm_basic_linalg.F:1882
cp_fm_basic_linalg::cp_fm_upper_to_full
subroutine, public cp_fm_upper_to_full(matrix, work)
given an upper triangular matrix computes the corresponding full matrix
Definition: cp_fm_basic_linalg.F:1759
cp_fm_basic_linalg::cp_fm_schur_product
subroutine, public cp_fm_schur_product(matrix_a, matrix_b, matrix_c)
computes the schur product of two matrices c_ij = a_ij * b_ij
Definition: cp_fm_basic_linalg.F:728
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition: cp_fm_basic_linalg.F:167
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition: cp_fm_basic_linalg.F:1680
cp_fm_cholesky
various cholesky decomposition related routines
Definition: cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_invert
subroutine, public cp_fm_cholesky_invert(matrix, n, info_out)
used to replace the cholesky decomposition by the inverse
Definition: cp_fm_cholesky.F:122
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition: cp_fm_cholesky.F:50
cp_fm_cholesky::cp_fm_cholesky_reduce
subroutine, public cp_fm_cholesky_reduce(matrix, matrixb, itype)
reduce a matrix pencil A,B to normal form B has to be cholesky decomposed with cp_fm_cholesky_decompo...
Definition: cp_fm_cholesky.F:192
cp_fm_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa)
...
Definition: cp_fm_cholesky.F:254
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
cp_fm_diag::cp_fm_syevd
subroutine, public cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
Computes all eigenvalues and vectors of a real symmetric matrix significantly faster than syevx,...
Definition: cp_fm_diag.F:413
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_start_copy_general
subroutine, public cp_fm_start_copy_general(source, destination, para_env, info)
Initiates the copy operation: get distribution data, post MPI isend and irecvs.
Definition: cp_fm_types.F:1568
cp_fm_types::cp_fm_cleanup_copy_general
subroutine, public cp_fm_cleanup_copy_general(info)
Completes the copy operation: wait for comms clean up MPI state.
Definition: cp_fm_types.F:1946
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_finish_copy_general
subroutine, public cp_fm_finish_copy_general(destination, info)
Completes the copy operation: wait for comms, unpack, clean up MPI state.
Definition: cp_fm_types.F:1882
cp_fm_types::cp_fm_set_element
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
Definition: cp_fm_types.F:700
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_admm_purify_mo_no_diag
integer, parameter, public do_admm_purify_mo_no_diag
Definition: input_constants.F:848
input_constants::do_admm_purify_none
integer, parameter, public do_admm_purify_none
Definition: input_constants.F:848
input_constants::do_admm_purify_cauchy_subspace
integer, parameter, public do_admm_purify_cauchy_subspace
Definition: input_constants.F:848
input_constants::do_admm_purify_cauchy
integer, parameter, public do_admm_purify_cauchy
Definition: input_constants.F:848
input_constants::do_admm_purify_mo_diag
integer, parameter, public do_admm_purify_mo_diag
Definition: input_constants.F:848
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
kpoint_methods
Routines needed for kpoint calculation.
Definition: kpoint_methods.F:15
kpoint_methods::rskp_transform
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
Definition: kpoint_methods.F:780
kpoint_methods::kpoint_density_transform
subroutine, public kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext)
generate real space density matrices in DBCSR format
Definition: kpoint_methods.F:1114
kpoint_methods::kpoint_density_matrices
subroutine, public kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
Calculate kpoint density matrices (rho(k), owned by kpoint groups)
Definition: kpoint_methods.F:990
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
kpoint_types::get_kpoint_env
subroutine, public get_kpoint_env(kpoint_env, nkpoint, wkp, xkp, is_local, mos)
Get information from a single kpoint environment.
Definition: kpoint_types.F:747
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition: kpoint_types.F:333
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition: mathconstants.F:143
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition: mathconstants.F:142
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition: mathconstants.F:143
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
pw_types
Definition: pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition: qs_collocate_density.F:38
qs_collocate_density::calculate_rho_elec
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Definition: qs_collocate_density.F:1710
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_force_types
Definition: qs_force_types.F:13
qs_force_types::add_qs_force
subroutine, public add_qs_force(force, qs_force, forcetype, atomic_kind_set)
Add force to a force_type variable.
Definition: qs_force_types.F:425
qs_gapw_densities
Definition: qs_gapw_densities.F:9
qs_gapw_densities::prepare_gapw_den
subroutine, public prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external)
...
Definition: qs_gapw_densities.F:49
qs_ks_atom
routines that build the Kohn-Sham matrix contributions coming from local atomic densities
Definition: qs_ks_atom.F:12
qs_ks_atom::update_ks_atom
subroutine, public update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, kind_set_external, oce_external, sab_external, kscale, kintegral, kforce, fscale)
The correction to the KS matrix due to the GAPW local terms to the hartree and XC contributions is he...
Definition: qs_ks_atom.F:110
qs_ks_types
Definition: qs_ks_types.F:15
qs_local_rho_types
Definition: qs_local_rho_types.F:8
qs_local_rho_types::local_rho_set_create
subroutine, public local_rho_set_create(local_rho_set)
...
Definition: qs_local_rho_types.F:194
qs_local_rho_types::local_rho_set_release
subroutine, public local_rho_set_release(local_rho_set)
...
Definition: qs_local_rho_types.F:213
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition: qs_overlap.F:19
qs_overlap::build_overlap_force
subroutine, public build_overlap_force(ks_env, force, basis_type_a, basis_type_b, sab_nl, matrix_p, matrixkp_p)
Calculation of the force contribution from an overlap matrix over Cartesian Gaussian functions.
Definition: qs_overlap.F:784
qs_rho_atom_methods
Definition: qs_rho_atom_methods.F:7
qs_rho_atom_methods::allocate_rho_atom_internals
subroutine, public allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
Definition: qs_rho_atom_methods.F:916
qs_rho_atom_methods::calculate_rho_atom_coeff
subroutine, public calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
...
Definition: qs_rho_atom_methods.F:420
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_set
subroutine, public qs_rho_set(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
...
Definition: qs_rho_types.F:308
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
qs_scf_types
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14
qs_vxc_atom
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
Definition: qs_vxc_atom.F:12
qs_vxc_atom::calculate_vxc_atom
subroutine, public calculate_vxc_atom(qs_env, energy_only, exc1, gradient_atom_set, adiabatic_rescale_factor, kind_set_external, rho_atom_set_external, xc_section_external)
...
Definition: qs_vxc_atom.F:87
qs_vxc
Definition: qs_vxc.F:16
qs_vxc::qs_vxc_create
subroutine, public qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, just_energy, edisp, dispersion_env, adiabatic_rescale_factor, pw_env_external)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition: qs_vxc.F:98
task_list_types
types for task lists
Definition: task_list_types.F:14