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