(git:ed6f26b)
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-2025 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 END DO
742 END IF
743 END IF
744
745 CALL timestop(handle)
746
747 END SUBROUTINE build_gfn0_xtb_matrices
748
749! **************************************************************************************************
750!> \brief ...
751!> \param qs_env ...
752!> \param calculate_forces ...
753! **************************************************************************************************
754 SUBROUTINE build_gfn1_xtb_matrices(qs_env, calculate_forces)
755
756 TYPE(qs_environment_type), POINTER :: qs_env
757 LOGICAL, INTENT(IN) :: calculate_forces
758
759 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_gfn1_xtb_matrices'
760
761 INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, iset, &
762 j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, lmaxa, lmaxb, maxder, n1, n2, na, &
763 natom, natorb_a, natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, &
764 nsetb, sgfa, sgfb, za, zb
765 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
766 INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
767 INTEGER, DIMENSION(3) :: cell
768 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
769 npgfb, nsgfa, nsgfb
770 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
771 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
772 LOGICAL :: defined, diagblock, do_nonbonded, found, &
773 use_virial, xb_inter
774 REAL(kind=dp) :: dfp, dhij, dr, drk, drx, enonbonded, &
775 enscale, erep, etaa, etab, exb, f0, &
776 f1, fhua, fhub, fhud, foab, hij, kf, &
777 rcova, rcovab, rcovb, rrab
778 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
779 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, huckel, owork
780 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
781 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
782 REAL(kind=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_ab, rij, rik
783 REAL(kind=dp), DIMENSION(5) :: dpia, dpib, kpolya, kpolyb, pia, pib
784 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
785 REAL(kind=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
786 scon_a, scon_b, wblock, zeta, zetb
787 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
788 TYPE(atprop_type), POINTER :: atprop
789 TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
790 TYPE(cp_logger_type), POINTER :: logger
791 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
792 TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
793 TYPE(dft_control_type), POINTER :: dft_control
794 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
795 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
796 TYPE(kpoint_type), POINTER :: kpoints
797 TYPE(mp_para_env_type), POINTER :: para_env
798 TYPE(neighbor_list_iterator_p_type), &
799 DIMENSION(:), POINTER :: nl_iterator
800 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
801 POINTER :: sab_orb, sab_xtb_nonbond
802 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
803 TYPE(qs_dispersion_type), POINTER :: dispersion_env
804 TYPE(qs_energy_type), POINTER :: energy
805 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
806 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
807 TYPE(qs_ks_env_type), POINTER :: ks_env
808 TYPE(qs_rho_type), POINTER :: rho
809 TYPE(virial_type), POINTER :: virial
810 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
811 TYPE(xtb_control_type), POINTER :: xtb_control
812
813 CALL timeset(routinen, handle)
814
815 NULLIFY (logger, virial, atprop)
816 logger => cp_get_default_logger()
817
818 NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
819 qs_kind_set, sab_orb, ks_env)
820
821 CALL get_qs_env(qs_env=qs_env, &
822 ks_env=ks_env, &
823 energy=energy, &
824 atomic_kind_set=atomic_kind_set, &
825 qs_kind_set=qs_kind_set, &
826 matrix_h_kp=matrix_h, &
827 matrix_s_kp=matrix_s, &
828 para_env=para_env, &
829 atprop=atprop, &
830 dft_control=dft_control, &
831 sab_orb=sab_orb)
832
833 nkind = SIZE(atomic_kind_set)
834 xtb_control => dft_control%qs_control%xtb_control
835 xb_inter = xtb_control%xb_interaction
836 do_nonbonded = xtb_control%do_nonbonded
837 nimg = dft_control%nimages
838 nderivatives = 0
839 IF (calculate_forces) nderivatives = 1
840 IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
841 maxder = ncoset(nderivatives)
842
843 NULLIFY (particle_set)
844 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
845 natom = SIZE(particle_set)
846 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
847 atom_of_kind=atom_of_kind, kind_of=kind_of)
848
849 IF (calculate_forces) THEN
850 NULLIFY (rho, force, matrix_w)
851 CALL get_qs_env(qs_env=qs_env, &
852 rho=rho, matrix_w_kp=matrix_w, &
853 virial=virial, force=force)
854 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
855
856 IF (SIZE(matrix_p, 1) == 2) THEN
857 DO img = 1, nimg
858 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
859 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
860 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
861 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
862 END DO
863 END IF
864 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
865 END IF
866 ! atomic energy decomposition
867 IF (atprop%energy) THEN
868 CALL atprop_array_init(atprop%atecc, natom)
869 END IF
870
871 NULLIFY (cell_to_index)
872 IF (nimg > 1) THEN
873 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
874 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
875 END IF
876
877 ! set up basis set lists
878 ALLOCATE (basis_set_list(nkind))
879 CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
880
881 ! allocate overlap matrix
882 CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
883 CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
884 sab_orb, .true.)
885 CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
886
887 ! initialize H matrix
888 CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
889 DO img = 1, nimg
890 ALLOCATE (matrix_h(1, img)%matrix)
891 CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
892 name="HAMILTONIAN MATRIX")
893 CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
894 END DO
895 CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
896
897 ! Calculate coordination numbers
898 ! needed for effective atomic energy levels (Eq. 12)
899 ! code taken from D3 dispersion energy
900 CALL cnumber_init(qs_env, cnumbers, dcnum, 1, calculate_forces)
901
902 ! vdW Potential
903 CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
904 CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
905 energy%dispersion, calculate_forces)
906
907 ! Calculate Huckel parameters
908 CALL gfn1_huckel(qs_env, cnumbers, huckel, dhuckel, calculate_forces)
909
910 ! Calculate KAB parameters and electronegativity correction
911 CALL gfn1_kpair(qs_env, kijab)
912
913 ! loop over all atom pairs with a non-zero overlap (sab_orb)
914 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
915 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
916 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
917 iatom=iatom, jatom=jatom, r=rij, cell=cell)
918 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
919 CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
920 IF (.NOT. defined .OR. natorb_a < 1) cycle
921 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
922 CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
923 IF (.NOT. defined .OR. natorb_b < 1) cycle
924
925 dr = sqrt(sum(rij(:)**2))
926
927 ! atomic parameters
928 CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
929 lmax=lmaxa, nshell=nsa, kpoly=kpolya)
930 CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
931 lmax=lmaxb, nshell=nsb, kpoly=kpolyb)
932
933 IF (nimg == 1) THEN
934 ic = 1
935 ELSE
936 ic = cell_to_index(cell(1), cell(2), cell(3))
937 cpassert(ic > 0)
938 END IF
939
940 icol = max(iatom, jatom)
941 irow = min(iatom, jatom)
942 NULLIFY (sblock, fblock)
943 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
944 row=irow, col=icol, block=sblock, found=found)
945 cpassert(found)
946 CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
947 row=irow, col=icol, block=fblock, found=found)
948 cpassert(found)
949
950 IF (calculate_forces) THEN
951 NULLIFY (pblock)
952 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
953 row=irow, col=icol, block=pblock, found=found)
954 cpassert(found)
955 NULLIFY (wblock)
956 CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
957 row=irow, col=icol, block=wblock, found=found)
958 cpassert(found)
959 DO i = 2, 4
960 NULLIFY (dsblocks(i)%block)
961 CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
962 row=irow, col=icol, block=dsblocks(i)%block, found=found)
963 cpassert(found)
964 END DO
965 END IF
966
967 ! overlap
968 basis_set_a => basis_set_list(ikind)%gto_basis_set
969 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
970 basis_set_b => basis_set_list(jkind)%gto_basis_set
971 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
972 atom_a = atom_of_kind(iatom)
973 atom_b = atom_of_kind(jatom)
974 ! basis ikind
975 first_sgfa => basis_set_a%first_sgf
976 la_max => basis_set_a%lmax
977 la_min => basis_set_a%lmin
978 npgfa => basis_set_a%npgf
979 nseta = basis_set_a%nset
980 nsgfa => basis_set_a%nsgf_set
981 rpgfa => basis_set_a%pgf_radius
982 set_radius_a => basis_set_a%set_radius
983 scon_a => basis_set_a%scon
984 zeta => basis_set_a%zet
985 ! basis jkind
986 first_sgfb => basis_set_b%first_sgf
987 lb_max => basis_set_b%lmax
988 lb_min => basis_set_b%lmin
989 npgfb => basis_set_b%npgf
990 nsetb = basis_set_b%nset
991 nsgfb => basis_set_b%nsgf_set
992 rpgfb => basis_set_b%pgf_radius
993 set_radius_b => basis_set_b%set_radius
994 scon_b => basis_set_b%scon
995 zetb => basis_set_b%zet
996
997 ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
998 ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
999 ALLOCATE (sint(natorb_a, natorb_b, maxder))
1000 sint = 0.0_dp
1001
1002 DO iset = 1, nseta
1003 ncoa = npgfa(iset)*ncoset(la_max(iset))
1004 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1005 sgfa = first_sgfa(1, iset)
1006 DO jset = 1, nsetb
1007 IF (set_radius_a(iset) + set_radius_b(jset) < dr) cycle
1008 ncob = npgfb(jset)*ncoset(lb_max(jset))
1009 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1010 sgfb = first_sgfb(1, jset)
1011 IF (calculate_forces) THEN
1012 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1013 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1014 rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
1015 ELSE
1016 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1017 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1018 rij, sab=oint(:, :, 1))
1019 END IF
1020 ! Contraction
1021 CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1022 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
1023 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.false.)
1024 IF (calculate_forces) THEN
1025 DO i = 2, 4
1026 CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1027 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
1028 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.false.)
1029 END DO
1030 END IF
1031 END DO
1032 END DO
1033 ! forces W matrix
1034 IF (calculate_forces) THEN
1035 DO i = 1, 3
1036 IF (iatom <= jatom) THEN
1037 force_ab(i) = sum(sint(:, :, i + 1)*wblock(:, :))
1038 ELSE
1039 force_ab(i) = sum(sint(:, :, i + 1)*transpose(wblock(:, :)))
1040 END IF
1041 END DO
1042 f1 = 2.0_dp
1043 force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
1044 force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
1045 IF (use_virial .AND. dr > 1.e-3_dp) THEN
1046 IF (iatom == jatom) f1 = 1.0_dp
1047 CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1048 END IF
1049 END IF
1050 ! update S matrix
1051 IF (iatom <= jatom) THEN
1052 sblock(:, :) = sblock(:, :) + sint(:, :, 1)
1053 ELSE
1054 sblock(:, :) = sblock(:, :) + transpose(sint(:, :, 1))
1055 END IF
1056 IF (calculate_forces) THEN
1057 DO i = 2, 4
1058 IF (iatom <= jatom) THEN
1059 dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
1060 ELSE
1061 dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - transpose(sint(:, :, i))
1062 END IF
1063 END DO
1064 END IF
1065
1066 ! Calculate Pi = Pia * Pib (Eq. 11)
1067 rcovab = rcova + rcovb
1068 rrab = sqrt(dr/rcovab)
1069 pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
1070 pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
1071 IF (calculate_forces) THEN
1072 IF (dr > 1.e-6_dp) THEN
1073 drx = 0.5_dp/rrab/rcovab
1074 ELSE
1075 drx = 0.0_dp
1076 END IF
1077 dpia(1:nsa) = drx*kpolya(1:nsa)
1078 dpib(1:nsb) = drx*kpolyb(1:nsb)
1079 END IF
1080
1081 ! diagonal block
1082 diagblock = .false.
1083 IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .true.
1084 !
1085 ! Eq. 10
1086 !
1087 IF (diagblock) THEN
1088 DO i = 1, natorb_a
1089 na = naoa(i)
1090 fblock(i, i) = fblock(i, i) + huckel(na, iatom)
1091 END DO
1092 ELSE
1093 DO j = 1, natorb_b
1094 nb = naob(j)
1095 DO i = 1, natorb_a
1096 na = naoa(i)
1097 hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1098 IF (iatom <= jatom) THEN
1099 fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1100 ELSE
1101 fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1102 END IF
1103 END DO
1104 END DO
1105 END IF
1106 IF (calculate_forces) THEN
1107 f0 = 1.0_dp
1108 IF (irow == iatom) f0 = -1.0_dp
1109 ! Derivative wrt coordination number
1110 fhua = 0.0_dp
1111 fhub = 0.0_dp
1112 fhud = 0.0_dp
1113 IF (diagblock) THEN
1114 DO i = 1, natorb_a
1115 la = laoa(i)
1116 na = naoa(i)
1117 fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
1118 END DO
1119 ELSE
1120 DO j = 1, natorb_b
1121 lb = laob(j)
1122 nb = naob(j)
1123 DO i = 1, natorb_a
1124 la = laoa(i)
1125 na = naoa(i)
1126 hij = 0.5_dp*pia(na)*pib(nb)
1127 IF (iatom <= jatom) THEN
1128 fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(na, iatom)
1129 fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(nb, jatom)
1130 ELSE
1131 fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(na, iatom)
1132 fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(nb, jatom)
1133 END IF
1134 END DO
1135 END DO
1136 IF (iatom /= jatom) THEN
1137 fhua = 2.0_dp*fhua
1138 fhub = 2.0_dp*fhub
1139 END IF
1140 END IF
1141 ! iatom
1142 atom_a = atom_of_kind(iatom)
1143 DO i = 1, dcnum(iatom)%neighbors
1144 katom = dcnum(iatom)%nlist(i)
1145 kkind = kind_of(katom)
1146 atom_c = atom_of_kind(katom)
1147 rik = dcnum(iatom)%rik(:, i)
1148 drk = sqrt(sum(rik(:)**2))
1149 IF (drk > 1.e-3_dp) THEN
1150 fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
1151 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
1152 force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
1153 fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
1154 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
1155 force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
1156 IF (use_virial) THEN
1157 fdik = fdika + fdikb
1158 CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
1159 END IF
1160 END IF
1161 END DO
1162 ! jatom
1163 atom_b = atom_of_kind(jatom)
1164 DO i = 1, dcnum(jatom)%neighbors
1165 katom = dcnum(jatom)%nlist(i)
1166 kkind = kind_of(katom)
1167 atom_c = atom_of_kind(katom)
1168 rik = dcnum(jatom)%rik(:, i)
1169 drk = sqrt(sum(rik(:)**2))
1170 IF (drk > 1.e-3_dp) THEN
1171 fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
1172 force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
1173 force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
1174 IF (use_virial) THEN
1175 CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
1176 END IF
1177 END IF
1178 END DO
1179 ! force from R dendent Huckel element: Pia*Pib
1180 IF (diagblock) THEN
1181 force_ab = 0._dp
1182 ELSE
1183 n1 = SIZE(fblock, 1)
1184 n2 = SIZE(fblock, 2)
1185 ALLOCATE (dfblock(n1, n2))
1186 dfblock = 0.0_dp
1187 DO j = 1, natorb_b
1188 lb = laob(j)
1189 nb = naob(j)
1190 DO i = 1, natorb_a
1191 la = laoa(i)
1192 na = naoa(i)
1193 dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
1194 IF (iatom <= jatom) THEN
1195 dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1196 ELSE
1197 dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
1198 END IF
1199 END DO
1200 END DO
1201 dfp = f0*sum(dfblock(:, :)*pblock(:, :))
1202 DO ir = 1, 3
1203 foab = 2.0_dp*dfp*rij(ir)/dr
1204 ! force from overlap matrix contribution to H
1205 DO j = 1, natorb_b
1206 lb = laob(j)
1207 nb = naob(j)
1208 DO i = 1, natorb_a
1209 la = laoa(i)
1210 na = naoa(i)
1211 hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
1212 IF (iatom <= jatom) THEN
1213 foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
1214 ELSE
1215 foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
1216 END IF
1217 END DO
1218 END DO
1219 force_ab(ir) = foab
1220 END DO
1221 DEALLOCATE (dfblock)
1222 END IF
1223 END IF
1224
1225 IF (calculate_forces) THEN
1226 atom_a = atom_of_kind(iatom)
1227 atom_b = atom_of_kind(jatom)
1228 IF (irow == iatom) force_ab = -force_ab
1229 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
1230 force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
1231 IF (use_virial) THEN
1232 f1 = 1.0_dp
1233 IF (iatom == jatom) f1 = 0.5_dp
1234 CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
1235 END IF
1236 END IF
1237
1238 DEALLOCATE (oint, owork, sint)
1239
1240 END DO
1241 CALL neighbor_list_iterator_release(nl_iterator)
1242
1243 DO i = 1, SIZE(matrix_h, 1)
1244 DO img = 1, nimg
1245 CALL dbcsr_finalize(matrix_h(i, img)%matrix)
1246 CALL dbcsr_finalize(matrix_s(i, img)%matrix)
1247 END DO
1248 END DO
1249
1250 kf = xtb_control%kf
1251 enscale = xtb_control%enscale
1252 erep = 0.0_dp
1253 CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
1254
1255 exb = 0.0_dp
1256 IF (xb_inter) THEN
1257 CALL xb_interaction(qs_env, exb, calculate_forces)
1258 END IF
1259
1260 enonbonded = 0.0_dp
1261 IF (do_nonbonded) THEN
1262 ! nonbonded interactions
1263 NULLIFY (sab_xtb_nonbond)
1264 CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond)
1265 CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
1266 atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
1267 END IF
1268
1269 ! set repulsive energy
1270 erep = erep + exb + enonbonded
1271 IF (xb_inter) THEN
1272 CALL para_env%sum(exb)
1273 energy%xtb_xb_inter = exb
1274 END IF
1275 IF (do_nonbonded) THEN
1276 CALL para_env%sum(enonbonded)
1277 energy%xtb_nonbonded = enonbonded
1278 END IF
1279 CALL para_env%sum(erep)
1280 energy%repulsive = erep
1281
1282 ! deallocate coordination numbers
1283 CALL cnumber_release(cnumbers, dcnum, calculate_forces)
1284
1285 ! deallocate Huckel parameters
1286 DEALLOCATE (huckel)
1287 IF (calculate_forces) THEN
1288 DEALLOCATE (dhuckel)
1289 END IF
1290 ! deallocate KAB parameters
1291 DEALLOCATE (kijab)
1292
1293 ! AO matrix outputs
1294 CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
1295
1296 DEALLOCATE (basis_set_list)
1297 IF (calculate_forces) THEN
1298 IF (SIZE(matrix_p, 1) == 2) THEN
1299 DO img = 1, nimg
1300 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
1301 beta_scalar=-1.0_dp)
1302 END DO
1303 END IF
1304 END IF
1305
1306 CALL timestop(handle)
1307
1308 END SUBROUTINE build_gfn1_xtb_matrices
1309
1310! **************************************************************************************************
1311!> \brief ...
1312!> \param qs_env ...
1313!> \param matrix_h ...
1314!> \param matrix_s ...
1315!> \param calculate_forces ...
1316! **************************************************************************************************
1317 SUBROUTINE ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
1318 TYPE(qs_environment_type), POINTER :: qs_env
1319 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s
1320 LOGICAL, INTENT(IN) :: calculate_forces
1321
1322 INTEGER :: after, i, img, iw, nimg
1323 LOGICAL :: norml1, norml2, omit_headers, use_arnoldi
1324 REAL(kind=dp), DIMENSION(2) :: condnum
1325 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1326 TYPE(cp_logger_type), POINTER :: logger
1327 TYPE(mp_para_env_type), POINTER :: para_env
1328
1329 logger => cp_get_default_logger()
1330
1331 CALL get_qs_env(qs_env, para_env=para_env)
1332 nimg = SIZE(matrix_h, 2)
1333 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
1334 IF (btest(cp_print_key_should_output(logger%iter_info, &
1335 qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
1336 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
1337 extension=".Log")
1338 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
1339 after = min(max(after, 1), 16)
1340 DO img = 1, nimg
1341 CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, &
1342 output_unit=iw, omit_headers=omit_headers)
1343 END DO
1344 CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
1345 END IF
1346
1347 IF (btest(cp_print_key_should_output(logger%iter_info, &
1348 qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
1349 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
1350 extension=".Log")
1351 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
1352 after = min(max(after, 1), 16)
1353 DO img = 1, nimg
1354 CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, &
1355 output_unit=iw, omit_headers=omit_headers)
1356 IF (btest(cp_print_key_should_output(logger%iter_info, &
1357 qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
1358 DO i = 2, SIZE(matrix_s, 1)
1359 CALL cp_dbcsr_write_sparse_matrix(matrix_s(i, img)%matrix, 4, after, qs_env, para_env, &
1360 output_unit=iw, omit_headers=omit_headers)
1361 END DO
1362 END IF
1363 END DO
1364 CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP")
1365 END IF
1366
1367 ! *** Overlap condition number
1368 IF (.NOT. calculate_forces) THEN
1369 IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
1370 "DFT%PRINT%OVERLAP_CONDITION") .NE. 0) THEN
1371 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
1372 extension=".Log")
1373 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
1374 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
1375 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
1376 CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
1377 CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
1378 END IF
1379 END IF
1380
1381 END SUBROUTINE ao_matrix_output
1382
1383END MODULE xtb_matrices
1384
ai_contraction::block_add
Definition ai_contraction.F:43
ai_contraction::contraction
Definition ai_contraction.F:31
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
qs_integral_utils::get_memory_usage
Definition qs_integral_utils.F:34
qs_overlap::create_sab_matrix
Definition qs_overlap.F:81
ai_contraction
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Definition ai_contraction.F:18
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
ai_overlap::overlap_ab
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:680
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
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:679
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:659
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_output
DBCSR output in CP2K.
Definition cp_dbcsr_output.F:17
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers)
...
Definition cp_dbcsr_output.F:163
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
cp_output_handling::cp_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:596
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)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:365
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:1426
qs_dispersion_cnum::cnumber_init
subroutine, public cnumber_init(qs_env, cnumbers, dcnum, ftype, derivatives, disp_env)
...
Definition qs_dispersion_cnum.F:1364
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, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, 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, 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)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1055
qs_force_types
Definition qs_force_types.F:13
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:572
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, 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, 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:333
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:570
cp_control_types::xtb_control_type
Definition cp_control_types.F:138
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
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:150
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_dispersion_cnum::dcnum_type
Definition qs_dispersion_cnum.F:48
qs_dispersion_types::qs_dispersion_type
Definition qs_dispersion_types.F:33
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:136
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
virial_types::virial_type
Definition virial_types.F:26
xtb_types::xtb_atom_type
Definition xtb_types.F:42