Fix for tblite to not use OT optimizer (#4804)

src/qs_environment.F

@@ -94,10 +94,10 @@ MODULE qs_environment
         do_method_mndod, do_method_ofgpw, do_method_pdg, do_method_pm3, do_method_pm6, &
         do_method_pm6fm, do_method_pnnl, do_method_rigpw, do_method_rm1, do_method_xtb, &
         do_qmmm_gauss, do_qmmm_swave, general_roks, hden_atomic, kg_tnadd_embed_ri, rel_none, &
-        rel_trans_atom, vdw_pairpot_dftd2, vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, &
-        vdw_pairpot_dftd4, wfi_aspc_nr, wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, &
-        wfi_ps_method_nr, wfi_use_guess_method_nr, xc_vdw_fun_none, xc_vdw_fun_nonloc, &
-        xc_vdw_fun_pairpot, xtb_vdw_type_d3, xtb_vdw_type_d4, xtb_vdw_type_none
+        rel_trans_atom, smear_fermi_dirac, vdw_pairpot_dftd2, vdw_pairpot_dftd3, &
+        vdw_pairpot_dftd3bj, vdw_pairpot_dftd4, wfi_aspc_nr, wfi_linear_ps_method_nr, &
+        wfi_linear_wf_method_nr, wfi_ps_method_nr, wfi_use_guess_method_nr, xc_vdw_fun_none, &
+        xc_vdw_fun_nonloc, xc_vdw_fun_pairpot, xtb_vdw_type_d3, xtb_vdw_type_d4, xtb_vdw_type_none
    USE input_section_types,             ONLY: section_vals_get,&
                                               section_vals_get_subs_vals,&
                                               section_vals_type,&
@@ -132,6 +132,7 @@ MODULE qs_environment
                                               write_qs_particle_coordinates,&
                                               write_structure_data
    USE particle_types,                  ONLY: particle_type
+   USE physcon,                         ONLY: kelvin
    USE pw_env_types,                    ONLY: pw_env_type
    USE qmmm_types_low,                  ONLY: qmmm_env_qm_type
    USE qs_basis_rotation_methods,       ONLY: qs_basis_rotation
@@ -1260,6 +1261,15 @@ CONTAINS
       END IF
 
       ! XXXXXXX this is backwards XXXXXXXX
+      IF (dft_control%qs_control%xtb_control%do_tblite) THEN
+         IF (.NOT. scf_control%smear%do_smear) THEN
+            ! set tblite default smearing
+            scf_control%smear%do_smear = .TRUE.
+            scf_control%smear%method = smear_fermi_dirac
+            scf_control%smear%electronic_temperature = 300._dp/kelvin
+            scf_control%smear%eps_fermi_dirac = 1.E-6_dp
+         END IF
+      END IF
       dft_control%smear = scf_control%smear%do_smear
 
       ! Periodic efield needs equal occupation and orbital gradients
@@ -1690,7 +1700,9 @@ CONTAINS
       ! However, if we try to proceed alpha orbitals first, this leads us n_mo(1:2) = (10,8)
       ! due to the following assignment instruction above:
       !   IF (scf_control%added_mos(2) > 0) THEN ... ELSE; n_mo_add = scf_control%added_mos(1); END IF
-      IF (scf_control%added_mos(1) < 0) THEN
+      IF (dft_control%qs_control%xtb_control%do_tblite) THEN
+         scf_control%added_mos(1) = n_ao - n_mo(1)  ! tblite needs all MO's
+      ELSEIF (scf_control%added_mos(1) < 0) THEN
          scf_control%added_mos(1) = n_ao - n_mo(1)  ! use all available MOs
       ELSEIF (scf_control%added_mos(1) > n_ao - n_mo(1)) THEN
          CALL cp_warn(__LOCATION__, &

src/qs_ks_methods.F

@@ -125,7 +125,6 @@ MODULE qs_ks_methods
    USE qs_rho_types,                    ONLY: qs_rho_get,&
                                               qs_rho_type
    USE qs_sccs,                         ONLY: sccs
-   USE qs_scf_types,                    ONLY: qs_scf_env_type
    USE qs_vxc,                          ONLY: qs_vxc_create
    USE qs_vxc_atom,                     ONLY: calculate_vxc_atom
    USE rtp_admm_methods,                ONLY: rtp_admm_calc_rho_aux,&
@@ -146,7 +145,7 @@ MODULE qs_ks_methods
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_methods'
 
    PUBLIC :: calc_rho_tot_gspace, qs_ks_update_qs_env, qs_ks_build_kohn_sham_matrix, &
-             qs_ks_allocate_basics, evaluate_core_matrix_p_mix_new, evaluate_core_matrix_traces
+             qs_ks_allocate_basics, evaluate_core_matrix_traces
 
 CONTAINS
 
@@ -1182,10 +1181,13 @@ CONTAINS
 ! **************************************************************************************************
 !> \brief Calculates the traces of the core matrices and the density matrix.
 !> \param qs_env ...
+!> \param rho_ao_ext ...
 !> \author Ole Schuett
 ! **************************************************************************************************
-   SUBROUTINE evaluate_core_matrix_traces(qs_env)
+   SUBROUTINE evaluate_core_matrix_traces(qs_env, rho_ao_ext)
       TYPE(qs_environment_type), POINTER                 :: qs_env
+      TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
+         POINTER                                         :: rho_ao_ext
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'evaluate_core_matrix_traces'
 
@@ -1206,7 +1208,11 @@ CONTAINS
                       kinetic_kp=matrixkp_t, &
                       matrix_h_kp=matrixkp_h)
 
-      CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
+      IF (PRESENT(rho_ao_ext)) THEN
+         rho_ao_kp => rho_ao_ext
+      ELSE
+         CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
+      END IF
 
       CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy%core, dft_control%nspins)
 
