Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for calculating local energy and stress tensor
10 : !> \author JGH
11 : !> \par History
12 : !> - 07.2019 created
13 : ! **************************************************************************************************
14 : MODULE qs_local_properties
15 : USE bibliography, ONLY: Cohen2000,&
16 : Filippetti2000,&
17 : Rogers2002,&
18 : cite_reference
19 : USE cp_control_types, ONLY: dft_control_type
20 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
21 : dbcsr_p_type,&
22 : dbcsr_set,&
23 : dbcsr_type
24 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
25 : USE cp_log_handling, ONLY: cp_get_default_logger,&
26 : cp_logger_get_default_io_unit,&
27 : cp_logger_type
28 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
29 : section_vals_type
30 : USE kinds, ONLY: dp
31 : USE mathlib, ONLY: det_3x3
32 : USE pw_env_types, ONLY: pw_env_get,&
33 : pw_env_type
34 : USE pw_methods, ONLY: pw_axpy,&
35 : pw_copy,&
36 : pw_derive,&
37 : pw_integrate_function,&
38 : pw_multiply,&
39 : pw_transfer,&
40 : pw_zero
41 : USE pw_pool_types, ONLY: pw_pool_type
42 : USE pw_types, ONLY: pw_c1d_gs_type,&
43 : pw_r3d_rs_type
44 : USE qs_collocate_density, ONLY: calculate_rho_elec
45 : USE qs_core_energies, ONLY: calculate_ptrace
46 : USE qs_energy_types, ONLY: qs_energy_type
47 : USE qs_environment_types, ONLY: get_qs_env,&
48 : qs_environment_type
49 : USE qs_ks_methods, ONLY: calc_rho_tot_gspace
50 : USE qs_ks_types, ONLY: qs_ks_env_type,&
51 : set_ks_env
52 : USE qs_matrix_w, ONLY: compute_matrix_w
53 : USE qs_rho_types, ONLY: qs_rho_get,&
54 : qs_rho_type
55 : USE qs_vxc, ONLY: qs_xc_density
56 : USE virial_types, ONLY: virial_type
57 : #include "./base/base_uses.f90"
58 :
59 : IMPLICIT NONE
60 :
61 : PRIVATE
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_local_properties'
64 :
65 : PUBLIC :: qs_local_energy, qs_local_stress
66 :
67 : ! **************************************************************************************************
68 :
69 : CONTAINS
70 :
71 : ! **************************************************************************************************
72 : !> \brief Routine to calculate the local energy
73 : !> \param qs_env the qs_env to update
74 : !> \param energy_density ...
75 : !> \par History
76 : !> 07.2019 created
77 : !> \author JGH
78 : ! **************************************************************************************************
79 34 : SUBROUTINE qs_local_energy(qs_env, energy_density)
80 : TYPE(qs_environment_type), POINTER :: qs_env
81 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: energy_density
82 :
83 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_local_energy'
84 :
85 : INTEGER :: handle, img, iounit, ispin, nimages, &
86 : nkind, nspins
87 : LOGICAL :: gapw, gapw_xc, local_soft_only
88 : REAL(KIND=dp) :: eban, eband, eh, exc, ovol
89 : TYPE(cp_logger_type), POINTER :: logger
90 34 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
91 34 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_w, rho_ao_kp
92 : TYPE(dbcsr_type), POINTER :: matrix
93 : TYPE(dft_control_type), POINTER :: dft_control
94 : TYPE(pw_c1d_gs_type) :: edens_g
95 : TYPE(pw_c1d_gs_type), POINTER :: rho_core, rho_tot_gspace
96 : TYPE(pw_env_type), POINTER :: pw_env
97 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
98 : TYPE(pw_r3d_rs_type) :: band_density, edens_r, hartree_density, &
99 : xc_density
100 34 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
101 : TYPE(pw_r3d_rs_type), POINTER :: rho_tot_rspace, v_hartree_rspace
102 : TYPE(qs_energy_type), POINTER :: energy
103 : TYPE(qs_ks_env_type), POINTER :: ks_env
104 : TYPE(qs_rho_type), POINTER :: rho, rho_struct
105 : TYPE(section_vals_type), POINTER :: input, xc_section
106 :
107 34 : CALL timeset(routineN, handle)
108 :
109 34 : CALL cite_reference(Cohen2000)
110 :
111 34 : CPASSERT(ASSOCIATED(qs_env))
112 34 : logger => cp_get_default_logger()
113 34 : iounit = cp_logger_get_default_io_unit()
114 :
115 : ! Check for GAPW method : additional terms for local densities
116 34 : CALL get_qs_env(qs_env, nkind=nkind, dft_control=dft_control)
117 34 : gapw = dft_control%qs_control%gapw
118 34 : gapw_xc = dft_control%qs_control%gapw_xc
119 34 : local_soft_only = dft_control%do_admm .OR. gapw .OR. gapw_xc
120 :
121 34 : nimages = dft_control%nimages
122 34 : nspins = dft_control%nspins
123 :
124 : ! get working arrays
125 34 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
126 34 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
127 34 : CALL auxbas_pw_pool%create_pw(band_density)
128 34 : CALL auxbas_pw_pool%create_pw(hartree_density)
129 34 : CALL auxbas_pw_pool%create_pw(xc_density)
130 :
131 : ! w matrix
132 34 : CALL get_qs_env(qs_env, matrix_w_kp=matrix_w)
133 34 : IF (.NOT. ASSOCIATED(matrix_w)) THEN
134 : CALL get_qs_env(qs_env, &
135 : ks_env=ks_env, &
136 4 : matrix_s_kp=matrix_s)
137 4 : matrix => matrix_s(1, 1)%matrix
138 4 : CALL dbcsr_allocate_matrix_set(matrix_w, nspins, nimages)
139 8 : DO ispin = 1, nspins
140 12 : DO img = 1, nimages
141 4 : ALLOCATE (matrix_w(ispin, img)%matrix)
142 4 : CALL dbcsr_copy(matrix_w(ispin, img)%matrix, matrix, name="W MATRIX")
143 8 : CALL dbcsr_set(matrix_w(ispin, img)%matrix, 0.0_dp)
144 : END DO
145 : END DO
146 4 : CALL set_ks_env(ks_env, matrix_w_kp=matrix_w)
147 : END IF
148 : ! band structure energy density
149 34 : CALL compute_matrix_w(qs_env, .TRUE.)
150 34 : CALL get_qs_env(qs_env, ks_env=ks_env, matrix_w_kp=matrix_w)
151 34 : CALL auxbas_pw_pool%create_pw(edens_r)
152 34 : CALL auxbas_pw_pool%create_pw(edens_g)
153 34 : CALL pw_zero(band_density)
154 68 : DO ispin = 1, nspins
155 34 : rho_ao => matrix_w(ispin, :)
156 : CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
157 : rho=edens_r, &
158 : rho_gspace=edens_g, &
159 34 : ks_env=ks_env, soft_valid=(gapw .OR. gapw_xc))
160 68 : CALL pw_axpy(edens_r, band_density)
161 : END DO
162 34 : CALL auxbas_pw_pool%give_back_pw(edens_r)
163 34 : CALL auxbas_pw_pool%give_back_pw(edens_g)
164 :
165 : ! Hartree energy density correction = -0.5 * V_H(r) * [rho(r) - rho_core(r)]
166 34 : ALLOCATE (rho_tot_gspace, rho_tot_rspace)
167 34 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
168 34 : CALL auxbas_pw_pool%create_pw(rho_tot_rspace)
169 34 : NULLIFY (rho_core)
170 : CALL get_qs_env(qs_env, &
171 : v_hartree_rspace=v_hartree_rspace, &
172 34 : rho_core=rho_core, rho=rho)
173 34 : CALL qs_rho_get(rho, rho_r=rho_r)
174 34 : IF (ASSOCIATED(rho_core)) THEN
175 32 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
176 32 : CALL pw_transfer(rho_core, rho_tot_rspace)
177 : ELSE
178 2 : CALL pw_zero(rho_tot_rspace)
179 : END IF
180 68 : DO ispin = 1, nspins
181 68 : CALL pw_axpy(rho_r(ispin), rho_tot_rspace, alpha=-1.0_dp)
182 : END DO
183 34 : CALL pw_zero(hartree_density)
184 34 : ovol = 0.5_dp/hartree_density%pw_grid%dvol
185 34 : CALL pw_multiply(hartree_density, v_hartree_rspace, rho_tot_rspace, alpha=ovol)
186 34 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
187 34 : CALL auxbas_pw_pool%give_back_pw(rho_tot_rspace)
188 34 : DEALLOCATE (rho_tot_gspace, rho_tot_rspace)
189 :
190 34 : IF (dft_control%do_admm) THEN
191 : CALL cp_warn(__LOCATION__, &
192 0 : "ADMM local energy density contains only the regular-grid contribution")
193 : END IF
194 34 : IF (gapw_xc .OR. gapw) THEN
195 : CALL cp_warn(__LOCATION__, &
196 2 : "GAPW/GAPW_XC local energy density contains only the soft regular-grid part")
197 : END IF
198 : ! XC energy density correction = E_xc(r) - V_xc(r)*rho(r)
199 34 : CALL get_qs_env(qs_env, input=input)
200 34 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
201 34 : CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
202 : !
203 34 : CALL qs_xc_density(ks_env, rho_struct, xc_section, xc_ener=xc_density)
204 : !
205 : ! energies
206 34 : CALL get_qs_env(qs_env=qs_env, energy=energy)
207 34 : eban = pw_integrate_function(band_density)
208 34 : eh = pw_integrate_function(hartree_density)
209 34 : exc = pw_integrate_function(xc_density)
210 :
211 : ! band energy
212 34 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks)
213 34 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
214 34 : CALL calculate_ptrace(matrix_ks, rho_ao_kp, eband, nspins, .TRUE.)
215 :
216 : ! get full density
217 34 : CALL pw_copy(band_density, energy_density)
218 34 : CALL pw_axpy(hartree_density, energy_density)
219 34 : CALL pw_axpy(xc_density, energy_density)
220 :
221 34 : IF (iounit > 0) THEN
222 17 : WRITE (UNIT=iounit, FMT="(/,T3,A)") REPEAT("=", 78)
223 17 : IF (local_soft_only) THEN
224 : WRITE (UNIT=iounit, FMT="(T4,A)") &
225 1 : "Local Energy Calculation (regular-grid contribution only)"
226 1 : WRITE (UNIT=iounit, FMT="(T4,A,T52,A,T68,A)") "Energy Component", "Reference", "Local Grid"
227 : ELSE
228 16 : WRITE (UNIT=iounit, FMT="(T4,A,T52,A,T75,A)") "Local Energy Calculation", "GPW/GAPW", "Local"
229 : END IF
230 17 : WRITE (UNIT=iounit, FMT="(T4,A,T45,F15.8,T65,F15.8)") "Band Energy", eband, eban
231 17 : WRITE (UNIT=iounit, FMT="(T4,A,T65,F15.8)") "Hartree Energy Correction", eh
232 17 : WRITE (UNIT=iounit, FMT="(T4,A,T65,F15.8)") "XC Energy Correction", exc
233 17 : IF (local_soft_only) THEN
234 1 : WRITE (UNIT=iounit, FMT="(T4,A,T45,F15.8)") "Reference Total Energy", energy%total
235 1 : WRITE (UNIT=iounit, FMT="(T4,A,T65,F15.8)") "Local Grid Energy Sum", &
236 2 : eban + eh + exc + energy%core_overlap + energy%core_self + energy%dispersion
237 : ELSE
238 16 : WRITE (UNIT=iounit, FMT="(T4,A,T45,F15.8,T65,F15.8)") "Total Energy", &
239 32 : energy%total, eban + eh + exc + energy%core_overlap + energy%core_self + energy%dispersion
240 : END IF
241 17 : WRITE (UNIT=iounit, FMT="(T3,A)") REPEAT("=", 78)
242 : END IF
243 :
244 : ! return temp arrays
245 34 : CALL auxbas_pw_pool%give_back_pw(band_density)
246 34 : CALL auxbas_pw_pool%give_back_pw(hartree_density)
247 34 : CALL auxbas_pw_pool%give_back_pw(xc_density)
248 :
249 34 : CALL timestop(handle)
250 :
251 34 : END SUBROUTINE qs_local_energy
252 :
253 : ! **************************************************************************************************
254 : !> \brief Routine to calculate the local stress
255 : !> \param qs_env the qs_env to update
256 : !> \param stress_tensor ...
257 : !> \param beta ...
258 : !> \par History
259 : !> 07.2019 created
260 : !> \author JGH
261 : ! **************************************************************************************************
262 30 : SUBROUTINE qs_local_stress(qs_env, stress_tensor, beta)
263 : TYPE(qs_environment_type), POINTER :: qs_env
264 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), &
265 : INTENT(INOUT), OPTIONAL :: stress_tensor
266 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: beta
267 :
268 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_local_stress'
269 :
270 : INTEGER :: handle, i, iounit, j, nimages, nkind, &
271 : nspins
272 : LOGICAL :: do_stress, gapw, gapw_xc, use_virial
273 : REAL(KIND=dp) :: my_beta
274 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc
275 : TYPE(cp_logger_type), POINTER :: logger
276 : TYPE(dft_control_type), POINTER :: dft_control
277 : TYPE(pw_c1d_gs_type) :: v_hartree_gspace
278 120 : TYPE(pw_c1d_gs_type), DIMENSION(3) :: efield
279 : TYPE(pw_env_type), POINTER :: pw_env
280 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
281 : TYPE(pw_r3d_rs_type) :: xc_density
282 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
283 : TYPE(qs_ks_env_type), POINTER :: ks_env
284 : TYPE(qs_rho_type), POINTER :: rho_struct
285 : TYPE(section_vals_type), POINTER :: input, xc_section
286 : TYPE(virial_type), POINTER :: virial
287 :
288 30 : CALL cp_warn(__LOCATION__, "Local Stress Tensor code is not working, skipping")
289 30 : RETURN
290 :
291 : CALL timeset(routineN, handle)
292 :
293 : CALL cite_reference(Filippetti2000)
294 : CALL cite_reference(Rogers2002)
295 :
296 : CPASSERT(ASSOCIATED(qs_env))
297 :
298 : IF (PRESENT(stress_tensor)) THEN
299 : do_stress = .TRUE.
300 : ELSE
301 : do_stress = .FALSE.
302 : END IF
303 : IF (PRESENT(beta)) THEN
304 : my_beta = beta
305 : ELSE
306 : my_beta = 0.0_dp
307 : END IF
308 :
309 : logger => cp_get_default_logger()
310 : iounit = cp_logger_get_default_io_unit()
311 :
312 : !!!!!!
313 : CALL cp_warn(__LOCATION__, "Local Stress Tensor code is not tested")
314 : !!!!!!
315 :
316 : ! Check for GAPW method : additional terms for local densities
317 : CALL get_qs_env(qs_env, nkind=nkind, dft_control=dft_control)
318 : gapw = dft_control%qs_control%gapw
319 : gapw_xc = dft_control%qs_control%gapw_xc
320 :
321 : nimages = dft_control%nimages
322 : nspins = dft_control%nspins
323 :
324 : ! get working arrays
325 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
326 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
327 : CALL auxbas_pw_pool%create_pw(xc_density)
328 :
329 : ! init local stress tensor
330 : IF (do_stress) THEN
331 : DO i = 1, 3
332 : DO j = 1, 3
333 : CALL pw_zero(stress_tensor(i, j))
334 : END DO
335 : END DO
336 : END IF
337 :
338 : IF (dft_control%do_admm) THEN
339 : CALL cp_warn(__LOCATION__, &
340 : "ADMM local stress density contains only the regular-grid contribution")
341 : END IF
342 : IF (gapw_xc .OR. gapw) THEN
343 : CALL cp_warn(__LOCATION__, &
344 : "GAPW/GAPW_XC local stress density contains only the soft regular-grid part")
345 : END IF
346 : ! XC energy density correction = E_xc(r) - V_xc(r)*rho(r)
347 : CALL get_qs_env(qs_env, ks_env=ks_env, input=input, rho=rho_struct)
348 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
349 : !
350 : CALL qs_xc_density(ks_env, rho_struct, xc_section, xc_ener=xc_density)
351 :
352 : ! Electrical field terms
353 : CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace)
354 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
355 : CALL pw_transfer(v_hartree_rspace, v_hartree_gspace)
356 : DO i = 1, 3
357 : CALL auxbas_pw_pool%create_pw(efield(i))
358 : CALL pw_copy(v_hartree_gspace, efield(i))
359 : END DO
360 : CALL pw_derive(efield(1), [1, 0, 0])
361 : CALL pw_derive(efield(2), [0, 1, 0])
362 : CALL pw_derive(efield(3), [0, 0, 1])
363 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
364 : DO i = 1, 3
365 : CALL auxbas_pw_pool%give_back_pw(efield(i))
366 : END DO
367 :
368 : pv_loc = 0.0_dp
369 :
370 : CALL get_qs_env(qs_env=qs_env, virial=virial)
371 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
372 : IF (.NOT. use_virial) THEN
373 : CALL cp_warn(__LOCATION__, "Local stress should be used with standard stress calculation.")
374 : END IF
375 : IF (iounit > 0 .AND. use_virial) THEN
376 : WRITE (UNIT=iounit, FMT="(/,T3,A)") REPEAT("=", 78)
377 : WRITE (UNIT=iounit, FMT="(T4,A)") "Local Stress Calculation"
378 : WRITE (UNIT=iounit, FMT="(T42,A,T64,A)") " 1/3 Trace", " Determinant"
379 : WRITE (UNIT=iounit, FMT="(T4,A,T42,F16.8,T64,F16.8)") "Total Stress", &
380 : (pv_loc(1, 1) + pv_loc(2, 2) + pv_loc(3, 3))/3.0_dp, det_3x3(pv_loc)
381 : WRITE (UNIT=iounit, FMT="(T3,A)") REPEAT("=", 78)
382 : END IF
383 :
384 : CALL timestop(handle)
385 :
386 30 : END SUBROUTINE qs_local_stress
387 :
388 : ! **************************************************************************************************
389 :
390 : END MODULE qs_local_properties
|
