1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of DFTB3 Terms 8!> \author JGH 9! ************************************************************************************************** 10MODULE qs_dftb3_methods 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 USE cp_para_types, ONLY: cp_para_env_type 16 USE dbcsr_api, ONLY: dbcsr_add,& 17 dbcsr_get_block_p,& 18 dbcsr_iterator_blocks_left,& 19 dbcsr_iterator_next_block,& 20 dbcsr_iterator_start,& 21 dbcsr_iterator_stop,& 22 dbcsr_iterator_type,& 23 dbcsr_p_type 24 USE distribution_1d_types, ONLY: distribution_1d_type 25 USE kinds, ONLY: dp 26 USE kpoint_types, ONLY: get_kpoint_info,& 27 kpoint_type 28 USE message_passing, ONLY: mp_sum 29 USE particle_types, ONLY: particle_type 30 USE qs_energy_types, ONLY: qs_energy_type 31 USE qs_environment_types, ONLY: get_qs_env,& 32 qs_environment_type 33 USE qs_force_types, ONLY: qs_force_type 34 USE qs_kind_types, ONLY: qs_kind_type 35 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 36 neighbor_list_iterate,& 37 neighbor_list_iterator_create,& 38 neighbor_list_iterator_p_type,& 39 neighbor_list_iterator_release,& 40 neighbor_list_set_p_type 41 USE qs_rho_types, ONLY: qs_rho_get,& 42 qs_rho_type 43 USE sap_kind_types, ONLY: sap_int_type 44 USE virial_methods, ONLY: virial_pair_force 45 USE virial_types, ONLY: virial_type 46#include "./base/base_uses.f90" 47 48 IMPLICIT NONE 49 50 PRIVATE 51 52 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb3_methods' 53 54 PUBLIC :: build_dftb3_diagonal 55 56CONTAINS 57 58! ************************************************************************************************** 59!> \brief ... 60!> \param qs_env ... 61!> \param ks_matrix ... 62!> \param rho ... 63!> \param mcharge ... 64!> \param energy ... 65!> \param xgamma ... 66!> \param zeff ... 67!> \param sap_int ... 68!> \param calculate_forces ... 69!> \param just_energy ... 70! ************************************************************************************************** 71 SUBROUTINE build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeff, & 72 sap_int, calculate_forces, just_energy) 73 74 TYPE(qs_environment_type), POINTER :: qs_env 75 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix 76 TYPE(qs_rho_type), POINTER :: rho 77 REAL(dp), DIMENSION(:), INTENT(in) :: mcharge 78 TYPE(qs_energy_type), POINTER :: energy 79 REAL(dp), DIMENSION(:), INTENT(in) :: xgamma, zeff 80 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int 81 LOGICAL, INTENT(in) :: calculate_forces, just_energy 82 83 CHARACTER(len=*), PARAMETER :: routineN = 'build_dftb3_diagonal', & 84 routineP = moduleN//':'//routineN 85 86 INTEGER :: atom_i, atom_j, blk, handle, i, ia, iac, & 87 iatom, ic, icol, ikind, irow, is, & 88 jatom, jkind, natom, nimg, nkind 89 INTEGER, DIMENSION(3) :: cellind 90 INTEGER, DIMENSION(:), POINTER :: atom_of_kind, kind_of 91 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 92 LOGICAL :: found, use_virial 93 REAL(KIND=dp) :: dr, eb3, eloc, fi, gmij, ua, ui, uj 94 REAL(KIND=dp), DIMENSION(3) :: fij, rij 95 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, ksblock, pblock, sblock 96 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint 97 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 98 TYPE(atprop_type), POINTER :: atprop 99 TYPE(cell_type), POINTER :: cell 100 TYPE(cp_para_env_type), POINTER :: para_env 101 TYPE(dbcsr_iterator_type) :: iter 102 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s 103 TYPE(distribution_1d_type), POINTER :: local_particles 104 TYPE(kpoint_type), POINTER :: kpoints 105 TYPE(neighbor_list_iterator_p_type), & 106 DIMENSION(:), POINTER :: nl_iterator 107 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 108 POINTER :: n_list 109 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 110 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 111 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 112 TYPE(virial_type), POINTER :: virial 113 114 CALL timeset(routineN, handle) 115 NULLIFY (atprop) 116 117 ! Energy 118 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, & 119 qs_kind_set=qs_kind_set, atprop=atprop) 120 121 eb3 = 0.0_dp 122 CALL get_qs_env(qs_env=qs_env, local_particles=local_particles) 123 DO ikind = 1, SIZE(local_particles%n_el) 124 ua = xgamma(ikind) 125 DO ia = 1, local_particles%n_el(ikind) 126 iatom = local_particles%list(ikind)%array(ia) 127 eloc = -1.0_dp/6.0_dp*ua*mcharge(iatom)**3 128 eb3 = eb3 + eloc 129 IF (atprop%energy) THEN 130 ! we have to add the part not covered by 0.5*Tr(FP) 131 eloc = -0.5_dp*eloc - 0.25_dp*ua*zeff(ikind)*mcharge(iatom)**2 132 atprop%atecoul(iatom) = atprop%atecoul(iatom) + eloc 133 END IF 134 END DO 135 END DO 136 CALL get_qs_env(qs_env=qs_env, para_env=para_env) 137 CALL mp_sum(eb3, para_env%group) 138 energy%dftb3 = eb3 139 140 ! Forces and Virial 141 IF (calculate_forces) THEN 142 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s, natom=natom, force=force, & 143 cell=cell, virial=virial, particle_set=particle_set) 144 ! virial 145 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 146 147 ALLOCATE (atom_of_kind(natom), kind_of(natom)) 148 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 149 kind_of=kind_of, atom_of_kind=atom_of_kind) 150 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 151 IF (SIZE(matrix_p, 1) == 2) THEN 152 DO ic = 1, SIZE(matrix_p, 2) 153 CALL dbcsr_add(matrix_p(1, ic)%matrix, matrix_p(2, ic)%matrix, & 154 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 155 END DO 156 END IF 157 ! 158 nimg = SIZE(matrix_p, 2) 159 NULLIFY (cell_to_index) 160 IF (nimg > 1) THEN 161 NULLIFY (kpoints) 162 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints) 163 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 164 END IF 165 IF (nimg == 1) THEN 166 ! no k-points; all matrices have been transformed to periodic bsf 167 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix) 168 DO WHILE (dbcsr_iterator_blocks_left(iter)) 169 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk) 170 ikind = kind_of(irow) 171 atom_i = atom_of_kind(irow) 172 ui = xgamma(ikind) 173 jkind = kind_of(icol) 174 atom_j = atom_of_kind(icol) 175 uj = xgamma(jkind) 176 ! 177 gmij = -0.5_dp*(ui*mcharge(irow)**2 + uj*mcharge(icol)**2) 178 ! 179 NULLIFY (pblock) 180 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, & 181 row=irow, col=icol, block=pblock, found=found) 182 CPASSERT(found) 183 DO i = 1, 3 184 NULLIFY (dsblock) 185 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, & 186 row=irow, col=icol, block=dsblock, found=found) 187 CPASSERT(found) 188 fi = -gmij*SUM(pblock*dsblock) 189 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi 190 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi 191 fij(i) = fi 192 END DO 193 ENDDO 194 CALL dbcsr_iterator_stop(iter) 195 ! use dsint list 196 IF (use_virial) THEN 197 CPASSERT(ASSOCIATED(sap_int)) 198 CALL get_qs_env(qs_env, nkind=nkind) 199 DO ikind = 1, nkind 200 DO jkind = 1, nkind 201 iac = ikind + nkind*(jkind - 1) 202 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE 203 ui = xgamma(ikind) 204 uj = xgamma(jkind) 205 DO ia = 1, sap_int(iac)%nalist 206 IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE 207 iatom = sap_int(iac)%alist(ia)%aatom 208 DO ic = 1, sap_int(iac)%alist(ia)%nclist 209 jatom = sap_int(iac)%alist(ia)%clist(ic)%catom 210 rij = sap_int(iac)%alist(ia)%clist(ic)%rac 211 dr = SQRT(SUM(rij(:)**2)) 212 IF (dr > 1.e-6_dp) THEN 213 dsint => sap_int(iac)%alist(ia)%clist(ic)%acint 214 gmij = -0.5_dp*(ui*mcharge(iatom)**2 + uj*mcharge(jatom)**2) 215 icol = MAX(iatom, jatom) 216 irow = MIN(iatom, jatom) 217 NULLIFY (pblock) 218 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, & 219 row=irow, col=icol, block=pblock, found=found) 220 CPASSERT(found) 221 DO i = 1, 3 222 IF (irow == iatom) THEN 223 fij(i) = -gmij*SUM(pblock*dsint(:, :, i)) 224 ELSE 225 fij(i) = -gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i)) 226 END IF 227 END DO 228 fi = 1.0_dp 229 IF (iatom == jatom) fi = 0.5_dp 230 CALL virial_pair_force(virial%pv_virial, fi, fij, rij) 231 IF (atprop%stress) THEN 232 CALL virial_pair_force(atprop%atstress(:, :, irow), fi*0.5_dp, fij, rij) 233 CALL virial_pair_force(atprop%atstress(:, :, icol), fi*0.5_dp, fij, rij) 234 END IF 235 END IF 236 END DO 237 END DO 238 END DO 239 END DO 240 END IF 241 ELSE 242 NULLIFY (n_list) 243 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) 244 CALL neighbor_list_iterator_create(nl_iterator, n_list) 245 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 246 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 247 iatom=iatom, jatom=jatom, r=rij, cell=cellind) 248 249 dr = SQRT(SUM(rij**2)) 250 IF (iatom == jatom .AND. dr < 1.0e-6_dp) CYCLE 251 252 icol = MAX(iatom, jatom) 253 irow = MIN(iatom, jatom) 254 255 ic = cell_to_index(cellind(1), cellind(2), cellind(3)) 256 CPASSERT(ic > 0) 257 258 ikind = kind_of(iatom) 259 atom_i = atom_of_kind(iatom) 260 ui = xgamma(ikind) 261 jkind = kind_of(jatom) 262 atom_j = atom_of_kind(jatom) 263 uj = xgamma(jkind) 264 ! 265 gmij = -0.5_dp*(ui*mcharge(iatom)**2 + uj*mcharge(jatom)**2) 266 ! 267 NULLIFY (pblock) 268 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, & 269 row=irow, col=icol, block=pblock, found=found) 270 CPASSERT(found) 271 DO i = 1, 3 272 NULLIFY (dsblock) 273 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, & 274 row=irow, col=icol, block=dsblock, found=found) 275 CPASSERT(found) 276 IF (irow == iatom) THEN 277 fi = -gmij*SUM(pblock*dsblock) 278 ELSE 279 fi = gmij*SUM(pblock*dsblock) 280 END IF 281 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi 282 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi 283 fij(i) = fi 284 END DO 285 IF (use_virial) THEN 286 fi = 1.0_dp 287 IF (iatom == jatom) fi = 0.5_dp 288 CALL virial_pair_force(virial%pv_virial, fi, fij, rij) 289 IF (atprop%stress) THEN 290 CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij) 291 CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij) 292 END IF 293 END IF 294 295 END DO 296 CALL neighbor_list_iterator_release(nl_iterator) 297 ! 298 END IF 299 300 DEALLOCATE (atom_of_kind, kind_of) 301 IF (SIZE(matrix_p, 1) == 2) THEN 302 DO ic = 1, SIZE(matrix_p, 2) 303 CALL dbcsr_add(matrix_p(1, ic)%matrix, matrix_p(2, ic)%matrix, & 304 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp) 305 END DO 306 END IF 307 END IF 308 309 ! KS matrix 310 IF (.NOT. just_energy) THEN 311 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s, natom=natom) 312 ALLOCATE (kind_of(natom)) 313 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of) 314 ! 315 nimg = SIZE(ks_matrix, 2) 316 NULLIFY (cell_to_index) 317 IF (nimg > 1) THEN 318 NULLIFY (kpoints) 319 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints) 320 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 321 END IF 322 323 IF (nimg == 1) THEN 324 ! no k-points; all matrices have been transformed to periodic bsf 325 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix) 326 DO WHILE (dbcsr_iterator_blocks_left(iter)) 327 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk) 328 ikind = kind_of(irow) 329 ui = xgamma(ikind) 330 jkind = kind_of(icol) 331 uj = xgamma(jkind) 332 gmij = -0.5_dp*(ui*mcharge(irow)**2 + uj*mcharge(icol)**2) 333 DO is = 1, SIZE(ks_matrix, 1) 334 NULLIFY (ksblock) 335 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, & 336 row=irow, col=icol, block=ksblock, found=found) 337 CPASSERT(found) 338 ksblock = ksblock - 0.5_dp*gmij*sblock 339 END DO 340 ENDDO 341 CALL dbcsr_iterator_stop(iter) 342 ELSE 343 NULLIFY (n_list) 344 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) 345 CALL neighbor_list_iterator_create(nl_iterator, n_list) 346 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 347 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 348 iatom=iatom, jatom=jatom, r=rij, cell=cellind) 349 350 icol = MAX(iatom, jatom) 351 irow = MIN(iatom, jatom) 352 353 ic = cell_to_index(cellind(1), cellind(2), cellind(3)) 354 CPASSERT(ic > 0) 355 356 ikind = kind_of(iatom) 357 ui = xgamma(ikind) 358 jkind = kind_of(jatom) 359 uj = xgamma(jkind) 360 gmij = -0.5_dp*(ui*mcharge(iatom)**2 + uj*mcharge(jatom)**2) 361 ! 362 NULLIFY (sblock) 363 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, & 364 row=irow, col=icol, block=sblock, found=found) 365 CPASSERT(found) 366 DO is = 1, SIZE(ks_matrix, 1) 367 NULLIFY (ksblock) 368 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, & 369 row=irow, col=icol, block=ksblock, found=found) 370 CPASSERT(found) 371 ksblock = ksblock - 0.5_dp*gmij*sblock 372 END DO 373 374 END DO 375 CALL neighbor_list_iterator_release(nl_iterator) 376 ! 377 END IF 378 DEALLOCATE (kind_of) 379 END IF 380 381 CALL timestop(handle) 382 383 END SUBROUTINE build_dftb3_diagonal 384 385END MODULE qs_dftb3_methods 386 387