(git:55939f2)
Loading...
Searching...
No Matches
qs_moments.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 Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
10!> \par History
11!> added angular moments (JGH 11.2012)
12!> \author JGH (20.07.2006)
13! **************************************************************************************************
14MODULE qs_moments
15 USE ai_angmom, ONLY: angmom
16 USE ai_moments, ONLY: contract_expikr, &
17 diff_momop, &
18 expikr, &
19 moment
20 USE atomic_kind_types, ONLY: atomic_kind_type, &
21 get_atomic_kind
22 USE basis_set_types, ONLY: gto_basis_set_p_type, &
23 gto_basis_set_type
24 USE bibliography, ONLY: mattiat2019, &
25 cite_reference
26 USE block_p_types, ONLY: block_p_type
27 USE cell_types, ONLY: cell_type, &
28 pbc, &
29 pbc_stable, &
30 get_cell
31 USE commutator_rpnl, ONLY: build_com_mom_nl
32 USE cp_blacs_env, ONLY: cp_blacs_env_type
33 USE cp_cfm_basic_linalg, ONLY: cp_cfm_det
34 USE cp_cfm_types, ONLY: cp_cfm_copy_general, &
35 cp_cfm_create, &
36 cp_cfm_get_element, &
37 cp_cfm_get_info, &
38 cp_cfm_release, &
39 cp_cfm_to_fm, &
40 cp_fm_to_cfm, &
41 cp_cfm_type
42 USE cp_control_types, ONLY: dft_control_type
43 USE cp_dbcsr_api, ONLY: &
44 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_type, &
45 dbcsr_get_block_p, dbcsr_get_info, dbcsr_get_matrix_type, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, &
46 dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
47 USE cp_dbcsr_contrib, ONLY: dbcsr_dot, &
48 dbcsr_trace
49 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
50 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm, &
51 cp_dbcsr_sm_fm_multiply, &
52 dbcsr_allocate_matrix_set, &
53 dbcsr_deallocate_matrix_set
54 USE cp_files, ONLY: close_file, &
55 open_file
56 USE cp_fm_struct, ONLY: cp_fm_struct_create, &
57 cp_fm_struct_double, &
58 cp_fm_struct_release, &
59 cp_fm_struct_type
60 USE cp_fm_types, ONLY: cp_fm_copy_general, &
61 cp_fm_create, &
62 cp_fm_get_info, &
63 cp_fm_release, &
64 cp_fm_set_all, &
65 cp_fm_type
66 USE cp_result_methods, ONLY: cp_results_erase, &
67 put_results
68 USE cp_result_types, ONLY: cp_result_type
69 USE cp_units, ONLY: cp_ukind_dipole, &
70 cp_unit_create, &
71 cp_unit_from_cp2k, &
72 cp_unit_from_cp2k1, &
73 cp_unit_release, &
74 cp_unit_type, &
75 cp_units_au, &
76 cp_units_c_m, &
77 cp_units_debye, &
78 cp_units_e_angstrom
79 USE distribution_1d_types, ONLY: distribution_1d_type
80 USE input_constants, ONLY: use_mom_ref_coac
81 USE kinds, ONLY: default_string_length, &
82 dp, &
83 max_line_length
84 USE kpoint_k_r_trafo_simple, ONLY: replicate_rs_matrices, &
85 rs_to_kp
86 USE kpoint_methods, ONLY: kpoint_init_cell_index, &
87 kpoint_initialize, &
88 rskp_transform
89 USE kpoint_types, ONLY: get_kpoint_info, &
90 kpoint_env_p_type, &
91 kpoint_env_type, &
92 kpoint_type, &
93 kpoint_create, &
94 kpoint_release, &
95 read_kpoint_section
96 USE mathconstants, ONLY: pi, &
97 twopi, &
98 gaussi, &
99 z_one, &
100 z_zero
101 USE message_passing, ONLY: mp_para_env_type
102 USE parallel_gemm_api, ONLY: parallel_gemm
103 USE moments_utils, ONLY: get_reference_point
104 USE orbital_pointers, ONLY: current_maxl, &
105 indco, &
106 ncoset
107 USE particle_methods, ONLY: get_particle_set
108 USE particle_types, ONLY: particle_type
109 USE physcon, ONLY: angstrom
110 USE qs_environment_types, ONLY: get_qs_env, &
111 qs_environment_type
112 USE qs_kind_types, ONLY: get_qs_kind, &
113 get_qs_kind_set, &
114 qs_kind_type
115 USE qs_ks_types, ONLY: get_ks_env, &
116 qs_ks_env_type
117 USE qs_mo_types, ONLY: get_mo_set, &
118 mo_set_type
119 USE qs_neighbor_list_types, ONLY: get_iterator_info, &
120 neighbor_list_iterate, &
121 neighbor_list_iterator_create, &
122 neighbor_list_iterator_p_type, &
123 neighbor_list_iterator_release, &
124 neighbor_list_set_p_type
125 USE qs_operators_ao, ONLY: build_exp_ikr_matrix, &
126 build_lin_mom_matrix
127 USE qs_overlap, ONLY: build_overlap_matrix
128 USE qs_rho_types, ONLY: qs_rho_get, &
129 qs_rho_type
130 USE rt_propagation_types, ONLY: get_rtp, &
131 rt_prop_type
132 USE soc_pseudopotential_methods, ONLY: v_soc_xyz_from_pseudopotential
133 USE cp_parser_methods, ONLY: read_float_object
134 USE input_section_types, ONLY: section_vals_get, &
135 section_vals_get_subs_vals, &
136 section_vals_type, &
137 section_vals_val_get
138 USE mathlib, ONLY: geeig_right, &
139 gemm_square
140 USE string_utilities, ONLY: uppercase
141
142#include "./base/base_uses.f90"
143
144 IMPLICIT NONE
145
146 PRIVATE
147
148 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_moments'
149
150 ! Public subroutines
151 PUBLIC :: build_berry_moment_matrix, build_local_moment_matrix, get_multipoles_ao
152 PUBLIC :: build_berry_kpoint_matrix, build_local_magmom_matrix
153 PUBLIC :: qs_moment_berry_phase, qs_moment_locop
154 PUBLIC :: qs_moment_kpoints, build_local_moment_matrix_rs_img
155 PUBLIC :: qs_moment_kpoints_deep, qs_moment_kpoints_deep_soc
156 PUBLIC :: qs_moment_kpoints_scf_mos
157 PUBLIC :: build_local_moments_der_matrix
158 PUBLIC :: calculate_commutator_nl_terms, op_orbbas, op_orbbas_rtp, &
159 print_moments, print_moments_nl, set_label
160
161CONTAINS
162
163! **************************************************************************************************
164!> \brief ...
165!> \param qs_env ...
166!> \param moments ...
167!> \param nmoments ...
168!> \param ref_point ...
169!> \param ref_points ...
170!> \param basis_type ...
171!> \param all_images accumulate all neighbor images belonging to one atom pair
172!> \param minimum_image apply the half-cell image filter
173!> \param neighbor_image preserve the neighbor-list image geometry
174!> \param first_component first Cartesian moment component to write
175! **************************************************************************************************
176 SUBROUTINE build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type, &
177 all_images, minimum_image, neighbor_image, first_component)
178
179 TYPE(qs_environment_type), POINTER :: qs_env
180 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moments
181 INTEGER, INTENT(IN) :: nmoments
182 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
183 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
184 OPTIONAL :: ref_points
185 CHARACTER(len=*), OPTIONAL :: basis_type
186 LOGICAL, OPTIONAL :: all_images, minimum_image, neighbor_image
187 INTEGER, OPTIONAL :: first_component
188
189 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moment_matrix'
190
191 INTEGER :: component_start, handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, &
192 last_jatom, maxco, maxsgf, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
193 LOGICAL :: found, my_all_images, my_minimum_image, &
194 my_neighbor_image, new_atom_pair
195 REAL(kind=dp) :: dab
196 REAL(kind=dp), DIMENSION(3) :: half_cell
197 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
198 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
199 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
200 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
201 TYPE(cell_type), POINTER :: cell
202 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
203 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
204 TYPE(neighbor_list_iterator_p_type), &
205 DIMENSION(:), POINTER :: nl_iterator
206 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
207 POINTER :: sab_orb
208 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
209 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
210 TYPE(qs_kind_type), POINTER :: qs_kind
211
212 IF (nmoments < 1) RETURN
213
214 CALL timeset(routinen, handle)
215
216 my_all_images = .false.
217 IF (PRESENT(all_images)) my_all_images = all_images
218 my_minimum_image = .false.
219 IF (PRESENT(minimum_image)) my_minimum_image = minimum_image
220 my_neighbor_image = .false.
221 IF (PRESENT(neighbor_image)) my_neighbor_image = neighbor_image
222
223 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
224 cpassert(SIZE(moments) >= nm)
225 component_start = 1
226 IF (PRESENT(first_component)) component_start = first_component
227 cpassert(component_start >= 1 .AND. component_start <= nm)
228
229 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
230 CALL get_qs_env(qs_env=qs_env, &
231 qs_kind_set=qs_kind_set, &
232 particle_set=particle_set, cell=cell, &
233 sab_orb=sab_orb)
234
235 IF (my_minimum_image) THEN
236 half_cell = 0.5_dp*norm2(cell%hmat, dim=1)
237 END IF
238
239 nkind = SIZE(qs_kind_set)
240 ! Allocate work storage
241 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
242 maxco=maxco, maxsgf=maxsgf, &
243 basis_type=basis_type)
244
245 ALLOCATE (mab(maxco, maxco, nm))
246 mab(:, :, :) = 0.0_dp
247
248 ALLOCATE (work(maxco, maxsgf))
249 work(:, :) = 0.0_dp
250
251 ALLOCATE (mint(nm))
252 DO i = component_start, nm
253 NULLIFY (mint(i)%block)
254 END DO
255
256 ALLOCATE (basis_set_list(nkind))
257 DO ikind = 1, nkind
258 qs_kind => qs_kind_set(ikind)
259 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
260 IF (ASSOCIATED(basis_set_a)) THEN
261 basis_set_list(ikind)%gto_basis_set => basis_set_a
262 ELSE
263 NULLIFY (basis_set_list(ikind)%gto_basis_set)
264 END IF
265 END DO
266 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
267 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
268 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
269 iatom=iatom, jatom=jatom, r=rab)
270 basis_set_a => basis_set_list(ikind)%gto_basis_set
271 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
272 basis_set_b => basis_set_list(jkind)%gto_basis_set
273 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
274 associate( &
275 ! basis ikind
276 first_sgfa => basis_set_a%first_sgf, &
277 la_max => basis_set_a%lmax, &
278 la_min => basis_set_a%lmin, &
279 npgfa => basis_set_a%npgf, &
280 nsgfa => basis_set_a%nsgf_set, &
281 rpgfa => basis_set_a%pgf_radius, &
282 set_radius_a => basis_set_a%set_radius, &
283 sphi_a => basis_set_a%sphi, &
284 zeta => basis_set_a%zet, &
285 ! basis jkind, &
286 first_sgfb => basis_set_b%first_sgf, &
287 lb_max => basis_set_b%lmax, &
288 lb_min => basis_set_b%lmin, &
289 npgfb => basis_set_b%npgf, &
290 nsgfb => basis_set_b%nsgf_set, &
291 rpgfb => basis_set_b%pgf_radius, &
292 set_radius_b => basis_set_b%set_radius, &
293 sphi_b => basis_set_b%sphi, &
294 zetb => basis_set_b%zet)
295
296 nseta = basis_set_a%nset
297 nsetb = basis_set_b%nset
298
299 IF (inode == 1) last_jatom = 0
300
301 IF (my_minimum_image) THEN
302 IF (any(abs(rab(:)) > half_cell(:))) cycle
303 END IF
304
305 new_atom_pair = jatom /= last_jatom
306 IF (.NOT. my_all_images .AND. .NOT. new_atom_pair) THEN
307 cycle
308 END IF
309
310 IF (new_atom_pair) THEN
311 last_jatom = jatom
312
313 IF (iatom <= jatom) THEN
314 irow = iatom
315 icol = jatom
316 ELSE
317 irow = jatom
318 icol = iatom
319 END IF
320
321 DO i = component_start, nm
322 NULLIFY (mint(i)%block)
323 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
324 row=irow, col=icol, block=mint(i)%block, found=found)
325 cpassert(found .AND. ASSOCIATED(mint(i)%block))
326 mint(i)%block = 0._dp
327 END DO
328 END IF
329
330 ! Reference point and image geometry.
331 IF (PRESENT(ref_points)) THEN
332 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
333 ELSE IF (PRESENT(ref_point)) THEN
334 rc(:) = ref_point(:)
335 ELSE
336 rc(:) = 0._dp
337 END IF
338
339 IF (my_all_images) THEN
340 ! Keep the image represented by the current neighbor-list node.
341 ra(:) = pbc(particle_set(iatom)%r(:), cell)
342 rb(:) = ra(:) + rab(:)
343 rac(:) = pbc(rc, ra, cell)
344 IF (my_neighbor_image) THEN
345 rbc(:) = rac(:) + rab(:)
346 ELSE
347 rbc(:) = pbc(rc, rb, cell)
348 END IF
349 dab = norm2(rab)
350 ELSE IF (my_neighbor_image) THEN
351 ! Keep the neighbor-list image while selecting one atom-pair representative.
352 ra(:) = particle_set(iatom)%r(:)
353 rb(:) = ra(:) + rab(:)
354 rac(:) = pbc(rc, ra, cell)
355 rbc(:) = rac(:) + rab(:)
356 dab = norm2(rab)
357 ELSE
358 ! Fold atomic positions around the reference point.
359 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
360 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
361 rab(:) = ra(:) - rb(:)
362 rac(:) = ra(:) - rc(:)
363 rbc(:) = rb(:) - rc(:)
364 dab = norm2(rab)
365 END IF
366
367 DO iset = 1, nseta
368
369 ncoa = npgfa(iset)*ncoset(la_max(iset))
370 sgfa = first_sgfa(1, iset)
371
372 DO jset = 1, nsetb
373
374 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
375
376 ncob = npgfb(jset)*ncoset(lb_max(jset))
377 sgfb = first_sgfb(1, jset)
378
379 ! Calculate the primitive integrals
380 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
381 rpgfa(:, iset), la_min(iset), &
382 lb_max(jset), npgfb(jset), zetb(:, jset), &
383 rpgfb(:, jset), nmoments, rac, rbc, mab)
384
385 ! Contraction step
386 DO i = component_start, nm
387
388 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
389 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
390 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
391 0.0_dp, work(1, 1), SIZE(work, 1))
392
393 IF (iatom <= jatom) THEN
394
395 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
396 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
397 work(1, 1), SIZE(work, 1), &
398 1.0_dp, mint(i)%block(sgfa, sgfb), &
399 SIZE(mint(i)%block, 1))
400
401 ELSE
402
403 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
404 1.0_dp, work(1, 1), SIZE(work, 1), &
405 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
406 1.0_dp, mint(i)%block(sgfb, sgfa), &
407 SIZE(mint(i)%block, 1))
408
409 END IF
410
411 END DO
412
413 END DO
414 END DO
415 END associate
416
417 END DO
418 CALL neighbor_list_iterator_release(nl_iterator)
419
420 ! Release work storage
421 DEALLOCATE (mab, basis_set_list)
422 DEALLOCATE (work)
423 DO i = component_start, nm
424 NULLIFY (mint(i)%block)
425 END DO
426 DEALLOCATE (mint)
427
428 CALL timestop(handle)
429
430 END SUBROUTINE build_local_moment_matrix
431
432! **************************************************************************************************
433!> \brief The AO multipoles M^k_µν = <φ_µ|(r - r_0)^k|φ_ν> of every order up to n_moments, on the
434!> block pattern of the overlap, about r_0, the centre of atomic charges;
435!> k = x, y, z for order 1, then xx, xy, xz, yy, yz, zz for order 2
436!> \param qs_env source of the basis, the geometry and the overlap's block pattern
437!> \param n_moments highest order, 1 for dipoles, 2 adds quadrupoles
438!> \param matrix_multipole M^k_µν, one entry per component k, allocated here; the caller releases it
439!> \param rpoint r_0
440! **************************************************************************************************
441 SUBROUTINE get_multipoles_ao(qs_env, n_moments, matrix_multipole, rpoint)
442
443 TYPE(qs_environment_type), POINTER :: qs_env
444 INTEGER, INTENT(IN) :: n_moments
445 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_multipole
446 REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: rpoint
447
448 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_multipoles_ao'
449
450 INTEGER :: handle, idir, n_multipole
451 REAL(kind=dp), DIMENSION(:), POINTER :: ref_point
452 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
453 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
454 POINTER :: sab_orb
455
456 CALL timeset(routinen, handle)
457
458 NULLIFY (matrix_s, sab_orb, ref_point)
459 CALL get_qs_env(qs_env, matrix_s=matrix_s, sab_orb=sab_orb)
460
461 ! components of every order up to n_moments: 3 dipoles, then 6 quadrupoles
462 n_multipole = (6 + 11*n_moments + 6*n_moments**2 + n_moments**3)/6 - 1
463 NULLIFY (matrix_multipole)
464 CALL dbcsr_allocate_matrix_set(matrix_multipole, n_multipole)
465 DO idir = 1, n_multipole
466 CALL dbcsr_init_p(matrix_multipole(idir)%matrix)
467 CALL dbcsr_create(matrix_multipole(idir)%matrix, name="ao_multipole", &
468 template=matrix_s(1)%matrix, matrix_type=dbcsr_type_symmetric)
469 CALL cp_dbcsr_alloc_block_from_nbl(matrix_multipole(idir)%matrix, sab_orb)
470 CALL dbcsr_set(matrix_multipole(idir)%matrix, 0._dp)
471 END DO
472
473 CALL get_reference_point(rpoint, qs_env=qs_env, reference=use_mom_ref_coac, ref_point=ref_point)
474 CALL build_local_moment_matrix(qs_env, matrix_multipole, n_moments, ref_point=rpoint)
475
476 CALL timestop(handle)
477
478 END SUBROUTINE get_multipoles_ao
479
480! **************************************************************************************************
481!> \brief Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b >
482!> Optionally stores the multipole moments themselves for free.
483!> Note that the multipole moments are symmetric while their derivatives are anti-symmetric
484!> Only first derivatives are performed, e. g. x d/dy
485!> \param qs_env ...
486!> \param moments_der will contain the derivatives of the multipole moments
487!> \param nmoments_der order of the moments with derivatives
488!> \param nmoments order of the multipole moments (no derivatives, same output as
489!> build_local_moment_matrix, needs moments as arguments to store results)
490!> \param ref_point ...
491!> \param moments contains the multipole moments, optionally for free, up to order nmoments
492!> \param basis_type basis set used for the integral evaluation
493!> \param minimum_image apply the half-cell image filter
494!> \param ordered use ordered neighbor blocks and, for ordinary derivatives, the full neighbor list
495!> \param lambda optional atom selector for the right-hand derivative
496!> \param deltaR optional weights for left and right nuclear-coordinate derivatives
497!> \param neighbor_image use the neighbor-list image for coordinate construction
498! **************************************************************************************************
499 SUBROUTINE build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, &
500 ref_point, moments, basis_type, minimum_image, ordered, lambda, &
501 deltaR, neighbor_image)
502 TYPE(qs_environment_type), POINTER :: qs_env
503 TYPE(dbcsr_p_type), DIMENSION(:, :), &
504 INTENT(INOUT), POINTER :: moments_der
505 INTEGER, INTENT(IN) :: nmoments_der, nmoments
506 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
507 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
508 OPTIONAL, POINTER :: moments
509 CHARACTER(len=*), OPTIONAL :: basis_type
510 LOGICAL, OPTIONAL :: minimum_image, ordered
511 INTEGER, OPTIONAL :: lambda
512 REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: deltar
513 LOGICAL, OPTIONAL :: neighbor_image
514
515 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moments_der_matrix'
516
517 INTEGER :: dimders, handle, i, iatom, icol, ider, ii, ikind, inode, ipgf, irow, iset, j, &
518 jatom, jkind, jpgf, jset, last_jatom, lda_min, m_dim, maxco, maxsgf, na, nb, ncoa, ncob, nda, ndb, &
519 nders, nkind, nm, nmom_build, nseta, nsetb, sgfa, sgfb
520 LOGICAL :: accumulate_images, found, my_minimum_image, my_neighbor_image, &
521 my_ordered, new_atom_pair
522 REAL(kind=dp) :: dab, derivative_sign
523 REAL(kind=dp), DIMENSION(3) :: half_cell
524 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
525 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
526 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
527 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
528 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab_tmp
529 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mom_block
530 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mom_block_der
531 TYPE(cell_type), POINTER :: cell
532 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
533 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
534 TYPE(neighbor_list_iterator_p_type), &
535 DIMENSION(:), POINTER :: nl_iterator
536 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
537 POINTER :: sab_all, sab_orb
538 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
539 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
540 TYPE(qs_kind_type), POINTER :: qs_kind
541
542 nmom_build = max(nmoments, nmoments_der) ! build moments up to order nmom_buiod
543 IF (nmom_build < 1) RETURN
544
545 CALL timeset(routinen, handle)
546
547 nders = 1 ! only first order derivatives
548 dimders = ncoset(nders) - 1
549
550 my_minimum_image = .false.
551 IF (PRESENT(minimum_image)) my_minimum_image = minimum_image
552 my_ordered = .false.
553 IF (PRESENT(ordered)) my_ordered = ordered
554 my_neighbor_image = my_ordered
555 IF (PRESENT(neighbor_image)) my_neighbor_image = neighbor_image
556 cpassert(.NOT. (PRESENT(deltar) .AND. PRESENT(lambda)))
557 IF (PRESENT(lambda) .OR. PRESENT(deltar)) THEN
558 cpassert(my_ordered)
559 END IF
560 accumulate_images = my_ordered .AND. my_neighbor_image
561 IF (PRESENT(deltar) .OR. PRESENT(lambda)) accumulate_images = .false.
562 derivative_sign = 1._dp
563 ! diff_momop returns the right nuclear-center derivative. The former
564 ! rRc_xyz_der_ao interface exposed the corresponding electronic derivative
565 ! on the ordinary ordered path.
566 IF (my_ordered .AND. .NOT. PRESENT(deltar) .AND. .NOT. PRESENT(lambda)) THEN
567 derivative_sign = -1._dp
568 END IF
569
570 NULLIFY (qs_kind_set, particle_set, sab_all, sab_orb, cell)
571 CALL get_qs_env(qs_env=qs_env, &
572 qs_kind_set=qs_kind_set, &
573 particle_set=particle_set, &
574 cell=cell, &
575 sab_orb=sab_orb, &
576 sab_all=sab_all)
577
578 IF (my_ordered) THEN
579 cpassert(nmoments_der > 0)
580 cpassert(dbcsr_get_matrix_type(moments_der(1, 1)%matrix) == dbcsr_type_no_symmetry)
581 IF (nmoments > 0) THEN
582 cpassert(PRESENT(moments))
583 cpassert(dbcsr_get_matrix_type(moments(1)%matrix) == dbcsr_type_no_symmetry)
584 END IF
585 END IF
586
587 IF (my_minimum_image) THEN
588 half_cell = 0.5_dp*norm2(cell%hmat, dim=1)
589 END IF
590
591 nkind = SIZE(qs_kind_set)
592
593 ! Work storage
594 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
595 maxco=maxco, maxsgf=maxsgf, &
596 basis_type=basis_type)
597
598 IF (nmoments > 0) THEN
599 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
600 cpassert(SIZE(moments) == nm)
601 ! storage for integrals
602 ALLOCATE (mab(maxco, maxco, nm))
603 ! blocks
604 mab(:, :, :) = 0.0_dp
605 ALLOCATE (mom_block(nm))
606 DO i = 1, nm
607 NULLIFY (mom_block(i)%block)
608 END DO
609 END IF
610
611 IF (nmoments_der > 0) THEN
612 m_dim = ncoset(nmoments_der) - 1
613 cpassert(SIZE(moments_der, dim=1) == m_dim)
614 cpassert(SIZE(moments_der, dim=2) == dimders)
615 ! storage for integrals
616 ALLOCATE (difmab(maxco, maxco, m_dim, dimders))
617 difmab(:, :, :, :) = 0.0_dp
618 ! blocks
619 ALLOCATE (mom_block_der(m_dim, dimders))
620 DO i = 1, m_dim
621 DO ider = 1, dimders
622 NULLIFY (mom_block_der(i, ider)%block)
623 END DO
624 END DO
625 END IF
626
627 ALLOCATE (work(maxco, maxsgf))
628 work(:, :) = 0.0_dp
629
630 NULLIFY (basis_set_a, basis_set_b, basis_set_list)
631 NULLIFY (qs_kind)
632 ALLOCATE (basis_set_list(nkind))
633 DO ikind = 1, nkind
634 qs_kind => qs_kind_set(ikind)
635 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
636 IF (ASSOCIATED(basis_set_a)) THEN
637 basis_set_list(ikind)%gto_basis_set => basis_set_a
638 ELSE
639 NULLIFY (basis_set_list(ikind)%gto_basis_set)
640 END IF
641 END DO
642
643 ! Calculate derivatives looping over neighbour list
644 NULLIFY (nl_iterator)
645 IF (my_ordered) THEN
646 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
647 ELSE
648 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
649 END IF
650 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
651 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
652 iatom=iatom, jatom=jatom, r=rab)
653 basis_set_a => basis_set_list(ikind)%gto_basis_set
654 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
655 basis_set_b => basis_set_list(jkind)%gto_basis_set
656 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
657 associate( &
658 ! basis ikind
659 first_sgfa => basis_set_a%first_sgf, &
660 la_max => basis_set_a%lmax, &
661 la_min => basis_set_a%lmin, &
662 npgfa => basis_set_a%npgf, &
663 nsgfa => basis_set_a%nsgf_set, &
664 rpgfa => basis_set_a%pgf_radius, &
665 set_radius_a => basis_set_a%set_radius, &
666 sphi_a => basis_set_a%sphi, &
667 zeta => basis_set_a%zet, &
668 ! basis jkind, &
669 first_sgfb => basis_set_b%first_sgf, &
670 lb_max => basis_set_b%lmax, &
671 lb_min => basis_set_b%lmin, &
672 npgfb => basis_set_b%npgf, &
673 nsgfb => basis_set_b%nsgf_set, &
674 rpgfb => basis_set_b%pgf_radius, &
675 set_radius_b => basis_set_b%set_radius, &
676 sphi_b => basis_set_b%sphi, &
677 zetb => basis_set_b%zet)
678
679 IF (SIZE(sphi_a, 1) == 0) cycle
680 IF (SIZE(sphi_b, 1) == 0) cycle
681
682 nseta = basis_set_a%nset
683 nsetb = basis_set_b%nset
684
685 IF (inode == 1) last_jatom = 0
686
687 IF (my_minimum_image) THEN
688 IF (any(abs(rab(:)) > half_cell(:))) cycle
689 END IF
690
691 new_atom_pair = jatom /= last_jatom
692 ! The ordinary ordered neighbor-image path accumulates all images for one atom pair.
693 ! Specialized derivative paths select one representative image, as in the former builders.
694 IF (.NOT. accumulate_images .AND. .NOT. new_atom_pair) cycle
695
696 IF (new_atom_pair) THEN
697 last_jatom = jatom
698
699 ! get blocks
700 IF (my_ordered) THEN
701 irow = iatom
702 icol = jatom
703 ELSE IF (iatom <= jatom) THEN
704 irow = iatom
705 icol = jatom
706 ELSE
707 irow = jatom
708 icol = iatom
709 END IF
710
711 IF (nmoments > 0) THEN
712 DO i = 1, nm
713 NULLIFY (mom_block(i)%block)
714 ! get block from pre calculated overlap matrix
715 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
716 row=irow, col=icol, block=mom_block(i)%block, found=found)
717 cpassert(found .AND. ASSOCIATED(mom_block(i)%block))
718 mom_block(i)%block = 0._dp
719 END DO
720 END IF
721 IF (nmoments_der > 0) THEN
722 DO i = 1, m_dim
723 DO ider = 1, dimders
724 NULLIFY (mom_block_der(i, ider)%block)
725 CALL dbcsr_get_block_p(matrix=moments_der(i, ider)%matrix, &
726 row=irow, col=icol, &
727 block=mom_block_der(i, ider)%block, &
728 found=found)
729 cpassert(found .AND. ASSOCIATED(mom_block_der(i, ider)%block))
730 mom_block_der(i, ider)%block = 0._dp
731 END DO
732 END DO
733 END IF
734 END IF
735
736 ! reference point
737 IF (PRESENT(ref_point)) THEN
738 rc(:) = ref_point(:)
739 ELSE
740 rc(:) = 0._dp
741 END IF
742 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
743 ! by folding around the center, such screwing can be avoided for a proper choice of center.
744 IF (my_neighbor_image) THEN
745 ra(:) = pbc(particle_set(iatom)%r(:), cell)
746 rb(:) = ra(:) + rab(:)
747 rac(:) = pbc(rc, ra, cell)
748 rbc(:) = rac(:) + rab(:)
749 ELSE IF (my_ordered) THEN
750 ra(:) = particle_set(iatom)%r(:)
751 rb(:) = particle_set(jatom)%r(:)
752 rab(:) = pbc(rb, ra, cell)
753 rac(:) = pbc(ra - rc, cell)
754 rbc(:) = pbc(rb - rc, cell)
755 ELSE
756 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
757 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
758 ! we dont use PBC at this point
759 rab(:) = ra(:) - rb(:)
760 rac(:) = ra(:) - rc(:)
761 rbc(:) = rb(:) - rc(:)
762 END IF
763 dab = norm2(rab)
764
765 DO iset = 1, nseta
766
767 ncoa = npgfa(iset)*ncoset(la_max(iset))
768 sgfa = first_sgfa(1, iset)
769
770 DO jset = 1, nsetb
771
772 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
773
774 ncob = npgfb(jset)*ncoset(lb_max(jset))
775 sgfb = first_sgfb(1, jset)
776
777 NULLIFY (mab_tmp)
778 ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
779 npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(nmom_build) - 1))
780
781 ! Calculate the primitive integrals (need l+1 for derivatives)
782 lda_min = max(0, la_min(iset) - 1)
783 CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), &
784 rpgfa(:, iset), lda_min, &
785 lb_max(jset) + 1, npgfb(jset), zetb(:, jset), &
786 rpgfb(:, jset), nmom_build, rac, rbc, mab_tmp)
787
788 IF (nmoments_der > 0) THEN
789 CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), &
790 rpgfa(:, iset), la_min(iset), &
791 lb_max(jset), npgfb(jset), zetb(:, jset), &
792 rpgfb(:, jset), lb_min(jset), &
793 nmoments_der, rac, rbc, difmab, mab_ext=mab_tmp, &
794 lambda=lambda, deltar=deltar, iatom=iatom, jatom=jatom)
795 END IF
796
797 IF (nmoments > 0) THEN
798 ! copy subset of mab_tmp (l+1) to mab (l)
799 mab = 0.0_dp
800 DO ii = 1, nm
801 na = 0
802 nda = 0
803 DO ipgf = 1, npgfa(iset)
804 nb = 0
805 ndb = 0
806 DO jpgf = 1, npgfb(jset)
807 DO j = 1, ncoset(lb_max(jset))
808 DO i = 1, ncoset(la_max(iset))
809 mab(i + na, j + nb, ii) = mab_tmp(i + nda, j + ndb, ii)
810 END DO ! i
811 END DO ! j
812 nb = nb + ncoset(lb_max(jset))
813 ndb = ndb + ncoset(lb_max(jset) + 1)
814 END DO ! jpgf
815 na = na + ncoset(la_max(iset))
816 nda = nda + ncoset(la_max(iset) + 1)
817 END DO ! ipgf
818 END DO
819 ! Contraction step
820 DO i = 1, nm
821
822 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
823 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
824 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
825 0.0_dp, work(1, 1), SIZE(work, 1))
826
827 IF (my_ordered .OR. iatom <= jatom) THEN
828 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
829 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
830 work(1, 1), SIZE(work, 1), &
831 1.0_dp, mom_block(i)%block(sgfa, sgfb), &
832 SIZE(mom_block(i)%block, 1))
833 ELSE
834 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
835 1.0_dp, work(1, 1), SIZE(work, 1), &
836 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
837 1.0_dp, mom_block(i)%block(sgfb, sgfa), &
838 SIZE(mom_block(i)%block, 1))
839 END IF
840 END DO
841 END IF
842
843 IF (nmoments_der > 0) THEN
844 DO i = 1, m_dim
845 DO ider = 1, dimders
846 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
847 derivative_sign, difmab(1, 1, i, ider), SIZE(difmab, 1), &
848 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
849 0._dp, work(1, 1), SIZE(work, 1))
850
851 IF (my_ordered .OR. iatom <= jatom) THEN
852 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
853 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
854 work(1, 1), SIZE(work, 1), &
855 1.0_dp, mom_block_der(i, ider)%block(sgfa, sgfb), &
856 SIZE(mom_block_der(i, ider)%block, 1))
857 ELSE
858 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
859 -1.0_dp, work(1, 1), SIZE(work, 1), &
860 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
861 1.0_dp, mom_block_der(i, ider)%block(sgfb, sgfa), &
862 SIZE(mom_block_der(i, ider)%block, 1))
863 END IF
864 END DO
865 END DO
866 END IF
867 DEALLOCATE (mab_tmp)
868 END DO
869 END DO
870 END associate
871 END DO
872 CALL neighbor_list_iterator_release(nl_iterator)
873
874 ! deallocations
875 DEALLOCATE (basis_set_list)
876 DEALLOCATE (work)
877 IF (nmoments > 0) THEN
878 DEALLOCATE (mab)
879 DO i = 1, nm
880 NULLIFY (mom_block(i)%block)
881 END DO
882 DEALLOCATE (mom_block)
883 END IF
884 IF (nmoments_der > 0) THEN
885 DEALLOCATE (difmab)
886 DO i = 1, m_dim
887 DO ider = 1, dimders
888 NULLIFY (mom_block_der(i, ider)%block)
889 END DO
890 END DO
891 DEALLOCATE (mom_block_der)
892 END IF
893
894 CALL timestop(handle)
895
896 END SUBROUTINE build_local_moments_der_matrix
897
898! **************************************************************************************************
899!> \brief ...
900!> \param qs_env ...
901!> \param magmom ...
902!> \param nmoments ...
903!> \param ref_point ...
904!> \param ref_points ...
905!> \param basis_type ...
906! **************************************************************************************************
907 SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
908
909 TYPE(qs_environment_type), POINTER :: qs_env
910 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom
911 INTEGER, INTENT(IN) :: nmoments
912 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
913 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
914 OPTIONAL :: ref_points
915 CHARACTER(len=*), OPTIONAL :: basis_type
916
917 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_magmom_matrix'
918
919 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, maxco, &
920 maxsgf, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
921 LOGICAL :: found
922 REAL(kind=dp) :: dab
923 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
924 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
925 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
926 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
927 TYPE(cell_type), POINTER :: cell
928 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
929 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
930 TYPE(neighbor_list_iterator_p_type), &
931 DIMENSION(:), POINTER :: nl_iterator
932 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
933 POINTER :: sab_orb
934 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
935 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
936 TYPE(qs_kind_type), POINTER :: qs_kind
937
938 IF (nmoments < 1) RETURN
939
940 CALL timeset(routinen, handle)
941
942 ! magnetic dipoles/angular moments only
943 nm = 3
944
945 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
946 CALL get_qs_env(qs_env=qs_env, &
947 qs_kind_set=qs_kind_set, &
948 particle_set=particle_set, cell=cell, &
949 sab_orb=sab_orb)
950
951 nkind = SIZE(qs_kind_set)
952
953 ! Allocate work storage
954 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
955 maxco=maxco, maxsgf=maxsgf)
956
957 ALLOCATE (mab(maxco, maxco, nm))
958 mab(:, :, :) = 0.0_dp
959
960 ALLOCATE (work(maxco, maxsgf))
961 work(:, :) = 0.0_dp
962
963 ALLOCATE (mint(nm))
964 DO i = 1, nm
965 NULLIFY (mint(i)%block)
966 END DO
967
968 ALLOCATE (basis_set_list(nkind))
969 DO ikind = 1, nkind
970 qs_kind => qs_kind_set(ikind)
971 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
972 IF (ASSOCIATED(basis_set_a)) THEN
973 basis_set_list(ikind)%gto_basis_set => basis_set_a
974 ELSE
975 NULLIFY (basis_set_list(ikind)%gto_basis_set)
976 END IF
977 END DO
978 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
979 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
980 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
981 iatom=iatom, jatom=jatom, r=rab)
982 basis_set_a => basis_set_list(ikind)%gto_basis_set
983 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
984 basis_set_b => basis_set_list(jkind)%gto_basis_set
985 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
986 associate( &
987 ! basis ikind
988 first_sgfa => basis_set_a%first_sgf, &
989 la_max => basis_set_a%lmax, &
990 la_min => basis_set_a%lmin, &
991 npgfa => basis_set_a%npgf, &
992 nsgfa => basis_set_a%nsgf_set, &
993 rpgfa => basis_set_a%pgf_radius, &
994 set_radius_a => basis_set_a%set_radius, &
995 sphi_a => basis_set_a%sphi, &
996 zeta => basis_set_a%zet, &
997 ! basis jkind, &
998 first_sgfb => basis_set_b%first_sgf, &
999 lb_max => basis_set_b%lmax, &
1000 lb_min => basis_set_b%lmin, &
1001 npgfb => basis_set_b%npgf, &
1002 nsgfb => basis_set_b%nsgf_set, &
1003 rpgfb => basis_set_b%pgf_radius, &
1004 set_radius_b => basis_set_b%set_radius, &
1005 sphi_b => basis_set_b%sphi, &
1006 zetb => basis_set_b%zet)
1007
1008 nseta = basis_set_a%nset
1009 nsetb = basis_set_b%nset
1010
1011 IF (iatom <= jatom) THEN
1012 irow = iatom
1013 icol = jatom
1014 ELSE
1015 irow = jatom
1016 icol = iatom
1017 END IF
1018
1019 DO i = 1, nm
1020 NULLIFY (mint(i)%block)
1021 CALL dbcsr_get_block_p(matrix=magmom(i)%matrix, &
1022 row=irow, col=icol, block=mint(i)%block, found=found)
1023 cpassert(found .AND. ASSOCIATED(mint(i)%block))
1024 mint(i)%block = 0._dp
1025 END DO
1026
1027 ! fold atomic position back into unit cell
1028 IF (PRESENT(ref_points)) THEN
1029 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
1030 ELSE IF (PRESENT(ref_point)) THEN
1031 rc(:) = ref_point(:)
1032 ELSE
1033 rc(:) = 0._dp
1034 END IF
1035 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
1036 ! by folding around the center, such screwing can be avoided for a proper choice of center.
1037 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
1038 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
1039 ! we dont use PBC at this point
1040 rab(:) = ra(:) - rb(:)
1041 rac(:) = ra(:) - rc(:)
1042 rbc(:) = rb(:) - rc(:)
1043 dab = norm2(rab)
1044
1045 DO iset = 1, nseta
1046
1047 ncoa = npgfa(iset)*ncoset(la_max(iset))
1048 sgfa = first_sgfa(1, iset)
1049
1050 DO jset = 1, nsetb
1051
1052 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1053
1054 ncob = npgfb(jset)*ncoset(lb_max(jset))
1055 sgfb = first_sgfb(1, jset)
1056
1057 ! Calculate the primitive integrals
1058 CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), &
1059 rpgfa(:, iset), la_min(iset), &
1060 lb_max(jset), npgfb(jset), zetb(:, jset), &
1061 rpgfb(:, jset), rac, rbc, mab)
1062
1063 ! Contraction step
1064 DO i = 1, nm
1065 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1066 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
1067 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1068 0.0_dp, work(1, 1), SIZE(work, 1))
1069
1070 IF (iatom <= jatom) THEN
1071 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1072 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1073 work(1, 1), SIZE(work, 1), &
1074 1.0_dp, mint(i)%block(sgfa, sgfb), &
1075 SIZE(mint(i)%block, 1))
1076 ELSE
1077 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1078 -1.0_dp, work(1, 1), SIZE(work, 1), &
1079 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1080 1.0_dp, mint(i)%block(sgfb, sgfa), &
1081 SIZE(mint(i)%block, 1))
1082 END IF
1083
1084 END DO
1085
1086 END DO
1087 END DO
1088 END associate
1089 END DO
1090 CALL neighbor_list_iterator_release(nl_iterator)
1091
1092 ! Release work storage
1093 DEALLOCATE (mab, basis_set_list)
1094 DEALLOCATE (work)
1095 DO i = 1, nm
1096 NULLIFY (mint(i)%block)
1097 END DO
1098 DEALLOCATE (mint)
1099
1100 CALL timestop(handle)
1101
1102 END SUBROUTINE build_local_magmom_matrix
1103
1104! **************************************************************************************************
1105!> \brief ...
1106!> \param qs_env ...
1107!> \param cosmat ...
1108!> \param sinmat ...
1109!> \param kvec ...
1110!> \param sab_orb_external ...
1111!> \param basis_type ...
1112! **************************************************************************************************
1113 SUBROUTINE build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
1114
1115 TYPE(qs_environment_type), POINTER :: qs_env
1116 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1117 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1118 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1119 OPTIONAL, POINTER :: sab_orb_external
1120 CHARACTER(len=*), OPTIONAL :: basis_type
1121
1122 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_moment_matrix'
1123
1124 INTEGER :: handle
1125 TYPE(dbcsr_p_type), DIMENSION(2, 1) :: op_sm_set
1126 REAL(kind=dp), DIMENSION(3, 1) :: kvec_batch
1127
1128 CALL timeset(routinen, handle)
1129
1130 op_sm_set(1, 1)%matrix => cosmat
1131 op_sm_set(2, 1)%matrix => sinmat
1132 kvec_batch(:, 1) = kvec(:)
1133
1134 CALL build_exp_ikr_matrix(qs_env, op_sm_set, kvec_batch, &
1135 sab_orb_external=sab_orb_external, basis_type=basis_type)
1136
1137 CALL timestop(handle)
1138
1139 END SUBROUTINE build_berry_moment_matrix
1140
1141! **************************************************************************************************
1142!> \brief ...
1143!> \param qs_env ...
1144!> \param cosmat ...
1145!> \param sinmat ...
1146!> \param kvec ...
1147!> \param basis_type ...
1148!> \param ordered use all ordered atom pairs for a non-Hermitian cross-k overlap
1149!> \param ordered_kpoints image-cell mapping for the ordered neighbour list
1150! **************************************************************************************************
1151 SUBROUTINE build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type, ordered, ordered_kpoints)
1152
1153 TYPE(qs_environment_type), POINTER :: qs_env
1154 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: cosmat, sinmat
1155 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1156 CHARACTER(len=*), OPTIONAL :: basis_type
1157 LOGICAL, INTENT(IN), OPTIONAL :: ordered
1158 TYPE(kpoint_type), POINTER, OPTIONAL :: ordered_kpoints
1159
1160 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_kpoint_matrix'
1161
1162 CHARACTER :: matrix_symmetry
1163 COMPLEX(dp), ALLOCATABLE :: expab(:, :)
1164 INTEGER :: handle, i, iatom, ic, icol, ikind, inode, irow, iset, jatom, jkind, jset, &
1165 ldsa, ldsb, ldwork, natom, ncoa, ncob, nimg, nkind, nseta, nsetb, sgfa, sgfb
1166 INTEGER, DIMENSION(3) :: icell
1167 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
1168 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1169 LOGICAL :: found, full_pairs, stable_images, use_cell_mapping
1170 REAL(dp), DIMENSION(:, :), POINTER :: cblock, sblock, work
1171 REAL(kind=dp) :: dab
1172 REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
1173 TYPE(cell_type), POINTER :: cell
1174 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1175 TYPE(dft_control_type), POINTER :: dft_control
1176 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1177 TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
1178 TYPE(kpoint_type), POINTER :: kpoints
1179 TYPE(neighbor_list_iterator_p_type), &
1180 DIMENSION(:), POINTER :: nl_iterator
1181 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1182 POINTER :: sab_orb
1183 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1184 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1185 TYPE(qs_kind_type), POINTER :: qs_kind
1186 TYPE(qs_ks_env_type), POINTER :: ks_env
1187
1188 CALL timeset(routinen, handle)
1189
1190 full_pairs = .false.
1191 IF (PRESENT(ordered)) full_pairs = ordered
1192 matrix_symmetry = dbcsr_type_symmetric
1193 IF (full_pairs) matrix_symmetry = dbcsr_type_no_symmetry
1194
1195 CALL get_qs_env(qs_env, &
1196 ks_env=ks_env, &
1197 dft_control=dft_control)
1198 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1199 stable_images = kpoints%nkp > 0 .AND. kpoints%symmetry
1200 nimg = dft_control%nimages
1201 IF (nimg > 1) THEN
1202 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1203 use_cell_mapping = .true.
1204 ELSE
1205 use_cell_mapping = .false.
1206 END IF
1207
1208 CALL get_qs_env(qs_env=qs_env, &
1209 qs_kind_set=qs_kind_set, &
1210 particle_set=particle_set, cell=cell, &
1211 sab_orb=sab_orb)
1212
1213 IF (full_pairs) THEN
1214 CALL get_qs_env(qs_env, sab_all=sab_orb)
1215 IF (.NOT. ASSOCIATED(sab_orb)) cpabort("Ordered Berry integrals require k-point neighbour lists.")
1216 cpassert(PRESENT(ordered_kpoints))
1217 CALL get_kpoint_info(ordered_kpoints, cell_to_index=cell_to_index)
1218 nimg = SIZE(ordered_kpoints%index_to_cell, 2)
1219 use_cell_mapping = .true.
1220 END IF
1221
1222 nkind = SIZE(qs_kind_set)
1223 natom = SIZE(particle_set)
1224 ALLOCATE (basis_set_list(nkind))
1225 DO ikind = 1, nkind
1226 qs_kind => qs_kind_set(ikind)
1227 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=basis_type)
1228 IF (ASSOCIATED(basis_set)) THEN
1229 basis_set_list(ikind)%gto_basis_set => basis_set
1230 ELSE
1231 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1232 END IF
1233 END DO
1234
1235 ALLOCATE (row_blk_sizes(natom))
1236 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1237 basis=basis_set_list)
1238 CALL get_ks_env(ks_env, dbcsr_dist=dbcsr_dist)
1239 ! (re)allocate matrix sets
1240 CALL dbcsr_allocate_matrix_set(sinmat, 1, nimg)
1241 CALL dbcsr_allocate_matrix_set(cosmat, 1, nimg)
1242 DO i = 1, nimg
1243 ! sin
1244 ALLOCATE (sinmat(1, i)%matrix)
1245 CALL dbcsr_create(matrix=sinmat(1, i)%matrix, &
1246 name="SINMAT", &
1247 dist=dbcsr_dist, matrix_type=matrix_symmetry, &
1248 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
1249 CALL cp_dbcsr_alloc_block_from_nbl(sinmat(1, i)%matrix, sab_orb)
1250 CALL dbcsr_set(sinmat(1, i)%matrix, 0.0_dp)
1251 ! cos
1252 ALLOCATE (cosmat(1, i)%matrix)
1253 CALL dbcsr_create(matrix=cosmat(1, i)%matrix, &
1254 name="COSMAT", &
1255 dist=dbcsr_dist, matrix_type=matrix_symmetry, &
1256 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
1257 CALL cp_dbcsr_alloc_block_from_nbl(cosmat(1, i)%matrix, sab_orb)
1258 CALL dbcsr_set(cosmat(1, i)%matrix, 0.0_dp)
1259 END DO
1260
1261 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
1262 ALLOCATE (expab(ldwork, ldwork))
1263 ALLOCATE (work(ldwork, ldwork))
1264
1265 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1266 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1267 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1268 iatom=iatom, jatom=jatom, r=rab, cell=icell)
1269 basis_set_a => basis_set_list(ikind)%gto_basis_set
1270 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1271 basis_set_b => basis_set_list(jkind)%gto_basis_set
1272 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1273 associate( &
1274 ! basis ikind
1275 first_sgfa => basis_set_a%first_sgf, &
1276 la_max => basis_set_a%lmax, &
1277 la_min => basis_set_a%lmin, &
1278 npgfa => basis_set_a%npgf, &
1279 nsgfa => basis_set_a%nsgf_set, &
1280 rpgfa => basis_set_a%pgf_radius, &
1281 set_radius_a => basis_set_a%set_radius, &
1282 sphi_a => basis_set_a%sphi, &
1283 zeta => basis_set_a%zet, &
1284 ! basis jkind, &
1285 first_sgfb => basis_set_b%first_sgf, &
1286 lb_max => basis_set_b%lmax, &
1287 lb_min => basis_set_b%lmin, &
1288 npgfb => basis_set_b%npgf, &
1289 nsgfb => basis_set_b%nsgf_set, &
1290 rpgfb => basis_set_b%pgf_radius, &
1291 set_radius_b => basis_set_b%set_radius, &
1292 sphi_b => basis_set_b%sphi, &
1293 zetb => basis_set_b%zet)
1294
1295 nseta = basis_set_a%nset
1296 nsetb = basis_set_b%nset
1297
1298 ldsa = SIZE(sphi_a, 1)
1299 ldsb = SIZE(sphi_b, 1)
1300
1301 IF (iatom <= jatom .OR. full_pairs) THEN
1302 irow = iatom
1303 icol = jatom
1304 ELSE
1305 irow = jatom
1306 icol = iatom
1307 END IF
1308
1309 IF (use_cell_mapping) THEN
1310 ic = cell_to_index(icell(1), icell(2), icell(3))
1311 cpassert(ic > 0)
1312 ELSE
1313 ic = 1
1314 END IF
1315
1316 NULLIFY (sblock)
1317 CALL dbcsr_get_block_p(matrix=sinmat(1, ic)%matrix, &
1318 row=irow, col=icol, block=sblock, found=found)
1319 cpassert(found)
1320 NULLIFY (cblock)
1321 CALL dbcsr_get_block_p(matrix=cosmat(1, ic)%matrix, &
1322 row=irow, col=icol, block=cblock, found=found)
1323 cpassert(found)
1324
1325 ! Absolute phases must use the same atom images as the orbital neighbour list.
1326 IF (stable_images) THEN
1327 ra(:) = pbc_stable(particle_set(iatom)%r(:), cell)
1328 ELSE
1329 ra(:) = pbc(particle_set(iatom)%r(:), cell)
1330 END IF
1331 rb(:) = ra + rab
1332 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1333
1334 DO iset = 1, nseta
1335
1336 ncoa = npgfa(iset)*ncoset(la_max(iset))
1337 sgfa = first_sgfa(1, iset)
1338
1339 DO jset = 1, nsetb
1340
1341 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1342
1343 ncob = npgfb(jset)*ncoset(lb_max(jset))
1344 sgfb = first_sgfb(1, jset)
1345
1346 ! Calculate the primitive integrals
1347 CALL expikr(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1348 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1349 ra, rb, kvec, expab)
1350 CALL contract_expikr(cblock, sblock, &
1351 iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1352 jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1353 expab, work, ldwork, ordered=full_pairs)
1354
1355 END DO
1356 END DO
1357 END associate
1358 END DO
1359 CALL neighbor_list_iterator_release(nl_iterator)
1360
1361 DEALLOCATE (expab)
1362 DEALLOCATE (work)
1363 DEALLOCATE (basis_set_list)
1364 DEALLOCATE (row_blk_sizes)
1365
1366 CALL timestop(handle)
1367
1368 END SUBROUTINE build_berry_kpoint_matrix
1369
1370! **************************************************************************************************
1371!> \brief ...
1372!> \param qs_env ...
1373!> \param magnetic ...
1374!> \param nmoments ...
1375!> \param reference ...
1376!> \param ref_point ...
1377!> \param unit_number ...
1378!> \param unit_str unit of the printed moments (dipole unit of cp_units), default debye
1379! **************************************************************************************************
1380 SUBROUTINE qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number, unit_str)
1381
1382 TYPE(qs_environment_type), POINTER :: qs_env
1383 LOGICAL, INTENT(IN) :: magnetic
1384 INTEGER, INTENT(IN) :: nmoments, reference
1385 REAL(dp), DIMENSION(:), POINTER :: ref_point
1386 INTEGER, INTENT(IN) :: unit_number
1387 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: unit_str
1388
1389 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_berry_phase'
1390 REAL(dp), PARAMETER :: eps_occ = 1.0e-6_dp
1391
1392 CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
1393 CHARACTER(LEN=default_string_length) :: description
1394 COMPLEX(dp) :: xphase(3), zdet, zdeta, zi(3), &
1395 zij(3, 3), zijk(3, 3, 3), &
1396 zijkl(3, 3, 3, 3), zphase(3), zz
1397 INTEGER :: handle, i, ia, idim, ikind, ispin, ix, &
1398 iy, iz, j, k, l, nao, nm, nmo, nmo_max, &
1399 nmom, tmp_dim
1400 LOGICAL :: floating, fractional, ghost
1401 LOGICAL, ALLOCATABLE, DIMENSION(:) :: all_occupied
1402 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: mo_occupied
1403 REAL(dp) :: charge, ci(3), cij(3, 3), dd, nel_det, &
1404 occ, trace
1405 REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom
1406 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
1407 REAL(dp), DIMENSION(3) :: kvec, qq, rcc, ria
1408 REAL(dp), DIMENSION(:), POINTER :: occupation_numbers
1409 TYPE(atomic_kind_type), POINTER :: atomic_kind
1410 TYPE(cell_type), POINTER :: cell
1411 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat
1412 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
1413 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: opvec
1414 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: op_fm_set
1415 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
1416 TYPE(cp_fm_type), POINTER :: mo_coeff
1417 TYPE(cp_result_type), POINTER :: results
1418 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
1419 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1420 TYPE(dft_control_type), POINTER :: dft_control
1421 TYPE(distribution_1d_type), POINTER :: local_particles
1422 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1423 TYPE(mp_para_env_type), POINTER :: para_env
1424 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1425 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1426 TYPE(qs_rho_type), POINTER :: rho
1427 TYPE(rt_prop_type), POINTER :: rtp
1428
1429 cpassert(ASSOCIATED(qs_env))
1430
1431 IF (ASSOCIATED(qs_env%ls_scf_env)) THEN
1432 IF (unit_number > 0) WRITE (unit_number, *) "Periodic moment calculation not implemented in linear scaling code"
1433 RETURN
1434 END IF
1435
1436 CALL timeset(routinen, handle)
1437
1438 ! restrict maximum moment available
1439 nmom = min(nmoments, 2)
1440
1441 nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
1442 ! rmom(:,1)=electronic
1443 ! rmom(:,2)=nuclear
1444 ! rmom(:,1)=total
1445 ALLOCATE (rmom(nm + 1, 3))
1446 ALLOCATE (rlab(nm + 1))
1447 rmom = 0.0_dp
1448 rlab = ""
1449 IF (magnetic) THEN
1450 nm = 3
1451 ALLOCATE (mmom(nm))
1452 mmom = 0._dp
1453 END IF
1454
1455 NULLIFY (dft_control, rho, cell, particle_set, results, para_env, &
1456 local_particles, matrix_s, mos, rho_ao)
1457
1458 CALL get_qs_env(qs_env, &
1459 dft_control=dft_control, &
1460 rho=rho, &
1461 cell=cell, &
1462 results=results, &
1463 particle_set=particle_set, &
1464 qs_kind_set=qs_kind_set, &
1465 para_env=para_env, &
1466 local_particles=local_particles, &
1467 matrix_s=matrix_s, &
1468 mos=mos)
1469
1470 CALL qs_rho_get(rho, rho_ao=rho_ao)
1471
1472 NULLIFY (cosmat, sinmat)
1473 ALLOCATE (cosmat, sinmat)
1474 CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
1475 CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
1476
1477 ALLOCATE (op_fm_set(2, dft_control%nspins))
1478 ALLOCATE (opvec(dft_control%nspins))
1479 ALLOCATE (eigrmat(dft_control%nspins))
1480 nmo_max = 0
1481 DO ispin = 1, dft_control%nspins
1482 NULLIFY (tmp_fm_struct, mo_coeff)
1483 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
1484 nmo_max = max(nmo_max, nmo)
1485 CALL cp_fm_create(opvec(ispin), mo_coeff%matrix_struct)
1486 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
1487 ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
1488 DO i = 1, SIZE(op_fm_set, 1)
1489 CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
1490 END DO
1491 CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
1492 CALL cp_fm_struct_release(tmp_fm_struct)
1493 END DO
1494
1495 ! occupation
1496 ALLOCATE (mo_occupied(nmo_max, dft_control%nspins), all_occupied(dft_control%nspins))
1497 mo_occupied = .false.
1498 fractional = .false.
1499 nel_det = 0.0_dp
1500 DO ispin = 1, dft_control%nspins
1501 NULLIFY (occupation_numbers)
1502 CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, nmo=nmo, occupation_numbers=occupation_numbers)
1503 DO i = 1, nmo
1504 IF (occupation_numbers(i) > eps_occ) THEN
1505 mo_occupied(i, ispin) = .true.
1506 nel_det = nel_det + occ
1507 IF (abs(occupation_numbers(i) - occ) > eps_occ) fractional = .true.
1508 END IF
1509 END DO
1510 all_occupied(ispin) = all(mo_occupied(1:nmo, ispin))
1511 END DO
1512 IF (fractional) THEN
1513 cpwarn("Berry phase moments for fractional MOs' occupation numbers not implemented")
1514 END IF
1515
1516 ! reference point
1517 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
1518 rcc = pbc(rcc, cell)
1519
1520 ! label
1521 DO l = 1, nm
1522 ix = indco(1, l + 1)
1523 iy = indco(2, l + 1)
1524 iz = indco(3, l + 1)
1525 CALL set_label(rlab(l + 1), ix, iy, iz)
1526 END DO
1527
1528 ! nuclear contribution
1529 DO ia = 1, SIZE(particle_set)
1530 atomic_kind => particle_set(ia)%atomic_kind
1531 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1532 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1533 IF (.NOT. ghost .AND. .NOT. floating) THEN
1534 rmom(1, 2) = rmom(1, 2) - charge
1535 END IF
1536 END DO
1537 ria = twopi*matmul(cell%h_inv, rcc)
1538 zphase = cmplx(cos(ria), sin(ria), dp)**rmom(1, 2)
1539
1540 zi = 0._dp
1541 zij = 0._dp
1542 zijk = 0._dp
1543 zijkl = 0._dp
1544
1545 DO l = 1, nmom
1546 SELECT CASE (l)
1547 CASE (1)
1548 ! Dipole
1549 zi(:) = cmplx(1._dp, 0._dp, dp)
1550 DO ia = 1, SIZE(particle_set)
1551 atomic_kind => particle_set(ia)%atomic_kind
1552 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1553 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1554 IF (.NOT. ghost .AND. .NOT. floating) THEN
1555 ria = particle_set(ia)%r
1556 ria = pbc(ria, cell)
1557 DO i = 1, 3
1558 kvec(:) = twopi*cell%h_inv(i, :)
1559 dd = sum(kvec(:)*ria(:))
1560 zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
1561 zi(i) = zi(i)*zdeta
1562 END DO
1563 END IF
1564 END DO
1565 zi = zi*zphase
1566 ci = aimag(log(zi))/twopi
1567 qq = aimag(log(zi))
1568 rmom(2:4, 2) = matmul(cell%hmat, ci)
1569 CASE (2)
1570 ! Quadrupole
1571 cpabort("Berry phase moments bigger than 1 not implemented")
1572 zij(:, :) = cmplx(1._dp, 0._dp, dp)
1573 DO ia = 1, SIZE(particle_set)
1574 atomic_kind => particle_set(ia)%atomic_kind
1575 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1576 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
1577 ria = particle_set(ia)%r
1578 ria = pbc(ria, cell)
1579 DO i = 1, 3
1580 DO j = i, 3
1581 kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1582 dd = sum(kvec(:)*ria(:))
1583 zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
1584 zij(i, j) = zij(i, j)*zdeta
1585 zij(j, i) = zij(i, j)
1586 END DO
1587 END DO
1588 END DO
1589 DO i = 1, 3
1590 DO j = 1, 3
1591 zij(i, j) = zij(i, j)*zphase(i)*zphase(j)
1592 zz = zij(i, j)/zi(i)/zi(j)
1593 cij(i, j) = aimag(log(zz))/twopi
1594 END DO
1595 END DO
1596 cij = 0.5_dp*cij/twopi/twopi
1597 cij = matmul(matmul(cell%hmat, cij), transpose(cell%hmat))
1598 DO k = 4, 9
1599 ix = indco(1, k + 1)
1600 iy = indco(2, k + 1)
1601 iz = indco(3, k + 1)
1602 IF (ix == 0) THEN
1603 rmom(k + 1, 2) = cij(iy, iz)
1604 ELSE IF (iy == 0) THEN
1605 rmom(k + 1, 2) = cij(ix, iz)
1606 ELSE IF (iz == 0) THEN
1607 rmom(k + 1, 2) = cij(ix, iy)
1608 END IF
1609 END DO
1610 CASE (3)
1611 ! Octapole
1612 cpabort("Berry phase moments bigger than 2 not implemented")
1613 CASE (4)
1614 ! Hexadecapole
1615 cpabort("Berry phase moments bigger than 3 not implemented")
1616 CASE DEFAULT
1617 cpabort("Berry phase moments bigger than 4 not implemented")
1618 END SELECT
1619 END DO
1620
1621 ! electronic contribution
1622
1623 ria = twopi*nel_det*matmul(cell%h_inv, rcc)
1624 xphase = cmplx(cos(ria), sin(ria), dp)
1625
1626 ! charge
1627 trace = 0.0_dp
1628 DO ispin = 1, dft_control%nspins
1629 CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
1630 rmom(1, 1) = rmom(1, 1) + trace
1631 END DO
1632
1633 zi = 0._dp
1634 zij = 0._dp
1635 zijk = 0._dp
1636 zijkl = 0._dp
1637
1638 DO l = 1, nmom
1639 SELECT CASE (l)
1640 CASE (1)
1641 ! Dipole
1642 DO i = 1, 3
1643 kvec(:) = twopi*cell%h_inv(i, :)
1644 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1645 IF (qs_env%run_rtp) THEN
1646 CALL get_qs_env(qs_env, rtp=rtp)
1647 CALL get_rtp(rtp, mos_new=mos_new)
1648 CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1649 ELSE
1650 CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1651 END IF
1652 zdet = cmplx(1._dp, 0._dp, dp)
1653 DO ispin = 1, dft_control%nspins
1654 CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
1655 DO idim = 1, tmp_dim
1656 eigrmat(ispin)%local_data(:, idim) = &
1657 cmplx(op_fm_set(1, ispin)%local_data(:, idim), &
1658 -op_fm_set(2, ispin)%local_data(:, idim), dp)
1659 END DO
1660 IF (.NOT. all_occupied(ispin)) THEN
1661 CALL berry_phase_remove_unocc(eigrmat(ispin), mo_occupied(:, ispin))
1662 END IF
1663 ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
1664 CALL cp_cfm_det(eigrmat(ispin), zdeta)
1665 zdet = zdet*zdeta
1666 IF (dft_control%nspins == 1) THEN
1667 zdet = zdet*zdeta
1668 END IF
1669 END DO
1670 zi(i) = zdet
1671 END DO
1672 zi = zi*xphase
1673 CASE (2)
1674 ! Quadrupole
1675 cpabort("Berry phase moments bigger than 1 not implemented")
1676 DO i = 1, 3
1677 DO j = i, 3
1678 kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1679 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1680 IF (qs_env%run_rtp) THEN
1681 CALL get_qs_env(qs_env, rtp=rtp)
1682 CALL get_rtp(rtp, mos_new=mos_new)
1683 CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1684 ELSE
1685 CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1686 END IF
1687 zdet = cmplx(1._dp, 0._dp, dp)
1688 DO ispin = 1, dft_control%nspins
1689 CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
1690 DO idim = 1, tmp_dim
1691 eigrmat(ispin)%local_data(:, idim) = &
1692 cmplx(op_fm_set(1, ispin)%local_data(:, idim), &
1693 -op_fm_set(2, ispin)%local_data(:, idim), dp)
1694 END DO
1695 IF (.NOT. all_occupied(ispin)) THEN
1696 CALL berry_phase_remove_unocc(eigrmat(ispin), mo_occupied(:, ispin))
1697 END IF
1698 ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
1699 CALL cp_cfm_det(eigrmat(ispin), zdeta)
1700 zdet = zdet*zdeta
1701 IF (dft_control%nspins == 1) THEN
1702 zdet = zdet*zdeta
1703 END IF
1704 END DO
1705 zij(i, j) = zdet*xphase(i)*xphase(j)
1706 zij(j, i) = zdet*xphase(i)*xphase(j)
1707 END DO
1708 END DO
1709 CASE (3)
1710 ! Octapole
1711 cpabort("Berry phase moments bigger than 2 not implemented")
1712 CASE (4)
1713 ! Hexadecapole
1714 cpabort("Berry phase moments bigger than 3 not implemented")
1715 CASE DEFAULT
1716 cpabort("Berry phase moments bigger than 4 not implemented")
1717 END SELECT
1718 END DO
1719 DO l = 1, nmom
1720 SELECT CASE (l)
1721 CASE (1)
1722 ! Dipole (apply periodic (2 Pi) boundary conditions)
1723 ci = aimag(log(zi))
1724 DO i = 1, 3
1725 IF (qq(i) + ci(i) > pi) ci(i) = ci(i) - twopi
1726 IF (qq(i) + ci(i) < -pi) ci(i) = ci(i) + twopi
1727 END DO
1728 rmom(2:4, 1) = matmul(cell%hmat, ci)/twopi
1729 CASE (2)
1730 ! Quadrupole
1731 cpabort("Berry phase moments bigger than 1 not implemented")
1732 DO i = 1, 3
1733 DO j = 1, 3
1734 zz = zij(i, j)/zi(i)/zi(j)
1735 cij(i, j) = aimag(log(zz))/twopi
1736 END DO
1737 END DO
1738 cij = 0.5_dp*cij/twopi/twopi
1739 cij = matmul(matmul(cell%hmat, cij), transpose(cell%hmat))
1740 DO k = 4, 9
1741 ix = indco(1, k + 1)
1742 iy = indco(2, k + 1)
1743 iz = indco(3, k + 1)
1744 IF (ix == 0) THEN
1745 rmom(k + 1, 1) = cij(iy, iz)
1746 ELSE IF (iy == 0) THEN
1747 rmom(k + 1, 1) = cij(ix, iz)
1748 ELSE IF (iz == 0) THEN
1749 rmom(k + 1, 1) = cij(ix, iy)
1750 END IF
1751 END DO
1752 CASE (3)
1753 ! Octapole
1754 cpabort("Berry phase moments bigger than 2 not implemented")
1755 CASE (4)
1756 ! Hexadecapole
1757 cpabort("Berry phase moments bigger than 3 not implemented")
1758 CASE DEFAULT
1759 cpabort("Berry phase moments bigger than 4 not implemented")
1760 END SELECT
1761 END DO
1762
1763 rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
1764 description = "[DIPOLE]"
1765 CALL cp_results_erase(results=results, description=description)
1766 CALL put_results(results=results, description=description, &
1767 values=rmom(2:4, 3))
1768 IF (magnetic) THEN
1769 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.true., mmom=mmom, &
1770 unit_str=unit_str)
1771 ELSE
1772 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.true., unit_str=unit_str)
1773 END IF
1774
1775 DEALLOCATE (rmom)
1776 DEALLOCATE (rlab)
1777 IF (magnetic) THEN
1778 DEALLOCATE (mmom)
1779 END IF
1780
1781 CALL dbcsr_deallocate_matrix(cosmat)
1782 CALL dbcsr_deallocate_matrix(sinmat)
1783
1784 CALL cp_fm_release(opvec)
1785 CALL cp_fm_release(op_fm_set)
1786 DO ispin = 1, dft_control%nspins
1787 CALL cp_cfm_release(eigrmat(ispin))
1788 END DO
1789 DEALLOCATE (eigrmat)
1790 DEALLOCATE (mo_occupied, all_occupied)
1791
1792 CALL timestop(handle)
1793
1794 END SUBROUTINE qs_moment_berry_phase
1795
1796! **************************************************************************************************
1797!> \brief Removes unoccupied MOs from the Berry phase matrix Z_ij = <psi_i|exp(i k.r)|psi_j> by
1798!> replacing their rows and columns with those of the unit matrix. The determinant of the
1799!> result equals the determinant of Z.
1800!> \param zmat nmo x nmo Berry phase matrix of one spin channel, modified on output
1801!> \param mo_occupied .TRUE. for the MOs that are part of the Slater determinant
1802! **************************************************************************************************
1803 SUBROUTINE berry_phase_remove_unocc(zmat, mo_occupied)
1804
1805 TYPE(cp_cfm_type), INTENT(IN) :: zmat
1806 LOGICAL, DIMENSION(:), INTENT(IN) :: mo_occupied
1807
1808 INTEGER :: icol, irow, ncol_local, nrow_local
1809 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1810
1811 NULLIFY (col_indices, row_indices)
1812 CALL cp_cfm_get_info(zmat, nrow_local=nrow_local, ncol_local=ncol_local, &
1813 row_indices=row_indices, col_indices=col_indices)
1814 DO icol = 1, ncol_local
1815 DO irow = 1, nrow_local
1816 IF (mo_occupied(row_indices(irow)) .AND. mo_occupied(col_indices(icol))) cycle
1817 IF (row_indices(irow) == col_indices(icol)) THEN
1818 zmat%local_data(irow, icol) = z_one
1819 ELSE
1820 zmat%local_data(irow, icol) = z_zero
1821 END IF
1822 END DO
1823 END DO
1824
1825 END SUBROUTINE berry_phase_remove_unocc
1826
1827! **************************************************************************************************
1828!> \brief ...
1829!> \param cosmat ...
1830!> \param sinmat ...
1831!> \param mos ...
1832!> \param op_fm_set ...
1833!> \param opvec ...
1834! **************************************************************************************************
1835 SUBROUTINE op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1836
1837 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1838 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1839 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
1840 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: opvec
1841
1842 INTEGER :: i, nao, nmo
1843 TYPE(cp_fm_type), POINTER :: mo_coeff
1844
1845 DO i = 1, SIZE(op_fm_set, 2) ! spin
1846 CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
1847 CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff, opvec(i), ncol=nmo)
1848 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
1849 op_fm_set(1, i))
1850 CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff, opvec(i), ncol=nmo)
1851 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
1852 op_fm_set(2, i))
1853 END DO
1854
1855 END SUBROUTINE op_orbbas
1856
1857! **************************************************************************************************
1858!> \brief ...
1859!> \param cosmat ...
1860!> \param sinmat ...
1861!> \param mos ...
1862!> \param op_fm_set ...
1863!> \param mos_new ...
1864! **************************************************************************************************
1865 SUBROUTINE op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1866
1867 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1868 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1869 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
1870 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
1871
1872 INTEGER :: i, icol, lcol, nao, newdim, nmo
1873 LOGICAL :: double_col, double_row
1874 TYPE(cp_fm_struct_type), POINTER :: newstruct, newstruct1
1875 TYPE(cp_fm_type) :: work, work1, work2
1876 TYPE(cp_fm_type), POINTER :: mo_coeff
1877
1878 DO i = 1, SIZE(op_fm_set, 2) ! spin
1879 CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
1880 CALL cp_fm_get_info(mos_new(2*i), ncol_local=lcol, ncol_global=nmo)
1881 double_col = .true.
1882 double_row = .false.
1883 CALL cp_fm_struct_double(newstruct, &
1884 mos_new(2*i)%matrix_struct, &
1885 mos_new(2*i)%matrix_struct%context, &
1886 double_col, &
1887 double_row)
1888
1889 CALL cp_fm_create(work, matrix_struct=newstruct)
1890 CALL cp_fm_create(work1, matrix_struct=newstruct)
1891 CALL cp_fm_create(work2, matrix_struct=newstruct)
1892 CALL cp_fm_get_info(work, ncol_global=newdim)
1893
1894 CALL cp_fm_set_all(work, 0.0_dp, 0.0_dp)
1895 DO icol = 1, lcol
1896 work%local_data(:, icol) = mos_new(2*i - 1)%local_data(:, icol)
1897 work%local_data(:, icol + lcol) = mos_new(2*i)%local_data(:, icol)
1898 END DO
1899
1900 CALL cp_dbcsr_sm_fm_multiply(cosmat, work, work1, ncol=newdim)
1901 CALL cp_dbcsr_sm_fm_multiply(sinmat, work, work2, ncol=newdim)
1902
1903 DO icol = 1, lcol
1904 work%local_data(:, icol) = work1%local_data(:, icol) - work2%local_data(:, icol + lcol)
1905 work%local_data(:, icol + lcol) = work1%local_data(:, icol + lcol) + work2%local_data(:, icol)
1906 END DO
1907
1908 CALL cp_fm_release(work1)
1909 CALL cp_fm_release(work2)
1910
1911 CALL cp_fm_struct_double(newstruct1, &
1912 op_fm_set(1, i)%matrix_struct, &
1913 op_fm_set(1, i)%matrix_struct%context, &
1914 double_col, &
1915 double_row)
1916
1917 CALL cp_fm_create(work1, matrix_struct=newstruct1)
1918
1919 CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i - 1), &
1920 work, 0.0_dp, work1)
1921
1922 DO icol = 1, lcol
1923 op_fm_set(1, i)%local_data(:, icol) = work1%local_data(:, icol)
1924 op_fm_set(2, i)%local_data(:, icol) = work1%local_data(:, icol + lcol)
1925 END DO
1926
1927 CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i), &
1928 work, 0.0_dp, work1)
1929
1930 DO icol = 1, lcol
1931 op_fm_set(1, i)%local_data(:, icol) = &
1932 op_fm_set(1, i)%local_data(:, icol) + work1%local_data(:, icol + lcol)
1933 op_fm_set(2, i)%local_data(:, icol) = &
1934 op_fm_set(2, i)%local_data(:, icol) - work1%local_data(:, icol)
1935 END DO
1936
1937 CALL cp_fm_release(work)
1938 CALL cp_fm_release(work1)
1939 CALL cp_fm_struct_release(newstruct)
1940 CALL cp_fm_struct_release(newstruct1)
1941
1942 END DO
1943
1944 END SUBROUTINE op_orbbas_rtp
1945
1946! **************************************************************************************************
1947!> \brief ...
1948!> \param qs_env ...
1949!> \param magnetic ...
1950!> \param nmoments ...
1951!> \param reference ...
1952!> \param ref_point ...
1953!> \param unit_number ...
1954!> \param vel_reprs ...
1955!> \param com_nl ...
1956!> \param unit_str unit of the printed moments (dipole unit of cp_units), default debye
1957! **************************************************************************************************
1958 SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl, &
1959 unit_str)
1960
1961 TYPE(qs_environment_type), POINTER :: qs_env
1962 LOGICAL, INTENT(IN) :: magnetic
1963 INTEGER, INTENT(IN) :: nmoments, reference
1964 REAL(dp), DIMENSION(:), INTENT(IN), POINTER :: ref_point
1965 INTEGER, INTENT(IN) :: unit_number
1966 LOGICAL, INTENT(IN), OPTIONAL :: vel_reprs, com_nl
1967 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: unit_str
1968
1969 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_locop'
1970
1971 CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
1972 CHARACTER(LEN=default_string_length) :: description
1973 INTEGER :: akind, handle, i, ia, iatom, idir, &
1974 ikind, ispin, ix, iy, iz, l, nm, nmom, &
1975 order
1976 LOGICAL :: my_com_nl, my_velreprs, use_periodic_images
1977 REAL(dp) :: charge, dd, strace, trace
1978 REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rrv_vrr, &
1979 nlcom_rv, nlcom_rvr, nlcom_rxrv, &
1980 qupole_der, rmom_vel
1981 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
1982 REAL(dp), DIMENSION(3) :: rcc, ria
1983 TYPE(atomic_kind_type), POINTER :: atomic_kind
1984 TYPE(cell_type), POINTER :: cell
1985 TYPE(cp_result_type), POINTER :: results
1986 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom, matrix_s, moments, momentum, &
1987 rho_ao
1988 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_der
1989 TYPE(dbcsr_type), POINTER :: tmp_ao
1990 TYPE(dft_control_type), POINTER :: dft_control
1991 TYPE(distribution_1d_type), POINTER :: local_particles
1992 TYPE(mp_para_env_type), POINTER :: para_env
1993 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1994 POINTER :: sab_all, sab_orb
1995 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1996 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1997 TYPE(qs_rho_type), POINTER :: rho
1998
1999 cpassert(ASSOCIATED(qs_env))
2000
2001 CALL timeset(routinen, handle)
2002
2003 my_velreprs = .false.
2004 IF (PRESENT(vel_reprs)) my_velreprs = vel_reprs
2005 IF (PRESENT(com_nl)) my_com_nl = com_nl
2006 IF (my_velreprs) CALL cite_reference(mattiat2019)
2007
2008 ! reference point
2009 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
2010
2011 ! only allow for moments up to maxl set by basis
2012 nmom = min(nmoments, current_maxl)
2013 ! electronic contribution
2014 NULLIFY (dft_control, rho, cell, particle_set, qs_kind_set, results, para_env, matrix_s, rho_ao, sab_all, sab_orb)
2015 CALL get_qs_env(qs_env, &
2016 dft_control=dft_control, &
2017 rho=rho, &
2018 cell=cell, &
2019 results=results, &
2020 particle_set=particle_set, &
2021 qs_kind_set=qs_kind_set, &
2022 para_env=para_env, &
2023 matrix_s=matrix_s, &
2024 sab_all=sab_all, &
2025 sab_orb=sab_orb)
2026
2027 IF (my_com_nl) THEN
2028 IF ((nmom >= 1) .AND. my_velreprs) THEN
2029 ALLOCATE (nlcom_rv(3))
2030 nlcom_rv(:) = 0._dp
2031 END IF
2032 IF ((nmom >= 2) .AND. my_velreprs) THEN
2033 ALLOCATE (nlcom_rrv(6))
2034 nlcom_rrv(:) = 0._dp
2035 ALLOCATE (nlcom_rvr(6))
2036 nlcom_rvr(:) = 0._dp
2037 ALLOCATE (nlcom_rrv_vrr(6))
2038 nlcom_rrv_vrr(:) = 0._dp
2039 END IF
2040 IF (magnetic) THEN
2041 ALLOCATE (nlcom_rxrv(3))
2042 nlcom_rxrv = 0._dp
2043 END IF
2044 ! Calculate non local correction terms
2045 CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr, rcc)
2046 END IF
2047
2048 NULLIFY (moments)
2049 nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
2050 CALL dbcsr_allocate_matrix_set(moments, nm)
2051 DO i = 1, nm
2052 ALLOCATE (moments(i)%matrix)
2053 IF (my_velreprs .AND. (nmom >= 2)) THEN
2054 CALL dbcsr_create(moments(i)%matrix, template=matrix_s(1)%matrix, &
2055 matrix_type=dbcsr_type_symmetric)
2056 CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_orb)
2057 ELSE
2058 CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
2059 END IF
2060 CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
2061 END DO
2062
2063 ! calculate derivatives if quadrupole in vel. reprs. is requested
2064 IF (my_velreprs .AND. (nmom >= 2)) THEN
2065 NULLIFY (moments_der)
2066 CALL dbcsr_allocate_matrix_set(moments_der, 3, 3)
2067 DO i = 1, 3 ! x, y, z
2068 DO idir = 1, 3 ! d/dx, d/dy, d/dz
2069 CALL dbcsr_init_p(moments_der(i, idir)%matrix)
2070 CALL dbcsr_create(moments_der(i, idir)%matrix, template=matrix_s(1)%matrix, &
2071 matrix_type=dbcsr_type_antisymmetric)
2072 CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_orb)
2073 CALL dbcsr_set(moments_der(i, idir)%matrix, 0.0_dp)
2074 END DO
2075 END DO
2076 CALL build_local_moments_der_matrix(qs_env, moments_der, 1, 2, ref_point=rcc, moments=moments)
2077 ELSE
2078 ! In non-periodic directions the moment must include the same images as the overlap.
2079 use_periodic_images = any(cell%perd /= 0) .AND. any(cell%perd == 0)
2080 CALL build_local_moment_matrix(qs_env, moments, nmom, ref_point=rcc, &
2081 all_images=use_periodic_images, neighbor_image=use_periodic_images)
2082 END IF
2083
2084 CALL qs_rho_get(rho, rho_ao=rho_ao)
2085
2086 ALLOCATE (rmom(nm + 1, 3))
2087 ALLOCATE (rlab(nm + 1))
2088 rmom = 0.0_dp
2089 rlab = ""
2090
2091 IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2092 ! Allocate matrix to store the matrix product to be traced (dbcsr_dot only works for products of
2093 ! symmetric matrices)
2094 NULLIFY (tmp_ao)
2095 CALL dbcsr_init_p(tmp_ao)
2096 CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
2097 CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
2098 CALL dbcsr_set(tmp_ao, 0.0_dp)
2099 END IF
2100
2101 trace = 0.0_dp
2102 DO ispin = 1, dft_control%nspins
2103 CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
2104 rmom(1, 1) = rmom(1, 1) + trace
2105 END DO
2106
2107 DO i = 1, SIZE(moments)
2108 strace = 0._dp
2109 DO ispin = 1, dft_control%nspins
2110 IF (my_velreprs .AND. nmoments >= 2) THEN
2111 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, moments(i)%matrix, &
2112 0.0_dp, tmp_ao)
2113 CALL dbcsr_trace(tmp_ao, trace)
2114 ELSE
2115 CALL dbcsr_dot(rho_ao(ispin)%matrix, moments(i)%matrix, trace)
2116 END IF
2117 strace = strace + trace
2118 END DO
2119 rmom(i + 1, 1) = strace
2120 END DO
2121
2122 CALL dbcsr_deallocate_matrix_set(moments)
2123
2124 ! nuclear contribution
2125 CALL get_qs_env(qs_env=qs_env, &
2126 local_particles=local_particles)
2127 DO ikind = 1, SIZE(local_particles%n_el)
2128 DO ia = 1, local_particles%n_el(ikind)
2129 iatom = local_particles%list(ikind)%array(ia)
2130 ! fold atomic positions back into unit cell
2131 ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
2132 ria = ria - rcc
2133 atomic_kind => particle_set(iatom)%atomic_kind
2134 CALL get_atomic_kind(atomic_kind, kind_number=akind)
2135 CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
2136 rmom(1, 2) = rmom(1, 2) - charge
2137 DO l = 1, nm
2138 ix = indco(1, l + 1)
2139 iy = indco(2, l + 1)
2140 iz = indco(3, l + 1)
2141 dd = 1._dp
2142 IF (ix > 0) dd = dd*ria(1)**ix
2143 IF (iy > 0) dd = dd*ria(2)**iy
2144 IF (iz > 0) dd = dd*ria(3)**iz
2145 rmom(l + 1, 2) = rmom(l + 1, 2) - charge*dd
2146 CALL set_label(rlab(l + 1), ix, iy, iz)
2147 END DO
2148 END DO
2149 END DO
2150 CALL para_env%sum(rmom(:, 2))
2151 rmom(:, :) = -rmom(:, :)
2152 rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
2153
2154 ! magnetic moments
2155 IF (magnetic) THEN
2156 NULLIFY (magmom)
2157 CALL dbcsr_allocate_matrix_set(magmom, 3)
2158 DO i = 1, SIZE(magmom)
2159 CALL dbcsr_init_p(magmom(i)%matrix)
2160 CALL dbcsr_create(magmom(i)%matrix, template=matrix_s(1)%matrix, &
2161 matrix_type=dbcsr_type_antisymmetric)
2162 CALL cp_dbcsr_alloc_block_from_nbl(magmom(i)%matrix, sab_orb)
2163 CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
2164 END DO
2165
2166 CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
2167
2168 ALLOCATE (mmom(SIZE(magmom)))
2169 mmom(:) = 0.0_dp
2170 IF (qs_env%run_rtp) THEN
2171 ! get imaginary part of the density in rho_ao (the real part is not needed since the trace of the product
2172 ! of a symmetric (REAL(rho_ao)) and an anti-symmetric (L_AO) matrix is zero)
2173 ! There may be other cases, where the imaginary part of the density is relevant
2174 NULLIFY (rho_ao)
2175 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2176 END IF
2177 ! if the density is purely real this is an expensive way to calculate zero
2178 DO i = 1, SIZE(magmom)
2179 strace = 0._dp
2180 DO ispin = 1, dft_control%nspins
2181 CALL dbcsr_set(tmp_ao, 0.0_dp)
2182 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, magmom(i)%matrix, &
2183 0.0_dp, tmp_ao)
2184 CALL dbcsr_trace(tmp_ao, trace)
2185 strace = strace + trace
2186 END DO
2187 mmom(i) = strace
2188 END DO
2189
2190 CALL dbcsr_deallocate_matrix_set(magmom)
2191 END IF
2192
2193 ! velocity representations
2194 IF (my_velreprs) THEN
2195 ALLOCATE (rmom_vel(nm))
2196 rmom_vel = 0.0_dp
2197
2198 DO order = 1, nmom
2199 SELECT CASE (order)
2200
2201 CASE (1) ! expectation value of momentum
2202 NULLIFY (momentum)
2203 CALL dbcsr_allocate_matrix_set(momentum, 3)
2204 DO i = 1, 3
2205 CALL dbcsr_init_p(momentum(i)%matrix)
2206 CALL dbcsr_create(momentum(i)%matrix, template=matrix_s(1)%matrix, &
2207 matrix_type=dbcsr_type_antisymmetric)
2208 CALL cp_dbcsr_alloc_block_from_nbl(momentum(i)%matrix, sab_orb)
2209 CALL dbcsr_set(momentum(i)%matrix, 0.0_dp)
2210 END DO
2211 CALL build_lin_mom_matrix(qs_env, momentum)
2212
2213 ! imaginary part of the density for RTP, real part gives 0 since momentum is antisymmetric
2214 IF (qs_env%run_rtp) THEN
2215 NULLIFY (rho_ao)
2216 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2217 DO idir = 1, SIZE(momentum)
2218 strace = 0._dp
2219 DO ispin = 1, dft_control%nspins
2220 CALL dbcsr_set(tmp_ao, 0.0_dp)
2221 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, momentum(idir)%matrix, &
2222 0.0_dp, tmp_ao)
2223 CALL dbcsr_trace(tmp_ao, trace)
2224 strace = strace + trace
2225 END DO
2226 rmom_vel(idir) = rmom_vel(idir) + strace
2227 END DO
2228 END IF
2229
2230 CALL dbcsr_deallocate_matrix_set(momentum)
2231
2232 CASE (2) ! expectation value of quadrupole moment in vel. reprs.
2233 ALLOCATE (qupole_der(9)) ! will contain the expectation values of r_\alpha * d/d r_\beta
2234 qupole_der = 0._dp
2235
2236 NULLIFY (rho_ao)
2237 CALL qs_rho_get(rho, rho_ao=rho_ao)
2238
2239 ! Calculate expectation value over real part
2240 trace = 0._dp
2241 DO i = 1, 3
2242 DO idir = 1, 3
2243 strace = 0._dp
2244 DO ispin = 1, dft_control%nspins
2245 CALL dbcsr_set(tmp_ao, 0._dp)
2246 CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2247 CALL dbcsr_trace(tmp_ao, trace)
2248 strace = strace + trace
2249 END DO
2250 qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2251 END DO
2252 END DO
2253
2254 IF (qs_env%run_rtp) THEN
2255 NULLIFY (rho_ao)
2256 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2257
2258 ! Calculate expectation value over imaginary part
2259 trace = 0._dp
2260 DO i = 1, 3
2261 DO idir = 1, 3
2262 strace = 0._dp
2263 DO ispin = 1, dft_control%nspins
2264 CALL dbcsr_set(tmp_ao, 0._dp)
2265 CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2266 CALL dbcsr_trace(tmp_ao, trace)
2267 strace = strace + trace
2268 END DO
2269 qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2270 END DO
2271 END DO
2272 END IF
2273
2274 ! calculate vel. reprs. of quadrupole moment from derivatives
2275 rmom_vel(4) = -2*qupole_der(1) - rmom(1, 1)
2276 rmom_vel(5) = -qupole_der(2) - qupole_der(4)
2277 rmom_vel(6) = -qupole_der(3) - qupole_der(7)
2278 rmom_vel(7) = -2*qupole_der(5) - rmom(1, 1)
2279 rmom_vel(8) = -qupole_der(6) - qupole_der(8)
2280 rmom_vel(9) = -2*qupole_der(9) - rmom(1, 1)
2281
2282 DEALLOCATE (qupole_der)
2283 CASE DEFAULT
2284 END SELECT
2285 END DO
2286 END IF
2287
2288 IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2289 CALL dbcsr_deallocate_matrix(tmp_ao)
2290 END IF
2291 IF (my_velreprs .AND. (nmoments >= 2)) THEN
2292 CALL dbcsr_deallocate_matrix_set(moments_der)
2293 END IF
2294
2295 description = "[DIPOLE]"
2296 CALL cp_results_erase(results=results, description=description)
2297 CALL put_results(results=results, description=description, &
2298 values=rmom(2:4, 3))
2299
2300 IF (unit_number > 0 .AND. any(cell%perd /= 0) .AND. any(cell%perd == 0)) THEN
2301 WRITE (unit_number, "(T3,A)") &
2302 "Only dipole components perpendicular to the periodic lattice vectors are well-defined."
2303 END IF
2304
2305 IF (magnetic .AND. my_velreprs) THEN
2306 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., mmom=mmom, rmom_vel=rmom_vel, &
2307 unit_str=unit_str)
2308 ELSE IF (magnetic .AND. .NOT. my_velreprs) THEN
2309 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., mmom=mmom, &
2310 unit_str=unit_str)
2311 ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2312 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., rmom_vel=rmom_vel, &
2313 unit_str=unit_str)
2314 ELSE
2315 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., unit_str=unit_str)
2316 END IF
2317
2318 IF (my_com_nl) THEN
2319 IF (magnetic) THEN
2320 mmom(:) = nlcom_rxrv(:)
2321 END IF
2322 IF (my_velreprs .AND. (nmom >= 1)) THEN
2323 DEALLOCATE (rmom_vel)
2324 ALLOCATE (rmom_vel(21))
2325 rmom_vel(1:3) = nlcom_rv
2326 END IF
2327 IF (my_velreprs .AND. (nmom >= 2)) THEN
2328 rmom_vel(4:9) = nlcom_rrv
2329 rmom_vel(10:15) = nlcom_rvr
2330 rmom_vel(16:21) = nlcom_rrv_vrr
2331 END IF
2332 IF (magnetic .AND. .NOT. my_velreprs) THEN
2333 CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom)
2334 ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2335 CALL print_moments_nl(unit_number, nmom, rlab, rmom_vel=rmom_vel)
2336 ELSE IF (my_velreprs .AND. magnetic) THEN
2337 CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom, rmom_vel=rmom_vel)
2338 END IF
2339
2340 END IF
2341
2342 IF (my_com_nl) THEN
2343 IF (nmom >= 1 .AND. my_velreprs) DEALLOCATE (nlcom_rv)
2344 IF (nmom >= 2 .AND. my_velreprs) THEN
2345 DEALLOCATE (nlcom_rrv)
2346 DEALLOCATE (nlcom_rvr)
2347 DEALLOCATE (nlcom_rrv_vrr)
2348 END IF
2349 IF (magnetic) DEALLOCATE (nlcom_rxrv)
2350 END IF
2351
2352 DEALLOCATE (rmom)
2353 DEALLOCATE (rlab)
2354 IF (magnetic) THEN
2355 DEALLOCATE (mmom)
2356 END IF
2357 IF (my_velreprs) THEN
2358 DEALLOCATE (rmom_vel)
2359 END IF
2360
2361 CALL timestop(handle)
2362
2363 END SUBROUTINE qs_moment_locop
2364
2365! **************************************************************************************************
2366!> \brief ...
2367!> \param label ...
2368!> \param ix ...
2369!> \param iy ...
2370!> \param iz ...
2371! **************************************************************************************************
2372 SUBROUTINE set_label(label, ix, iy, iz)
2373 CHARACTER(LEN=*), INTENT(OUT) :: label
2374 INTEGER, INTENT(IN) :: ix, iy, iz
2375
2376 INTEGER :: i
2377
2378 label = ""
2379 DO i = 1, ix
2380 WRITE (label(i:), "(A1)") "X"
2381 END DO
2382 DO i = ix + 1, ix + iy
2383 WRITE (label(i:), "(A1)") "Y"
2384 END DO
2385 DO i = ix + iy + 1, ix + iy + iz
2386 WRITE (label(i:), "(A1)") "Z"
2387 END DO
2388
2389 END SUBROUTINE set_label
2390
2391! **************************************************************************************************
2392!> \brief ...
2393!> \param unit_number ...
2394!> \param nmom ...
2395!> \param rmom ...
2396!> \param rlab ...
2397!> \param rcc ...
2398!> \param cell ...
2399!> \param periodic ...
2400!> \param mmom ...
2401!> \param rmom_vel ...
2402!> \param unit_str unit of the printed electric moments (a dipole unit of cp_units), default debye
2403! **************************************************************************************************
2404 SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmom, rmom_vel, unit_str)
2405 INTEGER, INTENT(IN) :: unit_number, nmom
2406 REAL(dp), DIMENSION(:, :), INTENT(IN) :: rmom
2407 CHARACTER(LEN=8), DIMENSION(:) :: rlab
2408 REAL(dp), DIMENSION(3), INTENT(IN) :: rcc
2409 TYPE(cell_type), POINTER :: cell
2410 LOGICAL :: periodic
2411 REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2412 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: unit_str
2413
2414 CHARACTER(LEN=16) :: dip_label, fnum, ftot, len_label
2415 INTEGER :: i, i0, i1, j, l
2416 LOGICAL :: use_sci
2417 REAL(dp) :: dd, dip_fac, len_fac
2418
2419 IF (unit_number > 0) THEN
2420 IF (PRESENT(unit_str)) THEN
2421 CALL get_moment_print_units(unit_str, dip_fac, len_fac, dip_label, len_label, use_sci)
2422 ELSE
2423 CALL get_moment_print_units("debye", dip_fac, len_fac, dip_label, len_label, use_sci)
2424 END IF
2425 ! fixed-point fields would print SI values (~1e-30) as zero
2426 fnum = "F14.8"
2427 ftot = "F13.7"
2428 IF (use_sci) THEN
2429 fnum = "ES14.6"
2430 ftot = "ES13.5"
2431 END IF
2432 DO l = 0, nmom
2433 SELECT CASE (l)
2434 CASE (0)
2435 WRITE (unit_number, "(T3,A,T33,3F16.8)") "Reference Point [Bohr]", rcc
2436 WRITE (unit_number, "(T3,A)") "Charges"
2437 WRITE (unit_number, "(T5,A,T18,F14.8,T36,A,T42,F14.8,T60,A,T67,F14.8)") &
2438 "Electronic=", rmom(1, 1), "Core=", rmom(1, 2), "Total=", rmom(1, 3)
2439 CASE (1)
2440 IF (periodic) THEN
2441 WRITE (unit_number, "(T3,A)") "Dipole vectors are based on the periodic (Berry phase) operator."
2442 WRITE (unit_number, "(T3,A)") "They are defined modulo integer multiples of the cell matrix ["// &
2443 trim(dip_label)//"]."
2444 WRITE (unit_number, "(T3,A,3("//trim(fnum)//",1X),A)") "[X] [", cell%hmat(1, :)*dip_fac, "] [i]"
2445 WRITE (unit_number, "(T3,A,3("//trim(fnum)//",1X),A)") "[Y]=[", cell%hmat(2, :)*dip_fac, "]*[j]"
2446 WRITE (unit_number, "(T3,A,3("//trim(fnum)//",1X),A)") "[Z] [", cell%hmat(3, :)*dip_fac, "] [k]"
2447 ELSE
2448 WRITE (unit_number, "(T3,A)") "Dipoles are based on the traditional operator."
2449 END IF
2450 dd = sqrt(sum(rmom(2:4, 3)**2))*dip_fac
2451 WRITE (unit_number, "(T3,A)") "Dipole moment ["//trim(dip_label)//"]"
2452 WRITE (unit_number, "(T5,3(A,A,E15.7,1X),T60,A,T68,"//trim(ftot)//")") &
2453 (trim(rlab(i)), "=", rmom(i, 3)*dip_fac, i=2, 4), "Total=", dd
2454 CASE (2)
2455 dd = dip_fac*len_fac
2456 WRITE (unit_number, "(T3,A)") "Quadrupole moment ["//trim(dip_label)//"*"//trim(len_label)//"]"
2457 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2458 (trim(rlab(i)), "=", rmom(i, 3)*dd, i=5, 7)
2459 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2460 (trim(rlab(i)), "=", rmom(i, 3)*dd, i=8, 10)
2461 CASE (3)
2462 dd = dip_fac*len_fac**2
2463 WRITE (unit_number, "(T3,A)") "Octapole moment ["//trim(dip_label)//"*"//trim(len_label)//"**2]"
2464 WRITE (unit_number, "(T7,4(A,A,"//trim(fnum)//",3X))") &
2465 (trim(rlab(i)), "=", rmom(i, 3)*dd, i=11, 14)
2466 WRITE (unit_number, "(T7,4(A,A,"//trim(fnum)//",3X))") &
2467 (trim(rlab(i)), "=", rmom(i, 3)*dd, i=15, 18)
2468 WRITE (unit_number, "(T7,4(A,A,"//trim(fnum)//",3X))") &
2469 (trim(rlab(i)), "=", rmom(i, 3)*dd, i=19, 20)
2470 CASE (4)
2471 dd = dip_fac*len_fac**3
2472 WRITE (unit_number, "(T3,A)") "Hexadecapole moment ["//trim(dip_label)//"*"//trim(len_label)//"**3]"
2473 WRITE (unit_number, "(T6,4(A,A,"//trim(fnum)//",2X))") &
2474 (trim(rlab(i)), "=", rmom(i, 3)*dd, i=21, 24)
2475 WRITE (unit_number, "(T6,4(A,A,"//trim(fnum)//",2X))") &
2476 (trim(rlab(i)), "=", rmom(i, 3)*dd, i=25, 28)
2477 WRITE (unit_number, "(T6,4(A,A,"//trim(fnum)//",2X))") &
2478 (trim(rlab(i)), "=", rmom(i, 3)*dd, i=29, 32)
2479 WRITE (unit_number, "(T6,4(A,A,"//trim(fnum)//",2X))") &
2480 (trim(rlab(i)), "=", rmom(i, 3)*dd, i=32, 35)
2481 CASE DEFAULT
2482 WRITE (unit_number, "(T3,A,A,I2)") "Higher moment ["//trim(dip_label)//"*"//trim(len_label)// &
2483 "**(L-1)]", " L=", l
2484 i0 = (6 + 11*(l - 1) + 6*(l - 1)**2 + (l - 1)**3)/6
2485 i1 = (6 + 11*l + 6*l**2 + l**3)/6 - 1
2486 dd = dip_fac*len_fac**(l - 1)
2487 DO i = i0, i1, 3
2488 WRITE (unit_number, "(T18,3(A,A,"//trim(fnum)//",4X))") &
2489 (trim(rlab(j + 1)), "=", rmom(j + 1, 3)*dd, j=i, min(i1, i + 2))
2490 END DO
2491 END SELECT
2492 END DO
2493 IF (PRESENT(mmom)) THEN
2494 IF (nmom >= 1) THEN
2495 dd = sqrt(sum(mmom(1:3)**2))
2496 WRITE (unit_number, "(T3,A)") "Orbital angular momentum [a. u.]"
2497 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2498 (trim(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2499 END IF
2500 END IF
2501 IF (PRESENT(rmom_vel)) THEN
2502 DO l = 1, nmom
2503 SELECT CASE (l)
2504 CASE (1)
2505 dd = sqrt(sum(rmom_vel(1:3)**2))
2506 WRITE (unit_number, "(T3,A)") "Expectation value of momentum operator [a. u.]"
2507 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2508 (trim(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2509 CASE (2)
2510 WRITE (unit_number, "(T3,A)") "Expectation value of quadrupole operator in vel. repr. [a. u.]"
2511 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2512 (trim(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2513 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2514 (trim(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2515 CASE DEFAULT
2516 END SELECT
2517 END DO
2518 END IF
2519 END IF
2520
2521 END SUBROUTINE print_moments
2522
2523! **************************************************************************************************
2524!> \brief Conversion factors and labels for printing electric moments in a chosen dipole unit.
2525!> Multipoles of order L are printed in (dipole unit)*(length unit)^(L-1), where the length
2526!> unit is paired with the dipole unit (debye, e_angstrom -> angstrom; au_dip -> bohr;
2527!> C_m -> m).
2528!> \param unit_str dipole unit understood by cp_units (debye, e_angstrom, au_dip, C_m)
2529!> \param dip_fac factor converting a dipole from a.u. (e*bohr) to unit_str
2530!> \param len_fac factor converting a length from bohr to the paired length unit
2531!> \param dip_label label of the dipole unit for the output headers
2532!> \param len_label label of the paired length unit for the output headers
2533!> \param use_sci .TRUE. if the values are too small for fixed-point output (SI units)
2534! **************************************************************************************************
2535 SUBROUTINE get_moment_print_units(unit_str, dip_fac, len_fac, dip_label, len_label, use_sci)
2536 CHARACTER(LEN=*), INTENT(IN) :: unit_str
2537 REAL(dp), INTENT(OUT) :: dip_fac, len_fac
2538 CHARACTER(LEN=*), INTENT(OUT) :: dip_label, len_label
2539 LOGICAL, INTENT(OUT) :: use_sci
2540
2541 CHARACTER(LEN=16) :: len_unit
2542 TYPE(cp_unit_type) :: unit
2543
2544 CALL cp_unit_create(unit, unit_str)
2545 IF (unit%n_kinds /= 1 .OR. unit%kind_id(1) /= cp_ukind_dipole .OR. unit%power(1) /= 1) THEN
2546 CALL cp_abort(__location__, "MOMENTS%UNIT must be a dipole moment unit "// &
2547 "(debye, e_angstrom, au_dip or C_m), got: "//trim(unit_str))
2548 END IF
2549 use_sci = .false.
2550 SELECT CASE (unit%unit_id(1))
2551 CASE (cp_units_debye)
2552 dip_label = "Debye"
2553 len_label = "Angstrom"
2554 len_unit = "angstrom"
2555 CASE (cp_units_e_angstrom)
2556 dip_label = "e*Angstrom"
2557 len_label = "Angstrom"
2558 len_unit = "angstrom"
2559 CASE (cp_units_au)
2560 dip_label = "a.u."
2561 len_label = "Bohr"
2562 len_unit = "bohr"
2563 CASE (cp_units_c_m)
2564 dip_label = "C*m"
2565 len_label = "m"
2566 len_unit = "m"
2567 use_sci = .true.
2568 CASE DEFAULT
2569 CALL cp_abort(__location__, "MOMENTS%UNIT: dipole unit not fully specified: "//trim(unit_str))
2570 END SELECT
2571 dip_fac = cp_unit_from_cp2k1(1.0_dp, unit)
2572 len_fac = cp_unit_from_cp2k(1.0_dp, trim(len_unit))
2573 CALL cp_unit_release(unit)
2574
2575 END SUBROUTINE get_moment_print_units
2576
2577! **************************************************************************************************
2578!> \brief ...
2579!> \param unit_number ...
2580!> \param nmom ...
2581!> \param rlab ...
2582!> \param mmom ...
2583!> \param rmom_vel ...
2584! **************************************************************************************************
2585 SUBROUTINE print_moments_nl(unit_number, nmom, rlab, mmom, rmom_vel)
2586 INTEGER, INTENT(IN) :: unit_number, nmom
2587 CHARACTER(LEN=8), DIMENSION(:) :: rlab
2588 REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2589
2590 INTEGER :: i, l
2591 REAL(dp) :: dd
2592
2593 IF (unit_number > 0) THEN
2594 IF (PRESENT(mmom)) THEN
2595 IF (nmom >= 1) THEN
2596 dd = sqrt(sum(mmom(1:3)**2))
2597 WRITE (unit_number, "(T3,A)") "Expectation value of rx[r,V_nl] [a. u.]"
2598 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2599 (trim(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2600 END IF
2601 END IF
2602 IF (PRESENT(rmom_vel)) THEN
2603 DO l = 1, nmom
2604 SELECT CASE (l)
2605 CASE (1)
2606 dd = sqrt(sum(rmom_vel(1:3)**2))
2607 WRITE (unit_number, "(T3,A)") "Expectation value of [r,V_nl] [a. u.]"
2608 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2609 (trim(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2610 CASE (2)
2611 WRITE (unit_number, "(T3,A)") "Expectation value of [rr,V_nl] [a. u.]"
2612 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2613 (trim(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2614 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2615 (trim(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2616 WRITE (unit_number, "(T3,A)") "Expectation value of r x V_nl x r [a. u.]"
2617 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2618 (trim(rlab(i + 1 - 6)), "=", rmom_vel(i), i=10, 12)
2619 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2620 (trim(rlab(i + 1 - 6)), "=", rmom_vel(i), i=13, 15)
2621 WRITE (unit_number, "(T3,A)") "Expectation value of r x r x V_nl + V_nl x r x r [a. u.]"
2622 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2623 (trim(rlab(i + 1 - 12)), "=", rmom_vel(i), i=16, 18)
2624 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2625 (trim(rlab(i + 1 - 12)), "=", rmom_vel(i), i=19, 21)
2626 CASE DEFAULT
2627 END SELECT
2628 END DO
2629 END IF
2630 END IF
2631
2632 END SUBROUTINE print_moments_nl
2633
2634! **************************************************************************************************
2635!> \brief Calculate the expectation value of operators related to non-local potential:
2636!> [r, Vnl], noted rv
2637!> r x [r,Vnl], noted rxrv
2638!> [rr,Vnl], noted rrv
2639!> r x Vnl x r, noted rvr
2640!> r x r x Vnl + Vnl x r x r, noted rrv_vrr
2641!> Note that the 3 first operator are commutator while the 2 last
2642!> are not. For reading clarity the same notation is used for all 5
2643!> operators.
2644!> \param qs_env ...
2645!> \param nlcom_rv ...
2646!> \param nlcom_rxrv ...
2647!> \param nlcom_rrv ...
2648!> \param nlcom_rvr ...
2649!> \param nlcom_rrv_vrr ...
2650!> \param ref_point ...
2651! **************************************************************************************************
2652 SUBROUTINE calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, &
2653 nlcom_rrv_vrr, ref_point)
2654
2655 TYPE(qs_environment_type), POINTER :: qs_env
2656 REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: nlcom_rv, nlcom_rxrv, nlcom_rrv, &
2657 nlcom_rvr, nlcom_rrv_vrr
2658 REAL(dp), DIMENSION(3) :: ref_point
2659
2660 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_commutator_nl_terms'
2661
2662 INTEGER :: handle, ind, ispin
2663 LOGICAL :: calc_rrv, calc_rrv_vrr, calc_rv, &
2664 calc_rvr, calc_rxrv
2665 REAL(dp) :: eps_ppnl, strace, trace
2666 TYPE(cell_type), POINTER :: cell
2667 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rrv, matrix_rrv_vrr, matrix_rv, &
2668 matrix_rvr, matrix_rxrv, matrix_s, &
2669 rho_ao
2670 TYPE(dbcsr_type), POINTER :: tmp_ao
2671 TYPE(dft_control_type), POINTER :: dft_control
2672 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2673 POINTER :: sab_all, sab_orb, sap_ppnl
2674 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2675 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2676 TYPE(qs_rho_type), POINTER :: rho
2677
2678 CALL timeset(routinen, handle)
2679
2680 calc_rv = .false.
2681 calc_rxrv = .false.
2682 calc_rrv = .false.
2683 calc_rvr = .false.
2684 calc_rrv_vrr = .false.
2685
2686 ! rv, rxrv and rrv are commutator matrices: anti-symmetric.
2687 ! The real part of the density matrix rho_ao is symmetric so that
2688 ! the expectation value of real density matrix is zero. Hence, if
2689 ! the density matrix is real, no need to compute these quantities.
2690 ! This is not the case for rvr and rrv_vrr which are symmetric.
2691
2692 IF (ALLOCATED(nlcom_rv)) THEN
2693 nlcom_rv(:) = 0._dp
2694 IF (qs_env%run_rtp) calc_rv = .true.
2695 END IF
2696 IF (ALLOCATED(nlcom_rxrv)) THEN
2697 nlcom_rxrv(:) = 0._dp
2698 IF (qs_env%run_rtp) calc_rxrv = .true.
2699 END IF
2700 IF (ALLOCATED(nlcom_rrv)) THEN
2701 nlcom_rrv(:) = 0._dp
2702 IF (qs_env%run_rtp) calc_rrv = .true.
2703 END IF
2704 IF (ALLOCATED(nlcom_rvr)) THEN
2705 nlcom_rvr(:) = 0._dp
2706 calc_rvr = .true.
2707 END IF
2708 IF (ALLOCATED(nlcom_rrv_vrr)) THEN
2709 nlcom_rrv_vrr(:) = 0._dp
2710 calc_rrv_vrr = .true.
2711 END IF
2712
2713 IF (.NOT. (calc_rv .OR. calc_rrv .OR. calc_rxrv .OR. calc_rvr .OR. calc_rrv_vrr)) THEN
2714 CALL timestop(handle)
2715 RETURN
2716 END IF
2717
2718 NULLIFY (cell, matrix_s, particle_set, qs_kind_set, rho, sab_all, sab_orb, sap_ppnl)
2719 CALL get_qs_env(qs_env, &
2720 cell=cell, &
2721 dft_control=dft_control, &
2722 matrix_s=matrix_s, &
2723 particle_set=particle_set, &
2724 qs_kind_set=qs_kind_set, &
2725 rho=rho, &
2726 sab_orb=sab_orb, &
2727 sab_all=sab_all, &
2728 sap_ppnl=sap_ppnl)
2729
2730 eps_ppnl = dft_control%qs_control%eps_ppnl
2731
2732 ! Allocate storage
2733 NULLIFY (matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr)
2734 IF (calc_rv) THEN
2735 CALL dbcsr_allocate_matrix_set(matrix_rv, 3)
2736 DO ind = 1, 3
2737 CALL dbcsr_init_p(matrix_rv(ind)%matrix)
2738 CALL dbcsr_create(matrix_rv(ind)%matrix, template=matrix_s(1)%matrix, &
2739 matrix_type=dbcsr_type_antisymmetric)
2740 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(ind)%matrix, sab_orb)
2741 CALL dbcsr_set(matrix_rv(ind)%matrix, 0._dp)
2742 END DO
2743 END IF
2744
2745 IF (calc_rxrv) THEN
2746 CALL dbcsr_allocate_matrix_set(matrix_rxrv, 3)
2747 DO ind = 1, 3
2748 CALL dbcsr_init_p(matrix_rxrv(ind)%matrix)
2749 CALL dbcsr_create(matrix_rxrv(ind)%matrix, template=matrix_s(1)%matrix, &
2750 matrix_type=dbcsr_type_antisymmetric)
2751 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rxrv(ind)%matrix, sab_orb)
2752 CALL dbcsr_set(matrix_rxrv(ind)%matrix, 0._dp)
2753 END DO
2754 END IF
2755
2756 IF (calc_rrv) THEN
2757 CALL dbcsr_allocate_matrix_set(matrix_rrv, 6)
2758 DO ind = 1, 6
2759 CALL dbcsr_init_p(matrix_rrv(ind)%matrix)
2760 CALL dbcsr_create(matrix_rrv(ind)%matrix, template=matrix_s(1)%matrix, &
2761 matrix_type=dbcsr_type_antisymmetric)
2762 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv(ind)%matrix, sab_orb)
2763 CALL dbcsr_set(matrix_rrv(ind)%matrix, 0._dp)
2764 END DO
2765 END IF
2766
2767 IF (calc_rvr) THEN
2768 CALL dbcsr_allocate_matrix_set(matrix_rvr, 6)
2769 DO ind = 1, 6
2770 CALL dbcsr_init_p(matrix_rvr(ind)%matrix)
2771 CALL dbcsr_create(matrix_rvr(ind)%matrix, template=matrix_s(1)%matrix, &
2772 matrix_type=dbcsr_type_symmetric)
2773 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rvr(ind)%matrix, sab_orb)
2774 CALL dbcsr_set(matrix_rvr(ind)%matrix, 0._dp)
2775 END DO
2776 END IF
2777 IF (calc_rrv_vrr) THEN
2778 CALL dbcsr_allocate_matrix_set(matrix_rrv_vrr, 6)
2779 DO ind = 1, 6
2780 CALL dbcsr_init_p(matrix_rrv_vrr(ind)%matrix)
2781 CALL dbcsr_create(matrix_rrv_vrr(ind)%matrix, template=matrix_s(1)%matrix, &
2782 matrix_type=dbcsr_type_symmetric)
2783 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv_vrr(ind)%matrix, sab_orb)
2784 CALL dbcsr_set(matrix_rrv_vrr(ind)%matrix, 0._dp)
2785 END DO
2786 END IF
2787
2788 ! calculate evaluation of operators in AO basis set
2789 CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
2790 matrix_rxrv=matrix_rxrv, matrix_rrv=matrix_rrv, matrix_rvr=matrix_rvr, &
2791 matrix_rrv_vrr=matrix_rrv_vrr, ref_point=ref_point)
2792
2793 ! Calculate expectation values
2794 ! Real part
2795 NULLIFY (tmp_ao)
2796 CALL dbcsr_init_p(tmp_ao)
2797 CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
2798 CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
2799 CALL dbcsr_set(tmp_ao, 0.0_dp)
2800
2801 IF (calc_rvr .OR. calc_rrv_vrr) THEN
2802 NULLIFY (rho_ao)
2803 CALL qs_rho_get(rho, rho_ao=rho_ao)
2804
2805 IF (calc_rvr) THEN
2806 trace = 0._dp
2807 DO ind = 1, SIZE(matrix_rvr)
2808 strace = 0._dp
2809 DO ispin = 1, dft_control%nspins
2810 CALL dbcsr_set(tmp_ao, 0.0_dp)
2811 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rvr(ind)%matrix, &
2812 0.0_dp, tmp_ao)
2813 CALL dbcsr_trace(tmp_ao, trace)
2814 strace = strace + trace
2815 END DO
2816 nlcom_rvr(ind) = nlcom_rvr(ind) + strace
2817 END DO
2818 END IF
2819
2820 IF (calc_rrv_vrr) THEN
2821 trace = 0._dp
2822 DO ind = 1, SIZE(matrix_rrv_vrr)
2823 strace = 0._dp
2824 DO ispin = 1, dft_control%nspins
2825 CALL dbcsr_set(tmp_ao, 0.0_dp)
2826 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv_vrr(ind)%matrix, &
2827 0.0_dp, tmp_ao)
2828 CALL dbcsr_trace(tmp_ao, trace)
2829 strace = strace + trace
2830 END DO
2831 nlcom_rrv_vrr(ind) = nlcom_rrv_vrr(ind) + strace
2832 END DO
2833 END IF
2834 END IF
2835
2836 ! imagninary part of the density matrix
2837 NULLIFY (rho_ao)
2838 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2839
2840 IF (calc_rv) THEN
2841 trace = 0._dp
2842 DO ind = 1, SIZE(matrix_rv)
2843 strace = 0._dp
2844 DO ispin = 1, dft_control%nspins
2845 CALL dbcsr_set(tmp_ao, 0.0_dp)
2846 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rv(ind)%matrix, &
2847 0.0_dp, tmp_ao)
2848 CALL dbcsr_trace(tmp_ao, trace)
2849 strace = strace + trace
2850 END DO
2851 nlcom_rv(ind) = nlcom_rv(ind) + strace
2852 END DO
2853 END IF
2854
2855 IF (calc_rrv) THEN
2856 trace = 0._dp
2857 DO ind = 1, SIZE(matrix_rrv)
2858 strace = 0._dp
2859 DO ispin = 1, dft_control%nspins
2860 CALL dbcsr_set(tmp_ao, 0.0_dp)
2861 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv(ind)%matrix, &
2862 0.0_dp, tmp_ao)
2863 CALL dbcsr_trace(tmp_ao, trace)
2864 strace = strace + trace
2865 END DO
2866 nlcom_rrv(ind) = nlcom_rrv(ind) + strace
2867 END DO
2868 END IF
2869
2870 IF (calc_rxrv) THEN
2871 trace = 0._dp
2872 DO ind = 1, SIZE(matrix_rxrv)
2873 strace = 0._dp
2874 DO ispin = 1, dft_control%nspins
2875 CALL dbcsr_set(tmp_ao, 0.0_dp)
2876 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rxrv(ind)%matrix, &
2877 0.0_dp, tmp_ao)
2878 CALL dbcsr_trace(tmp_ao, trace)
2879 strace = strace + trace
2880 END DO
2881 nlcom_rxrv(ind) = nlcom_rxrv(ind) + strace
2882 END DO
2883 END IF
2884 CALL dbcsr_deallocate_matrix(tmp_ao)
2885 IF (calc_rv) CALL dbcsr_deallocate_matrix_set(matrix_rv)
2886 IF (calc_rxrv) CALL dbcsr_deallocate_matrix_set(matrix_rxrv)
2887 IF (calc_rrv) CALL dbcsr_deallocate_matrix_set(matrix_rrv)
2888 IF (calc_rvr) CALL dbcsr_deallocate_matrix_set(matrix_rvr)
2889 IF (calc_rrv_vrr) CALL dbcsr_deallocate_matrix_set(matrix_rrv_vrr)
2890
2891 CALL timestop(handle)
2892 END SUBROUTINE calculate_commutator_nl_terms
2893
2894! **************************************************************************************************
2895!> \brief Get list of kpoints from input to compute dipole moment elements
2896!> \param qs_env ...
2897!> \param xkp ...
2898!> \param special_pnts ...
2899!> \author Shridhar Shanbhag
2900! **************************************************************************************************
2901 SUBROUTINE get_xkp_for_dipole_calc(qs_env, xkp, special_pnts)
2902 TYPE(qs_environment_type), POINTER :: qs_env
2903 TYPE(section_vals_type), POINTER :: kpnts, kpset
2904 REAL(kind=dp), DIMENSION(3, 3) :: cart_hmat, hmat
2905
2906 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_xkp_for_dipole_calc'
2907
2908 CHARACTER(LEN=default_string_length) :: ustr
2909 TYPE(kpoint_type), POINTER :: kpoint_work
2910 TYPE(cell_type), POINTER :: cell
2911 CHARACTER(LEN=default_string_length), &
2912 DIMENSION(:), POINTER :: strptr
2913 CHARACTER(LEN=default_string_length), &
2914 DIMENSION(:), POINTER :: special_pnts, spname
2915 CHARACTER(LEN=max_line_length) :: error_message
2916 INTEGER :: handle, i, ik, ikk, ip, &
2917 n_ptr, npline, nkp
2918 LOGICAL :: explicit_kpnts, explicit_kpset
2919 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: kspecial, xkp
2920 REAL(kind=dp), DIMENSION(3) :: kpptr
2921 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2922
2923 CALL timeset(routinen, handle)
2924 kpset => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINT_SET")
2925 kpnts => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINTS")
2926 CALL section_vals_get(kpset, explicit=explicit_kpset)
2927 CALL section_vals_get(kpnts, explicit=explicit_kpnts)
2928 IF (explicit_kpset .AND. explicit_kpnts) then
2929 cpabort("Both KPOINT_SET and KPOINTS present in MOMENTS section")
2930 end if
2931
2932 IF (explicit_kpset) THEN
2933 CALL get_qs_env(qs_env, cell=cell)
2934 CALL get_cell(cell, h=hmat)
2935 cart_hmat(:, :) = hmat(:, :)
2936 IF (cell%input_cell_canonicalized) cart_hmat(:, :) = cell%input_hmat(:, :)
2937 CALL section_vals_val_get(kpset, "NPOINTS", i_val=npline)
2938 CALL section_vals_val_get(kpset, "UNITS", c_val=ustr)
2939 CALL uppercase(ustr)
2940 CALL section_vals_val_get(kpset, "SPECIAL_POINT", n_rep_val=n_ptr)
2941 cpassert(n_ptr > 0)
2942 ALLOCATE (kspecial(3, n_ptr))
2943 ALLOCATE (spname(n_ptr))
2944 DO ip = 1, n_ptr
2945 CALL section_vals_val_get(kpset, "SPECIAL_POINT", i_rep_val=ip, c_vals=strptr)
2946 IF (SIZE(strptr(:), 1) == 4) THEN
2947 spname(ip) = strptr(1)
2948 DO i = 1, 3
2949 CALL read_float_object(strptr(i + 1), kpptr(i), error_message)
2950 IF (len_trim(error_message) > 0) cpabort(trim(error_message))
2951 END DO
2952 ELSE IF (SIZE(strptr(:), 1) == 3) THEN
2953 spname(ip) = "not specified"
2954 DO i = 1, 3
2955 CALL read_float_object(strptr(i), kpptr(i), error_message)
2956 IF (len_trim(error_message) > 0) cpabort(trim(error_message))
2957 END DO
2958 ELSE
2959 cpabort("Input SPECIAL_POINT invalid")
2960 END IF
2961 SELECT CASE (ustr)
2962 CASE ("B_VECTOR")
2963 kspecial(1:3, ip) = kpptr(1:3)
2964 CASE ("CART_ANGSTROM")
2965 kspecial(1:3, ip) = (kpptr(1)*cart_hmat(1, 1:3) + &
2966 kpptr(2)*cart_hmat(2, 1:3) + &
2967 kpptr(3)*cart_hmat(3, 1:3))/twopi*angstrom
2968 CASE ("CART_BOHR")
2969 kspecial(1:3, ip) = (kpptr(1)*cart_hmat(1, 1:3) + &
2970 kpptr(2)*cart_hmat(2, 1:3) + &
2971 kpptr(3)*cart_hmat(3, 1:3))/twopi
2972 CASE DEFAULT
2973 cpabort("Unknown unit <"//trim(ustr)//"> specified for k-point definition")
2974 END SELECT
2975 END DO
2976 nkp = (n_ptr - 1)*npline + 1
2977 cpassert(nkp >= 1)
2978
2979 ! Initialize environment and calculate MOs
2980 ALLOCATE (xkp(3, nkp))
2981 ALLOCATE (special_pnts(nkp))
2982 special_pnts(:) = ""
2983 xkp(1:3, 1) = kspecial(1:3, 1)
2984 ikk = 1
2985 special_pnts(ikk) = spname(1)
2986 DO ik = 2, n_ptr
2987 DO ip = 1, npline
2988 ikk = ikk + 1
2989 xkp(1:3, ikk) = kspecial(1:3, ik - 1) + &
2990 REAL(ip, kind=dp)/real(npline, kind=dp)* &
2991 (kspecial(1:3, ik) - kspecial(1:3, ik - 1))
2992 END DO
2993 special_pnts(ikk) = spname(ik)
2994 END DO
2995 DEALLOCATE (spname, kspecial)
2996 ELSE IF (explicit_kpnts) THEN
2997 CALL get_qs_env(qs_env, particle_set=particle_set, cell=cell)
2998 CALL get_cell(cell, h=hmat)
2999 NULLIFY (kpoint_work)
3000 CALL kpoint_create(kpoint_work)
3001 CALL read_kpoint_section(kpoint_work, kpnts, hmat, cell)
3002 CALL kpoint_initialize(kpoint_work, particle_set, cell)
3003 nkp = kpoint_work%nkp
3004 ALLOCATE (xkp(3, nkp))
3005 ALLOCATE (special_pnts(nkp))
3006 special_pnts(:) = ""
3007 xkp(1:3, :) = kpoint_work%xkp(1:3, :)
3008 CALL kpoint_release(kpoint_work)
3009 ELSE
3010 ! use k-point mesh from DFT calculation
3011 CALL get_qs_env(qs_env, kpoints=kpoint_work)
3012 nkp = kpoint_work%nkp
3013 nkp = kpoint_work%nkp
3014 ALLOCATE (xkp(3, nkp))
3015 ALLOCATE (special_pnts(nkp))
3016 special_pnts(:) = ""
3017 xkp(1:3, :) = kpoint_work%xkp(1:3, :)
3018 END IF
3019 CALL timestop(handle)
3020
3021 END SUBROUTINE get_xkp_for_dipole_calc
3022! **************************************************************************************************
3023!> \brief Calculate local moment matrix for a periodic system for all image cells
3024!> \param qs_env ...
3025!> \param moments_rs_img ...
3026!> \param rcc ...
3027!> \author Shridhar Shanbhag
3028! **************************************************************************************************
3029 SUBROUTINE build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc)
3030
3031 TYPE(qs_environment_type), POINTER :: qs_env
3032 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_rs_img
3033 REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
3034
3035 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moment_matrix_rs_img'
3036
3037 INTEGER :: handle, i_dir, iatom, ic, ikind, iset, j, jatom, jkind, jset, &
3038 ldsa, ldsb, ldwork, ncoa, ncob, nimg, nkind, nseta, nsetb, nsize, sgfa, sgfb
3039 INTEGER, DIMENSION(3) :: icell
3040 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
3041 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3042 LOGICAL :: found
3043 REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
3044 REAL(dp), DIMENSION(:, :), POINTER :: dblock, work
3045 REAL(dp), DIMENSION(:, :, :), POINTER :: dipab
3046 REAL(kind=dp) :: dab
3047 TYPE(cell_type), POINTER :: cell
3048 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
3049 TYPE(dft_control_type), POINTER :: dft_control
3050 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
3051 TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
3052 TYPE(kpoint_type), POINTER :: kpoints_all
3053 TYPE(mp_para_env_type), POINTER :: para_env
3054 TYPE(neighbor_list_iterator_p_type), &
3055 DIMENSION(:), POINTER :: nl_iterator
3056 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3057 POINTER :: sab_all
3058 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3059 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3060 TYPE(qs_kind_type), POINTER :: qs_kind
3061
3062 CALL timeset(routinen, handle)
3063
3064 CALL get_qs_env(qs_env=qs_env, &
3065 dft_control=dft_control, &
3066 qs_kind_set=qs_kind_set, &
3067 matrix_ks_kp=matrix_ks_kp, &
3068 particle_set=particle_set, &
3069 cell=cell, &
3070 para_env=para_env, &
3071 sab_all=sab_all)
3072
3073 NULLIFY (kpoints_all)
3074 CALL kpoint_create(kpoints_all)
3075 CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, nimg)
3076
3077 nkind = SIZE(qs_kind_set)
3078 ALLOCATE (basis_set_list(nkind))
3079 DO ikind = 1, nkind
3080 qs_kind => qs_kind_set(ikind)
3081 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set)
3082 IF (ASSOCIATED(basis_set)) THEN
3083 basis_set_list(ikind)%gto_basis_set => basis_set
3084 ELSE
3085 NULLIFY (basis_set_list(ikind)%gto_basis_set)
3086 END IF
3087 END DO
3088
3089 rc(:) = 0._dp
3090 IF (PRESENT(rcc)) rc(:) = rcc(:)
3091
3092 CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_list)
3093 CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
3094 nsize = SIZE(index_to_cell, 2)
3095 cpassert(SIZE(moments_rs_img, 2) == nsize)
3096 DO i_dir = 1, 3
3097 DO j = 1, nsize
3098 ALLOCATE (moments_rs_img(i_dir, j)%matrix)
3099 CALL dbcsr_create(matrix=moments_rs_img(i_dir, j)%matrix, &
3100 template=matrix_ks_kp(1, 1)%matrix, &
3101 matrix_type=dbcsr_type_no_symmetry, &
3102 name="DIPMAT")
3103 CALL cp_dbcsr_alloc_block_from_nbl(moments_rs_img(i_dir, j)%matrix, sab_all)
3104 CALL dbcsr_set(moments_rs_img(i_dir, j)%matrix, 0.0_dp)
3105 END DO
3106 END DO
3107
3108 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
3109 ALLOCATE (dipab(ldwork, ldwork, 3))
3110 ALLOCATE (work(ldwork, ldwork))
3111
3112 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
3113 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3114 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3115 iatom=iatom, jatom=jatom, r=rab, cell=icell)
3116
3117 basis_set_a => basis_set_list(ikind)%gto_basis_set
3118 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
3119 basis_set_b => basis_set_list(jkind)%gto_basis_set
3120 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
3121 associate( &
3122 ! basis ikind
3123 first_sgfa => basis_set_a%first_sgf, &
3124 la_max => basis_set_a%lmax, &
3125 la_min => basis_set_a%lmin, &
3126 npgfa => basis_set_a%npgf, &
3127 nsgfa => basis_set_a%nsgf_set, &
3128 rpgfa => basis_set_a%pgf_radius, &
3129 set_radius_a => basis_set_a%set_radius, &
3130 sphi_a => basis_set_a%sphi, &
3131 zeta => basis_set_a%zet, &
3132 ! basis jkind, &
3133 first_sgfb => basis_set_b%first_sgf, &
3134 lb_max => basis_set_b%lmax, &
3135 lb_min => basis_set_b%lmin, &
3136 npgfb => basis_set_b%npgf, &
3137 nsgfb => basis_set_b%nsgf_set, &
3138 rpgfb => basis_set_b%pgf_radius, &
3139 set_radius_b => basis_set_b%set_radius, &
3140 sphi_b => basis_set_b%sphi, &
3141 zetb => basis_set_b%zet)
3142
3143 nseta = basis_set_a%nset
3144 nsetb = basis_set_b%nset
3145
3146 ldsa = SIZE(sphi_a, 1)
3147 ldsb = SIZE(sphi_b, 1)
3148
3149 NULLIFY (dblock)
3150
3151 ra = pbc(particle_set(iatom)%r(:), cell)
3152 rb(:) = ra(:) + rab(:)
3153 rac = ra - rc
3154 rbc = rb - rc
3155 dab = norm2(rab)
3156
3157 ic = cell_to_index(icell(1), icell(2), icell(3))
3158
3159 DO iset = 1, nseta
3160
3161 ncoa = npgfa(iset)*ncoset(la_max(iset))
3162 sgfa = first_sgfa(1, iset)
3163
3164 DO jset = 1, nsetb
3165
3166 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
3167
3168 ncob = npgfb(jset)*ncoset(lb_max(jset))
3169 sgfb = first_sgfb(1, jset)
3170
3171 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
3172 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), 1, &
3173 rac, rbc, dipab)
3174 DO i_dir = 1, 3
3175 CALL dbcsr_get_block_p(matrix=moments_rs_img(i_dir, ic)%matrix, &
3176 row=iatom, col=jatom, block=dblock, found=found)
3177 cpassert(found)
3178 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
3179 1.0_dp, dipab(1, 1, i_dir), ldwork, &
3180 sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
3181
3182 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
3183 1.0_dp, sphi_a(1, sgfa), ldsa, &
3184 work(1, 1), ldwork, 1.0_dp, dblock(sgfa, sgfb), SIZE(dblock, 1))
3185 END DO
3186 END DO
3187 END DO
3188 END associate
3189 END DO
3190 CALL neighbor_list_iterator_release(nl_iterator)
3191 CALL kpoint_release(kpoints_all)
3192 DEALLOCATE (dipab, work, basis_set_list)
3193 CALL timestop(handle)
3194
3195 END SUBROUTINE build_local_moment_matrix_rs_img
3196
3197! **************************************************************************************************
3198!> \brief Image-cell index covering the images of two kpoint environments and their mirrors -R.
3199!> Needed when the SCF neighbor list reaches further than sab_all (e.g. with HFX), so
3200!> that real-space matrices of both kinds can be replicated into one set of arrays.
3201!> \param kp_a first kpoint environment (cell index initialised)
3202!> \param kp_b second kpoint environment (cell index initialised)
3203!> \param index_to_cell union image list, dim (3, n_images); allocated here, caller deallocates
3204!> \param cell_to_index mapping of cell coordinates to the union image index (0 = not present);
3205!> allocated here, caller deallocates
3206!> \author Shridhar Shanbhag
3207! **************************************************************************************************
3208 SUBROUTINE get_cell_index_union(kp_a, kp_b, index_to_cell, cell_to_index)
3209 TYPE(kpoint_type), POINTER :: kp_a, kp_b
3210 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
3211 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3212
3213 INTEGER :: i, ic, isgn, n
3214 INTEGER, DIMENSION(3) :: cell, itm
3215 INTEGER, DIMENSION(:, :), ALLOCATABLE :: list
3216 INTEGER, DIMENSION(:, :), POINTER :: itc_a, itc_b
3217
3218 NULLIFY (itc_a, itc_b)
3219 CALL get_kpoint_info(kp_a, index_to_cell=itc_a)
3220 CALL get_kpoint_info(kp_b, index_to_cell=itc_b)
3221 cpassert(ASSOCIATED(itc_a))
3222 cpassert(ASSOCIATED(itc_b))
3223
3224 DO i = 1, 3
3225 itm(i) = max(maxval(abs(itc_a(i, :))), maxval(abs(itc_b(i, :))))
3226 END DO
3227 ALLOCATE (cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
3228 cell_to_index(:, :, :) = 0
3229 ALLOCATE (list(3, 2*(SIZE(itc_a, 2) + SIZE(itc_b, 2))))
3230
3231 ! the images of kp_a keep their order, so index i of kp_a is index i of the union
3232 n = 0
3233 DO isgn = 1, -1, -2
3234 DO ic = 1, SIZE(itc_a, 2)
3235 cell(:) = isgn*itc_a(:, ic)
3236 IF (cell_to_index(cell(1), cell(2), cell(3)) /= 0) cycle
3237 n = n + 1
3238 cell_to_index(cell(1), cell(2), cell(3)) = n
3239 list(:, n) = cell(:)
3240 END DO
3241 END DO
3242 DO isgn = 1, -1, -2
3243 DO ic = 1, SIZE(itc_b, 2)
3244 cell(:) = isgn*itc_b(:, ic)
3245 IF (cell_to_index(cell(1), cell(2), cell(3)) /= 0) cycle
3246 n = n + 1
3247 cell_to_index(cell(1), cell(2), cell(3)) = n
3248 list(:, n) = cell(:)
3249 END DO
3250 END DO
3251
3252 ALLOCATE (index_to_cell(3, n))
3253 index_to_cell(:, :) = list(:, 1:n)
3254 DEALLOCATE (list)
3255
3256 END SUBROUTINE get_cell_index_union
3257
3258! **************************************************************************************************
3259!> \brief Calculate interband k-point dipole moments d_nm(k) = i<u_n|∇_k|u_m>, Berry curvature,
3260!> momentum and band energies for a periodic system for k-points
3261!> \param qs_env ...
3262!> \param xkp list of kpoints
3263!> \param dipole ...
3264!> \param rcc coordinates about which to calculate the dipole
3265!> \param berry_c berry curvature calculated using Ω^γ_n = -Σ_{m≠n} 2*Im[d^α_nm d^β_mn]
3266!> \param momentum ...
3267!> \param energy ...
3268!> \param do_parallel option to distribute the result in dipole across different MPI ranks
3269!> \author Shridhar Shanbhag
3270! **************************************************************************************************
3271 SUBROUTINE qs_moment_kpoints_deep(qs_env, xkp, dipole, rcc, berry_c, momentum, energy, do_parallel)
3272 TYPE(qs_environment_type), POINTER :: qs_env
3273 LOGICAL, OPTIONAL :: do_parallel
3274 LOGICAL :: my_do_parallel, calc_dipole, calc_bc, &
3275 calc_momentum, calc_energy
3276 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_kpoints_deep'
3277 COMPLEX(KIND=dp) :: phase, tmp_max
3278 COMPLEX(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: c_k, h_k, s_k, d_k, cdc, c_dh_c, &
3279 c_ds_c, dh_dk_i, ds_dk_i
3280 COMPLEX(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dip, mom
3281 COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
3282 ALLOCATABLE, OPTIONAL :: dipole, momentum
3283 INTEGER :: handle, i_dir, ikp, nkp, &
3284 n_img_scf, n_img_all, n_img_mom, nao, &
3285 num_pe, num_copy, mepos, n, m, mu, &
3286 ispin, nspin
3287 INTEGER, DIMENSION(3) :: periodic
3288 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_all, index_to_cell_mom
3289 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_all
3290 REAL(kind=dp), PARAMETER :: eps_degenerate = 1.0e-10_dp
3291 REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
3292 REAL(kind=dp), DIMENSION(3) :: my_rcc
3293 REAL(kind=dp), DIMENSION(3, 3) :: hmat
3294 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: eigenvals
3295 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: bc, xkp
3296 REAL(kind=dp), DIMENSION(:, :, :, :), &
3297 ALLOCATABLE, OPTIONAL :: berry_c
3298 REAL(kind=dp), DIMENSION(:, :, :), &
3299 ALLOCATABLE, OPTIONAL :: energy
3300 REAL(kind=dp), DIMENSION(:, :, :, :), &
3301 ALLOCATABLE :: d_rs, h_rs, s_rs
3302 TYPE(cell_type), POINTER :: cell
3303 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_rs_img, matrix_ks_kp, &
3304 matrix_s_kp
3305 TYPE(dft_control_type), POINTER :: dft_control
3306 TYPE(kpoint_type), POINTER :: kpoints_all, kpoints_scf
3307 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3308 TYPE(mp_para_env_type), POINTER :: para_env
3309 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3310 POINTER :: sab_all
3311
3312 CALL timeset(routinen, handle)
3313 calc_dipole = PRESENT(dipole)
3314 calc_bc = PRESENT(berry_c)
3315 calc_momentum = PRESENT(momentum)
3316 calc_energy = PRESENT(energy)
3317 my_do_parallel = .false.
3318 my_rcc = 0.0_dp
3319 IF (PRESENT(do_parallel)) my_do_parallel = do_parallel
3320 IF (PRESENT(rcc)) my_rcc = rcc
3321
3322 CALL get_qs_env(qs_env, &
3323 matrix_ks_kp=matrix_ks_kp, &
3324 matrix_s_kp=matrix_s_kp, &
3325 sab_all=sab_all, &
3326 cell=cell, &
3327 kpoints=kpoints_scf, &
3328 para_env=para_env, &
3329 dft_control=dft_control, &
3330 mos=mos)
3331
3332 CALL get_mo_set(mo_set=mos(1), nao=nao)
3333 CALL get_cell(cell=cell, h=hmat, periodic=periodic)
3334 nspin = SIZE(matrix_ks_kp, 1)
3335 nkp = SIZE(xkp, 2)
3336
3337 ! create kpoint environment kpoints_all which contains all neighbor cells R of sab_all
3338 ! (the images of the dipole operator) without considering any lattice symmetry
3339 NULLIFY (kpoints_all, index_to_cell_mom)
3340 CALL kpoint_create(kpoints_all)
3341 CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, n_img_scf)
3342 CALL get_kpoint_info(kpoints_all, index_to_cell=index_to_cell_mom)
3343 n_img_mom = SIZE(index_to_cell_mom, 2)
3344
3345 ! The SCF neighbor list (e.g. extended for HFX) can contain images not in sab_all, so the
3346 ! arrays S_rs, H_rs, D_rs are indexed by the union of both image sets (incl. mirrors -R)
3347 CALL get_cell_index_union(kpoints_all, kpoints_scf, index_to_cell_all, cell_to_index_all)
3348 n_img_all = SIZE(index_to_cell_all, 2)
3349
3350 NULLIFY (moments_rs_img)
3351 CALL dbcsr_allocate_matrix_set(moments_rs_img, 3, n_img_mom)
3352 ! D_μ,ν = <φ_μ|r|φ_ν>
3353 CALL build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc=my_rcc)
3354
3355 ALLOCATE (s_rs(1, nao, nao, n_img_all), h_rs(nspin, nao, nao, n_img_all), source=0.0_dp)
3356 ALLOCATE (d_rs(3, nao, nao, n_img_all), source=0.0_dp)
3357
3358 ! Convert real-space dbcsr matrices into arrays
3359 CALL replicate_rs_matrices(matrix_s_kp, kpoints_scf, s_rs, cell_to_index_all)
3360 CALL replicate_rs_matrices(matrix_ks_kp, kpoints_scf, h_rs, cell_to_index_all)
3361 CALL replicate_rs_matrices(moments_rs_img, kpoints_all, d_rs, cell_to_index_all)
3362
3363 mepos = 0
3364 num_pe = 1
3365 num_copy = nkp
3366 IF (my_do_parallel) THEN
3367 mepos = para_env%mepos
3368 num_pe = para_env%num_pe
3369 num_copy = ceiling(real(nkp)/num_pe)
3370 END IF
3371
3372 IF (calc_dipole) ALLOCATE (dipole(nspin, num_copy, 3, nao, nao), source=z_zero)
3373 IF (calc_bc) ALLOCATE (berry_c(nspin, num_copy, 3, nao), source=0.0_dp)
3374 IF (calc_momentum) ALLOCATE (momentum(nspin, num_copy, 3, nao, nao), source=z_zero)
3375 IF (calc_energy) ALLOCATE (energy(nspin, num_copy, nao), source=0.0_dp)
3376
3377!$OMP PARALLEL DEFAULT(NONE) PRIVATE(ikp, S_k, H_k, eigenvals, C_k, ispin, n, m, &
3378!$OMP i_dir, dS_dk_i, dH_dk_i, D_k, dip, bc, mom, C_dS_C, C_dH_C, CDC, tmp_max, phase) &
3379!$OMP SHARED(num_pe, mepos, dipole, berry_c, momentum, energy, nao, nspin, periodic, nkp, &
3380!$OMP xkp, S_rs, H_rs, D_rs, index_to_cell_all, hmat, calc_dipole, calc_bc, calc_momentum, calc_energy)
3381 ALLOCATE (ds_dk_i(nao, nao), c_ds_c(nao, nao), dh_dk_i(nao, nao), c_dh_c(nao, nao), source=z_zero)
3382 ALLOCATE (cdc(nao, nao), s_k(nao, nao), h_k(nao, nao), source=z_zero)
3383 ALLOCATE (c_k(nao, nao), d_k(nao, nao), source=z_zero)
3384 ALLOCATE (eigenvals(nao), source=0.0_dp)
3385 IF (calc_dipole .OR. calc_bc) ALLOCATE (dip(3, nao, nao), source=z_zero)
3386 IF (calc_bc) ALLOCATE (bc(3, nao), source=0.0_dp)
3387 IF (calc_momentum) ALLOCATE (mom(3, nao, nao), source=z_zero)
3388!$OMP DO COLLAPSE(2)
3389 DO ispin = 1, nspin
3390 DO ikp = 1, nkp
3391 IF (mod(ikp - 1, num_pe) /= mepos) cycle
3392
3393 ! S^R -> S(k), H^R -> H(k)
3394 s_k = 0
3395 h_k = 0
3396 CALL rs_to_kp(s_rs(1, :, :, :), s_k, index_to_cell_all, xkp(:, ikp))
3397 CALL rs_to_kp(h_rs(ispin, :, :, :), h_k, index_to_cell_all, xkp(:, ikp))
3398
3399 ! Diagonalize H(k)C(k) = S(k)C(k)ε(k)
3400 CALL geeig_right(h_k, s_k, eigenvals, c_k)
3401
3402 ! To have a smooth complex phase of C(k) as function of k, for every n, we force
3403 ! the largest C_μ,n(k) to be real.
3404 ! This is important to have a continuous dipole moment d_nm(k) as a function of k
3405 DO n = 1, nao
3406 tmp_max = c_k(1, n)
3407 DO mu = 1, nao
3408 IF (abs(c_k(mu, n)) < abs(tmp_max)) cycle
3409 tmp_max = c_k(mu, n)
3410 END DO
3411 phase = tmp_max/abs(tmp_max)
3412 c_k(:, n) = c_k(:, n)/phase
3413 END DO
3414
3415 ! reset per k-point/spin: degenerate pairs are skipped below and must not keep
3416 ! values from a previous iteration
3417 IF (calc_dipole .OR. calc_bc) dip = z_zero
3418 IF (calc_momentum) mom = z_zero
3419
3420 DO i_dir = 1, 3 ! d^x, d^y, d^z
3421
3422 IF (periodic(i_dir) == 0) cycle
3423 ! ∇ S(k) = Σ_R iR S^R e^(ikR), ∇ H(k) = Σ_R iR H^R e^(ikR)
3424 CALL rs_to_kp(s_rs(1, :, :, :), ds_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
3425 CALL rs_to_kp(h_rs(ispin, :, :, :), dh_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
3426
3427 ! Σ_R D^R e^(ikR) = D(k), D_μ,ν = <φ_μ|r|φ_ν>
3428 CALL rs_to_kp(d_rs(i_dir, :, :, :), d_k(:, :), index_to_cell_all, xkp(:, ikp))
3429
3430 ! Basis transform to Kohn-Sham basis: (C^H) ∇ S C, (C^H) ∇ H C, (C^H) D C
3431 CALL gemm_square(c_k, 'C', ds_dk_i, 'N', c_k, 'N', c_ds_c)
3432 CALL gemm_square(c_k, 'C', dh_dk_i, 'N', c_k, 'N', c_dh_c)
3433 CALL gemm_square(c_k, 'C', d_k, 'N', c_k, 'N', cdc)
3434
3435 ! Compute the dipole (n ≠ m)
3436 ! d_nm (k) = - i/(ε_n - ε_m) [ C^† dH(k)/dk C ]_nm
3437 ! + i ε_n/(ε_n - ε_m) [ C^† dS(k)/dk C ]_nm + [ C^†D(k)C ]_nm
3438 ! and the momentum
3439 ! p_nm = <ψ_n|-iħ∇_r|ψ_m> = <u_n|(ħk - iħ∇_r)|u_m>
3440 ! = i d_nm (ε_n - ε_m) + ∇_k ε_n δ_nm
3441 ! = [ C^† dH(k)/dk C ]_nm - ε_n [ C^† dS(k)/dk C ]_nm
3442 ! + i (ε_n - ε_m) [ (C^H)D(k)C ]_nm
3443 DO n = 1, nao
3444 DO m = 1, nao
3445 IF (n == m) THEN
3446 ! diagonal elements of the dipole would need to be computed from
3447 ! a numerical k-derivative which is not implemented
3448 IF (calc_momentum) THEN
3449 mom(i_dir, n, n) = dble(c_dh_c(n, n)) - dble(eigenvals(n)*c_ds_c(n, n))
3450 END IF
3451 ELSE
3452 IF ((calc_dipole .OR. calc_bc) .AND. (abs(eigenvals(n) - eigenvals(m)) >= eps_degenerate)) THEN
3453 dip(i_dir, n, m) = -gaussi*c_dh_c(n, m)/(eigenvals(n) - eigenvals(m)) &
3454 + gaussi*eigenvals(n)*c_ds_c(n, m)/(eigenvals(n) - eigenvals(m)) &
3455 + cdc(n, m)
3456 END IF
3457 IF (calc_momentum) THEN
3458 mom(i_dir, n, m) = c_dh_c(n, m) - eigenvals(n)*c_ds_c(n, m) &
3459 + gaussi*cdc(n, m)*(eigenvals(n) - eigenvals(m))
3460 END IF
3461 END IF
3462 END DO
3463 END DO
3464 END DO
3465 ! Compute the Berry curvature from the dipoles
3466 ! Ω^γ_n = -Σ_{m≠n} 2*Im[d^α_nm d^β_mn], where (α, β, γ) are cyclic in {x, y, z}
3467 IF (calc_bc) THEN
3468 bc = 0.0_dp
3469 DO i_dir = 1, 3
3470 DO n = 1, nao
3471 DO m = 1, nao
3472 IF (n == m) cycle
3473 bc(i_dir, n) = bc(i_dir, n) &
3474 - 2*aimag(dip(1 + mod(i_dir, 3), n, m)*dip(1 + mod(i_dir + 1, 3), m, n))
3475 END DO
3476 END DO
3477 END DO
3478 END IF
3479 ! Store the dipoles, momentum, berry curvature and energies for each MPI rank
3480 IF (calc_dipole) dipole(ispin, ceiling(real(ikp)/num_pe), :, :, :) = dip(:, :, :)
3481 IF (calc_momentum) momentum(ispin, ceiling(real(ikp)/num_pe), :, :, :) = mom(:, :, :)
3482 IF (calc_bc) berry_c(ispin, ceiling(real(ikp)/num_pe), :, :) = bc(:, :)
3483 IF (calc_energy) energy(ispin, ceiling(real(ikp)/num_pe), :) = eigenvals(:)
3484 END DO
3485 END DO
3486!$OMP END DO
3487 DEALLOCATE (ds_dk_i, c_ds_c, dh_dk_i, c_dh_c, cdc, s_k, h_k, c_k, d_k, eigenvals)
3488 IF (calc_dipole .OR. calc_bc) DEALLOCATE (dip)
3489 IF (calc_bc) DEALLOCATE (bc)
3490 IF (calc_momentum) DEALLOCATE (mom)
3491!$OMP END PARALLEL
3492 DEALLOCATE (s_rs, h_rs, d_rs)
3493 DEALLOCATE (index_to_cell_all, cell_to_index_all)
3494 CALL dbcsr_deallocate_matrix_set(moments_rs_img)
3495 CALL kpoint_release(kpoints_all)
3496 CALL timestop(handle)
3497 END SUBROUTINE qs_moment_kpoints_deep
3498
3499! **************************************************************************************************
3500!> \brief Add the spin-orbit coupling operator V^SOC(k) = sigma . V(k) (or its k-derivative) to a
3501!> 2*nao spinor matrix. The real-space components VSOC_rs store V^SOC without the imaginary
3502!> unit (real, antisymmetric), so the physical Hermitian component is W_c = i*V_c(k). The
3503!> Pauli placement mirrors H_KS_spinor_kp in post_scf_bandstructure_utils.F:
3504!> uu += W_z ; dd += -W_z ; ud += W_x + i*W_y ; du += W_x - i*W_y.
3505!> \param mat_spinor 2*nao spinor matrix to be updated in place
3506!> \param VSOC_rs real-space SOC images, dim (3=xyz, nao, nao, n_img)
3507!> \param index_to_cell image-cell coordinates
3508!> \param xkp k-point
3509!> \param nao number of atomic orbitals (block size)
3510!> \param deriv_direction optional cartesian direction for the k-derivative of V^SOC
3511!> \param hmat cell matrix, required for the derivative
3512!> \author Shridhar Shanbhag
3513! **************************************************************************************************
3514 SUBROUTINE add_vsoc_to_spinor(mat_spinor, VSOC_rs, index_to_cell, xkp, nao, &
3515 deriv_direction, hmat)
3516 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: mat_spinor
3517 REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: vsoc_rs
3518 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
3519 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xkp
3520 INTEGER, INTENT(IN) :: nao
3521 INTEGER, INTENT(IN), OPTIONAL :: deriv_direction
3522 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN), &
3523 OPTIONAL :: hmat
3524
3525 COMPLEX(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: vx, vy, vz
3526 INTEGER :: s
3527
3528 s = nao + 1
3529 ALLOCATE (vx(nao, nao), vy(nao, nao), vz(nao, nao), source=z_zero)
3530
3531 IF (PRESENT(deriv_direction)) THEN
3532 CALL rs_to_kp(vsoc_rs(1, :, :, :), vx, index_to_cell, xkp, deriv_direction, hmat)
3533 CALL rs_to_kp(vsoc_rs(2, :, :, :), vy, index_to_cell, xkp, deriv_direction, hmat)
3534 CALL rs_to_kp(vsoc_rs(3, :, :, :), vz, index_to_cell, xkp, deriv_direction, hmat)
3535 ELSE
3536 CALL rs_to_kp(vsoc_rs(1, :, :, :), vx, index_to_cell, xkp)
3537 CALL rs_to_kp(vsoc_rs(2, :, :, :), vy, index_to_cell, xkp)
3538 CALL rs_to_kp(vsoc_rs(3, :, :, :), vz, index_to_cell, xkp)
3539 END IF
3540
3541 ! physical component operators W_c = i * V_c (VSOC_rs stores V without the imaginary unit)
3542 vx(:, :) = gaussi*vx(:, :)
3543 vy(:, :) = gaussi*vy(:, :)
3544 vz(:, :) = gaussi*vz(:, :)
3545
3546 mat_spinor(1:nao, 1:nao) = mat_spinor(1:nao, 1:nao) + vz(:, :)
3547 mat_spinor(s:2*nao, s:2*nao) = mat_spinor(s:2*nao, s:2*nao) - vz(:, :)
3548 mat_spinor(1:nao, s:2*nao) = mat_spinor(1:nao, s:2*nao) + vx(:, :) + gaussi*vy(:, :)
3549 mat_spinor(s:2*nao, 1:nao) = mat_spinor(s:2*nao, 1:nao) + vx(:, :) - gaussi*vy(:, :)
3550
3551 DEALLOCATE (vx, vy, vz)
3552
3553 END SUBROUTINE add_vsoc_to_spinor
3554
3555! **************************************************************************************************
3556!> \brief Calculate interband k-point dipole moments d_nm(k) = i<u_n|∇_k|u_m>, Berry curvature,
3557!> momentum and band energies for a periodic system for k-points with spin-orbit coupling
3558!> (SOC) in the 2*nao spinor basis, via the length-gauge sum rule.
3559!> SOC enters both the Hamiltonian and the velocity operator (the k-derivative of V^SOC is
3560!> included). Mirrors qs_moment_kpoints_deep but promoted to spinor space; the two spin
3561!> channels are merged into 2*nao spinor bands (returned with spin dimension 1).
3562!> \param qs_env ...
3563!> \param xkp k-points at which the dipoles are evaluated
3564!> \param dipole optional output dipoles, dim (1, k copy, 3, 2*nao, 2*nao)
3565!> \param rcc reference point for the local dipole operator (interband dipoles are origin-independent)
3566!> \param berry_c optional SOC Berry curvature, dim (1, k copy, 3, 2*nao)
3567!> \param momentum optional SOC momentum matrix elements p_nm, dim (1, k copy, 3, 2*nao, 2*nao)
3568!> \param energy optional spinor band energies, dim (1, k copy, 2*nao)
3569!> \param do_parallel distribute k-points over MPI ranks
3570!> \author Shridhar Shanbhag
3571! **************************************************************************************************
3572 SUBROUTINE qs_moment_kpoints_deep_soc(qs_env, xkp, dipole, rcc, berry_c, momentum, energy, do_parallel)
3573 TYPE(qs_environment_type), POINTER :: qs_env
3574 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: xkp
3575 COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
3576 ALLOCATABLE, OPTIONAL :: dipole
3577 REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
3578 REAL(kind=dp), DIMENSION(:, :, :, :), &
3579 ALLOCATABLE, OPTIONAL :: berry_c
3580 COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
3581 ALLOCATABLE, OPTIONAL :: momentum
3582 REAL(kind=dp), DIMENSION(:, :, :), &
3583 ALLOCATABLE, OPTIONAL :: energy
3584 LOGICAL, OPTIONAL :: do_parallel
3585
3586 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_kpoints_deep_soc'
3587 REAL(kind=dp), PARAMETER :: eps_degenerate = 1.0e-10_dp
3588
3589 COMPLEX(KIND=dp) :: phase, tmp_max
3590 COMPLEX(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: d_k, dhd_k, dhu_k, ds_k, hd_k, hu_k, &
3591 s_k, d2, dh2, ds2, h2, s2, w, wd2w, &
3592 wdh2w, wds2w
3593 COMPLEX(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dip, mom
3594 INTEGER :: handle, i_dir, ikp, mepos, mu, m, n, &
3595 n_img_all, n_img_mom, nao, nao2, nimg, nkp, &
3596 nspin, num_copy, num_pe, s
3597 INTEGER, DIMENSION(3) :: periodic
3598 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_all, index_to_cell_mom
3599 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_all
3600 LOGICAL :: calc_bc, calc_dipole, calc_energy, &
3601 calc_momentum, my_do_parallel
3602 REAL(kind=dp), DIMENSION(3) :: my_rcc
3603 REAL(kind=dp), DIMENSION(3, 3) :: hmat
3604 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: eigenvals
3605 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: bc
3606 REAL(kind=dp), DIMENSION(:, :, :, :), ALLOCATABLE :: d_rs, h_rs, s_rs, vsoc_rs
3607 TYPE(cell_type), POINTER :: cell
3608 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_v_soc_xyz, matrix_ks_kp, &
3609 matrix_s_kp, moments_rs_img
3610 TYPE(dft_control_type), POINTER :: dft_control
3611 TYPE(kpoint_type), POINTER :: kpoints_all, kpoints_scf
3612 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3613 TYPE(mp_para_env_type), POINTER :: para_env
3614 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3615 POINTER :: sab_all
3616
3617 CALL timeset(routinen, handle)
3618 calc_dipole = PRESENT(dipole)
3619 calc_bc = PRESENT(berry_c)
3620 calc_momentum = PRESENT(momentum)
3621 calc_energy = PRESENT(energy)
3622 my_do_parallel = .false.
3623 my_rcc = 0.0_dp
3624 IF (PRESENT(do_parallel)) my_do_parallel = do_parallel
3625 IF (PRESENT(rcc)) my_rcc = rcc
3626
3627 NULLIFY (cell, cell_to_index_all, dft_control, index_to_cell_all, index_to_cell_mom, kpoints_all, &
3628 kpoints_scf, mat_v_soc_xyz, matrix_ks_kp, matrix_s_kp, moments_rs_img, mos, &
3629 para_env, sab_all)
3630
3631 CALL get_qs_env(qs_env, &
3632 matrix_ks_kp=matrix_ks_kp, &
3633 matrix_s_kp=matrix_s_kp, &
3634 sab_all=sab_all, &
3635 cell=cell, &
3636 kpoints=kpoints_scf, &
3637 para_env=para_env, &
3638 dft_control=dft_control, &
3639 mos=mos)
3640
3641 CALL get_mo_set(mo_set=mos(1), nao=nao)
3642 CALL get_cell(cell=cell, h=hmat, periodic=periodic)
3643 nspin = SIZE(matrix_ks_kp, 1)
3644 nkp = SIZE(xkp, 2)
3645 nao2 = 2*nao
3646 s = nao + 1
3647
3648 ! create kpoint environment kpoints_all which contains all neighbor cells R of sab_all
3649 ! (the images of the dipole operator)
3650 CALL kpoint_create(kpoints_all)
3651 CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, nimg)
3652 CALL get_kpoint_info(kpoints_all, index_to_cell=index_to_cell_mom)
3653 n_img_mom = SIZE(index_to_cell_mom, 2)
3654
3655 ! The SCF neighbor list (e.g. extended for HFX) can contain images not in sab_all, so the
3656 ! arrays S_rs, H_rs, D_rs, VSOC_rs are indexed by the union of both image sets (incl. -R)
3657 CALL get_cell_index_union(kpoints_all, kpoints_scf, index_to_cell_all, cell_to_index_all)
3658 n_img_all = SIZE(index_to_cell_all, 2)
3659
3660 NULLIFY (moments_rs_img)
3661 CALL dbcsr_allocate_matrix_set(moments_rs_img, 3, n_img_mom)
3662 ! D_mu,nu = <phi_mu|r|phi_nu>
3663 CALL build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc=my_rcc)
3664
3665 ! Real-space SOC operator V^SOC_c^R (c=x,y,z), stored real and antisymmetric on sab_orb.
3666 ! Aborts inside if the atoms carry no pseudopotential SOC parameters.
3667 CALL v_soc_xyz_from_pseudopotential(qs_env, mat_v_soc_xyz)
3668
3669 ALLOCATE (s_rs(1, nao, nao, n_img_all), h_rs(nspin, nao, nao, n_img_all), source=0.0_dp)
3670 ALLOCATE (d_rs(3, nao, nao, n_img_all), vsoc_rs(3, nao, nao, n_img_all), source=0.0_dp)
3671
3672 ! Convert real-space dbcsr matrices into arrays (S/H on the SCF kpoints, D on kpoints_all)
3673 CALL replicate_rs_matrices(matrix_s_kp, kpoints_scf, s_rs, cell_to_index_all)
3674 CALL replicate_rs_matrices(matrix_ks_kp, kpoints_scf, h_rs, cell_to_index_all)
3675 CALL replicate_rs_matrices(moments_rs_img, kpoints_all, d_rs, cell_to_index_all)
3676 CALL replicate_rs_matrices(mat_v_soc_xyz, kpoints_scf, vsoc_rs, cell_to_index_all)
3677
3678 mepos = 0
3679 num_pe = 1
3680 num_copy = nkp
3681 IF (my_do_parallel) THEN
3682 mepos = para_env%mepos
3683 num_pe = para_env%num_pe
3684 num_copy = ceiling(real(nkp)/num_pe)
3685 END IF
3686
3687 IF (calc_dipole) ALLOCATE (dipole(1, num_copy, 3, nao2, nao2), source=z_zero)
3688 IF (calc_bc) ALLOCATE (berry_c(1, num_copy, 3, nao2), source=0.0_dp)
3689 IF (calc_momentum) ALLOCATE (momentum(1, num_copy, 3, nao2, nao2), source=z_zero)
3690 IF (calc_energy) ALLOCATE (energy(1, num_copy, nao2), source=0.0_dp)
3691
3692!$OMP PARALLEL DEFAULT(NONE) PRIVATE(ikp, i_dir, n, m, mu, tmp_max, phase, eigenvals, &
3693!$OMP S_k, Hu_k, Hd_k, D_k, dS_k, dHu_k, dHd_k, S2, H2, D2, dS2, dH2, W, &
3694!$OMP WdH2W, WdS2W, WD2W, dip, bc, mom) &
3695!$OMP SHARED(num_pe, mepos, dipole, berry_c, momentum, energy, nao, nao2, s, nspin, periodic, &
3696!$OMP nkp, xkp, S_rs, H_rs, D_rs, VSOC_rs, index_to_cell_all, hmat, calc_dipole, calc_bc, &
3697!$OMP calc_momentum, calc_energy)
3698 ALLOCATE (s_k(nao, nao), hu_k(nao, nao), hd_k(nao, nao), d_k(nao, nao), source=z_zero)
3699 ALLOCATE (ds_k(nao, nao), dhu_k(nao, nao), dhd_k(nao, nao), source=z_zero)
3700 ALLOCATE (s2(nao2, nao2), h2(nao2, nao2), d2(nao2, nao2), ds2(nao2, nao2), &
3701 dh2(nao2, nao2), w(nao2, nao2), source=z_zero)
3702 ALLOCATE (wdh2w(nao2, nao2), wds2w(nao2, nao2), wd2w(nao2, nao2), source=z_zero)
3703 ALLOCATE (eigenvals(nao2), source=0.0_dp)
3704 IF (calc_dipole .OR. calc_bc) ALLOCATE (dip(3, nao2, nao2), source=z_zero)
3705 IF (calc_bc) ALLOCATE (bc(3, nao2), source=0.0_dp)
3706 IF (calc_momentum) ALLOCATE (mom(3, nao2, nao2), source=z_zero)
3707!$OMP DO
3708 DO ikp = 1, nkp
3709 IF (mod(ikp - 1, num_pe) /= mepos) cycle
3710
3711 ! S(k), H_up(k), H_dn(k)
3712 CALL rs_to_kp(s_rs(1, :, :, :), s_k, index_to_cell_all, xkp(:, ikp))
3713 CALL rs_to_kp(h_rs(1, :, :, :), hu_k, index_to_cell_all, xkp(:, ikp))
3714 CALL rs_to_kp(h_rs(nspin, :, :, :), hd_k, index_to_cell_all, xkp(:, ikp))
3715
3716 ! spinor overlap S2 = blockdiag(S, S)
3717 s2 = z_zero
3718 s2(1:nao, 1:nao) = s_k(:, :)
3719 s2(s:nao2, s:nao2) = s_k(:, :)
3720
3721 ! spinor Hamiltonian H2 = blockdiag(H_up, H_dn) + V^SOC(k)
3722 h2 = z_zero
3723 h2(1:nao, 1:nao) = hu_k(:, :)
3724 h2(s:nao2, s:nao2) = hd_k(:, :)
3725 CALL add_vsoc_to_spinor(h2, vsoc_rs, index_to_cell_all, xkp(:, ikp), nao)
3726
3727 ! Diagonalize H2 W = S2 W E -> spinor bands and AO-spinor eigenvectors
3728 CALL geeig_right(h2, s2, eigenvals, w)
3729
3730 ! To have a smooth complex phase of C(k) as function of k, for every n, we force
3731 ! the largest C_μ,n(k) to be real.
3732 ! This is important to have a continuous dipole moment d_nm(k) as a function of k
3733 DO n = 1, nao2
3734 tmp_max = w(1, n)
3735 DO mu = 1, nao2
3736 IF (abs(w(mu, n)) < abs(tmp_max)) cycle
3737 tmp_max = w(mu, n)
3738 END DO
3739 phase = tmp_max/abs(tmp_max)
3740 w(:, n) = w(:, n)/phase
3741 END DO
3742
3743 IF (calc_dipole .OR. calc_bc) dip = z_zero
3744 IF (calc_momentum) mom = z_zero
3745 DO i_dir = 1, 3 ! d^x, d^y, d^z
3746
3747 IF (periodic(i_dir) == 0) cycle
3748
3749 ! spinor overlap derivative dS2 = blockdiag(grad S, grad S)
3750 CALL rs_to_kp(s_rs(1, :, :, :), ds_k, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
3751 ds2 = z_zero
3752 ds2(1:nao, 1:nao) = ds_k(:, :)
3753 ds2(s:nao2, s:nao2) = ds_k(:, :)
3754
3755 ! spinor Hamiltonian derivative dH2 = blockdiag(grad H_up, grad H_dn) + grad V^SOC
3756 CALL rs_to_kp(h_rs(1, :, :, :), dhu_k, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
3757 CALL rs_to_kp(h_rs(nspin, :, :, :), dhd_k, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
3758 dh2 = z_zero
3759 dh2(1:nao, 1:nao) = dhu_k(:, :)
3760 dh2(s:nao2, s:nao2) = dhd_k(:, :)
3761 CALL add_vsoc_to_spinor(dh2, vsoc_rs, index_to_cell_all, xkp(:, ikp), nao, i_dir, hmat)
3762
3763 ! spinor dipole operator D2 = blockdiag(D(k), D(k))
3764 CALL rs_to_kp(d_rs(i_dir, :, :, :), d_k, index_to_cell_all, xkp(:, ikp))
3765 d2 = z_zero
3766 d2(1:nao, 1:nao) = d_k(:, :)
3767 d2(s:nao2, s:nao2) = d_k(:, :)
3768
3769 ! Basis transform to the spinor eigenbasis: W^H (grad S) W, W^H (grad H) W, W^H D W
3770 CALL gemm_square(w, 'C', dh2, 'N', w, 'N', wdh2w)
3771 CALL gemm_square(w, 'C', ds2, 'N', w, 'N', wds2w)
3772 CALL gemm_square(w, 'C', d2, 'N', w, 'N', wd2w)
3773
3774 ! d_nm(k) = -i/(E_n-E_m)[W^H grad H W]_nm + i E_n/(E_n-E_m)[W^H grad S W]_nm + [W^H D W]_nm
3775 ! Degenerate spinor pairs (|E_n-E_m| < eps_degenerate, e.g. Kramers) are skipped for the
3776 ! dipole (and hence the Berry curvature): the length-gauge sum rule needs degenerate
3777 ! perturbation theory there. The momentum below is well defined and is still computed.
3778 ! Momentum (grad H includes grad V^SOC, i.e. the SOC part of the velocity):
3779 ! p_nm = i d_nm (E_n - E_m) + grad_k E_n delta_nm
3780 ! = [W^H grad H W]_nm - E_n [W^H grad S W]_nm + i (E_n - E_m) [W^H D W]_nm
3781 DO n = 1, nao2
3782 DO m = 1, nao2
3783 IF (n == m) THEN
3784 ! diagonal elements of the dipole would need to be computed from
3785 ! a numerical k-derivative which is not implemented
3786 IF (calc_momentum) THEN
3787 mom(i_dir, n, n) = dble(wdh2w(n, n)) - dble(eigenvals(n)*wds2w(n, n))
3788 END IF
3789 ELSE
3790 IF ((calc_dipole .OR. calc_bc) .AND. (abs(eigenvals(n) - eigenvals(m)) >= eps_degenerate)) THEN
3791 dip(i_dir, n, m) = -gaussi*wdh2w(n, m)/(eigenvals(n) - eigenvals(m)) &
3792 + gaussi*eigenvals(n)*wds2w(n, m)/(eigenvals(n) - eigenvals(m)) &
3793 + wd2w(n, m)
3794 END IF
3795 IF (calc_momentum) THEN
3796 mom(i_dir, n, m) = wdh2w(n, m) - eigenvals(n)*wds2w(n, m) &
3797 + gaussi*wd2w(n, m)*(eigenvals(n) - eigenvals(m))
3798 END IF
3799 END IF
3800 END DO
3801 END DO
3802 END DO
3803
3804 ! Berry curvature Omega^gamma_n = -Sum_{m/=n} 2*Im[d^alpha_nm d^beta_mn], (alpha,beta,gamma) cyclic
3805 IF (calc_bc) THEN
3806 bc = 0.0_dp
3807 DO i_dir = 1, 3
3808 DO n = 1, nao2
3809 DO m = 1, nao2
3810 IF (n == m) cycle
3811 bc(i_dir, n) = bc(i_dir, n) &
3812 - 2*aimag(dip(1 + mod(i_dir, 3), n, m)*dip(1 + mod(i_dir + 1, 3), m, n))
3813 END DO
3814 END DO
3815 END DO
3816 END IF
3817
3818 ! Store the dipoles, momentum, berry curvature and energies for each MPI rank
3819 IF (calc_dipole) dipole(1, ceiling(real(ikp)/num_pe), :, :, :) = dip(:, :, :)
3820 IF (calc_momentum) momentum(1, ceiling(real(ikp)/num_pe), :, :, :) = mom(:, :, :)
3821 IF (calc_bc) berry_c(1, ceiling(real(ikp)/num_pe), :, :) = bc(:, :)
3822 IF (calc_energy) energy(1, ceiling(real(ikp)/num_pe), :) = eigenvals(:)
3823 END DO
3824!$OMP END DO
3825 DEALLOCATE (s_k, hu_k, hd_k, d_k, ds_k, dhu_k, dhd_k, s2, h2, d2, ds2, dh2, w, &
3826 wdh2w, wds2w, wd2w, eigenvals)
3827 IF (calc_dipole .OR. calc_bc) DEALLOCATE (dip)
3828 IF (calc_bc) DEALLOCATE (bc)
3829 IF (calc_momentum) DEALLOCATE (mom)
3830!$OMP END PARALLEL
3831 DEALLOCATE (s_rs, h_rs, d_rs, vsoc_rs)
3832 DEALLOCATE (index_to_cell_all, cell_to_index_all)
3833 CALL dbcsr_deallocate_matrix_set(moments_rs_img)
3834 CALL dbcsr_deallocate_matrix_set(mat_v_soc_xyz)
3835 CALL kpoint_release(kpoints_all)
3836 CALL timestop(handle)
3837 END SUBROUTINE qs_moment_kpoints_deep_soc
3838
3839! **************************************************************************************************
3840!> \brief Calculates interband k-point dipoles in the existing SCF MO basis.
3841!> \param qs_env ...
3842!> \param dipole ...
3843!> \param rcc retained for interface compatibility; interband dipoles are origin independent
3844!> \param nmo_spin_out number of SCF MOs available for each spin
3845! **************************************************************************************************
3846 SUBROUTINE qs_moment_kpoints_scf_mos(qs_env, dipole, rcc, nmo_spin_out)
3847 TYPE(qs_environment_type), POINTER :: qs_env
3848 COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
3849 ALLOCATABLE :: dipole
3850 REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
3851 INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT), &
3852 OPTIONAL :: nmo_spin_out
3853
3854 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_kpoints_scf_mos'
3855
3856 INTEGER :: handle, i_dir, ikp, ikp_local, ispin, &
3857 m, n, nao, nkp, nmo, nspin
3858 INTEGER, DIMENSION(:), ALLOCATABLE :: nmo_spin
3859 INTEGER, DIMENSION(2) :: kp_range
3860 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3861 LOGICAL :: my_kpgrp
3862 REAL(kind=dp), PARAMETER :: eps_degenerate = 1.0e-10_dp
3863 COMPLEX(KIND=dp) :: moment_element
3864 REAL(kind=dp) :: energy_diff
3865 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: eigenvalues_kp
3866 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvals
3867 TYPE(cp_blacs_env_type), POINTER :: blacs_env_all
3868 TYPE(cp_fm_struct_type), POINTER :: fm_struct
3869 TYPE(cp_cfm_type) :: cfm_empty, cfm_mo_coeff_global, &
3870 cfm_moment, cfm_operator, cfm_work
3871 TYPE(cp_fm_type) :: fm_dummy, fm_mo_coeff_global, &
3872 fm_operator_imag, fm_operator_real
3873 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: overlap_deriv
3874 TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix
3875 TYPE(dft_control_type), POINTER :: dft_control
3876 TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env
3877 TYPE(kpoint_env_type), POINTER :: kp
3878 TYPE(kpoint_type), POINTER :: kpoints_scf
3879 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_kp
3880 TYPE(cp_cfm_type), POINTER :: cmo_coeff_source
3881 TYPE(cp_fm_type), POINTER :: fm_mo_coeff_source
3882 TYPE(mp_para_env_type), POINTER :: para_env, para_env_kp
3883 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3884 POINTER :: sab_kp, sab_orb
3885 TYPE(qs_ks_env_type), POINTER :: ks_env
3886
3887 CALL timeset(routinen, handle)
3888
3889 NULLIFY (blacs_env_all, cell_to_index, cmo_coeff_source, cmatrix, dft_control, eigenvals, &
3890 fm_mo_coeff_source, &
3891 fm_struct, kp, kp_env, kpoints_scf, ks_env, mos_kp, &
3892 overlap_deriv, para_env, para_env_kp, rmatrix, sab_kp, sab_orb)
3893 IF (PRESENT(rcc)) THEN
3894 mark_used(rcc)
3895 END IF
3896
3897 CALL get_qs_env(qs_env, dft_control=dft_control, kpoints=kpoints_scf, ks_env=ks_env, &
3898 para_env=para_env, sab_orb=sab_orb)
3899 cpassert(ASSOCIATED(dft_control))
3900 cpassert(ASSOCIATED(kpoints_scf))
3901 cpassert(ASSOCIATED(ks_env))
3902 cpassert(ASSOCIATED(para_env))
3903 cpassert(ASSOCIATED(sab_orb))
3904
3905 CALL get_kpoint_info(kpoints_scf, nkp=nkp, kp_range=kp_range, kp_env=kp_env, &
3906 para_env_kp=para_env_kp, blacs_env_all=blacs_env_all, &
3907 cell_to_index=cell_to_index, sab_nl=sab_kp)
3908 IF (kp_range(2) >= kp_range(1)) THEN
3909 cpassert(ASSOCIATED(kp_env))
3910 END IF
3911 cpassert(ASSOCIATED(para_env_kp))
3912 cpassert(ASSOCIATED(blacs_env_all))
3913 cpassert(ASSOCIATED(cell_to_index))
3914 cpassert(ASSOCIATED(sab_kp))
3915
3916 CALL build_overlap_matrix(ks_env, matrixkp_s=overlap_deriv, nderivative=1, &
3917 basis_type_a="ORB", basis_type_b="ORB", sab_nl=sab_orb, &
3918 ext_kpoints=kpoints_scf)
3919
3920 nspin = dft_control%nspins
3921 CALL dbcsr_get_info(overlap_deriv(1, 1)%matrix, nfullrows_total=nao)
3922 ALLOCATE (nmo_spin(nspin), source=0)
3923 IF (kp_range(2) >= kp_range(1)) THEN
3924 kp => kp_env(1)%kpoint_env
3925 mos_kp => kp%mos
3926 cpassert(ASSOCIATED(mos_kp))
3927 DO ispin = 1, nspin
3928 CALL get_mo_set(mos_kp(ispin), nmo=nmo_spin(ispin))
3929 END DO
3930 END IF
3931 CALL para_env%max(nmo_spin)
3932 ALLOCATE (dipole(nspin, nkp, 3, maxval(nmo_spin), maxval(nmo_spin)), source=z_zero)
3933 IF (PRESENT(nmo_spin_out)) THEN
3934 ALLOCATE (nmo_spin_out(nspin))
3935 nmo_spin_out(:) = nmo_spin(:)
3936 END IF
3937
3938 ALLOCATE (rmatrix, cmatrix)
3939 CALL dbcsr_create(rmatrix, template=overlap_deriv(1, 1)%matrix, &
3940 matrix_type=dbcsr_type_antisymmetric)
3941 CALL dbcsr_create(cmatrix, template=overlap_deriv(1, 1)%matrix, &
3942 matrix_type=dbcsr_type_symmetric)
3943 CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_kp)
3944 CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_kp)
3945
3946 DO ikp = 1, nkp
3947 my_kpgrp = (ikp >= kp_range(1) .AND. ikp <= kp_range(2))
3948 IF (my_kpgrp) THEN
3949 ikp_local = ikp - kp_range(1) + 1
3950 kp => kp_env(ikp_local)%kpoint_env
3951 mos_kp => kp%mos
3952 ELSE
3953 NULLIFY (kp, mos_kp)
3954 END IF
3955 DO ispin = 1, nspin
3956 nmo = nmo_spin(ispin)
3957 ALLOCATE (eigenvalues_kp(nmo), source=0.0_dp)
3958
3959 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nmo, &
3960 para_env=para_env, context=blacs_env_all)
3961 CALL cp_cfm_create(cfm_mo_coeff_global, fm_struct)
3962 CALL cp_cfm_create(cfm_work, fm_struct)
3963 CALL cp_fm_struct_release(fm_struct)
3964 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo, ncol_global=nmo, &
3965 para_env=para_env, context=blacs_env_all)
3966 CALL cp_cfm_create(cfm_moment, fm_struct)
3967 CALL cp_fm_struct_release(fm_struct)
3968 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, &
3969 para_env=para_env, context=blacs_env_all)
3970 CALL cp_cfm_create(cfm_operator, fm_struct)
3971 CALL cp_fm_create(fm_operator_real, fm_struct)
3972 CALL cp_fm_create(fm_operator_imag, fm_struct)
3973 CALL cp_fm_struct_release(fm_struct)
3974
3975 IF (my_kpgrp) THEN
3976 CALL get_mo_set(mos_kp(ispin), eigenvalues=eigenvals)
3977 cpassert(ASSOCIATED(eigenvals))
3978 IF (kpoints_scf%use_real_wfn) THEN
3979 NULLIFY (cmo_coeff_source)
3980 fm_mo_coeff_source => mos_kp(ispin)%mo_coeff
3981 ELSE
3982 NULLIFY (fm_mo_coeff_source)
3983 cmo_coeff_source => mos_kp(ispin)%cmo_coeff
3984 END IF
3985 IF (para_env_kp%is_source()) eigenvalues_kp(1:nmo) = eigenvals(1:nmo)
3986 ELSE
3987 NULLIFY (cmo_coeff_source)
3988 NULLIFY (fm_mo_coeff_source)
3989 END IF
3990 IF (kpoints_scf%use_real_wfn) THEN
3991 CALL cp_fm_create(fm_mo_coeff_global, cfm_mo_coeff_global%matrix_struct)
3992 IF (ASSOCIATED(fm_mo_coeff_source)) THEN
3993 CALL cp_fm_copy_general(fm_mo_coeff_source, fm_mo_coeff_global, para_env)
3994 ELSE
3995 CALL cp_fm_copy_general(fm_dummy, fm_mo_coeff_global, para_env)
3996 END IF
3997 CALL cp_fm_to_cfm(msourcer=fm_mo_coeff_global, mtarget=cfm_mo_coeff_global)
3998 CALL cp_fm_release(fm_mo_coeff_global)
3999 ELSE
4000 IF (ASSOCIATED(cmo_coeff_source)) THEN
4001 CALL cp_cfm_copy_general(cmo_coeff_source, cfm_mo_coeff_global, para_env)
4002 ELSE
4003 CALL cp_cfm_copy_general(cfm_empty, cfm_mo_coeff_global, para_env)
4004 END IF
4005 END IF
4006 CALL para_env%sum(eigenvalues_kp)
4007
4008 DO i_dir = 1, 3
4009 CALL dbcsr_set(rmatrix, 0.0_dp)
4010 CALL dbcsr_set(cmatrix, 0.0_dp)
4011 CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=overlap_deriv, &
4012 ispin=i_dir + 1, xkp=kpoints_scf%xkp(:, ikp), &
4013 cell_to_index=cell_to_index, sab_nl=sab_kp)
4014
4015 ! Project A = rmatrix + i*cmatrix as C^H A C.
4016 CALL copy_dbcsr_to_fm(rmatrix, fm_operator_real)
4017 CALL copy_dbcsr_to_fm(cmatrix, fm_operator_imag)
4018 CALL cp_fm_to_cfm(fm_operator_real, fm_operator_imag, cfm_operator)
4019 CALL parallel_gemm("N", "N", nao, nmo, nao, z_one, cfm_operator, &
4020 cfm_mo_coeff_global, z_zero, cfm_work)
4021 CALL parallel_gemm("C", "N", nmo, nmo, nao, z_one, cfm_mo_coeff_global, &
4022 cfm_work, z_zero, cfm_moment)
4023
4024 DO n = 1, nmo
4025 DO m = 1, nmo
4026 IF (n == m) cycle
4027 energy_diff = eigenvalues_kp(m) - eigenvalues_kp(n)
4028 IF (abs(energy_diff) <= eps_degenerate) cycle
4029 CALL cp_cfm_get_element(cfm_moment, m, n, moment_element)
4030 IF (para_env%is_source()) THEN
4031 dipole(ispin, ikp, i_dir, n, m) = moment_element/energy_diff
4032 END IF
4033 END DO
4034 END DO
4035 END DO
4036 CALL cp_cfm_release(cfm_mo_coeff_global)
4037 CALL cp_cfm_release(cfm_work)
4038 CALL cp_cfm_release(cfm_moment)
4039 CALL cp_cfm_release(cfm_operator)
4040 CALL cp_fm_release(fm_operator_real)
4041 CALL cp_fm_release(fm_operator_imag)
4042 DEALLOCATE (eigenvalues_kp)
4043 END DO
4044 END DO
4045
4046 DO ispin = 1, nspin
4047 DO ikp = 1, nkp
4048 DO i_dir = 1, 3
4049 CALL para_env%sum(dipole(ispin, ikp, i_dir, :, :))
4050 END DO
4051 END DO
4052 END DO
4053
4054 CALL dbcsr_deallocate_matrix(cmatrix)
4055 CALL dbcsr_deallocate_matrix(rmatrix)
4056 CALL dbcsr_deallocate_matrix_set(overlap_deriv)
4057 DEALLOCATE (nmo_spin)
4058 CALL timestop(handle)
4059
4060 END SUBROUTINE qs_moment_kpoints_scf_mos
4061
4062! **************************************************************************************************
4063!> \brief Calculate and print dipole moment elements d_nm(k) = i<u_n|∇_k|u_m> for k-point
4064!> calculations
4065!> \param qs_env ...
4066!> \param nmoments ...
4067!> \param reference ...
4068!> \param ref_point ...
4069!> \param max_nmo ...
4070!> \param unit_number ...
4071!> \param do_soc include spin-orbit coupling (spinor sum-rule; requires an explicit KPOINT set)
4072!> \param unit_str unit of the printed dipoles (dipole unit of cp_units), default debye; the Berry
4073!> curvature is printed in the paired length unit squared
4074!> \param dip_file if present and non-empty, the dipole matrix elements are written to this file
4075!> instead of being printed to unit_number
4076!> \param bc_file if present and non-empty, the Berry curvature is written to this file instead of
4077!> being printed to unit_number
4078!> \author Shridhar Shanbhag
4079! **************************************************************************************************
4080 SUBROUTINE qs_moment_kpoints(qs_env, nmoments, reference, ref_point, max_nmo, unit_number, do_soc, &
4081 unit_str, dip_file, bc_file)
4082 TYPE(qs_environment_type), POINTER :: qs_env
4083 INTEGER, INTENT(IN) :: nmoments, reference, max_nmo
4084 REAL(dp), DIMENSION(:), INTENT(IN), POINTER :: ref_point
4085 INTEGER, INTENT(IN) :: unit_number
4086 LOGICAL, INTENT(IN) :: do_soc
4087 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: unit_str, dip_file, bc_file
4088 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_kpoints'
4089 CHARACTER(LEN=16) :: dip_label, len_label
4090 INTEGER :: bc_unit, dip_unit, iu
4091 INTEGER, DIMENSION(3) :: out_units
4092 LOGICAL :: use_sci
4093 LOGICAL :: to_bc_file, to_dip_file
4094 LOGICAL, DIMENSION(3) :: print_bc, print_dip
4095 REAL(kind=dp) :: dip_fac, len_fac
4096 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
4097 COMPLEX(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dipole_to_print
4098 COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
4099 ALLOCATABLE :: dipole
4100 INTEGER :: handle, i_dir, ikp, nmo_dim, nkp, nao, &
4101 num_pe, mepos, n, m, &
4102 ispin, nspin, nspin_dip, nmin, nmax, homo
4103 INTEGER, DIMENSION(:), ALLOCATABLE :: nmo_spin_scf
4104 LOGICAL :: explicit_kpnts, explicit_kpset, use_scf_mos
4105 REAL(kind=dp), DIMENSION(3) :: rcc
4106 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: xkp
4107 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: bc_to_print
4108 REAL(kind=dp), DIMENSION(:, :, :, :), ALLOCATABLE :: berry_c
4109 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
4110 TYPE(mp_para_env_type), POINTER :: para_env
4111 TYPE(section_vals_type), POINTER :: kpnts, kpset
4112 CHARACTER(LEN=default_string_length), &
4113 DIMENSION(:), POINTER :: special_pnts
4114
4115 CALL timeset(routinen, handle)
4116
4117 IF (nmoments > 1) cpabort("KPOINT quadrupole and higher moments not implemented.")
4118 IF (max_nmo < 0) cpabort("Negative maximum number of molecular orbitals max_nmo provided.")
4119
4120 CALL get_qs_env(qs_env, &
4121 para_env=para_env, &
4122 matrix_ks_kp=matrix_ks_kp, &
4123 mos=mos)
4124
4125 CALL get_mo_set(mo_set=mos(1), nao=nao)
4126 CALL get_xkp_for_dipole_calc(qs_env, xkp, special_pnts)
4127 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
4128 nspin = SIZE(matrix_ks_kp, 1)
4129 nkp = SIZE(xkp, 2)
4130
4131 kpset => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINT_SET")
4132 kpnts => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINTS")
4133 CALL section_vals_get(kpset, explicit=explicit_kpset)
4134 CALL section_vals_get(kpnts, explicit=explicit_kpnts)
4135 use_scf_mos = .NOT. explicit_kpset .AND. .NOT. explicit_kpnts
4136
4137 IF (do_soc .AND. use_scf_mos) THEN
4138 CALL cp_abort(__location__, "MOMENTS with SOC requires an explicit KPOINT_SET or "// &
4139 "KPOINTS section (the interband sum-rule / deep path).")
4140 END IF
4141
4142 ! output units: d_nm in the dipole unit, Berry curvature in the paired length unit squared
4143 IF (PRESENT(unit_str)) THEN
4144 CALL get_moment_print_units(unit_str, dip_fac, len_fac, dip_label, len_label, use_sci)
4145 ELSE
4146 CALL get_moment_print_units("debye", dip_fac, len_fac, dip_label, len_label, use_sci)
4147 END IF
4148 mark_used(use_sci)
4149
4150 IF (unit_number > 0) WRITE (unit_number, fmt="(/,T2,A)") &
4151 '!-----------------------------------------------------------------------------!'
4152 IF (unit_number > 0) WRITE (unit_number, "(T22,A)") "Periodic Dipole Matrix Elements ["// &
4153 trim(dip_label)//"]"
4154 IF (unit_number > 0 .AND. do_soc) WRITE (unit_number, "(T20,A)") &
4155 "with Spin-Orbit Coupling (spinor bands)"
4156
4157 ! optional separate files, written by the I/O rank only; data written to a file is not
4158 ! repeated in the main MOMENTS output
4159 to_dip_file = .false.
4160 to_bc_file = .false.
4161 IF (PRESENT(dip_file)) to_dip_file = (len_trim(dip_file) > 0)
4162 IF (PRESENT(bc_file)) to_bc_file = (len_trim(bc_file) > 0)
4163 IF (unit_number > 0 .AND. to_dip_file) WRITE (unit_number, "(/,T3,A)") &
4164 "Dipole matrix elements written to file: "//trim(dip_file)
4165 IF (unit_number > 0 .AND. to_bc_file) WRITE (unit_number, "(/,T3,A)") &
4166 "Berry curvature written to file: "//trim(bc_file)
4167 dip_unit = -1
4168 bc_unit = -1
4169 IF (para_env%is_source()) THEN
4170 IF (to_dip_file) THEN
4171 CALL open_file(file_name=trim(dip_file), file_status="REPLACE", file_action="WRITE", &
4172 unit_number=dip_unit)
4173 WRITE (dip_unit, fmt="(/,T2,A)") &
4174 '!-----------------------------------------------------------------------------!'
4175 WRITE (dip_unit, "(T22,A)") "Periodic Dipole Matrix Elements ["//trim(dip_label)//"]"
4176 IF (do_soc) WRITE (dip_unit, "(T20,A)") "with Spin-Orbit Coupling (spinor bands)"
4177 END IF
4178 IF (to_bc_file) THEN
4179 CALL open_file(file_name=trim(bc_file), file_status="REPLACE", file_action="WRITE", &
4180 unit_number=bc_unit)
4181 WRITE (bc_unit, fmt="(/,T2,A)") &
4182 '!-----------------------------------------------------------------------------!'
4183 WRITE (bc_unit, "(T22,A)") "Berry Curvature ["//trim(len_label)//"^2]"
4184 IF (do_soc) WRITE (bc_unit, "(T20,A)") "with Spin-Orbit Coupling (spinor bands)"
4185 END IF
4186 END IF
4187 ! output targets: main output, dipole file, Berry curvature file; each table goes
4188 ! to its file if one is given, otherwise to the main output
4189 out_units = [unit_number, dip_unit, bc_unit]
4190 print_dip = [.NOT. to_dip_file, to_dip_file, .false.]
4191 print_bc = [.NOT. to_bc_file, .false., to_bc_file]
4192 IF (to_dip_file .AND. to_bc_file) out_units(1) = -1
4193
4194 nspin_dip = nspin
4195 IF (do_soc) THEN
4196 ! spinor sum-rule: SOC merges the two spin channels into 2*nao spinor bands
4197 CALL qs_moment_kpoints_deep_soc(qs_env, &
4198 xkp, &
4199 dipole, &
4200 rcc, &
4201 berry_c, &
4202 do_parallel=.true.)
4203 nmo_dim = 2*nao
4204 nspin_dip = 1
4205 ELSE IF (use_scf_mos) THEN
4206 CALL qs_moment_kpoints_scf_mos(qs_env, dipole, rcc, nmo_spin_scf)
4207 nmo_dim = SIZE(dipole, 4)
4208 ALLOCATE (berry_c(nspin, nkp, 3, nmo_dim), source=0.0_dp)
4209 DO ispin = 1, nspin
4210 DO ikp = 1, nkp
4211 DO i_dir = 1, 3
4212 DO n = 1, nmo_dim
4213 DO m = 1, nmo_dim
4214 IF (n == m) cycle
4215 berry_c(ispin, ikp, i_dir, n) = berry_c(ispin, ikp, i_dir, n) &
4216 + 2*aimag(dipole(ispin, ikp, 1 + mod(i_dir, 3), n, m)* &
4217 dipole(ispin, ikp, 1 + mod(i_dir + 1, 3), m, n))
4218 END DO
4219 END DO
4220 END DO
4221 END DO
4222 END DO
4223 ELSE
4224 CALL qs_moment_kpoints_deep(qs_env, &
4225 xkp, &
4226 dipole, &
4227 rcc, &
4228 berry_c, &
4229 do_parallel=.true.)
4230 nmo_dim = nao
4231 END IF
4232
4233 ALLOCATE (dipole_to_print(3, nmo_dim, nmo_dim), source=z_zero)
4234 ALLOCATE (bc_to_print(3, nmo_dim), source=0.0_dp)
4235
4236 mepos = para_env%mepos
4237 num_pe = para_env%num_pe
4238
4239 DO ikp = 1, nkp
4240 DO ispin = 1, nspin_dip
4241 IF (do_soc) THEN
4242 ! spinor HOMO = number of occupied spin-up + spin-down states (mirrors SOC_ev)
4243 CALL get_mo_set(mo_set=mos(1), homo=homo)
4244 CALL get_mo_set(mo_set=mos(nspin), homo=n)
4245 homo = homo + n
4246 ELSE
4247 CALL get_mo_set(mo_set=mos(ispin), homo=homo)
4248 END IF
4249 nmin = max(1, homo - (max_nmo - 1)/2)
4250 nmax = min(nmo_dim, homo + max_nmo/2)
4251 IF (max_nmo == 0) THEN
4252 nmin = 1
4253 nmax = nmo_dim
4254 END IF
4255 IF (use_scf_mos) THEN
4256 nmax = min(nmax, nmo_spin_scf(ispin))
4257 END IF
4258 dipole_to_print = 0.0_dp
4259 bc_to_print = 0.0_dp
4260 IF (use_scf_mos) THEN
4261 dipole_to_print(:, :, :) = dipole(ispin, ikp, :, :, :)
4262 bc_to_print(:, :) = berry_c(ispin, ikp, :, :)
4263 ELSE IF (mod(ikp - 1, num_pe) == mepos) THEN
4264 dipole_to_print(:, :, :) = dipole(ispin, ceiling(real(ikp)/num_pe), :, :, :)
4265 bc_to_print(:, :) = berry_c(ispin, ceiling(real(ikp)/num_pe), :, :)
4266 END IF
4267 IF (.NOT. use_scf_mos) THEN
4268 CALL para_env%sum(dipole_to_print)
4269 CALL para_env%sum(bc_to_print)
4270 END IF
4271 DO iu = 1, SIZE(out_units)
4272 IF (out_units(iu) <= 0) cycle
4273 IF (special_pnts(ikp) /= "") WRITE (out_units(iu), "(/,2X,A,A)") &
4274 "Special point: ", adjustl(trim(special_pnts(ikp)))
4275 WRITE (out_units(iu), "(/,1X,A,I3,1X,3(A,1F12.6))") &
4276 "Kpoint:", ikp, ", kx:", xkp(1, ikp), ", ky:", xkp(2, ikp), ", kz:", xkp(3, ikp)
4277 IF (nspin > 1 .AND. .NOT. do_soc) WRITE (out_units(iu), "(/,2X,A,I2)") &
4278 "Open Shell System. Spin:", ispin
4279 IF (print_dip(iu)) THEN
4280 WRITE (out_units(iu), "(2X,A)") " kp n m Re(dx_nm) Im(dx_nm) &
4281 & Re(dy_nm) Im(dy_nm) Re(dz_nm) Im(dz_nm)"
4282 DO n = nmin, nmax
4283 DO m = nmin, nmax
4284 IF (n == m) cycle
4285 WRITE (out_units(iu), "(2X,I4,2I4,6(G11.3))") ikp, n, m, dipole_to_print(1:3, n, m)*dip_fac
4286 END DO
4287 END DO
4288 END IF
4289 IF (print_bc(iu)) THEN
4290 WRITE (out_units(iu), "(/,1X,A)") "Berry Curvature ["//trim(len_label)//"^2]"
4291 WRITE (out_units(iu), "(2X,A)") " kp n YZ ZX XY"
4292 DO n = nmin, nmax
4293 WRITE (out_units(iu), "(2X,2I5,3(1X,G11.3))") &
4294 ikp, n, bc_to_print(1:3, n)*len_fac**2
4295 END DO
4296 END IF
4297 END DO
4298 END DO
4299 END DO
4300 IF (dip_unit > 0) CALL close_file(unit_number=dip_unit)
4301 IF (bc_unit > 0) CALL close_file(unit_number=bc_unit)
4302 DEALLOCATE (dipole_to_print, bc_to_print, berry_c, dipole)
4303 IF (ALLOCATED(nmo_spin_scf)) DEALLOCATE (nmo_spin_scf)
4304 DEALLOCATE (special_pnts, xkp)
4305
4306 CALL timestop(handle)
4307
4308 END SUBROUTINE qs_moment_kpoints
4309
4310END MODULE qs_moments
ncoset
static GRID_HOST_DEVICE int ncoset(const int l)
Number of Cartesian orbitals up to given angular momentum quantum.
Definition grid_common.h:82
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
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:96
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:104
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_result_methods::put_results
Definition cp_result_methods.F:39
mathlib::gemm_square
Definition mathlib.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
ai_angmom
Calculation of the angular momentum integrals over Cartesian Gaussian-type functions.
Definition ai_angmom.F:17
ai_angmom::angmom
subroutine, public angmom(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, rac, rbc, angab)
Real angular-momentum kernel (r-C) cross grad, equal to i times L.
Definition ai_angmom.F:54
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
ai_moments::diff_momop
subroutine, public diff_momop(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, difmab, mab_ext, deltar, lambda, iatom, jatom)
This returns the derivative of the moment integrals [a|\mu|b]. By default, it differentiates the prim...
Definition ai_moments.F:390
ai_moments::contract_expikr
subroutine, public contract_expikr(cos_block, sin_block, iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, expab, work, ldwork, ordered)
Contract complex primitive integrals into real cosine and sine AO blocks.
Definition ai_moments.F:80
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:131
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::mattiat2019
integer, save, public mattiat2019
Definition bibliography.F:36
block_p_types
collect pointers to a block of reals
Definition block_p_types.F:14
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:234
cell_types::pbc_stable
real(kind=dp) function, dimension(3), public pbc_stable(r, cell)
Apply a stable periodic-image convention for k-point Bloch gauges.
Definition cell_types.F:458
commutator_rpnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition commutator_rpnl.F:14
commutator_rpnl::build_com_mom_nl
subroutine, public build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point)
Calculate [r,Vnl] (matrix_rv), r x [r,Vnl] (matrix_rxrv) or [rr,Vnl] (matrix_rrv) in AO basis....
Definition commutator_rpnl.F:76
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_det
subroutine, public cp_cfm_det(matrix_a, det_a)
Computes the determinant (with a correct sign even in parallel environment!) of a complex square matr...
Definition cp_cfm_basic_linalg.F:80
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_get_element
subroutine, public cp_cfm_get_element(matrix, irow_global, icol_global, alpha)
Get the matrix element by its global index.
Definition cp_cfm_types.F:256
cp_cfm_types::cp_cfm_copy_general
subroutine, public cp_cfm_copy_general(source, destination, para_env)
General copy of a cfm matrix to another cfm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition cp_cfm_types.F:905
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:845
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:632
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:793
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_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_get_matrix_type
character function, public dbcsr_get_matrix_type(matrix)
...
Definition cp_dbcsr_api.F:861
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
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_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
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_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_trace
subroutine, public dbcsr_trace(matrix, trace)
Computes the trace of the given matrix, also known as the sum of its diagonal elements.
Definition cp_dbcsr_contrib.F:458
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:372
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
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:940
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:344
cp_files
Utility routines to open and close files.
Definition cp_files.F:22
cp_files::open_file
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
Definition cp_files.F:197
cp_files::close_file
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
Definition cp_files.F:64
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_double
subroutine, public cp_fm_struct_double(fmstruct, struct, context, col, row)
creates a struct with twice the number of blocks on each core. If matrix A has to be multiplied with ...
Definition cp_fm_struct.F:527
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_copy_general
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition cp_fm_types.F:1587
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1046
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:536
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition cp_parser_methods.F:20
cp_parser_methods::read_float_object
elemental subroutine, public read_float_object(string, object, error_message)
Returns a floating point number read from a string including fraction like z1/z2.
Definition cp_parser_methods.F:1265
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition cp_result_methods.F:17
cp_result_methods::cp_results_erase
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
Definition cp_result_methods.F:365
cp_result_types
set of type/routines to handle the storage of results in force_envs
Definition cp_result_types.F:18
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_units_au
integer, parameter, public cp_units_au
Definition cp_units.F:66
cp_units::cp_unit_from_cp2k1
real(kind=dp) function, public cp_unit_from_cp2k1(value, unit, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1253
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1311
cp_units::cp_units_debye
integer, parameter, public cp_units_debye
Definition cp_units.F:124
cp_units::cp_units_c_m
integer, parameter, public cp_units_c_m
Definition cp_units.F:124
cp_units::cp_unit_create
subroutine, public cp_unit_create(unit, string)
creates a unit parsing a string
Definition cp_units.F:174
cp_units::cp_units_e_angstrom
integer, parameter, public cp_units_e_angstrom
Definition cp_units.F:124
cp_units::cp_ukind_dipole
integer, parameter, public cp_ukind_dipole
Definition cp_units.F:50
cp_units::cp_unit_release
elemental subroutine, public cp_unit_release(unit)
releases the given unit
Definition cp_units.F:587
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::use_mom_ref_coac
integer, parameter, public use_mom_ref_coac
Definition input_constants.F:520
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
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::max_line_length
integer, parameter, public max_line_length
Definition kinds.F:59
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_k_r_trafo_simple
Implements simple k-space/R-space transformations and Gamma-point MIC reconstruction.
Definition kpoint_k_r_trafo_simple.F:16
kpoint_k_r_trafo_simple::rs_to_kp
subroutine, public rs_to_kp(rs_real, ks_complex, index_to_cell, xkp, deriv_direction, hmat)
Integrate RS matrices (stored as Fortran array) into a kpoint matrix at given kp.
Definition kpoint_k_r_trafo_simple.F:241
kpoint_k_r_trafo_simple::replicate_rs_matrices
subroutine, public replicate_rs_matrices(rs_dbcsr_in, kpoint_in, rs_array_out, cell_to_index_out)
Convert dbcsr matrices representing operators in real-space image cells to arrays.
Definition kpoint_k_r_trafo_simple.F:70
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::rskp_transform
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
Definition kpoint_methods.F:1157
kpoint_methods::kpoint_init_cell_index
subroutine, public kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
Generates the mapping of cell indices and linear RS index CELL (0,0,0) is always mapped to index 1.
Definition kpoint_methods.F:1023
kpoint_methods::kpoint_initialize
subroutine, public kpoint_initialize(kpoint, particle_set, cell)
Generate the kpoints and initialize the kpoint environment.
Definition kpoint_methods.F:377
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::kpoint_release
subroutine, public kpoint_release(kpoint)
Release a kpoint environment, deallocate all data.
Definition kpoint_types.F:321
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:613
kpoint_types::kpoint_create
subroutine, public kpoint_create(kpoint)
Create a kpoint environment.
Definition kpoint_types.F:268
kpoint_types::read_kpoint_section
subroutine, public read_kpoint_section(kpoint, kpoint_section, a_vec, cell)
Read the kpoint input section.
Definition kpoint_types.F:820
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:143
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::geeig_right
subroutine, public geeig_right(a_in, b_in, eigenvalues, eigenvectors)
Solve the generalized eigenvalue equation for complex matrices A*v = B*v*λ
Definition mathlib.F:2203
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
moments_utils
Calculates the moment integrals <a|r^m|b>.
Definition moments_utils.F:15
moments_utils::get_reference_point
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Definition moments_utils.F:61
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::current_maxl
integer, save, public current_maxl
Definition orbital_pointers.F:40
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition orbital_pointers.F:43
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:123
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
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_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:483
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:976
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_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:440
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>.
Definition qs_moments.F:14
qs_moments::qs_moment_kpoints_deep
subroutine, public qs_moment_kpoints_deep(qs_env, xkp, dipole, rcc, berry_c, momentum, energy, do_parallel)
Calculate interband k-point dipole moments d_nm(k) = i<u_n|∇_k|u_m>, Berry curvature,...
Definition qs_moments.F:3272
qs_moments::build_local_moment_matrix
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type, all_images, minimum_image, neighbor_image, first_component)
...
Definition qs_moments.F:178
qs_moments::build_local_magmom_matrix
subroutine, public build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
...
Definition qs_moments.F:908
qs_moments::qs_moment_kpoints_deep_soc
subroutine, public qs_moment_kpoints_deep_soc(qs_env, xkp, dipole, rcc, berry_c, momentum, energy, do_parallel)
Calculate interband k-point dipole moments d_nm(k) = i<u_n|∇_k|u_m>, Berry curvature,...
Definition qs_moments.F:3573
qs_moments::print_moments_nl
subroutine, public print_moments_nl(unit_number, nmom, rlab, mmom, rmom_vel)
...
Definition qs_moments.F:2586
qs_moments::set_label
subroutine, public set_label(label, ix, iy, iz)
...
Definition qs_moments.F:2373
qs_moments::qs_moment_kpoints
subroutine, public qs_moment_kpoints(qs_env, nmoments, reference, ref_point, max_nmo, unit_number, do_soc, unit_str, dip_file, bc_file)
Calculate and print dipole moment elements d_nm(k) = i<u_n|∇_k|u_m> for k-point calculations.
Definition qs_moments.F:4082
qs_moments::build_berry_moment_matrix
subroutine, public build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
...
Definition qs_moments.F:1114
qs_moments::op_orbbas
subroutine, public op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
...
Definition qs_moments.F:1836
qs_moments::build_local_moment_matrix_rs_img
subroutine, public build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc)
Calculate local moment matrix for a periodic system for all image cells.
Definition qs_moments.F:3030
qs_moments::qs_moment_berry_phase
subroutine, public qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number, unit_str)
...
Definition qs_moments.F:1381
qs_moments::qs_moment_locop
subroutine, public qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl, unit_str)
...
Definition qs_moments.F:1960
qs_moments::calculate_commutator_nl_terms
subroutine, public calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr, ref_point)
Calculate the expectation value of operators related to non-local potential: [r, Vnl],...
Definition qs_moments.F:2654
qs_moments::qs_moment_kpoints_scf_mos
subroutine, public qs_moment_kpoints_scf_mos(qs_env, dipole, rcc, nmo_spin_out)
Calculates interband k-point dipoles in the existing SCF MO basis.
Definition qs_moments.F:3847
qs_moments::get_multipoles_ao
subroutine, public get_multipoles_ao(qs_env, n_moments, matrix_multipole, rpoint)
The AO multipoles M^k_µν = <φ_µ|(r - r_0)^k|φ_ν> of every order up to n_moments, on the block pattern...
Definition qs_moments.F:442
qs_moments::op_orbbas_rtp
subroutine, public op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
...
Definition qs_moments.F:1866
qs_moments::build_local_moments_der_matrix
subroutine, public build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, ref_point, moments, basis_type, minimum_image, ordered, lambda, deltar, neighbor_image)
Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b > Optionally sto...
Definition qs_moments.F:502
qs_moments::print_moments
subroutine, public print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmom, rmom_vel, unit_str)
...
Definition qs_moments.F:2405
qs_moments::build_berry_kpoint_matrix
subroutine, public build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type, ordered, ordered_kpoints)
...
Definition qs_moments.F:1152
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:550
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_operators_ao::build_exp_ikr_matrix
subroutine, public build_exp_ikr_matrix(qs_env, op_sm_set, kvec, sab_orb_external, basis_type, force_periodic, cell_external)
Build real and imaginary AO matrices for exp(i*k*r).
Definition qs_operators_ao.F:74
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
qs_overlap::build_overlap_matrix
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p, ext_kpoints)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:121
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_types
Types and set_get for real time propagation depending on runtype and diagonalization method different...
Definition rt_propagation_types.F:21
rt_propagation_types::get_rtp
subroutine, public get_rtp(rtp, exp_h_old, exp_h_new, h_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, s_inv, s_half, s_minus_half, b_mat, c_mat, propagator_matrix, mixing, mixing_factor, s_der, dt, nsteps, sinvh, sinvh_imag, sinvb, admm_mos)
...
Definition rt_propagation_types.F:416
soc_pseudopotential_methods
Definition soc_pseudopotential_methods.F:8
soc_pseudopotential_methods::v_soc_xyz_from_pseudopotential
subroutine, public v_soc_xyz_from_pseudopotential(qs_env, mat_v_soc_xyz)
V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = x,...
Definition soc_pseudopotential_methods.F:63
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
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
block_p_types::block_p_type
Definition block_p_types.F:23
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_control_types::dft_control_type
Definition cp_control_types.F:786
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_result_types::cp_result_type
contains arbitrary information which need to be stored
Definition cp_result_types.F:73
cp_units::cp_unit_type
stores a unit
Definition cp_units.F:143
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_env_p_type
Definition kpoint_types.F:119
kpoint_types::kpoint_env_type
Keeps information about a specific k-point.
Definition kpoint_types.F:94
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:196
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
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_mo_types::mo_set_type
Definition qs_mo_types.F:53
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
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
rt_propagation_types::rt_prop_type
Definition rt_propagation_types.F:76