18 USE dbcsr_api,
ONLY: dbcsr_get_block_p,&
19 dbcsr_iterator_blocks_left,&
20 dbcsr_iterator_next_block,&
21 dbcsr_iterator_start,&
41 neighbor_list_iterator_p_type,&
43 neighbor_list_set_p_type
53 #include "./base/base_uses.f90"
59 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'efield_tb_methods'
75 SUBROUTINE efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
77 TYPE(qs_environment_type),
POINTER :: qs_env
78 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
79 TYPE(qs_rho_type),
POINTER :: rho
80 REAL(
dp),
DIMENSION(:),
INTENT(in) :: mcharge
81 TYPE(qs_energy_type),
POINTER :: energy
82 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
84 CHARACTER(len=*),
PARAMETER :: routinen =
'efield_tb_matrix'
87 TYPE(dft_control_type),
POINTER :: dft_control
89 CALL timeset(routinen, handle)
91 energy%efield = 0.0_dp
92 CALL get_qs_env(qs_env, dft_control=dft_control)
93 IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb)
THEN
94 IF (dft_control%apply_period_efield)
THEN
95 IF (dft_control%period_efield%displacement_field)
THEN
96 CALL dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
98 CALL efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
100 ELSE IF (dft_control%apply_efield)
THEN
101 CALL efield_tb_local(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
102 ELSE IF (dft_control%apply_efield_field)
THEN
103 cpabort(
"efield_filed")
106 cpabort(
"This routine should only be called from TB")
109 CALL timestop(handle)
123 SUBROUTINE efield_tb_local(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
124 TYPE(qs_environment_type),
POINTER :: qs_env
125 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
126 TYPE(qs_rho_type),
POINTER :: rho
127 REAL(
dp),
DIMENSION(:),
INTENT(in) :: mcharge
128 TYPE(qs_energy_type),
POINTER :: energy
129 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
131 CHARACTER(LEN=*),
PARAMETER :: routinen =
'efield_tb_local'
133 INTEGER :: atom_a, atom_b, blk, handle, ia, icol, &
134 idir, ikind, irow, ispin, jkind, &
136 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
137 LOGICAL :: do_kpoints, found, use_virial
138 REAL(
dp) :: charge, fdir
139 REAL(
dp),
DIMENSION(3) :: ci, fieldpol, fij, ria, rib
140 REAL(
dp),
DIMENSION(:, :),
POINTER :: ds_block, ks_block, p_block, s_block
141 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
142 TYPE(cell_type),
POINTER :: cell
143 TYPE(dbcsr_iterator_type) :: iter
144 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_p, matrix_s
145 TYPE(dft_control_type),
POINTER :: dft_control
146 TYPE(mp_para_env_type),
POINTER :: para_env
147 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
148 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
149 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
150 TYPE(virial_type),
POINTER :: virial
152 CALL timeset(routinen, handle)
154 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, particle_set=particle_set)
155 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, energy=energy, para_env=para_env)
156 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints, virial=virial)
158 cpabort(
"Local electric field with kpoints not possible. Use Berry phase periodic version")
161 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
163 cpabort(
"Stress tensor for non-periodic E-field not possible")
166 fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
167 dft_control%efield_fields(1)%efield%strength
169 natom =
SIZE(particle_set)
173 ria = particle_set(ia)%r
175 ci(:) = ci(:) + charge*ria(:)
177 energy%efield = -sum(ci(:)*fieldpol(:))
179 IF (.NOT. just_energy)
THEN
181 IF (calculate_forces)
THEN
182 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
184 IF (para_env%mepos == 0)
THEN
188 atom_a = atom_of_kind(ia)
189 force(ikind)%efield(1:3, atom_a) = -charge*fieldpol(:)
194 atom_a = atom_of_kind(ia)
195 force(ikind)%efield(1:3, atom_a) = 0.0_dp
202 nspin =
SIZE(ks_matrix, 1)
204 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
205 CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
206 DO WHILE (dbcsr_iterator_blocks_left(iter))
207 NULLIFY (ks_block, s_block, p_block)
208 CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)
209 ria = particle_set(irow)%r
211 rib = particle_set(icol)%r
213 fdir = 0.5_dp*sum(fieldpol(1:3)*(ria(1:3) + rib(1:3)))
215 CALL dbcsr_get_block_p(matrix=ks_matrix(ispin, 1)%matrix, &
216 row=irow, col=icol, block=ks_block, found=found)
217 ks_block = ks_block + fdir*s_block
220 IF (calculate_forces)
THEN
221 ikind = kind_of(irow)
222 jkind = kind_of(icol)
223 atom_a = atom_of_kind(irow)
224 atom_b = atom_of_kind(icol)
227 CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
228 row=irow, col=icol, block=p_block, found=found)
231 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1)%matrix, &
232 row=irow, col=icol, block=ds_block, found=found)
234 fij(idir) = fij(idir) + sum(p_block*ds_block)
237 fdir = sum(ria(1:3)*fieldpol(1:3))
238 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
239 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
240 fdir = sum(rib(1:3)*fieldpol(1:3))
241 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
242 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
245 CALL dbcsr_iterator_stop(iter)
247 IF (calculate_forces)
THEN
248 DO ikind = 1,
SIZE(atomic_kind_set)
249 CALL para_env%sum(force(ikind)%efield)
255 CALL timestop(handle)
257 END SUBROUTINE efield_tb_local
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
277 CHARACTER(LEN=*),
PARAMETER :: routinen =
'efield_tb_berry'
279 COMPLEX(KIND=dp) :: zdeta
280 COMPLEX(KIND=dp),
DIMENSION(3) :: zi(3)
281 INTEGER :: atom_a, atom_b, blk, handle, ia, iac, &
282 iatom, ic, icol, idir, ikind, irow, &
283 is, ispin, jatom, jkind, natom, nimg, &
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, dr, fdir, fi
290 REAL(kind=
dp),
DIMENSION(3) :: fieldpol, fij, forcea, fpolvec, kvec, &
291 qi, rab, ria, rib, rij
292 REAL(kind=
dp),
DIMENSION(3, 3) :: hmat
293 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: ds_block, ks_block, p_block, s_block
294 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: dsint
295 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
296 TYPE(cell_type),
POINTER :: cell
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(mp_para_env_type),
POINTER :: para_env
302 TYPE(neighbor_list_iterator_p_type), &
303 DIMENSION(:),
POINTER :: nl_iterator
304 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
306 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
307 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
308 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
309 TYPE(sap_int_type),
DIMENSION(:),
POINTER :: sap_int
310 TYPE(virial_type),
POINTER :: virial
312 CALL timeset(routinen, handle)
314 NULLIFY (dft_control, cell, particle_set)
315 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
316 particle_set=particle_set, virial=virial)
317 NULLIFY (qs_kind_set, para_env, sab_orb)
318 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
319 energy=energy, para_env=para_env, sab_orb=sab_orb)
322 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
323 use_virial = use_virial .AND. calculate_forces
325 fieldpol = dft_control%period_efield%polarisation
326 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
327 fieldpol = -fieldpol*dft_control%period_efield%strength
328 hmat = cell%hmat(:, :)/
twopi
330 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir)
333 natom =
SIZE(particle_set)
334 nspin =
SIZE(ks_matrix, 1)
336 zi(:) = cmplx(1._dp, 0._dp,
dp)
339 ria = particle_set(ia)%r
341 kvec(:) =
twopi*cell%h_inv(idir, :)
342 dd = sum(kvec(:)*ria(:))
343 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)**charge
344 zi(idir) = zi(idir)*zdeta
348 energy%efield = -sum(fpolvec(:)*qi(:))
350 IF (.NOT. just_energy)
THEN
351 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
354 nimg = dft_control%nimages
355 NULLIFY (cell_to_index)
358 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
362 IF (calculate_forces)
THEN
363 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
365 IF (para_env%mepos == 0)
THEN
367 charge = -mcharge(ia)
368 iatom = atom_of_kind(ia)
370 force(ikind)%efield(:, iatom) = fieldpol(:)*charge
372 ria = particle_set(ia)%r
374 forcea(1:3) = fieldpol(1:3)*charge
381 iatom = atom_of_kind(ia)
383 force(ikind)%efield(:, iatom) = 0.0_dp
390 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
391 DO WHILE (dbcsr_iterator_blocks_left(iter))
392 CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)
395 ria = particle_set(irow)%r
396 rib = particle_set(icol)%r
398 kvec(:) =
twopi*cell%h_inv(idir, :)
399 dd = sum(kvec(:)*ria(:))
400 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
401 fdir = fdir + fpolvec(idir)*aimag(log(zdeta))
402 dd = sum(kvec(:)*rib(:))
403 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
404 fdir = fdir + fpolvec(idir)*aimag(log(zdeta))
409 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
410 row=irow, col=icol, block=ks_block, found=found)
412 ks_block = ks_block + 0.5_dp*fdir*s_block
414 IF (calculate_forces)
THEN
415 ikind = kind_of(irow)
416 jkind = kind_of(icol)
417 atom_a = atom_of_kind(irow)
418 atom_b = atom_of_kind(icol)
421 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
422 row=irow, col=icol, block=p_block, found=found)
425 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
426 row=irow, col=icol, block=ds_block, found=found)
428 fij(idir) = fij(idir) + sum(p_block*ds_block)
431 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
432 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
435 CALL dbcsr_iterator_stop(iter)
442 IF (dft_control%qs_control%dftb)
THEN
443 cpabort(
"DFTB stress tensor for periodic efield not implemented")
444 ELSEIF (dft_control%qs_control%xtb)
THEN
447 cpabort(
"TB method unknown")
453 iac = ikind + nkind*(jkind - 1)
454 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist)) cycle
455 DO ia = 1, sap_int(iac)%nalist
456 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist(ia)%clist)) cycle
457 iatom = sap_int(iac)%alist(ia)%aatom
458 DO ic = 1, sap_int(iac)%alist(ia)%nclist
459 jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
460 rij = sap_int(iac)%alist(ia)%clist(ic)%rac
461 dr = sqrt(sum(rij(:)**2))
462 IF (dr > 1.e-6_dp)
THEN
463 dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
464 icol = max(iatom, jatom)
465 irow = min(iatom, jatom)
466 IF (irow == iatom) rij = -rij
468 ria = particle_set(irow)%r
469 rib = particle_set(icol)%r
471 kvec(:) =
twopi*cell%h_inv(idir, :)
472 dd = sum(kvec(:)*ria(:))
473 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
474 fdir = fdir + fpolvec(idir)*aimag(log(zdeta))
475 dd = sum(kvec(:)*rib(:))
476 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
477 fdir = fdir + fpolvec(idir)*aimag(log(zdeta))
480 IF (iatom == jatom) fi = 0.5_dp
483 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
484 row=irow, col=icol, block=p_block, found=found)
488 IF (irow == iatom)
THEN
489 fij(idir) = sum(p_block*dsint(:, :, idir))
491 fij(idir) = sum(transpose(p_block)*dsint(:, :, idir))
494 IF (irow == iatom) fij = -fij
508 iatom=iatom, jatom=jatom, r=rab, cell=cellind)
510 icol = max(iatom, jatom)
511 irow = min(iatom, jatom)
513 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
517 ria = particle_set(irow)%r
518 rib = particle_set(icol)%r
520 kvec(:) =
twopi*cell%h_inv(idir, :)
521 dd = sum(kvec(:)*ria(:))
522 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
523 fdir = fdir + fpolvec(idir)*aimag(log(zdeta))
524 dd = sum(kvec(:)*rib(:))
525 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
526 fdir = fdir + fpolvec(idir)*aimag(log(zdeta))
530 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
531 row=irow, col=icol, block=s_block, found=found)
535 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
536 row=irow, col=icol, block=ks_block, found=found)
538 ks_block = ks_block + 0.5_dp*fdir*s_block
540 IF (calculate_forces)
THEN
541 atom_a = atom_of_kind(iatom)
542 atom_b = atom_of_kind(jatom)
545 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
546 row=irow, col=icol, block=p_block, found=found)
549 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
550 row=irow, col=icol, block=ds_block, found=found)
552 fij(idir) = fij(idir) + sum(p_block*ds_block)
555 IF (irow == iatom) fij = -fij
556 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3)
557 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) + fdir*fij(1:3)
559 dr = sqrt(sum(rab(:)**2))
560 IF (dr > 1.e-6_dp)
THEN
562 IF (iatom == jatom) fi = 0.5_dp
571 IF (calculate_forces)
THEN
572 DO ikind = 1,
SIZE(atomic_kind_set)
573 CALL para_env%sum(force(ikind)%efield)
579 CALL timestop(handle)
581 END SUBROUTINE efield_tb_berry
593 SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
594 TYPE(qs_environment_type),
POINTER :: qs_env
595 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
596 TYPE(qs_rho_type),
POINTER :: rho
597 REAL(
dp),
DIMENSION(:),
INTENT(in) :: mcharge
598 TYPE(qs_energy_type),
POINTER :: energy
599 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
601 CHARACTER(LEN=*),
PARAMETER :: routinen =
'dfield_tb_berry'
603 COMPLEX(KIND=dp) :: zdeta
604 COMPLEX(KIND=dp),
DIMENSION(3) :: zi(3)
605 INTEGER :: atom_a, atom_b, blk, handle, i, ia, &
606 iatom, ic, icol, idir, ikind, irow, &
607 is, ispin, jatom, jkind, natom, nimg, &
609 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
610 INTEGER,
DIMENSION(3) :: cellind
611 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
612 LOGICAL :: found, use_virial
613 REAL(kind=
dp) :: charge, dd, ener_field, fdir, omega
614 REAL(kind=
dp),
DIMENSION(3) :: ci, cqi, dfilter, di, fieldpol, fij, &
615 hdi, kvec, qi, rab, ria, rib
616 REAL(kind=
dp),
DIMENSION(3, 3) :: hmat
617 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: ds_block, ks_block, p_block, s_block
618 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
619 TYPE(cell_type),
POINTER :: cell
620 TYPE(dbcsr_iterator_type) :: iter
621 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_s
622 TYPE(dft_control_type),
POINTER :: dft_control
623 TYPE(efield_berry_type),
POINTER :: efield
624 TYPE(kpoint_type),
POINTER :: kpoints
625 TYPE(mp_para_env_type),
POINTER :: para_env
626 TYPE(neighbor_list_iterator_p_type), &
627 DIMENSION(:),
POINTER :: nl_iterator
628 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
630 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
631 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
632 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
633 TYPE(virial_type),
POINTER :: virial
635 CALL timeset(routinen, handle)
637 NULLIFY (dft_control, cell, particle_set)
638 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
639 particle_set=particle_set, virial=virial)
640 NULLIFY (qs_kind_set, para_env, sab_orb)
641 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
642 efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
649 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
650 use_virial = use_virial .AND. calculate_forces
653 cpabort(
"Stress tensor for periodic D-field not implemented")
656 dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
658 fieldpol = dft_control%period_efield%polarisation
659 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
660 fieldpol = fieldpol*dft_control%period_efield%strength
663 hmat = cell%hmat(:, :)/
twopi
665 natom =
SIZE(particle_set)
666 nspin =
SIZE(ks_matrix, 1)
668 zi(:) = cmplx(1._dp, 0._dp,
dp)
671 ria = particle_set(ia)%r
673 kvec(:) =
twopi*cell%h_inv(idir, :)
674 dd = sum(kvec(:)*ria(:))
675 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)**charge
676 zi(idir) = zi(idir)*zdeta
683 cqi(idir) = qi(idir)/omega
684 IF (cqi(idir) >
pi) cqi(idir) = cqi(idir) -
twopi
685 IF (cqi(idir) < -
pi) cqi(idir) = cqi(idir) +
twopi
687 IF (calculate_forces)
THEN
688 IF (abs(efield%polarisation(idir) - cqi(idir)) >
pi)
THEN
689 di(idir) = (efield%polarisation(idir) - cqi(idir))/
pi
691 cqi(idir) = cqi(idir) + sign(1.0_dp, di(idir))*
twopi
692 IF (abs(efield%polarisation(idir) - cqi(idir)) <
pi)
EXIT
696 cqi(idir) = cqi(idir)*omega
701 ci(idir) = ci(idir) + hmat(idir, i)*cqi(i)
705 IF (calculate_forces)
THEN
708 IF (abs(efield%field_energy - ener_field) >
pi*abs(sum(hmat)))
THEN
709 cpwarn(
"Large change of e-field energy detected. Correct for non-smooth energy surface")
711 efield%field_energy = ener_field
712 efield%polarisation(:) = cqi(:)/omega
718 ener_field = ener_field + dfilter(idir)*(fieldpol(idir) - 2._dp*
twopi/omega*ci(idir))**2
720 energy%efield = 0.25_dp/
twopi*ener_field
722 IF (.NOT. just_energy)
THEN
723 di(:) = -(fieldpol(:) - 2._dp*
twopi/omega*ci(:))*dfilter(:)/omega
725 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
728 nimg = dft_control%nimages
729 NULLIFY (cell_to_index)
732 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
736 IF (calculate_forces)
THEN
737 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
739 IF (para_env%mepos == 0)
THEN
742 iatom = atom_of_kind(ia)
744 force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + di(:)*charge
751 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
752 DO WHILE (dbcsr_iterator_blocks_left(iter))
753 CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)
756 hdi(idir) = -sum(di(1:3)*hmat(1:3, idir))
759 ria = particle_set(irow)%r
760 rib = particle_set(icol)%r
762 kvec(:) =
twopi*cell%h_inv(idir, :)
763 dd = sum(kvec(:)*ria(:))
764 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
765 fdir = fdir + hdi(idir)*aimag(log(zdeta))
766 dd = sum(kvec(:)*rib(:))
767 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
768 fdir = fdir + hdi(idir)*aimag(log(zdeta))
773 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
774 row=irow, col=icol, block=ks_block, found=found)
776 ks_block = ks_block + 0.5_dp*fdir*s_block
778 IF (calculate_forces)
THEN
779 ikind = kind_of(irow)
780 jkind = kind_of(icol)
781 atom_a = atom_of_kind(irow)
782 atom_b = atom_of_kind(icol)
785 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
786 row=irow, col=icol, block=p_block, found=found)
789 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
790 row=irow, col=icol, block=ds_block, found=found)
792 fij(idir) = fij(idir) + sum(p_block*ds_block)
795 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
796 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
800 CALL dbcsr_iterator_stop(iter)
805 iatom=iatom, jatom=jatom, r=rab, cell=cellind)
807 icol = max(iatom, jatom)
808 irow = min(iatom, jatom)
810 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
814 hdi(idir) = -sum(di(1:3)*hmat(1:3, idir))
817 ria = particle_set(irow)%r
818 rib = particle_set(icol)%r
820 kvec(:) =
twopi*cell%h_inv(idir, :)
821 dd = sum(kvec(:)*ria(:))
822 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
823 fdir = fdir + hdi(idir)*aimag(log(zdeta))
824 dd = sum(kvec(:)*rib(:))
825 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
826 fdir = fdir + hdi(idir)*aimag(log(zdeta))
830 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
831 row=irow, col=icol, block=s_block, found=found)
835 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
836 row=irow, col=icol, block=ks_block, found=found)
838 ks_block = ks_block + 0.5_dp*fdir*s_block
840 IF (calculate_forces)
THEN
841 atom_a = atom_of_kind(iatom)
842 atom_b = atom_of_kind(jatom)
845 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
846 row=irow, col=icol, block=p_block, found=found)
849 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
850 row=irow, col=icol, block=ds_block, found=found)
852 fij(idir) = fij(idir) + sum(p_block*ds_block)
855 IF (irow == iatom) fij = -fij
856 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3)
857 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) + fdir*fij(1:3)
866 CALL timestop(handle)
868 END SUBROUTINE dfield_tb_berry
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Handles all functions related to the CELL.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Calculation of electric field contributions in TB.
subroutine, public efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
...
Defines the basic variable types.
integer, parameter, public dp
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public twopi
Interface to the message passing library MPI.
Define the data structure for the particle information.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
type for berry phase efield matrices. At the moment only used for cosmat and sinmat
subroutine, public init_efield_matrices(efield)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
General overlap type integrals containers.
subroutine, public release_sap_int(sap_int)
...
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Calculation of Coulomb contributions in xTB.
subroutine, public xtb_dsint_list(qs_env, sap_int)
...