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 contributions in xTB
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE xtb_coulomb
13 : USE ai_contraction, ONLY: block_add,&
14 : contraction
15 : USE ai_overlap, ONLY: overlap_ab
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind_set
18 : USE atprop_types, ONLY: atprop_array_init,&
19 : atprop_type
20 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
21 : gto_basis_set_type
22 : USE cell_types, ONLY: cell_type,&
23 : get_cell,&
24 : pbc
25 : USE cp_control_types, ONLY: dft_control_type,&
26 : xtb_control_type
27 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
28 : dbcsr_get_block_p,&
29 : dbcsr_iterator_blocks_left,&
30 : dbcsr_iterator_next_block,&
31 : dbcsr_iterator_start,&
32 : dbcsr_iterator_stop,&
33 : dbcsr_iterator_type,&
34 : dbcsr_p_type
35 : USE distribution_1d_types, ONLY: distribution_1d_type
36 : USE ewald_environment_types, ONLY: ewald_env_get,&
37 : ewald_environment_type
38 : USE ewald_methods_tb, ONLY: tb_ewald_overlap,&
39 : tb_spme_evaluate
40 : USE ewald_pw_types, ONLY: ewald_pw_type
41 : USE kinds, ONLY: dp
42 : USE kpoint_types, ONLY: get_kpoint_info,&
43 : kpoint_type
44 : USE mathconstants, ONLY: oorootpi,&
45 : pi
46 : USE message_passing, ONLY: mp_para_env_type
47 : USE orbital_pointers, ONLY: ncoset
48 : USE particle_types, ONLY: particle_type
49 : USE pw_poisson_types, ONLY: do_ewald_ewald,&
50 : do_ewald_none,&
51 : do_ewald_pme,&
52 : do_ewald_spme
53 : USE qmmm_tb_coulomb, ONLY: build_tb_coulomb_qmqm
54 : USE qs_dftb3_methods, ONLY: build_dftb3_diagonal
55 : USE qs_energy_types, ONLY: qs_energy_type
56 : USE qs_environment_types, ONLY: get_qs_env,&
57 : qs_environment_type
58 : USE qs_force_types, ONLY: qs_force_type
59 : USE qs_integral_utils, ONLY: basis_set_list_setup,&
60 : get_memory_usage
61 : USE qs_kind_types, ONLY: get_qs_kind,&
62 : qs_kind_type
63 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
64 : neighbor_list_iterate,&
65 : neighbor_list_iterator_create,&
66 : neighbor_list_iterator_p_type,&
67 : neighbor_list_iterator_release,&
68 : neighbor_list_set_p_type
69 : USE qs_rho_types, ONLY: qs_rho_get,&
70 : qs_rho_type
71 : USE sap_kind_types, ONLY: clist_type,&
72 : release_sap_int,&
73 : sap_int_type
74 : USE virial_methods, ONLY: virial_pair_force
75 : USE virial_types, ONLY: virial_type
76 : USE xtb_spinpol, ONLY: build_xtb_spinpol
77 : USE xtb_types, ONLY: get_xtb_atom_param,&
78 : xtb_atom_type
79 : #include "./base/base_uses.f90"
80 :
81 : IMPLICIT NONE
82 :
83 : PRIVATE
84 :
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_coulomb'
86 :
87 : PUBLIC :: build_xtb_coulomb, gamma_rab_sr, dgamma_rab_sr, xtb_dsint_list
88 :
89 : CONTAINS
90 :
91 : ! **************************************************************************************************
92 : !> \brief ...
93 : !> \param qs_env ...
94 : !> \param ks_matrix ...
95 : !> \param rho ...
96 : !> \param charges ...
97 : !> \param mcharge ...
98 : !> \param energy ...
99 : !> \param calculate_forces ...
100 : !> \param just_energy ...
101 : ! **************************************************************************************************
102 31702 : SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
103 : calculate_forces, just_energy)
104 :
105 : TYPE(qs_environment_type), POINTER :: qs_env
106 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
107 : TYPE(qs_rho_type), POINTER :: rho
108 : REAL(dp), DIMENSION(:, :), INTENT(in) :: charges
109 : REAL(dp), DIMENSION(:), INTENT(in) :: mcharge
110 : TYPE(qs_energy_type), POINTER :: energy
111 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
112 :
113 : CHARACTER(len=*), PARAMETER :: routineN = 'build_xtb_coulomb'
114 :
115 : INTEGER :: atom_i, atom_j, ewald_type, handle, i, ia, iac, iatom, ic, icol, ikind, img, &
116 : irow, is, j, jatom, jkind, la, lb, lmaxa, lmaxb, natom, natorb_a, natorb_b, ni, nimg, nj, &
117 : nkind, nmat, za, zb
118 31702 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
119 : INTEGER, DIMENSION(25) :: laoa, laob
120 : INTEGER, DIMENSION(3) :: cellind, periodic
121 : INTEGER, DIMENSION(5) :: occ
122 31702 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
123 : LOGICAL :: defined, do_ewald, do_gamma_stress, &
124 : found, use_virial
125 : REAL(KIND=dp) :: alpha, deth, dr, ecsr, etaa, etab, f1, &
126 : f2, fi, gmij, kg, rcut, rcuta, rcutb, &
127 : zeff
128 31702 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xgamma, zeffk
129 31702 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gammab, gcij, gmcharge
130 31702 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gchrg
131 : REAL(KIND=dp), DIMENSION(25) :: gcint
132 : REAL(KIND=dp), DIMENSION(3) :: fij, rij
133 : REAL(KIND=dp), DIMENSION(5) :: kappaa, kappab
134 31702 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, ksblock, pblock, sblock
135 31702 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
136 31702 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
137 : TYPE(atprop_type), POINTER :: atprop
138 : TYPE(cell_type), POINTER :: cell
139 : TYPE(dbcsr_iterator_type) :: iter
140 31702 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
141 : TYPE(dft_control_type), POINTER :: dft_control
142 : TYPE(distribution_1d_type), POINTER :: local_particles
143 : TYPE(ewald_environment_type), POINTER :: ewald_env
144 : TYPE(ewald_pw_type), POINTER :: ewald_pw
145 : TYPE(kpoint_type), POINTER :: kpoints
146 : TYPE(mp_para_env_type), POINTER :: para_env
147 : TYPE(neighbor_list_iterator_p_type), &
148 31702 : DIMENSION(:), POINTER :: nl_iterator
149 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
150 31702 : POINTER :: n_list
151 31702 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
152 31702 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
153 31702 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
154 31702 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
155 : TYPE(virial_type), POINTER :: virial
156 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind
157 : TYPE(xtb_control_type), POINTER :: xtb_control
158 :
159 31702 : CALL timeset(routineN, handle)
160 :
161 31702 : NULLIFY (matrix_p, matrix_s, virial, atprop, dft_control)
162 :
163 : CALL get_qs_env(qs_env, &
164 : qs_kind_set=qs_kind_set, &
165 : particle_set=particle_set, &
166 : cell=cell, &
167 : virial=virial, &
168 : atprop=atprop, &
169 31702 : dft_control=dft_control)
170 :
171 31702 : xtb_control => dft_control%qs_control%xtb_control
172 :
173 31702 : use_virial = .FALSE.
174 31702 : IF (calculate_forces) THEN
175 930 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
176 : END IF
177 :
178 31702 : do_gamma_stress = .FALSE.
179 31702 : IF (.NOT. just_energy .AND. use_virial) THEN
180 138 : IF (dft_control%nimages == 1) do_gamma_stress = .TRUE.
181 : END IF
182 :
183 31702 : IF (atprop%energy) THEN
184 172 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
185 172 : natom = SIZE(particle_set)
186 172 : CALL atprop_array_init(atprop%atecoul, natom)
187 : END IF
188 :
189 31702 : IF (calculate_forces) THEN
190 : nmat = 4
191 : ELSE
192 31168 : nmat = 1
193 : END IF
194 :
195 31702 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
196 126808 : ALLOCATE (gchrg(natom, 5, nmat))
197 31702 : gchrg = 0._dp
198 126808 : ALLOCATE (gmcharge(natom, nmat))
199 31702 : gmcharge = 0._dp
200 :
201 : ! short range contribution (gamma)
202 : ! loop over all atom pairs (sab_xtbe)
203 31702 : kg = xtb_control%kg
204 31702 : NULLIFY (n_list)
205 31702 : CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
206 31702 : IF (.NOT. ASSOCIATED(n_list)) THEN
207 0 : CPABORT("sab_xtbe neighbor list is not associated in build_xtb_coulomb")
208 : END IF
209 31702 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
210 9193355 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
211 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
212 9161653 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
213 9161653 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
214 9161653 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
215 9161653 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
216 9161646 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
217 9161646 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
218 9161646 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
219 : ! atomic parameters
220 9161639 : CALL get_xtb_atom_param(xtb_atom_a, eta=etaa, lmax=lmaxa, kappa=kappaa, rcut=rcuta)
221 9161639 : CALL get_xtb_atom_param(xtb_atom_b, eta=etab, lmax=lmaxb, kappa=kappab, rcut=rcutb)
222 : ! gamma matrix
223 9161639 : ni = lmaxa + 1
224 9161639 : nj = lmaxb + 1
225 36646556 : ALLOCATE (gammab(ni, nj))
226 9161639 : rcut = rcuta + rcutb
227 36646556 : dr = SQRT(SUM(rij(:)**2))
228 9161639 : CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
229 110018116 : gchrg(iatom, 1:ni, 1) = gchrg(iatom, 1:ni, 1) + MATMUL(gammab, charges(jatom, 1:nj))
230 9161639 : IF (iatom /= jatom) THEN
231 88537924 : gchrg(jatom, 1:nj, 1) = gchrg(jatom, 1:nj, 1) + MATMUL(charges(iatom, 1:ni), gammab)
232 : END IF
233 9161639 : IF (calculate_forces) THEN
234 362407 : IF (dr > 1.e-6_dp) THEN
235 360139 : CALL dgamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
236 1440556 : DO i = 1, 3
237 : gchrg(iatom, 1:ni, i + 1) = gchrg(iatom, 1:ni, i + 1) &
238 14561808 : + MATMUL(gammab, charges(jatom, 1:nj))*rij(i)/dr
239 1440556 : IF (iatom /= jatom) THEN
240 : gchrg(jatom, 1:nj, i + 1) = gchrg(jatom, 1:nj, i + 1) &
241 11601735 : - MATMUL(charges(iatom, 1:ni), gammab)*rij(i)/dr
242 : END IF
243 : END DO
244 360139 : IF (use_virial) THEN
245 2848827 : gcint(1:ni) = MATMUL(gammab, charges(jatom, 1:nj))
246 834912 : DO i = 1, 3
247 2056449 : fij(i) = -SUM(charges(iatom, 1:ni)*gcint(1:ni))*rij(i)/dr
248 : END DO
249 208728 : fi = 1.0_dp
250 208728 : IF (iatom == jatom) fi = 0.5_dp
251 208728 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
252 : END IF
253 : END IF
254 : END IF
255 36646577 : DEALLOCATE (gammab)
256 : END DO
257 31702 : CALL neighbor_list_iterator_release(nl_iterator)
258 :
259 : ! 1/R contribution
260 :
261 31702 : IF (xtb_control%coulomb_lr) THEN
262 31702 : do_ewald = xtb_control%do_ewald
263 31702 : IF (do_ewald) THEN
264 : ! Ewald sum
265 15654 : NULLIFY (ewald_env, ewald_pw)
266 : CALL get_qs_env(qs_env=qs_env, &
267 15654 : ewald_env=ewald_env, ewald_pw=ewald_pw)
268 15654 : CALL get_cell(cell=cell, periodic=periodic, deth=deth)
269 15654 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
270 15654 : CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
271 15654 : CALL tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial)
272 0 : SELECT CASE (ewald_type)
273 : CASE DEFAULT
274 0 : CPABORT("Invalid Ewald type")
275 : CASE (do_ewald_none)
276 0 : CPABORT("Not allowed with xTB/DFTB")
277 : CASE (do_ewald_ewald)
278 0 : CPABORT("Standard Ewald not implemented in xTB/DFTB")
279 : CASE (do_ewald_pme)
280 0 : CPABORT("PME not implemented in xTB/DFTB")
281 : CASE (do_ewald_spme)
282 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
283 15654 : gmcharge, mcharge, calculate_forces, virial, use_virial)
284 : END SELECT
285 : ELSE
286 : ! direct sum
287 : CALL get_qs_env(qs_env=qs_env, &
288 16048 : local_particles=local_particles)
289 59642 : DO ikind = 1, SIZE(local_particles%n_el)
290 135931 : DO ia = 1, local_particles%n_el(ikind)
291 76289 : iatom = local_particles%list(ikind)%array(ia)
292 899995 : DO jatom = 1, iatom - 1
293 3120448 : rij = particle_set(iatom)%r - particle_set(jatom)%r
294 3120448 : rij = pbc(rij, cell)
295 3120448 : dr = SQRT(SUM(rij(:)**2))
296 856401 : IF (dr > 1.e-6_dp) THEN
297 780112 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)/dr
298 780112 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)/dr
299 785371 : DO i = 2, nmat
300 5259 : gmcharge(iatom, i) = gmcharge(iatom, i) + rij(i - 1)*mcharge(jatom)/dr**3
301 785371 : gmcharge(jatom, i) = gmcharge(jatom, i) - rij(i - 1)*mcharge(iatom)/dr**3
302 : END DO
303 780112 : IF (use_virial) THEN
304 24 : DO i = 1, 3
305 24 : fij(i) = mcharge(iatom)*mcharge(jatom)*rij(i)/dr**3
306 : END DO
307 6 : CALL virial_pair_force(virial%pv_virial, 1.0_dp, fij, rij)
308 : END IF
309 : END IF
310 : END DO
311 : END DO
312 : END DO
313 : END IF
314 : END IF
315 :
316 : ! global sum of gamma*p arrays
317 : CALL get_qs_env(qs_env=qs_env, &
318 : atomic_kind_set=atomic_kind_set, &
319 31702 : force=force, para_env=para_env)
320 31702 : CALL para_env%sum(gmcharge(:, 1))
321 31702 : CALL para_env%sum(gchrg(:, :, 1))
322 :
323 31702 : IF (xtb_control%coulomb_lr) THEN
324 31702 : IF (do_ewald) THEN
325 : ! add self charge interaction and background charge contribution
326 118920 : gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:)
327 18240 : IF (ANY(periodic(:) == 1)) THEN
328 117360 : gmcharge(:, 1) = gmcharge(:, 1) - pi/alpha**2/deth
329 : END IF
330 : END IF
331 : END IF
332 :
333 : ! energy
334 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
335 : kind_of=kind_of, &
336 31702 : atom_of_kind=atom_of_kind)
337 31702 : ecsr = 0.0_dp
338 284588 : DO iatom = 1, natom
339 252886 : ikind = kind_of(iatom)
340 252886 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
341 252886 : CALL get_xtb_atom_param(xtb_kind, lmax=ni)
342 252886 : ni = ni + 1
343 679872 : ecsr = ecsr + SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, 1))
344 : END DO
345 :
346 31702 : energy%hartree = energy%hartree + 0.5_dp*ecsr
347 284588 : energy%hartree = energy%hartree + 0.5_dp*SUM(mcharge(:)*gmcharge(:, 1))
348 :
349 31702 : IF (atprop%energy) THEN
350 172 : CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
351 748 : DO ikind = 1, SIZE(local_particles%n_el)
352 576 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
353 576 : CALL get_xtb_atom_param(xtb_kind, lmax=ni, occupation=occ)
354 576 : ni = ni + 1
355 3456 : zeff = SUM(REAL(occ, KIND=dp))
356 4360 : DO ia = 1, local_particles%n_el(ikind)
357 3036 : iatom = local_particles%list(ikind)%array(ia)
358 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
359 7258 : 0.5_dp*SUM(REAL(occ(1:ni), KIND=dp)*gchrg(iatom, 1:ni, 1))
360 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
361 3612 : 0.5_dp*zeff*gmcharge(iatom, 1)
362 : END DO
363 : END DO
364 : END IF
365 :
366 31702 : IF (calculate_forces) THEN
367 4764 : DO iatom = 1, natom
368 4230 : ikind = kind_of(iatom)
369 4230 : atom_i = atom_of_kind(iatom)
370 4230 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
371 4230 : CALL get_xtb_atom_param(xtb_kind, lmax=ni)
372 : ! short range
373 4230 : ni = ni + 1
374 16920 : DO i = 1, 3
375 36972 : fij(i) = SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, i + 1))
376 : END DO
377 4230 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
378 4230 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
379 4230 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
380 : ! long range
381 16920 : DO i = 1, 3
382 16920 : fij(i) = gmcharge(iatom, i + 1)*mcharge(iatom)
383 : END DO
384 4230 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
385 4230 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
386 8994 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
387 : END DO
388 : END IF
389 :
390 31702 : IF (.NOT. just_energy) THEN
391 31214 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
392 31214 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
393 :
394 31214 : nimg = dft_control%nimages
395 31214 : NULLIFY (cell_to_index)
396 31214 : IF (nimg > 1) THEN
397 6532 : NULLIFY (kpoints)
398 6532 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
399 6532 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
400 : END IF
401 :
402 31214 : IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
403 752 : DO img = 1, nimg
404 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
405 752 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
406 : END DO
407 : END IF
408 :
409 31214 : NULLIFY (sap_int)
410 31214 : IF (do_gamma_stress) THEN
411 : ! derivative overlap integral (non collapsed)
412 116 : CALL xtb_dsint_list(qs_env, sap_int)
413 : END IF
414 :
415 31214 : IF (nimg == 1) THEN
416 : ! no k-points; all matrices have been transformed to periodic bsf
417 24682 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
418 1484672 : DO WHILE (dbcsr_iterator_blocks_left(iter))
419 1459990 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
420 1459990 : ikind = kind_of(irow)
421 1459990 : jkind = kind_of(icol)
422 :
423 : ! atomic parameters
424 1459990 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
425 1459990 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
426 1459990 : CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
427 1459990 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
428 :
429 1459990 : ni = SIZE(sblock, 1)
430 1459990 : nj = SIZE(sblock, 2)
431 5839960 : ALLOCATE (gcij(ni, nj))
432 5999347 : DO i = 1, ni
433 20142765 : DO j = 1, nj
434 14143418 : la = laoa(i) + 1
435 14143418 : lb = laob(j) + 1
436 18682775 : gcij(i, j) = 0.5_dp*(gchrg(irow, la, 1) + gchrg(icol, lb, 1))
437 : END DO
438 : END DO
439 1459990 : gmij = 0.5_dp*(gmcharge(irow, 1) + gmcharge(icol, 1))
440 2977557 : DO is = 1, SIZE(ks_matrix, 1)
441 1517567 : NULLIFY (ksblock)
442 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
443 1517567 : row=irow, col=icol, block=ksblock, found=found)
444 1517567 : CPASSERT(found)
445 40208083 : ksblock = ksblock - gcij*sblock
446 43185640 : ksblock = ksblock - gmij*sblock
447 : END DO
448 1459990 : IF (calculate_forces) THEN
449 46503 : atom_i = atom_of_kind(irow)
450 46503 : atom_j = atom_of_kind(icol)
451 46503 : NULLIFY (pblock)
452 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
453 46503 : row=irow, col=icol, block=pblock, found=found)
454 46503 : CPASSERT(found)
455 186012 : DO i = 1, 3
456 139509 : NULLIFY (dsblock)
457 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
458 139509 : row=irow, col=icol, block=dsblock, found=found)
459 139509 : CPASSERT(found)
460 139509 : fij(i) = 0.0_dp
461 : ! short range
462 1687116 : fi = -2.0_dp*SUM(pblock*dsblock*gcij)
463 139509 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
464 139509 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
465 : fij(i) = fij(i) + fi
466 : ! long range
467 1687116 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
468 139509 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
469 139509 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
470 325521 : fij(i) = fij(i) + fi
471 : END DO
472 : END IF
473 4404652 : DEALLOCATE (gcij)
474 : END DO
475 24682 : CALL dbcsr_iterator_stop(iter)
476 : ! stress tensor (needs recalculation of overlap integrals)
477 24682 : IF (do_gamma_stress) THEN
478 342 : DO ikind = 1, nkind
479 860 : DO jkind = 1, nkind
480 518 : iac = ikind + nkind*(jkind - 1)
481 518 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
482 : ! atomic parameters
483 372 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
484 372 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
485 372 : CALL get_xtb_atom_param(xtb_atom_a, lao=laoa, natorb=ni)
486 372 : CALL get_xtb_atom_param(xtb_atom_b, lao=laob, natorb=nj)
487 2040 : DO ia = 1, sap_int(iac)%nalist
488 1442 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE
489 1422 : iatom = sap_int(iac)%alist(ia)%aatom
490 130827 : DO ic = 1, sap_int(iac)%alist(ia)%nclist
491 128887 : jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
492 515548 : rij = sap_int(iac)%alist(ia)%clist(ic)%rac
493 515548 : dr = SQRT(SUM(rij(:)**2))
494 130329 : IF (dr > 1.e-6_dp) THEN
495 128153 : dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
496 512612 : ALLOCATE (gcij(ni, nj))
497 810853 : DO i = 1, ni
498 5444921 : DO j = 1, nj
499 4634068 : la = laoa(i) + 1
500 4634068 : lb = laob(j) + 1
501 5316768 : gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
502 : END DO
503 : END DO
504 128153 : gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
505 128153 : icol = MAX(iatom, jatom)
506 128153 : irow = MIN(iatom, jatom)
507 128153 : NULLIFY (pblock)
508 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
509 128153 : row=irow, col=icol, block=pblock, found=found)
510 128153 : CPASSERT(found)
511 128153 : fij = 0.0_dp
512 512612 : DO i = 1, 3
513 : ! short/long range
514 384459 : IF (irow == iatom) THEN
515 9632850 : f1 = -2.0_dp*SUM(pblock*dsint(:, :, i)*gcij)
516 9632850 : f2 = -2.0_dp*gmij*SUM(pblock*dsint(:, :, i))
517 : ELSE
518 6701619 : f1 = -2.0_dp*SUM(TRANSPOSE(pblock)*dsint(:, :, i)*gcij)
519 6701619 : f2 = -2.0_dp*gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i))
520 : END IF
521 512612 : fij(i) = f1 + f2
522 : END DO
523 128153 : DEALLOCATE (gcij)
524 128153 : fi = 1.0_dp
525 128153 : IF (iatom == jatom) fi = 0.5_dp
526 256306 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
527 : END IF
528 : END DO
529 : END DO
530 : END DO
531 : END DO
532 : END IF
533 : ELSE
534 6532 : NULLIFY (n_list)
535 6532 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
536 6532 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
537 2184819 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
538 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
539 2178287 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
540 :
541 2178287 : icol = MAX(iatom, jatom)
542 2178287 : irow = MIN(iatom, jatom)
543 :
544 2178287 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
545 2178287 : CPASSERT(ic > 0)
546 :
547 2178287 : NULLIFY (sblock)
548 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
549 2178287 : row=irow, col=icol, block=sblock, found=found)
550 2178287 : CPASSERT(found)
551 :
552 : ! atomic parameters
553 2178287 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
554 2178287 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
555 2178287 : CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
556 2178287 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
557 :
558 2178287 : ni = SIZE(sblock, 1)
559 2178287 : nj = SIZE(sblock, 2)
560 8713148 : ALLOCATE (gcij(ni, nj))
561 12490737 : DO i = 1, ni
562 77609905 : DO j = 1, nj
563 75431618 : IF (irow == iatom) THEN
564 37397384 : la = laoa(i) + 1
565 37397384 : lb = laob(j) + 1
566 37397384 : gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
567 : ELSE
568 27721784 : la = laoa(j) + 1
569 27721784 : lb = laob(i) + 1
570 27721784 : gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
571 : END IF
572 : END DO
573 : END DO
574 2178287 : gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
575 4705046 : DO is = 1, SIZE(ks_matrix, 1)
576 2526759 : NULLIFY (ksblock)
577 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
578 2526759 : row=irow, col=icol, block=ksblock, found=found)
579 2526759 : CPASSERT(found)
580 159381979 : ksblock = ksblock - gcij*sblock
581 164087025 : ksblock = ksblock - gmij*sblock
582 : END DO
583 :
584 2178287 : IF (calculate_forces) THEN
585 33164 : atom_i = atom_of_kind(iatom)
586 33164 : atom_j = atom_of_kind(jatom)
587 33164 : IF (irow /= iatom) THEN
588 13579 : gmij = -gmij
589 878373 : gcij = -gcij
590 : END IF
591 33164 : NULLIFY (pblock)
592 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
593 33164 : row=irow, col=icol, block=pblock, found=found)
594 33164 : CPASSERT(found)
595 132656 : DO i = 1, 3
596 99492 : NULLIFY (dsblock)
597 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
598 99492 : row=irow, col=icol, block=dsblock, found=found)
599 99492 : CPASSERT(found)
600 99492 : fij(i) = 0.0_dp
601 : ! short range
602 6655500 : fi = -2.0_dp*SUM(pblock*dsblock*gcij)
603 99492 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
604 99492 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
605 : fij(i) = fij(i) + fi
606 : ! long range
607 6655500 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
608 99492 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
609 99492 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
610 232148 : fij(i) = fij(i) + fi
611 : END DO
612 33164 : IF (use_virial) THEN
613 87068 : dr = SQRT(SUM(rij(:)**2))
614 21767 : IF (dr > 1.e-6_dp) THEN
615 21685 : fi = 1.0_dp
616 21685 : IF (iatom == jatom) fi = 0.5_dp
617 21685 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
618 : END IF
619 : END IF
620 : END IF
621 6541393 : DEALLOCATE (gcij)
622 :
623 : END DO
624 6532 : CALL neighbor_list_iterator_release(nl_iterator)
625 : END IF
626 :
627 31214 : IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
628 752 : DO img = 1, nimg
629 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
630 752 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
631 : END DO
632 : END IF
633 : END IF
634 :
635 31702 : IF (xtb_control%tb3_interaction) THEN
636 31702 : CALL get_qs_env(qs_env, nkind=nkind)
637 126808 : ALLOCATE (zeffk(nkind), xgamma(nkind))
638 112112 : DO ikind = 1, nkind
639 80410 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
640 112112 : CALL get_xtb_atom_param(xtb_kind, xgamma=xgamma(ikind), zeff=zeffk(ikind))
641 : END DO
642 : ! Diagonal 3rd order correction (DFTB3)
643 : CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, &
644 31702 : sap_int, calculate_forces, just_energy)
645 31702 : DEALLOCATE (zeffk, xgamma)
646 : END IF
647 :
648 31702 : IF (xtb_control%do_spinpol) THEN
649 1678 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
650 : CALL build_xtb_spinpol(qs_env, ks_matrix, matrix_p, energy, &
651 1678 : sap_int, calculate_forces, just_energy)
652 : END IF
653 :
654 : ! QMMM
655 31702 : IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
656 : CALL build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, &
657 862 : calculate_forces, just_energy)
658 : END IF
659 :
660 31702 : IF (do_gamma_stress) THEN
661 116 : CALL release_sap_int(sap_int)
662 : END IF
663 :
664 31702 : CALL timestop(handle)
665 :
666 63404 : END SUBROUTINE build_xtb_coulomb
667 :
668 : ! **************************************************************************************************
669 : !> \brief Computes the short-range gamma parameter from
670 : !> Nataga-Mishimoto-Ohno-Klopman formula for xTB
671 : !> WARNING: The xTB function (gamma - 1/r) has still an l-dependent longrange (1/r^3)
672 : !> behaviour. We use a cutoff function to smoothly remove this part.
673 : !> However, this will change energies and effect final results.
674 : !>
675 : !> \param gmat ...
676 : !> \param rab ...
677 : !> \param nla ...
678 : !> \param kappaa ...
679 : !> \param etaa ...
680 : !> \param nlb ...
681 : !> \param kappab ...
682 : !> \param etab ...
683 : !> \param kg ...
684 : !> \param rcut ...
685 : !> \par History
686 : !> 10.2018 JGH
687 : !> \version 1.1
688 : ! **************************************************************************************************
689 9387525 : SUBROUTINE gamma_rab_sr(gmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
690 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: gmat
691 : REAL(dp), INTENT(IN) :: rab
692 : INTEGER, INTENT(IN) :: nla
693 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappaa
694 : REAL(dp), INTENT(IN) :: etaa
695 : INTEGER, INTENT(IN) :: nlb
696 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappab
697 : REAL(dp), INTENT(IN) :: etab, kg, rcut
698 :
699 : REAL(KIND=dp), PARAMETER :: rsmooth = 1.0_dp
700 :
701 : INTEGER :: i, j
702 : REAL(KIND=dp) :: fcut, r, rk, x
703 9387525 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eta
704 :
705 37550100 : ALLOCATE (eta(nla, nlb))
706 9387525 : eta = 0.0_dp
707 :
708 24435788 : DO j = 1, nlb
709 51313721 : DO i = 1, nla
710 26877933 : eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
711 41926196 : eta(i, j) = 2._dp/eta(i, j)
712 : END DO
713 : END DO
714 :
715 51313721 : gmat = 0.0_dp
716 9387525 : IF (rab < 1.e-6_dp) THEN
717 : ! on site terms
718 686579 : gmat(:, :) = eta(:, :)
719 9258304 : ELSE IF (rab > rcut) THEN
720 : ! do nothing
721 : ELSE
722 9258304 : rk = rab**kg
723 50627142 : eta = eta**(-kg)
724 9258304 : IF (rab < rcut - rsmooth) THEN
725 : fcut = 1.0_dp
726 : ELSE
727 1039040 : r = rab - (rcut - rsmooth)
728 1039040 : x = r/rsmooth
729 1039040 : fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
730 : END IF
731 50627142 : gmat(:, :) = fcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg) - fcut/rab
732 : END IF
733 :
734 9387525 : DEALLOCATE (eta)
735 :
736 9387525 : END SUBROUTINE gamma_rab_sr
737 :
738 : ! **************************************************************************************************
739 : !> \brief Computes the derivative of the short-range gamma parameter from
740 : !> Nataga-Mishimoto-Ohno-Klopman formula for xTB
741 : !> WARNING: The xTB function (gamma - 1/r) has still an l-dependent longrange (1/r^3)
742 : !> behaviour. We use a cutoff function to smoothly remove this part.
743 : !> However, this will change energies and effect final results.
744 : !>
745 : !> \param dgmat ...
746 : !> \param rab ...
747 : !> \param nla ...
748 : !> \param kappaa ...
749 : !> \param etaa ...
750 : !> \param nlb ...
751 : !> \param kappab ...
752 : !> \param etab ...
753 : !> \param kg ...
754 : !> \param rcut ...
755 : !> \par History
756 : !> 10.2018 JGH
757 : !> \version 1.1
758 : ! **************************************************************************************************
759 385005 : SUBROUTINE dgamma_rab_sr(dgmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
760 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: dgmat
761 : REAL(dp), INTENT(IN) :: rab
762 : INTEGER, INTENT(IN) :: nla
763 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappaa
764 : REAL(dp), INTENT(IN) :: etaa
765 : INTEGER, INTENT(IN) :: nlb
766 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappab
767 : REAL(dp), INTENT(IN) :: etab, kg, rcut
768 :
769 : REAL(KIND=dp), PARAMETER :: rsmooth = 1.0_dp
770 :
771 : INTEGER :: i, j
772 : REAL(KIND=dp) :: dfcut, fcut, r, rk, x
773 385005 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eta
774 :
775 1540020 : ALLOCATE (eta(nla, nlb))
776 :
777 1042054 : DO j = 1, nlb
778 2344643 : DO i = 1, nla
779 1302589 : eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
780 1959638 : eta(i, j) = 2._dp/eta(i, j)
781 : END DO
782 : END DO
783 :
784 385005 : IF (rab < 1.e-6) THEN
785 : ! on site terms
786 0 : dgmat(:, :) = 0.0_dp
787 385005 : ELSE IF (rab > rcut) THEN
788 0 : dgmat(:, :) = 0.0_dp
789 : ELSE
790 2344643 : eta = eta**(-kg)
791 385005 : rk = rab**kg
792 385005 : IF (rab < rcut - rsmooth) THEN
793 : fcut = 1.0_dp
794 : dfcut = 0.0_dp
795 : ELSE
796 48570 : r = rab - (rcut - rsmooth)
797 48570 : x = r/rsmooth
798 48570 : fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
799 48570 : dfcut = -30._dp*x**4 + 60._dp*x**3 - 30._dp*x**2
800 48570 : dfcut = dfcut/rsmooth
801 : END IF
802 2344643 : dgmat(:, :) = dfcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg)
803 2344643 : dgmat(:, :) = dgmat(:, :) - dfcut/rab + fcut/rab**2
804 2344643 : dgmat(:, :) = dgmat(:, :) - fcut/(rk + eta(:, :))*(1._dp/(rk + eta(:, :)))**(1._dp/kg)*rk/rab
805 : END IF
806 :
807 385005 : DEALLOCATE (eta)
808 :
809 385005 : END SUBROUTINE dgamma_rab_sr
810 :
811 : ! **************************************************************************************************
812 : !> \brief ...
813 : !> \param qs_env ...
814 : !> \param sap_int ...
815 : ! **************************************************************************************************
816 122 : SUBROUTINE xtb_dsint_list(qs_env, sap_int)
817 :
818 : TYPE(qs_environment_type), POINTER :: qs_env
819 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
820 :
821 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xtb_dsint_list'
822 :
823 : INTEGER :: handle, i, iac, iatom, ikind, ilist, iset, jatom, jkind, jneighbor, jset, ldsab, &
824 : n1, n2, natorb_a, natorb_b, ncoa, ncob, nkind, nlist, nneighbor, nseta, nsetb, sgfa, sgfb
825 : INTEGER, DIMENSION(3) :: cell
826 122 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
827 122 : npgfb, nsgfa, nsgfb
828 122 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
829 : LOGICAL :: defined
830 : REAL(KIND=dp) :: dr
831 122 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
832 122 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
833 : REAL(KIND=dp), DIMENSION(3) :: rij
834 122 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
835 122 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb
836 : TYPE(clist_type), POINTER :: clist
837 : TYPE(dft_control_type), POINTER :: dft_control
838 122 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
839 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
840 : TYPE(neighbor_list_iterator_p_type), &
841 122 : DIMENSION(:), POINTER :: nl_iterator
842 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
843 122 : POINTER :: sab_orb
844 122 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
845 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
846 :
847 122 : CALL timeset(routineN, handle)
848 :
849 122 : CALL get_qs_env(qs_env=qs_env, nkind=nkind)
850 122 : CPASSERT(.NOT. ASSOCIATED(sap_int))
851 938 : ALLOCATE (sap_int(nkind*nkind))
852 694 : DO i = 1, nkind*nkind
853 572 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
854 694 : sap_int(i)%nalist = 0
855 : END DO
856 :
857 : CALL get_qs_env(qs_env=qs_env, &
858 : qs_kind_set=qs_kind_set, &
859 : dft_control=dft_control, &
860 122 : sab_orb=sab_orb)
861 :
862 : ! set up basis set lists
863 610 : ALLOCATE (basis_set_list(nkind))
864 122 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
865 :
866 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
867 122 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
868 130391 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
869 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, &
870 : jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, &
871 130269 : inode=jneighbor, cell=cell, r=rij)
872 130269 : iac = ikind + nkind*(jkind - 1)
873 : !
874 130269 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
875 130269 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
876 130269 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
877 130269 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
878 130269 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
879 130269 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
880 :
881 521076 : dr = SQRT(SUM(rij(:)**2))
882 :
883 : ! integral list
884 130269 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
885 399 : sap_int(iac)%a_kind = ikind
886 399 : sap_int(iac)%p_kind = jkind
887 399 : sap_int(iac)%nalist = nlist
888 2666 : ALLOCATE (sap_int(iac)%alist(nlist))
889 1868 : DO i = 1, nlist
890 1469 : NULLIFY (sap_int(iac)%alist(i)%clist)
891 1469 : sap_int(iac)%alist(i)%aatom = 0
892 1868 : sap_int(iac)%alist(i)%nclist = 0
893 : END DO
894 : END IF
895 130269 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
896 1449 : sap_int(iac)%alist(ilist)%aatom = iatom
897 1449 : sap_int(iac)%alist(ilist)%nclist = nneighbor
898 143310 : ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
899 131718 : DO i = 1, nneighbor
900 131718 : sap_int(iac)%alist(ilist)%clist(i)%catom = 0
901 : END DO
902 : END IF
903 130269 : clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
904 130269 : clist%catom = jatom
905 521076 : clist%cell = cell
906 521076 : clist%rac = rij
907 651345 : ALLOCATE (clist%acint(natorb_a, natorb_b, 3))
908 130269 : NULLIFY (clist%achint)
909 16577580 : clist%acint = 0._dp
910 130269 : clist%nsgf_cnt = 0
911 130269 : NULLIFY (clist%sgf_list)
912 :
913 : ! overlap
914 130269 : basis_set_a => basis_set_list(ikind)%gto_basis_set
915 130269 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
916 130269 : basis_set_b => basis_set_list(jkind)%gto_basis_set
917 130269 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
918 : ! basis ikind
919 130269 : first_sgfa => basis_set_a%first_sgf
920 130269 : la_max => basis_set_a%lmax
921 130269 : la_min => basis_set_a%lmin
922 130269 : npgfa => basis_set_a%npgf
923 130269 : nseta = basis_set_a%nset
924 130269 : nsgfa => basis_set_a%nsgf_set
925 130269 : rpgfa => basis_set_a%pgf_radius
926 130269 : set_radius_a => basis_set_a%set_radius
927 130269 : scon_a => basis_set_a%scon
928 130269 : zeta => basis_set_a%zet
929 : ! basis jkind
930 130269 : first_sgfb => basis_set_b%first_sgf
931 130269 : lb_max => basis_set_b%lmax
932 130269 : lb_min => basis_set_b%lmin
933 130269 : npgfb => basis_set_b%npgf
934 130269 : nsetb = basis_set_b%nset
935 130269 : nsgfb => basis_set_b%nsgf_set
936 130269 : rpgfb => basis_set_b%pgf_radius
937 130269 : set_radius_b => basis_set_b%set_radius
938 130269 : scon_b => basis_set_b%scon
939 130269 : zetb => basis_set_b%zet
940 :
941 130269 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
942 1042152 : ALLOCATE (oint(ldsab, ldsab, 4), owork(ldsab, ldsab))
943 651345 : ALLOCATE (sint(natorb_a, natorb_b, 4))
944 130269 : sint = 0.0_dp
945 :
946 443485 : DO iset = 1, nseta
947 313216 : ncoa = npgfa(iset)*ncoset(la_max(iset))
948 313216 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
949 313216 : sgfa = first_sgfa(1, iset)
950 1218787 : DO jset = 1, nsetb
951 775302 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
952 521296 : ncob = npgfb(jset)*ncoset(lb_max(jset))
953 521296 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
954 521296 : sgfb = first_sgfb(1, jset)
955 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
956 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
957 521296 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
958 : ! Contraction
959 2919696 : DO i = 1, 4
960 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
961 2085184 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
962 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), &
963 2860486 : sgfa, sgfb, trans=.FALSE.)
964 : END DO
965 : END DO
966 : END DO
967 : ! update dS/dR matrix
968 16577580 : clist%acint(1:natorb_a, 1:natorb_b, 1:3) = sint(1:natorb_a, 1:natorb_b, 2:4)
969 :
970 390929 : DEALLOCATE (oint, owork, sint)
971 :
972 : END DO
973 122 : CALL neighbor_list_iterator_release(nl_iterator)
974 :
975 122 : DEALLOCATE (basis_set_list)
976 :
977 122 : CALL timestop(handle)
978 :
979 244 : END SUBROUTINE xtb_dsint_list
980 :
981 18821697 : END MODULE xtb_coulomb
|
