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_basis_type, ONLY: get_cutoff
21 : USE tblite_container, ONLY: container_cache
22 : USE tblite_integral_multipole, ONLY: multipole_cgto, multipole_grad_cgto, maxl, msao
23 : USE tblite_scf_info, ONLY: scf_info, atom_resolved, shell_resolved, &
24 : orbital_resolved, not_used
25 : USE tblite_scf_potential, ONLY: potential_type, new_potential
26 : USE tblite_wavefunction_type, ONLY: wavefunction_type, new_wavefunction
27 : USE tblite_xtb_calculator, ONLY: xtb_calculator, new_xtb_calculator
28 : USE tblite_xtb_gfn1, ONLY: new_gfn1_calculator
29 : USE tblite_xtb_gfn2, ONLY: new_gfn2_calculator
30 : USE tblite_xtb_h0, ONLY: get_selfenergy, get_hamiltonian, get_occupation, &
31 : get_hamiltonian_gradient, tb_hamiltonian
32 : USE tblite_xtb_ipea1, ONLY: new_ipea1_calculator
33 : #endif
34 : USE ai_contraction, ONLY: block_add, &
35 : contraction
36 : USE ai_overlap, ONLY: overlap_ab
37 : USE atomic_kind_types, ONLY: atomic_kind_type, get_atomic_kind, get_atomic_kind_set
38 : USE atprop_types, ONLY: atprop_type
39 : USE basis_set_types, ONLY: gto_basis_set_type, gto_basis_set_p_type, &
40 : & allocate_gto_basis_set, write_gto_basis_set, process_gto_basis
41 : USE cell_types, ONLY: cell_type, get_cell
42 : USE cp_blacs_env, ONLY: cp_blacs_env_type
43 : USE cp_control_types, ONLY: dft_control_type
44 : USE cp_dbcsr_api, ONLY: dbcsr_type, dbcsr_p_type, dbcsr_create, dbcsr_add, &
45 : dbcsr_get_block_p, dbcsr_finalize, &
46 : dbcsr_iterator_type, dbcsr_iterator_blocks_left, &
47 : dbcsr_iterator_start, dbcsr_iterator_stop, &
48 : dbcsr_iterator_next_block
49 : USE cp_dbcsr_contrib, ONLY: dbcsr_print
50 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
51 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
52 : USE cp_log_handling, ONLY: cp_get_default_logger, &
53 : cp_logger_type, cp_logger_get_default_io_unit
54 : USE cp_output_handling, ONLY: cp_print_key_should_output, &
55 : cp_print_key_unit_nr, cp_print_key_finished_output
56 : USE input_constants, ONLY: gfn1xtb, gfn2xtb, ipea1xtb
57 : USE input_section_types, ONLY: section_vals_val_get
58 : USE kinds, ONLY: dp, default_string_length
59 : USE kpoint_types, ONLY: get_kpoint_info, kpoint_type
60 : USE memory_utilities, ONLY: reallocate
61 : USE message_passing, ONLY: mp_para_env_type
62 : USE mulliken, ONLY: ao_charges
63 : USE orbital_pointers, ONLY: ncoset
64 : USE particle_types, ONLY: particle_type
65 : USE qs_charge_mixing, ONLY: charge_mixing
66 : USE qs_condnum, ONLY: overlap_condnum
67 : USE qs_energy_types, ONLY: qs_energy_type
68 : USE qs_environment_types, ONLY: get_qs_env, qs_environment_type
69 : USE qs_force_types, ONLY: qs_force_type
70 : USE qs_integral_utils, ONLY: basis_set_list_setup, get_memory_usage
71 : USE qs_kind_types, ONLY: get_qs_kind, qs_kind_type, get_qs_kind_set
72 : USE qs_ks_types, ONLY: qs_ks_env_type, set_ks_env
73 : USE qs_neighbor_list_types, ONLY: neighbor_list_iterator_create, neighbor_list_iterate, &
74 : get_iterator_info, neighbor_list_set_p_type, &
75 : neighbor_list_iterator_p_type, neighbor_list_iterator_release
76 : USE qs_overlap, ONLY: create_sab_matrix
77 : USE qs_rho_types, ONLY: qs_rho_get, qs_rho_type
78 : USE qs_scf_types, ONLY: qs_scf_env_type
79 : USE input_section_types, ONLY: section_vals_get_subs_vals, section_vals_type
80 : USE string_utilities, ONLY: integer_to_string
81 : USE tblite_types, ONLY: tblite_type, deallocate_tblite_type, allocate_tblite_type
82 : USE virial_types, ONLY: virial_type
83 : USE xtb_types, ONLY: get_xtb_atom_param, xtb_atom_type
84 : USE xtb_types, ONLY: xtb_atom_type
85 :
86 : #include "./base/base_uses.f90"
87 : IMPLICIT NONE
88 :
89 : PRIVATE
90 :
91 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tblite_interface'
92 :
93 : INTEGER, PARAMETER :: dip_n = 3
94 : INTEGER, PARAMETER :: quad_n = 6
95 : REAL(KIND=dp), PARAMETER :: same_atom = 0.00001_dp
96 :
97 : PUBLIC :: tb_set_calculator, tb_init_geometry, tb_init_wf
98 : PUBLIC :: tb_get_basis, build_tblite_matrices
99 : PUBLIC :: tb_get_energy, tb_update_charges, tb_ham_add_coulomb
100 : PUBLIC :: tb_get_multipole
101 : PUBLIC :: tb_derive_dH_diag, tb_derive_dH_off
102 :
103 : CONTAINS
104 :
105 : ! **************************************************************************************************
106 : !> \brief intialize geometry objects ...
107 : !> \param qs_env ...
108 : !> \param tb ...
109 : ! **************************************************************************************************
110 48 : SUBROUTINE tb_init_geometry(qs_env, tb)
111 :
112 : TYPE(qs_environment_type), POINTER :: qs_env
113 : TYPE(tblite_type), POINTER :: tb
114 :
115 : #if defined(__TBLITE)
116 :
117 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tblite_init_geometry'
118 :
119 : TYPE(cell_type), POINTER :: cell
120 48 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
121 48 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
122 : INTEGER :: iatom, natom
123 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xyz
124 : INTEGER :: handle, ikind
125 : INTEGER, DIMENSION(3) :: periodic
126 : LOGICAL, DIMENSION(3) :: lperiod
127 :
128 48 : CALL timeset(routineN, handle)
129 :
130 : !get info from environment vaiarable
131 48 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell, qs_kind_set=qs_kind_set)
132 :
133 : !get information about particles
134 48 : natom = SIZE(particle_set)
135 144 : ALLOCATE (xyz(3, natom))
136 48 : CALL allocate_tblite_type(tb)
137 144 : ALLOCATE (tb%el_num(natom))
138 264 : tb%el_num = -9
139 264 : DO iatom = 1, natom
140 864 : xyz(:, iatom) = particle_set(iatom)%r(:)
141 216 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
142 216 : CALL get_qs_kind(qs_kind_set(ikind), zatom=tb%el_num(iatom))
143 480 : IF (tb%el_num(iatom) < 1 .OR. tb%el_num(iatom) > 85) THEN
144 0 : CPABORT("only elements 1-85 are supported by tblite")
145 : END IF
146 : END DO
147 :
148 : !get information about cell / lattice
149 48 : CALL get_cell(cell=cell, periodic=periodic)
150 48 : lperiod(1) = periodic(1) == 1
151 48 : lperiod(2) = periodic(2) == 1
152 48 : lperiod(3) = periodic(3) == 1
153 :
154 : !prepare for the call to the dispersion function
155 48 : CALL new(tb%mol, tb%el_num, xyz, lattice=cell%hmat, periodic=lperiod)
156 :
157 48 : DEALLOCATE (xyz)
158 :
159 48 : CALL timestop(handle)
160 :
161 : #else
162 : MARK_USED(qs_env)
163 : MARK_USED(tb)
164 : CPABORT("Built without TBLITE")
165 : #endif
166 :
167 48 : END SUBROUTINE tb_init_geometry
168 :
169 : ! **************************************************************************************************
170 : !> \brief updating coordinates...
171 : !> \param qs_env ...
172 : !> \param tb ...
173 : ! **************************************************************************************************
174 928 : SUBROUTINE tb_update_geometry(qs_env, tb)
175 :
176 : TYPE(qs_environment_type) :: qs_env
177 : TYPE(tblite_type) :: tb
178 :
179 : #if defined(__TBLITE)
180 :
181 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tblite_update_geometry'
182 :
183 928 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
184 : INTEGER :: iatom, natom
185 928 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xyz
186 : INTEGER :: handle
187 :
188 928 : CALL timeset(routineN, handle)
189 :
190 : !get info from environment vaiarable
191 928 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
192 :
193 : !get information about particles
194 928 : natom = SIZE(particle_set)
195 2784 : ALLOCATE (xyz(3, natom))
196 5162 : DO iatom = 1, natom
197 17864 : xyz(:, iatom) = particle_set(iatom)%r(:)
198 : END DO
199 17864 : tb%mol%xyz(:, :) = xyz
200 :
201 928 : DEALLOCATE (xyz)
202 :
203 928 : CALL timestop(handle)
204 :
205 : #else
206 : MARK_USED(qs_env)
207 : MARK_USED(tb)
208 : CPABORT("Built without TBLITE")
209 : #endif
210 :
211 928 : END SUBROUTINE tb_update_geometry
212 :
213 : ! **************************************************************************************************
214 : !> \brief initialize wavefunction ...
215 : !> \param tb ...
216 : ! **************************************************************************************************
217 48 : SUBROUTINE tb_init_wf(tb)
218 :
219 : TYPE(tblite_type), POINTER :: tb
220 :
221 : #if defined(__TBLITE)
222 :
223 : INTEGER, PARAMETER :: nSpin = 1 !number of spin channels
224 :
225 : TYPE(scf_info) :: info
226 :
227 48 : info = tb%calc%variable_info()
228 0 : IF (info%charge > shell_resolved) CPABORT("tblite: no support for orbital resolved charge")
229 48 : IF (info%dipole > atom_resolved) CPABORT("tblite: no support for shell resolved dipole moment")
230 48 : IF (info%quadrupole > atom_resolved) &
231 0 : CPABORT("tblite: no support shell resolved quadrupole moment")
232 :
233 48 : CALL new_wavefunction(tb%wfn, tb%mol%nat, tb%calc%bas%nsh, tb%calc%bas%nao, nSpin, 0.0_dp)
234 :
235 48 : CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
236 :
237 : !allocate quantities later required
238 336 : ALLOCATE (tb%e_hal(tb%mol%nat), tb%e_rep(tb%mol%nat), tb%e_disp(tb%mol%nat))
239 240 : ALLOCATE (tb%e_scd(tb%mol%nat), tb%e_es(tb%mol%nat))
240 144 : ALLOCATE (tb%selfenergy(tb%calc%bas%nsh))
241 144 : IF (ALLOCATED(tb%calc%ncoord)) ALLOCATE (tb%cn(tb%mol%nat))
242 :
243 : #else
244 : MARK_USED(tb)
245 : CPABORT("Built without TBLITE")
246 : #endif
247 :
248 48 : END SUBROUTINE tb_init_wf
249 :
250 : ! **************************************************************************************************
251 : !> \brief ...
252 : !> \param tb ...
253 : !> \param typ ...
254 : ! **************************************************************************************************
255 48 : SUBROUTINE tb_set_calculator(tb, typ)
256 :
257 : TYPE(tblite_type), POINTER :: tb
258 : INTEGER :: typ
259 :
260 : #if defined(__TBLITE)
261 :
262 48 : TYPE(error_type), ALLOCATABLE :: error
263 :
264 48 : SELECT CASE (typ)
265 : CASE default
266 0 : CPABORT("Unknown xtb type")
267 : CASE (gfn1xtb)
268 26 : CALL new_gfn1_calculator(tb%calc, tb%mol, error)
269 : CASE (gfn2xtb)
270 12 : CALL new_gfn2_calculator(tb%calc, tb%mol, error)
271 : CASE (ipea1xtb)
272 48 : CALL new_ipea1_calculator(tb%calc, tb%mol, error)
273 : END SELECT
274 :
275 : #else
276 : MARK_USED(tb)
277 : MARK_USED(typ)
278 : CPABORT("Built without TBLITE")
279 : #endif
280 :
281 48 : END SUBROUTINE tb_set_calculator
282 :
283 : ! **************************************************************************************************
284 : !> \brief ...
285 : !> \param qs_env ...
286 : !> \param tb ...
287 : !> \param para_env ...
288 : ! **************************************************************************************************
289 928 : SUBROUTINE tb_init_ham(qs_env, tb, para_env)
290 :
291 : TYPE(qs_environment_type) :: qs_env
292 : TYPE(tblite_type) :: tb
293 : TYPE(mp_para_env_type) :: para_env
294 :
295 : #if defined(__TBLITE)
296 :
297 928 : TYPE(container_cache) :: hcache, rcache
298 :
299 5162 : tb%e_hal = 0.0_dp
300 5162 : tb%e_rep = 0.0_dp
301 5162 : tb%e_disp = 0.0_dp
302 928 : IF (ALLOCATED(tb%grad)) THEN
303 16952 : tb%grad = 0.0_dp
304 880 : CALL tb_zero_force(qs_env)
305 : END IF
306 12064 : tb%sigma = 0.0_dp
307 :
308 928 : IF (ALLOCATED(tb%calc%halogen)) THEN
309 694 : CALL tb%calc%halogen%update(tb%mol, hcache)
310 694 : IF (ALLOCATED(tb%grad)) THEN
311 13498 : tb%grad = 0.0_dp
312 : CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal, &
313 658 : & tb%grad, tb%sigma)
314 658 : CALL tb_grad2force(qs_env, tb, para_env, 0)
315 : ELSE
316 36 : CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal)
317 : END IF
318 : END IF
319 :
320 928 : IF (ALLOCATED(tb%calc%repulsion)) THEN
321 928 : CALL tb%calc%repulsion%update(tb%mol, rcache)
322 928 : IF (ALLOCATED(tb%grad)) THEN
323 16952 : tb%grad = 0.0_dp
324 : CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep, &
325 880 : & tb%grad, tb%sigma)
326 880 : CALL tb_grad2force(qs_env, tb, para_env, 1)
327 : ELSE
328 48 : CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep)
329 : END IF
330 : END IF
331 :
332 928 : IF (ALLOCATED(tb%calc%dispersion)) THEN
333 928 : CALL tb%calc%dispersion%update(tb%mol, tb%dcache)
334 928 : IF (ALLOCATED(tb%grad)) THEN
335 16952 : tb%grad = 0.0_dp
336 : CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp, &
337 880 : & tb%grad, tb%sigma)
338 880 : CALL tb_grad2force(qs_env, tb, para_env, 2)
339 : ELSE
340 48 : CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp)
341 : END IF
342 : END IF
343 :
344 928 : CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
345 928 : IF (ALLOCATED(tb%calc%coulomb)) THEN
346 928 : CALL tb%calc%coulomb%update(tb%mol, tb%cache)
347 : END IF
348 :
349 928 : IF (ALLOCATED(tb%grad)) THEN
350 880 : IF (ALLOCATED(tb%calc%ncoord)) THEN
351 880 : CALL tb%calc%ncoord%get_cn(tb%mol, tb%cn, tb%dcndr, tb%dcndL)
352 : END IF
353 : CALL get_selfenergy(tb%calc%h0, tb%mol%id, tb%calc%bas%ish_at, &
354 880 : & tb%calc%bas%nsh_id, cn=tb%cn, selfenergy=tb%selfenergy, dsedcn=tb%dsedcn)
355 : ELSE
356 48 : IF (ALLOCATED(tb%calc%ncoord)) THEN
357 48 : CALL tb%calc%ncoord%get_cn(tb%mol, tb%cn)
358 : END IF
359 : CALL get_selfenergy(tb%calc%h0, tb%mol%id, tb%calc%bas%ish_at, &
360 48 : & tb%calc%bas%nsh_id, cn=tb%cn, selfenergy=tb%selfenergy, dsedcn=tb%dsedcn)
361 : END IF
362 :
363 : #else
364 : MARK_USED(qs_env)
365 : MARK_USED(tb)
366 : MARK_USED(para_env)
367 : CPABORT("Built without TBLITE")
368 : #endif
369 :
370 928 : END SUBROUTINE tb_init_ham
371 :
372 : ! **************************************************************************************************
373 : !> \brief ...
374 : !> \param qs_env ...
375 : !> \param tb ...
376 : !> \param energy ...
377 : ! **************************************************************************************************
378 11320 : SUBROUTINE tb_get_energy(qs_env, tb, energy)
379 :
380 : TYPE(qs_environment_type), POINTER :: qs_env
381 : TYPE(tblite_type), POINTER :: tb
382 : TYPE(qs_energy_type), POINTER :: energy
383 :
384 : #if defined(__TBLITE)
385 :
386 : INTEGER :: iounit
387 : TYPE(cp_logger_type), POINTER :: logger
388 : TYPE(section_vals_type), POINTER :: scf_section
389 :
390 11320 : NULLIFY (scf_section, logger)
391 :
392 11320 : logger => cp_get_default_logger()
393 11320 : iounit = cp_logger_get_default_io_unit(logger)
394 11320 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
395 :
396 60984 : energy%repulsive = SUM(tb%e_rep)
397 60984 : energy%el_stat = SUM(tb%e_es)
398 60984 : energy%dispersion = SUM(tb%e_disp)
399 60984 : energy%dispersion_sc = SUM(tb%e_scd)
400 60984 : energy%xtb_xb_inter = SUM(tb%e_hal)
401 :
402 : energy%total = energy%core + energy%repulsive + energy%el_stat + energy%dispersion &
403 : + energy%dispersion_sc + energy%xtb_xb_inter + energy%kTS &
404 11320 : + energy%efield + energy%qmmm_el
405 :
406 : iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
407 11320 : extension=".scfLog")
408 11320 : IF (iounit > 0) THEN
409 : WRITE (UNIT=iounit, FMT="(/,(T9,A,T60,F20.10))") &
410 222 : "Repulsive pair potential energy: ", energy%repulsive, &
411 222 : "Zeroth order Hamiltonian energy: ", energy%core, &
412 222 : "Electrostatic energy: ", energy%el_stat, &
413 222 : "Self-consistent dispersion energy: ", energy%dispersion_sc, &
414 444 : "Non-self consistent dispersion energy: ", energy%dispersion
415 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
416 222 : "Correction for halogen bonding: ", energy%xtb_xb_inter
417 222 : IF (ABS(energy%efield) > 1.e-12_dp) THEN
418 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
419 0 : "Electric field interaction energy: ", energy%efield
420 : END IF
421 222 : IF (qs_env%qmmm) THEN
422 : WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") &
423 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
424 : END IF
425 : END IF
426 : CALL cp_print_key_finished_output(iounit, logger, scf_section, &
427 11320 : "PRINT%DETAILED_ENERGY")
428 :
429 : #else
430 : MARK_USED(qs_env)
431 : MARK_USED(tb)
432 : MARK_USED(energy)
433 : CPABORT("Built without TBLITE")
434 : #endif
435 :
436 11320 : END SUBROUTINE tb_get_energy
437 :
438 : ! **************************************************************************************************
439 : !> \brief ...
440 : !> \param tb ...
441 : !> \param gto_basis_set ...
442 : !> \param element_symbol ...
443 : !> \param param ...
444 : !> \param occ ...
445 : ! **************************************************************************************************
446 130 : SUBROUTINE tb_get_basis(tb, gto_basis_set, element_symbol, param, occ)
447 :
448 : TYPE(tblite_type), POINTER :: tb
449 : TYPE(gto_basis_set_type), POINTER :: gto_basis_set
450 : CHARACTER(len=2), INTENT(IN) :: element_symbol
451 : TYPE(xtb_atom_type), POINTER :: param
452 : INTEGER, DIMENSION(5), INTENT(out) :: occ
453 :
454 : #if defined(__TBLITE)
455 :
456 : CHARACTER(LEN=default_string_length) :: sng
457 : INTEGER :: ang, i_type, id_atom, ind_ao, ipgf, iSH, &
458 : ishell, ityp, maxl, mprim, natorb, &
459 : nset, nshell
460 : LOGICAL :: do_ortho
461 :
462 130 : CALL allocate_gto_basis_set(gto_basis_set)
463 :
464 : !identifying element in the bas data
465 130 : CALL symbol_to_number(i_type, element_symbol)
466 280 : DO id_atom = 1, tb%mol%nat
467 280 : IF (i_type == tb%el_num(id_atom)) EXIT
468 : END DO
469 130 : param%z = i_type
470 130 : param%symbol = element_symbol
471 130 : param%defined = .TRUE.
472 130 : ityp = tb%mol%id(id_atom)
473 :
474 : !getting size information
475 130 : nset = tb%calc%bas%nsh_id(ityp)
476 130 : nshell = 1
477 130 : mprim = 0
478 420 : DO ishell = 1, nset
479 420 : mprim = MAX(mprim, tb%calc%bas%cgto(ishell, ityp)%nprim)
480 : END DO
481 130 : param%nshell = nset
482 130 : natorb = 0
483 :
484 : !write basis set information
485 130 : CALL integer_to_string(mprim, sng)
486 130 : gto_basis_set%name = element_symbol//"_STO-"//TRIM(sng)//"G"
487 130 : gto_basis_set%nset = nset
488 130 : CALL reallocate(gto_basis_set%lmax, 1, nset)
489 130 : CALL reallocate(gto_basis_set%lmin, 1, nset)
490 130 : CALL reallocate(gto_basis_set%npgf, 1, nset)
491 130 : CALL reallocate(gto_basis_set%nshell, 1, nset)
492 130 : CALL reallocate(gto_basis_set%n, 1, 1, 1, nset)
493 130 : CALL reallocate(gto_basis_set%l, 1, 1, 1, nset)
494 130 : CALL reallocate(gto_basis_set%zet, 1, mprim, 1, nset)
495 130 : CALL reallocate(gto_basis_set%gcc, 1, mprim, 1, 1, 1, nset)
496 :
497 130 : ind_ao = 0
498 130 : maxl = 0
499 420 : DO ishell = 1, nset
500 290 : ang = tb%calc%bas%cgto(ishell, ityp)%ang
501 290 : natorb = natorb + (2*ang + 1)
502 290 : param%lval(ishell) = ang
503 290 : maxl = MAX(ang, maxl)
504 290 : gto_basis_set%lmax(ishell) = ang
505 290 : gto_basis_set%lmin(ishell) = ang
506 290 : gto_basis_set%npgf(ishell) = tb%calc%bas%cgto(ishell, ityp)%nprim
507 290 : gto_basis_set%nshell(ishell) = nshell
508 290 : gto_basis_set%n(1, ishell) = ang + 1
509 290 : gto_basis_set%l(1, ishell) = ang
510 1906 : DO ipgf = 1, gto_basis_set%npgf(ishell)
511 1616 : gto_basis_set%gcc(ipgf, 1, ishell) = tb%calc%bas%cgto(ishell, ityp)%coeff(ipgf)
512 1906 : gto_basis_set%zet(ipgf, ishell) = tb%calc%bas%cgto(ishell, ityp)%alpha(ipgf)
513 : END DO
514 998 : DO ipgf = 1, (2*ang + 1)
515 578 : ind_ao = ind_ao + 1
516 578 : param%lao(ind_ao) = ang
517 868 : param%nao(ind_ao) = ishell
518 : END DO
519 : END DO
520 :
521 130 : do_ortho = .FALSE.
522 130 : CALL process_gto_basis(gto_basis_set, do_ortho, nset, maxl)
523 :
524 : !setting additional values in parameter
525 130 : param%rcut = get_cutoff(tb%calc%bas)
526 130 : param%natorb = natorb
527 130 : param%lmax = maxl !max angular momentum
528 :
529 : !getting occupation
530 130 : occ = 0
531 420 : DO iSh = 1, MIN(tb%calc%bas%nsh_at(id_atom), 5)
532 290 : occ(iSh) = NINT(tb%calc%h0%refocc(iSh, ityp))
533 420 : param%occupation(iSh) = occ(iSh)
534 : END DO
535 780 : param%zeff = SUM(occ) !effective core charge
536 :
537 : !set normalization process
538 130 : gto_basis_set%norm_type = 3
539 :
540 : #else
541 : occ = 0
542 : MARK_USED(tb)
543 : MARK_USED(gto_basis_set)
544 : MARK_USED(element_symbol)
545 : MARK_USED(param)
546 : CPABORT("Built without TBLITE")
547 : #endif
548 :
549 130 : END SUBROUTINE tb_get_basis
550 :
551 : ! **************************************************************************************************
552 : !> \brief ...
553 : !> \param qs_env ...
554 : !> \param calculate_forces ...
555 : ! **************************************************************************************************
556 928 : SUBROUTINE build_tblite_matrices(qs_env, calculate_forces)
557 :
558 : TYPE(qs_environment_type), POINTER :: qs_env
559 : LOGICAL, INTENT(IN) :: calculate_forces
560 :
561 : #if defined(__TBLITE)
562 :
563 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_tblite_matrices'
564 :
565 : INTEGER :: handle, maxder, nderivatives, nimg, img, nkind, i, ic, iw, &
566 : iatom, jatom, ikind, jkind, iset, jset, n1, n2, icol, irow, &
567 : ia, ib, sgfa, sgfb, atom_a, atom_b, &
568 : ldsab, nseta, nsetb, natorb_a, natorb_b, sgfa0
569 : LOGICAL :: found, norml1, norml2, use_arnoldi
570 : REAL(KIND=dp) :: dr, rr
571 : INTEGER, DIMENSION(3) :: cell
572 : REAL(KIND=dp) :: hij, shpoly
573 : REAL(KIND=dp), DIMENSION(2) :: condnum
574 : REAL(KIND=dp), DIMENSION(3) :: rij
575 928 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
576 928 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
577 928 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint, hint
578 928 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min
579 928 : INTEGER, DIMENSION(:), POINTER :: npgfa, npgfb, nsgfa, nsgfb
580 928 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
581 928 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
582 928 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, zeta, zetb, scon_a, scon_b
583 928 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sblock, fblock
584 :
585 928 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
586 : TYPE(atprop_type), POINTER :: atprop
587 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
588 : TYPE(cp_logger_type), POINTER :: logger
589 928 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s, matrix_p, matrix_w
590 : TYPE(dft_control_type), POINTER :: dft_control
591 928 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
592 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
593 928 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
594 928 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
595 : TYPE(neighbor_list_iterator_p_type), &
596 928 : DIMENSION(:), POINTER :: nl_iterator
597 : TYPE(mp_para_env_type), POINTER :: para_env
598 : TYPE(qs_energy_type), POINTER :: energy
599 : TYPE(qs_ks_env_type), POINTER :: ks_env
600 928 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
601 : TYPE(qs_rho_type), POINTER :: rho
602 : TYPE(tblite_type), POINTER :: tb
603 : TYPE(tb_hamiltonian), POINTER :: h0
604 : TYPE(virial_type), POINTER :: virial
605 :
606 928 : CALL timeset(routineN, handle)
607 :
608 928 : NULLIFY (ks_env, energy, atomic_kind_set, qs_kind_set)
609 928 : NULLIFY (matrix_h, matrix_s, atprop, dft_control)
610 928 : NULLIFY (sab_orb, rho, tb)
611 :
612 : CALL get_qs_env(qs_env=qs_env, &
613 : ks_env=ks_env, para_env=para_env, &
614 : energy=energy, &
615 : atomic_kind_set=atomic_kind_set, &
616 : qs_kind_set=qs_kind_set, &
617 : matrix_h_kp=matrix_h, &
618 : matrix_s_kp=matrix_s, &
619 : atprop=atprop, &
620 : dft_control=dft_control, &
621 : sab_orb=sab_orb, &
622 928 : rho=rho, tb_tblite=tb)
623 928 : h0 => tb%calc%h0
624 :
625 : !update geometry (required for debug / geometry optimization)
626 928 : CALL tb_update_geometry(qs_env, tb)
627 :
628 928 : nkind = SIZE(atomic_kind_set)
629 928 : nderivatives = 0
630 928 : IF (calculate_forces) THEN
631 36 : nderivatives = 1
632 36 : IF (ALLOCATED(tb%grad)) DEALLOCATE (tb%grad)
633 108 : ALLOCATE (tb%grad(3, tb%mol%nat))
634 36 : IF (ALLOCATED(tb%dsedcn)) DEALLOCATE (tb%dsedcn)
635 108 : ALLOCATE (tb%dsedcn(tb%calc%bas%nsh))
636 36 : IF (ALLOCATED(tb%calc%ncoord)) THEN
637 36 : IF (ALLOCATED(tb%dcndr)) DEALLOCATE (tb%dcndr)
638 144 : ALLOCATE (tb%dcndr(3, tb%mol%nat, tb%mol%nat))
639 36 : IF (ALLOCATED(tb%dcndL)) DEALLOCATE (tb%dcndL)
640 108 : ALLOCATE (tb%dcndL(3, 3, tb%mol%nat))
641 : END IF
642 : END IF
643 928 : maxder = ncoset(nderivatives)
644 928 : nimg = dft_control%nimages
645 :
646 : !intialise hamiltonian
647 928 : CALL tb_init_ham(qs_env, tb, para_env)
648 :
649 : ! get density matrtix
650 928 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
651 :
652 : ! set up matrices for force calculations
653 928 : IF (calculate_forces) THEN
654 36 : NULLIFY (force, matrix_w, virial)
655 : CALL get_qs_env(qs_env=qs_env, &
656 : matrix_w_kp=matrix_w, &
657 36 : virial=virial, force=force)
658 :
659 36 : IF (SIZE(matrix_p, 1) == 2) THEN
660 0 : DO img = 1, nimg
661 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
662 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
663 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
664 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
665 : END DO
666 : END IF
667 50 : tb%use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
668 : END IF
669 :
670 928 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
671 :
672 : ! set up basis set lists
673 5216 : ALLOCATE (basis_set_list(nkind))
674 928 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
675 :
676 : ! allocate overlap matrix
677 928 : CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
678 : CALL create_sab_matrix(ks_env, matrix_s, "OVERLAP MATRIX", basis_set_list, basis_set_list, &
679 928 : sab_orb, .TRUE.)
680 928 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
681 :
682 : ! initialize H matrix
683 928 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
684 1856 : DO img = 1, nimg
685 928 : ALLOCATE (matrix_h(1, img)%matrix)
686 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
687 928 : name="HAMILTONIAN MATRIX")
688 1856 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
689 : END DO
690 928 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
691 :
692 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
693 928 : NULLIFY (nl_iterator)
694 928 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
695 166548 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
696 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
697 165620 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
698 :
699 165620 : icol = MAX(iatom, jatom)
700 165620 : irow = MIN(iatom, jatom)
701 165620 : IF (icol == jatom) THEN
702 377312 : rij = -rij
703 94328 : i = ikind
704 94328 : ikind = jkind
705 94328 : jkind = i
706 : END IF
707 :
708 662480 : dr = NORM2(rij(:))
709 :
710 165620 : ic = 1
711 165620 : NULLIFY (sblock)
712 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
713 165620 : row=irow, col=icol, BLOCK=sblock, found=found)
714 165620 : CPASSERT(found)
715 165620 : NULLIFY (fblock)
716 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
717 165620 : row=irow, col=icol, BLOCK=fblock, found=found)
718 165620 : CPASSERT(found)
719 :
720 : ! --------- Overlap
721 : !get basis information
722 165620 : basis_set_a => basis_set_list(ikind)%gto_basis_set
723 165620 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
724 165620 : basis_set_b => basis_set_list(jkind)%gto_basis_set
725 165620 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
726 165620 : atom_a = atom_of_kind(icol)
727 165620 : atom_b = atom_of_kind(irow)
728 : ! basis a
729 165620 : first_sgfa => basis_set_a%first_sgf
730 165620 : la_max => basis_set_a%lmax
731 165620 : la_min => basis_set_a%lmin
732 165620 : npgfa => basis_set_a%npgf
733 165620 : nseta = basis_set_a%nset
734 165620 : nsgfa => basis_set_a%nsgf_set
735 165620 : rpgfa => basis_set_a%pgf_radius
736 165620 : set_radius_a => basis_set_a%set_radius
737 165620 : scon_a => basis_set_a%scon
738 165620 : zeta => basis_set_a%zet
739 : ! basis b
740 165620 : first_sgfb => basis_set_b%first_sgf
741 165620 : lb_max => basis_set_b%lmax
742 165620 : lb_min => basis_set_b%lmin
743 165620 : npgfb => basis_set_b%npgf
744 165620 : nsetb = basis_set_b%nset
745 165620 : nsgfb => basis_set_b%nsgf_set
746 165620 : rpgfb => basis_set_b%pgf_radius
747 165620 : set_radius_b => basis_set_b%set_radius
748 165620 : scon_b => basis_set_b%scon
749 165620 : zetb => basis_set_b%zet
750 :
751 165620 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
752 1324960 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
753 165620 : natorb_a = 0
754 643037 : DO iset = 1, nseta
755 643037 : natorb_a = natorb_a + (2*basis_set_a%l(1, iset) + 1)
756 : END DO
757 : natorb_b = 0
758 643631 : DO iset = 1, nsetb
759 643631 : natorb_b = natorb_b + (2*basis_set_b%l(1, iset) + 1)
760 : END DO
761 828100 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
762 14628036 : sint = 0.0_dp
763 662480 : ALLOCATE (hint(natorb_a, natorb_b, maxder))
764 14628036 : hint = 0.0_dp
765 :
766 : !----------------- overlap integrals
767 643037 : DO iset = 1, nseta
768 477417 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
769 477417 : sgfa = first_sgfa(1, iset)
770 2037788 : DO jset = 1, nsetb
771 1394751 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
772 1021871 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
773 1021871 : sgfb = first_sgfb(1, jset)
774 1021871 : IF (calculate_forces) THEN
775 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
776 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
777 21607 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
778 : ELSE
779 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
780 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
781 1000264 : rij, sab=oint(:, :, 1))
782 : END IF
783 : ! Contraction
784 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
785 1021871 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
786 1021871 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
787 1499288 : IF (calculate_forces) THEN
788 86428 : DO i = 2, 4
789 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
790 64821 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
791 86428 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
792 : END DO
793 : END IF
794 : END DO
795 : END DO
796 :
797 : ! update S matrix
798 165620 : IF (icol <= irow) THEN
799 3189395 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
800 : ELSE
801 10436767 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
802 : END IF
803 :
804 : ! --------- Hamiltonian
805 165620 : IF (icol == irow .AND. dr < same_atom) THEN
806 : !get diagonal F matrix from selfenergy
807 2117 : n1 = tb%calc%bas%ish_at(icol)
808 7061 : DO iset = 1, nseta
809 4944 : sgfa = first_sgfa(1, iset)
810 4944 : hij = tb%selfenergy(n1 + iset)
811 17967 : DO ia = sgfa, sgfa + nsgfa(iset) - 1
812 15850 : hint(ia, ia, 1) = hij
813 : END DO
814 : END DO
815 : ELSE
816 : !get off-diagonal F matrix
817 163503 : rr = SQRT(dr/(h0%rad(jkind) + h0%rad(ikind)))
818 163503 : n1 = tb%calc%bas%ish_at(icol)
819 163503 : sgfa0 = 1
820 635976 : DO iset = 1, nseta
821 472473 : sgfa = first_sgfa(1, iset)
822 472473 : sgfa0 = sgfa0 + tb%calc%bas%cgto(iset, tb%mol%id(icol))%nprim
823 472473 : n2 = tb%calc%bas%ish_at(irow)
824 2018513 : DO jset = 1, nsetb
825 1382537 : sgfb = first_sgfb(1, jset)
826 : shpoly = (1.0_dp + h0%shpoly(iset, ikind)*rr) &
827 1382537 : *(1.0_dp + h0%shpoly(jset, jkind)*rr)
828 : hij = 0.5_dp*(tb%selfenergy(n1 + iset) + tb%selfenergy(n2 + jset)) &
829 1382537 : *h0%hscale(iset, jset, ikind, jkind)*shpoly
830 5902375 : DO ia = sgfa, sgfa + nsgfa(iset) - 1
831 17439883 : DO ib = sgfb, sgfb + nsgfb(jset) - 1
832 16057346 : hint(ia, ib, 1) = hij*sint(ia, ib, 1)
833 : END DO
834 : END DO
835 : END DO
836 : END DO
837 : END IF
838 :
839 : ! update F matrix
840 165620 : IF (icol <= irow) THEN
841 3189395 : fblock(:, :) = fblock(:, :) + hint(:, :, 1)
842 : ELSE
843 10436767 : fblock(:, :) = fblock(:, :) + TRANSPOSE(hint(:, :, 1))
844 : END IF
845 :
846 497788 : DEALLOCATE (oint, owork, sint, hint)
847 :
848 : END DO
849 928 : CALL neighbor_list_iterator_release(nl_iterator)
850 :
851 1856 : DO img = 1, nimg
852 1964 : DO i = 1, SIZE(matrix_s, 1)
853 1964 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
854 : END DO
855 2784 : DO i = 1, SIZE(matrix_h, 1)
856 1856 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
857 : END DO
858 : END DO
859 :
860 : !compute multipole moments for gfn2
861 928 : IF (dft_control%qs_control%xtb_control%tblite_method == gfn2xtb) &
862 234 : CALL tb_get_multipole(qs_env, tb)
863 :
864 : ! output overlap information
865 928 : NULLIFY (logger)
866 928 : logger => cp_get_default_logger()
867 928 : IF (.NOT. calculate_forces) THEN
868 892 : IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
869 : "DFT%PRINT%OVERLAP_CONDITION") /= 0) THEN
870 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
871 2 : extension=".Log")
872 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
873 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
874 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
875 2 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
876 2 : CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
877 : END IF
878 : END IF
879 :
880 928 : DEALLOCATE (basis_set_list)
881 :
882 928 : CALL timestop(handle)
883 :
884 : #else
885 : MARK_USED(qs_env)
886 : MARK_USED(calculate_forces)
887 : CPABORT("Built without TBLITE")
888 : #endif
889 :
890 1856 : END SUBROUTINE build_tblite_matrices
891 :
892 : ! **************************************************************************************************
893 : !> \brief ...
894 : !> \param qs_env ...
895 : !> \param dft_control ...
896 : !> \param tb ...
897 : !> \param calculate_forces ...
898 : !> \param use_rho ...
899 : ! **************************************************************************************************
900 23506 : SUBROUTINE tb_update_charges(qs_env, dft_control, tb, calculate_forces, use_rho)
901 :
902 : TYPE(qs_environment_type), POINTER :: qs_env
903 : TYPE(dft_control_type), POINTER :: dft_control
904 : TYPE(tblite_type), POINTER :: tb
905 : LOGICAL, INTENT(IN) :: calculate_forces
906 : LOGICAL, INTENT(IN) :: use_rho
907 :
908 : #if defined(__TBLITE)
909 :
910 : INTEGER :: iatom, ikind, is, ns, atom_a, ii, im
911 : INTEGER :: nimg, nkind, nsgf, natorb, na
912 : INTEGER :: n_atom, max_orb, max_shell
913 : LOGICAL :: do_dipole, do_quadrupole
914 : REAL(KIND=dp) :: norm
915 : INTEGER, DIMENSION(5) :: occ
916 : INTEGER, DIMENSION(25) :: lao
917 : INTEGER, DIMENSION(25) :: nao
918 23506 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ch_atom, ch_shell, ch_ref, ch_orb
919 23506 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, ao_dip, ao_quad
920 :
921 23506 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
922 23506 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_p
923 23506 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
924 : TYPE(dbcsr_type), POINTER :: s_matrix
925 : TYPE(mp_para_env_type), POINTER :: para_env
926 23506 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
927 23506 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
928 : TYPE(qs_rho_type), POINTER :: rho
929 : TYPE(qs_scf_env_type), POINTER :: scf_env
930 : TYPE(xtb_atom_type), POINTER :: xtb_kind
931 :
932 : ! also compute multipoles needed by GFN2
933 23506 : do_dipole = .FALSE.
934 23506 : do_quadrupole = .FALSE.
935 :
936 : ! compute mulliken charges required for charge update
937 23506 : NULLIFY (particle_set, qs_kind_set, atomic_kind_set)
938 : CALL get_qs_env(qs_env=qs_env, scf_env=scf_env, particle_set=particle_set, qs_kind_set=qs_kind_set, &
939 23506 : atomic_kind_set=atomic_kind_set, matrix_s_kp=matrix_s, rho=rho, para_env=para_env)
940 23506 : NULLIFY (matrix_p)
941 23506 : IF (use_rho) THEN
942 12186 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
943 12186 : IF (ASSOCIATED(tb%dipbra)) do_dipole = .TRUE.
944 12186 : IF (ASSOCIATED(tb%quadbra)) do_quadrupole = .TRUE.
945 : ELSE
946 11320 : matrix_p => scf_env%p_mix_new
947 : END IF
948 23506 : n_atom = SIZE(particle_set)
949 23506 : nkind = SIZE(atomic_kind_set)
950 23506 : nimg = dft_control%nimages
951 23506 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
952 94024 : ALLOCATE (aocg(nsgf, n_atom))
953 806956 : aocg = 0.0_dp
954 32906 : IF (do_dipole) ALLOCATE (ao_dip(n_atom, dip_n))
955 32906 : IF (do_quadrupole) ALLOCATE (ao_quad(n_atom, quad_n))
956 23506 : max_orb = 0
957 23506 : max_shell = 0
958 87544 : DO ikind = 1, nkind
959 64038 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
960 64038 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb)
961 87544 : max_orb = MAX(max_orb, natorb)
962 : END DO
963 126796 : DO is = 1, n_atom
964 126796 : max_shell = MAX(max_shell, tb%calc%bas%nsh_at(is))
965 : END DO
966 141036 : ALLOCATE (ch_atom(n_atom, 1), ch_shell(n_atom, max_shell))
967 141036 : ALLOCATE (ch_orb(max_orb, n_atom), ch_ref(max_orb, n_atom))
968 150302 : ch_atom = 0.0_dp
969 355218 : ch_shell = 0.0_dp
970 806956 : ch_orb = 0.0_dp
971 806956 : ch_ref = 0.0_dp
972 23506 : IF (nimg > 1) THEN
973 0 : CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
974 0 : IF (do_dipole .OR. do_quadrupole) THEN
975 0 : CPABORT("missing contraction with density matrix for multiple k-points")
976 : END IF
977 : ELSE
978 23506 : NULLIFY (p_matrix, s_matrix)
979 23506 : p_matrix => matrix_p(:, 1)
980 23506 : s_matrix => matrix_s(1, 1)%matrix
981 23506 : CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
982 23506 : IF (do_dipole) THEN
983 18800 : DO im = 1, dip_n
984 18800 : CALL contract_dens(p_matrix, tb%dipbra(im)%matrix, tb%dipket(im)%matrix, ao_dip(:, im), para_env)
985 : END DO
986 : END IF
987 23506 : IF (do_quadrupole) THEN
988 32900 : DO im = 1, quad_n
989 32900 : CALL contract_dens(p_matrix, tb%quadbra(im)%matrix, tb%quadket(im)%matrix, ao_quad(:, im), para_env)
990 : END DO
991 : END IF
992 : END IF
993 23506 : NULLIFY (xtb_kind)
994 87544 : DO ikind = 1, nkind
995 64038 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
996 64038 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
997 64038 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, nao=nao, occupation=occ)
998 254872 : DO iatom = 1, na
999 103290 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1000 603288 : DO is = 1, natorb
1001 499998 : ns = lao(is) + 1
1002 499998 : norm = 2*lao(is) + 1
1003 499998 : ch_ref(is, atom_a) = tb%calc%h0%refocc(nao(is), ikind)/norm
1004 499998 : ch_orb(is, atom_a) = aocg(is, atom_a) - ch_ref(is, atom_a)
1005 603288 : ch_shell(atom_a, ns) = ch_orb(is, atom_a) + ch_shell(atom_a, ns)
1006 : END DO
1007 847488 : ch_atom(atom_a, 1) = SUM(ch_orb(:, atom_a))
1008 : END DO
1009 : END DO
1010 23506 : DEALLOCATE (aocg)
1011 :
1012 : ! charge mixing
1013 23506 : IF (dft_control%qs_control%do_ls_scf) THEN
1014 : !
1015 : ELSE
1016 : CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
1017 23506 : ch_shell, para_env, scf_env%iter_count)
1018 : END IF
1019 :
1020 : !setting new wave function
1021 23506 : CALL tb%pot%reset
1022 126796 : tb%e_es = 0.0_dp
1023 126796 : tb%e_scd = 0.0_dp
1024 126796 : DO iatom = 1, n_atom
1025 103290 : ii = tb%calc%bas%ish_at(iatom)
1026 337932 : DO is = 1, tb%calc%bas%nsh_at(iatom)
1027 337932 : tb%wfn%qsh(ii + is, 1) = -ch_shell(iatom, is)
1028 : END DO
1029 126796 : tb%wfn%qat(iatom, 1) = -ch_atom(iatom, 1)
1030 : END DO
1031 :
1032 23506 : IF (do_dipole) THEN
1033 21318 : DO iatom = 1, n_atom
1034 71172 : DO im = 1, dip_n
1035 66472 : tb%wfn%dpat(im, iatom, 1) = -ao_dip(iatom, im)
1036 : END DO
1037 : END DO
1038 4700 : DEALLOCATE (ao_dip)
1039 : END IF
1040 23506 : IF (do_quadrupole) THEN
1041 21318 : DO iatom = 1, n_atom
1042 121026 : DO im = 1, quad_n
1043 116326 : tb%wfn%qpat(im, iatom, 1) = -ao_quad(iatom, im)
1044 : END DO
1045 : END DO
1046 4700 : DEALLOCATE (ao_quad)
1047 : END IF
1048 :
1049 23506 : IF (ALLOCATED(tb%calc%coulomb)) THEN
1050 23506 : CALL tb%calc%coulomb%get_potential(tb%mol, tb%cache, tb%wfn, tb%pot)
1051 23506 : CALL tb%calc%coulomb%get_energy(tb%mol, tb%cache, tb%wfn, tb%e_es)
1052 : END IF
1053 23506 : IF (ALLOCATED(tb%calc%dispersion)) THEN
1054 23506 : CALL tb%calc%dispersion%get_potential(tb%mol, tb%dcache, tb%wfn, tb%pot)
1055 23506 : CALL tb%calc%dispersion%get_energy(tb%mol, tb%dcache, tb%wfn, tb%e_scd)
1056 : END IF
1057 :
1058 23506 : IF (calculate_forces) THEN
1059 36 : IF (ALLOCATED(tb%calc%coulomb)) THEN
1060 636 : tb%grad = 0.0_dp
1061 36 : CALL tb%calc%coulomb%get_gradient(tb%mol, tb%cache, tb%wfn, tb%grad, tb%sigma)
1062 36 : CALL tb_grad2force(qs_env, tb, para_env, 3)
1063 : END IF
1064 :
1065 36 : IF (ALLOCATED(tb%calc%dispersion)) THEN
1066 636 : tb%grad = 0.0_dp
1067 36 : CALL tb%calc%dispersion%get_gradient(tb%mol, tb%dcache, tb%wfn, tb%grad, tb%sigma)
1068 36 : CALL tb_grad2force(qs_env, tb, para_env, 2)
1069 : END IF
1070 : END IF
1071 :
1072 23506 : DEALLOCATE (ch_atom, ch_shell, ch_orb, ch_ref)
1073 :
1074 : #else
1075 : MARK_USED(qs_env)
1076 : MARK_USED(tb)
1077 : MARK_USED(dft_control)
1078 : MARK_USED(calculate_forces)
1079 : MARK_USED(use_rho)
1080 : CPABORT("Built without TBLITE")
1081 : #endif
1082 :
1083 70518 : END SUBROUTINE tb_update_charges
1084 :
1085 : ! **************************************************************************************************
1086 : !> \brief ...
1087 : !> \param qs_env ...
1088 : !> \param tb ...
1089 : !> \param dft_control ...
1090 : ! **************************************************************************************************
1091 12186 : SUBROUTINE tb_ham_add_coulomb(qs_env, tb, dft_control)
1092 :
1093 : TYPE(qs_environment_type), POINTER :: qs_env
1094 : TYPE(tblite_type), POINTER :: tb
1095 : TYPE(dft_control_type), POINTER :: dft_control
1096 :
1097 : #if defined(__TBLITE)
1098 :
1099 : INTEGER :: ikind, jkind, iatom, jatom, icol, irow
1100 : INTEGER :: ic, is, nimg, ni, nj, i, j
1101 : INTEGER :: la, lb, za, zb
1102 : LOGICAL :: found
1103 : INTEGER, DIMENSION(3) :: cellind
1104 : INTEGER, DIMENSION(25) :: naoa, naob
1105 : REAL(KIND=dp), DIMENSION(3) :: rij
1106 12186 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, sum_shell
1107 12186 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ashift, bshift
1108 12186 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksblock, sblock
1109 12186 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1110 12186 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dip_ket1, dip_ket2, dip_ket3, &
1111 12186 : dip_bra1, dip_bra2, dip_bra3
1112 12186 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: quad_ket1, quad_ket2, quad_ket3, &
1113 12186 : quad_ket4, quad_ket5, quad_ket6, &
1114 12186 : quad_bra1, quad_bra2, quad_bra3, &
1115 12186 : quad_bra4, quad_bra5, quad_bra6
1116 :
1117 12186 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1118 : TYPE(dbcsr_iterator_type) :: iter
1119 12186 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1120 12186 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
1121 : TYPE(kpoint_type), POINTER :: kpoints
1122 : TYPE(neighbor_list_iterator_p_type), &
1123 12186 : DIMENSION(:), POINTER :: nl_iterator
1124 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1125 12186 : POINTER :: n_list
1126 12186 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1127 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
1128 :
1129 12186 : nimg = dft_control%nimages
1130 :
1131 12186 : NULLIFY (matrix_s, ks_matrix, n_list, qs_kind_set)
1132 12186 : 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)
1133 :
1134 : !creating sum of shell lists
1135 36558 : ALLOCATE (sum_shell(tb%mol%nat))
1136 65812 : i = 0
1137 65812 : DO j = 1, tb%mol%nat
1138 53626 : sum_shell(j) = i
1139 65812 : i = i + tb%calc%bas%nsh_at(j)
1140 : END DO
1141 :
1142 12186 : IF (nimg == 1) THEN
1143 : ! no k-points; all matrices have been transformed to periodic bsf
1144 12186 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1145 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1146 12186 : kind_of=kind_of)
1147 12186 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
1148 93572 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1149 81386 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1150 :
1151 81386 : ikind = kind_of(irow)
1152 81386 : jkind = kind_of(icol)
1153 :
1154 : ! atomic parameters
1155 81386 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
1156 81386 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
1157 81386 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa)
1158 81386 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob)
1159 :
1160 81386 : ni = SIZE(sblock, 1)
1161 325544 : ALLOCATE (ashift(ni, ni))
1162 3576932 : ashift = 0.0_dp
1163 533309 : DO i = 1, ni
1164 451923 : la = naoa(i) + sum_shell(irow)
1165 533309 : ashift(i, i) = tb%pot%vsh(la, 1)
1166 : END DO
1167 :
1168 81386 : nj = SIZE(sblock, 2)
1169 325544 : ALLOCATE (bshift(nj, nj))
1170 3112656 : bshift = 0.0_dp
1171 481341 : DO j = 1, nj
1172 399955 : lb = naob(j) + sum_shell(icol)
1173 481341 : bshift(j, j) = tb%pot%vsh(lb, 1)
1174 : END DO
1175 :
1176 162772 : DO is = 1, SIZE(ks_matrix, 1)
1177 81386 : NULLIFY (ksblock)
1178 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
1179 81386 : row=irow, col=icol, block=ksblock, found=found)
1180 81386 : CPASSERT(found)
1181 325544 : ksblock = ksblock - 0.5_dp*(MATMUL(ashift, sblock) &
1182 176780140 : + MATMUL(sblock, bshift))
1183 : ksblock = ksblock - 0.5_dp*(tb%pot%vat(irow, 1) &
1184 6472068 : + tb%pot%vat(icol, 1))*sblock
1185 : END DO
1186 256344 : DEALLOCATE (ashift, bshift)
1187 : END DO
1188 12186 : CALL dbcsr_iterator_stop(iter)
1189 :
1190 12186 : IF (ASSOCIATED(tb%dipbra)) THEN
1191 4700 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
1192 24010 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1193 19310 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1194 :
1195 19310 : NULLIFY (dip_bra1, dip_bra2, dip_bra3)
1196 : CALL dbcsr_get_block_p(matrix=tb%dipbra(1)%matrix, &
1197 19310 : row=irow, col=icol, BLOCK=dip_bra1, found=found)
1198 19310 : CPASSERT(found)
1199 : CALL dbcsr_get_block_p(matrix=tb%dipbra(2)%matrix, &
1200 19310 : row=irow, col=icol, BLOCK=dip_bra2, found=found)
1201 19310 : CPASSERT(found)
1202 : CALL dbcsr_get_block_p(matrix=tb%dipbra(3)%matrix, &
1203 19310 : row=irow, col=icol, BLOCK=dip_bra3, found=found)
1204 19310 : CPASSERT(found)
1205 19310 : NULLIFY (dip_ket1, dip_ket2, dip_ket3)
1206 : CALL dbcsr_get_block_p(matrix=tb%dipket(1)%matrix, &
1207 19310 : row=irow, col=icol, BLOCK=dip_ket1, found=found)
1208 19310 : CPASSERT(found)
1209 : CALL dbcsr_get_block_p(matrix=tb%dipket(2)%matrix, &
1210 19310 : row=irow, col=icol, BLOCK=dip_ket2, found=found)
1211 19310 : CPASSERT(found)
1212 : CALL dbcsr_get_block_p(matrix=tb%dipket(3)%matrix, &
1213 19310 : row=irow, col=icol, BLOCK=dip_ket3, found=found)
1214 19310 : CPASSERT(found)
1215 :
1216 43320 : DO is = 1, SIZE(ks_matrix, 1)
1217 19310 : NULLIFY (ksblock)
1218 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
1219 19310 : row=irow, col=icol, block=ksblock, found=found)
1220 19310 : CPASSERT(found)
1221 : ksblock = ksblock - 0.5_dp*(dip_ket1*tb%pot%vdp(1, irow, 1) &
1222 : + dip_ket2*tb%pot%vdp(2, irow, 1) &
1223 : + dip_ket3*tb%pot%vdp(3, irow, 1) &
1224 : + dip_bra1*tb%pot%vdp(1, icol, 1) &
1225 : + dip_bra2*tb%pot%vdp(2, icol, 1) &
1226 968084 : + dip_bra3*tb%pot%vdp(3, icol, 1))
1227 : END DO
1228 : END DO
1229 4700 : CALL dbcsr_iterator_stop(iter)
1230 : END IF
1231 :
1232 12186 : IF (ASSOCIATED(tb%quadbra)) THEN
1233 4700 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
1234 24010 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1235 19310 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1236 :
1237 19310 : NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
1238 : CALL dbcsr_get_block_p(matrix=tb%quadbra(1)%matrix, &
1239 19310 : row=irow, col=icol, BLOCK=quad_bra1, found=found)
1240 19310 : CPASSERT(found)
1241 : CALL dbcsr_get_block_p(matrix=tb%quadbra(2)%matrix, &
1242 19310 : row=irow, col=icol, BLOCK=quad_bra2, found=found)
1243 19310 : CPASSERT(found)
1244 : CALL dbcsr_get_block_p(matrix=tb%quadbra(3)%matrix, &
1245 19310 : row=irow, col=icol, BLOCK=quad_bra3, found=found)
1246 19310 : CPASSERT(found)
1247 : CALL dbcsr_get_block_p(matrix=tb%quadbra(4)%matrix, &
1248 19310 : row=irow, col=icol, BLOCK=quad_bra4, found=found)
1249 19310 : CPASSERT(found)
1250 : CALL dbcsr_get_block_p(matrix=tb%quadbra(5)%matrix, &
1251 19310 : row=irow, col=icol, BLOCK=quad_bra5, found=found)
1252 19310 : CPASSERT(found)
1253 : CALL dbcsr_get_block_p(matrix=tb%quadbra(6)%matrix, &
1254 19310 : row=irow, col=icol, BLOCK=quad_bra6, found=found)
1255 19310 : CPASSERT(found)
1256 :
1257 19310 : NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
1258 : CALL dbcsr_get_block_p(matrix=tb%quadket(1)%matrix, &
1259 19310 : row=irow, col=icol, BLOCK=quad_ket1, found=found)
1260 19310 : CPASSERT(found)
1261 : CALL dbcsr_get_block_p(matrix=tb%quadket(2)%matrix, &
1262 19310 : row=irow, col=icol, BLOCK=quad_ket2, found=found)
1263 19310 : CPASSERT(found)
1264 : CALL dbcsr_get_block_p(matrix=tb%quadket(3)%matrix, &
1265 19310 : row=irow, col=icol, BLOCK=quad_ket3, found=found)
1266 19310 : CPASSERT(found)
1267 : CALL dbcsr_get_block_p(matrix=tb%quadket(4)%matrix, &
1268 19310 : row=irow, col=icol, BLOCK=quad_ket4, found=found)
1269 19310 : CPASSERT(found)
1270 : CALL dbcsr_get_block_p(matrix=tb%quadket(5)%matrix, &
1271 19310 : row=irow, col=icol, BLOCK=quad_ket5, found=found)
1272 19310 : CPASSERT(found)
1273 : CALL dbcsr_get_block_p(matrix=tb%quadket(6)%matrix, &
1274 19310 : row=irow, col=icol, BLOCK=quad_ket6, found=found)
1275 19310 : CPASSERT(found)
1276 :
1277 43320 : DO is = 1, SIZE(ks_matrix, 1)
1278 19310 : NULLIFY (ksblock)
1279 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
1280 19310 : row=irow, col=icol, block=ksblock, found=found)
1281 19310 : CPASSERT(found)
1282 :
1283 : ksblock = ksblock - 0.5_dp*(quad_ket1*tb%pot%vqp(1, irow, 1) &
1284 : + quad_ket2*tb%pot%vqp(2, irow, 1) &
1285 : + quad_ket3*tb%pot%vqp(3, irow, 1) &
1286 : + quad_ket4*tb%pot%vqp(4, irow, 1) &
1287 : + quad_ket5*tb%pot%vqp(5, irow, 1) &
1288 : + quad_ket6*tb%pot%vqp(6, irow, 1) &
1289 : + quad_bra1*tb%pot%vqp(1, icol, 1) &
1290 : + quad_bra2*tb%pot%vqp(2, icol, 1) &
1291 : + quad_bra3*tb%pot%vqp(3, icol, 1) &
1292 : + quad_bra4*tb%pot%vqp(4, icol, 1) &
1293 : + quad_bra5*tb%pot%vqp(5, icol, 1) &
1294 968084 : + quad_bra6*tb%pot%vqp(6, icol, 1))
1295 : END DO
1296 : END DO
1297 4700 : CALL dbcsr_iterator_stop(iter)
1298 : END IF
1299 :
1300 : ELSE
1301 0 : CPABORT("GFN methods with k-points not tested")
1302 0 : NULLIFY (kpoints)
1303 0 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
1304 0 : NULLIFY (cell_to_index)
1305 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1306 :
1307 0 : NULLIFY (nl_iterator)
1308 0 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
1309 0 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1310 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
1311 0 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
1312 :
1313 0 : icol = MAX(iatom, jatom)
1314 0 : irow = MIN(iatom, jatom)
1315 :
1316 0 : IF (icol == jatom) THEN
1317 0 : rij = -rij
1318 0 : i = ikind
1319 0 : ikind = jkind
1320 0 : jkind = i
1321 : END IF
1322 :
1323 0 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
1324 0 : CPASSERT(ic > 0)
1325 :
1326 0 : NULLIFY (sblock)
1327 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
1328 0 : row=irow, col=icol, block=sblock, found=found)
1329 0 : CPASSERT(found)
1330 :
1331 : ! atomic parameters
1332 0 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
1333 0 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
1334 0 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa)
1335 0 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob)
1336 :
1337 0 : ni = SIZE(sblock, 1)
1338 0 : ALLOCATE (ashift(ni, ni))
1339 0 : ashift = 0.0_dp
1340 0 : DO i = 1, ni
1341 0 : la = naoa(i) + sum_shell(irow)
1342 0 : ashift(i, i) = tb%pot%vsh(la, 1)
1343 : END DO
1344 :
1345 0 : nj = SIZE(sblock, 2)
1346 0 : ALLOCATE (bshift(nj, nj))
1347 0 : bshift = 0.0_dp
1348 0 : DO j = 1, nj
1349 0 : lb = naob(j) + sum_shell(icol)
1350 0 : bshift(j, j) = tb%pot%vsh(lb, 1)
1351 : END DO
1352 :
1353 0 : DO is = 1, SIZE(ks_matrix, 1)
1354 0 : NULLIFY (ksblock)
1355 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
1356 0 : row=irow, col=icol, block=ksblock, found=found)
1357 0 : CPASSERT(found)
1358 0 : ksblock = ksblock - 0.5_dp*(MATMUL(ashift, sblock) &
1359 0 : + MATMUL(sblock, bshift))
1360 : ksblock = ksblock - 0.5_dp*(tb%pot%vat(irow, 1) &
1361 0 : + tb%pot%vat(icol, 1))*sblock
1362 : END DO
1363 : END DO
1364 0 : CALL neighbor_list_iterator_release(nl_iterator)
1365 :
1366 0 : IF (ASSOCIATED(tb%dipbra)) THEN
1367 0 : NULLIFY (nl_iterator)
1368 0 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
1369 0 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1370 : CALL get_iterator_info(nl_iterator, &
1371 0 : iatom=iatom, jatom=jatom, cell=cellind)
1372 0 : icol = MAX(iatom, jatom)
1373 0 : irow = MIN(iatom, jatom)
1374 :
1375 0 : NULLIFY (dip_bra1, dip_bra2, dip_bra3)
1376 : CALL dbcsr_get_block_p(matrix=tb%dipbra(1)%matrix, &
1377 0 : row=irow, col=icol, BLOCK=dip_bra1, found=found)
1378 0 : CPASSERT(found)
1379 : CALL dbcsr_get_block_p(matrix=tb%dipbra(2)%matrix, &
1380 0 : row=irow, col=icol, BLOCK=dip_bra2, found=found)
1381 0 : CPASSERT(found)
1382 : CALL dbcsr_get_block_p(matrix=tb%dipbra(3)%matrix, &
1383 0 : row=irow, col=icol, BLOCK=dip_bra3, found=found)
1384 0 : CPASSERT(found)
1385 0 : NULLIFY (dip_ket1, dip_ket2, dip_ket3)
1386 : CALL dbcsr_get_block_p(matrix=tb%dipket(1)%matrix, &
1387 0 : row=irow, col=icol, BLOCK=dip_ket1, found=found)
1388 0 : CPASSERT(found)
1389 : CALL dbcsr_get_block_p(matrix=tb%dipket(2)%matrix, &
1390 0 : row=irow, col=icol, BLOCK=dip_ket2, found=found)
1391 0 : CPASSERT(found)
1392 : CALL dbcsr_get_block_p(matrix=tb%dipket(3)%matrix, &
1393 0 : row=irow, col=icol, BLOCK=dip_ket3, found=found)
1394 0 : CPASSERT(found)
1395 :
1396 0 : DO is = 1, SIZE(ks_matrix, 1)
1397 0 : NULLIFY (ksblock)
1398 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
1399 0 : row=irow, col=icol, block=ksblock, found=found)
1400 0 : CPASSERT(found)
1401 : ksblock = ksblock - 0.5_dp*(dip_ket1*tb%pot%vdp(1, irow, 1) &
1402 : + dip_ket2*tb%pot%vdp(2, irow, 1) &
1403 : + dip_ket3*tb%pot%vdp(3, irow, 1) &
1404 : + dip_bra1*tb%pot%vdp(1, icol, 1) &
1405 : + dip_bra2*tb%pot%vdp(2, icol, 1) &
1406 0 : + dip_bra3*tb%pot%vdp(3, icol, 1))
1407 : END DO
1408 : END DO
1409 0 : CALL neighbor_list_iterator_release(nl_iterator)
1410 : END IF
1411 :
1412 0 : IF (ASSOCIATED(tb%quadbra)) THEN
1413 0 : NULLIFY (nl_iterator)
1414 0 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
1415 0 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1416 : CALL get_iterator_info(nl_iterator, &
1417 0 : iatom=iatom, jatom=jatom, cell=cellind)
1418 0 : icol = MAX(iatom, jatom)
1419 0 : irow = MIN(iatom, jatom)
1420 :
1421 0 : NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
1422 : CALL dbcsr_get_block_p(matrix=tb%quadbra(1)%matrix, &
1423 0 : row=irow, col=icol, BLOCK=quad_bra1, found=found)
1424 0 : CPASSERT(found)
1425 : CALL dbcsr_get_block_p(matrix=tb%quadbra(2)%matrix, &
1426 0 : row=irow, col=icol, BLOCK=quad_bra2, found=found)
1427 0 : CPASSERT(found)
1428 : CALL dbcsr_get_block_p(matrix=tb%quadbra(3)%matrix, &
1429 0 : row=irow, col=icol, BLOCK=quad_bra3, found=found)
1430 0 : CPASSERT(found)
1431 : CALL dbcsr_get_block_p(matrix=tb%quadbra(4)%matrix, &
1432 0 : row=irow, col=icol, BLOCK=quad_bra4, found=found)
1433 0 : CPASSERT(found)
1434 : CALL dbcsr_get_block_p(matrix=tb%quadbra(5)%matrix, &
1435 0 : row=irow, col=icol, BLOCK=quad_bra5, found=found)
1436 0 : CPASSERT(found)
1437 : CALL dbcsr_get_block_p(matrix=tb%quadbra(6)%matrix, &
1438 0 : row=irow, col=icol, BLOCK=quad_bra6, found=found)
1439 0 : CPASSERT(found)
1440 :
1441 0 : NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
1442 : CALL dbcsr_get_block_p(matrix=tb%quadket(1)%matrix, &
1443 0 : row=irow, col=icol, BLOCK=quad_ket1, found=found)
1444 0 : CPASSERT(found)
1445 : CALL dbcsr_get_block_p(matrix=tb%quadket(2)%matrix, &
1446 0 : row=irow, col=icol, BLOCK=quad_ket2, found=found)
1447 0 : CPASSERT(found)
1448 : CALL dbcsr_get_block_p(matrix=tb%quadket(3)%matrix, &
1449 0 : row=irow, col=icol, BLOCK=quad_ket3, found=found)
1450 0 : CPASSERT(found)
1451 : CALL dbcsr_get_block_p(matrix=tb%quadket(4)%matrix, &
1452 0 : row=irow, col=icol, BLOCK=quad_ket4, found=found)
1453 0 : CPASSERT(found)
1454 : CALL dbcsr_get_block_p(matrix=tb%quadket(5)%matrix, &
1455 0 : row=irow, col=icol, BLOCK=quad_ket5, found=found)
1456 0 : CPASSERT(found)
1457 : CALL dbcsr_get_block_p(matrix=tb%quadket(6)%matrix, &
1458 0 : row=irow, col=icol, BLOCK=quad_ket6, found=found)
1459 0 : CPASSERT(found)
1460 :
1461 0 : DO is = 1, SIZE(ks_matrix, 1)
1462 0 : NULLIFY (ksblock)
1463 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
1464 0 : row=irow, col=icol, block=ksblock, found=found)
1465 0 : CPASSERT(found)
1466 :
1467 : ksblock = ksblock - 0.5_dp*(quad_ket1*tb%pot%vqp(1, irow, 1) &
1468 : + quad_ket2*tb%pot%vqp(2, irow, 1) &
1469 : + quad_ket3*tb%pot%vqp(3, irow, 1) &
1470 : + quad_ket4*tb%pot%vqp(4, irow, 1) &
1471 : + quad_ket5*tb%pot%vqp(5, irow, 1) &
1472 : + quad_ket6*tb%pot%vqp(6, irow, 1) &
1473 : + quad_bra1*tb%pot%vqp(1, icol, 1) &
1474 : + quad_bra2*tb%pot%vqp(2, icol, 1) &
1475 : + quad_bra3*tb%pot%vqp(3, icol, 1) &
1476 : + quad_bra4*tb%pot%vqp(4, icol, 1) &
1477 : + quad_bra5*tb%pot%vqp(5, icol, 1) &
1478 0 : + quad_bra6*tb%pot%vqp(6, icol, 1))
1479 : END DO
1480 : END DO
1481 0 : CALL neighbor_list_iterator_release(nl_iterator)
1482 : END IF
1483 :
1484 : END IF
1485 :
1486 : #else
1487 : MARK_USED(qs_env)
1488 : MARK_USED(tb)
1489 : MARK_USED(dft_control)
1490 : CPABORT("Built without TBLITE")
1491 : #endif
1492 :
1493 24372 : END SUBROUTINE tb_ham_add_coulomb
1494 :
1495 : ! **************************************************************************************************
1496 : !> \brief ...
1497 : !> \param qs_env ...
1498 : !> \param tb ...
1499 : ! **************************************************************************************************
1500 234 : SUBROUTINE tb_get_multipole(qs_env, tb)
1501 :
1502 : TYPE(qs_environment_type), POINTER :: qs_env
1503 : TYPE(tblite_type), POINTER :: tb
1504 :
1505 : #if defined(__TBLITE)
1506 :
1507 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tb_get_multipole'
1508 :
1509 : INTEGER :: ikind, jkind, iatom, jatom, icol, irow, iset, jset, ityp, jtyp
1510 : INTEGER :: nkind, natom, handle, nimg, i, inda, indb
1511 : INTEGER :: atom_a, atom_b, nseta, nsetb, ia, ib, ij
1512 : LOGICAL :: found
1513 : REAL(KIND=dp) :: r2
1514 : REAL(KIND=dp), DIMENSION(3) :: rij
1515 234 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max
1516 234 : INTEGER, DIMENSION(:), POINTER :: nsgfa, nsgfb
1517 234 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1518 234 : INTEGER, ALLOCATABLE :: atom_of_kind(:)
1519 234 : REAL(KIND=dp), ALLOCATABLE :: stmp(:)
1520 234 : REAL(KIND=dp), ALLOCATABLE :: dtmp(:, :), qtmp(:, :), dtmpj(:, :), qtmpj(:, :)
1521 234 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dip_ket1, dip_ket2, dip_ket3, &
1522 234 : dip_bra1, dip_bra2, dip_bra3
1523 234 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: quad_ket1, quad_ket2, quad_ket3, &
1524 234 : quad_ket4, quad_ket5, quad_ket6, &
1525 234 : quad_bra1, quad_bra2, quad_bra3, &
1526 234 : quad_bra4, quad_bra5, quad_bra6
1527 :
1528 234 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1529 234 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1530 : TYPE(dft_control_type), POINTER :: dft_control
1531 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1532 234 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1533 234 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
1534 : TYPE(neighbor_list_iterator_p_type), &
1535 234 : DIMENSION(:), POINTER :: nl_iterator
1536 234 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1537 234 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1538 :
1539 234 : CALL timeset(routineN, handle)
1540 :
1541 : !get info from environment vaiarable
1542 234 : NULLIFY (atomic_kind_set, qs_kind_set, sab_orb, particle_set)
1543 234 : NULLIFY (dft_control, matrix_s)
1544 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
1545 : qs_kind_set=qs_kind_set, &
1546 : sab_orb=sab_orb, &
1547 : particle_set=particle_set, &
1548 : dft_control=dft_control, &
1549 234 : matrix_s_kp=matrix_s)
1550 234 : natom = SIZE(particle_set)
1551 234 : nkind = SIZE(atomic_kind_set)
1552 234 : nimg = dft_control%nimages
1553 :
1554 234 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
1555 :
1556 : !set up basis set lists
1557 1352 : ALLOCATE (basis_set_list(nkind))
1558 234 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
1559 :
1560 702 : ALLOCATE (stmp(msao(tb%calc%bas%maxl)**2))
1561 702 : ALLOCATE (dtmp(dip_n, msao(tb%calc%bas%maxl)**2))
1562 702 : ALLOCATE (qtmp(quad_n, msao(tb%calc%bas%maxl)**2))
1563 468 : ALLOCATE (dtmpj(dip_n, msao(tb%calc%bas%maxl)**2))
1564 468 : ALLOCATE (qtmpj(quad_n, msao(tb%calc%bas%maxl)**2))
1565 :
1566 : ! allocate dipole/quadrupole moment matrix elemnts
1567 234 : CALL dbcsr_allocate_matrix_set(tb%dipbra, dip_n)
1568 234 : CALL dbcsr_allocate_matrix_set(tb%dipket, dip_n)
1569 234 : CALL dbcsr_allocate_matrix_set(tb%quadbra, quad_n)
1570 234 : CALL dbcsr_allocate_matrix_set(tb%quadket, quad_n)
1571 936 : DO i = 1, dip_n
1572 702 : ALLOCATE (tb%dipbra(i)%matrix)
1573 702 : ALLOCATE (tb%dipket(i)%matrix)
1574 : CALL dbcsr_create(tb%dipbra(i)%matrix, template=matrix_s(1, 1)%matrix, &
1575 702 : name="DIPOLE BRAMATRIX")
1576 : CALL dbcsr_create(tb%dipket(i)%matrix, template=matrix_s(1, 1)%matrix, &
1577 702 : name="DIPOLE KETMATRIX")
1578 702 : CALL cp_dbcsr_alloc_block_from_nbl(tb%dipbra(i)%matrix, sab_orb)
1579 936 : CALL cp_dbcsr_alloc_block_from_nbl(tb%dipket(i)%matrix, sab_orb)
1580 : END DO
1581 1638 : DO i = 1, quad_n
1582 1404 : ALLOCATE (tb%quadbra(i)%matrix)
1583 1404 : ALLOCATE (tb%quadket(i)%matrix)
1584 : CALL dbcsr_create(tb%quadbra(i)%matrix, template=matrix_s(1, 1)%matrix, &
1585 1404 : name="QUADRUPOLE BRAMATRIX")
1586 : CALL dbcsr_create(tb%quadket(i)%matrix, template=matrix_s(1, 1)%matrix, &
1587 1404 : name="QUADRUPOLE KETMATRIX")
1588 1404 : CALL cp_dbcsr_alloc_block_from_nbl(tb%quadbra(i)%matrix, sab_orb)
1589 1638 : CALL cp_dbcsr_alloc_block_from_nbl(tb%quadket(i)%matrix, sab_orb)
1590 : END DO
1591 :
1592 : !loop over all atom pairs with a non-zero overlap (sab_orb)
1593 234 : NULLIFY (nl_iterator)
1594 234 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1595 1254 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1596 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
1597 1020 : iatom=iatom, jatom=jatom, r=rij)
1598 :
1599 4080 : r2 = NORM2(rij(:))**2
1600 :
1601 1020 : icol = MAX(iatom, jatom)
1602 1020 : irow = MIN(iatom, jatom)
1603 :
1604 1020 : IF (icol == jatom) THEN
1605 2512 : rij = -rij
1606 628 : i = ikind
1607 628 : ikind = jkind
1608 628 : jkind = i
1609 : END IF
1610 :
1611 1020 : ityp = tb%mol%id(icol)
1612 1020 : jtyp = tb%mol%id(irow)
1613 :
1614 1020 : NULLIFY (dip_bra1, dip_bra2, dip_bra3)
1615 : CALL dbcsr_get_block_p(matrix=tb%dipbra(1)%matrix, &
1616 1020 : row=irow, col=icol, BLOCK=dip_bra1, found=found)
1617 1020 : CPASSERT(found)
1618 : CALL dbcsr_get_block_p(matrix=tb%dipbra(2)%matrix, &
1619 1020 : row=irow, col=icol, BLOCK=dip_bra2, found=found)
1620 1020 : CPASSERT(found)
1621 : CALL dbcsr_get_block_p(matrix=tb%dipbra(3)%matrix, &
1622 1020 : row=irow, col=icol, BLOCK=dip_bra3, found=found)
1623 1020 : CPASSERT(found)
1624 :
1625 1020 : NULLIFY (dip_ket1, dip_ket2, dip_ket3)
1626 : CALL dbcsr_get_block_p(matrix=tb%dipket(1)%matrix, &
1627 1020 : row=irow, col=icol, BLOCK=dip_ket1, found=found)
1628 1020 : CPASSERT(found)
1629 : CALL dbcsr_get_block_p(matrix=tb%dipket(2)%matrix, &
1630 1020 : row=irow, col=icol, BLOCK=dip_ket2, found=found)
1631 1020 : CPASSERT(found)
1632 : CALL dbcsr_get_block_p(matrix=tb%dipket(3)%matrix, &
1633 1020 : row=irow, col=icol, BLOCK=dip_ket3, found=found)
1634 1020 : CPASSERT(found)
1635 :
1636 1020 : NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
1637 : CALL dbcsr_get_block_p(matrix=tb%quadbra(1)%matrix, &
1638 1020 : row=irow, col=icol, BLOCK=quad_bra1, found=found)
1639 1020 : CPASSERT(found)
1640 : CALL dbcsr_get_block_p(matrix=tb%quadbra(2)%matrix, &
1641 1020 : row=irow, col=icol, BLOCK=quad_bra2, found=found)
1642 1020 : CPASSERT(found)
1643 : CALL dbcsr_get_block_p(matrix=tb%quadbra(3)%matrix, &
1644 1020 : row=irow, col=icol, BLOCK=quad_bra3, found=found)
1645 1020 : CPASSERT(found)
1646 : CALL dbcsr_get_block_p(matrix=tb%quadbra(4)%matrix, &
1647 1020 : row=irow, col=icol, BLOCK=quad_bra4, found=found)
1648 1020 : CPASSERT(found)
1649 : CALL dbcsr_get_block_p(matrix=tb%quadbra(5)%matrix, &
1650 1020 : row=irow, col=icol, BLOCK=quad_bra5, found=found)
1651 1020 : CPASSERT(found)
1652 : CALL dbcsr_get_block_p(matrix=tb%quadbra(6)%matrix, &
1653 1020 : row=irow, col=icol, BLOCK=quad_bra6, found=found)
1654 1020 : CPASSERT(found)
1655 :
1656 1020 : NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
1657 : CALL dbcsr_get_block_p(matrix=tb%quadket(1)%matrix, &
1658 1020 : row=irow, col=icol, BLOCK=quad_ket1, found=found)
1659 1020 : CPASSERT(found)
1660 : CALL dbcsr_get_block_p(matrix=tb%quadket(2)%matrix, &
1661 1020 : row=irow, col=icol, BLOCK=quad_ket2, found=found)
1662 1020 : CPASSERT(found)
1663 : CALL dbcsr_get_block_p(matrix=tb%quadket(3)%matrix, &
1664 1020 : row=irow, col=icol, BLOCK=quad_ket3, found=found)
1665 1020 : CPASSERT(found)
1666 : CALL dbcsr_get_block_p(matrix=tb%quadket(4)%matrix, &
1667 1020 : row=irow, col=icol, BLOCK=quad_ket4, found=found)
1668 1020 : CPASSERT(found)
1669 : CALL dbcsr_get_block_p(matrix=tb%quadket(5)%matrix, &
1670 1020 : row=irow, col=icol, BLOCK=quad_ket5, found=found)
1671 1020 : CPASSERT(found)
1672 : CALL dbcsr_get_block_p(matrix=tb%quadket(6)%matrix, &
1673 1020 : row=irow, col=icol, BLOCK=quad_ket6, found=found)
1674 1020 : CPASSERT(found)
1675 :
1676 : !get basis information
1677 1020 : basis_set_a => basis_set_list(ikind)%gto_basis_set
1678 1020 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
1679 1020 : basis_set_b => basis_set_list(jkind)%gto_basis_set
1680 1020 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
1681 1020 : atom_a = atom_of_kind(icol)
1682 1020 : atom_b = atom_of_kind(irow)
1683 : ! basis a
1684 1020 : first_sgfa => basis_set_a%first_sgf
1685 1020 : la_max => basis_set_a%lmax
1686 1020 : nseta = basis_set_a%nset
1687 1020 : nsgfa => basis_set_a%nsgf_set
1688 : ! basis b
1689 1020 : first_sgfb => basis_set_b%first_sgf
1690 1020 : lb_max => basis_set_b%lmax
1691 1020 : nsetb = basis_set_b%nset
1692 1020 : nsgfb => basis_set_b%nsgf_set
1693 :
1694 : ! --------- Hamiltonian
1695 1254 : IF (icol == irow) THEN
1696 1242 : DO iset = 1, nseta
1697 2929 : DO jset = 1, nsetb
1698 : CALL multipole_cgto(tb%calc%bas%cgto(jset, ityp), tb%calc%bas%cgto(iset, ityp), &
1699 6748 : & r2, -rij, tb%calc%bas%intcut, stmp, dtmp, qtmp)
1700 :
1701 6253 : DO inda = 1, nsgfa(iset)
1702 3753 : ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
1703 14279 : DO indb = 1, nsgfb(jset)
1704 8839 : ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
1705 8839 : ij = indb + nsgfb(jset)*(inda - 1)
1706 :
1707 8839 : dip_ket1(ib, ia) = dtmp(1, ij)
1708 8839 : dip_ket2(ib, ia) = dtmp(2, ij)
1709 8839 : dip_ket3(ib, ia) = dtmp(3, ij)
1710 :
1711 8839 : quad_ket1(ib, ia) = qtmp(1, ij)
1712 8839 : quad_ket2(ib, ia) = qtmp(2, ij)
1713 8839 : quad_ket3(ib, ia) = qtmp(3, ij)
1714 8839 : quad_ket4(ib, ia) = qtmp(4, ij)
1715 8839 : quad_ket5(ib, ia) = qtmp(5, ij)
1716 8839 : quad_ket6(ib, ia) = qtmp(6, ij)
1717 :
1718 8839 : dip_bra1(ib, ia) = dtmp(1, ij)
1719 8839 : dip_bra2(ib, ia) = dtmp(2, ij)
1720 8839 : dip_bra3(ib, ia) = dtmp(3, ij)
1721 :
1722 8839 : quad_bra1(ib, ia) = qtmp(1, ij)
1723 8839 : quad_bra2(ib, ia) = qtmp(2, ij)
1724 8839 : quad_bra3(ib, ia) = qtmp(3, ij)
1725 8839 : quad_bra4(ib, ia) = qtmp(4, ij)
1726 8839 : quad_bra5(ib, ia) = qtmp(5, ij)
1727 12592 : quad_bra6(ib, ia) = qtmp(6, ij)
1728 : END DO
1729 : END DO
1730 : END DO
1731 : END DO
1732 : ELSE
1733 1652 : DO iset = 1, nseta
1734 3813 : DO jset = 1, nsetb
1735 : CALL multipole_cgto(tb%calc%bas%cgto(jset, jtyp), tb%calc%bas%cgto(iset, ityp), &
1736 8644 : & r2, -rij, tb%calc%bas%intcut, stmp, dtmp, qtmp)
1737 : CALL multipole_cgto(tb%calc%bas%cgto(iset, ityp), tb%calc%bas%cgto(jset, jtyp), &
1738 2161 : & r2, rij, tb%calc%bas%intcut, stmp, dtmpj, qtmpj)
1739 :
1740 7947 : DO inda = 1, nsgfa(iset)
1741 4725 : ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
1742 17043 : DO indb = 1, nsgfb(jset)
1743 10157 : ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
1744 :
1745 10157 : ij = indb + nsgfb(jset)*(inda - 1)
1746 :
1747 10157 : dip_bra1(ib, ia) = dtmp(1, ij)
1748 10157 : dip_bra2(ib, ia) = dtmp(2, ij)
1749 10157 : dip_bra3(ib, ia) = dtmp(3, ij)
1750 :
1751 10157 : quad_bra1(ib, ia) = qtmp(1, ij)
1752 10157 : quad_bra2(ib, ia) = qtmp(2, ij)
1753 10157 : quad_bra3(ib, ia) = qtmp(3, ij)
1754 10157 : quad_bra4(ib, ia) = qtmp(4, ij)
1755 10157 : quad_bra5(ib, ia) = qtmp(5, ij)
1756 10157 : quad_bra6(ib, ia) = qtmp(6, ij)
1757 :
1758 10157 : ij = inda + nsgfa(iset)*(indb - 1)
1759 :
1760 10157 : dip_ket1(ib, ia) = dtmpj(1, ij)
1761 10157 : dip_ket2(ib, ia) = dtmpj(2, ij)
1762 10157 : dip_ket3(ib, ia) = dtmpj(3, ij)
1763 :
1764 10157 : quad_ket1(ib, ia) = qtmpj(1, ij)
1765 10157 : quad_ket2(ib, ia) = qtmpj(2, ij)
1766 10157 : quad_ket3(ib, ia) = qtmpj(3, ij)
1767 10157 : quad_ket4(ib, ia) = qtmpj(4, ij)
1768 10157 : quad_ket5(ib, ia) = qtmpj(5, ij)
1769 14882 : quad_ket6(ib, ia) = qtmpj(6, ij)
1770 : END DO
1771 : END DO
1772 : END DO
1773 : END DO
1774 : END IF
1775 : END DO
1776 234 : CALL neighbor_list_iterator_release(nl_iterator)
1777 :
1778 936 : DO i = 1, dip_n
1779 702 : CALL dbcsr_finalize(tb%dipbra(i)%matrix)
1780 936 : CALL dbcsr_finalize(tb%dipket(i)%matrix)
1781 : END DO
1782 1638 : DO i = 1, quad_n
1783 1404 : CALL dbcsr_finalize(tb%quadbra(i)%matrix)
1784 1638 : CALL dbcsr_finalize(tb%quadket(i)%matrix)
1785 : END DO
1786 :
1787 234 : DEALLOCATE (basis_set_list)
1788 :
1789 234 : CALL timestop(handle)
1790 :
1791 : #else
1792 : MARK_USED(qs_env)
1793 : MARK_USED(tb)
1794 : CPABORT("Built without TBLITE")
1795 : #endif
1796 :
1797 468 : END SUBROUTINE tb_get_multipole
1798 :
1799 : ! **************************************************************************************************
1800 : !> \brief compute the mulliken properties (AO resolved)
1801 : !> \param p_matrix ...
1802 : !> \param bra_mat ...
1803 : !> \param ket_mat ...
1804 : !> \param output ...
1805 : !> \param para_env ...
1806 : !> \par History
1807 : !> adapted from ao_charges_2
1808 : !> \note
1809 : ! **************************************************************************************************
1810 42300 : SUBROUTINE contract_dens(p_matrix, bra_mat, ket_mat, output, para_env)
1811 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
1812 : TYPE(dbcsr_type), POINTER :: bra_mat, ket_mat
1813 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: output
1814 : TYPE(mp_para_env_type), POINTER :: para_env
1815 :
1816 : CHARACTER(len=*), PARAMETER :: routineN = 'contract_dens'
1817 :
1818 : INTEGER :: handle, i, iblock_col, iblock_row, &
1819 : ispin, j, nspin
1820 : LOGICAL :: found
1821 42300 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: bra, ket, p_block
1822 : TYPE(dbcsr_iterator_type) :: iter
1823 :
1824 42300 : CALL timeset(routineN, handle)
1825 :
1826 42300 : nspin = SIZE(p_matrix)
1827 191862 : output = 0.0_dp
1828 84600 : DO ispin = 1, nspin
1829 42300 : CALL dbcsr_iterator_start(iter, bra_mat)
1830 216090 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1831 173790 : NULLIFY (p_block, bra, ket)
1832 :
1833 173790 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, bra)
1834 : CALL dbcsr_get_block_p(matrix=p_matrix(ispin)%matrix, &
1835 173790 : row=iblock_row, col=iblock_col, BLOCK=p_block, found=found)
1836 173790 : IF (.NOT. found) CYCLE
1837 : CALL dbcsr_get_block_p(matrix=ket_mat, &
1838 173790 : row=iblock_row, col=iblock_col, BLOCK=ket, found=found)
1839 173790 : IF (.NOT. found) CPABORT("missing block")
1840 :
1841 173790 : IF (.NOT. (ASSOCIATED(bra) .AND. ASSOCIATED(p_block))) CYCLE
1842 880353 : DO j = 1, SIZE(p_block, 1)
1843 4192092 : DO i = 1, SIZE(p_block, 2)
1844 4018302 : output(iblock_row) = output(iblock_row) + p_block(j, i)*ket(j, i)
1845 : END DO
1846 : END DO
1847 216090 : IF (iblock_col /= iblock_row) THEN
1848 503631 : DO j = 1, SIZE(p_block, 1)
1849 2270349 : DO i = 1, SIZE(p_block, 2)
1850 2171340 : output(iblock_col) = output(iblock_col) + p_block(j, i)*bra(j, i)
1851 : END DO
1852 : END DO
1853 : END IF
1854 : END DO
1855 126900 : CALL dbcsr_iterator_stop(iter)
1856 : END DO
1857 :
1858 341424 : CALL para_env%sum(output)
1859 42300 : CALL timestop(handle)
1860 :
1861 42300 : END SUBROUTINE contract_dens
1862 :
1863 : ! **************************************************************************************************
1864 : !> \brief save gradient to force
1865 : !> \param qs_env ...
1866 : !> \param tb ...
1867 : !> \param para_env ...
1868 : !> \param ityp ...
1869 : !> \note
1870 : ! **************************************************************************************************
1871 2562 : SUBROUTINE tb_grad2force(qs_env, tb, para_env, ityp)
1872 :
1873 : TYPE(qs_environment_type) :: qs_env
1874 : TYPE(tblite_type) :: tb
1875 : TYPE(mp_para_env_type) :: para_env
1876 : INTEGER :: ityp
1877 :
1878 : CHARACTER(len=*), PARAMETER :: routineN = 'tb_grad2force'
1879 :
1880 : INTEGER :: atoma, handle, iatom, ikind, natom
1881 2562 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
1882 2562 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1883 2562 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1884 2562 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1885 :
1886 2562 : CALL timeset(routineN, handle)
1887 :
1888 2562 : NULLIFY (force, atomic_kind_set)
1889 : CALL get_qs_env(qs_env=qs_env, force=force, particle_set=particle_set, &
1890 2562 : atomic_kind_set=atomic_kind_set)
1891 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1892 2562 : atom_of_kind=atom_of_kind, kind_of=kind_of)
1893 :
1894 2562 : natom = SIZE(particle_set)
1895 :
1896 2562 : SELECT CASE (ityp)
1897 : CASE DEFAULT
1898 0 : CPABORT("unknown force type")
1899 : CASE (0)
1900 3868 : DO iatom = 1, natom
1901 3210 : ikind = kind_of(iatom)
1902 3210 : atoma = atom_of_kind(iatom)
1903 : force(ikind)%all_potential(:, atoma) = &
1904 13498 : force(ikind)%all_potential(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
1905 : END DO
1906 : CASE (1)
1907 4898 : DO iatom = 1, natom
1908 4018 : ikind = kind_of(iatom)
1909 4018 : atoma = atom_of_kind(iatom)
1910 : force(ikind)%repulsive(:, atoma) = &
1911 16952 : force(ikind)%repulsive(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
1912 : END DO
1913 : CASE (2)
1914 5084 : DO iatom = 1, natom
1915 4168 : ikind = kind_of(iatom)
1916 4168 : atoma = atom_of_kind(iatom)
1917 : force(ikind)%dispersion(:, atoma) = &
1918 17588 : force(ikind)%dispersion(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
1919 : END DO
1920 : CASE (3)
1921 186 : DO iatom = 1, natom
1922 150 : ikind = kind_of(iatom)
1923 150 : atoma = atom_of_kind(iatom)
1924 : force(ikind)%rho_elec(:, atoma) = &
1925 636 : force(ikind)%rho_elec(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
1926 : END DO
1927 : CASE (4)
1928 372 : DO iatom = 1, natom
1929 300 : ikind = kind_of(iatom)
1930 300 : atoma = atom_of_kind(iatom)
1931 : force(ikind)%overlap(:, atoma) = &
1932 1272 : force(ikind)%overlap(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
1933 : END DO
1934 : CASE (5)
1935 2562 : DO iatom = 1, natom
1936 0 : ikind = kind_of(iatom)
1937 0 : atoma = atom_of_kind(iatom)
1938 : force(ikind)%efield(:, atoma) = &
1939 0 : force(ikind)%efield(:, atoma) + tb%grad(:, iatom)/para_env%num_pe
1940 : END DO
1941 : END SELECT
1942 :
1943 2562 : CALL timestop(handle)
1944 :
1945 5124 : END SUBROUTINE tb_grad2force
1946 :
1947 : ! **************************************************************************************************
1948 : !> \brief set gradient to zero
1949 : !> \param qs_env ...
1950 : !> \note
1951 : ! **************************************************************************************************
1952 880 : SUBROUTINE tb_zero_force(qs_env)
1953 :
1954 : TYPE(qs_environment_type) :: qs_env
1955 :
1956 : INTEGER :: iatom, ikind, natom
1957 880 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
1958 880 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1959 880 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1960 880 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1961 :
1962 880 : NULLIFY (force, atomic_kind_set)
1963 : CALL get_qs_env(qs_env=qs_env, force=force, particle_set=particle_set, &
1964 880 : atomic_kind_set=atomic_kind_set)
1965 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1966 880 : kind_of=kind_of)
1967 :
1968 880 : natom = SIZE(particle_set)
1969 :
1970 4898 : DO iatom = 1, natom
1971 4018 : ikind = kind_of(iatom)
1972 53354 : force(ikind)%all_potential = 0.0_dp
1973 53354 : force(ikind)%repulsive = 0.0_dp
1974 53354 : force(ikind)%dispersion = 0.0_dp
1975 53354 : force(ikind)%rho_elec = 0.0_dp
1976 53354 : force(ikind)%overlap = 0.0_dp
1977 54234 : force(ikind)%efield = 0.0_dp
1978 : END DO
1979 :
1980 1760 : END SUBROUTINE tb_zero_force
1981 :
1982 : ! **************************************************************************************************
1983 : !> \brief compute the derivative of the energy
1984 : !> \param qs_env ...
1985 : !> \param use_rho ...
1986 : !> \param nimg ...
1987 : ! **************************************************************************************************
1988 36 : SUBROUTINE tb_derive_dH_diag(qs_env, use_rho, nimg)
1989 :
1990 : TYPE(qs_environment_type), POINTER :: qs_env
1991 : LOGICAL, INTENT(IN) :: use_rho
1992 : INTEGER, INTENT(IN) :: nimg
1993 :
1994 : #if defined(__TBLITE)
1995 : INTEGER :: i, iatom, ic, ikind, ind, is, ish, &
1996 : jatom, jkind
1997 : INTEGER, DIMENSION(3) :: cellind
1998 36 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1999 : LOGICAL :: found
2000 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dE
2001 36 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock
2002 36 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
2003 : TYPE(kpoint_type), POINTER :: kpoints
2004 : TYPE(mp_para_env_type), POINTER :: para_env
2005 : TYPE(neighbor_list_iterator_p_type), &
2006 36 : DIMENSION(:), POINTER :: nl_iterator
2007 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2008 36 : POINTER :: sab_orb
2009 : TYPE(qs_rho_type), POINTER :: rho
2010 : TYPE(qs_scf_env_type), POINTER :: scf_env
2011 : TYPE(tblite_type), POINTER :: tb
2012 :
2013 : ! compute mulliken charges required for charge update
2014 36 : NULLIFY (scf_env, rho, tb, sab_orb, para_env, kpoints)
2015 : CALL get_qs_env(qs_env=qs_env, scf_env=scf_env, rho=rho, tb_tblite=tb, &
2016 36 : sab_orb=sab_orb, para_env=para_env, kpoints=kpoints)
2017 :
2018 36 : NULLIFY (cell_to_index)
2019 36 : IF (nimg > 1) THEN
2020 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
2021 : END IF
2022 :
2023 36 : NULLIFY (matrix_p)
2024 36 : IF (use_rho) THEN
2025 36 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2026 : ELSE
2027 0 : matrix_p => scf_env%p_mix_new
2028 : END IF
2029 :
2030 108 : ALLOCATE (dE(tb%mol%nat))
2031 :
2032 186 : dE = 0.0_dp
2033 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
2034 36 : NULLIFY (nl_iterator)
2035 36 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
2036 3718 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2037 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2038 3682 : iatom=iatom, jatom=jatom, cell=cellind)
2039 :
2040 3682 : IF (iatom /= jatom) CYCLE
2041 :
2042 1045 : IF (ikind /= jkind) CPABORT("Type wrong")
2043 :
2044 1045 : is = tb%calc%bas%ish_at(iatom)
2045 :
2046 1045 : IF (nimg == 1) THEN
2047 : ic = 1
2048 : ELSE
2049 0 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2050 0 : CPASSERT(ic > 0)
2051 : END IF
2052 :
2053 1045 : IF (cellind(1) /= 0) CYCLE
2054 327 : IF (cellind(2) /= 0) CYCLE
2055 129 : IF (cellind(3) /= 0) CYCLE
2056 :
2057 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
2058 75 : row=iatom, col=jatom, BLOCK=pblock, found=found)
2059 :
2060 75 : IF (.NOT. found) CPABORT("block not found")
2061 :
2062 75 : ind = 0
2063 281 : DO ish = 1, tb%calc%bas%nsh_id(ikind)
2064 4208 : DO i = 1, msao(tb%calc%bas%cgto(ish, ikind)%ang)
2065 356 : ind = ind + 1
2066 526 : dE(iatom) = dE(iatom) + tb%dsedcn(is + ish)*pblock(ind, ind)
2067 : END DO
2068 : END DO
2069 : END DO
2070 36 : CALL neighbor_list_iterator_release(nl_iterator)
2071 36 : CALL para_env%sum(dE)
2072 :
2073 636 : tb%grad = 0.0_dp
2074 36 : CALL tb_add_grad(tb%grad, tb%dcndr, dE, tb%mol%nat)
2075 36 : IF (tb%use_virial) CALL tb_add_sig(tb%sigma, tb%dcndL, dE, tb%mol%nat)
2076 36 : CALL tb_grad2force(qs_env, tb, para_env, 4)
2077 36 : DEALLOCATE (dE)
2078 :
2079 : #else
2080 : MARK_USED(qs_env)
2081 : MARK_USED(use_rho)
2082 : MARK_USED(nimg)
2083 : CPABORT("Built without TBLITE")
2084 : #endif
2085 :
2086 36 : END SUBROUTINE tb_derive_dH_diag
2087 :
2088 : ! **************************************************************************************************
2089 : !> \brief compute the derivative of the energy
2090 : !> \param qs_env ...
2091 : !> \param use_rho ...
2092 : !> \param nimg ...
2093 : ! **************************************************************************************************
2094 36 : SUBROUTINE tb_derive_dH_off(qs_env, use_rho, nimg)
2095 :
2096 : TYPE(qs_environment_type), POINTER :: qs_env
2097 : LOGICAL, INTENT(IN) :: use_rho
2098 : INTEGER, INTENT(IN) :: nimg
2099 :
2100 : #if defined(__TBLITE)
2101 : INTEGER :: i, ij, iatom, ic, icol, ikind, ni, nj, nkind, nel, &
2102 : sgfi, sgfj, ityp, jatom, jkind, jrow, jtyp, iset, jset, &
2103 : nseti, nsetj, ia, ib, inda, indb
2104 : INTEGER, DIMENSION(3) :: cellind
2105 36 : INTEGER, DIMENSION(:), POINTER :: nsgfa, nsgfb
2106 36 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
2107 36 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2108 : LOGICAL :: found
2109 : REAL(KIND=dp) :: r2, dr, itemp, jtemp, rr, hij, shpoly, dshpoly, idHdc, jdHdc, &
2110 : scal, hp, i_a_shift, j_a_shift, ishift, jshift
2111 : REAL(KIND=dp), DIMENSION(3) :: rij, dgrad
2112 : REAL(KIND=dp), DIMENSION(3, 3) :: hsigma
2113 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dE, t_ov, idip, jdip, iquad, jquad
2114 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: t_dip, t_quad, t_d_ov
2115 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: t_i_dip, t_i_quad, t_j_dip, t_j_quad
2116 36 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock, wblock
2117 :
2118 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2119 36 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_w
2120 36 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
2121 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
2122 : TYPE(kpoint_type), POINTER :: kpoints
2123 : TYPE(mp_para_env_type), POINTER :: para_env
2124 : TYPE(neighbor_list_iterator_p_type), &
2125 36 : DIMENSION(:), POINTER :: nl_iterator
2126 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2127 36 : POINTER :: sab_orb
2128 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2129 : TYPE(qs_rho_type), POINTER :: rho
2130 : TYPE(qs_scf_env_type), POINTER :: scf_env
2131 : TYPE(tb_hamiltonian), POINTER :: h0
2132 : TYPE(tblite_type), POINTER :: tb
2133 :
2134 : ! compute mulliken charges required for charge update
2135 36 : NULLIFY (scf_env, rho, tb, sab_orb, para_env, kpoints, matrix_w)
2136 : CALL get_qs_env(qs_env=qs_env, &
2137 : atomic_kind_set=atomic_kind_set, &
2138 : scf_env=scf_env, &
2139 : rho=rho, &
2140 : tb_tblite=tb, &
2141 : sab_orb=sab_orb, &
2142 : para_env=para_env, &
2143 : qs_kind_set=qs_kind_set, &
2144 : kpoints=kpoints, &
2145 36 : matrix_w_kp=matrix_w)
2146 :
2147 36 : NULLIFY (cell_to_index)
2148 36 : IF (nimg > 1) THEN
2149 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
2150 : END IF
2151 :
2152 36 : h0 => tb%calc%h0
2153 :
2154 36 : NULLIFY (matrix_p)
2155 36 : IF (use_rho) THEN
2156 36 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2157 : ELSE
2158 0 : matrix_p => scf_env%p_mix_new
2159 : END IF
2160 :
2161 : ! set up basis set lists
2162 36 : nkind = SIZE(atomic_kind_set)
2163 206 : ALLOCATE (basis_set_list(nkind))
2164 36 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
2165 :
2166 108 : ALLOCATE (dE(tb%mol%nat))
2167 :
2168 36 : nel = msao(tb%calc%bas%maxl)**2
2169 108 : ALLOCATE (t_ov(nel))
2170 108 : ALLOCATE (t_d_ov(3, nel))
2171 72 : ALLOCATE (t_dip(dip_n, nel))
2172 144 : ALLOCATE (t_i_dip(3, dip_n, nel), t_j_dip(3, dip_n, nel))
2173 108 : ALLOCATE (t_quad(quad_n, nel))
2174 144 : ALLOCATE (t_i_quad(3, quad_n, nel), t_j_quad(3, quad_n, nel))
2175 :
2176 36 : ALLOCATE (idip(dip_n), jdip(dip_n))
2177 36 : ALLOCATE (iquad(quad_n), jquad(quad_n))
2178 :
2179 186 : dE = 0.0_dp
2180 636 : tb%grad = 0.0_dp
2181 36 : hsigma = 0.0_dp
2182 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
2183 36 : NULLIFY (nl_iterator)
2184 36 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
2185 3718 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2186 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2187 3682 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
2188 :
2189 3682 : icol = MAX(iatom, jatom)
2190 3682 : jrow = MIN(iatom, jatom)
2191 :
2192 3682 : IF (icol == jatom) THEN
2193 8576 : rij = -rij
2194 2144 : i = ikind
2195 2144 : ikind = jkind
2196 2144 : jkind = i
2197 : END IF
2198 :
2199 3682 : ityp = tb%mol%id(icol)
2200 3682 : jtyp = tb%mol%id(jrow)
2201 :
2202 14728 : r2 = DOT_PRODUCT(rij, rij)
2203 3682 : dr = SQRT(r2)
2204 3682 : IF (icol == jrow .AND. dr < same_atom) CYCLE
2205 3607 : rr = SQRT(dr/(h0%rad(ikind) + h0%rad(jkind)))
2206 :
2207 : !get basis information
2208 3607 : basis_set_a => basis_set_list(ikind)%gto_basis_set
2209 3607 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
2210 3607 : first_sgfa => basis_set_a%first_sgf
2211 3607 : nsgfa => basis_set_a%nsgf_set
2212 3607 : nseti = basis_set_a%nset
2213 3607 : basis_set_b => basis_set_list(jkind)%gto_basis_set
2214 3607 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
2215 3607 : first_sgfb => basis_set_b%first_sgf
2216 3607 : nsgfb => basis_set_b%nsgf_set
2217 3607 : nsetj = basis_set_b%nset
2218 :
2219 3607 : IF (nimg == 1) THEN
2220 : ic = 1
2221 : ELSE
2222 0 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
2223 0 : CPASSERT(ic > 0)
2224 : END IF
2225 :
2226 3607 : NULLIFY (pblock)
2227 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
2228 3607 : row=jrow, col=icol, BLOCK=pblock, found=found)
2229 3607 : IF (.NOT. found) CPABORT("pblock not found")
2230 :
2231 3607 : NULLIFY (wblock)
2232 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
2233 3607 : row=jrow, col=icol, block=wblock, found=found)
2234 3607 : CPASSERT(found)
2235 :
2236 3607 : i_a_shift = tb%pot%vat(icol, 1)
2237 3607 : j_a_shift = tb%pot%vat(jrow, 1)
2238 14428 : idip(:) = tb%pot%vdp(:, icol, 1)
2239 14428 : jdip(:) = tb%pot%vdp(:, jrow, 1)
2240 25249 : iquad(:) = tb%pot%vqp(:, icol, 1)
2241 25249 : jquad(:) = tb%pot%vqp(:, jrow, 1)
2242 :
2243 3607 : ni = tb%calc%bas%ish_at(icol)
2244 13714 : DO iset = 1, nseti
2245 10071 : sgfi = first_sgfa(1, iset)
2246 10071 : ishift = i_a_shift + tb%pot%vsh(ni + iset, 1)
2247 10071 : nj = tb%calc%bas%ish_at(jrow)
2248 42506 : DO jset = 1, nsetj
2249 28753 : sgfj = first_sgfb(1, jset)
2250 28753 : jshift = j_a_shift + tb%pot%vsh(nj + jset, 1)
2251 :
2252 : !get integrals and derivatives
2253 : CALL multipole_grad_cgto(tb%calc%bas%cgto(iset, ityp), tb%calc%bas%cgto(jset, jtyp), &
2254 : & r2, rij, tb%calc%bas%intcut, t_ov, t_dip, t_quad, t_d_ov, t_i_dip, t_i_quad, &
2255 28753 : & t_j_dip, t_j_quad)
2256 :
2257 : shpoly = (1.0_dp + h0%shpoly(iset, ikind)*rr) &
2258 28753 : *(1.0_dp + h0%shpoly(jset, jkind)*rr)
2259 : dshpoly = ((1.0_dp + h0%shpoly(iset, ikind)*rr)*h0%shpoly(jset, jkind)*rr &
2260 : + (1.0_dp + h0%shpoly(jset, jkind)*rr)*h0%shpoly(iset, ikind)*rr) &
2261 28753 : *0.5_dp/r2
2262 28753 : scal = h0%hscale(iset, jset, ikind, jkind)
2263 28753 : hij = 0.5_dp*(tb%selfenergy(ni + iset) + tb%selfenergy(nj + jset))*scal
2264 :
2265 28753 : idHdc = tb%dsedcn(ni + iset)*shpoly*scal
2266 28753 : jdHdc = tb%dsedcn(nj + jset)*shpoly*scal
2267 :
2268 28753 : itemp = 0.0_dp
2269 28753 : jtemp = 0.0_dp
2270 28753 : dgrad = 0.0_dp
2271 110992 : DO inda = 1, nsgfa(iset)
2272 82239 : ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
2273 352097 : DO indb = 1, nsgfb(jset)
2274 241105 : ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
2275 :
2276 241105 : ij = inda + nsgfa(iset)*(indb - 1)
2277 :
2278 241105 : itemp = itemp + idHdc*pblock(ib, ia)*t_ov(ij)
2279 241105 : jtemp = jtemp + jdHdc*pblock(ib, ia)*t_ov(ij)
2280 :
2281 241105 : hp = 2*hij*pblock(ib, ia)
2282 : dgrad(:) = dgrad(:) &
2283 : - (ishift + jshift)*pblock(ib, ia)*t_d_ov(:, ij) &
2284 : - 2*wblock(ib, ia)*t_d_ov(:, ij) &
2285 : + hp*shpoly*t_d_ov(:, ij) &
2286 : + hp*dshpoly*t_ov(ij)*rij &
2287 241105 : - pblock(ib, ia)*( &
2288 : MATMUL(t_i_dip(:, :, ij), idip) &
2289 241105 : + MATMUL(t_j_dip(:, :, ij), jdip) &
2290 241105 : + MATMUL(t_i_quad(:, :, ij), iquad) &
2291 22022794 : + MATMUL(t_j_quad(:, :, ij), jquad))
2292 :
2293 : END DO
2294 : END DO
2295 28753 : dE(icol) = dE(icol) + itemp
2296 28753 : dE(jrow) = dE(jrow) + jtemp
2297 115012 : tb%grad(:, icol) = tb%grad(:, icol) - dgrad
2298 115012 : tb%grad(:, jrow) = tb%grad(:, jrow) + dgrad
2299 38824 : IF (tb%use_virial) THEN
2300 28525 : IF (icol == jrow) THEN
2301 28320 : DO ia = 1, 3
2302 92040 : DO ib = 1, 3
2303 84960 : hsigma(ia, ib) = hsigma(ia, ib) + 0.25_dp*(rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
2304 : END DO
2305 : END DO
2306 : ELSE
2307 85780 : DO ia = 1, 3
2308 278785 : DO ib = 1, 3
2309 257340 : hsigma(ia, ib) = hsigma(ia, ib) + 0.50_dp*(rij(ia)*dgrad(ib) + rij(ib)*dgrad(ia))
2310 : END DO
2311 : END DO
2312 : END IF
2313 : END IF
2314 : END DO
2315 : END DO
2316 : END DO
2317 36 : CALL neighbor_list_iterator_release(nl_iterator)
2318 :
2319 36 : CALL para_env%sum(hsigma)
2320 36 : CALL para_env%sum(dE)
2321 36 : CALL para_env%sum(tb%grad)
2322 :
2323 36 : CALL tb_add_grad(tb%grad, tb%dcndr, dE, tb%mol%nat)
2324 36 : CALL tb_grad2force(qs_env, tb, para_env, 4)
2325 :
2326 468 : tb%sigma = tb%sigma + hsigma
2327 36 : IF (tb%use_virial) CALL tb_add_sig(tb%sigma, tb%dcndL, dE, tb%mol%nat)
2328 :
2329 36 : DEALLOCATE (dE)
2330 36 : DEALLOCATE (basis_set_list)
2331 36 : DEALLOCATE (t_ov, t_d_ov)
2332 36 : DEALLOCATE (t_dip, t_i_dip, t_j_dip)
2333 36 : DEALLOCATE (t_quad, t_i_quad, t_j_quad)
2334 36 : DEALLOCATE (idip, jdip, iquad, jquad)
2335 :
2336 36 : IF (tb%use_virial) CALL tb_add_stress(qs_env, tb, para_env)
2337 :
2338 : #else
2339 : MARK_USED(qs_env)
2340 : MARK_USED(use_rho)
2341 : MARK_USED(nimg)
2342 : CPABORT("Built without TBLITE")
2343 : #endif
2344 :
2345 36 : END SUBROUTINE tb_derive_dH_off
2346 :
2347 : ! **************************************************************************************************
2348 : !> \brief save stress tensor
2349 : !> \param qs_env ...
2350 : !> \param tb ...
2351 : !> \param para_env ...
2352 : ! **************************************************************************************************
2353 22 : SUBROUTINE tb_add_stress(qs_env, tb, para_env)
2354 :
2355 : TYPE(qs_environment_type) :: qs_env
2356 : TYPE(tblite_type) :: tb
2357 : TYPE(mp_para_env_type) :: para_env
2358 :
2359 : TYPE(cell_type), POINTER :: cell
2360 : TYPE(virial_type), POINTER :: virial
2361 :
2362 22 : NULLIFY (virial, cell)
2363 22 : CALL get_qs_env(qs_env=qs_env, virial=virial, cell=cell)
2364 :
2365 286 : virial%pv_virial = virial%pv_virial - tb%sigma/para_env%num_pe
2366 :
2367 22 : END SUBROUTINE tb_add_stress
2368 :
2369 : ! **************************************************************************************************
2370 : !> \brief add contrib. to gradient
2371 : !> \param grad ...
2372 : !> \param deriv ...
2373 : !> \param dE ...
2374 : !> \param natom ...
2375 : ! **************************************************************************************************
2376 72 : SUBROUTINE tb_add_grad(grad, deriv, dE, natom)
2377 :
2378 : REAL(KIND=dp), DIMENSION(:, :) :: grad
2379 : REAL(KIND=dp), DIMENSION(:, :, :) :: deriv
2380 : REAL(KIND=dp), DIMENSION(:) :: dE
2381 : INTEGER :: natom
2382 :
2383 : INTEGER :: i, j
2384 :
2385 372 : DO i = 1, natom
2386 1784 : DO j = 1, natom
2387 5948 : grad(:, i) = grad(:, i) + deriv(:, i, j)*dE(j)
2388 : END DO
2389 : END DO
2390 :
2391 72 : END SUBROUTINE tb_add_grad
2392 :
2393 : ! **************************************************************************************************
2394 : !> \brief add contrib. to sigma
2395 : !> \param sig ...
2396 : !> \param deriv ...
2397 : !> \param dE ...
2398 : !> \param natom ...
2399 : ! **************************************************************************************************
2400 44 : SUBROUTINE tb_add_sig(sig, deriv, dE, natom)
2401 :
2402 : REAL(KIND=dp), DIMENSION(:, :) :: sig
2403 : REAL(KIND=dp), DIMENSION(:, :, :) :: deriv
2404 : REAL(KIND=dp), DIMENSION(:) :: dE
2405 : INTEGER :: natom
2406 :
2407 : INTEGER :: i, j
2408 :
2409 176 : DO i = 1, 3
2410 740 : DO j = 1, natom
2411 2388 : sig(:, i) = sig(:, i) + deriv(:, i, j)*dE(j)
2412 : END DO
2413 : END DO
2414 :
2415 44 : END SUBROUTINE tb_add_sig
2416 :
2417 322491 : END MODULE tblite_interface
2418 :
|
