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 : !> <a|H|b> = <a|T|b> + <a|V|b>
13 : !>
14 : !> Kinetic energy:
15 : !>
16 : !> <a|T|b> = <a|-nabla**2/2|b>
17 : !> \_______________/
18 : !> |
19 : !> kinetic
20 : !>
21 : !> Nuclear potential energy:
22 : !>
23 : !> a) Allelectron calculation:
24 : !>
25 : !> erfc(r)
26 : !> <a|V|b> = -Z*<a|---------|b>
27 : !> r
28 : !>
29 : !> 1 - erf(r)
30 : !> = -Z*<a|------------|b>
31 : !> r
32 : !>
33 : !> 1 erf(r)
34 : !> = -Z*(<a|---|b> - <a|--------|b>)
35 : !> r r
36 : !>
37 : !> 1
38 : !> = -Z*(<a|---|b> - N*<ab||c>)
39 : !> r
40 : !>
41 : !> -Z
42 : !> = <a|---|b> + Z*N*<ab||c>
43 : !> r
44 : !> \_______/ \_____/
45 : !> | |
46 : !> nuclear coulomb
47 : !>
48 : !> b) Pseudopotential calculation (Goedecker, Teter and Hutter; GTH):
49 : !>
50 : !> <a|V|b> = <a|(V(local) + V(non-local))|b>
51 : !>
52 : !> = <a|(V(local)|b> + <a|V(non-local))|b>
53 : !>
54 : !> <a|V(local)|b> = <a|-Z(eff)*erf(SQRT(2)*alpha*r)/r +
55 : !> (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
56 : !> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
57 : !>
58 : !> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b>
59 : !> \par Literature
60 : !> S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B 54, 1703 (1996)
61 : !> C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B 58, 3641 (1998)
62 : !> M. Krack and M. Parrinello, Phys. Chem. Chem. Phys. 2, 2105 (2000)
63 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
64 : !> \par History
65 : !> - Joost VandeVondele (April 2003) : added LSD forces
66 : !> - Non-redundant calculation of the non-local part of the GTH PP
67 : !> (22.05.2003,MK)
68 : !> - New parallelization scheme (27.06.2003,MK)
69 : !> - OpenMP version (07.12.2003,JGH)
70 : !> - Binary search loop for VPPNL operators (09.01.2004,JGH,MK)
71 : !> - Refactoring of pseudopotential and nuclear attraction integrals (25.02.2009,JGH)
72 : !> - General refactoring (01.10.2010,JGH)
73 : !> - Refactoring related to the new kinetic energy and overlap routines (07.2014,JGH)
74 : !> - k-point functionality (07.2015,JGH)
75 : !> \author Matthias Krack (14.09.2000,21.03.02)
76 : ! **************************************************************************************************
77 : MODULE qs_core_hamiltonian
78 : USE atomic_kind_types, ONLY: atomic_kind_type
79 : USE cp_blacs_env, ONLY: cp_blacs_env_type
80 : USE cp_control_types, ONLY: dft_control_type
81 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
82 : dbcsr_copy,&
83 : dbcsr_create,&
84 : dbcsr_distribution_type,&
85 : dbcsr_p_type,&
86 : dbcsr_set,&
87 : dbcsr_type,&
88 : dbcsr_type_antisymmetric
89 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
90 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
91 : dbcsr_deallocate_matrix_set
92 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_matrix_dist,&
93 : cp_dbcsr_write_sparse_matrix
94 : USE cp_log_handling, ONLY: cp_get_default_logger,&
95 : cp_logger_type
96 : USE cp_output_handling, ONLY: cp_p_file,&
97 : cp_print_key_finished_output,&
98 : cp_print_key_should_output,&
99 : cp_print_key_unit_nr
100 : USE input_constants, ONLY: do_admm_purify_none,&
101 : kg_tnadd_atomic,&
102 : plus_u_tensorial
103 : USE input_section_types, ONLY: section_vals_val_get
104 : USE kg_environment_types, ONLY: kg_environment_type
105 : USE kg_tnadd_mat, ONLY: build_tnadd_mat
106 : USE kinds, ONLY: default_string_length,&
107 : dp
108 : USE message_passing, ONLY: mp_para_env_type
109 : USE particle_types, ONLY: particle_type
110 : USE qs_cneo_methods, ONLY: cneo_core_matrices
111 : USE qs_condnum, ONLY: overlap_condnum
112 : USE qs_core_matrices, ONLY: core_matrices,&
113 : kinetic_energy_matrix
114 : USE qs_environment_types, ONLY: get_qs_env,&
115 : qs_environment_type,&
116 : set_qs_env
117 : USE qs_force_types, ONLY: qs_force_type
118 : USE qs_kind_types, ONLY: get_qs_kind,&
119 : qs_kind_type
120 : USE qs_ks_types, ONLY: get_ks_env,&
121 : qs_ks_env_type,&
122 : set_ks_env
123 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
124 : USE qs_oce_methods, ONLY: build_oce_matrices
125 : USE qs_oce_types, ONLY: allocate_oce_set,&
126 : create_oce_set,&
127 : oce_matrix_type
128 : USE qs_overlap, ONLY: build_overlap_matrix
129 : USE qs_rho_types, ONLY: qs_rho_get,&
130 : qs_rho_type
131 : USE virial_types, ONLY: virial_type
132 : #include "./base/base_uses.f90"
133 :
134 : IMPLICIT NONE
135 :
136 : PRIVATE
137 :
138 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_hamiltonian'
139 :
140 : PUBLIC :: build_core_hamiltonian_matrix
141 : PUBLIC :: dump_info_core_hamiltonian, qs_matrix_h_allocate_imag_from_real
142 :
143 : CONTAINS
144 :
145 : ! **************************************************************************************************
146 : !> \brief Cosntruction of the QS Core Hamiltonian Matrix
147 : !> \param qs_env ...
148 : !> \param calculate_forces ...
149 : !> \author Creation (11.03.2002,MK)
150 : !> Non-redundant calculation of the non-local part of the GTH PP (22.05.2003,MK)
151 : !> New parallelization scheme (27.06.2003,MK)
152 : ! **************************************************************************************************
153 18933 : SUBROUTINE build_core_hamiltonian_matrix(qs_env, calculate_forces)
154 :
155 : TYPE(qs_environment_type), POINTER :: qs_env
156 : LOGICAL, INTENT(IN) :: calculate_forces
157 :
158 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_hamiltonian_matrix'
159 :
160 : INTEGER :: handle, ikind, img, iw, nder, nders, &
161 : nimages, nkind
162 : LOGICAL :: dft_plus_u_atom, h_is_complex, norml1, &
163 : norml2, ofdft, &
164 : tensorial_plus_u_forces, use_arnoldi, &
165 : use_virial
166 : REAL(KIND=dp) :: eps_filter, eps_fit
167 : REAL(KIND=dp), DIMENSION(2) :: condnum
168 18933 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
169 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
170 : TYPE(cp_logger_type), POINTER :: logger
171 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
172 18933 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_t, &
173 18933 : matrix_w
174 : TYPE(dft_control_type), POINTER :: dft_control
175 : TYPE(kg_environment_type), POINTER :: kg_env
176 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
177 18933 : POINTER :: sab_orb, sap_oce
178 : TYPE(oce_matrix_type), POINTER :: oce
179 18933 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
180 18933 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
181 18933 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
182 : TYPE(qs_ks_env_type), POINTER :: ks_env
183 : TYPE(qs_rho_type), POINTER :: rho
184 : TYPE(virial_type), POINTER :: virial
185 :
186 37866 : IF (calculate_forces) THEN
187 6341 : CALL timeset(routineN//"_forces", handle)
188 : ELSE
189 12592 : CALL timeset(routineN, handle)
190 : END IF
191 :
192 18933 : NULLIFY (logger)
193 18933 : logger => cp_get_default_logger()
194 :
195 18933 : NULLIFY (dft_control)
196 18933 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
197 :
198 18933 : NULLIFY (qs_kind_set)
199 18933 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
200 :
201 : ! is this a orbital-free method calculation
202 18933 : ofdft = dft_control%qs_control%ofgpw
203 :
204 18933 : nimages = dft_control%nimages
205 18933 : IF (ofdft) THEN
206 0 : CPASSERT(nimages == 1)
207 : END IF
208 :
209 18933 : nders = 0
210 18933 : IF (calculate_forces) THEN
211 6341 : nder = 1
212 : ELSE
213 12592 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
214 : "DFT%PRINT%AO_MATRICES/DERIVATIVES") /= 0) THEN
215 4 : nder = 1
216 : ELSE
217 12588 : nder = 0
218 : END IF
219 : END IF
220 :
221 18933 : IF ((cp_print_key_should_output(logger%iter_info, qs_env%input, &
222 : "DFT%PRINT%AO_MATRICES/OVERLAP") /= 0 .AND. &
223 : BTEST(cp_print_key_should_output(logger%iter_info, qs_env%input, &
224 : "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file))) THEN
225 4 : nders = 1
226 : END IF
227 :
228 : ! the delta pulse in the periodic case needs the momentum operator,
229 : ! which is equivalent to the derivative of the overlap matrix
230 18933 : IF (ASSOCIATED(dft_control%rtp_control)) THEN
231 1854 : IF (dft_control%rtp_control%apply_delta_pulse .AND. &
232 : dft_control%rtp_control%periodic) THEN
233 180 : nders = 1
234 : END IF
235 : END IF
236 :
237 18933 : tensorial_plus_u_forces = .FALSE.
238 18933 : IF (dft_control%plus_u_method_id == plus_u_tensorial) THEN
239 32 : IF (ASSOCIATED(qs_kind_set)) THEN
240 32 : DO ikind = 1, SIZE(qs_kind_set)
241 : CALL get_qs_kind(qs_kind_set(ikind), &
242 32 : dft_plus_u_atom=dft_plus_u_atom)
243 32 : IF (dft_plus_u_atom) THEN
244 : tensorial_plus_u_forces = .TRUE.
245 : EXIT
246 : END IF
247 : END DO
248 : END IF
249 : END IF
250 32 : IF (tensorial_plus_u_forces) THEN
251 32 : nders = 1
252 : END IF
253 :
254 18933 : IF (dft_control%tddfpt2_control%enabled) THEN
255 1810 : nders = 1
256 1810 : IF (dft_control%do_admm) THEN
257 390 : IF (dft_control%admm_control%purification_method /= do_admm_purify_none) THEN
258 : CALL cp_abort(__LOCATION__, &
259 0 : "Only purification method NONE is possible with TDDFT at the moment")
260 : END IF
261 : END IF
262 : END IF
263 :
264 : ! filter for new matrices
265 18933 : eps_filter = dft_control%qs_control%eps_filter_matrix
266 : !
267 18933 : NULLIFY (ks_env)
268 18933 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
269 18933 : NULLIFY (matrix_s, matrix_t)
270 18933 : CALL get_qs_env(qs_env=qs_env, kinetic_kp=matrix_t, matrix_s_kp=matrix_s)
271 18933 : NULLIFY (sab_orb)
272 18933 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
273 18933 : NULLIFY (rho, force, matrix_p, matrix_w)
274 18933 : IF (calculate_forces) THEN
275 6341 : CALL get_qs_env(qs_env=qs_env, force=force, matrix_w_kp=matrix_w)
276 6341 : CALL get_qs_env(qs_env=qs_env, rho=rho)
277 6341 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
278 : ! *** If LSD, then combine alpha density and beta density to
279 : ! *** total density: alpha <- alpha + beta and
280 : ! *** spin density: beta <- alpha - beta
281 : ! (since all things can be computed based on the sum of these matrices anyway)
282 : ! (matrix_p is restored at the end of the run, matrix_w is left in its modified state
283 : ! (as it should not be needed afterwards)
284 6341 : IF (SIZE(matrix_p, 1) == 2) THEN
285 2486 : DO img = 1, nimages
286 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
287 1624 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
288 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
289 1624 : alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
290 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
291 2486 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
292 : END DO
293 : END IF
294 : ELSE
295 : NULLIFY (matrix_p, matrix_w)
296 : END IF
297 :
298 : ! S matrix
299 : CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, &
300 : matrix_name="OVERLAP MATRIX", &
301 : basis_type_a="ORB", &
302 : basis_type_b="ORB", &
303 : sab_nl=sab_orb, calculate_forces=calculate_forces, &
304 18933 : matrixkp_p=matrix_w)
305 :
306 18933 : IF (calculate_forces) THEN
307 : ! *** If LSD, then recover alpha density and beta density ***
308 : ! *** from the total density (1) and the spin density (2) ***
309 : ! *** The W matrix is neglected, since it will be destroyed ***
310 : ! *** in the calling force routine after leaving this routine ***
311 6341 : IF (SIZE(matrix_p, 1) == 2) THEN
312 2486 : DO img = 1, nimages
313 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
314 1624 : alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
315 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
316 2486 : alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
317 : END DO
318 : END IF
319 : END IF
320 :
321 : ! T matrix
322 : CALL kinetic_energy_matrix(qs_env, matrixkp_t=matrix_t, &
323 : matrix_p=matrix_p, &
324 : matrix_name="KINETIC ENERGY MATRIX", &
325 : basis_type="ORB", &
326 : sab_orb=sab_orb, &
327 : calculate_forces=calculate_forces, &
328 18933 : eps_filter=eps_filter)
329 :
330 : ! (Re-)allocate H matrix based on overlap matrix
331 18933 : CALL get_ks_env(ks_env, complex_ks=h_is_complex)
332 18933 : CALL qs_matrix_h_allocate(qs_env, matrix_s(1, 1)%matrix, is_complex=h_is_complex)
333 :
334 18933 : NULLIFY (matrix_h)
335 18933 : CALL get_qs_env(qs_env, matrix_h_kp=matrix_h)
336 :
337 18933 : IF (.NOT. ofdft) THEN
338 127210 : DO img = 1, nimages
339 : CALL dbcsr_copy(matrix_h(1, img)%matrix, matrix_t(1, img)%matrix, &
340 127210 : keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX")
341 : END DO
342 : END IF
343 :
344 18933 : NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
345 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
346 18933 : particle_set=particle_set)
347 :
348 : ! *** core and pseudopotentials
349 18933 : CALL core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder)
350 :
351 : ! *** CNEO nuclear V_core
352 18933 : CALL cneo_core_matrices(qs_env, calculate_forces, nder)
353 :
354 : ! *** GAPW one-center-expansion (oce) matrices
355 18933 : NULLIFY (sap_oce)
356 18933 : CALL get_qs_env(qs_env=qs_env, sap_oce=sap_oce)
357 18933 : NULLIFY (oce)
358 18933 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
359 3372 : CALL get_qs_env(qs_env=qs_env, oce=oce)
360 3372 : CALL create_oce_set(oce)
361 3372 : nkind = SIZE(atomic_kind_set)
362 3372 : CALL allocate_oce_set(oce, nkind)
363 3372 : eps_fit = dft_control%qs_control%gapw_control%eps_fit
364 3372 : IF (ASSOCIATED(sap_oce)) THEN
365 : CALL build_oce_matrices(oce%intac, calculate_forces, nder, qs_kind_set, particle_set, &
366 3234 : sap_oce, eps_fit)
367 : END IF
368 : END IF
369 :
370 : ! *** KG atomic potentials for nonadditive kinetic energy
371 18933 : IF (dft_control%qs_control%do_kg) THEN
372 230 : IF (qs_env%kg_env%tnadd_method == kg_tnadd_atomic) THEN
373 42 : CALL get_qs_env(qs_env=qs_env, kg_env=kg_env, virial=virial, dbcsr_dist=dbcsr_dist)
374 42 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
375 : CALL build_tnadd_mat(kg_env, matrix_p, force, virial, calculate_forces, use_virial, &
376 42 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, dbcsr_dist)
377 : END IF
378 : END IF
379 :
380 : ! *** Put the core Hamiltonian matrix in the QS environment ***
381 18933 : CALL set_qs_env(qs_env, oce=oce)
382 18933 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s, kinetic_kp=matrix_t, matrix_h_kp=matrix_h)
383 :
384 : ! *** Print matrices if requested
385 18933 : CALL dump_info_core_hamiltonian(qs_env, calculate_forces)
386 :
387 : ! *** Overlap condition number
388 18933 : IF (.NOT. calculate_forces) THEN
389 12592 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
390 : "DFT%PRINT%OVERLAP_CONDITION") /= 0) THEN
391 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
392 38 : extension=".Log")
393 38 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
394 38 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
395 38 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
396 38 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
397 38 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
398 : END IF
399 : END IF
400 :
401 18933 : CALL timestop(handle)
402 :
403 18933 : END SUBROUTINE build_core_hamiltonian_matrix
404 :
405 : ! **************************************************************************************************
406 : !> \brief Possibly prints matrices after the construction of the Core
407 : !> Hamiltonian Matrix
408 : !> \param qs_env ...
409 : !> \param calculate_forces ...
410 : ! **************************************************************************************************
411 37866 : SUBROUTINE dump_info_core_hamiltonian(qs_env, calculate_forces)
412 : TYPE(qs_environment_type), POINTER :: qs_env
413 : LOGICAL, INTENT(IN) :: calculate_forces
414 :
415 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_info_core_hamiltonian'
416 :
417 : INTEGER :: after, handle, i, ic, iw, output_unit
418 : LOGICAL :: omit_headers
419 : TYPE(cp_logger_type), POINTER :: logger
420 18933 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_v
421 18933 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_h, matrixkp_s, matrixkp_t
422 : TYPE(mp_para_env_type), POINTER :: para_env
423 :
424 18933 : CALL timeset(routineN, handle)
425 :
426 18933 : NULLIFY (logger, matrix_v, para_env)
427 18933 : logger => cp_get_default_logger()
428 18933 : CALL get_qs_env(qs_env, para_env=para_env)
429 :
430 : ! Print the distribution of the overlap matrix blocks
431 : ! this duplicates causes duplicate printing at the force calc
432 18933 : IF (.NOT. calculate_forces) THEN
433 12592 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
434 : qs_env%input, "PRINT%DISTRIBUTION"), cp_p_file)) THEN
435 : output_unit = cp_print_key_unit_nr(logger, qs_env%input, "PRINT%DISTRIBUTION", &
436 92 : extension=".distribution")
437 92 : CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
438 92 : CALL cp_dbcsr_write_matrix_dist(matrixkp_s(1, 1)%matrix, output_unit, para_env)
439 92 : CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, "PRINT%DISTRIBUTION")
440 : END IF
441 : END IF
442 :
443 18933 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
444 : ! Print the overlap integral matrix, if requested
445 18933 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
446 : qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
447 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
448 6 : extension=".Log")
449 6 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
450 6 : after = MIN(MAX(after, 1), 16)
451 6 : CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
452 6 : IF (ASSOCIATED(matrixkp_s)) THEN
453 12 : DO ic = 1, SIZE(matrixkp_s, 2)
454 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(1, ic)%matrix, 4, after, qs_env, para_env, &
455 12 : output_unit=iw, omit_headers=omit_headers)
456 : END DO
457 6 : IF (BTEST(cp_print_key_should_output(logger%iter_info, qs_env%input, &
458 : "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
459 8 : DO ic = 1, SIZE(matrixkp_s, 2)
460 20 : DO i = 2, SIZE(matrixkp_s, 1)
461 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(i, ic)%matrix, 4, after, qs_env, para_env, &
462 16 : output_unit=iw, omit_headers=omit_headers)
463 : END DO
464 : END DO
465 : END IF
466 : END IF
467 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
468 6 : "DFT%PRINT%AO_MATRICES/OVERLAP")
469 : END IF
470 :
471 : ! Print the kinetic energy integral matrix, if requested
472 18933 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
473 : qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY"), cp_p_file)) THEN
474 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY", &
475 50 : extension=".Log")
476 50 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
477 50 : after = MIN(MAX(after, 1), 16)
478 50 : CALL get_qs_env(qs_env, kinetic_kp=matrixkp_t)
479 50 : IF (ASSOCIATED(matrixkp_t)) THEN
480 100 : DO ic = 1, SIZE(matrixkp_t, 2)
481 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_t(1, ic)%matrix, 4, after, qs_env, para_env, &
482 100 : output_unit=iw, omit_headers=omit_headers)
483 : END DO
484 : END IF
485 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
486 50 : "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY")
487 : END IF
488 :
489 : ! Print the potential energy matrix, if requested
490 18933 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
491 : qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY"), cp_p_file)) THEN
492 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY", &
493 50 : extension=".Log")
494 50 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
495 50 : after = MIN(MAX(after, 1), 16)
496 50 : CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h, kinetic_kp=matrixkp_t)
497 50 : IF (ASSOCIATED(matrixkp_h)) THEN
498 50 : IF (SIZE(matrixkp_h, 2) == 1) THEN
499 50 : CALL dbcsr_allocate_matrix_set(matrix_v, 1)
500 50 : ALLOCATE (matrix_v(1)%matrix)
501 50 : CALL dbcsr_copy(matrix_v(1)%matrix, matrixkp_h(1, 1)%matrix, name="POTENTIAL ENERGY MATRIX")
502 : CALL dbcsr_add(matrix_v(1)%matrix, matrixkp_t(1, 1)%matrix, &
503 50 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
504 : CALL cp_dbcsr_write_sparse_matrix(matrix_v(1)%matrix, 4, after, qs_env, &
505 50 : para_env, output_unit=iw, omit_headers=omit_headers)
506 50 : CALL dbcsr_deallocate_matrix_set(matrix_v)
507 : ELSE
508 0 : CPWARN("Printing of potential energy matrix not implemented for k-points")
509 : END IF
510 : END IF
511 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
512 50 : "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY")
513 : END IF
514 :
515 : ! Print the core Hamiltonian matrix, if requested
516 18933 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
517 : qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
518 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
519 50 : extension=".Log")
520 50 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
521 50 : after = MIN(MAX(after, 1), 16)
522 50 : CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h)
523 50 : IF (ASSOCIATED(matrixkp_h)) THEN
524 100 : DO ic = 1, SIZE(matrixkp_h, 2)
525 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_h(1, ic)%matrix, 4, after, qs_env, para_env, &
526 100 : output_unit=iw, omit_headers=omit_headers)
527 : END DO
528 : END IF
529 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
530 50 : "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
531 : END IF
532 :
533 18933 : CALL timestop(handle)
534 :
535 18933 : END SUBROUTINE dump_info_core_hamiltonian
536 :
537 : ! **************************************************************************************************
538 : !> \brief (Re-)allocate matrix_h based on the template (typically the overlap matrix)
539 : !> \param qs_env ...
540 : !> \param template ...
541 : !> \param is_complex ...
542 : ! **************************************************************************************************
543 18933 : SUBROUTINE qs_matrix_h_allocate(qs_env, template, is_complex)
544 : TYPE(qs_environment_type) :: qs_env
545 : TYPE(dbcsr_type), INTENT(in) :: template
546 : LOGICAL, INTENT(in) :: is_complex
547 :
548 : CHARACTER(LEN=default_string_length) :: headline
549 : INTEGER :: img, nimages
550 18933 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im
551 : TYPE(dft_control_type), POINTER :: dft_control
552 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
553 18933 : POINTER :: sab_orb
554 : TYPE(qs_ks_env_type), POINTER :: ks_env
555 :
556 18933 : NULLIFY (matrix_h, matrix_h_im, sab_orb, dft_control, ks_env)
557 : CALL get_qs_env(qs_env=qs_env, &
558 : matrix_h_kp=matrix_h, &
559 : matrix_h_im_kp=matrix_h_im, &
560 : sab_orb=sab_orb, &
561 : dft_control=dft_control, &
562 18933 : ks_env=ks_env)
563 :
564 18933 : nimages = dft_control%nimages
565 18933 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimages)
566 18933 : headline = "CORE HAMILTONIAN MATRIX"
567 127210 : DO img = 1, nimages
568 108277 : ALLOCATE (matrix_h(1, img)%matrix)
569 108277 : CALL dbcsr_create(matrix_h(1, img)%matrix, name=TRIM(headline), template=template)
570 108277 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
571 127210 : CALL dbcsr_set(matrix_h(1, img)%matrix, 0.0_dp)
572 : END DO
573 18933 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
574 :
575 18933 : IF (is_complex) THEN
576 366 : headline = "IMAGINARY PART OF CORE HAMILTONIAN MATRIX"
577 366 : CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
578 732 : DO img = 1, nimages
579 366 : ALLOCATE (matrix_h_im(1, img)%matrix)
580 : CALL dbcsr_create(matrix_h_im(1, img)%matrix, name=TRIM(headline), template=template, &
581 366 : matrix_type=dbcsr_type_antisymmetric)
582 366 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, img)%matrix, sab_orb)
583 732 : CALL dbcsr_set(matrix_h_im(1, img)%matrix, 0.0_dp)
584 : END DO
585 366 : CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
586 : END IF
587 :
588 18933 : END SUBROUTINE qs_matrix_h_allocate
589 :
590 : ! **************************************************************************************************
591 : !> \brief (Re-)allocates matrix_h_im from matrix_h
592 : !> \param qs_env ...
593 : ! **************************************************************************************************
594 8 : SUBROUTINE qs_matrix_h_allocate_imag_from_real(qs_env)
595 : TYPE(qs_environment_type) :: qs_env
596 :
597 : CHARACTER(LEN=default_string_length) :: headline
598 : INTEGER :: image, nimages
599 8 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im
600 : TYPE(dbcsr_type), POINTER :: template
601 : TYPE(dft_control_type), POINTER :: dft_control
602 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
603 8 : POINTER :: sab_orb
604 : TYPE(qs_ks_env_type), POINTER :: ks_env
605 :
606 8 : NULLIFY (matrix_h_im, matrix_h, dft_control, template, sab_orb, ks_env)
607 :
608 : CALL get_qs_env(qs_env, &
609 : matrix_h_im_kp=matrix_h_im, &
610 : matrix_h_kp=matrix_h, &
611 : dft_control=dft_control, &
612 : sab_orb=sab_orb, &
613 8 : ks_env=ks_env)
614 :
615 8 : nimages = dft_control%nimages
616 :
617 8 : CPASSERT(nimages == SIZE(matrix_h, 2))
618 :
619 8 : CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
620 :
621 16 : DO image = 1, nimages
622 8 : headline = "IMAGINARY CORE HAMILTONIAN MATRIX"
623 8 : ALLOCATE (matrix_h_im(1, image)%matrix)
624 8 : template => matrix_h(1, image)%matrix ! base on real part, but anti-symmetric
625 : CALL dbcsr_create(matrix=matrix_h_im(1, image)%matrix, template=template, &
626 8 : name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric)
627 8 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, image)%matrix, sab_orb)
628 16 : CALL dbcsr_set(matrix_h_im(1, image)%matrix, 0.0_dp)
629 : END DO
630 8 : CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
631 :
632 8 : END SUBROUTINE qs_matrix_h_allocate_imag_from_real
633 :
634 : END MODULE qs_core_hamiltonian
|
