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 Routines for calculating local energy and stress tensor
      10             : !> \author JGH
      11             : !> \par History
      12             : !>      - 07.2019 created
      13             : ! **************************************************************************************************
      14             : MODULE qs_local_properties
      15             :    USE bibliography,                    ONLY: Cohen2000,&
      16             :                                               Filippetti2000,&
      17             :                                               Rogers2002,&
      18             :                                               cite_reference
      19             :    USE cp_control_types,                ONLY: dft_control_type
      20             :    USE cp_dbcsr_api,                    ONLY: dbcsr_copy,&
      21             :                                               dbcsr_p_type,&
      22             :                                               dbcsr_set,&
      23             :                                               dbcsr_type
      24             :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set
      25             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      26             :                                               cp_logger_get_default_io_unit,&
      27             :                                               cp_logger_type
      28             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      29             :                                               section_vals_type
      30             :    USE kinds,                           ONLY: dp
      31             :    USE mathlib,                         ONLY: det_3x3
      32             :    USE pw_env_types,                    ONLY: pw_env_get,&
      33             :                                               pw_env_type
      34             :    USE pw_methods,                      ONLY: pw_axpy,&
      35             :                                               pw_copy,&
      36             :                                               pw_derive,&
      37             :                                               pw_integrate_function,&
      38             :                                               pw_multiply,&
      39             :                                               pw_transfer,&
      40             :                                               pw_zero
      41             :    USE pw_pool_types,                   ONLY: pw_pool_type
      42             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      43             :                                               pw_r3d_rs_type
      44             :    USE qs_collocate_density,            ONLY: calculate_rho_elec
      45             :    USE qs_core_energies,                ONLY: calculate_ptrace
      46             :    USE qs_energy_types,                 ONLY: qs_energy_type
      47             :    USE qs_environment_types,            ONLY: get_qs_env,&
      48             :                                               qs_environment_type
      49             :    USE qs_ks_methods,                   ONLY: calc_rho_tot_gspace
      50             :    USE qs_ks_types,                     ONLY: qs_ks_env_type,&
      51             :                                               set_ks_env
      52             :    USE qs_matrix_w,                     ONLY: compute_matrix_w
      53             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      54             :                                               qs_rho_type
      55             :    USE qs_vxc,                          ONLY: qs_xc_density
      56             :    USE virial_types,                    ONLY: virial_type
      57             : #include "./base/base_uses.f90"
      58             : 
      59             :    IMPLICIT NONE
      60             : 
      61             :    PRIVATE
      62             : 
      63             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_local_properties'
      64             : 
      65             :    PUBLIC :: qs_local_energy, qs_local_stress
      66             : 
      67             : ! **************************************************************************************************
      68             : 
      69             : CONTAINS
      70             : 
      71             : ! **************************************************************************************************
      72             : !> \brief Routine to calculate the local energy
      73             : !> \param qs_env the qs_env to update
      74             : !> \param energy_density ...
      75             : !> \par History
      76             : !>      07.2019 created
      77             : !> \author JGH
      78             : ! **************************************************************************************************
      79          64 :    SUBROUTINE qs_local_energy(qs_env, energy_density)
      80             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      81             :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: energy_density
      82             : 
      83             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'qs_local_energy'
      84             : 
      85             :       INTEGER                                            :: handle, img, iounit, ispin, nimages, &
      86             :                                                             nkind, nspins
      87             :       LOGICAL                                            :: gapw, gapw_xc
      88             :       REAL(KIND=dp)                                      :: eban, eband, eh, exc, ovol
      89             :       TYPE(cp_logger_type), POINTER                      :: logger
      90          32 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao
      91          32 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks, matrix_s, matrix_w, rho_ao_kp
      92             :       TYPE(dbcsr_type), POINTER                          :: matrix
      93             :       TYPE(dft_control_type), POINTER                    :: dft_control
      94             :       TYPE(pw_c1d_gs_type)                               :: edens_g
      95             :       TYPE(pw_c1d_gs_type), POINTER                      :: rho_core, rho_tot_gspace
      96             :       TYPE(pw_env_type), POINTER                         :: pw_env
      97             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
      98             :       TYPE(pw_r3d_rs_type)                               :: band_density, edens_r, hartree_density, &
      99             :                                                             xc_density
     100          32 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r
     101             :       TYPE(pw_r3d_rs_type), POINTER                      :: rho_tot_rspace, v_hartree_rspace
     102             :       TYPE(qs_energy_type), POINTER                      :: energy
     103             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     104             :       TYPE(qs_rho_type), POINTER                         :: rho, rho_struct
     105             :       TYPE(section_vals_type), POINTER                   :: input, xc_section
     106             : 
     107          32 :       CALL timeset(routineN, handle)
     108             : 
     109          32 :       CALL cite_reference(Cohen2000)
     110             : 
     111          32 :       CPASSERT(ASSOCIATED(qs_env))
     112          32 :       logger => cp_get_default_logger()
     113          32 :       iounit = cp_logger_get_default_io_unit()
     114             : 
     115             :       ! Check for GAPW method : additional terms for local densities
     116          32 :       CALL get_qs_env(qs_env, nkind=nkind, dft_control=dft_control)
     117          32 :       gapw = dft_control%qs_control%gapw
     118          32 :       gapw_xc = dft_control%qs_control%gapw_xc
     119             : 
     120          32 :       nimages = dft_control%nimages
     121          32 :       nspins = dft_control%nspins
     122             : 
     123             :       ! get working arrays
     124          32 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     125          32 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     126          32 :       CALL auxbas_pw_pool%create_pw(band_density)
     127          32 :       CALL auxbas_pw_pool%create_pw(hartree_density)
     128          32 :       CALL auxbas_pw_pool%create_pw(xc_density)
     129             : 
     130             :       ! w matrix
     131          32 :       CALL get_qs_env(qs_env, matrix_w_kp=matrix_w)
     132          32 :       IF (.NOT. ASSOCIATED(matrix_w)) THEN
     133             :          CALL get_qs_env(qs_env, &
     134             :                          ks_env=ks_env, &
     135           2 :                          matrix_s_kp=matrix_s)
     136           2 :          matrix => matrix_s(1, 1)%matrix
     137           2 :          CALL dbcsr_allocate_matrix_set(matrix_w, nspins, nimages)
     138           4 :          DO ispin = 1, nspins
     139           6 :             DO img = 1, nimages
     140           2 :                ALLOCATE (matrix_w(ispin, img)%matrix)
     141           2 :                CALL dbcsr_copy(matrix_w(ispin, img)%matrix, matrix, name="W MATRIX")
     142           4 :                CALL dbcsr_set(matrix_w(ispin, img)%matrix, 0.0_dp)
     143             :             END DO
     144             :          END DO
     145           2 :          CALL set_ks_env(ks_env, matrix_w_kp=matrix_w)
     146             :       END IF
     147             :       ! band structure energy density
     148          32 :       CALL compute_matrix_w(qs_env, .TRUE.)
     149          32 :       CALL get_qs_env(qs_env, ks_env=ks_env, matrix_w_kp=matrix_w)
     150          32 :       CALL auxbas_pw_pool%create_pw(edens_r)
     151          32 :       CALL auxbas_pw_pool%create_pw(edens_g)
     152          32 :       CALL pw_zero(band_density)
     153          64 :       DO ispin = 1, nspins
     154          32 :          rho_ao => matrix_w(ispin, :)
     155             :          CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
     156             :                                  rho=edens_r, &
     157             :                                  rho_gspace=edens_g, &
     158          32 :                                  ks_env=ks_env, soft_valid=(gapw .OR. gapw_xc))
     159          64 :          CALL pw_axpy(edens_r, band_density)
     160             :       END DO
     161          32 :       CALL auxbas_pw_pool%give_back_pw(edens_r)
     162          32 :       CALL auxbas_pw_pool%give_back_pw(edens_g)
     163             : 
     164             :       ! Hartree energy density correction = -0.5 * V_H(r) * [rho(r) - rho_core(r)]
     165          32 :       ALLOCATE (rho_tot_gspace, rho_tot_rspace)
     166          32 :       CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
     167          32 :       CALL auxbas_pw_pool%create_pw(rho_tot_rspace)
     168          32 :       NULLIFY (rho_core)
     169             :       CALL get_qs_env(qs_env, &
     170             :                       v_hartree_rspace=v_hartree_rspace, &
     171          32 :                       rho_core=rho_core, rho=rho)
     172          32 :       CALL qs_rho_get(rho, rho_r=rho_r)
     173          32 :       IF (ASSOCIATED(rho_core)) THEN
     174          32 :          CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
     175          32 :          CALL pw_transfer(rho_core, rho_tot_rspace)
     176             :       ELSE
     177           0 :          CALL pw_zero(rho_tot_rspace)
     178             :       END IF
     179          64 :       DO ispin = 1, nspins
     180          64 :          CALL pw_axpy(rho_r(ispin), rho_tot_rspace, alpha=-1.0_dp)
     181             :       END DO
     182          32 :       CALL pw_zero(hartree_density)
     183          32 :       ovol = 0.5_dp/hartree_density%pw_grid%dvol
     184          32 :       CALL pw_multiply(hartree_density, v_hartree_rspace, rho_tot_rspace, alpha=ovol)
     185          32 :       CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
     186          32 :       CALL auxbas_pw_pool%give_back_pw(rho_tot_rspace)
     187          32 :       DEALLOCATE (rho_tot_gspace, rho_tot_rspace)
     188             : 
     189          32 :       IF (dft_control%do_admm) THEN
     190           0 :          CALL cp_warn(__LOCATION__, "ADMM not supported for local energy calculation")
     191             :       END IF
     192          32 :       IF (gapw_xc .OR. gapw) THEN
     193           0 :          CALL cp_warn(__LOCATION__, "GAPW/GAPW_XC not supported for local energy calculation")
     194             :       END IF
     195             :       ! XC energy density correction = E_xc(r) - V_xc(r)*rho(r)
     196          32 :       CALL get_qs_env(qs_env, input=input)
     197          32 :       xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     198          32 :       CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
     199             :       !
     200          32 :       CALL qs_xc_density(ks_env, rho_struct, xc_section, xc_ener=xc_density)
     201             :       !
     202             :       ! energies
     203          32 :       CALL get_qs_env(qs_env=qs_env, energy=energy)
     204          32 :       eban = pw_integrate_function(band_density)
     205          32 :       eh = pw_integrate_function(hartree_density)
     206          32 :       exc = pw_integrate_function(xc_density)
     207             : 
     208             :       ! band energy
     209          32 :       CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks)
     210          32 :       CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
     211          32 :       CALL calculate_ptrace(matrix_ks, rho_ao_kp, eband, nspins)
     212             : 
     213             :       ! get full density
     214          32 :       CALL pw_copy(band_density, energy_density)
     215          32 :       CALL pw_axpy(hartree_density, energy_density)
     216          32 :       CALL pw_axpy(xc_density, energy_density)
     217             : 
     218          32 :       IF (iounit > 0) THEN
     219          16 :          WRITE (UNIT=iounit, FMT="(/,T3,A)") REPEAT("=", 78)
     220          16 :          WRITE (UNIT=iounit, FMT="(T4,A,T52,A,T75,A)") "Local Energy Calculation", "GPW/GAPW", "Local"
     221          16 :          WRITE (UNIT=iounit, FMT="(T4,A,T45,F15.8,T65,F15.8)") "Band Energy", eband, eban
     222          16 :          WRITE (UNIT=iounit, FMT="(T4,A,T65,F15.8)") "Hartree Energy Correction", eh
     223          16 :          WRITE (UNIT=iounit, FMT="(T4,A,T65,F15.8)") "XC Energy Correction", exc
     224          16 :          WRITE (UNIT=iounit, FMT="(T4,A,T45,F15.8,T65,F15.8)") "Total Energy", &
     225          32 :             energy%total, eban + eh + exc + energy%core_overlap + energy%core_self + energy%dispersion
     226          16 :          WRITE (UNIT=iounit, FMT="(T3,A)") REPEAT("=", 78)
     227             :       END IF
     228             : 
     229             :       ! return temp arrays
     230          32 :       CALL auxbas_pw_pool%give_back_pw(band_density)
     231          32 :       CALL auxbas_pw_pool%give_back_pw(hartree_density)
     232          32 :       CALL auxbas_pw_pool%give_back_pw(xc_density)
     233             : 
     234          32 :       CALL timestop(handle)
     235             : 
     236          32 :    END SUBROUTINE qs_local_energy
     237             : 
     238             : ! **************************************************************************************************
     239             : !> \brief Routine to calculate the local stress
     240             : !> \param qs_env the qs_env to update
     241             : !> \param stress_tensor ...
     242             : !> \param beta ...
     243             : !> \par History
     244             : !>      07.2019 created
     245             : !> \author JGH
     246             : ! **************************************************************************************************
     247          30 :    SUBROUTINE qs_local_stress(qs_env, stress_tensor, beta)
     248             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     249             :       TYPE(pw_r3d_rs_type), DIMENSION(:, :), &
     250             :          INTENT(INOUT), OPTIONAL                         :: stress_tensor
     251             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: beta
     252             : 
     253             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'qs_local_stress'
     254             : 
     255             :       INTEGER                                            :: handle, i, iounit, j, nimages, nkind, &
     256             :                                                             nspins
     257             :       LOGICAL                                            :: do_stress, gapw, gapw_xc, use_virial
     258             :       REAL(KIND=dp)                                      :: my_beta
     259             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_loc
     260             :       TYPE(cp_logger_type), POINTER                      :: logger
     261             :       TYPE(dft_control_type), POINTER                    :: dft_control
     262             :       TYPE(pw_c1d_gs_type)                               :: v_hartree_gspace
     263         120 :       TYPE(pw_c1d_gs_type), DIMENSION(3)                 :: efield
     264             :       TYPE(pw_env_type), POINTER                         :: pw_env
     265             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     266             :       TYPE(pw_r3d_rs_type)                               :: xc_density
     267             :       TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_rspace
     268             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     269             :       TYPE(qs_rho_type), POINTER                         :: rho_struct
     270             :       TYPE(section_vals_type), POINTER                   :: input, xc_section
     271             :       TYPE(virial_type), POINTER                         :: virial
     272             : 
     273          30 :       CALL cp_warn(__LOCATION__, "Local Stress Tensor code is not working, skipping")
     274          30 :       RETURN
     275             : 
     276             :       CALL timeset(routineN, handle)
     277             : 
     278             :       CALL cite_reference(Filippetti2000)
     279             :       CALL cite_reference(Rogers2002)
     280             : 
     281             :       CPASSERT(ASSOCIATED(qs_env))
     282             : 
     283             :       IF (PRESENT(stress_tensor)) THEN
     284             :          do_stress = .TRUE.
     285             :       ELSE
     286             :          do_stress = .FALSE.
     287             :       END IF
     288             :       IF (PRESENT(beta)) THEN
     289             :          my_beta = beta
     290             :       ELSE
     291             :          my_beta = 0.0_dp
     292             :       END IF
     293             : 
     294             :       logger => cp_get_default_logger()
     295             :       iounit = cp_logger_get_default_io_unit()
     296             : 
     297             :       !!!!!!
     298             :       CALL cp_warn(__LOCATION__, "Local Stress Tensor code is not tested")
     299             :       !!!!!!
     300             : 
     301             :       ! Check for GAPW method : additional terms for local densities
     302             :       CALL get_qs_env(qs_env, nkind=nkind, dft_control=dft_control)
     303             :       gapw = dft_control%qs_control%gapw
     304             :       gapw_xc = dft_control%qs_control%gapw_xc
     305             : 
     306             :       nimages = dft_control%nimages
     307             :       nspins = dft_control%nspins
     308             : 
     309             :       ! get working arrays
     310             :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     311             :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     312             :       CALL auxbas_pw_pool%create_pw(xc_density)
     313             : 
     314             :       ! init local stress tensor
     315             :       IF (do_stress) THEN
     316             :          DO i = 1, 3
     317             :             DO j = 1, 3
     318             :                CALL pw_zero(stress_tensor(i, j))
     319             :             END DO
     320             :          END DO
     321             :       END IF
     322             : 
     323             :       IF (dft_control%do_admm) THEN
     324             :          CALL cp_warn(__LOCATION__, "ADMM not supported for local energy calculation")
     325             :       END IF
     326             :       IF (gapw_xc .OR. gapw) THEN
     327             :          CALL cp_warn(__LOCATION__, "GAPW/GAPW_XC not supported for local energy calculation")
     328             :       END IF
     329             :       ! XC energy density correction = E_xc(r) - V_xc(r)*rho(r)
     330             :       CALL get_qs_env(qs_env, ks_env=ks_env, input=input, rho=rho_struct)
     331             :       xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     332             :       !
     333             :       CALL qs_xc_density(ks_env, rho_struct, xc_section, xc_ener=xc_density)
     334             : 
     335             :       ! Electrical field terms
     336             :       CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace)
     337             :       CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
     338             :       CALL pw_transfer(v_hartree_rspace, v_hartree_gspace)
     339             :       DO i = 1, 3
     340             :          CALL auxbas_pw_pool%create_pw(efield(i))
     341             :          CALL pw_copy(v_hartree_gspace, efield(i))
     342             :       END DO
     343             :       CALL pw_derive(efield(1), (/1, 0, 0/))
     344             :       CALL pw_derive(efield(2), (/0, 1, 0/))
     345             :       CALL pw_derive(efield(3), (/0, 0, 1/))
     346             :       CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
     347             :       DO i = 1, 3
     348             :          CALL auxbas_pw_pool%give_back_pw(efield(i))
     349             :       END DO
     350             : 
     351             :       pv_loc = 0.0_dp
     352             : 
     353             :       CALL get_qs_env(qs_env=qs_env, virial=virial)
     354             :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     355             :       IF (.NOT. use_virial) THEN
     356             :          CALL cp_warn(__LOCATION__, "Local stress should be used with standard stress calculation.")
     357             :       END IF
     358             :       IF (iounit > 0 .AND. use_virial) THEN
     359             :          WRITE (UNIT=iounit, FMT="(/,T3,A)") REPEAT("=", 78)
     360             :          WRITE (UNIT=iounit, FMT="(T4,A)") "Local Stress Calculation"
     361             :          WRITE (UNIT=iounit, FMT="(T42,A,T64,A)") "       1/3 Trace", "     Determinant"
     362             :          WRITE (UNIT=iounit, FMT="(T4,A,T42,F16.8,T64,F16.8)") "Total Stress", &
     363             :             (pv_loc(1, 1) + pv_loc(2, 2) + pv_loc(3, 3))/3.0_dp, det_3x3(pv_loc)
     364             :          WRITE (UNIT=iounit, FMT="(T3,A)") REPEAT("=", 78)
     365             :       END IF
     366             : 
     367             :       CALL timestop(handle)
     368             : 
     369          30 :    END SUBROUTINE qs_local_stress
     370             : 
     371             : ! **************************************************************************************************
     372             : 
     373             : END MODULE qs_local_properties