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 Calculation of the core Hamiltonian integral matrix <a|H|b> over
10 : !> Cartesian Gaussian-type functions.
11 : !>
12 : !> Nuclear potential energy:
13 : !>
14 : !> a) Allelectron calculation:
15 : !>
16 : !> erfc(r)
17 : !> <a|V|b> = -Z*<a|---------|b>
18 : !> r
19 : !>
20 : !> 1 - erf(r)
21 : !> = -Z*<a|------------|b>
22 : !> r
23 : !>
24 : !> 1 erf(r)
25 : !> = -Z*(<a|---|b> - <a|--------|b>)
26 : !> r r
27 : !>
28 : !> 1
29 : !> = -Z*(<a|---|b> - N*<ab||c>)
30 : !> r
31 : !>
32 : !> -Z
33 : !> = <a|---|b> + Z*N*<ab||c>
34 : !> r
35 : !> \_______/ \_____/
36 : !> | |
37 : !> nuclear coulomb
38 : !>
39 : !> b) Pseudopotential calculation (Goedecker, Teter and Hutter; GTH):
40 : !>
41 : !> <a|V|b> = <a|(V(local) + V(non-local))|b>
42 : !>
43 : !> = <a|(V(local)|b> + <a|V(non-local))|b>
44 : !>
45 : !> <a|V(local)|b> = <a|-Z(eff)*erf(SQRT(2)*alpha*r)/r +
46 : !> (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
47 : !> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
48 : !>
49 : !> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b>
50 : !> \par Literature
51 : !> S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B 54, 1703 (1996)
52 : !> C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B 58, 3641 (1998)
53 : !> M. Krack and M. Parrinello, Phys. Chem. Chem. Phys. 2, 2105 (2000)
54 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
55 : !> \par History
56 : !> - Joost VandeVondele (April 2003) : added LSD forces
57 : !> - Non-redundant calculation of the non-local part of the GTH PP
58 : !> (22.05.2003,MK)
59 : !> - New parallelization scheme (27.06.2003,MK)
60 : !> - OpenMP version (07.12.2003,JGH)
61 : !> - Binary search loop for VPPNL operators (09.01.2004,JGH,MK)
62 : !> - Refactoring of pseudopotential and nuclear attraction integrals (25.02.2009,JGH)
63 : !> - General refactoring (01.10.2010,JGH)
64 : !> - Refactoring related to the new kinetic energy and overlap routines (07.2014,JGH)
65 : !> - k-point functionality (07.2015,JGH)
66 : !> - Refactor for PP functionality (08.2025,JGH)
67 : !> \author Matthias Krack (14.09.2000,21.03.02)
68 : ! **************************************************************************************************
69 : MODULE qs_core_matrices
70 : USE atomic_kind_types, ONLY: atomic_kind_type,&
71 : get_atomic_kind_set
72 : USE cell_types, ONLY: cell_type
73 : USE core_ae, ONLY: build_core_ae
74 : USE core_ppl, ONLY: build_core_ppl
75 : USE core_ppnl, ONLY: build_core_ppnl
76 : USE cp_control_types, ONLY: dft_control_type
77 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
78 : dbcsr_iterator_blocks_left,&
79 : dbcsr_iterator_next_block,&
80 : dbcsr_iterator_start,&
81 : dbcsr_iterator_stop,&
82 : dbcsr_iterator_type,&
83 : dbcsr_p_type,&
84 : dbcsr_type
85 : USE cp_log_handling, ONLY: cp_get_default_logger,&
86 : cp_logger_get_default_unit_nr,&
87 : cp_logger_type
88 : USE ec_env_types, ONLY: energy_correction_type
89 : USE input_constants, ONLY: do_ppl_analytic,&
90 : rel_none,&
91 : rel_trans_atom
92 : USE kinds, ONLY: default_string_length,&
93 : dp
94 : USE kpoint_types, ONLY: get_kpoint_info,&
95 : kpoint_type
96 : USE lri_environment_types, ONLY: lri_environment_type
97 : USE mathlib, ONLY: det_3x3
98 : USE message_passing, ONLY: mp_para_env_type
99 : USE particle_types, ONLY: particle_type
100 : USE physcon, ONLY: pascal
101 : USE qs_environment_types, ONLY: get_qs_env,&
102 : qs_environment_type
103 : USE qs_force_types, ONLY: qs_force_type
104 : USE qs_kind_types, ONLY: get_qs_kind,&
105 : qs_kind_type
106 : USE qs_kinetic, ONLY: build_kinetic_matrix
107 : USE qs_ks_types, ONLY: get_ks_env,&
108 : qs_ks_env_type
109 : USE qs_linres_types, ONLY: dcdr_env_type
110 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
111 : USE virial_methods, ONLY: one_third_sum_diag
112 : USE virial_types, ONLY: virial_type
113 : #include "./base/base_uses.f90"
114 :
115 : IMPLICIT NONE
116 :
117 : PRIVATE
118 :
119 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_matrices'
120 :
121 : PUBLIC :: core_matrices, kinetic_energy_matrix
122 :
123 : CONTAINS
124 :
125 : ! **************************************************************************************************
126 : !> \brief ...
127 : !> \param qs_env ...
128 : !> \param matrix_h ...
129 : !> \param matrix_p ...
130 : !> \param calculate_forces ...
131 : !> \param nder ...
132 : !> \param ec_env ...
133 : !> \param dcdr_env ...
134 : !> \param ec_env_matrices ...
135 : !> \param basis_type ...
136 : !> \param debug_forces ...
137 : !> \param debug_stress ...
138 : !> \param atcore ...
139 : ! **************************************************************************************************
140 18931 : SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, ec_env, dcdr_env, &
141 66 : ec_env_matrices, basis_type, debug_forces, debug_stress, atcore)
142 :
143 : TYPE(qs_environment_type), POINTER :: qs_env
144 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
145 : LOGICAL, INTENT(IN) :: calculate_forces
146 : INTEGER, INTENT(IN) :: nder
147 : TYPE(energy_correction_type), OPTIONAL, POINTER :: ec_env
148 : TYPE(dcdr_env_type), OPTIONAL :: dcdr_env
149 : LOGICAL, INTENT(IN), OPTIONAL :: ec_env_matrices
150 : CHARACTER(LEN=*), OPTIONAL :: basis_type
151 : LOGICAL, INTENT(IN), OPTIONAL :: debug_forces, debug_stress
152 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atcore
153 :
154 : CHARACTER(LEN=default_string_length) :: my_basis_type
155 : INTEGER :: iounit, natom, nimages
156 18931 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
157 : LOGICAL :: all_present, debfor, debstr, my_gt_nl, &
158 : ppl_present, ppnl_present, use_lrigpw, &
159 : use_virial
160 : REAL(KIND=dp) :: eps_ppnl, fconv
161 : REAL(KIND=dp), DIMENSION(3) :: fodeb
162 : REAL(KIND=dp), DIMENSION(3, 3) :: stdeb
163 18931 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: deltaR
164 18931 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
165 : TYPE(cell_type), POINTER :: cell
166 : TYPE(cp_logger_type), POINTER :: logger
167 18931 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_hloc, matrix_ploc
168 : TYPE(dft_control_type), POINTER :: dft_control
169 : TYPE(kpoint_type), POINTER :: kpoints
170 : TYPE(lri_environment_type), POINTER :: lri_env
171 : TYPE(mp_para_env_type), POINTER :: para_env
172 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
173 18931 : POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
174 18931 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
175 18931 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
176 18931 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
177 : TYPE(qs_ks_env_type), POINTER :: ks_env
178 : TYPE(virial_type), POINTER :: virial
179 :
180 37862 : logger => cp_get_default_logger()
181 18931 : IF (logger%para_env%is_source()) THEN
182 9736 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
183 : ELSE
184 : iounit = -1
185 : END IF
186 :
187 18931 : NULLIFY (dft_control)
188 18931 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env)
189 :
190 : ! basis type
191 18931 : IF (PRESENT(basis_type)) THEN
192 802 : my_basis_type = basis_type
193 : ELSE
194 18129 : my_basis_type = "ORB"
195 : END IF
196 :
197 18931 : IF (PRESENT(dcdr_env) .AND. PRESENT(ec_env)) THEN
198 0 : CPABORT("core_matrices: conflicting options")
199 : END IF
200 :
201 : ! check size of atcore array
202 18931 : IF (PRESENT(atcore)) THEN
203 66 : CALL get_qs_env(qs_env=qs_env, natom=natom)
204 66 : CPASSERT(SIZE(atcore) >= natom)
205 : END IF
206 :
207 : ! check whether a gauge transformed version of the non-local potential part has to be used
208 18931 : my_gt_nl = .FALSE.
209 18931 : IF (qs_env%run_rtp) THEN
210 1238 : CPASSERT(ASSOCIATED(dft_control%rtp_control))
211 1238 : IF (dft_control%rtp_control%velocity_gauge) THEN
212 54 : my_gt_nl = dft_control%rtp_control%nl_gauge_transform
213 : END IF
214 : END IF
215 :
216 : ! prepare for k-points
217 18931 : nimages = dft_control%nimages
218 18931 : NULLIFY (cell_to_index)
219 18931 : IF (nimages > 1) THEN
220 310 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
221 310 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
222 310 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
223 310 : dft_control%qs_control%do_ppl_method = do_ppl_analytic
224 : END IF
225 :
226 : ! Possible dc/dr terms
227 18931 : IF (PRESENT(dcdr_env)) THEN
228 72 : deltaR => dcdr_env%delta_basis_function
229 72 : matrix_hloc(1:3, 1:1) => dcdr_env%matrix_hc(1:3)
230 : ELSE
231 18859 : NULLIFY (deltaR)
232 18859 : matrix_hloc => matrix_h
233 : END IF
234 18931 : matrix_ploc => matrix_p
235 :
236 : ! force analytic ppl calculation for GAPW methods
237 18931 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
238 2418 : dft_control%qs_control%do_ppl_method = do_ppl_analytic
239 : END IF
240 :
241 : ! force
242 18931 : NULLIFY (force)
243 18931 : IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force)
244 : ! virial
245 18931 : CALL get_qs_env(qs_env=qs_env, virial=virial)
246 18931 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
247 : ! lri/gpw
248 18931 : use_lrigpw = .FALSE.
249 :
250 : !debug
251 18931 : debfor = .FALSE.
252 18931 : IF (PRESENT(debug_forces)) debfor = debug_forces
253 1772 : debfor = debfor .AND. calculate_forces
254 18931 : debstr = .FALSE.
255 18931 : IF (PRESENT(debug_stress)) debstr = debug_stress
256 18931 : debstr = debstr .AND. use_virial
257 18931 : IF (debstr) THEN
258 50 : CALL get_qs_env(qs_env=qs_env, cell=cell)
259 50 : fconv = 1.0E-9_dp*pascal/cell%deth
260 : END IF
261 :
262 18931 : NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
263 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
264 18931 : particle_set=particle_set)
265 :
266 18931 : NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
267 18931 : IF (PRESENT(ec_env)) THEN
268 802 : sab_orb => ec_env%sab_orb
269 802 : sac_ae => ec_env%sac_ae
270 802 : sac_ppl => ec_env%sac_ppl
271 802 : sap_ppnl => ec_env%sap_ppnl
272 : ELSE
273 : CALL get_qs_env(qs_env=qs_env, &
274 : sab_orb=sab_orb, &
275 : sac_ae=sac_ae, &
276 : sac_ppl=sac_ppl, &
277 18129 : sap_ppnl=sap_ppnl)
278 : END IF
279 :
280 18931 : IF (PRESENT(ec_env) .AND. PRESENT(ec_env_matrices)) THEN
281 802 : IF (ec_env_matrices) THEN
282 346 : matrix_hloc => ec_env%matrix_h
283 346 : matrix_ploc => ec_env%matrix_p
284 : END IF
285 : END IF
286 :
287 : ! *** compute the nuclear attraction contribution to the core hamiltonian ***
288 18931 : all_present = ASSOCIATED(sac_ae)
289 18931 : IF (all_present) THEN
290 1076 : IF (PRESENT(dcdr_env)) THEN
291 0 : CPABORT("ECP/AE functionality for dcdr missing")
292 : END IF
293 1112 : IF (debfor) fodeb(1:3) = force(1)%all_potential(1:3, 1)
294 1076 : IF (debstr) stdeb = virial%pv_ppl
295 : CALL build_core_ae(matrix_hloc, matrix_ploc, force, virial, calculate_forces, use_virial, nder, &
296 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
297 2148 : nimages, cell_to_index, my_basis_type, atcore=atcore)
298 1076 : IF (debfor) THEN
299 48 : fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
300 12 : CALL para_env%sum(fodeb)
301 12 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHae ", fodeb
302 : END IF
303 1076 : IF (debstr) THEN
304 0 : stdeb = fconv*(virial%pv_ppl - stdeb)
305 0 : CALL para_env%sum(stdeb)
306 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
307 0 : 'STRESS| P*dHae ', one_third_sum_diag(stdeb), det_3x3(stdeb)
308 : END IF
309 : END IF
310 :
311 : ! *** compute the ppl contribution to the core hamiltonian ***
312 18931 : ppl_present = ASSOCIATED(sac_ppl)
313 18931 : IF (ppl_present) THEN
314 17925 : IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
315 17901 : CALL get_qs_env(qs_env, lri_env=lri_env)
316 17901 : IF (ASSOCIATED(lri_env)) THEN
317 86 : use_lrigpw = (dft_control%qs_control%lrigpw .AND. lri_env%ppl_ri)
318 : END IF
319 : IF (use_lrigpw) THEN
320 4 : IF (lri_env%exact_1c_terms) THEN
321 0 : CPABORT("not implemented")
322 : END IF
323 : ELSE
324 19067 : IF (debfor) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
325 18497 : IF (debstr) stdeb = virial%pv_ppl
326 : CALL build_core_ppl(matrix_hloc, matrix_ploc, force, &
327 : virial, calculate_forces, use_virial, nder, &
328 : qs_kind_set, atomic_kind_set, particle_set, &
329 : sab_orb, sac_ppl, nimages, cell_to_index, my_basis_type, &
330 35732 : deltaR=deltaR, atcore=atcore)
331 17897 : IF (debfor) THEN
332 1560 : fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
333 390 : CALL para_env%sum(fodeb)
334 390 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppl ", fodeb
335 : END IF
336 17897 : IF (debstr) THEN
337 650 : stdeb = fconv*(virial%pv_ppl - stdeb)
338 50 : CALL para_env%sum(stdeb)
339 50 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
340 25 : 'STRESS| P*dHppl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
341 : END IF
342 : END IF
343 : END IF
344 : END IF
345 :
346 : ! *** compute the ppnl contribution to the core hamiltonian ***
347 18931 : eps_ppnl = dft_control%qs_control%eps_ppnl
348 18931 : ppnl_present = ASSOCIATED(sap_ppnl)
349 18931 : IF (ppnl_present) THEN
350 15032 : IF (PRESENT(dcdr_env)) THEN
351 72 : matrix_hloc(1:3, 1:1) => dcdr_env%matrix_ppnl_1(1:3)
352 : END IF
353 15032 : IF (.NOT. my_gt_nl) THEN
354 16124 : IF (debfor) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
355 15596 : IF (debstr) stdeb = virial%pv_ppnl
356 : CALL build_core_ppnl(matrix_hloc, matrix_ploc, force, &
357 : virial, calculate_forces, use_virial, nder, &
358 : qs_kind_set, atomic_kind_set, particle_set, &
359 : sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, my_basis_type, &
360 29930 : deltaR=deltaR, atcore=atcore)
361 14996 : IF (debfor) THEN
362 1504 : fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
363 376 : CALL para_env%sum(fodeb)
364 376 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppnl ", fodeb
365 : END IF
366 14996 : IF (debstr) THEN
367 650 : stdeb = fconv*(virial%pv_ppnl - stdeb)
368 50 : CALL para_env%sum(stdeb)
369 50 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
370 25 : 'STRESS| P*dHppnl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
371 : END IF
372 : END IF
373 : END IF
374 :
375 18997 : END SUBROUTINE core_matrices
376 :
377 : ! **************************************************************************************************
378 : !> \brief Calculate kinetic energy matrix and possible relativistic correction
379 : !> \param qs_env ...
380 : !> \param matrixkp_t ...
381 : !> \param matrix_t ...
382 : !> \param matrix_p ...
383 : !> \param matrix_name ...
384 : !> \param calculate_forces ...
385 : !> \param nderivative ...
386 : !> \param sab_orb ...
387 : !> \param eps_filter ...
388 : !> \param basis_type ...
389 : !> \param debug_forces ...
390 : !> \param debug_stress ...
391 : !> \author Creation (21.08.2025,JGH)
392 : ! **************************************************************************************************
393 18471 : SUBROUTINE kinetic_energy_matrix(qs_env, matrixkp_t, matrix_t, matrix_p, &
394 : matrix_name, calculate_forces, nderivative, &
395 : sab_orb, eps_filter, basis_type, debug_forces, debug_stress)
396 :
397 : TYPE(qs_environment_type), POINTER :: qs_env
398 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
399 : POINTER :: matrixkp_t
400 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
401 : POINTER :: matrix_t
402 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
403 : POINTER :: matrix_p
404 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
405 : LOGICAL, INTENT(IN) :: calculate_forces
406 : INTEGER, INTENT(IN), OPTIONAL :: nderivative
407 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
408 : OPTIONAL, POINTER :: sab_orb
409 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter
410 : CHARACTER(LEN=*), OPTIONAL :: basis_type
411 : LOGICAL, INTENT(IN), OPTIONAL :: debug_forces, debug_stress
412 :
413 : CHARACTER(LEN=default_string_length) :: my_basis_type
414 : INTEGER :: ic, img, iounit, nimages
415 18471 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
416 : LOGICAL :: debfor, debstr, use_virial
417 : REAL(KIND=dp) :: fconv
418 : REAL(KIND=dp), DIMENSION(3) :: fodeb
419 : REAL(KIND=dp), DIMENSION(3, 3) :: stdeb
420 18471 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
421 : TYPE(cell_type), POINTER :: cell
422 : TYPE(cp_logger_type), POINTER :: logger
423 : TYPE(dft_control_type), POINTER :: dft_control
424 : TYPE(kpoint_type), POINTER :: kpoints
425 : TYPE(mp_para_env_type), POINTER :: para_env
426 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
427 18471 : POINTER :: sab_nl
428 18471 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
429 18471 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
430 : TYPE(qs_ks_env_type), POINTER :: ks_env
431 : TYPE(virial_type), POINTER :: virial
432 :
433 36942 : logger => cp_get_default_logger()
434 18471 : IF (logger%para_env%is_source()) THEN
435 9506 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
436 : ELSE
437 : iounit = -1
438 : END IF
439 :
440 18471 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env)
441 : ! is this a orbital-free method calculation
442 18471 : IF (dft_control%qs_control%ofgpw) RETURN
443 :
444 : ! Matrix images (kp)
445 18471 : nimages = dft_control%nimages
446 :
447 : ! basis type
448 18471 : IF (PRESENT(basis_type)) THEN
449 18187 : my_basis_type = basis_type
450 : ELSE
451 284 : my_basis_type = "ORB"
452 : END IF
453 :
454 18471 : debfor = .FALSE.
455 18471 : IF (PRESENT(debug_forces)) debfor = debug_forces
456 1772 : debfor = debfor .AND. calculate_forces
457 :
458 18471 : CALL get_qs_env(qs_env=qs_env, virial=virial)
459 18471 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
460 :
461 18471 : debstr = .FALSE.
462 18471 : IF (PRESENT(debug_stress)) debstr = debug_stress
463 20243 : debstr = debstr .AND. use_virial
464 1772 : IF (debstr) THEN
465 50 : CALL get_qs_env(qs_env=qs_env, cell=cell)
466 50 : fconv = 1.0E-9_dp*pascal/cell%deth
467 : END IF
468 :
469 18471 : NULLIFY (ks_env)
470 18471 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
471 18471 : NULLIFY (sab_nl)
472 18471 : IF (PRESENT(sab_orb)) THEN
473 18187 : sab_nl => sab_orb
474 : ELSE
475 284 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
476 : END IF
477 :
478 18471 : IF (calculate_forces) THEN
479 7697 : IF (SIZE(matrix_p, 1) == 2) THEN
480 2834 : DO img = 1, nimages
481 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
482 2834 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
483 : END DO
484 : END IF
485 : END IF
486 :
487 18471 : IF (debfor) THEN
488 402 : CALL get_qs_env(qs_env=qs_env, force=force)
489 1608 : fodeb(1:3) = force(1)%kinetic(1:3, 1)
490 : END IF
491 18471 : IF (debstr) THEN
492 650 : stdeb = virial%pv_ekinetic
493 : END IF
494 :
495 : ! T matrix
496 : CALL build_kinetic_matrix(ks_env, matrixkp_t=matrixkp_t, matrix_t=matrix_t, &
497 : matrix_name=matrix_name, &
498 : basis_type=my_basis_type, &
499 : sab_nl=sab_nl, &
500 : calculate_forces=calculate_forces, &
501 : matrixkp_p=matrix_p, &
502 : nderivative=nderivative, &
503 19211 : eps_filter=eps_filter)
504 :
505 18471 : IF (debfor) THEN
506 1608 : fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
507 402 : CALL para_env%sum(fodeb)
508 402 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dT ", fodeb
509 : END IF
510 18471 : IF (debstr) THEN
511 650 : stdeb = fconv*(virial%pv_ekinetic - stdeb)
512 50 : CALL para_env%sum(stdeb)
513 50 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
514 25 : 'STRESS| P*dT', one_third_sum_diag(stdeb), det_3x3(stdeb)
515 : END IF
516 :
517 18471 : IF (calculate_forces) THEN
518 7697 : IF (SIZE(matrix_p, 1) == 2) THEN
519 2834 : DO img = 1, nimages
520 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
521 2834 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
522 : END DO
523 : END IF
524 : END IF
525 :
526 : ! relativistic atomic correction to kinetic energy
527 18471 : IF (qs_env%rel_control%rel_method /= rel_none) THEN
528 58 : IF (qs_env%rel_control%rel_transformation == rel_trans_atom) THEN
529 58 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
530 58 : IF (nimages == 1) THEN
531 : ic = 1
532 : ELSE
533 4 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
534 4 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
535 4 : ic = cell_to_index(0, 0, 0)
536 : END IF
537 58 : IF (my_basis_type /= "ORB") THEN
538 0 : CPABORT("Basis mismatch for relativistic correction")
539 : END IF
540 58 : IF (PRESENT(matrixkp_t)) THEN
541 58 : CALL build_atomic_relmat(matrixkp_t(1, ic)%matrix, atomic_kind_set, qs_kind_set)
542 0 : ELSEIF (PRESENT(matrix_t)) THEN
543 0 : CALL build_atomic_relmat(matrix_t(1)%matrix, atomic_kind_set, qs_kind_set)
544 : END IF
545 : ELSE
546 0 : CPABORT("Relativistic corrections of this type are currently not implemented")
547 : END IF
548 : END IF
549 :
550 18471 : END SUBROUTINE kinetic_energy_matrix
551 :
552 : ! **************************************************************************************************
553 : !> \brief Adds atomic blocks of relativistic correction for the kinetic energy
554 : !> \param matrix_t ...
555 : !> \param atomic_kind_set ...
556 : !> \param qs_kind_set ...
557 : ! **************************************************************************************************
558 58 : SUBROUTINE build_atomic_relmat(matrix_t, atomic_kind_set, qs_kind_set)
559 : TYPE(dbcsr_type), POINTER :: matrix_t
560 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
561 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
562 :
563 : INTEGER :: iatom, ikind, jatom
564 58 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
565 58 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: reltmat, tblock
566 : TYPE(dbcsr_iterator_type) :: iter
567 :
568 58 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
569 :
570 58 : CALL dbcsr_iterator_start(iter, matrix_t)
571 221 : DO WHILE (dbcsr_iterator_blocks_left(iter))
572 163 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, tblock)
573 221 : IF (iatom == jatom) THEN
574 83 : ikind = kind_of(iatom)
575 83 : CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat)
576 192766 : IF (ASSOCIATED(reltmat)) tblock = tblock + reltmat
577 : END IF
578 : END DO
579 58 : CALL dbcsr_iterator_stop(iter)
580 :
581 116 : END SUBROUTINE build_atomic_relmat
582 :
583 : END MODULE qs_core_matrices
|
