Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief interface to tblite
10 : !> \author JVP
11 : !> \history creation 09.2024
12 : ! **************************************************************************************************
13 :
14 : MODULE tblite_interface
15 :
16 : #if defined(__TBLITE)
17 : USE mctc_env, ONLY: error_type
18 : USE mctc_io, ONLY: new
19 : USE mctc_io_symbols, ONLY: symbol_to_number
20 : USE tblite_basis_type, ONLY: get_cutoff
21 : USE tblite_container_cache, ONLY: container_cache
22 : USE tblite_container_type, ONLY: container_type
23 : USE tblite_data_spin, ONLY: get_spin_constant
24 : USE tblite_integral_multipole, ONLY: multipole_cgto, multipole_grad_cgto, msao
25 : USE tblite_scf, ONLY: get_mixer_dimension
26 : USE tblite_scf_info, ONLY: scf_info, atom_resolved, shell_resolved
27 : USE tblite_scf_potential, ONLY: new_potential
28 : USE tblite_spin, ONLY: spin_polarization, new_spin_polarization
29 : USE tblite_wavefunction_type, ONLY: new_wavefunction
30 : USE tblite_xtb_calculator, ONLY: new_xtb_calculator
31 : USE tblite_xtb_gfn1, ONLY: new_gfn1_calculator
32 : USE tblite_xtb_gfn2, ONLY: new_gfn2_calculator
33 : USE tblite_xtb_h0, ONLY: get_selfenergy, get_occupation, tb_hamiltonian
34 : USE tblite_xtb_ipea1, ONLY: new_ipea1_calculator
35 : #endif
36 : USE ai_contraction, ONLY: block_add, contraction
37 : USE ai_overlap, ONLY: overlap_ab
38 : USE atomic_kind_types, ONLY: atomic_kind_type, get_atomic_kind, get_atomic_kind_set
39 : USE atprop_types, ONLY: atprop_type
40 : USE basis_set_types, ONLY: gto_basis_set_type, gto_basis_set_p_type, allocate_gto_basis_set, &
41 : process_gto_basis
42 : USE block_p_types, ONLY: block_p_type
43 : USE cell_types, ONLY: cell_type, get_cell
44 : USE cp_blacs_env, ONLY: cp_blacs_env_type
45 : USE cp_control_types, ONLY: dft_control_type
46 : USE cp_dbcsr_api, ONLY: dbcsr_p_type, dbcsr_create, dbcsr_add, dbcsr_copy, &
47 : dbcsr_get_block_p, dbcsr_finalize, dbcsr_iterator_type, &
48 : dbcsr_iterator_blocks_left, dbcsr_iterator_start, dbcsr_iterator_stop, &
49 : dbcsr_iterator_next_block
50 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
51 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
52 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set, dbcsr_deallocate_matrix_set
53 : USE cp_log_handling, ONLY: cp_get_default_logger, cp_logger_type, cp_logger_get_default_io_unit
54 : USE cp_output_handling, ONLY: cp_print_key_should_output, cp_print_key_unit_nr, &
55 : cp_print_key_finished_output
56 : USE input_constants, ONLY: gfn1xtb, gfn2xtb, ipea1xtb, tblite_scc_mixer_auto, tblite_scc_mixer_cp2k, &
57 : tblite_scc_mixer_none, tblite_scc_mixer_tblite, tblite_solver_gvd, &
58 : tblite_solver_gvr
59 : USE input_section_types, ONLY: section_vals_val_get, section_vals_get_subs_vals, section_vals_type
60 : USE kinds, ONLY: default_string_length, dp, int_8
61 : USE kpoint_types, ONLY: get_kpoint_info, kpoint_type
62 : USE memory_utilities, ONLY: reallocate
63 : USE message_passing, ONLY: mp_para_env_type
64 : USE orbital_pointers, ONLY: ncoset
65 : USE particle_types, ONLY: particle_type
66 : USE qs_charge_mixing, ONLY: charge_mixing, tblite_scc_error_on_cp2k_scale, tblite_scc_pconv
67 : USE qs_condnum, ONLY: overlap_condnum
68 : USE qs_density_mixing_types, ONLY: modified_broyden_mixing_nr
69 : USE qs_energy_types, ONLY: qs_energy_type
70 : USE qs_environment_types, ONLY: get_qs_env, qs_environment_type
71 : USE qs_force_types, ONLY: qs_force_type
72 : USE qs_integral_utils, ONLY: basis_set_list_setup, get_memory_usage
73 : USE qs_kind_types, ONLY: get_qs_kind, qs_kind_type, get_qs_kind_set
74 : USE qs_ks_types, ONLY: get_ks_env, qs_ks_env_type, set_ks_env
75 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
76 : USE qs_overlap, ONLY: build_overlap_matrix, create_sab_matrix
77 : USE qs_rho_types, ONLY: qs_rho_get, qs_rho_type
78 : USE qs_scf_types, ONLY: qs_scf_env_type
79 : USE scf_control_types, ONLY: scf_control_type
80 : USE string_utilities, ONLY: integer_to_string
81 : #if defined(__TBLITE)
82 : USE tblite_scc_mixer, ONLY: new_cp2k_tblite_mixer
83 : #endif
84 : USE tblite_types, ONLY: tblite_type, allocate_tblite_type
85 : USE virial_types, ONLY: virial_type
86 : USE xtb_types, ONLY: get_xtb_atom_param, xtb_atom_type
87 :
88 : !$ USE OMP_LIB, ONLY: omp_destroy_lock, omp_init_lock, omp_set_lock, &
89 : !$ omp_unset_lock, omp_lock_kind
90 :
91 : #include "../base/base_uses.f90"
92 : IMPLICIT NONE
93 :
94 : PRIVATE
95 :
96 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tblite_interface'
97 :
98 : INTEGER, PARAMETER :: dip_n = 3
99 : INTEGER, PARAMETER :: quad_n = 6
100 : REAL(KIND=dp), PARAMETER :: same_atom = 0.00001_dp
101 :
102 : PUBLIC :: tb_set_calculator, tb_init_geometry, tb_init_wf
103 : PUBLIC :: tb_get_basis, build_tblite_matrices
104 : PUBLIC :: tb_get_energy, tb_update_charges, tb_ham_add_coulomb
105 : PUBLIC :: tb_native_scc_mixer_active
106 : PUBLIC :: tb_scf_mixer_error
107 : PUBLIC :: tb_get_multipole
108 : PUBLIC :: tb_derive_dH_off
109 :
110 : CONTAINS
111 :
112 : #if defined(__TBLITE)
113 : ! **************************************************************************************************
114 : !> \brief Project a tblite charge/magnetization quantity to a CP2K spin channel.
115 : !> \param values ...
116 : !> \param ispin ...
117 : !> \return ...
118 : ! **************************************************************************************************
119 1985072 : PURE FUNCTION tb_spin_project(values, ispin) RESULT(value)
120 :
121 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: values
122 : INTEGER, INTENT(IN) :: ispin
123 : REAL(KIND=dp) :: value
124 :
125 1985072 : value = values(1)
126 1985072 : IF (SIZE(values) > 1) THEN
127 256620 : SELECT CASE (ispin)
128 : CASE (1)
129 85540 : value = values(1) + values(2)
130 : CASE (2)
131 85540 : value = values(1) - values(2)
132 : CASE DEFAULT
133 171080 : value = values(1)
134 : END SELECT
135 : END IF
136 :
137 1985072 : END FUNCTION tb_spin_project
138 :
139 : ! **************************************************************************************************
140 : !> \brief Preserve spin-resolved density before force assembly combines spin channels.
141 : !> \param tb ...
142 : !> \param matrix_p ...
143 : ! **************************************************************************************************
144 22 : SUBROUTINE tb_store_density_ref(tb, matrix_p)
145 :
146 : TYPE(tblite_type), POINTER :: tb
147 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
148 :
149 : INTEGER :: img, ispin, nimg, nspin
150 :
151 22 : nspin = SIZE(matrix_p, 1)
152 22 : nimg = SIZE(matrix_p, 2)
153 22 : IF (ASSOCIATED(tb%rho_ao_kp_ref)) THEN
154 2 : IF (SIZE(tb%rho_ao_kp_ref, 1) /= nspin .OR. SIZE(tb%rho_ao_kp_ref, 2) /= nimg) THEN
155 0 : CALL dbcsr_deallocate_matrix_set(tb%rho_ao_kp_ref)
156 : END IF
157 : END IF
158 22 : IF (.NOT. ASSOCIATED(tb%rho_ao_kp_ref)) THEN
159 20 : CALL dbcsr_allocate_matrix_set(tb%rho_ao_kp_ref, nspin, nimg)
160 94 : DO img = 1, nimg
161 242 : DO ispin = 1, nspin
162 222 : ALLOCATE (tb%rho_ao_kp_ref(ispin, img)%matrix)
163 : END DO
164 : END DO
165 : END IF
166 152 : DO img = 1, nimg
167 412 : DO ispin = 1, nspin
168 390 : CALL dbcsr_copy(tb%rho_ao_kp_ref(ispin, img)%matrix, matrix_p(ispin, img)%matrix)
169 : END DO
170 : END DO
171 :
172 22 : END SUBROUTINE tb_store_density_ref
173 : #endif
174 :
175 : ! **************************************************************************************************
176 : !> \brief intialize geometry objects ...
177 : !> \param qs_env ...
178 : !> \param tb ...
179 : ! **************************************************************************************************
180 198 : SUBROUTINE tb_init_geometry(qs_env, tb)
181 :
182 : TYPE(qs_environment_type), POINTER :: qs_env
183 : TYPE(tblite_type), POINTER :: tb
184 :
185 : #if defined(__TBLITE)
186 :
187 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tblite_init_geometry'
188 :
189 : TYPE(cell_type), POINTER :: cell
190 198 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
191 198 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
192 : INTEGER :: iatom, natom
193 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xyz
194 : INTEGER :: handle, ikind
195 : INTEGER, DIMENSION(3) :: periodic
196 : LOGICAL, DIMENSION(3) :: lperiod
197 :
198 198 : CALL timeset(routineN, handle)
199 :
200 : !get info from environment vaiarable
201 198 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell, qs_kind_set=qs_kind_set)
202 :
203 : !get information about particles
204 198 : natom = SIZE(particle_set)
205 594 : ALLOCATE (xyz(3, natom))
206 198 : CALL allocate_tblite_type(tb)
207 594 : ALLOCATE (tb%el_num(natom))
208 1182 : tb%el_num = -9
209 1182 : DO iatom = 1, natom
210 3936 : xyz(:, iatom) = particle_set(iatom)%r(:)
211 984 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
212 984 : CALL get_qs_kind(qs_kind_set(ikind), zatom=tb%el_num(iatom))
213 2166 : IF (tb%el_num(iatom) < 1 .OR. tb%el_num(iatom) > 85) THEN
214 0 : CPABORT("only elements 1-85 are supported by tblite")
215 : END IF
216 : END DO
217 :
218 : !get information about cell / lattice
219 198 : CALL get_cell(cell=cell, periodic=periodic)
220 198 : lperiod(1) = periodic(1) == 1
221 198 : lperiod(2) = periodic(2) == 1
222 198 : lperiod(3) = periodic(3) == 1
223 :
224 : !prepare for the call to the dispersion function
225 198 : CALL new(tb%mol, tb%el_num, xyz, lattice=cell%hmat, periodic=lperiod)
226 :
227 198 : DEALLOCATE (xyz)
228 :
229 198 : CALL timestop(handle)
230 :
231 : #else
232 : MARK_USED(qs_env)
233 : MARK_USED(tb)
234 : CPABORT("Built without TBLITE")
235 : #endif
236 :
237 198 : END SUBROUTINE tb_init_geometry
238 :
239 : ! **************************************************************************************************
240 : !> \brief updating coordinates...
241 : !> \param qs_env ...
242 : !> \param tb ...
243 : ! **************************************************************************************************
244 2770 : SUBROUTINE tb_update_geometry(qs_env, tb)
245 :
246 : TYPE(qs_environment_type) :: qs_env
247 : TYPE(tblite_type) :: tb
248 :
249 : #if defined(__TBLITE)
250 :
251 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tblite_update_geometry'
252 :
253 : TYPE(cell_type), POINTER :: cell
254 2770 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
255 : INTEGER :: iatom, natom
256 : INTEGER :: handle
257 :
258 2770 : CALL timeset(routineN, handle)
259 :
260 : !get info from environment vaiarable
261 2770 : NULLIFY (cell)
262 2770 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell)
263 :
264 : !get information about particles
265 2770 : natom = SIZE(particle_set)
266 14982 : DO iatom = 1, natom
267 51618 : tb%mol%xyz(:, iatom) = particle_set(iatom)%r(:)
268 : END DO
269 36010 : tb%mol%lattice(:, :) = cell%hmat
270 :
271 2770 : CALL timestop(handle)
272 :
273 : #else
274 : MARK_USED(qs_env)
275 : MARK_USED(tb)
276 : CPABORT("Built without TBLITE")
277 : #endif
278 :
279 2770 : END SUBROUTINE tb_update_geometry
280 :
281 : ! **************************************************************************************************
282 : !> \brief initialize wavefunction ...
283 : !> \param tb ...
284 : !> \param dft_control ...
285 : ! **************************************************************************************************
286 198 : SUBROUTINE tb_init_wf(tb, dft_control)
287 :
288 : TYPE(tblite_type), POINTER :: tb
289 : TYPE(dft_control_type), POINTER :: dft_control
290 :
291 : #if defined(__TBLITE)
292 :
293 : INTEGER :: nSpin
294 :
295 : TYPE(scf_info) :: info
296 :
297 198 : nSpin = dft_control%nspins
298 0 : IF (nSpin /= 1 .AND. nSpin /= 2) CPABORT("tblite supports only one or two spin channels")
299 :
300 198 : tb%mol%charge = dft_control%charge
301 198 : tb%mol%uhf = MAX(0, dft_control%multiplicity - 1)
302 198 : IF (nSpin == 2) CALL tb_add_spin_polarization(tb)
303 :
304 198 : info = tb%calc%variable_info()
305 198 : IF (info%charge > shell_resolved) CPABORT("tblite: no support for orbital resolved charge")
306 198 : IF (info%dipole > atom_resolved) CPABORT("tblite: no support for shell resolved dipole moment")
307 198 : IF (info%quadrupole > atom_resolved) THEN
308 0 : CPABORT("tblite: no support shell resolved quadrupole moment")
309 : END IF
310 :
311 198 : CALL new_wavefunction(tb%wfn, tb%mol%nat, tb%calc%bas%nsh, tb%calc%bas%nao, nSpin, 0.0_dp)
312 198 : CALL get_occupation(tb%mol, tb%calc%bas, tb%calc%h0, tb%wfn%nocc, tb%wfn%n0at, tb%wfn%n0sh)
313 198 : CALL tb_reset_mixer(tb)
314 :
315 198 : CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
316 :
317 : !allocate quantities later required
318 1386 : ALLOCATE (tb%e_hal(tb%mol%nat), tb%e_rep(tb%mol%nat), tb%e_disp(tb%mol%nat))
319 1386 : ALLOCATE (tb%e_scd(tb%mol%nat), tb%e_es(tb%mol%nat), tb%e_int(tb%mol%nat))
320 594 : ALLOCATE (tb%selfenergy(tb%calc%bas%nsh))
321 594 : IF (ALLOCATED(tb%calc%ncoord)) ALLOCATE (tb%cn(tb%mol%nat))
322 :
323 : #else
324 : MARK_USED(tb)
325 : MARK_USED(dft_control)
326 : CPABORT("Built without TBLITE")
327 : #endif
328 :
329 198 : END SUBROUTINE tb_init_wf
330 :
331 : #if defined(__TBLITE)
332 : ! **************************************************************************************************
333 : !> \brief Add tblite's on-site spin-polarization interaction.
334 : !> \param tb ...
335 : ! **************************************************************************************************
336 30 : SUBROUTINE tb_add_spin_polarization(tb)
337 :
338 : TYPE(tblite_type), POINTER :: tb
339 :
340 30 : CLASS(container_type), ALLOCATABLE :: cont
341 30 : TYPE(spin_polarization), ALLOCATABLE :: spin
342 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wll
343 :
344 30 : ALLOCATE (spin)
345 30 : CALL tb_get_spin_constants(tb, wll)
346 30 : CALL new_spin_polarization(spin, tb%mol, wll, tb%calc%bas%nsh_id)
347 30 : CALL MOVE_ALLOC(spin, cont)
348 30 : CALL tb%calc%push_back(cont)
349 :
350 30 : END SUBROUTINE tb_add_spin_polarization
351 :
352 : ! **************************************************************************************************
353 : !> \brief Build tblite spin constants for the current basis.
354 : !> \param tb ...
355 : !> \param wll ...
356 : ! **************************************************************************************************
357 30 : SUBROUTINE tb_get_spin_constants(tb, wll)
358 :
359 : TYPE(tblite_type), POINTER :: tb
360 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
361 : INTENT(OUT) :: wll
362 :
363 : INTEGER :: il, ish, izp, jl, jsh, max_shell
364 :
365 68 : max_shell = MAXVAL(tb%calc%bas%nsh_id)
366 150 : ALLOCATE (wll(max_shell, max_shell, tb%mol%nid))
367 30 : wll = 0.0_dp
368 68 : DO izp = 1, tb%mol%nid
369 144 : DO ish = 1, tb%calc%bas%nsh_id(izp)
370 76 : il = tb%calc%bas%cgto(ish, izp)%ang
371 270 : DO jsh = 1, tb%calc%bas%nsh_id(izp)
372 156 : jl = tb%calc%bas%cgto(jsh, izp)%ang
373 232 : wll(jsh, ish, izp) = get_spin_constant(jl, il, tb%mol%num(izp))
374 : END DO
375 : END DO
376 : END DO
377 :
378 30 : END SUBROUTINE tb_get_spin_constants
379 : #endif
380 :
381 : ! **************************************************************************************************
382 : !> \brief Reset tblite's internal SCC mixer for a new CP2K SCF cycle.
383 : !> \param tb ...
384 : ! **************************************************************************************************
385 2656 : SUBROUTINE tb_reset_mixer(tb)
386 :
387 : TYPE(tblite_type), POINTER :: tb
388 :
389 : #if defined(__TBLITE)
390 :
391 : TYPE(scf_info) :: info
392 :
393 2656 : info = tb%calc%variable_info()
394 5114 : IF (ALLOCATED(tb%mixer)) DEALLOCATE (tb%mixer)
395 : CALL new_cp2k_tblite_mixer(tb%mixer, tb%mixer_memory, &
396 : tb%wfn%nspin*get_mixer_dimension(tb%mol, tb%calc%bas, info), &
397 : tb%mixer_damping, tb%mixer_omega0, tb%mixer_min_weight, &
398 2656 : tb%mixer_max_weight, tb%mixer_weight_factor)
399 :
400 : #else
401 : MARK_USED(tb)
402 : CPABORT("Built without TBLITE")
403 : #endif
404 :
405 2656 : END SUBROUTINE tb_reset_mixer
406 :
407 : ! **************************************************************************************************
408 : !> \brief Configure tblite's internal SCC mixer from CP2K input.
409 : !> \param tb ...
410 : !> \param iterations ...
411 : !> \param memory ...
412 : !> \param damping ...
413 : !> \param omega0 ...
414 : !> \param min_weight ...
415 : !> \param max_weight ...
416 : !> \param weight_factor ...
417 : !> \param solver ...
418 : ! **************************************************************************************************
419 2458 : SUBROUTINE tb_configure_mixer(tb, iterations, memory, damping, omega0, min_weight, max_weight, &
420 : weight_factor, solver)
421 :
422 : TYPE(tblite_type), POINTER :: tb
423 : INTEGER, INTENT(IN) :: iterations, memory, solver
424 : REAL(KIND=dp), INTENT(IN) :: damping, max_weight, min_weight, omega0, &
425 : weight_factor
426 :
427 : #if defined(__TBLITE)
428 :
429 2458 : IF (iterations < 1) CPABORT("tblite SCC mixer ITERATIONS must be positive")
430 2458 : IF (memory < 1) CPABORT("tblite SCC mixer MEMORY must be positive")
431 2458 : IF (damping <= 0.0_dp) CPABORT("tblite SCC mixer damping must be positive")
432 2458 : IF (omega0 <= 0.0_dp) CPABORT("tblite SCC mixer OMEGA0 must be positive")
433 2458 : IF (min_weight <= 0.0_dp) CPABORT("tblite SCC mixer MIN_WEIGHT must be positive")
434 2458 : IF (max_weight <= 0.0_dp) CPABORT("tblite SCC mixer MAX_WEIGHT must be positive")
435 2458 : IF (max_weight < min_weight) THEN
436 0 : CPABORT("tblite SCC mixer MAX_WEIGHT must not be smaller than MIN_WEIGHT")
437 : END IF
438 2458 : IF (weight_factor <= 0.0_dp) CPABORT("tblite SCC mixer WEIGHT_FACTOR must be positive")
439 2458 : SELECT CASE (solver)
440 : CASE (tblite_solver_gvd, tblite_solver_gvr)
441 : CASE DEFAULT
442 2458 : CPABORT("Unknown tblite SCC mixer SOLVER")
443 : END SELECT
444 :
445 2458 : tb%calc%max_iter = iterations
446 2458 : tb%mixer_memory = memory
447 2458 : tb%mixer_damping = damping
448 2458 : tb%calc%mixer_input%damping = damping
449 2458 : tb%mixer_omega0 = omega0
450 2458 : tb%mixer_min_weight = min_weight
451 2458 : tb%mixer_max_weight = max_weight
452 2458 : tb%mixer_weight_factor = weight_factor
453 :
454 : #else
455 : MARK_USED(tb)
456 : MARK_USED(iterations)
457 : MARK_USED(memory)
458 : MARK_USED(damping)
459 : MARK_USED(omega0)
460 : MARK_USED(min_weight)
461 : MARK_USED(max_weight)
462 : MARK_USED(weight_factor)
463 : MARK_USED(solver)
464 : CPABORT("Built without TBLITE")
465 : #endif
466 :
467 2458 : END SUBROUTINE tb_configure_mixer
468 :
469 : ! **************************************************************************************************
470 : !> \brief Return whether the tblite native SCC mixer is active for this run.
471 : !> \param dft_control ...
472 : !> \return ...
473 : ! **************************************************************************************************
474 94574 : FUNCTION tb_native_scc_mixer_active(dft_control) RESULT(use_native_mixer)
475 :
476 : TYPE(dft_control_type), POINTER :: dft_control
477 : LOGICAL :: use_native_mixer
478 :
479 94574 : use_native_mixer = .FALSE.
480 94574 : IF (.NOT. ASSOCIATED(dft_control)) RETURN
481 94574 : IF (dft_control%qs_control%do_ls_scf) RETURN
482 :
483 94574 : SELECT CASE (dft_control%qs_control%xtb_control%tblite_scc_mixer)
484 : CASE (tblite_scc_mixer_auto)
485 : use_native_mixer = .TRUE.
486 : CASE (tblite_scc_mixer_tblite)
487 5076 : use_native_mixer = .TRUE.
488 : CASE (tblite_scc_mixer_cp2k, tblite_scc_mixer_none)
489 5076 : use_native_mixer = .FALSE.
490 : CASE DEFAULT
491 94574 : CPABORT("Unknown tblite SCC mixer")
492 : END SELECT
493 :
494 : END FUNCTION tb_native_scc_mixer_active
495 :
496 : ! **************************************************************************************************
497 : !> \brief Return the native tblite SCC mixer residual on the CP2K iter_delta scale.
498 : !> \param dft_control ...
499 : !> \param tb ...
500 : !> \param eps_scf CP2K reporting scale for the native residual.
501 : !> \return ...
502 : ! **************************************************************************************************
503 24766 : FUNCTION tb_scf_mixer_error(dft_control, tb, eps_scf) RESULT(mixer_error)
504 :
505 : TYPE(dft_control_type), POINTER :: dft_control
506 : TYPE(tblite_type), POINTER :: tb
507 : REAL(KIND=dp), INTENT(IN) :: eps_scf
508 : REAL(KIND=dp) :: mixer_error
509 :
510 : #if defined(__TBLITE)
511 : REAL(KIND=dp) :: raw_error
512 : #endif
513 :
514 24766 : mixer_error = 0.0_dp
515 :
516 : #if defined(__TBLITE)
517 24766 : IF (.NOT. ASSOCIATED(tb)) RETURN
518 24766 : IF (.NOT. tb_native_scc_mixer_active(dft_control)) RETURN
519 22608 : IF (ALLOCATED(tb%mixer)) THEN
520 22608 : raw_error = REAL(tb%mixer%get_error(), KIND=dp)
521 : mixer_error = tblite_scc_error_on_cp2k_scale( &
522 : raw_error, eps_scf, &
523 22608 : tblite_scc_pconv*dft_control%qs_control%xtb_control%tblite_accuracy)
524 : END IF
525 : #else
526 : MARK_USED(dft_control)
527 : MARK_USED(tb)
528 : MARK_USED(eps_scf)
529 : #endif
530 :
531 : END FUNCTION tb_scf_mixer_error
532 :
533 : ! **************************************************************************************************
534 : !> \brief ...
535 : !> \param tb ...
536 : !> \param typ ...
537 : !> \param accuracy ...
538 : !> \param param_file ...
539 : ! **************************************************************************************************
540 198 : SUBROUTINE tb_set_calculator(tb, typ, accuracy, param_file)
541 :
542 : TYPE(tblite_type), POINTER :: tb
543 : INTEGER :: typ
544 : REAL(KIND=dp), INTENT(IN) :: accuracy
545 : CHARACTER(LEN=*), INTENT(IN) :: param_file
546 :
547 : #if defined(__TBLITE)
548 :
549 198 : TYPE(error_type), ALLOCATABLE :: error
550 :
551 198 : IF (ALLOCATED(tb%param)) DEALLOCATE (tb%param)
552 198 : IF (LEN_TRIM(param_file) > 0) THEN
553 0 : ALLOCATE (tb%param)
554 0 : CALL tb%param%load(TRIM(param_file), error)
555 0 : IF (ALLOCATED(error)) CPABORT("Could not load tblite PARAM file: "//TRIM(param_file))
556 0 : CALL new_xtb_calculator(tb%calc, tb%mol, tb%param, error)
557 : ELSE
558 198 : SELECT CASE (typ)
559 : CASE default
560 0 : CPABORT("Unknown xtb type")
561 : CASE (gfn1xtb)
562 72 : CALL new_gfn1_calculator(tb%calc, tb%mol, error)
563 : CASE (gfn2xtb)
564 110 : CALL new_gfn2_calculator(tb%calc, tb%mol, error)
565 : CASE (ipea1xtb)
566 198 : CALL new_ipea1_calculator(tb%calc, tb%mol, error)
567 : END SELECT
568 : END IF
569 198 : IF (ALLOCATED(error)) CPABORT("tblite calculator setup failed")
570 :
571 198 : tb%accuracy = accuracy
572 :
573 : #else
574 : MARK_USED(tb)
575 : MARK_USED(typ)
576 : MARK_USED(accuracy)
577 : MARK_USED(param_file)
578 : CPABORT("Built without TBLITE")
579 : #endif
580 :
581 198 : END SUBROUTINE tb_set_calculator
582 :
583 : ! **************************************************************************************************
584 : !> \brief ...
585 : !> \param qs_env ...
586 : !> \param tb ...
587 : !> \param para_env ...
588 : ! **************************************************************************************************
589 2770 : SUBROUTINE tb_init_ham(qs_env, tb, para_env)
590 :
591 : TYPE(qs_environment_type) :: qs_env
592 : TYPE(tblite_type) :: tb
593 : TYPE(mp_para_env_type) :: para_env
594 :
595 : #if defined(__TBLITE)
596 :
597 2770 : TYPE(container_cache) :: hcache, rcache
598 :
599 14982 : tb%e_hal = 0.0_dp
600 14982 : tb%e_rep = 0.0_dp
601 14982 : tb%e_disp = 0.0_dp
602 14982 : tb%e_int = 0.0_dp
603 2770 : IF (ALLOCATED(tb%grad)) THEN
604 3754 : tb%grad = 0.0_dp
605 154 : CALL tb_zero_force(qs_env)
606 : END IF
607 36010 : tb%sigma = 0.0_dp
608 :
609 2770 : IF (ALLOCATED(tb%calc%halogen)) THEN
610 1388 : CALL tb%calc%halogen%update(tb%mol, hcache)
611 1388 : IF (ALLOCATED(tb%grad)) THEN
612 2662 : tb%grad = 0.0_dp
613 : CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal, &
614 86 : & tb%grad, tb%sigma)
615 86 : CALL tb_grad2force(qs_env, tb, para_env, 0)
616 : ELSE
617 1302 : CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal)
618 : END IF
619 : END IF
620 :
621 2770 : IF (ALLOCATED(tb%calc%repulsion)) THEN
622 2770 : CALL tb%calc%repulsion%update(tb%mol, rcache)
623 2770 : IF (ALLOCATED(tb%grad)) THEN
624 3754 : tb%grad = 0.0_dp
625 : CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep, &
626 154 : & tb%grad, tb%sigma)
627 154 : CALL tb_grad2force(qs_env, tb, para_env, 1)
628 : ELSE
629 2616 : CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep)
630 : END IF
631 : END IF
632 :
633 2770 : IF (ALLOCATED(tb%calc%dispersion)) THEN
634 2770 : CALL tb%calc%dispersion%update(tb%mol, tb%dcache)
635 2770 : IF (ALLOCATED(tb%grad)) THEN
636 3754 : tb%grad = 0.0_dp
637 : CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp, &
638 154 : & tb%grad, tb%sigma)
639 154 : CALL tb_grad2force(qs_env, tb, para_env, 2)
640 : ELSE
641 2616 : CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp)
642 : END IF
643 : END IF
644 :
645 2770 : IF (ALLOCATED(tb%calc%interactions)) THEN
646 356 : CALL tb%calc%interactions%update(tb%mol, tb%icache)
647 : END IF
648 :
649 2770 : CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
650 2770 : IF (ALLOCATED(tb%calc%coulomb)) THEN
651 2770 : CALL tb%calc%coulomb%update(tb%mol, tb%cache)
652 : END IF
653 :
654 2770 : IF (ALLOCATED(tb%grad)) THEN
655 154 : IF (ALLOCATED(tb%calc%ncoord)) THEN
656 154 : CALL tb%calc%ncoord%get_cn(tb%mol, tb%cn, tb%dcndr, tb%dcndL)
657 : END IF
658 : ELSE
659 2616 : IF (ALLOCATED(tb%calc%ncoord)) THEN
660 2616 : CALL tb%calc%ncoord%get_cn(tb%mol, tb%cn)
661 : END IF
662 : END IF
663 :
664 : CALL get_selfenergy(tb%calc%h0, tb%mol%id, tb%calc%bas%ish_at, &
665 2770 : tb%calc%bas%nsh_id, cn=tb%cn, selfenergy=tb%selfenergy, dsedcn=tb%dsedcn)
666 :
667 : #else
668 : MARK_USED(qs_env)
669 : MARK_USED(tb)
670 : MARK_USED(para_env)
671 : CPABORT("Built without TBLITE")
672 : #endif
673 :
674 2770 : END SUBROUTINE tb_init_ham
675 :
676 : ! **************************************************************************************************
677 : !> \brief ...
678 : !> \param qs_env ...
679 : !> \param tb ...
680 : !> \param energy ...
681 : ! **************************************************************************************************
682 25680 : SUBROUTINE tb_get_energy(qs_env, tb, energy)
683 :
684 : TYPE(qs_environment_type), POINTER :: qs_env
685 : TYPE(tblite_type), POINTER :: tb
686 : TYPE(qs_energy_type), POINTER :: energy
687 :
688 : #if defined(__TBLITE)
689 :
690 : INTEGER :: iounit
691 : TYPE(cp_logger_type), POINTER :: logger
692 : TYPE(section_vals_type), POINTER :: scf_section
693 : REAL(KIND=dp) :: xtb_inter
694 25680 : NULLIFY (scf_section, logger)
695 :
696 25680 : logger => cp_get_default_logger()
697 25680 : iounit = cp_logger_get_default_io_unit(logger)
698 25680 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
699 :
700 148766 : energy%repulsive = SUM(tb%e_rep)
701 148766 : energy%el_stat = SUM(tb%e_es)
702 148766 : energy%dispersion = SUM(tb%e_disp)
703 148766 : energy%dispersion_sc = SUM(tb%e_scd)
704 148766 : energy%xtb_xb_inter = SUM(tb%e_hal)
705 148766 : xtb_inter = SUM(tb%e_int)
706 :
707 : energy%total = energy%core + energy%repulsive + energy%el_stat + energy%dispersion &
708 : + energy%dispersion_sc + energy%xtb_xb_inter + xtb_inter &
709 25680 : + energy%kTS + energy%efield + energy%qmmm_el
710 :
711 : iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
712 25680 : extension=".scfLog")
713 25680 : IF (iounit > 0) THEN
714 : WRITE (UNIT=iounit, FMT="(/,(T9,A,T60,F20.10))") &
715 161 : "Repulsive pair potential energy: ", energy%repulsive, &
716 161 : "Zeroth order Hamiltonian energy: ", energy%core, &
717 161 : "Electrostatic energy: ", energy%el_stat, &
718 161 : "Self-consistent dispersion energy: ", energy%dispersion_sc, &
719 322 : "Non-self consistent dispersion energy: ", energy%dispersion
720 161 : IF (ABS(energy%xtb_xb_inter) > 1.e-9_dp) THEN
721 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
722 0 : "Correction for halogen bonding: ", energy%xtb_xb_inter
723 : END IF
724 161 : IF (ABS(xtb_inter) > 1.e-9_dp) THEN
725 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
726 0 : "Additional interaction (e.g. spin): ", xtb_inter
727 : END IF
728 161 : IF (ABS(energy%efield) > 1.e-9_dp) THEN
729 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
730 0 : "Electric field interaction energy: ", energy%efield
731 : END IF
732 161 : IF (qs_env%qmmm) THEN
733 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
734 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
735 : END IF
736 : END IF
737 : CALL cp_print_key_finished_output(iounit, logger, scf_section, &
738 25680 : "PRINT%DETAILED_ENERGY")
739 :
740 : #else
741 : MARK_USED(qs_env)
742 : MARK_USED(tb)
743 : MARK_USED(energy)
744 : CPABORT("Built without TBLITE")
745 : #endif
746 :
747 25680 : END SUBROUTINE tb_get_energy
748 :
749 : ! **************************************************************************************************
750 : !> \brief ...
751 : !> \param tb ...
752 : !> \param gto_basis_set ...
753 : !> \param element_symbol ...
754 : !> \param param ...
755 : !> \param occ ...
756 : ! **************************************************************************************************
757 334 : SUBROUTINE tb_get_basis(tb, gto_basis_set, element_symbol, param, occ)
758 :
759 : TYPE(tblite_type), POINTER :: tb
760 : TYPE(gto_basis_set_type), POINTER :: gto_basis_set
761 : CHARACTER(len=2), INTENT(IN) :: element_symbol
762 : TYPE(xtb_atom_type), POINTER :: param
763 : INTEGER, DIMENSION(5), INTENT(out) :: occ
764 :
765 : #if defined(__TBLITE)
766 :
767 : REAL(KIND=dp) :: docc
768 : CHARACTER(LEN=default_string_length) :: sng
769 : INTEGER :: ang, i_type, id_atom, ind_ao, ipgf, iSH, &
770 : ishell, ityp, maxl, mprim, natorb, &
771 : nset, nshell
772 : LOGICAL :: do_ortho
773 :
774 334 : CALL allocate_gto_basis_set(gto_basis_set)
775 :
776 : !identifying element in the bas data
777 334 : CALL symbol_to_number(i_type, element_symbol)
778 684 : DO id_atom = 1, tb%mol%nat
779 684 : IF (i_type == tb%el_num(id_atom)) EXIT
780 : END DO
781 334 : param%z = i_type
782 334 : param%symbol = element_symbol
783 334 : param%defined = .TRUE.
784 334 : ityp = tb%mol%id(id_atom)
785 :
786 : !getting size information
787 334 : nset = tb%calc%bas%nsh_id(ityp)
788 334 : nshell = 1
789 334 : mprim = 0
790 1096 : DO ishell = 1, nset
791 1096 : mprim = MAX(mprim, tb%calc%bas%cgto(ishell, ityp)%nprim)
792 : END DO
793 334 : param%nshell = nset
794 334 : natorb = 0
795 :
796 : !write basis set information
797 334 : CALL integer_to_string(mprim, sng)
798 334 : gto_basis_set%name = element_symbol//"_STO-"//TRIM(sng)//"G"
799 334 : gto_basis_set%nset = nset
800 334 : CALL reallocate(gto_basis_set%lmax, 1, nset)
801 334 : CALL reallocate(gto_basis_set%lmin, 1, nset)
802 334 : CALL reallocate(gto_basis_set%npgf, 1, nset)
803 334 : CALL reallocate(gto_basis_set%nshell, 1, nset)
804 334 : CALL reallocate(gto_basis_set%n, 1, 1, 1, nset)
805 334 : CALL reallocate(gto_basis_set%l, 1, 1, 1, nset)
806 334 : CALL reallocate(gto_basis_set%zet, 1, mprim, 1, nset)
807 334 : CALL reallocate(gto_basis_set%gcc, 1, mprim, 1, 1, 1, nset)
808 :
809 334 : ind_ao = 0
810 334 : maxl = 0
811 1096 : DO ishell = 1, nset
812 762 : ang = tb%calc%bas%cgto(ishell, ityp)%ang
813 762 : natorb = natorb + (2*ang + 1)
814 762 : param%lval(ishell) = ang
815 762 : maxl = MAX(ang, maxl)
816 762 : gto_basis_set%lmax(ishell) = ang
817 762 : gto_basis_set%lmin(ishell) = ang
818 762 : gto_basis_set%npgf(ishell) = tb%calc%bas%cgto(ishell, ityp)%nprim
819 762 : gto_basis_set%nshell(ishell) = nshell
820 762 : gto_basis_set%n(1, ishell) = ang + 1
821 762 : gto_basis_set%l(1, ishell) = ang
822 4530 : DO ipgf = 1, gto_basis_set%npgf(ishell)
823 3768 : gto_basis_set%gcc(ipgf, 1, ishell) = tb%calc%bas%cgto(ishell, ityp)%coeff(ipgf)
824 4530 : gto_basis_set%zet(ipgf, ishell) = tb%calc%bas%cgto(ishell, ityp)%alpha(ipgf)
825 : END DO
826 2810 : DO ipgf = 1, (2*ang + 1)
827 1714 : ind_ao = ind_ao + 1
828 1714 : param%lao(ind_ao) = ang
829 2476 : param%nao(ind_ao) = ishell
830 : END DO
831 : END DO
832 :
833 334 : do_ortho = .FALSE.
834 334 : CALL process_gto_basis(gto_basis_set, do_ortho, nset, maxl)
835 :
836 : !setting additional values in parameter
837 334 : param%rcut = get_cutoff(tb%calc%bas, tb%accuracy)
838 334 : param%natorb = natorb
839 334 : param%lmax = maxl !max angular momentum
840 :
841 : !getting occupation
842 334 : occ = 0
843 334 : docc = 0.0_dp
844 334 : IF (tb%calc%bas%nsh_at(id_atom) > 5) CPABORT("too many shells in tblite")
845 1096 : DO iSh = 1, tb%calc%bas%nsh_at(id_atom)
846 762 : occ(iSh) = NINT(tb%calc%h0%refocc(iSh, ityp) + docc)
847 762 : docc = docc + tb%calc%h0%refocc(iSh, ityp) - REAL(occ(iSh))
848 1096 : param%occupation(iSh) = occ(iSh)
849 : END DO
850 334 : IF (ABS(docc) > 0.1_dp) CPABORT("Getting occupation numbers from tblite fails")
851 2004 : param%zeff = SUM(occ) !effective core charge
852 :
853 : !set normalization process
854 334 : gto_basis_set%norm_type = 3
855 :
856 : #else
857 : occ = 0
858 : MARK_USED(tb)
859 : MARK_USED(gto_basis_set)
860 : MARK_USED(element_symbol)
861 : MARK_USED(param)
862 : CPABORT("Built without TBLITE")
863 : #endif
864 :
865 334 : END SUBROUTINE tb_get_basis
866 :
867 : ! **************************************************************************************************
868 : !> \brief ...
869 : !> \param qs_env ...
870 : !> \param calculate_forces ...
871 : ! **************************************************************************************************
872 2770 : SUBROUTINE build_tblite_matrices(qs_env, calculate_forces)
873 :
874 : TYPE(qs_environment_type), POINTER :: qs_env
875 : LOGICAL, INTENT(IN) :: calculate_forces
876 :
877 : #if defined(__TBLITE)
878 :
879 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_tblite_matrices'
880 :
881 : INTEGER :: handle, maxder, maxsgf, nderivatives, nimg, img, nkind, i, &
882 : ic, iw, iatom, jatom, ikind, jkind, iset, jset, n1, n2, icol, &
883 : irow, ia, ib, sgfa, sgfb, ldsab, nseta, nsetb, &
884 : natorb_a, natorb_b, raw_iatom, raw_jatom, &
885 : slot, idim, jdim
886 : LOGICAL :: found, norml1, norml2, use_arnoldi
887 : REAL(KIND=dp) :: dr, dshpoly, ff, hij_base, r2, rr
888 : REAL(KIND=dp) :: native_dot_tmp, native_cn_icol, native_cn_irow, &
889 : native_dot_weight
890 : INTEGER, DIMENSION(3) :: cell
891 : REAL(KIND=dp) :: hij, shpoly
892 : REAL(KIND=dp), DIMENSION(2) :: condnum
893 : REAL(KIND=dp), DIMENSION(3) :: native_h0_overlap_force, native_radial_force, raw_rij, rij
894 : REAL(KIND=dp), DIMENSION(3, 3) :: native_radial_dot
895 2770 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
896 2770 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
897 2770 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
898 2770 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: native_cn_deriv_sum, local_cn_deriv
899 2770 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: native_radial_force_sum, local_radial_force
900 5540 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint, hint
901 2770 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: radial_hint
902 2770 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min
903 2770 : INTEGER, DIMENSION(:), POINTER :: npgfa, npgfb, nsgfa, nsgfb
904 2770 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
905 2770 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
906 2770 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, zeta, zetb, scon_a, scon_b
907 5540 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, sblock
908 36010 : TYPE(block_p_type), DIMENSION(3, 3) :: radial_blocks
909 : !$ INTEGER :: hash, hash1, lock_num
910 : !$ INTEGER(KIND=int_8) :: iatom8
911 2770 : !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
912 : INTEGER, PARAMETER :: nlock = 501
913 :
914 2770 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
915 : TYPE(atprop_type), POINTER :: atprop
916 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
917 : TYPE(cp_logger_type), POINTER :: logger
918 2770 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_dh_native, matrix_h, matrix_p, &
919 2770 : matrix_q_native, matrix_s, matrix_s_native, matrix_w
920 : TYPE(dft_control_type), POINTER :: dft_control
921 2770 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
922 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
923 2770 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
924 : TYPE(kpoint_type), POINTER :: kpoints
925 2770 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
926 2770 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_kp
927 : TYPE(mp_para_env_type), POINTER :: para_env
928 : TYPE(qs_energy_type), POINTER :: energy
929 : TYPE(qs_ks_env_type), POINTER :: ks_env
930 2770 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
931 : TYPE(qs_rho_type), POINTER :: rho
932 : TYPE(tblite_type), POINTER :: tb
933 : TYPE(tb_hamiltonian), POINTER :: h0
934 : TYPE(virial_type), POINTER :: virial
935 :
936 2770 : CALL timeset(routineN, handle)
937 :
938 2770 : NULLIFY (ks_env, energy, atomic_kind_set, qs_kind_set)
939 2770 : NULLIFY (matrix_dh_native, matrix_h, matrix_q_native, matrix_s, matrix_s_native, atprop, dft_control)
940 2770 : NULLIFY (sab_orb, sab_kp, rho, tb, kpoints, cell_to_index)
941 :
942 : CALL get_qs_env(qs_env=qs_env, &
943 : ks_env=ks_env, para_env=para_env, &
944 : energy=energy, &
945 : atomic_kind_set=atomic_kind_set, &
946 : qs_kind_set=qs_kind_set, &
947 : matrix_h_kp=matrix_h, &
948 : matrix_s_kp=matrix_s, &
949 : atprop=atprop, &
950 : dft_control=dft_control, &
951 : sab_orb=sab_orb, &
952 : sab_kp=sab_kp, &
953 2770 : rho=rho, tb_tblite=tb)
954 2770 : h0 => tb%calc%h0
955 :
956 : !update geometry (required for debug / geometry optimization)
957 2770 : CALL tb_update_geometry(qs_env, tb)
958 :
959 2770 : nkind = SIZE(atomic_kind_set)
960 2770 : nderivatives = 0
961 2770 : IF (calculate_forces) THEN
962 154 : nderivatives = 1
963 154 : IF (ALLOCATED(tb%grad)) DEALLOCATE (tb%grad)
964 462 : ALLOCATE (tb%grad(3, tb%mol%nat))
965 154 : IF (ALLOCATED(tb%dsedcn)) DEALLOCATE (tb%dsedcn)
966 462 : ALLOCATE (tb%dsedcn(tb%calc%bas%nsh))
967 154 : IF (ALLOCATED(tb%calc%ncoord)) THEN
968 154 : IF (ALLOCATED(tb%dcndr)) DEALLOCATE (tb%dcndr)
969 616 : ALLOCATE (tb%dcndr(3, tb%mol%nat, tb%mol%nat))
970 154 : IF (ALLOCATED(tb%dcndL)) DEALLOCATE (tb%dcndL)
971 462 : ALLOCATE (tb%dcndL(3, 3, tb%mol%nat))
972 : END IF
973 : ELSE
974 2616 : IF (ALLOCATED(tb%grad)) DEALLOCATE (tb%grad)
975 2616 : IF (ALLOCATED(tb%dcndr)) DEALLOCATE (tb%dcndr)
976 2616 : IF (ALLOCATED(tb%dcndL)) DEALLOCATE (tb%dcndL)
977 : END IF
978 2770 : maxder = ncoset(nderivatives)
979 2770 : nimg = dft_control%nimages
980 2770 : IF (nimg > 1) THEN
981 1296 : IF (.NOT. ASSOCIATED(sab_kp)) CPABORT("Missing k-point neighbor list for tblite")
982 1296 : sab_orb => sab_kp
983 1296 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
984 1296 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
985 : END IF
986 :
987 : !intialise hamiltonian
988 2770 : CALL tb_init_ham(qs_env, tb, para_env)
989 :
990 : ! get density matrtix
991 2770 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
992 :
993 : ! set up matrices for force calculations
994 2770 : IF (calculate_forces) THEN
995 154 : NULLIFY (force, matrix_w, virial)
996 : CALL get_qs_env(qs_env=qs_env, &
997 : matrix_w_kp=matrix_w, &
998 154 : virial=virial, force=force)
999 :
1000 154 : IF (SIZE(matrix_p, 1) == 2) THEN
1001 : ! SCF has restored the final density; preserve both spins before summation.
1002 22 : CALL tb_store_density_ref(tb, matrix_p)
1003 152 : DO img = 1, nimg
1004 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
1005 130 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1006 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
1007 152 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1008 : END DO
1009 : END IF
1010 222 : tb%use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1011 : END IF
1012 :
1013 2770 : IF (calculate_forces) THEN
1014 154 : CALL dbcsr_allocate_matrix_set(matrix_q_native, 1, nimg)
1015 5224 : DO img = 1, nimg
1016 5070 : ALLOCATE (matrix_q_native(1, img)%matrix)
1017 : CALL dbcsr_copy(matrix_q_native(1, img)%matrix, matrix_w(1, img)%matrix, &
1018 5224 : name="TBLITE NATIVE OVERLAP FORCE MATRIX")
1019 : END DO
1020 : END IF
1021 :
1022 2770 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
1023 2770 : IF (calculate_forces) THEN
1024 462 : ALLOCATE (native_radial_force_sum(3, SIZE(atom_of_kind)))
1025 462 : ALLOCATE (native_cn_deriv_sum(SIZE(atom_of_kind)))
1026 154 : native_cn_deriv_sum = 0.0_dp
1027 154 : native_radial_force_sum = 0.0_dp
1028 : END IF
1029 :
1030 : ! set up basis set lists
1031 12964 : ALLOCATE (basis_set_list(nkind))
1032 2770 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
1033 :
1034 : ! allocate overlap matrix
1035 2770 : CALL dbcsr_allocate_matrix_set(matrix_s, 1, nimg)
1036 : CALL create_sab_matrix(ks_env, matrix_s, "OVERLAP MATRIX", basis_set_list, basis_set_list, &
1037 2770 : sab_orb, .TRUE.)
1038 2770 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
1039 :
1040 : ! initialize H matrix
1041 2770 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
1042 83744 : DO img = 1, nimg
1043 80974 : ALLOCATE (matrix_h(1, img)%matrix)
1044 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, img)%matrix, &
1045 80974 : name="HAMILTONIAN MATRIX")
1046 83744 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
1047 : END DO
1048 2770 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
1049 2770 : IF (calculate_forces .AND. nimg > 1) THEN
1050 92 : CALL dbcsr_allocate_matrix_set(matrix_dh_native, 9, nimg)
1051 5100 : DO img = 1, nimg
1052 50172 : DO idim = 1, 9
1053 45072 : ALLOCATE (matrix_dh_native(idim, img)%matrix)
1054 : CALL dbcsr_create(matrix_dh_native(idim, img)%matrix, template=matrix_s(1, img)%matrix, &
1055 45072 : name="TBLITE H0 STRAIN DERIVATIVE")
1056 50080 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_dh_native(idim, img)%matrix, sab_orb)
1057 : END DO
1058 : END DO
1059 : END IF
1060 2770 : CALL get_qs_kind_set(qs_kind_set, maxsgf=maxsgf)
1061 2770 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
1062 :
1063 2770 : native_radial_dot = 0.0_dp
1064 :
1065 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
1066 : !$OMP PARALLEL DEFAULT(NONE) &
1067 : !$OMP SHARED (calculate_forces, basis_set_list, sab_orb, matrix_dh_native, matrix_s, matrix_h, matrix_p, &
1068 : !$OMP cell_to_index, nimg, atom_of_kind, qs_kind_set, tb, h0, ldsab, maxder, maxsgf, locks, &
1069 : !$OMP native_cn_deriv_sum, native_radial_force_sum, ncoset) &
1070 : !$OMP PRIVATE (slot, hash, hash1, iatom8, lock_num, ikind, jkind, iatom, jatom, cell, rij, ic, irow, &
1071 : !$OMP icol, dr, dshpoly, ff, hij_base, r2, raw_iatom, raw_jatom, &
1072 : !$OMP raw_rij, native_cn_icol, native_cn_irow, native_dot_weight, native_h0_overlap_force, &
1073 : !$OMP native_radial_force, local_cn_deriv, local_radial_force, &
1074 : !$OMP n1, n2, ia, ib, i, idim, jdim, iset, jset, sgfa, sgfb, &
1075 : !$OMP nseta, nsetb, natorb_a, &
1076 : !$OMP natorb_b, found, zeta, first_sgfa, la_max, la_min, npgfa, nsgfa, rpgfa, set_radius_a, &
1077 : !$OMP scon_a, first_sgfb, lb_max, lb_min, npgfb, nsgfb, rpgfb, set_radius_b, scon_b, zetb, rr, hij, shpoly, &
1078 : !$OMP owork, oint, sint, hint, radial_hint, basis_set_a, basis_set_b, sblock, fblock, &
1079 : !$OMP pblock, radial_blocks) &
1080 2770 : !$OMP REDUCTION (+ : native_radial_dot)
1081 :
1082 : !$OMP SINGLE
1083 : !$ ALLOCATE (locks(nlock))
1084 : !$OMP END SINGLE
1085 : !$OMP DO
1086 : !$ DO lock_num = 1, nlock
1087 : !$ CALL omp_init_lock(locks(lock_num))
1088 : !$ END DO
1089 : !$OMP END DO
1090 :
1091 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
1092 : ALLOCATE (sint(maxsgf, maxsgf, maxder), hint(maxsgf, maxsgf, 1))
1093 : IF (calculate_forces) THEN
1094 : ALLOCATE (local_cn_deriv(SIZE(atom_of_kind)), local_radial_force(3, SIZE(atom_of_kind)))
1095 : local_cn_deriv = 0.0_dp
1096 : local_radial_force = 0.0_dp
1097 : IF (nimg > 1) ALLOCATE (radial_hint(maxsgf, maxsgf, 3, 3))
1098 : END IF
1099 :
1100 : !$OMP DO SCHEDULE(GUIDED)
1101 : DO slot = 1, sab_orb(1)%nl_size
1102 :
1103 : ikind = sab_orb(1)%nlist_task(slot)%ikind
1104 : jkind = sab_orb(1)%nlist_task(slot)%jkind
1105 : iatom = sab_orb(1)%nlist_task(slot)%iatom
1106 : jatom = sab_orb(1)%nlist_task(slot)%jatom
1107 : cell(:) = sab_orb(1)%nlist_task(slot)%cell(:)
1108 : rij(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
1109 : raw_iatom = iatom
1110 : raw_jatom = jatom
1111 : raw_rij = rij
1112 : native_cn_icol = 0.0_dp
1113 : native_cn_irow = 0.0_dp
1114 : native_h0_overlap_force = 0.0_dp
1115 : native_radial_force = 0.0_dp
1116 :
1117 : ! canonicalize pair ordering as in current serial flow
1118 : icol = MAX(iatom, jatom)
1119 : irow = MIN(iatom, jatom)
1120 : IF (iatom < jatom) THEN
1121 : rij = -rij
1122 : i = ikind
1123 : ikind = jkind
1124 : jkind = i
1125 : END IF
1126 :
1127 : dr = NORM2(rij(:))
1128 : r2 = dr*dr
1129 :
1130 : IF (nimg == 1) THEN
1131 : ic = 1
1132 : ELSE
1133 : ic = cell_to_index(cell(1), cell(2), cell(3))
1134 : CPASSERT(ic > 0)
1135 : END IF
1136 :
1137 : NULLIFY (sblock)
1138 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
1139 : row=irow, col=icol, BLOCK=sblock, found=found)
1140 : CPASSERT(found)
1141 : NULLIFY (fblock)
1142 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
1143 : row=irow, col=icol, BLOCK=fblock, found=found)
1144 : CPASSERT(found)
1145 : IF (calculate_forces) THEN
1146 : NULLIFY (pblock)
1147 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
1148 : row=irow, col=icol, BLOCK=pblock, found=found)
1149 : CPASSERT(found)
1150 : IF (nimg > 1) THEN
1151 : DO jdim = 1, 3
1152 : DO idim = 1, 3
1153 : NULLIFY (radial_blocks(idim, jdim)%block)
1154 : CALL dbcsr_get_block_p(matrix=matrix_dh_native(idim + 3*(jdim - 1), ic)%matrix, &
1155 : row=irow, col=icol, BLOCK=radial_blocks(idim, jdim)%block, &
1156 : found=found)
1157 : CPASSERT(found)
1158 : END DO
1159 : END DO
1160 : END IF
1161 : END IF
1162 :
1163 : ! --------- Overlap
1164 : !get basis information
1165 : basis_set_a => basis_set_list(ikind)%gto_basis_set
1166 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
1167 : basis_set_b => basis_set_list(jkind)%gto_basis_set
1168 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
1169 : ! basis a
1170 : first_sgfa => basis_set_a%first_sgf
1171 : la_max => basis_set_a%lmax
1172 : la_min => basis_set_a%lmin
1173 : npgfa => basis_set_a%npgf
1174 : nseta = basis_set_a%nset
1175 : nsgfa => basis_set_a%nsgf_set
1176 : rpgfa => basis_set_a%pgf_radius
1177 : set_radius_a => basis_set_a%set_radius
1178 : scon_a => basis_set_a%scon
1179 : zeta => basis_set_a%zet
1180 : ! basis b
1181 : first_sgfb => basis_set_b%first_sgf
1182 : lb_max => basis_set_b%lmax
1183 : lb_min => basis_set_b%lmin
1184 : npgfb => basis_set_b%npgf
1185 : nsetb = basis_set_b%nset
1186 : nsgfb => basis_set_b%nsgf_set
1187 : rpgfb => basis_set_b%pgf_radius
1188 : set_radius_b => basis_set_b%set_radius
1189 : scon_b => basis_set_b%scon
1190 : zetb => basis_set_b%zet
1191 :
1192 : natorb_a = 0
1193 : DO iset = 1, nseta
1194 : natorb_a = natorb_a + (2*basis_set_a%l(1, iset) + 1)
1195 : END DO
1196 : natorb_b = 0
1197 : DO iset = 1, nsetb
1198 : natorb_b = natorb_b + (2*basis_set_b%l(1, iset) + 1)
1199 : END DO
1200 : sint(1:natorb_a, 1:natorb_b, :) = 0.0_dp
1201 : hint(1:natorb_a, 1:natorb_b, :) = 0.0_dp
1202 : IF (calculate_forces .AND. nimg > 1) THEN
1203 : radial_hint(1:natorb_a, 1:natorb_b, :, :) = 0.0_dp
1204 : END IF
1205 :
1206 : !----------------- overlap integrals
1207 : DO iset = 1, nseta
1208 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1209 : sgfa = first_sgfa(1, iset)
1210 : DO jset = 1, nsetb
1211 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
1212 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1213 : sgfb = first_sgfb(1, jset)
1214 : IF (calculate_forces) THEN
1215 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1216 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1217 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
1218 : ELSE
1219 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1220 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1221 : rij, sab=oint(:, :, 1))
1222 : END IF
1223 : ! Contraction
1224 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1225 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1226 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
1227 : IF (calculate_forces) THEN
1228 : DO i = 2, 4
1229 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1230 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1231 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
1232 : END DO
1233 : END IF
1234 : END DO
1235 : END DO
1236 :
1237 : !$ iatom8 = INT(irow - 1, int_8)*INT(SIZE(atom_of_kind), int_8) + INT(icol, int_8)
1238 : !$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
1239 : ! update S matrix
1240 : !$ hash = hash1
1241 : !$ CALL omp_set_lock(locks(hash))
1242 : IF (icol <= irow) THEN
1243 : sblock(:, :) = sblock(:, :) + sint(1:natorb_a, 1:natorb_b, 1)
1244 : ELSE
1245 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(1:natorb_a, 1:natorb_b, 1))
1246 : END IF
1247 : !$ CALL omp_unset_lock(locks(hash))
1248 :
1249 : ! --------- Hamiltonian
1250 : IF (icol == irow .AND. dr < same_atom) THEN
1251 : !get diagonal F matrix from selfenergy
1252 : n1 = tb%calc%bas%ish_at(icol)
1253 : DO iset = 1, nseta
1254 : sgfa = first_sgfa(1, iset)
1255 : hij = tb%selfenergy(n1 + iset)
1256 : DO ia = sgfa, sgfa + nsgfa(iset) - 1
1257 : hint(ia, ia, 1) = hij
1258 : IF (calculate_forces) THEN
1259 : native_cn_icol = native_cn_icol + tb%dsedcn(n1 + iset)*pblock(ia, ia)
1260 : DO jdim = 1, 3
1261 : DO idim = 1, 3
1262 : IF (nimg == 1) THEN
1263 : native_radial_dot(idim, jdim) = native_radial_dot(idim, jdim) + &
1264 : tb%dsedcn(n1 + iset)*tb%dcndL(idim, jdim, icol)*pblock(ia, ia)
1265 : ELSE
1266 : radial_hint(ia, ia, idim, jdim) = radial_hint(ia, ia, idim, jdim) + &
1267 : tb%dsedcn(n1 + iset)*tb%dcndL(idim, jdim, icol)
1268 : END IF
1269 : END DO
1270 : END DO
1271 : END IF
1272 : END DO
1273 : END DO
1274 : ELSE
1275 : !get off-diagonal F matrix
1276 : rr = SQRT(dr/(h0%rad(jkind) + h0%rad(ikind)))
1277 : n1 = tb%calc%bas%ish_at(icol)
1278 : DO iset = 1, nseta
1279 : sgfa = first_sgfa(1, iset)
1280 : n2 = tb%calc%bas%ish_at(irow)
1281 : DO jset = 1, nsetb
1282 : sgfb = first_sgfb(1, jset)
1283 : shpoly = (1.0_dp + h0%shpoly(iset, ikind)*rr) &
1284 : *(1.0_dp + h0%shpoly(jset, jkind)*rr)
1285 : dshpoly = ((1.0_dp + h0%shpoly(iset, ikind)*rr)*h0%shpoly(jset, jkind)*rr &
1286 : + (1.0_dp + h0%shpoly(jset, jkind)*rr)*h0%shpoly(iset, ikind)*rr) &
1287 : *0.5_dp/r2
1288 : hij_base = 0.5_dp*(tb%selfenergy(n1 + iset) + tb%selfenergy(n2 + jset)) &
1289 : *h0%hscale(iset, jset, ikind, jkind)
1290 : hij = hij_base*shpoly
1291 : DO ia = sgfa, sgfa + nsgfa(iset) - 1
1292 : DO ib = sgfb, sgfb + nsgfb(jset) - 1
1293 : hint(ia, ib, 1) = hij*sint(ia, ib, 1)
1294 : IF (calculate_forces) THEN
1295 : native_dot_weight = 2.0_dp
1296 : IF (icol == irow) native_dot_weight = 1.0_dp
1297 : native_cn_icol = native_cn_icol + native_dot_weight*0.5_dp* &
1298 : h0%hscale(iset, jset, ikind, jkind)*shpoly* &
1299 : tb%dsedcn(n1 + iset)*pblock(ib, ia)*sint(ia, ib, 1)
1300 : native_cn_irow = native_cn_irow + native_dot_weight*0.5_dp* &
1301 : h0%hscale(iset, jset, ikind, jkind)*shpoly* &
1302 : tb%dsedcn(n2 + jset)*pblock(ib, ia)*sint(ia, ib, 1)
1303 : native_radial_force = native_radial_force + &
1304 : hij_base*dshpoly*pblock(ib, ia)*sint(ia, ib, 1)*raw_rij
1305 : native_h0_overlap_force = native_h0_overlap_force + &
1306 : hij*pblock(ib, ia)*sint(ia, ib, 2:4)
1307 : DO jdim = 1, 3
1308 : DO idim = 1, 3
1309 : IF (nimg == 1) THEN
1310 : native_radial_dot(idim, jdim) = native_radial_dot(idim, jdim) + &
1311 : native_dot_weight*pblock(ib, ia)*( &
1312 : (-hij)*sint(ia, ib, idim + 1)*rij(jdim) + &
1313 : hij_base*dshpoly*sint(ia, ib, 1)*rij(idim)*rij(jdim) + &
1314 : 0.5_dp*h0%hscale(iset, jset, ikind, jkind)*shpoly*sint(ia, ib, 1)* &
1315 : (tb%dsedcn(n1 + iset)*tb%dcndL(idim, jdim, icol) + &
1316 : tb%dsedcn(n2 + jset)*tb%dcndL(idim, jdim, irow)))
1317 : ELSE
1318 : radial_hint(ia, ib, idim, jdim) = radial_hint(ia, ib, idim, jdim) + &
1319 : (-hij)*sint(ia, ib, idim + 1)*rij(jdim) + &
1320 : hij_base*dshpoly*sint(ia, ib, 1)*rij(idim)*rij(jdim) + &
1321 : 0.5_dp*h0%hscale(iset, jset, ikind, jkind)*shpoly*sint(ia, ib, 1)* &
1322 : (tb%dsedcn(n1 + iset)*tb%dcndL(idim, jdim, icol) + &
1323 : tb%dsedcn(n2 + jset)*tb%dcndL(idim, jdim, irow))
1324 : END IF
1325 : END DO
1326 : END DO
1327 : END IF
1328 : END DO
1329 : END DO
1330 : END DO
1331 : END DO
1332 : END IF
1333 :
1334 : ! update F matrix
1335 : !$ CALL omp_set_lock(locks(hash))
1336 : IF (icol <= irow) THEN
1337 : fblock(:, :) = fblock(:, :) + hint(1:natorb_a, 1:natorb_b, 1)
1338 : IF (calculate_forces .AND. nimg > 1) THEN
1339 : DO jdim = 1, 3
1340 : DO idim = 1, 3
1341 : radial_blocks(idim, jdim)%block(:, :) = radial_blocks(idim, jdim)%block(:, :) + &
1342 : radial_hint(1:natorb_a, 1:natorb_b, idim, jdim)
1343 : END DO
1344 : END DO
1345 : END IF
1346 : ELSE
1347 : fblock(:, :) = fblock(:, :) + TRANSPOSE(hint(1:natorb_a, 1:natorb_b, 1))
1348 : IF (calculate_forces .AND. nimg > 1) THEN
1349 : DO jdim = 1, 3
1350 : DO idim = 1, 3
1351 : radial_blocks(idim, jdim)%block(:, :) = radial_blocks(idim, jdim)%block(:, :) + &
1352 : TRANSPOSE(radial_hint(1:natorb_a, 1:natorb_b, idim, jdim))
1353 : END DO
1354 : END DO
1355 : END IF
1356 : END IF
1357 : !$ CALL omp_unset_lock(locks(hash))
1358 :
1359 : IF (calculate_forces) THEN
1360 : ff = 2.0_dp
1361 : local_radial_force(:, raw_iatom) = &
1362 : local_radial_force(:, raw_iatom) - ff*native_radial_force
1363 : local_radial_force(:, raw_jatom) = &
1364 : local_radial_force(:, raw_jatom) + ff*native_radial_force
1365 : local_radial_force(:, icol) = &
1366 : local_radial_force(:, icol) + ff*native_h0_overlap_force
1367 : local_radial_force(:, irow) = &
1368 : local_radial_force(:, irow) - ff*native_h0_overlap_force
1369 : local_cn_deriv(icol) = local_cn_deriv(icol) + native_cn_icol
1370 : local_cn_deriv(irow) = local_cn_deriv(irow) + native_cn_irow
1371 : END IF
1372 :
1373 : END DO
1374 : !$OMP END DO
1375 : IF (calculate_forces) THEN
1376 : !$OMP CRITICAL(tblite_native_radial_force)
1377 : native_radial_force_sum(:, :) = native_radial_force_sum(:, :) + local_radial_force(:, :)
1378 : native_cn_deriv_sum(:) = native_cn_deriv_sum(:) + local_cn_deriv(:)
1379 : !$OMP END CRITICAL(tblite_native_radial_force)
1380 : DEALLOCATE (local_radial_force, local_cn_deriv)
1381 : END IF
1382 : DEALLOCATE (oint, owork, sint, hint)
1383 : IF (ALLOCATED(radial_hint)) DEALLOCATE (radial_hint)
1384 : !$OMP DO
1385 : !$ DO lock_num = 1, nlock
1386 : !$ CALL omp_destroy_lock(locks(lock_num))
1387 : !$ END DO
1388 : !$OMP END DO
1389 : !$OMP SINGLE
1390 : !$ DEALLOCATE (locks)
1391 : !$OMP END SINGLE NOWAIT
1392 : !$OMP END PARALLEL
1393 :
1394 2770 : IF (calculate_forces) THEN
1395 154 : IF (nimg > 1) THEN
1396 92 : native_radial_dot = 0.0_dp
1397 5100 : DO img = 1, nimg
1398 20124 : DO jdim = 1, 3
1399 65104 : DO idim = 1, 3
1400 45072 : CALL dbcsr_finalize(matrix_dh_native(idim + 3*(jdim - 1), img)%matrix)
1401 : CALL dbcsr_dot(matrix_dh_native(idim + 3*(jdim - 1), img)%matrix, &
1402 45072 : matrix_p(1, img)%matrix, native_dot_tmp)
1403 60096 : native_radial_dot(idim, jdim) = native_radial_dot(idim, jdim) + native_dot_tmp
1404 : END DO
1405 : END DO
1406 : END DO
1407 92 : CALL dbcsr_deallocate_matrix_set(matrix_dh_native)
1408 : ELSE
1409 62 : CALL para_env%sum(native_radial_dot)
1410 : END IF
1411 154 : CALL para_env%sum(native_cn_deriv_sum)
1412 3754 : tb%grad = 0.0_dp
1413 154 : CALL tb_add_grad(tb%grad, tb%dcndr, native_cn_deriv_sum, tb%mol%nat)
1414 154 : CALL tb_grad2force(qs_env, tb, para_env, 4)
1415 1054 : DO iatom = 1, SIZE(atom_of_kind)
1416 900 : ikind = kind_of(iatom)
1417 : force(ikind)%overlap(:, atom_of_kind(iatom)) = &
1418 3754 : force(ikind)%overlap(:, atom_of_kind(iatom)) + native_radial_force_sum(:, iatom)
1419 : END DO
1420 154 : IF (tb%use_virial) THEN
1421 1118 : virial%pv_overlap = virial%pv_overlap - native_radial_dot/para_env%num_pe
1422 1118 : virial%pv_virial = virial%pv_virial - native_radial_dot/para_env%num_pe
1423 : END IF
1424 : END IF
1425 :
1426 83744 : DO img = 1, nimg
1427 161948 : DO i = 1, SIZE(matrix_s, 1)
1428 161948 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1429 : END DO
1430 164718 : DO i = 1, SIZE(matrix_h, 1)
1431 161948 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1432 : END DO
1433 : END DO
1434 :
1435 2770 : IF (calculate_forces) THEN
1436 : CALL build_overlap_matrix(ks_env, matrixkp_s=matrix_s_native, &
1437 : matrix_name="OVERLAP MATRIX", nderivative=nderivatives, &
1438 : basis_type_a="ORB", basis_type_b="ORB", sab_nl=sab_orb, &
1439 154 : calculate_forces=.TRUE., matrixkp_p=matrix_q_native)
1440 154 : CALL dbcsr_deallocate_matrix_set(matrix_s_native)
1441 154 : CALL dbcsr_deallocate_matrix_set(matrix_q_native)
1442 : END IF
1443 :
1444 : !compute multipole moments for gfn2
1445 2770 : IF (dft_control%qs_control%xtb_control%tblite_method == gfn2xtb) THEN
1446 1382 : CALL tb_get_multipole(qs_env, tb)
1447 : END IF
1448 :
1449 : ! output overlap information
1450 2770 : NULLIFY (logger)
1451 2770 : logger => cp_get_default_logger()
1452 2770 : IF (.NOT. calculate_forces) THEN
1453 2616 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
1454 : "DFT%PRINT%OVERLAP_CONDITION") /= 0) THEN
1455 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
1456 2 : extension=".Log")
1457 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
1458 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
1459 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
1460 2 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
1461 2 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
1462 : END IF
1463 : END IF
1464 :
1465 2770 : DEALLOCATE (basis_set_list)
1466 2770 : IF (ALLOCATED(native_cn_deriv_sum)) DEALLOCATE (native_cn_deriv_sum)
1467 2770 : IF (ALLOCATED(native_radial_force_sum)) DEALLOCATE (native_radial_force_sum)
1468 :
1469 2770 : CALL timestop(handle)
1470 :
1471 : #else
1472 : MARK_USED(qs_env)
1473 : MARK_USED(calculate_forces)
1474 : CPABORT("Built without TBLITE")
1475 : #endif
1476 :
1477 11080 : END SUBROUTINE build_tblite_matrices
1478 :
1479 : ! **************************************************************************************************
1480 : !> \brief ...
1481 : !> \param qs_env ...
1482 : !> \param dft_control ...
1483 : !> \param tb ...
1484 : !> \param calculate_forces ...
1485 : !> \param use_rho ...
1486 : ! **************************************************************************************************
1487 53808 : SUBROUTINE tb_update_charges(qs_env, dft_control, tb, calculate_forces, use_rho)
1488 :
1489 : TYPE(qs_environment_type), POINTER :: qs_env
1490 : TYPE(dft_control_type), POINTER :: dft_control
1491 : TYPE(tblite_type), POINTER :: tb
1492 : LOGICAL, INTENT(IN) :: calculate_forces
1493 : LOGICAL, INTENT(IN) :: use_rho
1494 :
1495 : #if defined(__TBLITE)
1496 :
1497 : INTEGER :: iatom, ikind, is, ns, atom_a, im
1498 : INTEGER :: ispin, nspin
1499 : INTEGER :: nkind, nsgf, natorb, na
1500 : INTEGER :: n_atom, max_orb, max_shell
1501 : LOGICAL :: advance_native_mixer, do_dipole, do_quadrupole, &
1502 : reuse_native_input, seed_native_from_rho, &
1503 : use_native_mixer, use_no_mixer
1504 : REAL(KIND=dp) :: native_seed_charge, norm, pao
1505 : INTEGER, DIMENSION(5) :: occ
1506 : INTEGER, DIMENSION(25) :: lao
1507 : INTEGER, DIMENSION(25) :: nao
1508 53808 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ch_atom
1509 53808 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg
1510 53808 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: aocg_spin, ao_dip_spin, ao_quad_spin, &
1511 53808 : ch_orb_spin, ch_shell_spin
1512 :
1513 53808 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1514 53808 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_p
1515 53808 : TYPE(error_type), ALLOCATABLE :: error
1516 : TYPE(mp_para_env_type), POINTER :: para_env
1517 53808 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1518 53808 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1519 : TYPE(qs_rho_type), POINTER :: rho
1520 : TYPE(qs_scf_env_type), POINTER :: scf_env
1521 : TYPE(scf_control_type), POINTER :: scf_control
1522 : TYPE(xtb_atom_type), POINTER :: xtb_kind
1523 :
1524 : ! compute mulliken charges required for charge update
1525 53808 : NULLIFY (particle_set, qs_kind_set, atomic_kind_set, scf_control, scf_env)
1526 : CALL get_qs_env(qs_env=qs_env, scf_env=scf_env, particle_set=particle_set, qs_kind_set=qs_kind_set, &
1527 : atomic_kind_set=atomic_kind_set, matrix_s_kp=matrix_s, rho=rho, para_env=para_env, &
1528 53808 : scf_control=scf_control)
1529 :
1530 : ! also compute multipoles needed by GFN2
1531 53808 : do_dipole = .FALSE.
1532 53808 : do_quadrupole = .FALSE.
1533 53808 : IF (dft_control%qs_control%do_ls_scf .OR. scf_control%use_ot) THEN
1534 : use_native_mixer = .FALSE.
1535 : use_no_mixer = .TRUE.
1536 : ELSE
1537 93674 : SELECT CASE (dft_control%qs_control%xtb_control%tblite_scc_mixer)
1538 : CASE (tblite_scc_mixer_auto)
1539 41618 : use_native_mixer = tb_native_scc_mixer_active(dft_control)
1540 41618 : use_no_mixer = .FALSE.
1541 : CASE (tblite_scc_mixer_tblite)
1542 : use_native_mixer = .TRUE.
1543 : use_no_mixer = .FALSE.
1544 : CASE (tblite_scc_mixer_cp2k)
1545 : use_native_mixer = .FALSE.
1546 : use_no_mixer = .FALSE.
1547 : CASE (tblite_scc_mixer_none)
1548 : use_native_mixer = .FALSE.
1549 0 : use_no_mixer = .TRUE.
1550 : CASE DEFAULT
1551 52056 : CPABORT("Unknown tblite SCC mixer")
1552 : END SELECT
1553 : END IF
1554 41618 : IF (use_native_mixer) THEN
1555 47596 : IF (.NOT. ASSOCIATED(scf_env)) CPABORT("tblite SCC mixer requires a QS SCF environment")
1556 47596 : IF (use_rho .AND. (.NOT. calculate_forces)) THEN
1557 24844 : IF (scf_env%iter_count > dft_control%qs_control%xtb_control%tblite_mixer_iterations) THEN
1558 0 : CPABORT("tblite SCC mixer exceeded TBLITE_MIXER/ITERATIONS")
1559 : END IF
1560 24844 : IF (scf_env%iter_count == 1) THEN
1561 : CALL tb_configure_mixer(tb, dft_control%qs_control%xtb_control%tblite_mixer_iterations, &
1562 : dft_control%qs_control%xtb_control%tblite_mixer_memory, &
1563 : dft_control%qs_control%xtb_control%tblite_mixer_damping, &
1564 : dft_control%qs_control%xtb_control%tblite_mixer_omega0, &
1565 : dft_control%qs_control%xtb_control%tblite_mixer_min_weight, &
1566 : dft_control%qs_control%xtb_control%tblite_mixer_max_weight, &
1567 : dft_control%qs_control%xtb_control%tblite_mixer_weight_factor, &
1568 2458 : dft_control%qs_control%xtb_control%tblite_mixer_solver)
1569 2458 : CALL tb_reset_mixer(tb)
1570 : END IF
1571 : END IF
1572 : END IF
1573 53808 : nspin = dft_control%nspins
1574 53808 : IF (nspin /= tb%wfn%nspin) CPABORT("CP2K/tblite spin channel mismatch")
1575 :
1576 53808 : NULLIFY (matrix_p)
1577 53808 : IF (use_rho) THEN
1578 29020 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1579 24788 : ELSE IF (calculate_forces .AND. nspin > 1) THEN
1580 22 : IF (.NOT. ASSOCIATED(tb%rho_ao_kp_ref)) THEN
1581 0 : CPABORT("Missing converged tblite density for UKS/LSD forces")
1582 : END IF
1583 22 : matrix_p => tb%rho_ao_kp_ref
1584 : ELSE
1585 24766 : matrix_p => scf_env%p_mix_new
1586 : END IF
1587 53808 : IF (ASSOCIATED(tb%dipbra)) do_dipole = .TRUE.
1588 53808 : IF (ASSOCIATED(tb%quadbra)) do_quadrupole = .TRUE.
1589 53808 : reuse_native_input = .FALSE.
1590 53808 : IF (use_native_mixer) THEN
1591 47596 : IF (scf_env%iter_count == 1) THEN
1592 8606 : reuse_native_input = ANY(ABS(tb%wfn%qsh) > 1.0E-14_dp)
1593 4918 : IF (do_dipole) reuse_native_input = reuse_native_input .OR. &
1594 5998 : ANY(ABS(tb%wfn%dpat) > 1.0E-14_dp)
1595 4918 : IF (do_quadrupole) reuse_native_input = reuse_native_input .OR. &
1596 8620 : ANY(ABS(tb%wfn%qpat) > 1.0E-14_dp)
1597 : END IF
1598 : END IF
1599 53808 : n_atom = SIZE(particle_set)
1600 53808 : nkind = SIZE(atomic_kind_set)
1601 53808 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
1602 215232 : ALLOCATE (aocg(nsgf, n_atom))
1603 269040 : ALLOCATE (aocg_spin(nsgf, n_atom, nspin))
1604 53808 : aocg = 0.0_dp
1605 53808 : IF (do_dipole) THEN
1606 149560 : ALLOCATE (ao_dip_spin(n_atom, dip_n, nspin))
1607 : END IF
1608 53808 : IF (do_quadrupole) THEN
1609 149560 : ALLOCATE (ao_quad_spin(n_atom, quad_n, nspin))
1610 : END IF
1611 : max_orb = 0
1612 135988 : max_shell = 0
1613 135988 : DO ikind = 1, nkind
1614 82180 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1615 82180 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb)
1616 135988 : max_orb = MAX(max_orb, natorb)
1617 : END DO
1618 310558 : DO is = 1, n_atom
1619 310558 : max_shell = MAX(max_shell, tb%calc%bas%nsh_at(is))
1620 : END DO
1621 215232 : ALLOCATE (ch_atom(n_atom, nspin))
1622 484272 : ALLOCATE (ch_orb_spin(max_orb, n_atom, nspin), ch_shell_spin(n_atom, max_shell, nspin))
1623 53808 : ch_atom = 0.0_dp
1624 53808 : ch_orb_spin = 0.0_dp
1625 53808 : ch_shell_spin = 0.0_dp
1626 53808 : CALL tb_contract_population(tb, matrix_p, matrix_s, para_env, aocg_spin, ao_dip_spin, ao_quad_spin)
1627 53808 : IF (nspin == 1) THEN
1628 2204332 : aocg(:, :) = aocg_spin(:, :, 1)
1629 : ELSE
1630 80214 : aocg(:, :) = aocg_spin(:, :, 1) + aocg_spin(:, :, 2)
1631 5502 : IF (do_dipole) THEN
1632 17040 : DO im = 1, dip_n
1633 49212 : DO iatom = 1, n_atom
1634 32172 : pao = ao_dip_spin(iatom, im, 1)
1635 32172 : ao_dip_spin(iatom, im, 1) = pao + ao_dip_spin(iatom, im, 2)
1636 44952 : ao_dip_spin(iatom, im, 2) = pao - ao_dip_spin(iatom, im, 2)
1637 : END DO
1638 : END DO
1639 : END IF
1640 5502 : IF (do_quadrupole) THEN
1641 29820 : DO im = 1, quad_n
1642 94164 : DO iatom = 1, n_atom
1643 64344 : pao = ao_quad_spin(iatom, im, 1)
1644 64344 : ao_quad_spin(iatom, im, 1) = pao + ao_quad_spin(iatom, im, 2)
1645 89904 : ao_quad_spin(iatom, im, 2) = pao - ao_quad_spin(iatom, im, 2)
1646 : END DO
1647 : END DO
1648 : END IF
1649 : END IF
1650 53808 : NULLIFY (xtb_kind)
1651 135988 : DO ikind = 1, nkind
1652 82180 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
1653 82180 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1654 82180 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, nao=nao, occupation=occ)
1655 474918 : DO iatom = 1, na
1656 256750 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1657 2096236 : DO is = 1, natorb
1658 : ! Preserve tblite's native shell order, including repeated angular momenta.
1659 1839486 : ns = nao(is)
1660 1839486 : norm = 2*lao(is) + 1
1661 1839486 : ch_orb_spin(is, atom_a, 1) = aocg(is, atom_a) - tb%calc%h0%refocc(nao(is), ikind)/norm
1662 1839486 : IF (nspin == 2) ch_orb_spin(is, atom_a, 2) = &
1663 49960 : aocg_spin(is, atom_a, 1) - aocg_spin(is, atom_a, 2)
1664 3985682 : DO ispin = 1, nspin
1665 : ch_shell_spin(atom_a, ns, ispin) = ch_orb_spin(is, atom_a, ispin) + &
1666 3728932 : ch_shell_spin(atom_a, ns, ispin)
1667 : END DO
1668 : END DO
1669 610532 : DO ispin = 1, nspin
1670 2562200 : ch_atom(atom_a, ispin) = SUM(ch_orb_spin(:, atom_a, ispin))
1671 : END DO
1672 : END DO
1673 : END DO
1674 310558 : native_seed_charge = -SUM(ch_atom(:, 1))
1675 : seed_native_from_rho = SUM(ABS(aocg_spin)) > 1.0E-10_dp .AND. &
1676 2418568 : ABS(native_seed_charge - REAL(dft_control%charge, dp)) < 1.0E-5_dp
1677 53808 : DEALLOCATE (aocg, aocg_spin)
1678 :
1679 : ! Population arrays use CP2K electron-excess signs. The spin axis already
1680 : ! contains total charge and magnetization, as required by tblite.
1681 53808 : IF (use_native_mixer) THEN
1682 47596 : IF (.NOT. ALLOCATED(tb%mixer)) CPABORT("tblite mixer not initialized")
1683 47596 : advance_native_mixer = .FALSE.
1684 47596 : IF (use_rho .AND. (.NOT. calculate_forces)) THEN
1685 24844 : advance_native_mixer = scf_env%iter_count > 1
1686 : END IF
1687 24844 : IF (advance_native_mixer) THEN
1688 : ! Build the next Hamiltonian from the previous SCC residual.
1689 22386 : CALL tb%mixer%next(error)
1690 22386 : IF (ALLOCATED(error)) CPABORT("tblite native mixer failed: "//error%message)
1691 22386 : CALL tb%mixer%get(tb%wfn%qsh)
1692 22386 : CALL tb_sum_shell_charges(tb)
1693 22386 : IF (do_dipole) CALL tb%mixer%get(tb%wfn%dpat)
1694 22386 : IF (do_quadrupole) CALL tb%mixer%get(tb%wfn%qpat)
1695 25210 : ELSE IF (use_rho .AND. (.NOT. calculate_forces)) THEN
1696 2458 : IF (.NOT. reuse_native_input) THEN
1697 158 : IF (seed_native_from_rho) THEN
1698 : ! SCF_GUESS/RESTART controls the initial SCC state.
1699 140 : CALL tb_set_population(tb, ch_shell_spin, ao_dip_spin, ao_quad_spin)
1700 : ELSE
1701 262 : tb%wfn%qsh(:, :) = 0.0_dp
1702 114 : tb%wfn%qat(:, :) = 0.0_dp
1703 280 : IF (do_dipole) tb%wfn%dpat(:, :, :) = 0.0_dp
1704 466 : IF (do_quadrupole) tb%wfn%qpat(:, :, :) = 0.0_dp
1705 : END IF
1706 : END IF
1707 : ELSE
1708 : ! Record input/output pairs only after diagonalization, never for forces.
1709 22752 : IF (.NOT. calculate_forces) THEN
1710 22608 : CALL tb%mixer%set(tb%wfn%qsh)
1711 22608 : IF (do_dipole) CALL tb%mixer%set(tb%wfn%dpat)
1712 22608 : IF (do_quadrupole) CALL tb%mixer%set(tb%wfn%qpat)
1713 : END IF
1714 22752 : CALL tb_set_population(tb, ch_shell_spin, ao_dip_spin, ao_quad_spin)
1715 22752 : IF (.NOT. calculate_forces) THEN
1716 22608 : CALL tb%mixer%diff(tb%wfn%qsh)
1717 22608 : IF (do_dipole) CALL tb%mixer%diff(tb%wfn%dpat)
1718 22608 : IF (do_quadrupole) CALL tb%mixer%diff(tb%wfn%qpat)
1719 : END IF
1720 : END IF
1721 : ELSE
1722 6212 : IF (.NOT. (dft_control%qs_control%do_ls_scf .OR. scf_control%use_ot)) THEN
1723 4460 : IF (.NOT. ASSOCIATED(scf_env)) CPABORT("CP2K SCC mixer requires a QS SCF environment")
1724 4460 : IF (.NOT. (do_dipole .OR. do_quadrupole) .AND. &
1725 : scf_env%mixing_method == modified_broyden_mixing_nr) THEN
1726 0 : CPABORT("MODIFIED_BROYDEN_MIXING with SCC_MIXER CP2K requires GFN2")
1727 : END IF
1728 : END IF
1729 6212 : IF (.NOT. use_no_mixer) THEN
1730 4460 : CALL tb_mix_population(scf_env, para_env, ch_shell_spin, ao_dip_spin, ao_quad_spin)
1731 : END IF
1732 6212 : CALL tb_set_population(tb, ch_shell_spin, ao_dip_spin, ao_quad_spin)
1733 6212 : IF (nspin == 1) THEN
1734 : IF (.NOT. (do_dipole .OR. do_quadrupole) .OR. &
1735 4584 : dft_control%qs_control%do_ls_scf .OR. scf_control%use_ot) THEN
1736 : ! Preserve the existing atom-charge convention for shell-only mixing
1737 : ! and direct density optimization: qat comes from the unmixed density.
1738 12238 : tb%wfn%qat(:, 1) = -ch_atom(:, 1)
1739 : END IF
1740 : END IF
1741 : END IF
1742 :
1743 53808 : CALL tb_update_scf_potential(qs_env, tb, para_env, calculate_forces)
1744 :
1745 53808 : IF (ALLOCATED(ao_dip_spin)) DEALLOCATE (ao_dip_spin)
1746 53808 : IF (ALLOCATED(ao_quad_spin)) DEALLOCATE (ao_quad_spin)
1747 53808 : DEALLOCATE (ch_atom, ch_orb_spin, ch_shell_spin)
1748 :
1749 : #else
1750 : MARK_USED(qs_env)
1751 : MARK_USED(tb)
1752 : MARK_USED(dft_control)
1753 : MARK_USED(calculate_forces)
1754 : MARK_USED(use_rho)
1755 : CPABORT("Built without TBLITE")
1756 : #endif
1757 :
1758 161424 : END SUBROUTINE tb_update_charges
1759 :
1760 : #if defined(__TBLITE)
1761 : ! **************************************************************************************************
1762 : !> \brief Convert CP2K shell/multipole populations to tblite charge/magnetization state.
1763 : !> \param tb ...
1764 : !> \param shell Electron excess per atom, shell and spin channel.
1765 : !> \param dipole Electron dipoles, allocated only when required.
1766 : !> \param quadrupole Electron quadrupoles, allocated only when required.
1767 : ! **************************************************************************************************
1768 29104 : SUBROUTINE tb_set_population(tb, shell, dipole, quadrupole)
1769 : TYPE(tblite_type), INTENT(INOUT) :: tb
1770 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: shell
1771 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1772 : INTENT(IN) :: dipole, quadrupole
1773 :
1774 : INTEGER :: iatom, ii, is, ispin
1775 :
1776 165136 : DO iatom = 1, tb%mol%nat
1777 136032 : ii = tb%calc%bas%ish_at(iatom)
1778 309932 : DO ispin = 1, tb%wfn%nspin
1779 516864 : DO is = 1, tb%calc%bas%nsh_at(iatom)
1780 516864 : tb%wfn%qsh(ii + is, ispin) = -shell(iatom, is, ispin)
1781 : END DO
1782 448228 : IF (ALLOCATED(dipole)) tb%wfn%dpat(:, iatom, ispin) = -dipole(iatom, :, ispin)
1783 887692 : IF (ALLOCATED(quadrupole)) tb%wfn%qpat(:, iatom, ispin) = -quadrupole(iatom, :, ispin)
1784 : END DO
1785 : END DO
1786 29104 : CALL tb_sum_shell_charges(tb)
1787 :
1788 29104 : END SUBROUTINE tb_set_population
1789 :
1790 : ! **************************************************************************************************
1791 : !> \brief Reconstruct tblite atomic charges from its current shell charges.
1792 : !> \param tb ...
1793 : ! **************************************************************************************************
1794 51490 : SUBROUTINE tb_sum_shell_charges(tb)
1795 : TYPE(tblite_type), INTENT(INOUT) :: tb
1796 :
1797 : INTEGER :: iatom, ii, ispin
1798 :
1799 298134 : DO iatom = 1, tb%mol%nat
1800 246644 : ii = tb%calc%bas%ish_at(iatom)
1801 558926 : DO ispin = 1, tb%wfn%nspin
1802 1186900 : tb%wfn%qat(iatom, ispin) = SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), ispin))
1803 : END DO
1804 : END DO
1805 :
1806 51490 : END SUBROUTINE tb_sum_shell_charges
1807 :
1808 : ! **************************************************************************************************
1809 : !> \brief Apply the CP2K mixer to shell and multipole populations for all spin channels.
1810 : !> \param scf_env ...
1811 : !> \param para_env ...
1812 : !> \param shell ...
1813 : !> \param dipole ...
1814 : !> \param quadrupole ...
1815 : ! **************************************************************************************************
1816 4460 : SUBROUTINE tb_mix_population(scf_env, para_env, shell, dipole, quadrupole)
1817 : TYPE(qs_scf_env_type), POINTER :: scf_env
1818 : TYPE(mp_para_env_type), POINTER :: para_env
1819 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: shell
1820 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1821 : INTENT(INOUT) :: dipole, quadrupole
1822 :
1823 : INTEGER :: ispin, ncols, nsh, nspin, offset
1824 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mix_vars
1825 :
1826 4460 : nspin = SIZE(shell, 3)
1827 4460 : IF (nspin == 1 .AND. .NOT. ALLOCATED(dipole) .AND. .NOT. ALLOCATED(quadrupole)) THEN
1828 : ! Retain CP2K's electron-excess convention for single-spin shell-only mixing.
1829 : CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
1830 0 : shell(:, :, 1), para_env, scf_env%iter_count)
1831 : RETURN
1832 : END IF
1833 :
1834 4460 : nsh = SIZE(shell, 2)
1835 4460 : ncols = nspin*nsh
1836 4460 : IF (ALLOCATED(dipole)) ncols = ncols + nspin*dip_n
1837 4460 : IF (ALLOCATED(quadrupole)) ncols = ncols + nspin*quad_n
1838 17840 : ALLOCATE (mix_vars(SIZE(shell, 1), ncols))
1839 :
1840 : ! Preserve the mixer history layout: all shell channels, then dipoles, then quadrupoles.
1841 4460 : offset = 0
1842 10490 : DO ispin = 1, nspin
1843 47866 : mix_vars(:, offset + 1:offset + nsh) = -shell(:, :, ispin)
1844 10490 : offset = offset + nsh
1845 : END DO
1846 4460 : IF (ALLOCATED(dipole)) THEN
1847 10490 : DO ispin = 1, nspin
1848 68784 : mix_vars(:, offset + 1:offset + dip_n) = -dipole(:, :, ispin)
1849 10490 : offset = offset + dip_n
1850 : END DO
1851 : END IF
1852 4460 : IF (ALLOCATED(quadrupole)) THEN
1853 10490 : DO ispin = 1, nspin
1854 131538 : mix_vars(:, offset + 1:offset + quad_n) = -quadrupole(:, :, ispin)
1855 10490 : offset = offset + quad_n
1856 : END DO
1857 : END IF
1858 : CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
1859 4460 : mix_vars, para_env, scf_env%iter_count)
1860 4460 : offset = 0
1861 10490 : DO ispin = 1, nspin
1862 47866 : shell(:, :, ispin) = -mix_vars(:, offset + 1:offset + nsh)
1863 10490 : offset = offset + nsh
1864 : END DO
1865 4460 : IF (ALLOCATED(dipole)) THEN
1866 10490 : DO ispin = 1, nspin
1867 68784 : dipole(:, :, ispin) = -mix_vars(:, offset + 1:offset + dip_n)
1868 10490 : offset = offset + dip_n
1869 : END DO
1870 : END IF
1871 4460 : IF (ALLOCATED(quadrupole)) THEN
1872 10490 : DO ispin = 1, nspin
1873 131538 : quadrupole(:, :, ispin) = -mix_vars(:, offset + 1:offset + quad_n)
1874 10490 : offset = offset + quad_n
1875 : END DO
1876 : END IF
1877 :
1878 4460 : END SUBROUTINE tb_mix_population
1879 :
1880 : ! **************************************************************************************************
1881 : !> \brief Evaluate SCC potentials, energies and optional gradients for the current tblite state.
1882 : !> \param qs_env ...
1883 : !> \param tb ...
1884 : !> \param para_env ...
1885 : !> \param calculate_forces ...
1886 : ! **************************************************************************************************
1887 53808 : SUBROUTINE tb_update_scf_potential(qs_env, tb, para_env, calculate_forces)
1888 : TYPE(qs_environment_type), POINTER :: qs_env
1889 : TYPE(tblite_type), POINTER :: tb
1890 : TYPE(mp_para_env_type), POINTER :: para_env
1891 : LOGICAL, INTENT(IN) :: calculate_forces
1892 :
1893 : ! These library results are replicated across MPI ranks. The force conversion
1894 : ! retains the existing normalization before CP2K's global force reduction.
1895 53808 : CALL tb%pot%reset
1896 310558 : tb%e_es = 0.0_dp
1897 310558 : tb%e_scd = 0.0_dp
1898 310558 : tb%e_int = 0.0_dp
1899 53808 : IF (ALLOCATED(tb%calc%coulomb)) THEN
1900 53808 : CALL tb%calc%coulomb%get_potential(tb%mol, tb%cache, tb%wfn, tb%pot)
1901 53808 : CALL tb%calc%coulomb%get_energy(tb%mol, tb%cache, tb%wfn, tb%e_es)
1902 : END IF
1903 53808 : IF (ALLOCATED(tb%calc%dispersion)) THEN
1904 53808 : CALL tb%calc%dispersion%get_potential(tb%mol, tb%dcache, tb%wfn, tb%pot)
1905 53808 : CALL tb%calc%dispersion%get_energy(tb%mol, tb%dcache, tb%wfn, tb%e_scd)
1906 : END IF
1907 53808 : IF (ALLOCATED(tb%calc%interactions)) THEN
1908 5502 : CALL tb%calc%interactions%get_potential(tb%mol, tb%icache, tb%wfn, tb%pot)
1909 5502 : CALL tb%calc%interactions%get_energy(tb%mol, tb%icache, tb%wfn, tb%e_int)
1910 : END IF
1911 :
1912 53808 : IF (calculate_forces) THEN
1913 154 : IF (ALLOCATED(tb%calc%coulomb)) THEN
1914 3754 : tb%grad = 0.0_dp
1915 154 : CALL tb%calc%coulomb%get_gradient(tb%mol, tb%cache, tb%wfn, tb%grad, tb%sigma)
1916 154 : CALL tb_grad2force(qs_env, tb, para_env, 3)
1917 : END IF
1918 :
1919 154 : IF (ALLOCATED(tb%calc%dispersion)) THEN
1920 3754 : tb%grad = 0.0_dp
1921 154 : CALL tb%calc%dispersion%get_gradient(tb%mol, tb%dcache, tb%wfn, tb%grad, tb%sigma)
1922 154 : CALL tb_grad2force(qs_env, tb, para_env, 2)
1923 : END IF
1924 :
1925 154 : IF (ALLOCATED(tb%calc%interactions)) THEN
1926 230 : tb%grad = 0.0_dp
1927 22 : CALL tb%calc%interactions%get_gradient(tb%mol, tb%icache, tb%wfn, tb%grad, tb%sigma)
1928 22 : CALL tb_grad2force(qs_env, tb, para_env, 3)
1929 : END IF
1930 : END IF
1931 :
1932 53808 : END SUBROUTINE tb_update_scf_potential
1933 : #endif
1934 :
1935 : ! **************************************************************************************************
1936 : !> \brief ...
1937 : !> \param qs_env ...
1938 : !> \param tb ...
1939 : !> \param dft_control ...
1940 : ! **************************************************************************************************
1941 28208 : SUBROUTINE tb_ham_add_coulomb(qs_env, tb, dft_control)
1942 :
1943 : TYPE(qs_environment_type), POINTER :: qs_env
1944 : TYPE(tblite_type), POINTER :: tb
1945 : TYPE(dft_control_type), POINTER :: dft_control
1946 :
1947 : #if defined(__TBLITE)
1948 : INTEGER :: i, im, ispin, ikind, iatom, jatom, irow, icol, img, slot
1949 : INTEGER :: nimg, nspin, nsgf, natorb, natom_kind, ndip, nquad
1950 : INTEGER, DIMENSION(3) :: cell
1951 : INTEGER, DIMENSION(25) :: nao
1952 28208 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1953 : LOGICAL :: found
1954 28208 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sblock
1955 28208 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: vat
1956 28208 : REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: vsh, vmp
1957 28208 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1958 : TYPE(dbcsr_iterator_type) :: iter
1959 28208 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, ks_matrix
1960 : TYPE(kpoint_type), POINTER :: kpoints
1961 28208 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: n_list
1962 28208 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1963 : TYPE(xtb_atom_type), POINTER :: xtb_atom
1964 : !$ INTEGER :: hash, lock_num
1965 : !$ INTEGER(KIND=int_8) :: pair_index
1966 28208 : !$ INTEGER(KIND=omp_lock_kind), DIMENSION(:), ALLOCATABLE :: locks
1967 : INTEGER, PARAMETER :: nlock = 501
1968 :
1969 28208 : nimg = dft_control%nimages
1970 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1971 28208 : matrix_s_kp=matrix_s, matrix_ks_kp=ks_matrix)
1972 28208 : nspin = SIZE(ks_matrix, 1)
1973 28208 : ndip = MERGE(dip_n, 0, ASSOCIATED(tb%dipbra))
1974 28208 : nquad = MERGE(quad_n, 0, ASSOCIATED(tb%quadbra))
1975 28208 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
1976 329654 : ALLOCATE (vsh(nsgf, tb%mol%nat, nspin), vat(tb%mol%nat, nspin), vmp(ndip + nquad, tb%mol%nat, nspin))
1977 :
1978 : ! Project each potential once per atom and spin, outside the matrix-block loop.
1979 28208 : vsh = 0.0_dp
1980 71440 : DO ikind = 1, SIZE(atomic_kind_set)
1981 43232 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_kind)
1982 43232 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom)
1983 43232 : CALL get_xtb_atom_param(xtb_atom, nao=nao, natorb=natorb)
1984 248536 : DO i = 1, natom_kind
1985 133864 : iatom = atomic_kind_set(ikind)%atom_list(i)
1986 318848 : DO ispin = 1, nspin
1987 1124618 : DO im = 1, natorb
1988 1124618 : vsh(im, iatom, ispin) = tb_spin_project(tb%pot%vsh(nao(im) + tb%calc%bas%ish_at(iatom), :), ispin)
1989 : END DO
1990 141752 : vat(iatom, ispin) = tb_spin_project(tb%pot%vat(iatom, :), ispin)
1991 428570 : DO im = 1, ndip
1992 428570 : vmp(im, iatom, ispin) = tb_spin_project(tb%pot%vdp(im, iatom, :), ispin)
1993 : END DO
1994 849252 : DO im = 1, nquad
1995 715388 : vmp(ndip + im, iatom, ispin) = tb_spin_project(tb%pot%vqp(im, iatom, :), ispin)
1996 : END DO
1997 : END DO
1998 : END DO
1999 : END DO
2000 :
2001 28208 : IF (nimg == 1) THEN
2002 15248 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tb, matrix_s, ks_matrix, vsh, vat, vmp) PRIVATE(iter, irow, icol, sblock)
2003 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
2004 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2005 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
2006 : CALL tb_ham_add_block(tb, ks_matrix, sblock, irow, icol, 1, vsh, vat, vmp)
2007 : END DO
2008 : CALL dbcsr_iterator_stop(iter)
2009 : !$OMP END PARALLEL
2010 : ELSE
2011 12960 : CALL get_qs_env(qs_env=qs_env, sab_kp=n_list, kpoints=kpoints)
2012 12960 : IF (.NOT. ASSOCIATED(n_list)) THEN
2013 0 : CPABORT("Missing k-point neighbor list for tblite Hamiltonian")
2014 : END IF
2015 12960 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
2016 12960 : !$ ALLOCATE (locks(nlock))
2017 6505920 : !$ DO lock_num = 1, nlock
2018 6505920 : !$ CALL omp_init_lock(locks(lock_num))
2019 : !$ END DO
2020 : ! Preserve neighbor-list multiplicities when several lattice cells share a block.
2021 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(GUIDED) &
2022 : !$OMP SHARED(tb, n_list, cell_to_index, matrix_s, ks_matrix, vsh, vat, vmp, locks) &
2023 12960 : !$OMP PRIVATE(slot, iatom, jatom, irow, icol, cell, img, sblock, found, pair_index, hash)
2024 : DO slot = 1, n_list(1)%nl_size
2025 : iatom = n_list(1)%nlist_task(slot)%iatom
2026 : jatom = n_list(1)%nlist_task(slot)%jatom
2027 : cell(:) = n_list(1)%nlist_task(slot)%cell(:)
2028 : irow = MIN(iatom, jatom)
2029 : icol = MAX(iatom, jatom)
2030 : img = cell_to_index(cell(1), cell(2), cell(3))
2031 : CPASSERT(img > 0)
2032 : CALL dbcsr_get_block_p(matrix_s(1, img)%matrix, irow, icol, BLOCK=sblock, found=found)
2033 : CPASSERT(found)
2034 : !$ pair_index = INT(irow - 1, int_8)*INT(tb%mol%nat, int_8) + INT(icol, int_8)
2035 : !$ hash = INT(MOD(pair_index, INT(nlock, int_8))) + 1
2036 : !$ CALL omp_set_lock(locks(hash))
2037 : CALL tb_ham_add_block(tb, ks_matrix, sblock, irow, icol, img, vsh, vat, vmp)
2038 : !$ CALL omp_unset_lock(locks(hash))
2039 : END DO
2040 : !$OMP END PARALLEL DO
2041 6505920 : !$ DO lock_num = 1, nlock
2042 6505920 : !$ CALL omp_destroy_lock(locks(lock_num))
2043 : !$ END DO
2044 12960 : !$ DEALLOCATE (locks)
2045 : END IF
2046 28208 : DEALLOCATE (vsh, vat, vmp)
2047 : #else
2048 : MARK_USED(qs_env)
2049 : MARK_USED(tb)
2050 : MARK_USED(dft_control)
2051 : CPABORT("Built without TBLITE")
2052 : #endif
2053 :
2054 84624 : END SUBROUTINE tb_ham_add_coulomb
2055 :
2056 : #if defined(__TBLITE)
2057 : ! **************************************************************************************************
2058 : !> \brief Apply shell, atomic and multipole potentials to one matrix block.
2059 : !> \param tb ...
2060 : !> \param ks_matrix ...
2061 : !> \param sblock ...
2062 : !> \param irow ...
2063 : !> \param icol ...
2064 : !> \param img ...
2065 : !> \param vsh ...
2066 : !> \param vat ...
2067 : !> \param vmp ...
2068 : ! **************************************************************************************************
2069 2790085 : SUBROUTINE tb_ham_add_block(tb, ks_matrix, sblock, irow, icol, img, vsh, vat, vmp)
2070 : TYPE(tblite_type), POINTER :: tb
2071 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
2072 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sblock
2073 : INTEGER, INTENT(IN) :: irow, icol, img
2074 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: vsh
2075 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: vat
2076 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: vmp
2077 :
2078 : INTEGER :: i, im, ispin, j
2079 : LOGICAL :: found
2080 : REAL(KIND=dp) :: vatom
2081 2790085 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksblock
2082 55801700 : TYPE(block_p_type), DIMENSION(dip_n+quad_n) :: bra, ket
2083 :
2084 2790085 : CALL tb_get_multipole_blocks(tb, irow, icol, img, bra, ket)
2085 5604985 : DO ispin = 1, SIZE(ks_matrix, 1)
2086 2814900 : CALL dbcsr_get_block_p(ks_matrix(ispin, img)%matrix, irow, icol, BLOCK=ksblock, found=found)
2087 2814900 : CPASSERT(found)
2088 2814900 : vatom = vat(irow, ispin) + vat(icol, ispin)
2089 22831875 : DO j = 1, SIZE(sblock, 2)
2090 189649234 : DO i = 1, SIZE(sblock, 1)
2091 : ksblock(i, j) = ksblock(i, j) - 0.5_dp*(vsh(i, irow, ispin)*sblock(i, j) &
2092 166817359 : + sblock(i, j)*vsh(j, icol, ispin))
2093 186834334 : ksblock(i, j) = ksblock(i, j) - 0.5_dp*vatom*sblock(i, j)
2094 : END DO
2095 : END DO
2096 19524913 : DO im = 1, SIZE(vmp, 1)
2097 : ksblock(:, :) = ksblock(:, :) - 0.5_dp*(ket(im)%block(:, :)*vmp(im, irow, ispin) &
2098 853586142 : + bra(im)%block(:, :)*vmp(im, icol, ispin))
2099 : END DO
2100 : END DO
2101 2790085 : END SUBROUTINE tb_ham_add_block
2102 : #endif
2103 :
2104 : ! **************************************************************************************************
2105 : !> \brief ...
2106 : !> \param qs_env ...
2107 : !> \param tb ...
2108 : ! **************************************************************************************************
2109 1382 : SUBROUTINE tb_get_multipole(qs_env, tb)
2110 :
2111 : TYPE(qs_environment_type), POINTER :: qs_env
2112 : TYPE(tblite_type), POINTER :: tb
2113 :
2114 : #if defined(__TBLITE)
2115 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tb_get_multipole'
2116 :
2117 : INTEGER :: handle, nimg, maxsgf, img, i, idx, slot, ic, im
2118 : INTEGER :: ikind, jkind, iatom, jatom, irow, icol, ityp, jtyp
2119 : INTEGER :: iset, jset, inda, indb, ia, ib, ij, ni, nj
2120 : INTEGER, DIMENSION(3) :: cell
2121 1382 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2122 : REAL(KIND=dp) :: r2
2123 : REAL(KIND=dp), DIMENSION(3) :: rij
2124 1382 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: stmp
2125 2764 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: dtmp, qtmp, dtmpj, qtmpj
2126 1382 : REAL(KIND=dp), DIMENSION(:, :, :, :), ALLOCATABLE :: pair_multipole
2127 27640 : TYPE(block_p_type), DIMENSION(dip_n + quad_n) :: bra, ket
2128 1382 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
2129 : TYPE(dft_control_type), POINTER :: dft_control
2130 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
2131 1382 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
2132 : TYPE(kpoint_type), POINTER :: kpoints
2133 1382 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb, sab_kp
2134 1382 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2135 : !$ INTEGER :: hash, lock_num
2136 : !$ INTEGER(KIND=int_8) :: pair_index
2137 1382 : !$ INTEGER(KIND=omp_lock_kind), DIMENSION(:), ALLOCATABLE :: locks
2138 : INTEGER, PARAMETER :: nlock = 501
2139 :
2140 1382 : CALL timeset(routineN, handle)
2141 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, sab_orb=sab_orb, sab_kp=sab_kp, &
2142 1382 : dft_control=dft_control, kpoints=kpoints, matrix_s_kp=matrix_s)
2143 1382 : nimg = dft_control%nimages
2144 1382 : IF (nimg > 1) THEN
2145 586 : IF (.NOT. ASSOCIATED(sab_kp)) THEN
2146 0 : CPABORT("Missing k-point neighbor list for tblite multipoles")
2147 : END IF
2148 586 : sab_orb => sab_kp
2149 586 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
2150 : END IF
2151 1382 : CALL get_qs_kind_set(qs_kind_set, maxsgf=maxsgf)
2152 6194 : ALLOCATE (basis_set_list(SIZE(qs_kind_set)))
2153 1382 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
2154 :
2155 : ! allocate dipole/quadrupole moment matrix elemnts
2156 1382 : CALL dbcsr_allocate_matrix_set(tb%dipbra, dip_n*nimg)
2157 1382 : CALL dbcsr_allocate_matrix_set(tb%dipket, dip_n*nimg)
2158 1382 : CALL dbcsr_allocate_matrix_set(tb%quadbra, quad_n*nimg)
2159 1382 : CALL dbcsr_allocate_matrix_set(tb%quadket, quad_n*nimg)
2160 26084 : DO img = 1, nimg
2161 98808 : DO i = 1, dip_n
2162 74106 : idx = i + dip_n*(img - 1)
2163 74106 : ALLOCATE (tb%dipbra(idx)%matrix)
2164 74106 : ALLOCATE (tb%dipket(idx)%matrix)
2165 : CALL dbcsr_create(tb%dipbra(idx)%matrix, template=matrix_s(1, img)%matrix, &
2166 74106 : name="DIPOLE BRAMATRIX")
2167 : CALL dbcsr_create(tb%dipket(idx)%matrix, template=matrix_s(1, img)%matrix, &
2168 74106 : name="DIPOLE KETMATRIX")
2169 74106 : CALL cp_dbcsr_alloc_block_from_nbl(tb%dipbra(idx)%matrix, sab_orb)
2170 98808 : CALL cp_dbcsr_alloc_block_from_nbl(tb%dipket(idx)%matrix, sab_orb)
2171 : END DO
2172 174296 : DO i = 1, quad_n
2173 148212 : idx = i + quad_n*(img - 1)
2174 148212 : ALLOCATE (tb%quadbra(idx)%matrix)
2175 148212 : ALLOCATE (tb%quadket(idx)%matrix)
2176 : CALL dbcsr_create(tb%quadbra(idx)%matrix, template=matrix_s(1, img)%matrix, &
2177 148212 : name="QUADRUPOLE BRAMATRIX")
2178 : CALL dbcsr_create(tb%quadket(idx)%matrix, template=matrix_s(1, img)%matrix, &
2179 148212 : name="QUADRUPOLE KETMATRIX")
2180 148212 : CALL cp_dbcsr_alloc_block_from_nbl(tb%quadbra(idx)%matrix, sab_orb)
2181 172914 : CALL cp_dbcsr_alloc_block_from_nbl(tb%quadket(idx)%matrix, sab_orb)
2182 : END DO
2183 : END DO
2184 :
2185 1382 : !$ ALLOCATE (locks(nlock))
2186 693764 : !$ DO lock_num = 1, nlock
2187 693764 : !$ CALL omp_init_lock(locks(lock_num))
2188 : !$ END DO
2189 :
2190 : !$OMP PARALLEL DEFAULT(NONE) &
2191 : !$OMP SHARED(tb, sab_orb, basis_set_list, nimg, cell_to_index, maxsgf, locks) &
2192 : !$OMP PRIVATE(stmp, dtmp, qtmp, dtmpj, qtmpj, pair_multipole, slot, ikind, jkind, iatom, jatom, cell, rij, r2, &
2193 : !$OMP irow, icol, i, ic, ityp, jtyp, basis_set_a, basis_set_b, iset, jset, inda, indb, ia, ib, ij, &
2194 1382 : !$OMP bra, ket, im, ni, nj, pair_index, hash)
2195 : ALLOCATE (stmp(msao(tb%calc%bas%maxl)**2))
2196 : ALLOCATE (dtmp(dip_n, SIZE(stmp)), qtmp(quad_n, SIZE(stmp)))
2197 : ALLOCATE (dtmpj(dip_n, SIZE(stmp)), qtmpj(quad_n, SIZE(stmp)))
2198 : ALLOCATE (pair_multipole(maxsgf, maxsgf, dip_n + quad_n, 2))
2199 : !$OMP DO SCHEDULE(GUIDED)
2200 : DO slot = 1, sab_orb(1)%nl_size
2201 : ikind = sab_orb(1)%nlist_task(slot)%ikind
2202 : jkind = sab_orb(1)%nlist_task(slot)%jkind
2203 : iatom = sab_orb(1)%nlist_task(slot)%iatom
2204 : jatom = sab_orb(1)%nlist_task(slot)%jatom
2205 : cell(:) = sab_orb(1)%nlist_task(slot)%cell(:)
2206 : rij(:) = sab_orb(1)%nlist_task(slot)%r(:)
2207 : r2 = NORM2(rij)**2
2208 : irow = MIN(iatom, jatom)
2209 : icol = MAX(iatom, jatom)
2210 : IF (iatom < jatom) THEN
2211 : rij = -rij
2212 : i = ikind
2213 : ikind = jkind
2214 : jkind = i
2215 : END IF
2216 : IF (nimg == 1) THEN
2217 : ic = 1
2218 : ELSE
2219 : ic = cell_to_index(cell(1), cell(2), cell(3))
2220 : CPASSERT(ic > 0)
2221 : END IF
2222 : ityp = tb%mol%id(icol)
2223 : jtyp = tb%mol%id(irow)
2224 : basis_set_a => basis_set_list(ikind)%gto_basis_set
2225 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
2226 : basis_set_b => basis_set_list(jkind)%gto_basis_set
2227 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
2228 : CALL tb_get_multipole_blocks(tb, irow, icol, ic, bra, ket)
2229 : ni = SIZE(bra(1)%block, 1)
2230 : nj = SIZE(bra(1)%block, 2)
2231 : pair_multipole(1:ni, 1:nj, :, :) = 0.0_dp
2232 :
2233 : DO iset = 1, basis_set_a%nset
2234 : DO jset = 1, basis_set_b%nset
2235 : CALL multipole_cgto(tb%calc%bas%cgto(jset, jtyp), tb%calc%bas%cgto(iset, ityp), &
2236 : r2, -rij, tb%calc%bas%intcut, stmp, dtmp, qtmp)
2237 : DO inda = 1, basis_set_a%nsgf_set(iset)
2238 : ia = basis_set_a%first_sgf(1, iset) - basis_set_a%first_sgf(1, 1) + inda
2239 : DO indb = 1, basis_set_b%nsgf_set(jset)
2240 : ib = basis_set_b%first_sgf(1, jset) - basis_set_b%first_sgf(1, 1) + indb
2241 : ij = indb + basis_set_b%nsgf_set(jset)*(inda - 1)
2242 : ! Periodic self-images still require the shift to the other center.
2243 : IF (icol == irow .AND. r2 < same_atom**2) THEN
2244 : dtmpj(:, ij) = dtmp(:, ij)
2245 : qtmpj(:, ij) = qtmp(:, ij)
2246 : ELSE
2247 : CALL tb_shift_multipole(-rij, stmp(ij), dtmp(:, ij), qtmp(:, ij), &
2248 : dtmpj(:, ij), qtmpj(:, ij))
2249 : END IF
2250 : pair_multipole(ib, ia, 1:dip_n, 1) = pair_multipole(ib, ia, 1:dip_n, 1) + dtmp(:, ij)
2251 : pair_multipole(ib, ia, dip_n + 1:, 1) = pair_multipole(ib, ia, dip_n + 1:, 1) + qtmp(:, ij)
2252 : pair_multipole(ib, ia, 1:dip_n, 2) = pair_multipole(ib, ia, 1:dip_n, 2) + dtmpj(:, ij)
2253 : pair_multipole(ib, ia, dip_n + 1:, 2) = pair_multipole(ib, ia, dip_n + 1:, 2) + qtmpj(:, ij)
2254 : END DO
2255 : END DO
2256 : END DO
2257 : END DO
2258 :
2259 : ! Integrals are thread-local; lock only the accumulation into shared blocks.
2260 : !$ pair_index = INT(irow - 1, int_8)*INT(tb%mol%nat, int_8) + INT(icol, int_8)
2261 : !$ hash = INT(MOD(pair_index, INT(nlock, int_8))) + 1
2262 : !$ CALL omp_set_lock(locks(hash))
2263 : DO im = 1, dip_n + quad_n
2264 : bra(im)%block(:, :) = bra(im)%block(:, :) + pair_multipole(1:ni, 1:nj, im, 1)
2265 : ket(im)%block(:, :) = ket(im)%block(:, :) + pair_multipole(1:ni, 1:nj, im, 2)
2266 : END DO
2267 : !$ CALL omp_unset_lock(locks(hash))
2268 : END DO
2269 : !$OMP END DO
2270 : DEALLOCATE (stmp, dtmp, qtmp, dtmpj, qtmpj, pair_multipole)
2271 : !$OMP END PARALLEL
2272 693764 : !$ DO lock_num = 1, nlock
2273 693764 : !$ CALL omp_destroy_lock(locks(lock_num))
2274 : !$ END DO
2275 1382 : !$ DEALLOCATE (locks)
2276 :
2277 75488 : DO i = 1, SIZE(tb%dipbra)
2278 74106 : CALL dbcsr_finalize(tb%dipbra(i)%matrix)
2279 75488 : CALL dbcsr_finalize(tb%dipket(i)%matrix)
2280 : END DO
2281 149594 : DO i = 1, SIZE(tb%quadbra)
2282 148212 : CALL dbcsr_finalize(tb%quadbra(i)%matrix)
2283 149594 : CALL dbcsr_finalize(tb%quadket(i)%matrix)
2284 : END DO
2285 1382 : DEALLOCATE (basis_set_list)
2286 1382 : CALL timestop(handle)
2287 : #else
2288 : MARK_USED(qs_env)
2289 : MARK_USED(tb)
2290 : CPABORT("Built without TBLITE")
2291 : #endif
2292 :
2293 6910 : END SUBROUTINE tb_get_multipole
2294 :
2295 : ! **************************************************************************************************
2296 : !> \brief Shift a multipole operator from one center to the other.
2297 : !> \param vec displacement vector between the two centers
2298 : !> \param s overlap integral
2299 : !> \param di dipole integral on the original center
2300 : !> \param qi quadrupole integral on the original center
2301 : !> \param dj dipole integral on the shifted center
2302 : !> \param qj quadrupole integral on the shifted center
2303 : ! **************************************************************************************************
2304 20169594 : PURE SUBROUTINE tb_shift_multipole(vec, s, di, qi, dj, qj)
2305 :
2306 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vec
2307 : REAL(KIND=dp), INTENT(IN) :: s
2308 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: di, qi
2309 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: dj, qj
2310 :
2311 : REAL(KIND=dp) :: tr
2312 :
2313 20169594 : dj(1) = di(1) + vec(1)*s
2314 20169594 : dj(2) = di(2) + vec(2)*s
2315 20169594 : dj(3) = di(3) + vec(3)*s
2316 :
2317 20169594 : qj(1) = 2*vec(1)*di(1) + vec(1)**2*s
2318 20169594 : qj(3) = 2*vec(2)*di(2) + vec(2)**2*s
2319 20169594 : qj(6) = 2*vec(3)*di(3) + vec(3)**2*s
2320 20169594 : qj(2) = vec(1)*di(2) + vec(2)*di(1) + vec(1)*vec(2)*s
2321 20169594 : qj(4) = vec(1)*di(3) + vec(3)*di(1) + vec(1)*vec(3)*s
2322 20169594 : qj(5) = vec(2)*di(3) + vec(3)*di(2) + vec(2)*vec(3)*s
2323 20169594 : tr = 0.5_dp*(qj(1) + qj(3) + qj(6))
2324 :
2325 20169594 : qj(1) = qi(1) + 1.5_dp*qj(1) - tr
2326 20169594 : qj(2) = qi(2) + 1.5_dp*qj(2)
2327 20169594 : qj(3) = qi(3) + 1.5_dp*qj(3) - tr
2328 20169594 : qj(4) = qi(4) + 1.5_dp*qj(4)
2329 20169594 : qj(5) = qi(5) + 1.5_dp*qj(5)
2330 20169594 : qj(6) = qi(6) + 1.5_dp*qj(6) - tr
2331 :
2332 20169594 : END SUBROUTINE tb_shift_multipole
2333 :
2334 : #if defined(__TBLITE)
2335 : ! **************************************************************************************************
2336 : !> \brief Fetch the multipole blocks in dipole, then quadrupole order.
2337 : !> \param tb ...
2338 : !> \param irow ...
2339 : !> \param icol ...
2340 : !> \param img ...
2341 : !> \param bra ...
2342 : !> \param ket ...
2343 : ! **************************************************************************************************
2344 324060600 : SUBROUTINE tb_get_multipole_blocks(tb, irow, icol, img, bra, ket)
2345 : TYPE(tblite_type), POINTER :: tb
2346 : INTEGER, INTENT(IN) :: irow, icol, img
2347 : TYPE(block_p_type), DIMENSION(:), INTENT(OUT) :: bra, ket
2348 :
2349 : INTEGER :: idx, im, ndip
2350 : LOGICAL :: found
2351 :
2352 16203030 : ndip = 0
2353 16203030 : IF (ASSOCIATED(tb%dipbra)) THEN
2354 29419912 : ndip = dip_n
2355 29419912 : DO im = 1, dip_n
2356 22064934 : idx = im + dip_n*(img - 1)
2357 22064934 : CALL dbcsr_get_block_p(tb%dipbra(idx)%matrix, irow, icol, BLOCK=bra(im)%block, found=found)
2358 22064934 : CPASSERT(found)
2359 22064934 : CALL dbcsr_get_block_p(tb%dipket(idx)%matrix, irow, icol, BLOCK=ket(im)%block, found=found)
2360 51484846 : CPASSERT(found)
2361 : END DO
2362 : END IF
2363 16203030 : IF (ASSOCIATED(tb%quadbra)) THEN
2364 51484846 : DO im = 1, quad_n
2365 44129868 : idx = im + quad_n*(img - 1)
2366 44129868 : CALL dbcsr_get_block_p(tb%quadbra(idx)%matrix, irow, icol, BLOCK=bra(ndip + im)%block, found=found)
2367 44129868 : CPASSERT(found)
2368 44129868 : CALL dbcsr_get_block_p(tb%quadket(idx)%matrix, irow, icol, BLOCK=ket(ndip + im)%block, found=found)
2369 95614714 : CPASSERT(found)
2370 : END DO
2371 : END IF
2372 16203030 : END SUBROUTINE tb_get_multipole_blocks
2373 :
2374 : ! **************************************************************************************************
2375 : !> \brief Contract all spin channels and multipoles, reducing all images together.
2376 : !> \param tb ...
2377 : !> \param matrix_p ...
2378 : !> \param matrix_s ...
2379 : !> \param para_env ...
2380 : !> \param charges ...
2381 : !> \param dipoles ...
2382 : !> \param quadrupoles ...
2383 : ! **************************************************************************************************
2384 53808 : SUBROUTINE tb_contract_population(tb, matrix_p, matrix_s, para_env, charges, dipoles, quadrupoles)
2385 : TYPE(tblite_type), POINTER :: tb
2386 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
2387 : TYPE(mp_para_env_type), POINTER :: para_env
2388 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: charges
2389 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
2390 : INTENT(INOUT) :: dipoles, quadrupoles
2391 :
2392 : INTEGER :: i, icol, im, img, irow, ispin, j, natom, &
2393 : ndip, nquad, nsgf, nspin
2394 : LOGICAL :: found
2395 : REAL(KIND=dp) :: pcol, prow, ps
2396 53808 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: local_population, population
2397 53808 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock, sblock
2398 1076160 : TYPE(block_p_type), DIMENSION(dip_n+quad_n) :: bra, ket
2399 : TYPE(dbcsr_iterator_type) :: iter
2400 :
2401 53808 : nsgf = SIZE(charges, 1)
2402 53808 : natom = SIZE(charges, 2)
2403 53808 : nspin = SIZE(charges, 3)
2404 16418 : ndip = MERGE(dip_n, 0, ALLOCATED(dipoles))
2405 53808 : nquad = MERGE(quad_n, 0, ALLOCATED(quadrupoles))
2406 269040 : ALLOCATE (population(nsgf + ndip + nquad, natom, nspin))
2407 53808 : population = 0.0_dp
2408 :
2409 : !$OMP PARALLEL DEFAULT(NONE) &
2410 : !$OMP SHARED(tb, matrix_p, matrix_s, population, nsgf, natom, nspin, ndip, nquad) &
2411 53808 : !$OMP PRIVATE(local_population, iter, img, irow, icol, sblock, pblock, found, bra, ket, ispin, i, j, im, ps, prow, pcol)
2412 : ALLOCATE (local_population(nsgf + ndip + nquad, natom, nspin))
2413 : local_population = 0.0_dp
2414 : DO img = 1, SIZE(matrix_s, 2)
2415 : CALL dbcsr_iterator_start(iter, matrix_s(1, img)%matrix)
2416 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2417 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
2418 : CALL tb_get_multipole_blocks(tb, irow, icol, img, bra, ket)
2419 : DO ispin = 1, nspin
2420 : CALL dbcsr_get_block_p(matrix_p(ispin, img)%matrix, irow, icol, BLOCK=pblock, found=found)
2421 : IF (.NOT. found) CYCLE
2422 : DO j = 1, SIZE(pblock, 2)
2423 : DO i = 1, SIZE(pblock, 1)
2424 : ps = pblock(i, j)*sblock(i, j)
2425 : local_population(i, irow, ispin) = local_population(i, irow, ispin) + ps
2426 : IF (irow /= icol) THEN
2427 : local_population(j, icol, ispin) = local_population(j, icol, ispin) + ps
2428 : END IF
2429 : END DO
2430 : END DO
2431 : DO im = 1, ndip + nquad
2432 : IF (irow == icol) THEN
2433 : prow = SUM(pblock*bra(im)%block)
2434 : ELSE
2435 : prow = SUM(pblock*ket(im)%block)
2436 : pcol = SUM(pblock*bra(im)%block)
2437 : local_population(nsgf + im, icol, ispin) = local_population(nsgf + im, icol, ispin) + pcol
2438 : END IF
2439 : local_population(nsgf + im, irow, ispin) = local_population(nsgf + im, irow, ispin) + prow
2440 : END DO
2441 : END DO
2442 : END DO
2443 : CALL dbcsr_iterator_stop(iter)
2444 : END DO
2445 : !$OMP CRITICAL(tblite_population)
2446 : population(:, :, :) = population(:, :, :) + local_population(:, :, :)
2447 : !$OMP END CRITICAL(tblite_population)
2448 : DEALLOCATE (local_population)
2449 : !$OMP END PARALLEL
2450 :
2451 : ! Real-space density images already contain the k-point weights.
2452 53808 : CALL para_env%sum(population)
2453 2418568 : charges(:, :, :) = population(1:nsgf, :, :)
2454 113118 : DO ispin = 1, nspin
2455 184260 : DO im = 1, ndip
2456 740028 : dipoles(:, im, ispin) = population(nsgf + im, :, ispin)
2457 : END DO
2458 363018 : DO im = 1, nquad
2459 1420746 : quadrupoles(:, im, ispin) = population(nsgf + ndip + im, :, ispin)
2460 : END DO
2461 : END DO
2462 53808 : DEALLOCATE (population)
2463 107616 : END SUBROUTINE tb_contract_population
2464 : #endif
2465 :
2466 : ! **************************************************************************************************
2467 : !> \brief save gradient to force
2468 : !> \param qs_env ...
2469 : !> \param tb ...
2470 : !> \param para_env ...
2471 : !> \param ityp ...
2472 : !> \note
2473 : ! **************************************************************************************************
2474 1032 : SUBROUTINE tb_grad2force(qs_env, tb, para_env, ityp)
2475 :
2476 : TYPE(qs_environment_type) :: qs_env
2477 : TYPE(tblite_type) :: tb
2478 : TYPE(mp_para_env_type) :: para_env
2479 : INTEGER :: ityp
2480 :
2481 : CHARACTER(len=*), PARAMETER :: routineN = 'tb_grad2force'
2482 :
2483 : INTEGER :: atoma, handle, iatom, ikind, natom
2484 1032 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
2485 1032 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2486 1032 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2487 1032 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2488 :
2489 1032 : CALL timeset(routineN, handle)
2490 :
2491 1032 : NULLIFY (force, atomic_kind_set)
2492 : CALL get_qs_env(qs_env=qs_env, force=force, particle_set=particle_set, &
2493 1032 : atomic_kind_set=atomic_kind_set)
2494 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
2495 1032 : atom_of_kind=atom_of_kind, kind_of=kind_of)
2496 :
2497 1032 : natom = SIZE(particle_set)
2498 :
2499 1032 : SELECT CASE (ityp)
2500 : CASE DEFAULT
2501 0 : CPABORT("unknown force type")
2502 : CASE (0)
2503 730 : DO iatom = 1, natom
2504 644 : ikind = kind_of(iatom)
2505 644 : atoma = atom_of_kind(iatom)
2506 : force(ikind)%all_potential(:, atoma) = &
2507 2662 : force(ikind)%all_potential(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
2508 : END DO
2509 : CASE (1)
2510 1054 : DO iatom = 1, natom
2511 900 : ikind = kind_of(iatom)
2512 900 : atoma = atom_of_kind(iatom)
2513 : force(ikind)%repulsive(:, atoma) = &
2514 3754 : force(ikind)%repulsive(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
2515 : END DO
2516 : CASE (2)
2517 2108 : DO iatom = 1, natom
2518 1800 : ikind = kind_of(iatom)
2519 1800 : atoma = atom_of_kind(iatom)
2520 : force(ikind)%dispersion(:, atoma) = &
2521 7508 : force(ikind)%dispersion(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
2522 : END DO
2523 : CASE (3)
2524 1128 : DO iatom = 1, natom
2525 952 : ikind = kind_of(iatom)
2526 952 : atoma = atom_of_kind(iatom)
2527 : force(ikind)%rho_elec(:, atoma) = &
2528 3984 : force(ikind)%rho_elec(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
2529 : END DO
2530 : CASE (4)
2531 2108 : DO iatom = 1, natom
2532 1800 : ikind = kind_of(iatom)
2533 1800 : atoma = atom_of_kind(iatom)
2534 : force(ikind)%overlap(:, atoma) = &
2535 7508 : force(ikind)%overlap(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
2536 : END DO
2537 : CASE (5)
2538 1032 : DO iatom = 1, natom
2539 0 : ikind = kind_of(iatom)
2540 0 : atoma = atom_of_kind(iatom)
2541 : force(ikind)%efield(:, atoma) = &
2542 0 : force(ikind)%efield(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
2543 : END DO
2544 : END SELECT
2545 :
2546 1032 : CALL timestop(handle)
2547 :
2548 2064 : END SUBROUTINE tb_grad2force
2549 :
2550 : ! **************************************************************************************************
2551 : !> \brief set gradient to zero
2552 : !> \param qs_env ...
2553 : !> \note
2554 : ! **************************************************************************************************
2555 154 : SUBROUTINE tb_zero_force(qs_env)
2556 :
2557 : TYPE(qs_environment_type) :: qs_env
2558 :
2559 : INTEGER :: iatom, ikind, natom
2560 154 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
2561 154 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2562 154 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2563 154 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2564 :
2565 154 : NULLIFY (force, atomic_kind_set)
2566 : CALL get_qs_env(qs_env=qs_env, force=force, particle_set=particle_set, &
2567 154 : atomic_kind_set=atomic_kind_set)
2568 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
2569 154 : kind_of=kind_of)
2570 :
2571 154 : natom = SIZE(particle_set)
2572 :
2573 1054 : DO iatom = 1, natom
2574 900 : ikind = kind_of(iatom)
2575 17780 : force(ikind)%all_potential = 0.0_dp
2576 17780 : force(ikind)%repulsive = 0.0_dp
2577 17780 : force(ikind)%dispersion = 0.0_dp
2578 17780 : force(ikind)%rho_elec = 0.0_dp
2579 17780 : force(ikind)%overlap = 0.0_dp
2580 17934 : force(ikind)%efield = 0.0_dp
2581 : END DO
2582 :
2583 308 : END SUBROUTINE tb_zero_force
2584 :
2585 : ! **************************************************************************************************
2586 : !> \brief Add SCC-overlap and direct multipole Hamiltonian derivatives.
2587 : !> \param qs_env ...
2588 : !> \param use_rho ...
2589 : !> \param nimg ...
2590 : ! **************************************************************************************************
2591 154 : SUBROUTINE tb_derive_dH_off(qs_env, use_rho, nimg)
2592 :
2593 : TYPE(qs_environment_type), POINTER :: qs_env
2594 : LOGICAL, INTENT(IN) :: use_rho
2595 : INTEGER, INTENT(IN) :: nimg
2596 :
2597 : #if defined(__TBLITE)
2598 : INTEGER :: i, idim, ij, iatom, ic, icol, ikind, img, ispin, &
2599 : jdim, ni, nj, nkind, nel, &
2600 : ityp, jatom, jkind, jrow, jtyp, iset, jset, nseti, nsetj, &
2601 : ia, ib, inda, indb, sampled_axes, sampled_even_axes, &
2602 : sampled_gamma_axes, nspin, ikp_axis, slot, maxsgf, natorb_a, natorb_b
2603 : INTEGER, DIMENSION(3) :: cellind, nkp_cellind, nkp_grid
2604 154 : INTEGER, DIMENSION(:), POINTER :: nsgfa, nsgfb
2605 154 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
2606 154 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2607 : LOGICAL :: found, gamma_centered, gamma_sampled_image_pair, &
2608 : has_multipole_response, sampled_image_pair, &
2609 : use_matrix_scc_stress
2610 : LOGICAL, DIMENSION(3) :: mesh_has_gamma
2611 : REAL(KIND=dp) :: r2, dr, i_a_shift, j_a_shift, i_a_shift_mag, j_a_shift_mag, &
2612 : ishift, jshift, ishift_mag, jshift_mag, pij_charge, &
2613 : pij_magnet, mp_pair_scale, kpoint_coordinate, native_dot_tmp
2614 : REAL(KIND=dp), DIMENSION(3) :: kp_shift
2615 : REAL(KIND=dp), DIMENSION(3) :: rij, dgrad, dhgrad_charge, dhgrad_magnet, &
2616 : mpgrad_charge, mpgrad_magnet
2617 : REAL(KIND=dp), DIMENSION(3, 3) :: hsigma
2618 : REAL(KIND=dp), DIMENSION(dip_n) :: idip, jdip, idip_mag, jdip_mag
2619 : REAL(KIND=dp), DIMENSION(quad_n) :: iquad, jquad, iquad_mag, jquad_mag
2620 154 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: t_ov
2621 154 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: local_grad
2622 154 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: t_dip, t_quad, t_d_ov
2623 154 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: t_i_dip, t_i_quad, t_j_dip, t_j_quad
2624 154 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :, :) :: scc_strain_hint
2625 154 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock, pblock_beta
2626 4158 : TYPE(block_p_type), DIMENSION(3, 3, 2) :: scc_strain_blocks
2627 :
2628 154 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2629 154 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_dh_scc, matrix_p
2630 154 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
2631 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
2632 : TYPE(kpoint_type), POINTER :: kpoints
2633 : TYPE(mp_para_env_type), POINTER :: para_env
2634 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2635 154 : POINTER :: sab_orb
2636 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2637 154 : POINTER :: sab_kp
2638 154 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2639 : TYPE(qs_rho_type), POINTER :: rho
2640 : TYPE(qs_scf_env_type), POINTER :: scf_env
2641 : TYPE(tblite_type), POINTER :: tb
2642 : TYPE(virial_type), POINTER :: virial
2643 :
2644 : !$ INTEGER :: hash, lock_num
2645 : !$ INTEGER(KIND=int_8) :: pair_index
2646 154 : !$ INTEGER(KIND=omp_lock_kind), DIMENSION(:), ALLOCATABLE :: locks
2647 : INTEGER, PARAMETER :: nlock = 501
2648 :
2649 : ! compute mulliken charges required for charge update
2650 154 : NULLIFY (scf_env, rho, tb, sab_orb, sab_kp, para_env, kpoints, matrix_dh_scc, virial)
2651 : CALL get_qs_env(qs_env=qs_env, &
2652 : atomic_kind_set=atomic_kind_set, &
2653 : scf_env=scf_env, &
2654 : rho=rho, &
2655 : tb_tblite=tb, &
2656 : sab_orb=sab_orb, &
2657 : sab_kp=sab_kp, &
2658 : para_env=para_env, &
2659 : qs_kind_set=qs_kind_set, &
2660 154 : kpoints=kpoints)
2661 :
2662 154 : NULLIFY (cell_to_index)
2663 154 : IF (nimg > 1) THEN
2664 92 : IF (.NOT. ASSOCIATED(sab_kp)) CPABORT("Missing tblite k-point neighbor list")
2665 92 : sab_orb => sab_kp
2666 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index, nkp_grid=nkp_grid, &
2667 92 : kp_shift=kp_shift, gamma_centered=gamma_centered)
2668 : ELSE
2669 248 : nkp_grid = 1
2670 62 : kp_shift = 0.0_dp
2671 62 : gamma_centered = .TRUE.
2672 : END IF
2673 154 : mesh_has_gamma = .FALSE.
2674 616 : DO i = 1, 3
2675 1312 : DO ikp_axis = 1, nkp_grid(i)
2676 696 : IF (gamma_centered .AND. MODULO(nkp_grid(i), 2) == 0) THEN
2677 : kpoint_coordinate = REAL(2*ikp_axis - nkp_grid(i), KIND=dp)/ &
2678 0 : REAL(2*nkp_grid(i), KIND=dp) + kp_shift(i)
2679 : ELSE
2680 : kpoint_coordinate = REAL(2*ikp_axis - nkp_grid(i) - 1, KIND=dp)/ &
2681 696 : REAL(2*nkp_grid(i), KIND=dp) + kp_shift(i)
2682 : END IF
2683 1158 : IF (ABS(kpoint_coordinate - ANINT(kpoint_coordinate)) < 1.0e-12_dp) THEN
2684 330 : mesh_has_gamma(i) = .TRUE.
2685 : END IF
2686 : END DO
2687 : END DO
2688 154 : has_multipole_response = ASSOCIATED(tb%dipbra) .OR. ASSOCIATED(tb%quadbra)
2689 :
2690 154 : NULLIFY (matrix_p)
2691 154 : IF (use_rho) THEN
2692 132 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2693 22 : ELSE IF (ASSOCIATED(tb%rho_ao_kp_ref)) THEN
2694 22 : matrix_p => tb%rho_ao_kp_ref
2695 : ELSE
2696 0 : matrix_p => scf_env%p_mix_new
2697 : END IF
2698 154 : nspin = SIZE(matrix_p, 1)
2699 400 : use_matrix_scc_stress = tb%use_virial .AND. nimg > 1 .AND. ALL(tb%mol%periodic)
2700 :
2701 : IF (use_matrix_scc_stress) THEN
2702 50 : CALL dbcsr_allocate_matrix_set(matrix_dh_scc, 9*nspin, nimg)
2703 3994 : DO img = 1, nimg
2704 7994 : DO ispin = 1, nspin
2705 19944 : DO jdim = 1, 3
2706 52000 : DO idim = 1, 3
2707 36000 : i = idim + 3*(jdim - 1) + 9*(ispin - 1)
2708 36000 : ALLOCATE (matrix_dh_scc(i, img)%matrix)
2709 : CALL dbcsr_create(matrix_dh_scc(i, img)%matrix, template=matrix_p(ispin, img)%matrix, &
2710 36000 : name="TBLITE SCC STRAIN DERIVATIVE")
2711 48000 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_dh_scc(i, img)%matrix, sab_orb)
2712 : END DO
2713 : END DO
2714 : END DO
2715 : END DO
2716 : END IF
2717 :
2718 : ! set up basis set lists
2719 154 : nkind = SIZE(atomic_kind_set)
2720 764 : ALLOCATE (basis_set_list(nkind))
2721 154 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
2722 :
2723 154 : nel = msao(tb%calc%bas%maxl)**2
2724 154 : CALL get_qs_kind_set(qs_kind_set, maxsgf=maxsgf)
2725 3754 : tb%grad = 0.0_dp
2726 154 : hsigma = 0.0_dp
2727 154 : !$ ALLOCATE (locks(nlock))
2728 77308 : !$ DO lock_num = 1, nlock
2729 77308 : !$ CALL omp_init_lock(locks(lock_num))
2730 : !$ END DO
2731 :
2732 : !$OMP PARALLEL DEFAULT(NONE) &
2733 : !$OMP SHARED(tb, sab_orb, basis_set_list, nimg, nspin, nel, maxsgf, cell_to_index, nkp_grid, mesh_has_gamma, &
2734 : !$OMP has_multipole_response, matrix_p, matrix_dh_scc, use_matrix_scc_stress, locks) &
2735 : !$OMP PRIVATE(local_grad, t_ov, t_d_ov, t_dip, t_i_dip, t_j_dip, t_quad, t_i_quad, t_j_quad, &
2736 : !$OMP idip, jdip, idip_mag, jdip_mag, iquad, jquad, iquad_mag, jquad_mag, scc_strain_hint, &
2737 : !$OMP slot, i, idim, jdim, ij, iatom, jatom, ikind, jkind, icol, jrow, ityp, jtyp, cellind, rij, r2, dr, &
2738 : !$OMP basis_set_a, basis_set_b, first_sgfa, first_sgfb, nsgfa, nsgfb, nseti, nsetj, ic, nkp_cellind, &
2739 : !$OMP sampled_axes, sampled_even_axes, sampled_gamma_axes, sampled_image_pair, gamma_sampled_image_pair, &
2740 : !$OMP mp_pair_scale, pblock, pblock_beta, found, ispin, scc_strain_blocks, i_a_shift, j_a_shift, &
2741 : !$OMP i_a_shift_mag, j_a_shift_mag, ni, nj, iset, jset, ishift, jshift, ishift_mag, jshift_mag, &
2742 : !$OMP dgrad, inda, indb, ia, ib, pij_charge, pij_magnet, mpgrad_charge, mpgrad_magnet, &
2743 154 : !$OMP dhgrad_charge, dhgrad_magnet, natorb_a, natorb_b, pair_index, hash) REDUCTION(+:hsigma)
2744 : ALLOCATE (t_ov(nel), t_d_ov(3, nel), t_dip(dip_n, nel), t_quad(quad_n, nel))
2745 : ALLOCATE (t_i_dip(3, dip_n, nel), t_j_dip(3, dip_n, nel))
2746 : ALLOCATE (t_i_quad(3, quad_n, nel), t_j_quad(3, quad_n, nel))
2747 : ALLOCATE (local_grad(3, tb%mol%nat))
2748 : local_grad = 0.0_dp
2749 : IF (use_matrix_scc_stress) ALLOCATE (scc_strain_hint(maxsgf, maxsgf, 3, 3, nspin))
2750 : !$OMP DO SCHEDULE(GUIDED)
2751 : DO slot = 1, sab_orb(1)%nl_size
2752 : ikind = sab_orb(1)%nlist_task(slot)%ikind
2753 : jkind = sab_orb(1)%nlist_task(slot)%jkind
2754 : iatom = sab_orb(1)%nlist_task(slot)%iatom
2755 : jatom = sab_orb(1)%nlist_task(slot)%jatom
2756 : cellind(:) = sab_orb(1)%nlist_task(slot)%cell(:)
2757 : rij(:) = sab_orb(1)%nlist_task(slot)%r(:)
2758 :
2759 : icol = MAX(iatom, jatom)
2760 : jrow = MIN(iatom, jatom)
2761 :
2762 : IF (iatom < jatom) THEN
2763 : rij = -rij
2764 : i = ikind
2765 : ikind = jkind
2766 : jkind = i
2767 : END IF
2768 :
2769 : ityp = tb%mol%id(icol)
2770 : jtyp = tb%mol%id(jrow)
2771 :
2772 : r2 = DOT_PRODUCT(rij, rij)
2773 : dr = SQRT(r2)
2774 : IF (icol == jrow .AND. dr < same_atom) CYCLE
2775 :
2776 : !get basis information
2777 : basis_set_a => basis_set_list(ikind)%gto_basis_set
2778 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
2779 : first_sgfa => basis_set_a%first_sgf
2780 : nsgfa => basis_set_a%nsgf_set
2781 : nseti = basis_set_a%nset
2782 : basis_set_b => basis_set_list(jkind)%gto_basis_set
2783 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
2784 : first_sgfb => basis_set_b%first_sgf
2785 : nsgfb => basis_set_b%nsgf_set
2786 : nsetj = basis_set_b%nset
2787 :
2788 : IF (nimg == 1) THEN
2789 : ic = 1
2790 : ELSE
2791 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2792 : CPASSERT(ic > 0)
2793 : END IF
2794 : nkp_cellind = 0
2795 : DO i = 1, 3
2796 : IF (nkp_grid(i) > 1) THEN
2797 : nkp_cellind(i) = MODULO(cellind(i), nkp_grid(i))
2798 : IF (2*nkp_cellind(i) > nkp_grid(i)) nkp_cellind(i) = nkp_cellind(i) - nkp_grid(i)
2799 : END IF
2800 : END DO
2801 : sampled_axes = 0
2802 : IF (nkp_cellind(1) /= 0) sampled_axes = sampled_axes + 1
2803 : IF (nkp_cellind(2) /= 0) sampled_axes = sampled_axes + 1
2804 : IF (nkp_cellind(3) /= 0) sampled_axes = sampled_axes + 1
2805 : sampled_even_axes = 0
2806 : sampled_gamma_axes = 0
2807 : DO i = 1, 3
2808 : IF (nkp_grid(i) > 1 .AND. MODULO(nkp_grid(i), 2) == 0 .AND. &
2809 : ABS(2*nkp_cellind(i)) == nkp_grid(i)) THEN
2810 : sampled_even_axes = sampled_even_axes + 1
2811 : IF (mesh_has_gamma(i)) sampled_gamma_axes = sampled_gamma_axes + 1
2812 : END IF
2813 : END DO
2814 : sampled_image_pair = sampled_axes > 0 .AND. sampled_even_axes > 0
2815 : gamma_sampled_image_pair = sampled_image_pair .AND. sampled_gamma_axes == sampled_even_axes
2816 : mp_pair_scale = 1.0_dp
2817 : IF (icol == jrow .AND. sampled_image_pair .AND. has_multipole_response .AND. &
2818 : .NOT. gamma_sampled_image_pair) mp_pair_scale = -1.0_dp
2819 :
2820 : NULLIFY (pblock, pblock_beta)
2821 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
2822 : row=jrow, col=icol, BLOCK=pblock, found=found)
2823 : IF (.NOT. found) CPABORT("pblock not found")
2824 : IF (nspin > 1) THEN
2825 : CALL dbcsr_get_block_p(matrix=matrix_p(2, ic)%matrix, &
2826 : row=jrow, col=icol, BLOCK=pblock_beta, found=found)
2827 : IF (.NOT. found) CPABORT("pblock beta not found")
2828 : END IF
2829 : IF (use_matrix_scc_stress) THEN
2830 : natorb_a = SIZE(pblock, 2)
2831 : natorb_b = SIZE(pblock, 1)
2832 : scc_strain_hint(1:natorb_a, 1:natorb_b, :, :, :) = 0.0_dp
2833 : DO ispin = 1, nspin
2834 : DO jdim = 1, 3
2835 : DO idim = 1, 3
2836 : NULLIFY (scc_strain_blocks(idim, jdim, ispin)%block)
2837 : i = idim + 3*(jdim - 1) + 9*(ispin - 1)
2838 : CALL dbcsr_get_block_p(matrix=matrix_dh_scc(i, ic)%matrix, &
2839 : row=jrow, col=icol, &
2840 : BLOCK=scc_strain_blocks(idim, jdim, ispin)%block, found=found)
2841 : IF (.NOT. found) CPABORT("SCC strain derivative block not found")
2842 : END DO
2843 : END DO
2844 : END DO
2845 : END IF
2846 : i_a_shift = tb%pot%vat(icol, 1)
2847 : j_a_shift = tb%pot%vat(jrow, 1)
2848 : i_a_shift_mag = 0.0_dp
2849 : j_a_shift_mag = 0.0_dp
2850 : IF (SIZE(tb%pot%vat, 2) > 1) THEN
2851 : i_a_shift_mag = tb%pot%vat(icol, 2)
2852 : j_a_shift_mag = tb%pot%vat(jrow, 2)
2853 : END IF
2854 : idip(:) = tb%pot%vdp(:, icol, 1)
2855 : jdip(:) = tb%pot%vdp(:, jrow, 1)
2856 : idip_mag(:) = 0.0_dp
2857 : jdip_mag(:) = 0.0_dp
2858 : IF (SIZE(tb%pot%vdp, 3) > 1) THEN
2859 : idip_mag(:) = tb%pot%vdp(:, icol, 2)
2860 : jdip_mag(:) = tb%pot%vdp(:, jrow, 2)
2861 : END IF
2862 : iquad(:) = tb%pot%vqp(:, icol, 1)
2863 : jquad(:) = tb%pot%vqp(:, jrow, 1)
2864 : iquad_mag(:) = 0.0_dp
2865 : jquad_mag(:) = 0.0_dp
2866 : IF (SIZE(tb%pot%vqp, 3) > 1) THEN
2867 : iquad_mag(:) = tb%pot%vqp(:, icol, 2)
2868 : jquad_mag(:) = tb%pot%vqp(:, jrow, 2)
2869 : END IF
2870 : ni = tb%calc%bas%ish_at(icol)
2871 : DO iset = 1, nseti
2872 : ishift = i_a_shift + tb%pot%vsh(ni + iset, 1)
2873 : ishift_mag = 0.0_dp
2874 : IF (SIZE(tb%pot%vsh, 2) > 1) ishift_mag = i_a_shift_mag + tb%pot%vsh(ni + iset, 2)
2875 : nj = tb%calc%bas%ish_at(jrow)
2876 : DO jset = 1, nsetj
2877 : jshift = j_a_shift + tb%pot%vsh(nj + jset, 1)
2878 : jshift_mag = 0.0_dp
2879 : IF (SIZE(tb%pot%vsh, 2) > 1) jshift_mag = j_a_shift_mag + tb%pot%vsh(nj + jset, 2)
2880 :
2881 : !get integrals and derivatives
2882 : CALL multipole_grad_cgto(tb%calc%bas%cgto(iset, ityp), tb%calc%bas%cgto(jset, jtyp), &
2883 : & r2, rij, tb%calc%bas%intcut, t_ov, t_dip, t_quad, t_d_ov, t_i_dip, t_i_quad, &
2884 : & t_j_dip, t_j_quad)
2885 :
2886 : dgrad = 0.0_dp
2887 : DO inda = 1, nsgfa(iset)
2888 : ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
2889 : DO indb = 1, nsgfb(jset)
2890 : ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
2891 :
2892 : ij = inda + nsgfa(iset)*(indb - 1)
2893 :
2894 : pij_charge = pblock(ib, ia)
2895 : pij_magnet = 0.0_dp
2896 : IF (nspin > 1) THEN
2897 : pij_charge = pij_charge + pblock_beta(ib, ia)
2898 : pij_magnet = pblock(ib, ia) - pblock_beta(ib, ia)
2899 : END IF
2900 : mpgrad_charge = MATMUL(t_i_dip(:, :, ij), idip) &
2901 : + MATMUL(t_j_dip(:, :, ij), jdip) &
2902 : + MATMUL(t_i_quad(:, :, ij), iquad) &
2903 : + MATMUL(t_j_quad(:, :, ij), jquad)
2904 : mpgrad_magnet = MATMUL(t_i_dip(:, :, ij), idip_mag) &
2905 : + MATMUL(t_j_dip(:, :, ij), jdip_mag) &
2906 : + MATMUL(t_i_quad(:, :, ij), iquad_mag) &
2907 : + MATMUL(t_j_quad(:, :, ij), jquad_mag)
2908 : dhgrad_charge = -(ishift + jshift)*t_d_ov(:, ij) - mp_pair_scale*mpgrad_charge
2909 : dhgrad_magnet = -(ishift_mag + jshift_mag)*t_d_ov(:, ij) - mp_pair_scale*mpgrad_magnet
2910 : dgrad(:) = dgrad(:) - &
2911 : ((ishift + jshift)*pij_charge + &
2912 : (ishift_mag + jshift_mag)*pij_magnet)*t_d_ov(:, ij) - &
2913 : mp_pair_scale*(pij_charge*mpgrad_charge + pij_magnet*mpgrad_magnet)
2914 :
2915 : IF (use_matrix_scc_stress) THEN
2916 : DO jdim = 1, 3
2917 : DO idim = 1, 3
2918 : scc_strain_hint(ia, ib, idim, jdim, 1) = &
2919 : scc_strain_hint(ia, ib, idim, jdim, 1) &
2920 : + (dhgrad_charge(idim) + MERGE(dhgrad_magnet(idim), 0.0_dp, nspin > 1))*rij(jdim)
2921 : IF (nspin > 1) THEN
2922 : scc_strain_hint(ia, ib, idim, jdim, 2) = &
2923 : scc_strain_hint(ia, ib, idim, jdim, 2) &
2924 : + (dhgrad_charge(idim) - dhgrad_magnet(idim))*rij(jdim)
2925 : END IF
2926 : END DO
2927 : END DO
2928 : END IF
2929 :
2930 : END DO
2931 : END DO
2932 : local_grad(:, icol) = local_grad(:, icol) - dgrad
2933 : local_grad(:, jrow) = local_grad(:, jrow) + dgrad
2934 : IF (tb%use_virial .AND. .NOT. use_matrix_scc_stress) THEN
2935 : IF (icol == jrow) THEN
2936 : DO ia = 1, 3
2937 : DO ib = 1, 3
2938 : IF (sampled_image_pair .AND. .NOT. gamma_sampled_image_pair) THEN
2939 : hsigma(ia, ib) = hsigma(ia, ib) - 0.25_dp* &
2940 : (rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
2941 : ELSE
2942 : hsigma(ia, ib) = hsigma(ia, ib) + 0.25_dp* &
2943 : (rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
2944 : END IF
2945 : END DO
2946 : END DO
2947 : ELSE
2948 : DO ia = 1, 3
2949 : DO ib = 1, 3
2950 : hsigma(ia, ib) = hsigma(ia, ib) + 0.50_dp*(rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
2951 : END DO
2952 : END DO
2953 : END IF
2954 : END IF
2955 : END DO
2956 : END DO
2957 : IF (use_matrix_scc_stress) THEN
2958 : !$ pair_index = INT(jrow - 1, int_8)*INT(tb%mol%nat, int_8) + INT(icol, int_8)
2959 : !$ hash = INT(MOD(pair_index, INT(nlock, int_8))) + 1
2960 : !$ CALL omp_set_lock(locks(hash))
2961 : DO ispin = 1, nspin
2962 : DO jdim = 1, 3
2963 : DO idim = 1, 3
2964 : IF (icol <= jrow) THEN
2965 : scc_strain_blocks(idim, jdim, ispin)%block(:, :) = &
2966 : scc_strain_blocks(idim, jdim, ispin)%block(:, :) &
2967 : + scc_strain_hint(1:natorb_a, 1:natorb_b, idim, jdim, ispin)
2968 : ELSE
2969 : scc_strain_blocks(idim, jdim, ispin)%block(:, :) = &
2970 : scc_strain_blocks(idim, jdim, ispin)%block(:, :) &
2971 : + TRANSPOSE(scc_strain_hint(1:natorb_a, 1:natorb_b, idim, jdim, ispin))
2972 : END IF
2973 : END DO
2974 : END DO
2975 : END DO
2976 : !$ CALL omp_unset_lock(locks(hash))
2977 : END IF
2978 : END DO
2979 : !$OMP END DO
2980 : !$OMP CRITICAL(tblite_scc_gradient)
2981 : tb%grad(:, :) = tb%grad(:, :) + local_grad(:, :)
2982 : !$OMP END CRITICAL(tblite_scc_gradient)
2983 : DEALLOCATE (local_grad, t_ov, t_d_ov, t_dip, t_quad, t_i_dip, t_j_dip, t_i_quad, t_j_quad)
2984 : IF (ALLOCATED(scc_strain_hint)) DEALLOCATE (scc_strain_hint)
2985 : !$OMP END PARALLEL
2986 77308 : !$ DO lock_num = 1, nlock
2987 77308 : !$ CALL omp_destroy_lock(locks(lock_num))
2988 : !$ END DO
2989 154 : !$ DEALLOCATE (locks)
2990 :
2991 154 : IF (use_matrix_scc_stress) THEN
2992 50 : hsigma = 0.0_dp
2993 3994 : DO img = 1, nimg
2994 7994 : DO ispin = 1, nspin
2995 19944 : DO jdim = 1, 3
2996 52000 : DO idim = 1, 3
2997 36000 : i = idim + 3*(jdim - 1) + 9*(ispin - 1)
2998 36000 : CALL dbcsr_finalize(matrix_dh_scc(i, img)%matrix)
2999 36000 : CALL dbcsr_dot(matrix_dh_scc(i, img)%matrix, matrix_p(ispin, img)%matrix, native_dot_tmp)
3000 : ! Symmetric DBCSR storage represents both orientations of each assembled pair.
3001 48000 : hsigma(idim, jdim) = hsigma(idim, jdim) + 0.5_dp*native_dot_tmp
3002 : END DO
3003 : END DO
3004 : END DO
3005 : END DO
3006 50 : CALL dbcsr_deallocate_matrix_set(matrix_dh_scc)
3007 : ELSE
3008 104 : CALL para_env%sum(hsigma)
3009 : END IF
3010 154 : CALL para_env%sum(tb%grad)
3011 154 : CALL tb_grad2force(qs_env, tb, para_env, 4)
3012 :
3013 1402 : IF (.NOT. use_matrix_scc_stress) tb%sigma = tb%sigma + hsigma
3014 :
3015 154 : DEALLOCATE (basis_set_list)
3016 :
3017 154 : IF (tb%use_virial) THEN
3018 86 : CALL tb_add_stress(qs_env, tb, para_env)
3019 86 : IF (use_matrix_scc_stress) THEN
3020 50 : CALL get_qs_env(qs_env=qs_env, virial=virial)
3021 650 : virial%pv_virial = virial%pv_virial - hsigma/para_env%num_pe
3022 : END IF
3023 : END IF
3024 :
3025 : #else
3026 : MARK_USED(qs_env)
3027 : MARK_USED(use_rho)
3028 : MARK_USED(nimg)
3029 : CPABORT("Built without TBLITE")
3030 : #endif
3031 :
3032 462 : END SUBROUTINE tb_derive_dH_off
3033 :
3034 : ! **************************************************************************************************
3035 : !> \brief save stress tensor
3036 : !> \param qs_env ...
3037 : !> \param tb ...
3038 : !> \param para_env ...
3039 : ! **************************************************************************************************
3040 86 : SUBROUTINE tb_add_stress(qs_env, tb, para_env)
3041 :
3042 : TYPE(qs_environment_type) :: qs_env
3043 : TYPE(tblite_type) :: tb
3044 : TYPE(mp_para_env_type) :: para_env
3045 :
3046 : INTEGER, DIMENSION(3) :: periodic
3047 : TYPE(cell_type), POINTER :: cell
3048 : TYPE(virial_type), POINTER :: virial
3049 :
3050 86 : NULLIFY (virial, cell)
3051 86 : CALL get_qs_env(qs_env=qs_env, virial=virial, cell=cell)
3052 86 : CALL get_cell(cell=cell, periodic=periodic)
3053 :
3054 108 : IF (ALL(periodic == 0)) THEN
3055 : CALL cp_warn(__LOCATION__, &
3056 : "tblite stress tensor requested for an isolated system. "// &
3057 : "The reported virial is useful for finite-difference checks, "// &
3058 4 : "but it is not a physically meaningful bulk stress for an isolated molecule.")
3059 : END IF
3060 :
3061 1118 : virial%pv_virial = virial%pv_virial - tb%sigma/para_env%num_pe
3062 :
3063 86 : END SUBROUTINE tb_add_stress
3064 :
3065 : ! **************************************************************************************************
3066 : !> \brief add contrib. to gradient
3067 : !> \param grad ...
3068 : !> \param deriv ...
3069 : !> \param dE ...
3070 : !> \param natom ...
3071 : ! **************************************************************************************************
3072 154 : SUBROUTINE tb_add_grad(grad, deriv, dE, natom)
3073 :
3074 : REAL(KIND=dp), DIMENSION(:, :) :: grad
3075 : REAL(KIND=dp), DIMENSION(:, :, :) :: deriv
3076 : REAL(KIND=dp), DIMENSION(:) :: dE
3077 : INTEGER :: natom
3078 :
3079 : INTEGER :: i, j
3080 :
3081 1054 : DO i = 1, natom
3082 9990 : DO j = 1, natom
3083 36644 : grad(:, i) = grad(:, i) + deriv(:, i, j)*dE(j)
3084 : END DO
3085 : END DO
3086 :
3087 154 : END SUBROUTINE tb_add_grad
3088 :
3089 : END MODULE tblite_interface
|
