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 charge response in xTB (EEQ only)
10 : !> Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov
11 : !> JCTC 13, 1989-2009, (2017)
12 : !> DOI: 10.1021/acs.jctc.7b00118
13 : !> \author JGH
14 : ! **************************************************************************************************
15 : MODULE xtb_qresp
16 : USE ai_contraction, ONLY: block_add,&
17 : contraction
18 : USE ai_overlap, ONLY: overlap_ab
19 : USE atomic_kind_types, ONLY: atomic_kind_type,&
20 : get_atomic_kind_set
21 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
22 : gto_basis_set_type
23 : USE cp_control_types, ONLY: dft_control_type,&
24 : xtb_control_type
25 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
26 : dbcsr_get_block_p,&
27 : dbcsr_p_type
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_type
30 : USE eeq_input, ONLY: eeq_solver_type
31 : USE input_constants, ONLY: vdw_pairpot_dftd4
32 : USE kinds, ONLY: dp
33 : USE kpoint_types, ONLY: get_kpoint_info,&
34 : kpoint_type
35 : USE message_passing, ONLY: mp_para_env_type
36 : USE orbital_pointers, ONLY: ncoset
37 : USE particle_types, ONLY: particle_type
38 : USE qs_dispersion_cnum, ONLY: cnumber_init,&
39 : cnumber_release,&
40 : dcnum_type
41 : USE qs_dispersion_types, ONLY: qs_dispersion_type
42 : USE qs_energy_types, ONLY: qs_energy_type
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_integral_utils, ONLY: basis_set_list_setup,&
46 : get_memory_usage
47 : USE qs_kind_types, ONLY: get_qs_kind,&
48 : qs_kind_type
49 : USE qs_ks_types, ONLY: get_ks_env,&
50 : qs_ks_env_type
51 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
52 : neighbor_list_iterate,&
53 : neighbor_list_iterator_create,&
54 : neighbor_list_iterator_p_type,&
55 : neighbor_list_iterator_release,&
56 : neighbor_list_set_p_type
57 : USE qs_rho_types, ONLY: qs_rho_get,&
58 : qs_rho_type
59 : USE xtb_eeq, ONLY: xtb_eeq_calculation,&
60 : xtb_eeq_lagrange
61 : USE xtb_hcore, ONLY: gfn0_huckel,&
62 : gfn0_kpair
63 : USE xtb_types, ONLY: get_xtb_atom_param,&
64 : xtb_atom_type
65 : #include "./base/base_uses.f90"
66 :
67 : IMPLICIT NONE
68 :
69 : PRIVATE
70 :
71 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_qresp'
72 :
73 : PUBLIC :: build_xtb_qresp
74 :
75 : CONTAINS
76 :
77 : ! **************************************************************************************************
78 : !> \brief ...
79 : !> \param qs_env ...
80 : !> \param qresp ...
81 : ! **************************************************************************************************
82 58 : SUBROUTINE build_xtb_qresp(qs_env, qresp)
83 :
84 : TYPE(qs_environment_type), POINTER :: qs_env
85 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: qresp
86 :
87 : INTEGER :: gfn_type
88 : TYPE(dft_control_type), POINTER :: dft_control
89 :
90 58 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
91 58 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
92 :
93 290 : qresp = 0.0_dp
94 58 : SELECT CASE (gfn_type)
95 : CASE (0)
96 58 : CALL build_gfn0_qresp(qs_env, qresp)
97 : CASE (1)
98 0 : CPABORT("QRESP: gfn_type = 1 not available")
99 : CASE (2)
100 0 : CPABORT("QRESP: gfn_type = 2 not available")
101 : CASE DEFAULT
102 58 : CPABORT("QRESP: Unknown gfn_type")
103 : END SELECT
104 :
105 58 : END SUBROUTINE build_xtb_qresp
106 :
107 : ! **************************************************************************************************
108 : !> \brief ...
109 : !> \param qs_env ...
110 : !> \param qresp ...
111 : ! **************************************************************************************************
112 58 : SUBROUTINE build_gfn0_qresp(qs_env, qresp)
113 :
114 : TYPE(qs_environment_type), POINTER :: qs_env
115 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: qresp
116 :
117 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_gfn0_qresp'
118 :
119 : INTEGER :: atom_a, atom_b, handle, i, iatom, ic, icol, ikind, img, irow, iset, j, jatom, &
120 : jkind, jset, ldsab, lmaxa, lmaxb, maxder, n1, n2, na, natom, natorb_a, natorb_b, nb, &
121 : ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, nsetb, sgfa, sgfb, za, zb
122 58 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
123 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
124 : INTEGER, DIMENSION(3) :: cell
125 58 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
126 58 : npgfb, nsgfa, nsgfb
127 58 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
128 58 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
129 : LOGICAL :: defined, diagblock, do_nonbonded, found
130 : REAL(KIND=dp) :: dr, drx, eeq_energy, ef_energy, etaa, &
131 : etab, f2, fqa, fqb, hij, qlambda, &
132 : rcova, rcovab, rcovb, rrab
133 58 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charges, cnumbers, dcharges
134 58 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dhuckel, dqhuckel, huckel, owork
135 58 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
136 58 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
137 : REAL(KIND=dp), DIMENSION(3) :: rij
138 : REAL(KIND=dp), DIMENSION(5) :: hena, henb, kpolya, kpolyb, pia, pib
139 58 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
140 58 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock, rpgfa, rpgfb, scon_a, scon_b, &
141 58 : zeta, zetb
142 58 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
143 : TYPE(cp_logger_type), POINTER :: logger
144 58 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
145 58 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
146 : TYPE(dft_control_type), POINTER :: dft_control
147 : TYPE(eeq_solver_type) :: eeq_sparam
148 58 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
149 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
150 : TYPE(kpoint_type), POINTER :: kpoints
151 : TYPE(mp_para_env_type), POINTER :: para_env
152 : TYPE(neighbor_list_iterator_p_type), &
153 58 : DIMENSION(:), POINTER :: nl_iterator
154 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
155 58 : POINTER :: sab_orb
156 58 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
157 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
158 : TYPE(qs_energy_type), POINTER :: energy
159 58 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
160 : TYPE(qs_ks_env_type), POINTER :: ks_env
161 : TYPE(qs_rho_type), POINTER :: rho
162 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
163 : TYPE(xtb_control_type), POINTER :: xtb_control
164 :
165 58 : CALL timeset(routineN, handle)
166 :
167 58 : NULLIFY (logger)
168 58 : logger => cp_get_default_logger()
169 :
170 58 : NULLIFY (matrix_p, atomic_kind_set, qs_kind_set, sab_orb, ks_env)
171 : CALL get_qs_env(qs_env=qs_env, &
172 : ks_env=ks_env, &
173 : energy=energy, &
174 : atomic_kind_set=atomic_kind_set, &
175 : qs_kind_set=qs_kind_set, &
176 : para_env=para_env, &
177 : dft_control=dft_control, &
178 58 : sab_orb=sab_orb)
179 :
180 58 : nkind = SIZE(atomic_kind_set)
181 58 : xtb_control => dft_control%qs_control%xtb_control
182 58 : do_nonbonded = xtb_control%do_nonbonded
183 58 : nimg = dft_control%nimages
184 58 : nderivatives = 0
185 58 : maxder = ncoset(nderivatives)
186 :
187 58 : NULLIFY (particle_set)
188 58 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
189 58 : natom = SIZE(particle_set)
190 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
191 58 : atom_of_kind=atom_of_kind, kind_of=kind_of)
192 :
193 58 : NULLIFY (rho)
194 58 : CALL get_qs_env(qs_env=qs_env, rho=rho)
195 58 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
196 :
197 58 : IF (SIZE(matrix_p, 1) == 2) THEN
198 4 : DO img = 1, nimg
199 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
200 4 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
201 : END DO
202 : END IF
203 :
204 58 : NULLIFY (cell_to_index)
205 58 : IF (nimg > 1) THEN
206 0 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
207 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
208 : END IF
209 :
210 : ! set up basis set lists
211 348 : ALLOCATE (basis_set_list(nkind))
212 58 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
213 :
214 : ! Calculate coordination numbers
215 : ! needed for effective atomic energy levels
216 : ! code taken from D3 dispersion energy
217 58 : CALL cnumber_init(qs_env, cnumbers, dcnum, 2, .TRUE.)
218 :
219 232 : ALLOCATE (charges(natom), dcharges(natom))
220 58 : charges = 0.0_dp
221 58 : CALL xtb_eeq_calculation(qs_env, charges, cnumbers, eeq_sparam, eeq_energy, ef_energy, qlambda)
222 290 : dcharges = qlambda/REAL(para_env%num_pe, KIND=dp)
223 :
224 58 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
225 58 : IF (dispersion_env%pp_type == vdw_pairpot_dftd4 .AND. dispersion_env%ext_charges) THEN
226 28 : IF (ASSOCIATED(dispersion_env%dcharges)) THEN
227 120 : dcharges(1:natom) = dcharges(1:natom) + dispersion_env%dcharges(1:natom)
228 : ELSE
229 4 : CPWARN("gfn0-xTB variational dipole DFTD4 contribution missing")
230 : END IF
231 : END IF
232 :
233 : ! Calculate Huckel parameters
234 58 : CALL gfn0_huckel(qs_env, cnumbers, charges, huckel, dhuckel, dqhuckel, .TRUE.)
235 :
236 : ! Calculate KAB parameters and electronegativity correction
237 58 : CALL gfn0_kpair(qs_env, kijab)
238 :
239 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
240 58 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
241 362 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
242 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
243 304 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
244 304 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
245 304 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
246 304 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
247 304 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
248 304 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
249 304 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
250 :
251 1216 : dr = SQRT(SUM(rij(:)**2))
252 :
253 : ! atomic parameters
254 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
255 304 : lmax=lmaxa, nshell=nsa, kpoly=kpolya, hen=hena)
256 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
257 304 : lmax=lmaxb, nshell=nsb, kpoly=kpolyb, hen=henb)
258 :
259 304 : IF (nimg == 1) THEN
260 : ic = 1
261 : ELSE
262 0 : ic = cell_to_index(cell(1), cell(2), cell(3))
263 0 : CPASSERT(ic > 0)
264 : END IF
265 :
266 304 : icol = MAX(iatom, jatom)
267 304 : irow = MIN(iatom, jatom)
268 :
269 304 : NULLIFY (pblock)
270 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
271 304 : row=irow, col=icol, block=pblock, found=found)
272 304 : CPASSERT(ASSOCIATED(pblock))
273 :
274 : ! overlap
275 304 : basis_set_a => basis_set_list(ikind)%gto_basis_set
276 304 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
277 304 : basis_set_b => basis_set_list(jkind)%gto_basis_set
278 304 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
279 304 : atom_a = atom_of_kind(iatom)
280 304 : atom_b = atom_of_kind(jatom)
281 : ! basis ikind
282 304 : first_sgfa => basis_set_a%first_sgf
283 304 : la_max => basis_set_a%lmax
284 304 : la_min => basis_set_a%lmin
285 304 : npgfa => basis_set_a%npgf
286 304 : nseta = basis_set_a%nset
287 304 : nsgfa => basis_set_a%nsgf_set
288 304 : rpgfa => basis_set_a%pgf_radius
289 304 : set_radius_a => basis_set_a%set_radius
290 304 : scon_a => basis_set_a%scon
291 304 : zeta => basis_set_a%zet
292 : ! basis jkind
293 304 : first_sgfb => basis_set_b%first_sgf
294 304 : lb_max => basis_set_b%lmax
295 304 : lb_min => basis_set_b%lmin
296 304 : npgfb => basis_set_b%npgf
297 304 : nsetb = basis_set_b%nset
298 304 : nsgfb => basis_set_b%nsgf_set
299 304 : rpgfb => basis_set_b%pgf_radius
300 304 : set_radius_b => basis_set_b%set_radius
301 304 : scon_b => basis_set_b%scon
302 304 : zetb => basis_set_b%zet
303 :
304 304 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
305 2432 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
306 1520 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
307 304 : sint = 0.0_dp
308 :
309 912 : DO iset = 1, nseta
310 608 : ncoa = npgfa(iset)*ncoset(la_max(iset))
311 608 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
312 608 : sgfa = first_sgfa(1, iset)
313 2128 : DO jset = 1, nsetb
314 1216 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
315 1216 : ncob = npgfb(jset)*ncoset(lb_max(jset))
316 1216 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
317 1216 : sgfb = first_sgfb(1, jset)
318 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
319 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
320 1216 : rij, sab=oint(:, :, 1))
321 : ! Contraction
322 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
323 1216 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
324 1824 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
325 : END DO
326 : END DO
327 :
328 : ! Calculate Pi = Pia * Pib (Eq. 11)
329 304 : rcovab = rcova + rcovb
330 304 : rrab = SQRT(dr/rcovab)
331 912 : pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
332 912 : pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
333 :
334 : ! diagonal block
335 304 : diagblock = .FALSE.
336 304 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
337 : !
338 304 : f2 = 1.0_dp
339 304 : IF (iatom /= jatom) f2 = 2.0_dp
340 : ! Derivative wrt coordination number
341 304 : fqa = 0.0_dp
342 304 : fqb = 0.0_dp
343 304 : IF (diagblock) THEN
344 464 : DO i = 1, natorb_a
345 348 : na = naoa(i)
346 464 : fqa = fqa + pblock(i, i)*dqhuckel(na, iatom)
347 : END DO
348 116 : dcharges(iatom) = dcharges(iatom) + fqa
349 : ELSE
350 738 : DO j = 1, natorb_b
351 550 : nb = naob(j)
352 2302 : DO i = 1, natorb_a
353 1564 : na = naoa(i)
354 1564 : hij = 0.5_dp*pia(na)*pib(nb)
355 1564 : drx = f2*hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)
356 2114 : IF (iatom <= jatom) THEN
357 464 : fqa = fqa + drx*pblock(i, j)*dqhuckel(na, iatom)
358 464 : fqb = fqb + drx*pblock(i, j)*dqhuckel(nb, jatom)
359 : ELSE
360 1100 : fqa = fqa + drx*pblock(j, i)*dqhuckel(na, iatom)
361 1100 : fqb = fqb + drx*pblock(j, i)*dqhuckel(nb, jatom)
362 : END IF
363 : END DO
364 : END DO
365 188 : dcharges(iatom) = dcharges(iatom) + fqa
366 188 : dcharges(jatom) = dcharges(jatom) + fqb
367 : END IF
368 :
369 1578 : DEALLOCATE (oint, owork, sint)
370 :
371 : END DO
372 58 : CALL neighbor_list_iterator_release(nl_iterator)
373 :
374 : ! EEQ response
375 58 : CALL para_env%sum(dcharges)
376 58 : CALL xtb_eeq_lagrange(qs_env, dcharges, qresp, eeq_sparam)
377 :
378 : ! deallocate coordination numbers
379 58 : CALL cnumber_release(cnumbers, dcnum, .TRUE.)
380 :
381 : ! deallocate Huckel parameters
382 58 : DEALLOCATE (huckel, dhuckel, dqhuckel)
383 : ! deallocate KAB parameters
384 58 : DEALLOCATE (kijab)
385 :
386 : ! deallocate charges
387 58 : DEALLOCATE (charges, dcharges)
388 :
389 58 : DEALLOCATE (basis_set_list)
390 58 : IF (SIZE(matrix_p, 1) == 2) THEN
391 4 : DO img = 1, nimg
392 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
393 4 : beta_scalar=-1.0_dp)
394 : END DO
395 : END IF
396 :
397 58 : CALL timestop(handle)
398 :
399 116 : END SUBROUTINE build_gfn0_qresp
400 :
401 : END MODULE xtb_qresp
402 :
|
