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 Hessian contributions in xTB 8!> \author JGH 9! ************************************************************************************************** 10MODULE xtb_ehess 11 USE atomic_kind_types, ONLY: atomic_kind_type,& 12 get_atomic_kind_set 13 USE atprop_types, ONLY: atprop_type 14 USE cell_types, ONLY: cell_type,& 15 get_cell,& 16 pbc 17 USE cp_control_types, ONLY: dft_control_type 18 USE cp_para_types, ONLY: cp_para_env_type 19 USE dbcsr_api, ONLY: dbcsr_get_block_p,& 20 dbcsr_iterator_blocks_left,& 21 dbcsr_iterator_next_block,& 22 dbcsr_iterator_start,& 23 dbcsr_iterator_stop,& 24 dbcsr_iterator_type,& 25 dbcsr_p_type 26 USE distribution_1d_types, ONLY: distribution_1d_type 27 USE ewald_environment_types, ONLY: ewald_env_get,& 28 ewald_environment_type 29 USE ewald_methods_tb, ONLY: tb_ewald_overlap,& 30 tb_spme_evaluate 31 USE ewald_pw_types, ONLY: ewald_pw_type 32 USE kinds, ONLY: dp 33 USE mathconstants, ONLY: oorootpi,& 34 pi 35 USE message_passing, ONLY: mp_sum 36 USE particle_types, ONLY: particle_type 37 USE pw_poisson_types, ONLY: do_ewald_ewald,& 38 do_ewald_none,& 39 do_ewald_pme,& 40 do_ewald_spme 41 USE qs_environment_types, ONLY: get_qs_env,& 42 qs_environment_type 43 USE qs_kind_types, ONLY: get_qs_kind,& 44 qs_kind_type 45 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 46 neighbor_list_iterate,& 47 neighbor_list_iterator_create,& 48 neighbor_list_iterator_p_type,& 49 neighbor_list_iterator_release,& 50 neighbor_list_set_p_type 51 USE virial_types, ONLY: virial_type 52 USE xtb_coulomb, ONLY: gamma_rab_sr 53 USE xtb_types, ONLY: get_xtb_atom_param,& 54 xtb_atom_type 55#include "./base/base_uses.f90" 56 57 IMPLICIT NONE 58 59 PRIVATE 60 61 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_ehess' 62 63 PUBLIC :: xtb_coulomb_hessian 64 65CONTAINS 66 67! ************************************************************************************************** 68!> \brief ... 69!> \param qs_env ... 70!> \param ks_matrix ... 71!> \param charges1 ... 72!> \param mcharge1 ... 73!> \param mcharge ... 74! ************************************************************************************************** 75 SUBROUTINE xtb_coulomb_hessian(qs_env, ks_matrix, charges1, mcharge1, mcharge) 76 77 TYPE(qs_environment_type), POINTER :: qs_env 78 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix 79 REAL(dp), DIMENSION(:, :) :: charges1 80 REAL(dp), DIMENSION(:) :: mcharge1, mcharge 81 82 CHARACTER(len=*), PARAMETER :: routineN = 'xtb_coulomb_hessian', & 83 routineP = moduleN//':'//routineN 84 85 INTEGER :: blk, ewald_type, handle, i, ia, iatom, icol, ikind, irow, is, j, jatom, jkind, & 86 la, lb, lmaxa, lmaxb, natom, natorb_a, natorb_b, ni, nj, nkind, nmat, za, zb 87 INTEGER, DIMENSION(25) :: laoa, laob 88 INTEGER, DIMENSION(3) :: cellind, periodic 89 INTEGER, DIMENSION(:), POINTER :: kind_of 90 LOGICAL :: defined, do_ewald, found 91 REAL(KIND=dp) :: alpha, deth, dr, etaa, etab, gmij, kg, & 92 rcut, rcuta, rcutb 93 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xgamma 94 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gammab, gcij, gmcharge 95 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gchrg 96 REAL(KIND=dp), DIMENSION(3) :: rij 97 REAL(KIND=dp), DIMENSION(5) :: kappaa, kappab 98 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksblock, sblock 99 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 100 TYPE(atprop_type), POINTER :: atprop 101 TYPE(cell_type), POINTER :: cell 102 TYPE(cp_para_env_type), POINTER :: para_env 103 TYPE(dbcsr_iterator_type) :: iter 104 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s 105 TYPE(dft_control_type), POINTER :: dft_control 106 TYPE(distribution_1d_type), POINTER :: local_particles 107 TYPE(ewald_environment_type), POINTER :: ewald_env 108 TYPE(ewald_pw_type), POINTER :: ewald_pw 109 TYPE(neighbor_list_iterator_p_type), & 110 DIMENSION(:), POINTER :: nl_iterator 111 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 112 POINTER :: n_list 113 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 114 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 115 TYPE(virial_type), POINTER :: virial 116 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind 117 118 CALL timeset(routineN, handle) 119 120 CALL get_qs_env(qs_env, & 121 matrix_s_kp=matrix_s, & 122 qs_kind_set=qs_kind_set, & 123 particle_set=particle_set, & 124 cell=cell, & 125 dft_control=dft_control) 126 127 IF (dft_control%nimages /= 1) THEN 128 CPABORT("No kpoints allowed in xTB response calculation") 129 END IF 130 131 CALL get_qs_env(qs_env, nkind=nkind, natom=natom) 132 nmat = 1 133 ALLOCATE (gchrg(natom, 5, nmat)) 134 gchrg = 0._dp 135 ALLOCATE (gmcharge(natom, nmat)) 136 gmcharge = 0._dp 137 138 ! short range contribution (gamma) 139 ! loop over all atom pairs with a non-zero overlap (sab_orb) 140 kg = dft_control%qs_control%xtb_control%kg 141 NULLIFY (n_list) 142 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) 143 CALL neighbor_list_iterator_create(nl_iterator, n_list) 144 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 145 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 146 iatom=iatom, jatom=jatom, r=rij, cell=cellind) 147 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a) 148 CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a) 149 IF (.NOT. defined .OR. natorb_a < 1) CYCLE 150 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b) 151 CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b) 152 IF (.NOT. defined .OR. natorb_b < 1) CYCLE 153 ! atomic parameters 154 CALL get_xtb_atom_param(xtb_atom_a, eta=etaa, lmax=lmaxa, kappa=kappaa, rcut=rcuta) 155 CALL get_xtb_atom_param(xtb_atom_b, eta=etab, lmax=lmaxb, kappa=kappab, rcut=rcutb) 156 ! gamma matrix 157 ni = lmaxa + 1 158 nj = lmaxb + 1 159 ALLOCATE (gammab(ni, nj)) 160 rcut = rcuta + rcutb 161 dr = SQRT(SUM(rij(:)**2)) 162 CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut) 163 gchrg(iatom, 1:ni, 1) = gchrg(iatom, 1:ni, 1) + MATMUL(gammab, charges1(jatom, 1:nj)) 164 IF (iatom /= jatom) THEN 165 gchrg(jatom, 1:nj, 1) = gchrg(jatom, 1:nj, 1) + MATMUL(charges1(iatom, 1:ni), gammab) 166 END IF 167 DEALLOCATE (gammab) 168 END DO 169 CALL neighbor_list_iterator_release(nl_iterator) 170 171 ! 1/R contribution 172 173 do_ewald = dft_control%qs_control%xtb_control%do_ewald 174 IF (do_ewald) THEN 175 ! Ewald sum 176 NULLIFY (ewald_env, ewald_pw) 177 NULLIFY (virial, atprop) 178 CALL get_qs_env(qs_env=qs_env, & 179 ewald_env=ewald_env, ewald_pw=ewald_pw) 180 CALL get_cell(cell=cell, periodic=periodic, deth=deth) 181 CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type) 182 CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list) 183 CALL tb_ewald_overlap(gmcharge, mcharge1, alpha, n_list, virial, .FALSE., atprop) 184 SELECT CASE (ewald_type) 185 CASE DEFAULT 186 CPABORT("Invalid Ewald type") 187 CASE (do_ewald_none) 188 CPABORT("Not allowed with DFTB") 189 CASE (do_ewald_ewald) 190 CPABORT("Standard Ewald not implemented in DFTB") 191 CASE (do_ewald_pme) 192 CPABORT("PME not implemented in DFTB") 193 CASE (do_ewald_spme) 194 CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, & 195 gmcharge, mcharge1, .FALSE., virial, .FALSE., atprop) 196 END SELECT 197 ELSE 198 ! direct sum 199 CALL get_qs_env(qs_env=qs_env, & 200 local_particles=local_particles) 201 DO ikind = 1, SIZE(local_particles%n_el) 202 DO ia = 1, local_particles%n_el(ikind) 203 iatom = local_particles%list(ikind)%array(ia) 204 DO jatom = 1, iatom - 1 205 rij = particle_set(iatom)%r - particle_set(jatom)%r 206 rij = pbc(rij, cell) 207 dr = SQRT(SUM(rij(:)**2)) 208 IF (dr > 1.e-6_dp) THEN 209 gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge1(jatom)/dr 210 gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge1(iatom)/dr 211 END IF 212 END DO 213 END DO 214 END DO 215 END IF 216 217 ! global sum of gamma*p arrays 218 CALL get_qs_env(qs_env=qs_env, para_env=para_env) 219 CALL mp_sum(gmcharge(:, 1), para_env%group) 220 CALL mp_sum(gchrg(:, :, 1), para_env%group) 221 222 IF (do_ewald) THEN 223 ! add self charge interaction and background charge contribution 224 gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:) 225 IF (ANY(periodic(:) == 1)) THEN 226 gmcharge(:, 1) = gmcharge(:, 1) - pi/alpha**2/deth 227 END IF 228 END IF 229 230 ALLOCATE (kind_of(natom)) 231 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set) 232 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of) 233 234 ! no k-points; all matrices have been transformed to periodic bsf 235 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix) 236 DO WHILE (dbcsr_iterator_blocks_left(iter)) 237 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk) 238 ikind = kind_of(irow) 239 jkind = kind_of(icol) 240 241 ! atomic parameters 242 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a) 243 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b) 244 CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa) 245 CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob) 246 247 ni = SIZE(sblock, 1) 248 nj = SIZE(sblock, 2) 249 ALLOCATE (gcij(ni, nj)) 250 DO i = 1, ni 251 DO j = 1, nj 252 la = laoa(i) + 1 253 lb = laob(j) + 1 254 gcij(i, j) = 0.5_dp*(gchrg(irow, la, 1) + gchrg(icol, lb, 1)) 255 END DO 256 END DO 257 gmij = 0.5_dp*(gmcharge(irow, 1) + gmcharge(icol, 1)) 258 DO is = 1, SIZE(ks_matrix) 259 NULLIFY (ksblock) 260 CALL dbcsr_get_block_p(matrix=ks_matrix(is)%matrix, & 261 row=irow, col=icol, block=ksblock, found=found) 262 CPASSERT(found) 263 ksblock = ksblock - gcij*sblock 264 ksblock = ksblock - gmij*sblock 265 END DO 266 DEALLOCATE (gcij) 267 ENDDO 268 CALL dbcsr_iterator_stop(iter) 269 270 IF (dft_control%qs_control%xtb_control%tb3_interaction) THEN 271 CALL get_qs_env(qs_env, nkind=nkind) 272 ALLOCATE (xgamma(nkind)) 273 DO ikind = 1, nkind 274 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind) 275 CALL get_xtb_atom_param(xtb_kind, xgamma=xgamma(ikind)) 276 END DO 277 ! Diagonal 3rd order correction (DFTB3) 278 CALL dftb3_diagonal_hessian(qs_env, ks_matrix, mcharge, mcharge1, xgamma) 279 DEALLOCATE (xgamma) 280 END IF 281 282 IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN 283 CPABORT("QMMM not available in xTB response calculations") 284 END IF 285 286 DEALLOCATE (gmcharge, gchrg, kind_of) 287 288 CALL timestop(handle) 289 290 END SUBROUTINE xtb_coulomb_hessian 291 292! ************************************************************************************************** 293!> \brief ... 294!> \param qs_env ... 295!> \param ks_matrix ... 296!> \param mcharge ... 297!> \param mcharge1 ... 298!> \param xgamma ... 299! ************************************************************************************************** 300 SUBROUTINE dftb3_diagonal_hessian(qs_env, ks_matrix, mcharge, mcharge1, xgamma) 301 302 TYPE(qs_environment_type), POINTER :: qs_env 303 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix 304 REAL(dp), DIMENSION(:) :: mcharge, mcharge1, xgamma 305 306 CHARACTER(len=*), PARAMETER :: routineN = 'dftb3_diagonal_hessian', & 307 routineP = moduleN//':'//routineN 308 309 INTEGER :: blk, handle, icol, ikind, irow, is, & 310 jkind, natom 311 INTEGER, DIMENSION(:), POINTER :: kind_of 312 LOGICAL :: found 313 REAL(KIND=dp) :: gmij, ui, uj 314 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksblock, sblock 315 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 316 TYPE(dbcsr_iterator_type) :: iter 317 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s 318 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 319 320 CALL timeset(routineN, handle) 321 322 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s, natom=natom) 323 ALLOCATE (kind_of(natom)) 324 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set) 325 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of) 326 ! no k-points; all matrices have been transformed to periodic bsf 327 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix) 328 DO WHILE (dbcsr_iterator_blocks_left(iter)) 329 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk) 330 ikind = kind_of(irow) 331 ui = xgamma(ikind) 332 jkind = kind_of(icol) 333 uj = xgamma(jkind) 334 gmij = ui*mcharge(irow)*mcharge1(irow) + uj*mcharge(icol)*mcharge1(icol) 335 DO is = 1, SIZE(ks_matrix) 336 NULLIFY (ksblock) 337 CALL dbcsr_get_block_p(matrix=ks_matrix(is)%matrix, & 338 row=irow, col=icol, block=ksblock, found=found) 339 CPASSERT(found) 340 ksblock = ksblock + 0.5_dp*gmij*sblock 341 END DO 342 ENDDO 343 CALL dbcsr_iterator_stop(iter) 344 DEALLOCATE (kind_of) 345 346 CALL timestop(handle) 347 348 END SUBROUTINE dftb3_diagonal_hessian 349 350END MODULE xtb_ehess 351 352