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