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