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 Routines to perform the RTP in the velocity gauge
10 : ! **************************************************************************************************
11 :
12 : MODULE rt_propagation_velocity_gauge
13 : USE ai_moments, ONLY: cossin
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind_set
16 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
17 : gto_basis_set_type
18 : USE bibliography, ONLY: Mattiat2022,&
19 : cite_reference
20 : USE cell_types, ONLY: cell_type,&
21 : pbc
22 : USE core_ppnl, ONLY: build_core_ppnl
23 : USE cp_control_types, ONLY: dft_control_type
24 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
25 : dbcsr_create,&
26 : dbcsr_get_block_p,&
27 : dbcsr_init_p,&
28 : dbcsr_p_type,&
29 : dbcsr_set,&
30 : dbcsr_type_antisymmetric,&
31 : dbcsr_type_symmetric
32 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
33 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
34 : dbcsr_deallocate_matrix_set
35 : USE efield_utils, ONLY: make_field
36 : USE external_potential_types, ONLY: gth_potential_p_type,&
37 : gth_potential_type,&
38 : sgp_potential_p_type,&
39 : sgp_potential_type
40 : USE input_section_types, ONLY: section_vals_type
41 : USE kinds, ONLY: dp,&
42 : int_8
43 : USE kpoint_types, ONLY: get_kpoint_info,&
44 : kpoint_type
45 : USE mathconstants, ONLY: one,&
46 : zero
47 : USE orbital_pointers, ONLY: coset,&
48 : init_orbital_pointers,&
49 : nco,&
50 : ncoset
51 : USE particle_types, ONLY: particle_type
52 : USE qs_environment_types, ONLY: get_qs_env,&
53 : qs_environment_type
54 : USE qs_force_types, ONLY: qs_force_type
55 : USE qs_kind_types, ONLY: get_qs_kind,&
56 : get_qs_kind_set,&
57 : qs_kind_type
58 : USE qs_ks_types, ONLY: get_ks_env,&
59 : qs_ks_env_type
60 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
61 : USE qs_operators_ao, ONLY: build_lin_mom_matrix
62 : USE qs_rho_types, ONLY: qs_rho_get,&
63 : qs_rho_type
64 : USE sap_kind_types, ONLY: alist_type,&
65 : clist_type,&
66 : get_alist,&
67 : release_sap_int,&
68 : sap_int_type,&
69 : sap_sort
70 : USE virial_types, ONLY: virial_type
71 :
72 : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
73 : !$ omp_init_lock, omp_set_lock, &
74 : !$ omp_unset_lock, omp_destroy_lock
75 :
76 : #include "./base/base_uses.f90"
77 :
78 : IMPLICIT NONE
79 :
80 : PRIVATE
81 :
82 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_velocity_gauge'
83 :
84 : PUBLIC :: velocity_gauge_ks_matrix, update_vector_potential, velocity_gauge_nl_force, &
85 : velocity_gauge_nl_commutator
86 :
87 : CONTAINS
88 :
89 : ! **************************************************************************************************
90 : !> \brief ...
91 : !> \param qs_env ...
92 : !> \param subtract_nl_term ...
93 : ! **************************************************************************************************
94 82 : SUBROUTINE velocity_gauge_ks_matrix(qs_env, subtract_nl_term)
95 : TYPE(qs_environment_type), POINTER :: qs_env
96 : LOGICAL, INTENT(IN), OPTIONAL :: subtract_nl_term
97 :
98 : CHARACTER(len=*), PARAMETER :: routineN = 'velocity_gauge_ks_matrix'
99 :
100 : INTEGER :: handle, idir, image, nder, nimages
101 82 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
102 : LOGICAL :: calculate_forces, my_subtract_nl_term, &
103 : ppnl_present, use_virial
104 : REAL(KIND=dp) :: eps_ppnl, factor
105 : REAL(KIND=dp), DIMENSION(3) :: vec_pot
106 82 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
107 : TYPE(cell_type), POINTER :: cell
108 82 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: momentum, nl_term
109 82 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im, matrix_nl, &
110 82 : matrix_p, matrix_s
111 : TYPE(dft_control_type), POINTER :: dft_control
112 : TYPE(kpoint_type), POINTER :: kpoints
113 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
114 82 : POINTER :: sab_orb, sap_ppnl
115 82 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
116 82 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
117 82 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
118 : TYPE(qs_ks_env_type), POINTER :: ks_env
119 : TYPE(qs_rho_type), POINTER :: rho
120 : TYPE(section_vals_type), POINTER :: input
121 : TYPE(virial_type), POINTER :: virial
122 :
123 82 : CALL timeset(routineN, handle)
124 :
125 82 : CALL cite_reference(Mattiat2022)
126 :
127 82 : my_subtract_nl_term = .FALSE.
128 82 : IF (PRESENT(subtract_nl_term)) my_subtract_nl_term = subtract_nl_term
129 :
130 82 : NULLIFY (dft_control, matrix_s, sab_orb, matrix_h, cell, input, matrix_h_im, kpoints, cell_to_index, &
131 82 : sap_ppnl, particle_set, qs_kind_set, atomic_kind_set, virial, force, matrix_p, rho, matrix_nl)
132 :
133 : CALL get_qs_env(qs_env, &
134 : rho=rho, &
135 : dft_control=dft_control, &
136 : sab_orb=sab_orb, &
137 : sap_ppnl=sap_ppnl, &
138 : matrix_s_kp=matrix_s, &
139 : matrix_h_kp=matrix_h, &
140 : cell=cell, &
141 : input=input, &
142 82 : matrix_h_im_kp=matrix_h_im)
143 :
144 82 : nimages = dft_control%nimages
145 82 : ppnl_present = ASSOCIATED(sap_ppnl)
146 :
147 82 : IF (nimages > 1) THEN
148 0 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
149 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
150 : END IF
151 :
152 82 : IF (my_subtract_nl_term) THEN
153 12 : IF (ppnl_present) THEN
154 : CALL get_qs_env(qs_env, &
155 : qs_kind_set=qs_kind_set, &
156 : particle_set=particle_set, &
157 : atomic_kind_set=atomic_kind_set, &
158 : virial=virial, &
159 : rho=rho, &
160 8 : force=force)
161 :
162 8 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
163 8 : calculate_forces = .FALSE.
164 8 : use_virial = .FALSE.
165 8 : nder = 1
166 8 : eps_ppnl = dft_control%qs_control%eps_ppnl
167 :
168 8 : CALL dbcsr_allocate_matrix_set(matrix_nl, 1, nimages)
169 16 : DO image = 1, nimages
170 8 : ALLOCATE (matrix_nl(1, image)%matrix)
171 8 : CALL dbcsr_create(matrix_nl(1, image)%matrix, template=matrix_s(1, 1)%matrix)
172 8 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_nl(1, image)%matrix, sab_orb)
173 16 : CALL dbcsr_set(matrix_nl(1, image)%matrix, zero)
174 : END DO
175 :
176 : CALL build_core_ppnl(matrix_nl, matrix_p, force, virial, calculate_forces, use_virial, nder, &
177 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
178 8 : nimages, cell_to_index, "ORB")
179 :
180 16 : DO image = 1, nimages
181 16 : CALL dbcsr_add(matrix_h(1, image)%matrix, matrix_nl(1, image)%matrix, one, -one)
182 : END DO
183 :
184 8 : CALL dbcsr_deallocate_matrix_set(matrix_nl)
185 : END IF
186 : END IF
187 :
188 : !get vector potential
189 328 : vec_pot = dft_control%rtp_control%vec_pot
190 :
191 : ! allocate and build matrices for linear momentum term
192 82 : NULLIFY (momentum)
193 82 : CALL dbcsr_allocate_matrix_set(momentum, 3)
194 328 : DO idir = 1, 3
195 246 : CALL dbcsr_init_p(momentum(idir)%matrix)
196 : CALL dbcsr_create(momentum(idir)%matrix, template=matrix_s(1, 1)%matrix, &
197 246 : matrix_type=dbcsr_type_antisymmetric)
198 246 : CALL cp_dbcsr_alloc_block_from_nbl(momentum(idir)%matrix, sab_orb)
199 328 : CALL dbcsr_set(momentum(idir)%matrix, zero)
200 : END DO
201 82 : CALL build_lin_mom_matrix(qs_env, momentum)
202 :
203 : ! set imaginary part of KS matrix to zero
204 164 : DO image = 1, nimages
205 164 : CALL dbcsr_set(matrix_h_im(1, image)%matrix, zero)
206 : END DO
207 :
208 : ! add linear term in vector potential to imaginary part of KS-matrix
209 164 : DO image = 1, nimages
210 410 : DO idir = 1, 3
211 328 : CALL dbcsr_add(matrix_h_im(1, image)%matrix, momentum(idir)%matrix, one, -vec_pot(idir))
212 : END DO
213 : END DO
214 :
215 82 : CALL dbcsr_deallocate_matrix_set(momentum)
216 :
217 : ! add quadratic term to real part of KS matrix
218 82 : factor = 0._dp
219 328 : DO idir = 1, 3
220 328 : factor = factor + vec_pot(idir)**2
221 : END DO
222 :
223 164 : DO image = 1, nimages
224 164 : CALL dbcsr_add(matrix_h(1, image)%matrix, matrix_s(1, image)%matrix, one, 0.5*factor)
225 : END DO
226 :
227 : ! add Non local term
228 82 : IF (ppnl_present) THEN
229 60 : IF (dft_control%rtp_control%nl_gauge_transform) THEN
230 60 : NULLIFY (nl_term)
231 60 : CALL dbcsr_allocate_matrix_set(nl_term, 2)
232 :
233 60 : CALL dbcsr_init_p(nl_term(1)%matrix)
234 : CALL dbcsr_create(nl_term(1)%matrix, template=matrix_s(1, 1)%matrix, &
235 60 : matrix_type=dbcsr_type_symmetric, name="nl gauge term real part")
236 60 : CALL cp_dbcsr_alloc_block_from_nbl(nl_term(1)%matrix, sab_orb)
237 60 : CALL dbcsr_set(nl_term(1)%matrix, zero)
238 :
239 60 : CALL dbcsr_init_p(nl_term(2)%matrix)
240 : CALL dbcsr_create(nl_term(2)%matrix, template=matrix_s(1, 1)%matrix, &
241 60 : matrix_type=dbcsr_type_antisymmetric, name="nl gauge term imaginary part")
242 60 : CALL cp_dbcsr_alloc_block_from_nbl(nl_term(2)%matrix, sab_orb)
243 60 : CALL dbcsr_set(nl_term(2)%matrix, zero)
244 :
245 60 : CALL velocity_gauge_nl_term(qs_env, nl_term, vec_pot)
246 :
247 120 : DO image = 1, nimages
248 60 : CALL dbcsr_add(matrix_h(1, image)%matrix, nl_term(1)%matrix, one, one)
249 120 : CALL dbcsr_add(matrix_h_im(1, image)%matrix, nl_term(2)%matrix, one, one)
250 : END DO
251 60 : CALL dbcsr_deallocate_matrix_set(nl_term)
252 : END IF
253 : END IF
254 :
255 82 : CALL timestop(handle)
256 :
257 82 : END SUBROUTINE velocity_gauge_ks_matrix
258 :
259 : ! **************************************************************************************************
260 : !> \brief Update the vector potential in the case where a time-dependant
261 : !> electric field is apply.
262 : !> \param qs_env ...
263 : !> \param dft_control ...
264 : ! **************************************************************************************************
265 36 : SUBROUTINE update_vector_potential(qs_env, dft_control)
266 : TYPE(qs_environment_type), INTENT(INOUT), POINTER :: qs_env
267 : TYPE(dft_control_type), INTENT(INOUT), POINTER :: dft_control
268 :
269 : REAL(kind=dp) :: field(3)
270 :
271 36 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
272 144 : dft_control%rtp_control%field = field
273 144 : dft_control%rtp_control%vec_pot = dft_control%rtp_control%vec_pot - field*qs_env%rtp%dt
274 : ! Update the vec_pot_initial value for RTP restart:
275 144 : dft_control%efield_fields(1)%efield%vec_pot_initial = dft_control%rtp_control%vec_pot
276 :
277 36 : END SUBROUTINE update_vector_potential
278 :
279 : ! **************************************************************************************************
280 : !> \brief ...
281 : !> \param qs_env ...
282 : !> \param nl_term ...
283 : !> \param vec_pot ...
284 : ! **************************************************************************************************
285 60 : SUBROUTINE velocity_gauge_nl_term(qs_env, nl_term, vec_pot)
286 : TYPE(qs_environment_type), POINTER :: qs_env
287 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
288 : POINTER :: nl_term
289 : REAL(KIND=dp), DIMENSION(3), INTENT(in) :: vec_pot
290 :
291 : CHARACTER(len=*), PARAMETER :: routiuneN = "velocity_gauge_nl_term"
292 :
293 : INTEGER :: handle, i, iac, iatom, ibc, icol, ikind, &
294 : irow, jatom, jkind, kac, kbc, kkind, &
295 : maxl, maxlgto, maxlppnl, na, natom, &
296 : nb, nkind, np, slot
297 : INTEGER, DIMENSION(3) :: cell_b
298 : LOGICAL :: found
299 : REAL(dp) :: eps_ppnl
300 : REAL(KIND=dp), DIMENSION(3) :: rab
301 60 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: imag_block, real_block
302 60 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint_cos, achint_sin, acint_cos, &
303 60 : acint_sin, bchint_cos, bchint_sin, &
304 60 : bcint_cos, bcint_sin
305 : TYPE(alist_type), POINTER :: alist_cos_ac, alist_cos_bc, &
306 : alist_sin_ac, alist_sin_bc
307 60 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
308 : TYPE(cell_type), POINTER :: cell
309 : TYPE(dft_control_type), POINTER :: dft_control
310 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
311 60 : DIMENSION(:) :: basis_set
312 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
313 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
314 60 : POINTER :: sab_orb, sap_ppnl
315 60 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
316 60 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
317 60 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int_cos, sap_int_sin
318 :
319 : !$ INTEGER(kind=omp_lock_kind), &
320 60 : !$ ALLOCATABLE, DIMENSION(:) :: locks
321 : !$ INTEGER(KIND=int_8) :: iatom8
322 : !$ INTEGER :: lock_num, hash
323 : !$ INTEGER, PARAMETER :: nlock = 501
324 :
325 : MARK_USED(int_8)
326 :
327 60 : CALL timeset(routiuneN, handle)
328 :
329 60 : NULLIFY (sap_ppnl, sab_orb)
330 : CALL get_qs_env(qs_env, &
331 : sap_ppnl=sap_ppnl, &
332 60 : sab_orb=sab_orb)
333 :
334 60 : IF (ASSOCIATED(sap_ppnl)) THEN
335 60 : NULLIFY (qs_kind_set, particle_set, cell, dft_control)
336 : CALL get_qs_env(qs_env, &
337 : dft_control=dft_control, &
338 : qs_kind_set=qs_kind_set, &
339 : particle_set=particle_set, &
340 : cell=cell, &
341 60 : atomic_kind_set=atomic_kind_set)
342 :
343 60 : nkind = SIZE(atomic_kind_set)
344 60 : natom = SIZE(particle_set)
345 60 : eps_ppnl = dft_control%qs_control%eps_ppnl
346 :
347 : CALL get_qs_kind_set(qs_kind_set, &
348 : maxlgto=maxlgto, &
349 60 : maxlppnl=maxlppnl)
350 :
351 60 : maxl = MAX(maxlppnl, maxlgto)
352 60 : CALL init_orbital_pointers(maxl + 1)
353 :
354 : ! initalize sab_int types to store the integrals
355 60 : NULLIFY (sap_int_cos, sap_int_sin)
356 780 : ALLOCATE (sap_int_cos(nkind*nkind), sap_int_sin(nkind*nkind))
357 300 : DO i = 1, SIZE(sap_int_cos)
358 240 : NULLIFY (sap_int_cos(i)%alist, sap_int_cos(i)%asort, sap_int_cos(i)%aindex)
359 240 : sap_int_cos(i)%nalist = 0
360 240 : NULLIFY (sap_int_sin(i)%alist, sap_int_sin(i)%asort, sap_int_sin(i)%aindex)
361 300 : sap_int_sin(i)%nalist = 0
362 : END DO
363 :
364 : ! get basis set
365 300 : ALLOCATE (basis_set(nkind))
366 180 : DO ikind = 1, nkind
367 120 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
368 180 : IF (ASSOCIATED(orb_basis_set)) THEN
369 120 : basis_set(ikind)%gto_basis_set => orb_basis_set
370 : ELSE
371 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
372 : END IF
373 : END DO
374 :
375 : ! calculate exponential integrals
376 : CALL build_sap_exp_ints(sap_int_cos, sap_int_sin, sap_ppnl, qs_kind_set, particle_set, &
377 : cell, kvec=vec_pot, basis_set=basis_set, nkind=nkind, &
378 60 : derivative=.FALSE.)
379 :
380 60 : CALL sap_sort(sap_int_cos)
381 60 : CALL sap_sort(sap_int_sin)
382 :
383 : ! assemble the integrals for the gauge term
384 : !$OMP PARALLEL &
385 : !$OMP DEFAULT (NONE) &
386 : !$OMP SHARED (basis_set, nl_term, sab_orb, sap_int_cos, sap_int_sin, eps_ppnl, locks, nkind, natom) &
387 : !$OMP PRIVATE (real_block, imag_block, acint_cos, achint_cos, bcint_cos, bchint_cos, acint_sin,&
388 : !$OMP achint_sin, bcint_sin, bchint_sin, slot, ikind, jkind, iatom, jatom, cell_b, rab, irow, icol,&
389 : !$OMP found, kkind, iac, ibc, alist_cos_ac, alist_cos_bc, alist_sin_ac, alist_sin_bc, kac, kbc, &
390 60 : !$OMP na, np, nb, iatom8, hash, lock_num)
391 :
392 : !$OMP SINGLE
393 : !$ ALLOCATE (locks(nlock))
394 : !$OMP END SINGLE
395 :
396 : !$OMP DO
397 : !$ DO lock_num = 1, nlock
398 : !$ call omp_init_lock(locks(lock_num))
399 : !$ END DO
400 : !$OMP END DO
401 :
402 : NULLIFY (real_block, imag_block)
403 : NULLIFY (acint_cos, bcint_cos, achint_cos, bchint_cos)
404 : NULLIFY (acint_sin, bcint_sin, achint_sin, bchint_sin)
405 :
406 : ! loop over atom pairs
407 : !$OMP DO SCHEDULE(GUIDED)
408 : DO slot = 1, sab_orb(1)%nl_size
409 : ikind = sab_orb(1)%nlist_task(slot)%ikind
410 : jkind = sab_orb(1)%nlist_task(slot)%jkind
411 : iatom = sab_orb(1)%nlist_task(slot)%iatom
412 : jatom = sab_orb(1)%nlist_task(slot)%jatom
413 : cell_b(:) = sab_orb(1)%nlist_task(slot)%cell
414 : rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
415 :
416 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
417 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
418 :
419 : IF (iatom <= jatom) THEN
420 : irow = iatom
421 : icol = jatom
422 : ELSE
423 : irow = jatom
424 : icol = iatom
425 : END IF
426 :
427 : CALL dbcsr_get_block_p(nl_term(1)%matrix, irow, icol, real_block, found)
428 : CALL dbcsr_get_block_p(nl_term(2)%matrix, irow, icol, imag_block, found)
429 :
430 : IF (ASSOCIATED(real_block) .AND. ASSOCIATED(imag_block)) THEN
431 : ! loop over the <gto_a|ppln_c>h_ij<ppnl_c|gto_b> pairs
432 : DO kkind = 1, nkind
433 : iac = ikind + nkind*(kkind - 1)
434 : ibc = jkind + nkind*(kkind - 1)
435 : IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist)) CYCLE
436 : IF (.NOT. ASSOCIATED(sap_int_cos(ibc)%alist)) CYCLE
437 : IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist)) CYCLE
438 : IF (.NOT. ASSOCIATED(sap_int_sin(ibc)%alist)) CYCLE
439 : CALL get_alist(sap_int_cos(iac), alist_cos_ac, iatom)
440 : CALL get_alist(sap_int_cos(ibc), alist_cos_bc, jatom)
441 : CALL get_alist(sap_int_sin(iac), alist_sin_ac, iatom)
442 : CALL get_alist(sap_int_sin(ibc), alist_sin_bc, jatom)
443 : IF (.NOT. ASSOCIATED(alist_cos_ac)) CYCLE
444 : IF (.NOT. ASSOCIATED(alist_cos_bc)) CYCLE
445 : IF (.NOT. ASSOCIATED(alist_sin_ac)) CYCLE
446 : IF (.NOT. ASSOCIATED(alist_sin_bc)) CYCLE
447 :
448 : ! only use cos for indexing, as cos and sin integrals are constructed by the same routine
449 : ! in the same way
450 : DO kac = 1, alist_cos_ac%nclist
451 : DO kbc = 1, alist_cos_bc%nclist
452 : ! the next two ifs should be the same for sine integrals
453 : IF (alist_cos_ac%clist(kac)%catom /= alist_cos_bc%clist(kbc)%catom) CYCLE
454 : IF (ALL(cell_b + alist_cos_bc%clist(kbc)%cell - alist_cos_ac%clist(kac)%cell == 0)) THEN
455 : ! screening
456 : IF (alist_cos_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
457 : .AND. alist_cos_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl &
458 : .AND. alist_sin_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
459 : .AND. alist_sin_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
460 :
461 : acint_cos => alist_cos_ac%clist(kac)%acint
462 : bcint_cos => alist_cos_bc%clist(kbc)%acint
463 : achint_cos => alist_cos_ac%clist(kac)%achint
464 : bchint_cos => alist_cos_bc%clist(kbc)%achint
465 : acint_sin => alist_sin_ac%clist(kac)%acint
466 : bcint_sin => alist_sin_bc%clist(kbc)%acint
467 : achint_sin => alist_sin_ac%clist(kac)%achint
468 : bchint_sin => alist_sin_bc%clist(kbc)%achint
469 :
470 : na = SIZE(acint_cos, 1)
471 : np = SIZE(acint_cos, 2)
472 : nb = SIZE(bcint_cos, 1)
473 : !$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
474 : !$ hash = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
475 : !$ CALL omp_set_lock(locks(hash))
476 : IF (iatom <= jatom) THEN
477 : ! cos*cos + sin*sin
478 : real_block(1:na, 1:nb) = real_block(1:na, 1:nb) + &
479 : MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, 1))) + &
480 : MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, 1)))
481 : ! sin * cos - cos * sin
482 : imag_block(1:na, 1:nb) = imag_block(1:na, 1:nb) - &
483 : MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, 1))) + &
484 : MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, 1)))
485 : ELSE
486 : ! cos*cos + sin*sin
487 : real_block(1:nb, 1:na) = real_block(1:nb, 1:na) + &
488 : MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, 1))) + &
489 : MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, 1)))
490 : ! sin * cos - cos * sin
491 : imag_block(1:nb, 1:na) = imag_block(1:nb, 1:na) - &
492 : MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, 1))) + &
493 : MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, 1)))
494 :
495 : END IF
496 : !$ CALL omp_unset_lock(locks(hash))
497 : EXIT
498 : END IF
499 : END DO
500 : END DO
501 : END DO
502 : END IF
503 :
504 : END DO
505 :
506 : !$OMP DO
507 : !$ DO lock_num = 1, nlock
508 : !$ call omp_destroy_lock(locks(lock_num))
509 : !$ END DO
510 : !$OMP END DO
511 :
512 : !$OMP SINGLE
513 : !$ DEALLOCATE (locks)
514 : !$OMP END SINGLE NOWAIT
515 :
516 : !$OMP END PARALLEL
517 60 : CALL release_sap_int(sap_int_cos)
518 60 : CALL release_sap_int(sap_int_sin)
519 :
520 120 : DEALLOCATE (basis_set)
521 : END IF
522 :
523 60 : CALL timestop(handle)
524 :
525 120 : END SUBROUTINE velocity_gauge_nl_term
526 :
527 : ! **************************************************************************************************
528 : !> \brief Calculate the commutator [Vnl~, r_d] (d=x,y,z) of the gauge-transformed nonlocal
529 : !> pseudopotential with the position operator, needed for the nonlocal correction
530 : !> to the integrated electronic current in RTP/EMD.
531 : !> The moment is taken relative to each projector center, which makes the result
532 : !> origin-independent (the constant shifts cancel exactly in the commutator).
533 : !> For vec_pot = 0 the plain commutator [Vnl, r_d] is obtained (matrix_cim = 0).
534 : !> \param qs_env ...
535 : !> \param matrix_cre real part of the commutator, antisymmetric, dimension(3)
536 : !> \param matrix_cim imaginary part of the commutator, symmetric, dimension(3)
537 : !> \param vec_pot vector potential defining the gauge phase, use 0 for no gauge transform
538 : ! **************************************************************************************************
539 92 : SUBROUTINE velocity_gauge_nl_commutator(qs_env, matrix_cre, matrix_cim, vec_pot)
540 : TYPE(qs_environment_type), POINTER :: qs_env
541 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
542 : POINTER :: matrix_cre, matrix_cim
543 : REAL(KIND=dp), DIMENSION(3), INTENT(in) :: vec_pot
544 :
545 : CHARACTER(len=*), PARAMETER :: routiuneN = "velocity_gauge_nl_commutator"
546 :
547 : INTEGER :: handle, i, iac, iatom, ibc, icol, idir, &
548 : ikind, irow, jatom, jkind, kac, kbc, &
549 : kkind, m, maxl, maxlgto, maxlppnl, na, &
550 : natom, nb, nkind, np, slot
551 : INTEGER, DIMENSION(3) :: cell_b
552 : LOGICAL :: found
553 : REAL(dp) :: eps_ppnl
554 : REAL(KIND=dp), DIMENSION(3) :: rab
555 92 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: imag_block, real_block
556 92 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint_cos, achint_sin, acint_cos, &
557 92 : acint_sin, bchint_cos, bchint_sin, &
558 92 : bcint_cos, bcint_sin
559 : TYPE(alist_type), POINTER :: alist_cos_ac, alist_cos_bc, &
560 : alist_sin_ac, alist_sin_bc
561 92 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
562 : TYPE(cell_type), POINTER :: cell
563 : TYPE(dft_control_type), POINTER :: dft_control
564 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
565 92 : DIMENSION(:) :: basis_set
566 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
567 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
568 92 : POINTER :: sab_orb, sap_ppnl
569 92 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
570 92 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
571 92 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int_cos, sap_int_sin
572 :
573 : !$ INTEGER(kind=omp_lock_kind), &
574 92 : !$ ALLOCATABLE, DIMENSION(:) :: locks
575 : !$ INTEGER(KIND=int_8) :: iatom8
576 : !$ INTEGER :: lock_num, hash
577 : !$ INTEGER, PARAMETER :: nlock = 501
578 :
579 : MARK_USED(int_8)
580 :
581 92 : CALL timeset(routiuneN, handle)
582 :
583 92 : NULLIFY (sap_ppnl, sab_orb)
584 : CALL get_qs_env(qs_env, &
585 : sap_ppnl=sap_ppnl, &
586 92 : sab_orb=sab_orb)
587 :
588 92 : IF (ASSOCIATED(sap_ppnl)) THEN
589 92 : NULLIFY (qs_kind_set, particle_set, cell, dft_control)
590 : CALL get_qs_env(qs_env, &
591 : dft_control=dft_control, &
592 : qs_kind_set=qs_kind_set, &
593 : particle_set=particle_set, &
594 : cell=cell, &
595 92 : atomic_kind_set=atomic_kind_set)
596 :
597 92 : nkind = SIZE(atomic_kind_set)
598 92 : natom = SIZE(particle_set)
599 92 : eps_ppnl = dft_control%qs_control%eps_ppnl
600 :
601 : CALL get_qs_kind_set(qs_kind_set, &
602 : maxlgto=maxlgto, &
603 92 : maxlppnl=maxlppnl)
604 :
605 : ! the moment integrals need one higher angular momentum on the basis side,
606 : ! and cossin internally uses one more
607 92 : maxl = MAX(maxlppnl, maxlgto)
608 92 : CALL init_orbital_pointers(maxl + 2)
609 :
610 : ! initalize sab_int types to store the integrals
611 92 : NULLIFY (sap_int_cos, sap_int_sin)
612 1292 : ALLOCATE (sap_int_cos(nkind*nkind), sap_int_sin(nkind*nkind))
613 508 : DO i = 1, SIZE(sap_int_cos)
614 416 : NULLIFY (sap_int_cos(i)%alist, sap_int_cos(i)%asort, sap_int_cos(i)%aindex)
615 416 : sap_int_cos(i)%nalist = 0
616 416 : NULLIFY (sap_int_sin(i)%alist, sap_int_sin(i)%asort, sap_int_sin(i)%aindex)
617 508 : sap_int_sin(i)%nalist = 0
618 : END DO
619 :
620 : ! get basis set
621 468 : ALLOCATE (basis_set(nkind))
622 284 : DO ikind = 1, nkind
623 192 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
624 284 : IF (ASSOCIATED(orb_basis_set)) THEN
625 192 : basis_set(ikind)%gto_basis_set => orb_basis_set
626 : ELSE
627 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
628 : END IF
629 : END DO
630 :
631 : ! calculate exponential integrals, plain (slot 1) and moment-weighted (slots 2:4)
632 : CALL build_sap_exp_ints(sap_int_cos, sap_int_sin, sap_ppnl, qs_kind_set, particle_set, &
633 : cell, kvec=vec_pot, basis_set=basis_set, nkind=nkind, &
634 92 : derivative=.FALSE., moment=.TRUE.)
635 :
636 92 : CALL sap_sort(sap_int_cos)
637 92 : CALL sap_sort(sap_int_sin)
638 :
639 : ! assemble the commutator matrices
640 : !$OMP PARALLEL &
641 : !$OMP DEFAULT (NONE) &
642 : !$OMP SHARED (basis_set, matrix_cre, matrix_cim, sab_orb, sap_int_cos, sap_int_sin, eps_ppnl, locks, nkind, natom) &
643 : !$OMP PRIVATE (real_block, imag_block, acint_cos, achint_cos, bcint_cos, bchint_cos, acint_sin,&
644 : !$OMP achint_sin, bcint_sin, bchint_sin, slot, ikind, jkind, iatom, jatom, cell_b, rab, irow, icol,&
645 : !$OMP found, kkind, iac, ibc, alist_cos_ac, alist_cos_bc, alist_sin_ac, alist_sin_bc, kac, kbc, &
646 92 : !$OMP na, np, nb, idir, m, iatom8, hash, lock_num)
647 :
648 : !$OMP SINGLE
649 : !$ ALLOCATE (locks(nlock))
650 : !$OMP END SINGLE
651 :
652 : !$OMP DO
653 : !$ DO lock_num = 1, nlock
654 : !$ call omp_init_lock(locks(lock_num))
655 : !$ END DO
656 : !$OMP END DO
657 :
658 : NULLIFY (real_block, imag_block)
659 : NULLIFY (acint_cos, bcint_cos, achint_cos, bchint_cos)
660 : NULLIFY (acint_sin, bcint_sin, achint_sin, bchint_sin)
661 :
662 : ! loop over atom pairs
663 : !$OMP DO SCHEDULE(GUIDED)
664 : DO slot = 1, sab_orb(1)%nl_size
665 : ikind = sab_orb(1)%nlist_task(slot)%ikind
666 : jkind = sab_orb(1)%nlist_task(slot)%jkind
667 : iatom = sab_orb(1)%nlist_task(slot)%iatom
668 : jatom = sab_orb(1)%nlist_task(slot)%jatom
669 : cell_b(:) = sab_orb(1)%nlist_task(slot)%cell
670 : rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
671 :
672 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
673 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
674 :
675 : IF (iatom <= jatom) THEN
676 : irow = iatom
677 : icol = jatom
678 : ELSE
679 : irow = jatom
680 : icol = iatom
681 : END IF
682 :
683 : ! loop over the <gto_a|ppln_c>h_ij<ppnl_c|gto_b> pairs
684 : DO kkind = 1, nkind
685 : iac = ikind + nkind*(kkind - 1)
686 : ibc = jkind + nkind*(kkind - 1)
687 : IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist)) CYCLE
688 : IF (.NOT. ASSOCIATED(sap_int_cos(ibc)%alist)) CYCLE
689 : IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist)) CYCLE
690 : IF (.NOT. ASSOCIATED(sap_int_sin(ibc)%alist)) CYCLE
691 : CALL get_alist(sap_int_cos(iac), alist_cos_ac, iatom)
692 : CALL get_alist(sap_int_cos(ibc), alist_cos_bc, jatom)
693 : CALL get_alist(sap_int_sin(iac), alist_sin_ac, iatom)
694 : CALL get_alist(sap_int_sin(ibc), alist_sin_bc, jatom)
695 : IF (.NOT. ASSOCIATED(alist_cos_ac)) CYCLE
696 : IF (.NOT. ASSOCIATED(alist_cos_bc)) CYCLE
697 : IF (.NOT. ASSOCIATED(alist_sin_ac)) CYCLE
698 : IF (.NOT. ASSOCIATED(alist_sin_bc)) CYCLE
699 :
700 : ! only use cos for indexing, as cos and sin integrals are constructed by the same routine
701 : ! in the same way
702 : DO kac = 1, alist_cos_ac%nclist
703 : DO kbc = 1, alist_cos_bc%nclist
704 : ! the next two ifs should be the same for sine integrals
705 : IF (alist_cos_ac%clist(kac)%catom /= alist_cos_bc%clist(kbc)%catom) CYCLE
706 : IF (ALL(cell_b + alist_cos_bc%clist(kbc)%cell - alist_cos_ac%clist(kac)%cell == 0)) THEN
707 : ! screening
708 : IF (alist_cos_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
709 : .AND. alist_cos_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl &
710 : .AND. alist_sin_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
711 : .AND. alist_sin_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
712 :
713 : acint_cos => alist_cos_ac%clist(kac)%acint
714 : bcint_cos => alist_cos_bc%clist(kbc)%acint
715 : achint_cos => alist_cos_ac%clist(kac)%achint
716 : bchint_cos => alist_cos_bc%clist(kbc)%achint
717 : acint_sin => alist_sin_ac%clist(kac)%acint
718 : bcint_sin => alist_sin_bc%clist(kbc)%acint
719 : achint_sin => alist_sin_ac%clist(kac)%achint
720 : bchint_sin => alist_sin_bc%clist(kbc)%achint
721 :
722 : na = SIZE(acint_cos, 1)
723 : np = SIZE(acint_cos, 2)
724 : nb = SIZE(bcint_cos, 1)
725 : !$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
726 : !$ hash = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
727 : !$ CALL omp_set_lock(locks(hash))
728 : DO idir = 1, 3
729 : m = 1 + idir
730 : CALL dbcsr_get_block_p(matrix_cre(idir)%matrix, irow, icol, real_block, found)
731 : CALL dbcsr_get_block_p(matrix_cim(idir)%matrix, irow, icol, imag_block, found)
732 : IF (.NOT. (ASSOCIATED(real_block) .AND. ASSOCIATED(imag_block))) CYCLE
733 : IF (iatom <= jatom) THEN
734 : ! <a|Vnl~ (r_d - R_c)|b> - <a|(r_d - R_c) Vnl~|b>
735 : real_block(1:na, 1:nb) = real_block(1:na, 1:nb) + &
736 : MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, m))) + &
737 : MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, m))) - &
738 : MATMUL(achint_cos(1:na, 1:np, m), TRANSPOSE(bcint_cos(1:nb, 1:np, 1))) - &
739 : MATMUL(achint_sin(1:na, 1:np, m), TRANSPOSE(bcint_sin(1:nb, 1:np, 1)))
740 : imag_block(1:na, 1:nb) = imag_block(1:na, 1:nb) + &
741 : MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, m))) - &
742 : MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, m))) - &
743 : MATMUL(achint_cos(1:na, 1:np, m), TRANSPOSE(bcint_sin(1:nb, 1:np, 1))) + &
744 : MATMUL(achint_sin(1:na, 1:np, m), TRANSPOSE(bcint_cos(1:nb, 1:np, 1)))
745 : ELSE
746 : ! stored block is (jatom, iatom): same formulas with the roles of a and b swapped
747 : real_block(1:nb, 1:na) = real_block(1:nb, 1:na) + &
748 : MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, m))) + &
749 : MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, m))) - &
750 : MATMUL(bchint_cos(1:nb, 1:np, m), TRANSPOSE(acint_cos(1:na, 1:np, 1))) - &
751 : MATMUL(bchint_sin(1:nb, 1:np, m), TRANSPOSE(acint_sin(1:na, 1:np, 1)))
752 : imag_block(1:nb, 1:na) = imag_block(1:nb, 1:na) + &
753 : MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, m))) - &
754 : MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, m))) - &
755 : MATMUL(bchint_cos(1:nb, 1:np, m), TRANSPOSE(acint_sin(1:na, 1:np, 1))) + &
756 : MATMUL(bchint_sin(1:nb, 1:np, m), TRANSPOSE(acint_cos(1:na, 1:np, 1)))
757 : END IF
758 : END DO
759 : !$ CALL omp_unset_lock(locks(hash))
760 : EXIT
761 : END IF
762 : END DO
763 : END DO
764 : END DO
765 :
766 : END DO
767 :
768 : !$OMP DO
769 : !$ DO lock_num = 1, nlock
770 : !$ call omp_destroy_lock(locks(lock_num))
771 : !$ END DO
772 : !$OMP END DO
773 :
774 : !$OMP SINGLE
775 : !$ DEALLOCATE (locks)
776 : !$OMP END SINGLE NOWAIT
777 :
778 : !$OMP END PARALLEL
779 92 : CALL release_sap_int(sap_int_cos)
780 92 : CALL release_sap_int(sap_int_sin)
781 :
782 184 : DEALLOCATE (basis_set)
783 : END IF
784 :
785 92 : CALL timestop(handle)
786 :
787 184 : END SUBROUTINE velocity_gauge_nl_commutator
788 :
789 : ! **************************************************************************************************
790 : !> \brief calculate <a|sin/cos|p> integrals and store in sap_int_type
791 : !> adapted from build_sap_ints
792 : !> Do this on each MPI task as the integrals need to be available globally.
793 : !> Might be faster than communicating as the integrals are obtained analytically.
794 : !> If asked, compute <da/dRa|sin/cos|p>
795 : !> If moment is requested, compute the moment-weighted integrals
796 : !> <a|sin/cos (r_d - R_c,d)|p> (d=x,y,z) in slots 2:4, with the moment taken
797 : !> relative to the projector center R_c (the frame origin used here)
798 : !> \param sap_int_cos ...
799 : !> \param sap_int_sin ...
800 : !> \param sap_ppnl ...
801 : !> \param qs_kind_set ...
802 : !> \param particle_set ...
803 : !> \param cell ...
804 : !> \param kvec ...
805 : !> \param basis_set ...
806 : !> \param nkind ...
807 : !> \param derivative ...
808 : !> \param moment ...
809 : ! **************************************************************************************************
810 184 : SUBROUTINE build_sap_exp_ints(sap_int_cos, sap_int_sin, sap_ppnl, qs_kind_set, particle_set, cell, &
811 184 : kvec, basis_set, nkind, derivative, moment)
812 : TYPE(sap_int_type), DIMENSION(:), INTENT(INOUT), &
813 : POINTER :: sap_int_cos, sap_int_sin
814 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
815 : INTENT(IN), POINTER :: sap_ppnl
816 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
817 : POINTER :: qs_kind_set
818 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
819 : POINTER :: particle_set
820 : TYPE(cell_type), INTENT(IN), POINTER :: cell
821 : REAL(KIND=dp), DIMENSION(3), INTENT(in) :: kvec
822 : TYPE(gto_basis_set_p_type), DIMENSION(:), &
823 : INTENT(IN) :: basis_set
824 : INTEGER, INTENT(IN) :: nkind
825 : LOGICAL, INTENT(IN) :: derivative
826 : LOGICAL, INTENT(IN), OPTIONAL :: moment
827 :
828 : CHARACTER(len=*), PARAMETER :: routiuneN = "build_sap_exp_ints"
829 :
830 : INTEGER :: ax, ay, az, co, cox, coy, coz, handle, i, iac, iatom, idir, ikind, ilist, iob, &
831 : ioe, ipgf, iset, jneighbor, katom, kkind, l, la, lc_max, lc_min, ldai, ldax, ldints, &
832 : lppnl, maxco, maxl, maxlgto, maxlppnl, maxppnl, maxsgf, na, nb, nco_b, nco_e, ncoa, ncoc, &
833 : nlist, nneighbor, np, nppnl, nprjc, nseta, nsgfa, nsl, prjc, sgfa, slot
834 : INTEGER, DIMENSION(3) :: cell_c
835 184 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
836 184 : nsgf_seta
837 184 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
838 : LOGICAL :: dogth, my_moment
839 : REAL(KIND=dp) :: dac, ppnl_radius
840 184 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ai_work_cos, ai_work_mcos, ai_work_msin, &
841 184 : ai_work_sin, work_cos, work_sin
842 184 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work_dcos, ai_work_dsin, work_dcos, &
843 184 : work_dsin
844 : REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc
845 : REAL(KIND=dp), DIMENSION(3) :: ra, rac, raf, rc, rcf
846 184 : REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_ppnl, set_radius_a
847 184 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, rpgfa, sphi_a, vprj_ppnl, zeta
848 : TYPE(clist_type), POINTER :: clist, clist_sin
849 184 : TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
850 : TYPE(gth_potential_type), POINTER :: gth_potential
851 184 : TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential
852 : TYPE(sgp_potential_type), POINTER :: sgp_potential
853 :
854 184 : CALL timeset(routiuneN, handle)
855 :
856 184 : my_moment = .FALSE.
857 184 : IF (PRESENT(moment)) my_moment = moment
858 : ! derivative and moment integrals use the same storage slots 2:4
859 184 : CPASSERT(.NOT. (derivative .AND. my_moment))
860 184 : nsl = 1
861 184 : IF (derivative .OR. my_moment) nsl = 4
862 :
863 : CALL get_qs_kind_set(qs_kind_set, &
864 : maxco=maxco, &
865 : maxlppnl=maxlppnl, &
866 : maxppnl=maxppnl, &
867 : maxsgf=maxsgf, &
868 184 : maxlgto=maxlgto)
869 :
870 : ! maximum dimensions for allocations
871 184 : maxl = MAX(maxlppnl, maxlgto)
872 184 : ldints = MAX(maxco, ncoset(maxlppnl), maxsgf, maxppnl)
873 184 : ldai = ncoset(maxl + 1)
874 :
875 : ! leading dimension of the la_max+1 integrals needed in moment mode
876 184 : ldax = 0
877 184 : IF (my_moment) THEN
878 284 : DO i = 1, SIZE(basis_set)
879 192 : IF (.NOT. ASSOCIATED(basis_set(i)%gto_basis_set)) CYCLE
880 572 : DO iset = 1, basis_set(i)%gto_basis_set%nset
881 : ldax = MAX(ldax, basis_set(i)%gto_basis_set%npgf(iset)* &
882 480 : ncoset(basis_set(i)%gto_basis_set%lmax(iset) + 1))
883 : END DO
884 : END DO
885 : END IF
886 :
887 : !set up direct access to basis and potential
888 184 : NULLIFY (gpotential, spotential)
889 1672 : ALLOCATE (gpotential(nkind), spotential(nkind))
890 560 : DO ikind = 1, nkind
891 376 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
892 376 : NULLIFY (gpotential(ikind)%gth_potential)
893 376 : NULLIFY (spotential(ikind)%sgp_potential)
894 560 : IF (ASSOCIATED(gth_potential)) THEN
895 376 : gpotential(ikind)%gth_potential => gth_potential
896 0 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
897 0 : spotential(ikind)%sgp_potential => sgp_potential
898 : END IF
899 : END DO
900 :
901 : !allocate sap int
902 184 : NULLIFY (clist)
903 2364 : DO slot = 1, sap_ppnl(1)%nl_size
904 :
905 2180 : ikind = sap_ppnl(1)%nlist_task(slot)%ikind
906 2180 : kkind = sap_ppnl(1)%nlist_task(slot)%jkind
907 2180 : iatom = sap_ppnl(1)%nlist_task(slot)%iatom
908 2180 : katom = sap_ppnl(1)%nlist_task(slot)%jatom
909 2180 : nlist = sap_ppnl(1)%nlist_task(slot)%nlist
910 2180 : ilist = sap_ppnl(1)%nlist_task(slot)%ilist
911 2180 : nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
912 :
913 2180 : iac = ikind + nkind*(kkind - 1)
914 2180 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
915 2180 : IF (.NOT. ASSOCIATED(gpotential(kkind)%gth_potential) .AND. &
916 : .NOT. ASSOCIATED(spotential(kkind)%sgp_potential)) CYCLE
917 2180 : IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist)) THEN
918 416 : sap_int_cos(iac)%a_kind = ikind
919 416 : sap_int_cos(iac)%p_kind = kkind
920 416 : sap_int_cos(iac)%nalist = nlist
921 2256 : ALLOCATE (sap_int_cos(iac)%alist(nlist))
922 1424 : DO i = 1, nlist
923 1008 : NULLIFY (sap_int_cos(iac)%alist(i)%clist)
924 1008 : sap_int_cos(iac)%alist(i)%aatom = 0
925 1424 : sap_int_cos(iac)%alist(i)%nclist = 0
926 : END DO
927 : END IF
928 2180 : IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist(ilist)%clist)) THEN
929 1008 : sap_int_cos(iac)%alist(ilist)%aatom = iatom
930 1008 : sap_int_cos(iac)%alist(ilist)%nclist = nneighbor
931 12260 : ALLOCATE (sap_int_cos(iac)%alist(ilist)%clist(nneighbor))
932 3188 : DO i = 1, nneighbor
933 2180 : clist => sap_int_cos(iac)%alist(ilist)%clist(i)
934 2180 : clist%catom = 0
935 2180 : NULLIFY (clist%acint)
936 2180 : NULLIFY (clist%achint)
937 3188 : NULLIFY (clist%sgf_list)
938 : END DO
939 : END IF
940 2180 : IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist)) THEN
941 416 : sap_int_sin(iac)%a_kind = ikind
942 416 : sap_int_sin(iac)%p_kind = kkind
943 416 : sap_int_sin(iac)%nalist = nlist
944 2256 : ALLOCATE (sap_int_sin(iac)%alist(nlist))
945 1424 : DO i = 1, nlist
946 1008 : NULLIFY (sap_int_sin(iac)%alist(i)%clist)
947 1008 : sap_int_sin(iac)%alist(i)%aatom = 0
948 1424 : sap_int_sin(iac)%alist(i)%nclist = 0
949 : END DO
950 : END IF
951 2364 : IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist(ilist)%clist)) THEN
952 1008 : sap_int_sin(iac)%alist(ilist)%aatom = iatom
953 1008 : sap_int_sin(iac)%alist(ilist)%nclist = nneighbor
954 12260 : ALLOCATE (sap_int_sin(iac)%alist(ilist)%clist(nneighbor))
955 3188 : DO i = 1, nneighbor
956 2180 : clist => sap_int_sin(iac)%alist(ilist)%clist(i)
957 2180 : clist%catom = 0
958 2180 : NULLIFY (clist%acint)
959 2180 : NULLIFY (clist%achint)
960 3188 : NULLIFY (clist%sgf_list)
961 : END DO
962 : END IF
963 : END DO
964 :
965 : ! actual calculation of the integrals <a|cos|p> and <a|sin|p>
966 : ! allocate temporary storage using maximum dimensions
967 :
968 : !$OMP PARALLEL &
969 : !$OMP DEFAULT (NONE) &
970 : !$OMP SHARED (basis_set, gpotential, ncoset, coset, sap_ppnl, sap_int_cos, sap_int_sin, nkind, &
971 : !$OMP ldints, ldax, maxco, nco, cell, particle_set, kvec, derivative, my_moment, nsl) &
972 : !$OMP PRIVATE (slot, ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, &
973 : !$OMP cell_c, rac, dac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta,&
974 : !$OMP rpgfa, set_radius_a, sphi_a, zeta, alpha_ppnl, cprj, lppnl, nppnl, nprj_ppnl,&
975 : !$OMP ppnl_radius, vprj_ppnl, clist, clist_sin, ra, rc, ncoa, sgfa, prjc, work_cos, work_sin,&
976 : !$OMP nprjc, rprjc, lc_max, lc_min, zetc, ncoc, ai_work_sin, ai_work_cos, na, nb, np, dogth, &
977 : !$OMP raf, rcf, work_dcos, work_dsin, ai_work_dcos, ai_work_dsin, idir, &
978 184 : !$OMP ai_work_mcos, ai_work_msin, ipgf, iob, ioe, nco_b, nco_e, la, ax, ay, az, co, cox, coy, coz)
979 :
980 : ALLOCATE (work_cos(ldints, ldints), work_sin(ldints, ldints))
981 : ALLOCATE (ai_work_cos(maxco, maxco), ai_work_sin(maxco, maxco))
982 : IF (derivative .OR. my_moment) THEN
983 : ALLOCATE (work_dcos(ldints, ldints, 3), work_dsin(ldints, ldints, 3))
984 : ALLOCATE (ai_work_dcos(maxco, maxco, 3), ai_work_dsin(maxco, maxco, 3))
985 : END IF
986 : IF (my_moment) THEN
987 : ALLOCATE (ai_work_mcos(ldax, maxco), ai_work_msin(ldax, maxco))
988 : ai_work_mcos = 0.0_dp
989 : ai_work_msin = 0.0_dp
990 : END IF
991 : work_cos = 0.0_dp
992 : work_sin = 0.0_dp
993 : ai_work_cos = 0.0_dp
994 : ai_work_sin = 0.0_dp
995 : IF (derivative .OR. my_moment) THEN
996 : ai_work_dcos = 0.0_dp
997 : ai_work_dsin = 0.0_dp
998 : END IF
999 : dogth = .FALSE.
1000 :
1001 : NULLIFY (first_sgfa, la_max, la_min, npgfa, nsgf_seta, rpgfa, set_radius_a, sphi_a, zeta)
1002 : NULLIFY (alpha_ppnl, cprj, nprj_ppnl, vprj_ppnl)
1003 : NULLIFY (clist, clist_sin)
1004 :
1005 : !$OMP DO SCHEDULE(GUIDED)
1006 : DO slot = 1, sap_ppnl(1)%nl_size
1007 : ikind = sap_ppnl(1)%nlist_task(slot)%ikind
1008 : kkind = sap_ppnl(1)%nlist_task(slot)%jkind
1009 : iatom = sap_ppnl(1)%nlist_task(slot)%iatom
1010 : katom = sap_ppnl(1)%nlist_task(slot)%jatom
1011 : nlist = sap_ppnl(1)%nlist_task(slot)%nlist
1012 : ilist = sap_ppnl(1)%nlist_task(slot)%ilist
1013 : nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
1014 : jneighbor = sap_ppnl(1)%nlist_task(slot)%inode
1015 : cell_c(:) = sap_ppnl(1)%nlist_task(slot)%cell(:)
1016 : rac(1:3) = sap_ppnl(1)%nlist_task(slot)%r(1:3)
1017 : dac = NORM2(rac)
1018 :
1019 : iac = ikind + nkind*(kkind - 1)
1020 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
1021 : ! get definition of gto basis set
1022 : first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf
1023 : la_max => basis_set(ikind)%gto_basis_set%lmax
1024 : la_min => basis_set(ikind)%gto_basis_set%lmin
1025 : npgfa => basis_set(ikind)%gto_basis_set%npgf
1026 : nseta = basis_set(ikind)%gto_basis_set%nset
1027 : nsgfa = basis_set(ikind)%gto_basis_set%nsgf
1028 : nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set
1029 : rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius
1030 : set_radius_a => basis_set(ikind)%gto_basis_set%set_radius
1031 : sphi_a => basis_set(ikind)%gto_basis_set%sphi
1032 : zeta => basis_set(ikind)%gto_basis_set%zet
1033 :
1034 : IF (ASSOCIATED(gpotential(kkind)%gth_potential)) THEN
1035 : ! GTH potential
1036 : dogth = .TRUE.
1037 : alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl
1038 : cprj => gpotential(kkind)%gth_potential%cprj
1039 : lppnl = gpotential(kkind)%gth_potential%lppnl
1040 : nppnl = gpotential(kkind)%gth_potential%nppnl
1041 : nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl
1042 : ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius
1043 : vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl
1044 : ELSE
1045 : CYCLE
1046 : END IF
1047 :
1048 : clist => sap_int_cos(iac)%alist(ilist)%clist(jneighbor)
1049 : clist_sin => sap_int_sin(iac)%alist(ilist)%clist(jneighbor)
1050 :
1051 : clist%catom = katom
1052 : clist%cell = cell_c
1053 : clist%rac = rac
1054 : clist_sin%catom = katom
1055 : clist_sin%cell = cell_c
1056 : clist_sin%rac = rac
1057 :
1058 : ALLOCATE (clist%acint(nsgfa, nppnl, nsl), clist%achint(nsgfa, nppnl, nsl))
1059 : clist%acint = 0.0_dp
1060 : clist%achint = 0.0_dp
1061 : clist%nsgf_cnt = 0
1062 :
1063 : ALLOCATE (clist_sin%acint(nsgfa, nppnl, nsl), clist_sin%achint(nsgfa, nppnl, nsl))
1064 : clist_sin%acint = 0.0_dp
1065 : clist_sin%achint = 0.0_dp
1066 : clist_sin%nsgf_cnt = 0
1067 :
1068 : ! reference point at zero
1069 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
1070 : rc(:) = ra + rac
1071 :
1072 : ! reference point at pseudized atom
1073 : raf(:) = ra - rc
1074 : rcf(:) = 0._dp
1075 :
1076 : DO iset = 1, nseta
1077 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1078 : sgfa = first_sgfa(1, iset)
1079 : IF (dogth) THEN
1080 : prjc = 1
1081 : work_cos = 0.0_dp
1082 : work_sin = 0.0_dp
1083 : DO l = 0, lppnl
1084 : nprjc = nprj_ppnl(l)*nco(l)
1085 : IF (nprjc == 0) CYCLE
1086 : rprjc(1) = ppnl_radius
1087 : IF (set_radius_a(iset) + rprjc(1) < dac) CYCLE
1088 : lc_max = l + 2*(nprj_ppnl(l) - 1)
1089 : lc_min = l
1090 : zetc(1) = alpha_ppnl(l)
1091 : ncoc = ncoset(lc_max)
1092 :
1093 : IF (derivative) THEN
1094 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1095 : lc_max, 1, zetc, rprjc, lc_min, raf, rcf, kvec, ai_work_cos, ai_work_sin, &
1096 : dcosab=ai_work_dcos, dsinab=ai_work_dsin)
1097 : ELSE IF (my_moment) THEN
1098 : ! moment-weighted integrals <a|cos/sin (r_d - R_c,d)|p> with the moment relative
1099 : ! to the projector center (the frame origin): apply the shift relation
1100 : ! (r_d - R_c,d) phi_a = phi_(a+1_d) + raf(d) phi_a
1101 : ! to the integrals computed with one higher angular momentum on the basis side
1102 : CALL cossin(la_max(iset) + 1, npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1103 : lc_max, 1, zetc, rprjc, lc_min, raf, rcf, kvec, ai_work_mcos, ai_work_msin)
1104 : nco_b = ncoset(la_max(iset))
1105 : nco_e = ncoset(la_max(iset) + 1)
1106 : DO ipgf = 1, npgfa(iset)
1107 : iob = (ipgf - 1)*nco_b
1108 : ioe = (ipgf - 1)*nco_e
1109 : ai_work_cos(iob + 1:iob + nco_b, 1:ncoc) = ai_work_mcos(ioe + 1:ioe + nco_b, 1:ncoc)
1110 : ai_work_sin(iob + 1:iob + nco_b, 1:ncoc) = ai_work_msin(ioe + 1:ioe + nco_b, 1:ncoc)
1111 : DO la = 0, la_max(iset)
1112 : DO ax = 0, la
1113 : DO ay = 0, la - ax
1114 : az = la - ax - ay
1115 : co = coset(ax, ay, az)
1116 : cox = coset(ax + 1, ay, az)
1117 : coy = coset(ax, ay + 1, az)
1118 : coz = coset(ax, ay, az + 1)
1119 : ai_work_dcos(iob + co, 1:ncoc, 1) = ai_work_mcos(ioe + cox, 1:ncoc) + &
1120 : raf(1)*ai_work_mcos(ioe + co, 1:ncoc)
1121 : ai_work_dcos(iob + co, 1:ncoc, 2) = ai_work_mcos(ioe + coy, 1:ncoc) + &
1122 : raf(2)*ai_work_mcos(ioe + co, 1:ncoc)
1123 : ai_work_dcos(iob + co, 1:ncoc, 3) = ai_work_mcos(ioe + coz, 1:ncoc) + &
1124 : raf(3)*ai_work_mcos(ioe + co, 1:ncoc)
1125 : ai_work_dsin(iob + co, 1:ncoc, 1) = ai_work_msin(ioe + cox, 1:ncoc) + &
1126 : raf(1)*ai_work_msin(ioe + co, 1:ncoc)
1127 : ai_work_dsin(iob + co, 1:ncoc, 2) = ai_work_msin(ioe + coy, 1:ncoc) + &
1128 : raf(2)*ai_work_msin(ioe + co, 1:ncoc)
1129 : ai_work_dsin(iob + co, 1:ncoc, 3) = ai_work_msin(ioe + coz, 1:ncoc) + &
1130 : raf(3)*ai_work_msin(ioe + co, 1:ncoc)
1131 : END DO
1132 : END DO
1133 : END DO
1134 : END DO
1135 : ELSE
1136 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1137 : lc_max, 1, zetc, rprjc, lc_min, raf, rcf, kvec, ai_work_cos, ai_work_sin)
1138 : END IF
1139 : ! projector functions: Cartesian -> spherical
1140 : na = ncoa
1141 : nb = nprjc
1142 : np = ncoc
1143 : work_cos(1:na, prjc:prjc + nb - 1) = &
1144 : MATMUL(ai_work_cos(1:na, 1:np), cprj(1:np, prjc:prjc + nb - 1))
1145 : work_sin(1:na, prjc:prjc + nb - 1) = &
1146 : MATMUL(ai_work_sin(1:na, 1:np), cprj(1:np, prjc:prjc + nb - 1))
1147 :
1148 : IF (derivative .OR. my_moment) THEN
1149 : DO idir = 1, 3
1150 : work_dcos(1:na, prjc:prjc + nb - 1, idir) = &
1151 : MATMUL(ai_work_dcos(1:na, 1:np, idir), cprj(1:np, prjc:prjc + nb - 1))
1152 : work_dsin(1:na, prjc:prjc + nb - 1, idir) = &
1153 : MATMUL(ai_work_dsin(1:na, 1:np, idir), cprj(1:np, prjc:prjc + nb - 1))
1154 : END DO
1155 : END IF
1156 :
1157 : prjc = prjc + nprjc
1158 : END DO
1159 :
1160 : ! contract gto basis set into acint
1161 : na = nsgf_seta(iset)
1162 : nb = nppnl
1163 : np = ncoa
1164 : clist%acint(sgfa:sgfa + na - 1, 1:nb, 1) = &
1165 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_cos(1:np, 1:nb))
1166 : clist_sin%acint(sgfa:sgfa + na - 1, 1:nb, 1) = &
1167 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_sin(1:np, 1:nb))
1168 : IF (derivative .OR. my_moment) THEN
1169 : DO idir = 1, 3
1170 : clist%acint(sgfa:sgfa + na - 1, 1:nb, 1 + idir) = &
1171 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_dcos(1:np, 1:nb, idir))
1172 : clist_sin%acint(sgfa:sgfa + na - 1, 1:nb, 1 + idir) = &
1173 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_dsin(1:np, 1:nb, idir))
1174 : END DO
1175 : END IF
1176 :
1177 : ! multiply with interaction matrix h_ij of the nl pp
1178 : clist%achint(sgfa:sgfa + na - 1, 1:nb, 1) = &
1179 : MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, 1), vprj_ppnl(1:nb, 1:nb))
1180 : clist_sin%achint(sgfa:sgfa + na - 1, 1:nb, 1) = &
1181 : MATMUL(clist_sin%acint(sgfa:sgfa + na - 1, 1:nb, 1), vprj_ppnl(1:nb, 1:nb))
1182 : IF (derivative .OR. my_moment) THEN
1183 : DO idir = 1, 3
1184 : clist%achint(sgfa:sgfa + na - 1, 1:nb, 1 + idir) = &
1185 : MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, 1 + idir), vprj_ppnl(1:nb, 1:nb))
1186 : clist_sin%achint(sgfa:sgfa + na - 1, 1:nb, 1 + idir) = &
1187 : MATMUL(clist_sin%acint(sgfa:sgfa + na - 1, 1:nb, 1 + idir), vprj_ppnl(1:nb, 1:nb))
1188 : END DO
1189 : END IF
1190 : END IF
1191 :
1192 : END DO
1193 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
1194 : clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1)))
1195 : clist_sin%maxac = MAXVAL(ABS(clist_sin%acint(:, :, 1)))
1196 : clist_sin%maxach = MAXVAL(ABS(clist_sin%achint(:, :, 1)))
1197 : END DO
1198 :
1199 : DEALLOCATE (work_cos, work_sin, ai_work_cos, ai_work_sin)
1200 : IF (derivative .OR. my_moment) DEALLOCATE (work_dcos, work_dsin, ai_work_dcos, ai_work_dsin)
1201 : IF (my_moment) DEALLOCATE (ai_work_mcos, ai_work_msin)
1202 :
1203 : !$OMP END PARALLEL
1204 :
1205 184 : DEALLOCATE (gpotential, spotential)
1206 :
1207 184 : CALL timestop(handle)
1208 :
1209 552 : END SUBROUTINE build_sap_exp_ints
1210 :
1211 : ! **************************************************************************************************
1212 : !> \brief Calculate the force associated to non-local pseudo potential in the velocity gauge
1213 : !> \param qs_env ...
1214 : !> \param particle_set ...
1215 : !> \date 09.2023
1216 : !> \author Guillaume Le Breton
1217 : ! **************************************************************************************************
1218 32 : SUBROUTINE velocity_gauge_nl_force(qs_env, particle_set)
1219 : TYPE(qs_environment_type), POINTER :: qs_env
1220 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1221 :
1222 : CHARACTER(len=*), PARAMETER :: routiuneN = "velocity_gauge_nl_force"
1223 :
1224 : INTEGER :: handle, i, iac, iatom, ibc, icol, idir, ikind, irow, jatom, jkind, kac, katom, &
1225 : kbc, kkind, maxl, maxlgto, maxlppnl, na, natom, nb, nkind, np, slot
1226 32 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
1227 : INTEGER, DIMENSION(3) :: cell_b
1228 : LOGICAL :: found_imag, found_real
1229 : REAL(dp) :: eps_ppnl, f0, sign_imag
1230 : REAL(KIND=dp), DIMENSION(3) :: fa, fb, rab, vec_pot
1231 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1232 32 : POINTER :: sab_orb, sap_ppnl
1233 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1234 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1235 32 : DIMENSION(:) :: basis_set
1236 : TYPE(dft_control_type), POINTER :: dft_control
1237 32 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_ao_im
1238 : TYPE(cell_type), POINTER :: cell
1239 32 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1240 : TYPE(alist_type), POINTER :: alist_cos_ac, alist_cos_bc, &
1241 : alist_sin_ac, alist_sin_bc
1242 32 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint_cos, achint_sin, acint_cos, &
1243 32 : acint_sin, bchint_cos, bchint_sin, &
1244 32 : bcint_cos, bcint_sin
1245 32 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix_p_imag, matrix_p_real
1246 64 : REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
1247 32 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1248 32 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1249 : TYPE(qs_rho_type), POINTER :: rho
1250 32 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int_cos, sap_int_sin
1251 :
1252 32 : CALL timeset(routiuneN, handle)
1253 :
1254 32 : NULLIFY (sap_ppnl)
1255 :
1256 : CALL get_qs_env(qs_env, &
1257 32 : sap_ppnl=sap_ppnl)
1258 :
1259 32 : IF (ASSOCIATED(sap_ppnl)) THEN
1260 32 : NULLIFY (qs_kind_set, cell, dft_control, force, sab_orb, atomic_kind_set, &
1261 32 : sap_int_cos, sap_int_sin)
1262 : ! Load and initialized the required quantities
1263 :
1264 : CALL get_qs_env(qs_env, &
1265 : sab_orb=sab_orb, &
1266 : force=force, &
1267 : dft_control=dft_control, &
1268 : qs_kind_set=qs_kind_set, &
1269 : cell=cell, &
1270 : atomic_kind_set=atomic_kind_set, &
1271 32 : rho=rho)
1272 :
1273 32 : nkind = SIZE(atomic_kind_set)
1274 32 : natom = SIZE(particle_set)
1275 32 : eps_ppnl = dft_control%qs_control%eps_ppnl
1276 :
1277 : CALL get_qs_kind_set(qs_kind_set, &
1278 : maxlgto=maxlgto, &
1279 32 : maxlppnl=maxlppnl)
1280 :
1281 32 : maxl = MAX(maxlppnl, maxlgto)
1282 32 : CALL init_orbital_pointers(maxl + 1)
1283 :
1284 : ! initalize sab_int types to store the integrals
1285 416 : ALLOCATE (sap_int_cos(nkind*nkind), sap_int_sin(nkind*nkind))
1286 160 : DO i = 1, SIZE(sap_int_cos)
1287 128 : NULLIFY (sap_int_cos(i)%alist, sap_int_cos(i)%asort, sap_int_cos(i)%aindex)
1288 128 : sap_int_cos(i)%nalist = 0
1289 128 : NULLIFY (sap_int_sin(i)%alist, sap_int_sin(i)%asort, sap_int_sin(i)%aindex)
1290 160 : sap_int_sin(i)%nalist = 0
1291 : END DO
1292 :
1293 : ! get basis set
1294 160 : ALLOCATE (basis_set(nkind))
1295 96 : DO ikind = 1, nkind
1296 64 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1297 96 : IF (ASSOCIATED(orb_basis_set)) THEN
1298 64 : basis_set(ikind)%gto_basis_set => orb_basis_set
1299 : ELSE
1300 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
1301 : END IF
1302 : END DO
1303 :
1304 : !get vector potential
1305 128 : vec_pot = dft_control%rtp_control%vec_pot
1306 :
1307 536 : force_thread = 0.0_dp
1308 :
1309 32 : CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao, rho_ao_im=rho_ao_im)
1310 : ! To avoid FOR loop over spin, sum the 2 spin into the first one directly. Undone later on
1311 32 : IF (SIZE(rho_ao) == 2) THEN
1312 : CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, &
1313 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1314 : CALL dbcsr_add(rho_ao_im(1)%matrix, rho_ao_im(2)%matrix, &
1315 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1316 : END IF
1317 :
1318 : ! Compute cosap = <a|cos kr|p>, sindap = <a|sin kr|p>, cosdap = <da/dRA|cos kr|p>, and sindap = <da/dRA|sin kr|p>
1319 : CALL build_sap_exp_ints(sap_int_cos, sap_int_sin, sap_ppnl, qs_kind_set, particle_set, &
1320 32 : cell, kvec=vec_pot, basis_set=basis_set, nkind=nkind, derivative=.TRUE.)
1321 32 : CALL sap_sort(sap_int_cos)
1322 32 : CALL sap_sort(sap_int_sin)
1323 :
1324 : ! Compute the force, on nuclei A it is given by: Re(P_ab) Re(dV_ab/dRA) - Im(P_ab) Im(dV_ab/dRA)
1325 :
1326 : !$OMP PARALLEL &
1327 : !$OMP DEFAULT (NONE) &
1328 : !$OMP SHARED (basis_set, sab_orb, sap_int_cos, sap_int_sin, eps_ppnl, nkind, natom,&
1329 : !$OMP rho_ao, rho_ao_im) &
1330 : !$OMP PRIVATE (matrix_p_real, matrix_p_imag, acint_cos, achint_cos, bcint_cos, bchint_cos, acint_sin,&
1331 : !$OMP achint_sin, bcint_sin, bchint_sin, slot, ikind, jkind, iatom, jatom,&
1332 : !$OMP cell_b, rab, irow, icol, fa, fb, f0, found_real, found_imag, sign_imag, &
1333 : !$OMP kkind, iac, ibc, alist_cos_ac, alist_cos_bc, alist_sin_ac, alist_sin_bc, kac, kbc,&
1334 : !$OMP na, np, nb, katom) &
1335 32 : !$OMP REDUCTION (+ : force_thread )
1336 :
1337 : NULLIFY (acint_cos, bcint_cos, achint_cos, bchint_cos)
1338 : NULLIFY (acint_sin, bcint_sin, achint_sin, bchint_sin)
1339 :
1340 : ! loop over atom pairs
1341 : !$OMP DO SCHEDULE(GUIDED)
1342 : DO slot = 1, sab_orb(1)%nl_size
1343 : ikind = sab_orb(1)%nlist_task(slot)%ikind
1344 : jkind = sab_orb(1)%nlist_task(slot)%jkind
1345 : iatom = sab_orb(1)%nlist_task(slot)%iatom
1346 : jatom = sab_orb(1)%nlist_task(slot)%jatom
1347 : cell_b(:) = sab_orb(1)%nlist_task(slot)%cell
1348 : rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
1349 :
1350 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
1351 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
1352 :
1353 : ! Use the symmetry of the first derivatives
1354 : IF (iatom == jatom) THEN
1355 : f0 = 1.0_dp
1356 : ELSE
1357 : f0 = 2.0_dp
1358 : END IF
1359 :
1360 : fa = 0.0_dp
1361 : fb = 0.0_dp
1362 :
1363 : IF (iatom <= jatom) THEN
1364 : irow = iatom
1365 : icol = jatom
1366 : sign_imag = +1.0_dp
1367 : ELSE
1368 : irow = jatom
1369 : icol = iatom
1370 : sign_imag = -1.0_dp
1371 : END IF
1372 : NULLIFY (matrix_p_real, matrix_p_imag)
1373 : CALL dbcsr_get_block_p(rho_ao(1)%matrix, irow, icol, matrix_p_real, found_real)
1374 : CALL dbcsr_get_block_p(rho_ao_im(1)%matrix, irow, icol, matrix_p_imag, found_imag)
1375 :
1376 : IF (found_real .OR. found_imag) THEN
1377 : ! loop over the <gto_a|ppln_c>h_ij<ppnl_c|gto_b> pairs
1378 : DO kkind = 1, nkind
1379 : iac = ikind + nkind*(kkind - 1)
1380 : ibc = jkind + nkind*(kkind - 1)
1381 : IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist)) CYCLE
1382 : IF (.NOT. ASSOCIATED(sap_int_cos(ibc)%alist)) CYCLE
1383 : IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist)) CYCLE
1384 : IF (.NOT. ASSOCIATED(sap_int_sin(ibc)%alist)) CYCLE
1385 : CALL get_alist(sap_int_cos(iac), alist_cos_ac, iatom)
1386 : CALL get_alist(sap_int_cos(ibc), alist_cos_bc, jatom)
1387 : CALL get_alist(sap_int_sin(iac), alist_sin_ac, iatom)
1388 : CALL get_alist(sap_int_sin(ibc), alist_sin_bc, jatom)
1389 : IF (.NOT. ASSOCIATED(alist_cos_ac)) CYCLE
1390 : IF (.NOT. ASSOCIATED(alist_cos_bc)) CYCLE
1391 : IF (.NOT. ASSOCIATED(alist_sin_ac)) CYCLE
1392 : IF (.NOT. ASSOCIATED(alist_sin_bc)) CYCLE
1393 :
1394 : ! only use cos for indexing, as cos and sin integrals are constructed by the same routine
1395 : ! in the same way
1396 : DO kac = 1, alist_cos_ac%nclist
1397 : DO kbc = 1, alist_cos_bc%nclist
1398 : ! the next two ifs should be the same for sine integrals
1399 : IF (alist_cos_ac%clist(kac)%catom /= alist_cos_bc%clist(kbc)%catom) CYCLE
1400 : IF (ALL(cell_b + alist_cos_bc%clist(kbc)%cell - alist_cos_ac%clist(kac)%cell == 0)) THEN
1401 : ! screening
1402 : IF (alist_cos_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
1403 : .AND. alist_cos_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl &
1404 : .AND. alist_sin_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
1405 : .AND. alist_sin_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
1406 :
1407 : acint_cos => alist_cos_ac%clist(kac)%acint
1408 : bcint_cos => alist_cos_bc%clist(kbc)%acint
1409 : achint_cos => alist_cos_ac%clist(kac)%achint
1410 : bchint_cos => alist_cos_bc%clist(kbc)%achint
1411 : acint_sin => alist_sin_ac%clist(kac)%acint
1412 : bcint_sin => alist_sin_bc%clist(kbc)%acint
1413 : achint_sin => alist_sin_ac%clist(kac)%achint
1414 : bchint_sin => alist_sin_bc%clist(kbc)%achint
1415 :
1416 : na = SIZE(acint_cos, 1)
1417 : np = SIZE(acint_cos, 2)
1418 : nb = SIZE(bcint_cos, 1)
1419 : ! Re(dV_ab/dRA) = <da/dRA|cos kr|p><p|cos kr|b> + <db/dRA|cos kr|p><p|cos kr|a>
1420 : ! + <da/dRA|sin kr|p><p|sin kr|b> + <db/dRA|sin kr|p><p|sin|a>
1421 : ! Im(dV_ab/dRA) = <da/dRA|sin kr|p><p|cos kr|b> - <db/dRA|sin kr|p><p|cos kr|a>
1422 : ! - <da/dRA|cos kr|p><p|sin kr|b> + <db/dRA|cos kr|p><p|sin|a>
1423 : katom = alist_cos_ac%clist(kac)%catom
1424 : DO idir = 1, 3
1425 : IF (iatom <= jatom) THEN
1426 : ! For fa:
1427 : IF (found_real) THEN
1428 : fa(idir) = SUM(matrix_p_real(1:na, 1:nb)* &
1429 : (+MATMUL(acint_cos(1:na, 1:np, 1 + idir), TRANSPOSE(bchint_cos(1:nb, 1:np, 1))) &
1430 : + MATMUL(acint_sin(1:na, 1:np, 1 + idir), TRANSPOSE(bchint_sin(1:nb, 1:np, 1)))))
1431 : END IF
1432 : IF (found_imag) THEN
1433 : fa(idir) = fa(idir) - sign_imag*SUM(matrix_p_imag(1:na, 1:nb)* &
1434 : (+MATMUL(acint_sin(1:na, 1:np, 1 + idir), TRANSPOSE(bchint_cos(1:nb, 1:np, 1))) &
1435 : - MATMUL(acint_cos(1:na, 1:np, 1 + idir), TRANSPOSE(bchint_sin(1:nb, 1:np, 1)))))
1436 : END IF
1437 : ! For fb:
1438 : IF (found_real) THEN
1439 : fb(idir) = SUM(matrix_p_real(1:na, 1:nb)* &
1440 : (+MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, 1 + idir))) &
1441 : + MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, 1 + idir)))))
1442 : END IF
1443 : IF (found_imag) THEN
1444 : fb(idir) = fb(idir) - sign_imag*SUM(matrix_p_imag(1:na, 1:nb)* &
1445 : (-MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, 1 + idir))) &
1446 : + MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, 1 + idir)))))
1447 : END IF
1448 : ELSE
1449 : ! For fa:
1450 : IF (found_real) THEN
1451 : fa(idir) = SUM(matrix_p_real(1:nb, 1:na)* &
1452 : (+MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, 1 + idir))) &
1453 : + MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, 1 + idir)))))
1454 : END IF
1455 : IF (found_imag) THEN
1456 : fa(idir) = fa(idir) - sign_imag*SUM(matrix_p_imag(1:nb, 1:na)* &
1457 : (+MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, 1 + idir))) &
1458 : - MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, 1 + idir)))))
1459 : END IF
1460 : ! For fb
1461 : IF (found_real) THEN
1462 : fb(idir) = SUM(matrix_p_real(1:nb, 1:na)* &
1463 : (+MATMUL(bcint_cos(1:nb, 1:np, 1 + idir), TRANSPOSE(achint_cos(1:na, 1:np, 1))) &
1464 : + MATMUL(bcint_sin(1:nb, 1:np, 1 + idir), TRANSPOSE(achint_sin(1:na, 1:np, 1)))))
1465 : END IF
1466 : IF (found_imag) THEN
1467 : fb(idir) = fb(idir) - sign_imag*SUM(matrix_p_imag(1:nb, 1:na)* &
1468 : (-MATMUL(bcint_cos(1:nb, 1:np, 1 + idir), TRANSPOSE(achint_sin(1:na, 1:np, 1))) &
1469 : + MATMUL(bcint_sin(1:nb, 1:np, 1 + idir), TRANSPOSE(achint_cos(1:na, 1:np, 1)))))
1470 : END IF
1471 : END IF
1472 : force_thread(idir, iatom) = force_thread(idir, iatom) + f0*fa(idir)
1473 : force_thread(idir, katom) = force_thread(idir, katom) - f0*fa(idir)
1474 : force_thread(idir, jatom) = force_thread(idir, jatom) + f0*fb(idir)
1475 : force_thread(idir, katom) = force_thread(idir, katom) - f0*fb(idir)
1476 : END DO
1477 : EXIT
1478 : END IF
1479 : END DO
1480 : END DO
1481 : END DO
1482 : END IF
1483 :
1484 : END DO
1485 :
1486 : !$OMP END PARALLEL
1487 :
1488 : ! Update the force
1489 32 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
1490 : !$OMP DO
1491 : DO iatom = 1, natom
1492 126 : i = atom_of_kind(iatom)
1493 126 : ikind = kind_of(iatom)
1494 504 : force(ikind)%gth_ppnl(:, i) = force(ikind)%gth_ppnl(:, i) + force_thread(:, iatom)
1495 : END DO
1496 : !$OMP END DO
1497 :
1498 : ! Clean up
1499 32 : IF (SIZE(rho_ao) == 2) THEN
1500 : CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, &
1501 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1502 : CALL dbcsr_add(rho_ao_im(1)%matrix, rho_ao_im(2)%matrix, &
1503 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1504 : END IF
1505 32 : CALL release_sap_int(sap_int_cos)
1506 32 : CALL release_sap_int(sap_int_sin)
1507 :
1508 64 : DEALLOCATE (basis_set, atom_of_kind, kind_of)
1509 :
1510 : END IF
1511 :
1512 32 : CALL timestop(handle)
1513 :
1514 64 : END SUBROUTINE velocity_gauge_nl_force
1515 :
1516 : END MODULE rt_propagation_velocity_gauge
|
