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