(git:b279b6b)
efield_tb_methods.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Calculation of electric field contributions in TB
10 !> \author JGH
11 ! **************************************************************************************************
13  USE atomic_kind_types, ONLY: atomic_kind_type,&
15  USE cell_types, ONLY: cell_type,&
16  pbc
17  USE cp_control_types, ONLY: dft_control_type
18  USE dbcsr_api, ONLY: dbcsr_get_block_p,&
19  dbcsr_iterator_blocks_left,&
20  dbcsr_iterator_next_block,&
21  dbcsr_iterator_start,&
22  dbcsr_iterator_stop,&
23  dbcsr_iterator_type,&
24  dbcsr_p_type
25  USE kinds, ONLY: dp
26  USE kpoint_types, ONLY: get_kpoint_info,&
27  kpoint_type
28  USE mathconstants, ONLY: pi,&
29  twopi
30  USE message_passing, ONLY: mp_para_env_type
31  USE particle_types, ONLY: particle_type
32  USE qs_energy_types, ONLY: qs_energy_type
33  USE qs_environment_types, ONLY: get_qs_env,&
34  qs_environment_type,&
36  USE qs_force_types, ONLY: qs_force_type
37  USE qs_kind_types, ONLY: qs_kind_type
41  neighbor_list_iterator_p_type,&
43  neighbor_list_set_p_type
44  USE qs_period_efield_types, ONLY: efield_berry_type,&
46  USE qs_rho_types, ONLY: qs_rho_get,&
47  qs_rho_type
48  USE sap_kind_types, ONLY: release_sap_int,&
49  sap_int_type
51  USE virial_types, ONLY: virial_type
52  USE xtb_coulomb, ONLY: xtb_dsint_list
53 #include "./base/base_uses.f90"
54 
55  IMPLICIT NONE
56 
57  PRIVATE
58 
59  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'efield_tb_methods'
60 
61  PUBLIC :: efield_tb_matrix
62 
63 CONTAINS
64 
65 ! **************************************************************************************************
66 !> \brief ...
67 !> \param qs_env ...
68 !> \param ks_matrix ...
69 !> \param rho ...
70 !> \param mcharge ...
71 !> \param energy ...
72 !> \param calculate_forces ...
73 !> \param just_energy ...
74 ! **************************************************************************************************
75  SUBROUTINE efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
76 
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
83 
84  CHARACTER(len=*), PARAMETER :: routinen = 'efield_tb_matrix'
85 
86  INTEGER :: handle
87  TYPE(dft_control_type), POINTER :: dft_control
88 
89  CALL timeset(routinen, handle)
90 
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)
97  ELSE
98  CALL efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
99  END IF
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")
104  END IF
105  ELSE
106  cpabort("This routine should only be called from TB")
107  END IF
108 
109  CALL timestop(handle)
110 
111  END SUBROUTINE efield_tb_matrix
112 
113 ! **************************************************************************************************
114 !> \brief ...
115 !> \param qs_env ...
116 !> \param ks_matrix ...
117 !> \param rho ...
118 !> \param mcharge ...
119 !> \param energy ...
120 !> \param calculate_forces ...
121 !> \param just_energy ...
122 ! **************************************************************************************************
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
130 
131  CHARACTER(LEN=*), PARAMETER :: routinen = 'efield_tb_local'
132 
133  INTEGER :: atom_a, atom_b, blk, handle, ia, icol, &
134  idir, ikind, irow, ispin, jkind, &
135  natom, nspin
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
151 
152  CALL timeset(routinen, handle)
153 
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)
157  IF (do_kpoints) THEN
158  cpabort("Local electric field with kpoints not possible. Use Berry phase periodic version")
159  END IF
160  ! disable stress calculation
161  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
162  IF (use_virial) THEN
163  cpabort("Stress tensor for non-periodic E-field not possible")
164  END IF
165 
166  fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
167  dft_control%efield_fields(1)%efield%strength
168 
169  natom = SIZE(particle_set)
170  ci = 0.0_dp
171  DO ia = 1, natom
172  charge = mcharge(ia)
173  ria = particle_set(ia)%r
174  ria = pbc(ria, cell)
175  ci(:) = ci(:) + charge*ria(:)
176  END DO
177  energy%efield = -sum(ci(:)*fieldpol(:))
178 
179  IF (.NOT. just_energy) THEN
180 
181  IF (calculate_forces) THEN
182  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
183  CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
184  IF (para_env%mepos == 0) THEN
185  DO ia = 1, natom
186  charge = mcharge(ia)
187  ikind = kind_of(ia)
188  atom_a = atom_of_kind(ia)
189  force(ikind)%efield(1:3, atom_a) = -charge*fieldpol(:)
190  END DO
191  ELSE
192  DO ia = 1, natom
193  ikind = kind_of(ia)
194  atom_a = atom_of_kind(ia)
195  force(ikind)%efield(1:3, atom_a) = 0.0_dp
196  END DO
197  END IF
198  CALL qs_rho_get(rho, rho_ao=matrix_p)
199  END IF
200 
201  ! Update KS matrix
202  nspin = SIZE(ks_matrix, 1)
203  NULLIFY (matrix_s)
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
210  ria = pbc(ria, cell)
211  rib = particle_set(icol)%r
212  rib = pbc(rib, cell)
213  fdir = 0.5_dp*sum(fieldpol(1:3)*(ria(1:3) + rib(1:3)))
214  DO ispin = 1, nspin
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
218  cpassert(found)
219  END DO
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)
225  fij = 0.0_dp
226  DO ispin = 1, nspin
227  CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
228  row=irow, col=icol, block=p_block, found=found)
229  cpassert(found)
230  DO idir = 1, 3
231  CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1)%matrix, &
232  row=irow, col=icol, block=ds_block, found=found)
233  cpassert(found)
234  fij(idir) = fij(idir) + sum(p_block*ds_block)
235  END DO
236  END DO
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)
243  END IF
244  END DO
245  CALL dbcsr_iterator_stop(iter)
246 
247  IF (calculate_forces) THEN
248  DO ikind = 1, SIZE(atomic_kind_set)
249  CALL para_env%sum(force(ikind)%efield)
250  END DO
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 
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, &
284  nkind, 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, 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(:), &
305  POINTER :: sab_orb
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
311 
312  CALL timeset(routinen, handle)
313 
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)
320 
321  ! calculate stress only if forces requested also
322  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
323  use_virial = use_virial .AND. calculate_forces
324 
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
329  DO idir = 1, 3
330  fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir)
331  END DO
332 
333  natom = SIZE(particle_set)
334  nspin = SIZE(ks_matrix, 1)
335 
336  zi(:) = cmplx(1._dp, 0._dp, dp)
337  DO ia = 1, natom
338  charge = mcharge(ia)
339  ria = particle_set(ia)%r
340  DO idir = 1, 3
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
345  END DO
346  END DO
347  qi = aimag(log(zi))
348  energy%efield = -sum(fpolvec(:)*qi(:))
349 
350  IF (.NOT. just_energy) THEN
351  CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
352  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
353 
354  nimg = dft_control%nimages
355  NULLIFY (cell_to_index)
356  IF (nimg > 1) THEN
357  NULLIFY (kpoints)
358  CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
359  CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
360  END IF
361 
362  IF (calculate_forces) THEN
363  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
364  CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
365  IF (para_env%mepos == 0) THEN
366  DO ia = 1, natom
367  charge = -mcharge(ia)
368  iatom = atom_of_kind(ia)
369  ikind = kind_of(ia)
370  force(ikind)%efield(:, iatom) = fieldpol(:)*charge
371  IF (use_virial) THEN
372  ria = particle_set(ia)%r
373  ria = pbc(ria, cell)
374  forcea(1:3) = fieldpol(1:3)*charge
375  CALL virial_pair_force(virial%pv_virial, -0.5_dp, forcea, ria)
376  CALL virial_pair_force(virial%pv_virial, -0.5_dp, ria, forcea)
377  END IF
378  END DO
379  ELSE
380  DO ia = 1, natom
381  iatom = atom_of_kind(ia)
382  ikind = kind_of(ia)
383  force(ikind)%efield(:, iatom) = 0.0_dp
384  END DO
385  END IF
386  END IF
387 
388  IF (nimg == 1) THEN
389  ! no k-points; all matrices have been transformed to periodic bsf
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)
393 
394  fdir = 0.0_dp
395  ria = particle_set(irow)%r
396  rib = particle_set(icol)%r
397  DO idir = 1, 3
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))
405  END DO
406 
407  DO is = 1, nspin
408  NULLIFY (ks_block)
409  CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
410  row=irow, col=icol, block=ks_block, found=found)
411  cpassert(found)
412  ks_block = ks_block + 0.5_dp*fdir*s_block
413  END DO
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)
419  fij = 0.0_dp
420  DO ispin = 1, nspin
421  CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
422  row=irow, col=icol, block=p_block, found=found)
423  cpassert(found)
424  DO idir = 1, 3
425  CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
426  row=irow, col=icol, block=ds_block, found=found)
427  cpassert(found)
428  fij(idir) = fij(idir) + sum(p_block*ds_block)
429  END DO
430  END DO
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)
433  END IF
434  END DO
435  CALL dbcsr_iterator_stop(iter)
436  !
437  ! stress tensor for Gamma point needs to recalculate overlap integral derivatives
438  !
439  IF (use_virial) THEN
440  ! derivative overlap integral (non collapsed)
441  NULLIFY (sap_int)
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
445  CALL xtb_dsint_list(qs_env, sap_int)
446  ELSE
447  cpabort("TB method unknown")
448  END IF
449  !
450  CALL get_qs_env(qs_env, nkind=nkind)
451  DO ikind = 1, nkind
452  DO jkind = 1, nkind
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
467  fdir = 0.0_dp
468  ria = particle_set(irow)%r
469  rib = particle_set(icol)%r
470  DO idir = 1, 3
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))
478  END DO
479  fi = 1.0_dp
480  IF (iatom == jatom) fi = 0.5_dp
481  DO ispin = 1, nspin
482  NULLIFY (p_block)
483  CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
484  row=irow, col=icol, block=p_block, found=found)
485  cpassert(found)
486  fij = 0.0_dp
487  DO idir = 1, 3
488  IF (irow == iatom) THEN
489  fij(idir) = sum(p_block*dsint(:, :, idir))
490  ELSE
491  fij(idir) = sum(transpose(p_block)*dsint(:, :, idir))
492  END IF
493  END DO
494  IF (irow == iatom) fij = -fij
495  CALL virial_pair_force(virial%pv_virial, fi, fdir*fij(1:3), rij)
496  END DO
497  END IF
498  END DO
499  END DO
500  END DO
501  END DO
502  CALL release_sap_int(sap_int)
503  END IF
504  ELSE
505  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
506  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
507  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
508  iatom=iatom, jatom=jatom, r=rab, cell=cellind)
509 
510  icol = max(iatom, jatom)
511  irow = min(iatom, jatom)
512 
513  ic = cell_to_index(cellind(1), cellind(2), cellind(3))
514  cpassert(ic > 0)
515 
516  fdir = 0.0_dp
517  ria = particle_set(irow)%r
518  rib = particle_set(icol)%r
519  DO idir = 1, 3
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))
527  END DO
528 
529  NULLIFY (s_block)
530  CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
531  row=irow, col=icol, block=s_block, found=found)
532  cpassert(found)
533  DO is = 1, nspin
534  NULLIFY (ks_block)
535  CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
536  row=irow, col=icol, block=ks_block, found=found)
537  cpassert(found)
538  ks_block = ks_block + 0.5_dp*fdir*s_block
539  END DO
540  IF (calculate_forces) THEN
541  atom_a = atom_of_kind(iatom)
542  atom_b = atom_of_kind(jatom)
543  fij = 0.0_dp
544  DO ispin = 1, nspin
545  CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
546  row=irow, col=icol, block=p_block, found=found)
547  cpassert(found)
548  DO idir = 1, 3
549  CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
550  row=irow, col=icol, block=ds_block, found=found)
551  cpassert(found)
552  fij(idir) = fij(idir) + sum(p_block*ds_block)
553  END DO
554  END DO
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)
558  IF (use_virial) THEN
559  dr = sqrt(sum(rab(:)**2))
560  IF (dr > 1.e-6_dp) THEN
561  fi = 1.0_dp
562  IF (iatom == jatom) fi = 0.5_dp
563  CALL virial_pair_force(virial%pv_virial, fi, -fdir*fij(1:3), rab)
564  END IF
565  END IF
566  END IF
567  END DO
568  CALL neighbor_list_iterator_release(nl_iterator)
569  END IF
570 
571  IF (calculate_forces) THEN
572  DO ikind = 1, SIZE(atomic_kind_set)
573  CALL para_env%sum(force(ikind)%efield)
574  END DO
575  END IF
576 
577  END IF
578 
579  CALL timestop(handle)
580 
581  END SUBROUTINE efield_tb_berry
582 
583 ! **************************************************************************************************
584 !> \brief ...
585 !> \param qs_env ...
586 !> \param ks_matrix ...
587 !> \param rho ...
588 !> \param mcharge ...
589 !> \param energy ...
590 !> \param calculate_forces ...
591 !> \param just_energy ...
592 ! **************************************************************************************************
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
600 
601  CHARACTER(LEN=*), PARAMETER :: routinen = 'dfield_tb_berry'
602 
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, &
608  nspin
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(:), &
629  POINTER :: sab_orb
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
634 
635  CALL timeset(routinen, handle)
636 
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)
643 
644  ! efield history
645  CALL init_efield_matrices(efield)
646  CALL set_qs_env(qs_env, efield=efield)
647 
648  ! calculate stress only if forces requested also
649  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
650  use_virial = use_virial .AND. calculate_forces
651  ! disable stress calculation
652  IF (use_virial) THEN
653  cpabort("Stress tensor for periodic D-field not implemented")
654  END IF
655 
656  dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
657 
658  fieldpol = dft_control%period_efield%polarisation
659  fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
660  fieldpol = fieldpol*dft_control%period_efield%strength
661 
662  omega = cell%deth
663  hmat = cell%hmat(:, :)/twopi
664 
665  natom = SIZE(particle_set)
666  nspin = SIZE(ks_matrix, 1)
667 
668  zi(:) = cmplx(1._dp, 0._dp, dp)
669  DO ia = 1, natom
670  charge = mcharge(ia)
671  ria = particle_set(ia)%r
672  DO idir = 1, 3
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
677  END DO
678  END DO
679  qi = aimag(log(zi))
680 
681  ! make sure the total normalized polarization is within [-1:1]
682  DO idir = 1, 3
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
686  ! now check for log branch
687  IF (calculate_forces) THEN
688  IF (abs(efield%polarisation(idir) - cqi(idir)) > pi) THEN
689  di(idir) = (efield%polarisation(idir) - cqi(idir))/pi
690  DO i = 1, 10
691  cqi(idir) = cqi(idir) + sign(1.0_dp, di(idir))*twopi
692  IF (abs(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
693  END DO
694  END IF
695  END IF
696  cqi(idir) = cqi(idir)*omega
697  END DO
698  DO idir = 1, 3
699  ci(idir) = 0.0_dp
700  DO i = 1, 3
701  ci(idir) = ci(idir) + hmat(idir, i)*cqi(i)
702  END DO
703  END DO
704  ! update the references
705  IF (calculate_forces) THEN
706  ener_field = sum(ci)
707  ! check for smoothness of energy surface
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")
710  END IF
711  efield%field_energy = ener_field
712  efield%polarisation(:) = cqi(:)/omega
713  END IF
714 
715  ! Energy
716  ener_field = 0.0_dp
717  DO idir = 1, 3
718  ener_field = ener_field + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi/omega*ci(idir))**2
719  END DO
720  energy%efield = 0.25_dp/twopi*ener_field
721 
722  IF (.NOT. just_energy) THEN
723  di(:) = -(fieldpol(:) - 2._dp*twopi/omega*ci(:))*dfilter(:)/omega
724 
725  CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
726  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
727 
728  nimg = dft_control%nimages
729  NULLIFY (cell_to_index)
730  IF (nimg > 1) THEN
731  NULLIFY (kpoints)
732  CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
733  CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
734  END IF
735 
736  IF (calculate_forces) THEN
737  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
738  CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
739  IF (para_env%mepos == 0) THEN
740  DO ia = 1, natom
741  charge = mcharge(ia)
742  iatom = atom_of_kind(ia)
743  ikind = kind_of(ia)
744  force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + di(:)*charge
745  END DO
746  END IF
747  END IF
748 
749  IF (nimg == 1) THEN
750  ! no k-points; all matrices have been transformed to periodic bsf
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)
754 
755  DO idir = 1, 3
756  hdi(idir) = -sum(di(1:3)*hmat(1:3, idir))
757  END DO
758  fdir = 0.0_dp
759  ria = particle_set(irow)%r
760  rib = particle_set(icol)%r
761  DO idir = 1, 3
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))
769  END DO
770 
771  DO is = 1, nspin
772  NULLIFY (ks_block)
773  CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
774  row=irow, col=icol, block=ks_block, found=found)
775  cpassert(found)
776  ks_block = ks_block + 0.5_dp*fdir*s_block
777  END DO
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)
783  fij = 0.0_dp
784  DO ispin = 1, nspin
785  CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
786  row=irow, col=icol, block=p_block, found=found)
787  cpassert(found)
788  DO idir = 1, 3
789  CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
790  row=irow, col=icol, block=ds_block, found=found)
791  cpassert(found)
792  fij(idir) = fij(idir) + sum(p_block*ds_block)
793  END DO
794  END DO
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)
797  END IF
798 
799  END DO
800  CALL dbcsr_iterator_stop(iter)
801  ELSE
802  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
803  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
804  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
805  iatom=iatom, jatom=jatom, r=rab, cell=cellind)
806 
807  icol = max(iatom, jatom)
808  irow = min(iatom, jatom)
809 
810  ic = cell_to_index(cellind(1), cellind(2), cellind(3))
811  cpassert(ic > 0)
812 
813  DO idir = 1, 3
814  hdi(idir) = -sum(di(1:3)*hmat(1:3, idir))
815  END DO
816  fdir = 0.0_dp
817  ria = particle_set(irow)%r
818  rib = particle_set(icol)%r
819  DO idir = 1, 3
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))
827  END DO
828 
829  NULLIFY (s_block)
830  CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
831  row=irow, col=icol, block=s_block, found=found)
832  cpassert(found)
833  DO is = 1, nspin
834  NULLIFY (ks_block)
835  CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
836  row=irow, col=icol, block=ks_block, found=found)
837  cpassert(found)
838  ks_block = ks_block + 0.5_dp*fdir*s_block
839  END DO
840  IF (calculate_forces) THEN
841  atom_a = atom_of_kind(iatom)
842  atom_b = atom_of_kind(jatom)
843  fij = 0.0_dp
844  DO ispin = 1, nspin
845  CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
846  row=irow, col=icol, block=p_block, found=found)
847  cpassert(found)
848  DO idir = 1, 3
849  CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
850  row=irow, col=icol, block=ds_block, found=found)
851  cpassert(found)
852  fij(idir) = fij(idir) + sum(p_block*ds_block)
853  END DO
854  END DO
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)
858  END IF
859 
860  END DO
861  CALL neighbor_list_iterator_release(nl_iterator)
862  END IF
863 
864  END IF
865 
866  CALL timestop(handle)
867 
868  END SUBROUTINE dfield_tb_berry
869 
870 END MODULE efield_tb_methods
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
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.
Definition: cell_types.F:15
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.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
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: kpoint_types.F:333
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
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.
Definition: qs_kind_types.F:23
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...
Definition: qs_rho_types.F:18
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...
Definition: qs_rho_types.F:229
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.
Definition: xtb_coulomb.F:12
subroutine, public xtb_dsint_list(qs_env, sap_int)
...
Definition: xtb_coulomb.F:820