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