(git:ed6f26b)
Loading...
Searching...
No Matches
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-2025 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! **************************************************************************************************
15 USE cell_types, ONLY: cell_type,&
16 pbc
25 USE kinds, ONLY: dp
26 USE kpoint_types, ONLY: get_kpoint_info,&
28 USE mathconstants, ONLY: pi,&
29 twopi
46 USE qs_rho_types, ONLY: qs_rho_get,&
51 USE virial_types, ONLY: virial_type
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
63CONTAINS
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 ! check if the periodic efield should be applied in the current step
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
98
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)
101 ELSE
102 CALL efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
103 END IF
104 END IF
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")
109 END IF
110 ELSE
111 cpabort("This routine should only be called from TB")
112 END IF
113
114 CALL timestop(handle)
115
116 END SUBROUTINE efield_tb_matrix
117
118! **************************************************************************************************
119!> \brief ...
120!> \param qs_env ...
121!> \param ks_matrix ...
122!> \param rho ...
123!> \param mcharge ...
124!> \param energy ...
125!> \param calculate_forces ...
126!> \param just_energy ...
127! **************************************************************************************************
128 SUBROUTINE efield_tb_local(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
129 TYPE(qs_environment_type), POINTER :: qs_env
130 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
131 TYPE(qs_rho_type), POINTER :: rho
132 REAL(dp), DIMENSION(:), INTENT(in) :: mcharge
133 TYPE(qs_energy_type), POINTER :: energy
134 LOGICAL, INTENT(in) :: calculate_forces, just_energy
135
136 CHARACTER(LEN=*), PARAMETER :: routinen = 'efield_tb_local'
137
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
145 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
146 TYPE(cell_type), POINTER :: cell
147 TYPE(dbcsr_iterator_type) :: iter
148 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s
149 TYPE(dft_control_type), POINTER :: dft_control
150 TYPE(mp_para_env_type), POINTER :: para_env
151 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
152 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
153 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
154 TYPE(virial_type), POINTER :: virial
155
156 CALL timeset(routinen, handle)
157
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)
161 IF (do_kpoints) THEN
162 cpabort("Local electric field with kpoints not possible. Use Berry phase periodic version")
163 END IF
164 ! disable stress calculation
165 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
166 IF (use_virial) THEN
167 cpabort("Stress tensor for non-periodic E-field not possible")
168 END IF
169
170 fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
171 dft_control%efield_fields(1)%efield%strength
172
173 natom = SIZE(particle_set)
174 ci = 0.0_dp
175 DO ia = 1, natom
176 charge = mcharge(ia)
177 ria = particle_set(ia)%r
178 ria = pbc(ria, cell)
179 ci(:) = ci(:) + charge*ria(:)
180 END DO
181 energy%efield = -sum(ci(:)*fieldpol(:))
182
183 IF (.NOT. just_energy) THEN
184
185 IF (calculate_forces) THEN
186 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
187 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
188 IF (para_env%mepos == 0) THEN
189 DO ia = 1, natom
190 charge = mcharge(ia)
191 ikind = kind_of(ia)
192 atom_a = atom_of_kind(ia)
193 force(ikind)%efield(1:3, atom_a) = -charge*fieldpol(:)
194 END DO
195 ELSE
196 DO ia = 1, natom
197 ikind = kind_of(ia)
198 atom_a = atom_of_kind(ia)
199 force(ikind)%efield(1:3, atom_a) = 0.0_dp
200 END DO
201 END IF
202 CALL qs_rho_get(rho, rho_ao=matrix_p)
203 END IF
204
205 ! Update KS matrix
206 nspin = SIZE(ks_matrix, 1)
207 NULLIFY (matrix_s)
208 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
209 CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
210 DO WHILE (dbcsr_iterator_blocks_left(iter))
211 NULLIFY (ks_block, s_block, p_block)
212 CALL dbcsr_iterator_next_block(iter, irow, icol, s_block)
213 ria = particle_set(irow)%r
214 ria = pbc(ria, cell)
215 rib = particle_set(icol)%r
216 rib = pbc(rib, cell)
217 fdir = 0.5_dp*sum(fieldpol(1:3)*(ria(1:3) + rib(1:3)))
218 DO ispin = 1, nspin
219 CALL dbcsr_get_block_p(matrix=ks_matrix(ispin, 1)%matrix, &
220 row=irow, col=icol, block=ks_block, found=found)
221 ks_block = ks_block + fdir*s_block
222 cpassert(found)
223 END DO
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)
229 fij = 0.0_dp
230 DO ispin = 1, nspin
231 CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
232 row=irow, col=icol, block=p_block, found=found)
233 cpassert(found)
234 DO idir = 1, 3
235 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1)%matrix, &
236 row=irow, col=icol, block=ds_block, found=found)
237 cpassert(found)
238 fij(idir) = fij(idir) + sum(p_block*ds_block)
239 END DO
240 END DO
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)
247 END IF
248 END DO
249 CALL dbcsr_iterator_stop(iter)
250
251 IF (calculate_forces) THEN
252 DO ikind = 1, SIZE(atomic_kind_set)
253 CALL para_env%sum(force(ikind)%efield)
254 END DO
255 END IF
256
257 END IF
258
259 CALL timestop(handle)
260
261 END SUBROUTINE efield_tb_local
262
263! **************************************************************************************************
264!> \brief ...
265!> \param qs_env ...
266!> \param ks_matrix ...
267!> \param rho ...
268!> \param mcharge ...
269!> \param energy ...
270!> \param calculate_forces ...
271!> \param just_energy ...
272! **************************************************************************************************
273 SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
274 TYPE(qs_environment_type), POINTER :: qs_env
275 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
276 TYPE(qs_rho_type), POINTER :: rho
277 REAL(dp), DIMENSION(:), INTENT(in) :: mcharge
278 TYPE(qs_energy_type), POINTER :: energy
279 LOGICAL, INTENT(in) :: calculate_forces, just_energy
280
281 CHARACTER(LEN=*), PARAMETER :: routinen = 'efield_tb_berry'
282
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, &
288 nkind, nspin
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
299 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
300 TYPE(cell_type), POINTER :: cell
301 TYPE(dbcsr_iterator_type) :: iter
302 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
303 TYPE(dft_control_type), POINTER :: dft_control
304 TYPE(kpoint_type), POINTER :: kpoints
305 TYPE(mp_para_env_type), POINTER :: para_env
307 DIMENSION(:), POINTER :: nl_iterator
308 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
309 POINTER :: sab_orb
310 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
311 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
312 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
313 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
314 TYPE(virial_type), POINTER :: virial
315
316 CALL timeset(routinen, handle)
317
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)
324
325 ! calculate stress only if forces requested also
326 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
327 use_virial = use_virial .AND. calculate_forces
328
329 ! if an intensities list is given, select the value for the current step
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)
334 END IF
335
336 fieldpol = dft_control%period_efield%polarisation
337 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
338 fieldpol = -fieldpol*strength
339 hmat = cell%hmat(:, :)/twopi
340 DO idir = 1, 3
341 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir)
342 END DO
343
344 natom = SIZE(particle_set)
345 nspin = SIZE(ks_matrix, 1)
346
347 zi(:) = cmplx(1._dp, 0._dp, dp)
348 DO ia = 1, natom
349 charge = mcharge(ia)
350 ria = particle_set(ia)%r
351 DO idir = 1, 3
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
356 END DO
357 END DO
358 qi = aimag(log(zi))
359 energy%efield = sum(fpolvec(:)*qi(:))
360
361 IF (.NOT. just_energy) THEN
362 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
363 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
364
365 nimg = dft_control%nimages
366 NULLIFY (cell_to_index)
367 IF (nimg > 1) THEN
368 NULLIFY (kpoints)
369 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
370 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
371 END IF
372
373 IF (calculate_forces) THEN
374 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
375 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
376 IF (para_env%mepos == 0) THEN
377 DO ia = 1, natom
378 charge = mcharge(ia)
379 iatom = atom_of_kind(ia)
380 ikind = kind_of(ia)
381 force(ikind)%efield(:, iatom) = fieldpol(:)*charge
382 IF (use_virial) THEN
383 ria = particle_set(ia)%r
384 ria = pbc(ria, cell)
385 forcea(1:3) = fieldpol(1:3)*charge
386 CALL virial_pair_force(virial%pv_virial, -0.5_dp, forcea, ria)
387 CALL virial_pair_force(virial%pv_virial, -0.5_dp, ria, forcea)
388 END IF
389 END DO
390 ELSE
391 DO ia = 1, natom
392 iatom = atom_of_kind(ia)
393 ikind = kind_of(ia)
394 force(ikind)%efield(:, iatom) = 0.0_dp
395 END DO
396 END IF
397 END IF
398
399 IF (nimg == 1) THEN
400 ! no k-points; all matrices have been transformed to periodic bsf
401 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
402 DO WHILE (dbcsr_iterator_blocks_left(iter))
403 CALL dbcsr_iterator_next_block(iter, irow, icol, s_block)
404
405 fdir = 0.0_dp
406 ria = particle_set(irow)%r
407 rib = particle_set(icol)%r
408 DO idir = 1, 3
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))
416 END DO
417
418 DO is = 1, nspin
419 NULLIFY (ks_block)
420 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
421 row=irow, col=icol, block=ks_block, found=found)
422 cpassert(found)
423 ks_block = ks_block - 0.5_dp*fdir*s_block
424 END DO
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)
430 fij = 0.0_dp
431 DO ispin = 1, nspin
432 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
433 row=irow, col=icol, block=p_block, found=found)
434 cpassert(found)
435 DO idir = 1, 3
436 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
437 row=irow, col=icol, block=ds_block, found=found)
438 cpassert(found)
439 fij(idir) = fij(idir) - sum(p_block*ds_block)
440 END DO
441 END DO
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)
444 END IF
445 END DO
446 CALL dbcsr_iterator_stop(iter)
447 !
448 ! stress tensor for Gamma point needs to recalculate overlap integral derivatives
449 !
450 IF (use_virial) THEN
451 ! derivative overlap integral (non collapsed)
452 NULLIFY (sap_int)
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
456 CALL xtb_dsint_list(qs_env, sap_int)
457 ELSE
458 cpabort("TB method unknown")
459 END IF
460 !
461 CALL get_qs_env(qs_env, nkind=nkind)
462 DO ikind = 1, nkind
463 DO jkind = 1, nkind
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
478 fdir = 0.0_dp
479 ria = particle_set(irow)%r
480 rib = particle_set(icol)%r
481 DO idir = 1, 3
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))
489 END DO
490 fi = 1.0_dp
491 IF (iatom == jatom) fi = 0.5_dp
492 DO ispin = 1, nspin
493 NULLIFY (p_block)
494 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
495 row=irow, col=icol, block=p_block, found=found)
496 cpassert(found)
497 fij = 0.0_dp
498 DO idir = 1, 3
499 IF (irow == iatom) THEN
500 fij(idir) = sum(p_block*dsint(:, :, idir))
501 ELSE
502 fij(idir) = sum(transpose(p_block)*dsint(:, :, idir))
503 END IF
504 END DO
505 IF (irow == iatom) fij = -fij
506 CALL virial_pair_force(virial%pv_virial, fi, fdir*fij(1:3), rij)
507 END DO
508 END IF
509 END DO
510 END DO
511 END DO
512 END DO
513 CALL release_sap_int(sap_int)
514 END IF
515 ELSE
516 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
517 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
518 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
519 iatom=iatom, jatom=jatom, r=rab, cell=cellind)
520
521 icol = max(iatom, jatom)
522 irow = min(iatom, jatom)
523
524 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
525 cpassert(ic > 0)
526
527 fdir = 0.0_dp
528 ria = particle_set(irow)%r
529 rib = particle_set(icol)%r
530 DO idir = 1, 3
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))
538 END DO
539
540 NULLIFY (s_block)
541 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
542 row=irow, col=icol, block=s_block, found=found)
543 cpassert(found)
544 DO is = 1, nspin
545 NULLIFY (ks_block)
546 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
547 row=irow, col=icol, block=ks_block, found=found)
548 cpassert(found)
549 ks_block = ks_block - 0.5_dp*fdir*s_block
550 END DO
551 IF (calculate_forces) THEN
552 atom_a = atom_of_kind(iatom)
553 atom_b = atom_of_kind(jatom)
554 fij = 0.0_dp
555 DO ispin = 1, nspin
556 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
557 row=irow, col=icol, block=p_block, found=found)
558 cpassert(found)
559 DO idir = 1, 3
560 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
561 row=irow, col=icol, block=ds_block, found=found)
562 cpassert(found)
563 fij(idir) = fij(idir) - sum(p_block*ds_block)
564 END DO
565 END DO
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)
569 IF (use_virial) THEN
570 dr = sqrt(sum(rab(:)**2))
571 IF (dr > 1.e-6_dp) THEN
572 fi = 1.0_dp
573 IF (iatom == jatom) fi = 0.5_dp
574 CALL virial_pair_force(virial%pv_virial, fi, -fdir*fij(1:3), rab)
575 END IF
576 END IF
577 END IF
578 END DO
579 CALL neighbor_list_iterator_release(nl_iterator)
580 END IF
581
582 IF (calculate_forces) THEN
583 DO ikind = 1, SIZE(atomic_kind_set)
584 CALL para_env%sum(force(ikind)%efield)
585 END DO
586 END IF
587
588 END IF
589
590 CALL timestop(handle)
591
592 END SUBROUTINE efield_tb_berry
593
594! **************************************************************************************************
595!> \brief ...
596!> \param qs_env ...
597!> \param ks_matrix ...
598!> \param rho ...
599!> \param mcharge ...
600!> \param energy ...
601!> \param calculate_forces ...
602!> \param just_energy ...
603! **************************************************************************************************
604 SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
605 TYPE(qs_environment_type), POINTER :: qs_env
606 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
607 TYPE(qs_rho_type), POINTER :: rho
608 REAL(dp), DIMENSION(:), INTENT(in) :: mcharge
609 TYPE(qs_energy_type), POINTER :: energy
610 LOGICAL, INTENT(in) :: calculate_forces, just_energy
611
612 CHARACTER(LEN=*), PARAMETER :: routinen = 'dfield_tb_berry'
613
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, &
624 strength
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
629 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
630 TYPE(cell_type), POINTER :: cell
631 TYPE(dbcsr_iterator_type) :: iter
632 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
633 TYPE(dft_control_type), POINTER :: dft_control
634 TYPE(efield_berry_type), POINTER :: efield
635 TYPE(kpoint_type), POINTER :: kpoints
636 TYPE(mp_para_env_type), POINTER :: para_env
638 DIMENSION(:), POINTER :: nl_iterator
639 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
640 POINTER :: sab_orb
641 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
642 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
643 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
644 TYPE(virial_type), POINTER :: virial
645
646 CALL timeset(routinen, handle)
647
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)
654
655 ! efield history
656 CALL init_efield_matrices(efield)
657 CALL set_qs_env(qs_env, efield=efield)
658
659 ! calculate stress only if forces requested also
660 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
661 use_virial = use_virial .AND. calculate_forces
662 ! disable stress calculation
663 IF (use_virial) THEN
664 cpabort("Stress tensor for periodic D-field not implemented")
665 END IF
666
667 dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
668
669 ! if an intensities list is given, select the value for the current step
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)
674 END IF
675
676 fieldpol = dft_control%period_efield%polarisation
677 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
678 fieldpol = fieldpol*strength
679
680 omega = cell%deth
681 hmat = cell%hmat(:, :)/twopi
682
683 natom = SIZE(particle_set)
684 nspin = SIZE(ks_matrix, 1)
685
686 zi(:) = cmplx(1._dp, 0._dp, dp)
687 DO ia = 1, natom
688 charge = mcharge(ia)
689 ria = particle_set(ia)%r
690 DO idir = 1, 3
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
695 END DO
696 END DO
697 qi = aimag(log(zi))
698
699 ! make sure the total normalized polarization is within [-1:1]
700 DO idir = 1, 3
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
704 ! now check for log branch
705 IF (calculate_forces) THEN
706 IF (abs(efield%polarisation(idir) - cqi(idir)) > pi) THEN
707 di(idir) = (efield%polarisation(idir) - cqi(idir))/pi
708 DO i = 1, 10
709 cqi(idir) = cqi(idir) + sign(1.0_dp, di(idir))*twopi
710 IF (abs(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
711 END DO
712 END IF
713 END IF
714 cqi(idir) = cqi(idir)*omega
715 END DO
716 DO idir = 1, 3
717 ci(idir) = 0.0_dp
718 DO i = 1, 3
719 ci(idir) = ci(idir) + hmat(idir, i)*cqi(i)
720 END DO
721 END DO
722 ! update the references
723 IF (calculate_forces) THEN
724 ener_field = sum(ci)
725 ! check for smoothness of energy surface
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")
728 END IF
729 efield%field_energy = ener_field
730 efield%polarisation(:) = cqi(:)/omega
731 END IF
732
733 ! Energy
734 ener_field = 0.0_dp
735 DO idir = 1, 3
736 ener_field = ener_field + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi/omega*ci(idir))**2
737 END DO
738 energy%efield = 0.25_dp/twopi*ener_field
739
740 IF (.NOT. just_energy) THEN
741 di(:) = -(fieldpol(:) - 2._dp*twopi/omega*ci(:))*dfilter(:)/omega
742
743 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
744 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
745
746 nimg = dft_control%nimages
747 NULLIFY (cell_to_index)
748 IF (nimg > 1) THEN
749 NULLIFY (kpoints)
750 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
751 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
752 END IF
753
754 IF (calculate_forces) THEN
755 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
756 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
757 IF (para_env%mepos == 0) THEN
758 DO ia = 1, natom
759 charge = mcharge(ia)
760 iatom = atom_of_kind(ia)
761 ikind = kind_of(ia)
762 force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + di(:)*charge
763 END DO
764 END IF
765 END IF
766
767 IF (nimg == 1) THEN
768 ! no k-points; all matrices have been transformed to periodic bsf
769 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
770 DO WHILE (dbcsr_iterator_blocks_left(iter))
771 CALL dbcsr_iterator_next_block(iter, irow, icol, s_block)
772
773 DO idir = 1, 3
774 hdi(idir) = -sum(di(1:3)*hmat(1:3, idir))
775 END DO
776 fdir = 0.0_dp
777 ria = particle_set(irow)%r
778 rib = particle_set(icol)%r
779 DO idir = 1, 3
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))
787 END DO
788
789 DO is = 1, nspin
790 NULLIFY (ks_block)
791 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
792 row=irow, col=icol, block=ks_block, found=found)
793 cpassert(found)
794 ks_block = ks_block + 0.5_dp*fdir*s_block
795 END DO
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)
801 fij = 0.0_dp
802 DO ispin = 1, nspin
803 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
804 row=irow, col=icol, block=p_block, found=found)
805 cpassert(found)
806 DO idir = 1, 3
807 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
808 row=irow, col=icol, block=ds_block, found=found)
809 cpassert(found)
810 fij(idir) = fij(idir) + sum(p_block*ds_block)
811 END DO
812 END DO
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)
815 END IF
816
817 END DO
818 CALL dbcsr_iterator_stop(iter)
819 ELSE
820 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
821 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
822 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
823 iatom=iatom, jatom=jatom, r=rab, cell=cellind)
824
825 icol = max(iatom, jatom)
826 irow = min(iatom, jatom)
827
828 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
829 cpassert(ic > 0)
830
831 DO idir = 1, 3
832 hdi(idir) = -sum(di(1:3)*hmat(1:3, idir))
833 END DO
834 fdir = 0.0_dp
835 ria = particle_set(irow)%r
836 rib = particle_set(icol)%r
837 DO idir = 1, 3
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))
845 END DO
846
847 NULLIFY (s_block)
848 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
849 row=irow, col=icol, block=s_block, found=found)
850 cpassert(found)
851 DO is = 1, nspin
852 NULLIFY (ks_block)
853 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
854 row=irow, col=icol, block=ks_block, found=found)
855 cpassert(found)
856 ks_block = ks_block + 0.5_dp*fdir*s_block
857 END DO
858 IF (calculate_forces) THEN
859 atom_a = atom_of_kind(iatom)
860 atom_b = atom_of_kind(jatom)
861 fij = 0.0_dp
862 DO ispin = 1, nspin
863 CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
864 row=irow, col=icol, block=p_block, found=found)
865 cpassert(found)
866 DO idir = 1, 3
867 CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
868 row=irow, col=icol, block=ds_block, found=found)
869 cpassert(found)
870 fij(idir) = fij(idir) + sum(p_block*ds_block)
871 END DO
872 END DO
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)
876 END IF
877
878 END DO
879 CALL neighbor_list_iterator_release(nl_iterator)
880 END IF
881
882 END IF
883
884 CALL timestop(handle)
885
886 END SUBROUTINE dfield_tb_berry
887
888END MODULE efield_tb_methods
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...
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.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
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.
Definition xtb_coulomb.F:12
subroutine, public xtb_dsint_list(qs_env, sap_int)
...
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
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.