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