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