(git:374b731)
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-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! **************************************************************************************************
15 USE cell_types, ONLY: cell_type,&
16 pbc
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,&
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 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
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
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
870END 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...
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 set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
type for berry phase efield matrices. At the moment only used for cosmat and sinmat
subroutine, public init_efield_matrices(efield)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
General overlap type integrals containers.
subroutine, public release_sap_int(sap_int)
...
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Calculation of Coulomb contributions in xTB.
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.