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