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 total energy and forces of excited states
      10             : !> \par History
      11             : !>       01.2020 created
      12             : !> \author JGH
      13             : ! **************************************************************************************************
      14             : MODULE excited_states
      15             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      16             :                                               get_atomic_kind_set
      17             :    USE cp_control_types,                ONLY: dft_control_type
      18             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      19             :                                               cp_logger_get_default_unit_nr,&
      20             :                                               cp_logger_type
      21             :    USE ex_property_calculation,         ONLY: ex_properties
      22             :    USE exstates_types,                  ONLY: excited_energy_type
      23             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      24             :                                               section_vals_type
      25             :    USE qs_energy_types,                 ONLY: qs_energy_type
      26             :    USE qs_environment_types,            ONLY: get_qs_env,&
      27             :                                               qs_environment_type,&
      28             :                                               set_qs_env
      29             :    USE qs_force_types,                  ONLY: allocate_qs_force,&
      30             :                                               deallocate_qs_force,&
      31             :                                               qs_force_type,&
      32             :                                               sum_qs_force,&
      33             :                                               zero_qs_force
      34             :    USE qs_p_env_types,                  ONLY: p_env_release,&
      35             :                                               qs_p_env_type
      36             :    USE response_solver,                 ONLY: response_equation,&
      37             :                                               response_force,&
      38             :                                               response_force_xtb
      39             : #include "./base/base_uses.f90"
      40             : 
      41             :    IMPLICIT NONE
      42             : 
      43             :    PRIVATE
      44             : 
      45             : ! *** Global parameters ***
      46             : 
      47             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'excited_states'
      48             : 
      49             :    PUBLIC :: excited_state_energy
      50             : 
      51             : CONTAINS
      52             : 
      53             : ! **************************************************************************************************
      54             : !> \brief Excited state energy and forces
      55             : !>
      56             : !> \param qs_env ...
      57             : !> \param calculate_forces ...
      58             : !> \par History
      59             : !>       03.2014 created
      60             : !> \author JGH
      61             : ! **************************************************************************************************
      62       30624 :    SUBROUTINE excited_state_energy(qs_env, calculate_forces)
      63             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      64             :       LOGICAL, INTENT(IN), OPTIONAL                      :: calculate_forces
      65             : 
      66             :       CHARACTER(len=*), PARAMETER :: routineN = 'excited_state_energy'
      67             : 
      68             :       INTEGER                                            :: handle, unit_nr
      69       30624 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: natom_of_kind
      70       30624 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      71             :       TYPE(cp_logger_type), POINTER                      :: logger
      72             :       TYPE(dft_control_type), POINTER                    :: dft_control
      73             :       TYPE(excited_energy_type), POINTER                 :: ex_env
      74             :       TYPE(qs_energy_type), POINTER                      :: energy
      75       30624 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: ks_force, lr_force
      76             :       TYPE(qs_p_env_type)                                :: p_env
      77             :       TYPE(section_vals_type), POINTER                   :: tdlr_section
      78             : 
      79       30624 :       CALL timeset(routineN, handle)
      80             : 
      81             :       ! Check for energy correction
      82       30624 :       IF (qs_env%excited_state) THEN
      83        1466 :          logger => cp_get_default_logger()
      84        1466 :          IF (logger%para_env%is_source()) THEN
      85         733 :             unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
      86             :          ELSE
      87         733 :             unit_nr = -1
      88             :          END IF
      89             : 
      90        1466 :          CALL get_qs_env(qs_env, exstate_env=ex_env, energy=energy)
      91             : 
      92        1466 :          energy%excited_state = ex_env%evalue
      93        1466 :          energy%total = energy%total + ex_env%evalue
      94        1466 :          IF (calculate_forces) THEN
      95         550 :             IF (unit_nr > 0) THEN
      96         275 :                WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 27), &
      97         550 :                   " Excited State Forces ", REPEAT("-", 28), "!"
      98             :             END IF
      99             :             ! prepare force array
     100         550 :             CALL get_qs_env(qs_env, force=ks_force, atomic_kind_set=atomic_kind_set)
     101         550 :             CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
     102         550 :             NULLIFY (lr_force)
     103         550 :             CALL allocate_qs_force(lr_force, natom_of_kind)
     104         550 :             DEALLOCATE (natom_of_kind)
     105         550 :             CALL zero_qs_force(lr_force)
     106         550 :             CALL set_qs_env(qs_env, force=lr_force)
     107             :             !
     108         550 :             tdlr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%LINRES")
     109         550 :             CALL response_equation(qs_env, p_env, ex_env%cpmos, unit_nr, tdlr_section)
     110             :             !
     111         550 :             CALL get_qs_env(qs_env, dft_control=dft_control)
     112         550 :             IF (dft_control%qs_control%semi_empirical) THEN
     113           0 :                CPABORT("Not available")
     114         550 :             ELSEIF (dft_control%qs_control%dftb) THEN
     115           0 :                CPABORT("Not available")
     116         550 :             ELSEIF (dft_control%qs_control%xtb) THEN
     117          16 :                CALL response_force_xtb(qs_env, p_env, ex_env%matrix_hz, ex_env, debug=ex_env%debug_forces)
     118             :             ELSE
     119             :                ! KS-DFT
     120             :                CALL response_force(qs_env=qs_env, vh_rspace=ex_env%vh_rspace, &
     121             :                                    vxc_rspace=ex_env%vxc_rspace, vtau_rspace=ex_env%vtau_rspace, &
     122             :                                    vadmm_rspace=ex_env%vadmm_rspace, matrix_hz=ex_env%matrix_hz, &
     123             :                                    matrix_pz=ex_env%matrix_px1, matrix_pz_admm=p_env%p1_admm, &
     124             :                                    matrix_wz=p_env%w1, &
     125             :                                    p_env=p_env, ex_env=ex_env, &
     126         534 :                                    debug=ex_env%debug_forces)
     127             :             END IF
     128             :             ! add TD and KS forces
     129         550 :             CALL get_qs_env(qs_env, force=lr_force)
     130         550 :             CALL sum_qs_force(ks_force, lr_force)
     131         550 :             CALL set_qs_env(qs_env, force=ks_force)
     132         550 :             CALL deallocate_qs_force(lr_force)
     133             :             !
     134         550 :             CALL ex_properties(qs_env, ex_env%matrix_pe, p_env)
     135             :             !
     136         550 :             CALL p_env_release(p_env)
     137             :             !
     138             :          ELSE
     139         916 :             IF (unit_nr > 0) THEN
     140         458 :                WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 27), &
     141         916 :                   " Excited State Energy ", REPEAT("-", 28), "!"
     142         458 :                WRITE (unit_nr, '(T2,A,T75,I6)') "Results for Excited State Nr.", ABS(ex_env%state)
     143         458 :                WRITE (unit_nr, '(T2,A,T65,F16.10)') "Excitation Energy [Hartree] ", ex_env%evalue
     144         458 :                WRITE (unit_nr, '(T2,A,T65,F16.10)') "Total Energy [Hartree]", energy%total
     145             :             END IF
     146             :          END IF
     147             : 
     148        1466 :          IF (unit_nr > 0) THEN
     149         733 :             WRITE (unit_nr, '(T2,A,A,A)') "!", REPEAT("-", 77), "!"
     150             :          END IF
     151             : 
     152             :       END IF
     153             : 
     154       30624 :       CALL timestop(handle)
     155             : 
     156      183744 :    END SUBROUTINE excited_state_energy
     157             : 
     158             : ! **************************************************************************************************
     159             : 
     160             : END MODULE excited_states