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(:, :) :: 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 600 : 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 9374 : 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 9374 : 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 154 : q_eigval(:) = 0.0_dp
547 66 : ALLOCATE (q_eigvec(n, n))
548 946 : q_eigvec(:, :) = 0.0_dp
549 22 : norb = SIZE(orbitals)
550 22 : CALL jacobi(q_matrix, q_eigval, q_eigvec)
551 946 : 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 176 : 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 31350 : 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) &
770 : CALL cp_warn(__LOCATION__, &
771 : "DFT+U energy contribution is negative possibly due "// &
772 0 : "to unphysical Lowdin charges!")
773 :
774 : ! Release (local) full matrices
775 150 : NULLIFY (fm_s_half)
776 150 : CALL cp_fm_release(matrix=fm_work1)
777 150 : CALL cp_fm_release(matrix=fm_work2)
778 :
779 : ! Release (local) sparse matrices
780 150 : CALL dbcsr_release(sm_q)
781 150 : CALL dbcsr_release(sm_v)
782 :
783 150 : CALL timestop(handle)
784 :
785 600 : END SUBROUTINE lowdin
786 :
787 : ! **************************************************************************************************
788 : !> \brief Add a DFT+U contribution to the Hamiltonian matrix\n
789 : !> using a method based on the Mulliken population analysis
790 : !> \f[q_{\mu\nu} = \frac{1}{2} (P_{\mu\nu} S_{\nu\mu} +
791 : !> S_{\mu\nu} P_{\nu\mu})\f]
792 : !> where \b P and \b S are the density and the
793 : !> overlap matrix, respectively.
794 : !> \param[in] qs_env Quickstep environment
795 : !> \param orthonormal_basis ...
796 : !> \param[in,out] matrix_h Hamiltonian matrices for each spin
797 : !> \param should_output ...
798 : !> \param output_unit ...
799 : !> \param print_level ...
800 : !> \date 03.07.2008
801 : !> \par
802 : !> \f{eqnarray*}{
803 : !> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}\\\
804 : !> & = & E^{\rm DFT} + \frac{1}{2}\sum_A(U_A - J_A)(Tr(q_A) - Tr(q^2_A))\\[1ex]
805 : !> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
806 : !> & = & \frac{\partial E^{\rm DFT}}
807 : !> {\partial P_{\mu\nu}} +
808 : !> \frac{\partial E^{\rm U}}
809 : !> {\partial P_{\mu\nu}}\\\
810 : !> & = & H_{\mu\nu} + \sum_A
811 : !> \frac{\partial E^{\rm U}}{\partial q_A}
812 : !> \frac{\partial q_A}{\partial P_{\mu\nu}}\\\
813 : !> \f}
814 : !> \author Matthias Krack (MK)
815 : !> \version 1.0
816 : ! **************************************************************************************************
817 1290 : SUBROUTINE mulliken(qs_env, orthonormal_basis, matrix_h, should_output, &
818 : output_unit, print_level)
819 :
820 : TYPE(qs_environment_type), POINTER :: qs_env
821 : LOGICAL, INTENT(IN) :: orthonormal_basis
822 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
823 : POINTER :: matrix_h
824 : LOGICAL, INTENT(IN) :: should_output
825 : INTEGER, INTENT(IN) :: output_unit, print_level
826 :
827 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mulliken'
828 :
829 : CHARACTER(LEN=10) :: spin_info
830 1290 : CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: symbol
831 : CHARACTER(LEN=default_string_length) :: atomic_kind_name
832 : INTEGER :: atom_a, handle, i, i0, iatom, ic, ikind, iorb, isb, iset, isgf, ishell, ispin, j, &
833 : jsb, jset, jsgf, jshell, lu, m, max_scf, n, natom, natom_of_kind, nimg, nkind, norb, nsb, &
834 : nsbsize, nset, nsgf_kind, nspin
835 : INTEGER, DIMENSION(1) :: iloc
836 1290 : INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals
837 1290 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
838 : LOGICAL :: debug, dft_plus_u_atom, found, &
839 : just_energy, occupation_enforced, smear
840 1290 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_plus_u_kind, orb_occ
841 : REAL(KIND=dp) :: eps_scf, eps_u_ramping, fspin, occ, trq, &
842 : trq2, u_minus_j, u_minus_j_target, &
843 : u_ramping
844 1290 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: q_eigval
845 1290 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: q_eigvec, q_matrix, q_work
846 1290 : REAL(KIND=dp), DIMENSION(:), POINTER :: nelec
847 1290 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, q_block, s_block, &
848 1290 : v_block
849 1290 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
850 : TYPE(atomic_kind_type), POINTER :: kind_a
851 1290 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
852 : TYPE(dbcsr_type), POINTER :: sm_h, sm_p, sm_q, sm_s, sm_v
853 : TYPE(dft_control_type), POINTER :: dft_control
854 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
855 : TYPE(mp_para_env_type), POINTER :: para_env
856 1290 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
857 : TYPE(qs_energy_type), POINTER :: energy
858 1290 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
859 : TYPE(qs_rho_type), POINTER :: rho
860 :
861 1290 : CALL timeset(routineN, handle)
862 :
863 1290 : debug = .FALSE. ! Set to .TRUE. to print debug information
864 :
865 1290 : NULLIFY (atom_list)
866 1290 : NULLIFY (atomic_kind_set)
867 1290 : NULLIFY (qs_kind_set)
868 1290 : NULLIFY (dft_control)
869 1290 : NULLIFY (energy)
870 1290 : NULLIFY (first_sgf)
871 1290 : NULLIFY (h_block)
872 1290 : NULLIFY (matrix_p)
873 1290 : NULLIFY (matrix_s)
874 1290 : NULLIFY (l)
875 1290 : NULLIFY (last_sgf)
876 1290 : NULLIFY (nelec)
877 1290 : NULLIFY (nshell)
878 1290 : NULLIFY (orb_basis_set)
879 1290 : NULLIFY (p_block)
880 1290 : NULLIFY (particle_set)
881 1290 : NULLIFY (q_block)
882 1290 : NULLIFY (rho)
883 1290 : NULLIFY (s_block)
884 1290 : NULLIFY (orbitals)
885 1290 : NULLIFY (sm_h)
886 1290 : NULLIFY (sm_p)
887 1290 : NULLIFY (sm_q)
888 1290 : NULLIFY (sm_s)
889 1290 : NULLIFY (sm_v)
890 1290 : NULLIFY (v_block)
891 1290 : NULLIFY (para_env)
892 :
893 1290 : smear = .FALSE.
894 1290 : max_scf = -1
895 1290 : eps_scf = 1.0E30_dp
896 1290 : occupation_enforced = .FALSE.
897 :
898 : CALL get_qs_env(qs_env=qs_env, &
899 : atomic_kind_set=atomic_kind_set, &
900 : qs_kind_set=qs_kind_set, &
901 : dft_control=dft_control, &
902 : energy=energy, &
903 : particle_set=particle_set, &
904 : rho=rho, &
905 1290 : para_env=para_env)
906 :
907 1290 : CPASSERT(ASSOCIATED(atomic_kind_set))
908 1290 : CPASSERT(ASSOCIATED(dft_control))
909 1290 : CPASSERT(ASSOCIATED(energy))
910 1290 : CPASSERT(ASSOCIATED(particle_set))
911 1290 : CPASSERT(ASSOCIATED(rho))
912 :
913 1290 : IF (orthonormal_basis) THEN
914 : NULLIFY (sm_s)
915 : ELSE
916 : ! Get overlap matrix in sparse format
917 : CALL get_qs_env(qs_env=qs_env, &
918 1290 : matrix_s_kp=matrix_s)
919 1290 : CPASSERT(ASSOCIATED(matrix_s))
920 : END IF
921 1290 : nimg = dft_control%nimages
922 :
923 : ! Get density matrices in sparse format
924 :
925 1290 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
926 :
927 1290 : energy%dft_plus_u = 0.0_dp
928 :
929 1290 : nspin = dft_control%nspins
930 :
931 1290 : IF (nspin == 2) THEN
932 : fspin = 1.0_dp
933 : ELSE
934 660 : fspin = 0.5_dp
935 : END IF
936 :
937 : ! Get the total number of atoms, contracted spherical Gaussian basis
938 : ! functions, and atomic kinds
939 :
940 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
941 1290 : natom=natom)
942 :
943 1290 : nkind = SIZE(atomic_kind_set)
944 :
945 3870 : ALLOCATE (is_plus_u_kind(nkind))
946 3870 : is_plus_u_kind(:) = .FALSE.
947 :
948 1290 : IF (PRESENT(matrix_h)) THEN
949 : just_energy = .FALSE.
950 : ELSE
951 570 : just_energy = .TRUE.
952 : END IF
953 :
954 : ! Loop over all spins
955 3210 : DO ispin = 1, nspin
956 :
957 : ! Loop over cell images
958 5130 : DO ic = 1, nimg
959 1920 : IF (.NOT. orthonormal_basis) THEN
960 1920 : sm_s => matrix_s(1, ic)%matrix
961 : END IF
962 :
963 1920 : IF (PRESENT(matrix_h)) THEN
964 : ! Hamiltonian matrix for spin ispin in sparse format
965 1072 : sm_h => matrix_h(ispin, ic)%matrix
966 : ELSE
967 : NULLIFY (sm_h)
968 : END IF
969 :
970 : ! Get density matrix for spin ispin in sparse format
971 :
972 1920 : sm_p => matrix_p(ispin, ic)%matrix
973 :
974 1920 : IF (.NOT. ASSOCIATED(sm_q)) THEN
975 1290 : ALLOCATE (sm_q)
976 1290 : CALL dbcsr_get_block_diag(sm_p, sm_q)
977 : END IF
978 1920 : CALL dbcsr_set(sm_q, 0.0_dp)
979 :
980 1920 : IF (.NOT. ASSOCIATED(sm_v)) THEN
981 1290 : ALLOCATE (sm_v)
982 1290 : CALL dbcsr_get_block_diag(sm_p, sm_v)
983 : END IF
984 1920 : CALL dbcsr_set(sm_v, 0.0_dp)
985 :
986 7680 : DO iatom = 1, natom
987 :
988 : CALL dbcsr_get_block_p(matrix=sm_p, &
989 : row=iatom, &
990 : col=iatom, &
991 : block=p_block, &
992 5760 : found=found)
993 :
994 5760 : IF (.NOT. ASSOCIATED(p_block)) CYCLE
995 :
996 : CALL dbcsr_get_block_p(matrix=sm_q, &
997 : row=iatom, &
998 : col=iatom, &
999 : block=q_block, &
1000 2880 : found=found)
1001 2880 : CPASSERT(ASSOCIATED(q_block))
1002 :
1003 13440 : IF (orthonormal_basis) THEN
1004 : ! S is the unit matrix
1005 0 : DO isgf = 1, SIZE(q_block, 1)
1006 0 : q_block(isgf, isgf) = p_block(isgf, isgf)
1007 : END DO
1008 : ELSE
1009 : CALL dbcsr_get_block_p(matrix=sm_s, &
1010 : row=iatom, &
1011 : col=iatom, &
1012 : block=s_block, &
1013 2880 : found=found)
1014 2880 : CPASSERT(ASSOCIATED(s_block))
1015 : ! Exploit that P and S are symmetric
1016 24960 : DO jsgf = 1, SIZE(p_block, 2)
1017 238080 : DO isgf = 1, SIZE(p_block, 1)
1018 232320 : q_block(isgf, jsgf) = p_block(isgf, jsgf)*s_block(isgf, jsgf)
1019 : END DO
1020 : END DO
1021 : END IF ! orthonormal basis set
1022 :
1023 : END DO ! next atom "iatom"
1024 :
1025 : ! E[DFT+U] = E[DFT] + E[U]
1026 : ! = E[DFT] + (U - J)*(Tr(q) - Tr(q*q))/2
1027 :
1028 : ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U]
1029 : ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j)
1030 : ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j))
1031 :
1032 : ! Loop over all atomic kinds
1033 :
1034 5760 : DO ikind = 1, nkind
1035 :
1036 : ! Load the required atomic kind data
1037 :
1038 : CALL get_atomic_kind(atomic_kind_set(ikind), &
1039 : atom_list=atom_list, &
1040 : name=atomic_kind_name, &
1041 3840 : natom=natom_of_kind)
1042 :
1043 : CALL get_qs_kind(qs_kind_set(ikind), &
1044 : dft_plus_u_atom=dft_plus_u_atom, &
1045 : l_of_dft_plus_u=lu, &
1046 : nsgf=nsgf_kind, &
1047 : basis_set=orb_basis_set, &
1048 : u_minus_j=u_minus_j, &
1049 : u_minus_j_target=u_minus_j_target, &
1050 : u_ramping=u_ramping, &
1051 : eps_u_ramping=eps_u_ramping, &
1052 : nelec=nelec, &
1053 : orbitals=orbitals, &
1054 : eps_scf=eps_scf, &
1055 : max_scf=max_scf, &
1056 3840 : smear=smear)
1057 :
1058 : ! Check, if the atoms of this atomic kind need a DFT+U correction
1059 :
1060 3840 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
1061 3840 : IF (.NOT. dft_plus_u_atom) CYCLE
1062 1920 : IF (lu < 0) CYCLE
1063 :
1064 : ! Apply U ramping if requested
1065 :
1066 1920 : IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
1067 976 : IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
1068 464 : u_minus_j = MIN(u_minus_j + u_ramping, u_minus_j_target)
1069 464 : CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
1070 : END IF
1071 976 : IF (should_output .AND. (output_unit > 0)) THEN
1072 : WRITE (UNIT=output_unit, FMT="(T3,A,3X,A,F0.3,A)") &
1073 476 : "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)), &
1074 952 : "U(eff) = ", u_minus_j*evolt, " eV"
1075 : END IF
1076 : END IF
1077 :
1078 1920 : IF (u_minus_j == 0.0_dp) CYCLE
1079 :
1080 1920 : is_plus_u_kind(ikind) = .TRUE.
1081 :
1082 : ! Load the required Gaussian basis set data
1083 :
1084 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1085 : first_sgf=first_sgf, &
1086 : l=l, &
1087 : last_sgf=last_sgf, &
1088 : nset=nset, &
1089 1920 : nshell=nshell)
1090 :
1091 : ! Count the relevant shell blocks of this atomic kind
1092 :
1093 1920 : nsb = 0
1094 5760 : DO iset = 1, nset
1095 15360 : DO ishell = 1, nshell(iset)
1096 13440 : IF (l(ishell, iset) == lu) nsb = nsb + 1
1097 : END DO
1098 : END DO
1099 :
1100 1920 : nsbsize = (2*lu + 1)
1101 1920 : n = nsb*nsbsize
1102 :
1103 7680 : ALLOCATE (q_matrix(n, n))
1104 82560 : q_matrix(:, :) = 0.0_dp
1105 :
1106 : ! Print headline if requested
1107 :
1108 1920 : IF (should_output .AND. (print_level > low_print_level)) THEN
1109 0 : IF (output_unit > 0) THEN
1110 0 : ALLOCATE (symbol(nsbsize))
1111 0 : DO m = -lu, lu
1112 0 : symbol(lu + m + 1) = sgf_symbol(0, lu, m)
1113 : END DO
1114 0 : IF (nspin > 1) THEN
1115 0 : WRITE (UNIT=spin_info, FMT="(A8,I2)") " of spin", ispin
1116 : ELSE
1117 0 : spin_info = ""
1118 : END IF
1119 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
1120 0 : "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ", ikind, &
1121 0 : ": "//TRIM(atomic_kind_name), &
1122 0 : "Atom Shell ", (ADJUSTR(symbol(i)), i=1, nsbsize), " Trace"
1123 0 : DEALLOCATE (symbol)
1124 : END IF
1125 : END IF
1126 :
1127 : ! Loop over all atoms of the current atomic kind
1128 :
1129 3840 : DO iatom = 1, natom_of_kind
1130 :
1131 1920 : atom_a = atom_list(iatom)
1132 :
1133 82560 : q_matrix(:, :) = 0.0_dp
1134 :
1135 : ! Get diagonal block
1136 :
1137 : CALL dbcsr_get_block_p(matrix=sm_q, &
1138 : row=atom_a, &
1139 : col=atom_a, &
1140 : block=q_block, &
1141 1920 : found=found)
1142 :
1143 : ! Calculate energy contribution to E(U)
1144 :
1145 1920 : IF (ASSOCIATED(q_block)) THEN
1146 :
1147 960 : i = 0
1148 2880 : DO iset = 1, nset
1149 7680 : DO ishell = 1, nshell(iset)
1150 4800 : IF (l(ishell, iset) /= lu) CYCLE
1151 9600 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1152 5760 : i = i + 1
1153 5760 : j = 0
1154 22080 : DO jset = 1, nset
1155 46080 : DO jshell = 1, nshell(jset)
1156 28800 : IF (l(jshell, jset) /= lu) CYCLE
1157 57600 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
1158 34560 : j = j + 1
1159 63360 : q_matrix(i, j) = q_block(isgf, jsgf)
1160 : END DO ! next contracted spherical Gaussian function "jsgf"
1161 : END DO ! next shell "jshell"
1162 : END DO ! next shell set "jset"
1163 : END DO ! next contracted spherical Gaussian function "isgf"
1164 : END DO ! next shell "ishell"
1165 : END DO ! next shell set "iset"
1166 :
1167 : ! Perform the requested manipulations of the (initial) orbital occupations
1168 :
1169 960 : IF (ASSOCIATED(orbitals)) THEN
1170 0 : IF ((qs_env%scf_env%iter_delta >= eps_scf) .OR. &
1171 : ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
1172 : (qs_env%scf_env%iter_count <= max_scf))) THEN
1173 0 : ALLOCATE (orb_occ(nsbsize))
1174 0 : ALLOCATE (q_eigval(n))
1175 0 : q_eigval(:) = 0.0_dp
1176 0 : ALLOCATE (q_eigvec(n, n))
1177 0 : q_eigvec(:, :) = 0.0_dp
1178 0 : norb = SIZE(orbitals)
1179 0 : CALL jacobi(q_matrix, q_eigval, q_eigvec)
1180 0 : q_matrix(:, :) = 0.0_dp
1181 0 : IF (nelec(ispin) >= 0.5_dp) THEN
1182 0 : trq = nelec(ispin)/SUM(q_eigval(1:n))
1183 0 : q_eigval(1:n) = trq*q_eigval(1:n)
1184 : END IF
1185 0 : DO isb = 1, nsb
1186 0 : trq = 0.0_dp
1187 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
1188 0 : trq = trq + q_eigval(i)
1189 : END DO
1190 0 : IF (smear) THEN
1191 0 : occ = trq/REAL(norb, KIND=dp)
1192 : ELSE
1193 0 : occ = 1.0_dp/fspin
1194 : END IF
1195 0 : orb_occ(:) = .FALSE.
1196 0 : iloc = MAXLOC(q_eigvec(:, isb*nsbsize))
1197 0 : jsb = INT((iloc(1) - 1)/nsbsize) + 1
1198 0 : i = 0
1199 0 : i0 = (jsb - 1)*nsbsize + 1
1200 0 : iorb = -1000
1201 0 : DO j = i0, jsb*nsbsize
1202 0 : i = i + 1
1203 0 : IF (i > norb) THEN
1204 0 : DO m = -lu, lu
1205 0 : IF (.NOT. orb_occ(lu + m + 1)) THEN
1206 0 : iorb = i0 + lu + m
1207 0 : orb_occ(lu + m + 1) = .TRUE.
1208 : END IF
1209 : END DO
1210 : ELSE
1211 0 : iorb = i0 + lu + orbitals(i)
1212 0 : orb_occ(lu + orbitals(i) + 1) = .TRUE.
1213 : END IF
1214 0 : CPASSERT(iorb /= -1000)
1215 0 : iloc = MAXLOC(q_eigvec(iorb, :))
1216 0 : q_eigval(iloc(1)) = MIN(occ, trq)
1217 0 : q_matrix(:, iloc(1)) = q_eigval(iloc(1))*q_eigvec(:, iloc(1)) ! backtransform left
1218 0 : trq = trq - q_eigval(iloc(1))
1219 : END DO
1220 : END DO
1221 0 : q_matrix(:, :) = MATMUL(q_matrix, TRANSPOSE(q_eigvec)) ! backtransform right
1222 0 : DEALLOCATE (orb_occ)
1223 0 : DEALLOCATE (q_eigval)
1224 0 : DEALLOCATE (q_eigvec)
1225 0 : occupation_enforced = .TRUE.
1226 : END IF
1227 : END IF ! orbitals associated
1228 :
1229 960 : trq = 0.0_dp
1230 960 : trq2 = 0.0_dp
1231 :
1232 6720 : DO i = 1, n
1233 5760 : trq = trq + q_matrix(i, i)
1234 41280 : DO j = 1, n
1235 40320 : trq2 = trq2 + q_matrix(i, j)*q_matrix(j, i)
1236 : END DO
1237 : END DO
1238 :
1239 960 : trq = fspin*trq
1240 960 : trq2 = fspin*fspin*trq2
1241 :
1242 : ! Calculate energy contribution to E(U)
1243 :
1244 960 : energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin
1245 :
1246 : ! Calculate potential V(U) = dE(U)/dq
1247 :
1248 960 : IF (.NOT. just_energy) THEN
1249 :
1250 : CALL dbcsr_get_block_p(matrix=sm_v, &
1251 : row=atom_a, &
1252 : col=atom_a, &
1253 : block=v_block, &
1254 536 : found=found)
1255 536 : CPASSERT(ASSOCIATED(v_block))
1256 :
1257 536 : i = 0
1258 1608 : DO iset = 1, nset
1259 4288 : DO ishell = 1, nshell(iset)
1260 2680 : IF (l(ishell, iset) /= lu) CYCLE
1261 5360 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1262 3216 : i = i + 1
1263 3216 : j = 0
1264 12328 : DO jset = 1, nset
1265 25728 : DO jshell = 1, nshell(jset)
1266 16080 : IF (l(jshell, jset) /= lu) CYCLE
1267 32160 : DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
1268 19296 : j = j + 1
1269 35376 : IF (isgf == jsgf) THEN
1270 3216 : v_block(isgf, isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i, i))
1271 : ELSE
1272 16080 : v_block(isgf, jsgf) = -u_minus_j*fspin*q_matrix(j, i)
1273 : END IF
1274 : END DO ! next contracted spherical Gaussian function "jsgf"
1275 : END DO ! next shell "jshell"
1276 : END DO ! next shell set "jset"
1277 : END DO ! next contracted spherical Gaussian function "isgf"
1278 : END DO ! next shell "ishell"
1279 : END DO ! next shell set "iset"
1280 :
1281 : END IF ! not just energy
1282 :
1283 : END IF ! q_block associated
1284 :
1285 : ! Consider print requests
1286 :
1287 5760 : IF (should_output .AND. (print_level > low_print_level)) THEN
1288 0 : CALL para_env%sum(q_matrix)
1289 0 : IF (output_unit > 0) THEN
1290 0 : ALLOCATE (q_work(nsb, nsbsize))
1291 0 : q_work(:, :) = 0.0_dp
1292 0 : DO isb = 1, nsb
1293 0 : j = 0
1294 0 : DO i = (isb - 1)*nsbsize + 1, isb*nsbsize
1295 0 : j = j + 1
1296 0 : q_work(isb, j) = q_matrix(i, i)
1297 : END DO
1298 : END DO
1299 0 : DO isb = 1, nsb
1300 : WRITE (UNIT=output_unit, FMT="(T3,I6,2X,I6,2X,10F8.3)") &
1301 0 : atom_a, isb, q_work(isb, :), SUM(q_work(isb, :))
1302 : END DO
1303 : WRITE (UNIT=output_unit, FMT="(T12,A,2X,10F8.3)") &
1304 0 : "Total", (SUM(q_work(:, i)), i=1, nsbsize), SUM(q_work)
1305 0 : WRITE (UNIT=output_unit, FMT="(A)") ""
1306 0 : DEALLOCATE (q_work)
1307 : IF (debug) THEN
1308 : ! Print the DFT+U occupation matrix
1309 : WRITE (UNIT=output_unit, FMT="(T9,70I10)") (i, i=1, n)
1310 : DO i = 1, n
1311 : WRITE (UNIT=output_unit, FMT="(T3,I6,70F10.6)") i, q_matrix(i, :)
1312 : END DO
1313 : ! Print the eigenvalues and eigenvectors of the occupation matrix
1314 : ALLOCATE (q_eigval(n))
1315 : q_eigval(:) = 0.0_dp
1316 : ALLOCATE (q_eigvec(n, n))
1317 : q_eigvec(:, :) = 0.0_dp
1318 : CALL jacobi(q_matrix, q_eigval, q_eigvec)
1319 : WRITE (UNIT=output_unit, FMT="(/,T9,70I10)") (i, i=1, n)
1320 : WRITE (UNIT=output_unit, FMT="(T9,71F10.6)") (q_eigval(i), i=1, n), &
1321 : SUM(q_eigval(1:n))
1322 : DO i = 1, n
1323 : WRITE (UNIT=output_unit, FMT="(T3,I6,70F10.6)") i, q_eigvec(i, :)
1324 : END DO
1325 : DEALLOCATE (q_eigval)
1326 : DEALLOCATE (q_eigvec)
1327 : END IF ! debug
1328 : END IF
1329 : IF (debug) THEN
1330 : ! Print the full atomic occupation matrix block
1331 : ALLOCATE (q_work(nsgf_kind, nsgf_kind))
1332 : q_work(:, :) = 0.0_dp
1333 : IF (ASSOCIATED(q_block)) q_work(:, :) = q_block(:, :)
1334 : CALL para_env%sum(q_work)
1335 : IF (output_unit > 0) THEN
1336 : norb = SIZE(q_work, 1)
1337 : WRITE (UNIT=output_unit, FMT="(/,T9,200I10)") (i, i=1, norb)
1338 : DO i = 1, norb
1339 : WRITE (UNIT=output_unit, FMT="(T3,I6,200F10.6)") i, q_work(i, :)
1340 : END DO
1341 : ALLOCATE (q_eigval(norb))
1342 : q_eigval(:) = 0.0_dp
1343 : ALLOCATE (q_eigvec(norb, norb))
1344 : q_eigvec(:, :) = 0.0_dp
1345 : CALL jacobi(q_work, q_eigval, q_eigvec)
1346 : WRITE (UNIT=output_unit, FMT="(/,T9,200I10)") (i, i=1, norb)
1347 : WRITE (UNIT=output_unit, FMT="(T9,201F10.6)") (q_eigval(i), i=1, norb), &
1348 : SUM(q_eigval(1:norb))
1349 : DO i = 1, norb
1350 : WRITE (UNIT=output_unit, FMT="(T3,I6,200F10.6)") i, q_eigvec(i, :)
1351 : END DO
1352 : DEALLOCATE (q_eigval)
1353 : DEALLOCATE (q_eigvec)
1354 : END IF
1355 : DEALLOCATE (q_work)
1356 : END IF ! debug
1357 : END IF ! should output
1358 :
1359 : END DO ! next atom "iatom" of atomic kind "ikind"
1360 :
1361 9600 : IF (ALLOCATED(q_matrix)) THEN
1362 1920 : DEALLOCATE (q_matrix)
1363 : END IF
1364 :
1365 : END DO ! next atomic kind "ikind"
1366 :
1367 : ! Add V(i,j)[U] to V(i,j)[DFT]
1368 :
1369 3840 : IF (ASSOCIATED(sm_h)) THEN
1370 :
1371 3216 : DO ikind = 1, nkind
1372 :
1373 2144 : IF (.NOT. is_plus_u_kind(ikind)) CYCLE
1374 :
1375 1072 : kind_a => atomic_kind_set(ikind)
1376 :
1377 : CALL get_atomic_kind(atomic_kind=kind_a, &
1378 : atom_list=atom_list, &
1379 1072 : natom=natom_of_kind)
1380 :
1381 3216 : DO iatom = 1, natom_of_kind
1382 :
1383 1072 : atom_a = atom_list(iatom)
1384 :
1385 : CALL dbcsr_get_block_p(matrix=sm_h, &
1386 : row=atom_a, &
1387 : col=atom_a, &
1388 : block=h_block, &
1389 1072 : found=found)
1390 :
1391 1072 : IF (.NOT. ASSOCIATED(h_block)) CYCLE
1392 :
1393 : CALL dbcsr_get_block_p(matrix=sm_v, &
1394 : row=atom_a, &
1395 : col=atom_a, &
1396 : block=v_block, &
1397 536 : found=found)
1398 536 : CPASSERT(ASSOCIATED(v_block))
1399 :
1400 4288 : IF (orthonormal_basis) THEN
1401 0 : DO isgf = 1, SIZE(h_block, 1)
1402 0 : h_block(isgf, isgf) = h_block(isgf, isgf) + v_block(isgf, isgf)
1403 : END DO
1404 : ELSE
1405 : CALL dbcsr_get_block_p(matrix=sm_s, &
1406 : row=atom_a, &
1407 : col=atom_a, &
1408 : block=s_block, &
1409 536 : found=found)
1410 536 : CPASSERT(ASSOCIATED(s_block))
1411 7504 : DO jsgf = 1, SIZE(h_block, 2)
1412 98624 : DO isgf = 1, SIZE(h_block, 1)
1413 97552 : h_block(isgf, jsgf) = h_block(isgf, jsgf) + v_block(isgf, jsgf)*s_block(isgf, jsgf)
1414 : END DO
1415 : END DO
1416 : END IF ! orthonormal basis set
1417 :
1418 : END DO ! next atom "iatom" of atomic kind "ikind"
1419 :
1420 : END DO ! Next atomic kind "ikind"
1421 :
1422 : END IF ! An update of the Hamiltonian matrix is requested
1423 :
1424 : END DO ! next cell image
1425 :
1426 : END DO ! next spin "ispin"
1427 :
1428 : ! Collect the energy contributions from all processes
1429 :
1430 1290 : CALL para_env%sum(energy%dft_plus_u)
1431 :
1432 1290 : IF (energy%dft_plus_u < 0.0_dp) THEN
1433 0 : IF (.NOT. occupation_enforced) THEN
1434 : CALL cp_warn(__LOCATION__, &
1435 : "DFT+U energy contribution is negative possibly due "// &
1436 0 : "to unphysical Mulliken charges!")
1437 : END IF
1438 : END IF
1439 :
1440 1290 : CALL dbcsr_deallocate_matrix(sm_q)
1441 1290 : CALL dbcsr_deallocate_matrix(sm_v)
1442 :
1443 1290 : CALL timestop(handle)
1444 :
1445 3870 : END SUBROUTINE mulliken
1446 :
1447 : ! **************************************************************************************************
1448 : !> \brief Add a DFT+U contribution to the Hamiltonian matrix\n
1449 : !> using a method based on Mulliken charges
1450 : !> \f[q_\mu = \sum_\nu \frac{1}{2}(P_{\mu\nu} S_{\nu\mu} +
1451 : !> S_{\mu\nu} P_{\nu\mu})
1452 : !> = \sum_\nu P_{\mu\nu} S_{\nu\mu}\f]
1453 : !> where \b P and \b S are the density and the
1454 : !> overlap matrix, respectively.
1455 : !> \param[in] qs_env Quickstep environment
1456 : !> \param orthonormal_basis ...
1457 : !> \param[in,out] matrix_h Hamiltonian matrices for each spin
1458 : !> \param[in,out] matrix_w Energy weighted density matrices for each spin
1459 : !> \param should_output ...
1460 : !> \param output_unit ...
1461 : !> \param print_level ...
1462 : !> \date 11.01.2008
1463 : !> \par
1464 : !> \f{eqnarray*}{
1465 : !> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}\\\
1466 : !> & = & E^{\rm DFT} + \frac{1}{2}(U - J)\sum_\mu (q_\mu - q_\mu^2)\\[1ex]
1467 : !> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
1468 : !> & = & \frac{\partial E^{\rm DFT}}
1469 : !> {\partial P_{\mu\nu}} +
1470 : !> \frac{\partial E^{\rm U}}
1471 : !> {\partial P_{\mu\nu}}\\\
1472 : !> & = & H_{\mu\nu} +
1473 : !> \frac{\partial E^{\rm U}}{\partial q_\mu}
1474 : !> \frac{\partial q_\mu}{\partial P_{\mu\nu}}\\\
1475 : !> & = & H_{\mu\nu} +
1476 : !> \frac{1}{2}(U - J)(1 - q_\mu - q_\nu) S_{\mu\nu}\\\
1477 : !> \f}
1478 : !> \author Matthias Krack (MK)
1479 : !> \version 1.0
1480 : !> \note The use of any full matrices was avoided. Thus no ScaLAPACK
1481 : !> calls are performed
1482 : ! **************************************************************************************************
1483 376 : SUBROUTINE mulliken_charges(qs_env, orthonormal_basis, matrix_h, matrix_w, &
1484 : should_output, output_unit, print_level)
1485 :
1486 : TYPE(qs_environment_type), POINTER :: qs_env
1487 : LOGICAL, INTENT(IN) :: orthonormal_basis
1488 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
1489 : POINTER :: matrix_h, matrix_w
1490 : LOGICAL, INTENT(IN) :: should_output
1491 : INTEGER, INTENT(IN) :: output_unit, print_level
1492 :
1493 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mulliken_charges'
1494 :
1495 : CHARACTER(LEN=10) :: spin_info
1496 376 : CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: symbol
1497 : CHARACTER(LEN=default_string_length) :: atomic_kind_name
1498 : INTEGER :: atom_a, handle, i, iatom, ic, ikind, isb, iset, isgf, ishell, ispin, jatom, jsgf, &
1499 : lu, m, natom, natom_of_kind, nimg, nkind, nsb, nset, nsgf, nspin, sgf
1500 376 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
1501 376 : INTEGER, DIMENSION(:), POINTER :: atom_list, nshell
1502 376 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf
1503 : LOGICAL :: dft_plus_u_atom, found, just_energy
1504 : REAL(KIND=dp) :: eps_u_ramping, fspin, q, u_minus_j, &
1505 : u_minus_j_target, u_ramping, v
1506 376 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dEdq, trps
1507 376 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: q_ii
1508 376 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, s_block, w_block
1509 376 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1510 : TYPE(dbcsr_iterator_type) :: iter
1511 376 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
1512 : TYPE(dbcsr_type), POINTER :: sm_h, sm_p, sm_s, sm_w
1513 : TYPE(dft_control_type), POINTER :: dft_control
1514 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1515 : TYPE(mp_para_env_type), POINTER :: para_env
1516 376 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1517 : TYPE(qs_energy_type), POINTER :: energy
1518 376 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1519 : TYPE(qs_rho_type), POINTER :: rho
1520 :
1521 376 : CALL timeset(routineN, handle)
1522 :
1523 376 : NULLIFY (atom_list)
1524 376 : NULLIFY (atomic_kind_set)
1525 376 : NULLIFY (qs_kind_set)
1526 376 : NULLIFY (dft_control)
1527 376 : NULLIFY (energy)
1528 376 : NULLIFY (first_sgf)
1529 376 : NULLIFY (h_block)
1530 376 : NULLIFY (matrix_p)
1531 376 : NULLIFY (matrix_s)
1532 376 : NULLIFY (l)
1533 376 : NULLIFY (last_sgf)
1534 376 : NULLIFY (nshell)
1535 376 : NULLIFY (orb_basis_set)
1536 376 : NULLIFY (p_block)
1537 376 : NULLIFY (particle_set)
1538 376 : NULLIFY (rho)
1539 376 : NULLIFY (s_block)
1540 376 : NULLIFY (sm_h)
1541 376 : NULLIFY (sm_p)
1542 376 : NULLIFY (sm_s)
1543 376 : NULLIFY (w_block)
1544 376 : NULLIFY (para_env)
1545 :
1546 : CALL get_qs_env(qs_env=qs_env, &
1547 : atomic_kind_set=atomic_kind_set, &
1548 : qs_kind_set=qs_kind_set, &
1549 : dft_control=dft_control, &
1550 : energy=energy, &
1551 : particle_set=particle_set, &
1552 : rho=rho, &
1553 376 : para_env=para_env)
1554 :
1555 376 : CPASSERT(ASSOCIATED(atomic_kind_set))
1556 376 : CPASSERT(ASSOCIATED(dft_control))
1557 376 : CPASSERT(ASSOCIATED(energy))
1558 376 : CPASSERT(ASSOCIATED(particle_set))
1559 376 : CPASSERT(ASSOCIATED(rho))
1560 :
1561 376 : IF (orthonormal_basis) THEN
1562 376 : NULLIFY (sm_s)
1563 : ELSE
1564 : ! Get overlap matrix in sparse format
1565 : CALL get_qs_env(qs_env=qs_env, &
1566 376 : matrix_s_kp=matrix_s)
1567 376 : CPASSERT(ASSOCIATED(matrix_s))
1568 : END IF
1569 :
1570 : ! Get density matrices in sparse format
1571 :
1572 376 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1573 :
1574 376 : energy%dft_plus_u = 0.0_dp
1575 :
1576 376 : nspin = dft_control%nspins
1577 376 : nimg = dft_control%nimages
1578 :
1579 376 : IF (nspin == 2) THEN
1580 : fspin = 1.0_dp
1581 : ELSE
1582 220 : fspin = 0.5_dp
1583 : END IF
1584 :
1585 : ! Get the total number of atoms, contracted spherical Gaussian basis
1586 : ! functions, and atomic kinds
1587 :
1588 376 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
1589 376 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
1590 :
1591 376 : nkind = SIZE(atomic_kind_set)
1592 :
1593 1128 : ALLOCATE (first_sgf_atom(natom))
1594 1504 : first_sgf_atom(:) = 0
1595 :
1596 : CALL get_particle_set(particle_set, qs_kind_set, &
1597 376 : first_sgf=first_sgf_atom)
1598 :
1599 1128 : ALLOCATE (trps(nsgf))
1600 9024 : trps(:) = 0.0_dp
1601 :
1602 376 : IF (PRESENT(matrix_h) .OR. PRESENT(matrix_w)) THEN
1603 604 : ALLOCATE (dEdq(nsgf))
1604 302 : just_energy = .FALSE.
1605 : ELSE
1606 : just_energy = .TRUE.
1607 : END IF
1608 :
1609 : ! Loop over all spins
1610 :
1611 908 : DO ispin = 1, nspin
1612 :
1613 10146 : IF (.NOT. just_energy) dEdq(:) = 0.0_dp
1614 :
1615 : ! Calculate Trace(P*S) assuming symmetric matrices
1616 :
1617 12768 : trps(:) = 0.0_dp
1618 :
1619 5264 : DO ic = 1, nimg
1620 4732 : IF (orthonormal_basis) THEN
1621 : NULLIFY (sm_s)
1622 : ELSE
1623 4732 : sm_s => matrix_s(1, ic)%matrix
1624 : END IF
1625 4732 : sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
1626 :
1627 4732 : CALL dbcsr_iterator_start(iter, sm_p)
1628 :
1629 18928 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1630 :
1631 14196 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
1632 :
1633 18928 : IF (orthonormal_basis) THEN
1634 :
1635 0 : IF (iatom /= jatom) CYCLE
1636 :
1637 0 : IF (ASSOCIATED(p_block)) THEN
1638 0 : sgf = first_sgf_atom(iatom)
1639 0 : DO isgf = 1, SIZE(p_block, 1)
1640 0 : trps(sgf) = trps(sgf) + p_block(isgf, isgf)
1641 0 : sgf = sgf + 1
1642 : END DO
1643 : END IF
1644 :
1645 : ELSE
1646 :
1647 : CALL dbcsr_get_block_p(matrix=sm_s, &
1648 : row=iatom, &
1649 : col=jatom, &
1650 : block=s_block, &
1651 14196 : found=found)
1652 14196 : CPASSERT(ASSOCIATED(s_block))
1653 :
1654 14196 : sgf = first_sgf_atom(jatom)
1655 104104 : DO jsgf = 1, SIZE(p_block, 2)
1656 974792 : DO isgf = 1, SIZE(p_block, 1)
1657 974792 : trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
1658 : END DO
1659 104104 : sgf = sgf + 1
1660 : END DO
1661 :
1662 14196 : IF (iatom /= jatom) THEN
1663 7098 : sgf = first_sgf_atom(iatom)
1664 80444 : DO isgf = 1, SIZE(p_block, 1)
1665 440076 : DO jsgf = 1, SIZE(p_block, 2)
1666 440076 : trps(sgf) = trps(sgf) + p_block(isgf, jsgf)*s_block(isgf, jsgf)
1667 : END DO
1668 80444 : sgf = sgf + 1
1669 : END DO
1670 : END IF
1671 :
1672 : END IF ! orthonormal basis set
1673 :
1674 : END DO ! next atom "iatom"
1675 :
1676 9996 : CALL dbcsr_iterator_stop(iter)
1677 :
1678 : END DO ! cell images
1679 :
1680 532 : CALL para_env%sum(trps)
1681 :
1682 : ! q <- Trace(PS)
1683 :
1684 : ! E[DFT+U] = E[DFT] + E[U]
1685 : ! = E[DFT] + (U - J)*(q - q**2))/2
1686 :
1687 : ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U]
1688 : ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j)
1689 : ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j))
1690 :
1691 : ! Loop over all atomic kinds
1692 :
1693 1596 : DO ikind = 1, nkind
1694 :
1695 : ! Load the required atomic kind data
1696 : CALL get_atomic_kind(atomic_kind_set(ikind), &
1697 : atom_list=atom_list, &
1698 : name=atomic_kind_name, &
1699 1064 : natom=natom_of_kind)
1700 :
1701 : CALL get_qs_kind(qs_kind_set(ikind), &
1702 : dft_plus_u_atom=dft_plus_u_atom, &
1703 : l_of_dft_plus_u=lu, &
1704 : basis_set=orb_basis_set, &
1705 : u_minus_j=u_minus_j, &
1706 : u_minus_j_target=u_minus_j_target, &
1707 : u_ramping=u_ramping, &
1708 1064 : eps_u_ramping=eps_u_ramping)
1709 :
1710 : ! Check, if this atom needs a DFT+U correction
1711 :
1712 1064 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
1713 1064 : IF (.NOT. dft_plus_u_atom) CYCLE
1714 532 : IF (lu < 0) CYCLE
1715 :
1716 : ! Apply U ramping if requested
1717 :
1718 532 : IF ((ispin == 1) .AND. (u_ramping > 0.0_dp)) THEN
1719 0 : IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN
1720 0 : u_minus_j = MIN(u_minus_j + u_ramping, u_minus_j_target)
1721 0 : CALL set_qs_kind(qs_kind_set(ikind), u_minus_j=u_minus_j)
1722 : END IF
1723 0 : IF (should_output .AND. (output_unit > 0)) THEN
1724 : WRITE (UNIT=output_unit, FMT="(T3,A,3X,A,F0.3,A)") &
1725 0 : "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)), &
1726 0 : "U(eff) = ", u_minus_j*evolt, " eV"
1727 : END IF
1728 : END IF
1729 :
1730 532 : IF (u_minus_j == 0.0_dp) CYCLE
1731 :
1732 : ! Load the required Gaussian basis set data
1733 :
1734 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1735 : first_sgf=first_sgf, &
1736 : l=l, &
1737 : last_sgf=last_sgf, &
1738 : nset=nset, &
1739 532 : nshell=nshell)
1740 :
1741 : ! Count the relevant shell blocks of this atomic kind
1742 :
1743 532 : nsb = 0
1744 1596 : DO iset = 1, nset
1745 4256 : DO ishell = 1, nshell(iset)
1746 3724 : IF (l(ishell, iset) == lu) nsb = nsb + 1
1747 : END DO
1748 : END DO
1749 :
1750 2128 : ALLOCATE (q_ii(nsb, 2*lu + 1))
1751 :
1752 : ! Print headline if requested
1753 :
1754 532 : IF (should_output .AND. (print_level > low_print_level)) THEN
1755 0 : IF (output_unit > 0) THEN
1756 0 : ALLOCATE (symbol(2*lu + 1))
1757 0 : DO m = -lu, lu
1758 0 : symbol(lu + m + 1) = sgf_symbol(0, lu, m)
1759 : END DO
1760 0 : IF (nspin > 1) THEN
1761 0 : WRITE (UNIT=spin_info, FMT="(A8,I2)") " of spin", ispin
1762 : ELSE
1763 0 : spin_info = ""
1764 : END IF
1765 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))") &
1766 0 : "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ", ikind, &
1767 0 : ": "//TRIM(atomic_kind_name), &
1768 0 : "Atom Shell ", (ADJUSTR(symbol(i)), i=1, 2*lu + 1), " Trace"
1769 0 : DEALLOCATE (symbol)
1770 : END IF
1771 : END IF
1772 :
1773 : ! Loop over all atoms of the current atomic kind
1774 :
1775 1064 : DO iatom = 1, natom_of_kind
1776 :
1777 532 : atom_a = atom_list(iatom)
1778 :
1779 5320 : q_ii(:, :) = 0.0_dp
1780 :
1781 : ! Get diagonal block
1782 :
1783 : CALL dbcsr_get_block_p(matrix=sm_p, &
1784 : row=atom_a, &
1785 : col=atom_a, &
1786 : block=p_block, &
1787 532 : found=found)
1788 :
1789 : ! Calculate E(U) and dE(U)/dq
1790 :
1791 532 : IF (ASSOCIATED(p_block)) THEN
1792 :
1793 266 : sgf = first_sgf_atom(atom_a)
1794 :
1795 266 : isb = 0
1796 798 : DO iset = 1, nset
1797 2128 : DO ishell = 1, nshell(iset)
1798 1862 : IF (l(ishell, iset) == lu) THEN
1799 532 : isb = isb + 1
1800 532 : i = 0
1801 2128 : DO isgf = first_sgf(ishell, iset), last_sgf(ishell, iset)
1802 1596 : q = fspin*trps(sgf)
1803 1596 : i = i + 1
1804 1596 : q_ii(isb, i) = q
1805 : energy%dft_plus_u = energy%dft_plus_u + &
1806 1596 : 0.5_dp*u_minus_j*(q - q**2)/fspin
1807 1596 : IF (.NOT. just_energy) THEN
1808 1254 : dEdq(sgf) = dEdq(sgf) + u_minus_j*(0.5_dp - q)
1809 : END IF
1810 2128 : sgf = sgf + 1
1811 : END DO ! next contracted spherical Gaussian function "isgf"
1812 : ELSE
1813 798 : sgf = sgf + last_sgf(ishell, iset) - first_sgf(ishell, iset) + 1
1814 : END IF ! angular momentum requested for DFT+U correction
1815 : END DO ! next shell "ishell"
1816 : END DO ! next shell set "iset"
1817 :
1818 : END IF ! this process is the owner of the sparse matrix block?
1819 :
1820 : ! Consider print requests
1821 :
1822 1596 : IF (should_output .AND. (print_level > low_print_level)) THEN
1823 0 : CALL para_env%sum(q_ii)
1824 0 : IF (output_unit > 0) THEN
1825 0 : DO isb = 1, nsb
1826 : WRITE (UNIT=output_unit, FMT="(T3,I6,2X,I6,2X,10F8.3)") &
1827 0 : atom_a, isb, q_ii(isb, :), SUM(q_ii(isb, :))
1828 : END DO
1829 : WRITE (UNIT=output_unit, FMT="(T12,A,2X,10F8.3)") &
1830 0 : "Total", (SUM(q_ii(:, i)), i=1, 2*lu + 1), SUM(q_ii)
1831 0 : WRITE (UNIT=output_unit, FMT="(A)") ""
1832 : END IF
1833 : END IF ! should output
1834 :
1835 : END DO ! next atom "iatom" of atomic kind "ikind"
1836 :
1837 2660 : IF (ALLOCATED(q_ii)) THEN
1838 532 : DEALLOCATE (q_ii)
1839 : END IF
1840 :
1841 : END DO ! next atomic kind "ikind"
1842 :
1843 532 : IF (.NOT. just_energy) THEN
1844 418 : CALL para_env%sum(dEdq)
1845 : END IF
1846 :
1847 : ! Add V(i,j)[U] to V(i,j)[DFT]
1848 :
1849 532 : IF (PRESENT(matrix_h)) THEN
1850 :
1851 4840 : DO ic = 1, nimg
1852 4460 : IF (orthonormal_basis) THEN
1853 : NULLIFY (sm_s)
1854 : ELSE
1855 4460 : sm_s => matrix_s(1, ic)%matrix
1856 : END IF
1857 4460 : sm_h => matrix_h(ispin, ic)%matrix
1858 :
1859 4460 : CALL dbcsr_iterator_start(iter, sm_h)
1860 :
1861 17840 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1862 :
1863 13380 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, h_block)
1864 :
1865 17840 : IF (orthonormal_basis) THEN
1866 :
1867 0 : IF (iatom /= jatom) CYCLE
1868 :
1869 0 : IF (ASSOCIATED(h_block)) THEN
1870 0 : sgf = first_sgf_atom(iatom)
1871 0 : DO isgf = 1, SIZE(h_block, 1)
1872 0 : h_block(isgf, isgf) = h_block(isgf, isgf) + dEdq(sgf)
1873 0 : sgf = sgf + 1
1874 : END DO
1875 : END IF
1876 :
1877 : ELSE
1878 :
1879 : ! Request katom just to check for consistent sparse matrix pattern
1880 :
1881 : CALL dbcsr_get_block_p(matrix=sm_s, &
1882 : row=iatom, &
1883 : col=jatom, &
1884 : block=s_block, &
1885 13380 : found=found)
1886 13380 : CPASSERT(ASSOCIATED(s_block))
1887 :
1888 : ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
1889 :
1890 13380 : sgf = first_sgf_atom(iatom)
1891 :
1892 133800 : DO isgf = 1, SIZE(h_block, 1)
1893 120420 : IF (dEdq(sgf) /= 0.0_dp) THEN
1894 40140 : v = 0.5_dp*dEdq(sgf)
1895 347880 : DO jsgf = 1, SIZE(h_block, 2)
1896 347880 : h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
1897 : END DO
1898 : END IF
1899 133800 : sgf = sgf + 1
1900 : END DO
1901 :
1902 13380 : sgf = first_sgf_atom(jatom)
1903 :
1904 98120 : DO jsgf = 1, SIZE(h_block, 2)
1905 84740 : IF (dEdq(sgf) /= 0.0_dp) THEN
1906 13380 : v = 0.5_dp*dEdq(sgf)
1907 187320 : DO isgf = 1, SIZE(h_block, 1)
1908 187320 : h_block(isgf, jsgf) = h_block(isgf, jsgf) + v*s_block(isgf, jsgf)
1909 : END DO
1910 : END IF
1911 98120 : sgf = sgf + 1
1912 : END DO
1913 :
1914 : END IF ! orthonormal basis set
1915 :
1916 : END DO ! Next atom "iatom"
1917 :
1918 9300 : CALL dbcsr_iterator_stop(iter)
1919 :
1920 : END DO
1921 :
1922 : END IF ! An update of the Hamiltonian matrix is requested
1923 :
1924 : ! Calculate the contribution (non-Pulay part) to the derivatives
1925 : ! w.r.t. the nuclear positions, which requires an update of the
1926 : ! energy weighted density W.
1927 :
1928 908 : IF (PRESENT(matrix_w) .AND. (.NOT. orthonormal_basis)) THEN
1929 :
1930 196 : DO ic = 1, nimg
1931 158 : sm_s => matrix_s(1, ic)%matrix
1932 158 : sm_p => matrix_p(ispin, ic)%matrix
1933 158 : sm_w => matrix_w(ispin, ic)%matrix
1934 :
1935 158 : CALL dbcsr_iterator_start(iter, sm_p)
1936 :
1937 632 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1938 :
1939 474 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, p_block)
1940 :
1941 : ! Skip the diagonal blocks of the W matrix
1942 :
1943 474 : IF (iatom == jatom) CYCLE
1944 :
1945 : ! Request katom just to check for consistent sparse matrix patterns
1946 :
1947 : CALL dbcsr_get_block_p(matrix=sm_w, &
1948 : row=iatom, &
1949 : col=jatom, &
1950 : block=w_block, &
1951 237 : found=found)
1952 237 : CPASSERT(ASSOCIATED(w_block))
1953 :
1954 : ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation
1955 :
1956 237 : sgf = first_sgf_atom(iatom)
1957 :
1958 2686 : DO isgf = 1, SIZE(w_block, 1)
1959 2449 : IF (dEdq(sgf) /= 0.0_dp) THEN
1960 948 : v = -0.5_dp*dEdq(sgf)
1961 5688 : DO jsgf = 1, SIZE(w_block, 2)
1962 5688 : w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
1963 : END DO
1964 : END IF
1965 2686 : sgf = sgf + 1
1966 : END DO
1967 :
1968 237 : sgf = first_sgf_atom(jatom)
1969 :
1970 1580 : DO jsgf = 1, SIZE(w_block, 2)
1971 1185 : IF (dEdq(sgf) /= 0.0_dp) THEN
1972 0 : v = -0.5_dp*dEdq(sgf)
1973 0 : DO isgf = 1, SIZE(w_block, 1)
1974 0 : w_block(isgf, jsgf) = w_block(isgf, jsgf) + v*p_block(isgf, jsgf)
1975 : END DO
1976 : END IF
1977 1659 : sgf = sgf + 1
1978 : END DO
1979 :
1980 : END DO ! next block node "jatom"
1981 :
1982 354 : CALL dbcsr_iterator_stop(iter)
1983 :
1984 : END DO
1985 :
1986 : END IF ! W matrix update requested
1987 :
1988 : END DO ! next spin "ispin"
1989 :
1990 : ! Collect the energy contributions from all processes
1991 :
1992 376 : CALL para_env%sum(energy%dft_plus_u)
1993 :
1994 376 : IF (energy%dft_plus_u < 0.0_dp) &
1995 : CALL cp_warn(__LOCATION__, &
1996 : "DFT+U energy contribution is negative possibly due "// &
1997 0 : "to unphysical Mulliken charges!")
1998 :
1999 : ! Release local work storage
2000 :
2001 376 : IF (ALLOCATED(first_sgf_atom)) THEN
2002 376 : DEALLOCATE (first_sgf_atom)
2003 : END IF
2004 :
2005 376 : IF (ALLOCATED(trps)) THEN
2006 376 : DEALLOCATE (trps)
2007 : END IF
2008 :
2009 376 : IF (ALLOCATED(dEdq)) THEN
2010 302 : DEALLOCATE (dEdq)
2011 : END IF
2012 :
2013 376 : CALL timestop(handle)
2014 :
2015 1128 : END SUBROUTINE mulliken_charges
2016 :
2017 : END MODULE dft_plus_u
|
