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