Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of electric field contributions in TB
10 : !> \author JGH
11 : ! **************************************************************************************************
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 19200 : 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 19200 : CALL timeset(routineN, handle)
90 :
91 19200 : energy%efield = 0.0_dp
92 19200 : CALL get_qs_env(qs_env, dft_control=dft_control)
93 19200 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
94 19200 : IF (dft_control%apply_period_efield) THEN
95 : ! check if the periodic efield should be applied in the current step
96 5116 : 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 5116 : 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 4784 : CALL efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
103 : END IF
104 : END IF
105 14084 : ELSE IF (dft_control%apply_efield) THEN
106 2742 : CALL efield_tb_local(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
107 11342 : 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 19200 : CALL timestop(handle)
115 :
116 19200 : 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 2742 : 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 2742 : 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 2742 : REAL(dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block
145 2742 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
146 : TYPE(cell_type), POINTER :: cell
147 : TYPE(dbcsr_iterator_type) :: iter
148 2742 : 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 2742 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
152 2742 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
153 2742 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
154 : TYPE(virial_type), POINTER :: virial
155 :
156 2742 : CALL timeset(routineN, handle)
157 :
158 2742 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, particle_set=particle_set)
159 2742 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, energy=energy, para_env=para_env)
160 2742 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints, virial=virial)
161 2742 : 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 2742 : 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 10968 : dft_control%efield_fields(1)%efield%strength
172 :
173 2742 : natom = SIZE(particle_set)
174 2742 : ci = 0.0_dp
175 13430 : DO ia = 1, natom
176 10688 : charge = mcharge(ia)
177 42752 : ria = particle_set(ia)%r
178 42752 : ria = pbc(ria, cell)
179 45494 : ci(:) = ci(:) + charge*ria(:)
180 : END DO
181 10968 : energy%efield = -SUM(ci(:)*fieldpol(:))
182 :
183 2742 : IF (.NOT. just_energy) THEN
184 :
185 2742 : IF (calculate_forces) THEN
186 36 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
187 36 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
188 36 : IF (para_env%mepos == 0) THEN
189 84 : DO ia = 1, natom
190 66 : charge = mcharge(ia)
191 66 : ikind = kind_of(ia)
192 66 : atom_a = atom_of_kind(ia)
193 282 : force(ikind)%efield(1:3, atom_a) = -charge*fieldpol(:)
194 : END DO
195 : ELSE
196 84 : DO ia = 1, natom
197 66 : ikind = kind_of(ia)
198 66 : atom_a = atom_of_kind(ia)
199 282 : force(ikind)%efield(1:3, atom_a) = 0.0_dp
200 : END DO
201 : END IF
202 36 : CALL qs_rho_get(rho, rho_ao=matrix_p)
203 : END IF
204 :
205 : ! Update KS matrix
206 2742 : nspin = SIZE(ks_matrix, 1)
207 2742 : NULLIFY (matrix_s)
208 2742 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
209 2742 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
210 16150 : DO WHILE (dbcsr_iterator_blocks_left(iter))
211 13408 : NULLIFY (ks_block, s_block, p_block)
212 13408 : CALL dbcsr_iterator_next_block(iter, irow, icol, s_block)
213 53632 : ria = particle_set(irow)%r
214 53632 : ria = pbc(ria, cell)
215 53632 : rib = particle_set(icol)%r
216 53632 : rib = pbc(rib, cell)
217 53632 : fdir = 0.5_dp*SUM(fieldpol(1:3)*(ria(1:3) + rib(1:3)))
218 26816 : DO ispin = 1, nspin
219 : CALL dbcsr_get_block_p(matrix=ks_matrix(ispin, 1)%matrix, &
220 13408 : row=irow, col=icol, BLOCK=ks_block, found=found)
221 359280 : ks_block = ks_block + fdir*s_block
222 40224 : CPASSERT(found)
223 : END DO
224 16150 : IF (calculate_forces) THEN
225 163 : ikind = kind_of(irow)
226 163 : jkind = kind_of(icol)
227 163 : atom_a = atom_of_kind(irow)
228 163 : atom_b = atom_of_kind(icol)
229 163 : fij = 0.0_dp
230 326 : DO ispin = 1, nspin
231 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
232 163 : row=irow, col=icol, BLOCK=p_block, found=found)
233 163 : CPASSERT(found)
234 978 : DO idir = 1, 3
235 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1)%matrix, &
236 489 : row=irow, col=icol, BLOCK=ds_block, found=found)
237 489 : CPASSERT(found)
238 7639 : fij(idir) = fij(idir) + SUM(p_block*ds_block)
239 : END DO
240 : END DO
241 652 : fdir = SUM(ria(1:3)*fieldpol(1:3))
242 652 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
243 652 : force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
244 652 : fdir = SUM(rib(1:3)*fieldpol(1:3))
245 652 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
246 652 : 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 2742 : CALL dbcsr_iterator_stop(iter)
250 :
251 2742 : IF (calculate_forces) THEN
252 138 : DO ikind = 1, SIZE(atomic_kind_set)
253 1194 : CALL para_env%sum(force(ikind)%efield)
254 : END DO
255 : END IF
256 :
257 : END IF
258 :
259 2742 : CALL timestop(handle)
260 :
261 5484 : 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 4784 : SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
274 : TYPE(qs_environment_type), POINTER :: qs_env
275 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
276 : TYPE(qs_rho_type), POINTER :: rho
277 : REAL(dp), DIMENSION(:), INTENT(in) :: mcharge
278 : TYPE(qs_energy_type), POINTER :: energy
279 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
280 :
281 : CHARACTER(LEN=*), PARAMETER :: routineN = 'efield_tb_berry'
282 :
283 : COMPLEX(KIND=dp) :: zdeta
284 : COMPLEX(KIND=dp), DIMENSION(3) :: zi(3)
285 : INTEGER :: atom_a, atom_b, handle, ia, iac, iatom, &
286 : ic, icol, idir, ikind, irow, is, &
287 : ispin, jatom, jkind, natom, nimg, &
288 : nkind, nspin
289 4784 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
290 : INTEGER, DIMENSION(3) :: cellind
291 4784 : 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 4784 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block
298 4784 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
299 4784 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
300 : TYPE(cell_type), POINTER :: cell
301 : TYPE(dbcsr_iterator_type) :: iter
302 4784 : 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 4784 : DIMENSION(:), POINTER :: nl_iterator
308 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
309 4784 : POINTER :: sab_orb
310 4784 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
311 4784 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
312 4784 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
313 4784 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
314 : TYPE(virial_type), POINTER :: virial
315 :
316 4784 : CALL timeset(routineN, handle)
317 :
318 4784 : NULLIFY (dft_control, cell, particle_set)
319 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
320 4784 : particle_set=particle_set, virial=virial)
321 4784 : NULLIFY (qs_kind_set, para_env, sab_orb)
322 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
323 4784 : energy=energy, para_env=para_env, sab_orb=sab_orb)
324 :
325 : ! calculate stress only if forces requested also
326 4920 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
327 146 : use_virial = use_virial .AND. calculate_forces
328 :
329 : ! if an intensities list is given, select the value for the current step
330 4784 : strength = dft_control%period_efield%strength
331 4784 : 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 19136 : fieldpol = dft_control%period_efield%polarisation
337 33488 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
338 19136 : fieldpol = -fieldpol*strength
339 62192 : hmat = cell%hmat(:, :)/twopi
340 19136 : DO idir = 1, 3
341 19136 : fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir)
342 : END DO
343 :
344 4784 : natom = SIZE(particle_set)
345 4784 : nspin = SIZE(ks_matrix, 1)
346 :
347 19136 : zi(:) = CMPLX(1._dp, 0._dp, dp)
348 22598 : DO ia = 1, natom
349 17814 : charge = mcharge(ia)
350 71256 : ria = particle_set(ia)%r
351 76040 : DO idir = 1, 3
352 213768 : kvec(:) = twopi*cell%h_inv(idir, :)
353 213768 : dd = SUM(kvec(:)*ria(:))
354 53442 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
355 71256 : zi(idir) = zi(idir)*zdeta
356 : END DO
357 : END DO
358 19136 : qi = AIMAG(LOG(zi))
359 19136 : energy%efield = SUM(fpolvec(:)*qi(:))
360 :
361 4784 : IF (.NOT. just_energy) THEN
362 4784 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
363 4784 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
364 :
365 4784 : nimg = dft_control%nimages
366 4784 : NULLIFY (cell_to_index)
367 4784 : IF (nimg > 1) THEN
368 1870 : NULLIFY (kpoints)
369 1870 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
370 1870 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
371 : END IF
372 :
373 4784 : IF (calculate_forces) THEN
374 54 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
375 54 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
376 54 : IF (para_env%mepos == 0) THEN
377 114 : DO ia = 1, natom
378 87 : charge = mcharge(ia)
379 87 : iatom = atom_of_kind(ia)
380 87 : ikind = kind_of(ia)
381 348 : force(ikind)%efield(:, iatom) = fieldpol(:)*charge
382 114 : 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 114 : DO ia = 1, natom
392 87 : iatom = atom_of_kind(ia)
393 87 : ikind = kind_of(ia)
394 375 : force(ikind)%efield(:, iatom) = 0.0_dp
395 : END DO
396 : END IF
397 : END IF
398 :
399 4784 : IF (nimg == 1) THEN
400 : ! no k-points; all matrices have been transformed to periodic bsf
401 2914 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
402 16173 : DO WHILE (dbcsr_iterator_blocks_left(iter))
403 13259 : CALL dbcsr_iterator_next_block(iter, irow, icol, s_block)
404 :
405 13259 : fdir = 0.0_dp
406 53036 : ria = particle_set(irow)%r
407 53036 : rib = particle_set(icol)%r
408 53036 : DO idir = 1, 3
409 159108 : kvec(:) = twopi*cell%h_inv(idir, :)
410 159108 : dd = SUM(kvec(:)*ria(:))
411 39777 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
412 39777 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
413 159108 : dd = SUM(kvec(:)*rib(:))
414 39777 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
415 53036 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
416 : END DO
417 :
418 30910 : DO is = 1, nspin
419 17651 : NULLIFY (ks_block)
420 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
421 17651 : row=irow, col=icol, block=ks_block, found=found)
422 17651 : CPASSERT(found)
423 464305 : ks_block = ks_block - 0.5_dp*fdir*s_block
424 : END DO
425 16173 : IF (calculate_forces .AND. irow /= icol) THEN
426 66 : ikind = kind_of(irow)
427 66 : jkind = kind_of(icol)
428 66 : atom_a = atom_of_kind(irow)
429 66 : atom_b = atom_of_kind(icol)
430 66 : fij = 0.0_dp
431 156 : DO ispin = 1, nspin
432 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
433 90 : row=irow, col=icol, BLOCK=p_block, found=found)
434 90 : CPASSERT(found)
435 516 : DO idir = 1, 3
436 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
437 270 : row=irow, col=icol, BLOCK=ds_block, found=found)
438 270 : CPASSERT(found)
439 3444 : fij(idir) = fij(idir) - SUM(p_block*ds_block)
440 : END DO
441 : END DO
442 264 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
443 264 : 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 2914 : CALL dbcsr_iterator_stop(iter)
447 : !
448 : ! stress tensor for Gamma point needs to recalculate overlap integral derivatives
449 : !
450 2914 : 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 : ELSEIF (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 1870 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
517 498343 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
518 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
519 496473 : iatom=iatom, jatom=jatom, r=rab, cell=cellind)
520 :
521 496473 : icol = MAX(iatom, jatom)
522 496473 : irow = MIN(iatom, jatom)
523 :
524 496473 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
525 496473 : CPASSERT(ic > 0)
526 :
527 496473 : fdir = 0.0_dp
528 1985892 : ria = particle_set(irow)%r
529 1985892 : rib = particle_set(icol)%r
530 1985892 : DO idir = 1, 3
531 5957676 : kvec(:) = twopi*cell%h_inv(idir, :)
532 5957676 : dd = SUM(kvec(:)*ria(:))
533 1489419 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
534 1489419 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
535 5957676 : dd = SUM(kvec(:)*rib(:))
536 1489419 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
537 1985892 : fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
538 : END DO
539 :
540 496473 : NULLIFY (s_block)
541 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
542 496473 : row=irow, col=icol, block=s_block, found=found)
543 496473 : CPASSERT(found)
544 1180546 : DO is = 1, nspin
545 684073 : NULLIFY (ks_block)
546 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
547 684073 : row=irow, col=icol, block=ks_block, found=found)
548 684073 : CPASSERT(found)
549 16026019 : ks_block = ks_block - 0.5_dp*fdir*s_block
550 : END DO
551 498343 : 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 1870 : CALL neighbor_list_iterator_release(nl_iterator)
580 : END IF
581 :
582 4784 : IF (calculate_forces) THEN
583 186 : DO ikind = 1, SIZE(atomic_kind_set)
584 1578 : CALL para_env%sum(force(ikind)%efield)
585 : END DO
586 : END IF
587 :
588 : END IF
589 :
590 4784 : CALL timestop(handle)
591 :
592 9568 : 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'
613 :
614 : COMPLEX(KIND=dp) :: zdeta
615 : COMPLEX(KIND=dp), DIMENSION(3) :: zi(3)
616 : INTEGER :: atom_a, atom_b, handle, i, ia, iatom, &
617 : ic, icol, idir, ikind, irow, is, &
618 : ispin, jatom, jkind, natom, nimg, nspin
619 332 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
620 : INTEGER, DIMENSION(3) :: cellind
621 332 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
622 : LOGICAL :: found, use_virial
623 : REAL(KIND=dp) :: charge, dd, ener_field, fdir, omega, &
624 : strength
625 : REAL(KIND=dp), DIMENSION(3) :: ci, cqi, dfilter, di, fieldpol, fij, &
626 : hdi, kvec, qi, rab, ria, rib
627 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
628 332 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block
629 332 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
630 : TYPE(cell_type), POINTER :: cell
631 : TYPE(dbcsr_iterator_type) :: iter
632 332 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
633 : TYPE(dft_control_type), POINTER :: dft_control
634 : TYPE(efield_berry_type), POINTER :: efield
635 : TYPE(kpoint_type), POINTER :: kpoints
636 : TYPE(mp_para_env_type), POINTER :: para_env
637 : TYPE(neighbor_list_iterator_p_type), &
638 332 : DIMENSION(:), POINTER :: nl_iterator
639 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
640 332 : POINTER :: sab_orb
641 332 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
642 332 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
643 332 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
644 : TYPE(virial_type), POINTER :: virial
645 :
646 332 : CALL timeset(routineN, handle)
647 :
648 332 : NULLIFY (dft_control, cell, particle_set)
649 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
650 332 : particle_set=particle_set, virial=virial)
651 332 : NULLIFY (qs_kind_set, para_env, sab_orb)
652 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
653 332 : efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
654 :
655 : ! efield history
656 332 : CALL init_efield_matrices(efield)
657 332 : CALL set_qs_env(qs_env, efield=efield)
658 :
659 : ! calculate stress only if forces requested also
660 332 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
661 0 : use_virial = use_virial .AND. calculate_forces
662 : ! disable stress calculation
663 : IF (use_virial) THEN
664 0 : CPABORT("Stress tensor for periodic D-field not implemented")
665 : END IF
666 :
667 1328 : dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
668 :
669 : ! if an intensities list is given, select the value for the current step
670 332 : strength = dft_control%period_efield%strength
671 332 : IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
672 : strength = dft_control%period_efield%strength_list(MOD(qs_env%sim_step &
673 0 : - dft_control%period_efield%start_frame, SIZE(dft_control%period_efield%strength_list)) + 1)
674 : END IF
675 :
676 1328 : fieldpol = dft_control%period_efield%polarisation
677 2324 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
678 1328 : fieldpol = fieldpol*strength
679 :
680 332 : omega = cell%deth
681 4316 : hmat = cell%hmat(:, :)/twopi
682 :
683 332 : natom = SIZE(particle_set)
684 332 : nspin = SIZE(ks_matrix, 1)
685 :
686 1328 : zi(:) = CMPLX(1._dp, 0._dp, dp)
687 996 : DO ia = 1, natom
688 664 : charge = mcharge(ia)
689 2656 : ria = particle_set(ia)%r
690 2988 : DO idir = 1, 3
691 7968 : kvec(:) = twopi*cell%h_inv(idir, :)
692 7968 : dd = SUM(kvec(:)*ria(:))
693 1992 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
694 2656 : zi(idir) = zi(idir)*zdeta
695 : END DO
696 : END DO
697 1328 : qi = AIMAG(LOG(zi))
698 :
699 : ! make sure the total normalized polarization is within [-1:1]
700 1328 : DO idir = 1, 3
701 996 : cqi(idir) = qi(idir)/omega
702 996 : IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi
703 996 : IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi
704 : ! now check for log branch
705 996 : IF (calculate_forces) THEN
706 30 : IF (ABS(efield%polarisation(idir) - cqi(idir)) > pi) THEN
707 0 : di(idir) = (efield%polarisation(idir) - cqi(idir))/pi
708 0 : DO i = 1, 10
709 0 : cqi(idir) = cqi(idir) + SIGN(1.0_dp, di(idir))*twopi
710 0 : IF (ABS(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
711 : END DO
712 : END IF
713 : END IF
714 1328 : cqi(idir) = cqi(idir)*omega
715 : END DO
716 1328 : DO idir = 1, 3
717 996 : ci(idir) = 0.0_dp
718 4316 : DO i = 1, 3
719 3984 : ci(idir) = ci(idir) + hmat(idir, i)*cqi(i)
720 : END DO
721 : END DO
722 : ! update the references
723 332 : IF (calculate_forces) THEN
724 40 : ener_field = SUM(ci)
725 : ! check for smoothness of energy surface
726 130 : IF (ABS(efield%field_energy - ener_field) > pi*ABS(SUM(hmat))) THEN
727 0 : CPWARN("Large change of e-field energy detected. Correct for non-smooth energy surface")
728 : END IF
729 10 : efield%field_energy = ener_field
730 40 : efield%polarisation(:) = cqi(:)/omega
731 : END IF
732 :
733 : ! Energy
734 332 : ener_field = 0.0_dp
735 1328 : DO idir = 1, 3
736 1328 : ener_field = ener_field + dfilter(idir)*(fieldpol(idir) - 2._dp*twopi/omega*ci(idir))**2
737 : END DO
738 332 : energy%efield = 0.25_dp/twopi*ener_field
739 :
740 332 : IF (.NOT. just_energy) THEN
741 1328 : di(:) = -(fieldpol(:) - 2._dp*twopi/omega*ci(:))*dfilter(:)/omega
742 :
743 332 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
744 332 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
745 :
746 332 : nimg = dft_control%nimages
747 332 : NULLIFY (cell_to_index)
748 332 : IF (nimg > 1) THEN
749 0 : NULLIFY (kpoints)
750 0 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
751 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
752 : END IF
753 :
754 332 : IF (calculate_forces) THEN
755 10 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
756 10 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
757 10 : IF (para_env%mepos == 0) THEN
758 15 : DO ia = 1, natom
759 10 : charge = mcharge(ia)
760 10 : iatom = atom_of_kind(ia)
761 10 : ikind = kind_of(ia)
762 45 : force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + di(:)*charge
763 : END DO
764 : END IF
765 : END IF
766 :
767 332 : IF (nimg == 1) THEN
768 : ! no k-points; all matrices have been transformed to periodic bsf
769 332 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
770 830 : DO WHILE (dbcsr_iterator_blocks_left(iter))
771 498 : CALL dbcsr_iterator_next_block(iter, irow, icol, s_block)
772 :
773 1992 : DO idir = 1, 3
774 6474 : hdi(idir) = -SUM(di(1:3)*hmat(1:3, idir))
775 : END DO
776 498 : fdir = 0.0_dp
777 1992 : ria = particle_set(irow)%r
778 1992 : rib = particle_set(icol)%r
779 1992 : DO idir = 1, 3
780 5976 : kvec(:) = twopi*cell%h_inv(idir, :)
781 5976 : dd = SUM(kvec(:)*ria(:))
782 1494 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
783 1494 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
784 5976 : dd = SUM(kvec(:)*rib(:))
785 1494 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
786 1992 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
787 : END DO
788 :
789 996 : DO is = 1, nspin
790 498 : NULLIFY (ks_block)
791 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
792 498 : row=irow, col=icol, block=ks_block, found=found)
793 498 : CPASSERT(found)
794 14110 : ks_block = ks_block + 0.5_dp*fdir*s_block
795 : END DO
796 830 : IF (calculate_forces) THEN
797 15 : ikind = kind_of(irow)
798 15 : jkind = kind_of(icol)
799 15 : atom_a = atom_of_kind(irow)
800 15 : atom_b = atom_of_kind(icol)
801 15 : fij = 0.0_dp
802 30 : DO ispin = 1, nspin
803 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
804 15 : row=irow, col=icol, BLOCK=p_block, found=found)
805 15 : CPASSERT(found)
806 90 : DO idir = 1, 3
807 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
808 45 : row=irow, col=icol, BLOCK=ds_block, found=found)
809 45 : CPASSERT(found)
810 675 : fij(idir) = fij(idir) + SUM(p_block*ds_block)
811 : END DO
812 : END DO
813 60 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
814 60 : force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
815 : END IF
816 :
817 : END DO
818 332 : CALL dbcsr_iterator_stop(iter)
819 : ELSE
820 0 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
821 0 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
822 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
823 0 : iatom=iatom, jatom=jatom, r=rab, cell=cellind)
824 :
825 0 : icol = MAX(iatom, jatom)
826 0 : irow = MIN(iatom, jatom)
827 :
828 0 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
829 0 : CPASSERT(ic > 0)
830 :
831 0 : DO idir = 1, 3
832 0 : hdi(idir) = -SUM(di(1:3)*hmat(1:3, idir))
833 : END DO
834 0 : fdir = 0.0_dp
835 0 : ria = particle_set(irow)%r
836 0 : rib = particle_set(icol)%r
837 0 : DO idir = 1, 3
838 0 : kvec(:) = twopi*cell%h_inv(idir, :)
839 0 : dd = SUM(kvec(:)*ria(:))
840 0 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
841 0 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
842 0 : dd = SUM(kvec(:)*rib(:))
843 0 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
844 0 : fdir = fdir + hdi(idir)*AIMAG(LOG(zdeta))
845 : END DO
846 :
847 0 : NULLIFY (s_block)
848 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
849 0 : row=irow, col=icol, block=s_block, found=found)
850 0 : CPASSERT(found)
851 0 : DO is = 1, nspin
852 0 : NULLIFY (ks_block)
853 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
854 0 : row=irow, col=icol, block=ks_block, found=found)
855 0 : CPASSERT(found)
856 0 : ks_block = ks_block + 0.5_dp*fdir*s_block
857 : END DO
858 0 : IF (calculate_forces) THEN
859 0 : atom_a = atom_of_kind(iatom)
860 0 : atom_b = atom_of_kind(jatom)
861 0 : fij = 0.0_dp
862 0 : DO ispin = 1, nspin
863 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
864 0 : row=irow, col=icol, BLOCK=p_block, found=found)
865 0 : CPASSERT(found)
866 0 : DO idir = 1, 3
867 : CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
868 0 : row=irow, col=icol, BLOCK=ds_block, found=found)
869 0 : CPASSERT(found)
870 0 : fij(idir) = fij(idir) + SUM(p_block*ds_block)
871 : END DO
872 : END DO
873 0 : IF (irow == iatom) fij = -fij
874 0 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3)
875 0 : force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) + fdir*fij(1:3)
876 : END IF
877 :
878 : END DO
879 0 : CALL neighbor_list_iterator_release(nl_iterator)
880 : END IF
881 :
882 : END IF
883 :
884 332 : CALL timestop(handle)
885 :
886 664 : END SUBROUTINE dfield_tb_berry
887 :
888 : END MODULE efield_tb_methods
|
