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 Calculates the energy contribution and the mo_derivative of
10 : !> a static electric field (nonperiodic)
11 : !> \par History
12 : !> none
13 : !> \author JGH (05.2015)
14 : ! **************************************************************************************************
15 : MODULE qs_efield_local
16 : USE ai_moments, ONLY: dipole_force
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind,&
19 : get_atomic_kind_set
20 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
21 : gto_basis_set_type
22 : USE cell_types, ONLY: cell_type,&
23 : pbc
24 : USE cp_control_types, ONLY: dft_control_type
25 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
26 : dbcsr_copy,&
27 : dbcsr_get_block_p,&
28 : dbcsr_p_type,&
29 : dbcsr_set
30 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
31 : USE kinds, ONLY: dp
32 : USE message_passing, ONLY: mp_para_env_type
33 : USE orbital_pointers, ONLY: ncoset
34 : USE particle_types, ONLY: particle_type
35 : USE qs_energy_types, ONLY: qs_energy_type
36 : USE qs_environment_types, ONLY: get_qs_env,&
37 : qs_environment_type,&
38 : set_qs_env
39 : USE qs_force_types, ONLY: qs_force_type
40 : USE qs_kind_types, ONLY: get_qs_kind,&
41 : qs_kind_type
42 : USE qs_moments, ONLY: build_local_moment_matrix
43 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
44 : neighbor_list_iterate,&
45 : neighbor_list_iterator_create,&
46 : neighbor_list_iterator_p_type,&
47 : neighbor_list_iterator_release,&
48 : neighbor_list_set_p_type
49 : USE qs_period_efield_types, ONLY: efield_berry_type,&
50 : init_efield_matrices,&
51 : set_efield_matrices
52 : USE qs_rho_types, ONLY: qs_rho_get,&
53 : qs_rho_type
54 : #include "./base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 :
58 : PRIVATE
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_efield_local'
61 :
62 : ! *** Public subroutines ***
63 :
64 : PUBLIC :: qs_efield_local_operator
65 :
66 : ! **************************************************************************************************
67 :
68 : CONTAINS
69 :
70 : ! **************************************************************************************************
71 :
72 : ! **************************************************************************************************
73 : !> \brief ...
74 : !> \param qs_env ...
75 : !> \param just_energy ...
76 : !> \param calculate_forces ...
77 : ! **************************************************************************************************
78 110283 : SUBROUTINE qs_efield_local_operator(qs_env, just_energy, calculate_forces)
79 :
80 : TYPE(qs_environment_type), POINTER :: qs_env
81 : LOGICAL, INTENT(IN) :: just_energy, calculate_forces
82 :
83 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_local_operator'
84 :
85 : INTEGER :: handle
86 : LOGICAL :: s_mstruct_changed
87 : REAL(dp), DIMENSION(3) :: rpoint
88 : TYPE(dft_control_type), POINTER :: dft_control
89 :
90 110283 : CALL timeset(routineN, handle)
91 :
92 110283 : NULLIFY (dft_control)
93 : CALL get_qs_env(qs_env, s_mstruct_changed=s_mstruct_changed, &
94 110283 : dft_control=dft_control)
95 :
96 110283 : IF (dft_control%apply_efield) THEN
97 9990 : rpoint = 0.0_dp
98 9990 : IF (s_mstruct_changed) CALL qs_efield_integrals(qs_env, rpoint)
99 9990 : CALL qs_efield_mo_derivatives(qs_env, rpoint, just_energy, calculate_forces)
100 : END IF
101 :
102 110283 : CALL timestop(handle)
103 :
104 110283 : END SUBROUTINE qs_efield_local_operator
105 :
106 : ! **************************************************************************************************
107 : !> \brief ...
108 : !> \param qs_env ...
109 : !> \param rpoint ...
110 : ! **************************************************************************************************
111 830 : SUBROUTINE qs_efield_integrals(qs_env, rpoint)
112 :
113 : TYPE(qs_environment_type), POINTER :: qs_env
114 : REAL(dp), DIMENSION(3), INTENT(IN) :: rpoint
115 :
116 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_integrals'
117 :
118 : INTEGER :: handle, i
119 830 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s
120 : TYPE(dft_control_type), POINTER :: dft_control
121 : TYPE(efield_berry_type), POINTER :: efield
122 :
123 830 : CALL timeset(routineN, handle)
124 830 : CPASSERT(ASSOCIATED(qs_env))
125 :
126 830 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
127 830 : NULLIFY (matrix_s)
128 830 : CALL get_qs_env(qs_env=qs_env, efield=efield, matrix_s=matrix_s)
129 830 : CALL init_efield_matrices(efield)
130 3320 : ALLOCATE (dipmat(3))
131 3320 : DO i = 1, 3
132 2490 : ALLOCATE (dipmat(i)%matrix)
133 2490 : CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'DIP MAT')
134 3320 : CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
135 : END DO
136 830 : CALL build_local_moment_matrix(qs_env, dipmat, 1, rpoint)
137 830 : CALL set_efield_matrices(efield=efield, dipmat=dipmat)
138 830 : CALL set_qs_env(qs_env=qs_env, efield=efield)
139 830 : CALL timestop(handle)
140 :
141 830 : END SUBROUTINE qs_efield_integrals
142 :
143 : ! **************************************************************************************************
144 : !> \brief ...
145 : !> \param qs_env ...
146 : !> \param rpoint ...
147 : !> \param just_energy ...
148 : !> \param calculate_forces ...
149 : ! **************************************************************************************************
150 9990 : SUBROUTINE qs_efield_mo_derivatives(qs_env, rpoint, just_energy, calculate_forces)
151 : TYPE(qs_environment_type), POINTER :: qs_env
152 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rpoint
153 : LOGICAL :: just_energy, calculate_forces
154 :
155 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_mo_derivatives'
156 :
157 : INTEGER :: atom_a, atom_b, handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, &
158 : jatom, jkind, jset, ldab, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
159 9990 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
160 9990 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
161 9990 : npgfb, nsgfa, nsgfb
162 9990 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
163 : LOGICAL :: found, trans
164 : REAL(dp) :: charge, ci(3), dab, ener_field, fdir, &
165 : fieldpol(3), tmp
166 : REAL(dp), DIMENSION(3) :: ra, rab, rac, rbc, ria
167 : REAL(dp), DIMENSION(3, 3) :: forcea, forceb
168 9990 : REAL(dp), DIMENSION(:, :), POINTER :: p_block_a, p_block_b, pblock, pmat, work
169 9990 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
170 9990 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
171 9990 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
172 : TYPE(cell_type), POINTER :: cell
173 9990 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_ks, matrix_p
174 : TYPE(dft_control_type), POINTER :: dft_control
175 : TYPE(efield_berry_type), POINTER :: efield
176 9990 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
177 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
178 : TYPE(mp_para_env_type), POINTER :: para_env
179 : TYPE(neighbor_list_iterator_p_type), &
180 9990 : DIMENSION(:), POINTER :: nl_iterator
181 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
182 9990 : POINTER :: sab_orb
183 9990 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
184 : TYPE(qs_energy_type), POINTER :: energy
185 9990 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
186 9990 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
187 : TYPE(qs_kind_type), POINTER :: qs_kind
188 : TYPE(qs_rho_type), POINTER :: rho
189 :
190 9990 : CALL timeset(routineN, handle)
191 :
192 9990 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, particle_set=particle_set)
193 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
194 9990 : efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
195 :
196 : fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
197 39960 : dft_control%efield_fields(1)%efield%strength
198 :
199 : ! nuclear contribution
200 9990 : natom = SIZE(particle_set)
201 9990 : IF (calculate_forces) THEN
202 170 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
203 170 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
204 : END IF
205 9990 : ci = 0.0_dp
206 41340 : DO ia = 1, natom
207 31350 : CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
208 31350 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
209 125400 : ria = particle_set(ia)%r - rpoint
210 125400 : ria = pbc(ria, cell)
211 125400 : ci(:) = ci(:) + charge*ria(:)
212 72690 : IF (calculate_forces) THEN
213 516 : IF (para_env%mepos == 0) THEN
214 258 : iatom = atom_of_kind(ia)
215 1032 : DO idir = 1, 3
216 1032 : force(ikind)%efield(idir, iatom) = force(ikind)%efield(idir, iatom) - fieldpol(idir)*charge
217 : END DO
218 : END IF
219 : END IF
220 : END DO
221 39960 : ener_field = -SUM(ci(:)*fieldpol(:))
222 :
223 : ! Energy
224 9990 : dipmat => efield%dipmat
225 9990 : NULLIFY (rho, matrix_p)
226 9990 : CALL get_qs_env(qs_env=qs_env, rho=rho)
227 9990 : CALL qs_rho_get(rho, rho_ao=matrix_p)
228 20894 : DO ispin = 1, SIZE(matrix_p)
229 53606 : DO idir = 1, 3
230 32712 : CALL dbcsr_dot(matrix_p(ispin)%matrix, dipmat(idir)%matrix, tmp)
231 43616 : ener_field = ener_field + fieldpol(idir)*tmp
232 : END DO
233 : END DO
234 9990 : energy%efield = ener_field
235 :
236 9990 : IF (.NOT. just_energy) THEN
237 :
238 : ! Update KS matrix
239 9990 : NULLIFY (matrix_ks)
240 9990 : CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks)
241 20894 : DO ispin = 1, SIZE(matrix_ks)
242 53606 : DO idir = 1, 3
243 : CALL dbcsr_add(matrix_ks(ispin)%matrix, dipmat(idir)%matrix, &
244 43616 : alpha_scalar=1.0_dp, beta_scalar=fieldpol(idir))
245 : END DO
246 : END DO
247 :
248 : ! forces from the efield contribution
249 9990 : IF (calculate_forces) THEN
250 170 : nkind = SIZE(qs_kind_set)
251 170 : natom = SIZE(particle_set)
252 :
253 870 : ALLOCATE (basis_set_list(nkind))
254 530 : DO ikind = 1, nkind
255 360 : qs_kind => qs_kind_set(ikind)
256 360 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
257 530 : IF (ASSOCIATED(basis_set_a)) THEN
258 360 : basis_set_list(ikind)%gto_basis_set => basis_set_a
259 : ELSE
260 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
261 : END IF
262 : END DO
263 : !
264 170 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
265 873 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
266 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
267 703 : iatom=iatom, jatom=jatom, r=rab)
268 703 : basis_set_a => basis_set_list(ikind)%gto_basis_set
269 703 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
270 703 : basis_set_b => basis_set_list(jkind)%gto_basis_set
271 703 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
272 : ! basis ikind
273 703 : first_sgfa => basis_set_a%first_sgf
274 703 : la_max => basis_set_a%lmax
275 703 : la_min => basis_set_a%lmin
276 703 : npgfa => basis_set_a%npgf
277 703 : nseta = basis_set_a%nset
278 703 : nsgfa => basis_set_a%nsgf_set
279 703 : rpgfa => basis_set_a%pgf_radius
280 703 : set_radius_a => basis_set_a%set_radius
281 703 : sphi_a => basis_set_a%sphi
282 703 : zeta => basis_set_a%zet
283 : ! basis jkind
284 703 : first_sgfb => basis_set_b%first_sgf
285 703 : lb_max => basis_set_b%lmax
286 703 : lb_min => basis_set_b%lmin
287 703 : npgfb => basis_set_b%npgf
288 703 : nsetb = basis_set_b%nset
289 703 : nsgfb => basis_set_b%nsgf_set
290 703 : rpgfb => basis_set_b%pgf_radius
291 703 : set_radius_b => basis_set_b%set_radius
292 703 : sphi_b => basis_set_b%sphi
293 703 : zetb => basis_set_b%zet
294 :
295 703 : atom_a = atom_of_kind(iatom)
296 703 : atom_b = atom_of_kind(jatom)
297 :
298 2812 : ra(:) = particle_set(iatom)%r(:) - rpoint(:)
299 703 : rac(:) = pbc(ra(:), cell)
300 2812 : rbc(:) = rac(:) + rab(:)
301 703 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
302 :
303 703 : IF (iatom <= jatom) THEN
304 454 : irow = iatom
305 454 : icol = jatom
306 454 : trans = .FALSE.
307 : ELSE
308 249 : irow = jatom
309 249 : icol = iatom
310 249 : trans = .TRUE.
311 : END IF
312 :
313 703 : fdir = 2.0_dp
314 703 : IF (iatom == jatom .AND. dab < 1.e-10_dp) fdir = 1.0_dp
315 :
316 : ! density matrix
317 703 : NULLIFY (p_block_a)
318 703 : CALL dbcsr_get_block_p(matrix_p(1)%matrix, irow, icol, p_block_a, found)
319 703 : IF (.NOT. found) CYCLE
320 703 : IF (SIZE(matrix_p) > 1) THEN
321 45 : NULLIFY (p_block_b)
322 45 : CALL dbcsr_get_block_p(matrix_p(2)%matrix, irow, icol, p_block_b, found)
323 45 : CPASSERT(found)
324 : END IF
325 703 : forcea = 0.0_dp
326 703 : forceb = 0.0_dp
327 :
328 1938 : DO iset = 1, nseta
329 1235 : ncoa = npgfa(iset)*ncoset(la_max(iset))
330 1235 : sgfa = first_sgfa(1, iset)
331 4194 : DO jset = 1, nsetb
332 2256 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
333 1870 : ncob = npgfb(jset)*ncoset(lb_max(jset))
334 1870 : sgfb = first_sgfb(1, jset)
335 : ! Calculate the primitive integrals (da|O|b) and (a|O|db)
336 1870 : ldab = MAX(ncoa, ncob)
337 13090 : ALLOCATE (work(ldab, ldab), pmat(ncoa, ncob))
338 : ! Decontract P matrix block
339 151299 : pmat = 0.0_dp
340 3824 : DO i = 1, SIZE(matrix_p)
341 1954 : IF (i == 1) THEN
342 1870 : pblock => p_block_a
343 : ELSE
344 84 : pblock => p_block_b
345 : END IF
346 1954 : IF (trans) THEN
347 : CALL dgemm("N", "T", ncoa, nsgfb(jset), nsgfa(iset), &
348 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
349 : pblock(sgfb, sgfa), SIZE(pblock, 1), &
350 647 : 0.0_dp, work(1, 1), ldab)
351 : ELSE
352 : CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
353 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
354 : pblock(sgfa, sgfb), SIZE(pblock, 1), &
355 1307 : 0.0_dp, work(1, 1), ldab)
356 : END IF
357 : CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
358 : 1.0_dp, work(1, 1), ldab, &
359 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
360 3824 : 1.0_dp, pmat(1, 1), ncoa)
361 : END DO
362 :
363 : CALL dipole_force(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
364 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
365 1870 : 1, rac, rbc, pmat, forcea, forceb)
366 :
367 3491 : DEALLOCATE (work, pmat)
368 : END DO
369 : END DO
370 :
371 2982 : DO idir = 1, 3
372 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) &
373 8436 : + fdir*fieldpol(idir)*forcea(idir, 1:3)
374 : force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) &
375 9139 : + fdir*fieldpol(idir)*forceb(idir, 1:3)
376 : END DO
377 :
378 : END DO
379 170 : CALL neighbor_list_iterator_release(nl_iterator)
380 170 : DEALLOCATE (basis_set_list)
381 : END IF
382 :
383 : END IF
384 :
385 9990 : IF (calculate_forces) THEN
386 530 : DO ikind = 1, SIZE(atomic_kind_set)
387 4658 : CALL para_env%sum(force(ikind)%efield)
388 : END DO
389 : END IF
390 :
391 9990 : CALL timestop(handle)
392 :
393 19980 : END SUBROUTINE qs_efield_mo_derivatives
394 :
395 : END MODULE qs_efield_local
|
