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 : !> \brief Simplified Tamm Dancoff approach (sTDA).
9 : ! **************************************************************************************************
10 : MODULE qs_tddfpt2_stda_types
11 :
12 : USE atomic_kind_types, ONLY: atomic_kind_type,&
13 : get_atomic_kind
14 : USE cp_control_types, ONLY: dft_control_type,&
15 : stda_control_type
16 : USE cp_log_handling, ONLY: cp_get_default_logger,&
17 : cp_logger_type
18 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
19 : cp_print_key_unit_nr
20 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
21 : section_vals_type
22 : USE kinds, ONLY: dp
23 : USE physcon, ONLY: evolt
24 : USE qs_environment_types, ONLY: get_qs_env,&
25 : qs_environment_type
26 : USE qs_kind_types, ONLY: qs_kind_type
27 : #include "./base/base_uses.f90"
28 :
29 : IMPLICIT NONE
30 :
31 : PRIVATE
32 :
33 : TYPE stda_kind_type
34 : CHARACTER(LEN=2) :: symbol = "" !element symbol
35 : INTEGER :: z = -1 !atomic_number
36 : INTEGER :: kind_number = -1 !kind number
37 : REAL(KIND=dp) :: hardness_param = -1.0_dp !hardness_parameter eta
38 : REAL(KIND=dp) :: rcut = -1.0_dp !cutoff radius for short-range Coulomb
39 : END TYPE
40 :
41 : TYPE stda_kind_p_type
42 : TYPE(stda_kind_type), POINTER :: kind_param => NULL()
43 : END TYPE
44 :
45 : TYPE stda_env_type
46 : !fraction of non-local HF exchange
47 : REAL(KIND=dp) :: hfx_fraction = -1.0_dp
48 : LOGICAL :: do_exchange = .FALSE.
49 : ! empirical parameters
50 : REAL(KIND=dp) :: alpha_param = -1.0_dp
51 : REAL(KIND=dp) :: beta_param = -1.0_dp
52 : ! filter for TD matrix
53 : REAL(KIND=dp) :: eps_td_filter = -1.0_dp
54 : TYPE(stda_kind_p_type), DIMENSION(:), POINTER:: kind_param_set => NULL()
55 : !number of atomic orbitals, number of occupied orbitals
56 : INTEGER, DIMENSION(2) :: n_ao = -1
57 : INTEGER, DIMENSION(2) :: nactive = -1
58 : END TYPE
59 :
60 : !PARAMETERS
61 : !&<
62 : INTEGER, PARAMETER, PRIVATE :: nelem = 103
63 : ! H He
64 : ! Li Be B C N O F Ne
65 : ! Na Mg Al Si P S Cl Ar
66 : ! K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr
67 : ! Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe
68 : ! Cs Ba La Ce-Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn
69 : ! Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr 103
70 :
71 : REAL(KIND=dp), DIMENSION(1:nelem), PARAMETER, PRIVATE:: hardness = &
72 : (/6.4299544220_dp, 12.544911890_dp, & ! 2 H-He
73 : 2.3745866560_dp, 3.4967633530_dp, 4.6190089720_dp, 5.7409789220_dp, &
74 : 6.8624665290_dp, 7.9854357010_dp, 9.1064753720_dp, 10.23034050_dp, & ! 8 Li-Ne
75 : 2.444141360_dp, 3.0146513830_dp, 3.5849070740_dp, 4.15513090_dp, &
76 : 4.7258039740_dp, 5.2959792410_dp, 5.8661864840_dp, 6.4366187140_dp, & ! 8 Na-Ar
77 : 2.3273178360_dp, 2.7587238140_dp, 2.8581921140_dp, 2.9578300430_dp, &
78 : 3.0573410060_dp, 3.1567254290_dp, 3.2563827230_dp, 3.3559314050_dp, &
79 : 3.4556091170_dp, 3.5550133130_dp, 3.6544183480_dp, 3.7541601450_dp, &
80 : 4.1855197930_dp, 4.6166272460_dp, 5.0662145070_dp, 5.4794960970_dp, &
81 : 5.9110996450_dp, 6.3418467680_dp, & ! 18 K-Kr
82 : 2.1204582570_dp, 2.5373700480_dp, 2.6335468980_dp, 2.7297528930_dp, &
83 : 2.8259738860_dp, 2.9221296040_dp, 3.0183708780_dp, 3.1145981770_dp, &
84 : 3.210756280_dp, 3.3069474480_dp, 3.4031948570_dp, 3.4993761390_dp, &
85 : 3.9163692460_dp, 4.3332332190_dp, 4.7500787860_dp, 5.1669793270_dp, &
86 : 5.5838871020_dp, 6.000897330_dp, & ! 18 Rb-Xe
87 : 0.6829150240_dp, 0.9200946840_dp, 1.1570887860_dp, 1.39427570_dp, &
88 : 1.6314731730_dp, 1.8684389980_dp, 2.1056577930_dp, 2.3426646420_dp, &
89 : 2.5798149820_dp, 2.8170264230_dp, 3.0540365330_dp, 3.2911692310_dp, &
90 : 3.5282971610_dp, 3.7655249290_dp, 4.0025547030_dp, 4.2394783410_dp, &
91 : 4.4765830210_dp, 4.7065224490_dp, 4.9508466940_dp, 5.1879311720_dp, &
92 : 5.4256076210_dp, 5.6619144310_dp, 5.900042920_dp, 6.1367145320_dp, &
93 : 6.3741299770_dp, 6.6102656130_dp, 1.7043485810_dp, 1.9413526120_dp, &
94 : 2.178491510_dp, 2.4158121060_dp, 2.6527780840_dp, 2.8899554570_dp, & ! 32 Cs-Rn
95 : 0.9882529880_dp, 1.2819499970_dp, 1.3497250380_dp, 1.4175257380_dp, &
96 : 1.9368567520_dp, 2.2305576050_dp, 2.5241204960_dp, 3.0436128480_dp, &
97 : 3.4168675260_dp, 3.4049844440_dp, 3.9244199680_dp, 4.2180813280_dp, &
98 : 4.5115926320_dp, 4.8050928950_dp, 5.0989816210_dp, 5.3926054620_dp, &
99 : 5.4606987930_dp/) ! 17 Fr-Lr
100 : !&>
101 :
102 : REAL(KIND=dp), DIMENSION(2), PARAMETER, PRIVATE:: alpha = (/1.420_dp, 0.480_dp/)
103 : REAL(KIND=dp), DIMENSION(2), PARAMETER, PRIVATE:: beta = (/0.200_dp, 1.830_dp/)
104 :
105 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_stda_types'
106 :
107 : PUBLIC:: stda_env_type, &
108 : allocate_stda_env, deallocate_stda_env, stda_init_param
109 :
110 : CONTAINS
111 :
112 : ! **************************************************************************************************
113 : !> \brief Get the parameters needed for an sTDA calculation
114 : !> \param qs_env ...
115 : !> \param stda_kernel ...
116 : !> \param stda_control ...
117 : ! **************************************************************************************************
118 402 : SUBROUTINE stda_init_param(qs_env, stda_kernel, stda_control)
119 :
120 : TYPE(qs_environment_type), POINTER :: qs_env
121 : TYPE(stda_env_type) :: stda_kernel
122 : TYPE(stda_control_type) :: stda_control
123 :
124 : INTEGER :: ikind, log_unit, nkind
125 : REAL(KIND=dp) :: eta, fxx, rcut
126 402 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
127 : TYPE(atomic_kind_type), POINTER :: atomic_kind
128 : TYPE(cp_logger_type), POINTER :: logger
129 : TYPE(dft_control_type), POINTER :: dft_control
130 : TYPE(section_vals_type), POINTER :: tddfpt_print_section
131 : TYPE(stda_kind_type), POINTER :: kind_param
132 :
133 402 : NULLIFY (logger)
134 804 : logger => cp_get_default_logger()
135 :
136 402 : CPASSERT(ASSOCIATED(stda_kernel%kind_param_set))
137 :
138 402 : NULLIFY (atomic_kind_set)
139 402 : CALL get_qs_env(qs_env, dft_control=dft_control, atomic_kind_set=atomic_kind_set)
140 402 : nkind = SIZE(atomic_kind_set)
141 :
142 402 : NULLIFY (tddfpt_print_section)
143 402 : tddfpt_print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
144 402 : log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, "PROGRAM_BANNER", extension=".tddfptLog")
145 :
146 1264 : DO ikind = 1, nkind
147 862 : atomic_kind => atomic_kind_set(ikind)
148 862 : kind_param => stda_kernel%kind_param_set(ikind)%kind_param
149 : ! element symbol, kind_number, atomic number
150 : CALL get_atomic_kind(atomic_kind, &
151 : element_symbol=stda_kernel%kind_param_set(ikind)%kind_param%symbol, &
152 : kind_number=stda_kernel%kind_param_set(ikind)%kind_param%kind_number, &
153 1264 : z=stda_kernel%kind_param_set(ikind)%kind_param%z)
154 : END DO
155 :
156 402 : IF (stda_control%do_exchange) THEN ! option to switch off exchange
157 : ! HFx_fraction
158 370 : stda_kernel%do_exchange = .TRUE.
159 370 : stda_kernel%hfx_fraction = stda_control%hfx_fraction
160 : ELSE
161 32 : stda_kernel%do_exchange = .FALSE.
162 32 : stda_kernel%hfx_fraction = 0.0_dp
163 : END IF
164 :
165 : ! alpha and beta parameter
166 402 : IF (stda_control%mn_alpha < -98.0_dp) THEN
167 400 : IF (dft_control%qs_control%xtb) THEN
168 200 : stda_kernel%alpha_param = 2.0_dp
169 : ELSE
170 200 : stda_kernel%alpha_param = alpha(1) + stda_kernel%hfx_fraction*alpha(2)
171 : END IF
172 : ELSE
173 2 : stda_kernel%alpha_param = stda_control%mn_alpha
174 : END IF
175 402 : IF (stda_control%mn_beta < -98.0_dp) THEN
176 400 : IF (dft_control%qs_control%xtb) THEN
177 200 : stda_kernel%beta_param = 4.0_dp
178 : ELSE
179 200 : stda_kernel%beta_param = beta(1) + stda_kernel%hfx_fraction*beta(2)
180 : END IF
181 : ELSE
182 2 : stda_kernel%beta_param = stda_control%mn_beta
183 : END IF
184 :
185 : ! TD Filter
186 402 : stda_kernel%eps_td_filter = stda_control%eps_td_filter
187 :
188 1264 : DO ikind = 1, nkind
189 : ! hardness parameter
190 : stda_kernel%kind_param_set(ikind)%kind_param%hardness_param = &
191 862 : hardness(stda_kernel%kind_param_set(ikind)%kind_param%z)*2.0_dp/evolt
192 : ! rcut parameter
193 862 : eta = stda_kernel%kind_param_set(ikind)%kind_param%hardness_param
194 862 : fxx = 2.0_dp*eta**2*stda_control%coulomb_sr_eps
195 862 : fxx = 0.5_dp*(1.0_dp/fxx)**0.33333_dp
196 862 : rcut = stda_control%coulomb_sr_cut
197 1264 : stda_kernel%kind_param_set(ikind)%kind_param%rcut = MIN(rcut, fxx)
198 : END DO
199 :
200 402 : IF (log_unit > 0) THEN
201 201 : IF (.NOT. stda_kernel%do_exchange) THEN
202 16 : WRITE (log_unit, "(T2,A,T78,A3)") "sTDA| Exchange term is not used!"
203 : END IF
204 201 : WRITE (log_unit, "(T2,A,T71,F10.4)") "sTDA| HFX Fraction", stda_kernel%hfx_fraction
205 201 : WRITE (log_unit, "(T2,A,T71,F10.4)") "sTDA| Mataga-Nishimoto exponent (C)", stda_kernel%alpha_param
206 201 : WRITE (log_unit, "(T2,A,T71,F10.4)") "sTDA| Mataga-Nishimoto exponent (X)", stda_kernel%beta_param
207 201 : WRITE (log_unit, "(T2,A,T61,E20.8)") "sTDA| TD matrix filter", stda_kernel%eps_td_filter
208 : END IF
209 402 : CALL cp_print_key_finished_output(log_unit, logger, tddfpt_print_section, "PROGRAM_BANNER")
210 :
211 402 : END SUBROUTINE stda_init_param
212 : ! **************************************************************************************************
213 : !> \brief Allocate the sTDA environment
214 : !> \param qs_env ...
215 : !> \param stda_kernel ...
216 : !> \param n_ao ...
217 : !> \param nactive ...
218 : ! **************************************************************************************************
219 402 : SUBROUTINE allocate_stda_env(qs_env, stda_kernel, n_ao, nactive)
220 :
221 : TYPE(qs_environment_type), POINTER :: qs_env
222 : TYPE(stda_env_type) :: stda_kernel
223 : INTEGER, INTENT(IN) :: n_ao
224 : INTEGER, DIMENSION(:), INTENT(IN) :: nactive
225 :
226 : INTEGER :: ii, nkind
227 402 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
228 : TYPE(stda_kind_type), POINTER :: kind_param
229 :
230 402 : stda_kernel%hfx_fraction = 0.0_dp
231 402 : stda_kernel%alpha_param = 0.0_dp
232 402 : stda_kernel%beta_param = 0.0_dp
233 1206 : stda_kernel%nactive = 0
234 1206 : stda_kernel%nactive(1:2) = nactive(1:2)
235 1206 : stda_kernel%n_ao = n_ao
236 402 : NULLIFY (stda_kernel%kind_param_set)
237 :
238 : ! initialize stda_kind_parameters
239 402 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
240 402 : nkind = SIZE(qs_kind_set)
241 :
242 2068 : ALLOCATE (stda_kernel%kind_param_set(nkind))
243 1264 : DO ii = 1, nkind
244 862 : NULLIFY (kind_param)
245 862 : CALL allocate_stda_kind_param(kind_param)
246 1264 : stda_kernel%kind_param_set(ii)%kind_param => kind_param
247 : END DO
248 :
249 402 : END SUBROUTINE allocate_stda_env
250 : ! **************************************************************************************************
251 : !> \brief Deallocate the sTDA environment
252 : !> \param stda_kernel ...
253 : ! **************************************************************************************************
254 402 : SUBROUTINE deallocate_stda_env(stda_kernel)
255 :
256 : TYPE(stda_env_type) :: stda_kernel
257 :
258 : INTEGER :: ii
259 : TYPE(stda_kind_type), POINTER :: kind_param
260 :
261 : ! deallocate stda_kind_parameters
262 402 : IF (ASSOCIATED(stda_kernel%kind_param_set)) THEN
263 1264 : DO ii = 1, SIZE(stda_kernel%kind_param_set)
264 862 : kind_param => stda_kernel%kind_param_set(ii)%kind_param
265 1264 : CALL deallocate_stda_kind_param(kind_param)
266 : END DO
267 402 : DEALLOCATE (stda_kernel%kind_param_set)
268 : NULLIFY (stda_kernel%kind_param_set)
269 : END IF
270 :
271 402 : END SUBROUTINE deallocate_stda_env
272 : ! **************************************************************************************************
273 : !> \brief Allocate sTDA kind parameter
274 : !> \param kind_param ...
275 : ! **************************************************************************************************
276 862 : SUBROUTINE allocate_stda_kind_param(kind_param)
277 :
278 : TYPE(stda_kind_type), POINTER :: kind_param
279 :
280 862 : IF (ASSOCIATED(kind_param)) &
281 0 : CALL deallocate_stda_kind_param(kind_param)
282 :
283 862 : ALLOCATE (kind_param)
284 :
285 862 : kind_param%symbol = ""
286 862 : kind_param%z = 0
287 862 : kind_param%kind_number = 0
288 862 : kind_param%hardness_param = 0.0_dp
289 862 : kind_param%rcut = 0.0_dp
290 :
291 862 : END SUBROUTINE allocate_stda_kind_param
292 : ! **************************************************************************************************
293 : !> \brief Deallocate sTDA kind parameter
294 : !> \param kind_param ...
295 : ! **************************************************************************************************
296 862 : SUBROUTINE deallocate_stda_kind_param(kind_param)
297 :
298 : TYPE(stda_kind_type), POINTER :: kind_param
299 :
300 862 : CPASSERT(ASSOCIATED(kind_param))
301 862 : DEALLOCATE (kind_param)
302 : NULLIFY (kind_param)
303 :
304 862 : END SUBROUTINE deallocate_stda_kind_param
305 :
306 0 : END MODULE qs_tddfpt2_stda_types
|
