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