(git:295f299)
Loading...
Searching...
No Matches
dft_plus_u.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7! **************************************************************************************************
8!> \brief Add the DFT+U contribution to the Hamiltonian matrix
9!> \details The implemented methods refers to:\n
10!> S. L. Dudarev, D. Nguyen Manh, and A. P. Sutton,
11!> Philos. Mag. B \b 75, 613 (1997)\n
12!> S. L. Dudarev et al.,
13!> Phys. Rev. B \b 57, 1505 (1998)
14!> \author Matthias Krack (MK)
15!> \date 14.01.2008
16!> \version 1.0
17! **************************************************************************************************
18MODULE dft_plus_u
19 USE 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
95CONTAINS
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 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 CALL timeset(routinen, handle)
121
122 cpassert(ASSOCIATED(qs_env))
123
124 NULLIFY (input, dft_control)
125
126 logger => cp_get_default_logger()
127
128 CALL get_qs_env(qs_env=qs_env, &
129 input=input, &
130 dft_control=dft_control)
131
132 CALL cite_reference(dudarev1997)
133 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 orthonormal_basis = .false.
139
140 ! Setup print control
141
142 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 (.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 log_filename=.false.)
150
151 ! Select DFT+U method
152
153 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 should_output, output_unit, print_level)
163 END IF
164 CASE (plus_u_mulliken)
165 CALL mulliken(qs_env, orthonormal_basis, matrix_h, &
166 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 should_output, output_unit, print_level)
170 CASE (plus_u_tensorial)
171 CALL tensorial(qs_env, matrix_h, matrix_w, &
172 should_output, output_unit, print_level)
173 CASE DEFAULT
174 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 ignore_should_output=should_output)
179
180 CALL timestop(handle)
181
182 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 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
231 INTEGER, DIMENSION(1) :: iloc
232 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals
233 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigval, q_eigval
240 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: orbq, q_eigvec, q_matrix, q_opp_matrix, &
241 q_work, slam
242 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
243 POINTER :: local_data
244 REAL(kind=dp), DIMENSION(:, :), POINTER :: q_block, q_opp_block, v_block, vhxc_block
245 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 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fm_wmat
249 TYPE(cp_fm_type), POINTER :: fm_s_half
250 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vhxc
251 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
259 TYPE(qs_energy_type), POINTER :: energy
260 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 CALL timeset(routinen, handle)
265
266 debug = .false. ! Set to .TRUE. to print debug information
267
268 NULLIFY (sm_h, sm_p, sm_p_opp, sm_s, sm_w, vhxc_block, sm_vhxc, matrix_vhxc)
269
270 smear = .false.
271 max_scf = -1
272 eps_scf = 1.0e30_dp
273 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 para_env=para_env)
288
289 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
290
291 energy%dft_plus_u = 0.0_dp
292
293 nspin = dft_control%nspins
294 nimg = dft_control%nimages
295 IF (dft_control%mtlr_dft_with_perturbation) THEN
296 IF (.NOT. ASSOCIATED(matrix_vhxc)) THEN
297 cpabort("MTLR requires the projected Hxc matrix for every spin channel.")
298 END IF
299 IF (SIZE(matrix_vhxc) /= nspin) THEN
300 cpabort("The number of projected Hxc matrices does not match the spin channels.")
301 END IF
302 END IF
303
304 IF (nspin == 2) THEN
305 fspin = 1.0_dp
306 ELSE
307 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 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
314 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
315
316 nkind = SIZE(atomic_kind_set)
317
318 ALLOCATE (first_sgf_atom(natom))
319 first_sgf_atom(:) = 0
320
321 CALL get_particle_set(particle_set, qs_kind_set, &
322 first_sgf=first_sgf_atom)
323
324 IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
325 just_energy = .false.
326 ELSE
327 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 DO ikind = 1, nkind
333 NULLIFY (orbitals)
334 CALL get_qs_kind(qs_kind_set(ikind), &
335 hund_j=hund_j, &
336 orbitals=orbitals)
337 IF (hund_j /= 0.0_dp) THEN
338 if_plus_j = .true.
339 IF (ASSOCIATED(orbitals)) THEN
340 CALL cp_abort(__location__, "ENFORCE_OCCUPATION is incompatible with a nonzero Hund J "// &
341 "in the Lowdin DFT+U+J method.")
342 END IF
343 END IF
344 END DO
345
346 NULLIFY (orbitals)
347
348 IF (do_kpoints) THEN
349 fm_wmat => scf_env%scf_work1
350 fmstruct => fm_wmat(1)%matrix_struct
351 ELSE
352 ! Retrieve S^(1/2) from the SCF environment
353 fm_s_half => scf_env%s_half
354 cpassert(ASSOCIATED(fm_s_half))
355 ! work matrices
356 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 name="FULL WORK MATRIX 1")
360 CALL cp_fm_create(matrix=fm_work2, matrix_struct=fmstruct, &
361 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 IF (PRESENT(matrix_w)) THEN
367 IF (do_kpoints) THEN
368 cpabort("Lowdin forces with k-points NYA in DFT+U")
369 END IF
370 CALL cp_fm_create(matrix=fm_sev, matrix_struct=fmstruct)
371 CALL cp_fm_create(matrix=slambda, matrix_struct=fmstruct)
372 ALLOCATE (eigval(nsgf), slam(nsgf, 1))
373 sm_s => matrix_s(1, 1)%matrix
374 CALL copy_dbcsr_to_fm(sm_s, fm_work1)
375 CALL choose_eigv_solver(fm_work1, fm_sev, eigval)
376 !
377 DO i = 1, nsgf
378 IF (eigval(i) > 0._dp) THEN
379 slam(i, 1) = sqrt(eigval(i))
380 ELSE
381 cpabort("S matrix not positive definit")
382 END IF
383 END DO
384 DO i = 1, nsgf
385 CALL cp_fm_set_submatrix(slambda, slam, 1, i, nsgf, 1, 1.0_dp, 0.0_dp)
386 END DO
387 DO i = 1, nsgf
388 CALL cp_fm_set_submatrix(slambda, slam, i, 1, 1, nsgf, 1.0_dp, 1.0_dp, .true.)
389 END DO
390 CALL cp_fm_get_info(slambda, local_data=local_data)
391 DO i = 1, SIZE(local_data, 2)
392 DO j = 1, SIZE(local_data, 1)
393 sij = local_data(j, i)
394 IF (sij > 0.0_dp) sij = 1.0_dp/sij
395 local_data(j, i) = sij
396 END DO
397 END DO
398 DEALLOCATE (eigval, slam)
399 END IF
400
401 ! Calculate S^(1/2)*P*S^(1/2)
402 IF (do_kpoints) THEN
403 cpabort("Lowdin option with k-points NYA in DFT+U")
404 ALLOCATE (orbq(nsgf, nspin))
405 CALL lowdin_kp_trans(kpoints, orbq)
406 DEALLOCATE (orbq)
407 END IF
408
409 ! Create local block diagonal matrices
410 sm_s => matrix_s(1, 1)%matrix
411 CALL dbcsr_get_block_diag(sm_s, sm_q)
412 IF (if_plus_j .AND. nspin == 2) THEN
413 CALL dbcsr_get_block_diag(sm_s, sm_q_opp)
414 END IF
415 CALL dbcsr_get_block_diag(sm_s, sm_v)
416
417 ! Loop over all spins
418 DO ispin = 1, nspin
419
420 CALL dbcsr_set(sm_q, 0.0_dp)
421 IF (if_plus_j .AND. nspin == 2) THEN
422 CALL dbcsr_set(sm_q_opp, 0.0_dp)
423 END IF
424 CALL dbcsr_set(sm_v, 0.0_dp)
425
426 IF (dft_control%mtlr_dft_with_perturbation) THEN
427 IF (ispin == 1) perturbation_strength = 0.9_dp*dft_control%perturbation_strength
428 IF (ispin == 2) perturbation_strength = 1.1_dp*dft_control%perturbation_strength
429 IF (.NOT. ASSOCIATED(matrix_vhxc(ispin)%matrix)) THEN
430 cpabort("MTLR projected Hxc matrix is not initialized.")
431 END IF
432 sm_vhxc => matrix_vhxc(ispin)%matrix
433 END IF
434
435 IF (do_kpoints) THEN
436 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 sm_p => matrix_p(ispin, 1)%matrix
440 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 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 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 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 sm_p_opp => matrix_p(3 - ispin, 1)%matrix
470 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 matrix_c=fm_work2)
481 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 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 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 smear=smear)
515
516 ! Check, if the atoms of this atomic kind need a DFT+U correction
517 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
518 IF (.NOT. dft_plus_u_atom) cycle
519 IF (lu < 0) cycle
520
521 ! Apply U ramping if requested
522 IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
523 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
524 u_minus_j = min(u_minus_j + u_ramping, u_minus_j_target)
525 CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
526 END IF
527 IF (should_output .AND. (output_unit > 0)) THEN
528 WRITE (unit=output_unit, fmt="(T3,A,3X,A,F0.3,A)") &
529 "Kind name: "//trim(adjustl(atomic_kind_name)), &
530 "U(eff) = ", u_minus_j*evolt, " eV"
531 END IF
532 END IF
533
534 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 nshell=nshell)
543
544 ! Count the relevant shell blocks of this atomic kind
545 nsb = 0
546 DO iset = 1, nset
547 DO ishell = 1, nshell(iset)
548 IF (l(ishell, iset) == lu) nsb = nsb + 1
549 END DO
550 END DO
551
552 nsbsize = (2*lu + 1)
553 n = nsb*nsbsize
554
555 ALLOCATE (q_matrix(n, n))
556 q_matrix(:, :) = 0.0_dp
557
558 IF (hund_j /= 0.0_dp .AND. nspin == 2) THEN
559 ALLOCATE (q_opp_matrix(n, n))
560 q_opp_matrix(:, :) = 0.0_dp
561 END IF
562
563 IF (dft_control%mtlr_dft_with_perturbation .AND. &
564 dft_control%mtlr_ikind == ikind) THEN
565 dft_control%vhxc(ispin) = 0.0_dp
566 dft_control%trq(ispin) = 0.0_dp
567 END IF
568
569 ! Print headline if requested
570 IF (should_output .AND. (print_level > low_print_level)) THEN
571 IF (output_unit > 0) THEN
572 ALLOCATE (symbol(nsbsize))
573 DO m = -lu, lu
574 symbol(lu + m + 1) = sgf_symbol(0, lu, m)
575 END DO
576 IF (nspin > 1) THEN
577 WRITE (unit=spin_info, fmt="(A8,I2)") " of spin", ispin
578 ELSE
579 spin_info = ""
580 END IF
581 WRITE (unit=output_unit, fmt="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
582 "DFT+U occupations"//trim(spin_info)//" for the atoms of atomic kind ", ikind, &
583 ": "//trim(atomic_kind_name), &
584 "Atom Shell ", (adjustr(symbol(i)), i=1, nsbsize), " Trace"
585 DEALLOCATE (symbol)
586 END IF
587 END IF
588
589 ! Loop over all atoms of the current atomic kind
590 DO iatom = 1, natom_of_kind
591 atom_a = atom_list(iatom)
592 q_matrix(:, :) = 0.0_dp
593 IF (ALLOCATED(q_opp_matrix)) THEN
594 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 found=found)
603 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 found=found)
609 END IF
610 IF (dft_control%mtlr_dft_with_perturbation .AND. &
611 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 found=found)
618 END IF
619
620 IF (ASSOCIATED(q_block)) THEN
621 ! Calculate energy contribution to E(U)
622 i = 0
623 DO iset = 1, nset
624 DO ishell = 1, nshell(iset)
625 IF (l(ishell, iset) /= lu) cycle
626 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
627 i = i + 1
628 j = 0
629 DO jset = 1, nset
630 DO jshell = 1, nshell(jset)
631 IF (l(jshell, jset) /= lu) cycle
632 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
633 j = j + 1
634 IF (isgf == jsgf) q_matrix(i, j) = q_block(isgf, jsgf)
635 IF (hund_j /= 0.0_dp .AND. nspin == 2 .AND. ASSOCIATED(q_opp_block)) THEN
636 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 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a .AND. &
642 isgf == jsgf) THEN
643 dft_control%vhxc(ispin) = dft_control%vhxc(ispin) + vhxc_block(isgf, isgf)
644 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 IF (ASSOCIATED(orbitals)) THEN
655 IF ((qs_env%scf_env%iter_delta >= eps_scf) .OR. &
656 ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
657 (qs_env%scf_env%iter_count <= max_scf))) THEN
658 ALLOCATE (orb_occ(nsbsize))
659 ALLOCATE (q_eigval(n))
660 q_eigval(:) = 0.0_dp
661 ALLOCATE (q_eigvec(n, n))
662 q_eigvec(:, :) = 0.0_dp
663 norb = SIZE(orbitals)
664 CALL jacobi(q_matrix, q_eigval, q_eigvec)
665 q_matrix(:, :) = 0.0_dp
666 DO isb = 1, nsb
667 trq = 0.0_dp
668 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
669 trq = trq + q_eigval(i)
670 END DO
671 IF (smear) THEN
672 occ = trq/real(norb, kind=dp)
673 ELSE
674 occ = 1.0_dp/fspin
675 END IF
676 orb_occ(:) = .false.
677 iloc = maxloc(q_eigvec(:, isb*nsbsize))
678 jsb = int((iloc(1) - 1)/nsbsize) + 1
679 i = 0
680 i0 = (jsb - 1)*nsbsize + 1
681 iorb = -1000
682 DO j = i0, jsb*nsbsize
683 i = i + 1
684 IF (i > norb) THEN
685 DO m = -lu, lu
686 IF (.NOT. orb_occ(lu + m + 1)) THEN
687 iorb = i0 + lu + m
688 orb_occ(lu + m + 1) = .true.
689 END IF
690 END DO
691 ELSE
692 iorb = i0 + lu + orbitals(i)
693 orb_occ(lu + orbitals(i) + 1) = .true.
694 END IF
695 cpassert(iorb /= -1000)
696 iloc = maxloc(q_eigvec(iorb, :))
697 q_eigval(iloc(1)) = min(occ, trq)
698 q_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1)) ! backtransform left
699 trq = trq - q_eigval(iloc(1))
700 END DO
701 END DO
702 q_matrix(:, :) = matmul(q_matrix, transpose(q_eigvec)) ! backtransform right
703 DEALLOCATE (orb_occ)
704 DEALLOCATE (q_eigval)
705 DEALLOCATE (q_eigvec)
706 END IF
707 END IF ! orbitals associated
708
709 trq = 0.0_dp
710 trq2 = 0.0_dp
711 trqxq_hund_j = 0.0_dp
712 DO i = 1, n
713 trq = trq + q_matrix(i, i)
714 DO j = 1, n
715 trq2 = trq2 + q_matrix(i, j)*q_matrix(j, i)
716 IF (hund_j /= 0.0_dp) THEN
717 IF (nspin == 2) THEN
718 trqxq_hund_j = trqxq_hund_j + q_matrix(i, j)*q_opp_matrix(j, i)
719 ELSE
720 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 trq = fspin*trq
726 trq2 = fspin*fspin*trq2
727
728 ! Calculate energy contribution to E(U)
729 energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
730 IF (hund_j /= 0.0_dp) THEN
731 trqxq_hund_j = fspin*fspin*trqxq_hund_j
732 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 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 found=found)
742 cpassert(ASSOCIATED(v_block))
743
744 i = 0
745 DO iset = 1, nset
746 DO ishell = 1, nshell(iset)
747 IF (l(ishell, iset) /= lu) cycle
748 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
749 i = i + 1
750 j = 0
751 DO jset = 1, nset
752 DO jshell = 1, nshell(jset)
753 IF (l(jshell, jset) /= lu) cycle
754 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
755 j = j + 1
756 IF (isgf == jsgf) THEN
757 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 dft_control%mtlr_ikind == ikind .AND. &
760 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
761 v_block(isgf, isgf) = v_block(isgf, isgf) + perturbation_strength
762 END IF
763 ELSE
764 cpassert(abs(q_matrix(j, i)) < 1.0e-14_dp)
765 IF (hund_j /= 0.0_dp .AND. nspin == 2) THEN
766 cpassert(abs(q_opp_matrix(j, i)) < 1.0e-14_dp)
767 END IF
768 v_block(isgf, jsgf) = -u_minus_j*fspin*q_matrix(j, i)
769 END IF
770 IF (hund_j /= 0.0_dp) THEN
771 IF (nspin == 2) THEN
772 v_block(isgf, jsgf) = v_block(isgf, jsgf) + hund_j*fspin*q_opp_matrix(j, i)
773 ELSE
774 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 IF (should_output .AND. (print_level > low_print_level)) THEN
789 CALL para_env%sum(q_matrix)
790 IF (output_unit > 0) THEN
791 ALLOCATE (q_work(nsb, nsbsize))
792 q_work(:, :) = 0.0_dp
793 DO isb = 1, nsb
794 j = 0
795 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
796 j = j + 1
797 q_work(isb, j) = q_matrix(i, i)
798 END DO
799 END DO
800 DO isb = 1, nsb
801 WRITE (unit=output_unit, fmt="(T3,I6,2X,I6,2X,10F8.3)") &
802 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 "Total", (sum(q_work(:, i)), i=1, nsbsize), sum(q_work)
806 WRITE (unit=output_unit, fmt="(A)") ""
807 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 IF (dft_control%mtlr_dft_with_perturbation .AND. &
863 dft_control%mtlr_ikind == ikind) THEN
864 dft_control%vhxc(ispin) = dft_control%vhxc(ispin)/real(n, dp)
865 CALL para_env%sum(dft_control%vhxc(ispin))
866 CALL para_env%sum(dft_control%trq(ispin))
867 END IF
868
869 IF (ALLOCATED(q_matrix)) THEN
870 DEALLOCATE (q_matrix)
871 END IF
872
873 IF (ALLOCATED(q_opp_matrix)) THEN
874 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 IF (PRESENT(matrix_h)) THEN
881 IF (do_kpoints) THEN
882 cpabort("Lowdin option with k-points NYA in DFT+U")
883 ELSE
884 sm_h => matrix_h(ispin, 1)%matrix
885 CALL cp_dbcsr_sm_fm_multiply(sm_v, fm_s_half, fm_work1, nsgf)
886 CALL cp_fm_transpose(fm_work1, fm_work2)
887 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 IF (PRESENT(matrix_w)) THEN
894
895 sm_p => matrix_p(ispin, 1)%matrix
896 sm_w => matrix_w(ispin, 1)%matrix
897
898 CALL cp_dbcsr_sm_fm_multiply(sm_v, fm_s_half, fm_work1, nsgf)
899 CALL cp_fm_transpose(fm_work1, fm_work2)
900 CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_work2, fm_work1, nsgf)
901 CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, 1.0_dp, fm_work1, fm_sev, 0.0_dp, fm_work2)
902 CALL parallel_gemm('T', 'N', nsgf, nsgf, nsgf, 1.0_dp, fm_sev, fm_work2, 0.0_dp, fm_work1)
903 CALL cp_fm_schur_product(fm_work1, slambda, fm_work2)
904 CALL cp_fm_transpose(fm_work2, fm_work1)
905 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_work1, matrix_b=fm_work2)
906 CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, -2.0_dp, fm_sev, fm_work2, 0.0_dp, fm_work1)
907 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 IF (PRESENT(matrix_w)) THEN
914 CALL cp_fm_release(matrix=fm_sev)
915 CALL cp_fm_release(matrix=slambda)
916 END IF
917
918 ! Collect the energy contributions from all processes
919
920 CALL para_env%sum(energy%dft_plus_u)
921
922 IF (energy%dft_plus_u < 0.0_dp) THEN
923 CALL cp_warn(__location__, &
924 "DFT+U energy contribution is negative possibly due "// &
925 "to unphysical Lowdin charges!")
926 END IF
927
928 ! Release (local) full matrices
929 NULLIFY (fm_s_half)
930 CALL cp_fm_release(matrix=fm_work1)
931 CALL cp_fm_release(matrix=fm_work2)
932
933 ! Release (local) sparse matrices
934 CALL dbcsr_release(sm_q)
935 CALL dbcsr_release(sm_v)
936
937 IF (if_plus_j .AND. nspin == 2) THEN
938 CALL dbcsr_release(sm_q_opp)
939 END IF
940
941 CALL timestop(handle)
942
943 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 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 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 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals
995 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
996 LOGICAL :: debug, dft_plus_u_atom, found, &
997 just_energy, occupation_enforced, smear
998 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: q_eigval
1003 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: q_eigvec, q_matrix, q_work
1004 REAL(kind=dp), DIMENSION(:), POINTER :: nelec
1005 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, p_block, q_block, s_block, &
1006 v_block
1007 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1008 TYPE(atomic_kind_type), POINTER :: kind_a
1009 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1015 TYPE(qs_energy_type), POINTER :: energy
1016 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1017 TYPE(qs_rho_type), POINTER :: rho
1018
1019 CALL timeset(routinen, handle)
1020
1021 debug = .false. ! Set to .TRUE. to print debug information
1022
1023 NULLIFY (atom_list)
1024 NULLIFY (atomic_kind_set)
1025 NULLIFY (qs_kind_set)
1026 NULLIFY (dft_control)
1027 NULLIFY (energy)
1028 NULLIFY (first_sgf)
1029 NULLIFY (h_block)
1030 NULLIFY (matrix_p)
1031 NULLIFY (matrix_s)
1032 NULLIFY (l)
1033 NULLIFY (last_sgf)
1034 NULLIFY (nelec)
1035 NULLIFY (nshell)
1036 NULLIFY (orb_basis_set)
1037 NULLIFY (p_block)
1038 NULLIFY (particle_set)
1039 NULLIFY (q_block)
1040 NULLIFY (rho)
1041 NULLIFY (s_block)
1042 NULLIFY (orbitals)
1043 NULLIFY (sm_h)
1044 NULLIFY (sm_p)
1045 NULLIFY (sm_q)
1046 NULLIFY (sm_s)
1047 NULLIFY (sm_v)
1048 NULLIFY (v_block)
1049 NULLIFY (para_env)
1050
1051 smear = .false.
1052 max_scf = -1
1053 eps_scf = 1.0e30_dp
1054 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 para_env=para_env)
1064
1065 cpassert(ASSOCIATED(atomic_kind_set))
1066 cpassert(ASSOCIATED(dft_control))
1067 cpassert(ASSOCIATED(energy))
1068 cpassert(ASSOCIATED(particle_set))
1069 cpassert(ASSOCIATED(rho))
1070
1071 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 matrix_s_kp=matrix_s)
1077 cpassert(ASSOCIATED(matrix_s))
1078 END IF
1079 nimg = dft_control%nimages
1080
1081 ! Get density matrices in sparse format
1082
1083 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1084
1085 energy%dft_plus_u = 0.0_dp
1086
1087 nspin = dft_control%nspins
1088
1089 IF (nspin == 2) THEN
1090 fspin = 1.0_dp
1091 ELSE
1092 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 natom=natom)
1100
1101 nkind = SIZE(atomic_kind_set)
1102
1103 ALLOCATE (is_plus_u_kind(nkind))
1104 is_plus_u_kind(:) = .false.
1105
1106 IF (PRESENT(matrix_h)) THEN
1107 just_energy = .false.
1108 ELSE
1109 just_energy = .true.
1110 END IF
1111
1112 ! Loop over all spins
1113 DO ispin = 1, nspin
1114
1115 ! Loop over cell images
1116 DO ic = 1, nimg
1117 IF (.NOT. orthonormal_basis) THEN
1118 sm_s => matrix_s(1, ic)%matrix
1119 END IF
1120
1121 IF (PRESENT(matrix_h)) THEN
1122 ! Hamiltonian matrix for spin ispin in sparse format
1123 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 sm_p => matrix_p(ispin, ic)%matrix
1131
1132 IF (.NOT. ASSOCIATED(sm_q)) THEN
1133 ALLOCATE (sm_q)
1134 CALL dbcsr_get_block_diag(sm_p, sm_q)
1135 END IF
1136 CALL dbcsr_set(sm_q, 0.0_dp)
1137
1138 IF (.NOT. ASSOCIATED(sm_v)) THEN
1139 ALLOCATE (sm_v)
1140 CALL dbcsr_get_block_diag(sm_p, sm_v)
1141 END IF
1142 CALL dbcsr_set(sm_v, 0.0_dp)
1143
1144 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 found=found)
1151
1152 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 found=found)
1159 cpassert(ASSOCIATED(q_block))
1160
1161 IF (orthonormal_basis) THEN
1162 ! S is the unit matrix
1163 DO isgf = 1, SIZE(q_block, 1)
1164 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 found=found)
1172 cpassert(ASSOCIATED(s_block))
1173 ! Exploit that P and S are symmetric
1174 DO jsgf = 1, SIZE(p_block, 2)
1175 DO isgf = 1, SIZE(p_block, 1)
1176 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 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 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 smear=smear)
1215
1216 ! Check, if the atoms of this atomic kind need a DFT+U correction
1217
1218 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
1219 IF (.NOT. dft_plus_u_atom) cycle
1220 IF (lu < 0) cycle
1221
1222 ! Apply U ramping if requested
1223
1224 IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
1225 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
1226 u_minus_j = min(u_minus_j + u_ramping, u_minus_j_target)
1227 CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
1228 END IF
1229 IF (should_output .AND. (output_unit > 0)) THEN
1230 WRITE (unit=output_unit, fmt="(T3,A,3X,A,F0.3,A)") &
1231 "Kind name: "//trim(adjustl(atomic_kind_name)), &
1232 "U(eff) = ", u_minus_j*evolt, " eV"
1233 END IF
1234 END IF
1235
1236 IF (u_minus_j == 0.0_dp) cycle
1237
1238 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 nshell=nshell)
1248
1249 ! Count the relevant shell blocks of this atomic kind
1250
1251 nsb = 0
1252 DO iset = 1, nset
1253 DO ishell = 1, nshell(iset)
1254 IF (l(ishell, iset) == lu) nsb = nsb + 1
1255 END DO
1256 END DO
1257
1258 nsbsize = (2*lu + 1)
1259 n = nsb*nsbsize
1260
1261 ALLOCATE (q_matrix(n, n))
1262 q_matrix(:, :) = 0.0_dp
1263
1264 ! Print headline if requested
1265
1266 IF (should_output .AND. (print_level > low_print_level)) THEN
1267 IF (output_unit > 0) THEN
1268 ALLOCATE (symbol(nsbsize))
1269 DO m = -lu, lu
1270 symbol(lu + m + 1) = sgf_symbol(0, lu, m)
1271 END DO
1272 IF (nspin > 1) THEN
1273 WRITE (unit=spin_info, fmt="(A8,I2)") " of spin", ispin
1274 ELSE
1275 spin_info = ""
1276 END IF
1277 WRITE (unit=output_unit, fmt="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
1278 "DFT+U occupations"//trim(spin_info)//" for the atoms of atomic kind ", ikind, &
1279 ": "//trim(atomic_kind_name), &
1280 "Atom Shell ", (adjustr(symbol(i)), i=1, nsbsize), " Trace"
1281 DEALLOCATE (symbol)
1282 END IF
1283 END IF
1284
1285 ! Loop over all atoms of the current atomic kind
1286
1287 DO iatom = 1, natom_of_kind
1288
1289 atom_a = atom_list(iatom)
1290
1291 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 found=found)
1300
1301 ! Calculate energy contribution to E(U)
1302
1303 IF (ASSOCIATED(q_block)) THEN
1304
1305 i = 0
1306 DO iset = 1, nset
1307 DO ishell = 1, nshell(iset)
1308 IF (l(ishell, iset) /= lu) cycle
1309 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1310 i = i + 1
1311 j = 0
1312 DO jset = 1, nset
1313 DO jshell = 1, nshell(jset)
1314 IF (l(jshell, jset) /= lu) cycle
1315 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
1316 j = j + 1
1317 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 IF (ASSOCIATED(orbitals)) THEN
1328 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 ALLOCATE (orb_occ(nsbsize))
1332 ALLOCATE (q_eigval(n))
1333 q_eigval(:) = 0.0_dp
1334 ALLOCATE (q_eigvec(n, n))
1335 q_eigvec(:, :) = 0.0_dp
1336 norb = SIZE(orbitals)
1337 CALL jacobi(q_matrix, q_eigval, q_eigvec)
1338 q_matrix(:, :) = 0.0_dp
1339 IF (nelec(ispin) >= 0.5_dp) THEN
1340 trq = nelec(ispin)/sum(q_eigval(1:n))
1341 q_eigval(1:n) = trq*q_eigval(1:n)
1342 END IF
1343 DO isb = 1, nsb
1344 trq = 0.0_dp
1345 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
1346 trq = trq + q_eigval(i)
1347 END DO
1348 IF (smear) THEN
1349 occ = trq/real(norb, kind=dp)
1350 ELSE
1351 occ = 1.0_dp/fspin
1352 END IF
1353 orb_occ(:) = .false.
1354 iloc = maxloc(q_eigvec(:, isb*nsbsize))
1355 jsb = int((iloc(1) - 1)/nsbsize) + 1
1356 i = 0
1357 i0 = (jsb - 1)*nsbsize + 1
1358 iorb = -1000
1359 DO j = i0, jsb*nsbsize
1360 i = i + 1
1361 IF (i > norb) THEN
1362 DO m = -lu, lu
1363 IF (.NOT. orb_occ(lu + m + 1)) THEN
1364 iorb = i0 + lu + m
1365 orb_occ(lu + m + 1) = .true.
1366 END IF
1367 END DO
1368 ELSE
1369 iorb = i0 + lu + orbitals(i)
1370 orb_occ(lu + orbitals(i) + 1) = .true.
1371 END IF
1372 cpassert(iorb /= -1000)
1373 iloc = maxloc(q_eigvec(iorb, :))
1374 q_eigval(iloc(1)) = min(occ, trq)
1375 q_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1)) ! backtransform left
1376 trq = trq - q_eigval(iloc(1))
1377 END DO
1378 END DO
1379 q_matrix(:, :) = matmul(q_matrix, transpose(q_eigvec)) ! backtransform right
1380 DEALLOCATE (orb_occ)
1381 DEALLOCATE (q_eigval)
1382 DEALLOCATE (q_eigvec)
1383 occupation_enforced = .true.
1384 END IF
1385 END IF ! orbitals associated
1386
1387 trq = 0.0_dp
1388 trq2 = 0.0_dp
1389
1390 DO i = 1, n
1391 trq = trq + q_matrix(i, i)
1392 DO j = 1, n
1393 trq2 = trq2 + q_matrix(i, j)*q_matrix(j, i)
1394 END DO
1395 END DO
1396
1397 trq = fspin*trq
1398 trq2 = fspin*fspin*trq2
1399
1400 ! Calculate energy contribution to E(U)
1401
1402 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 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 found=found)
1413 cpassert(ASSOCIATED(v_block))
1414
1415 i = 0
1416 DO iset = 1, nset
1417 DO ishell = 1, nshell(iset)
1418 IF (l(ishell, iset) /= lu) cycle
1419 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1420 i = i + 1
1421 j = 0
1422 DO jset = 1, nset
1423 DO jshell = 1, nshell(jset)
1424 IF (l(jshell, jset) /= lu) cycle
1425 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
1426 j = j + 1
1427 IF (isgf == jsgf) THEN
1428 v_block(isgf, isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i, i))
1429 ELSE
1430 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 IF (should_output .AND. (print_level > low_print_level)) THEN
1446 CALL para_env%sum(q_matrix)
1447 IF (output_unit > 0) THEN
1448 ALLOCATE (q_work(nsb, nsbsize))
1449 q_work(:, :) = 0.0_dp
1450 DO isb = 1, nsb
1451 j = 0
1452 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
1453 j = j + 1
1454 q_work(isb, j) = q_matrix(i, i)
1455 END DO
1456 END DO
1457 DO isb = 1, nsb
1458 WRITE (unit=output_unit, fmt="(T3,I6,2X,I6,2X,10F8.3)") &
1459 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 "Total", (sum(q_work(:, i)), i=1, nsbsize), sum(q_work)
1463 WRITE (unit=output_unit, fmt="(A)") ""
1464 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 IF (ALLOCATED(q_matrix)) THEN
1520 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 IF (ASSOCIATED(sm_h)) THEN
1528
1529 DO ikind = 1, nkind
1530
1531 IF (.NOT. is_plus_u_kind(ikind)) cycle
1532
1533 kind_a => atomic_kind_set(ikind)
1534
1535 CALL get_atomic_kind(atomic_kind=kind_a, &
1536 atom_list=atom_list, &
1537 natom=natom_of_kind)
1538
1539 DO iatom = 1, natom_of_kind
1540
1541 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 found=found)
1548
1549 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 found=found)
1556 cpassert(ASSOCIATED(v_block))
1557
1558 IF (orthonormal_basis) THEN
1559 DO isgf = 1, SIZE(h_block, 1)
1560 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 found=found)
1568 cpassert(ASSOCIATED(s_block))
1569 DO jsgf = 1, SIZE(h_block, 2)
1570 DO isgf = 1, SIZE(h_block, 1)
1571 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 CALL para_env%sum(energy%dft_plus_u)
1589
1590 IF (energy%dft_plus_u < 0.0_dp) THEN
1591 IF (.NOT. occupation_enforced) THEN
1592 CALL cp_warn(__location__, &
1593 "DFT+U energy contribution is negative possibly due "// &
1594 "to unphysical Mulliken charges!")
1595 END IF
1596 END IF
1597
1598 CALL dbcsr_deallocate_matrix(sm_q)
1599 CALL dbcsr_deallocate_matrix(sm_v)
1600
1601 CALL timestop(handle)
1602
1603 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 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
1659 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell
1660 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dedq, trps
1665 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: q_ii
1666 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, p_block, s_block, w_block
1667 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1668 TYPE(dbcsr_iterator_type) :: iter
1669 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1675 TYPE(qs_energy_type), POINTER :: energy
1676 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1677 TYPE(qs_rho_type), POINTER :: rho
1678
1679 CALL timeset(routinen, handle)
1680
1681 NULLIFY (atom_list)
1682 NULLIFY (atomic_kind_set)
1683 NULLIFY (qs_kind_set)
1684 NULLIFY (dft_control)
1685 NULLIFY (energy)
1686 NULLIFY (first_sgf)
1687 NULLIFY (h_block)
1688 NULLIFY (matrix_p)
1689 NULLIFY (matrix_s)
1690 NULLIFY (l)
1691 NULLIFY (last_sgf)
1692 NULLIFY (nshell)
1693 NULLIFY (orb_basis_set)
1694 NULLIFY (p_block)
1695 NULLIFY (particle_set)
1696 NULLIFY (rho)
1697 NULLIFY (s_block)
1698 NULLIFY (sm_h)
1699 NULLIFY (sm_p)
1700 NULLIFY (sm_s)
1701 NULLIFY (w_block)
1702 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 para_env=para_env)
1712
1713 cpassert(ASSOCIATED(atomic_kind_set))
1714 cpassert(ASSOCIATED(dft_control))
1715 cpassert(ASSOCIATED(energy))
1716 cpassert(ASSOCIATED(particle_set))
1717 cpassert(ASSOCIATED(rho))
1718
1719 IF (orthonormal_basis) THEN
1720 NULLIFY (sm_s)
1721 ELSE
1722 ! Get overlap matrix in sparse format
1723 CALL get_qs_env(qs_env=qs_env, &
1724 matrix_s_kp=matrix_s)
1725 cpassert(ASSOCIATED(matrix_s))
1726 END IF
1727
1728 ! Get density matrices in sparse format
1729
1730 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1731
1732 energy%dft_plus_u = 0.0_dp
1733
1734 nspin = dft_control%nspins
1735 nimg = dft_control%nimages
1736
1737 IF (nspin == 2) THEN
1738 fspin = 1.0_dp
1739 ELSE
1740 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 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
1747 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
1748
1749 nkind = SIZE(atomic_kind_set)
1750
1751 ALLOCATE (first_sgf_atom(natom))
1752 first_sgf_atom(:) = 0
1753
1754 CALL get_particle_set(particle_set, qs_kind_set, &
1755 first_sgf=first_sgf_atom)
1756
1757 ALLOCATE (trps(nsgf))
1758 trps(:) = 0.0_dp
1759
1760 IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
1761 ALLOCATE (dedq(nsgf))
1762 just_energy = .false.
1763 ELSE
1764 just_energy = .true.
1765 END IF
1766
1767 ! Loop over all spins
1768
1769 DO ispin = 1, nspin
1770
1771 IF (.NOT. just_energy) dedq(:) = 0.0_dp
1772
1773 ! Calculate Trace(P*S) assuming symmetric matrices
1774
1775 trps(:) = 0.0_dp
1776
1777 DO ic = 1, nimg
1778 IF (orthonormal_basis) THEN
1779 NULLIFY (sm_s)
1780 ELSE
1781 sm_s => matrix_s(1, ic)%matrix
1782 END IF
1783 sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
1784
1785 CALL dbcsr_iterator_start(iter, sm_p)
1786
1787 DO WHILE (dbcsr_iterator_blocks_left(iter))
1788
1789 CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
1790
1791 IF (orthonormal_basis) THEN
1792
1793 IF (iatom /= jatom) cycle
1794
1795 IF (ASSOCIATED(p_block)) THEN
1796 sgf = first_sgf_atom(iatom)
1797 DO isgf = 1, SIZE(p_block, 1)
1798 trps(sgf) = trps(sgf) + p_block(isgf, isgf)
1799 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 found=found)
1810 cpassert(ASSOCIATED(s_block))
1811
1812 sgf = first_sgf_atom(jatom)
1813 DO jsgf = 1, SIZE(p_block, 2)
1814 DO isgf = 1, SIZE(p_block, 1)
1815 trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
1816 END DO
1817 sgf = sgf + 1
1818 END DO
1819
1820 IF (iatom /= jatom) THEN
1821 sgf = first_sgf_atom(iatom)
1822 DO isgf = 1, SIZE(p_block, 1)
1823 DO jsgf = 1, SIZE(p_block, 2)
1824 trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
1825 END DO
1826 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 CALL dbcsr_iterator_stop(iter)
1835
1836 END DO ! cell images
1837
1838 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 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 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 eps_u_ramping=eps_u_ramping)
1867
1868 ! Check, if this atom needs a DFT+U correction
1869
1870 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
1871 IF (.NOT. dft_plus_u_atom) cycle
1872 IF (lu < 0) cycle
1873
1874 ! Apply U ramping if requested
1875
1876 IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
1877 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
1878 u_minus_j = min(u_minus_j + u_ramping, u_minus_j_target)
1879 CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
1880 END IF
1881 IF (should_output .AND. (output_unit > 0)) THEN
1882 WRITE (unit=output_unit, fmt="(T3,A,3X,A,F0.3,A)") &
1883 "Kind name: "//trim(adjustl(atomic_kind_name)), &
1884 "U(eff) = ", u_minus_j*evolt, " eV"
1885 END IF
1886 END IF
1887
1888 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 nshell=nshell)
1898
1899 ! Count the relevant shell blocks of this atomic kind
1900
1901 nsb = 0
1902 DO iset = 1, nset
1903 DO ishell = 1, nshell(iset)
1904 IF (l(ishell, iset) == lu) nsb = nsb + 1
1905 END DO
1906 END DO
1907
1908 ALLOCATE (q_ii(nsb, 2*lu + 1))
1909
1910 ! Print headline if requested
1911
1912 IF (should_output .AND. (print_level > low_print_level)) THEN
1913 IF (output_unit > 0) THEN
1914 ALLOCATE (symbol(2*lu + 1))
1915 DO m = -lu, lu
1916 symbol(lu + m + 1) = sgf_symbol(0, lu, m)
1917 END DO
1918 IF (nspin > 1) THEN
1919 WRITE (unit=spin_info, fmt="(A8,I2)") " of spin", ispin
1920 ELSE
1921 spin_info = ""
1922 END IF
1923 WRITE (unit=output_unit, fmt="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
1924 "DFT+U occupations"//trim(spin_info)//" for the atoms of atomic kind ", ikind, &
1925 ": "//trim(atomic_kind_name), &
1926 "Atom Shell ", (adjustr(symbol(i)), i=1, 2*lu + 1), " Trace"
1927 DEALLOCATE (symbol)
1928 END IF
1929 END IF
1930
1931 ! Loop over all atoms of the current atomic kind
1932
1933 DO iatom = 1, natom_of_kind
1934
1935 atom_a = atom_list(iatom)
1936
1937 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 found=found)
1946
1947 ! Calculate E(U) and dE(U)/dq
1948
1949 IF (ASSOCIATED(p_block)) THEN
1950
1951 sgf = first_sgf_atom(atom_a)
1952
1953 isb = 0
1954 DO iset = 1, nset
1955 DO ishell = 1, nshell(iset)
1956 IF (l(ishell, iset) == lu) THEN
1957 isb = isb + 1
1958 i = 0
1959 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1960 q = fspin*trps(sgf)
1961 i = i + 1
1962 q_ii(isb, i) = q
1963 energy%dft_plus_u = energy%dft_plus_u + &
1964 0.5_dp*u_minus_j*(q - q**2)/fspin
1965 IF (.NOT. just_energy) THEN
1966 dedq(sgf) = dedq(sgf) + u_minus_j*(0.5_dp - q)
1967 END IF
1968 sgf = sgf + 1
1969 END DO ! next contracted spherical Gaussian function "isgf"
1970 ELSE
1971 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 IF (should_output .AND. (print_level > low_print_level)) THEN
1981 CALL para_env%sum(q_ii)
1982 IF (output_unit > 0) THEN
1983 DO isb = 1, nsb
1984 WRITE (unit=output_unit, fmt="(T3,I6,2X,I6,2X,10F8.3)") &
1985 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 "Total", (sum(q_ii(:, i)), i=1, 2*lu + 1), sum(q_ii)
1989 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 IF (ALLOCATED(q_ii)) THEN
1996 DEALLOCATE (q_ii)
1997 END IF
1998
1999 END DO ! next atomic kind "ikind"
2000
2001 IF (.NOT. just_energy) THEN
2002 CALL para_env%sum(dedq)
2003 END IF
2004
2005 ! Add V(i,j)[U] to V(i,j)[DFT]
2006
2007 IF (PRESENT(matrix_h)) THEN
2008
2009 DO ic = 1, nimg
2010 IF (orthonormal_basis) THEN
2011 NULLIFY (sm_s)
2012 ELSE
2013 sm_s => matrix_s(1, ic)%matrix
2014 END IF
2015 sm_h => matrix_h(ispin, ic)%matrix
2016
2017 CALL dbcsr_iterator_start(iter, sm_h)
2018
2019 DO WHILE (dbcsr_iterator_blocks_left(iter))
2020
2021 CALL dbcsr_iterator_next_block(iter, iatom, jatom, h_block)
2022
2023 IF (orthonormal_basis) THEN
2024
2025 IF (iatom /= jatom) cycle
2026
2027 IF (ASSOCIATED(h_block)) THEN
2028 sgf = first_sgf_atom(iatom)
2029 DO isgf = 1, SIZE(h_block, 1)
2030 h_block(isgf, isgf) = h_block(isgf, isgf) + dedq(sgf)
2031 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 found=found)
2044 cpassert(ASSOCIATED(s_block))
2045
2046 ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
2047
2048 sgf = first_sgf_atom(iatom)
2049
2050 DO isgf = 1, SIZE(h_block, 1)
2051 IF (dedq(sgf) /= 0.0_dp) THEN
2052 v = 0.5_dp*dedq(sgf)
2053 DO jsgf = 1, SIZE(h_block, 2)
2054 h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
2055 END DO
2056 END IF
2057 sgf = sgf + 1
2058 END DO
2059
2060 sgf = first_sgf_atom(jatom)
2061
2062 DO jsgf = 1, SIZE(h_block, 2)
2063 IF (dedq(sgf) /= 0.0_dp) THEN
2064 v = 0.5_dp*dedq(sgf)
2065 DO isgf = 1, SIZE(h_block, 1)
2066 h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
2067 END DO
2068 END IF
2069 sgf = sgf + 1
2070 END DO
2071
2072 END IF ! orthonormal basis set
2073
2074 END DO ! Next atom "iatom"
2075
2076 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 IF (PRESENT(matrix_w) .AND. (.NOT. orthonormal_basis)) THEN
2087
2088 DO ic = 1, nimg
2089 sm_s => matrix_s(1, ic)%matrix
2090 sm_p => matrix_p(ispin, ic)%matrix
2091 sm_w => matrix_w(ispin, ic)%matrix
2092
2093 CALL dbcsr_iterator_start(iter, sm_p)
2094
2095 DO WHILE (dbcsr_iterator_blocks_left(iter))
2096
2097 CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
2098
2099 ! Skip the diagonal blocks of the W matrix
2100
2101 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 found=found)
2110 cpassert(ASSOCIATED(w_block))
2111
2112 ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
2113
2114 sgf = first_sgf_atom(iatom)
2115
2116 DO isgf = 1, SIZE(w_block, 1)
2117 IF (dedq(sgf) /= 0.0_dp) THEN
2118 v = -0.5_dp*dedq(sgf)
2119 DO jsgf = 1, SIZE(w_block, 2)
2120 w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
2121 END DO
2122 END IF
2123 sgf = sgf + 1
2124 END DO
2125
2126 sgf = first_sgf_atom(jatom)
2127
2128 DO jsgf = 1, SIZE(w_block, 2)
2129 IF (dedq(sgf) /= 0.0_dp) THEN
2130 v = -0.5_dp*dedq(sgf)
2131 DO isgf = 1, SIZE(w_block, 1)
2132 w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
2133 END DO
2134 END IF
2135 sgf = sgf + 1
2136 END DO
2137
2138 END DO ! next block node "jatom"
2139
2140 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 CALL para_env%sum(energy%dft_plus_u)
2151
2152 IF (energy%dft_plus_u < 0.0_dp) THEN
2153 CALL cp_warn(__location__, &
2154 "DFT+U energy contribution is negative possibly due "// &
2155 "to unphysical Mulliken charges!")
2156 END IF
2157
2158 ! Release local work storage
2159
2160 IF (ALLOCATED(first_sgf_atom)) THEN
2161 DEALLOCATE (first_sgf_atom)
2162 END IF
2163
2164 IF (ALLOCATED(trps)) THEN
2165 DEALLOCATE (trps)
2166 END IF
2167
2168 IF (ALLOCATED(dedq)) THEN
2169 DEALLOCATE (dedq)
2170 END IF
2171
2172 CALL timestop(handle)
2173
2174 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 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 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 INTEGER, ALLOCATABLE :: first_sgf_atom(:)
2206 INTEGER, DIMENSION(1) :: iloc
2207 INTEGER, DIMENSION(:), POINTER :: atom_list, atom_list2, nshell, orbitals
2208 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
2209 LOGICAL :: debug, dft_plus_u_atom, found, &
2210 just_energy, smear
2211 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 REAL(kind=dp), ALLOCATABLE :: a(:), q_a_matrix(:, :), q_b_matrix(:, :), q_eigval(:), &
2215 q_eigvec(:, :), q_work(:, :), s_inv_matrix(:, :), s_matrix(:, :), tmp_matrix1(:, :), &
2216 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 TYPE(block3d), ALLOCATABLE :: a_matrix(:, :), d_matrix(:, :), &
2225 e_matrix(:, :), i_matrix(:, :), &
2226 j_matrix(:, :)
2227 TYPE(block4d), ALLOCATABLE :: dd_matrix(:, :), ee_matrix(:, :), &
2228 ff_matrix(:, :)
2229 REAL(kind=dp), DIMENSION(:), POINTER :: ao_coef
2230 REAL(kind=dp), DIMENSION(:, :), POINTER :: a_block, d_block, e_block, f_block, &
2231 i_block, j_block, q_a_block, q_b_block, &
2232 s_block, tmp_block, v_block, vhxc_block
2233
2234 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2235 TYPE(mp_para_env_type), POINTER :: para_env
2236 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s, matrix_vhxc
2237 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2245 TYPE(qs_energy_type), POINTER :: energy
2246 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2247 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 CALL timeset(routinen, handle)
2252
2253 debug = .false.
2254 smear = .false.
2255 max_scf = -1
2256 eps_scf = 1.0e30_dp
2257
2258 NULLIFY (atom_list)
2259 NULLIFY (atom_list2)
2260 NULLIFY (atomic_kind_set)
2261 NULLIFY (qs_kind_set)
2262 NULLIFY (dft_control)
2263 NULLIFY (energy)
2264 NULLIFY (force)
2265 NULLIFY (first_sgf)
2266 NULLIFY (matrix_p)
2267 NULLIFY (matrix_s)
2268 NULLIFY (matrix_vhxc)
2269 NULLIFY (l)
2270 NULLIFY (last_sgf)
2271 NULLIFY (nshell)
2272 NULLIFY (orb_basis_set)
2273 NULLIFY (orbitals)
2274 NULLIFY (ao_coef)
2275 NULLIFY (particle_set)
2276 NULLIFY (q_a_block)
2277 NULLIFY (q_b_block)
2278 NULLIFY (vhxc_block)
2279 NULLIFY (rho)
2280 NULLIFY (scf_env)
2281 NULLIFY (sm_h)
2282 NULLIFY (sm_p)
2283 NULLIFY (sm_q_a)
2284 NULLIFY (sm_q_b)
2285 NULLIFY (sm_s)
2286 NULLIFY (sm_v)
2287 NULLIFY (sm_a)
2288 NULLIFY (sm_b)
2289 NULLIFY (sm_c)
2290 NULLIFY (sm_d)
2291 NULLIFY (sm_e)
2292 NULLIFY (sm_f)
2293 NULLIFY (sm_i)
2294 NULLIFY (sm_j)
2295 NULLIFY (sm_k)
2296 NULLIFY (sm_mid)
2297 NULLIFY (sm_vhxc)
2298 NULLIFY (v_block)
2299 NULLIFY (a_block)
2300 NULLIFY (d_block)
2301 NULLIFY (e_block)
2302 NULLIFY (f_block)
2303 NULLIFY (i_block)
2304 NULLIFY (j_block)
2305 NULLIFY (para_env)
2306 NULLIFY (s_block)
2307 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 para_env=para_env)
2321
2322 cpassert(ASSOCIATED(atomic_kind_set))
2323 cpassert(ASSOCIATED(dft_control))
2324 cpassert(ASSOCIATED(energy))
2325 cpassert(ASSOCIATED(matrix_s))
2326 cpassert(ASSOCIATED(particle_set))
2327 cpassert(ASSOCIATED(rho))
2328
2329 sm_s => matrix_s(1)%matrix
2330 CALL qs_rho_get(rho, rho_ao=matrix_p)
2331
2332 energy%dft_plus_u = 0.0_dp
2333
2334 nspin = dft_control%nspins
2335 IF (dft_control%mtlr_dft_with_perturbation) THEN
2336 IF (.NOT. ASSOCIATED(matrix_vhxc)) THEN
2337 cpabort("MTLR requires the projected Hxc matrix for every spin channel.")
2338 END IF
2339 IF (SIZE(matrix_vhxc) /= nspin) THEN
2340 cpabort("The number of projected Hxc matrices does not match the spin channels.")
2341 END IF
2342 END IF
2343
2344 IF (nspin == 2) THEN
2345 fspin = 1.0_dp
2346 ELSE
2347 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 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
2354 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
2355
2356 nkind = SIZE(atomic_kind_set)
2357
2358 ALLOCATE (first_sgf_atom(natom))
2359 first_sgf_atom(:) = 0
2360
2361 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf_atom)
2362
2363 IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
2364 just_energy = .false.
2365 ELSE
2366 just_energy = .true.
2367 END IF
2368
2369 CALL dbcsr_init_p(sm_mid)
2370 CALL dbcsr_create(sm_mid, template=sm_s, matrix_type="N")
2371
2372 ! Create local block diagonal matrices
2373
2374 ALLOCATE (sm_q_a)
2375 CALL dbcsr_get_block_diag(sm_s, sm_q_a)
2376
2377 ALLOCATE (sm_q_b)
2378 CALL dbcsr_get_block_diag(sm_s, sm_q_b)
2379
2380 CALL dbcsr_set(sm_q_a, 0.0_dp)
2381 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 DO ispin = 1, nspin
2386
2387 sm_p => matrix_p(ispin)%matrix
2388
2389 IF (nspin == 2) THEN
2390 IF (ispin == 1) THEN
2391 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2392 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_a)
2393 ELSE
2394 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2395 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_b)
2396 END IF
2397 ELSE
2398 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_p, 0.0_dp, sm_mid)
2399 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_mid, sm_s, 0.0_dp, sm_q_a)
2400 CALL dbcsr_copy(sm_q_b, sm_q_a)
2401 END IF
2402
2403 END DO
2404
2405 ALLOCATE (sm_v)
2406 CALL dbcsr_get_block_diag(sm_s, sm_v)
2407
2408 IF (PRESENT(matrix_w)) THEN
2409
2410 CALL dbcsr_init_p(sm_b)
2411 CALL dbcsr_create(sm_b, template=matrix_s(1)%matrix, matrix_type="N")
2412
2413 CALL dbcsr_init_p(sm_d)
2414 CALL dbcsr_create(sm_d, template=matrix_s(1)%matrix, matrix_type="N")
2415
2416 CALL dbcsr_init_p(sm_c)
2417 CALL dbcsr_create(sm_c, template=matrix_s(1)%matrix, matrix_type="N")
2418
2419 CALL dbcsr_init_p(sm_e)
2420 CALL dbcsr_create(sm_e, template=matrix_s(1)%matrix, matrix_type="N")
2421
2422 CALL dbcsr_init_p(sm_f)
2423 CALL dbcsr_create(sm_f, template=matrix_s(1)%matrix, matrix_type="N")
2424
2425 CALL dbcsr_init_p(sm_i)
2426 CALL dbcsr_create(sm_i, template=matrix_s(1)%matrix, matrix_type="N")
2427
2428 CALL dbcsr_init_p(sm_j)
2429 CALL dbcsr_create(sm_j, template=matrix_s(1)%matrix, matrix_type="N")
2430
2431 CALL dbcsr_init_p(sm_k)
2432 CALL dbcsr_create(sm_k, template=matrix_s(1)%matrix, matrix_type="N")
2433
2434 END IF
2435
2436 ALLOCATE (sm_a)
2437 CALL dbcsr_get_block_diag(sm_s, sm_a)
2438
2439 ! Loop over all spins
2440
2441 DO ispin = 1, nspin
2442
2443 IF (PRESENT(matrix_h)) THEN
2444 sm_h => matrix_h(ispin, 1)%matrix
2445 ELSE
2446 NULLIFY (sm_h)
2447 END IF
2448
2449 IF (PRESENT(matrix_w)) THEN
2450 sm_w => matrix_w(ispin, 1)%matrix
2451 ELSE
2452 NULLIFY (sm_w)
2453 END IF
2454
2455 CALL dbcsr_set(sm_v, 0.0_dp)
2456 CALL dbcsr_set(sm_a, 0.0_dp)
2457
2458 IF (dft_control%mtlr_dft_with_perturbation) THEN
2459 IF (ispin == 1) perturbation_strength = 0.9_dp*dft_control%perturbation_strength
2460 IF (ispin == 2) perturbation_strength = 1.1_dp*dft_control%perturbation_strength
2461 IF (.NOT. ASSOCIATED(matrix_vhxc(ispin)%matrix)) THEN
2462 cpabort("MTLR projected Hxc matrix is not initialized.")
2463 END IF
2464 sm_vhxc => matrix_vhxc(ispin)%matrix
2465 END IF
2466
2467 ! Loop over all atomic kinds
2468
2469 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 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 smear=smear)
2491
2492 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
2493 IF (.NOT. dft_plus_u_atom) cycle
2494 IF (lu < 0) cycle
2495 IF (u_minus_j == 0.0_dp .AND. hund_j == 0.0_dp) cycle
2496 IF (.NOT. ASSOCIATED(ao_coef)) THEN
2497 cpabort("Tensorial DFT+U requires associated ao_coef")
2498 END IF
2499
2500 ALLOCATE (a(SIZE(ao_coef)))
2501 a(:) = ao_coef(:)
2502
2503 ! Apply U ramping if requested
2504
2505 IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
2506 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
2507 u_minus_j = min(u_minus_j + u_ramping, u_minus_j_target)
2508 CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
2509 END IF
2510
2511 IF (should_output .AND. (output_unit > 0)) THEN
2512 WRITE (unit=output_unit, fmt="(T3,A,3X,A,F0.3,A)") &
2513 "Kind name: "//trim(adjustl(atomic_kind_name)), &
2514 "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 nshell=nshell)
2524
2525 nsb = 1
2526 nsbsize = (2*lu + 1)
2527 n = nsb*nsbsize
2528
2529 ALLOCATE (q_a_matrix(n, n))
2530 ALLOCATE (q_b_matrix(n, n))
2531 ALLOCATE (s_matrix(n, n))
2532 ALLOCATE (s_inv_matrix(n, n))
2533 ALLOCATE (v_matrix(n, n))
2534
2535 q_a_matrix(:, :) = 0.0_dp
2536 q_b_matrix(:, :) = 0.0_dp
2537 s_matrix(:, :) = 0.0_dp
2538 s_inv_matrix(:, :) = 0.0_dp
2539 v_matrix(:, :) = 0.0_dp
2540
2541 ! Print headline if requested
2542
2543 IF (should_output .AND. (print_level > low_print_level)) THEN
2544 IF (output_unit > 0) THEN
2545 ALLOCATE (symbol(nsbsize))
2546 DO m = -lu, lu
2547 symbol(lu + m + 1) = sgf_symbol(0, lu, m)
2548 END DO
2549
2550 IF (nspin > 1) THEN
2551 WRITE (unit=spin_info, fmt="(A8,I2)") " of spin", ispin
2552 ELSE
2553 spin_info = ""
2554 END IF
2555
2556 WRITE (unit=output_unit, fmt="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
2557 "DFT+U occupations"//trim(spin_info)//" for the atoms of atomic kind ", ikind, &
2558 ": "//trim(atomic_kind_name), &
2559 "Atom Shell ", (adjustr(symbol(i)), i=1, nsbsize), " Trace"
2560
2561 DEALLOCATE (symbol)
2562 END IF
2563 END IF
2564
2565 IF (dft_control%mtlr_dft_with_perturbation .AND. &
2566 dft_control%mtlr_ikind == ikind) THEN
2567 dft_control%vhxc(ispin) = 0.0_dp
2568 dft_control%trq(ispin) = 0.0_dp
2569 ALLOCATE (vhxc_matrix(n, n))
2570 vhxc_matrix(:, :) = 0.0_dp
2571 END IF
2572
2573 ! Loop over all atoms of the current atomic kind
2574
2575 DO iatom = 1, natom_of_kind
2576
2577 atom_a = atom_list(iatom)
2578
2579 q_a_matrix(:, :) = 0.0_dp
2580 q_b_matrix(:, :) = 0.0_dp
2581 s_matrix(:, :) = 0.0_dp
2582 s_inv_matrix(:, :) = 0.0_dp
2583 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 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 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 found=found)
2602
2603 IF (dft_control%mtlr_dft_with_perturbation .AND. &
2604 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 found=found)
2611 END IF
2612
2613 IF (ASSOCIATED(q_a_block) .AND. ASSOCIATED(q_b_block) .AND. ASSOCIATED(s_block)) THEN
2614
2615 isb = 0
2616 DO iset = 1, nset
2617 DO ishell = 1, nshell(iset)
2618 IF (l(ishell, iset) /= lu) cycle
2619 isb = isb + 1
2620 i = 0
2621 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
2622 i = i + 1
2623 jsb = 0
2624 DO jset = 1, nset
2625 DO jshell = 1, nshell(jset)
2626 IF (l(jshell, jset) /= lu) cycle
2627 jsb = jsb + 1
2628 j = 0
2629 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
2630 j = j + 1
2631 q_a_matrix(i, j) = q_a_matrix(i, j) + q_a_block(isgf, jsgf)*a(isb)*a(jsb)
2632 q_b_matrix(i, j) = q_b_matrix(i, j) + q_b_block(isgf, jsgf)*a(isb)*a(jsb)
2633 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 dft_control%mtlr_ikind == ikind .AND. &
2637 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
2638 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 CALL invert_matrix(s_matrix, s_inv_matrix, some_real)
2648
2649 q_a_matrix(:, :) = matmul(q_a_matrix, s_inv_matrix)
2650 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 dft_control%mtlr_ikind == ikind .AND. &
2655 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
2656 vhxc_matrix(:, :) = matmul(vhxc_matrix, s_inv_matrix)
2657 DO i = 1, n
2658 dft_control%vhxc(ispin) = dft_control%vhxc(ispin) + vhxc_matrix(i, i)
2659 IF (ispin == 1) dft_control%trq(1) = dft_control%trq(1) + q_a_matrix(i, i)
2660 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 IF (ASSOCIATED(orbitals)) THEN
2667 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 ALLOCATE (orb_occ(nsbsize))
2672 ALLOCATE (q_eigval(n))
2673 ALLOCATE (q_eigvec(n, n))
2674
2675 q_eigval(:) = 0.0_dp
2676 q_eigvec(:, :) = 0.0_dp
2677 norb = SIZE(orbitals)
2678
2679 IF (ispin == 1) THEN
2680 CALL jacobi(q_a_matrix, q_eigval, q_eigvec)
2681 q_a_matrix(:, :) = 0.0_dp
2682 ELSE
2683 CALL jacobi(q_b_matrix, q_eigval, q_eigvec)
2684 q_b_matrix(:, :) = 0.0_dp
2685 END IF
2686
2687 DO isb = 1, nsb
2688 trq = 0.0_dp
2689 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
2690 trq = trq + q_eigval(i)
2691 END DO
2692
2693 IF (smear) THEN
2694 occ = trq/real(norb, kind=dp)
2695 ELSE
2696 occ = 1.0_dp/fspin
2697 END IF
2698
2699 orb_occ(:) = .false.
2700 iloc = maxloc(q_eigvec(:, isb*nsbsize))
2701 jsb = int((iloc(1) - 1)/nsbsize) + 1
2702 i = 0
2703 i0 = (jsb - 1)*nsbsize + 1
2704 iorb = -1000
2705
2706 DO j = i0, jsb*nsbsize
2707 i = i + 1
2708 IF (i > norb) THEN
2709 DO m = -lu, lu
2710 IF (.NOT. orb_occ(lu + m + 1)) THEN
2711 iorb = i0 + lu + m
2712 orb_occ(lu + m + 1) = .true.
2713 END IF
2714 END DO
2715 ELSE
2716 iorb = i0 + lu + orbitals(i)
2717 orb_occ(lu + orbitals(i) + 1) = .true.
2718 END IF
2719
2720 cpassert(iorb /= -1000)
2721
2722 iloc = maxloc(q_eigvec(iorb, :))
2723 q_eigval(iloc(1)) = min(occ, trq)
2724
2725 IF (ispin == 1) THEN
2726 q_a_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1))
2727 ELSE
2728 q_b_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1))
2729 END IF
2730
2731 trq = trq - q_eigval(iloc(1))
2732 END DO
2733 END DO
2734
2735 IF (ispin == 1) THEN
2736 q_a_matrix(:, :) = matmul(q_a_matrix, transpose(q_eigvec))
2737 ELSE
2738 q_b_matrix(:, :) = matmul(q_b_matrix, transpose(q_eigvec))
2739 END IF
2740
2741 DEALLOCATE (orb_occ)
2742 DEALLOCATE (q_eigval)
2743 DEALLOCATE (q_eigvec)
2744 END IF
2745 END IF
2746
2747 trq = 0.0_dp
2748 trq2 = 0.0_dp
2749 trqxq_hund_j = 0.0_dp
2750
2751 IF (ispin == 1) THEN
2752 DO i = 1, n
2753 trq = trq + q_a_matrix(i, i)
2754 DO j = 1, n
2755 trq2 = trq2 + q_a_matrix(i, j)*q_a_matrix(j, i)
2756 trqxq_hund_j = trqxq_hund_j + q_a_matrix(i, j)*q_b_matrix(j, i)
2757
2758 IF (i == j) THEN
2759 v_matrix(i, i) = u_minus_j*(0.5_dp - fspin*q_a_matrix(i, i))
2760 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 dft_control%mtlr_ikind == ikind .AND. &
2763 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
2764 v_matrix(i, i) = v_matrix(i, i) + perturbation_strength
2765 END IF
2766 ELSE
2767 v_matrix(i, j) = -u_minus_j*fspin*q_a_matrix(j, i)
2768 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 DO i = 1, n
2774 trq = trq + q_b_matrix(i, i)
2775 DO j = 1, n
2776 trq2 = trq2 + q_b_matrix(i, j)*q_b_matrix(j, i)
2777 trqxq_hund_j = trqxq_hund_j + q_b_matrix(i, j)*q_a_matrix(j, i)
2778
2779 IF (i == j) THEN
2780 v_matrix(i, i) = u_minus_j*(0.5_dp - fspin*q_b_matrix(i, i))
2781 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 dft_control%mtlr_ikind == ikind .AND. &
2784 qs_kind_set(ikind)%dft_plus_u%lr_atom == atom_a) THEN
2785 v_matrix(i, i) = v_matrix(i, i) + perturbation_strength
2786 END IF
2787 ELSE
2788 v_matrix(i, j) = -u_minus_j*fspin*q_b_matrix(j, i)
2789 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 v_matrix(:, :) = matmul(s_inv_matrix, v_matrix)
2796
2797 trq = fspin*trq
2798 trq2 = fspin*fspin*trq2
2799 trqxq_hund_j = fspin*fspin*trqxq_hund_j
2800
2801 energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
2802 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 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 found=found)
2813 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 found=found)
2820 cpassert(ASSOCIATED(a_block))
2821
2822 isb = 0
2823 DO iset = 1, nset
2824 DO ishell = 1, nshell(iset)
2825 IF (l(ishell, iset) /= lu) cycle
2826 isb = isb + 1
2827 i = 0
2828 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
2829 i = i + 1
2830 jsb = 0
2831 DO jset = 1, nset
2832 DO jshell = 1, nshell(jset)
2833 IF (l(jshell, jset) /= lu) cycle
2834 jsb = jsb + 1
2835 j = 0
2836 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
2837 j = j + 1
2838 v_block(isgf, jsgf) = v_matrix(i, j)*a(isb)*a(jsb)
2839 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 IF (should_output .AND. (print_level > low_print_level)) THEN
2854 IF (ispin == 1) THEN
2855 CALL para_env%sum(q_a_matrix)
2856 ELSE
2857 CALL para_env%sum(q_b_matrix)
2858 END IF
2859
2860 IF (output_unit > 0) THEN
2861 ALLOCATE (q_work(nsb, nsbsize))
2862 q_work(:, :) = 0.0_dp
2863
2864 DO isb = 1, nsb
2865 j = 0
2866 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
2867 j = j + 1
2868 IF (ispin == 1) THEN
2869 q_work(isb, j) = q_a_matrix(i, i)
2870 ELSE
2871 q_work(isb, j) = q_b_matrix(i, i)
2872 END IF
2873 END DO
2874 END DO
2875
2876 DO isb = 1, nsb
2877 WRITE (unit=output_unit, fmt="(T3,I6,2X,I6,2X,10F8.3)") &
2878 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 "Total", (sum(q_work(:, i)), i=1, nsbsize), sum(q_work)
2883 WRITE (unit=output_unit, fmt="(A)") ""
2884
2885 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 IF (ALLOCATED(q_a_matrix)) DEALLOCATE (q_a_matrix)
2965 IF (ALLOCATED(q_b_matrix)) DEALLOCATE (q_b_matrix)
2966 IF (ALLOCATED(s_matrix)) DEALLOCATE (s_matrix)
2967 IF (ALLOCATED(s_inv_matrix)) DEALLOCATE (s_inv_matrix)
2968 IF (ALLOCATED(v_matrix)) DEALLOCATE (v_matrix)
2969
2970 DEALLOCATE (a)
2971
2972 IF (dft_control%mtlr_dft_with_perturbation .AND. &
2973 dft_control%mtlr_ikind == ikind) THEN
2974 dft_control%vhxc(ispin) = dft_control%vhxc(ispin)/real(n, dp)
2975 CALL para_env%sum(dft_control%vhxc(ispin))
2976 CALL para_env%sum(dft_control%trq(ispin))
2977 DEALLOCATE (vhxc_matrix)
2978 END IF
2979
2980 END DO
2981
2982 ! Add V(i,j)[U] to V(i,j)[DFT]
2983
2984 IF (ASSOCIATED(sm_h)) THEN
2985 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, sm_v, 0.0_dp, sm_mid)
2986 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 CALL dbcsr_deallocate_matrix(sm_mid)
2992
2993 IF (PRESENT(matrix_w)) THEN
2994
2995 CALL get_qs_env(qs_env=qs_env, force=force)
2996 DO ikind = 1, nkind
2997 force(ikind)%tensorial_u(:, :) = 0.0_dp
2998 END DO
2999 IF (nspin == 2) THEN
3000 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_b)
3001 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_s, matrix_p(2)%matrix, 0.0_dp, sm_c)
3002 ELSE
3003 CALL dbcsr_multiply("N", "N", 0.5_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_b)
3004 CALL dbcsr_multiply("N", "N", 0.5_dp, sm_s, matrix_p(1)%matrix, 0.0_dp, sm_c)
3005 END IF
3006 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_b, sm_s, 0.0_dp, sm_i)
3007 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_c, sm_s, 0.0_dp, sm_j)
3008
3009 DO jj = 1, 3
3010 CALL dbcsr_desymmetrize(matrixkp_s(jj + 1, 1)%matrix, sm_k)
3011
3012 CALL dbcsr_set(sm_f, 0.0_dp)
3013
3014 CALL dbcsr_iterator_start(iter, sm_k)
3015 DO WHILE (dbcsr_iterator_blocks_left(iter))
3016 NULLIFY (tmp_block)
3017 CALL dbcsr_iterator_next_block(iter, jatom, katom, tmp_block)
3018 IF (jatom == katom) cycle
3019 CALL dbcsr_put_block(matrix=sm_f, row=jatom, col=katom, block=tmp_block)
3020 END DO
3021 CALL dbcsr_iterator_stop(iter)
3022 CALL dbcsr_finalize(sm_f)
3023
3024 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_b, sm_f, 0.0_dp, sm_d)
3025 CALL dbcsr_multiply("N", "N", 1.0_dp, sm_c, sm_f, 0.0_dp, sm_e)
3026
3027 ALLOCATE (a_matrix(nkind, nkind))
3028 ALLOCATE (d_matrix(nkind, nkind))
3029 ALLOCATE (e_matrix(nkind, nkind))
3030 ALLOCATE (i_matrix(nkind, nkind))
3031 ALLOCATE (j_matrix(nkind, nkind))
3032 ALLOCATE (dd_matrix(nkind, nkind))
3033 ALLOCATE (ee_matrix(nkind, nkind))
3034 ALLOCATE (ff_matrix(nkind, nkind))
3035
3036 DO ikind = 1, nkind
3037
3038 CALL get_atomic_kind(atomic_kind_set(ikind), &
3039 atom_list=atom_list, &
3040 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 dft_plus_u_atom=dft_plus_u_atom)
3049
3050 IF (.NOT. dft_plus_u_atom) cycle
3051 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
3052 IF (lu < 0) cycle
3053 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 nshell=nshell)
3061
3062 nsb = 0
3063 DO iset = 1, nset
3064 DO ishell = 1, nshell(iset)
3065 IF (l(ishell, iset) == lu) nsb = nsb + 1
3066 END DO
3067 END DO
3068
3069 nsbsize = (2*lu + 1)
3070 n = nsb*nsbsize
3071
3072 ALLOCATE (a_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3073 ALLOCATE (d_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3074 ALLOCATE (e_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3075 ALLOCATE (i_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3076 ALLOCATE (j_matrix(ikind, ikind)%matrix(n, n, natom_of_kind))
3077
3078 a_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3079 d_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3080 e_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3081 i_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3082 j_matrix(ikind, ikind)%matrix(:, :, :) = 0.0_dp
3083
3084 DO iatom = 1, natom_of_kind
3085
3086 atom_aa = atom_list(iatom)
3087
3088 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 found=found)
3094
3095 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 found=found)
3101
3102 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 found=found)
3108
3109 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 found=found)
3115
3116 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 found=found)
3122
3123 IF (ASSOCIATED(a_block) .AND. ASSOCIATED(d_block) .AND. ASSOCIATED(e_block) &
3124 .AND. ASSOCIATED(i_block) .AND. ASSOCIATED(j_block)) THEN
3125
3126 i = 0
3127 DO iset = 1, nset
3128 DO ishell = 1, nshell(iset)
3129 IF (l(ishell, iset) /= lu) cycle
3130 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3131 i = i + 1
3132 j = 0
3133 DO jset = 1, nset
3134 DO jshell = 1, nshell(jset)
3135 IF (l(jshell, jset) /= lu) cycle
3136 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
3137 j = j + 1
3138 a_matrix(ikind, ikind)%matrix(i, j, iatom) = a_block(isgf, jsgf)
3139 d_matrix(ikind, ikind)%matrix(i, j, iatom) = -d_block(isgf, jsgf)
3140 e_matrix(ikind, ikind)%matrix(i, j, iatom) = -e_block(isgf, jsgf)
3141 i_matrix(ikind, ikind)%matrix(i, j, iatom) = i_block(isgf, jsgf)
3142 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 DO jkind = 1, nkind
3153
3154 CALL get_atomic_kind(atomic_kind_set(jkind), &
3155 atom_list=atom_list2, &
3156 natom=matom_of_kind)
3157
3158 CALL get_qs_kind(qs_kind_set(jkind), nsgf=nsgf_kind2)
3159
3160 IF (.NOT. ALLOCATED(dd_matrix(ikind, jkind)%matrix)) THEN
3161 ALLOCATE (dd_matrix(ikind, jkind)%matrix(n, nsgf_kind2, natom_of_kind, matom_of_kind))
3162 dd_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3163 END IF
3164
3165 IF (.NOT. ALLOCATED(ee_matrix(ikind, jkind)%matrix)) THEN
3166 ALLOCATE (ee_matrix(ikind, jkind)%matrix(n, nsgf_kind2, natom_of_kind, matom_of_kind))
3167 ee_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3168 END IF
3169
3170 IF (.NOT. ALLOCATED(ff_matrix(ikind, jkind)%matrix)) THEN
3171 ALLOCATE (ff_matrix(ikind, jkind)%matrix(nsgf_kind2, n, natom_of_kind, matom_of_kind))
3172 ff_matrix(ikind, jkind)%matrix(:, :, :, :) = 0.0_dp
3173 END IF
3174
3175 DO jatom = 1, matom_of_kind
3176
3177 atom_a = atom_list2(jatom)
3178
3179 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 found=found)
3185
3186 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 found=found)
3192
3193 IF (ASSOCIATED(d_block) .AND. ASSOCIATED(e_block)) THEN
3194 i = 0
3195 DO iset = 1, nset
3196 DO ishell = 1, nshell(iset)
3197 IF (l(ishell, iset) /= lu) cycle
3198 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3199 i = i + 1
3200 dd_matrix(ikind, jkind)%matrix(i, :, iatom, jatom) = d_block(isgf, :)
3201 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 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 found=found)
3213
3214 IF (ASSOCIATED(f_block)) THEN
3215 i = 0
3216 DO iset = 1, nset
3217 DO ishell = 1, nshell(iset)
3218 IF (l(ishell, iset) /= lu) cycle
3219 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
3220 i = i + 1
3221 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 DO ikind = 1, nkind
3233
3234 CALL get_atomic_kind(atomic_kind_set(ikind), &
3235 atom_list=atom_list, &
3236 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 dft_plus_u_atom=dft_plus_u_atom)
3245
3246 IF (.NOT. dft_plus_u_atom) cycle
3247 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
3248 IF (lu < 0) cycle
3249 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 nshell=nshell)
3257
3258 CALL para_env%sum(a_matrix(ikind, ikind)%matrix(:, :, :))
3259 CALL para_env%sum(d_matrix(ikind, ikind)%matrix(:, :, :))
3260 CALL para_env%sum(e_matrix(ikind, ikind)%matrix(:, :, :))
3261 CALL para_env%sum(i_matrix(ikind, ikind)%matrix(:, :, :))
3262 CALL para_env%sum(j_matrix(ikind, ikind)%matrix(:, :, :))
3263
3264 DO jkind = 1, nkind
3265 CALL para_env%sum(dd_matrix(ikind, jkind)%matrix(:, :, :, :))
3266 CALL para_env%sum(ee_matrix(ikind, jkind)%matrix(:, :, :, :))
3267 CALL para_env%sum(ff_matrix(ikind, jkind)%matrix(:, :, :, :))
3268 END DO
3269
3270 END DO
3271
3272 ! Force calculation
3273
3274 DO ikind = 1, nkind
3275
3276 CALL get_atomic_kind(atomic_kind_set(ikind), &
3277 atom_list=atom_list, &
3278 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 dft_plus_u_atom=dft_plus_u_atom)
3287
3288 IF (.NOT. dft_plus_u_atom) cycle
3289 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
3290 IF (lu < 0) cycle
3291 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 nshell=nshell)
3299
3300 nsb = 0
3301 DO iset = 1, nset
3302 DO ishell = 1, nshell(iset)
3303 IF (l(ishell, iset) == lu) nsb = nsb + 1
3304 END DO
3305 END DO
3306
3307 nsbsize = (2*lu + 1)
3308 n = nsb*nsbsize
3309
3310 ALLOCATE (tmp_matrix1(n, n))
3311 ALLOCATE (tmp_matrix2(n, n))
3312 ALLOCATE (tmp_matrix3(n, n))
3313
3314 DO iatom = 1, natom_of_kind
3315
3316 atom_aa = atom_list(iatom)
3317
3318 tmp_matrix1(:, :) = matmul(d_matrix(ikind, ikind)%matrix(:, :, iatom), &
3319 a_matrix(ikind, ikind)%matrix(:, :, iatom))
3320 tmp_matrix2(:, :) = matmul(e_matrix(ikind, ikind)%matrix(:, :, iatom), &
3321 a_matrix(ikind, ikind)%matrix(:, :, iatom))
3322
3323 i_matrix(ikind, ikind)%matrix(:, :, iatom) = matmul(i_matrix(ikind, ikind)%matrix(:, :, iatom), &
3324 a_matrix(ikind, ikind)%matrix(:, :, iatom))
3325 j_matrix(ikind, ikind)%matrix(:, :, iatom) = matmul(j_matrix(ikind, ikind)%matrix(:, :, iatom), &
3326 a_matrix(ikind, ikind)%matrix(:, :, iatom))
3327
3328 tmp_matrix3(:, :) = 0.5_dp*u_minus_j*tmp_matrix1 - &
3329 u_minus_j*matmul(tmp_matrix1, i_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3330 + 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 u_minus_j*matmul(tmp_matrix2, j_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3333 + hund_j*matmul(tmp_matrix1, j_matrix(ikind, ikind)%matrix(:, :, iatom))
3334
3335 DO i = 1, n
3336 force(ikind)%tensorial_u(jj, iatom) = force(ikind)%tensorial_u(jj, iatom) + tmp_matrix3(i, i)
3337 END DO
3338
3339 DO jkind = 1, nkind
3340
3341 CALL get_atomic_kind(atomic_kind_set(jkind), &
3342 atom_list=atom_list2, &
3343 natom=matom_of_kind)
3344
3345 CALL get_qs_kind(qs_kind_set(jkind), nsgf=nsgf_kind2)
3346
3347 DO jatom = 1, matom_of_kind
3348
3349 atom_a = atom_list2(jatom)
3350
3351 tmp_matrix1(:, :) = matmul(dd_matrix(ikind, jkind)%matrix(:, :, iatom, jatom), &
3352 ff_matrix(ikind, jkind)%matrix(:, :, iatom, jatom))
3353 tmp_matrix1(:, :) = matmul(tmp_matrix1, a_matrix(ikind, ikind)%matrix(:, :, iatom))
3354
3355 tmp_matrix2(:, :) = matmul(ee_matrix(ikind, jkind)%matrix(:, :, iatom, jatom), &
3356 ff_matrix(ikind, jkind)%matrix(:, :, iatom, jatom))
3357 tmp_matrix2(:, :) = matmul(tmp_matrix2, a_matrix(ikind, ikind)%matrix(:, :, iatom))
3358
3359 tmp_matrix3(:, :) = 0.5_dp*u_minus_j*tmp_matrix1 - &
3360 u_minus_j*matmul(tmp_matrix1, i_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3361 + 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 u_minus_j*matmul(tmp_matrix2, j_matrix(ikind, ikind)%matrix(:, :, iatom)) &
3364 + hund_j*matmul(tmp_matrix1, j_matrix(ikind, ikind)%matrix(:, :, iatom))
3365
3366 IF (atom_a /= atom_aa) THEN
3367 DO i = 1, n
3368 force(jkind)%tensorial_u(jj, jatom) = force(jkind)%tensorial_u(jj, jatom) &
3369 + tmp_matrix3(i, i)
3370 END DO
3371 END IF
3372
3373 END DO
3374 END DO
3375 END DO
3376
3377 DEALLOCATE (tmp_matrix1)
3378 DEALLOCATE (tmp_matrix2)
3379 DEALLOCATE (tmp_matrix3)
3380
3381 END DO
3382
3383 DO ikind = 1, nkind
3384
3385 IF (ALLOCATED(a_matrix(ikind, ikind)%matrix)) DEALLOCATE (a_matrix(ikind, ikind)%matrix)
3386 IF (ALLOCATED(d_matrix(ikind, ikind)%matrix)) DEALLOCATE (d_matrix(ikind, ikind)%matrix)
3387 IF (ALLOCATED(e_matrix(ikind, ikind)%matrix)) DEALLOCATE (e_matrix(ikind, ikind)%matrix)
3388 IF (ALLOCATED(i_matrix(ikind, ikind)%matrix)) DEALLOCATE (i_matrix(ikind, ikind)%matrix)
3389 IF (ALLOCATED(j_matrix(ikind, ikind)%matrix)) DEALLOCATE (j_matrix(ikind, ikind)%matrix)
3390
3391 DO jkind = 1, nkind
3392 IF (ALLOCATED(dd_matrix(ikind, jkind)%matrix)) DEALLOCATE (dd_matrix(ikind, jkind)%matrix)
3393 IF (ALLOCATED(ee_matrix(ikind, jkind)%matrix)) DEALLOCATE (ee_matrix(ikind, jkind)%matrix)
3394 IF (ALLOCATED(ff_matrix(ikind, jkind)%matrix)) DEALLOCATE (ff_matrix(ikind, jkind)%matrix)
3395 END DO
3396
3397 END DO
3398
3399 DEALLOCATE (a_matrix)
3400 DEALLOCATE (d_matrix)
3401 DEALLOCATE (e_matrix)
3402 DEALLOCATE (i_matrix)
3403 DEALLOCATE (j_matrix)
3404 DEALLOCATE (dd_matrix)
3405 DEALLOCATE (ee_matrix)
3406 DEALLOCATE (ff_matrix)
3407
3408 END DO
3409
3410 DO ikind = 1, nkind
3411 force(ikind)%tensorial_u(:, :) = force(ikind)%tensorial_u(:, :)*2.0_dp
3412 END DO
3413
3414 CALL date_and_time(values=timevalues)
3415
3416 END IF
3417
3418 ! Collect the energy contributions from all processes
3419
3420 CALL para_env%sum(energy%dft_plus_u)
3421
3422 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 "if this warning persists or try a different method!")
3427 END IF
3428
3429 ! Release local sparse matrices
3430
3431 CALL dbcsr_deallocate_matrix(sm_q_a)
3432 CALL dbcsr_deallocate_matrix(sm_q_b)
3433 CALL dbcsr_deallocate_matrix(sm_v)
3434 CALL dbcsr_deallocate_matrix(sm_a)
3435
3436 IF (PRESENT(matrix_w)) THEN
3437 CALL dbcsr_deallocate_matrix(sm_b)
3438 CALL dbcsr_deallocate_matrix(sm_c)
3439 CALL dbcsr_deallocate_matrix(sm_d)
3440 CALL dbcsr_deallocate_matrix(sm_e)
3441 CALL dbcsr_deallocate_matrix(sm_f)
3442 CALL dbcsr_deallocate_matrix(sm_k)
3443 CALL dbcsr_deallocate_matrix(sm_i)
3444 CALL dbcsr_deallocate_matrix(sm_j)
3445 END IF
3446
3447 CALL timestop(handle)
3448
3449 END SUBROUTINE tensorial
3450
3451END MODULE dft_plus_u
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
mathlib::invert_matrix
Definition mathlib.F:74
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::dudarev1997
integer, save, public dudarev1997
Definition bibliography.F:36
bibliography::dudarev1998
integer, save, public dudarev1998
Definition bibliography.F:36
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1125
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_get_block_diag
subroutine, public dbcsr_get_block_diag(matrix, diag)
Copies the diagonal blocks of matrix into diag.
Definition cp_dbcsr_contrib.F:450
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:548
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Definition cp_dbcsr_operations.F:693
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_dbcsr_output
DBCSR output in CP2K.
Definition cp_dbcsr_output.F:17
cp_dbcsr_output::write_fm_with_basis_info
subroutine, public write_fm_with_basis_info(blacs_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, output_unit, omit_headers)
Print a spherical matrix of blacs type.
Definition cp_dbcsr_output.F:79
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers, cartesian_basis)
...
Definition cp_dbcsr_output.F:165
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_transpose
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
Definition cp_fm_basic_linalg.F:1532
cp_fm_basic_linalg::cp_fm_schur_product
subroutine, public cp_fm_schur_product(matrix_a, matrix_b, matrix_c)
computes the schur product of two matrices c_ij = a_ij * b_ij
Definition cp_fm_basic_linalg.F:724
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:170
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:245
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:750
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::low_print_level
integer, parameter, public low_print_level
Definition cp_output_handling.F:73
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
dft_plus_u
Add the DFT+U contribution to the Hamiltonian matrix.
Definition dft_plus_u.F:18
dft_plus_u::plus_u
subroutine, public plus_u(qs_env, matrix_h, matrix_w)
Add the DFT+U contribution to the Hamiltonian matrix. Wrapper routine for all "+U" methods.
Definition dft_plus_u.F:107
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::plus_u_lowdin
integer, parameter, public plus_u_lowdin
Definition input_constants.F:924
input_constants::plus_u_mulliken_charges
integer, parameter, public plus_u_mulliken_charges
Definition input_constants.F:924
input_constants::plus_u_mulliken
integer, parameter, public plus_u_mulliken
Definition input_constants.F:924
input_constants::plus_u_tensorial
integer, parameter, public plus_u_tensorial
Definition input_constants.F:924
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::lowdin_kp_trans
subroutine, public lowdin_kp_trans(kpoint, pmat_diag)
Calculate Lowdin transformation of density matrix S^1/2 P S^1/2 Integrate diagonal elements over k-po...
Definition kpoint_methods.F:1610
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::jacobi
subroutine, public jacobi(a, d, v)
Jacobi matrix diagonalization. The eigenvalues are returned in vector d and the eigenvectors are retu...
Definition mathlib.F:1605
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
mulliken
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition mulliken.F:13
orbital_symbols
orbital_symbols
Definition orbital_symbols.F:14
orbital_symbols::sgf_symbol
character(len=6) function, public sgf_symbol(n, l, m)
Build a spherical orbital symbol (orbital labels for printing).
Definition orbital_symbols.F:107
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:121
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_force_types
Definition qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:486
qs_kind_types::set_qs_kind
subroutine, public set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, covalent_radius, vdw_radius, lmax_rho0, zeff, no_optimize, dispersion, u_minus_j, hund_j, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
Definition qs_kind_types.F:3222
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:979
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:175
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:746
particle_types::particle_type
Definition particle_types.F:35
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:185
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:97