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