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