Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief interface to tblite
10 : !> \author JVP
11 : !> \history creation 09.2024
12 : ! **************************************************************************************************
13 :
14 : MODULE tblite_interface
15 :
16 : #if defined(__TBLITE)
17 : USE mctc_env, ONLY: error_type
18 : USE mctc_io, ONLY: 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 :
138 : CONTAINS
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 46551350 : 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 46551350 : value = values(1)
154 46551350 : IF (SIZE(values) > 1) THEN
155 434982 : SELECT CASE (ispin)
156 : CASE (1)
157 144994 : value = values(1) + values(2)
158 : CASE (2)
159 144994 : value = values(1) - values(2)
160 : CASE DEFAULT
161 289988 : value = values(1)
162 : END SELECT
163 : END IF
164 :
165 46551350 : 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 4870 : 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 4870 : nspin = SIZE(matrix_p, 1)
180 4870 : nimg = SIZE(matrix_p, 2)
181 4870 : IF (ASSOCIATED(tb%rho_ao_kp_ref)) THEN
182 4842 : IF (SIZE(tb%rho_ao_kp_ref, 1) /= nspin .OR. SIZE(tb%rho_ao_kp_ref, 2) /= nimg) &
183 0 : CALL dbcsr_deallocate_matrix_set(tb%rho_ao_kp_ref)
184 : END IF
185 4870 : IF (.NOT. ASSOCIATED(tb%rho_ao_kp_ref)) THEN
186 28 : CALL dbcsr_allocate_matrix_set(tb%rho_ao_kp_ref, nspin, nimg)
187 56 : DO img = 1, nimg
188 112 : DO ispin = 1, nspin
189 84 : ALLOCATE (tb%rho_ao_kp_ref(ispin, img)%matrix)
190 : END DO
191 : END DO
192 : END IF
193 9740 : DO img = 1, nimg
194 19480 : DO ispin = 1, nspin
195 14610 : CALL dbcsr_copy(tb%rho_ao_kp_ref(ispin, img)%matrix, matrix_p(ispin, img)%matrix)
196 : END DO
197 : END DO
198 :
199 4870 : END SUBROUTINE tb_store_density_ref
200 : #endif
201 :
202 : ! **************************************************************************************************
203 : !> \brief intialize geometry objects ...
204 : !> \param qs_env ...
205 : !> \param tb ...
206 : ! **************************************************************************************************
207 152 : 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 152 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
218 152 : 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 152 : CALL timeset(routineN, handle)
226 :
227 : !get info from environment vaiarable
228 152 : 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 152 : natom = SIZE(particle_set)
232 456 : ALLOCATE (xyz(3, natom))
233 152 : CALL allocate_tblite_type(tb)
234 456 : ALLOCATE (tb%el_num(natom))
235 806 : tb%el_num = -9
236 806 : DO iatom = 1, natom
237 2616 : xyz(:, iatom) = particle_set(iatom)%r(:)
238 654 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
239 654 : CALL get_qs_kind(qs_kind_set(ikind), zatom=tb%el_num(iatom))
240 1460 : IF (tb%el_num(iatom) < 1 .OR. tb%el_num(iatom) > 85) THEN
241 0 : CPABORT("only elements 1-85 are supported by tblite")
242 : END IF
243 : END DO
244 :
245 : !get information about cell / lattice
246 152 : CALL get_cell(cell=cell, periodic=periodic)
247 152 : lperiod(1) = periodic(1) == 1
248 152 : lperiod(2) = periodic(2) == 1
249 152 : lperiod(3) = periodic(3) == 1
250 :
251 : !prepare for the call to the dispersion function
252 152 : CALL new(tb%mol, tb%el_num, xyz, lattice=cell%hmat, periodic=lperiod)
253 :
254 152 : DEALLOCATE (xyz)
255 :
256 152 : CALL timestop(handle)
257 :
258 : #else
259 : MARK_USED(qs_env)
260 : MARK_USED(tb)
261 : CPABORT("Built without TBLITE")
262 : #endif
263 :
264 152 : END SUBROUTINE tb_init_geometry
265 :
266 : ! **************************************************************************************************
267 : !> \brief updating coordinates...
268 : !> \param qs_env ...
269 : !> \param tb ...
270 : ! **************************************************************************************************
271 2396 : 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 2396 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
282 : INTEGER :: iatom, natom
283 2396 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xyz
284 : INTEGER :: handle
285 :
286 2396 : CALL timeset(routineN, handle)
287 :
288 : !get info from environment vaiarable
289 2396 : NULLIFY (cell)
290 2396 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell)
291 :
292 : !get information about particles
293 2396 : natom = SIZE(particle_set)
294 7188 : ALLOCATE (xyz(3, natom))
295 12920 : DO iatom = 1, natom
296 44492 : xyz(:, iatom) = particle_set(iatom)%r(:)
297 : END DO
298 44492 : tb%mol%xyz(:, :) = xyz
299 31148 : tb%mol%lattice(:, :) = cell%hmat
300 :
301 2396 : DEALLOCATE (xyz)
302 :
303 2396 : CALL timestop(handle)
304 :
305 : #else
306 : MARK_USED(qs_env)
307 : MARK_USED(tb)
308 : CPABORT("Built without TBLITE")
309 : #endif
310 :
311 2396 : END SUBROUTINE tb_update_geometry
312 :
313 : ! **************************************************************************************************
314 : !> \brief initialize wavefunction ...
315 : !> \param tb ...
316 : !> \param dft_control ...
317 : ! **************************************************************************************************
318 152 : 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 152 : nSpin = dft_control%nspins
330 0 : IF (nSpin /= 1 .AND. nSpin /= 2) CPABORT("tblite supports only one or two spin channels")
331 :
332 152 : tb%mol%charge = dft_control%charge
333 152 : tb%mol%uhf = MAX(0, dft_control%multiplicity - 1)
334 152 : IF (nSpin == 2) CALL tb_add_spin_polarization(tb)
335 :
336 152 : info = tb%calc%variable_info()
337 152 : IF (info%charge > shell_resolved) CPABORT("tblite: no support for orbital resolved charge")
338 152 : IF (info%dipole > atom_resolved) CPABORT("tblite: no support for shell resolved dipole moment")
339 152 : IF (info%quadrupole > atom_resolved) &
340 0 : CPABORT("tblite: no support shell resolved quadrupole moment")
341 :
342 152 : CALL new_wavefunction(tb%wfn, tb%mol%nat, tb%calc%bas%nsh, tb%calc%bas%nao, nSpin, 0.0_dp)
343 152 : CALL get_occupation(tb%mol, tb%calc%bas, tb%calc%h0, tb%wfn%nocc, tb%wfn%n0at, tb%wfn%n0sh)
344 152 : CALL tb_reset_mixer(tb)
345 :
346 152 : CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
347 :
348 : !allocate quantities later required
349 1064 : ALLOCATE (tb%e_hal(tb%mol%nat), tb%e_rep(tb%mol%nat), tb%e_disp(tb%mol%nat))
350 1064 : ALLOCATE (tb%e_scd(tb%mol%nat), tb%e_es(tb%mol%nat), tb%e_int(tb%mol%nat))
351 456 : ALLOCATE (tb%selfenergy(tb%calc%bas%nsh))
352 456 : 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 152 : 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 28 : SUBROUTINE tb_add_spin_polarization(tb)
368 :
369 : TYPE(tblite_type), POINTER :: tb
370 :
371 28 : CLASS(container_type), ALLOCATABLE :: cont
372 28 : TYPE(spin_polarization), ALLOCATABLE :: spin
373 28 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wll
374 :
375 28 : ALLOCATE (spin)
376 28 : CALL tb_get_spin_constants(tb, wll)
377 28 : CALL new_spin_polarization(spin, tb%mol, wll, tb%calc%bas%nsh_id)
378 28 : CALL MOVE_ALLOC(spin, cont)
379 28 : CALL tb%calc%push_back(cont)
380 :
381 28 : 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 28 : 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 140 : ALLOCATE (wll(tb%calc%bas%nsh, tb%calc%bas%nsh, tb%mol%nid))
397 28 : wll = 0.0_dp
398 64 : DO izp = 1, tb%mol%nid
399 136 : DO ish = 1, tb%calc%bas%nsh_id(izp)
400 72 : il = tb%calc%bas%cgto(ish, izp)%ang
401 256 : DO jsh = 1, tb%calc%bas%nsh_id(izp)
402 148 : jl = tb%calc%bas%cgto(jsh, izp)%ang
403 220 : wll(jsh, ish, izp) = get_spin_constant(jl, il, tb%mol%num(izp))
404 : END DO
405 : END DO
406 : END DO
407 :
408 28 : 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 2332 : 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 2332 : info = tb%calc%variable_info()
424 4512 : 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 2332 : tb%mixer_max_weight, tb%mixer_weight_factor)
429 :
430 : #else
431 : MARK_USED(tb)
432 : CPABORT("Built without TBLITE")
433 : #endif
434 :
435 2332 : 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 2180 : 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 2180 : IF (iterations < 1) CPABORT("tblite SCC mixer ITERATIONS must be positive")
460 2180 : IF (memory < 1) CPABORT("tblite SCC mixer MEMORY must be positive")
461 2180 : IF (damping <= 0.0_dp) CPABORT("tblite SCC mixer damping must be positive")
462 2180 : IF (omega0 <= 0.0_dp) CPABORT("tblite SCC mixer OMEGA0 must be positive")
463 2180 : IF (min_weight <= 0.0_dp) CPABORT("tblite SCC mixer MIN_WEIGHT must be positive")
464 2180 : IF (max_weight <= 0.0_dp) CPABORT("tblite SCC mixer MAX_WEIGHT must be positive")
465 2180 : IF (max_weight < min_weight) &
466 0 : CPABORT("tblite SCC mixer MAX_WEIGHT must not be smaller than MIN_WEIGHT")
467 2180 : IF (weight_factor <= 0.0_dp) CPABORT("tblite SCC mixer WEIGHT_FACTOR must be positive")
468 2180 : SELECT CASE (solver)
469 : CASE (tblite_solver_gvd, tblite_solver_gvr)
470 : CASE DEFAULT
471 2180 : CPABORT("Unknown tblite SCC mixer SOLVER")
472 : END SELECT
473 :
474 2180 : tb%calc%max_iter = iterations
475 2180 : tb%mixer_memory = memory
476 2180 : tb%mixer_solver = solver
477 2180 : tb%mixer_damping = damping
478 2180 : tb%calc%mixer_damping = damping
479 2180 : tb%mixer_omega0 = omega0
480 2180 : tb%mixer_min_weight = min_weight
481 2180 : tb%mixer_max_weight = max_weight
482 2180 : 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 2180 : 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 91612 : 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 91612 : use_native_mixer = .FALSE.
510 91612 : IF (.NOT. ASSOCIATED(dft_control)) RETURN
511 91612 : IF (dft_control%qs_control%do_ls_scf) RETURN
512 :
513 91612 : 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 3354 : use_native_mixer = .TRUE.
518 : CASE (tblite_scc_mixer_cp2k, tblite_scc_mixer_none)
519 3354 : use_native_mixer = .FALSE.
520 : CASE DEFAULT
521 91612 : 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 23376 : 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 23376 : mixer_error = 0.0_dp
545 :
546 : #if defined(__TBLITE)
547 23376 : IF (.NOT. ASSOCIATED(tb)) RETURN
548 23376 : IF (.NOT. tb_native_scc_mixer_active(dft_control)) RETURN
549 21752 : IF (ALLOCATED(tb%mixer)) THEN
550 21752 : raw_error = REAL(tb%mixer%get_error(), KIND=dp)
551 21752 : IF (eps_scf > 0.0_dp) THEN
552 : mixer_error = eps_scf*raw_error/ &
553 21752 : (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 152 : 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 152 : TYPE(error_type), ALLOCATABLE :: error
583 :
584 152 : IF (ALLOCATED(tb%param)) DEALLOCATE (tb%param)
585 152 : IF (LEN_TRIM(param_file) > 0) THEN
586 0 : ALLOCATE (tb%param)
587 0 : CALL tb%param%load(TRIM(param_file), error)
588 0 : IF (ALLOCATED(error)) CPABORT("Could not load tblite PARAM file: "//TRIM(param_file))
589 0 : CALL new_xtb_calculator(tb%calc, tb%mol, tb%param, error)
590 : ELSE
591 152 : SELECT CASE (typ)
592 : CASE default
593 0 : CPABORT("Unknown xtb type")
594 : CASE (gfn1xtb)
595 54 : CALL new_gfn1_calculator(tb%calc, tb%mol, error)
596 : CASE (gfn2xtb)
597 82 : CALL new_gfn2_calculator(tb%calc, tb%mol, error)
598 : CASE (ipea1xtb)
599 152 : CALL new_ipea1_calculator(tb%calc, tb%mol, error)
600 : END SELECT
601 : END IF
602 152 : IF (ALLOCATED(error)) CPABORT("tblite calculator setup failed")
603 :
604 152 : 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 152 : END SUBROUTINE tb_set_calculator
615 :
616 : ! **************************************************************************************************
617 : !> \brief ...
618 : !> \param qs_env ...
619 : !> \param tb ...
620 : !> \param para_env ...
621 : ! **************************************************************************************************
622 2396 : 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 2396 : TYPE(container_cache) :: hcache, rcache
631 :
632 12920 : tb%e_hal = 0.0_dp
633 12920 : tb%e_rep = 0.0_dp
634 12920 : tb%e_disp = 0.0_dp
635 12920 : tb%e_int = 0.0_dp
636 2396 : IF (ALLOCATED(tb%grad)) THEN
637 1946 : tb%grad = 0.0_dp
638 114 : CALL tb_zero_force(qs_env)
639 : END IF
640 31148 : tb%sigma = 0.0_dp
641 :
642 2396 : IF (ALLOCATED(tb%calc%halogen)) THEN
643 1172 : CALL tb%calc%halogen%update(tb%mol, hcache)
644 1172 : IF (ALLOCATED(tb%grad)) THEN
645 884 : tb%grad = 0.0_dp
646 : CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal, &
647 52 : & tb%grad, tb%sigma)
648 : CALL tb_dump_sigma_component("after_halogen", tb%sigma, para_env)
649 52 : CALL tb_grad2force(qs_env, tb, para_env, 0)
650 : ELSE
651 1120 : CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal)
652 : END IF
653 : END IF
654 :
655 2396 : IF (ALLOCATED(tb%calc%repulsion)) THEN
656 2396 : CALL tb%calc%repulsion%update(tb%mol, rcache)
657 2396 : IF (ALLOCATED(tb%grad)) THEN
658 1946 : tb%grad = 0.0_dp
659 : CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep, &
660 114 : & tb%grad, tb%sigma)
661 : CALL tb_dump_sigma_component("after_repulsion", tb%sigma, para_env)
662 114 : CALL tb_grad2force(qs_env, tb, para_env, 1)
663 : ELSE
664 2282 : CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep)
665 : END IF
666 : END IF
667 :
668 2396 : IF (ALLOCATED(tb%calc%dispersion)) THEN
669 2396 : CALL tb%calc%dispersion%update(tb%mol, tb%dcache)
670 2396 : IF (ALLOCATED(tb%grad)) THEN
671 1946 : tb%grad = 0.0_dp
672 : CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp, &
673 114 : & tb%grad, tb%sigma)
674 : CALL tb_dump_sigma_component("after_dispersion_static", tb%sigma, para_env)
675 114 : CALL tb_grad2force(qs_env, tb, para_env, 2)
676 : ELSE
677 2282 : CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp)
678 : END IF
679 : END IF
680 :
681 2396 : IF (ALLOCATED(tb%calc%interactions)) THEN
682 286 : CALL tb%calc%interactions%update(tb%mol, tb%icache)
683 : END IF
684 :
685 2396 : CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
686 2396 : IF (ALLOCATED(tb%calc%coulomb)) THEN
687 2396 : CALL tb%calc%coulomb%update(tb%mol, tb%cache)
688 : END IF
689 :
690 2396 : IF (ALLOCATED(tb%grad)) THEN
691 114 : IF (ALLOCATED(tb%calc%ncoord)) THEN
692 114 : 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 114 : & tb%calc%bas%nsh_id, cn=tb%cn, selfenergy=tb%selfenergy, dsedcn=tb%dsedcn)
696 : ELSE
697 2282 : IF (ALLOCATED(tb%calc%ncoord)) THEN
698 2282 : 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 2282 : & 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 2396 : END SUBROUTINE tb_init_ham
712 :
713 : ! **************************************************************************************************
714 : !> \brief ...
715 : !> \param qs_env ...
716 : !> \param tb ...
717 : !> \param energy ...
718 : ! **************************************************************************************************
719 23488 : 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 23488 : NULLIFY (scf_section, logger)
732 :
733 23488 : logger => cp_get_default_logger()
734 23488 : iounit = cp_logger_get_default_io_unit(logger)
735 23488 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
736 :
737 143954 : energy%repulsive = SUM(tb%e_rep)
738 143954 : energy%el_stat = SUM(tb%e_es)
739 143954 : energy%dispersion = SUM(tb%e_disp)
740 143954 : energy%dispersion_sc = SUM(tb%e_scd)
741 143954 : energy%xtb_xb_inter = SUM(tb%e_hal)
742 143954 : 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 23488 : + energy%efield + energy%qmmm_el
747 :
748 : iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
749 23488 : extension=".scfLog")
750 23488 : IF (iounit > 0) THEN
751 : WRITE (UNIT=iounit, FMT="(/,(T9,A,T60,F20.10))") &
752 153 : "Repulsive pair potential energy: ", energy%repulsive, &
753 153 : "Zeroth order Hamiltonian energy: ", energy%core, &
754 153 : "Electrostatic energy: ", energy%el_stat, &
755 153 : "Self-consistent dispersion energy: ", energy%dispersion_sc, &
756 306 : "Non-self consistent dispersion energy: ", energy%dispersion
757 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
758 153 : "Correction for halogen bonding: ", energy%xtb_xb_inter
759 153 : IF (ABS(energy%efield) > 1.e-12_dp) THEN
760 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
761 0 : "Electric field interaction energy: ", energy%efield
762 : END IF
763 153 : IF (qs_env%qmmm) THEN
764 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
765 0 : "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 23488 : "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 23488 : 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 264 : 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 264 : CALL allocate_gto_basis_set(gto_basis_set)
806 :
807 : !identifying element in the bas data
808 264 : CALL symbol_to_number(i_type, element_symbol)
809 462 : DO id_atom = 1, tb%mol%nat
810 462 : IF (i_type == tb%el_num(id_atom)) EXIT
811 : END DO
812 264 : param%z = i_type
813 264 : param%symbol = element_symbol
814 264 : param%defined = .TRUE.
815 264 : ityp = tb%mol%id(id_atom)
816 :
817 : !getting size information
818 264 : nset = tb%calc%bas%nsh_id(ityp)
819 264 : nshell = 1
820 264 : mprim = 0
821 870 : DO ishell = 1, nset
822 870 : mprim = MAX(mprim, tb%calc%bas%cgto(ishell, ityp)%nprim)
823 : END DO
824 264 : param%nshell = nset
825 264 : natorb = 0
826 :
827 : !write basis set information
828 264 : CALL integer_to_string(mprim, sng)
829 264 : gto_basis_set%name = element_symbol//"_STO-"//TRIM(sng)//"G"
830 264 : gto_basis_set%nset = nset
831 264 : CALL reallocate(gto_basis_set%lmax, 1, nset)
832 264 : CALL reallocate(gto_basis_set%lmin, 1, nset)
833 264 : CALL reallocate(gto_basis_set%npgf, 1, nset)
834 264 : CALL reallocate(gto_basis_set%nshell, 1, nset)
835 264 : CALL reallocate(gto_basis_set%n, 1, 1, 1, nset)
836 264 : CALL reallocate(gto_basis_set%l, 1, 1, 1, nset)
837 264 : CALL reallocate(gto_basis_set%zet, 1, mprim, 1, nset)
838 264 : CALL reallocate(gto_basis_set%gcc, 1, mprim, 1, 1, 1, nset)
839 :
840 264 : ind_ao = 0
841 264 : maxl = 0
842 870 : DO ishell = 1, nset
843 606 : ang = tb%calc%bas%cgto(ishell, ityp)%ang
844 606 : natorb = natorb + (2*ang + 1)
845 606 : param%lval(ishell) = ang
846 606 : maxl = MAX(ang, maxl)
847 606 : gto_basis_set%lmax(ishell) = ang
848 606 : gto_basis_set%lmin(ishell) = ang
849 606 : gto_basis_set%npgf(ishell) = tb%calc%bas%cgto(ishell, ityp)%nprim
850 606 : gto_basis_set%nshell(ishell) = nshell
851 606 : gto_basis_set%n(1, ishell) = ang + 1
852 606 : gto_basis_set%l(1, ishell) = ang
853 3632 : DO ipgf = 1, gto_basis_set%npgf(ishell)
854 3026 : gto_basis_set%gcc(ipgf, 1, ishell) = tb%calc%bas%cgto(ishell, ityp)%coeff(ipgf)
855 3632 : gto_basis_set%zet(ipgf, ishell) = tb%calc%bas%cgto(ishell, ityp)%alpha(ipgf)
856 : END DO
857 2224 : DO ipgf = 1, (2*ang + 1)
858 1354 : ind_ao = ind_ao + 1
859 1354 : param%lao(ind_ao) = ang
860 1960 : param%nao(ind_ao) = ishell
861 : END DO
862 : END DO
863 :
864 264 : do_ortho = .FALSE.
865 264 : CALL process_gto_basis(gto_basis_set, do_ortho, nset, maxl)
866 :
867 : !setting additional values in parameter
868 264 : param%rcut = get_cutoff(tb%calc%bas, tb%accuracy)
869 264 : param%natorb = natorb
870 264 : param%lmax = maxl !max angular momentum
871 :
872 : !getting occupation
873 264 : occ = 0
874 264 : docc = 0.0_dp
875 264 : IF (tb%calc%bas%nsh_at(id_atom) > 5) CPABORT("too many shells in tblite")
876 870 : DO iSh = 1, tb%calc%bas%nsh_at(id_atom)
877 606 : occ(iSh) = NINT(tb%calc%h0%refocc(iSh, ityp) + docc)
878 606 : docc = docc + tb%calc%h0%refocc(iSh, ityp) - REAL(occ(iSh))
879 870 : param%occupation(iSh) = occ(iSh)
880 : END DO
881 264 : IF (ABS(docc) > 0.1_dp) CPABORT("Getting occupation numbers from tblite fails")
882 1584 : param%zeff = SUM(occ) !effective core charge
883 :
884 : !set normalization process
885 264 : 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 264 : END SUBROUTINE tb_get_basis
897 :
898 : ! **************************************************************************************************
899 : !> \brief ...
900 : !> \param qs_env ...
901 : !> \param calculate_forces ...
902 : ! **************************************************************************************************
903 2396 : 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 2396 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
927 2396 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
928 2396 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
929 2396 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint, hint
930 2396 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min
931 2396 : INTEGER, DIMENSION(:), POINTER :: npgfa, npgfb, nsgfa, nsgfb
932 2396 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
933 2396 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
934 2396 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, zeta, zetb, scon_a, scon_b
935 2396 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sblock, fblock
936 : !$ INTEGER :: hash, hash1, lock_num
937 : !$ INTEGER(KIND=int_8) :: iatom8
938 2396 : !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
939 : INTEGER, PARAMETER :: nlock = 501
940 :
941 2396 : 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 2396 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s, matrix_p, matrix_w
946 : TYPE(dft_control_type), POINTER :: dft_control
947 2396 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
948 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
949 2396 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
950 : TYPE(kpoint_type), POINTER :: kpoints
951 2396 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
952 2396 : 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 2396 : 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 2396 : CALL timeset(routineN, handle)
963 :
964 2396 : NULLIFY (ks_env, energy, atomic_kind_set, qs_kind_set)
965 2396 : NULLIFY (matrix_h, matrix_s, atprop, dft_control)
966 2396 : 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 2396 : rho=rho, tb_tblite=tb)
980 2396 : h0 => tb%calc%h0
981 :
982 : !update geometry (required for debug / geometry optimization)
983 2396 : CALL tb_update_geometry(qs_env, tb)
984 :
985 2396 : nkind = SIZE(atomic_kind_set)
986 2396 : nderivatives = 0
987 2396 : IF (calculate_forces) THEN
988 114 : nderivatives = 1
989 114 : IF (ALLOCATED(tb%grad)) DEALLOCATE (tb%grad)
990 342 : ALLOCATE (tb%grad(3, tb%mol%nat))
991 114 : IF (ALLOCATED(tb%dsedcn)) DEALLOCATE (tb%dsedcn)
992 342 : ALLOCATE (tb%dsedcn(tb%calc%bas%nsh))
993 114 : IF (ALLOCATED(tb%calc%ncoord)) THEN
994 114 : IF (ALLOCATED(tb%dcndr)) DEALLOCATE (tb%dcndr)
995 456 : ALLOCATE (tb%dcndr(3, tb%mol%nat, tb%mol%nat))
996 114 : IF (ALLOCATED(tb%dcndL)) DEALLOCATE (tb%dcndL)
997 342 : ALLOCATE (tb%dcndL(3, 3, tb%mol%nat))
998 : END IF
999 : ELSE
1000 2282 : IF (ALLOCATED(tb%grad)) DEALLOCATE (tb%grad)
1001 2282 : IF (ALLOCATED(tb%dcndr)) DEALLOCATE (tb%dcndr)
1002 2282 : IF (ALLOCATED(tb%dcndL)) DEALLOCATE (tb%dcndL)
1003 : END IF
1004 2396 : maxder = ncoset(nderivatives)
1005 2396 : nimg = dft_control%nimages
1006 2396 : IF (nimg > 1) THEN
1007 1046 : IF (.NOT. ASSOCIATED(sab_kp)) CPABORT("Missing k-point neighbor list for tblite")
1008 1046 : sab_orb => sab_kp
1009 1046 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1010 1046 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1011 : END IF
1012 :
1013 : !intialise hamiltonian
1014 2396 : CALL tb_init_ham(qs_env, tb, para_env)
1015 :
1016 : ! get density matrtix
1017 2396 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1018 :
1019 : ! set up matrices for force calculations
1020 2396 : IF (calculate_forces) THEN
1021 114 : NULLIFY (force, matrix_w, virial)
1022 : CALL get_qs_env(qs_env=qs_env, &
1023 : matrix_w_kp=matrix_w, &
1024 114 : virial=virial, force=force)
1025 :
1026 114 : IF (SIZE(matrix_p, 1) == 2) THEN
1027 36 : DO img = 1, nimg
1028 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
1029 18 : 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 36 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1032 : END DO
1033 : END IF
1034 176 : tb%use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1035 : END IF
1036 :
1037 2396 : 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 11036 : ALLOCATE (basis_set_list(nkind))
1041 2396 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
1042 :
1043 : ! allocate overlap matrix
1044 2396 : 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 2396 : sab_orb, .TRUE.)
1047 2396 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
1048 :
1049 : ! initialize H matrix
1050 2396 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
1051 49394 : DO img = 1, nimg
1052 46998 : ALLOCATE (matrix_h(1, img)%matrix)
1053 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, img)%matrix, &
1054 46998 : name="HAMILTONIAN MATRIX")
1055 49394 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
1056 : END DO
1057 2396 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
1058 2396 : 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 2396 : !$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 49394 : DO img = 1, nimg
1265 100968 : DO i = 1, SIZE(matrix_s, 1)
1266 100968 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1267 : END DO
1268 96392 : DO i = 1, SIZE(matrix_h, 1)
1269 93996 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1270 : END DO
1271 : END DO
1272 :
1273 : !compute multipole moments for gfn2
1274 2396 : IF (dft_control%qs_control%xtb_control%tblite_method == gfn2xtb) &
1275 1224 : CALL tb_get_multipole(qs_env, tb)
1276 :
1277 : ! output overlap information
1278 2396 : NULLIFY (logger)
1279 2396 : logger => cp_get_default_logger()
1280 2396 : IF (.NOT. calculate_forces) THEN
1281 2282 : 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 2 : extension=".Log")
1285 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
1286 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
1287 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
1288 2 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
1289 2 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
1290 : END IF
1291 : END IF
1292 :
1293 2396 : DEALLOCATE (basis_set_list)
1294 :
1295 2396 : 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 4792 : 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 49086 : 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 49086 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ch_atom, ch_shell, ch_ref, ch_orb, mix_vars
1345 49086 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, ao_dip, ao_quad
1346 49086 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: aocg_spin, ao_dip_spin, ao_quad_spin, &
1347 49086 : ch_orb_spin, ch_shell_spin
1348 :
1349 49086 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1350 49086 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_p
1351 49086 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
1352 : TYPE(dbcsr_type), POINTER :: s_matrix
1353 49086 : TYPE(error_type), ALLOCATABLE :: error
1354 : TYPE(mp_para_env_type), POINTER :: para_env
1355 49086 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1356 49086 : 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 49086 : 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 49086 : scf_control=scf_control)
1367 :
1368 : ! also compute multipoles needed by GFN2
1369 49086 : do_dipole = .FALSE.
1370 49086 : do_quadrupole = .FALSE.
1371 49086 : 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 49086 : skip_scf_dispersion_energy = skip_scf_dispersion
1377 49086 : skip_scf_dispersion_gradient = skip_scf_dispersion
1378 49086 : 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 49086 : 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 91276 : SELECT CASE (dft_control%qs_control%xtb_control%tblite_scc_mixer)
1392 : CASE (tblite_scc_mixer_auto)
1393 42590 : use_native_mixer = tb_native_scc_mixer_active(dft_control)
1394 42590 : 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 0 : use_no_mixer = .TRUE.
1404 : CASE DEFAULT
1405 48686 : CPABORT("Unknown tblite SCC mixer")
1406 : END SELECT
1407 : END IF
1408 45334 : 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 0 : CPABORT("tblite SCC mixer exceeded TBLITE_MIXER/ITERATIONS")
1411 : END IF
1412 45622 : 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 2180 : dft_control%qs_control%xtb_control%tblite_mixer_solver)
1422 2180 : CALL tb_reset_mixer(tb)
1423 : END IF
1424 49086 : nspin = dft_control%nspins
1425 49086 : IF (nspin /= tb%wfn%nspin) CPABORT("CP2K/tblite spin channel mismatch")
1426 :
1427 49086 : NULLIFY (matrix_p)
1428 49086 : IF (use_rho) THEN
1429 25836 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1430 23250 : ELSEIF (calculate_forces .AND. nspin > 1) THEN
1431 18 : IF (.NOT. ASSOCIATED(tb%rho_ao_kp_ref)) &
1432 0 : CPABORT("Missing converged tblite density for UKS/LSD forces")
1433 18 : matrix_p => tb%rho_ao_kp_ref
1434 : ELSE
1435 23232 : matrix_p => scf_env%p_mix_new
1436 : END IF
1437 49086 : IF (nspin > 1 .AND. (.NOT. calculate_forces)) CALL tb_store_density_ref(tb, matrix_p)
1438 49086 : IF (ASSOCIATED(tb%dipbra)) do_dipole = .TRUE.
1439 49086 : IF (ASSOCIATED(tb%quadbra)) do_quadrupole = .TRUE.
1440 49086 : reuse_native_input = .FALSE.
1441 49086 : IF (use_native_mixer .AND. scf_env%iter_count == 1) THEN
1442 7284 : reuse_native_input = ANY(ABS(tb%wfn%qsh) > 1.0E-14_dp)
1443 4344 : IF (do_dipole) reuse_native_input = reuse_native_input .OR. &
1444 4442 : ANY(ABS(tb%wfn%dpat) > 1.0E-14_dp)
1445 4344 : IF (do_quadrupole) reuse_native_input = reuse_native_input .OR. &
1446 6148 : ANY(ABS(tb%wfn%qpat) > 1.0E-14_dp)
1447 : END IF
1448 49086 : n_atom = SIZE(particle_set)
1449 49086 : nkind = SIZE(atomic_kind_set)
1450 49086 : nimg = dft_control%nimages
1451 49086 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
1452 196344 : ALLOCATE (aocg(nsgf, n_atom))
1453 245430 : ALLOCATE (aocg_spin(nsgf, n_atom, nspin))
1454 49086 : aocg = 0.0_dp
1455 49086 : aocg_spin = 0.0_dp
1456 49086 : IF (do_dipole) THEN
1457 110370 : ALLOCATE (ao_dip(n_atom, dip_n))
1458 147160 : ALLOCATE (ao_dip_spin(n_atom, dip_n, nspin))
1459 36790 : ao_dip_spin = 0.0_dp
1460 : END IF
1461 49086 : IF (do_quadrupole) THEN
1462 110370 : ALLOCATE (ao_quad(n_atom, quad_n))
1463 147160 : ALLOCATE (ao_quad_spin(n_atom, quad_n, nspin))
1464 36790 : ao_quad_spin = 0.0_dp
1465 : END IF
1466 49086 : max_orb = 0
1467 49086 : max_shell = 0
1468 120080 : DO ikind = 1, nkind
1469 70994 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1470 70994 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb)
1471 120080 : max_orb = MAX(max_orb, natorb)
1472 : END DO
1473 299552 : DO is = 1, n_atom
1474 299552 : max_shell = MAX(max_shell, tb%calc%bas%nsh_at(is))
1475 : END DO
1476 343602 : ALLOCATE (ch_atom(n_atom, nspin), ch_shell(n_atom, max_shell))
1477 294516 : ALLOCATE (ch_orb(max_orb, n_atom), ch_ref(max_orb, n_atom))
1478 441774 : ALLOCATE (ch_orb_spin(max_orb, n_atom, nspin), ch_shell_spin(n_atom, max_shell, nspin))
1479 49086 : ch_atom = 0.0_dp
1480 49086 : ch_shell = 0.0_dp
1481 49086 : ch_orb = 0.0_dp
1482 49086 : ch_orb_spin = 0.0_dp
1483 49086 : ch_shell_spin = 0.0_dp
1484 49086 : ch_ref = 0.0_dp
1485 49086 : IF (nimg > 1) THEN
1486 37240 : DO ispin = 1, nspin
1487 18620 : CALL tb_ao_charges_kp_spin(matrix_p, matrix_s, aocg_spin(:, :, ispin), ispin, para_env)
1488 18620 : IF (do_dipole) THEN
1489 57824 : DO im = 1, dip_n
1490 : CALL tb_contract_dens_kp_spin(matrix_p, tb%dipbra, tb%dipket, im, dip_n, &
1491 57824 : ao_dip_spin(:, im, ispin), ispin, para_env)
1492 : END DO
1493 : END IF
1494 37240 : IF (do_quadrupole) THEN
1495 101192 : DO im = 1, quad_n
1496 : CALL tb_contract_dens_kp_spin(matrix_p, tb%quadbra, tb%quadket, im, quad_n, &
1497 101192 : 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 30466 : p_matrix => matrix_p(:, 1)
1504 30466 : s_matrix => matrix_s(1, 1)%matrix
1505 65820 : DO ispin = 1, nspin
1506 35354 : CALL tb_ao_charges_matrix(matrix_p(ispin, 1)%matrix, s_matrix, aocg_spin(:, :, ispin), para_env)
1507 35354 : IF (do_dipole) THEN
1508 104768 : DO im = 1, dip_n
1509 : CALL tb_contract_dens_matrix(matrix_p(ispin, 1)%matrix, tb%dipbra(im)%matrix, &
1510 104768 : tb%dipket(im)%matrix, ao_dip_spin(:, im, ispin), para_env)
1511 : END DO
1512 : END IF
1513 65820 : IF (do_quadrupole) THEN
1514 183344 : DO im = 1, quad_n
1515 : CALL tb_contract_dens_matrix(matrix_p(ispin, 1)%matrix, tb%quadbra(im)%matrix, &
1516 183344 : tb%quadket(im)%matrix, ao_quad_spin(:, im, ispin), para_env)
1517 : END DO
1518 : END IF
1519 : END DO
1520 : END IF
1521 49086 : IF (nspin == 1) THEN
1522 2290856 : aocg(:, :) = aocg_spin(:, :, 1)
1523 688280 : IF (do_dipole) ao_dip(:, :) = ao_dip_spin(:, :, 1)
1524 1332362 : IF (do_quadrupole) ao_quad(:, :) = ao_quad_spin(:, :, 1)
1525 : ELSE
1526 70188 : aocg(:, :) = aocg_spin(:, :, 1) + aocg_spin(:, :, 2)
1527 4888 : IF (do_dipole) THEN
1528 15432 : DO im = 1, dip_n
1529 44304 : DO iatom = 1, n_atom
1530 28872 : pao = ao_dip_spin(iatom, im, 1)
1531 28872 : ao_dip_spin(iatom, im, 1) = pao + ao_dip_spin(iatom, im, 2)
1532 40446 : ao_dip_spin(iatom, im, 2) = pao - ao_dip_spin(iatom, im, 2)
1533 : END DO
1534 : END DO
1535 44304 : ao_dip(:, :) = ao_dip_spin(:, :, 1)
1536 : END IF
1537 4888 : IF (do_quadrupole) THEN
1538 27006 : DO im = 1, quad_n
1539 84750 : DO iatom = 1, n_atom
1540 57744 : pao = ao_quad_spin(iatom, im, 1)
1541 57744 : ao_quad_spin(iatom, im, 1) = pao + ao_quad_spin(iatom, im, 2)
1542 80892 : ao_quad_spin(iatom, im, 2) = pao - ao_quad_spin(iatom, im, 2)
1543 : END DO
1544 : END DO
1545 84750 : ao_quad(:, :) = ao_quad_spin(:, :, 1)
1546 : END IF
1547 : END IF
1548 49086 : NULLIFY (xtb_kind)
1549 120080 : DO ikind = 1, nkind
1550 70994 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
1551 70994 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1552 70994 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, nao=nao, occupation=occ)
1553 441540 : DO iatom = 1, na
1554 250466 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1555 2210332 : DO is = 1, natorb
1556 1959866 : ns = lao(is) + 1
1557 1959866 : norm = 2*lao(is) + 1
1558 1959866 : ch_ref(is, atom_a) = tb%calc%h0%refocc(nao(is), ikind)/norm
1559 1959866 : ch_orb(is, atom_a) = aocg(is, atom_a) - ch_ref(is, atom_a)
1560 1959866 : ch_orb_spin(is, atom_a, 1) = ch_orb(is, atom_a)
1561 1959866 : IF (nspin == 2) ch_orb_spin(is, atom_a, 2) = &
1562 43992 : aocg_spin(is, atom_a, 1) - aocg_spin(is, atom_a, 2)
1563 1959866 : ch_shell(atom_a, ns) = ch_orb(is, atom_a) + ch_shell(atom_a, ns)
1564 4214190 : DO ispin = 1, nspin
1565 : ch_shell_spin(atom_a, ns, ispin) = ch_orb_spin(is, atom_a, ispin) + &
1566 3963724 : ch_shell_spin(atom_a, ns, ispin)
1567 : END DO
1568 : END DO
1569 584910 : DO ispin = 1, nspin
1570 2627724 : ch_atom(atom_a, ispin) = SUM(ch_orb_spin(:, atom_a, ispin))
1571 : END DO
1572 : END DO
1573 : END DO
1574 299552 : native_seed_charge = -SUM(ch_atom(:, 1))
1575 : seed_native_from_rho = SUM(ABS(aocg_spin)) > 1.0E-10_dp .AND. &
1576 2480318 : ABS(native_seed_charge - REAL(dft_control%charge, dp)) < 1.0E-5_dp
1577 49086 : DEALLOCATE (aocg, aocg_spin)
1578 :
1579 49086 : 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 49086 : IF (use_native_mixer) THEN
1611 45334 : IF (.NOT. ALLOCATED(tb%mixer)) CPABORT("tblite mixer not initialized")
1612 45334 : IF (use_rho .AND. (.NOT. calculate_forces) .AND. scf_env%iter_count > 1) THEN
1613 21440 : CALL tb%mixer%next(error)
1614 21440 : IF (ALLOCATED(error)) CPABORT("tblite native mixer failed")
1615 21440 : CALL tb%mixer%get(tb%wfn%qsh)
1616 162654 : tb%wfn%qat(:, :) = 0.0_dp
1617 136004 : DO iatom = 1, n_atom
1618 114564 : ii = tb%calc%bas%ish_at(iatom)
1619 254540 : DO ispin = 1, nspin
1620 : tb%wfn%qat(iatom, ispin) = &
1621 564030 : SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), ispin))
1622 : END DO
1623 : END DO
1624 21440 : IF (do_dipole) THEN
1625 15926 : CALL tb%mixer%get(tb%wfn%dpat)
1626 15926 : DEALLOCATE (ao_dip)
1627 : END IF
1628 21440 : IF (do_quadrupole) THEN
1629 15926 : CALL tb%mixer%get(tb%wfn%qpat)
1630 15926 : DEALLOCATE (ao_quad)
1631 : END IF
1632 : ELSE
1633 23894 : IF (use_rho .AND. (.NOT. calculate_forces)) THEN
1634 2180 : IF (.NOT. reuse_native_input) THEN
1635 130 : 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 610 : DO iatom = 1, n_atom
1639 508 : ii = tb%calc%bas%ish_at(iatom)
1640 1118 : DO ispin = 1, nspin
1641 1844 : DO is = 1, tb%calc%bas%nsh_at(iatom)
1642 1844 : tb%wfn%qsh(ii + is, ispin) = -ch_shell_spin(iatom, is, ispin)
1643 : END DO
1644 : tb%wfn%qat(iatom, ispin) = &
1645 2352 : SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), ispin))
1646 : END DO
1647 : END DO
1648 102 : IF (do_dipole) THEN
1649 300 : DO iatom = 1, n_atom
1650 548 : DO ispin = 1, nspin
1651 1240 : tb%wfn%dpat(:, iatom, ispin) = -ao_dip_spin(iatom, :, ispin)
1652 : END DO
1653 : END DO
1654 : END IF
1655 102 : IF (do_quadrupole) THEN
1656 300 : DO iatom = 1, n_atom
1657 548 : DO ispin = 1, nspin
1658 1984 : tb%wfn%qpat(:, iatom, ispin) = -ao_quad_spin(iatom, :, ispin)
1659 : END DO
1660 : END DO
1661 : END IF
1662 : ELSE
1663 364 : tb%wfn%qsh(:, :) = 0.0_dp
1664 200 : tb%wfn%qat(:, :) = 0.0_dp
1665 204 : IF (do_dipole) tb%wfn%dpat(:, :, :) = 0.0_dp
1666 324 : IF (do_quadrupole) tb%wfn%qpat(:, :, :) = 0.0_dp
1667 : END IF
1668 : END IF
1669 2180 : IF (do_dipole) DEALLOCATE (ao_dip)
1670 2180 : IF (do_quadrupole) DEALLOCATE (ao_quad)
1671 : ELSE
1672 21714 : IF ((.NOT. use_rho) .AND. (.NOT. calculate_forces)) THEN
1673 21608 : CALL tb%mixer%set(tb%wfn%qsh)
1674 21608 : IF (do_dipole) CALL tb%mixer%set(tb%wfn%dpat)
1675 21608 : IF (do_quadrupole) CALL tb%mixer%set(tb%wfn%qpat)
1676 : END IF
1677 137454 : DO iatom = 1, n_atom
1678 115740 : ii = tb%calc%bas%ish_at(iatom)
1679 257234 : DO ispin = 1, nspin
1680 453926 : DO is = 1, tb%calc%bas%nsh_at(iatom)
1681 453926 : tb%wfn%qsh(ii + is, ispin) = -ch_shell_spin(iatom, is, ispin)
1682 : END DO
1683 : tb%wfn%qat(iatom, ispin) = &
1684 569666 : SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), ispin))
1685 : END DO
1686 : END DO
1687 21714 : IF (do_dipole) THEN
1688 104480 : DO iatom = 1, n_atom
1689 195396 : DO ispin = 1, nspin
1690 452084 : tb%wfn%dpat(:, iatom, ispin) = -ao_dip_spin(iatom, :, ispin)
1691 : END DO
1692 : END DO
1693 16060 : DEALLOCATE (ao_dip)
1694 : END IF
1695 21714 : IF (do_quadrupole) THEN
1696 104480 : DO iatom = 1, n_atom
1697 195396 : DO ispin = 1, nspin
1698 724832 : tb%wfn%qpat(:, iatom, ispin) = -ao_quad_spin(iatom, :, ispin)
1699 : END DO
1700 : END DO
1701 16060 : DEALLOCATE (ao_quad)
1702 : END IF
1703 21714 : IF ((.NOT. use_rho) .AND. (.NOT. calculate_forces)) THEN
1704 21608 : CALL tb%mixer%diff(tb%wfn%qsh)
1705 21608 : IF (do_dipole) CALL tb%mixer%diff(tb%wfn%dpat)
1706 21608 : 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 3752 : native_sign_mixing = .FALSE.
1713 3752 : 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 3352 : ELSEIF (nspin > 1) THEN
1718 1938 : n_mix_cols = nspin*max_shell
1719 1938 : IF (do_dipole) n_mix_cols = n_mix_cols + nspin*dip_n
1720 1938 : IF (do_quadrupole) n_mix_cols = n_mix_cols + nspin*quad_n
1721 7752 : ALLOCATE (mix_vars(n_atom, n_mix_cols))
1722 1938 : mix_vars = 0.0_dp
1723 :
1724 1938 : mix_offset = 0
1725 5814 : DO ispin = 1, nspin
1726 : mix_vars(:, mix_offset + 1:mix_offset + max_shell) = &
1727 27132 : -ch_shell_spin(:, 1:max_shell, ispin)
1728 5814 : mix_offset = mix_offset + max_shell
1729 : END DO
1730 1938 : IF (do_dipole) THEN
1731 5814 : DO ispin = 1, nspin
1732 : mix_vars(:, mix_offset + 1:mix_offset + dip_n) = &
1733 38760 : -ao_dip_spin(:, 1:dip_n, ispin)
1734 5814 : mix_offset = mix_offset + dip_n
1735 : END DO
1736 : END IF
1737 1938 : IF (do_quadrupole) THEN
1738 5814 : DO ispin = 1, nspin
1739 : mix_vars(:, mix_offset + 1:mix_offset + quad_n) = &
1740 73644 : -ao_quad_spin(:, 1:quad_n, ispin)
1741 5814 : mix_offset = mix_offset + quad_n
1742 : END DO
1743 : END IF
1744 :
1745 1938 : IF (.NOT. use_no_mixer) THEN
1746 : CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
1747 1938 : mix_vars, para_env, scf_env%iter_count)
1748 : END IF
1749 :
1750 1938 : mix_offset = 0
1751 5814 : DO ispin = 1, nspin
1752 : ch_shell_spin(:, 1:max_shell, ispin) = &
1753 27132 : -mix_vars(:, mix_offset + 1:mix_offset + max_shell)
1754 5814 : mix_offset = mix_offset + max_shell
1755 : END DO
1756 13566 : ch_shell(:, 1:max_shell) = ch_shell_spin(:, 1:max_shell, 1)
1757 1938 : IF (do_dipole) THEN
1758 5814 : DO ispin = 1, nspin
1759 : ao_dip_spin(:, 1:dip_n, ispin) = &
1760 38760 : -mix_vars(:, mix_offset + 1:mix_offset + dip_n)
1761 5814 : mix_offset = mix_offset + dip_n
1762 : END DO
1763 19380 : ao_dip(:, 1:dip_n) = ao_dip_spin(:, 1:dip_n, 1)
1764 : END IF
1765 1938 : IF (do_quadrupole) THEN
1766 5814 : DO ispin = 1, nspin
1767 : ao_quad_spin(:, 1:quad_n, ispin) = &
1768 73644 : -mix_vars(:, mix_offset + 1:mix_offset + quad_n)
1769 5814 : mix_offset = mix_offset + quad_n
1770 : END DO
1771 36822 : ao_quad(:, 1:quad_n) = ao_quad_spin(:, 1:quad_n, 1)
1772 : END IF
1773 1938 : DEALLOCATE (mix_vars)
1774 : ELSE
1775 1414 : do_combined_mixing = do_dipole .OR. do_quadrupole
1776 1414 : native_sign_mixing = do_dipole .OR. do_quadrupole
1777 1414 : discard_mixed_output = .FALSE.
1778 1414 : skip_charge_mixing = use_no_mixer
1779 1414 : IF (skip_charge_mixing) THEN
1780 : !
1781 1414 : ELSEIF (do_combined_mixing) THEN
1782 1414 : n_mix_cols = max_shell
1783 1414 : IF (do_dipole) n_mix_cols = n_mix_cols + dip_n
1784 1414 : IF (do_quadrupole) n_mix_cols = n_mix_cols + quad_n
1785 5656 : ALLOCATE (mix_vars(n_atom, n_mix_cols))
1786 : IF (native_sign_mixing) THEN
1787 15554 : 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 1414 : mix_offset = max_shell
1792 1414 : IF (do_dipole) THEN
1793 : IF (native_sign_mixing) THEN
1794 22624 : 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 1414 : IF (do_quadrupole) THEN
1801 : IF (native_sign_mixing) THEN
1802 43834 : 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 1414 : mix_vars, para_env, scf_env%iter_count)
1809 : IF (.NOT. discard_mixed_output) THEN
1810 : IF (native_sign_mixing) THEN
1811 15554 : 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 1414 : mix_offset = max_shell
1816 1414 : IF (do_dipole) THEN
1817 : IF (native_sign_mixing) THEN
1818 22624 : 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 1414 : IF (do_quadrupole) THEN
1825 : IF (native_sign_mixing) THEN
1826 43834 : 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 1414 : DEALLOCATE (mix_vars)
1833 : ELSE
1834 : CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
1835 0 : ch_shell, para_env, scf_env%iter_count)
1836 : END IF
1837 : END IF
1838 :
1839 : !setting new wave function
1840 3752 : CALL tb%pot%reset
1841 14124 : tb%e_es = 0.0_dp
1842 14124 : tb%e_scd = 0.0_dp
1843 3752 : moment_alpha = scf_env%p_mix_alpha
1844 3752 : IF (nspin > 1) THEN
1845 6540 : DO iatom = 1, n_atom
1846 4360 : ii = tb%calc%bas%ish_at(iatom)
1847 15260 : DO ispin = 1, nspin
1848 26160 : DO is = 1, tb%calc%bas%nsh_at(iatom)
1849 26160 : tb%wfn%qsh(ii + is, ispin) = -ch_shell_spin(iatom, is, ispin)
1850 : END DO
1851 : tb%wfn%qat(iatom, ispin) = &
1852 30520 : SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), ispin))
1853 : END DO
1854 : END DO
1855 2180 : IF (do_dipole) THEN
1856 6540 : DO iatom = 1, n_atom
1857 15260 : DO ispin = 1, nspin
1858 39240 : tb%wfn%dpat(:, iatom, ispin) = -ao_dip_spin(iatom, :, ispin)
1859 : END DO
1860 : END DO
1861 2180 : DEALLOCATE (ao_dip)
1862 : END IF
1863 2180 : IF (do_quadrupole) THEN
1864 6540 : DO iatom = 1, n_atom
1865 15260 : DO ispin = 1, nspin
1866 65400 : tb%wfn%qpat(:, iatom, ispin) = -ao_quad_spin(iatom, :, ispin)
1867 : END DO
1868 : END DO
1869 2180 : DEALLOCATE (ao_quad)
1870 : END IF
1871 : ELSE
1872 7584 : DO iatom = 1, n_atom
1873 6012 : ii = tb%calc%bas%ish_at(iatom)
1874 15288 : DO is = 1, tb%calc%bas%nsh_at(iatom)
1875 9276 : new_charge = -ch_shell(iatom, is)
1876 15288 : tb%wfn%qsh(ii + is, 1) = new_charge
1877 : END DO
1878 7584 : IF (native_sign_mixing) THEN
1879 14140 : tb%wfn%qat(iatom, 1) = SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), 1))
1880 : ELSE
1881 356 : tb%wfn%qat(iatom, 1) = -ch_atom(iatom, 1)
1882 : END IF
1883 : END DO
1884 :
1885 1572 : IF (do_dipole) THEN
1886 7584 : DO iatom = 1, n_atom
1887 25620 : DO im = 1, dip_n
1888 24048 : tb%wfn%dpat(im, iatom, 1) = -ao_dip(iatom, im)
1889 : END DO
1890 : END DO
1891 1572 : DEALLOCATE (ao_dip)
1892 : END IF
1893 1572 : IF (do_quadrupole) THEN
1894 7584 : DO iatom = 1, n_atom
1895 43656 : DO im = 1, quad_n
1896 42084 : tb%wfn%qpat(im, iatom, 1) = -ao_quad(iatom, im)
1897 : END DO
1898 : END DO
1899 1572 : DEALLOCATE (ao_quad)
1900 : END IF
1901 : END IF
1902 : END IF
1903 :
1904 49086 : CALL tb%pot%reset
1905 299552 : tb%e_es = 0.0_dp
1906 299552 : tb%e_scd = 0.0_dp
1907 299552 : tb%e_int = 0.0_dp
1908 49086 : IF (ALLOCATED(tb%calc%coulomb)) THEN
1909 49086 : CALL tb%calc%coulomb%get_potential(tb%mol, tb%cache, tb%wfn, tb%pot)
1910 49086 : CALL tb%calc%coulomb%get_energy(tb%mol, tb%cache, tb%wfn, tb%e_es)
1911 : END IF
1912 49086 : IF (ALLOCATED(tb%calc%dispersion)) THEN
1913 : IF (.NOT. skip_scf_dispersion_potential) &
1914 49086 : CALL tb%calc%dispersion%get_potential(tb%mol, tb%dcache, tb%wfn, tb%pot)
1915 : IF (.NOT. skip_scf_dispersion_energy) &
1916 49086 : CALL tb%calc%dispersion%get_energy(tb%mol, tb%dcache, tb%wfn, tb%e_scd)
1917 : END IF
1918 49086 : IF (ALLOCATED(tb%calc%interactions)) THEN
1919 4888 : CALL tb%calc%interactions%get_potential(tb%mol, tb%icache, tb%wfn, tb%pot)
1920 4888 : CALL tb%calc%interactions%get_energy(tb%mol, tb%icache, tb%wfn, tb%e_int)
1921 : END IF
1922 :
1923 49086 : IF (calculate_forces) THEN
1924 114 : IF (ALLOCATED(tb%calc%coulomb)) THEN
1925 1946 : tb%grad = 0.0_dp
1926 114 : CALL tb%calc%coulomb%get_gradient(tb%mol, tb%cache, tb%wfn, tb%grad, tb%sigma)
1927 114 : CALL tb_dump_sigma_component("after_coulomb", tb%sigma, para_env)
1928 114 : CALL tb_grad2force(qs_env, tb, para_env, 3)
1929 : END IF
1930 :
1931 114 : IF (ALLOCATED(tb%calc%dispersion) .AND. .NOT. skip_scf_dispersion_gradient) THEN
1932 1946 : tb%grad = 0.0_dp
1933 114 : CALL tb%calc%dispersion%get_gradient(tb%mol, tb%dcache, tb%wfn, tb%grad, tb%sigma)
1934 114 : CALL tb_dump_sigma_component("after_dispersion_scf", tb%sigma, para_env)
1935 114 : CALL tb_grad2force(qs_env, tb, para_env, 2)
1936 : END IF
1937 :
1938 114 : IF (ALLOCATED(tb%calc%interactions)) THEN
1939 194 : tb%grad = 0.0_dp
1940 18 : CALL tb%calc%interactions%get_gradient(tb%mol, tb%icache, tb%wfn, tb%grad, tb%sigma)
1941 18 : CALL tb_dump_sigma_component("after_interactions_scf", tb%sigma, para_env)
1942 18 : CALL tb_grad2force(qs_env, tb, para_env, 3)
1943 : END IF
1944 : END IF
1945 :
1946 49086 : IF (ALLOCATED(ao_dip_spin)) DEALLOCATE (ao_dip_spin)
1947 49086 : IF (ALLOCATED(ao_quad_spin)) DEALLOCATE (ao_quad_spin)
1948 49086 : 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 147258 : END SUBROUTINE tb_update_charges
1960 :
1961 : ! **************************************************************************************************
1962 : !> \brief ...
1963 : !> \param qs_env ...
1964 : !> \param tb ...
1965 : !> \param dft_control ...
1966 : ! **************************************************************************************************
1967 25710 : 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 25710 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, sum_shell
1989 25710 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ashift, bshift
1990 25710 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksblock, sblock
1991 25710 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1992 25710 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dip_ket1, dip_ket2, dip_ket3, &
1993 25710 : dip_bra1, dip_bra2, dip_bra3
1994 25710 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: quad_ket1, quad_ket2, quad_ket3, &
1995 25710 : quad_ket4, quad_ket5, quad_ket6, &
1996 25710 : quad_bra1, quad_bra2, quad_bra3, &
1997 25710 : quad_bra4, quad_bra5, quad_bra6
1998 :
1999 25710 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2000 : TYPE(dbcsr_iterator_type) :: iter
2001 25710 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
2002 25710 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
2003 : TYPE(kpoint_type), POINTER :: kpoints
2004 : TYPE(neighbor_list_iterator_p_type), &
2005 25710 : DIMENSION(:), POINTER :: nl_iterator
2006 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2007 25710 : POINTER :: n_list
2008 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2009 25710 : POINTER :: kp_list
2010 25710 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2011 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
2012 :
2013 25710 : nimg = dft_control%nimages
2014 25710 : mpfac = -0.5_dp
2015 :
2016 25710 : 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 25710 : matrix_s_kp=matrix_s, matrix_ks_kp=ks_matrix, qs_kind_set=qs_kind_set)
2019 25710 : IF (nimg > 1) THEN
2020 9798 : IF (.NOT. ASSOCIATED(kp_list)) CPABORT("Missing k-point neighbor list for tblite Hamiltonian")
2021 9798 : n_list => kp_list
2022 : END IF
2023 25710 : nspin = SIZE(ks_matrix, 1)
2024 :
2025 : !creating sum of shell lists
2026 77130 : ALLOCATE (sum_shell(tb%mol%nat))
2027 155974 : i = 0
2028 155974 : DO j = 1, tb%mol%nat
2029 130264 : sum_shell(j) = i
2030 155974 : i = i + tb%calc%bas%nsh_at(j)
2031 : END DO
2032 :
2033 25710 : IF (nimg == 1) THEN
2034 : ! no k-points; all matrices have been transformed to periodic bsf
2035 15912 : 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 15912 : kind_of=kind_of)
2038 15912 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
2039 188401 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2040 172489 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
2041 :
2042 172489 : ikind = kind_of(irow)
2043 172489 : jkind = kind_of(icol)
2044 :
2045 : ! atomic parameters
2046 172489 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
2047 172489 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
2048 172489 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa)
2049 172489 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob)
2050 :
2051 172489 : ni = SIZE(sblock, 1)
2052 689956 : ALLOCATE (ashift(ni, ni))
2053 :
2054 172489 : nj = SIZE(sblock, 2)
2055 689956 : ALLOCATE (bshift(nj, nj))
2056 :
2057 351872 : DO is = 1, nspin
2058 179383 : ashift = 0.0_dp
2059 1569844 : DO i = 1, ni
2060 1390461 : la = naoa(i) + sum_shell(irow)
2061 1569844 : ashift(i, i) = tb_spin_project(tb%pot%vsh(la, :), is)
2062 : END DO
2063 179383 : bshift = 0.0_dp
2064 1533056 : DO j = 1, nj
2065 1353673 : lb = naob(j) + sum_shell(icol)
2066 1533056 : bshift(j, j) = tb_spin_project(tb%pot%vsh(lb, :), is)
2067 : END DO
2068 179383 : NULLIFY (ksblock)
2069 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
2070 179383 : row=irow, col=icol, block=ksblock, found=found)
2071 179383 : CPASSERT(found)
2072 717532 : ksblock = ksblock - 0.5_dp*(MATMUL(ashift, sblock) &
2073 444226831 : + MATMUL(sblock, bshift))
2074 : ksblock = ksblock - 0.5_dp*(tb_spin_project(tb%pot%vat(irow, :), is) &
2075 26302671 : + tb_spin_project(tb%pot%vat(icol, :), is))*sblock
2076 : END DO
2077 533379 : DEALLOCATE (ashift, bshift)
2078 : END DO
2079 15912 : CALL dbcsr_iterator_stop(iter)
2080 :
2081 15912 : IF (ASSOCIATED(tb%dipbra)) THEN
2082 11530 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
2083 144249 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2084 132719 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
2085 :
2086 132719 : NULLIFY (dip_bra1, dip_bra2, dip_bra3)
2087 : CALL dbcsr_get_block_p(matrix=tb%dipbra(1)%matrix, &
2088 132719 : row=irow, col=icol, BLOCK=dip_bra1, found=found)
2089 132719 : CPASSERT(found)
2090 : CALL dbcsr_get_block_p(matrix=tb%dipbra(2)%matrix, &
2091 132719 : row=irow, col=icol, BLOCK=dip_bra2, found=found)
2092 132719 : CPASSERT(found)
2093 : CALL dbcsr_get_block_p(matrix=tb%dipbra(3)%matrix, &
2094 132719 : row=irow, col=icol, BLOCK=dip_bra3, found=found)
2095 132719 : CPASSERT(found)
2096 132719 : NULLIFY (dip_ket1, dip_ket2, dip_ket3)
2097 : CALL dbcsr_get_block_p(matrix=tb%dipket(1)%matrix, &
2098 132719 : row=irow, col=icol, BLOCK=dip_ket1, found=found)
2099 132719 : CPASSERT(found)
2100 : CALL dbcsr_get_block_p(matrix=tb%dipket(2)%matrix, &
2101 132719 : row=irow, col=icol, BLOCK=dip_ket2, found=found)
2102 132719 : CPASSERT(found)
2103 : CALL dbcsr_get_block_p(matrix=tb%dipket(3)%matrix, &
2104 132719 : row=irow, col=icol, BLOCK=dip_ket3, found=found)
2105 132719 : CPASSERT(found)
2106 :
2107 281788 : DO is = 1, nspin
2108 137539 : NULLIFY (ksblock)
2109 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
2110 137539 : row=irow, col=icol, block=ksblock, found=found)
2111 137539 : 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 21994497 : + dip_bra3*tb_spin_project(tb%pot%vdp(3, icol, :), is))
2118 : END DO
2119 : END DO
2120 11530 : CALL dbcsr_iterator_stop(iter)
2121 : END IF
2122 :
2123 15912 : IF (ASSOCIATED(tb%quadbra)) THEN
2124 11530 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
2125 144249 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2126 132719 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
2127 :
2128 132719 : 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 132719 : row=irow, col=icol, BLOCK=quad_bra1, found=found)
2131 132719 : CPASSERT(found)
2132 : CALL dbcsr_get_block_p(matrix=tb%quadbra(2)%matrix, &
2133 132719 : row=irow, col=icol, BLOCK=quad_bra2, found=found)
2134 132719 : CPASSERT(found)
2135 : CALL dbcsr_get_block_p(matrix=tb%quadbra(3)%matrix, &
2136 132719 : row=irow, col=icol, BLOCK=quad_bra3, found=found)
2137 132719 : CPASSERT(found)
2138 : CALL dbcsr_get_block_p(matrix=tb%quadbra(4)%matrix, &
2139 132719 : row=irow, col=icol, BLOCK=quad_bra4, found=found)
2140 132719 : CPASSERT(found)
2141 : CALL dbcsr_get_block_p(matrix=tb%quadbra(5)%matrix, &
2142 132719 : row=irow, col=icol, BLOCK=quad_bra5, found=found)
2143 132719 : CPASSERT(found)
2144 : CALL dbcsr_get_block_p(matrix=tb%quadbra(6)%matrix, &
2145 132719 : row=irow, col=icol, BLOCK=quad_bra6, found=found)
2146 132719 : CPASSERT(found)
2147 :
2148 132719 : 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 132719 : row=irow, col=icol, BLOCK=quad_ket1, found=found)
2151 132719 : CPASSERT(found)
2152 : CALL dbcsr_get_block_p(matrix=tb%quadket(2)%matrix, &
2153 132719 : row=irow, col=icol, BLOCK=quad_ket2, found=found)
2154 132719 : CPASSERT(found)
2155 : CALL dbcsr_get_block_p(matrix=tb%quadket(3)%matrix, &
2156 132719 : row=irow, col=icol, BLOCK=quad_ket3, found=found)
2157 132719 : CPASSERT(found)
2158 : CALL dbcsr_get_block_p(matrix=tb%quadket(4)%matrix, &
2159 132719 : row=irow, col=icol, BLOCK=quad_ket4, found=found)
2160 132719 : CPASSERT(found)
2161 : CALL dbcsr_get_block_p(matrix=tb%quadket(5)%matrix, &
2162 132719 : row=irow, col=icol, BLOCK=quad_ket5, found=found)
2163 132719 : CPASSERT(found)
2164 : CALL dbcsr_get_block_p(matrix=tb%quadket(6)%matrix, &
2165 132719 : row=irow, col=icol, BLOCK=quad_ket6, found=found)
2166 132719 : CPASSERT(found)
2167 :
2168 281788 : DO is = 1, nspin
2169 137539 : NULLIFY (ksblock)
2170 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
2171 137539 : row=irow, col=icol, block=ksblock, found=found)
2172 137539 : 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 21994497 : + quad_bra6*tb_spin_project(tb%pot%vqp(6, icol, :), is))
2186 : END DO
2187 : END DO
2188 11530 : CALL dbcsr_iterator_stop(iter)
2189 : END IF
2190 :
2191 : ELSE
2192 9798 : NULLIFY (kpoints)
2193 9798 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
2194 9798 : NULLIFY (cell_to_index)
2195 9798 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
2196 :
2197 9798 : NULLIFY (nl_iterator)
2198 9798 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
2199 1310935 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2200 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2201 1301137 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
2202 :
2203 1301137 : icol = MAX(iatom, jatom)
2204 1301137 : irow = MIN(iatom, jatom)
2205 :
2206 1301137 : IF (iatom > jatom) THEN
2207 554345 : i = ikind
2208 554345 : ikind = jkind
2209 554345 : jkind = i
2210 : END IF
2211 :
2212 1301137 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2213 1301137 : CPASSERT(ic > 0)
2214 :
2215 1301137 : NULLIFY (sblock)
2216 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
2217 1301137 : row=irow, col=icol, block=sblock, found=found)
2218 1301137 : CPASSERT(found)
2219 :
2220 : ! atomic parameters
2221 1301137 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
2222 1301137 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
2223 1301137 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa)
2224 1301137 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob)
2225 :
2226 1301137 : ni = SIZE(sblock, 1)
2227 5204548 : ALLOCATE (ashift(ni, ni))
2228 :
2229 1301137 : nj = SIZE(sblock, 2)
2230 5204548 : ALLOCATE (bshift(nj, nj))
2231 :
2232 2602274 : DO is = 1, nspin
2233 1301137 : ashift = 0.0_dp
2234 12968215 : DO i = 1, ni
2235 11667078 : la = naoa(i) + sum_shell(irow)
2236 12968215 : ashift(i, i) = tb_spin_project(tb%pot%vsh(la, :), is)
2237 : END DO
2238 1301137 : bshift = 0.0_dp
2239 12968215 : DO j = 1, nj
2240 11667078 : lb = naob(j) + sum_shell(icol)
2241 12968215 : bshift(j, j) = tb_spin_project(tb%pot%vsh(lb, :), is)
2242 : END DO
2243 1301137 : NULLIFY (ksblock)
2244 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
2245 1301137 : row=irow, col=icol, block=ksblock, found=found)
2246 1301137 : CPASSERT(found)
2247 5204548 : ksblock = ksblock - 0.5_dp*(MATMUL(ashift, sblock) &
2248 4127157197 : + MATMUL(sblock, bshift))
2249 : ksblock = ksblock - 0.5_dp*(tb_spin_project(tb%pot%vat(irow, :), is) &
2250 236899731 : + tb_spin_project(tb%pot%vat(icol, :), is))*sblock
2251 : END DO
2252 3913209 : DEALLOCATE (ashift, bshift)
2253 : END DO
2254 9798 : CALL neighbor_list_iterator_release(nl_iterator)
2255 :
2256 9798 : IF (ASSOCIATED(tb%dipbra)) THEN
2257 7486 : NULLIFY (nl_iterator)
2258 7486 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
2259 842837 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2260 : CALL get_iterator_info(nl_iterator, &
2261 835351 : iatom=iatom, jatom=jatom, cell=cellind)
2262 835351 : icol = MAX(iatom, jatom)
2263 835351 : irow = MIN(iatom, jatom)
2264 835351 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2265 835351 : CPASSERT(ic > 0)
2266 835351 : id1 = 1 + dip_n*(ic - 1)
2267 835351 : id2 = 2 + dip_n*(ic - 1)
2268 835351 : id3 = 3 + dip_n*(ic - 1)
2269 :
2270 835351 : NULLIFY (dip_bra1, dip_bra2, dip_bra3)
2271 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id1)%matrix, &
2272 835351 : row=irow, col=icol, BLOCK=dip_bra1, found=found)
2273 835351 : CPASSERT(found)
2274 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id2)%matrix, &
2275 835351 : row=irow, col=icol, BLOCK=dip_bra2, found=found)
2276 835351 : CPASSERT(found)
2277 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id3)%matrix, &
2278 835351 : row=irow, col=icol, BLOCK=dip_bra3, found=found)
2279 835351 : CPASSERT(found)
2280 835351 : NULLIFY (dip_ket1, dip_ket2, dip_ket3)
2281 : CALL dbcsr_get_block_p(matrix=tb%dipket(id1)%matrix, &
2282 835351 : row=irow, col=icol, BLOCK=dip_ket1, found=found)
2283 835351 : CPASSERT(found)
2284 : CALL dbcsr_get_block_p(matrix=tb%dipket(id2)%matrix, &
2285 835351 : row=irow, col=icol, BLOCK=dip_ket2, found=found)
2286 835351 : CPASSERT(found)
2287 : CALL dbcsr_get_block_p(matrix=tb%dipket(id3)%matrix, &
2288 835351 : row=irow, col=icol, BLOCK=dip_ket3, found=found)
2289 835351 : CPASSERT(found)
2290 :
2291 1678188 : DO is = 1, nspin
2292 835351 : NULLIFY (ksblock)
2293 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
2294 835351 : row=irow, col=icol, block=ksblock, found=found)
2295 835351 : 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 151660893 : + dip_bra3*tb_spin_project(tb%pot%vdp(3, icol, :), is))
2302 : END DO
2303 : END DO
2304 7486 : CALL neighbor_list_iterator_release(nl_iterator)
2305 : END IF
2306 :
2307 9798 : IF (ASSOCIATED(tb%quadbra)) THEN
2308 7486 : NULLIFY (nl_iterator)
2309 7486 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
2310 842837 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2311 : CALL get_iterator_info(nl_iterator, &
2312 835351 : iatom=iatom, jatom=jatom, cell=cellind)
2313 835351 : icol = MAX(iatom, jatom)
2314 835351 : irow = MIN(iatom, jatom)
2315 835351 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2316 835351 : CPASSERT(ic > 0)
2317 835351 : iq1 = 1 + quad_n*(ic - 1)
2318 835351 : iq2 = 2 + quad_n*(ic - 1)
2319 835351 : iq3 = 3 + quad_n*(ic - 1)
2320 835351 : iq4 = 4 + quad_n*(ic - 1)
2321 835351 : iq5 = 5 + quad_n*(ic - 1)
2322 835351 : iq6 = 6 + quad_n*(ic - 1)
2323 :
2324 835351 : 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 835351 : row=irow, col=icol, BLOCK=quad_bra1, found=found)
2327 835351 : CPASSERT(found)
2328 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq2)%matrix, &
2329 835351 : row=irow, col=icol, BLOCK=quad_bra2, found=found)
2330 835351 : CPASSERT(found)
2331 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq3)%matrix, &
2332 835351 : row=irow, col=icol, BLOCK=quad_bra3, found=found)
2333 835351 : CPASSERT(found)
2334 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq4)%matrix, &
2335 835351 : row=irow, col=icol, BLOCK=quad_bra4, found=found)
2336 835351 : CPASSERT(found)
2337 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq5)%matrix, &
2338 835351 : row=irow, col=icol, BLOCK=quad_bra5, found=found)
2339 835351 : CPASSERT(found)
2340 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq6)%matrix, &
2341 835351 : row=irow, col=icol, BLOCK=quad_bra6, found=found)
2342 835351 : CPASSERT(found)
2343 :
2344 835351 : 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 835351 : row=irow, col=icol, BLOCK=quad_ket1, found=found)
2347 835351 : CPASSERT(found)
2348 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq2)%matrix, &
2349 835351 : row=irow, col=icol, BLOCK=quad_ket2, found=found)
2350 835351 : CPASSERT(found)
2351 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq3)%matrix, &
2352 835351 : row=irow, col=icol, BLOCK=quad_ket3, found=found)
2353 835351 : CPASSERT(found)
2354 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq4)%matrix, &
2355 835351 : row=irow, col=icol, BLOCK=quad_ket4, found=found)
2356 835351 : CPASSERT(found)
2357 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq5)%matrix, &
2358 835351 : row=irow, col=icol, BLOCK=quad_ket5, found=found)
2359 835351 : CPASSERT(found)
2360 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq6)%matrix, &
2361 835351 : row=irow, col=icol, BLOCK=quad_ket6, found=found)
2362 835351 : CPASSERT(found)
2363 :
2364 1678188 : DO is = 1, nspin
2365 835351 : NULLIFY (ksblock)
2366 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
2367 835351 : row=irow, col=icol, block=ksblock, found=found)
2368 835351 : 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 151660893 : + quad_bra6*tb_spin_project(tb%pot%vqp(6, icol, :), is))
2382 : END DO
2383 : END DO
2384 7486 : 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 51420 : END SUBROUTINE tb_ham_add_coulomb
2397 :
2398 : ! **************************************************************************************************
2399 : !> \brief ...
2400 : !> \param qs_env ...
2401 : !> \param tb ...
2402 : ! **************************************************************************************************
2403 1224 : 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 1224 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2420 : REAL(KIND=dp), DIMENSION(3) :: rij
2421 1224 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max
2422 1224 : INTEGER, DIMENSION(:), POINTER :: nsgfa, nsgfb
2423 1224 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
2424 1224 : INTEGER, ALLOCATABLE :: atom_of_kind(:)
2425 1224 : REAL(KIND=dp), ALLOCATABLE :: stmp(:)
2426 1224 : REAL(KIND=dp), ALLOCATABLE :: dtmp(:, :), qtmp(:, :), dtmpj(:, :), qtmpj(:, :)
2427 1224 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dip_ket1, dip_ket2, dip_ket3, &
2428 1224 : dip_bra1, dip_bra2, dip_bra3
2429 1224 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: quad_ket1, quad_ket2, quad_ket3, &
2430 1224 : quad_ket4, quad_ket5, quad_ket6, &
2431 1224 : quad_bra1, quad_bra2, quad_bra3, &
2432 1224 : quad_bra4, quad_bra5, quad_bra6
2433 :
2434 1224 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2435 1224 : 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 1224 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
2439 : TYPE(kpoint_type), POINTER :: kpoints
2440 1224 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
2441 1224 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_kp
2442 : TYPE(neighbor_list_iterator_p_type), &
2443 1224 : DIMENSION(:), POINTER :: nl_iterator
2444 1224 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2445 1224 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2446 :
2447 1224 : CALL timeset(routineN, handle)
2448 :
2449 : !get info from environment vaiarable
2450 1224 : NULLIFY (atomic_kind_set, qs_kind_set, sab_orb, sab_kp, particle_set)
2451 1224 : 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 1224 : matrix_s_kp=matrix_s)
2460 1224 : natom = SIZE(particle_set)
2461 1224 : nkind = SIZE(atomic_kind_set)
2462 1224 : nimg = dft_control%nimages
2463 1224 : IF (nimg > 1) THEN
2464 552 : IF (.NOT. ASSOCIATED(sab_kp)) CPABORT("Missing k-point neighbor list for tblite multipoles")
2465 552 : sab_orb => sab_kp
2466 552 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
2467 : END IF
2468 :
2469 1224 : 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 5446 : ALLOCATE (basis_set_list(nkind))
2473 1224 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
2474 :
2475 3672 : ALLOCATE (stmp(msao(tb%calc%bas%maxl)**2))
2476 3672 : ALLOCATE (dtmp(dip_n, msao(tb%calc%bas%maxl)**2))
2477 3672 : ALLOCATE (qtmp(quad_n, msao(tb%calc%bas%maxl)**2))
2478 2448 : ALLOCATE (dtmpj(dip_n, msao(tb%calc%bas%maxl)**2))
2479 2448 : ALLOCATE (qtmpj(quad_n, msao(tb%calc%bas%maxl)**2))
2480 :
2481 : ! allocate dipole/quadrupole moment matrix elemnts
2482 1224 : CALL dbcsr_allocate_matrix_set(tb%dipbra, dip_n*nimg)
2483 1224 : CALL dbcsr_allocate_matrix_set(tb%dipket, dip_n*nimg)
2484 1224 : CALL dbcsr_allocate_matrix_set(tb%quadbra, quad_n*nimg)
2485 1224 : CALL dbcsr_allocate_matrix_set(tb%quadket, quad_n*nimg)
2486 21988 : DO img = 1, nimg
2487 83056 : DO i = 1, dip_n
2488 62292 : idx = i + dip_n*(img - 1)
2489 62292 : ALLOCATE (tb%dipbra(idx)%matrix)
2490 62292 : ALLOCATE (tb%dipket(idx)%matrix)
2491 : CALL dbcsr_create(tb%dipbra(idx)%matrix, template=matrix_s(1, img)%matrix, &
2492 62292 : name="DIPOLE BRAMATRIX")
2493 : CALL dbcsr_create(tb%dipket(idx)%matrix, template=matrix_s(1, img)%matrix, &
2494 62292 : name="DIPOLE KETMATRIX")
2495 62292 : CALL cp_dbcsr_alloc_block_from_nbl(tb%dipbra(idx)%matrix, sab_orb)
2496 83056 : CALL cp_dbcsr_alloc_block_from_nbl(tb%dipket(idx)%matrix, sab_orb)
2497 : END DO
2498 146572 : DO i = 1, quad_n
2499 124584 : idx = i + quad_n*(img - 1)
2500 124584 : ALLOCATE (tb%quadbra(idx)%matrix)
2501 124584 : ALLOCATE (tb%quadket(idx)%matrix)
2502 : CALL dbcsr_create(tb%quadbra(idx)%matrix, template=matrix_s(1, img)%matrix, &
2503 124584 : name="QUADRUPOLE BRAMATRIX")
2504 : CALL dbcsr_create(tb%quadket(idx)%matrix, template=matrix_s(1, img)%matrix, &
2505 124584 : name="QUADRUPOLE KETMATRIX")
2506 124584 : CALL cp_dbcsr_alloc_block_from_nbl(tb%quadbra(idx)%matrix, sab_orb)
2507 145348 : 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 1224 : NULLIFY (nl_iterator)
2513 1224 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
2514 254049 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2515 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2516 252825 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
2517 :
2518 1011300 : r2 = NORM2(rij(:))**2
2519 :
2520 252825 : icol = MAX(iatom, jatom)
2521 252825 : irow = MIN(iatom, jatom)
2522 :
2523 252825 : IF (iatom < jatom) THEN
2524 311704 : rij = -rij
2525 77926 : i = ikind
2526 77926 : ikind = jkind
2527 77926 : jkind = i
2528 : END IF
2529 :
2530 252825 : IF (nimg == 1) THEN
2531 : ic = 1
2532 : ELSE
2533 60641 : ic = cell_to_index(cell(1), cell(2), cell(3))
2534 60641 : CPASSERT(ic > 0)
2535 : END IF
2536 252825 : id1 = 1 + dip_n*(ic - 1)
2537 252825 : id2 = 2 + dip_n*(ic - 1)
2538 252825 : id3 = 3 + dip_n*(ic - 1)
2539 252825 : iq1 = 1 + quad_n*(ic - 1)
2540 252825 : iq2 = 2 + quad_n*(ic - 1)
2541 252825 : iq3 = 3 + quad_n*(ic - 1)
2542 252825 : iq4 = 4 + quad_n*(ic - 1)
2543 252825 : iq5 = 5 + quad_n*(ic - 1)
2544 252825 : iq6 = 6 + quad_n*(ic - 1)
2545 :
2546 252825 : ityp = tb%mol%id(icol)
2547 252825 : jtyp = tb%mol%id(irow)
2548 :
2549 252825 : NULLIFY (dip_bra1, dip_bra2, dip_bra3)
2550 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id1)%matrix, &
2551 252825 : row=irow, col=icol, BLOCK=dip_bra1, found=found)
2552 252825 : CPASSERT(found)
2553 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id2)%matrix, &
2554 252825 : row=irow, col=icol, BLOCK=dip_bra2, found=found)
2555 252825 : CPASSERT(found)
2556 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id3)%matrix, &
2557 252825 : row=irow, col=icol, BLOCK=dip_bra3, found=found)
2558 252825 : CPASSERT(found)
2559 :
2560 252825 : NULLIFY (dip_ket1, dip_ket2, dip_ket3)
2561 : CALL dbcsr_get_block_p(matrix=tb%dipket(id1)%matrix, &
2562 252825 : row=irow, col=icol, BLOCK=dip_ket1, found=found)
2563 252825 : CPASSERT(found)
2564 : CALL dbcsr_get_block_p(matrix=tb%dipket(id2)%matrix, &
2565 252825 : row=irow, col=icol, BLOCK=dip_ket2, found=found)
2566 252825 : CPASSERT(found)
2567 : CALL dbcsr_get_block_p(matrix=tb%dipket(id3)%matrix, &
2568 252825 : row=irow, col=icol, BLOCK=dip_ket3, found=found)
2569 252825 : CPASSERT(found)
2570 :
2571 252825 : 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 252825 : row=irow, col=icol, BLOCK=quad_bra1, found=found)
2574 252825 : CPASSERT(found)
2575 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq2)%matrix, &
2576 252825 : row=irow, col=icol, BLOCK=quad_bra2, found=found)
2577 252825 : CPASSERT(found)
2578 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq3)%matrix, &
2579 252825 : row=irow, col=icol, BLOCK=quad_bra3, found=found)
2580 252825 : CPASSERT(found)
2581 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq4)%matrix, &
2582 252825 : row=irow, col=icol, BLOCK=quad_bra4, found=found)
2583 252825 : CPASSERT(found)
2584 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq5)%matrix, &
2585 252825 : row=irow, col=icol, BLOCK=quad_bra5, found=found)
2586 252825 : CPASSERT(found)
2587 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq6)%matrix, &
2588 252825 : row=irow, col=icol, BLOCK=quad_bra6, found=found)
2589 252825 : CPASSERT(found)
2590 :
2591 252825 : 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 252825 : row=irow, col=icol, BLOCK=quad_ket1, found=found)
2594 252825 : CPASSERT(found)
2595 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq2)%matrix, &
2596 252825 : row=irow, col=icol, BLOCK=quad_ket2, found=found)
2597 252825 : CPASSERT(found)
2598 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq3)%matrix, &
2599 252825 : row=irow, col=icol, BLOCK=quad_ket3, found=found)
2600 252825 : CPASSERT(found)
2601 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq4)%matrix, &
2602 252825 : row=irow, col=icol, BLOCK=quad_ket4, found=found)
2603 252825 : CPASSERT(found)
2604 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq5)%matrix, &
2605 252825 : row=irow, col=icol, BLOCK=quad_ket5, found=found)
2606 252825 : CPASSERT(found)
2607 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq6)%matrix, &
2608 252825 : row=irow, col=icol, BLOCK=quad_ket6, found=found)
2609 252825 : CPASSERT(found)
2610 :
2611 : !get basis information
2612 252825 : basis_set_a => basis_set_list(ikind)%gto_basis_set
2613 252825 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
2614 252825 : basis_set_b => basis_set_list(jkind)%gto_basis_set
2615 252825 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
2616 252825 : atom_a = atom_of_kind(icol)
2617 252825 : atom_b = atom_of_kind(irow)
2618 : ! basis a
2619 252825 : first_sgfa => basis_set_a%first_sgf
2620 252825 : la_max => basis_set_a%lmax
2621 252825 : nseta = basis_set_a%nset
2622 252825 : nsgfa => basis_set_a%nsgf_set
2623 : ! basis b
2624 252825 : first_sgfb => basis_set_b%first_sgf
2625 252825 : lb_max => basis_set_b%lmax
2626 252825 : nsetb = basis_set_b%nset
2627 252825 : nsgfb => basis_set_b%nsgf_set
2628 :
2629 : ! --------- Hamiltonian
2630 : ! Periodic self-images are off-diagonal lattice contributions, not the on-site block.
2631 254049 : IF (icol == irow .AND. r2 < same_atom**2) THEN
2632 9729 : DO iset = 1, nseta
2633 29631 : DO jset = 1, nsetb
2634 : CALL multipole_cgto(tb%calc%bas%cgto(jset, ityp), tb%calc%bas%cgto(iset, ityp), &
2635 79608 : & r2, -rij, tb%calc%bas%intcut, stmp, dtmp, qtmp)
2636 :
2637 84440 : DO inda = 1, nsgfa(iset)
2638 57460 : ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
2639 245668 : DO indb = 1, nsgfb(jset)
2640 168306 : ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
2641 168306 : ij = indb + nsgfb(jset)*(inda - 1)
2642 :
2643 168306 : dip_ket1(ib, ia) = dip_ket1(ib, ia) + dtmp(1, ij)
2644 168306 : dip_ket2(ib, ia) = dip_ket2(ib, ia) + dtmp(2, ij)
2645 168306 : dip_ket3(ib, ia) = dip_ket3(ib, ia) + dtmp(3, ij)
2646 :
2647 168306 : quad_ket1(ib, ia) = quad_ket1(ib, ia) + qtmp(1, ij)
2648 168306 : quad_ket2(ib, ia) = quad_ket2(ib, ia) + qtmp(2, ij)
2649 168306 : quad_ket3(ib, ia) = quad_ket3(ib, ia) + qtmp(3, ij)
2650 168306 : quad_ket4(ib, ia) = quad_ket4(ib, ia) + qtmp(4, ij)
2651 168306 : quad_ket5(ib, ia) = quad_ket5(ib, ia) + qtmp(5, ij)
2652 168306 : quad_ket6(ib, ia) = quad_ket6(ib, ia) + qtmp(6, ij)
2653 :
2654 168306 : dip_bra1(ib, ia) = dip_bra1(ib, ia) + dtmp(1, ij)
2655 168306 : dip_bra2(ib, ia) = dip_bra2(ib, ia) + dtmp(2, ij)
2656 168306 : dip_bra3(ib, ia) = dip_bra3(ib, ia) + dtmp(3, ij)
2657 :
2658 168306 : quad_bra1(ib, ia) = quad_bra1(ib, ia) + qtmp(1, ij)
2659 168306 : quad_bra2(ib, ia) = quad_bra2(ib, ia) + qtmp(2, ij)
2660 168306 : quad_bra3(ib, ia) = quad_bra3(ib, ia) + qtmp(3, ij)
2661 168306 : quad_bra4(ib, ia) = quad_bra4(ib, ia) + qtmp(4, ij)
2662 168306 : quad_bra5(ib, ia) = quad_bra5(ib, ia) + qtmp(5, ij)
2663 225766 : 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 998420 : DO iset = 1, nseta
2670 3239764 : DO jset = 1, nsetb
2671 : CALL multipole_cgto(tb%calc%bas%cgto(jset, jtyp), tb%calc%bas%cgto(iset, ityp), &
2672 8965376 : & r2, -rij, tb%calc%bas%intcut, stmp, dtmp, qtmp)
2673 :
2674 9706820 : DO inda = 1, nsgfa(iset)
2675 6717230 : ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
2676 29097556 : DO indb = 1, nsgfb(jset)
2677 20138982 : ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
2678 :
2679 20138982 : ij = indb + nsgfb(jset)*(inda - 1)
2680 : CALL tb_shift_multipole(-rij, stmp(ij), dtmp(:, ij), qtmp(:, ij), &
2681 80555928 : dtmpj(:, ij), qtmpj(:, ij))
2682 :
2683 20138982 : dip_bra1(ib, ia) = dip_bra1(ib, ia) + dtmp(1, ij)
2684 20138982 : dip_bra2(ib, ia) = dip_bra2(ib, ia) + dtmp(2, ij)
2685 20138982 : dip_bra3(ib, ia) = dip_bra3(ib, ia) + dtmp(3, ij)
2686 :
2687 20138982 : quad_bra1(ib, ia) = quad_bra1(ib, ia) + qtmp(1, ij)
2688 20138982 : quad_bra2(ib, ia) = quad_bra2(ib, ia) + qtmp(2, ij)
2689 20138982 : quad_bra3(ib, ia) = quad_bra3(ib, ia) + qtmp(3, ij)
2690 20138982 : quad_bra4(ib, ia) = quad_bra4(ib, ia) + qtmp(4, ij)
2691 20138982 : quad_bra5(ib, ia) = quad_bra5(ib, ia) + qtmp(5, ij)
2692 20138982 : quad_bra6(ib, ia) = quad_bra6(ib, ia) + qtmp(6, ij)
2693 :
2694 20138982 : dip_ket1(ib, ia) = dip_ket1(ib, ia) + dtmpj(1, ij)
2695 20138982 : dip_ket2(ib, ia) = dip_ket2(ib, ia) + dtmpj(2, ij)
2696 20138982 : dip_ket3(ib, ia) = dip_ket3(ib, ia) + dtmpj(3, ij)
2697 :
2698 20138982 : quad_ket1(ib, ia) = quad_ket1(ib, ia) + qtmpj(1, ij)
2699 20138982 : quad_ket2(ib, ia) = quad_ket2(ib, ia) + qtmpj(2, ij)
2700 20138982 : quad_ket3(ib, ia) = quad_ket3(ib, ia) + qtmpj(3, ij)
2701 20138982 : quad_ket4(ib, ia) = quad_ket4(ib, ia) + qtmpj(4, ij)
2702 20138982 : quad_ket5(ib, ia) = quad_ket5(ib, ia) + qtmpj(5, ij)
2703 26856212 : 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 1224 : CALL neighbor_list_iterator_release(nl_iterator)
2711 :
2712 63516 : DO i = 1, SIZE(tb%dipbra)
2713 62292 : CALL dbcsr_finalize(tb%dipbra(i)%matrix)
2714 63516 : CALL dbcsr_finalize(tb%dipket(i)%matrix)
2715 : END DO
2716 125808 : DO i = 1, SIZE(tb%quadbra)
2717 124584 : CALL dbcsr_finalize(tb%quadbra(i)%matrix)
2718 125808 : CALL dbcsr_finalize(tb%quadket(i)%matrix)
2719 : END DO
2720 :
2721 1224 : DEALLOCATE (basis_set_list)
2722 :
2723 1224 : CALL timestop(handle)
2724 :
2725 : #else
2726 : MARK_USED(qs_env)
2727 : MARK_USED(tb)
2728 : CPABORT("Built without TBLITE")
2729 : #endif
2730 :
2731 2448 : 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 20138982 : 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 20138982 : dj(1) = di(1) + vec(1)*s
2752 20138982 : dj(2) = di(2) + vec(2)*s
2753 20138982 : dj(3) = di(3) + vec(3)*s
2754 :
2755 20138982 : qj(1) = 2*vec(1)*di(1) + vec(1)**2*s
2756 20138982 : qj(3) = 2*vec(2)*di(2) + vec(2)**2*s
2757 20138982 : qj(6) = 2*vec(3)*di(3) + vec(3)**2*s
2758 20138982 : qj(2) = vec(1)*di(2) + vec(2)*di(1) + vec(1)*vec(2)*s
2759 20138982 : qj(4) = vec(1)*di(3) + vec(3)*di(1) + vec(1)*vec(3)*s
2760 20138982 : qj(5) = vec(2)*di(3) + vec(3)*di(2) + vec(2)*vec(3)*s
2761 20138982 : tr = 0.5_dp*(qj(1) + qj(3) + qj(6))
2762 :
2763 20138982 : qj(1) = qi(1) + 1.5_dp*qj(1) - tr
2764 20138982 : qj(2) = qi(2) + 1.5_dp*qj(2)
2765 20138982 : qj(3) = qi(3) + 1.5_dp*qj(3) - tr
2766 20138982 : qj(4) = qi(4) + 1.5_dp*qj(4)
2767 20138982 : qj(5) = qi(5) + 1.5_dp*qj(5)
2768 20138982 : qj(6) = qi(6) + 1.5_dp*qj(6) - tr
2769 :
2770 20138982 : 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 751246 : 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 751246 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, s_block
2787 : TYPE(dbcsr_iterator_type) :: iter
2788 :
2789 36445824 : charges = 0.0_dp
2790 751246 : CALL dbcsr_iterator_start(iter, s_matrix)
2791 6709689 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2792 5958443 : NULLIFY (s_block, p_block)
2793 5958443 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, s_block)
2794 5958443 : CALL dbcsr_get_block_p(matrix=p_mat, row=iblock_row, col=iblock_col, BLOCK=p_block, found=found)
2795 5958443 : IF (.NOT. found) CYCLE
2796 5958443 : IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) CYCLE
2797 :
2798 59006235 : DO j = 1, SIZE(p_block, 2)
2799 534867799 : DO i = 1, SIZE(p_block, 1)
2800 528909356 : charges(i, iblock_row) = charges(i, iblock_row) + p_block(i, j)*s_block(i, j)
2801 : END DO
2802 : END DO
2803 6709689 : IF (iblock_col /= iblock_row) THEN
2804 41213921 : DO j = 1, SIZE(p_block, 2)
2805 374104760 : DO i = 1, SIZE(p_block, 1)
2806 369952530 : 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 751246 : CALL dbcsr_iterator_stop(iter)
2812 72140402 : CALL para_env%sum(charges)
2813 :
2814 751246 : 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 18620 : 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 781680 : charges = 0.0_dp
2835 74480 : ALLOCATE (image_charges(SIZE(charges, 1), SIZE(charges, 2)))
2836 734512 : DO ic = 1, SIZE(s_matrix_kp, 2)
2837 : NULLIFY (p_mat, s_mat)
2838 715892 : p_mat => p_matrix_kp(ispin, ic)%matrix
2839 715892 : s_mat => s_matrix_kp(1, ic)%matrix
2840 734512 : IF (ASSOCIATED(p_mat) .AND. ASSOCIATED(s_mat)) THEN
2841 715892 : image_charges = 0.0_dp
2842 715892 : CALL tb_ao_charges_matrix(p_mat, s_mat, image_charges, para_env)
2843 34796272 : charges(:, :) = charges(:, :) + image_charges(:, :)
2844 : END IF
2845 : END DO
2846 18620 : DEALLOCATE (image_charges)
2847 :
2848 18620 : 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 4645152 : 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 4645152 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: bra, ket, p_block
2866 : TYPE(dbcsr_iterator_type) :: iter
2867 :
2868 26433090 : output = 0.0_dp
2869 4645152 : CALL dbcsr_iterator_start(iter, bra_mat)
2870 40250043 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2871 35604891 : NULLIFY (p_block, bra, ket)
2872 35604891 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, bra)
2873 35604891 : CALL dbcsr_get_block_p(matrix=p_mat, row=iblock_row, col=iblock_col, BLOCK=p_block, found=found)
2874 35604891 : IF (.NOT. found) CYCLE
2875 35604891 : CALL dbcsr_get_block_p(matrix=ket_mat, row=iblock_row, col=iblock_col, BLOCK=ket, found=found)
2876 35604891 : IF (.NOT. found) CPABORT("missing block")
2877 :
2878 35604891 : IF (.NOT. (ASSOCIATED(bra) .AND. ASSOCIATED(p_block))) CYCLE
2879 40250043 : IF (iblock_col == iblock_row) THEN
2880 107220069 : DO j = 1, SIZE(p_block, 1)
2881 968330133 : DO i = 1, SIZE(p_block, 2)
2882 957436164 : output(iblock_row) = output(iblock_row) + p_block(j, i)*bra(j, i)
2883 : END DO
2884 : END DO
2885 : ELSE
2886 246069432 : DO j = 1, SIZE(p_block, 1)
2887 2233890450 : DO i = 1, SIZE(p_block, 2)
2888 2209179528 : output(iblock_row) = output(iblock_row) + p_block(j, i)*ket(j, i)
2889 : END DO
2890 : END DO
2891 246069432 : DO j = 1, SIZE(p_block, 1)
2892 2233890450 : DO i = 1, SIZE(p_block, 2)
2893 2209179528 : 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 4645152 : CALL dbcsr_iterator_stop(iter)
2899 48221028 : CALL para_env%sum(output)
2900 :
2901 4645152 : 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 130104 : 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 130104 : nimg = SIZE(p_matrix, 2)
2927 665460 : output = 0.0_dp
2928 390312 : ALLOCATE (image_output(SIZE(output)))
2929 4539528 : DO ic = 1, nimg
2930 4409424 : idx = iop + nops*(ic - 1)
2931 4409424 : CPASSERT(idx <= SIZE(bra_mat))
2932 4409424 : CPASSERT(idx <= SIZE(ket_mat))
2933 : NULLIFY (p_mat)
2934 4409424 : p_mat => p_matrix(ispin, ic)%matrix
2935 4409424 : image_output = 0.0_dp
2936 4409424 : CALL tb_contract_dens_matrix(p_mat, bra_mat(idx)%matrix, ket_mat(idx)%matrix, image_output, para_env)
2937 25053732 : output = output + image_output
2938 : END DO
2939 130104 : DEALLOCATE (image_output)
2940 :
2941 130104 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: bra, ket, p_block
2966 : TYPE(dbcsr_iterator_type) :: iter
2967 :
2968 0 : CALL timeset(routineN, handle)
2969 :
2970 0 : nspin = SIZE(p_matrix)
2971 0 : output = 0.0_dp
2972 0 : DO ispin = 1, nspin
2973 0 : CALL dbcsr_iterator_start(iter, bra_mat)
2974 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2975 0 : NULLIFY (p_block, bra, ket)
2976 :
2977 0 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, bra)
2978 : CALL dbcsr_get_block_p(matrix=p_matrix(ispin)%matrix, &
2979 0 : row=iblock_row, col=iblock_col, BLOCK=p_block, found=found)
2980 0 : IF (.NOT. found) CYCLE
2981 : CALL dbcsr_get_block_p(matrix=ket_mat, &
2982 0 : row=iblock_row, col=iblock_col, BLOCK=ket, found=found)
2983 0 : IF (.NOT. found) CPABORT("missing block")
2984 :
2985 0 : IF (.NOT. (ASSOCIATED(bra) .AND. ASSOCIATED(p_block))) CYCLE
2986 0 : IF (iblock_col == iblock_row) THEN
2987 0 : DO j = 1, SIZE(p_block, 1)
2988 0 : DO i = 1, SIZE(p_block, 2)
2989 0 : output(iblock_row) = output(iblock_row) + p_block(j, i)*bra(j, i)
2990 : END DO
2991 : END DO
2992 : ELSE
2993 0 : DO j = 1, SIZE(p_block, 1)
2994 0 : DO i = 1, SIZE(p_block, 2)
2995 0 : output(iblock_row) = output(iblock_row) + p_block(j, i)*ket(j, i)
2996 : END DO
2997 : END DO
2998 0 : DO j = 1, SIZE(p_block, 1)
2999 0 : DO i = 1, SIZE(p_block, 2)
3000 0 : 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 0 : CALL dbcsr_iterator_stop(iter)
3006 : END DO
3007 :
3008 0 : CALL para_env%sum(output)
3009 0 : CALL timestop(handle)
3010 :
3011 0 : 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 0 : 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 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_image
3033 :
3034 0 : nimg = SIZE(p_matrix, 2)
3035 0 : output = 0.0_dp
3036 0 : ALLOCATE (image_output(SIZE(output)))
3037 0 : DO ic = 1, nimg
3038 0 : idx = iop + nops*(ic - 1)
3039 0 : CPASSERT(idx <= SIZE(bra_mat))
3040 0 : CPASSERT(idx <= SIZE(ket_mat))
3041 : NULLIFY (p_image)
3042 0 : p_image => p_matrix(:, ic)
3043 0 : image_output = 0.0_dp
3044 0 : CALL contract_dens(p_image, bra_mat(idx)%matrix, ket_mat(idx)%matrix, image_output, para_env)
3045 0 : output = output + image_output
3046 : END DO
3047 0 : DEALLOCATE (image_output)
3048 :
3049 0 : 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 754 : 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 754 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
3072 754 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3073 754 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3074 754 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3075 :
3076 754 : CALL timeset(routineN, handle)
3077 :
3078 754 : NULLIFY (force, atomic_kind_set)
3079 : CALL get_qs_env(qs_env=qs_env, force=force, particle_set=particle_set, &
3080 754 : atomic_kind_set=atomic_kind_set)
3081 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
3082 754 : atom_of_kind=atom_of_kind, kind_of=kind_of)
3083 :
3084 754 : natom = SIZE(particle_set)
3085 :
3086 754 : 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 754 : SELECT CASE (ityp)
3101 : CASE DEFAULT
3102 0 : CPABORT("unknown force type")
3103 : CASE (0)
3104 260 : DO iatom = 1, natom
3105 208 : ikind = kind_of(iatom)
3106 208 : atoma = atom_of_kind(iatom)
3107 : force(ikind)%all_potential(:, atoma) = &
3108 884 : force(ikind)%all_potential(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3109 : END DO
3110 : CASE (1)
3111 572 : DO iatom = 1, natom
3112 458 : ikind = kind_of(iatom)
3113 458 : atoma = atom_of_kind(iatom)
3114 : force(ikind)%repulsive(:, atoma) = &
3115 1946 : force(ikind)%repulsive(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3116 : END DO
3117 : CASE (2)
3118 1144 : DO iatom = 1, natom
3119 916 : ikind = kind_of(iatom)
3120 916 : atoma = atom_of_kind(iatom)
3121 : force(ikind)%dispersion(:, atoma) = &
3122 3892 : force(ikind)%dispersion(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3123 : END DO
3124 : CASE (3)
3125 634 : DO iatom = 1, natom
3126 502 : ikind = kind_of(iatom)
3127 502 : atoma = atom_of_kind(iatom)
3128 : force(ikind)%rho_elec(:, atoma) = &
3129 2140 : force(ikind)%rho_elec(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3130 : END DO
3131 : CASE (4)
3132 1144 : DO iatom = 1, natom
3133 916 : ikind = kind_of(iatom)
3134 916 : atoma = atom_of_kind(iatom)
3135 : force(ikind)%overlap(:, atoma) = &
3136 3892 : force(ikind)%overlap(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3137 : END DO
3138 : CASE (5)
3139 754 : DO iatom = 1, natom
3140 0 : ikind = kind_of(iatom)
3141 0 : atoma = atom_of_kind(iatom)
3142 : force(ikind)%efield(:, atoma) = &
3143 0 : force(ikind)%efield(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3144 : END DO
3145 : END SELECT
3146 :
3147 754 : CALL timestop(handle)
3148 :
3149 1508 : END SUBROUTINE tb_grad2force
3150 :
3151 : ! **************************************************************************************************
3152 : !> \brief set gradient to zero
3153 : !> \param qs_env ...
3154 : !> \note
3155 : ! **************************************************************************************************
3156 114 : SUBROUTINE tb_zero_force(qs_env)
3157 :
3158 : TYPE(qs_environment_type) :: qs_env
3159 :
3160 : INTEGER :: iatom, ikind, natom
3161 114 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
3162 114 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3163 114 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3164 114 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3165 :
3166 114 : NULLIFY (force, atomic_kind_set)
3167 : CALL get_qs_env(qs_env=qs_env, force=force, particle_set=particle_set, &
3168 114 : atomic_kind_set=atomic_kind_set)
3169 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
3170 114 : kind_of=kind_of)
3171 :
3172 114 : natom = SIZE(particle_set)
3173 :
3174 572 : DO iatom = 1, natom
3175 458 : ikind = kind_of(iatom)
3176 8162 : force(ikind)%all_potential = 0.0_dp
3177 8162 : force(ikind)%repulsive = 0.0_dp
3178 8162 : force(ikind)%dispersion = 0.0_dp
3179 8162 : force(ikind)%rho_elec = 0.0_dp
3180 8162 : force(ikind)%overlap = 0.0_dp
3181 8276 : force(ikind)%efield = 0.0_dp
3182 : END DO
3183 :
3184 228 : 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 114 : 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 114 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3203 : LOGICAL :: found
3204 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dE
3205 114 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock
3206 114 : 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 114 : DIMENSION(:), POINTER :: nl_iterator
3211 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3212 114 : POINTER :: sab_orb
3213 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3214 114 : 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 114 : 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 114 : sab_orb=sab_orb, sab_kp=sab_kp, para_env=para_env, kpoints=kpoints)
3223 :
3224 114 : NULLIFY (cell_to_index)
3225 114 : IF (nimg > 1) THEN
3226 56 : IF (.NOT. ASSOCIATED(sab_kp)) CPABORT("Missing tblite k-point neighbor list")
3227 56 : sab_orb => sab_kp
3228 56 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
3229 : END IF
3230 :
3231 114 : NULLIFY (matrix_p)
3232 114 : IF (use_rho) THEN
3233 96 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
3234 18 : ELSEIF (ASSOCIATED(tb%rho_ao_kp_ref)) THEN
3235 18 : matrix_p => tb%rho_ao_kp_ref
3236 : ELSE
3237 0 : matrix_p => scf_env%p_mix_new
3238 : END IF
3239 :
3240 342 : ALLOCATE (dE(tb%mol%nat))
3241 :
3242 114 : dE = 0.0_dp
3243 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
3244 114 : NULLIFY (nl_iterator)
3245 114 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
3246 22021 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3247 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3248 21907 : iatom=iatom, jatom=jatom, cell=cellind)
3249 :
3250 21907 : IF (iatom /= jatom) CYCLE
3251 :
3252 5807 : IF (ikind /= jkind) CPABORT("Type wrong")
3253 :
3254 5807 : is = tb%calc%bas%ish_at(iatom)
3255 :
3256 5807 : IF (nimg == 1) THEN
3257 : ic = 1
3258 : ELSE
3259 2468 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
3260 2468 : CPASSERT(ic > 0)
3261 : END IF
3262 :
3263 5807 : IF (cellind(1) /= 0) CYCLE
3264 2017 : IF (cellind(2) /= 0) CYCLE
3265 651 : IF (cellind(3) /= 0) CYCLE
3266 :
3267 594 : DO ispin = 1, SIZE(matrix_p, 1)
3268 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
3269 251 : row=iatom, col=jatom, BLOCK=pblock, found=found)
3270 :
3271 251 : IF (.NOT. found) CPABORT("block not found")
3272 :
3273 251 : ind = 0
3274 23064 : DO ish = 1, tb%calc%bas%nsh_id(ikind)
3275 2643 : DO i = 1, msao(tb%calc%bas%cgto(ish, ikind)%ang)
3276 1737 : ind = ind + 1
3277 2392 : dE(iatom) = dE(iatom) + tb%dsedcn(is + ish)*pblock(ind, ind)
3278 : END DO
3279 : END DO
3280 : END DO
3281 : END DO
3282 114 : CALL neighbor_list_iterator_release(nl_iterator)
3283 114 : CALL para_env%sum(dE)
3284 :
3285 1946 : tb%grad = 0.0_dp
3286 114 : CALL tb_add_grad(tb%grad, tb%dcndr, dE, tb%mol%nat)
3287 114 : IF (tb%use_virial) CALL tb_add_sig(tb%sigma, tb%dcndL, dE, tb%mol%nat)
3288 114 : CALL tb_grad2force(qs_env, tb, para_env, 4)
3289 114 : 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 114 : 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 114 : 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 114 : INTEGER, DIMENSION(:), POINTER :: nsgfa, nsgfb
3322 114 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
3323 114 : 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 114 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dE, t_ov, idip, jdip, idip_mag, jdip_mag, &
3347 114 : iquad, jquad, iquad_mag, jquad_mag
3348 114 : INTEGER, ALLOCATABLE, DIMENSION(:) :: non_nyquist_folded_seen
3349 114 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: t_dip, t_quad, t_d_ov
3350 114 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: t_i_dip, t_i_quad, t_j_dip, t_j_quad
3351 114 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock, pblock_beta, wblock, wblock_beta
3352 :
3353 114 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3354 114 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_w
3355 114 : 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 114 : DIMENSION(:), POINTER :: nl_iterator
3361 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3362 114 : POINTER :: sab_orb
3363 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3364 114 : POINTER :: sab_kp
3365 114 : 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 114 : 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 114 : matrix_w_kp=matrix_w)
3384 :
3385 114 : NULLIFY (cell_to_index)
3386 114 : IF (nimg > 1) THEN
3387 56 : IF (.NOT. ASSOCIATED(sab_kp)) CPABORT("Missing tblite k-point neighbor list")
3388 56 : sab_orb => sab_kp
3389 56 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index, nkp_grid=nkp_grid)
3390 : ELSE
3391 232 : nkp_grid = 1
3392 : END IF
3393 114 : sampled_odd_axes = 0
3394 114 : IF (nkp_grid(1) > 1 .AND. MODULO(nkp_grid(1), 2) == 1) sampled_odd_axes = sampled_odd_axes + 1
3395 114 : IF (nkp_grid(2) > 1 .AND. MODULO(nkp_grid(2), 2) == 1) sampled_odd_axes = sampled_odd_axes + 1
3396 114 : IF (nkp_grid(3) > 1 .AND. MODULO(nkp_grid(3), 2) == 1) sampled_odd_axes = sampled_odd_axes + 1
3397 456 : n_periodic = COUNT(tb%mol%periodic)
3398 : n_sampled_kpoints = 1
3399 456 : DO i = 1, 3
3400 456 : IF (tb%mol%periodic(i)) n_sampled_kpoints = n_sampled_kpoints*nkp_grid(i)
3401 : END DO
3402 :
3403 114 : h0 => tb%calc%h0
3404 114 : has_multipole_response = ASSOCIATED(tb%dipbra) .OR. ASSOCIATED(tb%quadbra)
3405 114 : dip_deriv_fac = -1.0_dp
3406 114 : quad_deriv_fac = -1.0_dp
3407 114 : pot_deriv_scale = 1.0_dp
3408 114 : w_deriv_scale = 1.0_dp
3409 114 : h0_deriv_scale = 1.0_dp
3410 114 : 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 114 : pot_image_deriv_scale = pot_deriv_scale
3422 114 : w_image_deriv_scale = w_deriv_scale
3423 114 : h0_image_deriv_scale = h0_deriv_scale
3424 114 : 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 114 : pot_self_image_deriv_scale = pot_image_deriv_scale
3436 114 : w_self_image_deriv_scale = w_image_deriv_scale
3437 114 : h0_self_image_deriv_scale = h0_image_deriv_scale
3438 114 : 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 304 : IF (nimg > 1 .AND. ANY(nkp_grid > 1) .AND. has_multipole_response) &
3445 30 : 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 114 : cn_image_scale = 1.0_dp
3457 256 : IF (nimg > 1 .AND. ANY(tb%mol%periodic)) cn_image_scale = -0.25_dp
3458 114 : 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 114 : NULLIFY (matrix_p)
3465 114 : IF (use_rho) THEN
3466 96 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
3467 18 : ELSEIF (ASSOCIATED(tb%rho_ao_kp_ref)) THEN
3468 18 : matrix_p => tb%rho_ao_kp_ref
3469 : ELSE
3470 0 : matrix_p => scf_env%p_mix_new
3471 : END IF
3472 114 : nspin = SIZE(matrix_p, 1)
3473 :
3474 : ! set up basis set lists
3475 114 : nkind = SIZE(atomic_kind_set)
3476 516 : ALLOCATE (basis_set_list(nkind))
3477 114 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
3478 :
3479 342 : ALLOCATE (dE(tb%mol%nat))
3480 :
3481 114 : nel = msao(tb%calc%bas%maxl)**2
3482 342 : ALLOCATE (t_ov(nel))
3483 342 : ALLOCATE (t_d_ov(3, nel))
3484 228 : ALLOCATE (t_dip(dip_n, nel))
3485 456 : ALLOCATE (t_i_dip(3, dip_n, nel), t_j_dip(3, dip_n, nel))
3486 342 : ALLOCATE (t_quad(quad_n, nel))
3487 456 : ALLOCATE (t_i_quad(3, quad_n, nel), t_j_quad(3, quad_n, nel))
3488 :
3489 114 : ALLOCATE (idip(dip_n), jdip(dip_n), idip_mag(dip_n), jdip_mag(dip_n))
3490 114 : ALLOCATE (iquad(quad_n), jquad(quad_n), iquad_mag(quad_n), jquad_mag(quad_n))
3491 :
3492 114 : dE = 0.0_dp
3493 1946 : tb%grad = 0.0_dp
3494 114 : hsigma = 0.0_dp
3495 456 : fold_shape = 2*nkp_grid + 1
3496 684 : ALLOCATE (non_nyquist_folded_seen(PRODUCT(fold_shape)))
3497 114 : non_nyquist_folded_seen = 0
3498 114 : has_nyquist_image = 0
3499 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
3500 114 : NULLIFY (nl_iterator)
3501 114 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
3502 22021 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3503 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3504 21907 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
3505 :
3506 21907 : icol = MAX(iatom, jatom)
3507 21907 : jrow = MIN(iatom, jatom)
3508 :
3509 21907 : IF (iatom < jatom) THEN
3510 26676 : rij = -rij
3511 6669 : i = ikind
3512 6669 : ikind = jkind
3513 6669 : jkind = i
3514 : END IF
3515 :
3516 21907 : ityp = tb%mol%id(icol)
3517 21907 : jtyp = tb%mol%id(jrow)
3518 :
3519 87628 : r2 = DOT_PRODUCT(rij, rij)
3520 21907 : dr = SQRT(r2)
3521 21907 : IF (icol == jrow .AND. dr < same_atom) CYCLE
3522 21678 : rr = SQRT(dr/(h0%rad(ikind) + h0%rad(jkind)))
3523 :
3524 : !get basis information
3525 21678 : basis_set_a => basis_set_list(ikind)%gto_basis_set
3526 21678 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
3527 21678 : first_sgfa => basis_set_a%first_sgf
3528 21678 : nsgfa => basis_set_a%nsgf_set
3529 21678 : nseti = basis_set_a%nset
3530 21678 : basis_set_b => basis_set_list(jkind)%gto_basis_set
3531 21678 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
3532 21678 : first_sgfb => basis_set_b%first_sgf
3533 21678 : nsgfb => basis_set_b%nsgf_set
3534 21678 : nsetj = basis_set_b%nset
3535 :
3536 21678 : IF (nimg == 1) THEN
3537 : ic = 1
3538 : ELSE
3539 8854 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
3540 8854 : CPASSERT(ic > 0)
3541 : END IF
3542 21678 : image_pair = (cellind(1) /= 0 .OR. cellind(2) /= 0 .OR. cellind(3) /= 0)
3543 21678 : nkp_cellind = 0
3544 86712 : DO i = 1, 3
3545 86712 : IF (nkp_grid(i) > 1) THEN
3546 25902 : nkp_cellind(i) = MODULO(cellind(i), nkp_grid(i))
3547 25902 : IF (2*nkp_cellind(i) > nkp_grid(i)) nkp_cellind(i) = nkp_cellind(i) - nkp_grid(i)
3548 : END IF
3549 : END DO
3550 21678 : sampled_axes = 0
3551 21678 : IF (nkp_cellind(1) /= 0) sampled_axes = sampled_axes + 1
3552 21678 : IF (nkp_cellind(2) /= 0) sampled_axes = sampled_axes + 1
3553 21678 : IF (nkp_cellind(3) /= 0) sampled_axes = sampled_axes + 1
3554 21678 : kpoint_image_pair = image_pair .AND. sampled_axes > 0
3555 21678 : sampled_even_axes = 0
3556 21678 : IF (nkp_grid(1) > 1 .AND. MODULO(nkp_grid(1), 2) == 0 .AND. &
3557 3669 : ABS(2*nkp_cellind(1)) == nkp_grid(1)) sampled_even_axes = sampled_even_axes + 1
3558 21678 : IF (nkp_grid(2) > 1 .AND. MODULO(nkp_grid(2), 2) == 0 .AND. &
3559 3830 : ABS(2*nkp_cellind(2)) == nkp_grid(2)) sampled_even_axes = sampled_even_axes + 1
3560 21678 : IF (nkp_grid(3) > 1 .AND. MODULO(nkp_grid(3), 2) == 0 .AND. &
3561 3889 : ABS(2*nkp_cellind(3)) == nkp_grid(3)) sampled_even_axes = sampled_even_axes + 1
3562 21678 : sampled_non_nyquist_axes = MAX(0, sampled_axes - sampled_even_axes)
3563 21678 : IF (kpoint_image_pair .AND. sampled_even_axes > 0) has_nyquist_image = 1
3564 21678 : 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 1406 : + fold_shape(1)*fold_shape(2)*(nkp_cellind(3) + nkp_grid(3))
3568 1406 : non_nyquist_folded_seen(fold_index) = 1
3569 : END IF
3570 21678 : 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 21678 : IF (icol == jrow .AND. sampled_image_pair) THEN
3576 2008 : pot_pair_scale = pot_self_image_deriv_scale
3577 2008 : 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 2008 : nyquist_self_image_deriv_scale = 1.0_dp
3582 2008 : IF (has_multipole_response) THEN
3583 712 : IF (n_periodic > 1) &
3584 : nyquist_self_image_deriv_scale = 1.0_dp + REAL(n_periodic, dp)/ &
3585 696 : REAL(n_sampled_kpoints*(1 + dip_n + quad_n), dp)
3586 1296 : ELSEIF (n_periodic == 3) THEN
3587 1216 : nyquist_self_image_deriv_scale = 1.0_dp + 1.0_dp/REAL(n_sampled_kpoints*n_periodic, dp)
3588 : END IF
3589 2008 : 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 2008 : w_pair_scale = w_pair_scale*(1.0_dp - 0.125_dp*REAL(sampled_odd_axes, dp))
3594 2008 : h0_pair_scale = h0_self_image_deriv_scale
3595 2008 : mp_pair_scale = mp_self_image_deriv_scale
3596 : END IF
3597 :
3598 21678 : NULLIFY (pblock, pblock_beta)
3599 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
3600 21678 : row=jrow, col=icol, BLOCK=pblock, found=found)
3601 21678 : IF (.NOT. found) CPABORT("pblock not found")
3602 21678 : IF (nspin > 1) THEN
3603 : CALL dbcsr_get_block_p(matrix=matrix_p(2, ic)%matrix, &
3604 19 : row=jrow, col=icol, BLOCK=pblock_beta, found=found)
3605 19 : IF (.NOT. found) CPABORT("pblock beta not found")
3606 : END IF
3607 :
3608 21678 : NULLIFY (wblock, wblock_beta)
3609 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
3610 21678 : row=jrow, col=icol, block=wblock, found=found)
3611 21678 : CPASSERT(found)
3612 21678 : IF (nspin > 1) THEN
3613 : CALL dbcsr_get_block_p(matrix=matrix_w(2, ic)%matrix, &
3614 19 : row=jrow, col=icol, block=wblock_beta, found=found)
3615 19 : CPASSERT(found)
3616 : END IF
3617 :
3618 21678 : i_a_shift = tb%pot%vat(icol, 1)
3619 21678 : j_a_shift = tb%pot%vat(jrow, 1)
3620 21678 : i_a_shift_mag = 0.0_dp
3621 21678 : j_a_shift_mag = 0.0_dp
3622 21678 : IF (SIZE(tb%pot%vat, 2) > 1) THEN
3623 19 : i_a_shift_mag = tb%pot%vat(icol, 2)
3624 19 : j_a_shift_mag = tb%pot%vat(jrow, 2)
3625 : END IF
3626 86712 : idip(:) = tb%pot%vdp(:, icol, 1)
3627 86712 : jdip(:) = tb%pot%vdp(:, jrow, 1)
3628 21678 : idip_mag(:) = 0.0_dp
3629 21678 : jdip_mag(:) = 0.0_dp
3630 21678 : IF (SIZE(tb%pot%vdp, 3) > 1) THEN
3631 76 : idip_mag(:) = tb%pot%vdp(:, icol, 2)
3632 76 : jdip_mag(:) = tb%pot%vdp(:, jrow, 2)
3633 : END IF
3634 151746 : iquad(:) = tb%pot%vqp(:, icol, 1)
3635 151746 : jquad(:) = tb%pot%vqp(:, jrow, 1)
3636 21678 : iquad_mag(:) = 0.0_dp
3637 21678 : jquad_mag(:) = 0.0_dp
3638 21678 : IF (SIZE(tb%pot%vqp, 3) > 1) THEN
3639 133 : iquad_mag(:) = tb%pot%vqp(:, icol, 2)
3640 133 : jquad_mag(:) = tb%pot%vqp(:, jrow, 2)
3641 : END IF
3642 21678 : dip_deriv_fac_pair = dip_deriv_fac
3643 21678 : quad_deriv_fac_pair = quad_deriv_fac
3644 21678 : overlap_deriv_fac = 1.0_dp
3645 :
3646 21678 : ni = tb%calc%bas%ish_at(icol)
3647 86060 : DO iset = 1, nseti
3648 64268 : sgfi = first_sgfa(1, iset)
3649 64268 : ishift = i_a_shift + tb%pot%vsh(ni + iset, 1)
3650 64268 : ishift_mag = 0.0_dp
3651 64268 : IF (SIZE(tb%pot%vsh, 2) > 1) ishift_mag = i_a_shift_mag + tb%pot%vsh(ni + iset, 2)
3652 64268 : nj = tb%calc%bas%ish_at(jrow)
3653 277517 : DO jset = 1, nsetj
3654 191342 : sgfj = first_sgfb(1, jset)
3655 191342 : jshift = j_a_shift + tb%pot%vsh(nj + jset, 1)
3656 191342 : jshift_mag = 0.0_dp
3657 191342 : 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 191342 : & t_j_dip, t_j_quad)
3663 :
3664 : shpoly = (1.0_dp + h0%shpoly(iset, ikind)*rr) &
3665 191342 : *(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 191342 : *0.5_dp/r2
3669 191342 : scal = h0%hscale(iset, jset, ikind, jkind)
3670 191342 : hij = 0.5_dp*(tb%selfenergy(ni + iset) + tb%selfenergy(nj + jset))*scal
3671 :
3672 191342 : idHdc = tb%dsedcn(ni + iset)*shpoly*scal
3673 191342 : jdHdc = tb%dsedcn(nj + jset)*shpoly*scal
3674 :
3675 191342 : itemp = 0.0_dp
3676 191342 : jtemp = 0.0_dp
3677 191342 : dgrad = 0.0_dp
3678 761422 : DO inda = 1, nsgfa(iset)
3679 570080 : ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
3680 2466192 : DO indb = 1, nsgfb(jset)
3681 1704770 : ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
3682 :
3683 1704770 : ij = inda + nsgfa(iset)*(indb - 1)
3684 :
3685 1704770 : pij_charge = pblock(ib, ia)
3686 1704770 : pij_magnet = 0.0_dp
3687 1704770 : wij_charge = wblock(ib, ia)
3688 1704770 : IF (nspin > 1) THEN
3689 206 : pij_charge = pij_charge + pblock_beta(ib, ia)
3690 206 : 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 1704770 : itemp = itemp + overlap_deriv_fac*idHdc*pij_charge*t_ov(ij)
3696 1704770 : jtemp = jtemp + overlap_deriv_fac*jdHdc*pij_charge*t_ov(ij)
3697 :
3698 1704770 : 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 1704770 : dip_deriv_fac_pair*(MATMUL(t_i_dip(:, :, ij), idip) &
3707 1704770 : + MATMUL(t_j_dip(:, :, ij), jdip)) &
3708 1704770 : + quad_deriv_fac_pair*(MATMUL(t_i_quad(:, :, ij), iquad) &
3709 1704770 : + MATMUL(t_j_quad(:, :, ij), jquad))) &
3710 : + pij_magnet*( &
3711 1704770 : dip_deriv_fac_pair*(MATMUL(t_i_dip(:, :, ij), idip_mag) &
3712 1704770 : + MATMUL(t_j_dip(:, :, ij), jdip_mag)) &
3713 1704770 : + quad_deriv_fac_pair*(MATMUL(t_i_quad(:, :, ij), iquad_mag) &
3714 264809430 : + 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 191342 : IF (icol == jrow) THEN
3720 48552 : IF (sampled_image_pair) THEN
3721 18072 : dE(icol) = dE(icol) + cn_image_scale*0.5_dp*(itemp + jtemp)
3722 : ELSE
3723 30480 : dE(icol) = dE(icol) + 0.5_dp*(itemp + jtemp)
3724 : END IF
3725 : ELSE
3726 142790 : dE(icol) = dE(icol) + itemp
3727 142790 : dE(jrow) = dE(jrow) + jtemp
3728 : END IF
3729 765368 : tb%grad(:, icol) = tb%grad(:, icol) - dgrad
3730 765368 : tb%grad(:, jrow) = tb%grad(:, jrow) + dgrad
3731 255610 : IF (tb%use_virial) THEN
3732 177679 : IF (icol == jrow) THEN
3733 173760 : DO ia = 1, 3
3734 564720 : DO ib = 1, 3
3735 521280 : IF (sampled_image_pair) THEN
3736 116640 : hsigma(ia, ib) = hsigma(ia, ib) - 0.25_dp*(rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
3737 : ELSE
3738 274320 : 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 536956 : DO ia = 1, 3
3744 1745107 : DO ib = 1, 3
3745 1610868 : 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 114 : CALL neighbor_list_iterator_release(nl_iterator)
3754 :
3755 114 : CALL para_env%sum(hsigma)
3756 114 : CALL para_env%sum(dE)
3757 114 : CALL para_env%max(non_nyquist_folded_seen)
3758 114 : CALL para_env%max(has_nyquist_image)
3759 7236 : n_non_nyquist_images = COUNT(non_nyquist_folded_seen > 0)
3760 114 : 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 114 : IF (has_multipole_response .AND. n_periodic == 3 .AND. n_non_nyquist_images > 0) &
3765 2 : cn_deriv_scale = 1.0_dp - 0.5_dp/REAL(n_cn_norm_images, dp)
3766 572 : dE = cn_deriv_scale*dE
3767 114 : CALL para_env%sum(tb%grad)
3768 :
3769 114 : CALL tb_add_grad(tb%grad, tb%dcndr, dE, tb%mol%nat)
3770 114 : CALL tb_grad2force(qs_env, tb, para_env, 4)
3771 :
3772 114 : CALL tb_dump_sigma_component("before_h0_off", tb%sigma, para_env)
3773 1482 : tb%sigma = tb%sigma + hsigma
3774 : CALL tb_dump_sigma_component("after_h0_off_direct", tb%sigma, para_env)
3775 114 : IF (tb%use_virial) THEN
3776 52 : 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 114 : DEALLOCATE (dE)
3781 114 : DEALLOCATE (non_nyquist_folded_seen)
3782 114 : DEALLOCATE (basis_set_list)
3783 114 : DEALLOCATE (t_ov, t_d_ov)
3784 114 : DEALLOCATE (t_dip, t_i_dip, t_j_dip)
3785 114 : DEALLOCATE (t_quad, t_i_quad, t_j_quad)
3786 114 : DEALLOCATE (idip, jdip, idip_mag, jdip_mag, iquad, jquad, iquad_mag, jquad_mag)
3787 :
3788 114 : 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 114 : 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 114 : 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 114 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cli_gradient, cli_virial, cp_gradient
3829 114 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
3830 114 : 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 114 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3837 : TYPE(qs_energy_type), POINTER :: energy
3838 114 : 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 114 : CALL timeset(routineN, handle)
3844 :
3845 114 : NULLIFY (atomic_kind_set, cell, dft_control, energy, force, kpoints, logger, para_env, particle_set, &
3846 114 : 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 114 : virial=virial, do_kpoints=do_kpoints, kpoints=kpoints)
3851 :
3852 114 : ref = dft_control%qs_control%xtb_control%reference_cli
3853 114 : IF (.NOT. ref%enabled) THEN
3854 112 : CALL timestop(handle)
3855 112 : RETURN
3856 : END IF
3857 2 : IF (.NOT. para_env%is_source()) THEN
3858 1 : CALL timestop(handle)
3859 1 : RETURN
3860 : END IF
3861 :
3862 1 : logger => cp_get_default_logger()
3863 1 : iounit = cp_logger_get_default_io_unit(logger)
3864 1 : verbosity_str = ""
3865 1 : IF (logger%iter_info%print_level == silent_print_level) THEN
3866 0 : verbosity_str = " --silent"
3867 1 : ELSE IF (logger%iter_info%print_level == high_print_level .OR. &
3868 : logger%iter_info%print_level == debug_print_level) THEN
3869 0 : verbosity_str = " --verbose"
3870 : END IF
3871 1 : IF (ref%solvation_active) THEN
3872 1 : periodic = 0
3873 1 : IF (ASSOCIATED(cell)) CALL get_cell(cell=cell, periodic=periodic)
3874 4 : n_periodic = COUNT(periodic == 1)
3875 1 : IF (n_periodic == 3) THEN
3876 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
3877 0 : "tblite reference CLI implicit solvation is not supported for PERIODIC XYZ."
3878 : WRITE (UNIT=iounit, FMT="(T2,A)") &
3879 0 : "Use PERIODIC NONE for molecular solvation diagnostics, or remove IMPLICIT_SOLVATION."
3880 0 : CPABORT("REFERENCE_CLI implicit solvation is incompatible with PERIODIC XYZ")
3881 1 : ELSE IF (n_periodic > 0) THEN
3882 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
3883 0 : "WARNING: tblite reference CLI implicit solvation with finite periodicity is diagnostic only."
3884 : WRITE (UNIT=iounit, FMT="(T2,A,I0,A)") &
3885 0 : "The generated native tblite reference geometry has ", n_periodic, &
3886 0 : " periodic direction(s); continuum-solvation conventions are primarily molecular."
3887 : END IF
3888 : END IF
3889 1 : unsupported_kpoints = .FALSE.
3890 1 : IF (do_kpoints .AND. ASSOCIATED(kpoints)) THEN
3891 0 : nkp = 0
3892 0 : CALL get_kpoint_info(kpoint=kpoints, nkp=nkp, xkp=xkp)
3893 0 : unsupported_kpoints = nkp > 1
3894 0 : IF (nkp == 1 .AND. ASSOCIATED(xkp)) unsupported_kpoints = ANY(ABS(xkp(:, 1)) > 1.0E-12_dp)
3895 : END IF
3896 0 : IF (unsupported_kpoints) THEN
3897 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
3898 0 : "tblite reference CLI check skipped: CP2K KPOINTS are active."
3899 : WRITE (UNIT=iounit, FMT="(T2,A)") &
3900 0 : "The native tblite CLI reference path does not reproduce CP2K multi-k-point sampling."
3901 0 : IF (ref%stop_on_error) CPABORT("tblite reference CLI cannot check CP2K k-point calculations")
3902 0 : CALL timestop(handle)
3903 0 : RETURN
3904 : END IF
3905 1 : IF (dft_control%qs_control%xtb_control%tblite_scc_mixer == tblite_scc_mixer_cp2k) THEN
3906 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
3907 0 : "WARNING: tblite reference CLI cannot reproduce XTB/SCC_MIXER CP2K."
3908 : WRITE (UNIT=iounit, FMT="(T2,A)") &
3909 0 : "The external native tblite run uses tblite's own SCC mixer; only the converged result is compared."
3910 0 : IF (ref%stop_on_error) CPABORT("tblite reference CLI cannot reproduce SCC_MIXER CP2K")
3911 : END IF
3912 1 : natom = SIZE(particle_set)
3913 1 : method = tb_reference_method_name(dft_control%qs_control%xtb_control%tblite_method)
3914 1 : guess = tb_reference_guess_name(ref%guess)
3915 1 : solver = tb_reference_solver_name(dft_control%qs_control%xtb_control%tblite_mixer_solver)
3916 1 : file_base = tb_join_path(ref%work_directory, ref%prefix)
3917 1 : gen_file = TRIM(file_base)//".gen"
3918 1 : grad_file = TRIM(file_base)//".grad"
3919 1 : json_file = TRIM(file_base)//".json"
3920 1 : log_file = TRIM(file_base)//".log"
3921 1 : post_processing_output_file = ""
3922 1 : IF (LEN_TRIM(ref%grad_file) > 0) grad_file = ref%grad_file
3923 1 : IF (LEN_TRIM(ref%json_file) > 0) json_file = ref%json_file
3924 1 : IF (LEN_TRIM(ref%post_processing_output_file) > 0) &
3925 1 : post_processing_output_file = ref%post_processing_output_file
3926 :
3927 1 : WRITE (charge_str, "(I0)") dft_control%charge
3928 1 : spin = MAX(0, dft_control%multiplicity - 1)
3929 1 : WRITE (spin_str, "(I0)") spin
3930 1 : WRITE (acc_str, "(ES16.8)") dft_control%qs_control%xtb_control%tblite_accuracy
3931 1 : reference_iterations = dft_control%qs_control%xtb_control%tblite_mixer_iterations
3932 1 : WRITE (iter_str, "(I0)") reference_iterations
3933 1 : efield_str = ""
3934 1 : IF (ref%efield_active) THEN
3935 1 : WRITE (efield_x_str, "(ES16.8)") ref%efield(1)
3936 1 : WRITE (efield_y_str, "(ES16.8)") ref%efield(2)
3937 1 : WRITE (efield_z_str, "(ES16.8)") ref%efield(3)
3938 : efield_str = " --efield "//TRIM(ADJUSTL(efield_x_str))//","// &
3939 1 : TRIM(ADJUSTL(efield_y_str))//","//TRIM(ADJUSTL(efield_z_str))
3940 : END IF
3941 1 : solvation_str = ""
3942 1 : solvation_model_name = ""
3943 1 : IF (ref%solvation_active) THEN
3944 2 : SELECT CASE (ref%solvation_model)
3945 : CASE (tblite_cli_solvation_alpb)
3946 1 : solvation_model_name = "ALPB"
3947 1 : solvation_str = " --alpb "
3948 : CASE (tblite_cli_solvation_gbsa)
3949 0 : solvation_model_name = "GBSA"
3950 0 : solvation_str = " --gbsa "
3951 : CASE (tblite_cli_solvation_gbe)
3952 0 : solvation_model_name = "GBE"
3953 0 : solvation_str = " --gbe "
3954 : CASE (tblite_cli_solvation_gb)
3955 0 : solvation_model_name = "GB"
3956 0 : solvation_str = " --gb "
3957 : CASE (tblite_cli_solvation_cpcm)
3958 0 : solvation_model_name = "CPCM"
3959 0 : solvation_str = " --cpcm "
3960 : CASE DEFAULT
3961 1 : CPABORT("Unknown tblite reference CLI implicit-solvation model")
3962 : END SELECT
3963 1 : solvation_str = TRIM(solvation_str)//" "//TRIM(tb_shell_quote(ref%solvation_solvent))
3964 2 : SELECT CASE (ref%solvation_born_kernel)
3965 : CASE (tblite_cli_born_kernel_auto)
3966 : CASE (tblite_cli_born_kernel_p16)
3967 1 : solvation_str = TRIM(solvation_str)//" --born-kernel p16"
3968 : CASE (tblite_cli_born_kernel_still)
3969 0 : solvation_str = TRIM(solvation_str)//" --born-kernel still"
3970 : CASE DEFAULT
3971 1 : CPABORT("Unknown tblite reference CLI Born kernel")
3972 : END SELECT
3973 2 : SELECT CASE (ref%solvation_state)
3974 : CASE (tblite_cli_solution_state_gsolv)
3975 : CASE (tblite_cli_solution_state_bar1mol)
3976 1 : solvation_str = TRIM(solvation_str)//" --solv-state bar1mol"
3977 : CASE (tblite_cli_solution_state_reference)
3978 0 : solvation_str = TRIM(solvation_str)//" --solv-state reference"
3979 : CASE DEFAULT
3980 1 : CPABORT("Unknown tblite reference CLI solution state")
3981 : END SELECT
3982 : END IF
3983 1 : IF (dft_control%qs_control%xtb_control%tblite_mixer_memory /= reference_iterations) THEN
3984 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
3985 0 : "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/MEMORY."
3986 : WRITE (UNIT=iounit, FMT="(T2,A,I0,A,I0,A)") &
3987 0 : "The native reference run uses tblite's internal mixer memory tied to --iterations (", &
3988 0 : reference_iterations, "), while CP2K uses MEMORY ", &
3989 0 : dft_control%qs_control%xtb_control%tblite_mixer_memory, "."
3990 0 : IF (ref%stop_on_error) CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/MEMORY")
3991 : END IF
3992 1 : 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 0 : "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/DAMPING."
3996 : WRITE (UNIT=iounit, FMT="(T2,A,F8.4,A,F8.4,A)") &
3997 0 : "The native reference run uses tblite's library default ", tblite_mixer_damping_default, &
3998 0 : ", while CP2K uses DAMPING ", dft_control%qs_control%xtb_control%tblite_mixer_damping, "."
3999 0 : IF (ref%stop_on_error) CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/DAMPING")
4000 : END IF
4001 1 : 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 0 : "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/OMEGA0."
4005 : WRITE (UNIT=iounit, FMT="(T2,A,ES12.4,A,ES12.4,A)") &
4006 0 : "The native reference run uses tblite's library default ", tblite_mixer_omega0_default, &
4007 0 : ", while CP2K uses OMEGA0 ", dft_control%qs_control%xtb_control%tblite_mixer_omega0, "."
4008 0 : IF (ref%stop_on_error) CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/OMEGA0")
4009 : END IF
4010 1 : 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 0 : "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 0 : "The native reference run uses tblite's library default ", tblite_mixer_min_weight_default, &
4016 0 : ", while CP2K uses MIN_WEIGHT ", dft_control%qs_control%xtb_control%tblite_mixer_min_weight, "."
4017 0 : IF (ref%stop_on_error) &
4018 0 : CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/MIN_WEIGHT")
4019 : END IF
4020 1 : 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 0 : "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 0 : "The native reference run uses tblite's library default ", tblite_mixer_max_weight_default, &
4026 0 : ", while CP2K uses MAX_WEIGHT ", dft_control%qs_control%xtb_control%tblite_mixer_max_weight, "."
4027 0 : IF (ref%stop_on_error) &
4028 0 : CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/MAX_WEIGHT")
4029 : END IF
4030 1 : 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 0 : "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 0 : "The native reference run uses tblite's library default ", tblite_mixer_weight_factor_default, &
4036 0 : ", while CP2K uses WEIGHT_FACTOR ", &
4037 0 : dft_control%qs_control%xtb_control%tblite_mixer_weight_factor, "."
4038 0 : IF (ref%stop_on_error) &
4039 0 : CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/WEIGHT_FACTOR")
4040 : END IF
4041 1 : etemp = 300.0_dp
4042 1 : IF (ASSOCIATED(scf_control) .AND. ASSOCIATED(scf_control%smear)) THEN
4043 1 : IF (scf_control%smear%do_smear) THEN
4044 1 : etemp = cp_unit_from_cp2k(scf_control%smear%electronic_temperature, "K")
4045 1 : IF (scf_control%smear%method /= smear_fermi_dirac) THEN
4046 : WRITE (UNIT=iounit, FMT="(/,T2,A,A,A)") &
4047 0 : "WARNING: tblite reference CLI cannot reproduce CP2K smearing method ", &
4048 0 : TRIM(tb_reference_smear_method_name(scf_control%smear%method)), "."
4049 : WRITE (UNIT=iounit, FMT="(T2,A,F12.3,A)") &
4050 0 : "The native reference run uses Fermi-Dirac electronic temperature ", etemp, " K instead."
4051 : END IF
4052 : END IF
4053 : END IF
4054 1 : WRITE (etemp_str, "(ES16.8)") etemp
4055 1 : etemp_guess = 0.0_dp
4056 1 : etemp_guess_str = ""
4057 1 : IF (ref%electronic_temperature_guess > 0.0_dp) THEN
4058 0 : etemp_guess = cp_unit_from_cp2k(ref%electronic_temperature_guess, "K")
4059 0 : WRITE (etemp_guess_val_str, "(ES16.8)") etemp_guess
4060 0 : etemp_guess_str = " --etemp-guess "//TRIM(ADJUSTL(etemp_guess_val_str))
4061 : END IF
4062 1 : param_str = ""
4063 1 : IF (LEN_TRIM(dft_control%qs_control%xtb_control%tblite_param_file) > 0) &
4064 0 : param_str = " --param "//TRIM(tb_shell_quote(dft_control%qs_control%xtb_control%tblite_param_file))
4065 1 : spinpol_str = ""
4066 1 : IF (dft_control%lsd) spinpol_str = " --spin-polarized"
4067 1 : post_processing_str = ""
4068 1 : IF (LEN_TRIM(ref%post_processing) > 0) &
4069 1 : post_processing_str = " --post-processing "//TRIM(tb_shell_quote(ref%post_processing))
4070 1 : post_processing_output_str = ""
4071 1 : IF (LEN_TRIM(post_processing_output_file) > 0) THEN
4072 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
4073 1 : "WARNING: tblite reference CLI POST_PROCESSING_OUTPUT was requested explicitly."
4074 : WRITE (UNIT=iounit, FMT="(T2,A)") &
4075 1 : "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 1 : TRIM(tb_shell_quote(post_processing_output_file))
4078 : END IF
4079 1 : restart_str = " --no-restart"
4080 1 : IF (LEN_TRIM(ref%restart_file) > 0) &
4081 0 : restart_str = " --restart "//TRIM(tb_shell_quote(ref%restart_file))
4082 :
4083 1 : 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 1 : " > "//TRIM(tb_shell_quote(log_file))//" 2>&1"
4105 :
4106 1 : cmdstat = 0
4107 1 : exitstat = 0
4108 1 : CALL execute_command_line(TRIM(command), exitstat=exitstat, cmdstat=cmdstat)
4109 1 : IF (cmdstat /= 0 .OR. exitstat /= 0) THEN
4110 1 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "tblite reference CLI check failed to run."
4111 1 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Command: ", TRIM(command)
4112 1 : WRITE (UNIT=iounit, FMT="(T2,A,I0,T32,A,I0)") "cmdstat:", cmdstat, "exitstat:", exitstat
4113 1 : IF (ref%stop_on_error) CPABORT("tblite reference CLI command failed")
4114 1 : CALL tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file, post_processing_output_file)
4115 1 : CALL timestop(handle)
4116 1 : RETURN
4117 : END IF
4118 :
4119 0 : 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 0 : have_energy, have_gradient, have_virial)
4122 :
4123 0 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "tblite reference CLI check"
4124 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Executable: ", TRIM(ref%program_name)
4125 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Method: ", TRIM(method)
4126 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Guess: ", TRIM(guess)
4127 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Solver: ", TRIM(solver)
4128 0 : WRITE (UNIT=iounit, FMT="(T2,A,L1)") "Spin-pol.: ", dft_control%lsd
4129 0 : IF (ref%efield_active) &
4130 0 : WRITE (UNIT=iounit, FMT="(T2,A,3ES16.8,A)") "Efield: ", ref%efield, " V/Angstrom"
4131 0 : IF (ref%solvation_active) &
4132 0 : WRITE (UNIT=iounit, FMT="(T2,A,A,1X,A)") "Solvation: ", TRIM(solvation_model_name), &
4133 0 : TRIM(ref%solvation_solvent)
4134 0 : IF (LEN_TRIM(ref%post_processing) > 0) &
4135 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Post proc.: ", TRIM(ref%post_processing)
4136 0 : IF (LEN_TRIM(post_processing_output_file) > 0) &
4137 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "PP output: ", TRIM(post_processing_output_file)
4138 0 : IF (ref%electronic_temperature_guess > 0.0_dp) &
4139 0 : WRITE (UNIT=iounit, FMT="(T2,A,F12.3,A)") "Guess etemp:", etemp_guess, " K"
4140 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Grad file: ", TRIM(grad_file)
4141 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "JSON file: ", TRIM(json_file)
4142 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Log file: ", TRIM(log_file)
4143 :
4144 0 : too_large = .FALSE.
4145 0 : IF (ref%check_energy) THEN
4146 0 : IF (have_energy) THEN
4147 0 : cp_energy = energy%total
4148 0 : ediff = ABS(cp_energy - cli_energy)
4149 : WRITE (UNIT=iounit, FMT="(T2,A,3ES22.12)") &
4150 0 : "Energy CP2K/CLI/absdiff:", cp_energy, cli_energy, ediff
4151 0 : too_large = too_large .OR. ediff > ref%error_limit
4152 : ELSE
4153 0 : WRITE (UNIT=iounit, FMT="(T2,A)") "Energy check skipped: no CLI energy found."
4154 : END IF
4155 : END IF
4156 :
4157 0 : IF (ref%check_forces) THEN
4158 0 : IF (have_gradient .AND. ASSOCIATED(force)) THEN
4159 0 : ALLOCATE (cp_gradient(3, natom))
4160 0 : CALL total_qs_force(cp_gradient, force, atomic_kind_set)
4161 0 : fsum = SUM(ABS(cp_gradient - cli_gradient))
4162 0 : fmax = MAXVAL(ABS(cp_gradient - cli_gradient))
4163 0 : WRITE (UNIT=iounit, FMT="(T2,A,2ES22.12)") "Gradient diff sum/max:", fsum, fmax
4164 0 : too_large = too_large .OR. fmax > ref%error_limit
4165 0 : DEALLOCATE (cp_gradient)
4166 : ELSE
4167 0 : WRITE (UNIT=iounit, FMT="(T2,A)") "Gradient check skipped: no CLI gradient or CP2K force found."
4168 : END IF
4169 : END IF
4170 :
4171 0 : IF (ref%check_virial) THEN
4172 0 : 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 0 : vsum = SUM(ABS(-virial%pv_virial - cli_virial))
4176 0 : vmax = MAXVAL(ABS(-virial%pv_virial - cli_virial))
4177 0 : WRITE (UNIT=iounit, FMT="(T2,A,2ES22.12)") "Virial diff sum/max:", vsum, vmax
4178 0 : too_large = too_large .OR. vmax > ref%error_limit
4179 : ELSE
4180 0 : WRITE (UNIT=iounit, FMT="(T2,A)") "Virial check skipped: no CLI virial or CP2K virial found."
4181 : END IF
4182 : END IF
4183 :
4184 0 : IF (too_large) THEN
4185 : WRITE (UNIT=iounit, FMT="(T2,A,ES12.4)") &
4186 0 : "tblite reference CLI deviation exceeded ERROR_LIMIT = ", ref%error_limit
4187 0 : IF (ref%stop_on_error) CPABORT("tblite reference CLI deviation exceeded ERROR_LIMIT")
4188 : END IF
4189 :
4190 0 : CALL tb_reference_cli_aux_commands(ref, dft_control, gen_file, file_base, verbosity_str, iounit)
4191 :
4192 0 : CALL tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file, post_processing_output_file)
4193 0 : DEALLOCATE (cli_gradient, cli_virial)
4194 :
4195 0 : CALL timestop(handle)
4196 :
4197 912 : 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 1 : FUNCTION tb_reference_method_name(method_id) RESULT(method)
4205 : INTEGER, INTENT(IN) :: method_id
4206 : CHARACTER(LEN=8) :: method
4207 :
4208 1 : SELECT CASE (method_id)
4209 : CASE (gfn1xtb)
4210 0 : method = "gfn1"
4211 : CASE (gfn2xtb)
4212 1 : method = "gfn2"
4213 : CASE (ipea1xtb)
4214 0 : method = "ipea1"
4215 : CASE DEFAULT
4216 1 : CPABORT("Unknown tblite reference CLI method")
4217 : END SELECT
4218 :
4219 1 : 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 1 : FUNCTION tb_reference_guess_name(guess_id) RESULT(guess)
4227 : INTEGER, INTENT(IN) :: guess_id
4228 : CHARACTER(LEN=8) :: guess
4229 :
4230 2 : SELECT CASE (guess_id)
4231 : CASE (tblite_guess_sad)
4232 1 : guess = "sad"
4233 : CASE (tblite_guess_eeq)
4234 0 : guess = "eeq"
4235 : CASE (tblite_guess_ceh)
4236 0 : guess = "ceh"
4237 : CASE DEFAULT
4238 1 : CPABORT("Unknown tblite reference CLI guess")
4239 : END SELECT
4240 :
4241 1 : 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 1 : FUNCTION tb_reference_solver_name(solver_id) RESULT(solver)
4249 : INTEGER, INTENT(IN) :: solver_id
4250 : CHARACTER(LEN=8) :: solver
4251 :
4252 2 : SELECT CASE (solver_id)
4253 : CASE (tblite_solver_gvd)
4254 1 : solver = "gvd"
4255 : CASE (tblite_solver_gvr)
4256 0 : solver = "gvr"
4257 : CASE DEFAULT
4258 1 : CPABORT("Unknown tblite reference CLI solver")
4259 : END SELECT
4260 :
4261 1 : 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 0 : FUNCTION tb_reference_smear_method_name(method_id) RESULT(method)
4269 : INTEGER, INTENT(IN) :: method_id
4270 : CHARACTER(LEN=24) :: method
4271 :
4272 0 : SELECT CASE (method_id)
4273 : CASE (smear_fermi_dirac)
4274 0 : method = "FERMI_DIRAC"
4275 : CASE (smear_energy_window)
4276 0 : method = "ENERGY_WINDOW"
4277 : CASE (smear_list)
4278 0 : method = "LIST"
4279 : CASE (smear_gaussian)
4280 0 : method = "GAUSSIAN"
4281 : CASE (smear_mp)
4282 0 : method = "METHFESSEL_PAXTON"
4283 : CASE (smear_mv)
4284 0 : method = "MARZARI_VANDERBILT"
4285 : CASE DEFAULT
4286 0 : method = "UNKNOWN"
4287 : END SELECT
4288 :
4289 0 : 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 0 : 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 0 : WRITE (charge_str, "(I0)") dft_control%charge
4319 0 : spin = MAX(0, dft_control%multiplicity - 1)
4320 0 : WRITE (spin_str, "(I0)") spin
4321 :
4322 0 : IF (ref%guess_cli%enabled) THEN
4323 0 : guess_input = ref%guess_cli%input_file
4324 0 : IF (LEN_TRIM(guess_input) == 0) guess_input = gen_file
4325 0 : etemp_guess_str = ""
4326 0 : IF (ref%guess_cli%electronic_temperature_guess > 0.0_dp) THEN
4327 0 : etemp_guess = cp_unit_from_cp2k(ref%guess_cli%electronic_temperature_guess, "K")
4328 0 : WRITE (etemp_guess_val_str, "(ES16.8)") etemp_guess
4329 0 : etemp_guess_str = " --etemp-guess "//TRIM(ADJUSTL(etemp_guess_val_str))
4330 : END IF
4331 0 : efield_str = ""
4332 0 : IF (ref%guess_cli%efield_active) THEN
4333 0 : WRITE (efield_x_str, "(ES16.8)") ref%guess_cli%efield(1)
4334 0 : WRITE (efield_y_str, "(ES16.8)") ref%guess_cli%efield(2)
4335 0 : WRITE (efield_z_str, "(ES16.8)") ref%guess_cli%efield(3)
4336 : efield_str = " --efield "//TRIM(ADJUSTL(efield_x_str))//","// &
4337 0 : TRIM(ADJUSTL(efield_y_str))//","//TRIM(ADJUSTL(efield_z_str))
4338 : END IF
4339 0 : grad_str = ""
4340 0 : IF (ref%guess_cli%grad) grad_str = " --grad"
4341 0 : json_str = ""
4342 0 : IF (LEN_TRIM(ref%guess_cli%json_file) > 0) &
4343 0 : json_str = " --json "//TRIM(tb_shell_quote(ref%guess_cli%json_file))
4344 0 : 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 0 : " > "//TRIM(tb_shell_quote(log_file))//" 2>&1"
4357 0 : CALL tb_reference_cli_execute(ref, "guess", command, log_file, iounit)
4358 0 : IF (.NOT. ref%keep_files .AND. LEN_TRIM(ref%guess_cli%json_file) > 0) &
4359 0 : CALL tb_delete_file(ref%guess_cli%json_file)
4360 : END IF
4361 :
4362 0 : IF (ref%param_cli%enabled) THEN
4363 0 : method_str = ""
4364 0 : 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 0 : ref%param_cli%method_explicit)))
4369 0 : output_str = ""
4370 0 : IF (LEN_TRIM(ref%param_cli%output_file) > 0) &
4371 0 : output_str = " --output "//TRIM(tb_shell_quote(ref%param_cli%output_file))
4372 0 : log_file = TRIM(file_base)//".param.log"
4373 : command = TRIM(tb_shell_quote(ref%program_name))//" param"// &
4374 : TRIM(method_str)// &
4375 0 : TRIM(output_str)
4376 0 : IF (LEN_TRIM(ref%param_cli%input_file) > 0) &
4377 0 : command = TRIM(command)//" "//TRIM(tb_shell_quote(ref%param_cli%input_file))
4378 0 : command = TRIM(command)//" > "//TRIM(tb_shell_quote(log_file))//" 2>&1"
4379 0 : CALL tb_reference_cli_execute(ref, "param", command, log_file, iounit)
4380 0 : IF (.NOT. ref%keep_files .AND. LEN_TRIM(ref%param_cli%output_file) > 0) &
4381 0 : CALL tb_delete_file(ref%param_cli%output_file)
4382 : END IF
4383 :
4384 0 : IF (ref%fit_cli%enabled) THEN
4385 0 : dry_run_str = ""
4386 0 : IF (ref%fit_cli%dry_run) dry_run_str = " --dry-run"
4387 0 : copy_str = ""
4388 0 : IF (LEN_TRIM(ref%fit_cli%copy_file) > 0) &
4389 0 : copy_str = " --copy "//TRIM(tb_shell_quote(ref%fit_cli%copy_file))
4390 0 : 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 0 : " > "//TRIM(tb_shell_quote(log_file))//" 2>&1"
4397 0 : CALL tb_reference_cli_execute(ref, "fit", command, log_file, iounit)
4398 0 : IF (.NOT. ref%keep_files .AND. LEN_TRIM(ref%fit_cli%copy_file) > 0) &
4399 0 : CALL tb_delete_file(ref%fit_cli%copy_file)
4400 : END IF
4401 :
4402 0 : IF (ref%tagdiff_cli%enabled) THEN
4403 0 : method_str = ""
4404 0 : IF (ref%tagdiff_cli%fit) method_str = " --fit"
4405 0 : 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 0 : " > "//TRIM(tb_shell_quote(log_file))//" 2>&1"
4410 0 : CALL tb_reference_cli_execute(ref, "tagdiff", command, log_file, iounit)
4411 : END IF
4412 :
4413 0 : 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 0 : 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 0 : cmdstat = 0
4432 0 : exitstat = 0
4433 0 : CALL execute_command_line(TRIM(command), exitstat=exitstat, cmdstat=cmdstat)
4434 0 : IF (cmdstat /= 0 .OR. exitstat /= 0) THEN
4435 0 : WRITE (UNIT=iounit, FMT="(/,T2,A,A)") "tblite reference CLI auxiliary command failed: ", &
4436 0 : TRIM(label)
4437 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Command: ", TRIM(command)
4438 0 : WRITE (UNIT=iounit, FMT="(T2,A,I0,T32,A,I0)") "cmdstat:", cmdstat, "exitstat:", exitstat
4439 0 : IF (ref%stop_on_error) CPABORT("tblite reference CLI auxiliary command failed")
4440 : ELSE
4441 0 : WRITE (UNIT=iounit, FMT="(/,T2,A,A)") "tblite reference CLI auxiliary command completed: ", &
4442 0 : TRIM(label)
4443 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Log file: ", TRIM(log_file)
4444 : END IF
4445 0 : IF (.NOT. ref%keep_files) CALL tb_delete_file(log_file)
4446 :
4447 0 : END SUBROUTINE tb_reference_cli_execute
4448 :
4449 : ! **************************************************************************************************
4450 : !> \brief Join directory and filename.
4451 : !> \param directory ...
4452 : !> \param filename ...
4453 : !> \return ...
4454 : ! **************************************************************************************************
4455 1 : FUNCTION tb_join_path(directory, filename) RESULT(path)
4456 : CHARACTER(LEN=*), INTENT(IN) :: directory, filename
4457 : CHARACTER(LEN=default_path_length) :: path
4458 :
4459 1 : IF (LEN_TRIM(directory) == 0 .OR. TRIM(directory) == ".") THEN
4460 1 : path = TRIM(filename)
4461 0 : ELSE IF (directory(LEN_TRIM(directory):LEN_TRIM(directory)) == "/") THEN
4462 0 : path = TRIM(directory)//TRIM(filename)
4463 : ELSE
4464 0 : path = TRIM(directory)//"/"//TRIM(filename)
4465 : END IF
4466 :
4467 1 : END FUNCTION tb_join_path
4468 :
4469 : ! **************************************************************************************************
4470 : !> \brief Shell-quote a filename or executable path.
4471 : !> \param text ...
4472 : !> \return ...
4473 : ! **************************************************************************************************
4474 9 : 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 9 : quoted = "'"
4481 150 : DO i = 1, LEN_TRIM(text)
4482 150 : IF (text(i:i) == "'") THEN
4483 0 : quoted = TRIM(quoted)//"'\\''"
4484 : ELSE
4485 141 : quoted = TRIM(quoted)//text(i:i)
4486 : END IF
4487 : END DO
4488 9 : quoted = TRIM(quoted)//"'"
4489 :
4490 9 : 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 1 : 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 1 : CHARACTER(LEN=2), ALLOCATABLE, DIMENSION(:) :: symbols, unique_symbols
4503 : INTEGER :: iatom, ikind, ios, natom, nuniq, unit_nr
4504 1 : 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 1 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
4510 :
4511 1 : NULLIFY (cell, particle_set)
4512 1 : CALL get_qs_env(qs_env=qs_env, cell=cell, particle_set=particle_set)
4513 :
4514 1 : natom = SIZE(particle_set)
4515 1 : to_angstrom = cp_unit_from_cp2k(1.0_dp, "angstrom")
4516 5 : ALLOCATE (symbols(natom), unique_symbols(natom), species(natom))
4517 1 : nuniq = 0
4518 5 : DO iatom = 1, natom
4519 4 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=symbols(iatom))
4520 4 : found = .FALSE.
4521 9 : DO ikind = 1, nuniq
4522 9 : IF (TRIM(unique_symbols(ikind)) == TRIM(symbols(iatom))) THEN
4523 : found = .TRUE.
4524 : EXIT
4525 : END IF
4526 : END DO
4527 4 : IF (.NOT. found) THEN
4528 3 : nuniq = nuniq + 1
4529 3 : unique_symbols(nuniq) = symbols(iatom)
4530 3 : ikind = nuniq
4531 : END IF
4532 5 : species(iatom) = ikind
4533 : END DO
4534 :
4535 : OPEN (NEWUNIT=unit_nr, FILE=TRIM(filename), STATUS="REPLACE", ACTION="WRITE", &
4536 1 : FORM="FORMATTED", IOSTAT=ios)
4537 1 : IF (ios /= 0) CPABORT("Could not open tblite reference CLI geometry file")
4538 :
4539 1 : CALL get_cell(cell=cell, periodic=periodic)
4540 4 : IF (ANY(periodic == 1)) THEN
4541 0 : WRITE (UNIT=unit_nr, FMT="(I0,1X,A)") natom, "S"
4542 : ELSE
4543 1 : WRITE (UNIT=unit_nr, FMT="(I0,1X,A)") natom, "C"
4544 : END IF
4545 4 : WRITE (UNIT=unit_nr, FMT="(*(A,1X))") (TRIM(unique_symbols(ikind)), ikind=1, nuniq)
4546 5 : DO iatom = 1, natom
4547 : WRITE (UNIT=unit_nr, FMT="(I0,1X,I0,3(1X,ES24.16))") &
4548 17 : iatom, species(iatom), particle_set(iatom)%r(:)*to_angstrom
4549 : END DO
4550 4 : IF (ANY(periodic == 1)) THEN
4551 0 : WRITE (UNIT=unit_nr, FMT="(3(1X,ES24.16))") 0.0_dp, 0.0_dp, 0.0_dp
4552 0 : DO ikind = 1, 3
4553 0 : WRITE (UNIT=unit_nr, FMT="(3(1X,ES24.16))") cell%hmat(:, ikind)*to_angstrom
4554 : END DO
4555 : END IF
4556 1 : CLOSE (unit_nr)
4557 :
4558 1 : DEALLOCATE (symbols, unique_symbols, species)
4559 :
4560 1 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: values
4586 :
4587 0 : have_energy = .FALSE.
4588 0 : have_gradient = .FALSE.
4589 0 : have_virial = .FALSE.
4590 0 : energy = 0.0_dp
4591 0 : gradient = 0.0_dp
4592 0 : virial = 0.0_dp
4593 :
4594 0 : INQUIRE (FILE=TRIM(filename), EXIST=exists)
4595 0 : IF (.NOT. exists) RETURN
4596 :
4597 : OPEN (NEWUNIT=unit_nr, FILE=TRIM(filename), STATUS="OLD", ACTION="READ", &
4598 0 : FORM="FORMATTED", IOSTAT=ios)
4599 0 : IF (ios /= 0) RETURN
4600 :
4601 : DO
4602 0 : READ (UNIT=unit_nr, FMT="(A)", IOSTAT=ios) line
4603 0 : IF (ios /= 0) EXIT
4604 0 : IF (INDEX(line, "energy :real:0:") > 0) THEN
4605 0 : READ (UNIT=unit_nr, FMT="(A)", IOSTAT=ios) line
4606 0 : IF (ios == 0) THEN
4607 0 : READ (line, *, IOSTAT=ios) energy
4608 0 : have_energy = ios == 0
4609 : END IF
4610 0 : ELSE IF (INDEX(line, "gradient :real:2:3,") > 0) THEN
4611 0 : ALLOCATE (values(3*natom))
4612 0 : CALL tb_read_real_values(unit_nr, values, nread)
4613 0 : IF (nread == 3*natom) THEN
4614 0 : CALL tb_values_to_matrix(values, gradient)
4615 0 : have_gradient = .TRUE.
4616 : END IF
4617 0 : DEALLOCATE (values)
4618 0 : ELSE IF (INDEX(line, "virial :real:2:3,3") > 0) THEN
4619 0 : ALLOCATE (values(9))
4620 0 : CALL tb_read_real_values(unit_nr, values, nread)
4621 0 : IF (nread == 9) THEN
4622 0 : CALL tb_values_to_matrix(values, virial)
4623 0 : have_virial = .TRUE.
4624 : END IF
4625 0 : DEALLOCATE (values)
4626 : END IF
4627 : END DO
4628 0 : CLOSE (unit_nr)
4629 :
4630 0 : 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 0 : 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 0 : nread = 0
4648 0 : DO WHILE (nread < SIZE(values))
4649 0 : READ (UNIT=unit_nr, FMT="(A)", IOSTAT=ios) line
4650 0 : IF (ios /= 0) EXIT
4651 0 : CALL tb_parse_real_line(line, values, nread)
4652 : END DO
4653 :
4654 0 : 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 0 : 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 0 : pos = 1
4672 0 : DO WHILE (pos <= LEN_TRIM(line) .AND. nread < SIZE(values))
4673 0 : DO WHILE (pos <= LEN_TRIM(line) .AND. INDEX(" ,[]", line(pos:pos)) > 0)
4674 0 : pos = pos + 1
4675 : END DO
4676 0 : IF (pos > LEN_TRIM(line)) EXIT
4677 : first = pos
4678 0 : DO WHILE (pos <= LEN_TRIM(line) .AND. INDEX(" ,[]", line(pos:pos)) == 0)
4679 0 : pos = pos + 1
4680 : END DO
4681 0 : last = pos - 1
4682 0 : token = line(first:last)
4683 0 : READ (token, *, IOSTAT=ios) values(nread + 1)
4684 0 : IF (ios == 0) nread = nread + 1
4685 : END DO
4686 :
4687 0 : 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 0 : 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 0 : n = 0
4702 0 : DO j = 1, SIZE(matrix, 2)
4703 0 : DO i = 1, SIZE(matrix, 1)
4704 0 : n = n + 1
4705 0 : matrix(i, j) = values(n)
4706 : END DO
4707 : END DO
4708 :
4709 0 : 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 1 : 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 1 : IF (ref%keep_files) RETURN
4727 1 : CALL tb_delete_file(gen_file)
4728 1 : CALL tb_delete_file(grad_file)
4729 1 : CALL tb_delete_file(json_file)
4730 1 : CALL tb_delete_file(log_file)
4731 1 : IF (LEN_TRIM(post_processing_output_file) > 0) &
4732 1 : 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 5 : SUBROUTINE tb_delete_file(filename)
4741 :
4742 : CHARACTER(LEN=*), INTENT(IN) :: filename
4743 :
4744 : INTEGER :: ios, unit_nr
4745 : LOGICAL :: exists
4746 :
4747 5 : INQUIRE (FILE=TRIM(filename), EXIST=exists)
4748 5 : IF (.NOT. exists) RETURN
4749 2 : OPEN (NEWUNIT=unit_nr, FILE=TRIM(filename), STATUS="OLD", IOSTAT=ios)
4750 2 : 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 0 : 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 0 : 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 52 : 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 52 : NULLIFY (virial, cell)
4804 52 : CALL get_qs_env(qs_env=qs_env, virial=virial, cell=cell)
4805 52 : CALL get_cell(cell=cell, periodic=periodic)
4806 :
4807 72 : 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 4 : "but it is not a physically meaningful bulk stress for an isolated molecule.")
4812 : END IF
4813 :
4814 52 : 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 676 : virial%pv_virial = virial%pv_virial - tb%sigma/para_env%num_pe
4829 :
4830 52 : 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 228 : 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 1144 : DO i = 1, natom
4849 5628 : DO j = 1, natom
4850 18852 : grad(:, i) = grad(:, i) + deriv(:, i, j)*dE(j)
4851 : END DO
4852 : END DO
4853 :
4854 228 : 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 104 : 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 416 : DO i = 1, 3
4873 1940 : DO j = 1, natom
4874 6408 : sig(:, i) = sig(:, i) + deriv(:, i, j)*dE(j)
4875 : END DO
4876 : END DO
4877 :
4878 104 : END SUBROUTINE tb_add_sig
4879 :
4880 3185290 : END MODULE tblite_interface
|
