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