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_cutoff, ONLY: get_lattice_points
24 : USE tblite_integral_multipole, ONLY: multipole_cgto, multipole_grad_cgto, maxl, msao
25 : USE tblite_integral_type, ONLY: integral_type, new_integral
26 : USE tblite_scf, ONLY: get_mixer_dimension, new_mixer
27 : USE tblite_scf_info, ONLY: scf_info, atom_resolved, shell_resolved, &
28 : orbital_resolved, not_used
29 : USE tblite_scf_potential, ONLY: potential_type, new_potential, add_pot_to_h1
30 : USE tblite_wavefunction_type, ONLY: wavefunction_type, new_wavefunction
31 : USE tblite_xtb_calculator, ONLY: xtb_calculator, new_xtb_calculator
32 : USE tblite_xtb_gfn1, ONLY: new_gfn1_calculator
33 : USE tblite_xtb_gfn2, ONLY: new_gfn2_calculator
34 : USE tblite_xtb_h0, ONLY: get_selfenergy, get_hamiltonian, get_occupation, &
35 : get_hamiltonian_gradient, tb_hamiltonian
36 : USE tblite_xtb_ipea1, ONLY: new_ipea1_calculator
37 : #endif
38 : USE ai_contraction, ONLY: block_add, &
39 : contraction
40 : USE ai_overlap, ONLY: overlap_ab
41 : USE atomic_kind_types, ONLY: atomic_kind_type, get_atomic_kind, get_atomic_kind_set
42 : USE atprop_types, ONLY: atprop_type
43 : USE basis_set_types, ONLY: gto_basis_set_type, gto_basis_set_p_type, &
44 : & allocate_gto_basis_set, write_gto_basis_set, process_gto_basis
45 : USE cell_types, ONLY: cell_type, get_cell
46 : USE cp_blacs_env, ONLY: cp_blacs_env_type
47 : USE cp_control_types, ONLY: dft_control_type, xtb_reference_cli_type
48 : USE cp_dbcsr_api, ONLY: dbcsr_type, dbcsr_p_type, dbcsr_create, dbcsr_add, &
49 : dbcsr_get_block_p, dbcsr_finalize, &
50 : dbcsr_iterator_type, dbcsr_iterator_blocks_left, &
51 : dbcsr_iterator_start, dbcsr_iterator_stop, &
52 : dbcsr_iterator_next_block
53 : USE cp_dbcsr_contrib, ONLY: dbcsr_print
54 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
55 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
56 : USE cp_log_handling, ONLY: cp_get_default_logger, &
57 : cp_logger_type, cp_logger_get_default_io_unit
58 : USE cp_output_handling, ONLY: cp_print_key_should_output, &
59 : cp_print_key_unit_nr, cp_print_key_finished_output
60 : USE cp_units, ONLY: cp_unit_from_cp2k
61 : USE input_constants, ONLY: gfn1xtb, gfn2xtb, ipea1xtb, &
62 : tblite_scc_mixer_auto, tblite_scc_mixer_cp2k, &
63 : tblite_scc_mixer_none, tblite_scc_mixer_tblite
64 : USE input_section_types, ONLY: section_vals_val_get
65 : USE kinds, ONLY: default_path_length, default_string_length, dp
66 : USE kpoint_types, ONLY: get_kpoint_info, kpoint_type
67 : USE memory_utilities, ONLY: reallocate
68 : USE message_passing, ONLY: mp_para_env_type
69 : USE mulliken, ONLY: ao_charges
70 : USE orbital_pointers, ONLY: ncoset
71 : USE particle_types, ONLY: particle_type
72 : USE qs_charge_mixing, ONLY: charge_mixing
73 : USE qs_condnum, ONLY: overlap_condnum
74 : USE qs_energy_types, ONLY: qs_energy_type
75 : USE qs_environment_types, ONLY: get_qs_env, qs_environment_type
76 : USE qs_force_types, ONLY: qs_force_type, total_qs_force
77 : USE qs_integral_utils, ONLY: basis_set_list_setup, get_memory_usage
78 : USE qs_kind_types, ONLY: get_qs_kind, qs_kind_type, get_qs_kind_set
79 : USE qs_ks_types, ONLY: get_ks_env, qs_ks_env_type, set_ks_env
80 : USE qs_neighbor_list_types, ONLY: neighbor_list_iterator_create, neighbor_list_iterate, &
81 : get_iterator_info, neighbor_list_set_p_type, &
82 : neighbor_list_iterator_p_type, neighbor_list_iterator_release
83 : USE qs_overlap, ONLY: create_sab_matrix
84 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
85 : USE qs_rho_types, ONLY: qs_rho_get, qs_rho_type
86 : USE qs_scf_types, ONLY: qs_scf_env_type
87 : USE scf_control_types, ONLY: scf_control_type
88 : USE input_section_types, ONLY: section_vals_get_subs_vals, section_vals_type
89 : USE string_utilities, ONLY: integer_to_string
90 : USE tblite_types, ONLY: tblite_type, deallocate_tblite_type, allocate_tblite_type
91 : USE virial_types, ONLY: virial_type
92 : USE xtb_types, ONLY: get_xtb_atom_param, xtb_atom_type
93 : USE xtb_types, ONLY: xtb_atom_type
94 :
95 : #include "./base/base_uses.f90"
96 : IMPLICIT NONE
97 :
98 : PRIVATE
99 :
100 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tblite_interface'
101 :
102 : INTEGER, PARAMETER :: dip_n = 3
103 : INTEGER, PARAMETER :: quad_n = 6
104 : REAL(KIND=dp), PARAMETER :: same_atom = 0.00001_dp
105 : REAL(KIND=dp), PARAMETER :: tblite_scc_pconv = 2.0E-5_dp
106 :
107 : PUBLIC :: tb_set_calculator, tb_init_geometry, tb_init_wf
108 : PUBLIC :: tb_get_basis, build_tblite_matrices
109 : PUBLIC :: tb_get_energy, tb_update_charges, tb_ham_add_coulomb
110 : PUBLIC :: tb_native_scc_mixer_active
111 : PUBLIC :: tb_scf_mixer_error
112 : PUBLIC :: tb_get_multipole
113 : PUBLIC :: tb_derive_dH_diag, tb_derive_dH_off
114 : PUBLIC :: tb_reference_cli_compare
115 :
116 : CONTAINS
117 :
118 : ! **************************************************************************************************
119 : !> \brief intialize geometry objects ...
120 : !> \param qs_env ...
121 : !> \param tb ...
122 : ! **************************************************************************************************
123 82 : SUBROUTINE tb_init_geometry(qs_env, tb)
124 :
125 : TYPE(qs_environment_type), POINTER :: qs_env
126 : TYPE(tblite_type), POINTER :: tb
127 :
128 : #if defined(__TBLITE)
129 :
130 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tblite_init_geometry'
131 :
132 : TYPE(cell_type), POINTER :: cell
133 82 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
134 82 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
135 : INTEGER :: iatom, natom
136 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xyz
137 : INTEGER :: handle, ikind
138 : INTEGER, DIMENSION(3) :: periodic
139 : LOGICAL, DIMENSION(3) :: lperiod
140 :
141 82 : CALL timeset(routineN, handle)
142 :
143 : !get info from environment vaiarable
144 82 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell, qs_kind_set=qs_kind_set)
145 :
146 : !get information about particles
147 82 : natom = SIZE(particle_set)
148 246 : ALLOCATE (xyz(3, natom))
149 82 : CALL allocate_tblite_type(tb)
150 246 : ALLOCATE (tb%el_num(natom))
151 490 : tb%el_num = -9
152 490 : DO iatom = 1, natom
153 1632 : xyz(:, iatom) = particle_set(iatom)%r(:)
154 408 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
155 408 : CALL get_qs_kind(qs_kind_set(ikind), zatom=tb%el_num(iatom))
156 898 : IF (tb%el_num(iatom) < 1 .OR. tb%el_num(iatom) > 85) THEN
157 0 : CPABORT("only elements 1-85 are supported by tblite")
158 : END IF
159 : END DO
160 :
161 : !get information about cell / lattice
162 82 : CALL get_cell(cell=cell, periodic=periodic)
163 82 : lperiod(1) = periodic(1) == 1
164 82 : lperiod(2) = periodic(2) == 1
165 82 : lperiod(3) = periodic(3) == 1
166 :
167 : !prepare for the call to the dispersion function
168 82 : CALL new(tb%mol, tb%el_num, xyz, lattice=cell%hmat, periodic=lperiod)
169 :
170 82 : DEALLOCATE (xyz)
171 :
172 82 : CALL timestop(handle)
173 :
174 : #else
175 : MARK_USED(qs_env)
176 : MARK_USED(tb)
177 : CPABORT("Built without TBLITE")
178 : #endif
179 :
180 82 : END SUBROUTINE tb_init_geometry
181 :
182 : ! **************************************************************************************************
183 : !> \brief updating coordinates...
184 : !> \param qs_env ...
185 : !> \param tb ...
186 : ! **************************************************************************************************
187 1396 : SUBROUTINE tb_update_geometry(qs_env, tb)
188 :
189 : TYPE(qs_environment_type) :: qs_env
190 : TYPE(tblite_type) :: tb
191 :
192 : #if defined(__TBLITE)
193 :
194 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tblite_update_geometry'
195 :
196 : TYPE(cell_type), POINTER :: cell
197 1396 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
198 : INTEGER :: iatom, natom
199 1396 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xyz
200 : INTEGER :: handle
201 :
202 1396 : CALL timeset(routineN, handle)
203 :
204 : !get info from environment vaiarable
205 1396 : NULLIFY (cell)
206 1396 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell)
207 :
208 : !get information about particles
209 1396 : natom = SIZE(particle_set)
210 4188 : ALLOCATE (xyz(3, natom))
211 8270 : DO iatom = 1, natom
212 28892 : xyz(:, iatom) = particle_set(iatom)%r(:)
213 : END DO
214 28892 : tb%mol%xyz(:, :) = xyz
215 18148 : tb%mol%lattice(:, :) = cell%hmat
216 :
217 1396 : DEALLOCATE (xyz)
218 :
219 1396 : CALL timestop(handle)
220 :
221 : #else
222 : MARK_USED(qs_env)
223 : MARK_USED(tb)
224 : CPABORT("Built without TBLITE")
225 : #endif
226 :
227 1396 : END SUBROUTINE tb_update_geometry
228 :
229 : ! **************************************************************************************************
230 : !> \brief initialize wavefunction ...
231 : !> \param tb ...
232 : ! **************************************************************************************************
233 82 : SUBROUTINE tb_init_wf(tb)
234 :
235 : TYPE(tblite_type), POINTER :: tb
236 :
237 : #if defined(__TBLITE)
238 :
239 : INTEGER, PARAMETER :: nSpin = 1 !number of spin channels
240 :
241 : TYPE(scf_info) :: info
242 :
243 82 : info = tb%calc%variable_info()
244 0 : IF (info%charge > shell_resolved) CPABORT("tblite: no support for orbital resolved charge")
245 82 : IF (info%dipole > atom_resolved) CPABORT("tblite: no support for shell resolved dipole moment")
246 82 : IF (info%quadrupole > atom_resolved) &
247 0 : CPABORT("tblite: no support shell resolved quadrupole moment")
248 :
249 82 : CALL new_wavefunction(tb%wfn, tb%mol%nat, tb%calc%bas%nsh, tb%calc%bas%nao, nSpin, 0.0_dp)
250 82 : CALL tb_reset_mixer(tb)
251 :
252 82 : CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
253 :
254 : !allocate quantities later required
255 574 : ALLOCATE (tb%e_hal(tb%mol%nat), tb%e_rep(tb%mol%nat), tb%e_disp(tb%mol%nat))
256 410 : ALLOCATE (tb%e_scd(tb%mol%nat), tb%e_es(tb%mol%nat))
257 246 : ALLOCATE (tb%selfenergy(tb%calc%bas%nsh))
258 246 : IF (ALLOCATED(tb%calc%ncoord)) ALLOCATE (tb%cn(tb%mol%nat))
259 :
260 : #else
261 : MARK_USED(tb)
262 : CPABORT("Built without TBLITE")
263 : #endif
264 :
265 82 : END SUBROUTINE tb_init_wf
266 :
267 : ! **************************************************************************************************
268 : !> \brief Reset tblite's internal SCC mixer for a new CP2K SCF cycle.
269 : !> \param tb ...
270 : ! **************************************************************************************************
271 1382 : SUBROUTINE tb_reset_mixer(tb)
272 :
273 : TYPE(tblite_type), POINTER :: tb
274 :
275 : #if defined(__TBLITE)
276 :
277 : TYPE(scf_info) :: info
278 :
279 1382 : info = tb%calc%variable_info()
280 2682 : IF (ALLOCATED(tb%mixer)) DEALLOCATE (tb%mixer)
281 : CALL new_mixer(tb%mixer, tb%calc%max_iter, &
282 : tb%wfn%nspin*get_mixer_dimension(tb%mol, tb%calc%bas, info), &
283 1382 : tb%calc%mixer_damping)
284 :
285 : #else
286 : MARK_USED(tb)
287 : CPABORT("Built without TBLITE")
288 : #endif
289 :
290 1382 : END SUBROUTINE tb_reset_mixer
291 :
292 : ! **************************************************************************************************
293 : !> \brief Configure tblite's internal SCC mixer from CP2K input.
294 : !> \param tb ...
295 : !> \param max_scf ...
296 : !> \param damping ...
297 : ! **************************************************************************************************
298 1300 : SUBROUTINE tb_configure_mixer(tb, max_scf, damping)
299 :
300 : TYPE(tblite_type), POINTER :: tb
301 : INTEGER, INTENT(IN) :: max_scf
302 : REAL(KIND=dp), INTENT(IN) :: damping
303 :
304 : #if defined(__TBLITE)
305 :
306 1300 : IF (max_scf < 1) CPABORT("tblite SCC mixer MAX_SCF must be positive")
307 1300 : IF (damping <= 0.0_dp) CPABORT("tblite SCC mixer damping must be positive")
308 :
309 1300 : tb%calc%max_iter = max_scf
310 1300 : tb%calc%mixer_damping = damping
311 :
312 : #else
313 : MARK_USED(tb)
314 : MARK_USED(max_scf)
315 : MARK_USED(damping)
316 : CPABORT("Built without TBLITE")
317 : #endif
318 :
319 1300 : END SUBROUTINE tb_configure_mixer
320 :
321 : ! **************************************************************************************************
322 : !> \brief Return whether the tblite native SCC mixer is active for this run.
323 : !> \param dft_control ...
324 : !> \return ...
325 : ! **************************************************************************************************
326 58938 : FUNCTION tb_native_scc_mixer_active(dft_control) RESULT(use_native_mixer)
327 :
328 : TYPE(dft_control_type), POINTER :: dft_control
329 : LOGICAL :: use_native_mixer
330 :
331 58938 : use_native_mixer = .FALSE.
332 58938 : IF (.NOT. ASSOCIATED(dft_control)) RETURN
333 :
334 58938 : SELECT CASE (dft_control%qs_control%xtb_control%tblite_scc_mixer)
335 : CASE (tblite_scc_mixer_auto)
336 : use_native_mixer = .TRUE.
337 : CASE (tblite_scc_mixer_tblite)
338 1414 : use_native_mixer = .TRUE.
339 : CASE (tblite_scc_mixer_cp2k, tblite_scc_mixer_none)
340 1414 : use_native_mixer = .FALSE.
341 : CASE DEFAULT
342 58938 : CPABORT("Unknown tblite SCC mixer")
343 : END SELECT
344 :
345 : END FUNCTION tb_native_scc_mixer_active
346 :
347 : ! **************************************************************************************************
348 : !> \brief Return the native tblite SCC mixer residual for CP2K SCF convergence.
349 : !> \param dft_control ...
350 : !> \param tb ...
351 : !> \param eps_scf ...
352 : !> \return ...
353 : ! **************************************************************************************************
354 14426 : FUNCTION tb_scf_mixer_error(dft_control, tb, eps_scf) RESULT(mixer_error)
355 :
356 : TYPE(dft_control_type), POINTER :: dft_control
357 : TYPE(tblite_type), POINTER :: tb
358 : REAL(KIND=dp), INTENT(IN) :: eps_scf
359 : REAL(KIND=dp) :: mixer_error
360 :
361 : REAL(KIND=dp) :: raw_error
362 :
363 14426 : mixer_error = 0.0_dp
364 :
365 : #if defined(__TBLITE)
366 14426 : IF (.NOT. ASSOCIATED(tb)) RETURN
367 14426 : IF (.NOT. tb_native_scc_mixer_active(dft_control)) RETURN
368 13744 : IF (ALLOCATED(tb%mixer)) THEN
369 13744 : raw_error = REAL(tb%mixer%get_error(), KIND=dp)
370 13744 : IF (eps_scf > 0.0_dp) THEN
371 13744 : mixer_error = eps_scf*raw_error/tblite_scc_pconv
372 : ELSE
373 : mixer_error = raw_error
374 : END IF
375 : END IF
376 : #else
377 : MARK_USED(dft_control)
378 : MARK_USED(tb)
379 : MARK_USED(eps_scf)
380 : #endif
381 :
382 : END FUNCTION tb_scf_mixer_error
383 :
384 : ! **************************************************************************************************
385 : !> \brief ...
386 : !> \param tb ...
387 : !> \param typ ...
388 : ! **************************************************************************************************
389 82 : SUBROUTINE tb_set_calculator(tb, typ)
390 :
391 : TYPE(tblite_type), POINTER :: tb
392 : INTEGER :: typ
393 :
394 : #if defined(__TBLITE)
395 :
396 82 : TYPE(error_type), ALLOCATABLE :: error
397 :
398 82 : SELECT CASE (typ)
399 : CASE default
400 0 : CPABORT("Unknown xtb type")
401 : CASE (gfn1xtb)
402 34 : CALL new_gfn1_calculator(tb%calc, tb%mol, error)
403 : CASE (gfn2xtb)
404 34 : CALL new_gfn2_calculator(tb%calc, tb%mol, error)
405 : CASE (ipea1xtb)
406 82 : CALL new_ipea1_calculator(tb%calc, tb%mol, error)
407 : END SELECT
408 :
409 : #else
410 : MARK_USED(tb)
411 : MARK_USED(typ)
412 : CPABORT("Built without TBLITE")
413 : #endif
414 :
415 82 : END SUBROUTINE tb_set_calculator
416 :
417 : ! **************************************************************************************************
418 : !> \brief ...
419 : !> \param qs_env ...
420 : !> \param tb ...
421 : !> \param para_env ...
422 : ! **************************************************************************************************
423 1396 : SUBROUTINE tb_init_ham(qs_env, tb, para_env)
424 :
425 : TYPE(qs_environment_type) :: qs_env
426 : TYPE(tblite_type) :: tb
427 : TYPE(mp_para_env_type) :: para_env
428 :
429 : #if defined(__TBLITE)
430 :
431 1396 : TYPE(container_cache) :: hcache, rcache
432 :
433 8270 : tb%e_hal = 0.0_dp
434 8270 : tb%e_rep = 0.0_dp
435 8270 : tb%e_disp = 0.0_dp
436 1396 : IF (ALLOCATED(tb%grad)) THEN
437 1038 : tb%grad = 0.0_dp
438 54 : CALL tb_zero_force(qs_env)
439 : END IF
440 18148 : tb%sigma = 0.0_dp
441 :
442 1396 : IF (ALLOCATED(tb%calc%halogen)) THEN
443 720 : CALL tb%calc%halogen%update(tb%mol, hcache)
444 720 : IF (ALLOCATED(tb%grad)) THEN
445 506 : tb%grad = 0.0_dp
446 : CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal, &
447 26 : & tb%grad, tb%sigma)
448 : CALL tb_dump_sigma_component("after_halogen", tb%sigma, para_env)
449 26 : CALL tb_grad2force(qs_env, tb, para_env, 0)
450 : ELSE
451 694 : CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal)
452 : END IF
453 : END IF
454 :
455 1396 : IF (ALLOCATED(tb%calc%repulsion)) THEN
456 1396 : CALL tb%calc%repulsion%update(tb%mol, rcache)
457 1396 : IF (ALLOCATED(tb%grad)) THEN
458 1038 : tb%grad = 0.0_dp
459 : CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep, &
460 54 : & tb%grad, tb%sigma)
461 : CALL tb_dump_sigma_component("after_repulsion", tb%sigma, para_env)
462 54 : CALL tb_grad2force(qs_env, tb, para_env, 1)
463 : ELSE
464 1342 : CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep)
465 : END IF
466 : END IF
467 :
468 1396 : IF (ALLOCATED(tb%calc%dispersion)) THEN
469 1396 : CALL tb%calc%dispersion%update(tb%mol, tb%dcache)
470 1396 : IF (ALLOCATED(tb%grad)) THEN
471 1038 : tb%grad = 0.0_dp
472 : CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp, &
473 54 : & tb%grad, tb%sigma)
474 : CALL tb_dump_sigma_component("after_dispersion_static", tb%sigma, para_env)
475 54 : CALL tb_grad2force(qs_env, tb, para_env, 2)
476 : ELSE
477 1342 : CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp)
478 : END IF
479 : END IF
480 :
481 1396 : CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
482 1396 : IF (ALLOCATED(tb%calc%coulomb)) THEN
483 1396 : CALL tb%calc%coulomb%update(tb%mol, tb%cache)
484 : END IF
485 :
486 1396 : IF (ALLOCATED(tb%grad)) THEN
487 54 : IF (ALLOCATED(tb%calc%ncoord)) THEN
488 54 : CALL tb%calc%ncoord%get_cn(tb%mol, tb%cn, tb%dcndr, tb%dcndL)
489 : END IF
490 : CALL get_selfenergy(tb%calc%h0, tb%mol%id, tb%calc%bas%ish_at, &
491 54 : & tb%calc%bas%nsh_id, cn=tb%cn, selfenergy=tb%selfenergy, dsedcn=tb%dsedcn)
492 : ELSE
493 1342 : IF (ALLOCATED(tb%calc%ncoord)) THEN
494 1342 : CALL tb%calc%ncoord%get_cn(tb%mol, tb%cn)
495 : END IF
496 : CALL get_selfenergy(tb%calc%h0, tb%mol%id, tb%calc%bas%ish_at, &
497 1342 : & tb%calc%bas%nsh_id, cn=tb%cn, selfenergy=tb%selfenergy, dsedcn=tb%dsedcn)
498 : END IF
499 :
500 : #else
501 : MARK_USED(qs_env)
502 : MARK_USED(tb)
503 : MARK_USED(para_env)
504 : CPABORT("Built without TBLITE")
505 : #endif
506 :
507 1396 : END SUBROUTINE tb_init_ham
508 :
509 : ! **************************************************************************************************
510 : !> \brief ...
511 : !> \param qs_env ...
512 : !> \param tb ...
513 : !> \param energy ...
514 : ! **************************************************************************************************
515 14426 : SUBROUTINE tb_get_energy(qs_env, tb, energy)
516 :
517 : TYPE(qs_environment_type), POINTER :: qs_env
518 : TYPE(tblite_type), POINTER :: tb
519 : TYPE(qs_energy_type), POINTER :: energy
520 :
521 : #if defined(__TBLITE)
522 :
523 : INTEGER :: iounit
524 : TYPE(cp_logger_type), POINTER :: logger
525 : TYPE(section_vals_type), POINTER :: scf_section
526 :
527 14426 : NULLIFY (scf_section, logger)
528 :
529 14426 : logger => cp_get_default_logger()
530 14426 : iounit = cp_logger_get_default_io_unit(logger)
531 14426 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
532 :
533 96154 : energy%repulsive = SUM(tb%e_rep)
534 96154 : energy%el_stat = SUM(tb%e_es)
535 96154 : energy%dispersion = SUM(tb%e_disp)
536 96154 : energy%dispersion_sc = SUM(tb%e_scd)
537 96154 : energy%xtb_xb_inter = SUM(tb%e_hal)
538 :
539 : energy%total = energy%core + energy%repulsive + energy%el_stat + energy%dispersion &
540 : + energy%dispersion_sc + energy%xtb_xb_inter + energy%kTS &
541 14426 : + energy%efield + energy%qmmm_el
542 :
543 : iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
544 14426 : extension=".scfLog")
545 14426 : IF (iounit > 0) THEN
546 : WRITE (UNIT=iounit, FMT="(/,(T9,A,T60,F20.10))") &
547 147 : "Repulsive pair potential energy: ", energy%repulsive, &
548 147 : "Zeroth order Hamiltonian energy: ", energy%core, &
549 147 : "Electrostatic energy: ", energy%el_stat, &
550 147 : "Self-consistent dispersion energy: ", energy%dispersion_sc, &
551 294 : "Non-self consistent dispersion energy: ", energy%dispersion
552 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
553 147 : "Correction for halogen bonding: ", energy%xtb_xb_inter
554 147 : IF (ABS(energy%efield) > 1.e-12_dp) THEN
555 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
556 0 : "Electric field interaction energy: ", energy%efield
557 : END IF
558 147 : IF (qs_env%qmmm) THEN
559 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
560 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
561 : END IF
562 : END IF
563 : CALL cp_print_key_finished_output(iounit, logger, scf_section, &
564 14426 : "PRINT%DETAILED_ENERGY")
565 :
566 : #else
567 : MARK_USED(qs_env)
568 : MARK_USED(tb)
569 : MARK_USED(energy)
570 : CPABORT("Built without TBLITE")
571 : #endif
572 :
573 14426 : END SUBROUTINE tb_get_energy
574 :
575 : ! **************************************************************************************************
576 : !> \brief ...
577 : !> \param tb ...
578 : !> \param gto_basis_set ...
579 : !> \param element_symbol ...
580 : !> \param param ...
581 : !> \param occ ...
582 : ! **************************************************************************************************
583 178 : SUBROUTINE tb_get_basis(tb, gto_basis_set, element_symbol, param, occ)
584 :
585 : TYPE(tblite_type), POINTER :: tb
586 : TYPE(gto_basis_set_type), POINTER :: gto_basis_set
587 : CHARACTER(len=2), INTENT(IN) :: element_symbol
588 : TYPE(xtb_atom_type), POINTER :: param
589 : INTEGER, DIMENSION(5), INTENT(out) :: occ
590 :
591 : #if defined(__TBLITE)
592 :
593 : REAL(KIND=dp) :: docc
594 : CHARACTER(LEN=default_string_length) :: sng
595 : INTEGER :: ang, i_type, id_atom, ind_ao, ipgf, iSH, &
596 : ishell, ityp, maxl, mprim, natorb, &
597 : nset, nshell
598 : LOGICAL :: do_ortho
599 :
600 178 : CALL allocate_gto_basis_set(gto_basis_set)
601 :
602 : !identifying element in the bas data
603 178 : CALL symbol_to_number(i_type, element_symbol)
604 352 : DO id_atom = 1, tb%mol%nat
605 352 : IF (i_type == tb%el_num(id_atom)) EXIT
606 : END DO
607 178 : param%z = i_type
608 178 : param%symbol = element_symbol
609 178 : param%defined = .TRUE.
610 178 : ityp = tb%mol%id(id_atom)
611 :
612 : !getting size information
613 178 : nset = tb%calc%bas%nsh_id(ityp)
614 178 : nshell = 1
615 178 : mprim = 0
616 582 : DO ishell = 1, nset
617 582 : mprim = MAX(mprim, tb%calc%bas%cgto(ishell, ityp)%nprim)
618 : END DO
619 178 : param%nshell = nset
620 178 : natorb = 0
621 :
622 : !write basis set information
623 178 : CALL integer_to_string(mprim, sng)
624 178 : gto_basis_set%name = element_symbol//"_STO-"//TRIM(sng)//"G"
625 178 : gto_basis_set%nset = nset
626 178 : CALL reallocate(gto_basis_set%lmax, 1, nset)
627 178 : CALL reallocate(gto_basis_set%lmin, 1, nset)
628 178 : CALL reallocate(gto_basis_set%npgf, 1, nset)
629 178 : CALL reallocate(gto_basis_set%nshell, 1, nset)
630 178 : CALL reallocate(gto_basis_set%n, 1, 1, 1, nset)
631 178 : CALL reallocate(gto_basis_set%l, 1, 1, 1, nset)
632 178 : CALL reallocate(gto_basis_set%zet, 1, mprim, 1, nset)
633 178 : CALL reallocate(gto_basis_set%gcc, 1, mprim, 1, 1, 1, nset)
634 :
635 178 : ind_ao = 0
636 178 : maxl = 0
637 582 : DO ishell = 1, nset
638 404 : ang = tb%calc%bas%cgto(ishell, ityp)%ang
639 404 : natorb = natorb + (2*ang + 1)
640 404 : param%lval(ishell) = ang
641 404 : maxl = MAX(ang, maxl)
642 404 : gto_basis_set%lmax(ishell) = ang
643 404 : gto_basis_set%lmin(ishell) = ang
644 404 : gto_basis_set%npgf(ishell) = tb%calc%bas%cgto(ishell, ityp)%nprim
645 404 : gto_basis_set%nshell(ishell) = nshell
646 404 : gto_basis_set%n(1, ishell) = ang + 1
647 404 : gto_basis_set%l(1, ishell) = ang
648 2502 : DO ipgf = 1, gto_basis_set%npgf(ishell)
649 2098 : gto_basis_set%gcc(ipgf, 1, ishell) = tb%calc%bas%cgto(ishell, ityp)%coeff(ipgf)
650 2502 : gto_basis_set%zet(ipgf, ishell) = tb%calc%bas%cgto(ishell, ityp)%alpha(ipgf)
651 : END DO
652 1458 : DO ipgf = 1, (2*ang + 1)
653 876 : ind_ao = ind_ao + 1
654 876 : param%lao(ind_ao) = ang
655 1280 : param%nao(ind_ao) = ishell
656 : END DO
657 : END DO
658 :
659 178 : do_ortho = .FALSE.
660 178 : CALL process_gto_basis(gto_basis_set, do_ortho, nset, maxl)
661 :
662 : !setting additional values in parameter
663 178 : param%rcut = get_cutoff(tb%calc%bas)
664 178 : param%natorb = natorb
665 178 : param%lmax = maxl !max angular momentum
666 :
667 : !getting occupation
668 178 : occ = 0
669 178 : docc = 0.0_dp
670 178 : IF (tb%calc%bas%nsh_at(id_atom) > 5) CPABORT("too many shells in tblite")
671 582 : DO iSh = 1, tb%calc%bas%nsh_at(id_atom)
672 404 : occ(iSh) = NINT(tb%calc%h0%refocc(iSh, ityp) + docc)
673 404 : docc = docc + tb%calc%h0%refocc(iSh, ityp) - REAL(occ(iSh))
674 582 : param%occupation(iSh) = occ(iSh)
675 : END DO
676 178 : IF (ABS(docc) > 0.1_dp) CPABORT("Getting occupation numbers from tblite fails")
677 1068 : param%zeff = SUM(occ) !effective core charge
678 :
679 : !set normalization process
680 178 : gto_basis_set%norm_type = 3
681 :
682 : #else
683 : occ = 0
684 : MARK_USED(tb)
685 : MARK_USED(gto_basis_set)
686 : MARK_USED(element_symbol)
687 : MARK_USED(param)
688 : CPABORT("Built without TBLITE")
689 : #endif
690 :
691 178 : END SUBROUTINE tb_get_basis
692 :
693 : ! **************************************************************************************************
694 : !> \brief ...
695 : !> \param qs_env ...
696 : !> \param calculate_forces ...
697 : ! **************************************************************************************************
698 1396 : SUBROUTINE build_tblite_matrices(qs_env, calculate_forces)
699 :
700 : TYPE(qs_environment_type), POINTER :: qs_env
701 : LOGICAL, INTENT(IN) :: calculate_forces
702 :
703 : #if defined(__TBLITE)
704 :
705 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_tblite_matrices'
706 :
707 : INTEGER :: handle, maxder, nderivatives, nimg, img, nkind, i, j, &
708 : ic, iw, iatom, jatom, ikind, jkind, iset, jset, n1, n2, icol, &
709 : irow, ia, ib, sgfa, sgfb, atom_a, atom_b, ldsab, nseta, nsetb, &
710 : natorb_a, natorb_b, sgfa0, mp_test_mode, i0, j0
711 : LOGICAL :: found, norml1, norml2, use_arnoldi
712 : REAL(KIND=dp) :: dr, rr
713 : INTEGER, DIMENSION(3) :: cell
714 : REAL(KIND=dp) :: hij, shpoly
715 : REAL(KIND=dp), DIMENSION(2) :: condnum
716 : REAL(KIND=dp), DIMENSION(3) :: rij
717 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
718 : INTEGER :: debug_status
719 : CHARACTER(LEN=32) :: debug_value
720 : #endif
721 1396 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
722 1396 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
723 1396 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork, native_s, native_h, native_lattr
724 1396 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint, hint, native_dip, native_quad
725 1396 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min
726 1396 : INTEGER, DIMENSION(:), POINTER :: npgfa, npgfb, nsgfa, nsgfb
727 1396 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
728 1396 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
729 1396 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, zeta, zetb, scon_a, scon_b
730 1396 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sblock, fblock
731 :
732 1396 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
733 1396 : TYPE(adjacency_list) :: native_list
734 : TYPE(atprop_type), POINTER :: atprop
735 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
736 : TYPE(cp_logger_type), POINTER :: logger
737 1396 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s, matrix_p, matrix_w
738 : TYPE(dft_control_type), POINTER :: dft_control
739 1396 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
740 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
741 1396 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
742 : TYPE(kpoint_type), POINTER :: kpoints
743 1396 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
744 : TYPE(neighbor_list_iterator_p_type), &
745 1396 : DIMENSION(:), POINTER :: nl_iterator
746 : TYPE(mp_para_env_type), POINTER :: para_env
747 : TYPE(qs_energy_type), POINTER :: energy
748 : TYPE(qs_ks_env_type), POINTER :: ks_env
749 1396 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
750 : TYPE(qs_rho_type), POINTER :: rho
751 : TYPE(tblite_type), POINTER :: tb
752 : TYPE(tb_hamiltonian), POINTER :: h0
753 : TYPE(virial_type), POINTER :: virial
754 :
755 1396 : CALL timeset(routineN, handle)
756 :
757 1396 : mp_test_mode = 0
758 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
759 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_MULTIPOLE_TEST", debug_value, STATUS=debug_status)
760 : IF (debug_status == 0) READ (debug_value, *, IOSTAT=debug_status) mp_test_mode
761 : #endif
762 :
763 1396 : NULLIFY (ks_env, energy, atomic_kind_set, qs_kind_set)
764 1396 : NULLIFY (matrix_h, matrix_s, atprop, dft_control)
765 1396 : NULLIFY (sab_orb, rho, tb, kpoints, cell_to_index)
766 :
767 : CALL get_qs_env(qs_env=qs_env, &
768 : ks_env=ks_env, para_env=para_env, &
769 : energy=energy, &
770 : atomic_kind_set=atomic_kind_set, &
771 : qs_kind_set=qs_kind_set, &
772 : matrix_h_kp=matrix_h, &
773 : matrix_s_kp=matrix_s, &
774 : atprop=atprop, &
775 : dft_control=dft_control, &
776 : sab_orb=sab_orb, &
777 1396 : rho=rho, tb_tblite=tb)
778 1396 : h0 => tb%calc%h0
779 :
780 : !update geometry (required for debug / geometry optimization)
781 1396 : CALL tb_update_geometry(qs_env, tb)
782 :
783 1396 : nkind = SIZE(atomic_kind_set)
784 1396 : nderivatives = 0
785 1396 : IF (calculate_forces) THEN
786 54 : nderivatives = 1
787 54 : IF (ALLOCATED(tb%grad)) DEALLOCATE (tb%grad)
788 162 : ALLOCATE (tb%grad(3, tb%mol%nat))
789 54 : IF (ALLOCATED(tb%dsedcn)) DEALLOCATE (tb%dsedcn)
790 162 : ALLOCATE (tb%dsedcn(tb%calc%bas%nsh))
791 54 : IF (ALLOCATED(tb%calc%ncoord)) THEN
792 54 : IF (ALLOCATED(tb%dcndr)) DEALLOCATE (tb%dcndr)
793 216 : ALLOCATE (tb%dcndr(3, tb%mol%nat, tb%mol%nat))
794 54 : IF (ALLOCATED(tb%dcndL)) DEALLOCATE (tb%dcndL)
795 162 : ALLOCATE (tb%dcndL(3, 3, tb%mol%nat))
796 : END IF
797 : ELSE
798 1342 : IF (ALLOCATED(tb%grad)) DEALLOCATE (tb%grad)
799 1342 : IF (ALLOCATED(tb%dcndr)) DEALLOCATE (tb%dcndr)
800 1342 : IF (ALLOCATED(tb%dcndL)) DEALLOCATE (tb%dcndL)
801 : END IF
802 1396 : maxder = ncoset(nderivatives)
803 1396 : nimg = dft_control%nimages
804 1396 : IF (nimg > 1) THEN
805 426 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
806 426 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
807 : END IF
808 :
809 : !intialise hamiltonian
810 1396 : CALL tb_init_ham(qs_env, tb, para_env)
811 :
812 : ! get density matrtix
813 1396 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
814 :
815 : ! set up matrices for force calculations
816 1396 : IF (calculate_forces) THEN
817 54 : NULLIFY (force, matrix_w, virial)
818 : CALL get_qs_env(qs_env=qs_env, &
819 : matrix_w_kp=matrix_w, &
820 54 : virial=virial, force=force)
821 :
822 54 : IF (SIZE(matrix_p, 1) == 2) THEN
823 0 : DO img = 1, nimg
824 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
825 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
826 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
827 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
828 : END DO
829 : END IF
830 82 : tb%use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
831 : END IF
832 :
833 1396 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
834 :
835 : ! set up basis set lists
836 6740 : ALLOCATE (basis_set_list(nkind))
837 1396 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
838 :
839 : ! allocate overlap matrix
840 1396 : CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
841 : CALL create_sab_matrix(ks_env, matrix_s, "OVERLAP MATRIX", basis_set_list, basis_set_list, &
842 1396 : sab_orb, .TRUE.)
843 1396 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
844 :
845 : ! initialize H matrix
846 1396 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
847 24028 : DO img = 1, nimg
848 22632 : ALLOCATE (matrix_h(1, img)%matrix)
849 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, img)%matrix, &
850 22632 : name="HAMILTONIAN MATRIX")
851 24028 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
852 : END DO
853 1396 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
854 :
855 : IF (mp_test_mode == 81 .OR. mp_test_mode == 82 .OR. &
856 : mp_test_mode == 232 .OR. mp_test_mode == 233) THEN
857 : IF (nimg > 1) CPABORT("Native tblite matrix test with k-points not implemented")
858 : ALLOCATE (native_s(tb%calc%bas%nao, tb%calc%bas%nao))
859 : ALLOCATE (native_h(tb%calc%bas%nao, tb%calc%bas%nao))
860 : ALLOCATE (native_dip(dip_n, tb%calc%bas%nao, tb%calc%bas%nao))
861 : ALLOCATE (native_quad(quad_n, tb%calc%bas%nao, tb%calc%bas%nao))
862 : CALL get_lattice_points(tb%mol%periodic, tb%mol%lattice, get_cutoff(tb%calc%bas, 1.0_dp), native_lattr)
863 : CALL new_adjacency_list(native_list, tb%mol, native_lattr, get_cutoff(tb%calc%bas, 1.0_dp))
864 : CALL get_hamiltonian(tb%mol, native_lattr, native_list, tb%calc%bas, tb%calc%h0, tb%selfenergy, &
865 : native_s, native_dip, native_quad, native_h)
866 :
867 : NULLIFY (nl_iterator)
868 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
869 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
870 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom)
871 : icol = MAX(iatom, jatom)
872 : irow = MIN(iatom, jatom)
873 : i0 = tb%calc%bas%iao_sh(tb%calc%bas%ish_at(irow) + 1)
874 : j0 = tb%calc%bas%iao_sh(tb%calc%bas%ish_at(icol) + 1)
875 :
876 : NULLIFY (sblock, fblock)
877 : CALL dbcsr_get_block_p(matrix=matrix_s(1, 1)%matrix, row=irow, col=icol, BLOCK=sblock, found=found)
878 : CPASSERT(found)
879 : CALL dbcsr_get_block_p(matrix=matrix_h(1, 1)%matrix, row=irow, col=icol, BLOCK=fblock, found=found)
880 : CPASSERT(found)
881 :
882 : DO j = 1, SIZE(sblock, 1)
883 : DO i = 1, SIZE(sblock, 2)
884 : sblock(j, i) = native_s(i0 + j, j0 + i)
885 : fblock(j, i) = native_h(i0 + j, j0 + i)
886 : END DO
887 : END DO
888 : END DO
889 : CALL neighbor_list_iterator_release(nl_iterator)
890 :
891 : DO img = 1, nimg
892 : DO i = 1, SIZE(matrix_s, 1)
893 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
894 : END DO
895 : DO i = 1, SIZE(matrix_h, 1)
896 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
897 : END DO
898 : END DO
899 :
900 : IF (dft_control%qs_control%xtb_control%tblite_method == gfn2xtb) &
901 : CALL tb_get_multipole(qs_env, tb)
902 :
903 : DEALLOCATE (native_s, native_h, native_dip, native_quad, native_lattr)
904 : DEALLOCATE (basis_set_list)
905 : CALL timestop(handle)
906 : RETURN
907 : END IF
908 :
909 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
910 1396 : NULLIFY (nl_iterator)
911 1396 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
912 418446 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
913 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
914 417050 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
915 :
916 417050 : icol = MAX(iatom, jatom)
917 417050 : irow = MIN(iatom, jatom)
918 417050 : IF (icol == jatom) THEN
919 952636 : rij = -rij
920 238159 : i = ikind
921 238159 : ikind = jkind
922 238159 : jkind = i
923 : END IF
924 :
925 1668200 : dr = NORM2(rij(:))
926 :
927 417050 : IF (nimg == 1) THEN
928 : ic = 1
929 : ELSE
930 53098 : ic = cell_to_index(cell(1), cell(2), cell(3))
931 53098 : CPASSERT(ic > 0)
932 : END IF
933 417050 : NULLIFY (sblock)
934 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
935 417050 : row=irow, col=icol, BLOCK=sblock, found=found)
936 417050 : CPASSERT(found)
937 417050 : NULLIFY (fblock)
938 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
939 417050 : row=irow, col=icol, BLOCK=fblock, found=found)
940 417050 : CPASSERT(found)
941 :
942 : ! --------- Overlap
943 : !get basis information
944 417050 : basis_set_a => basis_set_list(ikind)%gto_basis_set
945 417050 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
946 417050 : basis_set_b => basis_set_list(jkind)%gto_basis_set
947 417050 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
948 417050 : atom_a = atom_of_kind(icol)
949 417050 : atom_b = atom_of_kind(irow)
950 : ! basis a
951 417050 : first_sgfa => basis_set_a%first_sgf
952 417050 : la_max => basis_set_a%lmax
953 417050 : la_min => basis_set_a%lmin
954 417050 : npgfa => basis_set_a%npgf
955 417050 : nseta = basis_set_a%nset
956 417050 : nsgfa => basis_set_a%nsgf_set
957 417050 : rpgfa => basis_set_a%pgf_radius
958 417050 : set_radius_a => basis_set_a%set_radius
959 417050 : scon_a => basis_set_a%scon
960 417050 : zeta => basis_set_a%zet
961 : ! basis b
962 417050 : first_sgfb => basis_set_b%first_sgf
963 417050 : lb_max => basis_set_b%lmax
964 417050 : lb_min => basis_set_b%lmin
965 417050 : npgfb => basis_set_b%npgf
966 417050 : nsetb = basis_set_b%nset
967 417050 : nsgfb => basis_set_b%nsgf_set
968 417050 : rpgfb => basis_set_b%pgf_radius
969 417050 : set_radius_b => basis_set_b%set_radius
970 417050 : scon_b => basis_set_b%scon
971 417050 : zetb => basis_set_b%zet
972 :
973 417050 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
974 3336400 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
975 417050 : natorb_a = 0
976 1649357 : DO iset = 1, nseta
977 1649357 : natorb_a = natorb_a + (2*basis_set_a%l(1, iset) + 1)
978 : END DO
979 : natorb_b = 0
980 1649892 : DO iset = 1, nsetb
981 1649892 : natorb_b = natorb_b + (2*basis_set_b%l(1, iset) + 1)
982 : END DO
983 2085250 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
984 40187137 : sint = 0.0_dp
985 1668200 : ALLOCATE (hint(natorb_a, natorb_b, maxder))
986 40187137 : hint = 0.0_dp
987 :
988 : !----------------- overlap integrals
989 1649357 : DO iset = 1, nseta
990 1232307 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
991 1232307 : sgfa = first_sgfa(1, iset)
992 5310236 : DO jset = 1, nsetb
993 3660879 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
994 2614166 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
995 2614166 : sgfb = first_sgfb(1, jset)
996 2614166 : IF (calculate_forces) THEN
997 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
998 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
999 76580 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
1000 : ELSE
1001 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1002 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1003 2537586 : rij, sab=oint(:, :, 1))
1004 : END IF
1005 : ! Contraction
1006 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1007 2614166 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1008 2614166 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
1009 3846473 : IF (calculate_forces) THEN
1010 306320 : DO i = 2, 4
1011 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1012 229740 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
1013 306320 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
1014 : END DO
1015 : END IF
1016 : END DO
1017 : END DO
1018 :
1019 : ! update S matrix
1020 417050 : IF (icol <= irow) THEN
1021 9306963 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
1022 : ELSE
1023 27265741 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
1024 : END IF
1025 :
1026 : ! --------- Hamiltonian
1027 417050 : IF (icol == irow .AND. dr < same_atom) THEN
1028 : !get diagonal F matrix from selfenergy
1029 3437 : n1 = tb%calc%bas%ish_at(icol)
1030 12523 : DO iset = 1, nseta
1031 9086 : sgfa = first_sgfa(1, iset)
1032 9086 : hij = tb%selfenergy(n1 + iset)
1033 36841 : DO ia = sgfa, sgfa + nsgfa(iset) - 1
1034 33404 : hint(ia, ia, 1) = hij
1035 : END DO
1036 : END DO
1037 : ELSE
1038 : !get off-diagonal F matrix
1039 413613 : rr = SQRT(dr/(h0%rad(jkind) + h0%rad(ikind)))
1040 413613 : n1 = tb%calc%bas%ish_at(icol)
1041 413613 : sgfa0 = 1
1042 1636834 : DO iset = 1, nseta
1043 1223221 : sgfa = first_sgfa(1, iset)
1044 1223221 : sgfa0 = sgfa0 + tb%calc%bas%cgto(iset, tb%mol%id(icol))%nprim
1045 1223221 : n2 = tb%calc%bas%ish_at(irow)
1046 5272591 : DO jset = 1, nsetb
1047 3635757 : sgfb = first_sgfb(1, jset)
1048 : shpoly = (1.0_dp + h0%shpoly(iset, ikind)*rr) &
1049 3635757 : *(1.0_dp + h0%shpoly(jset, jkind)*rr)
1050 : hij = 0.5_dp*(tb%selfenergy(n1 + iset) + tb%selfenergy(n2 + jset)) &
1051 3635757 : *h0%hscale(iset, jset, ikind, jkind)*shpoly
1052 15670919 : DO ia = sgfa, sgfa + nsgfa(iset) - 1
1053 46757761 : DO ib = sgfb, sgfb + nsgfb(jset) - 1
1054 43122004 : hint(ia, ib, 1) = hij*sint(ia, ib, 1)
1055 : END DO
1056 : END DO
1057 : END DO
1058 : END DO
1059 : END IF
1060 :
1061 : ! update F matrix
1062 417050 : IF (icol <= irow) THEN
1063 9306963 : fblock(:, :) = fblock(:, :) + hint(:, :, 1)
1064 : ELSE
1065 27265741 : fblock(:, :) = fblock(:, :) + TRANSPOSE(hint(:, :, 1))
1066 : END IF
1067 :
1068 1252546 : DEALLOCATE (oint, owork, sint, hint)
1069 :
1070 : END DO
1071 1396 : CALL neighbor_list_iterator_release(nl_iterator)
1072 :
1073 24028 : DO img = 1, nimg
1074 48006 : DO i = 1, SIZE(matrix_s, 1)
1075 48006 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1076 : END DO
1077 46660 : DO i = 1, SIZE(matrix_h, 1)
1078 45264 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1079 : END DO
1080 : END DO
1081 :
1082 : !compute multipole moments for gfn2
1083 1396 : IF (dft_control%qs_control%xtb_control%tblite_method == gfn2xtb) &
1084 676 : CALL tb_get_multipole(qs_env, tb)
1085 :
1086 : ! output overlap information
1087 1396 : NULLIFY (logger)
1088 1396 : logger => cp_get_default_logger()
1089 1396 : IF (.NOT. calculate_forces) THEN
1090 1342 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
1091 : "DFT%PRINT%OVERLAP_CONDITION") /= 0) THEN
1092 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
1093 2 : extension=".Log")
1094 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
1095 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
1096 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
1097 2 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
1098 2 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
1099 : END IF
1100 : END IF
1101 :
1102 1396 : DEALLOCATE (basis_set_list)
1103 :
1104 1396 : CALL timestop(handle)
1105 :
1106 : #else
1107 : MARK_USED(qs_env)
1108 : MARK_USED(calculate_forces)
1109 : CPABORT("Built without TBLITE")
1110 : #endif
1111 :
1112 2792 : END SUBROUTINE build_tblite_matrices
1113 :
1114 : ! **************************************************************************************************
1115 : !> \brief ...
1116 : !> \param qs_env ...
1117 : !> \param dft_control ...
1118 : !> \param tb ...
1119 : !> \param calculate_forces ...
1120 : !> \param use_rho ...
1121 : ! **************************************************************************************************
1122 30158 : SUBROUTINE tb_update_charges(qs_env, dft_control, tb, calculate_forces, use_rho)
1123 :
1124 : TYPE(qs_environment_type), POINTER :: qs_env
1125 : TYPE(dft_control_type), POINTER :: dft_control
1126 : TYPE(tblite_type), POINTER :: tb
1127 : LOGICAL, INTENT(IN) :: calculate_forces
1128 : LOGICAL, INTENT(IN) :: use_rho
1129 :
1130 : #if defined(__TBLITE)
1131 :
1132 : INTEGER :: iatom, ikind, is, ns, atom_a, ii, im, irow, icol, i, j
1133 : INTEGER :: iao, jao, ish, ispin, i0, j0, nspin
1134 : INTEGER :: nimg, nkind, nsgf, natorb, na, n_mix_cols, mix_offset
1135 : INTEGER :: n_atom, max_orb, max_shell
1136 : INTEGER :: mp_test_mode, raw_state_status, raw_state_unit, &
1137 : state_status, state_unit, last_mix_slot, ia_local, ns_mix
1138 : LOGICAL :: discard_mixed_output, do_combined_mixing, do_dipole, do_quadrupole, &
1139 : native_sign_mixing, phased_combined_mixing, skip_charge_mixing, &
1140 : reuse_native_input, skip_scf_dispersion, skip_scf_dispersion_energy, &
1141 : skip_scf_dispersion_gradient, skip_scf_dispersion_potential, &
1142 : use_native_mixer, use_no_mixer
1143 : REAL(KIND=dp) :: mix_alpha, norm, moment_alpha, new_charge, new_moment, pao
1144 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
1145 : INTEGER :: debug_status
1146 : CHARACTER(LEN=32) :: debug_value
1147 : #endif
1148 : CHARACTER(LEN=default_path_length) :: raw_state_file, state_file
1149 : INTEGER, DIMENSION(5) :: occ
1150 : INTEGER, DIMENSION(25) :: lao
1151 : INTEGER, DIMENSION(25) :: nao
1152 30158 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ch_atom, ch_shell, ch_ref, ch_orb, mix_vars, prev_vars
1153 30158 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, ao_dip, ao_quad, native_lattr
1154 30158 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: native_p
1155 :
1156 30158 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1157 30158 : TYPE(adjacency_list) :: native_list
1158 : TYPE(dbcsr_iterator_type) :: iter
1159 30158 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_p
1160 30158 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
1161 : TYPE(dbcsr_type), POINTER :: s_matrix
1162 30158 : TYPE(error_type), ALLOCATABLE :: error
1163 30158 : TYPE(integral_type) :: native_ints
1164 : TYPE(mp_para_env_type), POINTER :: para_env
1165 30158 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1166 30158 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block
1167 30158 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1168 : TYPE(qs_rho_type), POINTER :: rho
1169 : TYPE(qs_scf_env_type), POINTER :: scf_env
1170 : TYPE(scf_control_type), POINTER :: scf_control
1171 : TYPE(xtb_atom_type), POINTER :: xtb_kind
1172 :
1173 : ! also compute multipoles needed by GFN2
1174 30158 : do_dipole = .FALSE.
1175 30158 : do_quadrupole = .FALSE.
1176 30158 : mp_test_mode = 0
1177 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
1178 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_MULTIPOLE_TEST", debug_value, STATUS=debug_status)
1179 : IF (debug_status == 0) READ (debug_value, *, IOSTAT=debug_status) mp_test_mode
1180 : #endif
1181 30158 : skip_scf_dispersion = .FALSE.
1182 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
1183 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DEBUG_SKIP_SCF_DISPERSION", debug_value, STATUS=debug_status)
1184 : IF (debug_status == 0) READ (debug_value, *, IOSTAT=debug_status) skip_scf_dispersion
1185 : #endif
1186 30158 : skip_scf_dispersion_energy = skip_scf_dispersion
1187 30158 : skip_scf_dispersion_gradient = skip_scf_dispersion
1188 30158 : skip_scf_dispersion_potential = skip_scf_dispersion
1189 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
1190 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DEBUG_SKIP_SCF_DISPERSION_ENERGY", debug_value, STATUS=debug_status)
1191 : IF (debug_status == 0) READ (debug_value, *, IOSTAT=debug_status) skip_scf_dispersion_energy
1192 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DEBUG_SKIP_SCF_DISPERSION_GRADIENT", debug_value, STATUS=debug_status)
1193 : IF (debug_status == 0) READ (debug_value, *, IOSTAT=debug_status) skip_scf_dispersion_gradient
1194 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DEBUG_SKIP_SCF_DISPERSION_POTENTIAL", debug_value, STATUS=debug_status)
1195 : IF (debug_status == 0) READ (debug_value, *, IOSTAT=debug_status) skip_scf_dispersion_potential
1196 : #endif
1197 58902 : SELECT CASE (dft_control%qs_control%xtb_control%tblite_scc_mixer)
1198 : CASE (tblite_scc_mixer_auto)
1199 28744 : use_native_mixer = tb_native_scc_mixer_active(dft_control)
1200 28744 : use_no_mixer = .FALSE.
1201 : CASE (tblite_scc_mixer_tblite)
1202 : use_native_mixer = .TRUE.
1203 1414 : use_no_mixer = .FALSE.
1204 : CASE (tblite_scc_mixer_cp2k)
1205 1414 : use_native_mixer = .FALSE.
1206 1414 : use_no_mixer = .FALSE.
1207 : CASE (tblite_scc_mixer_none)
1208 0 : use_native_mixer = .FALSE.
1209 0 : use_no_mixer = .TRUE.
1210 : CASE DEFAULT
1211 0 : CPABORT("Unknown tblite SCC mixer")
1212 : END SELECT
1213 :
1214 : ! compute mulliken charges required for charge update
1215 30158 : NULLIFY (particle_set, qs_kind_set, atomic_kind_set, scf_control)
1216 : CALL get_qs_env(qs_env=qs_env, scf_env=scf_env, particle_set=particle_set, qs_kind_set=qs_kind_set, &
1217 : atomic_kind_set=atomic_kind_set, matrix_s_kp=matrix_s, rho=rho, para_env=para_env, &
1218 30158 : scf_control=scf_control)
1219 30158 : IF (use_native_mixer .AND. use_rho .AND. (.NOT. calculate_forces) .AND. &
1220 : scf_env%iter_count == 1) THEN
1221 : CALL tb_configure_mixer(tb, scf_control%max_scf, &
1222 1300 : dft_control%qs_control%xtb_control%tblite_mixer_damping)
1223 1300 : CALL tb_reset_mixer(tb)
1224 : END IF
1225 30158 : NULLIFY (matrix_p)
1226 30158 : IF (use_rho) THEN
1227 15732 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1228 : ELSE
1229 14426 : matrix_p => scf_env%p_mix_new
1230 : END IF
1231 30158 : IF (ASSOCIATED(tb%dipbra)) do_dipole = .TRUE.
1232 30158 : IF (ASSOCIATED(tb%quadbra)) do_quadrupole = .TRUE.
1233 30158 : reuse_native_input = .FALSE.
1234 30158 : IF (use_native_mixer .AND. scf_env%iter_count == 1) THEN
1235 4772 : reuse_native_input = ANY(ABS(tb%wfn%qsh) > 1.0E-14_dp)
1236 2592 : IF (do_dipole) reuse_native_input = reuse_native_input .OR. &
1237 2500 : ANY(ABS(tb%wfn%dpat) > 1.0E-14_dp)
1238 2592 : IF (do_quadrupole) reuse_native_input = reuse_native_input .OR. &
1239 3424 : ANY(ABS(tb%wfn%qpat) > 1.0E-14_dp)
1240 : END IF
1241 30158 : n_atom = SIZE(particle_set)
1242 30158 : nkind = SIZE(atomic_kind_set)
1243 30158 : nimg = dft_control%nimages
1244 30158 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
1245 120632 : ALLOCATE (aocg(nsgf, n_atom))
1246 1611896 : aocg = 0.0_dp
1247 74194 : IF (do_dipole) ALLOCATE (ao_dip(n_atom, dip_n))
1248 74194 : IF (do_quadrupole) ALLOCATE (ao_quad(n_atom, quad_n))
1249 30158 : max_orb = 0
1250 30158 : max_shell = 0
1251 78080 : DO ikind = 1, nkind
1252 47922 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1253 47922 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb)
1254 78080 : max_orb = MAX(max_orb, natorb)
1255 : END DO
1256 200048 : DO is = 1, n_atom
1257 200048 : max_shell = MAX(max_shell, tb%calc%bas%nsh_at(is))
1258 : END DO
1259 180948 : ALLOCATE (ch_atom(n_atom, 1), ch_shell(n_atom, max_shell))
1260 180948 : ALLOCATE (ch_orb(max_orb, n_atom), ch_ref(max_orb, n_atom))
1261 230206 : ch_atom = 0.0_dp
1262 604734 : ch_shell = 0.0_dp
1263 1611896 : ch_orb = 0.0_dp
1264 1611896 : ch_ref = 0.0_dp
1265 30158 : IF (nimg > 1) THEN
1266 7872 : CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
1267 7872 : IF (do_dipole) THEN
1268 25504 : DO im = 1, dip_n
1269 25504 : CALL contract_dens_kp(matrix_p, tb%dipbra, tb%dipket, im, dip_n, ao_dip(:, im), para_env)
1270 : END DO
1271 : END IF
1272 7872 : IF (do_quadrupole) THEN
1273 44632 : DO im = 1, quad_n
1274 44632 : CALL contract_dens_kp(matrix_p, tb%quadbra, tb%quadket, im, quad_n, ao_quad(:, im), para_env)
1275 : END DO
1276 : END IF
1277 : ELSE
1278 22286 : NULLIFY (p_matrix, s_matrix)
1279 22286 : p_matrix => matrix_p(:, 1)
1280 22286 : s_matrix => matrix_s(1, 1)%matrix
1281 22286 : CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
1282 22286 : IF (do_dipole) THEN
1283 62568 : DO im = 1, dip_n
1284 62568 : CALL contract_dens(p_matrix, tb%dipbra(im)%matrix, tb%dipket(im)%matrix, ao_dip(:, im), para_env)
1285 : END DO
1286 : END IF
1287 22286 : IF (do_quadrupole) THEN
1288 109494 : DO im = 1, quad_n
1289 109494 : CALL contract_dens(p_matrix, tb%quadbra(im)%matrix, tb%quadket(im)%matrix, ao_quad(:, im), para_env)
1290 : END DO
1291 : END IF
1292 : END IF
1293 30158 : NULLIFY (xtb_kind)
1294 78080 : DO ikind = 1, nkind
1295 47922 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
1296 47922 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1297 47922 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, nao=nao, occupation=occ)
1298 295892 : DO iatom = 1, na
1299 169890 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1300 1490088 : DO is = 1, natorb
1301 1320198 : ns = lao(is) + 1
1302 1320198 : norm = 2*lao(is) + 1
1303 1320198 : ch_ref(is, atom_a) = tb%calc%h0%refocc(nao(is), ikind)/norm
1304 1320198 : ch_orb(is, atom_a) = aocg(is, atom_a) - ch_ref(is, atom_a)
1305 1490088 : ch_shell(atom_a, ns) = ch_orb(is, atom_a) + ch_shell(atom_a, ns)
1306 : END DO
1307 1629660 : ch_atom(atom_a, 1) = SUM(ch_orb(:, atom_a))
1308 : END DO
1309 : END DO
1310 30158 : DEALLOCATE (aocg)
1311 :
1312 30158 : raw_state_status = 1
1313 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
1314 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_RAW_STATE_DUMP", raw_state_file, STATUS=raw_state_status)
1315 : #endif
1316 : IF (raw_state_status == 0) THEN
1317 : OPEN (NEWUNIT=raw_state_unit, FILE=TRIM(raw_state_file), STATUS="REPLACE", ACTION="WRITE")
1318 : WRITE (raw_state_unit, *) "qat"
1319 : DO iatom = 1, n_atom
1320 : WRITE (raw_state_unit, "(I0,1X,ES24.16)") iatom, -ch_atom(iatom, 1)
1321 : END DO
1322 : WRITE (raw_state_unit, *) "qsh"
1323 : DO iatom = 1, n_atom
1324 : DO is = 1, tb%calc%bas%nsh_at(iatom)
1325 : WRITE (raw_state_unit, "(I0,1X,ES24.16)") tb%calc%bas%ish_at(iatom) + is, -ch_shell(iatom, is)
1326 : END DO
1327 : END DO
1328 : IF (do_dipole) THEN
1329 : WRITE (raw_state_unit, *) "dpat"
1330 : DO iatom = 1, n_atom
1331 : WRITE (raw_state_unit, "(I0,3(1X,ES24.16))") iatom, -ao_dip(iatom, :)
1332 : END DO
1333 : END IF
1334 : IF (do_quadrupole) THEN
1335 : WRITE (raw_state_unit, *) "qpat"
1336 : DO iatom = 1, n_atom
1337 : WRITE (raw_state_unit, "(I0,6(1X,ES24.16))") iatom, -ao_quad(iatom, :)
1338 : END DO
1339 : END IF
1340 : CLOSE (raw_state_unit)
1341 : END IF
1342 :
1343 : IF (mp_test_mode == 150) THEN
1344 : n_mix_cols = max_shell
1345 : IF (do_dipole) n_mix_cols = n_mix_cols + dip_n
1346 : IF (do_quadrupole) n_mix_cols = n_mix_cols + quad_n
1347 : ALLOCATE (mix_vars(n_atom, n_mix_cols))
1348 : mix_vars = 0.0_dp
1349 :
1350 : mix_vars(:, 1:max_shell) = -ch_shell(:, 1:max_shell)
1351 : mix_offset = max_shell
1352 : IF (do_dipole) THEN
1353 : mix_vars(:, mix_offset + 1:mix_offset + dip_n) = -ao_dip(:, 1:dip_n)
1354 : mix_offset = mix_offset + dip_n
1355 : END IF
1356 : IF (do_quadrupole) THEN
1357 : mix_vars(:, mix_offset + 1:mix_offset + quad_n) = -ao_quad(:, 1:quad_n)
1358 : END IF
1359 :
1360 : IF (use_rho .AND. (.NOT. calculate_forces) .AND. scf_env%mixing_store%ncall > 0) THEN
1361 : mix_vars = 0.0_dp
1362 : ns_mix = MIN(n_mix_cols, scf_env%mixing_store%max_shell)
1363 : last_mix_slot = MOD(scf_env%mixing_store%ncall - 1, scf_env%mixing_store%nbuffer) + 1
1364 : DO ia_local = 1, scf_env%mixing_store%nat_local
1365 : iatom = scf_env%mixing_store%atlist(ia_local)
1366 : mix_vars(iatom, 1:ns_mix) = scf_env%mixing_store%acharge(ia_local, 1:ns_mix, last_mix_slot)
1367 : END DO
1368 : CALL para_env%sum(mix_vars)
1369 : ELSEIF ((.NOT. use_rho) .AND. (.NOT. calculate_forces)) THEN
1370 : ALLOCATE (prev_vars(n_atom, n_mix_cols))
1371 : prev_vars(:, :) = mix_vars
1372 : IF (scf_env%mixing_store%ncall > 0) THEN
1373 : prev_vars = 0.0_dp
1374 : ns_mix = MIN(n_mix_cols, scf_env%mixing_store%max_shell)
1375 : last_mix_slot = MOD(scf_env%mixing_store%ncall - 1, scf_env%mixing_store%nbuffer) + 1
1376 : DO ia_local = 1, scf_env%mixing_store%nat_local
1377 : iatom = scf_env%mixing_store%atlist(ia_local)
1378 : prev_vars(iatom, 1:ns_mix) = scf_env%mixing_store%acharge(ia_local, 1:ns_mix, last_mix_slot)
1379 : END DO
1380 : CALL para_env%sum(prev_vars)
1381 : END IF
1382 : mix_alpha = scf_env%mixing_store%alpha
1383 : prev_vars(:, :) = (1.0_dp - mix_alpha)*prev_vars + mix_alpha*mix_vars
1384 : scf_env%mixing_store%ncall = scf_env%mixing_store%ncall + 1
1385 : ns_mix = MIN(n_mix_cols, scf_env%mixing_store%max_shell)
1386 : last_mix_slot = MOD(scf_env%mixing_store%ncall - 1, scf_env%mixing_store%nbuffer) + 1
1387 : DO ia_local = 1, scf_env%mixing_store%nat_local
1388 : iatom = scf_env%mixing_store%atlist(ia_local)
1389 : scf_env%mixing_store%acharge(ia_local, 1:ns_mix, last_mix_slot) = prev_vars(iatom, 1:ns_mix)
1390 : END DO
1391 : DEALLOCATE (prev_vars)
1392 : END IF
1393 :
1394 : DO iatom = 1, n_atom
1395 : ii = tb%calc%bas%ish_at(iatom)
1396 : DO is = 1, tb%calc%bas%nsh_at(iatom)
1397 : tb%wfn%qsh(ii + is, 1) = mix_vars(iatom, is)
1398 : END DO
1399 : tb%wfn%qat(iatom, 1) = SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), 1))
1400 : END DO
1401 : mix_offset = max_shell
1402 : IF (do_dipole) THEN
1403 : DO iatom = 1, n_atom
1404 : tb%wfn%dpat(:, iatom, 1) = mix_vars(iatom, mix_offset + 1:mix_offset + dip_n)
1405 : END DO
1406 : mix_offset = mix_offset + dip_n
1407 : DEALLOCATE (ao_dip)
1408 : END IF
1409 : IF (do_quadrupole) THEN
1410 : DO iatom = 1, n_atom
1411 : tb%wfn%qpat(:, iatom, 1) = mix_vars(iatom, mix_offset + 1:mix_offset + quad_n)
1412 : END DO
1413 : DEALLOCATE (ao_quad)
1414 : END IF
1415 : DEALLOCATE (mix_vars)
1416 : ELSEIF (use_native_mixer .OR. mp_test_mode == 120 .OR. mp_test_mode == 220 .OR. &
1417 : mp_test_mode == 230 .OR. mp_test_mode == 231 .OR. &
1418 30158 : mp_test_mode == 232 .OR. mp_test_mode == 233) THEN
1419 28744 : IF (.NOT. ALLOCATED(tb%mixer)) CPABORT("tblite mixer not initialized")
1420 28744 : IF (use_rho .AND. (.NOT. calculate_forces) .AND. scf_env%iter_count > 1) THEN
1421 13648 : CALL tb%mixer%next(error)
1422 13648 : IF (ALLOCATED(error)) CPABORT("tblite native mixer failed")
1423 13648 : CALL tb%mixer%get(tb%wfn%qsh)
1424 92184 : tb%wfn%qat(:, 1) = 0.0_dp
1425 92184 : DO iatom = 1, n_atom
1426 78536 : ii = tb%calc%bas%ish_at(iatom)
1427 313364 : tb%wfn%qat(iatom, 1) = SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), 1))
1428 : END DO
1429 13648 : IF (do_dipole) THEN
1430 9964 : CALL tb%mixer%get(tb%wfn%dpat)
1431 9964 : DEALLOCATE (ao_dip)
1432 : END IF
1433 13648 : IF (do_quadrupole) THEN
1434 9964 : CALL tb%mixer%get(tb%wfn%qpat)
1435 9964 : DEALLOCATE (ao_quad)
1436 : END IF
1437 : ELSE
1438 : IF ((use_native_mixer .OR. mp_test_mode == 230 .OR. mp_test_mode == 231 .OR. &
1439 : mp_test_mode == 232 .OR. mp_test_mode == 233) .AND. &
1440 15096 : use_rho .AND. (.NOT. calculate_forces)) THEN
1441 1300 : IF (.NOT. reuse_native_input) THEN
1442 1170 : tb%wfn%qsh(:, :) = 0.0_dp
1443 560 : tb%wfn%qat(:, :) = 0.0_dp
1444 728 : IF (do_dipole) tb%wfn%dpat(:, :, :) = 0.0_dp
1445 1190 : IF (do_quadrupole) tb%wfn%qpat(:, :, :) = 0.0_dp
1446 : END IF
1447 1300 : IF (do_dipole) DEALLOCATE (ao_dip)
1448 1300 : IF (do_quadrupole) DEALLOCATE (ao_quad)
1449 : ELSE
1450 13796 : IF ((.NOT. use_rho) .AND. (.NOT. calculate_forces)) THEN
1451 13744 : CALL tb%mixer%set(tb%wfn%qsh)
1452 13744 : IF (do_dipole) CALL tb%mixer%set(tb%wfn%dpat)
1453 13744 : IF (do_quadrupole) CALL tb%mixer%set(tb%wfn%qpat)
1454 : END IF
1455 93034 : DO iatom = 1, n_atom
1456 79238 : ii = tb%calc%bas%ish_at(iatom)
1457 302198 : DO is = 1, tb%calc%bas%nsh_at(iatom)
1458 302198 : tb%wfn%qsh(ii + is, 1) = -ch_shell(iatom, is)
1459 : END DO
1460 315994 : tb%wfn%qat(iatom, 1) = SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), 1))
1461 : END DO
1462 13796 : IF (do_dipole) THEN
1463 69064 : DO iatom = 1, n_atom
1464 246154 : tb%wfn%dpat(:, iatom, 1) = -ao_dip(iatom, :)
1465 : END DO
1466 10034 : DEALLOCATE (ao_dip)
1467 : END IF
1468 13796 : IF (do_quadrupole) THEN
1469 69064 : DO iatom = 1, n_atom
1470 423244 : tb%wfn%qpat(:, iatom, 1) = -ao_quad(iatom, :)
1471 : END DO
1472 10034 : DEALLOCATE (ao_quad)
1473 : END IF
1474 13796 : IF ((.NOT. use_rho) .AND. (.NOT. calculate_forces)) THEN
1475 13744 : CALL tb%mixer%diff(tb%wfn%qsh)
1476 13744 : IF (do_dipole) CALL tb%mixer%diff(tb%wfn%dpat)
1477 13744 : IF (do_quadrupole) CALL tb%mixer%diff(tb%wfn%qpat)
1478 : END IF
1479 : END IF
1480 : END IF
1481 : ELSE
1482 : ! charge mixing
1483 1414 : IF (dft_control%qs_control%do_ls_scf) THEN
1484 : !
1485 : ELSE
1486 : ! CP2K SCC mixing has to see all tblite SCC variables, not only shell charges.
1487 : do_combined_mixing = do_dipole .OR. do_quadrupole .OR. &
1488 : (mp_test_mode == 110 .OR. mp_test_mode == 111 .OR. &
1489 : mp_test_mode == 112 .OR. mp_test_mode == 114 .OR. &
1490 : mp_test_mode == 115 .OR. &
1491 : mp_test_mode == 116 .OR. mp_test_mode == 117 .OR. &
1492 : mp_test_mode == 118 .OR. mp_test_mode == 122 .OR. &
1493 : mp_test_mode == 130 .OR. mp_test_mode == 131 .OR. &
1494 : mp_test_mode == 140 .OR. mp_test_mode == 141 .OR. &
1495 1414 : mp_test_mode == 142 .OR. mp_test_mode == 183)
1496 : native_sign_mixing = ((do_dipole .OR. do_quadrupole) .AND. mp_test_mode == 0) .OR. &
1497 : (mp_test_mode == 116 .OR. mp_test_mode == 117 .OR. &
1498 1414 : mp_test_mode == 118)
1499 1414 : phased_combined_mixing = (mp_test_mode == 130 .OR. mp_test_mode == 131)
1500 : discard_mixed_output = (mp_test_mode == 114 .AND. .NOT. use_rho) .OR. &
1501 : (mp_test_mode == 115 .AND. (.NOT. use_rho .OR. calculate_forces)) .OR. &
1502 : (mp_test_mode == 117 .AND. .NOT. use_rho) .OR. &
1503 : (mp_test_mode == 118 .AND. (.NOT. use_rho .OR. calculate_forces)) .OR. &
1504 1414 : (mp_test_mode == 131 .AND. (.NOT. use_rho .OR. calculate_forces))
1505 : skip_charge_mixing = use_no_mixer .OR. &
1506 : (mp_test_mode == 111 .AND. use_rho) .OR. &
1507 : (mp_test_mode == 112 .AND. .NOT. use_rho) .OR. &
1508 1414 : (phased_combined_mixing .AND. use_rho)
1509 : IF (phased_combined_mixing .AND. use_rho .AND. scf_env%mixing_store%ncall > 0 .AND. &
1510 : (mp_test_mode /= 131 .OR. .NOT. calculate_forces)) THEN
1511 : n_mix_cols = max_shell
1512 : IF (do_dipole) n_mix_cols = n_mix_cols + dip_n
1513 : IF (do_quadrupole) n_mix_cols = n_mix_cols + quad_n
1514 : ALLOCATE (mix_vars(n_atom, n_mix_cols))
1515 : mix_vars = 0.0_dp
1516 : ns_mix = MIN(n_mix_cols, scf_env%mixing_store%max_shell)
1517 : last_mix_slot = MOD(scf_env%mixing_store%ncall - 1, scf_env%mixing_store%nbuffer) + 1
1518 : DO ia_local = 1, scf_env%mixing_store%nat_local
1519 : iatom = scf_env%mixing_store%atlist(ia_local)
1520 : mix_vars(iatom, 1:ns_mix) = scf_env%mixing_store%acharge(ia_local, 1:ns_mix, last_mix_slot)
1521 : END DO
1522 : CALL para_env%sum(mix_vars)
1523 : ch_shell(:, 1:max_shell) = mix_vars(:, 1:max_shell)
1524 : mix_offset = max_shell
1525 : IF (do_dipole) THEN
1526 : ao_dip(:, 1:dip_n) = mix_vars(:, mix_offset + 1:mix_offset + dip_n)
1527 : mix_offset = mix_offset + dip_n
1528 : END IF
1529 : IF (do_quadrupole) THEN
1530 : ao_quad(:, 1:quad_n) = mix_vars(:, mix_offset + 1:mix_offset + quad_n)
1531 : END IF
1532 : DEALLOCATE (mix_vars)
1533 : END IF
1534 : IF (skip_charge_mixing) THEN
1535 : !
1536 1414 : ELSEIF (do_combined_mixing) THEN
1537 1414 : n_mix_cols = max_shell
1538 1414 : IF (do_dipole) n_mix_cols = n_mix_cols + dip_n
1539 1414 : IF (do_quadrupole) n_mix_cols = n_mix_cols + quad_n
1540 5656 : ALLOCATE (mix_vars(n_atom, n_mix_cols))
1541 : IF (native_sign_mixing) THEN
1542 15554 : mix_vars(:, 1:max_shell) = -ch_shell(:, 1:max_shell)
1543 : ELSE
1544 : mix_vars(:, 1:max_shell) = ch_shell(:, 1:max_shell)
1545 : END IF
1546 1414 : mix_offset = max_shell
1547 1414 : IF (do_dipole) THEN
1548 : IF (native_sign_mixing) THEN
1549 22624 : mix_vars(:, mix_offset + 1:mix_offset + dip_n) = -ao_dip(:, 1:dip_n)
1550 : ELSE
1551 : mix_vars(:, mix_offset + 1:mix_offset + dip_n) = ao_dip(:, 1:dip_n)
1552 : END IF
1553 : mix_offset = mix_offset + dip_n
1554 : END IF
1555 1414 : IF (do_quadrupole) THEN
1556 : IF (native_sign_mixing) THEN
1557 43834 : mix_vars(:, mix_offset + 1:mix_offset + quad_n) = -ao_quad(:, 1:quad_n)
1558 : ELSE
1559 : mix_vars(:, mix_offset + 1:mix_offset + quad_n) = ao_quad(:, 1:quad_n)
1560 : END IF
1561 : END IF
1562 : CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
1563 1414 : mix_vars, para_env, scf_env%iter_count)
1564 : IF (.NOT. discard_mixed_output) THEN
1565 : IF (native_sign_mixing) THEN
1566 15554 : ch_shell(:, 1:max_shell) = -mix_vars(:, 1:max_shell)
1567 : ELSE
1568 : ch_shell(:, 1:max_shell) = mix_vars(:, 1:max_shell)
1569 : END IF
1570 1414 : mix_offset = max_shell
1571 1414 : IF (do_dipole) THEN
1572 : IF (native_sign_mixing) THEN
1573 22624 : ao_dip(:, 1:dip_n) = -mix_vars(:, mix_offset + 1:mix_offset + dip_n)
1574 : ELSE
1575 : ao_dip(:, 1:dip_n) = mix_vars(:, mix_offset + 1:mix_offset + dip_n)
1576 : END IF
1577 : mix_offset = mix_offset + dip_n
1578 : END IF
1579 1414 : IF (do_quadrupole) THEN
1580 : IF (native_sign_mixing) THEN
1581 43834 : ao_quad(:, 1:quad_n) = -mix_vars(:, mix_offset + 1:mix_offset + quad_n)
1582 : ELSE
1583 : ao_quad(:, 1:quad_n) = mix_vars(:, mix_offset + 1:mix_offset + quad_n)
1584 : END IF
1585 : END IF
1586 : END IF
1587 1414 : DEALLOCATE (mix_vars)
1588 : ELSE
1589 : CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
1590 0 : ch_shell, para_env, scf_env%iter_count)
1591 : END IF
1592 : END IF
1593 :
1594 : !setting new wave function
1595 1414 : CALL tb%pot%reset
1596 7070 : tb%e_es = 0.0_dp
1597 7070 : tb%e_scd = 0.0_dp
1598 7070 : moment_alpha = scf_env%p_mix_alpha
1599 : IF (mp_test_mode == 61 .OR. mp_test_mode == 91) moment_alpha = 0.05_dp
1600 : IF (mp_test_mode == 62 .OR. mp_test_mode == 92) moment_alpha = 0.5_dp
1601 7070 : DO iatom = 1, n_atom
1602 5656 : ii = tb%calc%bas%ish_at(iatom)
1603 14140 : DO is = 1, tb%calc%bas%nsh_at(iatom)
1604 8484 : new_charge = -ch_shell(iatom, is)
1605 : IF ((mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. &
1606 : scf_env%iter_count > 1) THEN
1607 : new_charge = moment_alpha*new_charge + (1.0_dp - moment_alpha)*tb%wfn%qsh(ii + is, 1)
1608 : END IF
1609 14140 : tb%wfn%qsh(ii + is, 1) = new_charge
1610 : END DO
1611 7070 : IF (native_sign_mixing .OR. mp_test_mode == 122) THEN
1612 14140 : tb%wfn%qat(iatom, 1) = SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), 1))
1613 : ELSE
1614 0 : new_charge = -ch_atom(iatom, 1)
1615 : IF ((mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. &
1616 : scf_env%iter_count > 1) THEN
1617 : new_charge = moment_alpha*new_charge + (1.0_dp - moment_alpha)*tb%wfn%qat(iatom, 1)
1618 : END IF
1619 0 : tb%wfn%qat(iatom, 1) = new_charge
1620 : END IF
1621 : END DO
1622 :
1623 1414 : IF (do_dipole) THEN
1624 7070 : DO iatom = 1, n_atom
1625 24038 : DO im = 1, dip_n
1626 16968 : new_moment = -ao_dip(iatom, im)
1627 : IF ((mp_test_mode == 60 .OR. mp_test_mode == 61 .OR. mp_test_mode == 62 .OR. &
1628 : mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. &
1629 : scf_env%iter_count > 1) THEN
1630 : new_moment = moment_alpha*new_moment + (1.0_dp - moment_alpha)*tb%wfn%dpat(im, iatom, 1)
1631 : END IF
1632 22624 : tb%wfn%dpat(im, iatom, 1) = new_moment
1633 : END DO
1634 : END DO
1635 1414 : DEALLOCATE (ao_dip)
1636 : END IF
1637 1414 : IF (do_quadrupole) THEN
1638 7070 : DO iatom = 1, n_atom
1639 41006 : DO im = 1, quad_n
1640 33936 : new_moment = -ao_quad(iatom, im)
1641 : IF ((mp_test_mode == 60 .OR. mp_test_mode == 61 .OR. mp_test_mode == 62 .OR. &
1642 : mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. &
1643 : scf_env%iter_count > 1) THEN
1644 : new_moment = moment_alpha*new_moment + (1.0_dp - moment_alpha)*tb%wfn%qpat(im, iatom, 1)
1645 : END IF
1646 39592 : tb%wfn%qpat(im, iatom, 1) = new_moment
1647 : END DO
1648 : END DO
1649 1414 : DEALLOCATE (ao_quad)
1650 : END IF
1651 : END IF
1652 : IF (mp_test_mode == 2 .OR. mp_test_mode == 7 .OR. mp_test_mode == 8 .OR. &
1653 : mp_test_mode == 13 .OR. mp_test_mode == 14 .OR. mp_test_mode == 15 .OR. &
1654 : mp_test_mode == 16 .OR. mp_test_mode == 17 .OR. mp_test_mode == 18 .OR. &
1655 : mp_test_mode == 20 .OR. mp_test_mode == 21 .OR. &
1656 : mp_test_mode == 32 .OR. mp_test_mode == 33 .OR. &
1657 : mp_test_mode == 34 .OR. mp_test_mode == 35 .OR. &
1658 : mp_test_mode == 36 .OR. mp_test_mode == 37 .OR. &
1659 : mp_test_mode == 10 .OR. mp_test_mode == 11 .OR. mp_test_mode == 12) THEN
1660 : tb%wfn%dpat(:, :, :) = -tb%wfn%dpat(:, :, :)
1661 : tb%wfn%qpat(:, :, :) = -tb%wfn%qpat(:, :, :)
1662 : ELSEIF (mp_test_mode == 4) THEN
1663 : tb%wfn%dpat(:, :, :) = 0.0_dp
1664 : tb%wfn%qpat(:, :, :) = 0.0_dp
1665 : END IF
1666 :
1667 : IF (mp_test_mode == 100 .OR. mp_test_mode == 101 .OR. mp_test_mode == 102) THEN
1668 : nspin = SIZE(p_matrix)
1669 : CALL new_integral(native_ints, tb%calc%bas%nao)
1670 : CALL get_lattice_points(tb%mol%periodic, tb%mol%lattice, get_cutoff(tb%calc%bas, 1.0_dp), native_lattr)
1671 : CALL new_adjacency_list(native_list, tb%mol, native_lattr, get_cutoff(tb%calc%bas, 1.0_dp))
1672 : CALL get_hamiltonian(tb%mol, native_lattr, native_list, tb%calc%bas, tb%calc%h0, tb%selfenergy, &
1673 : native_ints%overlap, native_ints%dipole, native_ints%quadrupole, &
1674 : native_ints%hamiltonian)
1675 : ALLOCATE (native_p(tb%calc%bas%nao, tb%calc%bas%nao, nspin))
1676 : native_p = 0.0_dp
1677 : DO ispin = 1, nspin
1678 : CALL dbcsr_iterator_start(iter, p_matrix(ispin)%matrix)
1679 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1680 : NULLIFY (p_block)
1681 : CALL dbcsr_iterator_next_block(iter, irow, icol, p_block)
1682 : i0 = tb%calc%bas%iao_sh(tb%calc%bas%ish_at(irow) + 1)
1683 : j0 = tb%calc%bas%iao_sh(tb%calc%bas%ish_at(icol) + 1)
1684 : DO j = 1, SIZE(p_block, 1)
1685 : DO i = 1, SIZE(p_block, 2)
1686 : IF (mp_test_mode == 102 .AND. irow /= icol) THEN
1687 : native_p(i0 + j, j0 + i, ispin) = 0.5_dp*p_block(j, i)
1688 : native_p(j0 + i, i0 + j, ispin) = 0.5_dp*p_block(j, i)
1689 : ELSE
1690 : native_p(i0 + j, j0 + i, ispin) = p_block(j, i)
1691 : IF (mp_test_mode == 100 .AND. irow /= icol) &
1692 : native_p(j0 + i, i0 + j, ispin) = p_block(j, i)
1693 : END IF
1694 : END DO
1695 : END DO
1696 : END DO
1697 : CALL dbcsr_iterator_stop(iter)
1698 : END DO
1699 :
1700 : tb%wfn%qsh(:, :) = 0.0_dp
1701 : DO ispin = 1, nspin
1702 : DO iao = 1, tb%calc%bas%nao
1703 : pao = 0.0_dp
1704 : DO jao = 1, tb%calc%bas%nao
1705 : pao = pao + native_p(jao, iao, ispin)*native_ints%overlap(jao, iao)
1706 : END DO
1707 : ish = tb%calc%bas%ao2sh(iao)
1708 : tb%wfn%qsh(ish, ispin) = tb%wfn%qsh(ish, ispin) - pao
1709 : END DO
1710 : END DO
1711 : tb%wfn%qsh(:, 1) = tb%wfn%qsh(:, 1) + tb%wfn%n0sh(:)
1712 :
1713 : tb%wfn%qat(:, :) = 0.0_dp
1714 : DO ispin = 1, nspin
1715 : DO ish = 1, tb%calc%bas%nsh
1716 : tb%wfn%qat(tb%calc%bas%sh2at(ish), ispin) = &
1717 : tb%wfn%qat(tb%calc%bas%sh2at(ish), ispin) + tb%wfn%qsh(ish, ispin)
1718 : END DO
1719 : END DO
1720 :
1721 : IF (do_dipole) THEN
1722 : tb%wfn%dpat(:, :, :) = 0.0_dp
1723 : DO ispin = 1, nspin
1724 : DO iao = 1, tb%calc%bas%nao
1725 : atom_a = tb%calc%bas%ao2at(iao)
1726 : DO jao = 1, tb%calc%bas%nao
1727 : DO im = 1, dip_n
1728 : tb%wfn%dpat(im, atom_a, ispin) = tb%wfn%dpat(im, atom_a, ispin) &
1729 : - native_p(jao, iao, ispin)*native_ints%dipole(im, jao, iao)
1730 : END DO
1731 : END DO
1732 : END DO
1733 : END DO
1734 : END IF
1735 : IF (do_quadrupole) THEN
1736 : tb%wfn%qpat(:, :, :) = 0.0_dp
1737 : DO ispin = 1, nspin
1738 : DO iao = 1, tb%calc%bas%nao
1739 : atom_a = tb%calc%bas%ao2at(iao)
1740 : DO jao = 1, tb%calc%bas%nao
1741 : DO im = 1, quad_n
1742 : tb%wfn%qpat(im, atom_a, ispin) = tb%wfn%qpat(im, atom_a, ispin) &
1743 : - native_p(jao, iao, ispin)*native_ints%quadrupole(im, jao, iao)
1744 : END DO
1745 : END DO
1746 : END DO
1747 : END DO
1748 : END IF
1749 : DEALLOCATE (native_p, native_lattr)
1750 : END IF
1751 :
1752 30158 : CALL tb%pot%reset
1753 200048 : tb%e_es = 0.0_dp
1754 200048 : tb%e_scd = 0.0_dp
1755 30158 : IF (ALLOCATED(tb%calc%coulomb)) THEN
1756 30158 : CALL tb%calc%coulomb%get_potential(tb%mol, tb%cache, tb%wfn, tb%pot)
1757 30158 : CALL tb%calc%coulomb%get_energy(tb%mol, tb%cache, tb%wfn, tb%e_es)
1758 : END IF
1759 30158 : IF (ALLOCATED(tb%calc%dispersion)) THEN
1760 : IF (.NOT. skip_scf_dispersion_potential) &
1761 30158 : CALL tb%calc%dispersion%get_potential(tb%mol, tb%dcache, tb%wfn, tb%pot)
1762 : IF (.NOT. skip_scf_dispersion_energy) &
1763 30158 : CALL tb%calc%dispersion%get_energy(tb%mol, tb%dcache, tb%wfn, tb%e_scd)
1764 : END IF
1765 :
1766 30158 : state_status = 1
1767 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
1768 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_STATE_DUMP", state_file, STATUS=state_status)
1769 : #endif
1770 : IF (state_status == 0) THEN
1771 : OPEN (NEWUNIT=state_unit, FILE=TRIM(state_file), STATUS="REPLACE", ACTION="WRITE")
1772 : WRITE (state_unit, *) "qat"
1773 : DO iatom = 1, n_atom
1774 : WRITE (state_unit, "(I0,1X,ES24.16)") iatom, tb%wfn%qat(iatom, 1)
1775 : END DO
1776 : WRITE (state_unit, *) "qsh"
1777 : DO iatom = 1, SIZE(tb%wfn%qsh, 1)
1778 : WRITE (state_unit, "(I0,1X,ES24.16)") iatom, tb%wfn%qsh(iatom, 1)
1779 : END DO
1780 : IF (do_dipole) THEN
1781 : WRITE (state_unit, *) "dpat"
1782 : DO iatom = 1, n_atom
1783 : WRITE (state_unit, "(I0,3(1X,ES24.16))") iatom, tb%wfn%dpat(:, iatom, 1)
1784 : END DO
1785 : END IF
1786 : IF (do_quadrupole) THEN
1787 : WRITE (state_unit, *) "qpat"
1788 : DO iatom = 1, n_atom
1789 : WRITE (state_unit, "(I0,6(1X,ES24.16))") iatom, tb%wfn%qpat(:, iatom, 1)
1790 : END DO
1791 : END IF
1792 : CLOSE (state_unit)
1793 : END IF
1794 :
1795 30158 : IF (calculate_forces) THEN
1796 54 : IF (ALLOCATED(tb%calc%coulomb)) THEN
1797 1038 : tb%grad = 0.0_dp
1798 54 : CALL tb%calc%coulomb%get_gradient(tb%mol, tb%cache, tb%wfn, tb%grad, tb%sigma)
1799 54 : CALL tb_dump_sigma_component("after_coulomb", tb%sigma, para_env)
1800 54 : CALL tb_grad2force(qs_env, tb, para_env, 3)
1801 : END IF
1802 :
1803 54 : IF (ALLOCATED(tb%calc%dispersion) .AND. .NOT. skip_scf_dispersion_gradient) THEN
1804 1038 : tb%grad = 0.0_dp
1805 54 : CALL tb%calc%dispersion%get_gradient(tb%mol, tb%dcache, tb%wfn, tb%grad, tb%sigma)
1806 54 : CALL tb_dump_sigma_component("after_dispersion_scf", tb%sigma, para_env)
1807 54 : CALL tb_grad2force(qs_env, tb, para_env, 2)
1808 : END IF
1809 : END IF
1810 :
1811 30158 : DEALLOCATE (ch_atom, ch_shell, ch_orb, ch_ref)
1812 :
1813 : #else
1814 : MARK_USED(qs_env)
1815 : MARK_USED(tb)
1816 : MARK_USED(dft_control)
1817 : MARK_USED(calculate_forces)
1818 : MARK_USED(use_rho)
1819 : CPABORT("Built without TBLITE")
1820 : #endif
1821 :
1822 90474 : END SUBROUTINE tb_update_charges
1823 :
1824 : ! **************************************************************************************************
1825 : !> \brief ...
1826 : !> \param qs_env ...
1827 : !> \param tb ...
1828 : !> \param dft_control ...
1829 : ! **************************************************************************************************
1830 15732 : SUBROUTINE tb_ham_add_coulomb(qs_env, tb, dft_control)
1831 :
1832 : TYPE(qs_environment_type), POINTER :: qs_env
1833 : TYPE(tblite_type), POINTER :: tb
1834 : TYPE(dft_control_type), POINTER :: dft_control
1835 :
1836 : #if defined(__TBLITE)
1837 :
1838 : INTEGER :: ikind, jkind, iatom, jatom, icol, irow
1839 : INTEGER :: ic, id1, id2, id3, iq1, iq2, iq3, iq4, iq5, iq6, &
1840 : is, nimg, ni, nj, i, j, i0, j0
1841 : INTEGER :: la, lb, za, zb
1842 : INTEGER :: mp_test_mode
1843 : LOGICAL :: found
1844 : INTEGER, DIMENSION(3) :: cellind
1845 : INTEGER, DIMENSION(25) :: naoa, naob
1846 : REAL(KIND=dp), DIMENSION(3) :: rij
1847 : REAL(KIND=dp) :: mpfac
1848 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
1849 : INTEGER :: debug_status
1850 : CHARACTER(LEN=32) :: debug_value
1851 : #endif
1852 15732 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, sum_shell
1853 15732 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ashift, bshift, native_lattr
1854 15732 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: native_h1
1855 15732 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksblock, sblock
1856 15732 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1857 15732 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dip_ket1, dip_ket2, dip_ket3, &
1858 15732 : dip_bra1, dip_bra2, dip_bra3
1859 15732 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: quad_ket1, quad_ket2, quad_ket3, &
1860 15732 : quad_ket4, quad_ket5, quad_ket6, &
1861 15732 : quad_bra1, quad_bra2, quad_bra3, &
1862 15732 : quad_bra4, quad_bra5, quad_bra6
1863 :
1864 15732 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1865 15732 : TYPE(adjacency_list) :: native_list
1866 : TYPE(dbcsr_iterator_type) :: iter
1867 15732 : TYPE(integral_type) :: native_ints
1868 15732 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1869 15732 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
1870 : TYPE(kpoint_type), POINTER :: kpoints
1871 : TYPE(neighbor_list_iterator_p_type), &
1872 15732 : DIMENSION(:), POINTER :: nl_iterator
1873 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1874 15732 : POINTER :: n_list
1875 15732 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1876 15732 : TYPE(potential_type) :: native_pot
1877 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
1878 :
1879 15732 : nimg = dft_control%nimages
1880 15732 : mp_test_mode = 0
1881 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
1882 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_MULTIPOLE_TEST", debug_value, STATUS=debug_status)
1883 : IF (debug_status == 0) READ (debug_value, *, IOSTAT=debug_status) mp_test_mode
1884 : #endif
1885 : mpfac = -0.5_dp
1886 : IF (mp_test_mode == 3 .OR. mp_test_mode == 10 .OR. mp_test_mode == 11 .OR. &
1887 : mp_test_mode == 12) mpfac = 0.5_dp
1888 : IF (mp_test_mode == 5) mpfac = 0.0_dp
1889 :
1890 15732 : NULLIFY (matrix_s, ks_matrix, n_list, qs_kind_set)
1891 15732 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list, matrix_s_kp=matrix_s, matrix_ks_kp=ks_matrix, qs_kind_set=qs_kind_set)
1892 :
1893 : IF (mp_test_mode == 83 .OR. mp_test_mode == 183 .OR. &
1894 : mp_test_mode == 220 .OR. mp_test_mode == 231 .OR. &
1895 : mp_test_mode == 233) THEN
1896 : IF (nimg > 1) CPABORT("Native tblite Hamiltonian test with k-points not implemented")
1897 : CALL new_integral(native_ints, tb%calc%bas%nao)
1898 : CALL get_lattice_points(tb%mol%periodic, tb%mol%lattice, get_cutoff(tb%calc%bas, 1.0_dp), native_lattr)
1899 : CALL new_adjacency_list(native_list, tb%mol, native_lattr, get_cutoff(tb%calc%bas, 1.0_dp))
1900 : CALL get_hamiltonian(tb%mol, native_lattr, native_list, tb%calc%bas, tb%calc%h0, tb%selfenergy, &
1901 : native_ints%overlap, native_ints%dipole, native_ints%quadrupole, &
1902 : native_ints%hamiltonian)
1903 : ALLOCATE (native_h1(tb%calc%bas%nao, tb%calc%bas%nao, 1))
1904 : native_pot = tb%pot
1905 : CALL add_pot_to_h1(tb%calc%bas, native_ints, native_pot, native_h1)
1906 :
1907 : CALL dbcsr_iterator_start(iter, ks_matrix(1, 1)%matrix)
1908 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1909 : CALL dbcsr_iterator_next_block(iter, irow, icol, ksblock)
1910 : i0 = tb%calc%bas%iao_sh(tb%calc%bas%ish_at(irow) + 1)
1911 : j0 = tb%calc%bas%iao_sh(tb%calc%bas%ish_at(icol) + 1)
1912 : DO j = 1, SIZE(ksblock, 1)
1913 : DO i = 1, SIZE(ksblock, 2)
1914 : ksblock(j, i) = native_h1(i0 + j, j0 + i, 1)
1915 : END DO
1916 : END DO
1917 : END DO
1918 : CALL dbcsr_iterator_stop(iter)
1919 : DEALLOCATE (native_h1, native_lattr)
1920 : RETURN
1921 : END IF
1922 :
1923 : !creating sum of shell lists
1924 47196 : ALLOCATE (sum_shell(tb%mol%nat))
1925 103894 : i = 0
1926 103894 : DO j = 1, tb%mol%nat
1927 88162 : sum_shell(j) = i
1928 103894 : i = i + tb%calc%bas%nsh_at(j)
1929 : END DO
1930 :
1931 15732 : IF (nimg == 1) THEN
1932 : ! no k-points; all matrices have been transformed to periodic bsf
1933 11596 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1934 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1935 11596 : kind_of=kind_of)
1936 11596 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
1937 151674 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1938 140078 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1939 :
1940 140078 : ikind = kind_of(irow)
1941 140078 : jkind = kind_of(icol)
1942 :
1943 : ! atomic parameters
1944 140078 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
1945 140078 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
1946 140078 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa)
1947 140078 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob)
1948 :
1949 140078 : ni = SIZE(sblock, 1)
1950 560312 : ALLOCATE (ashift(ni, ni))
1951 10850692 : ashift = 0.0_dp
1952 1259542 : DO i = 1, ni
1953 1119464 : la = naoa(i) + sum_shell(irow)
1954 1259542 : ashift(i, i) = tb%pot%vsh(la, 1)
1955 : END DO
1956 :
1957 140078 : nj = SIZE(sblock, 2)
1958 560312 : ALLOCATE (bshift(nj, nj))
1959 10628294 : bshift = 0.0_dp
1960 1232515 : DO j = 1, nj
1961 1092437 : lb = naob(j) + sum_shell(icol)
1962 1232515 : bshift(j, j) = tb%pot%vsh(lb, 1)
1963 : END DO
1964 :
1965 280156 : DO is = 1, SIZE(ks_matrix, 1)
1966 140078 : NULLIFY (ksblock)
1967 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
1968 140078 : row=irow, col=icol, block=ksblock, found=found)
1969 140078 : CPASSERT(found)
1970 560312 : ksblock = ksblock - 0.5_dp*(MATMUL(ashift, sblock) &
1971 716968082 : + MATMUL(sblock, bshift))
1972 : ksblock = ksblock - 0.5_dp*(tb%pot%vat(irow, 1) &
1973 21562920 : + tb%pot%vat(icol, 1))*sblock
1974 : END DO
1975 431830 : DEALLOCATE (ashift, bshift)
1976 : END DO
1977 11596 : CALL dbcsr_iterator_stop(iter)
1978 :
1979 11596 : IF (ASSOCIATED(tb%dipbra)) THEN
1980 8030 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
1981 111536 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1982 103506 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1983 :
1984 103506 : NULLIFY (dip_bra1, dip_bra2, dip_bra3)
1985 : CALL dbcsr_get_block_p(matrix=tb%dipbra(1)%matrix, &
1986 103506 : row=irow, col=icol, BLOCK=dip_bra1, found=found)
1987 103506 : CPASSERT(found)
1988 : CALL dbcsr_get_block_p(matrix=tb%dipbra(2)%matrix, &
1989 103506 : row=irow, col=icol, BLOCK=dip_bra2, found=found)
1990 103506 : CPASSERT(found)
1991 : CALL dbcsr_get_block_p(matrix=tb%dipbra(3)%matrix, &
1992 103506 : row=irow, col=icol, BLOCK=dip_bra3, found=found)
1993 103506 : CPASSERT(found)
1994 103506 : NULLIFY (dip_ket1, dip_ket2, dip_ket3)
1995 : CALL dbcsr_get_block_p(matrix=tb%dipket(1)%matrix, &
1996 103506 : row=irow, col=icol, BLOCK=dip_ket1, found=found)
1997 103506 : CPASSERT(found)
1998 : CALL dbcsr_get_block_p(matrix=tb%dipket(2)%matrix, &
1999 103506 : row=irow, col=icol, BLOCK=dip_ket2, found=found)
2000 103506 : CPASSERT(found)
2001 : CALL dbcsr_get_block_p(matrix=tb%dipket(3)%matrix, &
2002 103506 : row=irow, col=icol, BLOCK=dip_ket3, found=found)
2003 103506 : CPASSERT(found)
2004 :
2005 215042 : DO is = 1, SIZE(ks_matrix, 1)
2006 103506 : NULLIFY (ksblock)
2007 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
2008 103506 : row=irow, col=icol, block=ksblock, found=found)
2009 103506 : CPASSERT(found)
2010 : ksblock = ksblock + mpfac*(dip_ket1*tb%pot%vdp(1, irow, 1) &
2011 : + dip_ket2*tb%pot%vdp(2, irow, 1) &
2012 : + dip_ket3*tb%pot%vdp(3, irow, 1) &
2013 : + dip_bra1*tb%pot%vdp(1, icol, 1) &
2014 : + dip_bra2*tb%pot%vdp(2, icol, 1) &
2015 17357732 : + dip_bra3*tb%pot%vdp(3, icol, 1))
2016 : END DO
2017 : END DO
2018 8030 : CALL dbcsr_iterator_stop(iter)
2019 : END IF
2020 :
2021 11596 : IF (ASSOCIATED(tb%quadbra)) THEN
2022 8030 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
2023 111536 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2024 103506 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
2025 :
2026 103506 : NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
2027 : CALL dbcsr_get_block_p(matrix=tb%quadbra(1)%matrix, &
2028 103506 : row=irow, col=icol, BLOCK=quad_bra1, found=found)
2029 103506 : CPASSERT(found)
2030 : CALL dbcsr_get_block_p(matrix=tb%quadbra(2)%matrix, &
2031 103506 : row=irow, col=icol, BLOCK=quad_bra2, found=found)
2032 103506 : CPASSERT(found)
2033 : CALL dbcsr_get_block_p(matrix=tb%quadbra(3)%matrix, &
2034 103506 : row=irow, col=icol, BLOCK=quad_bra3, found=found)
2035 103506 : CPASSERT(found)
2036 : CALL dbcsr_get_block_p(matrix=tb%quadbra(4)%matrix, &
2037 103506 : row=irow, col=icol, BLOCK=quad_bra4, found=found)
2038 103506 : CPASSERT(found)
2039 : CALL dbcsr_get_block_p(matrix=tb%quadbra(5)%matrix, &
2040 103506 : row=irow, col=icol, BLOCK=quad_bra5, found=found)
2041 103506 : CPASSERT(found)
2042 : CALL dbcsr_get_block_p(matrix=tb%quadbra(6)%matrix, &
2043 103506 : row=irow, col=icol, BLOCK=quad_bra6, found=found)
2044 103506 : CPASSERT(found)
2045 :
2046 103506 : NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
2047 : CALL dbcsr_get_block_p(matrix=tb%quadket(1)%matrix, &
2048 103506 : row=irow, col=icol, BLOCK=quad_ket1, found=found)
2049 103506 : CPASSERT(found)
2050 : CALL dbcsr_get_block_p(matrix=tb%quadket(2)%matrix, &
2051 103506 : row=irow, col=icol, BLOCK=quad_ket2, found=found)
2052 103506 : CPASSERT(found)
2053 : CALL dbcsr_get_block_p(matrix=tb%quadket(3)%matrix, &
2054 103506 : row=irow, col=icol, BLOCK=quad_ket3, found=found)
2055 103506 : CPASSERT(found)
2056 : CALL dbcsr_get_block_p(matrix=tb%quadket(4)%matrix, &
2057 103506 : row=irow, col=icol, BLOCK=quad_ket4, found=found)
2058 103506 : CPASSERT(found)
2059 : CALL dbcsr_get_block_p(matrix=tb%quadket(5)%matrix, &
2060 103506 : row=irow, col=icol, BLOCK=quad_ket5, found=found)
2061 103506 : CPASSERT(found)
2062 : CALL dbcsr_get_block_p(matrix=tb%quadket(6)%matrix, &
2063 103506 : row=irow, col=icol, BLOCK=quad_ket6, found=found)
2064 103506 : CPASSERT(found)
2065 :
2066 215042 : DO is = 1, SIZE(ks_matrix, 1)
2067 103506 : NULLIFY (ksblock)
2068 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
2069 103506 : row=irow, col=icol, block=ksblock, found=found)
2070 103506 : CPASSERT(found)
2071 :
2072 : ksblock = ksblock + mpfac*(quad_ket1*tb%pot%vqp(1, irow, 1) &
2073 : + quad_ket2*tb%pot%vqp(2, irow, 1) &
2074 : + quad_ket3*tb%pot%vqp(3, irow, 1) &
2075 : + quad_ket4*tb%pot%vqp(4, irow, 1) &
2076 : + quad_ket5*tb%pot%vqp(5, irow, 1) &
2077 : + quad_ket6*tb%pot%vqp(6, irow, 1) &
2078 : + quad_bra1*tb%pot%vqp(1, icol, 1) &
2079 : + quad_bra2*tb%pot%vqp(2, icol, 1) &
2080 : + quad_bra3*tb%pot%vqp(3, icol, 1) &
2081 : + quad_bra4*tb%pot%vqp(4, icol, 1) &
2082 : + quad_bra5*tb%pot%vqp(5, icol, 1) &
2083 17357732 : + quad_bra6*tb%pot%vqp(6, icol, 1))
2084 : END DO
2085 : END DO
2086 8030 : CALL dbcsr_iterator_stop(iter)
2087 : END IF
2088 :
2089 : ELSE
2090 4136 : NULLIFY (kpoints)
2091 4136 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
2092 4136 : NULLIFY (cell_to_index)
2093 4136 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
2094 :
2095 4136 : NULLIFY (nl_iterator)
2096 4136 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
2097 414954 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2098 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2099 410818 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
2100 :
2101 410818 : icol = MAX(iatom, jatom)
2102 410818 : irow = MIN(iatom, jatom)
2103 :
2104 410818 : IF (icol == jatom) THEN
2105 982660 : rij = -rij
2106 245665 : i = ikind
2107 245665 : ikind = jkind
2108 245665 : jkind = i
2109 : END IF
2110 :
2111 410818 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2112 410818 : CPASSERT(ic > 0)
2113 :
2114 410818 : NULLIFY (sblock)
2115 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
2116 410818 : row=irow, col=icol, block=sblock, found=found)
2117 410818 : CPASSERT(found)
2118 :
2119 : ! atomic parameters
2120 410818 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
2121 410818 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
2122 410818 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa)
2123 410818 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob)
2124 :
2125 410818 : ni = SIZE(sblock, 1)
2126 1643272 : ALLOCATE (ashift(ni, ni))
2127 37384438 : ashift = 0.0_dp
2128 4108180 : DO i = 1, ni
2129 3697362 : la = naoa(i) + sum_shell(irow)
2130 4108180 : ashift(i, i) = tb%pot%vsh(la, 1)
2131 : END DO
2132 :
2133 410818 : nj = SIZE(sblock, 2)
2134 1643272 : ALLOCATE (bshift(nj, nj))
2135 37384438 : bshift = 0.0_dp
2136 4108180 : DO j = 1, nj
2137 3697362 : lb = naob(j) + sum_shell(icol)
2138 4108180 : bshift(j, j) = tb%pot%vsh(lb, 1)
2139 : END DO
2140 :
2141 821636 : DO is = 1, SIZE(ks_matrix, 1)
2142 410818 : NULLIFY (ksblock)
2143 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
2144 410818 : row=irow, col=icol, block=ksblock, found=found)
2145 410818 : CPASSERT(found)
2146 1643272 : ksblock = ksblock - 0.5_dp*(MATMUL(ashift, sblock) &
2147 2583634402 : + MATMUL(sblock, bshift))
2148 : ksblock = ksblock - 0.5_dp*(tb%pot%vat(irow, 1) &
2149 75590512 : + tb%pot%vat(icol, 1))*sblock
2150 : END DO
2151 1236590 : DEALLOCATE (ashift, bshift)
2152 : END DO
2153 4136 : CALL neighbor_list_iterator_release(nl_iterator)
2154 :
2155 4136 : IF (ASSOCIATED(tb%dipbra)) THEN
2156 3298 : NULLIFY (nl_iterator)
2157 3298 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
2158 231392 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2159 : CALL get_iterator_info(nl_iterator, &
2160 228094 : iatom=iatom, jatom=jatom, cell=cellind)
2161 228094 : icol = MAX(iatom, jatom)
2162 228094 : irow = MIN(iatom, jatom)
2163 228094 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2164 228094 : CPASSERT(ic > 0)
2165 228094 : id1 = 1 + dip_n*(ic - 1)
2166 228094 : id2 = 2 + dip_n*(ic - 1)
2167 228094 : id3 = 3 + dip_n*(ic - 1)
2168 :
2169 228094 : NULLIFY (dip_bra1, dip_bra2, dip_bra3)
2170 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id1)%matrix, &
2171 228094 : row=irow, col=icol, BLOCK=dip_bra1, found=found)
2172 228094 : CPASSERT(found)
2173 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id2)%matrix, &
2174 228094 : row=irow, col=icol, BLOCK=dip_bra2, found=found)
2175 228094 : CPASSERT(found)
2176 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id3)%matrix, &
2177 228094 : row=irow, col=icol, BLOCK=dip_bra3, found=found)
2178 228094 : CPASSERT(found)
2179 228094 : NULLIFY (dip_ket1, dip_ket2, dip_ket3)
2180 : CALL dbcsr_get_block_p(matrix=tb%dipket(id1)%matrix, &
2181 228094 : row=irow, col=icol, BLOCK=dip_ket1, found=found)
2182 228094 : CPASSERT(found)
2183 : CALL dbcsr_get_block_p(matrix=tb%dipket(id2)%matrix, &
2184 228094 : row=irow, col=icol, BLOCK=dip_ket2, found=found)
2185 228094 : CPASSERT(found)
2186 : CALL dbcsr_get_block_p(matrix=tb%dipket(id3)%matrix, &
2187 228094 : row=irow, col=icol, BLOCK=dip_ket3, found=found)
2188 228094 : CPASSERT(found)
2189 :
2190 459486 : DO is = 1, SIZE(ks_matrix, 1)
2191 228094 : NULLIFY (ksblock)
2192 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
2193 228094 : row=irow, col=icol, block=ksblock, found=found)
2194 228094 : CPASSERT(found)
2195 : ksblock = ksblock + mpfac*(dip_ket1*tb%pot%vdp(1, irow, 1) &
2196 : + dip_ket2*tb%pot%vdp(2, irow, 1) &
2197 : + dip_ket3*tb%pot%vdp(3, irow, 1) &
2198 : + dip_bra1*tb%pot%vdp(1, icol, 1) &
2199 : + dip_bra2*tb%pot%vdp(2, icol, 1) &
2200 41969296 : + dip_bra3*tb%pot%vdp(3, icol, 1))
2201 : END DO
2202 : END DO
2203 3298 : CALL neighbor_list_iterator_release(nl_iterator)
2204 : END IF
2205 :
2206 4136 : IF (ASSOCIATED(tb%quadbra)) THEN
2207 3298 : NULLIFY (nl_iterator)
2208 3298 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
2209 231392 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2210 : CALL get_iterator_info(nl_iterator, &
2211 228094 : iatom=iatom, jatom=jatom, cell=cellind)
2212 228094 : icol = MAX(iatom, jatom)
2213 228094 : irow = MIN(iatom, jatom)
2214 228094 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2215 228094 : CPASSERT(ic > 0)
2216 228094 : iq1 = 1 + quad_n*(ic - 1)
2217 228094 : iq2 = 2 + quad_n*(ic - 1)
2218 228094 : iq3 = 3 + quad_n*(ic - 1)
2219 228094 : iq4 = 4 + quad_n*(ic - 1)
2220 228094 : iq5 = 5 + quad_n*(ic - 1)
2221 228094 : iq6 = 6 + quad_n*(ic - 1)
2222 :
2223 228094 : NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
2224 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq1)%matrix, &
2225 228094 : row=irow, col=icol, BLOCK=quad_bra1, found=found)
2226 228094 : CPASSERT(found)
2227 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq2)%matrix, &
2228 228094 : row=irow, col=icol, BLOCK=quad_bra2, found=found)
2229 228094 : CPASSERT(found)
2230 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq3)%matrix, &
2231 228094 : row=irow, col=icol, BLOCK=quad_bra3, found=found)
2232 228094 : CPASSERT(found)
2233 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq4)%matrix, &
2234 228094 : row=irow, col=icol, BLOCK=quad_bra4, found=found)
2235 228094 : CPASSERT(found)
2236 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq5)%matrix, &
2237 228094 : row=irow, col=icol, BLOCK=quad_bra5, found=found)
2238 228094 : CPASSERT(found)
2239 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq6)%matrix, &
2240 228094 : row=irow, col=icol, BLOCK=quad_bra6, found=found)
2241 228094 : CPASSERT(found)
2242 :
2243 228094 : NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
2244 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq1)%matrix, &
2245 228094 : row=irow, col=icol, BLOCK=quad_ket1, found=found)
2246 228094 : CPASSERT(found)
2247 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq2)%matrix, &
2248 228094 : row=irow, col=icol, BLOCK=quad_ket2, found=found)
2249 228094 : CPASSERT(found)
2250 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq3)%matrix, &
2251 228094 : row=irow, col=icol, BLOCK=quad_ket3, found=found)
2252 228094 : CPASSERT(found)
2253 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq4)%matrix, &
2254 228094 : row=irow, col=icol, BLOCK=quad_ket4, found=found)
2255 228094 : CPASSERT(found)
2256 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq5)%matrix, &
2257 228094 : row=irow, col=icol, BLOCK=quad_ket5, found=found)
2258 228094 : CPASSERT(found)
2259 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq6)%matrix, &
2260 228094 : row=irow, col=icol, BLOCK=quad_ket6, found=found)
2261 228094 : CPASSERT(found)
2262 :
2263 459486 : DO is = 1, SIZE(ks_matrix, 1)
2264 228094 : NULLIFY (ksblock)
2265 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
2266 228094 : row=irow, col=icol, block=ksblock, found=found)
2267 228094 : CPASSERT(found)
2268 :
2269 : ksblock = ksblock + mpfac*(quad_ket1*tb%pot%vqp(1, irow, 1) &
2270 : + quad_ket2*tb%pot%vqp(2, irow, 1) &
2271 : + quad_ket3*tb%pot%vqp(3, irow, 1) &
2272 : + quad_ket4*tb%pot%vqp(4, irow, 1) &
2273 : + quad_ket5*tb%pot%vqp(5, irow, 1) &
2274 : + quad_ket6*tb%pot%vqp(6, irow, 1) &
2275 : + quad_bra1*tb%pot%vqp(1, icol, 1) &
2276 : + quad_bra2*tb%pot%vqp(2, icol, 1) &
2277 : + quad_bra3*tb%pot%vqp(3, icol, 1) &
2278 : + quad_bra4*tb%pot%vqp(4, icol, 1) &
2279 : + quad_bra5*tb%pot%vqp(5, icol, 1) &
2280 41969296 : + quad_bra6*tb%pot%vqp(6, icol, 1))
2281 : END DO
2282 : END DO
2283 3298 : CALL neighbor_list_iterator_release(nl_iterator)
2284 : END IF
2285 :
2286 : END IF
2287 :
2288 : #else
2289 : MARK_USED(qs_env)
2290 : MARK_USED(tb)
2291 : MARK_USED(dft_control)
2292 : CPABORT("Built without TBLITE")
2293 : #endif
2294 :
2295 31464 : END SUBROUTINE tb_ham_add_coulomb
2296 :
2297 : ! **************************************************************************************************
2298 : !> \brief ...
2299 : !> \param qs_env ...
2300 : !> \param tb ...
2301 : ! **************************************************************************************************
2302 676 : SUBROUTINE tb_get_multipole(qs_env, tb)
2303 :
2304 : TYPE(qs_environment_type), POINTER :: qs_env
2305 : TYPE(tblite_type), POINTER :: tb
2306 :
2307 : #if defined(__TBLITE)
2308 :
2309 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tb_get_multipole'
2310 :
2311 : INTEGER :: ikind, jkind, iatom, jatom, icol, irow, iset, jset, ityp, jtyp
2312 : INTEGER :: ic, idx, id1, id2, id3, img, iq1, iq2, iq3, iq4, iq5, iq6
2313 : INTEGER :: nkind, natom, handle, nimg, i, j, inda, indb
2314 : INTEGER :: mp_test_mode
2315 : INTEGER :: atom_a, atom_b, nseta, nsetb, ia, ib, ij, i0, j0
2316 : LOGICAL :: found
2317 : REAL(KIND=dp) :: r2
2318 : INTEGER, DIMENSION(3) :: cell
2319 676 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2320 : REAL(KIND=dp), DIMENSION(3) :: rij, mp_int_vec, mp_shift_vec
2321 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
2322 : INTEGER :: debug_status
2323 : CHARACTER(LEN=32) :: debug_value
2324 : #endif
2325 676 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max
2326 676 : INTEGER, DIMENSION(:), POINTER :: nsgfa, nsgfb
2327 676 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
2328 676 : INTEGER, ALLOCATABLE :: atom_of_kind(:)
2329 676 : REAL(KIND=dp), ALLOCATABLE :: stmp(:), native_s(:, :), native_h(:, :), native_lattr(:, :)
2330 676 : REAL(KIND=dp), ALLOCATABLE :: dtmp(:, :), qtmp(:, :), dtmpj(:, :), qtmpj(:, :)
2331 676 : REAL(KIND=dp), ALLOCATABLE :: native_dip(:, :, :), native_quad(:, :, :)
2332 676 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dip_ket1, dip_ket2, dip_ket3, &
2333 676 : dip_bra1, dip_bra2, dip_bra3
2334 676 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: quad_ket1, quad_ket2, quad_ket3, &
2335 676 : quad_ket4, quad_ket5, quad_ket6, &
2336 676 : quad_bra1, quad_bra2, quad_bra3, &
2337 676 : quad_bra4, quad_bra5, quad_bra6
2338 :
2339 676 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2340 676 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
2341 : TYPE(dft_control_type), POINTER :: dft_control
2342 : TYPE(dbcsr_iterator_type) :: iter
2343 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
2344 676 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
2345 676 : TYPE(adjacency_list) :: native_list
2346 : TYPE(kpoint_type), POINTER :: kpoints
2347 676 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
2348 : TYPE(neighbor_list_iterator_p_type), &
2349 676 : DIMENSION(:), POINTER :: nl_iterator
2350 676 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2351 676 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2352 :
2353 676 : CALL timeset(routineN, handle)
2354 :
2355 : !get info from environment vaiarable
2356 676 : NULLIFY (atomic_kind_set, qs_kind_set, sab_orb, particle_set)
2357 676 : NULLIFY (dft_control, matrix_s, kpoints, cell_to_index)
2358 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
2359 : qs_kind_set=qs_kind_set, &
2360 : sab_orb=sab_orb, &
2361 : particle_set=particle_set, &
2362 : dft_control=dft_control, &
2363 : kpoints=kpoints, &
2364 676 : matrix_s_kp=matrix_s)
2365 676 : natom = SIZE(particle_set)
2366 676 : nkind = SIZE(atomic_kind_set)
2367 676 : nimg = dft_control%nimages
2368 676 : IF (nimg > 1) THEN
2369 232 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
2370 : END IF
2371 676 : mp_test_mode = 0
2372 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
2373 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_MULTIPOLE_TEST", debug_value, STATUS=debug_status)
2374 : IF (debug_status == 0) READ (debug_value, *, IOSTAT=debug_status) mp_test_mode
2375 : #endif
2376 :
2377 676 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
2378 :
2379 : !set up basis set lists
2380 3150 : ALLOCATE (basis_set_list(nkind))
2381 676 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
2382 :
2383 2028 : ALLOCATE (stmp(msao(tb%calc%bas%maxl)**2))
2384 2028 : ALLOCATE (dtmp(dip_n, msao(tb%calc%bas%maxl)**2))
2385 2028 : ALLOCATE (qtmp(quad_n, msao(tb%calc%bas%maxl)**2))
2386 1352 : ALLOCATE (dtmpj(dip_n, msao(tb%calc%bas%maxl)**2))
2387 1352 : ALLOCATE (qtmpj(quad_n, msao(tb%calc%bas%maxl)**2))
2388 :
2389 : ! allocate dipole/quadrupole moment matrix elemnts
2390 676 : CALL dbcsr_allocate_matrix_set(tb%dipbra, dip_n*nimg)
2391 676 : CALL dbcsr_allocate_matrix_set(tb%dipket, dip_n*nimg)
2392 676 : CALL dbcsr_allocate_matrix_set(tb%quadbra, quad_n*nimg)
2393 676 : CALL dbcsr_allocate_matrix_set(tb%quadket, quad_n*nimg)
2394 9354 : DO img = 1, nimg
2395 34712 : DO i = 1, dip_n
2396 26034 : idx = i + dip_n*(img - 1)
2397 26034 : ALLOCATE (tb%dipbra(idx)%matrix)
2398 26034 : ALLOCATE (tb%dipket(idx)%matrix)
2399 : CALL dbcsr_create(tb%dipbra(idx)%matrix, template=matrix_s(1, img)%matrix, &
2400 26034 : name="DIPOLE BRAMATRIX")
2401 : CALL dbcsr_create(tb%dipket(idx)%matrix, template=matrix_s(1, img)%matrix, &
2402 26034 : name="DIPOLE KETMATRIX")
2403 26034 : CALL cp_dbcsr_alloc_block_from_nbl(tb%dipbra(idx)%matrix, sab_orb)
2404 34712 : CALL cp_dbcsr_alloc_block_from_nbl(tb%dipket(idx)%matrix, sab_orb)
2405 : END DO
2406 61422 : DO i = 1, quad_n
2407 52068 : idx = i + quad_n*(img - 1)
2408 52068 : ALLOCATE (tb%quadbra(idx)%matrix)
2409 52068 : ALLOCATE (tb%quadket(idx)%matrix)
2410 : CALL dbcsr_create(tb%quadbra(idx)%matrix, template=matrix_s(1, img)%matrix, &
2411 52068 : name="QUADRUPOLE BRAMATRIX")
2412 : CALL dbcsr_create(tb%quadket(idx)%matrix, template=matrix_s(1, img)%matrix, &
2413 52068 : name="QUADRUPOLE KETMATRIX")
2414 52068 : CALL cp_dbcsr_alloc_block_from_nbl(tb%quadbra(idx)%matrix, sab_orb)
2415 60746 : CALL cp_dbcsr_alloc_block_from_nbl(tb%quadket(idx)%matrix, sab_orb)
2416 : END DO
2417 : END DO
2418 :
2419 : IF (mp_test_mode == 80 .OR. mp_test_mode == 82) THEN
2420 : IF (nimg > 1) CPABORT("Native tblite multipole test with k-points not implemented")
2421 : ALLOCATE (native_s(tb%calc%bas%nao, tb%calc%bas%nao))
2422 : ALLOCATE (native_h(tb%calc%bas%nao, tb%calc%bas%nao))
2423 : ALLOCATE (native_dip(dip_n, tb%calc%bas%nao, tb%calc%bas%nao))
2424 : ALLOCATE (native_quad(quad_n, tb%calc%bas%nao, tb%calc%bas%nao))
2425 : CALL get_lattice_points(tb%mol%periodic, tb%mol%lattice, get_cutoff(tb%calc%bas, 1.0_dp), native_lattr)
2426 : CALL new_adjacency_list(native_list, tb%mol, native_lattr, get_cutoff(tb%calc%bas, 1.0_dp))
2427 : CALL get_hamiltonian(tb%mol, native_lattr, native_list, tb%calc%bas, tb%calc%h0, tb%selfenergy, &
2428 : native_s, native_dip, native_quad, native_h)
2429 :
2430 : CALL dbcsr_iterator_start(iter, tb%dipbra(1)%matrix)
2431 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2432 : CALL dbcsr_iterator_next_block(iter, irow, icol, dip_bra1)
2433 : CALL dbcsr_get_block_p(matrix=tb%dipbra(2)%matrix, row=irow, col=icol, BLOCK=dip_bra2, found=found)
2434 : CPASSERT(found)
2435 : CALL dbcsr_get_block_p(matrix=tb%dipbra(3)%matrix, row=irow, col=icol, BLOCK=dip_bra3, found=found)
2436 : CPASSERT(found)
2437 : CALL dbcsr_get_block_p(matrix=tb%dipket(1)%matrix, row=irow, col=icol, BLOCK=dip_ket1, found=found)
2438 : CPASSERT(found)
2439 : CALL dbcsr_get_block_p(matrix=tb%dipket(2)%matrix, row=irow, col=icol, BLOCK=dip_ket2, found=found)
2440 : CPASSERT(found)
2441 : CALL dbcsr_get_block_p(matrix=tb%dipket(3)%matrix, row=irow, col=icol, BLOCK=dip_ket3, found=found)
2442 : CPASSERT(found)
2443 :
2444 : CALL dbcsr_get_block_p(matrix=tb%quadbra(1)%matrix, row=irow, col=icol, BLOCK=quad_bra1, found=found)
2445 : CPASSERT(found)
2446 : CALL dbcsr_get_block_p(matrix=tb%quadbra(2)%matrix, row=irow, col=icol, BLOCK=quad_bra2, found=found)
2447 : CPASSERT(found)
2448 : CALL dbcsr_get_block_p(matrix=tb%quadbra(3)%matrix, row=irow, col=icol, BLOCK=quad_bra3, found=found)
2449 : CPASSERT(found)
2450 : CALL dbcsr_get_block_p(matrix=tb%quadbra(4)%matrix, row=irow, col=icol, BLOCK=quad_bra4, found=found)
2451 : CPASSERT(found)
2452 : CALL dbcsr_get_block_p(matrix=tb%quadbra(5)%matrix, row=irow, col=icol, BLOCK=quad_bra5, found=found)
2453 : CPASSERT(found)
2454 : CALL dbcsr_get_block_p(matrix=tb%quadbra(6)%matrix, row=irow, col=icol, BLOCK=quad_bra6, found=found)
2455 : CPASSERT(found)
2456 : CALL dbcsr_get_block_p(matrix=tb%quadket(1)%matrix, row=irow, col=icol, BLOCK=quad_ket1, found=found)
2457 : CPASSERT(found)
2458 : CALL dbcsr_get_block_p(matrix=tb%quadket(2)%matrix, row=irow, col=icol, BLOCK=quad_ket2, found=found)
2459 : CPASSERT(found)
2460 : CALL dbcsr_get_block_p(matrix=tb%quadket(3)%matrix, row=irow, col=icol, BLOCK=quad_ket3, found=found)
2461 : CPASSERT(found)
2462 : CALL dbcsr_get_block_p(matrix=tb%quadket(4)%matrix, row=irow, col=icol, BLOCK=quad_ket4, found=found)
2463 : CPASSERT(found)
2464 : CALL dbcsr_get_block_p(matrix=tb%quadket(5)%matrix, row=irow, col=icol, BLOCK=quad_ket5, found=found)
2465 : CPASSERT(found)
2466 : CALL dbcsr_get_block_p(matrix=tb%quadket(6)%matrix, row=irow, col=icol, BLOCK=quad_ket6, found=found)
2467 : CPASSERT(found)
2468 :
2469 : i0 = tb%calc%bas%iao_sh(tb%calc%bas%ish_at(irow) + 1)
2470 : j0 = tb%calc%bas%iao_sh(tb%calc%bas%ish_at(icol) + 1)
2471 : DO j = 1, SIZE(dip_bra1, 1)
2472 : DO i = 1, SIZE(dip_bra1, 2)
2473 : dip_bra1(j, i) = native_dip(1, i0 + j, j0 + i)
2474 : dip_bra2(j, i) = native_dip(2, i0 + j, j0 + i)
2475 : dip_bra3(j, i) = native_dip(3, i0 + j, j0 + i)
2476 : dip_ket1(j, i) = native_dip(1, j0 + i, i0 + j)
2477 : dip_ket2(j, i) = native_dip(2, j0 + i, i0 + j)
2478 : dip_ket3(j, i) = native_dip(3, j0 + i, i0 + j)
2479 :
2480 : quad_bra1(j, i) = native_quad(1, i0 + j, j0 + i)
2481 : quad_bra2(j, i) = native_quad(2, i0 + j, j0 + i)
2482 : quad_bra3(j, i) = native_quad(3, i0 + j, j0 + i)
2483 : quad_bra4(j, i) = native_quad(4, i0 + j, j0 + i)
2484 : quad_bra5(j, i) = native_quad(5, i0 + j, j0 + i)
2485 : quad_bra6(j, i) = native_quad(6, i0 + j, j0 + i)
2486 : quad_ket1(j, i) = native_quad(1, j0 + i, i0 + j)
2487 : quad_ket2(j, i) = native_quad(2, j0 + i, i0 + j)
2488 : quad_ket3(j, i) = native_quad(3, j0 + i, i0 + j)
2489 : quad_ket4(j, i) = native_quad(4, j0 + i, i0 + j)
2490 : quad_ket5(j, i) = native_quad(5, j0 + i, i0 + j)
2491 : quad_ket6(j, i) = native_quad(6, j0 + i, i0 + j)
2492 : END DO
2493 : END DO
2494 : END DO
2495 : CALL dbcsr_iterator_stop(iter)
2496 :
2497 : DO i = 1, dip_n
2498 : CALL dbcsr_finalize(tb%dipbra(i)%matrix)
2499 : CALL dbcsr_finalize(tb%dipket(i)%matrix)
2500 : END DO
2501 : DO i = 1, quad_n
2502 : CALL dbcsr_finalize(tb%quadbra(i)%matrix)
2503 : CALL dbcsr_finalize(tb%quadket(i)%matrix)
2504 : END DO
2505 :
2506 : DEALLOCATE (native_s, native_h, native_dip, native_quad, native_lattr, basis_set_list)
2507 : CALL timestop(handle)
2508 : RETURN
2509 : END IF
2510 :
2511 : !loop over all atom pairs with a non-zero overlap (sab_orb)
2512 676 : NULLIFY (nl_iterator)
2513 676 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
2514 159898 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2515 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2516 159222 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
2517 :
2518 636888 : r2 = NORM2(rij(:))**2
2519 :
2520 159222 : icol = MAX(iatom, jatom)
2521 159222 : irow = MIN(iatom, jatom)
2522 :
2523 159222 : IF (icol == jatom) THEN
2524 366364 : rij = -rij
2525 91591 : i = ikind
2526 91591 : ikind = jkind
2527 91591 : jkind = i
2528 : END IF
2529 :
2530 159222 : IF (nimg == 1) THEN
2531 : ic = 1
2532 : ELSE
2533 18586 : ic = cell_to_index(cell(1), cell(2), cell(3))
2534 18586 : CPASSERT(ic > 0)
2535 : END IF
2536 159222 : id1 = 1 + dip_n*(ic - 1)
2537 159222 : id2 = 2 + dip_n*(ic - 1)
2538 159222 : id3 = 3 + dip_n*(ic - 1)
2539 159222 : iq1 = 1 + quad_n*(ic - 1)
2540 159222 : iq2 = 2 + quad_n*(ic - 1)
2541 159222 : iq3 = 3 + quad_n*(ic - 1)
2542 159222 : iq4 = 4 + quad_n*(ic - 1)
2543 159222 : iq5 = 5 + quad_n*(ic - 1)
2544 159222 : iq6 = 6 + quad_n*(ic - 1)
2545 :
2546 159222 : ityp = tb%mol%id(icol)
2547 159222 : jtyp = tb%mol%id(irow)
2548 :
2549 159222 : NULLIFY (dip_bra1, dip_bra2, dip_bra3)
2550 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id1)%matrix, &
2551 159222 : row=irow, col=icol, BLOCK=dip_bra1, found=found)
2552 159222 : CPASSERT(found)
2553 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id2)%matrix, &
2554 159222 : row=irow, col=icol, BLOCK=dip_bra2, found=found)
2555 159222 : CPASSERT(found)
2556 : CALL dbcsr_get_block_p(matrix=tb%dipbra(id3)%matrix, &
2557 159222 : row=irow, col=icol, BLOCK=dip_bra3, found=found)
2558 159222 : CPASSERT(found)
2559 :
2560 159222 : NULLIFY (dip_ket1, dip_ket2, dip_ket3)
2561 : CALL dbcsr_get_block_p(matrix=tb%dipket(id1)%matrix, &
2562 159222 : row=irow, col=icol, BLOCK=dip_ket1, found=found)
2563 159222 : CPASSERT(found)
2564 : CALL dbcsr_get_block_p(matrix=tb%dipket(id2)%matrix, &
2565 159222 : row=irow, col=icol, BLOCK=dip_ket2, found=found)
2566 159222 : CPASSERT(found)
2567 : CALL dbcsr_get_block_p(matrix=tb%dipket(id3)%matrix, &
2568 159222 : row=irow, col=icol, BLOCK=dip_ket3, found=found)
2569 159222 : CPASSERT(found)
2570 :
2571 159222 : NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
2572 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq1)%matrix, &
2573 159222 : row=irow, col=icol, BLOCK=quad_bra1, found=found)
2574 159222 : CPASSERT(found)
2575 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq2)%matrix, &
2576 159222 : row=irow, col=icol, BLOCK=quad_bra2, found=found)
2577 159222 : CPASSERT(found)
2578 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq3)%matrix, &
2579 159222 : row=irow, col=icol, BLOCK=quad_bra3, found=found)
2580 159222 : CPASSERT(found)
2581 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq4)%matrix, &
2582 159222 : row=irow, col=icol, BLOCK=quad_bra4, found=found)
2583 159222 : CPASSERT(found)
2584 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq5)%matrix, &
2585 159222 : row=irow, col=icol, BLOCK=quad_bra5, found=found)
2586 159222 : CPASSERT(found)
2587 : CALL dbcsr_get_block_p(matrix=tb%quadbra(iq6)%matrix, &
2588 159222 : row=irow, col=icol, BLOCK=quad_bra6, found=found)
2589 159222 : CPASSERT(found)
2590 :
2591 159222 : NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
2592 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq1)%matrix, &
2593 159222 : row=irow, col=icol, BLOCK=quad_ket1, found=found)
2594 159222 : CPASSERT(found)
2595 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq2)%matrix, &
2596 159222 : row=irow, col=icol, BLOCK=quad_ket2, found=found)
2597 159222 : CPASSERT(found)
2598 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq3)%matrix, &
2599 159222 : row=irow, col=icol, BLOCK=quad_ket3, found=found)
2600 159222 : CPASSERT(found)
2601 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq4)%matrix, &
2602 159222 : row=irow, col=icol, BLOCK=quad_ket4, found=found)
2603 159222 : CPASSERT(found)
2604 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq5)%matrix, &
2605 159222 : row=irow, col=icol, BLOCK=quad_ket5, found=found)
2606 159222 : CPASSERT(found)
2607 : CALL dbcsr_get_block_p(matrix=tb%quadket(iq6)%matrix, &
2608 159222 : row=irow, col=icol, BLOCK=quad_ket6, found=found)
2609 159222 : CPASSERT(found)
2610 :
2611 : !get basis information
2612 159222 : basis_set_a => basis_set_list(ikind)%gto_basis_set
2613 159222 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
2614 159222 : basis_set_b => basis_set_list(jkind)%gto_basis_set
2615 159222 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
2616 159222 : atom_a = atom_of_kind(icol)
2617 159222 : atom_b = atom_of_kind(irow)
2618 : ! basis a
2619 159222 : first_sgfa => basis_set_a%first_sgf
2620 159222 : la_max => basis_set_a%lmax
2621 159222 : nseta = basis_set_a%nset
2622 159222 : nsgfa => basis_set_a%nsgf_set
2623 : ! basis b
2624 159222 : first_sgfb => basis_set_b%first_sgf
2625 159222 : lb_max => basis_set_b%lmax
2626 159222 : nsetb = basis_set_b%nset
2627 159222 : nsgfb => basis_set_b%nsgf_set
2628 :
2629 : ! --------- Hamiltonian
2630 : ! Periodic self-images are off-diagonal lattice contributions, not the on-site block.
2631 159898 : IF (icol == irow .AND. r2 < same_atom**2) THEN
2632 5902 : DO iset = 1, nseta
2633 17884 : DO jset = 1, nsetb
2634 : CALL multipole_cgto(tb%calc%bas%cgto(jset, ityp), tb%calc%bas%cgto(iset, ityp), &
2635 47928 : & r2, -rij, tb%calc%bas%intcut, stmp, dtmp, qtmp)
2636 :
2637 50806 : DO inda = 1, nsgfa(iset)
2638 34544 : ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
2639 147668 : DO indb = 1, nsgfb(jset)
2640 101142 : ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
2641 101142 : ij = indb + nsgfb(jset)*(inda - 1)
2642 :
2643 101142 : dip_ket1(ib, ia) = dip_ket1(ib, ia) + dtmp(1, ij)
2644 101142 : dip_ket2(ib, ia) = dip_ket2(ib, ia) + dtmp(2, ij)
2645 101142 : dip_ket3(ib, ia) = dip_ket3(ib, ia) + dtmp(3, ij)
2646 :
2647 101142 : quad_ket1(ib, ia) = quad_ket1(ib, ia) + qtmp(1, ij)
2648 101142 : quad_ket2(ib, ia) = quad_ket2(ib, ia) + qtmp(2, ij)
2649 101142 : quad_ket3(ib, ia) = quad_ket3(ib, ia) + qtmp(3, ij)
2650 101142 : quad_ket4(ib, ia) = quad_ket4(ib, ia) + qtmp(4, ij)
2651 101142 : quad_ket5(ib, ia) = quad_ket5(ib, ia) + qtmp(5, ij)
2652 101142 : quad_ket6(ib, ia) = quad_ket6(ib, ia) + qtmp(6, ij)
2653 :
2654 101142 : dip_bra1(ib, ia) = dip_bra1(ib, ia) + dtmp(1, ij)
2655 101142 : dip_bra2(ib, ia) = dip_bra2(ib, ia) + dtmp(2, ij)
2656 101142 : dip_bra3(ib, ia) = dip_bra3(ib, ia) + dtmp(3, ij)
2657 :
2658 101142 : quad_bra1(ib, ia) = quad_bra1(ib, ia) + qtmp(1, ij)
2659 101142 : quad_bra2(ib, ia) = quad_bra2(ib, ia) + qtmp(2, ij)
2660 101142 : quad_bra3(ib, ia) = quad_bra3(ib, ia) + qtmp(3, ij)
2661 101142 : quad_bra4(ib, ia) = quad_bra4(ib, ia) + qtmp(4, ij)
2662 101142 : quad_bra5(ib, ia) = quad_bra5(ib, ia) + qtmp(5, ij)
2663 135686 : quad_bra6(ib, ia) = quad_bra6(ib, ia) + qtmp(6, ij)
2664 : END DO
2665 : END DO
2666 : END DO
2667 : END DO
2668 : ELSE
2669 629434 : DO iset = 1, nseta
2670 2043798 : DO jset = 1, nsetb
2671 5657456 : mp_int_vec = -rij
2672 5657456 : mp_shift_vec = -rij
2673 : IF (mp_test_mode == 14) THEN
2674 : mp_int_vec = rij
2675 : mp_shift_vec = rij
2676 : ELSEIF (mp_test_mode == 15) THEN
2677 : mp_shift_vec = rij
2678 : ELSEIF (mp_test_mode == 16) THEN
2679 : mp_int_vec = rij
2680 : END IF
2681 : CALL multipole_cgto(tb%calc%bas%cgto(jset, jtyp), tb%calc%bas%cgto(iset, ityp), &
2682 1414364 : & r2, mp_int_vec, tb%calc%bas%intcut, stmp, dtmp, qtmp)
2683 :
2684 6127156 : DO inda = 1, nsgfa(iset)
2685 4240958 : ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
2686 18374148 : DO indb = 1, nsgfb(jset)
2687 12718826 : ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
2688 :
2689 12718826 : ij = indb + nsgfb(jset)*(inda - 1)
2690 : CALL tb_shift_multipole(mp_shift_vec, stmp(ij), dtmp(:, ij), qtmp(:, ij), &
2691 12718826 : dtmpj(:, ij), qtmpj(:, ij))
2692 12718826 : IF (icol == irow .AND. r2 >= same_atom**2) THEN
2693 3282768 : IF (nimg > 1) THEN
2694 : ! Match the k-point image block orientation used for self-image quadrupoles.
2695 3882816 : qtmp(:, ij) = -qtmp(:, ij)
2696 3882816 : qtmpj(:, ij) = -qtmpj(:, ij)
2697 : END IF
2698 : IF (mp_test_mode == 30 .OR. mp_test_mode == 32 .OR. mp_test_mode == 33) THEN
2699 : dtmpj(:, ij) = dtmp(:, ij)
2700 : qtmpj(:, ij) = qtmp(:, ij)
2701 : ELSEIF (mp_test_mode == 31) THEN
2702 : CALL tb_shift_multipole(-mp_shift_vec, stmp(ij), dtmp(:, ij), qtmp(:, ij), &
2703 : dtmpj(:, ij), qtmpj(:, ij))
2704 : END IF
2705 : IF (mp_test_mode == 40 .OR. mp_test_mode == 41 .OR. mp_test_mode == 44) THEN
2706 : dtmp(:, ij) = -dtmp(:, ij)
2707 : END IF
2708 : IF (mp_test_mode == 40 .OR. mp_test_mode == 42 .OR. mp_test_mode == 44) THEN
2709 : dtmpj(:, ij) = -dtmpj(:, ij)
2710 : END IF
2711 : IF (mp_test_mode == 40 .OR. mp_test_mode == 41 .OR. mp_test_mode == 45) THEN
2712 : qtmp(:, ij) = -qtmp(:, ij)
2713 : END IF
2714 : IF (mp_test_mode == 40 .OR. mp_test_mode == 42 .OR. mp_test_mode == 45) THEN
2715 : qtmpj(:, ij) = -qtmpj(:, ij)
2716 : END IF
2717 : IF (mp_test_mode == 43) THEN
2718 : dtmp(:, ij) = 0.0_dp
2719 : qtmp(:, ij) = 0.0_dp
2720 : dtmpj(:, ij) = 0.0_dp
2721 : qtmpj(:, ij) = 0.0_dp
2722 : END IF
2723 : END IF
2724 : IF (r2 >= same_atom**2) THEN
2725 : IF (mp_test_mode == 46 .OR. mp_test_mode == 48) THEN
2726 : dtmp(:, ij) = -dtmp(:, ij)
2727 : dtmpj(:, ij) = -dtmpj(:, ij)
2728 : END IF
2729 : IF (mp_test_mode == 46 .OR. mp_test_mode == 49) THEN
2730 : qtmp(:, ij) = -qtmp(:, ij)
2731 : qtmpj(:, ij) = -qtmpj(:, ij)
2732 : END IF
2733 : IF (mp_test_mode == 47) THEN
2734 : dtmp(:, ij) = 0.0_dp
2735 : qtmp(:, ij) = 0.0_dp
2736 : dtmpj(:, ij) = 0.0_dp
2737 : qtmpj(:, ij) = 0.0_dp
2738 : END IF
2739 : IF (mp_test_mode == 50) THEN
2740 : dtmp(:, ij) = 0.25_dp*dtmp(:, ij)
2741 : qtmp(:, ij) = 0.25_dp*qtmp(:, ij)
2742 : dtmpj(:, ij) = 0.25_dp*dtmpj(:, ij)
2743 : qtmpj(:, ij) = 0.25_dp*qtmpj(:, ij)
2744 : ELSEIF (mp_test_mode == 51) THEN
2745 : dtmp(:, ij) = 0.5_dp*dtmp(:, ij)
2746 : qtmp(:, ij) = 0.5_dp*qtmp(:, ij)
2747 : dtmpj(:, ij) = 0.5_dp*dtmpj(:, ij)
2748 : qtmpj(:, ij) = 0.5_dp*qtmpj(:, ij)
2749 : ELSEIF (mp_test_mode == 52) THEN
2750 : dtmp(:, ij) = 0.75_dp*dtmp(:, ij)
2751 : qtmp(:, ij) = 0.75_dp*qtmp(:, ij)
2752 : dtmpj(:, ij) = 0.75_dp*dtmpj(:, ij)
2753 : qtmpj(:, ij) = 0.75_dp*qtmpj(:, ij)
2754 : ELSEIF (mp_test_mode == 53) THEN
2755 : dtmp(:, ij) = -0.25_dp*dtmp(:, ij)
2756 : qtmp(:, ij) = -0.25_dp*qtmp(:, ij)
2757 : dtmpj(:, ij) = -0.25_dp*dtmpj(:, ij)
2758 : qtmpj(:, ij) = -0.25_dp*qtmpj(:, ij)
2759 : ELSEIF (mp_test_mode == 54) THEN
2760 : dtmp(:, ij) = -0.5_dp*dtmp(:, ij)
2761 : qtmp(:, ij) = -0.5_dp*qtmp(:, ij)
2762 : dtmpj(:, ij) = -0.5_dp*dtmpj(:, ij)
2763 : qtmpj(:, ij) = -0.5_dp*qtmpj(:, ij)
2764 : ELSEIF (mp_test_mode == 55) THEN
2765 : dtmp(:, ij) = -0.75_dp*dtmp(:, ij)
2766 : qtmp(:, ij) = -0.75_dp*qtmp(:, ij)
2767 : dtmpj(:, ij) = -0.75_dp*dtmpj(:, ij)
2768 : qtmpj(:, ij) = -0.75_dp*qtmpj(:, ij)
2769 : END IF
2770 : END IF
2771 :
2772 12718826 : dip_bra1(ib, ia) = dip_bra1(ib, ia) + dtmp(1, ij)
2773 12718826 : dip_bra2(ib, ia) = dip_bra2(ib, ia) + dtmp(2, ij)
2774 12718826 : dip_bra3(ib, ia) = dip_bra3(ib, ia) + dtmp(3, ij)
2775 :
2776 12718826 : quad_bra1(ib, ia) = quad_bra1(ib, ia) + qtmp(1, ij)
2777 12718826 : quad_bra2(ib, ia) = quad_bra2(ib, ia) + qtmp(2, ij)
2778 12718826 : quad_bra3(ib, ia) = quad_bra3(ib, ia) + qtmp(3, ij)
2779 12718826 : quad_bra4(ib, ia) = quad_bra4(ib, ia) + qtmp(4, ij)
2780 12718826 : quad_bra5(ib, ia) = quad_bra5(ib, ia) + qtmp(5, ij)
2781 12718826 : quad_bra6(ib, ia) = quad_bra6(ib, ia) + qtmp(6, ij)
2782 :
2783 12718826 : dip_ket1(ib, ia) = dip_ket1(ib, ia) + dtmpj(1, ij)
2784 12718826 : dip_ket2(ib, ia) = dip_ket2(ib, ia) + dtmpj(2, ij)
2785 12718826 : dip_ket3(ib, ia) = dip_ket3(ib, ia) + dtmpj(3, ij)
2786 :
2787 12718826 : quad_ket1(ib, ia) = quad_ket1(ib, ia) + qtmpj(1, ij)
2788 12718826 : quad_ket2(ib, ia) = quad_ket2(ib, ia) + qtmpj(2, ij)
2789 12718826 : quad_ket3(ib, ia) = quad_ket3(ib, ia) + qtmpj(3, ij)
2790 12718826 : quad_ket4(ib, ia) = quad_ket4(ib, ia) + qtmpj(4, ij)
2791 12718826 : quad_ket5(ib, ia) = quad_ket5(ib, ia) + qtmpj(5, ij)
2792 16959784 : quad_ket6(ib, ia) = quad_ket6(ib, ia) + qtmpj(6, ij)
2793 : END DO
2794 : END DO
2795 : END DO
2796 : END DO
2797 : END IF
2798 : END DO
2799 676 : CALL neighbor_list_iterator_release(nl_iterator)
2800 :
2801 26710 : DO i = 1, SIZE(tb%dipbra)
2802 26034 : CALL dbcsr_finalize(tb%dipbra(i)%matrix)
2803 26710 : CALL dbcsr_finalize(tb%dipket(i)%matrix)
2804 : END DO
2805 52744 : DO i = 1, SIZE(tb%quadbra)
2806 52068 : CALL dbcsr_finalize(tb%quadbra(i)%matrix)
2807 52744 : CALL dbcsr_finalize(tb%quadket(i)%matrix)
2808 : END DO
2809 :
2810 676 : DEALLOCATE (basis_set_list)
2811 :
2812 676 : CALL timestop(handle)
2813 :
2814 : #else
2815 : MARK_USED(qs_env)
2816 : MARK_USED(tb)
2817 : CPABORT("Built without TBLITE")
2818 : #endif
2819 :
2820 1352 : END SUBROUTINE tb_get_multipole
2821 :
2822 : ! **************************************************************************************************
2823 : !> \brief Shift a multipole operator from one center to the other.
2824 : !> \param vec displacement vector between the two centers
2825 : !> \param s overlap integral
2826 : !> \param di dipole integral on the original center
2827 : !> \param qi quadrupole integral on the original center
2828 : !> \param dj dipole integral on the shifted center
2829 : !> \param qj quadrupole integral on the shifted center
2830 : ! **************************************************************************************************
2831 12718826 : PURE SUBROUTINE tb_shift_multipole(vec, s, di, qi, dj, qj)
2832 :
2833 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vec
2834 : REAL(KIND=dp), INTENT(IN) :: s
2835 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: di, qi
2836 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: dj, qj
2837 :
2838 : REAL(KIND=dp) :: tr
2839 :
2840 12718826 : dj(1) = di(1) + vec(1)*s
2841 12718826 : dj(2) = di(2) + vec(2)*s
2842 12718826 : dj(3) = di(3) + vec(3)*s
2843 :
2844 12718826 : qj(1) = 2*vec(1)*di(1) + vec(1)**2*s
2845 12718826 : qj(3) = 2*vec(2)*di(2) + vec(2)**2*s
2846 12718826 : qj(6) = 2*vec(3)*di(3) + vec(3)**2*s
2847 12718826 : qj(2) = vec(1)*di(2) + vec(2)*di(1) + vec(1)*vec(2)*s
2848 12718826 : qj(4) = vec(1)*di(3) + vec(3)*di(1) + vec(1)*vec(3)*s
2849 12718826 : qj(5) = vec(2)*di(3) + vec(3)*di(2) + vec(2)*vec(3)*s
2850 12718826 : tr = 0.5_dp*(qj(1) + qj(3) + qj(6))
2851 :
2852 12718826 : qj(1) = qi(1) + 1.5_dp*qj(1) - tr
2853 12718826 : qj(2) = qi(2) + 1.5_dp*qj(2)
2854 12718826 : qj(3) = qi(3) + 1.5_dp*qj(3) - tr
2855 12718826 : qj(4) = qi(4) + 1.5_dp*qj(4)
2856 12718826 : qj(5) = qi(5) + 1.5_dp*qj(5)
2857 12718826 : qj(6) = qi(6) + 1.5_dp*qj(6) - tr
2858 :
2859 12718826 : END SUBROUTINE tb_shift_multipole
2860 :
2861 : ! **************************************************************************************************
2862 : !> \brief compute the mulliken properties (AO resolved)
2863 : !> \param p_matrix ...
2864 : !> \param bra_mat ...
2865 : !> \param ket_mat ...
2866 : !> \param output ...
2867 : !> \param para_env ...
2868 : !> \par History
2869 : !> adapted from ao_charges_2
2870 : !> \note
2871 : ! **************************************************************************************************
2872 1886886 : SUBROUTINE contract_dens(p_matrix, bra_mat, ket_mat, output, para_env)
2873 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
2874 : TYPE(dbcsr_type), POINTER :: bra_mat, ket_mat
2875 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: output
2876 : TYPE(mp_para_env_type), POINTER :: para_env
2877 :
2878 : CHARACTER(len=*), PARAMETER :: routineN = 'contract_dens'
2879 :
2880 : INTEGER :: handle, i, iblock_col, &
2881 : iblock_row, ispin, j, mp_test_mode, &
2882 : nspin
2883 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
2884 : INTEGER :: debug_status
2885 : CHARACTER(LEN=32) :: debug_value
2886 : #endif
2887 : LOGICAL :: found
2888 1886886 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: bra, ket, p_block
2889 : TYPE(dbcsr_iterator_type) :: iter
2890 :
2891 1886886 : CALL timeset(routineN, handle)
2892 :
2893 1886886 : mp_test_mode = 0
2894 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
2895 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_MULTIPOLE_TEST", debug_value, STATUS=debug_status)
2896 : IF (debug_status == 0) READ (debug_value, *, IOSTAT=debug_status) mp_test_mode
2897 : #endif
2898 :
2899 1886886 : nspin = SIZE(p_matrix)
2900 9918288 : output = 0.0_dp
2901 3773772 : DO ispin = 1, nspin
2902 1886886 : CALL dbcsr_iterator_start(iter, bra_mat)
2903 12903030 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2904 11016144 : NULLIFY (p_block, bra, ket)
2905 :
2906 11016144 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, bra)
2907 : CALL dbcsr_get_block_p(matrix=p_matrix(ispin)%matrix, &
2908 11016144 : row=iblock_row, col=iblock_col, BLOCK=p_block, found=found)
2909 11016144 : IF (.NOT. found) CYCLE
2910 : CALL dbcsr_get_block_p(matrix=ket_mat, &
2911 11016144 : row=iblock_row, col=iblock_col, BLOCK=ket, found=found)
2912 11016144 : IF (.NOT. found) CPABORT("missing block")
2913 :
2914 11016144 : IF (.NOT. (ASSOCIATED(bra) .AND. ASSOCIATED(p_block))) CYCLE
2915 12903030 : IF (iblock_col == iblock_row) THEN
2916 : IF (mp_test_mode == 70) THEN
2917 : DO j = 1, SIZE(p_block, 1)
2918 : DO i = 1, SIZE(p_block, 2)
2919 : output(iblock_row) = output(iblock_row) + p_block(j, i)*ket(j, i)
2920 : END DO
2921 : END DO
2922 : ELSEIF (mp_test_mode == 71) THEN
2923 : DO j = 1, SIZE(p_block, 1)
2924 : DO i = 1, SIZE(p_block, 2)
2925 : output(iblock_row) = output(iblock_row) + 0.5_dp*p_block(j, i)*(bra(j, i) + ket(j, i))
2926 : END DO
2927 : END DO
2928 : ELSE
2929 39706083 : DO j = 1, SIZE(p_block, 1)
2930 359758629 : DO i = 1, SIZE(p_block, 2)
2931 355742928 : output(iblock_row) = output(iblock_row) + p_block(j, i)*bra(j, i)
2932 : END DO
2933 : END DO
2934 : END IF
2935 : ELSE
2936 69381711 : DO j = 1, SIZE(p_block, 1)
2937 628284645 : DO i = 1, SIZE(p_block, 2)
2938 621284202 : output(iblock_row) = output(iblock_row) + p_block(j, i)*ket(j, i)
2939 : END DO
2940 : END DO
2941 69381711 : DO j = 1, SIZE(p_block, 1)
2942 628284645 : DO i = 1, SIZE(p_block, 2)
2943 621284202 : output(iblock_col) = output(iblock_col) + p_block(j, i)*bra(j, i)
2944 : END DO
2945 : END DO
2946 : END IF
2947 : END DO
2948 5660658 : CALL dbcsr_iterator_stop(iter)
2949 : END DO
2950 :
2951 17949690 : CALL para_env%sum(output)
2952 1886886 : CALL timestop(handle)
2953 :
2954 1886886 : END SUBROUTINE contract_dens
2955 :
2956 : ! **************************************************************************************************
2957 : !> \brief compute the AO-resolved density contraction for real-space k-point image matrices
2958 : !> \param p_matrix ...
2959 : !> \param bra_mat ...
2960 : !> \param ket_mat ...
2961 : !> \param iop ...
2962 : !> \param nops ...
2963 : !> \param output ...
2964 : !> \param para_env ...
2965 : ! **************************************************************************************************
2966 57384 : SUBROUTINE contract_dens_kp(p_matrix, bra_mat, ket_mat, iop, nops, output, para_env)
2967 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_matrix
2968 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: bra_mat, ket_mat
2969 : INTEGER, INTENT(IN) :: iop, nops
2970 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: output
2971 : TYPE(mp_para_env_type), POINTER :: para_env
2972 :
2973 : INTEGER :: ic, idx, nimg
2974 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: image_output
2975 57384 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_image
2976 :
2977 57384 : nimg = SIZE(p_matrix, 2)
2978 288072 : output = 0.0_dp
2979 172152 : ALLOCATE (image_output(SIZE(output)))
2980 1803492 : DO ic = 1, nimg
2981 1746108 : idx = iop + nops*(ic - 1)
2982 1746108 : CPASSERT(idx <= SIZE(bra_mat))
2983 1746108 : CPASSERT(idx <= SIZE(ket_mat))
2984 : NULLIFY (p_image)
2985 1746108 : p_image => p_matrix(:, ic)
2986 8871084 : image_output = 0.0_dp
2987 1746108 : CALL contract_dens(p_image, bra_mat(idx)%matrix, ket_mat(idx)%matrix, image_output, para_env)
2988 8928468 : output = output + image_output
2989 : END DO
2990 57384 : DEALLOCATE (image_output)
2991 :
2992 57384 : END SUBROUTINE contract_dens_kp
2993 :
2994 : ! **************************************************************************************************
2995 : !> \brief save gradient to force
2996 : !> \param qs_env ...
2997 : !> \param tb ...
2998 : !> \param para_env ...
2999 : !> \param ityp ...
3000 : !> \note
3001 : ! **************************************************************************************************
3002 350 : SUBROUTINE tb_grad2force(qs_env, tb, para_env, ityp)
3003 :
3004 : TYPE(qs_environment_type) :: qs_env
3005 : TYPE(tblite_type) :: tb
3006 : TYPE(mp_para_env_type) :: para_env
3007 : INTEGER :: ityp
3008 :
3009 : CHARACTER(len=*), PARAMETER :: routineN = 'tb_grad2force'
3010 :
3011 : CHARACTER(LEN=default_path_length) :: dump_file
3012 : INTEGER :: atoma, dump_status, dump_unit, handle, &
3013 : iatom, ikind, natom
3014 350 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
3015 350 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3016 350 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3017 350 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3018 :
3019 350 : CALL timeset(routineN, handle)
3020 :
3021 350 : NULLIFY (force, atomic_kind_set)
3022 : CALL get_qs_env(qs_env=qs_env, force=force, particle_set=particle_set, &
3023 350 : atomic_kind_set=atomic_kind_set)
3024 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
3025 350 : atom_of_kind=atom_of_kind, kind_of=kind_of)
3026 :
3027 350 : natom = SIZE(particle_set)
3028 :
3029 350 : dump_status = 1
3030 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
3031 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_FORCE_DUMP", dump_file, STATUS=dump_status)
3032 : #endif
3033 : IF (dump_status == 0) THEN
3034 : OPEN (NEWUNIT=dump_unit, FILE=TRIM(dump_file), STATUS="UNKNOWN", &
3035 : POSITION="APPEND", ACTION="WRITE")
3036 : WRITE (dump_unit, "(A,1X,I0)") "component", ityp
3037 : DO iatom = 1, natom
3038 : WRITE (dump_unit, "(I0,3(1X,ES24.16))") iatom, tb%grad(:, iatom)/para_env%num_pe
3039 : END DO
3040 : CLOSE (dump_unit)
3041 : END IF
3042 :
3043 350 : SELECT CASE (ityp)
3044 : CASE DEFAULT
3045 0 : CPABORT("unknown force type")
3046 : CASE (0)
3047 146 : DO iatom = 1, natom
3048 120 : ikind = kind_of(iatom)
3049 120 : atoma = atom_of_kind(iatom)
3050 : force(ikind)%all_potential(:, atoma) = &
3051 506 : force(ikind)%all_potential(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3052 : END DO
3053 : CASE (1)
3054 300 : DO iatom = 1, natom
3055 246 : ikind = kind_of(iatom)
3056 246 : atoma = atom_of_kind(iatom)
3057 : force(ikind)%repulsive(:, atoma) = &
3058 1038 : force(ikind)%repulsive(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3059 : END DO
3060 : CASE (2)
3061 600 : DO iatom = 1, natom
3062 492 : ikind = kind_of(iatom)
3063 492 : atoma = atom_of_kind(iatom)
3064 : force(ikind)%dispersion(:, atoma) = &
3065 2076 : force(ikind)%dispersion(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3066 : END DO
3067 : CASE (3)
3068 300 : DO iatom = 1, natom
3069 246 : ikind = kind_of(iatom)
3070 246 : atoma = atom_of_kind(iatom)
3071 : force(ikind)%rho_elec(:, atoma) = &
3072 1038 : force(ikind)%rho_elec(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3073 : END DO
3074 : CASE (4)
3075 600 : DO iatom = 1, natom
3076 492 : ikind = kind_of(iatom)
3077 492 : atoma = atom_of_kind(iatom)
3078 : force(ikind)%overlap(:, atoma) = &
3079 2076 : force(ikind)%overlap(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3080 : END DO
3081 : CASE (5)
3082 350 : DO iatom = 1, natom
3083 0 : ikind = kind_of(iatom)
3084 0 : atoma = atom_of_kind(iatom)
3085 : force(ikind)%efield(:, atoma) = &
3086 0 : force(ikind)%efield(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
3087 : END DO
3088 : END SELECT
3089 :
3090 350 : CALL timestop(handle)
3091 :
3092 700 : END SUBROUTINE tb_grad2force
3093 :
3094 : ! **************************************************************************************************
3095 : !> \brief set gradient to zero
3096 : !> \param qs_env ...
3097 : !> \note
3098 : ! **************************************************************************************************
3099 54 : SUBROUTINE tb_zero_force(qs_env)
3100 :
3101 : TYPE(qs_environment_type) :: qs_env
3102 :
3103 : INTEGER :: iatom, ikind, natom
3104 54 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
3105 54 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3106 54 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3107 54 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3108 :
3109 54 : NULLIFY (force, atomic_kind_set)
3110 : CALL get_qs_env(qs_env=qs_env, force=force, particle_set=particle_set, &
3111 54 : atomic_kind_set=atomic_kind_set)
3112 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
3113 54 : kind_of=kind_of)
3114 :
3115 54 : natom = SIZE(particle_set)
3116 :
3117 300 : DO iatom = 1, natom
3118 246 : ikind = kind_of(iatom)
3119 4254 : force(ikind)%all_potential = 0.0_dp
3120 4254 : force(ikind)%repulsive = 0.0_dp
3121 4254 : force(ikind)%dispersion = 0.0_dp
3122 4254 : force(ikind)%rho_elec = 0.0_dp
3123 4254 : force(ikind)%overlap = 0.0_dp
3124 4308 : force(ikind)%efield = 0.0_dp
3125 : END DO
3126 :
3127 108 : END SUBROUTINE tb_zero_force
3128 :
3129 : ! **************************************************************************************************
3130 : !> \brief compute the derivative of the energy
3131 : !> \param qs_env ...
3132 : !> \param use_rho ...
3133 : !> \param nimg ...
3134 : ! **************************************************************************************************
3135 54 : SUBROUTINE tb_derive_dH_diag(qs_env, use_rho, nimg)
3136 :
3137 : TYPE(qs_environment_type), POINTER :: qs_env
3138 : LOGICAL, INTENT(IN) :: use_rho
3139 : INTEGER, INTENT(IN) :: nimg
3140 :
3141 : #if defined(__TBLITE)
3142 : INTEGER :: i, iatom, ic, ikind, ind, is, ish, &
3143 : jatom, jkind
3144 : INTEGER, DIMENSION(3) :: cellind
3145 54 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3146 : LOGICAL :: found
3147 54 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dE
3148 54 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock
3149 54 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
3150 : TYPE(kpoint_type), POINTER :: kpoints
3151 : TYPE(mp_para_env_type), POINTER :: para_env
3152 : TYPE(neighbor_list_iterator_p_type), &
3153 54 : DIMENSION(:), POINTER :: nl_iterator
3154 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3155 54 : POINTER :: sab_orb
3156 : TYPE(qs_rho_type), POINTER :: rho
3157 : TYPE(qs_scf_env_type), POINTER :: scf_env
3158 : TYPE(tblite_type), POINTER :: tb
3159 :
3160 : ! compute mulliken charges required for charge update
3161 54 : NULLIFY (scf_env, rho, tb, sab_orb, para_env, kpoints)
3162 : CALL get_qs_env(qs_env=qs_env, scf_env=scf_env, rho=rho, tb_tblite=tb, &
3163 54 : sab_orb=sab_orb, para_env=para_env, kpoints=kpoints)
3164 :
3165 54 : NULLIFY (cell_to_index)
3166 54 : IF (nimg > 1) THEN
3167 18 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
3168 : END IF
3169 :
3170 54 : NULLIFY (matrix_p)
3171 54 : IF (use_rho) THEN
3172 54 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
3173 : ELSE
3174 0 : matrix_p => scf_env%p_mix_new
3175 : END IF
3176 :
3177 162 : ALLOCATE (dE(tb%mol%nat))
3178 :
3179 300 : dE = 0.0_dp
3180 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
3181 54 : NULLIFY (nl_iterator)
3182 54 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
3183 12358 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3184 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3185 12304 : iatom=iatom, jatom=jatom, cell=cellind)
3186 :
3187 12304 : IF (iatom /= jatom) CYCLE
3188 :
3189 3397 : IF (ikind /= jkind) CPABORT("Type wrong")
3190 :
3191 3397 : is = tb%calc%bas%ish_at(iatom)
3192 :
3193 3397 : IF (nimg == 1) THEN
3194 : ic = 1
3195 : ELSE
3196 732 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
3197 732 : CPASSERT(ic > 0)
3198 : END IF
3199 :
3200 3397 : IF (cellind(1) /= 0) CYCLE
3201 1127 : IF (cellind(2) /= 0) CYCLE
3202 353 : IF (cellind(3) /= 0) CYCLE
3203 :
3204 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
3205 123 : row=iatom, col=jatom, BLOCK=pblock, found=found)
3206 :
3207 123 : IF (.NOT. found) CPABORT("block not found")
3208 :
3209 123 : ind = 0
3210 496 : DO ish = 1, tb%calc%bas%nsh_id(ikind)
3211 13468 : DO i = 1, msao(tb%calc%bas%cgto(ish, ikind)%ang)
3212 845 : ind = ind + 1
3213 1164 : dE(iatom) = dE(iatom) + tb%dsedcn(is + ish)*pblock(ind, ind)
3214 : END DO
3215 : END DO
3216 : END DO
3217 54 : CALL neighbor_list_iterator_release(nl_iterator)
3218 54 : CALL para_env%sum(dE)
3219 :
3220 1038 : tb%grad = 0.0_dp
3221 54 : CALL tb_add_grad(tb%grad, tb%dcndr, dE, tb%mol%nat)
3222 54 : IF (tb%use_virial) CALL tb_add_sig(tb%sigma, tb%dcndL, dE, tb%mol%nat)
3223 54 : CALL tb_grad2force(qs_env, tb, para_env, 4)
3224 54 : DEALLOCATE (dE)
3225 :
3226 : #else
3227 : MARK_USED(qs_env)
3228 : MARK_USED(use_rho)
3229 : MARK_USED(nimg)
3230 : CPABORT("Built without TBLITE")
3231 : #endif
3232 :
3233 54 : END SUBROUTINE tb_derive_dH_diag
3234 :
3235 : ! **************************************************************************************************
3236 : !> \brief compute the derivative of the energy
3237 : !> \param qs_env ...
3238 : !> \param use_rho ...
3239 : !> \param nimg ...
3240 : ! **************************************************************************************************
3241 54 : SUBROUTINE tb_derive_dH_off(qs_env, use_rho, nimg)
3242 :
3243 : TYPE(qs_environment_type), POINTER :: qs_env
3244 : LOGICAL, INTENT(IN) :: use_rho
3245 : INTEGER, INTENT(IN) :: nimg
3246 :
3247 : #if defined(__TBLITE)
3248 : INTEGER :: i, ij, iatom, ic, icol, ikind, ni, nj, nkind, nel, &
3249 : sgfi, sgfj, ityp, jatom, jkind, jrow, jtyp, iset, jset, &
3250 : nseti, nsetj, ia, ib, inda, indb
3251 : INTEGER, DIMENSION(3) :: cellind
3252 54 : INTEGER, DIMENSION(:), POINTER :: nsgfa, nsgfb
3253 54 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
3254 54 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3255 : LOGICAL :: found, image_pair
3256 : REAL(KIND=dp) :: r2, dr, itemp, jtemp, rr, hij, shpoly, dshpoly, idHdc, jdHdc, &
3257 : scal, hp, i_a_shift, j_a_shift, ishift, jshift, &
3258 : dip_deriv_fac, quad_deriv_fac, dip_deriv_fac_pair, &
3259 : quad_deriv_fac_pair, overlap_deriv_fac, pot_deriv_scale, &
3260 : w_deriv_scale, h0_deriv_scale, mp_deriv_scale, &
3261 : pot_image_deriv_scale, w_image_deriv_scale, &
3262 : h0_image_deriv_scale, mp_image_deriv_scale, &
3263 : pot_pair_scale, w_pair_scale, h0_pair_scale, mp_pair_scale, &
3264 : pot_self_image_deriv_scale, w_self_image_deriv_scale, &
3265 : h0_self_image_deriv_scale, mp_self_image_deriv_scale, &
3266 : cn_image_scale
3267 : INTEGER :: mp_test_mode
3268 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
3269 : INTEGER :: debug_status
3270 : CHARACTER(LEN=32) :: debug_value, scale_value
3271 : #endif
3272 : REAL(KIND=dp), DIMENSION(3) :: rij, dgrad
3273 : REAL(KIND=dp), DIMENSION(3, 3) :: hsigma
3274 54 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dE, t_ov, idip, jdip, iquad, jquad
3275 54 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: t_dip, t_quad, t_d_ov
3276 54 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: t_i_dip, t_i_quad, t_j_dip, t_j_quad
3277 54 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock, wblock
3278 :
3279 54 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3280 54 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_w
3281 54 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
3282 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
3283 : TYPE(kpoint_type), POINTER :: kpoints
3284 : TYPE(mp_para_env_type), POINTER :: para_env
3285 : TYPE(neighbor_list_iterator_p_type), &
3286 54 : DIMENSION(:), POINTER :: nl_iterator
3287 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3288 54 : POINTER :: sab_orb
3289 54 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3290 : TYPE(qs_rho_type), POINTER :: rho
3291 : TYPE(qs_scf_env_type), POINTER :: scf_env
3292 : TYPE(tb_hamiltonian), POINTER :: h0
3293 : TYPE(tblite_type), POINTER :: tb
3294 :
3295 : ! compute mulliken charges required for charge update
3296 54 : NULLIFY (scf_env, rho, tb, sab_orb, para_env, kpoints, matrix_w)
3297 : CALL get_qs_env(qs_env=qs_env, &
3298 : atomic_kind_set=atomic_kind_set, &
3299 : scf_env=scf_env, &
3300 : rho=rho, &
3301 : tb_tblite=tb, &
3302 : sab_orb=sab_orb, &
3303 : para_env=para_env, &
3304 : qs_kind_set=qs_kind_set, &
3305 : kpoints=kpoints, &
3306 54 : matrix_w_kp=matrix_w)
3307 :
3308 54 : NULLIFY (cell_to_index)
3309 54 : IF (nimg > 1) THEN
3310 18 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
3311 : END IF
3312 :
3313 54 : h0 => tb%calc%h0
3314 54 : mp_test_mode = 0
3315 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
3316 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_MULTIPOLE_TEST", debug_value, STATUS=debug_status)
3317 : IF (debug_status == 0) READ (debug_value, *, IOSTAT=debug_status) mp_test_mode
3318 : #endif
3319 : dip_deriv_fac = -1.0_dp
3320 : quad_deriv_fac = -1.0_dp
3321 : IF (mp_test_mode == 7 .OR. mp_test_mode == 11) THEN
3322 : dip_deriv_fac = 1.0_dp
3323 : quad_deriv_fac = 1.0_dp
3324 : END IF
3325 : IF (mp_test_mode == 8 .OR. mp_test_mode == 12) THEN
3326 : dip_deriv_fac = 0.0_dp
3327 : quad_deriv_fac = 0.0_dp
3328 : END IF
3329 : IF (mp_test_mode == 13 .OR. mp_test_mode == 14 .OR. mp_test_mode == 15 .OR. &
3330 : mp_test_mode == 16 .OR. mp_test_mode == 21 .OR. mp_test_mode == 33 .OR. &
3331 : mp_test_mode == 34 .OR. mp_test_mode == 37) THEN
3332 : dip_deriv_fac = -0.5_dp
3333 : quad_deriv_fac = -0.5_dp
3334 : END IF
3335 : IF (mp_test_mode == 17) THEN
3336 : dip_deriv_fac = -1.0_dp
3337 : quad_deriv_fac = 0.0_dp
3338 : END IF
3339 : IF (mp_test_mode == 18) THEN
3340 : dip_deriv_fac = 0.0_dp
3341 : quad_deriv_fac = -1.0_dp
3342 : END IF
3343 54 : pot_deriv_scale = 1.0_dp
3344 54 : w_deriv_scale = 1.0_dp
3345 54 : h0_deriv_scale = 1.0_dp
3346 54 : mp_deriv_scale = 1.0_dp
3347 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
3348 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_POT_SCALE", scale_value, STATUS=debug_status)
3349 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) pot_deriv_scale
3350 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_W_SCALE", scale_value, STATUS=debug_status)
3351 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) w_deriv_scale
3352 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_H0_SCALE", scale_value, STATUS=debug_status)
3353 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) h0_deriv_scale
3354 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_MP_SCALE", scale_value, STATUS=debug_status)
3355 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) mp_deriv_scale
3356 : #endif
3357 54 : pot_image_deriv_scale = pot_deriv_scale
3358 54 : w_image_deriv_scale = w_deriv_scale
3359 54 : h0_image_deriv_scale = h0_deriv_scale
3360 54 : mp_image_deriv_scale = mp_deriv_scale
3361 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
3362 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_POT_IMAGE_SCALE", scale_value, STATUS=debug_status)
3363 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) pot_image_deriv_scale
3364 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_W_IMAGE_SCALE", scale_value, STATUS=debug_status)
3365 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) w_image_deriv_scale
3366 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_H0_IMAGE_SCALE", scale_value, STATUS=debug_status)
3367 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) h0_image_deriv_scale
3368 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_MP_IMAGE_SCALE", scale_value, STATUS=debug_status)
3369 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) mp_image_deriv_scale
3370 : #endif
3371 54 : pot_self_image_deriv_scale = pot_image_deriv_scale
3372 54 : w_self_image_deriv_scale = w_image_deriv_scale
3373 54 : h0_self_image_deriv_scale = h0_image_deriv_scale
3374 54 : mp_self_image_deriv_scale = mp_image_deriv_scale
3375 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
3376 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_POT_SELF_IMAGE_SCALE", scale_value, STATUS=debug_status)
3377 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) pot_self_image_deriv_scale
3378 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_W_SELF_IMAGE_SCALE", scale_value, STATUS=debug_status)
3379 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) w_self_image_deriv_scale
3380 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_H0_SELF_IMAGE_SCALE", scale_value, STATUS=debug_status)
3381 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) h0_self_image_deriv_scale
3382 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_MP_SELF_IMAGE_SCALE", scale_value, STATUS=debug_status)
3383 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) mp_self_image_deriv_scale
3384 : #endif
3385 54 : cn_image_scale = 1.0_dp
3386 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
3387 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_DH_CN_IMAGE_SCALE", scale_value, STATUS=debug_status)
3388 : IF (debug_status == 0) READ (scale_value, *, IOSTAT=debug_status) cn_image_scale
3389 : #endif
3390 :
3391 54 : NULLIFY (matrix_p)
3392 54 : IF (use_rho) THEN
3393 54 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
3394 : ELSE
3395 0 : matrix_p => scf_env%p_mix_new
3396 : END IF
3397 :
3398 : ! set up basis set lists
3399 54 : nkind = SIZE(atomic_kind_set)
3400 264 : ALLOCATE (basis_set_list(nkind))
3401 54 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
3402 :
3403 162 : ALLOCATE (dE(tb%mol%nat))
3404 :
3405 54 : nel = msao(tb%calc%bas%maxl)**2
3406 162 : ALLOCATE (t_ov(nel))
3407 162 : ALLOCATE (t_d_ov(3, nel))
3408 108 : ALLOCATE (t_dip(dip_n, nel))
3409 216 : ALLOCATE (t_i_dip(3, dip_n, nel), t_j_dip(3, dip_n, nel))
3410 162 : ALLOCATE (t_quad(quad_n, nel))
3411 216 : ALLOCATE (t_i_quad(3, quad_n, nel), t_j_quad(3, quad_n, nel))
3412 :
3413 54 : ALLOCATE (idip(dip_n), jdip(dip_n))
3414 54 : ALLOCATE (iquad(quad_n), jquad(quad_n))
3415 :
3416 300 : dE = 0.0_dp
3417 1038 : tb%grad = 0.0_dp
3418 54 : hsigma = 0.0_dp
3419 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
3420 54 : NULLIFY (nl_iterator)
3421 54 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
3422 12358 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3423 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3424 12304 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
3425 :
3426 12304 : icol = MAX(iatom, jatom)
3427 12304 : jrow = MIN(iatom, jatom)
3428 :
3429 12304 : IF (icol == jatom) THEN
3430 28328 : rij = -rij
3431 7082 : i = ikind
3432 7082 : ikind = jkind
3433 7082 : jkind = i
3434 : END IF
3435 :
3436 12304 : ityp = tb%mol%id(icol)
3437 12304 : jtyp = tb%mol%id(jrow)
3438 :
3439 49216 : r2 = DOT_PRODUCT(rij, rij)
3440 12304 : dr = SQRT(r2)
3441 12304 : IF (icol == jrow .AND. dr < same_atom) CYCLE
3442 : IF (mp_test_mode == 140 .AND. icol == jrow) CYCLE
3443 12181 : rr = SQRT(dr/(h0%rad(ikind) + h0%rad(jkind)))
3444 :
3445 : !get basis information
3446 12181 : basis_set_a => basis_set_list(ikind)%gto_basis_set
3447 12181 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
3448 12181 : first_sgfa => basis_set_a%first_sgf
3449 12181 : nsgfa => basis_set_a%nsgf_set
3450 12181 : nseti = basis_set_a%nset
3451 12181 : basis_set_b => basis_set_list(jkind)%gto_basis_set
3452 12181 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
3453 12181 : first_sgfb => basis_set_b%first_sgf
3454 12181 : nsgfb => basis_set_b%nsgf_set
3455 12181 : nsetj = basis_set_b%nset
3456 :
3457 12181 : IF (nimg == 1) THEN
3458 : ic = 1
3459 : ELSE
3460 1798 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
3461 1798 : CPASSERT(ic > 0)
3462 : END IF
3463 12181 : image_pair = (cellind(1) /= 0 .OR. cellind(2) /= 0 .OR. cellind(3) /= 0)
3464 : pot_pair_scale = MERGE(pot_image_deriv_scale, pot_deriv_scale, image_pair)
3465 : w_pair_scale = MERGE(w_image_deriv_scale, w_deriv_scale, image_pair)
3466 : h0_pair_scale = MERGE(h0_image_deriv_scale, h0_deriv_scale, image_pair)
3467 : mp_pair_scale = MERGE(mp_image_deriv_scale, mp_deriv_scale, image_pair)
3468 : IF (icol == jrow .AND. image_pair) THEN
3469 : pot_pair_scale = pot_self_image_deriv_scale
3470 : w_pair_scale = w_self_image_deriv_scale
3471 : h0_pair_scale = h0_self_image_deriv_scale
3472 : mp_pair_scale = mp_self_image_deriv_scale
3473 : END IF
3474 :
3475 12181 : NULLIFY (pblock)
3476 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
3477 12181 : row=jrow, col=icol, BLOCK=pblock, found=found)
3478 12181 : IF (.NOT. found) CPABORT("pblock not found")
3479 :
3480 12181 : NULLIFY (wblock)
3481 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
3482 12181 : row=jrow, col=icol, block=wblock, found=found)
3483 12181 : CPASSERT(found)
3484 :
3485 12181 : i_a_shift = tb%pot%vat(icol, 1)
3486 12181 : j_a_shift = tb%pot%vat(jrow, 1)
3487 48724 : idip(:) = tb%pot%vdp(:, icol, 1)
3488 48724 : jdip(:) = tb%pot%vdp(:, jrow, 1)
3489 85267 : iquad(:) = tb%pot%vqp(:, icol, 1)
3490 85267 : jquad(:) = tb%pot%vqp(:, jrow, 1)
3491 : dip_deriv_fac_pair = dip_deriv_fac
3492 : quad_deriv_fac_pair = quad_deriv_fac
3493 : overlap_deriv_fac = 1.0_dp
3494 : IF (icol == jrow) THEN
3495 : IF (mp_test_mode == 141) THEN
3496 : dip_deriv_fac_pair = 0.0_dp
3497 : quad_deriv_fac_pair = 0.0_dp
3498 : ELSEIF (mp_test_mode == 142) THEN
3499 : overlap_deriv_fac = 0.0_dp
3500 : END IF
3501 : END IF
3502 :
3503 12181 : ni = tb%calc%bas%ish_at(icol)
3504 48041 : DO iset = 1, nseti
3505 35806 : sgfi = first_sgfa(1, iset)
3506 35806 : ishift = i_a_shift + tb%pot%vsh(ni + iset, 1)
3507 35806 : nj = tb%calc%bas%ish_at(jrow)
3508 154110 : DO jset = 1, nsetj
3509 106000 : sgfj = first_sgfb(1, jset)
3510 106000 : jshift = j_a_shift + tb%pot%vsh(nj + jset, 1)
3511 :
3512 : !get integrals and derivatives
3513 : IF (mp_test_mode == 20 .OR. mp_test_mode == 21) THEN
3514 : CALL multipole_grad_cgto(tb%calc%bas%cgto(jset, jtyp), tb%calc%bas%cgto(iset, ityp), &
3515 : & r2, rij, tb%calc%bas%intcut, t_ov, t_dip, t_quad, t_d_ov, t_j_dip, t_j_quad, &
3516 : & t_i_dip, t_i_quad)
3517 : ELSEIF (mp_test_mode == 34 .OR. mp_test_mode == 35) THEN
3518 : CALL multipole_grad_cgto(tb%calc%bas%cgto(jset, jtyp), tb%calc%bas%cgto(iset, ityp), &
3519 : & r2, -rij, tb%calc%bas%intcut, t_ov, t_dip, t_quad, t_d_ov, t_j_dip, t_j_quad, &
3520 : & t_i_dip, t_i_quad)
3521 : ELSEIF (mp_test_mode == 36 .OR. mp_test_mode == 37) THEN
3522 : CALL multipole_grad_cgto(tb%calc%bas%cgto(iset, ityp), tb%calc%bas%cgto(jset, jtyp), &
3523 : & r2, -rij, tb%calc%bas%intcut, t_ov, t_dip, t_quad, t_d_ov, t_i_dip, t_i_quad, &
3524 : & t_j_dip, t_j_quad)
3525 : ELSE
3526 : CALL multipole_grad_cgto(tb%calc%bas%cgto(iset, ityp), tb%calc%bas%cgto(jset, jtyp), &
3527 : & r2, rij, tb%calc%bas%intcut, t_ov, t_dip, t_quad, t_d_ov, t_i_dip, t_i_quad, &
3528 106000 : & t_j_dip, t_j_quad)
3529 : END IF
3530 :
3531 : shpoly = (1.0_dp + h0%shpoly(iset, ikind)*rr) &
3532 106000 : *(1.0_dp + h0%shpoly(jset, jkind)*rr)
3533 : dshpoly = ((1.0_dp + h0%shpoly(iset, ikind)*rr)*h0%shpoly(jset, jkind)*rr &
3534 : + (1.0_dp + h0%shpoly(jset, jkind)*rr)*h0%shpoly(iset, ikind)*rr) &
3535 106000 : *0.5_dp/r2
3536 106000 : scal = h0%hscale(iset, jset, ikind, jkind)
3537 106000 : hij = 0.5_dp*(tb%selfenergy(ni + iset) + tb%selfenergy(nj + jset))*scal
3538 :
3539 106000 : idHdc = tb%dsedcn(ni + iset)*shpoly*scal
3540 106000 : jdHdc = tb%dsedcn(nj + jset)*shpoly*scal
3541 :
3542 106000 : itemp = 0.0_dp
3543 106000 : jtemp = 0.0_dp
3544 106000 : dgrad = 0.0_dp
3545 420200 : DO inda = 1, nsgfa(iset)
3546 314200 : ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
3547 1357480 : DO indb = 1, nsgfb(jset)
3548 937280 : ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
3549 :
3550 : IF (mp_test_mode == 20 .OR. mp_test_mode == 21 .OR. &
3551 : mp_test_mode == 34 .OR. mp_test_mode == 35) THEN
3552 : ij = indb + nsgfb(jset)*(inda - 1)
3553 : ELSE
3554 937280 : ij = inda + nsgfa(iset)*(indb - 1)
3555 : END IF
3556 :
3557 937280 : itemp = itemp + overlap_deriv_fac*idHdc*pblock(ib, ia)*t_ov(ij)
3558 937280 : jtemp = jtemp + overlap_deriv_fac*jdHdc*pblock(ib, ia)*t_ov(ij)
3559 :
3560 937280 : hp = 2*hij*pblock(ib, ia)
3561 : dgrad(:) = dgrad(:) &
3562 : + overlap_deriv_fac*(-pot_pair_scale*(ishift + jshift)*pblock(ib, ia)*t_d_ov(:, ij) &
3563 : - 2*w_pair_scale*wblock(ib, ia)*t_d_ov(:, ij) &
3564 : + h0_pair_scale*(hp*shpoly*t_d_ov(:, ij) &
3565 : + hp*dshpoly*t_ov(ij)*rij)) &
3566 : + mp_pair_scale*pblock(ib, ia)*( &
3567 937280 : dip_deriv_fac_pair*(MATMUL(t_i_dip(:, :, ij), idip) &
3568 937280 : + MATMUL(t_j_dip(:, :, ij), jdip)) &
3569 937280 : + quad_deriv_fac_pair*(MATMUL(t_i_quad(:, :, ij), iquad) &
3570 85606680 : + MATMUL(t_j_quad(:, :, ij), jquad)))
3571 :
3572 : END DO
3573 : END DO
3574 : ! Periodic self-image pairs share one atom in the sorted CP2K block representation.
3575 106000 : IF (icol == jrow) THEN
3576 27816 : IF (image_pair) THEN
3577 27816 : dE(icol) = dE(icol) + cn_image_scale*0.5_dp*(itemp + jtemp)
3578 : ELSE
3579 0 : dE(icol) = dE(icol) + 0.5_dp*(itemp + jtemp)
3580 : END IF
3581 : ELSE
3582 78184 : dE(icol) = dE(icol) + itemp
3583 78184 : dE(jrow) = dE(jrow) + jtemp
3584 : END IF
3585 424000 : tb%grad(:, icol) = tb%grad(:, icol) - dgrad
3586 424000 : tb%grad(:, jrow) = tb%grad(:, jrow) + dgrad
3587 141806 : IF (tb%use_virial) THEN
3588 98437 : IF (icol == jrow) THEN
3589 100032 : DO ia = 1, 3
3590 325104 : DO ib = 1, 3
3591 300096 : IF (image_pair .AND. nimg > 1) THEN
3592 31104 : hsigma(ia, ib) = hsigma(ia, ib) - 0.25_dp*(rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
3593 : ELSE
3594 193968 : hsigma(ia, ib) = hsigma(ia, ib) + 0.25_dp*(rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
3595 : END IF
3596 : END DO
3597 : END DO
3598 : ELSE
3599 293716 : DO ia = 1, 3
3600 954577 : DO ib = 1, 3
3601 881148 : hsigma(ia, ib) = hsigma(ia, ib) + 0.50_dp*(rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
3602 : END DO
3603 : END DO
3604 : END IF
3605 : END IF
3606 : END DO
3607 : END DO
3608 : END DO
3609 54 : CALL neighbor_list_iterator_release(nl_iterator)
3610 :
3611 54 : CALL para_env%sum(hsigma)
3612 54 : CALL para_env%sum(dE)
3613 54 : CALL para_env%sum(tb%grad)
3614 :
3615 54 : CALL tb_add_grad(tb%grad, tb%dcndr, dE, tb%mol%nat)
3616 54 : CALL tb_grad2force(qs_env, tb, para_env, 4)
3617 :
3618 54 : CALL tb_dump_sigma_component("before_h0_off", tb%sigma, para_env)
3619 702 : tb%sigma = tb%sigma + hsigma
3620 : CALL tb_dump_sigma_component("after_h0_off_direct", tb%sigma, para_env)
3621 54 : IF (tb%use_virial) THEN
3622 26 : CALL tb_add_sig(tb%sigma, tb%dcndL, dE, tb%mol%nat)
3623 : CALL tb_dump_sigma_component("after_h0_off_cn", tb%sigma, para_env)
3624 : END IF
3625 :
3626 54 : DEALLOCATE (dE)
3627 54 : DEALLOCATE (basis_set_list)
3628 54 : DEALLOCATE (t_ov, t_d_ov)
3629 54 : DEALLOCATE (t_dip, t_i_dip, t_j_dip)
3630 54 : DEALLOCATE (t_quad, t_i_quad, t_j_quad)
3631 54 : DEALLOCATE (idip, jdip, iquad, jquad)
3632 :
3633 54 : IF (tb%use_virial) CALL tb_add_stress(qs_env, tb, para_env)
3634 :
3635 : #else
3636 : MARK_USED(qs_env)
3637 : MARK_USED(use_rho)
3638 : MARK_USED(nimg)
3639 : CPABORT("Built without TBLITE")
3640 : #endif
3641 :
3642 54 : END SUBROUTINE tb_derive_dH_off
3643 :
3644 : ! **************************************************************************************************
3645 : !> \brief Run native tblite CLI and compare against CP2K/tblite.
3646 : !> \param qs_env ...
3647 : ! **************************************************************************************************
3648 54 : SUBROUTINE tb_reference_cli_compare(qs_env)
3649 :
3650 : TYPE(qs_environment_type), POINTER :: qs_env
3651 :
3652 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tb_reference_cli_compare'
3653 :
3654 : CHARACTER(LEN=32) :: acc_str, charge_str, iter_str, spin_str
3655 : CHARACTER(LEN=8) :: method
3656 : CHARACTER(LEN=8*default_path_length) :: command
3657 : CHARACTER(LEN=default_path_length) :: file_base, gen_file, grad_file, &
3658 : json_file, log_file
3659 : INTEGER :: cmdstat, exitstat, handle, iounit, &
3660 : natom, nkp, spin
3661 : LOGICAL :: do_kpoints, have_energy, have_gradient, &
3662 : have_virial, too_large, &
3663 : unsupported_kpoints
3664 : REAL(KIND=dp) :: cli_energy, cp_energy, ediff, fmax, &
3665 : fsum, vmax, vsum
3666 54 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cli_gradient, cli_virial, cp_gradient
3667 54 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
3668 54 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3669 : TYPE(cp_logger_type), POINTER :: logger
3670 : TYPE(dft_control_type), POINTER :: dft_control
3671 : TYPE(kpoint_type), POINTER :: kpoints
3672 : TYPE(mp_para_env_type), POINTER :: para_env
3673 54 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3674 : TYPE(qs_energy_type), POINTER :: energy
3675 54 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3676 : TYPE(virial_type), POINTER :: virial
3677 : TYPE(xtb_reference_cli_type) :: ref
3678 :
3679 54 : CALL timeset(routineN, handle)
3680 :
3681 54 : NULLIFY (atomic_kind_set, dft_control, energy, force, kpoints, logger, para_env, particle_set, &
3682 54 : virial, xkp)
3683 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
3684 : dft_control=dft_control, energy=energy, force=force, &
3685 : para_env=para_env, particle_set=particle_set, virial=virial, &
3686 54 : do_kpoints=do_kpoints, kpoints=kpoints)
3687 :
3688 54 : ref = dft_control%qs_control%xtb_control%reference_cli
3689 54 : IF (.NOT. ref%enabled) THEN
3690 52 : CALL timestop(handle)
3691 52 : RETURN
3692 : END IF
3693 2 : IF (.NOT. para_env%is_source()) THEN
3694 1 : CALL timestop(handle)
3695 1 : RETURN
3696 : END IF
3697 :
3698 1 : logger => cp_get_default_logger()
3699 1 : iounit = cp_logger_get_default_io_unit(logger)
3700 1 : unsupported_kpoints = .FALSE.
3701 1 : IF (do_kpoints .AND. ASSOCIATED(kpoints)) THEN
3702 0 : nkp = 0
3703 0 : CALL get_kpoint_info(kpoint=kpoints, nkp=nkp, xkp=xkp)
3704 0 : unsupported_kpoints = nkp > 1
3705 0 : IF (nkp == 1 .AND. ASSOCIATED(xkp)) unsupported_kpoints = ANY(ABS(xkp(:, 1)) > 1.0E-12_dp)
3706 : END IF
3707 0 : IF (unsupported_kpoints) THEN
3708 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
3709 0 : "tblite reference CLI check skipped: CP2K KPOINTS are active."
3710 : WRITE (UNIT=iounit, FMT="(T2,A)") &
3711 0 : "The native tblite CLI reference path does not reproduce CP2K multi-k-point sampling."
3712 0 : IF (ref%stop_on_error) CPABORT("tblite reference CLI cannot check CP2K k-point calculations")
3713 0 : CALL timestop(handle)
3714 0 : RETURN
3715 : END IF
3716 1 : natom = SIZE(particle_set)
3717 1 : method = tb_reference_method_name(dft_control%qs_control%xtb_control%tblite_method)
3718 1 : file_base = tb_join_path(ref%work_directory, ref%prefix)
3719 1 : gen_file = TRIM(file_base)//".gen"
3720 1 : grad_file = TRIM(file_base)//".grad"
3721 1 : json_file = TRIM(file_base)//".json"
3722 1 : log_file = TRIM(file_base)//".log"
3723 :
3724 1 : WRITE (charge_str, "(I0)") dft_control%charge
3725 1 : spin = MAX(0, dft_control%multiplicity - 1)
3726 1 : WRITE (spin_str, "(I0)") spin
3727 1 : WRITE (acc_str, "(ES16.8)") ref%accuracy
3728 1 : WRITE (iter_str, "(I0)") ref%iterations
3729 :
3730 1 : CALL tb_write_reference_gen(qs_env, TRIM(gen_file))
3731 :
3732 : command = TRIM(tb_shell_quote(ref%program_name))//" run --method "//TRIM(method)// &
3733 : " --charge "//TRIM(ADJUSTL(charge_str))// &
3734 : " --spin "//TRIM(ADJUSTL(spin_str))// &
3735 : " --acc "//TRIM(ADJUSTL(acc_str))// &
3736 : " --iterations "//TRIM(ADJUSTL(iter_str))// &
3737 : " --no-restart --input gen --grad "//TRIM(tb_shell_quote(grad_file))// &
3738 : " --json "//TRIM(tb_shell_quote(json_file))//" "//TRIM(tb_shell_quote(gen_file))// &
3739 1 : " > "//TRIM(tb_shell_quote(log_file))//" 2>&1"
3740 :
3741 1 : cmdstat = 0
3742 1 : exitstat = 0
3743 1 : CALL execute_command_line(TRIM(command), exitstat=exitstat, cmdstat=cmdstat)
3744 1 : IF (cmdstat /= 0 .OR. exitstat /= 0) THEN
3745 1 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "tblite reference CLI check failed to run."
3746 1 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Command: ", TRIM(command)
3747 1 : WRITE (UNIT=iounit, FMT="(T2,A,I0,T32,A,I0)") "cmdstat:", cmdstat, "exitstat:", exitstat
3748 1 : IF (ref%stop_on_error) CPABORT("tblite reference CLI command failed")
3749 1 : CALL tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file)
3750 1 : CALL timestop(handle)
3751 1 : RETURN
3752 : END IF
3753 :
3754 0 : ALLOCATE (cli_gradient(3, natom), cli_virial(3, 3))
3755 : CALL tb_read_reference_grad(TRIM(grad_file), natom, cli_energy, cli_gradient, cli_virial, &
3756 0 : have_energy, have_gradient, have_virial)
3757 :
3758 0 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "tblite reference CLI check"
3759 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Executable: ", TRIM(ref%program_name)
3760 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Method: ", TRIM(method)
3761 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Log file: ", TRIM(log_file)
3762 :
3763 0 : too_large = .FALSE.
3764 0 : IF (ref%check_energy) THEN
3765 0 : IF (have_energy) THEN
3766 0 : cp_energy = energy%total
3767 0 : ediff = ABS(cp_energy - cli_energy)
3768 : WRITE (UNIT=iounit, FMT="(T2,A,3ES22.12)") &
3769 0 : "Energy CP2K/CLI/absdiff:", cp_energy, cli_energy, ediff
3770 0 : too_large = too_large .OR. ediff > ref%error_limit
3771 : ELSE
3772 0 : WRITE (UNIT=iounit, FMT="(T2,A)") "Energy check skipped: no CLI energy found."
3773 : END IF
3774 : END IF
3775 :
3776 0 : IF (ref%check_forces) THEN
3777 0 : IF (have_gradient .AND. ASSOCIATED(force)) THEN
3778 0 : ALLOCATE (cp_gradient(3, natom))
3779 0 : CALL total_qs_force(cp_gradient, force, atomic_kind_set)
3780 0 : fsum = SUM(ABS(cp_gradient - cli_gradient))
3781 0 : fmax = MAXVAL(ABS(cp_gradient - cli_gradient))
3782 0 : WRITE (UNIT=iounit, FMT="(T2,A,2ES22.12)") "Gradient diff sum/max:", fsum, fmax
3783 0 : too_large = too_large .OR. fmax > ref%error_limit
3784 0 : DEALLOCATE (cp_gradient)
3785 : ELSE
3786 0 : WRITE (UNIT=iounit, FMT="(T2,A)") "Gradient check skipped: no CLI gradient or CP2K force found."
3787 : END IF
3788 : END IF
3789 :
3790 0 : IF (ref%check_virial) THEN
3791 0 : IF (have_virial .AND. ASSOCIATED(virial)) THEN
3792 : ! Native tblite prints the positive cell derivative; CP2K stores the PV virial
3793 : ! with the opposite sign.
3794 0 : vsum = SUM(ABS(-virial%pv_virial - cli_virial))
3795 0 : vmax = MAXVAL(ABS(-virial%pv_virial - cli_virial))
3796 0 : WRITE (UNIT=iounit, FMT="(T2,A,2ES22.12)") "Virial diff sum/max:", vsum, vmax
3797 0 : too_large = too_large .OR. vmax > ref%error_limit
3798 : ELSE
3799 0 : WRITE (UNIT=iounit, FMT="(T2,A)") "Virial check skipped: no CLI virial or CP2K virial found."
3800 : END IF
3801 : END IF
3802 :
3803 0 : IF (too_large) THEN
3804 : WRITE (UNIT=iounit, FMT="(T2,A,ES12.4)") &
3805 0 : "tblite reference CLI deviation exceeded ERROR_LIMIT = ", ref%error_limit
3806 0 : IF (ref%stop_on_error) CPABORT("tblite reference CLI deviation exceeded ERROR_LIMIT")
3807 : END IF
3808 :
3809 0 : CALL tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file)
3810 0 : DEALLOCATE (cli_gradient, cli_virial)
3811 :
3812 0 : CALL timestop(handle)
3813 :
3814 108 : END SUBROUTINE tb_reference_cli_compare
3815 :
3816 : ! **************************************************************************************************
3817 : !> \brief Map CP2K tblite method id to native tblite CLI method name.
3818 : !> \param method_id ...
3819 : !> \return ...
3820 : ! **************************************************************************************************
3821 1 : FUNCTION tb_reference_method_name(method_id) RESULT(method)
3822 : INTEGER, INTENT(IN) :: method_id
3823 : CHARACTER(LEN=8) :: method
3824 :
3825 1 : SELECT CASE (method_id)
3826 : CASE (gfn1xtb)
3827 0 : method = "gfn1"
3828 : CASE (gfn2xtb)
3829 1 : method = "gfn2"
3830 : CASE (ipea1xtb)
3831 0 : method = "ipea1"
3832 : CASE DEFAULT
3833 1 : CPABORT("Unknown tblite reference CLI method")
3834 : END SELECT
3835 :
3836 1 : END FUNCTION tb_reference_method_name
3837 :
3838 : ! **************************************************************************************************
3839 : !> \brief Join directory and filename.
3840 : !> \param directory ...
3841 : !> \param filename ...
3842 : !> \return ...
3843 : ! **************************************************************************************************
3844 1 : FUNCTION tb_join_path(directory, filename) RESULT(path)
3845 : CHARACTER(LEN=*), INTENT(IN) :: directory, filename
3846 : CHARACTER(LEN=default_path_length) :: path
3847 :
3848 1 : IF (LEN_TRIM(directory) == 0 .OR. TRIM(directory) == ".") THEN
3849 1 : path = TRIM(filename)
3850 0 : ELSE IF (directory(LEN_TRIM(directory):LEN_TRIM(directory)) == "/") THEN
3851 0 : path = TRIM(directory)//TRIM(filename)
3852 : ELSE
3853 0 : path = TRIM(directory)//"/"//TRIM(filename)
3854 : END IF
3855 :
3856 1 : END FUNCTION tb_join_path
3857 :
3858 : ! **************************************************************************************************
3859 : !> \brief Shell-quote a filename or executable path.
3860 : !> \param text ...
3861 : !> \return ...
3862 : ! **************************************************************************************************
3863 5 : FUNCTION tb_shell_quote(text) RESULT(quoted)
3864 : CHARACTER(LEN=*), INTENT(IN) :: text
3865 : CHARACTER(LEN=4*default_path_length) :: quoted
3866 :
3867 : INTEGER :: i
3868 :
3869 5 : quoted = "'"
3870 97 : DO i = 1, LEN_TRIM(text)
3871 97 : IF (text(i:i) == "'") THEN
3872 0 : quoted = TRIM(quoted)//"'\\''"
3873 : ELSE
3874 92 : quoted = TRIM(quoted)//text(i:i)
3875 : END IF
3876 : END DO
3877 5 : quoted = TRIM(quoted)//"'"
3878 :
3879 5 : END FUNCTION tb_shell_quote
3880 :
3881 : ! **************************************************************************************************
3882 : !> \brief Write current CP2K geometry as DFTB+ gen format for native tblite.
3883 : !> \param qs_env ...
3884 : !> \param filename ...
3885 : ! **************************************************************************************************
3886 1 : SUBROUTINE tb_write_reference_gen(qs_env, filename)
3887 :
3888 : TYPE(qs_environment_type), POINTER :: qs_env
3889 : CHARACTER(LEN=*), INTENT(IN) :: filename
3890 :
3891 1 : CHARACTER(LEN=2), ALLOCATABLE, DIMENSION(:) :: symbols, unique_symbols
3892 : INTEGER :: iatom, ikind, ios, natom, nuniq, unit_nr
3893 1 : INTEGER, ALLOCATABLE, DIMENSION(:) :: species
3894 : INTEGER, DIMENSION(3) :: periodic
3895 : LOGICAL :: found
3896 : REAL(KIND=dp) :: to_angstrom
3897 : TYPE(cell_type), POINTER :: cell
3898 1 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3899 :
3900 1 : NULLIFY (cell, particle_set)
3901 1 : CALL get_qs_env(qs_env=qs_env, cell=cell, particle_set=particle_set)
3902 :
3903 1 : natom = SIZE(particle_set)
3904 1 : to_angstrom = cp_unit_from_cp2k(1.0_dp, "angstrom")
3905 5 : ALLOCATE (symbols(natom), unique_symbols(natom), species(natom))
3906 1 : nuniq = 0
3907 5 : DO iatom = 1, natom
3908 4 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=symbols(iatom))
3909 4 : found = .FALSE.
3910 9 : DO ikind = 1, nuniq
3911 9 : IF (TRIM(unique_symbols(ikind)) == TRIM(symbols(iatom))) THEN
3912 : found = .TRUE.
3913 : EXIT
3914 : END IF
3915 : END DO
3916 4 : IF (.NOT. found) THEN
3917 3 : nuniq = nuniq + 1
3918 3 : unique_symbols(nuniq) = symbols(iatom)
3919 3 : ikind = nuniq
3920 : END IF
3921 5 : species(iatom) = ikind
3922 : END DO
3923 :
3924 : OPEN (NEWUNIT=unit_nr, FILE=TRIM(filename), STATUS="REPLACE", ACTION="WRITE", &
3925 1 : FORM="FORMATTED", IOSTAT=ios)
3926 1 : IF (ios /= 0) CPABORT("Could not open tblite reference CLI geometry file")
3927 :
3928 1 : CALL get_cell(cell=cell, periodic=periodic)
3929 4 : IF (ANY(periodic == 1)) THEN
3930 0 : WRITE (UNIT=unit_nr, FMT="(I0,1X,A)") natom, "S"
3931 : ELSE
3932 1 : WRITE (UNIT=unit_nr, FMT="(I0,1X,A)") natom, "C"
3933 : END IF
3934 4 : WRITE (UNIT=unit_nr, FMT="(*(A,1X))") (TRIM(unique_symbols(ikind)), ikind=1, nuniq)
3935 5 : DO iatom = 1, natom
3936 : WRITE (UNIT=unit_nr, FMT="(I0,1X,I0,3(1X,ES24.16))") &
3937 17 : iatom, species(iatom), particle_set(iatom)%r(:)*to_angstrom
3938 : END DO
3939 4 : IF (ANY(periodic == 1)) THEN
3940 0 : WRITE (UNIT=unit_nr, FMT="(3(1X,ES24.16))") 0.0_dp, 0.0_dp, 0.0_dp
3941 0 : DO ikind = 1, 3
3942 0 : WRITE (UNIT=unit_nr, FMT="(3(1X,ES24.16))") cell%hmat(:, ikind)*to_angstrom
3943 : END DO
3944 : END IF
3945 1 : CLOSE (unit_nr)
3946 :
3947 1 : DEALLOCATE (symbols, unique_symbols, species)
3948 :
3949 1 : END SUBROUTINE tb_write_reference_gen
3950 :
3951 : ! **************************************************************************************************
3952 : !> \brief Read native tblite gradient file.
3953 : !> \param filename ...
3954 : !> \param natom ...
3955 : !> \param energy ...
3956 : !> \param gradient ...
3957 : !> \param virial ...
3958 : !> \param have_energy ...
3959 : !> \param have_gradient ...
3960 : !> \param have_virial ...
3961 : ! **************************************************************************************************
3962 0 : SUBROUTINE tb_read_reference_grad(filename, natom, energy, gradient, virial, &
3963 : have_energy, have_gradient, have_virial)
3964 :
3965 : CHARACTER(LEN=*), INTENT(IN) :: filename
3966 : INTEGER, INTENT(IN) :: natom
3967 : REAL(KIND=dp), INTENT(OUT) :: energy
3968 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: gradient, virial
3969 : LOGICAL, INTENT(OUT) :: have_energy, have_gradient, have_virial
3970 :
3971 : CHARACTER(LEN=1024) :: line
3972 : INTEGER :: ios, nread, unit_nr
3973 : LOGICAL :: exists
3974 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: values
3975 :
3976 0 : have_energy = .FALSE.
3977 0 : have_gradient = .FALSE.
3978 0 : have_virial = .FALSE.
3979 0 : energy = 0.0_dp
3980 0 : gradient = 0.0_dp
3981 0 : virial = 0.0_dp
3982 :
3983 0 : INQUIRE (FILE=TRIM(filename), EXIST=exists)
3984 0 : IF (.NOT. exists) RETURN
3985 :
3986 : OPEN (NEWUNIT=unit_nr, FILE=TRIM(filename), STATUS="OLD", ACTION="READ", &
3987 0 : FORM="FORMATTED", IOSTAT=ios)
3988 0 : IF (ios /= 0) RETURN
3989 :
3990 : DO
3991 0 : READ (UNIT=unit_nr, FMT="(A)", IOSTAT=ios) line
3992 0 : IF (ios /= 0) EXIT
3993 0 : IF (INDEX(line, "energy :real:0:") > 0) THEN
3994 0 : READ (UNIT=unit_nr, FMT="(A)", IOSTAT=ios) line
3995 0 : IF (ios == 0) THEN
3996 0 : READ (line, *, IOSTAT=ios) energy
3997 0 : have_energy = ios == 0
3998 : END IF
3999 0 : ELSE IF (INDEX(line, "gradient :real:2:3,") > 0) THEN
4000 0 : ALLOCATE (values(3*natom))
4001 0 : CALL tb_read_real_values(unit_nr, values, nread)
4002 0 : IF (nread == 3*natom) THEN
4003 0 : CALL tb_values_to_matrix(values, gradient)
4004 0 : have_gradient = .TRUE.
4005 : END IF
4006 0 : DEALLOCATE (values)
4007 0 : ELSE IF (INDEX(line, "virial :real:2:3,3") > 0) THEN
4008 0 : ALLOCATE (values(9))
4009 0 : CALL tb_read_real_values(unit_nr, values, nread)
4010 0 : IF (nread == 9) THEN
4011 0 : CALL tb_values_to_matrix_rowmajor(values, virial)
4012 0 : have_virial = .TRUE.
4013 : END IF
4014 0 : DEALLOCATE (values)
4015 : END IF
4016 : END DO
4017 0 : CLOSE (unit_nr)
4018 :
4019 0 : END SUBROUTINE tb_read_reference_grad
4020 :
4021 : ! **************************************************************************************************
4022 : !> \brief Read a fixed number of real values from following lines.
4023 : !> \param unit_nr ...
4024 : !> \param values ...
4025 : !> \param nread ...
4026 : ! **************************************************************************************************
4027 0 : SUBROUTINE tb_read_real_values(unit_nr, values, nread)
4028 :
4029 : INTEGER, INTENT(IN) :: unit_nr
4030 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: values
4031 : INTEGER, INTENT(OUT) :: nread
4032 :
4033 : CHARACTER(LEN=1024) :: line
4034 : INTEGER :: ios
4035 :
4036 0 : nread = 0
4037 0 : DO WHILE (nread < SIZE(values))
4038 0 : READ (UNIT=unit_nr, FMT="(A)", IOSTAT=ios) line
4039 0 : IF (ios /= 0) EXIT
4040 0 : CALL tb_parse_real_line(line, values, nread)
4041 : END DO
4042 :
4043 0 : END SUBROUTINE tb_read_real_values
4044 :
4045 : ! **************************************************************************************************
4046 : !> \brief Parse real values from one text line.
4047 : !> \param line ...
4048 : !> \param values ...
4049 : !> \param nread ...
4050 : ! **************************************************************************************************
4051 0 : SUBROUTINE tb_parse_real_line(line, values, nread)
4052 :
4053 : CHARACTER(LEN=*), INTENT(IN) :: line
4054 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: values
4055 : INTEGER, INTENT(INOUT) :: nread
4056 :
4057 : CHARACTER(LEN=128) :: token
4058 : INTEGER :: first, ios, last, pos
4059 :
4060 0 : pos = 1
4061 0 : DO WHILE (pos <= LEN_TRIM(line) .AND. nread < SIZE(values))
4062 0 : DO WHILE (pos <= LEN_TRIM(line) .AND. INDEX(" ,[]", line(pos:pos)) > 0)
4063 0 : pos = pos + 1
4064 : END DO
4065 0 : IF (pos > LEN_TRIM(line)) EXIT
4066 : first = pos
4067 0 : DO WHILE (pos <= LEN_TRIM(line) .AND. INDEX(" ,[]", line(pos:pos)) == 0)
4068 0 : pos = pos + 1
4069 : END DO
4070 0 : last = pos - 1
4071 0 : token = line(first:last)
4072 0 : READ (token, *, IOSTAT=ios) values(nread + 1)
4073 0 : IF (ios == 0) nread = nread + 1
4074 : END DO
4075 :
4076 0 : END SUBROUTINE tb_parse_real_line
4077 :
4078 : ! **************************************************************************************************
4079 : !> \brief Convert flat native tblite values to CP2K atom-major matrix layout.
4080 : !> \param values ...
4081 : !> \param matrix ...
4082 : ! **************************************************************************************************
4083 0 : SUBROUTINE tb_values_to_matrix(values, matrix)
4084 :
4085 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: values
4086 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: matrix
4087 :
4088 : INTEGER :: i, j, n
4089 :
4090 0 : n = 0
4091 0 : DO j = 1, SIZE(matrix, 2)
4092 0 : DO i = 1, SIZE(matrix, 1)
4093 0 : n = n + 1
4094 0 : matrix(i, j) = values(n)
4095 : END DO
4096 : END DO
4097 :
4098 0 : END SUBROUTINE tb_values_to_matrix
4099 :
4100 : ! **************************************************************************************************
4101 : !> \brief Convert flat native tblite row-major values to a 2D matrix.
4102 : !> \param values ...
4103 : !> \param matrix ...
4104 : ! **************************************************************************************************
4105 0 : SUBROUTINE tb_values_to_matrix_rowmajor(values, matrix)
4106 :
4107 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: values
4108 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: matrix
4109 :
4110 : INTEGER :: i, j, n
4111 :
4112 0 : n = 0
4113 0 : DO i = 1, SIZE(matrix, 1)
4114 0 : DO j = 1, SIZE(matrix, 2)
4115 0 : n = n + 1
4116 0 : matrix(i, j) = values(n)
4117 : END DO
4118 : END DO
4119 :
4120 0 : END SUBROUTINE tb_values_to_matrix_rowmajor
4121 :
4122 : ! **************************************************************************************************
4123 : !> \brief Delete temporary files unless requested otherwise.
4124 : !> \param ref ...
4125 : !> \param gen_file ...
4126 : !> \param grad_file ...
4127 : !> \param json_file ...
4128 : !> \param log_file ...
4129 : ! **************************************************************************************************
4130 1 : SUBROUTINE tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file)
4131 :
4132 : TYPE(xtb_reference_cli_type), INTENT(IN) :: ref
4133 : CHARACTER(LEN=*), INTENT(IN) :: gen_file, grad_file, json_file, log_file
4134 :
4135 1 : IF (ref%keep_files) RETURN
4136 1 : CALL tb_delete_file(gen_file)
4137 1 : CALL tb_delete_file(grad_file)
4138 1 : CALL tb_delete_file(json_file)
4139 1 : CALL tb_delete_file(log_file)
4140 :
4141 : END SUBROUTINE tb_reference_cleanup
4142 :
4143 : ! **************************************************************************************************
4144 : !> \brief Delete a file if it exists.
4145 : !> \param filename ...
4146 : ! **************************************************************************************************
4147 4 : SUBROUTINE tb_delete_file(filename)
4148 :
4149 : CHARACTER(LEN=*), INTENT(IN) :: filename
4150 :
4151 : INTEGER :: ios, unit_nr
4152 : LOGICAL :: exists
4153 :
4154 4 : INQUIRE (FILE=TRIM(filename), EXIST=exists)
4155 4 : IF (.NOT. exists) RETURN
4156 2 : OPEN (NEWUNIT=unit_nr, FILE=TRIM(filename), STATUS="OLD", IOSTAT=ios)
4157 2 : IF (ios == 0) CLOSE (unit_nr, STATUS="DELETE")
4158 :
4159 : END SUBROUTINE tb_delete_file
4160 :
4161 : ! **************************************************************************************************
4162 : !> \brief Dump cumulative tblite virial pieces for local debugging.
4163 : !> \param label ...
4164 : !> \param sigma ...
4165 : !> \param para_env ...
4166 : ! **************************************************************************************************
4167 0 : SUBROUTINE tb_dump_sigma_component(label, sigma, para_env)
4168 :
4169 : CHARACTER(LEN=*), INTENT(IN) :: label
4170 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sigma
4171 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
4172 :
4173 : CHARACTER(LEN=default_path_length) :: dump_file
4174 : INTEGER :: dump_status, dump_unit, i
4175 :
4176 0 : dump_status = 1
4177 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
4178 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_SIGMA_COMPONENT_DUMP", dump_file, STATUS=dump_status)
4179 : #endif
4180 : IF (dump_status /= 0) RETURN
4181 :
4182 : OPEN (NEWUNIT=dump_unit, FILE=TRIM(dump_file), STATUS="UNKNOWN", &
4183 : POSITION="APPEND", ACTION="WRITE")
4184 : WRITE (dump_unit, "(A)") TRIM(label)
4185 : DO i = 1, 3
4186 : WRITE (dump_unit, "(3(1X,ES24.16))") sigma(i, :)/para_env%num_pe
4187 : END DO
4188 : CLOSE (dump_unit)
4189 :
4190 : END SUBROUTINE tb_dump_sigma_component
4191 :
4192 : ! **************************************************************************************************
4193 : !> \brief save stress tensor
4194 : !> \param qs_env ...
4195 : !> \param tb ...
4196 : !> \param para_env ...
4197 : ! **************************************************************************************************
4198 26 : SUBROUTINE tb_add_stress(qs_env, tb, para_env)
4199 :
4200 : TYPE(qs_environment_type) :: qs_env
4201 : TYPE(tblite_type) :: tb
4202 : TYPE(mp_para_env_type) :: para_env
4203 :
4204 : CHARACTER(LEN=default_path_length) :: dump_file
4205 : INTEGER :: dump_status, dump_unit, i
4206 : INTEGER, DIMENSION(3) :: periodic
4207 : TYPE(cell_type), POINTER :: cell
4208 : TYPE(virial_type), POINTER :: virial
4209 :
4210 26 : NULLIFY (virial, cell)
4211 26 : CALL get_qs_env(qs_env=qs_env, virial=virial, cell=cell)
4212 26 : CALL get_cell(cell=cell, periodic=periodic)
4213 :
4214 32 : IF (ALL(periodic == 0)) THEN
4215 : CALL cp_warn(__LOCATION__, &
4216 : "tblite stress tensor requested for an isolated system. "// &
4217 : "The reported virial is useful for finite-difference checks, "// &
4218 2 : "but it is not a physically meaningful bulk stress for an isolated molecule.")
4219 : END IF
4220 :
4221 26 : dump_status = 1
4222 : #if defined(__TBLITE_DEBUG_DIAGNOSTICS)
4223 : CALL GET_ENVIRONMENT_VARIABLE("CP2K_TBLITE_VIRIAL_DUMP", dump_file, STATUS=dump_status)
4224 : #endif
4225 : IF (dump_status == 0) THEN
4226 : OPEN (NEWUNIT=dump_unit, FILE=TRIM(dump_file), STATUS="UNKNOWN", &
4227 : POSITION="APPEND", ACTION="WRITE")
4228 : WRITE (dump_unit, "(A)") "sigma"
4229 : DO i = 1, 3
4230 : WRITE (dump_unit, "(3(1X,ES24.16))") tb%sigma(i, :)/para_env%num_pe
4231 : END DO
4232 : CLOSE (dump_unit)
4233 : END IF
4234 :
4235 338 : virial%pv_virial = virial%pv_virial - tb%sigma/para_env%num_pe
4236 :
4237 26 : END SUBROUTINE tb_add_stress
4238 :
4239 : ! **************************************************************************************************
4240 : !> \brief add contrib. to gradient
4241 : !> \param grad ...
4242 : !> \param deriv ...
4243 : !> \param dE ...
4244 : !> \param natom ...
4245 : ! **************************************************************************************************
4246 108 : SUBROUTINE tb_add_grad(grad, deriv, dE, natom)
4247 :
4248 : REAL(KIND=dp), DIMENSION(:, :) :: grad
4249 : REAL(KIND=dp), DIMENSION(:, :, :) :: deriv
4250 : REAL(KIND=dp), DIMENSION(:) :: dE
4251 : INTEGER :: natom
4252 :
4253 : INTEGER :: i, j
4254 :
4255 600 : DO i = 1, natom
4256 3116 : DO j = 1, natom
4257 10556 : grad(:, i) = grad(:, i) + deriv(:, i, j)*dE(j)
4258 : END DO
4259 : END DO
4260 :
4261 108 : END SUBROUTINE tb_add_grad
4262 :
4263 : ! **************************************************************************************************
4264 : !> \brief add contrib. to sigma
4265 : !> \param sig ...
4266 : !> \param deriv ...
4267 : !> \param dE ...
4268 : !> \param natom ...
4269 : ! **************************************************************************************************
4270 52 : SUBROUTINE tb_add_sig(sig, deriv, dE, natom)
4271 :
4272 : REAL(KIND=dp), DIMENSION(:, :) :: sig
4273 : REAL(KIND=dp), DIMENSION(:, :, :) :: deriv
4274 : REAL(KIND=dp), DIMENSION(:) :: dE
4275 : INTEGER :: natom
4276 :
4277 : INTEGER :: i, j
4278 :
4279 208 : DO i = 1, 3
4280 1012 : DO j = 1, natom
4281 3372 : sig(:, i) = sig(:, i) + deriv(:, i, j)*dE(j)
4282 : END DO
4283 : END DO
4284 :
4285 52 : END SUBROUTINE tb_add_sig
4286 :
4287 1488176 : END MODULE tblite_interface
|
