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 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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
25 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
26 : dbcsr_deallocate_matrix_set
27 : USE dbcsr_api, ONLY: dbcsr_add,&
28 : dbcsr_create,&
29 : dbcsr_get_block_p,&
30 : dbcsr_init_p,&
31 : dbcsr_p_type,&
32 : dbcsr_set,&
33 : dbcsr_type_antisymmetric,&
34 : dbcsr_type_symmetric
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: init_orbital_pointers,&
48 : nco,&
49 : ncoset
50 : USE particle_types, ONLY: particle_type
51 : USE qs_environment_types, ONLY: get_qs_env,&
52 : qs_environment_type
53 : USE qs_force_types, ONLY: qs_force_type
54 : USE qs_kind_types, ONLY: get_qs_kind,&
55 : get_qs_kind_set,&
56 : qs_kind_type
57 : USE qs_ks_types, ONLY: get_ks_env,&
58 : qs_ks_env_type
59 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
60 : USE qs_operators_ao, ONLY: build_lin_mom_matrix
61 : USE qs_rho_types, ONLY: qs_rho_get,&
62 : qs_rho_type
63 : USE sap_kind_types, ONLY: alist_type,&
64 : clist_type,&
65 : get_alist,&
66 : release_sap_int,&
67 : sap_int_type,&
68 : sap_sort
69 : USE virial_types, ONLY: virial_type
70 :
71 : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
72 : !$ omp_init_lock, omp_set_lock, &
73 : !$ omp_unset_lock, omp_destroy_lock
74 :
75 : #include "./base/base_uses.f90"
76 :
77 : IMPLICIT NONE
78 :
79 : PRIVATE
80 :
81 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_velocity_gauge'
82 :
83 : PUBLIC :: velocity_gauge_ks_matrix, update_vector_potential, velocity_gauge_nl_force
84 :
85 : CONTAINS
86 :
87 : ! **************************************************************************************************
88 : !> \brief ...
89 : !> \param qs_env ...
90 : !> \param subtract_nl_term ...
91 : ! **************************************************************************************************
92 62 : SUBROUTINE velocity_gauge_ks_matrix(qs_env, subtract_nl_term)
93 : TYPE(qs_environment_type), POINTER :: qs_env
94 : LOGICAL, INTENT(IN), OPTIONAL :: subtract_nl_term
95 :
96 : CHARACTER(len=*), PARAMETER :: routineN = 'velocity_gauge_ks_matrix'
97 :
98 : INTEGER :: handle, idir, image, nder, nimages
99 62 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
100 : LOGICAL :: calculate_forces, my_subtract_nl_term, &
101 : ppnl_present, use_virial
102 : REAL(KIND=dp) :: eps_ppnl, factor
103 : REAL(KIND=dp), DIMENSION(3) :: vec_pot
104 62 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
105 : TYPE(cell_type), POINTER :: cell
106 62 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: momentum, nl_term
107 62 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im, matrix_nl, &
108 62 : matrix_p, matrix_s
109 : TYPE(dft_control_type), POINTER :: dft_control
110 : TYPE(kpoint_type), POINTER :: kpoints
111 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
112 62 : POINTER :: sab_orb, sap_ppnl
113 62 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
114 62 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
115 62 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
116 : TYPE(qs_ks_env_type), POINTER :: ks_env
117 : TYPE(qs_rho_type), POINTER :: rho
118 : TYPE(section_vals_type), POINTER :: input
119 : TYPE(virial_type), POINTER :: virial
120 :
121 62 : CALL timeset(routineN, handle)
122 :
123 62 : CALL cite_reference(Mattiat2022)
124 :
125 62 : my_subtract_nl_term = .FALSE.
126 62 : IF (PRESENT(subtract_nl_term)) my_subtract_nl_term = subtract_nl_term
127 :
128 62 : NULLIFY (dft_control, matrix_s, sab_orb, matrix_h, cell, input, matrix_h_im, kpoints, cell_to_index, &
129 62 : sap_ppnl, particle_set, qs_kind_set, atomic_kind_set, virial, force, matrix_p, rho, matrix_nl)
130 :
131 : CALL get_qs_env(qs_env, &
132 : rho=rho, &
133 : dft_control=dft_control, &
134 : sab_orb=sab_orb, &
135 : sap_ppnl=sap_ppnl, &
136 : matrix_s_kp=matrix_s, &
137 : matrix_h_kp=matrix_h, &
138 : cell=cell, &
139 : input=input, &
140 62 : matrix_h_im_kp=matrix_h_im)
141 :
142 62 : nimages = dft_control%nimages
143 62 : ppnl_present = ASSOCIATED(sap_ppnl)
144 :
145 62 : IF (nimages > 1) THEN
146 0 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
147 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
148 : END IF
149 :
150 62 : IF (my_subtract_nl_term) THEN
151 8 : IF (ppnl_present) THEN
152 : CALL get_qs_env(qs_env, &
153 : qs_kind_set=qs_kind_set, &
154 : particle_set=particle_set, &
155 : atomic_kind_set=atomic_kind_set, &
156 : virial=virial, &
157 : rho=rho, &
158 4 : force=force)
159 :
160 4 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
161 4 : calculate_forces = .FALSE.
162 4 : use_virial = .FALSE.
163 4 : nder = 1
164 4 : eps_ppnl = dft_control%qs_control%eps_ppnl
165 :
166 4 : CALL dbcsr_allocate_matrix_set(matrix_nl, 1, nimages)
167 8 : DO image = 1, nimages
168 4 : ALLOCATE (matrix_nl(1, image)%matrix)
169 4 : CALL dbcsr_create(matrix_nl(1, image)%matrix, template=matrix_s(1, 1)%matrix)
170 4 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_nl(1, image)%matrix, sab_orb)
171 8 : CALL dbcsr_set(matrix_nl(1, image)%matrix, zero)
172 : END DO
173 :
174 : CALL build_core_ppnl(matrix_nl, matrix_p, force, virial, calculate_forces, use_virial, nder, &
175 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
176 4 : nimages, cell_to_index, "ORB")
177 :
178 8 : DO image = 1, nimages
179 8 : CALL dbcsr_add(matrix_h(1, image)%matrix, matrix_nl(1, image)%matrix, one, -one)
180 : END DO
181 :
182 4 : CALL dbcsr_deallocate_matrix_set(matrix_nl)
183 : END IF
184 : END IF
185 :
186 : !get vector potential
187 248 : vec_pot = dft_control%rtp_control%vec_pot
188 :
189 : ! allocate and build matrices for linear momentum term
190 62 : NULLIFY (momentum)
191 62 : CALL dbcsr_allocate_matrix_set(momentum, 3)
192 248 : DO idir = 1, 3
193 186 : CALL dbcsr_init_p(momentum(idir)%matrix)
194 : CALL dbcsr_create(momentum(idir)%matrix, template=matrix_s(1, 1)%matrix, &
195 186 : matrix_type=dbcsr_type_antisymmetric)
196 186 : CALL cp_dbcsr_alloc_block_from_nbl(momentum(idir)%matrix, sab_orb)
197 248 : CALL dbcsr_set(momentum(idir)%matrix, zero)
198 : END DO
199 62 : CALL build_lin_mom_matrix(qs_env, momentum)
200 :
201 : ! set imaginary part of KS matrix to zero
202 124 : DO image = 1, nimages
203 124 : CALL dbcsr_set(matrix_h_im(1, image)%matrix, zero)
204 : END DO
205 :
206 : ! add linear term in vector potential to imaginary part of KS-matrix
207 124 : DO image = 1, nimages
208 310 : DO idir = 1, 3
209 248 : CALL dbcsr_add(matrix_h_im(1, image)%matrix, momentum(idir)%matrix, one, -vec_pot(idir))
210 : END DO
211 : END DO
212 :
213 62 : CALL dbcsr_deallocate_matrix_set(momentum)
214 :
215 : ! add quadratic term to real part of KS matrix
216 62 : factor = 0._dp
217 248 : DO idir = 1, 3
218 248 : factor = factor + vec_pot(idir)**2
219 : END DO
220 :
221 124 : DO image = 1, nimages
222 124 : CALL dbcsr_add(matrix_h(1, image)%matrix, matrix_s(1, image)%matrix, one, 0.5*factor)
223 : END DO
224 :
225 : ! add Non local term
226 62 : IF (ppnl_present) THEN
227 40 : IF (dft_control%rtp_control%nl_gauge_transform) THEN
228 40 : NULLIFY (nl_term)
229 40 : CALL dbcsr_allocate_matrix_set(nl_term, 2)
230 :
231 40 : CALL dbcsr_init_p(nl_term(1)%matrix)
232 : CALL dbcsr_create(nl_term(1)%matrix, template=matrix_s(1, 1)%matrix, &
233 40 : matrix_type=dbcsr_type_symmetric, name="nl gauge term real part")
234 40 : CALL cp_dbcsr_alloc_block_from_nbl(nl_term(1)%matrix, sab_orb)
235 40 : CALL dbcsr_set(nl_term(1)%matrix, zero)
236 :
237 40 : CALL dbcsr_init_p(nl_term(2)%matrix)
238 : CALL dbcsr_create(nl_term(2)%matrix, template=matrix_s(1, 1)%matrix, &
239 40 : matrix_type=dbcsr_type_antisymmetric, name="nl gauge term imaginary part")
240 40 : CALL cp_dbcsr_alloc_block_from_nbl(nl_term(2)%matrix, sab_orb)
241 40 : CALL dbcsr_set(nl_term(2)%matrix, zero)
242 :
243 40 : CALL velocity_gauge_nl_term(qs_env, nl_term, vec_pot)
244 :
245 80 : DO image = 1, nimages
246 40 : CALL dbcsr_add(matrix_h(1, image)%matrix, nl_term(1)%matrix, one, one)
247 80 : CALL dbcsr_add(matrix_h_im(1, image)%matrix, nl_term(2)%matrix, one, one)
248 : END DO
249 40 : CALL dbcsr_deallocate_matrix_set(nl_term)
250 : END IF
251 : END IF
252 :
253 62 : CALL timestop(handle)
254 :
255 62 : END SUBROUTINE velocity_gauge_ks_matrix
256 :
257 : ! **************************************************************************************************
258 : !> \brief Update the vector potential in the case where a time-dependant
259 : !> electric field is apply.
260 : !> \param qs_env ...
261 : !> \param dft_control ...
262 : ! **************************************************************************************************
263 32 : SUBROUTINE update_vector_potential(qs_env, dft_control)
264 : TYPE(qs_environment_type), INTENT(INOUT), POINTER :: qs_env
265 : TYPE(dft_control_type), INTENT(INOUT), POINTER :: dft_control
266 :
267 : REAL(kind=dp) :: field(3)
268 :
269 32 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
270 128 : dft_control%rtp_control%field = field
271 128 : dft_control%rtp_control%vec_pot = dft_control%rtp_control%vec_pot - field*qs_env%rtp%dt
272 : ! Update the vec_pot_initial value for RTP restart:
273 128 : dft_control%efield_fields(1)%efield%vec_pot_initial = dft_control%rtp_control%vec_pot
274 :
275 32 : END SUBROUTINE update_vector_potential
276 :
277 : ! **************************************************************************************************
278 : !> \brief ...
279 : !> \param qs_env ...
280 : !> \param nl_term ...
281 : !> \param vec_pot ...
282 : ! **************************************************************************************************
283 40 : SUBROUTINE velocity_gauge_nl_term(qs_env, nl_term, vec_pot)
284 : TYPE(qs_environment_type), POINTER :: qs_env
285 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
286 : POINTER :: nl_term
287 : REAL(KIND=dp), DIMENSION(3), INTENT(in) :: vec_pot
288 :
289 : CHARACTER(len=*), PARAMETER :: routiuneN = "velocity_gauge_nl_term"
290 :
291 : INTEGER :: handle, i, iac, iatom, ibc, icol, ikind, &
292 : irow, jatom, jkind, kac, kbc, kkind, &
293 : maxl, maxlgto, maxlppnl, na, natom, &
294 : nb, nkind, np, slot
295 : INTEGER, DIMENSION(3) :: cell_b
296 : LOGICAL :: found
297 : REAL(dp) :: eps_ppnl
298 : REAL(KIND=dp), DIMENSION(3) :: rab
299 40 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: imag_block, real_block
300 40 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint_cos, achint_sin, acint_cos, &
301 40 : acint_sin, bchint_cos, bchint_sin, &
302 40 : bcint_cos, bcint_sin
303 : TYPE(alist_type), POINTER :: alist_cos_ac, alist_cos_bc, &
304 : alist_sin_ac, alist_sin_bc
305 40 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
306 : TYPE(cell_type), POINTER :: cell
307 : TYPE(dft_control_type), POINTER :: dft_control
308 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
309 40 : DIMENSION(:) :: basis_set
310 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
311 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
312 40 : POINTER :: sab_orb, sap_ppnl
313 40 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
314 40 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
315 40 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int_cos, sap_int_sin
316 :
317 : !$ INTEGER(kind=omp_lock_kind), &
318 40 : !$ ALLOCATABLE, DIMENSION(:) :: locks
319 : !$ INTEGER(KIND=int_8) :: iatom8
320 : !$ INTEGER :: lock_num, hash
321 : !$ INTEGER, PARAMETER :: nlock = 501
322 :
323 : MARK_USED(int_8)
324 :
325 40 : CALL timeset(routiuneN, handle)
326 :
327 40 : NULLIFY (sap_ppnl, sab_orb)
328 : CALL get_qs_env(qs_env, &
329 : sap_ppnl=sap_ppnl, &
330 40 : sab_orb=sab_orb)
331 :
332 40 : IF (ASSOCIATED(sap_ppnl)) THEN
333 40 : NULLIFY (qs_kind_set, particle_set, cell, dft_control)
334 : CALL get_qs_env(qs_env, &
335 : dft_control=dft_control, &
336 : qs_kind_set=qs_kind_set, &
337 : particle_set=particle_set, &
338 : cell=cell, &
339 40 : atomic_kind_set=atomic_kind_set)
340 :
341 40 : nkind = SIZE(atomic_kind_set)
342 40 : natom = SIZE(particle_set)
343 40 : eps_ppnl = dft_control%qs_control%eps_ppnl
344 :
345 : CALL get_qs_kind_set(qs_kind_set, &
346 : maxlgto=maxlgto, &
347 40 : maxlppnl=maxlppnl)
348 :
349 40 : maxl = MAX(maxlppnl, maxlgto)
350 40 : CALL init_orbital_pointers(maxl + 1)
351 :
352 : ! initalize sab_int types to store the integrals
353 40 : NULLIFY (sap_int_cos, sap_int_sin)
354 520 : ALLOCATE (sap_int_cos(nkind*nkind), sap_int_sin(nkind*nkind))
355 200 : DO i = 1, SIZE(sap_int_cos)
356 160 : NULLIFY (sap_int_cos(i)%alist, sap_int_cos(i)%asort, sap_int_cos(i)%aindex)
357 160 : sap_int_cos(i)%nalist = 0
358 160 : NULLIFY (sap_int_sin(i)%alist, sap_int_sin(i)%asort, sap_int_sin(i)%aindex)
359 200 : sap_int_sin(i)%nalist = 0
360 : END DO
361 :
362 : ! get basis set
363 200 : ALLOCATE (basis_set(nkind))
364 120 : DO ikind = 1, nkind
365 80 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
366 120 : IF (ASSOCIATED(orb_basis_set)) THEN
367 80 : basis_set(ikind)%gto_basis_set => orb_basis_set
368 : ELSE
369 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
370 : END IF
371 : END DO
372 :
373 : ! calculate exponential integrals
374 : CALL build_sap_exp_ints(sap_int_cos, sap_int_sin, sap_ppnl, qs_kind_set, particle_set, &
375 : cell, kvec=vec_pot, basis_set=basis_set, nkind=nkind, &
376 40 : derivative=.FALSE.)
377 :
378 40 : CALL sap_sort(sap_int_cos)
379 40 : CALL sap_sort(sap_int_sin)
380 :
381 : ! assemble the integrals for the gauge term
382 : !$OMP PARALLEL &
383 : !$OMP DEFAULT (NONE) &
384 : !$OMP SHARED (basis_set, nl_term, sab_orb, sap_int_cos, sap_int_sin, eps_ppnl, locks, nkind, natom) &
385 : !$OMP PRIVATE (real_block, imag_block, acint_cos, achint_cos, bcint_cos, bchint_cos, acint_sin,&
386 : !$OMP achint_sin, bcint_sin, bchint_sin, slot, ikind, jkind, iatom, jatom, cell_b, rab, irow, icol,&
387 : !$OMP found, kkind, iac, ibc, alist_cos_ac, alist_cos_bc, alist_sin_ac, alist_sin_bc, kac, kbc, &
388 40 : !$OMP na, np, nb, iatom8, hash)
389 :
390 : !$OMP SINGLE
391 : !$ ALLOCATE (locks(nlock))
392 : !$OMP END SINGLE
393 :
394 : !$OMP DO
395 : !$ DO lock_num = 1, nlock
396 : !$ call omp_init_lock(locks(lock_num))
397 : !$ END DO
398 : !$OMP END DO
399 :
400 : NULLIFY (real_block, imag_block)
401 : NULLIFY (acint_cos, bcint_cos, achint_cos, bchint_cos)
402 : NULLIFY (acint_sin, bcint_sin, achint_sin, bchint_sin)
403 :
404 : ! loop over atom pairs
405 : !$OMP DO SCHEDULE(GUIDED)
406 : DO slot = 1, sab_orb(1)%nl_size
407 : ikind = sab_orb(1)%nlist_task(slot)%ikind
408 : jkind = sab_orb(1)%nlist_task(slot)%jkind
409 : iatom = sab_orb(1)%nlist_task(slot)%iatom
410 : jatom = sab_orb(1)%nlist_task(slot)%jatom
411 : cell_b(:) = sab_orb(1)%nlist_task(slot)%cell
412 : rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
413 :
414 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
415 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
416 :
417 : IF (iatom <= jatom) THEN
418 : irow = iatom
419 : icol = jatom
420 : ELSE
421 : irow = jatom
422 : icol = iatom
423 : END IF
424 :
425 : CALL dbcsr_get_block_p(nl_term(1)%matrix, irow, icol, real_block, found)
426 : CALL dbcsr_get_block_p(nl_term(2)%matrix, irow, icol, imag_block, found)
427 :
428 : IF (ASSOCIATED(real_block) .AND. ASSOCIATED(imag_block)) THEN
429 : ! loop over the <gto_a|ppln_c>h_ij<ppnl_c|gto_b> pairs
430 : DO kkind = 1, nkind
431 : iac = ikind + nkind*(kkind - 1)
432 : ibc = jkind + nkind*(kkind - 1)
433 : IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist)) CYCLE
434 : IF (.NOT. ASSOCIATED(sap_int_cos(ibc)%alist)) CYCLE
435 : IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist)) CYCLE
436 : IF (.NOT. ASSOCIATED(sap_int_sin(ibc)%alist)) CYCLE
437 : CALL get_alist(sap_int_cos(iac), alist_cos_ac, iatom)
438 : CALL get_alist(sap_int_cos(ibc), alist_cos_bc, jatom)
439 : CALL get_alist(sap_int_sin(iac), alist_sin_ac, iatom)
440 : CALL get_alist(sap_int_sin(ibc), alist_sin_bc, jatom)
441 : IF (.NOT. ASSOCIATED(alist_cos_ac)) CYCLE
442 : IF (.NOT. ASSOCIATED(alist_cos_bc)) CYCLE
443 : IF (.NOT. ASSOCIATED(alist_sin_ac)) CYCLE
444 : IF (.NOT. ASSOCIATED(alist_sin_bc)) CYCLE
445 :
446 : ! only use cos for indexing, as cos and sin integrals are constructed by the same routine
447 : ! in the same way
448 : DO kac = 1, alist_cos_ac%nclist
449 : DO kbc = 1, alist_cos_bc%nclist
450 : ! the next two ifs should be the same for sine integrals
451 : IF (alist_cos_ac%clist(kac)%catom /= alist_cos_bc%clist(kbc)%catom) CYCLE
452 : IF (ALL(cell_b + alist_cos_bc%clist(kbc)%cell - alist_cos_ac%clist(kac)%cell == 0)) THEN
453 : ! screening
454 : IF (alist_cos_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
455 : .AND. alist_cos_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl &
456 : .AND. alist_sin_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
457 : .AND. alist_sin_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
458 :
459 : acint_cos => alist_cos_ac%clist(kac)%acint
460 : bcint_cos => alist_cos_bc%clist(kbc)%acint
461 : achint_cos => alist_cos_ac%clist(kac)%achint
462 : bchint_cos => alist_cos_bc%clist(kbc)%achint
463 : acint_sin => alist_sin_ac%clist(kac)%acint
464 : bcint_sin => alist_sin_bc%clist(kbc)%acint
465 : achint_sin => alist_sin_ac%clist(kac)%achint
466 : bchint_sin => alist_sin_bc%clist(kbc)%achint
467 :
468 : na = SIZE(acint_cos, 1)
469 : np = SIZE(acint_cos, 2)
470 : nb = SIZE(bcint_cos, 1)
471 : !$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
472 : !$ hash = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
473 : !$ CALL omp_set_lock(locks(hash))
474 : IF (iatom <= jatom) THEN
475 : ! cos*cos + sin*sin
476 : real_block(1:na, 1:nb) = real_block(1:na, 1:nb) + &
477 : MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, 1))) + &
478 : MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, 1)))
479 : ! sin * cos - cos * sin
480 : imag_block(1:na, 1:nb) = imag_block(1:na, 1:nb) - &
481 : MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, 1))) + &
482 : MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, 1)))
483 : ELSE
484 : ! cos*cos + sin*sin
485 : real_block(1:nb, 1:na) = real_block(1:nb, 1:na) + &
486 : MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, 1))) + &
487 : MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, 1)))
488 : ! sin * cos - cos * sin
489 : imag_block(1:nb, 1:na) = imag_block(1:nb, 1:na) - &
490 : MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, 1))) + &
491 : MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, 1)))
492 :
493 : END IF
494 : !$ CALL omp_unset_lock(locks(hash))
495 : EXIT
496 : END IF
497 : END DO
498 : END DO
499 : END DO
500 : END IF
501 :
502 : END DO
503 :
504 : !$OMP DO
505 : !$ DO lock_num = 1, nlock
506 : !$ call omp_destroy_lock(locks(lock_num))
507 : !$ END DO
508 : !$OMP END DO
509 :
510 : !$OMP SINGLE
511 : !$ DEALLOCATE (locks)
512 : !$OMP END SINGLE NOWAIT
513 :
514 : !$OMP END PARALLEL
515 40 : CALL release_sap_int(sap_int_cos)
516 40 : CALL release_sap_int(sap_int_sin)
517 :
518 80 : DEALLOCATE (basis_set)
519 : END IF
520 :
521 40 : CALL timestop(handle)
522 :
523 80 : END SUBROUTINE velocity_gauge_nl_term
524 :
525 : ! **************************************************************************************************
526 : !> \brief calculate <a|sin/cos|p> integrals and store in sap_int_type
527 : !> adapted from build_sap_ints
528 : !> Do this on each MPI task as the integrals need to be available globally.
529 : !> Might be faster than communicating as the integrals are obtained analytically.
530 : !> If asked, compute <da/dRa|sin/cos|p>
531 : !> \param sap_int_cos ...
532 : !> \param sap_int_sin ...
533 : !> \param sap_ppnl ...
534 : !> \param qs_kind_set ...
535 : !> \param particle_set ...
536 : !> \param cell ...
537 : !> \param kvec ...
538 : !> \param basis_set ...
539 : !> \param nkind ...
540 : !> \param derivative ...
541 : ! **************************************************************************************************
542 62 : SUBROUTINE build_sap_exp_ints(sap_int_cos, sap_int_sin, sap_ppnl, qs_kind_set, particle_set, cell, &
543 62 : kvec, basis_set, nkind, derivative)
544 : TYPE(sap_int_type), DIMENSION(:), INTENT(INOUT), &
545 : POINTER :: sap_int_cos, sap_int_sin
546 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
547 : INTENT(IN), POINTER :: sap_ppnl
548 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
549 : POINTER :: qs_kind_set
550 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
551 : POINTER :: particle_set
552 : TYPE(cell_type), INTENT(IN), POINTER :: cell
553 : REAL(KIND=dp), DIMENSION(3), INTENT(in) :: kvec
554 : TYPE(gto_basis_set_p_type), DIMENSION(:), &
555 : INTENT(IN) :: basis_set
556 : INTEGER, INTENT(IN) :: nkind
557 : LOGICAL, INTENT(IN) :: derivative
558 :
559 : CHARACTER(len=*), PARAMETER :: routiuneN = "build_sap_exp_ints"
560 :
561 : INTEGER :: handle, i, iac, iatom, idir, ikind, ilist, iset, jneighbor, katom, kkind, l, &
562 : lc_max, lc_min, ldai, ldints, lppnl, maxco, maxl, maxlgto, maxlppnl, maxppnl, maxsgf, na, &
563 : nb, ncoa, ncoc, nlist, nneighbor, np, nppnl, nprjc, nseta, nsgfa, prjc, sgfa, slot
564 : INTEGER, DIMENSION(3) :: cell_c
565 62 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
566 62 : nsgf_seta
567 62 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
568 : LOGICAL :: dogth
569 : REAL(KIND=dp) :: dac, ppnl_radius
570 62 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ai_work_cos, ai_work_sin, work_cos, &
571 62 : work_sin
572 62 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work_dcos, ai_work_dsin, work_dcos, &
573 62 : work_dsin
574 : REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc
575 : REAL(KIND=dp), DIMENSION(3) :: ra, rac, raf, rc, rcf
576 62 : REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_ppnl, set_radius_a
577 62 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, rpgfa, sphi_a, vprj_ppnl, zeta
578 : TYPE(clist_type), POINTER :: clist, clist_sin
579 62 : TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
580 : TYPE(gth_potential_type), POINTER :: gth_potential
581 62 : TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential
582 : TYPE(sgp_potential_type), POINTER :: sgp_potential
583 :
584 62 : CALL timeset(routiuneN, handle)
585 :
586 : CALL get_qs_kind_set(qs_kind_set, &
587 : maxco=maxco, &
588 : maxlppnl=maxlppnl, &
589 : maxppnl=maxppnl, &
590 : maxsgf=maxsgf, &
591 62 : maxlgto=maxlgto)
592 :
593 : ! maximum dimensions for allocations
594 62 : maxl = MAX(maxlppnl, maxlgto)
595 62 : ldints = MAX(maxco, ncoset(maxlppnl), maxsgf, maxppnl)
596 62 : ldai = ncoset(maxl + 1)
597 :
598 : !set up direct access to basis and potential
599 62 : NULLIFY (gpotential, spotential)
600 620 : ALLOCATE (gpotential(nkind), spotential(nkind))
601 186 : DO ikind = 1, nkind
602 124 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
603 124 : NULLIFY (gpotential(ikind)%gth_potential)
604 124 : NULLIFY (spotential(ikind)%sgp_potential)
605 186 : IF (ASSOCIATED(gth_potential)) THEN
606 124 : gpotential(ikind)%gth_potential => gth_potential
607 0 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
608 0 : spotential(ikind)%sgp_potential => sgp_potential
609 : END IF
610 : END DO
611 :
612 : !allocate sap int
613 62 : NULLIFY (clist)
614 248 : DO slot = 1, sap_ppnl(1)%nl_size
615 :
616 186 : ikind = sap_ppnl(1)%nlist_task(slot)%ikind
617 186 : kkind = sap_ppnl(1)%nlist_task(slot)%jkind
618 186 : iatom = sap_ppnl(1)%nlist_task(slot)%iatom
619 186 : katom = sap_ppnl(1)%nlist_task(slot)%jatom
620 186 : nlist = sap_ppnl(1)%nlist_task(slot)%nlist
621 186 : ilist = sap_ppnl(1)%nlist_task(slot)%ilist
622 186 : nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
623 :
624 186 : iac = ikind + nkind*(kkind - 1)
625 186 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
626 186 : IF (.NOT. ASSOCIATED(gpotential(kkind)%gth_potential) .AND. &
627 : .NOT. ASSOCIATED(spotential(kkind)%sgp_potential)) CYCLE
628 186 : IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist)) THEN
629 124 : sap_int_cos(iac)%a_kind = ikind
630 124 : sap_int_cos(iac)%p_kind = kkind
631 124 : sap_int_cos(iac)%nalist = nlist
632 558 : ALLOCATE (sap_int_cos(iac)%alist(nlist))
633 310 : DO i = 1, nlist
634 186 : NULLIFY (sap_int_cos(iac)%alist(i)%clist)
635 186 : sap_int_cos(iac)%alist(i)%aatom = 0
636 310 : sap_int_cos(iac)%alist(i)%nclist = 0
637 : END DO
638 : END IF
639 186 : IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist(ilist)%clist)) THEN
640 186 : sap_int_cos(iac)%alist(ilist)%aatom = iatom
641 186 : sap_int_cos(iac)%alist(ilist)%nclist = nneighbor
642 2046 : ALLOCATE (sap_int_cos(iac)%alist(ilist)%clist(nneighbor))
643 372 : DO i = 1, nneighbor
644 186 : clist => sap_int_cos(iac)%alist(ilist)%clist(i)
645 186 : clist%catom = 0
646 186 : NULLIFY (clist%acint)
647 186 : NULLIFY (clist%achint)
648 372 : NULLIFY (clist%sgf_list)
649 : END DO
650 : END IF
651 186 : IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist)) THEN
652 124 : sap_int_sin(iac)%a_kind = ikind
653 124 : sap_int_sin(iac)%p_kind = kkind
654 124 : sap_int_sin(iac)%nalist = nlist
655 558 : ALLOCATE (sap_int_sin(iac)%alist(nlist))
656 310 : DO i = 1, nlist
657 186 : NULLIFY (sap_int_sin(iac)%alist(i)%clist)
658 186 : sap_int_sin(iac)%alist(i)%aatom = 0
659 310 : sap_int_sin(iac)%alist(i)%nclist = 0
660 : END DO
661 : END IF
662 248 : IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist(ilist)%clist)) THEN
663 186 : sap_int_sin(iac)%alist(ilist)%aatom = iatom
664 186 : sap_int_sin(iac)%alist(ilist)%nclist = nneighbor
665 2046 : ALLOCATE (sap_int_sin(iac)%alist(ilist)%clist(nneighbor))
666 372 : DO i = 1, nneighbor
667 186 : clist => sap_int_sin(iac)%alist(ilist)%clist(i)
668 186 : clist%catom = 0
669 186 : NULLIFY (clist%acint)
670 186 : NULLIFY (clist%achint)
671 372 : NULLIFY (clist%sgf_list)
672 : END DO
673 : END IF
674 : END DO
675 :
676 : ! actual calculation of the integrals <a|cos|p> and <a|sin|p>
677 : ! allocate temporary storage using maximum dimensions
678 :
679 : !$OMP PARALLEL &
680 : !$OMP DEFAULT (NONE) &
681 : !$OMP SHARED (basis_set, gpotential, ncoset, sap_ppnl, sap_int_cos, sap_int_sin, nkind, &
682 : !$OMP ldints, maxco, nco, cell, particle_set, kvec, derivative) &
683 : !$OMP PRIVATE (slot, ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, &
684 : !$OMP cell_c, rac, dac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta,&
685 : !$OMP rpgfa, set_radius_a, sphi_a, zeta, alpha_ppnl, cprj, lppnl, nppnl, nprj_ppnl,&
686 : !$OMP ppnl_radius, vprj_ppnl, clist, clist_sin, ra, rc, ncoa, sgfa, prjc, work_cos, work_sin,&
687 : !$OMP nprjc, rprjc, lc_max, lc_min, zetc, ncoc, ai_work_sin, ai_work_cos, na, nb, np, dogth, &
688 62 : !$OMP raf, rcf, work_dcos, work_dsin, ai_work_dcos, ai_work_dsin, idir)
689 :
690 : ALLOCATE (work_cos(ldints, ldints), work_sin(ldints, ldints))
691 : ALLOCATE (ai_work_cos(maxco, maxco), ai_work_sin(maxco, maxco))
692 : IF (derivative) THEN
693 : ALLOCATE (work_dcos(ldints, ldints, 3), work_dsin(ldints, ldints, 3))
694 : ALLOCATE (ai_work_dcos(maxco, maxco, 3), ai_work_dsin(maxco, maxco, 3))
695 : END IF
696 : work_cos = 0.0_dp
697 : work_sin = 0.0_dp
698 : ai_work_cos = 0.0_dp
699 : ai_work_sin = 0.0_dp
700 : IF (derivative) THEN
701 : ai_work_dcos = 0.0_dp
702 : ai_work_dsin = 0.0_dp
703 : END IF
704 : dogth = .FALSE.
705 :
706 : NULLIFY (first_sgfa, la_max, la_min, npgfa, nsgf_seta, rpgfa, set_radius_a, sphi_a, zeta)
707 : NULLIFY (alpha_ppnl, cprj, nprj_ppnl, vprj_ppnl)
708 : NULLIFY (clist, clist_sin)
709 :
710 : !$OMP DO SCHEDULE(GUIDED)
711 : DO slot = 1, sap_ppnl(1)%nl_size
712 : ikind = sap_ppnl(1)%nlist_task(slot)%ikind
713 : kkind = sap_ppnl(1)%nlist_task(slot)%jkind
714 : iatom = sap_ppnl(1)%nlist_task(slot)%iatom
715 : katom = sap_ppnl(1)%nlist_task(slot)%jatom
716 : nlist = sap_ppnl(1)%nlist_task(slot)%nlist
717 : ilist = sap_ppnl(1)%nlist_task(slot)%ilist
718 : nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
719 : jneighbor = sap_ppnl(1)%nlist_task(slot)%inode
720 : cell_c(:) = sap_ppnl(1)%nlist_task(slot)%cell(:)
721 : rac(1:3) = sap_ppnl(1)%nlist_task(slot)%r(1:3)
722 : dac = NORM2(rac)
723 :
724 : iac = ikind + nkind*(kkind - 1)
725 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
726 : ! get definition of gto basis set
727 : first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf
728 : la_max => basis_set(ikind)%gto_basis_set%lmax
729 : la_min => basis_set(ikind)%gto_basis_set%lmin
730 : npgfa => basis_set(ikind)%gto_basis_set%npgf
731 : nseta = basis_set(ikind)%gto_basis_set%nset
732 : nsgfa = basis_set(ikind)%gto_basis_set%nsgf
733 : nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set
734 : rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius
735 : set_radius_a => basis_set(ikind)%gto_basis_set%set_radius
736 : sphi_a => basis_set(ikind)%gto_basis_set%sphi
737 : zeta => basis_set(ikind)%gto_basis_set%zet
738 :
739 : IF (ASSOCIATED(gpotential(kkind)%gth_potential)) THEN
740 : ! GTH potential
741 : dogth = .TRUE.
742 : alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl
743 : cprj => gpotential(kkind)%gth_potential%cprj
744 : lppnl = gpotential(kkind)%gth_potential%lppnl
745 : nppnl = gpotential(kkind)%gth_potential%nppnl
746 : nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl
747 : ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius
748 : vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl
749 : ELSE
750 : CYCLE
751 : END IF
752 :
753 : clist => sap_int_cos(iac)%alist(ilist)%clist(jneighbor)
754 : clist_sin => sap_int_sin(iac)%alist(ilist)%clist(jneighbor)
755 :
756 : clist%catom = katom
757 : clist%cell = cell_c
758 : clist%rac = rac
759 : clist_sin%catom = katom
760 : clist_sin%cell = cell_c
761 : clist_sin%rac = rac
762 :
763 : IF (.NOT. derivative) THEN
764 : ALLOCATE (clist%acint(nsgfa, nppnl, 1), clist%achint(nsgfa, nppnl, 1))
765 : ELSE
766 : ALLOCATE (clist%acint(nsgfa, nppnl, 4), clist%achint(nsgfa, nppnl, 4))
767 : END IF
768 : clist%acint = 0.0_dp
769 : clist%achint = 0.0_dp
770 : clist%nsgf_cnt = 0
771 :
772 : IF (.NOT. derivative) THEN
773 : ALLOCATE (clist_sin%acint(nsgfa, nppnl, 1), clist_sin%achint(nsgfa, nppnl, 1))
774 : ELSE
775 : ALLOCATE (clist_sin%acint(nsgfa, nppnl, 4), clist_sin%achint(nsgfa, nppnl, 4))
776 : END IF
777 : clist_sin%acint = 0.0_dp
778 : clist_sin%achint = 0.0_dp
779 : clist_sin%nsgf_cnt = 0
780 :
781 : ! reference point at zero
782 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
783 : rc(:) = ra + rac
784 :
785 : ! reference point at pseudized atom
786 : raf(:) = ra - rc
787 : rcf(:) = 0._dp
788 :
789 : DO iset = 1, nseta
790 : ncoa = npgfa(iset)*ncoset(la_max(iset))
791 : sgfa = first_sgfa(1, iset)
792 : IF (dogth) THEN
793 : prjc = 1
794 : work_cos = 0.0_dp
795 : work_sin = 0.0_dp
796 : DO l = 0, lppnl
797 : nprjc = nprj_ppnl(l)*nco(l)
798 : IF (nprjc == 0) CYCLE
799 : rprjc(1) = ppnl_radius
800 : IF (set_radius_a(iset) + rprjc(1) < dac) CYCLE
801 : lc_max = l + 2*(nprj_ppnl(l) - 1)
802 : lc_min = l
803 : zetc(1) = alpha_ppnl(l)
804 : ncoc = ncoset(lc_max)
805 :
806 : IF (.NOT. derivative) THEN
807 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
808 : lc_max, 1, zetc, rprjc, lc_min, raf, rcf, kvec, ai_work_cos, ai_work_sin)
809 : ELSE
810 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
811 : lc_max, 1, zetc, rprjc, lc_min, raf, rcf, kvec, ai_work_cos, ai_work_sin, &
812 : dcosab=ai_work_dcos, dsinab=ai_work_dsin)
813 : END IF
814 : ! projector functions: Cartesian -> spherical
815 : na = ncoa
816 : nb = nprjc
817 : np = ncoc
818 : work_cos(1:na, prjc:prjc + nb - 1) = &
819 : MATMUL(ai_work_cos(1:na, 1:np), cprj(1:np, prjc:prjc + nb - 1))
820 : work_sin(1:na, prjc:prjc + nb - 1) = &
821 : MATMUL(ai_work_sin(1:na, 1:np), cprj(1:np, prjc:prjc + nb - 1))
822 :
823 : IF (derivative) THEN
824 : DO idir = 1, 3
825 : work_dcos(1:na, prjc:prjc + nb - 1, idir) = &
826 : MATMUL(ai_work_dcos(1:na, 1:np, idir), cprj(1:np, prjc:prjc + nb - 1))
827 : work_dsin(1:na, prjc:prjc + nb - 1, idir) = &
828 : MATMUL(ai_work_dsin(1:na, 1:np, idir), cprj(1:np, prjc:prjc + nb - 1))
829 : END DO
830 : END IF
831 :
832 : prjc = prjc + nprjc
833 : END DO
834 :
835 : ! contract gto basis set into acint
836 : na = nsgf_seta(iset)
837 : nb = nppnl
838 : np = ncoa
839 : clist%acint(sgfa:sgfa + na - 1, 1:nb, 1) = &
840 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_cos(1:np, 1:nb))
841 : clist_sin%acint(sgfa:sgfa + na - 1, 1:nb, 1) = &
842 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_sin(1:np, 1:nb))
843 : IF (derivative) THEN
844 : DO idir = 1, 3
845 : clist%acint(sgfa:sgfa + na - 1, 1:nb, 1 + idir) = &
846 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_dcos(1:np, 1:nb, idir))
847 : clist_sin%acint(sgfa:sgfa + na - 1, 1:nb, 1 + idir) = &
848 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_dsin(1:np, 1:nb, idir))
849 : END DO
850 : END IF
851 :
852 : ! multiply with interaction matrix h_ij of the nl pp
853 : clist%achint(sgfa:sgfa + na - 1, 1:nb, 1) = &
854 : MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, 1), vprj_ppnl(1:nb, 1:nb))
855 : clist_sin%achint(sgfa:sgfa + na - 1, 1:nb, 1) = &
856 : MATMUL(clist_sin%acint(sgfa:sgfa + na - 1, 1:nb, 1), vprj_ppnl(1:nb, 1:nb))
857 : IF (derivative) THEN
858 : DO idir = 1, 3
859 : clist%achint(sgfa:sgfa + na - 1, 1:nb, 1 + idir) = &
860 : MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, 1 + idir), vprj_ppnl(1:nb, 1:nb))
861 : clist_sin%achint(sgfa:sgfa + na - 1, 1:nb, 1 + idir) = &
862 : MATMUL(clist_sin%acint(sgfa:sgfa + na - 1, 1:nb, 1 + idir), vprj_ppnl(1:nb, 1:nb))
863 : END DO
864 : END IF
865 : END IF
866 :
867 : END DO
868 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
869 : clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1)))
870 : clist_sin%maxac = MAXVAL(ABS(clist_sin%acint(:, :, 1)))
871 : clist_sin%maxach = MAXVAL(ABS(clist_sin%achint(:, :, 1)))
872 : END DO
873 :
874 : DEALLOCATE (work_cos, work_sin, ai_work_cos, ai_work_sin)
875 : IF (derivative) DEALLOCATE (work_dcos, work_dsin, ai_work_dcos, ai_work_dsin)
876 :
877 : !$OMP END PARALLEL
878 :
879 62 : DEALLOCATE (gpotential, spotential)
880 :
881 62 : CALL timestop(handle)
882 :
883 124 : END SUBROUTINE build_sap_exp_ints
884 :
885 : ! **************************************************************************************************
886 : !> \brief Calculate the force associated to non-local pseudo potential in the velocity gauge
887 : !> \param qs_env ...
888 : !> \param particle_set ...
889 : !> \date 09.2023
890 : !> \author Guillaume Le Breton
891 : ! **************************************************************************************************
892 22 : SUBROUTINE velocity_gauge_nl_force(qs_env, particle_set)
893 : TYPE(qs_environment_type), POINTER :: qs_env
894 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
895 :
896 : CHARACTER(len=*), PARAMETER :: routiuneN = "velocity_gauge_nl_force"
897 :
898 : INTEGER :: handle, i, iac, iatom, ibc, icol, idir, ikind, irow, jatom, jkind, kac, katom, &
899 : kbc, kkind, maxl, maxlgto, maxlppnl, na, natom, nb, nkind, np, slot
900 22 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
901 : INTEGER, DIMENSION(3) :: cell_b
902 : LOGICAL :: found_imag, found_real
903 : REAL(dp) :: eps_ppnl, f0, sign_imag
904 : REAL(KIND=dp), DIMENSION(3) :: fa, fb, rab, vec_pot
905 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
906 22 : POINTER :: sab_orb, sap_ppnl
907 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
908 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
909 22 : DIMENSION(:) :: basis_set
910 : TYPE(dft_control_type), POINTER :: dft_control
911 22 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_ao_im
912 : TYPE(cell_type), POINTER :: cell
913 22 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
914 : TYPE(alist_type), POINTER :: alist_cos_ac, alist_cos_bc, &
915 : alist_sin_ac, alist_sin_bc
916 22 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint_cos, achint_sin, acint_cos, &
917 22 : acint_sin, bchint_cos, bchint_sin, &
918 22 : bcint_cos, bcint_sin
919 22 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix_p_imag, matrix_p_real
920 44 : REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
921 22 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
922 22 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
923 : TYPE(qs_rho_type), POINTER :: rho
924 22 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int_cos, sap_int_sin
925 :
926 22 : CALL timeset(routiuneN, handle)
927 :
928 22 : NULLIFY (sap_ppnl)
929 :
930 : CALL get_qs_env(qs_env, &
931 22 : sap_ppnl=sap_ppnl)
932 :
933 22 : IF (ASSOCIATED(sap_ppnl)) THEN
934 22 : NULLIFY (qs_kind_set, cell, dft_control, force, sab_orb, atomic_kind_set, &
935 22 : sap_int_cos, sap_int_sin)
936 : ! Load and initialized the required quantities
937 :
938 : CALL get_qs_env(qs_env, &
939 : sab_orb=sab_orb, &
940 : force=force, &
941 : dft_control=dft_control, &
942 : qs_kind_set=qs_kind_set, &
943 : cell=cell, &
944 : atomic_kind_set=atomic_kind_set, &
945 22 : rho=rho)
946 :
947 22 : nkind = SIZE(atomic_kind_set)
948 22 : natom = SIZE(particle_set)
949 22 : eps_ppnl = dft_control%qs_control%eps_ppnl
950 :
951 : CALL get_qs_kind_set(qs_kind_set, &
952 : maxlgto=maxlgto, &
953 22 : maxlppnl=maxlppnl)
954 :
955 22 : maxl = MAX(maxlppnl, maxlgto)
956 22 : CALL init_orbital_pointers(maxl + 1)
957 :
958 : ! initalize sab_int types to store the integrals
959 286 : ALLOCATE (sap_int_cos(nkind*nkind), sap_int_sin(nkind*nkind))
960 110 : DO i = 1, SIZE(sap_int_cos)
961 88 : NULLIFY (sap_int_cos(i)%alist, sap_int_cos(i)%asort, sap_int_cos(i)%aindex)
962 88 : sap_int_cos(i)%nalist = 0
963 88 : NULLIFY (sap_int_sin(i)%alist, sap_int_sin(i)%asort, sap_int_sin(i)%aindex)
964 110 : sap_int_sin(i)%nalist = 0
965 : END DO
966 :
967 : ! get basis set
968 110 : ALLOCATE (basis_set(nkind))
969 66 : DO ikind = 1, nkind
970 44 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
971 66 : IF (ASSOCIATED(orb_basis_set)) THEN
972 44 : basis_set(ikind)%gto_basis_set => orb_basis_set
973 : ELSE
974 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
975 : END IF
976 : END DO
977 :
978 : !get vector potential
979 88 : vec_pot = dft_control%rtp_control%vec_pot
980 :
981 286 : force_thread = 0.0_dp
982 :
983 22 : CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao, rho_ao_im=rho_ao_im)
984 : ! To avoid FOR loop over spin, sum the 2 spin into the first one directly. Undone later on
985 22 : IF (SIZE(rho_ao) == 2) THEN
986 : CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, &
987 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
988 : CALL dbcsr_add(rho_ao_im(1)%matrix, rho_ao_im(2)%matrix, &
989 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
990 : END IF
991 :
992 : ! Compute cosap = <a|cos kr|p>, sindap = <a|sin kr|p>, cosdap = <da/dRA|cos kr|p>, and sindap = <da/dRA|sin kr|p>
993 : CALL build_sap_exp_ints(sap_int_cos, sap_int_sin, sap_ppnl, qs_kind_set, particle_set, &
994 22 : cell, kvec=vec_pot, basis_set=basis_set, nkind=nkind, derivative=.TRUE.)
995 22 : CALL sap_sort(sap_int_cos)
996 22 : CALL sap_sort(sap_int_sin)
997 :
998 : ! Compute the force, on nuclei A it is given by: Re(P_ab) Re(dV_ab/dRA) - Im(P_ab) Im(dV_ab/dRA)
999 :
1000 : !$OMP PARALLEL &
1001 : !$OMP DEFAULT (NONE) &
1002 : !$OMP SHARED (basis_set, sab_orb, sap_int_cos, sap_int_sin, eps_ppnl, nkind, natom,&
1003 : !$OMP rho_ao, rho_ao_im) &
1004 : !$OMP PRIVATE (matrix_p_real, matrix_p_imag, acint_cos, achint_cos, bcint_cos, bchint_cos, acint_sin,&
1005 : !$OMP achint_sin, bcint_sin, bchint_sin, slot, ikind, jkind, iatom, jatom,&
1006 : !$OMP cell_b, rab, irow, icol, fa, fb, f0, found_real, found_imag, sign_imag, &
1007 : !$OMP kkind, iac, ibc, alist_cos_ac, alist_cos_bc, alist_sin_ac, alist_sin_bc, kac, kbc,&
1008 : !$OMP na, np, nb, katom) &
1009 22 : !$OMP REDUCTION (+ : force_thread )
1010 :
1011 : NULLIFY (acint_cos, bcint_cos, achint_cos, bchint_cos)
1012 : NULLIFY (acint_sin, bcint_sin, achint_sin, bchint_sin)
1013 :
1014 : ! loop over atom pairs
1015 : !$OMP DO SCHEDULE(GUIDED)
1016 : DO slot = 1, sab_orb(1)%nl_size
1017 : ikind = sab_orb(1)%nlist_task(slot)%ikind
1018 : jkind = sab_orb(1)%nlist_task(slot)%jkind
1019 : iatom = sab_orb(1)%nlist_task(slot)%iatom
1020 : jatom = sab_orb(1)%nlist_task(slot)%jatom
1021 : cell_b(:) = sab_orb(1)%nlist_task(slot)%cell
1022 : rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
1023 :
1024 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
1025 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
1026 :
1027 : ! Use the symmetry of the first derivatives
1028 : IF (iatom == jatom) THEN
1029 : f0 = 1.0_dp
1030 : ELSE
1031 : f0 = 2.0_dp
1032 : END IF
1033 :
1034 : fa = 0.0_dp
1035 : fb = 0.0_dp
1036 :
1037 : IF (iatom <= jatom) THEN
1038 : irow = iatom
1039 : icol = jatom
1040 : sign_imag = +1.0_dp
1041 : ELSE
1042 : irow = jatom
1043 : icol = iatom
1044 : sign_imag = -1.0_dp
1045 : END IF
1046 : NULLIFY (matrix_p_real, matrix_p_imag)
1047 : CALL dbcsr_get_block_p(rho_ao(1)%matrix, irow, icol, matrix_p_real, found_real)
1048 : CALL dbcsr_get_block_p(rho_ao_im(1)%matrix, irow, icol, matrix_p_imag, found_imag)
1049 :
1050 : IF (found_real .OR. found_imag) THEN
1051 : ! loop over the <gto_a|ppln_c>h_ij<ppnl_c|gto_b> pairs
1052 : DO kkind = 1, nkind
1053 : iac = ikind + nkind*(kkind - 1)
1054 : ibc = jkind + nkind*(kkind - 1)
1055 : IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist)) CYCLE
1056 : IF (.NOT. ASSOCIATED(sap_int_cos(ibc)%alist)) CYCLE
1057 : IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist)) CYCLE
1058 : IF (.NOT. ASSOCIATED(sap_int_sin(ibc)%alist)) CYCLE
1059 : CALL get_alist(sap_int_cos(iac), alist_cos_ac, iatom)
1060 : CALL get_alist(sap_int_cos(ibc), alist_cos_bc, jatom)
1061 : CALL get_alist(sap_int_sin(iac), alist_sin_ac, iatom)
1062 : CALL get_alist(sap_int_sin(ibc), alist_sin_bc, jatom)
1063 : IF (.NOT. ASSOCIATED(alist_cos_ac)) CYCLE
1064 : IF (.NOT. ASSOCIATED(alist_cos_bc)) CYCLE
1065 : IF (.NOT. ASSOCIATED(alist_sin_ac)) CYCLE
1066 : IF (.NOT. ASSOCIATED(alist_sin_bc)) CYCLE
1067 :
1068 : ! only use cos for indexing, as cos and sin integrals are constructed by the same routine
1069 : ! in the same way
1070 : DO kac = 1, alist_cos_ac%nclist
1071 : DO kbc = 1, alist_cos_bc%nclist
1072 : ! the next two ifs should be the same for sine integrals
1073 : IF (alist_cos_ac%clist(kac)%catom /= alist_cos_bc%clist(kbc)%catom) CYCLE
1074 : IF (ALL(cell_b + alist_cos_bc%clist(kbc)%cell - alist_cos_ac%clist(kac)%cell == 0)) THEN
1075 : ! screening
1076 : IF (alist_cos_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
1077 : .AND. alist_cos_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl &
1078 : .AND. alist_sin_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
1079 : .AND. alist_sin_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
1080 :
1081 : acint_cos => alist_cos_ac%clist(kac)%acint
1082 : bcint_cos => alist_cos_bc%clist(kbc)%acint
1083 : achint_cos => alist_cos_ac%clist(kac)%achint
1084 : bchint_cos => alist_cos_bc%clist(kbc)%achint
1085 : acint_sin => alist_sin_ac%clist(kac)%acint
1086 : bcint_sin => alist_sin_bc%clist(kbc)%acint
1087 : achint_sin => alist_sin_ac%clist(kac)%achint
1088 : bchint_sin => alist_sin_bc%clist(kbc)%achint
1089 :
1090 : na = SIZE(acint_cos, 1)
1091 : np = SIZE(acint_cos, 2)
1092 : nb = SIZE(bcint_cos, 1)
1093 : ! Re(dV_ab/dRA) = <da/dRA|cos kr|p><p|cos kr|b> + <db/dRA|cos kr|p><p|cos kr|a> + <da/dRA|sin kr|p><p|sin kr|b> + <db/dRA|sin kr|p><p|sin|a>
1094 : ! Im(dV_ab/dRA) = <da/dRA|sin kr|p><p|cos kr|b> - <db/dRA|sin kr|p><p|cos kr|a> - <da/dRA|cos kr|p><p|sin kr|b> + <db/dRA|cos kr|p><p|sin|a>
1095 : katom = alist_cos_ac%clist(kac)%catom
1096 : DO idir = 1, 3
1097 : IF (iatom <= jatom) THEN
1098 : ! For fa:
1099 : IF (found_real) &
1100 : fa(idir) = SUM(matrix_p_real(1:na, 1:nb)* &
1101 : (+MATMUL(acint_cos(1:na, 1:np, 1 + idir), TRANSPOSE(bchint_cos(1:nb, 1:np, 1))) &
1102 : + MATMUL(acint_sin(1:na, 1:np, 1 + idir), TRANSPOSE(bchint_sin(1:nb, 1:np, 1)))))
1103 : IF (found_imag) &
1104 : fa(idir) = fa(idir) - sign_imag*SUM(matrix_p_imag(1:na, 1:nb)* &
1105 : (+MATMUL(acint_sin(1:na, 1:np, 1 + idir), TRANSPOSE(bchint_cos(1:nb, 1:np, 1))) &
1106 : - MATMUL(acint_cos(1:na, 1:np, 1 + idir), TRANSPOSE(bchint_sin(1:nb, 1:np, 1)))))
1107 : ! For fb:
1108 : IF (found_real) &
1109 : fb(idir) = SUM(matrix_p_real(1:na, 1:nb)* &
1110 : (+MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, 1 + idir))) &
1111 : + MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, 1 + idir)))))
1112 : IF (found_imag) &
1113 : fb(idir) = fb(idir) - sign_imag*SUM(matrix_p_imag(1:na, 1:nb)* &
1114 : (-MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, 1 + idir))) &
1115 : + MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, 1 + idir)))))
1116 : ELSE
1117 : ! For fa:
1118 : IF (found_real) &
1119 : fa(idir) = SUM(matrix_p_real(1:nb, 1:na)* &
1120 : (+MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, 1 + idir))) &
1121 : + MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, 1 + idir)))))
1122 : IF (found_imag) &
1123 : fa(idir) = fa(idir) - sign_imag*SUM(matrix_p_imag(1:nb, 1:na)* &
1124 : (+MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, 1 + idir))) &
1125 : - MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, 1 + idir)))))
1126 : ! For fb
1127 : IF (found_real) &
1128 : fb(idir) = SUM(matrix_p_real(1:nb, 1:na)* &
1129 : (+MATMUL(bcint_cos(1:nb, 1:np, 1 + idir), TRANSPOSE(achint_cos(1:na, 1:np, 1))) &
1130 : + MATMUL(bcint_sin(1:nb, 1:np, 1 + idir), TRANSPOSE(achint_sin(1:na, 1:np, 1)))))
1131 : IF (found_imag) &
1132 : fb(idir) = fb(idir) - sign_imag*SUM(matrix_p_imag(1:nb, 1:na)* &
1133 : (-MATMUL(bcint_cos(1:nb, 1:np, 1 + idir), TRANSPOSE(achint_sin(1:na, 1:np, 1))) &
1134 : + MATMUL(bcint_sin(1:nb, 1:np, 1 + idir), TRANSPOSE(achint_cos(1:na, 1:np, 1)))))
1135 : END IF
1136 : force_thread(idir, iatom) = force_thread(idir, iatom) + f0*fa(idir)
1137 : force_thread(idir, katom) = force_thread(idir, katom) - f0*fa(idir)
1138 : force_thread(idir, jatom) = force_thread(idir, jatom) + f0*fb(idir)
1139 : force_thread(idir, katom) = force_thread(idir, katom) - f0*fb(idir)
1140 : END DO
1141 : EXIT
1142 : END IF
1143 : END DO
1144 : END DO
1145 : END DO
1146 : END IF
1147 :
1148 : END DO
1149 :
1150 : !$OMP END PARALLEL
1151 :
1152 : ! Update the force
1153 22 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
1154 : !$OMP DO
1155 : DO iatom = 1, natom
1156 66 : i = atom_of_kind(iatom)
1157 66 : ikind = kind_of(iatom)
1158 264 : force(ikind)%gth_ppnl(:, i) = force(ikind)%gth_ppnl(:, i) + force_thread(:, iatom)
1159 : END DO
1160 : !$OMP END DO
1161 :
1162 : ! Clean up
1163 22 : IF (SIZE(rho_ao) == 2) THEN
1164 : CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, &
1165 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1166 : CALL dbcsr_add(rho_ao_im(1)%matrix, rho_ao_im(2)%matrix, &
1167 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1168 : END IF
1169 22 : CALL release_sap_int(sap_int_cos)
1170 22 : CALL release_sap_int(sap_int_sin)
1171 :
1172 44 : DEALLOCATE (basis_set, atom_of_kind, kind_of)
1173 :
1174 : END IF
1175 :
1176 22 : CALL timestop(handle)
1177 :
1178 44 : END SUBROUTINE velocity_gauge_nl_force
1179 :
1180 : END MODULE rt_propagation_velocity_gauge
|
