(git:374b731)
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-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
12MODULE rt_propagation_velocity_gauge
13 USE ai_moments, ONLY: cossin
14 USE atomic_kind_types, ONLY: atomic_kind_type,&
15 get_atomic_kind_set
16 USE basis_set_types, ONLY: gto_basis_set_p_type,&
17 gto_basis_set_type
18 USE bibliography, ONLY: mattiat2022,&
19 cite_reference
20 USE cell_types, ONLY: cell_type,&
21 pbc
22 USE core_ppnl, ONLY: build_core_ppnl
23 USE cp_control_types, ONLY: dft_control_type
24 USE cp_dbcsr_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
85CONTAINS
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
1180END MODULE rt_propagation_velocity_gauge
cell_types::pbc
Definition cell_types.F:92
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
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:93
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:559
external_potential_types::gth_potential_p_type
Definition external_potential_types.F:227
external_potential_types::gth_potential_type
Definition external_potential_types.F:116
external_potential_types::sgp_potential_p_type
Definition external_potential_types.F:235
external_potential_types::sgp_potential_type
Definition external_potential_types.F:172
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:138
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
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:171
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:135
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
sap_kind_types::alist_type
Definition sap_kind_types.F:58
sap_kind_types::clist_type
Definition sap_kind_types.F:43
sap_kind_types::sap_int_type
Definition sap_kind_types.F:64
virial_types::virial_type
Definition virial_types.F:26