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