(git:ed6f26b)
Loading...
Searching...
No Matches
admm_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 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! **************************************************************************************************
15MODULE admm_methods
16 USE admm_types, ONLY: admm_gapw_r3d_rs_type,&
17 admm_type,&
18 get_admm_env
19 USE atomic_kind_types, ONLY: atomic_kind_type
20 USE bibliography, ONLY: merlot2014,&
21 cite_reference
22 USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale,&
23 cp_cfm_scale_and_add,&
24 cp_cfm_scale_and_add_fm,&
25 cp_cfm_uplo_to_full
26 USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose,&
27 cp_cfm_cholesky_invert
28 USE cp_cfm_types, ONLY: cp_cfm_create,&
29 cp_cfm_release,&
30 cp_cfm_to_fm,&
31 cp_cfm_type,&
32 cp_fm_to_cfm
33 USE cp_control_types, ONLY: dft_control_type
34 USE cp_dbcsr_api, ONLY: &
35 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
36 dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
37 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
38 dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
39 dbcsr_type_no_symmetry, dbcsr_type_symmetric
40 USE cp_dbcsr_contrib, ONLY: dbcsr_dot
41 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
42 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
43 copy_fm_to_dbcsr,&
44 cp_dbcsr_plus_fm_fm_t,&
45 dbcsr_allocate_matrix_set,&
46 dbcsr_deallocate_matrix_set
47 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
48 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
49 cp_fm_scale,&
50 cp_fm_scale_and_add,&
51 cp_fm_schur_product,&
52 cp_fm_uplo_to_full
53 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
54 cp_fm_cholesky_invert,&
55 cp_fm_cholesky_reduce,&
56 cp_fm_cholesky_restore
57 USE cp_fm_diag, ONLY: cp_fm_syevd
58 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
59 cp_fm_struct_release,&
60 cp_fm_struct_type
61 USE cp_fm_types, ONLY: &
62 copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, &
63 cp_fm_get_info, cp_fm_release, cp_fm_set_all, cp_fm_set_element, cp_fm_start_copy_general, &
64 cp_fm_to_fm, cp_fm_type
65 USE cp_log_handling, ONLY: cp_get_default_logger,&
66 cp_logger_type,&
67 cp_to_string
68 USE cp_output_handling, ONLY: cp_p_file,&
69 cp_print_key_finished_output,&
70 cp_print_key_should_output,&
71 cp_print_key_unit_nr
72 USE input_constants, ONLY: do_admm_purify_cauchy,&
73 do_admm_purify_cauchy_subspace,&
74 do_admm_purify_mo_diag,&
75 do_admm_purify_mo_no_diag,&
76 do_admm_purify_none
77 USE input_section_types, ONLY: section_vals_type,&
78 section_vals_val_get
79 USE kinds, ONLY: default_string_length,&
80 dp
81 USE kpoint_methods, ONLY: kpoint_density_matrices,&
82 kpoint_density_transform,&
83 rskp_transform
84 USE kpoint_types, ONLY: get_kpoint_env,&
85 get_kpoint_info,&
86 kpoint_env_type,&
87 kpoint_type
88 USE mathconstants, ONLY: gaussi,&
89 z_one,&
90 z_zero
91 USE message_passing, ONLY: mp_para_env_type
92 USE parallel_gemm_api, ONLY: parallel_gemm
93 USE pw_types, ONLY: pw_c1d_gs_type,&
94 pw_r3d_rs_type
95 USE qs_collocate_density, ONLY: calculate_rho_elec
96 USE qs_energy_types, ONLY: qs_energy_type
97 USE qs_environment_types, ONLY: get_qs_env,&
98 qs_environment_type
99 USE qs_force_types, ONLY: add_qs_force,&
100 qs_force_type
101 USE qs_gapw_densities, ONLY: prepare_gapw_den
102 USE qs_ks_atom, ONLY: update_ks_atom
103 USE qs_ks_types, ONLY: qs_ks_env_type
104 USE qs_local_rho_types, ONLY: local_rho_set_create,&
105 local_rho_set_release,&
106 local_rho_type
107 USE qs_mo_types, ONLY: get_mo_set,&
108 mo_set_type
109 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
110 USE qs_overlap, ONLY: build_overlap_force
111 USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,&
112 calculate_rho_atom_coeff
113 USE qs_rho_types, ONLY: qs_rho_get,&
114 qs_rho_set,&
115 qs_rho_type
116 USE qs_scf_types, ONLY: qs_scf_env_type
117 USE qs_vxc, ONLY: qs_vxc_create
118 USE qs_vxc_atom, ONLY: calculate_vxc_atom
119 USE task_list_types, ONLY: task_list_type
120#include "./base/base_uses.f90"
121
122 IMPLICIT NONE
123 PRIVATE
124
125 PUBLIC :: admm_mo_calc_rho_aux, &
126 admm_mo_calc_rho_aux_kp, &
127 admm_mo_merge_ks_matrix, &
128 admm_mo_merge_derivs, &
129 admm_aux_response_density, &
130 calc_mixed_overlap_force, &
131 scale_dm, &
132 admm_fit_mo_coeffs, &
133 admm_update_ks_atom, &
134 calc_admm_mo_derivatives, &
135 calc_admm_ovlp_forces, &
136 calc_admm_ovlp_forces_kp, &
137 admm_projection_derivative, &
138 kpoint_calc_admm_matrices
139
140 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'admm_methods'
141
142CONTAINS
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 :: 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)
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_uplo_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_uplo_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_uplo_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_uplo_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_uplo_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 :: 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_uplo_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_uplo_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)
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_uplo_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_uplo_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_uplo_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_uplo_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 :: 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)
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_uplo_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_uplo_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)
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_uplo_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_uplo_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_uplo_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 cp_cfm_uplo_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 :: 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)
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_uplo_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_uplo_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 cp_cfm_uplo_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
4129END MODULE admm_methods
cp_cfm_basic_linalg::cp_cfm_scale
Definition cp_cfm_basic_linalg.F:60
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:82
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
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:21
bibliography::merlot2014
integer, save, public merlot2014
Definition bibliography.F:36
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:240
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:351
cp_cfm_basic_linalg::cp_cfm_uplo_to_full
subroutine, public cp_cfm_uplo_to_full(matrix, workspace, uplo)
...
Definition cp_cfm_basic_linalg.F:1347
cp_cfm_cholesky
various cholesky decomposition related routines
Definition cp_cfm_cholesky.F:13
cp_cfm_cholesky::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_cholesky.F:51
cp_cfm_cholesky::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_cholesky.F:112
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:155
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:813
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:761
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_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1165
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:233
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
Definition cp_dbcsr_api.F:956
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:930
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1027
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:517
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:679
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:989
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
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:214
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:695
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:111
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:1885
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:713
cp_fm_basic_linalg::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
...
Definition cp_fm_basic_linalg.F:3189
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_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1747
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:1665
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_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:434
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:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
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:1567
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:1945
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:1009
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:528
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:1881
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:693
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:164
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:853
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:1063
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:854
input_constants::do_admm_purify_none
integer, parameter, public do_admm_purify_none
Definition input_constants.F:854
input_constants::do_admm_purify_cauchy_subspace
integer, parameter, public do_admm_purify_cauchy_subspace
Definition input_constants.F:854
input_constants::do_admm_purify_cauchy
integer, parameter, public do_admm_purify_cauchy
Definition input_constants.F:854
input_constants::do_admm_purify_mo_diag
integer, parameter, public do_admm_purify_mo_diag
Definition input_constants.F:854
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:1047
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:809
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:1143
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:1019
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:779
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:365
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::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
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:1711
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_pp, 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, harris_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, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, pw_env_sub)
...
Definition qs_gapw_densities.F:55
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:85
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
admm_types::admm_gapw_r3d_rs_type
A subtype of the admm_env that contains the extra data needed for an ADMM GAPW calculation.
Definition admm_types.F:85
admm_types::admm_type
stores some data used in wavefunction fitting
Definition admm_types.F:120
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:68
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:186
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::copy_info_type
Stores the state of a copy between cp_fm_start_copy_general and cp_fm_finish_copy_general.
Definition cp_fm_types.F:139
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_env_type
Keeps information about a specific k-point.
Definition kpoint_types.F:73
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:150
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:136
qs_local_rho_types::local_rho_type
Definition qs_local_rho_types.F:43
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:97
task_list_types::task_list_type
Definition task_list_types.F:58