Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 : get_atomic_kind_set
80 : USE core_ae, ONLY: build_core_ae
81 : USE core_ppl, ONLY: build_core_ppl
82 : USE core_ppnl, ONLY: build_core_ppnl
83 : USE cp_blacs_env, ONLY: cp_blacs_env_type
84 : USE cp_control_types, ONLY: dft_control_type
85 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
86 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
87 : dbcsr_deallocate_matrix_set
88 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_matrix_dist,&
89 : cp_dbcsr_write_sparse_matrix
90 : USE cp_log_handling, ONLY: cp_get_default_logger,&
91 : cp_logger_type
92 : USE cp_output_handling, ONLY: cp_p_file,&
93 : cp_print_key_finished_output,&
94 : cp_print_key_should_output,&
95 : cp_print_key_unit_nr
96 : USE dbcsr_api, ONLY: &
97 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_iterator_blocks_left, &
98 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
99 : dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric
100 : USE input_constants, ONLY: do_admm_purify_none,&
101 : do_ppl_analytic,&
102 : kg_tnadd_atomic,&
103 : rel_none,&
104 : rel_trans_atom
105 : USE input_section_types, ONLY: section_vals_val_get
106 : USE kg_environment_types, ONLY: kg_environment_type
107 : USE kg_tnadd_mat, ONLY: build_tnadd_mat
108 : USE kinds, ONLY: default_string_length,&
109 : dp
110 : USE kpoint_types, ONLY: get_kpoint_info,&
111 : kpoint_type
112 : USE lri_environment_types, ONLY: lri_environment_type
113 : USE message_passing, ONLY: mp_para_env_type
114 : USE particle_types, ONLY: particle_type
115 : USE qs_condnum, ONLY: overlap_condnum
116 : USE qs_environment_types, ONLY: get_qs_env,&
117 : qs_environment_type,&
118 : set_qs_env
119 : USE qs_force_types, ONLY: qs_force_type
120 : USE qs_kind_types, ONLY: get_qs_kind,&
121 : qs_kind_type
122 : USE qs_kinetic, ONLY: build_kinetic_matrix
123 : USE qs_ks_types, ONLY: get_ks_env,&
124 : qs_ks_env_type,&
125 : set_ks_env
126 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
127 : USE qs_oce_methods, ONLY: build_oce_matrices
128 : USE qs_oce_types, ONLY: allocate_oce_set,&
129 : create_oce_set,&
130 : oce_matrix_type
131 : USE qs_overlap, ONLY: build_overlap_matrix
132 : USE qs_rho_types, ONLY: qs_rho_get,&
133 : qs_rho_type
134 : USE virial_types, ONLY: virial_type
135 : #include "./base/base_uses.f90"
136 :
137 : IMPLICIT NONE
138 :
139 : PRIVATE
140 :
141 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_hamiltonian'
142 :
143 : PUBLIC :: build_core_hamiltonian_matrix
144 : PUBLIC :: dump_info_core_hamiltonian, qs_matrix_h_allocate_imag_from_real
145 :
146 : CONTAINS
147 :
148 : ! **************************************************************************************************
149 : !> \brief Cosntruction of the QS Core Hamiltonian Matrix
150 : !> \param qs_env ...
151 : !> \param calculate_forces ...
152 : !> \author Creation (11.03.2002,MK)
153 : !> Non-redundant calculation of the non-local part of the GTH PP (22.05.2003,MK)
154 : !> New parallelization scheme (27.06.2003,MK)
155 : ! **************************************************************************************************
156 15891 : SUBROUTINE build_core_hamiltonian_matrix(qs_env, calculate_forces)
157 :
158 : TYPE(qs_environment_type), POINTER :: qs_env
159 : LOGICAL, INTENT(IN) :: calculate_forces
160 :
161 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_hamiltonian_matrix'
162 :
163 : INTEGER :: handle, ic, img, iw, nder, nders, &
164 : nimages, nkind
165 15891 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
166 : LOGICAL :: h_is_complex, norml1, norml2, ofdft, &
167 : use_arnoldi, use_virial
168 : REAL(KIND=dp) :: eps_filter, eps_fit
169 : REAL(KIND=dp), DIMENSION(2) :: condnum
170 15891 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
171 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
172 : TYPE(cp_logger_type), POINTER :: logger
173 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
174 15891 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_t, &
175 15891 : matrix_w
176 : TYPE(dft_control_type), POINTER :: dft_control
177 : TYPE(kg_environment_type), POINTER :: kg_env
178 : TYPE(kpoint_type), POINTER :: kpoints
179 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
180 15891 : POINTER :: sab_orb, sap_oce
181 : TYPE(oce_matrix_type), POINTER :: oce
182 15891 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
183 15891 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
184 15891 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
185 : TYPE(qs_ks_env_type), POINTER :: ks_env
186 : TYPE(qs_rho_type), POINTER :: rho
187 : TYPE(virial_type), POINTER :: virial
188 :
189 31782 : IF (calculate_forces) THEN
190 5645 : CALL timeset(routineN//"_forces", handle)
191 : ELSE
192 10246 : CALL timeset(routineN, handle)
193 : END IF
194 :
195 15891 : NULLIFY (logger)
196 15891 : logger => cp_get_default_logger()
197 :
198 15891 : NULLIFY (dft_control)
199 15891 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
200 :
201 : ! is this a orbital-free method calculation
202 15891 : ofdft = dft_control%qs_control%ofgpw
203 :
204 15891 : nimages = dft_control%nimages
205 15891 : IF (ofdft) THEN
206 0 : CPASSERT(nimages == 1)
207 : END IF
208 :
209 15891 : nders = 0
210 15891 : IF (calculate_forces) THEN
211 5645 : nder = 1
212 : ELSE
213 10246 : 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 10242 : nder = 0
218 : END IF
219 : END IF
220 :
221 15891 : 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 15891 : IF (ASSOCIATED(dft_control%rtp_control)) THEN
231 1782 : IF (dft_control%rtp_control%apply_delta_pulse .AND. &
232 : dft_control%rtp_control%periodic) THEN
233 116 : nders = 1
234 : END IF
235 : END IF
236 :
237 15891 : IF (dft_control%tddfpt2_control%enabled) THEN
238 1396 : nders = 1
239 1396 : IF (dft_control%do_admm) THEN
240 314 : IF (dft_control%admm_control%purification_method /= do_admm_purify_none) &
241 : CALL cp_abort(__LOCATION__, &
242 0 : "Only purification method NONE is possible with TDDFT at the moment")
243 : END IF
244 : END IF
245 :
246 : ! filter for new matrices
247 15891 : eps_filter = dft_control%qs_control%eps_filter_matrix
248 : !
249 15891 : NULLIFY (ks_env)
250 15891 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
251 15891 : NULLIFY (matrix_s, matrix_t)
252 15891 : CALL get_qs_env(qs_env=qs_env, kinetic_kp=matrix_t, matrix_s_kp=matrix_s)
253 15891 : NULLIFY (sab_orb)
254 15891 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
255 15891 : NULLIFY (rho, force, matrix_p, matrix_w)
256 15891 : IF (calculate_forces) THEN
257 5645 : CALL get_qs_env(qs_env=qs_env, force=force, matrix_w_kp=matrix_w)
258 5645 : CALL get_qs_env(qs_env=qs_env, rho=rho)
259 5645 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
260 : ! *** If LSD, then combine alpha density and beta density to
261 : ! *** total density: alpha <- alpha + beta and
262 : ! *** spin density: beta <- alpha - beta
263 : ! (since all things can be computed based on the sum of these matrices anyway)
264 : ! (matrix_p is restored at the end of the run, matrix_w is left in its modified state
265 : ! (as it should not be needed afterwards)
266 5645 : IF (SIZE(matrix_p, 1) == 2) THEN
267 2346 : DO img = 1, nimages
268 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
269 1536 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
270 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
271 1536 : alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
272 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
273 2346 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
274 : END DO
275 : END IF
276 : ! S matrix
277 : CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, &
278 : matrix_name="OVERLAP MATRIX", &
279 : basis_type_a="ORB", &
280 : basis_type_b="ORB", &
281 : sab_nl=sab_orb, calculate_forces=.TRUE., &
282 5645 : matrixkp_p=matrix_w)
283 : ! T matrix
284 5645 : IF (.NOT. ofdft) &
285 : CALL build_kinetic_matrix(ks_env, matrixkp_t=matrix_t, &
286 : matrix_name="KINETIC ENERGY MATRIX", &
287 : basis_type="ORB", &
288 : sab_nl=sab_orb, calculate_forces=.TRUE., &
289 : matrixkp_p=matrix_p, &
290 5645 : eps_filter=eps_filter)
291 : ! *** If LSD, then recover alpha density and beta density ***
292 : ! *** from the total density (1) and the spin density (2) ***
293 : ! *** The W matrix is neglected, since it will be destroyed ***
294 : ! *** in the calling force routine after leaving this routine ***
295 5645 : IF (SIZE(matrix_p, 1) == 2) THEN
296 2346 : DO img = 1, nimages
297 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
298 1536 : alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
299 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
300 2346 : alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
301 : END DO
302 : END IF
303 : ELSE
304 : ! S matrix
305 : CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, &
306 : matrix_name="OVERLAP MATRIX", &
307 : basis_type_a="ORB", &
308 : basis_type_b="ORB", &
309 10246 : sab_nl=sab_orb)
310 : ! T matrix
311 10246 : IF (.NOT. ofdft) &
312 : CALL build_kinetic_matrix(ks_env, matrixkp_t=matrix_t, &
313 : matrix_name="KINETIC ENERGY MATRIX", &
314 : basis_type="ORB", &
315 : sab_nl=sab_orb, &
316 10246 : eps_filter=eps_filter)
317 : END IF
318 :
319 : ! (Re-)allocate H matrix based on overlap matrix
320 15891 : CALL get_ks_env(ks_env, complex_ks=h_is_complex)
321 15891 : CALL qs_matrix_h_allocate(qs_env, matrix_s(1, 1)%matrix, is_complex=h_is_complex)
322 :
323 15891 : NULLIFY (matrix_h)
324 15891 : CALL get_qs_env(qs_env, matrix_h_kp=matrix_h)
325 :
326 15891 : IF (.NOT. ofdft) THEN
327 63066 : DO img = 1, nimages
328 : CALL dbcsr_copy(matrix_h(1, img)%matrix, matrix_t(1, img)%matrix, &
329 63066 : keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX")
330 : END DO
331 : END IF
332 :
333 15891 : NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
334 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
335 15891 : particle_set=particle_set)
336 :
337 15891 : IF (.NOT. ofdft) THEN
338 : ! relativistic atomic correction to kinetic energy
339 15891 : IF (qs_env%rel_control%rel_method /= rel_none) THEN
340 58 : IF (qs_env%rel_control%rel_transformation == rel_trans_atom) THEN
341 58 : IF (nimages == 1) THEN
342 58 : ic = 1
343 : ELSE
344 4 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
345 4 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
346 4 : ic = cell_to_index(0, 0, 0)
347 : END IF
348 58 : CALL build_atomic_relmat(matrix_h(1, 1)%matrix, atomic_kind_set, qs_kind_set)
349 : ELSE
350 0 : CPABORT("Relativistic corrections of this type are currently not implemented")
351 : END IF
352 : END IF
353 : END IF
354 :
355 : ! *** core and pseudopotentials
356 15891 : CALL core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder)
357 :
358 : ! *** GAPW one-center-expansion (oce) matrices
359 15891 : NULLIFY (sap_oce)
360 15891 : CALL get_qs_env(qs_env=qs_env, sap_oce=sap_oce)
361 15891 : NULLIFY (oce)
362 15891 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
363 2086 : CALL get_qs_env(qs_env=qs_env, oce=oce)
364 2086 : CALL create_oce_set(oce)
365 2086 : nkind = SIZE(atomic_kind_set)
366 2086 : CALL allocate_oce_set(oce, nkind)
367 2086 : eps_fit = dft_control%qs_control%gapw_control%eps_fit
368 2086 : IF (ASSOCIATED(sap_oce)) &
369 : CALL build_oce_matrices(oce%intac, calculate_forces, nder, qs_kind_set, particle_set, &
370 2054 : sap_oce, eps_fit)
371 : END IF
372 :
373 : ! *** KG atomic potentials for nonadditive kinetic energy
374 15891 : IF (dft_control%qs_control%do_kg) THEN
375 182 : IF (qs_env%kg_env%tnadd_method == kg_tnadd_atomic) THEN
376 30 : CALL get_qs_env(qs_env=qs_env, kg_env=kg_env, virial=virial, dbcsr_dist=dbcsr_dist)
377 30 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
378 : CALL build_tnadd_mat(kg_env, matrix_p, force, virial, calculate_forces, use_virial, &
379 30 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, dbcsr_dist)
380 : END IF
381 : END IF
382 :
383 : ! *** Put the core Hamiltonian matrix in the QS environment ***
384 15891 : CALL set_qs_env(qs_env, oce=oce)
385 15891 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s, kinetic_kp=matrix_t, matrix_h_kp=matrix_h)
386 :
387 : ! *** Print matrices if requested
388 15891 : CALL dump_info_core_hamiltonian(qs_env, calculate_forces)
389 :
390 : ! *** Overlap condition number
391 15891 : IF (.NOT. calculate_forces) THEN
392 10246 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
393 : "DFT%PRINT%OVERLAP_CONDITION") .NE. 0) THEN
394 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
395 36 : extension=".Log")
396 36 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
397 36 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
398 36 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
399 36 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
400 36 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
401 : END IF
402 : END IF
403 :
404 15891 : CALL timestop(handle)
405 :
406 15891 : END SUBROUTINE build_core_hamiltonian_matrix
407 :
408 : ! **************************************************************************************************
409 : !> \brief ...
410 : !> \param qs_env ...
411 : !> \param matrix_h ...
412 : !> \param matrix_p ...
413 : !> \param calculate_forces ...
414 : !> \param nder ...
415 : ! **************************************************************************************************
416 15891 : SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder)
417 :
418 : TYPE(qs_environment_type), POINTER :: qs_env
419 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
420 : LOGICAL, INTENT(IN) :: calculate_forces
421 : INTEGER, INTENT(IN) :: nder
422 :
423 : INTEGER :: nimages
424 15891 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
425 : LOGICAL :: all_present, my_gt_nl, ppl_present, &
426 : ppnl_present, use_virial
427 : REAL(KIND=dp) :: eps_ppnl
428 15891 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
429 : TYPE(dft_control_type), POINTER :: dft_control
430 : TYPE(kpoint_type), POINTER :: kpoints
431 : TYPE(lri_environment_type), POINTER :: lri_env
432 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
433 15891 : POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
434 15891 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
435 15891 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
436 15891 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
437 : TYPE(qs_ks_env_type), POINTER :: ks_env
438 : TYPE(virial_type), POINTER :: virial
439 :
440 15891 : NULLIFY (dft_control)
441 15891 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control)
442 15891 : nimages = dft_control%nimages
443 :
444 : ! check whether a gauge transformed version of the non-local potential part has to be used
445 15891 : my_gt_nl = .FALSE.
446 15891 : IF (qs_env%run_rtp) THEN
447 1244 : CPASSERT(ASSOCIATED(dft_control%rtp_control))
448 1244 : IF (dft_control%rtp_control%velocity_gauge) THEN
449 54 : my_gt_nl = dft_control%rtp_control%nl_gauge_transform
450 : END IF
451 : END IF
452 :
453 : ! prepare for k-points
454 15891 : NULLIFY (cell_to_index)
455 15891 : IF (nimages > 1) THEN
456 280 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
457 280 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
458 280 : dft_control%qs_control%do_ppl_method = do_ppl_analytic
459 : END IF
460 : ! force analytic ppl calculation for GAPW methods
461 15891 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
462 2086 : dft_control%qs_control%do_ppl_method = do_ppl_analytic
463 : END IF
464 :
465 : ! force
466 15891 : NULLIFY (force)
467 15891 : IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force)
468 : ! virial
469 15891 : CALL get_qs_env(qs_env=qs_env, virial=virial)
470 15891 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
471 :
472 15891 : NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
473 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
474 15891 : particle_set=particle_set)
475 :
476 15891 : NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
477 : CALL get_qs_env(qs_env=qs_env, &
478 : sab_orb=sab_orb, &
479 : sac_ae=sac_ae, &
480 : sac_ppl=sac_ppl, &
481 15891 : sap_ppnl=sap_ppnl)
482 :
483 : ! *** compute the nuclear attraction contribution to the core hamiltonian ***
484 15891 : all_present = ASSOCIATED(sac_ae)
485 15891 : IF (all_present) THEN
486 : CALL build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
487 980 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index)
488 : END IF
489 :
490 : ! *** compute the ppl contribution to the core hamiltonian ***
491 15891 : ppl_present = ASSOCIATED(sac_ppl)
492 15891 : IF (ppl_present) THEN
493 14967 : IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
494 14943 : IF (dft_control%qs_control%lrigpw) THEN
495 80 : CALL get_qs_env(qs_env, lri_env=lri_env)
496 80 : IF (lri_env%ppl_ri) THEN
497 4 : IF (lri_env%exact_1c_terms) THEN
498 0 : CPABORT("not implemented")
499 : END IF
500 : ELSE
501 : CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
502 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
503 76 : nimages, cell_to_index, "ORB")
504 : END IF
505 : ELSE
506 : CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
507 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
508 14863 : nimages, cell_to_index, "ORB")
509 : END IF
510 : END IF
511 : END IF
512 :
513 : ! *** compute the ppnl contribution to the core hamiltonian ***
514 15891 : eps_ppnl = dft_control%qs_control%eps_ppnl
515 15891 : ppnl_present = ASSOCIATED(sap_ppnl)
516 15891 : IF (ppnl_present) THEN
517 12100 : IF (.NOT. my_gt_nl) THEN
518 : CALL build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
519 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
520 12064 : nimages, cell_to_index, "ORB")
521 : END IF
522 : END IF
523 :
524 15891 : END SUBROUTINE core_matrices
525 :
526 : ! **************************************************************************************************
527 : !> \brief Adds atomic blocks of relativistic correction for the kinetic energy
528 : !> \param matrix_h ...
529 : !> \param atomic_kind_set ...
530 : !> \param qs_kind_set ...
531 : ! **************************************************************************************************
532 58 : SUBROUTINE build_atomic_relmat(matrix_h, atomic_kind_set, qs_kind_set)
533 : TYPE(dbcsr_type), POINTER :: matrix_h
534 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
535 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
536 :
537 : INTEGER :: blk, iatom, ikind, jatom
538 58 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
539 58 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hblock, reltmat
540 : TYPE(dbcsr_iterator_type) :: iter
541 :
542 58 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
543 :
544 58 : CALL dbcsr_iterator_start(iter, matrix_h)
545 221 : DO WHILE (dbcsr_iterator_blocks_left(iter))
546 163 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, hblock, blk)
547 221 : IF (iatom == jatom) THEN
548 83 : ikind = kind_of(iatom)
549 83 : CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat)
550 192683 : IF (ASSOCIATED(reltmat)) hblock = hblock + reltmat
551 : END IF
552 : END DO
553 58 : CALL dbcsr_iterator_stop(iter)
554 :
555 116 : END SUBROUTINE build_atomic_relmat
556 :
557 : ! **************************************************************************************************
558 : !> \brief Possibly prints matrices after the construction of the Core
559 : !> Hamiltonian Matrix
560 : !> \param qs_env ...
561 : !> \param calculate_forces ...
562 : ! **************************************************************************************************
563 31782 : SUBROUTINE dump_info_core_hamiltonian(qs_env, calculate_forces)
564 : TYPE(qs_environment_type), POINTER :: qs_env
565 : LOGICAL, INTENT(IN) :: calculate_forces
566 :
567 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_info_core_hamiltonian'
568 :
569 : INTEGER :: after, handle, i, ic, iw, output_unit
570 : LOGICAL :: omit_headers
571 : TYPE(cp_logger_type), POINTER :: logger
572 15891 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_v
573 15891 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_h, matrixkp_s, matrixkp_t
574 : TYPE(mp_para_env_type), POINTER :: para_env
575 :
576 15891 : CALL timeset(routineN, handle)
577 :
578 15891 : NULLIFY (logger, matrix_v, para_env)
579 15891 : logger => cp_get_default_logger()
580 15891 : CALL get_qs_env(qs_env, para_env=para_env)
581 :
582 : ! Print the distribution of the overlap matrix blocks
583 : ! this duplicates causes duplicate printing at the force calc
584 15891 : IF (.NOT. calculate_forces) THEN
585 10246 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
586 : qs_env%input, "PRINT%DISTRIBUTION"), cp_p_file)) THEN
587 : output_unit = cp_print_key_unit_nr(logger, qs_env%input, "PRINT%DISTRIBUTION", &
588 92 : extension=".distribution")
589 92 : CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
590 92 : CALL cp_dbcsr_write_matrix_dist(matrixkp_s(1, 1)%matrix, output_unit, para_env)
591 92 : CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, "PRINT%DISTRIBUTION")
592 : END IF
593 : END IF
594 :
595 15891 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
596 : ! Print the overlap integral matrix, if requested
597 15891 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
598 : qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
599 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
600 4 : extension=".Log")
601 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
602 4 : after = MIN(MAX(after, 1), 16)
603 4 : CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
604 4 : IF (ASSOCIATED(matrixkp_s)) THEN
605 8 : DO ic = 1, SIZE(matrixkp_s, 2)
606 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(1, ic)%matrix, 4, after, qs_env, para_env, &
607 8 : output_unit=iw, omit_headers=omit_headers)
608 : END DO
609 4 : IF (BTEST(cp_print_key_should_output(logger%iter_info, qs_env%input, &
610 : "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
611 8 : DO ic = 1, SIZE(matrixkp_s, 2)
612 20 : DO i = 2, SIZE(matrixkp_s, 1)
613 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(i, ic)%matrix, 4, after, qs_env, para_env, &
614 16 : output_unit=iw, omit_headers=omit_headers)
615 : END DO
616 : END DO
617 : END IF
618 : END IF
619 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
620 4 : "DFT%PRINT%AO_MATRICES/OVERLAP")
621 : END IF
622 :
623 : ! Print the kinetic energy integral matrix, if requested
624 15891 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
625 : qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY"), cp_p_file)) THEN
626 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY", &
627 92 : extension=".Log")
628 92 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
629 92 : after = MIN(MAX(after, 1), 16)
630 92 : CALL get_qs_env(qs_env, kinetic_kp=matrixkp_t)
631 92 : IF (ASSOCIATED(matrixkp_t)) THEN
632 184 : DO ic = 1, SIZE(matrixkp_t, 2)
633 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_t(1, ic)%matrix, 4, after, qs_env, para_env, &
634 184 : output_unit=iw, omit_headers=omit_headers)
635 : END DO
636 : END IF
637 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
638 92 : "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY")
639 : END IF
640 :
641 : ! Print the potential energy matrix, if requested
642 15891 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
643 : qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY"), cp_p_file)) THEN
644 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY", &
645 92 : extension=".Log")
646 92 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
647 92 : after = MIN(MAX(after, 1), 16)
648 92 : CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h, kinetic_kp=matrixkp_t)
649 92 : IF (ASSOCIATED(matrixkp_h)) THEN
650 92 : IF (SIZE(matrixkp_h, 2) == 1) THEN
651 92 : CALL dbcsr_allocate_matrix_set(matrix_v, 1)
652 92 : ALLOCATE (matrix_v(1)%matrix)
653 92 : CALL dbcsr_copy(matrix_v(1)%matrix, matrixkp_h(1, 1)%matrix, name="POTENTIAL ENERGY MATRIX")
654 : CALL dbcsr_add(matrix_v(1)%matrix, matrixkp_t(1, 1)%matrix, &
655 92 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
656 : CALL cp_dbcsr_write_sparse_matrix(matrix_v(1)%matrix, 4, after, qs_env, &
657 92 : para_env, output_unit=iw, omit_headers=omit_headers)
658 92 : CALL dbcsr_deallocate_matrix_set(matrix_v)
659 : ELSE
660 0 : CPWARN("Printing of potential energy matrix not implemented for k-points")
661 : END IF
662 : END IF
663 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
664 92 : "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY")
665 : END IF
666 :
667 : ! Print the core Hamiltonian matrix, if requested
668 15891 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
669 : qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
670 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
671 92 : extension=".Log")
672 92 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
673 92 : after = MIN(MAX(after, 1), 16)
674 92 : CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h)
675 92 : IF (ASSOCIATED(matrixkp_h)) THEN
676 184 : DO ic = 1, SIZE(matrixkp_h, 2)
677 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_h(1, ic)%matrix, 4, after, qs_env, para_env, &
678 184 : output_unit=iw, omit_headers=omit_headers)
679 : END DO
680 : END IF
681 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
682 92 : "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
683 : END IF
684 :
685 15891 : CALL timestop(handle)
686 :
687 15891 : END SUBROUTINE dump_info_core_hamiltonian
688 :
689 : ! **************************************************************************************************
690 : !> \brief (Re-)allocate matrix_h based on the template (typically the overlap matrix)
691 : !> \param qs_env ...
692 : !> \param template ...
693 : !> \param is_complex ...
694 : ! **************************************************************************************************
695 15891 : SUBROUTINE qs_matrix_h_allocate(qs_env, template, is_complex)
696 : TYPE(qs_environment_type) :: qs_env
697 : TYPE(dbcsr_type), INTENT(in) :: template
698 : LOGICAL, INTENT(in) :: is_complex
699 :
700 : CHARACTER(LEN=default_string_length) :: headline
701 : INTEGER :: img, nimages
702 15891 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im
703 : TYPE(dft_control_type), POINTER :: dft_control
704 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
705 15891 : POINTER :: sab_orb
706 : TYPE(qs_ks_env_type), POINTER :: ks_env
707 :
708 15891 : NULLIFY (matrix_h, matrix_h_im, sab_orb, dft_control, ks_env)
709 : CALL get_qs_env(qs_env=qs_env, &
710 : matrix_h_kp=matrix_h, &
711 : matrix_h_im_kp=matrix_h_im, &
712 : sab_orb=sab_orb, &
713 : dft_control=dft_control, &
714 15891 : ks_env=ks_env)
715 :
716 15891 : nimages = dft_control%nimages
717 15891 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimages)
718 15891 : headline = "CORE HAMILTONIAN MATRIX"
719 63066 : DO img = 1, nimages
720 47175 : ALLOCATE (matrix_h(1, img)%matrix)
721 47175 : CALL dbcsr_create(matrix_h(1, img)%matrix, name=TRIM(headline), template=template)
722 47175 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
723 63066 : CALL dbcsr_set(matrix_h(1, img)%matrix, 0.0_dp)
724 : END DO
725 15891 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
726 :
727 15891 : IF (is_complex) THEN
728 352 : headline = "IMAGINARY PART OF CORE HAMILTONIAN MATRIX"
729 352 : CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
730 704 : DO img = 1, nimages
731 352 : ALLOCATE (matrix_h_im(1, img)%matrix)
732 : CALL dbcsr_create(matrix_h_im(1, img)%matrix, name=TRIM(headline), template=template, &
733 352 : matrix_type=dbcsr_type_antisymmetric)
734 352 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, img)%matrix, sab_orb)
735 704 : CALL dbcsr_set(matrix_h_im(1, img)%matrix, 0.0_dp)
736 : END DO
737 352 : CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
738 : END IF
739 :
740 15891 : END SUBROUTINE qs_matrix_h_allocate
741 :
742 : ! **************************************************************************************************
743 : !> \brief (Re-)allocates matrix_h_im from matrix_h
744 : !> \param qs_env ...
745 : ! **************************************************************************************************
746 8 : SUBROUTINE qs_matrix_h_allocate_imag_from_real(qs_env)
747 : TYPE(qs_environment_type) :: qs_env
748 :
749 : CHARACTER(LEN=default_string_length) :: headline
750 : INTEGER :: image, nimages
751 8 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im
752 : TYPE(dbcsr_type), POINTER :: template
753 : TYPE(dft_control_type), POINTER :: dft_control
754 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
755 8 : POINTER :: sab_orb
756 : TYPE(qs_ks_env_type), POINTER :: ks_env
757 :
758 8 : NULLIFY (matrix_h_im, matrix_h, dft_control, template, sab_orb, ks_env)
759 :
760 : CALL get_qs_env(qs_env, &
761 : matrix_h_im_kp=matrix_h_im, &
762 : matrix_h_kp=matrix_h, &
763 : dft_control=dft_control, &
764 : sab_orb=sab_orb, &
765 8 : ks_env=ks_env)
766 :
767 8 : nimages = dft_control%nimages
768 :
769 8 : CPASSERT(nimages .EQ. SIZE(matrix_h, 2))
770 :
771 8 : CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
772 :
773 16 : DO image = 1, nimages
774 8 : headline = "IMAGINARY CORE HAMILTONIAN MATRIX"
775 8 : ALLOCATE (matrix_h_im(1, image)%matrix)
776 8 : template => matrix_h(1, image)%matrix ! base on real part, but anti-symmetric
777 : CALL dbcsr_create(matrix=matrix_h_im(1, image)%matrix, template=template, &
778 8 : name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric, nze=0)
779 8 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, image)%matrix, sab_orb)
780 16 : CALL dbcsr_set(matrix_h_im(1, image)%matrix, 0.0_dp)
781 : END DO
782 8 : CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
783 :
784 8 : END SUBROUTINE qs_matrix_h_allocate_imag_from_real
785 :
786 : END MODULE qs_core_hamiltonian
|