@@ -1227,51 +1233,6 @@ CONTAINS
       CALL timestop(handle)
    END SUBROUTINE evaluate_core_matrix_traces
 
-! **************************************************************************************************
-!> \brief Calculates the traces of the core matrices and the density matrix.
-!> \param qs_env ...
-!> \author Johann Pototschnig
-! **************************************************************************************************
-   SUBROUTINE evaluate_core_matrix_p_mix_new(qs_env)
-      TYPE(qs_environment_type), POINTER                 :: qs_env
-
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'evaluate_core_matrix_p_mix_new'
-
-      INTEGER                                            :: handle
-      REAL(KIND=dp)                                      :: energy_core_im
-      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrixkp_h, matrixkp_t
-      TYPE(dft_control_type), POINTER                    :: dft_control
-      TYPE(qs_energy_type), POINTER                      :: energy
-      TYPE(qs_scf_env_type), POINTER                     :: scf_env
-
-      CALL timeset(routineN, handle)
-      NULLIFY (energy, dft_control, matrixkp_t, matrixkp_h)
-
-      CALL get_qs_env(qs_env, &
-                      energy=energy, &
-                      dft_control=dft_control, &
-                      kinetic_kp=matrixkp_t, &
-                      matrix_h_kp=matrixkp_h, &
-                      scf_env=scf_env)
-
-      CALL calculate_ptrace(matrixkp_h, scf_env%p_mix_new, energy%core, dft_control%nspins)
-
-      ! Add the imaginary part in the RTP case
-      IF (qs_env%run_rtp) THEN
-         IF (dft_control%rtp_control%velocity_gauge) THEN
-            CALL get_qs_env(qs_env, matrix_h_im_kp=matrixkp_h)
-            CALL calculate_ptrace(matrixkp_h, scf_env%p_mix_new, energy_core_im, dft_control%nspins)
-            energy%core = energy%core - energy_core_im
-         END IF
-      END IF
-
-      ! kinetic energy
-      IF (ASSOCIATED(matrixkp_t)) &
-         CALL calculate_ptrace(matrixkp_t, scf_env%p_mix_new, energy%kinetic, dft_control%nspins)
-
-      CALL timestop(handle)
-   END SUBROUTINE evaluate_core_matrix_p_mix_new
-
 ! **************************************************************************************************
 !> \brief Constructs a new Khon-Sham matrix
 !> \param qs_env ...

src/qs_scf.F

@@ -123,7 +123,7 @@ MODULE qs_scf
                                               set_qs_env
    USE qs_integrate_potential,          ONLY: integrate_v_rspace
    USE qs_kind_types,                   ONLY: qs_kind_type
-   USE qs_ks_methods,                   ONLY: evaluate_core_matrix_p_mix_new,&
+   USE qs_ks_methods,                   ONLY: evaluate_core_matrix_traces,&
                                               qs_ks_update_qs_env
    USE qs_ks_types,                     ONLY: get_ks_env,&
                                               qs_ks_did_change,&
@@ -524,8 +524,9 @@ CONTAINS
             CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.FALSE.)
 
             IF (dft_control%qs_control%xtb_control%do_tblite) THEN
+               CPASSERT(scf_env%mixing_method > 0)
                CALL tb_update_charges(qs_env, dft_control, qs_env%tb_tblite, .FALSE., .FALSE.)
-               CALL evaluate_core_matrix_p_mix_new(qs_env)
+               CALL evaluate_core_matrix_traces(qs_env, rho_ao_ext=scf_env%p_mix_new)
                CALL tb_get_energy(qs_env, qs_env%tb_tblite, energy)
             END IF
 

src/qs_scf_initialization.F

@@ -945,6 +945,16 @@ CONTAINS
                       para_env=para_env, &
                       xas_env=xas_env)
 
+      !Check restricted optimizers available for tblite library
+      IF (dft_control%qs_control%xtb_control%do_tblite) THEN
+         IF (dft_control%lsd) THEN
+            CPABORT("LSD option not compatible with tblite library.")
+         END IF
+         IF (scf_env%method == ot_method_nr) THEN
+            CPABORT("OT SCF option not compatible with tblite library.")
+         END IF
+      END IF
+
       ! Calculate ortho matrix
       ndep = 0
       IF (scf_env%needs_ortho) THEN

src/tblite_interface.F

@@ -74,6 +74,7 @@ MODULE tblite_interface
                                      get_iterator_info, neighbor_list_set_p_type, &
                                      neighbor_list_iterator_p_type, neighbor_list_iterator_release
    USE qs_overlap, ONLY: create_sab_matrix
+   USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
    USE qs_rho_types, ONLY: qs_rho_get, qs_rho_type
    USE qs_scf_types, ONLY: qs_scf_env_type
    USE input_section_types, ONLY: section_vals_get_subs_vals, section_vals_type