(git:c09f6e5)
Loading...
Searching...
No Matches
qs_core_matrices.F
Go to the documentation of this file.
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! **************************************************************************************************
69MODULE 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
123CONTAINS
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 SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, ec_env, dcdr_env, &
141 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: deltar
164 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
165 TYPE(cell_type), POINTER :: cell
166 TYPE(cp_logger_type), POINTER :: logger
167 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 POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
174 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
175 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
176 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 logger => cp_get_default_logger()
181 IF (logger%para_env%is_source()) THEN
182 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
183 ELSE
184 iounit = -1
185 END IF
186
187 NULLIFY (dft_control)
188 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env)
189
190 ! basis type
191 IF (PRESENT(basis_type)) THEN
192 my_basis_type = basis_type
193 ELSE
194 my_basis_type = "ORB"
195 END IF
196
197 IF (PRESENT(dcdr_env) .AND. PRESENT(ec_env)) THEN
198 cpabort("core_matrices: conflicting options")
199 END IF
200
201 ! check size of atcore array
202 IF (PRESENT(atcore)) THEN
203 CALL get_qs_env(qs_env=qs_env, natom=natom)
204 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 my_gt_nl = .false.
209 IF (qs_env%run_rtp) THEN
210 cpassert(ASSOCIATED(dft_control%rtp_control))
211 IF (dft_control%rtp_control%velocity_gauge) THEN
212 my_gt_nl = dft_control%rtp_control%nl_gauge_transform
213 END IF
214 END IF
215
216 ! prepare for k-points
217 nimages = dft_control%nimages
218 NULLIFY (cell_to_index)
219 IF (nimages > 1) THEN
220 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
221 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
222 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
223 dft_control%qs_control%do_ppl_method = do_ppl_analytic
224 END IF
225
226 ! Possible dc/dr terms
227 IF (PRESENT(dcdr_env)) THEN
228 deltar => dcdr_env%delta_basis_function
229 matrix_hloc(1:3, 1:1) => dcdr_env%matrix_hc(1:3)
230 ELSE
231 NULLIFY (deltar)
232 matrix_hloc => matrix_h
233 END IF
234 matrix_ploc => matrix_p
235
236 ! force analytic ppl calculation for GAPW methods
237 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
238 dft_control%qs_control%do_ppl_method = do_ppl_analytic
239 END IF
240
241 ! force
242 NULLIFY (force)
243 IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force)
244 ! virial
245 CALL get_qs_env(qs_env=qs_env, virial=virial)
246 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
247 ! lri/gpw
248 use_lrigpw = .false.
249
250 !debug
251 debfor = .false.
252 IF (PRESENT(debug_forces)) debfor = debug_forces
253 debfor = debfor .AND. calculate_forces
254 debstr = .false.
255 IF (PRESENT(debug_stress)) debstr = debug_stress
256 debstr = debstr .AND. use_virial
257 IF (debstr) THEN
258 CALL get_qs_env(qs_env=qs_env, cell=cell)
259 fconv = 1.0e-9_dp*pascal/cell%deth
260 END IF
261
262 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 particle_set=particle_set)
265
266 NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
267 IF (PRESENT(ec_env)) THEN
268 sab_orb => ec_env%sab_orb
269 sac_ae => ec_env%sac_ae
270 sac_ppl => ec_env%sac_ppl
271 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 sap_ppnl=sap_ppnl)
278 END IF
279
280 IF (PRESENT(ec_env) .AND. PRESENT(ec_env_matrices)) THEN
281 IF (ec_env_matrices) THEN
282 matrix_hloc => ec_env%matrix_h
283 matrix_ploc => ec_env%matrix_p
284 END IF
285 END IF
286
287 ! *** compute the nuclear attraction contribution to the core hamiltonian ***
288 all_present = ASSOCIATED(sac_ae)
289 IF (all_present) THEN
290 IF (PRESENT(dcdr_env)) THEN
291 cpabort("ECP/AE functionality for dcdr missing")
292 END IF
293 IF (debfor) fodeb(1:3) = force(1)%all_potential(1:3, 1)
294 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 nimages, cell_to_index, my_basis_type, atcore=atcore)
298 IF (debfor) THEN
299 fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
300 CALL para_env%sum(fodeb)
301 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHae ", fodeb
302 END IF
303 IF (debstr) THEN
304 stdeb = fconv*(virial%pv_ppl - stdeb)
305 CALL para_env%sum(stdeb)
306 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
307 '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 ppl_present = ASSOCIATED(sac_ppl)
313 IF (ppl_present) THEN
314 IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
315 CALL get_qs_env(qs_env, lri_env=lri_env)
316 IF (ASSOCIATED(lri_env)) THEN
317 use_lrigpw = (dft_control%qs_control%lrigpw .AND. lri_env%ppl_ri)
318 END IF
319 IF (use_lrigpw) THEN
320 IF (lri_env%exact_1c_terms) THEN
321 cpabort("not implemented")
322 END IF
323 ELSE
324 IF (debfor) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
325 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 deltar=deltar, atcore=atcore)
331 IF (debfor) THEN
332 fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
333 CALL para_env%sum(fodeb)
334 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppl ", fodeb
335 END IF
336 IF (debstr) THEN
337 stdeb = fconv*(virial%pv_ppl - stdeb)
338 CALL para_env%sum(stdeb)
339 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
340 '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 eps_ppnl = dft_control%qs_control%eps_ppnl
348 ppnl_present = ASSOCIATED(sap_ppnl)
349 IF (ppnl_present) THEN
350 IF (PRESENT(dcdr_env)) THEN
351 matrix_hloc(1:3, 1:1) => dcdr_env%matrix_ppnl_1(1:3)
352 END IF
353 IF (.NOT. my_gt_nl) THEN
354 IF (debfor) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
355 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 deltar=deltar, atcore=atcore)
361 IF (debfor) THEN
362 fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
363 CALL para_env%sum(fodeb)
364 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppnl ", fodeb
365 END IF
366 IF (debstr) THEN
367 stdeb = fconv*(virial%pv_ppnl - stdeb)
368 CALL para_env%sum(stdeb)
369 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
370 'STRESS| P*dHppnl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
371 END IF
372 END IF
373 END IF
374
375 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 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 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 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 POINTER :: sab_nl
428 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
429 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 logger => cp_get_default_logger()
434 IF (logger%para_env%is_source()) THEN
435 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
436 ELSE
437 iounit = -1
438 END IF
439
440 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 IF (dft_control%qs_control%ofgpw) RETURN
443
444 ! Matrix images (kp)
445 nimages = dft_control%nimages
446
447 ! basis type
448 IF (PRESENT(basis_type)) THEN
449 my_basis_type = basis_type
450 ELSE
451 my_basis_type = "ORB"
452 END IF
453
454 debfor = .false.
455 IF (PRESENT(debug_forces)) debfor = debug_forces
456 debfor = debfor .AND. calculate_forces
457
458 CALL get_qs_env(qs_env=qs_env, virial=virial)
459 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
460
461 debstr = .false.
462 IF (PRESENT(debug_stress)) debstr = debug_stress
463 debstr = debstr .AND. use_virial
464 IF (debstr) THEN
465 CALL get_qs_env(qs_env=qs_env, cell=cell)
466 fconv = 1.0e-9_dp*pascal/cell%deth
467 END IF
468
469 NULLIFY (ks_env)
470 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
471 NULLIFY (sab_nl)
472 IF (PRESENT(sab_orb)) THEN
473 sab_nl => sab_orb
474 ELSE
475 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
476 END IF
477
478 IF (calculate_forces) THEN
479 IF (SIZE(matrix_p, 1) == 2) THEN
480 DO img = 1, nimages
481 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
482 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
483 END DO
484 END IF
485 END IF
486
487 IF (debfor) THEN
488 CALL get_qs_env(qs_env=qs_env, force=force)
489 fodeb(1:3) = force(1)%kinetic(1:3, 1)
490 END IF
491 IF (debstr) THEN
492 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 eps_filter=eps_filter)
504
505 IF (debfor) THEN
506 fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
507 CALL para_env%sum(fodeb)
508 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dT ", fodeb
509 END IF
510 IF (debstr) THEN
511 stdeb = fconv*(virial%pv_ekinetic - stdeb)
512 CALL para_env%sum(stdeb)
513 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
514 'STRESS| P*dT', one_third_sum_diag(stdeb), det_3x3(stdeb)
515 END IF
516
517 IF (calculate_forces) THEN
518 IF (SIZE(matrix_p, 1) == 2) THEN
519 DO img = 1, nimages
520 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
521 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 IF (qs_env%rel_control%rel_method /= rel_none) THEN
528 IF (qs_env%rel_control%rel_transformation == rel_trans_atom) THEN
529 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
530 IF (nimages == 1) THEN
531 ic = 1
532 ELSE
533 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
534 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
535 ic = cell_to_index(0, 0, 0)
536 END IF
537 IF (my_basis_type /= "ORB") THEN
538 cpabort("Basis mismatch for relativistic correction")
539 END IF
540 IF (PRESENT(matrixkp_t)) THEN
541 CALL build_atomic_relmat(matrixkp_t(1, ic)%matrix, atomic_kind_set, qs_kind_set)
542 ELSEIF (PRESENT(matrix_t)) THEN
543 CALL build_atomic_relmat(matrix_t(1)%matrix, atomic_kind_set, qs_kind_set)
544 END IF
545 ELSE
546 cpabort("Relativistic corrections of this type are currently not implemented")
547 END IF
548 END IF
549
550 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
565 REAL(kind=dp), DIMENSION(:, :), POINTER :: reltmat, tblock
566 TYPE(dbcsr_iterator_type) :: iter
567
568 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
569
570 CALL dbcsr_iterator_start(iter, matrix_t)
571 DO WHILE (dbcsr_iterator_blocks_left(iter))
572 CALL dbcsr_iterator_next_block(iter, iatom, jatom, tblock)
573 IF (iatom == jatom) THEN
574 ikind = kind_of(iatom)
575 CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat)
576 IF (ASSOCIATED(reltmat)) tblock = tblock + reltmat
577 END IF
578 END DO
579 CALL dbcsr_iterator_stop(iter)
580
581 END SUBROUTINE build_atomic_relmat
582
583END MODULE qs_core_matrices
mathlib::det_3x3
Definition mathlib.F:66
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
core_ae
Calculation of the nuclear attraction contribution to the core Hamiltonian <a|erfc|b> :we only calcul...
Definition core_ae.F:14
core_ae::build_core_ae
subroutine, public build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index, basis_type, atcore)
...
Definition core_ae.F:87
core_ppl
Calculation of the local pseudopotential contribution to the core Hamiltonian <a|V(local)|b> = <a|Sum...
Definition core_ppl.F:18
core_ppl::build_core_ppl
subroutine, public build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, nimages, cell_to_index, basis_type, deltar, atcore)
...
Definition core_ppl.F:97
core_ppnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition core_ppnl.F:15
core_ppnl::build_core_ppnl
subroutine, public build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, basis_type, deltar, matrix_l, atcore)
...
Definition core_ppnl.F:90
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
Definition cp_dbcsr_api.F:969
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:943
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1040
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1002
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:253
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
ec_env_types
Types needed for a for a Energy Correction.
Definition ec_env_types.F:14
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::rel_trans_atom
integer, parameter, public rel_trans_atom
Definition input_constants.F:615
input_constants::rel_none
integer, parameter, public rel_none
Definition input_constants.F:609
input_constants::do_ppl_analytic
integer, parameter, public do_ppl_analytic
Definition input_constants.F:275
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:365
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition lri_environment_types.F:18
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::pascal
real(kind=dp), parameter, public pascal
Definition physcon.F:174
qs_core_matrices
Calculation of the core Hamiltonian integral matrix <a|H|b> over Cartesian Gaussian-type functions.
Definition qs_core_matrices.F:69
qs_core_matrices::kinetic_energy_matrix
subroutine, public kinetic_energy_matrix(qs_env, matrixkp_t, matrix_t, matrix_p, matrix_name, calculate_forces, nderivative, sab_orb, eps_filter, basis_type, debug_forces, debug_stress)
Calculate kinetic energy matrix and possible relativistic correction.
Definition qs_core_matrices.F:396
qs_core_matrices::core_matrices
subroutine, public core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, ec_env, dcdr_env, ec_env_matrices, basis_type, debug_forces, debug_stress, atcore)
...
Definition qs_core_matrices.F:142
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:526
qs_force_types
Definition qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:466
qs_kinetic
Calculation of kinetic energy matrix and forces.
Definition qs_kinetic.F:15
qs_kinetic::build_kinetic_matrix
subroutine, public build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
Definition qs_kinetic.F:101
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition qs_ks_types.F:336
qs_linres_types
Type definitiona for linear response calculations.
Definition qs_linres_types.F:12
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
virial_methods
Definition virial_methods.F:12
virial_methods::one_third_sum_diag
pure real(kind=dp) function, public one_third_sum_diag(a)
...
Definition virial_methods.F:244
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:602
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
ec_env_types::energy_correction_type
Contains information on the energy correction functional for KG.
Definition ec_env_types.F:55
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:150
lri_environment_types::lri_environment_type
Definition lri_environment_types.F:271
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_linres_types::dcdr_env_type
Definition qs_linres_types.F:247
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
virial_types::virial_type
Definition virial_types.F:26