Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 cp_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 36966 : 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 36966 : CALL timeset(routineN, handle)
90 :
91 36966 : energy%efield = 0.0_dp
92 36966 : CALL get_qs_env(qs_env, dft_control=dft_control)
93 36966 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
94 36966 : IF (dft_control%apply_period_efield) THEN
95 : ! check if the periodic efield should be applied in the current step
96 4736 : 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 4736 : IF (dft_control%period_efield%displacement_field) THEN
100 332 : CALL dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
101 : ELSE
102 4404 : CALL efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
103 : END IF
104 : END IF
105 32230 : ELSE IF (dft_control%apply_efield) THEN
106 3042 : CALL efield_tb_local(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
107 29188 : ELSE IF (dft_control%apply_efield_field) THEN
108 0 : CPABORT("efield_filed")
109 : END IF
110 : ELSE
111 0 : CPABORT("This routine should only be called from TB")
112 : END IF
113 :
114 36966 : CALL timestop(handle)
115 :
116 36966 : 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 3042 : 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 3042 : 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 3042 : REAL(dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block
145 3042 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
146 : TYPE(cell_type), POINTER :: cell
147 : TYPE(dbcsr_iterator_type) :: iter
148 3042 : 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 3042 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
152 3042 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
153 3042 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
154 : TYPE(virial_type), POINTER :: virial
155 :
156 3042 : CALL timeset(routineN, handle)
157 :
158 3042 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, particle_set=particle_set)
159 3042 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, energy=energy, para_env=para_env)
160 3042 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints, virial=virial)
161 3042 : IF (do_kpoints) THEN
162 0 : CPABORT("Local electric field with kpoints not possible. Use Berry phase periodic version")
163 : END IF
164 : ! disable stress calculation
165 3042 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
166 : IF (use_virial) THEN
167 0 : CPABORT("Stress tensor for non-periodic E-field not possible")
168 : END IF
169 :
170 : fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
171 12168 : dft_control%efield_fields(1)%efield%strength
172 :
173 3042 : natom = SIZE(particle_set)
174 3042 : ci = 0.0_dp
175 14942 : DO ia = 1, natom
176 11900 : charge = mcharge(ia)
177 47600 : ria = particle_set(ia)%r
178 47600 : ria = pbc(ria, cell)
179 50642 : ci(:) = ci(:) + charge*ria(:)
180 : END DO
181 12168 : energy%efield = -SUM(ci(:)*fieldpol(:))
182 :
183 3042 : IF (.NOT. just_energy) THEN
184 :
185 3042 : IF (calculate_forces) THEN
186 34 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
187 34 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
188 34 : IF (para_env%mepos == 0) THEN
189 79 : DO ia = 1, natom
190 62 : charge = mcharge(ia)
191 62 : ikind = kind_of(ia)
192 62 : atom_a = atom_of_kind(ia)
193 265 : force(ikind)%efield(1:3, atom_a) = -charge*fieldpol(:)
194 : END DO
195 : ELSE
196 79 : DO ia = 1, natom
197 62 : ikind = kind_of(ia)
198 62 : atom_a = atom_of_kind(ia)
199 265 : force(ikind)%efield(1:3, atom_a) = 0.0_dp
200 : END DO
201 : END IF
202 34 : CALL qs_rho_get(rho, rho_ao=matrix_p)
203 : END IF
204 :
205 : ! Update KS matrix
206 3042 : nspin = SIZE(ks_matrix, 1)
207 3042 : NULLIFY (matrix_s)
208 3042 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
209 3042 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
210 17971 : DO WHILE (dbcsr_iterator_blocks_left(iter))
211 14929 : NULLIFY (ks_block, s_block, p_block)
212 14929 : CALL dbcsr_iterator_next_block(iter, irow, icol, s_block)
213 59716 : ria = particle_set(irow)%r
214 59716 : ria = pbc(ria, cell)
215 59716 : rib = particle_set(icol)%r
216 59716 : rib = pbc(rib, cell)
217 59716 : fdir = 0.5_dp*SUM(fieldpol(1:3)*(ria(1:3) + rib(1:3)))
218 31908 : DO ispin = 1, nspin
219 : CALL dbcsr_get_block_p(matrix=ks_matrix(ispin, 1)%matrix, &
220 16979 : row=irow, col=icol, BLOCK=ks_block, found=found)
221 450682 : ks_block = ks_block + fdir*s_block
222 48887 : CPASSERT(found)
223 : END DO
224 17971 : IF (calculate_forces) THEN
225 153 : ikind = kind_of(irow)
226 153 : jkind = kind_of(icol)
227 153 : atom_a = atom_of_kind(irow)
228 153 : atom_b = atom_of_kind(icol)
229 153 : fij = 0.0_dp
230 306 : DO ispin = 1, nspin
231 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
232 153 : row=irow, col=icol, BLOCK=p_block, found=found)
233 153 : CPASSERT(found)
234 918 : DO idir = 1, 3
235 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1)%matrix, &
236 459 : row=irow, col=icol, BLOCK=ds_block, found=found)
237 459 : CPASSERT(found)
238 7215 : fij(idir) = fij(idir) + SUM(p_block*ds_block)
239 : END DO
240 : END DO
241 612 : fdir = SUM(ria(1:3)*fieldpol(1:3))
242 612 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
243 612 : force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
244 612 : fdir = SUM(rib(1:3)*fieldpol(1:3))
245 612 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
246 612 : 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 3042 : CALL dbcsr_iterator_stop(iter)
250 :
251 3042 : IF (calculate_forces) THEN
252 130 : DO ikind = 1, SIZE(atomic_kind_set)
253 1122 : CALL para_env%sum(force(ikind)%efield)
254 : END DO
255 : END IF
256 :
257 : END IF
258 :
259 3042 : CALL timestop(handle)
260 :
261 6084 : 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 4404 : 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
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 4404 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
290 : INTEGER, DIMENSION(3) :: cellind
291 4404 : 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 4404 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block
298 4404 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
299 4404 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
300 : TYPE(cell_type), POINTER :: cell
301 : TYPE(dbcsr_iterator_type) :: iter
302 4404 : 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
306 : TYPE(neighbor_list_iterator_p_type), &
307 4404 : DIMENSION(:), POINTER :: nl_iterator
308 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
309 4404 : POINTER :: sab_orb
310 4404 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
311 4404 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
312 4404 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
313 4404 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
314 : TYPE(virial_type), POINTER :: virial
315 :
316 4404 : CALL timeset(routineN, handle)
317 :
318 4404 : NULLIFY (dft_control, cell, particle_set)
319 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
320 4404 : particle_set=particle_set, virial=virial)
321 4404 : NULLIFY (qs_kind_set, para_env, sab_orb)
322 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
323 4404 : energy=energy, para_env=para_env, sab_orb=sab_orb)
324 :
325 : ! calculate stress only if forces requested also
326 4544 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
327 150 : use_virial = use_virial .AND. calculate_forces
328 :
329 : ! if an intensities list is given, select the value for the current step
330 4404 : strength = dft_control%period_efield%strength
331 4404 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
332 : strength = dft_control%period_efield%strength_list(MOD(qs_env%sim_step &
333 0 : - dft_control%period_efield%start_frame, SIZE(dft_control%period_efield%strength_list)) + 1)
334 : END IF
335 :
336 17616 : fieldpol = dft_control%period_efield%polarisation
337 30828 : fieldpol = fieldpol/NORM2(fieldpol)
338 17616 : fieldpol = -fieldpol*strength
339 57252 : hmat = cell%hmat(:, :)/twopi
340 17616 : DO idir = 1, 3
341 17616 : fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir)
342 : END DO
343 :
344 4404 : natom = SIZE(particle_set)
345 4404 : nspin = SIZE(ks_matrix, 1)
346 :
347 17616 : zi(:) = CMPLX(1._dp, 0._dp, dp)
348 20952 : DO ia = 1, natom
349 16548 : charge = mcharge(ia)
350 66192 : ria = particle_set(ia)%r
351 70596 : DO idir = 1, 3
352 198576 : kvec(:) = twopi*cell%h_inv(idir, :)
353 198576 : dd = SUM(kvec(:)*ria(:))
354 49644 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
355 66192 : zi(idir) = zi(idir)*zdeta
356 : END DO
357 : END DO
358 17616 : qi = AIMAG(LOG(zi))
359 17616 : energy%efield = SUM(fpolvec(:)*qi(:))
360 :
361 4404 : IF (.NOT. just_energy) THEN
362 4404 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
363 4404 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
364 :
365 4404 : nimg = dft_control%nimages
366 4404 : NULLIFY (cell_to_index)
367 4404 : IF (nimg > 1) THEN
368 1646 : NULLIFY (kpoints)
369 1646 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
370 1646 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
371 : END IF
372 :
373 4404 : IF (calculate_forces) THEN
374 52 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
375 52 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
376 52 : IF (para_env%mepos == 0) THEN
377 109 : DO ia = 1, natom
378 83 : charge = mcharge(ia)
379 83 : iatom = atom_of_kind(ia)
380 83 : ikind = kind_of(ia)
381 332 : force(ikind)%efield(:, iatom) = fieldpol(:)*charge
382 109 : IF (use_virial) THEN
383 80 : ria = particle_set(ia)%r
384 80 : ria = pbc(ria, cell)
385 80 : forcea(1:3) = fieldpol(1:3)*charge
386 20 : CALL virial_pair_force(virial%pv_virial, -0.5_dp, forcea, ria)
387 20 : CALL virial_pair_force(virial%pv_virial, -0.5_dp, ria, forcea)
388 : END IF
389 : END DO
390 : ELSE
391 109 : DO ia = 1, natom
392 83 : iatom = atom_of_kind(ia)
393 83 : ikind = kind_of(ia)
394 358 : force(ikind)%efield(:, iatom) = 0.0_dp
395 : END DO
396 : END IF
397 : END IF
398 :
399 4404 : IF (nimg == 1) THEN
400 : ! no k-points; all matrices have been transformed to periodic bsf
401 2758 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
402 15244 : DO WHILE (dbcsr_iterator_blocks_left(iter))
403 12486 : CALL dbcsr_iterator_next_block(iter, irow, icol, s_block)
404 :
405 12486 : fdir = 0.0_dp
406 49944 : ria = particle_set(irow)%r
407 49944 : rib = particle_set(icol)%r
408 49944 : DO idir = 1, 3
409 149832 : kvec(:) = twopi*cell%h_inv(idir, :)
410 149832 : dd = SUM(kvec(:)*ria(:))
411 37458 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
412 37458 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
413 149832 : dd = SUM(kvec(:)*rib(:))
414 37458 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
415 49944 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
416 : END DO
417 :
418 29254 : DO is = 1, nspin
419 16768 : NULLIFY (ks_block)
420 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
421 16768 : row=irow, col=icol, block=ks_block, found=found)
422 16768 : CPASSERT(found)
423 440922 : ks_block = ks_block - 0.5_dp*fdir*s_block
424 : END DO
425 15244 : IF (calculate_forces .AND. irow /= icol) THEN
426 60 : ikind = kind_of(irow)
427 60 : jkind = kind_of(icol)
428 60 : atom_a = atom_of_kind(irow)
429 60 : atom_b = atom_of_kind(icol)
430 60 : fij = 0.0_dp
431 144 : DO ispin = 1, nspin
432 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
433 84 : row=irow, col=icol, BLOCK=p_block, found=found)
434 84 : CPASSERT(found)
435 480 : DO idir = 1, 3
436 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
437 252 : row=irow, col=icol, BLOCK=ds_block, found=found)
438 252 : CPASSERT(found)
439 3204 : fij(idir) = fij(idir) - SUM(p_block*ds_block)
440 : END DO
441 : END DO
442 240 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
443 240 : 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 2758 : CALL dbcsr_iterator_stop(iter)
447 : !
448 : ! stress tensor for Gamma point needs to recalculate overlap integral derivatives
449 : !
450 2758 : IF (use_virial) THEN
451 : ! derivative overlap integral (non collapsed)
452 6 : NULLIFY (sap_int)
453 6 : IF (dft_control%qs_control%dftb) THEN
454 0 : CPABORT("DFTB stress tensor for periodic efield not implemented")
455 6 : ELSE IF (dft_control%qs_control%xtb) THEN
456 6 : CALL xtb_dsint_list(qs_env, sap_int)
457 : ELSE
458 0 : CPABORT("TB method unknown")
459 : END IF
460 : !
461 6 : CALL get_qs_env(qs_env, nkind=nkind)
462 24 : DO ikind = 1, nkind
463 78 : DO jkind = 1, nkind
464 54 : iac = ikind + nkind*(jkind - 1)
465 54 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
466 72 : DO ia = 1, sap_int(iac)%nalist
467 27 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE
468 27 : iatom = sap_int(iac)%alist(ia)%aatom
469 1463 : DO ic = 1, sap_int(iac)%alist(ia)%nclist
470 1382 : jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
471 5528 : rij = sap_int(iac)%alist(ia)%clist(ic)%rac
472 5528 : dr = SQRT(SUM(rij(:)**2))
473 1409 : IF (dr > 1.e-6_dp) THEN
474 1370 : dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
475 1370 : icol = MAX(iatom, jatom)
476 1370 : irow = MIN(iatom, jatom)
477 3737 : IF (irow == iatom) rij = -rij
478 1370 : fdir = 0.0_dp
479 5480 : ria = particle_set(irow)%r
480 5480 : rib = particle_set(icol)%r
481 5480 : DO idir = 1, 3
482 16440 : kvec(:) = twopi*cell%h_inv(idir, :)
483 16440 : dd = SUM(kvec(:)*ria(:))
484 4110 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
485 4110 : fdir = fdir - fpolvec(idir)*AIMAG(LOG(zdeta))
486 16440 : dd = SUM(kvec(:)*rib(:))
487 4110 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
488 5480 : fdir = fdir - fpolvec(idir)*AIMAG(LOG(zdeta))
489 : END DO
490 1370 : fi = 1.0_dp
491 1370 : IF (iatom == jatom) fi = 0.5_dp
492 3406 : DO ispin = 1, nspin
493 2036 : NULLIFY (p_block)
494 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
495 2036 : row=irow, col=icol, block=p_block, found=found)
496 2036 : CPASSERT(found)
497 2036 : fij = 0.0_dp
498 8144 : DO idir = 1, 3
499 8144 : IF (irow == iatom) THEN
500 40725 : fij(idir) = SUM(p_block*dsint(:, :, idir))
501 : ELSE
502 32265 : fij(idir) = SUM(TRANSPOSE(p_block)*dsint(:, :, idir))
503 : END IF
504 : END DO
505 5555 : IF (irow == iatom) fij = -fij
506 11550 : 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 6 : CALL release_sap_int(sap_int)
514 : END IF
515 : ELSE
516 1646 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
517 478599 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
518 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
519 476953 : iatom=iatom, jatom=jatom, r=rab, cell=cellind)
520 :
521 476953 : icol = MAX(iatom, jatom)
522 476953 : irow = MIN(iatom, jatom)
523 :
524 476953 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
525 476953 : CPASSERT(ic > 0)
526 :
527 476953 : fdir = 0.0_dp
528 1907812 : ria = particle_set(irow)%r
529 1907812 : rib = particle_set(icol)%r
530 1907812 : DO idir = 1, 3
531 5723436 : kvec(:) = twopi*cell%h_inv(idir, :)
532 5723436 : dd = SUM(kvec(:)*ria(:))
533 1430859 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
534 1430859 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
535 5723436 : dd = SUM(kvec(:)*rib(:))
536 1430859 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
537 1907812 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
538 : END DO
539 :
540 476953 : NULLIFY (s_block)
541 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
542 476953 : row=irow, col=icol, block=s_block, found=found)
543 476953 : CPASSERT(found)
544 1148206 : DO is = 1, nspin
545 671253 : NULLIFY (ks_block)
546 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
547 671253 : row=irow, col=icol, block=ks_block, found=found)
548 671253 : CPASSERT(found)
549 15852175 : ks_block = ks_block - 0.5_dp*fdir*s_block
550 : END DO
551 478599 : IF (calculate_forces) THEN
552 3323 : atom_a = atom_of_kind(iatom)
553 3323 : atom_b = atom_of_kind(jatom)
554 3323 : fij = 0.0_dp
555 7986 : DO ispin = 1, nspin
556 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
557 4663 : row=irow, col=icol, BLOCK=p_block, found=found)
558 4663 : CPASSERT(found)
559 26638 : DO idir = 1, 3
560 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
561 13989 : row=irow, col=icol, BLOCK=ds_block, found=found)
562 13989 : CPASSERT(found)
563 176149 : fij(idir) = fij(idir) - SUM(p_block*ds_block)
564 : END DO
565 : END DO
566 9197 : IF (irow == iatom) fij = -fij
567 13292 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3)
568 13292 : force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) + fdir*fij(1:3)
569 3323 : IF (use_virial) THEN
570 5360 : dr = SQRT(SUM(rab(:)**2))
571 1340 : IF (dr > 1.e-6_dp) THEN
572 1332 : fi = 1.0_dp
573 1332 : IF (iatom == jatom) fi = 0.5_dp
574 5328 : 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 1646 : CALL neighbor_list_iterator_release(nl_iterator)
580 : END IF
581 :
582 4404 : IF (calculate_forces) THEN
583 178 : DO ikind = 1, SIZE(atomic_kind_set)
584 1506 : CALL para_env%sum(force(ikind)%efield)
585 : END DO
586 : END IF
587 :
588 : END IF
589 :
590 4404 : CALL timestop(handle)
591 :
592 8808 : 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 332 : 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', smoothness_warning = &
613 : "Large change of e-field energy detected. Correct for non-smooth energy surface"
614 :
615 : COMPLEX(KIND=dp) :: zdeta
616 : COMPLEX(KIND=dp), DIMENSION(3) :: zi
617 : INTEGER :: atom_a, atom_b, handle, i, ia, iatom, &
618 : ic, icol, idir, ikind, irow, is, &
619 : ispin, jatom, jkind, natom, nimg, nspin
620 332 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
621 : INTEGER, DIMENSION(3) :: cellind
622 332 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
623 : LOGICAL :: found, use_virial
624 : REAL(KIND=dp) :: charge, dd, ener_field, fdir, omega, &
625 : strength
626 : REAL(KIND=dp), DIMENSION(3) :: ci, cqi, dfilter, di, fieldpol, fij, &
627 : hdi, kvec, qi, rab, ria, rib
628 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
629 332 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block
630 332 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
631 : TYPE(cell_type), POINTER :: cell
632 : TYPE(dbcsr_iterator_type) :: iter
633 332 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
634 : TYPE(dft_control_type), POINTER :: dft_control
635 : TYPE(efield_berry_type), POINTER :: efield
636 : TYPE(kpoint_type), POINTER :: kpoints
637 : TYPE(mp_para_env_type), POINTER :: para_env
638 : TYPE(neighbor_list_iterator_p_type), &
639 332 : DIMENSION(:), POINTER :: nl_iterator
640 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
641 332 : POINTER :: sab_orb
642 332 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
643 332 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
644 332 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
645 : TYPE(virial_type), POINTER :: virial
646 :
647 332 : CALL timeset(routineN, handle)
648 :
649 332 : NULLIFY (dft_control, cell, particle_set)
650 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
651 332 : particle_set=particle_set, virial=virial)
652 332 : NULLIFY (qs_kind_set, para_env, sab_orb)
653 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
654 332 : efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
655 :
656 : ! efield history
657 332 : CALL init_efield_matrices(efield)
658 332 : CALL set_qs_env(qs_env, efield=efield)
659 :
660 : ! calculate stress only if forces requested also
661 332 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
662 0 : use_virial = use_virial .AND. calculate_forces
663 : ! disable stress calculation
664 : IF (use_virial) THEN
665 0 : CPABORT("Stress tensor for periodic D-field not implemented")
666 : END IF
667 :
668 1328 : dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
669 :
670 : ! if an intensities list is given, select the value for the current step
671 332 : strength = dft_control%period_efield%strength
672 332 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
673 : strength = dft_control%period_efield%strength_list(MOD(qs_env%sim_step &
674 0 : - dft_control%period_efield%start_frame, SIZE(dft_control%period_efield%strength_list)) + 1)
675 : END IF
676 :
677 1328 : fieldpol = dft_control%period_efield%polarisation
678 2324 : fieldpol = fieldpol/NORM2(fieldpol)
679 1328 : fieldpol = fieldpol*strength
680 :
681 332 : omega = cell%deth
682 4316 : hmat = cell%hmat(:, :)/twopi
683 :
684 332 : natom = SIZE(particle_set)
685 332 : nspin = SIZE(ks_matrix, 1)
686 :
687 1328 : zi(:) = CMPLX(1._dp, 0._dp, dp)
688 996 : DO ia = 1, natom
689 664 : charge = mcharge(ia)
690 2656 : ria = particle_set(ia)%r
691 2988 : DO idir = 1, 3
692 7968 : kvec(:) = twopi*cell%h_inv(idir, :)
693 7968 : dd = SUM(kvec(:)*ria(:))
694 1992 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
695 2656 : zi(idir) = zi(idir)*zdeta
696 : END DO
697 : END DO
698 1328 : qi = AIMAG(LOG(zi))
699 :
700 : ! make sure the total normalized polarization is within [-1:1]
701 1328 : DO idir = 1, 3
702 996 : cqi(idir) = qi(idir)/omega
703 996 : IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi
704 996 : IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi
705 : ! now check for log branch
706 996 : IF (calculate_forces) THEN
707 30 : IF (ABS(efield%polarisation(idir) - cqi(idir)) > pi) THEN
708 0 : di(idir) = (efield%polarisation(idir) - cqi(idir))/pi
709 0 : DO i = 1, 10
710 0 : cqi(idir) = cqi(idir) + SIGN(1.0_dp, di(idir))*twopi
711 0 : IF (ABS(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
712 : END DO
713 : END IF
714 : END IF
715 1328 : cqi(idir) = cqi(idir)*omega
716 : END DO
717 1328 : DO idir = 1, 3
718 996 : ci(idir) = 0.0_dp
719 4316 : DO i = 1, 3
720 3984 : ci(idir) = ci(idir) + hmat(idir, i)*cqi(i)
721 : END DO
722 : END DO
723 : ! update the references
724 332 : IF (calculate_forces) THEN
725 40 : ener_field = SUM(ci)
726 : ! check for smoothness of energy surface
727 130 : IF (ABS(efield%field_energy - ener_field) > pi*ABS(SUM(hmat))) THEN
728 0 : CPWARN(smoothness_warning)
729 : END IF
730 10 : efield%field_energy = ener_field
731 40 : efield%polarisation(:) = cqi(:)/omega
732 : END IF
733 :
734 : ! Energy
735 332 : ener_field = 0.0_dp
736 1328 : DO idir = 1, 3
737 1328 : ener_field = ener_field + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi/omega*ci(idir))**2
738 : END DO
739 332 : energy%efield = 0.25_dp/twopi*ener_field
740 :
741 332 : IF (.NOT. just_energy) THEN
742 1328 : di(:) = -(fieldpol(:) - 2._dp*twopi/omega*ci(:))*dfilter(:)/omega
743 :
744 332 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
745 332 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
746 :
747 332 : nimg = dft_control%nimages
748 332 : NULLIFY (cell_to_index)
749 332 : IF (nimg > 1) THEN
750 0 : NULLIFY (kpoints)
751 0 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
752 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
753 : END IF
754 :
755 332 : IF (calculate_forces) THEN
756 10 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
757 10 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
758 10 : IF (para_env%mepos == 0) THEN
759 15 : DO ia = 1, natom
760 10 : charge = mcharge(ia)
761 10 : iatom = atom_of_kind(ia)
762 10 : ikind = kind_of(ia)
763 45 : force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + di(:)*charge
764 : END DO
765 : END IF
766 : END IF
767 :
768 332 : IF (nimg == 1) THEN
769 : ! no k-points; all matrices have been transformed to periodic bsf
770 332 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
771 830 : DO WHILE (dbcsr_iterator_blocks_left(iter))
772 498 : CALL dbcsr_iterator_next_block(iter, irow, icol, s_block)
773 :
774 1992 : DO idir = 1, 3
775 6474 : hdi(idir) = -SUM(di(1:3)*hmat(1:3, idir))
776 : END DO
777 498 : fdir = 0.0_dp
778 1992 : ria = particle_set(irow)%r
779 1992 : rib = particle_set(icol)%r
780 1992 : DO idir = 1, 3
781 5976 : kvec(:) = twopi*cell%h_inv(idir, :)
782 5976 : dd = SUM(kvec(:)*ria(:))
783 1494 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
784 1494 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
785 5976 : dd = SUM(kvec(:)*rib(:))
786 1494 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
787 1992 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
788 : END DO
789 :
790 996 : DO is = 1, nspin
791 498 : NULLIFY (ks_block)
792 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
793 498 : row=irow, col=icol, block=ks_block, found=found)
794 498 : CPASSERT(found)
795 14110 : ks_block = ks_block + 0.5_dp*fdir*s_block
796 : END DO
797 830 : IF (calculate_forces) THEN
798 15 : ikind = kind_of(irow)
799 15 : jkind = kind_of(icol)
800 15 : atom_a = atom_of_kind(irow)
801 15 : atom_b = atom_of_kind(icol)
802 15 : fij = 0.0_dp
803 30 : DO ispin = 1, nspin
804 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
805 15 : row=irow, col=icol, BLOCK=p_block, found=found)
806 15 : CPASSERT(found)
807 90 : DO idir = 1, 3
808 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
809 45 : row=irow, col=icol, BLOCK=ds_block, found=found)
810 45 : CPASSERT(found)
811 675 : fij(idir) = fij(idir) + SUM(p_block*ds_block)
812 : END DO
813 : END DO
814 60 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
815 60 : force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
816 : END IF
817 :
818 : END DO
819 332 : CALL dbcsr_iterator_stop(iter)
820 : ELSE
821 0 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
822 0 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
823 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
824 0 : iatom=iatom, jatom=jatom, r=rab, cell=cellind)
825 :
826 0 : icol = MAX(iatom, jatom)
827 0 : irow = MIN(iatom, jatom)
828 :
829 0 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
830 0 : CPASSERT(ic > 0)
831 :
832 0 : DO idir = 1, 3
833 0 : hdi(idir) = -SUM(di(1:3)*hmat(1:3, idir))
834 : END DO
835 0 : fdir = 0.0_dp
836 0 : ria = particle_set(irow)%r
837 0 : rib = particle_set(icol)%r
838 0 : DO idir = 1, 3
839 0 : kvec(:) = twopi*cell%h_inv(idir, :)
840 0 : dd = SUM(kvec(:)*ria(:))
841 0 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
842 0 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
843 0 : dd = SUM(kvec(:)*rib(:))
844 0 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
845 0 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
846 : END DO
847 :
848 0 : NULLIFY (s_block)
849 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
850 0 : row=irow, col=icol, block=s_block, found=found)
851 0 : CPASSERT(found)
852 0 : DO is = 1, nspin
853 0 : NULLIFY (ks_block)
854 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
855 0 : row=irow, col=icol, block=ks_block, found=found)
856 0 : CPASSERT(found)
857 0 : ks_block = ks_block + 0.5_dp*fdir*s_block
858 : END DO
859 0 : IF (calculate_forces) THEN
860 0 : atom_a = atom_of_kind(iatom)
861 0 : atom_b = atom_of_kind(jatom)
862 0 : fij = 0.0_dp
863 0 : DO ispin = 1, nspin
864 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
865 0 : row=irow, col=icol, BLOCK=p_block, found=found)
866 0 : CPASSERT(found)
867 0 : DO idir = 1, 3
868 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
869 0 : row=irow, col=icol, BLOCK=ds_block, found=found)
870 0 : CPASSERT(found)
871 0 : fij(idir) = fij(idir) + SUM(p_block*ds_block)
872 : END DO
873 : END DO
874 0 : IF (irow == iatom) fij = -fij
875 0 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3)
876 0 : force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) + fdir*fij(1:3)
877 : END IF
878 :
879 : END DO
880 0 : CALL neighbor_list_iterator_release(nl_iterator)
881 : END IF
882 :
883 : END IF
884 :
885 332 : CALL timestop(handle)
886 :
887 664 : END SUBROUTINE dfield_tb_berry
888 :
889 : END MODULE efield_tb_methods
|
