(git:6ba6522)
Loading...
Searching...
No Matches
tblite_interface.F
Go to the documentation of this file.
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
14MODULE tblite_interface
15
16#if defined(__TBLITE)
17 USE mctc_env, ONLY: error_type
18 USE mctc_io, ONLY: structure_type, new
19 USE mctc_io_symbols, ONLY: symbol_to_number
20 USE tblite_adjlist, ONLY: adjacency_list, new_adjacency_list
21 USE tblite_basis_type, ONLY: get_cutoff
22 USE tblite_container, ONLY: container_cache
23 USE tblite_container_type, ONLY: container_type
24 USE tblite_cutoff, ONLY: get_lattice_points
25 USE tblite_data_spin, ONLY: get_spin_constant
26 USE tblite_integral_multipole, ONLY: multipole_cgto, multipole_grad_cgto, maxl, msao
27 USE tblite_integral_type, ONLY: integral_type, new_integral
28 USE tblite_scf, ONLY: get_mixer_dimension
29 USE tblite_scf_info, ONLY: scf_info, atom_resolved, shell_resolved, &
30 orbital_resolved, not_used
31 USE tblite_scf_potential, ONLY: potential_type, new_potential, add_pot_to_h1
32 USE tblite_spin, ONLY: spin_polarization, new_spin_polarization
33 USE tblite_wavefunction_type, ONLY: wavefunction_type, new_wavefunction
34 USE tblite_xtb_calculator, ONLY: xtb_calculator, new_xtb_calculator
35 USE tblite_xtb_gfn1, ONLY: new_gfn1_calculator
36 USE tblite_xtb_gfn2, ONLY: new_gfn2_calculator
37 USE tblite_xtb_h0, ONLY: get_selfenergy, get_hamiltonian, get_occupation, &
38 get_hamiltonian_gradient, tb_hamiltonian
39 USE tblite_xtb_ipea1, ONLY: new_ipea1_calculator
40#endif
41 USE ai_contraction, ONLY: block_add, &
42 contraction
43 USE ai_overlap, ONLY: overlap_ab
44 USE atomic_kind_types, ONLY: atomic_kind_type, get_atomic_kind, get_atomic_kind_set
45 USE atprop_types, ONLY: atprop_type
46 USE basis_set_types, ONLY: gto_basis_set_type, gto_basis_set_p_type, &
47 & allocate_gto_basis_set, write_gto_basis_set, process_gto_basis
48 USE block_p_types, ONLY: block_p_type
49 USE cell_types, ONLY: cell_type, get_cell
50 USE cp_blacs_env, ONLY: cp_blacs_env_type
51 USE cp_control_types, ONLY: dft_control_type, xtb_reference_cli_type
52 USE cp_dbcsr_api, ONLY: dbcsr_type, dbcsr_p_type, dbcsr_create, dbcsr_add, dbcsr_copy, &
53 dbcsr_get_block_p, dbcsr_finalize, &
54 dbcsr_iterator_type, dbcsr_iterator_blocks_left, &
55 dbcsr_iterator_start, dbcsr_iterator_stop, &
56 dbcsr_iterator_next_block
57 USE cp_dbcsr_contrib, ONLY: dbcsr_dot, dbcsr_print
58 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
59 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set, dbcsr_deallocate_matrix_set
60 USE cp_log_handling, ONLY: cp_get_default_logger, &
61 cp_logger_type, cp_logger_get_default_io_unit
62 USE cp_output_handling, ONLY: cp_print_key_should_output, &
63 cp_print_key_unit_nr, cp_print_key_finished_output, &
64 debug_print_level, high_print_level, silent_print_level
65 USE cp_units, ONLY: cp_unit_from_cp2k
66 USE input_constants, ONLY: gfn1xtb, gfn2xtb, ipea1xtb, smear_energy_window, &
67 smear_fermi_dirac, smear_gaussian, smear_list, smear_mp, smear_mv, &
68 tblite_cli_born_kernel_auto, tblite_cli_born_kernel_p16, &
69 tblite_cli_born_kernel_still, tblite_cli_solution_state_bar1mol, &
70 tblite_cli_solution_state_gsolv, tblite_cli_solution_state_reference, &
71 tblite_cli_solvation_alpb, tblite_cli_solvation_cpcm, &
72 tblite_cli_solvation_gb, tblite_cli_solvation_gbe, &
73 tblite_cli_solvation_gbsa, &
74 tblite_guess_ceh, tblite_guess_eeq, tblite_guess_sad, &
75 tblite_mixer_damping_default, &
76 tblite_mixer_max_weight_default, tblite_mixer_min_weight_default, &
77 tblite_mixer_omega0_default, tblite_mixer_weight_factor_default, &
78 tblite_scc_mixer_auto, tblite_scc_mixer_cp2k, &
79 tblite_scc_mixer_none, tblite_scc_mixer_tblite, &
80 tblite_solver_gvd, tblite_solver_gvr
81 USE input_section_types, ONLY: section_vals_val_get
82 USE kinds, ONLY: default_path_length, default_string_length, dp, int_8
83 USE kpoint_types, ONLY: get_kpoint_info, kpoint_type
84 USE memory_utilities, ONLY: reallocate
85 USE message_passing, ONLY: mp_para_env_type
86 USE mulliken, ONLY: ao_charges
87 USE orbital_pointers, ONLY: ncoset
88 USE particle_types, ONLY: particle_type
89 USE qs_charge_mixing, ONLY: charge_mixing, tblite_scc_error_on_cp2k_scale, tblite_scc_pconv
90 USE qs_condnum, ONLY: overlap_condnum
91 USE qs_density_mixing_types, ONLY: modified_broyden_mixing_nr
92 USE qs_energy_types, ONLY: qs_energy_type
93 USE qs_environment_types, ONLY: get_qs_env, qs_environment_type
94 USE qs_force_types, ONLY: qs_force_type, total_qs_force
95 USE qs_integral_utils, ONLY: basis_set_list_setup, get_memory_usage
96 USE qs_kind_types, ONLY: get_qs_kind, qs_kind_type, get_qs_kind_set
97 USE qs_ks_types, ONLY: get_ks_env, qs_ks_env_type, set_ks_env
98 USE qs_neighbor_list_types, ONLY: neighbor_list_iterator_create, neighbor_list_iterate, &
99 get_iterator_info, neighbor_list_set_p_type, &
100 neighbor_list_iterator_p_type, neighbor_list_iterator_release
101 USE qs_overlap, ONLY: build_overlap_matrix, create_sab_matrix
102 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
103 USE qs_rho_types, ONLY: qs_rho_get, qs_rho_type
104 USE qs_scf_types, ONLY: qs_scf_env_type
105 USE scf_control_types, ONLY: scf_control_type
106 USE input_section_types, ONLY: section_vals_get_subs_vals, section_vals_type
107 USE string_utilities, ONLY: integer_to_string
108#if defined(__TBLITE)
109 USE tblite_scc_mixer, ONLY: new_cp2k_tblite_mixer
110#endif
111 USE tblite_types, ONLY: tblite_type, deallocate_tblite_type, allocate_tblite_type
112 USE virial_types, ONLY: virial_type
113 USE xtb_types, ONLY: get_xtb_atom_param, xtb_atom_type
114 USE xtb_types, ONLY: xtb_atom_type
115
116!$ USE OMP_LIB, ONLY: omp_destroy_lock, omp_init_lock, omp_set_lock, &
117!$ omp_unset_lock, omp_lock_kind
118
119#include "./base/base_uses.f90"
120 IMPLICIT NONE
121
122 PRIVATE
123
124 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tblite_interface'
125
126 INTEGER, PARAMETER :: dip_n = 3
127 INTEGER, PARAMETER :: quad_n = 6
128 REAL(KIND=dp), PARAMETER :: same_atom = 0.00001_dp
129
130 PUBLIC :: tb_set_calculator, tb_init_geometry, tb_init_wf
131 PUBLIC :: tb_get_basis, build_tblite_matrices
132 PUBLIC :: tb_get_energy, tb_update_charges, tb_ham_add_coulomb
133 PUBLIC :: tb_native_scc_mixer_active
134 PUBLIC :: tb_scf_mixer_error
135 PUBLIC :: tb_get_multipole
136 PUBLIC :: tb_derive_dh_off
137 PUBLIC :: tb_reference_cli_compare
138
139CONTAINS
140
141#if defined(__TBLITE)
142! **************************************************************************************************
143!> \brief Project a tblite charge/magnetization quantity to a CP2K spin channel.
144!> \param values ...
145!> \param ispin ...
146!> \return ...
147! **************************************************************************************************
148 PURE FUNCTION tb_spin_project(values, ispin) RESULT(value)
149
150 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: values
151 INTEGER, INTENT(IN) :: ispin
152 REAL(KIND=dp) :: value
153
154 value = values(1)
155 IF (SIZE(values) > 1) THEN
156 SELECT CASE (ispin)
157 CASE (1)
158 value = values(1) + values(2)
159 CASE (2)
160 value = values(1) - values(2)
161 CASE DEFAULT
162 value = values(1)
163 END SELECT
164 END IF
165
166 END FUNCTION tb_spin_project
167
168! **************************************************************************************************
169!> \brief Store a private copy of the converged density for tblite force derivatives.
170!> \param tb ...
171!> \param matrix_p ...
172! **************************************************************************************************
173 SUBROUTINE tb_store_density_ref(tb, matrix_p)
174
175 TYPE(tblite_type), POINTER :: tb
176 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
177
178 INTEGER :: img, ispin, nimg, nspin
179
180 nspin = SIZE(matrix_p, 1)
181 nimg = SIZE(matrix_p, 2)
182 IF (ASSOCIATED(tb%rho_ao_kp_ref)) THEN
183 IF (SIZE(tb%rho_ao_kp_ref, 1) /= nspin .OR. SIZE(tb%rho_ao_kp_ref, 2) /= nimg) THEN
184 CALL dbcsr_deallocate_matrix_set(tb%rho_ao_kp_ref)
185 END IF
186 END IF
187 IF (.NOT. ASSOCIATED(tb%rho_ao_kp_ref)) THEN
188 CALL dbcsr_allocate_matrix_set(tb%rho_ao_kp_ref, nspin, nimg)
189 DO img = 1, nimg
190 DO ispin = 1, nspin
191 ALLOCATE (tb%rho_ao_kp_ref(ispin, img)%matrix)
192 END DO
193 END DO
194 END IF
195 DO img = 1, nimg
196 DO ispin = 1, nspin
197 CALL dbcsr_copy(tb%rho_ao_kp_ref(ispin, img)%matrix, matrix_p(ispin, img)%matrix)
198 END DO
199 END DO
200
201 END SUBROUTINE tb_store_density_ref
202#endif
203
204! **************************************************************************************************
205!> \brief intialize geometry objects ...
206!> \param qs_env ...
207!> \param tb ...
208! **************************************************************************************************
209 SUBROUTINE tb_init_geometry(qs_env, tb)
210
211 TYPE(qs_environment_type), POINTER :: qs_env
212 TYPE(tblite_type), POINTER :: tb
213
214#if defined(__TBLITE)
215
216 CHARACTER(LEN=*), PARAMETER :: routinen = 'tblite_init_geometry'
217
218 TYPE(cell_type), POINTER :: cell
219 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
220 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
221 INTEGER :: iatom, natom
222 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: xyz
223 INTEGER :: handle, ikind
224 INTEGER, DIMENSION(3) :: periodic
225 LOGICAL, DIMENSION(3) :: lperiod
226
227 CALL timeset(routinen, handle)
228
229 !get info from environment vaiarable
230 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell, qs_kind_set=qs_kind_set)
231
232 !get information about particles
233 natom = SIZE(particle_set)
234 ALLOCATE (xyz(3, natom))
235 CALL allocate_tblite_type(tb)
236 ALLOCATE (tb%el_num(natom))
237 tb%el_num = -9
238 DO iatom = 1, natom
239 xyz(:, iatom) = particle_set(iatom)%r(:)
240 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
241 CALL get_qs_kind(qs_kind_set(ikind), zatom=tb%el_num(iatom))
242 IF (tb%el_num(iatom) < 1 .OR. tb%el_num(iatom) > 85) THEN
243 cpabort("only elements 1-85 are supported by tblite")
244 END IF
245 END DO
246
247 !get information about cell / lattice
248 CALL get_cell(cell=cell, periodic=periodic)
249 lperiod(1) = periodic(1) == 1
250 lperiod(2) = periodic(2) == 1
251 lperiod(3) = periodic(3) == 1
252
253 !prepare for the call to the dispersion function
254 CALL new(tb%mol, tb%el_num, xyz, lattice=cell%hmat, periodic=lperiod)
255
256 DEALLOCATE (xyz)
257
258 CALL timestop(handle)
259
260#else
261 mark_used(qs_env)
262 mark_used(tb)
263 cpabort("Built without TBLITE")
264#endif
265
266 END SUBROUTINE tb_init_geometry
267
268! **************************************************************************************************
269!> \brief updating coordinates...
270!> \param qs_env ...
271!> \param tb ...
272! **************************************************************************************************
273 SUBROUTINE tb_update_geometry(qs_env, tb)
274
275 TYPE(qs_environment_type) :: qs_env
276 TYPE(tblite_type) :: tb
277
278#if defined(__TBLITE)
279
280 CHARACTER(LEN=*), PARAMETER :: routinen = 'tblite_update_geometry'
281
282 TYPE(cell_type), POINTER :: cell
283 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
284 INTEGER :: iatom, natom
285 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: xyz
286 INTEGER :: handle
287
288 CALL timeset(routinen, handle)
289
290 !get info from environment vaiarable
291 NULLIFY (cell)
292 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell)
293
294 !get information about particles
295 natom = SIZE(particle_set)
296 ALLOCATE (xyz(3, natom))
297 DO iatom = 1, natom
298 xyz(:, iatom) = particle_set(iatom)%r(:)
299 END DO
300 tb%mol%xyz(:, :) = xyz
301 tb%mol%lattice(:, :) = cell%hmat
302
303 DEALLOCATE (xyz)
304
305 CALL timestop(handle)
306
307#else
308 mark_used(qs_env)
309 mark_used(tb)
310 cpabort("Built without TBLITE")
311#endif
312
313 END SUBROUTINE tb_update_geometry
314
315! **************************************************************************************************
316!> \brief initialize wavefunction ...
317!> \param tb ...
318!> \param dft_control ...
319! **************************************************************************************************
320 SUBROUTINE tb_init_wf(tb, dft_control)
321
322 TYPE(tblite_type), POINTER :: tb
323 TYPE(dft_control_type), POINTER :: dft_control
324
325#if defined(__TBLITE)
326
327 INTEGER :: nspin
328
329 TYPE(scf_info) :: info
330
331 nspin = dft_control%nspins
332 IF (nspin /= 1 .AND. nspin /= 2) cpabort("tblite supports only one or two spin channels")
333
334 tb%mol%charge = dft_control%charge
335 tb%mol%uhf = max(0, dft_control%multiplicity - 1)
336 IF (nspin == 2) CALL tb_add_spin_polarization(tb)
337
338 info = tb%calc%variable_info()
339 IF (info%charge > shell_resolved) cpabort("tblite: no support for orbital resolved charge")
340 IF (info%dipole > atom_resolved) cpabort("tblite: no support for shell resolved dipole moment")
341 IF (info%quadrupole > atom_resolved) THEN
342 cpabort("tblite: no support shell resolved quadrupole moment")
343 END IF
344
345 CALL new_wavefunction(tb%wfn, tb%mol%nat, tb%calc%bas%nsh, tb%calc%bas%nao, nspin, 0.0_dp)
346 CALL get_occupation(tb%mol, tb%calc%bas, tb%calc%h0, tb%wfn%nocc, tb%wfn%n0at, tb%wfn%n0sh)
347 CALL tb_reset_mixer(tb)
348
349 CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
350
351 !allocate quantities later required
352 ALLOCATE (tb%e_hal(tb%mol%nat), tb%e_rep(tb%mol%nat), tb%e_disp(tb%mol%nat))
353 ALLOCATE (tb%e_scd(tb%mol%nat), tb%e_es(tb%mol%nat), tb%e_int(tb%mol%nat))
354 ALLOCATE (tb%selfenergy(tb%calc%bas%nsh))
355 IF (ALLOCATED(tb%calc%ncoord)) ALLOCATE (tb%cn(tb%mol%nat))
356
357#else
358 mark_used(tb)
359 mark_used(dft_control)
360 cpabort("Built without TBLITE")
361#endif
362
363 END SUBROUTINE tb_init_wf
364
365#if defined(__TBLITE)
366! **************************************************************************************************
367!> \brief Add tblite's on-site spin-polarization interaction.
368!> \param tb ...
369! **************************************************************************************************
370 SUBROUTINE tb_add_spin_polarization(tb)
371
372 TYPE(tblite_type), POINTER :: tb
373
374 CLASS(container_type), ALLOCATABLE :: cont
375 TYPE(spin_polarization), ALLOCATABLE :: spin
376 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wll
377
378 ALLOCATE (spin)
379 CALL tb_get_spin_constants(tb, wll)
380 CALL new_spin_polarization(spin, tb%mol, wll, tb%calc%bas%nsh_id)
381 CALL move_alloc(spin, cont)
382 CALL tb%calc%push_back(cont)
383
384 END SUBROUTINE tb_add_spin_polarization
385
386! **************************************************************************************************
387!> \brief Build tblite spin constants for the current basis.
388!> \param tb ...
389!> \param wll ...
390! **************************************************************************************************
391 SUBROUTINE tb_get_spin_constants(tb, wll)
392
393 TYPE(tblite_type), POINTER :: tb
394 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
395 INTENT(OUT) :: wll
396
397 INTEGER :: il, ish, izp, jl, jsh
398
399 ALLOCATE (wll(tb%calc%bas%nsh, tb%calc%bas%nsh, tb%mol%nid))
400 wll = 0.0_dp
401 DO izp = 1, tb%mol%nid
402 DO ish = 1, tb%calc%bas%nsh_id(izp)
403 il = tb%calc%bas%cgto(ish, izp)%ang
404 DO jsh = 1, tb%calc%bas%nsh_id(izp)
405 jl = tb%calc%bas%cgto(jsh, izp)%ang
406 wll(jsh, ish, izp) = get_spin_constant(jl, il, tb%mol%num(izp))
407 END DO
408 END DO
409 END DO
410
411 END SUBROUTINE tb_get_spin_constants
412#endif
413
414! **************************************************************************************************
415!> \brief Reset tblite's internal SCC mixer for a new CP2K SCF cycle.
416!> \param tb ...
417! **************************************************************************************************
418 SUBROUTINE tb_reset_mixer(tb)
419
420 TYPE(tblite_type), POINTER :: tb
421
422#if defined(__TBLITE)
423
424 TYPE(scf_info) :: info
425
426 info = tb%calc%variable_info()
427 IF (ALLOCATED(tb%mixer)) DEALLOCATE (tb%mixer)
428 CALL new_cp2k_tblite_mixer(tb%mixer, tb%mixer_memory, &
429 tb%wfn%nspin*get_mixer_dimension(tb%mol, tb%calc%bas, info), &
430 tb%mixer_damping, tb%mixer_omega0, tb%mixer_min_weight, &
431 tb%mixer_max_weight, tb%mixer_weight_factor)
432
433#else
434 mark_used(tb)
435 cpabort("Built without TBLITE")
436#endif
437
438 END SUBROUTINE tb_reset_mixer
439
440! **************************************************************************************************
441!> \brief Configure tblite's internal SCC mixer from CP2K input.
442!> \param tb ...
443!> \param iterations ...
444!> \param memory ...
445!> \param damping ...
446!> \param omega0 ...
447!> \param min_weight ...
448!> \param max_weight ...
449!> \param weight_factor ...
450!> \param solver ...
451! **************************************************************************************************
452 SUBROUTINE tb_configure_mixer(tb, iterations, memory, damping, omega0, min_weight, max_weight, &
453 weight_factor, solver)
454
455 TYPE(tblite_type), POINTER :: tb
456 INTEGER, INTENT(IN) :: iterations, memory, solver
457 REAL(kind=dp), INTENT(IN) :: damping, max_weight, min_weight, omega0, &
458 weight_factor
459
460#if defined(__TBLITE)
461
462 IF (iterations < 1) cpabort("tblite SCC mixer ITERATIONS must be positive")
463 IF (memory < 1) cpabort("tblite SCC mixer MEMORY must be positive")
464 IF (damping <= 0.0_dp) cpabort("tblite SCC mixer damping must be positive")
465 IF (omega0 <= 0.0_dp) cpabort("tblite SCC mixer OMEGA0 must be positive")
466 IF (min_weight <= 0.0_dp) cpabort("tblite SCC mixer MIN_WEIGHT must be positive")
467 IF (max_weight <= 0.0_dp) cpabort("tblite SCC mixer MAX_WEIGHT must be positive")
468 IF (max_weight < min_weight) THEN
469 cpabort("tblite SCC mixer MAX_WEIGHT must not be smaller than MIN_WEIGHT")
470 END IF
471 IF (weight_factor <= 0.0_dp) cpabort("tblite SCC mixer WEIGHT_FACTOR must be positive")
472 SELECT CASE (solver)
473 CASE (tblite_solver_gvd, tblite_solver_gvr)
474 CASE DEFAULT
475 cpabort("Unknown tblite SCC mixer SOLVER")
476 END SELECT
477
478 tb%calc%max_iter = iterations
479 tb%mixer_memory = memory
480 tb%mixer_solver = solver
481 tb%mixer_damping = damping
482 tb%calc%mixer_input%damping = damping
483 tb%mixer_omega0 = omega0
484 tb%mixer_min_weight = min_weight
485 tb%mixer_max_weight = max_weight
486 tb%mixer_weight_factor = weight_factor
487
488#else
489 mark_used(tb)
490 mark_used(iterations)
491 mark_used(memory)
492 mark_used(damping)
493 mark_used(omega0)
494 mark_used(min_weight)
495 mark_used(max_weight)
496 mark_used(weight_factor)
497 mark_used(solver)
498 cpabort("Built without TBLITE")
499#endif
500
501 END SUBROUTINE tb_configure_mixer
502
503! **************************************************************************************************
504!> \brief Return whether the tblite native SCC mixer is active for this run.
505!> \param dft_control ...
506!> \return ...
507! **************************************************************************************************
508 FUNCTION tb_native_scc_mixer_active(dft_control) RESULT(use_native_mixer)
509
510 TYPE(dft_control_type), POINTER :: dft_control
511 LOGICAL :: use_native_mixer
512
513 use_native_mixer = .false.
514 IF (.NOT. ASSOCIATED(dft_control)) RETURN
515 IF (dft_control%qs_control%do_ls_scf) RETURN
516
517 SELECT CASE (dft_control%qs_control%xtb_control%tblite_scc_mixer)
518 CASE (tblite_scc_mixer_auto)
519 use_native_mixer = .true.
520 CASE (tblite_scc_mixer_tblite)
521 use_native_mixer = .true.
522 CASE (tblite_scc_mixer_cp2k, tblite_scc_mixer_none)
523 use_native_mixer = .false.
524 CASE DEFAULT
525 cpabort("Unknown tblite SCC mixer")
526 END SELECT
527
528 END FUNCTION tb_native_scc_mixer_active
529
530! **************************************************************************************************
531!> \brief Return the native tblite SCC mixer residual on the CP2K iter_delta scale.
532!> \param dft_control ...
533!> \param tb ...
534!> \param eps_scf CP2K reporting scale for the native residual.
535!> \return ...
536! **************************************************************************************************
537 FUNCTION tb_scf_mixer_error(dft_control, tb, eps_scf) RESULT(mixer_error)
538
539 TYPE(dft_control_type), POINTER :: dft_control
540 TYPE(tblite_type), POINTER :: tb
541 REAL(kind=dp), INTENT(IN) :: eps_scf
542 REAL(kind=dp) :: mixer_error
543
544#if defined(__TBLITE)
545 REAL(kind=dp) :: raw_error
546#endif
547
548 mixer_error = 0.0_dp
549
550#if defined(__TBLITE)
551 IF (.NOT. ASSOCIATED(tb)) RETURN
552 IF (.NOT. tb_native_scc_mixer_active(dft_control)) RETURN
553 IF (ALLOCATED(tb%mixer)) THEN
554 raw_error = real(tb%mixer%get_error(), kind=dp)
555 mixer_error = tblite_scc_error_on_cp2k_scale( &
556 raw_error, eps_scf, &
557 tblite_scc_pconv*dft_control%qs_control%xtb_control%tblite_accuracy)
558 END IF
559#else
560 mark_used(dft_control)
561 mark_used(tb)
562 mark_used(eps_scf)
563#endif
564
565 END FUNCTION tb_scf_mixer_error
566
567! **************************************************************************************************
568!> \brief ...
569!> \param tb ...
570!> \param typ ...
571!> \param accuracy ...
572!> \param param_file ...
573! **************************************************************************************************
574 SUBROUTINE tb_set_calculator(tb, typ, accuracy, param_file)
575
576 TYPE(tblite_type), POINTER :: tb
577 INTEGER :: typ
578 REAL(kind=dp), INTENT(IN) :: accuracy
579 CHARACTER(LEN=*), INTENT(IN) :: param_file
580
581#if defined(__TBLITE)
582
583 TYPE(error_type), ALLOCATABLE :: error
584
585 IF (ALLOCATED(tb%param)) DEALLOCATE (tb%param)
586 IF (len_trim(param_file) > 0) THEN
587 ALLOCATE (tb%param)
588 CALL tb%param%load(trim(param_file), error)
589 IF (ALLOCATED(error)) cpabort("Could not load tblite PARAM file: "//trim(param_file))
590 CALL new_xtb_calculator(tb%calc, tb%mol, tb%param, error)
591 ELSE
592 SELECT CASE (typ)
593 CASE default
594 cpabort("Unknown xtb type")
595 CASE (gfn1xtb)
596 CALL new_gfn1_calculator(tb%calc, tb%mol, error)
597 CASE (gfn2xtb)
598 CALL new_gfn2_calculator(tb%calc, tb%mol, error)
599 CASE (ipea1xtb)
600 CALL new_ipea1_calculator(tb%calc, tb%mol, error)
601 END SELECT
602 END IF
603 IF (ALLOCATED(error)) cpabort("tblite calculator setup failed")
604
605 tb%accuracy = accuracy
606
607#else
608 mark_used(tb)
609 mark_used(typ)
610 mark_used(accuracy)
611 mark_used(param_file)
612 cpabort("Built without TBLITE")
613#endif
614
615 END SUBROUTINE tb_set_calculator
616
617! **************************************************************************************************
618!> \brief ...
619!> \param qs_env ...
620!> \param tb ...
621!> \param para_env ...
622! **************************************************************************************************
623 SUBROUTINE tb_init_ham(qs_env, tb, para_env)
624
625 TYPE(qs_environment_type) :: qs_env
626 TYPE(tblite_type) :: tb
627 TYPE(mp_para_env_type) :: para_env
628
629#if defined(__TBLITE)
630
631 TYPE(container_cache) :: hcache, rcache
632
633 tb%e_hal = 0.0_dp
634 tb%e_rep = 0.0_dp
635 tb%e_disp = 0.0_dp
636 tb%e_int = 0.0_dp
637 IF (ALLOCATED(tb%grad)) THEN
638 tb%grad = 0.0_dp
639 CALL tb_zero_force(qs_env)
640 END IF
641 tb%sigma = 0.0_dp
642
643 IF (ALLOCATED(tb%calc%halogen)) THEN
644 CALL tb%calc%halogen%update(tb%mol, hcache)
645 IF (ALLOCATED(tb%grad)) THEN
646 tb%grad = 0.0_dp
647 CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal, &
648 & tb%grad, tb%sigma)
649 CALL tb_dump_sigma_component("after_halogen", tb%sigma, para_env)
650 CALL tb_grad2force(qs_env, tb, para_env, 0)
651 ELSE
652 CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal)
653 END IF
654 END IF
655
656 IF (ALLOCATED(tb%calc%repulsion)) THEN
657 CALL tb%calc%repulsion%update(tb%mol, rcache)
658 IF (ALLOCATED(tb%grad)) THEN
659 tb%grad = 0.0_dp
660 CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep, &
661 & tb%grad, tb%sigma)
662 CALL tb_dump_sigma_component("after_repulsion", tb%sigma, para_env)
663 CALL tb_grad2force(qs_env, tb, para_env, 1)
664 ELSE
665 CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep)
666 END IF
667 END IF
668
669 IF (ALLOCATED(tb%calc%dispersion)) THEN
670 CALL tb%calc%dispersion%update(tb%mol, tb%dcache)
671 IF (ALLOCATED(tb%grad)) THEN
672 tb%grad = 0.0_dp
673 CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp, &
674 & tb%grad, tb%sigma)
675 CALL tb_dump_sigma_component("after_dispersion_static", tb%sigma, para_env)
676 CALL tb_grad2force(qs_env, tb, para_env, 2)
677 ELSE
678 CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp)
679 END IF
680 END IF
681
682 IF (ALLOCATED(tb%calc%interactions)) THEN
683 CALL tb%calc%interactions%update(tb%mol, tb%icache)
684 END IF
685
686 CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
687 IF (ALLOCATED(tb%calc%coulomb)) THEN
688 CALL tb%calc%coulomb%update(tb%mol, tb%cache)
689 END IF
690
691 IF (ALLOCATED(tb%grad)) THEN
692 IF (ALLOCATED(tb%calc%ncoord)) THEN
693 CALL tb%calc%ncoord%get_cn(tb%mol, tb%cn, tb%dcndr, tb%dcndL)
694 END IF
695 CALL get_selfenergy(tb%calc%h0, tb%mol%id, tb%calc%bas%ish_at, &
696 & tb%calc%bas%nsh_id, cn=tb%cn, selfenergy=tb%selfenergy, dsedcn=tb%dsedcn)
697 ELSE
698 IF (ALLOCATED(tb%calc%ncoord)) THEN
699 CALL tb%calc%ncoord%get_cn(tb%mol, tb%cn)
700 END IF
701 CALL get_selfenergy(tb%calc%h0, tb%mol%id, tb%calc%bas%ish_at, &
702 & tb%calc%bas%nsh_id, cn=tb%cn, selfenergy=tb%selfenergy, dsedcn=tb%dsedcn)
703 END IF
704
705#else
706 mark_used(qs_env)
707 mark_used(tb)
708 mark_used(para_env)
709 cpabort("Built without TBLITE")
710#endif
711
712 END SUBROUTINE tb_init_ham
713
714! **************************************************************************************************
715!> \brief ...
716!> \param qs_env ...
717!> \param tb ...
718!> \param energy ...
719! **************************************************************************************************
720 SUBROUTINE tb_get_energy(qs_env, tb, energy)
721
722 TYPE(qs_environment_type), POINTER :: qs_env
723 TYPE(tblite_type), POINTER :: tb
724 TYPE(qs_energy_type), POINTER :: energy
725
726#if defined(__TBLITE)
727
728 INTEGER :: iounit
729 TYPE(cp_logger_type), POINTER :: logger
730 TYPE(section_vals_type), POINTER :: scf_section
731 REAL(kind=dp) :: xtb_inter
732 NULLIFY (scf_section, logger)
733
734 logger => cp_get_default_logger()
735 iounit = cp_logger_get_default_io_unit(logger)
736 scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
737
738 energy%repulsive = sum(tb%e_rep)
739 energy%el_stat = sum(tb%e_es)
740 energy%dispersion = sum(tb%e_disp)
741 energy%dispersion_sc = sum(tb%e_scd)
742 energy%xtb_xb_inter = sum(tb%e_hal)
743 xtb_inter = sum(tb%e_int)
744
745 energy%total = energy%core + energy%repulsive + energy%el_stat + energy%dispersion &
746 + energy%dispersion_sc + energy%xtb_xb_inter + xtb_inter &
747 + energy%kTS + energy%efield + energy%qmmm_el
748
749 iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
750 extension=".scfLog")
751 IF (iounit > 0) THEN
752 WRITE (unit=iounit, fmt="(/,(T9,A,T60,F20.10))") &
753 "Repulsive pair potential energy: ", energy%repulsive, &
754 "Zeroth order Hamiltonian energy: ", energy%core, &
755 "Electrostatic energy: ", energy%el_stat, &
756 "Self-consistent dispersion energy: ", energy%dispersion_sc, &
757 "Non-self consistent dispersion energy: ", energy%dispersion
758 IF (abs(energy%xtb_xb_inter) > 1.e-9_dp) THEN
759 WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
760 "Correction for halogen bonding: ", energy%xtb_xb_inter
761 END IF
762 IF (abs(xtb_inter) > 1.e-9_dp) THEN
763 WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
764 "Additional interaction (e.g. spin): ", xtb_inter
765 END IF
766 IF (abs(energy%efield) > 1.e-9_dp) THEN
767 WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
768 "Electric field interaction energy: ", energy%efield
769 END IF
770 IF (qs_env%qmmm) THEN
771 WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
772 "QM/MM Electrostatic energy: ", energy%qmmm_el
773 END IF
774 END IF
775 CALL cp_print_key_finished_output(iounit, logger, scf_section, &
776 "PRINT%DETAILED_ENERGY")
777
778#else
779 mark_used(qs_env)
780 mark_used(tb)
781 mark_used(energy)
782 cpabort("Built without TBLITE")
783#endif
784
785 END SUBROUTINE tb_get_energy
786
787! **************************************************************************************************
788!> \brief ...
789!> \param tb ...
790!> \param gto_basis_set ...
791!> \param element_symbol ...
792!> \param param ...
793!> \param occ ...
794! **************************************************************************************************
795 SUBROUTINE tb_get_basis(tb, gto_basis_set, element_symbol, param, occ)
796
797 TYPE(tblite_type), POINTER :: tb
798 TYPE(gto_basis_set_type), POINTER :: gto_basis_set
799 CHARACTER(len=2), INTENT(IN) :: element_symbol
800 TYPE(xtb_atom_type), POINTER :: param
801 INTEGER, DIMENSION(5), INTENT(out) :: occ
802
803#if defined(__TBLITE)
804
805 REAL(kind=dp) :: docc
806 CHARACTER(LEN=default_string_length) :: sng
807 INTEGER :: ang, i_type, id_atom, ind_ao, ipgf, ish, &
808 ishell, ityp, maxl, mprim, natorb, &
809 nset, nshell
810 LOGICAL :: do_ortho
811
812 CALL allocate_gto_basis_set(gto_basis_set)
813
814 !identifying element in the bas data
815 CALL symbol_to_number(i_type, element_symbol)
816 DO id_atom = 1, tb%mol%nat
817 IF (i_type == tb%el_num(id_atom)) EXIT
818 END DO
819 param%z = i_type
820 param%symbol = element_symbol
821 param%defined = .true.
822 ityp = tb%mol%id(id_atom)
823
824 !getting size information
825 nset = tb%calc%bas%nsh_id(ityp)
826 nshell = 1
827 mprim = 0
828 DO ishell = 1, nset
829 mprim = max(mprim, tb%calc%bas%cgto(ishell, ityp)%nprim)
830 END DO
831 param%nshell = nset
832 natorb = 0
833
834 !write basis set information
835 CALL integer_to_string(mprim, sng)
836 gto_basis_set%name = element_symbol//"_STO-"//trim(sng)//"G"
837 gto_basis_set%nset = nset
838 CALL reallocate(gto_basis_set%lmax, 1, nset)
839 CALL reallocate(gto_basis_set%lmin, 1, nset)
840 CALL reallocate(gto_basis_set%npgf, 1, nset)
841 CALL reallocate(gto_basis_set%nshell, 1, nset)
842 CALL reallocate(gto_basis_set%n, 1, 1, 1, nset)
843 CALL reallocate(gto_basis_set%l, 1, 1, 1, nset)
844 CALL reallocate(gto_basis_set%zet, 1, mprim, 1, nset)
845 CALL reallocate(gto_basis_set%gcc, 1, mprim, 1, 1, 1, nset)
846
847 ind_ao = 0
848 maxl = 0
849 DO ishell = 1, nset
850 ang = tb%calc%bas%cgto(ishell, ityp)%ang
851 natorb = natorb + (2*ang + 1)
852 param%lval(ishell) = ang
853 maxl = max(ang, maxl)
854 gto_basis_set%lmax(ishell) = ang
855 gto_basis_set%lmin(ishell) = ang
856 gto_basis_set%npgf(ishell) = tb%calc%bas%cgto(ishell, ityp)%nprim
857 gto_basis_set%nshell(ishell) = nshell
858 gto_basis_set%n(1, ishell) = ang + 1
859 gto_basis_set%l(1, ishell) = ang
860 DO ipgf = 1, gto_basis_set%npgf(ishell)
861 gto_basis_set%gcc(ipgf, 1, ishell) = tb%calc%bas%cgto(ishell, ityp)%coeff(ipgf)
862 gto_basis_set%zet(ipgf, ishell) = tb%calc%bas%cgto(ishell, ityp)%alpha(ipgf)
863 END DO
864 DO ipgf = 1, (2*ang + 1)
865 ind_ao = ind_ao + 1
866 param%lao(ind_ao) = ang
867 param%nao(ind_ao) = ishell
868 END DO
869 END DO
870
871 do_ortho = .false.
872 CALL process_gto_basis(gto_basis_set, do_ortho, nset, maxl)
873
874 !setting additional values in parameter
875 param%rcut = get_cutoff(tb%calc%bas, tb%accuracy)
876 param%natorb = natorb
877 param%lmax = maxl !max angular momentum
878
879 !getting occupation
880 occ = 0
881 docc = 0.0_dp
882 IF (tb%calc%bas%nsh_at(id_atom) > 5) cpabort("too many shells in tblite")
883 DO ish = 1, tb%calc%bas%nsh_at(id_atom)
884 occ(ish) = nint(tb%calc%h0%refocc(ish, ityp) + docc)
885 docc = docc + tb%calc%h0%refocc(ish, ityp) - real(occ(ish))
886 param%occupation(ish) = occ(ish)
887 END DO
888 IF (abs(docc) > 0.1_dp) cpabort("Getting occupation numbers from tblite fails")
889 param%zeff = sum(occ) !effective core charge
890
891 !set normalization process
892 gto_basis_set%norm_type = 3
893
894#else
895 occ = 0
896 mark_used(tb)
897 mark_used(gto_basis_set)
898 mark_used(element_symbol)
899 mark_used(param)
900 cpabort("Built without TBLITE")
901#endif
902
903 END SUBROUTINE tb_get_basis
904
905! **************************************************************************************************
906!> \brief ...
907!> \param qs_env ...
908!> \param calculate_forces ...
909! **************************************************************************************************
910 SUBROUTINE build_tblite_matrices(qs_env, calculate_forces)
911
912 TYPE(qs_environment_type), POINTER :: qs_env
913 LOGICAL, INTENT(IN) :: calculate_forces
914
915#if defined(__TBLITE)
916
917 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_tblite_matrices'
918
919 INTEGER :: handle, maxder, nderivatives, nimg, img, nkind, i, &
920 ic, iw, iatom, jatom, ikind, jkind, iset, jset, n1, n2, icol, &
921 irow, ia, ib, sgfa, sgfb, ldsab, nseta, nsetb, &
922 natorb_a, natorb_b, raw_iatom, raw_jatom, &
923 slot, idim, jdim
924 LOGICAL :: found, norml1, norml2, use_arnoldi
925 REAL(kind=dp) :: dr, dshpoly, ff, hij_base, r2, rr
926 REAL(kind=dp) :: native_dot_tmp, native_cn_icol, native_cn_irow, &
927 native_dot_weight
928 INTEGER, DIMENSION(3) :: cell
929 REAL(kind=dp) :: hij, shpoly
930 REAL(kind=dp), DIMENSION(2) :: condnum
931 REAL(kind=dp), DIMENSION(3) :: native_h0_overlap_force, native_radial_force, raw_rij, rij
932 REAL(kind=dp), DIMENSION(3, 3) :: native_radial_dot
933 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
934 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
935 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
936 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: native_cn_deriv_thread
937 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: native_radial_force_thread
938 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint, hint
939 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: radial_hint
940 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min
941 INTEGER, DIMENSION(:), POINTER :: npgfa, npgfb, nsgfa, nsgfb
942 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
943 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
944 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, zeta, zetb, scon_a, scon_b
945 REAL(kind=dp), DIMENSION(:, :), POINTER :: fblock, pblock, sblock
946 TYPE(block_p_type), DIMENSION(3, 3) :: radial_blocks
947!$ INTEGER :: hash, hash1, lock_num
948!$ INTEGER(KIND=int_8) :: iatom8
949!$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
950 INTEGER, PARAMETER :: nlock = 501
951
952 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
953 TYPE(atprop_type), POINTER :: atprop
954 TYPE(cp_blacs_env_type), POINTER :: blacs_env
955 TYPE(cp_logger_type), POINTER :: logger
956 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_dh_native, matrix_h, matrix_p, &
957 matrix_q_native, matrix_s, matrix_s_native, matrix_w
958 TYPE(dft_control_type), POINTER :: dft_control
959 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
960 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
961 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
962 TYPE(kpoint_type), POINTER :: kpoints
963 TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
964 TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_kp
965 TYPE(mp_para_env_type), POINTER :: para_env
966 TYPE(qs_energy_type), POINTER :: energy
967 TYPE(qs_ks_env_type), POINTER :: ks_env
968 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
969 TYPE(qs_rho_type), POINTER :: rho
970 TYPE(tblite_type), POINTER :: tb
971 TYPE(tb_hamiltonian), POINTER :: h0
972 TYPE(virial_type), POINTER :: virial
973
974 CALL timeset(routinen, handle)
975
976 NULLIFY (ks_env, energy, atomic_kind_set, qs_kind_set)
977 NULLIFY (matrix_dh_native, matrix_h, matrix_q_native, matrix_s, matrix_s_native, atprop, dft_control)
978 NULLIFY (sab_orb, sab_kp, rho, tb, kpoints, cell_to_index)
979
980 CALL get_qs_env(qs_env=qs_env, &
981 ks_env=ks_env, para_env=para_env, &
982 energy=energy, &
983 atomic_kind_set=atomic_kind_set, &
984 qs_kind_set=qs_kind_set, &
985 matrix_h_kp=matrix_h, &
986 matrix_s_kp=matrix_s, &
987 atprop=atprop, &
988 dft_control=dft_control, &
989 sab_orb=sab_orb, &
990 sab_kp=sab_kp, &
991 rho=rho, tb_tblite=tb)
992 h0 => tb%calc%h0
993
994 !update geometry (required for debug / geometry optimization)
995 CALL tb_update_geometry(qs_env, tb)
996
997 nkind = SIZE(atomic_kind_set)
998 nderivatives = 0
999 IF (calculate_forces) THEN
1000 nderivatives = 1
1001 IF (ALLOCATED(tb%grad)) DEALLOCATE (tb%grad)
1002 ALLOCATE (tb%grad(3, tb%mol%nat))
1003 IF (ALLOCATED(tb%dsedcn)) DEALLOCATE (tb%dsedcn)
1004 ALLOCATE (tb%dsedcn(tb%calc%bas%nsh))
1005 IF (ALLOCATED(tb%calc%ncoord)) THEN
1006 IF (ALLOCATED(tb%dcndr)) DEALLOCATE (tb%dcndr)
1007 ALLOCATE (tb%dcndr(3, tb%mol%nat, tb%mol%nat))
1008 IF (ALLOCATED(tb%dcndL)) DEALLOCATE (tb%dcndL)
1009 ALLOCATE (tb%dcndL(3, 3, tb%mol%nat))
1010 END IF
1011 ELSE
1012 IF (ALLOCATED(tb%grad)) DEALLOCATE (tb%grad)
1013 IF (ALLOCATED(tb%dcndr)) DEALLOCATE (tb%dcndr)
1014 IF (ALLOCATED(tb%dcndL)) DEALLOCATE (tb%dcndL)
1015 END IF
1016 maxder = ncoset(nderivatives)
1017 nimg = dft_control%nimages
1018 IF (nimg > 1) THEN
1019 IF (.NOT. ASSOCIATED(sab_kp)) cpabort("Missing k-point neighbor list for tblite")
1020 sab_orb => sab_kp
1021 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1022 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1023 END IF
1024
1025 !intialise hamiltonian
1026 CALL tb_init_ham(qs_env, tb, para_env)
1027
1028 ! get density matrtix
1029 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1030
1031 ! set up matrices for force calculations
1032 IF (calculate_forces) THEN
1033 NULLIFY (force, matrix_w, virial)
1034 CALL get_qs_env(qs_env=qs_env, &
1035 matrix_w_kp=matrix_w, &
1036 virial=virial, force=force)
1037
1038 IF (SIZE(matrix_p, 1) == 2) THEN
1039 DO img = 1, nimg
1040 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
1041 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1042 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
1043 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1044 END DO
1045 END IF
1046 tb%use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1047 END IF
1048
1049 IF (calculate_forces) THEN
1050 CALL dbcsr_allocate_matrix_set(matrix_q_native, 1, nimg)
1051 DO img = 1, nimg
1052 ALLOCATE (matrix_q_native(1, img)%matrix)
1053 CALL dbcsr_copy(matrix_q_native(1, img)%matrix, matrix_w(1, img)%matrix, &
1054 name="TBLITE NATIVE OVERLAP FORCE MATRIX")
1055 END DO
1056 END IF
1057
1058 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
1059 IF (calculate_forces) THEN
1060 ALLOCATE (native_radial_force_thread(3, SIZE(atom_of_kind)))
1061 ALLOCATE (native_cn_deriv_thread(SIZE(atom_of_kind)))
1062 native_cn_deriv_thread = 0.0_dp
1063 native_radial_force_thread = 0.0_dp
1064 END IF
1065
1066 ! set up basis set lists
1067 ALLOCATE (basis_set_list(nkind))
1068 CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
1069
1070 ! allocate overlap matrix
1071 CALL dbcsr_allocate_matrix_set(matrix_s, 1, nimg)
1072 CALL create_sab_matrix(ks_env, matrix_s, "OVERLAP MATRIX", basis_set_list, basis_set_list, &
1073 sab_orb, .true.)
1074 CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
1075
1076 ! initialize H matrix
1077 CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
1078 DO img = 1, nimg
1079 ALLOCATE (matrix_h(1, img)%matrix)
1080 CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, img)%matrix, &
1081 name="HAMILTONIAN MATRIX")
1082 CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
1083 END DO
1084 CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
1085 IF (calculate_forces .AND. nimg > 1) THEN
1086 CALL dbcsr_allocate_matrix_set(matrix_dh_native, 9, nimg)
1087 DO img = 1, nimg
1088 DO idim = 1, 9
1089 ALLOCATE (matrix_dh_native(idim, img)%matrix)
1090 CALL dbcsr_create(matrix_dh_native(idim, img)%matrix, template=matrix_s(1, img)%matrix, &
1091 name="TBLITE H0 STRAIN DERIVATIVE")
1092 CALL cp_dbcsr_alloc_block_from_nbl(matrix_dh_native(idim, img)%matrix, sab_orb)
1093 END DO
1094 END DO
1095 END IF
1096 ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
1097
1098 native_radial_dot = 0.0_dp
1099
1100 ! loop over all atom pairs with a non-zero overlap (sab_orb)
1101!$OMP PARALLEL DEFAULT(NONE) &
1102!$OMP SHARED (calculate_forces, basis_set_list, sab_orb, matrix_dh_native, matrix_s, matrix_h, matrix_p, &
1103!$OMP cell_to_index, nimg, atom_of_kind, qs_kind_set, tb, h0, ldsab, maxder, locks, &
1104!$OMP native_cn_deriv_thread, native_radial_force_thread, ncoset) &
1105!$OMP PRIVATE (slot, hash, hash1, iatom8, lock_num, ikind, jkind, iatom, jatom, cell, rij, ic, irow, &
1106!$OMP icol, dr, dshpoly, ff, hij_base, r2, raw_iatom, raw_jatom, &
1107!$OMP raw_rij, native_cn_icol, native_cn_irow, native_dot_weight, native_h0_overlap_force, &
1108!$OMP native_radial_force, &
1109!$OMP n1, n2, ia, ib, i, idim, jdim, iset, jset, sgfa, sgfb, &
1110!$OMP nseta, nsetb, natorb_a, &
1111!$OMP natorb_b, found, zeta, first_sgfa, la_max, la_min, npgfa, nsgfa, rpgfa, set_radius_a, &
1112!$OMP scon_a, first_sgfb, lb_max, lb_min, npgfb, nsgfb, rpgfb, set_radius_b, scon_b, zetb, rr, hij, shpoly, &
1113!$OMP owork, oint, sint, hint, radial_hint, basis_set_a, basis_set_b, sblock, fblock, &
1114!$OMP pblock, radial_blocks) &
1115!$OMP REDUCTION (+ : native_radial_dot)
1116
1117!$OMP SINGLE
1118!$ ALLOCATE (locks(nlock))
1119!$OMP END SINGLE
1120!$OMP DO
1121!$ DO lock_num = 1, nlock
1122!$ CALL omp_init_lock(locks(lock_num))
1123!$ END DO
1124!$OMP END DO
1125
1126 ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
1127
1128!$OMP DO SCHEDULE(GUIDED)
1129 DO slot = 1, sab_orb(1)%nl_size
1130
1131 ikind = sab_orb(1)%nlist_task(slot)%ikind
1132 jkind = sab_orb(1)%nlist_task(slot)%jkind
1133 iatom = sab_orb(1)%nlist_task(slot)%iatom
1134 jatom = sab_orb(1)%nlist_task(slot)%jatom
1135 cell(:) = sab_orb(1)%nlist_task(slot)%cell(:)
1136 rij(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
1137 raw_iatom = iatom
1138 raw_jatom = jatom
1139 raw_rij = rij
1140 native_cn_icol = 0.0_dp
1141 native_cn_irow = 0.0_dp
1142 native_h0_overlap_force = 0.0_dp
1143 native_radial_force = 0.0_dp
1144
1145 ! canonicalize pair ordering as in current serial flow
1146 icol = max(iatom, jatom)
1147 irow = min(iatom, jatom)
1148 IF (iatom < jatom) THEN
1149 rij = -rij
1150 i = ikind
1151 ikind = jkind
1152 jkind = i
1153 END IF
1154
1155 dr = norm2(rij(:))
1156 r2 = dr*dr
1157
1158 IF (nimg == 1) THEN
1159 ic = 1
1160 ELSE
1161 ic = cell_to_index(cell(1), cell(2), cell(3))
1162 cpassert(ic > 0)
1163 END IF
1164
1165 NULLIFY (sblock)
1166 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
1167 row=irow, col=icol, block=sblock, found=found)
1168 cpassert(found)
1169 NULLIFY (fblock)
1170 CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
1171 row=irow, col=icol, block=fblock, found=found)
1172 cpassert(found)
1173 IF (calculate_forces) THEN
1174 NULLIFY (pblock)
1175 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
1176 row=irow, col=icol, block=pblock, found=found)
1177 cpassert(found)
1178 IF (nimg > 1) THEN
1179 DO jdim = 1, 3
1180 DO idim = 1, 3
1181 NULLIFY (radial_blocks(idim, jdim)%block)
1182 CALL dbcsr_get_block_p(matrix=matrix_dh_native(idim + 3*(jdim - 1), ic)%matrix, &
1183 row=irow, col=icol, block=radial_blocks(idim, jdim)%block, &
1184 found=found)
1185 cpassert(found)
1186 END DO
1187 END DO
1188 END IF
1189 END IF
1190
1191 ! --------- Overlap
1192 !get basis information
1193 basis_set_a => basis_set_list(ikind)%gto_basis_set
1194 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1195 basis_set_b => basis_set_list(jkind)%gto_basis_set
1196 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1197 ! basis a
1198 first_sgfa => basis_set_a%first_sgf
1199 la_max => basis_set_a%lmax
1200 la_min => basis_set_a%lmin
1201 npgfa => basis_set_a%npgf
1202 nseta = basis_set_a%nset
1203 nsgfa => basis_set_a%nsgf_set
1204 rpgfa => basis_set_a%pgf_radius
1205 set_radius_a => basis_set_a%set_radius
1206 scon_a => basis_set_a%scon
1207 zeta => basis_set_a%zet
1208 ! basis b
1209 first_sgfb => basis_set_b%first_sgf
1210 lb_max => basis_set_b%lmax
1211 lb_min => basis_set_b%lmin
1212 npgfb => basis_set_b%npgf
1213 nsetb = basis_set_b%nset
1214 nsgfb => basis_set_b%nsgf_set
1215 rpgfb => basis_set_b%pgf_radius
1216 set_radius_b => basis_set_b%set_radius
1217 scon_b => basis_set_b%scon
1218 zetb => basis_set_b%zet
1219
1220 natorb_a = 0
1221 DO iset = 1, nseta
1222 natorb_a = natorb_a + (2*basis_set_a%l(1, iset) + 1)
1223 END DO
1224 natorb_b = 0
1225 DO iset = 1, nsetb
1226 natorb_b = natorb_b + (2*basis_set_b%l(1, iset) + 1)
1227 END DO
1228 ALLOCATE (sint(natorb_a, natorb_b, maxder))
1229 sint = 0.0_dp
1230 ALLOCATE (hint(natorb_a, natorb_b, maxder))
1231 hint = 0.0_dp
1232 IF (calculate_forces .AND. nimg > 1) THEN
1233 ALLOCATE (radial_hint(natorb_a, natorb_b, 3, 3))
1234 radial_hint = 0.0_dp
1235 END IF
1236
1237 !----------------- overlap integrals
1238 DO iset = 1, nseta
1239 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1240 sgfa = first_sgfa(1, iset)
1241 DO jset = 1, nsetb
1242 IF (set_radius_a(iset) + set_radius_b(jset) < dr) cycle
1243 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1244 sgfb = first_sgfb(1, jset)
1245 IF (calculate_forces) THEN
1246 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1247 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1248 rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
1249 ELSE
1250 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1251 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1252 rij, sab=oint(:, :, 1))
1253 END IF
1254 ! Contraction
1255 CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1256 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
1257 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.false.)
1258 IF (calculate_forces) THEN
1259 DO i = 2, 4
1260 CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1261 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
1262 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.false.)
1263 END DO
1264 END IF
1265 END DO
1266 END DO
1267
1268!$ iatom8 = INT(iatom - 1, int_8)*INT(SIZE(atom_of_kind), int_8) + INT(jatom, int_8)
1269!$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
1270 ! update S matrix
1271!$ hash = hash1
1272!$ CALL omp_set_lock(locks(hash))
1273 IF (icol <= irow) THEN
1274 sblock(:, :) = sblock(:, :) + sint(:, :, 1)
1275 ELSE
1276 sblock(:, :) = sblock(:, :) + transpose(sint(:, :, 1))
1277 END IF
1278!$ CALL omp_unset_lock(locks(hash))
1279
1280 ! --------- Hamiltonian
1281 IF (icol == irow .AND. dr < same_atom) THEN
1282 !get diagonal F matrix from selfenergy
1283 n1 = tb%calc%bas%ish_at(icol)
1284 DO iset = 1, nseta
1285 sgfa = first_sgfa(1, iset)
1286 hij = tb%selfenergy(n1 + iset)
1287 DO ia = sgfa, sgfa + nsgfa(iset) - 1
1288 hint(ia, ia, 1) = hij
1289 IF (calculate_forces) THEN
1290 native_cn_icol = native_cn_icol + tb%dsedcn(n1 + iset)*pblock(ia, ia)
1291 DO jdim = 1, 3
1292 DO idim = 1, 3
1293 IF (nimg == 1) THEN
1294 native_radial_dot(idim, jdim) = native_radial_dot(idim, jdim) + &
1295 tb%dsedcn(n1 + iset)*tb%dcndL(idim, jdim, icol)*pblock(ia, ia)
1296 ELSE
1297 radial_hint(ia, ia, idim, jdim) = radial_hint(ia, ia, idim, jdim) + &
1298 tb%dsedcn(n1 + iset)*tb%dcndL(idim, jdim, icol)
1299 END IF
1300 END DO
1301 END DO
1302 END IF
1303 END DO
1304 END DO
1305 ELSE
1306 !get off-diagonal F matrix
1307 rr = sqrt(dr/(h0%rad(jkind) + h0%rad(ikind)))
1308 n1 = tb%calc%bas%ish_at(icol)
1309 DO iset = 1, nseta
1310 sgfa = first_sgfa(1, iset)
1311 n2 = tb%calc%bas%ish_at(irow)
1312 DO jset = 1, nsetb
1313 sgfb = first_sgfb(1, jset)
1314 shpoly = (1.0_dp + h0%shpoly(iset, ikind)*rr) &
1315 *(1.0_dp + h0%shpoly(jset, jkind)*rr)
1316 dshpoly = ((1.0_dp + h0%shpoly(iset, ikind)*rr)*h0%shpoly(jset, jkind)*rr &
1317 + (1.0_dp + h0%shpoly(jset, jkind)*rr)*h0%shpoly(iset, ikind)*rr) &
1318 *0.5_dp/r2
1319 hij_base = 0.5_dp*(tb%selfenergy(n1 + iset) + tb%selfenergy(n2 + jset)) &
1320 *h0%hscale(iset, jset, ikind, jkind)
1321 hij = hij_base*shpoly
1322 DO ia = sgfa, sgfa + nsgfa(iset) - 1
1323 DO ib = sgfb, sgfb + nsgfb(jset) - 1
1324 hint(ia, ib, 1) = hij*sint(ia, ib, 1)
1325 IF (calculate_forces) THEN
1326 native_dot_weight = 2.0_dp
1327 IF (icol == irow) native_dot_weight = 1.0_dp
1328 native_cn_icol = native_cn_icol + native_dot_weight*0.5_dp* &
1329 h0%hscale(iset, jset, ikind, jkind)*shpoly* &
1330 tb%dsedcn(n1 + iset)*pblock(ib, ia)*sint(ia, ib, 1)
1331 native_cn_irow = native_cn_irow + native_dot_weight*0.5_dp* &
1332 h0%hscale(iset, jset, ikind, jkind)*shpoly* &
1333 tb%dsedcn(n2 + jset)*pblock(ib, ia)*sint(ia, ib, 1)
1334 native_radial_force = native_radial_force + &
1335 hij_base*dshpoly*pblock(ib, ia)*sint(ia, ib, 1)*raw_rij
1336 native_h0_overlap_force = native_h0_overlap_force + &
1337 hij*pblock(ib, ia)*sint(ia, ib, 2:4)
1338 DO jdim = 1, 3
1339 DO idim = 1, 3
1340 IF (nimg == 1) THEN
1341 native_radial_dot(idim, jdim) = native_radial_dot(idim, jdim) + &
1342 native_dot_weight*pblock(ib, ia)*( &
1343 (-hij)*sint(ia, ib, idim + 1)*rij(jdim) + &
1344 hij_base*dshpoly*sint(ia, ib, 1)*rij(idim)*rij(jdim) + &
1345 0.5_dp*h0%hscale(iset, jset, ikind, jkind)*shpoly*sint(ia, ib, 1)* &
1346 (tb%dsedcn(n1 + iset)*tb%dcndL(idim, jdim, icol) + &
1347 tb%dsedcn(n2 + jset)*tb%dcndL(idim, jdim, irow)))
1348 ELSE
1349 radial_hint(ia, ib, idim, jdim) = radial_hint(ia, ib, idim, jdim) + &
1350 (-hij)*sint(ia, ib, idim + 1)*rij(jdim) + &
1351 hij_base*dshpoly*sint(ia, ib, 1)*rij(idim)*rij(jdim) + &
1352 0.5_dp*h0%hscale(iset, jset, ikind, jkind)*shpoly*sint(ia, ib, 1)* &
1353 (tb%dsedcn(n1 + iset)*tb%dcndL(idim, jdim, icol) + &
1354 tb%dsedcn(n2 + jset)*tb%dcndL(idim, jdim, irow))
1355 END IF
1356 END DO
1357 END DO
1358 END IF
1359 END DO
1360 END DO
1361 END DO
1362 END DO
1363 END IF
1364
1365 ! update F matrix
1366!$ CALL omp_set_lock(locks(hash))
1367 IF (icol <= irow) THEN
1368 fblock(:, :) = fblock(:, :) + hint(:, :, 1)
1369 IF (calculate_forces .AND. nimg > 1) THEN
1370 DO jdim = 1, 3
1371 DO idim = 1, 3
1372 radial_blocks(idim, jdim)%block(:, :) = radial_blocks(idim, jdim)%block(:, :) + &
1373 radial_hint(:, :, idim, jdim)
1374 END DO
1375 END DO
1376 END IF
1377 ELSE
1378 fblock(:, :) = fblock(:, :) + transpose(hint(:, :, 1))
1379 IF (calculate_forces .AND. nimg > 1) THEN
1380 DO jdim = 1, 3
1381 DO idim = 1, 3
1382 radial_blocks(idim, jdim)%block(:, :) = radial_blocks(idim, jdim)%block(:, :) + &
1383 transpose(radial_hint(:, :, idim, jdim))
1384 END DO
1385 END DO
1386 END IF
1387 END IF
1388!$ CALL omp_unset_lock(locks(hash))
1389
1390 IF (calculate_forces) THEN
1391 ff = 2.0_dp
1392!$OMP CRITICAL(tblite_native_radial_force)
1393 native_radial_force_thread(:, raw_iatom) = &
1394 native_radial_force_thread(:, raw_iatom) - ff*native_radial_force
1395 native_radial_force_thread(:, raw_jatom) = &
1396 native_radial_force_thread(:, raw_jatom) + ff*native_radial_force
1397 native_radial_force_thread(:, icol) = &
1398 native_radial_force_thread(:, icol) + ff*native_h0_overlap_force
1399 native_radial_force_thread(:, irow) = &
1400 native_radial_force_thread(:, irow) - ff*native_h0_overlap_force
1401 native_cn_deriv_thread(icol) = native_cn_deriv_thread(icol) + native_cn_icol
1402 native_cn_deriv_thread(irow) = native_cn_deriv_thread(irow) + native_cn_irow
1403!$OMP END CRITICAL(tblite_native_radial_force)
1404 END IF
1405
1406 DEALLOCATE (sint, hint)
1407 IF (ALLOCATED(radial_hint)) DEALLOCATE (radial_hint)
1408
1409 END DO
1410!$OMP END DO
1411 DEALLOCATE (oint, owork)
1412!$OMP DO
1413!$ DO lock_num = 1, nlock
1414!$ CALL omp_destroy_lock(locks(lock_num))
1415!$ END DO
1416!$OMP END DO
1417!$OMP SINGLE
1418!$ DEALLOCATE (locks)
1419!$OMP END SINGLE NOWAIT
1420!$OMP END PARALLEL
1421
1422 IF (calculate_forces) THEN
1423 IF (nimg > 1) THEN
1424 native_radial_dot = 0.0_dp
1425 DO img = 1, nimg
1426 DO jdim = 1, 3
1427 DO idim = 1, 3
1428 CALL dbcsr_finalize(matrix_dh_native(idim + 3*(jdim - 1), img)%matrix)
1429 CALL dbcsr_dot(matrix_dh_native(idim + 3*(jdim - 1), img)%matrix, &
1430 matrix_p(1, img)%matrix, native_dot_tmp)
1431 native_radial_dot(idim, jdim) = native_radial_dot(idim, jdim) + native_dot_tmp
1432 END DO
1433 END DO
1434 END DO
1435 CALL dbcsr_deallocate_matrix_set(matrix_dh_native)
1436 ELSE
1437 CALL para_env%sum(native_radial_dot)
1438 END IF
1439 CALL para_env%sum(native_cn_deriv_thread)
1440 tb%grad = 0.0_dp
1441 CALL tb_add_grad(tb%grad, tb%dcndr, native_cn_deriv_thread, tb%mol%nat)
1442 CALL tb_grad2force(qs_env, tb, para_env, 4)
1443 DO iatom = 1, SIZE(atom_of_kind)
1444 ikind = kind_of(iatom)
1445 force(ikind)%overlap(:, atom_of_kind(iatom)) = &
1446 force(ikind)%overlap(:, atom_of_kind(iatom)) + native_radial_force_thread(:, iatom)
1447 END DO
1448 IF (tb%use_virial) THEN
1449 virial%pv_overlap = virial%pv_overlap - native_radial_dot/para_env%num_pe
1450 virial%pv_virial = virial%pv_virial - native_radial_dot/para_env%num_pe
1451 END IF
1452 END IF
1453
1454 DO img = 1, nimg
1455 DO i = 1, SIZE(matrix_s, 1)
1456 CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1457 END DO
1458 DO i = 1, SIZE(matrix_h, 1)
1459 CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1460 END DO
1461 END DO
1462
1463 IF (calculate_forces) THEN
1464 CALL build_overlap_matrix(ks_env, matrixkp_s=matrix_s_native, &
1465 matrix_name="OVERLAP MATRIX", nderivative=nderivatives, &
1466 basis_type_a="ORB", basis_type_b="ORB", sab_nl=sab_orb, &
1467 calculate_forces=.true., matrixkp_p=matrix_q_native)
1468 CALL dbcsr_deallocate_matrix_set(matrix_s_native)
1469 CALL dbcsr_deallocate_matrix_set(matrix_q_native)
1470 END IF
1471
1472 !compute multipole moments for gfn2
1473 IF (dft_control%qs_control%xtb_control%tblite_method == gfn2xtb) THEN
1474 CALL tb_get_multipole(qs_env, tb)
1475 END IF
1476
1477 ! output overlap information
1478 NULLIFY (logger)
1479 logger => cp_get_default_logger()
1480 IF (.NOT. calculate_forces) THEN
1481 IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
1482 "DFT%PRINT%OVERLAP_CONDITION") /= 0) THEN
1483 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
1484 extension=".Log")
1485 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
1486 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
1487 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
1488 CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
1489 CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
1490 END IF
1491 END IF
1492
1493 DEALLOCATE (basis_set_list)
1494 IF (ALLOCATED(native_cn_deriv_thread)) DEALLOCATE (native_cn_deriv_thread)
1495 IF (ALLOCATED(native_radial_force_thread)) DEALLOCATE (native_radial_force_thread)
1496
1497 CALL timestop(handle)
1498
1499#else
1500 mark_used(qs_env)
1501 mark_used(calculate_forces)
1502 cpabort("Built without TBLITE")
1503#endif
1504
1505 END SUBROUTINE build_tblite_matrices
1506
1507! **************************************************************************************************
1508!> \brief ...
1509!> \param qs_env ...
1510!> \param dft_control ...
1511!> \param tb ...
1512!> \param calculate_forces ...
1513!> \param use_rho ...
1514! **************************************************************************************************
1515 SUBROUTINE tb_update_charges(qs_env, dft_control, tb, calculate_forces, use_rho)
1516
1517 TYPE(qs_environment_type), POINTER :: qs_env
1518 TYPE(dft_control_type), POINTER :: dft_control
1519 TYPE(tblite_type), POINTER :: tb
1520 LOGICAL, INTENT(IN) :: calculate_forces
1521 LOGICAL, INTENT(IN) :: use_rho
1522
1523#if defined(__TBLITE)
1524
1525 INTEGER :: iatom, ikind, is, ns, atom_a, ii, im
1526 INTEGER :: ispin, nspin
1527 INTEGER :: nimg, nkind, nsgf, natorb, na, n_mix_cols, mix_offset
1528 INTEGER :: n_atom, max_orb, max_shell
1529 INTEGER :: raw_state_status, raw_state_unit
1530 LOGICAL :: advance_native_mixer, discard_mixed_output, do_combined_mixing, &
1531 do_dipole, do_quadrupole, native_sign_mixing, &
1532 skip_charge_mixing, reuse_native_input, skip_scf_dispersion, &
1533 seed_native_from_rho, use_native_mixer, use_no_mixer
1534 REAL(kind=dp) :: native_seed_charge, norm, new_charge, pao
1535#if defined(__TBLITE_DEBUG_DIAGNOSTICS)
1536 INTEGER :: debug_status
1537 CHARACTER(LEN=32) :: debug_value
1538#endif
1539 CHARACTER(LEN=default_path_length) :: raw_state_file
1540 INTEGER, DIMENSION(5) :: occ
1541 INTEGER, DIMENSION(25) :: lao
1542 INTEGER, DIMENSION(25) :: nao
1543 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ch_atom, ch_shell, ch_ref, ch_orb, mix_vars
1544 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, ao_dip, ao_quad
1545 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: aocg_spin, ao_dip_spin, ao_quad_spin, &
1546 ch_orb_spin, ch_shell_spin
1547
1548 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1549 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_p
1550 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
1551 TYPE(dbcsr_type), POINTER :: s_matrix
1552 TYPE(error_type), ALLOCATABLE :: error
1553 TYPE(mp_para_env_type), POINTER :: para_env
1554 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1555 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1556 TYPE(qs_rho_type), POINTER :: rho
1557 TYPE(qs_scf_env_type), POINTER :: scf_env
1558 TYPE(scf_control_type), POINTER :: scf_control
1559 TYPE(xtb_atom_type), POINTER :: xtb_kind
1560
1561 ! compute mulliken charges required for charge update
1562 NULLIFY (particle_set, qs_kind_set, atomic_kind_set, scf_control, scf_env)
1563 CALL get_qs_env(qs_env=qs_env, scf_env=scf_env, particle_set=particle_set, qs_kind_set=qs_kind_set, &
1564 atomic_kind_set=atomic_kind_set, matrix_s_kp=matrix_s, rho=rho, para_env=para_env, &
1565 scf_control=scf_control)
1566
1567 ! also compute multipoles needed by GFN2
1568 do_dipole = .false.
1569 do_quadrupole = .false.
1570 skip_scf_dispersion = .false.
1571#if defined(__TBLITE_DEBUG_DIAGNOSTICS)
1572 CALL get_environment_variable("CP2K_TBLITE_DEBUG_SKIP_SCF_DISPERSION", debug_value, status=debug_status)
1573 IF (debug_status == 0) THEN
1574 READ (debug_value, *, iostat=debug_status) skip_scf_dispersion
1575 IF (debug_status /= 0) skip_scf_dispersion = .false.
1576 END IF
1577#endif
1578 IF (dft_control%qs_control%do_ls_scf .OR. scf_control%use_ot) THEN
1579 use_native_mixer = .false.
1580 use_no_mixer = .true.
1581 ELSE
1582 SELECT CASE (dft_control%qs_control%xtb_control%tblite_scc_mixer)
1583 CASE (tblite_scc_mixer_auto)
1584 use_native_mixer = tb_native_scc_mixer_active(dft_control)
1585 use_no_mixer = .false.
1586 CASE (tblite_scc_mixer_tblite)
1587 use_native_mixer = .true.
1588 use_no_mixer = .false.
1589 CASE (tblite_scc_mixer_cp2k)
1590 use_native_mixer = .false.
1591 use_no_mixer = .false.
1592 CASE (tblite_scc_mixer_none)
1593 use_native_mixer = .false.
1594 use_no_mixer = .true.
1595 CASE DEFAULT
1596 cpabort("Unknown tblite SCC mixer")
1597 END SELECT
1598 END IF
1599 IF (use_native_mixer) THEN
1600 IF (.NOT. ASSOCIATED(scf_env)) cpabort("tblite SCC mixer requires a QS SCF environment")
1601 IF (use_rho .AND. (.NOT. calculate_forces)) THEN
1602 IF (scf_env%iter_count > dft_control%qs_control%xtb_control%tblite_mixer_iterations) THEN
1603 cpabort("tblite SCC mixer exceeded TBLITE_MIXER/ITERATIONS")
1604 END IF
1605 IF (scf_env%iter_count == 1) THEN
1606 CALL tb_configure_mixer(tb, dft_control%qs_control%xtb_control%tblite_mixer_iterations, &
1607 dft_control%qs_control%xtb_control%tblite_mixer_memory, &
1608 dft_control%qs_control%xtb_control%tblite_mixer_damping, &
1609 dft_control%qs_control%xtb_control%tblite_mixer_omega0, &
1610 dft_control%qs_control%xtb_control%tblite_mixer_min_weight, &
1611 dft_control%qs_control%xtb_control%tblite_mixer_max_weight, &
1612 dft_control%qs_control%xtb_control%tblite_mixer_weight_factor, &
1613 dft_control%qs_control%xtb_control%tblite_mixer_solver)
1614 CALL tb_reset_mixer(tb)
1615 END IF
1616 END IF
1617 END IF
1618 nspin = dft_control%nspins
1619 IF (nspin /= tb%wfn%nspin) cpabort("CP2K/tblite spin channel mismatch")
1620
1621 NULLIFY (matrix_p)
1622 IF (use_rho) THEN
1623 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1624 ELSE IF (calculate_forces .AND. nspin > 1) THEN
1625 IF (.NOT. ASSOCIATED(tb%rho_ao_kp_ref)) THEN
1626 cpabort("Missing converged tblite density for UKS/LSD forces")
1627 END IF
1628 matrix_p => tb%rho_ao_kp_ref
1629 ELSE
1630 matrix_p => scf_env%p_mix_new
1631 END IF
1632 IF (nspin > 1 .AND. (.NOT. calculate_forces)) CALL tb_store_density_ref(tb, matrix_p)
1633 IF (ASSOCIATED(tb%dipbra)) do_dipole = .true.
1634 IF (ASSOCIATED(tb%quadbra)) do_quadrupole = .true.
1635 reuse_native_input = .false.
1636 IF (use_native_mixer) THEN
1637 IF (scf_env%iter_count == 1) THEN
1638 reuse_native_input = any(abs(tb%wfn%qsh) > 1.0e-14_dp)
1639 IF (do_dipole) reuse_native_input = reuse_native_input .OR. &
1640 any(abs(tb%wfn%dpat) > 1.0e-14_dp)
1641 IF (do_quadrupole) reuse_native_input = reuse_native_input .OR. &
1642 any(abs(tb%wfn%qpat) > 1.0e-14_dp)
1643 END IF
1644 END IF
1645 n_atom = SIZE(particle_set)
1646 nkind = SIZE(atomic_kind_set)
1647 nimg = dft_control%nimages
1648 CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
1649 ALLOCATE (aocg(nsgf, n_atom))
1650 ALLOCATE (aocg_spin(nsgf, n_atom, nspin))
1651 aocg = 0.0_dp
1652 aocg_spin = 0.0_dp
1653 IF (do_dipole) THEN
1654 ALLOCATE (ao_dip(n_atom, dip_n))
1655 ALLOCATE (ao_dip_spin(n_atom, dip_n, nspin))
1656 ao_dip_spin = 0.0_dp
1657 END IF
1658 IF (do_quadrupole) THEN
1659 ALLOCATE (ao_quad(n_atom, quad_n))
1660 ALLOCATE (ao_quad_spin(n_atom, quad_n, nspin))
1661 ao_quad_spin = 0.0_dp
1662 END IF
1663 max_orb = 0
1664 max_shell = 0
1665 DO ikind = 1, nkind
1666 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1667 CALL get_xtb_atom_param(xtb_kind, natorb=natorb)
1668 max_orb = max(max_orb, natorb)
1669 END DO
1670 DO is = 1, n_atom
1671 max_shell = max(max_shell, tb%calc%bas%nsh_at(is))
1672 END DO
1673 ALLOCATE (ch_atom(n_atom, nspin), ch_shell(n_atom, max_shell))
1674 ALLOCATE (ch_orb(max_orb, n_atom), ch_ref(max_orb, n_atom))
1675 ALLOCATE (ch_orb_spin(max_orb, n_atom, nspin), ch_shell_spin(n_atom, max_shell, nspin))
1676 ch_atom = 0.0_dp
1677 ch_shell = 0.0_dp
1678 ch_orb = 0.0_dp
1679 ch_orb_spin = 0.0_dp
1680 ch_shell_spin = 0.0_dp
1681 ch_ref = 0.0_dp
1682 IF (nimg > 1) THEN
1683 DO ispin = 1, nspin
1684 CALL tb_ao_charges_kp_spin(matrix_p, matrix_s, aocg_spin(:, :, ispin), ispin, para_env)
1685 IF (do_dipole) THEN
1686 DO im = 1, dip_n
1687 CALL tb_contract_dens_kp_spin(matrix_p, tb%dipbra, tb%dipket, im, dip_n, &
1688 ao_dip_spin(:, im, ispin), ispin, para_env)
1689 END DO
1690 END IF
1691 IF (do_quadrupole) THEN
1692 DO im = 1, quad_n
1693 CALL tb_contract_dens_kp_spin(matrix_p, tb%quadbra, tb%quadket, im, quad_n, &
1694 ao_quad_spin(:, im, ispin), ispin, para_env)
1695 END DO
1696 END IF
1697 END DO
1698 ELSE
1699 NULLIFY (p_matrix, s_matrix)
1700 p_matrix => matrix_p(:, 1)
1701 s_matrix => matrix_s(1, 1)%matrix
1702 DO ispin = 1, nspin
1703 CALL tb_ao_charges_matrix(matrix_p(ispin, 1)%matrix, s_matrix, aocg_spin(:, :, ispin), para_env)
1704 IF (do_dipole) THEN
1705 DO im = 1, dip_n
1706 CALL tb_contract_dens_matrix(matrix_p(ispin, 1)%matrix, tb%dipbra(im)%matrix, &
1707 tb%dipket(im)%matrix, ao_dip_spin(:, im, ispin), para_env)
1708 END DO
1709 END IF
1710 IF (do_quadrupole) THEN
1711 DO im = 1, quad_n
1712 CALL tb_contract_dens_matrix(matrix_p(ispin, 1)%matrix, tb%quadbra(im)%matrix, &
1713 tb%quadket(im)%matrix, ao_quad_spin(:, im, ispin), para_env)
1714 END DO
1715 END IF
1716 END DO
1717 END IF
1718 IF (nspin == 1) THEN
1719 aocg(:, :) = aocg_spin(:, :, 1)
1720 IF (do_dipole) ao_dip(:, :) = ao_dip_spin(:, :, 1)
1721 IF (do_quadrupole) ao_quad(:, :) = ao_quad_spin(:, :, 1)
1722 ELSE
1723 aocg(:, :) = aocg_spin(:, :, 1) + aocg_spin(:, :, 2)
1724 IF (do_dipole) THEN
1725 DO im = 1, dip_n
1726 DO iatom = 1, n_atom
1727 pao = ao_dip_spin(iatom, im, 1)
1728 ao_dip_spin(iatom, im, 1) = pao + ao_dip_spin(iatom, im, 2)
1729 ao_dip_spin(iatom, im, 2) = pao - ao_dip_spin(iatom, im, 2)
1730 END DO
1731 END DO
1732 ao_dip(:, :) = ao_dip_spin(:, :, 1)
1733 END IF
1734 IF (do_quadrupole) THEN
1735 DO im = 1, quad_n
1736 DO iatom = 1, n_atom
1737 pao = ao_quad_spin(iatom, im, 1)
1738 ao_quad_spin(iatom, im, 1) = pao + ao_quad_spin(iatom, im, 2)
1739 ao_quad_spin(iatom, im, 2) = pao - ao_quad_spin(iatom, im, 2)
1740 END DO
1741 END DO
1742 ao_quad(:, :) = ao_quad_spin(:, :, 1)
1743 END IF
1744 END IF
1745 NULLIFY (xtb_kind)
1746 DO ikind = 1, nkind
1747 CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
1748 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1749 CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, nao=nao, occupation=occ)
1750 DO iatom = 1, na
1751 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1752 DO is = 1, natorb
1753 ns = lao(is) + 1
1754 norm = 2*lao(is) + 1
1755 ch_ref(is, atom_a) = tb%calc%h0%refocc(nao(is), ikind)/norm
1756 ch_orb(is, atom_a) = aocg(is, atom_a) - ch_ref(is, atom_a)
1757 ch_orb_spin(is, atom_a, 1) = ch_orb(is, atom_a)
1758 IF (nspin == 2) ch_orb_spin(is, atom_a, 2) = &
1759 aocg_spin(is, atom_a, 1) - aocg_spin(is, atom_a, 2)
1760 ch_shell(atom_a, ns) = ch_orb(is, atom_a) + ch_shell(atom_a, ns)
1761 DO ispin = 1, nspin
1762 ch_shell_spin(atom_a, ns, ispin) = ch_orb_spin(is, atom_a, ispin) + &
1763 ch_shell_spin(atom_a, ns, ispin)
1764 END DO
1765 END DO
1766 DO ispin = 1, nspin
1767 ch_atom(atom_a, ispin) = sum(ch_orb_spin(:, atom_a, ispin))
1768 END DO
1769 END DO
1770 END DO
1771 native_seed_charge = -sum(ch_atom(:, 1))
1772 seed_native_from_rho = sum(abs(aocg_spin)) > 1.0e-10_dp .AND. &
1773 abs(native_seed_charge - real(dft_control%charge, dp)) < 1.0e-5_dp
1774 DEALLOCATE (aocg, aocg_spin)
1775
1776 raw_state_status = 1
1777#if defined(__TBLITE_DEBUG_DIAGNOSTICS)
1778 CALL get_environment_variable("CP2K_TBLITE_RAW_STATE_DUMP", raw_state_file, status=raw_state_status)
1779#endif
1780 IF (raw_state_status == 0) THEN
1781 OPEN (newunit=raw_state_unit, file=trim(raw_state_file), status="REPLACE", action="WRITE")
1782 WRITE (raw_state_unit, *) "qat"
1783 DO iatom = 1, n_atom
1784 WRITE (raw_state_unit, "(I0,1X,ES24.16)") iatom, -ch_atom(iatom, 1)
1785 END DO
1786 WRITE (raw_state_unit, *) "qsh"
1787 DO iatom = 1, n_atom
1788 DO is = 1, tb%calc%bas%nsh_at(iatom)
1789 WRITE (raw_state_unit, "(I0,1X,ES24.16)") tb%calc%bas%ish_at(iatom) + is, -ch_shell(iatom, is)
1790 END DO
1791 END DO
1792 IF (do_dipole) THEN
1793 WRITE (raw_state_unit, *) "dpat"
1794 DO iatom = 1, n_atom
1795 WRITE (raw_state_unit, "(I0,3(1X,ES24.16))") iatom, -ao_dip(iatom, :)
1796 END DO
1797 END IF
1798 IF (do_quadrupole) THEN
1799 WRITE (raw_state_unit, *) "qpat"
1800 DO iatom = 1, n_atom
1801 WRITE (raw_state_unit, "(I0,6(1X,ES24.16))") iatom, -ao_quad(iatom, :)
1802 END DO
1803 END IF
1804 CLOSE (raw_state_unit)
1805 END IF
1806
1807 IF (use_native_mixer) THEN
1808 IF (.NOT. ALLOCATED(tb%mixer)) cpabort("tblite mixer not initialized")
1809 advance_native_mixer = .false.
1810 IF (use_rho .AND. (.NOT. calculate_forces)) THEN
1811 advance_native_mixer = scf_env%iter_count > 1
1812 END IF
1813 IF (advance_native_mixer) THEN
1814 CALL tb%mixer%next(error)
1815 IF (ALLOCATED(error)) cpabort("tblite native mixer failed")
1816 CALL tb%mixer%get(tb%wfn%qsh)
1817 tb%wfn%qat(:, :) = 0.0_dp
1818 DO iatom = 1, n_atom
1819 ii = tb%calc%bas%ish_at(iatom)
1820 DO ispin = 1, nspin
1821 tb%wfn%qat(iatom, ispin) = &
1822 sum(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), ispin))
1823 END DO
1824 END DO
1825 IF (do_dipole) THEN
1826 CALL tb%mixer%get(tb%wfn%dpat)
1827 DEALLOCATE (ao_dip)
1828 END IF
1829 IF (do_quadrupole) THEN
1830 CALL tb%mixer%get(tb%wfn%qpat)
1831 DEALLOCATE (ao_quad)
1832 END IF
1833 ELSE
1834 IF (use_rho .AND. (.NOT. calculate_forces)) THEN
1835 IF (.NOT. reuse_native_input) THEN
1836 IF (seed_native_from_rho) THEN
1837 ! Seed the native tblite mixer from CP2K's current density so SCF_GUESS/RESTART
1838 ! controls the initial SCC state.
1839 DO iatom = 1, n_atom
1840 ii = tb%calc%bas%ish_at(iatom)
1841 DO ispin = 1, nspin
1842 DO is = 1, tb%calc%bas%nsh_at(iatom)
1843 tb%wfn%qsh(ii + is, ispin) = -ch_shell_spin(iatom, is, ispin)
1844 END DO
1845 tb%wfn%qat(iatom, ispin) = &
1846 sum(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), ispin))
1847 END DO
1848 END DO
1849 IF (do_dipole) THEN
1850 DO iatom = 1, n_atom
1851 DO ispin = 1, nspin
1852 tb%wfn%dpat(:, iatom, ispin) = -ao_dip_spin(iatom, :, ispin)
1853 END DO
1854 END DO
1855 END IF
1856 IF (do_quadrupole) THEN
1857 DO iatom = 1, n_atom
1858 DO ispin = 1, nspin
1859 tb%wfn%qpat(:, iatom, ispin) = -ao_quad_spin(iatom, :, ispin)
1860 END DO
1861 END DO
1862 END IF
1863 ELSE
1864 tb%wfn%qsh(:, :) = 0.0_dp
1865 tb%wfn%qat(:, :) = 0.0_dp
1866 IF (do_dipole) tb%wfn%dpat(:, :, :) = 0.0_dp
1867 IF (do_quadrupole) tb%wfn%qpat(:, :, :) = 0.0_dp
1868 END IF
1869 END IF
1870 IF (do_dipole) DEALLOCATE (ao_dip)
1871 IF (do_quadrupole) DEALLOCATE (ao_quad)
1872 ELSE
1873 IF ((.NOT. use_rho) .AND. (.NOT. calculate_forces)) THEN
1874 CALL tb%mixer%set(tb%wfn%qsh)
1875 IF (do_dipole) CALL tb%mixer%set(tb%wfn%dpat)
1876 IF (do_quadrupole) CALL tb%mixer%set(tb%wfn%qpat)
1877 END IF
1878 DO iatom = 1, n_atom
1879 ii = tb%calc%bas%ish_at(iatom)
1880 DO ispin = 1, nspin
1881 DO is = 1, tb%calc%bas%nsh_at(iatom)
1882 tb%wfn%qsh(ii + is, ispin) = -ch_shell_spin(iatom, is, ispin)
1883 END DO
1884 tb%wfn%qat(iatom, ispin) = &
1885 sum(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), ispin))
1886 END DO
1887 END DO
1888 IF (do_dipole) THEN
1889 DO iatom = 1, n_atom
1890 DO ispin = 1, nspin
1891 tb%wfn%dpat(:, iatom, ispin) = -ao_dip_spin(iatom, :, ispin)
1892 END DO
1893 END DO
1894 DEALLOCATE (ao_dip)
1895 END IF
1896 IF (do_quadrupole) THEN
1897 DO iatom = 1, n_atom
1898 DO ispin = 1, nspin
1899 tb%wfn%qpat(:, iatom, ispin) = -ao_quad_spin(iatom, :, ispin)
1900 END DO
1901 END DO
1902 DEALLOCATE (ao_quad)
1903 END IF
1904 IF ((.NOT. use_rho) .AND. (.NOT. calculate_forces)) THEN
1905 CALL tb%mixer%diff(tb%wfn%qsh)
1906 IF (do_dipole) CALL tb%mixer%diff(tb%wfn%dpat)
1907 IF (do_quadrupole) CALL tb%mixer%diff(tb%wfn%qpat)
1908 END IF
1909 END IF
1910 END IF
1911 ELSE
1912 ! charge mixing
1913 native_sign_mixing = .false.
1914 IF (.NOT. (dft_control%qs_control%do_ls_scf .OR. scf_control%use_ot)) THEN
1915 IF (.NOT. ASSOCIATED(scf_env)) cpabort("CP2K SCC mixer requires a QS SCF environment")
1916 IF (.NOT. (do_dipole .OR. do_quadrupole) .AND. &
1917 scf_env%mixing_method == modified_broyden_mixing_nr) THEN
1918 cpabort("MODIFIED_BROYDEN_MIXING with SCC_MIXER CP2K requires GFN2")
1919 END IF
1920 END IF
1921 IF (dft_control%qs_control%do_ls_scf .OR. scf_control%use_ot) THEN
1922 ! LS_SCF and OT optimize the electronic variables directly. Use the
1923 ! current shell/multipole response from that density without an
1924 ! extra SCC variable mixing step.
1925 ELSE IF (nspin > 1) THEN
1926 n_mix_cols = nspin*max_shell
1927 IF (do_dipole) n_mix_cols = n_mix_cols + nspin*dip_n
1928 IF (do_quadrupole) n_mix_cols = n_mix_cols + nspin*quad_n
1929 ALLOCATE (mix_vars(n_atom, n_mix_cols))
1930 mix_vars = 0.0_dp
1931
1932 mix_offset = 0
1933 DO ispin = 1, nspin
1934 mix_vars(:, mix_offset + 1:mix_offset + max_shell) = &
1935 -ch_shell_spin(:, 1:max_shell, ispin)
1936 mix_offset = mix_offset + max_shell
1937 END DO
1938 IF (do_dipole) THEN
1939 DO ispin = 1, nspin
1940 mix_vars(:, mix_offset + 1:mix_offset + dip_n) = &
1941 -ao_dip_spin(:, 1:dip_n, ispin)
1942 mix_offset = mix_offset + dip_n
1943 END DO
1944 END IF
1945 IF (do_quadrupole) THEN
1946 DO ispin = 1, nspin
1947 mix_vars(:, mix_offset + 1:mix_offset + quad_n) = &
1948 -ao_quad_spin(:, 1:quad_n, ispin)
1949 mix_offset = mix_offset + quad_n
1950 END DO
1951 END IF
1952
1953 IF (.NOT. use_no_mixer) THEN
1954 CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
1955 mix_vars, para_env, scf_env%iter_count)
1956 END IF
1957
1958 mix_offset = 0
1959 DO ispin = 1, nspin
1960 ch_shell_spin(:, 1:max_shell, ispin) = &
1961 -mix_vars(:, mix_offset + 1:mix_offset + max_shell)
1962 mix_offset = mix_offset + max_shell
1963 END DO
1964 ch_shell(:, 1:max_shell) = ch_shell_spin(:, 1:max_shell, 1)
1965 IF (do_dipole) THEN
1966 DO ispin = 1, nspin
1967 ao_dip_spin(:, 1:dip_n, ispin) = &
1968 -mix_vars(:, mix_offset + 1:mix_offset + dip_n)
1969 mix_offset = mix_offset + dip_n
1970 END DO
1971 ao_dip(:, 1:dip_n) = ao_dip_spin(:, 1:dip_n, 1)
1972 END IF
1973 IF (do_quadrupole) THEN
1974 DO ispin = 1, nspin
1975 ao_quad_spin(:, 1:quad_n, ispin) = &
1976 -mix_vars(:, mix_offset + 1:mix_offset + quad_n)
1977 mix_offset = mix_offset + quad_n
1978 END DO
1979 ao_quad(:, 1:quad_n) = ao_quad_spin(:, 1:quad_n, 1)
1980 END IF
1981 DEALLOCATE (mix_vars)
1982 ELSE
1983 do_combined_mixing = do_dipole .OR. do_quadrupole
1984 native_sign_mixing = do_dipole .OR. do_quadrupole
1985 discard_mixed_output = .false.
1986 skip_charge_mixing = use_no_mixer
1987 IF (skip_charge_mixing) THEN
1988 !
1989 ELSE IF (do_combined_mixing) THEN
1990 n_mix_cols = max_shell
1991 IF (do_dipole) n_mix_cols = n_mix_cols + dip_n
1992 IF (do_quadrupole) n_mix_cols = n_mix_cols + quad_n
1993 ALLOCATE (mix_vars(n_atom, n_mix_cols))
1994 IF (native_sign_mixing) THEN
1995 mix_vars(:, 1:max_shell) = -ch_shell(:, 1:max_shell)
1996 ELSE
1997 mix_vars(:, 1:max_shell) = ch_shell(:, 1:max_shell)
1998 END IF
1999 mix_offset = max_shell
2000 IF (do_dipole) THEN
2001 IF (native_sign_mixing) THEN
2002 mix_vars(:, mix_offset + 1:mix_offset + dip_n) = -ao_dip(:, 1:dip_n)
2003 ELSE
2004 mix_vars(:, mix_offset + 1:mix_offset + dip_n) = ao_dip(:, 1:dip_n)
2005 END IF
2006 mix_offset = mix_offset + dip_n
2007 END IF
2008 IF (do_quadrupole) THEN
2009 IF (native_sign_mixing) THEN
2010 mix_vars(:, mix_offset + 1:mix_offset + quad_n) = -ao_quad(:, 1:quad_n)
2011 ELSE
2012 mix_vars(:, mix_offset + 1:mix_offset + quad_n) = ao_quad(:, 1:quad_n)
2013 END IF
2014 END IF
2015 CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
2016 mix_vars, para_env, scf_env%iter_count)
2017 IF (.NOT. discard_mixed_output) THEN
2018 IF (native_sign_mixing) THEN
2019 ch_shell(:, 1:max_shell) = -mix_vars(:, 1:max_shell)
2020 ELSE
2021 ch_shell(:, 1:max_shell) = mix_vars(:, 1:max_shell)
2022 END IF
2023 mix_offset = max_shell
2024 IF (do_dipole) THEN
2025 IF (native_sign_mixing) THEN
2026 ao_dip(:, 1:dip_n) = -mix_vars(:, mix_offset + 1:mix_offset + dip_n)
2027 ELSE
2028 ao_dip(:, 1:dip_n) = mix_vars(:, mix_offset + 1:mix_offset + dip_n)
2029 END IF
2030 mix_offset = mix_offset + dip_n
2031 END IF
2032 IF (do_quadrupole) THEN
2033 IF (native_sign_mixing) THEN
2034 ao_quad(:, 1:quad_n) = -mix_vars(:, mix_offset + 1:mix_offset + quad_n)
2035 ELSE
2036 ao_quad(:, 1:quad_n) = mix_vars(:, mix_offset + 1:mix_offset + quad_n)
2037 END IF
2038 END IF
2039 END IF
2040 DEALLOCATE (mix_vars)
2041 ELSE
2042 CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
2043 ch_shell, para_env, scf_env%iter_count)
2044 END IF
2045 END IF
2046
2047 !setting new wave function
2048 CALL tb%pot%reset
2049 tb%e_es = 0.0_dp
2050 tb%e_scd = 0.0_dp
2051 IF (nspin > 1) THEN
2052 DO iatom = 1, n_atom
2053 ii = tb%calc%bas%ish_at(iatom)
2054 DO ispin = 1, nspin
2055 DO is = 1, tb%calc%bas%nsh_at(iatom)
2056 tb%wfn%qsh(ii + is, ispin) = -ch_shell_spin(iatom, is, ispin)
2057 END DO
2058 tb%wfn%qat(iatom, ispin) = &
2059 sum(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), ispin))
2060 END DO
2061 END DO
2062 IF (do_dipole) THEN
2063 DO iatom = 1, n_atom
2064 DO ispin = 1, nspin
2065 tb%wfn%dpat(:, iatom, ispin) = -ao_dip_spin(iatom, :, ispin)
2066 END DO
2067 END DO
2068 DEALLOCATE (ao_dip)
2069 END IF
2070 IF (do_quadrupole) THEN
2071 DO iatom = 1, n_atom
2072 DO ispin = 1, nspin
2073 tb%wfn%qpat(:, iatom, ispin) = -ao_quad_spin(iatom, :, ispin)
2074 END DO
2075 END DO
2076 DEALLOCATE (ao_quad)
2077 END IF
2078 ELSE
2079 DO iatom = 1, n_atom
2080 ii = tb%calc%bas%ish_at(iatom)
2081 DO is = 1, tb%calc%bas%nsh_at(iatom)
2082 new_charge = -ch_shell(iatom, is)
2083 tb%wfn%qsh(ii + is, 1) = new_charge
2084 END DO
2085 IF (native_sign_mixing) THEN
2086 tb%wfn%qat(iatom, 1) = sum(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), 1))
2087 ELSE
2088 tb%wfn%qat(iatom, 1) = -ch_atom(iatom, 1)
2089 END IF
2090 END DO
2091
2092 IF (do_dipole) THEN
2093 DO iatom = 1, n_atom
2094 DO im = 1, dip_n
2095 tb%wfn%dpat(im, iatom, 1) = -ao_dip(iatom, im)
2096 END DO
2097 END DO
2098 DEALLOCATE (ao_dip)
2099 END IF
2100 IF (do_quadrupole) THEN
2101 DO iatom = 1, n_atom
2102 DO im = 1, quad_n
2103 tb%wfn%qpat(im, iatom, 1) = -ao_quad(iatom, im)
2104 END DO
2105 END DO
2106 DEALLOCATE (ao_quad)
2107 END IF
2108 END IF
2109 END IF
2110
2111 CALL tb%pot%reset
2112 tb%e_es = 0.0_dp
2113 tb%e_scd = 0.0_dp
2114 tb%e_int = 0.0_dp
2115 IF (ALLOCATED(tb%calc%coulomb)) THEN
2116 CALL tb%calc%coulomb%get_potential(tb%mol, tb%cache, tb%wfn, tb%pot)
2117 CALL tb%calc%coulomb%get_energy(tb%mol, tb%cache, tb%wfn, tb%e_es)
2118 END IF
2119 IF (ALLOCATED(tb%calc%dispersion)) THEN
2120 IF (.NOT. skip_scf_dispersion) THEN
2121 CALL tb%calc%dispersion%get_potential(tb%mol, tb%dcache, tb%wfn, tb%pot)
2122 CALL tb%calc%dispersion%get_energy(tb%mol, tb%dcache, tb%wfn, tb%e_scd)
2123 END IF
2124 END IF
2125 IF (ALLOCATED(tb%calc%interactions)) THEN
2126 CALL tb%calc%interactions%get_potential(tb%mol, tb%icache, tb%wfn, tb%pot)
2127 CALL tb%calc%interactions%get_energy(tb%mol, tb%icache, tb%wfn, tb%e_int)
2128 END IF
2129
2130 IF (calculate_forces) THEN
2131 IF (ALLOCATED(tb%calc%coulomb)) THEN
2132 tb%grad = 0.0_dp
2133 CALL tb%calc%coulomb%get_gradient(tb%mol, tb%cache, tb%wfn, tb%grad, tb%sigma)
2134 CALL tb_dump_sigma_component("after_coulomb", tb%sigma, para_env)
2135 CALL tb_grad2force(qs_env, tb, para_env, 3)
2136 END IF
2137
2138 IF (ALLOCATED(tb%calc%dispersion) .AND. .NOT. skip_scf_dispersion) THEN
2139 tb%grad = 0.0_dp
2140 CALL tb%calc%dispersion%get_gradient(tb%mol, tb%dcache, tb%wfn, tb%grad, tb%sigma)
2141 CALL tb_dump_sigma_component("after_dispersion_scf", tb%sigma, para_env)
2142 CALL tb_grad2force(qs_env, tb, para_env, 2)
2143 END IF
2144
2145 IF (ALLOCATED(tb%calc%interactions)) THEN
2146 tb%grad = 0.0_dp
2147 CALL tb%calc%interactions%get_gradient(tb%mol, tb%icache, tb%wfn, tb%grad, tb%sigma)
2148 CALL tb_dump_sigma_component("after_interactions_scf", tb%sigma, para_env)
2149 CALL tb_grad2force(qs_env, tb, para_env, 3)
2150 END IF
2151 END IF
2152
2153 IF (ALLOCATED(ao_dip_spin)) DEALLOCATE (ao_dip_spin)
2154 IF (ALLOCATED(ao_quad_spin)) DEALLOCATE (ao_quad_spin)
2155 DEALLOCATE (ch_atom, ch_shell, ch_orb, ch_ref, ch_orb_spin, ch_shell_spin)
2156
2157#else
2158 mark_used(qs_env)
2159 mark_used(tb)
2160 mark_used(dft_control)
2161 mark_used(calculate_forces)
2162 mark_used(use_rho)
2163 cpabort("Built without TBLITE")
2164#endif
2165
2166 END SUBROUTINE tb_update_charges
2167
2168! **************************************************************************************************
2169!> \brief ...
2170!> \param qs_env ...
2171!> \param tb ...
2172!> \param dft_control ...
2173! **************************************************************************************************
2174 SUBROUTINE tb_ham_add_coulomb(qs_env, tb, dft_control)
2175
2176 TYPE(qs_environment_type), POINTER :: qs_env
2177 TYPE(tblite_type), POINTER :: tb
2178 TYPE(dft_control_type), POINTER :: dft_control
2179
2180#if defined(__TBLITE)
2181
2182 INTEGER :: ikind, jkind, iatom, jatom, icol, irow
2183 INTEGER :: ic, id1, id2, id3, iq1, iq2, iq3, iq4, iq5, iq6, &
2184 is, nimg, ni, nj, i, j, nspin
2185 INTEGER :: la, lb, za, zb
2186 LOGICAL :: found
2187 INTEGER, DIMENSION(3) :: cellind
2188 INTEGER, DIMENSION(25) :: naoa, naob
2189 REAL(kind=dp), DIMENSION(3) :: rij
2190 REAL(kind=dp) :: mpfac
2191#if defined(__TBLITE_DEBUG_DIAGNOSTICS)
2192 INTEGER :: debug_status
2193 CHARACTER(LEN=32) :: debug_value
2194#endif
2195 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, sum_shell
2196 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ashift, bshift
2197 REAL(kind=dp), DIMENSION(:, :), POINTER :: ksblock, sblock
2198 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2199 REAL(kind=dp), DIMENSION(:, :), POINTER :: dip_ket1, dip_ket2, dip_ket3, &
2200 dip_bra1, dip_bra2, dip_bra3
2201 REAL(kind=dp), DIMENSION(:, :), POINTER :: quad_ket1, quad_ket2, quad_ket3, &
2202 quad_ket4, quad_ket5, quad_ket6, &
2203 quad_bra1, quad_bra2, quad_bra3, &
2204 quad_bra4, quad_bra5, quad_bra6
2205
2206 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2207 TYPE(dbcsr_iterator_type) :: iter
2208 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
2209 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
2210 TYPE(kpoint_type), POINTER :: kpoints
2211 TYPE(neighbor_list_iterator_p_type), &
2212 DIMENSION(:), POINTER :: nl_iterator
2213 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2214 POINTER :: n_list
2215 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2216 POINTER :: kp_list
2217 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2218 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
2219
2220 nimg = dft_control%nimages
2221 mpfac = -0.5_dp
2222
2223 NULLIFY (matrix_s, ks_matrix, n_list, kp_list, qs_kind_set)
2224 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list, sab_kp=kp_list, &
2225 matrix_s_kp=matrix_s, matrix_ks_kp=ks_matrix, qs_kind_set=qs_kind_set)
2226 IF (nimg > 1) THEN
2227 IF (.NOT. ASSOCIATED(kp_list)) cpabort("Missing k-point neighbor list for tblite Hamiltonian")
2228 n_list => kp_list
2229 END IF
2230 nspin = SIZE(ks_matrix, 1)
2231
2232 !creating sum of shell lists
2233 ALLOCATE (sum_shell(tb%mol%nat))
2234 i = 0
2235 DO j = 1, tb%mol%nat
2236 sum_shell(j) = i
2237 i = i + tb%calc%bas%nsh_at(j)
2238 END DO
2239
2240 IF (nimg == 1) THEN
2241 ! no k-points; all matrices have been transformed to periodic bsf
2242 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
2243 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
2244 kind_of=kind_of)
2245 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
2246 DO WHILE (dbcsr_iterator_blocks_left(iter))
2247 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
2248
2249 ikind = kind_of(irow)
2250 jkind = kind_of(icol)
2251
2252 ! atomic parameters
2253 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
2254 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
2255 CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa)
2256 CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob)
2257
2258 ni = SIZE(sblock, 1)
2259 ALLOCATE (ashift(ni, ni))
2260
2261 nj = SIZE(sblock, 2)
2262 ALLOCATE (bshift(nj, nj))
2263
2264 DO is = 1, nspin
2265 ashift = 0.0_dp
2266 DO i = 1, ni
2267 la = naoa(i) + sum_shell(irow)
2268 ashift(i, i) = tb_spin_project(tb%pot%vsh(la, :), is)
2269 END DO
2270 bshift = 0.0_dp
2271 DO j = 1, nj
2272 lb = naob(j) + sum_shell(icol)
2273 bshift(j, j) = tb_spin_project(tb%pot%vsh(lb, :), is)
2274 END DO
2275 NULLIFY (ksblock)
2276 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
2277 row=irow, col=icol, block=ksblock, found=found)
2278 cpassert(found)
2279 ksblock = ksblock - 0.5_dp*(matmul(ashift, sblock) &
2280 + matmul(sblock, bshift))
2281 ksblock = ksblock - 0.5_dp*(tb_spin_project(tb%pot%vat(irow, :), is) &
2282 + tb_spin_project(tb%pot%vat(icol, :), is))*sblock
2283 END DO
2284 DEALLOCATE (ashift, bshift)
2285 END DO
2286 CALL dbcsr_iterator_stop(iter)
2287
2288 IF (ASSOCIATED(tb%dipbra)) THEN
2289 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
2290 DO WHILE (dbcsr_iterator_blocks_left(iter))
2291 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
2292
2293 NULLIFY (dip_bra1, dip_bra2, dip_bra3)
2294 CALL dbcsr_get_block_p(matrix=tb%dipbra(1)%matrix, &
2295 row=irow, col=icol, block=dip_bra1, found=found)
2296 cpassert(found)
2297 CALL dbcsr_get_block_p(matrix=tb%dipbra(2)%matrix, &
2298 row=irow, col=icol, block=dip_bra2, found=found)
2299 cpassert(found)
2300 CALL dbcsr_get_block_p(matrix=tb%dipbra(3)%matrix, &
2301 row=irow, col=icol, block=dip_bra3, found=found)
2302 cpassert(found)
2303 NULLIFY (dip_ket1, dip_ket2, dip_ket3)
2304 CALL dbcsr_get_block_p(matrix=tb%dipket(1)%matrix, &
2305 row=irow, col=icol, block=dip_ket1, found=found)
2306 cpassert(found)
2307 CALL dbcsr_get_block_p(matrix=tb%dipket(2)%matrix, &
2308 row=irow, col=icol, block=dip_ket2, found=found)
2309 cpassert(found)
2310 CALL dbcsr_get_block_p(matrix=tb%dipket(3)%matrix, &
2311 row=irow, col=icol, block=dip_ket3, found=found)
2312 cpassert(found)
2313
2314 DO is = 1, nspin
2315 NULLIFY (ksblock)
2316 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
2317 row=irow, col=icol, block=ksblock, found=found)
2318 cpassert(found)
2319 ksblock = ksblock + mpfac*(dip_ket1*tb_spin_project(tb%pot%vdp(1, irow, :), is) &
2320 + dip_ket2*tb_spin_project(tb%pot%vdp(2, irow, :), is) &
2321 + dip_ket3*tb_spin_project(tb%pot%vdp(3, irow, :), is) &
2322 + dip_bra1*tb_spin_project(tb%pot%vdp(1, icol, :), is) &
2323 + dip_bra2*tb_spin_project(tb%pot%vdp(2, icol, :), is) &
2324 + dip_bra3*tb_spin_project(tb%pot%vdp(3, icol, :), is))
2325 END DO
2326 END DO
2327 CALL dbcsr_iterator_stop(iter)
2328 END IF
2329
2330 IF (ASSOCIATED(tb%quadbra)) THEN
2331 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
2332 DO WHILE (dbcsr_iterator_blocks_left(iter))
2333 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
2334
2335 NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
2336 CALL dbcsr_get_block_p(matrix=tb%quadbra(1)%matrix, &
2337 row=irow, col=icol, block=quad_bra1, found=found)
2338 cpassert(found)
2339 CALL dbcsr_get_block_p(matrix=tb%quadbra(2)%matrix, &
2340 row=irow, col=icol, block=quad_bra2, found=found)
2341 cpassert(found)
2342 CALL dbcsr_get_block_p(matrix=tb%quadbra(3)%matrix, &
2343 row=irow, col=icol, block=quad_bra3, found=found)
2344 cpassert(found)
2345 CALL dbcsr_get_block_p(matrix=tb%quadbra(4)%matrix, &
2346 row=irow, col=icol, block=quad_bra4, found=found)
2347 cpassert(found)
2348 CALL dbcsr_get_block_p(matrix=tb%quadbra(5)%matrix, &
2349 row=irow, col=icol, block=quad_bra5, found=found)
2350 cpassert(found)
2351 CALL dbcsr_get_block_p(matrix=tb%quadbra(6)%matrix, &
2352 row=irow, col=icol, block=quad_bra6, found=found)
2353 cpassert(found)
2354
2355 NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
2356 CALL dbcsr_get_block_p(matrix=tb%quadket(1)%matrix, &
2357 row=irow, col=icol, block=quad_ket1, found=found)
2358 cpassert(found)
2359 CALL dbcsr_get_block_p(matrix=tb%quadket(2)%matrix, &
2360 row=irow, col=icol, block=quad_ket2, found=found)
2361 cpassert(found)
2362 CALL dbcsr_get_block_p(matrix=tb%quadket(3)%matrix, &
2363 row=irow, col=icol, block=quad_ket3, found=found)
2364 cpassert(found)
2365 CALL dbcsr_get_block_p(matrix=tb%quadket(4)%matrix, &
2366 row=irow, col=icol, block=quad_ket4, found=found)
2367 cpassert(found)
2368 CALL dbcsr_get_block_p(matrix=tb%quadket(5)%matrix, &
2369 row=irow, col=icol, block=quad_ket5, found=found)
2370 cpassert(found)
2371 CALL dbcsr_get_block_p(matrix=tb%quadket(6)%matrix, &
2372 row=irow, col=icol, block=quad_ket6, found=found)
2373 cpassert(found)
2374
2375 DO is = 1, nspin
2376 NULLIFY (ksblock)
2377 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
2378 row=irow, col=icol, block=ksblock, found=found)
2379 cpassert(found)
2380
2381 ksblock = ksblock + mpfac*(quad_ket1*tb_spin_project(tb%pot%vqp(1, irow, :), is) &
2382 + quad_ket2*tb_spin_project(tb%pot%vqp(2, irow, :), is) &
2383 + quad_ket3*tb_spin_project(tb%pot%vqp(3, irow, :), is) &
2384 + quad_ket4*tb_spin_project(tb%pot%vqp(4, irow, :), is) &
2385 + quad_ket5*tb_spin_project(tb%pot%vqp(5, irow, :), is) &
2386 + quad_ket6*tb_spin_project(tb%pot%vqp(6, irow, :), is) &
2387 + quad_bra1*tb_spin_project(tb%pot%vqp(1, icol, :), is) &
2388 + quad_bra2*tb_spin_project(tb%pot%vqp(2, icol, :), is) &
2389 + quad_bra3*tb_spin_project(tb%pot%vqp(3, icol, :), is) &
2390 + quad_bra4*tb_spin_project(tb%pot%vqp(4, icol, :), is) &
2391 + quad_bra5*tb_spin_project(tb%pot%vqp(5, icol, :), is) &
2392 + quad_bra6*tb_spin_project(tb%pot%vqp(6, icol, :), is))
2393 END DO
2394 END DO
2395 CALL dbcsr_iterator_stop(iter)
2396 END IF
2397
2398 ELSE
2399 NULLIFY (kpoints)
2400 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
2401 NULLIFY (cell_to_index)
2402 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
2403
2404 NULLIFY (nl_iterator)
2405 CALL neighbor_list_iterator_create(nl_iterator, n_list)
2406 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2407 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2408 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
2409
2410 icol = max(iatom, jatom)
2411 irow = min(iatom, jatom)
2412
2413 IF (iatom > jatom) THEN
2414 i = ikind
2415 ikind = jkind
2416 jkind = i
2417 END IF
2418
2419 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2420 cpassert(ic > 0)
2421
2422 NULLIFY (sblock)
2423 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
2424 row=irow, col=icol, block=sblock, found=found)
2425 cpassert(found)
2426
2427 ! atomic parameters
2428 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
2429 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
2430 CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa)
2431 CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob)
2432
2433 ni = SIZE(sblock, 1)
2434 ALLOCATE (ashift(ni, ni))
2435
2436 nj = SIZE(sblock, 2)
2437 ALLOCATE (bshift(nj, nj))
2438
2439 DO is = 1, nspin
2440 ashift = 0.0_dp
2441 DO i = 1, ni
2442 la = naoa(i) + sum_shell(irow)
2443 ashift(i, i) = tb_spin_project(tb%pot%vsh(la, :), is)
2444 END DO
2445 bshift = 0.0_dp
2446 DO j = 1, nj
2447 lb = naob(j) + sum_shell(icol)
2448 bshift(j, j) = tb_spin_project(tb%pot%vsh(lb, :), is)
2449 END DO
2450 NULLIFY (ksblock)
2451 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
2452 row=irow, col=icol, block=ksblock, found=found)
2453 cpassert(found)
2454 ksblock = ksblock - 0.5_dp*(matmul(ashift, sblock) &
2455 + matmul(sblock, bshift))
2456 ksblock = ksblock - 0.5_dp*(tb_spin_project(tb%pot%vat(irow, :), is) &
2457 + tb_spin_project(tb%pot%vat(icol, :), is))*sblock
2458 END DO
2459 DEALLOCATE (ashift, bshift)
2460 END DO
2461 CALL neighbor_list_iterator_release(nl_iterator)
2462
2463 IF (ASSOCIATED(tb%dipbra)) THEN
2464 NULLIFY (nl_iterator)
2465 CALL neighbor_list_iterator_create(nl_iterator, n_list)
2466 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2467 CALL get_iterator_info(nl_iterator, &
2468 iatom=iatom, jatom=jatom, cell=cellind)
2469 icol = max(iatom, jatom)
2470 irow = min(iatom, jatom)
2471 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2472 cpassert(ic > 0)
2473 id1 = 1 + dip_n*(ic - 1)
2474 id2 = 2 + dip_n*(ic - 1)
2475 id3 = 3 + dip_n*(ic - 1)
2476
2477 NULLIFY (dip_bra1, dip_bra2, dip_bra3)
2478 CALL dbcsr_get_block_p(matrix=tb%dipbra(id1)%matrix, &
2479 row=irow, col=icol, block=dip_bra1, found=found)
2480 cpassert(found)
2481 CALL dbcsr_get_block_p(matrix=tb%dipbra(id2)%matrix, &
2482 row=irow, col=icol, block=dip_bra2, found=found)
2483 cpassert(found)
2484 CALL dbcsr_get_block_p(matrix=tb%dipbra(id3)%matrix, &
2485 row=irow, col=icol, block=dip_bra3, found=found)
2486 cpassert(found)
2487 NULLIFY (dip_ket1, dip_ket2, dip_ket3)
2488 CALL dbcsr_get_block_p(matrix=tb%dipket(id1)%matrix, &
2489 row=irow, col=icol, block=dip_ket1, found=found)
2490 cpassert(found)
2491 CALL dbcsr_get_block_p(matrix=tb%dipket(id2)%matrix, &
2492 row=irow, col=icol, block=dip_ket2, found=found)
2493 cpassert(found)
2494 CALL dbcsr_get_block_p(matrix=tb%dipket(id3)%matrix, &
2495 row=irow, col=icol, block=dip_ket3, found=found)
2496 cpassert(found)
2497
2498 DO is = 1, nspin
2499 NULLIFY (ksblock)
2500 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
2501 row=irow, col=icol, block=ksblock, found=found)
2502 cpassert(found)
2503 ksblock = ksblock + mpfac*(dip_ket1*tb_spin_project(tb%pot%vdp(1, irow, :), is) &
2504 + dip_ket2*tb_spin_project(tb%pot%vdp(2, irow, :), is) &
2505 + dip_ket3*tb_spin_project(tb%pot%vdp(3, irow, :), is) &
2506 + dip_bra1*tb_spin_project(tb%pot%vdp(1, icol, :), is) &
2507 + dip_bra2*tb_spin_project(tb%pot%vdp(2, icol, :), is) &
2508 + dip_bra3*tb_spin_project(tb%pot%vdp(3, icol, :), is))
2509 END DO
2510 END DO
2511 CALL neighbor_list_iterator_release(nl_iterator)
2512 END IF
2513
2514 IF (ASSOCIATED(tb%quadbra)) THEN
2515 NULLIFY (nl_iterator)
2516 CALL neighbor_list_iterator_create(nl_iterator, n_list)
2517 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2518 CALL get_iterator_info(nl_iterator, &
2519 iatom=iatom, jatom=jatom, cell=cellind)
2520 icol = max(iatom, jatom)
2521 irow = min(iatom, jatom)
2522 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2523 cpassert(ic > 0)
2524 iq1 = 1 + quad_n*(ic - 1)
2525 iq2 = 2 + quad_n*(ic - 1)
2526 iq3 = 3 + quad_n*(ic - 1)
2527 iq4 = 4 + quad_n*(ic - 1)
2528 iq5 = 5 + quad_n*(ic - 1)
2529 iq6 = 6 + quad_n*(ic - 1)
2530
2531 NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
2532 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq1)%matrix, &
2533 row=irow, col=icol, block=quad_bra1, found=found)
2534 cpassert(found)
2535 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq2)%matrix, &
2536 row=irow, col=icol, block=quad_bra2, found=found)
2537 cpassert(found)
2538 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq3)%matrix, &
2539 row=irow, col=icol, block=quad_bra3, found=found)
2540 cpassert(found)
2541 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq4)%matrix, &
2542 row=irow, col=icol, block=quad_bra4, found=found)
2543 cpassert(found)
2544 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq5)%matrix, &
2545 row=irow, col=icol, block=quad_bra5, found=found)
2546 cpassert(found)
2547 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq6)%matrix, &
2548 row=irow, col=icol, block=quad_bra6, found=found)
2549 cpassert(found)
2550
2551 NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
2552 CALL dbcsr_get_block_p(matrix=tb%quadket(iq1)%matrix, &
2553 row=irow, col=icol, block=quad_ket1, found=found)
2554 cpassert(found)
2555 CALL dbcsr_get_block_p(matrix=tb%quadket(iq2)%matrix, &
2556 row=irow, col=icol, block=quad_ket2, found=found)
2557 cpassert(found)
2558 CALL dbcsr_get_block_p(matrix=tb%quadket(iq3)%matrix, &
2559 row=irow, col=icol, block=quad_ket3, found=found)
2560 cpassert(found)
2561 CALL dbcsr_get_block_p(matrix=tb%quadket(iq4)%matrix, &
2562 row=irow, col=icol, block=quad_ket4, found=found)
2563 cpassert(found)
2564 CALL dbcsr_get_block_p(matrix=tb%quadket(iq5)%matrix, &
2565 row=irow, col=icol, block=quad_ket5, found=found)
2566 cpassert(found)
2567 CALL dbcsr_get_block_p(matrix=tb%quadket(iq6)%matrix, &
2568 row=irow, col=icol, block=quad_ket6, found=found)
2569 cpassert(found)
2570
2571 DO is = 1, nspin
2572 NULLIFY (ksblock)
2573 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
2574 row=irow, col=icol, block=ksblock, found=found)
2575 cpassert(found)
2576
2577 ksblock = ksblock + mpfac*(quad_ket1*tb_spin_project(tb%pot%vqp(1, irow, :), is) &
2578 + quad_ket2*tb_spin_project(tb%pot%vqp(2, irow, :), is) &
2579 + quad_ket3*tb_spin_project(tb%pot%vqp(3, irow, :), is) &
2580 + quad_ket4*tb_spin_project(tb%pot%vqp(4, irow, :), is) &
2581 + quad_ket5*tb_spin_project(tb%pot%vqp(5, irow, :), is) &
2582 + quad_ket6*tb_spin_project(tb%pot%vqp(6, irow, :), is) &
2583 + quad_bra1*tb_spin_project(tb%pot%vqp(1, icol, :), is) &
2584 + quad_bra2*tb_spin_project(tb%pot%vqp(2, icol, :), is) &
2585 + quad_bra3*tb_spin_project(tb%pot%vqp(3, icol, :), is) &
2586 + quad_bra4*tb_spin_project(tb%pot%vqp(4, icol, :), is) &
2587 + quad_bra5*tb_spin_project(tb%pot%vqp(5, icol, :), is) &
2588 + quad_bra6*tb_spin_project(tb%pot%vqp(6, icol, :), is))
2589 END DO
2590 END DO
2591 CALL neighbor_list_iterator_release(nl_iterator)
2592 END IF
2593
2594 END IF
2595
2596#else
2597 mark_used(qs_env)
2598 mark_used(tb)
2599 mark_used(dft_control)
2600 cpabort("Built without TBLITE")
2601#endif
2602
2603 END SUBROUTINE tb_ham_add_coulomb
2604
2605! **************************************************************************************************
2606!> \brief ...
2607!> \param qs_env ...
2608!> \param tb ...
2609! **************************************************************************************************
2610 SUBROUTINE tb_get_multipole(qs_env, tb)
2611
2612 TYPE(qs_environment_type), POINTER :: qs_env
2613 TYPE(tblite_type), POINTER :: tb
2614
2615#if defined(__TBLITE)
2616
2617 CHARACTER(LEN=*), PARAMETER :: routinen = 'tb_get_multipole'
2618
2619 INTEGER :: ikind, jkind, iatom, jatom, icol, irow, iset, jset, ityp, jtyp
2620 INTEGER :: ic, idx, id1, id2, id3, img, iq1, iq2, iq3, iq4, iq5, iq6
2621 INTEGER :: nkind, natom, handle, nimg, i, inda, indb
2622 INTEGER :: atom_a, atom_b, nseta, nsetb, ia, ib, ij
2623 LOGICAL :: found
2624 REAL(kind=dp) :: r2
2625 INTEGER, DIMENSION(3) :: cell
2626 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2627 REAL(kind=dp), DIMENSION(3) :: rij
2628 INTEGER, DIMENSION(:), POINTER :: la_max, lb_max
2629 INTEGER, DIMENSION(:), POINTER :: nsgfa, nsgfb
2630 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
2631 INTEGER, ALLOCATABLE :: atom_of_kind(:)
2632 REAL(kind=dp), ALLOCATABLE :: stmp(:)
2633 REAL(kind=dp), ALLOCATABLE :: dtmp(:, :), qtmp(:, :), dtmpj(:, :), qtmpj(:, :)
2634 REAL(kind=dp), DIMENSION(:, :), POINTER :: dip_ket1, dip_ket2, dip_ket3, &
2635 dip_bra1, dip_bra2, dip_bra3
2636 REAL(kind=dp), DIMENSION(:, :), POINTER :: quad_ket1, quad_ket2, quad_ket3, &
2637 quad_ket4, quad_ket5, quad_ket6, &
2638 quad_bra1, quad_bra2, quad_bra3, &
2639 quad_bra4, quad_bra5, quad_bra6
2640
2641 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2642 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
2643 TYPE(dft_control_type), POINTER :: dft_control
2644 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
2645 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
2646 TYPE(kpoint_type), POINTER :: kpoints
2647 TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
2648 TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_kp
2649 TYPE(neighbor_list_iterator_p_type), &
2650 DIMENSION(:), POINTER :: nl_iterator
2651 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2652 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2653
2654 CALL timeset(routinen, handle)
2655
2656 !get info from environment vaiarable
2657 NULLIFY (atomic_kind_set, qs_kind_set, sab_orb, sab_kp, particle_set)
2658 NULLIFY (dft_control, matrix_s, kpoints, cell_to_index)
2659 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
2660 qs_kind_set=qs_kind_set, &
2661 sab_orb=sab_orb, &
2662 sab_kp=sab_kp, &
2663 particle_set=particle_set, &
2664 dft_control=dft_control, &
2665 kpoints=kpoints, &
2666 matrix_s_kp=matrix_s)
2667 natom = SIZE(particle_set)
2668 nkind = SIZE(atomic_kind_set)
2669 nimg = dft_control%nimages
2670 IF (nimg > 1) THEN
2671 IF (.NOT. ASSOCIATED(sab_kp)) cpabort("Missing k-point neighbor list for tblite multipoles")
2672 sab_orb => sab_kp
2673 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
2674 END IF
2675
2676 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
2677
2678 !set up basis set lists
2679 ALLOCATE (basis_set_list(nkind))
2680 CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
2681
2682 ALLOCATE (stmp(msao(tb%calc%bas%maxl)**2))
2683 ALLOCATE (dtmp(dip_n, msao(tb%calc%bas%maxl)**2))
2684 ALLOCATE (qtmp(quad_n, msao(tb%calc%bas%maxl)**2))
2685 ALLOCATE (dtmpj(dip_n, msao(tb%calc%bas%maxl)**2))
2686 ALLOCATE (qtmpj(quad_n, msao(tb%calc%bas%maxl)**2))
2687
2688 ! allocate dipole/quadrupole moment matrix elemnts
2689 CALL dbcsr_allocate_matrix_set(tb%dipbra, dip_n*nimg)
2690 CALL dbcsr_allocate_matrix_set(tb%dipket, dip_n*nimg)
2691 CALL dbcsr_allocate_matrix_set(tb%quadbra, quad_n*nimg)
2692 CALL dbcsr_allocate_matrix_set(tb%quadket, quad_n*nimg)
2693 DO img = 1, nimg
2694 DO i = 1, dip_n
2695 idx = i + dip_n*(img - 1)
2696 ALLOCATE (tb%dipbra(idx)%matrix)
2697 ALLOCATE (tb%dipket(idx)%matrix)
2698 CALL dbcsr_create(tb%dipbra(idx)%matrix, template=matrix_s(1, img)%matrix, &
2699 name="DIPOLE BRAMATRIX")
2700 CALL dbcsr_create(tb%dipket(idx)%matrix, template=matrix_s(1, img)%matrix, &
2701 name="DIPOLE KETMATRIX")
2702 CALL cp_dbcsr_alloc_block_from_nbl(tb%dipbra(idx)%matrix, sab_orb)
2703 CALL cp_dbcsr_alloc_block_from_nbl(tb%dipket(idx)%matrix, sab_orb)
2704 END DO
2705 DO i = 1, quad_n
2706 idx = i + quad_n*(img - 1)
2707 ALLOCATE (tb%quadbra(idx)%matrix)
2708 ALLOCATE (tb%quadket(idx)%matrix)
2709 CALL dbcsr_create(tb%quadbra(idx)%matrix, template=matrix_s(1, img)%matrix, &
2710 name="QUADRUPOLE BRAMATRIX")
2711 CALL dbcsr_create(tb%quadket(idx)%matrix, template=matrix_s(1, img)%matrix, &
2712 name="QUADRUPOLE KETMATRIX")
2713 CALL cp_dbcsr_alloc_block_from_nbl(tb%quadbra(idx)%matrix, sab_orb)
2714 CALL cp_dbcsr_alloc_block_from_nbl(tb%quadket(idx)%matrix, sab_orb)
2715 END DO
2716 END DO
2717
2718 !loop over all atom pairs with a non-zero overlap (sab_orb)
2719 NULLIFY (nl_iterator)
2720 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
2721 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2722 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2723 iatom=iatom, jatom=jatom, r=rij, cell=cell)
2724
2725 r2 = norm2(rij(:))**2
2726
2727 icol = max(iatom, jatom)
2728 irow = min(iatom, jatom)
2729
2730 IF (iatom < jatom) THEN
2731 rij = -rij
2732 i = ikind
2733 ikind = jkind
2734 jkind = i
2735 END IF
2736
2737 IF (nimg == 1) THEN
2738 ic = 1
2739 ELSE
2740 ic = cell_to_index(cell(1), cell(2), cell(3))
2741 cpassert(ic > 0)
2742 END IF
2743 id1 = 1 + dip_n*(ic - 1)
2744 id2 = 2 + dip_n*(ic - 1)
2745 id3 = 3 + dip_n*(ic - 1)
2746 iq1 = 1 + quad_n*(ic - 1)
2747 iq2 = 2 + quad_n*(ic - 1)
2748 iq3 = 3 + quad_n*(ic - 1)
2749 iq4 = 4 + quad_n*(ic - 1)
2750 iq5 = 5 + quad_n*(ic - 1)
2751 iq6 = 6 + quad_n*(ic - 1)
2752
2753 ityp = tb%mol%id(icol)
2754 jtyp = tb%mol%id(irow)
2755
2756 NULLIFY (dip_bra1, dip_bra2, dip_bra3)
2757 CALL dbcsr_get_block_p(matrix=tb%dipbra(id1)%matrix, &
2758 row=irow, col=icol, block=dip_bra1, found=found)
2759 cpassert(found)
2760 CALL dbcsr_get_block_p(matrix=tb%dipbra(id2)%matrix, &
2761 row=irow, col=icol, block=dip_bra2, found=found)
2762 cpassert(found)
2763 CALL dbcsr_get_block_p(matrix=tb%dipbra(id3)%matrix, &
2764 row=irow, col=icol, block=dip_bra3, found=found)
2765 cpassert(found)
2766
2767 NULLIFY (dip_ket1, dip_ket2, dip_ket3)
2768 CALL dbcsr_get_block_p(matrix=tb%dipket(id1)%matrix, &
2769 row=irow, col=icol, block=dip_ket1, found=found)
2770 cpassert(found)
2771 CALL dbcsr_get_block_p(matrix=tb%dipket(id2)%matrix, &
2772 row=irow, col=icol, block=dip_ket2, found=found)
2773 cpassert(found)
2774 CALL dbcsr_get_block_p(matrix=tb%dipket(id3)%matrix, &
2775 row=irow, col=icol, block=dip_ket3, found=found)
2776 cpassert(found)
2777
2778 NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
2779 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq1)%matrix, &
2780 row=irow, col=icol, block=quad_bra1, found=found)
2781 cpassert(found)
2782 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq2)%matrix, &
2783 row=irow, col=icol, block=quad_bra2, found=found)
2784 cpassert(found)
2785 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq3)%matrix, &
2786 row=irow, col=icol, block=quad_bra3, found=found)
2787 cpassert(found)
2788 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq4)%matrix, &
2789 row=irow, col=icol, block=quad_bra4, found=found)
2790 cpassert(found)
2791 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq5)%matrix, &
2792 row=irow, col=icol, block=quad_bra5, found=found)
2793 cpassert(found)
2794 CALL dbcsr_get_block_p(matrix=tb%quadbra(iq6)%matrix, &
2795 row=irow, col=icol, block=quad_bra6, found=found)
2796 cpassert(found)
2797
2798 NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
2799 CALL dbcsr_get_block_p(matrix=tb%quadket(iq1)%matrix, &
2800 row=irow, col=icol, block=quad_ket1, found=found)
2801 cpassert(found)
2802 CALL dbcsr_get_block_p(matrix=tb%quadket(iq2)%matrix, &
2803 row=irow, col=icol, block=quad_ket2, found=found)
2804 cpassert(found)
2805 CALL dbcsr_get_block_p(matrix=tb%quadket(iq3)%matrix, &
2806 row=irow, col=icol, block=quad_ket3, found=found)
2807 cpassert(found)
2808 CALL dbcsr_get_block_p(matrix=tb%quadket(iq4)%matrix, &
2809 row=irow, col=icol, block=quad_ket4, found=found)
2810 cpassert(found)
2811 CALL dbcsr_get_block_p(matrix=tb%quadket(iq5)%matrix, &
2812 row=irow, col=icol, block=quad_ket5, found=found)
2813 cpassert(found)
2814 CALL dbcsr_get_block_p(matrix=tb%quadket(iq6)%matrix, &
2815 row=irow, col=icol, block=quad_ket6, found=found)
2816 cpassert(found)
2817
2818 !get basis information
2819 basis_set_a => basis_set_list(ikind)%gto_basis_set
2820 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
2821 basis_set_b => basis_set_list(jkind)%gto_basis_set
2822 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
2823 atom_a = atom_of_kind(icol)
2824 atom_b = atom_of_kind(irow)
2825 ! basis a
2826 first_sgfa => basis_set_a%first_sgf
2827 la_max => basis_set_a%lmax
2828 nseta = basis_set_a%nset
2829 nsgfa => basis_set_a%nsgf_set
2830 ! basis b
2831 first_sgfb => basis_set_b%first_sgf
2832 lb_max => basis_set_b%lmax
2833 nsetb = basis_set_b%nset
2834 nsgfb => basis_set_b%nsgf_set
2835
2836 ! --------- Hamiltonian
2837 ! Periodic self-images are off-diagonal lattice contributions, not the on-site block.
2838 IF (icol == irow .AND. r2 < same_atom**2) THEN
2839 DO iset = 1, nseta
2840 DO jset = 1, nsetb
2841 CALL multipole_cgto(tb%calc%bas%cgto(jset, ityp), tb%calc%bas%cgto(iset, ityp), &
2842 & r2, -rij, tb%calc%bas%intcut, stmp, dtmp, qtmp)
2843
2844 DO inda = 1, nsgfa(iset)
2845 ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
2846 DO indb = 1, nsgfb(jset)
2847 ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
2848 ij = indb + nsgfb(jset)*(inda - 1)
2849
2850 dip_ket1(ib, ia) = dip_ket1(ib, ia) + dtmp(1, ij)
2851 dip_ket2(ib, ia) = dip_ket2(ib, ia) + dtmp(2, ij)
2852 dip_ket3(ib, ia) = dip_ket3(ib, ia) + dtmp(3, ij)
2853
2854 quad_ket1(ib, ia) = quad_ket1(ib, ia) + qtmp(1, ij)
2855 quad_ket2(ib, ia) = quad_ket2(ib, ia) + qtmp(2, ij)
2856 quad_ket3(ib, ia) = quad_ket3(ib, ia) + qtmp(3, ij)
2857 quad_ket4(ib, ia) = quad_ket4(ib, ia) + qtmp(4, ij)
2858 quad_ket5(ib, ia) = quad_ket5(ib, ia) + qtmp(5, ij)
2859 quad_ket6(ib, ia) = quad_ket6(ib, ia) + qtmp(6, ij)
2860
2861 dip_bra1(ib, ia) = dip_bra1(ib, ia) + dtmp(1, ij)
2862 dip_bra2(ib, ia) = dip_bra2(ib, ia) + dtmp(2, ij)
2863 dip_bra3(ib, ia) = dip_bra3(ib, ia) + dtmp(3, ij)
2864
2865 quad_bra1(ib, ia) = quad_bra1(ib, ia) + qtmp(1, ij)
2866 quad_bra2(ib, ia) = quad_bra2(ib, ia) + qtmp(2, ij)
2867 quad_bra3(ib, ia) = quad_bra3(ib, ia) + qtmp(3, ij)
2868 quad_bra4(ib, ia) = quad_bra4(ib, ia) + qtmp(4, ij)
2869 quad_bra5(ib, ia) = quad_bra5(ib, ia) + qtmp(5, ij)
2870 quad_bra6(ib, ia) = quad_bra6(ib, ia) + qtmp(6, ij)
2871 END DO
2872 END DO
2873 END DO
2874 END DO
2875 ELSE
2876 DO iset = 1, nseta
2877 DO jset = 1, nsetb
2878 CALL multipole_cgto(tb%calc%bas%cgto(jset, jtyp), tb%calc%bas%cgto(iset, ityp), &
2879 & r2, -rij, tb%calc%bas%intcut, stmp, dtmp, qtmp)
2880
2881 DO inda = 1, nsgfa(iset)
2882 ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
2883 DO indb = 1, nsgfb(jset)
2884 ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
2885
2886 ij = indb + nsgfb(jset)*(inda - 1)
2887 CALL tb_shift_multipole(-rij, stmp(ij), dtmp(:, ij), qtmp(:, ij), &
2888 dtmpj(:, ij), qtmpj(:, ij))
2889
2890 dip_bra1(ib, ia) = dip_bra1(ib, ia) + dtmp(1, ij)
2891 dip_bra2(ib, ia) = dip_bra2(ib, ia) + dtmp(2, ij)
2892 dip_bra3(ib, ia) = dip_bra3(ib, ia) + dtmp(3, ij)
2893
2894 quad_bra1(ib, ia) = quad_bra1(ib, ia) + qtmp(1, ij)
2895 quad_bra2(ib, ia) = quad_bra2(ib, ia) + qtmp(2, ij)
2896 quad_bra3(ib, ia) = quad_bra3(ib, ia) + qtmp(3, ij)
2897 quad_bra4(ib, ia) = quad_bra4(ib, ia) + qtmp(4, ij)
2898 quad_bra5(ib, ia) = quad_bra5(ib, ia) + qtmp(5, ij)
2899 quad_bra6(ib, ia) = quad_bra6(ib, ia) + qtmp(6, ij)
2900
2901 dip_ket1(ib, ia) = dip_ket1(ib, ia) + dtmpj(1, ij)
2902 dip_ket2(ib, ia) = dip_ket2(ib, ia) + dtmpj(2, ij)
2903 dip_ket3(ib, ia) = dip_ket3(ib, ia) + dtmpj(3, ij)
2904
2905 quad_ket1(ib, ia) = quad_ket1(ib, ia) + qtmpj(1, ij)
2906 quad_ket2(ib, ia) = quad_ket2(ib, ia) + qtmpj(2, ij)
2907 quad_ket3(ib, ia) = quad_ket3(ib, ia) + qtmpj(3, ij)
2908 quad_ket4(ib, ia) = quad_ket4(ib, ia) + qtmpj(4, ij)
2909 quad_ket5(ib, ia) = quad_ket5(ib, ia) + qtmpj(5, ij)
2910 quad_ket6(ib, ia) = quad_ket6(ib, ia) + qtmpj(6, ij)
2911 END DO
2912 END DO
2913 END DO
2914 END DO
2915 END IF
2916 END DO
2917 CALL neighbor_list_iterator_release(nl_iterator)
2918
2919 DO i = 1, SIZE(tb%dipbra)
2920 CALL dbcsr_finalize(tb%dipbra(i)%matrix)
2921 CALL dbcsr_finalize(tb%dipket(i)%matrix)
2922 END DO
2923 DO i = 1, SIZE(tb%quadbra)
2924 CALL dbcsr_finalize(tb%quadbra(i)%matrix)
2925 CALL dbcsr_finalize(tb%quadket(i)%matrix)
2926 END DO
2927
2928 DEALLOCATE (basis_set_list)
2929
2930 CALL timestop(handle)
2931
2932#else
2933 mark_used(qs_env)
2934 mark_used(tb)
2935 cpabort("Built without TBLITE")
2936#endif
2937
2938 END SUBROUTINE tb_get_multipole
2939
2940! **************************************************************************************************
2941!> \brief Shift a multipole operator from one center to the other.
2942!> \param vec displacement vector between the two centers
2943!> \param s overlap integral
2944!> \param di dipole integral on the original center
2945!> \param qi quadrupole integral on the original center
2946!> \param dj dipole integral on the shifted center
2947!> \param qj quadrupole integral on the shifted center
2948! **************************************************************************************************
2949 PURE SUBROUTINE tb_shift_multipole(vec, s, di, qi, dj, qj)
2950
2951 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: vec
2952 REAL(kind=dp), INTENT(IN) :: s
2953 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: di, qi
2954 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: dj, qj
2955
2956 REAL(kind=dp) :: tr
2957
2958 dj(1) = di(1) + vec(1)*s
2959 dj(2) = di(2) + vec(2)*s
2960 dj(3) = di(3) + vec(3)*s
2961
2962 qj(1) = 2*vec(1)*di(1) + vec(1)**2*s
2963 qj(3) = 2*vec(2)*di(2) + vec(2)**2*s
2964 qj(6) = 2*vec(3)*di(3) + vec(3)**2*s
2965 qj(2) = vec(1)*di(2) + vec(2)*di(1) + vec(1)*vec(2)*s
2966 qj(4) = vec(1)*di(3) + vec(3)*di(1) + vec(1)*vec(3)*s
2967 qj(5) = vec(2)*di(3) + vec(3)*di(2) + vec(2)*vec(3)*s
2968 tr = 0.5_dp*(qj(1) + qj(3) + qj(6))
2969
2970 qj(1) = qi(1) + 1.5_dp*qj(1) - tr
2971 qj(2) = qi(2) + 1.5_dp*qj(2)
2972 qj(3) = qi(3) + 1.5_dp*qj(3) - tr
2973 qj(4) = qi(4) + 1.5_dp*qj(4)
2974 qj(5) = qi(5) + 1.5_dp*qj(5)
2975 qj(6) = qi(6) + 1.5_dp*qj(6) - tr
2976
2977 END SUBROUTINE tb_shift_multipole
2978
2979! **************************************************************************************************
2980!> \brief compute the mulliken properties (AO resolved)
2981!> \param p_mat ...
2982!> \param s_matrix ...
2983!> \param charges ...
2984!> \param para_env ...
2985! **************************************************************************************************
2986 SUBROUTINE tb_ao_charges_matrix(p_mat, s_matrix, charges, para_env)
2987 TYPE(dbcsr_type), POINTER :: p_mat, s_matrix
2988 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
2989 TYPE(mp_para_env_type), POINTER :: para_env
2990
2991 INTEGER :: i, iblock_col, iblock_row, j
2992 LOGICAL :: found
2993 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, s_block
2994 TYPE(dbcsr_iterator_type) :: iter
2995
2996 charges = 0.0_dp
2997 CALL dbcsr_iterator_start(iter, s_matrix)
2998 DO WHILE (dbcsr_iterator_blocks_left(iter))
2999 NULLIFY (s_block, p_block)
3000 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, s_block)
3001 CALL dbcsr_get_block_p(matrix=p_mat, row=iblock_row, col=iblock_col, block=p_block, found=found)
3002 IF (.NOT. found) cycle
3003 IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) cycle
3004
3005 DO j = 1, SIZE(p_block, 2)
3006 DO i = 1, SIZE(p_block, 1)
3007 charges(i, iblock_row) = charges(i, iblock_row) + p_block(i, j)*s_block(i, j)
3008 END DO
3009 END DO
3010 IF (iblock_col /= iblock_row) THEN
3011 DO j = 1, SIZE(p_block, 2)
3012 DO i = 1, SIZE(p_block, 1)
3013 charges(j, iblock_col) = charges(j, iblock_col) + p_block(i, j)*s_block(i, j)
3014 END DO
3015 END DO
3016 END IF
3017 END DO
3018 CALL dbcsr_iterator_stop(iter)
3019 CALL para_env%sum(charges)
3020
3021 END SUBROUTINE tb_ao_charges_matrix
3022
3023! **************************************************************************************************
3024!> \brief compute the AO-resolved Mulliken charges for one k-point spin channel.
3025!> \param p_matrix_kp ...
3026!> \param s_matrix_kp ...
3027!> \param charges ...
3028!> \param ispin ...
3029!> \param para_env ...
3030! **************************************************************************************************
3031 SUBROUTINE tb_ao_charges_kp_spin(p_matrix_kp, s_matrix_kp, charges, ispin, para_env)
3032 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_matrix_kp, s_matrix_kp
3033 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
3034 INTEGER, INTENT(IN) :: ispin
3035 TYPE(mp_para_env_type), POINTER :: para_env
3036
3037 INTEGER :: ic
3038 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: image_charges
3039 TYPE(dbcsr_type), POINTER :: p_mat, s_mat
3040
3041 charges = 0.0_dp
3042 ALLOCATE (image_charges(SIZE(charges, 1), SIZE(charges, 2)))
3043 DO ic = 1, SIZE(s_matrix_kp, 2)
3044 NULLIFY (p_mat, s_mat)
3045 p_mat => p_matrix_kp(ispin, ic)%matrix
3046 s_mat => s_matrix_kp(1, ic)%matrix
3047 IF (ASSOCIATED(p_mat) .AND. ASSOCIATED(s_mat)) THEN
3048 image_charges = 0.0_dp
3049 CALL tb_ao_charges_matrix(p_mat, s_mat, image_charges, para_env)
3050 charges(:, :) = charges(:, :) + image_charges(:, :)
3051 END IF
3052 END DO
3053 DEALLOCATE (image_charges)
3054
3055 END SUBROUTINE tb_ao_charges_kp_spin
3056
3057! **************************************************************************************************
3058!> \brief compute the mulliken properties (AO resolved)
3059!> \param p_mat ...
3060!> \param bra_mat ...
3061!> \param ket_mat ...
3062!> \param output ...
3063!> \param para_env ...
3064! **************************************************************************************************
3065 SUBROUTINE tb_contract_dens_matrix(p_mat, bra_mat, ket_mat, output, para_env)
3066 TYPE(dbcsr_type), POINTER :: p_mat, bra_mat, ket_mat
3067 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: output
3068 TYPE(mp_para_env_type), POINTER :: para_env
3069
3070 INTEGER :: i, iblock_col, iblock_row, j
3071 LOGICAL :: found
3072 REAL(kind=dp), DIMENSION(:, :), POINTER :: bra, ket, p_block
3073 TYPE(dbcsr_iterator_type) :: iter
3074
3075 output = 0.0_dp
3076 CALL dbcsr_iterator_start(iter, bra_mat)
3077 DO WHILE (dbcsr_iterator_blocks_left(iter))
3078 NULLIFY (p_block, bra, ket)
3079 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, bra)
3080 CALL dbcsr_get_block_p(matrix=p_mat, row=iblock_row, col=iblock_col, block=p_block, found=found)
3081 IF (.NOT. found) cycle
3082 CALL dbcsr_get_block_p(matrix=ket_mat, row=iblock_row, col=iblock_col, block=ket, found=found)
3083 IF (.NOT. found) cpabort("missing block")
3084
3085 IF (.NOT. (ASSOCIATED(bra) .AND. ASSOCIATED(p_block))) cycle
3086 IF (iblock_col == iblock_row) THEN
3087 DO j = 1, SIZE(p_block, 1)
3088 DO i = 1, SIZE(p_block, 2)
3089 output(iblock_row) = output(iblock_row) + p_block(j, i)*bra(j, i)
3090 END DO
3091 END DO
3092 ELSE
3093 DO j = 1, SIZE(p_block, 1)
3094 DO i = 1, SIZE(p_block, 2)
3095 output(iblock_row) = output(iblock_row) + p_block(j, i)*ket(j, i)
3096 END DO
3097 END DO
3098 DO j = 1, SIZE(p_block, 1)
3099 DO i = 1, SIZE(p_block, 2)
3100 output(iblock_col) = output(iblock_col) + p_block(j, i)*bra(j, i)
3101 END DO
3102 END DO
3103 END IF
3104 END DO
3105 CALL dbcsr_iterator_stop(iter)
3106 CALL para_env%sum(output)
3107
3108 END SUBROUTINE tb_contract_dens_matrix
3109
3110! **************************************************************************************************
3111!> \brief compute the AO-resolved density contraction for one k-point spin channel.
3112!> \param p_matrix ...
3113!> \param bra_mat ...
3114!> \param ket_mat ...
3115!> \param iop ...
3116!> \param nops ...
3117!> \param output ...
3118!> \param ispin ...
3119!> \param para_env ...
3120! **************************************************************************************************
3121 SUBROUTINE tb_contract_dens_kp_spin(p_matrix, bra_mat, ket_mat, iop, nops, output, ispin, para_env)
3122 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_matrix
3123 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: bra_mat, ket_mat
3124 INTEGER, INTENT(IN) :: iop, nops
3125 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: output
3126 INTEGER, INTENT(IN) :: ispin
3127 TYPE(mp_para_env_type), POINTER :: para_env
3128
3129 INTEGER :: ic, idx, nimg
3130 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: image_output
3131 TYPE(dbcsr_type), POINTER :: p_mat
3132
3133 nimg = SIZE(p_matrix, 2)
3134 output = 0.0_dp
3135 ALLOCATE (image_output(SIZE(output)))
3136 DO ic = 1, nimg
3137 idx = iop + nops*(ic - 1)
3138 cpassert(idx <= SIZE(bra_mat))
3139 cpassert(idx <= SIZE(ket_mat))
3140 NULLIFY (p_mat)
3141 p_mat => p_matrix(ispin, ic)%matrix
3142 image_output = 0.0_dp
3143 CALL tb_contract_dens_matrix(p_mat, bra_mat(idx)%matrix, ket_mat(idx)%matrix, image_output, para_env)
3144 output = output + image_output
3145 END DO
3146 DEALLOCATE (image_output)
3147
3148 END SUBROUTINE tb_contract_dens_kp_spin
3149
3150! **************************************************************************************************
3151!> \brief compute the mulliken properties (AO resolved)
3152!> \param p_matrix ...
3153!> \param bra_mat ...
3154!> \param ket_mat ...
3155!> \param output ...
3156!> \param para_env ...
3157!> \par History
3158!> adapted from ao_charges_2
3159!> \note
3160! **************************************************************************************************
3161 SUBROUTINE contract_dens(p_matrix, bra_mat, ket_mat, output, para_env)
3162 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
3163 TYPE(dbcsr_type), POINTER :: bra_mat, ket_mat
3164 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: output
3165 TYPE(mp_para_env_type), POINTER :: para_env
3166
3167 CHARACTER(len=*), PARAMETER :: routinen = 'contract_dens'
3168
3169 INTEGER :: handle, i, iblock_col, iblock_row, &
3170 ispin, j, nspin
3171 LOGICAL :: found
3172 REAL(kind=dp), DIMENSION(:, :), POINTER :: bra, ket, p_block
3173 TYPE(dbcsr_iterator_type) :: iter
3174
3175 CALL timeset(routinen, handle)
3176
3177 nspin = SIZE(p_matrix)
3178 output = 0.0_dp
3179 DO ispin = 1, nspin
3180 CALL dbcsr_iterator_start(iter, bra_mat)
3181 DO WHILE (dbcsr_iterator_blocks_left(iter))
3182 NULLIFY (p_block, bra, ket)
3183
3184 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, bra)
3185 CALL dbcsr_get_block_p(matrix=p_matrix(ispin)%matrix, &
3186 row=iblock_row, col=iblock_col, block=p_block, found=found)
3187 IF (.NOT. found) cycle
3188 CALL dbcsr_get_block_p(matrix=ket_mat, &
3189 row=iblock_row, col=iblock_col, block=ket, found=found)
3190 IF (.NOT. found) cpabort("missing block")
3191
3192 IF (.NOT. (ASSOCIATED(bra) .AND. ASSOCIATED(p_block))) cycle
3193 IF (iblock_col == iblock_row) THEN
3194 DO j = 1, SIZE(p_block, 1)
3195 DO i = 1, SIZE(p_block, 2)
3196 output(iblock_row) = output(iblock_row) + p_block(j, i)*bra(j, i)
3197 END DO
3198 END DO
3199 ELSE
3200 DO j = 1, SIZE(p_block, 1)
3201 DO i = 1, SIZE(p_block, 2)
3202 output(iblock_row) = output(iblock_row) + p_block(j, i)*ket(j, i)
3203 END DO
3204 END DO
3205 DO j = 1, SIZE(p_block, 1)
3206 DO i = 1, SIZE(p_block, 2)
3207 output(iblock_col) = output(iblock_col) + p_block(j, i)*bra(j, i)
3208 END DO
3209 END DO
3210 END IF
3211 END DO
3212 CALL dbcsr_iterator_stop(iter)
3213 END DO
3214
3215 CALL para_env%sum(output)
3216 CALL timestop(handle)
3217
3218 END SUBROUTINE contract_dens
3219
3220! **************************************************************************************************
3221!> \brief compute the AO-resolved density contraction for real-space k-point image matrices
3222!> \param p_matrix ...
3223!> \param bra_mat ...
3224!> \param ket_mat ...
3225!> \param iop ...
3226!> \param nops ...
3227!> \param output ...
3228!> \param para_env ...
3229! **************************************************************************************************
3230 SUBROUTINE contract_dens_kp(p_matrix, bra_mat, ket_mat, iop, nops, output, para_env)
3231 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_matrix
3232 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: bra_mat, ket_mat
3233 INTEGER, INTENT(IN) :: iop, nops
3234 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: output
3235 TYPE(mp_para_env_type), POINTER :: para_env
3236
3237 INTEGER :: ic, idx, nimg
3238 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: image_output
3239 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_image
3240
3241 nimg = SIZE(p_matrix, 2)
3242 output = 0.0_dp
3243 ALLOCATE (image_output(SIZE(output)))
3244 DO ic = 1, nimg
3245 idx = iop + nops*(ic - 1)
3246 cpassert(idx <= SIZE(bra_mat))
3247 cpassert(idx <= SIZE(ket_mat))
3248 NULLIFY (p_image)
3249 p_image => p_matrix(:, ic)
3250 image_output = 0.0_dp
3251 CALL contract_dens(p_image, bra_mat(idx)%matrix, ket_mat(idx)%matrix, image_output, para_env)
3252 output = output + image_output
3253 END DO
3254 DEALLOCATE (image_output)
3255
3256 END SUBROUTINE contract_dens_kp
3257
3258! **************************************************************************************************
3259!> \brief save gradient to force
3260!> \param qs_env ...
3261!> \param tb ...
3262!> \param para_env ...
3263!> \param ityp ...
3264!> \note
3265! **************************************************************************************************
3266 SUBROUTINE tb_grad2force(qs_env, tb, para_env, ityp)
3267
3268 TYPE(qs_environment_type) :: qs_env
3269 TYPE(tblite_type) :: tb
3270 TYPE(mp_para_env_type) :: para_env
3271 INTEGER :: ityp
3272
3273 CHARACTER(len=*), PARAMETER :: routinen = 'tb_grad2force'
3274
3275 CHARACTER(LEN=default_path_length) :: dump_file
3276 INTEGER :: atoma, dump_status, dump_unit, handle, &
3277 iatom, ikind, natom
3278 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
3279 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3280 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3281 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3282
3283 CALL timeset(routinen, handle)
3284
3285 NULLIFY (force, atomic_kind_set)
3286 CALL get_qs_env(qs_env=qs_env, force=force, particle_set=particle_set, &
3287 atomic_kind_set=atomic_kind_set)
3288 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
3289 atom_of_kind=atom_of_kind, kind_of=kind_of)
3290
3291 natom = SIZE(particle_set)
3292
3293 dump_status = 1
3294#if defined(__TBLITE_DEBUG_DIAGNOSTICS)
3295 CALL get_environment_variable("CP2K_TBLITE_FORCE_DUMP", dump_file, status=dump_status)
3296#endif
3297 IF (dump_status == 0) THEN
3298 OPEN (newunit=dump_unit, file=trim(dump_file), status="UNKNOWN", &
3299 position="APPEND", action="WRITE")
3300 WRITE (dump_unit, "(A,1X,I0)") "component", ityp
3301 DO iatom = 1, natom
3302 WRITE (dump_unit, "(I0,3(1X,ES24.16))") iatom, tb%grad(:, iatom)/para_env%num_pe
3303 END DO
3304 CLOSE (dump_unit)
3305 END IF
3306
3307 SELECT CASE (ityp)
3308 CASE DEFAULT
3309 cpabort("unknown force type")
3310 CASE (0)
3311 DO iatom = 1, natom
3312 ikind = kind_of(iatom)
3313 atoma = atom_of_kind(iatom)
3314 force(ikind)%all_potential(:, atoma) = &
3315 force(ikind)%all_potential(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3316 END DO
3317 CASE (1)
3318 DO iatom = 1, natom
3319 ikind = kind_of(iatom)
3320 atoma = atom_of_kind(iatom)
3321 force(ikind)%repulsive(:, atoma) = &
3322 force(ikind)%repulsive(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3323 END DO
3324 CASE (2)
3325 DO iatom = 1, natom
3326 ikind = kind_of(iatom)
3327 atoma = atom_of_kind(iatom)
3328 force(ikind)%dispersion(:, atoma) = &
3329 force(ikind)%dispersion(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3330 END DO
3331 CASE (3)
3332 DO iatom = 1, natom
3333 ikind = kind_of(iatom)
3334 atoma = atom_of_kind(iatom)
3335 force(ikind)%rho_elec(:, atoma) = &
3336 force(ikind)%rho_elec(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3337 END DO
3338 CASE (4)
3339 DO iatom = 1, natom
3340 ikind = kind_of(iatom)
3341 atoma = atom_of_kind(iatom)
3342 force(ikind)%overlap(:, atoma) = &
3343 force(ikind)%overlap(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3344 END DO
3345 CASE (5)
3346 DO iatom = 1, natom
3347 ikind = kind_of(iatom)
3348 atoma = atom_of_kind(iatom)
3349 force(ikind)%efield(:, atoma) = &
3350 force(ikind)%efield(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3351 END DO
3352 END SELECT
3353
3354 CALL timestop(handle)
3355
3356 END SUBROUTINE tb_grad2force
3357
3358! **************************************************************************************************
3359!> \brief set gradient to zero
3360!> \param qs_env ...
3361!> \note
3362! **************************************************************************************************
3363 SUBROUTINE tb_zero_force(qs_env)
3364
3365 TYPE(qs_environment_type) :: qs_env
3366
3367 INTEGER :: iatom, ikind, natom
3368 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
3369 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3370 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3371 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3372
3373 NULLIFY (force, atomic_kind_set)
3374 CALL get_qs_env(qs_env=qs_env, force=force, particle_set=particle_set, &
3375 atomic_kind_set=atomic_kind_set)
3376 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
3377 kind_of=kind_of)
3378
3379 natom = SIZE(particle_set)
3380
3381 DO iatom = 1, natom
3382 ikind = kind_of(iatom)
3383 force(ikind)%all_potential = 0.0_dp
3384 force(ikind)%repulsive = 0.0_dp
3385 force(ikind)%dispersion = 0.0_dp
3386 force(ikind)%rho_elec = 0.0_dp
3387 force(ikind)%overlap = 0.0_dp
3388 force(ikind)%efield = 0.0_dp
3389 END DO
3390
3391 END SUBROUTINE tb_zero_force
3392
3393! **************************************************************************************************
3394!> \brief Add SCC-overlap and direct multipole Hamiltonian derivatives.
3395!> \param qs_env ...
3396!> \param use_rho ...
3397!> \param nimg ...
3398! **************************************************************************************************
3399 SUBROUTINE tb_derive_dh_off(qs_env, use_rho, nimg)
3400
3401 TYPE(qs_environment_type), POINTER :: qs_env
3402 LOGICAL, INTENT(IN) :: use_rho
3403 INTEGER, INTENT(IN) :: nimg
3404
3405#if defined(__TBLITE)
3406 INTEGER :: i, idim, ij, iatom, ic, icol, ikind, img, ispin, &
3407 jdim, ni, nj, nkind, nel, &
3408 ityp, jatom, jkind, jrow, jtyp, iset, jset, nseti, nsetj, &
3409 ia, ib, inda, indb, sampled_axes, sampled_even_axes, &
3410 sampled_gamma_axes, nspin, ikp_axis
3411 INTEGER, DIMENSION(3) :: cellind, nkp_cellind, nkp_grid
3412 INTEGER, DIMENSION(:), POINTER :: nsgfa, nsgfb
3413 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
3414 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3415 LOGICAL :: found, gamma_centered, gamma_sampled_image_pair, &
3416 has_multipole_response, sampled_image_pair, &
3417 use_matrix_scc_stress
3418 LOGICAL, DIMENSION(3) :: mesh_has_gamma
3419 REAL(kind=dp) :: r2, dr, i_a_shift, j_a_shift, i_a_shift_mag, j_a_shift_mag, &
3420 ishift, jshift, ishift_mag, jshift_mag, pij_charge, &
3421 pij_magnet, mp_pair_scale, kpoint_coordinate, native_dot_tmp
3422 REAL(kind=dp), DIMENSION(3) :: kp_shift
3423 REAL(kind=dp), DIMENSION(3) :: rij, dgrad, dhgrad_charge, dhgrad_magnet, &
3424 mpgrad_charge, mpgrad_magnet
3425 REAL(kind=dp), DIMENSION(3, 3) :: hsigma
3426 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: t_ov, idip, jdip, idip_mag, jdip_mag, &
3427 iquad, jquad, iquad_mag, jquad_mag
3428 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: t_dip, t_quad, t_d_ov
3429 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: t_i_dip, t_i_quad, t_j_dip, t_j_quad
3430 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :, :) :: scc_strain_hint
3431 REAL(kind=dp), DIMENSION(:, :), POINTER :: pblock, pblock_beta
3432 TYPE(block_p_type), DIMENSION(3, 3, 2) :: scc_strain_blocks
3433
3434 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3435 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_dh_scc, matrix_p
3436 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
3437 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
3438 TYPE(kpoint_type), POINTER :: kpoints
3439 TYPE(mp_para_env_type), POINTER :: para_env
3440 TYPE(neighbor_list_iterator_p_type), &
3441 DIMENSION(:), POINTER :: nl_iterator
3442 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3443 POINTER :: sab_orb
3444 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3445 POINTER :: sab_kp
3446 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3447 TYPE(qs_rho_type), POINTER :: rho
3448 TYPE(qs_scf_env_type), POINTER :: scf_env
3449 TYPE(tblite_type), POINTER :: tb
3450 TYPE(virial_type), POINTER :: virial
3451
3452 ! compute mulliken charges required for charge update
3453 NULLIFY (scf_env, rho, tb, sab_orb, sab_kp, para_env, kpoints, matrix_dh_scc, virial)
3454 CALL get_qs_env(qs_env=qs_env, &
3455 atomic_kind_set=atomic_kind_set, &
3456 scf_env=scf_env, &
3457 rho=rho, &
3458 tb_tblite=tb, &
3459 sab_orb=sab_orb, &
3460 sab_kp=sab_kp, &
3461 para_env=para_env, &
3462 qs_kind_set=qs_kind_set, &
3463 kpoints=kpoints)
3464
3465 NULLIFY (cell_to_index)
3466 IF (nimg > 1) THEN
3467 IF (.NOT. ASSOCIATED(sab_kp)) cpabort("Missing tblite k-point neighbor list")
3468 sab_orb => sab_kp
3469 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index, nkp_grid=nkp_grid, &
3470 kp_shift=kp_shift, gamma_centered=gamma_centered)
3471 ELSE
3472 nkp_grid = 1
3473 kp_shift = 0.0_dp
3474 gamma_centered = .true.
3475 END IF
3476 mesh_has_gamma = .false.
3477 DO i = 1, 3
3478 DO ikp_axis = 1, nkp_grid(i)
3479 IF (gamma_centered .AND. modulo(nkp_grid(i), 2) == 0) THEN
3480 kpoint_coordinate = real(2*ikp_axis - nkp_grid(i), kind=dp)/ &
3481 REAL(2*nkp_grid(i), kind=dp) + kp_shift(i)
3482 ELSE
3483 kpoint_coordinate = real(2*ikp_axis - nkp_grid(i) - 1, kind=dp)/ &
3484 REAL(2*nkp_grid(i), kind=dp) + kp_shift(i)
3485 END IF
3486 IF (abs(kpoint_coordinate - anint(kpoint_coordinate)) < 1.0e-12_dp) THEN
3487 mesh_has_gamma(i) = .true.
3488 END IF
3489 END DO
3490 END DO
3491 has_multipole_response = ASSOCIATED(tb%dipbra) .OR. ASSOCIATED(tb%quadbra)
3492
3493 NULLIFY (matrix_p)
3494 IF (use_rho) THEN
3495 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
3496 ELSE IF (ASSOCIATED(tb%rho_ao_kp_ref)) THEN
3497 matrix_p => tb%rho_ao_kp_ref
3498 ELSE
3499 matrix_p => scf_env%p_mix_new
3500 END IF
3501 nspin = SIZE(matrix_p, 1)
3502 use_matrix_scc_stress = tb%use_virial .AND. nimg > 1 .AND. all(tb%mol%periodic)
3503
3504 IF (use_matrix_scc_stress) THEN
3505 CALL dbcsr_allocate_matrix_set(matrix_dh_scc, 9*nspin, nimg)
3506 DO img = 1, nimg
3507 DO ispin = 1, nspin
3508 DO jdim = 1, 3
3509 DO idim = 1, 3
3510 i = idim + 3*(jdim - 1) + 9*(ispin - 1)
3511 ALLOCATE (matrix_dh_scc(i, img)%matrix)
3512 CALL dbcsr_create(matrix_dh_scc(i, img)%matrix, template=matrix_p(ispin, img)%matrix, &
3513 name="TBLITE SCC STRAIN DERIVATIVE")
3514 CALL cp_dbcsr_alloc_block_from_nbl(matrix_dh_scc(i, img)%matrix, sab_orb)
3515 END DO
3516 END DO
3517 END DO
3518 END DO
3519 END IF
3520
3521 ! set up basis set lists
3522 nkind = SIZE(atomic_kind_set)
3523 ALLOCATE (basis_set_list(nkind))
3524 CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
3525
3526 nel = msao(tb%calc%bas%maxl)**2
3527 ALLOCATE (t_ov(nel))
3528 ALLOCATE (t_d_ov(3, nel))
3529 ALLOCATE (t_dip(dip_n, nel))
3530 ALLOCATE (t_i_dip(3, dip_n, nel), t_j_dip(3, dip_n, nel))
3531 ALLOCATE (t_quad(quad_n, nel))
3532 ALLOCATE (t_i_quad(3, quad_n, nel), t_j_quad(3, quad_n, nel))
3533
3534 ALLOCATE (idip(dip_n), jdip(dip_n), idip_mag(dip_n), jdip_mag(dip_n))
3535 ALLOCATE (iquad(quad_n), jquad(quad_n), iquad_mag(quad_n), jquad_mag(quad_n))
3536
3537 tb%grad = 0.0_dp
3538 hsigma = 0.0_dp
3539 ! loop over all atom pairs with a non-zero overlap (sab_orb)
3540 NULLIFY (nl_iterator)
3541 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
3542 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3543 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3544 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
3545
3546 icol = max(iatom, jatom)
3547 jrow = min(iatom, jatom)
3548
3549 IF (iatom < jatom) THEN
3550 rij = -rij
3551 i = ikind
3552 ikind = jkind
3553 jkind = i
3554 END IF
3555
3556 ityp = tb%mol%id(icol)
3557 jtyp = tb%mol%id(jrow)
3558
3559 r2 = dot_product(rij, rij)
3560 dr = sqrt(r2)
3561 IF (icol == jrow .AND. dr < same_atom) cycle
3562
3563 !get basis information
3564 basis_set_a => basis_set_list(ikind)%gto_basis_set
3565 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
3566 first_sgfa => basis_set_a%first_sgf
3567 nsgfa => basis_set_a%nsgf_set
3568 nseti = basis_set_a%nset
3569 basis_set_b => basis_set_list(jkind)%gto_basis_set
3570 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
3571 first_sgfb => basis_set_b%first_sgf
3572 nsgfb => basis_set_b%nsgf_set
3573 nsetj = basis_set_b%nset
3574
3575 IF (nimg == 1) THEN
3576 ic = 1
3577 ELSE
3578 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
3579 cpassert(ic > 0)
3580 END IF
3581 nkp_cellind = 0
3582 DO i = 1, 3
3583 IF (nkp_grid(i) > 1) THEN
3584 nkp_cellind(i) = modulo(cellind(i), nkp_grid(i))
3585 IF (2*nkp_cellind(i) > nkp_grid(i)) nkp_cellind(i) = nkp_cellind(i) - nkp_grid(i)
3586 END IF
3587 END DO
3588 sampled_axes = 0
3589 IF (nkp_cellind(1) /= 0) sampled_axes = sampled_axes + 1
3590 IF (nkp_cellind(2) /= 0) sampled_axes = sampled_axes + 1
3591 IF (nkp_cellind(3) /= 0) sampled_axes = sampled_axes + 1
3592 sampled_even_axes = 0
3593 sampled_gamma_axes = 0
3594 DO i = 1, 3
3595 IF (nkp_grid(i) > 1 .AND. modulo(nkp_grid(i), 2) == 0 .AND. &
3596 abs(2*nkp_cellind(i)) == nkp_grid(i)) THEN
3597 sampled_even_axes = sampled_even_axes + 1
3598 IF (mesh_has_gamma(i)) sampled_gamma_axes = sampled_gamma_axes + 1
3599 END IF
3600 END DO
3601 sampled_image_pair = sampled_axes > 0 .AND. sampled_even_axes > 0
3602 gamma_sampled_image_pair = sampled_image_pair .AND. sampled_gamma_axes == sampled_even_axes
3603 mp_pair_scale = 1.0_dp
3604 IF (icol == jrow .AND. sampled_image_pair .AND. has_multipole_response .AND. &
3605 .NOT. gamma_sampled_image_pair) mp_pair_scale = -1.0_dp
3606
3607 NULLIFY (pblock, pblock_beta)
3608 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
3609 row=jrow, col=icol, block=pblock, found=found)
3610 IF (.NOT. found) cpabort("pblock not found")
3611 IF (nspin > 1) THEN
3612 CALL dbcsr_get_block_p(matrix=matrix_p(2, ic)%matrix, &
3613 row=jrow, col=icol, block=pblock_beta, found=found)
3614 IF (.NOT. found) cpabort("pblock beta not found")
3615 END IF
3616 IF (use_matrix_scc_stress) THEN
3617 ALLOCATE (scc_strain_hint(SIZE(pblock, 2), SIZE(pblock, 1), 3, 3, nspin))
3618 scc_strain_hint = 0.0_dp
3619 DO ispin = 1, nspin
3620 DO jdim = 1, 3
3621 DO idim = 1, 3
3622 NULLIFY (scc_strain_blocks(idim, jdim, ispin)%block)
3623 i = idim + 3*(jdim - 1) + 9*(ispin - 1)
3624 CALL dbcsr_get_block_p(matrix=matrix_dh_scc(i, ic)%matrix, &
3625 row=jrow, col=icol, &
3626 block=scc_strain_blocks(idim, jdim, ispin)%block, found=found)
3627 IF (.NOT. found) cpabort("SCC strain derivative block not found")
3628 END DO
3629 END DO
3630 END DO
3631 END IF
3632 i_a_shift = tb%pot%vat(icol, 1)
3633 j_a_shift = tb%pot%vat(jrow, 1)
3634 i_a_shift_mag = 0.0_dp
3635 j_a_shift_mag = 0.0_dp
3636 IF (SIZE(tb%pot%vat, 2) > 1) THEN
3637 i_a_shift_mag = tb%pot%vat(icol, 2)
3638 j_a_shift_mag = tb%pot%vat(jrow, 2)
3639 END IF
3640 idip(:) = tb%pot%vdp(:, icol, 1)
3641 jdip(:) = tb%pot%vdp(:, jrow, 1)
3642 idip_mag(:) = 0.0_dp
3643 jdip_mag(:) = 0.0_dp
3644 IF (SIZE(tb%pot%vdp, 3) > 1) THEN
3645 idip_mag(:) = tb%pot%vdp(:, icol, 2)
3646 jdip_mag(:) = tb%pot%vdp(:, jrow, 2)
3647 END IF
3648 iquad(:) = tb%pot%vqp(:, icol, 1)
3649 jquad(:) = tb%pot%vqp(:, jrow, 1)
3650 iquad_mag(:) = 0.0_dp
3651 jquad_mag(:) = 0.0_dp
3652 IF (SIZE(tb%pot%vqp, 3) > 1) THEN
3653 iquad_mag(:) = tb%pot%vqp(:, icol, 2)
3654 jquad_mag(:) = tb%pot%vqp(:, jrow, 2)
3655 END IF
3656 ni = tb%calc%bas%ish_at(icol)
3657 DO iset = 1, nseti
3658 ishift = i_a_shift + tb%pot%vsh(ni + iset, 1)
3659 ishift_mag = 0.0_dp
3660 IF (SIZE(tb%pot%vsh, 2) > 1) ishift_mag = i_a_shift_mag + tb%pot%vsh(ni + iset, 2)
3661 nj = tb%calc%bas%ish_at(jrow)
3662 DO jset = 1, nsetj
3663 jshift = j_a_shift + tb%pot%vsh(nj + jset, 1)
3664 jshift_mag = 0.0_dp
3665 IF (SIZE(tb%pot%vsh, 2) > 1) jshift_mag = j_a_shift_mag + tb%pot%vsh(nj + jset, 2)
3666
3667 !get integrals and derivatives
3668 CALL multipole_grad_cgto(tb%calc%bas%cgto(iset, ityp), tb%calc%bas%cgto(jset, jtyp), &
3669 & r2, rij, tb%calc%bas%intcut, t_ov, t_dip, t_quad, t_d_ov, t_i_dip, t_i_quad, &
3670 & t_j_dip, t_j_quad)
3671
3672 dgrad = 0.0_dp
3673 DO inda = 1, nsgfa(iset)
3674 ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
3675 DO indb = 1, nsgfb(jset)
3676 ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
3677
3678 ij = inda + nsgfa(iset)*(indb - 1)
3679
3680 pij_charge = pblock(ib, ia)
3681 pij_magnet = 0.0_dp
3682 IF (nspin > 1) THEN
3683 pij_charge = pij_charge + pblock_beta(ib, ia)
3684 pij_magnet = pblock(ib, ia) - pblock_beta(ib, ia)
3685 END IF
3686 mpgrad_charge = matmul(t_i_dip(:, :, ij), idip) &
3687 + matmul(t_j_dip(:, :, ij), jdip) &
3688 + matmul(t_i_quad(:, :, ij), iquad) &
3689 + matmul(t_j_quad(:, :, ij), jquad)
3690 mpgrad_magnet = matmul(t_i_dip(:, :, ij), idip_mag) &
3691 + matmul(t_j_dip(:, :, ij), jdip_mag) &
3692 + matmul(t_i_quad(:, :, ij), iquad_mag) &
3693 + matmul(t_j_quad(:, :, ij), jquad_mag)
3694 dhgrad_charge = -(ishift + jshift)*t_d_ov(:, ij) - mp_pair_scale*mpgrad_charge
3695 dhgrad_magnet = -(ishift_mag + jshift_mag)*t_d_ov(:, ij) - mp_pair_scale*mpgrad_magnet
3696 dgrad(:) = dgrad(:) - &
3697 ((ishift + jshift)*pij_charge + &
3698 (ishift_mag + jshift_mag)*pij_magnet)*t_d_ov(:, ij) - &
3699 mp_pair_scale*(pij_charge*mpgrad_charge + pij_magnet*mpgrad_magnet)
3700
3701 IF (ALLOCATED(scc_strain_hint)) THEN
3702 DO jdim = 1, 3
3703 DO idim = 1, 3
3704 scc_strain_hint(ia, ib, idim, jdim, 1) = &
3705 scc_strain_hint(ia, ib, idim, jdim, 1) &
3706 + (dhgrad_charge(idim) + merge(dhgrad_magnet(idim), 0.0_dp, nspin > 1))*rij(jdim)
3707 IF (nspin > 1) THEN
3708 scc_strain_hint(ia, ib, idim, jdim, 2) = &
3709 scc_strain_hint(ia, ib, idim, jdim, 2) &
3710 + (dhgrad_charge(idim) - dhgrad_magnet(idim))*rij(jdim)
3711 END IF
3712 END DO
3713 END DO
3714 END IF
3715
3716 END DO
3717 END DO
3718 tb%grad(:, icol) = tb%grad(:, icol) - dgrad
3719 tb%grad(:, jrow) = tb%grad(:, jrow) + dgrad
3720 IF (tb%use_virial .AND. .NOT. use_matrix_scc_stress) THEN
3721 IF (icol == jrow) THEN
3722 DO ia = 1, 3
3723 DO ib = 1, 3
3724 IF (sampled_image_pair .AND. .NOT. gamma_sampled_image_pair) THEN
3725 hsigma(ia, ib) = hsigma(ia, ib) - 0.25_dp* &
3726 (rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
3727 ELSE
3728 hsigma(ia, ib) = hsigma(ia, ib) + 0.25_dp* &
3729 (rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
3730 END IF
3731 END DO
3732 END DO
3733 ELSE
3734 DO ia = 1, 3
3735 DO ib = 1, 3
3736 hsigma(ia, ib) = hsigma(ia, ib) + 0.50_dp*(rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
3737 END DO
3738 END DO
3739 END IF
3740 END IF
3741 END DO
3742 END DO
3743 IF (ALLOCATED(scc_strain_hint)) THEN
3744 DO ispin = 1, nspin
3745 DO jdim = 1, 3
3746 DO idim = 1, 3
3747 IF (icol <= jrow) THEN
3748 scc_strain_blocks(idim, jdim, ispin)%block(:, :) = &
3749 scc_strain_blocks(idim, jdim, ispin)%block(:, :) &
3750 + scc_strain_hint(:, :, idim, jdim, ispin)
3751 ELSE
3752 scc_strain_blocks(idim, jdim, ispin)%block(:, :) = &
3753 scc_strain_blocks(idim, jdim, ispin)%block(:, :) &
3754 + transpose(scc_strain_hint(:, :, idim, jdim, ispin))
3755 END IF
3756 END DO
3757 END DO
3758 END DO
3759 DEALLOCATE (scc_strain_hint)
3760 END IF
3761 END DO
3762 CALL neighbor_list_iterator_release(nl_iterator)
3763
3764 IF (use_matrix_scc_stress) THEN
3765 hsigma = 0.0_dp
3766 DO img = 1, nimg
3767 DO ispin = 1, nspin
3768 DO jdim = 1, 3
3769 DO idim = 1, 3
3770 i = idim + 3*(jdim - 1) + 9*(ispin - 1)
3771 CALL dbcsr_finalize(matrix_dh_scc(i, img)%matrix)
3772 CALL dbcsr_dot(matrix_dh_scc(i, img)%matrix, matrix_p(ispin, img)%matrix, native_dot_tmp)
3773 ! Symmetric DBCSR storage represents both orientations of each assembled pair.
3774 hsigma(idim, jdim) = hsigma(idim, jdim) + 0.5_dp*native_dot_tmp
3775 END DO
3776 END DO
3777 END DO
3778 END DO
3779 CALL dbcsr_deallocate_matrix_set(matrix_dh_scc)
3780 ELSE
3781 CALL para_env%sum(hsigma)
3782 END IF
3783 CALL para_env%sum(tb%grad)
3784 CALL tb_grad2force(qs_env, tb, para_env, 4)
3785
3786 IF (.NOT. use_matrix_scc_stress) tb%sigma = tb%sigma + hsigma
3787
3788 DEALLOCATE (basis_set_list)
3789 DEALLOCATE (t_ov, t_d_ov)
3790 DEALLOCATE (t_dip, t_i_dip, t_j_dip)
3791 DEALLOCATE (t_quad, t_i_quad, t_j_quad)
3792 DEALLOCATE (idip, jdip, idip_mag, jdip_mag, iquad, jquad, iquad_mag, jquad_mag)
3793
3794 IF (tb%use_virial) THEN
3795 CALL tb_add_stress(qs_env, tb, para_env)
3796 IF (use_matrix_scc_stress) THEN
3797 CALL get_qs_env(qs_env=qs_env, virial=virial)
3798 virial%pv_virial = virial%pv_virial - hsigma/para_env%num_pe
3799 END IF
3800 END IF
3801
3802#else
3803 mark_used(qs_env)
3804 mark_used(use_rho)
3805 mark_used(nimg)
3806 cpabort("Built without TBLITE")
3807#endif
3808
3809 END SUBROUTINE tb_derive_dh_off
3810
3811! **************************************************************************************************
3812!> \brief Run native tblite CLI and compare against CP2K/tblite.
3813!> \param qs_env ...
3814! **************************************************************************************************
3815 SUBROUTINE tb_reference_cli_compare(qs_env)
3816
3817 TYPE(qs_environment_type), POINTER :: qs_env
3818
3819 CHARACTER(LEN=*), PARAMETER :: routinen = 'tb_reference_cli_compare'
3820
3821 CHARACTER(LEN=16) :: solvation_model_name
3822 CHARACTER(LEN=32) :: acc_str, charge_str, efield_x_str, efield_y_str, efield_z_str, &
3823 etemp_guess_val_str, etemp_str, iter_str, spin_str, spinpol_str
3824 CHARACTER(LEN=4*default_path_length+16) :: efield_str, etemp_guess_str, param_str, &
3825 post_processing_output_str, post_processing_str, restart_str, solvation_str, verbosity_str
3826 CHARACTER(LEN=8) :: guess, method, solver
3827 CHARACTER(LEN=8*default_path_length) :: command
3828 CHARACTER(LEN=default_path_length) :: file_base, gen_file, grad_file, &
3829 json_file, log_file, &
3830 post_processing_output_file
3831 INTEGER :: cmdstat, exitstat, handle, iounit, &
3832 n_periodic, natom, nkp, &
3833 reference_iterations, spin
3834 INTEGER, DIMENSION(3) :: periodic
3835 LOGICAL :: do_kpoints, have_energy, have_gradient, &
3836 have_virial, too_large, &
3837 unsupported_kpoints
3838 REAL(kind=dp) :: cli_energy, cp_energy, ediff, etemp, &
3839 etemp_guess, fmax, fsum, vmax, vsum
3840 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cli_gradient, cli_virial, cp_gradient
3841 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
3842 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3843 TYPE(cell_type), POINTER :: cell
3844 TYPE(cp_logger_type), POINTER :: logger
3845 TYPE(dft_control_type), POINTER :: dft_control
3846 TYPE(kpoint_type), POINTER :: kpoints
3847 TYPE(mp_para_env_type), POINTER :: para_env
3848 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3849 TYPE(qs_energy_type), POINTER :: energy
3850 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3851 TYPE(scf_control_type), POINTER :: scf_control
3852 TYPE(virial_type), POINTER :: virial
3853 TYPE(xtb_reference_cli_type) :: ref
3854
3855 CALL timeset(routinen, handle)
3856
3857 NULLIFY (atomic_kind_set, cell, dft_control, energy, force, kpoints, logger, para_env, particle_set, &
3858 scf_control, virial, xkp)
3859 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell, &
3860 dft_control=dft_control, energy=energy, force=force, &
3861 para_env=para_env, particle_set=particle_set, scf_control=scf_control, &
3862 virial=virial, do_kpoints=do_kpoints, kpoints=kpoints)
3863
3864 ref = dft_control%qs_control%xtb_control%reference_cli
3865 IF (.NOT. ref%enabled) THEN
3866 CALL timestop(handle)
3867 RETURN
3868 END IF
3869 IF (.NOT. para_env%is_source()) THEN
3870 CALL timestop(handle)
3871 RETURN
3872 END IF
3873
3874 logger => cp_get_default_logger()
3875 iounit = cp_logger_get_default_io_unit(logger)
3876 verbosity_str = ""
3877 IF (logger%iter_info%print_level == silent_print_level) THEN
3878 verbosity_str = " --silent"
3879 ELSE IF (logger%iter_info%print_level == high_print_level .OR. &
3880 logger%iter_info%print_level == debug_print_level) THEN
3881 verbosity_str = " --verbose"
3882 END IF
3883 IF (ref%solvation_active) THEN
3884 periodic = 0
3885 IF (ASSOCIATED(cell)) CALL get_cell(cell=cell, periodic=periodic)
3886 n_periodic = count(periodic == 1)
3887 IF (n_periodic == 3) THEN
3888 WRITE (unit=iounit, fmt="(/,T2,A)") &
3889 "tblite reference CLI implicit solvation is not supported for PERIODIC XYZ."
3890 WRITE (unit=iounit, fmt="(T2,A)") &
3891 "Use PERIODIC NONE for molecular solvation diagnostics, or remove IMPLICIT_SOLVATION."
3892 cpabort("REFERENCE_CLI implicit solvation is incompatible with PERIODIC XYZ")
3893 ELSE IF (n_periodic > 0) THEN
3894 WRITE (unit=iounit, fmt="(/,T2,A)") &
3895 "WARNING: tblite reference CLI implicit solvation with finite periodicity is diagnostic only."
3896 WRITE (unit=iounit, fmt="(T2,A,I0,A)") &
3897 "The generated native tblite reference geometry has ", n_periodic, &
3898 " periodic direction(s); continuum-solvation conventions are primarily molecular."
3899 END IF
3900 END IF
3901 unsupported_kpoints = .false.
3902 IF (do_kpoints .AND. ASSOCIATED(kpoints)) THEN
3903 nkp = 0
3904 CALL get_kpoint_info(kpoint=kpoints, nkp=nkp, xkp=xkp)
3905 unsupported_kpoints = nkp > 1
3906 IF (nkp == 1 .AND. ASSOCIATED(xkp)) unsupported_kpoints = any(abs(xkp(:, 1)) > 1.0e-12_dp)
3907 END IF
3908 IF (unsupported_kpoints) THEN
3909 WRITE (unit=iounit, fmt="(/,T2,A)") &
3910 "tblite reference CLI check skipped: CP2K KPOINTS are active."
3911 WRITE (unit=iounit, fmt="(T2,A)") &
3912 "The native tblite CLI reference path does not reproduce CP2K multi-k-point sampling."
3913 IF (ref%stop_on_error) cpabort("tblite reference CLI cannot check CP2K k-point calculations")
3914 CALL timestop(handle)
3915 RETURN
3916 END IF
3917 IF (dft_control%qs_control%xtb_control%tblite_scc_mixer == tblite_scc_mixer_cp2k) THEN
3918 WRITE (unit=iounit, fmt="(/,T2,A)") &
3919 "WARNING: tblite reference CLI cannot reproduce XTB/SCC_MIXER CP2K."
3920 WRITE (unit=iounit, fmt="(T2,A)") &
3921 "The external native tblite run uses tblite's own SCC mixer; only the converged result is compared."
3922 IF (ref%stop_on_error) cpabort("tblite reference CLI cannot reproduce SCC_MIXER CP2K")
3923 END IF
3924 natom = SIZE(particle_set)
3925 method = tb_reference_method_name(dft_control%qs_control%xtb_control%tblite_method)
3926 guess = tb_reference_guess_name(ref%guess)
3927 solver = tb_reference_solver_name(dft_control%qs_control%xtb_control%tblite_mixer_solver)
3928 file_base = tb_join_path(ref%work_directory, ref%prefix)
3929 gen_file = trim(file_base)//".gen"
3930 grad_file = trim(file_base)//".grad"
3931 json_file = trim(file_base)//".json"
3932 log_file = trim(file_base)//".log"
3933 post_processing_output_file = ""
3934 IF (len_trim(ref%grad_file) > 0) grad_file = ref%grad_file
3935 IF (len_trim(ref%json_file) > 0) json_file = ref%json_file
3936 IF (len_trim(ref%post_processing_output_file) > 0) THEN
3937 post_processing_output_file = ref%post_processing_output_file
3938 END IF
3939
3940 WRITE (charge_str, "(I0)") dft_control%charge
3941 spin = max(0, dft_control%multiplicity - 1)
3942 WRITE (spin_str, "(I0)") spin
3943 WRITE (acc_str, "(ES16.8)") dft_control%qs_control%xtb_control%tblite_accuracy
3944 reference_iterations = dft_control%qs_control%xtb_control%tblite_mixer_iterations
3945 WRITE (iter_str, "(I0)") reference_iterations
3946 efield_str = ""
3947 IF (ref%efield_active) THEN
3948 WRITE (efield_x_str, "(ES16.8)") ref%efield(1)
3949 WRITE (efield_y_str, "(ES16.8)") ref%efield(2)
3950 WRITE (efield_z_str, "(ES16.8)") ref%efield(3)
3951 efield_str = " --efield "//trim(adjustl(efield_x_str))//","// &
3952 trim(adjustl(efield_y_str))//","//trim(adjustl(efield_z_str))
3953 END IF
3954 solvation_str = ""
3955 solvation_model_name = ""
3956 IF (ref%solvation_active) THEN
3957 SELECT CASE (ref%solvation_model)
3958 CASE (tblite_cli_solvation_alpb)
3959 solvation_model_name = "ALPB"
3960 solvation_str = " --alpb "
3961 CASE (tblite_cli_solvation_gbsa)
3962 solvation_model_name = "GBSA"
3963 solvation_str = " --gbsa "
3964 CASE (tblite_cli_solvation_gbe)
3965 solvation_model_name = "GBE"
3966 solvation_str = " --gbe "
3967 CASE (tblite_cli_solvation_gb)
3968 solvation_model_name = "GB"
3969 solvation_str = " --gb "
3970 CASE (tblite_cli_solvation_cpcm)
3971 solvation_model_name = "CPCM"
3972 solvation_str = " --cpcm "
3973 CASE DEFAULT
3974 cpabort("Unknown tblite reference CLI implicit-solvation model")
3975 END SELECT
3976 solvation_str = trim(solvation_str)//" "//trim(tb_shell_quote(ref%solvation_solvent))
3977 SELECT CASE (ref%solvation_born_kernel)
3978 CASE (tblite_cli_born_kernel_auto)
3979 CASE (tblite_cli_born_kernel_p16)
3980 solvation_str = trim(solvation_str)//" --born-kernel p16"
3981 CASE (tblite_cli_born_kernel_still)
3982 solvation_str = trim(solvation_str)//" --born-kernel still"
3983 CASE DEFAULT
3984 cpabort("Unknown tblite reference CLI Born kernel")
3985 END SELECT
3986 SELECT CASE (ref%solvation_state)
3987 CASE (tblite_cli_solution_state_gsolv)
3988 CASE (tblite_cli_solution_state_bar1mol)
3989 solvation_str = trim(solvation_str)//" --solv-state bar1mol"
3990 CASE (tblite_cli_solution_state_reference)
3991 solvation_str = trim(solvation_str)//" --solv-state reference"
3992 CASE DEFAULT
3993 cpabort("Unknown tblite reference CLI solution state")
3994 END SELECT
3995 END IF
3996 IF (dft_control%qs_control%xtb_control%tblite_mixer_memory /= reference_iterations) THEN
3997 WRITE (unit=iounit, fmt="(/,T2,A)") &
3998 "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/MEMORY."
3999 WRITE (unit=iounit, fmt="(T2,A,I0,A,I0,A)") &
4000 "The native reference run uses tblite's internal mixer memory tied to --iterations (", &
4001 reference_iterations, "), while CP2K uses MEMORY ", &
4002 dft_control%qs_control%xtb_control%tblite_mixer_memory, "."
4003 IF (ref%stop_on_error) cpabort("tblite reference CLI cannot reproduce TBLITE_MIXER/MEMORY")
4004 END IF
4005 IF (abs(dft_control%qs_control%xtb_control%tblite_mixer_damping - &
4006 tblite_mixer_damping_default) > 10.0_dp*epsilon(1.0_dp)) THEN
4007 WRITE (unit=iounit, fmt="(/,T2,A)") &
4008 "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/DAMPING."
4009 WRITE (unit=iounit, fmt="(T2,A,F8.4,A,F8.4,A)") &
4010 "The native reference run uses tblite's library default ", tblite_mixer_damping_default, &
4011 ", while CP2K uses DAMPING ", dft_control%qs_control%xtb_control%tblite_mixer_damping, "."
4012 IF (ref%stop_on_error) cpabort("tblite reference CLI cannot reproduce TBLITE_MIXER/DAMPING")
4013 END IF
4014 IF (abs(dft_control%qs_control%xtb_control%tblite_mixer_omega0 - &
4015 tblite_mixer_omega0_default) > 10.0_dp*epsilon(1.0_dp)) THEN
4016 WRITE (unit=iounit, fmt="(/,T2,A)") &
4017 "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/OMEGA0."
4018 WRITE (unit=iounit, fmt="(T2,A,ES12.4,A,ES12.4,A)") &
4019 "The native reference run uses tblite's library default ", tblite_mixer_omega0_default, &
4020 ", while CP2K uses OMEGA0 ", dft_control%qs_control%xtb_control%tblite_mixer_omega0, "."
4021 IF (ref%stop_on_error) cpabort("tblite reference CLI cannot reproduce TBLITE_MIXER/OMEGA0")
4022 END IF
4023 IF (abs(dft_control%qs_control%xtb_control%tblite_mixer_min_weight - &
4024 tblite_mixer_min_weight_default) > 10.0_dp*epsilon(1.0_dp)) THEN
4025 WRITE (unit=iounit, fmt="(/,T2,A)") &
4026 "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/MIN_WEIGHT."
4027 WRITE (unit=iounit, fmt="(T2,A,ES12.4,A,ES12.4,A)") &
4028 "The native reference run uses tblite's library default ", tblite_mixer_min_weight_default, &
4029 ", while CP2K uses MIN_WEIGHT ", dft_control%qs_control%xtb_control%tblite_mixer_min_weight, "."
4030 IF (ref%stop_on_error) THEN
4031 cpabort("tblite reference CLI cannot reproduce TBLITE_MIXER/MIN_WEIGHT")
4032 END IF
4033 END IF
4034 IF (abs(dft_control%qs_control%xtb_control%tblite_mixer_max_weight - &
4035 tblite_mixer_max_weight_default) > 10.0_dp*epsilon(1.0_dp)) THEN
4036 WRITE (unit=iounit, fmt="(/,T2,A)") &
4037 "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/MAX_WEIGHT."
4038 WRITE (unit=iounit, fmt="(T2,A,ES12.4,A,ES12.4,A)") &
4039 "The native reference run uses tblite's library default ", tblite_mixer_max_weight_default, &
4040 ", while CP2K uses MAX_WEIGHT ", dft_control%qs_control%xtb_control%tblite_mixer_max_weight, "."
4041 IF (ref%stop_on_error) THEN
4042 cpabort("tblite reference CLI cannot reproduce TBLITE_MIXER/MAX_WEIGHT")
4043 END IF
4044 END IF
4045 IF (abs(dft_control%qs_control%xtb_control%tblite_mixer_weight_factor - &
4046 tblite_mixer_weight_factor_default) > 10.0_dp*epsilon(1.0_dp)) THEN
4047 WRITE (unit=iounit, fmt="(/,T2,A)") &
4048 "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/WEIGHT_FACTOR."
4049 WRITE (unit=iounit, fmt="(T2,A,ES12.4,A,ES12.4,A)") &
4050 "The native reference run uses tblite's library default ", tblite_mixer_weight_factor_default, &
4051 ", while CP2K uses WEIGHT_FACTOR ", &
4052 dft_control%qs_control%xtb_control%tblite_mixer_weight_factor, "."
4053 IF (ref%stop_on_error) THEN
4054 cpabort("tblite reference CLI cannot reproduce TBLITE_MIXER/WEIGHT_FACTOR")
4055 END IF
4056 END IF
4057 etemp = 300.0_dp
4058 IF (ASSOCIATED(scf_control)) THEN
4059 IF (ASSOCIATED(scf_control%smear)) THEN
4060 IF (scf_control%smear%do_smear) THEN
4061 etemp = cp_unit_from_cp2k(scf_control%smear%electronic_temperature, "K")
4062 IF (scf_control%smear%method /= smear_fermi_dirac) THEN
4063 WRITE (unit=iounit, fmt="(/,T2,A,A,A)") &
4064 "WARNING: tblite reference CLI cannot reproduce CP2K smearing method ", &
4065 trim(tb_reference_smear_method_name(scf_control%smear%method)), "."
4066 WRITE (unit=iounit, fmt="(T2,A,F12.3,A)") &
4067 "The native reference run uses Fermi-Dirac electronic temperature ", etemp, " K instead."
4068 END IF
4069 END IF
4070 END IF
4071 END IF
4072 WRITE (etemp_str, "(ES16.8)") etemp
4073 etemp_guess = 0.0_dp
4074 etemp_guess_str = ""
4075 IF (ref%electronic_temperature_guess > 0.0_dp) THEN
4076 etemp_guess = cp_unit_from_cp2k(ref%electronic_temperature_guess, "K")
4077 WRITE (etemp_guess_val_str, "(ES16.8)") etemp_guess
4078 etemp_guess_str = " --etemp-guess "//trim(adjustl(etemp_guess_val_str))
4079 END IF
4080 param_str = ""
4081 IF (len_trim(dft_control%qs_control%xtb_control%tblite_param_file) > 0) THEN
4082 param_str = " --param "//trim(tb_shell_quote(dft_control%qs_control%xtb_control%tblite_param_file))
4083 END IF
4084 spinpol_str = ""
4085 IF (dft_control%lsd) spinpol_str = " --spin-polarized"
4086 post_processing_str = ""
4087 IF (len_trim(ref%post_processing) > 0) THEN
4088 post_processing_str = " --post-processing "//trim(tb_shell_quote(ref%post_processing))
4089 END IF
4090 post_processing_output_str = ""
4091 IF (len_trim(post_processing_output_file) > 0) THEN
4092 WRITE (unit=iounit, fmt="(/,T2,A)") &
4093 "WARNING: tblite reference CLI POST_PROCESSING_OUTPUT was requested explicitly."
4094 WRITE (unit=iounit, fmt="(T2,A)") &
4095 "Some tblite 0.5.0 command-line builds document --post-processing-output but do not parse it."
4096 post_processing_output_str = " --post-processing-output "// &
4097 trim(tb_shell_quote(post_processing_output_file))
4098 END IF
4099 restart_str = " --no-restart"
4100 IF (len_trim(ref%restart_file) > 0) THEN
4101 restart_str = " --restart "//trim(tb_shell_quote(ref%restart_file))
4102 END IF
4103
4104 CALL tb_write_reference_gen(qs_env, trim(gen_file))
4105
4106 command = trim(tb_shell_quote(ref%program_name))//" run --method "//trim(method)// &
4107 trim(param_str)// &
4108 trim(spinpol_str)// &
4109 " --charge "//trim(adjustl(charge_str))// &
4110 " --spin "//trim(adjustl(spin_str))// &
4111 " --acc "//trim(adjustl(acc_str))// &
4112 " --guess "//trim(guess)// &
4113 " --solver "//trim(solver)// &
4114 " --iterations "//trim(adjustl(iter_str))// &
4115 " --etemp "//trim(adjustl(etemp_str))// &
4116 trim(etemp_guess_str)// &
4117 trim(efield_str)// &
4118 trim(solvation_str)// &
4119 trim(post_processing_str)// &
4120 trim(post_processing_output_str)// &
4121 trim(restart_str)// &
4122 trim(verbosity_str)//" --input "//trim(tb_shell_quote(ref%input_format))// &
4123 " --grad "//trim(tb_shell_quote(grad_file))// &
4124 " --json "//trim(tb_shell_quote(json_file))//" "//trim(tb_shell_quote(gen_file))// &
4125 " > "//trim(tb_shell_quote(log_file))//" 2>&1"
4126
4127 cmdstat = 0
4128 exitstat = 0
4129 CALL execute_command_line(trim(command), exitstat=exitstat, cmdstat=cmdstat)
4130 IF (cmdstat /= 0 .OR. exitstat /= 0) THEN
4131 WRITE (unit=iounit, fmt="(/,T2,A)") "tblite reference CLI check failed to run."
4132 WRITE (unit=iounit, fmt="(T2,A,A)") "Command: ", trim(command)
4133 WRITE (unit=iounit, fmt="(T2,A,I0,T32,A,I0)") "cmdstat:", cmdstat, "exitstat:", exitstat
4134 IF (ref%stop_on_error) cpabort("tblite reference CLI command failed")
4135 CALL tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file, post_processing_output_file)
4136 CALL timestop(handle)
4137 RETURN
4138 END IF
4139
4140 ALLOCATE (cli_gradient(3, natom), cli_virial(3, 3))
4141 CALL tb_read_reference_grad(trim(grad_file), natom, cli_energy, cli_gradient, cli_virial, &
4142 have_energy, have_gradient, have_virial)
4143
4144 WRITE (unit=iounit, fmt="(/,T2,A)") "tblite reference CLI check"
4145 WRITE (unit=iounit, fmt="(T2,A,A)") "Executable: ", trim(ref%program_name)
4146 WRITE (unit=iounit, fmt="(T2,A,A)") "Method: ", trim(method)
4147 WRITE (unit=iounit, fmt="(T2,A,A)") "Guess: ", trim(guess)
4148 WRITE (unit=iounit, fmt="(T2,A,A)") "Solver: ", trim(solver)
4149 WRITE (unit=iounit, fmt="(T2,A,L1)") "Spin-pol.: ", dft_control%lsd
4150 IF (ref%efield_active) THEN
4151 WRITE (unit=iounit, fmt="(T2,A,3ES16.8,A)") "Efield: ", ref%efield, " V/Angstrom"
4152 END IF
4153 IF (ref%solvation_active) THEN
4154 WRITE (unit=iounit, fmt="(T2,A,A,1X,A)") "Solvation: ", trim(solvation_model_name), &
4155 trim(ref%solvation_solvent)
4156 END IF
4157 IF (len_trim(ref%post_processing) > 0) THEN
4158 WRITE (unit=iounit, fmt="(T2,A,A)") "Post proc.: ", trim(ref%post_processing)
4159 END IF
4160 IF (len_trim(post_processing_output_file) > 0) THEN
4161 WRITE (unit=iounit, fmt="(T2,A,A)") "PP output: ", trim(post_processing_output_file)
4162 END IF
4163 IF (ref%electronic_temperature_guess > 0.0_dp) THEN
4164 WRITE (unit=iounit, fmt="(T2,A,F12.3,A)") "Guess etemp:", etemp_guess, " K"
4165 END IF
4166 WRITE (unit=iounit, fmt="(T2,A,A)") "Grad file: ", trim(grad_file)
4167 WRITE (unit=iounit, fmt="(T2,A,A)") "JSON file: ", trim(json_file)
4168 WRITE (unit=iounit, fmt="(T2,A,A)") "Log file: ", trim(log_file)
4169
4170 too_large = .false.
4171 IF (ref%check_energy) THEN
4172 IF (have_energy) THEN
4173 cp_energy = energy%total
4174 ediff = abs(cp_energy - cli_energy)
4175 WRITE (unit=iounit, fmt="(T2,A,3ES22.12)") &
4176 "Energy CP2K/CLI/absdiff:", cp_energy, cli_energy, ediff
4177 too_large = too_large .OR. ediff > ref%error_limit
4178 ELSE
4179 WRITE (unit=iounit, fmt="(T2,A)") "Energy check skipped: no CLI energy found."
4180 END IF
4181 END IF
4182
4183 IF (ref%check_forces) THEN
4184 IF (have_gradient .AND. ASSOCIATED(force)) THEN
4185 ALLOCATE (cp_gradient(3, natom))
4186 CALL total_qs_force(cp_gradient, force, atomic_kind_set)
4187 fsum = sum(abs(cp_gradient - cli_gradient))
4188 fmax = maxval(abs(cp_gradient - cli_gradient))
4189 WRITE (unit=iounit, fmt="(T2,A,2ES22.12)") "Gradient diff sum/max:", fsum, fmax
4190 too_large = too_large .OR. fmax > ref%error_limit
4191 DEALLOCATE (cp_gradient)
4192 ELSE
4193 WRITE (unit=iounit, fmt="(T2,A)") "Gradient check skipped: no CLI gradient or CP2K force found."
4194 END IF
4195 END IF
4196
4197 IF (ref%check_virial) THEN
4198 IF (have_virial .AND. ASSOCIATED(virial)) THEN
4199 ! Native tblite prints the positive cell derivative; CP2K stores the PV virial
4200 ! with the opposite sign.
4201 vsum = sum(abs(-virial%pv_virial - cli_virial))
4202 vmax = maxval(abs(-virial%pv_virial - cli_virial))
4203 WRITE (unit=iounit, fmt="(T2,A,2ES22.12)") "Virial diff sum/max:", vsum, vmax
4204 too_large = too_large .OR. vmax > ref%error_limit
4205 ELSE
4206 WRITE (unit=iounit, fmt="(T2,A)") "Virial check skipped: no CLI virial or CP2K virial found."
4207 END IF
4208 END IF
4209
4210 IF (too_large) THEN
4211 WRITE (unit=iounit, fmt="(T2,A,ES12.4)") &
4212 "tblite reference CLI deviation exceeded ERROR_LIMIT = ", ref%error_limit
4213 IF (ref%stop_on_error) cpabort("tblite reference CLI deviation exceeded ERROR_LIMIT")
4214 END IF
4215
4216 CALL tb_reference_cli_aux_commands(ref, dft_control, gen_file, file_base, verbosity_str, iounit)
4217
4218 CALL tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file, post_processing_output_file)
4219 DEALLOCATE (cli_gradient, cli_virial)
4220
4221 CALL timestop(handle)
4222
4223 END SUBROUTINE tb_reference_cli_compare
4224
4225! **************************************************************************************************
4226!> \brief Map CP2K tblite method id to native tblite CLI method name.
4227!> \param method_id ...
4228!> \return ...
4229! **************************************************************************************************
4230 FUNCTION tb_reference_method_name(method_id) RESULT(method)
4231 INTEGER, INTENT(IN) :: method_id
4232 CHARACTER(LEN=8) :: method
4233
4234 SELECT CASE (method_id)
4235 CASE (gfn1xtb)
4236 method = "gfn1"
4237 CASE (gfn2xtb)
4238 method = "gfn2"
4239 CASE (ipea1xtb)
4240 method = "ipea1"
4241 CASE DEFAULT
4242 cpabort("Unknown tblite reference CLI method")
4243 END SELECT
4244
4245 END FUNCTION tb_reference_method_name
4246
4247! **************************************************************************************************
4248!> \brief Map CP2K tblite reference CLI guess id to native tblite CLI guess name.
4249!> \param guess_id ...
4250!> \return ...
4251! **************************************************************************************************
4252 FUNCTION tb_reference_guess_name(guess_id) RESULT(guess)
4253 INTEGER, INTENT(IN) :: guess_id
4254 CHARACTER(LEN=8) :: guess
4255
4256 SELECT CASE (guess_id)
4257 CASE (tblite_guess_sad)
4258 guess = "sad"
4259 CASE (tblite_guess_eeq)
4260 guess = "eeq"
4261 CASE (tblite_guess_ceh)
4262 guess = "ceh"
4263 CASE DEFAULT
4264 cpabort("Unknown tblite reference CLI guess")
4265 END SELECT
4266
4267 END FUNCTION tb_reference_guess_name
4268
4269! **************************************************************************************************
4270!> \brief Map CP2K tblite solver id to native tblite CLI solver name.
4271!> \param solver_id ...
4272!> \return ...
4273! **************************************************************************************************
4274 FUNCTION tb_reference_solver_name(solver_id) RESULT(solver)
4275 INTEGER, INTENT(IN) :: solver_id
4276 CHARACTER(LEN=8) :: solver
4277
4278 SELECT CASE (solver_id)
4279 CASE (tblite_solver_gvd)
4280 solver = "gvd"
4281 CASE (tblite_solver_gvr)
4282 solver = "gvr"
4283 CASE DEFAULT
4284 cpabort("Unknown tblite reference CLI solver")
4285 END SELECT
4286
4287 END FUNCTION tb_reference_solver_name
4288
4289! **************************************************************************************************
4290!> \brief Map CP2K smearing method id to a diagnostic label.
4291!> \param method_id ...
4292!> \return ...
4293! **************************************************************************************************
4294 FUNCTION tb_reference_smear_method_name(method_id) RESULT(method)
4295 INTEGER, INTENT(IN) :: method_id
4296 CHARACTER(LEN=24) :: method
4297
4298 SELECT CASE (method_id)
4299 CASE (smear_fermi_dirac)
4300 method = "FERMI_DIRAC"
4301 CASE (smear_energy_window)
4302 method = "ENERGY_WINDOW"
4303 CASE (smear_list)
4304 method = "LIST"
4305 CASE (smear_gaussian)
4306 method = "GAUSSIAN"
4307 CASE (smear_mp)
4308 method = "METHFESSEL_PAXTON"
4309 CASE (smear_mv)
4310 method = "MARZARI_VANDERBILT"
4311 CASE DEFAULT
4312 method = "UNKNOWN"
4313 END SELECT
4314
4315 END FUNCTION tb_reference_smear_method_name
4316
4317! **************************************************************************************************
4318!> \brief Run optional native tblite auxiliary subcommands.
4319!> \param ref ...
4320!> \param dft_control ...
4321!> \param gen_file ...
4322!> \param file_base ...
4323!> \param verbosity_str ...
4324!> \param iounit ...
4325! **************************************************************************************************
4326 SUBROUTINE tb_reference_cli_aux_commands(ref, dft_control, gen_file, file_base, verbosity_str, iounit)
4327
4328 TYPE(xtb_reference_cli_type), INTENT(IN) :: ref
4329 TYPE(dft_control_type), INTENT(IN) :: dft_control
4330 CHARACTER(LEN=*), INTENT(IN) :: gen_file, file_base, verbosity_str
4331 INTEGER, INTENT(IN) :: iounit
4332
4333 CHARACTER(LEN=32) :: charge_str, efield_x_str, efield_y_str, &
4334 efield_z_str, etemp_guess_val_str, &
4335 spin_str
4336 CHARACTER(LEN=4*default_path_length+16) :: copy_str, dry_run_str, efield_str, &
4337 etemp_guess_str, grad_str, json_str, &
4338 method_str, output_str
4339 CHARACTER(LEN=8*default_path_length) :: command
4340 CHARACTER(LEN=default_path_length) :: guess_input, log_file
4341 INTEGER :: spin
4342 REAL(kind=dp) :: etemp_guess
4343
4344 WRITE (charge_str, "(I0)") dft_control%charge
4345 spin = max(0, dft_control%multiplicity - 1)
4346 WRITE (spin_str, "(I0)") spin
4347
4348 IF (ref%guess_cli%enabled) THEN
4349 guess_input = ref%guess_cli%input_file
4350 IF (len_trim(guess_input) == 0) guess_input = gen_file
4351 etemp_guess_str = ""
4352 IF (ref%guess_cli%electronic_temperature_guess > 0.0_dp) THEN
4353 etemp_guess = cp_unit_from_cp2k(ref%guess_cli%electronic_temperature_guess, "K")
4354 WRITE (etemp_guess_val_str, "(ES16.8)") etemp_guess
4355 etemp_guess_str = " --etemp-guess "//trim(adjustl(etemp_guess_val_str))
4356 END IF
4357 efield_str = ""
4358 IF (ref%guess_cli%efield_active) THEN
4359 WRITE (efield_x_str, "(ES16.8)") ref%guess_cli%efield(1)
4360 WRITE (efield_y_str, "(ES16.8)") ref%guess_cli%efield(2)
4361 WRITE (efield_z_str, "(ES16.8)") ref%guess_cli%efield(3)
4362 efield_str = " --efield "//trim(adjustl(efield_x_str))//","// &
4363 trim(adjustl(efield_y_str))//","//trim(adjustl(efield_z_str))
4364 END IF
4365 grad_str = ""
4366 IF (ref%guess_cli%grad) grad_str = " --grad"
4367 json_str = ""
4368 IF (len_trim(ref%guess_cli%json_file) > 0) THEN
4369 json_str = " --json "//trim(tb_shell_quote(ref%guess_cli%json_file))
4370 END IF
4371 log_file = trim(file_base)//".guess.log"
4372 command = trim(tb_shell_quote(ref%program_name))// &
4373 " guess --charge "//trim(adjustl(charge_str))// &
4374 " --spin "//trim(adjustl(spin_str))// &
4375 " --method "//trim(tb_reference_guess_name(ref%guess_cli%method))// &
4376 " --solver "//trim(tb_reference_solver_name(ref%guess_cli%solver))// &
4377 trim(etemp_guess_str)// &
4378 trim(efield_str)// &
4379 trim(grad_str)// &
4380 trim(json_str)// &
4381 trim(verbosity_str)//" --input "//trim(tb_shell_quote(ref%guess_cli%input_format))// &
4382 " "//trim(tb_shell_quote(guess_input))// &
4383 " > "//trim(tb_shell_quote(log_file))//" 2>&1"
4384 CALL tb_reference_cli_execute(ref, "guess", command, log_file, iounit)
4385 IF (.NOT. ref%keep_files .AND. len_trim(ref%guess_cli%json_file) > 0) THEN
4386 CALL tb_delete_file(ref%guess_cli%json_file)
4387 END IF
4388 END IF
4389
4390 IF (ref%param_cli%enabled) THEN
4391 method_str = ""
4392 IF (ref%param_cli%method_explicit .OR. len_trim(ref%param_cli%input_file) == 0) THEN
4393 method_str = " --method "// &
4394 trim(tb_reference_method_name(merge(ref%param_cli%method, &
4395 dft_control%qs_control%xtb_control%tblite_method, &
4396 ref%param_cli%method_explicit)))
4397 END IF
4398 output_str = ""
4399 IF (len_trim(ref%param_cli%output_file) > 0) THEN
4400 output_str = " --output "//trim(tb_shell_quote(ref%param_cli%output_file))
4401 END IF
4402 log_file = trim(file_base)//".param.log"
4403 command = trim(tb_shell_quote(ref%program_name))//" param"// &
4404 trim(method_str)// &
4405 trim(output_str)
4406 IF (len_trim(ref%param_cli%input_file) > 0) THEN
4407 command = trim(command)//" "//trim(tb_shell_quote(ref%param_cli%input_file))
4408 END IF
4409 command = trim(command)//" > "//trim(tb_shell_quote(log_file))//" 2>&1"
4410 CALL tb_reference_cli_execute(ref, "param", command, log_file, iounit)
4411 IF (.NOT. ref%keep_files .AND. len_trim(ref%param_cli%output_file) > 0) THEN
4412 CALL tb_delete_file(ref%param_cli%output_file)
4413 END IF
4414 END IF
4415
4416 IF (ref%fit_cli%enabled) THEN
4417 dry_run_str = ""
4418 IF (ref%fit_cli%dry_run) dry_run_str = " --dry-run"
4419 copy_str = ""
4420 IF (len_trim(ref%fit_cli%copy_file) > 0) THEN
4421 copy_str = " --copy "//trim(tb_shell_quote(ref%fit_cli%copy_file))
4422 END IF
4423 log_file = trim(file_base)//".fit.log"
4424 command = trim(tb_shell_quote(ref%program_name))//" fit"// &
4425 trim(dry_run_str)// &
4426 trim(copy_str)// &
4427 trim(verbosity_str)//" "//trim(tb_shell_quote(ref%fit_cli%param_file))// &
4428 " "//trim(tb_shell_quote(ref%fit_cli%input_file))// &
4429 " > "//trim(tb_shell_quote(log_file))//" 2>&1"
4430 CALL tb_reference_cli_execute(ref, "fit", command, log_file, iounit)
4431 IF (.NOT. ref%keep_files .AND. len_trim(ref%fit_cli%copy_file) > 0) THEN
4432 CALL tb_delete_file(ref%fit_cli%copy_file)
4433 END IF
4434 END IF
4435
4436 IF (ref%tagdiff_cli%enabled) THEN
4437 method_str = ""
4438 IF (ref%tagdiff_cli%fit) method_str = " --fit"
4439 log_file = trim(file_base)//".tagdiff.log"
4440 command = trim(tb_shell_quote(ref%program_name))//" tagdiff"// &
4441 trim(method_str)//" "//trim(tb_shell_quote(ref%tagdiff_cli%actual_file))// &
4442 " "//trim(tb_shell_quote(ref%tagdiff_cli%reference_file))// &
4443 " > "//trim(tb_shell_quote(log_file))//" 2>&1"
4444 CALL tb_reference_cli_execute(ref, "tagdiff", command, log_file, iounit)
4445 END IF
4446
4447 END SUBROUTINE tb_reference_cli_aux_commands
4448
4449! **************************************************************************************************
4450!> \brief Execute one native tblite REFERENCE_CLI auxiliary command.
4451!> \param ref ...
4452!> \param label ...
4453!> \param command ...
4454!> \param log_file ...
4455!> \param iounit ...
4456! **************************************************************************************************
4457 SUBROUTINE tb_reference_cli_execute(ref, label, command, log_file, iounit)
4458
4459 TYPE(xtb_reference_cli_type), INTENT(IN) :: ref
4460 CHARACTER(LEN=*), INTENT(IN) :: label, command, log_file
4461 INTEGER, INTENT(IN) :: iounit
4462
4463 INTEGER :: cmdstat, exitstat
4464
4465 cmdstat = 0
4466 exitstat = 0
4467 CALL execute_command_line(trim(command), exitstat=exitstat, cmdstat=cmdstat)
4468 IF (cmdstat /= 0 .OR. exitstat /= 0) THEN
4469 WRITE (unit=iounit, fmt="(/,T2,A,A)") "tblite reference CLI auxiliary command failed: ", &
4470 trim(label)
4471 WRITE (unit=iounit, fmt="(T2,A,A)") "Command: ", trim(command)
4472 WRITE (unit=iounit, fmt="(T2,A,I0,T32,A,I0)") "cmdstat:", cmdstat, "exitstat:", exitstat
4473 IF (ref%stop_on_error) cpabort("tblite reference CLI auxiliary command failed")
4474 ELSE
4475 WRITE (unit=iounit, fmt="(/,T2,A,A)") "tblite reference CLI auxiliary command completed: ", &
4476 trim(label)
4477 WRITE (unit=iounit, fmt="(T2,A,A)") "Log file: ", trim(log_file)
4478 END IF
4479 IF (.NOT. ref%keep_files) CALL tb_delete_file(log_file)
4480
4481 END SUBROUTINE tb_reference_cli_execute
4482
4483! **************************************************************************************************
4484!> \brief Join directory and filename.
4485!> \param directory ...
4486!> \param filename ...
4487!> \return ...
4488! **************************************************************************************************
4489 FUNCTION tb_join_path(directory, filename) RESULT(path)
4490 CHARACTER(LEN=*), INTENT(IN) :: directory, filename
4491 CHARACTER(LEN=default_path_length) :: path
4492
4493 IF (len_trim(directory) == 0 .OR. trim(directory) == ".") THEN
4494 path = trim(filename)
4495 ELSE IF (directory(len_trim(directory):len_trim(directory)) == "/") THEN
4496 path = trim(directory)//trim(filename)
4497 ELSE
4498 path = trim(directory)//"/"//trim(filename)
4499 END IF
4500
4501 END FUNCTION tb_join_path
4502
4503! **************************************************************************************************
4504!> \brief Shell-quote a filename or executable path.
4505!> \param text ...
4506!> \return ...
4507! **************************************************************************************************
4508 FUNCTION tb_shell_quote(text) RESULT(quoted)
4509 CHARACTER(LEN=*), INTENT(IN) :: text
4510 CHARACTER(LEN=4*default_path_length) :: quoted
4511
4512 INTEGER :: i
4513
4514 quoted = "'"
4515 DO i = 1, len_trim(text)
4516 IF (text(i:i) == "'") THEN
4517 quoted = trim(quoted)//"'\\''"
4518 ELSE
4519 quoted = trim(quoted)//text(i:i)
4520 END IF
4521 END DO
4522 quoted = trim(quoted)//"'"
4523
4524 END FUNCTION tb_shell_quote
4525
4526! **************************************************************************************************
4527!> \brief Write current CP2K geometry as DFTB+ gen format for native tblite.
4528!> \param qs_env ...
4529!> \param filename ...
4530! **************************************************************************************************
4531 SUBROUTINE tb_write_reference_gen(qs_env, filename)
4532
4533 TYPE(qs_environment_type), POINTER :: qs_env
4534 CHARACTER(LEN=*), INTENT(IN) :: filename
4535
4536 CHARACTER(LEN=2), ALLOCATABLE, DIMENSION(:) :: symbols, unique_symbols
4537 INTEGER :: iatom, ikind, ios, natom, nuniq, unit_nr
4538 INTEGER, ALLOCATABLE, DIMENSION(:) :: species
4539 INTEGER, DIMENSION(3) :: periodic
4540 LOGICAL :: found
4541 REAL(kind=dp) :: to_angstrom
4542 TYPE(cell_type), POINTER :: cell
4543 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
4544
4545 NULLIFY (cell, particle_set)
4546 CALL get_qs_env(qs_env=qs_env, cell=cell, particle_set=particle_set)
4547
4548 natom = SIZE(particle_set)
4549 to_angstrom = cp_unit_from_cp2k(1.0_dp, "angstrom")
4550 ALLOCATE (symbols(natom), unique_symbols(natom), species(natom))
4551 nuniq = 0
4552 DO iatom = 1, natom
4553 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=symbols(iatom))
4554 found = .false.
4555 DO ikind = 1, nuniq
4556 IF (trim(unique_symbols(ikind)) == trim(symbols(iatom))) THEN
4557 found = .true.
4558 EXIT
4559 END IF
4560 END DO
4561 IF (.NOT. found) THEN
4562 nuniq = nuniq + 1
4563 unique_symbols(nuniq) = symbols(iatom)
4564 ikind = nuniq
4565 END IF
4566 species(iatom) = ikind
4567 END DO
4568
4569 OPEN (newunit=unit_nr, file=trim(filename), status="REPLACE", action="WRITE", &
4570 form="FORMATTED", iostat=ios)
4571 IF (ios /= 0) cpabort("Could not open tblite reference CLI geometry file")
4572
4573 CALL get_cell(cell=cell, periodic=periodic)
4574 IF (any(periodic == 1)) THEN
4575 WRITE (unit=unit_nr, fmt="(I0,1X,A)") natom, "S"
4576 ELSE
4577 WRITE (unit=unit_nr, fmt="(I0,1X,A)") natom, "C"
4578 END IF
4579 WRITE (unit=unit_nr, fmt="(*(A,1X))") (trim(unique_symbols(ikind)), ikind=1, nuniq)
4580 DO iatom = 1, natom
4581 WRITE (unit=unit_nr, fmt="(I0,1X,I0,3(1X,ES24.16))") &
4582 iatom, species(iatom), particle_set(iatom)%r(:)*to_angstrom
4583 END DO
4584 IF (any(periodic == 1)) THEN
4585 WRITE (unit=unit_nr, fmt="(3(1X,ES24.16))") 0.0_dp, 0.0_dp, 0.0_dp
4586 DO ikind = 1, 3
4587 WRITE (unit=unit_nr, fmt="(3(1X,ES24.16))") cell%hmat(:, ikind)*to_angstrom
4588 END DO
4589 END IF
4590 CLOSE (unit_nr)
4591
4592 DEALLOCATE (symbols, unique_symbols, species)
4593
4594 END SUBROUTINE tb_write_reference_gen
4595
4596! **************************************************************************************************
4597!> \brief Read native tblite gradient file.
4598!> \param filename ...
4599!> \param natom ...
4600!> \param energy ...
4601!> \param gradient ...
4602!> \param virial ...
4603!> \param have_energy ...
4604!> \param have_gradient ...
4605!> \param have_virial ...
4606! **************************************************************************************************
4607 SUBROUTINE tb_read_reference_grad(filename, natom, energy, gradient, virial, &
4608 have_energy, have_gradient, have_virial)
4609
4610 CHARACTER(LEN=*), INTENT(IN) :: filename
4611 INTEGER, INTENT(IN) :: natom
4612 REAL(kind=dp), INTENT(OUT) :: energy
4613 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: gradient, virial
4614 LOGICAL, INTENT(OUT) :: have_energy, have_gradient, have_virial
4615
4616 CHARACTER(LEN=1024) :: line
4617 INTEGER :: ios, nread, unit_nr
4618 LOGICAL :: exists
4619 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: values
4620
4621 have_energy = .false.
4622 have_gradient = .false.
4623 have_virial = .false.
4624 energy = 0.0_dp
4625 gradient = 0.0_dp
4626 virial = 0.0_dp
4627
4628 INQUIRE (file=trim(filename), exist=exists)
4629 IF (.NOT. exists) RETURN
4630
4631 OPEN (newunit=unit_nr, file=trim(filename), status="OLD", action="READ", &
4632 form="FORMATTED", iostat=ios)
4633 IF (ios /= 0) RETURN
4634
4635 DO
4636 READ (unit=unit_nr, fmt="(A)", iostat=ios) line
4637 IF (ios /= 0) EXIT
4638 IF (index(line, "energy :real:0:") > 0) THEN
4639 READ (unit=unit_nr, fmt="(A)", iostat=ios) line
4640 IF (ios == 0) THEN
4641 READ (line, *, iostat=ios) energy
4642 have_energy = ios == 0
4643 END IF
4644 ELSE IF (index(line, "gradient :real:2:3,") > 0) THEN
4645 ALLOCATE (values(3*natom))
4646 CALL tb_read_real_values(unit_nr, values, nread)
4647 IF (nread == 3*natom) THEN
4648 CALL tb_values_to_matrix(values, gradient)
4649 have_gradient = .true.
4650 END IF
4651 DEALLOCATE (values)
4652 ELSE IF (index(line, "virial :real:2:3,3") > 0) THEN
4653 ALLOCATE (values(9))
4654 CALL tb_read_real_values(unit_nr, values, nread)
4655 IF (nread == 9) THEN
4656 CALL tb_values_to_matrix(values, virial)
4657 have_virial = .true.
4658 END IF
4659 DEALLOCATE (values)
4660 END IF
4661 END DO
4662 CLOSE (unit_nr)
4663
4664 END SUBROUTINE tb_read_reference_grad
4665
4666! **************************************************************************************************
4667!> \brief Read a fixed number of real values from following lines.
4668!> \param unit_nr ...
4669!> \param values ...
4670!> \param nread ...
4671! **************************************************************************************************
4672 SUBROUTINE tb_read_real_values(unit_nr, values, nread)
4673
4674 INTEGER, INTENT(IN) :: unit_nr
4675 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: values
4676 INTEGER, INTENT(OUT) :: nread
4677
4678 CHARACTER(LEN=1024) :: line
4679 INTEGER :: ios
4680
4681 nread = 0
4682 DO WHILE (nread < SIZE(values))
4683 READ (unit=unit_nr, fmt="(A)", iostat=ios) line
4684 IF (ios /= 0) EXIT
4685 CALL tb_parse_real_line(line, values, nread)
4686 END DO
4687
4688 END SUBROUTINE tb_read_real_values
4689
4690! **************************************************************************************************
4691!> \brief Parse real values from one text line.
4692!> \param line ...
4693!> \param values ...
4694!> \param nread ...
4695! **************************************************************************************************
4696 SUBROUTINE tb_parse_real_line(line, values, nread)
4697
4698 CHARACTER(LEN=*), INTENT(IN) :: line
4699 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: values
4700 INTEGER, INTENT(INOUT) :: nread
4701
4702 CHARACTER(LEN=128) :: token
4703 INTEGER :: first, ios, last, pos
4704
4705 pos = 1
4706 DO WHILE (pos <= len_trim(line) .AND. nread < SIZE(values))
4707 DO WHILE (pos <= len_trim(line) .AND. index(" ,[]", line(pos:pos)) > 0)
4708 pos = pos + 1
4709 END DO
4710 IF (pos > len_trim(line)) EXIT
4711 first = pos
4712 DO WHILE (pos <= len_trim(line) .AND. index(" ,[]", line(pos:pos)) == 0)
4713 pos = pos + 1
4714 END DO
4715 last = pos - 1
4716 token = line(first:last)
4717 READ (token, *, iostat=ios) values(nread + 1)
4718 IF (ios == 0) nread = nread + 1
4719 END DO
4720
4721 END SUBROUTINE tb_parse_real_line
4722
4723! **************************************************************************************************
4724!> \brief Convert flat native tblite values to CP2K atom-major matrix layout.
4725!> \param values ...
4726!> \param matrix ...
4727! **************************************************************************************************
4728 SUBROUTINE tb_values_to_matrix(values, matrix)
4729
4730 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: values
4731 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: matrix
4732
4733 INTEGER :: i, j, n
4734
4735 n = 0
4736 DO j = 1, SIZE(matrix, 2)
4737 DO i = 1, SIZE(matrix, 1)
4738 n = n + 1
4739 matrix(i, j) = values(n)
4740 END DO
4741 END DO
4742
4743 END SUBROUTINE tb_values_to_matrix
4744
4745! **************************************************************************************************
4746!> \brief Delete temporary files unless requested otherwise.
4747!> \param ref ...
4748!> \param gen_file ...
4749!> \param grad_file ...
4750!> \param json_file ...
4751!> \param log_file ...
4752!> \param post_processing_output_file ...
4753! **************************************************************************************************
4754 SUBROUTINE tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file, post_processing_output_file)
4755
4756 TYPE(xtb_reference_cli_type), INTENT(IN) :: ref
4757 CHARACTER(LEN=*), INTENT(IN) :: gen_file, grad_file, json_file, &
4758 log_file, post_processing_output_file
4759
4760 IF (ref%keep_files) RETURN
4761 CALL tb_delete_file(gen_file)
4762 CALL tb_delete_file(grad_file)
4763 CALL tb_delete_file(json_file)
4764 CALL tb_delete_file(log_file)
4765 IF (len_trim(post_processing_output_file) > 0) THEN
4766 CALL tb_delete_file(post_processing_output_file)
4767 END IF
4768
4769 END SUBROUTINE tb_reference_cleanup
4770
4771! **************************************************************************************************
4772!> \brief Delete a file if it exists.
4773!> \param filename ...
4774! **************************************************************************************************
4775 SUBROUTINE tb_delete_file(filename)
4776
4777 CHARACTER(LEN=*), INTENT(IN) :: filename
4778
4779 INTEGER :: ios, unit_nr
4780 LOGICAL :: exists
4781
4782 INQUIRE (file=trim(filename), exist=exists)
4783 IF (.NOT. exists) RETURN
4784 OPEN (newunit=unit_nr, file=trim(filename), status="OLD", iostat=ios)
4785 IF (ios == 0) CLOSE (unit_nr, status="DELETE")
4786
4787 END SUBROUTINE tb_delete_file
4788
4789! **************************************************************************************************
4790!> \brief Dump cumulative tblite virial pieces for local debugging.
4791!> \param label ...
4792!> \param sigma ...
4793!> \param para_env ...
4794! **************************************************************************************************
4795 SUBROUTINE tb_dump_sigma_component(label, sigma, para_env)
4796
4797 CHARACTER(LEN=*), INTENT(IN) :: label
4798 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: sigma
4799 TYPE(mp_para_env_type), INTENT(IN) :: para_env
4800
4801 CHARACTER(LEN=default_path_length) :: dump_file
4802 INTEGER :: dump_status, dump_unit, i
4803
4804 dump_status = 1
4805#if defined(__TBLITE_DEBUG_DIAGNOSTICS)
4806 CALL get_environment_variable("CP2K_TBLITE_SIGMA_COMPONENT_DUMP", dump_file, status=dump_status)
4807#endif
4808 IF (dump_status /= 0) RETURN
4809
4810 OPEN (newunit=dump_unit, file=trim(dump_file), status="UNKNOWN", &
4811 position="APPEND", action="WRITE")
4812 WRITE (dump_unit, "(A)") trim(label)
4813 DO i = 1, 3
4814 WRITE (dump_unit, "(3(1X,ES24.16))") sigma(i, :)/para_env%num_pe
4815 END DO
4816 CLOSE (dump_unit)
4817
4818 END SUBROUTINE tb_dump_sigma_component
4819
4820! **************************************************************************************************
4821!> \brief save stress tensor
4822!> \param qs_env ...
4823!> \param tb ...
4824!> \param para_env ...
4825! **************************************************************************************************
4826 SUBROUTINE tb_add_stress(qs_env, tb, para_env)
4827
4828 TYPE(qs_environment_type) :: qs_env
4829 TYPE(tblite_type) :: tb
4830 TYPE(mp_para_env_type) :: para_env
4831
4832 CHARACTER(LEN=default_path_length) :: dump_file
4833 INTEGER :: dump_status, dump_unit, i
4834 INTEGER, DIMENSION(3) :: periodic
4835 TYPE(cell_type), POINTER :: cell
4836 TYPE(virial_type), POINTER :: virial
4837
4838 NULLIFY (virial, cell)
4839 CALL get_qs_env(qs_env=qs_env, virial=virial, cell=cell)
4840 CALL get_cell(cell=cell, periodic=periodic)
4841
4842 IF (all(periodic == 0)) THEN
4843 CALL cp_warn(__location__, &
4844 "tblite stress tensor requested for an isolated system. "// &
4845 "The reported virial is useful for finite-difference checks, "// &
4846 "but it is not a physically meaningful bulk stress for an isolated molecule.")
4847 END IF
4848
4849 dump_status = 1
4850#if defined(__TBLITE_DEBUG_DIAGNOSTICS)
4851 CALL get_environment_variable("CP2K_TBLITE_VIRIAL_DUMP", dump_file, status=dump_status)
4852#endif
4853 IF (dump_status == 0) THEN
4854 OPEN (newunit=dump_unit, file=trim(dump_file), status="UNKNOWN", &
4855 position="APPEND", action="WRITE")
4856 WRITE (dump_unit, "(A)") "sigma"
4857 DO i = 1, 3
4858 WRITE (dump_unit, "(3(1X,ES24.16))") tb%sigma(i, :)/para_env%num_pe
4859 END DO
4860 CLOSE (dump_unit)
4861 END IF
4862
4863 virial%pv_virial = virial%pv_virial - tb%sigma/para_env%num_pe
4864
4865 END SUBROUTINE tb_add_stress
4866
4867! **************************************************************************************************
4868!> \brief add contrib. to gradient
4869!> \param grad ...
4870!> \param deriv ...
4871!> \param dE ...
4872!> \param natom ...
4873! **************************************************************************************************
4874 SUBROUTINE tb_add_grad(grad, deriv, dE, natom)
4875
4876 REAL(kind=dp), DIMENSION(:, :) :: grad
4877 REAL(kind=dp), DIMENSION(:, :, :) :: deriv
4878 REAL(kind=dp), DIMENSION(:) :: de
4879 INTEGER :: natom
4880
4881 INTEGER :: i, j
4882
4883 DO i = 1, natom
4884 DO j = 1, natom
4885 grad(:, i) = grad(:, i) + deriv(:, i, j)*de(j)
4886 END DO
4887 END DO
4888
4889 END SUBROUTINE tb_add_grad
4890
4891! **************************************************************************************************
4892!> \brief add contrib. to sigma
4893!> \param sig ...
4894!> \param deriv ...
4895!> \param dE ...
4896!> \param natom ...
4897! **************************************************************************************************
4898 SUBROUTINE tb_add_sig(sig, deriv, dE, natom)
4899
4900 REAL(kind=dp), DIMENSION(:, :) :: sig
4901 REAL(kind=dp), DIMENSION(:, :, :) :: deriv
4902 REAL(kind=dp), DIMENSION(:) :: de
4903 INTEGER :: natom
4904
4905 INTEGER :: i, j
4906
4907 DO i = 1, 3
4908 DO j = 1, natom
4909 sig(:, i) = sig(:, i) + deriv(:, i, j)*de(j)
4910 END DO
4911 END DO
4912
4913 END SUBROUTINE tb_add_sig
4914
4915END MODULE tblite_interface
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:161
ai_contraction::block_add
Definition ai_contraction.F:43
ai_contraction::contraction
Definition ai_contraction.F:31
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:95
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:103
memory_utilities::reallocate
Definition memory_utilities.F:27
mulliken::ao_charges
Definition mulliken.F:42
qs_integral_utils::get_memory_usage
Definition qs_integral_utils.F:34
qs_overlap::create_sab_matrix
Definition qs_overlap.F:81
ai_contraction
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Definition ai_contraction.F:18
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
ai_overlap::overlap_ab
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab, rr_work)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:273
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
basis_set_types
Definition basis_set_types.F:15
basis_set_types::process_gto_basis
subroutine, public process_gto_basis(gto_basis_set, do_ortho, nset, maxl)
...
Definition basis_set_types.F:2512
basis_set_types::allocate_gto_basis_set
subroutine, public allocate_gto_basis_set(gto_basis_set)
...
Definition basis_set_types.F:152
basis_set_types::write_gto_basis_set
subroutine, public write_gto_basis_set(gto_basis_set, output_unit, header)
Write a Gaussian-type orbital (GTO) basis set data set to the output unit.
Definition basis_set_types.F:1924
block_p_types
collect pointers to a block of reals
Definition block_p_types.F:14
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:233
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_print
subroutine, public dbcsr_print(matrix, variable_name, unit_nr)
Prints given matrix in matlab format (only present blocks).
Definition cp_dbcsr_contrib.F:682
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_output
DBCSR output in CP2K.
Definition cp_dbcsr_output.F:17
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers, cartesian_basis)
...
Definition cp_dbcsr_output.F:165
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::debug_print_level
integer, parameter, public debug_print_level
Definition cp_output_handling.F:73
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_output_handling::high_print_level
integer, parameter, public high_print_level
Definition cp_output_handling.F:73
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_output_handling::silent_print_level
integer, parameter, public silent_print_level
Definition cp_output_handling.F:73
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1251
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::smear_fermi_dirac
integer, parameter, public smear_fermi_dirac
Definition input_constants.F:950
input_constants::gfn1xtb
integer, parameter, public gfn1xtb
Definition input_constants.F:1429
input_constants::smear_energy_window
integer, parameter, public smear_energy_window
Definition input_constants.F:950
input_constants::tblite_scc_mixer_cp2k
integer, parameter, public tblite_scc_mixer_cp2k
Definition input_constants.F:1435
input_constants::tblite_scc_mixer_none
integer, parameter, public tblite_scc_mixer_none
Definition input_constants.F:1435
input_constants::tblite_guess_sad
integer, parameter, public tblite_guess_sad
Definition input_constants.F:1450
input_constants::tblite_cli_born_kernel_p16
integer, parameter, public tblite_cli_born_kernel_p16
Definition input_constants.F:1466
input_constants::tblite_solver_gvd
integer, parameter, public tblite_solver_gvd
Definition input_constants.F:1455
input_constants::tblite_mixer_damping_default
real(kind=dp), parameter, public tblite_mixer_damping_default
Definition input_constants.F:1443
input_constants::tblite_cli_solution_state_gsolv
integer, parameter, public tblite_cli_solution_state_gsolv
Definition input_constants.F:1470
input_constants::tblite_guess_eeq
integer, parameter, public tblite_guess_eeq
Definition input_constants.F:1450
input_constants::tblite_cli_solution_state_reference
integer, parameter, public tblite_cli_solution_state_reference
Definition input_constants.F:1470
input_constants::tblite_cli_born_kernel_still
integer, parameter, public tblite_cli_born_kernel_still
Definition input_constants.F:1466
input_constants::smear_list
integer, parameter, public smear_list
Definition input_constants.F:950
input_constants::ipea1xtb
integer, parameter, public ipea1xtb
Definition input_constants.F:1429
input_constants::tblite_cli_solvation_gb
integer, parameter, public tblite_cli_solvation_gb
Definition input_constants.F:1459
input_constants::tblite_solver_gvr
integer, parameter, public tblite_solver_gvr
Definition input_constants.F:1455
input_constants::tblite_scc_mixer_tblite
integer, parameter, public tblite_scc_mixer_tblite
Definition input_constants.F:1435
input_constants::tblite_cli_solution_state_bar1mol
integer, parameter, public tblite_cli_solution_state_bar1mol
Definition input_constants.F:1470
input_constants::tblite_cli_solvation_alpb
integer, parameter, public tblite_cli_solvation_alpb
Definition input_constants.F:1459
input_constants::smear_gaussian
integer, parameter, public smear_gaussian
Definition input_constants.F:950
input_constants::tblite_mixer_max_weight_default
real(kind=dp), parameter, public tblite_mixer_max_weight_default
Definition input_constants.F:1443
input_constants::tblite_cli_born_kernel_auto
integer, parameter, public tblite_cli_born_kernel_auto
Definition input_constants.F:1466
input_constants::tblite_mixer_omega0_default
real(kind=dp), parameter, public tblite_mixer_omega0_default
Definition input_constants.F:1443
input_constants::smear_mv
integer, parameter, public smear_mv
Definition input_constants.F:950
input_constants::tblite_scc_mixer_auto
integer, parameter, public tblite_scc_mixer_auto
Definition input_constants.F:1435
input_constants::tblite_cli_solvation_cpcm
integer, parameter, public tblite_cli_solvation_cpcm
Definition input_constants.F:1459
input_constants::tblite_cli_solvation_gbe
integer, parameter, public tblite_cli_solvation_gbe
Definition input_constants.F:1459
input_constants::tblite_cli_solvation_gbsa
integer, parameter, public tblite_cli_solvation_gbsa
Definition input_constants.F:1459
input_constants::smear_mp
integer, parameter, public smear_mp
Definition input_constants.F:950
input_constants::gfn2xtb
integer, parameter, public gfn2xtb
Definition input_constants.F:1429
input_constants::tblite_guess_ceh
integer, parameter, public tblite_guess_ceh
Definition input_constants.F:1450
input_constants::tblite_mixer_weight_factor_default
real(kind=dp), parameter, public tblite_mixer_weight_factor_default
Definition input_constants.F:1443
input_constants::tblite_mixer_min_weight_default
real(kind=dp), parameter, public tblite_mixer_min_weight_default
Definition input_constants.F:1443
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:585
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
mulliken
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii}...
Definition mulliken.F:13
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_charge_mixing
Definition qs_charge_mixing.F:9
qs_charge_mixing::charge_mixing
subroutine, public charge_mixing(mixing_method, mixing_store, charges, para_env, iter_count, scc_mixer, tblite_mixer_iterations, tblite_mixer_damping, tblite_mixer_memory, tblite_mixer_omega0, tblite_mixer_min_weight, tblite_mixer_max_weight, tblite_mixer_weight_factor)
Driver for TB SCC variable mixing, calls the requested method.
Definition qs_charge_mixing.F:71
qs_charge_mixing::tblite_scc_error_on_cp2k_scale
pure real(kind=dp) function, public tblite_scc_error_on_cp2k_scale(raw_error, eps_scf, pconv)
Map a raw tblite SCC residual to CP2K's EPS_SCF reporting scale.
Definition qs_charge_mixing.F:226
qs_charge_mixing::tblite_scc_pconv
real(kind=dp), parameter, public tblite_scc_pconv
Definition qs_charge_mixing.F:45
qs_condnum
Calculation of overlap matrix condition numbers.
Definition qs_condnum.F:13
qs_condnum::overlap_condnum
subroutine, public overlap_condnum(matrixkp_s, condnum, iunit, norml1, norml2, use_arnoldi, blacs_env)
Calculation of the overlap matrix Condition Number.
Definition qs_condnum.F:66
qs_density_mixing_types
module that contains the definitions of the scf types
Definition qs_density_mixing_types.F:14
qs_density_mixing_types::modified_broyden_mixing_nr
integer, parameter, public modified_broyden_mixing_nr
Definition qs_density_mixing_types.F:50
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:545
qs_force_types
Definition qs_force_types.F:13
qs_force_types::total_qs_force
subroutine, public total_qs_force(force, qs_force, atomic_kind_set)
Get current total force.
Definition qs_force_types.F:581
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:984
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:591
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition qs_ks_types.F:343
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:550
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
qs_overlap::build_overlap_matrix
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p, ext_kpoints)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:121
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::integer_to_string
subroutine, public integer_to_string(inumber, string)
Converts an integer number to a string. The WRITE statement will return an error message,...
Definition string_utilities.F:3239
tblite_interface
interface to tblite
Definition tblite_interface.F:14
tblite_interface::tb_native_scc_mixer_active
logical function, public tb_native_scc_mixer_active(dft_control)
Return whether the tblite native SCC mixer is active for this run.
Definition tblite_interface.F:509
tblite_interface::tb_ham_add_coulomb
subroutine, public tb_ham_add_coulomb(qs_env, tb, dft_control)
...
Definition tblite_interface.F:2175
tblite_interface::tb_init_wf
subroutine, public tb_init_wf(tb, dft_control)
initialize wavefunction ...
Definition tblite_interface.F:321
tblite_interface::tb_update_charges
subroutine, public tb_update_charges(qs_env, dft_control, tb, calculate_forces, use_rho)
...
Definition tblite_interface.F:1516
tblite_interface::tb_init_geometry
subroutine, public tb_init_geometry(qs_env, tb)
intialize geometry objects ...
Definition tblite_interface.F:210
tblite_interface::build_tblite_matrices
subroutine, public build_tblite_matrices(qs_env, calculate_forces)
...
Definition tblite_interface.F:911
tblite_interface::tb_get_energy
subroutine, public tb_get_energy(qs_env, tb, energy)
...
Definition tblite_interface.F:721
tblite_interface::tb_set_calculator
subroutine, public tb_set_calculator(tb, typ, accuracy, param_file)
...
Definition tblite_interface.F:575
tblite_interface::tb_derive_dh_off
subroutine, public tb_derive_dh_off(qs_env, use_rho, nimg)
Add SCC-overlap and direct multipole Hamiltonian derivatives.
Definition tblite_interface.F:3400
tblite_interface::tb_get_multipole
subroutine, public tb_get_multipole(qs_env, tb)
...
Definition tblite_interface.F:2611
tblite_interface::tb_reference_cli_compare
subroutine, public tb_reference_cli_compare(qs_env)
Run native tblite CLI and compare against CP2K/tblite.
Definition tblite_interface.F:3816
tblite_interface::tb_scf_mixer_error
real(kind=dp) function, public tb_scf_mixer_error(dft_control, tb, eps_scf)
Return the native tblite SCC mixer residual on the CP2K iter_delta scale.
Definition tblite_interface.F:538
tblite_interface::tb_get_basis
subroutine, public tb_get_basis(tb, gto_basis_set, element_symbol, param, occ)
...
Definition tblite_interface.F:796
tblite_scc_mixer
CP2K-side tblite-compatible SCC Broyden mixer.
Definition tblite_scc_mixer.F:11
tblite_types
types for tblite
Definition tblite_types.F:14
tblite_types::allocate_tblite_type
subroutine, public allocate_tblite_type(tb_tblite)
...
Definition tblite_types.F:102
tblite_types::deallocate_tblite_type
subroutine, public deallocate_tblite_type(tb_tblite)
...
Definition tblite_types.F:116
virial_types
Definition virial_types.F:13
xtb_types
Definition of the xTB parameter types.
Definition xtb_types.F:20
xtb_types::get_xtb_atom_param
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, wall, hen, zeta, xi, kappa0, alpg, occupation, ngauss, electronegativity, chmax, en, kqat2, kcn, kq)
...
Definition xtb_types.F:206
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
block_p_types::block_p_type
Definition block_p_types.F:23
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:780
cp_control_types::xtb_reference_cli_type
Definition cp_control_types.F:224
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:196
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
particle_types::particle_type
Definition particle_types.F:35
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:224
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:99
scf_control_types::scf_control_type
Definition scf_control_types.F:134
tblite_types::tblite_type
Definition tblite_types.F:46
virial_types::virial_type
Definition virial_types.F:26
xtb_types::xtb_atom_type
Definition xtb_types.F:42