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