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 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 32282 : 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 32282 : INTEGER, ALLOCATABLE, DIMENSION(:) :: natom_of_kind
70 32282 : 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 32282 : 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 32282 : CALL timeset(routineN, handle)
80 :
81 : ! Check for energy correction
82 32282 : IF (qs_env%excited_state) THEN
83 1470 : logger => cp_get_default_logger()
84 1470 : IF (logger%para_env%is_source()) THEN
85 735 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
86 : ELSE
87 735 : unit_nr = -1
88 : END IF
89 :
90 1470 : CALL get_qs_env(qs_env, exstate_env=ex_env, energy=energy)
91 :
92 1470 : energy%excited_state = ex_env%evalue
93 1470 : energy%total = energy%total + ex_env%evalue
94 1470 : IF (calculate_forces) THEN
95 552 : IF (unit_nr > 0) THEN
96 276 : WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 27), &
97 552 : " Excited State Forces ", REPEAT("-", 28), "!"
98 : END IF
99 : ! prepare force array
100 552 : CALL get_qs_env(qs_env, force=ks_force, atomic_kind_set=atomic_kind_set)
101 552 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
102 552 : NULLIFY (lr_force)
103 552 : CALL allocate_qs_force(lr_force, natom_of_kind)
104 552 : DEALLOCATE (natom_of_kind)
105 552 : CALL zero_qs_force(lr_force)
106 552 : CALL set_qs_env(qs_env, force=lr_force)
107 : !
108 552 : tdlr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%LINRES")
109 552 : CALL response_equation(qs_env, p_env, ex_env%cpmos, unit_nr, tdlr_section)
110 : !
111 552 : CALL get_qs_env(qs_env, dft_control=dft_control)
112 552 : IF (dft_control%qs_control%semi_empirical) THEN
113 0 : CPABORT("Not available")
114 552 : ELSEIF (dft_control%qs_control%dftb) THEN
115 0 : CPABORT("Not available")
116 552 : 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 536 : debug=ex_env%debug_forces)
127 : END IF
128 : ! add TD and KS forces
129 552 : CALL get_qs_env(qs_env, force=lr_force)
130 552 : CALL sum_qs_force(ks_force, lr_force)
131 552 : CALL set_qs_env(qs_env, force=ks_force)
132 552 : CALL deallocate_qs_force(lr_force)
133 : !
134 552 : CALL ex_properties(qs_env, ex_env%matrix_pe, p_env)
135 : !
136 552 : CALL p_env_release(p_env)
137 : !
138 : ELSE
139 918 : IF (unit_nr > 0) THEN
140 459 : WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 27), &
141 918 : " Excited State Energy ", REPEAT("-", 28), "!"
142 459 : WRITE (unit_nr, '(T2,A,T75,I6)') "Results for Excited State Nr.", ABS(ex_env%state)
143 459 : WRITE (unit_nr, '(T2,A,T65,F16.10)') "Excitation Energy [Hartree] ", ex_env%evalue
144 459 : WRITE (unit_nr, '(T2,A,T65,F16.10)') "Total Energy [Hartree]", energy%total
145 : END IF
146 : END IF
147 :
148 1470 : IF (unit_nr > 0) THEN
149 735 : WRITE (unit_nr, '(T2,A,A,A)') "!", REPEAT("-", 77), "!"
150 : END IF
151 :
152 : END IF
153 :
154 32282 : CALL timestop(handle)
155 :
156 193692 : END SUBROUTINE excited_state_energy
157 :
158 : ! **************************************************************************************************
159 :
160 : END MODULE excited_states
|
