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