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