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