(git:37370f3)
Loading...
Searching...
No Matches
dft_plus_u.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!> \brief Add the DFT+U contribution to the Hamiltonian matrix
9!> \details The implemented methods refers to:\n
10!> S. L. Dudarev, D. Nguyen Manh, and A. P. Sutton,
11!> Philos. Mag. B \b 75, 613 (1997)\n
12!> S. L. Dudarev et al.,
13!> Phys. Rev. B \b 57, 1505 (1998)
14!> \author Matthias Krack (MK)
15!> \date 14.01.2008
16!> \version 1.0
17! **************************************************************************************************
18MODULE dft_plus_u
19 USE omp_lib, ONLY: omp_get_max_threads,&
20 omp_get_thread_num
21 USE atomic_kind_types, ONLY: atomic_kind_type,&
22 get_atomic_kind,&
23 get_atomic_kind_set
24 USE basis_set_types, ONLY: get_gto_basis_set,&
25 gto_basis_set_type
26 USE bibliography, ONLY: dudarev1997,&
27 dudarev1998,&
28 cite_reference
29 USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
30 cp_cfm_transpose
31 USE cp_cfm_diag, ONLY: cp_cfm_heevd
32 USE cp_cfm_types, ONLY: cp_cfm_create,&
33 cp_cfm_get_info,&
34 cp_cfm_release,&
35 cp_cfm_to_cfm,&
36 cp_cfm_type,&
37 cp_fm_to_cfm
38 USE cp_control_types, ONLY: dft_control_type
39 USE cp_dbcsr_api, ONLY: &
40 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_finalize, &
41 dbcsr_get_block_p, dbcsr_init_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
42 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
43 dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_set, dbcsr_type
44 USE cp_dbcsr_contrib, ONLY: dbcsr_get_block_diag,&
45 dbcsr_get_diag,&
46 dbcsr_set_diag
47 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
48 copy_fm_to_dbcsr,&
49 cp_dbcsr_plus_fm_fm_t,&
50 cp_dbcsr_sm_fm_multiply
51 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix,&
52 write_fm_with_basis_info
53 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
54 cp_fm_schur_product,&
55 cp_fm_transpose
56 USE cp_fm_diag, ONLY: choose_eigv_solver
57 USE cp_fm_struct, ONLY: cp_fm_struct_type
58 USE cp_fm_types, ONLY: cp_fm_create,&
59 cp_fm_get_info,&
60 cp_fm_release,&
61 cp_fm_set_submatrix,&
62 cp_fm_to_fm,&
63 cp_fm_type
64 USE cp_log_handling, ONLY: cp_get_default_logger,&
65 cp_logger_type
66 USE cp_output_handling, ONLY: cp_p_file,&
67 cp_print_key_finished_output,&
68 cp_print_key_should_output,&
69 cp_print_key_unit_nr,&
70 low_print_level
71 USE input_constants, ONLY: diag_update_method_adiis,&
72 plus_u_lowdin,&
73 plus_u_mulliken,&
74 plus_u_mulliken_charges,&
75 plus_u_tensorial
76 USE input_section_types, ONLY: section_vals_type
77 USE kinds, ONLY: default_string_length,&
78 dp
79 USE kpoint_types, ONLY: kpoint_env_type,&
80 kpoint_type
81 USE local_gemm_api, ONLY: local_gemm_pu_gpu,&
82 local_gemm_ctxt_type
83 USE mathconstants, ONLY: z_one,&
84 z_zero
85 USE mathlib, ONLY: invert_matrix,&
86 jacobi
87 USE message_passing, ONLY: mp_para_env_type
88 USE orbital_symbols, ONLY: sgf_symbol
89 USE parallel_gemm_api, ONLY: parallel_gemm
90 USE particle_methods, ONLY: get_particle_set
91 USE particle_types, ONLY: particle_type
92 USE physcon, ONLY: evolt
93 USE qs_energy_types, ONLY: qs_energy_type
94 USE qs_environment_types, ONLY: get_qs_env,&
95 qs_environment_type
96 USE qs_force_types, ONLY: qs_force_type
97 USE qs_kind_types, ONLY: get_qs_kind,&
98 get_qs_kind_set,&
99 qs_kind_type,&
100 set_qs_kind
101 USE qs_kpoint_operators, ONLY: kpoint_operator_context_create,&
102 kpoint_operator_context_release,&
103 kpoint_operator_context_type,&
104 kpoint_operator_get
105 USE qs_rho_types, ONLY: qs_rho_get,&
106 qs_rho_type
107 USE qs_scf_types, ONLY: qs_scf_env_type
108#include "./base/base_uses.f90"
109
110 IMPLICIT NONE
111
112 PRIVATE
113
114 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dft_plus_u'
115
116 PUBLIC :: plus_u
117 PUBLIC :: lowdin_kp_init, lowdin_kp_initial_density, lowdin_kp_population, lowdin_kp_apply
118 PUBLIC :: lowdin_kp_forces
119 PUBLIC :: lowdin_kp_apply_mo, lowdin_kp_kernel
120 PUBLIC :: lowdin_kp_prepare
121
122CONTAINS
123! **************************************************************************************************
124!> \brief Add the DFT+U contribution to the Hamiltonian matrix.\n
125!> Wrapper routine for all "+U" methods
126!> \param[in] qs_env Quickstep environment
127!> \param[in,out] matrix_h Hamiltonian matrices for each spin
128!> \param[in,out] matrix_w Energy weighted density matrices for each spin
129!> \date 14.01.2008
130!> \author Matthias Krack (MK)
131!> \version 1.0
132! **************************************************************************************************
133 SUBROUTINE plus_u(qs_env, matrix_h, matrix_w)
134
135 TYPE(qs_environment_type), POINTER :: qs_env
136 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
137 POINTER :: matrix_h, matrix_w
138
139 CHARACTER(LEN=*), PARAMETER :: routinen = 'plus_u'
140
141 INTEGER :: handle, output_unit, print_level
142 LOGICAL :: orthonormal_basis, should_output
143 TYPE(cp_logger_type), POINTER :: logger
144 TYPE(dft_control_type), POINTER :: dft_control
145 TYPE(section_vals_type), POINTER :: input
146
147 CALL timeset(routinen, handle)
148
149 cpassert(ASSOCIATED(qs_env))
150
151 NULLIFY (input, dft_control)
152
153 logger => cp_get_default_logger()
154
155 CALL get_qs_env(qs_env=qs_env, &
156 input=input, &
157 dft_control=dft_control)
158
159 CALL cite_reference(dudarev1997)
160 CALL cite_reference(dudarev1998)
161
162 ! Later we could save here some time, if the method in use has this property
163 ! which then has to be figured out here
164
165 orthonormal_basis = .false.
166
167 ! Setup print control
168
169 print_level = logger%iter_info%print_level
170 should_output = (btest(cp_print_key_should_output(logger%iter_info, input, &
171 "DFT%PRINT%PLUS_U"), cp_p_file) .AND. &
172 (.NOT. PRESENT(matrix_w)))
173 output_unit = cp_print_key_unit_nr(logger, input, "DFT%PRINT%PLUS_U", &
174 extension=".plus_u", &
175 ignore_should_output=should_output, &
176 log_filename=.false.)
177
178 ! Select DFT+U method
179
180 SELECT CASE (dft_control%plus_u_method_id)
181 CASE (plus_u_lowdin)
182 IF (orthonormal_basis) THEN
183 ! For an orthonormal basis the Lowdin method and the Mulliken method
184 ! are equivalent
185 CALL mulliken(qs_env, orthonormal_basis, matrix_h, &
186 should_output, output_unit, print_level)
187 ELSE
188 CALL lowdin(qs_env, matrix_h, matrix_w, &
189 should_output, output_unit, print_level)
190 END IF
191 CASE (plus_u_mulliken)
192 CALL mulliken(qs_env, orthonormal_basis, matrix_h, &
193 should_output, output_unit, print_level)
194 CASE (plus_u_mulliken_charges)
195 CALL mulliken_charges(qs_env, orthonormal_basis, matrix_h, matrix_w, &
196 should_output, output_unit, print_level)
197 CASE (plus_u_tensorial)
198 CALL tensorial(qs_env, matrix_h, matrix_w, &
199 should_output, output_unit, print_level)
200 CASE DEFAULT
201 cpabort("Invalid DFT+U method requested")
202 END SELECT
203
204 CALL cp_print_key_finished_output(output_unit, logger, input, "DFT%PRINT%PLUS_U", &
205 ignore_should_output=should_output)
206
207 CALL timestop(handle)
208
209 END SUBROUTINE plus_u
210
211! **************************************************************************************************
212!> \brief Add a DFT+U contribution to the Hamiltonian matrix\n
213!> using a method based on Lowdin charges
214!> \f[Q = S^{1/2} P S^{1/2}\f]
215!> where \b P and \b S are the density and the
216!> overlap matrix, respectively.
217!> \param[in] qs_env Quickstep environment
218!> \param[in,out] matrix_h Hamiltonian matrices for each spin
219!> \param[in,out] matrix_w Energy weighted density matrices for each spin
220!> \param should_output ...
221!> \param output_unit ...
222!> \param print_level ...
223!> \date 02.07.2008
224!> \par
225!> \f{eqnarray*}{
226!> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}
227!> & = & E^{\rm DFT} + \frac{1}{2}(U - J)\sum_\mu (q_\mu - q_\mu^2)\\‍[1ex]
228!> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
229!> & = & \frac{\partial E^{\rm DFT}}
230!> {\partial P_{\mu\nu}} +
231!> \frac{\partial E^{\rm U}}
232!> {\partial P_{\mu\nu}}\\\
233!> & = & H_{\mu\nu} +
234!> \frac{\partial E^{\rm U}}{\partial q_\mu}
235!> \frac{\partial q_\mu}{\partial P_{\mu\nu}}\\\
236!> \f}
237!> \author Matthias Krack (MK)
238!> \version 1.0
239! **************************************************************************************************
240 SUBROUTINE lowdin(qs_env, matrix_h, matrix_w, should_output, output_unit, &
241 print_level)
242
243 TYPE(qs_environment_type), POINTER :: qs_env
244 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
245 POINTER :: matrix_h, matrix_w
246 LOGICAL, INTENT(IN) :: should_output
247 INTEGER, INTENT(IN) :: output_unit, print_level
248
249 CHARACTER(LEN=*), PARAMETER :: routinen = 'lowdin'
250
251 CHARACTER(LEN=10) :: spin_info
252 CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: symbol
253 CHARACTER(LEN=default_string_length) :: atomic_kind_name
254 INTEGER :: atom_a, handle, i, i0, iatom, ikind, iorb, isb, iset, isgf, ishell, ispin, j, &
255 jsb, jsgf, lu, m, max_scf, n, natom_of_kind, nkind, norb, nsb, nsbsize, nset, nsgf, &
256 nsgf_kind, nspin
257 INTEGER, ALLOCATABLE, DIMENSION(:) :: shell_sgf
258 INTEGER, DIMENSION(1) :: iloc
259 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals
260 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
261 LOGICAL :: debug, dft_plus_u_atom, do_kpoints, &
262 found, if_plus_j, just_energy, smear
263 LOGICAL, ALLOCATABLE, DIMENSION(:) :: orb_occ
264 REAL(kind=dp) :: eps_scf, eps_u_ramping, fspin, hund_j, occ, perturbation_strength, sij, &
265 trq, trq2, trqxq_hund_j, u_minus_j, u_minus_j_target, u_ramping
266 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigval, q_eigval
267 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: orbq, q_eigvec, q_matrix, q_opp_matrix, &
268 q_work, slam
269 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
270 POINTER :: local_data
271 REAL(kind=dp), DIMENSION(:, :), POINTER :: q_block, q_opp_block, v_block, vhxc_block
272 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
273 TYPE(cp_fm_struct_type), POINTER :: fmstruct
274 TYPE(cp_fm_type) :: fm_sev, fm_work1, fm_work2, slambda
275 TYPE(cp_fm_type), POINTER :: fm_s_half
276 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vhxc
277 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
278 TYPE(dbcsr_type) :: sm_q, sm_q_opp, sm_v
279 TYPE(dbcsr_type), POINTER :: sm_h, sm_p, sm_p_opp, sm_s, sm_vhxc, sm_w
280 TYPE(dft_control_type), POINTER :: dft_control
281 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
282 TYPE(kpoint_type), POINTER :: kpoints
283 TYPE(mp_para_env_type), POINTER :: para_env
284 TYPE(qs_energy_type), POINTER :: energy
285 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
286 TYPE(qs_rho_type), POINTER :: rho
287 TYPE(qs_scf_env_type), POINTER :: scf_env
288
289 CALL timeset(routinen, handle)
290
291 debug = .false. ! Set to .TRUE. to print debug information
292
293 NULLIFY (sm_h, sm_p, sm_p_opp, sm_s, sm_w, vhxc_block, sm_vhxc, matrix_vhxc)
294
295 smear = .false.
296 max_scf = -1
297 eps_scf = 1.0e30_dp
298 if_plus_j = .false.
299
300 CALL get_qs_env(qs_env=qs_env, &
301 atomic_kind_set=atomic_kind_set, &
302 qs_kind_set=qs_kind_set, &
303 dft_control=dft_control, &
304 do_kpoints=do_kpoints, &
305 kpoints=kpoints, &
306 energy=energy, &
307 matrix_s_kp=matrix_s, &
308 matrix_vhxc=matrix_vhxc, &
309 rho=rho, &
310 scf_env=scf_env, &
311 para_env=para_env)
312
313 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
314
315 energy%dft_plus_u = 0.0_dp
316
317 nspin = dft_control%nspins
318 IF (dft_control%mtlr_dft_with_perturbation) THEN
319 IF (.NOT. ASSOCIATED(matrix_vhxc)) THEN
320 cpabort("MTLR requires the projected Hxc matrix for every spin channel.")
321 END IF
322 IF (SIZE(matrix_vhxc) /= nspin) THEN
323 cpabort("The number of projected Hxc matrices does not match the spin channels.")
324 END IF
325 END IF
326
327 IF (nspin == 2) THEN
328 fspin = 1.0_dp
329 ELSE
330 fspin = 0.5_dp
331 END IF
332
333 ! Get the number of basis functions and atomic kinds.
334 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
335
336 nkind = SIZE(atomic_kind_set)
337
338 IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
339 just_energy = .false.
340 ELSE
341 just_energy = .true.
342 END IF
343
344 ! Check if unlike spin +J correction is needed and reject an
345 ! incompatible enforced occupation setup before any expensive work.
346 DO ikind = 1, nkind
347 NULLIFY (orbitals)
348 CALL get_qs_kind(qs_kind_set(ikind), &
349 hund_j=hund_j, &
350 orbitals=orbitals)
351 IF (hund_j /= 0.0_dp) THEN
352 if_plus_j = .true.
353 IF (ASSOCIATED(orbitals)) THEN
354 CALL cp_abort(__location__, "ENFORCE_OCCUPATION is incompatible with a nonzero Hund J "// &
355 "in the Lowdin DFT+U+J method.")
356 END IF
357 END IF
358 END DO
359
360 NULLIFY (orbitals)
361
362 IF (.NOT. do_kpoints) THEN
363 ! Retrieve S^(1/2) from the SCF environment
364 fm_s_half => scf_env%s_half
365 cpassert(ASSOCIATED(fm_s_half))
366 ! work matrices
367 CALL cp_fm_get_info(fm_s_half, matrix_struct=fmstruct)
368 CALL cp_fm_create(matrix=fm_work1, matrix_struct=fmstruct, &
369 name="FULL WORK MATRIX 1")
370 CALL cp_fm_create(matrix=fm_work2, matrix_struct=fmstruct, &
371 name="FULL WORK MATRIX 2")
372 END IF
373
374 ! Calculate S eigenvectors and Lambda matrix for forces
375 ! See sTDA forces (get_lowdin_mo_coefficients in qs_tddfpt2_stda_utils
376 ! A. Hehn et al JCTC 2022, 18, 4186
377 IF (PRESENT(matrix_w) .AND. .NOT. do_kpoints) THEN
378 CALL cp_fm_create(matrix=fm_sev, matrix_struct=fmstruct)
379 CALL cp_fm_create(matrix=slambda, matrix_struct=fmstruct)
380 ALLOCATE (eigval(nsgf), slam(nsgf, 1))
381 sm_s => matrix_s(1, 1)%matrix
382 CALL copy_dbcsr_to_fm(sm_s, fm_work1)
383 CALL choose_eigv_solver(fm_work1, fm_sev, eigval)
384 !
385 DO i = 1, nsgf
386 IF (eigval(i) > 0._dp) THEN
387 slam(i, 1) = sqrt(eigval(i))
388 ELSE
389 cpabort("S matrix not positive definit")
390 END IF
391 END DO
392 DO i = 1, nsgf
393 CALL cp_fm_set_submatrix(slambda, slam, 1, i, nsgf, 1, 1.0_dp, 0.0_dp)
394 END DO
395 DO i = 1, nsgf
396 CALL cp_fm_set_submatrix(slambda, slam, i, 1, 1, nsgf, 1.0_dp, 1.0_dp, .true.)
397 END DO
398 CALL cp_fm_get_info(slambda, local_data=local_data)
399 DO i = 1, SIZE(local_data, 2)
400 DO j = 1, SIZE(local_data, 1)
401 sij = local_data(j, i)
402 IF (sij > 0.0_dp) sij = 1.0_dp/sij
403 local_data(j, i) = sij
404 END DO
405 END DO
406 DEALLOCATE (eigval, slam)
407 END IF
408
409 ! Calculate S^(1/2)*P*S^(1/2)
410 IF (do_kpoints) THEN
411 ! OT and ADIIS evaluate the current trial density, including energy-only line searches.
412 ! Density mixing instead owns separate accepted coordinates.
413 IF (qs_env%scf_control%use_ot .OR. &
414 qs_env%scf_control%diagonalization%update_method == diag_update_method_adiis .OR. PRESENT(matrix_w)) THEN
415 CALL lowdin_kp_population(kpoints, alpha=1.0_dp)
416 END IF
417 ALLOCATE (orbq(nsgf, nspin))
418 orbq = 0.0_dp
419 orbq(kpoints%lowdin_orbitals, :) = kpoints%lowdin_q
420 END IF
421
422 ! Create local block diagonal matrices
423 sm_s => matrix_s(1, 1)%matrix
424 CALL dbcsr_get_block_diag(sm_s, sm_q)
425 IF (if_plus_j .AND. nspin == 2) THEN
426 CALL dbcsr_get_block_diag(sm_s, sm_q_opp)
427 END IF
428 CALL dbcsr_get_block_diag(sm_s, sm_v)
429
430 ! Loop over all spins
431 DO ispin = 1, nspin
432
433 CALL dbcsr_set(sm_q, 0.0_dp)
434 IF (if_plus_j .AND. nspin == 2) THEN
435 CALL dbcsr_set(sm_q_opp, 0.0_dp)
436 END IF
437 CALL dbcsr_set(sm_v, 0.0_dp)
438
439 IF (dft_control%mtlr_dft_with_perturbation) THEN
440 IF (ispin == 1) perturbation_strength = 0.9_dp*dft_control%perturbation_strength
441 IF (ispin == 2) perturbation_strength = 1.1_dp*dft_control%perturbation_strength
442 IF (.NOT. ASSOCIATED(matrix_vhxc(ispin)%matrix)) THEN
443 cpabort("MTLR projected Hxc matrix is not initialized.")
444 END IF
445 sm_vhxc => matrix_vhxc(ispin)%matrix
446 END IF
447
448 IF (do_kpoints) THEN
449 CALL dbcsr_set_diag(sm_q, orbq(:, ispin))
450 IF (if_plus_j .AND. nspin == 2) CALL dbcsr_set_diag(sm_q_opp, orbq(:, 3 - ispin))
451 ELSE
452 ! Calculate S^(1/2)*P*S^(1/2) as a full matrix (Lowdin)
453 sm_p => matrix_p(ispin, 1)%matrix
454 CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_s_half, fm_work1, nsgf)
455 CALL parallel_gemm(transa="N", &
456 transb="N", &
457 m=nsgf, &
458 n=nsgf, &
459 k=nsgf, &
460 alpha=1.0_dp, &
461 matrix_a=fm_s_half, &
462 matrix_b=fm_work1, &
463 beta=0.0_dp, &
464 matrix_c=fm_work2)
465 IF (debug) THEN
466 CALL cp_dbcsr_write_sparse_matrix(sm_p, 4, 6, qs_env, para_env, &
467 output_unit=output_unit)
468 CALL write_fm_with_basis_info(fm_s_half, 4, 6, qs_env, para_env, &
469 output_unit=output_unit)
470 CALL write_fm_with_basis_info(fm_work2, 4, 6, qs_env, para_env, &
471 output_unit=output_unit)
472 END IF ! debug
473 ! Copy occupation matrix to sparse matrix format, finally we are only
474 ! interested in the diagonal (atomic) blocks, i.e. the previous full
475 ! matrix product is not the most efficient choice, anyway.
476 CALL copy_fm_to_dbcsr(fm_work2, sm_q, keep_sparsity=.true.)
477
478 ! If nspin == 1 and one would like to add "unlike" spin +J corrections,
479 ! sm_p is already enough and sm_p_opp is not needed
480 IF (if_plus_j .AND. nspin == 2) THEN
481 ! Calculate S^(1/2)*P*S^(1/2) for the other spin channel if
482 ! unlike spin +J correction is activated.
483 sm_p_opp => matrix_p(3 - ispin, 1)%matrix
484 CALL cp_dbcsr_sm_fm_multiply(sm_p_opp, fm_s_half, fm_work1, nsgf)
485 CALL parallel_gemm(transa="N", &
486 transb="N", &
487 m=nsgf, &
488 n=nsgf, &
489 k=nsgf, &
490 alpha=1.0_dp, &
491 matrix_a=fm_s_half, &
492 matrix_b=fm_work1, &
493 beta=0.0_dp, &
494 matrix_c=fm_work2)
495 CALL copy_fm_to_dbcsr(fm_work2, sm_q_opp, keep_sparsity=.true.)
496 END IF
497 END IF
498
499 ! E[DFT+U] = E[DFT] + E[U]
500 ! = E[DFT] + (U - J)*(Tr(q) - Tr(q*q))/2
501
502 ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U]
503 ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j)
504 ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j))
505
506 ! Loop over all atomic kinds
507 DO ikind = 1, nkind
508
509 ! Load the required atomic kind data
510 CALL get_atomic_kind(atomic_kind_set(ikind), &
511 atom_list=atom_list, &
512 name=atomic_kind_name, &
513 natom=natom_of_kind)
514
515 CALL get_qs_kind(qs_kind_set(ikind), &
516 dft_plus_u_atom=dft_plus_u_atom, &
517 l_of_dft_plus_u=lu, &
518 nsgf=nsgf_kind, &
519 basis_set=orb_basis_set, &
520 u_minus_j=u_minus_j, &
521 hund_j=hund_j, &
522 u_minus_j_target=u_minus_j_target, &
523 u_ramping=u_ramping, &
524 eps_u_ramping=eps_u_ramping, &
525 orbitals=orbitals, &
526 eps_scf=eps_scf, &
527 max_scf=max_scf, &
528 smear=smear)
529
530 ! Check, if the atoms of this atomic kind need a DFT+U correction
531 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
532 IF (.NOT. dft_plus_u_atom) cycle
533 IF (lu < 0) cycle
534
535 ! Apply U ramping if requested
536 IF ((ispin == 1) .AND. (u_ramping > 0.0_dp) .AND. .NOT. do_kpoints) THEN
537 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
538 u_minus_j = min(u_minus_j + u_ramping, u_minus_j_target)
539 CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
540 END IF
541 IF (should_output .AND. (output_unit > 0)) THEN
542 WRITE (unit=output_unit, fmt="(T3,A,3X,A,F0.3,A)") &
543 "Kind name: "//trim(adjustl(atomic_kind_name)), &
544 "U(eff) = ", u_minus_j*evolt, " eV"
545 END IF
546 END IF
547
548 IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) cycle
549
550 ! Load the required Gaussian basis set data
551 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
552 first_sgf=first_sgf, &
553 l=l, &
554 last_sgf=last_sgf, &
555 nset=nset, &
556 nshell=nshell)
557
558 ! Select this kind's shell once; occupations and potentials use the same AO order.
559 ALLOCATE (shell_sgf(nsgf_kind))
560 n = 0
561 DO iset = 1, nset
562 DO ishell = 1, nshell(iset)
563 IF (l(ishell, iset) /= lu) cycle
564 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
565 n = n + 1
566 shell_sgf(n) = isgf
567 END DO
568 END DO
569 END DO
570 nsbsize = 2*lu + 1
571 nsb = n/nsbsize
572
573 ALLOCATE (q_matrix(n, n))
574 q_matrix(:, :) = 0.0_dp
575
576 IF (hund_j /= 0.0_dp .AND. nspin == 2) THEN
577 ALLOCATE (q_opp_matrix(n, n))
578 q_opp_matrix(:, :) = 0.0_dp
579 END IF
580
581 IF (dft_control%mtlr_dft_with_perturbation .AND. &
582 dft_control%mtlr_ikind == ikind) THEN
583 dft_control%vhxc(ispin) = 0.0_dp
584 dft_control%trq(ispin) = 0.0_dp
585 END IF
586
587 ! Print headline if requested
588 IF (should_output .AND. (print_level > low_print_level)) THEN
589 IF (output_unit > 0) THEN
590 ALLOCATE (symbol(nsbsize))
591 DO m = -lu, lu
592 symbol(lu + m + 1) = sgf_symbol(0, lu, m)
593 END DO
594 IF (nspin > 1) THEN
595 WRITE (unit=spin_info, fmt="(A8,I2)") " of spin", ispin
596 ELSE
597 spin_info = ""
598 END IF
599 WRITE (unit=output_unit, fmt="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
600 "DFT+U occupations"//trim(spin_info)//" for the atoms of atomic kind ", ikind, &
601 ": "//trim(atomic_kind_name), &
602 "Atom Shell ", (adjustr(symbol(i)), i=1, nsbsize), " Trace"
603 DEALLOCATE (symbol)
604 END IF
605 END IF
606
607 ! Loop over all atoms of the current atomic kind
608 DO iatom = 1, natom_of_kind
609 atom_a = atom_list(iatom)
610 q_matrix(:, :) = 0.0_dp
611 IF (ALLOCATED(q_opp_matrix)) THEN
612 q_opp_matrix(:, :) = 0.0_dp
613 END IF
614
615 ! Get diagonal block
616 CALL dbcsr_get_block_p(matrix=sm_q, &
617 row=atom_a, &
618 col=atom_a, &
619 block=q_block, &
620 found=found)
621 IF (hund_j /= 0.0_dp .AND. nspin == 2) THEN
622 CALL dbcsr_get_block_p(matrix=sm_q_opp, &
623 row=atom_a, &
624 col=atom_a, &
625 block=q_opp_block, &
626 found=found)
627 END IF
628 IF (dft_control%mtlr_dft_with_perturbation .AND. &
629 dft_control%mtlr_ikind == ikind .AND. &
630 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
631 CALL dbcsr_get_block_p(matrix=sm_vhxc, &
632 row=atom_a, &
633 col=atom_a, &
634 block=vhxc_block, &
635 found=found)
636 END IF
637
638 IF (ASSOCIATED(q_block)) THEN
639 ! Calculate energy contribution to E(U)
640 DO i = 1, n
641 isgf = shell_sgf(i)
642 q_matrix(i, i) = q_block(isgf, isgf)
643 IF (hund_j /= 0.0_dp .AND. nspin == 2 .AND. ASSOCIATED(q_opp_block)) THEN
644 q_opp_matrix(i, i) = q_opp_block(isgf, isgf)
645 END IF
646 IF (ASSOCIATED(vhxc_block) .AND. dft_control%mtlr_dft_with_perturbation .AND. &
647 dft_control%mtlr_ikind == ikind .AND. qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
648 dft_control%vhxc(ispin) = dft_control%vhxc(ispin) + vhxc_block(isgf, isgf)
649 dft_control%trq(ispin) = dft_control%trq(ispin) + q_block(isgf, isgf)
650 END IF
651 END DO
652
653 ! Perform the requested manipulations of the (initial) orbital occupations
654 IF (ASSOCIATED(orbitals)) THEN
655 IF ((qs_env%scf_env%iter_delta >= eps_scf) .OR. &
656 ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
657 (qs_env%scf_env%iter_count <= max_scf))) THEN
658 ALLOCATE (orb_occ(nsbsize))
659 ALLOCATE (q_eigval(n))
660 q_eigval(:) = 0.0_dp
661 ALLOCATE (q_eigvec(n, n))
662 q_eigvec(:, :) = 0.0_dp
663 norb = SIZE(orbitals)
664 CALL jacobi(q_matrix, q_eigval, q_eigvec)
665 q_matrix(:, :) = 0.0_dp
666 DO isb = 1, nsb
667 trq = 0.0_dp
668 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
669 trq = trq + q_eigval(i)
670 END DO
671 IF (smear) THEN
672 occ = trq/real(norb, kind=dp)
673 ELSE
674 occ = 1.0_dp/fspin
675 END IF
676 orb_occ(:) = .false.
677 iloc = maxloc(q_eigvec(:, isb*nsbsize))
678 jsb = int((iloc(1) - 1)/nsbsize) + 1
679 i = 0
680 i0 = (jsb - 1)*nsbsize + 1
681 iorb = -1000
682 DO j = i0, jsb*nsbsize
683 i = i + 1
684 IF (i > norb) THEN
685 DO m = -lu, lu
686 IF (.NOT. orb_occ(lu + m + 1)) THEN
687 iorb = i0 + lu + m
688 orb_occ(lu + m + 1) = .true.
689 END IF
690 END DO
691 ELSE
692 iorb = i0 + lu + orbitals(i)
693 orb_occ(lu + orbitals(i) + 1) = .true.
694 END IF
695 cpassert(iorb /= -1000)
696 iloc = maxloc(q_eigvec(iorb, :))
697 q_eigval(iloc(1)) = min(occ, trq)
698 q_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1)) ! backtransform left
699 trq = trq - q_eigval(iloc(1))
700 END DO
701 END DO
702 q_matrix(:, :) = matmul(q_matrix, transpose(q_eigvec)) ! backtransform right
703 DEALLOCATE (orb_occ)
704 DEALLOCATE (q_eigval)
705 DEALLOCATE (q_eigvec)
706 END IF
707 END IF ! orbitals associated
708
709 trq = fspin*sum([(q_matrix(i, i), i=1, n)])
710 trq2 = fspin*fspin*sum(q_matrix*transpose(q_matrix))
711 trqxq_hund_j = 0.0_dp
712 IF (hund_j /= 0.0_dp) THEN
713 IF (nspin == 2) THEN
714 trqxq_hund_j = sum(q_matrix*transpose(q_opp_matrix))
715 ELSE
716 trqxq_hund_j = sum(q_matrix*transpose(q_matrix))
717 END IF
718 END IF
719
720 ! Calculate energy contribution to E(U)
721 energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
722 IF (hund_j /= 0.0_dp) THEN
723 trqxq_hund_j = fspin*fspin*trqxq_hund_j
724 energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*hund_j*trqxq_hund_j/fspin
725 END IF
726
727 ! Calculate potential V(U) = dE(U)/dq
728 IF (.NOT. just_energy) THEN
729 CALL dbcsr_get_block_p(matrix=sm_v, &
730 row=atom_a, &
731 col=atom_a, &
732 block=v_block, &
733 found=found)
734 cpassert(ASSOCIATED(v_block))
735
736 DO j = 1, n
737 jsgf = shell_sgf(j)
738 DO i = 1, n
739 isgf = shell_sgf(i)
740 IF (i == j) THEN
741 v_block(isgf, isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i, i))
742 IF (dft_control%mtlr_dft_with_perturbation .AND. &
743 dft_control%mtlr_ikind == ikind .AND. &
744 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
745 v_block(isgf, isgf) = v_block(isgf, isgf) + perturbation_strength
746 END IF
747 ELSE
748 cpassert(abs(q_matrix(j, i)) < 1.0e-14_dp)
749 IF (hund_j /= 0.0_dp .AND. nspin == 2) THEN
750 cpassert(abs(q_opp_matrix(j, i)) < 1.0e-14_dp)
751 END IF
752 v_block(isgf, jsgf) = -u_minus_j*fspin*q_matrix(j, i)
753 END IF
754 IF (hund_j /= 0.0_dp) THEN
755 IF (nspin == 2) THEN
756 v_block(isgf, jsgf) = v_block(isgf, jsgf) + hund_j*fspin*q_opp_matrix(j, i)
757 ELSE
758 v_block(isgf, jsgf) = v_block(isgf, jsgf) + hund_j*fspin*q_matrix(j, i)
759 END IF
760 END IF
761 END DO
762 END DO
763 END IF ! not just energy
764
765 END IF ! q_block associated
766
767 ! Consider print requests
768 IF (should_output .AND. (print_level > low_print_level)) THEN
769 CALL para_env%sum(q_matrix)
770 IF (output_unit > 0) THEN
771 ALLOCATE (q_work(nsb, nsbsize))
772 q_work(:, :) = 0.0_dp
773 DO isb = 1, nsb
774 j = 0
775 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
776 j = j + 1
777 q_work(isb, j) = q_matrix(i, i)
778 END DO
779 END DO
780 DO isb = 1, nsb
781 WRITE (unit=output_unit, fmt="(T3,I6,2X,I6,2X,10F8.3)") &
782 atom_a, isb, q_work(isb, :), sum(q_work(isb, :))
783 END DO
784 WRITE (unit=output_unit, fmt="(T12,A,2X,10F8.3)") &
785 "Total", (sum(q_work(:, i)), i=1, nsbsize), sum(q_work)
786 WRITE (unit=output_unit, fmt="(A)") ""
787 DEALLOCATE (q_work)
788 IF (debug) THEN
789 ! Print the DFT+U occupation matrix
790 WRITE (unit=output_unit, fmt="(T9,70I10)") (i, i=1, n)
791 DO i = 1, n
792 WRITE (unit=output_unit, fmt="(T3,I6,70F10.6)") i, q_matrix(i, :)
793 END DO
794 ! Print the eigenvalues and eigenvectors of the occupation matrix
795 ALLOCATE (q_eigval(n))
796 q_eigval(:) = 0.0_dp
797 ALLOCATE (q_eigvec(n, n))
798 q_eigvec(:, :) = 0.0_dp
799 CALL jacobi(q_matrix, q_eigval, q_eigvec)
800 WRITE (unit=output_unit, fmt="(/,T9,70I10)") (i, i=1, n)
801 WRITE (unit=output_unit, fmt="(T9,71F10.6)") (q_eigval(i), i=1, n), &
802 sum(q_eigval(1:n))
803 DO i = 1, n
804 WRITE (unit=output_unit, fmt="(T3,I6,70F10.6)") i, q_eigvec(i, :)
805 END DO
806 DEALLOCATE (q_eigval)
807 DEALLOCATE (q_eigvec)
808 END IF ! debug
809 END IF
810 IF (debug) THEN
811 ! Print the full atomic occupation matrix block
812 ALLOCATE (q_work(nsgf_kind, nsgf_kind))
813 q_work(:, :) = 0.0_dp
814 IF (ASSOCIATED(q_block)) q_work(:, :) = q_block(:, :)
815 CALL para_env%sum(q_work)
816 IF (output_unit > 0) THEN
817 norb = SIZE(q_work, 1)
818 WRITE (unit=output_unit, fmt="(/,T9,200I10)") (i, i=1, norb)
819 DO i = 1, norb
820 WRITE (unit=output_unit, fmt="(T3,I6,200F10.6)") i, q_work(i, :)
821 END DO
822 ALLOCATE (q_eigval(norb))
823 q_eigval(:) = 0.0_dp
824 ALLOCATE (q_eigvec(norb, norb))
825 q_eigvec(:, :) = 0.0_dp
826 CALL jacobi(q_work, q_eigval, q_eigvec)
827 WRITE (unit=output_unit, fmt="(/,T9,200I10)") (i, i=1, norb)
828 WRITE (unit=output_unit, fmt="(T9,201F10.6)") (q_eigval(i), i=1, norb), &
829 sum(q_eigval(1:norb))
830 DO i = 1, norb
831 WRITE (unit=output_unit, fmt="(T3,I6,200F10.6)") i, q_eigvec(i, :)
832 END DO
833 DEALLOCATE (q_eigval)
834 DEALLOCATE (q_eigvec)
835 END IF
836 DEALLOCATE (q_work)
837 END IF ! debug
838 END IF ! should output
839
840 END DO ! next atom "iatom" of atomic kind "ikind"
841
842 IF (dft_control%mtlr_dft_with_perturbation .AND. &
843 dft_control%mtlr_ikind == ikind) THEN
844 dft_control%vhxc(ispin) = dft_control%vhxc(ispin)/real(n, dp)
845 CALL para_env%sum(dft_control%vhxc(ispin))
846 CALL para_env%sum(dft_control%trq(ispin))
847 END IF
848
849 DEALLOCATE (shell_sgf, q_matrix)
850
851 IF (ALLOCATED(q_opp_matrix)) THEN
852 DEALLOCATE (q_opp_matrix)
853 END IF
854
855 END DO ! next atomic kind "ikind"
856
857 ! Add V(i,j)[U] to V(i,j)[DFT]
858 IF (PRESENT(matrix_h) .OR. (PRESENT(matrix_w) .AND. do_kpoints)) THEN
859 IF (do_kpoints) THEN
860 block
861 REAL(dp) :: diagonal(nsgf)
862 CALL dbcsr_get_diag(sm_v, diagonal)
863 CALL para_env%sum(diagonal)
864 kpoints%lowdin_v(:, ispin) = diagonal(kpoints%lowdin_orbitals)
865 END block
866 ELSE
867 sm_h => matrix_h(ispin, 1)%matrix
868 CALL cp_dbcsr_sm_fm_multiply(sm_v, fm_s_half, fm_work1, nsgf)
869 CALL cp_fm_transpose(fm_work1, fm_work2)
870 CALL cp_dbcsr_plus_fm_fm_t(sm_h, fm_s_half, fm_work2, nsgf)
871 END IF
872 END IF ! An update of the Hamiltonian matrix is requested
873
874 ! Calculate the contribution (non-Pulay part) to the derivatives
875 ! w.r.t. the nuclear positions
876 IF (PRESENT(matrix_w) .AND. .NOT. do_kpoints) THEN
877
878 sm_p => matrix_p(ispin, 1)%matrix
879 sm_w => matrix_w(ispin, 1)%matrix
880
881 CALL cp_dbcsr_sm_fm_multiply(sm_v, fm_s_half, fm_work1, nsgf)
882 CALL cp_fm_transpose(fm_work1, fm_work2)
883 CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_work2, fm_work1, nsgf)
884 CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, 1.0_dp, fm_work1, fm_sev, 0.0_dp, fm_work2)
885 CALL parallel_gemm('T', 'N', nsgf, nsgf, nsgf, 1.0_dp, fm_sev, fm_work2, 0.0_dp, fm_work1)
886 CALL cp_fm_schur_product(fm_work1, slambda, fm_work2)
887 CALL cp_fm_transpose(fm_work2, fm_work1)
888 ! Symmetrize before projecting onto the stored triangular sparse blocks.
889 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_work2, matrix_b=fm_work1)
890 CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, -1.0_dp, fm_sev, fm_work2, 0.0_dp, fm_work1)
891 CALL cp_dbcsr_plus_fm_fm_t(sm_w, fm_work1, fm_sev, nsgf)
892
893 END IF ! W matrix update requested
894
895 END DO ! next spin "ispin"
896
897 IF (PRESENT(matrix_w)) THEN
898 IF (do_kpoints) CALL lowdin_kp_forces(kpoints)
899 CALL cp_fm_release(matrix=fm_sev)
900 CALL cp_fm_release(matrix=slambda)
901 END IF
902
903 ! Collect the energy contributions from all processes
904
905 CALL para_env%sum(energy%dft_plus_u)
906
907 IF (energy%dft_plus_u < 0.0_dp) THEN
908 CALL cp_warn(__location__, &
909 "DFT+U energy contribution is negative possibly due "// &
910 "to unphysical Lowdin charges!")
911 END IF
912
913 ! Release (local) full matrices
914 NULLIFY (fm_s_half)
915 CALL cp_fm_release(matrix=fm_work1)
916 CALL cp_fm_release(matrix=fm_work2)
917
918 ! Release (local) sparse matrices
919 CALL dbcsr_release(sm_q)
920 CALL dbcsr_release(sm_v)
921
922 IF (if_plus_j .AND. nspin == 2) THEN
923 CALL dbcsr_release(sm_q_opp)
924 END IF
925
926 CALL timestop(handle)
927
928 END SUBROUTINE lowdin
929
930! **************************************************************************************************
931!> \brief Select the AO columns used by the existing orbital-diagonal Lowdin functional.
932!> \param qs_env Quickstep environment, after S(k)**(1/2) has been prepared.
933! **************************************************************************************************
934 SUBROUTINE lowdin_kp_init(qs_env)
935 TYPE(qs_environment_type), POINTER :: qs_env
936
937 INTEGER :: i, ia, ik, ikind, iset, ish, lu, m, nao, &
938 nset, nspin
939 INTEGER, ALLOCATABLE :: first_atom(:)
940 INTEGER, POINTER :: atoms(:), first(:, :), l(:, :), &
941 last(:, :), nshell(:), orbitals(:)
942 LOGICAL :: active
943 LOGICAL, ALLOCATABLE :: selected(:)
944 REAL(kind=dp) :: j, u
945 TYPE(atomic_kind_type), POINTER :: kinds_set(:)
946 TYPE(dft_control_type), POINTER :: dft_control
947 TYPE(gto_basis_set_type), POINTER :: basis
948 TYPE(kpoint_env_type), POINTER :: kp
949 TYPE(kpoint_type), POINTER :: kpoints
950 TYPE(particle_type), POINTER :: particles(:)
951 TYPE(qs_kind_type), POINTER :: qs_kinds(:)
952
953 CALL get_qs_env(qs_env, kpoints=kpoints, atomic_kind_set=kinds_set, qs_kind_set=qs_kinds, &
954 particle_set=particles, dft_control=dft_control)
955 IF (dft_control%mtlr_dft_with_perturbation) THEN
956 cpabort("MTLR with Lowdin DFT+U and k-points is not implemented")
957 END IF
958 DO ik = 1, SIZE(kpoints%kp_sym)
959 IF (kpoints%kp_sym(ik)%kpoint_sym%apply_symmetry) THEN
960 cpabort("Lowdin DFT+U requires a full k-point grid or time-reversal-only reduction")
961 END IF
962 END DO
963 CALL get_qs_kind_set(qs_kinds, nsgf=nao)
964 nspin = dft_control%nspins
965 ALLOCATE (first_atom(SIZE(particles)), selected(nao))
966 CALL get_particle_set(particles, qs_kinds, first_sgf=first_atom)
967 selected = .false.
968 DO ikind = 1, SIZE(qs_kinds)
969 NULLIFY (orbitals)
970 CALL get_qs_kind(qs_kinds(ikind), dft_plus_u_atom=active, l_of_dft_plus_u=lu, &
971 basis_set=basis, u_minus_j_target=u, hund_j=j, orbitals=orbitals)
972 IF (.NOT. active) cycle
973 IF (ASSOCIATED(orbitals)) THEN
974 cpabort("Lowdin DFT+U with k-points does not support ENFORCE_OCCUPATION")
975 END IF
976 IF (u == 0.0_dp .AND. j == 0.0_dp) cycle
977 CALL get_atomic_kind(kinds_set(ikind), atom_list=atoms)
978 CALL get_gto_basis_set(basis, first_sgf=first, last_sgf=last, l=l, nset=nset, nshell=nshell)
979 DO ia = 1, SIZE(atoms)
980 DO iset = 1, nset
981 DO ish = 1, nshell(iset)
982 IF (l(ish, iset) /= lu) cycle
983 selected(first_atom(atoms(ia)) + first(ish, iset) - 1: &
984 first_atom(atoms(ia)) + last(ish, iset) - 1) = .true.
985 END DO
986 END DO
987 END DO
988 END DO
989 IF (ALLOCATED(kpoints%lowdin_orbitals)) THEN
990 DEALLOCATE (kpoints%lowdin_orbitals, kpoints%lowdin_q, kpoints%lowdin_q_new, kpoints%lowdin_v)
991 END IF
992 m = count(selected)
993 ALLOCATE (kpoints%lowdin_orbitals(m))
994 kpoints%lowdin_orbitals(:) = pack([(i, i=1, nao)], selected)
995 ALLOCATE (kpoints%lowdin_q(m, nspin), kpoints%lowdin_q_new(m, nspin), kpoints%lowdin_v(m, nspin))
996 kpoints%lowdin_q(:, :) = 0.0_dp
997 kpoints%lowdin_q_new(:, :) = 0.0_dp
998 kpoints%lowdin_v(:, :) = 0.0_dp
999 kpoints%lowdin_density_ready = .false.
1000 kpoints%lowdin_population_ready = .false.
1001 DO ik = 1, SIZE(kpoints%kp_env)
1002 kp => kpoints%kp_env(ik)%kpoint_env
1003 CALL lowdin_kp_prepare(kp, kpoints%lowdin_orbitals, kpoints%use_real_wfn)
1004 END DO
1005 END SUBROUTINE lowdin_kp_init
1006
1007! **************************************************************************************************
1008!> \brief Prepare selected Lowdin projectors and their reusable product workspace.
1009!> \param kp One k point; SCF already provides its overlap square root.
1010!> \param orbitals Selected global AO indices.
1011!> \param use_real_wfn Whether to use the real overlap square root.
1012!> \param smat Optional complex overlap for post-SCF reconstruction; remains unchanged.
1013! **************************************************************************************************
1014 SUBROUTINE lowdin_kp_prepare(kp, orbitals, use_real_wfn, smat)
1015 TYPE(kpoint_env_type), POINTER :: kp
1016 INTEGER, INTENT(IN) :: orbitals(:)
1017 LOGICAL, INTENT(IN) :: use_real_wfn
1018 TYPE(cp_cfm_type), INTENT(IN), OPTIONAL :: smat
1019
1020 INTEGER :: i, m, n
1021 REAL(dp), ALLOCATABLE :: eigenvalues(:)
1022 TYPE(cp_cfm_type) :: evec, work
1023
1024 CALL cp_fm_release(kp%lowdin_b)
1025 CALL cp_fm_release(kp%lowdin_work)
1026 CALL cp_cfm_release(kp%lowdin_cb)
1027 CALL cp_cfm_release(kp%lowdin_cwork)
1028 CALL kp%lowdin_gemm%destroy()
1029 m = SIZE(orbitals)
1030 IF (m == 0) RETURN
1031 IF (PRESENT(smat)) THEN
1032 cpassert(.NOT. use_real_wfn)
1033 CALL cp_cfm_get_info(smat, nrow_global=n)
1034 ALLOCATE (eigenvalues(n))
1035 CALL cp_cfm_create(evec, smat%matrix_struct)
1036 CALL cp_cfm_create(work, smat%matrix_struct)
1037 CALL cp_cfm_release(kp%cshalf)
1038 CALL cp_cfm_create(kp%cshalf, smat%matrix_struct)
1039 CALL cp_cfm_to_cfm(smat, work)
1040 CALL cp_cfm_heevd(work, evec, eigenvalues)
1041 cpassert(all(eigenvalues >= 0.0_dp))
1042 CALL cp_cfm_to_cfm(evec, work)
1043 CALL cp_cfm_column_scale(work, cmplx(sqrt(eigenvalues), 0.0_dp, dp))
1044 CALL parallel_gemm('N', 'C', n, n, n, z_one, work, evec, z_zero, kp%cshalf)
1045 CALL cp_cfm_release(evec)
1046 CALL cp_cfm_release(work)
1047 END IF
1048 CALL kp%lowdin_gemm%create(local_gemm_pu_gpu, timing=.false.)
1049 IF (use_real_wfn) THEN
1050 CALL cp_fm_create(kp%lowdin_b, kp%shalf%matrix_struct, ncol=m)
1051 CALL cp_fm_create(kp%lowdin_work, kp%lowdin_b%matrix_struct)
1052 DO i = 1, m
1053 CALL cp_fm_to_fm(kp%shalf, kp%lowdin_b, 1, orbitals(i), i)
1054 END DO
1055 ! The real driver also uses complex storage during its ADIIS/CDIIS handover.
1056 CALL cp_cfm_create(kp%lowdin_cb, kp%lowdin_b%matrix_struct)
1057 CALL cp_fm_to_cfm(kp%lowdin_b, mtarget=kp%lowdin_cb)
1058 ELSE
1059 CALL cp_cfm_create(kp%lowdin_cb, kp%cshalf%matrix_struct, ncol=m)
1060 DO i = 1, m
1061 CALL cp_cfm_to_cfm(kp%cshalf, kp%lowdin_cb, 1, orbitals(i), i)
1062 END DO
1063 END IF
1064 CALL cp_cfm_create(kp%lowdin_cwork, kp%lowdin_cb%matrix_struct)
1065 IF (PRESENT(smat)) CALL cp_cfm_release(kp%cshalf)
1066 END SUBROUTINE lowdin_kp_prepare
1067
1068! **************************************************************************************************
1069!> \brief Complete the initial density's Lowdin coordinates, including atomic and density restarts.
1070!> \param qs_env Quickstep environment with an initialized real-space density.
1071!> \param fmwork Global AO matrix workspace from the SCF environment being initialized.
1072! **************************************************************************************************
1073 SUBROUTINE lowdin_kp_initial_density(qs_env, fmwork)
1074 TYPE(qs_environment_type), POINTER :: qs_env
1075 TYPE(cp_fm_type), INTENT(IN) :: fmwork
1076
1077 INTEGER :: ik, ispin, kplocal, kpmax
1078 TYPE(cp_cfm_type) :: cpmat
1079 TYPE(cp_fm_type) :: pmat
1080 TYPE(dbcsr_p_type), POINTER :: density(:, :)
1081 TYPE(kpoint_env_type), POINTER :: kp
1082 TYPE(kpoint_operator_context_type) :: ctx
1083 TYPE(kpoint_type), POINTER :: kpoints
1084 TYPE(qs_rho_type), POINTER :: rho
1085
1086 CALL get_qs_env(qs_env, kpoints=kpoints, rho=rho)
1087 IF (.NOT. ALLOCATED(kpoints%lowdin_orbitals)) RETURN
1088 IF (.NOT. kpoints%lowdin_density_ready .AND. SIZE(kpoints%lowdin_orbitals) > 0) THEN
1089 ! The initial real-space guess defines a finite Fourier series. Subsequent
1090 ! accepted densities use their full Bloch matrices, never this sparse reconstruction.
1091 CALL qs_rho_get(rho, rho_ao_kp=density)
1092 CALL kpoint_operator_context_create(ctx, kpoints, fmwork, SIZE(density, 1))
1093 kp => kpoints%kp_env(1)%kpoint_env
1094 IF (kpoints%use_real_wfn) THEN
1095 CALL cp_fm_create(pmat, kp%shalf%matrix_struct)
1096 ELSE
1097 CALL cp_cfm_create(cpmat, kp%cshalf%matrix_struct)
1098 END IF
1099 kplocal = SIZE(kpoints%kp_env)
1100 kpmax = maxval(kpoints%kp_dist(2, :) - kpoints%kp_dist(1, :) + 1)
1101 kpoints%lowdin_q_new(:, :) = 0.0_dp
1102 DO ik = 1, kpmax
1103 DO ispin = 1, SIZE(density, 1)
1104 IF (kpoints%use_real_wfn) THEN
1105 CALL kpoint_operator_get(ctx, ik, ispin, density, fm_re=pmat, matrix_row=ispin)
1106 ELSE
1107 CALL kpoint_operator_get(ctx, ik, ispin, density, cmat=cpmat, matrix_row=ispin)
1108 END IF
1109 IF (ik > kplocal) cycle
1110 kp => kpoints%kp_env(ik)%kpoint_env
1111 IF (kpoints%use_real_wfn) THEN
1112 CALL lowdin_kp_accumulate(kp, kpoints%lowdin_q_new(:, ispin), pmat=pmat)
1113 ELSE
1114 CALL lowdin_kp_accumulate(kp, kpoints%lowdin_q_new(:, ispin), cpmat=cpmat)
1115 END IF
1116 END DO
1117 END DO
1118 CALL kpoints%para_env%sum(kpoints%lowdin_q_new)
1119 CALL kpoint_operator_context_release(ctx)
1120 CALL cp_fm_release(pmat)
1121 CALL cp_cfm_release(cpmat)
1122 kpoints%lowdin_population_ready = .true.
1123 END IF
1124 CALL lowdin_kp_population(kpoints, alpha=1.0_dp)
1125 END SUBROUTINE lowdin_kp_initial_density
1126
1127! **************************************************************************************************
1128!> \brief Cache the current occupations, optionally measure their residual and accept or mix them.
1129!> \param kpoints Optional k-point state; absent or inactive Lowdin state is a no-op.
1130!> \param delta Include the unmixed occupation residual in this convergence measure.
1131!> \param alpha Weight of the current occupations in the accepted state; absent leaves it unchanged.
1132! **************************************************************************************************
1133 SUBROUTINE lowdin_kp_population(kpoints, delta, alpha)
1134 TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
1135 REAL(kind=dp), INTENT(INOUT), OPTIONAL :: delta
1136 REAL(kind=dp), INTENT(IN), OPTIONAL :: alpha
1137
1138 IF (.NOT. PRESENT(kpoints)) RETURN
1139 IF (.NOT. ASSOCIATED(kpoints)) RETURN
1140 IF (.NOT. ALLOCATED(kpoints%lowdin_q)) RETURN
1141 IF (.NOT. kpoints%lowdin_population_ready) THEN
1142 IF (SIZE(kpoints%lowdin_q) > 0) THEN
1143 cpassert(kpoints%lowdin_density_ready)
1144 CALL lowdin_kp_project(kpoints, kpoints%lowdin_q_new)
1145 END IF
1146 kpoints%lowdin_population_ready = .true.
1147 END IF
1148 IF (PRESENT(delta)) delta = max(delta, maxval(abs(kpoints%lowdin_q_new - kpoints%lowdin_q)))
1149 IF (PRESENT(alpha)) THEN
1150 IF (alpha == 1.0_dp) THEN
1151 kpoints%lowdin_q(:, :) = kpoints%lowdin_q_new
1152 ELSE IF (alpha /= 0.0_dp) THEN
1153 kpoints%lowdin_q(:, :) = (1.0_dp - alpha)*kpoints%lowdin_q + alpha*kpoints%lowdin_q_new
1154 END IF
1155 END IF
1156 END SUBROUTINE lowdin_kp_population
1157
1158! **************************************************************************************************
1159!> \brief Project and reduce a Bloch density with one accumulator and matrix workspace per worker.
1160!> \param kpoints Selected overlap columns and parallel distribution.
1161!> \param q Weighted occupations; does not change the ground-state occupation caches.
1162!> \param density Optional split-complex response; absent uses the current Bloch density.
1163! **************************************************************************************************
1164 SUBROUTINE lowdin_kp_project(kpoints, q, density)
1165 TYPE(kpoint_type), POINTER :: kpoints
1166 REAL(kind=dp), CONTIGUOUS, INTENT(OUT) :: q(:, :)
1167 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: density(:, :, :)
1168
1169 INTEGER :: handle, ik, ispin, kplocal, m, nspin, &
1170 nworkers, thread
1171 TYPE(cp_cfm_type), ALLOCATABLE :: work(:)
1172 TYPE(kpoint_env_type), POINTER :: kp
1173
1174 m = SIZE(q, 1)
1175 IF (m == 0) RETURN
1176 CALL timeset("lowdin_kp_population", handle)
1177 nspin = SIZE(q, 2)
1178 kplocal = SIZE(kpoints%kp_env)
1179 kp => kpoints%kp_env(1)%kpoint_env
1180 nworkers = 1
1181!$ IF (kpoints%para_env_kp%num_pe == 1) nworkers = MIN(omp_get_max_threads(), kplocal)
1182 ALLOCATE (work(nworkers))
1183 q(:, :) = 0.0_dp
1184 IF (PRESENT(density) .OR. .NOT. kpoints%use_real_wfn) THEN
1185 DO thread = 1, nworkers
1186 CALL cp_cfm_create(work(thread), kp%pmat(1, 1)%matrix_struct)
1187 END DO
1188 END IF
1189!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) NUM_THREADS(nworkers) IF(nworkers > 1) &
1190!$OMP SHARED(kpoints,density,work,nspin,kplocal,nworkers) PRIVATE(ik,ispin,kp,thread) REDUCTION(+:q)
1191 DO ik = 1, kplocal
1192 thread = 1
1193!$ thread = omp_get_thread_num() + 1
1194 kp => kpoints%kp_env(ik)%kpoint_env
1195 DO ispin = 1, nspin
1196 IF (PRESENT(density)) THEN
1197 CALL cp_fm_to_cfm(density(ik, 1, ispin), density(ik, 2, ispin), work(thread))
1198 CALL lowdin_kp_accumulate(kp, q(:, ispin), cpmat=work(thread))
1199 ELSE IF (kpoints%use_real_wfn) THEN
1200 CALL lowdin_kp_accumulate(kp, q(:, ispin), pmat=kp%pmat(1, ispin))
1201 ELSE
1202 CALL cp_fm_to_cfm(kp%pmat(1, ispin), kp%pmat(2, ispin), work(thread))
1203 CALL lowdin_kp_accumulate(kp, q(:, ispin), cpmat=work(thread))
1204 END IF
1205 END DO
1206 END DO
1207!$OMP END PARALLEL DO
1208 ! Each entry is a local tile contraction: reduce once across groups AND group ranks.
1209 CALL kpoints%para_env%sum(q)
1210 DO thread = 1, nworkers
1211 CALL cp_cfm_release(work(thread))
1212 END DO
1213 CALL timestop(handle)
1214 END SUBROUTINE lowdin_kp_project
1215
1216! **************************************************************************************************
1217!> \brief Accumulate weighted local entries of diag(B**H P B), without a full transformed density.
1218!> \param kp One k point, with a private reusable product workspace.
1219!> \param q Accumulated local contributions; the caller zeroes and reduces the worker buffers.
1220!> \param pmat Real Bloch density.
1221!> \param cpmat Complex Bloch density (exclusive with pmat).
1222! **************************************************************************************************
1223 SUBROUTINE lowdin_kp_accumulate(kp, q, pmat, cpmat)
1224 TYPE(kpoint_env_type), POINTER :: kp
1225 REAL(kind=dp), INTENT(INOUT) :: q(:)
1226 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: pmat
1227 TYPE(cp_cfm_type), INTENT(IN), OPTIONAL :: cpmat
1228
1229 INTEGER :: j, m, n, nc, nr
1230 INTEGER, POINTER :: columns(:)
1231 LOGICAL :: local
1232
1233 m = SIZE(q)
1234 cpassert(PRESENT(pmat) .NEQV. PRESENT(cpmat))
1235 IF (PRESENT(pmat)) THEN
1236 CALL cp_fm_get_info(kp%lowdin_b, nrow_global=n, nrow_local=nr, ncol_local=nc, col_indices=columns)
1237 local = product(kp%lowdin_b%matrix_struct%context%num_pe) == 1
1238 IF (local) THEN
1239 CALL kp%lowdin_gemm%gemm('N', 'N', n, m, n, 1.0_dp, pmat%local_data, SIZE(pmat%local_data, 1), &
1240 kp%lowdin_b%local_data, SIZE(kp%lowdin_b%local_data, 1), 0.0_dp, &
1241 kp%lowdin_work%local_data, SIZE(kp%lowdin_work%local_data, 1))
1242 ELSE
1243 CALL parallel_gemm('N', 'N', n, m, n, 1.0_dp, pmat, kp%lowdin_b, 0.0_dp, kp%lowdin_work)
1244 END IF
1245 DO j = 1, nc
1246 q(columns(j)) = q(columns(j)) + kp%wkp* &
1247 sum(kp%lowdin_b%local_data(:nr, j)*kp%lowdin_work%local_data(:nr, j))
1248 END DO
1249 ELSE
1250 CALL cp_cfm_get_info(kp%lowdin_cb, nrow_global=n, nrow_local=nr, ncol_local=nc, col_indices=columns)
1251 CALL lowdin_kp_gemm('N', 'N', n, m, n, cpmat, kp%lowdin_cb, z_zero, kp%lowdin_cwork, kp%lowdin_gemm)
1252 DO j = 1, nc
1253 q(columns(j)) = q(columns(j)) + kp%wkp*real(sum(conjg(kp%lowdin_cb%local_data(:nr, j))* &
1254 kp%lowdin_cwork%local_data(:nr, j)), dp)
1255 END DO
1256 END IF
1257 END SUBROUTINE lowdin_kp_accumulate
1258
1259! **************************************************************************************************
1260!> \brief Add B diag(v) B**H before any solver, DIIS residual, or history storage uses H(k).
1261!> \param kp One k point, with its private workspace.
1262!> \param v Selected orbital potentials for the same raw/mixed state as the supplied Hamiltonian.
1263!> \param rmat Real Hamiltonian.
1264!> \param cmat Complex Hamiltonian (exclusive with rmat).
1265! **************************************************************************************************
1266 SUBROUTINE lowdin_kp_apply(kp, v, rmat, cmat)
1267 TYPE(kpoint_env_type), POINTER :: kp
1268 REAL(kind=dp), INTENT(IN) :: v(:)
1269 TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: rmat
1270 TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: cmat
1271
1272 INTEGER :: j, m, n, nc, nr
1273 INTEGER, POINTER :: columns(:)
1274 LOGICAL :: local
1275
1276 m = SIZE(v)
1277 IF (m == 0) RETURN
1278 IF (all(v == 0.0_dp)) RETURN
1279 cpassert(PRESENT(rmat) .NEQV. PRESENT(cmat))
1280 IF (PRESENT(rmat)) THEN
1281 CALL cp_fm_get_info(kp%lowdin_b, nrow_global=n, nrow_local=nr, ncol_local=nc, col_indices=columns)
1282 DO j = 1, nc
1283 kp%lowdin_work%local_data(:nr, j) = v(columns(j))*kp%lowdin_b%local_data(:nr, j)
1284 END DO
1285 local = product(kp%lowdin_b%matrix_struct%context%num_pe) == 1
1286 IF (local) THEN
1287 CALL kp%lowdin_gemm%gemm('N', 'T', n, n, m, 1.0_dp, &
1288 kp%lowdin_work%local_data, SIZE(kp%lowdin_work%local_data, 1), &
1289 kp%lowdin_b%local_data, SIZE(kp%lowdin_b%local_data, 1), 1.0_dp, &
1290 rmat%local_data, SIZE(rmat%local_data, 1))
1291 ELSE
1292 CALL parallel_gemm('N', 'T', n, n, m, 1.0_dp, kp%lowdin_work, kp%lowdin_b, 1.0_dp, rmat)
1293 END IF
1294 ELSE
1295 CALL cp_cfm_get_info(kp%lowdin_cb, nrow_global=n, nrow_local=nr, ncol_local=nc, col_indices=columns)
1296 DO j = 1, nc
1297 kp%lowdin_cwork%local_data(:nr, j) = v(columns(j))*kp%lowdin_cb%local_data(:nr, j)
1298 END DO
1299 CALL lowdin_kp_gemm('N', 'C', n, n, m, kp%lowdin_cwork, kp%lowdin_cb, z_one, cmat, kp%lowdin_gemm)
1300 END IF
1301 END SUBROUTINE lowdin_kp_apply
1302
1303! **************************************************************************************************
1304!> \brief Add B diag(v) B**H C to the complex orbital gradient, without forming an AO matrix.
1305!> \param kp K point owning the selected overlap columns and a private GEMM context.
1306!> \param v Orbital potential, without k weights or occupations (applied by the caller).
1307!> \param coeff Active complex orbitals.
1308!> \param hc H*C, updated in place.
1309! **************************************************************************************************
1310 SUBROUTINE lowdin_kp_apply_mo(kp, v, coeff, hc)
1311 TYPE(kpoint_env_type), POINTER :: kp
1312 REAL(kind=dp), INTENT(IN) :: v(:)
1313 TYPE(cp_cfm_type), INTENT(IN) :: coeff
1314 TYPE(cp_cfm_type), INTENT(INOUT) :: hc
1315
1316 INTEGER :: i, m, n, nc, nmo, nr
1317 INTEGER, POINTER :: rows(:)
1318 TYPE(cp_cfm_type) :: projected
1319
1320 m = SIZE(v)
1321 IF (m == 0 .OR. all(v == 0.0_dp)) RETURN
1322 CALL cp_cfm_get_info(coeff, nrow_global=n, ncol_global=nmo)
1323 CALL cp_cfm_create(projected, coeff%matrix_struct, nrow=m)
1324 CALL lowdin_kp_gemm('C', 'N', m, nmo, n, kp%lowdin_cb, coeff, z_zero, projected, kp%lowdin_gemm)
1325 CALL cp_cfm_get_info(projected, nrow_local=nr, ncol_local=nc, row_indices=rows)
1326 DO i = 1, nr
1327 projected%local_data(i, :nc) = v(rows(i))*projected%local_data(i, :nc)
1328 END DO
1329 CALL lowdin_kp_gemm('N', 'N', n, nmo, m, kp%lowdin_cb, projected, z_one, hc, kp%lowdin_gemm)
1330 CALL cp_cfm_release(projected)
1331 END SUBROUTINE lowdin_kp_apply_mo
1332
1333! **************************************************************************************************
1334!> \brief Project a Bloch density response and return its orbital potential response, including +J.
1335!> \param qs_env Quickstep environment with selected Lowdin projectors and current U/J.
1336!> \param density Unweighted split-complex Bloch density response (local k point, component, spin).
1337!> \param response Globally summed orbital potential response; ground-state caches are preserved.
1338! **************************************************************************************************
1339 SUBROUTINE lowdin_kp_kernel(qs_env, density, response)
1340 TYPE(qs_environment_type), POINTER :: qs_env
1341 TYPE(cp_fm_type), INTENT(IN) :: density(:, :, :)
1342 REAL(kind=dp), ALLOCATABLE, INTENT(OUT) :: response(:, :)
1343
1344 INTEGER :: i, iatom, ikind, m, nspin
1345 INTEGER, ALLOCATABLE :: last(:)
1346 REAL(kind=dp) :: j, u
1347 TYPE(kpoint_type), POINTER :: kpoints
1348 TYPE(particle_type), POINTER :: particles(:)
1349 TYPE(qs_kind_type), POINTER :: qs_kinds(:)
1350
1351 CALL get_qs_env(qs_env, kpoints=kpoints, particle_set=particles, qs_kind_set=qs_kinds)
1352 m = SIZE(kpoints%lowdin_orbitals)
1353 nspin = SIZE(density, 3)
1354 ALLOCATE (response(m, nspin))
1355 response(:, :) = 0.0_dp
1356 IF (m == 0) RETURN
1357 CALL lowdin_kp_project(kpoints, response, density)
1358 ALLOCATE (last(SIZE(particles)))
1359 CALL get_particle_set(particles, qs_kinds, last_sgf=last)
1360 iatom = 1
1361 DO i = 1, m
1362 DO WHILE (kpoints%lowdin_orbitals(i) > last(iatom))
1363 iatom = iatom + 1
1364 END DO
1365 CALL get_atomic_kind(particles(iatom)%atomic_kind, kind_number=ikind)
1366 CALL get_qs_kind(qs_kinds(ikind), u_minus_j=u, hund_j=j)
1367 response(i, :) = 0.5_dp*real(nspin, dp)*(j*response(i, nspin:1:-1) - u*response(i, :))
1368 END DO
1369 END SUBROUTINE lowdin_kp_kernel
1370
1371! **************************************************************************************************
1372!> \brief Add -dE_U/dS(k) to W(k), before its ordinary inverse Fourier transform.
1373!> If A=S**(1/2), A*dA+dA*A=dS. In the overlap eigenbasis its adjoint is
1374!> G_ij=(U**H*(P*A*V+V*A*P)*U)_ij/(sqrt(s_i)+sqrt(s_j)).
1375!> W <- W-U*G*U**H; the existing overlap derivatives give forces and virial.
1376!> k-point weights belong to the Fourier transform and are not applied here.
1377!> \param kpoints Final Bloch densities, orbital potentials and cached overlap spectra.
1378! **************************************************************************************************
1379 SUBROUTINE lowdin_kp_forces(kpoints)
1380 TYPE(kpoint_type), POINTER :: kpoints
1381
1382 INTEGER :: handle, ik, ispin, iw, nworkers, thread
1383 TYPE(cp_cfm_type), ALLOCATABLE :: work(:, :)
1384 TYPE(kpoint_env_type), POINTER :: kp
1385 TYPE(local_gemm_ctxt_type), ALLOCATABLE :: gemm(:)
1386
1387 IF (SIZE(kpoints%lowdin_orbitals) == 0) RETURN
1388 IF (all(kpoints%lowdin_v == 0.0_dp)) RETURN
1389 CALL timeset("lowdin_kp_forces", handle)
1390 nworkers = 1
1391!$ IF (kpoints%para_env_kp%num_pe == 1) &
1392!$ nworkers = MIN(omp_get_max_threads(), SIZE(kpoints%kp_env)*SIZE(kpoints%lowdin_v, 2))
1393 ALLOCATE (work(3, nworkers), gemm(nworkers))
1394 kp => kpoints%kp_env(1)%kpoint_env
1395 DO thread = 1, nworkers
1396 DO iw = 1, 2
1397 CALL cp_cfm_create(work(iw, thread), kp%pmat(1, 1)%matrix_struct)
1398 END DO
1399 CALL cp_cfm_create(work(3, thread), kp%lowdin_cb%matrix_struct)
1400 CALL gemm(thread)%create(local_gemm_pu_gpu, timing=.false.)
1401 END DO
1402!$OMP PARALLEL DO COLLAPSE(2) DEFAULT(NONE) SCHEDULE(DYNAMIC, 1) NUM_THREADS(nworkers) IF(nworkers > 1) &
1403!$OMP SHARED(kpoints,work,gemm,nworkers) PRIVATE(ik,ispin,thread)
1404 DO ispin = 1, SIZE(kpoints%lowdin_v, 2)
1405 DO ik = 1, SIZE(kpoints%kp_env)
1406 thread = 1
1407!$ thread = omp_get_thread_num() + 1
1408 CALL lowdin_kp_force_matrix(kpoints%kp_env(ik)%kpoint_env, kpoints%lowdin_orbitals, ispin, &
1409 kpoints%lowdin_v(:, ispin), work(:, thread), gemm(thread))
1410 END DO
1411 END DO
1412!$OMP END PARALLEL DO
1413 DO thread = 1, nworkers
1414 CALL gemm(thread)%destroy()
1415 DO iw = 1, 3
1416 CALL cp_cfm_release(work(iw, thread))
1417 END DO
1418 END DO
1419 CALL timestop(handle)
1420 END SUBROUTINE lowdin_kp_forces
1421
1422! **************************************************************************************************
1423!> \brief One spin/k-point's Hermitian overlap derivative, factoring through the selected AO columns.
1424!> \param kp K-point state; W is updated, P and the overlap spectrum are preserved.
1425!> \param orbitals Selected AO columns.
1426!> \param ispin Spin channel to update.
1427!> \param v Derivatives of the common U/J functional with respect to this spin's q.
1428!> \param work Two square matrices and selected-column scratch, private to this worker.
1429!> \param gemm Private local GEMM context; concurrent spins must not share the k-point context.
1430! **************************************************************************************************
1431 SUBROUTINE lowdin_kp_force_matrix(kp, orbitals, ispin, v, work, gemm)
1432 TYPE(kpoint_env_type), POINTER :: kp
1433 INTEGER, INTENT(IN) :: orbitals(:), ispin
1434 REAL(kind=dp), INTENT(IN) :: v(:)
1435 TYPE(cp_cfm_type), INTENT(INOUT) :: work(3)
1436 TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: gemm
1437
1438 COMPLEX(KIND=dp) :: value
1439 INTEGER :: i, j, m, n, nc, nr
1440 INTEGER, POINTER :: columns(:), rows(:)
1441 LOGICAL :: local, real_only
1442
1443 IF (all(v == 0.0_dp)) RETURN
1444 cpassert(ALLOCATED(kp%lowdin_sqrt))
1445 cpassert(ASSOCIATED(kp%lowdin_evec%local_data))
1446 cpassert(ASSOCIATED(kp%wmat))
1447 CALL cp_cfm_get_info(work(1), nrow_global=n, nrow_local=nr, ncol_local=nc, &
1448 row_indices=rows, col_indices=columns)
1449 m = SIZE(orbitals)
1450 local = product(work(1)%matrix_struct%context%num_pe) == 1
1451 real_only = SIZE(kp%pmat, 1) == 1
1452 IF (real_only) THEN
1453 CALL cp_fm_to_cfm(kp%pmat(1, ispin), mtarget=work(1))
1454 ELSE
1455 CALL cp_fm_to_cfm(kp%pmat(1, ispin), kp%pmat(2, ispin), work(1))
1456 END IF
1457 ! U**H * P * A * V * U = (U**H * P * B) * diag(v) * U(orbitals,:).
1458 ! Keep these products thin: only the final backtransform needs square GEMMs.
1459 CALL lowdin_kp_gemm('N', 'N', n, m, n, work(1), kp%lowdin_cb, z_zero, work(2), gemm)
1460 CALL lowdin_kp_gemm('C', 'N', n, m, n, kp%lowdin_evec, work(2), z_zero, work(1), gemm)
1461 DO j = 1, nc
1462 IF (columns(j) <= m) work(1)%local_data(:nr, j) = v(columns(j))*work(1)%local_data(:nr, j)
1463 END DO
1464 IF (local) THEN
1465 DO j = 1, m
1466 work(3)%local_data(:n, j) = conjg(kp%lowdin_evec%local_data(orbitals(j), :n))
1467 END DO
1468 ELSE
1469 CALL cp_cfm_transpose(kp%lowdin_evec, 'C', work(2))
1470 DO j = 1, m
1471 CALL cp_cfm_to_cfm(work(2), work(3), 1, orbitals(j), j)
1472 END DO
1473 END IF
1474 CALL lowdin_kp_gemm('N', 'C', n, n, m, work(1), work(3), z_zero, work(2), gemm)
1475 ! Eigenvalue sums remain regular for equal or nearly equal overlap eigenvalues.
1476 DO j = 1, nc
1477 DO i = 1, nr
1478 work(2)%local_data(i, j) = work(2)%local_data(i, j)/ &
1479 (kp%lowdin_sqrt(rows(i)) + kp%lowdin_sqrt(columns(j)))
1480 END DO
1481 END DO
1482 CALL lowdin_kp_gemm('N', 'N', n, n, n, kp%lowdin_evec, work(2), z_zero, work(1), gemm)
1483 CALL lowdin_kp_gemm('N', 'C', n, n, n, work(1), kp%lowdin_evec, z_zero, work(2), gemm)
1484 IF (.NOT. local) CALL cp_cfm_transpose(work(2), 'C', work(1))
1485 DO j = 1, nc
1486 DO i = 1, nr
1487 IF (local) THEN
1488 value = work(2)%local_data(i, j) + conjg(work(2)%local_data(j, i))
1489 ELSE
1490 value = work(2)%local_data(i, j) + work(1)%local_data(i, j)
1491 END IF
1492 kp%wmat(1, ispin)%local_data(i, j) = kp%wmat(1, ispin)%local_data(i, j) - real(value, dp)
1493 IF (.NOT. real_only) THEN
1494 kp%wmat(2, ispin)%local_data(i, j) = kp%wmat(2, ispin)%local_data(i, j) - aimag(value)
1495 END IF
1496 END DO
1497 END DO
1498 END SUBROUTINE lowdin_kp_force_matrix
1499
1500! **************************************************************************************************
1501!> \brief C = op(A)*op(B) + beta*C, with a private local context or collective parallel GEMM.
1502!> \param ta Left transpose flag.
1503!> \param tb Right transpose flag.
1504!> \param nrow Output rows.
1505!> \param ncol Output columns.
1506!> \param inner Contraction dimension; may use only the selected columns of a square workspace.
1507!> \param a Left operand.
1508!> \param b Right operand.
1509!> \param beta Scale of the existing output.
1510!> \param c Output, distinct from both inputs.
1511!> \param gemm Context owned by the current k-point or force worker.
1512! **************************************************************************************************
1513 SUBROUTINE lowdin_kp_gemm(ta, tb, nrow, ncol, inner, a, b, beta, c, gemm)
1514 CHARACTER, INTENT(IN) :: ta, tb
1515 INTEGER, INTENT(IN) :: nrow, ncol, inner
1516 TYPE(cp_cfm_type), INTENT(IN) :: a, b
1517 COMPLEX(KIND=dp), INTENT(IN) :: beta
1518 TYPE(cp_cfm_type), INTENT(INOUT) :: c
1519 TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: gemm
1520
1521 IF (product(c%matrix_struct%context%num_pe) == 1) THEN
1522 CALL gemm%gemm(ta, tb, nrow, ncol, inner, z_one, &
1523 a%local_data, SIZE(a%local_data, 1), b%local_data, SIZE(b%local_data, 1), &
1524 beta, c%local_data, SIZE(c%local_data, 1))
1525 ELSE
1526 CALL parallel_gemm(ta, tb, nrow, ncol, inner, z_one, a, b, beta, c)
1527 END IF
1528 END SUBROUTINE lowdin_kp_gemm
1529
1530! **************************************************************************************************
1531!> \brief Add a DFT+U contribution to the Hamiltonian matrix\n
1532!> using a method based on the Mulliken population analysis
1533!> \f[q_{\mu\nu} = \frac{1}{2} (P_{\mu\nu} S_{\nu\mu} +
1534!> S_{\mu\nu} P_{\nu\mu})\f]
1535!> where \b P and \b S are the density and the
1536!> overlap matrix, respectively.
1537!> \param[in] qs_env Quickstep environment
1538!> \param orthonormal_basis ...
1539!> \param[in,out] matrix_h Hamiltonian matrices for each spin
1540!> \param should_output ...
1541!> \param output_unit ...
1542!> \param print_level ...
1543!> \date 03.07.2008
1544!> \par
1545!> \f{eqnarray*}{
1546!> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}\\\
1547!> & = & E^{\rm DFT} + \frac{1}{2}\sum_A(U_A - J_A)(Tr(q_A) - Tr(q^2_A))\\‍[1ex]
1548!> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
1549!> & = & \frac{\partial E^{\rm DFT}}
1550!> {\partial P_{\mu\nu}} +
1551!> \frac{\partial E^{\rm U}}
1552!> {\partial P_{\mu\nu}}\\\
1553!> & = & H_{\mu\nu} + \sum_A
1554!> \frac{\partial E^{\rm U}}{\partial q_A}
1555!> \frac{\partial q_A}{\partial P_{\mu\nu}}\\\
1556!> \f}
1557!> \author Matthias Krack (MK)
1558!> \version 1.0
1559! **************************************************************************************************
1560 SUBROUTINE mulliken(qs_env, orthonormal_basis, matrix_h, should_output, &
1561 output_unit, print_level)
1562
1563 TYPE(qs_environment_type), POINTER :: qs_env
1564 LOGICAL, INTENT(IN) :: orthonormal_basis
1565 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
1566 POINTER :: matrix_h
1567 LOGICAL, INTENT(IN) :: should_output
1568 INTEGER, INTENT(IN) :: output_unit, print_level
1569
1570 CHARACTER(LEN=*), PARAMETER :: routinen = 'mulliken'
1571
1572 CHARACTER(LEN=10) :: spin_info
1573 CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: symbol
1574 CHARACTER(LEN=default_string_length) :: atomic_kind_name
1575 INTEGER :: atom_a, handle, i, i0, iatom, ic, ikind, iorb, isb, iset, isgf, ishell, ispin, j, &
1576 jsb, jset, jsgf, jshell, lu, m, max_scf, n, natom, natom_of_kind, nimg, nkind, norb, nsb, &
1577 nsbsize, nset, nsgf_kind, nspin
1578 INTEGER, DIMENSION(1) :: iloc
1579 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals
1580 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
1581 LOGICAL :: debug, dft_plus_u_atom, found, &
1582 just_energy, occupation_enforced, smear
1583 LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_plus_u_kind, orb_occ
1584 REAL(kind=dp) :: eps_scf, eps_u_ramping, fspin, occ, trq, &
1585 trq2, u_minus_j, u_minus_j_target, &
1586 u_ramping
1587 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: q_eigval
1588 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: q_eigvec, q_matrix, q_work
1589 REAL(kind=dp), DIMENSION(:), POINTER :: nelec
1590 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, p_block, q_block, s_block, &
1591 v_block
1592 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1593 TYPE(atomic_kind_type), POINTER :: kind_a
1594 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
1595 TYPE(dbcsr_type), POINTER :: sm_h, sm_p, sm_q, sm_s, sm_v
1596 TYPE(dft_control_type), POINTER :: dft_control
1597 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1598 TYPE(mp_para_env_type), POINTER :: para_env
1599 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1600 TYPE(qs_energy_type), POINTER :: energy
1601 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1602 TYPE(qs_rho_type), POINTER :: rho
1603
1604 CALL timeset(routinen, handle)
1605
1606 debug = .false. ! Set to .TRUE. to print debug information
1607
1608 NULLIFY (atom_list)
1609 NULLIFY (atomic_kind_set)
1610 NULLIFY (qs_kind_set)
1611 NULLIFY (dft_control)
1612 NULLIFY (energy)
1613 NULLIFY (first_sgf)
1614 NULLIFY (h_block)
1615 NULLIFY (matrix_p)
1616 NULLIFY (matrix_s)
1617 NULLIFY (l)
1618 NULLIFY (last_sgf)
1619 NULLIFY (nelec)
1620 NULLIFY (nshell)
1621 NULLIFY (orb_basis_set)
1622 NULLIFY (p_block)
1623 NULLIFY (particle_set)
1624 NULLIFY (q_block)
1625 NULLIFY (rho)
1626 NULLIFY (s_block)
1627 NULLIFY (orbitals)
1628 NULLIFY (sm_h)
1629 NULLIFY (sm_p)
1630 NULLIFY (sm_q)
1631 NULLIFY (sm_s)
1632 NULLIFY (sm_v)
1633 NULLIFY (v_block)
1634 NULLIFY (para_env)
1635
1636 smear = .false.
1637 max_scf = -1
1638 eps_scf = 1.0e30_dp
1639 occupation_enforced = .false.
1640
1641 CALL get_qs_env(qs_env=qs_env, &
1642 atomic_kind_set=atomic_kind_set, &
1643 qs_kind_set=qs_kind_set, &
1644 dft_control=dft_control, &
1645 energy=energy, &
1646 particle_set=particle_set, &
1647 rho=rho, &
1648 para_env=para_env)
1649
1650 cpassert(ASSOCIATED(atomic_kind_set))
1651 cpassert(ASSOCIATED(dft_control))
1652 cpassert(ASSOCIATED(energy))
1653 cpassert(ASSOCIATED(particle_set))
1654 cpassert(ASSOCIATED(rho))
1655
1656 IF (orthonormal_basis) THEN
1657 NULLIFY (sm_s)
1658 ELSE
1659 ! Get overlap matrix in sparse format
1660 CALL get_qs_env(qs_env=qs_env, &
1661 matrix_s_kp=matrix_s)
1662 cpassert(ASSOCIATED(matrix_s))
1663 END IF
1664 nimg = dft_control%nimages
1665
1666 ! Get density matrices in sparse format
1667
1668 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1669
1670 energy%dft_plus_u = 0.0_dp
1671
1672 nspin = dft_control%nspins
1673
1674 IF (nspin == 2) THEN
1675 fspin = 1.0_dp
1676 ELSE
1677 fspin = 0.5_dp
1678 END IF
1679
1680 ! Get the total number of atoms, contracted spherical Gaussian basis
1681 ! functions, and atomic kinds
1682
1683 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1684 natom=natom)
1685
1686 nkind = SIZE(atomic_kind_set)
1687
1688 ALLOCATE (is_plus_u_kind(nkind))
1689 is_plus_u_kind(:) = .false.
1690
1691 IF (PRESENT(matrix_h)) THEN
1692 just_energy = .false.
1693 ELSE
1694 just_energy = .true.
1695 END IF
1696
1697 ! Loop over all spins
1698 DO ispin = 1, nspin
1699
1700 ! Loop over cell images
1701 DO ic = 1, nimg
1702 IF (.NOT. orthonormal_basis) THEN
1703 sm_s => matrix_s(1, ic)%matrix
1704 END IF
1705
1706 IF (PRESENT(matrix_h)) THEN
1707 ! Hamiltonian matrix for spin ispin in sparse format
1708 sm_h => matrix_h(ispin, ic)%matrix
1709 ELSE
1710 NULLIFY (sm_h)
1711 END IF
1712
1713 ! Get density matrix for spin ispin in sparse format
1714
1715 sm_p => matrix_p(ispin, ic)%matrix
1716
1717 IF (.NOT. ASSOCIATED(sm_q)) THEN
1718 ALLOCATE (sm_q)
1719 CALL dbcsr_get_block_diag(sm_p, sm_q)
1720 END IF
1721 CALL dbcsr_set(sm_q, 0.0_dp)
1722
1723 IF (.NOT. ASSOCIATED(sm_v)) THEN
1724 ALLOCATE (sm_v)
1725 CALL dbcsr_get_block_diag(sm_p, sm_v)
1726 END IF
1727 CALL dbcsr_set(sm_v, 0.0_dp)
1728
1729 DO iatom = 1, natom
1730
1731 CALL dbcsr_get_block_p(matrix=sm_p, &
1732 row=iatom, &
1733 col=iatom, &
1734 block=p_block, &
1735 found=found)
1736
1737 IF (.NOT. ASSOCIATED(p_block)) cycle
1738
1739 CALL dbcsr_get_block_p(matrix=sm_q, &
1740 row=iatom, &
1741 col=iatom, &
1742 block=q_block, &
1743 found=found)
1744 cpassert(ASSOCIATED(q_block))
1745
1746 IF (orthonormal_basis) THEN
1747 ! S is the unit matrix
1748 DO isgf = 1, SIZE(q_block, 1)
1749 q_block(isgf, isgf) = p_block(isgf, isgf)
1750 END DO
1751 ELSE
1752 CALL dbcsr_get_block_p(matrix=sm_s, &
1753 row=iatom, &
1754 col=iatom, &
1755 block=s_block, &
1756 found=found)
1757 cpassert(ASSOCIATED(s_block))
1758 ! Exploit that P and S are symmetric
1759 DO jsgf = 1, SIZE(p_block, 2)
1760 DO isgf = 1, SIZE(p_block, 1)
1761 q_block(isgf, jsgf) = p_block(isgf, jsgf)*s_block(isgf, jsgf)
1762 END DO
1763 END DO
1764 END IF ! orthonormal basis set
1765
1766 END DO ! next atom "iatom"
1767
1768 ! E[DFT+U] = E[DFT] + E[U]
1769 ! = E[DFT] + (U - J)*(Tr(q) - Tr(q*q))/2
1770
1771 ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U]
1772 ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j)
1773 ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j))
1774
1775 ! Loop over all atomic kinds
1776
1777 DO ikind = 1, nkind
1778
1779 ! Load the required atomic kind data
1780
1781 CALL get_atomic_kind(atomic_kind_set(ikind), &
1782 atom_list=atom_list, &
1783 name=atomic_kind_name, &
1784 natom=natom_of_kind)
1785
1786 CALL get_qs_kind(qs_kind_set(ikind), &
1787 dft_plus_u_atom=dft_plus_u_atom, &
1788 l_of_dft_plus_u=lu, &
1789 nsgf=nsgf_kind, &
1790 basis_set=orb_basis_set, &
1791 u_minus_j=u_minus_j, &
1792 u_minus_j_target=u_minus_j_target, &
1793 u_ramping=u_ramping, &
1794 eps_u_ramping=eps_u_ramping, &
1795 nelec=nelec, &
1796 orbitals=orbitals, &
1797 eps_scf=eps_scf, &
1798 max_scf=max_scf, &
1799 smear=smear)
1800
1801 ! Check, if the atoms of this atomic kind need a DFT+U correction
1802
1803 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
1804 IF (.NOT. dft_plus_u_atom) cycle
1805 IF (lu < 0) cycle
1806
1807 ! Apply U ramping if requested
1808
1809 IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
1810 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
1811 u_minus_j = min(u_minus_j + u_ramping, u_minus_j_target)
1812 CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
1813 END IF
1814 IF (should_output .AND. (output_unit > 0)) THEN
1815 WRITE (unit=output_unit, fmt="(T3,A,3X,A,F0.3,A)") &
1816 "Kind name: "//trim(adjustl(atomic_kind_name)), &
1817 "U(eff) = ", u_minus_j*evolt, " eV"
1818 END IF
1819 END IF
1820
1821 IF (u_minus_j == 0.0_dp) cycle
1822
1823 is_plus_u_kind(ikind) = .true.
1824
1825 ! Load the required Gaussian basis set data
1826
1827 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1828 first_sgf=first_sgf, &
1829 l=l, &
1830 last_sgf=last_sgf, &
1831 nset=nset, &
1832 nshell=nshell)
1833
1834 ! Count the relevant shell blocks of this atomic kind
1835
1836 nsb = 0
1837 DO iset = 1, nset
1838 DO ishell = 1, nshell(iset)
1839 IF (l(ishell, iset) == lu) nsb = nsb + 1
1840 END DO
1841 END DO
1842
1843 nsbsize = (2*lu + 1)
1844 n = nsb*nsbsize
1845
1846 ALLOCATE (q_matrix(n, n))
1847 q_matrix(:, :) = 0.0_dp
1848
1849 ! Print headline if requested
1850
1851 IF (should_output .AND. (print_level > low_print_level)) THEN
1852 IF (output_unit > 0) THEN
1853 ALLOCATE (symbol(nsbsize))
1854 DO m = -lu, lu
1855 symbol(lu + m + 1) = sgf_symbol(0, lu, m)
1856 END DO
1857 IF (nspin > 1) THEN
1858 WRITE (unit=spin_info, fmt="(A8,I2)") " of spin", ispin
1859 ELSE
1860 spin_info = ""
1861 END IF
1862 WRITE (unit=output_unit, fmt="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
1863 "DFT+U occupations"//trim(spin_info)//" for the atoms of atomic kind ", ikind, &
1864 ": "//trim(atomic_kind_name), &
1865 "Atom Shell ", (adjustr(symbol(i)), i=1, nsbsize), " Trace"
1866 DEALLOCATE (symbol)
1867 END IF
1868 END IF
1869
1870 ! Loop over all atoms of the current atomic kind
1871
1872 DO iatom = 1, natom_of_kind
1873
1874 atom_a = atom_list(iatom)
1875
1876 q_matrix(:, :) = 0.0_dp
1877
1878 ! Get diagonal block
1879
1880 CALL dbcsr_get_block_p(matrix=sm_q, &
1881 row=atom_a, &
1882 col=atom_a, &
1883 block=q_block, &
1884 found=found)
1885
1886 ! Calculate energy contribution to E(U)
1887
1888 IF (ASSOCIATED(q_block)) THEN
1889
1890 i = 0
1891 DO iset = 1, nset
1892 DO ishell = 1, nshell(iset)
1893 IF (l(ishell, iset) /= lu) cycle
1894 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1895 i = i + 1
1896 j = 0
1897 DO jset = 1, nset
1898 DO jshell = 1, nshell(jset)
1899 IF (l(jshell, jset) /= lu) cycle
1900 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
1901 j = j + 1
1902 q_matrix(i, j) = q_block(isgf, jsgf)
1903 END DO ! next contracted spherical Gaussian function "jsgf"
1904 END DO ! next shell "jshell"
1905 END DO ! next shell set "jset"
1906 END DO ! next contracted spherical Gaussian function "isgf"
1907 END DO ! next shell "ishell"
1908 END DO ! next shell set "iset"
1909
1910 ! Perform the requested manipulations of the (initial) orbital occupations
1911
1912 IF (ASSOCIATED(orbitals)) THEN
1913 IF ((qs_env%scf_env%iter_delta >= eps_scf) .OR. &
1914 ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
1915 (qs_env%scf_env%iter_count <= max_scf))) THEN
1916 ALLOCATE (orb_occ(nsbsize))
1917 ALLOCATE (q_eigval(n))
1918 q_eigval(:) = 0.0_dp
1919 ALLOCATE (q_eigvec(n, n))
1920 q_eigvec(:, :) = 0.0_dp
1921 norb = SIZE(orbitals)
1922 CALL jacobi(q_matrix, q_eigval, q_eigvec)
1923 q_matrix(:, :) = 0.0_dp
1924 IF (nelec(ispin) >= 0.5_dp) THEN
1925 trq = nelec(ispin)/sum(q_eigval(1:n))
1926 q_eigval(1:n) = trq*q_eigval(1:n)
1927 END IF
1928 DO isb = 1, nsb
1929 trq = 0.0_dp
1930 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
1931 trq = trq + q_eigval(i)
1932 END DO
1933 IF (smear) THEN
1934 occ = trq/real(norb, kind=dp)
1935 ELSE
1936 occ = 1.0_dp/fspin
1937 END IF
1938 orb_occ(:) = .false.
1939 iloc = maxloc(q_eigvec(:, isb*nsbsize))
1940 jsb = int((iloc(1) - 1)/nsbsize) + 1
1941 i = 0
1942 i0 = (jsb - 1)*nsbsize + 1
1943 iorb = -1000
1944 DO j = i0, jsb*nsbsize
1945 i = i + 1
1946 IF (i > norb) THEN
1947 DO m = -lu, lu
1948 IF (.NOT. orb_occ(lu + m + 1)) THEN
1949 iorb = i0 + lu + m
1950 orb_occ(lu + m + 1) = .true.
1951 END IF
1952 END DO
1953 ELSE
1954 iorb = i0 + lu + orbitals(i)
1955 orb_occ(lu + orbitals(i) + 1) = .true.
1956 END IF
1957 cpassert(iorb /= -1000)
1958 iloc = maxloc(q_eigvec(iorb, :))
1959 q_eigval(iloc(1)) = min(occ, trq)
1960 q_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1)) ! backtransform left
1961 trq = trq - q_eigval(iloc(1))
1962 END DO
1963 END DO
1964 q_matrix(:, :) = matmul(q_matrix, transpose(q_eigvec)) ! backtransform right
1965 DEALLOCATE (orb_occ)
1966 DEALLOCATE (q_eigval)
1967 DEALLOCATE (q_eigvec)
1968 occupation_enforced = .true.
1969 END IF
1970 END IF ! orbitals associated
1971
1972 trq = 0.0_dp
1973 trq2 = 0.0_dp
1974
1975 DO i = 1, n
1976 trq = trq + q_matrix(i, i)
1977 DO j = 1, n
1978 trq2 = trq2 + q_matrix(i, j)*q_matrix(j, i)
1979 END DO
1980 END DO
1981
1982 trq = fspin*trq
1983 trq2 = fspin*fspin*trq2
1984
1985 ! Calculate energy contribution to E(U)
1986
1987 energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
1988
1989 ! Calculate potential V(U) = dE(U)/dq
1990
1991 IF (.NOT. just_energy) THEN
1992
1993 CALL dbcsr_get_block_p(matrix=sm_v, &
1994 row=atom_a, &
1995 col=atom_a, &
1996 block=v_block, &
1997 found=found)
1998 cpassert(ASSOCIATED(v_block))
1999
2000 i = 0
2001 DO iset = 1, nset
2002 DO ishell = 1, nshell(iset)
2003 IF (l(ishell, iset) /= lu) cycle
2004 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
2005 i = i + 1
2006 j = 0
2007 DO jset = 1, nset
2008 DO jshell = 1, nshell(jset)
2009 IF (l(jshell, jset) /= lu) cycle
2010 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
2011 j = j + 1
2012 IF (isgf == jsgf) THEN
2013 v_block(isgf, isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i, i))
2014 ELSE
2015 v_block(isgf, jsgf) = -u_minus_j*fspin*q_matrix(j, i)
2016 END IF
2017 END DO ! next contracted spherical Gaussian function "jsgf"
2018 END DO ! next shell "jshell"
2019 END DO ! next shell set "jset"
2020 END DO ! next contracted spherical Gaussian function "isgf"
2021 END DO ! next shell "ishell"
2022 END DO ! next shell set "iset"
2023
2024 END IF ! not just energy
2025
2026 END IF ! q_block associated
2027
2028 ! Consider print requests
2029
2030 IF (should_output .AND. (print_level > low_print_level)) THEN
2031 CALL para_env%sum(q_matrix)
2032 IF (output_unit > 0) THEN
2033 ALLOCATE (q_work(nsb, nsbsize))
2034 q_work(:, :) = 0.0_dp
2035 DO isb = 1, nsb
2036 j = 0
2037 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
2038 j = j + 1
2039 q_work(isb, j) = q_matrix(i, i)
2040 END DO
2041 END DO
2042 DO isb = 1, nsb
2043 WRITE (unit=output_unit, fmt="(T3,I6,2X,I6,2X,10F8.3)") &
2044 atom_a, isb, q_work(isb, :), sum(q_work(isb, :))
2045 END DO
2046 WRITE (unit=output_unit, fmt="(T12,A,2X,10F8.3)") &
2047 "Total", (sum(q_work(:, i)), i=1, nsbsize), sum(q_work)
2048 WRITE (unit=output_unit, fmt="(A)") ""
2049 DEALLOCATE (q_work)
2050 IF (debug) THEN
2051 ! Print the DFT+U occupation matrix
2052 WRITE (unit=output_unit, fmt="(T9,70I10)") (i, i=1, n)
2053 DO i = 1, n
2054 WRITE (unit=output_unit, fmt="(T3,I6,70F10.6)") i, q_matrix(i, :)
2055 END DO
2056 ! Print the eigenvalues and eigenvectors of the occupation matrix
2057 ALLOCATE (q_eigval(n))
2058 q_eigval(:) = 0.0_dp
2059 ALLOCATE (q_eigvec(n, n))
2060 q_eigvec(:, :) = 0.0_dp
2061 CALL jacobi(q_matrix, q_eigval, q_eigvec)
2062 WRITE (unit=output_unit, fmt="(/,T9,70I10)") (i, i=1, n)
2063 WRITE (unit=output_unit, fmt="(T9,71F10.6)") (q_eigval(i), i=1, n), &
2064 sum(q_eigval(1:n))
2065 DO i = 1, n
2066 WRITE (unit=output_unit, fmt="(T3,I6,70F10.6)") i, q_eigvec(i, :)
2067 END DO
2068 DEALLOCATE (q_eigval)
2069 DEALLOCATE (q_eigvec)
2070 END IF ! debug
2071 END IF
2072 IF (debug) THEN
2073 ! Print the full atomic occupation matrix block
2074 ALLOCATE (q_work(nsgf_kind, nsgf_kind))
2075 q_work(:, :) = 0.0_dp
2076 IF (ASSOCIATED(q_block)) q_work(:, :) = q_block(:, :)
2077 CALL para_env%sum(q_work)
2078 IF (output_unit > 0) THEN
2079 norb = SIZE(q_work, 1)
2080 WRITE (unit=output_unit, fmt="(/,T9,200I10)") (i, i=1, norb)
2081 DO i = 1, norb
2082 WRITE (unit=output_unit, fmt="(T3,I6,200F10.6)") i, q_work(i, :)
2083 END DO
2084 ALLOCATE (q_eigval(norb))
2085 q_eigval(:) = 0.0_dp
2086 ALLOCATE (q_eigvec(norb, norb))
2087 q_eigvec(:, :) = 0.0_dp
2088 CALL jacobi(q_work, q_eigval, q_eigvec)
2089 WRITE (unit=output_unit, fmt="(/,T9,200I10)") (i, i=1, norb)
2090 WRITE (unit=output_unit, fmt="(T9,201F10.6)") (q_eigval(i), i=1, norb), &
2091 sum(q_eigval(1:norb))
2092 DO i = 1, norb
2093 WRITE (unit=output_unit, fmt="(T3,I6,200F10.6)") i, q_eigvec(i, :)
2094 END DO
2095 DEALLOCATE (q_eigval)
2096 DEALLOCATE (q_eigvec)
2097 END IF
2098 DEALLOCATE (q_work)
2099 END IF ! debug
2100 END IF ! should output
2101
2102 END DO ! next atom "iatom" of atomic kind "ikind"
2103
2104 IF (ALLOCATED(q_matrix)) THEN
2105 DEALLOCATE (q_matrix)
2106 END IF
2107
2108 END DO ! next atomic kind "ikind"
2109
2110 ! Add V(i,j)[U] to V(i,j)[DFT]
2111
2112 IF (ASSOCIATED(sm_h)) THEN
2113
2114 DO ikind = 1, nkind
2115
2116 IF (.NOT. is_plus_u_kind(ikind)) cycle
2117
2118 kind_a => atomic_kind_set(ikind)
2119
2120 CALL get_atomic_kind(atomic_kind=kind_a, &
2121 atom_list=atom_list, &
2122 natom=natom_of_kind)
2123
2124 DO iatom = 1, natom_of_kind
2125
2126 atom_a = atom_list(iatom)
2127
2128 CALL dbcsr_get_block_p(matrix=sm_h, &
2129 row=atom_a, &
2130 col=atom_a, &
2131 block=h_block, &
2132 found=found)
2133
2134 IF (.NOT. ASSOCIATED(h_block)) cycle
2135
2136 CALL dbcsr_get_block_p(matrix=sm_v, &
2137 row=atom_a, &
2138 col=atom_a, &
2139 block=v_block, &
2140 found=found)
2141 cpassert(ASSOCIATED(v_block))
2142
2143 IF (orthonormal_basis) THEN
2144 DO isgf = 1, SIZE(h_block, 1)
2145 h_block(isgf, isgf) = h_block(isgf, isgf) + v_block(isgf, isgf)
2146 END DO
2147 ELSE
2148 CALL dbcsr_get_block_p(matrix=sm_s, &
2149 row=atom_a, &
2150 col=atom_a, &
2151 block=s_block, &
2152 found=found)
2153 cpassert(ASSOCIATED(s_block))
2154 DO jsgf = 1, SIZE(h_block, 2)
2155 DO isgf = 1, SIZE(h_block, 1)
2156 h_block(isgf, jsgf) = h_block(isgf, jsgf) + v_block(isgf, jsgf)*s_block(isgf, jsgf)
2157 END DO
2158 END DO
2159 END IF ! orthonormal basis set
2160
2161 END DO ! next atom "iatom" of atomic kind "ikind"
2162
2163 END DO ! Next atomic kind "ikind"
2164
2165 END IF ! An update of the Hamiltonian matrix is requested
2166
2167 END DO ! next cell image
2168
2169 END DO ! next spin "ispin"
2170
2171 ! Collect the energy contributions from all processes
2172
2173 CALL para_env%sum(energy%dft_plus_u)
2174
2175 IF (energy%dft_plus_u < 0.0_dp) THEN
2176 IF (.NOT. occupation_enforced) THEN
2177 CALL cp_warn(__location__, &
2178 "DFT+U energy contribution is negative possibly due "// &
2179 "to unphysical Mulliken charges!")
2180 END IF
2181 END IF
2182
2183 CALL dbcsr_deallocate_matrix(sm_q)
2184 CALL dbcsr_deallocate_matrix(sm_v)
2185
2186 CALL timestop(handle)
2187
2188 END SUBROUTINE mulliken
2189
2190! **************************************************************************************************
2191!> \brief Add a DFT+U contribution to the Hamiltonian matrix\n
2192!> using a method based on Mulliken charges
2193!> \f[q_\mu = \sum_\nu \frac{1}{2}(P_{\mu\nu} S_{\nu\mu} +
2194!> S_{\mu\nu} P_{\nu\mu})
2195!> = \sum_\nu P_{\mu\nu} S_{\nu\mu}\f]
2196!> where \b P and \b S are the density and the
2197!> overlap matrix, respectively.
2198!> \param[in] qs_env Quickstep environment
2199!> \param orthonormal_basis ...
2200!> \param[in,out] matrix_h Hamiltonian matrices for each spin
2201!> \param[in,out] matrix_w Energy weighted density matrices for each spin
2202!> \param should_output ...
2203!> \param output_unit ...
2204!> \param print_level ...
2205!> \date 11.01.2008
2206!> \par
2207!> \f{eqnarray*}{
2208!> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}\\\
2209!> & = & E^{\rm DFT} + \frac{1}{2}(U - J)\sum_\mu (q_\mu - q_\mu^2)\\‍[1ex]
2210!> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
2211!> & = & \frac{\partial E^{\rm DFT}}
2212!> {\partial P_{\mu\nu}} +
2213!> \frac{\partial E^{\rm U}}
2214!> {\partial P_{\mu\nu}}\\\
2215!> & = & H_{\mu\nu} +
2216!> \frac{\partial E^{\rm U}}{\partial q_\mu}
2217!> \frac{\partial q_\mu}{\partial P_{\mu\nu}}\\\
2218!> & = & H_{\mu\nu} +
2219!> \frac{1}{2}(U - J)(1 - q_\mu - q_\nu) S_{\mu\nu}\\\
2220!> \f}
2221!> \author Matthias Krack (MK)
2222!> \version 1.0
2223!> \note The use of any full matrices was avoided. Thus no ScaLAPACK
2224!> calls are performed
2225! **************************************************************************************************
2226 SUBROUTINE mulliken_charges(qs_env, orthonormal_basis, matrix_h, matrix_w, &
2227 should_output, output_unit, print_level)
2228
2229 TYPE(qs_environment_type), POINTER :: qs_env
2230 LOGICAL, INTENT(IN) :: orthonormal_basis
2231 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
2232 POINTER :: matrix_h, matrix_w
2233 LOGICAL, INTENT(IN) :: should_output
2234 INTEGER, INTENT(IN) :: output_unit, print_level
2235
2236 CHARACTER(LEN=*), PARAMETER :: routinen = 'mulliken_charges'
2237
2238 CHARACTER(LEN=10) :: spin_info
2239 CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: symbol
2240 CHARACTER(LEN=default_string_length) :: atomic_kind_name
2241 INTEGER :: atom_a, handle, i, iatom, ic, ikind, isb, iset, isgf, ishell, ispin, jatom, jsgf, &
2242 lu, m, natom, natom_of_kind, nimg, nkind, nsb, nset, nsgf, nspin, sgf
2243 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
2244 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell
2245 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
2246 LOGICAL :: dft_plus_u_atom, found, just_energy
2247 REAL(kind=dp) :: eps_u_ramping, fspin, q, u_minus_j, &
2248 u_minus_j_target, u_ramping, v
2249 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dedq, trps
2250 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: q_ii
2251 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, p_block, s_block, w_block
2252 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2253 TYPE(dbcsr_iterator_type) :: iter
2254 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
2255 TYPE(dbcsr_type), POINTER :: sm_h, sm_p, sm_s, sm_w
2256 TYPE(dft_control_type), POINTER :: dft_control
2257 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
2258 TYPE(mp_para_env_type), POINTER :: para_env
2259 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2260 TYPE(qs_energy_type), POINTER :: energy
2261 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2262 TYPE(qs_rho_type), POINTER :: rho
2263
2264 CALL timeset(routinen, handle)
2265
2266 NULLIFY (atom_list)
2267 NULLIFY (atomic_kind_set)
2268 NULLIFY (qs_kind_set)
2269 NULLIFY (dft_control)
2270 NULLIFY (energy)
2271 NULLIFY (first_sgf)
2272 NULLIFY (h_block)
2273 NULLIFY (matrix_p)
2274 NULLIFY (matrix_s)
2275 NULLIFY (l)
2276 NULLIFY (last_sgf)
2277 NULLIFY (nshell)
2278 NULLIFY (orb_basis_set)
2279 NULLIFY (p_block)
2280 NULLIFY (particle_set)
2281 NULLIFY (rho)
2282 NULLIFY (s_block)
2283 NULLIFY (sm_h)
2284 NULLIFY (sm_p)
2285 NULLIFY (sm_s)
2286 NULLIFY (w_block)
2287 NULLIFY (para_env)
2288
2289 CALL get_qs_env(qs_env=qs_env, &
2290 atomic_kind_set=atomic_kind_set, &
2291 qs_kind_set=qs_kind_set, &
2292 dft_control=dft_control, &
2293 energy=energy, &
2294 particle_set=particle_set, &
2295 rho=rho, &
2296 para_env=para_env)
2297
2298 cpassert(ASSOCIATED(atomic_kind_set))
2299 cpassert(ASSOCIATED(dft_control))
2300 cpassert(ASSOCIATED(energy))
2301 cpassert(ASSOCIATED(particle_set))
2302 cpassert(ASSOCIATED(rho))
2303
2304 IF (orthonormal_basis) THEN
2305 NULLIFY (sm_s)
2306 ELSE
2307 ! Get overlap matrix in sparse format
2308 CALL get_qs_env(qs_env=qs_env, &
2309 matrix_s_kp=matrix_s)
2310 cpassert(ASSOCIATED(matrix_s))
2311 END IF
2312
2313 ! Get density matrices in sparse format
2314
2315 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2316
2317 energy%dft_plus_u = 0.0_dp
2318
2319 nspin = dft_control%nspins
2320 nimg = dft_control%nimages
2321
2322 IF (nspin == 2) THEN
2323 fspin = 1.0_dp
2324 ELSE
2325 fspin = 0.5_dp
2326 END IF
2327
2328 ! Get the total number of atoms, contracted spherical Gaussian basis
2329 ! functions, and atomic kinds
2330
2331 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
2332 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
2333
2334 nkind = SIZE(atomic_kind_set)
2335
2336 ALLOCATE (first_sgf_atom(natom))
2337 first_sgf_atom(:) = 0
2338
2339 CALL get_particle_set(particle_set, qs_kind_set, &
2340 first_sgf=first_sgf_atom)
2341
2342 ALLOCATE (trps(nsgf))
2343 trps(:) = 0.0_dp
2344
2345 IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
2346 ALLOCATE (dedq(nsgf))
2347 just_energy = .false.
2348 ELSE
2349 just_energy = .true.
2350 END IF
2351
2352 ! Loop over all spins
2353
2354 DO ispin = 1, nspin
2355
2356 IF (.NOT. just_energy) dedq(:) = 0.0_dp
2357
2358 ! Calculate Trace(P*S) assuming symmetric matrices
2359
2360 trps(:) = 0.0_dp
2361
2362 DO ic = 1, nimg
2363 IF (orthonormal_basis) THEN
2364 NULLIFY (sm_s)
2365 ELSE
2366 sm_s => matrix_s(1, ic)%matrix
2367 END IF
2368 sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
2369
2370 CALL dbcsr_iterator_start(iter, sm_p)
2371
2372 DO WHILE (dbcsr_iterator_blocks_left(iter))
2373
2374 CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
2375
2376 IF (orthonormal_basis) THEN
2377
2378 IF (iatom /= jatom) cycle
2379
2380 IF (ASSOCIATED(p_block)) THEN
2381 sgf = first_sgf_atom(iatom)
2382 DO isgf = 1, SIZE(p_block, 1)
2383 trps(sgf) = trps(sgf) + p_block(isgf, isgf)
2384 sgf = sgf + 1
2385 END DO
2386 END IF
2387
2388 ELSE
2389
2390 CALL dbcsr_get_block_p(matrix=sm_s, &
2391 row=iatom, &
2392 col=jatom, &
2393 block=s_block, &
2394 found=found)
2395 cpassert(ASSOCIATED(s_block))
2396
2397 sgf = first_sgf_atom(jatom)
2398 DO jsgf = 1, SIZE(p_block, 2)
2399 DO isgf = 1, SIZE(p_block, 1)
2400 trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
2401 END DO
2402 sgf = sgf + 1
2403 END DO
2404
2405 IF (iatom /= jatom) THEN
2406 sgf = first_sgf_atom(iatom)
2407 DO isgf = 1, SIZE(p_block, 1)
2408 DO jsgf = 1, SIZE(p_block, 2)
2409 trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
2410 END DO
2411 sgf = sgf + 1
2412 END DO
2413 END IF
2414
2415 END IF ! orthonormal basis set
2416
2417 END DO ! next atom "iatom"
2418
2419 CALL dbcsr_iterator_stop(iter)
2420
2421 END DO ! cell images
2422
2423 CALL para_env%sum(trps)
2424
2425 ! q <- Trace(PS)
2426
2427 ! E[DFT+U] = E[DFT] + E[U]
2428 ! = E[DFT] + (U - J)*(q - q**2))/2
2429
2430 ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U]
2431 ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j)
2432 ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j))
2433
2434 ! Loop over all atomic kinds
2435
2436 DO ikind = 1, nkind
2437
2438 ! Load the required atomic kind data
2439 CALL get_atomic_kind(atomic_kind_set(ikind), &
2440 atom_list=atom_list, &
2441 name=atomic_kind_name, &
2442 natom=natom_of_kind)
2443
2444 CALL get_qs_kind(qs_kind_set(ikind), &
2445 dft_plus_u_atom=dft_plus_u_atom, &
2446 l_of_dft_plus_u=lu, &
2447 basis_set=orb_basis_set, &
2448 u_minus_j=u_minus_j, &
2449 u_minus_j_target=u_minus_j_target, &
2450 u_ramping=u_ramping, &
2451 eps_u_ramping=eps_u_ramping)
2452
2453 ! Check, if this atom needs a DFT+U correction
2454
2455 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
2456 IF (.NOT. dft_plus_u_atom) cycle
2457 IF (lu < 0) cycle
2458
2459 ! Apply U ramping if requested
2460
2461 IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
2462 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
2463 u_minus_j = min(u_minus_j + u_ramping, u_minus_j_target)
2464 CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
2465 END IF
2466 IF (should_output .AND. (output_unit > 0)) THEN
2467 WRITE (unit=output_unit, fmt="(T3,A,3X,A,F0.3,A)") &
2468 "Kind name: "//trim(adjustl(atomic_kind_name)), &
2469 "U(eff) = ", u_minus_j*evolt, " eV"
2470 END IF
2471 END IF
2472
2473 IF (u_minus_j == 0.0_dp) cycle
2474
2475 ! Load the required Gaussian basis set data
2476
2477 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2478 first_sgf=first_sgf, &
2479 l=l, &
2480 last_sgf=last_sgf, &
2481 nset=nset, &
2482 nshell=nshell)
2483
2484 ! Count the relevant shell blocks of this atomic kind
2485
2486 nsb = 0
2487 DO iset = 1, nset
2488 DO ishell = 1, nshell(iset)
2489 IF (l(ishell, iset) == lu) nsb = nsb + 1
2490 END DO
2491 END DO
2492
2493 ALLOCATE (q_ii(nsb, 2*lu + 1))
2494
2495 ! Print headline if requested
2496
2497 IF (should_output .AND. (print_level > low_print_level)) THEN
2498 IF (output_unit > 0) THEN
2499 ALLOCATE (symbol(2*lu + 1))
2500 DO m = -lu, lu
2501 symbol(lu + m + 1) = sgf_symbol(0, lu, m)
2502 END DO
2503 IF (nspin > 1) THEN
2504 WRITE (unit=spin_info, fmt="(A8,I2)") " of spin", ispin
2505 ELSE
2506 spin_info = ""
2507 END IF
2508 WRITE (unit=output_unit, fmt="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
2509 "DFT+U occupations"//trim(spin_info)//" for the atoms of atomic kind ", ikind, &
2510 ": "//trim(atomic_kind_name), &
2511 "Atom Shell ", (adjustr(symbol(i)), i=1, 2*lu + 1), " Trace"
2512 DEALLOCATE (symbol)
2513 END IF
2514 END IF
2515
2516 ! Loop over all atoms of the current atomic kind
2517
2518 DO iatom = 1, natom_of_kind
2519
2520 atom_a = atom_list(iatom)
2521
2522 q_ii(:, :) = 0.0_dp
2523
2524 ! Get diagonal block
2525
2526 CALL dbcsr_get_block_p(matrix=sm_p, &
2527 row=atom_a, &
2528 col=atom_a, &
2529 block=p_block, &
2530 found=found)
2531
2532 ! Calculate E(U) and dE(U)/dq
2533
2534 IF (ASSOCIATED(p_block)) THEN
2535
2536 sgf = first_sgf_atom(atom_a)
2537
2538 isb = 0
2539 DO iset = 1, nset
2540 DO ishell = 1, nshell(iset)
2541 IF (l(ishell, iset) == lu) THEN
2542 isb = isb + 1
2543 i = 0
2544 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
2545 q = fspin*trps(sgf)
2546 i = i + 1
2547 q_ii(isb, i) = q
2548 energy%dft_plus_u = energy%dft_plus_u + &
2549 0.5_dp*u_minus_j*(q - q**2)/fspin
2550 IF (.NOT. just_energy) THEN
2551 dedq(sgf) = dedq(sgf) + u_minus_j*(0.5_dp - q)
2552 END IF
2553 sgf = sgf + 1
2554 END DO ! next contracted spherical Gaussian function "isgf"
2555 ELSE
2556 sgf = sgf + last_sgf(ishell, iset) - first_sgf(ishell, iset) + 1
2557 END IF ! angular momentum requested for DFT+U correction
2558 END DO ! next shell "ishell"
2559 END DO ! next shell set "iset"
2560
2561 END IF ! this process is the owner of the sparse matrix block?
2562
2563 ! Consider print requests
2564
2565 IF (should_output .AND. (print_level > low_print_level)) THEN
2566 CALL para_env%sum(q_ii)
2567 IF (output_unit > 0) THEN
2568 DO isb = 1, nsb
2569 WRITE (unit=output_unit, fmt="(T3,I6,2X,I6,2X,10F8.3)") &
2570 atom_a, isb, q_ii(isb, :), sum(q_ii(isb, :))
2571 END DO
2572 WRITE (unit=output_unit, fmt="(T12,A,2X,10F8.3)") &
2573 "Total", (sum(q_ii(:, i)), i=1, 2*lu + 1), sum(q_ii)
2574 WRITE (unit=output_unit, fmt="(A)") ""
2575 END IF
2576 END IF ! should output
2577
2578 END DO ! next atom "iatom" of atomic kind "ikind"
2579
2580 IF (ALLOCATED(q_ii)) THEN
2581 DEALLOCATE (q_ii)
2582 END IF
2583
2584 END DO ! next atomic kind "ikind"
2585
2586 IF (.NOT. just_energy) THEN
2587 CALL para_env%sum(dedq)
2588 END IF
2589
2590 ! Add V(i,j)[U] to V(i,j)[DFT]
2591
2592 IF (PRESENT(matrix_h)) THEN
2593
2594 DO ic = 1, nimg
2595 IF (orthonormal_basis) THEN
2596 NULLIFY (sm_s)
2597 ELSE
2598 sm_s => matrix_s(1, ic)%matrix
2599 END IF
2600 sm_h => matrix_h(ispin, ic)%matrix
2601
2602 CALL dbcsr_iterator_start(iter, sm_h)
2603
2604 DO WHILE (dbcsr_iterator_blocks_left(iter))
2605
2606 CALL dbcsr_iterator_next_block(iter, iatom, jatom, h_block)
2607
2608 IF (orthonormal_basis) THEN
2609
2610 IF (iatom /= jatom) cycle
2611
2612 IF (ASSOCIATED(h_block)) THEN
2613 sgf = first_sgf_atom(iatom)
2614 DO isgf = 1, SIZE(h_block, 1)
2615 h_block(isgf, isgf) = h_block(isgf, isgf) + dedq(sgf)
2616 sgf = sgf + 1
2617 END DO
2618 END IF
2619
2620 ELSE
2621
2622 ! Request katom just to check for consistent sparse matrix pattern
2623
2624 CALL dbcsr_get_block_p(matrix=sm_s, &
2625 row=iatom, &
2626 col=jatom, &
2627 block=s_block, &
2628 found=found)
2629 cpassert(ASSOCIATED(s_block))
2630
2631 ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
2632
2633 sgf = first_sgf_atom(iatom)
2634
2635 DO isgf = 1, SIZE(h_block, 1)
2636 IF (dedq(sgf) /= 0.0_dp) THEN
2637 v = 0.5_dp*dedq(sgf)
2638 DO jsgf = 1, SIZE(h_block, 2)
2639 h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
2640 END DO
2641 END IF
2642 sgf = sgf + 1
2643 END DO
2644
2645 sgf = first_sgf_atom(jatom)
2646
2647 DO jsgf = 1, SIZE(h_block, 2)
2648 IF (dedq(sgf) /= 0.0_dp) THEN
2649 v = 0.5_dp*dedq(sgf)
2650 DO isgf = 1, SIZE(h_block, 1)
2651 h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
2652 END DO
2653 END IF
2654 sgf = sgf + 1
2655 END DO
2656
2657 END IF ! orthonormal basis set
2658
2659 END DO ! Next atom "iatom"
2660
2661 CALL dbcsr_iterator_stop(iter)
2662
2663 END DO
2664
2665 END IF ! An update of the Hamiltonian matrix is requested
2666
2667 ! Calculate the contribution (non-Pulay part) to the derivatives
2668 ! w.r.t. the nuclear positions, which requires an update of the
2669 ! energy weighted density W.
2670
2671 IF (PRESENT(matrix_w) .AND. (.NOT. orthonormal_basis)) THEN
2672
2673 DO ic = 1, nimg
2674 sm_s => matrix_s(1, ic)%matrix
2675 sm_p => matrix_p(ispin, ic)%matrix
2676 sm_w => matrix_w(ispin, ic)%matrix
2677
2678 CALL dbcsr_iterator_start(iter, sm_p)
2679
2680 DO WHILE (dbcsr_iterator_blocks_left(iter))
2681
2682 CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
2683
2684 ! Skip the diagonal blocks of the W matrix
2685
2686 IF (iatom == jatom) cycle
2687
2688 ! Request katom just to check for consistent sparse matrix patterns
2689
2690 CALL dbcsr_get_block_p(matrix=sm_w, &
2691 row=iatom, &
2692 col=jatom, &
2693 block=w_block, &
2694 found=found)
2695 cpassert(ASSOCIATED(w_block))
2696
2697 ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
2698
2699 sgf = first_sgf_atom(iatom)
2700
2701 DO isgf = 1, SIZE(w_block, 1)
2702 IF (dedq(sgf) /= 0.0_dp) THEN
2703 v = -0.5_dp*dedq(sgf)
2704 DO jsgf = 1, SIZE(w_block, 2)
2705 w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
2706 END DO
2707 END IF
2708 sgf = sgf + 1
2709 END DO
2710
2711 sgf = first_sgf_atom(jatom)
2712
2713 DO jsgf = 1, SIZE(w_block, 2)
2714 IF (dedq(sgf) /= 0.0_dp) THEN
2715 v = -0.5_dp*dedq(sgf)
2716 DO isgf = 1, SIZE(w_block, 1)
2717 w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
2718 END DO
2719 END IF
2720 sgf = sgf + 1
2721 END DO
2722
2723 END DO ! next block node "jatom"
2724
2725 CALL dbcsr_iterator_stop(iter)
2726
2727 END DO
2728
2729 END IF ! W matrix update requested
2730
2731 END DO ! next spin "ispin"
2732
2733 ! Collect the energy contributions from all processes
2734
2735 CALL para_env%sum(energy%dft_plus_u)
2736
2737 IF (energy%dft_plus_u < 0.0_dp) THEN
2738 CALL cp_warn(__location__, &
2739 "DFT+U energy contribution is negative possibly due "// &
2740 "to unphysical Mulliken charges!")
2741 END IF
2742
2743 ! Release local work storage
2744
2745 IF (ALLOCATED(first_sgf_atom)) THEN
2746 DEALLOCATE (first_sgf_atom)
2747 END IF
2748
2749 IF (ALLOCATED(trps)) THEN
2750 DEALLOCATE (trps)
2751 END IF
2752
2753 IF (ALLOCATED(dedq)) THEN
2754 DEALLOCATE (dedq)
2755 END IF
2756
2757 CALL timestop(handle)
2758
2759 END SUBROUTINE mulliken_charges
2760
2761! **************************************************************************************************
2762!> \brief Calculate the tensorial DFT+U+J energy contribution, Hamiltonian correction,
2763!> and optional force contribution.
2764!> \param[in] qs_env Quickstep environment
2765!> \param[in,out] matrix_h Hamiltonian matrices
2766!> \param[in,out] matrix_w Energy-weighted density matrices for each spin channel
2767!> \param[in] should_output Whether detailed output should be printed
2768!> \param[in] output_unit ...
2769!> \param[in] print_level ...
2770!> \date 14.03.2026
2771!> \author Ziwei Chai (ZC)
2772! **************************************************************************************************
2773 SUBROUTINE tensorial(qs_env, matrix_h, matrix_w, should_output, output_unit, print_level)
2774
2775 TYPE(qs_environment_type), POINTER :: qs_env
2776 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
2777 POINTER :: matrix_h, matrix_w
2778 LOGICAL, INTENT(IN) :: should_output
2779 INTEGER, INTENT(IN) :: output_unit, print_level
2780
2781 CHARACTER(LEN=*), PARAMETER :: routinen = 'tensorial'
2782
2783 CHARACTER(LEN=10) :: spin_info
2784 CHARACTER(LEN=6), ALLOCATABLE :: symbol(:)
2785 CHARACTER(LEN=default_string_length) :: atomic_kind_name
2786 INTEGER :: atom_a, atom_aa, handle, i, i0, iatom, ikind, iorb, isb, iset, isgf, ishell, &
2787 ispin, j, jatom, jj, jkind, jsb, jset, jsgf, jshell, katom, lu, m, matom_of_kind, &
2788 max_scf, n, natom, natom_of_kind, nkind, norb, nsb, nsbsize, nset, nsgf, nsgf_kind, &
2789 nsgf_kind2, nspin, timevalues(8)
2790 INTEGER, ALLOCATABLE :: first_sgf_atom(:)
2791 INTEGER, DIMENSION(1) :: iloc
2792 INTEGER, DIMENSION(:), POINTER :: atom_list, atom_list2, nshell, orbitals
2793 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
2794 LOGICAL :: debug, dft_plus_u_atom, found, &
2795 just_energy, smear
2796 LOGICAL, ALLOCATABLE :: orb_occ(:)
2797 REAL(kind=dp) :: eps_scf, eps_u_ramping, fspin, hund_j, occ, perturbation_strength, &
2798 some_real, trq, trq2, trqxq_hund_j, u_minus_j, u_minus_j_target, u_ramping
2799 REAL(kind=dp), ALLOCATABLE :: a(:), q_a_matrix(:, :), q_b_matrix(:, :), q_eigval(:), &
2800 q_eigvec(:, :), q_work(:, :), s_inv_matrix(:, :), s_matrix(:, :), tmp_matrix1(:, :), &
2801 tmp_matrix2(:, :), tmp_matrix3(:, :), v_matrix(:, :), vhxc_matrix(:, :)
2802
2803 TYPE :: block3d
2804 REAL(kind=dp), ALLOCATABLE :: matrix(:, :, :)
2805 END TYPE block3d
2806 TYPE :: block4d
2807 REAL(kind=dp), ALLOCATABLE :: matrix(:, :, :, :)
2808 END TYPE block4d
2809 TYPE(block3d), ALLOCATABLE :: a_matrix(:, :), d_matrix(:, :), &
2810 e_matrix(:, :), i_matrix(:, :), &
2811 j_matrix(:, :)
2812 TYPE(block4d), ALLOCATABLE :: dd_matrix(:, :), ee_matrix(:, :), &
2813 ff_matrix(:, :)
2814 REAL(kind=dp), DIMENSION(:), POINTER :: ao_coef
2815 REAL(kind=dp), DIMENSION(:, :), POINTER :: a_block, d_block, e_block, f_block, &
2816 i_block, j_block, q_a_block, q_b_block, &
2817 s_block, tmp_block, v_block, vhxc_block
2818
2819 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2820 TYPE(mp_para_env_type), POINTER :: para_env
2821 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s, matrix_vhxc
2822 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_s
2823 TYPE(dbcsr_type), POINTER :: sm_a, sm_b, sm_c, sm_d, sm_e, sm_f, &
2824 sm_h, sm_i, sm_j, sm_k, sm_q_a, sm_q_b, &
2825 sm_s, sm_v, sm_w, sm_mid, sm_vhxc, sm_p
2826 TYPE(dbcsr_iterator_type) :: iter
2827 TYPE(dft_control_type), POINTER :: dft_control
2828 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
2829 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2830 TYPE(qs_energy_type), POINTER :: energy
2831 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2832 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2833 TYPE(qs_rho_type), POINTER :: rho
2834 TYPE(qs_scf_env_type), POINTER :: scf_env
2835
2836 CALL timeset(routinen, handle)
2837
2838 debug = .false.
2839 smear = .false.
2840 max_scf = -1
2841 eps_scf = 1.0e30_dp
2842
2843 NULLIFY (atom_list)
2844 NULLIFY (atom_list2)
2845 NULLIFY (atomic_kind_set)
2846 NULLIFY (qs_kind_set)
2847 NULLIFY (dft_control)
2848 NULLIFY (energy)
2849 NULLIFY (force)
2850 NULLIFY (first_sgf)
2851 NULLIFY (matrix_p)
2852 NULLIFY (matrix_s)
2853 NULLIFY (matrix_vhxc)
2854 NULLIFY (l)
2855 NULLIFY (last_sgf)
2856 NULLIFY (nshell)
2857 NULLIFY (orb_basis_set)
2858 NULLIFY (orbitals)
2859 NULLIFY (ao_coef)
2860 NULLIFY (particle_set)
2861 NULLIFY (q_a_block)
2862 NULLIFY (q_b_block)
2863 NULLIFY (vhxc_block)
2864 NULLIFY (rho)
2865 NULLIFY (scf_env)
2866 NULLIFY (sm_h)
2867 NULLIFY (sm_p)
2868 NULLIFY (sm_q_a)
2869 NULLIFY (sm_q_b)
2870 NULLIFY (sm_s)
2871 NULLIFY (sm_v)
2872 NULLIFY (sm_a)
2873 NULLIFY (sm_b)
2874 NULLIFY (sm_c)
2875 NULLIFY (sm_d)
2876 NULLIFY (sm_e)
2877 NULLIFY (sm_f)
2878 NULLIFY (sm_i)
2879 NULLIFY (sm_j)
2880 NULLIFY (sm_k)
2881 NULLIFY (sm_mid)
2882 NULLIFY (sm_vhxc)
2883 NULLIFY (v_block)
2884 NULLIFY (a_block)
2885 NULLIFY (d_block)
2886 NULLIFY (e_block)
2887 NULLIFY (f_block)
2888 NULLIFY (i_block)
2889 NULLIFY (j_block)
2890 NULLIFY (para_env)
2891 NULLIFY (s_block)
2892 NULLIFY (tmp_block)
2893
2894 CALL get_qs_env(qs_env=qs_env, &
2895 atomic_kind_set=atomic_kind_set, &
2896 qs_kind_set=qs_kind_set, &
2897 dft_control=dft_control, &
2898 energy=energy, &
2899 matrix_s=matrix_s, &
2900 matrix_vhxc=matrix_vhxc, &
2901 matrix_s_kp=matrixkp_s, &
2902 particle_set=particle_set, &
2903 rho=rho, &
2904 scf_env=scf_env, &
2905 para_env=para_env)
2906
2907 cpassert(ASSOCIATED(atomic_kind_set))
2908 cpassert(ASSOCIATED(dft_control))
2909 cpassert(ASSOCIATED(energy))
2910 cpassert(ASSOCIATED(matrix_s))
2911 cpassert(ASSOCIATED(particle_set))
2912 cpassert(ASSOCIATED(rho))
2913
2914 sm_s => matrix_s(1)%matrix
2915 CALL qs_rho_get(rho, rho_ao=matrix_p)
2916
2917 energy%dft_plus_u = 0.0_dp
2918
2919 nspin = dft_control%nspins
2920 IF (dft_control%mtlr_dft_with_perturbation) THEN
2921 IF (.NOT. ASSOCIATED(matrix_vhxc)) THEN
2922 cpabort("MTLR requires the projected Hxc matrix for every spin channel.")
2923 END IF
2924 IF (SIZE(matrix_vhxc) /= nspin) THEN
2925 cpabort("The number of projected Hxc matrices does not match the spin channels.")
2926 END IF
2927 END IF
2928
2929 IF (nspin == 2) THEN
2930 fspin = 1.0_dp
2931 ELSE
2932 fspin = 0.5_dp
2933 END IF
2934
2935 ! Get the total number of atoms, contracted spherical Gaussian basis
2936 ! functions, and atomic kinds
2937
2938 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
2939 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
2940
2941 nkind = SIZE(atomic_kind_set)
2942
2943 ALLOCATE (first_sgf_atom(natom))
2944 first_sgf_atom(:) = 0
2945
2946 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf_atom)
2947
2948 IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
2949 just_energy = .false.
2950 ELSE
2951 just_energy = .true.
2952 END IF
2953
2954 CALL dbcsr_init_p(sm_mid)
2955 CALL dbcsr_create(sm_mid, template=sm_s, matrix_type="N")
2956
2957 ! Create local block diagonal matrices
2958
2959 ALLOCATE (sm_q_a)
2960 CALL dbcsr_get_block_diag(sm_s, sm_q_a)
2961
2962 ALLOCATE (sm_q_b)
2963 CALL dbcsr_get_block_diag(sm_s, sm_q_b)
2964
2965 CALL dbcsr_set(sm_q_a, 0.0_dp)
2966 CALL dbcsr_set(sm_q_b, 0.0_dp)
2967
2968 ! Loop over all spins for calculating sm_q_a and sm_q_b
2969
2970 DO ispin = 1, nspin
2971
2972 sm_p => matrix_p(ispin)%matrix
2973
2974 IF (nspin == 2) THEN
2975 IF (ispin == 1) THEN
2976 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2977 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_a)
2978 ELSE
2979 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2980 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_b)
2981 END IF
2982 ELSE
2983 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2984 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_a)
2985 CALL dbcsr_copy(sm_q_b, sm_q_a)
2986 END IF
2987
2988 END DO
2989
2990 ALLOCATE (sm_v)
2991 CALL dbcsr_get_block_diag(sm_s, sm_v)
2992
2993 IF (PRESENT(matrix_w)) THEN
2994
2995 CALL dbcsr_init_p(sm_b)
2996 CALL dbcsr_create(sm_b, template=matrix_s(1)%matrix, matrix_type="N")
2997
2998 CALL dbcsr_init_p(sm_d)
2999 CALL dbcsr_create(sm_d, template=matrix_s(1)%matrix, matrix_type="N")
3000
3001 CALL dbcsr_init_p(sm_c)
3002 CALL dbcsr_create(sm_c, template=matrix_s(1)%matrix, matrix_type="N")
3003
3004 CALL dbcsr_init_p(sm_e)
3005 CALL dbcsr_create(sm_e, template=matrix_s(1)%matrix, matrix_type="N")
3006
3007 CALL dbcsr_init_p(sm_f)
3008 CALL dbcsr_create(sm_f, template=matrix_s(1)%matrix, matrix_type="N")
3009
3010 CALL dbcsr_init_p(sm_i)
3011 CALL dbcsr_create(sm_i, template=matrix_s(1)%matrix, matrix_type="N")
3012
3013 CALL dbcsr_init_p(sm_j)
3014 CALL dbcsr_create(sm_j, template=matrix_s(1)%matrix, matrix_type="N")
3015
3016 CALL dbcsr_init_p(sm_k)
3017 CALL dbcsr_create(sm_k, template=matrix_s(1)%matrix, matrix_type="N")
3018
3019 END IF
3020
3021 ALLOCATE (sm_a)
3022 CALL dbcsr_get_block_diag(sm_s, sm_a)
3023
3024 ! Loop over all spins
3025
3026 DO ispin = 1, nspin
3027
3028 IF (PRESENT(matrix_h)) THEN
3029 sm_h => matrix_h(ispin, 1)%matrix
3030 ELSE
3031 NULLIFY (sm_h)
3032 END IF
3033
3034 IF (PRESENT(matrix_w)) THEN
3035 sm_w => matrix_w(ispin, 1)%matrix
3036 ELSE
3037 NULLIFY (sm_w)
3038 END IF
3039
3040 CALL dbcsr_set(sm_v, 0.0_dp)
3041 CALL dbcsr_set(sm_a, 0.0_dp)
3042
3043 IF (dft_control%mtlr_dft_with_perturbation) THEN
3044 IF (ispin == 1) perturbation_strength = 0.9_dp*dft_control%perturbation_strength
3045 IF (ispin == 2) perturbation_strength = 1.1_dp*dft_control%perturbation_strength
3046 IF (.NOT. ASSOCIATED(matrix_vhxc(ispin)%matrix)) THEN
3047 cpabort("MTLR projected Hxc matrix is not initialized.")
3048 END IF
3049 sm_vhxc => matrix_vhxc(ispin)%matrix
3050 END IF
3051
3052 ! Loop over all atomic kinds
3053
3054 DO ikind = 1, nkind
3055
3056 CALL get_atomic_kind(atomic_kind_set(ikind), &
3057 atom_list=atom_list, &
3058 name=atomic_kind_name, &
3059 natom=natom_of_kind)
3060
3061 CALL get_qs_kind(qs_kind_set(ikind), &
3062 dft_plus_u_atom=dft_plus_u_atom, &
3063 l_of_dft_plus_u=lu, &
3064 nsgf=nsgf_kind, &
3065 basis_set=orb_basis_set, &
3066 u_minus_j=u_minus_j, &
3067 hund_j=hund_j, &
3068 ao_coef=ao_coef, &
3069 u_minus_j_target=u_minus_j_target, &
3070 u_ramping=u_ramping, &
3071 eps_u_ramping=eps_u_ramping, &
3072 orbitals=orbitals, &
3073 eps_scf=eps_scf, &
3074 max_scf=max_scf, &
3075 smear=smear)
3076
3077 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
3078 IF (.NOT. dft_plus_u_atom) cycle
3079 IF (lu < 0) cycle
3080 IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) cycle
3081 IF (.NOT. ASSOCIATED(ao_coef)) THEN
3082 cpabort("Tensorial DFT+U requires associated ao_coef")
3083 END IF
3084
3085 ALLOCATE (a(SIZE(ao_coef)))
3086 a(:) = ao_coef(:)
3087
3088 ! Apply U ramping if requested
3089
3090 IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
3091 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
3092 u_minus_j = min(u_minus_j + u_ramping, u_minus_j_target)
3093 CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
3094 END IF
3095
3096 IF (should_output .AND. (output_unit > 0)) THEN
3097 WRITE (unit=output_unit, fmt="(T3,A,3X,A,F0.3,A)") &
3098 "Kind name: "//trim(adjustl(atomic_kind_name)), &
3099 "U(eff) = ", u_minus_j*evolt, " eV"
3100 END IF
3101 END IF
3102
3103 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
3104 first_sgf=first_sgf, &
3105 l=l, &
3106 last_sgf=last_sgf, &
3107 nset=nset, &
3108 nshell=nshell)
3109
3110 nsb = 1
3111 nsbsize = (2*lu + 1)
3112 n = nsb*nsbsize
3113
3114 ALLOCATE (q_a_matrix(n, n))
3115 ALLOCATE (q_b_matrix(n, n))
3116 ALLOCATE (s_matrix(n, n))
3117 ALLOCATE (s_inv_matrix(n, n))
3118 ALLOCATE (v_matrix(n, n))
3119
3120 q_a_matrix(:, :) = 0.0_dp
3121 q_b_matrix(:, :) = 0.0_dp
3122 s_matrix(:, :) = 0.0_dp
3123 s_inv_matrix(:, :) = 0.0_dp
3124 v_matrix(:, :) = 0.0_dp
3125
3126 ! Print headline if requested
3127
3128 IF (should_output .AND. (print_level > low_print_level)) THEN
3129 IF (output_unit > 0) THEN
3130 ALLOCATE (symbol(nsbsize))
3131 DO m = -lu, lu
3132 symbol(lu + m + 1) = sgf_symbol(0, lu, m)
3133 END DO
3134
3135 IF (nspin > 1) THEN
3136 WRITE (unit=spin_info, fmt="(A8,I2)") " of spin", ispin
3137 ELSE
3138 spin_info = ""
3139 END IF
3140
3141 WRITE (unit=output_unit, fmt="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
3142 "DFT+U occupations"//trim(spin_info)//" for the atoms of atomic kind ", ikind, &
3143 ": "//trim(atomic_kind_name), &
3144 "Atom Shell ", (adjustr(symbol(i)), i=1, nsbsize), " Trace"
3145
3146 DEALLOCATE (symbol)
3147 END IF
3148 END IF
3149
3150 IF (dft_control%mtlr_dft_with_perturbation .AND. &
3151 dft_control%mtlr_ikind == ikind) THEN
3152 dft_control%vhxc(ispin) = 0.0_dp
3153 dft_control%trq(ispin) = 0.0_dp
3154 ALLOCATE (vhxc_matrix(n, n))
3155 vhxc_matrix(:, :) = 0.0_dp
3156 END IF
3157
3158 ! Loop over all atoms of the current atomic kind
3159
3160 DO iatom = 1, natom_of_kind
3161
3162 atom_a = atom_list(iatom)
3163
3164 q_a_matrix(:, :) = 0.0_dp
3165 q_b_matrix(:, :) = 0.0_dp
3166 s_matrix(:, :) = 0.0_dp
3167 s_inv_matrix(:, :) = 0.0_dp
3168 v_matrix(:, :) = 0.0_dp
3169
3170 CALL dbcsr_get_block_p(matrix=sm_q_a, &
3171 row=atom_a, &
3172 col=atom_a, &
3173 block=q_a_block, &
3174 found=found)
3175
3176 CALL dbcsr_get_block_p(matrix=sm_q_b, &
3177 row=atom_a, &
3178 col=atom_a, &
3179 block=q_b_block, &
3180 found=found)
3181
3182 CALL dbcsr_get_block_p(matrix=sm_s, &
3183 row=atom_a, &
3184 col=atom_a, &
3185 block=s_block, &
3186 found=found)
3187
3188 IF (dft_control%mtlr_dft_with_perturbation .AND. &
3189 dft_control%mtlr_ikind == ikind .AND. &
3190 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
3191 CALL dbcsr_get_block_p(matrix=sm_vhxc, &
3192 row=atom_a, &
3193 col=atom_a, &
3194 block=vhxc_block, &
3195 found=found)
3196 END IF
3197
3198 IF (ASSOCIATED(q_a_block) .AND. ASSOCIATED(q_b_block) .AND. ASSOCIATED(s_block)) THEN
3199
3200 isb = 0
3201 DO iset = 1, nset
3202 DO ishell = 1, nshell(iset)
3203 IF (l(ishell, iset) /= lu) cycle
3204 isb = isb + 1
3205 i = 0
3206 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3207 i = i + 1
3208 jsb = 0
3209 DO jset = 1, nset
3210 DO jshell = 1, nshell(jset)
3211 IF (l(jshell, jset) /= lu) cycle
3212 jsb = jsb + 1
3213 j = 0
3214 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
3215 j = j + 1
3216 q_a_matrix(i, j) = q_a_matrix(i, j) + q_a_block(isgf, jsgf)*a(isb)*a(jsb)
3217 q_b_matrix(i, j) = q_b_matrix(i, j) + q_b_block(isgf, jsgf)*a(isb)*a(jsb)
3218 s_matrix(i, j) = s_matrix(i, j) + s_block(isgf, jsgf)*a(isb)*a(jsb)
3219 IF (ASSOCIATED(vhxc_block) .AND. &
3220 dft_control%mtlr_dft_with_perturbation .AND. &
3221 dft_control%mtlr_ikind == ikind .AND. &
3222 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
3223 vhxc_matrix(i, j) = vhxc_matrix(i, j) + vhxc_block(isgf, jsgf)*a(isb)*a(jsb)
3224 END IF
3225 END DO
3226 END DO
3227 END DO
3228 END DO
3229 END DO
3230 END DO
3231
3232 CALL invert_matrix(s_matrix, s_inv_matrix, some_real)
3233
3234 q_a_matrix(:, :) = matmul(q_a_matrix, s_inv_matrix)
3235 q_b_matrix(:, :) = matmul(q_b_matrix, s_inv_matrix)
3236
3237 IF (ASSOCIATED(vhxc_block) .AND. &
3238 dft_control%mtlr_dft_with_perturbation .AND. &
3239 dft_control%mtlr_ikind == ikind .AND. &
3240 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
3241 vhxc_matrix(:, :) = matmul(vhxc_matrix, s_inv_matrix)
3242 DO i = 1, n
3243 dft_control%vhxc(ispin) = dft_control%vhxc(ispin) + vhxc_matrix(i, i)
3244 IF (ispin == 1) dft_control%trq(1) = dft_control%trq(1) + q_a_matrix(i, i)
3245 IF (ispin == 2) dft_control%trq(2) = dft_control%trq(2) + q_b_matrix(i, i)
3246 END DO
3247 END IF
3248
3249 ! Perform the requested manipulations of the (initial) orbital occupations
3250
3251 IF (ASSOCIATED(orbitals)) THEN
3252 IF ((qs_env%scf_env%iter_delta >= eps_scf) .OR. &
3253 ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
3254 (qs_env%scf_env%iter_count <= max_scf))) THEN
3255
3256 ALLOCATE (orb_occ(nsbsize))
3257 ALLOCATE (q_eigval(n))
3258 ALLOCATE (q_eigvec(n, n))
3259
3260 q_eigval(:) = 0.0_dp
3261 q_eigvec(:, :) = 0.0_dp
3262 norb = SIZE(orbitals)
3263
3264 IF (ispin == 1) THEN
3265 CALL jacobi(q_a_matrix, q_eigval, q_eigvec)
3266 q_a_matrix(:, :) = 0.0_dp
3267 ELSE
3268 CALL jacobi(q_b_matrix, q_eigval, q_eigvec)
3269 q_b_matrix(:, :) = 0.0_dp
3270 END IF
3271
3272 DO isb = 1, nsb
3273 trq = 0.0_dp
3274 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
3275 trq = trq + q_eigval(i)
3276 END DO
3277
3278 IF (smear) THEN
3279 occ = trq/real(norb, kind=dp)
3280 ELSE
3281 occ = 1.0_dp/fspin
3282 END IF
3283
3284 orb_occ(:) = .false.
3285 iloc = maxloc(q_eigvec(:, isb*nsbsize))
3286 jsb = int((iloc(1) - 1)/nsbsize) + 1
3287 i = 0
3288 i0 = (jsb - 1)*nsbsize + 1
3289 iorb = -1000
3290
3291 DO j = i0, jsb*nsbsize
3292 i = i + 1
3293 IF (i > norb) THEN
3294 DO m = -lu, lu
3295 IF (.NOT. orb_occ(lu + m + 1)) THEN
3296 iorb = i0 + lu + m
3297 orb_occ(lu + m + 1) = .true.
3298 END IF
3299 END DO
3300 ELSE
3301 iorb = i0 + lu + orbitals(i)
3302 orb_occ(lu + orbitals(i) + 1) = .true.
3303 END IF
3304
3305 cpassert(iorb /= -1000)
3306
3307 iloc = maxloc(q_eigvec(iorb, :))
3308 q_eigval(iloc(1)) = min(occ, trq)
3309
3310 IF (ispin == 1) THEN
3311 q_a_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1))
3312 ELSE
3313 q_b_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1))
3314 END IF
3315
3316 trq = trq - q_eigval(iloc(1))
3317 END DO
3318 END DO
3319
3320 IF (ispin == 1) THEN
3321 q_a_matrix(:, :) = matmul(q_a_matrix, transpose(q_eigvec))
3322 ELSE
3323 q_b_matrix(:, :) = matmul(q_b_matrix, transpose(q_eigvec))
3324 END IF
3325
3326 DEALLOCATE (orb_occ)
3327 DEALLOCATE (q_eigval)
3328 DEALLOCATE (q_eigvec)
3329 END IF
3330 END IF
3331
3332 trq = 0.0_dp
3333 trq2 = 0.0_dp
3334 trqxq_hund_j = 0.0_dp
3335
3336 IF (ispin == 1) THEN
3337 DO i = 1, n
3338 trq = trq + q_a_matrix(i, i)
3339 DO j = 1, n
3340 trq2 = trq2 + q_a_matrix(i, j)*q_a_matrix(j, i)
3341 trqxq_hund_j = trqxq_hund_j + q_a_matrix(i, j)*q_b_matrix(j, i)
3342
3343 IF (i == j) THEN
3344 v_matrix(i, i) = u_minus_j*(0.5_dp - fspin*q_a_matrix(i, i))
3345 v_matrix(i, i) = v_matrix(i, i) + hund_j*fspin*q_b_matrix(i, i)
3346 IF (dft_control%mtlr_dft_with_perturbation .AND. &
3347 dft_control%mtlr_ikind == ikind .AND. &
3348 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
3349 v_matrix(i, i) = v_matrix(i, i) + perturbation_strength
3350 END IF
3351 ELSE
3352 v_matrix(i, j) = -u_minus_j*fspin*q_a_matrix(j, i)
3353 v_matrix(i, j) = v_matrix(i, j) + hund_j*fspin*q_b_matrix(j, i)
3354 END IF
3355 END DO
3356 END DO
3357 ELSE
3358 DO i = 1, n
3359 trq = trq + q_b_matrix(i, i)
3360 DO j = 1, n
3361 trq2 = trq2 + q_b_matrix(i, j)*q_b_matrix(j, i)
3362 trqxq_hund_j = trqxq_hund_j + q_b_matrix(i, j)*q_a_matrix(j, i)
3363
3364 IF (i == j) THEN
3365 v_matrix(i, i) = u_minus_j*(0.5_dp - fspin*q_b_matrix(i, i))
3366 v_matrix(i, i) = v_matrix(i, i) + hund_j*fspin*q_a_matrix(i, i)
3367 IF (dft_control%mtlr_dft_with_perturbation .AND. &
3368 dft_control%mtlr_ikind == ikind .AND. &
3369 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
3370 v_matrix(i, i) = v_matrix(i, i) + perturbation_strength
3371 END IF
3372 ELSE
3373 v_matrix(i, j) = -u_minus_j*fspin*q_b_matrix(j, i)
3374 v_matrix(i, j) = v_matrix(i, j) + hund_j*fspin*q_a_matrix(j, i)
3375 END IF
3376 END DO
3377 END DO
3378 END IF
3379
3380 v_matrix(:, :) = matmul(s_inv_matrix, v_matrix)
3381
3382 trq = fspin*trq
3383 trq2 = fspin*fspin*trq2
3384 trqxq_hund_j = fspin*fspin*trqxq_hund_j
3385
3386 energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
3387 energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*hund_j*trqxq_hund_j/fspin
3388
3389 ! Calculate potential V(U) = dE(U)/dq
3390
3391 IF (.NOT. just_energy) THEN
3392
3393 CALL dbcsr_get_block_p(matrix=sm_v, &
3394 row=atom_a, &
3395 col=atom_a, &
3396 block=v_block, &
3397 found=found)
3398 cpassert(ASSOCIATED(v_block))
3399
3400 CALL dbcsr_get_block_p(matrix=sm_a, &
3401 row=atom_a, &
3402 col=atom_a, &
3403 block=a_block, &
3404 found=found)
3405 cpassert(ASSOCIATED(a_block))
3406
3407 isb = 0
3408 DO iset = 1, nset
3409 DO ishell = 1, nshell(iset)
3410 IF (l(ishell, iset) /= lu) cycle
3411 isb = isb + 1
3412 i = 0
3413 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3414 i = i + 1
3415 jsb = 0
3416 DO jset = 1, nset
3417 DO jshell = 1, nshell(jset)
3418 IF (l(jshell, jset) /= lu) cycle
3419 jsb = jsb + 1
3420 j = 0
3421 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
3422 j = j + 1
3423 v_block(isgf, jsgf) = v_matrix(i, j)*a(isb)*a(jsb)
3424 a_block(isgf, jsgf) = s_inv_matrix(i, j)*a(isb)*a(jsb)
3425 END DO
3426 END DO
3427 END DO
3428 END DO
3429 END DO
3430 END DO
3431
3432 END IF
3433
3434 END IF
3435
3436 ! Consider print requests
3437
3438 IF (should_output .AND. (print_level > low_print_level)) THEN
3439 IF (ispin == 1) THEN
3440 CALL para_env%sum(q_a_matrix)
3441 ELSE
3442 CALL para_env%sum(q_b_matrix)
3443 END IF
3444
3445 IF (output_unit > 0) THEN
3446 ALLOCATE (q_work(nsb, nsbsize))
3447 q_work(:, :) = 0.0_dp
3448
3449 DO isb = 1, nsb
3450 j = 0
3451 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
3452 j = j + 1
3453 IF (ispin == 1) THEN
3454 q_work(isb, j) = q_a_matrix(i, i)
3455 ELSE
3456 q_work(isb, j) = q_b_matrix(i, i)
3457 END IF
3458 END DO
3459 END DO
3460
3461 DO isb = 1, nsb
3462 WRITE (unit=output_unit, fmt="(T3,I6,2X,I6,2X,10F8.3)") &
3463 atom_a, isb, q_work(isb, :), sum(q_work(isb, :))
3464 END DO
3465
3466 WRITE (unit=output_unit, fmt="(T12,A,2X,10F8.3)") &
3467 "Total", (sum(q_work(:, i)), i=1, nsbsize), sum(q_work)
3468 WRITE (unit=output_unit, fmt="(A)") ""
3469
3470 DEALLOCATE (q_work)
3471
3472 IF (debug) THEN
3473 WRITE (unit=output_unit, fmt="(T9,70I10)") (i, i=1, n)
3474 DO i = 1, n
3475 IF (ispin == 1) THEN
3476 WRITE (unit=output_unit, fmt="(T3,I6,70F10.6)") i, q_a_matrix(i, :)
3477 ELSE
3478 WRITE (unit=output_unit, fmt="(T3,I6,70F10.6)") i, q_b_matrix(i, :)
3479 END IF
3480 END DO
3481
3482 ALLOCATE (q_eigval(n))
3483 ALLOCATE (q_eigvec(n, n))
3484 q_eigval(:) = 0.0_dp
3485 q_eigvec(:, :) = 0.0_dp
3486
3487 IF (ispin == 1) THEN
3488 CALL jacobi(q_a_matrix, q_eigval, q_eigvec)
3489 ELSE
3490 CALL jacobi(q_b_matrix, q_eigval, q_eigvec)
3491 END IF
3492
3493 WRITE (unit=output_unit, fmt="(/,T9,70I10)") (i, i=1, n)
3494 WRITE (unit=output_unit, fmt="(T9,71F10.6)") (q_eigval(i), i=1, n), sum(q_eigval(1:n))
3495
3496 DO i = 1, n
3497 WRITE (unit=output_unit, fmt="(T3,I6,70F10.6)") i, q_eigvec(i, :)
3498 END DO
3499
3500 DEALLOCATE (q_eigval)
3501 DEALLOCATE (q_eigvec)
3502 END IF
3503 END IF
3504
3505 IF (debug) THEN
3506 ALLOCATE (q_work(nsgf_kind, nsgf_kind))
3507 q_work(:, :) = 0.0_dp
3508
3509 IF (ispin == 1) THEN
3510 IF (ASSOCIATED(q_a_block)) q_work(:, :) = q_a_block(:, :)
3511 ELSE
3512 IF (ASSOCIATED(q_b_block)) q_work(:, :) = q_b_block(:, :)
3513 END IF
3514
3515 CALL para_env%sum(q_work)
3516
3517 IF (output_unit > 0) THEN
3518 norb = SIZE(q_work, 1)
3519
3520 WRITE (unit=output_unit, fmt="(/,T9,200I10)") (i, i=1, norb)
3521 DO i = 1, norb
3522 WRITE (unit=output_unit, fmt="(T3,I6,200F10.6)") i, q_work(i, :)
3523 END DO
3524
3525 ALLOCATE (q_eigval(norb))
3526 ALLOCATE (q_eigvec(norb, norb))
3527 q_eigval(:) = 0.0_dp
3528 q_eigvec(:, :) = 0.0_dp
3529
3530 CALL jacobi(q_work, q_eigval, q_eigvec)
3531
3532 WRITE (unit=output_unit, fmt="(/,T9,200I10)") (i, i=1, norb)
3533 WRITE (unit=output_unit, fmt="(T9,201F10.6)") (q_eigval(i), i=1, norb), sum(q_eigval(1:norb))
3534
3535 DO i = 1, norb
3536 WRITE (unit=output_unit, fmt="(T3,I6,200F10.6)") i, q_eigvec(i, :)
3537 END DO
3538
3539 DEALLOCATE (q_eigval)
3540 DEALLOCATE (q_eigvec)
3541 END IF
3542
3543 DEALLOCATE (q_work)
3544 END IF
3545 END IF
3546
3547 END DO
3548
3549 IF (ALLOCATED(q_a_matrix)) DEALLOCATE (q_a_matrix)
3550 IF (ALLOCATED(q_b_matrix)) DEALLOCATE (q_b_matrix)
3551 IF (ALLOCATED(s_matrix)) DEALLOCATE (s_matrix)
3552 IF (ALLOCATED(s_inv_matrix)) DEALLOCATE (s_inv_matrix)
3553 IF (ALLOCATED(v_matrix)) DEALLOCATE (v_matrix)
3554
3555 DEALLOCATE (a)
3556
3557 IF (dft_control%mtlr_dft_with_perturbation .AND. &
3558 dft_control%mtlr_ikind == ikind) THEN
3559 dft_control%vhxc(ispin) = dft_control%vhxc(ispin)/real(n, dp)
3560 CALL para_env%sum(dft_control%vhxc(ispin))
3561 CALL para_env%sum(dft_control%trq(ispin))
3562 DEALLOCATE (vhxc_matrix)
3563 END IF
3564
3565 END DO
3566
3567 ! Add V(i,j)[U] to V(i,j)[DFT]
3568
3569 IF (ASSOCIATED(sm_h)) THEN
3570 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_v, 0.0_dp, sm_mid)
3571 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 1.0_dp, sm_h)
3572 END IF
3573
3574 END DO
3575
3576 CALL dbcsr_deallocate_matrix(sm_mid)
3577
3578 IF (PRESENT(matrix_w)) THEN
3579
3580 CALL get_qs_env(qs_env=qs_env, force=force)
3581 DO ikind = 1, nkind
3582 force(ikind)%tensorial_u(:, :) = 0.0_dp
3583 END DO
3584 IF (nspin == 2) THEN
3585 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_b)
3586 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, matrix_p(2)%matrix, 0.0_dp, sm_c)
3587 ELSE
3588 CALL dbcsr_multiply("N", "N", 0.5_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_b)
3589 CALL dbcsr_multiply("N", "N", 0.5_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_c)
3590 END IF
3591 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_b, sm_s, 0.0_dp, sm_i)
3592 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_c, sm_s, 0.0_dp, sm_j)
3593
3594 DO jj = 1, 3
3595 CALL dbcsr_desymmetrize(matrixkp_s(jj + 1, 1)%matrix, sm_k)
3596
3597 CALL dbcsr_set(sm_f, 0.0_dp)
3598
3599 CALL dbcsr_iterator_start(iter, sm_k)
3600 DO WHILE (dbcsr_iterator_blocks_left(iter))
3601 NULLIFY (tmp_block)
3602 CALL dbcsr_iterator_next_block(iter, jatom, katom, tmp_block)
3603 IF (jatom == katom) cycle
3604 CALL dbcsr_put_block(matrix=sm_f, row=jatom, col=katom, block=tmp_block)
3605 END DO
3606 CALL dbcsr_iterator_stop(iter)
3607 CALL dbcsr_finalize(sm_f)
3608
3609 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_b, sm_f, 0.0_dp, sm_d)
3610 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_c, sm_f, 0.0_dp, sm_e)
3611
3612 ALLOCATE (a_matrix(nkind, nkind))
3613 ALLOCATE (d_matrix(nkind, nkind))
3614 ALLOCATE (e_matrix(nkind, nkind))
3615 ALLOCATE (i_matrix(nkind, nkind))
3616 ALLOCATE (j_matrix(nkind, nkind))
3617 ALLOCATE (dd_matrix(nkind, nkind))
3618 ALLOCATE (ee_matrix(nkind, nkind))
3619 ALLOCATE (ff_matrix(nkind, nkind))
3620
3621 DO ikind = 1, nkind
3622
3623 CALL get_atomic_kind(atomic_kind_set(ikind), &
3624 atom_list=atom_list, &
3625 natom=natom_of_kind)
3626
3627 CALL get_qs_kind(qs_kind_set(ikind), &
3628 u_minus_j=u_minus_j, &
3629 l_of_dft_plus_u=lu, &
3630 basis_set=orb_basis_set, &
3631 hund_j=hund_j, &
3632 nsgf=nsgf_kind, &
3633 dft_plus_u_atom=dft_plus_u_atom)
3634
3635 IF (.NOT. dft_plus_u_atom) cycle
3636 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
3637 IF (lu < 0) cycle
3638 IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) cycle
3639
3640 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
3641 first_sgf=first_sgf, &
3642 l=l, &
3643 last_sgf=last_sgf, &
3644 nset=nset, &
3645 nshell=nshell)
3646
3647 nsb = 0
3648 DO iset = 1, nset
3649 DO ishell = 1, nshell(iset)
3650 IF (l(ishell, iset) == lu) nsb = nsb + 1
3651 END DO
3652 END DO
3653
3654 nsbsize = (2*lu + 1)
3655 n = nsb*nsbsize
3656
3657 ALLOCATE (a_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3658 ALLOCATE (d_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3659 ALLOCATE (e_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3660 ALLOCATE (i_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3661 ALLOCATE (j_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3662
3663 a_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3664 d_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3665 e_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3666 i_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3667 j_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3668
3669 DO iatom = 1, natom_of_kind
3670
3671 atom_aa = atom_list(iatom)
3672
3673 NULLIFY (a_block)
3674 CALL dbcsr_get_block_p(matrix=sm_a, &
3675 row=atom_aa, &
3676 col=atom_aa, &
3677 block=a_block, &
3678 found=found)
3679
3680 NULLIFY (d_block)
3681 CALL dbcsr_get_block_p(matrix=sm_d, &
3682 row=atom_aa, &
3683 col=atom_aa, &
3684 block=d_block, &
3685 found=found)
3686
3687 NULLIFY (e_block)
3688 CALL dbcsr_get_block_p(matrix=sm_e, &
3689 row=atom_aa, &
3690 col=atom_aa, &
3691 block=e_block, &
3692 found=found)
3693
3694 NULLIFY (i_block)
3695 CALL dbcsr_get_block_p(matrix=sm_i, &
3696 row=atom_aa, &
3697 col=atom_aa, &
3698 block=i_block, &
3699 found=found)
3700
3701 NULLIFY (j_block)
3702 CALL dbcsr_get_block_p(matrix=sm_j, &
3703 row=atom_aa, &
3704 col=atom_aa, &
3705 block=j_block, &
3706 found=found)
3707
3708 IF (ASSOCIATED(a_block) .AND. ASSOCIATED(d_block) .AND. ASSOCIATED(e_block) &
3709 .AND. ASSOCIATED(i_block) .AND. ASSOCIATED(j_block)) THEN
3710
3711 i = 0
3712 DO iset = 1, nset
3713 DO ishell = 1, nshell(iset)
3714 IF (l(ishell, iset) /= lu) cycle
3715 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3716 i = i + 1
3717 j = 0
3718 DO jset = 1, nset
3719 DO jshell = 1, nshell(jset)
3720 IF (l(jshell, jset) /= lu) cycle
3721 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
3722 j = j + 1
3723 a_matrix(ikind, ikind)%matrix(i, j, iatom) = a_block(isgf, jsgf)
3724 d_matrix(ikind, ikind)%matrix(i, j, iatom) = -d_block(isgf, jsgf)
3725 e_matrix(ikind, ikind)%matrix(i, j, iatom) = -e_block(isgf, jsgf)
3726 i_matrix(ikind, ikind)%matrix(i, j, iatom) = i_block(isgf, jsgf)
3727 j_matrix(ikind, ikind)%matrix(i, j, iatom) = j_block(isgf, jsgf)
3728 END DO
3729 END DO
3730 END DO
3731 END DO
3732 END DO
3733 END DO
3734
3735 END IF
3736
3737 DO jkind = 1, nkind
3738
3739 CALL get_atomic_kind(atomic_kind_set(jkind), &
3740 atom_list=atom_list2, &
3741 natom=matom_of_kind)
3742
3743 CALL get_qs_kind(qs_kind_set(jkind), nsgf=nsgf_kind2)
3744
3745 IF (.NOT. ALLOCATED(dd_matrix(ikind, jkind)%matrix)) THEN
3746 ALLOCATE (dd_matrix(ikind, jkind)%matrix(n, nsgf_kind2, natom_of_kind, matom_of_kind))
3747 dd_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3748 END IF
3749
3750 IF (.NOT. ALLOCATED(ee_matrix(ikind, jkind)%matrix)) THEN
3751 ALLOCATE (ee_matrix(ikind, jkind)%matrix(n, nsgf_kind2, natom_of_kind, matom_of_kind))
3752 ee_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3753 END IF
3754
3755 IF (.NOT. ALLOCATED(ff_matrix(ikind, jkind)%matrix)) THEN
3756 ALLOCATE (ff_matrix(ikind, jkind)%matrix(nsgf_kind2, n, natom_of_kind, matom_of_kind))
3757 ff_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3758 END IF
3759
3760 DO jatom = 1, matom_of_kind
3761
3762 atom_a = atom_list2(jatom)
3763
3764 NULLIFY (d_block)
3765 CALL dbcsr_get_block_p(matrix=sm_b, &
3766 row=atom_aa, &
3767 col=atom_a, &
3768 block=d_block, &
3769 found=found)
3770
3771 NULLIFY (e_block)
3772 CALL dbcsr_get_block_p(matrix=sm_c, &
3773 row=atom_aa, &
3774 col=atom_a, &
3775 block=e_block, &
3776 found=found)
3777
3778 IF (ASSOCIATED(d_block) .AND. ASSOCIATED(e_block)) THEN
3779 i = 0
3780 DO iset = 1, nset
3781 DO ishell = 1, nshell(iset)
3782 IF (l(ishell, iset) /= lu) cycle
3783 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3784 i = i + 1
3785 dd_matrix(ikind, jkind)%matrix(i, :, iatom, jatom) = d_block(isgf, :)
3786 ee_matrix(ikind, jkind)%matrix(i, :, iatom, jatom) = e_block(isgf, :)
3787 END DO
3788 END DO
3789 END DO
3790 END IF
3791
3792 NULLIFY (f_block)
3793 CALL dbcsr_get_block_p(matrix=sm_f, &
3794 row=atom_a, &
3795 col=atom_aa, &
3796 block=f_block, &
3797 found=found)
3798
3799 IF (ASSOCIATED(f_block)) THEN
3800 i = 0
3801 DO iset = 1, nset
3802 DO ishell = 1, nshell(iset)
3803 IF (l(ishell, iset) /= lu) cycle
3804 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3805 i = i + 1
3806 ff_matrix(ikind, jkind)%matrix(:, i, iatom, jatom) = f_block(:, isgf)
3807 END DO
3808 END DO
3809 END DO
3810 END IF
3811
3812 END DO
3813 END DO
3814 END DO
3815 END DO
3816
3817 DO ikind = 1, nkind
3818
3819 CALL get_atomic_kind(atomic_kind_set(ikind), &
3820 atom_list=atom_list, &
3821 natom=natom_of_kind)
3822
3823 CALL get_qs_kind(qs_kind_set(ikind), &
3824 u_minus_j=u_minus_j, &
3825 l_of_dft_plus_u=lu, &
3826 basis_set=orb_basis_set, &
3827 hund_j=hund_j, &
3828 nsgf=nsgf_kind, &
3829 dft_plus_u_atom=dft_plus_u_atom)
3830
3831 IF (.NOT. dft_plus_u_atom) cycle
3832 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
3833 IF (lu < 0) cycle
3834 IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) cycle
3835
3836 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
3837 first_sgf=first_sgf, &
3838 l=l, &
3839 last_sgf=last_sgf, &
3840 nset=nset, &
3841 nshell=nshell)
3842
3843 CALL para_env%sum(a_matrix(ikind, ikind)%matrix(:, :, :))
3844 CALL para_env%sum(d_matrix(ikind, ikind)%matrix(:, :, :))
3845 CALL para_env%sum(e_matrix(ikind, ikind)%matrix(:, :, :))
3846 CALL para_env%sum(i_matrix(ikind, ikind)%matrix(:, :, :))
3847 CALL para_env%sum(j_matrix(ikind, ikind)%matrix(:, :, :))
3848
3849 DO jkind = 1, nkind
3850 CALL para_env%sum(dd_matrix(ikind, jkind)%matrix(:, :, :, :))
3851 CALL para_env%sum(ee_matrix(ikind, jkind)%matrix(:, :, :, :))
3852 CALL para_env%sum(ff_matrix(ikind, jkind)%matrix(:, :, :, :))
3853 END DO
3854
3855 END DO
3856
3857 ! Force calculation
3858
3859 DO ikind = 1, nkind
3860
3861 CALL get_atomic_kind(atomic_kind_set(ikind), &
3862 atom_list=atom_list, &
3863 natom=natom_of_kind)
3864
3865 CALL get_qs_kind(qs_kind_set(ikind), &
3866 u_minus_j=u_minus_j, &
3867 l_of_dft_plus_u=lu, &
3868 basis_set=orb_basis_set, &
3869 hund_j=hund_j, &
3870 nsgf=nsgf_kind, &
3871 dft_plus_u_atom=dft_plus_u_atom)
3872
3873 IF (.NOT. dft_plus_u_atom) cycle
3874 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
3875 IF (lu < 0) cycle
3876 IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) cycle
3877
3878 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
3879 first_sgf=first_sgf, &
3880 l=l, &
3881 last_sgf=last_sgf, &
3882 nset=nset, &
3883 nshell=nshell)
3884
3885 nsb = 0
3886 DO iset = 1, nset
3887 DO ishell = 1, nshell(iset)
3888 IF (l(ishell, iset) == lu) nsb = nsb + 1
3889 END DO
3890 END DO
3891
3892 nsbsize = (2*lu + 1)
3893 n = nsb*nsbsize
3894
3895 ALLOCATE (tmp_matrix1(n, n))
3896 ALLOCATE (tmp_matrix2(n, n))
3897 ALLOCATE (tmp_matrix3(n, n))
3898
3899 DO iatom = 1, natom_of_kind
3900
3901 atom_aa = atom_list(iatom)
3902
3903 tmp_matrix1(:, :) = matmul(d_matrix(ikind, ikind)%matrix(:, :, iatom), &
3904 a_matrix(ikind, ikind)%matrix(:, :, iatom))
3905 tmp_matrix2(:, :) = matmul(e_matrix(ikind, ikind)%matrix(:, :, iatom), &
3906 a_matrix(ikind, ikind)%matrix(:, :, iatom))
3907
3908 i_matrix(ikind, ikind)%matrix(:, :, iatom) = matmul(i_matrix(ikind, ikind)%matrix(:, :, iatom), &
3909 a_matrix(ikind, ikind)%matrix(:, :, iatom))
3910 j_matrix(ikind, ikind)%matrix(:, :, iatom) = matmul(j_matrix(ikind, ikind)%matrix(:, :, iatom), &
3911 a_matrix(ikind, ikind)%matrix(:, :, iatom))
3912
3913 tmp_matrix3(:, :) = 0.5_dp*u_minus_j*tmp_matrix1 - &
3914 u_minus_j*matmul(tmp_matrix1, i_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3915 + hund_j*matmul(tmp_matrix2, i_matrix(ikind, ikind)%matrix(:, :, iatom))
3916 tmp_matrix3(:, :) = tmp_matrix3 + 0.5_dp*u_minus_j*tmp_matrix2 - &
3917 u_minus_j*matmul(tmp_matrix2, j_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3918 + hund_j*matmul(tmp_matrix1, j_matrix(ikind, ikind)%matrix(:, :, iatom))
3919
3920 DO i = 1, n
3921 force(ikind)%tensorial_u(jj, iatom) = force(ikind)%tensorial_u(jj, iatom) + tmp_matrix3(i, i)
3922 END DO
3923
3924 DO jkind = 1, nkind
3925
3926 CALL get_atomic_kind(atomic_kind_set(jkind), &
3927 atom_list=atom_list2, &
3928 natom=matom_of_kind)
3929
3930 CALL get_qs_kind(qs_kind_set(jkind), nsgf=nsgf_kind2)
3931
3932 DO jatom = 1, matom_of_kind
3933
3934 atom_a = atom_list2(jatom)
3935
3936 tmp_matrix1(:, :) = matmul(dd_matrix(ikind, jkind)%matrix(:, :, iatom, jatom), &
3937 ff_matrix(ikind, jkind)%matrix(:, :, iatom, jatom))
3938 tmp_matrix1(:, :) = matmul(tmp_matrix1, a_matrix(ikind, ikind)%matrix(:, :, iatom))
3939
3940 tmp_matrix2(:, :) = matmul(ee_matrix(ikind, jkind)%matrix(:, :, iatom, jatom), &
3941 ff_matrix(ikind, jkind)%matrix(:, :, iatom, jatom))
3942 tmp_matrix2(:, :) = matmul(tmp_matrix2, a_matrix(ikind, ikind)%matrix(:, :, iatom))
3943
3944 tmp_matrix3(:, :) = 0.5_dp*u_minus_j*tmp_matrix1 - &
3945 u_minus_j*matmul(tmp_matrix1, i_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3946 + hund_j*matmul(tmp_matrix2, i_matrix(ikind, ikind)%matrix(:, :, iatom))
3947 tmp_matrix3(:, :) = tmp_matrix3 + 0.5_dp*u_minus_j*tmp_matrix2 - &
3948 u_minus_j*matmul(tmp_matrix2, j_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3949 + hund_j*matmul(tmp_matrix1, j_matrix(ikind, ikind)%matrix(:, :, iatom))
3950
3951 IF (atom_a /= atom_aa) THEN
3952 DO i = 1, n
3953 force(jkind)%tensorial_u(jj, jatom) = force(jkind)%tensorial_u(jj, jatom) &
3954 + tmp_matrix3(i, i)
3955 END DO
3956 END IF
3957
3958 END DO
3959 END DO
3960 END DO
3961
3962 DEALLOCATE (tmp_matrix1)
3963 DEALLOCATE (tmp_matrix2)
3964 DEALLOCATE (tmp_matrix3)
3965
3966 END DO
3967
3968 DO ikind = 1, nkind
3969
3970 IF (ALLOCATED(a_matrix(ikind, ikind)%matrix)) DEALLOCATE (a_matrix(ikind, ikind)%matrix)
3971 IF (ALLOCATED(d_matrix(ikind, ikind)%matrix)) DEALLOCATE (d_matrix(ikind, ikind)%matrix)
3972 IF (ALLOCATED(e_matrix(ikind, ikind)%matrix)) DEALLOCATE (e_matrix(ikind, ikind)%matrix)
3973 IF (ALLOCATED(i_matrix(ikind, ikind)%matrix)) DEALLOCATE (i_matrix(ikind, ikind)%matrix)
3974 IF (ALLOCATED(j_matrix(ikind, ikind)%matrix)) DEALLOCATE (j_matrix(ikind, ikind)%matrix)
3975
3976 DO jkind = 1, nkind
3977 IF (ALLOCATED(dd_matrix(ikind, jkind)%matrix)) DEALLOCATE (dd_matrix(ikind, jkind)%matrix)
3978 IF (ALLOCATED(ee_matrix(ikind, jkind)%matrix)) DEALLOCATE (ee_matrix(ikind, jkind)%matrix)
3979 IF (ALLOCATED(ff_matrix(ikind, jkind)%matrix)) DEALLOCATE (ff_matrix(ikind, jkind)%matrix)
3980 END DO
3981
3982 END DO
3983
3984 DEALLOCATE (a_matrix)
3985 DEALLOCATE (d_matrix)
3986 DEALLOCATE (e_matrix)
3987 DEALLOCATE (i_matrix)
3988 DEALLOCATE (j_matrix)
3989 DEALLOCATE (dd_matrix)
3990 DEALLOCATE (ee_matrix)
3991 DEALLOCATE (ff_matrix)
3992
3993 END DO
3994
3995 DO ikind = 1, nkind
3996 force(ikind)%tensorial_u(:, :) = force(ikind)%tensorial_u(:, :)*2.0_dp
3997 END DO
3998
3999 CALL date_and_time(values=timevalues)
4000
4001 END IF
4002
4003 ! Collect the energy contributions from all processes
4004
4005 CALL para_env%sum(energy%dft_plus_u)
4006
4007 IF (energy%dft_plus_u < 0.0_dp) THEN
4008 CALL cp_warn(__location__, &
4009 "DFT+U energy contibution is negative possibly due "// &
4010 "to unphysical Lowdin charges. Check your input, "// &
4011 "if this warning persists or try a different method!")
4012 END IF
4013
4014 ! Release local sparse matrices
4015
4016 CALL dbcsr_deallocate_matrix(sm_q_a)
4017 CALL dbcsr_deallocate_matrix(sm_q_b)
4018 CALL dbcsr_deallocate_matrix(sm_v)
4019 CALL dbcsr_deallocate_matrix(sm_a)
4020
4021 IF (PRESENT(matrix_w)) THEN
4022 CALL dbcsr_deallocate_matrix(sm_b)
4023 CALL dbcsr_deallocate_matrix(sm_c)
4024 CALL dbcsr_deallocate_matrix(sm_d)
4025 CALL dbcsr_deallocate_matrix(sm_e)
4026 CALL dbcsr_deallocate_matrix(sm_f)
4027 CALL dbcsr_deallocate_matrix(sm_k)
4028 CALL dbcsr_deallocate_matrix(sm_i)
4029 CALL dbcsr_deallocate_matrix(sm_j)
4030 END IF
4031
4032 CALL timestop(handle)
4033
4034 END SUBROUTINE tensorial
4035
4036END MODULE dft_plus_u
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_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
mathlib::invert_matrix
Definition mathlib.F:74
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:216
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:131
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::dudarev1997
integer, save, public dudarev1997
Definition bibliography.F:36
bibliography::dudarev1998
integer, save, public dudarev1998
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_transpose
subroutine, public cp_cfm_transpose(matrix, trans, matrixt)
Transposes a BLACS distributed complex matrix.
Definition cp_cfm_basic_linalg.F:1117
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling, first_col, source)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:635
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition cp_cfm_diag.F:94
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_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:632
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_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_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
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_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
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_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1125
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_set_diag
subroutine, public dbcsr_set_diag(matrix, diag)
Copies the diagonal elements from the given array into the given matrix.
Definition cp_dbcsr_contrib.F:634
cp_dbcsr_contrib::dbcsr_get_diag
subroutine, public dbcsr_get_diag(matrix, diag)
Copies the diagonal elements from the given matrix into the given array.
Definition cp_dbcsr_contrib.F:593
cp_dbcsr_contrib::dbcsr_get_block_diag
subroutine, public dbcsr_get_block_diag(matrix, diag)
Copies the diagonal blocks of matrix into diag.
Definition cp_dbcsr_contrib.F:496
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:940
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::write_fm_with_basis_info
subroutine, public write_fm_with_basis_info(blacs_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, output_unit, omit_headers)
Print a spherical matrix of blacs type.
Definition cp_dbcsr_output.F:79
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_transpose
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
Definition cp_fm_basic_linalg.F:1533
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_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::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:262
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
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_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:758
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_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::low_print_level
integer, parameter, public low_print_level
Definition cp_output_handling.F:77
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
dft_plus_u
Add the DFT+U contribution to the Hamiltonian matrix.
Definition dft_plus_u.F:18
dft_plus_u::lowdin_kp_population
subroutine, public lowdin_kp_population(kpoints, delta, alpha)
Cache the current occupations, optionally measure their residual and accept or mix them.
Definition dft_plus_u.F:1134
dft_plus_u::lowdin_kp_initial_density
subroutine, public lowdin_kp_initial_density(qs_env, fmwork)
Complete the initial density's Lowdin coordinates, including atomic and density restarts.
Definition dft_plus_u.F:1074
dft_plus_u::lowdin_kp_kernel
subroutine, public lowdin_kp_kernel(qs_env, density, response)
Project a Bloch density response and return its orbital potential response, including +J.
Definition dft_plus_u.F:1340
dft_plus_u::lowdin_kp_prepare
subroutine, public lowdin_kp_prepare(kp, orbitals, use_real_wfn, smat)
Prepare selected Lowdin projectors and their reusable product workspace.
Definition dft_plus_u.F:1015
dft_plus_u::lowdin_kp_apply
subroutine, public lowdin_kp_apply(kp, v, rmat, cmat)
Add B diag(v) B**H before any solver, DIIS residual, or history storage uses H(k).
Definition dft_plus_u.F:1267
dft_plus_u::lowdin_kp_apply_mo
subroutine, public lowdin_kp_apply_mo(kp, v, coeff, hc)
Add B diag(v) B**H C to the complex orbital gradient, without forming an AO matrix.
Definition dft_plus_u.F:1311
dft_plus_u::plus_u
subroutine, public plus_u(qs_env, matrix_h, matrix_w)
Add the DFT+U contribution to the Hamiltonian matrix. Wrapper routine for all "+U" methods.
Definition dft_plus_u.F:134
dft_plus_u::lowdin_kp_init
subroutine, public lowdin_kp_init(qs_env)
Select the AO columns used by the existing orbital-diagonal Lowdin functional.
Definition dft_plus_u.F:935
dft_plus_u::lowdin_kp_forces
subroutine, public lowdin_kp_forces(kpoints)
Add -dE_U/dS(k) to W(k), before its ordinary inverse Fourier transform. If A=S**(1/2),...
Definition dft_plus_u.F:1380
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::plus_u_lowdin
integer, parameter, public plus_u_lowdin
Definition input_constants.F:945
input_constants::plus_u_mulliken_charges
integer, parameter, public plus_u_mulliken_charges
Definition input_constants.F:945
input_constants::plus_u_mulliken
integer, parameter, public plus_u_mulliken
Definition input_constants.F:945
input_constants::plus_u_tensorial
integer, parameter, public plus_u_tensorial
Definition input_constants.F:945
input_constants::diag_update_method_adiis
integer, parameter, public diag_update_method_adiis
Definition input_constants.F:970
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
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_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
local_gemm_api
Definition local_gemm_api.F:8
local_gemm_api::local_gemm_pu_gpu
integer, parameter, public local_gemm_pu_gpu
Definition local_gemm_api.F:41
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::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::jacobi
subroutine, public jacobi(a, d, v)
Jacobi matrix diagonalization. The eigenvalues are returned in vector d and the eigenvectors are retu...
Definition mathlib.F:1630
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
mulliken
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii}...
Definition mulliken.F:13
orbital_symbols
orbital_symbols
Definition orbital_symbols.F:14
orbital_symbols::sgf_symbol
character(len=6) function, public sgf_symbol(n, l, m)
Build a spherical orbital symbol (orbital labels for printing).
Definition orbital_symbols.F:107
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:123
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
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_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:483
qs_kind_types::set_qs_kind
subroutine, public set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, covalent_radius, vdw_radius, lmax_rho0, zeff, no_optimize, dispersion, u_minus_j, hund_j, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
Definition qs_kind_types.F:3320
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:976
qs_kpoint_operators
Assembly of real and complex k-point operators from real-space DBCSR matrices. Outputs use the distri...
Definition qs_kpoint_operators.F:26
qs_kpoint_operators::kpoint_operator_context_create
subroutine, public kpoint_operator_context_create(ctx, kpoints, ao_ao_fm, nspin)
Create an empty operator context. The context owns no buffers after this call; they are created on fi...
Definition qs_kpoint_operators.F:390
qs_kpoint_operators::kpoint_operator_get
subroutine, public kpoint_operator_get(ctx, ikp, ispin, rsmat, cmat, fm_re, fm_im, mat_re, mat_im, matrix_row, slot, grid)
Assemble one operator for one global round and spin in one call: run the single-group transfer of the...
Definition qs_kpoint_operators.F:1868
qs_kpoint_operators::kpoint_operator_context_release
subroutine, public kpoint_operator_context_release(ctx)
Release the context together with every buffer it created. Call once, at the end of the fixed-layout ...
Definition qs_kpoint_operators.F:427
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_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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
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::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
local_gemm_api::local_gemm_ctxt_type
Definition local_gemm_api.F:47
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:774
particle_types::particle_type
Definition particle_types.F:35
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_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_kpoint_operators::kpoint_operator_context_type
Caller-owned assembly context for a fixed geometry and parallel layout. OT retains it for one inner S...
Definition qs_kpoint_operators.F:167
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