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
96 IF (dft_control%period_efield%start_frame <= qs_env%sim_step .AND. &
97 (dft_control%period_efield%end_frame == -1 .OR. dft_control%period_efield%end_frame >= qs_env%sim_step))
THEN
99 IF (dft_control%period_efield%displacement_field)
THEN
100 CALL dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
102 CALL efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
105 ELSE IF (dft_control%apply_efield)
THEN
106 CALL efield_tb_local(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
107 ELSE IF (dft_control%apply_efield_field)
THEN
108 cpabort(
"efield_filed")
111 cpabort(
"This routine should only be called from TB")
114 CALL timestop(handle)
128 SUBROUTINE efield_tb_local(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
130 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
132 REAL(
dp),
DIMENSION(:),
INTENT(in) :: mcharge
134 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
136 CHARACTER(LEN=*),
PARAMETER :: routinen =
'efield_tb_local'
138 INTEGER :: atom_a, atom_b, handle, ia, icol, idir, &
139 ikind, irow, ispin, jkind, natom, nspin
140 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
141 LOGICAL :: do_kpoints, found, use_virial
142 REAL(
dp) :: charge, fdir
143 REAL(
dp),
DIMENSION(3) :: ci, fieldpol, fij, ria, rib
144 REAL(
dp),
DIMENSION(:, :),
POINTER :: ds_block, ks_block, p_block, s_block
148 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_p, matrix_s
153 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
156 CALL timeset(routinen, handle)
158 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, particle_set=particle_set)
159 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, energy=energy, para_env=para_env)
160 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints, virial=virial)
162 cpabort(
"Local electric field with kpoints not possible. Use Berry phase periodic version")
165 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
167 cpabort(
"Stress tensor for non-periodic E-field not possible")
170 fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
171 dft_control%efield_fields(1)%efield%strength
173 natom =
SIZE(particle_set)
177 ria = particle_set(ia)%r
179 ci(:) = ci(:) + charge*ria(:)
181 energy%efield = -sum(ci(:)*fieldpol(:))
183 IF (.NOT. just_energy)
THEN
185 IF (calculate_forces)
THEN
186 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
188 IF (para_env%mepos == 0)
THEN
192 atom_a = atom_of_kind(ia)
193 force(ikind)%efield(1:3, atom_a) = -charge*fieldpol(:)
198 atom_a = atom_of_kind(ia)
199 force(ikind)%efield(1:3, atom_a) = 0.0_dp
206 nspin =
SIZE(ks_matrix, 1)
208 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
211 NULLIFY (ks_block, s_block, p_block)
213 ria = particle_set(irow)%r
215 rib = particle_set(icol)%r
217 fdir = 0.5_dp*sum(fieldpol(1:3)*(ria(1:3) + rib(1:3)))
220 row=irow, col=icol, block=ks_block, found=found)
221 ks_block = ks_block + fdir*s_block
224 IF (calculate_forces)
THEN
225 ikind = kind_of(irow)
226 jkind = kind_of(icol)
227 atom_a = atom_of_kind(irow)
228 atom_b = atom_of_kind(icol)
232 row=irow, col=icol, block=p_block, found=found)
236 row=irow, col=icol, block=ds_block, found=found)
238 fij(idir) = fij(idir) + sum(p_block*ds_block)
241 fdir = sum(ria(1:3)*fieldpol(1:3))
242 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
243 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
244 fdir = sum(rib(1:3)*fieldpol(1:3))
245 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
246 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
251 IF (calculate_forces)
THEN
252 DO ikind = 1,
SIZE(atomic_kind_set)
253 CALL para_env%sum(force(ikind)%efield)
259 CALL timestop(handle)
261 END SUBROUTINE efield_tb_local
273 SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
275 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
277 REAL(
dp),
DIMENSION(:),
INTENT(in) :: mcharge
279 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
281 CHARACTER(LEN=*),
PARAMETER :: routinen =
'efield_tb_berry'
283 COMPLEX(KIND=dp) :: zdeta
284 COMPLEX(KIND=dp),
DIMENSION(3) :: zi(3)
285 INTEGER :: atom_a, atom_b, handle, ia, iac, iatom, &
286 ic, icol, idir, ikind, irow, is, &
287 ispin, jatom, jkind, natom, nimg, &
289 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
290 INTEGER,
DIMENSION(3) :: cellind
291 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
292 LOGICAL :: found, use_virial
293 REAL(kind=
dp) :: charge, dd, dr, fdir, fi, strength
294 REAL(kind=
dp),
DIMENSION(3) :: fieldpol, fij, forcea, fpolvec, kvec, &
295 qi, rab, ria, rib, rij
296 REAL(kind=
dp),
DIMENSION(3, 3) :: hmat
297 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: ds_block, ks_block, p_block, s_block
298 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: dsint
302 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_s
307 DIMENSION(:),
POINTER :: nl_iterator
312 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
316 CALL timeset(routinen, handle)
318 NULLIFY (dft_control, cell, particle_set)
319 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
320 particle_set=particle_set, virial=virial)
321 NULLIFY (qs_kind_set, para_env, sab_orb)
322 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
323 energy=energy, para_env=para_env, sab_orb=sab_orb)
326 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
327 use_virial = use_virial .AND. calculate_forces
330 strength = dft_control%period_efield%strength
331 IF (
ALLOCATED(dft_control%period_efield%strength_list))
THEN
332 strength = dft_control%period_efield%strength_list(mod(qs_env%sim_step &
333 - dft_control%period_efield%start_frame,
SIZE(dft_control%period_efield%strength_list)) + 1)
336 fieldpol = dft_control%period_efield%polarisation
337 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
338 fieldpol = -fieldpol*strength
339 hmat = cell%hmat(:, :)/
twopi
341 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir)
344 natom =
SIZE(particle_set)
345 nspin =
SIZE(ks_matrix, 1)
347 zi(:) = cmplx(1._dp, 0._dp,
dp)
350 ria = particle_set(ia)%r
352 kvec(:) =
twopi*cell%h_inv(idir, :)
353 dd = sum(kvec(:)*ria(:))
354 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)**charge
355 zi(idir) = zi(idir)*zdeta
359 energy%efield = sum(fpolvec(:)*qi(:))
361 IF (.NOT. just_energy)
THEN
362 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
365 nimg = dft_control%nimages
366 NULLIFY (cell_to_index)
369 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
373 IF (calculate_forces)
THEN
374 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
376 IF (para_env%mepos == 0)
THEN
379 iatom = atom_of_kind(ia)
381 force(ikind)%efield(:, iatom) = fieldpol(:)*charge
383 ria = particle_set(ia)%r
385 forcea(1:3) = fieldpol(1:3)*charge
392 iatom = atom_of_kind(ia)
394 force(ikind)%efield(:, iatom) = 0.0_dp
406 ria = particle_set(irow)%r
407 rib = particle_set(icol)%r
409 kvec(:) =
twopi*cell%h_inv(idir, :)
410 dd = sum(kvec(:)*ria(:))
411 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
412 fdir = fdir + fpolvec(idir)*aimag(log(zdeta))
413 dd = sum(kvec(:)*rib(:))
414 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
415 fdir = fdir + fpolvec(idir)*aimag(log(zdeta))
421 row=irow, col=icol, block=ks_block, found=found)
423 ks_block = ks_block - 0.5_dp*fdir*s_block
425 IF (calculate_forces .AND. irow /= icol)
THEN
426 ikind = kind_of(irow)
427 jkind = kind_of(icol)
428 atom_a = atom_of_kind(irow)
429 atom_b = atom_of_kind(icol)
433 row=irow, col=icol, block=p_block, found=found)
437 row=irow, col=icol, block=ds_block, found=found)
439 fij(idir) = fij(idir) - sum(p_block*ds_block)
442 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
443 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
453 IF (dft_control%qs_control%dftb)
THEN
454 cpabort(
"DFTB stress tensor for periodic efield not implemented")
455 ELSEIF (dft_control%qs_control%xtb)
THEN
458 cpabort(
"TB method unknown")
464 iac = ikind + nkind*(jkind - 1)
465 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist)) cycle
466 DO ia = 1, sap_int(iac)%nalist
467 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist(ia)%clist)) cycle
468 iatom = sap_int(iac)%alist(ia)%aatom
469 DO ic = 1, sap_int(iac)%alist(ia)%nclist
470 jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
471 rij = sap_int(iac)%alist(ia)%clist(ic)%rac
472 dr = sqrt(sum(rij(:)**2))
473 IF (dr > 1.e-6_dp)
THEN
474 dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
475 icol = max(iatom, jatom)
476 irow = min(iatom, jatom)
477 IF (irow == iatom) rij = -rij
479 ria = particle_set(irow)%r
480 rib = particle_set(icol)%r
482 kvec(:) =
twopi*cell%h_inv(idir, :)
483 dd = sum(kvec(:)*ria(:))
484 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
485 fdir = fdir - fpolvec(idir)*aimag(log(zdeta))
486 dd = sum(kvec(:)*rib(:))
487 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
488 fdir = fdir - fpolvec(idir)*aimag(log(zdeta))
491 IF (iatom == jatom) fi = 0.5_dp
495 row=irow, col=icol, block=p_block, found=found)
499 IF (irow == iatom)
THEN
500 fij(idir) = sum(p_block*dsint(:, :, idir))
502 fij(idir) = sum(transpose(p_block)*dsint(:, :, idir))
505 IF (irow == iatom) fij = -fij
519 iatom=iatom, jatom=jatom, r=rab, cell=cellind)
521 icol = max(iatom, jatom)
522 irow = min(iatom, jatom)
524 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
528 ria = particle_set(irow)%r
529 rib = particle_set(icol)%r
531 kvec(:) =
twopi*cell%h_inv(idir, :)
532 dd = sum(kvec(:)*ria(:))
533 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
534 fdir = fdir + fpolvec(idir)*aimag(log(zdeta))
535 dd = sum(kvec(:)*rib(:))
536 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
537 fdir = fdir + fpolvec(idir)*aimag(log(zdeta))
542 row=irow, col=icol, block=s_block, found=found)
547 row=irow, col=icol, block=ks_block, found=found)
549 ks_block = ks_block - 0.5_dp*fdir*s_block
551 IF (calculate_forces)
THEN
552 atom_a = atom_of_kind(iatom)
553 atom_b = atom_of_kind(jatom)
557 row=irow, col=icol, block=p_block, found=found)
561 row=irow, col=icol, block=ds_block, found=found)
563 fij(idir) = fij(idir) - sum(p_block*ds_block)
566 IF (irow == iatom) fij = -fij
567 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3)
568 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) + fdir*fij(1:3)
570 dr = sqrt(sum(rab(:)**2))
571 IF (dr > 1.e-6_dp)
THEN
573 IF (iatom == jatom) fi = 0.5_dp
582 IF (calculate_forces)
THEN
583 DO ikind = 1,
SIZE(atomic_kind_set)
584 CALL para_env%sum(force(ikind)%efield)
590 CALL timestop(handle)
592 END SUBROUTINE efield_tb_berry
604 SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
606 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
608 REAL(
dp),
DIMENSION(:),
INTENT(in) :: mcharge
610 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
612 CHARACTER(LEN=*),
PARAMETER :: routinen =
'dfield_tb_berry'
614 COMPLEX(KIND=dp) :: zdeta
615 COMPLEX(KIND=dp),
DIMENSION(3) :: zi(3)
616 INTEGER :: atom_a, atom_b, handle, i, ia, iatom, &
617 ic, icol, idir, ikind, irow, is, &
618 ispin, jatom, jkind, natom, nimg, nspin
619 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
620 INTEGER,
DIMENSION(3) :: cellind
621 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
622 LOGICAL :: found, use_virial
623 REAL(kind=
dp) :: charge, dd, ener_field, fdir, omega, &
625 REAL(kind=
dp),
DIMENSION(3) :: ci, cqi, dfilter, di, fieldpol, fij, &
626 hdi, kvec, qi, rab, ria, rib
627 REAL(kind=
dp),
DIMENSION(3, 3) :: hmat
628 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: ds_block, ks_block, p_block, s_block
632 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_s
638 DIMENSION(:),
POINTER :: nl_iterator
643 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
646 CALL timeset(routinen, handle)
648 NULLIFY (dft_control, cell, particle_set)
649 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
650 particle_set=particle_set, virial=virial)
651 NULLIFY (qs_kind_set, para_env, sab_orb)
652 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
653 efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
660 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
661 use_virial = use_virial .AND. calculate_forces
664 cpabort(
"Stress tensor for periodic D-field not implemented")
667 dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
670 strength = dft_control%period_efield%strength
671 IF (
ALLOCATED(dft_control%period_efield%strength_list))
THEN
672 strength = dft_control%period_efield%strength_list(mod(qs_env%sim_step &
673 - dft_control%period_efield%start_frame,
SIZE(dft_control%period_efield%strength_list)) + 1)
676 fieldpol = dft_control%period_efield%polarisation
677 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
678 fieldpol = fieldpol*strength
681 hmat = cell%hmat(:, :)/
twopi
683 natom =
SIZE(particle_set)
684 nspin =
SIZE(ks_matrix, 1)
686 zi(:) = cmplx(1._dp, 0._dp,
dp)
689 ria = particle_set(ia)%r
691 kvec(:) =
twopi*cell%h_inv(idir, :)
692 dd = sum(kvec(:)*ria(:))
693 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)**charge
694 zi(idir) = zi(idir)*zdeta
701 cqi(idir) = qi(idir)/omega
702 IF (cqi(idir) >
pi) cqi(idir) = cqi(idir) -
twopi
703 IF (cqi(idir) < -
pi) cqi(idir) = cqi(idir) +
twopi
705 IF (calculate_forces)
THEN
706 IF (abs(efield%polarisation(idir) - cqi(idir)) >
pi)
THEN
707 di(idir) = (efield%polarisation(idir) - cqi(idir))/
pi
709 cqi(idir) = cqi(idir) + sign(1.0_dp, di(idir))*
twopi
710 IF (abs(efield%polarisation(idir) - cqi(idir)) <
pi)
EXIT
714 cqi(idir) = cqi(idir)*omega
719 ci(idir) = ci(idir) + hmat(idir, i)*cqi(i)
723 IF (calculate_forces)
THEN
726 IF (abs(efield%field_energy - ener_field) >
pi*abs(sum(hmat)))
THEN
727 cpwarn(
"Large change of e-field energy detected. Correct for non-smooth energy surface")
729 efield%field_energy = ener_field
730 efield%polarisation(:) = cqi(:)/omega
736 ener_field = ener_field + dfilter(idir)*(fieldpol(idir) - 2._dp*
twopi/omega*ci(idir))**2
738 energy%efield = 0.25_dp/
twopi*ener_field
740 IF (.NOT. just_energy)
THEN
741 di(:) = -(fieldpol(:) - 2._dp*
twopi/omega*ci(:))*dfilter(:)/omega
743 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
746 nimg = dft_control%nimages
747 NULLIFY (cell_to_index)
750 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
754 IF (calculate_forces)
THEN
755 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
757 IF (para_env%mepos == 0)
THEN
760 iatom = atom_of_kind(ia)
762 force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + di(:)*charge
774 hdi(idir) = -sum(di(1:3)*hmat(1:3, idir))
777 ria = particle_set(irow)%r
778 rib = particle_set(icol)%r
780 kvec(:) =
twopi*cell%h_inv(idir, :)
781 dd = sum(kvec(:)*ria(:))
782 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
783 fdir = fdir + hdi(idir)*aimag(log(zdeta))
784 dd = sum(kvec(:)*rib(:))
785 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
786 fdir = fdir + hdi(idir)*aimag(log(zdeta))
792 row=irow, col=icol, block=ks_block, found=found)
794 ks_block = ks_block + 0.5_dp*fdir*s_block
796 IF (calculate_forces)
THEN
797 ikind = kind_of(irow)
798 jkind = kind_of(icol)
799 atom_a = atom_of_kind(irow)
800 atom_b = atom_of_kind(icol)
804 row=irow, col=icol, block=p_block, found=found)
808 row=irow, col=icol, block=ds_block, found=found)
810 fij(idir) = fij(idir) + sum(p_block*ds_block)
813 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
814 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
823 iatom=iatom, jatom=jatom, r=rab, cell=cellind)
825 icol = max(iatom, jatom)
826 irow = min(iatom, jatom)
828 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
832 hdi(idir) = -sum(di(1:3)*hmat(1:3, idir))
835 ria = particle_set(irow)%r
836 rib = particle_set(icol)%r
838 kvec(:) =
twopi*cell%h_inv(idir, :)
839 dd = sum(kvec(:)*ria(:))
840 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
841 fdir = fdir + hdi(idir)*aimag(log(zdeta))
842 dd = sum(kvec(:)*rib(:))
843 zdeta = cmplx(cos(dd), sin(dd), kind=
dp)
844 fdir = fdir + hdi(idir)*aimag(log(zdeta))
849 row=irow, col=icol, block=s_block, found=found)
854 row=irow, col=icol, block=ks_block, found=found)
856 ks_block = ks_block + 0.5_dp*fdir*s_block
858 IF (calculate_forces)
THEN
859 atom_a = atom_of_kind(iatom)
860 atom_b = atom_of_kind(jatom)
864 row=irow, col=icol, block=p_block, found=found)
868 row=irow, col=icol, block=ds_block, found=found)
870 fij(idir) = fij(idir) + sum(p_block*ds_block)
873 IF (irow == iatom) fij = -fij
874 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3)
875 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) + fdir*fij(1:3)
884 CALL timestop(handle)
886 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...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
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_pp, 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, harris_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, eeq, 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, harris_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, eeq, 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)
...
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.