(git:591cf04)
Loading...
Searching...
No Matches
xtb_matrices.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculation of Overlap and Hamiltonian matrices in xTB
10!> Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov
11!> JCTC 13, 1989-2009, (2017)
12!> DOI: 10.1021/acs.jctc.7b00118
13!> \author JGH
14! **************************************************************************************************
15MODULE xtb_matrices
16 USE ai_contraction, ONLY: block_add,&
17 contraction
18 USE ai_overlap, ONLY: overlap_ab
19 USE atomic_kind_types, ONLY: atomic_kind_type,&
20 get_atomic_kind_set
21 USE atprop_types, ONLY: atprop_array_init,&
22 atprop_type
23 USE basis_set_types, ONLY: gto_basis_set_p_type,&
24 gto_basis_set_type
25 USE block_p_types, ONLY: block_p_type
26 USE cp_blacs_env, ONLY: cp_blacs_env_type
27 USE cp_control_types, ONLY: dft_control_type,&
28 xtb_control_type
29 USE cp_dbcsr_api, ONLY: dbcsr_add,&
30 dbcsr_create,&
31 dbcsr_finalize,&
32 dbcsr_get_block_p,&
33 dbcsr_p_type
34 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
35 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
36 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
37 USE cp_log_handling, ONLY: cp_get_default_logger,&
38 cp_logger_type
39 USE cp_output_handling, ONLY: cp_p_file,&
40 cp_print_key_finished_output,&
41 cp_print_key_should_output,&
42 cp_print_key_unit_nr
43 USE eeq_input, ONLY: eeq_solver_type
44 USE input_constants, ONLY: vdw_pairpot_dftd4
45 USE input_section_types, ONLY: section_vals_val_get
46 USE kinds, ONLY: dp
47 USE kpoint_types, ONLY: get_kpoint_info,&
48 kpoint_type
49 USE message_passing, ONLY: mp_para_env_type
50 USE orbital_pointers, ONLY: ncoset
51 USE particle_types, ONLY: particle_type
52 USE qs_condnum, ONLY: overlap_condnum
53 USE qs_dispersion_cnum, ONLY: cnumber_init,&
54 cnumber_release,&
55 dcnum_type
56 USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
57 USE qs_dispersion_types, ONLY: qs_dispersion_type
58 USE qs_energy_types, ONLY: qs_energy_type
59 USE qs_environment_types, ONLY: get_qs_env,&
60 qs_environment_type,&
61 set_qs_env
62 USE qs_force_types, ONLY: qs_force_type
63 USE qs_integral_utils, ONLY: basis_set_list_setup,&
64 get_memory_usage
65 USE qs_kind_types, ONLY: get_qs_kind,&
66 qs_kind_type
67 USE qs_ks_types, ONLY: get_ks_env,&
68 qs_ks_env_type,&
69 set_ks_env
70 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
71 neighbor_list_iterate,&
72 neighbor_list_iterator_create,&
73 neighbor_list_iterator_p_type,&
74 neighbor_list_iterator_release,&
75 neighbor_list_set_p_type
76 USE qs_overlap, ONLY: create_sab_matrix
77 USE qs_rho_types, ONLY: qs_rho_get,&
78 qs_rho_type
79 USE virial_methods, ONLY: virial_pair_force
80 USE virial_types, ONLY: virial_type
81 USE xtb_eeq, ONLY: xtb_eeq_calculation,&
82 xtb_eeq_forces
83 USE xtb_hcore, ONLY: gfn0_huckel,&
84 gfn0_kpair,&
85 gfn1_huckel,&
86 gfn1_kpair
87 USE xtb_potentials, ONLY: nonbonded_correction,&
88 repulsive_potential,&
89 srb_potential,&
90 xb_interaction
91 USE xtb_types, ONLY: get_xtb_atom_param,&
92 xtb_atom_type
93#include "./base/base_uses.f90"
94
95 IMPLICIT NONE
96
97 PRIVATE
98
99 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_matrices'
100
101 PUBLIC :: build_xtb_matrices
102
103CONTAINS
104
105! **************************************************************************************************
106!> \brief ...
107!> \param qs_env ...
108!> \param calculate_forces ...
109! **************************************************************************************************
110 SUBROUTINE build_xtb_matrices(qs_env, calculate_forces)
111
112 TYPE(qs_environment_type), POINTER :: qs_env
113 LOGICAL, INTENT(IN) :: calculate_forces
114
115 INTEGER :: gfn_type
116 TYPE(dft_control_type), POINTER :: dft_control
117
118 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
119 gfn_type = dft_control%qs_control%xtb_control%gfn_type
120
121 SELECT CASE (gfn_type)
122 CASE (0)
123 CALL build_gfn0_xtb_matrices(qs_env, calculate_forces)
124 CASE (1)
125 CALL build_gfn1_xtb_matrices(qs_env, calculate_forces)
126 CASE (2)
127 cpabort("gfn_type = 2 not yet available")
128 CASE DEFAULT
129 cpabort("Unknown gfn_type")
130 END SELECT
131
132 END SUBROUTINE build_xtb_matrices
133
134! **************************************************************************************************
135!> \brief ...
136!> \param qs_env ...
137!> \param calculate_forces ...
138! **************************************************************************************************
139 SUBROUTINE build_gfn0_xtb_matrices(qs_env, calculate_forces)
140
141 TYPE(qs_environment_type), POINTER :: qs_env
142 LOGICAL, INTENT(IN) :: calculate_forces
143
144 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_gfn0_xtb_matrices'
145
146 INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, iset, &
147 j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, lmaxa, lmaxb, maxder, n1, n2, na, &
148 natom, natorb_a, natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, &
149 nsetb, sgfa, sgfb, za, zb
150 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
151 INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
152 INTEGER, DIMENSION(3) :: cell
153 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
154 npgfb, nsgfa, nsgfb
155 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
156 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
157 LOGICAL :: defined, diagblock, do_nonbonded, found, &
158 use_virial
159 REAL(kind=dp) :: dfp, dhij, dr, drk, drx, eeq_energy, ef_energy, enonbonded, enscale, erep, &
160 esrb, etaa, etab, f0, f1, f2, fhua, fhub, fhud, foab, fqa, fqb, hij, kf, qlambda, rcova, &
161 rcovab, rcovb, rrab
162 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: charges, cnumbers, dcharges, qlagrange
163 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, dqhuckel, huckel, owork
164 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
165 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
166 REAL(kind=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_ab, rij, rik
167 REAL(kind=dp), DIMENSION(5) :: dpia, dpib, hena, henb, kpolya, kpolyb, &
168 pia, pib
169 REAL(kind=dp), DIMENSION(:), POINTER :: eeq_q, set_radius_a, set_radius_b
170 REAL(kind=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
171 scon_a, scon_b, wblock, zeta, zetb
172 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
173 TYPE(atprop_type), POINTER :: atprop
174 TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
175 TYPE(cp_logger_type), POINTER :: logger
176 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
177 TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
178 TYPE(dft_control_type), POINTER :: dft_control
179 TYPE(eeq_solver_type) :: eeq_sparam
180 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
181 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
182 TYPE(kpoint_type), POINTER :: kpoints
183 TYPE(mp_para_env_type), POINTER :: para_env
184 TYPE(neighbor_list_iterator_p_type), &
185 DIMENSION(:), POINTER :: nl_iterator
186 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
187 POINTER :: sab_orb, sab_xtb_nonbond
188 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
189 TYPE(qs_dispersion_type), POINTER :: dispersion_env
190 TYPE(qs_energy_type), POINTER :: energy
191 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
192 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
193 TYPE(qs_ks_env_type), POINTER :: ks_env
194 TYPE(qs_rho_type), POINTER :: rho
195 TYPE(virial_type), POINTER :: virial
196 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
197 TYPE(xtb_control_type), POINTER :: xtb_control
198
199 CALL timeset(routinen, handle)
200
201 NULLIFY (logger, virial, atprop)
202 logger => cp_get_default_logger()
203
204 NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
205 qs_kind_set, sab_orb, ks_env)
206 CALL get_qs_env(qs_env=qs_env, &
207 ks_env=ks_env, &
208 energy=energy, &
209 atomic_kind_set=atomic_kind_set, &
210 qs_kind_set=qs_kind_set, &
211 matrix_h_kp=matrix_h, &
212 matrix_s_kp=matrix_s, &
213 para_env=para_env, &
214 atprop=atprop, &
215 dft_control=dft_control, &
216 sab_orb=sab_orb)
217
218 nkind = SIZE(atomic_kind_set)
219 xtb_control => dft_control%qs_control%xtb_control
220 eeq_sparam = xtb_control%eeq_sparam
221 do_nonbonded = xtb_control%do_nonbonded
222 nimg = dft_control%nimages
223 nderivatives = 0
224 IF (calculate_forces) nderivatives = 1
225 IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
226 maxder = ncoset(nderivatives)
227
228 NULLIFY (particle_set)
229 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
230 natom = SIZE(particle_set)
231 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
232 atom_of_kind=atom_of_kind, kind_of=kind_of)
233
234 IF (calculate_forces) THEN
235 NULLIFY (rho, force, matrix_w)
236 CALL get_qs_env(qs_env=qs_env, &
237 rho=rho, matrix_w_kp=matrix_w, &
238 virial=virial, force=force)
239 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
240
241 IF (SIZE(matrix_p, 1) == 2) THEN
242 DO img = 1, nimg
243 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
244 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
245 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
246 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
247 END DO
248 END IF
249 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
250 END IF
251 ! atomic energy decomposition
252 IF (atprop%energy) THEN
253 CALL atprop_array_init(atprop%atecc, natom)
254 END IF
255
256 NULLIFY (cell_to_index)
257 IF (nimg > 1) THEN
258 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
259 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
260 END IF
261
262 ! set up basis set lists
263 ALLOCATE (basis_set_list(nkind))
264 CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
265
266 ! allocate overlap matrix
267 CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
268 CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
269 sab_orb, .true.)
270 CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
271
272 ! initialize H matrix
273 CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
274 DO img = 1, nimg
275 ALLOCATE (matrix_h(1, img)%matrix)
276 CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
277 name="HAMILTONIAN MATRIX")
278 CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
279 END DO
280 CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
281
282 ! Calculate coordination numbers
283 ! needed for effective atomic energy levels
284 ! code taken from D3 dispersion energy
285 CALL cnumber_init(qs_env, cnumbers, dcnum, 2, calculate_forces)
286
287 ALLOCATE (charges(natom))
288 charges = 0.0_dp
289 CALL xtb_eeq_calculation(qs_env, charges, cnumbers, eeq_sparam, eeq_energy, ef_energy, qlambda)
290 IF (calculate_forces) THEN
291 ALLOCATE (dcharges(natom))
292 dcharges = qlambda/real(para_env%num_pe, kind=dp)
293 END IF
294 energy%eeq = eeq_energy
295 energy%efield = ef_energy
296
297 CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
298 ! prepare charges (needed for D4)
299 IF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
300 dispersion_env%ext_charges = .true.
301 IF (ASSOCIATED(dispersion_env%charges)) DEALLOCATE (dispersion_env%charges)
302 ALLOCATE (dispersion_env%charges(natom))
303 dispersion_env%charges = charges
304 IF (calculate_forces) THEN
305 IF (ASSOCIATED(dispersion_env%dcharges)) DEALLOCATE (dispersion_env%dcharges)
306 ALLOCATE (dispersion_env%dcharges(natom))
307 dispersion_env%dcharges = 0.0_dp
308 END IF
309 END IF
310 CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
311 energy%dispersion, calculate_forces)
312 IF (calculate_forces) THEN
313 IF (dispersion_env%pp_type == vdw_pairpot_dftd4 .AND. dispersion_env%ext_charges) THEN
314 dcharges(1:natom) = dcharges(1:natom) + dispersion_env%dcharges(1:natom)
315 END IF
316 END IF
317
318 ! Calculate Huckel parameters
319 CALL gfn0_huckel(qs_env, cnumbers, charges, huckel, dhuckel, dqhuckel, calculate_forces)
320
321 ! Calculate KAB parameters and electronegativity correction
322 CALL gfn0_kpair(qs_env, kijab)
323
324 ! loop over all atom pairs with a non-zero overlap (sab_orb)
325 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
326 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
327 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
328 iatom=iatom, jatom=jatom, r=rij, cell=cell)
329 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
330 CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
331 IF (.NOT. defined .OR. natorb_a < 1) cycle
332 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
333 CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
334 IF (.NOT. defined .OR. natorb_b < 1) cycle
335
336 dr = sqrt(sum(rij(:)**2))
337
338 ! atomic parameters
339 CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
340 lmax=lmaxa, nshell=nsa, kpoly=kpolya, hen=hena)
341 CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
342 lmax=lmaxb, nshell=nsb, kpoly=kpolyb, hen=henb)
343
344 IF (nimg == 1) THEN
345 ic = 1
346 ELSE
347 ic = cell_to_index(cell(1), cell(2), cell(3))
348 cpassert(ic > 0)
349 END IF
350
351 icol = max(iatom, jatom)
352 irow = min(iatom, jatom)
353 NULLIFY (sblock, fblock)
354 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
355 row=irow, col=icol, block=sblock, found=found)
356 cpassert(found)
357 CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
358 row=irow, col=icol, block=fblock, found=found)
359 cpassert(found)
360
361 IF (calculate_forces) THEN
362 NULLIFY (pblock)
363 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
364 row=irow, col=icol, block=pblock, found=found)
365 cpassert(ASSOCIATED(pblock))
366 NULLIFY (wblock)
367 CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
368 row=irow, col=icol, block=wblock, found=found)
369 cpassert(ASSOCIATED(wblock))
370 DO i = 2, 4
371 NULLIFY (dsblocks(i)%block)
372 CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
373 row=irow, col=icol, block=dsblocks(i)%block, found=found)
374 cpassert(found)
375 END DO
376 END IF
377
378 ! overlap
379 basis_set_a => basis_set_list(ikind)%gto_basis_set
380 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
381 basis_set_b => basis_set_list(jkind)%gto_basis_set
382 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
383 atom_a = atom_of_kind(iatom)
384 atom_b = atom_of_kind(jatom)
385 ! basis ikind
386 first_sgfa => basis_set_a%first_sgf
387 la_max => basis_set_a%lmax
388 la_min => basis_set_a%lmin
389 npgfa => basis_set_a%npgf
390 nseta = basis_set_a%nset
391 nsgfa => basis_set_a%nsgf_set
392 rpgfa => basis_set_a%pgf_radius
393 set_radius_a => basis_set_a%set_radius
394 scon_a => basis_set_a%scon
395 zeta => basis_set_a%zet
396 ! basis jkind
397 first_sgfb => basis_set_b%first_sgf
398 lb_max => basis_set_b%lmax
399 lb_min => basis_set_b%lmin
400 npgfb => basis_set_b%npgf
401 nsetb = basis_set_b%nset
402 nsgfb => basis_set_b%nsgf_set
403 rpgfb => basis_set_b%pgf_radius
404 set_radius_b => basis_set_b%set_radius
405 scon_b => basis_set_b%scon
406 zetb => basis_set_b%zet
407
408 ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
409 ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
410 ALLOCATE (sint(natorb_a, natorb_b, maxder))
411 sint = 0.0_dp
412
413 DO iset = 1, nseta
414 ncoa = npgfa(iset)*ncoset(la_max(iset))
415 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
416 sgfa = first_sgfa(1, iset)
417 DO jset = 1, nsetb
418 IF (set_radius_a(iset) + set_radius_b(jset) < dr) cycle
419 ncob = npgfb(jset)*ncoset(lb_max(jset))
420 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
421 sgfb = first_sgfb(1, jset)
422 IF (calculate_forces) THEN
423 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
424 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
425 rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
426 ELSE
427 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
428 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
429 rij, sab=oint(:, :, 1))
430 END IF
431 ! Contraction
432 CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
433 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
434 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.false.)
435 IF (calculate_forces) THEN
436 DO i = 2, 4
437 CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
438 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
439 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.false.)
440 END DO
441 END IF
442 END DO
443 END DO
444 ! forces W matrix
445 IF (calculate_forces) THEN
446 DO i = 1, 3
447 IF (iatom <= jatom) THEN
448 force_ab(i) = sum(sint(:, :, i + 1)*wblock(:, :))
449 ELSE
450 force_ab(i) = sum(sint(:, :, i + 1)*transpose(wblock(:, :)))
451 END IF
452 END DO
453 f1 = 2.0_dp
454 force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
455 force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
456 IF (use_virial .AND. dr > 1.e-3_dp) THEN
457 IF (iatom == jatom) f1 = 1.0_dp
458 CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
459 END IF
460 END IF
461 ! update S matrix
462 IF (iatom <= jatom) THEN
463 sblock(:, :) = sblock(:, :) + sint(:, :, 1)
464 ELSE
465 sblock(:, :) = sblock(:, :) + transpose(sint(:, :, 1))
466 END IF
467 IF (calculate_forces) THEN
468 DO i = 2, 4
469 IF (iatom <= jatom) THEN
470 dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
471 ELSE
472 dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - transpose(sint(:, :, i))
473 END IF
474 END DO
475 END IF
476
477 ! Calculate Pi = Pia * Pib (Eq. 11)
478 rcovab = rcova + rcovb
479 rrab = sqrt(dr/rcovab)
480 pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
481 pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
482 IF (calculate_forces) THEN
483 IF (dr > 1.e-6_dp) THEN
484 drx = 0.5_dp/rrab/rcovab
485 ELSE
486 drx = 0.0_dp
487 END IF
488 dpia(1:nsa) = drx*kpolya(1:nsa)
489 dpib(1:nsb) = drx*kpolyb(1:nsb)
490 END IF
491
492 ! diagonal block
493 diagblock = .false.
494 IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .true.
495 !
496 ! Eq. 10
497 !
498 IF (diagblock) THEN
499 DO i = 1, natorb_a
500 na = naoa(i)
501 fblock(i, i) = fblock(i, i) + huckel(na, iatom)
502 END DO
503 ELSE
504 DO j = 1, natorb_b
505 nb = naob(j)
506 DO i = 1, natorb_a
507 na = naoa(i)
508 hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
509 IF (iatom <= jatom) THEN
510 fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
511 ELSE
512 fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
513 END IF
514 END DO
515 END DO
516 END IF
517 IF (calculate_forces) THEN
518 f0 = 1.0_dp
519 IF (irow == iatom) f0 = -1.0_dp
520 f2 = 1.0_dp
521 IF (iatom /= jatom) f2 = 2.0_dp
522 ! Derivative wrt coordination number
523 fhua = 0.0_dp
524 fhub = 0.0_dp
525 fhud = 0.0_dp
526 fqa = 0.0_dp
527 fqb = 0.0_dp
528 IF (diagblock) THEN
529 DO i = 1, natorb_a
530 la = laoa(i)
531 na = naoa(i)
532 fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
533 fqa = fqa + pblock(i, i)*dqhuckel(na, iatom)
534 END DO
535 dcharges(iatom) = dcharges(iatom) + fqa
536 ELSE
537 DO j = 1, natorb_b
538 lb = laob(j)
539 nb = naob(j)
540 DO i = 1, natorb_a
541 la = laoa(i)
542 na = naoa(i)
543 hij = 0.5_dp*pia(na)*pib(nb)
544 drx = f2*hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)
545 IF (iatom <= jatom) THEN
546 fhua = fhua + drx*pblock(i, j)*dhuckel(na, iatom)
547 fhub = fhub + drx*pblock(i, j)*dhuckel(nb, jatom)
548 fqa = fqa + drx*pblock(i, j)*dqhuckel(na, iatom)
549 fqb = fqb + drx*pblock(i, j)*dqhuckel(nb, jatom)
550 ELSE
551 fhua = fhua + drx*pblock(j, i)*dhuckel(na, iatom)
552 fhub = fhub + drx*pblock(j, i)*dhuckel(nb, jatom)
553 fqa = fqa + drx*pblock(j, i)*dqhuckel(na, iatom)
554 fqb = fqb + drx*pblock(j, i)*dqhuckel(nb, jatom)
555 END IF
556 END DO
557 END DO
558 dcharges(iatom) = dcharges(iatom) + fqa
559 dcharges(jatom) = dcharges(jatom) + fqb
560 END IF
561 ! iatom
562 atom_a = atom_of_kind(iatom)
563 DO i = 1, dcnum(iatom)%neighbors
564 katom = dcnum(iatom)%nlist(i)
565 kkind = kind_of(katom)
566 atom_c = atom_of_kind(katom)
567 rik = dcnum(iatom)%rik(:, i)
568 drk = sqrt(sum(rik(:)**2))
569 IF (drk > 1.e-3_dp) THEN
570 fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
571 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
572 force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
573 fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
574 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
575 force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
576 IF (use_virial) THEN
577 fdik = fdika + fdikb
578 CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
579 END IF
580 END IF
581 END DO
582 ! jatom
583 atom_b = atom_of_kind(jatom)
584 DO i = 1, dcnum(jatom)%neighbors
585 katom = dcnum(jatom)%nlist(i)
586 kkind = kind_of(katom)
587 atom_c = atom_of_kind(katom)
588 rik = dcnum(jatom)%rik(:, i)
589 drk = sqrt(sum(rik(:)**2))
590 IF (drk > 1.e-3_dp) THEN
591 fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
592 force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
593 force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
594 IF (use_virial) THEN
595 CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
596 END IF
597 END IF
598 END DO
599 ! force from R dendent Huckel element: Pia*Pib
600 IF (diagblock) THEN
601 force_ab = 0._dp
602 ELSE
603 n1 = SIZE(fblock, 1)
604 n2 = SIZE(fblock, 2)
605 ALLOCATE (dfblock(n1, n2))
606 dfblock = 0.0_dp
607 DO j = 1, natorb_b
608 lb = laob(j)
609 nb = naob(j)
610 DO i = 1, natorb_a
611 la = laoa(i)
612 na = naoa(i)
613 dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
614 IF (iatom <= jatom) THEN
615 dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
616 ELSE
617 dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
618 END IF
619 END DO
620 END DO
621 dfp = f0*sum(dfblock(:, :)*pblock(:, :))
622 DO ir = 1, 3
623 foab = 2.0_dp*dfp*rij(ir)/dr
624 ! force from overlap matrix contribution to H
625 DO j = 1, natorb_b
626 lb = laob(j)
627 nb = naob(j)
628 DO i = 1, natorb_a
629 la = laoa(i)
630 na = naoa(i)
631 hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
632 IF (iatom <= jatom) THEN
633 foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
634 ELSE
635 foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
636 END IF
637 END DO
638 END DO
639 force_ab(ir) = foab
640 END DO
641 DEALLOCATE (dfblock)
642 END IF
643 END IF
644
645 IF (calculate_forces) THEN
646 atom_a = atom_of_kind(iatom)
647 atom_b = atom_of_kind(jatom)
648 IF (irow == iatom) force_ab = -force_ab
649 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
650 force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
651 IF (use_virial) THEN
652 f1 = 1.0_dp
653 IF (iatom == jatom) f1 = 0.5_dp
654 CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
655 END IF
656 END IF
657
658 DEALLOCATE (oint, owork, sint)
659
660 END DO
661 CALL neighbor_list_iterator_release(nl_iterator)
662
663 DO i = 1, SIZE(matrix_h, 1)
664 DO img = 1, nimg
665 CALL dbcsr_finalize(matrix_h(i, img)%matrix)
666 CALL dbcsr_finalize(matrix_s(i, img)%matrix)
667 END DO
668 END DO
669
670 ! EEQ forces (response and direct)
671 IF (calculate_forces) THEN
672 CALL para_env%sum(dcharges)
673 ALLOCATE (qlagrange(natom))
674 CALL xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
675 END IF
676
677 kf = xtb_control%kf
678 enscale = xtb_control%enscale
679 erep = 0.0_dp
680 CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
681
682 esrb = 0.0_dp
683 CALL srb_potential(qs_env, esrb, calculate_forces, xtb_control, cnumbers, dcnum)
684
685 enonbonded = 0.0_dp
686 IF (do_nonbonded) THEN
687 ! nonbonded interactions
688 NULLIFY (sab_xtb_nonbond)
689 CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond)
690 CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
691 atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
692 END IF
693
694 ! set repulsive energy
695 erep = erep + esrb + enonbonded
696 IF (do_nonbonded) THEN
697 CALL para_env%sum(enonbonded)
698 energy%xtb_nonbonded = enonbonded
699 END IF
700 CALL para_env%sum(esrb)
701 energy%srb = esrb
702 CALL para_env%sum(erep)
703 energy%repulsive = erep
704
705 ! save EEQ charges
706 NULLIFY (eeq_q)
707 CALL get_qs_env(qs_env, eeq=eeq_q)
708 IF (ASSOCIATED(eeq_q)) THEN
709 cpassert(SIZE(eeq_q) == natom)
710 ELSE
711 ALLOCATE (eeq_q(natom))
712 eeq_q(1:natom) = charges(1:natom)
713 END IF
714 CALL set_qs_env(qs_env, eeq=eeq_q)
715
716 ! deallocate coordination numbers
717 CALL cnumber_release(cnumbers, dcnum, calculate_forces)
718
719 ! deallocate Huckel parameters
720 DEALLOCATE (huckel)
721 IF (calculate_forces) THEN
722 DEALLOCATE (dhuckel, dqhuckel)
723 END IF
724 ! deallocate KAB parameters
725 DEALLOCATE (kijab)
726
727 ! deallocate charges
728 DEALLOCATE (charges)
729 IF (calculate_forces) THEN
730 DEALLOCATE (dcharges, qlagrange)
731 END IF
732
733 ! AO matrix outputs
734 CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
735
736 DEALLOCATE (basis_set_list)
737 IF (calculate_forces) THEN
738 IF (SIZE(matrix_p, 1) == 2) THEN
739 DO img = 1, nimg
740 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
741 beta_scalar=-1.0_dp)
742 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, alpha_scalar=1.0_dp, &
743 beta_scalar=-1.0_dp)
744 END DO
745 END IF
746 END IF
747
748 CALL timestop(handle)
749
750 END SUBROUTINE build_gfn0_xtb_matrices
751
752! **************************************************************************************************
753!> \brief ...
754!> \param qs_env ...
755!> \param calculate_forces ...
756! **************************************************************************************************
757 SUBROUTINE build_gfn1_xtb_matrices(qs_env, calculate_forces)
758
759 TYPE(qs_environment_type), POINTER :: qs_env
760 LOGICAL, INTENT(IN) :: calculate_forces
761
762 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_gfn1_xtb_matrices'
763
764 INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, iset, &
765 j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, lmaxa, lmaxb, maxder, n1, n2, na, &
766 natom, natorb_a, natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, &
767 nsetb, sgfa, sgfb, za, zb
768 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
769 INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
770 INTEGER, DIMENSION(3) :: cell
771 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
772 npgfb, nsgfa, nsgfb
773 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
774 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
775 LOGICAL :: defined, diagblock, do_nonbonded, found, &
776 use_virial, xb_inter
777 REAL(kind=dp) :: dfp, dhij, dr, drk, drx, enonbonded, &
778 enscale, erep, etaa, etab, exb, f0, &
779 f1, fhua, fhub, fhud, foab, hij, kf, &
780 rcova, rcovab, rcovb, rrab
781 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
782 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, huckel, owork
783 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
784 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
785 REAL(kind=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_ab, rij, rik
786 REAL(kind=dp), DIMENSION(5) :: dpia, dpib, kpolya, kpolyb, pia, pib
787 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
788 REAL(kind=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
789 scon_a, scon_b, wblock, zeta, zetb
790 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
791 TYPE(atprop_type), POINTER :: atprop
792 TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
793 TYPE(cp_logger_type), POINTER :: logger
794 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
795 TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
796 TYPE(dft_control_type), POINTER :: dft_control
797 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
798 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
799 TYPE(kpoint_type), POINTER :: kpoints
800 TYPE(mp_para_env_type), POINTER :: para_env
801 TYPE(neighbor_list_iterator_p_type), &
802 DIMENSION(:), POINTER :: nl_iterator
803 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
804 POINTER :: sab_orb, sab_xtb_nonbond
805 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
806 TYPE(qs_dispersion_type), POINTER :: dispersion_env
807 TYPE(qs_energy_type), POINTER :: energy
808 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
809 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
810 TYPE(qs_ks_env_type), POINTER :: ks_env
811 TYPE(qs_rho_type), POINTER :: rho
812 TYPE(virial_type), POINTER :: virial
813 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
814 TYPE(xtb_control_type), POINTER :: xtb_control
815
816 CALL timeset(routinen, handle)
817
818 NULLIFY (logger, virial, atprop)
819 logger => cp_get_default_logger()
820
821 NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
822 qs_kind_set, sab_orb, ks_env)
823
824 CALL get_qs_env(qs_env=qs_env, &
825 ks_env=ks_env, &
826 energy=energy, &
827 atomic_kind_set=atomic_kind_set, &
828 qs_kind_set=qs_kind_set, &
829 matrix_h_kp=matrix_h, &
830 matrix_s_kp=matrix_s, &
831 para_env=para_env, &
832 atprop=atprop, &
833 dft_control=dft_control, &
834 sab_orb=sab_orb)
835
836 nkind = SIZE(atomic_kind_set)
837 xtb_control => dft_control%qs_control%xtb_control
838 xb_inter = xtb_control%xb_interaction
839 do_nonbonded = xtb_control%do_nonbonded
840 nimg = dft_control%nimages
841 nderivatives = 0
842 IF (calculate_forces) nderivatives = 1
843 IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
844 maxder = ncoset(nderivatives)
845
846 NULLIFY (particle_set)
847 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
848 natom = SIZE(particle_set)
849 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
850 atom_of_kind=atom_of_kind, kind_of=kind_of)
851
852 IF (calculate_forces) THEN
853 NULLIFY (rho, force, matrix_w)
854 CALL get_qs_env(qs_env=qs_env, &
855 rho=rho, matrix_w_kp=matrix_w, &
856 virial=virial, force=force)
857 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
858
859 IF (SIZE(matrix_p, 1) == 2) THEN
860 DO img = 1, nimg
861 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
862 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
863 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
864 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
865 END DO
866 END IF
867 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
868 END IF
869 ! atomic energy decomposition
870 IF (atprop%energy) THEN
871 CALL atprop_array_init(atprop%atecc, natom)
872 END IF
873
874 NULLIFY (cell_to_index)
875 IF (nimg > 1) THEN
876 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
877 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
878 END IF
879
880 ! set up basis set lists
881 ALLOCATE (basis_set_list(nkind))
882 CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
883
884 ! allocate overlap matrix
885 CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
886 CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
887 sab_orb, .true.)
888 CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
889
890 ! initialize H matrix
891 CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
892 DO img = 1, nimg
893 ALLOCATE (matrix_h(1, img)%matrix)
894 CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
895 name="HAMILTONIAN MATRIX")
896 CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
897 END DO
898 CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
899
900 ! Calculate coordination numbers
901 ! needed for effective atomic energy levels (Eq. 12)
902 ! code taken from D3 dispersion energy
903 CALL cnumber_init(qs_env, cnumbers, dcnum, 1, calculate_forces)
904
905 ! vdW Potential
906 CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
907 CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
908 energy%dispersion, calculate_forces)
909
910 ! Calculate Huckel parameters
911 CALL gfn1_huckel(qs_env, cnumbers, huckel, dhuckel, calculate_forces)
912
913 ! Calculate KAB parameters and electronegativity correction
914 CALL gfn1_kpair(qs_env, kijab)
915
916 ! loop over all atom pairs with a non-zero overlap (sab_orb)
917 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
918 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
919 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
920 iatom=iatom, jatom=jatom, r=rij, cell=cell)
921 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
922 CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
923 IF (.NOT. defined .OR. natorb_a < 1) cycle
924 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
925 CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
926 IF (.NOT. defined .OR. natorb_b < 1) cycle
927
928 dr = sqrt(sum(rij(:)**2))
929
930 ! atomic parameters
931 CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
932 lmax=lmaxa, nshell=nsa, kpoly=kpolya)
933 CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
934 lmax=lmaxb, nshell=nsb, kpoly=kpolyb)
935
936 IF (nimg == 1) THEN
937 ic = 1
938 ELSE
939 ic = cell_to_index(cell(1), cell(2), cell(3))
940 cpassert(ic > 0)
941 END IF
942
943 icol = max(iatom, jatom)
944 irow = min(iatom, jatom)
945 NULLIFY (sblock, fblock)
946 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
947 row=irow, col=icol, block=sblock, found=found)
948 cpassert(found)
949 CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
950 row=irow, col=icol, block=fblock, found=found)
951 cpassert(found)
952
953 IF (calculate_forces) THEN
954 NULLIFY (pblock)
955 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
956 row=irow, col=icol, block=pblock, found=found)
957 cpassert(found)
958 NULLIFY (wblock)
959 CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
960 row=irow, col=icol, block=wblock, found=found)
961 cpassert(found)
962 DO i = 2, 4
963 NULLIFY (dsblocks(i)%block)
964 CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
965 row=irow, col=icol, block=dsblocks(i)%block, found=found)
966 cpassert(found)
967 END DO
968 END IF
969
970 ! overlap
971 basis_set_a => basis_set_list(ikind)%gto_basis_set
972 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
973 basis_set_b => basis_set_list(jkind)%gto_basis_set
974 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
975 atom_a = atom_of_kind(iatom)
976 atom_b = atom_of_kind(jatom)
977 ! basis ikind
978 first_sgfa => basis_set_a%first_sgf
979 la_max => basis_set_a%lmax
980 la_min => basis_set_a%lmin
981 npgfa => basis_set_a%npgf
982 nseta = basis_set_a%nset
983 nsgfa => basis_set_a%nsgf_set
984 rpgfa => basis_set_a%pgf_radius
985 set_radius_a => basis_set_a%set_radius
986 scon_a => basis_set_a%scon
987 zeta => basis_set_a%zet
988 ! basis jkind
989 first_sgfb => basis_set_b%first_sgf
990 lb_max => basis_set_b%lmax
991 lb_min => basis_set_b%lmin
992 npgfb => basis_set_b%npgf
993 nsetb = basis_set_b%nset
994 nsgfb => basis_set_b%nsgf_set
995 rpgfb => basis_set_b%pgf_radius
996 set_radius_b => basis_set_b%set_radius
997 scon_b => basis_set_b%scon
998 zetb => basis_set_b%zet
999
1000 ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
1001 ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
1002 ALLOCATE (sint(natorb_a, natorb_b, maxder))
1003 sint = 0.0_dp
1004
1005 DO iset = 1, nseta
1006 ncoa = npgfa(iset)*ncoset(la_max(iset))
1007 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1008 sgfa = first_sgfa(1, iset)
1009 DO jset = 1, nsetb
1010 IF (set_radius_a(iset) + set_radius_b(jset) < dr) cycle
1011 ncob = npgfb(jset)*ncoset(lb_max(jset))
1012 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1013 sgfb = first_sgfb(1, jset)
1014 IF (calculate_forces) THEN
1015 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1016 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1017 rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
1018 ELSE
1019 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1020 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1021 rij, sab=oint(:, :, 1))
1022 END IF
1023 ! Contraction
1024 CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1025 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
1026 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.false.)
1027 IF (calculate_forces) THEN
1028 DO i = 2, 4
1029 CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1030 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
1031 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.false.)
1032 END DO
1033 END IF
1034 END DO
1035 END DO
1036 ! forces W matrix
1037 IF (calculate_forces) THEN
1038 DO i = 1, 3
1039 IF (iatom <= jatom) THEN
1040 force_ab(i) = sum(sint(:, :, i + 1)*wblock(:, :))
1041 ELSE
1042 force_ab(i) = sum(sint(:, :, i + 1)*transpose(wblock(:, :)))
1043 END IF
1044 END DO
1045 f1 = 2.0_dp
1046 force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
1047 force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
1048 IF (use_virial .AND. dr > 1.e-3_dp) THEN
1049 IF (iatom == jatom) f1 = 1.0_dp
1050 CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1051 END IF
1052 END IF
1053 ! update S matrix
1054 IF (iatom <= jatom) THEN
1055 sblock(:, :) = sblock(:, :) + sint(:, :, 1)
1056 ELSE
1057 sblock(:, :) = sblock(:, :) + transpose(sint(:, :, 1))
1058 END IF
1059 IF (calculate_forces) THEN
1060 DO i = 2, 4
1061 IF (iatom <= jatom) THEN
1062 dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
1063 ELSE
1064 dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - transpose(sint(:, :, i))
1065 END IF
1066 END DO
1067 END IF
1068
1069 ! Calculate Pi = Pia * Pib (Eq. 11)
1070 rcovab = rcova + rcovb
1071 rrab = sqrt(dr/rcovab)
1072 pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
1073 pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
1074 IF (calculate_forces) THEN
1075 IF (dr > 1.e-6_dp) THEN
1076 drx = 0.5_dp/rrab/rcovab
1077 ELSE
1078 drx = 0.0_dp
1079 END IF
1080 dpia(1:nsa) = drx*kpolya(1:nsa)
1081 dpib(1:nsb) = drx*kpolyb(1:nsb)
1082 END IF
1083
1084 ! diagonal block
1085 diagblock = .false.
1086 IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .true.
1087 !
1088 ! Eq. 10
1089 !
1090 IF (diagblock) THEN
1091 DO i = 1, natorb_a
1092 na = naoa(i)
1093 fblock(i, i) = fblock(i, i) + huckel(na, iatom)
1094 END DO
1095 ELSE
1096 DO j = 1, natorb_b
1097 nb = naob(j)
1098 DO i = 1, natorb_a
1099 na = naoa(i)
1100 hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1101 IF (iatom <= jatom) THEN
1102 fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1103 ELSE
1104 fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1105 END IF
1106 END DO
1107 END DO
1108 END IF
1109 IF (calculate_forces) THEN
1110 f0 = 1.0_dp
1111 IF (irow == iatom) f0 = -1.0_dp
1112 ! Derivative wrt coordination number
1113 fhua = 0.0_dp
1114 fhub = 0.0_dp
1115 fhud = 0.0_dp
1116 IF (diagblock) THEN
1117 DO i = 1, natorb_a
1118 la = laoa(i)
1119 na = naoa(i)
1120 fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
1121 END DO
1122 ELSE
1123 DO j = 1, natorb_b
1124 lb = laob(j)
1125 nb = naob(j)
1126 DO i = 1, natorb_a
1127 la = laoa(i)
1128 na = naoa(i)
1129 hij = 0.5_dp*pia(na)*pib(nb)
1130 IF (iatom <= jatom) THEN
1131 fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(na, iatom)
1132 fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(nb, jatom)
1133 ELSE
1134 fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(na, iatom)
1135 fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(nb, jatom)
1136 END IF
1137 END DO
1138 END DO
1139 IF (iatom /= jatom) THEN
1140 fhua = 2.0_dp*fhua
1141 fhub = 2.0_dp*fhub
1142 END IF
1143 END IF
1144 ! iatom
1145 atom_a = atom_of_kind(iatom)
1146 DO i = 1, dcnum(iatom)%neighbors
1147 katom = dcnum(iatom)%nlist(i)
1148 kkind = kind_of(katom)
1149 atom_c = atom_of_kind(katom)
1150 rik = dcnum(iatom)%rik(:, i)
1151 drk = sqrt(sum(rik(:)**2))
1152 IF (drk > 1.e-3_dp) THEN
1153 fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
1154 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
1155 force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
1156 fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
1157 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
1158 force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
1159 IF (use_virial) THEN
1160 fdik = fdika + fdikb
1161 CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
1162 END IF
1163 END IF
1164 END DO
1165 ! jatom
1166 atom_b = atom_of_kind(jatom)
1167 DO i = 1, dcnum(jatom)%neighbors
1168 katom = dcnum(jatom)%nlist(i)
1169 kkind = kind_of(katom)
1170 atom_c = atom_of_kind(katom)
1171 rik = dcnum(jatom)%rik(:, i)
1172 drk = sqrt(sum(rik(:)**2))
1173 IF (drk > 1.e-3_dp) THEN
1174 fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
1175 force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
1176 force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
1177 IF (use_virial) THEN
1178 CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
1179 END IF
1180 END IF
1181 END DO
1182 ! force from R dendent Huckel element: Pia*Pib
1183 IF (diagblock) THEN
1184 force_ab = 0._dp
1185 ELSE
1186 n1 = SIZE(fblock, 1)
1187 n2 = SIZE(fblock, 2)
1188 ALLOCATE (dfblock(n1, n2))
1189 dfblock = 0.0_dp
1190 DO j = 1, natorb_b
1191 lb = laob(j)
1192 nb = naob(j)
1193 DO i = 1, natorb_a
1194 la = laoa(i)
1195 na = naoa(i)
1196 dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
1197 IF (iatom <= jatom) THEN
1198 dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1199 ELSE
1200 dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1201 END IF
1202 END DO
1203 END DO
1204 dfp = f0*sum(dfblock(:, :)*pblock(:, :))
1205 DO ir = 1, 3
1206 foab = 2.0_dp*dfp*rij(ir)/dr
1207 ! force from overlap matrix contribution to H
1208 DO j = 1, natorb_b
1209 lb = laob(j)
1210 nb = naob(j)
1211 DO i = 1, natorb_a
1212 la = laoa(i)
1213 na = naoa(i)
1214 hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1215 IF (iatom <= jatom) THEN
1216 foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
1217 ELSE
1218 foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
1219 END IF
1220 END DO
1221 END DO
1222 force_ab(ir) = foab
1223 END DO
1224 DEALLOCATE (dfblock)
1225 END IF
1226 END IF
1227
1228 IF (calculate_forces) THEN
1229 atom_a = atom_of_kind(iatom)
1230 atom_b = atom_of_kind(jatom)
1231 IF (irow == iatom) force_ab = -force_ab
1232 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
1233 force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
1234 IF (use_virial) THEN
1235 f1 = 1.0_dp
1236 IF (iatom == jatom) f1 = 0.5_dp
1237 CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1238 END IF
1239 END IF
1240
1241 DEALLOCATE (oint, owork, sint)
1242
1243 END DO
1244 CALL neighbor_list_iterator_release(nl_iterator)
1245
1246 DO i = 1, SIZE(matrix_h, 1)
1247 DO img = 1, nimg
1248 CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1249 CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1250 END DO
1251 END DO
1252
1253 kf = xtb_control%kf
1254 enscale = xtb_control%enscale
1255 erep = 0.0_dp
1256 CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
1257
1258 exb = 0.0_dp
1259 IF (xb_inter) THEN
1260 CALL xb_interaction(qs_env, exb, calculate_forces)
1261 END IF
1262
1263 enonbonded = 0.0_dp
1264 IF (do_nonbonded) THEN
1265 ! nonbonded interactions
1266 NULLIFY (sab_xtb_nonbond)
1267 CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond)
1268 CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
1269 atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
1270 END IF
1271
1272 ! set repulsive energy
1273 erep = erep + exb + enonbonded
1274 IF (xb_inter) THEN
1275 CALL para_env%sum(exb)
1276 energy%xtb_xb_inter = exb
1277 END IF
1278 IF (do_nonbonded) THEN
1279 CALL para_env%sum(enonbonded)
1280 energy%xtb_nonbonded = enonbonded
1281 END IF
1282 CALL para_env%sum(erep)
1283 energy%repulsive = erep
1284
1285 ! deallocate coordination numbers
1286 CALL cnumber_release(cnumbers, dcnum, calculate_forces)
1287
1288 ! deallocate Huckel parameters
1289 DEALLOCATE (huckel)
1290 IF (calculate_forces) THEN
1291 DEALLOCATE (dhuckel)
1292 END IF
1293 ! deallocate KAB parameters
1294 DEALLOCATE (kijab)
1295
1296 ! AO matrix outputs
1297 CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
1298
1299 DEALLOCATE (basis_set_list)
1300 IF (calculate_forces) THEN
1301 IF (SIZE(matrix_p, 1) == 2) THEN
1302 DO img = 1, nimg
1303 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
1304 beta_scalar=-1.0_dp)
1305 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, alpha_scalar=1.0_dp, &
1306 beta_scalar=-1.0_dp)
1307 END DO
1308 END IF
1309 END IF
1310
1311 CALL timestop(handle)
1312
1313 END SUBROUTINE build_gfn1_xtb_matrices
1314
1315! **************************************************************************************************
1316!> \brief ...
1317!> \param qs_env ...
1318!> \param matrix_h ...
1319!> \param matrix_s ...
1320!> \param calculate_forces ...
1321! **************************************************************************************************
1322 SUBROUTINE ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
1323 TYPE(qs_environment_type), POINTER :: qs_env
1324 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s
1325 LOGICAL, INTENT(IN) :: calculate_forces
1326
1327 INTEGER :: after, i, img, iw, nimg
1328 LOGICAL :: norml1, norml2, omit_headers, use_arnoldi
1329 REAL(kind=dp), DIMENSION(2) :: condnum
1330 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1331 TYPE(cp_logger_type), POINTER :: logger
1332 TYPE(mp_para_env_type), POINTER :: para_env
1333
1334 logger => cp_get_default_logger()
1335
1336 CALL get_qs_env(qs_env, para_env=para_env)
1337 nimg = SIZE(matrix_h, 2)
1338 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
1339 IF (btest(cp_print_key_should_output(logger%iter_info, &
1340 qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
1341 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
1342 extension=".Log")
1343 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
1344 after = min(max(after, 1), 16)
1345 DO img = 1, nimg
1346 CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, &
1347 output_unit=iw, omit_headers=omit_headers)
1348 END DO
1349 CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
1350 END IF
1351
1352 IF (btest(cp_print_key_should_output(logger%iter_info, &
1353 qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
1354 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
1355 extension=".Log")
1356 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
1357 after = min(max(after, 1), 16)
1358 DO img = 1, nimg
1359 CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, &
1360 output_unit=iw, omit_headers=omit_headers)
1361 IF (btest(cp_print_key_should_output(logger%iter_info, &
1362 qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
1363 DO i = 2, SIZE(matrix_s, 1)
1364 CALL cp_dbcsr_write_sparse_matrix(matrix_s(i, img)%matrix, 4, after, qs_env, para_env, &
1365 output_unit=iw, omit_headers=omit_headers)
1366 END DO
1367 END IF
1368 END DO
1369 CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP")
1370 END IF
1371
1372 ! *** Overlap condition number
1373 IF (.NOT. calculate_forces) THEN
1374 IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
1375 "DFT%PRINT%OVERLAP_CONDITION") /= 0) THEN
1376 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
1377 extension=".Log")
1378 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
1379 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
1380 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
1381 CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
1382 CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
1383 END IF
1384 END IF
1385
1386 END SUBROUTINE ao_matrix_output
1387
1388END MODULE xtb_matrices
ai_contraction::block_add
Definition ai_contraction.F:43
ai_contraction::contraction
Definition ai_contraction.F:31
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:95
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, rr_work)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:273
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
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
atprop_types::atprop_array_init
subroutine, public atprop_array_init(atarray, natom)
...
Definition atprop_types.F:87
basis_set_types
Definition basis_set_types.F:15
block_p_types
collect pointers to a block of reals
Definition block_p_types.F:14
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_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_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_output
DBCSR output in CP2K.
Definition cp_dbcsr_output.F:17
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers, cartesian_basis)
...
Definition cp_dbcsr_output.F:165
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_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:1073
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
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
eeq_input
Input definition and setup for EEQ model.
Definition eeq_input.F:12
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::vdw_pairpot_dftd4
integer, parameter, public vdw_pairpot_dftd4
Definition input_constants.F:640
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_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:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
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, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:568
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
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_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_dispersion_cnum
Coordination number routines for dispersion pairpotentials.
Definition qs_dispersion_cnum.F:12
qs_dispersion_cnum::cnumber_release
subroutine, public cnumber_release(cnumbers, dcnum, derivatives)
...
Definition qs_dispersion_cnum.F:1341
qs_dispersion_cnum::cnumber_init
subroutine, public cnumber_init(qs_env, cnumbers, dcnum, ftype, derivatives, disp_env)
...
Definition qs_dispersion_cnum.F:1279
qs_dispersion_pairpot
Calculation of dispersion using pair potentials.
Definition qs_dispersion_pairpot.F:12
qs_dispersion_pairpot::calculate_dispersion_pairpot
subroutine, public calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
...
Definition qs_dispersion_pairpot.F:393
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition qs_dispersion_types.F:12
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, mimic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1134
qs_force_types
Definition qs_force_types.F:13
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:591
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, 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, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition qs_ks_types.F:343
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:550
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
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F:142
virial_types
Definition virial_types.F:13
xtb_eeq
Calculation of charge equilibration in xTB.
Definition xtb_eeq.F:12
xtb_eeq::xtb_eeq_calculation
subroutine, public xtb_eeq_calculation(qs_env, charges, cnumbers, eeq_sparam, eeq_energy, ef_energy, lambda)
...
Definition xtb_eeq.F:79
xtb_eeq::xtb_eeq_forces
subroutine, public xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
...
Definition xtb_eeq.F:229
xtb_hcore
Calculation of EHT matrix elements in xTB Reference: Stefan Grimme, Christoph Bannwarth,...
Definition xtb_hcore.F:15
xtb_hcore::gfn0_kpair
subroutine, public gfn0_kpair(qs_env, kijab)
...
Definition xtb_hcore.F:203
xtb_hcore::gfn0_huckel
subroutine, public gfn0_huckel(qs_env, cnumbers, charges, huckel, dhuckel, dqhuckel, calculate_forces)
...
Definition xtb_hcore.F:57
xtb_hcore::gfn1_kpair
subroutine, public gfn1_kpair(qs_env, kijab)
...
Definition xtb_hcore.F:322
xtb_hcore::gfn1_huckel
subroutine, public gfn1_huckel(qs_env, cnumbers, huckel, dhuckel, calculate_forces)
...
Definition xtb_hcore.F:122
xtb_matrices
Calculation of Overlap and Hamiltonian matrices in xTB Reference: Stefan Grimme, Christoph Bannwarth,...
Definition xtb_matrices.F:15
xtb_matrices::build_xtb_matrices
subroutine, public build_xtb_matrices(qs_env, calculate_forces)
...
Definition xtb_matrices.F:111
xtb_potentials
xTB (repulsive) pair potentials Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov JCTC 1...
Definition xtb_potentials.F:15
xtb_potentials::nonbonded_correction
subroutine, public nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
Computes a correction for nonbonded interactions based on a generic potential.
Definition xtb_potentials.F:623
xtb_potentials::srb_potential
subroutine, public srb_potential(qs_env, esrb, calculate_forces, xtb_control, cnumbers, dcnum)
...
Definition xtb_potentials.F:202
xtb_potentials::xb_interaction
subroutine, public xb_interaction(qs_env, exb, calculate_forces)
...
Definition xtb_potentials.F:399
xtb_potentials::repulsive_potential
subroutine, public repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
...
Definition xtb_potentials.F:94
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, wall, hen, zeta, xi, kappa0, alpg, occupation, ngauss, electronegativity, chmax, en, kqat2, kcn, kq)
...
Definition xtb_types.F:206
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
block_p_types::block_p_type
Definition block_p_types.F:23
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:761
cp_control_types::xtb_control_type
Definition cp_control_types.F:256
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
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
eeq_input::eeq_solver_type
Definition eeq_input.F:27
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:188
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
particle_types::particle_type
Definition particle_types.F:35
qs_dispersion_cnum::dcnum_type
Definition qs_dispersion_cnum.F:48
qs_dispersion_types::qs_dispersion_type
Definition qs_dispersion_types.F:34
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:224
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:190
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
virial_types::virial_type
Definition virial_types.F:26
xtb_types::xtb_atom_type
Definition xtb_types.F:42