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 atomic_kind_types, ONLY: atomic_kind_type,&
20 : get_atomic_kind,&
21 : get_atomic_kind_set
22 : USE basis_set_types, ONLY: get_gto_basis_set,&
23 : gto_basis_set_type
24 : USE bibliography, ONLY: Dudarev1997,&
25 : Dudarev1998,&
26 : cite_reference
27 : USE cp_control_types, ONLY: dft_control_type
28 : USE cp_dbcsr_api, ONLY: &
29 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_finalize, &
30 : dbcsr_get_block_p, dbcsr_init_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
31 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
32 : dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_set, dbcsr_type
33 : USE cp_dbcsr_contrib, ONLY: dbcsr_get_block_diag
34 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
35 : copy_fm_to_dbcsr,&
36 : cp_dbcsr_plus_fm_fm_t,&
37 : cp_dbcsr_sm_fm_multiply
38 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix,&
39 : write_fm_with_basis_info
40 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
41 : cp_fm_schur_product,&
42 : cp_fm_transpose
43 : USE cp_fm_diag, ONLY: choose_eigv_solver
44 : USE cp_fm_struct, ONLY: cp_fm_struct_type
45 : USE cp_fm_types, ONLY: cp_fm_create,&
46 : cp_fm_get_info,&
47 : cp_fm_release,&
48 : cp_fm_set_submatrix,&
49 : cp_fm_type
50 : USE cp_log_handling, ONLY: cp_get_default_logger,&
51 : cp_logger_type
52 : USE cp_output_handling, ONLY: cp_p_file,&
53 : cp_print_key_finished_output,&
54 : cp_print_key_should_output,&
55 : cp_print_key_unit_nr,&
56 : low_print_level
57 : USE input_constants, ONLY: plus_u_lowdin,&
58 : plus_u_mulliken,&
59 : plus_u_mulliken_charges,&
60 : plus_u_tensorial
61 : USE input_section_types, ONLY: section_vals_type
62 : USE kinds, ONLY: default_string_length,&
63 : dp
64 : USE kpoint_methods, ONLY: lowdin_kp_trans
65 : USE kpoint_types, ONLY: kpoint_type
66 : USE mathlib, ONLY: invert_matrix,&
67 : jacobi
68 : USE message_passing, ONLY: mp_para_env_type
69 : USE orbital_symbols, ONLY: sgf_symbol
70 : USE parallel_gemm_api, ONLY: parallel_gemm
71 : USE particle_methods, ONLY: get_particle_set
72 : USE particle_types, ONLY: particle_type
73 : USE physcon, ONLY: evolt
74 : USE qs_energy_types, ONLY: qs_energy_type
75 : USE qs_environment_types, ONLY: get_qs_env,&
76 : qs_environment_type
77 : USE qs_force_types, ONLY: qs_force_type
78 : USE qs_kind_types, ONLY: get_qs_kind,&
79 : get_qs_kind_set,&
80 : qs_kind_type,&
81 : set_qs_kind
82 : USE qs_rho_types, ONLY: qs_rho_get,&
83 : qs_rho_type
84 : USE qs_scf_types, ONLY: qs_scf_env_type
85 : #include "./base/base_uses.f90"
86 :
87 : IMPLICIT NONE
88 :
89 : PRIVATE
90 :
91 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dft_plus_u'
92 :
93 : PUBLIC :: plus_u
94 :
95 : CONTAINS
96 : ! **************************************************************************************************
97 : !> \brief Add the DFT+U contribution to the Hamiltonian matrix.\n
98 : !> Wrapper routine for all "+U" methods
99 : !> \param[in] qs_env Quickstep environment
100 : !> \param[in,out] matrix_h Hamiltonian matrices for each spin
101 : !> \param[in,out] matrix_w Energy weighted density matrices for each spin
102 : !> \date 14.01.2008
103 : !> \author Matthias Krack (MK)
104 : !> \version 1.0
105 : ! **************************************************************************************************
106 2682 : SUBROUTINE plus_u(qs_env, matrix_h, matrix_w)
107 :
108 : TYPE(qs_environment_type), POINTER :: qs_env
109 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
110 : POINTER :: matrix_h, matrix_w
111 :
112 : CHARACTER(LEN=*), PARAMETER :: routineN = 'plus_u'
113 :
114 : INTEGER :: handle, output_unit, print_level
115 : LOGICAL :: orthonormal_basis, should_output
116 : TYPE(cp_logger_type), POINTER :: logger
117 : TYPE(dft_control_type), POINTER :: dft_control
118 : TYPE(section_vals_type), POINTER :: input
119 :
120 2682 : CALL timeset(routineN, handle)
121 :
122 2682 : CPASSERT(ASSOCIATED(qs_env))
123 :
124 2682 : NULLIFY (input, dft_control)
125 :
126 2682 : logger => cp_get_default_logger()
127 :
128 : CALL get_qs_env(qs_env=qs_env, &
129 : input=input, &
130 2682 : dft_control=dft_control)
131 :
132 2682 : CALL cite_reference(Dudarev1997)
133 2682 : CALL cite_reference(Dudarev1998)
134 :
135 : ! Later we could save here some time, if the method in use has this property
136 : ! which then has to be figured out here
137 :
138 2682 : orthonormal_basis = .FALSE.
139 :
140 : ! Setup print control
141 :
142 2682 : print_level = logger%iter_info%print_level
143 : should_output = (BTEST(cp_print_key_should_output(logger%iter_info, input, &
144 : "DFT%PRINT%PLUS_U"), cp_p_file) .AND. &
145 2682 : (.NOT. PRESENT(matrix_w)))
146 : output_unit = cp_print_key_unit_nr(logger, input, "DFT%PRINT%PLUS_U", &
147 : extension=".plus_u", &
148 : ignore_should_output=should_output, &
149 2682 : log_filename=.FALSE.)
150 :
151 : ! Select DFT+U method
152 :
153 2682 : SELECT CASE (dft_control%plus_u_method_id)
154 : CASE (plus_u_lowdin)
155 : IF (orthonormal_basis) THEN
156 : ! For an orthonormal basis the Lowdin method and the Mulliken method
157 : ! are equivalent
158 : CALL mulliken(qs_env, orthonormal_basis, matrix_h, &
159 : should_output, output_unit, print_level)
160 : ELSE
161 : CALL lowdin(qs_env, matrix_h, matrix_w, &
162 430 : should_output, output_unit, print_level)
163 : END IF
164 : CASE (plus_u_mulliken)
165 : CALL mulliken(qs_env, orthonormal_basis, matrix_h, &
166 1290 : should_output, output_unit, print_level)
167 : CASE (plus_u_mulliken_charges)
168 : CALL mulliken_charges(qs_env, orthonormal_basis, matrix_h, matrix_w, &
169 404 : should_output, output_unit, print_level)
170 : CASE (plus_u_tensorial)
171 : CALL tensorial(qs_env, matrix_h, matrix_w, &
172 558 : should_output, output_unit, print_level)
173 : CASE DEFAULT
174 2682 : CPABORT("Invalid DFT+U method requested")
175 : END SELECT
176 :
177 : CALL cp_print_key_finished_output(output_unit, logger, input, "DFT%PRINT%PLUS_U", &
178 2682 : ignore_should_output=should_output)
179 :
180 2682 : CALL timestop(handle)
181 :
182 2682 : END SUBROUTINE plus_u
183 :
184 : ! **************************************************************************************************
185 : !> \brief Add a DFT+U contribution to the Hamiltonian matrix\n
186 : !> using a method based on Lowdin charges
187 : !> \f[Q = S^{1/2} P S^{1/2}\f]
188 : !> where \b P and \b S are the density and the
189 : !> overlap matrix, respectively.
190 : !> \param[in] qs_env Quickstep environment
191 : !> \param[in,out] matrix_h Hamiltonian matrices for each spin
192 : !> \param[in,out] matrix_w Energy weighted density matrices for each spin
193 : !> \param should_output ...
194 : !> \param output_unit ...
195 : !> \param print_level ...
196 : !> \date 02.07.2008
197 : !> \par
198 : !> \f{eqnarray*}{
199 : !> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}
200 : !> & = & E^{\rm DFT} + \frac{1}{2}(U - J)\sum_\mu (q_\mu - q_\mu^2)\\[1ex]
201 : !> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
202 : !> & = & \frac{\partial E^{\rm DFT}}
203 : !> {\partial P_{\mu\nu}} +
204 : !> \frac{\partial E^{\rm U}}
205 : !> {\partial P_{\mu\nu}}\\\
206 : !> & = & H_{\mu\nu} +
207 : !> \frac{\partial E^{\rm U}}{\partial q_\mu}
208 : !> \frac{\partial q_\mu}{\partial P_{\mu\nu}}\\\
209 : !> \f}
210 : !> \author Matthias Krack (MK)
211 : !> \version 1.0
212 : ! **************************************************************************************************
213 430 : SUBROUTINE lowdin(qs_env, matrix_h, matrix_w, should_output, output_unit, &
214 : print_level)
215 :
216 : TYPE(qs_environment_type), POINTER :: qs_env
217 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
218 : POINTER :: matrix_h, matrix_w
219 : LOGICAL, INTENT(IN) :: should_output
220 : INTEGER, INTENT(IN) :: output_unit, print_level
221 :
222 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lowdin'
223 :
224 : CHARACTER(LEN=10) :: spin_info
225 430 : CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: symbol
226 : CHARACTER(LEN=default_string_length) :: atomic_kind_name
227 : INTEGER :: atom_a, handle, i, i0, iatom, ikind, iorb, isb, iset, isgf, ishell, ispin, j, &
228 : jsb, jset, jsgf, jshell, lu, m, max_scf, n, natom, natom_of_kind, nimg, nkind, norb, nsb, &
229 : nsbsize, nset, nsgf, nsgf_kind, nspin
230 430 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
231 : INTEGER, DIMENSION(1) :: iloc
232 430 : INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals
233 430 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
234 : LOGICAL :: debug, dft_plus_u_atom, do_kpoints, &
235 : found, if_plus_j, just_energy, smear
236 430 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: orb_occ
237 : REAL(KIND=dp) :: eps_scf, eps_u_ramping, fspin, hund_j, occ, perturbation_strength, sij, &
238 : trq, trq2, trqxq_hund_j, u_minus_j, u_minus_j_target, u_ramping
239 430 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigval, q_eigval
240 430 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orbq, q_eigvec, q_matrix, q_opp_matrix, &
241 430 : q_work, slam
242 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
243 430 : POINTER :: local_data
244 430 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: q_block, q_opp_block, v_block, vhxc_block
245 430 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
246 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
247 : TYPE(cp_fm_type) :: fm_sev, fm_work1, fm_work2, slambda
248 430 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fm_wmat
249 : TYPE(cp_fm_type), POINTER :: fm_s_half
250 430 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vhxc
251 430 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
252 : TYPE(dbcsr_type) :: sm_q, sm_q_opp, sm_v
253 : TYPE(dbcsr_type), POINTER :: sm_h, sm_p, sm_p_opp, sm_s, sm_vhxc, sm_w
254 : TYPE(dft_control_type), POINTER :: dft_control
255 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
256 : TYPE(kpoint_type), POINTER :: kpoints
257 : TYPE(mp_para_env_type), POINTER :: para_env
258 430 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
259 : TYPE(qs_energy_type), POINTER :: energy
260 430 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
261 : TYPE(qs_rho_type), POINTER :: rho
262 : TYPE(qs_scf_env_type), POINTER :: scf_env
263 :
264 430 : CALL timeset(routineN, handle)
265 :
266 430 : debug = .FALSE. ! Set to .TRUE. to print debug information
267 :
268 430 : NULLIFY (sm_h, sm_p, sm_p_opp, sm_s, sm_w, vhxc_block, sm_vhxc, matrix_vhxc)
269 :
270 430 : smear = .FALSE.
271 430 : max_scf = -1
272 430 : eps_scf = 1.0E30_dp
273 430 : if_plus_j = .FALSE.
274 :
275 : CALL get_qs_env(qs_env=qs_env, &
276 : atomic_kind_set=atomic_kind_set, &
277 : qs_kind_set=qs_kind_set, &
278 : dft_control=dft_control, &
279 : do_kpoints=do_kpoints, &
280 : kpoints=kpoints, &
281 : energy=energy, &
282 : matrix_s_kp=matrix_s, &
283 : matrix_vhxc=matrix_vhxc, &
284 : particle_set=particle_set, &
285 : rho=rho, &
286 : scf_env=scf_env, &
287 430 : para_env=para_env)
288 :
289 430 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
290 :
291 430 : energy%dft_plus_u = 0.0_dp
292 :
293 430 : nspin = dft_control%nspins
294 430 : nimg = dft_control%nimages
295 430 : IF (dft_control%mtlr_dft_with_perturbation) THEN
296 222 : IF (.NOT. ASSOCIATED(matrix_vhxc)) THEN
297 0 : CPABORT("MTLR requires the projected Hxc matrix for every spin channel.")
298 : END IF
299 222 : IF (SIZE(matrix_vhxc) /= nspin) THEN
300 0 : CPABORT("The number of projected Hxc matrices does not match the spin channels.")
301 : END IF
302 : END IF
303 :
304 430 : IF (nspin == 2) THEN
305 : fspin = 1.0_dp
306 : ELSE
307 82 : fspin = 0.5_dp
308 : END IF
309 :
310 : ! Get the total number of atoms, contracted spherical Gaussian basis
311 : ! functions, and atomic kinds
312 :
313 430 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
314 430 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
315 :
316 430 : nkind = SIZE(atomic_kind_set)
317 :
318 1290 : ALLOCATE (first_sgf_atom(natom))
319 430 : first_sgf_atom(:) = 0
320 :
321 : CALL get_particle_set(particle_set, qs_kind_set, &
322 430 : first_sgf=first_sgf_atom)
323 :
324 430 : IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
325 : just_energy = .FALSE.
326 : ELSE
327 34 : just_energy = .TRUE.
328 : END IF
329 :
330 : ! Check if unlike spin +J correction is needed and reject an
331 : ! incompatible enforced occupation setup before any expensive work.
332 1040 : DO ikind = 1, nkind
333 610 : NULLIFY (orbitals)
334 : CALL get_qs_kind(qs_kind_set(ikind), &
335 : hund_j=hund_j, &
336 610 : orbitals=orbitals)
337 1040 : IF (hund_j /= 0.0_dp) THEN
338 280 : if_plus_j = .TRUE.
339 280 : IF (ASSOCIATED(orbitals)) THEN
340 : CALL cp_abort(__LOCATION__, "ENFORCE_OCCUPATION is incompatible with a nonzero Hund J "// &
341 0 : "in the Lowdin DFT+U+J method.")
342 : END IF
343 : END IF
344 : END DO
345 :
346 430 : NULLIFY (orbitals)
347 :
348 430 : IF (do_kpoints) THEN
349 0 : fm_wmat => scf_env%scf_work1
350 0 : fmstruct => fm_wmat(1)%matrix_struct
351 : ELSE
352 : ! Retrieve S^(1/2) from the SCF environment
353 430 : fm_s_half => scf_env%s_half
354 430 : CPASSERT(ASSOCIATED(fm_s_half))
355 : ! work matrices
356 430 : CALL cp_fm_get_info(fm_s_half, matrix_struct=fmstruct)
357 : END IF
358 : CALL cp_fm_create(matrix=fm_work1, matrix_struct=fmstruct, &
359 430 : name="FULL WORK MATRIX 1")
360 : CALL cp_fm_create(matrix=fm_work2, matrix_struct=fmstruct, &
361 430 : name="FULL WORK MATRIX 2")
362 :
363 : ! Calculate S eigenvectors and Lambda matrix for forces
364 : ! See sTDA forces (get_lowdin_mo_coefficients in qs_tddfpt2_stda_utils
365 : ! A. Hehn et al JCTC 2022, 18, 4186
366 430 : IF (PRESENT(matrix_w)) THEN
367 2 : IF (do_kpoints) THEN
368 0 : CPABORT("Lowdin forces with k-points NYA in DFT+U")
369 : END IF
370 2 : CALL cp_fm_create(matrix=fm_sev, matrix_struct=fmstruct)
371 2 : CALL cp_fm_create(matrix=slambda, matrix_struct=fmstruct)
372 8 : ALLOCATE (eigval(nsgf), slam(nsgf, 1))
373 2 : sm_s => matrix_s(1, 1)%matrix
374 2 : CALL copy_dbcsr_to_fm(sm_s, fm_work1)
375 2 : CALL choose_eigv_solver(fm_work1, fm_sev, eigval)
376 : !
377 48 : DO i = 1, nsgf
378 48 : IF (eigval(i) > 0._dp) THEN
379 46 : slam(i, 1) = SQRT(eigval(i))
380 : ELSE
381 0 : CPABORT("S matrix not positive definit")
382 : END IF
383 : END DO
384 48 : DO i = 1, nsgf
385 48 : CALL cp_fm_set_submatrix(slambda, slam, 1, i, nsgf, 1, 1.0_dp, 0.0_dp)
386 : END DO
387 48 : DO i = 1, nsgf
388 48 : CALL cp_fm_set_submatrix(slambda, slam, i, 1, 1, nsgf, 1.0_dp, 1.0_dp, .TRUE.)
389 : END DO
390 2 : CALL cp_fm_get_info(slambda, local_data=local_data)
391 48 : DO i = 1, SIZE(local_data, 2)
392 577 : DO j = 1, SIZE(local_data, 1)
393 529 : sij = local_data(j, i)
394 529 : IF (sij > 0.0_dp) sij = 1.0_dp/sij
395 575 : local_data(j, i) = sij
396 : END DO
397 : END DO
398 4 : DEALLOCATE (eigval, slam)
399 : END IF
400 :
401 : ! Calculate S^(1/2)*P*S^(1/2)
402 430 : IF (do_kpoints) THEN
403 0 : CPABORT("Lowdin option with k-points NYA in DFT+U")
404 0 : ALLOCATE (orbq(nsgf, nspin))
405 0 : CALL lowdin_kp_trans(kpoints, orbq)
406 0 : DEALLOCATE (orbq)
407 : END IF
408 :
409 : ! Create local block diagonal matrices
410 430 : sm_s => matrix_s(1, 1)%matrix
411 430 : CALL dbcsr_get_block_diag(sm_s, sm_q)
412 430 : IF (if_plus_j .AND. nspin == 2) THEN
413 280 : CALL dbcsr_get_block_diag(sm_s, sm_q_opp)
414 : END IF
415 430 : CALL dbcsr_get_block_diag(sm_s, sm_v)
416 :
417 : ! Loop over all spins
418 1208 : DO ispin = 1, nspin
419 :
420 778 : CALL dbcsr_set(sm_q, 0.0_dp)
421 778 : IF (if_plus_j .AND. nspin == 2) THEN
422 560 : CALL dbcsr_set(sm_q_opp, 0.0_dp)
423 : END IF
424 778 : CALL dbcsr_set(sm_v, 0.0_dp)
425 :
426 778 : IF (dft_control%mtlr_dft_with_perturbation) THEN
427 444 : IF (ispin == 1) perturbation_strength = 0.9_dp*dft_control%perturbation_strength
428 444 : IF (ispin == 2) perturbation_strength = 1.1_dp*dft_control%perturbation_strength
429 444 : IF (.NOT. ASSOCIATED(matrix_vhxc(ispin)%matrix)) THEN
430 0 : CPABORT("MTLR projected Hxc matrix is not initialized.")
431 : END IF
432 444 : sm_vhxc => matrix_vhxc(ispin)%matrix
433 : END IF
434 :
435 778 : IF (do_kpoints) THEN
436 0 : CPABORT("Lowdin option with k-points NYA in DFT+U")
437 : ELSE
438 : ! Calculate S^(1/2)*P*S^(1/2) as a full matrix (Lowdin)
439 778 : sm_p => matrix_p(ispin, 1)%matrix
440 778 : CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_s_half, fm_work1, nsgf)
441 : CALL parallel_gemm(transa="N", &
442 : transb="N", &
443 : m=nsgf, &
444 : n=nsgf, &
445 : k=nsgf, &
446 : alpha=1.0_dp, &
447 : matrix_a=fm_s_half, &
448 : matrix_b=fm_work1, &
449 : beta=0.0_dp, &
450 778 : matrix_c=fm_work2)
451 : IF (debug) THEN
452 : CALL cp_dbcsr_write_sparse_matrix(sm_p, 4, 6, qs_env, para_env, &
453 : output_unit=output_unit)
454 : CALL write_fm_with_basis_info(fm_s_half, 4, 6, qs_env, para_env, &
455 : output_unit=output_unit)
456 : CALL write_fm_with_basis_info(fm_work2, 4, 6, qs_env, para_env, &
457 : output_unit=output_unit)
458 : END IF ! debug
459 : ! Copy occupation matrix to sparse matrix format, finally we are only
460 : ! interested in the diagonal (atomic) blocks, i.e. the previous full
461 : ! matrix product is not the most efficient choice, anyway.
462 778 : CALL copy_fm_to_dbcsr(fm_work2, sm_q, keep_sparsity=.TRUE.)
463 :
464 : ! If nspin == 1 and one would like to add "unlike" spin +J corrections,
465 : ! sm_p is already enough and sm_p_opp is not needed
466 778 : IF (if_plus_j .AND. nspin == 2) THEN
467 : ! Calculate S^(1/2)*P*S^(1/2) for the other spin channel if
468 : ! unlike spin +J correction is activated.
469 560 : sm_p_opp => matrix_p(3 - ispin, 1)%matrix
470 560 : CALL cp_dbcsr_sm_fm_multiply(sm_p_opp, fm_s_half, fm_work1, nsgf)
471 : CALL parallel_gemm(transa="N", &
472 : transb="N", &
473 : m=nsgf, &
474 : n=nsgf, &
475 : k=nsgf, &
476 : alpha=1.0_dp, &
477 : matrix_a=fm_s_half, &
478 : matrix_b=fm_work1, &
479 : beta=0.0_dp, &
480 560 : matrix_c=fm_work2)
481 560 : CALL copy_fm_to_dbcsr(fm_work2, sm_q_opp, keep_sparsity=.TRUE.)
482 : END IF
483 : END IF
484 :
485 : ! E[DFT+U] = E[DFT] + E[U]
486 : ! = E[DFT] + (U - J)*(Tr(q) - Tr(q*q))/2
487 :
488 : ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U]
489 : ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j)
490 : ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j))
491 :
492 : ! Loop over all atomic kinds
493 1834 : DO ikind = 1, nkind
494 :
495 : ! Load the required atomic kind data
496 : CALL get_atomic_kind(atomic_kind_set(ikind), &
497 : atom_list=atom_list, &
498 : name=atomic_kind_name, &
499 1056 : natom=natom_of_kind)
500 :
501 : CALL get_qs_kind(qs_kind_set(ikind), &
502 : dft_plus_u_atom=dft_plus_u_atom, &
503 : l_of_dft_plus_u=lu, &
504 : nsgf=nsgf_kind, &
505 : basis_set=orb_basis_set, &
506 : u_minus_j=u_minus_j, &
507 : hund_j=hund_j, &
508 : u_minus_j_target=u_minus_j_target, &
509 : u_ramping=u_ramping, &
510 : eps_u_ramping=eps_u_ramping, &
511 : orbitals=orbitals, &
512 : eps_scf=eps_scf, &
513 : max_scf=max_scf, &
514 1056 : smear=smear)
515 :
516 : ! Check, if the atoms of this atomic kind need a DFT+U correction
517 1056 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
518 1056 : IF (.NOT. dft_plus_u_atom) CYCLE
519 778 : IF (lu < 0) CYCLE
520 :
521 : ! Apply U ramping if requested
522 778 : IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
523 0 : IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
524 0 : u_minus_j = MIN(u_minus_j + u_ramping, u_minus_j_target)
525 0 : CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
526 : END IF
527 0 : IF (should_output .AND. (output_unit > 0)) THEN
528 : WRITE (UNIT=output_unit, FMT="(T3,A,3X,A,F0.3,A)") &
529 0 : "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)), &
530 0 : "U(eff) = ", u_minus_j*evolt, " eV"
531 : END IF
532 : END IF
533 :
534 778 : IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) CYCLE
535 :
536 : ! Load the required Gaussian basis set data
537 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
538 : first_sgf=first_sgf, &
539 : l=l, &
540 : last_sgf=last_sgf, &
541 : nset=nset, &
542 778 : nshell=nshell)
543 :
544 : ! Count the relevant shell blocks of this atomic kind
545 778 : nsb = 0
546 1834 : DO iset = 1, nset
547 4224 : DO ishell = 1, nshell(iset)
548 3446 : IF (l(ishell, iset) == lu) nsb = nsb + 1
549 : END DO
550 : END DO
551 :
552 778 : nsbsize = (2*lu + 1)
553 778 : n = nsb*nsbsize
554 :
555 3112 : ALLOCATE (q_matrix(n, n))
556 778 : q_matrix(:, :) = 0.0_dp
557 :
558 778 : IF (hund_j /= 0.0_dp .AND. nspin == 2) THEN
559 1680 : ALLOCATE (q_opp_matrix(n, n))
560 560 : q_opp_matrix(:, :) = 0.0_dp
561 : END IF
562 :
563 778 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
564 : dft_control%mtlr_ikind == ikind) THEN
565 444 : dft_control%vhxc(ispin) = 0.0_dp
566 444 : dft_control%trq(ispin) = 0.0_dp
567 : END IF
568 :
569 : ! Print headline if requested
570 778 : IF (should_output .AND. (print_level > low_print_level)) THEN
571 0 : IF (output_unit > 0) THEN
572 0 : ALLOCATE (symbol(nsbsize))
573 0 : DO m = -lu, lu
574 0 : symbol(lu + m + 1) = sgf_symbol(0, lu, m)
575 : END DO
576 0 : IF (nspin > 1) THEN
577 0 : WRITE (UNIT=spin_info, FMT="(A8,I2)") " of spin", ispin
578 : ELSE
579 0 : spin_info = ""
580 : END IF
581 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
582 0 : "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ", ikind, &
583 0 : ": "//TRIM(atomic_kind_name), &
584 0 : "Atom Shell ", (ADJUSTR(symbol(i)), i=1, nsbsize), " Trace"
585 0 : DEALLOCATE (symbol)
586 : END IF
587 : END IF
588 :
589 : ! Loop over all atoms of the current atomic kind
590 2056 : DO iatom = 1, natom_of_kind
591 1278 : atom_a = atom_list(iatom)
592 1278 : q_matrix(:, :) = 0.0_dp
593 1278 : IF (ALLOCATED(q_opp_matrix)) THEN
594 1060 : q_opp_matrix(:, :) = 0.0_dp
595 : END IF
596 :
597 : ! Get diagonal block
598 : CALL dbcsr_get_block_p(matrix=sm_q, &
599 : row=atom_a, &
600 : col=atom_a, &
601 : block=q_block, &
602 1278 : found=found)
603 1278 : IF (hund_j /= 0.0_dp .AND. nspin == 2) THEN
604 : CALL dbcsr_get_block_p(matrix=sm_q_opp, &
605 : row=atom_a, &
606 : col=atom_a, &
607 : block=q_opp_block, &
608 1060 : found=found)
609 : END IF
610 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
611 1278 : dft_control%mtlr_ikind == ikind .AND. &
612 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
613 : CALL dbcsr_get_block_p(matrix=sm_vhxc, &
614 : row=atom_a, &
615 : col=atom_a, &
616 : block=vhxc_block, &
617 444 : found=found)
618 : END IF
619 :
620 1278 : IF (ASSOCIATED(q_block)) THEN
621 : ! Calculate energy contribution to E(U)
622 639 : i = 0
623 1417 : DO iset = 1, nset
624 3112 : DO ishell = 1, nshell(iset)
625 1695 : IF (l(ishell, iset) /= lu) CYCLE
626 3890 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
627 1834 : i = i + 1
628 1834 : j = 0
629 6197 : DO jset = 1, nset
630 10672 : DO jshell = 1, nshell(jset)
631 6170 : IF (l(jshell, jset) /= lu) CYCLE
632 13340 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
633 7004 : j = j + 1
634 7004 : IF (isgf == jsgf) q_matrix(i, j) = q_block(isgf, jsgf)
635 7004 : IF (hund_j /= 0.0_dp .AND. nspin == 2 .AND. ASSOCIATED(q_opp_block)) THEN
636 3080 : IF (isgf == jsgf) q_opp_matrix(i, j) = q_opp_block(isgf, jsgf)
637 : END IF
638 : IF (ASSOCIATED(vhxc_block) .AND. &
639 : dft_control%mtlr_dft_with_perturbation .AND. &
640 : dft_control%mtlr_ikind == ikind .AND. &
641 7004 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a .AND. &
642 6170 : isgf == jsgf) THEN
643 444 : dft_control%vhxc(ispin) = dft_control%vhxc(ispin) + vhxc_block(isgf, isgf)
644 444 : dft_control%trq(ispin) = dft_control%trq(ispin) + q_block(isgf, isgf)
645 : END IF
646 : END DO ! next contracted spherical Gaussian function "jsgf"
647 : END DO ! next shell "jshell"
648 : END DO ! next shell set "jset"
649 : END DO ! next contracted spherical Gaussian function "isgf"
650 : END DO ! next shell "ishell"
651 : END DO ! next shell set "iset"
652 :
653 : ! Perform the requested manipulations of the (initial) orbital occupations
654 639 : 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 639 : trq = 0.0_dp
710 639 : trq2 = 0.0_dp
711 639 : trqxq_hund_j = 0.0_dp
712 2473 : DO i = 1, n
713 1834 : trq = trq + q_matrix(i, i)
714 9477 : DO j = 1, n
715 7004 : trq2 = trq2 + q_matrix(i, j)*q_matrix(j, i)
716 8838 : IF (hund_j /= 0.0_dp) THEN
717 3080 : IF (nspin == 2) THEN
718 3080 : trqxq_hund_j = trqxq_hund_j + q_matrix(i, j)*q_opp_matrix(j, i)
719 : ELSE
720 0 : trqxq_hund_j = trqxq_hund_j + q_matrix(i, j)*q_matrix(j, i)
721 : END IF
722 : END IF
723 : END DO
724 : END DO
725 639 : trq = fspin*trq
726 639 : trq2 = fspin*fspin*trq2
727 :
728 : ! Calculate energy contribution to E(U)
729 639 : energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
730 639 : IF (hund_j /= 0.0_dp) THEN
731 530 : trqxq_hund_j = fspin*fspin*trqxq_hund_j
732 530 : energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*hund_j*trqxq_hund_j/fspin
733 : END IF
734 :
735 : ! Calculate potential V(U) = dE(U)/dq
736 639 : IF (.NOT. just_energy) THEN
737 : CALL dbcsr_get_block_p(matrix=sm_v, &
738 : row=atom_a, &
739 : col=atom_a, &
740 : block=v_block, &
741 605 : found=found)
742 605 : CPASSERT(ASSOCIATED(v_block))
743 :
744 605 : i = 0
745 1315 : DO iset = 1, nset
746 2840 : DO ishell = 1, nshell(iset)
747 1525 : IF (l(ishell, iset) /= lu) CYCLE
748 3550 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
749 1630 : i = i + 1
750 1630 : j = 0
751 5415 : DO jset = 1, nset
752 9040 : DO jshell = 1, nshell(jset)
753 5150 : IF (l(jshell, jset) /= lu) CYCLE
754 11300 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
755 5780 : j = j + 1
756 5780 : IF (isgf == jsgf) THEN
757 1630 : v_block(isgf, isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i, i))
758 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
759 1630 : dft_control%mtlr_ikind == ikind .AND. &
760 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
761 444 : v_block(isgf, isgf) = v_block(isgf, isgf) + perturbation_strength
762 : END IF
763 : ELSE
764 4150 : CPASSERT(ABS(q_matrix(j, i)) < 1.0E-14_dp)
765 4150 : IF (hund_j /= 0.0_dp .AND. nspin == 2) THEN
766 1480 : CPASSERT(ABS(q_opp_matrix(j, i)) < 1.0E-14_dp)
767 : END IF
768 4150 : v_block(isgf, jsgf) = -u_minus_j*fspin*q_matrix(j, i)
769 : END IF
770 10930 : IF (hund_j /= 0.0_dp) THEN
771 2576 : IF (nspin == 2) THEN
772 2576 : v_block(isgf, jsgf) = v_block(isgf, jsgf) + hund_j*fspin*q_opp_matrix(j, i)
773 : ELSE
774 0 : v_block(isgf, jsgf) = v_block(isgf, jsgf) + hund_j*fspin*q_matrix(j, i)
775 : END IF
776 : END IF
777 : END DO ! next contracted spherical Gaussian function "jsgf"
778 : END DO ! next shell "jshell"
779 : END DO ! next shell set "jset"
780 : END DO ! next contracted spherical Gaussian function "isgf"
781 : END DO ! next shell "ishell"
782 : END DO ! next shell set "iset"
783 : END IF ! not just energy
784 :
785 : END IF ! q_block associated
786 :
787 : ! Consider print requests
788 3334 : IF (should_output .AND. (print_level > low_print_level)) THEN
789 0 : CALL para_env%sum(q_matrix)
790 0 : IF (output_unit > 0) THEN
791 0 : ALLOCATE (q_work(nsb, nsbsize))
792 0 : q_work(:, :) = 0.0_dp
793 0 : DO isb = 1, nsb
794 0 : j = 0
795 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
796 0 : j = j + 1
797 0 : q_work(isb, j) = q_matrix(i, i)
798 : END DO
799 : END DO
800 0 : DO isb = 1, nsb
801 : WRITE (UNIT=output_unit, FMT="(T3,I6,2X,I6,2X,10F8.3)") &
802 0 : atom_a, isb, q_work(isb, :), SUM(q_work(isb, :))
803 : END DO
804 : WRITE (UNIT=output_unit, FMT="(T12,A,2X,10F8.3)") &
805 0 : "Total", (SUM(q_work(:, i)), i=1, nsbsize), SUM(q_work)
806 0 : WRITE (UNIT=output_unit, FMT="(A)") ""
807 0 : DEALLOCATE (q_work)
808 : IF (debug) THEN
809 : ! Print the DFT+U occupation matrix
810 : WRITE (UNIT=output_unit, FMT="(T9,70I10)") (i, i=1, n)
811 : DO i = 1, n
812 : WRITE (UNIT=output_unit, FMT="(T3,I6,70F10.6)") i, q_matrix(i, :)
813 : END DO
814 : ! Print the eigenvalues and eigenvectors of the occupation matrix
815 : ALLOCATE (q_eigval(n))
816 : q_eigval(:) = 0.0_dp
817 : ALLOCATE (q_eigvec(n, n))
818 : q_eigvec(:, :) = 0.0_dp
819 : CALL jacobi(q_matrix, q_eigval, q_eigvec)
820 : WRITE (UNIT=output_unit, FMT="(/,T9,70I10)") (i, i=1, n)
821 : WRITE (UNIT=output_unit, FMT="(T9,71F10.6)") (q_eigval(i), i=1, n), &
822 : SUM(q_eigval(1:n))
823 : DO i = 1, n
824 : WRITE (UNIT=output_unit, FMT="(T3,I6,70F10.6)") i, q_eigvec(i, :)
825 : END DO
826 : DEALLOCATE (q_eigval)
827 : DEALLOCATE (q_eigvec)
828 : END IF ! debug
829 : END IF
830 : IF (debug) THEN
831 : ! Print the full atomic occupation matrix block
832 : ALLOCATE (q_work(nsgf_kind, nsgf_kind))
833 : q_work(:, :) = 0.0_dp
834 : IF (ASSOCIATED(q_block)) q_work(:, :) = q_block(:, :)
835 : CALL para_env%sum(q_work)
836 : IF (output_unit > 0) THEN
837 : norb = SIZE(q_work, 1)
838 : WRITE (UNIT=output_unit, FMT="(/,T9,200I10)") (i, i=1, norb)
839 : DO i = 1, norb
840 : WRITE (UNIT=output_unit, FMT="(T3,I6,200F10.6)") i, q_work(i, :)
841 : END DO
842 : ALLOCATE (q_eigval(norb))
843 : q_eigval(:) = 0.0_dp
844 : ALLOCATE (q_eigvec(norb, norb))
845 : q_eigvec(:, :) = 0.0_dp
846 : CALL jacobi(q_work, q_eigval, q_eigvec)
847 : WRITE (UNIT=output_unit, FMT="(/,T9,200I10)") (i, i=1, norb)
848 : WRITE (UNIT=output_unit, FMT="(T9,201F10.6)") (q_eigval(i), i=1, norb), &
849 : SUM(q_eigval(1:norb))
850 : DO i = 1, norb
851 : WRITE (UNIT=output_unit, FMT="(T3,I6,200F10.6)") i, q_eigvec(i, :)
852 : END DO
853 : DEALLOCATE (q_eigval)
854 : DEALLOCATE (q_eigvec)
855 : END IF
856 : DEALLOCATE (q_work)
857 : END IF ! debug
858 : END IF ! should output
859 :
860 : END DO ! next atom "iatom" of atomic kind "ikind"
861 :
862 778 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
863 : dft_control%mtlr_ikind == ikind) THEN
864 444 : dft_control%vhxc(ispin) = dft_control%vhxc(ispin)/REAL(n, dp)
865 444 : CALL para_env%sum(dft_control%vhxc(ispin))
866 444 : CALL para_env%sum(dft_control%trq(ispin))
867 : END IF
868 :
869 778 : IF (ALLOCATED(q_matrix)) THEN
870 778 : DEALLOCATE (q_matrix)
871 : END IF
872 :
873 3390 : IF (ALLOCATED(q_opp_matrix)) THEN
874 560 : DEALLOCATE (q_opp_matrix)
875 : END IF
876 :
877 : END DO ! next atomic kind "ikind"
878 :
879 : ! Add V(i,j)[U] to V(i,j)[DFT]
880 778 : IF (PRESENT(matrix_h)) THEN
881 708 : IF (do_kpoints) THEN
882 0 : CPABORT("Lowdin option with k-points NYA in DFT+U")
883 : ELSE
884 708 : sm_h => matrix_h(ispin, 1)%matrix
885 708 : CALL cp_dbcsr_sm_fm_multiply(sm_v, fm_s_half, fm_work1, nsgf)
886 708 : CALL cp_fm_transpose(fm_work1, fm_work2)
887 708 : CALL cp_dbcsr_plus_fm_fm_t(sm_h, fm_s_half, fm_work2, nsgf)
888 : END IF
889 : END IF ! An update of the Hamiltonian matrix is requested
890 :
891 : ! Calculate the contribution (non-Pulay part) to the derivatives
892 : ! w.r.t. the nuclear positions
893 1208 : IF (PRESENT(matrix_w)) THEN
894 :
895 2 : sm_p => matrix_p(ispin, 1)%matrix
896 2 : sm_w => matrix_w(ispin, 1)%matrix
897 :
898 2 : CALL cp_dbcsr_sm_fm_multiply(sm_v, fm_s_half, fm_work1, nsgf)
899 2 : CALL cp_fm_transpose(fm_work1, fm_work2)
900 2 : CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_work2, fm_work1, nsgf)
901 2 : CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, 1.0_dp, fm_work1, fm_sev, 0.0_dp, fm_work2)
902 2 : CALL parallel_gemm('T', 'N', nsgf, nsgf, nsgf, 1.0_dp, fm_sev, fm_work2, 0.0_dp, fm_work1)
903 2 : CALL cp_fm_schur_product(fm_work1, slambda, fm_work2)
904 2 : CALL cp_fm_transpose(fm_work2, fm_work1)
905 2 : CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_work1, matrix_b=fm_work2)
906 2 : CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, -2.0_dp, fm_sev, fm_work2, 0.0_dp, fm_work1)
907 2 : CALL cp_dbcsr_plus_fm_fm_t(sm_w, fm_work1, fm_sev, nsgf)
908 :
909 : END IF ! W matrix update requested
910 :
911 : END DO ! next spin "ispin"
912 :
913 430 : IF (PRESENT(matrix_w)) THEN
914 2 : CALL cp_fm_release(matrix=fm_sev)
915 2 : CALL cp_fm_release(matrix=slambda)
916 : END IF
917 :
918 : ! Collect the energy contributions from all processes
919 :
920 430 : CALL para_env%sum(energy%dft_plus_u)
921 :
922 430 : IF (energy%dft_plus_u < 0.0_dp) THEN
923 : CALL cp_warn(__LOCATION__, &
924 : "DFT+U energy contribution is negative possibly due "// &
925 0 : "to unphysical Lowdin charges!")
926 : END IF
927 :
928 : ! Release (local) full matrices
929 430 : NULLIFY (fm_s_half)
930 430 : CALL cp_fm_release(matrix=fm_work1)
931 430 : CALL cp_fm_release(matrix=fm_work2)
932 :
933 : ! Release (local) sparse matrices
934 430 : CALL dbcsr_release(sm_q)
935 430 : CALL dbcsr_release(sm_v)
936 :
937 430 : IF (if_plus_j .AND. nspin == 2) THEN
938 280 : CALL dbcsr_release(sm_q_opp)
939 : END IF
940 :
941 430 : CALL timestop(handle)
942 :
943 1720 : END SUBROUTINE lowdin
944 :
945 : ! **************************************************************************************************
946 : !> \brief Add a DFT+U contribution to the Hamiltonian matrix\n
947 : !> using a method based on the Mulliken population analysis
948 : !> \f[q_{\mu\nu} = \frac{1}{2} (P_{\mu\nu} S_{\nu\mu} +
949 : !> S_{\mu\nu} P_{\nu\mu})\f]
950 : !> where \b P and \b S are the density and the
951 : !> overlap matrix, respectively.
952 : !> \param[in] qs_env Quickstep environment
953 : !> \param orthonormal_basis ...
954 : !> \param[in,out] matrix_h Hamiltonian matrices for each spin
955 : !> \param should_output ...
956 : !> \param output_unit ...
957 : !> \param print_level ...
958 : !> \date 03.07.2008
959 : !> \par
960 : !> \f{eqnarray*}{
961 : !> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}\\\
962 : !> & = & E^{\rm DFT} + \frac{1}{2}\sum_A(U_A - J_A)(Tr(q_A) - Tr(q^2_A))\\[1ex]
963 : !> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
964 : !> & = & \frac{\partial E^{\rm DFT}}
965 : !> {\partial P_{\mu\nu}} +
966 : !> \frac{\partial E^{\rm U}}
967 : !> {\partial P_{\mu\nu}}\\\
968 : !> & = & H_{\mu\nu} + \sum_A
969 : !> \frac{\partial E^{\rm U}}{\partial q_A}
970 : !> \frac{\partial q_A}{\partial P_{\mu\nu}}\\\
971 : !> \f}
972 : !> \author Matthias Krack (MK)
973 : !> \version 1.0
974 : ! **************************************************************************************************
975 1290 : SUBROUTINE mulliken(qs_env, orthonormal_basis, matrix_h, should_output, &
976 : output_unit, print_level)
977 :
978 : TYPE(qs_environment_type), POINTER :: qs_env
979 : LOGICAL, INTENT(IN) :: orthonormal_basis
980 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
981 : POINTER :: matrix_h
982 : LOGICAL, INTENT(IN) :: should_output
983 : INTEGER, INTENT(IN) :: output_unit, print_level
984 :
985 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mulliken'
986 :
987 : CHARACTER(LEN=10) :: spin_info
988 1290 : CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: symbol
989 : CHARACTER(LEN=default_string_length) :: atomic_kind_name
990 : INTEGER :: atom_a, handle, i, i0, iatom, ic, ikind, iorb, isb, iset, isgf, ishell, ispin, j, &
991 : jsb, jset, jsgf, jshell, lu, m, max_scf, n, natom, natom_of_kind, nimg, nkind, norb, nsb, &
992 : nsbsize, nset, nsgf_kind, nspin
993 : INTEGER, DIMENSION(1) :: iloc
994 1290 : INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals
995 1290 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
996 : LOGICAL :: debug, dft_plus_u_atom, found, &
997 : just_energy, occupation_enforced, smear
998 1290 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_plus_u_kind, orb_occ
999 : REAL(KIND=dp) :: eps_scf, eps_u_ramping, fspin, occ, trq, &
1000 : trq2, u_minus_j, u_minus_j_target, &
1001 : u_ramping
1002 1290 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: q_eigval
1003 1290 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: q_eigvec, q_matrix, q_work
1004 1290 : REAL(KIND=dp), DIMENSION(:), POINTER :: nelec
1005 1290 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, q_block, s_block, &
1006 1290 : v_block
1007 1290 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1008 : TYPE(atomic_kind_type), POINTER :: kind_a
1009 1290 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
1010 : TYPE(dbcsr_type), POINTER :: sm_h, sm_p, sm_q, sm_s, sm_v
1011 : TYPE(dft_control_type), POINTER :: dft_control
1012 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1013 : TYPE(mp_para_env_type), POINTER :: para_env
1014 1290 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1015 : TYPE(qs_energy_type), POINTER :: energy
1016 1290 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1017 : TYPE(qs_rho_type), POINTER :: rho
1018 :
1019 1290 : CALL timeset(routineN, handle)
1020 :
1021 1290 : debug = .FALSE. ! Set to .TRUE. to print debug information
1022 :
1023 1290 : NULLIFY (atom_list)
1024 1290 : NULLIFY (atomic_kind_set)
1025 1290 : NULLIFY (qs_kind_set)
1026 1290 : NULLIFY (dft_control)
1027 1290 : NULLIFY (energy)
1028 1290 : NULLIFY (first_sgf)
1029 1290 : NULLIFY (h_block)
1030 1290 : NULLIFY (matrix_p)
1031 1290 : NULLIFY (matrix_s)
1032 1290 : NULLIFY (l)
1033 1290 : NULLIFY (last_sgf)
1034 1290 : NULLIFY (nelec)
1035 1290 : NULLIFY (nshell)
1036 1290 : NULLIFY (orb_basis_set)
1037 1290 : NULLIFY (p_block)
1038 1290 : NULLIFY (particle_set)
1039 1290 : NULLIFY (q_block)
1040 1290 : NULLIFY (rho)
1041 1290 : NULLIFY (s_block)
1042 1290 : NULLIFY (orbitals)
1043 1290 : NULLIFY (sm_h)
1044 1290 : NULLIFY (sm_p)
1045 1290 : NULLIFY (sm_q)
1046 1290 : NULLIFY (sm_s)
1047 1290 : NULLIFY (sm_v)
1048 1290 : NULLIFY (v_block)
1049 1290 : NULLIFY (para_env)
1050 :
1051 1290 : smear = .FALSE.
1052 1290 : max_scf = -1
1053 1290 : eps_scf = 1.0E30_dp
1054 1290 : occupation_enforced = .FALSE.
1055 :
1056 : CALL get_qs_env(qs_env=qs_env, &
1057 : atomic_kind_set=atomic_kind_set, &
1058 : qs_kind_set=qs_kind_set, &
1059 : dft_control=dft_control, &
1060 : energy=energy, &
1061 : particle_set=particle_set, &
1062 : rho=rho, &
1063 1290 : para_env=para_env)
1064 :
1065 1290 : CPASSERT(ASSOCIATED(atomic_kind_set))
1066 1290 : CPASSERT(ASSOCIATED(dft_control))
1067 1290 : CPASSERT(ASSOCIATED(energy))
1068 1290 : CPASSERT(ASSOCIATED(particle_set))
1069 1290 : CPASSERT(ASSOCIATED(rho))
1070 :
1071 1290 : IF (orthonormal_basis) THEN
1072 : NULLIFY (sm_s)
1073 : ELSE
1074 : ! Get overlap matrix in sparse format
1075 : CALL get_qs_env(qs_env=qs_env, &
1076 1290 : matrix_s_kp=matrix_s)
1077 1290 : CPASSERT(ASSOCIATED(matrix_s))
1078 : END IF
1079 1290 : nimg = dft_control%nimages
1080 :
1081 : ! Get density matrices in sparse format
1082 :
1083 1290 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1084 :
1085 1290 : energy%dft_plus_u = 0.0_dp
1086 :
1087 1290 : nspin = dft_control%nspins
1088 :
1089 1290 : IF (nspin == 2) THEN
1090 : fspin = 1.0_dp
1091 : ELSE
1092 660 : fspin = 0.5_dp
1093 : END IF
1094 :
1095 : ! Get the total number of atoms, contracted spherical Gaussian basis
1096 : ! functions, and atomic kinds
1097 :
1098 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1099 1290 : natom=natom)
1100 :
1101 1290 : nkind = SIZE(atomic_kind_set)
1102 :
1103 3870 : ALLOCATE (is_plus_u_kind(nkind))
1104 1290 : is_plus_u_kind(:) = .FALSE.
1105 :
1106 1290 : IF (PRESENT(matrix_h)) THEN
1107 : just_energy = .FALSE.
1108 : ELSE
1109 570 : just_energy = .TRUE.
1110 : END IF
1111 :
1112 : ! Loop over all spins
1113 3210 : DO ispin = 1, nspin
1114 :
1115 : ! Loop over cell images
1116 5130 : DO ic = 1, nimg
1117 1920 : IF (.NOT. orthonormal_basis) THEN
1118 1920 : sm_s => matrix_s(1, ic)%matrix
1119 : END IF
1120 :
1121 1920 : IF (PRESENT(matrix_h)) THEN
1122 : ! Hamiltonian matrix for spin ispin in sparse format
1123 1072 : sm_h => matrix_h(ispin, ic)%matrix
1124 : ELSE
1125 : NULLIFY (sm_h)
1126 : END IF
1127 :
1128 : ! Get density matrix for spin ispin in sparse format
1129 :
1130 1920 : sm_p => matrix_p(ispin, ic)%matrix
1131 :
1132 1920 : IF (.NOT. ASSOCIATED(sm_q)) THEN
1133 1290 : ALLOCATE (sm_q)
1134 1290 : CALL dbcsr_get_block_diag(sm_p, sm_q)
1135 : END IF
1136 1920 : CALL dbcsr_set(sm_q, 0.0_dp)
1137 :
1138 1920 : IF (.NOT. ASSOCIATED(sm_v)) THEN
1139 1290 : ALLOCATE (sm_v)
1140 1290 : CALL dbcsr_get_block_diag(sm_p, sm_v)
1141 : END IF
1142 1920 : CALL dbcsr_set(sm_v, 0.0_dp)
1143 :
1144 7680 : DO iatom = 1, natom
1145 :
1146 : CALL dbcsr_get_block_p(matrix=sm_p, &
1147 : row=iatom, &
1148 : col=iatom, &
1149 : block=p_block, &
1150 5760 : found=found)
1151 :
1152 5760 : IF (.NOT. ASSOCIATED(p_block)) CYCLE
1153 :
1154 : CALL dbcsr_get_block_p(matrix=sm_q, &
1155 : row=iatom, &
1156 : col=iatom, &
1157 : block=q_block, &
1158 2880 : found=found)
1159 2880 : CPASSERT(ASSOCIATED(q_block))
1160 :
1161 13440 : IF (orthonormal_basis) THEN
1162 : ! S is the unit matrix
1163 0 : DO isgf = 1, SIZE(q_block, 1)
1164 0 : q_block(isgf, isgf) = p_block(isgf, isgf)
1165 : END DO
1166 : ELSE
1167 : CALL dbcsr_get_block_p(matrix=sm_s, &
1168 : row=iatom, &
1169 : col=iatom, &
1170 : block=s_block, &
1171 2880 : found=found)
1172 2880 : CPASSERT(ASSOCIATED(s_block))
1173 : ! Exploit that P and S are symmetric
1174 24960 : DO jsgf = 1, SIZE(p_block, 2)
1175 238080 : DO isgf = 1, SIZE(p_block, 1)
1176 232320 : q_block(isgf, jsgf) = p_block(isgf, jsgf)*s_block(isgf, jsgf)
1177 : END DO
1178 : END DO
1179 : END IF ! orthonormal basis set
1180 :
1181 : END DO ! next atom "iatom"
1182 :
1183 : ! E[DFT+U] = E[DFT] + E[U]
1184 : ! = E[DFT] + (U - J)*(Tr(q) - Tr(q*q))/2
1185 :
1186 : ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U]
1187 : ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j)
1188 : ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j))
1189 :
1190 : ! Loop over all atomic kinds
1191 :
1192 5760 : DO ikind = 1, nkind
1193 :
1194 : ! Load the required atomic kind data
1195 :
1196 : CALL get_atomic_kind(atomic_kind_set(ikind), &
1197 : atom_list=atom_list, &
1198 : name=atomic_kind_name, &
1199 3840 : natom=natom_of_kind)
1200 :
1201 : CALL get_qs_kind(qs_kind_set(ikind), &
1202 : dft_plus_u_atom=dft_plus_u_atom, &
1203 : l_of_dft_plus_u=lu, &
1204 : nsgf=nsgf_kind, &
1205 : basis_set=orb_basis_set, &
1206 : u_minus_j=u_minus_j, &
1207 : u_minus_j_target=u_minus_j_target, &
1208 : u_ramping=u_ramping, &
1209 : eps_u_ramping=eps_u_ramping, &
1210 : nelec=nelec, &
1211 : orbitals=orbitals, &
1212 : eps_scf=eps_scf, &
1213 : max_scf=max_scf, &
1214 3840 : smear=smear)
1215 :
1216 : ! Check, if the atoms of this atomic kind need a DFT+U correction
1217 :
1218 3840 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
1219 3840 : IF (.NOT. dft_plus_u_atom) CYCLE
1220 1920 : IF (lu < 0) CYCLE
1221 :
1222 : ! Apply U ramping if requested
1223 :
1224 1920 : IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
1225 976 : IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
1226 464 : u_minus_j = MIN(u_minus_j + u_ramping, u_minus_j_target)
1227 464 : CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
1228 : END IF
1229 976 : IF (should_output .AND. (output_unit > 0)) THEN
1230 : WRITE (UNIT=output_unit, FMT="(T3,A,3X,A,F0.3,A)") &
1231 476 : "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)), &
1232 952 : "U(eff) = ", u_minus_j*evolt, " eV"
1233 : END IF
1234 : END IF
1235 :
1236 1920 : IF (u_minus_j == 0.0_dp) CYCLE
1237 :
1238 1920 : is_plus_u_kind(ikind) = .TRUE.
1239 :
1240 : ! Load the required Gaussian basis set data
1241 :
1242 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1243 : first_sgf=first_sgf, &
1244 : l=l, &
1245 : last_sgf=last_sgf, &
1246 : nset=nset, &
1247 1920 : nshell=nshell)
1248 :
1249 : ! Count the relevant shell blocks of this atomic kind
1250 :
1251 1920 : nsb = 0
1252 5760 : DO iset = 1, nset
1253 15360 : DO ishell = 1, nshell(iset)
1254 13440 : IF (l(ishell, iset) == lu) nsb = nsb + 1
1255 : END DO
1256 : END DO
1257 :
1258 1920 : nsbsize = (2*lu + 1)
1259 1920 : n = nsb*nsbsize
1260 :
1261 7680 : ALLOCATE (q_matrix(n, n))
1262 1920 : q_matrix(:, :) = 0.0_dp
1263 :
1264 : ! Print headline if requested
1265 :
1266 1920 : IF (should_output .AND. (print_level > low_print_level)) THEN
1267 0 : IF (output_unit > 0) THEN
1268 0 : ALLOCATE (symbol(nsbsize))
1269 0 : DO m = -lu, lu
1270 0 : symbol(lu + m + 1) = sgf_symbol(0, lu, m)
1271 : END DO
1272 0 : IF (nspin > 1) THEN
1273 0 : WRITE (UNIT=spin_info, FMT="(A8,I2)") " of spin", ispin
1274 : ELSE
1275 0 : spin_info = ""
1276 : END IF
1277 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
1278 0 : "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ", ikind, &
1279 0 : ": "//TRIM(atomic_kind_name), &
1280 0 : "Atom Shell ", (ADJUSTR(symbol(i)), i=1, nsbsize), " Trace"
1281 0 : DEALLOCATE (symbol)
1282 : END IF
1283 : END IF
1284 :
1285 : ! Loop over all atoms of the current atomic kind
1286 :
1287 3840 : DO iatom = 1, natom_of_kind
1288 :
1289 1920 : atom_a = atom_list(iatom)
1290 :
1291 1920 : q_matrix(:, :) = 0.0_dp
1292 :
1293 : ! Get diagonal block
1294 :
1295 : CALL dbcsr_get_block_p(matrix=sm_q, &
1296 : row=atom_a, &
1297 : col=atom_a, &
1298 : block=q_block, &
1299 1920 : found=found)
1300 :
1301 : ! Calculate energy contribution to E(U)
1302 :
1303 1920 : IF (ASSOCIATED(q_block)) THEN
1304 :
1305 960 : i = 0
1306 2880 : DO iset = 1, nset
1307 7680 : DO ishell = 1, nshell(iset)
1308 4800 : IF (l(ishell, iset) /= lu) CYCLE
1309 9600 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1310 5760 : i = i + 1
1311 5760 : j = 0
1312 22080 : DO jset = 1, nset
1313 46080 : DO jshell = 1, nshell(jset)
1314 28800 : IF (l(jshell, jset) /= lu) CYCLE
1315 57600 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
1316 34560 : j = j + 1
1317 63360 : q_matrix(i, j) = q_block(isgf, jsgf)
1318 : END DO ! next contracted spherical Gaussian function "jsgf"
1319 : END DO ! next shell "jshell"
1320 : END DO ! next shell set "jset"
1321 : END DO ! next contracted spherical Gaussian function "isgf"
1322 : END DO ! next shell "ishell"
1323 : END DO ! next shell set "iset"
1324 :
1325 : ! Perform the requested manipulations of the (initial) orbital occupations
1326 :
1327 960 : IF (ASSOCIATED(orbitals)) THEN
1328 0 : IF ((qs_env%scf_env%iter_delta >= eps_scf) .OR. &
1329 : ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
1330 : (qs_env%scf_env%iter_count <= max_scf))) THEN
1331 0 : ALLOCATE (orb_occ(nsbsize))
1332 0 : ALLOCATE (q_eigval(n))
1333 0 : q_eigval(:) = 0.0_dp
1334 0 : ALLOCATE (q_eigvec(n, n))
1335 0 : q_eigvec(:, :) = 0.0_dp
1336 0 : norb = SIZE(orbitals)
1337 0 : CALL jacobi(q_matrix, q_eigval, q_eigvec)
1338 0 : q_matrix(:, :) = 0.0_dp
1339 0 : IF (nelec(ispin) >= 0.5_dp) THEN
1340 0 : trq = nelec(ispin)/SUM(q_eigval(1:n))
1341 0 : q_eigval(1:n) = trq*q_eigval(1:n)
1342 : END IF
1343 0 : DO isb = 1, nsb
1344 0 : trq = 0.0_dp
1345 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
1346 0 : trq = trq + q_eigval(i)
1347 : END DO
1348 0 : IF (smear) THEN
1349 0 : occ = trq/REAL(norb, KIND=dp)
1350 : ELSE
1351 0 : occ = 1.0_dp/fspin
1352 : END IF
1353 0 : orb_occ(:) = .FALSE.
1354 0 : iloc = MAXLOC(q_eigvec(:, isb*nsbsize))
1355 0 : jsb = INT((iloc(1) - 1)/nsbsize) + 1
1356 0 : i = 0
1357 0 : i0 = (jsb - 1)*nsbsize + 1
1358 0 : iorb = -1000
1359 0 : DO j = i0, jsb*nsbsize
1360 0 : i = i + 1
1361 0 : IF (i > norb) THEN
1362 0 : DO m = -lu, lu
1363 0 : IF (.NOT. orb_occ(lu + m + 1)) THEN
1364 0 : iorb = i0 + lu + m
1365 0 : orb_occ(lu + m + 1) = .TRUE.
1366 : END IF
1367 : END DO
1368 : ELSE
1369 0 : iorb = i0 + lu + orbitals(i)
1370 0 : orb_occ(lu + orbitals(i) + 1) = .TRUE.
1371 : END IF
1372 0 : CPASSERT(iorb /= -1000)
1373 0 : iloc = MAXLOC(q_eigvec(iorb, :))
1374 0 : q_eigval(iloc(1)) = MIN(occ, trq)
1375 0 : q_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1)) ! backtransform left
1376 0 : trq = trq - q_eigval(iloc(1))
1377 : END DO
1378 : END DO
1379 0 : q_matrix(:, :) = MATMUL(q_matrix, TRANSPOSE(q_eigvec)) ! backtransform right
1380 0 : DEALLOCATE (orb_occ)
1381 0 : DEALLOCATE (q_eigval)
1382 0 : DEALLOCATE (q_eigvec)
1383 0 : occupation_enforced = .TRUE.
1384 : END IF
1385 : END IF ! orbitals associated
1386 :
1387 960 : trq = 0.0_dp
1388 960 : trq2 = 0.0_dp
1389 :
1390 6720 : DO i = 1, n
1391 5760 : trq = trq + q_matrix(i, i)
1392 41280 : DO j = 1, n
1393 40320 : trq2 = trq2 + q_matrix(i, j)*q_matrix(j, i)
1394 : END DO
1395 : END DO
1396 :
1397 960 : trq = fspin*trq
1398 960 : trq2 = fspin*fspin*trq2
1399 :
1400 : ! Calculate energy contribution to E(U)
1401 :
1402 960 : energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
1403 :
1404 : ! Calculate potential V(U) = dE(U)/dq
1405 :
1406 960 : IF (.NOT. just_energy) THEN
1407 :
1408 : CALL dbcsr_get_block_p(matrix=sm_v, &
1409 : row=atom_a, &
1410 : col=atom_a, &
1411 : block=v_block, &
1412 536 : found=found)
1413 536 : CPASSERT(ASSOCIATED(v_block))
1414 :
1415 536 : i = 0
1416 1608 : DO iset = 1, nset
1417 4288 : DO ishell = 1, nshell(iset)
1418 2680 : IF (l(ishell, iset) /= lu) CYCLE
1419 5360 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1420 3216 : i = i + 1
1421 3216 : j = 0
1422 12328 : DO jset = 1, nset
1423 25728 : DO jshell = 1, nshell(jset)
1424 16080 : IF (l(jshell, jset) /= lu) CYCLE
1425 32160 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
1426 19296 : j = j + 1
1427 35376 : IF (isgf == jsgf) THEN
1428 3216 : v_block(isgf, isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i, i))
1429 : ELSE
1430 16080 : v_block(isgf, jsgf) = -u_minus_j*fspin*q_matrix(j, i)
1431 : END IF
1432 : END DO ! next contracted spherical Gaussian function "jsgf"
1433 : END DO ! next shell "jshell"
1434 : END DO ! next shell set "jset"
1435 : END DO ! next contracted spherical Gaussian function "isgf"
1436 : END DO ! next shell "ishell"
1437 : END DO ! next shell set "iset"
1438 :
1439 : END IF ! not just energy
1440 :
1441 : END IF ! q_block associated
1442 :
1443 : ! Consider print requests
1444 :
1445 5760 : IF (should_output .AND. (print_level > low_print_level)) THEN
1446 0 : CALL para_env%sum(q_matrix)
1447 0 : IF (output_unit > 0) THEN
1448 0 : ALLOCATE (q_work(nsb, nsbsize))
1449 0 : q_work(:, :) = 0.0_dp
1450 0 : DO isb = 1, nsb
1451 0 : j = 0
1452 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
1453 0 : j = j + 1
1454 0 : q_work(isb, j) = q_matrix(i, i)
1455 : END DO
1456 : END DO
1457 0 : DO isb = 1, nsb
1458 : WRITE (UNIT=output_unit, FMT="(T3,I6,2X,I6,2X,10F8.3)") &
1459 0 : atom_a, isb, q_work(isb, :), SUM(q_work(isb, :))
1460 : END DO
1461 : WRITE (UNIT=output_unit, FMT="(T12,A,2X,10F8.3)") &
1462 0 : "Total", (SUM(q_work(:, i)), i=1, nsbsize), SUM(q_work)
1463 0 : WRITE (UNIT=output_unit, FMT="(A)") ""
1464 0 : DEALLOCATE (q_work)
1465 : IF (debug) THEN
1466 : ! Print the DFT+U occupation matrix
1467 : WRITE (UNIT=output_unit, FMT="(T9,70I10)") (i, i=1, n)
1468 : DO i = 1, n
1469 : WRITE (UNIT=output_unit, FMT="(T3,I6,70F10.6)") i, q_matrix(i, :)
1470 : END DO
1471 : ! Print the eigenvalues and eigenvectors of the occupation matrix
1472 : ALLOCATE (q_eigval(n))
1473 : q_eigval(:) = 0.0_dp
1474 : ALLOCATE (q_eigvec(n, n))
1475 : q_eigvec(:, :) = 0.0_dp
1476 : CALL jacobi(q_matrix, q_eigval, q_eigvec)
1477 : WRITE (UNIT=output_unit, FMT="(/,T9,70I10)") (i, i=1, n)
1478 : WRITE (UNIT=output_unit, FMT="(T9,71F10.6)") (q_eigval(i), i=1, n), &
1479 : SUM(q_eigval(1:n))
1480 : DO i = 1, n
1481 : WRITE (UNIT=output_unit, FMT="(T3,I6,70F10.6)") i, q_eigvec(i, :)
1482 : END DO
1483 : DEALLOCATE (q_eigval)
1484 : DEALLOCATE (q_eigvec)
1485 : END IF ! debug
1486 : END IF
1487 : IF (debug) THEN
1488 : ! Print the full atomic occupation matrix block
1489 : ALLOCATE (q_work(nsgf_kind, nsgf_kind))
1490 : q_work(:, :) = 0.0_dp
1491 : IF (ASSOCIATED(q_block)) q_work(:, :) = q_block(:, :)
1492 : CALL para_env%sum(q_work)
1493 : IF (output_unit > 0) THEN
1494 : norb = SIZE(q_work, 1)
1495 : WRITE (UNIT=output_unit, FMT="(/,T9,200I10)") (i, i=1, norb)
1496 : DO i = 1, norb
1497 : WRITE (UNIT=output_unit, FMT="(T3,I6,200F10.6)") i, q_work(i, :)
1498 : END DO
1499 : ALLOCATE (q_eigval(norb))
1500 : q_eigval(:) = 0.0_dp
1501 : ALLOCATE (q_eigvec(norb, norb))
1502 : q_eigvec(:, :) = 0.0_dp
1503 : CALL jacobi(q_work, q_eigval, q_eigvec)
1504 : WRITE (UNIT=output_unit, FMT="(/,T9,200I10)") (i, i=1, norb)
1505 : WRITE (UNIT=output_unit, FMT="(T9,201F10.6)") (q_eigval(i), i=1, norb), &
1506 : SUM(q_eigval(1:norb))
1507 : DO i = 1, norb
1508 : WRITE (UNIT=output_unit, FMT="(T3,I6,200F10.6)") i, q_eigvec(i, :)
1509 : END DO
1510 : DEALLOCATE (q_eigval)
1511 : DEALLOCATE (q_eigvec)
1512 : END IF
1513 : DEALLOCATE (q_work)
1514 : END IF ! debug
1515 : END IF ! should output
1516 :
1517 : END DO ! next atom "iatom" of atomic kind "ikind"
1518 :
1519 9600 : IF (ALLOCATED(q_matrix)) THEN
1520 1920 : DEALLOCATE (q_matrix)
1521 : END IF
1522 :
1523 : END DO ! next atomic kind "ikind"
1524 :
1525 : ! Add V(i,j)[U] to V(i,j)[DFT]
1526 :
1527 3840 : IF (ASSOCIATED(sm_h)) THEN
1528 :
1529 3216 : DO ikind = 1, nkind
1530 :
1531 2144 : IF (.NOT. is_plus_u_kind(ikind)) CYCLE
1532 :
1533 1072 : kind_a => atomic_kind_set(ikind)
1534 :
1535 : CALL get_atomic_kind(atomic_kind=kind_a, &
1536 : atom_list=atom_list, &
1537 1072 : natom=natom_of_kind)
1538 :
1539 3216 : DO iatom = 1, natom_of_kind
1540 :
1541 1072 : atom_a = atom_list(iatom)
1542 :
1543 : CALL dbcsr_get_block_p(matrix=sm_h, &
1544 : row=atom_a, &
1545 : col=atom_a, &
1546 : block=h_block, &
1547 1072 : found=found)
1548 :
1549 1072 : IF (.NOT. ASSOCIATED(h_block)) CYCLE
1550 :
1551 : CALL dbcsr_get_block_p(matrix=sm_v, &
1552 : row=atom_a, &
1553 : col=atom_a, &
1554 : block=v_block, &
1555 536 : found=found)
1556 536 : CPASSERT(ASSOCIATED(v_block))
1557 :
1558 4288 : IF (orthonormal_basis) THEN
1559 0 : DO isgf = 1, SIZE(h_block, 1)
1560 0 : h_block(isgf, isgf) = h_block(isgf, isgf) + v_block(isgf, isgf)
1561 : END DO
1562 : ELSE
1563 : CALL dbcsr_get_block_p(matrix=sm_s, &
1564 : row=atom_a, &
1565 : col=atom_a, &
1566 : block=s_block, &
1567 536 : found=found)
1568 536 : CPASSERT(ASSOCIATED(s_block))
1569 7504 : DO jsgf = 1, SIZE(h_block, 2)
1570 98624 : DO isgf = 1, SIZE(h_block, 1)
1571 97552 : h_block(isgf, jsgf) = h_block(isgf, jsgf) + v_block(isgf, jsgf)*s_block(isgf, jsgf)
1572 : END DO
1573 : END DO
1574 : END IF ! orthonormal basis set
1575 :
1576 : END DO ! next atom "iatom" of atomic kind "ikind"
1577 :
1578 : END DO ! Next atomic kind "ikind"
1579 :
1580 : END IF ! An update of the Hamiltonian matrix is requested
1581 :
1582 : END DO ! next cell image
1583 :
1584 : END DO ! next spin "ispin"
1585 :
1586 : ! Collect the energy contributions from all processes
1587 :
1588 1290 : CALL para_env%sum(energy%dft_plus_u)
1589 :
1590 1290 : IF (energy%dft_plus_u < 0.0_dp) THEN
1591 0 : IF (.NOT. occupation_enforced) THEN
1592 : CALL cp_warn(__LOCATION__, &
1593 : "DFT+U energy contribution is negative possibly due "// &
1594 0 : "to unphysical Mulliken charges!")
1595 : END IF
1596 : END IF
1597 :
1598 1290 : CALL dbcsr_deallocate_matrix(sm_q)
1599 1290 : CALL dbcsr_deallocate_matrix(sm_v)
1600 :
1601 1290 : CALL timestop(handle)
1602 :
1603 3870 : END SUBROUTINE mulliken
1604 :
1605 : ! **************************************************************************************************
1606 : !> \brief Add a DFT+U contribution to the Hamiltonian matrix\n
1607 : !> using a method based on Mulliken charges
1608 : !> \f[q_\mu = \sum_\nu \frac{1}{2}(P_{\mu\nu} S_{\nu\mu} +
1609 : !> S_{\mu\nu} P_{\nu\mu})
1610 : !> = \sum_\nu P_{\mu\nu} S_{\nu\mu}\f]
1611 : !> where \b P and \b S are the density and the
1612 : !> overlap matrix, respectively.
1613 : !> \param[in] qs_env Quickstep environment
1614 : !> \param orthonormal_basis ...
1615 : !> \param[in,out] matrix_h Hamiltonian matrices for each spin
1616 : !> \param[in,out] matrix_w Energy weighted density matrices for each spin
1617 : !> \param should_output ...
1618 : !> \param output_unit ...
1619 : !> \param print_level ...
1620 : !> \date 11.01.2008
1621 : !> \par
1622 : !> \f{eqnarray*}{
1623 : !> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}\\\
1624 : !> & = & E^{\rm DFT} + \frac{1}{2}(U - J)\sum_\mu (q_\mu - q_\mu^2)\\[1ex]
1625 : !> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
1626 : !> & = & \frac{\partial E^{\rm DFT}}
1627 : !> {\partial P_{\mu\nu}} +
1628 : !> \frac{\partial E^{\rm U}}
1629 : !> {\partial P_{\mu\nu}}\\\
1630 : !> & = & H_{\mu\nu} +
1631 : !> \frac{\partial E^{\rm U}}{\partial q_\mu}
1632 : !> \frac{\partial q_\mu}{\partial P_{\mu\nu}}\\\
1633 : !> & = & H_{\mu\nu} +
1634 : !> \frac{1}{2}(U - J)(1 - q_\mu - q_\nu) S_{\mu\nu}\\\
1635 : !> \f}
1636 : !> \author Matthias Krack (MK)
1637 : !> \version 1.0
1638 : !> \note The use of any full matrices was avoided. Thus no ScaLAPACK
1639 : !> calls are performed
1640 : ! **************************************************************************************************
1641 404 : SUBROUTINE mulliken_charges(qs_env, orthonormal_basis, matrix_h, matrix_w, &
1642 : should_output, output_unit, print_level)
1643 :
1644 : TYPE(qs_environment_type), POINTER :: qs_env
1645 : LOGICAL, INTENT(IN) :: orthonormal_basis
1646 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
1647 : POINTER :: matrix_h, matrix_w
1648 : LOGICAL, INTENT(IN) :: should_output
1649 : INTEGER, INTENT(IN) :: output_unit, print_level
1650 :
1651 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mulliken_charges'
1652 :
1653 : CHARACTER(LEN=10) :: spin_info
1654 404 : CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: symbol
1655 : CHARACTER(LEN=default_string_length) :: atomic_kind_name
1656 : INTEGER :: atom_a, handle, i, iatom, ic, ikind, isb, iset, isgf, ishell, ispin, jatom, jsgf, &
1657 : lu, m, natom, natom_of_kind, nimg, nkind, nsb, nset, nsgf, nspin, sgf
1658 404 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
1659 404 : INTEGER, DIMENSION(:), POINTER :: atom_list, nshell
1660 404 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
1661 : LOGICAL :: dft_plus_u_atom, found, just_energy
1662 : REAL(KIND=dp) :: eps_u_ramping, fspin, q, u_minus_j, &
1663 : u_minus_j_target, u_ramping, v
1664 404 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dEdq, trps
1665 404 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: q_ii
1666 404 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, s_block, w_block
1667 404 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1668 : TYPE(dbcsr_iterator_type) :: iter
1669 404 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
1670 : TYPE(dbcsr_type), POINTER :: sm_h, sm_p, sm_s, sm_w
1671 : TYPE(dft_control_type), POINTER :: dft_control
1672 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1673 : TYPE(mp_para_env_type), POINTER :: para_env
1674 404 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1675 : TYPE(qs_energy_type), POINTER :: energy
1676 404 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1677 : TYPE(qs_rho_type), POINTER :: rho
1678 :
1679 404 : CALL timeset(routineN, handle)
1680 :
1681 404 : NULLIFY (atom_list)
1682 404 : NULLIFY (atomic_kind_set)
1683 404 : NULLIFY (qs_kind_set)
1684 404 : NULLIFY (dft_control)
1685 404 : NULLIFY (energy)
1686 404 : NULLIFY (first_sgf)
1687 404 : NULLIFY (h_block)
1688 404 : NULLIFY (matrix_p)
1689 404 : NULLIFY (matrix_s)
1690 404 : NULLIFY (l)
1691 404 : NULLIFY (last_sgf)
1692 404 : NULLIFY (nshell)
1693 404 : NULLIFY (orb_basis_set)
1694 404 : NULLIFY (p_block)
1695 404 : NULLIFY (particle_set)
1696 404 : NULLIFY (rho)
1697 404 : NULLIFY (s_block)
1698 404 : NULLIFY (sm_h)
1699 404 : NULLIFY (sm_p)
1700 404 : NULLIFY (sm_s)
1701 404 : NULLIFY (w_block)
1702 404 : NULLIFY (para_env)
1703 :
1704 : CALL get_qs_env(qs_env=qs_env, &
1705 : atomic_kind_set=atomic_kind_set, &
1706 : qs_kind_set=qs_kind_set, &
1707 : dft_control=dft_control, &
1708 : energy=energy, &
1709 : particle_set=particle_set, &
1710 : rho=rho, &
1711 404 : para_env=para_env)
1712 :
1713 404 : CPASSERT(ASSOCIATED(atomic_kind_set))
1714 404 : CPASSERT(ASSOCIATED(dft_control))
1715 404 : CPASSERT(ASSOCIATED(energy))
1716 404 : CPASSERT(ASSOCIATED(particle_set))
1717 404 : CPASSERT(ASSOCIATED(rho))
1718 :
1719 404 : IF (orthonormal_basis) THEN
1720 404 : NULLIFY (sm_s)
1721 : ELSE
1722 : ! Get overlap matrix in sparse format
1723 : CALL get_qs_env(qs_env=qs_env, &
1724 404 : matrix_s_kp=matrix_s)
1725 404 : CPASSERT(ASSOCIATED(matrix_s))
1726 : END IF
1727 :
1728 : ! Get density matrices in sparse format
1729 :
1730 404 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1731 :
1732 404 : energy%dft_plus_u = 0.0_dp
1733 :
1734 404 : nspin = dft_control%nspins
1735 404 : nimg = dft_control%nimages
1736 :
1737 404 : IF (nspin == 2) THEN
1738 : fspin = 1.0_dp
1739 : ELSE
1740 248 : fspin = 0.5_dp
1741 : END IF
1742 :
1743 : ! Get the total number of atoms, contracted spherical Gaussian basis
1744 : ! functions, and atomic kinds
1745 :
1746 404 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
1747 404 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
1748 :
1749 404 : nkind = SIZE(atomic_kind_set)
1750 :
1751 1212 : ALLOCATE (first_sgf_atom(natom))
1752 404 : first_sgf_atom(:) = 0
1753 :
1754 : CALL get_particle_set(particle_set, qs_kind_set, &
1755 404 : first_sgf=first_sgf_atom)
1756 :
1757 1212 : ALLOCATE (trps(nsgf))
1758 404 : trps(:) = 0.0_dp
1759 :
1760 404 : IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
1761 990 : ALLOCATE (dEdq(nsgf))
1762 330 : just_energy = .FALSE.
1763 : ELSE
1764 : just_energy = .TRUE.
1765 : END IF
1766 :
1767 : ! Loop over all spins
1768 :
1769 964 : DO ispin = 1, nspin
1770 :
1771 560 : IF (.NOT. just_energy) dEdq(:) = 0.0_dp
1772 :
1773 : ! Calculate Trace(P*S) assuming symmetric matrices
1774 :
1775 560 : trps(:) = 0.0_dp
1776 :
1777 7000 : DO ic = 1, nimg
1778 6440 : IF (orthonormal_basis) THEN
1779 : NULLIFY (sm_s)
1780 : ELSE
1781 6440 : sm_s => matrix_s(1, ic)%matrix
1782 : END IF
1783 6440 : sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
1784 :
1785 6440 : CALL dbcsr_iterator_start(iter, sm_p)
1786 :
1787 25760 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1788 :
1789 19320 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
1790 :
1791 25760 : IF (orthonormal_basis) THEN
1792 :
1793 0 : IF (iatom /= jatom) CYCLE
1794 :
1795 0 : IF (ASSOCIATED(p_block)) THEN
1796 0 : sgf = first_sgf_atom(iatom)
1797 0 : DO isgf = 1, SIZE(p_block, 1)
1798 0 : trps(sgf) = trps(sgf) + p_block(isgf, isgf)
1799 0 : sgf = sgf + 1
1800 : END DO
1801 : END IF
1802 :
1803 : ELSE
1804 :
1805 : CALL dbcsr_get_block_p(matrix=sm_s, &
1806 : row=iatom, &
1807 : col=jatom, &
1808 : block=s_block, &
1809 19320 : found=found)
1810 19320 : CPASSERT(ASSOCIATED(s_block))
1811 :
1812 19320 : sgf = first_sgf_atom(jatom)
1813 141680 : DO jsgf = 1, SIZE(p_block, 2)
1814 1326640 : DO isgf = 1, SIZE(p_block, 1)
1815 1326640 : trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
1816 : END DO
1817 141680 : sgf = sgf + 1
1818 : END DO
1819 :
1820 19320 : IF (iatom /= jatom) THEN
1821 9660 : sgf = first_sgf_atom(iatom)
1822 109480 : DO isgf = 1, SIZE(p_block, 1)
1823 598920 : DO jsgf = 1, SIZE(p_block, 2)
1824 598920 : trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
1825 : END DO
1826 109480 : sgf = sgf + 1
1827 : END DO
1828 : END IF
1829 :
1830 : END IF ! orthonormal basis set
1831 :
1832 : END DO ! next atom "iatom"
1833 :
1834 13440 : CALL dbcsr_iterator_stop(iter)
1835 :
1836 : END DO ! cell images
1837 :
1838 560 : CALL para_env%sum(trps)
1839 :
1840 : ! q <- Trace(PS)
1841 :
1842 : ! E[DFT+U] = E[DFT] + E[U]
1843 : ! = E[DFT] + (U - J)*(q - q**2))/2
1844 :
1845 : ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U]
1846 : ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j)
1847 : ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j))
1848 :
1849 : ! Loop over all atomic kinds
1850 :
1851 1680 : DO ikind = 1, nkind
1852 :
1853 : ! Load the required atomic kind data
1854 : CALL get_atomic_kind(atomic_kind_set(ikind), &
1855 : atom_list=atom_list, &
1856 : name=atomic_kind_name, &
1857 1120 : natom=natom_of_kind)
1858 :
1859 : CALL get_qs_kind(qs_kind_set(ikind), &
1860 : dft_plus_u_atom=dft_plus_u_atom, &
1861 : l_of_dft_plus_u=lu, &
1862 : basis_set=orb_basis_set, &
1863 : u_minus_j=u_minus_j, &
1864 : u_minus_j_target=u_minus_j_target, &
1865 : u_ramping=u_ramping, &
1866 1120 : eps_u_ramping=eps_u_ramping)
1867 :
1868 : ! Check, if this atom needs a DFT+U correction
1869 :
1870 1120 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
1871 1120 : IF (.NOT. dft_plus_u_atom) CYCLE
1872 560 : IF (lu < 0) CYCLE
1873 :
1874 : ! Apply U ramping if requested
1875 :
1876 560 : IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
1877 0 : IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
1878 0 : u_minus_j = MIN(u_minus_j + u_ramping, u_minus_j_target)
1879 0 : CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
1880 : END IF
1881 0 : IF (should_output .AND. (output_unit > 0)) THEN
1882 : WRITE (UNIT=output_unit, FMT="(T3,A,3X,A,F0.3,A)") &
1883 0 : "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)), &
1884 0 : "U(eff) = ", u_minus_j*evolt, " eV"
1885 : END IF
1886 : END IF
1887 :
1888 560 : IF (u_minus_j == 0.0_dp) CYCLE
1889 :
1890 : ! Load the required Gaussian basis set data
1891 :
1892 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1893 : first_sgf=first_sgf, &
1894 : l=l, &
1895 : last_sgf=last_sgf, &
1896 : nset=nset, &
1897 560 : nshell=nshell)
1898 :
1899 : ! Count the relevant shell blocks of this atomic kind
1900 :
1901 560 : nsb = 0
1902 1680 : DO iset = 1, nset
1903 4480 : DO ishell = 1, nshell(iset)
1904 3920 : IF (l(ishell, iset) == lu) nsb = nsb + 1
1905 : END DO
1906 : END DO
1907 :
1908 2240 : ALLOCATE (q_ii(nsb, 2*lu + 1))
1909 :
1910 : ! Print headline if requested
1911 :
1912 560 : IF (should_output .AND. (print_level > low_print_level)) THEN
1913 0 : IF (output_unit > 0) THEN
1914 0 : ALLOCATE (symbol(2*lu + 1))
1915 0 : DO m = -lu, lu
1916 0 : symbol(lu + m + 1) = sgf_symbol(0, lu, m)
1917 : END DO
1918 0 : IF (nspin > 1) THEN
1919 0 : WRITE (UNIT=spin_info, FMT="(A8,I2)") " of spin", ispin
1920 : ELSE
1921 0 : spin_info = ""
1922 : END IF
1923 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
1924 0 : "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ", ikind, &
1925 0 : ": "//TRIM(atomic_kind_name), &
1926 0 : "Atom Shell ", (ADJUSTR(symbol(i)), i=1, 2*lu + 1), " Trace"
1927 0 : DEALLOCATE (symbol)
1928 : END IF
1929 : END IF
1930 :
1931 : ! Loop over all atoms of the current atomic kind
1932 :
1933 1120 : DO iatom = 1, natom_of_kind
1934 :
1935 560 : atom_a = atom_list(iatom)
1936 :
1937 560 : q_ii(:, :) = 0.0_dp
1938 :
1939 : ! Get diagonal block
1940 :
1941 : CALL dbcsr_get_block_p(matrix=sm_p, &
1942 : row=atom_a, &
1943 : col=atom_a, &
1944 : block=p_block, &
1945 560 : found=found)
1946 :
1947 : ! Calculate E(U) and dE(U)/dq
1948 :
1949 560 : IF (ASSOCIATED(p_block)) THEN
1950 :
1951 280 : sgf = first_sgf_atom(atom_a)
1952 :
1953 280 : isb = 0
1954 840 : DO iset = 1, nset
1955 2240 : DO ishell = 1, nshell(iset)
1956 1960 : IF (l(ishell, iset) == lu) THEN
1957 560 : isb = isb + 1
1958 560 : i = 0
1959 2240 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1960 1680 : q = fspin*trps(sgf)
1961 1680 : i = i + 1
1962 1680 : q_ii(isb, i) = q
1963 : energy%dft_plus_u = energy%dft_plus_u + &
1964 1680 : 0.5_dp*u_minus_j*(q - q**2)/fspin
1965 1680 : IF (.NOT. just_energy) THEN
1966 1338 : dEdq(sgf) = dEdq(sgf) + u_minus_j*(0.5_dp - q)
1967 : END IF
1968 2240 : sgf = sgf + 1
1969 : END DO ! next contracted spherical Gaussian function "isgf"
1970 : ELSE
1971 840 : sgf = sgf + last_sgf(ishell, iset) - first_sgf(ishell, iset) + 1
1972 : END IF ! angular momentum requested for DFT+U correction
1973 : END DO ! next shell "ishell"
1974 : END DO ! next shell set "iset"
1975 :
1976 : END IF ! this process is the owner of the sparse matrix block?
1977 :
1978 : ! Consider print requests
1979 :
1980 1680 : IF (should_output .AND. (print_level > low_print_level)) THEN
1981 0 : CALL para_env%sum(q_ii)
1982 0 : IF (output_unit > 0) THEN
1983 0 : DO isb = 1, nsb
1984 : WRITE (UNIT=output_unit, FMT="(T3,I6,2X,I6,2X,10F8.3)") &
1985 0 : atom_a, isb, q_ii(isb, :), SUM(q_ii(isb, :))
1986 : END DO
1987 : WRITE (UNIT=output_unit, FMT="(T12,A,2X,10F8.3)") &
1988 0 : "Total", (SUM(q_ii(:, i)), i=1, 2*lu + 1), SUM(q_ii)
1989 0 : WRITE (UNIT=output_unit, FMT="(A)") ""
1990 : END IF
1991 : END IF ! should output
1992 :
1993 : END DO ! next atom "iatom" of atomic kind "ikind"
1994 :
1995 2800 : IF (ALLOCATED(q_ii)) THEN
1996 560 : DEALLOCATE (q_ii)
1997 : END IF
1998 :
1999 : END DO ! next atomic kind "ikind"
2000 :
2001 560 : IF (.NOT. just_energy) THEN
2002 446 : CALL para_env%sum(dEdq)
2003 : END IF
2004 :
2005 : ! Add V(i,j)[U] to V(i,j)[DFT]
2006 :
2007 560 : IF (PRESENT(matrix_h)) THEN
2008 :
2009 6576 : DO ic = 1, nimg
2010 6168 : IF (orthonormal_basis) THEN
2011 : NULLIFY (sm_s)
2012 : ELSE
2013 6168 : sm_s => matrix_s(1, ic)%matrix
2014 : END IF
2015 6168 : sm_h => matrix_h(ispin, ic)%matrix
2016 :
2017 6168 : CALL dbcsr_iterator_start(iter, sm_h)
2018 :
2019 24672 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2020 :
2021 18504 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, h_block)
2022 :
2023 24672 : IF (orthonormal_basis) THEN
2024 :
2025 0 : IF (iatom /= jatom) CYCLE
2026 :
2027 0 : IF (ASSOCIATED(h_block)) THEN
2028 0 : sgf = first_sgf_atom(iatom)
2029 0 : DO isgf = 1, SIZE(h_block, 1)
2030 0 : h_block(isgf, isgf) = h_block(isgf, isgf) + dEdq(sgf)
2031 0 : sgf = sgf + 1
2032 : END DO
2033 : END IF
2034 :
2035 : ELSE
2036 :
2037 : ! Request katom just to check for consistent sparse matrix pattern
2038 :
2039 : CALL dbcsr_get_block_p(matrix=sm_s, &
2040 : row=iatom, &
2041 : col=jatom, &
2042 : block=s_block, &
2043 18504 : found=found)
2044 18504 : CPASSERT(ASSOCIATED(s_block))
2045 :
2046 : ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
2047 :
2048 18504 : sgf = first_sgf_atom(iatom)
2049 :
2050 185040 : DO isgf = 1, SIZE(h_block, 1)
2051 166536 : IF (dEdq(sgf) /= 0.0_dp) THEN
2052 55512 : v = 0.5_dp*dEdq(sgf)
2053 481104 : DO jsgf = 1, SIZE(h_block, 2)
2054 481104 : h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
2055 : END DO
2056 : END IF
2057 185040 : sgf = sgf + 1
2058 : END DO
2059 :
2060 18504 : sgf = first_sgf_atom(jatom)
2061 :
2062 135696 : DO jsgf = 1, SIZE(h_block, 2)
2063 117192 : IF (dEdq(sgf) /= 0.0_dp) THEN
2064 18504 : v = 0.5_dp*dEdq(sgf)
2065 259056 : DO isgf = 1, SIZE(h_block, 1)
2066 259056 : h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
2067 : END DO
2068 : END IF
2069 135696 : sgf = sgf + 1
2070 : END DO
2071 :
2072 : END IF ! orthonormal basis set
2073 :
2074 : END DO ! Next atom "iatom"
2075 :
2076 12744 : CALL dbcsr_iterator_stop(iter)
2077 :
2078 : END DO
2079 :
2080 : END IF ! An update of the Hamiltonian matrix is requested
2081 :
2082 : ! Calculate the contribution (non-Pulay part) to the derivatives
2083 : ! w.r.t. the nuclear positions, which requires an update of the
2084 : ! energy weighted density W.
2085 :
2086 964 : IF (PRESENT(matrix_w) .AND. (.NOT. orthonormal_basis)) THEN
2087 :
2088 196 : DO ic = 1, nimg
2089 158 : sm_s => matrix_s(1, ic)%matrix
2090 158 : sm_p => matrix_p(ispin, ic)%matrix
2091 158 : sm_w => matrix_w(ispin, ic)%matrix
2092 :
2093 158 : CALL dbcsr_iterator_start(iter, sm_p)
2094 :
2095 632 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2096 :
2097 474 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
2098 :
2099 : ! Skip the diagonal blocks of the W matrix
2100 :
2101 474 : IF (iatom == jatom) CYCLE
2102 :
2103 : ! Request katom just to check for consistent sparse matrix patterns
2104 :
2105 : CALL dbcsr_get_block_p(matrix=sm_w, &
2106 : row=iatom, &
2107 : col=jatom, &
2108 : block=w_block, &
2109 237 : found=found)
2110 237 : CPASSERT(ASSOCIATED(w_block))
2111 :
2112 : ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
2113 :
2114 237 : sgf = first_sgf_atom(iatom)
2115 :
2116 2686 : DO isgf = 1, SIZE(w_block, 1)
2117 2449 : IF (dEdq(sgf) /= 0.0_dp) THEN
2118 948 : v = -0.5_dp*dEdq(sgf)
2119 5688 : DO jsgf = 1, SIZE(w_block, 2)
2120 5688 : w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
2121 : END DO
2122 : END IF
2123 2686 : sgf = sgf + 1
2124 : END DO
2125 :
2126 237 : sgf = first_sgf_atom(jatom)
2127 :
2128 1580 : DO jsgf = 1, SIZE(w_block, 2)
2129 1185 : IF (dEdq(sgf) /= 0.0_dp) THEN
2130 0 : v = -0.5_dp*dEdq(sgf)
2131 0 : DO isgf = 1, SIZE(w_block, 1)
2132 0 : w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
2133 : END DO
2134 : END IF
2135 1659 : sgf = sgf + 1
2136 : END DO
2137 :
2138 : END DO ! next block node "jatom"
2139 :
2140 354 : CALL dbcsr_iterator_stop(iter)
2141 :
2142 : END DO
2143 :
2144 : END IF ! W matrix update requested
2145 :
2146 : END DO ! next spin "ispin"
2147 :
2148 : ! Collect the energy contributions from all processes
2149 :
2150 404 : CALL para_env%sum(energy%dft_plus_u)
2151 :
2152 404 : IF (energy%dft_plus_u < 0.0_dp) THEN
2153 : CALL cp_warn(__LOCATION__, &
2154 : "DFT+U energy contribution is negative possibly due "// &
2155 0 : "to unphysical Mulliken charges!")
2156 : END IF
2157 :
2158 : ! Release local work storage
2159 :
2160 404 : IF (ALLOCATED(first_sgf_atom)) THEN
2161 404 : DEALLOCATE (first_sgf_atom)
2162 : END IF
2163 :
2164 404 : IF (ALLOCATED(trps)) THEN
2165 404 : DEALLOCATE (trps)
2166 : END IF
2167 :
2168 404 : IF (ALLOCATED(dEdq)) THEN
2169 330 : DEALLOCATE (dEdq)
2170 : END IF
2171 :
2172 404 : CALL timestop(handle)
2173 :
2174 1212 : END SUBROUTINE mulliken_charges
2175 :
2176 : ! **************************************************************************************************
2177 : !> \brief Calculate the tensorial DFT+U+J energy contribution, Hamiltonian correction,
2178 : !> and optional force contribution.
2179 : !> \param[in] qs_env Quickstep environment
2180 : !> \param[in,out] matrix_h Hamiltonian matrices
2181 : !> \param[in,out] matrix_w Energy-weighted density matrices for each spin channel
2182 : !> \param[in] should_output Whether detailed output should be printed
2183 : !> \param[in] output_unit ...
2184 : !> \param[in] print_level ...
2185 : !> \date 14.03.2026
2186 : !> \author Ziwei Chai (ZC)
2187 : ! **************************************************************************************************
2188 558 : SUBROUTINE tensorial(qs_env, matrix_h, matrix_w, should_output, output_unit, print_level)
2189 :
2190 : TYPE(qs_environment_type), POINTER :: qs_env
2191 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
2192 : POINTER :: matrix_h, matrix_w
2193 : LOGICAL, INTENT(IN) :: should_output
2194 : INTEGER, INTENT(IN) :: output_unit, print_level
2195 :
2196 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tensorial'
2197 :
2198 : CHARACTER(LEN=10) :: spin_info
2199 558 : CHARACTER(LEN=6), ALLOCATABLE :: symbol(:)
2200 : CHARACTER(LEN=default_string_length) :: atomic_kind_name
2201 : INTEGER :: atom_a, atom_aa, handle, i, i0, iatom, ikind, iorb, isb, iset, isgf, ishell, &
2202 : ispin, j, jatom, jj, jkind, jsb, jset, jsgf, jshell, katom, lu, m, matom_of_kind, &
2203 : max_scf, n, natom, natom_of_kind, nkind, norb, nsb, nsbsize, nset, nsgf, nsgf_kind, &
2204 : nsgf_kind2, nspin, timevalues(8)
2205 558 : INTEGER, ALLOCATABLE :: first_sgf_atom(:)
2206 : INTEGER, DIMENSION(1) :: iloc
2207 558 : INTEGER, DIMENSION(:), POINTER :: atom_list, atom_list2, nshell, orbitals
2208 558 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
2209 : LOGICAL :: debug, dft_plus_u_atom, found, &
2210 : just_energy, smear
2211 558 : LOGICAL, ALLOCATABLE :: orb_occ(:)
2212 : REAL(KIND=dp) :: eps_scf, eps_u_ramping, fspin, hund_j, occ, perturbation_strength, &
2213 : some_real, trq, trq2, trqxq_hund_j, u_minus_j, u_minus_j_target, u_ramping
2214 558 : REAL(KIND=dp), ALLOCATABLE :: a(:), q_a_matrix(:, :), q_b_matrix(:, :), q_eigval(:), &
2215 558 : q_eigvec(:, :), q_work(:, :), s_inv_matrix(:, :), s_matrix(:, :), tmp_matrix1(:, :), &
2216 558 : tmp_matrix2(:, :), tmp_matrix3(:, :), v_matrix(:, :), vhxc_matrix(:, :)
2217 :
2218 : TYPE :: block3d
2219 : REAL(KIND=dp), ALLOCATABLE :: matrix(:, :, :)
2220 : END TYPE block3d
2221 : TYPE :: block4d
2222 : REAL(KIND=dp), ALLOCATABLE :: matrix(:, :, :, :)
2223 : END TYPE block4d
2224 558 : TYPE(block3d), ALLOCATABLE :: a_matrix(:, :), d_matrix(:, :), &
2225 558 : e_matrix(:, :), i_matrix(:, :), &
2226 558 : j_matrix(:, :)
2227 558 : TYPE(block4d), ALLOCATABLE :: dd_matrix(:, :), ee_matrix(:, :), &
2228 558 : ff_matrix(:, :)
2229 558 : REAL(KIND=dp), DIMENSION(:), POINTER :: ao_coef
2230 558 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a_block, d_block, e_block, f_block, &
2231 558 : i_block, j_block, q_a_block, q_b_block, &
2232 558 : s_block, tmp_block, v_block, vhxc_block
2233 :
2234 558 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2235 : TYPE(mp_para_env_type), POINTER :: para_env
2236 558 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s, matrix_vhxc
2237 558 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_s
2238 : TYPE(dbcsr_type), POINTER :: sm_a, sm_b, sm_c, sm_d, sm_e, sm_f, &
2239 : sm_h, sm_i, sm_j, sm_k, sm_q_a, sm_q_b, &
2240 : sm_s, sm_v, sm_w, sm_mid, sm_vhxc, sm_p
2241 : TYPE(dbcsr_iterator_type) :: iter
2242 : TYPE(dft_control_type), POINTER :: dft_control
2243 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
2244 558 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2245 : TYPE(qs_energy_type), POINTER :: energy
2246 558 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2247 558 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2248 : TYPE(qs_rho_type), POINTER :: rho
2249 : TYPE(qs_scf_env_type), POINTER :: scf_env
2250 :
2251 558 : CALL timeset(routineN, handle)
2252 :
2253 558 : debug = .FALSE.
2254 558 : smear = .FALSE.
2255 558 : max_scf = -1
2256 558 : eps_scf = 1.0E30_dp
2257 :
2258 558 : NULLIFY (atom_list)
2259 558 : NULLIFY (atom_list2)
2260 558 : NULLIFY (atomic_kind_set)
2261 558 : NULLIFY (qs_kind_set)
2262 558 : NULLIFY (dft_control)
2263 558 : NULLIFY (energy)
2264 558 : NULLIFY (force)
2265 558 : NULLIFY (first_sgf)
2266 558 : NULLIFY (matrix_p)
2267 558 : NULLIFY (matrix_s)
2268 558 : NULLIFY (matrix_vhxc)
2269 558 : NULLIFY (l)
2270 558 : NULLIFY (last_sgf)
2271 558 : NULLIFY (nshell)
2272 558 : NULLIFY (orb_basis_set)
2273 558 : NULLIFY (orbitals)
2274 558 : NULLIFY (ao_coef)
2275 558 : NULLIFY (particle_set)
2276 558 : NULLIFY (q_a_block)
2277 558 : NULLIFY (q_b_block)
2278 558 : NULLIFY (vhxc_block)
2279 558 : NULLIFY (rho)
2280 558 : NULLIFY (scf_env)
2281 558 : NULLIFY (sm_h)
2282 558 : NULLIFY (sm_p)
2283 558 : NULLIFY (sm_q_a)
2284 558 : NULLIFY (sm_q_b)
2285 558 : NULLIFY (sm_s)
2286 558 : NULLIFY (sm_v)
2287 558 : NULLIFY (sm_a)
2288 558 : NULLIFY (sm_b)
2289 558 : NULLIFY (sm_c)
2290 558 : NULLIFY (sm_d)
2291 558 : NULLIFY (sm_e)
2292 558 : NULLIFY (sm_f)
2293 558 : NULLIFY (sm_i)
2294 558 : NULLIFY (sm_j)
2295 558 : NULLIFY (sm_k)
2296 558 : NULLIFY (sm_mid)
2297 558 : NULLIFY (sm_vhxc)
2298 558 : NULLIFY (v_block)
2299 558 : NULLIFY (a_block)
2300 558 : NULLIFY (d_block)
2301 558 : NULLIFY (e_block)
2302 558 : NULLIFY (f_block)
2303 558 : NULLIFY (i_block)
2304 558 : NULLIFY (j_block)
2305 558 : NULLIFY (para_env)
2306 558 : NULLIFY (s_block)
2307 558 : NULLIFY (tmp_block)
2308 :
2309 : CALL get_qs_env(qs_env=qs_env, &
2310 : atomic_kind_set=atomic_kind_set, &
2311 : qs_kind_set=qs_kind_set, &
2312 : dft_control=dft_control, &
2313 : energy=energy, &
2314 : matrix_s=matrix_s, &
2315 : matrix_vhxc=matrix_vhxc, &
2316 : matrix_s_kp=matrixkp_s, &
2317 : particle_set=particle_set, &
2318 : rho=rho, &
2319 : scf_env=scf_env, &
2320 558 : para_env=para_env)
2321 :
2322 558 : CPASSERT(ASSOCIATED(atomic_kind_set))
2323 558 : CPASSERT(ASSOCIATED(dft_control))
2324 558 : CPASSERT(ASSOCIATED(energy))
2325 558 : CPASSERT(ASSOCIATED(matrix_s))
2326 558 : CPASSERT(ASSOCIATED(particle_set))
2327 558 : CPASSERT(ASSOCIATED(rho))
2328 :
2329 558 : sm_s => matrix_s(1)%matrix
2330 558 : CALL qs_rho_get(rho, rho_ao=matrix_p)
2331 :
2332 558 : energy%dft_plus_u = 0.0_dp
2333 :
2334 558 : nspin = dft_control%nspins
2335 558 : IF (dft_control%mtlr_dft_with_perturbation) THEN
2336 198 : IF (.NOT. ASSOCIATED(matrix_vhxc)) THEN
2337 0 : CPABORT("MTLR requires the projected Hxc matrix for every spin channel.")
2338 : END IF
2339 198 : IF (SIZE(matrix_vhxc) /= nspin) THEN
2340 0 : CPABORT("The number of projected Hxc matrices does not match the spin channels.")
2341 : END IF
2342 : END IF
2343 :
2344 558 : IF (nspin == 2) THEN
2345 : fspin = 1.0_dp
2346 : ELSE
2347 0 : fspin = 0.5_dp
2348 : END IF
2349 :
2350 : ! Get the total number of atoms, contracted spherical Gaussian basis
2351 : ! functions, and atomic kinds
2352 :
2353 558 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
2354 558 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
2355 :
2356 558 : nkind = SIZE(atomic_kind_set)
2357 :
2358 1674 : ALLOCATE (first_sgf_atom(natom))
2359 558 : first_sgf_atom(:) = 0
2360 :
2361 558 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf_atom)
2362 :
2363 558 : IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
2364 : just_energy = .FALSE.
2365 : ELSE
2366 116 : just_energy = .TRUE.
2367 : END IF
2368 :
2369 558 : CALL dbcsr_init_p(sm_mid)
2370 558 : CALL dbcsr_create(sm_mid, template=sm_s, matrix_type="N")
2371 :
2372 : ! Create local block diagonal matrices
2373 :
2374 558 : ALLOCATE (sm_q_a)
2375 558 : CALL dbcsr_get_block_diag(sm_s, sm_q_a)
2376 :
2377 558 : ALLOCATE (sm_q_b)
2378 558 : CALL dbcsr_get_block_diag(sm_s, sm_q_b)
2379 :
2380 558 : CALL dbcsr_set(sm_q_a, 0.0_dp)
2381 558 : CALL dbcsr_set(sm_q_b, 0.0_dp)
2382 :
2383 : ! Loop over all spins for calculating sm_q_a and sm_q_b
2384 :
2385 1674 : DO ispin = 1, nspin
2386 :
2387 1116 : sm_p => matrix_p(ispin)%matrix
2388 :
2389 1674 : IF (nspin == 2) THEN
2390 1116 : IF (ispin == 1) THEN
2391 558 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2392 558 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_a)
2393 : ELSE
2394 558 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2395 558 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_b)
2396 : END IF
2397 : ELSE
2398 0 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2399 0 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_a)
2400 0 : CALL dbcsr_copy(sm_q_b, sm_q_a)
2401 : END IF
2402 :
2403 : END DO
2404 :
2405 558 : ALLOCATE (sm_v)
2406 558 : CALL dbcsr_get_block_diag(sm_s, sm_v)
2407 :
2408 558 : IF (PRESENT(matrix_w)) THEN
2409 :
2410 2 : CALL dbcsr_init_p(sm_b)
2411 2 : CALL dbcsr_create(sm_b, template=matrix_s(1)%matrix, matrix_type="N")
2412 :
2413 2 : CALL dbcsr_init_p(sm_d)
2414 2 : CALL dbcsr_create(sm_d, template=matrix_s(1)%matrix, matrix_type="N")
2415 :
2416 2 : CALL dbcsr_init_p(sm_c)
2417 2 : CALL dbcsr_create(sm_c, template=matrix_s(1)%matrix, matrix_type="N")
2418 :
2419 2 : CALL dbcsr_init_p(sm_e)
2420 2 : CALL dbcsr_create(sm_e, template=matrix_s(1)%matrix, matrix_type="N")
2421 :
2422 2 : CALL dbcsr_init_p(sm_f)
2423 2 : CALL dbcsr_create(sm_f, template=matrix_s(1)%matrix, matrix_type="N")
2424 :
2425 2 : CALL dbcsr_init_p(sm_i)
2426 2 : CALL dbcsr_create(sm_i, template=matrix_s(1)%matrix, matrix_type="N")
2427 :
2428 2 : CALL dbcsr_init_p(sm_j)
2429 2 : CALL dbcsr_create(sm_j, template=matrix_s(1)%matrix, matrix_type="N")
2430 :
2431 2 : CALL dbcsr_init_p(sm_k)
2432 2 : CALL dbcsr_create(sm_k, template=matrix_s(1)%matrix, matrix_type="N")
2433 :
2434 : END IF
2435 :
2436 558 : ALLOCATE (sm_a)
2437 558 : CALL dbcsr_get_block_diag(sm_s, sm_a)
2438 :
2439 : ! Loop over all spins
2440 :
2441 1674 : DO ispin = 1, nspin
2442 :
2443 1116 : IF (PRESENT(matrix_h)) THEN
2444 880 : sm_h => matrix_h(ispin, 1)%matrix
2445 : ELSE
2446 : NULLIFY (sm_h)
2447 : END IF
2448 :
2449 : IF (PRESENT(matrix_w)) THEN
2450 1116 : sm_w => matrix_w(ispin, 1)%matrix
2451 : ELSE
2452 1116 : NULLIFY (sm_w)
2453 : END IF
2454 :
2455 1116 : CALL dbcsr_set(sm_v, 0.0_dp)
2456 1116 : CALL dbcsr_set(sm_a, 0.0_dp)
2457 :
2458 1116 : IF (dft_control%mtlr_dft_with_perturbation) THEN
2459 396 : IF (ispin == 1) perturbation_strength = 0.9_dp*dft_control%perturbation_strength
2460 396 : IF (ispin == 2) perturbation_strength = 1.1_dp*dft_control%perturbation_strength
2461 396 : IF (.NOT. ASSOCIATED(matrix_vhxc(ispin)%matrix)) THEN
2462 0 : CPABORT("MTLR projected Hxc matrix is not initialized.")
2463 : END IF
2464 396 : sm_vhxc => matrix_vhxc(ispin)%matrix
2465 : END IF
2466 :
2467 : ! Loop over all atomic kinds
2468 :
2469 2712 : DO ikind = 1, nkind
2470 :
2471 : CALL get_atomic_kind(atomic_kind_set(ikind), &
2472 : atom_list=atom_list, &
2473 : name=atomic_kind_name, &
2474 1596 : natom=natom_of_kind)
2475 :
2476 : CALL get_qs_kind(qs_kind_set(ikind), &
2477 : dft_plus_u_atom=dft_plus_u_atom, &
2478 : l_of_dft_plus_u=lu, &
2479 : nsgf=nsgf_kind, &
2480 : basis_set=orb_basis_set, &
2481 : u_minus_j=u_minus_j, &
2482 : hund_j=hund_j, &
2483 : ao_coef=ao_coef, &
2484 : u_minus_j_target=u_minus_j_target, &
2485 : u_ramping=u_ramping, &
2486 : eps_u_ramping=eps_u_ramping, &
2487 : orbitals=orbitals, &
2488 : eps_scf=eps_scf, &
2489 : max_scf=max_scf, &
2490 1596 : smear=smear)
2491 :
2492 1596 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
2493 1596 : IF (.NOT. dft_plus_u_atom) CYCLE
2494 1116 : IF (lu < 0) CYCLE
2495 1116 : IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) CYCLE
2496 1116 : IF (.NOT. ASSOCIATED(ao_coef)) THEN
2497 0 : CPABORT("Tensorial DFT+U requires associated ao_coef")
2498 : END IF
2499 :
2500 3348 : ALLOCATE (a(SIZE(ao_coef)))
2501 3348 : a(:) = ao_coef(:)
2502 :
2503 : ! Apply U ramping if requested
2504 :
2505 1116 : IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
2506 0 : IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
2507 0 : u_minus_j = MIN(u_minus_j + u_ramping, u_minus_j_target)
2508 0 : CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
2509 : END IF
2510 :
2511 0 : IF (should_output .AND. (output_unit > 0)) THEN
2512 : WRITE (UNIT=output_unit, FMT="(T3,A,3X,A,F0.3,A)") &
2513 0 : "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)), &
2514 0 : "U(eff) = ", u_minus_j*evolt, " eV"
2515 : END IF
2516 : END IF
2517 :
2518 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
2519 : first_sgf=first_sgf, &
2520 : l=l, &
2521 : last_sgf=last_sgf, &
2522 : nset=nset, &
2523 1116 : nshell=nshell)
2524 :
2525 1116 : nsb = 1
2526 1116 : nsbsize = (2*lu + 1)
2527 1116 : n = nsb*nsbsize
2528 :
2529 4464 : ALLOCATE (q_a_matrix(n, n))
2530 3348 : ALLOCATE (q_b_matrix(n, n))
2531 3348 : ALLOCATE (s_matrix(n, n))
2532 3348 : ALLOCATE (s_inv_matrix(n, n))
2533 3348 : ALLOCATE (v_matrix(n, n))
2534 :
2535 1116 : q_a_matrix(:, :) = 0.0_dp
2536 1116 : q_b_matrix(:, :) = 0.0_dp
2537 1116 : s_matrix(:, :) = 0.0_dp
2538 1116 : s_inv_matrix(:, :) = 0.0_dp
2539 1116 : v_matrix(:, :) = 0.0_dp
2540 :
2541 : ! Print headline if requested
2542 :
2543 1116 : IF (should_output .AND. (print_level > low_print_level)) THEN
2544 0 : IF (output_unit > 0) THEN
2545 0 : ALLOCATE (symbol(nsbsize))
2546 0 : DO m = -lu, lu
2547 0 : symbol(lu + m + 1) = sgf_symbol(0, lu, m)
2548 : END DO
2549 :
2550 0 : IF (nspin > 1) THEN
2551 0 : WRITE (UNIT=spin_info, FMT="(A8,I2)") " of spin", ispin
2552 : ELSE
2553 0 : spin_info = ""
2554 : END IF
2555 :
2556 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
2557 0 : "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ", ikind, &
2558 0 : ": "//TRIM(atomic_kind_name), &
2559 0 : "Atom Shell ", (ADJUSTR(symbol(i)), i=1, nsbsize), " Trace"
2560 :
2561 0 : DEALLOCATE (symbol)
2562 : END IF
2563 : END IF
2564 :
2565 1116 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
2566 : dft_control%mtlr_ikind == ikind) THEN
2567 396 : dft_control%vhxc(ispin) = 0.0_dp
2568 396 : dft_control%trq(ispin) = 0.0_dp
2569 1188 : ALLOCATE (vhxc_matrix(n, n))
2570 396 : vhxc_matrix(:, :) = 0.0_dp
2571 : END IF
2572 :
2573 : ! Loop over all atoms of the current atomic kind
2574 :
2575 2868 : DO iatom = 1, natom_of_kind
2576 :
2577 1752 : atom_a = atom_list(iatom)
2578 :
2579 1752 : q_a_matrix(:, :) = 0.0_dp
2580 1752 : q_b_matrix(:, :) = 0.0_dp
2581 1752 : s_matrix(:, :) = 0.0_dp
2582 1752 : s_inv_matrix(:, :) = 0.0_dp
2583 1752 : v_matrix(:, :) = 0.0_dp
2584 :
2585 : CALL dbcsr_get_block_p(matrix=sm_q_a, &
2586 : row=atom_a, &
2587 : col=atom_a, &
2588 : block=q_a_block, &
2589 1752 : found=found)
2590 :
2591 : CALL dbcsr_get_block_p(matrix=sm_q_b, &
2592 : row=atom_a, &
2593 : col=atom_a, &
2594 : block=q_b_block, &
2595 1752 : found=found)
2596 :
2597 : CALL dbcsr_get_block_p(matrix=sm_s, &
2598 : row=atom_a, &
2599 : col=atom_a, &
2600 : block=s_block, &
2601 1752 : found=found)
2602 :
2603 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
2604 1752 : dft_control%mtlr_ikind == ikind .AND. &
2605 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
2606 : CALL dbcsr_get_block_p(matrix=sm_vhxc, &
2607 : row=atom_a, &
2608 : col=atom_a, &
2609 : block=vhxc_block, &
2610 396 : found=found)
2611 : END IF
2612 :
2613 1752 : IF (ASSOCIATED(q_a_block) .AND. ASSOCIATED(q_b_block) .AND. ASSOCIATED(s_block)) THEN
2614 :
2615 876 : isb = 0
2616 1992 : DO iset = 1, nset
2617 4464 : DO ishell = 1, nshell(iset)
2618 2472 : IF (l(ishell, iset) /= lu) CYCLE
2619 1752 : isb = isb + 1
2620 1752 : i = 0
2621 5580 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
2622 2712 : i = i + 1
2623 2712 : jsb = 0
2624 9336 : DO jset = 1, nset
2625 16608 : DO jshell = 1, nshell(jset)
2626 9744 : IF (l(jshell, jset) /= lu) CYCLE
2627 5424 : jsb = jsb + 1
2628 5424 : j = 0
2629 20760 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
2630 11184 : j = j + 1
2631 11184 : q_a_matrix(i, j) = q_a_matrix(i, j) + q_a_block(isgf, jsgf)*a(isb)*a(jsb)
2632 11184 : q_b_matrix(i, j) = q_b_matrix(i, j) + q_b_block(isgf, jsgf)*a(isb)*a(jsb)
2633 11184 : s_matrix(i, j) = s_matrix(i, j) + s_block(isgf, jsgf)*a(isb)*a(jsb)
2634 : IF (ASSOCIATED(vhxc_block) .AND. &
2635 : dft_control%mtlr_dft_with_perturbation .AND. &
2636 11184 : dft_control%mtlr_ikind == ikind .AND. &
2637 9744 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
2638 792 : vhxc_matrix(i, j) = vhxc_matrix(i, j) + vhxc_block(isgf, jsgf)*a(isb)*a(jsb)
2639 : END IF
2640 : END DO
2641 : END DO
2642 : END DO
2643 : END DO
2644 : END DO
2645 : END DO
2646 :
2647 876 : CALL invert_matrix(s_matrix, s_inv_matrix, some_real)
2648 :
2649 53628 : q_a_matrix(:, :) = MATMUL(q_a_matrix, s_inv_matrix)
2650 53628 : q_b_matrix(:, :) = MATMUL(q_b_matrix, s_inv_matrix)
2651 :
2652 : IF (ASSOCIATED(vhxc_block) .AND. &
2653 : dft_control%mtlr_dft_with_perturbation .AND. &
2654 876 : dft_control%mtlr_ikind == ikind .AND. &
2655 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
2656 2574 : vhxc_matrix(:, :) = MATMUL(vhxc_matrix, s_inv_matrix)
2657 396 : DO i = 1, n
2658 198 : dft_control%vhxc(ispin) = dft_control%vhxc(ispin) + vhxc_matrix(i, i)
2659 198 : IF (ispin == 1) dft_control%trq(1) = dft_control%trq(1) + q_a_matrix(i, i)
2660 396 : IF (ispin == 2) dft_control%trq(2) = dft_control%trq(2) + q_b_matrix(i, i)
2661 : END DO
2662 : END IF
2663 :
2664 : ! Perform the requested manipulations of the (initial) orbital occupations
2665 :
2666 876 : IF (ASSOCIATED(orbitals)) THEN
2667 0 : IF ((qs_env%scf_env%iter_delta >= eps_scf) .OR. &
2668 : ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
2669 : (qs_env%scf_env%iter_count <= max_scf))) THEN
2670 :
2671 0 : ALLOCATE (orb_occ(nsbsize))
2672 0 : ALLOCATE (q_eigval(n))
2673 0 : ALLOCATE (q_eigvec(n, n))
2674 :
2675 0 : q_eigval(:) = 0.0_dp
2676 0 : q_eigvec(:, :) = 0.0_dp
2677 0 : norb = SIZE(orbitals)
2678 :
2679 0 : IF (ispin == 1) THEN
2680 0 : CALL jacobi(q_a_matrix, q_eigval, q_eigvec)
2681 0 : q_a_matrix(:, :) = 0.0_dp
2682 : ELSE
2683 0 : CALL jacobi(q_b_matrix, q_eigval, q_eigvec)
2684 0 : q_b_matrix(:, :) = 0.0_dp
2685 : END IF
2686 :
2687 0 : DO isb = 1, nsb
2688 0 : trq = 0.0_dp
2689 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
2690 0 : trq = trq + q_eigval(i)
2691 : END DO
2692 :
2693 0 : IF (smear) THEN
2694 0 : occ = trq/REAL(norb, KIND=dp)
2695 : ELSE
2696 0 : occ = 1.0_dp/fspin
2697 : END IF
2698 :
2699 0 : orb_occ(:) = .FALSE.
2700 0 : iloc = MAXLOC(q_eigvec(:, isb*nsbsize))
2701 0 : jsb = INT((iloc(1) - 1)/nsbsize) + 1
2702 0 : i = 0
2703 0 : i0 = (jsb - 1)*nsbsize + 1
2704 0 : iorb = -1000
2705 :
2706 0 : DO j = i0, jsb*nsbsize
2707 0 : i = i + 1
2708 0 : IF (i > norb) THEN
2709 0 : DO m = -lu, lu
2710 0 : IF (.NOT. orb_occ(lu + m + 1)) THEN
2711 0 : iorb = i0 + lu + m
2712 0 : orb_occ(lu + m + 1) = .TRUE.
2713 : END IF
2714 : END DO
2715 : ELSE
2716 0 : iorb = i0 + lu + orbitals(i)
2717 0 : orb_occ(lu + orbitals(i) + 1) = .TRUE.
2718 : END IF
2719 :
2720 0 : CPASSERT(iorb /= -1000)
2721 :
2722 0 : iloc = MAXLOC(q_eigvec(iorb, :))
2723 0 : q_eigval(iloc(1)) = MIN(occ, trq)
2724 :
2725 0 : IF (ispin == 1) THEN
2726 0 : q_a_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1))
2727 : ELSE
2728 0 : q_b_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1))
2729 : END IF
2730 :
2731 0 : trq = trq - q_eigval(iloc(1))
2732 : END DO
2733 : END DO
2734 :
2735 0 : IF (ispin == 1) THEN
2736 0 : q_a_matrix(:, :) = MATMUL(q_a_matrix, TRANSPOSE(q_eigvec))
2737 : ELSE
2738 0 : q_b_matrix(:, :) = MATMUL(q_b_matrix, TRANSPOSE(q_eigvec))
2739 : END IF
2740 :
2741 0 : DEALLOCATE (orb_occ)
2742 0 : DEALLOCATE (q_eigval)
2743 0 : DEALLOCATE (q_eigvec)
2744 : END IF
2745 : END IF
2746 :
2747 876 : trq = 0.0_dp
2748 876 : trq2 = 0.0_dp
2749 876 : trqxq_hund_j = 0.0_dp
2750 :
2751 876 : IF (ispin == 1) THEN
2752 1116 : DO i = 1, n
2753 678 : trq = trq + q_a_matrix(i, i)
2754 2514 : DO j = 1, n
2755 1398 : trq2 = trq2 + q_a_matrix(i, j)*q_a_matrix(j, i)
2756 1398 : trqxq_hund_j = trqxq_hund_j + q_a_matrix(i, j)*q_b_matrix(j, i)
2757 :
2758 2076 : IF (i == j) THEN
2759 678 : v_matrix(i, i) = u_minus_j*(0.5_dp - fspin*q_a_matrix(i, i))
2760 678 : v_matrix(i, i) = v_matrix(i, i) + hund_j*fspin*q_b_matrix(i, i)
2761 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
2762 678 : dft_control%mtlr_ikind == ikind .AND. &
2763 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
2764 99 : v_matrix(i, i) = v_matrix(i, i) + perturbation_strength
2765 : END IF
2766 : ELSE
2767 720 : v_matrix(i, j) = -u_minus_j*fspin*q_a_matrix(j, i)
2768 720 : v_matrix(i, j) = v_matrix(i, j) + hund_j*fspin*q_b_matrix(j, i)
2769 : END IF
2770 : END DO
2771 : END DO
2772 : ELSE
2773 1116 : DO i = 1, n
2774 678 : trq = trq + q_b_matrix(i, i)
2775 2514 : DO j = 1, n
2776 1398 : trq2 = trq2 + q_b_matrix(i, j)*q_b_matrix(j, i)
2777 1398 : trqxq_hund_j = trqxq_hund_j + q_b_matrix(i, j)*q_a_matrix(j, i)
2778 :
2779 2076 : IF (i == j) THEN
2780 678 : v_matrix(i, i) = u_minus_j*(0.5_dp - fspin*q_b_matrix(i, i))
2781 678 : v_matrix(i, i) = v_matrix(i, i) + hund_j*fspin*q_a_matrix(i, i)
2782 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
2783 678 : dft_control%mtlr_ikind == ikind .AND. &
2784 : qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
2785 99 : v_matrix(i, i) = v_matrix(i, i) + perturbation_strength
2786 : END IF
2787 : ELSE
2788 720 : v_matrix(i, j) = -u_minus_j*fspin*q_b_matrix(j, i)
2789 720 : v_matrix(i, j) = v_matrix(i, j) + hund_j*fspin*q_a_matrix(j, i)
2790 : END IF
2791 : END DO
2792 : END DO
2793 : END IF
2794 :
2795 37644 : v_matrix(:, :) = MATMUL(s_inv_matrix, v_matrix)
2796 :
2797 876 : trq = fspin*trq
2798 876 : trq2 = fspin*fspin*trq2
2799 876 : trqxq_hund_j = fspin*fspin*trqxq_hund_j
2800 :
2801 876 : energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
2802 876 : energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*hund_j*trqxq_hund_j/fspin
2803 :
2804 : ! Calculate potential V(U) = dE(U)/dq
2805 :
2806 1752 : IF (.NOT. just_energy) THEN
2807 :
2808 : CALL dbcsr_get_block_p(matrix=sm_v, &
2809 : row=atom_a, &
2810 : col=atom_a, &
2811 : block=v_block, &
2812 760 : found=found)
2813 760 : CPASSERT(ASSOCIATED(v_block))
2814 :
2815 : CALL dbcsr_get_block_p(matrix=sm_a, &
2816 : row=atom_a, &
2817 : col=atom_a, &
2818 : block=a_block, &
2819 760 : found=found)
2820 760 : CPASSERT(ASSOCIATED(a_block))
2821 :
2822 760 : isb = 0
2823 1644 : DO iset = 1, nset
2824 3536 : DO ishell = 1, nshell(iset)
2825 1892 : IF (l(ishell, iset) /= lu) CYCLE
2826 1520 : isb = isb + 1
2827 1520 : i = 0
2828 4420 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
2829 2016 : i = i + 1
2830 2016 : jsb = 0
2831 6668 : DO jset = 1, nset
2832 11040 : DO jshell = 1, nshell(jset)
2833 6264 : IF (l(jshell, jset) /= lu) CYCLE
2834 4032 : jsb = jsb + 1
2835 4032 : j = 0
2836 13800 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
2837 7008 : j = j + 1
2838 7008 : v_block(isgf, jsgf) = v_matrix(i, j)*a(isb)*a(jsb)
2839 13272 : a_block(isgf, jsgf) = s_inv_matrix(i, j)*a(isb)*a(jsb)
2840 : END DO
2841 : END DO
2842 : END DO
2843 : END DO
2844 : END DO
2845 : END DO
2846 :
2847 : END IF
2848 :
2849 : END IF
2850 :
2851 : ! Consider print requests
2852 :
2853 4620 : IF (should_output .AND. (print_level > low_print_level)) THEN
2854 0 : IF (ispin == 1) THEN
2855 0 : CALL para_env%sum(q_a_matrix)
2856 : ELSE
2857 0 : CALL para_env%sum(q_b_matrix)
2858 : END IF
2859 :
2860 0 : IF (output_unit > 0) THEN
2861 0 : ALLOCATE (q_work(nsb, nsbsize))
2862 0 : q_work(:, :) = 0.0_dp
2863 :
2864 0 : DO isb = 1, nsb
2865 0 : j = 0
2866 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
2867 0 : j = j + 1
2868 0 : IF (ispin == 1) THEN
2869 0 : q_work(isb, j) = q_a_matrix(i, i)
2870 : ELSE
2871 0 : q_work(isb, j) = q_b_matrix(i, i)
2872 : END IF
2873 : END DO
2874 : END DO
2875 :
2876 0 : DO isb = 1, nsb
2877 : WRITE (UNIT=output_unit, FMT="(T3,I6,2X,I6,2X,10F8.3)") &
2878 0 : atom_a, isb, q_work(isb, :), SUM(q_work(isb, :))
2879 : END DO
2880 :
2881 : WRITE (UNIT=output_unit, FMT="(T12,A,2X,10F8.3)") &
2882 0 : "Total", (SUM(q_work(:, i)), i=1, nsbsize), SUM(q_work)
2883 0 : WRITE (UNIT=output_unit, FMT="(A)") ""
2884 :
2885 0 : DEALLOCATE (q_work)
2886 :
2887 : IF (debug) THEN
2888 : WRITE (UNIT=output_unit, FMT="(T9,70I10)") (i, i=1, n)
2889 : DO i = 1, n
2890 : IF (ispin == 1) THEN
2891 : WRITE (UNIT=output_unit, FMT="(T3,I6,70F10.6)") i, q_a_matrix(i, :)
2892 : ELSE
2893 : WRITE (UNIT=output_unit, FMT="(T3,I6,70F10.6)") i, q_b_matrix(i, :)
2894 : END IF
2895 : END DO
2896 :
2897 : ALLOCATE (q_eigval(n))
2898 : ALLOCATE (q_eigvec(n, n))
2899 : q_eigval(:) = 0.0_dp
2900 : q_eigvec(:, :) = 0.0_dp
2901 :
2902 : IF (ispin == 1) THEN
2903 : CALL jacobi(q_a_matrix, q_eigval, q_eigvec)
2904 : ELSE
2905 : CALL jacobi(q_b_matrix, q_eigval, q_eigvec)
2906 : END IF
2907 :
2908 : WRITE (UNIT=output_unit, FMT="(/,T9,70I10)") (i, i=1, n)
2909 : WRITE (UNIT=output_unit, FMT="(T9,71F10.6)") (q_eigval(i), i=1, n), SUM(q_eigval(1:n))
2910 :
2911 : DO i = 1, n
2912 : WRITE (UNIT=output_unit, FMT="(T3,I6,70F10.6)") i, q_eigvec(i, :)
2913 : END DO
2914 :
2915 : DEALLOCATE (q_eigval)
2916 : DEALLOCATE (q_eigvec)
2917 : END IF
2918 : END IF
2919 :
2920 : IF (debug) THEN
2921 : ALLOCATE (q_work(nsgf_kind, nsgf_kind))
2922 : q_work(:, :) = 0.0_dp
2923 :
2924 : IF (ispin == 1) THEN
2925 : IF (ASSOCIATED(q_a_block)) q_work(:, :) = q_a_block(:, :)
2926 : ELSE
2927 : IF (ASSOCIATED(q_b_block)) q_work(:, :) = q_b_block(:, :)
2928 : END IF
2929 :
2930 : CALL para_env%sum(q_work)
2931 :
2932 : IF (output_unit > 0) THEN
2933 : norb = SIZE(q_work, 1)
2934 :
2935 : WRITE (UNIT=output_unit, FMT="(/,T9,200I10)") (i, i=1, norb)
2936 : DO i = 1, norb
2937 : WRITE (UNIT=output_unit, FMT="(T3,I6,200F10.6)") i, q_work(i, :)
2938 : END DO
2939 :
2940 : ALLOCATE (q_eigval(norb))
2941 : ALLOCATE (q_eigvec(norb, norb))
2942 : q_eigval(:) = 0.0_dp
2943 : q_eigvec(:, :) = 0.0_dp
2944 :
2945 : CALL jacobi(q_work, q_eigval, q_eigvec)
2946 :
2947 : WRITE (UNIT=output_unit, FMT="(/,T9,200I10)") (i, i=1, norb)
2948 : WRITE (UNIT=output_unit, FMT="(T9,201F10.6)") (q_eigval(i), i=1, norb), SUM(q_eigval(1:norb))
2949 :
2950 : DO i = 1, norb
2951 : WRITE (UNIT=output_unit, FMT="(T3,I6,200F10.6)") i, q_eigvec(i, :)
2952 : END DO
2953 :
2954 : DEALLOCATE (q_eigval)
2955 : DEALLOCATE (q_eigvec)
2956 : END IF
2957 :
2958 : DEALLOCATE (q_work)
2959 : END IF
2960 : END IF
2961 :
2962 : END DO
2963 :
2964 1116 : IF (ALLOCATED(q_a_matrix)) DEALLOCATE (q_a_matrix)
2965 1116 : IF (ALLOCATED(q_b_matrix)) DEALLOCATE (q_b_matrix)
2966 1116 : IF (ALLOCATED(s_matrix)) DEALLOCATE (s_matrix)
2967 1116 : IF (ALLOCATED(s_inv_matrix)) DEALLOCATE (s_inv_matrix)
2968 1116 : IF (ALLOCATED(v_matrix)) DEALLOCATE (v_matrix)
2969 :
2970 1116 : DEALLOCATE (a)
2971 :
2972 1116 : IF (dft_control%mtlr_dft_with_perturbation .AND. &
2973 3828 : dft_control%mtlr_ikind == ikind) THEN
2974 396 : dft_control%vhxc(ispin) = dft_control%vhxc(ispin)/REAL(n, dp)
2975 396 : CALL para_env%sum(dft_control%vhxc(ispin))
2976 396 : CALL para_env%sum(dft_control%trq(ispin))
2977 396 : DEALLOCATE (vhxc_matrix)
2978 : END IF
2979 :
2980 : END DO
2981 :
2982 : ! Add V(i,j)[U] to V(i,j)[DFT]
2983 :
2984 1674 : IF (ASSOCIATED(sm_h)) THEN
2985 880 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_v, 0.0_dp, sm_mid)
2986 880 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 1.0_dp, sm_h)
2987 : END IF
2988 :
2989 : END DO
2990 :
2991 558 : CALL dbcsr_deallocate_matrix(sm_mid)
2992 :
2993 558 : IF (PRESENT(matrix_w)) THEN
2994 :
2995 2 : CALL get_qs_env(qs_env=qs_env, force=force)
2996 6 : DO ikind = 1, nkind
2997 30 : force(ikind)%tensorial_u(:, :) = 0.0_dp
2998 : END DO
2999 2 : IF (nspin == 2) THEN
3000 2 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_b)
3001 2 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, matrix_p(2)%matrix, 0.0_dp, sm_c)
3002 : ELSE
3003 0 : CALL dbcsr_multiply("N", "N", 0.5_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_b)
3004 0 : CALL dbcsr_multiply("N", "N", 0.5_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_c)
3005 : END IF
3006 2 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_b, sm_s, 0.0_dp, sm_i)
3007 2 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_c, sm_s, 0.0_dp, sm_j)
3008 :
3009 8 : DO jj = 1, 3
3010 6 : CALL dbcsr_desymmetrize(matrixkp_s(jj + 1, 1)%matrix, sm_k)
3011 :
3012 6 : CALL dbcsr_set(sm_f, 0.0_dp)
3013 :
3014 6 : CALL dbcsr_iterator_start(iter, sm_k)
3015 33 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3016 27 : NULLIFY (tmp_block)
3017 27 : CALL dbcsr_iterator_next_block(iter, jatom, katom, tmp_block)
3018 27 : IF (jatom == katom) CYCLE
3019 27 : CALL dbcsr_put_block(matrix=sm_f, row=jatom, col=katom, block=tmp_block)
3020 : END DO
3021 6 : CALL dbcsr_iterator_stop(iter)
3022 6 : CALL dbcsr_finalize(sm_f)
3023 :
3024 6 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_b, sm_f, 0.0_dp, sm_d)
3025 6 : CALL dbcsr_multiply("N", "N", 1.0_dp, sm_c, sm_f, 0.0_dp, sm_e)
3026 :
3027 60 : ALLOCATE (a_matrix(nkind, nkind))
3028 54 : ALLOCATE (d_matrix(nkind, nkind))
3029 54 : ALLOCATE (e_matrix(nkind, nkind))
3030 54 : ALLOCATE (i_matrix(nkind, nkind))
3031 54 : ALLOCATE (j_matrix(nkind, nkind))
3032 60 : ALLOCATE (dd_matrix(nkind, nkind))
3033 54 : ALLOCATE (ee_matrix(nkind, nkind))
3034 54 : ALLOCATE (ff_matrix(nkind, nkind))
3035 :
3036 18 : DO ikind = 1, nkind
3037 :
3038 : CALL get_atomic_kind(atomic_kind_set(ikind), &
3039 : atom_list=atom_list, &
3040 12 : natom=natom_of_kind)
3041 :
3042 : CALL get_qs_kind(qs_kind_set(ikind), &
3043 : u_minus_j=u_minus_j, &
3044 : l_of_dft_plus_u=lu, &
3045 : basis_set=orb_basis_set, &
3046 : hund_j=hund_j, &
3047 : nsgf=nsgf_kind, &
3048 12 : dft_plus_u_atom=dft_plus_u_atom)
3049 :
3050 12 : IF (.NOT. dft_plus_u_atom) CYCLE
3051 6 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
3052 6 : IF (lu < 0) CYCLE
3053 6 : IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) CYCLE
3054 :
3055 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
3056 : first_sgf=first_sgf, &
3057 : l=l, &
3058 : last_sgf=last_sgf, &
3059 : nset=nset, &
3060 6 : nshell=nshell)
3061 :
3062 6 : nsb = 0
3063 18 : DO iset = 1, nset
3064 48 : DO ishell = 1, nshell(iset)
3065 42 : IF (l(ishell, iset) == lu) nsb = nsb + 1
3066 : END DO
3067 : END DO
3068 :
3069 6 : nsbsize = (2*lu + 1)
3070 6 : n = nsb*nsbsize
3071 :
3072 30 : ALLOCATE (a_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3073 24 : ALLOCATE (d_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3074 24 : ALLOCATE (e_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3075 24 : ALLOCATE (i_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3076 24 : ALLOCATE (j_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3077 :
3078 264 : a_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3079 264 : d_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3080 264 : e_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3081 264 : i_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3082 264 : j_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3083 :
3084 36 : DO iatom = 1, natom_of_kind
3085 :
3086 6 : atom_aa = atom_list(iatom)
3087 :
3088 6 : NULLIFY (a_block)
3089 : CALL dbcsr_get_block_p(matrix=sm_a, &
3090 : row=atom_aa, &
3091 : col=atom_aa, &
3092 : block=a_block, &
3093 6 : found=found)
3094 :
3095 6 : NULLIFY (d_block)
3096 : CALL dbcsr_get_block_p(matrix=sm_d, &
3097 : row=atom_aa, &
3098 : col=atom_aa, &
3099 : block=d_block, &
3100 6 : found=found)
3101 :
3102 6 : NULLIFY (e_block)
3103 : CALL dbcsr_get_block_p(matrix=sm_e, &
3104 : row=atom_aa, &
3105 : col=atom_aa, &
3106 : block=e_block, &
3107 6 : found=found)
3108 :
3109 6 : NULLIFY (i_block)
3110 : CALL dbcsr_get_block_p(matrix=sm_i, &
3111 : row=atom_aa, &
3112 : col=atom_aa, &
3113 : block=i_block, &
3114 6 : found=found)
3115 :
3116 6 : NULLIFY (j_block)
3117 : CALL dbcsr_get_block_p(matrix=sm_j, &
3118 : row=atom_aa, &
3119 : col=atom_aa, &
3120 : block=j_block, &
3121 6 : found=found)
3122 :
3123 : IF (ASSOCIATED(a_block) .AND. ASSOCIATED(d_block) .AND. ASSOCIATED(e_block) &
3124 6 : .AND. ASSOCIATED(i_block) .AND. ASSOCIATED(j_block)) THEN
3125 :
3126 3 : i = 0
3127 9 : DO iset = 1, nset
3128 24 : DO ishell = 1, nshell(iset)
3129 15 : IF (l(ishell, iset) /= lu) CYCLE
3130 30 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3131 18 : i = i + 1
3132 18 : j = 0
3133 69 : DO jset = 1, nset
3134 144 : DO jshell = 1, nshell(jset)
3135 90 : IF (l(jshell, jset) /= lu) CYCLE
3136 180 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
3137 108 : j = j + 1
3138 108 : a_matrix(ikind, ikind)%matrix(i, j, iatom) = a_block(isgf, jsgf)
3139 108 : d_matrix(ikind, ikind)%matrix(i, j, iatom) = -d_block(isgf, jsgf)
3140 108 : e_matrix(ikind, ikind)%matrix(i, j, iatom) = -e_block(isgf, jsgf)
3141 108 : i_matrix(ikind, ikind)%matrix(i, j, iatom) = i_block(isgf, jsgf)
3142 198 : j_matrix(ikind, ikind)%matrix(i, j, iatom) = j_block(isgf, jsgf)
3143 : END DO
3144 : END DO
3145 : END DO
3146 : END DO
3147 : END DO
3148 : END DO
3149 :
3150 : END IF
3151 :
3152 36 : DO jkind = 1, nkind
3153 :
3154 : CALL get_atomic_kind(atomic_kind_set(jkind), &
3155 : atom_list=atom_list2, &
3156 12 : natom=matom_of_kind)
3157 :
3158 12 : CALL get_qs_kind(qs_kind_set(jkind), nsgf=nsgf_kind2)
3159 :
3160 12 : IF (.NOT. ALLOCATED(dd_matrix(ikind, jkind)%matrix)) THEN
3161 72 : ALLOCATE (dd_matrix(ikind, jkind)%matrix(n, nsgf_kind2, natom_of_kind, matom_of_kind))
3162 1014 : dd_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3163 : END IF
3164 :
3165 12 : IF (.NOT. ALLOCATED(ee_matrix(ikind, jkind)%matrix)) THEN
3166 72 : ALLOCATE (ee_matrix(ikind, jkind)%matrix(n, nsgf_kind2, natom_of_kind, matom_of_kind))
3167 1014 : ee_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3168 : END IF
3169 :
3170 12 : IF (.NOT. ALLOCATED(ff_matrix(ikind, jkind)%matrix)) THEN
3171 72 : ALLOCATE (ff_matrix(ikind, jkind)%matrix(nsgf_kind2, n, natom_of_kind, matom_of_kind))
3172 984 : ff_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3173 : END IF
3174 :
3175 48 : DO jatom = 1, matom_of_kind
3176 :
3177 18 : atom_a = atom_list2(jatom)
3178 :
3179 18 : NULLIFY (d_block)
3180 : CALL dbcsr_get_block_p(matrix=sm_b, &
3181 : row=atom_aa, &
3182 : col=atom_a, &
3183 : block=d_block, &
3184 18 : found=found)
3185 :
3186 18 : NULLIFY (e_block)
3187 : CALL dbcsr_get_block_p(matrix=sm_c, &
3188 : row=atom_aa, &
3189 : col=atom_a, &
3190 : block=e_block, &
3191 18 : found=found)
3192 :
3193 18 : IF (ASSOCIATED(d_block) .AND. ASSOCIATED(e_block)) THEN
3194 9 : i = 0
3195 27 : DO iset = 1, nset
3196 72 : DO ishell = 1, nshell(iset)
3197 45 : IF (l(ishell, iset) /= lu) CYCLE
3198 90 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3199 54 : i = i + 1
3200 468 : dd_matrix(ikind, jkind)%matrix(i, :, iatom, jatom) = d_block(isgf, :)
3201 513 : ee_matrix(ikind, jkind)%matrix(i, :, iatom, jatom) = e_block(isgf, :)
3202 : END DO
3203 : END DO
3204 : END DO
3205 : END IF
3206 :
3207 18 : NULLIFY (f_block)
3208 : CALL dbcsr_get_block_p(matrix=sm_f, &
3209 : row=atom_a, &
3210 : col=atom_aa, &
3211 : block=f_block, &
3212 18 : found=found)
3213 :
3214 48 : IF (ASSOCIATED(f_block)) THEN
3215 6 : i = 0
3216 18 : DO iset = 1, nset
3217 48 : DO ishell = 1, nshell(iset)
3218 30 : IF (l(ishell, iset) /= lu) CYCLE
3219 60 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3220 36 : i = i + 1
3221 246 : ff_matrix(ikind, jkind)%matrix(:, i, iatom, jatom) = f_block(:, isgf)
3222 : END DO
3223 : END DO
3224 : END DO
3225 : END IF
3226 :
3227 : END DO
3228 : END DO
3229 : END DO
3230 : END DO
3231 :
3232 18 : DO ikind = 1, nkind
3233 :
3234 : CALL get_atomic_kind(atomic_kind_set(ikind), &
3235 : atom_list=atom_list, &
3236 12 : natom=natom_of_kind)
3237 :
3238 : CALL get_qs_kind(qs_kind_set(ikind), &
3239 : u_minus_j=u_minus_j, &
3240 : l_of_dft_plus_u=lu, &
3241 : basis_set=orb_basis_set, &
3242 : hund_j=hund_j, &
3243 : nsgf=nsgf_kind, &
3244 12 : dft_plus_u_atom=dft_plus_u_atom)
3245 :
3246 12 : IF (.NOT. dft_plus_u_atom) CYCLE
3247 6 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
3248 6 : IF (lu < 0) CYCLE
3249 6 : IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) CYCLE
3250 :
3251 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
3252 : first_sgf=first_sgf, &
3253 : l=l, &
3254 : last_sgf=last_sgf, &
3255 : nset=nset, &
3256 6 : nshell=nshell)
3257 :
3258 6 : CALL para_env%sum(a_matrix(ikind, ikind)%matrix(:, :, :))
3259 6 : CALL para_env%sum(d_matrix(ikind, ikind)%matrix(:, :, :))
3260 6 : CALL para_env%sum(e_matrix(ikind, ikind)%matrix(:, :, :))
3261 6 : CALL para_env%sum(i_matrix(ikind, ikind)%matrix(:, :, :))
3262 6 : CALL para_env%sum(j_matrix(ikind, ikind)%matrix(:, :, :))
3263 :
3264 42 : DO jkind = 1, nkind
3265 12 : CALL para_env%sum(dd_matrix(ikind, jkind)%matrix(:, :, :, :))
3266 12 : CALL para_env%sum(ee_matrix(ikind, jkind)%matrix(:, :, :, :))
3267 24 : CALL para_env%sum(ff_matrix(ikind, jkind)%matrix(:, :, :, :))
3268 : END DO
3269 :
3270 : END DO
3271 :
3272 : ! Force calculation
3273 :
3274 18 : DO ikind = 1, nkind
3275 :
3276 : CALL get_atomic_kind(atomic_kind_set(ikind), &
3277 : atom_list=atom_list, &
3278 12 : natom=natom_of_kind)
3279 :
3280 : CALL get_qs_kind(qs_kind_set(ikind), &
3281 : u_minus_j=u_minus_j, &
3282 : l_of_dft_plus_u=lu, &
3283 : basis_set=orb_basis_set, &
3284 : hund_j=hund_j, &
3285 : nsgf=nsgf_kind, &
3286 12 : dft_plus_u_atom=dft_plus_u_atom)
3287 :
3288 12 : IF (.NOT. dft_plus_u_atom) CYCLE
3289 6 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
3290 6 : IF (lu < 0) CYCLE
3291 6 : IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) CYCLE
3292 :
3293 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
3294 : first_sgf=first_sgf, &
3295 : l=l, &
3296 : last_sgf=last_sgf, &
3297 : nset=nset, &
3298 6 : nshell=nshell)
3299 :
3300 6 : nsb = 0
3301 18 : DO iset = 1, nset
3302 48 : DO ishell = 1, nshell(iset)
3303 42 : IF (l(ishell, iset) == lu) nsb = nsb + 1
3304 : END DO
3305 : END DO
3306 :
3307 6 : nsbsize = (2*lu + 1)
3308 6 : n = nsb*nsbsize
3309 :
3310 24 : ALLOCATE (tmp_matrix1(n, n))
3311 18 : ALLOCATE (tmp_matrix2(n, n))
3312 18 : ALLOCATE (tmp_matrix3(n, n))
3313 :
3314 12 : DO iatom = 1, natom_of_kind
3315 :
3316 6 : atom_aa = atom_list(iatom)
3317 :
3318 6 : tmp_matrix1(:, :) = MATMUL(d_matrix(ikind, ikind)%matrix(:, :, iatom), &
3319 8034 : a_matrix(ikind, ikind)%matrix(:, :, iatom))
3320 6 : tmp_matrix2(:, :) = MATMUL(e_matrix(ikind, ikind)%matrix(:, :, iatom), &
3321 8034 : a_matrix(ikind, ikind)%matrix(:, :, iatom))
3322 :
3323 6 : i_matrix(ikind, ikind)%matrix(:, :, iatom) = MATMUL(i_matrix(ikind, ikind)%matrix(:, :, iatom), &
3324 8544 : a_matrix(ikind, ikind)%matrix(:, :, iatom))
3325 6 : j_matrix(ikind, ikind)%matrix(:, :, iatom) = MATMUL(j_matrix(ikind, ikind)%matrix(:, :, iatom), &
3326 8544 : a_matrix(ikind, ikind)%matrix(:, :, iatom))
3327 :
3328 : tmp_matrix3(:, :) = 0.5_dp*u_minus_j*tmp_matrix1 - &
3329 6 : u_minus_j*MATMUL(tmp_matrix1, i_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3330 11154 : + hund_j*MATMUL(tmp_matrix2, i_matrix(ikind, ikind)%matrix(:, :, iatom))
3331 : tmp_matrix3(:, :) = tmp_matrix3 + 0.5_dp*u_minus_j*tmp_matrix2 - &
3332 24 : u_minus_j*MATMUL(tmp_matrix2, j_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3333 11154 : + hund_j*MATMUL(tmp_matrix1, j_matrix(ikind, ikind)%matrix(:, :, iatom))
3334 :
3335 42 : DO i = 1, n
3336 42 : force(ikind)%tensorial_u(jj, iatom) = force(ikind)%tensorial_u(jj, iatom) + tmp_matrix3(i, i)
3337 : END DO
3338 :
3339 24 : DO jkind = 1, nkind
3340 :
3341 : CALL get_atomic_kind(atomic_kind_set(jkind), &
3342 : atom_list=atom_list2, &
3343 12 : natom=matom_of_kind)
3344 :
3345 12 : CALL get_qs_kind(qs_kind_set(jkind), nsgf=nsgf_kind2)
3346 :
3347 36 : DO jatom = 1, matom_of_kind
3348 :
3349 18 : atom_a = atom_list2(jatom)
3350 :
3351 18 : tmp_matrix1(:, :) = MATMUL(dd_matrix(ikind, jkind)%matrix(:, :, iatom, jatom), &
3352 30762 : ff_matrix(ikind, jkind)%matrix(:, :, iatom, jatom))
3353 24894 : tmp_matrix1(:, :) = MATMUL(tmp_matrix1, a_matrix(ikind, ikind)%matrix(:, :, iatom))
3354 :
3355 18 : tmp_matrix2(:, :) = MATMUL(ee_matrix(ikind, jkind)%matrix(:, :, iatom, jatom), &
3356 30762 : ff_matrix(ikind, jkind)%matrix(:, :, iatom, jatom))
3357 24894 : tmp_matrix2(:, :) = MATMUL(tmp_matrix2, a_matrix(ikind, ikind)%matrix(:, :, iatom))
3358 :
3359 : tmp_matrix3(:, :) = 0.5_dp*u_minus_j*tmp_matrix1 - &
3360 18 : u_minus_j*MATMUL(tmp_matrix1, i_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3361 33462 : + hund_j*MATMUL(tmp_matrix2, i_matrix(ikind, ikind)%matrix(:, :, iatom))
3362 : tmp_matrix3(:, :) = tmp_matrix3 + 0.5_dp*u_minus_j*tmp_matrix2 - &
3363 72 : u_minus_j*MATMUL(tmp_matrix2, j_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3364 33462 : + hund_j*MATMUL(tmp_matrix1, j_matrix(ikind, ikind)%matrix(:, :, iatom))
3365 :
3366 30 : IF (atom_a /= atom_aa) THEN
3367 84 : DO i = 1, n
3368 : force(jkind)%tensorial_u(jj, jatom) = force(jkind)%tensorial_u(jj, jatom) &
3369 84 : + tmp_matrix3(i, i)
3370 : END DO
3371 : END IF
3372 :
3373 : END DO
3374 : END DO
3375 : END DO
3376 :
3377 6 : DEALLOCATE (tmp_matrix1)
3378 6 : DEALLOCATE (tmp_matrix2)
3379 36 : DEALLOCATE (tmp_matrix3)
3380 :
3381 : END DO
3382 :
3383 18 : DO ikind = 1, nkind
3384 :
3385 12 : IF (ALLOCATED(a_matrix(ikind, ikind)%matrix)) DEALLOCATE (a_matrix(ikind, ikind)%matrix)
3386 12 : IF (ALLOCATED(d_matrix(ikind, ikind)%matrix)) DEALLOCATE (d_matrix(ikind, ikind)%matrix)
3387 12 : IF (ALLOCATED(e_matrix(ikind, ikind)%matrix)) DEALLOCATE (e_matrix(ikind, ikind)%matrix)
3388 12 : IF (ALLOCATED(i_matrix(ikind, ikind)%matrix)) DEALLOCATE (i_matrix(ikind, ikind)%matrix)
3389 12 : IF (ALLOCATED(j_matrix(ikind, ikind)%matrix)) DEALLOCATE (j_matrix(ikind, ikind)%matrix)
3390 :
3391 42 : DO jkind = 1, nkind
3392 24 : IF (ALLOCATED(dd_matrix(ikind, jkind)%matrix)) DEALLOCATE (dd_matrix(ikind, jkind)%matrix)
3393 24 : IF (ALLOCATED(ee_matrix(ikind, jkind)%matrix)) DEALLOCATE (ee_matrix(ikind, jkind)%matrix)
3394 36 : IF (ALLOCATED(ff_matrix(ikind, jkind)%matrix)) DEALLOCATE (ff_matrix(ikind, jkind)%matrix)
3395 : END DO
3396 :
3397 : END DO
3398 :
3399 30 : DEALLOCATE (a_matrix)
3400 30 : DEALLOCATE (d_matrix)
3401 30 : DEALLOCATE (e_matrix)
3402 30 : DEALLOCATE (i_matrix)
3403 30 : DEALLOCATE (j_matrix)
3404 30 : DEALLOCATE (dd_matrix)
3405 30 : DEALLOCATE (ee_matrix)
3406 38 : DEALLOCATE (ff_matrix)
3407 :
3408 : END DO
3409 :
3410 6 : DO ikind = 1, nkind
3411 30 : force(ikind)%tensorial_u(:, :) = force(ikind)%tensorial_u(:, :)*2.0_dp
3412 : END DO
3413 :
3414 2 : CALL DATE_AND_TIME(values=timevalues)
3415 :
3416 : END IF
3417 :
3418 : ! Collect the energy contributions from all processes
3419 :
3420 558 : CALL para_env%sum(energy%dft_plus_u)
3421 :
3422 558 : IF (energy%dft_plus_u < 0.0_dp) THEN
3423 : CALL cp_warn(__LOCATION__, &
3424 : "DFT+U energy contibution is negative possibly due "// &
3425 : "to unphysical Lowdin charges. Check your input, "// &
3426 0 : "if this warning persists or try a different method!")
3427 : END IF
3428 :
3429 : ! Release local sparse matrices
3430 :
3431 558 : CALL dbcsr_deallocate_matrix(sm_q_a)
3432 558 : CALL dbcsr_deallocate_matrix(sm_q_b)
3433 558 : CALL dbcsr_deallocate_matrix(sm_v)
3434 558 : CALL dbcsr_deallocate_matrix(sm_a)
3435 :
3436 558 : IF (PRESENT(matrix_w)) THEN
3437 2 : CALL dbcsr_deallocate_matrix(sm_b)
3438 2 : CALL dbcsr_deallocate_matrix(sm_c)
3439 2 : CALL dbcsr_deallocate_matrix(sm_d)
3440 2 : CALL dbcsr_deallocate_matrix(sm_e)
3441 2 : CALL dbcsr_deallocate_matrix(sm_f)
3442 2 : CALL dbcsr_deallocate_matrix(sm_k)
3443 2 : CALL dbcsr_deallocate_matrix(sm_i)
3444 2 : CALL dbcsr_deallocate_matrix(sm_j)
3445 : END IF
3446 :
3447 558 : CALL timestop(handle)
3448 :
3449 1674 : END SUBROUTINE tensorial
3450 :
3451 60 : END MODULE dft_plus_u
|
