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