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