(git:5c1df3d)
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_cossin, &
17 cossin, &
18 diff_momop, &
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 get_cell
30 USE commutator_rpnl, ONLY: build_com_mom_nl
31 USE cp_blacs_env, ONLY: cp_blacs_env_type
32 USE cp_cfm_basic_linalg, ONLY: cp_cfm_det
33 USE cp_cfm_types, ONLY: cp_cfm_create, &
34 cp_cfm_get_info, &
35 cp_cfm_release, &
36 cp_cfm_type
37 USE cp_control_types, ONLY: dft_control_type
38 USE cp_dbcsr_api, ONLY: &
39 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_type, &
40 dbcsr_get_block_p, dbcsr_get_info, dbcsr_get_matrix_type, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, &
41 dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
42 USE cp_dbcsr_contrib, ONLY: dbcsr_dot, &
43 dbcsr_trace
44 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
45 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply, &
46 dbcsr_allocate_matrix_set, &
47 dbcsr_deallocate_matrix_set
48 USE cp_fm_struct, ONLY: cp_fm_struct_create, &
49 cp_fm_struct_double, &
50 cp_fm_struct_release, &
51 cp_fm_struct_type
52 USE cp_fm_types, ONLY: cp_fm_copy_general, &
53 cp_fm_create, &
54 cp_fm_get_element, &
55 cp_fm_get_info, &
56 cp_fm_release, &
57 cp_fm_set_all, &
58 cp_fm_type
59 USE cp_result_methods, ONLY: cp_results_erase, &
60 put_results
61 USE cp_result_types, ONLY: cp_result_type
62 USE distribution_1d_types, ONLY: distribution_1d_type
63 USE kinds, ONLY: default_string_length, &
64 dp, &
65 max_line_length
66 USE kpoint_k_r_trafo_simple, ONLY: replicate_rs_matrices, &
67 rs_to_kp
68 USE kpoint_methods, ONLY: kpoint_init_cell_index, &
69 kpoint_initialize, &
70 rskp_transform
71 USE kpoint_types, ONLY: get_kpoint_info, &
72 kpoint_env_p_type, &
73 kpoint_env_type, &
74 kpoint_type, &
75 kpoint_create, &
76 kpoint_release, &
77 read_kpoint_section
78 USE mathconstants, ONLY: pi, &
79 twopi, &
80 gaussi, &
81 z_zero
82 USE message_passing, ONLY: mp_para_env_type
83 USE moments_utils, ONLY: get_reference_point
84 USE orbital_pointers, ONLY: current_maxl, &
85 indco, &
86 ncoset
87 USE parallel_gemm_api, ONLY: parallel_gemm
88 USE particle_methods, ONLY: get_particle_set
89 USE particle_types, ONLY: particle_type
90 USE physcon, ONLY: bohr, &
91 debye, &
92 angstrom
93 USE qs_environment_types, ONLY: get_qs_env, &
94 qs_environment_type
95 USE qs_kind_types, ONLY: get_qs_kind, &
96 get_qs_kind_set, &
97 qs_kind_type
98 USE qs_ks_types, ONLY: get_ks_env, &
99 qs_ks_env_type
100 USE qs_mo_types, ONLY: get_mo_set, &
101 mo_set_type
102 USE qs_neighbor_list_types, ONLY: get_iterator_info, &
103 neighbor_list_iterate, &
104 neighbor_list_iterator_create, &
105 neighbor_list_iterator_p_type, &
106 neighbor_list_iterator_release, &
107 neighbor_list_set_p_type
108 USE qs_operators_ao, ONLY: build_exp_ikr_matrix, &
109 build_lin_mom_matrix
110 USE qs_overlap, ONLY: build_overlap_matrix
111 USE qs_rho_types, ONLY: qs_rho_get, &
112 qs_rho_type
113 USE rt_propagation_types, ONLY: get_rtp, &
114 rt_prop_type
115 USE cp_parser_methods, ONLY: read_float_object
116 USE input_section_types, ONLY: section_vals_get, &
117 section_vals_get_subs_vals, &
118 section_vals_type, &
119 section_vals_val_get
120 USE mathlib, ONLY: geeig_right, &
121 gemm_square
122 USE string_utilities, ONLY: uppercase
123
124#include "./base/base_uses.f90"
125
126 IMPLICIT NONE
127
128 PRIVATE
129
130 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_moments'
131
132 ! Public subroutines
133 PUBLIC :: build_berry_moment_matrix, build_local_moment_matrix
134 PUBLIC :: build_berry_kpoint_matrix, build_local_magmom_matrix
135 PUBLIC :: qs_moment_berry_phase, qs_moment_locop
136 PUBLIC :: qs_moment_kpoints, build_local_moment_matrix_rs_img
137 PUBLIC :: qs_moment_kpoints_deep
138 PUBLIC :: qs_moment_kpoints_scf_mos
139 PUBLIC :: build_local_moments_der_matrix
140 PUBLIC :: calculate_commutator_nl_terms, op_orbbas, op_orbbas_rtp, &
141 print_moments, print_moments_nl, set_label
142
143CONTAINS
144
145! **************************************************************************************************
146!> \brief ...
147!> \param qs_env ...
148!> \param moments ...
149!> \param nmoments ...
150!> \param ref_point ...
151!> \param ref_points ...
152!> \param basis_type ...
153!> \param all_images accumulate all neighbor images belonging to one atom pair
154!> \param minimum_image apply the half-cell image filter
155!> \param neighbor_image use the neighbor-list image vector for one atom-pair representative
156!> \param first_component first Cartesian moment component to write
157! **************************************************************************************************
158 SUBROUTINE build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type, &
159 all_images, minimum_image, neighbor_image, first_component)
160
161 TYPE(qs_environment_type), POINTER :: qs_env
162 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moments
163 INTEGER, INTENT(IN) :: nmoments
164 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
165 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
166 OPTIONAL :: ref_points
167 CHARACTER(len=*), OPTIONAL :: basis_type
168 LOGICAL, OPTIONAL :: all_images, minimum_image, neighbor_image
169 INTEGER, OPTIONAL :: first_component
170
171 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moment_matrix'
172
173 INTEGER :: component_start, handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, &
174 last_jatom, maxco, maxsgf, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
175 LOGICAL :: found, my_all_images, my_minimum_image, &
176 my_neighbor_image, new_atom_pair
177 REAL(kind=dp) :: dab
178 REAL(kind=dp), DIMENSION(3) :: half_cell
179 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
180 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
181 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
182 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
183 TYPE(cell_type), POINTER :: cell
184 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
185 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
186 TYPE(neighbor_list_iterator_p_type), &
187 DIMENSION(:), POINTER :: nl_iterator
188 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
189 POINTER :: sab_orb
190 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
191 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
192 TYPE(qs_kind_type), POINTER :: qs_kind
193
194 IF (nmoments < 1) RETURN
195
196 CALL timeset(routinen, handle)
197
198 my_all_images = .false.
199 IF (PRESENT(all_images)) my_all_images = all_images
200 my_minimum_image = .false.
201 IF (PRESENT(minimum_image)) my_minimum_image = minimum_image
202 my_neighbor_image = .false.
203 IF (PRESENT(neighbor_image)) my_neighbor_image = neighbor_image
204 cpassert(.NOT. (my_all_images .AND. my_neighbor_image))
205
206 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
207 cpassert(SIZE(moments) >= nm)
208 component_start = 1
209 IF (PRESENT(first_component)) component_start = first_component
210 cpassert(component_start >= 1 .AND. component_start <= nm)
211
212 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
213 CALL get_qs_env(qs_env=qs_env, &
214 qs_kind_set=qs_kind_set, &
215 particle_set=particle_set, cell=cell, &
216 sab_orb=sab_orb)
217
218 IF (my_minimum_image) THEN
219 half_cell = 0.5_dp*norm2(cell%hmat, dim=1)
220 END IF
221
222 nkind = SIZE(qs_kind_set)
223 ! Allocate work storage
224 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
225 maxco=maxco, maxsgf=maxsgf, &
226 basis_type=basis_type)
227
228 ALLOCATE (mab(maxco, maxco, nm))
229 mab(:, :, :) = 0.0_dp
230
231 ALLOCATE (work(maxco, maxsgf))
232 work(:, :) = 0.0_dp
233
234 ALLOCATE (mint(nm))
235 DO i = component_start, nm
236 NULLIFY (mint(i)%block)
237 END DO
238
239 ALLOCATE (basis_set_list(nkind))
240 DO ikind = 1, nkind
241 qs_kind => qs_kind_set(ikind)
242 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
243 IF (ASSOCIATED(basis_set_a)) THEN
244 basis_set_list(ikind)%gto_basis_set => basis_set_a
245 ELSE
246 NULLIFY (basis_set_list(ikind)%gto_basis_set)
247 END IF
248 END DO
249 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
250 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
251 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
252 iatom=iatom, jatom=jatom, r=rab)
253 basis_set_a => basis_set_list(ikind)%gto_basis_set
254 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
255 basis_set_b => basis_set_list(jkind)%gto_basis_set
256 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
257 associate( &
258 ! basis ikind
259 first_sgfa => basis_set_a%first_sgf, &
260 la_max => basis_set_a%lmax, &
261 la_min => basis_set_a%lmin, &
262 npgfa => basis_set_a%npgf, &
263 nsgfa => basis_set_a%nsgf_set, &
264 rpgfa => basis_set_a%pgf_radius, &
265 set_radius_a => basis_set_a%set_radius, &
266 sphi_a => basis_set_a%sphi, &
267 zeta => basis_set_a%zet, &
268 ! basis jkind, &
269 first_sgfb => basis_set_b%first_sgf, &
270 lb_max => basis_set_b%lmax, &
271 lb_min => basis_set_b%lmin, &
272 npgfb => basis_set_b%npgf, &
273 nsgfb => basis_set_b%nsgf_set, &
274 rpgfb => basis_set_b%pgf_radius, &
275 set_radius_b => basis_set_b%set_radius, &
276 sphi_b => basis_set_b%sphi, &
277 zetb => basis_set_b%zet)
278
279 nseta = basis_set_a%nset
280 nsetb = basis_set_b%nset
281
282 IF (inode == 1) last_jatom = 0
283
284 IF (my_minimum_image) THEN
285 IF (any(abs(rab(:)) > half_cell(:))) cycle
286 END IF
287
288 new_atom_pair = jatom /= last_jatom
289 IF (.NOT. my_all_images .AND. .NOT. new_atom_pair) THEN
290 cycle
291 END IF
292
293 IF (new_atom_pair) THEN
294 last_jatom = jatom
295
296 IF (iatom <= jatom) THEN
297 irow = iatom
298 icol = jatom
299 ELSE
300 irow = jatom
301 icol = iatom
302 END IF
303
304 DO i = component_start, nm
305 NULLIFY (mint(i)%block)
306 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
307 row=irow, col=icol, block=mint(i)%block, found=found)
308 cpassert(found .AND. ASSOCIATED(mint(i)%block))
309 mint(i)%block = 0._dp
310 END DO
311 END IF
312
313 ! Reference point and image geometry.
314 IF (PRESENT(ref_points)) THEN
315 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
316 ELSE IF (PRESENT(ref_point)) THEN
317 rc(:) = ref_point(:)
318 ELSE
319 rc(:) = 0._dp
320 END IF
321
322 IF (my_all_images) THEN
323 ! Keep the image represented by the current neighbor-list node.
324 ra(:) = pbc(particle_set(iatom)%r(:), cell)
325 rb(:) = ra(:) + rab(:)
326 rac(:) = pbc(rc, ra, cell)
327 rbc(:) = pbc(rc, rb, cell)
328 dab = norm2(rab)
329 ELSE IF (my_neighbor_image) THEN
330 ! Keep the neighbor-list image while selecting one atom-pair representative.
331 ra(:) = particle_set(iatom)%r(:)
332 rb(:) = ra(:) + rab(:)
333 rac(:) = pbc(rc, ra, cell)
334 rbc(:) = rac(:) + rab(:)
335 dab = norm2(rab)
336 ELSE
337 ! Fold atomic positions around the reference point.
338 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
339 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
340 rab(:) = ra(:) - rb(:)
341 rac(:) = ra(:) - rc(:)
342 rbc(:) = rb(:) - rc(:)
343 dab = norm2(rab)
344 END IF
345
346 DO iset = 1, nseta
347
348 ncoa = npgfa(iset)*ncoset(la_max(iset))
349 sgfa = first_sgfa(1, iset)
350
351 DO jset = 1, nsetb
352
353 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
354
355 ncob = npgfb(jset)*ncoset(lb_max(jset))
356 sgfb = first_sgfb(1, jset)
357
358 ! Calculate the primitive integrals
359 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
360 rpgfa(:, iset), la_min(iset), &
361 lb_max(jset), npgfb(jset), zetb(:, jset), &
362 rpgfb(:, jset), nmoments, rac, rbc, mab)
363
364 ! Contraction step
365 DO i = component_start, nm
366
367 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
368 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
369 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
370 0.0_dp, work(1, 1), SIZE(work, 1))
371
372 IF (iatom <= jatom) THEN
373
374 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
375 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
376 work(1, 1), SIZE(work, 1), &
377 1.0_dp, mint(i)%block(sgfa, sgfb), &
378 SIZE(mint(i)%block, 1))
379
380 ELSE
381
382 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
383 1.0_dp, work(1, 1), SIZE(work, 1), &
384 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
385 1.0_dp, mint(i)%block(sgfb, sgfa), &
386 SIZE(mint(i)%block, 1))
387
388 END IF
389
390 END DO
391
392 END DO
393 END DO
394 END associate
395
396 END DO
397 CALL neighbor_list_iterator_release(nl_iterator)
398
399 ! Release work storage
400 DEALLOCATE (mab, basis_set_list)
401 DEALLOCATE (work)
402 DO i = component_start, nm
403 NULLIFY (mint(i)%block)
404 END DO
405 DEALLOCATE (mint)
406
407 CALL timestop(handle)
408
409 END SUBROUTINE build_local_moment_matrix
410
411! **************************************************************************************************
412!> \brief Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b >
413!> Optionally stores the multipole moments themselves for free.
414!> Note that the multipole moments are symmetric while their derivatives are anti-symmetric
415!> Only first derivatives are performed, e. g. x d/dy
416!> \param qs_env ...
417!> \param moments_der will contain the derivatives of the multipole moments
418!> \param nmoments_der order of the moments with derivatives
419!> \param nmoments order of the multipole moments (no derivatives, same output as
420!> build_local_moment_matrix, needs moments as arguments to store results)
421!> \param ref_point ...
422!> \param moments contains the multipole moments, optionally for free, up to order nmoments
423!> \param basis_type basis set used for the integral evaluation
424!> \param minimum_image apply the half-cell image filter
425!> \param ordered use ordered neighbor blocks and, for ordinary derivatives, the full neighbor list
426!> \param lambda optional atom selector for the right-hand derivative
427!> \param deltaR optional weights for left and right nuclear-coordinate derivatives
428!> \param neighbor_image use the neighbor-list image for coordinate construction
429! **************************************************************************************************
430 SUBROUTINE build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, &
431 ref_point, moments, basis_type, minimum_image, ordered, lambda, &
432 deltaR, neighbor_image)
433 TYPE(qs_environment_type), POINTER :: qs_env
434 TYPE(dbcsr_p_type), DIMENSION(:, :), &
435 INTENT(INOUT), POINTER :: moments_der
436 INTEGER, INTENT(IN) :: nmoments_der, nmoments
437 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
438 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
439 OPTIONAL, POINTER :: moments
440 CHARACTER(len=*), OPTIONAL :: basis_type
441 LOGICAL, OPTIONAL :: minimum_image, ordered
442 INTEGER, OPTIONAL :: lambda
443 REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: deltar
444 LOGICAL, OPTIONAL :: neighbor_image
445
446 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moments_der_matrix'
447
448 INTEGER :: dimders, handle, i, iatom, icol, ider, ii, ikind, inode, ipgf, irow, iset, j, &
449 jatom, jkind, jpgf, jset, last_jatom, lda_min, m_dim, maxco, maxsgf, na, nb, ncoa, ncob, nda, ndb, &
450 nders, nkind, nm, nmom_build, nseta, nsetb, sgfa, sgfb
451 LOGICAL :: accumulate_images, found, my_minimum_image, my_neighbor_image, &
452 my_ordered, new_atom_pair
453 REAL(kind=dp) :: dab, derivative_sign
454 REAL(kind=dp), DIMENSION(3) :: half_cell
455 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
456 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
457 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
458 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
459 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab_tmp
460 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mom_block
461 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mom_block_der
462 TYPE(cell_type), POINTER :: cell
463 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
464 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
465 TYPE(neighbor_list_iterator_p_type), &
466 DIMENSION(:), POINTER :: nl_iterator
467 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
468 POINTER :: sab_all, sab_orb
469 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
470 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
471 TYPE(qs_kind_type), POINTER :: qs_kind
472
473 nmom_build = max(nmoments, nmoments_der) ! build moments up to order nmom_buiod
474 IF (nmom_build < 1) RETURN
475
476 CALL timeset(routinen, handle)
477
478 nders = 1 ! only first order derivatives
479 dimders = ncoset(nders) - 1
480
481 my_minimum_image = .false.
482 IF (PRESENT(minimum_image)) my_minimum_image = minimum_image
483 my_ordered = .false.
484 IF (PRESENT(ordered)) my_ordered = ordered
485 my_neighbor_image = my_ordered
486 IF (PRESENT(neighbor_image)) my_neighbor_image = neighbor_image
487 cpassert(.NOT. (PRESENT(deltar) .AND. PRESENT(lambda)))
488 IF (PRESENT(lambda) .OR. PRESENT(deltar)) THEN
489 cpassert(my_ordered)
490 END IF
491 accumulate_images = my_ordered .AND. my_neighbor_image
492 IF (PRESENT(deltar) .OR. PRESENT(lambda)) accumulate_images = .false.
493 derivative_sign = 1._dp
494 ! diff_momop returns the right nuclear-center derivative. The former
495 ! rRc_xyz_der_ao interface exposed the corresponding electronic derivative
496 ! on the ordinary ordered path.
497 IF (my_ordered .AND. .NOT. PRESENT(deltar) .AND. .NOT. PRESENT(lambda)) THEN
498 derivative_sign = -1._dp
499 END IF
500
501 NULLIFY (qs_kind_set, particle_set, sab_all, sab_orb, cell)
502 CALL get_qs_env(qs_env=qs_env, &
503 qs_kind_set=qs_kind_set, &
504 particle_set=particle_set, &
505 cell=cell, &
506 sab_orb=sab_orb, &
507 sab_all=sab_all)
508
509 IF (my_ordered) THEN
510 cpassert(nmoments_der > 0)
511 cpassert(dbcsr_get_matrix_type(moments_der(1, 1)%matrix) == dbcsr_type_no_symmetry)
512 IF (nmoments > 0) THEN
513 cpassert(PRESENT(moments))
514 cpassert(dbcsr_get_matrix_type(moments(1)%matrix) == dbcsr_type_no_symmetry)
515 END IF
516 END IF
517
518 IF (my_minimum_image) THEN
519 half_cell = 0.5_dp*norm2(cell%hmat, dim=1)
520 END IF
521
522 nkind = SIZE(qs_kind_set)
523
524 ! Work storage
525 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
526 maxco=maxco, maxsgf=maxsgf, &
527 basis_type=basis_type)
528
529 IF (nmoments > 0) THEN
530 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
531 cpassert(SIZE(moments) == nm)
532 ! storage for integrals
533 ALLOCATE (mab(maxco, maxco, nm))
534 ! blocks
535 mab(:, :, :) = 0.0_dp
536 ALLOCATE (mom_block(nm))
537 DO i = 1, nm
538 NULLIFY (mom_block(i)%block)
539 END DO
540 END IF
541
542 IF (nmoments_der > 0) THEN
543 m_dim = ncoset(nmoments_der) - 1
544 cpassert(SIZE(moments_der, dim=1) == m_dim)
545 cpassert(SIZE(moments_der, dim=2) == dimders)
546 ! storage for integrals
547 ALLOCATE (difmab(maxco, maxco, m_dim, dimders))
548 difmab(:, :, :, :) = 0.0_dp
549 ! blocks
550 ALLOCATE (mom_block_der(m_dim, dimders))
551 DO i = 1, m_dim
552 DO ider = 1, dimders
553 NULLIFY (mom_block_der(i, ider)%block)
554 END DO
555 END DO
556 END IF
557
558 ALLOCATE (work(maxco, maxsgf))
559 work(:, :) = 0.0_dp
560
561 NULLIFY (basis_set_a, basis_set_b, basis_set_list)
562 NULLIFY (qs_kind)
563 ALLOCATE (basis_set_list(nkind))
564 DO ikind = 1, nkind
565 qs_kind => qs_kind_set(ikind)
566 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
567 IF (ASSOCIATED(basis_set_a)) THEN
568 basis_set_list(ikind)%gto_basis_set => basis_set_a
569 ELSE
570 NULLIFY (basis_set_list(ikind)%gto_basis_set)
571 END IF
572 END DO
573
574 ! Calculate derivatives looping over neighbour list
575 NULLIFY (nl_iterator)
576 IF (my_ordered) THEN
577 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
578 ELSE
579 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
580 END IF
581 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
582 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
583 iatom=iatom, jatom=jatom, r=rab)
584 basis_set_a => basis_set_list(ikind)%gto_basis_set
585 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
586 basis_set_b => basis_set_list(jkind)%gto_basis_set
587 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
588 associate( &
589 ! basis ikind
590 first_sgfa => basis_set_a%first_sgf, &
591 la_max => basis_set_a%lmax, &
592 la_min => basis_set_a%lmin, &
593 npgfa => basis_set_a%npgf, &
594 nsgfa => basis_set_a%nsgf_set, &
595 rpgfa => basis_set_a%pgf_radius, &
596 set_radius_a => basis_set_a%set_radius, &
597 sphi_a => basis_set_a%sphi, &
598 zeta => basis_set_a%zet, &
599 ! basis jkind, &
600 first_sgfb => basis_set_b%first_sgf, &
601 lb_max => basis_set_b%lmax, &
602 lb_min => basis_set_b%lmin, &
603 npgfb => basis_set_b%npgf, &
604 nsgfb => basis_set_b%nsgf_set, &
605 rpgfb => basis_set_b%pgf_radius, &
606 set_radius_b => basis_set_b%set_radius, &
607 sphi_b => basis_set_b%sphi, &
608 zetb => basis_set_b%zet)
609
610 IF (SIZE(sphi_a, 1) == 0) cycle
611 IF (SIZE(sphi_b, 1) == 0) cycle
612
613 nseta = basis_set_a%nset
614 nsetb = basis_set_b%nset
615
616 IF (inode == 1) last_jatom = 0
617
618 IF (my_minimum_image) THEN
619 IF (any(abs(rab(:)) > half_cell(:))) cycle
620 END IF
621
622 new_atom_pair = jatom /= last_jatom
623 ! The ordinary ordered neighbor-image path accumulates all images for one atom pair.
624 ! Specialized derivative paths select one representative image, as in the former builders.
625 IF (.NOT. accumulate_images .AND. .NOT. new_atom_pair) cycle
626
627 IF (new_atom_pair) THEN
628 last_jatom = jatom
629
630 ! get blocks
631 IF (my_ordered) THEN
632 irow = iatom
633 icol = jatom
634 ELSE IF (iatom <= jatom) THEN
635 irow = iatom
636 icol = jatom
637 ELSE
638 irow = jatom
639 icol = iatom
640 END IF
641
642 IF (nmoments > 0) THEN
643 DO i = 1, nm
644 NULLIFY (mom_block(i)%block)
645 ! get block from pre calculated overlap matrix
646 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
647 row=irow, col=icol, block=mom_block(i)%block, found=found)
648 cpassert(found .AND. ASSOCIATED(mom_block(i)%block))
649 mom_block(i)%block = 0._dp
650 END DO
651 END IF
652 IF (nmoments_der > 0) THEN
653 DO i = 1, m_dim
654 DO ider = 1, dimders
655 NULLIFY (mom_block_der(i, ider)%block)
656 CALL dbcsr_get_block_p(matrix=moments_der(i, ider)%matrix, &
657 row=irow, col=icol, &
658 block=mom_block_der(i, ider)%block, &
659 found=found)
660 cpassert(found .AND. ASSOCIATED(mom_block_der(i, ider)%block))
661 mom_block_der(i, ider)%block = 0._dp
662 END DO
663 END DO
664 END IF
665 END IF
666
667 ! reference point
668 IF (PRESENT(ref_point)) THEN
669 rc(:) = ref_point(:)
670 ELSE
671 rc(:) = 0._dp
672 END IF
673 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
674 ! by folding around the center, such screwing can be avoided for a proper choice of center.
675 IF (my_neighbor_image) THEN
676 ra(:) = pbc(particle_set(iatom)%r(:), cell)
677 rb(:) = ra(:) + rab(:)
678 rac(:) = pbc(rc, ra, cell)
679 rbc(:) = rac(:) + rab(:)
680 ELSE IF (my_ordered) THEN
681 ra(:) = particle_set(iatom)%r(:)
682 rb(:) = particle_set(jatom)%r(:)
683 rab(:) = pbc(rb, ra, cell)
684 rac(:) = pbc(ra - rc, cell)
685 rbc(:) = pbc(rb - rc, cell)
686 ELSE
687 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
688 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
689 ! we dont use PBC at this point
690 rab(:) = ra(:) - rb(:)
691 rac(:) = ra(:) - rc(:)
692 rbc(:) = rb(:) - rc(:)
693 END IF
694 dab = norm2(rab)
695
696 DO iset = 1, nseta
697
698 ncoa = npgfa(iset)*ncoset(la_max(iset))
699 sgfa = first_sgfa(1, iset)
700
701 DO jset = 1, nsetb
702
703 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
704
705 ncob = npgfb(jset)*ncoset(lb_max(jset))
706 sgfb = first_sgfb(1, jset)
707
708 NULLIFY (mab_tmp)
709 ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
710 npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(nmom_build) - 1))
711
712 ! Calculate the primitive integrals (need l+1 for derivatives)
713 lda_min = max(0, la_min(iset) - 1)
714 CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), &
715 rpgfa(:, iset), lda_min, &
716 lb_max(jset) + 1, npgfb(jset), zetb(:, jset), &
717 rpgfb(:, jset), nmom_build, rac, rbc, mab_tmp)
718
719 IF (nmoments_der > 0) THEN
720 CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), &
721 rpgfa(:, iset), la_min(iset), &
722 lb_max(jset), npgfb(jset), zetb(:, jset), &
723 rpgfb(:, jset), lb_min(jset), &
724 nmoments_der, rac, rbc, difmab, mab_ext=mab_tmp, &
725 lambda=lambda, deltar=deltar, iatom=iatom, jatom=jatom)
726 END IF
727
728 IF (nmoments > 0) THEN
729 ! copy subset of mab_tmp (l+1) to mab (l)
730 mab = 0.0_dp
731 DO ii = 1, nm
732 na = 0
733 nda = 0
734 DO ipgf = 1, npgfa(iset)
735 nb = 0
736 ndb = 0
737 DO jpgf = 1, npgfb(jset)
738 DO j = 1, ncoset(lb_max(jset))
739 DO i = 1, ncoset(la_max(iset))
740 mab(i + na, j + nb, ii) = mab_tmp(i + nda, j + ndb, ii)
741 END DO ! i
742 END DO ! j
743 nb = nb + ncoset(lb_max(jset))
744 ndb = ndb + ncoset(lb_max(jset) + 1)
745 END DO ! jpgf
746 na = na + ncoset(la_max(iset))
747 nda = nda + ncoset(la_max(iset) + 1)
748 END DO ! ipgf
749 END DO
750 ! Contraction step
751 DO i = 1, nm
752
753 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
754 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
755 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
756 0.0_dp, work(1, 1), SIZE(work, 1))
757
758 IF (my_ordered .OR. iatom <= jatom) THEN
759 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
760 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
761 work(1, 1), SIZE(work, 1), &
762 1.0_dp, mom_block(i)%block(sgfa, sgfb), &
763 SIZE(mom_block(i)%block, 1))
764 ELSE
765 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
766 1.0_dp, work(1, 1), SIZE(work, 1), &
767 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
768 1.0_dp, mom_block(i)%block(sgfb, sgfa), &
769 SIZE(mom_block(i)%block, 1))
770 END IF
771 END DO
772 END IF
773
774 IF (nmoments_der > 0) THEN
775 DO i = 1, m_dim
776 DO ider = 1, dimders
777 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
778 derivative_sign, difmab(1, 1, i, ider), SIZE(difmab, 1), &
779 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
780 0._dp, work(1, 1), SIZE(work, 1))
781
782 IF (my_ordered .OR. iatom <= jatom) THEN
783 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
784 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
785 work(1, 1), SIZE(work, 1), &
786 1.0_dp, mom_block_der(i, ider)%block(sgfa, sgfb), &
787 SIZE(mom_block_der(i, ider)%block, 1))
788 ELSE
789 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
790 -1.0_dp, work(1, 1), SIZE(work, 1), &
791 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
792 1.0_dp, mom_block_der(i, ider)%block(sgfb, sgfa), &
793 SIZE(mom_block_der(i, ider)%block, 1))
794 END IF
795 END DO
796 END DO
797 END IF
798 DEALLOCATE (mab_tmp)
799 END DO
800 END DO
801 END associate
802 END DO
803 CALL neighbor_list_iterator_release(nl_iterator)
804
805 ! deallocations
806 DEALLOCATE (basis_set_list)
807 DEALLOCATE (work)
808 IF (nmoments > 0) THEN
809 DEALLOCATE (mab)
810 DO i = 1, nm
811 NULLIFY (mom_block(i)%block)
812 END DO
813 DEALLOCATE (mom_block)
814 END IF
815 IF (nmoments_der > 0) THEN
816 DEALLOCATE (difmab)
817 DO i = 1, m_dim
818 DO ider = 1, dimders
819 NULLIFY (mom_block_der(i, ider)%block)
820 END DO
821 END DO
822 DEALLOCATE (mom_block_der)
823 END IF
824
825 CALL timestop(handle)
826
827 END SUBROUTINE build_local_moments_der_matrix
828
829! **************************************************************************************************
830!> \brief ...
831!> \param qs_env ...
832!> \param magmom ...
833!> \param nmoments ...
834!> \param ref_point ...
835!> \param ref_points ...
836!> \param basis_type ...
837! **************************************************************************************************
838 SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
839
840 TYPE(qs_environment_type), POINTER :: qs_env
841 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom
842 INTEGER, INTENT(IN) :: nmoments
843 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
844 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
845 OPTIONAL :: ref_points
846 CHARACTER(len=*), OPTIONAL :: basis_type
847
848 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_magmom_matrix'
849
850 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, maxco, &
851 maxsgf, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
852 LOGICAL :: found
853 REAL(kind=dp) :: dab
854 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
855 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
856 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
857 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
858 TYPE(cell_type), POINTER :: cell
859 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
860 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
861 TYPE(neighbor_list_iterator_p_type), &
862 DIMENSION(:), POINTER :: nl_iterator
863 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
864 POINTER :: sab_orb
865 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
866 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
867 TYPE(qs_kind_type), POINTER :: qs_kind
868
869 IF (nmoments < 1) RETURN
870
871 CALL timeset(routinen, handle)
872
873 ! magnetic dipoles/angular moments only
874 nm = 3
875
876 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
877 CALL get_qs_env(qs_env=qs_env, &
878 qs_kind_set=qs_kind_set, &
879 particle_set=particle_set, cell=cell, &
880 sab_orb=sab_orb)
881
882 nkind = SIZE(qs_kind_set)
883
884 ! Allocate work storage
885 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
886 maxco=maxco, maxsgf=maxsgf)
887
888 ALLOCATE (mab(maxco, maxco, nm))
889 mab(:, :, :) = 0.0_dp
890
891 ALLOCATE (work(maxco, maxsgf))
892 work(:, :) = 0.0_dp
893
894 ALLOCATE (mint(nm))
895 DO i = 1, nm
896 NULLIFY (mint(i)%block)
897 END DO
898
899 ALLOCATE (basis_set_list(nkind))
900 DO ikind = 1, nkind
901 qs_kind => qs_kind_set(ikind)
902 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
903 IF (ASSOCIATED(basis_set_a)) THEN
904 basis_set_list(ikind)%gto_basis_set => basis_set_a
905 ELSE
906 NULLIFY (basis_set_list(ikind)%gto_basis_set)
907 END IF
908 END DO
909 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
910 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
911 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
912 iatom=iatom, jatom=jatom, r=rab)
913 basis_set_a => basis_set_list(ikind)%gto_basis_set
914 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
915 basis_set_b => basis_set_list(jkind)%gto_basis_set
916 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
917 associate( &
918 ! basis ikind
919 first_sgfa => basis_set_a%first_sgf, &
920 la_max => basis_set_a%lmax, &
921 la_min => basis_set_a%lmin, &
922 npgfa => basis_set_a%npgf, &
923 nsgfa => basis_set_a%nsgf_set, &
924 rpgfa => basis_set_a%pgf_radius, &
925 set_radius_a => basis_set_a%set_radius, &
926 sphi_a => basis_set_a%sphi, &
927 zeta => basis_set_a%zet, &
928 ! basis jkind, &
929 first_sgfb => basis_set_b%first_sgf, &
930 lb_max => basis_set_b%lmax, &
931 lb_min => basis_set_b%lmin, &
932 npgfb => basis_set_b%npgf, &
933 nsgfb => basis_set_b%nsgf_set, &
934 rpgfb => basis_set_b%pgf_radius, &
935 set_radius_b => basis_set_b%set_radius, &
936 sphi_b => basis_set_b%sphi, &
937 zetb => basis_set_b%zet)
938
939 nseta = basis_set_a%nset
940 nsetb = basis_set_b%nset
941
942 IF (iatom <= jatom) THEN
943 irow = iatom
944 icol = jatom
945 ELSE
946 irow = jatom
947 icol = iatom
948 END IF
949
950 DO i = 1, nm
951 NULLIFY (mint(i)%block)
952 CALL dbcsr_get_block_p(matrix=magmom(i)%matrix, &
953 row=irow, col=icol, block=mint(i)%block, found=found)
954 cpassert(found .AND. ASSOCIATED(mint(i)%block))
955 mint(i)%block = 0._dp
956 END DO
957
958 ! fold atomic position back into unit cell
959 IF (PRESENT(ref_points)) THEN
960 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
961 ELSE IF (PRESENT(ref_point)) THEN
962 rc(:) = ref_point(:)
963 ELSE
964 rc(:) = 0._dp
965 END IF
966 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
967 ! by folding around the center, such screwing can be avoided for a proper choice of center.
968 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
969 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
970 ! we dont use PBC at this point
971 rab(:) = ra(:) - rb(:)
972 rac(:) = ra(:) - rc(:)
973 rbc(:) = rb(:) - rc(:)
974 dab = norm2(rab)
975
976 DO iset = 1, nseta
977
978 ncoa = npgfa(iset)*ncoset(la_max(iset))
979 sgfa = first_sgfa(1, iset)
980
981 DO jset = 1, nsetb
982
983 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
984
985 ncob = npgfb(jset)*ncoset(lb_max(jset))
986 sgfb = first_sgfb(1, jset)
987
988 ! Calculate the primitive integrals
989 CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), &
990 rpgfa(:, iset), la_min(iset), &
991 lb_max(jset), npgfb(jset), zetb(:, jset), &
992 rpgfb(:, jset), rac, rbc, mab)
993
994 ! Contraction step
995 DO i = 1, nm
996 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
997 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
998 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
999 0.0_dp, work(1, 1), SIZE(work, 1))
1000
1001 IF (iatom <= jatom) THEN
1002 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1003 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1004 work(1, 1), SIZE(work, 1), &
1005 1.0_dp, mint(i)%block(sgfa, sgfb), &
1006 SIZE(mint(i)%block, 1))
1007 ELSE
1008 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1009 -1.0_dp, work(1, 1), SIZE(work, 1), &
1010 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1011 1.0_dp, mint(i)%block(sgfb, sgfa), &
1012 SIZE(mint(i)%block, 1))
1013 END IF
1014
1015 END DO
1016
1017 END DO
1018 END DO
1019 END associate
1020 END DO
1021 CALL neighbor_list_iterator_release(nl_iterator)
1022
1023 ! Release work storage
1024 DEALLOCATE (mab, basis_set_list)
1025 DEALLOCATE (work)
1026 DO i = 1, nm
1027 NULLIFY (mint(i)%block)
1028 END DO
1029 DEALLOCATE (mint)
1030
1031 CALL timestop(handle)
1032
1033 END SUBROUTINE build_local_magmom_matrix
1034
1035! **************************************************************************************************
1036!> \brief ...
1037!> \param qs_env ...
1038!> \param cosmat ...
1039!> \param sinmat ...
1040!> \param kvec ...
1041!> \param sab_orb_external ...
1042!> \param basis_type ...
1043! **************************************************************************************************
1044 SUBROUTINE build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
1045
1046 TYPE(qs_environment_type), POINTER :: qs_env
1047 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1048 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1049 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1050 OPTIONAL, POINTER :: sab_orb_external
1051 CHARACTER(len=*), OPTIONAL :: basis_type
1052
1053 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_moment_matrix'
1054
1055 INTEGER :: handle
1056 TYPE(dbcsr_p_type), DIMENSION(2, 1) :: op_sm_set
1057 REAL(kind=dp), DIMENSION(3, 1) :: kvec_batch
1058
1059 CALL timeset(routinen, handle)
1060
1061 op_sm_set(1, 1)%matrix => cosmat
1062 op_sm_set(2, 1)%matrix => sinmat
1063 kvec_batch(:, 1) = kvec(:)
1064
1065 CALL build_exp_ikr_matrix(qs_env, op_sm_set, kvec_batch, &
1066 sab_orb_external=sab_orb_external, basis_type=basis_type)
1067
1068 CALL timestop(handle)
1069
1070 END SUBROUTINE build_berry_moment_matrix
1071
1072! **************************************************************************************************
1073!> \brief ...
1074!> \param qs_env ...
1075!> \param cosmat ...
1076!> \param sinmat ...
1077!> \param kvec ...
1078!> \param basis_type ...
1079! **************************************************************************************************
1080 SUBROUTINE build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
1081
1082 TYPE(qs_environment_type), POINTER :: qs_env
1083 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: cosmat, sinmat
1084 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1085 CHARACTER(len=*), OPTIONAL :: basis_type
1086
1087 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_kpoint_matrix'
1088
1089 INTEGER :: handle, i, iatom, ic, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, &
1090 ldsa, ldsb, ldwork, natom, ncoa, ncob, nimg, nkind, nseta, nsetb, sgfa, sgfb
1091 INTEGER, DIMENSION(3) :: icell
1092 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
1093 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1094 LOGICAL :: found, use_cell_mapping
1095 REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1096 REAL(kind=dp) :: dab
1097 REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
1098 TYPE(cell_type), POINTER :: cell
1099 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1100 TYPE(dft_control_type), POINTER :: dft_control
1101 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1102 TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
1103 TYPE(kpoint_type), POINTER :: kpoints
1104 TYPE(neighbor_list_iterator_p_type), &
1105 DIMENSION(:), POINTER :: nl_iterator
1106 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1107 POINTER :: sab_orb
1108 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1109 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1110 TYPE(qs_kind_type), POINTER :: qs_kind
1111 TYPE(qs_ks_env_type), POINTER :: ks_env
1112
1113 CALL timeset(routinen, handle)
1114
1115 CALL get_qs_env(qs_env, &
1116 ks_env=ks_env, &
1117 dft_control=dft_control)
1118 nimg = dft_control%nimages
1119 IF (nimg > 1) THEN
1120 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1121 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1122 use_cell_mapping = .true.
1123 ELSE
1124 use_cell_mapping = .false.
1125 END IF
1126
1127 CALL get_qs_env(qs_env=qs_env, &
1128 qs_kind_set=qs_kind_set, &
1129 particle_set=particle_set, cell=cell, &
1130 sab_orb=sab_orb)
1131
1132 nkind = SIZE(qs_kind_set)
1133 natom = SIZE(particle_set)
1134 ALLOCATE (basis_set_list(nkind))
1135 DO ikind = 1, nkind
1136 qs_kind => qs_kind_set(ikind)
1137 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=basis_type)
1138 IF (ASSOCIATED(basis_set)) THEN
1139 basis_set_list(ikind)%gto_basis_set => basis_set
1140 ELSE
1141 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1142 END IF
1143 END DO
1144
1145 ALLOCATE (row_blk_sizes(natom))
1146 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1147 basis=basis_set_list)
1148 CALL get_ks_env(ks_env, dbcsr_dist=dbcsr_dist)
1149 ! (re)allocate matrix sets
1150 CALL dbcsr_allocate_matrix_set(sinmat, 1, nimg)
1151 CALL dbcsr_allocate_matrix_set(cosmat, 1, nimg)
1152 DO i = 1, nimg
1153 ! sin
1154 ALLOCATE (sinmat(1, i)%matrix)
1155 CALL dbcsr_create(matrix=sinmat(1, i)%matrix, &
1156 name="SINMAT", &
1157 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1158 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
1159 CALL cp_dbcsr_alloc_block_from_nbl(sinmat(1, i)%matrix, sab_orb)
1160 CALL dbcsr_set(sinmat(1, i)%matrix, 0.0_dp)
1161 ! cos
1162 ALLOCATE (cosmat(1, i)%matrix)
1163 CALL dbcsr_create(matrix=cosmat(1, i)%matrix, &
1164 name="COSMAT", &
1165 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1166 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
1167 CALL cp_dbcsr_alloc_block_from_nbl(cosmat(1, i)%matrix, sab_orb)
1168 CALL dbcsr_set(cosmat(1, i)%matrix, 0.0_dp)
1169 END DO
1170
1171 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
1172 ldab = ldwork
1173 ALLOCATE (cosab(ldab, ldab))
1174 ALLOCATE (sinab(ldab, ldab))
1175 ALLOCATE (work(ldwork, ldwork))
1176
1177 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1178 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1179 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1180 iatom=iatom, jatom=jatom, r=rab, cell=icell)
1181 basis_set_a => basis_set_list(ikind)%gto_basis_set
1182 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1183 basis_set_b => basis_set_list(jkind)%gto_basis_set
1184 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1185 associate( &
1186 ! basis ikind
1187 first_sgfa => basis_set_a%first_sgf, &
1188 la_max => basis_set_a%lmax, &
1189 la_min => basis_set_a%lmin, &
1190 npgfa => basis_set_a%npgf, &
1191 nsgfa => basis_set_a%nsgf_set, &
1192 rpgfa => basis_set_a%pgf_radius, &
1193 set_radius_a => basis_set_a%set_radius, &
1194 sphi_a => basis_set_a%sphi, &
1195 zeta => basis_set_a%zet, &
1196 ! basis jkind, &
1197 first_sgfb => basis_set_b%first_sgf, &
1198 lb_max => basis_set_b%lmax, &
1199 lb_min => basis_set_b%lmin, &
1200 npgfb => basis_set_b%npgf, &
1201 nsgfb => basis_set_b%nsgf_set, &
1202 rpgfb => basis_set_b%pgf_radius, &
1203 set_radius_b => basis_set_b%set_radius, &
1204 sphi_b => basis_set_b%sphi, &
1205 zetb => basis_set_b%zet)
1206
1207 nseta = basis_set_a%nset
1208 nsetb = basis_set_b%nset
1209
1210 ldsa = SIZE(sphi_a, 1)
1211 ldsb = SIZE(sphi_b, 1)
1212
1213 IF (iatom <= jatom) THEN
1214 irow = iatom
1215 icol = jatom
1216 ELSE
1217 irow = jatom
1218 icol = iatom
1219 END IF
1220
1221 IF (use_cell_mapping) THEN
1222 ic = cell_to_index(icell(1), icell(2), icell(3))
1223 cpassert(ic > 0)
1224 ELSE
1225 ic = 1
1226 END IF
1227
1228 NULLIFY (sblock)
1229 CALL dbcsr_get_block_p(matrix=sinmat(1, ic)%matrix, &
1230 row=irow, col=icol, block=sblock, found=found)
1231 cpassert(found)
1232 NULLIFY (cblock)
1233 CALL dbcsr_get_block_p(matrix=cosmat(1, ic)%matrix, &
1234 row=irow, col=icol, block=cblock, found=found)
1235 cpassert(found)
1236
1237 ra(:) = pbc(particle_set(iatom)%r(:), cell)
1238 rb(:) = ra + rab
1239 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1240
1241 DO iset = 1, nseta
1242
1243 ncoa = npgfa(iset)*ncoset(la_max(iset))
1244 sgfa = first_sgfa(1, iset)
1245
1246 DO jset = 1, nsetb
1247
1248 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1249
1250 ncob = npgfb(jset)*ncoset(lb_max(jset))
1251 sgfb = first_sgfb(1, jset)
1252
1253 ! Calculate the primitive integrals
1254 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1255 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1256 ra, rb, kvec, cosab, sinab)
1257 CALL contract_cossin(cblock, sblock, &
1258 iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1259 jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1260 cosab, sinab, ldab, work, ldwork)
1261
1262 END DO
1263 END DO
1264 END associate
1265 END DO
1266 CALL neighbor_list_iterator_release(nl_iterator)
1267
1268 DEALLOCATE (cosab)
1269 DEALLOCATE (sinab)
1270 DEALLOCATE (work)
1271 DEALLOCATE (basis_set_list)
1272 DEALLOCATE (row_blk_sizes)
1273
1274 CALL timestop(handle)
1275
1276 END SUBROUTINE build_berry_kpoint_matrix
1277
1278! **************************************************************************************************
1279!> \brief ...
1280!> \param qs_env ...
1281!> \param magnetic ...
1282!> \param nmoments ...
1283!> \param reference ...
1284!> \param ref_point ...
1285!> \param unit_number ...
1286! **************************************************************************************************
1287 SUBROUTINE qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number)
1288
1289 TYPE(qs_environment_type), POINTER :: qs_env
1290 LOGICAL, INTENT(IN) :: magnetic
1291 INTEGER, INTENT(IN) :: nmoments, reference
1292 REAL(dp), DIMENSION(:), POINTER :: ref_point
1293 INTEGER, INTENT(IN) :: unit_number
1294
1295 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_berry_phase'
1296
1297 CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
1298 CHARACTER(LEN=default_string_length) :: description
1299 COMPLEX(dp) :: xphase(3), zdet, zdeta, zi(3), &
1300 zij(3, 3), zijk(3, 3, 3), &
1301 zijkl(3, 3, 3, 3), zphase(3), zz
1302 INTEGER :: handle, i, ia, idim, ikind, ispin, ix, &
1303 iy, iz, j, k, l, nao, nm, nmo, nmom, &
1304 nmotot, tmp_dim
1305 LOGICAL :: floating, ghost, uniform
1306 REAL(dp) :: charge, ci(3), cij(3, 3), dd, occ, trace
1307 REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom
1308 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
1309 REAL(dp), DIMENSION(3) :: kvec, qq, rcc, ria
1310 TYPE(atomic_kind_type), POINTER :: atomic_kind
1311 TYPE(cell_type), POINTER :: cell
1312 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat
1313 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
1314 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: opvec
1315 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: op_fm_set
1316 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
1317 TYPE(cp_fm_type), POINTER :: mo_coeff
1318 TYPE(cp_result_type), POINTER :: results
1319 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
1320 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1321 TYPE(dft_control_type), POINTER :: dft_control
1322 TYPE(distribution_1d_type), POINTER :: local_particles
1323 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1324 TYPE(mp_para_env_type), POINTER :: para_env
1325 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1326 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1327 TYPE(qs_rho_type), POINTER :: rho
1328 TYPE(rt_prop_type), POINTER :: rtp
1329
1330 cpassert(ASSOCIATED(qs_env))
1331
1332 IF (ASSOCIATED(qs_env%ls_scf_env)) THEN
1333 IF (unit_number > 0) WRITE (unit_number, *) "Periodic moment calculation not implemented in linear scaling code"
1334 RETURN
1335 END IF
1336
1337 CALL timeset(routinen, handle)
1338
1339 ! restrict maximum moment available
1340 nmom = min(nmoments, 2)
1341
1342 nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
1343 ! rmom(:,1)=electronic
1344 ! rmom(:,2)=nuclear
1345 ! rmom(:,1)=total
1346 ALLOCATE (rmom(nm + 1, 3))
1347 ALLOCATE (rlab(nm + 1))
1348 rmom = 0.0_dp
1349 rlab = ""
1350 IF (magnetic) THEN
1351 nm = 3
1352 ALLOCATE (mmom(nm))
1353 mmom = 0._dp
1354 END IF
1355
1356 NULLIFY (dft_control, rho, cell, particle_set, results, para_env, &
1357 local_particles, matrix_s, mos, rho_ao)
1358
1359 CALL get_qs_env(qs_env, &
1360 dft_control=dft_control, &
1361 rho=rho, &
1362 cell=cell, &
1363 results=results, &
1364 particle_set=particle_set, &
1365 qs_kind_set=qs_kind_set, &
1366 para_env=para_env, &
1367 local_particles=local_particles, &
1368 matrix_s=matrix_s, &
1369 mos=mos)
1370
1371 CALL qs_rho_get(rho, rho_ao=rho_ao)
1372
1373 NULLIFY (cosmat, sinmat)
1374 ALLOCATE (cosmat, sinmat)
1375 CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
1376 CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
1377
1378 ALLOCATE (op_fm_set(2, dft_control%nspins))
1379 ALLOCATE (opvec(dft_control%nspins))
1380 ALLOCATE (eigrmat(dft_control%nspins))
1381 nmotot = 0
1382 DO ispin = 1, dft_control%nspins
1383 NULLIFY (tmp_fm_struct, mo_coeff)
1384 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
1385 nmotot = nmotot + nmo
1386 CALL cp_fm_create(opvec(ispin), mo_coeff%matrix_struct)
1387 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
1388 ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
1389 DO i = 1, SIZE(op_fm_set, 1)
1390 CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
1391 END DO
1392 CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
1393 CALL cp_fm_struct_release(tmp_fm_struct)
1394 END DO
1395
1396 ! occupation
1397 DO ispin = 1, dft_control%nspins
1398 CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
1399 IF (.NOT. uniform) THEN
1400 cpwarn("Berry phase moments for non uniform MOs' occupation numbers not implemented")
1401 END IF
1402 END DO
1403
1404 ! reference point
1405 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
1406 rcc = pbc(rcc, cell)
1407
1408 ! label
1409 DO l = 1, nm
1410 ix = indco(1, l + 1)
1411 iy = indco(2, l + 1)
1412 iz = indco(3, l + 1)
1413 CALL set_label(rlab(l + 1), ix, iy, iz)
1414 END DO
1415
1416 ! nuclear contribution
1417 DO ia = 1, SIZE(particle_set)
1418 atomic_kind => particle_set(ia)%atomic_kind
1419 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1420 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1421 IF (.NOT. ghost .AND. .NOT. floating) THEN
1422 rmom(1, 2) = rmom(1, 2) - charge
1423 END IF
1424 END DO
1425 ria = twopi*matmul(cell%h_inv, rcc)
1426 zphase = cmplx(cos(ria), sin(ria), dp)**rmom(1, 2)
1427
1428 zi = 0._dp
1429 zij = 0._dp
1430 zijk = 0._dp
1431 zijkl = 0._dp
1432
1433 DO l = 1, nmom
1434 SELECT CASE (l)
1435 CASE (1)
1436 ! Dipole
1437 zi(:) = cmplx(1._dp, 0._dp, dp)
1438 DO ia = 1, SIZE(particle_set)
1439 atomic_kind => particle_set(ia)%atomic_kind
1440 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1441 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1442 IF (.NOT. ghost .AND. .NOT. floating) THEN
1443 ria = particle_set(ia)%r
1444 ria = pbc(ria, cell)
1445 DO i = 1, 3
1446 kvec(:) = twopi*cell%h_inv(i, :)
1447 dd = sum(kvec(:)*ria(:))
1448 zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
1449 zi(i) = zi(i)*zdeta
1450 END DO
1451 END IF
1452 END DO
1453 zi = zi*zphase
1454 ci = aimag(log(zi))/twopi
1455 qq = aimag(log(zi))
1456 rmom(2:4, 2) = matmul(cell%hmat, ci)
1457 CASE (2)
1458 ! Quadrupole
1459 cpabort("Berry phase moments bigger than 1 not implemented")
1460 zij(:, :) = cmplx(1._dp, 0._dp, dp)
1461 DO ia = 1, SIZE(particle_set)
1462 atomic_kind => particle_set(ia)%atomic_kind
1463 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1464 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
1465 ria = particle_set(ia)%r
1466 ria = pbc(ria, cell)
1467 DO i = 1, 3
1468 DO j = i, 3
1469 kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1470 dd = sum(kvec(:)*ria(:))
1471 zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
1472 zij(i, j) = zij(i, j)*zdeta
1473 zij(j, i) = zij(i, j)
1474 END DO
1475 END DO
1476 END DO
1477 DO i = 1, 3
1478 DO j = 1, 3
1479 zij(i, j) = zij(i, j)*zphase(i)*zphase(j)
1480 zz = zij(i, j)/zi(i)/zi(j)
1481 cij(i, j) = aimag(log(zz))/twopi
1482 END DO
1483 END DO
1484 cij = 0.5_dp*cij/twopi/twopi
1485 cij = matmul(matmul(cell%hmat, cij), transpose(cell%hmat))
1486 DO k = 4, 9
1487 ix = indco(1, k + 1)
1488 iy = indco(2, k + 1)
1489 iz = indco(3, k + 1)
1490 IF (ix == 0) THEN
1491 rmom(k + 1, 2) = cij(iy, iz)
1492 ELSE IF (iy == 0) THEN
1493 rmom(k + 1, 2) = cij(ix, iz)
1494 ELSE IF (iz == 0) THEN
1495 rmom(k + 1, 2) = cij(ix, iy)
1496 END IF
1497 END DO
1498 CASE (3)
1499 ! Octapole
1500 cpabort("Berry phase moments bigger than 2 not implemented")
1501 CASE (4)
1502 ! Hexadecapole
1503 cpabort("Berry phase moments bigger than 3 not implemented")
1504 CASE DEFAULT
1505 cpabort("Berry phase moments bigger than 4 not implemented")
1506 END SELECT
1507 END DO
1508
1509 ! electronic contribution
1510
1511 ria = twopi*real(nmotot, dp)*occ*matmul(cell%h_inv, rcc)
1512 xphase = cmplx(cos(ria), sin(ria), dp)
1513
1514 ! charge
1515 trace = 0.0_dp
1516 DO ispin = 1, dft_control%nspins
1517 CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
1518 rmom(1, 1) = rmom(1, 1) + trace
1519 END DO
1520
1521 zi = 0._dp
1522 zij = 0._dp
1523 zijk = 0._dp
1524 zijkl = 0._dp
1525
1526 DO l = 1, nmom
1527 SELECT CASE (l)
1528 CASE (1)
1529 ! Dipole
1530 DO i = 1, 3
1531 kvec(:) = twopi*cell%h_inv(i, :)
1532 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1533 IF (qs_env%run_rtp) THEN
1534 CALL get_qs_env(qs_env, rtp=rtp)
1535 CALL get_rtp(rtp, mos_new=mos_new)
1536 CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1537 ELSE
1538 CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1539 END IF
1540 zdet = cmplx(1._dp, 0._dp, dp)
1541 DO ispin = 1, dft_control%nspins
1542 CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
1543 DO idim = 1, tmp_dim
1544 eigrmat(ispin)%local_data(:, idim) = &
1545 cmplx(op_fm_set(1, ispin)%local_data(:, idim), &
1546 -op_fm_set(2, ispin)%local_data(:, idim), dp)
1547 END DO
1548 ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
1549 CALL cp_cfm_det(eigrmat(ispin), zdeta)
1550 zdet = zdet*zdeta
1551 IF (dft_control%nspins == 1) THEN
1552 zdet = zdet*zdeta
1553 END IF
1554 END DO
1555 zi(i) = zdet
1556 END DO
1557 zi = zi*xphase
1558 CASE (2)
1559 ! Quadrupole
1560 cpabort("Berry phase moments bigger than 1 not implemented")
1561 DO i = 1, 3
1562 DO j = i, 3
1563 kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1564 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1565 IF (qs_env%run_rtp) THEN
1566 CALL get_qs_env(qs_env, rtp=rtp)
1567 CALL get_rtp(rtp, mos_new=mos_new)
1568 CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1569 ELSE
1570 CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1571 END IF
1572 zdet = cmplx(1._dp, 0._dp, dp)
1573 DO ispin = 1, dft_control%nspins
1574 CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
1575 DO idim = 1, tmp_dim
1576 eigrmat(ispin)%local_data(:, idim) = &
1577 cmplx(op_fm_set(1, ispin)%local_data(:, idim), &
1578 -op_fm_set(2, ispin)%local_data(:, idim), dp)
1579 END DO
1580 ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
1581 CALL cp_cfm_det(eigrmat(ispin), zdeta)
1582 zdet = zdet*zdeta
1583 IF (dft_control%nspins == 1) THEN
1584 zdet = zdet*zdeta
1585 END IF
1586 END DO
1587 zij(i, j) = zdet*xphase(i)*xphase(j)
1588 zij(j, i) = zdet*xphase(i)*xphase(j)
1589 END DO
1590 END DO
1591 CASE (3)
1592 ! Octapole
1593 cpabort("Berry phase moments bigger than 2 not implemented")
1594 CASE (4)
1595 ! Hexadecapole
1596 cpabort("Berry phase moments bigger than 3 not implemented")
1597 CASE DEFAULT
1598 cpabort("Berry phase moments bigger than 4 not implemented")
1599 END SELECT
1600 END DO
1601 DO l = 1, nmom
1602 SELECT CASE (l)
1603 CASE (1)
1604 ! Dipole (apply periodic (2 Pi) boundary conditions)
1605 ci = aimag(log(zi))
1606 DO i = 1, 3
1607 IF (qq(i) + ci(i) > pi) ci(i) = ci(i) - twopi
1608 IF (qq(i) + ci(i) < -pi) ci(i) = ci(i) + twopi
1609 END DO
1610 rmom(2:4, 1) = matmul(cell%hmat, ci)/twopi
1611 CASE (2)
1612 ! Quadrupole
1613 cpabort("Berry phase moments bigger than 1 not implemented")
1614 DO i = 1, 3
1615 DO j = 1, 3
1616 zz = zij(i, j)/zi(i)/zi(j)
1617 cij(i, j) = aimag(log(zz))/twopi
1618 END DO
1619 END DO
1620 cij = 0.5_dp*cij/twopi/twopi
1621 cij = matmul(matmul(cell%hmat, cij), transpose(cell%hmat))
1622 DO k = 4, 9
1623 ix = indco(1, k + 1)
1624 iy = indco(2, k + 1)
1625 iz = indco(3, k + 1)
1626 IF (ix == 0) THEN
1627 rmom(k + 1, 1) = cij(iy, iz)
1628 ELSE IF (iy == 0) THEN
1629 rmom(k + 1, 1) = cij(ix, iz)
1630 ELSE IF (iz == 0) THEN
1631 rmom(k + 1, 1) = cij(ix, iy)
1632 END IF
1633 END DO
1634 CASE (3)
1635 ! Octapole
1636 cpabort("Berry phase moments bigger than 2 not implemented")
1637 CASE (4)
1638 ! Hexadecapole
1639 cpabort("Berry phase moments bigger than 3 not implemented")
1640 CASE DEFAULT
1641 cpabort("Berry phase moments bigger than 4 not implemented")
1642 END SELECT
1643 END DO
1644
1645 rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
1646 description = "[DIPOLE]"
1647 CALL cp_results_erase(results=results, description=description)
1648 CALL put_results(results=results, description=description, &
1649 values=rmom(2:4, 3))
1650 IF (magnetic) THEN
1651 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.true., mmom=mmom)
1652 ELSE
1653 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.true.)
1654 END IF
1655
1656 DEALLOCATE (rmom)
1657 DEALLOCATE (rlab)
1658 IF (magnetic) THEN
1659 DEALLOCATE (mmom)
1660 END IF
1661
1662 CALL dbcsr_deallocate_matrix(cosmat)
1663 CALL dbcsr_deallocate_matrix(sinmat)
1664
1665 CALL cp_fm_release(opvec)
1666 CALL cp_fm_release(op_fm_set)
1667 DO ispin = 1, dft_control%nspins
1668 CALL cp_cfm_release(eigrmat(ispin))
1669 END DO
1670 DEALLOCATE (eigrmat)
1671
1672 CALL timestop(handle)
1673
1674 END SUBROUTINE qs_moment_berry_phase
1675
1676! **************************************************************************************************
1677!> \brief ...
1678!> \param cosmat ...
1679!> \param sinmat ...
1680!> \param mos ...
1681!> \param op_fm_set ...
1682!> \param opvec ...
1683! **************************************************************************************************
1684 SUBROUTINE op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1685
1686 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1687 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1688 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
1689 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: opvec
1690
1691 INTEGER :: i, nao, nmo
1692 TYPE(cp_fm_type), POINTER :: mo_coeff
1693
1694 DO i = 1, SIZE(op_fm_set, 2) ! spin
1695 CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
1696 CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff, opvec(i), ncol=nmo)
1697 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
1698 op_fm_set(1, i))
1699 CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff, opvec(i), ncol=nmo)
1700 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
1701 op_fm_set(2, i))
1702 END DO
1703
1704 END SUBROUTINE op_orbbas
1705
1706! **************************************************************************************************
1707!> \brief ...
1708!> \param cosmat ...
1709!> \param sinmat ...
1710!> \param mos ...
1711!> \param op_fm_set ...
1712!> \param mos_new ...
1713! **************************************************************************************************
1714 SUBROUTINE op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1715
1716 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1717 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1718 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
1719 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
1720
1721 INTEGER :: i, icol, lcol, nao, newdim, nmo
1722 LOGICAL :: double_col, double_row
1723 TYPE(cp_fm_struct_type), POINTER :: newstruct, newstruct1
1724 TYPE(cp_fm_type) :: work, work1, work2
1725 TYPE(cp_fm_type), POINTER :: mo_coeff
1726
1727 DO i = 1, SIZE(op_fm_set, 2) ! spin
1728 CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
1729 CALL cp_fm_get_info(mos_new(2*i), ncol_local=lcol, ncol_global=nmo)
1730 double_col = .true.
1731 double_row = .false.
1732 CALL cp_fm_struct_double(newstruct, &
1733 mos_new(2*i)%matrix_struct, &
1734 mos_new(2*i)%matrix_struct%context, &
1735 double_col, &
1736 double_row)
1737
1738 CALL cp_fm_create(work, matrix_struct=newstruct)
1739 CALL cp_fm_create(work1, matrix_struct=newstruct)
1740 CALL cp_fm_create(work2, matrix_struct=newstruct)
1741 CALL cp_fm_get_info(work, ncol_global=newdim)
1742
1743 CALL cp_fm_set_all(work, 0.0_dp, 0.0_dp)
1744 DO icol = 1, lcol
1745 work%local_data(:, icol) = mos_new(2*i - 1)%local_data(:, icol)
1746 work%local_data(:, icol + lcol) = mos_new(2*i)%local_data(:, icol)
1747 END DO
1748
1749 CALL cp_dbcsr_sm_fm_multiply(cosmat, work, work1, ncol=newdim)
1750 CALL cp_dbcsr_sm_fm_multiply(sinmat, work, work2, ncol=newdim)
1751
1752 DO icol = 1, lcol
1753 work%local_data(:, icol) = work1%local_data(:, icol) - work2%local_data(:, icol + lcol)
1754 work%local_data(:, icol + lcol) = work1%local_data(:, icol + lcol) + work2%local_data(:, icol)
1755 END DO
1756
1757 CALL cp_fm_release(work1)
1758 CALL cp_fm_release(work2)
1759
1760 CALL cp_fm_struct_double(newstruct1, &
1761 op_fm_set(1, i)%matrix_struct, &
1762 op_fm_set(1, i)%matrix_struct%context, &
1763 double_col, &
1764 double_row)
1765
1766 CALL cp_fm_create(work1, matrix_struct=newstruct1)
1767
1768 CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i - 1), &
1769 work, 0.0_dp, work1)
1770
1771 DO icol = 1, lcol
1772 op_fm_set(1, i)%local_data(:, icol) = work1%local_data(:, icol)
1773 op_fm_set(2, i)%local_data(:, icol) = work1%local_data(:, icol + lcol)
1774 END DO
1775
1776 CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i), &
1777 work, 0.0_dp, work1)
1778
1779 DO icol = 1, lcol
1780 op_fm_set(1, i)%local_data(:, icol) = &
1781 op_fm_set(1, i)%local_data(:, icol) + work1%local_data(:, icol + lcol)
1782 op_fm_set(2, i)%local_data(:, icol) = &
1783 op_fm_set(2, i)%local_data(:, icol) - work1%local_data(:, icol)
1784 END DO
1785
1786 CALL cp_fm_release(work)
1787 CALL cp_fm_release(work1)
1788 CALL cp_fm_struct_release(newstruct)
1789 CALL cp_fm_struct_release(newstruct1)
1790
1791 END DO
1792
1793 END SUBROUTINE op_orbbas_rtp
1794
1795! **************************************************************************************************
1796!> \brief ...
1797!> \param qs_env ...
1798!> \param magnetic ...
1799!> \param nmoments ...
1800!> \param reference ...
1801!> \param ref_point ...
1802!> \param unit_number ...
1803!> \param vel_reprs ...
1804!> \param com_nl ...
1805! **************************************************************************************************
1806 SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
1807
1808 TYPE(qs_environment_type), POINTER :: qs_env
1809 LOGICAL, INTENT(IN) :: magnetic
1810 INTEGER, INTENT(IN) :: nmoments, reference
1811 REAL(dp), DIMENSION(:), INTENT(IN), POINTER :: ref_point
1812 INTEGER, INTENT(IN) :: unit_number
1813 LOGICAL, INTENT(IN), OPTIONAL :: vel_reprs, com_nl
1814
1815 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_locop'
1816
1817 CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
1818 CHARACTER(LEN=default_string_length) :: description
1819 INTEGER :: akind, handle, i, ia, iatom, idir, &
1820 ikind, ispin, ix, iy, iz, l, nm, nmom, &
1821 order
1822 LOGICAL :: my_com_nl, my_velreprs
1823 REAL(dp) :: charge, dd, strace, trace
1824 REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rrv_vrr, &
1825 nlcom_rv, nlcom_rvr, nlcom_rxrv, &
1826 qupole_der, rmom_vel
1827 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
1828 REAL(dp), DIMENSION(3) :: rcc, ria
1829 TYPE(atomic_kind_type), POINTER :: atomic_kind
1830 TYPE(cell_type), POINTER :: cell
1831 TYPE(cp_result_type), POINTER :: results
1832 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom, matrix_s, moments, momentum, &
1833 rho_ao
1834 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_der
1835 TYPE(dbcsr_type), POINTER :: tmp_ao
1836 TYPE(dft_control_type), POINTER :: dft_control
1837 TYPE(distribution_1d_type), POINTER :: local_particles
1838 TYPE(mp_para_env_type), POINTER :: para_env
1839 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1840 POINTER :: sab_all, sab_orb
1841 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1842 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1843 TYPE(qs_rho_type), POINTER :: rho
1844
1845 cpassert(ASSOCIATED(qs_env))
1846
1847 CALL timeset(routinen, handle)
1848
1849 my_velreprs = .false.
1850 IF (PRESENT(vel_reprs)) my_velreprs = vel_reprs
1851 IF (PRESENT(com_nl)) my_com_nl = com_nl
1852 IF (my_velreprs) CALL cite_reference(mattiat2019)
1853
1854 ! reference point
1855 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
1856
1857 ! only allow for moments up to maxl set by basis
1858 nmom = min(nmoments, current_maxl)
1859 ! electronic contribution
1860 NULLIFY (dft_control, rho, cell, particle_set, qs_kind_set, results, para_env, matrix_s, rho_ao, sab_all, sab_orb)
1861 CALL get_qs_env(qs_env, &
1862 dft_control=dft_control, &
1863 rho=rho, &
1864 cell=cell, &
1865 results=results, &
1866 particle_set=particle_set, &
1867 qs_kind_set=qs_kind_set, &
1868 para_env=para_env, &
1869 matrix_s=matrix_s, &
1870 sab_all=sab_all, &
1871 sab_orb=sab_orb)
1872
1873 IF (my_com_nl) THEN
1874 IF ((nmom >= 1) .AND. my_velreprs) THEN
1875 ALLOCATE (nlcom_rv(3))
1876 nlcom_rv(:) = 0._dp
1877 END IF
1878 IF ((nmom >= 2) .AND. my_velreprs) THEN
1879 ALLOCATE (nlcom_rrv(6))
1880 nlcom_rrv(:) = 0._dp
1881 ALLOCATE (nlcom_rvr(6))
1882 nlcom_rvr(:) = 0._dp
1883 ALLOCATE (nlcom_rrv_vrr(6))
1884 nlcom_rrv_vrr(:) = 0._dp
1885 END IF
1886 IF (magnetic) THEN
1887 ALLOCATE (nlcom_rxrv(3))
1888 nlcom_rxrv = 0._dp
1889 END IF
1890 ! Calculate non local correction terms
1891 CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr, rcc)
1892 END IF
1893
1894 NULLIFY (moments)
1895 nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
1896 CALL dbcsr_allocate_matrix_set(moments, nm)
1897 DO i = 1, nm
1898 ALLOCATE (moments(i)%matrix)
1899 IF (my_velreprs .AND. (nmom >= 2)) THEN
1900 CALL dbcsr_create(moments(i)%matrix, template=matrix_s(1)%matrix, &
1901 matrix_type=dbcsr_type_symmetric)
1902 CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_orb)
1903 ELSE
1904 CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
1905 END IF
1906 CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
1907 END DO
1908
1909 ! calculate derivatives if quadrupole in vel. reprs. is requested
1910 IF (my_velreprs .AND. (nmom >= 2)) THEN
1911 NULLIFY (moments_der)
1912 CALL dbcsr_allocate_matrix_set(moments_der, 3, 3)
1913 DO i = 1, 3 ! x, y, z
1914 DO idir = 1, 3 ! d/dx, d/dy, d/dz
1915 CALL dbcsr_init_p(moments_der(i, idir)%matrix)
1916 CALL dbcsr_create(moments_der(i, idir)%matrix, template=matrix_s(1)%matrix, &
1917 matrix_type=dbcsr_type_antisymmetric)
1918 CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_orb)
1919 CALL dbcsr_set(moments_der(i, idir)%matrix, 0.0_dp)
1920 END DO
1921 END DO
1922 CALL build_local_moments_der_matrix(qs_env, moments_der, 1, 2, ref_point=rcc, moments=moments)
1923 ELSE
1924 CALL build_local_moment_matrix(qs_env, moments, nmom, ref_point=rcc)
1925 END IF
1926
1927 CALL qs_rho_get(rho, rho_ao=rho_ao)
1928
1929 ALLOCATE (rmom(nm + 1, 3))
1930 ALLOCATE (rlab(nm + 1))
1931 rmom = 0.0_dp
1932 rlab = ""
1933
1934 IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
1935 ! Allocate matrix to store the matrix product to be traced (dbcsr_dot only works for products of
1936 ! symmetric matrices)
1937 NULLIFY (tmp_ao)
1938 CALL dbcsr_init_p(tmp_ao)
1939 CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
1940 CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
1941 CALL dbcsr_set(tmp_ao, 0.0_dp)
1942 END IF
1943
1944 trace = 0.0_dp
1945 DO ispin = 1, dft_control%nspins
1946 CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
1947 rmom(1, 1) = rmom(1, 1) + trace
1948 END DO
1949
1950 DO i = 1, SIZE(moments)
1951 strace = 0._dp
1952 DO ispin = 1, dft_control%nspins
1953 IF (my_velreprs .AND. nmoments >= 2) THEN
1954 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, moments(i)%matrix, &
1955 0.0_dp, tmp_ao)
1956 CALL dbcsr_trace(tmp_ao, trace)
1957 ELSE
1958 CALL dbcsr_dot(rho_ao(ispin)%matrix, moments(i)%matrix, trace)
1959 END IF
1960 strace = strace + trace
1961 END DO
1962 rmom(i + 1, 1) = strace
1963 END DO
1964
1965 CALL dbcsr_deallocate_matrix_set(moments)
1966
1967 ! nuclear contribution
1968 CALL get_qs_env(qs_env=qs_env, &
1969 local_particles=local_particles)
1970 DO ikind = 1, SIZE(local_particles%n_el)
1971 DO ia = 1, local_particles%n_el(ikind)
1972 iatom = local_particles%list(ikind)%array(ia)
1973 ! fold atomic positions back into unit cell
1974 ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
1975 ria = ria - rcc
1976 atomic_kind => particle_set(iatom)%atomic_kind
1977 CALL get_atomic_kind(atomic_kind, kind_number=akind)
1978 CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
1979 rmom(1, 2) = rmom(1, 2) - charge
1980 DO l = 1, nm
1981 ix = indco(1, l + 1)
1982 iy = indco(2, l + 1)
1983 iz = indco(3, l + 1)
1984 dd = 1._dp
1985 IF (ix > 0) dd = dd*ria(1)**ix
1986 IF (iy > 0) dd = dd*ria(2)**iy
1987 IF (iz > 0) dd = dd*ria(3)**iz
1988 rmom(l + 1, 2) = rmom(l + 1, 2) - charge*dd
1989 CALL set_label(rlab(l + 1), ix, iy, iz)
1990 END DO
1991 END DO
1992 END DO
1993 CALL para_env%sum(rmom(:, 2))
1994 rmom(:, :) = -rmom(:, :)
1995 rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
1996
1997 ! magnetic moments
1998 IF (magnetic) THEN
1999 NULLIFY (magmom)
2000 CALL dbcsr_allocate_matrix_set(magmom, 3)
2001 DO i = 1, SIZE(magmom)
2002 CALL dbcsr_init_p(magmom(i)%matrix)
2003 CALL dbcsr_create(magmom(i)%matrix, template=matrix_s(1)%matrix, &
2004 matrix_type=dbcsr_type_antisymmetric)
2005 CALL cp_dbcsr_alloc_block_from_nbl(magmom(i)%matrix, sab_orb)
2006 CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
2007 END DO
2008
2009 CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
2010
2011 ALLOCATE (mmom(SIZE(magmom)))
2012 mmom(:) = 0.0_dp
2013 IF (qs_env%run_rtp) THEN
2014 ! get imaginary part of the density in rho_ao (the real part is not needed since the trace of the product
2015 ! of a symmetric (REAL(rho_ao)) and an anti-symmetric (L_AO) matrix is zero)
2016 ! There may be other cases, where the imaginary part of the density is relevant
2017 NULLIFY (rho_ao)
2018 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2019 END IF
2020 ! if the density is purely real this is an expensive way to calculate zero
2021 DO i = 1, SIZE(magmom)
2022 strace = 0._dp
2023 DO ispin = 1, dft_control%nspins
2024 CALL dbcsr_set(tmp_ao, 0.0_dp)
2025 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, magmom(i)%matrix, &
2026 0.0_dp, tmp_ao)
2027 CALL dbcsr_trace(tmp_ao, trace)
2028 strace = strace + trace
2029 END DO
2030 mmom(i) = strace
2031 END DO
2032
2033 CALL dbcsr_deallocate_matrix_set(magmom)
2034 END IF
2035
2036 ! velocity representations
2037 IF (my_velreprs) THEN
2038 ALLOCATE (rmom_vel(nm))
2039 rmom_vel = 0.0_dp
2040
2041 DO order = 1, nmom
2042 SELECT CASE (order)
2043
2044 CASE (1) ! expectation value of momentum
2045 NULLIFY (momentum)
2046 CALL dbcsr_allocate_matrix_set(momentum, 3)
2047 DO i = 1, 3
2048 CALL dbcsr_init_p(momentum(i)%matrix)
2049 CALL dbcsr_create(momentum(i)%matrix, template=matrix_s(1)%matrix, &
2050 matrix_type=dbcsr_type_antisymmetric)
2051 CALL cp_dbcsr_alloc_block_from_nbl(momentum(i)%matrix, sab_orb)
2052 CALL dbcsr_set(momentum(i)%matrix, 0.0_dp)
2053 END DO
2054 CALL build_lin_mom_matrix(qs_env, momentum)
2055
2056 ! imaginary part of the density for RTP, real part gives 0 since momentum is antisymmetric
2057 IF (qs_env%run_rtp) THEN
2058 NULLIFY (rho_ao)
2059 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2060 DO idir = 1, SIZE(momentum)
2061 strace = 0._dp
2062 DO ispin = 1, dft_control%nspins
2063 CALL dbcsr_set(tmp_ao, 0.0_dp)
2064 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, momentum(idir)%matrix, &
2065 0.0_dp, tmp_ao)
2066 CALL dbcsr_trace(tmp_ao, trace)
2067 strace = strace + trace
2068 END DO
2069 rmom_vel(idir) = rmom_vel(idir) + strace
2070 END DO
2071 END IF
2072
2073 CALL dbcsr_deallocate_matrix_set(momentum)
2074
2075 CASE (2) ! expectation value of quadrupole moment in vel. reprs.
2076 ALLOCATE (qupole_der(9)) ! will contain the expectation values of r_\alpha * d/d r_\beta
2077 qupole_der = 0._dp
2078
2079 NULLIFY (rho_ao)
2080 CALL qs_rho_get(rho, rho_ao=rho_ao)
2081
2082 ! Calculate expectation value over real part
2083 trace = 0._dp
2084 DO i = 1, 3
2085 DO idir = 1, 3
2086 strace = 0._dp
2087 DO ispin = 1, dft_control%nspins
2088 CALL dbcsr_set(tmp_ao, 0._dp)
2089 CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2090 CALL dbcsr_trace(tmp_ao, trace)
2091 strace = strace + trace
2092 END DO
2093 qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2094 END DO
2095 END DO
2096
2097 IF (qs_env%run_rtp) THEN
2098 NULLIFY (rho_ao)
2099 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2100
2101 ! Calculate expectation value over imaginary part
2102 trace = 0._dp
2103 DO i = 1, 3
2104 DO idir = 1, 3
2105 strace = 0._dp
2106 DO ispin = 1, dft_control%nspins
2107 CALL dbcsr_set(tmp_ao, 0._dp)
2108 CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2109 CALL dbcsr_trace(tmp_ao, trace)
2110 strace = strace + trace
2111 END DO
2112 qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2113 END DO
2114 END DO
2115 END IF
2116
2117 ! calculate vel. reprs. of quadrupole moment from derivatives
2118 rmom_vel(4) = -2*qupole_der(1) - rmom(1, 1)
2119 rmom_vel(5) = -qupole_der(2) - qupole_der(4)
2120 rmom_vel(6) = -qupole_der(3) - qupole_der(7)
2121 rmom_vel(7) = -2*qupole_der(5) - rmom(1, 1)
2122 rmom_vel(8) = -qupole_der(6) - qupole_der(8)
2123 rmom_vel(9) = -2*qupole_der(9) - rmom(1, 1)
2124
2125 DEALLOCATE (qupole_der)
2126 CASE DEFAULT
2127 END SELECT
2128 END DO
2129 END IF
2130
2131 IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2132 CALL dbcsr_deallocate_matrix(tmp_ao)
2133 END IF
2134 IF (my_velreprs .AND. (nmoments >= 2)) THEN
2135 CALL dbcsr_deallocate_matrix_set(moments_der)
2136 END IF
2137
2138 description = "[DIPOLE]"
2139 CALL cp_results_erase(results=results, description=description)
2140 CALL put_results(results=results, description=description, &
2141 values=rmom(2:4, 3))
2142
2143 IF (magnetic .AND. my_velreprs) THEN
2144 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., mmom=mmom, rmom_vel=rmom_vel)
2145 ELSE IF (magnetic .AND. .NOT. my_velreprs) THEN
2146 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., mmom=mmom)
2147 ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2148 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., rmom_vel=rmom_vel)
2149 ELSE
2150 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false.)
2151 END IF
2152
2153 IF (my_com_nl) THEN
2154 IF (magnetic) THEN
2155 mmom(:) = nlcom_rxrv(:)
2156 END IF
2157 IF (my_velreprs .AND. (nmom >= 1)) THEN
2158 DEALLOCATE (rmom_vel)
2159 ALLOCATE (rmom_vel(21))
2160 rmom_vel(1:3) = nlcom_rv
2161 END IF
2162 IF (my_velreprs .AND. (nmom >= 2)) THEN
2163 rmom_vel(4:9) = nlcom_rrv
2164 rmom_vel(10:15) = nlcom_rvr
2165 rmom_vel(16:21) = nlcom_rrv_vrr
2166 END IF
2167 IF (magnetic .AND. .NOT. my_velreprs) THEN
2168 CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom)
2169 ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2170 CALL print_moments_nl(unit_number, nmom, rlab, rmom_vel=rmom_vel)
2171 ELSE IF (my_velreprs .AND. magnetic) THEN
2172 CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom, rmom_vel=rmom_vel)
2173 END IF
2174
2175 END IF
2176
2177 IF (my_com_nl) THEN
2178 IF (nmom >= 1 .AND. my_velreprs) DEALLOCATE (nlcom_rv)
2179 IF (nmom >= 2 .AND. my_velreprs) THEN
2180 DEALLOCATE (nlcom_rrv)
2181 DEALLOCATE (nlcom_rvr)
2182 DEALLOCATE (nlcom_rrv_vrr)
2183 END IF
2184 IF (magnetic) DEALLOCATE (nlcom_rxrv)
2185 END IF
2186
2187 DEALLOCATE (rmom)
2188 DEALLOCATE (rlab)
2189 IF (magnetic) THEN
2190 DEALLOCATE (mmom)
2191 END IF
2192 IF (my_velreprs) THEN
2193 DEALLOCATE (rmom_vel)
2194 END IF
2195
2196 CALL timestop(handle)
2197
2198 END SUBROUTINE qs_moment_locop
2199
2200! **************************************************************************************************
2201!> \brief ...
2202!> \param label ...
2203!> \param ix ...
2204!> \param iy ...
2205!> \param iz ...
2206! **************************************************************************************************
2207 SUBROUTINE set_label(label, ix, iy, iz)
2208 CHARACTER(LEN=*), INTENT(OUT) :: label
2209 INTEGER, INTENT(IN) :: ix, iy, iz
2210
2211 INTEGER :: i
2212
2213 label = ""
2214 DO i = 1, ix
2215 WRITE (label(i:), "(A1)") "X"
2216 END DO
2217 DO i = ix + 1, ix + iy
2218 WRITE (label(i:), "(A1)") "Y"
2219 END DO
2220 DO i = ix + iy + 1, ix + iy + iz
2221 WRITE (label(i:), "(A1)") "Z"
2222 END DO
2223
2224 END SUBROUTINE set_label
2225
2226! **************************************************************************************************
2227!> \brief ...
2228!> \param unit_number ...
2229!> \param nmom ...
2230!> \param rmom ...
2231!> \param rlab ...
2232!> \param rcc ...
2233!> \param cell ...
2234!> \param periodic ...
2235!> \param mmom ...
2236!> \param rmom_vel ...
2237! **************************************************************************************************
2238 SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmom, rmom_vel)
2239 INTEGER, INTENT(IN) :: unit_number, nmom
2240 REAL(dp), DIMENSION(:, :), INTENT(IN) :: rmom
2241 CHARACTER(LEN=8), DIMENSION(:) :: rlab
2242 REAL(dp), DIMENSION(3), INTENT(IN) :: rcc
2243 TYPE(cell_type), POINTER :: cell
2244 LOGICAL :: periodic
2245 REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2246
2247 INTEGER :: i, i0, i1, j, l
2248 REAL(dp) :: dd
2249
2250 IF (unit_number > 0) THEN
2251 DO l = 0, nmom
2252 SELECT CASE (l)
2253 CASE (0)
2254 WRITE (unit_number, "(T3,A,T33,3F16.8)") "Reference Point [Bohr]", rcc
2255 WRITE (unit_number, "(T3,A)") "Charges"
2256 WRITE (unit_number, "(T5,A,T18,F14.8,T36,A,T42,F14.8,T60,A,T67,F14.8)") &
2257 "Electronic=", rmom(1, 1), "Core=", rmom(1, 2), "Total=", rmom(1, 3)
2258 CASE (1)
2259 IF (periodic) THEN
2260 WRITE (unit_number, "(T3,A)") "Dipole vectors are based on the periodic (Berry phase) operator."
2261 WRITE (unit_number, "(T3,A)") "They are defined modulo integer multiples of the cell matrix [Debye]."
2262 WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[X] [", cell%hmat(1, :)*debye, "] [i]"
2263 WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Y]=[", cell%hmat(2, :)*debye, "]*[j]"
2264 WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Z] [", cell%hmat(3, :)*debye, "] [k]"
2265 ELSE
2266 WRITE (unit_number, "(T3,A)") "Dipoles are based on the traditional operator."
2267 END IF
2268 dd = sqrt(sum(rmom(2:4, 3)**2))*debye
2269 WRITE (unit_number, "(T3,A)") "Dipole moment [Debye]"
2270 WRITE (unit_number, "(T5,3(A,A,E15.7,1X),T60,A,T68,F13.7)") &
2271 (trim(rlab(i)), "=", rmom(i, 3)*debye, i=2, 4), "Total=", dd
2272 CASE (2)
2273 WRITE (unit_number, "(T3,A)") "Quadrupole moment [Debye*Angstrom]"
2274 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2275 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=5, 7)
2276 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2277 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=8, 10)
2278 CASE (3)
2279 WRITE (unit_number, "(T3,A)") "Octapole moment [Debye*Angstrom**2]"
2280 WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2281 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=11, 14)
2282 WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2283 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=15, 18)
2284 WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2285 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=19, 20)
2286 CASE (4)
2287 WRITE (unit_number, "(T3,A)") "Hexadecapole moment [Debye*Angstrom**3]"
2288 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2289 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=21, 24)
2290 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2291 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=25, 28)
2292 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2293 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=29, 32)
2294 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2295 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=32, 35)
2296 CASE DEFAULT
2297 WRITE (unit_number, "(T3,A,A,I2)") "Higher moment [Debye*Angstrom**(L-1)]", &
2298 " L=", l
2299 i0 = (6 + 11*(l - 1) + 6*(l - 1)**2 + (l - 1)**3)/6
2300 i1 = (6 + 11*l + 6*l**2 + l**3)/6 - 1
2301 dd = debye/(bohr)**(l - 1)
2302 DO i = i0, i1, 3
2303 WRITE (unit_number, "(T18,3(A,A,F14.8,4X))") &
2304 (trim(rlab(j + 1)), "=", rmom(j + 1, 3)*dd, j=i, min(i1, i + 2))
2305 END DO
2306 END SELECT
2307 END DO
2308 IF (PRESENT(mmom)) THEN
2309 IF (nmom >= 1) THEN
2310 dd = sqrt(sum(mmom(1:3)**2))
2311 WRITE (unit_number, "(T3,A)") "Orbital angular momentum [a. u.]"
2312 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2313 (trim(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2314 END IF
2315 END IF
2316 IF (PRESENT(rmom_vel)) THEN
2317 DO l = 1, nmom
2318 SELECT CASE (l)
2319 CASE (1)
2320 dd = sqrt(sum(rmom_vel(1:3)**2))
2321 WRITE (unit_number, "(T3,A)") "Expectation value of momentum operator [a. u.]"
2322 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2323 (trim(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2324 CASE (2)
2325 WRITE (unit_number, "(T3,A)") "Expectation value of quadrupole operator in vel. repr. [a. u.]"
2326 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2327 (trim(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2328 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2329 (trim(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2330 CASE DEFAULT
2331 END SELECT
2332 END DO
2333 END IF
2334 END IF
2335
2336 END SUBROUTINE print_moments
2337
2338! **************************************************************************************************
2339!> \brief ...
2340!> \param unit_number ...
2341!> \param nmom ...
2342!> \param rlab ...
2343!> \param mmom ...
2344!> \param rmom_vel ...
2345! **************************************************************************************************
2346 SUBROUTINE print_moments_nl(unit_number, nmom, rlab, mmom, rmom_vel)
2347 INTEGER, INTENT(IN) :: unit_number, nmom
2348 CHARACTER(LEN=8), DIMENSION(:) :: rlab
2349 REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2350
2351 INTEGER :: i, l
2352 REAL(dp) :: dd
2353
2354 IF (unit_number > 0) THEN
2355 IF (PRESENT(mmom)) THEN
2356 IF (nmom >= 1) THEN
2357 dd = sqrt(sum(mmom(1:3)**2))
2358 WRITE (unit_number, "(T3,A)") "Expectation value of rx[r,V_nl] [a. u.]"
2359 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2360 (trim(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2361 END IF
2362 END IF
2363 IF (PRESENT(rmom_vel)) THEN
2364 DO l = 1, nmom
2365 SELECT CASE (l)
2366 CASE (1)
2367 dd = sqrt(sum(rmom_vel(1:3)**2))
2368 WRITE (unit_number, "(T3,A)") "Expectation value of [r,V_nl] [a. u.]"
2369 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2370 (trim(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2371 CASE (2)
2372 WRITE (unit_number, "(T3,A)") "Expectation value of [rr,V_nl] [a. u.]"
2373 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2374 (trim(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2375 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2376 (trim(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2377 WRITE (unit_number, "(T3,A)") "Expectation value of r x V_nl x r [a. u.]"
2378 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2379 (trim(rlab(i + 1 - 6)), "=", rmom_vel(i), i=10, 12)
2380 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2381 (trim(rlab(i + 1 - 6)), "=", rmom_vel(i), i=13, 15)
2382 WRITE (unit_number, "(T3,A)") "Expectation value of r x r x V_nl + V_nl x r x r [a. u.]"
2383 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2384 (trim(rlab(i + 1 - 12)), "=", rmom_vel(i), i=16, 18)
2385 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2386 (trim(rlab(i + 1 - 12)), "=", rmom_vel(i), i=19, 21)
2387 CASE DEFAULT
2388 END SELECT
2389 END DO
2390 END IF
2391 END IF
2392
2393 END SUBROUTINE print_moments_nl
2394
2395! **************************************************************************************************
2396!> \brief Calculate the expectation value of operators related to non-local potential:
2397!> [r, Vnl], noted rv
2398!> r x [r,Vnl], noted rxrv
2399!> [rr,Vnl], noted rrv
2400!> r x Vnl x r, noted rvr
2401!> r x r x Vnl + Vnl x r x r, noted rrv_vrr
2402!> Note that the 3 first operator are commutator while the 2 last
2403!> are not. For reading clarity the same notation is used for all 5
2404!> operators.
2405!> \param qs_env ...
2406!> \param nlcom_rv ...
2407!> \param nlcom_rxrv ...
2408!> \param nlcom_rrv ...
2409!> \param nlcom_rvr ...
2410!> \param nlcom_rrv_vrr ...
2411!> \param ref_point ...
2412! **************************************************************************************************
2413 SUBROUTINE calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, &
2414 nlcom_rrv_vrr, ref_point)
2415
2416 TYPE(qs_environment_type), POINTER :: qs_env
2417 REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: nlcom_rv, nlcom_rxrv, nlcom_rrv, &
2418 nlcom_rvr, nlcom_rrv_vrr
2419 REAL(dp), DIMENSION(3) :: ref_point
2420
2421 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_commutator_nl_terms'
2422
2423 INTEGER :: handle, ind, ispin
2424 LOGICAL :: calc_rrv, calc_rrv_vrr, calc_rv, &
2425 calc_rvr, calc_rxrv
2426 REAL(dp) :: eps_ppnl, strace, trace
2427 TYPE(cell_type), POINTER :: cell
2428 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rrv, matrix_rrv_vrr, matrix_rv, &
2429 matrix_rvr, matrix_rxrv, matrix_s, &
2430 rho_ao
2431 TYPE(dbcsr_type), POINTER :: tmp_ao
2432 TYPE(dft_control_type), POINTER :: dft_control
2433 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2434 POINTER :: sab_all, sab_orb, sap_ppnl
2435 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2436 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2437 TYPE(qs_rho_type), POINTER :: rho
2438
2439 CALL timeset(routinen, handle)
2440
2441 calc_rv = .false.
2442 calc_rxrv = .false.
2443 calc_rrv = .false.
2444 calc_rvr = .false.
2445 calc_rrv_vrr = .false.
2446
2447 ! rv, rxrv and rrv are commutator matrices: anti-symmetric.
2448 ! The real part of the density matrix rho_ao is symmetric so that
2449 ! the expectation value of real density matrix is zero. Hence, if
2450 ! the density matrix is real, no need to compute these quantities.
2451 ! This is not the case for rvr and rrv_vrr which are symmetric.
2452
2453 IF (ALLOCATED(nlcom_rv)) THEN
2454 nlcom_rv(:) = 0._dp
2455 IF (qs_env%run_rtp) calc_rv = .true.
2456 END IF
2457 IF (ALLOCATED(nlcom_rxrv)) THEN
2458 nlcom_rxrv(:) = 0._dp
2459 IF (qs_env%run_rtp) calc_rxrv = .true.
2460 END IF
2461 IF (ALLOCATED(nlcom_rrv)) THEN
2462 nlcom_rrv(:) = 0._dp
2463 IF (qs_env%run_rtp) calc_rrv = .true.
2464 END IF
2465 IF (ALLOCATED(nlcom_rvr)) THEN
2466 nlcom_rvr(:) = 0._dp
2467 calc_rvr = .true.
2468 END IF
2469 IF (ALLOCATED(nlcom_rrv_vrr)) THEN
2470 nlcom_rrv_vrr(:) = 0._dp
2471 calc_rrv_vrr = .true.
2472 END IF
2473
2474 IF (.NOT. (calc_rv .OR. calc_rrv .OR. calc_rxrv .OR. calc_rvr .OR. calc_rrv_vrr)) THEN
2475 CALL timestop(handle)
2476 RETURN
2477 END IF
2478
2479 NULLIFY (cell, matrix_s, particle_set, qs_kind_set, rho, sab_all, sab_orb, sap_ppnl)
2480 CALL get_qs_env(qs_env, &
2481 cell=cell, &
2482 dft_control=dft_control, &
2483 matrix_s=matrix_s, &
2484 particle_set=particle_set, &
2485 qs_kind_set=qs_kind_set, &
2486 rho=rho, &
2487 sab_orb=sab_orb, &
2488 sab_all=sab_all, &
2489 sap_ppnl=sap_ppnl)
2490
2491 eps_ppnl = dft_control%qs_control%eps_ppnl
2492
2493 ! Allocate storage
2494 NULLIFY (matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr)
2495 IF (calc_rv) THEN
2496 CALL dbcsr_allocate_matrix_set(matrix_rv, 3)
2497 DO ind = 1, 3
2498 CALL dbcsr_init_p(matrix_rv(ind)%matrix)
2499 CALL dbcsr_create(matrix_rv(ind)%matrix, template=matrix_s(1)%matrix, &
2500 matrix_type=dbcsr_type_antisymmetric)
2501 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(ind)%matrix, sab_orb)
2502 CALL dbcsr_set(matrix_rv(ind)%matrix, 0._dp)
2503 END DO
2504 END IF
2505
2506 IF (calc_rxrv) THEN
2507 CALL dbcsr_allocate_matrix_set(matrix_rxrv, 3)
2508 DO ind = 1, 3
2509 CALL dbcsr_init_p(matrix_rxrv(ind)%matrix)
2510 CALL dbcsr_create(matrix_rxrv(ind)%matrix, template=matrix_s(1)%matrix, &
2511 matrix_type=dbcsr_type_antisymmetric)
2512 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rxrv(ind)%matrix, sab_orb)
2513 CALL dbcsr_set(matrix_rxrv(ind)%matrix, 0._dp)
2514 END DO
2515 END IF
2516
2517 IF (calc_rrv) THEN
2518 CALL dbcsr_allocate_matrix_set(matrix_rrv, 6)
2519 DO ind = 1, 6
2520 CALL dbcsr_init_p(matrix_rrv(ind)%matrix)
2521 CALL dbcsr_create(matrix_rrv(ind)%matrix, template=matrix_s(1)%matrix, &
2522 matrix_type=dbcsr_type_antisymmetric)
2523 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv(ind)%matrix, sab_orb)
2524 CALL dbcsr_set(matrix_rrv(ind)%matrix, 0._dp)
2525 END DO
2526 END IF
2527
2528 IF (calc_rvr) THEN
2529 CALL dbcsr_allocate_matrix_set(matrix_rvr, 6)
2530 DO ind = 1, 6
2531 CALL dbcsr_init_p(matrix_rvr(ind)%matrix)
2532 CALL dbcsr_create(matrix_rvr(ind)%matrix, template=matrix_s(1)%matrix, &
2533 matrix_type=dbcsr_type_symmetric)
2534 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rvr(ind)%matrix, sab_orb)
2535 CALL dbcsr_set(matrix_rvr(ind)%matrix, 0._dp)
2536 END DO
2537 END IF
2538 IF (calc_rrv_vrr) THEN
2539 CALL dbcsr_allocate_matrix_set(matrix_rrv_vrr, 6)
2540 DO ind = 1, 6
2541 CALL dbcsr_init_p(matrix_rrv_vrr(ind)%matrix)
2542 CALL dbcsr_create(matrix_rrv_vrr(ind)%matrix, template=matrix_s(1)%matrix, &
2543 matrix_type=dbcsr_type_symmetric)
2544 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv_vrr(ind)%matrix, sab_orb)
2545 CALL dbcsr_set(matrix_rrv_vrr(ind)%matrix, 0._dp)
2546 END DO
2547 END IF
2548
2549 ! calculate evaluation of operators in AO basis set
2550 CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
2551 matrix_rxrv=matrix_rxrv, matrix_rrv=matrix_rrv, matrix_rvr=matrix_rvr, &
2552 matrix_rrv_vrr=matrix_rrv_vrr, ref_point=ref_point)
2553
2554 ! Calculate expectation values
2555 ! Real part
2556 NULLIFY (tmp_ao)
2557 CALL dbcsr_init_p(tmp_ao)
2558 CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
2559 CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
2560 CALL dbcsr_set(tmp_ao, 0.0_dp)
2561
2562 IF (calc_rvr .OR. calc_rrv_vrr) THEN
2563 NULLIFY (rho_ao)
2564 CALL qs_rho_get(rho, rho_ao=rho_ao)
2565
2566 IF (calc_rvr) THEN
2567 trace = 0._dp
2568 DO ind = 1, SIZE(matrix_rvr)
2569 strace = 0._dp
2570 DO ispin = 1, dft_control%nspins
2571 CALL dbcsr_set(tmp_ao, 0.0_dp)
2572 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rvr(ind)%matrix, &
2573 0.0_dp, tmp_ao)
2574 CALL dbcsr_trace(tmp_ao, trace)
2575 strace = strace + trace
2576 END DO
2577 nlcom_rvr(ind) = nlcom_rvr(ind) + strace
2578 END DO
2579 END IF
2580
2581 IF (calc_rrv_vrr) THEN
2582 trace = 0._dp
2583 DO ind = 1, SIZE(matrix_rrv_vrr)
2584 strace = 0._dp
2585 DO ispin = 1, dft_control%nspins
2586 CALL dbcsr_set(tmp_ao, 0.0_dp)
2587 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv_vrr(ind)%matrix, &
2588 0.0_dp, tmp_ao)
2589 CALL dbcsr_trace(tmp_ao, trace)
2590 strace = strace + trace
2591 END DO
2592 nlcom_rrv_vrr(ind) = nlcom_rrv_vrr(ind) + strace
2593 END DO
2594 END IF
2595 END IF
2596
2597 ! imagninary part of the density matrix
2598 NULLIFY (rho_ao)
2599 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2600
2601 IF (calc_rv) THEN
2602 trace = 0._dp
2603 DO ind = 1, SIZE(matrix_rv)
2604 strace = 0._dp
2605 DO ispin = 1, dft_control%nspins
2606 CALL dbcsr_set(tmp_ao, 0.0_dp)
2607 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rv(ind)%matrix, &
2608 0.0_dp, tmp_ao)
2609 CALL dbcsr_trace(tmp_ao, trace)
2610 strace = strace + trace
2611 END DO
2612 nlcom_rv(ind) = nlcom_rv(ind) + strace
2613 END DO
2614 END IF
2615
2616 IF (calc_rrv) THEN
2617 trace = 0._dp
2618 DO ind = 1, SIZE(matrix_rrv)
2619 strace = 0._dp
2620 DO ispin = 1, dft_control%nspins
2621 CALL dbcsr_set(tmp_ao, 0.0_dp)
2622 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv(ind)%matrix, &
2623 0.0_dp, tmp_ao)
2624 CALL dbcsr_trace(tmp_ao, trace)
2625 strace = strace + trace
2626 END DO
2627 nlcom_rrv(ind) = nlcom_rrv(ind) + strace
2628 END DO
2629 END IF
2630
2631 IF (calc_rxrv) THEN
2632 trace = 0._dp
2633 DO ind = 1, SIZE(matrix_rxrv)
2634 strace = 0._dp
2635 DO ispin = 1, dft_control%nspins
2636 CALL dbcsr_set(tmp_ao, 0.0_dp)
2637 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rxrv(ind)%matrix, &
2638 0.0_dp, tmp_ao)
2639 CALL dbcsr_trace(tmp_ao, trace)
2640 strace = strace + trace
2641 END DO
2642 nlcom_rxrv(ind) = nlcom_rxrv(ind) + strace
2643 END DO
2644 END IF
2645 CALL dbcsr_deallocate_matrix(tmp_ao)
2646 IF (calc_rv) CALL dbcsr_deallocate_matrix_set(matrix_rv)
2647 IF (calc_rxrv) CALL dbcsr_deallocate_matrix_set(matrix_rxrv)
2648 IF (calc_rrv) CALL dbcsr_deallocate_matrix_set(matrix_rrv)
2649 IF (calc_rvr) CALL dbcsr_deallocate_matrix_set(matrix_rvr)
2650 IF (calc_rrv_vrr) CALL dbcsr_deallocate_matrix_set(matrix_rrv_vrr)
2651
2652 CALL timestop(handle)
2653 END SUBROUTINE calculate_commutator_nl_terms
2654
2655! **************************************************************************************************
2656!> \brief Get list of kpoints from input to compute dipole moment elements
2657!> \param qs_env ...
2658!> \param xkp ...
2659!> \param special_pnts ...
2660!> \author Shridhar Shanbhag
2661! **************************************************************************************************
2662 SUBROUTINE get_xkp_for_dipole_calc(qs_env, xkp, special_pnts)
2663 TYPE(qs_environment_type), POINTER :: qs_env
2664 TYPE(section_vals_type), POINTER :: kpnts, kpset
2665 REAL(kind=dp), DIMENSION(3, 3) :: cart_hmat, hmat
2666
2667 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_xkp_for_dipole_calc'
2668
2669 CHARACTER(LEN=default_string_length) :: ustr
2670 TYPE(kpoint_type), POINTER :: kpoint_work
2671 TYPE(cell_type), POINTER :: cell
2672 CHARACTER(LEN=default_string_length), &
2673 DIMENSION(:), POINTER :: strptr
2674 CHARACTER(LEN=default_string_length), &
2675 DIMENSION(:), POINTER :: special_pnts, spname
2676 CHARACTER(LEN=max_line_length) :: error_message
2677 INTEGER :: handle, i, ik, ikk, ip, &
2678 n_ptr, npline, nkp
2679 LOGICAL :: explicit_kpnts, explicit_kpset
2680 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: kspecial, xkp
2681 REAL(kind=dp), DIMENSION(3) :: kpptr
2682 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2683
2684 CALL timeset(routinen, handle)
2685 kpset => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINT_SET")
2686 kpnts => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINTS")
2687 CALL section_vals_get(kpset, explicit=explicit_kpset)
2688 CALL section_vals_get(kpnts, explicit=explicit_kpnts)
2689 IF (explicit_kpset .AND. explicit_kpnts) then
2690 cpabort("Both KPOINT_SET and KPOINTS present in MOMENTS section")
2691 end if
2692
2693 IF (explicit_kpset) THEN
2694 CALL get_qs_env(qs_env, cell=cell)
2695 CALL get_cell(cell, h=hmat)
2696 cart_hmat(:, :) = hmat(:, :)
2697 IF (cell%input_cell_canonicalized) cart_hmat(:, :) = cell%input_hmat(:, :)
2698 CALL section_vals_val_get(kpset, "NPOINTS", i_val=npline)
2699 CALL section_vals_val_get(kpset, "UNITS", c_val=ustr)
2700 CALL uppercase(ustr)
2701 CALL section_vals_val_get(kpset, "SPECIAL_POINT", n_rep_val=n_ptr)
2702 cpassert(n_ptr > 0)
2703 ALLOCATE (kspecial(3, n_ptr))
2704 ALLOCATE (spname(n_ptr))
2705 DO ip = 1, n_ptr
2706 CALL section_vals_val_get(kpset, "SPECIAL_POINT", i_rep_val=ip, c_vals=strptr)
2707 IF (SIZE(strptr(:), 1) == 4) THEN
2708 spname(ip) = strptr(1)
2709 DO i = 1, 3
2710 CALL read_float_object(strptr(i + 1), kpptr(i), error_message)
2711 IF (len_trim(error_message) > 0) cpabort(trim(error_message))
2712 END DO
2713 ELSE IF (SIZE(strptr(:), 1) == 3) THEN
2714 spname(ip) = "not specified"
2715 DO i = 1, 3
2716 CALL read_float_object(strptr(i), kpptr(i), error_message)
2717 IF (len_trim(error_message) > 0) cpabort(trim(error_message))
2718 END DO
2719 ELSE
2720 cpabort("Input SPECIAL_POINT invalid")
2721 END IF
2722 SELECT CASE (ustr)
2723 CASE ("B_VECTOR")
2724 kspecial(1:3, ip) = kpptr(1:3)
2725 CASE ("CART_ANGSTROM")
2726 kspecial(1:3, ip) = (kpptr(1)*cart_hmat(1, 1:3) + &
2727 kpptr(2)*cart_hmat(2, 1:3) + &
2728 kpptr(3)*cart_hmat(3, 1:3))/twopi*angstrom
2729 CASE ("CART_BOHR")
2730 kspecial(1:3, ip) = (kpptr(1)*cart_hmat(1, 1:3) + &
2731 kpptr(2)*cart_hmat(2, 1:3) + &
2732 kpptr(3)*cart_hmat(3, 1:3))/twopi
2733 CASE DEFAULT
2734 cpabort("Unknown unit <"//trim(ustr)//"> specified for k-point definition")
2735 END SELECT
2736 END DO
2737 nkp = (n_ptr - 1)*npline + 1
2738 cpassert(nkp >= 1)
2739
2740 ! Initialize environment and calculate MOs
2741 ALLOCATE (xkp(3, nkp))
2742 ALLOCATE (special_pnts(nkp))
2743 special_pnts(:) = ""
2744 xkp(1:3, 1) = kspecial(1:3, 1)
2745 ikk = 1
2746 special_pnts(ikk) = spname(1)
2747 DO ik = 2, n_ptr
2748 DO ip = 1, npline
2749 ikk = ikk + 1
2750 xkp(1:3, ikk) = kspecial(1:3, ik - 1) + &
2751 REAL(ip, kind=dp)/real(npline, kind=dp)* &
2752 (kspecial(1:3, ik) - kspecial(1:3, ik - 1))
2753 END DO
2754 special_pnts(ikk) = spname(ik)
2755 END DO
2756 DEALLOCATE (spname, kspecial)
2757 ELSE IF (explicit_kpnts) THEN
2758 CALL get_qs_env(qs_env, particle_set=particle_set, cell=cell)
2759 CALL get_cell(cell, h=hmat)
2760 NULLIFY (kpoint_work)
2761 CALL kpoint_create(kpoint_work)
2762 CALL read_kpoint_section(kpoint_work, kpnts, hmat, cell)
2763 CALL kpoint_initialize(kpoint_work, particle_set, cell)
2764 nkp = kpoint_work%nkp
2765 ALLOCATE (xkp(3, nkp))
2766 ALLOCATE (special_pnts(nkp))
2767 special_pnts(:) = ""
2768 xkp(1:3, :) = kpoint_work%xkp(1:3, :)
2769 CALL kpoint_release(kpoint_work)
2770 ELSE
2771 ! use k-point mesh from DFT calculation
2772 CALL get_qs_env(qs_env, kpoints=kpoint_work)
2773 nkp = kpoint_work%nkp
2774 nkp = kpoint_work%nkp
2775 ALLOCATE (xkp(3, nkp))
2776 ALLOCATE (special_pnts(nkp))
2777 special_pnts(:) = ""
2778 xkp(1:3, :) = kpoint_work%xkp(1:3, :)
2779 END IF
2780 CALL timestop(handle)
2781
2782 END SUBROUTINE get_xkp_for_dipole_calc
2783! **************************************************************************************************
2784!> \brief Calculate local moment matrix for a periodic system for all image cells
2785!> \param qs_env ...
2786!> \param moments_rs_img ...
2787!> \param rcc ...
2788!> \author Shridhar Shanbhag
2789! **************************************************************************************************
2790 SUBROUTINE build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc)
2791
2792 TYPE(qs_environment_type), POINTER :: qs_env
2793 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_rs_img
2794 REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
2795
2796 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moment_matrix_rs_img'
2797
2798 INTEGER :: handle, i_dir, iatom, ic, ikind, iset, j, jatom, jkind, jset, &
2799 ldsa, ldsb, ldwork, ncoa, ncob, nimg, nkind, nseta, nsetb, nsize, sgfa, sgfb
2800 INTEGER, DIMENSION(3) :: icell
2801 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
2802 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2803 LOGICAL :: found
2804 REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
2805 REAL(dp), DIMENSION(:, :), POINTER :: dblock, work
2806 REAL(dp), DIMENSION(:, :, :), POINTER :: dipab
2807 REAL(kind=dp) :: dab
2808 TYPE(cell_type), POINTER :: cell
2809 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
2810 TYPE(dft_control_type), POINTER :: dft_control
2811 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
2812 TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
2813 TYPE(kpoint_type), POINTER :: kpoints_all
2814 TYPE(mp_para_env_type), POINTER :: para_env
2815 TYPE(neighbor_list_iterator_p_type), &
2816 DIMENSION(:), POINTER :: nl_iterator
2817 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2818 POINTER :: sab_all
2819 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2820 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2821 TYPE(qs_kind_type), POINTER :: qs_kind
2822
2823 CALL timeset(routinen, handle)
2824
2825 CALL get_qs_env(qs_env=qs_env, &
2826 dft_control=dft_control, &
2827 qs_kind_set=qs_kind_set, &
2828 matrix_ks_kp=matrix_ks_kp, &
2829 particle_set=particle_set, &
2830 cell=cell, &
2831 para_env=para_env, &
2832 sab_all=sab_all)
2833
2834 NULLIFY (kpoints_all)
2835 CALL kpoint_create(kpoints_all)
2836 CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, nimg)
2837
2838 nkind = SIZE(qs_kind_set)
2839 ALLOCATE (basis_set_list(nkind))
2840 DO ikind = 1, nkind
2841 qs_kind => qs_kind_set(ikind)
2842 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set)
2843 IF (ASSOCIATED(basis_set)) THEN
2844 basis_set_list(ikind)%gto_basis_set => basis_set
2845 ELSE
2846 NULLIFY (basis_set_list(ikind)%gto_basis_set)
2847 END IF
2848 END DO
2849
2850 rc(:) = 0._dp
2851 IF (PRESENT(rcc)) rc(:) = rcc(:)
2852
2853 CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_list)
2854 CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
2855 nsize = SIZE(index_to_cell, 2)
2856 cpassert(SIZE(moments_rs_img, 2) == nsize)
2857 DO i_dir = 1, 3
2858 DO j = 1, nsize
2859 ALLOCATE (moments_rs_img(i_dir, j)%matrix)
2860 CALL dbcsr_create(matrix=moments_rs_img(i_dir, j)%matrix, &
2861 template=matrix_ks_kp(1, 1)%matrix, &
2862 matrix_type=dbcsr_type_no_symmetry, &
2863 name="DIPMAT")
2864 CALL cp_dbcsr_alloc_block_from_nbl(moments_rs_img(i_dir, j)%matrix, sab_all)
2865 CALL dbcsr_set(moments_rs_img(i_dir, j)%matrix, 0.0_dp)
2866 END DO
2867 END DO
2868
2869 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
2870 ALLOCATE (dipab(ldwork, ldwork, 3))
2871 ALLOCATE (work(ldwork, ldwork))
2872
2873 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
2874 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2875 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2876 iatom=iatom, jatom=jatom, r=rab, cell=icell)
2877
2878 basis_set_a => basis_set_list(ikind)%gto_basis_set
2879 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
2880 basis_set_b => basis_set_list(jkind)%gto_basis_set
2881 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
2882 associate( &
2883 ! basis ikind
2884 first_sgfa => basis_set_a%first_sgf, &
2885 la_max => basis_set_a%lmax, &
2886 la_min => basis_set_a%lmin, &
2887 npgfa => basis_set_a%npgf, &
2888 nsgfa => basis_set_a%nsgf_set, &
2889 rpgfa => basis_set_a%pgf_radius, &
2890 set_radius_a => basis_set_a%set_radius, &
2891 sphi_a => basis_set_a%sphi, &
2892 zeta => basis_set_a%zet, &
2893 ! basis jkind, &
2894 first_sgfb => basis_set_b%first_sgf, &
2895 lb_max => basis_set_b%lmax, &
2896 lb_min => basis_set_b%lmin, &
2897 npgfb => basis_set_b%npgf, &
2898 nsgfb => basis_set_b%nsgf_set, &
2899 rpgfb => basis_set_b%pgf_radius, &
2900 set_radius_b => basis_set_b%set_radius, &
2901 sphi_b => basis_set_b%sphi, &
2902 zetb => basis_set_b%zet)
2903
2904 nseta = basis_set_a%nset
2905 nsetb = basis_set_b%nset
2906
2907 ldsa = SIZE(sphi_a, 1)
2908 ldsb = SIZE(sphi_b, 1)
2909
2910 NULLIFY (dblock)
2911
2912 ra = pbc(particle_set(iatom)%r(:), cell)
2913 rb(:) = ra(:) + rab(:)
2914 rac = ra - rc
2915 rbc = rb - rc
2916 dab = norm2(rab)
2917
2918 ic = cell_to_index(icell(1), icell(2), icell(3))
2919
2920 DO iset = 1, nseta
2921
2922 ncoa = npgfa(iset)*ncoset(la_max(iset))
2923 sgfa = first_sgfa(1, iset)
2924
2925 DO jset = 1, nsetb
2926
2927 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
2928
2929 ncob = npgfb(jset)*ncoset(lb_max(jset))
2930 sgfb = first_sgfb(1, jset)
2931
2932 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
2933 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), 1, &
2934 rac, rbc, dipab)
2935 DO i_dir = 1, 3
2936 CALL dbcsr_get_block_p(matrix=moments_rs_img(i_dir, ic)%matrix, &
2937 row=iatom, col=jatom, block=dblock, found=found)
2938 cpassert(found)
2939 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
2940 1.0_dp, dipab(1, 1, i_dir), ldwork, &
2941 sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
2942
2943 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
2944 1.0_dp, sphi_a(1, sgfa), ldsa, &
2945 work(1, 1), ldwork, 1.0_dp, dblock(1, 1), SIZE(dblock, 1))
2946 END DO
2947 END DO
2948 END DO
2949 END associate
2950 END DO
2951 CALL neighbor_list_iterator_release(nl_iterator)
2952 CALL kpoint_release(kpoints_all)
2953 DEALLOCATE (dipab, work, basis_set_list)
2954 CALL timestop(handle)
2955
2956 END SUBROUTINE build_local_moment_matrix_rs_img
2957
2958! **************************************************************************************************
2959!> \brief Calculates the dipole moments and berry curvature for periodic systems for kpoints
2960!> \param qs_env ...
2961!> \param xkp list of kpoints
2962!> \param dipole ...
2963!> \param rcc coordinates about which to calculate the dipole
2964!> \param berry_c berry curvature calculated using Ω^γ_n = Σ_m 2*Im[d^α_nm (d^β_mn)*]
2965!> \param do_parallel option to distribute the result in dipole across
2966!> different MPI ranks
2967!> \author Shridhar Shanbhag
2968! **************************************************************************************************
2969 SUBROUTINE qs_moment_kpoints_deep(qs_env, xkp, dipole, rcc, berry_c, do_parallel)
2970 TYPE(qs_environment_type), POINTER :: qs_env
2971 LOGICAL, OPTIONAL :: do_parallel
2972 LOGICAL :: my_do_parallel, calc_bc
2973 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_kpoints_deep'
2974 COMPLEX(KIND=dp) :: phase, tmp_max
2975 COMPLEX(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: c_k, h_k, s_k, d_k, cdc, c_dh_c, &
2976 c_ds_c, dh_dk_i, ds_dk_i
2977 COMPLEX(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dip
2978 COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
2979 ALLOCATABLE :: dipole
2980 INTEGER :: handle, i_dir, ikp, nkp, &
2981 n_img_scf, n_img_all, nao, &
2982 num_pe, num_copy, mepos, n, m, mu, &
2983 ispin, nspin
2984 INTEGER, DIMENSION(3) :: periodic
2985 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_all
2986 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_all
2987 REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
2988 REAL(kind=dp), DIMENSION(3) :: my_rcc
2989 REAL(kind=dp), DIMENSION(3, 3) :: hmat
2990 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: eigenvals
2991 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: bc, xkp
2992 REAL(kind=dp), DIMENSION(:, :, :, :), &
2993 ALLOCATABLE, OPTIONAL :: berry_c
2994 REAL(kind=dp), DIMENSION(:, :, :, :), &
2995 ALLOCATABLE :: d_rs, h_rs, s_rs
2996 TYPE(cell_type), POINTER :: cell
2997 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_rs_img, matrix_ks_kp, &
2998 matrix_s_kp
2999 TYPE(dft_control_type), POINTER :: dft_control
3000 TYPE(kpoint_type), POINTER :: kpoints_all, kpoints_scf
3001 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3002 TYPE(mp_para_env_type), POINTER :: para_env
3003 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3004 POINTER :: sab_all
3005
3006 CALL timeset(routinen, handle)
3007 calc_bc = PRESENT(berry_c)
3008 my_do_parallel = .false.
3009 my_rcc = 0.0_dp
3010 IF (PRESENT(do_parallel)) my_do_parallel = do_parallel
3011 IF (PRESENT(rcc)) my_rcc = rcc
3012
3013 CALL get_qs_env(qs_env, &
3014 matrix_ks_kp=matrix_ks_kp, &
3015 matrix_s_kp=matrix_s_kp, &
3016 sab_all=sab_all, &
3017 cell=cell, &
3018 kpoints=kpoints_scf, &
3019 para_env=para_env, &
3020 dft_control=dft_control, &
3021 mos=mos)
3022
3023 CALL get_mo_set(mo_set=mos(1), nao=nao)
3024 CALL get_cell(cell=cell, h=hmat, periodic=periodic)
3025 nspin = SIZE(matrix_ks_kp, 1)
3026 nkp = SIZE(xkp, 2)
3027
3028 ! create kpoint environment kpoints_all which contains all neighbor cells R
3029 ! without considering any lattice symmetry
3030 NULLIFY (kpoints_all)
3031 CALL kpoint_create(kpoints_all)
3032 CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, n_img_scf)
3033
3034 CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index_all, &
3035 index_to_cell=index_to_cell_all)
3036 n_img_all = SIZE(index_to_cell_all, 2)
3037
3038 NULLIFY (moments_rs_img)
3039 CALL dbcsr_allocate_matrix_set(moments_rs_img, 3, n_img_all)
3040 ! D_μ,ν = <φ_μ|r|φ_ν>
3041 CALL build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc=my_rcc)
3042
3043 ALLOCATE (s_rs(1, nao, nao, n_img_all), h_rs(nspin, nao, nao, n_img_all), source=0.0_dp)
3044 ALLOCATE (d_rs(3, nao, nao, n_img_all), source=0.0_dp)
3045
3046 ! Convert real-space dbcsr matrices into arrays
3047 CALL replicate_rs_matrices(matrix_s_kp, kpoints_scf, s_rs, cell_to_index_all)
3048 CALL replicate_rs_matrices(matrix_ks_kp, kpoints_scf, h_rs, cell_to_index_all)
3049 CALL replicate_rs_matrices(moments_rs_img, kpoints_all, d_rs, cell_to_index_all)
3050
3051 mepos = 0
3052 num_pe = 1
3053 num_copy = nkp
3054 IF (my_do_parallel) THEN
3055 mepos = para_env%mepos
3056 num_pe = para_env%num_pe
3057 num_copy = ceiling(real(nkp)/num_pe)
3058 END IF
3059
3060 ALLOCATE (dipole(nspin, num_copy, 3, nao, nao), source=z_zero)
3061 IF (calc_bc) ALLOCATE (berry_c(nspin, num_copy, 3, nao), source=0.0_dp)
3062
3063!$OMP PARALLEL DEFAULT(NONE) PRIVATE(ikp, S_k, H_k, eigenvals, C_k, ispin, n, m, &
3064!$OMP i_dir, dS_dk_i, dH_dk_i, D_k, dip, bc, C_dS_C, C_dH_C, CDC, tmp_max, phase) &
3065!$OMP SHARED(num_pe, mepos, dipole, berry_c, nao, nspin, periodic, &
3066!$OMP nkp, xkp, S_rs, H_rs, D_rs, index_to_cell_all, hmat, calc_bc)
3067 ALLOCATE (ds_dk_i(nao, nao), c_ds_c(nao, nao), dh_dk_i(nao, nao), c_dh_c(nao, nao), source=z_zero)
3068 ALLOCATE (cdc(nao, nao), dip(3, nao, nao), s_k(nao, nao), h_k(nao, nao), source=z_zero)
3069 ALLOCATE (c_k(nao, nao), d_k(nao, nao), source=z_zero)
3070 ALLOCATE (eigenvals(nao), source=0.0_dp)
3071 IF (calc_bc) ALLOCATE (bc(3, nao), source=0.0_dp)
3072!$OMP DO COLLAPSE(2)
3073 DO ispin = 1, nspin
3074 DO ikp = 1, nkp
3075 IF (mod(ikp - 1, num_pe) /= mepos) cycle
3076
3077 ! S^R -> S(k), H^R -> H(k)
3078 s_k = 0
3079 h_k = 0
3080 CALL rs_to_kp(s_rs(1, :, :, :), s_k, index_to_cell_all, xkp(:, ikp))
3081 CALL rs_to_kp(h_rs(ispin, :, :, :), h_k, index_to_cell_all, xkp(:, ikp))
3082
3083 ! Diagonalize H(k)C(k) = S(k)C(k)ε(k)
3084 CALL geeig_right(h_k, s_k, eigenvals, c_k)
3085
3086 ! To have a smooth complex phase of C(k) as function of k, for every n, we force
3087 ! the largest C_μ,n(k) to be real.
3088 ! This is important to have a continuous dipole moment d_nm(k) as a function of k
3089 DO n = 1, nao
3090 tmp_max = c_k(1, n)
3091 DO mu = 1, nao
3092 IF (abs(c_k(mu, n)) < abs(tmp_max)) cycle
3093 tmp_max = c_k(mu, n)
3094 END DO
3095 phase = tmp_max/abs(tmp_max)
3096 c_k(:, n) = c_k(:, n)/phase
3097 END DO
3098
3099 DO i_dir = 1, 3 ! d^x, d^y, d^z
3100
3101 IF (periodic(i_dir) == 0) cycle
3102 ! ∇ S(k) = Σ_R iR S^R e^(ikR), ∇ H(k) = Σ_R iR H^R e^(ikR)
3103 CALL rs_to_kp(s_rs(1, :, :, :), ds_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
3104 CALL rs_to_kp(h_rs(ispin, :, :, :), dh_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
3105
3106 ! Σ_R D^R e^(ikR) = D(k), D_μ,ν = <φ_μ|r|φ_ν>
3107 CALL rs_to_kp(d_rs(i_dir, :, :, :), d_k(:, :), index_to_cell_all, xkp(:, ikp))
3108
3109 ! Basis transform to Kohn-Sham basis: (C^H) ∇ S C, (C^H) ∇ H C, (C^H) D C
3110 CALL gemm_square(c_k, 'C', ds_dk_i, 'N', c_k, 'N', c_ds_c)
3111 CALL gemm_square(c_k, 'C', dh_dk_i, 'N', c_k, 'N', c_dh_c)
3112 CALL gemm_square(c_k, 'C', d_k, 'N', c_k, 'N', cdc)
3113
3114 ! Compute the dipole
3115 ! d_nm (k) = - i/(ε(n)-ε(m)) [ (C^H)(dH(k)/dk)C ]_nm
3116 ! + i ε(n)/(ε(n)-ε(m)) [ (C^H)(dS(k)/dk)C ]_nm + [ (C^H)D(k)C ]_nm
3117 DO n = 1, nao
3118 DO m = 1, nao
3119 IF (n == m) cycle ! diagonal elements would need to be computed from
3120 ! a numerical k-derivative which is not implemented
3121 dip(i_dir, n, m) = -gaussi*c_dh_c(n, m)/(eigenvals(n) - eigenvals(m)) &
3122 + gaussi*eigenvals(n)*c_ds_c(n, m)/(eigenvals(n) - eigenvals(m)) &
3123 + cdc(n, m)
3124 END DO
3125 END DO
3126 END DO
3127 ! Compute the Berry curvature from the dipoles
3128 ! Ω^γ_n = Σ_m 2*Im[d^α_nm d^β_mn], where, α, β, γ belong to {x, y, z}
3129 IF (calc_bc) THEN
3130 bc = 0.0_dp
3131 DO i_dir = 1, 3
3132 DO n = 1, nao
3133 DO m = 1, nao
3134 IF (n == m) cycle
3135 bc(i_dir, n) = bc(i_dir, n) &
3136 + 2*aimag(dip(1 + mod(i_dir, 3), n, m)*dip(1 + mod(i_dir + 1, 3), m, n))
3137 END DO
3138 END DO
3139 END DO
3140 END IF
3141 ! Store the dipoles and berry curvature for each MPI rank
3142 dipole(ispin, ceiling(real(ikp)/num_pe), :, :, :) = dip(:, :, :)
3143 IF (calc_bc) berry_c(ispin, ceiling(real(ikp)/num_pe), :, :) = bc(:, :)
3144 END DO
3145 END DO
3146!$OMP END DO
3147 DEALLOCATE (ds_dk_i, c_ds_c, dh_dk_i, c_dh_c, cdc, dip, s_k, h_k, c_k, d_k, eigenvals)
3148 IF (calc_bc) DEALLOCATE (bc)
3149!$OMP END PARALLEL
3150 DEALLOCATE (s_rs, h_rs, d_rs)
3151 CALL dbcsr_deallocate_matrix_set(moments_rs_img)
3152 CALL kpoint_release(kpoints_all)
3153 CALL timestop(handle)
3154 END SUBROUTINE qs_moment_kpoints_deep
3155
3156! **************************************************************************************************
3157!> \brief Calculates interband k-point dipoles in the existing SCF MO basis.
3158!> \param qs_env ...
3159!> \param dipole ...
3160!> \param rcc retained for interface compatibility; interband dipoles are origin independent
3161!> \param nmo_spin_out number of SCF MOs available for each spin
3162! **************************************************************************************************
3163 SUBROUTINE qs_moment_kpoints_scf_mos(qs_env, dipole, rcc, nmo_spin_out)
3164 TYPE(qs_environment_type), POINTER :: qs_env
3165 COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
3166 ALLOCATABLE :: dipole
3167 REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
3168 INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT), &
3169 OPTIONAL :: nmo_spin_out
3170
3171 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_kpoints_scf_mos'
3172
3173 INTEGER :: handle, i_dir, ikp, ikp_local, ispin, &
3174 m, n, nao, nkp, nmo, nspin
3175 INTEGER, DIMENSION(:), ALLOCATABLE :: nmo_spin
3176 INTEGER, DIMENSION(2) :: kp_range
3177 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3178 LOGICAL :: my_kpgrp
3179 REAL(kind=dp), PARAMETER :: eps_degenerate = 1.0e-10_dp
3180 REAL(kind=dp) :: cimag, creal, energy_diff
3181 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: eigenvalues_kp
3182 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvals
3183 TYPE(cp_blacs_env_type), POINTER :: blacs_env_all
3184 TYPE(cp_fm_struct_type), POINTER :: moment_struct
3185 TYPE(cp_fm_struct_type), POINTER :: fm_struct
3186 TYPE(cp_fm_type) :: fm_dummy, fm_tmp, mo_coeff_im_global, &
3187 mo_coeff_re_global, moment_im, &
3188 moment_re
3189 TYPE(cp_fm_type), POINTER :: mo_coeff_im, mo_coeff_re
3190 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: overlap_deriv
3191 TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix
3192 TYPE(dft_control_type), POINTER :: dft_control
3193 TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env
3194 TYPE(kpoint_env_type), POINTER :: kp
3195 TYPE(kpoint_type), POINTER :: kpoints_scf
3196 TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp
3197 TYPE(mp_para_env_type), POINTER :: para_env, para_env_kp
3198 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3199 POINTER :: sab_kp, sab_orb
3200 TYPE(qs_ks_env_type), POINTER :: ks_env
3201
3202 CALL timeset(routinen, handle)
3203
3204 NULLIFY (blacs_env_all, cell_to_index, cmatrix, dft_control, eigenvals, fm_struct, kp, &
3205 kp_env, kpoints_scf, ks_env, mo_coeff_im, mo_coeff_re, moment_struct, mos_kp, &
3206 overlap_deriv, para_env, para_env_kp, rmatrix, sab_kp, sab_orb)
3207 IF (PRESENT(rcc)) THEN
3208 mark_used(rcc)
3209 END IF
3210
3211 CALL get_qs_env(qs_env, dft_control=dft_control, kpoints=kpoints_scf, ks_env=ks_env, &
3212 para_env=para_env, sab_orb=sab_orb)
3213 cpassert(ASSOCIATED(dft_control))
3214 cpassert(ASSOCIATED(kpoints_scf))
3215 cpassert(ASSOCIATED(ks_env))
3216 cpassert(ASSOCIATED(para_env))
3217 cpassert(ASSOCIATED(sab_orb))
3218
3219 CALL get_kpoint_info(kpoints_scf, nkp=nkp, kp_range=kp_range, kp_env=kp_env, &
3220 para_env_kp=para_env_kp, blacs_env_all=blacs_env_all, &
3221 cell_to_index=cell_to_index, sab_nl=sab_kp)
3222 IF (kp_range(2) >= kp_range(1)) THEN
3223 cpassert(ASSOCIATED(kp_env))
3224 END IF
3225 cpassert(ASSOCIATED(para_env_kp))
3226 cpassert(ASSOCIATED(blacs_env_all))
3227 cpassert(ASSOCIATED(cell_to_index))
3228 cpassert(ASSOCIATED(sab_kp))
3229
3230 CALL build_overlap_matrix(ks_env, matrixkp_s=overlap_deriv, nderivative=1, &
3231 basis_type_a="ORB", basis_type_b="ORB", sab_nl=sab_orb, &
3232 ext_kpoints=kpoints_scf)
3233
3234 nspin = dft_control%nspins
3235 CALL dbcsr_get_info(overlap_deriv(1, 1)%matrix, nfullrows_total=nao)
3236 ALLOCATE (nmo_spin(nspin), source=0)
3237 IF (kp_range(2) >= kp_range(1)) THEN
3238 kp => kp_env(1)%kpoint_env
3239 mos_kp => kp%mos
3240 cpassert(ASSOCIATED(mos_kp))
3241 DO ispin = 1, nspin
3242 CALL get_mo_set(mos_kp(1, ispin), nmo=nmo_spin(ispin))
3243 END DO
3244 END IF
3245 CALL para_env%max(nmo_spin)
3246 ALLOCATE (dipole(nspin, nkp, 3, maxval(nmo_spin), maxval(nmo_spin)), source=z_zero)
3247 IF (PRESENT(nmo_spin_out)) THEN
3248 ALLOCATE (nmo_spin_out(nspin))
3249 nmo_spin_out(:) = nmo_spin(:)
3250 END IF
3251
3252 ALLOCATE (rmatrix, cmatrix)
3253 CALL dbcsr_create(rmatrix, template=overlap_deriv(1, 1)%matrix, &
3254 matrix_type=dbcsr_type_antisymmetric)
3255 CALL dbcsr_create(cmatrix, template=overlap_deriv(1, 1)%matrix, &
3256 matrix_type=dbcsr_type_symmetric)
3257 CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_kp)
3258 CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_kp)
3259
3260 DO ikp = 1, nkp
3261 my_kpgrp = (ikp >= kp_range(1) .AND. ikp <= kp_range(2))
3262 IF (my_kpgrp) THEN
3263 ikp_local = ikp - kp_range(1) + 1
3264 kp => kp_env(ikp_local)%kpoint_env
3265 mos_kp => kp%mos
3266 ELSE
3267 NULLIFY (kp, mos_kp)
3268 END IF
3269 DO ispin = 1, nspin
3270 nmo = nmo_spin(ispin)
3271 ALLOCATE (eigenvalues_kp(nmo), source=0.0_dp)
3272
3273 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nmo, &
3274 para_env=para_env, context=blacs_env_all)
3275 CALL cp_fm_create(mo_coeff_re_global, fm_struct)
3276 CALL cp_fm_create(mo_coeff_im_global, fm_struct)
3277 CALL cp_fm_create(fm_tmp, fm_struct)
3278 CALL cp_fm_struct_release(fm_struct)
3279 CALL cp_fm_struct_create(moment_struct, nrow_global=nmo, ncol_global=nmo, &
3280 para_env=para_env, context=blacs_env_all)
3281 CALL cp_fm_create(moment_re, moment_struct)
3282 CALL cp_fm_create(moment_im, moment_struct)
3283 CALL cp_fm_struct_release(moment_struct)
3284
3285 IF (my_kpgrp) THEN
3286 CALL get_mo_set(mos_kp(1, ispin), eigenvalues=eigenvals, mo_coeff=mo_coeff_re)
3287 CALL get_mo_set(mos_kp(2, ispin), mo_coeff=mo_coeff_im)
3288 cpassert(ASSOCIATED(eigenvals))
3289 cpassert(ASSOCIATED(mo_coeff_re))
3290 cpassert(ASSOCIATED(mo_coeff_im))
3291 IF (para_env_kp%is_source()) eigenvalues_kp(1:nmo) = eigenvals(1:nmo)
3292 CALL cp_fm_copy_general(mo_coeff_re, mo_coeff_re_global, para_env)
3293 CALL cp_fm_copy_general(mo_coeff_im, mo_coeff_im_global, para_env)
3294 ELSE
3295 CALL cp_fm_copy_general(fm_dummy, mo_coeff_re_global, para_env)
3296 CALL cp_fm_copy_general(fm_dummy, mo_coeff_im_global, para_env)
3297 END IF
3298 CALL para_env%sum(eigenvalues_kp)
3299
3300 DO i_dir = 1, 3
3301 CALL dbcsr_set(rmatrix, 0.0_dp)
3302 CALL dbcsr_set(cmatrix, 0.0_dp)
3303 CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=overlap_deriv, &
3304 ispin=i_dir + 1, xkp=kpoints_scf%xkp(:, ikp), &
3305 cell_to_index=cell_to_index, sab_nl=sab_kp)
3306
3307 ! Project the complex AO derivative operator as C^H A C. The
3308 ! off-diagonal length-gauge dipoles follow from the energy-gap relation.
3309 CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_re_global, fm_tmp, nmo)
3310 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3311 1.0_dp, mo_coeff_re_global, fm_tmp, 0.0_dp, moment_re)
3312 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3313 -1.0_dp, mo_coeff_im_global, fm_tmp, 0.0_dp, moment_im)
3314
3315 CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_im_global, fm_tmp, nmo)
3316 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3317 1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
3318 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3319 1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)
3320
3321 CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_re_global, fm_tmp, nmo)
3322 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3323 1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
3324 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3325 1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)
3326
3327 CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_im_global, fm_tmp, nmo)
3328 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3329 -1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_re)
3330 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3331 1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_im)
3332
3333 DO n = 1, nmo
3334 DO m = 1, nmo
3335 IF (n == m) cycle
3336 energy_diff = eigenvalues_kp(m) - eigenvalues_kp(n)
3337 IF (abs(energy_diff) <= eps_degenerate) cycle
3338 CALL cp_fm_get_element(moment_re, m, n, creal)
3339 CALL cp_fm_get_element(moment_im, m, n, cimag)
3340 IF (para_env%is_source()) then
3341 dipole(ispin, ikp, i_dir, n, m) = cmplx(creal, cimag, kind=dp)/energy_diff
3342 end if
3343 END DO
3344 END DO
3345 END DO
3346 CALL cp_fm_release(mo_coeff_im_global)
3347 CALL cp_fm_release(mo_coeff_re_global)
3348 CALL cp_fm_release(moment_im)
3349 CALL cp_fm_release(moment_re)
3350 CALL cp_fm_release(fm_tmp)
3351 DEALLOCATE (eigenvalues_kp)
3352 END DO
3353 END DO
3354
3355 DO ispin = 1, nspin
3356 DO ikp = 1, nkp
3357 DO i_dir = 1, 3
3358 CALL para_env%sum(dipole(ispin, ikp, i_dir, :, :))
3359 END DO
3360 END DO
3361 END DO
3362
3363 CALL dbcsr_deallocate_matrix(cmatrix)
3364 CALL dbcsr_deallocate_matrix(rmatrix)
3365 CALL dbcsr_deallocate_matrix_set(overlap_deriv)
3366 DEALLOCATE (nmo_spin)
3367 CALL timestop(handle)
3368
3369 END SUBROUTINE qs_moment_kpoints_scf_mos
3370
3371! **************************************************************************************************
3372!> \brief Calculate and print dipole moment elements d_nm(k) for k-point calculations
3373!> \param qs_env ...
3374!> \param nmoments ...
3375!> \param reference ...
3376!> \param ref_point ...
3377!> \param max_nmo ...
3378!> \param unit_number ...
3379!> \author Shridhar Shanbhag
3380! **************************************************************************************************
3381 SUBROUTINE qs_moment_kpoints(qs_env, nmoments, reference, ref_point, max_nmo, unit_number)
3382 TYPE(qs_environment_type), POINTER :: qs_env
3383 INTEGER, INTENT(IN) :: nmoments, reference, max_nmo
3384 REAL(dp), DIMENSION(:), INTENT(IN), POINTER :: ref_point
3385 INTEGER, INTENT(IN) :: unit_number
3386 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_kpoints'
3387 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
3388 COMPLEX(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dipole_to_print
3389 COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
3390 ALLOCATABLE :: dipole
3391 INTEGER :: handle, i_dir, ikp, nmo_dim, nkp, nao, &
3392 num_pe, mepos, n, m, &
3393 ispin, nspin, nmin, nmax, homo
3394 INTEGER, DIMENSION(:), ALLOCATABLE :: nmo_spin_scf
3395 LOGICAL :: explicit_kpnts, explicit_kpset, use_scf_mos
3396 REAL(kind=dp), DIMENSION(3) :: rcc
3397 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: xkp
3398 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: bc_to_print
3399 REAL(kind=dp), DIMENSION(:, :, :, :), ALLOCATABLE :: berry_c
3400 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3401 TYPE(mp_para_env_type), POINTER :: para_env
3402 TYPE(section_vals_type), POINTER :: kpnts, kpset
3403 CHARACTER(LEN=default_string_length), &
3404 DIMENSION(:), POINTER :: special_pnts
3405
3406 CALL timeset(routinen, handle)
3407
3408 IF (nmoments > 1) cpabort("KPOINT quadrupole and higher moments not implemented.")
3409 IF (max_nmo < 0) cpabort("Negative maximum number of molecular orbitals max_nmo provided.")
3410
3411 CALL get_qs_env(qs_env, &
3412 para_env=para_env, &
3413 matrix_ks_kp=matrix_ks_kp, &
3414 mos=mos)
3415
3416 CALL get_mo_set(mo_set=mos(1), nao=nao)
3417 CALL get_xkp_for_dipole_calc(qs_env, xkp, special_pnts)
3418 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
3419 nspin = SIZE(matrix_ks_kp, 1)
3420 nkp = SIZE(xkp, 2)
3421
3422 kpset => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINT_SET")
3423 kpnts => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINTS")
3424 CALL section_vals_get(kpset, explicit=explicit_kpset)
3425 CALL section_vals_get(kpnts, explicit=explicit_kpnts)
3426 use_scf_mos = .NOT. explicit_kpset .AND. .NOT. explicit_kpnts
3427
3428 IF (unit_number > 0) WRITE (unit_number, fmt="(/,T2,A)") &
3429 '!-----------------------------------------------------------------------------!'
3430 IF (unit_number > 0) WRITE (unit_number, "(T22,A)") "Periodic Dipole Matrix Elements"
3431
3432 IF (use_scf_mos) THEN
3433 CALL qs_moment_kpoints_scf_mos(qs_env, dipole, rcc, nmo_spin_scf)
3434 nmo_dim = SIZE(dipole, 4)
3435 ALLOCATE (berry_c(nspin, nkp, 3, nmo_dim), source=0.0_dp)
3436 DO ispin = 1, nspin
3437 DO ikp = 1, nkp
3438 DO i_dir = 1, 3
3439 DO n = 1, nmo_dim
3440 DO m = 1, nmo_dim
3441 IF (n == m) cycle
3442 berry_c(ispin, ikp, i_dir, n) = berry_c(ispin, ikp, i_dir, n) &
3443 + 2*aimag(dipole(ispin, ikp, 1 + mod(i_dir, 3), n, m)* &
3444 dipole(ispin, ikp, 1 + mod(i_dir + 1, 3), m, n))
3445 END DO
3446 END DO
3447 END DO
3448 END DO
3449 END DO
3450 ELSE
3451 CALL qs_moment_kpoints_deep(qs_env, &
3452 xkp, &
3453 dipole, &
3454 rcc, &
3455 berry_c, &
3456 do_parallel=.true.)
3457 nmo_dim = nao
3458 END IF
3459
3460 ALLOCATE (dipole_to_print(3, nmo_dim, nmo_dim), source=z_zero)
3461 ALLOCATE (bc_to_print(3, nmo_dim), source=0.0_dp)
3462
3463 mepos = para_env%mepos
3464 num_pe = para_env%num_pe
3465
3466 DO ikp = 1, nkp
3467 DO ispin = 1, nspin
3468 CALL get_mo_set(mo_set=mos(ispin), homo=homo)
3469 nmin = max(1, homo - (max_nmo - 1)/2)
3470 nmax = min(nao, homo + max_nmo/2)
3471 IF (max_nmo == 0) THEN
3472 nmin = 1
3473 nmax = nao
3474 END IF
3475 IF (use_scf_mos) THEN
3476 nmax = min(nmax, nmo_spin_scf(ispin))
3477 END IF
3478 dipole_to_print = 0.0_dp
3479 bc_to_print = 0.0_dp
3480 IF (use_scf_mos) THEN
3481 dipole_to_print(:, :, :) = dipole(ispin, ikp, :, :, :)
3482 bc_to_print(:, :) = berry_c(ispin, ikp, :, :)
3483 ELSE IF (mod(ikp - 1, num_pe) == mepos) THEN
3484 dipole_to_print(:, :, :) = dipole(ispin, ceiling(real(ikp)/num_pe), :, :, :)
3485 bc_to_print(:, :) = berry_c(ispin, ceiling(real(ikp)/num_pe), :, :)
3486 END IF
3487 IF (.NOT. use_scf_mos) THEN
3488 CALL para_env%sum(dipole_to_print)
3489 CALL para_env%sum(bc_to_print)
3490 END IF
3491 IF (unit_number > 0) THEN
3492 IF (special_pnts(ikp) /= "") WRITE (unit_number, "(/,2X,A,A)") &
3493 "Special point: ", adjustl(trim(special_pnts(ikp)))
3494 WRITE (unit_number, "(/,1X,A,I3,1X,3(A,1F12.6))") &
3495 "Kpoint:", ikp, ", kx:", xkp(1, ikp), ", ky:", xkp(2, ikp), ", kz:", xkp(3, ikp)
3496 IF (nspin > 1) WRITE (unit_number, "(/,2X,A,I2)") "Open Shell System. Spin:", ispin
3497 WRITE (unit_number, "(2X,A)") " kp n m Re(dx_nm) Im(dx_nm) &
3498 & Re(dy_nm) Im(dy_nm) Re(dz_nm) Im(dz_nm)"
3499 DO n = nmin, nmax
3500 DO m = nmin, nmax
3501 IF (n == m) cycle
3502 WRITE (unit_number, "(2X,I4,2I4,6(G11.3))") ikp, n, m, dipole_to_print(1:3, n, m)
3503 END DO
3504 END DO
3505 WRITE (unit_number, "(/,1X,A)") "Berry Curvature"
3506 WRITE (unit_number, "(2X,A)") " kp n YZ ZX XY"
3507 DO n = nmin, nmax
3508 WRITE (unit_number, "(2X,2I5,3(1X,G11.3))") &
3509 ikp, n, bc_to_print(1, n), bc_to_print(2, n), bc_to_print(3, n)
3510 END DO
3511 END IF
3512 END DO
3513 END DO
3514 DEALLOCATE (dipole_to_print, bc_to_print, berry_c, dipole)
3515 IF (ALLOCATED(nmo_spin_scf)) DEALLOCATE (nmo_spin_scf)
3516 DEALLOCATE (special_pnts, xkp)
3517
3518 CALL timestop(handle)
3519
3520 END SUBROUTINE qs_moment_kpoints
3521
3522END 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:81
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:104
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
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)
...
Definition ai_angmom.F:52
ai_moments
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition ai_moments.F:17
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:1348
ai_moments::contract_cossin
subroutine, public contract_cossin(cos_block, sin_block, iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, cosab, sinab, ldab, work, ldwork)
...
Definition ai_moments.F:77
ai_moments::cossin
subroutine, public cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max, npgfb, zetb, rpgfb, lb_min, rac, rbc, kvec, cosab, sinab, dcosab, dsinab)
...
Definition ai_moments.F:155
ai_moments::moment
subroutine, public moment(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lc_max, rac, rbc, mab)
...
Definition ai_moments.F:802
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
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:138
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:233
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:77
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
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:125
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:631
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:412
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:367
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:548
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:1524
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:1038
cp_fm_types::cp_fm_get_element
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
Definition cp_fm_types.F:627
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:534
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:165
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:1259
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
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_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 transformations from k-space to R-space for Fortran array matrices.
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:230
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:62
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:1175
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:1041
kpoint_methods::kpoint_initialize
subroutine, public kpoint_initialize(kpoint, particle_set, cell)
Generate the kpoints and initialize the kpoint environment.
Definition kpoint_methods.F:156
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:302
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:544
kpoint_types::kpoint_create
subroutine, public kpoint_create(kpoint)
Create a kpoint environment.
Definition kpoint_types.F:250
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:751
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:142
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:143
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:2178
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
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
physcon::debye
real(kind=dp), parameter, public debye
Definition physcon.F:201
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:542
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:984
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition qs_ks_types.F:343
qs_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)
Get the components of a MO set data structure.
Definition qs_mo_types.F:401
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>.
Definition qs_moments.F:14
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:160
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:839
qs_moments::print_moments_nl
subroutine, public print_moments_nl(unit_number, nmom, rlab, mmom, rmom_vel)
...
Definition qs_moments.F:2347
qs_moments::set_label
subroutine, public set_label(label, ix, iy, iz)
...
Definition qs_moments.F:2208
qs_moments::qs_moment_locop
subroutine, public qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
...
Definition qs_moments.F:1807
qs_moments::print_moments
subroutine, public print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmom, rmom_vel)
...
Definition qs_moments.F:2239
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:1045
qs_moments::op_orbbas
subroutine, public op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
...
Definition qs_moments.F:1685
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:2791
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:2415
qs_moments::qs_moment_kpoints_deep
subroutine, public qs_moment_kpoints_deep(qs_env, xkp, dipole, rcc, berry_c, do_parallel)
Calculates the dipole moments and berry curvature for periodic systems for kpoints.
Definition qs_moments.F:2970
qs_moments::qs_moment_kpoints
subroutine, public qs_moment_kpoints(qs_env, nmoments, reference, ref_point, max_nmo, unit_number)
Calculate and print dipole moment elements d_nm(k) for k-point calculations.
Definition qs_moments.F:3382
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:3164
qs_moments::build_berry_kpoint_matrix
subroutine, public build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
...
Definition qs_moments.F:1081
qs_moments::op_orbbas_rtp
subroutine, public op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
...
Definition qs_moments.F:1715
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:433
qs_moments::qs_moment_berry_phase
subroutine, public qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number)
...
Definition qs_moments.F:1288
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:275
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
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:69
cp_control_types::dft_control_type
Definition cp_control_types.F:758
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:114
cp_result_types::cp_result_type
contains arbitrary information which need to be stored
Definition cp_result_types.F:73
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:108
kpoint_types::kpoint_env_type
Keeps information about a specific k-point.
Definition kpoint_types.F:90
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:187
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:767
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:224
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:47
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