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