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