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 31800 : 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 31800 : INTEGER, ALLOCATABLE, DIMENSION(:) :: natom_of_kind
70 31800 : 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 31800 : 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 31800 : CALL timeset(routineN, handle)
80 :
81 : ! Check for energy correction
82 31800 : IF (qs_env%excited_state) THEN
83 1496 : logger => cp_get_default_logger()
84 1496 : IF (logger%para_env%is_source()) THEN
85 748 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
86 : ELSE
87 748 : unit_nr = -1
88 : END IF
89 :
90 1496 : CALL get_qs_env(qs_env, exstate_env=ex_env, energy=energy)
91 :
92 1496 : energy%excited_state = ex_env%evalue
93 1496 : energy%total = energy%total + ex_env%evalue
94 1496 : IF (calculate_forces) THEN
95 564 : IF (unit_nr > 0) THEN
96 282 : WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 27), &
97 564 : " Excited State Forces ", REPEAT("-", 28), "!"
98 : END IF
99 : ! prepare force array
100 564 : CALL get_qs_env(qs_env, force=ks_force, atomic_kind_set=atomic_kind_set)
101 564 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
102 564 : NULLIFY (lr_force)
103 564 : CALL allocate_qs_force(lr_force, natom_of_kind)
104 564 : DEALLOCATE (natom_of_kind)
105 564 : CALL zero_qs_force(lr_force)
106 564 : CALL set_qs_env(qs_env, force=lr_force)
107 : !
108 564 : tdlr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%LINRES")
109 : ! Solve the Z-vector linear equation system
110 564 : CALL response_equation(qs_env, p_env, ex_env%cpmos, unit_nr, tdlr_section)
111 : !
112 564 : CALL get_qs_env(qs_env, dft_control=dft_control)
113 564 : IF (dft_control%qs_control%semi_empirical) THEN
114 0 : CPABORT("Not available")
115 564 : ELSEIF (dft_control%qs_control%dftb) THEN
116 0 : CPABORT("Not available")
117 564 : ELSEIF (dft_control%qs_control%xtb) THEN
118 16 : CALL response_force_xtb(qs_env, p_env, ex_env%matrix_hz, ex_env, debug=ex_env%debug_forces)
119 : ELSE
120 : ! KS-DFT
121 : CALL response_force(qs_env=qs_env, vh_rspace=ex_env%vh_rspace, &
122 : vxc_rspace=ex_env%vxc_rspace, vtau_rspace=ex_env%vtau_rspace, &
123 : vadmm_rspace=ex_env%vadmm_rspace, matrix_hz=ex_env%matrix_hz, &
124 : matrix_pz=p_env%p1, matrix_pz_admm=p_env%p1_admm, &
125 : matrix_wz=p_env%w1, &
126 : p_env=p_env, ex_env=ex_env, &
127 548 : debug=ex_env%debug_forces)
128 : END IF
129 : ! add TD and KS forces
130 564 : CALL get_qs_env(qs_env, force=lr_force)
131 564 : CALL sum_qs_force(ks_force, lr_force)
132 564 : CALL set_qs_env(qs_env, force=ks_force)
133 564 : CALL deallocate_qs_force(lr_force)
134 : !
135 564 : CALL ex_properties(qs_env, ex_env%matrix_pe, p_env)
136 : !
137 564 : CALL p_env_release(p_env)
138 : !
139 : ELSE
140 932 : IF (unit_nr > 0) THEN
141 466 : WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 27), &
142 932 : " Excited State Energy ", REPEAT("-", 28), "!"
143 466 : WRITE (unit_nr, '(T2,A,T75,I6)') "Results for Excited State Nr.", ABS(ex_env%state)
144 466 : WRITE (unit_nr, '(T2,A,T65,F16.10)') "Excitation Energy [Hartree] ", ex_env%evalue
145 466 : WRITE (unit_nr, '(T2,A,T65,F16.10)') "Total Energy [Hartree]", energy%total
146 : END IF
147 : END IF
148 :
149 1496 : IF (unit_nr > 0) THEN
150 748 : WRITE (unit_nr, '(T2,A,A,A)') "!", REPEAT("-", 77), "!"
151 : END IF
152 :
153 : END IF
154 :
155 31800 : CALL timestop(handle)
156 :
157 190800 : END SUBROUTINE excited_state_energy
158 :
159 : ! **************************************************************************************************
160 :
161 : END MODULE excited_states
|
