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_types, ONLY: get_xtb_atom_param,&
77 : xtb_atom_type
78 : #include "./base/base_uses.f90"
79 :
80 : IMPLICIT NONE
81 :
82 : PRIVATE
83 :
84 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_coulomb'
85 :
86 : PUBLIC :: build_xtb_coulomb, gamma_rab_sr, dgamma_rab_sr, xtb_dsint_list
87 :
88 : CONTAINS
89 :
90 : ! **************************************************************************************************
91 : !> \brief ...
92 : !> \param qs_env ...
93 : !> \param ks_matrix ...
94 : !> \param rho ...
95 : !> \param charges ...
96 : !> \param mcharge ...
97 : !> \param energy ...
98 : !> \param calculate_forces ...
99 : !> \param just_energy ...
100 : ! **************************************************************************************************
101 29564 : SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
102 : calculate_forces, just_energy)
103 :
104 : TYPE(qs_environment_type), POINTER :: qs_env
105 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
106 : TYPE(qs_rho_type), POINTER :: rho
107 : REAL(dp), DIMENSION(:, :), INTENT(in) :: charges
108 : REAL(dp), DIMENSION(:), INTENT(in) :: mcharge
109 : TYPE(qs_energy_type), POINTER :: energy
110 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
111 :
112 : CHARACTER(len=*), PARAMETER :: routineN = 'build_xtb_coulomb'
113 :
114 : INTEGER :: atom_i, atom_j, ewald_type, handle, i, ia, iac, iatom, ic, icol, ikind, img, &
115 : irow, is, j, jatom, jkind, la, lb, lmaxa, lmaxb, natom, natorb_a, natorb_b, ni, nimg, nj, &
116 : nkind, nmat, za, zb
117 29564 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
118 : INTEGER, DIMENSION(25) :: laoa, laob
119 : INTEGER, DIMENSION(3) :: cellind, periodic
120 : INTEGER, DIMENSION(5) :: occ
121 29564 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
122 : LOGICAL :: defined, do_ewald, do_gamma_stress, &
123 : found, use_virial
124 : REAL(KIND=dp) :: alpha, deth, dr, ecsr, etaa, etab, f1, &
125 : f2, fi, gmij, kg, rcut, rcuta, rcutb, &
126 : zeff
127 29564 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xgamma, zeffk
128 29564 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gammab, gcij, gmcharge
129 29564 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gchrg
130 : REAL(KIND=dp), DIMENSION(25) :: gcint
131 : REAL(KIND=dp), DIMENSION(3) :: fij, rij
132 : REAL(KIND=dp), DIMENSION(5) :: kappaa, kappab
133 29564 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, ksblock, pblock, sblock
134 29564 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
135 29564 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
136 : TYPE(atprop_type), POINTER :: atprop
137 : TYPE(cell_type), POINTER :: cell
138 : TYPE(dbcsr_iterator_type) :: iter
139 29564 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
140 : TYPE(dft_control_type), POINTER :: dft_control
141 : TYPE(distribution_1d_type), POINTER :: local_particles
142 : TYPE(ewald_environment_type), POINTER :: ewald_env
143 : TYPE(ewald_pw_type), POINTER :: ewald_pw
144 : TYPE(kpoint_type), POINTER :: kpoints
145 : TYPE(mp_para_env_type), POINTER :: para_env
146 : TYPE(neighbor_list_iterator_p_type), &
147 29564 : DIMENSION(:), POINTER :: nl_iterator
148 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
149 29564 : POINTER :: n_list
150 29564 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
151 29564 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
152 29564 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
153 29564 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
154 : TYPE(virial_type), POINTER :: virial
155 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind
156 : TYPE(xtb_control_type), POINTER :: xtb_control
157 :
158 29564 : CALL timeset(routineN, handle)
159 :
160 29564 : NULLIFY (matrix_p, matrix_s, virial, atprop, dft_control)
161 :
162 : CALL get_qs_env(qs_env, &
163 : qs_kind_set=qs_kind_set, &
164 : particle_set=particle_set, &
165 : cell=cell, &
166 : virial=virial, &
167 : atprop=atprop, &
168 29564 : dft_control=dft_control)
169 :
170 29564 : xtb_control => dft_control%qs_control%xtb_control
171 :
172 29564 : use_virial = .FALSE.
173 29564 : IF (calculate_forces) THEN
174 896 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
175 : END IF
176 :
177 29564 : do_gamma_stress = .FALSE.
178 29564 : IF (.NOT. just_energy .AND. use_virial) THEN
179 132 : IF (dft_control%nimages == 1) do_gamma_stress = .TRUE.
180 : END IF
181 :
182 29564 : IF (atprop%energy) THEN
183 172 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
184 172 : natom = SIZE(particle_set)
185 172 : CALL atprop_array_init(atprop%atecoul, natom)
186 : END IF
187 :
188 29564 : IF (calculate_forces) THEN
189 : nmat = 4
190 : ELSE
191 29050 : nmat = 1
192 : END IF
193 :
194 29564 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
195 118256 : ALLOCATE (gchrg(natom, 5, nmat))
196 29564 : gchrg = 0._dp
197 118256 : ALLOCATE (gmcharge(natom, nmat))
198 29564 : gmcharge = 0._dp
199 :
200 : ! short range contribution (gamma)
201 : ! loop over all atom pairs (sab_xtbe)
202 29564 : kg = xtb_control%kg
203 29564 : NULLIFY (n_list)
204 29564 : CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
205 29564 : IF (.NOT. ASSOCIATED(n_list)) THEN
206 0 : CPABORT("sab_xtbe neighbor list is not associated in build_xtb_coulomb")
207 : END IF
208 29564 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
209 8810894 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
210 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
211 8781330 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
212 8781330 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
213 8781330 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
214 8781330 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
215 8781323 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
216 8781323 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
217 8781323 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
218 : ! atomic parameters
219 8781316 : CALL get_xtb_atom_param(xtb_atom_a, eta=etaa, lmax=lmaxa, kappa=kappaa, rcut=rcuta)
220 8781316 : CALL get_xtb_atom_param(xtb_atom_b, eta=etab, lmax=lmaxb, kappa=kappab, rcut=rcutb)
221 : ! gamma matrix
222 8781316 : ni = lmaxa + 1
223 8781316 : nj = lmaxb + 1
224 35125264 : ALLOCATE (gammab(ni, nj))
225 8781316 : rcut = rcuta + rcutb
226 35125264 : dr = SQRT(SUM(rij(:)**2))
227 8781316 : CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
228 106951727 : gchrg(iatom, 1:ni, 1) = gchrg(iatom, 1:ni, 1) + MATMUL(gammab, charges(jatom, 1:nj))
229 8781316 : IF (iatom /= jatom) THEN
230 86812743 : gchrg(jatom, 1:nj, 1) = gchrg(jatom, 1:nj, 1) + MATMUL(charges(iatom, 1:ni), gammab)
231 : END IF
232 8781316 : IF (calculate_forces) THEN
233 359783 : IF (dr > 1.e-6_dp) THEN
234 357550 : CALL dgamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
235 1430200 : DO i = 1, 3
236 : gchrg(iatom, 1:ni, i + 1) = gchrg(iatom, 1:ni, i + 1) &
237 14482533 : + MATMUL(gammab, charges(jatom, 1:nj))*rij(i)/dr
238 1430200 : IF (iatom /= jatom) THEN
239 : gchrg(jatom, 1:nj, i + 1) = gchrg(jatom, 1:nj, i + 1) &
240 11553132 : - MATMUL(charges(iatom, 1:ni), gammab)*rij(i)/dr
241 : END IF
242 : END DO
243 357550 : IF (use_virial) THEN
244 2839192 : gcint(1:ni) = MATMUL(gammab, charges(jatom, 1:nj))
245 830300 : DO i = 1, 3
246 2046743 : fij(i) = -SUM(charges(iatom, 1:ni)*gcint(1:ni))*rij(i)/dr
247 : END DO
248 207575 : fi = 1.0_dp
249 207575 : IF (iatom == jatom) fi = 0.5_dp
250 207575 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
251 : END IF
252 : END IF
253 : END IF
254 35125285 : DEALLOCATE (gammab)
255 : END DO
256 29564 : CALL neighbor_list_iterator_release(nl_iterator)
257 :
258 : ! 1/R contribution
259 :
260 29564 : IF (xtb_control%coulomb_lr) THEN
261 29564 : do_ewald = xtb_control%do_ewald
262 29564 : IF (do_ewald) THEN
263 : ! Ewald sum
264 14314 : NULLIFY (ewald_env, ewald_pw)
265 : CALL get_qs_env(qs_env=qs_env, &
266 14314 : ewald_env=ewald_env, ewald_pw=ewald_pw)
267 14314 : CALL get_cell(cell=cell, periodic=periodic, deth=deth)
268 14314 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
269 14314 : CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
270 14314 : CALL tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial)
271 0 : SELECT CASE (ewald_type)
272 : CASE DEFAULT
273 0 : CPABORT("Invalid Ewald type")
274 : CASE (do_ewald_none)
275 0 : CPABORT("Not allowed with xTB/DFTB")
276 : CASE (do_ewald_ewald)
277 0 : CPABORT("Standard Ewald not implemented in xTB/DFTB")
278 : CASE (do_ewald_pme)
279 0 : CPABORT("PME not implemented in xTB/DFTB")
280 : CASE (do_ewald_spme)
281 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
282 14314 : gmcharge, mcharge, calculate_forces, virial, use_virial)
283 : END SELECT
284 : ELSE
285 : ! direct sum
286 : CALL get_qs_env(qs_env=qs_env, &
287 15250 : local_particles=local_particles)
288 56600 : DO ikind = 1, SIZE(local_particles%n_el)
289 127730 : DO ia = 1, local_particles%n_el(ikind)
290 71130 : iatom = local_particles%list(ikind)%array(ia)
291 846767 : DO jatom = 1, iatom - 1
292 2937148 : rij = particle_set(iatom)%r - particle_set(jatom)%r
293 2937148 : rij = pbc(rij, cell)
294 2937148 : dr = SQRT(SUM(rij(:)**2))
295 805417 : IF (dr > 1.e-6_dp) THEN
296 734287 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)/dr
297 734287 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)/dr
298 739519 : DO i = 2, nmat
299 5232 : gmcharge(iatom, i) = gmcharge(iatom, i) + rij(i - 1)*mcharge(jatom)/dr**3
300 739519 : gmcharge(jatom, i) = gmcharge(jatom, i) - rij(i - 1)*mcharge(iatom)/dr**3
301 : END DO
302 734287 : IF (use_virial) THEN
303 24 : DO i = 1, 3
304 24 : fij(i) = mcharge(iatom)*mcharge(jatom)*rij(i)/dr**3
305 : END DO
306 6 : CALL virial_pair_force(virial%pv_virial, 1.0_dp, fij, rij)
307 : END IF
308 : END IF
309 : END DO
310 : END DO
311 : END DO
312 : END IF
313 : END IF
314 :
315 : ! global sum of gamma*p arrays
316 : CALL get_qs_env(qs_env=qs_env, &
317 : atomic_kind_set=atomic_kind_set, &
318 29564 : force=force, para_env=para_env)
319 29564 : CALL para_env%sum(gmcharge(:, 1))
320 29564 : CALL para_env%sum(gchrg(:, :, 1))
321 :
322 29564 : IF (xtb_control%coulomb_lr) THEN
323 29564 : IF (do_ewald) THEN
324 : ! add self charge interaction and background charge contribution
325 112570 : gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:)
326 16932 : IF (ANY(periodic(:) == 1)) THEN
327 111010 : gmcharge(:, 1) = gmcharge(:, 1) - pi/alpha**2/deth
328 : END IF
329 : END IF
330 : END IF
331 :
332 : ! energy
333 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
334 : kind_of=kind_of, &
335 29564 : atom_of_kind=atom_of_kind)
336 29564 : ecsr = 0.0_dp
337 267122 : DO iatom = 1, natom
338 237558 : ikind = kind_of(iatom)
339 237558 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
340 237558 : CALL get_xtb_atom_param(xtb_kind, lmax=ni)
341 237558 : ni = ni + 1
342 639824 : ecsr = ecsr + SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, 1))
343 : END DO
344 :
345 29564 : energy%hartree = energy%hartree + 0.5_dp*ecsr
346 267122 : energy%hartree = energy%hartree + 0.5_dp*SUM(mcharge(:)*gmcharge(:, 1))
347 :
348 29564 : IF (atprop%energy) THEN
349 172 : CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
350 748 : DO ikind = 1, SIZE(local_particles%n_el)
351 576 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
352 576 : CALL get_xtb_atom_param(xtb_kind, lmax=ni, occupation=occ)
353 576 : ni = ni + 1
354 3456 : zeff = SUM(REAL(occ, KIND=dp))
355 4360 : DO ia = 1, local_particles%n_el(ikind)
356 3036 : iatom = local_particles%list(ikind)%array(ia)
357 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
358 7258 : 0.5_dp*SUM(REAL(occ(1:ni), KIND=dp)*gchrg(iatom, 1:ni, 1))
359 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
360 3612 : 0.5_dp*zeff*gmcharge(iatom, 1)
361 : END DO
362 : END DO
363 : END IF
364 :
365 29564 : IF (calculate_forces) THEN
366 4674 : DO iatom = 1, natom
367 4160 : ikind = kind_of(iatom)
368 4160 : atom_i = atom_of_kind(iatom)
369 4160 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
370 4160 : CALL get_xtb_atom_param(xtb_kind, lmax=ni)
371 : ! short range
372 4160 : ni = ni + 1
373 16640 : DO i = 1, 3
374 36392 : fij(i) = SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, i + 1))
375 : END DO
376 4160 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
377 4160 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
378 4160 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
379 : ! long range
380 16640 : DO i = 1, 3
381 16640 : fij(i) = gmcharge(iatom, i + 1)*mcharge(iatom)
382 : END DO
383 4160 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
384 4160 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
385 8834 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
386 : END DO
387 : END IF
388 :
389 29564 : IF (.NOT. just_energy) THEN
390 29076 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
391 29076 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
392 :
393 29076 : nimg = dft_control%nimages
394 29076 : NULLIFY (cell_to_index)
395 29076 : IF (nimg > 1) THEN
396 5876 : NULLIFY (kpoints)
397 5876 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
398 5876 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
399 : END IF
400 :
401 29076 : IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
402 440 : DO img = 1, nimg
403 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
404 440 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
405 : END DO
406 : END IF
407 :
408 29076 : NULLIFY (sap_int)
409 29076 : IF (do_gamma_stress) THEN
410 : ! derivative overlap integral (non collapsed)
411 112 : CALL xtb_dsint_list(qs_env, sap_int)
412 : END IF
413 :
414 29076 : IF (nimg == 1) THEN
415 : ! no k-points; all matrices have been transformed to periodic bsf
416 23200 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
417 1428973 : DO WHILE (dbcsr_iterator_blocks_left(iter))
418 1405773 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
419 1405773 : ikind = kind_of(irow)
420 1405773 : jkind = kind_of(icol)
421 :
422 : ! atomic parameters
423 1405773 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
424 1405773 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
425 1405773 : CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
426 1405773 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
427 :
428 1405773 : ni = SIZE(sblock, 1)
429 1405773 : nj = SIZE(sblock, 2)
430 5623092 : ALLOCATE (gcij(ni, nj))
431 5745790 : DO i = 1, ni
432 19321154 : DO j = 1, nj
433 13575364 : la = laoa(i) + 1
434 13575364 : lb = laob(j) + 1
435 17915381 : gcij(i, j) = 0.5_dp*(gchrg(irow, la, 1) + gchrg(icol, lb, 1))
436 : END DO
437 : END DO
438 1405773 : gmij = 0.5_dp*(gmcharge(irow, 1) + gmcharge(icol, 1))
439 2823342 : DO is = 1, SIZE(ks_matrix, 1)
440 1417569 : NULLIFY (ksblock)
441 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
442 1417569 : row=irow, col=icol, block=ksblock, found=found)
443 1417569 : CPASSERT(found)
444 37433361 : ksblock = ksblock - gcij*sblock
445 40256703 : ksblock = ksblock - gmij*sblock
446 : END DO
447 1405773 : IF (calculate_forces) THEN
448 46443 : atom_i = atom_of_kind(irow)
449 46443 : atom_j = atom_of_kind(icol)
450 46443 : NULLIFY (pblock)
451 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
452 46443 : row=irow, col=icol, block=pblock, found=found)
453 46443 : CPASSERT(found)
454 185772 : DO i = 1, 3
455 139329 : NULLIFY (dsblock)
456 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
457 139329 : row=irow, col=icol, block=dsblock, found=found)
458 139329 : CPASSERT(found)
459 139329 : fij(i) = 0.0_dp
460 : ! short range
461 1685004 : fi = -2.0_dp*SUM(pblock*dsblock*gcij)
462 139329 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
463 139329 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
464 : fij(i) = fij(i) + fi
465 : ! long range
466 1685004 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
467 139329 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
468 139329 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
469 325101 : fij(i) = fij(i) + fi
470 : END DO
471 : END IF
472 4240519 : DEALLOCATE (gcij)
473 : END DO
474 23200 : CALL dbcsr_iterator_stop(iter)
475 : ! stress tensor (needs recalculation of overlap integrals)
476 23200 : IF (do_gamma_stress) THEN
477 326 : DO ikind = 1, nkind
478 808 : DO jkind = 1, nkind
479 482 : iac = ikind + nkind*(jkind - 1)
480 482 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
481 : ! atomic parameters
482 354 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
483 354 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
484 354 : CALL get_xtb_atom_param(xtb_atom_a, lao=laoa, natorb=ni)
485 354 : CALL get_xtb_atom_param(xtb_atom_b, lao=laob, natorb=nj)
486 1992 : DO ia = 1, sap_int(iac)%nalist
487 1424 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE
488 1404 : iatom = sap_int(iac)%alist(ia)%aatom
489 130562 : DO ic = 1, sap_int(iac)%alist(ia)%nclist
490 128676 : jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
491 514704 : rij = sap_int(iac)%alist(ia)%clist(ic)%rac
492 514704 : dr = SQRT(SUM(rij(:)**2))
493 130100 : IF (dr > 1.e-6_dp) THEN
494 127950 : dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
495 511800 : ALLOCATE (gcij(ni, nj))
496 810078 : DO i = 1, ni
497 5442510 : DO j = 1, nj
498 4632432 : la = laoa(i) + 1
499 4632432 : lb = laob(j) + 1
500 5314560 : gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
501 : END DO
502 : END DO
503 127950 : gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
504 127950 : icol = MAX(iatom, jatom)
505 127950 : irow = MIN(iatom, jatom)
506 127950 : NULLIFY (pblock)
507 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
508 127950 : row=irow, col=icol, block=pblock, found=found)
509 127950 : CPASSERT(found)
510 127950 : fij = 0.0_dp
511 511800 : DO i = 1, 3
512 : ! short/long range
513 383850 : IF (irow == iatom) THEN
514 9628827 : f1 = -2.0_dp*SUM(pblock*dsint(:, :, i)*gcij)
515 9628827 : f2 = -2.0_dp*gmij*SUM(pblock*dsint(:, :, i))
516 : ELSE
517 6698409 : f1 = -2.0_dp*SUM(TRANSPOSE(pblock)*dsint(:, :, i)*gcij)
518 6698409 : f2 = -2.0_dp*gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i))
519 : END IF
520 511800 : fij(i) = f1 + f2
521 : END DO
522 127950 : DEALLOCATE (gcij)
523 127950 : fi = 1.0_dp
524 127950 : IF (iatom == jatom) fi = 0.5_dp
525 255900 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
526 : END IF
527 : END DO
528 : END DO
529 : END DO
530 : END DO
531 : END IF
532 : ELSE
533 5876 : NULLIFY (n_list)
534 5876 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
535 5876 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
536 2134139 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
537 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
538 2128263 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
539 :
540 2128263 : icol = MAX(iatom, jatom)
541 2128263 : irow = MIN(iatom, jatom)
542 :
543 2128263 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
544 2128263 : CPASSERT(ic > 0)
545 :
546 2128263 : NULLIFY (sblock)
547 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
548 2128263 : row=irow, col=icol, block=sblock, found=found)
549 2128263 : CPASSERT(found)
550 :
551 : ! atomic parameters
552 2128263 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
553 2128263 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
554 2128263 : CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
555 2128263 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
556 :
557 2128263 : ni = SIZE(sblock, 1)
558 2128263 : nj = SIZE(sblock, 2)
559 8513052 : ALLOCATE (gcij(ni, nj))
560 12831593 : DO i = 1, ni
561 84376021 : DO j = 1, nj
562 82247758 : IF (irow == iatom) THEN
563 40922826 : la = laoa(i) + 1
564 40922826 : lb = laob(j) + 1
565 40922826 : gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
566 : ELSE
567 30621602 : la = laoa(j) + 1
568 30621602 : lb = laob(i) + 1
569 30621602 : gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
570 : END IF
571 : END DO
572 : END DO
573 2128263 : gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
574 4572930 : DO is = 1, SIZE(ks_matrix, 1)
575 2444667 : NULLIFY (ksblock)
576 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
577 2444667 : row=irow, col=icol, block=ksblock, found=found)
578 2444667 : CPASSERT(found)
579 190339231 : ksblock = ksblock - gcij*sblock
580 194912161 : ksblock = ksblock - gmij*sblock
581 : END DO
582 :
583 2128263 : IF (calculate_forces) THEN
584 32376 : atom_i = atom_of_kind(iatom)
585 32376 : atom_j = atom_of_kind(jatom)
586 32376 : IF (irow /= iatom) THEN
587 13261 : gmij = -gmij
588 874875 : gcij = -gcij
589 : END IF
590 32376 : NULLIFY (pblock)
591 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
592 32376 : row=irow, col=icol, block=pblock, found=found)
593 32376 : CPASSERT(found)
594 129504 : DO i = 1, 3
595 97128 : NULLIFY (dsblock)
596 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
597 97128 : row=irow, col=icol, block=dsblock, found=found)
598 97128 : CPASSERT(found)
599 97128 : fij(i) = 0.0_dp
600 : ! short range
601 6629472 : fi = -2.0_dp*SUM(pblock*dsblock*gcij)
602 97128 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
603 97128 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
604 : fij(i) = fij(i) + fi
605 : ! long range
606 6629472 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
607 97128 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
608 97128 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
609 226632 : fij(i) = fij(i) + fi
610 : END DO
611 32376 : IF (use_virial) THEN
612 85492 : dr = SQRT(SUM(rij(:)**2))
613 21373 : IF (dr > 1.e-6_dp) THEN
614 21295 : fi = 1.0_dp
615 21295 : IF (iatom == jatom) fi = 0.5_dp
616 21295 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
617 : END IF
618 : END IF
619 : END IF
620 6390665 : DEALLOCATE (gcij)
621 :
622 : END DO
623 5876 : CALL neighbor_list_iterator_release(nl_iterator)
624 : END IF
625 :
626 29076 : IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
627 440 : DO img = 1, nimg
628 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
629 440 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
630 : END DO
631 : END IF
632 : END IF
633 :
634 29564 : IF (xtb_control%tb3_interaction) THEN
635 29564 : CALL get_qs_env(qs_env, nkind=nkind)
636 118256 : ALLOCATE (zeffk(nkind), xgamma(nkind))
637 103654 : DO ikind = 1, nkind
638 74090 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
639 103654 : CALL get_xtb_atom_param(xtb_kind, xgamma=xgamma(ikind), zeff=zeffk(ikind))
640 : END DO
641 : ! Diagonal 3rd order correction (DFTB3)
642 : CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, &
643 29564 : sap_int, calculate_forces, just_energy)
644 29564 : DEALLOCATE (zeffk, xgamma)
645 : END IF
646 :
647 : ! QMMM
648 29564 : IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
649 : CALL build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, &
650 862 : calculate_forces, just_energy)
651 : END IF
652 :
653 29564 : IF (do_gamma_stress) THEN
654 112 : CALL release_sap_int(sap_int)
655 : END IF
656 :
657 29564 : CALL timestop(handle)
658 :
659 59128 : END SUBROUTINE build_xtb_coulomb
660 :
661 : ! **************************************************************************************************
662 : !> \brief Computes the short-range gamma parameter from
663 : !> Nataga-Mishimoto-Ohno-Klopman formula for xTB
664 : !> WARNING: The xTB function (gamma - 1/r) has still an l-dependent longrange (1/r^3)
665 : !> behaviour. We use a cutoff function to smoothly remove this part.
666 : !> However, this will change energies and effect final results.
667 : !>
668 : !> \param gmat ...
669 : !> \param rab ...
670 : !> \param nla ...
671 : !> \param kappaa ...
672 : !> \param etaa ...
673 : !> \param nlb ...
674 : !> \param kappab ...
675 : !> \param etab ...
676 : !> \param kg ...
677 : !> \param rcut ...
678 : !> \par History
679 : !> 10.2018 JGH
680 : !> \version 1.1
681 : ! **************************************************************************************************
682 9003843 : SUBROUTINE gamma_rab_sr(gmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
683 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: gmat
684 : REAL(dp), INTENT(IN) :: rab
685 : INTEGER, INTENT(IN) :: nla
686 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappaa
687 : REAL(dp), INTENT(IN) :: etaa
688 : INTEGER, INTENT(IN) :: nlb
689 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappab
690 : REAL(dp), INTENT(IN) :: etab, kg, rcut
691 :
692 : REAL(KIND=dp), PARAMETER :: rsmooth = 1.0_dp
693 :
694 : INTEGER :: i, j
695 : REAL(KIND=dp) :: fcut, r, rk, x
696 9003843 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eta
697 :
698 36015372 : ALLOCATE (eta(nla, nlb))
699 9003843 : eta = 0.0_dp
700 :
701 23567154 : DO j = 1, nlb
702 49956633 : DO i = 1, nla
703 26389479 : eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
704 40952790 : eta(i, j) = 2._dp/eta(i, j)
705 : END DO
706 : END DO
707 :
708 49956633 : gmat = 0.0_dp
709 9003843 : IF (rab < 1.e-6_dp) THEN
710 : ! on site terms
711 648275 : gmat(:, :) = eta(:, :)
712 8882498 : ELSEIF (rab > rcut) THEN
713 : ! do nothing
714 : ELSE
715 8882498 : rk = rab**kg
716 49308358 : eta = eta**(-kg)
717 8882498 : IF (rab < rcut - rsmooth) THEN
718 : fcut = 1.0_dp
719 : ELSE
720 1004419 : r = rab - (rcut - rsmooth)
721 1004419 : x = r/rsmooth
722 1004419 : fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
723 : END IF
724 49308358 : gmat(:, :) = fcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg) - fcut/rab
725 : END IF
726 :
727 9003843 : DEALLOCATE (eta)
728 :
729 9003843 : END SUBROUTINE gamma_rab_sr
730 :
731 : ! **************************************************************************************************
732 : !> \brief Computes the derivative of the short-range gamma parameter from
733 : !> Nataga-Mishimoto-Ohno-Klopman formula for xTB
734 : !> WARNING: The xTB function (gamma - 1/r) has still an l-dependent longrange (1/r^3)
735 : !> behaviour. We use a cutoff function to smoothly remove this part.
736 : !> However, this will change energies and effect final results.
737 : !>
738 : !> \param dgmat ...
739 : !> \param rab ...
740 : !> \param nla ...
741 : !> \param kappaa ...
742 : !> \param etaa ...
743 : !> \param nlb ...
744 : !> \param kappab ...
745 : !> \param etab ...
746 : !> \param kg ...
747 : !> \param rcut ...
748 : !> \par History
749 : !> 10.2018 JGH
750 : !> \version 1.1
751 : ! **************************************************************************************************
752 382127 : SUBROUTINE dgamma_rab_sr(dgmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
753 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: dgmat
754 : REAL(dp), INTENT(IN) :: rab
755 : INTEGER, INTENT(IN) :: nla
756 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappaa
757 : REAL(dp), INTENT(IN) :: etaa
758 : INTEGER, INTENT(IN) :: nlb
759 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappab
760 : REAL(dp), INTENT(IN) :: etab, kg, rcut
761 :
762 : REAL(KIND=dp), PARAMETER :: rsmooth = 1.0_dp
763 :
764 : INTEGER :: i, j
765 : REAL(KIND=dp) :: dfcut, fcut, r, rk, x
766 382127 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eta
767 :
768 1528508 : ALLOCATE (eta(nla, nlb))
769 :
770 1035041 : DO j = 1, nlb
771 2331567 : DO i = 1, nla
772 1296526 : eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
773 1949440 : eta(i, j) = 2._dp/eta(i, j)
774 : END DO
775 : END DO
776 :
777 382127 : IF (rab < 1.e-6) THEN
778 : ! on site terms
779 0 : dgmat(:, :) = 0.0_dp
780 382127 : ELSEIF (rab > rcut) THEN
781 0 : dgmat(:, :) = 0.0_dp
782 : ELSE
783 2331567 : eta = eta**(-kg)
784 382127 : rk = rab**kg
785 382127 : IF (rab < rcut - rsmooth) THEN
786 : fcut = 1.0_dp
787 : dfcut = 0.0_dp
788 : ELSE
789 48311 : r = rab - (rcut - rsmooth)
790 48311 : x = r/rsmooth
791 48311 : fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
792 48311 : dfcut = -30._dp*x**4 + 60._dp*x**3 - 30._dp*x**2
793 48311 : dfcut = dfcut/rsmooth
794 : END IF
795 2331567 : dgmat(:, :) = dfcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg)
796 2331567 : dgmat(:, :) = dgmat(:, :) - dfcut/rab + fcut/rab**2
797 2331567 : dgmat(:, :) = dgmat(:, :) - fcut/(rk + eta(:, :))*(1._dp/(rk + eta(:, :)))**(1._dp/kg)*rk/rab
798 : END IF
799 :
800 382127 : DEALLOCATE (eta)
801 :
802 382127 : END SUBROUTINE dgamma_rab_sr
803 :
804 : ! **************************************************************************************************
805 : !> \brief ...
806 : !> \param qs_env ...
807 : !> \param sap_int ...
808 : ! **************************************************************************************************
809 118 : SUBROUTINE xtb_dsint_list(qs_env, sap_int)
810 :
811 : TYPE(qs_environment_type), POINTER :: qs_env
812 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
813 :
814 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xtb_dsint_list'
815 :
816 : INTEGER :: handle, i, iac, iatom, ikind, ilist, iset, jatom, jkind, jneighbor, jset, ldsab, &
817 : n1, n2, natorb_a, natorb_b, ncoa, ncob, nkind, nlist, nneighbor, nseta, nsetb, sgfa, sgfb
818 : INTEGER, DIMENSION(3) :: cell
819 118 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
820 118 : npgfb, nsgfa, nsgfb
821 118 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
822 : LOGICAL :: defined
823 : REAL(KIND=dp) :: dr
824 118 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
825 118 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
826 : REAL(KIND=dp), DIMENSION(3) :: rij
827 118 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
828 118 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb
829 : TYPE(clist_type), POINTER :: clist
830 : TYPE(dft_control_type), POINTER :: dft_control
831 118 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
832 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
833 : TYPE(neighbor_list_iterator_p_type), &
834 118 : DIMENSION(:), POINTER :: nl_iterator
835 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
836 118 : POINTER :: sab_orb
837 118 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
838 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
839 :
840 118 : CALL timeset(routineN, handle)
841 :
842 118 : CALL get_qs_env(qs_env=qs_env, nkind=nkind)
843 118 : CPASSERT(.NOT. ASSOCIATED(sap_int))
844 890 : ALLOCATE (sap_int(nkind*nkind))
845 654 : DO i = 1, nkind*nkind
846 536 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
847 654 : sap_int(i)%nalist = 0
848 : END DO
849 :
850 : CALL get_qs_env(qs_env=qs_env, &
851 : qs_kind_set=qs_kind_set, &
852 : dft_control=dft_control, &
853 118 : sab_orb=sab_orb)
854 :
855 : ! set up basis set lists
856 586 : ALLOCATE (basis_set_list(nkind))
857 118 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
858 :
859 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
860 118 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
861 130176 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
862 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, &
863 : jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, &
864 130058 : inode=jneighbor, cell=cell, r=rij)
865 130058 : iac = ikind + nkind*(jkind - 1)
866 : !
867 130058 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
868 130058 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
869 130058 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
870 130058 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
871 130058 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
872 130058 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
873 :
874 520232 : dr = SQRT(SUM(rij(:)**2))
875 :
876 : ! integral list
877 130058 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
878 381 : sap_int(iac)%a_kind = ikind
879 381 : sap_int(iac)%p_kind = jkind
880 381 : sap_int(iac)%nalist = nlist
881 2594 : ALLOCATE (sap_int(iac)%alist(nlist))
882 1832 : DO i = 1, nlist
883 1451 : NULLIFY (sap_int(iac)%alist(i)%clist)
884 1451 : sap_int(iac)%alist(i)%aatom = 0
885 1832 : sap_int(iac)%alist(i)%nclist = 0
886 : END DO
887 : END IF
888 130058 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
889 1431 : sap_int(iac)%alist(ilist)%aatom = iatom
890 1431 : sap_int(iac)%alist(ilist)%nclist = nneighbor
891 142937 : ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
892 131489 : DO i = 1, nneighbor
893 131489 : sap_int(iac)%alist(ilist)%clist(i)%catom = 0
894 : END DO
895 : END IF
896 130058 : clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
897 130058 : clist%catom = jatom
898 520232 : clist%cell = cell
899 520232 : clist%rac = rij
900 650290 : ALLOCATE (clist%acint(natorb_a, natorb_b, 3))
901 130058 : NULLIFY (clist%achint)
902 16569800 : clist%acint = 0._dp
903 130058 : clist%nsgf_cnt = 0
904 130058 : NULLIFY (clist%sgf_list)
905 :
906 : ! overlap
907 130058 : basis_set_a => basis_set_list(ikind)%gto_basis_set
908 130058 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
909 130058 : basis_set_b => basis_set_list(jkind)%gto_basis_set
910 130058 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
911 : ! basis ikind
912 130058 : first_sgfa => basis_set_a%first_sgf
913 130058 : la_max => basis_set_a%lmax
914 130058 : la_min => basis_set_a%lmin
915 130058 : npgfa => basis_set_a%npgf
916 130058 : nseta = basis_set_a%nset
917 130058 : nsgfa => basis_set_a%nsgf_set
918 130058 : rpgfa => basis_set_a%pgf_radius
919 130058 : set_radius_a => basis_set_a%set_radius
920 130058 : scon_a => basis_set_a%scon
921 130058 : zeta => basis_set_a%zet
922 : ! basis jkind
923 130058 : first_sgfb => basis_set_b%first_sgf
924 130058 : lb_max => basis_set_b%lmax
925 130058 : lb_min => basis_set_b%lmin
926 130058 : npgfb => basis_set_b%npgf
927 130058 : nsetb = basis_set_b%nset
928 130058 : nsgfb => basis_set_b%nsgf_set
929 130058 : rpgfb => basis_set_b%pgf_radius
930 130058 : set_radius_b => basis_set_b%set_radius
931 130058 : scon_b => basis_set_b%scon
932 130058 : zetb => basis_set_b%zet
933 :
934 130058 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
935 1040464 : ALLOCATE (oint(ldsab, ldsab, 4), owork(ldsab, ldsab))
936 650290 : ALLOCATE (sint(natorb_a, natorb_b, 4))
937 130058 : sint = 0.0_dp
938 :
939 442852 : DO iset = 1, nseta
940 312794 : ncoa = npgfa(iset)*ncoset(la_max(iset))
941 312794 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
942 312794 : sgfa = first_sgfa(1, iset)
943 1217310 : DO jset = 1, nsetb
944 774458 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
945 520512 : ncob = npgfb(jset)*ncoset(lb_max(jset))
946 520512 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
947 520512 : sgfb = first_sgfb(1, jset)
948 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
949 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
950 520512 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
951 : ! Contraction
952 2915354 : DO i = 1, 4
953 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
954 2082048 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
955 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), &
956 2856506 : sgfa, sgfb, trans=.FALSE.)
957 : END DO
958 : END DO
959 : END DO
960 : ! update dS/dR matrix
961 16569800 : clist%acint(1:natorb_a, 1:natorb_b, 1:3) = sint(1:natorb_a, 1:natorb_b, 2:4)
962 :
963 390292 : DEALLOCATE (oint, owork, sint)
964 :
965 : END DO
966 118 : CALL neighbor_list_iterator_release(nl_iterator)
967 :
968 118 : DEALLOCATE (basis_set_list)
969 :
970 118 : CALL timestop(handle)
971 :
972 236 : END SUBROUTINE xtb_dsint_list
973 :
974 18187793 : END MODULE xtb_coulomb
|
