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