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             : ! **************************************************************************************************
       9             : !> \brief tblite matrix build
      10             : !> \author JVP
      11             : !> \history creation 09.2024
      12             : ! **************************************************************************************************
      13             : 
      14             : MODULE tblite_ks_matrix
      15             : 
      16             :    USE cp_control_types,                ONLY: dft_control_type
      17             :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      18             :                                               dbcsr_copy,&
      19             :                                               dbcsr_multiply,&
      20             :                                               dbcsr_p_type,&
      21             :                                               dbcsr_type
      22             :    USE cp_dbcsr_contrib,                ONLY: dbcsr_dot
      23             :    USE kinds,                           ONLY: dp
      24             :    USE message_passing,                 ONLY: mp_para_env_type
      25             :    USE qs_energy_types,                 ONLY: qs_energy_type
      26             :    USE qs_environment_types,            ONLY: get_qs_env,&
      27             :                                               qs_environment_type
      28             :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      29             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      30             :                                               mo_set_type
      31             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      32             :                                               qs_rho_type
      33             :    USE tblite_interface,                ONLY: tb_derive_dH_diag,&
      34             :                                               tb_derive_dH_off,&
      35             :                                               tb_ham_add_coulomb,&
      36             :                                               tb_update_charges
      37             : #include "./base/base_uses.f90"
      38             : 
      39             :    IMPLICIT NONE
      40             : 
      41             :    PRIVATE
      42             : 
      43             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tblite_ks_matrix'
      44             : 
      45             :    PUBLIC :: build_tblite_ks_matrix
      46             : 
      47             : CONTAINS
      48             : 
      49             : ! **************************************************************************************************
      50             : !> \brief ...
      51             : !> \param qs_env ...
      52             : !> \param calculate_forces ...
      53             : !> \param just_energy ...
      54             : !> \param ext_ks_matrix ...
      55             : ! **************************************************************************************************
      56        9102 :    SUBROUTINE build_tblite_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
      57             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      58             :       LOGICAL, INTENT(in)                                :: calculate_forces, just_energy
      59             :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
      60             :          POINTER                                         :: ext_ks_matrix
      61             : 
      62             :       CHARACTER(len=*), PARAMETER :: routineN = 'build_tblite_ks_matrix'
      63             : 
      64             :       INTEGER                                            :: handle, img, ispin, nimg, ns, nspins
      65             :       LOGICAL                                            :: do_efield
      66             :       REAL(KIND=dp)                                      :: pc_ener, qmmm_el
      67        9102 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p1, mo_derivs
      68        9102 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: ks_matrix, matrix_h
      69             :       TYPE(dbcsr_type), POINTER                          :: mo_coeff
      70             :       TYPE(dft_control_type), POINTER                    :: dft_control
      71             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      72             :       TYPE(qs_energy_type), POINTER                      :: energy
      73             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
      74             :       TYPE(qs_rho_type), POINTER                         :: rho
      75             : 
      76        9102 :       CALL timeset(routineN, handle)
      77             : 
      78        9102 :       NULLIFY (dft_control, ks_env, ks_matrix, rho, energy)
      79        9102 :       CPASSERT(ASSOCIATED(qs_env))
      80             : 
      81             :       CALL get_qs_env(qs_env, &
      82             :                       dft_control=dft_control, &
      83             :                       matrix_h_kp=matrix_h, &
      84             :                       para_env=para_env, &
      85             :                       ks_env=ks_env, &
      86             :                       matrix_ks_kp=ks_matrix, &
      87             :                       rho=rho, &
      88        9102 :                       energy=energy)
      89             : 
      90        9102 :       IF (PRESENT(ext_ks_matrix)) THEN
      91             :          ! remap pointer to allow for non-kpoint external ks matrix
      92             :          ! ext_ks_matrix is used in linear response code
      93           0 :          ns = SIZE(ext_ks_matrix)
      94           0 :          ks_matrix(1:ns, 1:1) => ext_ks_matrix(1:ns)
      95             :       END IF
      96             : 
      97        9102 :       energy%qmmm_el = 0.0_dp
      98             : 
      99        9102 :       nspins = dft_control%nspins
     100        9102 :       nimg = dft_control%nimages
     101        9102 :       CPASSERT(ASSOCIATED(matrix_h))
     102        9102 :       CPASSERT(ASSOCIATED(rho))
     103       27306 :       CPASSERT(SIZE(ks_matrix) > 0)
     104             : 
     105       18204 :       DO ispin = 1, nspins
     106       27306 :          DO img = 1, nimg
     107             :             ! copy the core matrix into the fock matrix
     108       18204 :             CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
     109             :          END DO
     110             :       END DO
     111             : 
     112        9102 :       IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
     113             :           dft_control%apply_efield_field) THEN
     114           0 :          do_efield = .TRUE.
     115           0 :          CPABORT("Not implemented yet. Use CP2K routines for GFN1")
     116             :       ELSE
     117        9102 :          do_efield = .FALSE.
     118             :       END IF
     119             : 
     120        9102 :       CALL tb_update_charges(qs_env, dft_control, qs_env%tb_tblite, calculate_forces, .TRUE.)
     121             : 
     122        9102 :       CALL tb_ham_add_coulomb(qs_env, qs_env%tb_tblite, dft_control)
     123             : 
     124        9102 :       IF (qs_env%qmmm) THEN
     125           0 :          CPASSERT(SIZE(ks_matrix, 2) == 1)
     126           0 :          DO ispin = 1, nspins
     127             :             ! If QM/MM sumup the 1el Hamiltonian
     128             :             CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
     129           0 :                            1.0_dp, 1.0_dp)
     130           0 :             CALL qs_rho_get(rho, rho_ao=matrix_p1)
     131             :             ! Compute QM/MM Energy
     132             :             CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
     133           0 :                            matrix_p1(ispin)%matrix, qmmm_el)
     134           0 :             energy%qmmm_el = energy%qmmm_el + qmmm_el
     135             :          END DO
     136           0 :          pc_ener = qs_env%ks_qmmm_env%pc_ener
     137           0 :          energy%qmmm_el = energy%qmmm_el + pc_ener
     138             :       END IF
     139             : 
     140        9102 :       IF (calculate_forces) THEN
     141          28 :          CALL tb_derive_dH_diag(qs_env, .TRUE.)
     142          28 :          CALL tb_derive_dH_off(qs_env, .TRUE.)
     143             :       END IF
     144             : 
     145             :       ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
     146        9102 :       IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
     147           0 :          CPASSERT(SIZE(ks_matrix, 2) == 1)
     148             :          BLOCK
     149           0 :             TYPE(mo_set_type), DIMENSION(:), POINTER         :: mo_array
     150           0 :             CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
     151           0 :             DO ispin = 1, SIZE(mo_derivs)
     152           0 :                CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
     153           0 :                IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
     154           0 :                   CPABORT("")
     155             :                END IF
     156             :                CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
     157           0 :                                    0.0_dp, mo_derivs(ispin)%matrix)
     158             :             END DO
     159             :          END BLOCK
     160             :       END IF
     161             : 
     162        9102 :       CALL timestop(handle)
     163             : 
     164        9102 :    END SUBROUTINE build_tblite_ks_matrix
     165             : 
     166             : END MODULE tblite_ks_matrix