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 cp_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 23990 : 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 23990 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
90 : INTEGER, DIMENSION(3) :: cellind
91 23990 : 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 23990 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, ksblock, pblock, sblock
96 23990 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
97 23990 : 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 23990 : 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 23990 : DIMENSION(:), POINTER :: nl_iterator
107 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
108 23990 : POINTER :: n_list
109 23990 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
110 23990 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
111 23990 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
112 : TYPE(virial_type), POINTER :: virial
113 :
114 23990 : CALL timeset(routineN, handle)
115 23990 : NULLIFY (atprop)
116 :
117 : ! Energy
118 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
119 23990 : qs_kind_set=qs_kind_set, atprop=atprop)
120 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
121 23990 : kind_of=kind_of, atom_of_kind=atom_of_kind)
122 :
123 23990 : eb3 = 0.0_dp
124 23990 : CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
125 85444 : DO ikind = 1, SIZE(local_particles%n_el)
126 61454 : ua = xgamma(ikind)
127 213434 : DO ia = 1, local_particles%n_el(ikind)
128 127990 : iatom = local_particles%list(ikind)%array(ia)
129 127990 : eloc = -1.0_dp/6.0_dp*ua*mcharge(iatom)**3
130 127990 : eb3 = eb3 + eloc
131 189444 : 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 23990 : CALL get_qs_env(qs_env=qs_env, para_env=para_env)
139 23990 : CALL para_env%sum(eb3)
140 23990 : energy%dftb3 = eb3
141 :
142 : ! Forces and Virial
143 23990 : 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 : END IF
230 : END DO
231 : END DO
232 : END DO
233 : END DO
234 : END IF
235 : ELSE
236 58 : NULLIFY (n_list)
237 58 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
238 58 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
239 33640 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
240 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
241 33582 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
242 :
243 134328 : dr = SQRT(SUM(rij**2))
244 33582 : IF (iatom == jatom .AND. dr < 1.0e-6_dp) CYCLE
245 :
246 33410 : icol = MAX(iatom, jatom)
247 33410 : irow = MIN(iatom, jatom)
248 :
249 33410 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
250 33410 : CPASSERT(ic > 0)
251 :
252 33410 : ikind = kind_of(iatom)
253 33410 : atom_i = atom_of_kind(iatom)
254 33410 : ui = xgamma(ikind)
255 33410 : jkind = kind_of(jatom)
256 33410 : atom_j = atom_of_kind(jatom)
257 33410 : uj = xgamma(jkind)
258 : !
259 33410 : gmij = -0.5_dp*(ui*mcharge(iatom)**2 + uj*mcharge(jatom)**2)
260 : !
261 33410 : NULLIFY (pblock)
262 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
263 33410 : row=irow, col=icol, block=pblock, found=found)
264 33410 : CPASSERT(found)
265 133640 : DO i = 1, 3
266 100230 : NULLIFY (dsblock)
267 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
268 100230 : row=irow, col=icol, block=dsblock, found=found)
269 100230 : CPASSERT(found)
270 100230 : IF (irow == iatom) THEN
271 3677568 : fi = -gmij*SUM(pblock*dsblock)
272 : ELSE
273 2379069 : fi = gmij*SUM(pblock*dsblock)
274 : END IF
275 100230 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
276 100230 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
277 233870 : fij(i) = fi
278 : END DO
279 33468 : IF (use_virial) THEN
280 22502 : fi = 1.0_dp
281 22502 : IF (iatom == jatom) fi = 0.5_dp
282 22502 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
283 : END IF
284 :
285 : END DO
286 58 : CALL neighbor_list_iterator_release(nl_iterator)
287 : !
288 : END IF
289 :
290 472 : IF (SIZE(matrix_p, 1) == 2) THEN
291 432 : DO ic = 1, SIZE(matrix_p, 2)
292 : CALL dbcsr_add(matrix_p(1, ic)%matrix, matrix_p(2, ic)%matrix, &
293 432 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
294 : END DO
295 : END IF
296 : END IF
297 :
298 : ! KS matrix
299 23990 : IF (.NOT. just_energy) THEN
300 23888 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s, natom=natom)
301 : !
302 23888 : nimg = SIZE(ks_matrix, 2)
303 23888 : NULLIFY (cell_to_index)
304 23888 : IF (nimg > 1) THEN
305 4614 : NULLIFY (kpoints)
306 4614 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
307 4614 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
308 : END IF
309 :
310 23888 : IF (nimg == 1) THEN
311 : ! no k-points; all matrices have been transformed to periodic bsf
312 19274 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
313 1868987 : DO WHILE (dbcsr_iterator_blocks_left(iter))
314 1849713 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
315 1849713 : ikind = kind_of(irow)
316 1849713 : ui = xgamma(ikind)
317 1849713 : jkind = kind_of(icol)
318 1849713 : uj = xgamma(jkind)
319 1849713 : gmij = -0.5_dp*(ui*mcharge(irow)**2 + uj*mcharge(icol)**2)
320 3727076 : DO is = 1, SIZE(ks_matrix, 1)
321 1858089 : NULLIFY (ksblock)
322 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
323 1858089 : row=irow, col=icol, block=ksblock, found=found)
324 1858089 : CPASSERT(found)
325 45806073 : ksblock = ksblock - 0.5_dp*gmij*sblock
326 : END DO
327 : END DO
328 19274 : CALL dbcsr_iterator_stop(iter)
329 : ELSE
330 4614 : NULLIFY (n_list)
331 4614 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
332 4614 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
333 2199332 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
334 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
335 2194718 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
336 :
337 2194718 : icol = MAX(iatom, jatom)
338 2194718 : irow = MIN(iatom, jatom)
339 :
340 2194718 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
341 2194718 : CPASSERT(ic > 0)
342 :
343 2194718 : ikind = kind_of(iatom)
344 2194718 : ui = xgamma(ikind)
345 2194718 : jkind = kind_of(jatom)
346 2194718 : uj = xgamma(jkind)
347 2194718 : gmij = -0.5_dp*(ui*mcharge(iatom)**2 + uj*mcharge(jatom)**2)
348 : !
349 2194718 : NULLIFY (sblock)
350 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
351 2194718 : row=irow, col=icol, block=sblock, found=found)
352 2194718 : CPASSERT(found)
353 4578300 : DO is = 1, SIZE(ks_matrix, 1)
354 2378968 : NULLIFY (ksblock)
355 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
356 2378968 : row=irow, col=icol, block=ksblock, found=found)
357 2378968 : CPASSERT(found)
358 155645036 : ksblock = ksblock - 0.5_dp*gmij*sblock
359 : END DO
360 :
361 : END DO
362 4614 : CALL neighbor_list_iterator_release(nl_iterator)
363 : !
364 : END IF
365 : END IF
366 :
367 23990 : CALL timestop(handle)
368 :
369 47980 : END SUBROUTINE build_dftb3_diagonal
370 :
371 : END MODULE qs_dftb3_methods
372 :
|
