Line data Source code
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 : ! **************************************************************************************************
12 : MODULE efield_tb_methods
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind_set
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,&
35 : set_qs_env
36 : USE qs_force_types, ONLY: qs_force_type
37 : USE qs_kind_types, ONLY: qs_kind_type
38 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
39 : neighbor_list_iterate,&
40 : neighbor_list_iterator_create,&
41 : neighbor_list_iterator_p_type,&
42 : neighbor_list_iterator_release,&
43 : neighbor_list_set_p_type
44 : USE qs_period_efield_types, ONLY: efield_berry_type,&
45 : init_efield_matrices
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
50 : USE virial_methods, ONLY: virial_pair_force
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 17810 : 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 17810 : CALL timeset(routineN, handle)
90 :
91 17810 : energy%efield = 0.0_dp
92 17810 : CALL get_qs_env(qs_env, dft_control=dft_control)
93 17810 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
94 17810 : IF (dft_control%apply_period_efield) THEN
95 4292 : IF (dft_control%period_efield%displacement_field) THEN
96 332 : CALL dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
97 : ELSE
98 3960 : CALL efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
99 : END IF
100 13518 : ELSE IF (dft_control%apply_efield) THEN
101 2176 : CALL efield_tb_local(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
102 11342 : ELSE IF (dft_control%apply_efield_field) THEN
103 0 : CPABORT("efield_filed")
104 : END IF
105 : ELSE
106 0 : CPABORT("This routine should only be called from TB")
107 : END IF
108 :
109 17810 : CALL timestop(handle)
110 :
111 17810 : 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 2176 : 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 2176 : 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 2176 : REAL(dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block
141 2176 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
142 : TYPE(cell_type), POINTER :: cell
143 : TYPE(dbcsr_iterator_type) :: iter
144 2176 : 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 2176 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
148 2176 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
149 2176 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
150 : TYPE(virial_type), POINTER :: virial
151 :
152 2176 : CALL timeset(routineN, handle)
153 :
154 2176 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, particle_set=particle_set)
155 2176 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, energy=energy, para_env=para_env)
156 2176 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints, virial=virial)
157 2176 : IF (do_kpoints) THEN
158 0 : CPABORT("Local electric field with kpoints not possible. Use Berry phase periodic version")
159 : END IF
160 : ! disable stress calculation
161 2176 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
162 : IF (use_virial) THEN
163 0 : CPABORT("Stress tensor for non-periodic E-field not possible")
164 : END IF
165 :
166 : fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
167 8704 : dft_control%efield_fields(1)%efield%strength
168 :
169 2176 : natom = SIZE(particle_set)
170 2176 : ci = 0.0_dp
171 10556 : DO ia = 1, natom
172 8380 : charge = mcharge(ia)
173 33520 : ria = particle_set(ia)%r
174 33520 : ria = pbc(ria, cell)
175 35696 : ci(:) = ci(:) + charge*ria(:)
176 : END DO
177 8704 : energy%efield = -SUM(ci(:)*fieldpol(:))
178 :
179 2176 : IF (.NOT. just_energy) THEN
180 :
181 2176 : IF (calculate_forces) THEN
182 32 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
183 32 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
184 32 : IF (para_env%mepos == 0) THEN
185 72 : DO ia = 1, natom
186 56 : charge = mcharge(ia)
187 56 : ikind = kind_of(ia)
188 56 : atom_a = atom_of_kind(ia)
189 240 : force(ikind)%efield(1:3, atom_a) = -charge*fieldpol(:)
190 : END DO
191 : ELSE
192 72 : DO ia = 1, natom
193 56 : ikind = kind_of(ia)
194 56 : atom_a = atom_of_kind(ia)
195 240 : force(ikind)%efield(1:3, atom_a) = 0.0_dp
196 : END DO
197 : END IF
198 32 : CALL qs_rho_get(rho, rho_ao=matrix_p)
199 : END IF
200 :
201 : ! Update KS matrix
202 2176 : nspin = SIZE(ks_matrix, 1)
203 2176 : NULLIFY (matrix_s)
204 2176 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
205 2176 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
206 12677 : DO WHILE (dbcsr_iterator_blocks_left(iter))
207 10501 : NULLIFY (ks_block, s_block, p_block)
208 10501 : CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)
209 42004 : ria = particle_set(irow)%r
210 42004 : ria = pbc(ria, cell)
211 42004 : rib = particle_set(icol)%r
212 42004 : rib = pbc(rib, cell)
213 42004 : fdir = 0.5_dp*SUM(fieldpol(1:3)*(ria(1:3) + rib(1:3)))
214 21002 : DO ispin = 1, nspin
215 : CALL dbcsr_get_block_p(matrix=ks_matrix(ispin, 1)%matrix, &
216 10501 : row=irow, col=icol, BLOCK=ks_block, found=found)
217 274417 : ks_block = ks_block + fdir*s_block
218 31503 : CPASSERT(found)
219 : END DO
220 12677 : IF (calculate_forces) THEN
221 136 : ikind = kind_of(irow)
222 136 : jkind = kind_of(icol)
223 136 : atom_a = atom_of_kind(irow)
224 136 : atom_b = atom_of_kind(icol)
225 136 : fij = 0.0_dp
226 272 : DO ispin = 1, nspin
227 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
228 136 : row=irow, col=icol, BLOCK=p_block, found=found)
229 136 : CPASSERT(found)
230 816 : DO idir = 1, 3
231 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1)%matrix, &
232 408 : row=irow, col=icol, BLOCK=ds_block, found=found)
233 408 : CPASSERT(found)
234 6502 : fij(idir) = fij(idir) + SUM(p_block*ds_block)
235 : END DO
236 : END DO
237 544 : fdir = SUM(ria(1:3)*fieldpol(1:3))
238 544 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
239 544 : force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
240 544 : fdir = SUM(rib(1:3)*fieldpol(1:3))
241 544 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
242 544 : 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 2176 : CALL dbcsr_iterator_stop(iter)
246 :
247 2176 : IF (calculate_forces) THEN
248 120 : DO ikind = 1, SIZE(atomic_kind_set)
249 1016 : CALL para_env%sum(force(ikind)%efield)
250 : END DO
251 : END IF
252 :
253 : END IF
254 :
255 2176 : CALL timestop(handle)
256 :
257 4352 : 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 3960 : 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 3960 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
286 : INTEGER, DIMENSION(3) :: cellind
287 3960 : 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 3960 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block
294 3960 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
295 3960 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
296 : TYPE(cell_type), POINTER :: cell
297 : TYPE(dbcsr_iterator_type) :: iter
298 3960 : 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 3960 : DIMENSION(:), POINTER :: nl_iterator
304 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
305 3960 : POINTER :: sab_orb
306 3960 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
307 3960 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
308 3960 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
309 3960 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
310 : TYPE(virial_type), POINTER :: virial
311 :
312 3960 : CALL timeset(routineN, handle)
313 :
314 3960 : NULLIFY (dft_control, cell, particle_set)
315 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
316 3960 : particle_set=particle_set, virial=virial)
317 3960 : NULLIFY (qs_kind_set, para_env, sab_orb)
318 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
319 3960 : energy=energy, para_env=para_env, sab_orb=sab_orb)
320 :
321 : ! calculate stress only if forces requested also
322 4076 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
323 124 : use_virial = use_virial .AND. calculate_forces
324 :
325 15840 : fieldpol = dft_control%period_efield%polarisation
326 27720 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
327 15840 : fieldpol = -fieldpol*dft_control%period_efield%strength
328 51480 : hmat = cell%hmat(:, :)/twopi
329 15840 : DO idir = 1, 3
330 15840 : fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir)
331 : END DO
332 :
333 3960 : natom = SIZE(particle_set)
334 3960 : nspin = SIZE(ks_matrix, 1)
335 :
336 15840 : zi(:) = CMPLX(1._dp, 0._dp, dp)
337 18474 : DO ia = 1, natom
338 14514 : charge = mcharge(ia)
339 58056 : ria = particle_set(ia)%r
340 62016 : DO idir = 1, 3
341 174168 : kvec(:) = twopi*cell%h_inv(idir, :)
342 174168 : dd = SUM(kvec(:)*ria(:))
343 43542 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
344 58056 : zi(idir) = zi(idir)*zdeta
345 : END DO
346 : END DO
347 15840 : qi = AIMAG(LOG(zi))
348 15840 : energy%efield = -SUM(fpolvec(:)*qi(:))
349 :
350 3960 : IF (.NOT. just_energy) THEN
351 3960 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
352 3960 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
353 :
354 3960 : nimg = dft_control%nimages
355 3960 : NULLIFY (cell_to_index)
356 3960 : IF (nimg > 1) THEN
357 1858 : NULLIFY (kpoints)
358 1858 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
359 1858 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
360 : END IF
361 :
362 3960 : IF (calculate_forces) THEN
363 46 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
364 46 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
365 46 : IF (para_env%mepos == 0) THEN
366 94 : DO ia = 1, natom
367 71 : charge = -mcharge(ia)
368 71 : iatom = atom_of_kind(ia)
369 71 : ikind = kind_of(ia)
370 284 : force(ikind)%efield(:, iatom) = fieldpol(:)*charge
371 94 : IF (use_virial) THEN
372 64 : ria = particle_set(ia)%r
373 64 : ria = pbc(ria, cell)
374 64 : forcea(1:3) = fieldpol(1:3)*charge
375 16 : CALL virial_pair_force(virial%pv_virial, -0.5_dp, forcea, ria)
376 16 : CALL virial_pair_force(virial%pv_virial, -0.5_dp, ria, forcea)
377 : END IF
378 : END DO
379 : ELSE
380 94 : DO ia = 1, natom
381 71 : iatom = atom_of_kind(ia)
382 71 : ikind = kind_of(ia)
383 307 : force(ikind)%efield(:, iatom) = 0.0_dp
384 : END DO
385 : END IF
386 : END IF
387 :
388 3960 : IF (nimg == 1) THEN
389 : ! no k-points; all matrices have been transformed to periodic bsf
390 2102 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
391 11294 : DO WHILE (dbcsr_iterator_blocks_left(iter))
392 9192 : CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)
393 :
394 9192 : fdir = 0.0_dp
395 36768 : ria = particle_set(irow)%r
396 36768 : rib = particle_set(icol)%r
397 36768 : DO idir = 1, 3
398 110304 : kvec(:) = twopi*cell%h_inv(idir, :)
399 110304 : dd = SUM(kvec(:)*ria(:))
400 27576 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
401 27576 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
402 110304 : dd = SUM(kvec(:)*rib(:))
403 27576 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
404 36768 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
405 : END DO
406 :
407 22756 : DO is = 1, nspin
408 13564 : NULLIFY (ks_block)
409 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
410 13564 : row=irow, col=icol, block=ks_block, found=found)
411 13564 : CPASSERT(found)
412 355608 : ks_block = ks_block + 0.5_dp*fdir*s_block
413 : END DO
414 11294 : IF (calculate_forces) THEN
415 82 : ikind = kind_of(irow)
416 82 : jkind = kind_of(icol)
417 82 : atom_a = atom_of_kind(irow)
418 82 : atom_b = atom_of_kind(icol)
419 82 : fij = 0.0_dp
420 208 : DO ispin = 1, nspin
421 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
422 126 : row=irow, col=icol, BLOCK=p_block, found=found)
423 126 : CPASSERT(found)
424 712 : DO idir = 1, 3
425 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
426 378 : row=irow, col=icol, BLOCK=ds_block, found=found)
427 378 : CPASSERT(found)
428 5076 : fij(idir) = fij(idir) + SUM(p_block*ds_block)
429 : END DO
430 : END DO
431 328 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
432 328 : 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 2102 : CALL dbcsr_iterator_stop(iter)
436 : !
437 : ! stress tensor for Gamma point needs to recalculate overlap integral derivatives
438 : !
439 2102 : IF (use_virial) THEN
440 : ! derivative overlap integral (non collapsed)
441 4 : NULLIFY (sap_int)
442 4 : IF (dft_control%qs_control%dftb) THEN
443 0 : CPABORT("DFTB stress tensor for periodic efield not implemented")
444 4 : ELSEIF (dft_control%qs_control%xtb) THEN
445 4 : CALL xtb_dsint_list(qs_env, sap_int)
446 : ELSE
447 0 : CPABORT("TB method unknown")
448 : END IF
449 : !
450 4 : CALL get_qs_env(qs_env, nkind=nkind)
451 16 : DO ikind = 1, nkind
452 52 : DO jkind = 1, nkind
453 36 : iac = ikind + nkind*(jkind - 1)
454 36 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
455 48 : DO ia = 1, sap_int(iac)%nalist
456 18 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE
457 18 : iatom = sap_int(iac)%alist(ia)%aatom
458 1394 : DO ic = 1, sap_int(iac)%alist(ia)%nclist
459 1340 : jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
460 5360 : rij = sap_int(iac)%alist(ia)%clist(ic)%rac
461 5360 : dr = SQRT(SUM(rij(:)**2))
462 1358 : IF (dr > 1.e-6_dp) THEN
463 1332 : dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
464 1332 : icol = MAX(iatom, jatom)
465 1332 : irow = MIN(iatom, jatom)
466 3636 : IF (irow == iatom) rij = -rij
467 1332 : fdir = 0.0_dp
468 5328 : ria = particle_set(irow)%r
469 5328 : rib = particle_set(icol)%r
470 5328 : DO idir = 1, 3
471 15984 : kvec(:) = twopi*cell%h_inv(idir, :)
472 15984 : dd = SUM(kvec(:)*ria(:))
473 3996 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
474 3996 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
475 15984 : dd = SUM(kvec(:)*rib(:))
476 3996 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
477 5328 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
478 : END DO
479 1332 : fi = 1.0_dp
480 1332 : IF (iatom == jatom) fi = 0.5_dp
481 3330 : DO ispin = 1, nspin
482 1998 : NULLIFY (p_block)
483 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
484 1998 : row=irow, col=icol, block=p_block, found=found)
485 1998 : CPASSERT(found)
486 1998 : fij = 0.0_dp
487 7992 : DO idir = 1, 3
488 7992 : IF (irow == iatom) THEN
489 40176 : fij(idir) = SUM(p_block*dsint(:, :, idir))
490 : ELSE
491 31662 : fij(idir) = SUM(TRANSPOSE(p_block)*dsint(:, :, idir))
492 : END IF
493 : END DO
494 5454 : IF (irow == iatom) fij = -fij
495 11322 : 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 4 : CALL release_sap_int(sap_int)
503 : END IF
504 : ELSE
505 1858 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
506 494311 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
507 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
508 492453 : iatom=iatom, jatom=jatom, r=rab, cell=cellind)
509 :
510 492453 : icol = MAX(iatom, jatom)
511 492453 : irow = MIN(iatom, jatom)
512 :
513 492453 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
514 492453 : CPASSERT(ic > 0)
515 :
516 492453 : fdir = 0.0_dp
517 1969812 : ria = particle_set(irow)%r
518 1969812 : rib = particle_set(icol)%r
519 1969812 : DO idir = 1, 3
520 5909436 : kvec(:) = twopi*cell%h_inv(idir, :)
521 5909436 : dd = SUM(kvec(:)*ria(:))
522 1477359 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
523 1477359 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
524 5909436 : dd = SUM(kvec(:)*rib(:))
525 1477359 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
526 1969812 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
527 : END DO
528 :
529 492453 : NULLIFY (s_block)
530 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
531 492453 : row=irow, col=icol, block=s_block, found=found)
532 492453 : CPASSERT(found)
533 1169156 : DO is = 1, nspin
534 676703 : NULLIFY (ks_block)
535 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
536 676703 : row=irow, col=icol, block=ks_block, found=found)
537 676703 : CPASSERT(found)
538 15851147 : ks_block = ks_block + 0.5_dp*fdir*s_block
539 : END DO
540 494311 : IF (calculate_forces) THEN
541 3323 : atom_a = atom_of_kind(iatom)
542 3323 : atom_b = atom_of_kind(jatom)
543 3323 : fij = 0.0_dp
544 7986 : DO ispin = 1, nspin
545 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
546 4663 : row=irow, col=icol, BLOCK=p_block, found=found)
547 4663 : CPASSERT(found)
548 26638 : DO idir = 1, 3
549 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
550 13989 : row=irow, col=icol, BLOCK=ds_block, found=found)
551 13989 : CPASSERT(found)
552 176149 : fij(idir) = fij(idir) + SUM(p_block*ds_block)
553 : END DO
554 : END DO
555 9197 : IF (irow == iatom) fij = -fij
556 13292 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3)
557 13292 : force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) + fdir*fij(1:3)
558 3323 : IF (use_virial) THEN
559 5360 : dr = SQRT(SUM(rab(:)**2))
560 1340 : IF (dr > 1.e-6_dp) THEN
561 1332 : fi = 1.0_dp
562 1332 : IF (iatom == jatom) fi = 0.5_dp
563 5328 : 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 1858 : CALL neighbor_list_iterator_release(nl_iterator)
569 : END IF
570 :
571 3960 : IF (calculate_forces) THEN
572 154 : DO ikind = 1, SIZE(atomic_kind_set)
573 1290 : CALL para_env%sum(force(ikind)%efield)
574 : END DO
575 : END IF
576 :
577 : END IF
578 :
579 3960 : CALL timestop(handle)
580 :
581 7920 : 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 332 : 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 332 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
610 : INTEGER, DIMENSION(3) :: cellind
611 332 : 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 332 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block
618 332 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
619 : TYPE(cell_type), POINTER :: cell
620 : TYPE(dbcsr_iterator_type) :: iter
621 332 : 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 332 : DIMENSION(:), POINTER :: nl_iterator
628 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
629 332 : POINTER :: sab_orb
630 332 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
631 332 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
632 332 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
633 : TYPE(virial_type), POINTER :: virial
634 :
635 332 : CALL timeset(routineN, handle)
636 :
637 332 : NULLIFY (dft_control, cell, particle_set)
638 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
639 332 : particle_set=particle_set, virial=virial)
640 332 : NULLIFY (qs_kind_set, para_env, sab_orb)
641 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
642 332 : efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
643 :
644 : ! efield history
645 332 : CALL init_efield_matrices(efield)
646 332 : CALL set_qs_env(qs_env, efield=efield)
647 :
648 : ! calculate stress only if forces requested also
649 332 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
650 0 : use_virial = use_virial .AND. calculate_forces
651 : ! disable stress calculation
652 : IF (use_virial) THEN
653 0 : CPABORT("Stress tensor for periodic D-field not implemented")
654 : END IF
655 :
656 1328 : dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
657 :
658 1328 : fieldpol = dft_control%period_efield%polarisation
659 2324 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
660 1328 : fieldpol = fieldpol*dft_control%period_efield%strength
661 :
662 332 : omega = cell%deth
663 4316 : hmat = cell%hmat(:, :)/twopi
664 :
665 332 : natom = SIZE(particle_set)
666 332 : nspin = SIZE(ks_matrix, 1)
667 :
668 1328 : zi(:) = CMPLX(1._dp, 0._dp, dp)
669 996 : DO ia = 1, natom
670 664 : charge = mcharge(ia)
671 2656 : ria = particle_set(ia)%r
672 2988 : DO idir = 1, 3
673 7968 : kvec(:) = twopi*cell%h_inv(idir, :)
674 7968 : dd = SUM(kvec(:)*ria(:))
675 1992 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
676 2656 : zi(idir) = zi(idir)*zdeta
677 : END DO
678 : END DO
679 1328 : qi = AIMAG(LOG(zi))
680 :
681 : ! make sure the total normalized polarization is within [-1:1]
682 1328 : DO idir = 1, 3
683 996 : cqi(idir) = qi(idir)/omega
684 996 : IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi
685 996 : IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi
686 : ! now check for log branch
687 996 : IF (calculate_forces) THEN
688 30 : IF (ABS(efield%polarisation(idir) - cqi(idir)) > pi) THEN
689 0 : di(idir) = (efield%polarisation(idir) - cqi(idir))/pi
690 0 : DO i = 1, 10
691 0 : cqi(idir) = cqi(idir) + SIGN(1.0_dp, di(idir))*twopi
692 0 : IF (ABS(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
693 : END DO
694 : END IF
695 : END IF
696 1328 : cqi(idir) = cqi(idir)*omega
697 : END DO
698 1328 : DO idir = 1, 3
699 996 : ci(idir) = 0.0_dp
700 4316 : DO i = 1, 3
701 3984 : ci(idir) = ci(idir) + hmat(idir, i)*cqi(i)
702 : END DO
703 : END DO
704 : ! update the references
705 332 : IF (calculate_forces) THEN
706 40 : ener_field = SUM(ci)
707 : ! check for smoothness of energy surface
708 130 : IF (ABS(efield%field_energy - ener_field) > pi*ABS(SUM(hmat))) THEN
709 0 : CPWARN("Large change of e-field energy detected. Correct for non-smooth energy surface")
710 : END IF
711 10 : efield%field_energy = ener_field
712 40 : efield%polarisation(:) = cqi(:)/omega
713 : END IF
714 :
715 : ! Energy
716 332 : ener_field = 0.0_dp
717 1328 : DO idir = 1, 3
718 1328 : ener_field = ener_field + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi/omega*ci(idir))**2
719 : END DO
720 332 : energy%efield = 0.25_dp/twopi*ener_field
721 :
722 332 : IF (.NOT. just_energy) THEN
723 1328 : di(:) = -(fieldpol(:) - 2._dp*twopi/omega*ci(:))*dfilter(:)/omega
724 :
725 332 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
726 332 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
727 :
728 332 : nimg = dft_control%nimages
729 332 : NULLIFY (cell_to_index)
730 332 : IF (nimg > 1) THEN
731 0 : NULLIFY (kpoints)
732 0 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
733 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
734 : END IF
735 :
736 332 : IF (calculate_forces) THEN
737 10 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
738 10 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
739 10 : IF (para_env%mepos == 0) THEN
740 15 : DO ia = 1, natom
741 10 : charge = mcharge(ia)
742 10 : iatom = atom_of_kind(ia)
743 10 : ikind = kind_of(ia)
744 45 : force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + di(:)*charge
745 : END DO
746 : END IF
747 : END IF
748 :
749 332 : IF (nimg == 1) THEN
750 : ! no k-points; all matrices have been transformed to periodic bsf
751 332 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
752 830 : DO WHILE (dbcsr_iterator_blocks_left(iter))
753 498 : CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)
754 :
755 1992 : DO idir = 1, 3
756 6474 : hdi(idir) = -SUM(di(1:3)*hmat(1:3, idir))
757 : END DO
758 498 : fdir = 0.0_dp
759 1992 : ria = particle_set(irow)%r
760 1992 : rib = particle_set(icol)%r
761 1992 : DO idir = 1, 3
762 5976 : kvec(:) = twopi*cell%h_inv(idir, :)
763 5976 : dd = SUM(kvec(:)*ria(:))
764 1494 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
765 1494 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
766 5976 : dd = SUM(kvec(:)*rib(:))
767 1494 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
768 1992 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
769 : END DO
770 :
771 996 : DO is = 1, nspin
772 498 : NULLIFY (ks_block)
773 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
774 498 : row=irow, col=icol, block=ks_block, found=found)
775 498 : CPASSERT(found)
776 14110 : ks_block = ks_block + 0.5_dp*fdir*s_block
777 : END DO
778 830 : IF (calculate_forces) THEN
779 15 : ikind = kind_of(irow)
780 15 : jkind = kind_of(icol)
781 15 : atom_a = atom_of_kind(irow)
782 15 : atom_b = atom_of_kind(icol)
783 15 : fij = 0.0_dp
784 30 : DO ispin = 1, nspin
785 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
786 15 : row=irow, col=icol, BLOCK=p_block, found=found)
787 15 : CPASSERT(found)
788 90 : DO idir = 1, 3
789 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
790 45 : row=irow, col=icol, BLOCK=ds_block, found=found)
791 45 : CPASSERT(found)
792 675 : fij(idir) = fij(idir) + SUM(p_block*ds_block)
793 : END DO
794 : END DO
795 60 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
796 60 : 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 332 : CALL dbcsr_iterator_stop(iter)
801 : ELSE
802 0 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
803 0 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
804 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
805 0 : iatom=iatom, jatom=jatom, r=rab, cell=cellind)
806 :
807 0 : icol = MAX(iatom, jatom)
808 0 : irow = MIN(iatom, jatom)
809 :
810 0 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
811 0 : CPASSERT(ic > 0)
812 :
813 0 : DO idir = 1, 3
814 0 : hdi(idir) = -SUM(di(1:3)*hmat(1:3, idir))
815 : END DO
816 0 : fdir = 0.0_dp
817 0 : ria = particle_set(irow)%r
818 0 : rib = particle_set(icol)%r
819 0 : DO idir = 1, 3
820 0 : kvec(:) = twopi*cell%h_inv(idir, :)
821 0 : dd = SUM(kvec(:)*ria(:))
822 0 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
823 0 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
824 0 : dd = SUM(kvec(:)*rib(:))
825 0 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
826 0 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
827 : END DO
828 :
829 0 : NULLIFY (s_block)
830 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
831 0 : row=irow, col=icol, block=s_block, found=found)
832 0 : CPASSERT(found)
833 0 : DO is = 1, nspin
834 0 : NULLIFY (ks_block)
835 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
836 0 : row=irow, col=icol, block=ks_block, found=found)
837 0 : CPASSERT(found)
838 0 : ks_block = ks_block + 0.5_dp*fdir*s_block
839 : END DO
840 0 : IF (calculate_forces) THEN
841 0 : atom_a = atom_of_kind(iatom)
842 0 : atom_b = atom_of_kind(jatom)
843 0 : fij = 0.0_dp
844 0 : DO ispin = 1, nspin
845 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
846 0 : row=irow, col=icol, BLOCK=p_block, found=found)
847 0 : CPASSERT(found)
848 0 : DO idir = 1, 3
849 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
850 0 : row=irow, col=icol, BLOCK=ds_block, found=found)
851 0 : CPASSERT(found)
852 0 : fij(idir) = fij(idir) + SUM(p_block*ds_block)
853 : END DO
854 : END DO
855 0 : IF (irow == iatom) fij = -fij
856 0 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3)
857 0 : 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 0 : CALL neighbor_list_iterator_release(nl_iterator)
862 : END IF
863 :
864 : END IF
865 :
866 332 : CALL timestop(handle)
867 :
868 664 : END SUBROUTINE dfield_tb_berry
869 :
870 : END MODULE efield_tb_methods
|
