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 : USE input_section_types, ONLY: section_vals_val_get
103 : USE kg_environment_types, ONLY: kg_environment_type
104 : USE kg_tnadd_mat, ONLY: build_tnadd_mat
105 : USE kinds, ONLY: default_string_length,&
106 : dp
107 : USE message_passing, ONLY: mp_para_env_type
108 : USE particle_types, ONLY: particle_type
109 : USE qs_cneo_methods, ONLY: cneo_core_matrices
110 : USE qs_condnum, ONLY: overlap_condnum
111 : USE qs_core_matrices, ONLY: core_matrices,&
112 : kinetic_energy_matrix
113 : USE qs_environment_types, ONLY: get_qs_env,&
114 : qs_environment_type,&
115 : set_qs_env
116 : USE qs_force_types, ONLY: qs_force_type
117 : USE qs_kind_types, ONLY: qs_kind_type
118 : USE qs_ks_types, ONLY: get_ks_env,&
119 : qs_ks_env_type,&
120 : set_ks_env
121 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
122 : USE qs_oce_methods, ONLY: build_oce_matrices
123 : USE qs_oce_types, ONLY: allocate_oce_set,&
124 : create_oce_set,&
125 : oce_matrix_type
126 : USE qs_overlap, ONLY: build_overlap_matrix
127 : USE qs_rho_types, ONLY: qs_rho_get,&
128 : qs_rho_type
129 : USE virial_types, ONLY: virial_type
130 : #include "./base/base_uses.f90"
131 :
132 : IMPLICIT NONE
133 :
134 : PRIVATE
135 :
136 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_hamiltonian'
137 :
138 : PUBLIC :: build_core_hamiltonian_matrix
139 : PUBLIC :: dump_info_core_hamiltonian, qs_matrix_h_allocate_imag_from_real
140 :
141 : CONTAINS
142 :
143 : ! **************************************************************************************************
144 : !> \brief Cosntruction of the QS Core Hamiltonian Matrix
145 : !> \param qs_env ...
146 : !> \param calculate_forces ...
147 : !> \author Creation (11.03.2002,MK)
148 : !> Non-redundant calculation of the non-local part of the GTH PP (22.05.2003,MK)
149 : !> New parallelization scheme (27.06.2003,MK)
150 : ! **************************************************************************************************
151 18707 : SUBROUTINE build_core_hamiltonian_matrix(qs_env, calculate_forces)
152 :
153 : TYPE(qs_environment_type), POINTER :: qs_env
154 : LOGICAL, INTENT(IN) :: calculate_forces
155 :
156 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_hamiltonian_matrix'
157 :
158 : INTEGER :: handle, img, iw, nder, nders, nimages, &
159 : nkind
160 : LOGICAL :: h_is_complex, norml1, norml2, ofdft, &
161 : use_arnoldi, use_virial
162 : REAL(KIND=dp) :: eps_filter, eps_fit
163 : REAL(KIND=dp), DIMENSION(2) :: condnum
164 18707 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
165 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
166 : TYPE(cp_logger_type), POINTER :: logger
167 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
168 18707 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_t, &
169 18707 : matrix_w
170 : TYPE(dft_control_type), POINTER :: dft_control
171 : TYPE(kg_environment_type), POINTER :: kg_env
172 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
173 18707 : POINTER :: sab_orb, sap_oce
174 : TYPE(oce_matrix_type), POINTER :: oce
175 18707 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
176 18707 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
177 18707 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
178 : TYPE(qs_ks_env_type), POINTER :: ks_env
179 : TYPE(qs_rho_type), POINTER :: rho
180 : TYPE(virial_type), POINTER :: virial
181 :
182 37414 : IF (calculate_forces) THEN
183 6309 : CALL timeset(routineN//"_forces", handle)
184 : ELSE
185 12398 : CALL timeset(routineN, handle)
186 : END IF
187 :
188 18707 : NULLIFY (logger)
189 18707 : logger => cp_get_default_logger()
190 :
191 18707 : NULLIFY (dft_control)
192 18707 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
193 :
194 : ! is this a orbital-free method calculation
195 18707 : ofdft = dft_control%qs_control%ofgpw
196 :
197 18707 : nimages = dft_control%nimages
198 18707 : IF (ofdft) THEN
199 0 : CPASSERT(nimages == 1)
200 : END IF
201 :
202 18707 : nders = 0
203 18707 : IF (calculate_forces) THEN
204 6309 : nder = 1
205 : ELSE
206 12398 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
207 : "DFT%PRINT%AO_MATRICES/DERIVATIVES") /= 0) THEN
208 4 : nder = 1
209 : ELSE
210 12394 : nder = 0
211 : END IF
212 : END IF
213 :
214 18707 : IF ((cp_print_key_should_output(logger%iter_info, qs_env%input, &
215 : "DFT%PRINT%AO_MATRICES/OVERLAP") /= 0 .AND. &
216 : BTEST(cp_print_key_should_output(logger%iter_info, qs_env%input, &
217 : "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file))) THEN
218 4 : nders = 1
219 : END IF
220 :
221 : ! the delta pulse in the periodic case needs the momentum operator,
222 : ! which is equivalent to the derivative of the overlap matrix
223 18707 : IF (ASSOCIATED(dft_control%rtp_control)) THEN
224 1802 : IF (dft_control%rtp_control%apply_delta_pulse .AND. &
225 : dft_control%rtp_control%periodic) THEN
226 128 : nders = 1
227 : END IF
228 : END IF
229 :
230 18707 : IF (dft_control%tddfpt2_control%enabled) THEN
231 1810 : nders = 1
232 1810 : IF (dft_control%do_admm) THEN
233 390 : IF (dft_control%admm_control%purification_method /= do_admm_purify_none) THEN
234 : CALL cp_abort(__LOCATION__, &
235 0 : "Only purification method NONE is possible with TDDFT at the moment")
236 : END IF
237 : END IF
238 : END IF
239 :
240 : ! filter for new matrices
241 18707 : eps_filter = dft_control%qs_control%eps_filter_matrix
242 : !
243 18707 : NULLIFY (ks_env)
244 18707 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
245 18707 : NULLIFY (matrix_s, matrix_t)
246 18707 : CALL get_qs_env(qs_env=qs_env, kinetic_kp=matrix_t, matrix_s_kp=matrix_s)
247 18707 : NULLIFY (sab_orb)
248 18707 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
249 18707 : NULLIFY (rho, force, matrix_p, matrix_w)
250 18707 : IF (calculate_forces) THEN
251 6309 : CALL get_qs_env(qs_env=qs_env, force=force, matrix_w_kp=matrix_w)
252 6309 : CALL get_qs_env(qs_env=qs_env, rho=rho)
253 6309 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
254 : ! *** If LSD, then combine alpha density and beta density to
255 : ! *** total density: alpha <- alpha + beta and
256 : ! *** spin density: beta <- alpha - beta
257 : ! (since all things can be computed based on the sum of these matrices anyway)
258 : ! (matrix_p is restored at the end of the run, matrix_w is left in its modified state
259 : ! (as it should not be needed afterwards)
260 6309 : IF (SIZE(matrix_p, 1) == 2) THEN
261 2470 : DO img = 1, nimages
262 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
263 1616 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
264 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
265 1616 : alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
266 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
267 2470 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
268 : END DO
269 : END IF
270 : ELSE
271 : NULLIFY (matrix_p, matrix_w)
272 : END IF
273 :
274 : ! S matrix
275 : CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, &
276 : matrix_name="OVERLAP MATRIX", &
277 : basis_type_a="ORB", &
278 : basis_type_b="ORB", &
279 : sab_nl=sab_orb, calculate_forces=calculate_forces, &
280 18707 : matrixkp_p=matrix_w)
281 :
282 18707 : IF (calculate_forces) THEN
283 : ! *** If LSD, then recover alpha density and beta density ***
284 : ! *** from the total density (1) and the spin density (2) ***
285 : ! *** The W matrix is neglected, since it will be destroyed ***
286 : ! *** in the calling force routine after leaving this routine ***
287 6309 : IF (SIZE(matrix_p, 1) == 2) THEN
288 2470 : DO img = 1, nimages
289 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
290 1616 : alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
291 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
292 2470 : alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
293 : END DO
294 : END IF
295 : END IF
296 :
297 : ! T matrix
298 : CALL kinetic_energy_matrix(qs_env, matrixkp_t=matrix_t, &
299 : matrix_p=matrix_p, &
300 : matrix_name="KINETIC ENERGY MATRIX", &
301 : basis_type="ORB", &
302 : sab_orb=sab_orb, &
303 : calculate_forces=calculate_forces, &
304 18707 : eps_filter=eps_filter)
305 :
306 : ! (Re-)allocate H matrix based on overlap matrix
307 18707 : CALL get_ks_env(ks_env, complex_ks=h_is_complex)
308 18707 : CALL qs_matrix_h_allocate(qs_env, matrix_s(1, 1)%matrix, is_complex=h_is_complex)
309 :
310 18707 : NULLIFY (matrix_h)
311 18707 : CALL get_qs_env(qs_env, matrix_h_kp=matrix_h)
312 :
313 18707 : IF (.NOT. ofdft) THEN
314 126448 : DO img = 1, nimages
315 : CALL dbcsr_copy(matrix_h(1, img)%matrix, matrix_t(1, img)%matrix, &
316 126448 : keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX")
317 : END DO
318 : END IF
319 :
320 18707 : NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
321 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
322 18707 : particle_set=particle_set)
323 :
324 : ! *** core and pseudopotentials
325 18707 : CALL core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder)
326 :
327 : ! *** CNEO nuclear V_core
328 18707 : CALL cneo_core_matrices(qs_env, calculate_forces, nder)
329 :
330 : ! *** GAPW one-center-expansion (oce) matrices
331 18707 : NULLIFY (sap_oce)
332 18707 : CALL get_qs_env(qs_env=qs_env, sap_oce=sap_oce)
333 18707 : NULLIFY (oce)
334 18707 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
335 3262 : CALL get_qs_env(qs_env=qs_env, oce=oce)
336 3262 : CALL create_oce_set(oce)
337 3262 : nkind = SIZE(atomic_kind_set)
338 3262 : CALL allocate_oce_set(oce, nkind)
339 3262 : eps_fit = dft_control%qs_control%gapw_control%eps_fit
340 3262 : IF (ASSOCIATED(sap_oce)) THEN
341 : CALL build_oce_matrices(oce%intac, calculate_forces, nder, qs_kind_set, particle_set, &
342 3202 : sap_oce, eps_fit)
343 : END IF
344 : END IF
345 :
346 : ! *** KG atomic potentials for nonadditive kinetic energy
347 18707 : IF (dft_control%qs_control%do_kg) THEN
348 230 : IF (qs_env%kg_env%tnadd_method == kg_tnadd_atomic) THEN
349 42 : CALL get_qs_env(qs_env=qs_env, kg_env=kg_env, virial=virial, dbcsr_dist=dbcsr_dist)
350 42 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
351 : CALL build_tnadd_mat(kg_env, matrix_p, force, virial, calculate_forces, use_virial, &
352 42 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, dbcsr_dist)
353 : END IF
354 : END IF
355 :
356 : ! *** Put the core Hamiltonian matrix in the QS environment ***
357 18707 : CALL set_qs_env(qs_env, oce=oce)
358 18707 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s, kinetic_kp=matrix_t, matrix_h_kp=matrix_h)
359 :
360 : ! *** Print matrices if requested
361 18707 : CALL dump_info_core_hamiltonian(qs_env, calculate_forces)
362 :
363 : ! *** Overlap condition number
364 18707 : IF (.NOT. calculate_forces) THEN
365 12398 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
366 : "DFT%PRINT%OVERLAP_CONDITION") /= 0) THEN
367 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
368 38 : extension=".Log")
369 38 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
370 38 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
371 38 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
372 38 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
373 38 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
374 : END IF
375 : END IF
376 :
377 18707 : CALL timestop(handle)
378 :
379 18707 : END SUBROUTINE build_core_hamiltonian_matrix
380 :
381 : ! **************************************************************************************************
382 : !> \brief Possibly prints matrices after the construction of the Core
383 : !> Hamiltonian Matrix
384 : !> \param qs_env ...
385 : !> \param calculate_forces ...
386 : ! **************************************************************************************************
387 37414 : SUBROUTINE dump_info_core_hamiltonian(qs_env, calculate_forces)
388 : TYPE(qs_environment_type), POINTER :: qs_env
389 : LOGICAL, INTENT(IN) :: calculate_forces
390 :
391 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_info_core_hamiltonian'
392 :
393 : INTEGER :: after, handle, i, ic, iw, output_unit
394 : LOGICAL :: omit_headers
395 : TYPE(cp_logger_type), POINTER :: logger
396 18707 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_v
397 18707 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_h, matrixkp_s, matrixkp_t
398 : TYPE(mp_para_env_type), POINTER :: para_env
399 :
400 18707 : CALL timeset(routineN, handle)
401 :
402 18707 : NULLIFY (logger, matrix_v, para_env)
403 18707 : logger => cp_get_default_logger()
404 18707 : CALL get_qs_env(qs_env, para_env=para_env)
405 :
406 : ! Print the distribution of the overlap matrix blocks
407 : ! this duplicates causes duplicate printing at the force calc
408 18707 : IF (.NOT. calculate_forces) THEN
409 12398 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
410 : qs_env%input, "PRINT%DISTRIBUTION"), cp_p_file)) THEN
411 : output_unit = cp_print_key_unit_nr(logger, qs_env%input, "PRINT%DISTRIBUTION", &
412 92 : extension=".distribution")
413 92 : CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
414 92 : CALL cp_dbcsr_write_matrix_dist(matrixkp_s(1, 1)%matrix, output_unit, para_env)
415 92 : CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, "PRINT%DISTRIBUTION")
416 : END IF
417 : END IF
418 :
419 18707 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
420 : ! Print the overlap integral matrix, if requested
421 18707 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
422 : qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
423 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
424 6 : extension=".Log")
425 6 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
426 6 : after = MIN(MAX(after, 1), 16)
427 6 : CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
428 6 : IF (ASSOCIATED(matrixkp_s)) THEN
429 12 : DO ic = 1, SIZE(matrixkp_s, 2)
430 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(1, ic)%matrix, 4, after, qs_env, para_env, &
431 12 : output_unit=iw, omit_headers=omit_headers)
432 : END DO
433 6 : IF (BTEST(cp_print_key_should_output(logger%iter_info, qs_env%input, &
434 : "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
435 8 : DO ic = 1, SIZE(matrixkp_s, 2)
436 20 : DO i = 2, SIZE(matrixkp_s, 1)
437 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(i, ic)%matrix, 4, after, qs_env, para_env, &
438 16 : output_unit=iw, omit_headers=omit_headers)
439 : END DO
440 : END DO
441 : END IF
442 : END IF
443 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
444 6 : "DFT%PRINT%AO_MATRICES/OVERLAP")
445 : END IF
446 :
447 : ! Print the kinetic energy integral matrix, if requested
448 18707 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
449 : qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY"), cp_p_file)) THEN
450 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY", &
451 50 : extension=".Log")
452 50 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
453 50 : after = MIN(MAX(after, 1), 16)
454 50 : CALL get_qs_env(qs_env, kinetic_kp=matrixkp_t)
455 50 : IF (ASSOCIATED(matrixkp_t)) THEN
456 100 : DO ic = 1, SIZE(matrixkp_t, 2)
457 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_t(1, ic)%matrix, 4, after, qs_env, para_env, &
458 100 : output_unit=iw, omit_headers=omit_headers)
459 : END DO
460 : END IF
461 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
462 50 : "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY")
463 : END IF
464 :
465 : ! Print the potential energy matrix, if requested
466 18707 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
467 : qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY"), cp_p_file)) THEN
468 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY", &
469 50 : extension=".Log")
470 50 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
471 50 : after = MIN(MAX(after, 1), 16)
472 50 : CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h, kinetic_kp=matrixkp_t)
473 50 : IF (ASSOCIATED(matrixkp_h)) THEN
474 50 : IF (SIZE(matrixkp_h, 2) == 1) THEN
475 50 : CALL dbcsr_allocate_matrix_set(matrix_v, 1)
476 50 : ALLOCATE (matrix_v(1)%matrix)
477 50 : CALL dbcsr_copy(matrix_v(1)%matrix, matrixkp_h(1, 1)%matrix, name="POTENTIAL ENERGY MATRIX")
478 : CALL dbcsr_add(matrix_v(1)%matrix, matrixkp_t(1, 1)%matrix, &
479 50 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
480 : CALL cp_dbcsr_write_sparse_matrix(matrix_v(1)%matrix, 4, after, qs_env, &
481 50 : para_env, output_unit=iw, omit_headers=omit_headers)
482 50 : CALL dbcsr_deallocate_matrix_set(matrix_v)
483 : ELSE
484 0 : CPWARN("Printing of potential energy matrix not implemented for k-points")
485 : END IF
486 : END IF
487 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
488 50 : "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY")
489 : END IF
490 :
491 : ! Print the core Hamiltonian matrix, if requested
492 18707 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
493 : qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
494 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
495 50 : extension=".Log")
496 50 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
497 50 : after = MIN(MAX(after, 1), 16)
498 50 : CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h)
499 50 : IF (ASSOCIATED(matrixkp_h)) THEN
500 100 : DO ic = 1, SIZE(matrixkp_h, 2)
501 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_h(1, ic)%matrix, 4, after, qs_env, para_env, &
502 100 : output_unit=iw, omit_headers=omit_headers)
503 : END DO
504 : END IF
505 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
506 50 : "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
507 : END IF
508 :
509 18707 : CALL timestop(handle)
510 :
511 18707 : END SUBROUTINE dump_info_core_hamiltonian
512 :
513 : ! **************************************************************************************************
514 : !> \brief (Re-)allocate matrix_h based on the template (typically the overlap matrix)
515 : !> \param qs_env ...
516 : !> \param template ...
517 : !> \param is_complex ...
518 : ! **************************************************************************************************
519 18707 : SUBROUTINE qs_matrix_h_allocate(qs_env, template, is_complex)
520 : TYPE(qs_environment_type) :: qs_env
521 : TYPE(dbcsr_type), INTENT(in) :: template
522 : LOGICAL, INTENT(in) :: is_complex
523 :
524 : CHARACTER(LEN=default_string_length) :: headline
525 : INTEGER :: img, nimages
526 18707 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im
527 : TYPE(dft_control_type), POINTER :: dft_control
528 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
529 18707 : POINTER :: sab_orb
530 : TYPE(qs_ks_env_type), POINTER :: ks_env
531 :
532 18707 : NULLIFY (matrix_h, matrix_h_im, sab_orb, dft_control, ks_env)
533 : CALL get_qs_env(qs_env=qs_env, &
534 : matrix_h_kp=matrix_h, &
535 : matrix_h_im_kp=matrix_h_im, &
536 : sab_orb=sab_orb, &
537 : dft_control=dft_control, &
538 18707 : ks_env=ks_env)
539 :
540 18707 : nimages = dft_control%nimages
541 18707 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimages)
542 18707 : headline = "CORE HAMILTONIAN MATRIX"
543 126448 : DO img = 1, nimages
544 107741 : ALLOCATE (matrix_h(1, img)%matrix)
545 107741 : CALL dbcsr_create(matrix_h(1, img)%matrix, name=TRIM(headline), template=template)
546 107741 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
547 126448 : CALL dbcsr_set(matrix_h(1, img)%matrix, 0.0_dp)
548 : END DO
549 18707 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
550 :
551 18707 : IF (is_complex) THEN
552 366 : headline = "IMAGINARY PART OF CORE HAMILTONIAN MATRIX"
553 366 : CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
554 732 : DO img = 1, nimages
555 366 : ALLOCATE (matrix_h_im(1, img)%matrix)
556 : CALL dbcsr_create(matrix_h_im(1, img)%matrix, name=TRIM(headline), template=template, &
557 366 : matrix_type=dbcsr_type_antisymmetric)
558 366 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, img)%matrix, sab_orb)
559 732 : CALL dbcsr_set(matrix_h_im(1, img)%matrix, 0.0_dp)
560 : END DO
561 366 : CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
562 : END IF
563 :
564 18707 : END SUBROUTINE qs_matrix_h_allocate
565 :
566 : ! **************************************************************************************************
567 : !> \brief (Re-)allocates matrix_h_im from matrix_h
568 : !> \param qs_env ...
569 : ! **************************************************************************************************
570 8 : SUBROUTINE qs_matrix_h_allocate_imag_from_real(qs_env)
571 : TYPE(qs_environment_type) :: qs_env
572 :
573 : CHARACTER(LEN=default_string_length) :: headline
574 : INTEGER :: image, nimages
575 8 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im
576 : TYPE(dbcsr_type), POINTER :: template
577 : TYPE(dft_control_type), POINTER :: dft_control
578 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
579 8 : POINTER :: sab_orb
580 : TYPE(qs_ks_env_type), POINTER :: ks_env
581 :
582 8 : NULLIFY (matrix_h_im, matrix_h, dft_control, template, sab_orb, ks_env)
583 :
584 : CALL get_qs_env(qs_env, &
585 : matrix_h_im_kp=matrix_h_im, &
586 : matrix_h_kp=matrix_h, &
587 : dft_control=dft_control, &
588 : sab_orb=sab_orb, &
589 8 : ks_env=ks_env)
590 :
591 8 : nimages = dft_control%nimages
592 :
593 8 : CPASSERT(nimages == SIZE(matrix_h, 2))
594 :
595 8 : CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
596 :
597 16 : DO image = 1, nimages
598 8 : headline = "IMAGINARY CORE HAMILTONIAN MATRIX"
599 8 : ALLOCATE (matrix_h_im(1, image)%matrix)
600 8 : template => matrix_h(1, image)%matrix ! base on real part, but anti-symmetric
601 : CALL dbcsr_create(matrix=matrix_h_im(1, image)%matrix, template=template, &
602 8 : name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric)
603 8 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, image)%matrix, sab_orb)
604 16 : CALL dbcsr_set(matrix_h_im(1, image)%matrix, 0.0_dp)
605 : END DO
606 8 : CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
607 :
608 8 : END SUBROUTINE qs_matrix_h_allocate_imag_from_real
609 :
610 : END MODULE qs_core_hamiltonian
|
