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