Line data Source code
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 : ! **************************************************************************************************
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_api, ONLY: &
22 : dbcsr_add, dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, &
23 : dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_block_p, dbcsr_multiply, dbcsr_p_type, &
24 : dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_symmetric
25 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
26 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
27 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
28 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
29 : USE cp_log_handling, ONLY: cp_get_default_logger,&
30 : cp_logger_type
31 : USE cp_output_handling, ONLY: cp_p_file,&
32 : cp_print_key_finished_output,&
33 : cp_print_key_should_output,&
34 : cp_print_key_unit_nr
35 : USE efield_tb_methods, ONLY: efield_tb_matrix
36 : USE input_constants, ONLY: tblite_scc_mixer_auto
37 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
38 : section_vals_type,&
39 : section_vals_val_get
40 : USE kinds, ONLY: default_string_length,&
41 : dp
42 : USE kpoint_types, ONLY: get_kpoint_info,&
43 : kpoint_type
44 : USE message_passing, ONLY: mp_para_env_type
45 : USE mulliken, ONLY: mulliken_charges
46 : USE particle_methods, ONLY: get_particle_set
47 : USE particle_types, ONLY: particle_type
48 : USE qs_charge_mixing, ONLY: charge_mixing
49 : USE qs_dftb_coulomb, ONLY: build_dftb_coulomb
50 : USE qs_dftb_types, ONLY: qs_dftb_atom_type,&
51 : qs_dftb_pairpot_type
52 : USE qs_dftb_utils, ONLY: compute_block_sk,&
53 : get_dftb_atom_param,&
54 : iptr,&
55 : urep_egr
56 : USE qs_energy_types, ONLY: qs_energy_type
57 : USE qs_environment_types, ONLY: get_qs_env,&
58 : qs_environment_type
59 : USE qs_force_types, ONLY: qs_force_type
60 : USE qs_kind_types, ONLY: get_qs_kind,&
61 : get_qs_kind_set,&
62 : qs_kind_type
63 : USE qs_ks_types, ONLY: get_ks_env,&
64 : qs_ks_env_type,&
65 : set_ks_env
66 : USE qs_mo_types, ONLY: get_mo_set,&
67 : mo_set_type
68 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
69 : neighbor_list_iterate,&
70 : neighbor_list_iterator_create,&
71 : neighbor_list_iterator_p_type,&
72 : neighbor_list_iterator_release,&
73 : neighbor_list_set_p_type
74 : USE qs_rho_types, ONLY: qs_rho_get,&
75 : qs_rho_type
76 : USE qs_scf_types, ONLY: qs_scf_env_type
77 : USE virial_methods, ONLY: virial_pair_force
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 :
91 : PUBLIC :: build_dftb_matrices, build_dftb_ks_matrix, build_dftb_overlap
92 :
93 : CONTAINS
94 :
95 : ! **************************************************************************************************
96 : !> \brief ...
97 : !> \param qs_env ...
98 : !> \param para_env ...
99 : !> \param calculate_forces ...
100 : ! **************************************************************************************************
101 4092 : 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 4092 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
113 : INTEGER, DIMENSION(3) :: cell
114 4092 : 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 8184 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dfblock, dsblock, fblock, fmatij, &
124 4092 : fmatji, pblock, sblock, scoeff, &
125 4092 : smatij, smatji, spxr, wblock
126 4092 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
127 : TYPE(atprop_type), POINTER :: atprop
128 16368 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
129 : TYPE(cp_logger_type), POINTER :: logger
130 4092 : 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
134 : TYPE(neighbor_list_iterator_p_type), &
135 4092 : DIMENSION(:), POINTER :: nl_iterator
136 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
137 4092 : POINTER :: sab_orb
138 4092 : 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 4092 : POINTER :: dftb_potential
142 : TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji
143 : TYPE(qs_energy_type), POINTER :: energy
144 4092 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
145 4092 : 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 4092 : CALL timeset(routineN, handle)
151 :
152 : ! set pointers
153 4092 : iptr = 0
154 20460 : DO la = 0, 3
155 85932 : DO lb = 0, 3
156 65472 : llm = 0
157 286440 : DO l1 = 0, MAX(la, lb)
158 597432 : DO l2 = 0, MIN(l1, la, lb)
159 1018908 : DO m = 0, l2
160 486948 : llm = llm + 1
161 814308 : iptr(l1, l2, m, la, lb) = llm
162 : END DO
163 : END DO
164 : END DO
165 : END DO
166 : END DO
167 :
168 4092 : NULLIFY (logger, virial, atprop)
169 4092 : logger => cp_get_default_logger()
170 :
171 4092 : NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
172 4092 : 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 4092 : ks_env=ks_env)
183 :
184 4092 : dftb_control => dft_control%qs_control%dftb_control
185 4092 : nimg = dft_control%nimages
186 : ! Allocate the overlap and Hamiltonian matrix
187 4092 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
188 4092 : nderivatives = 0
189 4092 : IF (dftb_control%self_consistent .AND. calculate_forces) nderivatives = 1
190 4092 : CALL setup_matrices2(qs_env, nderivatives, nimg, matrix_s, "OVERLAP", sab_orb)
191 4092 : CALL setup_matrices2(qs_env, 0, nimg, matrix_h, "CORE HAMILTONIAN", sab_orb)
192 4092 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
193 4092 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
194 :
195 4092 : NULLIFY (dftb_potential)
196 4092 : CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential)
197 4092 : NULLIFY (particle_set)
198 4092 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
199 :
200 4092 : IF (calculate_forces) THEN
201 786 : 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 786 : force=force)
207 786 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
208 :
209 786 : IF (SIZE(matrix_p, 1) == 2) THEN
210 340 : DO img = 1, nimg
211 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
212 170 : 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 340 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
215 : END DO
216 : END IF
217 786 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
218 786 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
219 : END IF
220 : ! atomic energy decomposition
221 4092 : IF (atprop%energy) THEN
222 162 : CALL atprop_array_init(atprop%atecc, natom=SIZE(particle_set))
223 : END IF
224 :
225 4092 : NULLIFY (cell_to_index)
226 4092 : IF (nimg > 1) THEN
227 960 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
228 960 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
229 : END IF
230 :
231 4092 : erep = 0._dp
232 :
233 4092 : nkind = SIZE(atomic_kind_set)
234 :
235 4092 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
236 773367 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
237 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
238 769275 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
239 769275 : 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 769275 : eta=eta_a, natorb=natorb_a)
243 769275 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
244 769275 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
245 : CALL get_dftb_atom_param(dftb_kind_b, &
246 769275 : defined=defined, lmax=lmaxj, eta=eta_b, natorb=natorb_b)
247 :
248 769275 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
249 :
250 : ! retrieve information on F and S matrix
251 769275 : dftb_param_ij => dftb_potential(ikind, jkind)
252 769275 : dftb_param_ji => dftb_potential(jkind, ikind)
253 : ! assume table size and type is symmetric
254 769275 : ngrd = dftb_param_ij%ngrd
255 769275 : ngrdcut = dftb_param_ij%ngrdcut
256 769275 : dgrd = dftb_param_ij%dgrd
257 769275 : ddr = dgrd*dftb_fd_deriv_step
258 769275 : CPASSERT(dftb_param_ij%llm == dftb_param_ji%llm)
259 769275 : llm = dftb_param_ij%llm
260 769275 : fmatij => dftb_param_ij%fmat
261 769275 : smatij => dftb_param_ij%smat
262 769275 : fmatji => dftb_param_ji%fmat
263 769275 : smatji => dftb_param_ji%smat
264 : ! repulsive pair potential
265 769275 : n_urpoly = dftb_param_ij%n_urpoly
266 769275 : urep_cut = dftb_param_ij%urep_cut
267 8462025 : urep = dftb_param_ij%urep
268 769275 : spxr => dftb_param_ij%spxr
269 769275 : scoeff => dftb_param_ij%scoeff
270 769275 : spdim = dftb_param_ij%spdim
271 769275 : s_cut = dftb_param_ij%s_cut
272 3077100 : srep = dftb_param_ij%srep
273 2307825 : surr = dftb_param_ij%surr
274 :
275 3077100 : dr = SQRT(SUM(rij(:)**2))
276 769275 : IF (NINT(dr/dgrd) <= ngrdcut) THEN
277 :
278 767613 : IF (nimg == 1) THEN
279 : ic = 1
280 : ELSE
281 74862 : ic = cell_to_index(cell(1), cell(2), cell(3))
282 74862 : CPASSERT(ic > 0)
283 : END IF
284 :
285 767613 : icol = MAX(iatom, jatom)
286 767613 : irow = MIN(iatom, jatom)
287 767613 : NULLIFY (sblock, fblock)
288 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
289 767613 : row=irow, col=icol, BLOCK=sblock, found=found)
290 767613 : CPASSERT(found)
291 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
292 767613 : row=irow, col=icol, BLOCK=fblock, found=found)
293 767613 : CPASSERT(found)
294 :
295 767613 : IF (calculate_forces) THEN
296 314670 : NULLIFY (pblock)
297 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
298 314670 : row=irow, col=icol, block=pblock, found=found)
299 314670 : CPASSERT(ASSOCIATED(pblock))
300 314670 : NULLIFY (wblock)
301 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
302 314670 : row=irow, col=icol, block=wblock, found=found)
303 314670 : CPASSERT(ASSOCIATED(wblock))
304 314670 : IF (dftb_control%self_consistent) THEN
305 1047392 : DO i = 2, 4
306 785544 : NULLIFY (dsblocks(i)%block)
307 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
308 785544 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
309 1047392 : CPASSERT(found)
310 : END DO
311 : END IF
312 : END IF
313 :
314 767613 : IF (iatom == jatom .AND. dr < 0.001_dp) THEN
315 : ! diagonal block
316 82055 : DO i = 1, natorb_a
317 57115 : sblock(i, i) = sblock(i, i) + 1._dp
318 82055 : 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 742673 : llm, lmaxi, lmaxj, irow, iatom)
324 : CALL compute_block_sk(fblock, fmatij, fmatji, rij, ngrd, ngrdcut, dgrd, &
325 742673 : llm, lmaxi, lmaxj, irow, iatom)
326 742673 : IF (calculate_forces) THEN
327 306072 : force_ab = 0._dp
328 306072 : force_w = 0._dp
329 306072 : n1 = SIZE(fblock, 1)
330 306072 : 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 306072 : f0 = 1.0_dp
334 306072 : IF (irow == iatom) f0 = -1.0_dp
335 :
336 1836432 : ALLOCATE (dfblock(n1, n2), dsblock(n1, n2))
337 :
338 1224288 : DO i = 1, 3
339 918216 : drij = rij
340 12505326 : dfblock = 0._dp; dsblock = 0._dp
341 :
342 918216 : drij(i) = rij(i) - ddr*f0
343 : CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
344 918216 : llm, lmaxi, lmaxj, irow, iatom)
345 : CALL compute_block_sk(dfblock, fmatij, fmatji, drij, ngrd, ngrdcut, dgrd, &
346 918216 : llm, lmaxi, lmaxj, irow, iatom)
347 :
348 6711771 : dsblock = -dsblock
349 6711771 : dfblock = -dfblock
350 :
351 918216 : drij(i) = rij(i) + ddr*f0
352 : CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
353 918216 : llm, lmaxi, lmaxj, irow, iatom)
354 : CALL compute_block_sk(dfblock, fmatij, fmatji, drij, ngrd, ngrdcut, dgrd, &
355 918216 : llm, lmaxi, lmaxj, irow, iatom)
356 :
357 6711771 : dfblock = dfblock/(2.0_dp*ddr)
358 6711771 : dsblock = dsblock/(2.0_dp*ddr)
359 :
360 6711771 : foab = 2.0_dp*SUM(dfblock*pblock)
361 6711771 : fow = -2.0_dp*SUM(dsblock*wblock)
362 :
363 918216 : force_ab(i) = force_ab(i) + foab
364 918216 : force_w(i) = force_w(i) + fow
365 1224288 : IF (dftb_control%self_consistent) THEN
366 763953 : CPASSERT(ASSOCIATED(dsblocks(i + 1)%block))
367 5559807 : dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock
368 : END IF
369 : END DO
370 306072 : IF (use_virial) THEN
371 190621 : IF (iatom == jatom) f0 = 0.5_dp*f0
372 190621 : CALL virial_pair_force(virial%pv_virial, -f0, force_ab, rij)
373 190621 : CALL virial_pair_force(virial%pv_virial, -f0, force_w, rij)
374 : END IF
375 306072 : DEALLOCATE (dfblock, dsblock)
376 : END IF
377 : END IF
378 :
379 767613 : IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
380 306072 : atom_a = atom_of_kind(iatom)
381 306072 : atom_b = atom_of_kind(jatom)
382 763029 : IF (irow == iatom) force_ab = -force_ab
383 763029 : IF (irow == iatom) force_w = -force_w
384 1224288 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
385 1224288 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
386 1224288 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - force_w(:)
387 1232886 : force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + force_w(:)
388 : END IF
389 :
390 : END IF
391 :
392 : ! repulsive potential
393 773367 : IF ((dr <= urep_cut .OR. spdim > 0) .AND. dr > 0.001_dp) THEN
394 581197 : erepij = 0._dp
395 : CALL urep_egr(rij, dr, erepij, force_rr, &
396 581197 : n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, calculate_forces)
397 581197 : erep = erep + erepij
398 581197 : IF (atprop%energy) THEN
399 228665 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*erepij
400 228665 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*erepij
401 : END IF
402 581197 : IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
403 257129 : atom_a = atom_of_kind(iatom)
404 257129 : atom_b = atom_of_kind(jatom)
405 : force(ikind)%repulsive(:, atom_a) = &
406 1028516 : force(ikind)%repulsive(:, atom_a) - force_rr(:)
407 : force(jkind)%repulsive(:, atom_b) = &
408 1028516 : force(jkind)%repulsive(:, atom_b) + force_rr(:)
409 257129 : IF (use_virial) THEN
410 171808 : f0 = -1.0_dp
411 171808 : IF (iatom == jatom) f0 = -0.5_dp
412 171808 : 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 4092 : CALL neighbor_list_iterator_release(nl_iterator)
419 :
420 10176 : DO i = 1, SIZE(matrix_s, 1)
421 36150 : DO img = 1, nimg
422 32058 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
423 : END DO
424 : END DO
425 8184 : DO i = 1, SIZE(matrix_h, 1)
426 28488 : DO img = 1, nimg
427 24396 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
428 : END DO
429 : END DO
430 :
431 : ! set repulsive energy
432 4092 : CALL para_env%sum(erep)
433 4092 : energy%repulsive = erep
434 :
435 4092 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
436 4092 : 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 0 : extension=".Log")
440 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
441 0 : after = MIN(MAX(after, 1), 16)
442 0 : DO img = 1, nimg
443 : CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, &
444 0 : output_unit=iw, omit_headers=omit_headers)
445 : END DO
446 :
447 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
448 0 : "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
449 : END IF
450 :
451 4092 : 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 0 : extension=".Log")
455 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
456 0 : after = MIN(MAX(after, 1), 16)
457 0 : DO img = 1, nimg
458 : CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, &
459 0 : output_unit=iw, omit_headers=omit_headers)
460 :
461 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
462 0 : qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
463 0 : 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 0 : 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 0 : "DFT%PRINT%AO_MATRICES/OVERLAP")
472 : END IF
473 :
474 4092 : IF (calculate_forces) THEN
475 786 : IF (SIZE(matrix_p, 1) == 2) THEN
476 340 : DO img = 1, nimg
477 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
478 170 : 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 340 : beta_scalar=-1.0_dp)
481 : END DO
482 : END IF
483 : END IF
484 :
485 4092 : CALL timestop(handle)
486 :
487 12276 : END SUBROUTINE build_dftb_matrices
488 :
489 : ! **************************************************************************************************
490 : !> \brief ...
491 : !> \param qs_env ...
492 : !> \param calculate_forces ...
493 : !> \param just_energy ...
494 : ! **************************************************************************************************
495 29784 : 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 29784 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mix_charge
506 29784 : REAL(KIND=dp), DIMENSION(:), POINTER :: mcharge, occupation_numbers
507 29784 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges
508 29784 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
509 : TYPE(cp_logger_type), POINTER :: logger
510 29784 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs
511 29784 : 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 29784 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
516 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
517 : TYPE(qs_energy_type), POINTER :: energy
518 29784 : 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 29784 : CALL timeset(routineN, handle)
525 29784 : NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, &
526 29784 : ks_matrix, rho, energy, scf_env)
527 29784 : logger => cp_get_default_logger()
528 29784 : 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 29784 : energy=energy)
540 :
541 29784 : energy%hartree = 0.0_dp
542 29784 : energy%qmmm_el = 0.0_dp
543 :
544 29784 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
545 29784 : nspins = dft_control%nspins
546 29784 : CPASSERT(ASSOCIATED(matrix_h))
547 29784 : CPASSERT(ASSOCIATED(rho))
548 89352 : CPASSERT(SIZE(ks_matrix) > 0)
549 :
550 61612 : DO ispin = 1, nspins
551 268994 : DO img = 1, SIZE(ks_matrix, 2)
552 : ! copy the core matrix into the fock matrix
553 239210 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
554 : END DO
555 : END DO
556 :
557 29784 : 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 28054 : matrix_s_kp=matrix_s)
561 28054 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
562 28054 : natom = SIZE(particle_set)
563 112216 : ALLOCATE (charges(natom, nspins))
564 : !
565 28054 : CALL mulliken_charges(matrix_p, matrix_s, para_env, charges)
566 : !
567 84162 : ALLOCATE (mcharge(natom))
568 28054 : nkind = SIZE(atomic_kind_set)
569 87464 : DO ikind = 1, nkind
570 59410 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
571 59410 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
572 59410 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
573 365596 : DO iatom = 1, natom
574 218722 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
575 504182 : mcharge(atom_a) = zeff - SUM(charges(atom_a, 1:nspins))
576 : END DO
577 : END DO
578 28054 : DEALLOCATE (charges)
579 :
580 28054 : 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 1944 : CALL get_qs_env(qs_env=qs_env, scf_env=scf_env)
583 5832 : ALLOCATE (mix_charge(SIZE(mcharge), 1))
584 16832 : 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 1944 : dft_control%qs_control%dftb_control%tblite_mixer_weight_factor)
599 16832 : mcharge(:) = mix_charge(:, 1)
600 1944 : DEALLOCATE (mix_charge)
601 : END IF
602 :
603 : CALL build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, &
604 28054 : calculate_forces, just_energy)
605 :
606 : CALL efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, &
607 28054 : calculate_forces, just_energy)
608 :
609 28054 : DEALLOCATE (mcharge)
610 :
611 : END IF
612 :
613 29784 : IF (qs_env%qmmm) THEN
614 4808 : CPASSERT(SIZE(ks_matrix, 2) == 1)
615 9616 : 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 4808 : 1.0_dp, 1.0_dp)
619 4808 : 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 4808 : matrix_p1(ispin)%matrix, qmmm_el)
623 9616 : energy%qmmm_el = energy%qmmm_el + qmmm_el
624 : END DO
625 4808 : pc_ener = qs_env%ks_qmmm_env%pc_ener
626 4808 : 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 29784 : energy%repulsive + energy%dispersion + energy%dftb3
631 :
632 29784 : 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 28054 : extension=".scfLog")
635 28054 : IF (output_unit > 0) THEN
636 : WRITE (UNIT=output_unit, FMT="(/,(T9,A,T60,F20.10))") &
637 15 : "Repulsive pair potential energy: ", energy%repulsive, &
638 15 : "Zeroth order Hamiltonian energy: ", energy%core, &
639 15 : "Charge fluctuation energy: ", energy%hartree, &
640 30 : "London dispersion energy: ", energy%dispersion
641 15 : IF (ABS(energy%efield) > 1.e-12_dp) THEN
642 : WRITE (UNIT=output_unit, FMT="(T9,A,T60,F20.10)") &
643 0 : "Electric field interaction energy: ", energy%efield
644 : END IF
645 15 : IF (dft_control%qs_control%dftb_control%dftb3_diagonal) THEN
646 : WRITE (UNIT=output_unit, FMT="(T9,A,T60,F20.10)") &
647 0 : "DFTB3 3rd Order Energy Correction ", energy%dftb3
648 : END IF
649 15 : IF (qs_env%qmmm) THEN
650 : WRITE (UNIT=output_unit, FMT="(T9,A,T60,F20.10)") &
651 0 : "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 28054 : "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 29784 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
659 386 : CPASSERT(SIZE(ks_matrix, 2) == 1)
660 : BLOCK
661 386 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
662 386 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
663 800 : DO ispin = 1, SIZE(mo_derivs)
664 : CALL get_mo_set(mo_set=mo_array(ispin), &
665 414 : mo_coeff_b=mo_coeff, occupation_numbers=occupation_numbers)
666 414 : IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
667 0 : CPABORT("")
668 : END IF
669 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
670 800 : 0.0_dp, mo_derivs(ispin)%matrix)
671 : END DO
672 : END BLOCK
673 : END IF
674 :
675 29784 : CALL timestop(handle)
676 :
677 59568 : END SUBROUTINE build_dftb_ks_matrix
678 :
679 : ! **************************************************************************************************
680 : !> \brief ...
681 : !> \param qs_env ...
682 : !> \param nderivative ...
683 : !> \param matrix_s ...
684 : ! **************************************************************************************************
685 1792 : 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 1792 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, dsblockm, sblock, smatij, smatji
700 896 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsblock1
701 896 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
702 8960 : 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 896 : DIMENSION(:), POINTER :: nl_iterator
708 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
709 896 : POINTER :: sab_orb
710 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a, dftb_kind_b
711 : TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
712 896 : POINTER :: dftb_potential
713 : TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji
714 896 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
715 :
716 896 : CALL timeset(routineN, handle)
717 :
718 : ! set pointers
719 896 : iptr = 0
720 4480 : DO la = 0, 3
721 18816 : DO lb = 0, 3
722 14336 : llm = 0
723 62720 : DO l1 = 0, MAX(la, lb)
724 130816 : DO l2 = 0, MIN(l1, la, lb)
725 223104 : DO m = 0, l2
726 106624 : llm = llm + 1
727 178304 : iptr(l1, l2, m, la, lb) = llm
728 : END DO
729 : END DO
730 : END DO
731 : END DO
732 : END DO
733 :
734 896 : NULLIFY (logger)
735 896 : logger => cp_get_default_logger()
736 :
737 896 : 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 896 : dft_control=dft_control)
742 :
743 896 : dftb_control => dft_control%qs_control%dftb_control
744 :
745 896 : NULLIFY (dftb_potential)
746 : CALL get_qs_env(qs_env=qs_env, &
747 896 : dftb_potential=dftb_potential)
748 :
749 896 : nkind = SIZE(atomic_kind_set)
750 :
751 : ! Allocate the overlap matrix
752 896 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
753 896 : CALL setup_matrices1(qs_env, nderivative, matrix_s, 'OVERLAP', sab_orb)
754 :
755 896 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
756 3668 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
757 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
758 2772 : iatom=iatom, jatom=jatom, r=rij)
759 :
760 2772 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
761 2772 : 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 2772 : natorb=natorb_a)
765 :
766 2772 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
767 :
768 2772 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
769 : CALL get_dftb_atom_param(dftb_kind_b, &
770 2772 : defined=defined, lmax=lmaxj, natorb=natorb_b)
771 :
772 2772 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
773 :
774 : ! retrieve information on F and S matrix
775 2772 : dftb_param_ij => dftb_potential(ikind, jkind)
776 2772 : dftb_param_ji => dftb_potential(jkind, ikind)
777 : ! assume table size and type is symmetric
778 2772 : ngrd = dftb_param_ij%ngrd
779 2772 : ngrdcut = dftb_param_ij%ngrdcut
780 2772 : dgrd = dftb_param_ij%dgrd
781 2772 : ddr = dgrd*dftb_fd_deriv_step
782 2772 : CPASSERT(dftb_param_ij%llm == dftb_param_ji%llm)
783 2772 : llm = dftb_param_ij%llm
784 2772 : smatij => dftb_param_ij%smat
785 2772 : smatji => dftb_param_ji%smat
786 :
787 11088 : dr = SQRT(SUM(rij(:)**2))
788 3668 : IF (NINT(dr/dgrd) <= ngrdcut) THEN
789 :
790 2772 : icol = MAX(iatom, jatom); irow = MIN(iatom, jatom)
791 :
792 2772 : NULLIFY (sblock)
793 : CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
794 2772 : row=irow, col=icol, BLOCK=sblock, found=found)
795 2772 : CPASSERT(found)
796 :
797 2772 : IF (nderivative > 0) THEN
798 3672 : DO i = 2, SIZE(matrix_s, 1)
799 3258 : NULLIFY (dsblocks(i)%block)
800 : CALL dbcsr_get_block_p(matrix=matrix_s(i)%matrix, &
801 3672 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
802 : END DO
803 : END IF
804 :
805 2772 : IF (iatom == jatom .AND. dr < 0.001_dp) THEN
806 : ! diagonal block
807 4137 : DO i = 1, natorb_a
808 4137 : 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 1407 : llm, lmaxi, lmaxj, irow, iatom)
814 :
815 1407 : IF (nderivative >= 1) THEN
816 228 : n1 = SIZE(sblock, 1); n2 = SIZE(sblock, 2)
817 228 : indder = 1 ! used to put the 2nd derivatives in the correct matric (5=xx,8=yy,10=zz)
818 :
819 2280 : ALLOCATE (dsblock1(n1, n2, 3), dsblock(n1, n2), dsblockm(n1, n2))
820 5394 : dsblock1 = 0.0_dp
821 912 : DO i = 1, 3
822 9648 : dsblock = 0._dp; dsblockm = 0.0_dp
823 684 : drij = rij
824 684 : f0 = 1.0_dp; IF (irow == iatom) f0 = -1.0_dp
825 :
826 684 : drij(i) = rij(i) - ddr*f0
827 : CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
828 684 : llm, lmaxi, lmaxj, irow, iatom)
829 :
830 684 : drij(i) = rij(i) + ddr*f0
831 : CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
832 684 : llm, lmaxi, lmaxj, irow, iatom)
833 :
834 9648 : dsblock1(:, :, i) = (dsblock + dsblockm)
835 9648 : dsblock = dsblock - dsblockm
836 5166 : dsblock = dsblock/(2.0_dp*ddr)
837 :
838 684 : CPASSERT(ASSOCIATED(dsblocks(i + 1)%block))
839 5166 : dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock
840 912 : IF (nderivative > 1) THEN
841 567 : indder = indder + 5 - i
842 4347 : dsblocks(indder)%block = 0.0_dp
843 : dsblocks(indder)%block = dsblocks(indder)%block + &
844 4347 : (dsblock1(:, :, i) - 2.0_dp*sblock)/ddr**2
845 : END IF
846 : END DO
847 :
848 228 : IF (nderivative > 1) THEN
849 567 : DO i = 1, 2
850 1134 : DO j = i + 1, 3
851 8127 : dsblock = 0._dp; dsblockm = 0.0_dp
852 567 : drij = rij
853 567 : f0 = 1.0_dp; IF (irow == iatom) f0 = -1.0_dp
854 :
855 567 : 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 567 : llm, lmaxi, lmaxj, irow, iatom)
858 :
859 567 : 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 567 : llm, lmaxi, lmaxj, irow, iatom)
862 :
863 567 : indder = 2 + 2*i + j
864 4347 : dsblocks(indder)%block = 0.0_dp
865 : dsblocks(indder)%block = &
866 : dsblocks(indder)%block + ( &
867 4725 : 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 228 : DEALLOCATE (dsblock1, dsblock, dsblockm)
873 : END IF
874 : END IF
875 : END IF
876 : END DO
877 896 : CALL neighbor_list_iterator_release(nl_iterator)
878 :
879 2626 : DO i = 1, SIZE(matrix_s, 1)
880 2626 : CALL dbcsr_finalize(matrix_s(i)%matrix)
881 : END DO
882 :
883 896 : CALL timestop(handle)
884 :
885 1792 : 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 896 : 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 896 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
910 896 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
911 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
912 896 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
913 896 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
914 :
915 896 : 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 896 : neighbor_list_id=neighbor_list_id)
923 :
924 896 : nkind = SIZE(atomic_kind_set)
925 896 : natom = SIZE(particle_set)
926 :
927 896 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
928 :
929 2688 : ALLOCATE (first_sgf(natom))
930 1792 : ALLOCATE (last_sgf(natom))
931 :
932 : CALL get_particle_set(particle_set, qs_kind_set, &
933 : first_sgf=first_sgf, &
934 896 : last_sgf=last_sgf)
935 :
936 896 : nmat = 0
937 896 : IF (nderivative == 0) nmat = 1
938 896 : IF (nderivative == 1) nmat = 4
939 896 : IF (nderivative == 2) nmat = 10
940 896 : CPASSERT(nmat > 0)
941 :
942 1792 : ALLOCATE (row_blk_sizes(natom))
943 896 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
944 :
945 896 : CALL dbcsr_allocate_matrix_set(matrices, nmat)
946 :
947 : ! Up to 2nd derivative take care to get the symmetries correct
948 2626 : DO i = 1, nmat
949 1730 : IF (i > 1) THEN
950 834 : matnames = TRIM(mnames)//" DERIVATIVE MATRIX DFTB"
951 834 : symmetry_type = dbcsr_type_antisymmetric
952 834 : IF (i > 4) symmetry_type = dbcsr_type_symmetric
953 : ELSE
954 896 : symmetry_type = dbcsr_type_symmetric
955 896 : matnames = TRIM(mnames)//" MATRIX DFTB"
956 : END IF
957 1730 : 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 1730 : mutable_work=.TRUE.)
963 2626 : CALL cp_dbcsr_alloc_block_from_nbl(matrices(i)%matrix, sab_nl)
964 : END DO
965 :
966 896 : DEALLOCATE (first_sgf)
967 896 : DEALLOCATE (last_sgf)
968 :
969 896 : DEALLOCATE (row_blk_sizes)
970 :
971 896 : 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 8184 : 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 8184 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
997 8184 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
998 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
999 8184 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1000 8184 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1001 :
1002 8184 : 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 8184 : neighbor_list_id=neighbor_list_id)
1010 :
1011 8184 : nkind = SIZE(atomic_kind_set)
1012 8184 : natom = SIZE(particle_set)
1013 :
1014 8184 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
1015 :
1016 24552 : ALLOCATE (first_sgf(natom))
1017 16368 : ALLOCATE (last_sgf(natom))
1018 :
1019 : CALL get_particle_set(particle_set, qs_kind_set, &
1020 : first_sgf=first_sgf, &
1021 8184 : last_sgf=last_sgf)
1022 :
1023 8184 : nmat = 0
1024 8184 : IF (nderivative == 0) nmat = 1
1025 8184 : IF (nderivative == 1) nmat = 4
1026 8184 : IF (nderivative == 2) nmat = 10
1027 8184 : CPASSERT(nmat > 0)
1028 :
1029 16368 : ALLOCATE (row_blk_sizes(natom))
1030 8184 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
1031 :
1032 8184 : CALL dbcsr_allocate_matrix_set(matrices, nmat, nimg)
1033 :
1034 : ! Up to 2nd derivative take care to get the symmetries correct
1035 48792 : DO img = 1, nimg
1036 95070 : DO i = 1, nmat
1037 46278 : IF (i > 1) THEN
1038 5670 : matnames = TRIM(mnames)//" DERIVATIVE MATRIX DFTB"
1039 5670 : symmetry_type = dbcsr_type_antisymmetric
1040 5670 : IF (i > 4) symmetry_type = dbcsr_type_symmetric
1041 : ELSE
1042 40608 : symmetry_type = dbcsr_type_symmetric
1043 40608 : matnames = TRIM(mnames)//" MATRIX DFTB"
1044 : END IF
1045 46278 : 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 46278 : mutable_work=.TRUE.)
1051 86886 : CALL cp_dbcsr_alloc_block_from_nbl(matrices(i, img)%matrix, sab_nl)
1052 : END DO
1053 : END DO
1054 :
1055 8184 : DEALLOCATE (first_sgf)
1056 8184 : DEALLOCATE (last_sgf)
1057 :
1058 8184 : DEALLOCATE (row_blk_sizes)
1059 :
1060 8184 : END SUBROUTINE setup_matrices2
1061 :
1062 : END MODULE qs_dftb_matrices
|
