(git:d43dca4)
Loading...
Searching...
No Matches
qs_dftb_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 DFTB
10!> \author JGH
11! **************************************************************************************************
21 USE cp_dbcsr_api, ONLY: &
22 dbcsr_add, dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, &
24 dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_symmetric
31 USE cp_output_handling, ONLY: cp_p_file,&
40 USE kinds, ONLY: default_string_length,&
41 dp
42 USE kpoint_types, ONLY: get_kpoint_info,&
45 USE mulliken, ONLY: mulliken_charges
54 iptr,&
60 USE qs_kind_types, ONLY: get_qs_kind,&
63 USE qs_ks_types, ONLY: get_ks_env,&
66 USE qs_mo_types, ONLY: get_mo_set,&
74 USE qs_rho_types, ONLY: qs_rho_get,&
78 USE virial_types, ONLY: virial_type
79#include "./base/base_uses.f90"
80
81 IMPLICIT NONE
82
83 INTEGER, DIMENSION(16), PARAMETER :: orbptr = [0, 1, 1, 1, &
84 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3]
85
86 PRIVATE
87
88 CHARACTER(len=*), PARAMETER, PRIVATE :: modulen = 'qs_dftb_matrices'
89 REAL(kind=dp), PARAMETER, PRIVATE :: dftb_fd_deriv_step = 1.0e-3_dp
90
92
93CONTAINS
94
95! **************************************************************************************************
96!> \brief ...
97!> \param qs_env ...
98!> \param para_env ...
99!> \param calculate_forces ...
100! **************************************************************************************************
101 SUBROUTINE build_dftb_matrices(qs_env, para_env, calculate_forces)
102
103 TYPE(qs_environment_type), POINTER :: qs_env
104 TYPE(mp_para_env_type), POINTER :: para_env
105 LOGICAL, INTENT(IN) :: calculate_forces
106
107 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_dftb_matrices'
108
109 INTEGER :: after, atom_a, atom_b, handle, i, iatom, ic, icol, ikind, img, irow, iw, jatom, &
110 jkind, l1, l2, la, lb, llm, lmaxi, lmaxj, m, n1, n2, n_urpoly, natorb_a, natorb_b, &
111 nderivatives, ngrd, ngrdcut, nimg, nkind, spdim
112 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
113 INTEGER, DIMENSION(3) :: cell
114 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
115 LOGICAL :: defined, found, omit_headers, use_virial
116 REAL(kind=dp) :: ddr, dgrd, dr, erep, erepij, f0, foab, &
117 fow, s_cut, urep_cut
118 REAL(kind=dp), DIMENSION(0:3) :: eta_a, eta_b, skself
119 REAL(kind=dp), DIMENSION(10) :: urep
120 REAL(kind=dp), DIMENSION(2) :: surr
121 REAL(kind=dp), DIMENSION(3) :: drij, force_ab, force_rr, force_w, rij, &
122 srep
123 REAL(kind=dp), DIMENSION(:, :), POINTER :: dfblock, dsblock, fblock, fmatij, &
124 fmatji, pblock, sblock, scoeff, &
125 smatij, smatji, spxr, wblock
126 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
127 TYPE(atprop_type), POINTER :: atprop
128 TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
129 TYPE(cp_logger_type), POINTER :: logger
130 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
131 TYPE(dft_control_type), POINTER :: dft_control
132 TYPE(dftb_control_type), POINTER :: dftb_control
133 TYPE(kpoint_type), POINTER :: kpoints
135 DIMENSION(:), POINTER :: nl_iterator
136 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
137 POINTER :: sab_orb
138 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
139 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a, dftb_kind_b
140 TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
141 POINTER :: dftb_potential
142 TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji
143 TYPE(qs_energy_type), POINTER :: energy
144 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
145 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
146 TYPE(qs_ks_env_type), POINTER :: ks_env
147 TYPE(qs_rho_type), POINTER :: rho
148 TYPE(virial_type), POINTER :: virial
149
150 CALL timeset(routinen, handle)
151
152 ! set pointers
153 iptr = 0
154 DO la = 0, 3
155 DO lb = 0, 3
156 llm = 0
157 DO l1 = 0, max(la, lb)
158 DO l2 = 0, min(l1, la, lb)
159 DO m = 0, l2
160 llm = llm + 1
161 iptr(l1, l2, m, la, lb) = llm
162 END DO
163 END DO
164 END DO
165 END DO
166 END DO
167
168 NULLIFY (logger, virial, atprop)
169 logger => cp_get_default_logger()
170
171 NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
172 qs_kind_set, sab_orb, ks_env)
173
174 CALL get_qs_env(qs_env=qs_env, &
175 energy=energy, &
176 atomic_kind_set=atomic_kind_set, &
177 qs_kind_set=qs_kind_set, &
178 matrix_h_kp=matrix_h, &
179 matrix_s_kp=matrix_s, &
180 atprop=atprop, &
181 dft_control=dft_control, &
182 ks_env=ks_env)
183
184 dftb_control => dft_control%qs_control%dftb_control
185 nimg = dft_control%nimages
186 ! Allocate the overlap and Hamiltonian matrix
187 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
188 nderivatives = 0
189 IF (dftb_control%self_consistent .AND. calculate_forces) nderivatives = 1
190 CALL setup_matrices2(qs_env, nderivatives, nimg, matrix_s, "OVERLAP", sab_orb)
191 CALL setup_matrices2(qs_env, 0, nimg, matrix_h, "CORE HAMILTONIAN", sab_orb)
192 CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
193 CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
194
195 NULLIFY (dftb_potential)
196 CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential)
197 NULLIFY (particle_set)
198 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
199
200 IF (calculate_forces) THEN
201 NULLIFY (rho, force, matrix_w)
202 CALL get_qs_env(qs_env=qs_env, &
203 rho=rho, &
204 matrix_w_kp=matrix_w, &
205 virial=virial, &
206 force=force)
207 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
208
209 IF (SIZE(matrix_p, 1) == 2) THEN
210 DO img = 1, nimg
211 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
212 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
213 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
214 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
215 END DO
216 END IF
217 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
218 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
219 END IF
220 ! atomic energy decomposition
221 IF (atprop%energy) THEN
222 CALL atprop_array_init(atprop%atecc, natom=SIZE(particle_set))
223 END IF
224
225 NULLIFY (cell_to_index)
226 IF (nimg > 1) THEN
227 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
228 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
229 END IF
230
231 erep = 0._dp
232
233 nkind = SIZE(atomic_kind_set)
234
235 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
236 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
237 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
238 iatom=iatom, jatom=jatom, r=rij, cell=cell)
239 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
240 CALL get_dftb_atom_param(dftb_kind_a, &
241 defined=defined, lmax=lmaxi, skself=skself, &
242 eta=eta_a, natorb=natorb_a)
243 IF (.NOT. defined .OR. natorb_a < 1) cycle
244 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
245 CALL get_dftb_atom_param(dftb_kind_b, &
246 defined=defined, lmax=lmaxj, eta=eta_b, natorb=natorb_b)
247
248 IF (.NOT. defined .OR. natorb_b < 1) cycle
249
250 ! retrieve information on F and S matrix
251 dftb_param_ij => dftb_potential(ikind, jkind)
252 dftb_param_ji => dftb_potential(jkind, ikind)
253 ! assume table size and type is symmetric
254 ngrd = dftb_param_ij%ngrd
255 ngrdcut = dftb_param_ij%ngrdcut
256 dgrd = dftb_param_ij%dgrd
257 ddr = dgrd*dftb_fd_deriv_step
258 cpassert(dftb_param_ij%llm == dftb_param_ji%llm)
259 llm = dftb_param_ij%llm
260 fmatij => dftb_param_ij%fmat
261 smatij => dftb_param_ij%smat
262 fmatji => dftb_param_ji%fmat
263 smatji => dftb_param_ji%smat
264 ! repulsive pair potential
265 n_urpoly = dftb_param_ij%n_urpoly
266 urep_cut = dftb_param_ij%urep_cut
267 urep = dftb_param_ij%urep
268 spxr => dftb_param_ij%spxr
269 scoeff => dftb_param_ij%scoeff
270 spdim = dftb_param_ij%spdim
271 s_cut = dftb_param_ij%s_cut
272 srep = dftb_param_ij%srep
273 surr = dftb_param_ij%surr
274
275 dr = sqrt(sum(rij(:)**2))
276 IF (nint(dr/dgrd) <= ngrdcut) THEN
277
278 IF (nimg == 1) THEN
279 ic = 1
280 ELSE
281 ic = cell_to_index(cell(1), cell(2), cell(3))
282 cpassert(ic > 0)
283 END IF
284
285 icol = max(iatom, jatom)
286 irow = min(iatom, jatom)
287 NULLIFY (sblock, fblock)
288 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
289 row=irow, col=icol, block=sblock, found=found)
290 cpassert(found)
291 CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
292 row=irow, col=icol, block=fblock, found=found)
293 cpassert(found)
294
295 IF (calculate_forces) THEN
296 NULLIFY (pblock)
297 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
298 row=irow, col=icol, block=pblock, found=found)
299 cpassert(ASSOCIATED(pblock))
300 NULLIFY (wblock)
301 CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
302 row=irow, col=icol, block=wblock, found=found)
303 cpassert(ASSOCIATED(wblock))
304 IF (dftb_control%self_consistent) THEN
305 DO i = 2, 4
306 NULLIFY (dsblocks(i)%block)
307 CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
308 row=irow, col=icol, block=dsblocks(i)%block, found=found)
309 cpassert(found)
310 END DO
311 END IF
312 END IF
313
314 IF (iatom == jatom .AND. dr < 0.001_dp) THEN
315 ! diagonal block
316 DO i = 1, natorb_a
317 sblock(i, i) = sblock(i, i) + 1._dp
318 fblock(i, i) = fblock(i, i) + skself(orbptr(i))
319 END DO
320 ELSE
321 ! off-diagonal block
322 CALL compute_block_sk(sblock, smatij, smatji, rij, ngrd, ngrdcut, dgrd, &
323 llm, lmaxi, lmaxj, irow, iatom)
324 CALL compute_block_sk(fblock, fmatij, fmatji, rij, ngrd, ngrdcut, dgrd, &
325 llm, lmaxi, lmaxj, irow, iatom)
326 IF (calculate_forces) THEN
327 force_ab = 0._dp
328 force_w = 0._dp
329 n1 = SIZE(fblock, 1)
330 n2 = SIZE(fblock, 2)
331 ! make sure that displacement is in the correct direction depending on the position
332 ! of the block (upper or lower triangle)
333 f0 = 1.0_dp
334 IF (irow == iatom) f0 = -1.0_dp
335
336 ALLOCATE (dfblock(n1, n2), dsblock(n1, n2))
337
338 DO i = 1, 3
339 drij = rij
340 dfblock = 0._dp; dsblock = 0._dp
341
342 drij(i) = rij(i) - ddr*f0
343 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
344 llm, lmaxi, lmaxj, irow, iatom)
345 CALL compute_block_sk(dfblock, fmatij, fmatji, drij, ngrd, ngrdcut, dgrd, &
346 llm, lmaxi, lmaxj, irow, iatom)
347
348 dsblock = -dsblock
349 dfblock = -dfblock
350
351 drij(i) = rij(i) + ddr*f0
352 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
353 llm, lmaxi, lmaxj, irow, iatom)
354 CALL compute_block_sk(dfblock, fmatij, fmatji, drij, ngrd, ngrdcut, dgrd, &
355 llm, lmaxi, lmaxj, irow, iatom)
356
357 dfblock = dfblock/(2.0_dp*ddr)
358 dsblock = dsblock/(2.0_dp*ddr)
359
360 foab = 2.0_dp*sum(dfblock*pblock)
361 fow = -2.0_dp*sum(dsblock*wblock)
362
363 force_ab(i) = force_ab(i) + foab
364 force_w(i) = force_w(i) + fow
365 IF (dftb_control%self_consistent) THEN
366 cpassert(ASSOCIATED(dsblocks(i + 1)%block))
367 dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock
368 END IF
369 END DO
370 IF (use_virial) THEN
371 IF (iatom == jatom) f0 = 0.5_dp*f0
372 CALL virial_pair_force(virial%pv_virial, -f0, force_ab, rij)
373 CALL virial_pair_force(virial%pv_virial, -f0, force_w, rij)
374 END IF
375 DEALLOCATE (dfblock, dsblock)
376 END IF
377 END IF
378
379 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
380 atom_a = atom_of_kind(iatom)
381 atom_b = atom_of_kind(jatom)
382 IF (irow == iatom) force_ab = -force_ab
383 IF (irow == iatom) force_w = -force_w
384 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
385 force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
386 force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - force_w(:)
387 force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + force_w(:)
388 END IF
389
390 END IF
391
392 ! repulsive potential
393 IF ((dr <= urep_cut .OR. spdim > 0) .AND. dr > 0.001_dp) THEN
394 erepij = 0._dp
395 CALL urep_egr(rij, dr, erepij, force_rr, &
396 n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, calculate_forces)
397 erep = erep + erepij
398 IF (atprop%energy) THEN
399 atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*erepij
400 atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*erepij
401 END IF
402 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
403 atom_a = atom_of_kind(iatom)
404 atom_b = atom_of_kind(jatom)
405 force(ikind)%repulsive(:, atom_a) = &
406 force(ikind)%repulsive(:, atom_a) - force_rr(:)
407 force(jkind)%repulsive(:, atom_b) = &
408 force(jkind)%repulsive(:, atom_b) + force_rr(:)
409 IF (use_virial) THEN
410 f0 = -1.0_dp
411 IF (iatom == jatom) f0 = -0.5_dp
412 CALL virial_pair_force(virial%pv_virial, f0, force_rr, rij)
413 END IF
414 END IF
415 END IF
416
417 END DO
418 CALL neighbor_list_iterator_release(nl_iterator)
419
420 DO i = 1, SIZE(matrix_s, 1)
421 DO img = 1, nimg
422 CALL dbcsr_finalize(matrix_s(i, img)%matrix)
423 END DO
424 END DO
425 DO i = 1, SIZE(matrix_h, 1)
426 DO img = 1, nimg
427 CALL dbcsr_finalize(matrix_h(i, img)%matrix)
428 END DO
429 END DO
430
431 ! set repulsive energy
432 CALL para_env%sum(erep)
433 energy%repulsive = erep
434
435 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
436 IF (btest(cp_print_key_should_output(logger%iter_info, &
437 qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
438 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
439 extension=".Log")
440 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
441 after = min(max(after, 1), 16)
442 DO img = 1, nimg
443 CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, &
444 output_unit=iw, omit_headers=omit_headers)
445 END DO
446
447 CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
448 "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
449 END IF
450
451 IF (btest(cp_print_key_should_output(logger%iter_info, &
452 qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
453 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
454 extension=".Log")
455 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
456 after = min(max(after, 1), 16)
457 DO img = 1, nimg
458 CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, &
459 output_unit=iw, omit_headers=omit_headers)
460
461 IF (btest(cp_print_key_should_output(logger%iter_info, &
462 qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
463 DO i = 2, SIZE(matrix_s, 1)
464 CALL cp_dbcsr_write_sparse_matrix(matrix_s(i, img)%matrix, 4, after, qs_env, para_env, &
465 output_unit=iw, omit_headers=omit_headers)
466 END DO
467 END IF
468 END DO
469
470 CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
471 "DFT%PRINT%AO_MATRICES/OVERLAP")
472 END IF
473
474 IF (calculate_forces) THEN
475 IF (SIZE(matrix_p, 1) == 2) THEN
476 DO img = 1, nimg
477 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
478 beta_scalar=-1.0_dp)
479 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, alpha_scalar=1.0_dp, &
480 beta_scalar=-1.0_dp)
481 END DO
482 END IF
483 END IF
484
485 CALL timestop(handle)
486
487 END SUBROUTINE build_dftb_matrices
488
489! **************************************************************************************************
490!> \brief ...
491!> \param qs_env ...
492!> \param calculate_forces ...
493!> \param just_energy ...
494! **************************************************************************************************
495 SUBROUTINE build_dftb_ks_matrix(qs_env, calculate_forces, just_energy)
496 TYPE(qs_environment_type), POINTER :: qs_env
497 LOGICAL, INTENT(in) :: calculate_forces, just_energy
498
499 CHARACTER(len=*), PARAMETER :: routinen = 'build_dftb_ks_matrix'
500
501 INTEGER :: atom_a, handle, iatom, ikind, img, &
502 ispin, natom, nkind, nspins, &
503 output_unit
504 REAL(kind=dp) :: pc_ener, qmmm_el, zeff
505 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mix_charge
506 REAL(kind=dp), DIMENSION(:), POINTER :: mcharge, occupation_numbers
507 REAL(kind=dp), DIMENSION(:, :), POINTER :: charges
508 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
509 TYPE(cp_logger_type), POINTER :: logger
510 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs
511 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_h, matrix_p, matrix_s
512 TYPE(dbcsr_type), POINTER :: mo_coeff
513 TYPE(dft_control_type), POINTER :: dft_control
514 TYPE(mp_para_env_type), POINTER :: para_env
515 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
516 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
517 TYPE(qs_energy_type), POINTER :: energy
518 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
519 TYPE(qs_ks_env_type), POINTER :: ks_env
520 TYPE(qs_rho_type), POINTER :: rho
521 TYPE(qs_scf_env_type), POINTER :: scf_env
522 TYPE(section_vals_type), POINTER :: scf_section
523
524 CALL timeset(routinen, handle)
525 NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, &
526 ks_matrix, rho, energy, scf_env)
527 logger => cp_get_default_logger()
528 cpassert(ASSOCIATED(qs_env))
529
530 CALL get_qs_env(qs_env, &
531 dft_control=dft_control, &
532 atomic_kind_set=atomic_kind_set, &
533 qs_kind_set=qs_kind_set, &
534 matrix_h_kp=matrix_h, &
535 para_env=para_env, &
536 ks_env=ks_env, &
537 matrix_ks_kp=ks_matrix, &
538 rho=rho, &
539 energy=energy)
540
541 energy%hartree = 0.0_dp
542 energy%qmmm_el = 0.0_dp
543
544 scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
545 nspins = dft_control%nspins
546 cpassert(ASSOCIATED(matrix_h))
547 cpassert(ASSOCIATED(rho))
548 cpassert(SIZE(ks_matrix) > 0)
549
550 DO ispin = 1, nspins
551 DO img = 1, SIZE(ks_matrix, 2)
552 ! copy the core matrix into the fock matrix
553 CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
554 END DO
555 END DO
556
557 IF (dft_control%qs_control%dftb_control%self_consistent) THEN
558 ! Mulliken charges
559 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
560 matrix_s_kp=matrix_s)
561 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
562 natom = SIZE(particle_set)
563 ALLOCATE (charges(natom, nspins))
564 !
565 CALL mulliken_charges(matrix_p, matrix_s, para_env, charges)
566 !
567 ALLOCATE (mcharge(natom))
568 nkind = SIZE(atomic_kind_set)
569 DO ikind = 1, nkind
570 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
571 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
572 CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
573 DO iatom = 1, natom
574 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
575 mcharge(atom_a) = zeff - sum(charges(atom_a, 1:nspins))
576 END DO
577 END DO
578 DEALLOCATE (charges)
579
580 IF ((.NOT. dft_control%qs_control%do_ls_scf) .AND. &
581 (dft_control%qs_control%dftb_control%tblite_scc_mixer /= tblite_scc_mixer_auto)) THEN
582 CALL get_qs_env(qs_env=qs_env, scf_env=scf_env)
583 ALLOCATE (mix_charge(SIZE(mcharge), 1))
584 mix_charge(:, 1) = mcharge(:)
585 CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
586 mix_charge, para_env, scf_env%iter_count, &
587 scc_mixer=dft_control%qs_control%dftb_control%tblite_scc_mixer, &
588 tblite_mixer_iterations= &
589 dft_control%qs_control%dftb_control%tblite_mixer_iterations, &
590 tblite_mixer_damping=dft_control%qs_control%dftb_control%tblite_mixer_damping, &
591 tblite_mixer_memory=dft_control%qs_control%dftb_control%tblite_mixer_memory, &
592 tblite_mixer_omega0=dft_control%qs_control%dftb_control%tblite_mixer_omega0, &
593 tblite_mixer_min_weight= &
594 dft_control%qs_control%dftb_control%tblite_mixer_min_weight, &
595 tblite_mixer_max_weight= &
596 dft_control%qs_control%dftb_control%tblite_mixer_max_weight, &
597 tblite_mixer_weight_factor= &
598 dft_control%qs_control%dftb_control%tblite_mixer_weight_factor)
599 mcharge(:) = mix_charge(:, 1)
600 DEALLOCATE (mix_charge)
601 END IF
602
603 CALL build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, &
604 calculate_forces, just_energy)
605
606 CALL efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, &
607 calculate_forces, just_energy)
608
609 DEALLOCATE (mcharge)
610
611 END IF
612
613 IF (qs_env%qmmm) THEN
614 cpassert(SIZE(ks_matrix, 2) == 1)
615 DO ispin = 1, nspins
616 ! If QM/MM sumup the 1el Hamiltonian
617 CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
618 1.0_dp, 1.0_dp)
619 CALL qs_rho_get(rho, rho_ao=matrix_p1)
620 ! Compute QM/MM Energy
621 CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
622 matrix_p1(ispin)%matrix, qmmm_el)
623 energy%qmmm_el = energy%qmmm_el + qmmm_el
624 END DO
625 pc_ener = qs_env%ks_qmmm_env%pc_ener
626 energy%qmmm_el = energy%qmmm_el + pc_ener
627 END IF
628
629 energy%total = energy%core + energy%hartree + energy%qmmm_el + energy%efield + &
630 energy%repulsive + energy%dispersion + energy%dftb3
631
632 IF (dft_control%qs_control%dftb_control%self_consistent) THEN
633 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
634 extension=".scfLog")
635 IF (output_unit > 0) THEN
636 WRITE (unit=output_unit, fmt="(/,(T9,A,T60,F20.10))") &
637 "Repulsive pair potential energy: ", energy%repulsive, &
638 "Zeroth order Hamiltonian energy: ", energy%core, &
639 "Charge fluctuation energy: ", energy%hartree, &
640 "London dispersion energy: ", energy%dispersion
641 IF (abs(energy%efield) > 1.e-12_dp) THEN
642 WRITE (unit=output_unit, fmt="(T9,A,T60,F20.10)") &
643 "Electric field interaction energy: ", energy%efield
644 END IF
645 IF (dft_control%qs_control%dftb_control%dftb3_diagonal) THEN
646 WRITE (unit=output_unit, fmt="(T9,A,T60,F20.10)") &
647 "DFTB3 3rd Order Energy Correction ", energy%dftb3
648 END IF
649 IF (qs_env%qmmm) THEN
650 WRITE (unit=output_unit, fmt="(T9,A,T60,F20.10)") &
651 "QM/MM Electrostatic energy: ", energy%qmmm_el
652 END IF
653 END IF
654 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
655 "PRINT%DETAILED_ENERGY")
656 END IF
657 ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C (plus occupation numbers)
658 IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
659 cpassert(SIZE(ks_matrix, 2) == 1)
660 block
661 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
662 CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
663 DO ispin = 1, SIZE(mo_derivs)
664 CALL get_mo_set(mo_set=mo_array(ispin), &
665 mo_coeff_b=mo_coeff, occupation_numbers=occupation_numbers)
666 IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
667 cpabort("")
668 END IF
669 CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
670 0.0_dp, mo_derivs(ispin)%matrix)
671 END DO
672 END block
673 END IF
674
675 CALL timestop(handle)
676
677 END SUBROUTINE build_dftb_ks_matrix
678
679! **************************************************************************************************
680!> \brief ...
681!> \param qs_env ...
682!> \param nderivative ...
683!> \param matrix_s ...
684! **************************************************************************************************
685 SUBROUTINE build_dftb_overlap(qs_env, nderivative, matrix_s)
686
687 TYPE(qs_environment_type), POINTER :: qs_env
688 INTEGER, INTENT(IN) :: nderivative
689 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
690
691 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_dftb_overlap'
692
693 INTEGER :: handle, i, iatom, icol, ikind, indder, irow, j, jatom, jkind, l1, l2, la, lb, &
694 llm, lmaxi, lmaxj, m, n1, n2, natom, natorb_a, natorb_b, ngrd, ngrdcut, nkind
695 LOGICAL :: defined, found
696 REAL(kind=dp) :: ddr, dgrd, dr, f0
697 REAL(kind=dp), DIMENSION(0:3) :: skself
698 REAL(kind=dp), DIMENSION(3) :: drij, rij
699 REAL(kind=dp), DIMENSION(:, :), POINTER :: dsblock, dsblockm, sblock, smatij, smatji
700 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dsblock1
701 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
702 TYPE(block_p_type), DIMENSION(2:10) :: dsblocks
703 TYPE(cp_logger_type), POINTER :: logger
704 TYPE(dft_control_type), POINTER :: dft_control
705 TYPE(dftb_control_type), POINTER :: dftb_control
706 TYPE(neighbor_list_iterator_p_type), &
707 DIMENSION(:), POINTER :: nl_iterator
708 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
709 POINTER :: sab_orb
710 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a, dftb_kind_b
711 TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
712 POINTER :: dftb_potential
713 TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji
714 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
715
716 CALL timeset(routinen, handle)
717
718 ! set pointers
719 iptr = 0
720 DO la = 0, 3
721 DO lb = 0, 3
722 llm = 0
723 DO l1 = 0, max(la, lb)
724 DO l2 = 0, min(l1, la, lb)
725 DO m = 0, l2
726 llm = llm + 1
727 iptr(l1, l2, m, la, lb) = llm
728 END DO
729 END DO
730 END DO
731 END DO
732 END DO
733
734 NULLIFY (logger)
735 logger => cp_get_default_logger()
736
737 NULLIFY (atomic_kind_set, qs_kind_set, sab_orb)
738
739 CALL get_qs_env(qs_env=qs_env, &
740 atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
741 dft_control=dft_control)
742
743 dftb_control => dft_control%qs_control%dftb_control
744
745 NULLIFY (dftb_potential)
746 CALL get_qs_env(qs_env=qs_env, &
747 dftb_potential=dftb_potential)
748
749 nkind = SIZE(atomic_kind_set)
750
751 ! Allocate the overlap matrix
752 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
753 CALL setup_matrices1(qs_env, nderivative, matrix_s, 'OVERLAP', sab_orb)
754
755 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
756 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
757 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
758 iatom=iatom, jatom=jatom, r=rij)
759
760 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
761 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
762 CALL get_dftb_atom_param(dftb_kind_a, &
763 defined=defined, lmax=lmaxi, skself=skself, &
764 natorb=natorb_a)
765
766 IF (.NOT. defined .OR. natorb_a < 1) cycle
767
768 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
769 CALL get_dftb_atom_param(dftb_kind_b, &
770 defined=defined, lmax=lmaxj, natorb=natorb_b)
771
772 IF (.NOT. defined .OR. natorb_b < 1) cycle
773
774 ! retrieve information on F and S matrix
775 dftb_param_ij => dftb_potential(ikind, jkind)
776 dftb_param_ji => dftb_potential(jkind, ikind)
777 ! assume table size and type is symmetric
778 ngrd = dftb_param_ij%ngrd
779 ngrdcut = dftb_param_ij%ngrdcut
780 dgrd = dftb_param_ij%dgrd
781 ddr = dgrd*dftb_fd_deriv_step
782 cpassert(dftb_param_ij%llm == dftb_param_ji%llm)
783 llm = dftb_param_ij%llm
784 smatij => dftb_param_ij%smat
785 smatji => dftb_param_ji%smat
786
787 dr = sqrt(sum(rij(:)**2))
788 IF (nint(dr/dgrd) <= ngrdcut) THEN
789
790 icol = max(iatom, jatom); irow = min(iatom, jatom)
791
792 NULLIFY (sblock)
793 CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
794 row=irow, col=icol, block=sblock, found=found)
795 cpassert(found)
796
797 IF (nderivative > 0) THEN
798 DO i = 2, SIZE(matrix_s, 1)
799 NULLIFY (dsblocks(i)%block)
800 CALL dbcsr_get_block_p(matrix=matrix_s(i)%matrix, &
801 row=irow, col=icol, block=dsblocks(i)%block, found=found)
802 END DO
803 END IF
804
805 IF (iatom == jatom .AND. dr < 0.001_dp) THEN
806 ! diagonal block
807 DO i = 1, natorb_a
808 sblock(i, i) = sblock(i, i) + 1._dp
809 END DO
810 ELSE
811 ! off-diagonal block
812 CALL compute_block_sk(sblock, smatij, smatji, rij, ngrd, ngrdcut, dgrd, &
813 llm, lmaxi, lmaxj, irow, iatom)
814
815 IF (nderivative >= 1) THEN
816 n1 = SIZE(sblock, 1); n2 = SIZE(sblock, 2)
817 indder = 1 ! used to put the 2nd derivatives in the correct matric (5=xx,8=yy,10=zz)
818
819 ALLOCATE (dsblock1(n1, n2, 3), dsblock(n1, n2), dsblockm(n1, n2))
820 dsblock1 = 0.0_dp
821 DO i = 1, 3
822 dsblock = 0._dp; dsblockm = 0.0_dp
823 drij = rij
824 f0 = 1.0_dp; IF (irow == iatom) f0 = -1.0_dp
825
826 drij(i) = rij(i) - ddr*f0
827 CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
828 llm, lmaxi, lmaxj, irow, iatom)
829
830 drij(i) = rij(i) + ddr*f0
831 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
832 llm, lmaxi, lmaxj, irow, iatom)
833
834 dsblock1(:, :, i) = (dsblock + dsblockm)
835 dsblock = dsblock - dsblockm
836 dsblock = dsblock/(2.0_dp*ddr)
837
838 cpassert(ASSOCIATED(dsblocks(i + 1)%block))
839 dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock
840 IF (nderivative > 1) THEN
841 indder = indder + 5 - i
842 dsblocks(indder)%block = 0.0_dp
843 dsblocks(indder)%block = dsblocks(indder)%block + &
844 (dsblock1(:, :, i) - 2.0_dp*sblock)/ddr**2
845 END IF
846 END DO
847
848 IF (nderivative > 1) THEN
849 DO i = 1, 2
850 DO j = i + 1, 3
851 dsblock = 0._dp; dsblockm = 0.0_dp
852 drij = rij
853 f0 = 1.0_dp; IF (irow == iatom) f0 = -1.0_dp
854
855 drij(i) = rij(i) - ddr*f0; drij(j) = rij(j) - ddr*f0
856 CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
857 llm, lmaxi, lmaxj, irow, iatom)
858
859 drij(i) = rij(i) + ddr*f0; drij(j) = rij(j) + ddr*f0
860 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
861 llm, lmaxi, lmaxj, irow, iatom)
862
863 indder = 2 + 2*i + j
864 dsblocks(indder)%block = 0.0_dp
865 dsblocks(indder)%block = &
866 dsblocks(indder)%block + ( &
867 dsblock + dsblockm - dsblock1(:, :, i) - dsblock1(:, :, j) + 2.0_dp*sblock)/(2.0_dp*ddr**2)
868 END DO
869 END DO
870 END IF
871
872 DEALLOCATE (dsblock1, dsblock, dsblockm)
873 END IF
874 END IF
875 END IF
876 END DO
877 CALL neighbor_list_iterator_release(nl_iterator)
878
879 DO i = 1, SIZE(matrix_s, 1)
880 CALL dbcsr_finalize(matrix_s(i)%matrix)
881 END DO
882
883 CALL timestop(handle)
884
885 END SUBROUTINE build_dftb_overlap
886
887! **************************************************************************************************
888!> \brief ...
889!> \param qs_env ...
890!> \param nderivative ...
891!> \param matrices ...
892!> \param mnames ...
893!> \param sab_nl ...
894! **************************************************************************************************
895 SUBROUTINE setup_matrices1(qs_env, nderivative, matrices, mnames, sab_nl)
896
897 TYPE(qs_environment_type), POINTER :: qs_env
898 INTEGER, INTENT(IN) :: nderivative
899 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
900 CHARACTER(LEN=*) :: mnames
901 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
902 POINTER :: sab_nl
903
904 CHARACTER(1) :: symmetry_type
905 CHARACTER(LEN=default_string_length) :: matnames
906 INTEGER :: i, natom, neighbor_list_id, nkind, nmat, &
907 nsgf
908 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
909 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
910 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
911 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
912 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
913 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
914
915 NULLIFY (particle_set, atomic_kind_set)
916
917 CALL get_qs_env(qs_env=qs_env, &
918 atomic_kind_set=atomic_kind_set, &
919 qs_kind_set=qs_kind_set, &
920 particle_set=particle_set, &
921 dbcsr_dist=dbcsr_dist, &
922 neighbor_list_id=neighbor_list_id)
923
924 nkind = SIZE(atomic_kind_set)
925 natom = SIZE(particle_set)
926
927 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
928
929 ALLOCATE (first_sgf(natom))
930 ALLOCATE (last_sgf(natom))
931
932 CALL get_particle_set(particle_set, qs_kind_set, &
933 first_sgf=first_sgf, &
934 last_sgf=last_sgf)
935
936 nmat = 0
937 IF (nderivative == 0) nmat = 1
938 IF (nderivative == 1) nmat = 4
939 IF (nderivative == 2) nmat = 10
940 cpassert(nmat > 0)
941
942 ALLOCATE (row_blk_sizes(natom))
943 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
944
945 CALL dbcsr_allocate_matrix_set(matrices, nmat)
946
947 ! Up to 2nd derivative take care to get the symmetries correct
948 DO i = 1, nmat
949 IF (i > 1) THEN
950 matnames = trim(mnames)//" DERIVATIVE MATRIX DFTB"
951 symmetry_type = dbcsr_type_antisymmetric
952 IF (i > 4) symmetry_type = dbcsr_type_symmetric
953 ELSE
954 symmetry_type = dbcsr_type_symmetric
955 matnames = trim(mnames)//" MATRIX DFTB"
956 END IF
957 ALLOCATE (matrices(i)%matrix)
958 CALL dbcsr_create(matrix=matrices(i)%matrix, &
959 name=trim(matnames), &
960 dist=dbcsr_dist, matrix_type=symmetry_type, &
961 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
962 mutable_work=.true.)
963 CALL cp_dbcsr_alloc_block_from_nbl(matrices(i)%matrix, sab_nl)
964 END DO
965
966 DEALLOCATE (first_sgf)
967 DEALLOCATE (last_sgf)
968
969 DEALLOCATE (row_blk_sizes)
970
971 END SUBROUTINE setup_matrices1
972
973! **************************************************************************************************
974!> \brief ...
975!> \param qs_env ...
976!> \param nderivative ...
977!> \param nimg ...
978!> \param matrices ...
979!> \param mnames ...
980!> \param sab_nl ...
981! **************************************************************************************************
982 SUBROUTINE setup_matrices2(qs_env, nderivative, nimg, matrices, mnames, sab_nl)
983
984 TYPE(qs_environment_type), POINTER :: qs_env
985 INTEGER, INTENT(IN) :: nderivative, nimg
986 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrices
987 CHARACTER(LEN=*) :: mnames
988 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
989 POINTER :: sab_nl
990
991 CHARACTER(1) :: symmetry_type
992 CHARACTER(LEN=default_string_length) :: matnames
993 INTEGER :: i, img, natom, neighbor_list_id, nkind, &
994 nmat, nsgf
995 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
996 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
997 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
998 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
999 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1000 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1001
1002 NULLIFY (particle_set, atomic_kind_set)
1003
1004 CALL get_qs_env(qs_env=qs_env, &
1005 atomic_kind_set=atomic_kind_set, &
1006 qs_kind_set=qs_kind_set, &
1007 particle_set=particle_set, &
1008 dbcsr_dist=dbcsr_dist, &
1009 neighbor_list_id=neighbor_list_id)
1010
1011 nkind = SIZE(atomic_kind_set)
1012 natom = SIZE(particle_set)
1013
1014 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
1015
1016 ALLOCATE (first_sgf(natom))
1017 ALLOCATE (last_sgf(natom))
1018
1019 CALL get_particle_set(particle_set, qs_kind_set, &
1020 first_sgf=first_sgf, &
1021 last_sgf=last_sgf)
1022
1023 nmat = 0
1024 IF (nderivative == 0) nmat = 1
1025 IF (nderivative == 1) nmat = 4
1026 IF (nderivative == 2) nmat = 10
1027 cpassert(nmat > 0)
1028
1029 ALLOCATE (row_blk_sizes(natom))
1030 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
1031
1032 CALL dbcsr_allocate_matrix_set(matrices, nmat, nimg)
1033
1034 ! Up to 2nd derivative take care to get the symmetries correct
1035 DO img = 1, nimg
1036 DO i = 1, nmat
1037 IF (i > 1) THEN
1038 matnames = trim(mnames)//" DERIVATIVE MATRIX DFTB"
1039 symmetry_type = dbcsr_type_antisymmetric
1040 IF (i > 4) symmetry_type = dbcsr_type_symmetric
1041 ELSE
1042 symmetry_type = dbcsr_type_symmetric
1043 matnames = trim(mnames)//" MATRIX DFTB"
1044 END IF
1045 ALLOCATE (matrices(i, img)%matrix)
1046 CALL dbcsr_create(matrix=matrices(i, img)%matrix, &
1047 name=trim(matnames), &
1048 dist=dbcsr_dist, matrix_type=symmetry_type, &
1049 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
1050 mutable_work=.true.)
1051 CALL cp_dbcsr_alloc_block_from_nbl(matrices(i, img)%matrix, sab_nl)
1052 END DO
1053 END DO
1054
1055 DEALLOCATE (first_sgf)
1056 DEALLOCATE (last_sgf)
1057
1058 DEALLOCATE (row_blk_sizes)
1059
1060 END SUBROUTINE setup_matrices2
1061
1062END MODULE qs_dftb_matrices
Define the atomic kind types and their sub types.
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.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Holds information on atomic properties.
subroutine, public atprop_array_init(atarray, natom)
...
collect pointers to a block of reals
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_finalize(matrix)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
DBCSR operations in CP2K.
DBCSR output in CP2K.
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)
...
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
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)
...
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,...
integer, parameter, public cp_p_file
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...
Calculation of electric field contributions in TB.
subroutine, public efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public tblite_scc_mixer_auto
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
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
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Types and basic routines needed for a kpoint calculation.
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.
Interface to the message passing library MPI.
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition mulliken.F:13
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
subroutine, public charge_mixing(mixing_method, mixing_store, charges, para_env, iter_count, scc_mixer, tblite_mixer_iterations, tblite_mixer_damping, tblite_mixer_memory, tblite_mixer_omega0, tblite_mixer_min_weight, tblite_mixer_max_weight, tblite_mixer_weight_factor)
Driver for TB SCC variable mixing, calls the requested method.
Calculation of Coulomb contributions in DFTB.
subroutine, public build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
...
Calculation of Overlap and Hamiltonian matrices in DFTB.
subroutine, public build_dftb_ks_matrix(qs_env, calculate_forces, just_energy)
...
integer, dimension(16), parameter orbptr
subroutine, public build_dftb_matrices(qs_env, para_env, calculate_forces)
...
subroutine, public build_dftb_overlap(qs_env, nderivative, matrix_s)
...
Definition of the DFTB parameter types.
Working with the DFTB parameter types.
subroutine, public urep_egr(rv, r, erep, derep, n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, dograd)
...
subroutine, public compute_block_sk(block, smatij, smatji, rij, ngrd, ngrdcut, dgrd, llm, lmaxi, lmaxj, irow, iatom)
...
integer, dimension(0:3, 0:3, 0:3, 0:3, 0:3), public iptr
subroutine, public get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
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.
Define the quickstep kind type and their sub types.
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.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
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)
...
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 and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
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...
module that contains the definitions of the scf types
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Provides all information about an atomic kind.
type for the atomic properties
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.