(git:e5b1968)
Loading...
Searching...
No Matches
tblite_interface.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief interface to tblite
10!> \author JVP
11!> \history creation 09.2024
12! **************************************************************************************************
13
14MODULE tblite_interface
15
16#if defined(__TBLITE)
17 USE mctc_env, ONLY: error_type
18 USE mctc_io, ONLY: structure_type, new
19 USE mctc_io_symbols, ONLY: symbol_to_number
20 USE tblite_basis_type, ONLY: get_cutoff
21 USE tblite_container, ONLY: container_cache
22 USE tblite_integral_multipole, ONLY: multipole_cgto, multipole_grad_cgto, maxl, msao
23 USE tblite_scf_info, ONLY: scf_info, atom_resolved, shell_resolved, &
24 orbital_resolved, not_used
25 USE tblite_scf_potential, ONLY: potential_type, new_potential
26 USE tblite_wavefunction_type, ONLY: wavefunction_type, new_wavefunction
27 USE tblite_xtb_calculator, ONLY: xtb_calculator, new_xtb_calculator
28 USE tblite_xtb_gfn1, ONLY: new_gfn1_calculator
29 USE tblite_xtb_gfn2, ONLY: new_gfn2_calculator
30 USE tblite_xtb_h0, ONLY: get_selfenergy, get_hamiltonian, get_occupation, &
31 get_hamiltonian_gradient, tb_hamiltonian
32 USE tblite_xtb_ipea1, ONLY: new_ipea1_calculator
33#endif
34 USE ai_contraction, ONLY: block_add, &
35 contraction
36 USE ai_overlap, ONLY: overlap_ab
37 USE atomic_kind_types, ONLY: atomic_kind_type, get_atomic_kind, get_atomic_kind_set
38 USE atprop_types, ONLY: atprop_type
39 USE basis_set_types, ONLY: gto_basis_set_type, gto_basis_set_p_type, &
40 & allocate_gto_basis_set, write_gto_basis_set, process_gto_basis
41 USE cell_types, ONLY: cell_type, get_cell
42 USE cp_blacs_env, ONLY: cp_blacs_env_type
43 USE cp_control_types, ONLY: dft_control_type
44 USE cp_dbcsr_api, ONLY: dbcsr_type, dbcsr_p_type, dbcsr_create, dbcsr_add, &
45 dbcsr_get_block_p, dbcsr_finalize, &
46 dbcsr_iterator_type, dbcsr_iterator_blocks_left, &
47 dbcsr_iterator_start, dbcsr_iterator_stop, &
48 dbcsr_iterator_next_block
49 USE cp_dbcsr_contrib, ONLY: dbcsr_print
50 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
51 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
52 USE cp_log_handling, ONLY: cp_get_default_logger, &
53 cp_logger_type, cp_logger_get_default_io_unit
54 USE cp_output_handling, ONLY: cp_print_key_should_output, &
55 cp_print_key_unit_nr, cp_print_key_finished_output
56 USE input_constants, ONLY: gfn1xtb, gfn2xtb, ipea1xtb
57 USE input_section_types, ONLY: section_vals_val_get
58 USE kinds, ONLY: dp, default_string_length
59 USE kpoint_types, ONLY: get_kpoint_info, kpoint_type
60 USE memory_utilities, ONLY: reallocate
61 USE message_passing, ONLY: mp_para_env_type
62 USE mulliken, ONLY: ao_charges
63 USE orbital_pointers, ONLY: ncoset
64 USE particle_types, ONLY: particle_type
65 USE qs_charge_mixing, ONLY: charge_mixing
66 USE qs_condnum, ONLY: overlap_condnum
67 USE qs_energy_types, ONLY: qs_energy_type
68 USE qs_environment_types, ONLY: get_qs_env, qs_environment_type
69 USE qs_force_types, ONLY: qs_force_type
70 USE qs_integral_utils, ONLY: basis_set_list_setup, get_memory_usage
71 USE qs_kind_types, ONLY: get_qs_kind, qs_kind_type, get_qs_kind_set
72 USE qs_ks_types, ONLY: qs_ks_env_type, set_ks_env
73 USE qs_neighbor_list_types, ONLY: neighbor_list_iterator_create, neighbor_list_iterate, &
74 get_iterator_info, neighbor_list_set_p_type, &
75 neighbor_list_iterator_p_type, neighbor_list_iterator_release
76 USE qs_overlap, ONLY: create_sab_matrix
77 USE qs_rho_types, ONLY: qs_rho_get, qs_rho_type
78 USE qs_scf_types, ONLY: qs_scf_env_type
79 USE input_section_types, ONLY: section_vals_get_subs_vals, section_vals_type
80 USE string_utilities, ONLY: integer_to_string
81 USE tblite_types, ONLY: tblite_type, deallocate_tblite_type, allocate_tblite_type
82 USE virial_types, ONLY: virial_type
83 USE xtb_types, ONLY: get_xtb_atom_param, xtb_atom_type
84 USE xtb_types, ONLY: xtb_atom_type
85
86#include "./base/base_uses.f90"
87 IMPLICIT NONE
88
89 PRIVATE
90
91 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tblite_interface'
92
93 INTEGER, PARAMETER :: dip_n = 3
94 INTEGER, PARAMETER :: quad_n = 6
95
96 PUBLIC :: tb_init_geometry, tb_init_ham, tb_init_wf
97 PUBLIC :: tb_set_calculator, tb_get_basis, build_tblite_matrices
98 PUBLIC :: tb_get_energy, tb_update_charges, tb_ham_add_coulomb
99 PUBLIC :: tb_get_multipole
100
101CONTAINS
102
103! **************************************************************************************************
104!> \brief ...
105!> \param qs_env ...
106!> \param tb ...
107! **************************************************************************************************
108 SUBROUTINE tb_init_geometry(qs_env, tb)
109
110 TYPE(qs_environment_type), POINTER :: qs_env
111 TYPE(tblite_type), POINTER :: tb
112
113#if defined(__TBLITE)
114
115 CHARACTER(LEN=*), PARAMETER :: routinen = 'tblite_init_geometry'
116
117 TYPE(cell_type), POINTER :: cell
118 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
119 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
120 INTEGER :: iatom, natom
121 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: xyz
122 INTEGER :: handle, ikind
123 INTEGER, DIMENSION(3) :: periodic
124 LOGICAL, DIMENSION(3) :: lperiod
125
126 CALL timeset(routinen, handle)
127
128 !get info from environment vaiarable
129 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell, qs_kind_set=qs_kind_set)
130
131 !get information about particles
132 natom = SIZE(particle_set)
133 ALLOCATE (xyz(3, natom))
134 CALL allocate_tblite_type(tb)
135 ALLOCATE (tb%el_num(natom))
136 tb%el_num = -9
137 DO iatom = 1, natom
138 xyz(:, iatom) = particle_set(iatom)%r(:)
139 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
140 CALL get_qs_kind(qs_kind_set(ikind), zatom=tb%el_num(iatom))
141 IF (tb%el_num(iatom) < 1 .OR. tb%el_num(iatom) > 85) THEN
142 cpabort("only elements 1-85 are supported by tblite")
143 END IF
144 END DO
145
146 !get information about cell / lattice
147 CALL get_cell(cell=cell, periodic=periodic)
148 lperiod(1) = periodic(1) == 1
149 lperiod(2) = periodic(2) == 1
150 lperiod(3) = periodic(3) == 1
151
152 !prepare for the call to the dispersion function
153 CALL new(tb%mol, tb%el_num, xyz, lattice=cell%hmat, periodic=lperiod)
154
155 DEALLOCATE (xyz)
156
157 CALL timestop(handle)
158
159#else
160 mark_used(qs_env)
161 mark_used(tb)
162 cpabort("Built without TBLITE")
163#endif
164
165 END SUBROUTINE tb_init_geometry
166
167 ! **************************************************************************************************
168!> \brief ...
169!> \param tb ...
170!> \param do_grad ...
171! **************************************************************************************************
172 SUBROUTINE tb_init_wf(tb, do_grad)
173
174 TYPE(tblite_type), POINTER :: tb
175 LOGICAL, INTENT(in) :: do_grad
176
177#if defined(__TBLITE)
178
179 INTEGER, PARAMETER :: nspin = 1 !number of spin channels
180
181 TYPE(scf_info) :: info
182
183 info = tb%calc%variable_info()
184 IF (info%charge > shell_resolved) cpabort("tblite: no support for orbital resolved charge")
185 IF (info%dipole > atom_resolved) cpabort("tblite: no support for shell resolved dipole moment")
186 IF (info%quadrupole > atom_resolved) &
187 cpabort("tblite: no support shell resolved quadrupole moment")
188
189 CALL new_wavefunction(tb%wfn, tb%mol%nat, tb%calc%bas%nsh, tb%calc%bas%nao, nspin, 0.0_dp)
190
191 CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
192
193 !allocate quantities later required
194 ALLOCATE (tb%e_hal(tb%mol%nat), tb%e_rep(tb%mol%nat), tb%e_disp(tb%mol%nat))
195 ALLOCATE (tb%e_scd(tb%mol%nat), tb%e_es(tb%mol%nat))
196 ALLOCATE (tb%selfenergy(tb%calc%bas%nsh))
197
198 IF (ALLOCATED(tb%calc%ncoord)) ALLOCATE (tb%cn(tb%mol%nat))
199
200 IF (do_grad) ALLOCATE (tb%grad(3, tb%mol%nat))
201
202 IF (ALLOCATED(tb%grad)) THEN
203 IF (ALLOCATED(tb%calc%ncoord)) THEN
204 ALLOCATE (tb%dcndr(3, tb%mol%nat, tb%mol%nat), tb%dcndL(3, 3, tb%mol%nat))
205 END IF
206 ALLOCATE (tb%dsedcn(tb%calc%bas%nsh))
207 END IF
208
209#else
210 mark_used(tb)
211 mark_used(do_grad)
212 cpabort("Built without TBLITE")
213#endif
214
215 END SUBROUTINE tb_init_wf
216
217! **************************************************************************************************
218!> \brief ...
219!> \param tb ...
220!> \param typ ...
221! **************************************************************************************************
222 SUBROUTINE tb_set_calculator(tb, typ)
223
224 TYPE(tblite_type), POINTER :: tb
225 INTEGER :: typ
226
227#if defined(__TBLITE)
228
229 TYPE(error_type), ALLOCATABLE :: error
230
231 SELECT CASE (typ)
232 CASE default
233 cpabort("Unknown xtb type")
234 CASE (gfn1xtb)
235 CALL new_gfn1_calculator(tb%calc, tb%mol, error)
236 CASE (gfn2xtb)
237 CALL new_gfn2_calculator(tb%calc, tb%mol, error)
238 CASE (ipea1xtb)
239 CALL new_ipea1_calculator(tb%calc, tb%mol, error)
240 END SELECT
241
242#else
243 mark_used(tb)
244 mark_used(typ)
245 cpabort("Built without TBLITE")
246#endif
247
248 END SUBROUTINE tb_set_calculator
249
250! **************************************************************************************************
251!> \brief ...
252!> \param tb ...
253! **************************************************************************************************
254 SUBROUTINE tb_init_ham(tb)
255
256 TYPE(tblite_type), POINTER :: tb
257
258#if defined(__TBLITE)
259
260 TYPE(container_cache) :: hcache, rcache
261
262 tb%e_hal = 0.0_dp
263 tb%e_rep = 0.0_dp
264 tb%e_disp = 0.0_dp
265 IF (ALLOCATED(tb%grad)) THEN
266 tb%grad = 0.0_dp
267 tb%sigma = 0.0_dp
268 END IF
269
270 IF (ALLOCATED(tb%calc%halogen)) THEN
271 CALL tb%calc%halogen%update(tb%mol, hcache)
272 IF (ALLOCATED(tb%grad)) THEN
273 CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal, &
274 & tb%grad, tb%sigma)
275 ELSE
276 CALL tb%calc%halogen%get_engrad(tb%mol, hcache, tb%e_hal)
277 END IF
278 END IF
279
280 IF (ALLOCATED(tb%calc%repulsion)) THEN
281 CALL tb%calc%repulsion%update(tb%mol, rcache)
282 IF (ALLOCATED(tb%grad)) THEN
283 CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep, &
284 & tb%grad, tb%sigma)
285 ELSE
286 CALL tb%calc%repulsion%get_engrad(tb%mol, rcache, tb%e_rep)
287 END IF
288 END IF
289
290 IF (ALLOCATED(tb%calc%dispersion)) THEN
291 CALL tb%calc%dispersion%update(tb%mol, tb%dcache)
292 IF (ALLOCATED(tb%grad)) THEN
293 CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp, &
294 & tb%grad, tb%sigma)
295 ELSE
296 CALL tb%calc%dispersion%get_engrad(tb%mol, tb%dcache, tb%e_disp)
297 END IF
298 END IF
299
300 CALL new_potential(tb%pot, tb%mol, tb%calc%bas, tb%wfn%nspin)
301 IF (ALLOCATED(tb%calc%coulomb)) THEN
302 CALL tb%calc%coulomb%update(tb%mol, tb%cache)
303 END IF
304
305 IF (ALLOCATED(tb%grad)) THEN
306 IF (ALLOCATED(tb%calc%ncoord)) THEN
307 CALL tb%calc%ncoord%get_cn(tb%mol, tb%cn, tb%dcndr, tb%dcndL)
308 END IF
309 CALL get_selfenergy(tb%calc%h0, tb%mol%id, tb%calc%bas%ish_at, &
310 & tb%calc%bas%nsh_id, cn=tb%cn, selfenergy=tb%selfenergy, dsedcn=tb%dsedcn)
311 ELSE
312 IF (ALLOCATED(tb%calc%ncoord)) THEN
313 CALL tb%calc%ncoord%get_cn(tb%mol, tb%cn)
314 END IF
315 CALL get_selfenergy(tb%calc%h0, tb%mol%id, tb%calc%bas%ish_at, &
316 & tb%calc%bas%nsh_id, cn=tb%cn, selfenergy=tb%selfenergy, dsedcn=tb%dsedcn)
317 END IF
318
319#else
320 mark_used(tb)
321 cpabort("Built without TBLITE")
322#endif
323
324 END SUBROUTINE tb_init_ham
325
326! **************************************************************************************************
327!> \brief ...
328!> \param qs_env ...
329!> \param tb ...
330!> \param energy ...
331! **************************************************************************************************
332 SUBROUTINE tb_get_energy(qs_env, tb, energy)
333
334 TYPE(qs_environment_type), POINTER :: qs_env
335 TYPE(tblite_type), POINTER :: tb
336 TYPE(qs_energy_type), POINTER :: energy
337
338#if defined(__TBLITE)
339
340 INTEGER :: iounit
341 TYPE(cp_logger_type), POINTER :: logger
342 TYPE(section_vals_type), POINTER :: scf_section
343
344 NULLIFY (scf_section, logger)
345
346 logger => cp_get_default_logger()
347 iounit = cp_logger_get_default_io_unit(logger)
348 scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
349
350 energy%repulsive = sum(tb%e_rep)
351 energy%el_stat = sum(tb%e_es)
352 energy%dispersion = sum(tb%e_disp)
353 energy%dispersion_sc = sum(tb%e_scd)
354 energy%xtb_xb_inter = sum(tb%e_hal)
355
356 energy%total = energy%core + energy%repulsive + energy%el_stat + energy%dispersion &
357 + energy%dispersion_sc + energy%xtb_xb_inter + energy%kTS &
358 + energy%efield + energy%qmmm_el
359
360 iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
361 extension=".scfLog")
362 IF (iounit > 0) THEN
363 WRITE (unit=iounit, fmt="(/,(T9,A,T60,F20.10))") &
364 "Repulsive pair potential energy: ", energy%repulsive, &
365 "Zeroth order Hamiltonian energy: ", energy%core, &
366 "Electrostatic energy: ", energy%el_stat, &
367 "Non-self consistent dispersion energy: ", energy%dispersion
368 WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
369 "Self-consistent dispersion energy: ", energy%dispersion_sc
370 WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
371 "Correction for halogen bonding: ", energy%xtb_xb_inter
372 IF (abs(energy%efield) > 1.e-12_dp) THEN
373 WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
374 "Electric field interaction energy: ", energy%efield
375 END IF
376 IF (qs_env%qmmm) THEN
377 WRITE (unit=iounit, fmt="(T9,A,T60,F20.10)") &
378 "QM/MM Electrostatic energy: ", energy%qmmm_el
379 END IF
380 END IF
381 CALL cp_print_key_finished_output(iounit, logger, scf_section, &
382 "PRINT%DETAILED_ENERGY")
383
384#else
385 mark_used(qs_env)
386 mark_used(tb)
387 mark_used(energy)
388 cpabort("Built without TBLITE")
389#endif
390
391 END SUBROUTINE tb_get_energy
392
393! **************************************************************************************************
394!> \brief ...
395!> \param tb ...
396!> \param gto_basis_set ...
397!> \param element_symbol ...
398!> \param param ...
399!> \param occ ...
400! **************************************************************************************************
401 SUBROUTINE tb_get_basis(tb, gto_basis_set, element_symbol, param, occ)
402
403 TYPE(tblite_type), POINTER :: tb
404 TYPE(gto_basis_set_type), POINTER :: gto_basis_set
405 CHARACTER(len=2), INTENT(IN) :: element_symbol
406 TYPE(xtb_atom_type), POINTER :: param
407 INTEGER, DIMENSION(5), INTENT(out) :: occ
408
409#if defined(__TBLITE)
410
411 CHARACTER(LEN=default_string_length) :: sng
412 INTEGER :: ang, i_type, id_atom, ind_ao, ipgf, ish, &
413 ishell, ityp, maxl, mprim, natorb, &
414 nset, nshell
415 LOGICAL :: do_ortho
416
417 CALL allocate_gto_basis_set(gto_basis_set)
418
419 !identifying element in the bas data
420 CALL symbol_to_number(i_type, element_symbol)
421 DO id_atom = 1, tb%mol%nat
422 IF (i_type == tb%el_num(id_atom)) EXIT
423 END DO
424 param%z = i_type
425 param%symbol = element_symbol
426 param%defined = .true.
427 ityp = tb%mol%id(id_atom)
428
429 !getting size information
430 nset = tb%calc%bas%nsh_id(ityp)
431 nshell = 1
432 mprim = 0
433 DO ishell = 1, nset
434 mprim = max(mprim, tb%calc%bas%cgto(ishell, ityp)%nprim)
435 END DO
436 param%nshell = nset
437 natorb = 0
438
439 !write basis set information
440 CALL integer_to_string(mprim, sng)
441 gto_basis_set%name = element_symbol//"_STO-"//trim(sng)//"G"
442 gto_basis_set%nset = nset
443 CALL reallocate(gto_basis_set%lmax, 1, nset)
444 CALL reallocate(gto_basis_set%lmin, 1, nset)
445 CALL reallocate(gto_basis_set%npgf, 1, nset)
446 CALL reallocate(gto_basis_set%nshell, 1, nset)
447 CALL reallocate(gto_basis_set%n, 1, 1, 1, nset)
448 CALL reallocate(gto_basis_set%l, 1, 1, 1, nset)
449 CALL reallocate(gto_basis_set%zet, 1, mprim, 1, nset)
450 CALL reallocate(gto_basis_set%gcc, 1, mprim, 1, 1, 1, nset)
451
452 ind_ao = 0
453 maxl = 0
454 DO ishell = 1, nset
455 ang = tb%calc%bas%cgto(ishell, ityp)%ang
456 natorb = natorb + (2*ang + 1)
457 param%lval(ishell) = ang
458 maxl = max(ang, maxl)
459 gto_basis_set%lmax(ishell) = ang
460 gto_basis_set%lmin(ishell) = ang
461 gto_basis_set%npgf(ishell) = tb%calc%bas%cgto(ishell, ityp)%nprim
462 gto_basis_set%nshell(ishell) = nshell
463 gto_basis_set%n(1, ishell) = ang + 1
464 gto_basis_set%l(1, ishell) = ang
465 DO ipgf = 1, gto_basis_set%npgf(ishell)
466 gto_basis_set%gcc(ipgf, 1, ishell) = tb%calc%bas%cgto(ishell, ityp)%coeff(ipgf)
467 gto_basis_set%zet(ipgf, ishell) = tb%calc%bas%cgto(ishell, ityp)%alpha(ipgf)
468 END DO
469 DO ipgf = 1, (2*ang + 1)
470 ind_ao = ind_ao + 1
471 param%lao(ind_ao) = ang
472 param%nao(ind_ao) = ishell
473 END DO
474 END DO
475
476 do_ortho = .false.
477 CALL process_gto_basis(gto_basis_set, do_ortho, nset, maxl)
478
479 !setting additional values in parameter
480 param%rcut = get_cutoff(tb%calc%bas)
481 param%natorb = natorb
482 param%lmax = maxl !max angular momentum
483
484 !getting occupation
485 occ = 0
486 DO ish = 1, min(tb%calc%bas%nsh_at(id_atom), 5)
487 occ(ish) = nint(tb%calc%h0%refocc(ish, ityp))
488 param%occupation(ish) = occ(ish)
489 END DO
490 param%zeff = sum(occ) !effective core charge
491
492 !set normalization process
493 gto_basis_set%norm_type = 3
494
495#else
496 mark_used(tb)
497 mark_used(gto_basis_set)
498 mark_used(element_symbol)
499 mark_used(param)
500 mark_used(occ)
501 cpabort("Built without TBLITE")
502#endif
503
504 END SUBROUTINE tb_get_basis
505
506 ! **************************************************************************************************
507!> \brief ...
508!> \param qs_env ...
509!> \param calculate_forces ...
510! **************************************************************************************************
511 SUBROUTINE build_tblite_matrices(qs_env, calculate_forces)
512
513 TYPE(qs_environment_type), POINTER :: qs_env
514 LOGICAL, INTENT(IN) :: calculate_forces
515
516#if defined(__TBLITE)
517
518 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_tblite_matrices'
519
520 INTEGER :: handle, maxder, nderivatives, nimg, img, nkind, i, ic, iw, &
521 iatom, jatom, ikind, jkind, iset, jset, n1, n2, icol, irow, &
522 ishell, jshell, ia, ib, sgfa, sgfb, atom_a, atom_b, &
523 ldsab, nseta, nsetb, natorb_a, natorb_b
524 LOGICAL :: found, norml1, norml2, use_arnoldi, use_virial
525 REAL(kind=dp) :: dr, rr
526 INTEGER, DIMENSION(3) :: cell
527 REAL(kind=dp) :: hij, shpoly
528 REAL(kind=dp), DIMENSION(2) :: condnum
529 REAL(kind=dp), DIMENSION(3) :: rij
530 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
531 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
532 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint, hint
533 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min
534 INTEGER, DIMENSION(:), POINTER :: npgfa, npgfb, nsgfa, nsgfb
535 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
536 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
537 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, zeta, zetb, scon_a, scon_b
538 REAL(kind=dp), DIMENSION(:, :), POINTER :: sblock, fblock
539
540 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
541 TYPE(atprop_type), POINTER :: atprop
542 TYPE(cp_blacs_env_type), POINTER :: blacs_env
543 TYPE(cp_logger_type), POINTER :: logger
544 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s, matrix_p, matrix_w
545 TYPE(dft_control_type), POINTER :: dft_control
546 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
547 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
548 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
549 TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
550 TYPE(neighbor_list_iterator_p_type), &
551 DIMENSION(:), POINTER :: nl_iterator
552 TYPE(qs_energy_type), POINTER :: energy
553 TYPE(qs_ks_env_type), POINTER :: ks_env
554 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
555 TYPE(qs_rho_type), POINTER :: rho
556 TYPE(tblite_type), POINTER :: tb
557 TYPE(tb_hamiltonian), POINTER :: h0
558 TYPE(virial_type), POINTER :: virial
559
560 CALL timeset(routinen, handle)
561
562 NULLIFY (ks_env, energy, atomic_kind_set, qs_kind_set)
563 NULLIFY (matrix_h, matrix_s, atprop, dft_control)
564 NULLIFY (sab_orb, rho, tb)
565
566 CALL get_qs_env(qs_env=qs_env, &
567 ks_env=ks_env, &
568 energy=energy, &
569 atomic_kind_set=atomic_kind_set, &
570 qs_kind_set=qs_kind_set, &
571 matrix_h_kp=matrix_h, &
572 matrix_s_kp=matrix_s, &
573 atprop=atprop, &
574 dft_control=dft_control, &
575 sab_orb=sab_orb, &
576 rho=rho, tb_tblite=tb)
577 h0 => tb%calc%h0
578
579 nkind = SIZE(atomic_kind_set)
580 nderivatives = 0
581 IF (calculate_forces) nderivatives = 1
582 maxder = ncoset(nderivatives)
583 nimg = dft_control%nimages
584
585 ! get density matrtix
586 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
587
588 ! set up matrices for force calculations
589 IF (calculate_forces) THEN
590 NULLIFY (force, matrix_w)
591 CALL get_qs_env(qs_env=qs_env, &
592 matrix_w_kp=matrix_w, &
593 virial=virial, force=force)
594
595 IF (SIZE(matrix_p, 1) == 2) THEN
596 DO img = 1, nimg
597 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
598 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
599 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
600 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
601 END DO
602 END IF
603 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
604 END IF
605
606 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
607
608 ! set up basis set lists
609 ALLOCATE (basis_set_list(nkind))
610 CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
611
612 ! allocate overlap matrix
613 CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
614 CALL create_sab_matrix(ks_env, matrix_s, "OVERLAP MATRIX", basis_set_list, basis_set_list, &
615 sab_orb, .true.)
616 CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
617
618 ! initialize H matrix
619 CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
620 DO img = 1, nimg
621 ALLOCATE (matrix_h(1, img)%matrix)
622 CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
623 name="HAMILTONIAN MATRIX")
624 CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
625 END DO
626 CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
627
628 ! loop over all atom pairs with a non-zero overlap (sab_orb)
629 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
630 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
631 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
632 iatom=iatom, jatom=jatom, r=rij, cell=cell)
633
634 icol = max(iatom, jatom)
635 irow = min(iatom, jatom)
636 IF (icol == jatom) THEN
637 rij = -rij
638 i = ikind
639 ikind = jkind
640 jkind = i
641 END IF
642
643 dr = norm2(rij(:))
644
645 ic = 1
646 NULLIFY (sblock)
647 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
648 row=irow, col=icol, block=sblock, found=found)
649 cpassert(found)
650 NULLIFY (fblock)
651 CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
652 row=irow, col=icol, block=fblock, found=found)
653 cpassert(found)
654
655 ! --------- Overlap
656 !get basis information
657 basis_set_a => basis_set_list(ikind)%gto_basis_set
658 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
659 basis_set_b => basis_set_list(jkind)%gto_basis_set
660 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
661 atom_a = atom_of_kind(icol)
662 atom_b = atom_of_kind(irow)
663 ! basis a
664 first_sgfa => basis_set_a%first_sgf
665 la_max => basis_set_a%lmax
666 la_min => basis_set_a%lmin
667 npgfa => basis_set_a%npgf
668 nseta = basis_set_a%nset
669 nsgfa => basis_set_a%nsgf_set
670 rpgfa => basis_set_a%pgf_radius
671 set_radius_a => basis_set_a%set_radius
672 scon_a => basis_set_a%scon
673 zeta => basis_set_a%zet
674 ! basis b
675 first_sgfb => basis_set_b%first_sgf
676 lb_max => basis_set_b%lmax
677 lb_min => basis_set_b%lmin
678 npgfb => basis_set_b%npgf
679 nsetb = basis_set_b%nset
680 nsgfb => basis_set_b%nsgf_set
681 rpgfb => basis_set_b%pgf_radius
682 set_radius_b => basis_set_b%set_radius
683 scon_b => basis_set_b%scon
684 zetb => basis_set_b%zet
685
686 ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
687 ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
688 natorb_a = 0
689 DO iset = 1, nseta
690 natorb_a = natorb_a + (2*basis_set_a%l(1, iset) + 1)
691 END DO
692 natorb_b = 0
693 DO iset = 1, nsetb
694 natorb_b = natorb_b + (2*basis_set_b%l(1, iset) + 1)
695 END DO
696 ALLOCATE (sint(natorb_a, natorb_b, maxder))
697 sint = 0.0_dp
698 ALLOCATE (hint(natorb_a, natorb_b, maxder))
699 hint = 0.0_dp
700
701 !----------------- overlap integrals
702 DO iset = 1, nseta
703 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
704 sgfa = first_sgfa(1, iset)
705 DO jset = 1, nsetb
706 IF (set_radius_a(iset) + set_radius_b(jset) < dr) cycle
707 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
708 sgfb = first_sgfb(1, jset)
709 IF (calculate_forces) THEN
710 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
711 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
712 rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
713 ELSE
714 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
715 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
716 rij, sab=oint(:, :, 1))
717 END IF
718 ! Contraction
719 CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
720 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
721 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.false.)
722 IF (calculate_forces) THEN
723 DO i = 2, 4
724 CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
725 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
726 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.false.)
727 END DO
728 END IF
729 END DO
730 END DO
731
732 ! update S matrix
733 IF (icol <= irow) THEN
734 sblock(:, :) = sblock(:, :) + sint(:, :, 1)
735 ELSE
736 sblock(:, :) = sblock(:, :) + transpose(sint(:, :, 1))
737 END IF
738
739 ! --------- Hamiltonian
740 IF (icol == irow .AND. dr < 0.001_dp) THEN
741 !get diagonal F matrix from selfenergy
742 n1 = 0
743 DO ishell = 1, icol - 1
744 n1 = n1 + tb%calc%bas%nsh_at(ishell)
745 END DO
746 DO iset = 1, nseta
747 sgfa = first_sgfa(1, iset)
748 hij = tb%selfenergy(n1 + iset)
749 DO ia = sgfa, sgfa + nsgfa(iset) - 1
750 hint(ia, ia, 1) = hij
751 END DO
752 END DO
753 ELSE
754 !get off-diagonal F matrix
755 rr = sqrt(dr/(h0%rad(jkind) + h0%rad(ikind)))
756 n1 = 0
757 DO ishell = 1, icol - 1
758 n1 = n1 + tb%calc%bas%nsh_at(ishell)
759 END DO
760 DO iset = 1, nseta
761 sgfa = first_sgfa(1, iset)
762 n2 = 0
763 DO jshell = 1, irow - 1
764 n2 = n2 + tb%calc%bas%nsh_at(jshell)
765 END DO
766 DO jset = 1, nsetb
767 sgfb = first_sgfb(1, jset)
768 shpoly = (1.0_dp + h0%shpoly(iset, ikind)*rr) &
769 *(1.0_dp + h0%shpoly(jset, jkind)*rr)
770 hij = 0.5_dp*(tb%selfenergy(n1 + iset) + tb%selfenergy(n2 + jset)) &
771 *h0%hscale(iset, jset, ikind, jkind)*shpoly
772 DO ia = sgfa, sgfa + nsgfa(iset) - 1
773 DO ib = sgfb, sgfb + nsgfb(jset) - 1
774 hint(ia, ib, 1) = hij*sint(ia, ib, 1)
775 END DO
776 END DO
777 END DO
778 END DO
779 END IF
780
781 ! update F matrix
782 IF (icol <= irow) THEN
783 fblock(:, :) = fblock(:, :) + hint(:, :, 1)
784 ELSE
785 fblock(:, :) = fblock(:, :) + transpose(hint(:, :, 1))
786 END IF
787
788 DEALLOCATE (oint, owork, sint, hint)
789
790 END DO
791 CALL neighbor_list_iterator_release(nl_iterator)
792
793 DO i = 1, SIZE(matrix_s, 1)
794 DO img = 1, nimg
795 CALL dbcsr_finalize(matrix_s(i, img)%matrix)
796 CALL dbcsr_finalize(matrix_h(i, img)%matrix)
797 END DO
798 END DO
799
800 !compute multipole moments for gfn2
801 IF (dft_control%qs_control%xtb_control%tblite_method == gfn2xtb) &
802 CALL tb_get_multipole(qs_env, tb)
803
804 ! output overlap information
805 NULLIFY (logger)
806 logger => cp_get_default_logger()
807 IF (.NOT. calculate_forces) THEN
808 IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
809 "DFT%PRINT%OVERLAP_CONDITION") .NE. 0) THEN
810 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
811 extension=".Log")
812 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
813 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
814 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
815 CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
816 CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
817 END IF
818 END IF
819
820 DEALLOCATE (basis_set_list)
821
822 CALL timestop(handle)
823
824#else
825 mark_used(qs_env)
826 mark_used(calculate_forces)
827 cpabort("Built without TBLITE")
828#endif
829
830 END SUBROUTINE build_tblite_matrices
831
832 ! **************************************************************************************************
833!> \brief ...
834!> \param qs_env ...
835!> \param dft_control ...
836!> \param tb ...
837!> \param calculate_forces ...
838!> \param use_rho ...
839! **************************************************************************************************
840 SUBROUTINE tb_update_charges(qs_env, dft_control, tb, calculate_forces, use_rho)
841
842 TYPE(qs_environment_type), POINTER :: qs_env
843 TYPE(dft_control_type), POINTER :: dft_control
844 TYPE(tblite_type), POINTER :: tb
845 LOGICAL, INTENT(IN) :: calculate_forces
846 LOGICAL, INTENT(IN) :: use_rho
847
848#if defined(__TBLITE)
849
850 INTEGER :: iatom, ikind, is, ns, atom_a, ii, im
851 INTEGER :: nimg, nkind, nsgf, natorb, na
852 INTEGER :: n_atom, max_orb, max_shell
853 LOGICAL :: do_dipole, do_quadrupole
854 REAL(kind=dp) :: norm
855 INTEGER, DIMENSION(5) :: occ
856 INTEGER, DIMENSION(25) :: lao
857 INTEGER, DIMENSION(25) :: nao
858 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ch_atom, ch_shell, ch_ref, ch_orb
859 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, ao_dip, ao_quad
860
861 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
862 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_p
863 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
864 TYPE(dbcsr_type), POINTER :: s_matrix
865 TYPE(mp_para_env_type), POINTER :: para_env
866 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
867 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
868 TYPE(qs_rho_type), POINTER :: rho
869 TYPE(qs_scf_env_type), POINTER :: scf_env
870 TYPE(xtb_atom_type), POINTER :: xtb_kind
871
872 IF (calculate_forces) cpabort("tblite: forces not yet available")
873
874 ! also compute multipoles needed by GFN2
875 do_dipole = .false.
876 do_quadrupole = .false.
877
878 ! compute mulliken charges required for charge update
879 NULLIFY (particle_set, qs_kind_set, atomic_kind_set)
880 CALL get_qs_env(qs_env=qs_env, scf_env=scf_env, particle_set=particle_set, qs_kind_set=qs_kind_set, &
881 atomic_kind_set=atomic_kind_set, matrix_s_kp=matrix_s, rho=rho, para_env=para_env)
882 NULLIFY (matrix_p)
883 IF (use_rho) THEN
884 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
885 IF (ASSOCIATED(tb%dipbra)) do_dipole = .true.
886 IF (ASSOCIATED(tb%quadbra)) do_quadrupole = .true.
887 ELSE
888 matrix_p => scf_env%p_mix_new
889 END IF
890 n_atom = SIZE(particle_set)
891 nkind = SIZE(atomic_kind_set)
892 nimg = dft_control%nimages
893 CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
894 ALLOCATE (aocg(nsgf, n_atom))
895 aocg = 0.0_dp
896 IF (do_dipole) ALLOCATE (ao_dip(n_atom, dip_n))
897 IF (do_quadrupole) ALLOCATE (ao_quad(n_atom, quad_n))
898 max_orb = 0
899 max_shell = 0
900 DO ikind = 1, nkind
901 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
902 CALL get_xtb_atom_param(xtb_kind, natorb=natorb)
903 max_orb = max(max_orb, natorb)
904 END DO
905 DO is = 1, n_atom
906 max_shell = max(max_shell, tb%calc%bas%nsh_at(is))
907 END DO
908 ALLOCATE (ch_atom(n_atom, 1), ch_shell(n_atom, max_shell))
909 ALLOCATE (ch_orb(max_orb, n_atom), ch_ref(max_orb, n_atom))
910 ch_atom = 0.0_dp
911 ch_shell = 0.0_dp
912 ch_orb = 0.0_dp
913 ch_ref = 0.0_dp
914 IF (nimg > 1) THEN
915 CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
916 IF (do_dipole .OR. do_quadrupole) THEN
917 cpabort("missing contraction with density matrix for multiple k-points")
918 END IF
919 ELSE
920 NULLIFY (p_matrix, s_matrix)
921 p_matrix => matrix_p(:, 1)
922 s_matrix => matrix_s(1, 1)%matrix
923 CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
924 IF (do_dipole) THEN
925 DO im = 1, dip_n
926 CALL contract_dens(p_matrix, tb%dipbra(im)%matrix, tb%dipket(im)%matrix, ao_dip(:, im), para_env)
927 END DO
928 END IF
929 IF (do_quadrupole) THEN
930 DO im = 1, quad_n
931 CALL contract_dens(p_matrix, tb%quadbra(im)%matrix, tb%quadket(im)%matrix, ao_quad(:, im), para_env)
932 END DO
933 END IF
934 END IF
935 NULLIFY (xtb_kind)
936 DO ikind = 1, nkind
937 CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
938 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
939 CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, nao=nao, occupation=occ)
940 DO iatom = 1, na
941 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
942 DO is = 1, natorb
943 ns = lao(is) + 1
944 norm = 2*lao(is) + 1
945 ch_ref(is, atom_a) = tb%calc%h0%refocc(nao(is), ikind)/norm
946 ch_orb(is, atom_a) = aocg(is, atom_a) - ch_ref(is, atom_a)
947 ch_shell(atom_a, ns) = ch_orb(is, atom_a) + ch_shell(atom_a, ns)
948 END DO
949 ch_atom(atom_a, 1) = sum(ch_orb(:, atom_a))
950 END DO
951 END DO
952 DEALLOCATE (aocg)
953
954 ! charge mixing
955 IF (dft_control%qs_control%do_ls_scf) THEN
956 !
957 ELSE
958 CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
959 ch_shell, para_env, scf_env%iter_count)
960 END IF
961
962 !setting new wave function
963 CALL tb%pot%reset
964 tb%e_es = 0.0_dp
965 tb%e_scd = 0.0_dp
966 DO iatom = 1, n_atom
967 ii = tb%calc%bas%ish_at(iatom)
968 DO is = 1, tb%calc%bas%nsh_at(iatom)
969 tb%wfn%qsh(ii + is, 1) = -ch_shell(iatom, is)
970 END DO
971 tb%wfn%qat(iatom, 1) = -ch_atom(iatom, 1)
972 END DO
973
974 IF (do_dipole) THEN
975 DO iatom = 1, n_atom
976 DO im = 1, dip_n
977 tb%wfn%dpat(im, iatom, 1) = -ao_dip(iatom, im)
978 END DO
979 END DO
980 DEALLOCATE (ao_dip)
981 END IF
982 IF (do_quadrupole) THEN
983 DO iatom = 1, n_atom
984 DO im = 1, quad_n
985 tb%wfn%qpat(im, iatom, 1) = -ao_quad(iatom, im)
986 END DO
987 END DO
988 DEALLOCATE (ao_quad)
989 END IF
990
991 IF (ALLOCATED(tb%calc%coulomb)) THEN
992 CALL tb%calc%coulomb%get_potential(tb%mol, tb%cache, tb%wfn, tb%pot)
993 END IF
994 IF (ALLOCATED(tb%calc%dispersion)) THEN
995 CALL tb%calc%dispersion%get_potential(tb%mol, tb%dcache, tb%wfn, tb%pot)
996 END IF
997
998 IF (ALLOCATED(tb%calc%coulomb)) THEN
999 CALL tb%calc%coulomb%get_energy(tb%mol, tb%cache, tb%wfn, tb%e_es)
1000 END IF
1001 IF (ALLOCATED(tb%calc%dispersion)) THEN
1002 CALL tb%calc%dispersion%get_energy(tb%mol, tb%dcache, tb%wfn, tb%e_scd)
1003 END IF
1004
1005 DEALLOCATE (ch_atom, ch_shell, ch_orb, ch_ref)
1006
1007#else
1008 mark_used(qs_env)
1009 mark_used(tb)
1010 mark_used(dft_control)
1011 mark_used(calculate_forces)
1012 mark_used(use_rho)
1013 cpabort("Built without TBLITE")
1014#endif
1015
1016 END SUBROUTINE tb_update_charges
1017
1018 ! **************************************************************************************************
1019!> \brief ...
1020!> \param qs_env ...
1021!> \param tb ...
1022!> \param dft_control ...
1023!> \param calculate_forces ...
1024! **************************************************************************************************
1025 SUBROUTINE tb_ham_add_coulomb(qs_env, tb, dft_control, calculate_forces)
1026
1027 TYPE(qs_environment_type), POINTER :: qs_env
1028 TYPE(tblite_type), POINTER :: tb
1029 TYPE(dft_control_type), POINTER :: dft_control
1030 LOGICAL, INTENT(IN) :: calculate_forces
1031
1032#if defined(__TBLITE)
1033
1034 INTEGER :: ikind, jkind, iatom, jatom, icol, irow
1035 INTEGER :: ic, is, nimg, ni, nj, i, j
1036 INTEGER :: la, lb, za, zb
1037 LOGICAL :: found
1038 INTEGER, DIMENSION(3) :: cellind
1039 INTEGER, DIMENSION(25) :: naoa, naob
1040 REAL(kind=dp), DIMENSION(3) :: rij
1041 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, sum_shell
1042 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ashift, bshift
1043 REAL(kind=dp), DIMENSION(:, :), POINTER :: ksblock, sblock
1044 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1045 REAL(kind=dp), DIMENSION(:, :), POINTER :: dip_ket1, dip_ket2, dip_ket3, &
1046 dip_bra1, dip_bra2, dip_bra3
1047 REAL(kind=dp), DIMENSION(:, :), POINTER :: quad_ket1, quad_ket2, quad_ket3, &
1048 quad_ket4, quad_ket5, quad_ket6, &
1049 quad_bra1, quad_bra2, quad_bra3, &
1050 quad_bra4, quad_bra5, quad_bra6
1051
1052 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1053 TYPE(dbcsr_iterator_type) :: iter
1054 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1055 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
1056 TYPE(kpoint_type), POINTER :: kpoints
1057 TYPE(neighbor_list_iterator_p_type), &
1058 DIMENSION(:), POINTER :: nl_iterator
1059 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1060 POINTER :: n_list
1061 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1062 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
1063
1064 IF (calculate_forces) cpabort("tblite: forces not yet available")
1065
1066 nimg = dft_control%nimages
1067
1068 NULLIFY (matrix_s, ks_matrix, n_list, qs_kind_set)
1069 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)
1070
1071 !creating sum of shell lists
1072 ALLOCATE (sum_shell(tb%mol%nat))
1073 i = 0
1074 DO j = 1, tb%mol%nat
1075 sum_shell(j) = i
1076 i = i + tb%calc%bas%nsh_at(j)
1077 END DO
1078
1079 IF (nimg == 1) THEN
1080 ! no k-points; all matrices have been transformed to periodic bsf
1081 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1082 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1083 kind_of=kind_of)
1084 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
1085 DO WHILE (dbcsr_iterator_blocks_left(iter))
1086 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1087
1088 ikind = kind_of(irow)
1089 jkind = kind_of(icol)
1090
1091 ! atomic parameters
1092 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
1093 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
1094 CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa)
1095 CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob)
1096
1097 ni = SIZE(sblock, 1)
1098 ALLOCATE (ashift(ni, ni))
1099 ashift = 0.0_dp
1100 DO i = 1, ni
1101 la = naoa(i) + sum_shell(irow)
1102 ashift(i, i) = tb%pot%vsh(la, 1)
1103 END DO
1104
1105 nj = SIZE(sblock, 2)
1106 ALLOCATE (bshift(nj, nj))
1107 bshift = 0.0_dp
1108 DO j = 1, nj
1109 lb = naob(j) + sum_shell(icol)
1110 bshift(j, j) = tb%pot%vsh(lb, 1)
1111 END DO
1112
1113 DO is = 1, SIZE(ks_matrix, 1)
1114 NULLIFY (ksblock)
1115 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
1116 row=irow, col=icol, block=ksblock, found=found)
1117 cpassert(found)
1118 ksblock = ksblock - 0.5_dp*(matmul(ashift, sblock) &
1119 + matmul(sblock, bshift))
1120 ksblock = ksblock - 0.5_dp*(tb%pot%vat(irow, 1) &
1121 + tb%pot%vat(icol, 1))*sblock
1122 END DO
1123 DEALLOCATE (ashift, bshift)
1124 END DO
1125 CALL dbcsr_iterator_stop(iter)
1126
1127 IF (ASSOCIATED(tb%dipbra)) THEN
1128 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
1129 DO WHILE (dbcsr_iterator_blocks_left(iter))
1130 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1131
1132 NULLIFY (dip_bra1, dip_bra2, dip_bra3)
1133 CALL dbcsr_get_block_p(matrix=tb%dipbra(1)%matrix, &
1134 row=irow, col=icol, block=dip_bra1, found=found)
1135 cpassert(found)
1136 CALL dbcsr_get_block_p(matrix=tb%dipbra(2)%matrix, &
1137 row=irow, col=icol, block=dip_bra2, found=found)
1138 cpassert(found)
1139 CALL dbcsr_get_block_p(matrix=tb%dipbra(3)%matrix, &
1140 row=irow, col=icol, block=dip_bra3, found=found)
1141 cpassert(found)
1142 NULLIFY (dip_ket1, dip_ket2, dip_ket3)
1143 CALL dbcsr_get_block_p(matrix=tb%dipket(1)%matrix, &
1144 row=irow, col=icol, block=dip_ket1, found=found)
1145 cpassert(found)
1146 CALL dbcsr_get_block_p(matrix=tb%dipket(2)%matrix, &
1147 row=irow, col=icol, block=dip_ket2, found=found)
1148 cpassert(found)
1149 CALL dbcsr_get_block_p(matrix=tb%dipket(3)%matrix, &
1150 row=irow, col=icol, block=dip_ket3, found=found)
1151 cpassert(found)
1152
1153 DO is = 1, SIZE(ks_matrix, 1)
1154 NULLIFY (ksblock)
1155 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
1156 row=irow, col=icol, block=ksblock, found=found)
1157 cpassert(found)
1158 ksblock = ksblock - 0.5_dp*(dip_ket1*tb%pot%vdp(1, irow, 1) &
1159 + dip_ket2*tb%pot%vdp(2, irow, 1) &
1160 + dip_ket3*tb%pot%vdp(3, irow, 1) &
1161 + dip_bra1*tb%pot%vdp(1, icol, 1) &
1162 + dip_bra2*tb%pot%vdp(2, icol, 1) &
1163 + dip_bra3*tb%pot%vdp(3, icol, 1))
1164 END DO
1165 END DO
1166 CALL dbcsr_iterator_stop(iter)
1167 END IF
1168
1169 IF (ASSOCIATED(tb%quadbra)) THEN
1170 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
1171 DO WHILE (dbcsr_iterator_blocks_left(iter))
1172 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
1173
1174 NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
1175 CALL dbcsr_get_block_p(matrix=tb%quadbra(1)%matrix, &
1176 row=irow, col=icol, block=quad_bra1, found=found)
1177 cpassert(found)
1178 CALL dbcsr_get_block_p(matrix=tb%quadbra(2)%matrix, &
1179 row=irow, col=icol, block=quad_bra2, found=found)
1180 cpassert(found)
1181 CALL dbcsr_get_block_p(matrix=tb%quadbra(3)%matrix, &
1182 row=irow, col=icol, block=quad_bra3, found=found)
1183 cpassert(found)
1184 CALL dbcsr_get_block_p(matrix=tb%quadbra(4)%matrix, &
1185 row=irow, col=icol, block=quad_bra4, found=found)
1186 cpassert(found)
1187 CALL dbcsr_get_block_p(matrix=tb%quadbra(5)%matrix, &
1188 row=irow, col=icol, block=quad_bra5, found=found)
1189 cpassert(found)
1190 CALL dbcsr_get_block_p(matrix=tb%quadbra(6)%matrix, &
1191 row=irow, col=icol, block=quad_bra6, found=found)
1192 cpassert(found)
1193
1194 NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
1195 CALL dbcsr_get_block_p(matrix=tb%quadket(1)%matrix, &
1196 row=irow, col=icol, block=quad_ket1, found=found)
1197 cpassert(found)
1198 CALL dbcsr_get_block_p(matrix=tb%quadket(2)%matrix, &
1199 row=irow, col=icol, block=quad_ket2, found=found)
1200 cpassert(found)
1201 CALL dbcsr_get_block_p(matrix=tb%quadket(3)%matrix, &
1202 row=irow, col=icol, block=quad_ket3, found=found)
1203 cpassert(found)
1204 CALL dbcsr_get_block_p(matrix=tb%quadket(4)%matrix, &
1205 row=irow, col=icol, block=quad_ket4, found=found)
1206 cpassert(found)
1207 CALL dbcsr_get_block_p(matrix=tb%quadket(5)%matrix, &
1208 row=irow, col=icol, block=quad_ket5, found=found)
1209 cpassert(found)
1210 CALL dbcsr_get_block_p(matrix=tb%quadket(6)%matrix, &
1211 row=irow, col=icol, block=quad_ket6, found=found)
1212 cpassert(found)
1213
1214 DO is = 1, SIZE(ks_matrix, 1)
1215 NULLIFY (ksblock)
1216 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
1217 row=irow, col=icol, block=ksblock, found=found)
1218 cpassert(found)
1219
1220 ksblock = ksblock - 0.5_dp*(quad_ket1*tb%pot%vqp(1, irow, 1) &
1221 + quad_ket2*tb%pot%vqp(2, irow, 1) &
1222 + quad_ket3*tb%pot%vqp(3, irow, 1) &
1223 + quad_ket4*tb%pot%vqp(4, irow, 1) &
1224 + quad_ket5*tb%pot%vqp(5, irow, 1) &
1225 + quad_ket6*tb%pot%vqp(6, irow, 1) &
1226 + quad_bra1*tb%pot%vqp(1, icol, 1) &
1227 + quad_bra2*tb%pot%vqp(2, icol, 1) &
1228 + quad_bra3*tb%pot%vqp(3, icol, 1) &
1229 + quad_bra4*tb%pot%vqp(4, icol, 1) &
1230 + quad_bra5*tb%pot%vqp(5, icol, 1) &
1231 + quad_bra6*tb%pot%vqp(6, icol, 1))
1232 END DO
1233 END DO
1234 CALL dbcsr_iterator_stop(iter)
1235 END IF
1236
1237 ELSE
1238 cpabort("GFN methods with k-points not tested")
1239 NULLIFY (kpoints)
1240 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
1241 NULLIFY (cell_to_index)
1242 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1243
1244 CALL neighbor_list_iterator_create(nl_iterator, n_list)
1245 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1246 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
1247 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
1248
1249 icol = max(iatom, jatom)
1250 irow = min(iatom, jatom)
1251
1252 IF (icol == jatom) THEN
1253 rij = -rij
1254 i = ikind
1255 ikind = jkind
1256 jkind = i
1257 END IF
1258
1259 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
1260 cpassert(ic > 0)
1261
1262 NULLIFY (sblock)
1263 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
1264 row=irow, col=icol, block=sblock, found=found)
1265 cpassert(found)
1266
1267 ! atomic parameters
1268 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
1269 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
1270 CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa)
1271 CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob)
1272
1273 ni = SIZE(sblock, 1)
1274 ALLOCATE (ashift(ni, ni))
1275 ashift = 0.0_dp
1276 DO i = 1, ni
1277 la = naoa(i) + sum_shell(irow)
1278 ashift(i, i) = tb%pot%vsh(la, 1)
1279 END DO
1280
1281 nj = SIZE(sblock, 2)
1282 ALLOCATE (bshift(nj, nj))
1283 bshift = 0.0_dp
1284 DO j = 1, nj
1285 lb = naob(j) + sum_shell(icol)
1286 bshift(j, j) = tb%pot%vsh(lb, 1)
1287 END DO
1288
1289 DO is = 1, SIZE(ks_matrix, 1)
1290 NULLIFY (ksblock)
1291 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
1292 row=irow, col=icol, block=ksblock, found=found)
1293 cpassert(found)
1294 ksblock = ksblock - 0.5_dp*(matmul(ashift, sblock) &
1295 + matmul(sblock, bshift))
1296 ksblock = ksblock - 0.5_dp*(tb%pot%vat(irow, 1) &
1297 + tb%pot%vat(icol, 1))*sblock
1298 END DO
1299 END DO
1300 CALL neighbor_list_iterator_release(nl_iterator)
1301
1302 IF (ASSOCIATED(tb%dipbra)) THEN
1303 CALL neighbor_list_iterator_create(nl_iterator, n_list)
1304 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1305 CALL get_iterator_info(nl_iterator, &
1306 iatom=iatom, jatom=jatom, cell=cellind)
1307 icol = max(iatom, jatom)
1308 irow = min(iatom, jatom)
1309
1310 NULLIFY (dip_bra1, dip_bra2, dip_bra3)
1311 CALL dbcsr_get_block_p(matrix=tb%dipbra(1)%matrix, &
1312 row=irow, col=icol, block=dip_bra1, found=found)
1313 cpassert(found)
1314 CALL dbcsr_get_block_p(matrix=tb%dipbra(2)%matrix, &
1315 row=irow, col=icol, block=dip_bra2, found=found)
1316 cpassert(found)
1317 CALL dbcsr_get_block_p(matrix=tb%dipbra(3)%matrix, &
1318 row=irow, col=icol, block=dip_bra3, found=found)
1319 cpassert(found)
1320 NULLIFY (dip_ket1, dip_ket2, dip_ket3)
1321 CALL dbcsr_get_block_p(matrix=tb%dipket(1)%matrix, &
1322 row=irow, col=icol, block=dip_ket1, found=found)
1323 cpassert(found)
1324 CALL dbcsr_get_block_p(matrix=tb%dipket(2)%matrix, &
1325 row=irow, col=icol, block=dip_ket2, found=found)
1326 cpassert(found)
1327 CALL dbcsr_get_block_p(matrix=tb%dipket(3)%matrix, &
1328 row=irow, col=icol, block=dip_ket3, found=found)
1329 cpassert(found)
1330
1331 DO is = 1, SIZE(ks_matrix, 1)
1332 NULLIFY (ksblock)
1333 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
1334 row=irow, col=icol, block=ksblock, found=found)
1335 cpassert(found)
1336 ksblock = ksblock - 0.5_dp*(dip_ket1*tb%pot%vdp(1, irow, 1) &
1337 + dip_ket2*tb%pot%vdp(2, irow, 1) &
1338 + dip_ket3*tb%pot%vdp(3, irow, 1) &
1339 + dip_bra1*tb%pot%vdp(1, icol, 1) &
1340 + dip_bra2*tb%pot%vdp(2, icol, 1) &
1341 + dip_bra3*tb%pot%vdp(3, icol, 1))
1342 END DO
1343 END DO
1344 CALL neighbor_list_iterator_release(nl_iterator)
1345 END IF
1346
1347 IF (ASSOCIATED(tb%quadbra)) THEN
1348 CALL neighbor_list_iterator_create(nl_iterator, n_list)
1349 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1350 CALL get_iterator_info(nl_iterator, &
1351 iatom=iatom, jatom=jatom, cell=cellind)
1352 icol = max(iatom, jatom)
1353 irow = min(iatom, jatom)
1354
1355 NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
1356 CALL dbcsr_get_block_p(matrix=tb%quadbra(1)%matrix, &
1357 row=irow, col=icol, block=quad_bra1, found=found)
1358 cpassert(found)
1359 CALL dbcsr_get_block_p(matrix=tb%quadbra(2)%matrix, &
1360 row=irow, col=icol, block=quad_bra2, found=found)
1361 cpassert(found)
1362 CALL dbcsr_get_block_p(matrix=tb%quadbra(3)%matrix, &
1363 row=irow, col=icol, block=quad_bra3, found=found)
1364 cpassert(found)
1365 CALL dbcsr_get_block_p(matrix=tb%quadbra(4)%matrix, &
1366 row=irow, col=icol, block=quad_bra4, found=found)
1367 cpassert(found)
1368 CALL dbcsr_get_block_p(matrix=tb%quadbra(5)%matrix, &
1369 row=irow, col=icol, block=quad_bra5, found=found)
1370 cpassert(found)
1371 CALL dbcsr_get_block_p(matrix=tb%quadbra(6)%matrix, &
1372 row=irow, col=icol, block=quad_bra6, found=found)
1373 cpassert(found)
1374
1375 NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
1376 CALL dbcsr_get_block_p(matrix=tb%quadket(1)%matrix, &
1377 row=irow, col=icol, block=quad_ket1, found=found)
1378 cpassert(found)
1379 CALL dbcsr_get_block_p(matrix=tb%quadket(2)%matrix, &
1380 row=irow, col=icol, block=quad_ket2, found=found)
1381 cpassert(found)
1382 CALL dbcsr_get_block_p(matrix=tb%quadket(3)%matrix, &
1383 row=irow, col=icol, block=quad_ket3, found=found)
1384 cpassert(found)
1385 CALL dbcsr_get_block_p(matrix=tb%quadket(4)%matrix, &
1386 row=irow, col=icol, block=quad_ket4, found=found)
1387 cpassert(found)
1388 CALL dbcsr_get_block_p(matrix=tb%quadket(5)%matrix, &
1389 row=irow, col=icol, block=quad_ket5, found=found)
1390 cpassert(found)
1391 CALL dbcsr_get_block_p(matrix=tb%quadket(6)%matrix, &
1392 row=irow, col=icol, block=quad_ket6, found=found)
1393 cpassert(found)
1394
1395 DO is = 1, SIZE(ks_matrix, 1)
1396 NULLIFY (ksblock)
1397 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
1398 row=irow, col=icol, block=ksblock, found=found)
1399 cpassert(found)
1400
1401 ksblock = ksblock - 0.5_dp*(quad_ket1*tb%pot%vqp(1, irow, 1) &
1402 + quad_ket2*tb%pot%vqp(2, irow, 1) &
1403 + quad_ket3*tb%pot%vqp(3, irow, 1) &
1404 + quad_ket4*tb%pot%vqp(4, irow, 1) &
1405 + quad_ket5*tb%pot%vqp(5, irow, 1) &
1406 + quad_ket6*tb%pot%vqp(6, irow, 1) &
1407 + quad_bra1*tb%pot%vqp(1, icol, 1) &
1408 + quad_bra2*tb%pot%vqp(2, icol, 1) &
1409 + quad_bra3*tb%pot%vqp(3, icol, 1) &
1410 + quad_bra4*tb%pot%vqp(4, icol, 1) &
1411 + quad_bra5*tb%pot%vqp(5, icol, 1) &
1412 + quad_bra6*tb%pot%vqp(6, icol, 1))
1413 END DO
1414 END DO
1415 CALL neighbor_list_iterator_release(nl_iterator)
1416 END IF
1417
1418 END IF
1419
1420#else
1421 mark_used(qs_env)
1422 mark_used(tb)
1423 mark_used(dft_control)
1424 mark_used(calculate_forces)
1425 cpabort("Built without TBLITE")
1426#endif
1427
1428 END SUBROUTINE tb_ham_add_coulomb
1429
1430! **************************************************************************************************
1431!> \brief ...
1432!> \param qs_env ...
1433!> \param tb ...
1434! **************************************************************************************************
1435 SUBROUTINE tb_get_multipole(qs_env, tb)
1436
1437 TYPE(qs_environment_type), POINTER :: qs_env
1438 TYPE(tblite_type), POINTER :: tb
1439
1440#if defined(__TBLITE)
1441
1442 CHARACTER(LEN=*), PARAMETER :: routinen = 'tb_get_multipole'
1443
1444 INTEGER :: ikind, jkind, iatom, jatom, icol, irow, iset, jset, ityp, jtyp
1445 INTEGER :: nkind, natom, handle, nimg, i, inda, indb
1446 INTEGER :: atom_a, atom_b, nseta, nsetb, ia, ib, ij
1447 LOGICAL :: found
1448 REAL(kind=dp) :: r2
1449 REAL(kind=dp), DIMENSION(3) :: rij
1450 INTEGER, DIMENSION(:), POINTER :: la_max, lb_max
1451 INTEGER, DIMENSION(:), POINTER :: nsgfa, nsgfb
1452 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1453 INTEGER, ALLOCATABLE :: atom_of_kind(:)
1454 REAL(kind=dp), ALLOCATABLE :: stmp(:)
1455 REAL(kind=dp), ALLOCATABLE :: dtmp(:, :), qtmp(:, :), dtmpj(:, :), qtmpj(:, :)
1456 REAL(kind=dp), DIMENSION(:, :), POINTER :: dip_ket1, dip_ket2, dip_ket3, &
1457 dip_bra1, dip_bra2, dip_bra3
1458 REAL(kind=dp), DIMENSION(:, :), POINTER :: quad_ket1, quad_ket2, quad_ket3, &
1459 quad_ket4, quad_ket5, quad_ket6, &
1460 quad_bra1, quad_bra2, quad_bra3, &
1461 quad_bra4, quad_bra5, quad_bra6
1462
1463 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1464 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1465 TYPE(dft_control_type), POINTER :: dft_control
1466 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1467 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1468 TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
1469 TYPE(neighbor_list_iterator_p_type), &
1470 DIMENSION(:), POINTER :: nl_iterator
1471 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1472 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1473
1474 CALL timeset(routinen, handle)
1475
1476 !get info from environment vaiarable
1477 NULLIFY (atomic_kind_set, qs_kind_set, sab_orb, particle_set)
1478 NULLIFY (dft_control, matrix_s)
1479 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
1480 qs_kind_set=qs_kind_set, &
1481 sab_orb=sab_orb, &
1482 particle_set=particle_set, &
1483 dft_control=dft_control, &
1484 matrix_s_kp=matrix_s)
1485 natom = SIZE(particle_set)
1486 nkind = SIZE(atomic_kind_set)
1487 nimg = dft_control%nimages
1488
1489 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
1490
1491 !set up basis set lists
1492 ALLOCATE (basis_set_list(nkind))
1493 CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
1494
1495 ALLOCATE (stmp(msao(tb%calc%bas%maxl)**2))
1496 ALLOCATE (dtmp(dip_n, msao(tb%calc%bas%maxl)**2))
1497 ALLOCATE (qtmp(quad_n, msao(tb%calc%bas%maxl)**2))
1498 ALLOCATE (dtmpj(dip_n, msao(tb%calc%bas%maxl)**2))
1499 ALLOCATE (qtmpj(quad_n, msao(tb%calc%bas%maxl)**2))
1500
1501 ! allocate dipole/quadrupole moment matrix elemnts
1502 CALL dbcsr_allocate_matrix_set(tb%dipbra, dip_n)
1503 CALL dbcsr_allocate_matrix_set(tb%dipket, dip_n)
1504 CALL dbcsr_allocate_matrix_set(tb%quadbra, quad_n)
1505 CALL dbcsr_allocate_matrix_set(tb%quadket, quad_n)
1506 DO i = 1, dip_n
1507 ALLOCATE (tb%dipbra(i)%matrix)
1508 ALLOCATE (tb%dipket(i)%matrix)
1509 CALL dbcsr_create(tb%dipbra(i)%matrix, template=matrix_s(1, 1)%matrix, &
1510 name="DIPOLE BRAMATRIX")
1511 CALL dbcsr_create(tb%dipket(i)%matrix, template=matrix_s(1, 1)%matrix, &
1512 name="DIPOLE KETMATRIX")
1513 CALL cp_dbcsr_alloc_block_from_nbl(tb%dipbra(i)%matrix, sab_orb)
1514 CALL cp_dbcsr_alloc_block_from_nbl(tb%dipket(i)%matrix, sab_orb)
1515 END DO
1516 DO i = 1, quad_n
1517 ALLOCATE (tb%quadbra(i)%matrix)
1518 ALLOCATE (tb%quadket(i)%matrix)
1519 CALL dbcsr_create(tb%quadbra(i)%matrix, template=matrix_s(1, 1)%matrix, &
1520 name="QUADRUPOLE BRAMATRIX")
1521 CALL dbcsr_create(tb%quadket(i)%matrix, template=matrix_s(1, 1)%matrix, &
1522 name="QUADRUPOLE KETMATRIX")
1523 CALL cp_dbcsr_alloc_block_from_nbl(tb%quadbra(i)%matrix, sab_orb)
1524 CALL cp_dbcsr_alloc_block_from_nbl(tb%quadket(i)%matrix, sab_orb)
1525 END DO
1526
1527 !loop over all atom pairs with a non-zero overlap (sab_orb)
1528 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1529 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1530 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
1531 iatom=iatom, jatom=jatom, r=rij)
1532
1533 r2 = norm2(rij(:))**2
1534
1535 icol = max(iatom, jatom)
1536 irow = min(iatom, jatom)
1537
1538 IF (icol == jatom) THEN
1539 rij = -rij
1540 i = ikind
1541 ikind = jkind
1542 jkind = i
1543 END IF
1544
1545 ityp = tb%mol%id(icol)
1546 jtyp = tb%mol%id(irow)
1547
1548 NULLIFY (dip_bra1, dip_bra2, dip_bra3)
1549 CALL dbcsr_get_block_p(matrix=tb%dipbra(1)%matrix, &
1550 row=irow, col=icol, block=dip_bra1, found=found)
1551 cpassert(found)
1552 CALL dbcsr_get_block_p(matrix=tb%dipbra(2)%matrix, &
1553 row=irow, col=icol, block=dip_bra2, found=found)
1554 cpassert(found)
1555 CALL dbcsr_get_block_p(matrix=tb%dipbra(3)%matrix, &
1556 row=irow, col=icol, block=dip_bra3, found=found)
1557 cpassert(found)
1558
1559 NULLIFY (dip_ket1, dip_ket2, dip_ket3)
1560 CALL dbcsr_get_block_p(matrix=tb%dipket(1)%matrix, &
1561 row=irow, col=icol, block=dip_ket1, found=found)
1562 cpassert(found)
1563 CALL dbcsr_get_block_p(matrix=tb%dipket(2)%matrix, &
1564 row=irow, col=icol, block=dip_ket2, found=found)
1565 cpassert(found)
1566 CALL dbcsr_get_block_p(matrix=tb%dipket(3)%matrix, &
1567 row=irow, col=icol, block=dip_ket3, found=found)
1568 cpassert(found)
1569
1570 NULLIFY (quad_bra1, quad_bra2, quad_bra3, quad_bra4, quad_bra5, quad_bra6)
1571 CALL dbcsr_get_block_p(matrix=tb%quadbra(1)%matrix, &
1572 row=irow, col=icol, block=quad_bra1, found=found)
1573 cpassert(found)
1574 CALL dbcsr_get_block_p(matrix=tb%quadbra(2)%matrix, &
1575 row=irow, col=icol, block=quad_bra2, found=found)
1576 cpassert(found)
1577 CALL dbcsr_get_block_p(matrix=tb%quadbra(3)%matrix, &
1578 row=irow, col=icol, block=quad_bra3, found=found)
1579 cpassert(found)
1580 CALL dbcsr_get_block_p(matrix=tb%quadbra(4)%matrix, &
1581 row=irow, col=icol, block=quad_bra4, found=found)
1582 cpassert(found)
1583 CALL dbcsr_get_block_p(matrix=tb%quadbra(5)%matrix, &
1584 row=irow, col=icol, block=quad_bra5, found=found)
1585 cpassert(found)
1586 CALL dbcsr_get_block_p(matrix=tb%quadbra(6)%matrix, &
1587 row=irow, col=icol, block=quad_bra6, found=found)
1588 cpassert(found)
1589
1590 NULLIFY (quad_ket1, quad_ket2, quad_ket3, quad_ket4, quad_ket5, quad_ket6)
1591 CALL dbcsr_get_block_p(matrix=tb%quadket(1)%matrix, &
1592 row=irow, col=icol, block=quad_ket1, found=found)
1593 cpassert(found)
1594 CALL dbcsr_get_block_p(matrix=tb%quadket(2)%matrix, &
1595 row=irow, col=icol, block=quad_ket2, found=found)
1596 cpassert(found)
1597 CALL dbcsr_get_block_p(matrix=tb%quadket(3)%matrix, &
1598 row=irow, col=icol, block=quad_ket3, found=found)
1599 cpassert(found)
1600 CALL dbcsr_get_block_p(matrix=tb%quadket(4)%matrix, &
1601 row=irow, col=icol, block=quad_ket4, found=found)
1602 cpassert(found)
1603 CALL dbcsr_get_block_p(matrix=tb%quadket(5)%matrix, &
1604 row=irow, col=icol, block=quad_ket5, found=found)
1605 cpassert(found)
1606 CALL dbcsr_get_block_p(matrix=tb%quadket(6)%matrix, &
1607 row=irow, col=icol, block=quad_ket6, found=found)
1608 cpassert(found)
1609
1610 !get basis information
1611 basis_set_a => basis_set_list(ikind)%gto_basis_set
1612 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1613 basis_set_b => basis_set_list(jkind)%gto_basis_set
1614 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1615 atom_a = atom_of_kind(icol)
1616 atom_b = atom_of_kind(irow)
1617 ! basis a
1618 first_sgfa => basis_set_a%first_sgf
1619 la_max => basis_set_a%lmax
1620 nseta = basis_set_a%nset
1621 nsgfa => basis_set_a%nsgf_set
1622 ! basis b
1623 first_sgfb => basis_set_b%first_sgf
1624 lb_max => basis_set_b%lmax
1625 nsetb = basis_set_b%nset
1626 nsgfb => basis_set_b%nsgf_set
1627
1628 ! --------- Hamiltonian
1629 IF (icol == irow) THEN
1630 DO iset = 1, nseta
1631 DO jset = 1, nsetb
1632 CALL multipole_cgto(tb%calc%bas%cgto(jset, ityp), tb%calc%bas%cgto(iset, ityp), &
1633 & r2, -rij, tb%calc%bas%intcut, stmp, dtmp, qtmp)
1634
1635 DO inda = 1, nsgfa(iset)
1636 ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
1637 DO indb = 1, nsgfb(jset)
1638 ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
1639 ij = indb + nsgfb(jset)*(inda - 1)
1640
1641 dip_ket1(ib, ia) = dtmp(1, ij)
1642 dip_ket2(ib, ia) = dtmp(2, ij)
1643 dip_ket3(ib, ia) = dtmp(3, ij)
1644
1645 quad_ket1(ib, ia) = qtmp(1, ij)
1646 quad_ket2(ib, ia) = qtmp(2, ij)
1647 quad_ket3(ib, ia) = qtmp(3, ij)
1648 quad_ket4(ib, ia) = qtmp(4, ij)
1649 quad_ket5(ib, ia) = qtmp(5, ij)
1650 quad_ket6(ib, ia) = qtmp(6, ij)
1651
1652 dip_bra1(ib, ia) = dtmp(1, ij)
1653 dip_bra2(ib, ia) = dtmp(2, ij)
1654 dip_bra3(ib, ia) = dtmp(3, ij)
1655
1656 quad_bra1(ib, ia) = qtmp(1, ij)
1657 quad_bra2(ib, ia) = qtmp(2, ij)
1658 quad_bra3(ib, ia) = qtmp(3, ij)
1659 quad_bra4(ib, ia) = qtmp(4, ij)
1660 quad_bra5(ib, ia) = qtmp(5, ij)
1661 quad_bra6(ib, ia) = qtmp(6, ij)
1662 END DO
1663 END DO
1664 END DO
1665 END DO
1666 ELSE
1667 DO iset = 1, nseta
1668 DO jset = 1, nsetb
1669 CALL multipole_cgto(tb%calc%bas%cgto(jset, jtyp), tb%calc%bas%cgto(iset, ityp), &
1670 & r2, -rij, tb%calc%bas%intcut, stmp, dtmp, qtmp)
1671 CALL multipole_cgto(tb%calc%bas%cgto(iset, ityp), tb%calc%bas%cgto(jset, jtyp), &
1672 & r2, rij, tb%calc%bas%intcut, stmp, dtmpj, qtmpj)
1673
1674 DO inda = 1, nsgfa(iset)
1675 ia = first_sgfa(1, iset) - first_sgfa(1, 1) + inda
1676 DO indb = 1, nsgfb(jset)
1677 ib = first_sgfb(1, jset) - first_sgfb(1, 1) + indb
1678
1679 ij = indb + nsgfb(jset)*(inda - 1)
1680
1681 dip_bra1(ib, ia) = dtmp(1, ij)
1682 dip_bra2(ib, ia) = dtmp(2, ij)
1683 dip_bra3(ib, ia) = dtmp(3, ij)
1684
1685 quad_bra1(ib, ia) = qtmp(1, ij)
1686 quad_bra2(ib, ia) = qtmp(2, ij)
1687 quad_bra3(ib, ia) = qtmp(3, ij)
1688 quad_bra4(ib, ia) = qtmp(4, ij)
1689 quad_bra5(ib, ia) = qtmp(5, ij)
1690 quad_bra6(ib, ia) = qtmp(6, ij)
1691
1692 ij = inda + nsgfa(iset)*(indb - 1)
1693
1694 dip_ket1(ib, ia) = dtmpj(1, ij)
1695 dip_ket2(ib, ia) = dtmpj(2, ij)
1696 dip_ket3(ib, ia) = dtmpj(3, ij)
1697
1698 quad_ket1(ib, ia) = qtmpj(1, ij)
1699 quad_ket2(ib, ia) = qtmpj(2, ij)
1700 quad_ket3(ib, ia) = qtmpj(3, ij)
1701 quad_ket4(ib, ia) = qtmpj(4, ij)
1702 quad_ket5(ib, ia) = qtmpj(5, ij)
1703 quad_ket6(ib, ia) = qtmpj(6, ij)
1704 END DO
1705 END DO
1706 END DO
1707 END DO
1708 END IF
1709 END DO
1710 CALL neighbor_list_iterator_release(nl_iterator)
1711
1712 DO i = 1, dip_n
1713 CALL dbcsr_finalize(tb%dipbra(i)%matrix)
1714 CALL dbcsr_finalize(tb%dipket(i)%matrix)
1715 END DO
1716 DO i = 1, quad_n
1717 CALL dbcsr_finalize(tb%quadbra(i)%matrix)
1718 CALL dbcsr_finalize(tb%quadket(i)%matrix)
1719 END DO
1720
1721 DEALLOCATE (basis_set_list)
1722
1723 CALL timestop(handle)
1724
1725#else
1726 mark_used(qs_env)
1727 mark_used(tb)
1728 cpabort("Built without TBLITE")
1729#endif
1730
1731 END SUBROUTINE tb_get_multipole
1732
1733! **************************************************************************************************
1734!> \brief compute the mulliken properties (AO resolved)
1735!> \param p_matrix ...
1736!> \param bra_mat ...
1737!> \param ket_mat ...
1738!> \param output ...
1739!> \param para_env ...
1740!> \par History
1741!> adapted from ao_charges_2
1742!> \note
1743! **************************************************************************************************
1744 SUBROUTINE contract_dens(p_matrix, bra_mat, ket_mat, output, para_env)
1745 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
1746 TYPE(dbcsr_type), POINTER :: bra_mat, ket_mat
1747 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: output
1748 TYPE(mp_para_env_type), POINTER :: para_env
1749
1750 CHARACTER(len=*), PARAMETER :: routinen = 'contract_dens'
1751
1752 INTEGER :: handle, i, iblock_col, iblock_row, &
1753 ispin, j, nspin
1754 LOGICAL :: found
1755 REAL(kind=dp), DIMENSION(:, :), POINTER :: bra, ket, p_block
1756 TYPE(dbcsr_iterator_type) :: iter
1757
1758 CALL timeset(routinen, handle)
1759
1760 nspin = SIZE(p_matrix)
1761 output = 0.0_dp
1762 DO ispin = 1, nspin
1763 CALL dbcsr_iterator_start(iter, bra_mat)
1764 DO WHILE (dbcsr_iterator_blocks_left(iter))
1765 NULLIFY (p_block, bra, ket)
1766
1767 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, bra)
1768 CALL dbcsr_get_block_p(matrix=p_matrix(ispin)%matrix, &
1769 row=iblock_row, col=iblock_col, block=p_block, found=found)
1770 IF (.NOT. found) cycle
1771 CALL dbcsr_get_block_p(matrix=ket_mat, &
1772 row=iblock_row, col=iblock_col, block=ket, found=found)
1773 IF (.NOT. found) cpabort("missing block")
1774
1775 IF (.NOT. (ASSOCIATED(bra) .AND. ASSOCIATED(p_block))) cycle
1776 DO j = 1, SIZE(p_block, 1)
1777 DO i = 1, SIZE(p_block, 2)
1778 output(iblock_row) = output(iblock_row) + p_block(j, i)*ket(j, i)
1779 END DO
1780 END DO
1781 IF (iblock_col /= iblock_row) THEN
1782 DO j = 1, SIZE(p_block, 1)
1783 DO i = 1, SIZE(p_block, 2)
1784 output(iblock_col) = output(iblock_col) + p_block(j, i)*bra(j, i)
1785 END DO
1786 END DO
1787 END IF
1788 END DO
1789 CALL dbcsr_iterator_stop(iter)
1790 END DO
1791
1792 CALL para_env%sum(output)
1793 CALL timestop(handle)
1794
1795 END SUBROUTINE contract_dens
1796
1797END MODULE tblite_interface
1798
ai_contraction::block_add
Definition ai_contraction.F:43
ai_contraction::contraction
Definition ai_contraction.F:31
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
memory_utilities::reallocate
Definition memory_utilities.F:27
mulliken::ao_charges
Definition mulliken.F:42
qs_integral_utils::get_memory_usage
Definition qs_integral_utils.F:34
qs_overlap::create_sab_matrix
Definition qs_overlap.F:81
ai_contraction
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Definition ai_contraction.F:18
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
ai_overlap::overlap_ab
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:680
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
basis_set_types
Definition basis_set_types.F:15
basis_set_types::process_gto_basis
subroutine, public process_gto_basis(gto_basis_set, do_ortho, nset, maxl)
...
Definition basis_set_types.F:2350
basis_set_types::allocate_gto_basis_set
subroutine, public allocate_gto_basis_set(gto_basis_set)
...
Definition basis_set_types.F:151
basis_set_types::write_gto_basis_set
subroutine, public write_gto_basis_set(gto_basis_set, output_unit, header)
Write a Gaussian-type orbital (GTO) basis set data set to the output unit.
Definition basis_set_types.F:1860
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:195
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
Definition cp_dbcsr_api.F:956
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:930
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1027
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:679
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:659
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:989
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_print
subroutine, public dbcsr_print(matrix, variable_name, unit_nr)
Prints given matrix in matlab format (only present blocks).
Definition cp_dbcsr_contrib.F:682
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::gfn1xtb
integer, parameter, public gfn1xtb
Definition input_constants.F:1341
input_constants::ipea1xtb
integer, parameter, public ipea1xtb
Definition input_constants.F:1341
input_constants::gfn2xtb
integer, parameter, public gfn2xtb
Definition input_constants.F:1341
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:365
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
mulliken
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition mulliken.F:13
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_charge_mixing
Definition qs_charge_mixing.F:9
qs_charge_mixing::charge_mixing
subroutine, public charge_mixing(mixing_method, mixing_store, charges, para_env, iter_count)
Driver for the charge mixing, calls the proper routine given the requested method.
Definition qs_charge_mixing.F:42
qs_condnum
Calculation of overlap matrix condition numbers.
Definition qs_condnum.F:13
qs_condnum::overlap_condnum
subroutine, public overlap_condnum(matrixkp_s, condnum, iunit, norml1, norml2, use_arnoldi, blacs_env)
Calculation of the overlap matrix Condition Number.
Definition qs_condnum.F:66
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:518
qs_force_types
Definition qs_force_types.F:13
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:879
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:572
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:549
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::integer_to_string
subroutine, public integer_to_string(inumber, string)
Converts an integer number to a string. The WRITE statement will return an error message,...
Definition string_utilities.F:3239
tblite_interface
interface to tblite
Definition tblite_interface.F:14
tblite_interface::tb_set_calculator
subroutine, public tb_set_calculator(tb, typ)
...
Definition tblite_interface.F:223
tblite_interface::tb_init_ham
subroutine, public tb_init_ham(tb)
...
Definition tblite_interface.F:255
tblite_interface::tb_update_charges
subroutine, public tb_update_charges(qs_env, dft_control, tb, calculate_forces, use_rho)
...
Definition tblite_interface.F:841
tblite_interface::tb_init_geometry
subroutine, public tb_init_geometry(qs_env, tb)
...
Definition tblite_interface.F:109
tblite_interface::tb_ham_add_coulomb
subroutine, public tb_ham_add_coulomb(qs_env, tb, dft_control, calculate_forces)
...
Definition tblite_interface.F:1026
tblite_interface::build_tblite_matrices
subroutine, public build_tblite_matrices(qs_env, calculate_forces)
...
Definition tblite_interface.F:512
tblite_interface::tb_get_energy
subroutine, public tb_get_energy(qs_env, tb, energy)
...
Definition tblite_interface.F:333
tblite_interface::tb_get_multipole
subroutine, public tb_get_multipole(qs_env, tb)
...
Definition tblite_interface.F:1436
tblite_interface::tb_init_wf
subroutine, public tb_init_wf(tb, do_grad)
...
Definition tblite_interface.F:173
tblite_interface::tb_get_basis
subroutine, public tb_get_basis(tb, gto_basis_set, element_symbol, param, occ)
...
Definition tblite_interface.F:402
tblite_types
types for tblite
Definition tblite_types.F:14
tblite_types::allocate_tblite_type
subroutine, public allocate_tblite_type(tb_tblite)
...
Definition tblite_types.F:78
tblite_types::deallocate_tblite_type
subroutine, public deallocate_tblite_type(tb_tblite)
...
Definition tblite_types.F:92
virial_types
Definition virial_types.F:13
xtb_types
Definition of the xTB parameter types.
Definition xtb_types.F:20
xtb_types::get_xtb_atom_param
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, xi, kappa0, alpg, occupation, electronegativity, chmax, en, kqat2, kcn, kq)
...
Definition xtb_types.F:199
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:586
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:186
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:150
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:219
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:136
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:97
tblite_types::tblite_type
Definition tblite_types.F:37
virial_types::virial_type
Definition virial_types.F:26
xtb_types::xtb_atom_type
Definition xtb_types.F:42