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