(git:b3d087c)
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-2026 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 ext_kpoints ...
136!> \param basis_type ...
137!> \param debug_forces ...
138!> \param debug_stress ...
139!> \param atcore ...
140! **************************************************************************************************
141 SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, ec_env, dcdr_env, &
142 ec_env_matrices, ext_kpoints, basis_type, &
143 debug_forces, debug_stress, atcore)
144
145 TYPE(qs_environment_type), POINTER :: qs_env
146 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
147 LOGICAL, INTENT(IN) :: calculate_forces
148 INTEGER, INTENT(IN) :: nder
149 TYPE(energy_correction_type), OPTIONAL, POINTER :: ec_env
150 TYPE(dcdr_env_type), OPTIONAL :: dcdr_env
151 LOGICAL, INTENT(IN), OPTIONAL :: ec_env_matrices
152 TYPE(kpoint_type), OPTIONAL, POINTER :: ext_kpoints
153 CHARACTER(LEN=*), OPTIONAL :: basis_type
154 LOGICAL, INTENT(IN), OPTIONAL :: debug_forces, debug_stress
155 REAL(kind=dp), DIMENSION(:), OPTIONAL :: atcore
156
157 CHARACTER(LEN=default_string_length) :: my_basis_type
158 INTEGER :: iounit, natom, nimages
159 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
160 LOGICAL :: all_present, debfor, debstr, ext_kp, &
161 my_gt_nl, ppl_present, ppnl_present, &
162 use_lrigpw, use_virial
163 REAL(kind=dp) :: eps_ppnl, fconv
164 REAL(kind=dp), DIMENSION(3) :: fodeb
165 REAL(kind=dp), DIMENSION(3, 3) :: stdeb
166 REAL(kind=dp), DIMENSION(:, :), POINTER :: deltar
167 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
168 TYPE(cell_type), POINTER :: cell
169 TYPE(cp_logger_type), POINTER :: logger
170 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_hloc, matrix_ploc
171 TYPE(dft_control_type), POINTER :: dft_control
172 TYPE(kpoint_type), POINTER :: kpoints
173 TYPE(lri_environment_type), POINTER :: lri_env
174 TYPE(mp_para_env_type), POINTER :: para_env
175 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
176 POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
177 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
178 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
179 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
180 TYPE(qs_ks_env_type), POINTER :: ks_env
181 TYPE(virial_type), POINTER :: virial
182
183 logger => cp_get_default_logger()
184 IF (logger%para_env%is_source()) THEN
185 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
186 ELSE
187 iounit = -1
188 END IF
189
190 ext_kp = .false.
191 IF (PRESENT(ext_kpoints)) THEN
192 IF (ASSOCIATED(ext_kpoints)) ext_kp = .true.
193 END IF
194
195 NULLIFY (dft_control)
196 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env)
197
198 ! basis type
199 IF (PRESENT(basis_type)) THEN
200 my_basis_type = basis_type
201 ELSE
202 my_basis_type = "ORB"
203 END IF
204
205 IF (PRESENT(dcdr_env) .AND. PRESENT(ec_env)) THEN
206 cpabort("core_matrices: conflicting options")
207 END IF
208
209 ! check size of atcore array
210 IF (PRESENT(atcore)) THEN
211 CALL get_qs_env(qs_env=qs_env, natom=natom)
212 cpassert(SIZE(atcore) >= natom)
213 END IF
214
215 ! check whether a gauge transformed version of the non-local potential part has to be used
216 my_gt_nl = .false.
217 IF (qs_env%run_rtp) THEN
218 cpassert(ASSOCIATED(dft_control%rtp_control))
219 IF (dft_control%rtp_control%velocity_gauge) THEN
220 my_gt_nl = dft_control%rtp_control%nl_gauge_transform
221 END IF
222 END IF
223
224 ! prepare for k-points
225 NULLIFY (cell_to_index)
226 IF (ext_kp) THEN
227 nimages = SIZE(ext_kpoints%index_to_cell, 2)
228 CALL get_kpoint_info(kpoint=ext_kpoints, cell_to_index=cell_to_index)
229 dft_control%qs_control%do_ppl_method = do_ppl_analytic
230 ELSE
231 nimages = dft_control%nimages
232 IF (nimages > 1) THEN
233 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
234 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
235 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
236 dft_control%qs_control%do_ppl_method = do_ppl_analytic
237 END IF
238 END IF
239
240 ! Possible dc/dr terms
241 IF (PRESENT(dcdr_env)) THEN
242 deltar => dcdr_env%delta_basis_function
243 matrix_hloc(1:3, 1:1) => dcdr_env%matrix_hc(1:3)
244 ELSE
245 NULLIFY (deltar)
246 matrix_hloc => matrix_h
247 END IF
248 matrix_ploc => matrix_p
249
250 ! force analytic ppl calculation for GAPW methods
251 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
252 dft_control%qs_control%do_ppl_method = do_ppl_analytic
253 END IF
254
255 ! force
256 NULLIFY (force)
257 IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force)
258 ! virial
259 CALL get_qs_env(qs_env=qs_env, virial=virial)
260 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
261 ! lri/gpw
262 use_lrigpw = .false.
263
264 !debug
265 debfor = .false.
266 IF (PRESENT(debug_forces)) debfor = debug_forces
267 debfor = debfor .AND. calculate_forces
268 debstr = .false.
269 IF (PRESENT(debug_stress)) debstr = debug_stress
270 debstr = debstr .AND. use_virial
271 IF (debstr) THEN
272 CALL get_qs_env(qs_env=qs_env, cell=cell)
273 fconv = 1.0e-9_dp*pascal/cell%deth
274 END IF
275
276 NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
277 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
278 particle_set=particle_set)
279
280 NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
281 IF (PRESENT(ec_env)) THEN
282 sab_orb => ec_env%sab_orb
283 sac_ae => ec_env%sac_ae
284 sac_ppl => ec_env%sac_ppl
285 sap_ppnl => ec_env%sap_ppnl
286 ELSE
287 CALL get_qs_env(qs_env=qs_env, &
288 sab_orb=sab_orb, &
289 sac_ae=sac_ae, &
290 sac_ppl=sac_ppl, &
291 sap_ppnl=sap_ppnl)
292 END IF
293
294 IF (PRESENT(ec_env) .AND. PRESENT(ec_env_matrices)) THEN
295 IF (ec_env_matrices) THEN
296 matrix_hloc => ec_env%matrix_h
297 matrix_ploc => ec_env%matrix_p
298 END IF
299 END IF
300
301 ! *** compute the nuclear attraction contribution to the core hamiltonian ***
302 all_present = ASSOCIATED(sac_ae)
303 IF (all_present) THEN
304 IF (PRESENT(dcdr_env)) THEN
305 cpabort("ECP/AE functionality for dcdr missing")
306 END IF
307 IF (debfor) fodeb(1:3) = force(1)%all_potential(1:3, 1)
308 IF (debstr) stdeb = virial%pv_ppl
309 CALL build_core_ae(matrix_hloc, matrix_ploc, force, virial, calculate_forces, use_virial, nder, &
310 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
311 nimages, cell_to_index, my_basis_type, atcore=atcore)
312 IF (debfor) THEN
313 fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
314 CALL para_env%sum(fodeb)
315 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHae ", fodeb
316 END IF
317 IF (debstr) THEN
318 stdeb = fconv*(virial%pv_ppl - stdeb)
319 CALL para_env%sum(stdeb)
320 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
321 'STRESS| P*dHae ', one_third_sum_diag(stdeb), det_3x3(stdeb)
322 END IF
323 END IF
324
325 ! *** compute the ppl contribution to the core hamiltonian ***
326 ppl_present = ASSOCIATED(sac_ppl)
327 IF (ppl_present) THEN
328 IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
329 CALL get_qs_env(qs_env, lri_env=lri_env)
330 IF (ASSOCIATED(lri_env)) THEN
331 use_lrigpw = (dft_control%qs_control%lrigpw .AND. lri_env%ppl_ri)
332 END IF
333 IF (use_lrigpw) THEN
334 IF (lri_env%exact_1c_terms) THEN
335 cpabort("not implemented")
336 END IF
337 ELSE
338 IF (debfor) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
339 IF (debstr) stdeb = virial%pv_ppl
340 CALL build_core_ppl(matrix_hloc, matrix_ploc, force, &
341 virial, calculate_forces, use_virial, nder, &
342 qs_kind_set, atomic_kind_set, particle_set, &
343 sab_orb, sac_ppl, nimages, cell_to_index, my_basis_type, &
344 deltar=deltar, atcore=atcore)
345 IF (debfor) THEN
346 fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
347 CALL para_env%sum(fodeb)
348 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppl ", fodeb
349 END IF
350 IF (debstr) THEN
351 stdeb = fconv*(virial%pv_ppl - stdeb)
352 CALL para_env%sum(stdeb)
353 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
354 'STRESS| P*dHppl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
355 END IF
356 END IF
357 END IF
358 END IF
359
360 ! *** compute the ppnl contribution to the core hamiltonian ***
361 eps_ppnl = dft_control%qs_control%eps_ppnl
362 ppnl_present = ASSOCIATED(sap_ppnl)
363 IF (ppnl_present) THEN
364 IF (PRESENT(dcdr_env)) THEN
365 matrix_hloc(1:3, 1:1) => dcdr_env%matrix_ppnl_1(1:3)
366 END IF
367 IF (.NOT. my_gt_nl) THEN
368 IF (debfor) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
369 IF (debstr) stdeb = virial%pv_ppnl
370 CALL build_core_ppnl(matrix_hloc, matrix_ploc, force, &
371 virial, calculate_forces, use_virial, nder, &
372 qs_kind_set, atomic_kind_set, particle_set, &
373 sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, my_basis_type, &
374 deltar=deltar, atcore=atcore)
375 IF (debfor) THEN
376 fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
377 CALL para_env%sum(fodeb)
378 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppnl ", fodeb
379 END IF
380 IF (debstr) THEN
381 stdeb = fconv*(virial%pv_ppnl - stdeb)
382 CALL para_env%sum(stdeb)
383 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
384 'STRESS| P*dHppnl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
385 END IF
386 END IF
387 END IF
388
389 END SUBROUTINE core_matrices
390
391! **************************************************************************************************
392!> \brief Calculate kinetic energy matrix and possible relativistic correction
393!> \param qs_env ...
394!> \param matrixkp_t ...
395!> \param matrix_t ...
396!> \param matrix_p ...
397!> \param ext_kpoints ...
398!> \param matrix_name ...
399!> \param calculate_forces ...
400!> \param nderivative ...
401!> \param sab_orb ...
402!> \param eps_filter ...
403!> \param basis_type ...
404!> \param debug_forces ...
405!> \param debug_stress ...
406!> \author Creation (21.08.2025,JGH)
407! **************************************************************************************************
408 SUBROUTINE kinetic_energy_matrix(qs_env, matrixkp_t, matrix_t, matrix_p, ext_kpoints, &
409 matrix_name, calculate_forces, nderivative, &
410 sab_orb, eps_filter, basis_type, debug_forces, debug_stress)
411
412 TYPE(qs_environment_type), POINTER :: qs_env
413 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
414 POINTER :: matrixkp_t
415 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
416 POINTER :: matrix_t
417 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
418 POINTER :: matrix_p
419 TYPE(kpoint_type), OPTIONAL, POINTER :: ext_kpoints
420 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
421 LOGICAL, INTENT(IN) :: calculate_forces
422 INTEGER, INTENT(IN), OPTIONAL :: nderivative
423 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
424 OPTIONAL, POINTER :: sab_orb
425 REAL(kind=dp), INTENT(IN), OPTIONAL :: eps_filter
426 CHARACTER(LEN=*), OPTIONAL :: basis_type
427 LOGICAL, INTENT(IN), OPTIONAL :: debug_forces, debug_stress
428
429 CHARACTER(LEN=default_string_length) :: my_basis_type
430 INTEGER :: ic, img, iounit, nimages
431 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
432 LOGICAL :: debfor, debstr, ext_kp, use_virial
433 REAL(kind=dp) :: fconv
434 REAL(kind=dp), DIMENSION(3) :: fodeb
435 REAL(kind=dp), DIMENSION(3, 3) :: stdeb
436 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
437 TYPE(cell_type), POINTER :: cell
438 TYPE(cp_logger_type), POINTER :: logger
439 TYPE(dft_control_type), POINTER :: dft_control
440 TYPE(kpoint_type), POINTER :: kpoints
441 TYPE(mp_para_env_type), POINTER :: para_env
442 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
443 POINTER :: sab_nl
444 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
445 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
446 TYPE(qs_ks_env_type), POINTER :: ks_env
447 TYPE(virial_type), POINTER :: virial
448
449 logger => cp_get_default_logger()
450 IF (logger%para_env%is_source()) THEN
451 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
452 ELSE
453 iounit = -1
454 END IF
455
456 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env)
457 ! is this a orbital-free method calculation
458 IF (dft_control%qs_control%ofgpw) RETURN
459
460 ext_kp = .false.
461 IF (PRESENT(ext_kpoints)) THEN
462 IF (ASSOCIATED(ext_kpoints)) ext_kp = .true.
463 END IF
464
465 ! basis type
466 IF (PRESENT(basis_type)) THEN
467 my_basis_type = basis_type
468 ELSE
469 my_basis_type = "ORB"
470 END IF
471
472 debfor = .false.
473 IF (PRESENT(debug_forces)) debfor = debug_forces
474 debfor = debfor .AND. calculate_forces
475
476 CALL get_qs_env(qs_env=qs_env, virial=virial)
477 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
478
479 debstr = .false.
480 IF (PRESENT(debug_stress)) debstr = debug_stress
481 debstr = debstr .AND. use_virial
482 IF (debstr) THEN
483 CALL get_qs_env(qs_env=qs_env, cell=cell)
484 fconv = 1.0e-9_dp*pascal/cell%deth
485 END IF
486
487 NULLIFY (ks_env)
488 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
489 NULLIFY (sab_nl)
490 IF (PRESENT(sab_orb)) THEN
491 sab_nl => sab_orb
492 ELSE
493 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
494 END IF
495
496 IF (calculate_forces) THEN
497 IF (SIZE(matrix_p, 1) == 2) THEN
498 DO img = 1, SIZE(matrix_p, 2)
499 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
500 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
501 END DO
502 END IF
503 END IF
504
505 IF (debfor) THEN
506 CALL get_qs_env(qs_env=qs_env, force=force)
507 fodeb(1:3) = force(1)%kinetic(1:3, 1)
508 END IF
509 IF (debstr) THEN
510 stdeb = virial%pv_ekinetic
511 END IF
512
513 ! T matrix
514 IF (ext_kp) THEN
515 CALL build_kinetic_matrix(ks_env, matrixkp_t=matrixkp_t, matrix_t=matrix_t, &
516 matrix_name=matrix_name, basis_type=my_basis_type, &
517 sab_nl=sab_nl, calculate_forces=calculate_forces, &
518 matrixkp_p=matrix_p, ext_kpoints=ext_kpoints, &
519 nderivative=nderivative, eps_filter=eps_filter)
520 ELSE
521 CALL build_kinetic_matrix(ks_env, matrixkp_t=matrixkp_t, matrix_t=matrix_t, &
522 matrix_name=matrix_name, basis_type=my_basis_type, &
523 sab_nl=sab_nl, calculate_forces=calculate_forces, &
524 matrixkp_p=matrix_p, nderivative=nderivative, &
525 eps_filter=eps_filter)
526 END IF
527
528 IF (debfor) THEN
529 fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
530 CALL para_env%sum(fodeb)
531 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dT ", fodeb
532 END IF
533 IF (debstr) THEN
534 stdeb = fconv*(virial%pv_ekinetic - stdeb)
535 CALL para_env%sum(stdeb)
536 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
537 'STRESS| P*dT', one_third_sum_diag(stdeb), det_3x3(stdeb)
538 END IF
539
540 IF (calculate_forces) THEN
541 IF (SIZE(matrix_p, 1) == 2) THEN
542 DO img = 1, SIZE(matrix_p, 2)
543 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
544 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
545 END DO
546 END IF
547 END IF
548
549 ! relativistic atomic correction to kinetic energy
550 IF (qs_env%rel_control%rel_method /= rel_none) THEN
551 IF (qs_env%rel_control%rel_transformation == rel_trans_atom) THEN
552 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
553 IF (ext_kp) THEN
554 ic = ext_kpoints%cell_to_index(0, 0, 0)
555 ELSE
556 nimages = dft_control%nimages
557 IF (nimages == 1) THEN
558 ic = 1
559 ELSE
560 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
561 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
562 ic = cell_to_index(0, 0, 0)
563 END IF
564 END IF
565 IF (my_basis_type /= "ORB") THEN
566 cpabort("Basis mismatch for relativistic correction")
567 END IF
568 IF (PRESENT(matrixkp_t)) THEN
569 CALL build_atomic_relmat(matrixkp_t(1, ic)%matrix, atomic_kind_set, qs_kind_set)
570 ELSEIF (PRESENT(matrix_t)) THEN
571 CALL build_atomic_relmat(matrix_t(1)%matrix, atomic_kind_set, qs_kind_set)
572 END IF
573 ELSE
574 cpabort("Relativistic corrections of this type are currently not implemented")
575 END IF
576 END IF
577
578 END SUBROUTINE kinetic_energy_matrix
579
580! **************************************************************************************************
581!> \brief Adds atomic blocks of relativistic correction for the kinetic energy
582!> \param matrix_t ...
583!> \param atomic_kind_set ...
584!> \param qs_kind_set ...
585! **************************************************************************************************
586 SUBROUTINE build_atomic_relmat(matrix_t, atomic_kind_set, qs_kind_set)
587 TYPE(dbcsr_type), POINTER :: matrix_t
588 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
589 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
590
591 INTEGER :: iatom, ikind, jatom
592 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
593 REAL(kind=dp), DIMENSION(:, :), POINTER :: reltmat, tblock
594 TYPE(dbcsr_iterator_type) :: iter
595
596 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
597
598 CALL dbcsr_iterator_start(iter, matrix_t)
599 DO WHILE (dbcsr_iterator_blocks_left(iter))
600 CALL dbcsr_iterator_next_block(iter, iatom, jatom, tblock)
601 IF (iatom == jatom) THEN
602 ikind = kind_of(iatom)
603 CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat)
604 IF (ASSOCIATED(reltmat)) tblock = tblock + reltmat
605 END IF
606 END DO
607 CALL dbcsr_iterator_stop(iter)
608
609 END SUBROUTINE build_atomic_relmat
610
611END MODULE qs_core_matrices
mathlib::det_3x3
Definition mathlib.F:70
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:622
input_constants::rel_none
integer, parameter, public rel_none
Definition input_constants.F:616
input_constants::do_ppl_analytic
integer, parameter, public do_ppl_analytic
Definition input_constants.F:280
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:372
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::core_matrices
subroutine, public core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, ec_env, dcdr_env, ec_env_matrices, ext_kpoints, basis_type, debug_forces, debug_stress, atcore)
...
Definition qs_core_matrices.F:144
qs_core_matrices::kinetic_energy_matrix
subroutine, public kinetic_energy_matrix(qs_env, matrixkp_t, matrix_t, matrix_p, ext_kpoints, 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:411
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, mimic, 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, xcint_weights, 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:530
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, monovalent, 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:468
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, ext_kpoints, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
Definition qs_kinetic.F:102
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, exc_accint, 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, xcint_weights, 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:341
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:60
cp_control_types::dft_control_type
Definition cp_control_types.F:633
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:63
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:157
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:743
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