Line data Source code
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 : ! **************************************************************************************************
18 : MODULE 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 :
122 : CONTAINS
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 6386 : 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 6386 : CALL timeset(routineN, handle)
148 :
149 6386 : CPASSERT(ASSOCIATED(qs_env))
150 :
151 6386 : NULLIFY (input, dft_control)
152 :
153 6386 : logger => cp_get_default_logger()
154 :
155 : CALL get_qs_env(qs_env=qs_env, &
156 : input=input, &
157 6386 : dft_control=dft_control)
158 :
159 6386 : CALL cite_reference(Dudarev1997)
160 6386 : 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 6386 : orthonormal_basis = .FALSE.
166 :
167 : ! Setup print control
168 :
169 6386 : 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 6386 : (.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 6386 : log_filename=.FALSE.)
177 :
178 : ! Select DFT+U method
179 :
180 6386 : 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 2556 : should_output, output_unit, print_level)
190 : END IF
191 : CASE (plus_u_mulliken)
192 : CALL mulliken(qs_env, orthonormal_basis, matrix_h, &
193 2096 : 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 962 : should_output, output_unit, print_level)
197 : CASE (plus_u_tensorial)
198 : CALL tensorial(qs_env, matrix_h, matrix_w, &
199 772 : should_output, output_unit, print_level)
200 : CASE DEFAULT
201 6386 : 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 6386 : ignore_should_output=should_output)
206 :
207 6386 : CALL timestop(handle)
208 :
209 6386 : 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 2556 : 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 2556 : 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 2556 : INTEGER, ALLOCATABLE, DIMENSION(:) :: shell_sgf
258 : INTEGER, DIMENSION(1) :: iloc
259 2556 : INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals
260 2556 : 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 2556 : 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 2556 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigval, q_eigval
267 2556 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orbq, q_eigvec, q_matrix, q_opp_matrix, &
268 2556 : q_work, slam
269 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
270 2556 : POINTER :: local_data
271 2556 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: q_block, q_opp_block, v_block, vhxc_block
272 2556 : 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 2556 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vhxc
277 2556 : 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 2556 : 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 2556 : CALL timeset(routineN, handle)
290 :
291 2556 : debug = .FALSE. ! Set to .TRUE. to print debug information
292 :
293 2556 : NULLIFY (sm_h, sm_p, sm_p_opp, sm_s, sm_w, vhxc_block, sm_vhxc, matrix_vhxc)
294 :
295 2556 : smear = .FALSE.
296 2556 : max_scf = -1
297 2556 : eps_scf = 1.0E30_dp
298 2556 : 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 2556 : para_env=para_env)
312 :
313 2556 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
314 :
315 2556 : energy%dft_plus_u = 0.0_dp
316 :
317 2556 : nspin = dft_control%nspins
318 2556 : IF (dft_control%mtlr_dft_with_perturbation) THEN
319 1448 : IF (.NOT. ASSOCIATED(matrix_vhxc)) THEN
320 0 : CPABORT("MTLR requires the projected Hxc matrix for every spin channel.")
321 : END IF
322 1448 : IF (SIZE(matrix_vhxc) /= nspin) THEN
323 0 : CPABORT("The number of projected Hxc matrices does not match the spin channels.")
324 : END IF
325 : END IF
326 :
327 2556 : IF (nspin == 2) THEN
328 : fspin = 1.0_dp
329 : ELSE
330 488 : fspin = 0.5_dp
331 : END IF
332 :
333 : ! Get the number of basis functions and atomic kinds.
334 2556 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
335 :
336 2556 : nkind = SIZE(atomic_kind_set)
337 :
338 2556 : IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
339 : just_energy = .FALSE.
340 : ELSE
341 34 : 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 6472 : DO ikind = 1, nkind
347 3916 : NULLIFY (orbitals)
348 : CALL get_qs_kind(qs_kind_set(ikind), &
349 : hund_j=hund_j, &
350 3916 : orbitals=orbitals)
351 6472 : IF (hund_j /= 0.0_dp) THEN
352 1774 : if_plus_j = .TRUE.
353 1774 : IF (ASSOCIATED(orbitals)) THEN
354 : CALL cp_abort(__LOCATION__, "ENFORCE_OCCUPATION is incompatible with a nonzero Hund J "// &
355 0 : "in the Lowdin DFT+U+J method.")
356 : END IF
357 : END IF
358 : END DO
359 :
360 2556 : NULLIFY (orbitals)
361 :
362 2556 : IF (.NOT. do_kpoints) THEN
363 : ! Retrieve S^(1/2) from the SCF environment
364 1924 : fm_s_half => scf_env%s_half
365 1924 : CPASSERT(ASSOCIATED(fm_s_half))
366 : ! work matrices
367 1924 : CALL cp_fm_get_info(fm_s_half, matrix_struct=fmstruct)
368 : CALL cp_fm_create(matrix=fm_work1, matrix_struct=fmstruct, &
369 1924 : name="FULL WORK MATRIX 1")
370 : CALL cp_fm_create(matrix=fm_work2, matrix_struct=fmstruct, &
371 1924 : 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 2556 : IF (PRESENT(matrix_w) .AND. .NOT. do_kpoints) THEN
378 2 : CALL cp_fm_create(matrix=fm_sev, matrix_struct=fmstruct)
379 2 : CALL cp_fm_create(matrix=slambda, matrix_struct=fmstruct)
380 8 : ALLOCATE (eigval(nsgf), slam(nsgf, 1))
381 2 : sm_s => matrix_s(1, 1)%matrix
382 2 : CALL copy_dbcsr_to_fm(sm_s, fm_work1)
383 2 : CALL choose_eigv_solver(fm_work1, fm_sev, eigval)
384 : !
385 48 : DO i = 1, nsgf
386 48 : IF (eigval(i) > 0._dp) THEN
387 46 : slam(i, 1) = SQRT(eigval(i))
388 : ELSE
389 0 : CPABORT("S matrix not positive definit")
390 : END IF
391 : END DO
392 48 : DO i = 1, nsgf
393 48 : CALL cp_fm_set_submatrix(slambda, slam, 1, i, nsgf, 1, 1.0_dp, 0.0_dp)
394 : END DO
395 48 : DO i = 1, nsgf
396 48 : CALL cp_fm_set_submatrix(slambda, slam, i, 1, 1, nsgf, 1.0_dp, 1.0_dp, .TRUE.)
397 : END DO
398 2 : CALL cp_fm_get_info(slambda, local_data=local_data)
399 48 : DO i = 1, SIZE(local_data, 2)
400 577 : DO j = 1, SIZE(local_data, 1)
401 529 : sij = local_data(j, i)
402 529 : IF (sij > 0.0_dp) sij = 1.0_dp/sij
403 575 : local_data(j, i) = sij
404 : END DO
405 : END DO
406 4 : DEALLOCATE (eigval, slam)
407 : END IF
408 :
409 : ! Calculate S^(1/2)*P*S^(1/2)
410 2556 : 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 632 : qs_env%scf_control%diagonalization%update_method == diag_update_method_adiis .OR. PRESENT(matrix_w)) THEN
415 262 : CALL lowdin_kp_population(kpoints, alpha=1.0_dp)
416 : END IF
417 2528 : ALLOCATE (orbq(nsgf, nspin))
418 632 : orbq = 0.0_dp
419 10430 : orbq(kpoints%lowdin_orbitals, :) = kpoints%lowdin_q
420 : END IF
421 :
422 : ! Create local block diagonal matrices
423 2556 : sm_s => matrix_s(1, 1)%matrix
424 2556 : CALL dbcsr_get_block_diag(sm_s, sm_q)
425 2556 : IF (if_plus_j .AND. nspin == 2) THEN
426 1774 : CALL dbcsr_get_block_diag(sm_s, sm_q_opp)
427 : END IF
428 2556 : CALL dbcsr_get_block_diag(sm_s, sm_v)
429 :
430 : ! Loop over all spins
431 7180 : DO ispin = 1, nspin
432 :
433 4624 : CALL dbcsr_set(sm_q, 0.0_dp)
434 4624 : IF (if_plus_j .AND. nspin == 2) THEN
435 3548 : CALL dbcsr_set(sm_q_opp, 0.0_dp)
436 : END IF
437 4624 : CALL dbcsr_set(sm_v, 0.0_dp)
438 :
439 4624 : IF (dft_control%mtlr_dft_with_perturbation) THEN
440 2896 : IF (ispin == 1) perturbation_strength = 0.9_dp*dft_control%perturbation_strength
441 2896 : IF (ispin == 2) perturbation_strength = 1.1_dp*dft_control%perturbation_strength
442 2896 : IF (.NOT. ASSOCIATED(matrix_vhxc(ispin)%matrix)) THEN
443 0 : CPABORT("MTLR projected Hxc matrix is not initialized.")
444 : END IF
445 2896 : sm_vhxc => matrix_vhxc(ispin)%matrix
446 : END IF
447 :
448 4624 : IF (do_kpoints) THEN
449 858 : CALL dbcsr_set_diag(sm_q, orbq(:, ispin))
450 858 : 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 3766 : sm_p => matrix_p(ispin, 1)%matrix
454 3766 : 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 3766 : 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 3766 : 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 3766 : 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 3548 : sm_p_opp => matrix_p(3 - ispin, 1)%matrix
484 3548 : 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 3548 : matrix_c=fm_work2)
495 3548 : 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 11480 : 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 6856 : 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 6856 : smear=smear)
529 :
530 : ! Check, if the atoms of this atomic kind need a DFT+U correction
531 6856 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
532 6856 : IF (.NOT. dft_plus_u_atom) CYCLE
533 4624 : IF (lu < 0) CYCLE
534 :
535 : ! Apply U ramping if requested
536 4624 : IF ((ispin == 1) .AND. (u_ramping > 0.0_dp) .AND. .NOT. do_kpoints) THEN
537 0 : IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
538 0 : u_minus_j = MIN(u_minus_j + u_ramping, u_minus_j_target)
539 0 : CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
540 : END IF
541 0 : IF (should_output .AND. (output_unit > 0)) THEN
542 : WRITE (UNIT=output_unit, FMT="(T3,A,3X,A,F0.3,A)") &
543 0 : "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)), &
544 0 : "U(eff) = ", u_minus_j*evolt, " eV"
545 : END IF
546 : END IF
547 :
548 4624 : 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 4624 : nshell=nshell)
557 :
558 : ! Select this kind's shell once; occupations and potentials use the same AO order.
559 13872 : ALLOCATE (shell_sgf(nsgf_kind))
560 4624 : n = 0
561 10622 : DO iset = 1, nset
562 26566 : DO ishell = 1, nshell(iset)
563 15944 : IF (l(ishell, iset) /= lu) CYCLE
564 33422 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
565 18176 : n = n + 1
566 34120 : shell_sgf(n) = isgf
567 : END DO
568 : END DO
569 : END DO
570 4624 : nsbsize = 2*lu + 1
571 4624 : nsb = n/nsbsize
572 :
573 18496 : ALLOCATE (q_matrix(n, n))
574 4624 : q_matrix(:, :) = 0.0_dp
575 :
576 4624 : IF (hund_j /= 0.0_dp .AND. nspin == 2) THEN
577 10644 : ALLOCATE (q_opp_matrix(n, n))
578 3548 : q_opp_matrix(:, :) = 0.0_dp
579 : END IF
580 :
581 4624 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
582 : dft_control%mtlr_ikind == ikind) THEN
583 2896 : dft_control%vhxc(ispin) = 0.0_dp
584 2896 : dft_control%trq(ispin) = 0.0_dp
585 : END IF
586 :
587 : ! Print headline if requested
588 4624 : IF (should_output .AND. (print_level > low_print_level)) THEN
589 0 : IF (output_unit > 0) THEN
590 0 : ALLOCATE (symbol(nsbsize))
591 0 : DO m = -lu, lu
592 0 : symbol(lu + m + 1) = sgf_symbol(0, lu, m)
593 : END DO
594 0 : IF (nspin > 1) THEN
595 0 : WRITE (UNIT=spin_info, FMT="(A8,I2)") " of spin", ispin
596 : ELSE
597 0 : spin_info = ""
598 : END IF
599 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
600 0 : "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ", ikind, &
601 0 : ": "//TRIM(atomic_kind_name), &
602 0 : "Atom Shell ", (ADJUSTR(symbol(i)), i=1, nsbsize), " Trace"
603 0 : DEALLOCATE (symbol)
604 : END IF
605 : END IF
606 :
607 : ! Loop over all atoms of the current atomic kind
608 11640 : DO iatom = 1, natom_of_kind
609 7016 : atom_a = atom_list(iatom)
610 7016 : q_matrix(:, :) = 0.0_dp
611 7016 : IF (ALLOCATED(q_opp_matrix)) THEN
612 5940 : 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 7016 : found=found)
621 7016 : 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 5940 : found=found)
627 : END IF
628 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
629 7016 : 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 2896 : found=found)
636 : END IF
637 :
638 7016 : IF (ASSOCIATED(q_block)) THEN
639 : ! Calculate energy contribution to E(U)
640 14988 : DO i = 1, n
641 11480 : isgf = shell_sgf(i)
642 11480 : q_matrix(i, i) = q_block(isgf, isgf)
643 11480 : IF (hund_j /= 0.0_dp .AND. nspin == 2 .AND. ASSOCIATED(q_opp_block)) THEN
644 8252 : 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 14988 : dft_control%mtlr_ikind == ikind .AND. qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
648 4792 : dft_control%vhxc(ispin) = dft_control%vhxc(ispin) + vhxc_block(isgf, isgf)
649 4792 : 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 3508 : IF (ASSOCIATED(orbitals)) THEN
655 68 : 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 66 : ALLOCATE (orb_occ(nsbsize))
659 66 : ALLOCATE (q_eigval(n))
660 22 : q_eigval(:) = 0.0_dp
661 66 : ALLOCATE (q_eigvec(n, n))
662 22 : q_eigvec(:, :) = 0.0_dp
663 22 : norb = SIZE(orbitals)
664 22 : CALL jacobi(q_matrix, q_eigval, q_eigvec)
665 22 : q_matrix(:, :) = 0.0_dp
666 66 : DO isb = 1, nsb
667 44 : trq = 0.0_dp
668 176 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
669 176 : trq = trq + q_eigval(i)
670 : END DO
671 44 : IF (smear) THEN
672 44 : occ = trq/REAL(norb, KIND=dp)
673 : ELSE
674 0 : occ = 1.0_dp/fspin
675 : END IF
676 44 : orb_occ(:) = .FALSE.
677 352 : iloc = MAXLOC(q_eigvec(:, isb*nsbsize))
678 44 : jsb = INT((iloc(1) - 1)/nsbsize) + 1
679 44 : i = 0
680 44 : i0 = (jsb - 1)*nsbsize + 1
681 44 : iorb = -1000
682 198 : DO j = i0, jsb*nsbsize
683 132 : i = i + 1
684 132 : IF (i > norb) THEN
685 0 : DO m = -lu, lu
686 0 : IF (.NOT. orb_occ(lu + m + 1)) THEN
687 0 : iorb = i0 + lu + m
688 0 : orb_occ(lu + m + 1) = .TRUE.
689 : END IF
690 : END DO
691 : ELSE
692 132 : iorb = i0 + lu + orbitals(i)
693 132 : orb_occ(lu + orbitals(i) + 1) = .TRUE.
694 : END IF
695 132 : CPASSERT(iorb /= -1000)
696 1056 : iloc = MAXLOC(q_eigvec(iorb, :))
697 132 : q_eigval(iloc(1)) = MIN(occ, trq)
698 924 : q_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1)) ! backtransform left
699 176 : trq = trq - q_eigval(iloc(1))
700 : END DO
701 : END DO
702 30426 : q_matrix(:, :) = MATMUL(q_matrix, TRANSPOSE(q_eigvec)) ! backtransform right
703 22 : DEALLOCATE (orb_occ)
704 22 : DEALLOCATE (q_eigval)
705 22 : DEALLOCATE (q_eigvec)
706 : END IF
707 : END IF ! orbitals associated
708 :
709 26468 : trq = fspin*SUM([(q_matrix(i, i), i=1, n)])
710 64732 : trq2 = fspin*fspin*SUM(q_matrix*TRANSPOSE(q_matrix))
711 3508 : trqxq_hund_j = 0.0_dp
712 3508 : IF (hund_j /= 0.0_dp) THEN
713 2970 : IF (nspin == 2) THEN
714 41598 : trqxq_hund_j = SUM(q_matrix*TRANSPOSE(q_opp_matrix))
715 : ELSE
716 0 : trqxq_hund_j = SUM(q_matrix*TRANSPOSE(q_matrix))
717 : END IF
718 : END IF
719 :
720 : ! Calculate energy contribution to E(U)
721 3508 : energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
722 3508 : IF (hund_j /= 0.0_dp) THEN
723 2970 : trqxq_hund_j = fspin*fspin*trqxq_hund_j
724 2970 : 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 3508 : 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 3474 : found=found)
734 3474 : CPASSERT(ASSOCIATED(v_block))
735 :
736 14750 : DO j = 1, n
737 11276 : jsgf = shell_sgf(j)
738 63270 : DO i = 1, n
739 48520 : isgf = shell_sgf(i)
740 48520 : IF (i == j) THEN
741 11276 : 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 11276 : dft_control%mtlr_ikind == ikind .AND. &
744 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
745 4792 : v_block(isgf, isgf) = v_block(isgf, isgf) + perturbation_strength
746 : END IF
747 : ELSE
748 37244 : CPASSERT(ABS(q_matrix(j, i)) < 1.0E-14_dp)
749 37244 : IF (hund_j /= 0.0_dp .AND. nspin == 2) THEN
750 21704 : CPASSERT(ABS(q_opp_matrix(j, i)) < 1.0E-14_dp)
751 : END IF
752 37244 : v_block(isgf, jsgf) = -u_minus_j*fspin*q_matrix(j, i)
753 : END IF
754 59796 : IF (hund_j /= 0.0_dp) THEN
755 29872 : IF (nspin == 2) THEN
756 29872 : v_block(isgf, jsgf) = v_block(isgf, jsgf) + hund_j*fspin*q_opp_matrix(j, i)
757 : ELSE
758 0 : 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 18656 : IF (should_output .AND. (print_level > low_print_level)) THEN
769 0 : CALL para_env%sum(q_matrix)
770 0 : IF (output_unit > 0) THEN
771 0 : ALLOCATE (q_work(nsb, nsbsize))
772 0 : q_work(:, :) = 0.0_dp
773 0 : DO isb = 1, nsb
774 0 : j = 0
775 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
776 0 : j = j + 1
777 0 : q_work(isb, j) = q_matrix(i, i)
778 : END DO
779 : END DO
780 0 : DO isb = 1, nsb
781 : WRITE (UNIT=output_unit, FMT="(T3,I6,2X,I6,2X,10F8.3)") &
782 0 : 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 0 : "Total", (SUM(q_work(:, i)), i=1, nsbsize), SUM(q_work)
786 0 : WRITE (UNIT=output_unit, FMT="(A)") ""
787 0 : 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 4624 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
843 : dft_control%mtlr_ikind == ikind) THEN
844 2896 : dft_control%vhxc(ispin) = dft_control%vhxc(ispin)/REAL(n, dp)
845 2896 : CALL para_env%sum(dft_control%vhxc(ispin))
846 2896 : CALL para_env%sum(dft_control%trq(ispin))
847 : END IF
848 :
849 4624 : DEALLOCATE (shell_sgf, q_matrix)
850 :
851 20728 : IF (ALLOCATED(q_opp_matrix)) THEN
852 3548 : 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 4624 : IF (PRESENT(matrix_h) .OR. (PRESENT(matrix_w) .AND. do_kpoints)) THEN
859 4554 : IF (do_kpoints) THEN
860 : BLOCK
861 1716 : REAL(dp) :: diagonal(nsgf)
862 858 : CALL dbcsr_get_diag(sm_v, diagonal)
863 858 : CALL para_env%sum(diagonal)
864 11154 : kpoints%lowdin_v(:, ispin) = diagonal(kpoints%lowdin_orbitals)
865 : END BLOCK
866 : ELSE
867 3696 : sm_h => matrix_h(ispin, 1)%matrix
868 3696 : CALL cp_dbcsr_sm_fm_multiply(sm_v, fm_s_half, fm_work1, nsgf)
869 3696 : CALL cp_fm_transpose(fm_work1, fm_work2)
870 3696 : 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 7112 : IF (PRESENT(matrix_w) .AND. .NOT. do_kpoints) THEN
877 :
878 2 : sm_p => matrix_p(ispin, 1)%matrix
879 2 : sm_w => matrix_w(ispin, 1)%matrix
880 :
881 2 : CALL cp_dbcsr_sm_fm_multiply(sm_v, fm_s_half, fm_work1, nsgf)
882 2 : CALL cp_fm_transpose(fm_work1, fm_work2)
883 2 : CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_work2, fm_work1, nsgf)
884 2 : CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, 1.0_dp, fm_work1, fm_sev, 0.0_dp, fm_work2)
885 2 : CALL parallel_gemm('T', 'N', nsgf, nsgf, nsgf, 1.0_dp, fm_sev, fm_work2, 0.0_dp, fm_work1)
886 2 : CALL cp_fm_schur_product(fm_work1, slambda, fm_work2)
887 2 : CALL cp_fm_transpose(fm_work2, fm_work1)
888 : ! Symmetrize before projecting onto the stored triangular sparse blocks.
889 2 : CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_work2, matrix_b=fm_work1)
890 2 : CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, -1.0_dp, fm_sev, fm_work2, 0.0_dp, fm_work1)
891 2 : 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 2556 : IF (PRESENT(matrix_w)) THEN
898 10 : IF (do_kpoints) CALL lowdin_kp_forces(kpoints)
899 10 : CALL cp_fm_release(matrix=fm_sev)
900 10 : CALL cp_fm_release(matrix=slambda)
901 : END IF
902 :
903 : ! Collect the energy contributions from all processes
904 :
905 2556 : CALL para_env%sum(energy%dft_plus_u)
906 :
907 2556 : IF (energy%dft_plus_u < 0.0_dp) THEN
908 : CALL cp_warn(__LOCATION__, &
909 : "DFT+U energy contribution is negative possibly due "// &
910 0 : "to unphysical Lowdin charges!")
911 : END IF
912 :
913 : ! Release (local) full matrices
914 2556 : NULLIFY (fm_s_half)
915 2556 : CALL cp_fm_release(matrix=fm_work1)
916 2556 : CALL cp_fm_release(matrix=fm_work2)
917 :
918 : ! Release (local) sparse matrices
919 2556 : CALL dbcsr_release(sm_q)
920 2556 : CALL dbcsr_release(sm_v)
921 :
922 2556 : IF (if_plus_j .AND. nspin == 2) THEN
923 1774 : CALL dbcsr_release(sm_q_opp)
924 : END IF
925 :
926 2556 : CALL timestop(handle)
927 :
928 5112 : 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 120 : 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 120 : INTEGER, POINTER :: atoms(:), first(:, :), l(:, :), &
941 120 : last(:, :), nshell(:), orbitals(:)
942 : LOGICAL :: active
943 : LOGICAL, ALLOCATABLE :: selected(:)
944 : REAL(KIND=dp) :: j, u
945 120 : 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 120 : TYPE(particle_type), POINTER :: particles(:)
951 120 : 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 120 : particle_set=particles, dft_control=dft_control)
955 120 : IF (dft_control%mtlr_dft_with_perturbation) THEN
956 0 : CPABORT("MTLR with Lowdin DFT+U and k-points is not implemented")
957 : END IF
958 600 : DO ik = 1, SIZE(kpoints%kp_sym)
959 600 : IF (kpoints%kp_sym(ik)%kpoint_sym%apply_symmetry) THEN
960 0 : CPABORT("Lowdin DFT+U requires a full k-point grid or time-reversal-only reduction")
961 : END IF
962 : END DO
963 120 : CALL get_qs_kind_set(qs_kinds, nsgf=nao)
964 120 : nspin = dft_control%nspins
965 600 : ALLOCATE (first_atom(SIZE(particles)), selected(nao))
966 120 : CALL get_particle_set(particles, qs_kinds, first_sgf=first_atom)
967 120 : selected = .FALSE.
968 360 : DO ikind = 1, SIZE(qs_kinds)
969 240 : NULLIFY (orbitals)
970 : CALL get_qs_kind(qs_kinds(ikind), dft_plus_u_atom=active, l_of_dft_plus_u=lu, &
971 240 : basis_set=basis, u_minus_j_target=u, hund_j=j, orbitals=orbitals)
972 240 : IF (.NOT. active) CYCLE
973 120 : IF (ASSOCIATED(orbitals)) THEN
974 0 : CPABORT("Lowdin DFT+U with k-points does not support ENFORCE_OCCUPATION")
975 : END IF
976 120 : IF (u == 0.0_dp .AND. j == 0.0_dp) CYCLE
977 120 : CALL get_atomic_kind(kinds_set(ikind), atom_list=atoms)
978 120 : CALL get_gto_basis_set(basis, first_sgf=first, last_sgf=last, l=l, nset=nset, nshell=nshell)
979 600 : DO ia = 1, SIZE(atoms)
980 480 : DO iset = 1, nset
981 840 : DO ish = 1, nshell(iset)
982 600 : IF (l(ish, iset) /= lu) CYCLE
983 : selected(first_atom(atoms(ia)) + first(ish, iset) - 1: &
984 1440 : first_atom(atoms(ia)) + last(ish, iset) - 1) = .TRUE.
985 : END DO
986 : END DO
987 : END DO
988 : END DO
989 120 : IF (ALLOCATED(kpoints%lowdin_orbitals)) THEN
990 112 : DEALLOCATE (kpoints%lowdin_orbitals, kpoints%lowdin_q, kpoints%lowdin_q_new, kpoints%lowdin_v)
991 : END IF
992 2880 : m = COUNT(selected)
993 360 : ALLOCATE (kpoints%lowdin_orbitals(m))
994 6480 : kpoints%lowdin_orbitals(:) = PACK([(i, i=1, nao)], selected)
995 960 : ALLOCATE (kpoints%lowdin_q(m, nspin), kpoints%lowdin_q_new(m, nspin), kpoints%lowdin_v(m, nspin))
996 1268 : kpoints%lowdin_q(:, :) = 0.0_dp
997 1268 : kpoints%lowdin_q_new(:, :) = 0.0_dp
998 1268 : kpoints%lowdin_v(:, :) = 0.0_dp
999 120 : kpoints%lowdin_density_ready = .FALSE.
1000 120 : kpoints%lowdin_population_ready = .FALSE.
1001 360 : DO ik = 1, SIZE(kpoints%kp_env)
1002 240 : kp => kpoints%kp_env(ik)%kpoint_env
1003 360 : CALL lowdin_kp_prepare(kp, kpoints%lowdin_orbitals, kpoints%use_real_wfn)
1004 : END DO
1005 360 : 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 240 : 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 240 : REAL(dp), ALLOCATABLE :: eigenvalues(:)
1022 : TYPE(cp_cfm_type) :: evec, work
1023 :
1024 240 : CALL cp_fm_release(kp%lowdin_b)
1025 240 : CALL cp_fm_release(kp%lowdin_work)
1026 240 : CALL cp_cfm_release(kp%lowdin_cb)
1027 240 : CALL cp_cfm_release(kp%lowdin_cwork)
1028 240 : CALL kp%lowdin_gemm%destroy()
1029 240 : m = SIZE(orbitals)
1030 240 : IF (m == 0) RETURN
1031 240 : IF (PRESENT(smat)) THEN
1032 0 : CPASSERT(.NOT. use_real_wfn)
1033 0 : CALL cp_cfm_get_info(smat, nrow_global=n)
1034 0 : ALLOCATE (eigenvalues(n))
1035 0 : CALL cp_cfm_create(evec, smat%matrix_struct)
1036 0 : CALL cp_cfm_create(work, smat%matrix_struct)
1037 0 : CALL cp_cfm_release(kp%cshalf)
1038 0 : CALL cp_cfm_create(kp%cshalf, smat%matrix_struct)
1039 0 : CALL cp_cfm_to_cfm(smat, work)
1040 0 : CALL cp_cfm_heevd(work, evec, eigenvalues)
1041 0 : CPASSERT(ALL(eigenvalues >= 0.0_dp))
1042 0 : CALL cp_cfm_to_cfm(evec, work)
1043 0 : CALL cp_cfm_column_scale(work, CMPLX(SQRT(eigenvalues), 0.0_dp, dp))
1044 0 : CALL parallel_gemm('N', 'C', n, n, n, z_one, work, evec, z_zero, kp%cshalf)
1045 0 : CALL cp_cfm_release(evec)
1046 0 : CALL cp_cfm_release(work)
1047 : END IF
1048 240 : CALL kp%lowdin_gemm%create(LOCAL_GEMM_PU_GPU, timing=.FALSE.)
1049 240 : IF (use_real_wfn) THEN
1050 0 : CALL cp_fm_create(kp%lowdin_b, kp%shalf%matrix_struct, ncol=m)
1051 0 : CALL cp_fm_create(kp%lowdin_work, kp%lowdin_b%matrix_struct)
1052 0 : DO i = 1, m
1053 0 : 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 0 : CALL cp_cfm_create(kp%lowdin_cb, kp%lowdin_b%matrix_struct)
1057 0 : CALL cp_fm_to_cfm(kp%lowdin_b, mtarget=kp%lowdin_cb)
1058 : ELSE
1059 240 : CALL cp_cfm_create(kp%lowdin_cb, kp%cshalf%matrix_struct, ncol=m)
1060 1680 : DO i = 1, m
1061 1680 : CALL cp_cfm_to_cfm(kp%cshalf, kp%lowdin_cb, 1, orbitals(i), i)
1062 : END DO
1063 : END IF
1064 240 : CALL cp_cfm_create(kp%lowdin_cwork, kp%lowdin_cb%matrix_struct)
1065 240 : IF (PRESENT(smat)) CALL cp_cfm_release(kp%cshalf)
1066 480 : 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 120 : 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 120 : TYPE(dbcsr_p_type), POINTER :: density(:, :)
1081 : TYPE(kpoint_env_type), POINTER :: kp
1082 120 : TYPE(kpoint_operator_context_type) :: ctx
1083 : TYPE(kpoint_type), POINTER :: kpoints
1084 : TYPE(qs_rho_type), POINTER :: rho
1085 :
1086 120 : CALL get_qs_env(qs_env, kpoints=kpoints, rho=rho)
1087 120 : IF (.NOT. ALLOCATED(kpoints%lowdin_orbitals)) RETURN
1088 120 : 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 8 : CALL qs_rho_get(rho, rho_ao_kp=density)
1092 8 : CALL kpoint_operator_context_create(ctx, kpoints, fmwork, SIZE(density, 1))
1093 8 : kp => kpoints%kp_env(1)%kpoint_env
1094 8 : IF (kpoints%use_real_wfn) THEN
1095 0 : CALL cp_fm_create(pmat, kp%shalf%matrix_struct)
1096 : ELSE
1097 8 : CALL cp_cfm_create(cpmat, kp%cshalf%matrix_struct)
1098 : END IF
1099 8 : kplocal = SIZE(kpoints%kp_env)
1100 24 : kpmax = MAXVAL(kpoints%kp_dist(2, :) - kpoints%kp_dist(1, :) + 1)
1101 92 : kpoints%lowdin_q_new(:, :) = 0.0_dp
1102 24 : DO ik = 1, kpmax
1103 48 : DO ispin = 1, SIZE(density, 1)
1104 24 : IF (kpoints%use_real_wfn) THEN
1105 0 : CALL kpoint_operator_get(ctx, ik, ispin, density, fm_re=pmat, matrix_row=ispin)
1106 : ELSE
1107 24 : CALL kpoint_operator_get(ctx, ik, ispin, density, cmat=cpmat, matrix_row=ispin)
1108 : END IF
1109 24 : IF (ik > kplocal) CYCLE
1110 24 : kp => kpoints%kp_env(ik)%kpoint_env
1111 40 : IF (kpoints%use_real_wfn) THEN
1112 0 : CALL lowdin_kp_accumulate(kp, kpoints%lowdin_q_new(:, ispin), pmat=pmat)
1113 : ELSE
1114 24 : CALL lowdin_kp_accumulate(kp, kpoints%lowdin_q_new(:, ispin), cpmat=cpmat)
1115 : END IF
1116 : END DO
1117 : END DO
1118 8 : CALL kpoints%para_env%sum(kpoints%lowdin_q_new)
1119 8 : CALL kpoint_operator_context_release(ctx)
1120 8 : CALL cp_fm_release(pmat)
1121 8 : CALL cp_cfm_release(cpmat)
1122 8 : kpoints%lowdin_population_ready = .TRUE.
1123 : END IF
1124 120 : CALL lowdin_kp_population(kpoints, alpha=1.0_dp)
1125 240 : 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 265240 : 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 265240 : IF (.NOT. PRESENT(kpoints)) RETURN
1139 265240 : IF (.NOT. ASSOCIATED(kpoints)) RETURN
1140 265240 : IF (.NOT. ALLOCATED(kpoints%lowdin_q)) RETURN
1141 1222 : IF (.NOT. kpoints%lowdin_population_ready) THEN
1142 1968 : IF (SIZE(kpoints%lowdin_q) > 0) THEN
1143 656 : CPASSERT(kpoints%lowdin_density_ready)
1144 656 : CALL lowdin_kp_project(kpoints, kpoints%lowdin_q_new)
1145 : END IF
1146 656 : kpoints%lowdin_population_ready = .TRUE.
1147 : END IF
1148 5954 : IF (PRESENT(delta)) delta = MAX(delta, MAXVAL(ABS(kpoints%lowdin_q_new - kpoints%lowdin_q)))
1149 1222 : IF (PRESENT(alpha)) THEN
1150 1222 : IF (alpha == 1.0_dp) THEN
1151 9664 : kpoints%lowdin_q(:, :) = kpoints%lowdin_q_new
1152 294 : ELSE IF (alpha /= 0.0_dp) THEN
1153 0 : 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 656 : 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 656 : TYPE(cp_cfm_type), ALLOCATABLE :: work(:)
1172 : TYPE(kpoint_env_type), POINTER :: kp
1173 :
1174 656 : m = SIZE(q, 1)
1175 656 : IF (m == 0) RETURN
1176 656 : CALL timeset("lowdin_kp_population", handle)
1177 656 : nspin = SIZE(q, 2)
1178 656 : kplocal = SIZE(kpoints%kp_env)
1179 656 : kp => kpoints%kp_env(1)%kpoint_env
1180 656 : nworkers = 1
1181 656 : !$ IF (kpoints%para_env_kp%num_pe == 1) nworkers = MIN(omp_get_max_threads(), kplocal)
1182 2624 : ALLOCATE (work(nworkers))
1183 6802 : q(:, :) = 0.0_dp
1184 656 : IF (PRESENT(density) .OR. .NOT. kpoints%use_real_wfn) THEN
1185 1312 : DO thread = 1, nworkers
1186 1312 : 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 656 : !$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 656 : CALL kpoints%para_env%sum(q)
1210 1312 : DO thread = 1, nworkers
1211 1312 : CALL cp_cfm_release(work(thread))
1212 : END DO
1213 656 : CALL timestop(handle)
1214 1312 : 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 1780 : 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 1780 : INTEGER, POINTER :: columns(:)
1231 : LOGICAL :: local
1232 :
1233 1780 : m = SIZE(q)
1234 0 : CPASSERT(PRESENT(pmat) .NEQV. PRESENT(cpmat))
1235 1780 : IF (PRESENT(pmat)) THEN
1236 0 : CALL cp_fm_get_info(kp%lowdin_b, nrow_global=n, nrow_local=nr, ncol_local=nc, col_indices=columns)
1237 0 : local = PRODUCT(kp%lowdin_b%matrix_struct%context%num_pe) == 1
1238 0 : 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 0 : kp%lowdin_work%local_data, SIZE(kp%lowdin_work%local_data, 1))
1242 : ELSE
1243 0 : 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 0 : DO j = 1, nc
1246 : q(columns(j)) = q(columns(j)) + kp%wkp* &
1247 0 : SUM(kp%lowdin_b%local_data(:nr, j)*kp%lowdin_work%local_data(:nr, j))
1248 : END DO
1249 : ELSE
1250 1780 : CALL cp_cfm_get_info(kp%lowdin_cb, nrow_global=n, nrow_local=nr, ncol_local=nc, col_indices=columns)
1251 1780 : CALL lowdin_kp_gemm('N', 'N', n, m, n, cpmat, kp%lowdin_cb, z_zero, kp%lowdin_cwork, kp%lowdin_gemm)
1252 12460 : DO j = 1, nc
1253 : q(columns(j)) = q(columns(j)) + kp%wkp*REAL(SUM(CONJG(kp%lowdin_cb%local_data(:nr, j))* &
1254 258100 : kp%lowdin_cwork%local_data(:nr, j)), dp)
1255 : END DO
1256 : END IF
1257 1780 : 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 792 : 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 792 : INTEGER, POINTER :: columns(:)
1274 : LOGICAL :: local
1275 :
1276 792 : m = SIZE(v)
1277 0 : IF (m == 0) RETURN
1278 792 : IF (ALL(v == 0.0_dp)) RETURN
1279 792 : CPASSERT(PRESENT(rmat) .NEQV. PRESENT(cmat))
1280 792 : IF (PRESENT(rmat)) THEN
1281 0 : CALL cp_fm_get_info(kp%lowdin_b, nrow_global=n, nrow_local=nr, ncol_local=nc, col_indices=columns)
1282 0 : DO j = 1, nc
1283 0 : kp%lowdin_work%local_data(:nr, j) = v(columns(j))*kp%lowdin_b%local_data(:nr, j)
1284 : END DO
1285 0 : local = PRODUCT(kp%lowdin_b%matrix_struct%context%num_pe) == 1
1286 0 : 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 0 : rmat%local_data, SIZE(rmat%local_data, 1))
1291 : ELSE
1292 0 : 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 792 : CALL cp_cfm_get_info(kp%lowdin_cb, nrow_global=n, nrow_local=nr, ncol_local=nc, col_indices=columns)
1296 5544 : DO j = 1, nc
1297 114840 : kp%lowdin_cwork%local_data(:nr, j) = v(columns(j))*kp%lowdin_cb%local_data(:nr, j)
1298 : END DO
1299 792 : 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 792 : 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 584 : 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 584 : INTEGER, POINTER :: rows(:)
1318 : TYPE(cp_cfm_type) :: projected
1319 :
1320 584 : m = SIZE(v)
1321 584 : IF (m == 0 .OR. ALL(v == 0.0_dp)) RETURN
1322 584 : CALL cp_cfm_get_info(coeff, nrow_global=n, ncol_global=nmo)
1323 584 : CALL cp_cfm_create(projected, coeff%matrix_struct, nrow=m)
1324 584 : CALL lowdin_kp_gemm('C', 'N', m, nmo, n, kp%lowdin_cb, coeff, z_zero, projected, kp%lowdin_gemm)
1325 584 : CALL cp_cfm_get_info(projected, nrow_local=nr, ncol_local=nc, row_indices=rows)
1326 4088 : DO i = 1, nr
1327 18104 : projected%local_data(i, :nc) = v(rows(i))*projected%local_data(i, :nc)
1328 : END DO
1329 584 : CALL lowdin_kp_gemm('N', 'N', n, nmo, m, kp%lowdin_cb, projected, z_one, hc, kp%lowdin_gemm)
1330 584 : CALL cp_cfm_release(projected)
1331 584 : 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 0 : 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 0 : INTEGER, ALLOCATABLE :: last(:)
1346 : REAL(KIND=dp) :: j, u
1347 : TYPE(kpoint_type), POINTER :: kpoints
1348 0 : TYPE(particle_type), POINTER :: particles(:)
1349 0 : TYPE(qs_kind_type), POINTER :: qs_kinds(:)
1350 :
1351 0 : CALL get_qs_env(qs_env, kpoints=kpoints, particle_set=particles, qs_kind_set=qs_kinds)
1352 0 : m = SIZE(kpoints%lowdin_orbitals)
1353 0 : nspin = SIZE(density, 3)
1354 0 : ALLOCATE (response(m, nspin))
1355 0 : response(:, :) = 0.0_dp
1356 0 : IF (m == 0) RETURN
1357 0 : CALL lowdin_kp_project(kpoints, response, density)
1358 0 : ALLOCATE (last(SIZE(particles)))
1359 0 : CALL get_particle_set(particles, qs_kinds, last_sgf=last)
1360 0 : iatom = 1
1361 0 : DO i = 1, m
1362 0 : DO WHILE (kpoints%lowdin_orbitals(i) > last(iatom))
1363 0 : iatom = iatom + 1
1364 : END DO
1365 0 : CALL get_atomic_kind(particles(iatom)%atomic_kind, kind_number=ikind)
1366 0 : CALL get_qs_kind(qs_kinds(ikind), u_minus_j=u, hund_j=j)
1367 0 : response(i, :) = 0.5_dp*REAL(nspin, dp)*(j*response(i, nspin:1:-1) - u*response(i, :))
1368 : END DO
1369 0 : 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 8 : SUBROUTINE lowdin_kp_forces(kpoints)
1380 : TYPE(kpoint_type), POINTER :: kpoints
1381 :
1382 : INTEGER :: handle, ik, ispin, iw, nworkers, thread
1383 8 : TYPE(cp_cfm_type), ALLOCATABLE :: work(:, :)
1384 : TYPE(kpoint_env_type), POINTER :: kp
1385 8 : TYPE(local_gemm_ctxt_type), ALLOCATABLE :: gemm(:)
1386 :
1387 8 : IF (SIZE(kpoints%lowdin_orbitals) == 0) RETURN
1388 8 : IF (ALL(kpoints%lowdin_v == 0.0_dp)) RETURN
1389 8 : CALL timeset("lowdin_kp_forces", handle)
1390 8 : nworkers = 1
1391 8 : !$ IF (kpoints%para_env_kp%num_pe == 1) &
1392 8 : !$ nworkers = MIN(omp_get_max_threads(), SIZE(kpoints%kp_env)*SIZE(kpoints%lowdin_v, 2))
1393 80 : ALLOCATE (work(3, nworkers), gemm(nworkers))
1394 8 : kp => kpoints%kp_env(1)%kpoint_env
1395 16 : DO thread = 1, nworkers
1396 24 : DO iw = 1, 2
1397 24 : CALL cp_cfm_create(work(iw, thread), kp%pmat(1, 1)%matrix_struct)
1398 : END DO
1399 8 : CALL cp_cfm_create(work(3, thread), kp%lowdin_cb%matrix_struct)
1400 16 : 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 8 : !$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 16 : DO thread = 1, nworkers
1414 8 : CALL gemm(thread)%destroy()
1415 40 : DO iw = 1, 3
1416 32 : CALL cp_cfm_release(work(iw, thread))
1417 : END DO
1418 : END DO
1419 8 : CALL timestop(handle)
1420 16 : 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 24 : 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 24 : INTEGER, POINTER :: columns(:), rows(:)
1441 : LOGICAL :: local, real_only
1442 :
1443 24 : IF (ALL(v == 0.0_dp)) RETURN
1444 24 : CPASSERT(ALLOCATED(kp%lowdin_sqrt))
1445 24 : CPASSERT(ASSOCIATED(kp%lowdin_evec%local_data))
1446 24 : CPASSERT(ASSOCIATED(kp%wmat))
1447 : CALL cp_cfm_get_info(work(1), nrow_global=n, nrow_local=nr, ncol_local=nc, &
1448 24 : row_indices=rows, col_indices=columns)
1449 24 : m = SIZE(orbitals)
1450 72 : local = PRODUCT(work(1)%matrix_struct%context%num_pe) == 1
1451 24 : real_only = SIZE(kp%pmat, 1) == 1
1452 24 : IF (real_only) THEN
1453 0 : CALL cp_fm_to_cfm(kp%pmat(1, ispin), mtarget=work(1))
1454 : ELSE
1455 24 : 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 24 : CALL lowdin_kp_gemm('N', 'N', n, m, n, work(1), kp%lowdin_cb, z_zero, work(2), gemm)
1460 24 : CALL lowdin_kp_gemm('C', 'N', n, m, n, kp%lowdin_evec, work(2), z_zero, work(1), gemm)
1461 576 : DO j = 1, nc
1462 3888 : IF (columns(j) <= m) work(1)%local_data(:nr, j) = v(columns(j))*work(1)%local_data(:nr, j)
1463 : END DO
1464 24 : IF (local) THEN
1465 168 : DO j = 1, m
1466 3480 : work(3)%local_data(:n, j) = CONJG(kp%lowdin_evec%local_data(orbitals(j), :n))
1467 : END DO
1468 : ELSE
1469 0 : CALL cp_cfm_transpose(kp%lowdin_evec, 'C', work(2))
1470 0 : DO j = 1, m
1471 0 : CALL cp_cfm_to_cfm(work(2), work(3), 1, orbitals(j), j)
1472 : END DO
1473 : END IF
1474 24 : 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 576 : DO j = 1, nc
1477 13272 : DO i = 1, nr
1478 : work(2)%local_data(i, j) = work(2)%local_data(i, j)/ &
1479 13248 : (kp%lowdin_sqrt(rows(i)) + kp%lowdin_sqrt(columns(j)))
1480 : END DO
1481 : END DO
1482 24 : CALL lowdin_kp_gemm('N', 'N', n, n, n, kp%lowdin_evec, work(2), z_zero, work(1), gemm)
1483 24 : CALL lowdin_kp_gemm('N', 'C', n, n, n, work(1), kp%lowdin_evec, z_zero, work(2), gemm)
1484 24 : IF (.NOT. local) CALL cp_cfm_transpose(work(2), 'C', work(1))
1485 576 : DO j = 1, nc
1486 13272 : DO i = 1, nr
1487 12696 : IF (local) THEN
1488 12696 : value = work(2)%local_data(i, j) + CONJG(work(2)%local_data(j, i))
1489 : ELSE
1490 0 : value = work(2)%local_data(i, j) + work(1)%local_data(i, j)
1491 : END IF
1492 12696 : kp%wmat(1, ispin)%local_data(i, j) = kp%wmat(1, ispin)%local_data(i, j) - REAL(value, dp)
1493 13248 : IF (.NOT. real_only) THEN
1494 12696 : 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 24 : 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 3860 : 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 11580 : 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 3860 : beta, c%local_data, SIZE(c%local_data, 1))
1525 : ELSE
1526 0 : CALL parallel_gemm(ta, tb, nrow, ncol, inner, z_one, a, b, beta, c)
1527 : END IF
1528 3860 : 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 2096 : 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 2096 : 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 2096 : INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals
1580 2096 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
1581 : LOGICAL :: debug, dft_plus_u_atom, found, &
1582 : just_energy, occupation_enforced, smear
1583 2096 : 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 2096 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: q_eigval
1588 2096 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: q_eigvec, q_matrix, q_work
1589 2096 : REAL(KIND=dp), DIMENSION(:), POINTER :: nelec
1590 2096 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, q_block, s_block, &
1591 2096 : v_block
1592 2096 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1593 : TYPE(atomic_kind_type), POINTER :: kind_a
1594 2096 : 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 2096 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1600 : TYPE(qs_energy_type), POINTER :: energy
1601 2096 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1602 : TYPE(qs_rho_type), POINTER :: rho
1603 :
1604 2096 : CALL timeset(routineN, handle)
1605 :
1606 2096 : debug = .FALSE. ! Set to .TRUE. to print debug information
1607 :
1608 2096 : NULLIFY (atom_list)
1609 2096 : NULLIFY (atomic_kind_set)
1610 2096 : NULLIFY (qs_kind_set)
1611 2096 : NULLIFY (dft_control)
1612 2096 : NULLIFY (energy)
1613 2096 : NULLIFY (first_sgf)
1614 2096 : NULLIFY (h_block)
1615 2096 : NULLIFY (matrix_p)
1616 2096 : NULLIFY (matrix_s)
1617 2096 : NULLIFY (l)
1618 2096 : NULLIFY (last_sgf)
1619 2096 : NULLIFY (nelec)
1620 2096 : NULLIFY (nshell)
1621 2096 : NULLIFY (orb_basis_set)
1622 2096 : NULLIFY (p_block)
1623 2096 : NULLIFY (particle_set)
1624 2096 : NULLIFY (q_block)
1625 2096 : NULLIFY (rho)
1626 2096 : NULLIFY (s_block)
1627 2096 : NULLIFY (orbitals)
1628 2096 : NULLIFY (sm_h)
1629 2096 : NULLIFY (sm_p)
1630 2096 : NULLIFY (sm_q)
1631 2096 : NULLIFY (sm_s)
1632 2096 : NULLIFY (sm_v)
1633 2096 : NULLIFY (v_block)
1634 2096 : NULLIFY (para_env)
1635 :
1636 2096 : smear = .FALSE.
1637 2096 : max_scf = -1
1638 2096 : eps_scf = 1.0E30_dp
1639 2096 : 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 2096 : para_env=para_env)
1649 :
1650 2096 : CPASSERT(ASSOCIATED(atomic_kind_set))
1651 2096 : CPASSERT(ASSOCIATED(dft_control))
1652 2096 : CPASSERT(ASSOCIATED(energy))
1653 2096 : CPASSERT(ASSOCIATED(particle_set))
1654 2096 : CPASSERT(ASSOCIATED(rho))
1655 :
1656 2096 : 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 2096 : matrix_s_kp=matrix_s)
1662 2096 : CPASSERT(ASSOCIATED(matrix_s))
1663 : END IF
1664 2096 : nimg = dft_control%nimages
1665 :
1666 : ! Get density matrices in sparse format
1667 :
1668 2096 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1669 :
1670 2096 : energy%dft_plus_u = 0.0_dp
1671 :
1672 2096 : nspin = dft_control%nspins
1673 :
1674 2096 : IF (nspin == 2) THEN
1675 : fspin = 1.0_dp
1676 : ELSE
1677 1114 : 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 2096 : natom=natom)
1685 :
1686 2096 : nkind = SIZE(atomic_kind_set)
1687 :
1688 6288 : ALLOCATE (is_plus_u_kind(nkind))
1689 2096 : is_plus_u_kind(:) = .FALSE.
1690 :
1691 2096 : IF (PRESENT(matrix_h)) THEN
1692 : just_energy = .FALSE.
1693 : ELSE
1694 578 : just_energy = .TRUE.
1695 : END IF
1696 :
1697 : ! Loop over all spins
1698 5174 : DO ispin = 1, nspin
1699 :
1700 : ! Loop over cell images
1701 101670 : DO ic = 1, nimg
1702 96496 : IF (.NOT. orthonormal_basis) THEN
1703 96496 : sm_s => matrix_s(1, ic)%matrix
1704 : END IF
1705 :
1706 96496 : IF (PRESENT(matrix_h)) THEN
1707 : ! Hamiltonian matrix for spin ispin in sparse format
1708 94616 : 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 96496 : sm_p => matrix_p(ispin, ic)%matrix
1716 :
1717 96496 : IF (.NOT. ASSOCIATED(sm_q)) THEN
1718 2096 : ALLOCATE (sm_q)
1719 2096 : CALL dbcsr_get_block_diag(sm_p, sm_q)
1720 : END IF
1721 96496 : CALL dbcsr_set(sm_q, 0.0_dp)
1722 :
1723 96496 : IF (.NOT. ASSOCIATED(sm_v)) THEN
1724 2096 : ALLOCATE (sm_v)
1725 2096 : CALL dbcsr_get_block_diag(sm_p, sm_v)
1726 : END IF
1727 96496 : CALL dbcsr_set(sm_v, 0.0_dp)
1728 :
1729 385984 : 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 289488 : found=found)
1736 :
1737 289488 : 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 144744 : found=found)
1744 144744 : CPASSERT(ASSOCIATED(q_block))
1745 :
1746 675472 : IF (orthonormal_basis) THEN
1747 : ! S is the unit matrix
1748 0 : DO isgf = 1, SIZE(q_block, 1)
1749 0 : 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 144744 : found=found)
1757 144744 : CPASSERT(ASSOCIATED(s_block))
1758 : ! Exploit that P and S are symmetric
1759 1219678 : DO jsgf = 1, SIZE(p_block, 2)
1760 11325736 : DO isgf = 1, SIZE(p_block, 1)
1761 11036248 : 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 289488 : 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 192992 : 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 192992 : smear=smear)
1800 :
1801 : ! Check, if the atoms of this atomic kind need a DFT+U correction
1802 :
1803 192992 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
1804 192992 : IF (.NOT. dft_plus_u_atom) CYCLE
1805 96496 : IF (lu < 0) CYCLE
1806 :
1807 : ! Apply U ramping if requested
1808 :
1809 96496 : IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
1810 976 : IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
1811 464 : u_minus_j = MIN(u_minus_j + u_ramping, u_minus_j_target)
1812 464 : CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
1813 : END IF
1814 976 : IF (should_output .AND. (output_unit > 0)) THEN
1815 : WRITE (UNIT=output_unit, FMT="(T3,A,3X,A,F0.3,A)") &
1816 476 : "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)), &
1817 952 : "U(eff) = ", u_minus_j*evolt, " eV"
1818 : END IF
1819 : END IF
1820 :
1821 96496 : IF (u_minus_j == 0.0_dp) CYCLE
1822 :
1823 96496 : 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 96496 : nshell=nshell)
1833 :
1834 : ! Count the relevant shell blocks of this atomic kind
1835 :
1836 96496 : nsb = 0
1837 194912 : DO iset = 1, nset
1838 663484 : DO ishell = 1, nshell(iset)
1839 566988 : IF (l(ishell, iset) == lu) nsb = nsb + 1
1840 : END DO
1841 : END DO
1842 :
1843 96496 : nsbsize = (2*lu + 1)
1844 96496 : n = nsb*nsbsize
1845 :
1846 385984 : ALLOCATE (q_matrix(n, n))
1847 96496 : q_matrix(:, :) = 0.0_dp
1848 :
1849 : ! Print headline if requested
1850 :
1851 96496 : IF (should_output .AND. (print_level > low_print_level)) THEN
1852 0 : IF (output_unit > 0) THEN
1853 0 : ALLOCATE (symbol(nsbsize))
1854 0 : DO m = -lu, lu
1855 0 : symbol(lu + m + 1) = sgf_symbol(0, lu, m)
1856 : END DO
1857 0 : IF (nspin > 1) THEN
1858 0 : WRITE (UNIT=spin_info, FMT="(A8,I2)") " of spin", ispin
1859 : ELSE
1860 0 : spin_info = ""
1861 : END IF
1862 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
1863 0 : "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ", ikind, &
1864 0 : ": "//TRIM(atomic_kind_name), &
1865 0 : "Atom Shell ", (ADJUSTR(symbol(i)), i=1, nsbsize), " Trace"
1866 0 : DEALLOCATE (symbol)
1867 : END IF
1868 : END IF
1869 :
1870 : ! Loop over all atoms of the current atomic kind
1871 :
1872 192992 : DO iatom = 1, natom_of_kind
1873 :
1874 96496 : atom_a = atom_list(iatom)
1875 :
1876 96496 : 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 96496 : found=found)
1885 :
1886 : ! Calculate energy contribution to E(U)
1887 :
1888 96496 : IF (ASSOCIATED(q_block)) THEN
1889 :
1890 48248 : i = 0
1891 97456 : DO iset = 1, nset
1892 331742 : DO ishell = 1, nshell(iset)
1893 234286 : IF (l(ishell, iset) /= lu) CYCLE
1894 421284 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1895 282534 : i = i + 1
1896 282534 : j = 0
1897 805114 : DO jset = 1, nset
1898 1948728 : DO jshell = 1, nshell(jset)
1899 1377900 : IF (l(jshell, jset) /= lu) CYCLE
1900 2479026 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
1901 1660434 : j = j + 1
1902 3038334 : 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 48248 : IF (ASSOCIATED(orbitals)) THEN
1913 0 : 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 0 : ALLOCATE (orb_occ(nsbsize))
1917 0 : ALLOCATE (q_eigval(n))
1918 0 : q_eigval(:) = 0.0_dp
1919 0 : ALLOCATE (q_eigvec(n, n))
1920 0 : q_eigvec(:, :) = 0.0_dp
1921 0 : norb = SIZE(orbitals)
1922 0 : CALL jacobi(q_matrix, q_eigval, q_eigvec)
1923 0 : q_matrix(:, :) = 0.0_dp
1924 0 : IF (nelec(ispin) >= 0.5_dp) THEN
1925 0 : trq = nelec(ispin)/SUM(q_eigval(1:n))
1926 0 : q_eigval(1:n) = trq*q_eigval(1:n)
1927 : END IF
1928 0 : DO isb = 1, nsb
1929 0 : trq = 0.0_dp
1930 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
1931 0 : trq = trq + q_eigval(i)
1932 : END DO
1933 0 : IF (smear) THEN
1934 0 : occ = trq/REAL(norb, KIND=dp)
1935 : ELSE
1936 0 : occ = 1.0_dp/fspin
1937 : END IF
1938 0 : orb_occ(:) = .FALSE.
1939 0 : iloc = MAXLOC(q_eigvec(:, isb*nsbsize))
1940 0 : jsb = INT((iloc(1) - 1)/nsbsize) + 1
1941 0 : i = 0
1942 0 : i0 = (jsb - 1)*nsbsize + 1
1943 0 : iorb = -1000
1944 0 : DO j = i0, jsb*nsbsize
1945 0 : i = i + 1
1946 0 : IF (i > norb) THEN
1947 0 : DO m = -lu, lu
1948 0 : IF (.NOT. orb_occ(lu + m + 1)) THEN
1949 0 : iorb = i0 + lu + m
1950 0 : orb_occ(lu + m + 1) = .TRUE.
1951 : END IF
1952 : END DO
1953 : ELSE
1954 0 : iorb = i0 + lu + orbitals(i)
1955 0 : orb_occ(lu + orbitals(i) + 1) = .TRUE.
1956 : END IF
1957 0 : CPASSERT(iorb /= -1000)
1958 0 : iloc = MAXLOC(q_eigvec(iorb, :))
1959 0 : q_eigval(iloc(1)) = MIN(occ, trq)
1960 0 : q_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1)) ! backtransform left
1961 0 : trq = trq - q_eigval(iloc(1))
1962 : END DO
1963 : END DO
1964 0 : q_matrix(:, :) = MATMUL(q_matrix, TRANSPOSE(q_eigvec)) ! backtransform right
1965 0 : DEALLOCATE (orb_occ)
1966 0 : DEALLOCATE (q_eigval)
1967 0 : DEALLOCATE (q_eigvec)
1968 0 : occupation_enforced = .TRUE.
1969 : END IF
1970 : END IF ! orbitals associated
1971 :
1972 48248 : trq = 0.0_dp
1973 48248 : trq2 = 0.0_dp
1974 :
1975 330782 : DO i = 1, n
1976 282534 : trq = trq + q_matrix(i, i)
1977 1991216 : DO j = 1, n
1978 1942968 : trq2 = trq2 + q_matrix(i, j)*q_matrix(j, i)
1979 : END DO
1980 : END DO
1981 :
1982 48248 : trq = fspin*trq
1983 48248 : trq2 = fspin*fspin*trq2
1984 :
1985 : ! Calculate energy contribution to E(U)
1986 :
1987 48248 : 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 48248 : 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 47308 : found=found)
1998 47308 : CPASSERT(ASSOCIATED(v_block))
1999 :
2000 47308 : i = 0
2001 95152 : DO iset = 1, nset
2002 324738 : DO ishell = 1, nshell(iset)
2003 229586 : IF (l(ishell, iset) /= lu) CYCLE
2004 412400 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
2005 276894 : i = i + 1
2006 276894 : j = 0
2007 786590 : DO jset = 1, nset
2008 1906704 : DO jshell = 1, nshell(jset)
2009 1349700 : IF (l(jshell, jset) /= lu) CYCLE
2010 2425722 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
2011 1626594 : j = j + 1
2012 2976294 : IF (isgf == jsgf) THEN
2013 276894 : v_block(isgf, isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i, i))
2014 : ELSE
2015 1349700 : 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 289488 : IF (should_output .AND. (print_level > low_print_level)) THEN
2031 0 : CALL para_env%sum(q_matrix)
2032 0 : IF (output_unit > 0) THEN
2033 0 : ALLOCATE (q_work(nsb, nsbsize))
2034 0 : q_work(:, :) = 0.0_dp
2035 0 : DO isb = 1, nsb
2036 0 : j = 0
2037 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
2038 0 : j = j + 1
2039 0 : q_work(isb, j) = q_matrix(i, i)
2040 : END DO
2041 : END DO
2042 0 : DO isb = 1, nsb
2043 : WRITE (UNIT=output_unit, FMT="(T3,I6,2X,I6,2X,10F8.3)") &
2044 0 : 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 0 : "Total", (SUM(q_work(:, i)), i=1, nsbsize), SUM(q_work)
2048 0 : WRITE (UNIT=output_unit, FMT="(A)") ""
2049 0 : 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 482480 : IF (ALLOCATED(q_matrix)) THEN
2105 96496 : 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 99574 : IF (ASSOCIATED(sm_h)) THEN
2113 :
2114 283848 : DO ikind = 1, nkind
2115 :
2116 189232 : IF (.NOT. is_plus_u_kind(ikind)) CYCLE
2117 :
2118 94616 : kind_a => atomic_kind_set(ikind)
2119 :
2120 : CALL get_atomic_kind(atomic_kind=kind_a, &
2121 : atom_list=atom_list, &
2122 94616 : natom=natom_of_kind)
2123 :
2124 283848 : DO iatom = 1, natom_of_kind
2125 :
2126 94616 : 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 94616 : found=found)
2133 :
2134 94616 : 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 47308 : found=found)
2141 47308 : CPASSERT(ASSOCIATED(v_block))
2142 :
2143 378464 : IF (orthonormal_basis) THEN
2144 0 : DO isgf = 1, SIZE(h_block, 1)
2145 0 : 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 47308 : found=found)
2153 47308 : CPASSERT(ASSOCIATED(s_block))
2154 641450 : DO jsgf = 1, SIZE(h_block, 2)
2155 8203984 : DO isgf = 1, SIZE(h_block, 1)
2156 8109368 : 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 2096 : CALL para_env%sum(energy%dft_plus_u)
2174 :
2175 2096 : IF (energy%dft_plus_u < 0.0_dp) THEN
2176 2 : IF (.NOT. occupation_enforced) THEN
2177 : CALL cp_warn(__LOCATION__, &
2178 : "DFT+U energy contribution is negative possibly due "// &
2179 2 : "to unphysical Mulliken charges!")
2180 : END IF
2181 : END IF
2182 :
2183 2096 : CALL dbcsr_deallocate_matrix(sm_q)
2184 2096 : CALL dbcsr_deallocate_matrix(sm_v)
2185 :
2186 2096 : CALL timestop(handle)
2187 :
2188 6288 : 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 962 : 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 962 : 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 962 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
2244 962 : INTEGER, DIMENSION(:), POINTER :: atom_list, nshell
2245 962 : 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 962 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dEdq, trps
2250 962 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: q_ii
2251 962 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, s_block, w_block
2252 962 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2253 : TYPE(dbcsr_iterator_type) :: iter
2254 962 : 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 962 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2260 : TYPE(qs_energy_type), POINTER :: energy
2261 962 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2262 : TYPE(qs_rho_type), POINTER :: rho
2263 :
2264 962 : CALL timeset(routineN, handle)
2265 :
2266 962 : NULLIFY (atom_list)
2267 962 : NULLIFY (atomic_kind_set)
2268 962 : NULLIFY (qs_kind_set)
2269 962 : NULLIFY (dft_control)
2270 962 : NULLIFY (energy)
2271 962 : NULLIFY (first_sgf)
2272 962 : NULLIFY (h_block)
2273 962 : NULLIFY (matrix_p)
2274 962 : NULLIFY (matrix_s)
2275 962 : NULLIFY (l)
2276 962 : NULLIFY (last_sgf)
2277 962 : NULLIFY (nshell)
2278 962 : NULLIFY (orb_basis_set)
2279 962 : NULLIFY (p_block)
2280 962 : NULLIFY (particle_set)
2281 962 : NULLIFY (rho)
2282 962 : NULLIFY (s_block)
2283 962 : NULLIFY (sm_h)
2284 962 : NULLIFY (sm_p)
2285 962 : NULLIFY (sm_s)
2286 962 : NULLIFY (w_block)
2287 962 : 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 962 : para_env=para_env)
2297 :
2298 962 : CPASSERT(ASSOCIATED(atomic_kind_set))
2299 962 : CPASSERT(ASSOCIATED(dft_control))
2300 962 : CPASSERT(ASSOCIATED(energy))
2301 962 : CPASSERT(ASSOCIATED(particle_set))
2302 962 : CPASSERT(ASSOCIATED(rho))
2303 :
2304 962 : IF (orthonormal_basis) THEN
2305 962 : NULLIFY (sm_s)
2306 : ELSE
2307 : ! Get overlap matrix in sparse format
2308 : CALL get_qs_env(qs_env=qs_env, &
2309 962 : matrix_s_kp=matrix_s)
2310 962 : CPASSERT(ASSOCIATED(matrix_s))
2311 : END IF
2312 :
2313 : ! Get density matrices in sparse format
2314 :
2315 962 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2316 :
2317 962 : energy%dft_plus_u = 0.0_dp
2318 :
2319 962 : nspin = dft_control%nspins
2320 962 : nimg = dft_control%nimages
2321 :
2322 962 : IF (nspin == 2) THEN
2323 : fspin = 1.0_dp
2324 : ELSE
2325 576 : 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 962 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
2332 962 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
2333 :
2334 962 : nkind = SIZE(atomic_kind_set)
2335 :
2336 2886 : ALLOCATE (first_sgf_atom(natom))
2337 962 : first_sgf_atom(:) = 0
2338 :
2339 : CALL get_particle_set(particle_set, qs_kind_set, &
2340 962 : first_sgf=first_sgf_atom)
2341 :
2342 2886 : ALLOCATE (trps(nsgf))
2343 962 : trps(:) = 0.0_dp
2344 :
2345 962 : IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
2346 2664 : ALLOCATE (dEdq(nsgf))
2347 888 : just_energy = .FALSE.
2348 : ELSE
2349 : just_energy = .TRUE.
2350 : END IF
2351 :
2352 : ! Loop over all spins
2353 :
2354 2310 : DO ispin = 1, nspin
2355 :
2356 1348 : IF (.NOT. just_energy) dEdq(:) = 0.0_dp
2357 :
2358 : ! Calculate Trace(P*S) assuming symmetric matrices
2359 :
2360 1348 : trps(:) = 0.0_dp
2361 :
2362 78006 : DO ic = 1, nimg
2363 76658 : IF (orthonormal_basis) THEN
2364 : NULLIFY (sm_s)
2365 : ELSE
2366 76658 : sm_s => matrix_s(1, ic)%matrix
2367 : END IF
2368 76658 : sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
2369 :
2370 76658 : CALL dbcsr_iterator_start(iter, sm_p)
2371 :
2372 306632 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2373 :
2374 229974 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
2375 :
2376 306632 : IF (orthonormal_basis) THEN
2377 :
2378 0 : IF (iatom /= jatom) CYCLE
2379 :
2380 0 : IF (ASSOCIATED(p_block)) THEN
2381 0 : sgf = first_sgf_atom(iatom)
2382 0 : DO isgf = 1, SIZE(p_block, 1)
2383 0 : trps(sgf) = trps(sgf) + p_block(isgf, isgf)
2384 0 : 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 229974 : found=found)
2395 229974 : CPASSERT(ASSOCIATED(s_block))
2396 :
2397 229974 : sgf = first_sgf_atom(jatom)
2398 1686476 : DO jsgf = 1, SIZE(p_block, 2)
2399 15791548 : DO isgf = 1, SIZE(p_block, 1)
2400 15791548 : trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
2401 : END DO
2402 1686476 : sgf = sgf + 1
2403 : END DO
2404 :
2405 229974 : IF (iatom /= jatom) THEN
2406 114987 : sgf = first_sgf_atom(iatom)
2407 1303186 : DO isgf = 1, SIZE(p_block, 1)
2408 7129194 : DO jsgf = 1, SIZE(p_block, 2)
2409 7129194 : trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
2410 : END DO
2411 1303186 : 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 154664 : CALL dbcsr_iterator_stop(iter)
2420 :
2421 : END DO ! cell images
2422 :
2423 1348 : 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 4044 : 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 2696 : 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 2696 : eps_u_ramping=eps_u_ramping)
2452 :
2453 : ! Check, if this atom needs a DFT+U correction
2454 :
2455 2696 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
2456 2696 : IF (.NOT. dft_plus_u_atom) CYCLE
2457 1348 : IF (lu < 0) CYCLE
2458 :
2459 : ! Apply U ramping if requested
2460 :
2461 1348 : IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
2462 0 : IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
2463 0 : u_minus_j = MIN(u_minus_j + u_ramping, u_minus_j_target)
2464 0 : CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
2465 : END IF
2466 0 : IF (should_output .AND. (output_unit > 0)) THEN
2467 : WRITE (UNIT=output_unit, FMT="(T3,A,3X,A,F0.3,A)") &
2468 0 : "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)), &
2469 0 : "U(eff) = ", u_minus_j*evolt, " eV"
2470 : END IF
2471 : END IF
2472 :
2473 1348 : 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 1348 : nshell=nshell)
2483 :
2484 : ! Count the relevant shell blocks of this atomic kind
2485 :
2486 1348 : nsb = 0
2487 3158 : DO iset = 1, nset
2488 9898 : DO ishell = 1, nshell(iset)
2489 8550 : IF (l(ishell, iset) == lu) nsb = nsb + 1
2490 : END DO
2491 : END DO
2492 :
2493 5392 : ALLOCATE (q_ii(nsb, 2*lu + 1))
2494 :
2495 : ! Print headline if requested
2496 :
2497 1348 : IF (should_output .AND. (print_level > low_print_level)) THEN
2498 0 : IF (output_unit > 0) THEN
2499 0 : ALLOCATE (symbol(2*lu + 1))
2500 0 : DO m = -lu, lu
2501 0 : symbol(lu + m + 1) = sgf_symbol(0, lu, m)
2502 : END DO
2503 0 : IF (nspin > 1) THEN
2504 0 : WRITE (UNIT=spin_info, FMT="(A8,I2)") " of spin", ispin
2505 : ELSE
2506 0 : spin_info = ""
2507 : END IF
2508 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
2509 0 : "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ", ikind, &
2510 0 : ": "//TRIM(atomic_kind_name), &
2511 0 : "Atom Shell ", (ADJUSTR(symbol(i)), i=1, 2*lu + 1), " Trace"
2512 0 : DEALLOCATE (symbol)
2513 : END IF
2514 : END IF
2515 :
2516 : ! Loop over all atoms of the current atomic kind
2517 :
2518 2696 : DO iatom = 1, natom_of_kind
2519 :
2520 1348 : atom_a = atom_list(iatom)
2521 :
2522 1348 : 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 1348 : found=found)
2531 :
2532 : ! Calculate E(U) and dE(U)/dq
2533 :
2534 1348 : IF (ASSOCIATED(p_block)) THEN
2535 :
2536 674 : sgf = first_sgf_atom(atom_a)
2537 :
2538 674 : isb = 0
2539 1579 : DO iset = 1, nset
2540 4949 : DO ishell = 1, nshell(iset)
2541 4275 : IF (l(ishell, iset) == lu) THEN
2542 1348 : isb = isb + 1
2543 1348 : i = 0
2544 5392 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
2545 4044 : q = fspin*trps(sgf)
2546 4044 : i = i + 1
2547 4044 : q_ii(isb, i) = q
2548 : energy%dft_plus_u = energy%dft_plus_u + &
2549 4044 : 0.5_dp*u_minus_j*(q - q**2)/fspin
2550 4044 : IF (.NOT. just_energy) THEN
2551 3702 : dEdq(sgf) = dEdq(sgf) + u_minus_j*(0.5_dp - q)
2552 : END IF
2553 5392 : sgf = sgf + 1
2554 : END DO ! next contracted spherical Gaussian function "isgf"
2555 : ELSE
2556 2022 : 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 4044 : IF (should_output .AND. (print_level > low_print_level)) THEN
2566 0 : CALL para_env%sum(q_ii)
2567 0 : IF (output_unit > 0) THEN
2568 0 : DO isb = 1, nsb
2569 : WRITE (UNIT=output_unit, FMT="(T3,I6,2X,I6,2X,10F8.3)") &
2570 0 : 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 0 : "Total", (SUM(q_ii(:, i)), i=1, 2*lu + 1), SUM(q_ii)
2574 0 : 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 6740 : IF (ALLOCATED(q_ii)) THEN
2581 1348 : DEALLOCATE (q_ii)
2582 : END IF
2583 :
2584 : END DO ! next atomic kind "ikind"
2585 :
2586 1348 : IF (.NOT. just_energy) THEN
2587 1234 : CALL para_env%sum(dEdq)
2588 : END IF
2589 :
2590 : ! Add V(i,j)[U] to V(i,j)[DFT]
2591 :
2592 1348 : IF (PRESENT(matrix_h)) THEN
2593 :
2594 76662 : DO ic = 1, nimg
2595 75476 : IF (orthonormal_basis) THEN
2596 : NULLIFY (sm_s)
2597 : ELSE
2598 75476 : sm_s => matrix_s(1, ic)%matrix
2599 : END IF
2600 75476 : sm_h => matrix_h(ispin, ic)%matrix
2601 :
2602 75476 : CALL dbcsr_iterator_start(iter, sm_h)
2603 :
2604 301904 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2605 :
2606 226428 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, h_block)
2607 :
2608 301904 : IF (orthonormal_basis) THEN
2609 :
2610 0 : IF (iatom /= jatom) CYCLE
2611 :
2612 0 : IF (ASSOCIATED(h_block)) THEN
2613 0 : sgf = first_sgf_atom(iatom)
2614 0 : DO isgf = 1, SIZE(h_block, 1)
2615 0 : h_block(isgf, isgf) = h_block(isgf, isgf) + dEdq(sgf)
2616 0 : 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 226428 : found=found)
2629 226428 : CPASSERT(ASSOCIATED(s_block))
2630 :
2631 : ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
2632 :
2633 226428 : sgf = first_sgf_atom(iatom)
2634 :
2635 2264280 : DO isgf = 1, SIZE(h_block, 1)
2636 2037852 : IF (dEdq(sgf) /= 0.0_dp) THEN
2637 679284 : v = 0.5_dp*dEdq(sgf)
2638 5887128 : DO jsgf = 1, SIZE(h_block, 2)
2639 5887128 : h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
2640 : END DO
2641 : END IF
2642 2264280 : sgf = sgf + 1
2643 : END DO
2644 :
2645 226428 : sgf = first_sgf_atom(jatom)
2646 :
2647 1660472 : DO jsgf = 1, SIZE(h_block, 2)
2648 1434044 : IF (dEdq(sgf) /= 0.0_dp) THEN
2649 226428 : v = 0.5_dp*dEdq(sgf)
2650 3169992 : DO isgf = 1, SIZE(h_block, 1)
2651 3169992 : h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
2652 : END DO
2653 : END IF
2654 1660472 : sgf = sgf + 1
2655 : END DO
2656 :
2657 : END IF ! orthonormal basis set
2658 :
2659 : END DO ! Next atom "iatom"
2660 :
2661 152138 : 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 2310 : IF (PRESENT(matrix_w) .AND. (.NOT. orthonormal_basis)) THEN
2672 :
2673 1116 : DO ic = 1, nimg
2674 1068 : sm_s => matrix_s(1, ic)%matrix
2675 1068 : sm_p => matrix_p(ispin, ic)%matrix
2676 1068 : sm_w => matrix_w(ispin, ic)%matrix
2677 :
2678 1068 : CALL dbcsr_iterator_start(iter, sm_p)
2679 :
2680 4272 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2681 :
2682 3204 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
2683 :
2684 : ! Skip the diagonal blocks of the W matrix
2685 :
2686 3204 : 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 1602 : found=found)
2695 1602 : CPASSERT(ASSOCIATED(w_block))
2696 :
2697 : ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
2698 :
2699 1602 : sgf = first_sgf_atom(iatom)
2700 :
2701 18156 : DO isgf = 1, SIZE(w_block, 1)
2702 16554 : IF (dEdq(sgf) /= 0.0_dp) THEN
2703 6408 : v = -0.5_dp*dEdq(sgf)
2704 38448 : DO jsgf = 1, SIZE(w_block, 2)
2705 38448 : w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
2706 : END DO
2707 : END IF
2708 18156 : sgf = sgf + 1
2709 : END DO
2710 :
2711 1602 : sgf = first_sgf_atom(jatom)
2712 :
2713 10680 : DO jsgf = 1, SIZE(w_block, 2)
2714 8010 : IF (dEdq(sgf) /= 0.0_dp) THEN
2715 0 : v = -0.5_dp*dEdq(sgf)
2716 0 : DO isgf = 1, SIZE(w_block, 1)
2717 0 : w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
2718 : END DO
2719 : END IF
2720 11214 : sgf = sgf + 1
2721 : END DO
2722 :
2723 : END DO ! next block node "jatom"
2724 :
2725 2184 : 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 962 : CALL para_env%sum(energy%dft_plus_u)
2736 :
2737 962 : IF (energy%dft_plus_u < 0.0_dp) THEN
2738 : CALL cp_warn(__LOCATION__, &
2739 : "DFT+U energy contribution is negative possibly due "// &
2740 4 : "to unphysical Mulliken charges!")
2741 : END IF
2742 :
2743 : ! Release local work storage
2744 :
2745 962 : IF (ALLOCATED(first_sgf_atom)) THEN
2746 962 : DEALLOCATE (first_sgf_atom)
2747 : END IF
2748 :
2749 962 : IF (ALLOCATED(trps)) THEN
2750 962 : DEALLOCATE (trps)
2751 : END IF
2752 :
2753 962 : IF (ALLOCATED(dEdq)) THEN
2754 888 : DEALLOCATE (dEdq)
2755 : END IF
2756 :
2757 962 : CALL timestop(handle)
2758 :
2759 2886 : 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 772 : 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 772 : 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 772 : INTEGER, ALLOCATABLE :: first_sgf_atom(:)
2791 : INTEGER, DIMENSION(1) :: iloc
2792 772 : INTEGER, DIMENSION(:), POINTER :: atom_list, atom_list2, nshell, orbitals
2793 772 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
2794 : LOGICAL :: debug, dft_plus_u_atom, found, &
2795 : just_energy, smear
2796 772 : 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 772 : REAL(KIND=dp), ALLOCATABLE :: a(:), q_a_matrix(:, :), q_b_matrix(:, :), q_eigval(:), &
2800 772 : q_eigvec(:, :), q_work(:, :), s_inv_matrix(:, :), s_matrix(:, :), tmp_matrix1(:, :), &
2801 772 : 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 772 : TYPE(block3d), ALLOCATABLE :: a_matrix(:, :), d_matrix(:, :), &
2810 772 : e_matrix(:, :), i_matrix(:, :), &
2811 772 : j_matrix(:, :)
2812 772 : TYPE(block4d), ALLOCATABLE :: dd_matrix(:, :), ee_matrix(:, :), &
2813 772 : ff_matrix(:, :)
2814 772 : REAL(KIND=dp), DIMENSION(:), POINTER :: ao_coef
2815 772 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a_block, d_block, e_block, f_block, &
2816 772 : i_block, j_block, q_a_block, q_b_block, &
2817 772 : s_block, tmp_block, v_block, vhxc_block
2818 :
2819 772 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2820 : TYPE(mp_para_env_type), POINTER :: para_env
2821 772 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s, matrix_vhxc
2822 772 : 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 772 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2830 : TYPE(qs_energy_type), POINTER :: energy
2831 772 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2832 772 : 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 772 : CALL timeset(routineN, handle)
2837 :
2838 772 : debug = .FALSE.
2839 772 : smear = .FALSE.
2840 772 : max_scf = -1
2841 772 : eps_scf = 1.0E30_dp
2842 :
2843 772 : NULLIFY (atom_list)
2844 772 : NULLIFY (atom_list2)
2845 772 : NULLIFY (atomic_kind_set)
2846 772 : NULLIFY (qs_kind_set)
2847 772 : NULLIFY (dft_control)
2848 772 : NULLIFY (energy)
2849 772 : NULLIFY (force)
2850 772 : NULLIFY (first_sgf)
2851 772 : NULLIFY (matrix_p)
2852 772 : NULLIFY (matrix_s)
2853 772 : NULLIFY (matrix_vhxc)
2854 772 : NULLIFY (l)
2855 772 : NULLIFY (last_sgf)
2856 772 : NULLIFY (nshell)
2857 772 : NULLIFY (orb_basis_set)
2858 772 : NULLIFY (orbitals)
2859 772 : NULLIFY (ao_coef)
2860 772 : NULLIFY (particle_set)
2861 772 : NULLIFY (q_a_block)
2862 772 : NULLIFY (q_b_block)
2863 772 : NULLIFY (vhxc_block)
2864 772 : NULLIFY (rho)
2865 772 : NULLIFY (scf_env)
2866 772 : NULLIFY (sm_h)
2867 772 : NULLIFY (sm_p)
2868 772 : NULLIFY (sm_q_a)
2869 772 : NULLIFY (sm_q_b)
2870 772 : NULLIFY (sm_s)
2871 772 : NULLIFY (sm_v)
2872 772 : NULLIFY (sm_a)
2873 772 : NULLIFY (sm_b)
2874 772 : NULLIFY (sm_c)
2875 772 : NULLIFY (sm_d)
2876 772 : NULLIFY (sm_e)
2877 772 : NULLIFY (sm_f)
2878 772 : NULLIFY (sm_i)
2879 772 : NULLIFY (sm_j)
2880 772 : NULLIFY (sm_k)
2881 772 : NULLIFY (sm_mid)
2882 772 : NULLIFY (sm_vhxc)
2883 772 : NULLIFY (v_block)
2884 772 : NULLIFY (a_block)
2885 772 : NULLIFY (d_block)
2886 772 : NULLIFY (e_block)
2887 772 : NULLIFY (f_block)
2888 772 : NULLIFY (i_block)
2889 772 : NULLIFY (j_block)
2890 772 : NULLIFY (para_env)
2891 772 : NULLIFY (s_block)
2892 772 : 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 772 : para_env=para_env)
2906 :
2907 772 : CPASSERT(ASSOCIATED(atomic_kind_set))
2908 772 : CPASSERT(ASSOCIATED(dft_control))
2909 772 : CPASSERT(ASSOCIATED(energy))
2910 772 : CPASSERT(ASSOCIATED(matrix_s))
2911 772 : CPASSERT(ASSOCIATED(particle_set))
2912 772 : CPASSERT(ASSOCIATED(rho))
2913 :
2914 772 : sm_s => matrix_s(1)%matrix
2915 772 : CALL qs_rho_get(rho, rho_ao=matrix_p)
2916 :
2917 772 : energy%dft_plus_u = 0.0_dp
2918 :
2919 772 : nspin = dft_control%nspins
2920 772 : IF (dft_control%mtlr_dft_with_perturbation) THEN
2921 396 : IF (.NOT. ASSOCIATED(matrix_vhxc)) THEN
2922 0 : CPABORT("MTLR requires the projected Hxc matrix for every spin channel.")
2923 : END IF
2924 396 : IF (SIZE(matrix_vhxc) /= nspin) THEN
2925 0 : CPABORT("The number of projected Hxc matrices does not match the spin channels.")
2926 : END IF
2927 : END IF
2928 :
2929 772 : IF (nspin == 2) THEN
2930 : fspin = 1.0_dp
2931 : ELSE
2932 0 : 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 772 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
2939 772 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
2940 :
2941 772 : nkind = SIZE(atomic_kind_set)
2942 :
2943 2316 : ALLOCATE (first_sgf_atom(natom))
2944 772 : first_sgf_atom(:) = 0
2945 :
2946 772 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf_atom)
2947 :
2948 772 : IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
2949 : just_energy = .FALSE.
2950 : ELSE
2951 116 : just_energy = .TRUE.
2952 : END IF
2953 :
2954 772 : CALL dbcsr_init_p(sm_mid)
2955 772 : CALL dbcsr_create(sm_mid, template=sm_s, matrix_type="N")
2956 :
2957 : ! Create local block diagonal matrices
2958 :
2959 772 : ALLOCATE (sm_q_a)
2960 772 : CALL dbcsr_get_block_diag(sm_s, sm_q_a)
2961 :
2962 772 : ALLOCATE (sm_q_b)
2963 772 : CALL dbcsr_get_block_diag(sm_s, sm_q_b)
2964 :
2965 772 : CALL dbcsr_set(sm_q_a, 0.0_dp)
2966 772 : 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 2316 : DO ispin = 1, nspin
2971 :
2972 1544 : sm_p => matrix_p(ispin)%matrix
2973 :
2974 2316 : IF (nspin == 2) THEN
2975 1544 : IF (ispin == 1) THEN
2976 772 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2977 772 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_a)
2978 : ELSE
2979 772 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2980 772 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_b)
2981 : END IF
2982 : ELSE
2983 0 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2984 0 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_a)
2985 0 : CALL dbcsr_copy(sm_q_b, sm_q_a)
2986 : END IF
2987 :
2988 : END DO
2989 :
2990 772 : ALLOCATE (sm_v)
2991 772 : CALL dbcsr_get_block_diag(sm_s, sm_v)
2992 :
2993 772 : IF (PRESENT(matrix_w)) THEN
2994 :
2995 2 : CALL dbcsr_init_p(sm_b)
2996 2 : CALL dbcsr_create(sm_b, template=matrix_s(1)%matrix, matrix_type="N")
2997 :
2998 2 : CALL dbcsr_init_p(sm_d)
2999 2 : CALL dbcsr_create(sm_d, template=matrix_s(1)%matrix, matrix_type="N")
3000 :
3001 2 : CALL dbcsr_init_p(sm_c)
3002 2 : CALL dbcsr_create(sm_c, template=matrix_s(1)%matrix, matrix_type="N")
3003 :
3004 2 : CALL dbcsr_init_p(sm_e)
3005 2 : CALL dbcsr_create(sm_e, template=matrix_s(1)%matrix, matrix_type="N")
3006 :
3007 2 : CALL dbcsr_init_p(sm_f)
3008 2 : CALL dbcsr_create(sm_f, template=matrix_s(1)%matrix, matrix_type="N")
3009 :
3010 2 : CALL dbcsr_init_p(sm_i)
3011 2 : CALL dbcsr_create(sm_i, template=matrix_s(1)%matrix, matrix_type="N")
3012 :
3013 2 : CALL dbcsr_init_p(sm_j)
3014 2 : CALL dbcsr_create(sm_j, template=matrix_s(1)%matrix, matrix_type="N")
3015 :
3016 2 : CALL dbcsr_init_p(sm_k)
3017 2 : CALL dbcsr_create(sm_k, template=matrix_s(1)%matrix, matrix_type="N")
3018 :
3019 : END IF
3020 :
3021 772 : ALLOCATE (sm_a)
3022 772 : CALL dbcsr_get_block_diag(sm_s, sm_a)
3023 :
3024 : ! Loop over all spins
3025 :
3026 2316 : DO ispin = 1, nspin
3027 :
3028 1544 : IF (PRESENT(matrix_h)) THEN
3029 1308 : sm_h => matrix_h(ispin, 1)%matrix
3030 : ELSE
3031 : NULLIFY (sm_h)
3032 : END IF
3033 :
3034 : IF (PRESENT(matrix_w)) THEN
3035 1544 : sm_w => matrix_w(ispin, 1)%matrix
3036 : ELSE
3037 1544 : NULLIFY (sm_w)
3038 : END IF
3039 :
3040 1544 : CALL dbcsr_set(sm_v, 0.0_dp)
3041 1544 : CALL dbcsr_set(sm_a, 0.0_dp)
3042 :
3043 1544 : IF (dft_control%mtlr_dft_with_perturbation) THEN
3044 792 : IF (ispin == 1) perturbation_strength = 0.9_dp*dft_control%perturbation_strength
3045 792 : IF (ispin == 2) perturbation_strength = 1.1_dp*dft_control%perturbation_strength
3046 792 : IF (.NOT. ASSOCIATED(matrix_vhxc(ispin)%matrix)) THEN
3047 0 : CPABORT("MTLR projected Hxc matrix is not initialized.")
3048 : END IF
3049 792 : sm_vhxc => matrix_vhxc(ispin)%matrix
3050 : END IF
3051 :
3052 : ! Loop over all atomic kinds
3053 :
3054 3568 : 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 2024 : 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 2024 : smear=smear)
3076 :
3077 2024 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
3078 2024 : IF (.NOT. dft_plus_u_atom) CYCLE
3079 1544 : IF (lu < 0) CYCLE
3080 1544 : IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) CYCLE
3081 1544 : IF (.NOT. ASSOCIATED(ao_coef)) THEN
3082 0 : CPABORT("Tensorial DFT+U requires associated ao_coef")
3083 : END IF
3084 :
3085 4632 : ALLOCATE (a(SIZE(ao_coef)))
3086 4632 : a(:) = ao_coef(:)
3087 :
3088 : ! Apply U ramping if requested
3089 :
3090 1544 : IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
3091 0 : IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
3092 0 : u_minus_j = MIN(u_minus_j + u_ramping, u_minus_j_target)
3093 0 : CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
3094 : END IF
3095 :
3096 0 : IF (should_output .AND. (output_unit > 0)) THEN
3097 : WRITE (UNIT=output_unit, FMT="(T3,A,3X,A,F0.3,A)") &
3098 0 : "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)), &
3099 0 : "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 1544 : nshell=nshell)
3109 :
3110 1544 : nsb = 1
3111 1544 : nsbsize = (2*lu + 1)
3112 1544 : n = nsb*nsbsize
3113 :
3114 6176 : ALLOCATE (q_a_matrix(n, n))
3115 4632 : ALLOCATE (q_b_matrix(n, n))
3116 4632 : ALLOCATE (s_matrix(n, n))
3117 4632 : ALLOCATE (s_inv_matrix(n, n))
3118 4632 : ALLOCATE (v_matrix(n, n))
3119 :
3120 1544 : q_a_matrix(:, :) = 0.0_dp
3121 1544 : q_b_matrix(:, :) = 0.0_dp
3122 1544 : s_matrix(:, :) = 0.0_dp
3123 1544 : s_inv_matrix(:, :) = 0.0_dp
3124 1544 : v_matrix(:, :) = 0.0_dp
3125 :
3126 : ! Print headline if requested
3127 :
3128 1544 : IF (should_output .AND. (print_level > low_print_level)) THEN
3129 0 : IF (output_unit > 0) THEN
3130 0 : ALLOCATE (symbol(nsbsize))
3131 0 : DO m = -lu, lu
3132 0 : symbol(lu + m + 1) = sgf_symbol(0, lu, m)
3133 : END DO
3134 :
3135 0 : IF (nspin > 1) THEN
3136 0 : WRITE (UNIT=spin_info, FMT="(A8,I2)") " of spin", ispin
3137 : ELSE
3138 0 : spin_info = ""
3139 : END IF
3140 :
3141 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
3142 0 : "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ", ikind, &
3143 0 : ": "//TRIM(atomic_kind_name), &
3144 0 : "Atom Shell ", (ADJUSTR(symbol(i)), i=1, nsbsize), " Trace"
3145 :
3146 0 : DEALLOCATE (symbol)
3147 : END IF
3148 : END IF
3149 :
3150 1544 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
3151 : dft_control%mtlr_ikind == ikind) THEN
3152 792 : dft_control%vhxc(ispin) = 0.0_dp
3153 792 : dft_control%trq(ispin) = 0.0_dp
3154 2376 : ALLOCATE (vhxc_matrix(n, n))
3155 792 : vhxc_matrix(:, :) = 0.0_dp
3156 : END IF
3157 :
3158 : ! Loop over all atoms of the current atomic kind
3159 :
3160 4152 : DO iatom = 1, natom_of_kind
3161 :
3162 2608 : atom_a = atom_list(iatom)
3163 :
3164 2608 : q_a_matrix(:, :) = 0.0_dp
3165 2608 : q_b_matrix(:, :) = 0.0_dp
3166 2608 : s_matrix(:, :) = 0.0_dp
3167 2608 : s_inv_matrix(:, :) = 0.0_dp
3168 2608 : 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 2608 : 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 2608 : 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 2608 : found=found)
3187 :
3188 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
3189 2608 : 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 792 : found=found)
3196 : END IF
3197 :
3198 2608 : IF (ASSOCIATED(q_a_block) .AND. ASSOCIATED(q_b_block) .AND. ASSOCIATED(s_block)) THEN
3199 :
3200 1304 : isb = 0
3201 2848 : DO iset = 1, nset
3202 6176 : DO ishell = 1, nshell(iset)
3203 3328 : IF (l(ishell, iset) /= lu) CYCLE
3204 2608 : isb = isb + 1
3205 2608 : i = 0
3206 7720 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3207 3568 : i = i + 1
3208 3568 : jsb = 0
3209 11904 : DO jset = 1, nset
3210 20032 : DO jshell = 1, nshell(jset)
3211 11456 : IF (l(jshell, jset) /= lu) CYCLE
3212 7136 : jsb = jsb + 1
3213 7136 : j = 0
3214 25040 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
3215 12896 : j = j + 1
3216 12896 : q_a_matrix(i, j) = q_a_matrix(i, j) + q_a_block(isgf, jsgf)*a(isb)*a(jsb)
3217 12896 : q_b_matrix(i, j) = q_b_matrix(i, j) + q_b_block(isgf, jsgf)*a(isb)*a(jsb)
3218 12896 : 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 12896 : dft_control%mtlr_ikind == ikind .AND. &
3222 11456 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
3223 1584 : 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 1304 : CALL invert_matrix(s_matrix, s_inv_matrix, some_real)
3233 :
3234 59192 : q_a_matrix(:, :) = MATMUL(q_a_matrix, s_inv_matrix)
3235 59192 : 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 1304 : dft_control%mtlr_ikind == ikind .AND. &
3240 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
3241 5148 : vhxc_matrix(:, :) = MATMUL(vhxc_matrix, s_inv_matrix)
3242 792 : DO i = 1, n
3243 396 : dft_control%vhxc(ispin) = dft_control%vhxc(ispin) + vhxc_matrix(i, i)
3244 396 : IF (ispin == 1) dft_control%trq(1) = dft_control%trq(1) + q_a_matrix(i, i)
3245 792 : 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 1304 : IF (ASSOCIATED(orbitals)) THEN
3252 0 : 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 0 : ALLOCATE (orb_occ(nsbsize))
3257 0 : ALLOCATE (q_eigval(n))
3258 0 : ALLOCATE (q_eigvec(n, n))
3259 :
3260 0 : q_eigval(:) = 0.0_dp
3261 0 : q_eigvec(:, :) = 0.0_dp
3262 0 : norb = SIZE(orbitals)
3263 :
3264 0 : IF (ispin == 1) THEN
3265 0 : CALL jacobi(q_a_matrix, q_eigval, q_eigvec)
3266 0 : q_a_matrix(:, :) = 0.0_dp
3267 : ELSE
3268 0 : CALL jacobi(q_b_matrix, q_eigval, q_eigvec)
3269 0 : q_b_matrix(:, :) = 0.0_dp
3270 : END IF
3271 :
3272 0 : DO isb = 1, nsb
3273 0 : trq = 0.0_dp
3274 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
3275 0 : trq = trq + q_eigval(i)
3276 : END DO
3277 :
3278 0 : IF (smear) THEN
3279 0 : occ = trq/REAL(norb, KIND=dp)
3280 : ELSE
3281 0 : occ = 1.0_dp/fspin
3282 : END IF
3283 :
3284 0 : orb_occ(:) = .FALSE.
3285 0 : iloc = MAXLOC(q_eigvec(:, isb*nsbsize))
3286 0 : jsb = INT((iloc(1) - 1)/nsbsize) + 1
3287 0 : i = 0
3288 0 : i0 = (jsb - 1)*nsbsize + 1
3289 0 : iorb = -1000
3290 :
3291 0 : DO j = i0, jsb*nsbsize
3292 0 : i = i + 1
3293 0 : IF (i > norb) THEN
3294 0 : DO m = -lu, lu
3295 0 : IF (.NOT. orb_occ(lu + m + 1)) THEN
3296 0 : iorb = i0 + lu + m
3297 0 : orb_occ(lu + m + 1) = .TRUE.
3298 : END IF
3299 : END DO
3300 : ELSE
3301 0 : iorb = i0 + lu + orbitals(i)
3302 0 : orb_occ(lu + orbitals(i) + 1) = .TRUE.
3303 : END IF
3304 :
3305 0 : CPASSERT(iorb /= -1000)
3306 :
3307 0 : iloc = MAXLOC(q_eigvec(iorb, :))
3308 0 : q_eigval(iloc(1)) = MIN(occ, trq)
3309 :
3310 0 : IF (ispin == 1) THEN
3311 0 : q_a_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1))
3312 : ELSE
3313 0 : q_b_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1))
3314 : END IF
3315 :
3316 0 : trq = trq - q_eigval(iloc(1))
3317 : END DO
3318 : END DO
3319 :
3320 0 : IF (ispin == 1) THEN
3321 0 : q_a_matrix(:, :) = MATMUL(q_a_matrix, TRANSPOSE(q_eigvec))
3322 : ELSE
3323 0 : q_b_matrix(:, :) = MATMUL(q_b_matrix, TRANSPOSE(q_eigvec))
3324 : END IF
3325 :
3326 0 : DEALLOCATE (orb_occ)
3327 0 : DEALLOCATE (q_eigval)
3328 0 : DEALLOCATE (q_eigvec)
3329 : END IF
3330 : END IF
3331 :
3332 1304 : trq = 0.0_dp
3333 1304 : trq2 = 0.0_dp
3334 1304 : trqxq_hund_j = 0.0_dp
3335 :
3336 1304 : IF (ispin == 1) THEN
3337 1544 : DO i = 1, n
3338 892 : trq = trq + q_a_matrix(i, i)
3339 3156 : DO j = 1, n
3340 1612 : trq2 = trq2 + q_a_matrix(i, j)*q_a_matrix(j, i)
3341 1612 : trqxq_hund_j = trqxq_hund_j + q_a_matrix(i, j)*q_b_matrix(j, i)
3342 :
3343 2504 : IF (i == j) THEN
3344 892 : v_matrix(i, i) = u_minus_j*(0.5_dp - fspin*q_a_matrix(i, i))
3345 892 : 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 892 : dft_control%mtlr_ikind == ikind .AND. &
3348 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
3349 198 : v_matrix(i, i) = v_matrix(i, i) + perturbation_strength
3350 : END IF
3351 : ELSE
3352 720 : v_matrix(i, j) = -u_minus_j*fspin*q_a_matrix(j, i)
3353 720 : 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 1544 : DO i = 1, n
3359 892 : trq = trq + q_b_matrix(i, i)
3360 3156 : DO j = 1, n
3361 1612 : trq2 = trq2 + q_b_matrix(i, j)*q_b_matrix(j, i)
3362 1612 : trqxq_hund_j = trqxq_hund_j + q_b_matrix(i, j)*q_a_matrix(j, i)
3363 :
3364 2504 : IF (i == j) THEN
3365 892 : v_matrix(i, i) = u_minus_j*(0.5_dp - fspin*q_b_matrix(i, i))
3366 892 : 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 892 : dft_control%mtlr_ikind == ikind .AND. &
3369 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
3370 198 : v_matrix(i, i) = v_matrix(i, i) + perturbation_strength
3371 : END IF
3372 : ELSE
3373 720 : v_matrix(i, j) = -u_minus_j*fspin*q_b_matrix(j, i)
3374 720 : 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 41496 : v_matrix(:, :) = MATMUL(s_inv_matrix, v_matrix)
3381 :
3382 1304 : trq = fspin*trq
3383 1304 : trq2 = fspin*fspin*trq2
3384 1304 : trqxq_hund_j = fspin*fspin*trqxq_hund_j
3385 :
3386 1304 : energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
3387 1304 : 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 2608 : 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 1188 : found=found)
3398 1188 : 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 1188 : found=found)
3405 1188 : CPASSERT(ASSOCIATED(a_block))
3406 :
3407 1188 : isb = 0
3408 2500 : DO iset = 1, nset
3409 5248 : DO ishell = 1, nshell(iset)
3410 2748 : IF (l(ishell, iset) /= lu) CYCLE
3411 2376 : isb = isb + 1
3412 2376 : i = 0
3413 6560 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3414 2872 : i = i + 1
3415 2872 : jsb = 0
3416 9236 : DO jset = 1, nset
3417 14464 : DO jshell = 1, nshell(jset)
3418 7976 : IF (l(jshell, jset) /= lu) CYCLE
3419 5744 : jsb = jsb + 1
3420 5744 : j = 0
3421 18080 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
3422 8720 : j = j + 1
3423 8720 : v_block(isgf, jsgf) = v_matrix(i, j)*a(isb)*a(jsb)
3424 16696 : 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 6760 : IF (should_output .AND. (print_level > low_print_level)) THEN
3439 0 : IF (ispin == 1) THEN
3440 0 : CALL para_env%sum(q_a_matrix)
3441 : ELSE
3442 0 : CALL para_env%sum(q_b_matrix)
3443 : END IF
3444 :
3445 0 : IF (output_unit > 0) THEN
3446 0 : ALLOCATE (q_work(nsb, nsbsize))
3447 0 : q_work(:, :) = 0.0_dp
3448 :
3449 0 : DO isb = 1, nsb
3450 0 : j = 0
3451 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
3452 0 : j = j + 1
3453 0 : IF (ispin == 1) THEN
3454 0 : q_work(isb, j) = q_a_matrix(i, i)
3455 : ELSE
3456 0 : q_work(isb, j) = q_b_matrix(i, i)
3457 : END IF
3458 : END DO
3459 : END DO
3460 :
3461 0 : DO isb = 1, nsb
3462 : WRITE (UNIT=output_unit, FMT="(T3,I6,2X,I6,2X,10F8.3)") &
3463 0 : 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 0 : "Total", (SUM(q_work(:, i)), i=1, nsbsize), SUM(q_work)
3468 0 : WRITE (UNIT=output_unit, FMT="(A)") ""
3469 :
3470 0 : 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 1544 : IF (ALLOCATED(q_a_matrix)) DEALLOCATE (q_a_matrix)
3550 1544 : IF (ALLOCATED(q_b_matrix)) DEALLOCATE (q_b_matrix)
3551 1544 : IF (ALLOCATED(s_matrix)) DEALLOCATE (s_matrix)
3552 1544 : IF (ALLOCATED(s_inv_matrix)) DEALLOCATE (s_inv_matrix)
3553 1544 : IF (ALLOCATED(v_matrix)) DEALLOCATE (v_matrix)
3554 :
3555 1544 : DEALLOCATE (a)
3556 :
3557 1544 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
3558 5112 : dft_control%mtlr_ikind == ikind) THEN
3559 792 : dft_control%vhxc(ispin) = dft_control%vhxc(ispin)/REAL(n, dp)
3560 792 : CALL para_env%sum(dft_control%vhxc(ispin))
3561 792 : CALL para_env%sum(dft_control%trq(ispin))
3562 792 : DEALLOCATE (vhxc_matrix)
3563 : END IF
3564 :
3565 : END DO
3566 :
3567 : ! Add V(i,j)[U] to V(i,j)[DFT]
3568 :
3569 2316 : IF (ASSOCIATED(sm_h)) THEN
3570 1308 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_v, 0.0_dp, sm_mid)
3571 1308 : 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 772 : CALL dbcsr_deallocate_matrix(sm_mid)
3577 :
3578 772 : IF (PRESENT(matrix_w)) THEN
3579 :
3580 2 : CALL get_qs_env(qs_env=qs_env, force=force)
3581 6 : DO ikind = 1, nkind
3582 30 : force(ikind)%tensorial_u(:, :) = 0.0_dp
3583 : END DO
3584 2 : IF (nspin == 2) THEN
3585 2 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_b)
3586 2 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, matrix_p(2)%matrix, 0.0_dp, sm_c)
3587 : ELSE
3588 0 : CALL dbcsr_multiply("N", "N", 0.5_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_b)
3589 0 : CALL dbcsr_multiply("N", "N", 0.5_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_c)
3590 : END IF
3591 2 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_b, sm_s, 0.0_dp, sm_i)
3592 2 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_c, sm_s, 0.0_dp, sm_j)
3593 :
3594 8 : DO jj = 1, 3
3595 6 : CALL dbcsr_desymmetrize(matrixkp_s(jj + 1, 1)%matrix, sm_k)
3596 :
3597 6 : CALL dbcsr_set(sm_f, 0.0_dp)
3598 :
3599 6 : CALL dbcsr_iterator_start(iter, sm_k)
3600 33 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3601 27 : NULLIFY (tmp_block)
3602 27 : CALL dbcsr_iterator_next_block(iter, jatom, katom, tmp_block)
3603 27 : IF (jatom == katom) CYCLE
3604 27 : CALL dbcsr_put_block(matrix=sm_f, row=jatom, col=katom, block=tmp_block)
3605 : END DO
3606 6 : CALL dbcsr_iterator_stop(iter)
3607 6 : CALL dbcsr_finalize(sm_f)
3608 :
3609 6 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_b, sm_f, 0.0_dp, sm_d)
3610 6 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_c, sm_f, 0.0_dp, sm_e)
3611 :
3612 60 : ALLOCATE (a_matrix(nkind, nkind))
3613 54 : ALLOCATE (d_matrix(nkind, nkind))
3614 54 : ALLOCATE (e_matrix(nkind, nkind))
3615 54 : ALLOCATE (i_matrix(nkind, nkind))
3616 54 : ALLOCATE (j_matrix(nkind, nkind))
3617 60 : ALLOCATE (dd_matrix(nkind, nkind))
3618 54 : ALLOCATE (ee_matrix(nkind, nkind))
3619 54 : ALLOCATE (ff_matrix(nkind, nkind))
3620 :
3621 18 : DO ikind = 1, nkind
3622 :
3623 : CALL get_atomic_kind(atomic_kind_set(ikind), &
3624 : atom_list=atom_list, &
3625 12 : 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 12 : dft_plus_u_atom=dft_plus_u_atom)
3634 :
3635 12 : IF (.NOT. dft_plus_u_atom) CYCLE
3636 6 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
3637 6 : IF (lu < 0) CYCLE
3638 6 : 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 6 : nshell=nshell)
3646 :
3647 6 : nsb = 0
3648 18 : DO iset = 1, nset
3649 48 : DO ishell = 1, nshell(iset)
3650 42 : IF (l(ishell, iset) == lu) nsb = nsb + 1
3651 : END DO
3652 : END DO
3653 :
3654 6 : nsbsize = (2*lu + 1)
3655 6 : n = nsb*nsbsize
3656 :
3657 30 : ALLOCATE (a_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3658 24 : ALLOCATE (d_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3659 24 : ALLOCATE (e_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3660 24 : ALLOCATE (i_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3661 24 : ALLOCATE (j_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3662 :
3663 264 : a_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3664 264 : d_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3665 264 : e_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3666 264 : i_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3667 264 : j_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3668 :
3669 36 : DO iatom = 1, natom_of_kind
3670 :
3671 6 : atom_aa = atom_list(iatom)
3672 :
3673 6 : 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 6 : found=found)
3679 :
3680 6 : 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 6 : found=found)
3686 :
3687 6 : 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 6 : found=found)
3693 :
3694 6 : 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 6 : found=found)
3700 :
3701 6 : 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 6 : found=found)
3707 :
3708 : IF (ASSOCIATED(a_block) .AND. ASSOCIATED(d_block) .AND. ASSOCIATED(e_block) &
3709 6 : .AND. ASSOCIATED(i_block) .AND. ASSOCIATED(j_block)) THEN
3710 :
3711 3 : i = 0
3712 9 : DO iset = 1, nset
3713 24 : DO ishell = 1, nshell(iset)
3714 15 : IF (l(ishell, iset) /= lu) CYCLE
3715 30 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3716 18 : i = i + 1
3717 18 : j = 0
3718 69 : DO jset = 1, nset
3719 144 : DO jshell = 1, nshell(jset)
3720 90 : IF (l(jshell, jset) /= lu) CYCLE
3721 180 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
3722 108 : j = j + 1
3723 108 : a_matrix(ikind, ikind)%matrix(i, j, iatom) = a_block(isgf, jsgf)
3724 108 : d_matrix(ikind, ikind)%matrix(i, j, iatom) = -d_block(isgf, jsgf)
3725 108 : e_matrix(ikind, ikind)%matrix(i, j, iatom) = -e_block(isgf, jsgf)
3726 108 : i_matrix(ikind, ikind)%matrix(i, j, iatom) = i_block(isgf, jsgf)
3727 198 : 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 36 : DO jkind = 1, nkind
3738 :
3739 : CALL get_atomic_kind(atomic_kind_set(jkind), &
3740 : atom_list=atom_list2, &
3741 12 : natom=matom_of_kind)
3742 :
3743 12 : CALL get_qs_kind(qs_kind_set(jkind), nsgf=nsgf_kind2)
3744 :
3745 12 : IF (.NOT. ALLOCATED(dd_matrix(ikind, jkind)%matrix)) THEN
3746 72 : ALLOCATE (dd_matrix(ikind, jkind)%matrix(n, nsgf_kind2, natom_of_kind, matom_of_kind))
3747 1014 : dd_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3748 : END IF
3749 :
3750 12 : IF (.NOT. ALLOCATED(ee_matrix(ikind, jkind)%matrix)) THEN
3751 72 : ALLOCATE (ee_matrix(ikind, jkind)%matrix(n, nsgf_kind2, natom_of_kind, matom_of_kind))
3752 1014 : ee_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3753 : END IF
3754 :
3755 12 : IF (.NOT. ALLOCATED(ff_matrix(ikind, jkind)%matrix)) THEN
3756 72 : ALLOCATE (ff_matrix(ikind, jkind)%matrix(nsgf_kind2, n, natom_of_kind, matom_of_kind))
3757 984 : ff_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3758 : END IF
3759 :
3760 48 : DO jatom = 1, matom_of_kind
3761 :
3762 18 : atom_a = atom_list2(jatom)
3763 :
3764 18 : 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 18 : found=found)
3770 :
3771 18 : 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 18 : found=found)
3777 :
3778 18 : IF (ASSOCIATED(d_block) .AND. ASSOCIATED(e_block)) THEN
3779 9 : i = 0
3780 27 : DO iset = 1, nset
3781 72 : DO ishell = 1, nshell(iset)
3782 45 : IF (l(ishell, iset) /= lu) CYCLE
3783 90 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3784 54 : i = i + 1
3785 468 : dd_matrix(ikind, jkind)%matrix(i, :, iatom, jatom) = d_block(isgf, :)
3786 513 : 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 18 : 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 18 : found=found)
3798 :
3799 48 : IF (ASSOCIATED(f_block)) THEN
3800 6 : i = 0
3801 18 : DO iset = 1, nset
3802 48 : DO ishell = 1, nshell(iset)
3803 30 : IF (l(ishell, iset) /= lu) CYCLE
3804 60 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3805 36 : i = i + 1
3806 246 : 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 18 : DO ikind = 1, nkind
3818 :
3819 : CALL get_atomic_kind(atomic_kind_set(ikind), &
3820 : atom_list=atom_list, &
3821 12 : 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 12 : dft_plus_u_atom=dft_plus_u_atom)
3830 :
3831 12 : IF (.NOT. dft_plus_u_atom) CYCLE
3832 6 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
3833 6 : IF (lu < 0) CYCLE
3834 6 : 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 6 : nshell=nshell)
3842 :
3843 6 : CALL para_env%sum(a_matrix(ikind, ikind)%matrix(:, :, :))
3844 6 : CALL para_env%sum(d_matrix(ikind, ikind)%matrix(:, :, :))
3845 6 : CALL para_env%sum(e_matrix(ikind, ikind)%matrix(:, :, :))
3846 6 : CALL para_env%sum(i_matrix(ikind, ikind)%matrix(:, :, :))
3847 6 : CALL para_env%sum(j_matrix(ikind, ikind)%matrix(:, :, :))
3848 :
3849 42 : DO jkind = 1, nkind
3850 12 : CALL para_env%sum(dd_matrix(ikind, jkind)%matrix(:, :, :, :))
3851 12 : CALL para_env%sum(ee_matrix(ikind, jkind)%matrix(:, :, :, :))
3852 24 : CALL para_env%sum(ff_matrix(ikind, jkind)%matrix(:, :, :, :))
3853 : END DO
3854 :
3855 : END DO
3856 :
3857 : ! Force calculation
3858 :
3859 18 : DO ikind = 1, nkind
3860 :
3861 : CALL get_atomic_kind(atomic_kind_set(ikind), &
3862 : atom_list=atom_list, &
3863 12 : 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 12 : dft_plus_u_atom=dft_plus_u_atom)
3872 :
3873 12 : IF (.NOT. dft_plus_u_atom) CYCLE
3874 6 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
3875 6 : IF (lu < 0) CYCLE
3876 6 : 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 6 : nshell=nshell)
3884 :
3885 6 : nsb = 0
3886 18 : DO iset = 1, nset
3887 48 : DO ishell = 1, nshell(iset)
3888 42 : IF (l(ishell, iset) == lu) nsb = nsb + 1
3889 : END DO
3890 : END DO
3891 :
3892 6 : nsbsize = (2*lu + 1)
3893 6 : n = nsb*nsbsize
3894 :
3895 24 : ALLOCATE (tmp_matrix1(n, n))
3896 18 : ALLOCATE (tmp_matrix2(n, n))
3897 18 : ALLOCATE (tmp_matrix3(n, n))
3898 :
3899 12 : DO iatom = 1, natom_of_kind
3900 :
3901 6 : atom_aa = atom_list(iatom)
3902 :
3903 6 : tmp_matrix1(:, :) = MATMUL(d_matrix(ikind, ikind)%matrix(:, :, iatom), &
3904 8034 : a_matrix(ikind, ikind)%matrix(:, :, iatom))
3905 6 : tmp_matrix2(:, :) = MATMUL(e_matrix(ikind, ikind)%matrix(:, :, iatom), &
3906 8034 : a_matrix(ikind, ikind)%matrix(:, :, iatom))
3907 :
3908 6 : i_matrix(ikind, ikind)%matrix(:, :, iatom) = MATMUL(i_matrix(ikind, ikind)%matrix(:, :, iatom), &
3909 8544 : a_matrix(ikind, ikind)%matrix(:, :, iatom))
3910 6 : j_matrix(ikind, ikind)%matrix(:, :, iatom) = MATMUL(j_matrix(ikind, ikind)%matrix(:, :, iatom), &
3911 8544 : a_matrix(ikind, ikind)%matrix(:, :, iatom))
3912 :
3913 : tmp_matrix3(:, :) = 0.5_dp*u_minus_j*tmp_matrix1 - &
3914 6 : u_minus_j*MATMUL(tmp_matrix1, i_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3915 11154 : + 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 24 : u_minus_j*MATMUL(tmp_matrix2, j_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3918 11154 : + hund_j*MATMUL(tmp_matrix1, j_matrix(ikind, ikind)%matrix(:, :, iatom))
3919 :
3920 42 : DO i = 1, n
3921 42 : force(ikind)%tensorial_u(jj, iatom) = force(ikind)%tensorial_u(jj, iatom) + tmp_matrix3(i, i)
3922 : END DO
3923 :
3924 24 : DO jkind = 1, nkind
3925 :
3926 : CALL get_atomic_kind(atomic_kind_set(jkind), &
3927 : atom_list=atom_list2, &
3928 12 : natom=matom_of_kind)
3929 :
3930 12 : CALL get_qs_kind(qs_kind_set(jkind), nsgf=nsgf_kind2)
3931 :
3932 36 : DO jatom = 1, matom_of_kind
3933 :
3934 18 : atom_a = atom_list2(jatom)
3935 :
3936 18 : tmp_matrix1(:, :) = MATMUL(dd_matrix(ikind, jkind)%matrix(:, :, iatom, jatom), &
3937 30762 : ff_matrix(ikind, jkind)%matrix(:, :, iatom, jatom))
3938 24894 : tmp_matrix1(:, :) = MATMUL(tmp_matrix1, a_matrix(ikind, ikind)%matrix(:, :, iatom))
3939 :
3940 18 : tmp_matrix2(:, :) = MATMUL(ee_matrix(ikind, jkind)%matrix(:, :, iatom, jatom), &
3941 30762 : ff_matrix(ikind, jkind)%matrix(:, :, iatom, jatom))
3942 24894 : tmp_matrix2(:, :) = MATMUL(tmp_matrix2, a_matrix(ikind, ikind)%matrix(:, :, iatom))
3943 :
3944 : tmp_matrix3(:, :) = 0.5_dp*u_minus_j*tmp_matrix1 - &
3945 18 : u_minus_j*MATMUL(tmp_matrix1, i_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3946 33462 : + 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 72 : u_minus_j*MATMUL(tmp_matrix2, j_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3949 33462 : + hund_j*MATMUL(tmp_matrix1, j_matrix(ikind, ikind)%matrix(:, :, iatom))
3950 :
3951 30 : IF (atom_a /= atom_aa) THEN
3952 84 : DO i = 1, n
3953 : force(jkind)%tensorial_u(jj, jatom) = force(jkind)%tensorial_u(jj, jatom) &
3954 84 : + tmp_matrix3(i, i)
3955 : END DO
3956 : END IF
3957 :
3958 : END DO
3959 : END DO
3960 : END DO
3961 :
3962 6 : DEALLOCATE (tmp_matrix1)
3963 6 : DEALLOCATE (tmp_matrix2)
3964 36 : DEALLOCATE (tmp_matrix3)
3965 :
3966 : END DO
3967 :
3968 18 : DO ikind = 1, nkind
3969 :
3970 12 : IF (ALLOCATED(a_matrix(ikind, ikind)%matrix)) DEALLOCATE (a_matrix(ikind, ikind)%matrix)
3971 12 : IF (ALLOCATED(d_matrix(ikind, ikind)%matrix)) DEALLOCATE (d_matrix(ikind, ikind)%matrix)
3972 12 : IF (ALLOCATED(e_matrix(ikind, ikind)%matrix)) DEALLOCATE (e_matrix(ikind, ikind)%matrix)
3973 12 : IF (ALLOCATED(i_matrix(ikind, ikind)%matrix)) DEALLOCATE (i_matrix(ikind, ikind)%matrix)
3974 12 : IF (ALLOCATED(j_matrix(ikind, ikind)%matrix)) DEALLOCATE (j_matrix(ikind, ikind)%matrix)
3975 :
3976 42 : DO jkind = 1, nkind
3977 24 : IF (ALLOCATED(dd_matrix(ikind, jkind)%matrix)) DEALLOCATE (dd_matrix(ikind, jkind)%matrix)
3978 24 : IF (ALLOCATED(ee_matrix(ikind, jkind)%matrix)) DEALLOCATE (ee_matrix(ikind, jkind)%matrix)
3979 36 : IF (ALLOCATED(ff_matrix(ikind, jkind)%matrix)) DEALLOCATE (ff_matrix(ikind, jkind)%matrix)
3980 : END DO
3981 :
3982 : END DO
3983 :
3984 30 : DEALLOCATE (a_matrix)
3985 30 : DEALLOCATE (d_matrix)
3986 30 : DEALLOCATE (e_matrix)
3987 30 : DEALLOCATE (i_matrix)
3988 30 : DEALLOCATE (j_matrix)
3989 30 : DEALLOCATE (dd_matrix)
3990 30 : DEALLOCATE (ee_matrix)
3991 38 : DEALLOCATE (ff_matrix)
3992 :
3993 : END DO
3994 :
3995 6 : DO ikind = 1, nkind
3996 30 : force(ikind)%tensorial_u(:, :) = force(ikind)%tensorial_u(:, :)*2.0_dp
3997 : END DO
3998 :
3999 2 : CALL DATE_AND_TIME(values=timevalues)
4000 :
4001 : END IF
4002 :
4003 : ! Collect the energy contributions from all processes
4004 :
4005 772 : CALL para_env%sum(energy%dft_plus_u)
4006 :
4007 772 : 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 0 : "if this warning persists or try a different method!")
4012 : END IF
4013 :
4014 : ! Release local sparse matrices
4015 :
4016 772 : CALL dbcsr_deallocate_matrix(sm_q_a)
4017 772 : CALL dbcsr_deallocate_matrix(sm_q_b)
4018 772 : CALL dbcsr_deallocate_matrix(sm_v)
4019 772 : CALL dbcsr_deallocate_matrix(sm_a)
4020 :
4021 772 : IF (PRESENT(matrix_w)) THEN
4022 2 : CALL dbcsr_deallocate_matrix(sm_b)
4023 2 : CALL dbcsr_deallocate_matrix(sm_c)
4024 2 : CALL dbcsr_deallocate_matrix(sm_d)
4025 2 : CALL dbcsr_deallocate_matrix(sm_e)
4026 2 : CALL dbcsr_deallocate_matrix(sm_f)
4027 2 : CALL dbcsr_deallocate_matrix(sm_k)
4028 2 : CALL dbcsr_deallocate_matrix(sm_i)
4029 2 : CALL dbcsr_deallocate_matrix(sm_j)
4030 : END IF
4031 :
4032 772 : CALL timestop(handle)
4033 :
4034 2316 : END SUBROUTINE tensorial
4035 :
4036 60 : END MODULE dft_plus_u
|
