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 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 22802 : 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, blk, 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 22802 : 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 22802 : 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 22802 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xgamma, zeffk
128 22802 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gammab, gcij, gmcharge
129 22802 : 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 22802 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, ksblock, pblock, sblock
134 22802 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
135 22802 : 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 22802 : 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 22802 : DIMENSION(:), POINTER :: nl_iterator
148 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
149 22802 : POINTER :: n_list
150 22802 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
151 22802 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
152 22802 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
153 22802 : 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 22802 : CALL timeset(routineN, handle)
159 :
160 22802 : 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 22802 : dft_control=dft_control)
169 :
170 22802 : xtb_control => dft_control%qs_control%xtb_control
171 :
172 22802 : use_virial = .FALSE.
173 22802 : IF (calculate_forces) THEN
174 816 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
175 : END IF
176 :
177 22802 : do_gamma_stress = .FALSE.
178 22802 : IF (.NOT. just_energy .AND. use_virial) THEN
179 52 : IF (dft_control%nimages == 1) do_gamma_stress = .TRUE.
180 : END IF
181 :
182 22802 : 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 22802 : IF (calculate_forces) THEN
189 : nmat = 4
190 : ELSE
191 22368 : nmat = 1
192 : END IF
193 :
194 22802 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
195 91208 : ALLOCATE (gchrg(natom, 5, nmat))
196 1275736 : gchrg = 0._dp
197 91208 : ALLOCATE (gmcharge(natom, nmat))
198 268568 : gmcharge = 0._dp
199 :
200 : ! short range contribution (gamma)
201 : ! loop over all atom pairs (sab_xtbe)
202 22802 : kg = xtb_control%kg
203 22802 : NULLIFY (n_list)
204 22802 : CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
205 22802 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
206 7526496 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
207 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
208 7503694 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
209 7503694 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
210 7503694 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
211 7503694 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
212 7503687 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
213 7503687 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
214 7503687 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
215 : ! atomic parameters
216 7503680 : CALL get_xtb_atom_param(xtb_atom_a, eta=etaa, lmax=lmaxa, kappa=kappaa, rcut=rcuta)
217 7503680 : CALL get_xtb_atom_param(xtb_atom_b, eta=etab, lmax=lmaxb, kappa=kappab, rcut=rcutb)
218 : ! gamma matrix
219 7503680 : ni = lmaxa + 1
220 7503680 : nj = lmaxb + 1
221 30014720 : ALLOCATE (gammab(ni, nj))
222 7503680 : rcut = rcuta + rcutb
223 30014720 : dr = SQRT(SUM(rij(:)**2))
224 7503680 : CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
225 91525824 : gchrg(iatom, 1:ni, 1) = gchrg(iatom, 1:ni, 1) + MATMUL(gammab, charges(jatom, 1:nj))
226 7503680 : IF (iatom /= jatom) THEN
227 91378626 : gchrg(jatom, 1:nj, 1) = gchrg(jatom, 1:nj, 1) + MATMUL(charges(iatom, 1:ni), gammab)
228 : END IF
229 7503680 : IF (calculate_forces) THEN
230 281544 : IF (dr > 1.e-6_dp) THEN
231 279572 : CALL dgamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
232 1118288 : DO i = 1, 3
233 : gchrg(iatom, 1:ni, i + 1) = gchrg(iatom, 1:ni, i + 1) &
234 10220643 : + MATMUL(gammab, charges(jatom, 1:nj))*rij(i)/dr
235 1118288 : IF (iatom /= jatom) THEN
236 : gchrg(jatom, 1:nj, i + 1) = gchrg(jatom, 1:nj, i + 1) &
237 10554975 : - MATMUL(charges(iatom, 1:ni), gammab)*rij(i)/dr
238 : END IF
239 : END DO
240 279572 : IF (use_virial) THEN
241 1162578 : gcint(1:ni) = MATMUL(gammab, charges(jatom, 1:nj))
242 516772 : DO i = 1, 3
243 1094410 : fij(i) = -SUM(charges(iatom, 1:ni)*gcint(1:ni))*rij(i)/dr
244 : END DO
245 129193 : fi = 1.0_dp
246 129193 : IF (iatom == jatom) fi = 0.5_dp
247 129193 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
248 : END IF
249 : END IF
250 : END IF
251 30014741 : DEALLOCATE (gammab)
252 : END DO
253 22802 : CALL neighbor_list_iterator_release(nl_iterator)
254 :
255 : ! 1/R contribution
256 :
257 22802 : IF (xtb_control%coulomb_lr) THEN
258 22802 : do_ewald = xtb_control%do_ewald
259 22802 : IF (do_ewald) THEN
260 : ! Ewald sum
261 10022 : NULLIFY (ewald_env, ewald_pw)
262 : CALL get_qs_env(qs_env=qs_env, &
263 10022 : ewald_env=ewald_env, ewald_pw=ewald_pw)
264 10022 : CALL get_cell(cell=cell, periodic=periodic, deth=deth)
265 10022 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
266 10022 : CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
267 10022 : CALL tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial)
268 0 : SELECT CASE (ewald_type)
269 : CASE DEFAULT
270 0 : CPABORT("Invalid Ewald type")
271 : CASE (do_ewald_none)
272 0 : CPABORT("Not allowed with xTB/DFTB")
273 : CASE (do_ewald_ewald)
274 0 : CPABORT("Standard Ewald not implemented in xTB/DFTB")
275 : CASE (do_ewald_pme)
276 0 : CPABORT("PME not implemented in xTB/DFTB")
277 : CASE (do_ewald_spme)
278 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
279 10022 : gmcharge, mcharge, calculate_forces, virial, use_virial)
280 : END SELECT
281 : ELSE
282 : ! direct sum
283 : CALL get_qs_env(qs_env=qs_env, &
284 12780 : local_particles=local_particles)
285 49190 : DO ikind = 1, SIZE(local_particles%n_el)
286 116615 : DO ia = 1, local_particles%n_el(ikind)
287 67425 : iatom = local_particles%list(ikind)%array(ia)
288 834417 : DO jatom = 1, iatom - 1
289 2922328 : rij = particle_set(iatom)%r - particle_set(jatom)%r
290 2922328 : rij = pbc(rij, cell)
291 2922328 : dr = SQRT(SUM(rij(:)**2))
292 798007 : IF (dr > 1.e-6_dp) THEN
293 730582 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)/dr
294 730582 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)/dr
295 735778 : DO i = 2, nmat
296 5196 : gmcharge(iatom, i) = gmcharge(iatom, i) + rij(i - 1)*mcharge(jatom)/dr**3
297 735778 : gmcharge(jatom, i) = gmcharge(jatom, i) - rij(i - 1)*mcharge(iatom)/dr**3
298 : END DO
299 : END IF
300 : END DO
301 : END DO
302 : END DO
303 12780 : CPASSERT(.NOT. use_virial)
304 : END IF
305 : END IF
306 :
307 : ! global sum of gamma*p arrays
308 : CALL get_qs_env(qs_env=qs_env, &
309 : atomic_kind_set=atomic_kind_set, &
310 22802 : force=force, para_env=para_env)
311 22802 : CALL para_env%sum(gmcharge(:, 1))
312 22802 : CALL para_env%sum(gchrg(:, :, 1))
313 :
314 22802 : IF (xtb_control%coulomb_lr) THEN
315 22802 : IF (do_ewald) THEN
316 : ! add self charge interaction and background charge contribution
317 88878 : gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:)
318 10958 : IF (ANY(periodic(:) == 1)) THEN
319 87318 : gmcharge(:, 1) = gmcharge(:, 1) - pi/alpha**2/deth
320 : END IF
321 : END IF
322 : END IF
323 :
324 : ! energy
325 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
326 : kind_of=kind_of, &
327 22802 : atom_of_kind=atom_of_kind)
328 22802 : ecsr = 0.0_dp
329 233550 : DO iatom = 1, natom
330 210748 : ikind = kind_of(iatom)
331 210748 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
332 210748 : CALL get_xtb_atom_param(xtb_kind, lmax=ni)
333 210748 : ni = ni + 1
334 559506 : ecsr = ecsr + SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, 1))
335 : END DO
336 :
337 22802 : energy%hartree = energy%hartree + 0.5_dp*ecsr
338 233550 : energy%hartree = energy%hartree + 0.5_dp*SUM(mcharge(:)*gmcharge(:, 1))
339 :
340 22802 : IF (atprop%energy) THEN
341 172 : CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
342 748 : DO ikind = 1, SIZE(local_particles%n_el)
343 576 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
344 576 : CALL get_xtb_atom_param(xtb_kind, lmax=ni, occupation=occ)
345 576 : ni = ni + 1
346 3456 : zeff = SUM(REAL(occ, KIND=dp))
347 4360 : DO ia = 1, local_particles%n_el(ikind)
348 3036 : iatom = local_particles%list(ikind)%array(ia)
349 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
350 7258 : 0.5_dp*SUM(REAL(occ(1:ni), KIND=dp)*gchrg(iatom, 1:ni, 1))
351 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
352 3612 : 0.5_dp*zeff*gmcharge(iatom, 1)
353 : END DO
354 : END DO
355 : END IF
356 :
357 22802 : IF (calculate_forces) THEN
358 4072 : DO iatom = 1, natom
359 3638 : ikind = kind_of(iatom)
360 3638 : atom_i = atom_of_kind(iatom)
361 3638 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
362 3638 : CALL get_xtb_atom_param(xtb_kind, lmax=ni)
363 : ! short range
364 3638 : ni = ni + 1
365 14552 : DO i = 1, 3
366 30458 : fij(i) = SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, i + 1))
367 : END DO
368 3638 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
369 3638 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
370 3638 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
371 : ! long range
372 14552 : DO i = 1, 3
373 14552 : fij(i) = gmcharge(iatom, i + 1)*mcharge(iatom)
374 : END DO
375 3638 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
376 3638 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
377 7710 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
378 : END DO
379 : END IF
380 :
381 22802 : IF (.NOT. just_energy) THEN
382 22760 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
383 22760 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
384 :
385 22760 : nimg = dft_control%nimages
386 22760 : NULLIFY (cell_to_index)
387 22760 : IF (nimg > 1) THEN
388 2656 : NULLIFY (kpoints)
389 2656 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
390 2656 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
391 : END IF
392 :
393 22760 : IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
394 432 : DO img = 1, nimg
395 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
396 432 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
397 : END DO
398 : END IF
399 :
400 22760 : NULLIFY (sap_int)
401 22760 : IF (do_gamma_stress) THEN
402 : ! derivative overlap integral (non collapsed)
403 38 : CALL xtb_dsint_list(qs_env, sap_int)
404 : END IF
405 :
406 22760 : IF (nimg == 1) THEN
407 : ! no k-points; all matrices have been transformed to periodic bsf
408 20104 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
409 1405465 : DO WHILE (dbcsr_iterator_blocks_left(iter))
410 1385361 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
411 1385361 : ikind = kind_of(irow)
412 1385361 : jkind = kind_of(icol)
413 :
414 : ! atomic parameters
415 1385361 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
416 1385361 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
417 1385361 : CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
418 1385361 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
419 :
420 1385361 : ni = SIZE(sblock, 1)
421 1385361 : nj = SIZE(sblock, 2)
422 5541444 : ALLOCATE (gcij(ni, nj))
423 5613830 : DO i = 1, ni
424 18669412 : DO j = 1, nj
425 13055582 : la = laoa(i) + 1
426 13055582 : lb = laob(j) + 1
427 17284051 : gcij(i, j) = 0.5_dp*(gchrg(irow, la, 1) + gchrg(icol, lb, 1))
428 : END DO
429 : END DO
430 1385361 : gmij = 0.5_dp*(gmcharge(irow, 1) + gmcharge(icol, 1))
431 2779118 : DO is = 1, SIZE(ks_matrix, 1)
432 1393757 : NULLIFY (ksblock)
433 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
434 1393757 : row=irow, col=icol, block=ksblock, found=found)
435 1393757 : CPASSERT(found)
436 36153145 : ksblock = ksblock - gcij*sblock
437 38932263 : ksblock = ksblock - gmij*sblock
438 : END DO
439 1385361 : IF (calculate_forces) THEN
440 45390 : atom_i = atom_of_kind(irow)
441 45390 : atom_j = atom_of_kind(icol)
442 45390 : NULLIFY (pblock)
443 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
444 45390 : row=irow, col=icol, block=pblock, found=found)
445 45390 : CPASSERT(found)
446 181560 : DO i = 1, 3
447 136170 : NULLIFY (dsblock)
448 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
449 136170 : row=irow, col=icol, block=dsblock, found=found)
450 136170 : CPASSERT(found)
451 136170 : fij(i) = 0.0_dp
452 : ! short range
453 1540485 : fi = -2.0_dp*SUM(pblock*dsblock*gcij)
454 136170 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
455 136170 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
456 136170 : fij(i) = fij(i) + fi
457 : ! long range
458 1540485 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
459 136170 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
460 136170 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
461 317730 : fij(i) = fij(i) + fi
462 : END DO
463 : END IF
464 4176187 : DEALLOCATE (gcij)
465 : END DO
466 20104 : CALL dbcsr_iterator_stop(iter)
467 : ! stress tensor (needs recalculation of overlap integrals)
468 20104 : IF (do_gamma_stress) THEN
469 134 : DO ikind = 1, nkind
470 410 : DO jkind = 1, nkind
471 276 : iac = ikind + nkind*(jkind - 1)
472 276 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
473 : ! atomic parameters
474 156 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
475 156 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
476 156 : CALL get_xtb_atom_param(xtb_atom_a, lao=laoa, natorb=ni)
477 156 : CALL get_xtb_atom_param(xtb_atom_b, lao=laob, natorb=nj)
478 1244 : DO ia = 1, sap_int(iac)%nalist
479 992 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE
480 974 : iatom = sap_int(iac)%alist(ia)%aatom
481 70609 : DO ic = 1, sap_int(iac)%alist(ia)%nclist
482 69359 : jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
483 277436 : rij = sap_int(iac)%alist(ia)%clist(ic)%rac
484 277436 : dr = SQRT(SUM(rij(:)**2))
485 70351 : IF (dr > 1.e-6_dp) THEN
486 68879 : dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
487 275516 : ALLOCATE (gcij(ni, nj))
488 283876 : DO i = 1, ni
489 1178463 : DO j = 1, nj
490 894587 : la = laoa(i) + 1
491 894587 : lb = laob(j) + 1
492 1109584 : gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
493 : END DO
494 : END DO
495 68879 : gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
496 68879 : icol = MAX(iatom, jatom)
497 68879 : irow = MIN(iatom, jatom)
498 68879 : NULLIFY (pblock)
499 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
500 68879 : row=irow, col=icol, block=pblock, found=found)
501 68879 : CPASSERT(found)
502 68879 : fij = 0.0_dp
503 275516 : DO i = 1, 3
504 : ! short/long range
505 206637 : IF (irow == iatom) THEN
506 1957914 : f1 = -2.0_dp*SUM(pblock*dsint(:, :, i)*gcij)
507 1957914 : f2 = -2.0_dp*gmij*SUM(pblock*dsint(:, :, i))
508 : ELSE
509 1577181 : f1 = -2.0_dp*SUM(TRANSPOSE(pblock)*dsint(:, :, i)*gcij)
510 1577181 : f2 = -2.0_dp*gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i))
511 : END IF
512 275516 : fij(i) = f1 + f2
513 : END DO
514 68879 : DEALLOCATE (gcij)
515 68879 : fi = 1.0_dp
516 68879 : IF (iatom == jatom) fi = 0.5_dp
517 137758 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
518 : END IF
519 : END DO
520 : END DO
521 : END DO
522 : END DO
523 : END IF
524 : ELSE
525 2656 : NULLIFY (n_list)
526 2656 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
527 2656 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
528 1963764 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
529 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
530 1961108 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
531 :
532 1961108 : icol = MAX(iatom, jatom)
533 1961108 : irow = MIN(iatom, jatom)
534 :
535 1961108 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
536 1961108 : CPASSERT(ic > 0)
537 :
538 1961108 : NULLIFY (sblock)
539 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
540 1961108 : row=irow, col=icol, block=sblock, found=found)
541 1961108 : CPASSERT(found)
542 :
543 : ! atomic parameters
544 1961108 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
545 1961108 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
546 1961108 : CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
547 1961108 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
548 :
549 1961108 : ni = SIZE(sblock, 1)
550 1961108 : nj = SIZE(sblock, 2)
551 7844432 : ALLOCATE (gcij(ni, nj))
552 11442495 : DO i = 1, ni
553 72925194 : DO j = 1, nj
554 70964086 : IF (irow == iatom) THEN
555 34957817 : la = laoa(i) + 1
556 34957817 : lb = laob(j) + 1
557 34957817 : gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
558 : ELSE
559 26524882 : la = laoa(j) + 1
560 26524882 : lb = laob(i) + 1
561 26524882 : gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
562 : END IF
563 : END DO
564 : END DO
565 1961108 : gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
566 4109816 : DO is = 1, SIZE(ks_matrix, 1)
567 2148708 : NULLIFY (ksblock)
568 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
569 2148708 : row=irow, col=icol, block=ksblock, found=found)
570 2148708 : CPASSERT(found)
571 147117272 : ksblock = ksblock - gcij*sblock
572 151227088 : ksblock = ksblock - gmij*sblock
573 : END DO
574 :
575 1961108 : IF (calculate_forces) THEN
576 29400 : atom_i = atom_of_kind(iatom)
577 29400 : atom_j = atom_of_kind(jatom)
578 29400 : IF (irow /= iatom) THEN
579 11957 : gmij = -gmij
580 759507 : gcij = -gcij
581 : END IF
582 29400 : NULLIFY (pblock)
583 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
584 29400 : row=irow, col=icol, block=pblock, found=found)
585 29400 : CPASSERT(found)
586 117600 : DO i = 1, 3
587 88200 : NULLIFY (dsblock)
588 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
589 88200 : row=irow, col=icol, block=dsblock, found=found)
590 88200 : CPASSERT(found)
591 88200 : fij(i) = 0.0_dp
592 : ! short range
593 5846208 : fi = -2.0_dp*SUM(pblock*dsblock*gcij)
594 88200 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
595 88200 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
596 88200 : fij(i) = fij(i) + fi
597 : ! long range
598 5846208 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
599 88200 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
600 88200 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
601 205800 : fij(i) = fij(i) + fi
602 : END DO
603 29400 : IF (use_virial) THEN
604 73828 : dr = SQRT(SUM(rij(:)**2))
605 18457 : IF (dr > 1.e-6_dp) THEN
606 18393 : fi = 1.0_dp
607 18393 : IF (iatom == jatom) fi = 0.5_dp
608 18393 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
609 : END IF
610 : END IF
611 : END IF
612 5885980 : DEALLOCATE (gcij)
613 :
614 : END DO
615 2656 : CALL neighbor_list_iterator_release(nl_iterator)
616 : END IF
617 :
618 22760 : IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
619 432 : DO img = 1, nimg
620 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
621 432 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
622 : END DO
623 : END IF
624 : END IF
625 :
626 22802 : IF (xtb_control%tb3_interaction) THEN
627 22802 : CALL get_qs_env(qs_env, nkind=nkind)
628 91208 : ALLOCATE (zeffk(nkind), xgamma(nkind))
629 83610 : DO ikind = 1, nkind
630 60808 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
631 83610 : CALL get_xtb_atom_param(xtb_kind, xgamma=xgamma(ikind), zeff=zeffk(ikind))
632 : END DO
633 : ! Diagonal 3rd order correction (DFTB3)
634 : CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, &
635 22802 : sap_int, calculate_forces, just_energy)
636 22802 : DEALLOCATE (zeffk, xgamma)
637 : END IF
638 :
639 : ! QMMM
640 22802 : IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
641 : CALL build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, &
642 862 : calculate_forces, just_energy)
643 : END IF
644 :
645 22802 : IF (do_gamma_stress) THEN
646 38 : CALL release_sap_int(sap_int)
647 : END IF
648 :
649 22802 : CALL timestop(handle)
650 :
651 45604 : END SUBROUTINE build_xtb_coulomb
652 :
653 : ! **************************************************************************************************
654 : !> \brief Computes the short-range gamma parameter from
655 : !> Nataga-Mishimoto-Ohno-Klopman formula for xTB
656 : !> WARNING: The xTB function (gamma - 1/r) has still an l-dependent longrange (1/r^3)
657 : !> behaviour. We use a cutoff function to smoothly remove this part.
658 : !> However, this will change energies and effect final results.
659 : !>
660 : !> \param gmat ...
661 : !> \param rab ...
662 : !> \param nla ...
663 : !> \param kappaa ...
664 : !> \param etaa ...
665 : !> \param nlb ...
666 : !> \param kappab ...
667 : !> \param etab ...
668 : !> \param kg ...
669 : !> \param rcut ...
670 : !> \par History
671 : !> 10.2018 JGH
672 : !> \version 1.1
673 : ! **************************************************************************************************
674 7726207 : SUBROUTINE gamma_rab_sr(gmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
675 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: gmat
676 : REAL(dp), INTENT(IN) :: rab
677 : INTEGER, INTENT(IN) :: nla
678 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappaa
679 : REAL(dp), INTENT(IN) :: etaa
680 : INTEGER, INTENT(IN) :: nlb
681 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappab
682 : REAL(dp), INTENT(IN) :: etab, kg, rcut
683 :
684 : REAL(KIND=dp), PARAMETER :: rsmooth = 1.0_dp
685 :
686 : INTEGER :: i, j
687 : REAL(KIND=dp) :: fcut, r, rk, x
688 7726207 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eta
689 :
690 30904828 : ALLOCATE (eta(nla, nlb))
691 37335179 : eta = 0.0_dp
692 :
693 19096703 : DO j = 1, nlb
694 37335179 : DO i = 1, nla
695 18238476 : eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
696 29608972 : eta(i, j) = 2._dp/eta(i, j)
697 : END DO
698 : END DO
699 :
700 37335179 : gmat = 0.0_dp
701 7726207 : IF (rab < 1.e-6_dp) THEN
702 : ! on site terms
703 564002 : gmat(:, :) = eta(:, :)
704 7618267 : ELSEIF (rab > rcut) THEN
705 : ! do nothing
706 : ELSE
707 7618267 : rk = rab**kg
708 36771177 : eta = eta**(-kg)
709 7618267 : IF (rab < rcut - rsmooth) THEN
710 : fcut = 1.0_dp
711 : ELSE
712 847948 : r = rab - (rcut - rsmooth)
713 847948 : x = r/rsmooth
714 847948 : fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
715 : END IF
716 36771177 : gmat(:, :) = fcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg) - fcut/rab
717 : END IF
718 :
719 7726207 : DEALLOCATE (eta)
720 :
721 7726207 : END SUBROUTINE gamma_rab_sr
722 :
723 : ! **************************************************************************************************
724 : !> \brief Computes the derivative of the short-range gamma parameter from
725 : !> Nataga-Mishimoto-Ohno-Klopman formula for xTB
726 : !> WARNING: The xTB function (gamma - 1/r) has still an l-dependent longrange (1/r^3)
727 : !> behaviour. We use a cutoff function to smoothly remove this part.
728 : !> However, this will change energies and effect final results.
729 : !>
730 : !> \param dgmat ...
731 : !> \param rab ...
732 : !> \param nla ...
733 : !> \param kappaa ...
734 : !> \param etaa ...
735 : !> \param nlb ...
736 : !> \param kappab ...
737 : !> \param etab ...
738 : !> \param kg ...
739 : !> \param rcut ...
740 : !> \par History
741 : !> 10.2018 JGH
742 : !> \version 1.1
743 : ! **************************************************************************************************
744 304149 : SUBROUTINE dgamma_rab_sr(dgmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
745 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: dgmat
746 : REAL(dp), INTENT(IN) :: rab
747 : INTEGER, INTENT(IN) :: nla
748 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappaa
749 : REAL(dp), INTENT(IN) :: etaa
750 : INTEGER, INTENT(IN) :: nlb
751 : REAL(dp), DIMENSION(:), INTENT(IN) :: kappab
752 : REAL(dp), INTENT(IN) :: etab, kg, rcut
753 :
754 : REAL(KIND=dp), PARAMETER :: rsmooth = 1.0_dp
755 :
756 : INTEGER :: i, j
757 : REAL(KIND=dp) :: dfcut, fcut, r, rk, x
758 304149 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eta
759 :
760 1216596 : ALLOCATE (eta(nla, nlb))
761 :
762 744637 : DO j = 1, nlb
763 1455665 : DO i = 1, nla
764 711028 : eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
765 1151516 : eta(i, j) = 2._dp/eta(i, j)
766 : END DO
767 : END DO
768 :
769 304149 : IF (rab < 1.e-6) THEN
770 : ! on site terms
771 0 : dgmat(:, :) = 0.0_dp
772 304149 : ELSEIF (rab > rcut) THEN
773 0 : dgmat(:, :) = 0.0_dp
774 : ELSE
775 1455665 : eta = eta**(-kg)
776 304149 : rk = rab**kg
777 304149 : IF (rab < rcut - rsmooth) THEN
778 : fcut = 1.0_dp
779 : dfcut = 0.0_dp
780 : ELSE
781 40748 : r = rab - (rcut - rsmooth)
782 40748 : x = r/rsmooth
783 40748 : fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
784 40748 : dfcut = -30._dp*x**4 + 60._dp*x**3 - 30._dp*x**2
785 40748 : dfcut = dfcut/rsmooth
786 : END IF
787 1455665 : dgmat(:, :) = dfcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg)
788 1455665 : dgmat(:, :) = dgmat(:, :) - dfcut/rab + fcut/rab**2
789 1455665 : dgmat(:, :) = dgmat(:, :) - fcut/(rk + eta(:, :))*(1._dp/(rk + eta(:, :)))**(1._dp/kg)*rk/rab
790 : END IF
791 :
792 304149 : DEALLOCATE (eta)
793 :
794 304149 : END SUBROUTINE dgamma_rab_sr
795 :
796 : ! **************************************************************************************************
797 : !> \brief ...
798 : !> \param qs_env ...
799 : !> \param sap_int ...
800 : ! **************************************************************************************************
801 44 : SUBROUTINE xtb_dsint_list(qs_env, sap_int)
802 :
803 : TYPE(qs_environment_type), POINTER :: qs_env
804 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
805 :
806 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xtb_dsint_list'
807 :
808 : INTEGER :: handle, i, iac, iatom, ikind, ilist, iset, jatom, jkind, jneighbor, jset, ldsab, &
809 : n1, n2, natorb_a, natorb_b, ncoa, ncob, nkind, nlist, nneighbor, nseta, nsetb, sgfa, sgfb
810 : INTEGER, DIMENSION(3) :: cell
811 44 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
812 44 : npgfb, nsgfa, nsgfb
813 44 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
814 : LOGICAL :: defined
815 : REAL(KIND=dp) :: dr
816 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
817 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
818 : REAL(KIND=dp), DIMENSION(3) :: rij
819 44 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
820 44 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb
821 : TYPE(clist_type), POINTER :: clist
822 : TYPE(dft_control_type), POINTER :: dft_control
823 44 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
824 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
825 : TYPE(neighbor_list_iterator_p_type), &
826 44 : DIMENSION(:), POINTER :: nl_iterator
827 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
828 44 : POINTER :: sab_orb
829 44 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
830 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
831 :
832 44 : CALL timeset(routineN, handle)
833 :
834 44 : CALL get_qs_env(qs_env=qs_env, nkind=nkind)
835 44 : CPASSERT(.NOT. ASSOCIATED(sap_int))
836 462 : ALLOCATE (sap_int(nkind*nkind))
837 374 : DO i = 1, nkind*nkind
838 330 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
839 374 : sap_int(i)%nalist = 0
840 : END DO
841 :
842 : CALL get_qs_env(qs_env=qs_env, &
843 : qs_kind_set=qs_kind_set, &
844 : dft_control=dft_control, &
845 44 : sab_orb=sab_orb)
846 :
847 : ! set up basis set lists
848 246 : ALLOCATE (basis_set_list(nkind))
849 44 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
850 :
851 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
852 44 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
853 70785 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
854 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, &
855 : jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, &
856 70741 : inode=jneighbor, cell=cell, r=rij)
857 70741 : iac = ikind + nkind*(jkind - 1)
858 : !
859 70741 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
860 70741 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
861 70741 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
862 70741 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
863 70741 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
864 70741 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
865 :
866 282964 : dr = SQRT(SUM(rij(:)**2))
867 :
868 : ! integral list
869 70741 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
870 183 : sap_int(iac)%a_kind = ikind
871 183 : sap_int(iac)%p_kind = jkind
872 183 : sap_int(iac)%nalist = nlist
873 1568 : ALLOCATE (sap_int(iac)%alist(nlist))
874 1202 : DO i = 1, nlist
875 1019 : NULLIFY (sap_int(iac)%alist(i)%clist)
876 1019 : sap_int(iac)%alist(i)%aatom = 0
877 1202 : sap_int(iac)%alist(i)%nclist = 0
878 : END DO
879 : END IF
880 70741 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
881 1001 : sap_int(iac)%alist(ilist)%aatom = iatom
882 1001 : sap_int(iac)%alist(ilist)%nclist = nneighbor
883 79750 : ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
884 71742 : DO i = 1, nneighbor
885 71742 : sap_int(iac)%alist(ilist)%clist(i)%catom = 0
886 : END DO
887 : END IF
888 70741 : clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
889 70741 : clist%catom = jatom
890 282964 : clist%cell = cell
891 282964 : clist%rac = rij
892 353705 : ALLOCATE (clist%acint(natorb_a, natorb_b, 3))
893 70741 : NULLIFY (clist%achint)
894 3677812 : clist%acint = 0._dp
895 70741 : clist%nsgf_cnt = 0
896 70741 : NULLIFY (clist%sgf_list)
897 :
898 : ! overlap
899 70741 : basis_set_a => basis_set_list(ikind)%gto_basis_set
900 70741 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
901 70741 : basis_set_b => basis_set_list(jkind)%gto_basis_set
902 70741 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
903 : ! basis ikind
904 70741 : first_sgfa => basis_set_a%first_sgf
905 70741 : la_max => basis_set_a%lmax
906 70741 : la_min => basis_set_a%lmin
907 70741 : npgfa => basis_set_a%npgf
908 70741 : nseta = basis_set_a%nset
909 70741 : nsgfa => basis_set_a%nsgf_set
910 70741 : rpgfa => basis_set_a%pgf_radius
911 70741 : set_radius_a => basis_set_a%set_radius
912 70741 : scon_a => basis_set_a%scon
913 70741 : zeta => basis_set_a%zet
914 : ! basis jkind
915 70741 : first_sgfb => basis_set_b%first_sgf
916 70741 : lb_max => basis_set_b%lmax
917 70741 : lb_min => basis_set_b%lmin
918 70741 : npgfb => basis_set_b%npgf
919 70741 : nsetb = basis_set_b%nset
920 70741 : nsgfb => basis_set_b%nsgf_set
921 70741 : rpgfb => basis_set_b%pgf_radius
922 70741 : set_radius_b => basis_set_b%set_radius
923 70741 : scon_b => basis_set_b%scon
924 70741 : zetb => basis_set_b%zet
925 :
926 70741 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
927 565928 : ALLOCATE (oint(ldsab, ldsab, 4), owork(ldsab, ldsab))
928 353705 : ALLOCATE (sint(natorb_a, natorb_b, 4))
929 4880169 : sint = 0.0_dp
930 :
931 218610 : DO iset = 1, nseta
932 147869 : ncoa = npgfa(iset)*ncoset(la_max(iset))
933 147869 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
934 147869 : sgfa = first_sgfa(1, iset)
935 532679 : DO jset = 1, nsetb
936 314069 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
937 258157 : ncob = npgfb(jset)*ncoset(lb_max(jset))
938 258157 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
939 258157 : sgfb = first_sgfb(1, jset)
940 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
941 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
942 258157 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
943 : ! Contraction
944 1438654 : DO i = 1, 4
945 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
946 1032628 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
947 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), &
948 1346697 : sgfa, sgfb, trans=.FALSE.)
949 : END DO
950 : END DO
951 : END DO
952 : ! update dS/dR matrix
953 3677812 : clist%acint(1:natorb_a, 1:natorb_b, 1:3) = sint(1:natorb_a, 1:natorb_b, 2:4)
954 :
955 212267 : DEALLOCATE (oint, owork, sint)
956 :
957 : END DO
958 44 : CALL neighbor_list_iterator_release(nl_iterator)
959 :
960 44 : DEALLOCATE (basis_set_list)
961 :
962 44 : CALL timestop(handle)
963 :
964 88 : END SUBROUTINE xtb_dsint_list
965 :
966 15552436 : END MODULE xtb_coulomb
967 :
|
