Line data Source code
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 : !> <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 16625 : 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 16625 : 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 16625 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_t, &
169 16625 : 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 16625 : POINTER :: sab_orb, sap_oce
174 : TYPE(oce_matrix_type), POINTER :: oce
175 16625 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
176 16625 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
177 16625 : 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 33250 : IF (calculate_forces) THEN
183 5875 : CALL timeset(routineN//"_forces", handle)
184 : ELSE
185 10750 : CALL timeset(routineN, handle)
186 : END IF
187 :
188 16625 : NULLIFY (logger)
189 16625 : logger => cp_get_default_logger()
190 :
191 16625 : NULLIFY (dft_control)
192 16625 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
193 :
194 : ! is this a orbital-free method calculation
195 16625 : ofdft = dft_control%qs_control%ofgpw
196 :
197 16625 : nimages = dft_control%nimages
198 16625 : IF (ofdft) THEN
199 0 : CPASSERT(nimages == 1)
200 : END IF
201 :
202 16625 : nders = 0
203 16625 : IF (calculate_forces) THEN
204 5875 : nder = 1
205 : ELSE
206 10750 : 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 10746 : nder = 0
211 : END IF
212 : END IF
213 :
214 16625 : 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 16625 : IF (ASSOCIATED(dft_control%rtp_control)) THEN
224 1788 : IF (dft_control%rtp_control%apply_delta_pulse .AND. &
225 : dft_control%rtp_control%periodic) THEN
226 116 : nders = 1
227 : END IF
228 : END IF
229 :
230 16625 : IF (dft_control%tddfpt2_control%enabled) THEN
231 1458 : nders = 1
232 1458 : IF (dft_control%do_admm) THEN
233 326 : IF (dft_control%admm_control%purification_method /= do_admm_purify_none) &
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 :
239 : ! filter for new matrices
240 16625 : eps_filter = dft_control%qs_control%eps_filter_matrix
241 : !
242 16625 : NULLIFY (ks_env)
243 16625 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
244 16625 : NULLIFY (matrix_s, matrix_t)
245 16625 : CALL get_qs_env(qs_env=qs_env, kinetic_kp=matrix_t, matrix_s_kp=matrix_s)
246 16625 : NULLIFY (sab_orb)
247 16625 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
248 16625 : NULLIFY (rho, force, matrix_p, matrix_w)
249 16625 : IF (calculate_forces) THEN
250 5875 : CALL get_qs_env(qs_env=qs_env, force=force, matrix_w_kp=matrix_w)
251 5875 : CALL get_qs_env(qs_env=qs_env, rho=rho)
252 5875 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
253 : ! *** If LSD, then combine alpha density and beta density to
254 : ! *** total density: alpha <- alpha + beta and
255 : ! *** spin density: beta <- alpha - beta
256 : ! (since all things can be computed based on the sum of these matrices anyway)
257 : ! (matrix_p is restored at the end of the run, matrix_w is left in its modified state
258 : ! (as it should not be needed afterwards)
259 5875 : IF (SIZE(matrix_p, 1) == 2) THEN
260 2422 : DO img = 1, nimages
261 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
262 1592 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
263 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
264 1592 : alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
265 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
266 2422 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
267 : END DO
268 : END IF
269 : ELSE
270 : NULLIFY (matrix_p, matrix_w)
271 : END IF
272 :
273 : ! S matrix
274 : CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, &
275 : matrix_name="OVERLAP MATRIX", &
276 : basis_type_a="ORB", &
277 : basis_type_b="ORB", &
278 : sab_nl=sab_orb, calculate_forces=calculate_forces, &
279 16625 : matrixkp_p=matrix_w)
280 :
281 16625 : IF (calculate_forces) THEN
282 : ! *** If LSD, then recover alpha density and beta density ***
283 : ! *** from the total density (1) and the spin density (2) ***
284 : ! *** The W matrix is neglected, since it will be destroyed ***
285 : ! *** in the calling force routine after leaving this routine ***
286 5875 : IF (SIZE(matrix_p, 1) == 2) THEN
287 2422 : DO img = 1, nimages
288 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
289 1592 : alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
290 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
291 2422 : alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
292 : END DO
293 : END IF
294 : END IF
295 :
296 : ! T matrix
297 : CALL kinetic_energy_matrix(qs_env, matrixkp_t=matrix_t, &
298 : matrix_p=matrix_p, &
299 : matrix_name="KINETIC ENERGY MATRIX", &
300 : basis_type="ORB", &
301 : sab_orb=sab_orb, &
302 : calculate_forces=calculate_forces, &
303 16625 : eps_filter=eps_filter)
304 :
305 : ! (Re-)allocate H matrix based on overlap matrix
306 16625 : CALL get_ks_env(ks_env, complex_ks=h_is_complex)
307 16625 : CALL qs_matrix_h_allocate(qs_env, matrix_s(1, 1)%matrix, is_complex=h_is_complex)
308 :
309 16625 : NULLIFY (matrix_h)
310 16625 : CALL get_qs_env(qs_env, matrix_h_kp=matrix_h)
311 :
312 16625 : IF (.NOT. ofdft) THEN
313 66254 : DO img = 1, nimages
314 : CALL dbcsr_copy(matrix_h(1, img)%matrix, matrix_t(1, img)%matrix, &
315 66254 : keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX")
316 : END DO
317 : END IF
318 :
319 16625 : NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
320 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
321 16625 : particle_set=particle_set)
322 :
323 : ! *** core and pseudopotentials
324 16625 : CALL core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder)
325 :
326 : ! *** CNEO nuclear V_core
327 16625 : CALL cneo_core_matrices(qs_env, calculate_forces, nder)
328 :
329 : ! *** GAPW one-center-expansion (oce) matrices
330 16625 : NULLIFY (sap_oce)
331 16625 : CALL get_qs_env(qs_env=qs_env, sap_oce=sap_oce)
332 16625 : NULLIFY (oce)
333 16625 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
334 2256 : CALL get_qs_env(qs_env=qs_env, oce=oce)
335 2256 : CALL create_oce_set(oce)
336 2256 : nkind = SIZE(atomic_kind_set)
337 2256 : CALL allocate_oce_set(oce, nkind)
338 2256 : eps_fit = dft_control%qs_control%gapw_control%eps_fit
339 2256 : IF (ASSOCIATED(sap_oce)) &
340 : CALL build_oce_matrices(oce%intac, calculate_forces, nder, qs_kind_set, particle_set, &
341 2222 : sap_oce, eps_fit)
342 : END IF
343 :
344 : ! *** KG atomic potentials for nonadditive kinetic energy
345 16625 : IF (dft_control%qs_control%do_kg) THEN
346 182 : IF (qs_env%kg_env%tnadd_method == kg_tnadd_atomic) THEN
347 30 : CALL get_qs_env(qs_env=qs_env, kg_env=kg_env, virial=virial, dbcsr_dist=dbcsr_dist)
348 30 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
349 : CALL build_tnadd_mat(kg_env, matrix_p, force, virial, calculate_forces, use_virial, &
350 30 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, dbcsr_dist)
351 : END IF
352 : END IF
353 :
354 : ! *** Put the core Hamiltonian matrix in the QS environment ***
355 16625 : CALL set_qs_env(qs_env, oce=oce)
356 16625 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s, kinetic_kp=matrix_t, matrix_h_kp=matrix_h)
357 :
358 : ! *** Print matrices if requested
359 16625 : CALL dump_info_core_hamiltonian(qs_env, calculate_forces)
360 :
361 : ! *** Overlap condition number
362 16625 : IF (.NOT. calculate_forces) THEN
363 10750 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
364 : "DFT%PRINT%OVERLAP_CONDITION") /= 0) THEN
365 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
366 36 : extension=".Log")
367 36 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
368 36 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
369 36 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
370 36 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
371 36 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
372 : END IF
373 : END IF
374 :
375 16625 : CALL timestop(handle)
376 :
377 16625 : END SUBROUTINE build_core_hamiltonian_matrix
378 :
379 : ! **************************************************************************************************
380 : !> \brief Possibly prints matrices after the construction of the Core
381 : !> Hamiltonian Matrix
382 : !> \param qs_env ...
383 : !> \param calculate_forces ...
384 : ! **************************************************************************************************
385 33250 : SUBROUTINE dump_info_core_hamiltonian(qs_env, calculate_forces)
386 : TYPE(qs_environment_type), POINTER :: qs_env
387 : LOGICAL, INTENT(IN) :: calculate_forces
388 :
389 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_info_core_hamiltonian'
390 :
391 : INTEGER :: after, handle, i, ic, iw, output_unit
392 : LOGICAL :: omit_headers
393 : TYPE(cp_logger_type), POINTER :: logger
394 16625 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_v
395 16625 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_h, matrixkp_s, matrixkp_t
396 : TYPE(mp_para_env_type), POINTER :: para_env
397 :
398 16625 : CALL timeset(routineN, handle)
399 :
400 16625 : NULLIFY (logger, matrix_v, para_env)
401 16625 : logger => cp_get_default_logger()
402 16625 : CALL get_qs_env(qs_env, para_env=para_env)
403 :
404 : ! Print the distribution of the overlap matrix blocks
405 : ! this duplicates causes duplicate printing at the force calc
406 16625 : IF (.NOT. calculate_forces) THEN
407 10750 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
408 : qs_env%input, "PRINT%DISTRIBUTION"), cp_p_file)) THEN
409 : output_unit = cp_print_key_unit_nr(logger, qs_env%input, "PRINT%DISTRIBUTION", &
410 92 : extension=".distribution")
411 92 : CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
412 92 : CALL cp_dbcsr_write_matrix_dist(matrixkp_s(1, 1)%matrix, output_unit, para_env)
413 92 : CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, "PRINT%DISTRIBUTION")
414 : END IF
415 : END IF
416 :
417 16625 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
418 : ! Print the overlap integral matrix, if requested
419 16625 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
420 : qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
421 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
422 4 : extension=".Log")
423 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
424 4 : after = MIN(MAX(after, 1), 16)
425 4 : CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
426 4 : IF (ASSOCIATED(matrixkp_s)) THEN
427 8 : DO ic = 1, SIZE(matrixkp_s, 2)
428 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(1, ic)%matrix, 4, after, qs_env, para_env, &
429 8 : output_unit=iw, omit_headers=omit_headers)
430 : END DO
431 4 : IF (BTEST(cp_print_key_should_output(logger%iter_info, qs_env%input, &
432 : "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
433 8 : DO ic = 1, SIZE(matrixkp_s, 2)
434 20 : DO i = 2, SIZE(matrixkp_s, 1)
435 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(i, ic)%matrix, 4, after, qs_env, para_env, &
436 16 : output_unit=iw, omit_headers=omit_headers)
437 : END DO
438 : END DO
439 : END IF
440 : END IF
441 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
442 4 : "DFT%PRINT%AO_MATRICES/OVERLAP")
443 : END IF
444 :
445 : ! Print the kinetic energy integral matrix, if requested
446 16625 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
447 : qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY"), cp_p_file)) THEN
448 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY", &
449 48 : extension=".Log")
450 48 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
451 48 : after = MIN(MAX(after, 1), 16)
452 48 : CALL get_qs_env(qs_env, kinetic_kp=matrixkp_t)
453 48 : IF (ASSOCIATED(matrixkp_t)) THEN
454 96 : DO ic = 1, SIZE(matrixkp_t, 2)
455 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_t(1, ic)%matrix, 4, after, qs_env, para_env, &
456 96 : output_unit=iw, omit_headers=omit_headers)
457 : END DO
458 : END IF
459 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
460 48 : "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY")
461 : END IF
462 :
463 : ! Print the potential energy matrix, if requested
464 16625 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
465 : qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY"), cp_p_file)) THEN
466 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY", &
467 48 : extension=".Log")
468 48 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
469 48 : after = MIN(MAX(after, 1), 16)
470 48 : CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h, kinetic_kp=matrixkp_t)
471 48 : IF (ASSOCIATED(matrixkp_h)) THEN
472 48 : IF (SIZE(matrixkp_h, 2) == 1) THEN
473 48 : CALL dbcsr_allocate_matrix_set(matrix_v, 1)
474 48 : ALLOCATE (matrix_v(1)%matrix)
475 48 : CALL dbcsr_copy(matrix_v(1)%matrix, matrixkp_h(1, 1)%matrix, name="POTENTIAL ENERGY MATRIX")
476 : CALL dbcsr_add(matrix_v(1)%matrix, matrixkp_t(1, 1)%matrix, &
477 48 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
478 : CALL cp_dbcsr_write_sparse_matrix(matrix_v(1)%matrix, 4, after, qs_env, &
479 48 : para_env, output_unit=iw, omit_headers=omit_headers)
480 48 : CALL dbcsr_deallocate_matrix_set(matrix_v)
481 : ELSE
482 0 : CPWARN("Printing of potential energy matrix not implemented for k-points")
483 : END IF
484 : END IF
485 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
486 48 : "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY")
487 : END IF
488 :
489 : ! Print the core Hamiltonian matrix, if requested
490 16625 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
491 : qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
492 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
493 48 : extension=".Log")
494 48 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
495 48 : after = MIN(MAX(after, 1), 16)
496 48 : CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h)
497 48 : IF (ASSOCIATED(matrixkp_h)) THEN
498 96 : DO ic = 1, SIZE(matrixkp_h, 2)
499 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_h(1, ic)%matrix, 4, after, qs_env, para_env, &
500 96 : output_unit=iw, omit_headers=omit_headers)
501 : END DO
502 : END IF
503 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
504 48 : "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
505 : END IF
506 :
507 16625 : CALL timestop(handle)
508 :
509 16625 : END SUBROUTINE dump_info_core_hamiltonian
510 :
511 : ! **************************************************************************************************
512 : !> \brief (Re-)allocate matrix_h based on the template (typically the overlap matrix)
513 : !> \param qs_env ...
514 : !> \param template ...
515 : !> \param is_complex ...
516 : ! **************************************************************************************************
517 16625 : SUBROUTINE qs_matrix_h_allocate(qs_env, template, is_complex)
518 : TYPE(qs_environment_type) :: qs_env
519 : TYPE(dbcsr_type), INTENT(in) :: template
520 : LOGICAL, INTENT(in) :: is_complex
521 :
522 : CHARACTER(LEN=default_string_length) :: headline
523 : INTEGER :: img, nimages
524 16625 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im
525 : TYPE(dft_control_type), POINTER :: dft_control
526 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
527 16625 : POINTER :: sab_orb
528 : TYPE(qs_ks_env_type), POINTER :: ks_env
529 :
530 16625 : NULLIFY (matrix_h, matrix_h_im, sab_orb, dft_control, ks_env)
531 : CALL get_qs_env(qs_env=qs_env, &
532 : matrix_h_kp=matrix_h, &
533 : matrix_h_im_kp=matrix_h_im, &
534 : sab_orb=sab_orb, &
535 : dft_control=dft_control, &
536 16625 : ks_env=ks_env)
537 :
538 16625 : nimages = dft_control%nimages
539 16625 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimages)
540 16625 : headline = "CORE HAMILTONIAN MATRIX"
541 66254 : DO img = 1, nimages
542 49629 : ALLOCATE (matrix_h(1, img)%matrix)
543 49629 : CALL dbcsr_create(matrix_h(1, img)%matrix, name=TRIM(headline), template=template)
544 49629 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
545 66254 : CALL dbcsr_set(matrix_h(1, img)%matrix, 0.0_dp)
546 : END DO
547 16625 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
548 :
549 16625 : IF (is_complex) THEN
550 352 : headline = "IMAGINARY PART OF CORE HAMILTONIAN MATRIX"
551 352 : CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
552 704 : DO img = 1, nimages
553 352 : ALLOCATE (matrix_h_im(1, img)%matrix)
554 : CALL dbcsr_create(matrix_h_im(1, img)%matrix, name=TRIM(headline), template=template, &
555 352 : matrix_type=dbcsr_type_antisymmetric)
556 352 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, img)%matrix, sab_orb)
557 704 : CALL dbcsr_set(matrix_h_im(1, img)%matrix, 0.0_dp)
558 : END DO
559 352 : CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
560 : END IF
561 :
562 16625 : END SUBROUTINE qs_matrix_h_allocate
563 :
564 : ! **************************************************************************************************
565 : !> \brief (Re-)allocates matrix_h_im from matrix_h
566 : !> \param qs_env ...
567 : ! **************************************************************************************************
568 8 : SUBROUTINE qs_matrix_h_allocate_imag_from_real(qs_env)
569 : TYPE(qs_environment_type) :: qs_env
570 :
571 : CHARACTER(LEN=default_string_length) :: headline
572 : INTEGER :: image, nimages
573 8 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im
574 : TYPE(dbcsr_type), POINTER :: template
575 : TYPE(dft_control_type), POINTER :: dft_control
576 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
577 8 : POINTER :: sab_orb
578 : TYPE(qs_ks_env_type), POINTER :: ks_env
579 :
580 8 : NULLIFY (matrix_h_im, matrix_h, dft_control, template, sab_orb, ks_env)
581 :
582 : CALL get_qs_env(qs_env, &
583 : matrix_h_im_kp=matrix_h_im, &
584 : matrix_h_kp=matrix_h, &
585 : dft_control=dft_control, &
586 : sab_orb=sab_orb, &
587 8 : ks_env=ks_env)
588 :
589 8 : nimages = dft_control%nimages
590 :
591 8 : CPASSERT(nimages == SIZE(matrix_h, 2))
592 :
593 8 : CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
594 :
595 16 : DO image = 1, nimages
596 8 : headline = "IMAGINARY CORE HAMILTONIAN MATRIX"
597 8 : ALLOCATE (matrix_h_im(1, image)%matrix)
598 8 : template => matrix_h(1, image)%matrix ! base on real part, but anti-symmetric
599 : CALL dbcsr_create(matrix=matrix_h_im(1, image)%matrix, template=template, &
600 8 : name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric)
601 8 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, image)%matrix, sab_orb)
602 16 : CALL dbcsr_set(matrix_h_im(1, image)%matrix, 0.0_dp)
603 : END DO
604 8 : CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
605 :
606 8 : END SUBROUTINE qs_matrix_h_allocate_imag_from_real
607 :
608 : END MODULE qs_core_hamiltonian
|
