Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Calculation of KS matrix in xTB
      10             : !>        Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov
      11             : !>                   JCTC 13, 1989-2009, (2017)
      12             : !>                   DOI: 10.1021/acs.jctc.7b00118
      13             : !> \author JGH
      14             : ! **************************************************************************************************
      15             : MODULE xtb_ks_matrix
      16             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      17             :                                               get_atomic_kind
      18             :    USE cp_control_types,                ONLY: dft_control_type
      19             :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      20             :                                               dbcsr_copy,&
      21             :                                               dbcsr_dot,&
      22             :                                               dbcsr_multiply,&
      23             :                                               dbcsr_p_type,&
      24             :                                               dbcsr_type
      25             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      26             :                                               cp_logger_get_default_io_unit,&
      27             :                                               cp_logger_type
      28             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      29             :                                               cp_print_key_unit_nr
      30             :    USE efield_tb_methods,               ONLY: efield_tb_matrix
      31             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      32             :                                               section_vals_type
      33             :    USE kinds,                           ONLY: dp
      34             :    USE message_passing,                 ONLY: mp_para_env_type
      35             :    USE mulliken,                        ONLY: ao_charges
      36             :    USE particle_types,                  ONLY: particle_type
      37             :    USE qs_charge_mixing,                ONLY: charge_mixing
      38             :    USE qs_energy_types,                 ONLY: qs_energy_type
      39             :    USE qs_environment_types,            ONLY: get_qs_env,&
      40             :                                               qs_environment_type
      41             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      42             :                                               get_qs_kind_set,&
      43             :                                               qs_kind_type
      44             :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      45             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      46             :                                               mo_set_type
      47             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      48             :                                               qs_rho_type
      49             :    USE qs_scf_types,                    ONLY: qs_scf_env_type
      50             :    USE xtb_coulomb,                     ONLY: build_xtb_coulomb
      51             :    USE xtb_types,                       ONLY: get_xtb_atom_param,&
      52             :                                               xtb_atom_type
      53             : #include "./base/base_uses.f90"
      54             : 
      55             :    IMPLICIT NONE
      56             : 
      57             :    PRIVATE
      58             : 
      59             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_ks_matrix'
      60             : 
      61             :    PUBLIC :: build_xtb_ks_matrix
      62             : 
      63             : CONTAINS
      64             : 
      65             : ! **************************************************************************************************
      66             : !> \brief ...
      67             : !> \param qs_env ...
      68             : !> \param calculate_forces ...
      69             : !> \param just_energy ...
      70             : !> \param ext_ks_matrix ...
      71             : ! **************************************************************************************************
      72       25240 :    SUBROUTINE build_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
      73             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      74             :       LOGICAL, INTENT(in)                                :: calculate_forces, just_energy
      75             :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
      76             :          POINTER                                         :: ext_ks_matrix
      77             : 
      78             :       INTEGER                                            :: gfn_type
      79             :       TYPE(dft_control_type), POINTER                    :: dft_control
      80             : 
      81       25240 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
      82       25240 :       gfn_type = dft_control%qs_control%xtb_control%gfn_type
      83             : 
      84        2184 :       SELECT CASE (gfn_type)
      85             :       CASE (0)
      86        2184 :          CPASSERT(.NOT. PRESENT(ext_ks_matrix))
      87        2184 :          CALL build_gfn0_xtb_ks_matrix(qs_env, calculate_forces, just_energy)
      88             :       CASE (1)
      89       23056 :          CALL build_gfn1_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
      90             :       CASE (2)
      91           0 :          CPABORT("gfn_type = 2 not yet available")
      92             :       CASE DEFAULT
      93       25240 :          CPABORT("Unknown gfn_type")
      94             :       END SELECT
      95             : 
      96       25240 :    END SUBROUTINE build_xtb_ks_matrix
      97             : 
      98             : ! **************************************************************************************************
      99             : !> \brief ...
     100             : !> \param qs_env ...
     101             : !> \param calculate_forces ...
     102             : !> \param just_energy ...
     103             : ! **************************************************************************************************
     104        2184 :    SUBROUTINE build_gfn0_xtb_ks_matrix(qs_env, calculate_forces, just_energy)
     105             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     106             :       LOGICAL, INTENT(in)                                :: calculate_forces, just_energy
     107             : 
     108             :       CHARACTER(len=*), PARAMETER :: routineN = 'build_gfn0_xtb_ks_matrix'
     109             : 
     110             :       INTEGER                                            :: handle, img, iounit, ispin, natom, nimg, &
     111             :                                                             nspins
     112             :       REAL(KIND=dp)                                      :: pc_ener, qmmm_el
     113        2184 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     114             :       TYPE(cp_logger_type), POINTER                      :: logger
     115        2184 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p1, mo_derivs
     116        2184 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: ks_matrix, matrix_h
     117             :       TYPE(dbcsr_type), POINTER                          :: mo_coeff
     118             :       TYPE(dft_control_type), POINTER                    :: dft_control
     119             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     120             :       TYPE(qs_energy_type), POINTER                      :: energy
     121        2184 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     122             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     123             :       TYPE(qs_rho_type), POINTER                         :: rho
     124             :       TYPE(section_vals_type), POINTER                   :: scf_section
     125             : 
     126        2184 :       CALL timeset(routineN, handle)
     127             : 
     128             :       MARK_USED(calculate_forces)
     129             : 
     130        2184 :       NULLIFY (dft_control, logger, scf_section, ks_env, ks_matrix, rho, &
     131        2184 :                energy)
     132        2184 :       CPASSERT(ASSOCIATED(qs_env))
     133             : 
     134        2184 :       logger => cp_get_default_logger()
     135        2184 :       iounit = cp_logger_get_default_io_unit(logger)
     136             : 
     137             :       CALL get_qs_env(qs_env, &
     138             :                       dft_control=dft_control, &
     139             :                       atomic_kind_set=atomic_kind_set, &
     140             :                       qs_kind_set=qs_kind_set, &
     141             :                       matrix_h_kp=matrix_h, &
     142             :                       para_env=para_env, &
     143             :                       ks_env=ks_env, &
     144             :                       matrix_ks_kp=ks_matrix, &
     145        2184 :                       energy=energy)
     146             : 
     147        2184 :       energy%hartree = 0.0_dp
     148        2184 :       energy%qmmm_el = 0.0_dp
     149             : 
     150        2184 :       scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
     151        2184 :       nspins = dft_control%nspins
     152        2184 :       nimg = dft_control%nimages
     153        2184 :       CPASSERT(ASSOCIATED(matrix_h))
     154        6552 :       CPASSERT(SIZE(ks_matrix) > 0)
     155             : 
     156        4696 :       DO ispin = 1, nspins
     157       60136 :          DO img = 1, nimg
     158             :             ! copy the core matrix into the fock matrix
     159       57952 :             CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
     160             :          END DO
     161             :       END DO
     162             : 
     163        2184 :       IF (qs_env%qmmm) THEN
     164           0 :          CPABORT("gfn0 QMMM NYA")
     165           0 :          CALL get_qs_env(qs_env=qs_env, rho=rho, natom=natom)
     166           0 :          CPASSERT(SIZE(ks_matrix, 2) == 1)
     167           0 :          DO ispin = 1, nspins
     168             :             ! If QM/MM sumup the 1el Hamiltonian
     169             :             CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
     170           0 :                            1.0_dp, 1.0_dp)
     171           0 :             CALL qs_rho_get(rho, rho_ao=matrix_p1)
     172             :             ! Compute QM/MM Energy
     173             :             CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
     174           0 :                            matrix_p1(ispin)%matrix, qmmm_el)
     175           0 :             energy%qmmm_el = energy%qmmm_el + qmmm_el
     176             :          END DO
     177           0 :          pc_ener = qs_env%ks_qmmm_env%pc_ener
     178           0 :          energy%qmmm_el = energy%qmmm_el + pc_ener
     179             :       END IF
     180             : 
     181             :       energy%total = energy%core + energy%eeq + energy%efield + energy%qmmm_el + &
     182        2184 :                      energy%repulsive + energy%dispersion + energy%kTS
     183             : 
     184             :       iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
     185        2184 :                                     extension=".scfLog")
     186        2184 :       IF (iounit > 0) THEN
     187             :          WRITE (UNIT=iounit, FMT="(/,(T9,A,T60,F20.10))") &
     188           0 :             "Repulsive pair potential energy:               ", energy%repulsive, &
     189           0 :             "SRB Correction energy:                         ", energy%srb, &
     190           0 :             "Zeroth order Hamiltonian energy:               ", energy%core, &
     191           0 :             "Charge equilibration energy:                   ", energy%eeq, &
     192           0 :             "London dispersion energy:                      ", energy%dispersion
     193           0 :          IF (dft_control%qs_control%xtb_control%do_nonbonded) &
     194             :             WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
     195           0 :             "Correction for nonbonded interactions:         ", energy%xtb_nonbonded
     196           0 :          IF (ABS(energy%efield) > 1.e-12_dp) THEN
     197             :             WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
     198           0 :                "Electric field interaction energy:             ", energy%efield
     199             :          END IF
     200           0 :          IF (qs_env%qmmm) THEN
     201             :             WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
     202           0 :                "QM/MM Electrostatic energy:                    ", energy%qmmm_el
     203             :          END IF
     204             :       END IF
     205             :       CALL cp_print_key_finished_output(iounit, logger, scf_section, &
     206        2184 :                                         "PRINT%DETAILED_ENERGY")
     207             :       ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
     208        2184 :       IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
     209           0 :          CPASSERT(SIZE(ks_matrix, 2) == 1)
     210             :          BLOCK
     211           0 :             TYPE(mo_set_type), DIMENSION(:), POINTER         :: mo_array
     212           0 :             CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
     213           0 :             DO ispin = 1, SIZE(mo_derivs)
     214           0 :                CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
     215           0 :                IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
     216           0 :                   CPABORT("")
     217             :                END IF
     218             :                CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
     219           0 :                                    0.0_dp, mo_derivs(ispin)%matrix)
     220             :             END DO
     221             :          END BLOCK
     222             :       END IF
     223             : 
     224        2184 :       CALL timestop(handle)
     225             : 
     226        2184 :    END SUBROUTINE build_gfn0_xtb_ks_matrix
     227             : 
     228             : ! **************************************************************************************************
     229             : !> \brief ...
     230             : !> \param qs_env ...
     231             : !> \param calculate_forces ...
     232             : !> \param just_energy ...
     233             : !> \param ext_ks_matrix ...
     234             : ! **************************************************************************************************
     235       23056 :    SUBROUTINE build_gfn1_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
     236             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     237             :       LOGICAL, INTENT(in)                                :: calculate_forces, just_energy
     238             :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
     239             :          POINTER                                         :: ext_ks_matrix
     240             : 
     241             :       CHARACTER(len=*), PARAMETER :: routineN = 'build_gfn1_xtb_ks_matrix'
     242             : 
     243             :       INTEGER                                            :: atom_a, handle, iatom, ikind, img, &
     244             :                                                             iounit, is, ispin, na, natom, natorb, &
     245             :                                                             nimg, nkind, ns, nsgf, nspins
     246             :       INTEGER, DIMENSION(25)                             :: lao
     247             :       INTEGER, DIMENSION(5)                              :: occ
     248             :       LOGICAL                                            :: do_efield, pass_check
     249             :       REAL(KIND=dp)                                      :: achrg, chmax, pc_ener, qmmm_el
     250       23056 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: mcharge
     251       23056 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: aocg, charges
     252       23056 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     253             :       TYPE(cp_logger_type), POINTER                      :: logger
     254       23056 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p1, mo_derivs, p_matrix
     255       23056 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: ks_matrix, matrix_h, matrix_p, matrix_s
     256             :       TYPE(dbcsr_type), POINTER                          :: mo_coeff, s_matrix
     257             :       TYPE(dft_control_type), POINTER                    :: dft_control
     258             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     259       23056 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     260             :       TYPE(qs_energy_type), POINTER                      :: energy
     261       23056 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     262             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     263             :       TYPE(qs_rho_type), POINTER                         :: rho
     264             :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     265             :       TYPE(section_vals_type), POINTER                   :: scf_section
     266             :       TYPE(xtb_atom_type), POINTER                       :: xtb_kind
     267             : 
     268       23056 :       CALL timeset(routineN, handle)
     269             : 
     270       23056 :       NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, &
     271       23056 :                ks_matrix, rho, energy)
     272       23056 :       CPASSERT(ASSOCIATED(qs_env))
     273             : 
     274       23056 :       logger => cp_get_default_logger()
     275       23056 :       iounit = cp_logger_get_default_io_unit(logger)
     276             : 
     277             :       CALL get_qs_env(qs_env, &
     278             :                       dft_control=dft_control, &
     279             :                       atomic_kind_set=atomic_kind_set, &
     280             :                       qs_kind_set=qs_kind_set, &
     281             :                       matrix_h_kp=matrix_h, &
     282             :                       para_env=para_env, &
     283             :                       ks_env=ks_env, &
     284             :                       matrix_ks_kp=ks_matrix, &
     285             :                       rho=rho, &
     286       23056 :                       energy=energy)
     287             : 
     288       23056 :       IF (PRESENT(ext_ks_matrix)) THEN
     289             :          ! remap pointer to allow for non-kpoint external ks matrix
     290             :          ! ext_ks_matrix is used in linear response code
     291          16 :          ns = SIZE(ext_ks_matrix)
     292          16 :          ks_matrix(1:ns, 1:1) => ext_ks_matrix(1:ns)
     293             :       END IF
     294             : 
     295       23056 :       energy%hartree = 0.0_dp
     296       23056 :       energy%qmmm_el = 0.0_dp
     297       23056 :       energy%efield = 0.0_dp
     298             : 
     299       23056 :       scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
     300       23056 :       nspins = dft_control%nspins
     301       23056 :       nimg = dft_control%nimages
     302       23056 :       CPASSERT(ASSOCIATED(matrix_h))
     303       23056 :       CPASSERT(ASSOCIATED(rho))
     304       69168 :       CPASSERT(SIZE(ks_matrix) > 0)
     305             : 
     306       48264 :       DO ispin = 1, nspins
     307      434282 :          DO img = 1, nimg
     308             :             ! copy the core matrix into the fock matrix
     309      411226 :             CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
     310             :          END DO
     311             :       END DO
     312             : 
     313       23056 :       IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
     314             :           dft_control%apply_efield_field) THEN
     315             :          do_efield = .TRUE.
     316             :       ELSE
     317       17160 :          do_efield = .FALSE.
     318             :       END IF
     319             : 
     320       23056 :       IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
     321             :          ! Mulliken charges
     322       21174 :          CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, matrix_s_kp=matrix_s)
     323       21174 :          CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
     324       21174 :          natom = SIZE(particle_set)
     325      105870 :          ALLOCATE (mcharge(natom), charges(natom, 5))
     326     1148204 :          charges = 0.0_dp
     327       21174 :          nkind = SIZE(atomic_kind_set)
     328       21174 :          CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
     329       84696 :          ALLOCATE (aocg(nsgf, natom))
     330     1093224 :          aocg = 0.0_dp
     331       21174 :          IF (nimg > 1) THEN
     332        2644 :             CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
     333             :          ELSE
     334       18530 :             p_matrix => matrix_p(:, 1)
     335       18530 :             s_matrix => matrix_s(1, 1)%matrix
     336       18530 :             CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
     337             :          END IF
     338       77096 :          DO ikind = 1, nkind
     339       55922 :             CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
     340       55922 :             CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
     341       55922 :             CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
     342      337250 :             DO iatom = 1, na
     343      204232 :                atom_a = atomic_kind_set(ikind)%atom_list(iatom)
     344     1225392 :                charges(atom_a, :) = REAL(occ(:), KIND=dp)
     345      911216 :                DO is = 1, natorb
     346      651062 :                   ns = lao(is) + 1
     347      855294 :                   charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
     348             :                END DO
     349             :             END DO
     350             :          END DO
     351       21174 :          DEALLOCATE (aocg)
     352             :          ! charge mixing
     353       21174 :          IF (dft_control%qs_control%do_ls_scf) THEN
     354             :             !
     355             :          ELSE
     356       20962 :             CALL get_qs_env(qs_env=qs_env, scf_env=scf_env)
     357             :             CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
     358       20962 :                                charges, para_env, scf_env%iter_count)
     359             :          END IF
     360             : 
     361      246580 :          DO iatom = 1, natom
     362     1248448 :             mcharge(iatom) = SUM(charges(iatom, :))
     363             :          END DO
     364             :       END IF
     365             : 
     366       23056 :       IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
     367             :          CALL build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
     368       21174 :                                 calculate_forces, just_energy)
     369             :       END IF
     370             : 
     371       23056 :       IF (do_efield) THEN
     372        5896 :          CALL efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
     373             :       END IF
     374             : 
     375       23056 :       IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
     376       21174 :          IF (dft_control%qs_control%xtb_control%check_atomic_charges) THEN
     377       20316 :             pass_check = .TRUE.
     378       73540 :             DO ikind = 1, nkind
     379       53224 :                CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
     380       53224 :                CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
     381       53224 :                CALL get_xtb_atom_param(xtb_kind, chmax=chmax)
     382      326672 :                DO iatom = 1, na
     383      199908 :                   atom_a = atomic_kind_set(ikind)%atom_list(iatom)
     384      199908 :                   achrg = mcharge(atom_a)
     385      253132 :                   IF (ABS(achrg) > chmax) THEN
     386           0 :                      IF (iounit > 0) THEN
     387           0 :                         WRITE (iounit, "(A,A,I3,I6,A,F4.2,A,F6.2)") " Charge outside chemical range:", &
     388           0 :                            "  Kind Atom=", ikind, atom_a, "   Limit=", chmax, "  Charge=", achrg
     389             :                      END IF
     390             :                      pass_check = .FALSE.
     391             :                   END IF
     392             :                END DO
     393             :             END DO
     394       20316 :             IF (.NOT. pass_check) THEN
     395             :                CALL cp_warn(__LOCATION__, "Atomic charges outside chemical range were detected."// &
     396             :                             " Switch-off CHECK_ATOMIC_CHARGES keyword in the &xTB section"// &
     397           0 :                             " if you want to force to continue the calculation.")
     398           0 :                CPABORT("xTB Charges")
     399             :             END IF
     400             :          END IF
     401             :       END IF
     402             : 
     403       23056 :       IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
     404       21174 :          DEALLOCATE (mcharge, charges)
     405             :       END IF
     406             : 
     407       23056 :       IF (qs_env%qmmm) THEN
     408        5146 :          CPASSERT(SIZE(ks_matrix, 2) == 1)
     409       10292 :          DO ispin = 1, nspins
     410             :             ! If QM/MM sumup the 1el Hamiltonian
     411             :             CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
     412        5146 :                            1.0_dp, 1.0_dp)
     413        5146 :             CALL qs_rho_get(rho, rho_ao=matrix_p1)
     414             :             ! Compute QM/MM Energy
     415             :             CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
     416        5146 :                            matrix_p1(ispin)%matrix, qmmm_el)
     417       10292 :             energy%qmmm_el = energy%qmmm_el + qmmm_el
     418             :          END DO
     419        5146 :          pc_ener = qs_env%ks_qmmm_env%pc_ener
     420        5146 :          energy%qmmm_el = energy%qmmm_el + pc_ener
     421             :       END IF
     422             : 
     423             :       energy%total = energy%core + energy%hartree + energy%efield + energy%qmmm_el + &
     424       23056 :                      energy%repulsive + energy%dispersion + energy%dftb3 + energy%kTS
     425             : 
     426             :       iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
     427       23056 :                                     extension=".scfLog")
     428       23056 :       IF (iounit > 0) THEN
     429             :          WRITE (UNIT=iounit, FMT="(/,(T9,A,T60,F20.10))") &
     430           0 :             "Repulsive pair potential energy:               ", energy%repulsive, &
     431           0 :             "Zeroth order Hamiltonian energy:               ", energy%core, &
     432           0 :             "Charge fluctuation energy:                     ", energy%hartree, &
     433           0 :             "London dispersion energy:                      ", energy%dispersion
     434           0 :          IF (dft_control%qs_control%xtb_control%xb_interaction) &
     435             :             WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
     436           0 :             "Correction for halogen bonding:                ", energy%xtb_xb_inter
     437           0 :          IF (dft_control%qs_control%xtb_control%do_nonbonded) &
     438             :             WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
     439           0 :             "Correction for nonbonded interactions:         ", energy%xtb_nonbonded
     440           0 :          IF (ABS(energy%efield) > 1.e-12_dp) THEN
     441             :             WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
     442           0 :                "Electric field interaction energy:             ", energy%efield
     443             :          END IF
     444             :          WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
     445           0 :             "DFTB3 3rd Order Energy Correction              ", energy%dftb3
     446           0 :          IF (qs_env%qmmm) THEN
     447             :             WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
     448           0 :                "QM/MM Electrostatic energy:                    ", energy%qmmm_el
     449             :          END IF
     450             :       END IF
     451             :       CALL cp_print_key_finished_output(iounit, logger, scf_section, &
     452       23056 :                                         "PRINT%DETAILED_ENERGY")
     453             :       ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
     454       23056 :       IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
     455        7238 :          CPASSERT(SIZE(ks_matrix, 2) == 1)
     456             :          BLOCK
     457        7238 :             TYPE(mo_set_type), DIMENSION(:), POINTER         :: mo_array
     458        7238 :             CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
     459       14524 :             DO ispin = 1, SIZE(mo_derivs)
     460        7286 :                CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
     461        7286 :                IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
     462           0 :                   CPABORT("")
     463             :                END IF
     464             :                CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
     465       14524 :                                    0.0_dp, mo_derivs(ispin)%matrix)
     466             :             END DO
     467             :          END BLOCK
     468             :       END IF
     469             : 
     470       23056 :       CALL timestop(handle)
     471             : 
     472       46112 :    END SUBROUTINE build_gfn1_xtb_ks_matrix
     473             : 
     474             : END MODULE xtb_ks_matrix
     475             :