Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 56 : 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 56 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
91 56 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
92 :
93 280 : qresp = 0.0_dp
94 56 : SELECT CASE (gfn_type)
95 : CASE (0)
96 56 : 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 56 : CPABORT("QRESP: Unknown gfn_type")
103 : END SELECT
104 :
105 56 : END SUBROUTINE build_xtb_qresp
106 :
107 : ! **************************************************************************************************
108 : !> \brief ...
109 : !> \param qs_env ...
110 : !> \param qresp ...
111 : ! **************************************************************************************************
112 56 : 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 56 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
123 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
124 : INTEGER, DIMENSION(3) :: cell
125 56 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
126 56 : npgfb, nsgfa, nsgfb
127 56 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
128 56 : 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 56 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charges, cnumbers, dcharges
134 56 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dhuckel, dqhuckel, huckel, owork
135 56 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
136 56 : 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 56 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
140 56 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock, rpgfa, rpgfb, scon_a, scon_b, &
141 56 : zeta, zetb
142 56 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
143 : TYPE(cp_logger_type), POINTER :: logger
144 56 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_w
145 56 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
146 : TYPE(dft_control_type), POINTER :: dft_control
147 : TYPE(eeq_solver_type) :: eeq_sparam
148 56 : 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 56 : DIMENSION(:), POINTER :: nl_iterator
154 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
155 56 : POINTER :: sab_orb
156 56 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
157 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
158 : TYPE(qs_energy_type), POINTER :: energy
159 56 : 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 56 : CALL timeset(routineN, handle)
166 :
167 56 : NULLIFY (logger)
168 56 : logger => cp_get_default_logger()
169 :
170 56 : 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 56 : sab_orb=sab_orb)
179 :
180 56 : nkind = SIZE(atomic_kind_set)
181 56 : xtb_control => dft_control%qs_control%xtb_control
182 56 : do_nonbonded = xtb_control%do_nonbonded
183 56 : nimg = dft_control%nimages
184 56 : nderivatives = 0
185 56 : maxder = ncoset(nderivatives)
186 :
187 56 : NULLIFY (particle_set)
188 56 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
189 56 : natom = SIZE(particle_set)
190 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
191 56 : atom_of_kind=atom_of_kind, kind_of=kind_of)
192 :
193 56 : NULLIFY (rho, matrix_w)
194 56 : CALL get_qs_env(qs_env=qs_env, rho=rho)
195 56 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
196 :
197 56 : IF (SIZE(matrix_p, 1) == 2) THEN
198 0 : DO img = 1, nimg
199 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
200 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
201 : CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
202 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
203 : END DO
204 : END IF
205 :
206 56 : NULLIFY (cell_to_index)
207 56 : IF (nimg > 1) THEN
208 0 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
209 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
210 : END IF
211 :
212 : ! set up basis set lists
213 336 : ALLOCATE (basis_set_list(nkind))
214 56 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
215 :
216 : ! Calculate coordination numbers
217 : ! needed for effective atomic energy levels
218 : ! code taken from D3 dispersion energy
219 56 : CALL cnumber_init(qs_env, cnumbers, dcnum, 2, .TRUE.)
220 :
221 224 : ALLOCATE (charges(natom), dcharges(natom))
222 280 : charges = 0.0_dp
223 56 : CALL xtb_eeq_calculation(qs_env, charges, cnumbers, eeq_sparam, eeq_energy, ef_energy, qlambda)
224 280 : dcharges = qlambda/REAL(para_env%num_pe, KIND=dp)
225 :
226 56 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
227 56 : IF (dispersion_env%pp_type == vdw_pairpot_dftd4 .AND. dispersion_env%ext_charges) THEN
228 28 : IF (ASSOCIATED(dispersion_env%dcharges)) THEN
229 120 : dcharges(1:natom) = dcharges(1:natom) + dispersion_env%dcharges(1:natom)
230 : ELSE
231 4 : CPWARN("gfn0-xTB variational dipole DFTD4 contribution missing")
232 : END IF
233 : END IF
234 :
235 : ! Calculate Huckel parameters
236 56 : CALL gfn0_huckel(qs_env, cnumbers, charges, huckel, dhuckel, dqhuckel, .TRUE.)
237 :
238 : ! Calculate KAB parameters and electronegativity correction
239 56 : CALL gfn0_kpair(qs_env, kijab)
240 :
241 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
242 56 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
243 350 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
244 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
245 294 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
246 294 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
247 294 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
248 294 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
249 294 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
250 294 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
251 294 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
252 :
253 1176 : dr = SQRT(SUM(rij(:)**2))
254 :
255 : ! atomic parameters
256 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
257 294 : lmax=lmaxa, nshell=nsa, kpoly=kpolya, hen=hena)
258 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
259 294 : lmax=lmaxb, nshell=nsb, kpoly=kpolyb, hen=henb)
260 :
261 294 : IF (nimg == 1) THEN
262 : ic = 1
263 : ELSE
264 0 : ic = cell_to_index(cell(1), cell(2), cell(3))
265 0 : CPASSERT(ic > 0)
266 : END IF
267 :
268 294 : icol = MAX(iatom, jatom)
269 294 : irow = MIN(iatom, jatom)
270 :
271 294 : NULLIFY (pblock)
272 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
273 294 : row=irow, col=icol, block=pblock, found=found)
274 294 : CPASSERT(ASSOCIATED(pblock))
275 :
276 : ! overlap
277 294 : basis_set_a => basis_set_list(ikind)%gto_basis_set
278 294 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
279 294 : basis_set_b => basis_set_list(jkind)%gto_basis_set
280 294 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
281 294 : atom_a = atom_of_kind(iatom)
282 294 : atom_b = atom_of_kind(jatom)
283 : ! basis ikind
284 294 : first_sgfa => basis_set_a%first_sgf
285 294 : la_max => basis_set_a%lmax
286 294 : la_min => basis_set_a%lmin
287 294 : npgfa => basis_set_a%npgf
288 294 : nseta = basis_set_a%nset
289 294 : nsgfa => basis_set_a%nsgf_set
290 294 : rpgfa => basis_set_a%pgf_radius
291 294 : set_radius_a => basis_set_a%set_radius
292 294 : scon_a => basis_set_a%scon
293 294 : zeta => basis_set_a%zet
294 : ! basis jkind
295 294 : first_sgfb => basis_set_b%first_sgf
296 294 : lb_max => basis_set_b%lmax
297 294 : lb_min => basis_set_b%lmin
298 294 : npgfb => basis_set_b%npgf
299 294 : nsetb = basis_set_b%nset
300 294 : nsgfb => basis_set_b%nsgf_set
301 294 : rpgfb => basis_set_b%pgf_radius
302 294 : set_radius_b => basis_set_b%set_radius
303 294 : scon_b => basis_set_b%scon
304 294 : zetb => basis_set_b%zet
305 :
306 294 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
307 2352 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
308 1470 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
309 4088 : sint = 0.0_dp
310 :
311 882 : DO iset = 1, nseta
312 588 : ncoa = npgfa(iset)*ncoset(la_max(iset))
313 588 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
314 588 : sgfa = first_sgfa(1, iset)
315 2058 : DO jset = 1, nsetb
316 1176 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
317 1176 : ncob = npgfb(jset)*ncoset(lb_max(jset))
318 1176 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
319 1176 : sgfb = first_sgfb(1, jset)
320 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
321 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
322 1176 : rij, sab=oint(:, :, 1))
323 : ! Contraction
324 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
325 1176 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
326 1764 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
327 : END DO
328 : END DO
329 :
330 : ! Calculate Pi = Pia * Pib (Eq. 11)
331 294 : rcovab = rcova + rcovb
332 294 : rrab = SQRT(dr/rcovab)
333 882 : pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
334 882 : pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
335 :
336 : ! diagonal block
337 294 : diagblock = .FALSE.
338 294 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
339 : !
340 294 : f2 = 1.0_dp
341 294 : IF (iatom /= jatom) f2 = 2.0_dp
342 : ! Derivative wrt coordination number
343 294 : fqa = 0.0_dp
344 294 : fqb = 0.0_dp
345 294 : IF (diagblock) THEN
346 448 : DO i = 1, natorb_a
347 336 : na = naoa(i)
348 448 : fqa = fqa + pblock(i, i)*dqhuckel(na, iatom)
349 : END DO
350 112 : dcharges(iatom) = dcharges(iatom) + fqa
351 : ELSE
352 714 : DO j = 1, natorb_b
353 532 : nb = naob(j)
354 2226 : DO i = 1, natorb_a
355 1512 : na = naoa(i)
356 1512 : hij = 0.5_dp*pia(na)*pib(nb)
357 1512 : drx = f2*hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)
358 2044 : IF (iatom <= jatom) THEN
359 448 : fqa = fqa + drx*pblock(i, j)*dqhuckel(na, iatom)
360 448 : fqb = fqb + drx*pblock(i, j)*dqhuckel(nb, jatom)
361 : ELSE
362 1064 : fqa = fqa + drx*pblock(j, i)*dqhuckel(na, iatom)
363 1064 : fqb = fqb + drx*pblock(j, i)*dqhuckel(nb, jatom)
364 : END IF
365 : END DO
366 : END DO
367 182 : dcharges(iatom) = dcharges(iatom) + fqa
368 182 : dcharges(jatom) = dcharges(jatom) + fqb
369 : END IF
370 :
371 1526 : DEALLOCATE (oint, owork, sint)
372 :
373 : END DO
374 56 : CALL neighbor_list_iterator_release(nl_iterator)
375 :
376 : ! EEQ response
377 56 : CALL para_env%sum(dcharges)
378 56 : CALL xtb_eeq_lagrange(qs_env, dcharges, qresp, eeq_sparam)
379 :
380 : ! deallocate coordination numbers
381 56 : CALL cnumber_release(cnumbers, dcnum, .TRUE.)
382 :
383 : ! deallocate Huckel parameters
384 56 : DEALLOCATE (huckel, dhuckel, dqhuckel)
385 : ! deallocate KAB parameters
386 56 : DEALLOCATE (kijab)
387 :
388 : ! deallocate charges
389 56 : DEALLOCATE (charges, dcharges)
390 :
391 56 : DEALLOCATE (basis_set_list)
392 56 : IF (SIZE(matrix_p, 1) == 2) THEN
393 0 : DO img = 1, nimg
394 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
395 0 : beta_scalar=-1.0_dp)
396 : END DO
397 : END IF
398 :
399 56 : CALL timestop(handle)
400 :
401 112 : END SUBROUTINE build_gfn0_qresp
402 :
403 : END MODULE xtb_qresp
404 :
|
