Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of Coulomb Hessian contributions in xTB
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE xtb_ehess
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind_set
15 : USE cell_types, ONLY: cell_type,&
16 : get_cell,&
17 : pbc
18 : USE cp_control_types, ONLY: dft_control_type,&
19 : xtb_control_type
20 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
21 : dbcsr_iterator_blocks_left,&
22 : dbcsr_iterator_next_block,&
23 : dbcsr_iterator_start,&
24 : dbcsr_iterator_stop,&
25 : dbcsr_iterator_type,&
26 : dbcsr_p_type
27 : USE distribution_1d_types, ONLY: distribution_1d_type
28 : USE ewald_environment_types, ONLY: ewald_env_get,&
29 : ewald_environment_type
30 : USE ewald_methods_tb, ONLY: tb_ewald_overlap,&
31 : tb_spme_evaluate
32 : USE ewald_pw_types, ONLY: ewald_pw_type
33 : USE kinds, ONLY: dp
34 : USE mathconstants, ONLY: oorootpi,&
35 : pi
36 : USE message_passing, ONLY: mp_para_env_type
37 : USE particle_types, ONLY: particle_type
38 : USE pw_poisson_types, ONLY: do_ewald_ewald,&
39 : do_ewald_none,&
40 : do_ewald_pme,&
41 : do_ewald_spme
42 : USE qs_environment_types, ONLY: get_qs_env,&
43 : qs_environment_type
44 : USE qs_kind_types, ONLY: get_qs_kind,&
45 : qs_kind_type
46 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
47 : neighbor_list_iterate,&
48 : neighbor_list_iterator_create,&
49 : neighbor_list_iterator_p_type,&
50 : neighbor_list_iterator_release,&
51 : neighbor_list_set_p_type
52 : USE virial_types, ONLY: virial_type
53 : USE xtb_coulomb, ONLY: gamma_rab_sr
54 : USE xtb_spinpol, ONLY: xtb_spinpol_hessian
55 : USE xtb_types, ONLY: get_xtb_atom_param,&
56 : xtb_atom_type
57 : #include "./base/base_uses.f90"
58 :
59 : IMPLICIT NONE
60 :
61 : PRIVATE
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_ehess'
64 :
65 : PUBLIC :: xtb_coulomb_hessian
66 :
67 : CONTAINS
68 :
69 : ! **************************************************************************************************
70 : !> \brief ...
71 : !> \param qs_env ...
72 : !> \param ks_matrix ...
73 : !> \param charges1 ...
74 : !> \param mcharge1 ...
75 : !> \param mcharge ...
76 : !> \param matrix_p1 ...
77 : ! **************************************************************************************************
78 250 : SUBROUTINE xtb_coulomb_hessian(qs_env, ks_matrix, charges1, mcharge1, mcharge, matrix_p1)
79 :
80 : TYPE(qs_environment_type), POINTER :: qs_env
81 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix
82 : REAL(dp), DIMENSION(:, :) :: charges1
83 : REAL(dp), DIMENSION(:) :: mcharge1, mcharge
84 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1
85 :
86 : CHARACTER(len=*), PARAMETER :: routineN = 'xtb_coulomb_hessian'
87 :
88 : INTEGER :: ewald_type, handle, i, ia, iatom, icol, ikind, irow, is, j, jatom, jkind, la, lb, &
89 : lmaxa, lmaxb, natom, natorb_a, natorb_b, ni, nj, nkind, nmat, za, zb
90 250 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
91 : INTEGER, DIMENSION(25) :: laoa, laob
92 : INTEGER, DIMENSION(3) :: cellind, periodic
93 : LOGICAL :: defined, do_ewald, found
94 : REAL(KIND=dp) :: alpha, deth, dr, etaa, etab, gmij, kg, &
95 : rcut, rcuta, rcutb
96 250 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xgamma
97 250 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gammab, gcij, gmcharge
98 250 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gchrg
99 : REAL(KIND=dp), DIMENSION(3) :: rij
100 : REAL(KIND=dp), DIMENSION(5) :: kappaa, kappab
101 250 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksblock, sblock
102 250 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
103 : TYPE(cell_type), POINTER :: cell
104 : TYPE(dbcsr_iterator_type) :: iter
105 250 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
106 : TYPE(dft_control_type), POINTER :: dft_control
107 : TYPE(distribution_1d_type), POINTER :: local_particles
108 : TYPE(ewald_environment_type), POINTER :: ewald_env
109 : TYPE(ewald_pw_type), POINTER :: ewald_pw
110 : TYPE(mp_para_env_type), POINTER :: para_env
111 : TYPE(neighbor_list_iterator_p_type), &
112 250 : DIMENSION(:), POINTER :: nl_iterator
113 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
114 250 : POINTER :: n_list
115 250 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
116 250 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
117 : TYPE(virial_type), POINTER :: virial
118 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind
119 : TYPE(xtb_control_type), POINTER :: xtb_control
120 :
121 250 : CALL timeset(routineN, handle)
122 :
123 : CALL get_qs_env(qs_env, &
124 : matrix_s_kp=matrix_s, &
125 : qs_kind_set=qs_kind_set, &
126 : particle_set=particle_set, &
127 : cell=cell, &
128 250 : dft_control=dft_control)
129 :
130 250 : xtb_control => dft_control%qs_control%xtb_control
131 :
132 250 : IF (dft_control%nimages /= 1) THEN
133 0 : CPABORT("No kpoints allowed in xTB response calculation")
134 : END IF
135 :
136 250 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
137 250 : nmat = 1
138 750 : ALLOCATE (gchrg(natom, 5, nmat))
139 250 : gchrg = 0._dp
140 750 : ALLOCATE (gmcharge(natom, nmat))
141 250 : gmcharge = 0._dp
142 :
143 : ! short range contribution (gamma)
144 : ! loop over all atom pairs (sab_xtbe)
145 250 : kg = xtb_control%kg
146 250 : NULLIFY (n_list)
147 250 : CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
148 250 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
149 201138 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
150 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
151 200888 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
152 200888 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
153 200888 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
154 200888 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
155 200888 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
156 200888 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
157 200888 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
158 : ! atomic parameters
159 200888 : CALL get_xtb_atom_param(xtb_atom_a, eta=etaa, lmax=lmaxa, kappa=kappaa, rcut=rcuta)
160 200888 : CALL get_xtb_atom_param(xtb_atom_b, eta=etab, lmax=lmaxb, kappa=kappab, rcut=rcutb)
161 : ! gamma matrix
162 200888 : ni = lmaxa + 1
163 200888 : nj = lmaxb + 1
164 803552 : ALLOCATE (gammab(ni, nj))
165 200888 : rcut = rcuta + rcutb
166 803552 : dr = SQRT(SUM(rij(:)**2))
167 200888 : CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
168 1784616 : gchrg(iatom, 1:ni, 1) = gchrg(iatom, 1:ni, 1) + MATMUL(gammab, charges1(jatom, 1:nj))
169 200888 : IF (iatom /= jatom) THEN
170 1734442 : gchrg(jatom, 1:nj, 1) = gchrg(jatom, 1:nj, 1) + MATMUL(charges1(iatom, 1:ni), gammab)
171 : END IF
172 803552 : DEALLOCATE (gammab)
173 : END DO
174 250 : CALL neighbor_list_iterator_release(nl_iterator)
175 :
176 : ! 1/R contribution
177 :
178 250 : IF (xtb_control%coulomb_lr) THEN
179 250 : do_ewald = xtb_control%do_ewald
180 250 : IF (do_ewald) THEN
181 : ! Ewald sum
182 98 : NULLIFY (ewald_env, ewald_pw)
183 98 : NULLIFY (virial)
184 : CALL get_qs_env(qs_env=qs_env, &
185 98 : ewald_env=ewald_env, ewald_pw=ewald_pw)
186 98 : CALL get_cell(cell=cell, periodic=periodic, deth=deth)
187 98 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
188 98 : CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
189 98 : CALL tb_ewald_overlap(gmcharge, mcharge1, alpha, n_list, virial, .FALSE.)
190 0 : SELECT CASE (ewald_type)
191 : CASE DEFAULT
192 0 : CPABORT("Invalid Ewald type")
193 : CASE (do_ewald_none)
194 0 : CPABORT("Not allowed with DFTB")
195 : CASE (do_ewald_ewald)
196 0 : CPABORT("Standard Ewald not implemented in DFTB")
197 : CASE (do_ewald_pme)
198 0 : CPABORT("PME not implemented in DFTB")
199 : CASE (do_ewald_spme)
200 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
201 98 : gmcharge, mcharge1, .FALSE., virial, .FALSE.)
202 : END SELECT
203 : ELSE
204 : ! direct sum
205 : CALL get_qs_env(qs_env=qs_env, &
206 152 : local_particles=local_particles)
207 528 : DO ikind = 1, SIZE(local_particles%n_el)
208 792 : DO ia = 1, local_particles%n_el(ikind)
209 264 : iatom = local_particles%list(ikind)%array(ia)
210 976 : DO jatom = 1, iatom - 1
211 1344 : rij = particle_set(iatom)%r - particle_set(jatom)%r
212 1344 : rij = pbc(rij, cell)
213 1344 : dr = SQRT(SUM(rij(:)**2))
214 600 : IF (dr > 1.e-6_dp) THEN
215 336 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge1(jatom)/dr
216 336 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge1(iatom)/dr
217 : END IF
218 : END DO
219 : END DO
220 : END DO
221 : END IF
222 : END IF
223 :
224 : ! global sum of gamma*p arrays
225 250 : CALL get_qs_env(qs_env=qs_env, para_env=para_env)
226 250 : CALL para_env%sum(gmcharge(:, 1))
227 250 : CALL para_env%sum(gchrg(:, :, 1))
228 :
229 250 : IF (xtb_control%coulomb_lr) THEN
230 250 : IF (do_ewald) THEN
231 : ! add self charge interaction and background charge contribution
232 1918 : gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge1(:)
233 146 : IF (ANY(periodic(:) == 1)) THEN
234 1838 : gmcharge(:, 1) = gmcharge(:, 1) - pi/alpha**2/deth
235 : END IF
236 : END IF
237 : END IF
238 :
239 250 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
240 250 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
241 :
242 : ! no k-points; all matrices have been transformed to periodic bsf
243 250 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
244 38236 : DO WHILE (dbcsr_iterator_blocks_left(iter))
245 37986 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
246 37986 : ikind = kind_of(irow)
247 37986 : jkind = kind_of(icol)
248 :
249 : ! atomic parameters
250 37986 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
251 37986 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
252 37986 : CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
253 37986 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
254 :
255 37986 : ni = SIZE(sblock, 1)
256 37986 : nj = SIZE(sblock, 2)
257 151944 : ALLOCATE (gcij(ni, nj))
258 155948 : DO i = 1, ni
259 426008 : DO j = 1, nj
260 270060 : la = laoa(i) + 1
261 270060 : lb = laob(j) + 1
262 388022 : gcij(i, j) = gchrg(irow, la, 1) + gchrg(icol, lb, 1)
263 : END DO
264 : END DO
265 37986 : gmij = gmcharge(irow, 1) + gmcharge(icol, 1)
266 76412 : DO is = 1, SIZE(ks_matrix)
267 38426 : NULLIFY (ksblock)
268 : CALL dbcsr_get_block_p(matrix=ks_matrix(is)%matrix, &
269 38426 : row=irow, col=icol, block=ksblock, found=found)
270 38426 : CPASSERT(found)
271 757046 : ksblock = ksblock - gcij*sblock
272 833458 : ksblock = ksblock - gmij*sblock
273 : END DO
274 114208 : DEALLOCATE (gcij)
275 : END DO
276 250 : CALL dbcsr_iterator_stop(iter)
277 :
278 250 : IF (xtb_control%tb3_interaction) THEN
279 250 : CALL get_qs_env(qs_env, nkind=nkind)
280 750 : ALLOCATE (xgamma(nkind))
281 860 : DO ikind = 1, nkind
282 610 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
283 860 : CALL get_xtb_atom_param(xtb_kind, xgamma=xgamma(ikind))
284 : END DO
285 : ! Diagonal 3rd order correction (DFTB3)
286 250 : CALL dftb3_diagonal_hessian(qs_env, ks_matrix, mcharge, mcharge1, xgamma)
287 250 : DEALLOCATE (xgamma)
288 : END IF
289 :
290 250 : IF (xtb_control%do_spinpol) THEN
291 34 : CALL xtb_spinpol_hessian(qs_env, ks_matrix, matrix_p1)
292 : END IF
293 :
294 250 : IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
295 0 : CPABORT("QMMM not available in xTB response calculations")
296 : END IF
297 :
298 250 : DEALLOCATE (gmcharge, gchrg)
299 :
300 250 : CALL timestop(handle)
301 :
302 750 : END SUBROUTINE xtb_coulomb_hessian
303 :
304 : ! **************************************************************************************************
305 : !> \brief ...
306 : !> \param qs_env ...
307 : !> \param ks_matrix ...
308 : !> \param mcharge ...
309 : !> \param mcharge1 ...
310 : !> \param xgamma ...
311 : ! **************************************************************************************************
312 250 : SUBROUTINE dftb3_diagonal_hessian(qs_env, ks_matrix, mcharge, mcharge1, xgamma)
313 :
314 : TYPE(qs_environment_type), POINTER :: qs_env
315 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix
316 : REAL(dp), DIMENSION(:) :: mcharge, mcharge1, xgamma
317 :
318 : CHARACTER(len=*), PARAMETER :: routineN = 'dftb3_diagonal_hessian'
319 :
320 : INTEGER :: handle, icol, ikind, irow, is, jkind
321 250 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
322 : LOGICAL :: found
323 : REAL(KIND=dp) :: gmij, ui, uj
324 250 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksblock, sblock
325 250 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
326 : TYPE(dbcsr_iterator_type) :: iter
327 250 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
328 250 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
329 :
330 250 : CALL timeset(routineN, handle)
331 :
332 250 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
333 250 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
334 250 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
335 : ! no k-points; all matrices have been transformed to periodic bsf
336 250 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
337 38236 : DO WHILE (dbcsr_iterator_blocks_left(iter))
338 37986 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
339 37986 : ikind = kind_of(irow)
340 37986 : ui = xgamma(ikind)
341 37986 : jkind = kind_of(icol)
342 37986 : uj = xgamma(jkind)
343 37986 : gmij = ui*mcharge(irow)*mcharge1(irow) + uj*mcharge(icol)*mcharge1(icol)
344 76662 : DO is = 1, SIZE(ks_matrix)
345 38426 : NULLIFY (ksblock)
346 : CALL dbcsr_get_block_p(matrix=ks_matrix(is)%matrix, &
347 38426 : row=irow, col=icol, block=ksblock, found=found)
348 38426 : CPASSERT(found)
349 871884 : ksblock = ksblock + gmij*sblock
350 : END DO
351 : END DO
352 250 : CALL dbcsr_iterator_stop(iter)
353 :
354 250 : CALL timestop(handle)
355 :
356 500 : END SUBROUTINE dftb3_diagonal_hessian
357 :
358 395726 : END MODULE xtb_ehess
359 :
|
