LCOV - code coverage report
Current view: top level - src - qs_local_properties.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.9 % 96 93
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 2

            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 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
        

Generated by: LCOV version 2.0-1