LCOV - code coverage report
Current view: top level - src - dft_plus_u.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:3130539) Lines: 607 767 79.1 %
Date: 2025-05-14 06:57:18 Functions: 4 4 100.0 %

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

Generated by: LCOV version 1.15