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