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