18 USE dbcsr_api,
ONLY: dbcsr_get_block_p,&
19 dbcsr_iterator_blocks_left,&
20 dbcsr_iterator_next_block,&
21 dbcsr_iterator_start,&
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)
78 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
80 REAL(
dp),
DIMENSION(:),
INTENT(in) :: mcharge
82 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
84 CHARACTER(len=*),
PARAMETER :: routinen =
'efield_tb_matrix'
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)
125 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
127 REAL(
dp),
DIMENSION(:),
INTENT(in) :: mcharge
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
143 TYPE(dbcsr_iterator_type) :: iter
144 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_p, matrix_s
149 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
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)
271 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
273 REAL(
dp),
DIMENSION(:),
INTENT(in) :: mcharge
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
297 TYPE(dbcsr_iterator_type) :: iter
298 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_s
303 DIMENSION(:),
POINTER :: nl_iterator
308 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
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)
595 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
597 REAL(
dp),
DIMENSION(:),
INTENT(in) :: mcharge
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
620 TYPE(dbcsr_iterator_type) :: iter
621 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_s
627 DIMENSION(:),
POINTER :: nl_iterator
632 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
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
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 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.
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.
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)
...
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.
keeps the density in various representations, keeping track of which ones are valid.