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