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