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