(git:01d6ff2)
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_deallocate_matrix, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
30 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
31 dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type
32 USE cp_dbcsr_contrib, ONLY: dbcsr_get_block_diag
33 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
34 copy_fm_to_dbcsr,&
35 cp_dbcsr_plus_fm_fm_t,&
36 cp_dbcsr_sm_fm_multiply
37 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix,&
38 write_fm_with_basis_info
39 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
40 cp_fm_schur_product,&
41 cp_fm_transpose
42 USE cp_fm_diag, ONLY: choose_eigv_solver
43 USE cp_fm_struct, ONLY: cp_fm_struct_type
44 USE cp_fm_types, ONLY: cp_fm_create,&
45 cp_fm_get_info,&
46 cp_fm_release,&
47 cp_fm_set_submatrix,&
48 cp_fm_type
49 USE cp_log_handling, ONLY: cp_get_default_logger,&
50 cp_logger_type
51 USE cp_output_handling, ONLY: cp_p_file,&
52 cp_print_key_finished_output,&
53 cp_print_key_should_output,&
54 cp_print_key_unit_nr,&
55 low_print_level
56 USE input_constants, ONLY: plus_u_lowdin,&
57 plus_u_mulliken,&
58 plus_u_mulliken_charges
59 USE input_section_types, ONLY: section_vals_type
60 USE kinds, ONLY: default_string_length,&
61 dp
62 USE kpoint_methods, ONLY: lowdin_kp_trans
63 USE kpoint_types, ONLY: kpoint_type
64 USE mathlib, ONLY: jacobi
65 USE message_passing, ONLY: mp_para_env_type
66 USE orbital_symbols, ONLY: sgf_symbol
67 USE parallel_gemm_api, ONLY: parallel_gemm
68 USE particle_methods, ONLY: get_particle_set
69 USE particle_types, ONLY: particle_type
70 USE physcon, ONLY: evolt
71 USE qs_energy_types, ONLY: qs_energy_type
72 USE qs_environment_types, ONLY: get_qs_env,&
73 qs_environment_type
74 USE qs_kind_types, ONLY: get_qs_kind,&
75 get_qs_kind_set,&
76 qs_kind_type,&
77 set_qs_kind
78 USE qs_rho_types, ONLY: qs_rho_get,&
79 qs_rho_type
80 USE qs_scf_types, ONLY: qs_scf_env_type
81#include "./base/base_uses.f90"
82
83 IMPLICIT NONE
84
85 PRIVATE
86
87 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dft_plus_u'
88
89 PUBLIC :: plus_u
90
91CONTAINS
92! **************************************************************************************************
93!> \brief Add the DFT+U contribution to the Hamiltonian matrix.\n
94!> Wrapper routine for all "+U" methods
95!> \param[in] qs_env Quickstep environment
96!> \param[in,out] matrix_h Hamiltonian matrices for each spin
97!> \param[in,out] matrix_w Energy weighted density matrices for each spin
98!> \date 14.01.2008
99!> \author Matthias Krack (MK)
100!> \version 1.0
101! **************************************************************************************************
102 SUBROUTINE plus_u(qs_env, matrix_h, matrix_w)
103
104 TYPE(qs_environment_type), POINTER :: qs_env
105 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
106 POINTER :: matrix_h, matrix_w
107
108 CHARACTER(LEN=*), PARAMETER :: routinen = 'plus_u'
109
110 INTEGER :: handle, output_unit, print_level
111 LOGICAL :: orthonormal_basis, should_output
112 TYPE(cp_logger_type), POINTER :: logger
113 TYPE(dft_control_type), POINTER :: dft_control
114 TYPE(section_vals_type), POINTER :: input
115
116 CALL timeset(routinen, handle)
117
118 cpassert(ASSOCIATED(qs_env))
119
120 NULLIFY (input, dft_control)
121
122 logger => cp_get_default_logger()
123
124 CALL get_qs_env(qs_env=qs_env, &
125 input=input, &
126 dft_control=dft_control)
127
128 CALL cite_reference(dudarev1997)
129 CALL cite_reference(dudarev1998)
130
131 ! Later we could save here some time, if the method in use has this property
132 ! which then has to be figured out here
133
134 orthonormal_basis = .false.
135
136 ! Setup print control
137
138 print_level = logger%iter_info%print_level
139 should_output = (btest(cp_print_key_should_output(logger%iter_info, input, &
140 "DFT%PRINT%PLUS_U"), cp_p_file) .AND. &
141 (.NOT. PRESENT(matrix_w)))
142 output_unit = cp_print_key_unit_nr(logger, input, "DFT%PRINT%PLUS_U", &
143 extension=".plus_u", &
144 ignore_should_output=should_output, &
145 log_filename=.false.)
146
147 ! Select DFT+U method
148
149 SELECT CASE (dft_control%plus_u_method_id)
150 CASE (plus_u_lowdin)
151 IF (orthonormal_basis) THEN
152 ! For an orthonormal basis the Lowdin method and the Mulliken method
153 ! are equivalent
154 CALL mulliken(qs_env, orthonormal_basis, matrix_h, &
155 should_output, output_unit, print_level)
156 ELSE
157 CALL lowdin(qs_env, matrix_h, matrix_w, &
158 should_output, output_unit, print_level)
159 END IF
160 CASE (plus_u_mulliken)
161 CALL mulliken(qs_env, orthonormal_basis, matrix_h, &
162 should_output, output_unit, print_level)
163 CASE (plus_u_mulliken_charges)
164 CALL mulliken_charges(qs_env, orthonormal_basis, matrix_h, matrix_w, &
165 should_output, output_unit, print_level)
166 CASE DEFAULT
167 cpabort("Invalid DFT+U method requested")
168 END SELECT
169
170 CALL cp_print_key_finished_output(output_unit, logger, input, "DFT%PRINT%PLUS_U", &
171 ignore_should_output=should_output)
172
173 CALL timestop(handle)
174
175 END SUBROUTINE plus_u
176
177! **************************************************************************************************
178!> \brief Add a DFT+U contribution to the Hamiltonian matrix\n
179!> using a method based on Lowdin charges
180!> \f[Q = S^{1/2} P S^{1/2}\f]
181!> where \b P and \b S are the density and the
182!> overlap matrix, respectively.
183!> \param[in] qs_env Quickstep environment
184!> \param[in,out] matrix_h Hamiltonian matrices for each spin
185!> \param[in,out] matrix_w Energy weighted density matrices for each spin
186!> \param should_output ...
187!> \param output_unit ...
188!> \param print_level ...
189!> \date 02.07.2008
190!> \par
191!> \f{eqnarray*}{
192!> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}\\\
193!> & = & E^{\rm DFT} + \frac{1}{2}(U - J)\sum_\mu (q_\mu - q_\mu^2)\\[1ex]
194!> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
195!> & = & \frac{\partial E^{\rm DFT}}
196!> {\partial P_{\mu\nu}} +
197!> \frac{\partial E^{\rm U}}
198!> {\partial P_{\mu\nu}}\\\
199!> & = & H_{\mu\nu} +
200!> \frac{\partial E^{\rm U}}{\partial q_\mu}
201!> \frac{\partial q_\mu}{\partial P_{\mu\nu}}\\\
202!> \f}
203!> \author Matthias Krack (MK)
204!> \version 1.0
205! **************************************************************************************************
206 SUBROUTINE lowdin(qs_env, matrix_h, matrix_w, should_output, output_unit, &
207 print_level)
208
209 TYPE(qs_environment_type), POINTER :: qs_env
210 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
211 POINTER :: matrix_h, matrix_w
212 LOGICAL, INTENT(IN) :: should_output
213 INTEGER, INTENT(IN) :: output_unit, print_level
214
215 CHARACTER(LEN=*), PARAMETER :: routinen = 'lowdin'
216
217 CHARACTER(LEN=10) :: spin_info
218 CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: symbol
219 CHARACTER(LEN=default_string_length) :: atomic_kind_name
220 INTEGER :: atom_a, handle, i, i0, iatom, ikind, iorb, isb, iset, isgf, ishell, ispin, j, &
221 jsb, jset, jsgf, jshell, lu, m, max_scf, n, natom, natom_of_kind, nimg, nkind, norb, nsb, &
222 nsbsize, nset, nsgf, nsgf_kind, nspin
223 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
224 INTEGER, DIMENSION(1) :: iloc
225 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals
226 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
227 LOGICAL :: debug, dft_plus_u_atom, do_kpoints, &
228 found, just_energy, smear
229 LOGICAL, ALLOCATABLE, DIMENSION(:) :: orb_occ
230 REAL(kind=dp) :: eps_scf, eps_u_ramping, fspin, occ, sij, &
231 trq, trq2, u_minus_j, &
232 u_minus_j_target, u_ramping
233 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigval, q_eigval
234 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: q_eigvec, q_matrix, q_work, slam
235 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
236 POINTER :: local_data
237 REAL(kind=dp), DIMENSION(:, :), POINTER :: q_block, v_block
238 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
239 TYPE(cp_fm_struct_type), POINTER :: fmstruct
240 TYPE(cp_fm_type) :: fm_sev, fm_work1, fm_work2, slambda
241 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fm_wmat
242 TYPE(cp_fm_type), POINTER :: fm_s_half
243 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
244 TYPE(dbcsr_type) :: sm_q, sm_v
245 TYPE(dbcsr_type), POINTER :: sm_h, sm_p, sm_s, sm_w
246 TYPE(dft_control_type), POINTER :: dft_control
247 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
248 TYPE(kpoint_type), POINTER :: kpoints
249 TYPE(mp_para_env_type), POINTER :: para_env
250 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
251 TYPE(qs_energy_type), POINTER :: energy
252 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
253 TYPE(qs_rho_type), POINTER :: rho
254 TYPE(qs_scf_env_type), POINTER :: scf_env
255
256 CALL timeset(routinen, handle)
257
258 debug = .false. ! Set to .TRUE. to print debug information
259
260 NULLIFY (atom_list)
261 NULLIFY (atomic_kind_set)
262 NULLIFY (qs_kind_set)
263 NULLIFY (dft_control)
264 NULLIFY (energy)
265 NULLIFY (first_sgf)
266 NULLIFY (fm_s_half)
267 NULLIFY (matrix_p)
268 NULLIFY (matrix_s)
269 NULLIFY (l)
270 NULLIFY (last_sgf)
271 NULLIFY (nshell)
272 NULLIFY (orb_basis_set)
273 NULLIFY (orbitals)
274 NULLIFY (particle_set)
275 NULLIFY (q_block)
276 NULLIFY (rho)
277 NULLIFY (scf_env)
278 NULLIFY (sm_h)
279 NULLIFY (sm_p)
280 NULLIFY (sm_s)
281 NULLIFY (sm_w)
282 NULLIFY (v_block)
283 NULLIFY (para_env)
284
285 smear = .false.
286 max_scf = -1
287 eps_scf = 1.0e30_dp
288
289 CALL get_qs_env(qs_env=qs_env, &
290 atomic_kind_set=atomic_kind_set, &
291 qs_kind_set=qs_kind_set, &
292 dft_control=dft_control, &
293 do_kpoints=do_kpoints, &
294 kpoints=kpoints, &
295 energy=energy, &
296 matrix_s_kp=matrix_s, &
297 particle_set=particle_set, &
298 rho=rho, &
299 scf_env=scf_env, &
300 para_env=para_env)
301
302 cpassert(ASSOCIATED(atomic_kind_set))
303 cpassert(ASSOCIATED(dft_control))
304 cpassert(ASSOCIATED(energy))
305 cpassert(ASSOCIATED(matrix_s))
306 cpassert(ASSOCIATED(particle_set))
307 cpassert(ASSOCIATED(rho))
308
309 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
310
311 energy%dft_plus_u = 0.0_dp
312
313 nspin = dft_control%nspins
314 nimg = dft_control%nimages
315
316 IF (nspin == 2) THEN
317 fspin = 1.0_dp
318 ELSE
319 fspin = 0.5_dp
320 END IF
321
322 ! Get the total number of atoms, contracted spherical Gaussian basis
323 ! functions, and atomic kinds
324
325 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
326 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
327
328 nkind = SIZE(atomic_kind_set)
329
330 ALLOCATE (first_sgf_atom(natom))
331 first_sgf_atom(:) = 0
332
333 CALL get_particle_set(particle_set, qs_kind_set, &
334 first_sgf=first_sgf_atom)
335
336 IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
337 just_energy = .false.
338 ELSE
339 just_energy = .true.
340 END IF
341
342 IF (do_kpoints) THEN
343 fm_wmat => scf_env%scf_work1
344 fmstruct => fm_wmat(1)%matrix_struct
345 ELSE
346 ! Retrieve S^(1/2) from the SCF environment
347 fm_s_half => scf_env%s_half
348 cpassert(ASSOCIATED(fm_s_half))
349 ! work matrices
350 CALL cp_fm_get_info(fm_s_half, matrix_struct=fmstruct)
351 END IF
352 CALL cp_fm_create(matrix=fm_work1, matrix_struct=fmstruct, &
353 name="FULL WORK MATRIX 1")
354 CALL cp_fm_create(matrix=fm_work2, matrix_struct=fmstruct, &
355 name="FULL WORK MATRIX 2")
356
357 ! Calculate S eigenvectors and Lambda matrix for forces
358 ! See sTDA forces (get_lowdin_mo_coefficients in qs_tddfpt2_stda_utils
359 ! A. Hehn et al JCTC 2022, 18, 4186
360 IF (PRESENT(matrix_w)) THEN
361 CALL cp_fm_create(matrix=fm_sev, matrix_struct=fmstruct)
362 CALL cp_fm_create(matrix=slambda, matrix_struct=fmstruct)
363 ALLOCATE (eigval(nsgf), slam(nsgf, 1))
364 sm_s => matrix_s(1, 1)%matrix
365 CALL copy_dbcsr_to_fm(sm_s, fm_work1)
366 CALL choose_eigv_solver(fm_work1, fm_sev, eigval)
367 !
368 DO i = 1, nsgf
369 IF (eigval(i) > 0._dp) THEN
370 slam(i, 1) = sqrt(eigval(i))
371 ELSE
372 cpabort("S matrix not positive definit")
373 END IF
374 END DO
375 DO i = 1, nsgf
376 CALL cp_fm_set_submatrix(slambda, slam, 1, i, nsgf, 1, 1.0_dp, 0.0_dp)
377 END DO
378 DO i = 1, nsgf
379 CALL cp_fm_set_submatrix(slambda, slam, i, 1, 1, nsgf, 1.0_dp, 1.0_dp, .true.)
380 END DO
381 CALL cp_fm_get_info(slambda, local_data=local_data)
382 DO i = 1, SIZE(local_data, 2)
383 DO j = 1, SIZE(local_data, 1)
384 sij = local_data(j, i)
385 IF (sij > 0.0_dp) sij = 1.0_dp/sij
386 local_data(j, i) = sij
387 END DO
388 END DO
389 DEALLOCATE (eigval, slam)
390 END IF
391
392 ! Calculate S^(1/2)*P*S^(1/2) as a full matrix (Lowdin)
393 IF (do_kpoints) THEN
394 ! stored in kp%pmat
395 DO ispin = 1, nspin
396 CALL lowdin_kp_trans(matrix_p, kpoints, ispin, fm_wmat)
397 END DO
398 END IF
399
400 ! Create local block diagonal matrices
401 sm_s => matrix_s(1, 1)%matrix
402 CALL dbcsr_get_block_diag(sm_s, sm_q)
403 CALL dbcsr_get_block_diag(sm_s, sm_v)
404
405 ! Loop over all spins
406 DO ispin = 1, nspin
407
408 CALL dbcsr_set(sm_q, 0.0_dp)
409 CALL dbcsr_set(sm_v, 0.0_dp)
410
411 ! Calculate S^(1/2)*P*S^(1/2) as a full matrix (Lowdin)
412 IF (do_kpoints) THEN
413 ! assemble results in sm_q from kp%pmat
414 cpabort("Lowdin option with k-points NYA in DFT+U")
415 ELSE
416 sm_p => matrix_p(ispin, 1)%matrix
417 CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_s_half, fm_work1, nsgf)
418 CALL parallel_gemm(transa="N", &
419 transb="N", &
420 m=nsgf, &
421 n=nsgf, &
422 k=nsgf, &
423 alpha=1.0_dp, &
424 matrix_a=fm_s_half, &
425 matrix_b=fm_work1, &
426 beta=0.0_dp, &
427 matrix_c=fm_work2)
428 IF (debug) THEN
429 CALL cp_dbcsr_write_sparse_matrix(sm_p, 4, 6, qs_env, para_env, &
430 output_unit=output_unit)
431 CALL write_fm_with_basis_info(fm_s_half, 4, 6, qs_env, para_env, &
432 output_unit=output_unit)
433 CALL write_fm_with_basis_info(fm_work2, 4, 6, qs_env, para_env, &
434 output_unit=output_unit)
435 END IF ! debug
436 ! Copy occupation matrix to sparse matrix format, finally we are only
437 ! interested in the diagonal (atomic) blocks, i.e. the previous full
438 ! matrix product is not the most efficient choice, anyway.
439 CALL copy_fm_to_dbcsr(fm_work2, sm_q, keep_sparsity=.true.)
440 END IF
441
442 ! E[DFT+U] = E[DFT] + E[U]
443 ! = E[DFT] + (U - J)*(Tr(q) - Tr(q*q))/2
444
445 ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U]
446 ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j)
447 ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j))
448
449 ! Loop over all atomic kinds
450 DO ikind = 1, nkind
451
452 ! Load the required atomic kind data
453 CALL get_atomic_kind(atomic_kind_set(ikind), &
454 atom_list=atom_list, &
455 name=atomic_kind_name, &
456 natom=natom_of_kind)
457
458 CALL get_qs_kind(qs_kind_set(ikind), &
459 dft_plus_u_atom=dft_plus_u_atom, &
460 l_of_dft_plus_u=lu, &
461 nsgf=nsgf_kind, &
462 basis_set=orb_basis_set, &
463 u_minus_j=u_minus_j, &
464 u_minus_j_target=u_minus_j_target, &
465 u_ramping=u_ramping, &
466 eps_u_ramping=eps_u_ramping, &
467 orbitals=orbitals, &
468 eps_scf=eps_scf, &
469 max_scf=max_scf, &
470 smear=smear)
471
472 ! Check, if the atoms of this atomic kind need a DFT+U correction
473 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
474 IF (.NOT. dft_plus_u_atom) cycle
475 IF (lu < 0) cycle
476
477 ! Apply U ramping if requested
478 IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
479 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
480 u_minus_j = min(u_minus_j + u_ramping, u_minus_j_target)
481 CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
482 END IF
483 IF (should_output .AND. (output_unit > 0)) THEN
484 WRITE (unit=output_unit, fmt="(T3,A,3X,A,F0.3,A)") &
485 "Kind name: "//trim(adjustl(atomic_kind_name)), &
486 "U(eff) = ", u_minus_j*evolt, " eV"
487 END IF
488 END IF
489
490 IF (u_minus_j == 0.0_dp) cycle
491
492 ! Load the required Gaussian basis set data
493 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
494 first_sgf=first_sgf, &
495 l=l, &
496 last_sgf=last_sgf, &
497 nset=nset, &
498 nshell=nshell)
499
500 ! Count the relevant shell blocks of this atomic kind
501 nsb = 0
502 DO iset = 1, nset
503 DO ishell = 1, nshell(iset)
504 IF (l(ishell, iset) == lu) nsb = nsb + 1
505 END DO
506 END DO
507
508 nsbsize = (2*lu + 1)
509 n = nsb*nsbsize
510
511 ALLOCATE (q_matrix(n, n))
512 q_matrix(:, :) = 0.0_dp
513
514 ! Print headline if requested
515 IF (should_output .AND. (print_level > low_print_level)) THEN
516 IF (output_unit > 0) THEN
517 ALLOCATE (symbol(nsbsize))
518 DO m = -lu, lu
519 symbol(lu + m + 1) = sgf_symbol(0, lu, m)
520 END DO
521 IF (nspin > 1) THEN
522 WRITE (unit=spin_info, fmt="(A8,I2)") " of spin", ispin
523 ELSE
524 spin_info = ""
525 END IF
526 WRITE (unit=output_unit, fmt="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
527 "DFT+U occupations"//trim(spin_info)//" for the atoms of atomic kind ", ikind, &
528 ": "//trim(atomic_kind_name), &
529 "Atom Shell ", (adjustr(symbol(i)), i=1, nsbsize), " Trace"
530 DEALLOCATE (symbol)
531 END IF
532 END IF
533
534 ! Loop over all atoms of the current atomic kind
535 DO iatom = 1, natom_of_kind
536 atom_a = atom_list(iatom)
537 q_matrix(:, :) = 0.0_dp
538
539 ! Get diagonal block
540 CALL dbcsr_get_block_p(matrix=sm_q, &
541 row=atom_a, &
542 col=atom_a, &
543 block=q_block, &
544 found=found)
545
546 IF (ASSOCIATED(q_block)) THEN
547 ! Calculate energy contribution to E(U)
548 i = 0
549 DO iset = 1, nset
550 DO ishell = 1, nshell(iset)
551 IF (l(ishell, iset) /= lu) cycle
552 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
553 i = i + 1
554 j = 0
555 DO jset = 1, nset
556 DO jshell = 1, nshell(jset)
557 IF (l(jshell, jset) /= lu) cycle
558 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
559 j = j + 1
560 IF (isgf == jsgf) q_matrix(i, j) = q_block(isgf, jsgf)
561 END DO ! next contracted spherical Gaussian function "jsgf"
562 END DO ! next shell "jshell"
563 END DO ! next shell set "jset"
564 END DO ! next contracted spherical Gaussian function "isgf"
565 END DO ! next shell "ishell"
566 END DO ! next shell set "iset"
567
568 ! Perform the requested manipulations of the (initial) orbital occupations
569 IF (ASSOCIATED(orbitals)) THEN
570 IF ((qs_env%scf_env%iter_delta >= eps_scf) .OR. &
571 ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
572 (qs_env%scf_env%iter_count <= max_scf))) THEN
573 ALLOCATE (orb_occ(nsbsize))
574 ALLOCATE (q_eigval(n))
575 q_eigval(:) = 0.0_dp
576 ALLOCATE (q_eigvec(n, n))
577 q_eigvec(:, :) = 0.0_dp
578 norb = SIZE(orbitals)
579 CALL jacobi(q_matrix, q_eigval, q_eigvec)
580 q_matrix(:, :) = 0.0_dp
581 DO isb = 1, nsb
582 trq = 0.0_dp
583 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
584 trq = trq + q_eigval(i)
585 END DO
586 IF (smear) THEN
587 occ = trq/real(norb, kind=dp)
588 ELSE
589 occ = 1.0_dp/fspin
590 END IF
591 orb_occ(:) = .false.
592 iloc = maxloc(q_eigvec(:, isb*nsbsize))
593 jsb = int((iloc(1) - 1)/nsbsize) + 1
594 i = 0
595 i0 = (jsb - 1)*nsbsize + 1
596 iorb = -1000
597 DO j = i0, jsb*nsbsize
598 i = i + 1
599 IF (i > norb) THEN
600 DO m = -lu, lu
601 IF (.NOT. orb_occ(lu + m + 1)) THEN
602 iorb = i0 + lu + m
603 orb_occ(lu + m + 1) = .true.
604 END IF
605 END DO
606 ELSE
607 iorb = i0 + lu + orbitals(i)
608 orb_occ(lu + orbitals(i) + 1) = .true.
609 END IF
610 cpassert(iorb /= -1000)
611 iloc = maxloc(q_eigvec(iorb, :))
612 q_eigval(iloc(1)) = min(occ, trq)
613 q_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1)) ! backtransform left
614 trq = trq - q_eigval(iloc(1))
615 END DO
616 END DO
617 q_matrix(:, :) = matmul(q_matrix, transpose(q_eigvec)) ! backtransform right
618 DEALLOCATE (orb_occ)
619 DEALLOCATE (q_eigval)
620 DEALLOCATE (q_eigvec)
621 END IF
622 END IF ! orbitals associated
623
624 trq = 0.0_dp
625 trq2 = 0.0_dp
626 DO i = 1, n
627 trq = trq + q_matrix(i, i)
628 DO j = 1, n
629 trq2 = trq2 + q_matrix(i, j)*q_matrix(j, i)
630 END DO
631 END DO
632 trq = fspin*trq
633 trq2 = fspin*fspin*trq2
634
635 ! Calculate energy contribution to E(U)
636 energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
637
638 ! Calculate potential V(U) = dE(U)/dq
639 IF (.NOT. just_energy) THEN
640 CALL dbcsr_get_block_p(matrix=sm_v, &
641 row=atom_a, &
642 col=atom_a, &
643 block=v_block, &
644 found=found)
645 cpassert(ASSOCIATED(v_block))
646
647 i = 0
648 DO iset = 1, nset
649 DO ishell = 1, nshell(iset)
650 IF (l(ishell, iset) /= lu) cycle
651 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
652 i = i + 1
653 j = 0
654 DO jset = 1, nset
655 DO jshell = 1, nshell(jset)
656 IF (l(jshell, jset) /= lu) cycle
657 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
658 j = j + 1
659 IF (isgf == jsgf) THEN
660 v_block(isgf, isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i, i))
661 ELSE
662 cpassert(abs(q_matrix(j, i)) < 1.0e-14_dp)
663 v_block(isgf, jsgf) = -u_minus_j*fspin*q_matrix(j, i)
664 END IF
665 END DO ! next contracted spherical Gaussian function "jsgf"
666 END DO ! next shell "jshell"
667 END DO ! next shell set "jset"
668 END DO ! next contracted spherical Gaussian function "isgf"
669 END DO ! next shell "ishell"
670 END DO ! next shell set "iset"
671 END IF ! not just energy
672
673 END IF ! q_block associated
674
675 ! Consider print requests
676 IF (should_output .AND. (print_level > low_print_level)) THEN
677 CALL para_env%sum(q_matrix)
678 IF (output_unit > 0) THEN
679 ALLOCATE (q_work(nsb, nsbsize))
680 q_work(:, :) = 0.0_dp
681 DO isb = 1, nsb
682 j = 0
683 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
684 j = j + 1
685 q_work(isb, j) = q_matrix(i, i)
686 END DO
687 END DO
688 DO isb = 1, nsb
689 WRITE (unit=output_unit, fmt="(T3,I6,2X,I6,2X,10F8.3)") &
690 atom_a, isb, q_work(isb, :), sum(q_work(isb, :))
691 END DO
692 WRITE (unit=output_unit, fmt="(T12,A,2X,10F8.3)") &
693 "Total", (sum(q_work(:, i)), i=1, nsbsize), sum(q_work)
694 WRITE (unit=output_unit, fmt="(A)") ""
695 DEALLOCATE (q_work)
696 IF (debug) THEN
697 ! Print the DFT+U occupation matrix
698 WRITE (unit=output_unit, fmt="(T9,70I10)") (i, i=1, n)
699 DO i = 1, n
700 WRITE (unit=output_unit, fmt="(T3,I6,70F10.6)") i, q_matrix(i, :)
701 END DO
702 ! Print the eigenvalues and eigenvectors of the occupation matrix
703 ALLOCATE (q_eigval(n))
704 q_eigval(:) = 0.0_dp
705 ALLOCATE (q_eigvec(n, n))
706 q_eigvec(:, :) = 0.0_dp
707 CALL jacobi(q_matrix, q_eigval, q_eigvec)
708 WRITE (unit=output_unit, fmt="(/,T9,70I10)") (i, i=1, n)
709 WRITE (unit=output_unit, fmt="(T9,71F10.6)") (q_eigval(i), i=1, n), &
710 sum(q_eigval(1:n))
711 DO i = 1, n
712 WRITE (unit=output_unit, fmt="(T3,I6,70F10.6)") i, q_eigvec(i, :)
713 END DO
714 DEALLOCATE (q_eigval)
715 DEALLOCATE (q_eigvec)
716 END IF ! debug
717 END IF
718 IF (debug) THEN
719 ! Print the full atomic occupation matrix block
720 ALLOCATE (q_work(nsgf_kind, nsgf_kind))
721 q_work(:, :) = 0.0_dp
722 IF (ASSOCIATED(q_block)) q_work(:, :) = q_block(:, :)
723 CALL para_env%sum(q_work)
724 IF (output_unit > 0) THEN
725 norb = SIZE(q_work, 1)
726 WRITE (unit=output_unit, fmt="(/,T9,200I10)") (i, i=1, norb)
727 DO i = 1, norb
728 WRITE (unit=output_unit, fmt="(T3,I6,200F10.6)") i, q_work(i, :)
729 END DO
730 ALLOCATE (q_eigval(norb))
731 q_eigval(:) = 0.0_dp
732 ALLOCATE (q_eigvec(norb, norb))
733 q_eigvec(:, :) = 0.0_dp
734 CALL jacobi(q_work, q_eigval, q_eigvec)
735 WRITE (unit=output_unit, fmt="(/,T9,200I10)") (i, i=1, norb)
736 WRITE (unit=output_unit, fmt="(T9,201F10.6)") (q_eigval(i), i=1, norb), &
737 sum(q_eigval(1:norb))
738 DO i = 1, norb
739 WRITE (unit=output_unit, fmt="(T3,I6,200F10.6)") i, q_eigvec(i, :)
740 END DO
741 DEALLOCATE (q_eigval)
742 DEALLOCATE (q_eigvec)
743 END IF
744 DEALLOCATE (q_work)
745 END IF ! debug
746 END IF ! should output
747
748 END DO ! next atom "iatom" of atomic kind "ikind"
749
750 IF (ALLOCATED(q_matrix)) THEN
751 DEALLOCATE (q_matrix)
752 END IF
753 END DO ! next atomic kind "ikind"
754
755 ! Add V(i,j)[U] to V(i,j)[DFT]
756 IF (PRESENT(matrix_h)) THEN
757 IF (do_kpoints) THEN
758 cpabort("Lowdin option with k-points NYA in DFT+U")
759 ELSE
760 sm_h => matrix_h(ispin, 1)%matrix
761 CALL cp_dbcsr_sm_fm_multiply(sm_v, fm_s_half, fm_work1, nsgf)
762 CALL cp_fm_transpose(fm_work1, fm_work2)
763 CALL cp_dbcsr_plus_fm_fm_t(sm_h, fm_s_half, fm_work2, nsgf)
764 END IF
765 END IF ! An update of the Hamiltonian matrix is requested
766
767 ! Calculate the contribution (non-Pulay part) to the derivatives
768 ! w.r.t. the nuclear positions
769 IF (PRESENT(matrix_w)) THEN
770
771 sm_p => matrix_p(ispin, 1)%matrix
772 sm_w => matrix_w(ispin, 1)%matrix
773
774 CALL cp_dbcsr_sm_fm_multiply(sm_v, fm_s_half, fm_work1, nsgf)
775 CALL cp_fm_transpose(fm_work1, fm_work2)
776 CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_work2, fm_work1, nsgf)
777 CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, 1.0_dp, fm_work1, fm_sev, 0.0_dp, fm_work2)
778 CALL parallel_gemm('T', 'N', nsgf, nsgf, nsgf, 1.0_dp, fm_sev, fm_work2, 0.0_dp, fm_work1)
779 CALL cp_fm_schur_product(fm_work1, slambda, fm_work2)
780 CALL cp_fm_transpose(fm_work2, fm_work1)
781 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_work1, matrix_b=fm_work2)
782 CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, -2.0_dp, fm_sev, fm_work2, 0.0_dp, fm_work1)
783 CALL cp_dbcsr_plus_fm_fm_t(sm_w, fm_work1, fm_sev, nsgf)
784
785 END IF ! W matrix update requested
786
787 END DO ! next spin "ispin"
788
789 IF (PRESENT(matrix_w)) THEN
790 CALL cp_fm_release(matrix=fm_sev)
791 CALL cp_fm_release(matrix=slambda)
792 END IF
793
794 ! Collect the energy contributions from all processes
795
796 CALL para_env%sum(energy%dft_plus_u)
797
798 IF (energy%dft_plus_u < 0.0_dp) &
799 CALL cp_warn(__location__, &
800 "DFT+U energy contribution is negative possibly due "// &
801 "to unphysical Lowdin charges!")
802
803 ! Release (local) full matrices
804 NULLIFY (fm_s_half)
805 CALL cp_fm_release(matrix=fm_work1)
806 CALL cp_fm_release(matrix=fm_work2)
807
808 ! Release (local) sparse matrices
809 CALL dbcsr_release(sm_q)
810 CALL dbcsr_release(sm_v)
811
812 CALL timestop(handle)
813
814 END SUBROUTINE lowdin
815
816! **************************************************************************************************
817!> \brief Add a DFT+U contribution to the Hamiltonian matrix\n
818!> using a method based on the Mulliken population analysis
819!> \f[q_{\mu\nu} = \frac{1}{2} (P_{\mu\nu} S_{\nu\mu} +
820!> S_{\mu\nu} P_{\nu\mu})\f]
821!> where \b P and \b S are the density and the
822!> overlap matrix, respectively.
823!> \param[in] qs_env Quickstep environment
824!> \param orthonormal_basis ...
825!> \param[in,out] matrix_h Hamiltonian matrices for each spin
826!> \param should_output ...
827!> \param output_unit ...
828!> \param print_level ...
829!> \date 03.07.2008
830!> \par
831!> \f{eqnarray*}{
832!> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}\\\
833!> & = & E^{\rm DFT} + \frac{1}{2}\sum_A(U_A - J_A)(Tr(q_A) - Tr(q^2_A))\\[1ex]
834!> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
835!> & = & \frac{\partial E^{\rm DFT}}
836!> {\partial P_{\mu\nu}} +
837!> \frac{\partial E^{\rm U}}
838!> {\partial P_{\mu\nu}}\\\
839!> & = & H_{\mu\nu} + \sum_A
840!> \frac{\partial E^{\rm U}}{\partial q_A}
841!> \frac{\partial q_A}{\partial P_{\mu\nu}}\\\
842!> \f}
843!> \author Matthias Krack (MK)
844!> \version 1.0
845! **************************************************************************************************
846 SUBROUTINE mulliken(qs_env, orthonormal_basis, matrix_h, should_output, &
847 output_unit, print_level)
848
849 TYPE(qs_environment_type), POINTER :: qs_env
850 LOGICAL, INTENT(IN) :: orthonormal_basis
851 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
852 POINTER :: matrix_h
853 LOGICAL, INTENT(IN) :: should_output
854 INTEGER, INTENT(IN) :: output_unit, print_level
855
856 CHARACTER(LEN=*), PARAMETER :: routinen = 'mulliken'
857
858 CHARACTER(LEN=10) :: spin_info
859 CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: symbol
860 CHARACTER(LEN=default_string_length) :: atomic_kind_name
861 INTEGER :: atom_a, handle, i, i0, iatom, ic, ikind, iorb, isb, iset, isgf, ishell, ispin, j, &
862 jsb, jset, jsgf, jshell, lu, m, max_scf, n, natom, natom_of_kind, nimg, nkind, norb, nsb, &
863 nsbsize, nset, nsgf_kind, nspin
864 INTEGER, DIMENSION(1) :: iloc
865 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals
866 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
867 LOGICAL :: debug, dft_plus_u_atom, found, &
868 just_energy, occupation_enforced, smear
869 LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_plus_u_kind, orb_occ
870 REAL(kind=dp) :: eps_scf, eps_u_ramping, fspin, occ, trq, &
871 trq2, u_minus_j, u_minus_j_target, &
872 u_ramping
873 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: q_eigval
874 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: q_eigvec, q_matrix, q_work
875 REAL(kind=dp), DIMENSION(:), POINTER :: nelec
876 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, p_block, q_block, s_block, &
877 v_block
878 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
879 TYPE(atomic_kind_type), POINTER :: kind_a
880 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
881 TYPE(dbcsr_type), POINTER :: sm_h, sm_p, sm_q, sm_s, sm_v
882 TYPE(dft_control_type), POINTER :: dft_control
883 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
884 TYPE(mp_para_env_type), POINTER :: para_env
885 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
886 TYPE(qs_energy_type), POINTER :: energy
887 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
888 TYPE(qs_rho_type), POINTER :: rho
889
890 CALL timeset(routinen, handle)
891
892 debug = .false. ! Set to .TRUE. to print debug information
893
894 NULLIFY (atom_list)
895 NULLIFY (atomic_kind_set)
896 NULLIFY (qs_kind_set)
897 NULLIFY (dft_control)
898 NULLIFY (energy)
899 NULLIFY (first_sgf)
900 NULLIFY (h_block)
901 NULLIFY (matrix_p)
902 NULLIFY (matrix_s)
903 NULLIFY (l)
904 NULLIFY (last_sgf)
905 NULLIFY (nelec)
906 NULLIFY (nshell)
907 NULLIFY (orb_basis_set)
908 NULLIFY (p_block)
909 NULLIFY (particle_set)
910 NULLIFY (q_block)
911 NULLIFY (rho)
912 NULLIFY (s_block)
913 NULLIFY (orbitals)
914 NULLIFY (sm_h)
915 NULLIFY (sm_p)
916 NULLIFY (sm_q)
917 NULLIFY (sm_s)
918 NULLIFY (sm_v)
919 NULLIFY (v_block)
920 NULLIFY (para_env)
921
922 smear = .false.
923 max_scf = -1
924 eps_scf = 1.0e30_dp
925 occupation_enforced = .false.
926
927 CALL get_qs_env(qs_env=qs_env, &
928 atomic_kind_set=atomic_kind_set, &
929 qs_kind_set=qs_kind_set, &
930 dft_control=dft_control, &
931 energy=energy, &
932 particle_set=particle_set, &
933 rho=rho, &
934 para_env=para_env)
935
936 cpassert(ASSOCIATED(atomic_kind_set))
937 cpassert(ASSOCIATED(dft_control))
938 cpassert(ASSOCIATED(energy))
939 cpassert(ASSOCIATED(particle_set))
940 cpassert(ASSOCIATED(rho))
941
942 IF (orthonormal_basis) THEN
943 NULLIFY (sm_s)
944 ELSE
945 ! Get overlap matrix in sparse format
946 CALL get_qs_env(qs_env=qs_env, &
947 matrix_s_kp=matrix_s)
948 cpassert(ASSOCIATED(matrix_s))
949 END IF
950 nimg = dft_control%nimages
951
952 ! Get density matrices in sparse format
953
954 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
955
956 energy%dft_plus_u = 0.0_dp
957
958 nspin = dft_control%nspins
959
960 IF (nspin == 2) THEN
961 fspin = 1.0_dp
962 ELSE
963 fspin = 0.5_dp
964 END IF
965
966 ! Get the total number of atoms, contracted spherical Gaussian basis
967 ! functions, and atomic kinds
968
969 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
970 natom=natom)
971
972 nkind = SIZE(atomic_kind_set)
973
974 ALLOCATE (is_plus_u_kind(nkind))
975 is_plus_u_kind(:) = .false.
976
977 IF (PRESENT(matrix_h)) THEN
978 just_energy = .false.
979 ELSE
980 just_energy = .true.
981 END IF
982
983 ! Loop over all spins
984 DO ispin = 1, nspin
985
986 ! Loop over cell images
987 DO ic = 1, nimg
988 IF (.NOT. orthonormal_basis) THEN
989 sm_s => matrix_s(1, ic)%matrix
990 END IF
991
992 IF (PRESENT(matrix_h)) THEN
993 ! Hamiltonian matrix for spin ispin in sparse format
994 sm_h => matrix_h(ispin, ic)%matrix
995 ELSE
996 NULLIFY (sm_h)
997 END IF
998
999 ! Get density matrix for spin ispin in sparse format
1000
1001 sm_p => matrix_p(ispin, ic)%matrix
1002
1003 IF (.NOT. ASSOCIATED(sm_q)) THEN
1004 ALLOCATE (sm_q)
1005 CALL dbcsr_get_block_diag(sm_p, sm_q)
1006 END IF
1007 CALL dbcsr_set(sm_q, 0.0_dp)
1008
1009 IF (.NOT. ASSOCIATED(sm_v)) THEN
1010 ALLOCATE (sm_v)
1011 CALL dbcsr_get_block_diag(sm_p, sm_v)
1012 END IF
1013 CALL dbcsr_set(sm_v, 0.0_dp)
1014
1015 DO iatom = 1, natom
1016
1017 CALL dbcsr_get_block_p(matrix=sm_p, &
1018 row=iatom, &
1019 col=iatom, &
1020 block=p_block, &
1021 found=found)
1022
1023 IF (.NOT. ASSOCIATED(p_block)) cycle
1024
1025 CALL dbcsr_get_block_p(matrix=sm_q, &
1026 row=iatom, &
1027 col=iatom, &
1028 block=q_block, &
1029 found=found)
1030 cpassert(ASSOCIATED(q_block))
1031
1032 IF (orthonormal_basis) THEN
1033 ! S is the unit matrix
1034 DO isgf = 1, SIZE(q_block, 1)
1035 q_block(isgf, isgf) = p_block(isgf, isgf)
1036 END DO
1037 ELSE
1038 CALL dbcsr_get_block_p(matrix=sm_s, &
1039 row=iatom, &
1040 col=iatom, &
1041 block=s_block, &
1042 found=found)
1043 cpassert(ASSOCIATED(s_block))
1044 ! Exploit that P and S are symmetric
1045 DO jsgf = 1, SIZE(p_block, 2)
1046 DO isgf = 1, SIZE(p_block, 1)
1047 q_block(isgf, jsgf) = p_block(isgf, jsgf)*s_block(isgf, jsgf)
1048 END DO
1049 END DO
1050 END IF ! orthonormal basis set
1051
1052 END DO ! next atom "iatom"
1053
1054 ! E[DFT+U] = E[DFT] + E[U]
1055 ! = E[DFT] + (U - J)*(Tr(q) - Tr(q*q))/2
1056
1057 ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U]
1058 ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j)
1059 ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j))
1060
1061 ! Loop over all atomic kinds
1062
1063 DO ikind = 1, nkind
1064
1065 ! Load the required atomic kind data
1066
1067 CALL get_atomic_kind(atomic_kind_set(ikind), &
1068 atom_list=atom_list, &
1069 name=atomic_kind_name, &
1070 natom=natom_of_kind)
1071
1072 CALL get_qs_kind(qs_kind_set(ikind), &
1073 dft_plus_u_atom=dft_plus_u_atom, &
1074 l_of_dft_plus_u=lu, &
1075 nsgf=nsgf_kind, &
1076 basis_set=orb_basis_set, &
1077 u_minus_j=u_minus_j, &
1078 u_minus_j_target=u_minus_j_target, &
1079 u_ramping=u_ramping, &
1080 eps_u_ramping=eps_u_ramping, &
1081 nelec=nelec, &
1082 orbitals=orbitals, &
1083 eps_scf=eps_scf, &
1084 max_scf=max_scf, &
1085 smear=smear)
1086
1087 ! Check, if the atoms of this atomic kind need a DFT+U correction
1088
1089 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
1090 IF (.NOT. dft_plus_u_atom) cycle
1091 IF (lu < 0) cycle
1092
1093 ! Apply U ramping if requested
1094
1095 IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
1096 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
1097 u_minus_j = min(u_minus_j + u_ramping, u_minus_j_target)
1098 CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
1099 END IF
1100 IF (should_output .AND. (output_unit > 0)) THEN
1101 WRITE (unit=output_unit, fmt="(T3,A,3X,A,F0.3,A)") &
1102 "Kind name: "//trim(adjustl(atomic_kind_name)), &
1103 "U(eff) = ", u_minus_j*evolt, " eV"
1104 END IF
1105 END IF
1106
1107 IF (u_minus_j == 0.0_dp) cycle
1108
1109 is_plus_u_kind(ikind) = .true.
1110
1111 ! Load the required Gaussian basis set data
1112
1113 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1114 first_sgf=first_sgf, &
1115 l=l, &
1116 last_sgf=last_sgf, &
1117 nset=nset, &
1118 nshell=nshell)
1119
1120 ! Count the relevant shell blocks of this atomic kind
1121
1122 nsb = 0
1123 DO iset = 1, nset
1124 DO ishell = 1, nshell(iset)
1125 IF (l(ishell, iset) == lu) nsb = nsb + 1
1126 END DO
1127 END DO
1128
1129 nsbsize = (2*lu + 1)
1130 n = nsb*nsbsize
1131
1132 ALLOCATE (q_matrix(n, n))
1133 q_matrix(:, :) = 0.0_dp
1134
1135 ! Print headline if requested
1136
1137 IF (should_output .AND. (print_level > low_print_level)) THEN
1138 IF (output_unit > 0) THEN
1139 ALLOCATE (symbol(nsbsize))
1140 DO m = -lu, lu
1141 symbol(lu + m + 1) = sgf_symbol(0, lu, m)
1142 END DO
1143 IF (nspin > 1) THEN
1144 WRITE (unit=spin_info, fmt="(A8,I2)") " of spin", ispin
1145 ELSE
1146 spin_info = ""
1147 END IF
1148 WRITE (unit=output_unit, fmt="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
1149 "DFT+U occupations"//trim(spin_info)//" for the atoms of atomic kind ", ikind, &
1150 ": "//trim(atomic_kind_name), &
1151 "Atom Shell ", (adjustr(symbol(i)), i=1, nsbsize), " Trace"
1152 DEALLOCATE (symbol)
1153 END IF
1154 END IF
1155
1156 ! Loop over all atoms of the current atomic kind
1157
1158 DO iatom = 1, natom_of_kind
1159
1160 atom_a = atom_list(iatom)
1161
1162 q_matrix(:, :) = 0.0_dp
1163
1164 ! Get diagonal block
1165
1166 CALL dbcsr_get_block_p(matrix=sm_q, &
1167 row=atom_a, &
1168 col=atom_a, &
1169 block=q_block, &
1170 found=found)
1171
1172 ! Calculate energy contribution to E(U)
1173
1174 IF (ASSOCIATED(q_block)) THEN
1175
1176 i = 0
1177 DO iset = 1, nset
1178 DO ishell = 1, nshell(iset)
1179 IF (l(ishell, iset) /= lu) cycle
1180 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1181 i = i + 1
1182 j = 0
1183 DO jset = 1, nset
1184 DO jshell = 1, nshell(jset)
1185 IF (l(jshell, jset) /= lu) cycle
1186 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
1187 j = j + 1
1188 q_matrix(i, j) = q_block(isgf, jsgf)
1189 END DO ! next contracted spherical Gaussian function "jsgf"
1190 END DO ! next shell "jshell"
1191 END DO ! next shell set "jset"
1192 END DO ! next contracted spherical Gaussian function "isgf"
1193 END DO ! next shell "ishell"
1194 END DO ! next shell set "iset"
1195
1196 ! Perform the requested manipulations of the (initial) orbital occupations
1197
1198 IF (ASSOCIATED(orbitals)) THEN
1199 IF ((qs_env%scf_env%iter_delta >= eps_scf) .OR. &
1200 ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
1201 (qs_env%scf_env%iter_count <= max_scf))) THEN
1202 ALLOCATE (orb_occ(nsbsize))
1203 ALLOCATE (q_eigval(n))
1204 q_eigval(:) = 0.0_dp
1205 ALLOCATE (q_eigvec(n, n))
1206 q_eigvec(:, :) = 0.0_dp
1207 norb = SIZE(orbitals)
1208 CALL jacobi(q_matrix, q_eigval, q_eigvec)
1209 q_matrix(:, :) = 0.0_dp
1210 IF (nelec(ispin) >= 0.5_dp) THEN
1211 trq = nelec(ispin)/sum(q_eigval(1:n))
1212 q_eigval(1:n) = trq*q_eigval(1:n)
1213 END IF
1214 DO isb = 1, nsb
1215 trq = 0.0_dp
1216 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
1217 trq = trq + q_eigval(i)
1218 END DO
1219 IF (smear) THEN
1220 occ = trq/real(norb, kind=dp)
1221 ELSE
1222 occ = 1.0_dp/fspin
1223 END IF
1224 orb_occ(:) = .false.
1225 iloc = maxloc(q_eigvec(:, isb*nsbsize))
1226 jsb = int((iloc(1) - 1)/nsbsize) + 1
1227 i = 0
1228 i0 = (jsb - 1)*nsbsize + 1
1229 iorb = -1000
1230 DO j = i0, jsb*nsbsize
1231 i = i + 1
1232 IF (i > norb) THEN
1233 DO m = -lu, lu
1234 IF (.NOT. orb_occ(lu + m + 1)) THEN
1235 iorb = i0 + lu + m
1236 orb_occ(lu + m + 1) = .true.
1237 END IF
1238 END DO
1239 ELSE
1240 iorb = i0 + lu + orbitals(i)
1241 orb_occ(lu + orbitals(i) + 1) = .true.
1242 END IF
1243 cpassert(iorb /= -1000)
1244 iloc = maxloc(q_eigvec(iorb, :))
1245 q_eigval(iloc(1)) = min(occ, trq)
1246 q_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1)) ! backtransform left
1247 trq = trq - q_eigval(iloc(1))
1248 END DO
1249 END DO
1250 q_matrix(:, :) = matmul(q_matrix, transpose(q_eigvec)) ! backtransform right
1251 DEALLOCATE (orb_occ)
1252 DEALLOCATE (q_eigval)
1253 DEALLOCATE (q_eigvec)
1254 occupation_enforced = .true.
1255 END IF
1256 END IF ! orbitals associated
1257
1258 trq = 0.0_dp
1259 trq2 = 0.0_dp
1260
1261 DO i = 1, n
1262 trq = trq + q_matrix(i, i)
1263 DO j = 1, n
1264 trq2 = trq2 + q_matrix(i, j)*q_matrix(j, i)
1265 END DO
1266 END DO
1267
1268 trq = fspin*trq
1269 trq2 = fspin*fspin*trq2
1270
1271 ! Calculate energy contribution to E(U)
1272
1273 energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
1274
1275 ! Calculate potential V(U) = dE(U)/dq
1276
1277 IF (.NOT. just_energy) THEN
1278
1279 CALL dbcsr_get_block_p(matrix=sm_v, &
1280 row=atom_a, &
1281 col=atom_a, &
1282 block=v_block, &
1283 found=found)
1284 cpassert(ASSOCIATED(v_block))
1285
1286 i = 0
1287 DO iset = 1, nset
1288 DO ishell = 1, nshell(iset)
1289 IF (l(ishell, iset) /= lu) cycle
1290 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1291 i = i + 1
1292 j = 0
1293 DO jset = 1, nset
1294 DO jshell = 1, nshell(jset)
1295 IF (l(jshell, jset) /= lu) cycle
1296 DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
1297 j = j + 1
1298 IF (isgf == jsgf) THEN
1299 v_block(isgf, isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i, i))
1300 ELSE
1301 v_block(isgf, jsgf) = -u_minus_j*fspin*q_matrix(j, i)
1302 END IF
1303 END DO ! next contracted spherical Gaussian function "jsgf"
1304 END DO ! next shell "jshell"
1305 END DO ! next shell set "jset"
1306 END DO ! next contracted spherical Gaussian function "isgf"
1307 END DO ! next shell "ishell"
1308 END DO ! next shell set "iset"
1309
1310 END IF ! not just energy
1311
1312 END IF ! q_block associated
1313
1314 ! Consider print requests
1315
1316 IF (should_output .AND. (print_level > low_print_level)) THEN
1317 CALL para_env%sum(q_matrix)
1318 IF (output_unit > 0) THEN
1319 ALLOCATE (q_work(nsb, nsbsize))
1320 q_work(:, :) = 0.0_dp
1321 DO isb = 1, nsb
1322 j = 0
1323 DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
1324 j = j + 1
1325 q_work(isb, j) = q_matrix(i, i)
1326 END DO
1327 END DO
1328 DO isb = 1, nsb
1329 WRITE (unit=output_unit, fmt="(T3,I6,2X,I6,2X,10F8.3)") &
1330 atom_a, isb, q_work(isb, :), sum(q_work(isb, :))
1331 END DO
1332 WRITE (unit=output_unit, fmt="(T12,A,2X,10F8.3)") &
1333 "Total", (sum(q_work(:, i)), i=1, nsbsize), sum(q_work)
1334 WRITE (unit=output_unit, fmt="(A)") ""
1335 DEALLOCATE (q_work)
1336 IF (debug) THEN
1337 ! Print the DFT+U occupation matrix
1338 WRITE (unit=output_unit, fmt="(T9,70I10)") (i, i=1, n)
1339 DO i = 1, n
1340 WRITE (unit=output_unit, fmt="(T3,I6,70F10.6)") i, q_matrix(i, :)
1341 END DO
1342 ! Print the eigenvalues and eigenvectors of the occupation matrix
1343 ALLOCATE (q_eigval(n))
1344 q_eigval(:) = 0.0_dp
1345 ALLOCATE (q_eigvec(n, n))
1346 q_eigvec(:, :) = 0.0_dp
1347 CALL jacobi(q_matrix, q_eigval, q_eigvec)
1348 WRITE (unit=output_unit, fmt="(/,T9,70I10)") (i, i=1, n)
1349 WRITE (unit=output_unit, fmt="(T9,71F10.6)") (q_eigval(i), i=1, n), &
1350 sum(q_eigval(1:n))
1351 DO i = 1, n
1352 WRITE (unit=output_unit, fmt="(T3,I6,70F10.6)") i, q_eigvec(i, :)
1353 END DO
1354 DEALLOCATE (q_eigval)
1355 DEALLOCATE (q_eigvec)
1356 END IF ! debug
1357 END IF
1358 IF (debug) THEN
1359 ! Print the full atomic occupation matrix block
1360 ALLOCATE (q_work(nsgf_kind, nsgf_kind))
1361 q_work(:, :) = 0.0_dp
1362 IF (ASSOCIATED(q_block)) q_work(:, :) = q_block(:, :)
1363 CALL para_env%sum(q_work)
1364 IF (output_unit > 0) THEN
1365 norb = SIZE(q_work, 1)
1366 WRITE (unit=output_unit, fmt="(/,T9,200I10)") (i, i=1, norb)
1367 DO i = 1, norb
1368 WRITE (unit=output_unit, fmt="(T3,I6,200F10.6)") i, q_work(i, :)
1369 END DO
1370 ALLOCATE (q_eigval(norb))
1371 q_eigval(:) = 0.0_dp
1372 ALLOCATE (q_eigvec(norb, norb))
1373 q_eigvec(:, :) = 0.0_dp
1374 CALL jacobi(q_work, q_eigval, q_eigvec)
1375 WRITE (unit=output_unit, fmt="(/,T9,200I10)") (i, i=1, norb)
1376 WRITE (unit=output_unit, fmt="(T9,201F10.6)") (q_eigval(i), i=1, norb), &
1377 sum(q_eigval(1:norb))
1378 DO i = 1, norb
1379 WRITE (unit=output_unit, fmt="(T3,I6,200F10.6)") i, q_eigvec(i, :)
1380 END DO
1381 DEALLOCATE (q_eigval)
1382 DEALLOCATE (q_eigvec)
1383 END IF
1384 DEALLOCATE (q_work)
1385 END IF ! debug
1386 END IF ! should output
1387
1388 END DO ! next atom "iatom" of atomic kind "ikind"
1389
1390 IF (ALLOCATED(q_matrix)) THEN
1391 DEALLOCATE (q_matrix)
1392 END IF
1393
1394 END DO ! next atomic kind "ikind"
1395
1396 ! Add V(i,j)[U] to V(i,j)[DFT]
1397
1398 IF (ASSOCIATED(sm_h)) THEN
1399
1400 DO ikind = 1, nkind
1401
1402 IF (.NOT. is_plus_u_kind(ikind)) cycle
1403
1404 kind_a => atomic_kind_set(ikind)
1405
1406 CALL get_atomic_kind(atomic_kind=kind_a, &
1407 atom_list=atom_list, &
1408 natom=natom_of_kind)
1409
1410 DO iatom = 1, natom_of_kind
1411
1412 atom_a = atom_list(iatom)
1413
1414 CALL dbcsr_get_block_p(matrix=sm_h, &
1415 row=atom_a, &
1416 col=atom_a, &
1417 block=h_block, &
1418 found=found)
1419
1420 IF (.NOT. ASSOCIATED(h_block)) cycle
1421
1422 CALL dbcsr_get_block_p(matrix=sm_v, &
1423 row=atom_a, &
1424 col=atom_a, &
1425 block=v_block, &
1426 found=found)
1427 cpassert(ASSOCIATED(v_block))
1428
1429 IF (orthonormal_basis) THEN
1430 DO isgf = 1, SIZE(h_block, 1)
1431 h_block(isgf, isgf) = h_block(isgf, isgf) + v_block(isgf, isgf)
1432 END DO
1433 ELSE
1434 CALL dbcsr_get_block_p(matrix=sm_s, &
1435 row=atom_a, &
1436 col=atom_a, &
1437 block=s_block, &
1438 found=found)
1439 cpassert(ASSOCIATED(s_block))
1440 DO jsgf = 1, SIZE(h_block, 2)
1441 DO isgf = 1, SIZE(h_block, 1)
1442 h_block(isgf, jsgf) = h_block(isgf, jsgf) + v_block(isgf, jsgf)*s_block(isgf, jsgf)
1443 END DO
1444 END DO
1445 END IF ! orthonormal basis set
1446
1447 END DO ! next atom "iatom" of atomic kind "ikind"
1448
1449 END DO ! Next atomic kind "ikind"
1450
1451 END IF ! An update of the Hamiltonian matrix is requested
1452
1453 END DO ! next cell image
1454
1455 END DO ! next spin "ispin"
1456
1457 ! Collect the energy contributions from all processes
1458
1459 CALL para_env%sum(energy%dft_plus_u)
1460
1461 IF (energy%dft_plus_u < 0.0_dp) THEN
1462 IF (.NOT. occupation_enforced) THEN
1463 CALL cp_warn(__location__, &
1464 "DFT+U energy contribution is negative possibly due "// &
1465 "to unphysical Mulliken charges!")
1466 END IF
1467 END IF
1468
1469 CALL dbcsr_deallocate_matrix(sm_q)
1470 CALL dbcsr_deallocate_matrix(sm_v)
1471
1472 CALL timestop(handle)
1473
1474 END SUBROUTINE mulliken
1475
1476! **************************************************************************************************
1477!> \brief Add a DFT+U contribution to the Hamiltonian matrix\n
1478!> using a method based on Mulliken charges
1479!> \f[q_\mu = \sum_\nu \frac{1}{2}(P_{\mu\nu} S_{\nu\mu} +
1480!> S_{\mu\nu} P_{\nu\mu})
1481!> = \sum_\nu P_{\mu\nu} S_{\nu\mu}\f]
1482!> where \b P and \b S are the density and the
1483!> overlap matrix, respectively.
1484!> \param[in] qs_env Quickstep environment
1485!> \param orthonormal_basis ...
1486!> \param[in,out] matrix_h Hamiltonian matrices for each spin
1487!> \param[in,out] matrix_w Energy weighted density matrices for each spin
1488!> \param should_output ...
1489!> \param output_unit ...
1490!> \param print_level ...
1491!> \date 11.01.2008
1492!> \par
1493!> \f{eqnarray*}{
1494!> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}\\\
1495!> & = & E^{\rm DFT} + \frac{1}{2}(U - J)\sum_\mu (q_\mu - q_\mu^2)\\[1ex]
1496!> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
1497!> & = & \frac{\partial E^{\rm DFT}}
1498!> {\partial P_{\mu\nu}} +
1499!> \frac{\partial E^{\rm U}}
1500!> {\partial P_{\mu\nu}}\\\
1501!> & = & H_{\mu\nu} +
1502!> \frac{\partial E^{\rm U}}{\partial q_\mu}
1503!> \frac{\partial q_\mu}{\partial P_{\mu\nu}}\\\
1504!> & = & H_{\mu\nu} +
1505!> \frac{1}{2}(U - J)(1 - q_\mu - q_\nu) S_{\mu\nu}\\\
1506!> \f}
1507!> \author Matthias Krack (MK)
1508!> \version 1.0
1509!> \note The use of any full matrices was avoided. Thus no ScaLAPACK
1510!> calls are performed
1511! **************************************************************************************************
1512 SUBROUTINE mulliken_charges(qs_env, orthonormal_basis, matrix_h, matrix_w, &
1513 should_output, output_unit, print_level)
1514
1515 TYPE(qs_environment_type), POINTER :: qs_env
1516 LOGICAL, INTENT(IN) :: orthonormal_basis
1517 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
1518 POINTER :: matrix_h, matrix_w
1519 LOGICAL, INTENT(IN) :: should_output
1520 INTEGER, INTENT(IN) :: output_unit, print_level
1521
1522 CHARACTER(LEN=*), PARAMETER :: routinen = 'mulliken_charges'
1523
1524 CHARACTER(LEN=10) :: spin_info
1525 CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: symbol
1526 CHARACTER(LEN=default_string_length) :: atomic_kind_name
1527 INTEGER :: atom_a, handle, i, iatom, ic, ikind, isb, iset, isgf, ishell, ispin, jatom, jsgf, &
1528 lu, m, natom, natom_of_kind, nimg, nkind, nsb, nset, nsgf, nspin, sgf
1529 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
1530 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell
1531 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
1532 LOGICAL :: dft_plus_u_atom, found, just_energy
1533 REAL(kind=dp) :: eps_u_ramping, fspin, q, u_minus_j, &
1534 u_minus_j_target, u_ramping, v
1535 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dedq, trps
1536 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: q_ii
1537 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, p_block, s_block, w_block
1538 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1539 TYPE(dbcsr_iterator_type) :: iter
1540 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
1541 TYPE(dbcsr_type), POINTER :: sm_h, sm_p, sm_s, sm_w
1542 TYPE(dft_control_type), POINTER :: dft_control
1543 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1544 TYPE(mp_para_env_type), POINTER :: para_env
1545 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1546 TYPE(qs_energy_type), POINTER :: energy
1547 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1548 TYPE(qs_rho_type), POINTER :: rho
1549
1550 CALL timeset(routinen, handle)
1551
1552 NULLIFY (atom_list)
1553 NULLIFY (atomic_kind_set)
1554 NULLIFY (qs_kind_set)
1555 NULLIFY (dft_control)
1556 NULLIFY (energy)
1557 NULLIFY (first_sgf)
1558 NULLIFY (h_block)
1559 NULLIFY (matrix_p)
1560 NULLIFY (matrix_s)
1561 NULLIFY (l)
1562 NULLIFY (last_sgf)
1563 NULLIFY (nshell)
1564 NULLIFY (orb_basis_set)
1565 NULLIFY (p_block)
1566 NULLIFY (particle_set)
1567 NULLIFY (rho)
1568 NULLIFY (s_block)
1569 NULLIFY (sm_h)
1570 NULLIFY (sm_p)
1571 NULLIFY (sm_s)
1572 NULLIFY (w_block)
1573 NULLIFY (para_env)
1574
1575 CALL get_qs_env(qs_env=qs_env, &
1576 atomic_kind_set=atomic_kind_set, &
1577 qs_kind_set=qs_kind_set, &
1578 dft_control=dft_control, &
1579 energy=energy, &
1580 particle_set=particle_set, &
1581 rho=rho, &
1582 para_env=para_env)
1583
1584 cpassert(ASSOCIATED(atomic_kind_set))
1585 cpassert(ASSOCIATED(dft_control))
1586 cpassert(ASSOCIATED(energy))
1587 cpassert(ASSOCIATED(particle_set))
1588 cpassert(ASSOCIATED(rho))
1589
1590 IF (orthonormal_basis) THEN
1591 NULLIFY (sm_s)
1592 ELSE
1593 ! Get overlap matrix in sparse format
1594 CALL get_qs_env(qs_env=qs_env, &
1595 matrix_s_kp=matrix_s)
1596 cpassert(ASSOCIATED(matrix_s))
1597 END IF
1598
1599 ! Get density matrices in sparse format
1600
1601 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1602
1603 energy%dft_plus_u = 0.0_dp
1604
1605 nspin = dft_control%nspins
1606 nimg = dft_control%nimages
1607
1608 IF (nspin == 2) THEN
1609 fspin = 1.0_dp
1610 ELSE
1611 fspin = 0.5_dp
1612 END IF
1613
1614 ! Get the total number of atoms, contracted spherical Gaussian basis
1615 ! functions, and atomic kinds
1616
1617 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
1618 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
1619
1620 nkind = SIZE(atomic_kind_set)
1621
1622 ALLOCATE (first_sgf_atom(natom))
1623 first_sgf_atom(:) = 0
1624
1625 CALL get_particle_set(particle_set, qs_kind_set, &
1626 first_sgf=first_sgf_atom)
1627
1628 ALLOCATE (trps(nsgf))
1629 trps(:) = 0.0_dp
1630
1631 IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
1632 ALLOCATE (dedq(nsgf))
1633 just_energy = .false.
1634 ELSE
1635 just_energy = .true.
1636 END IF
1637
1638 ! Loop over all spins
1639
1640 DO ispin = 1, nspin
1641
1642 IF (.NOT. just_energy) dedq(:) = 0.0_dp
1643
1644 ! Calculate Trace(P*S) assuming symmetric matrices
1645
1646 trps(:) = 0.0_dp
1647
1648 DO ic = 1, nimg
1649 IF (orthonormal_basis) THEN
1650 NULLIFY (sm_s)
1651 ELSE
1652 sm_s => matrix_s(1, ic)%matrix
1653 END IF
1654 sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
1655
1656 CALL dbcsr_iterator_start(iter, sm_p)
1657
1658 DO WHILE (dbcsr_iterator_blocks_left(iter))
1659
1660 CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
1661
1662 IF (orthonormal_basis) THEN
1663
1664 IF (iatom /= jatom) cycle
1665
1666 IF (ASSOCIATED(p_block)) THEN
1667 sgf = first_sgf_atom(iatom)
1668 DO isgf = 1, SIZE(p_block, 1)
1669 trps(sgf) = trps(sgf) + p_block(isgf, isgf)
1670 sgf = sgf + 1
1671 END DO
1672 END IF
1673
1674 ELSE
1675
1676 CALL dbcsr_get_block_p(matrix=sm_s, &
1677 row=iatom, &
1678 col=jatom, &
1679 block=s_block, &
1680 found=found)
1681 cpassert(ASSOCIATED(s_block))
1682
1683 sgf = first_sgf_atom(jatom)
1684 DO jsgf = 1, SIZE(p_block, 2)
1685 DO isgf = 1, SIZE(p_block, 1)
1686 trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
1687 END DO
1688 sgf = sgf + 1
1689 END DO
1690
1691 IF (iatom /= jatom) THEN
1692 sgf = first_sgf_atom(iatom)
1693 DO isgf = 1, SIZE(p_block, 1)
1694 DO jsgf = 1, SIZE(p_block, 2)
1695 trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
1696 END DO
1697 sgf = sgf + 1
1698 END DO
1699 END IF
1700
1701 END IF ! orthonormal basis set
1702
1703 END DO ! next atom "iatom"
1704
1705 CALL dbcsr_iterator_stop(iter)
1706
1707 END DO ! cell images
1708
1709 CALL para_env%sum(trps)
1710
1711 ! q <- Trace(PS)
1712
1713 ! E[DFT+U] = E[DFT] + E[U]
1714 ! = E[DFT] + (U - J)*(q - q**2))/2
1715
1716 ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U]
1717 ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j)
1718 ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j))
1719
1720 ! Loop over all atomic kinds
1721
1722 DO ikind = 1, nkind
1723
1724 ! Load the required atomic kind data
1725 CALL get_atomic_kind(atomic_kind_set(ikind), &
1726 atom_list=atom_list, &
1727 name=atomic_kind_name, &
1728 natom=natom_of_kind)
1729
1730 CALL get_qs_kind(qs_kind_set(ikind), &
1731 dft_plus_u_atom=dft_plus_u_atom, &
1732 l_of_dft_plus_u=lu, &
1733 basis_set=orb_basis_set, &
1734 u_minus_j=u_minus_j, &
1735 u_minus_j_target=u_minus_j_target, &
1736 u_ramping=u_ramping, &
1737 eps_u_ramping=eps_u_ramping)
1738
1739 ! Check, if this atom needs a DFT+U correction
1740
1741 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
1742 IF (.NOT. dft_plus_u_atom) cycle
1743 IF (lu < 0) cycle
1744
1745 ! Apply U ramping if requested
1746
1747 IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
1748 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
1749 u_minus_j = min(u_minus_j + u_ramping, u_minus_j_target)
1750 CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
1751 END IF
1752 IF (should_output .AND. (output_unit > 0)) THEN
1753 WRITE (unit=output_unit, fmt="(T3,A,3X,A,F0.3,A)") &
1754 "Kind name: "//trim(adjustl(atomic_kind_name)), &
1755 "U(eff) = ", u_minus_j*evolt, " eV"
1756 END IF
1757 END IF
1758
1759 IF (u_minus_j == 0.0_dp) cycle
1760
1761 ! Load the required Gaussian basis set data
1762
1763 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1764 first_sgf=first_sgf, &
1765 l=l, &
1766 last_sgf=last_sgf, &
1767 nset=nset, &
1768 nshell=nshell)
1769
1770 ! Count the relevant shell blocks of this atomic kind
1771
1772 nsb = 0
1773 DO iset = 1, nset
1774 DO ishell = 1, nshell(iset)
1775 IF (l(ishell, iset) == lu) nsb = nsb + 1
1776 END DO
1777 END DO
1778
1779 ALLOCATE (q_ii(nsb, 2*lu + 1))
1780
1781 ! Print headline if requested
1782
1783 IF (should_output .AND. (print_level > low_print_level)) THEN
1784 IF (output_unit > 0) THEN
1785 ALLOCATE (symbol(2*lu + 1))
1786 DO m = -lu, lu
1787 symbol(lu + m + 1) = sgf_symbol(0, lu, m)
1788 END DO
1789 IF (nspin > 1) THEN
1790 WRITE (unit=spin_info, fmt="(A8,I2)") " of spin", ispin
1791 ELSE
1792 spin_info = ""
1793 END IF
1794 WRITE (unit=output_unit, fmt="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
1795 "DFT+U occupations"//trim(spin_info)//" for the atoms of atomic kind ", ikind, &
1796 ": "//trim(atomic_kind_name), &
1797 "Atom Shell ", (adjustr(symbol(i)), i=1, 2*lu + 1), " Trace"
1798 DEALLOCATE (symbol)
1799 END IF
1800 END IF
1801
1802 ! Loop over all atoms of the current atomic kind
1803
1804 DO iatom = 1, natom_of_kind
1805
1806 atom_a = atom_list(iatom)
1807
1808 q_ii(:, :) = 0.0_dp
1809
1810 ! Get diagonal block
1811
1812 CALL dbcsr_get_block_p(matrix=sm_p, &
1813 row=atom_a, &
1814 col=atom_a, &
1815 block=p_block, &
1816 found=found)
1817
1818 ! Calculate E(U) and dE(U)/dq
1819
1820 IF (ASSOCIATED(p_block)) THEN
1821
1822 sgf = first_sgf_atom(atom_a)
1823
1824 isb = 0
1825 DO iset = 1, nset
1826 DO ishell = 1, nshell(iset)
1827 IF (l(ishell, iset) == lu) THEN
1828 isb = isb + 1
1829 i = 0
1830 DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1831 q = fspin*trps(sgf)
1832 i = i + 1
1833 q_ii(isb, i) = q
1834 energy%dft_plus_u = energy%dft_plus_u + &
1835 0.5_dp*u_minus_j*(q - q**2)/fspin
1836 IF (.NOT. just_energy) THEN
1837 dedq(sgf) = dedq(sgf) + u_minus_j*(0.5_dp - q)
1838 END IF
1839 sgf = sgf + 1
1840 END DO ! next contracted spherical Gaussian function "isgf"
1841 ELSE
1842 sgf = sgf + last_sgf(ishell, iset) - first_sgf(ishell, iset) + 1
1843 END IF ! angular momentum requested for DFT+U correction
1844 END DO ! next shell "ishell"
1845 END DO ! next shell set "iset"
1846
1847 END IF ! this process is the owner of the sparse matrix block?
1848
1849 ! Consider print requests
1850
1851 IF (should_output .AND. (print_level > low_print_level)) THEN
1852 CALL para_env%sum(q_ii)
1853 IF (output_unit > 0) THEN
1854 DO isb = 1, nsb
1855 WRITE (unit=output_unit, fmt="(T3,I6,2X,I6,2X,10F8.3)") &
1856 atom_a, isb, q_ii(isb, :), sum(q_ii(isb, :))
1857 END DO
1858 WRITE (unit=output_unit, fmt="(T12,A,2X,10F8.3)") &
1859 "Total", (sum(q_ii(:, i)), i=1, 2*lu + 1), sum(q_ii)
1860 WRITE (unit=output_unit, fmt="(A)") ""
1861 END IF
1862 END IF ! should output
1863
1864 END DO ! next atom "iatom" of atomic kind "ikind"
1865
1866 IF (ALLOCATED(q_ii)) THEN
1867 DEALLOCATE (q_ii)
1868 END IF
1869
1870 END DO ! next atomic kind "ikind"
1871
1872 IF (.NOT. just_energy) THEN
1873 CALL para_env%sum(dedq)
1874 END IF
1875
1876 ! Add V(i,j)[U] to V(i,j)[DFT]
1877
1878 IF (PRESENT(matrix_h)) THEN
1879
1880 DO ic = 1, nimg
1881 IF (orthonormal_basis) THEN
1882 NULLIFY (sm_s)
1883 ELSE
1884 sm_s => matrix_s(1, ic)%matrix
1885 END IF
1886 sm_h => matrix_h(ispin, ic)%matrix
1887
1888 CALL dbcsr_iterator_start(iter, sm_h)
1889
1890 DO WHILE (dbcsr_iterator_blocks_left(iter))
1891
1892 CALL dbcsr_iterator_next_block(iter, iatom, jatom, h_block)
1893
1894 IF (orthonormal_basis) THEN
1895
1896 IF (iatom /= jatom) cycle
1897
1898 IF (ASSOCIATED(h_block)) THEN
1899 sgf = first_sgf_atom(iatom)
1900 DO isgf = 1, SIZE(h_block, 1)
1901 h_block(isgf, isgf) = h_block(isgf, isgf) + dedq(sgf)
1902 sgf = sgf + 1
1903 END DO
1904 END IF
1905
1906 ELSE
1907
1908 ! Request katom just to check for consistent sparse matrix pattern
1909
1910 CALL dbcsr_get_block_p(matrix=sm_s, &
1911 row=iatom, &
1912 col=jatom, &
1913 block=s_block, &
1914 found=found)
1915 cpassert(ASSOCIATED(s_block))
1916
1917 ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
1918
1919 sgf = first_sgf_atom(iatom)
1920
1921 DO isgf = 1, SIZE(h_block, 1)
1922 IF (dedq(sgf) /= 0.0_dp) THEN
1923 v = 0.5_dp*dedq(sgf)
1924 DO jsgf = 1, SIZE(h_block, 2)
1925 h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
1926 END DO
1927 END IF
1928 sgf = sgf + 1
1929 END DO
1930
1931 sgf = first_sgf_atom(jatom)
1932
1933 DO jsgf = 1, SIZE(h_block, 2)
1934 IF (dedq(sgf) /= 0.0_dp) THEN
1935 v = 0.5_dp*dedq(sgf)
1936 DO isgf = 1, SIZE(h_block, 1)
1937 h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
1938 END DO
1939 END IF
1940 sgf = sgf + 1
1941 END DO
1942
1943 END IF ! orthonormal basis set
1944
1945 END DO ! Next atom "iatom"
1946
1947 CALL dbcsr_iterator_stop(iter)
1948
1949 END DO
1950
1951 END IF ! An update of the Hamiltonian matrix is requested
1952
1953 ! Calculate the contribution (non-Pulay part) to the derivatives
1954 ! w.r.t. the nuclear positions, which requires an update of the
1955 ! energy weighted density W.
1956
1957 IF (PRESENT(matrix_w) .AND. (.NOT. orthonormal_basis)) THEN
1958
1959 DO ic = 1, nimg
1960 sm_s => matrix_s(1, ic)%matrix
1961 sm_p => matrix_p(ispin, ic)%matrix
1962 sm_w => matrix_w(ispin, ic)%matrix
1963
1964 CALL dbcsr_iterator_start(iter, sm_p)
1965
1966 DO WHILE (dbcsr_iterator_blocks_left(iter))
1967
1968 CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
1969
1970 ! Skip the diagonal blocks of the W matrix
1971
1972 IF (iatom == jatom) cycle
1973
1974 ! Request katom just to check for consistent sparse matrix patterns
1975
1976 CALL dbcsr_get_block_p(matrix=sm_w, &
1977 row=iatom, &
1978 col=jatom, &
1979 block=w_block, &
1980 found=found)
1981 cpassert(ASSOCIATED(w_block))
1982
1983 ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
1984
1985 sgf = first_sgf_atom(iatom)
1986
1987 DO isgf = 1, SIZE(w_block, 1)
1988 IF (dedq(sgf) /= 0.0_dp) THEN
1989 v = -0.5_dp*dedq(sgf)
1990 DO jsgf = 1, SIZE(w_block, 2)
1991 w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
1992 END DO
1993 END IF
1994 sgf = sgf + 1
1995 END DO
1996
1997 sgf = first_sgf_atom(jatom)
1998
1999 DO jsgf = 1, SIZE(w_block, 2)
2000 IF (dedq(sgf) /= 0.0_dp) THEN
2001 v = -0.5_dp*dedq(sgf)
2002 DO isgf = 1, SIZE(w_block, 1)
2003 w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
2004 END DO
2005 END IF
2006 sgf = sgf + 1
2007 END DO
2008
2009 END DO ! next block node "jatom"
2010
2011 CALL dbcsr_iterator_stop(iter)
2012
2013 END DO
2014
2015 END IF ! W matrix update requested
2016
2017 END DO ! next spin "ispin"
2018
2019 ! Collect the energy contributions from all processes
2020
2021 CALL para_env%sum(energy%dft_plus_u)
2022
2023 IF (energy%dft_plus_u < 0.0_dp) &
2024 CALL cp_warn(__location__, &
2025 "DFT+U energy contribution is negative possibly due "// &
2026 "to unphysical Mulliken charges!")
2027
2028 ! Release local work storage
2029
2030 IF (ALLOCATED(first_sgf_atom)) THEN
2031 DEALLOCATE (first_sgf_atom)
2032 END IF
2033
2034 IF (ALLOCATED(trps)) THEN
2035 DEALLOCATE (trps)
2036 END IF
2037
2038 IF (ALLOCATED(dedq)) THEN
2039 DEALLOCATE (dedq)
2040 END IF
2041
2042 CALL timestop(handle)
2043
2044 END SUBROUTINE mulliken_charges
2045
2046END MODULE dft_plus_u
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
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)
...
Definition basis_set_types.F:631
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_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)
...
Definition cp_dbcsr_api.F:969
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:943
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1040
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:692
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1002
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1194
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1132
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:552
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:214
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:695
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:78
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)
...
Definition cp_dbcsr_output.F:163
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:1760
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:781
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:172
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:239
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:1087
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:169
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:789
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:1063
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:103
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::plus_u_lowdin
integer, parameter, public plus_u_lowdin
Definition input_constants.F:898
input_constants::plus_u_mulliken_charges
integer, parameter, public plus_u_mulliken_charges
Definition input_constants.F:898
input_constants::plus_u_mulliken
integer, parameter, public plus_u_mulliken
Definition input_constants.F:898
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(matrix_p, kpoints, ispin, fmwork)
Kpoint transformation of density matrix for Lowdin population analysis Transforms matrices to kpoint,...
Definition kpoint_methods.F:1258
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:1476
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)
Get the components of a particle set.
Definition particle_methods.F:99
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, 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:530
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, 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, 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:468
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, 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:913
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, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
Definition qs_kind_types.F:3026
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:631
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:115
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:156
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
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:220
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
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