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 DFTB3 Terms
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE qs_dftb3_methods
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind_set
15 : USE atprop_types, ONLY: atprop_type
16 : USE cell_types, ONLY: cell_type
17 : USE dbcsr_api, ONLY: dbcsr_add,&
18 : dbcsr_get_block_p,&
19 : dbcsr_iterator_blocks_left,&
20 : dbcsr_iterator_next_block,&
21 : dbcsr_iterator_start,&
22 : dbcsr_iterator_stop,&
23 : dbcsr_iterator_type,&
24 : dbcsr_p_type
25 : USE distribution_1d_types, ONLY: distribution_1d_type
26 : USE kinds, ONLY: dp
27 : USE kpoint_types, ONLY: get_kpoint_info,&
28 : kpoint_type
29 : USE message_passing, ONLY: mp_para_env_type
30 : USE particle_types, ONLY: particle_type
31 : USE qs_energy_types, ONLY: qs_energy_type
32 : USE qs_environment_types, ONLY: get_qs_env,&
33 : qs_environment_type
34 : USE qs_force_types, ONLY: qs_force_type
35 : USE qs_kind_types, ONLY: qs_kind_type
36 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
37 : neighbor_list_iterate,&
38 : neighbor_list_iterator_create,&
39 : neighbor_list_iterator_p_type,&
40 : neighbor_list_iterator_release,&
41 : neighbor_list_set_p_type
42 : USE qs_rho_types, ONLY: qs_rho_get,&
43 : qs_rho_type
44 : USE sap_kind_types, ONLY: sap_int_type
45 : USE virial_methods, ONLY: virial_pair_force
46 : USE virial_types, ONLY: virial_type
47 : #include "./base/base_uses.f90"
48 :
49 : IMPLICIT NONE
50 :
51 : PRIVATE
52 :
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb3_methods'
54 :
55 : PUBLIC :: build_dftb3_diagonal
56 :
57 : CONTAINS
58 :
59 : ! **************************************************************************************************
60 : !> \brief ...
61 : !> \param qs_env ...
62 : !> \param ks_matrix ...
63 : !> \param rho ...
64 : !> \param mcharge ...
65 : !> \param energy ...
66 : !> \param xgamma ...
67 : !> \param zeff ...
68 : !> \param sap_int ...
69 : !> \param calculate_forces ...
70 : !> \param just_energy ...
71 : ! **************************************************************************************************
72 23956 : SUBROUTINE build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeff, &
73 : sap_int, calculate_forces, just_energy)
74 :
75 : TYPE(qs_environment_type), POINTER :: qs_env
76 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
77 : TYPE(qs_rho_type), POINTER :: rho
78 : REAL(dp), DIMENSION(:), INTENT(in) :: mcharge
79 : TYPE(qs_energy_type), POINTER :: energy
80 : REAL(dp), DIMENSION(:), INTENT(in) :: xgamma, zeff
81 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
82 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
83 :
84 : CHARACTER(len=*), PARAMETER :: routineN = 'build_dftb3_diagonal'
85 :
86 : INTEGER :: atom_i, atom_j, blk, handle, i, ia, iac, &
87 : iatom, ic, icol, ikind, irow, is, &
88 : jatom, jkind, natom, nimg, nkind
89 23956 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
90 : INTEGER, DIMENSION(3) :: cellind
91 23956 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
92 : LOGICAL :: found, use_virial
93 : REAL(KIND=dp) :: dr, eb3, eloc, fi, gmij, ua, ui, uj
94 : REAL(KIND=dp), DIMENSION(3) :: fij, rij
95 23956 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, ksblock, pblock, sblock
96 23956 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
97 23956 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
98 : TYPE(atprop_type), POINTER :: atprop
99 : TYPE(cell_type), POINTER :: cell
100 : TYPE(dbcsr_iterator_type) :: iter
101 23956 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
102 : TYPE(distribution_1d_type), POINTER :: local_particles
103 : TYPE(kpoint_type), POINTER :: kpoints
104 : TYPE(mp_para_env_type), POINTER :: para_env
105 : TYPE(neighbor_list_iterator_p_type), &
106 23956 : DIMENSION(:), POINTER :: nl_iterator
107 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
108 23956 : POINTER :: n_list
109 23956 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
110 23956 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
111 23956 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
112 : TYPE(virial_type), POINTER :: virial
113 :
114 23956 : CALL timeset(routineN, handle)
115 23956 : NULLIFY (atprop)
116 :
117 : ! Energy
118 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
119 23956 : qs_kind_set=qs_kind_set, atprop=atprop)
120 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
121 23956 : kind_of=kind_of, atom_of_kind=atom_of_kind)
122 :
123 23956 : eb3 = 0.0_dp
124 23956 : CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
125 85314 : DO ikind = 1, SIZE(local_particles%n_el)
126 61358 : ua = xgamma(ikind)
127 213150 : DO ia = 1, local_particles%n_el(ikind)
128 127836 : iatom = local_particles%list(ikind)%array(ia)
129 127836 : eloc = -1.0_dp/6.0_dp*ua*mcharge(iatom)**3
130 127836 : eb3 = eb3 + eloc
131 189194 : IF (atprop%energy) THEN
132 : ! we have to add the part not covered by 0.5*Tr(FP)
133 9756 : eloc = -0.5_dp*eloc - 0.25_dp*ua*zeff(ikind)*mcharge(iatom)**2
134 9756 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + eloc
135 : END IF
136 : END DO
137 : END DO
138 23956 : CALL get_qs_env(qs_env=qs_env, para_env=para_env)
139 23956 : CALL para_env%sum(eb3)
140 23956 : energy%dftb3 = eb3
141 :
142 : ! Forces and Virial
143 23956 : IF (calculate_forces) THEN
144 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s, natom=natom, force=force, &
145 472 : cell=cell, virial=virial, particle_set=particle_set)
146 : ! virial
147 472 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
148 472 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
149 472 : IF (SIZE(matrix_p, 1) == 2) THEN
150 432 : DO ic = 1, SIZE(matrix_p, 2)
151 : CALL dbcsr_add(matrix_p(1, ic)%matrix, matrix_p(2, ic)%matrix, &
152 432 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
153 : END DO
154 : END IF
155 : !
156 472 : nimg = SIZE(matrix_p, 2)
157 472 : NULLIFY (cell_to_index)
158 472 : IF (nimg > 1) THEN
159 58 : NULLIFY (kpoints)
160 58 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
161 58 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
162 : END IF
163 472 : IF (nimg == 1) THEN
164 : ! no k-points; all matrices have been transformed to periodic bsf
165 414 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
166 103689 : DO WHILE (dbcsr_iterator_blocks_left(iter))
167 103275 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
168 103275 : ikind = kind_of(irow)
169 103275 : atom_i = atom_of_kind(irow)
170 103275 : ui = xgamma(ikind)
171 103275 : jkind = kind_of(icol)
172 103275 : atom_j = atom_of_kind(icol)
173 103275 : uj = xgamma(jkind)
174 : !
175 103275 : gmij = -0.5_dp*(ui*mcharge(irow)**2 + uj*mcharge(icol)**2)
176 : !
177 103275 : NULLIFY (pblock)
178 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
179 103275 : row=irow, col=icol, block=pblock, found=found)
180 103275 : CPASSERT(found)
181 413514 : DO i = 1, 3
182 309825 : NULLIFY (dsblock)
183 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
184 309825 : row=irow, col=icol, block=dsblock, found=found)
185 309825 : CPASSERT(found)
186 2755113 : fi = -gmij*SUM(pblock*dsblock)
187 309825 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
188 309825 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
189 722925 : fij(i) = fi
190 : END DO
191 : END DO
192 414 : CALL dbcsr_iterator_stop(iter)
193 : ! use dsint list
194 414 : IF (use_virial) THEN
195 70 : CPASSERT(ASSOCIATED(sap_int))
196 70 : CALL get_qs_env(qs_env, nkind=nkind)
197 226 : DO ikind = 1, nkind
198 610 : DO jkind = 1, nkind
199 384 : iac = ikind + nkind*(jkind - 1)
200 384 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
201 274 : ui = xgamma(ikind)
202 274 : uj = xgamma(jkind)
203 4488 : DO ia = 1, sap_int(iac)%nalist
204 4058 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE
205 4038 : iatom = sap_int(iac)%alist(ia)%aatom
206 138257 : DO ic = 1, sap_int(iac)%alist(ia)%nclist
207 133835 : jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
208 535340 : rij = sap_int(iac)%alist(ia)%clist(ic)%rac
209 535340 : dr = SQRT(SUM(rij(:)**2))
210 137893 : IF (dr > 1.e-6_dp) THEN
211 131821 : dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
212 131821 : gmij = -0.5_dp*(ui*mcharge(iatom)**2 + uj*mcharge(jatom)**2)
213 131821 : icol = MAX(iatom, jatom)
214 131821 : irow = MIN(iatom, jatom)
215 131821 : NULLIFY (pblock)
216 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
217 131821 : row=irow, col=icol, block=pblock, found=found)
218 131821 : CPASSERT(found)
219 527284 : DO i = 1, 3
220 527284 : IF (irow == iatom) THEN
221 2614437 : fij(i) = -gmij*SUM(pblock*dsint(:, :, i))
222 : ELSE
223 2236134 : fij(i) = -gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i))
224 : END IF
225 : END DO
226 131821 : fi = 1.0_dp
227 131821 : IF (iatom == jatom) fi = 0.5_dp
228 131821 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
229 131821 : IF (atprop%stress) THEN
230 98393 : CALL virial_pair_force(atprop%atstress(:, :, irow), fi*0.5_dp, fij, rij)
231 98393 : CALL virial_pair_force(atprop%atstress(:, :, icol), fi*0.5_dp, fij, rij)
232 : END IF
233 : END IF
234 : END DO
235 : END DO
236 : END DO
237 : END DO
238 : END IF
239 : ELSE
240 58 : NULLIFY (n_list)
241 58 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
242 58 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
243 33640 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
244 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
245 33582 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
246 :
247 134328 : dr = SQRT(SUM(rij**2))
248 33582 : IF (iatom == jatom .AND. dr < 1.0e-6_dp) CYCLE
249 :
250 33410 : icol = MAX(iatom, jatom)
251 33410 : irow = MIN(iatom, jatom)
252 :
253 33410 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
254 33410 : CPASSERT(ic > 0)
255 :
256 33410 : ikind = kind_of(iatom)
257 33410 : atom_i = atom_of_kind(iatom)
258 33410 : ui = xgamma(ikind)
259 33410 : jkind = kind_of(jatom)
260 33410 : atom_j = atom_of_kind(jatom)
261 33410 : uj = xgamma(jkind)
262 : !
263 33410 : gmij = -0.5_dp*(ui*mcharge(iatom)**2 + uj*mcharge(jatom)**2)
264 : !
265 33410 : NULLIFY (pblock)
266 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
267 33410 : row=irow, col=icol, block=pblock, found=found)
268 33410 : CPASSERT(found)
269 133640 : DO i = 1, 3
270 100230 : NULLIFY (dsblock)
271 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
272 100230 : row=irow, col=icol, block=dsblock, found=found)
273 100230 : CPASSERT(found)
274 100230 : IF (irow == iatom) THEN
275 3677568 : fi = -gmij*SUM(pblock*dsblock)
276 : ELSE
277 2379069 : fi = gmij*SUM(pblock*dsblock)
278 : END IF
279 100230 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
280 100230 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
281 233870 : fij(i) = fi
282 : END DO
283 33468 : IF (use_virial) THEN
284 22502 : fi = 1.0_dp
285 22502 : IF (iatom == jatom) fi = 0.5_dp
286 22502 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
287 22502 : IF (atprop%stress) THEN
288 0 : CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij)
289 0 : CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij)
290 : END IF
291 : END IF
292 :
293 : END DO
294 58 : CALL neighbor_list_iterator_release(nl_iterator)
295 : !
296 : END IF
297 :
298 472 : IF (SIZE(matrix_p, 1) == 2) THEN
299 432 : DO ic = 1, SIZE(matrix_p, 2)
300 : CALL dbcsr_add(matrix_p(1, ic)%matrix, matrix_p(2, ic)%matrix, &
301 432 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
302 : END DO
303 : END IF
304 : END IF
305 :
306 : ! KS matrix
307 23956 : IF (.NOT. just_energy) THEN
308 23854 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s, natom=natom)
309 : !
310 23854 : nimg = SIZE(ks_matrix, 2)
311 23854 : NULLIFY (cell_to_index)
312 23854 : IF (nimg > 1) THEN
313 4614 : NULLIFY (kpoints)
314 4614 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
315 4614 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
316 : END IF
317 :
318 23854 : IF (nimg == 1) THEN
319 : ! no k-points; all matrices have been transformed to periodic bsf
320 19240 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
321 1868170 : DO WHILE (dbcsr_iterator_blocks_left(iter))
322 1848930 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
323 1848930 : ikind = kind_of(irow)
324 1848930 : ui = xgamma(ikind)
325 1848930 : jkind = kind_of(icol)
326 1848930 : uj = xgamma(jkind)
327 1848930 : gmij = -0.5_dp*(ui*mcharge(irow)**2 + uj*mcharge(icol)**2)
328 3725476 : DO is = 1, SIZE(ks_matrix, 1)
329 1857306 : NULLIFY (ksblock)
330 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
331 1857306 : row=irow, col=icol, block=ksblock, found=found)
332 1857306 : CPASSERT(found)
333 45789208 : ksblock = ksblock - 0.5_dp*gmij*sblock
334 : END DO
335 : END DO
336 19240 : CALL dbcsr_iterator_stop(iter)
337 : ELSE
338 4614 : NULLIFY (n_list)
339 4614 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
340 4614 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
341 2199332 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
342 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
343 2194718 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
344 :
345 2194718 : icol = MAX(iatom, jatom)
346 2194718 : irow = MIN(iatom, jatom)
347 :
348 2194718 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
349 2194718 : CPASSERT(ic > 0)
350 :
351 2194718 : ikind = kind_of(iatom)
352 2194718 : ui = xgamma(ikind)
353 2194718 : jkind = kind_of(jatom)
354 2194718 : uj = xgamma(jkind)
355 2194718 : gmij = -0.5_dp*(ui*mcharge(iatom)**2 + uj*mcharge(jatom)**2)
356 : !
357 2194718 : NULLIFY (sblock)
358 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
359 2194718 : row=irow, col=icol, block=sblock, found=found)
360 2194718 : CPASSERT(found)
361 4578300 : DO is = 1, SIZE(ks_matrix, 1)
362 2378968 : NULLIFY (ksblock)
363 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
364 2378968 : row=irow, col=icol, block=ksblock, found=found)
365 2378968 : CPASSERT(found)
366 155645036 : ksblock = ksblock - 0.5_dp*gmij*sblock
367 : END DO
368 :
369 : END DO
370 4614 : CALL neighbor_list_iterator_release(nl_iterator)
371 : !
372 : END IF
373 : END IF
374 :
375 23956 : CALL timestop(handle)
376 :
377 47912 : END SUBROUTINE build_dftb3_diagonal
378 :
379 : END MODULE qs_dftb3_methods
380 :
|
