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 electric field contributions in TB 8!> \author JGH 9! ************************************************************************************************** 10MODULE efield_tb_methods 11 USE atomic_kind_types, ONLY: atomic_kind_type,& 12 get_atomic_kind_set 13 USE cell_types, ONLY: cell_type,& 14 pbc 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_para_types, ONLY: cp_para_env_type 17 USE dbcsr_api, ONLY: 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 kinds, ONLY: dp 25 USE kpoint_types, ONLY: get_kpoint_info,& 26 kpoint_type 27 USE mathconstants, ONLY: pi,& 28 twopi 29 USE message_passing, ONLY: mp_sum 30 USE particle_types, ONLY: particle_type 31 USE qs_energy_types, ONLY: qs_energy_type 32 USE qs_environment_types, ONLY: get_qs_env,& 33 qs_environment_type,& 34 set_qs_env 35 USE qs_force_types, ONLY: qs_force_type 36 USE qs_kind_types, ONLY: qs_kind_type 37 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 38 neighbor_list_iterate,& 39 neighbor_list_iterator_create,& 40 neighbor_list_iterator_p_type,& 41 neighbor_list_iterator_release,& 42 neighbor_list_set_p_type 43 USE qs_period_efield_types, ONLY: efield_berry_type,& 44 init_efield_matrices 45 USE qs_rho_types, ONLY: qs_rho_get,& 46 qs_rho_type 47 USE virial_methods, ONLY: virial_pair_force 48 USE virial_types, ONLY: virial_type 49#include "./base/base_uses.f90" 50 51 IMPLICIT NONE 52 53 PRIVATE 54 55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'efield_tb_methods' 56 57 PUBLIC :: efield_tb_matrix 58 59CONTAINS 60 61! ************************************************************************************************** 62!> \brief ... 63!> \param qs_env ... 64!> \param ks_matrix ... 65!> \param rho ... 66!> \param mcharge ... 67!> \param energy ... 68!> \param calculate_forces ... 69!> \param just_energy ... 70! ************************************************************************************************** 71 SUBROUTINE efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) 72 73 TYPE(qs_environment_type), POINTER :: qs_env 74 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix 75 TYPE(qs_rho_type), POINTER :: rho 76 REAL(dp), DIMENSION(:), INTENT(in) :: mcharge 77 TYPE(qs_energy_type), POINTER :: energy 78 LOGICAL, INTENT(in) :: calculate_forces, just_energy 79 80 CHARACTER(len=*), PARAMETER :: routineN = 'efield_tb_matrix', & 81 routineP = moduleN//':'//routineN 82 83 INTEGER :: handle 84 TYPE(dft_control_type), POINTER :: dft_control 85 86 CALL timeset(routineN, handle) 87 88 energy%efield = 0.0_dp 89 CALL get_qs_env(qs_env, dft_control=dft_control) 90 IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN 91 IF (dft_control%apply_period_efield) THEN 92 IF (dft_control%period_efield%displacement_field) THEN 93 CALL dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) 94 ELSE 95 CALL efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) 96 END IF 97 ELSE IF (dft_control%apply_efield) THEN 98 CALL efield_tb_local(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) 99 ELSE IF (dft_control%apply_efield_field) THEN 100 CPABORT("efield_filed") 101 END IF 102 ELSE 103 CPABORT("This routine should only be called from TB") 104 END IF 105 106 CALL timestop(handle) 107 108 END SUBROUTINE efield_tb_matrix 109 110! ************************************************************************************************** 111!> \brief ... 112!> \param qs_env ... 113!> \param ks_matrix ... 114!> \param rho ... 115!> \param mcharge ... 116!> \param energy ... 117!> \param calculate_forces ... 118!> \param just_energy ... 119! ************************************************************************************************** 120 SUBROUTINE efield_tb_local(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) 121 TYPE(qs_environment_type), POINTER :: qs_env 122 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix 123 TYPE(qs_rho_type), POINTER :: rho 124 REAL(dp), DIMENSION(:), INTENT(in) :: mcharge 125 TYPE(qs_energy_type), POINTER :: energy 126 LOGICAL, INTENT(in) :: calculate_forces, just_energy 127 128 CHARACTER(LEN=*), PARAMETER :: routineN = 'efield_tb_local', & 129 routineP = moduleN//':'//routineN 130 131 INTEGER :: atom_a, atom_b, blk, handle, ia, icol, & 132 idir, ikind, irow, ispin, jkind, & 133 natom, nspin 134 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of 135 LOGICAL :: do_kpoints, found, use_virial 136 REAL(dp) :: charge, fdir 137 REAL(dp), DIMENSION(3) :: ci, fieldpol, fij, ria, rib 138 REAL(dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block 139 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 140 TYPE(cell_type), POINTER :: cell 141 TYPE(cp_para_env_type), POINTER :: para_env 142 TYPE(dbcsr_iterator_type) :: iter 143 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s 144 TYPE(dft_control_type), POINTER :: dft_control 145 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 146 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 147 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 148 TYPE(virial_type), POINTER :: virial 149 150 CALL timeset(routineN, handle) 151 152 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, particle_set=particle_set) 153 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, energy=energy, para_env=para_env) 154 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints, virial=virial) 155 IF (do_kpoints) THEN 156 CPABORT("Local electric field with kpoints not possible. Use Berry phase periodic version") 157 END IF 158 ! disable stress calculation 159 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 160 IF (use_virial) THEN 161 CPABORT("Stress tensor for non-periodic E-field not possible") 162 END IF 163 164 fieldpol = dft_control%efield_fields(1)%efield%polarisation* & 165 dft_control%efield_fields(1)%efield%strength 166 167 natom = SIZE(particle_set) 168 ci = 0.0_dp 169 DO ia = 1, natom 170 charge = mcharge(ia) 171 ria = particle_set(ia)%r 172 ria = pbc(ria, cell) 173 ci(:) = ci(:) + charge*ria(:) 174 END DO 175 energy%efield = -SUM(ci(:)*fieldpol(:)) 176 177 IF (.NOT. just_energy) THEN 178 179 IF (calculate_forces) THEN 180 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force) 181 ALLOCATE (atom_of_kind(natom), kind_of(natom)) 182 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of) 183 IF (para_env%mepos == 0) THEN 184 DO ia = 1, natom 185 charge = mcharge(ia) 186 ikind = kind_of(ia) 187 atom_a = atom_of_kind(ia) 188 force(ikind)%efield(1:3, atom_a) = -charge*fieldpol(:) 189 END DO 190 ELSE 191 DO ia = 1, natom 192 ikind = kind_of(ia) 193 atom_a = atom_of_kind(ia) 194 force(ikind)%efield(1:3, atom_a) = 0.0_dp 195 END DO 196 END IF 197 CALL qs_rho_get(rho, rho_ao=matrix_p) 198 END IF 199 200 ! Update KS matrix 201 nspin = SIZE(ks_matrix, 1) 202 NULLIFY (matrix_s) 203 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s) 204 CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix) 205 DO WHILE (dbcsr_iterator_blocks_left(iter)) 206 NULLIFY (ks_block, s_block, p_block) 207 CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk) 208 ria = particle_set(irow)%r 209 ria = pbc(ria, cell) 210 rib = particle_set(icol)%r 211 rib = pbc(rib, cell) 212 fdir = 0.5_dp*SUM(fieldpol(1:3)*(ria(1:3) + rib(1:3))) 213 DO ispin = 1, nspin 214 CALL dbcsr_get_block_p(matrix=ks_matrix(ispin, 1)%matrix, & 215 row=irow, col=icol, BLOCK=ks_block, found=found) 216 ks_block = ks_block + fdir*s_block 217 CPASSERT(found) 218 END DO 219 IF (calculate_forces) THEN 220 ikind = kind_of(irow) 221 jkind = kind_of(icol) 222 atom_a = atom_of_kind(irow) 223 atom_b = atom_of_kind(icol) 224 fij = 0.0_dp 225 DO ispin = 1, nspin 226 CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, & 227 row=irow, col=icol, BLOCK=p_block, found=found) 228 CPASSERT(found) 229 DO idir = 1, 3 230 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1)%matrix, & 231 row=irow, col=icol, BLOCK=ds_block, found=found) 232 CPASSERT(found) 233 fij(idir) = fij(idir) + SUM(p_block*ds_block) 234 END DO 235 END DO 236 fdir = SUM(ria(1:3)*fieldpol(1:3)) 237 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3) 238 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3) 239 fdir = SUM(rib(1:3)*fieldpol(1:3)) 240 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3) 241 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3) 242 END IF 243 ENDDO 244 CALL dbcsr_iterator_stop(iter) 245 246 IF (calculate_forces) THEN 247 DO ikind = 1, SIZE(atomic_kind_set) 248 CALL mp_sum(force(ikind)%efield, para_env%group) 249 END DO 250 DEALLOCATE (atom_of_kind, kind_of) 251 END IF 252 253 END IF 254 255 CALL timestop(handle) 256 257 END SUBROUTINE efield_tb_local 258 259! ************************************************************************************************** 260!> \brief ... 261!> \param qs_env ... 262!> \param ks_matrix ... 263!> \param rho ... 264!> \param mcharge ... 265!> \param energy ... 266!> \param calculate_forces ... 267!> \param just_energy ... 268! ************************************************************************************************** 269 SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) 270 TYPE(qs_environment_type), POINTER :: qs_env 271 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix 272 TYPE(qs_rho_type), POINTER :: rho 273 REAL(dp), DIMENSION(:), INTENT(in) :: mcharge 274 TYPE(qs_energy_type), POINTER :: energy 275 LOGICAL, INTENT(in) :: calculate_forces, just_energy 276 277 CHARACTER(LEN=*), PARAMETER :: routineN = 'efield_tb_berry', & 278 routineP = moduleN//':'//routineN 279 280 COMPLEX(KIND=dp) :: zdeta 281 COMPLEX(KIND=dp), DIMENSION(3) :: zi(3) 282 INTEGER :: atom_a, atom_b, blk, handle, ia, iatom, & 283 ic, icol, idir, ikind, irow, is, & 284 ispin, jatom, jkind, natom, nimg, nspin 285 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of 286 INTEGER, DIMENSION(3) :: cellind 287 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 288 LOGICAL :: found, use_virial 289 REAL(KIND=dp) :: charge, dd, fdir 290 REAL(KIND=dp), DIMENSION(3) :: fieldpol, fij, forcea, fpolvec, kvec, & 291 qi, rab, ria, rib 292 REAL(KIND=dp), DIMENSION(3, 3) :: hmat 293 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block 294 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 295 TYPE(cell_type), POINTER :: cell 296 TYPE(cp_para_env_type), POINTER :: para_env 297 TYPE(dbcsr_iterator_type) :: iter 298 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s 299 TYPE(dft_control_type), POINTER :: dft_control 300 TYPE(kpoint_type), POINTER :: kpoints 301 TYPE(neighbor_list_iterator_p_type), & 302 DIMENSION(:), POINTER :: nl_iterator 303 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 304 POINTER :: sab_orb 305 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 306 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 307 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 308 TYPE(virial_type), POINTER :: virial 309 310 CALL timeset(routineN, handle) 311 312 NULLIFY (dft_control, cell, particle_set) 313 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, & 314 particle_set=particle_set, virial=virial) 315 NULLIFY (qs_kind_set, para_env, sab_orb) 316 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, & 317 energy=energy, para_env=para_env, sab_orb=sab_orb) 318 319 ! calculate stress only if forces requested also 320 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 321 use_virial = use_virial .AND. calculate_forces 322 ! disable stress calculation 323 IF (use_virial) THEN 324 CPABORT("Stress tensor for periodic E-field not implemented") 325 END IF 326 327 fieldpol = dft_control%period_efield%polarisation 328 fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) 329 fieldpol = -fieldpol*dft_control%period_efield%strength 330 hmat = cell%hmat(:, :)/twopi 331 DO idir = 1, 3 332 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir) 333 END DO 334 335 natom = SIZE(particle_set) 336 nspin = SIZE(ks_matrix, 1) 337 338 zi(:) = CMPLX(1._dp, 0._dp, dp) 339 DO ia = 1, natom 340 charge = mcharge(ia) 341 ria = particle_set(ia)%r 342 DO idir = 1, 3 343 kvec(:) = twopi*cell%h_inv(idir, :) 344 dd = SUM(kvec(:)*ria(:)) 345 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge 346 zi(idir) = zi(idir)*zdeta 347 END DO 348 END DO 349 qi = AIMAG(LOG(zi)) 350 energy%efield = -SUM(fpolvec(:)*qi(:)) 351 352 IF (.NOT. just_energy) THEN 353 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s) 354 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 355 356 nimg = dft_control%nimages 357 NULLIFY (cell_to_index) 358 IF (nimg > 1) THEN 359 NULLIFY (kpoints) 360 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints) 361 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 362 END IF 363 364 IF (calculate_forces) THEN 365 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force) 366 ALLOCATE (atom_of_kind(natom), kind_of(natom)) 367 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of) 368 IF (para_env%mepos == 0) THEN 369 DO ia = 1, natom 370 charge = -mcharge(ia) 371 iatom = atom_of_kind(ia) 372 ikind = kind_of(ia) 373 force(ikind)%efield(:, iatom) = fieldpol(:)*charge 374 IF (use_virial) THEN 375 ria = particle_set(ia)%r 376 ria = pbc(ria, cell) 377 CALL virial_pair_force(virial%pv_virial, 1.0_dp, forcea, ria) 378 END IF 379 END DO 380 ELSE 381 DO ia = 1, natom 382 iatom = atom_of_kind(ia) 383 ikind = kind_of(ia) 384 force(ikind)%efield(:, iatom) = 0.0_dp 385 END DO 386 END IF 387 END IF 388 389 IF (nimg == 1) THEN 390 ! no k-points; all matrices have been transformed to periodic bsf 391 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix) 392 DO WHILE (dbcsr_iterator_blocks_left(iter)) 393 CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk) 394 395 fdir = 0.0_dp 396 ria = particle_set(irow)%r 397 rib = particle_set(icol)%r 398 DO idir = 1, 3 399 kvec(:) = twopi*cell%h_inv(idir, :) 400 dd = SUM(kvec(:)*ria(:)) 401 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp) 402 fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta)) 403 dd = SUM(kvec(:)*rib(:)) 404 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp) 405 fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta)) 406 END DO 407 408 DO is = 1, nspin 409 NULLIFY (ks_block) 410 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, & 411 row=irow, col=icol, block=ks_block, found=found) 412 CPASSERT(found) 413 ks_block = ks_block + 0.5_dp*fdir*s_block 414 END DO 415 IF (calculate_forces) THEN 416 ikind = kind_of(irow) 417 jkind = kind_of(icol) 418 atom_a = atom_of_kind(irow) 419 atom_b = atom_of_kind(icol) 420 fij = 0.0_dp 421 DO ispin = 1, nspin 422 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, & 423 row=irow, col=icol, BLOCK=p_block, found=found) 424 CPASSERT(found) 425 DO idir = 1, 3 426 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, & 427 row=irow, col=icol, BLOCK=ds_block, found=found) 428 CPASSERT(found) 429 fij(idir) = fij(idir) + SUM(p_block*ds_block) 430 END DO 431 END DO 432 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3) 433 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3) 434 END IF 435 436 ENDDO 437 CALL dbcsr_iterator_stop(iter) 438 ELSE 439 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 440 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 441 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 442 iatom=iatom, jatom=jatom, r=rab, cell=cellind) 443 444 icol = MAX(iatom, jatom) 445 irow = MIN(iatom, jatom) 446 447 ic = cell_to_index(cellind(1), cellind(2), cellind(3)) 448 CPASSERT(ic > 0) 449 450 fdir = 0.0_dp 451 ria = particle_set(irow)%r 452 rib = particle_set(icol)%r 453 DO idir = 1, 3 454 kvec(:) = twopi*cell%h_inv(idir, :) 455 dd = SUM(kvec(:)*ria(:)) 456 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp) 457 fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta)) 458 dd = SUM(kvec(:)*rib(:)) 459 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp) 460 fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta)) 461 END DO 462 463 NULLIFY (s_block) 464 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, & 465 row=irow, col=icol, block=s_block, found=found) 466 CPASSERT(found) 467 DO is = 1, nspin 468 NULLIFY (ks_block) 469 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, & 470 row=irow, col=icol, block=ks_block, found=found) 471 CPASSERT(found) 472 ks_block = ks_block + 0.5_dp*fdir*s_block 473 END DO 474 IF (calculate_forces) THEN 475 atom_a = atom_of_kind(iatom) 476 atom_b = atom_of_kind(jatom) 477 fij = 0.0_dp 478 DO ispin = 1, nspin 479 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, & 480 row=irow, col=icol, BLOCK=p_block, found=found) 481 CPASSERT(found) 482 DO idir = 1, 3 483 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, & 484 row=irow, col=icol, BLOCK=ds_block, found=found) 485 CPASSERT(found) 486 fij(idir) = fij(idir) + SUM(p_block*ds_block) 487 END DO 488 END DO 489 IF (irow == iatom) fij = -fij 490 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3) 491 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) + fdir*fij(1:3) 492 END IF 493 494 END DO 495 CALL neighbor_list_iterator_release(nl_iterator) 496 END IF 497 498 IF (calculate_forces) THEN 499 DO ikind = 1, SIZE(atomic_kind_set) 500 CALL mp_sum(force(ikind)%efield, para_env%group) 501 END DO 502 DEALLOCATE (atom_of_kind, kind_of) 503 END IF 504 505 END IF 506 507 CALL timestop(handle) 508 509 END SUBROUTINE efield_tb_berry 510 511! ************************************************************************************************** 512!> \brief ... 513!> \param qs_env ... 514!> \param ks_matrix ... 515!> \param rho ... 516!> \param mcharge ... 517!> \param energy ... 518!> \param calculate_forces ... 519!> \param just_energy ... 520! ************************************************************************************************** 521 SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) 522 TYPE(qs_environment_type), POINTER :: qs_env 523 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix 524 TYPE(qs_rho_type), POINTER :: rho 525 REAL(dp), DIMENSION(:), INTENT(in) :: mcharge 526 TYPE(qs_energy_type), POINTER :: energy 527 LOGICAL, INTENT(in) :: calculate_forces, just_energy 528 529 CHARACTER(LEN=*), PARAMETER :: routineN = 'dfield_tb_berry', & 530 routineP = moduleN//':'//routineN 531 532 COMPLEX(KIND=dp) :: zdeta 533 COMPLEX(KIND=dp), DIMENSION(3) :: zi(3) 534 INTEGER :: atom_a, atom_b, blk, handle, i, ia, & 535 iatom, ic, icol, idir, ikind, irow, & 536 is, ispin, jatom, jkind, natom, nimg, & 537 nspin 538 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of 539 INTEGER, DIMENSION(3) :: cellind 540 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 541 LOGICAL :: found, use_virial 542 REAL(KIND=dp) :: charge, dd, ener_field, fdir, omega 543 REAL(KIND=dp), DIMENSION(3) :: ci, cqi, dfilter, di, fieldpol, fij, & 544 hdi, kvec, qi, rab, ria, rib 545 REAL(KIND=dp), DIMENSION(3, 3) :: hmat 546 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block 547 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 548 TYPE(cell_type), POINTER :: cell 549 TYPE(cp_para_env_type), POINTER :: para_env 550 TYPE(dbcsr_iterator_type) :: iter 551 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s 552 TYPE(dft_control_type), POINTER :: dft_control 553 TYPE(efield_berry_type), POINTER :: efield 554 TYPE(kpoint_type), POINTER :: kpoints 555 TYPE(neighbor_list_iterator_p_type), & 556 DIMENSION(:), POINTER :: nl_iterator 557 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 558 POINTER :: sab_orb 559 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 560 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 561 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 562 TYPE(virial_type), POINTER :: virial 563 564 CALL timeset(routineN, handle) 565 566 NULLIFY (dft_control, cell, particle_set) 567 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, & 568 particle_set=particle_set, virial=virial) 569 NULLIFY (qs_kind_set, para_env, sab_orb) 570 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, & 571 efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb) 572 573 ! efield history 574 CALL init_efield_matrices(efield) 575 CALL set_qs_env(qs_env, efield=efield) 576 577 ! calculate stress only if forces requested also 578 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 579 use_virial = use_virial .AND. calculate_forces 580 ! disable stress calculation 581 IF (use_virial) THEN 582 CPABORT("Stress tensor for periodic E-field not implemented") 583 END IF 584 585 dfilter(1:3) = dft_control%period_efield%d_filter(1:3) 586 587 fieldpol = dft_control%period_efield%polarisation 588 fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) 589 fieldpol = fieldpol*dft_control%period_efield%strength 590 591 omega = cell%deth 592 hmat = cell%hmat(:, :)/twopi 593 594 natom = SIZE(particle_set) 595 nspin = SIZE(ks_matrix, 1) 596 597 zi(:) = CMPLX(1._dp, 0._dp, dp) 598 DO ia = 1, natom 599 charge = mcharge(ia) 600 ria = particle_set(ia)%r 601 DO idir = 1, 3 602 kvec(:) = twopi*cell%h_inv(idir, :) 603 dd = SUM(kvec(:)*ria(:)) 604 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge 605 zi(idir) = zi(idir)*zdeta 606 END DO 607 END DO 608 qi = AIMAG(LOG(zi)) 609 610 ! make sure the total normalized polarization is within [-1:1] 611 DO idir = 1, 3 612 cqi(idir) = qi(idir)/omega 613 IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi 614 IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi 615 ! now check for log branch 616 IF (calculate_forces) THEN 617 IF (ABS(efield%polarisation(idir) - cqi(idir)) > pi) THEN 618 di(idir) = (efield%polarisation(idir) - cqi(idir))/pi 619 DO i = 1, 10 620 cqi(idir) = cqi(idir) + SIGN(1.0_dp, di(idir))*twopi 621 IF (ABS(efield%polarisation(idir) - cqi(idir)) < pi) EXIT 622 END DO 623 END IF 624 END IF 625 cqi(idir) = cqi(idir)*omega 626 END DO 627 DO idir = 1, 3 628 ci(idir) = 0.0_dp 629 DO i = 1, 3 630 ci(idir) = ci(idir) + hmat(idir, i)*cqi(i) 631 END DO 632 END DO 633 ! update the references 634 IF (calculate_forces) THEN 635 ener_field = SUM(ci) 636 ! check for smoothness of energy surface 637 IF (ABS(efield%field_energy - ener_field) > pi*ABS(SUM(hmat))) THEN 638 CPWARN("Large change of e-field energy detected. Correct for non-smooth energy surface") 639 END IF 640 efield%field_energy = ener_field 641 efield%polarisation(:) = cqi(:)/omega 642 END IF 643 644 ! Energy 645 ener_field = 0.0_dp 646 DO idir = 1, 3 647 ener_field = ener_field + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi/omega*ci(idir))**2 648 END DO 649 energy%efield = 0.25_dp/twopi*ener_field 650 651 IF (.NOT. just_energy) THEN 652 di(:) = -(fieldpol(:) - 2._dp*twopi/omega*ci(:))*dfilter(:)/omega 653 654 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s) 655 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 656 657 nimg = dft_control%nimages 658 NULLIFY (cell_to_index) 659 IF (nimg > 1) THEN 660 NULLIFY (kpoints) 661 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints) 662 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 663 END IF 664 665 IF (calculate_forces) THEN 666 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force) 667 ALLOCATE (atom_of_kind(natom), kind_of(natom)) 668 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of) 669 IF (para_env%mepos == 0) THEN 670 DO ia = 1, natom 671 charge = mcharge(ia) 672 iatom = atom_of_kind(ia) 673 ikind = kind_of(ia) 674 force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + di(:)*charge 675 IF (use_virial) THEN 676 ria = particle_set(ia)%r 677 ria = pbc(ria, cell) 678 CALL virial_pair_force(virial%pv_virial, 1.0_dp, di*charge, ria) 679 END IF 680 END DO 681 END IF 682 END IF 683 684 IF (nimg == 1) THEN 685 ! no k-points; all matrices have been transformed to periodic bsf 686 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix) 687 DO WHILE (dbcsr_iterator_blocks_left(iter)) 688 CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk) 689 690 DO idir = 1, 3 691 hdi(idir) = -SUM(di(1:3)*hmat(1:3, idir)) 692 END DO 693 fdir = 0.0_dp 694 ria = particle_set(irow)%r 695 rib = particle_set(icol)%r 696 DO idir = 1, 3 697 kvec(:) = twopi*cell%h_inv(idir, :) 698 dd = SUM(kvec(:)*ria(:)) 699 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp) 700 fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta)) 701 dd = SUM(kvec(:)*rib(:)) 702 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp) 703 fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta)) 704 END DO 705 706 DO is = 1, nspin 707 NULLIFY (ks_block) 708 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, & 709 row=irow, col=icol, block=ks_block, found=found) 710 CPASSERT(found) 711 ks_block = ks_block + 0.5_dp*fdir*s_block 712 END DO 713 IF (calculate_forces) THEN 714 ikind = kind_of(irow) 715 jkind = kind_of(icol) 716 atom_a = atom_of_kind(irow) 717 atom_b = atom_of_kind(icol) 718 fij = 0.0_dp 719 DO ispin = 1, nspin 720 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, & 721 row=irow, col=icol, BLOCK=p_block, found=found) 722 CPASSERT(found) 723 DO idir = 1, 3 724 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, & 725 row=irow, col=icol, BLOCK=ds_block, found=found) 726 CPASSERT(found) 727 fij(idir) = fij(idir) + SUM(p_block*ds_block) 728 END DO 729 END DO 730 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3) 731 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3) 732 END IF 733 734 ENDDO 735 CALL dbcsr_iterator_stop(iter) 736 ELSE 737 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 738 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 739 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 740 iatom=iatom, jatom=jatom, r=rab, cell=cellind) 741 742 icol = MAX(iatom, jatom) 743 irow = MIN(iatom, jatom) 744 745 ic = cell_to_index(cellind(1), cellind(2), cellind(3)) 746 CPASSERT(ic > 0) 747 748 DO idir = 1, 3 749 hdi(idir) = -SUM(di(1:3)*hmat(1:3, idir)) 750 END DO 751 fdir = 0.0_dp 752 ria = particle_set(irow)%r 753 rib = particle_set(icol)%r 754 DO idir = 1, 3 755 kvec(:) = twopi*cell%h_inv(idir, :) 756 dd = SUM(kvec(:)*ria(:)) 757 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp) 758 fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta)) 759 dd = SUM(kvec(:)*rib(:)) 760 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp) 761 fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta)) 762 END DO 763 764 NULLIFY (s_block) 765 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, & 766 row=irow, col=icol, block=s_block, found=found) 767 CPASSERT(found) 768 DO is = 1, nspin 769 NULLIFY (ks_block) 770 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, & 771 row=irow, col=icol, block=ks_block, found=found) 772 CPASSERT(found) 773 ks_block = ks_block + 0.5_dp*fdir*s_block 774 END DO 775 IF (calculate_forces) THEN 776 atom_a = atom_of_kind(iatom) 777 atom_b = atom_of_kind(jatom) 778 fij = 0.0_dp 779 DO ispin = 1, nspin 780 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, & 781 row=irow, col=icol, BLOCK=p_block, found=found) 782 CPASSERT(found) 783 DO idir = 1, 3 784 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, & 785 row=irow, col=icol, BLOCK=ds_block, found=found) 786 CPASSERT(found) 787 fij(idir) = fij(idir) + SUM(p_block*ds_block) 788 END DO 789 END DO 790 IF (irow == iatom) fij = -fij 791 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3) 792 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) + fdir*fij(1:3) 793 END IF 794 795 END DO 796 CALL neighbor_list_iterator_release(nl_iterator) 797 END IF 798 799 IF (calculate_forces) THEN 800 DEALLOCATE (atom_of_kind, kind_of) 801 END IF 802 803 END IF 804 805 CALL timestop(handle) 806 807 END SUBROUTINE dfield_tb_berry 808 809END MODULE efield_tb_methods 810