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