Line data Source code
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_api, ONLY: &
22 : dbcsr_add, dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, &
23 : dbcsr_distribution_type, dbcsr_dot, dbcsr_finalize, dbcsr_get_block_p, dbcsr_multiply, &
24 : dbcsr_p_type, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_symmetric
25 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
26 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
27 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_type
30 : USE cp_output_handling, ONLY: cp_p_file,&
31 : cp_print_key_finished_output,&
32 : cp_print_key_should_output,&
33 : cp_print_key_unit_nr
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 2822 : 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 2822 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
108 : INTEGER, DIMENSION(3) :: cell
109 2822 : 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 5644 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dfblock, dsblock, fblock, fmatij, &
119 2822 : fmatji, pblock, sblock, scoeff, &
120 2822 : smatij, smatji, spxr, wblock
121 2822 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
122 : TYPE(atprop_type), POINTER :: atprop
123 11288 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
124 : TYPE(cp_logger_type), POINTER :: logger
125 2822 : 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 2822 : DIMENSION(:), POINTER :: nl_iterator
131 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
132 2822 : POINTER :: sab_orb
133 2822 : 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 2822 : POINTER :: dftb_potential
137 : TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji
138 : TYPE(qs_energy_type), POINTER :: energy
139 2822 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
140 2822 : 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 2822 : CALL timeset(routineN, handle)
146 :
147 : ! set pointers
148 2822 : iptr = 0
149 14110 : DO la = 0, 3
150 59262 : DO lb = 0, 3
151 45152 : llm = 0
152 197540 : DO l1 = 0, MAX(la, lb)
153 412012 : DO l2 = 0, MIN(l1, la, lb)
154 702678 : DO m = 0, l2
155 335818 : llm = llm + 1
156 561578 : iptr(l1, l2, m, la, lb) = llm
157 : END DO
158 : END DO
159 : END DO
160 : END DO
161 : END DO
162 :
163 2822 : NULLIFY (logger, virial, atprop)
164 2822 : logger => cp_get_default_logger()
165 :
166 2822 : NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
167 2822 : 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 2822 : ks_env=ks_env)
178 :
179 2822 : dftb_control => dft_control%qs_control%dftb_control
180 2822 : nimg = dft_control%nimages
181 : ! Allocate the overlap and Hamiltonian matrix
182 2822 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
183 2822 : nderivatives = 0
184 2822 : IF (dftb_control%self_consistent .AND. calculate_forces) nderivatives = 1
185 2822 : CALL setup_matrices2(qs_env, nderivatives, nimg, matrix_s, "OVERLAP", sab_orb)
186 2822 : CALL setup_matrices2(qs_env, 0, nimg, matrix_h, "CORE HAMILTONIAN", sab_orb)
187 2822 : CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
188 2822 : CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
189 :
190 2822 : NULLIFY (dftb_potential)
191 2822 : CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential)
192 2822 : NULLIFY (particle_set)
193 2822 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
194 :
195 2822 : IF (calculate_forces) THEN
196 724 : 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 724 : force=force)
202 724 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
203 :
204 724 : IF (SIZE(matrix_p, 1) == 2) THEN
205 328 : DO img = 1, nimg
206 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
207 164 : 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 328 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
210 : END DO
211 : END IF
212 724 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
213 724 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
214 : END IF
215 : ! atomic energy decomposition
216 2822 : IF (atprop%energy) THEN
217 162 : CALL atprop_array_init(atprop%atecc, natom=SIZE(particle_set))
218 : END IF
219 :
220 2822 : NULLIFY (cell_to_index)
221 2822 : IF (nimg > 1) THEN
222 290 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
223 290 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
224 : END IF
225 :
226 2822 : erep = 0._dp
227 :
228 2822 : nkind = SIZE(atomic_kind_set)
229 :
230 2822 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
231 726673 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
232 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
233 723851 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
234 723851 : 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 723851 : eta=eta_a, natorb=natorb_a)
238 723851 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
239 723851 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
240 : CALL get_dftb_atom_param(dftb_kind_b, &
241 723851 : defined=defined, lmax=lmaxj, eta=eta_b, natorb=natorb_b)
242 :
243 723851 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
244 :
245 : ! retrieve information on F and S matrix
246 723851 : dftb_param_ij => dftb_potential(ikind, jkind)
247 723851 : dftb_param_ji => dftb_potential(jkind, ikind)
248 : ! assume table size and type is symmetric
249 723851 : ngrd = dftb_param_ij%ngrd
250 723851 : ngrdcut = dftb_param_ij%ngrdcut
251 723851 : dgrd = dftb_param_ij%dgrd
252 723851 : ddr = dgrd*0.1_dp
253 723851 : CPASSERT(dftb_param_ij%llm == dftb_param_ji%llm)
254 723851 : llm = dftb_param_ij%llm
255 723851 : fmatij => dftb_param_ij%fmat
256 723851 : smatij => dftb_param_ij%smat
257 723851 : fmatji => dftb_param_ji%fmat
258 723851 : smatji => dftb_param_ji%smat
259 : ! repulsive pair potential
260 723851 : n_urpoly = dftb_param_ij%n_urpoly
261 723851 : urep_cut = dftb_param_ij%urep_cut
262 7962361 : urep = dftb_param_ij%urep
263 723851 : spxr => dftb_param_ij%spxr
264 723851 : scoeff => dftb_param_ij%scoeff
265 723851 : spdim = dftb_param_ij%spdim
266 723851 : s_cut = dftb_param_ij%s_cut
267 2895404 : srep = dftb_param_ij%srep
268 2171553 : surr = dftb_param_ij%surr
269 :
270 2895404 : dr = SQRT(SUM(rij(:)**2))
271 723851 : IF (NINT(dr/dgrd) <= ngrdcut) THEN
272 :
273 722189 : IF (nimg == 1) THEN
274 : ic = 1
275 : ELSE
276 34363 : ic = cell_to_index(cell(1), cell(2), cell(3))
277 34363 : CPASSERT(ic > 0)
278 : END IF
279 :
280 722189 : icol = MAX(iatom, jatom)
281 722189 : irow = MIN(iatom, jatom)
282 722189 : NULLIFY (sblock, fblock)
283 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
284 722189 : row=irow, col=icol, BLOCK=sblock, found=found)
285 722189 : CPASSERT(found)
286 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
287 722189 : row=irow, col=icol, BLOCK=fblock, found=found)
288 722189 : CPASSERT(found)
289 :
290 722189 : IF (calculate_forces) THEN
291 312241 : NULLIFY (pblock)
292 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
293 312241 : row=irow, col=icol, block=pblock, found=found)
294 312241 : CPASSERT(ASSOCIATED(pblock))
295 312241 : NULLIFY (wblock)
296 : CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
297 312241 : row=irow, col=icol, block=wblock, found=found)
298 312241 : CPASSERT(ASSOCIATED(wblock))
299 312241 : IF (dftb_control%self_consistent) THEN
300 1038628 : DO i = 2, 4
301 778971 : NULLIFY (dsblocks(i)%block)
302 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
303 778971 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
304 1038628 : CPASSERT(found)
305 : END DO
306 : END IF
307 : END IF
308 :
309 722189 : IF (iatom == jatom .AND. dr < 0.001_dp) THEN
310 : ! diagonal block
311 71100 : DO i = 1, natorb_a
312 48649 : sblock(i, i) = sblock(i, i) + 1._dp
313 71100 : 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 699738 : llm, lmaxi, lmaxj, irow, iatom)
319 : CALL compute_block_sk(fblock, fmatij, fmatji, rij, ngrd, ngrdcut, dgrd, &
320 699738 : llm, lmaxi, lmaxj, irow, iatom)
321 699738 : IF (calculate_forces) THEN
322 303765 : force_ab = 0._dp
323 303765 : force_w = 0._dp
324 303765 : n1 = SIZE(fblock, 1)
325 303765 : 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 303765 : f0 = 1.0_dp
329 303765 : IF (irow == iatom) f0 = -1.0_dp
330 :
331 1822590 : ALLOCATE (dfblock(n1, n2), dsblock(n1, n2))
332 :
333 1215060 : DO i = 1, 3
334 911295 : drij = rij
335 12094617 : dfblock = 0._dp; dsblock = 0._dp
336 :
337 911295 : drij(i) = rij(i) - ddr*f0
338 : CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
339 911295 : llm, lmaxi, lmaxj, irow, iatom)
340 : CALL compute_block_sk(dfblock, fmatij, fmatji, drij, ngrd, ngrdcut, dgrd, &
341 911295 : llm, lmaxi, lmaxj, irow, iatom)
342 :
343 6502956 : dsblock = -dsblock
344 6502956 : dfblock = -dfblock
345 :
346 911295 : drij(i) = rij(i) + ddr*f0
347 : CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
348 911295 : llm, lmaxi, lmaxj, irow, iatom)
349 : CALL compute_block_sk(dfblock, fmatij, fmatji, drij, ngrd, ngrdcut, dgrd, &
350 911295 : llm, lmaxi, lmaxj, irow, iatom)
351 :
352 6502956 : dfblock = dfblock/(2.0_dp*ddr)
353 6502956 : dsblock = dsblock/(2.0_dp*ddr)
354 :
355 6502956 : foab = 2.0_dp*SUM(dfblock*pblock)
356 6502956 : fow = -2.0_dp*SUM(dsblock*wblock)
357 :
358 911295 : force_ab(i) = force_ab(i) + foab
359 911295 : force_w(i) = force_w(i) + fow
360 1215060 : IF (dftb_control%self_consistent) THEN
361 757650 : CPASSERT(ASSOCIATED(dsblocks(i + 1)%block))
362 5395806 : dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock
363 : END IF
364 : END DO
365 303765 : IF (use_virial) THEN
366 188575 : IF (iatom == jatom) f0 = 0.5_dp*f0
367 188575 : CALL virial_pair_force(virial%pv_virial, -f0, force_ab, rij)
368 188575 : CALL virial_pair_force(virial%pv_virial, -f0, force_w, rij)
369 188575 : IF (atprop%stress) THEN
370 82273 : f1 = 0.5_dp*f0
371 82273 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -f1, force_ab, rij)
372 82273 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -f1, force_w, rij)
373 82273 : CALL virial_pair_force(atprop%atstress(:, :, jatom), -f1, force_ab, rij)
374 82273 : CALL virial_pair_force(atprop%atstress(:, :, jatom), -f1, force_w, rij)
375 : END IF
376 : END IF
377 303765 : DEALLOCATE (dfblock, dsblock)
378 : END IF
379 : END IF
380 :
381 722189 : IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
382 303765 : atom_a = atom_of_kind(iatom)
383 303765 : atom_b = atom_of_kind(jatom)
384 756729 : IF (irow == iatom) force_ab = -force_ab
385 756729 : IF (irow == iatom) force_w = -force_w
386 1215060 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
387 1215060 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
388 1215060 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - force_w(:)
389 1223536 : force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + force_w(:)
390 : END IF
391 :
392 : END IF
393 :
394 : ! repulsive potential
395 726673 : IF ((dr <= urep_cut .OR. spdim > 0) .AND. dr > 0.001_dp) THEN
396 542295 : erepij = 0._dp
397 : CALL urep_egr(rij, dr, erepij, force_rr, &
398 542295 : n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, calculate_forces)
399 542295 : erep = erep + erepij
400 542295 : IF (atprop%energy) THEN
401 228665 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*erepij
402 228665 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*erepij
403 : END IF
404 542295 : IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
405 254994 : atom_a = atom_of_kind(iatom)
406 254994 : atom_b = atom_of_kind(jatom)
407 : force(ikind)%repulsive(:, atom_a) = &
408 1019976 : force(ikind)%repulsive(:, atom_a) - force_rr(:)
409 : force(jkind)%repulsive(:, atom_b) = &
410 1019976 : force(jkind)%repulsive(:, atom_b) + force_rr(:)
411 254994 : IF (use_virial) THEN
412 169917 : f0 = -1.0_dp
413 169917 : IF (iatom == jatom) f0 = -0.5_dp
414 169917 : CALL virial_pair_force(virial%pv_virial, f0, force_rr, rij)
415 169917 : IF (atprop%stress) THEN
416 71092 : CALL virial_pair_force(atprop%atstress(:, :, iatom), f0*0.5_dp, force_rr, rij)
417 71092 : 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 2822 : CALL neighbor_list_iterator_release(nl_iterator)
425 :
426 7492 : DO i = 1, SIZE(matrix_s, 1)
427 20984 : DO img = 1, nimg
428 18162 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
429 : END DO
430 : END DO
431 5644 : DO i = 1, SIZE(matrix_h, 1)
432 15026 : DO img = 1, nimg
433 12204 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
434 : END DO
435 : END DO
436 :
437 : ! set repulsive energy
438 2822 : CALL para_env%sum(erep)
439 2822 : energy%repulsive = erep
440 :
441 2822 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
442 2822 : 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 0 : extension=".Log")
446 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
447 0 : after = MIN(MAX(after, 1), 16)
448 0 : DO img = 1, nimg
449 : CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, &
450 0 : output_unit=iw, omit_headers=omit_headers)
451 : END DO
452 :
453 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
454 0 : "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
455 : END IF
456 :
457 2822 : 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 0 : extension=".Log")
461 0 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
462 0 : after = MIN(MAX(after, 1), 16)
463 0 : DO img = 1, nimg
464 : CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, &
465 0 : output_unit=iw, omit_headers=omit_headers)
466 :
467 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
468 0 : qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
469 0 : 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 0 : 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 0 : "DFT%PRINT%AO_MATRICES/OVERLAP")
478 : END IF
479 :
480 2822 : IF (calculate_forces) THEN
481 724 : IF (SIZE(matrix_p, 1) == 2) THEN
482 328 : DO img = 1, nimg
483 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
484 328 : beta_scalar=-1.0_dp)
485 : END DO
486 : END IF
487 : END IF
488 :
489 2822 : CALL timestop(handle)
490 :
491 8466 : END SUBROUTINE build_dftb_matrices
492 :
493 : ! **************************************************************************************************
494 : !> \brief ...
495 : !> \param qs_env ...
496 : !> \param calculate_forces ...
497 : !> \param just_energy ...
498 : ! **************************************************************************************************
499 13050 : 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 13050 : REAL(KIND=dp), DIMENSION(:), POINTER :: mcharge, occupation_numbers
510 13050 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges
511 13050 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
512 : TYPE(cp_logger_type), POINTER :: logger
513 13050 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs
514 13050 : 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 13050 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
519 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
520 : TYPE(qs_energy_type), POINTER :: energy
521 13050 : 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 13050 : CALL timeset(routineN, handle)
527 13050 : NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, &
528 13050 : ks_matrix, rho, energy)
529 13050 : logger => cp_get_default_logger()
530 13050 : 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 13050 : energy=energy)
542 :
543 13050 : energy%hartree = 0.0_dp
544 13050 : energy%qmmm_el = 0.0_dp
545 :
546 13050 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
547 13050 : nspins = dft_control%nspins
548 13050 : CPASSERT(ASSOCIATED(matrix_h))
549 13050 : CPASSERT(ASSOCIATED(rho))
550 39150 : CPASSERT(SIZE(ks_matrix) > 0)
551 :
552 26380 : DO ispin = 1, nspins
553 85856 : DO img = 1, SIZE(ks_matrix, 2)
554 : ! copy the core matrix into the fock matrix
555 72806 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
556 : END DO
557 : END DO
558 :
559 13050 : 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 11918 : matrix_s_kp=matrix_s)
563 11918 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
564 11918 : natom = SIZE(particle_set)
565 47672 : ALLOCATE (charges(natom, nspins))
566 : !
567 11918 : CALL mulliken_charges(matrix_p, matrix_s, para_env, charges)
568 : !
569 35754 : ALLOCATE (mcharge(natom))
570 11918 : nkind = SIZE(atomic_kind_set)
571 37402 : DO ikind = 1, nkind
572 25484 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
573 25484 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
574 25484 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
575 218168 : DO iatom = 1, natom
576 155282 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
577 337168 : mcharge(atom_a) = zeff - SUM(charges(atom_a, 1:nspins))
578 : END DO
579 : END DO
580 11918 : DEALLOCATE (charges)
581 :
582 : CALL build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, &
583 11918 : calculate_forces, just_energy)
584 :
585 : CALL efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, &
586 11918 : calculate_forces, just_energy)
587 :
588 11918 : DEALLOCATE (mcharge)
589 :
590 : END IF
591 :
592 13050 : IF (qs_env%qmmm) THEN
593 4808 : CPASSERT(SIZE(ks_matrix, 2) == 1)
594 9616 : 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 4808 : 1.0_dp, 1.0_dp)
598 4808 : 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 4808 : matrix_p1(ispin)%matrix, qmmm_el)
602 9616 : energy%qmmm_el = energy%qmmm_el + qmmm_el
603 : END DO
604 4808 : pc_ener = qs_env%ks_qmmm_env%pc_ener
605 4808 : 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 13050 : energy%repulsive + energy%dispersion + energy%dftb3
610 :
611 13050 : IF (dft_control%qs_control%dftb_control%self_consistent) THEN
612 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
613 11918 : extension=".scfLog")
614 11918 : IF (output_unit > 0) THEN
615 : WRITE (UNIT=output_unit, FMT="(/,(T9,A,T60,F20.10))") &
616 15 : "Repulsive pair potential energy: ", energy%repulsive, &
617 15 : "Zeroth order Hamiltonian energy: ", energy%core, &
618 15 : "Charge fluctuation energy: ", energy%hartree, &
619 30 : "London dispersion energy: ", energy%dispersion
620 15 : IF (ABS(energy%efield) > 1.e-12_dp) THEN
621 : WRITE (UNIT=output_unit, FMT="(T9,A,T60,F20.10)") &
622 0 : "Electric field interaction energy: ", energy%efield
623 : END IF
624 15 : IF (dft_control%qs_control%dftb_control%dftb3_diagonal) THEN
625 : WRITE (UNIT=output_unit, FMT="(T9,A,T60,F20.10)") &
626 0 : "DFTB3 3rd Order Energy Correction ", energy%dftb3
627 : END IF
628 15 : IF (qs_env%qmmm) THEN
629 : WRITE (UNIT=output_unit, FMT="(T9,A,T60,F20.10)") &
630 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
631 : END IF
632 : END IF
633 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
634 11918 : "PRINT%DETAILED_ENERGY")
635 : END IF
636 : ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C (plus occupation numbers)
637 13050 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
638 336 : CPASSERT(SIZE(ks_matrix, 2) == 1)
639 : BLOCK
640 336 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
641 336 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
642 700 : DO ispin = 1, SIZE(mo_derivs)
643 : CALL get_mo_set(mo_set=mo_array(ispin), &
644 364 : mo_coeff_b=mo_coeff, occupation_numbers=occupation_numbers)
645 364 : IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
646 0 : CPABORT("")
647 : END IF
648 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
649 700 : 0.0_dp, mo_derivs(ispin)%matrix)
650 : END DO
651 : END BLOCK
652 : END IF
653 :
654 13050 : CALL timestop(handle)
655 :
656 13050 : END SUBROUTINE build_dftb_ks_matrix
657 :
658 : ! **************************************************************************************************
659 : !> \brief ...
660 : !> \param qs_env ...
661 : !> \param nderivative ...
662 : !> \param matrix_s ...
663 : ! **************************************************************************************************
664 1792 : SUBROUTINE build_dftb_overlap(qs_env, nderivative, matrix_s)
665 :
666 : TYPE(qs_environment_type), POINTER :: qs_env
667 : INTEGER, INTENT(IN) :: nderivative
668 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
669 :
670 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_dftb_overlap'
671 :
672 : INTEGER :: handle, i, iatom, icol, ikind, indder, irow, j, jatom, jkind, l1, l2, la, lb, &
673 : llm, lmaxi, lmaxj, m, n1, n2, natom, natorb_a, natorb_b, ngrd, ngrdcut, nkind
674 : LOGICAL :: defined, found
675 : REAL(KIND=dp) :: ddr, dgrd, dr, f0
676 : REAL(KIND=dp), DIMENSION(0:3) :: skself
677 : REAL(KIND=dp), DIMENSION(3) :: drij, rij
678 1792 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, dsblockm, sblock, smatij, smatji
679 896 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsblock1
680 896 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
681 8960 : TYPE(block_p_type), DIMENSION(2:10) :: dsblocks
682 : TYPE(cp_logger_type), POINTER :: logger
683 : TYPE(dft_control_type), POINTER :: dft_control
684 : TYPE(dftb_control_type), POINTER :: dftb_control
685 : TYPE(neighbor_list_iterator_p_type), &
686 896 : DIMENSION(:), POINTER :: nl_iterator
687 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
688 896 : POINTER :: sab_orb
689 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a, dftb_kind_b
690 : TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
691 896 : POINTER :: dftb_potential
692 : TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji
693 896 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
694 :
695 896 : CALL timeset(routineN, handle)
696 :
697 : ! set pointers
698 896 : iptr = 0
699 4480 : DO la = 0, 3
700 18816 : DO lb = 0, 3
701 14336 : llm = 0
702 62720 : DO l1 = 0, MAX(la, lb)
703 130816 : DO l2 = 0, MIN(l1, la, lb)
704 223104 : DO m = 0, l2
705 106624 : llm = llm + 1
706 178304 : iptr(l1, l2, m, la, lb) = llm
707 : END DO
708 : END DO
709 : END DO
710 : END DO
711 : END DO
712 :
713 896 : NULLIFY (logger)
714 896 : logger => cp_get_default_logger()
715 :
716 896 : NULLIFY (atomic_kind_set, qs_kind_set, sab_orb)
717 :
718 : CALL get_qs_env(qs_env=qs_env, &
719 : atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
720 896 : dft_control=dft_control)
721 :
722 896 : dftb_control => dft_control%qs_control%dftb_control
723 :
724 896 : NULLIFY (dftb_potential)
725 : CALL get_qs_env(qs_env=qs_env, &
726 896 : dftb_potential=dftb_potential)
727 :
728 896 : nkind = SIZE(atomic_kind_set)
729 :
730 : ! Allocate the overlap matrix
731 896 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
732 896 : CALL setup_matrices1(qs_env, nderivative, matrix_s, 'OVERLAP', sab_orb)
733 :
734 896 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
735 3668 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
736 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
737 2772 : iatom=iatom, jatom=jatom, r=rij)
738 :
739 2772 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
740 2772 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
741 : CALL get_dftb_atom_param(dftb_kind_a, &
742 : defined=defined, lmax=lmaxi, skself=skself, &
743 2772 : natorb=natorb_a)
744 :
745 2772 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
746 :
747 2772 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
748 : CALL get_dftb_atom_param(dftb_kind_b, &
749 2772 : defined=defined, lmax=lmaxj, natorb=natorb_b)
750 :
751 2772 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
752 :
753 : ! retrieve information on F and S matrix
754 2772 : dftb_param_ij => dftb_potential(ikind, jkind)
755 2772 : dftb_param_ji => dftb_potential(jkind, ikind)
756 : ! assume table size and type is symmetric
757 2772 : ngrd = dftb_param_ij%ngrd
758 2772 : ngrdcut = dftb_param_ij%ngrdcut
759 2772 : dgrd = dftb_param_ij%dgrd
760 2772 : ddr = dgrd*0.1_dp
761 2772 : CPASSERT(dftb_param_ij%llm == dftb_param_ji%llm)
762 2772 : llm = dftb_param_ij%llm
763 2772 : smatij => dftb_param_ij%smat
764 2772 : smatji => dftb_param_ji%smat
765 :
766 11088 : dr = SQRT(SUM(rij(:)**2))
767 3668 : IF (NINT(dr/dgrd) <= ngrdcut) THEN
768 :
769 2772 : icol = MAX(iatom, jatom); irow = MIN(iatom, jatom)
770 :
771 2772 : NULLIFY (sblock)
772 : CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
773 2772 : row=irow, col=icol, BLOCK=sblock, found=found)
774 2772 : CPASSERT(found)
775 :
776 2772 : IF (nderivative .GT. 0) THEN
777 3672 : DO i = 2, SIZE(matrix_s, 1)
778 3258 : NULLIFY (dsblocks(i)%block)
779 : CALL dbcsr_get_block_p(matrix=matrix_s(i)%matrix, &
780 3672 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
781 : END DO
782 : END IF
783 :
784 2772 : IF (iatom == jatom .AND. dr < 0.001_dp) THEN
785 : ! diagonal block
786 4137 : DO i = 1, natorb_a
787 4137 : sblock(i, i) = sblock(i, i) + 1._dp
788 : END DO
789 : ELSE
790 : ! off-diagonal block
791 : CALL compute_block_sk(sblock, smatij, smatji, rij, ngrd, ngrdcut, dgrd, &
792 1407 : llm, lmaxi, lmaxj, irow, iatom)
793 :
794 1407 : IF (nderivative .GE. 1) THEN
795 228 : n1 = SIZE(sblock, 1); n2 = SIZE(sblock, 2)
796 228 : indder = 1 ! used to put the 2nd derivatives in the correct matric (5=xx,8=yy,10=zz)
797 :
798 2280 : ALLOCATE (dsblock1(n1, n2, 3), dsblock(n1, n2), dsblockm(n1, n2))
799 5394 : dsblock1 = 0.0_dp
800 912 : DO i = 1, 3
801 9648 : dsblock = 0._dp; dsblockm = 0.0_dp
802 684 : drij = rij
803 684 : f0 = 1.0_dp; IF (irow == iatom) f0 = -1.0_dp
804 :
805 684 : drij(i) = rij(i) - ddr*f0
806 : CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
807 684 : llm, lmaxi, lmaxj, irow, iatom)
808 :
809 684 : drij(i) = rij(i) + ddr*f0
810 : CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
811 684 : llm, lmaxi, lmaxj, irow, iatom)
812 :
813 9648 : dsblock1(:, :, i) = (dsblock + dsblockm)
814 9648 : dsblock = dsblock - dsblockm
815 5166 : dsblock = dsblock/(2.0_dp*ddr)
816 :
817 684 : CPASSERT(ASSOCIATED(dsblocks(i + 1)%block))
818 5166 : dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock
819 912 : IF (nderivative .GT. 1) THEN
820 567 : indder = indder + 5 - i
821 4347 : dsblocks(indder)%block = 0.0_dp
822 : dsblocks(indder)%block = dsblocks(indder)%block + &
823 4347 : (dsblock1(:, :, i) - 2.0_dp*sblock)/ddr**2
824 : END IF
825 : END DO
826 :
827 228 : IF (nderivative .GT. 1) THEN
828 567 : DO i = 1, 2
829 1134 : DO j = i + 1, 3
830 8127 : dsblock = 0._dp; dsblockm = 0.0_dp
831 567 : drij = rij
832 567 : f0 = 1.0_dp; IF (irow == iatom) f0 = -1.0_dp
833 :
834 567 : drij(i) = rij(i) - ddr*f0; drij(j) = rij(j) - ddr*f0
835 : CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
836 567 : llm, lmaxi, lmaxj, irow, iatom)
837 :
838 567 : drij(i) = rij(i) + ddr*f0; drij(j) = rij(j) + ddr*f0
839 : CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
840 567 : llm, lmaxi, lmaxj, irow, iatom)
841 :
842 567 : indder = 2 + 2*i + j
843 4347 : dsblocks(indder)%block = 0.0_dp
844 : dsblocks(indder)%block = &
845 : dsblocks(indder)%block + ( &
846 4725 : dsblock + dsblockm - dsblock1(:, :, i) - dsblock1(:, :, j) + 2.0_dp*sblock)/(2.0_dp*ddr**2)
847 : END DO
848 : END DO
849 : END IF
850 :
851 228 : DEALLOCATE (dsblock1, dsblock, dsblockm)
852 : END IF
853 : END IF
854 : END IF
855 : END DO
856 896 : CALL neighbor_list_iterator_release(nl_iterator)
857 :
858 2626 : DO i = 1, SIZE(matrix_s, 1)
859 2626 : CALL dbcsr_finalize(matrix_s(i)%matrix)
860 : END DO
861 :
862 896 : CALL timestop(handle)
863 :
864 1792 : END SUBROUTINE build_dftb_overlap
865 :
866 : ! **************************************************************************************************
867 : !> \brief ...
868 : !> \param qs_env ...
869 : !> \param nderivative ...
870 : !> \param matrices ...
871 : !> \param mnames ...
872 : !> \param sab_nl ...
873 : ! **************************************************************************************************
874 896 : SUBROUTINE setup_matrices1(qs_env, nderivative, matrices, mnames, sab_nl)
875 :
876 : TYPE(qs_environment_type), POINTER :: qs_env
877 : INTEGER, INTENT(IN) :: nderivative
878 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
879 : CHARACTER(LEN=*) :: mnames
880 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
881 : POINTER :: sab_nl
882 :
883 : CHARACTER(1) :: symmetry_type
884 : CHARACTER(LEN=default_string_length) :: matnames
885 : INTEGER :: i, natom, neighbor_list_id, nkind, nmat, &
886 : nsgf
887 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
888 896 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
889 896 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
890 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
891 896 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
892 896 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
893 :
894 896 : NULLIFY (particle_set, atomic_kind_set)
895 :
896 : CALL get_qs_env(qs_env=qs_env, &
897 : atomic_kind_set=atomic_kind_set, &
898 : qs_kind_set=qs_kind_set, &
899 : particle_set=particle_set, &
900 : dbcsr_dist=dbcsr_dist, &
901 896 : neighbor_list_id=neighbor_list_id)
902 :
903 896 : nkind = SIZE(atomic_kind_set)
904 896 : natom = SIZE(particle_set)
905 :
906 896 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
907 :
908 2688 : ALLOCATE (first_sgf(natom))
909 1792 : ALLOCATE (last_sgf(natom))
910 :
911 : CALL get_particle_set(particle_set, qs_kind_set, &
912 : first_sgf=first_sgf, &
913 896 : last_sgf=last_sgf)
914 :
915 896 : nmat = 0
916 896 : IF (nderivative == 0) nmat = 1
917 896 : IF (nderivative == 1) nmat = 4
918 896 : IF (nderivative == 2) nmat = 10
919 896 : CPASSERT(nmat > 0)
920 :
921 1792 : ALLOCATE (row_blk_sizes(natom))
922 896 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
923 :
924 896 : CALL dbcsr_allocate_matrix_set(matrices, nmat)
925 :
926 : ! Up to 2nd derivative take care to get the symmetries correct
927 2626 : DO i = 1, nmat
928 1730 : IF (i .GT. 1) THEN
929 834 : matnames = TRIM(mnames)//" DERIVATIVE MATRIX DFTB"
930 834 : symmetry_type = dbcsr_type_antisymmetric
931 834 : IF (i .GT. 4) symmetry_type = dbcsr_type_symmetric
932 : ELSE
933 896 : symmetry_type = dbcsr_type_symmetric
934 896 : matnames = TRIM(mnames)//" MATRIX DFTB"
935 : END IF
936 1730 : ALLOCATE (matrices(i)%matrix)
937 : CALL dbcsr_create(matrix=matrices(i)%matrix, &
938 : name=TRIM(matnames), &
939 : dist=dbcsr_dist, matrix_type=symmetry_type, &
940 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
941 1730 : nze=0, mutable_work=.TRUE.)
942 2626 : CALL cp_dbcsr_alloc_block_from_nbl(matrices(i)%matrix, sab_nl)
943 : END DO
944 :
945 896 : DEALLOCATE (first_sgf)
946 896 : DEALLOCATE (last_sgf)
947 :
948 896 : DEALLOCATE (row_blk_sizes)
949 :
950 896 : END SUBROUTINE setup_matrices1
951 :
952 : ! **************************************************************************************************
953 : !> \brief ...
954 : !> \param qs_env ...
955 : !> \param nderivative ...
956 : !> \param nimg ...
957 : !> \param matrices ...
958 : !> \param mnames ...
959 : !> \param sab_nl ...
960 : ! **************************************************************************************************
961 5644 : SUBROUTINE setup_matrices2(qs_env, nderivative, nimg, matrices, mnames, sab_nl)
962 :
963 : TYPE(qs_environment_type), POINTER :: qs_env
964 : INTEGER, INTENT(IN) :: nderivative, nimg
965 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrices
966 : CHARACTER(LEN=*) :: mnames
967 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
968 : POINTER :: sab_nl
969 :
970 : CHARACTER(1) :: symmetry_type
971 : CHARACTER(LEN=default_string_length) :: matnames
972 : INTEGER :: i, img, natom, neighbor_list_id, nkind, &
973 : nmat, nsgf
974 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
975 5644 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
976 5644 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
977 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
978 5644 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
979 5644 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
980 :
981 5644 : NULLIFY (particle_set, atomic_kind_set)
982 :
983 : CALL get_qs_env(qs_env=qs_env, &
984 : atomic_kind_set=atomic_kind_set, &
985 : qs_kind_set=qs_kind_set, &
986 : particle_set=particle_set, &
987 : dbcsr_dist=dbcsr_dist, &
988 5644 : neighbor_list_id=neighbor_list_id)
989 :
990 5644 : nkind = SIZE(atomic_kind_set)
991 5644 : natom = SIZE(particle_set)
992 :
993 5644 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
994 :
995 16932 : ALLOCATE (first_sgf(natom))
996 11288 : ALLOCATE (last_sgf(natom))
997 :
998 : CALL get_particle_set(particle_set, qs_kind_set, &
999 : first_sgf=first_sgf, &
1000 5644 : last_sgf=last_sgf)
1001 :
1002 5644 : nmat = 0
1003 5644 : IF (nderivative == 0) nmat = 1
1004 5644 : IF (nderivative == 1) nmat = 4
1005 5644 : IF (nderivative == 2) nmat = 10
1006 5644 : CPASSERT(nmat > 0)
1007 :
1008 11288 : ALLOCATE (row_blk_sizes(natom))
1009 5644 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
1010 :
1011 5644 : CALL dbcsr_allocate_matrix_set(matrices, nmat, nimg)
1012 :
1013 : ! Up to 2nd derivative take care to get the symmetries correct
1014 24408 : DO img = 1, nimg
1015 47282 : DO i = 1, nmat
1016 22874 : IF (i .GT. 1) THEN
1017 4110 : matnames = TRIM(mnames)//" DERIVATIVE MATRIX DFTB"
1018 4110 : symmetry_type = dbcsr_type_antisymmetric
1019 4110 : IF (i .GT. 4) symmetry_type = dbcsr_type_symmetric
1020 : ELSE
1021 18764 : symmetry_type = dbcsr_type_symmetric
1022 18764 : matnames = TRIM(mnames)//" MATRIX DFTB"
1023 : END IF
1024 22874 : ALLOCATE (matrices(i, img)%matrix)
1025 : CALL dbcsr_create(matrix=matrices(i, img)%matrix, &
1026 : name=TRIM(matnames), &
1027 : dist=dbcsr_dist, matrix_type=symmetry_type, &
1028 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
1029 22874 : nze=0, mutable_work=.TRUE.)
1030 41638 : CALL cp_dbcsr_alloc_block_from_nbl(matrices(i, img)%matrix, sab_nl)
1031 : END DO
1032 : END DO
1033 :
1034 5644 : DEALLOCATE (first_sgf)
1035 5644 : DEALLOCATE (last_sgf)
1036 :
1037 5644 : DEALLOCATE (row_blk_sizes)
1038 :
1039 5644 : END SUBROUTINE setup_matrices2
1040 :
1041 : END MODULE qs_dftb_matrices
1042 :
|