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_api, ONLY: &
86 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_iterator_blocks_left, &
87 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
88 : dbcsr_p_type, dbcsr_set, dbcsr_type, 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 : 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, core_matrices
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 16083 : 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 16083 : 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 16083 : 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 16083 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_t, &
175 16083 : 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 16083 : POINTER :: sab_orb, sap_oce
181 : TYPE(oce_matrix_type), POINTER :: oce
182 16083 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
183 16083 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
184 16083 : 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 32166 : IF (calculate_forces) THEN
190 5671 : CALL timeset(routineN//"_forces", handle)
191 : ELSE
192 10412 : CALL timeset(routineN, handle)
193 : END IF
194 :
195 16083 : NULLIFY (logger)
196 16083 : logger => cp_get_default_logger()
197 :
198 16083 : NULLIFY (dft_control)
199 16083 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
200 :
201 : ! is this a orbital-free method calculation
202 16083 : ofdft = dft_control%qs_control%ofgpw
203 :
204 16083 : nimages = dft_control%nimages
205 16083 : IF (ofdft) THEN
206 0 : CPASSERT(nimages == 1)
207 : END IF
208 :
209 16083 : nders = 0
210 16083 : IF (calculate_forces) THEN
211 5671 : nder = 1
212 : ELSE
213 10412 : 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 10408 : nder = 0
218 : END IF
219 : END IF
220 :
221 16083 : 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 16083 : IF (ASSOCIATED(dft_control%rtp_control)) THEN
231 1794 : 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 16083 : 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 16083 : eps_filter = dft_control%qs_control%eps_filter_matrix
248 : !
249 16083 : NULLIFY (ks_env)
250 16083 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
251 16083 : NULLIFY (matrix_s, matrix_t)
252 16083 : CALL get_qs_env(qs_env=qs_env, kinetic_kp=matrix_t, matrix_s_kp=matrix_s)
253 16083 : NULLIFY (sab_orb)
254 16083 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
255 16083 : NULLIFY (rho, force, matrix_p, matrix_w)
256 16083 : IF (calculate_forces) THEN
257 5671 : CALL get_qs_env(qs_env=qs_env, force=force, matrix_w_kp=matrix_w)
258 5671 : CALL get_qs_env(qs_env=qs_env, rho=rho)
259 5671 : 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 5671 : IF (SIZE(matrix_p, 1) == 2) THEN
267 2354 : DO img = 1, nimages
268 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
269 1540 : 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 1540 : 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 2354 : 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 5671 : matrixkp_p=matrix_w)
283 : ! T matrix
284 5671 : 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 5671 : 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 5671 : IF (SIZE(matrix_p, 1) == 2) THEN
296 2354 : DO img = 1, nimages
297 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
298 1540 : 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 2354 : 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 10412 : sab_nl=sab_orb)
310 : ! T matrix
311 10412 : 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 10412 : eps_filter=eps_filter)
317 : END IF
318 :
319 : ! (Re-)allocate H matrix based on overlap matrix
320 16083 : CALL get_ks_env(ks_env, complex_ks=h_is_complex)
321 16083 : CALL qs_matrix_h_allocate(qs_env, matrix_s(1, 1)%matrix, is_complex=h_is_complex)
322 :
323 16083 : NULLIFY (matrix_h)
324 16083 : CALL get_qs_env(qs_env, matrix_h_kp=matrix_h)
325 :
326 16083 : IF (.NOT. ofdft) THEN
327 64758 : DO img = 1, nimages
328 : CALL dbcsr_copy(matrix_h(1, img)%matrix, matrix_t(1, img)%matrix, &
329 64758 : keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX")
330 : END DO
331 : END IF
332 :
333 16083 : 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 16083 : particle_set=particle_set)
336 :
337 16083 : IF (.NOT. ofdft) THEN
338 : ! relativistic atomic correction to kinetic energy
339 16083 : 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 16083 : CALL core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder)
357 :
358 : ! *** GAPW one-center-expansion (oce) matrices
359 16083 : NULLIFY (sap_oce)
360 16083 : CALL get_qs_env(qs_env=qs_env, sap_oce=sap_oce)
361 16083 : NULLIFY (oce)
362 16083 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
363 2090 : CALL get_qs_env(qs_env=qs_env, oce=oce)
364 2090 : CALL create_oce_set(oce)
365 2090 : nkind = SIZE(atomic_kind_set)
366 2090 : CALL allocate_oce_set(oce, nkind)
367 2090 : eps_fit = dft_control%qs_control%gapw_control%eps_fit
368 2090 : IF (ASSOCIATED(sap_oce)) &
369 : CALL build_oce_matrices(oce%intac, calculate_forces, nder, qs_kind_set, particle_set, &
370 2058 : sap_oce, eps_fit)
371 : END IF
372 :
373 : ! *** KG atomic potentials for nonadditive kinetic energy
374 16083 : 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 16083 : CALL set_qs_env(qs_env, oce=oce)
385 16083 : 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 16083 : CALL dump_info_core_hamiltonian(qs_env, calculate_forces)
389 :
390 : ! *** Overlap condition number
391 16083 : IF (.NOT. calculate_forces) THEN
392 10412 : 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 16083 : CALL timestop(handle)
405 :
406 16083 : 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 : !> \param atcore ...
416 : ! **************************************************************************************************
417 16149 : SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, atcore)
418 :
419 : TYPE(qs_environment_type), POINTER :: qs_env
420 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
421 : LOGICAL, INTENT(IN) :: calculate_forces
422 : INTEGER, INTENT(IN) :: nder
423 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atcore
424 :
425 : INTEGER :: natom, nimages
426 16149 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
427 : LOGICAL :: all_present, my_gt_nl, ppl_present, &
428 : ppnl_present, use_virial
429 : REAL(KIND=dp) :: eps_ppnl
430 16149 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
431 : TYPE(dft_control_type), POINTER :: dft_control
432 : TYPE(kpoint_type), POINTER :: kpoints
433 : TYPE(lri_environment_type), POINTER :: lri_env
434 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
435 16149 : POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
436 16149 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
437 16149 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
438 16149 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
439 : TYPE(qs_ks_env_type), POINTER :: ks_env
440 : TYPE(virial_type), POINTER :: virial
441 :
442 16149 : NULLIFY (dft_control)
443 16149 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control, natom=natom)
444 16149 : nimages = dft_control%nimages
445 16149 : IF (PRESENT(atcore)) THEN
446 66 : CPASSERT(SIZE(atcore) >= natom)
447 : END IF
448 :
449 : ! check whether a gauge transformed version of the non-local potential part has to be used
450 16149 : my_gt_nl = .FALSE.
451 16149 : IF (qs_env%run_rtp) THEN
452 1244 : CPASSERT(ASSOCIATED(dft_control%rtp_control))
453 1244 : IF (dft_control%rtp_control%velocity_gauge) THEN
454 54 : my_gt_nl = dft_control%rtp_control%nl_gauge_transform
455 : END IF
456 : END IF
457 :
458 : ! prepare for k-points
459 16149 : NULLIFY (cell_to_index)
460 16149 : IF (nimages > 1) THEN
461 292 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
462 292 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
463 292 : dft_control%qs_control%do_ppl_method = do_ppl_analytic
464 : END IF
465 : ! force analytic ppl calculation for GAPW methods
466 16149 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
467 2110 : dft_control%qs_control%do_ppl_method = do_ppl_analytic
468 : END IF
469 :
470 : ! force
471 16149 : NULLIFY (force)
472 16149 : IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force)
473 : ! virial
474 16149 : CALL get_qs_env(qs_env=qs_env, virial=virial)
475 16149 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
476 :
477 16149 : NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
478 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
479 16149 : particle_set=particle_set)
480 :
481 16149 : NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
482 : CALL get_qs_env(qs_env=qs_env, &
483 : sab_orb=sab_orb, &
484 : sac_ae=sac_ae, &
485 : sac_ppl=sac_ppl, &
486 16149 : sap_ppnl=sap_ppnl)
487 :
488 : ! *** compute the nuclear attraction contribution to the core hamiltonian ***
489 16149 : all_present = ASSOCIATED(sac_ae)
490 16149 : IF (all_present) THEN
491 : CALL build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
492 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
493 1956 : nimages, cell_to_index, atcore=atcore)
494 : END IF
495 :
496 : ! *** compute the ppl contribution to the core hamiltonian ***
497 16149 : ppl_present = ASSOCIATED(sac_ppl)
498 16149 : IF (ppl_present) THEN
499 15219 : IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
500 15195 : IF (dft_control%qs_control%lrigpw) THEN
501 80 : CALL get_qs_env(qs_env, lri_env=lri_env)
502 80 : IF (lri_env%ppl_ri) THEN
503 4 : IF (lri_env%exact_1c_terms) THEN
504 0 : CPABORT("not implemented")
505 : END IF
506 : ELSE
507 : CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
508 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
509 152 : nimages, cell_to_index, "ORB", atcore=atcore)
510 : END IF
511 : ELSE
512 : CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
513 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
514 30168 : nimages, cell_to_index, "ORB", atcore=atcore)
515 : END IF
516 : END IF
517 : END IF
518 :
519 : ! *** compute the ppnl contribution to the core hamiltonian ***
520 16149 : eps_ppnl = dft_control%qs_control%eps_ppnl
521 16149 : ppnl_present = ASSOCIATED(sap_ppnl)
522 16149 : IF (ppnl_present) THEN
523 12348 : IF (.NOT. my_gt_nl) THEN
524 : CALL build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
525 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
526 24562 : nimages, cell_to_index, "ORB", atcore=atcore)
527 : END IF
528 : END IF
529 :
530 16149 : END SUBROUTINE core_matrices
531 :
532 : ! **************************************************************************************************
533 : !> \brief Adds atomic blocks of relativistic correction for the kinetic energy
534 : !> \param matrix_h ...
535 : !> \param atomic_kind_set ...
536 : !> \param qs_kind_set ...
537 : ! **************************************************************************************************
538 58 : SUBROUTINE build_atomic_relmat(matrix_h, atomic_kind_set, qs_kind_set)
539 : TYPE(dbcsr_type), POINTER :: matrix_h
540 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
541 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
542 :
543 : INTEGER :: blk, iatom, ikind, jatom
544 58 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
545 58 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hblock, reltmat
546 : TYPE(dbcsr_iterator_type) :: iter
547 :
548 58 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
549 :
550 58 : CALL dbcsr_iterator_start(iter, matrix_h)
551 221 : DO WHILE (dbcsr_iterator_blocks_left(iter))
552 163 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, hblock, blk)
553 221 : IF (iatom == jatom) THEN
554 83 : ikind = kind_of(iatom)
555 83 : CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat)
556 192683 : IF (ASSOCIATED(reltmat)) hblock = hblock + reltmat
557 : END IF
558 : END DO
559 58 : CALL dbcsr_iterator_stop(iter)
560 :
561 116 : END SUBROUTINE build_atomic_relmat
562 :
563 : ! **************************************************************************************************
564 : !> \brief Possibly prints matrices after the construction of the Core
565 : !> Hamiltonian Matrix
566 : !> \param qs_env ...
567 : !> \param calculate_forces ...
568 : ! **************************************************************************************************
569 32166 : SUBROUTINE dump_info_core_hamiltonian(qs_env, calculate_forces)
570 : TYPE(qs_environment_type), POINTER :: qs_env
571 : LOGICAL, INTENT(IN) :: calculate_forces
572 :
573 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_info_core_hamiltonian'
574 :
575 : INTEGER :: after, handle, i, ic, iw, output_unit
576 : LOGICAL :: omit_headers
577 : TYPE(cp_logger_type), POINTER :: logger
578 16083 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_v
579 16083 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_h, matrixkp_s, matrixkp_t
580 : TYPE(mp_para_env_type), POINTER :: para_env
581 :
582 16083 : CALL timeset(routineN, handle)
583 :
584 16083 : NULLIFY (logger, matrix_v, para_env)
585 16083 : logger => cp_get_default_logger()
586 16083 : CALL get_qs_env(qs_env, para_env=para_env)
587 :
588 : ! Print the distribution of the overlap matrix blocks
589 : ! this duplicates causes duplicate printing at the force calc
590 16083 : IF (.NOT. calculate_forces) THEN
591 10412 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
592 : qs_env%input, "PRINT%DISTRIBUTION"), cp_p_file)) THEN
593 : output_unit = cp_print_key_unit_nr(logger, qs_env%input, "PRINT%DISTRIBUTION", &
594 92 : extension=".distribution")
595 92 : CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
596 92 : CALL cp_dbcsr_write_matrix_dist(matrixkp_s(1, 1)%matrix, output_unit, para_env)
597 92 : CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, "PRINT%DISTRIBUTION")
598 : END IF
599 : END IF
600 :
601 16083 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
602 : ! Print the overlap integral matrix, if requested
603 16083 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
604 : qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
605 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
606 4 : extension=".Log")
607 4 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
608 4 : after = MIN(MAX(after, 1), 16)
609 4 : CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
610 4 : IF (ASSOCIATED(matrixkp_s)) THEN
611 8 : DO ic = 1, SIZE(matrixkp_s, 2)
612 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(1, ic)%matrix, 4, after, qs_env, para_env, &
613 8 : output_unit=iw, omit_headers=omit_headers)
614 : END DO
615 4 : IF (BTEST(cp_print_key_should_output(logger%iter_info, qs_env%input, &
616 : "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
617 8 : DO ic = 1, SIZE(matrixkp_s, 2)
618 20 : DO i = 2, SIZE(matrixkp_s, 1)
619 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(i, ic)%matrix, 4, after, qs_env, para_env, &
620 16 : output_unit=iw, omit_headers=omit_headers)
621 : END DO
622 : END DO
623 : END IF
624 : END IF
625 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
626 4 : "DFT%PRINT%AO_MATRICES/OVERLAP")
627 : END IF
628 :
629 : ! Print the kinetic energy integral matrix, if requested
630 16083 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
631 : qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY"), cp_p_file)) THEN
632 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY", &
633 92 : extension=".Log")
634 92 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
635 92 : after = MIN(MAX(after, 1), 16)
636 92 : CALL get_qs_env(qs_env, kinetic_kp=matrixkp_t)
637 92 : IF (ASSOCIATED(matrixkp_t)) THEN
638 184 : DO ic = 1, SIZE(matrixkp_t, 2)
639 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_t(1, ic)%matrix, 4, after, qs_env, para_env, &
640 184 : output_unit=iw, omit_headers=omit_headers)
641 : END DO
642 : END IF
643 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
644 92 : "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY")
645 : END IF
646 :
647 : ! Print the potential energy matrix, if requested
648 16083 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
649 : qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY"), cp_p_file)) THEN
650 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY", &
651 92 : extension=".Log")
652 92 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
653 92 : after = MIN(MAX(after, 1), 16)
654 92 : CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h, kinetic_kp=matrixkp_t)
655 92 : IF (ASSOCIATED(matrixkp_h)) THEN
656 92 : IF (SIZE(matrixkp_h, 2) == 1) THEN
657 92 : CALL dbcsr_allocate_matrix_set(matrix_v, 1)
658 92 : ALLOCATE (matrix_v(1)%matrix)
659 92 : CALL dbcsr_copy(matrix_v(1)%matrix, matrixkp_h(1, 1)%matrix, name="POTENTIAL ENERGY MATRIX")
660 : CALL dbcsr_add(matrix_v(1)%matrix, matrixkp_t(1, 1)%matrix, &
661 92 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
662 : CALL cp_dbcsr_write_sparse_matrix(matrix_v(1)%matrix, 4, after, qs_env, &
663 92 : para_env, output_unit=iw, omit_headers=omit_headers)
664 92 : CALL dbcsr_deallocate_matrix_set(matrix_v)
665 : ELSE
666 0 : CPWARN("Printing of potential energy matrix not implemented for k-points")
667 : END IF
668 : END IF
669 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
670 92 : "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY")
671 : END IF
672 :
673 : ! Print the core Hamiltonian matrix, if requested
674 16083 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
675 : qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
676 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
677 92 : extension=".Log")
678 92 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
679 92 : after = MIN(MAX(after, 1), 16)
680 92 : CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h)
681 92 : IF (ASSOCIATED(matrixkp_h)) THEN
682 184 : DO ic = 1, SIZE(matrixkp_h, 2)
683 : CALL cp_dbcsr_write_sparse_matrix(matrixkp_h(1, ic)%matrix, 4, after, qs_env, para_env, &
684 184 : output_unit=iw, omit_headers=omit_headers)
685 : END DO
686 : END IF
687 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
688 92 : "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
689 : END IF
690 :
691 16083 : CALL timestop(handle)
692 :
693 16083 : END SUBROUTINE dump_info_core_hamiltonian
694 :
695 : ! **************************************************************************************************
696 : !> \brief (Re-)allocate matrix_h based on the template (typically the overlap matrix)
697 : !> \param qs_env ...
698 : !> \param template ...
699 : !> \param is_complex ...
700 : ! **************************************************************************************************
701 16083 : SUBROUTINE qs_matrix_h_allocate(qs_env, template, is_complex)
702 : TYPE(qs_environment_type) :: qs_env
703 : TYPE(dbcsr_type), INTENT(in) :: template
704 : LOGICAL, INTENT(in) :: is_complex
705 :
706 : CHARACTER(LEN=default_string_length) :: headline
707 : INTEGER :: img, nimages
708 16083 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im
709 : TYPE(dft_control_type), POINTER :: dft_control
710 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
711 16083 : POINTER :: sab_orb
712 : TYPE(qs_ks_env_type), POINTER :: ks_env
713 :
714 16083 : NULLIFY (matrix_h, matrix_h_im, sab_orb, dft_control, ks_env)
715 : CALL get_qs_env(qs_env=qs_env, &
716 : matrix_h_kp=matrix_h, &
717 : matrix_h_im_kp=matrix_h_im, &
718 : sab_orb=sab_orb, &
719 : dft_control=dft_control, &
720 16083 : ks_env=ks_env)
721 :
722 16083 : nimages = dft_control%nimages
723 16083 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimages)
724 16083 : headline = "CORE HAMILTONIAN MATRIX"
725 64758 : DO img = 1, nimages
726 48675 : ALLOCATE (matrix_h(1, img)%matrix)
727 48675 : CALL dbcsr_create(matrix_h(1, img)%matrix, name=TRIM(headline), template=template)
728 48675 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
729 64758 : CALL dbcsr_set(matrix_h(1, img)%matrix, 0.0_dp)
730 : END DO
731 16083 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
732 :
733 16083 : IF (is_complex) THEN
734 352 : headline = "IMAGINARY PART OF CORE HAMILTONIAN MATRIX"
735 352 : CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
736 704 : DO img = 1, nimages
737 352 : ALLOCATE (matrix_h_im(1, img)%matrix)
738 : CALL dbcsr_create(matrix_h_im(1, img)%matrix, name=TRIM(headline), template=template, &
739 352 : matrix_type=dbcsr_type_antisymmetric)
740 352 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, img)%matrix, sab_orb)
741 704 : CALL dbcsr_set(matrix_h_im(1, img)%matrix, 0.0_dp)
742 : END DO
743 352 : CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
744 : END IF
745 :
746 16083 : END SUBROUTINE qs_matrix_h_allocate
747 :
748 : ! **************************************************************************************************
749 : !> \brief (Re-)allocates matrix_h_im from matrix_h
750 : !> \param qs_env ...
751 : ! **************************************************************************************************
752 8 : SUBROUTINE qs_matrix_h_allocate_imag_from_real(qs_env)
753 : TYPE(qs_environment_type) :: qs_env
754 :
755 : CHARACTER(LEN=default_string_length) :: headline
756 : INTEGER :: image, nimages
757 8 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im
758 : TYPE(dbcsr_type), POINTER :: template
759 : TYPE(dft_control_type), POINTER :: dft_control
760 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
761 8 : POINTER :: sab_orb
762 : TYPE(qs_ks_env_type), POINTER :: ks_env
763 :
764 8 : NULLIFY (matrix_h_im, matrix_h, dft_control, template, sab_orb, ks_env)
765 :
766 : CALL get_qs_env(qs_env, &
767 : matrix_h_im_kp=matrix_h_im, &
768 : matrix_h_kp=matrix_h, &
769 : dft_control=dft_control, &
770 : sab_orb=sab_orb, &
771 8 : ks_env=ks_env)
772 :
773 8 : nimages = dft_control%nimages
774 :
775 8 : CPASSERT(nimages .EQ. SIZE(matrix_h, 2))
776 :
777 8 : CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
778 :
779 16 : DO image = 1, nimages
780 8 : headline = "IMAGINARY CORE HAMILTONIAN MATRIX"
781 8 : ALLOCATE (matrix_h_im(1, image)%matrix)
782 8 : template => matrix_h(1, image)%matrix ! base on real part, but anti-symmetric
783 : CALL dbcsr_create(matrix=matrix_h_im(1, image)%matrix, template=template, &
784 8 : name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric, nze=0)
785 8 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, image)%matrix, sab_orb)
786 16 : CALL dbcsr_set(matrix_h_im(1, image)%matrix, 0.0_dp)
787 : END DO
788 8 : CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
789 :
790 8 : END SUBROUTINE qs_matrix_h_allocate_imag_from_real
791 :
792 : END MODULE qs_core_hamiltonian
|
