Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of the core Hamiltonian integral matrix <a|H|b> over
10 : !> Cartesian Gaussian-type functions.
11 : !>
12 : !> Nuclear potential energy:
13 : !>
14 : !> a) Allelectron calculation:
15 : !>
16 : !> erfc(r)
17 : !> <a|V|b> = -Z*<a|---------|b>
18 : !> r
19 : !>
20 : !> 1 - erf(r)
21 : !> = -Z*<a|------------|b>
22 : !> r
23 : !>
24 : !> 1 erf(r)
25 : !> = -Z*(<a|---|b> - <a|--------|b>)
26 : !> r r
27 : !>
28 : !> 1
29 : !> = -Z*(<a|---|b> - N*<ab||c>)
30 : !> r
31 : !>
32 : !> -Z
33 : !> = <a|---|b> + Z*N*<ab||c>
34 : !> r
35 : !> \_______/ \_____/
36 : !> | |
37 : !> nuclear coulomb
38 : !>
39 : !> b) Pseudopotential calculation (Goedecker, Teter and Hutter; GTH):
40 : !>
41 : !> <a|V|b> = <a|(V(local) + V(non-local))|b>
42 : !>
43 : !> = <a|(V(local)|b> + <a|V(non-local))|b>
44 : !>
45 : !> <a|V(local)|b> = <a|-Z(eff)*erf(SQRT(2)*alpha*r)/r +
46 : !> (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
47 : !> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
48 : !>
49 : !> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b>
50 : !> \par Literature
51 : !> S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B 54, 1703 (1996)
52 : !> C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B 58, 3641 (1998)
53 : !> M. Krack and M. Parrinello, Phys. Chem. Chem. Phys. 2, 2105 (2000)
54 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
55 : !> \par History
56 : !> - Joost VandeVondele (April 2003) : added LSD forces
57 : !> - Non-redundant calculation of the non-local part of the GTH PP
58 : !> (22.05.2003,MK)
59 : !> - New parallelization scheme (27.06.2003,MK)
60 : !> - OpenMP version (07.12.2003,JGH)
61 : !> - Binary search loop for VPPNL operators (09.01.2004,JGH,MK)
62 : !> - Refactoring of pseudopotential and nuclear attraction integrals (25.02.2009,JGH)
63 : !> - General refactoring (01.10.2010,JGH)
64 : !> - Refactoring related to the new kinetic energy and overlap routines (07.2014,JGH)
65 : !> - k-point functionality (07.2015,JGH)
66 : !> - Refactor for PP functionality (08.2025,JGH)
67 : !> \author Matthias Krack (14.09.2000,21.03.02)
68 : ! **************************************************************************************************
69 : MODULE qs_core_matrices
70 : USE atomic_kind_types, ONLY: atomic_kind_type,&
71 : get_atomic_kind_set
72 : USE cell_types, ONLY: cell_type
73 : USE core_ae, ONLY: build_core_ae
74 : USE core_ppl, ONLY: build_core_ppl
75 : USE core_ppnl, ONLY: build_core_ppnl
76 : USE cp_control_types, ONLY: dft_control_type
77 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
78 : dbcsr_iterator_blocks_left,&
79 : dbcsr_iterator_next_block,&
80 : dbcsr_iterator_start,&
81 : dbcsr_iterator_stop,&
82 : dbcsr_iterator_type,&
83 : dbcsr_p_type,&
84 : dbcsr_type
85 : USE cp_log_handling, ONLY: cp_get_default_logger,&
86 : cp_logger_get_default_unit_nr,&
87 : cp_logger_type
88 : USE ec_env_types, ONLY: energy_correction_type
89 : USE input_constants, ONLY: do_ppl_analytic,&
90 : rel_none,&
91 : rel_trans_atom
92 : USE kinds, ONLY: default_string_length,&
93 : dp
94 : USE kpoint_types, ONLY: get_kpoint_info,&
95 : kpoint_type
96 : USE lri_environment_types, ONLY: lri_environment_type
97 : USE mathlib, ONLY: det_3x3
98 : USE message_passing, ONLY: mp_para_env_type
99 : USE particle_types, ONLY: particle_type
100 : USE physcon, ONLY: pascal
101 : USE qs_environment_types, ONLY: get_qs_env,&
102 : qs_environment_type
103 : USE qs_force_types, ONLY: qs_force_type
104 : USE qs_kind_types, ONLY: get_qs_kind,&
105 : qs_kind_type
106 : USE qs_kinetic, ONLY: build_kinetic_matrix
107 : USE qs_ks_types, ONLY: get_ks_env,&
108 : qs_ks_env_type
109 : USE qs_linres_types, ONLY: dcdr_env_type
110 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
111 : USE virial_methods, ONLY: one_third_sum_diag
112 : USE virial_types, ONLY: virial_type
113 : #include "./base/base_uses.f90"
114 :
115 : IMPLICIT NONE
116 :
117 : PRIVATE
118 :
119 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_matrices'
120 :
121 : PUBLIC :: core_matrices, kinetic_energy_matrix
122 :
123 : CONTAINS
124 :
125 : ! **************************************************************************************************
126 : !> \brief ...
127 : !> \param qs_env ...
128 : !> \param matrix_h ...
129 : !> \param matrix_p ...
130 : !> \param calculate_forces ...
131 : !> \param nder ...
132 : !> \param ec_env ...
133 : !> \param dcdr_env ...
134 : !> \param ec_env_matrices ...
135 : !> \param ext_kpoints ...
136 : !> \param basis_type ...
137 : !> \param debug_forces ...
138 : !> \param debug_stress ...
139 : !> \param atcore ...
140 : ! **************************************************************************************************
141 19685 : 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 66 : 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 19685 : 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 19685 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: deltaR
167 19685 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
168 : TYPE(cell_type), POINTER :: cell
169 : TYPE(cp_logger_type), POINTER :: logger
170 19685 : 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 19685 : POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
177 19685 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
178 19685 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
179 19685 : 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 39370 : logger => cp_get_default_logger()
184 19685 : IF (logger%para_env%is_source()) THEN
185 10113 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
186 : ELSE
187 : iounit = -1
188 : END IF
189 :
190 19685 : ext_kp = .FALSE.
191 19685 : IF (PRESENT(ext_kpoints)) THEN
192 0 : IF (ASSOCIATED(ext_kpoints)) ext_kp = .TRUE.
193 : END IF
194 :
195 19685 : NULLIFY (dft_control)
196 19685 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env)
197 :
198 : ! basis type
199 19685 : IF (PRESENT(basis_type)) THEN
200 816 : my_basis_type = basis_type
201 : ELSE
202 18869 : my_basis_type = "ORB"
203 : END IF
204 :
205 19685 : IF (PRESENT(dcdr_env) .AND. PRESENT(ec_env)) THEN
206 0 : CPABORT("core_matrices: conflicting options")
207 : END IF
208 :
209 : ! check size of atcore array
210 19685 : IF (PRESENT(atcore)) THEN
211 66 : CALL get_qs_env(qs_env=qs_env, natom=natom)
212 66 : 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 19685 : my_gt_nl = .FALSE.
217 19685 : IF (qs_env%run_rtp) THEN
218 1252 : CPASSERT(ASSOCIATED(dft_control%rtp_control))
219 1252 : IF (dft_control%rtp_control%velocity_gauge) THEN
220 54 : my_gt_nl = dft_control%rtp_control%nl_gauge_transform
221 : END IF
222 : END IF
223 :
224 : ! prepare for k-points
225 19685 : NULLIFY (cell_to_index)
226 19685 : IF (ext_kp) THEN
227 0 : nimages = SIZE(ext_kpoints%index_to_cell, 2)
228 0 : CALL get_kpoint_info(kpoint=ext_kpoints, cell_to_index=cell_to_index)
229 0 : dft_control%qs_control%do_ppl_method = do_ppl_analytic
230 : ELSE
231 19685 : nimages = dft_control%nimages
232 19685 : IF (nimages > 1) THEN
233 390 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
234 390 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
235 390 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
236 390 : dft_control%qs_control%do_ppl_method = do_ppl_analytic
237 : END IF
238 : END IF
239 :
240 : ! Possible dc/dr terms
241 19685 : IF (PRESENT(dcdr_env)) THEN
242 72 : deltaR => dcdr_env%delta_basis_function
243 72 : matrix_hloc(1:3, 1:1) => dcdr_env%matrix_hc(1:3)
244 : ELSE
245 19613 : NULLIFY (deltaR)
246 19613 : matrix_hloc => matrix_h
247 : END IF
248 19685 : matrix_ploc => matrix_p
249 :
250 : ! force analytic ppl calculation for GAPW methods
251 19685 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
252 3048 : dft_control%qs_control%do_ppl_method = do_ppl_analytic
253 : END IF
254 :
255 : ! force
256 19685 : NULLIFY (force)
257 19685 : IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force)
258 : ! virial
259 19685 : CALL get_qs_env(qs_env=qs_env, virial=virial)
260 19685 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
261 : ! lri/gpw
262 19685 : use_lrigpw = .FALSE.
263 :
264 : !debug
265 19685 : debfor = .FALSE.
266 19685 : IF (PRESENT(debug_forces)) debfor = debug_forces
267 1890 : debfor = debfor .AND. calculate_forces
268 19685 : debstr = .FALSE.
269 19685 : IF (PRESENT(debug_stress)) debstr = debug_stress
270 19685 : debstr = debstr .AND. use_virial
271 19685 : IF (debstr) THEN
272 50 : CALL get_qs_env(qs_env=qs_env, cell=cell)
273 50 : fconv = 1.0E-9_dp*pascal/cell%deth
274 : END IF
275 :
276 19685 : 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 19685 : particle_set=particle_set)
279 :
280 19685 : NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
281 19685 : IF (PRESENT(ec_env)) THEN
282 816 : sab_orb => ec_env%sab_orb
283 816 : sac_ae => ec_env%sac_ae
284 816 : sac_ppl => ec_env%sac_ppl
285 816 : 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 18869 : sap_ppnl=sap_ppnl)
292 : END IF
293 :
294 19685 : IF (PRESENT(ec_env) .AND. PRESENT(ec_env_matrices)) THEN
295 816 : IF (ec_env_matrices) THEN
296 346 : matrix_hloc => ec_env%matrix_h
297 346 : matrix_ploc => ec_env%matrix_p
298 : END IF
299 : END IF
300 :
301 : ! *** compute the nuclear attraction contribution to the core hamiltonian ***
302 19685 : all_present = ASSOCIATED(sac_ae)
303 19685 : IF (all_present) THEN
304 1232 : IF (PRESENT(dcdr_env)) THEN
305 0 : CPABORT("ECP/AE functionality for dcdr missing")
306 : END IF
307 1286 : IF (debfor) fodeb(1:3) = force(1)%all_potential(1:3, 1)
308 1232 : 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 2460 : nimages, cell_to_index, my_basis_type, atcore=atcore)
312 1232 : IF (debfor) THEN
313 72 : fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
314 18 : CALL para_env%sum(fodeb)
315 18 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHae ", fodeb
316 : END IF
317 1232 : IF (debstr) THEN
318 0 : stdeb = fconv*(virial%pv_ppl - stdeb)
319 0 : CALL para_env%sum(stdeb)
320 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
321 0 : '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 19685 : ppl_present = ASSOCIATED(sac_ppl)
327 19685 : IF (ppl_present) THEN
328 18523 : IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
329 18499 : CALL get_qs_env(qs_env, lri_env=lri_env)
330 18499 : IF (ASSOCIATED(lri_env)) THEN
331 90 : use_lrigpw = (dft_control%qs_control%lrigpw .AND. lri_env%ppl_ri)
332 : END IF
333 : IF (use_lrigpw) THEN
334 4 : IF (lri_env%exact_1c_terms) THEN
335 0 : CPABORT("not implemented")
336 : END IF
337 : ELSE
338 19881 : IF (debfor) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
339 19095 : 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 36928 : deltaR=deltaR, atcore=atcore)
345 18495 : IF (debfor) THEN
346 1848 : fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
347 462 : CALL para_env%sum(fodeb)
348 462 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppl ", fodeb
349 : END IF
350 18495 : IF (debstr) THEN
351 650 : stdeb = fconv*(virial%pv_ppl - stdeb)
352 50 : CALL para_env%sum(stdeb)
353 50 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
354 25 : '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 19685 : eps_ppnl = dft_control%qs_control%eps_ppnl
362 19685 : ppnl_present = ASSOCIATED(sap_ppnl)
363 19685 : IF (ppnl_present) THEN
364 15620 : IF (PRESENT(dcdr_env)) THEN
365 72 : matrix_hloc(1:3, 1:1) => dcdr_env%matrix_ppnl_1(1:3)
366 : END IF
367 15620 : IF (.NOT. my_gt_nl) THEN
368 16928 : IF (debfor) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
369 16184 : 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 31106 : deltaR=deltaR, atcore=atcore)
375 15584 : IF (debfor) THEN
376 1792 : fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
377 448 : CALL para_env%sum(fodeb)
378 448 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppnl ", fodeb
379 : END IF
380 15584 : IF (debstr) THEN
381 650 : stdeb = fconv*(virial%pv_ppnl - stdeb)
382 50 : CALL para_env%sum(stdeb)
383 50 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
384 25 : 'STRESS| P*dHppnl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
385 : END IF
386 : END IF
387 : END IF
388 :
389 19751 : 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 19225 : 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 19225 : 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 19225 : 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 19225 : POINTER :: sab_nl
444 19225 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
445 19225 : 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 38450 : logger => cp_get_default_logger()
450 19225 : IF (logger%para_env%is_source()) THEN
451 9883 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
452 : ELSE
453 : iounit = -1
454 : END IF
455 :
456 19225 : 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 19225 : IF (dft_control%qs_control%ofgpw) RETURN
459 :
460 19225 : ext_kp = .FALSE.
461 19225 : IF (PRESENT(ext_kpoints)) THEN
462 0 : IF (ASSOCIATED(ext_kpoints)) ext_kp = .TRUE.
463 : END IF
464 :
465 : ! basis type
466 19225 : IF (PRESENT(basis_type)) THEN
467 18937 : my_basis_type = basis_type
468 : ELSE
469 288 : my_basis_type = "ORB"
470 : END IF
471 :
472 19225 : debfor = .FALSE.
473 19225 : IF (PRESENT(debug_forces)) debfor = debug_forces
474 1890 : debfor = debfor .AND. calculate_forces
475 :
476 19225 : CALL get_qs_env(qs_env=qs_env, virial=virial)
477 19225 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
478 :
479 19225 : debstr = .FALSE.
480 19225 : IF (PRESENT(debug_stress)) debstr = debug_stress
481 21115 : debstr = debstr .AND. use_virial
482 1890 : IF (debstr) THEN
483 50 : CALL get_qs_env(qs_env=qs_env, cell=cell)
484 50 : fconv = 1.0E-9_dp*pascal/cell%deth
485 : END IF
486 :
487 19225 : NULLIFY (ks_env)
488 19225 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
489 19225 : NULLIFY (sab_nl)
490 19225 : IF (PRESENT(sab_orb)) THEN
491 18937 : sab_nl => sab_orb
492 : ELSE
493 288 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
494 : END IF
495 :
496 19225 : IF (calculate_forces) THEN
497 7973 : IF (SIZE(matrix_p, 1) == 2) THEN
498 2862 : DO img = 1, SIZE(matrix_p, 2)
499 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
500 2862 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
501 : END DO
502 : END IF
503 : END IF
504 :
505 19225 : IF (debfor) THEN
506 480 : CALL get_qs_env(qs_env=qs_env, force=force)
507 1920 : fodeb(1:3) = force(1)%kinetic(1:3, 1)
508 : END IF
509 19225 : IF (debstr) THEN
510 650 : stdeb = virial%pv_ekinetic
511 : END IF
512 :
513 : ! T matrix
514 19225 : 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 0 : 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 19983 : eps_filter=eps_filter)
526 : END IF
527 :
528 19225 : IF (debfor) THEN
529 1920 : fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
530 480 : CALL para_env%sum(fodeb)
531 480 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dT ", fodeb
532 : END IF
533 19225 : IF (debstr) THEN
534 650 : stdeb = fconv*(virial%pv_ekinetic - stdeb)
535 50 : CALL para_env%sum(stdeb)
536 50 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
537 25 : 'STRESS| P*dT', one_third_sum_diag(stdeb), det_3x3(stdeb)
538 : END IF
539 :
540 19225 : IF (calculate_forces) THEN
541 7973 : IF (SIZE(matrix_p, 1) == 2) THEN
542 2862 : DO img = 1, SIZE(matrix_p, 2)
543 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
544 2862 : 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 19225 : IF (qs_env%rel_control%rel_method /= rel_none) THEN
551 58 : IF (qs_env%rel_control%rel_transformation == rel_trans_atom) THEN
552 58 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
553 58 : IF (ext_kp) THEN
554 0 : ic = ext_kpoints%cell_to_index(0, 0, 0)
555 : ELSE
556 58 : nimages = dft_control%nimages
557 58 : IF (nimages == 1) THEN
558 : ic = 1
559 : ELSE
560 4 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
561 4 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
562 4 : ic = cell_to_index(0, 0, 0)
563 : END IF
564 : END IF
565 58 : IF (my_basis_type /= "ORB") THEN
566 0 : CPABORT("Basis mismatch for relativistic correction")
567 : END IF
568 58 : IF (PRESENT(matrixkp_t)) THEN
569 58 : CALL build_atomic_relmat(matrixkp_t(1, ic)%matrix, atomic_kind_set, qs_kind_set)
570 0 : ELSEIF (PRESENT(matrix_t)) THEN
571 0 : CALL build_atomic_relmat(matrix_t(1)%matrix, atomic_kind_set, qs_kind_set)
572 : END IF
573 : ELSE
574 0 : CPABORT("Relativistic corrections of this type are currently not implemented")
575 : END IF
576 : END IF
577 :
578 19225 : 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 58 : 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 58 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
593 58 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: reltmat, tblock
594 : TYPE(dbcsr_iterator_type) :: iter
595 :
596 58 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
597 :
598 58 : CALL dbcsr_iterator_start(iter, matrix_t)
599 221 : DO WHILE (dbcsr_iterator_blocks_left(iter))
600 163 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, tblock)
601 221 : IF (iatom == jatom) THEN
602 83 : ikind = kind_of(iatom)
603 83 : CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat)
604 192766 : IF (ASSOCIATED(reltmat)) tblock = tblock + reltmat
605 : END IF
606 : END DO
607 58 : CALL dbcsr_iterator_stop(iter)
608 :
609 116 : END SUBROUTINE build_atomic_relmat
610 :
611 : END MODULE qs_core_matrices
|
