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 Utility subroutine for qs energy calculation
      10             : !> \par History
      11             : !>      none
      12             : !> \author MK (29.10.2002)
      13             : ! **************************************************************************************************
      14             : MODULE qs_matrix_w
      15             :    USE cp_control_types,                ONLY: dft_control_type
      16             :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type,&
      17             :                                               dbcsr_set
      18             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      19             :                                               cp_fm_struct_release,&
      20             :                                               cp_fm_struct_type
      21             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      22             :                                               cp_fm_release,&
      23             :                                               cp_fm_type
      24             :    USE kinds,                           ONLY: dp
      25             :    USE kpoint_methods,                  ONLY: kpoint_density_matrices,&
      26             :                                               kpoint_density_transform
      27             :    USE kpoint_types,                    ONLY: kpoint_type
      28             :    USE qs_density_matrices,             ONLY: calculate_w_matrix,&
      29             :                                               calculate_w_matrix_ot
      30             :    USE qs_environment_types,            ONLY: get_qs_env,&
      31             :                                               qs_environment_type
      32             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      33             :                                               mo_set_type
      34             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      35             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      36             :                                               qs_rho_type
      37             :    USE scf_control_types,               ONLY: scf_control_type
      38             : #include "./base/base_uses.f90"
      39             : 
      40             :    IMPLICIT NONE
      41             : 
      42             :    PRIVATE
      43             : 
      44             : ! *** Global parameters ***
      45             : 
      46             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_matrix_w'
      47             : 
      48             :    PUBLIC :: compute_matrix_w
      49             : 
      50             : CONTAINS
      51             : 
      52             : ! **************************************************************************************************
      53             : !> \brief Refactoring of qs_energies_scf. Moves computation of matrix_w
      54             : !>        into separate subroutine
      55             : !> \param qs_env ...
      56             : !> \param calc_forces ...
      57             : !> \par History
      58             : !>      05.2013 created [Florian Schiffmann]
      59             : ! **************************************************************************************************
      60             : 
      61       18935 :    SUBROUTINE compute_matrix_w(qs_env, calc_forces)
      62             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      63             :       LOGICAL, INTENT(IN)                                :: calc_forces
      64             : 
      65             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'compute_matrix_w'
      66             : 
      67             :       INTEGER                                            :: handle, is, ispin, nao, nspin
      68             :       LOGICAL                                            :: do_kpoints, has_unit_metric
      69       18935 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_s, matrix_w, &
      70       18935 :                                                             mo_derivs, rho_ao
      71             :       TYPE(dft_control_type), POINTER                    :: dft_control
      72       18935 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
      73             :       TYPE(mo_set_type), POINTER                         :: mo_set
      74             :       TYPE(qs_rho_type), POINTER                         :: rho
      75             :       TYPE(scf_control_type), POINTER                    :: scf_control
      76             : 
      77       18935 :       CALL timeset(routineN, handle)
      78             : 
      79             :       ! if calculate forces, time to compute the w matrix
      80       18935 :       CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
      81             : 
      82       18935 :       IF (calc_forces .AND. .NOT. has_unit_metric) THEN
      83        5593 :          CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
      84             : 
      85        5593 :          IF (do_kpoints) THEN
      86         150 :             BLOCK
      87         450 :                TYPE(cp_fm_type), DIMENSION(2)                   :: fmwork
      88             :                TYPE(cp_fm_struct_type), POINTER                   :: ao_ao_fmstruct
      89             :                TYPE(cp_fm_type), POINTER                          :: mo_coeff
      90         150 :                TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s_kp, matrix_w_kp
      91             :                TYPE(kpoint_type), POINTER                         :: kpoints
      92             :                TYPE(neighbor_list_set_p_type), DIMENSION(:), &
      93         150 :                   POINTER                                         :: sab_nl
      94             : 
      95             :                CALL get_qs_env(qs_env, &
      96             :                                matrix_w_kp=matrix_w_kp, &
      97             :                                matrix_s_kp=matrix_s_kp, &
      98             :                                sab_orb=sab_nl, &
      99             :                                mos=mos, &
     100         150 :                                kpoints=kpoints)
     101             : 
     102         150 :                CALL get_mo_set(mos(1), mo_coeff=mo_coeff, nao=nao)
     103             :                CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, nrow_global=nao, ncol_global=nao, &
     104         150 :                                         template_fmstruct=mo_coeff%matrix_struct)
     105             : 
     106         450 :                DO is = 1, SIZE(fmwork)
     107         450 :                   CALL cp_fm_create(fmwork(is), matrix_struct=ao_ao_fmstruct)
     108             :                END DO
     109         150 :                CALL cp_fm_struct_release(ao_ao_fmstruct)
     110             : 
     111             :                ! energy weighted density matrices in k-space
     112         150 :                CALL kpoint_density_matrices(kpoints, energy_weighted=.TRUE.)
     113             :                ! energy weighted density matrices in real space
     114             :                CALL kpoint_density_transform(kpoints, matrix_w_kp, .TRUE., &
     115         150 :                                              matrix_s_kp(1, 1)%matrix, sab_nl, fmwork)
     116             : 
     117         600 :                DO is = 1, SIZE(fmwork)
     118         450 :                   CALL cp_fm_release(fmwork(is))
     119             :                END DO
     120             : 
     121             :             END BLOCK
     122             :          ELSE
     123             : 
     124        5443 :             NULLIFY (dft_control, rho_ao)
     125             :             CALL get_qs_env(qs_env, &
     126             :                             matrix_w=matrix_w, &
     127             :                             matrix_ks=matrix_ks, &
     128             :                             matrix_s=matrix_s, &
     129             :                             mo_derivs=mo_derivs, &
     130             :                             scf_control=scf_control, &
     131             :                             mos=mos, &
     132             :                             rho=rho, &
     133        5443 :                             dft_control=dft_control)
     134             : 
     135        5443 :             CALL qs_rho_get(rho, rho_ao=rho_ao)
     136             : 
     137        5443 :             nspin = SIZE(mos)
     138       11522 :             DO ispin = 1, nspin
     139        6079 :                mo_set => mos(ispin)
     140       11522 :                IF (dft_control%roks) THEN
     141         168 :                   IF (scf_control%use_ot) THEN
     142         116 :                      IF (ispin > 1) THEN
     143             :                         ! not very elegant, indeed ...
     144          58 :                         CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
     145             :                      ELSE
     146             :                         CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, &
     147          58 :                                                    matrix_w(ispin)%matrix, matrix_s(1)%matrix)
     148             :                      END IF
     149             :                   ELSE
     150             :                      CALL calculate_w_matrix(mo_set=mo_set, &
     151             :                                              matrix_ks=matrix_ks(ispin)%matrix, &
     152             :                                              matrix_p=rho_ao(ispin)%matrix, &
     153          52 :                                              matrix_w=matrix_w(ispin)%matrix)
     154             :                   END IF
     155             :                ELSE
     156        5911 :                   IF (scf_control%use_ot) THEN
     157             :                      CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, &
     158        2413 :                                                 matrix_w(ispin)%matrix, matrix_s(1)%matrix)
     159             :                   ELSE
     160        3498 :                      CALL calculate_w_matrix(mo_set, matrix_w(ispin)%matrix)
     161             :                   END IF
     162             :                END IF
     163             :             END DO
     164             : 
     165             :          END IF
     166             : 
     167             :       END IF
     168             : 
     169       18935 :       CALL timestop(handle)
     170             : 
     171       18935 :    END SUBROUTINE compute_matrix_w
     172             : 
     173             : END MODULE qs_matrix_w