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