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 DFTB 8!> \author JGH 9! ************************************************************************************************** 10MODULE qs_dftb_coulomb 11 USE atomic_kind_types, ONLY: atomic_kind_type,& 12 get_atomic_kind_set,& 13 is_hydrogen 14 USE atprop_types, ONLY: atprop_array_init,& 15 atprop_type 16 USE cell_types, ONLY: cell_type,& 17 get_cell,& 18 pbc 19 USE cp_control_types, ONLY: dft_control_type,& 20 dftb_control_type 21 USE cp_para_types, ONLY: cp_para_env_type 22 USE dbcsr_api, ONLY: dbcsr_add,& 23 dbcsr_get_block_p,& 24 dbcsr_iterator_blocks_left,& 25 dbcsr_iterator_next_block,& 26 dbcsr_iterator_start,& 27 dbcsr_iterator_stop,& 28 dbcsr_iterator_type,& 29 dbcsr_p_type 30 USE distribution_1d_types, ONLY: distribution_1d_type 31 USE ewald_environment_types, ONLY: ewald_env_get,& 32 ewald_environment_type 33 USE ewald_methods_tb, ONLY: tb_ewald_overlap,& 34 tb_spme_evaluate 35 USE ewald_pw_types, ONLY: ewald_pw_type 36 USE kinds, ONLY: dp 37 USE kpoint_types, ONLY: get_kpoint_info,& 38 kpoint_type 39 USE mathconstants, ONLY: oorootpi,& 40 pi 41 USE message_passing, ONLY: mp_sum 42 USE particle_types, ONLY: particle_type 43 USE pw_poisson_types, ONLY: do_ewald_ewald,& 44 do_ewald_none,& 45 do_ewald_pme,& 46 do_ewald_spme 47 USE qmmm_tb_coulomb, ONLY: build_tb_coulomb_qmqm 48 USE qs_dftb3_methods, ONLY: build_dftb3_diagonal 49 USE qs_dftb_types, ONLY: qs_dftb_atom_type,& 50 qs_dftb_pairpot_type 51 USE qs_dftb_utils, ONLY: compute_block_sk,& 52 get_dftb_atom_param 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_kind_types, ONLY: get_qs_kind,& 58 qs_kind_type 59 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 60 neighbor_list_iterate,& 61 neighbor_list_iterator_create,& 62 neighbor_list_iterator_p_type,& 63 neighbor_list_iterator_release,& 64 neighbor_list_set_p_type 65 USE qs_rho_types, ONLY: qs_rho_get,& 66 qs_rho_type 67 USE sap_kind_types, ONLY: clist_type,& 68 release_sap_int,& 69 sap_int_type 70 USE virial_methods, ONLY: virial_pair_force 71 USE virial_types, ONLY: virial_type 72#include "./base/base_uses.f90" 73 74 IMPLICIT NONE 75 76 PRIVATE 77 78 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_coulomb' 79 80 ! screening for gamma function 81 REAL(dp), PARAMETER :: tol_gamma = 1.e-4_dp 82 ! small real number 83 REAL(dp), PARAMETER :: rtiny = 1.e-10_dp 84 85 PUBLIC :: build_dftb_coulomb, gamma_rab_sr 86 87CONTAINS 88 89! ************************************************************************************************** 90!> \brief ... 91!> \param qs_env ... 92!> \param ks_matrix ... 93!> \param rho ... 94!> \param mcharge ... 95!> \param energy ... 96!> \param calculate_forces ... 97!> \param just_energy ... 98! ************************************************************************************************** 99 SUBROUTINE build_dftb_coulomb(qs_env, ks_matrix, rho, 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(:) :: mcharge 106 TYPE(qs_energy_type), POINTER :: energy 107 LOGICAL, INTENT(in) :: calculate_forces, just_energy 108 109 CHARACTER(len=*), PARAMETER :: routineN = 'build_dftb_coulomb', & 110 routineP = moduleN//':'//routineN 111 112 INTEGER :: atom_i, atom_j, blk, ewald_type, handle, i, ia, iac, iatom, ic, icol, ikind, img, & 113 irow, is, jatom, jkind, natom, natorb_a, natorb_b, nimg, nkind, nmat 114 INTEGER, DIMENSION(3) :: cellind, periodic 115 INTEGER, DIMENSION(:), POINTER :: atom_of_kind, kind_of 116 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 117 LOGICAL :: defined, do_ewald, do_gamma_stress, & 118 found, hb_sr_damp, use_virial 119 REAL(KIND=dp) :: alpha, ddr, deth, dgam, dr, drm, drp, & 120 fi, ga, gb, gmat, gmij, hb_para, zeff 121 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xgamma, zeffk 122 REAL(KIND=dp), DIMENSION(0:3) :: eta_a, eta_b 123 REAL(KIND=dp), DIMENSION(3) :: fij, rij 124 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, gmcharge, ksblock, pblock, & 125 sblock 126 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint 127 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 128 TYPE(atprop_type), POINTER :: atprop 129 TYPE(cell_type), POINTER :: cell 130 TYPE(cp_para_env_type), POINTER :: para_env 131 TYPE(dbcsr_iterator_type) :: iter 132 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s 133 TYPE(dft_control_type), POINTER :: dft_control 134 TYPE(distribution_1d_type), POINTER :: local_particles 135 TYPE(ewald_environment_type), POINTER :: ewald_env 136 TYPE(ewald_pw_type), POINTER :: ewald_pw 137 TYPE(kpoint_type), POINTER :: kpoints 138 TYPE(neighbor_list_iterator_p_type), & 139 DIMENSION(:), POINTER :: nl_iterator 140 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 141 POINTER :: n_list 142 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 143 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind, dftb_kind_a, dftb_kind_b 144 TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), & 145 POINTER :: dftb_potential 146 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 147 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 148 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int 149 TYPE(virial_type), POINTER :: virial 150 151 CALL timeset(routineN, handle) 152 153 NULLIFY (matrix_p, matrix_s, virial, atprop, dft_control) 154 155 use_virial = .FALSE. 156 157 IF (calculate_forces) THEN 158 nmat = 4 159 ELSE 160 nmat = 1 161 END IF 162 163 CALL get_qs_env(qs_env, nkind=nkind, natom=natom) 164 ALLOCATE (gmcharge(natom, nmat)) 165 gmcharge = 0._dp 166 167 CALL get_qs_env(qs_env, & 168 particle_set=particle_set, & 169 cell=cell, & 170 virial=virial, & 171 atprop=atprop, & 172 dft_control=dft_control, & 173 atomic_kind_set=atomic_kind_set, & 174 qs_kind_set=qs_kind_set, & 175 force=force, para_env=para_env) 176 177 IF (calculate_forces) THEN 178 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 179 END IF 180 181 NULLIFY (dftb_potential) 182 CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential) 183 NULLIFY (n_list) 184 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) 185 CALL neighbor_list_iterator_create(nl_iterator, n_list) 186 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 187 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 188 iatom=iatom, jatom=jatom, r=rij, cell=cellind) 189 190 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a) 191 CALL get_dftb_atom_param(dftb_kind_a, & 192 defined=defined, eta=eta_a, natorb=natorb_a) 193 IF (.NOT. defined .OR. natorb_a < 1) CYCLE 194 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b) 195 CALL get_dftb_atom_param(dftb_kind_b, & 196 defined=defined, eta=eta_b, natorb=natorb_b) 197 IF (.NOT. defined .OR. natorb_b < 1) CYCLE 198 199 ! gamma matrix 200 hb_sr_damp = dft_control%qs_control%dftb_control%hb_sr_damp 201 IF (hb_sr_damp) THEN 202 ! short range correction enabled only when iatom XOR jatom are hydrogens 203 hb_sr_damp = is_hydrogen(particle_set(iatom)%atomic_kind) .NEQV. & 204 is_hydrogen(particle_set(jatom)%atomic_kind) 205 END IF 206 IF (hb_sr_damp) THEN 207 hb_para = dft_control%qs_control%dftb_control%hb_sr_para 208 ELSE 209 hb_para = 0.0_dp 210 END IF 211 ga = eta_a(0) 212 gb = eta_b(0) 213 dr = SQRT(SUM(rij(:)**2)) 214 gmat = gamma_rab_sr(dr, ga, gb, hb_para) 215 gmcharge(jatom, 1) = gmcharge(jatom, 1) + gmat*mcharge(iatom) 216 IF (iatom /= jatom) THEN 217 gmcharge(iatom, 1) = gmcharge(iatom, 1) + gmat*mcharge(jatom) 218 END IF 219 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN 220 ddr = 0.1_dp*dftb_potential(ikind, jkind)%dgrd 221 drp = dr + ddr 222 drm = dr - ddr 223 dgam = 0.5_dp*(gamma_rab_sr(drp, ga, gb, hb_para) - gamma_rab_sr(drm, ga, gb, hb_para))/ddr 224 DO i = 1, 3 225 gmcharge(jatom, i + 1) = gmcharge(jatom, i + 1) - dgam*mcharge(iatom)*rij(i)/dr 226 IF (dr > 0.001_dp) THEN 227 gmcharge(iatom, i + 1) = gmcharge(iatom, i + 1) + dgam*mcharge(jatom)*rij(i)/dr 228 END IF 229 IF (use_virial) THEN 230 fij(i) = -mcharge(iatom)*mcharge(jatom)*dgam*rij(i)/dr 231 END IF 232 END DO 233 IF (use_virial) THEN 234 fi = 1.0_dp 235 IF (iatom == jatom) fi = 0.5_dp 236 CALL virial_pair_force(virial%pv_virial, fi, fij, rij) 237 IF (atprop%stress) THEN 238 CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij) 239 CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij) 240 END IF 241 END IF 242 END IF 243 END DO 244 CALL neighbor_list_iterator_release(nl_iterator) 245 246 IF (atprop%energy) THEN 247 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set) 248 natom = SIZE(particle_set) 249 CALL atprop_array_init(atprop%atecoul, natom) 250 END IF 251 252 ! 1/R contribution 253 do_ewald = dft_control%qs_control%dftb_control%do_ewald 254 IF (do_ewald) THEN 255 ! Ewald sum 256 NULLIFY (ewald_env, ewald_pw) 257 CALL get_qs_env(qs_env=qs_env, & 258 ewald_env=ewald_env, ewald_pw=ewald_pw) 259 CALL get_cell(cell=cell, periodic=periodic, deth=deth) 260 CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type) 261 CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list) 262 CALL tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, & 263 virial, use_virial, atprop) 264 SELECT CASE (ewald_type) 265 CASE DEFAULT 266 CPABORT("Invalid Ewald type") 267 CASE (do_ewald_none) 268 CPABORT("Not allowed with DFTB") 269 CASE (do_ewald_ewald) 270 CPABORT("Standard Ewald not implemented in DFTB") 271 CASE (do_ewald_pme) 272 CPABORT("PME not implemented in DFTB") 273 CASE (do_ewald_spme) 274 CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, & 275 gmcharge, mcharge, calculate_forces, virial, use_virial, atprop) 276 END SELECT 277 ELSE 278 ! direct sum 279 CALL get_qs_env(qs_env=qs_env, & 280 local_particles=local_particles) 281 DO ikind = 1, SIZE(local_particles%n_el) 282 DO ia = 1, local_particles%n_el(ikind) 283 iatom = local_particles%list(ikind)%array(ia) 284 DO jatom = 1, iatom - 1 285 rij = particle_set(iatom)%r - particle_set(jatom)%r 286 rij = pbc(rij, cell) 287 dr = SQRT(SUM(rij(:)**2)) 288 gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)/dr 289 gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)/dr 290 DO i = 2, nmat 291 gmcharge(iatom, i) = gmcharge(iatom, i) + rij(i - 1)*mcharge(jatom)/dr**3 292 gmcharge(jatom, i) = gmcharge(jatom, i) - rij(i - 1)*mcharge(iatom)/dr**3 293 END DO 294 END DO 295 END DO 296 END DO 297 CPASSERT(.NOT. use_virial) 298 END IF 299 300 CALL mp_sum(gmcharge(:, 1), para_env%group) 301 302 IF (do_ewald) THEN 303 ! add self charge interaction and background charge contribution 304 gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:) 305 IF (ANY(periodic(:) == 1)) THEN 306 gmcharge(:, 1) = gmcharge(:, 1) - pi/alpha**2/deth 307 END IF 308 END IF 309 310 energy%hartree = energy%hartree + 0.5_dp*SUM(mcharge(:)*gmcharge(:, 1)) 311 IF (atprop%energy) THEN 312 CALL get_qs_env(qs_env=qs_env, & 313 local_particles=local_particles) 314 DO ikind = 1, SIZE(local_particles%n_el) 315 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind) 316 CALL get_dftb_atom_param(dftb_kind, zeff=zeff) 317 DO ia = 1, local_particles%n_el(ikind) 318 iatom = local_particles%list(ikind)%array(ia) 319 atprop%atecoul(iatom) = atprop%atecoul(iatom) + & 320 0.5_dp*zeff*gmcharge(iatom, 1) 321 END DO 322 END DO 323 END IF 324 325 IF (calculate_forces) THEN 326 ALLOCATE (atom_of_kind(natom), kind_of(natom)) 327 328 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 329 kind_of=kind_of, & 330 atom_of_kind=atom_of_kind) 331 332 gmcharge(:, 2) = gmcharge(:, 2)*mcharge(:) 333 gmcharge(:, 3) = gmcharge(:, 3)*mcharge(:) 334 gmcharge(:, 4) = gmcharge(:, 4)*mcharge(:) 335 DO iatom = 1, natom 336 ikind = kind_of(iatom) 337 atom_i = atom_of_kind(iatom) 338 force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - gmcharge(iatom, 2) 339 force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - gmcharge(iatom, 3) 340 force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - gmcharge(iatom, 4) 341 END DO 342 END IF 343 344 do_gamma_stress = .FALSE. 345 IF (.NOT. just_energy .AND. use_virial) THEN 346 IF (dft_control%nimages == 1) do_gamma_stress = .TRUE. 347 END IF 348 349 IF (.NOT. just_energy) THEN 350 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s) 351 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 352 353 nimg = dft_control%nimages 354 NULLIFY (cell_to_index) 355 IF (nimg > 1) THEN 356 NULLIFY (kpoints) 357 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints) 358 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 359 END IF 360 361 IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN 362 DO img = 1, nimg 363 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, & 364 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 365 END DO 366 END IF 367 368 NULLIFY (sap_int) 369 IF (do_gamma_stress) THEN 370 ! derivative overlap integral (non collapsed) 371 CALL dftb_dsint_list(qs_env, sap_int) 372 END IF 373 374 IF (nimg == 1) THEN 375 ! no k-points; all matrices have been transformed to periodic bsf 376 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix) 377 DO WHILE (dbcsr_iterator_blocks_left(iter)) 378 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk) 379 gmij = 0.5_dp*(gmcharge(irow, 1) + gmcharge(icol, 1)) 380 DO is = 1, SIZE(ks_matrix, 1) 381 NULLIFY (ksblock) 382 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, & 383 row=irow, col=icol, block=ksblock, found=found) 384 CPASSERT(found) 385 ksblock = ksblock - gmij*sblock 386 END DO 387 IF (calculate_forces) THEN 388 ikind = kind_of(irow) 389 atom_i = atom_of_kind(irow) 390 jkind = kind_of(icol) 391 atom_j = atom_of_kind(icol) 392 NULLIFY (pblock) 393 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, & 394 row=irow, col=icol, block=pblock, found=found) 395 CPASSERT(found) 396 DO i = 1, 3 397 NULLIFY (dsblock) 398 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, & 399 row=irow, col=icol, block=dsblock, found=found) 400 CPASSERT(found) 401 fi = -2.0_dp*gmij*SUM(pblock*dsblock) 402 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi 403 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi 404 fij(i) = fi 405 END DO 406 END IF 407 ENDDO 408 CALL dbcsr_iterator_stop(iter) 409 ! 410 IF (do_gamma_stress) THEN 411 DO ikind = 1, nkind 412 DO jkind = 1, nkind 413 iac = ikind + nkind*(jkind - 1) 414 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE 415 DO ia = 1, sap_int(iac)%nalist 416 IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE 417 iatom = sap_int(iac)%alist(ia)%aatom 418 DO ic = 1, sap_int(iac)%alist(ia)%nclist 419 jatom = sap_int(iac)%alist(ia)%clist(ic)%catom 420 rij = sap_int(iac)%alist(ia)%clist(ic)%rac 421 dr = SQRT(SUM(rij(:)**2)) 422 IF (dr > 1.e-6_dp) THEN 423 dsint => sap_int(iac)%alist(ia)%clist(ic)%acint 424 gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1)) 425 icol = MAX(iatom, jatom) 426 irow = MIN(iatom, jatom) 427 NULLIFY (pblock) 428 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, & 429 row=irow, col=icol, block=pblock, found=found) 430 CPASSERT(found) 431 DO i = 1, 3 432 IF (irow == iatom) THEN 433 fij(i) = -2.0_dp*gmij*SUM(pblock*dsint(:, :, i)) 434 ELSE 435 fij(i) = -2.0_dp*gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i)) 436 END IF 437 END DO 438 fi = 1.0_dp 439 IF (iatom == jatom) fi = 0.5_dp 440 CALL virial_pair_force(virial%pv_virial, fi, fij, rij) 441 IF (atprop%stress) THEN 442 CALL virial_pair_force(atprop%atstress(:, :, iatom), 0.5_dp*fi, fij, rij) 443 CALL virial_pair_force(atprop%atstress(:, :, jatom), 0.5_dp*fi, fij, rij) 444 END IF 445 END IF 446 END DO 447 END DO 448 END DO 449 END DO 450 END IF 451 ELSE 452 NULLIFY (n_list) 453 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) 454 CALL neighbor_list_iterator_create(nl_iterator, n_list) 455 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 456 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 457 iatom=iatom, jatom=jatom, r=rij, cell=cellind) 458 459 icol = MAX(iatom, jatom) 460 irow = MIN(iatom, jatom) 461 ic = cell_to_index(cellind(1), cellind(2), cellind(3)) 462 CPASSERT(ic > 0) 463 464 gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1)) 465 466 NULLIFY (sblock) 467 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, & 468 row=irow, col=icol, block=sblock, found=found) 469 CPASSERT(found) 470 DO is = 1, SIZE(ks_matrix, 1) 471 NULLIFY (ksblock) 472 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, & 473 row=irow, col=icol, block=ksblock, found=found) 474 CPASSERT(found) 475 ksblock = ksblock - gmij*sblock 476 END DO 477 478 IF (calculate_forces) THEN 479 ikind = kind_of(iatom) 480 atom_i = atom_of_kind(iatom) 481 jkind = kind_of(jatom) 482 atom_j = atom_of_kind(jatom) 483 IF (irow == jatom) gmij = -gmij 484 NULLIFY (pblock) 485 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, & 486 row=irow, col=icol, block=pblock, found=found) 487 CPASSERT(found) 488 DO i = 1, 3 489 NULLIFY (dsblock) 490 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, & 491 row=irow, col=icol, block=dsblock, found=found) 492 CPASSERT(found) 493 fi = -2.0_dp*gmij*SUM(pblock*dsblock) 494 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi 495 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi 496 fij(i) = fi 497 END DO 498 IF (use_virial) THEN 499 fi = 1.0_dp 500 IF (iatom == jatom) fi = 0.5_dp 501 CALL virial_pair_force(virial%pv_virial, fi, fij, rij) 502 IF (atprop%stress) THEN 503 CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij) 504 CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij) 505 END IF 506 END IF 507 END IF 508 END DO 509 CALL neighbor_list_iterator_release(nl_iterator) 510 END IF 511 512 IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN 513 DO img = 1, nimg 514 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, & 515 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp) 516 END DO 517 END IF 518 END IF 519 520 IF (dft_control%qs_control%dftb_control%dftb3_diagonal) THEN 521 CALL get_qs_env(qs_env, nkind=nkind) 522 ALLOCATE (zeffk(nkind), xgamma(nkind)) 523 DO ikind = 1, nkind 524 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind) 525 CALL get_dftb_atom_param(dftb_kind, dudq=xgamma(ikind), zeff=zeffk(ikind)) 526 END DO 527 ! Diagonal 3rd order correction (DFTB3) 528 CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, & 529 sap_int, calculate_forces, just_energy) 530 DEALLOCATE (zeffk, xgamma) 531 END IF 532 533 ! QMMM 534 IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN 535 CALL build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, & 536 calculate_forces, just_energy) 537 END IF 538 539 IF (do_gamma_stress) THEN 540 CALL release_sap_int(sap_int) 541 END IF 542 543 IF (calculate_forces) THEN 544 DEALLOCATE (atom_of_kind, kind_of) 545 END IF 546 DEALLOCATE (gmcharge) 547 548 CALL timestop(handle) 549 550 END SUBROUTINE build_dftb_coulomb 551 552! ************************************************************************************************** 553!> \brief Computes the short-range gamma parameter from exact Coulomb 554!> interaction of normalized exp(-a*r) charge distribution - 1/r 555!> \param r ... 556!> \param ga ... 557!> \param gb ... 558!> \param hb_para ... 559!> \return ... 560!> \par Literature 561!> Elstner et al, PRB 58 (1998) 7260 562!> \par History 563!> 10.2008 Axel Kohlmeyer - adding sr_damp 564!> 08.2014 JGH - adding flexible exponent for damping 565!> \version 1.1 566! ************************************************************************************************** 567 FUNCTION gamma_rab_sr(r, ga, gb, hb_para) RESULT(gamma) 568 REAL(dp), INTENT(in) :: r, ga, gb, hb_para 569 REAL(dp) :: gamma 570 571 REAL(dp) :: a, b, fac, g_sum 572 573 gamma = 0.0_dp 574 a = 3.2_dp*ga ! 3.2 = 16/5 in Eq. 18 and ff. 575 b = 3.2_dp*gb 576 g_sum = a + b 577 IF (g_sum < tol_gamma) RETURN ! hardness screening 578 IF (r < rtiny) THEN ! This is for short distances but non-onsite terms 579 ! This gives also correct diagonal elements (a=b, r=0) 580 gamma = 0.5_dp*(a*b/g_sum + (a*b)**2/g_sum**3) 581 RETURN 582 END IF 583 ! 584 ! distinguish two cases: Gamma's are very close, e.g. for the same atom type, 585 ! and Gamma's are different 586 ! 587 IF (ABS(a - b) < rtiny) THEN 588 fac = 1.6_dp*r*a*b/g_sum*(1.0_dp + a*b/g_sum**2) 589 gamma = -(48.0_dp + 33._dp*fac + (9.0_dp + fac)*fac**2)*EXP(-fac)/(48._dp*r) 590 ELSE 591 gamma = -EXP(-a*r)*(0.5_dp*a*b**4/(a**2 - b**2)**2 - & 592 (b**6 - 3._dp*a**2*b**4)/(r*(a**2 - b**2)**3)) - & ! a-> b 593 EXP(-b*r)*(0.5_dp*b*a**4/(b**2 - a**2)**2 - & 594 (a**6 - 3._dp*b**2*a**4)/(r*(b**2 - a**2)**3)) ! b-> a 595 END IF 596 ! 597 ! damping function for better short range hydrogen bonds. 598 ! functional form from Hu H. et al., J. Phys. Chem. A 2007, 111, 5685-5691 599 ! according to Elstner M, Theor. Chem. Acc. 2006, 116, 316-325, 600 ! this should only be applied to a-b pairs involving hydrogen. 601 IF (hb_para > 0.0_dp) THEN 602 gamma = gamma*EXP(-(0.5_dp*(ga + gb))**hb_para*r*r) 603 END IF 604 END FUNCTION gamma_rab_sr 605 606! ************************************************************************************************** 607!> \brief ... 608!> \param qs_env ... 609!> \param sap_int ... 610! ************************************************************************************************** 611 SUBROUTINE dftb_dsint_list(qs_env, sap_int) 612 613 TYPE(qs_environment_type), POINTER :: qs_env 614 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int 615 616 CHARACTER(LEN=*), PARAMETER :: routineN = 'dftb_dsint_list', & 617 routineP = moduleN//':'//routineN 618 619 INTEGER :: handle, i, iac, iatom, ikind, ilist, jatom, jkind, jneighbor, llm, lmaxi, lmaxj, & 620 n1, n2, natorb_a, natorb_b, ngrd, ngrdcut, nkind, nlist, nneighbor 621 INTEGER, DIMENSION(3) :: cell 622 LOGICAL :: defined 623 REAL(KIND=dp) :: ddr, dgrd, dr 624 REAL(KIND=dp), DIMENSION(3) :: drij, rij 625 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, dsblockm, smatij, smatji 626 TYPE(clist_type), POINTER :: clist 627 TYPE(dft_control_type), POINTER :: dft_control 628 TYPE(dftb_control_type), POINTER :: dftb_control 629 TYPE(neighbor_list_iterator_p_type), & 630 DIMENSION(:), POINTER :: nl_iterator 631 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 632 POINTER :: sab_orb 633 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a, dftb_kind_b 634 TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), & 635 POINTER :: dftb_potential 636 TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji 637 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 638 639 CALL timeset(routineN, handle) 640 641 CALL get_qs_env(qs_env=qs_env, nkind=nkind) 642 CPASSERT(.NOT. ASSOCIATED(sap_int)) 643 ALLOCATE (sap_int(nkind*nkind)) 644 DO i = 1, nkind*nkind 645 NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex) 646 sap_int(i)%nalist = 0 647 END DO 648 649 CALL get_qs_env(qs_env=qs_env, & 650 qs_kind_set=qs_kind_set, & 651 dft_control=dft_control, & 652 sab_orb=sab_orb) 653 654 dftb_control => dft_control%qs_control%dftb_control 655 656 NULLIFY (dftb_potential) 657 CALL get_qs_env(qs_env=qs_env, & 658 dftb_potential=dftb_potential) 659 660 ! loop over all atom pairs with a non-zero overlap (sab_orb) 661 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 662 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 663 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, & 664 jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, & 665 inode=jneighbor, cell=cell, r=rij) 666 iac = ikind + nkind*(jkind - 1) 667 ! 668 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a) 669 CALL get_dftb_atom_param(dftb_kind_a, & 670 defined=defined, lmax=lmaxi, natorb=natorb_a) 671 IF (.NOT. defined .OR. natorb_a < 1) CYCLE 672 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b) 673 CALL get_dftb_atom_param(dftb_kind_b, & 674 defined=defined, lmax=lmaxj, natorb=natorb_b) 675 IF (.NOT. defined .OR. natorb_b < 1) CYCLE 676 677 dr = SQRT(SUM(rij(:)**2)) 678 679 ! integral list 680 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN 681 sap_int(iac)%a_kind = ikind 682 sap_int(iac)%p_kind = jkind 683 sap_int(iac)%nalist = nlist 684 ALLOCATE (sap_int(iac)%alist(nlist)) 685 DO i = 1, nlist 686 NULLIFY (sap_int(iac)%alist(i)%clist) 687 sap_int(iac)%alist(i)%aatom = 0 688 sap_int(iac)%alist(i)%nclist = 0 689 END DO 690 END IF 691 IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN 692 sap_int(iac)%alist(ilist)%aatom = iatom 693 sap_int(iac)%alist(ilist)%nclist = nneighbor 694 ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor)) 695 DO i = 1, nneighbor 696 sap_int(iac)%alist(ilist)%clist(i)%catom = 0 697 END DO 698 END IF 699 clist => sap_int(iac)%alist(ilist)%clist(jneighbor) 700 clist%catom = jatom 701 clist%cell = cell 702 clist%rac = rij 703 ALLOCATE (clist%acint(natorb_a, natorb_b, 3)) 704 NULLIFY (clist%achint) 705 clist%acint = 0._dp 706 clist%nsgf_cnt = 0 707 NULLIFY (clist%sgf_list) 708 709 ! retrieve information on S matrix 710 dftb_param_ij => dftb_potential(ikind, jkind) 711 dftb_param_ji => dftb_potential(jkind, ikind) 712 ! assume table size and type is symmetric 713 ngrd = dftb_param_ij%ngrd 714 ngrdcut = dftb_param_ij%ngrdcut 715 dgrd = dftb_param_ij%dgrd 716 ddr = dgrd*0.1_dp 717 CPASSERT(dftb_param_ij%llm == dftb_param_ji%llm) 718 llm = dftb_param_ij%llm 719 smatij => dftb_param_ij%smat 720 smatji => dftb_param_ji%smat 721 722 IF (NINT(dr/dgrd) <= ngrdcut) THEN 723 IF (iatom == jatom .AND. dr < 0.001_dp) THEN 724 ! diagonal block 725 ELSE 726 ! off-diagonal block 727 n1 = natorb_a 728 n2 = natorb_b 729 ALLOCATE (dsblock(n1, n2), dsblockm(n1, n2)) 730 DO i = 1, 3 731 dsblock = 0._dp 732 dsblockm = 0.0_dp 733 drij = rij 734 735 drij(i) = rij(i) - ddr 736 CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, & 737 llm, lmaxi, lmaxj, iatom, iatom) 738 739 drij(i) = rij(i) + ddr 740 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, & 741 llm, lmaxi, lmaxj, iatom, iatom) 742 743 dsblock = dsblock - dsblockm 744 dsblock = dsblock/(2.0_dp*ddr) 745 746 clist%acint(1:n1, 1:n2, i) = -dsblock(1:n1, 1:n2) 747 ENDDO 748 DEALLOCATE (dsblock, dsblockm) 749 END IF 750 END IF 751 752 END DO 753 CALL neighbor_list_iterator_release(nl_iterator) 754 755 CALL timestop(handle) 756 757 END SUBROUTINE dftb_dsint_list 758 759END MODULE qs_dftb_coulomb 760 761