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 Performs density functional perturbation theory (tddfpt) calculations.
      10             : !>      Uses the self consistent approach. The tddfpt calculation uses the ground
      11             : !>      state of the unperturbed system as the initial state.
      12             : ! **************************************************************************************************
      13             : MODULE qs_tddfpt_module
      14             :    USE bibliography,                    ONLY: Iannuzzi2005,&
      15             :                                               cite_reference
      16             :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type
      17             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      18             :                                               cp_logger_get_default_io_unit,&
      19             :                                               cp_logger_type
      20             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      21             :                                               cp_print_key_unit_nr
      22             :    USE header,                          ONLY: tddfpt_header
      23             :    USE input_constants,                 ONLY: oe_gllb,&
      24             :                                               oe_lb,&
      25             :                                               oe_none,&
      26             :                                               oe_saop,&
      27             :                                               oe_sic,&
      28             :                                               tddfpt_excitations
      29             :    USE input_section_types,             ONLY: section_get_ival,&
      30             :                                               section_vals_create,&
      31             :                                               section_vals_get_subs_vals,&
      32             :                                               section_vals_release,&
      33             :                                               section_vals_retain,&
      34             :                                               section_vals_set_subs_vals,&
      35             :                                               section_vals_type,&
      36             :                                               section_vals_val_get
      37             :    USE qs_energy_types,                 ONLY: qs_energy_type
      38             :    USE qs_environment_types,            ONLY: get_qs_env,&
      39             :                                               qs_environment_type
      40             :    USE qs_ks_methods,                   ONLY: qs_ks_build_kohn_sham_matrix
      41             :    USE qs_ks_types,                     ONLY: qs_ks_env_type,&
      42             :                                               set_ks_env
      43             :    USE qs_p_env_types,                  ONLY: qs_p_env_type
      44             :    USE qs_rho_types,                    ONLY: qs_rho_type
      45             :    USE qs_scf_types,                    ONLY: qs_scf_env_type
      46             :    USE qs_tddfpt_eigensolver,           ONLY: eigensolver
      47             :    USE qs_tddfpt_types,                 ONLY: tddfpt_env_type
      48             :    USE qs_tddfpt_utils,                 ONLY: find_contributions,&
      49             :                                               tddfpt_cleanup,&
      50             :                                               tddfpt_init
      51             :    USE xc_pot_saop,                     ONLY: add_saop_pot
      52             : #include "./base/base_uses.f90"
      53             : 
      54             :    IMPLICIT NONE
      55             : 
      56             :    PRIVATE
      57             :    PUBLIC :: tddfpt_calculation
      58             : 
      59             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt_module'
      60             : 
      61             : CONTAINS
      62             : 
      63             : ! **************************************************************************************************
      64             : !> \brief Performs the perturbation calculation
      65             : !> \param qs_env ...
      66             : ! **************************************************************************************************
      67          12 :    SUBROUTINE tddfpt_calculation(qs_env)
      68             : 
      69             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      70             : 
      71             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_calculation'
      72             : 
      73             :       INTEGER                                            :: handle, iw
      74             :       TYPE(cp_logger_type), POINTER                      :: logger
      75             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
      76             :       TYPE(qs_p_env_type)                                :: p_env
      77             :       TYPE(section_vals_type), POINTER                   :: dft_section, input
      78             :       TYPE(tddfpt_env_type)                              :: t_env
      79             : 
      80          12 :       NULLIFY (logger)
      81          24 :       logger => cp_get_default_logger()
      82          12 :       NULLIFY (input, ks_env)
      83          12 :       CALL get_qs_env(qs_env, ks_env=ks_env, input=input)
      84          12 :       dft_section => section_vals_get_subs_vals(input, "DFT")
      85             : 
      86          12 :       IF (section_get_ival(dft_section, "EXCITATIONS") /= tddfpt_excitations) RETURN
      87          12 :       CALL cite_reference(Iannuzzi2005)
      88             : 
      89          12 :       CALL timeset(routineN, handle)
      90             : 
      91          12 :       IF (section_get_ival(dft_section, "TDDFPT%OE_CORR") /= oe_none) THEN
      92           2 :          CALL orbital_eigenvalue_correction(qs_env)
      93             :       END IF
      94             : 
      95             :       iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PROGRAM_BANNER", &
      96          12 :                                 extension=".Log")
      97          12 :       CALL tddfpt_header(iw)
      98             :       CALL cp_print_key_finished_output(iw, logger, dft_section, &
      99          12 :                                         "PRINT%PROGRAM_BANNER")
     100             : 
     101             :       !---------------------------------------!
     102             :       ! we don't want to update the KS matrix !
     103             :       !---------------------------------------!
     104          12 :       CALL set_ks_env(ks_env, rho_changed=.FALSE.)
     105             : 
     106          12 :       CALL tddfpt_init(p_env, t_env, qs_env)
     107             : 
     108          12 :       CALL eigensolver(p_env, qs_env, t_env)
     109             : 
     110          12 :       CALL find_contributions(qs_env, t_env)
     111             : 
     112          12 :       CALL tddfpt_cleanup(t_env, p_env)
     113             : 
     114          12 :       CALL timestop(handle)
     115             : 
     116          60 :    END SUBROUTINE tddfpt_calculation
     117             : 
     118             : ! **************************************************************************************************
     119             : !> \brief Apply a special potential to obtain better
     120             : !>       orbital eigenvalues.
     121             : !> \param qs_env ...
     122             : ! **************************************************************************************************
     123           4 :    SUBROUTINE orbital_eigenvalue_correction(qs_env)
     124             : 
     125             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     126             : 
     127             :       INTEGER                                            :: oe_corr, output_unit
     128             :       TYPE(cp_logger_type), POINTER                      :: logger
     129           2 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
     130             :       TYPE(qs_energy_type), POINTER                      :: energy
     131             :       TYPE(qs_rho_type), POINTER                         :: rho
     132             :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     133             :       TYPE(section_vals_type), POINTER                   :: input, xc_fun_orig, xc_fun_tmp
     134             : 
     135           0 :       CPASSERT(ASSOCIATED(qs_env))
     136             : 
     137           2 :       NULLIFY (logger, scf_env, input, energy, matrix_ks, rho)
     138           2 :       logger => cp_get_default_logger()
     139           2 :       output_unit = cp_logger_get_default_io_unit(logger)
     140             : 
     141           2 :       IF (output_unit > 0) THEN
     142           1 :          WRITE (output_unit, "(2X,A)") "", &
     143           1 :             "-----------------------------------------------------------------------------", &
     144           1 :             "-                    orbital eigenvalue correction started                  -", &
     145           1 :             "-----------------------------------------------------------------------------", &
     146           2 :             ""
     147             :       END IF
     148             : 
     149             :       CALL get_qs_env(qs_env, &
     150             :                       scf_env=scf_env, &
     151             :                       input=input, &
     152             :                       matrix_ks=matrix_ks, &
     153           2 :                       rho=rho)
     154             : 
     155             :       !----------------------!
     156             :       ! KS matrix without XC !
     157             :       !----------------------!
     158           2 :       xc_fun_orig => section_vals_get_subs_vals(input, "DFT%XC%XC_FUNCTIONAL")
     159           2 :       CALL section_vals_retain(xc_fun_orig)
     160           2 :       NULLIFY (xc_fun_tmp)
     161           2 :       CALL section_vals_create(xc_fun_tmp, xc_fun_orig%section)
     162           2 :       CALL section_vals_set_subs_vals(input, "DFT%XC%XC_FUNCTIONAL", xc_fun_tmp)
     163           2 :       CALL section_vals_release(xc_fun_tmp)
     164             : 
     165           2 :       CALL get_qs_env(qs_env, energy=energy)
     166             :       CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., &
     167           2 :                                         just_energy=.FALSE.)
     168             : 
     169           2 :       CALL section_vals_set_subs_vals(input, "DFT%XC%XC_FUNCTIONAL", xc_fun_orig)
     170           2 :       CALL section_vals_release(xc_fun_orig)
     171             : 
     172           2 :       CALL section_vals_val_get(input, "DFT%TDDFPT%OE_CORR", i_val=oe_corr)
     173             :       IF (oe_corr == oe_saop .OR. &
     174           2 :           oe_corr == oe_lb .OR. &
     175             :           oe_corr == oe_gllb) THEN
     176           2 :          CALL add_saop_pot(matrix_ks, qs_env, oe_corr)
     177             :       ELSE IF (oe_corr == oe_sic) THEN
     178             :       END IF
     179             : 
     180           2 :    END SUBROUTINE orbital_eigenvalue_correction
     181             : 
     182             : END MODULE qs_tddfpt_module