Line data Source code
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 : ! **************************************************************************************************
14 : MODULE 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 :
143 : CONTAINS
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 3428 : 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 3428 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
180 3428 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
181 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
182 3428 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
183 : TYPE(cell_type), POINTER :: cell
184 3428 : 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 3428 : DIMENSION(:), POINTER :: nl_iterator
188 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
189 3428 : POINTER :: sab_orb
190 3428 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
191 3428 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
192 : TYPE(qs_kind_type), POINTER :: qs_kind
193 :
194 3428 : IF (nmoments < 1) RETURN
195 :
196 3428 : CALL timeset(routineN, handle)
197 :
198 3428 : my_all_images = .FALSE.
199 3428 : IF (PRESENT(all_images)) my_all_images = all_images
200 3428 : my_minimum_image = .FALSE.
201 3428 : IF (PRESENT(minimum_image)) my_minimum_image = minimum_image
202 3428 : my_neighbor_image = .FALSE.
203 3428 : IF (PRESENT(neighbor_image)) my_neighbor_image = neighbor_image
204 6 : CPASSERT(.NOT. (my_all_images .AND. my_neighbor_image))
205 :
206 3428 : nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
207 3428 : CPASSERT(SIZE(moments) >= nm)
208 3428 : component_start = 1
209 3428 : IF (PRESENT(first_component)) component_start = first_component
210 3428 : CPASSERT(component_start >= 1 .AND. component_start <= nm)
211 :
212 3428 : 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 3428 : sab_orb=sab_orb)
217 :
218 3428 : IF (my_minimum_image) THEN
219 150 : half_cell = 0.5_dp*NORM2(cell%hmat, DIM=1)
220 : END IF
221 :
222 3428 : 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 6750 : basis_type=basis_type)
227 :
228 17140 : ALLOCATE (mab(maxco, maxco, nm))
229 3428 : mab(:, :, :) = 0.0_dp
230 :
231 13712 : ALLOCATE (work(maxco, maxsgf))
232 3428 : work(:, :) = 0.0_dp
233 :
234 20738 : ALLOCATE (mint(nm))
235 13882 : DO i = component_start, nm
236 13882 : NULLIFY (mint(i)%block)
237 : END DO
238 :
239 17166 : ALLOCATE (basis_set_list(nkind))
240 10310 : DO ikind = 1, nkind
241 6882 : qs_kind => qs_kind_set(ikind)
242 6882 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
243 10310 : IF (ASSOCIATED(basis_set_a)) THEN
244 6882 : basis_set_list(ikind)%gto_basis_set => basis_set_a
245 : ELSE
246 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
247 : END IF
248 : END DO
249 3428 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
250 39522 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
251 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
252 36094 : iatom=iatom, jatom=jatom, r=rab)
253 36094 : basis_set_a => basis_set_list(ikind)%gto_basis_set
254 36094 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
255 36094 : basis_set_b => basis_set_list(jkind)%gto_basis_set
256 36094 : 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 36094 : nseta = basis_set_a%nset
280 36094 : nsetb = basis_set_b%nset
281 :
282 9506 : IF (inode == 1) last_jatom = 0
283 :
284 36094 : IF (my_minimum_image) THEN
285 482 : IF (ANY(ABS(rab(:)) > half_cell(:))) CYCLE
286 : END IF
287 :
288 36092 : new_atom_pair = jatom /= last_jatom
289 36092 : IF (.NOT. my_all_images .AND. .NOT. new_atom_pair) THEN
290 : CYCLE
291 : END IF
292 :
293 13612 : IF (new_atom_pair) THEN
294 13363 : last_jatom = jatom
295 :
296 13363 : IF (iatom <= jatom) THEN
297 8323 : irow = iatom
298 8323 : icol = jatom
299 : ELSE
300 5040 : irow = jatom
301 5040 : icol = iatom
302 : END IF
303 :
304 54322 : DO i = component_start, nm
305 40959 : NULLIFY (mint(i)%block)
306 : CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
307 40959 : row=irow, col=icol, BLOCK=mint(i)%block, found=found)
308 40959 : CPASSERT(found .AND. ASSOCIATED(mint(i)%block))
309 2009108 : mint(i)%block = 0._dp
310 : END DO
311 : END IF
312 :
313 : ! Reference point and image geometry.
314 13612 : IF (PRESENT(ref_points)) THEN
315 0 : rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
316 13612 : ELSE IF (PRESENT(ref_point)) THEN
317 51368 : rc(:) = ref_point(:)
318 : ELSE
319 770 : rc(:) = 0._dp
320 : END IF
321 :
322 13612 : IF (my_all_images) THEN
323 : ! Keep the image represented by the current neighbor-list node.
324 400 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
325 1600 : rb(:) = ra(:) + rab(:)
326 400 : rac(:) = pbc(rc, ra, cell)
327 400 : rbc(:) = pbc(rc, rb, cell)
328 1600 : dab = NORM2(rab)
329 13212 : ELSE IF (my_neighbor_image) THEN
330 : ! Keep the neighbor-list image while selecting one atom-pair representative.
331 72 : ra(:) = particle_set(iatom)%r(:)
332 72 : rb(:) = ra(:) + rab(:)
333 18 : rac(:) = pbc(rc, ra, cell)
334 72 : rbc(:) = rac(:) + rab(:)
335 72 : dab = NORM2(rab)
336 : ELSE
337 : ! Fold atomic positions around the reference point.
338 105552 : ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
339 105552 : rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
340 52776 : rab(:) = ra(:) - rb(:)
341 52776 : rac(:) = ra(:) - rc(:)
342 52776 : rbc(:) = rb(:) - rc(:)
343 52776 : dab = NORM2(rab)
344 : END IF
345 :
346 70760 : DO iset = 1, nseta
347 :
348 21054 : ncoa = npgfa(iset)*ncoset(la_max(iset))
349 21054 : sgfa = first_sgfa(1, iset)
350 :
351 73709 : DO jset = 1, nsetb
352 :
353 39043 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
354 :
355 38161 : ncob = npgfb(jset)*ncoset(lb_max(jset))
356 38161 : 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 38161 : rpgfb(:, jset), nmoments, rac, rbc, mab)
363 :
364 : ! Contraction step
365 179128 : 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 119913 : 0.0_dp, work(1, 1), SIZE(work, 1))
371 :
372 158956 : 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 78094 : 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 41819 : 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 3428 : CALL neighbor_list_iterator_release(nl_iterator)
398 :
399 : ! Release work storage
400 3428 : DEALLOCATE (mab, basis_set_list)
401 3428 : DEALLOCATE (work)
402 13882 : DO i = component_start, nm
403 13882 : NULLIFY (mint(i)%block)
404 : END DO
405 3428 : DEALLOCATE (mint)
406 :
407 3428 : CALL timestop(handle)
408 :
409 6856 : 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 4302 : SUBROUTINE build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, &
431 4302 : 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 4302 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
456 4302 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
457 4302 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
458 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
459 4302 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab_tmp
460 4302 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mom_block
461 4302 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mom_block_der
462 : TYPE(cell_type), POINTER :: cell
463 4302 : 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 4302 : DIMENSION(:), POINTER :: nl_iterator
467 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
468 4302 : POINTER :: sab_all, sab_orb
469 4302 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
470 4302 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
471 : TYPE(qs_kind_type), POINTER :: qs_kind
472 :
473 4302 : nmom_build = MAX(nmoments, nmoments_der) ! build moments up to order nmom_buiod
474 4302 : IF (nmom_build < 1) RETURN
475 :
476 4302 : CALL timeset(routineN, handle)
477 :
478 4302 : nders = 1 ! only first order derivatives
479 4302 : dimders = ncoset(nders) - 1
480 :
481 4302 : my_minimum_image = .FALSE.
482 4302 : IF (PRESENT(minimum_image)) my_minimum_image = minimum_image
483 4302 : my_ordered = .FALSE.
484 4302 : IF (PRESENT(ordered)) my_ordered = ordered
485 4302 : my_neighbor_image = my_ordered
486 4302 : IF (PRESENT(neighbor_image)) my_neighbor_image = neighbor_image
487 4302 : CPASSERT(.NOT. (PRESENT(deltaR) .AND. PRESENT(lambda)))
488 4302 : IF (PRESENT(lambda) .OR. PRESENT(deltaR)) THEN
489 510 : CPASSERT(my_ordered)
490 : END IF
491 4302 : accumulate_images = my_ordered .AND. my_neighbor_image
492 4302 : IF (PRESENT(deltaR) .OR. PRESENT(lambda)) accumulate_images = .FALSE.
493 4302 : 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 4302 : IF (my_ordered .AND. .NOT. PRESENT(deltaR) .AND. .NOT. PRESENT(lambda)) THEN
498 3750 : derivative_sign = -1._dp
499 : END IF
500 :
501 4302 : 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 4302 : sab_all=sab_all)
508 :
509 4302 : IF (my_ordered) THEN
510 4260 : CPASSERT(nmoments_der > 0)
511 4260 : CPASSERT(dbcsr_get_matrix_type(moments_der(1, 1)%matrix) == dbcsr_type_no_symmetry)
512 4260 : IF (nmoments > 0) THEN
513 3750 : CPASSERT(PRESENT(moments))
514 3750 : CPASSERT(dbcsr_get_matrix_type(moments(1)%matrix) == dbcsr_type_no_symmetry)
515 : END IF
516 : END IF
517 :
518 4302 : IF (my_minimum_image) THEN
519 0 : half_cell = 0.5_dp*NORM2(cell%hmat, DIM=1)
520 : END IF
521 :
522 4302 : 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 4848 : basis_type=basis_type)
528 :
529 4302 : IF (nmoments > 0) THEN
530 3790 : nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
531 3790 : CPASSERT(SIZE(moments) == nm)
532 : ! storage for integrals
533 18938 : ALLOCATE (mab(maxco, maxco, nm))
534 : ! blocks
535 3790 : mab(:, :, :) = 0.0_dp
536 45480 : ALLOCATE (mom_block(nm))
537 37900 : DO i = 1, nm
538 37900 : NULLIFY (mom_block(i)%block)
539 : END DO
540 : END IF
541 :
542 4302 : IF (nmoments_der > 0) THEN
543 4302 : M_dim = ncoset(nmoments_der) - 1
544 4302 : CPASSERT(SIZE(moments_der, dim=1) == M_dim)
545 4302 : CPASSERT(SIZE(moments_der, dim=2) == dimders)
546 : ! storage for integrals
547 25800 : ALLOCATE (difmab(maxco, maxco, M_dim, dimders))
548 4302 : difmab(:, :, :, :) = 0.0_dp
549 : ! blocks
550 136368 : ALLOCATE (mom_block_der(M_dim, dimders))
551 39720 : DO i = 1, M_dim
552 145974 : DO ider = 1, dimders
553 141672 : NULLIFY (mom_block_der(i, ider)%block)
554 : END DO
555 : END DO
556 : END IF
557 :
558 17202 : ALLOCATE (work(maxco, maxsgf))
559 4302 : work(:, :) = 0.0_dp
560 :
561 4302 : NULLIFY (basis_set_a, basis_set_b, basis_set_list)
562 4302 : NULLIFY (qs_kind)
563 20138 : ALLOCATE (basis_set_list(nkind))
564 11534 : DO ikind = 1, nkind
565 7232 : qs_kind => qs_kind_set(ikind)
566 7232 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
567 11534 : IF (ASSOCIATED(basis_set_a)) THEN
568 7232 : basis_set_list(ikind)%gto_basis_set => basis_set_a
569 : ELSE
570 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
571 : END IF
572 : END DO
573 :
574 : ! Calculate derivatives looping over neighbour list
575 4302 : NULLIFY (nl_iterator)
576 4302 : IF (my_ordered) THEN
577 4260 : CALL neighbor_list_iterator_create(nl_iterator, sab_all)
578 : ELSE
579 42 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
580 : END IF
581 308016 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
582 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
583 303714 : iatom=iatom, jatom=jatom, r=rab)
584 303714 : basis_set_a => basis_set_list(ikind)%gto_basis_set
585 303714 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
586 303714 : basis_set_b => basis_set_list(jkind)%gto_basis_set
587 303714 : 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 303714 : IF (SIZE(sphi_a, 1) == 0) CYCLE
611 303693 : IF (SIZE(sphi_b, 1) == 0) CYCLE
612 :
613 303693 : nseta = basis_set_a%nset
614 303693 : nsetb = basis_set_b%nset
615 :
616 303693 : IF (inode == 1) last_jatom = 0
617 :
618 303693 : IF (my_minimum_image) THEN
619 0 : IF (ANY(ABS(rab(:)) > half_cell(:))) CYCLE
620 : END IF
621 :
622 303693 : 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 303693 : IF (.NOT. accumulate_images .AND. .NOT. new_atom_pair) CYCLE
626 :
627 284811 : IF (new_atom_pair) THEN
628 21597 : last_jatom = jatom
629 :
630 : ! get blocks
631 21597 : IF (my_ordered) THEN
632 20253 : irow = iatom
633 20253 : icol = jatom
634 1344 : ELSE IF (iatom <= jatom) THEN
635 710 : irow = iatom
636 710 : icol = jatom
637 : ELSE
638 634 : irow = jatom
639 634 : icol = iatom
640 : END IF
641 :
642 21597 : IF (nmoments > 0) THEN
643 192960 : DO i = 1, nm
644 173664 : NULLIFY (mom_block(i)%block)
645 : ! get block from pre calculated overlap matrix
646 : CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
647 173664 : row=irow, col=icol, BLOCK=mom_block(i)%block, found=found)
648 173664 : CPASSERT(found .AND. ASSOCIATED(mom_block(i)%block))
649 13565520 : mom_block(i)%block = 0._dp
650 : END DO
651 : END IF
652 21597 : IF (nmoments_der > 0) THEN
653 194172 : DO i = 1, M_dim
654 711897 : DO ider = 1, dimders
655 517725 : 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 517725 : found=found)
660 517725 : CPASSERT(found .AND. ASSOCIATED(mom_block_der(i, ider)%block))
661 40655619 : 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 284811 : IF (PRESENT(ref_point)) THEN
669 1139244 : rc(:) = ref_point(:)
670 : ELSE
671 0 : 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 284811 : IF (my_neighbor_image) THEN
676 281199 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
677 1124796 : rb(:) = ra(:) + rab(:)
678 281199 : rac(:) = pbc(rc, ra, cell)
679 1124796 : rbc(:) = rac(:) + rab(:)
680 3612 : ELSE IF (my_ordered) THEN
681 9072 : ra(:) = particle_set(iatom)%r(:)
682 9072 : rb(:) = particle_set(jatom)%r(:)
683 2268 : rab(:) = pbc(rb, ra, cell)
684 9072 : rac(:) = pbc(ra - rc, cell)
685 9072 : rbc(:) = pbc(rb - rc, cell)
686 : ELSE
687 10752 : ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
688 10752 : rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
689 : ! we dont use PBC at this point
690 5376 : rab(:) = ra(:) - rb(:)
691 5376 : rac(:) = ra(:) - rc(:)
692 5376 : rbc(:) = rb(:) - rc(:)
693 : END IF
694 1139244 : dab = NORM2(rab)
695 :
696 1243365 : DO iset = 1, nseta
697 :
698 654840 : ncoa = npgfa(iset)*ncoset(la_max(iset))
699 654840 : sgfa = first_sgfa(1, iset)
700 :
701 2522697 : DO jset = 1, nsetb
702 :
703 1583046 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
704 :
705 600231 : ncob = npgfb(jset)*ncoset(lb_max(jset))
706 600231 : sgfb = first_sgfb(1, jset)
707 :
708 600231 : NULLIFY (mab_tmp)
709 : ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
710 2974323 : npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(nmom_build) - 1))
711 :
712 : ! Calculate the primitive integrals (need l+1 for derivatives)
713 600231 : 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 600231 : rpgfb(:, jset), nmom_build, rac, rbc, mab_tmp)
718 :
719 600231 : 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 600231 : lambda=lambda, deltaR=deltaR, iatom=iatom, jatom=jatom)
726 : END IF
727 :
728 600231 : IF (nmoments > 0) THEN
729 : ! copy subset of mab_tmp (l+1) to mab (l)
730 597930 : mab = 0.0_dp
731 5979300 : DO ii = 1, nm
732 5381370 : na = 0
733 5381370 : nda = 0
734 17745306 : DO ipgf = 1, npgfa(iset)
735 11766006 : nb = 0
736 11766006 : ndb = 0
737 43098966 : DO jpgf = 1, npgfb(jset)
738 128571651 : DO j = 1, ncoset(lb_max(jset))
739 488676321 : DO i = 1, ncoset(la_max(iset))
740 457343361 : mab(i + na, j + nb, ii) = mab_tmp(i + nda, j + ndb, ii)
741 : END DO ! i
742 : END DO ! j
743 31332960 : nb = nb + ncoset(lb_max(jset))
744 43098966 : ndb = ndb + ncoset(lb_max(jset) + 1)
745 : END DO ! jpgf
746 11766006 : na = na + ncoset(la_max(iset))
747 17147376 : nda = nda + ncoset(la_max(iset) + 1)
748 : END DO ! ipgf
749 : END DO
750 : ! Contraction step
751 5979300 : 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 5381370 : 0.0_dp, work(1, 1), SIZE(work, 1))
757 :
758 5979300 : 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 5377653 : 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 3717 : SIZE(mom_block(i)%block, 1))
770 : END IF
771 : END DO
772 : END IF
773 :
774 600231 : IF (nmoments_der > 0) THEN
775 5982852 : DO i = 1, M_dim
776 22130715 : 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 16147863 : 0._dp, work(1, 1), SIZE(work, 1))
781 :
782 21530484 : 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 16144092 : 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 3771 : SIZE(mom_block_der(i, ider)%block, 1))
794 : END IF
795 : END DO
796 : END DO
797 : END IF
798 2237886 : DEALLOCATE (mab_tmp)
799 : END DO
800 : END DO
801 : END ASSOCIATE
802 : END DO
803 4302 : CALL neighbor_list_iterator_release(nl_iterator)
804 :
805 : ! deallocations
806 4302 : DEALLOCATE (basis_set_list)
807 4302 : DEALLOCATE (work)
808 4302 : IF (nmoments > 0) THEN
809 3790 : DEALLOCATE (mab)
810 37900 : DO i = 1, nm
811 37900 : NULLIFY (mom_block(i)%block)
812 : END DO
813 3790 : DEALLOCATE (mom_block)
814 : END IF
815 4302 : IF (nmoments_der > 0) THEN
816 4302 : DEALLOCATE (difmab)
817 39720 : DO i = 1, M_dim
818 145974 : DO ider = 1, dimders
819 141672 : NULLIFY (mom_block_der(i, ider)%block)
820 : END DO
821 : END DO
822 4302 : DEALLOCATE (mom_block_der)
823 : END IF
824 :
825 4302 : CALL timestop(handle)
826 :
827 8604 : 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 64 : 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 64 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
855 64 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
856 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
857 64 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
858 : TYPE(cell_type), POINTER :: cell
859 64 : 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 64 : DIMENSION(:), POINTER :: nl_iterator
863 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
864 64 : POINTER :: sab_orb
865 64 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
866 64 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
867 : TYPE(qs_kind_type), POINTER :: qs_kind
868 :
869 64 : IF (nmoments < 1) RETURN
870 :
871 64 : CALL timeset(routineN, handle)
872 :
873 : ! magnetic dipoles/angular moments only
874 64 : nm = 3
875 :
876 64 : 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 64 : sab_orb=sab_orb)
881 :
882 64 : nkind = SIZE(qs_kind_set)
883 :
884 : ! Allocate work storage
885 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
886 64 : maxco=maxco, maxsgf=maxsgf)
887 :
888 320 : ALLOCATE (mab(maxco, maxco, nm))
889 64 : mab(:, :, :) = 0.0_dp
890 :
891 256 : ALLOCATE (work(maxco, maxsgf))
892 64 : work(:, :) = 0.0_dp
893 :
894 256 : ALLOCATE (mint(nm))
895 256 : DO i = 1, nm
896 256 : NULLIFY (mint(i)%block)
897 : END DO
898 :
899 356 : ALLOCATE (basis_set_list(nkind))
900 228 : DO ikind = 1, nkind
901 164 : qs_kind => qs_kind_set(ikind)
902 328 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
903 228 : IF (ASSOCIATED(basis_set_a)) THEN
904 164 : basis_set_list(ikind)%gto_basis_set => basis_set_a
905 : ELSE
906 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
907 : END IF
908 : END DO
909 64 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
910 7696 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
911 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
912 7632 : iatom=iatom, jatom=jatom, r=rab)
913 7632 : basis_set_a => basis_set_list(ikind)%gto_basis_set
914 7632 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
915 7632 : basis_set_b => basis_set_list(jkind)%gto_basis_set
916 7632 : 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 7632 : nseta = basis_set_a%nset
940 7632 : nsetb = basis_set_b%nset
941 :
942 7632 : IF (iatom <= jatom) THEN
943 4216 : irow = iatom
944 4216 : icol = jatom
945 : ELSE
946 3416 : irow = jatom
947 3416 : icol = iatom
948 : END IF
949 :
950 30528 : DO i = 1, nm
951 22896 : NULLIFY (mint(i)%block)
952 : CALL dbcsr_get_block_p(matrix=magmom(i)%matrix, &
953 22896 : row=irow, col=icol, BLOCK=mint(i)%block, found=found)
954 22896 : CPASSERT(found .AND. ASSOCIATED(mint(i)%block))
955 384762 : mint(i)%block = 0._dp
956 : END DO
957 :
958 : ! fold atomic position back into unit cell
959 7632 : IF (PRESENT(ref_points)) THEN
960 0 : rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
961 7632 : ELSE IF (PRESENT(ref_point)) THEN
962 30528 : rc(:) = ref_point(:)
963 : ELSE
964 0 : 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 61056 : ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
969 61056 : rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
970 : ! we dont use PBC at this point
971 30528 : rab(:) = ra(:) - rb(:)
972 30528 : rac(:) = ra(:) - rc(:)
973 30528 : rbc(:) = rb(:) - rc(:)
974 30528 : dab = NORM2(rab)
975 :
976 23013 : DO iset = 1, nseta
977 :
978 7749 : ncoa = npgfa(iset)*ncoset(la_max(iset))
979 7749 : sgfa = first_sgfa(1, iset)
980 :
981 23364 : DO jset = 1, nsetb
982 :
983 7983 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
984 :
985 6990 : ncob = npgfb(jset)*ncoset(lb_max(jset))
986 6990 : 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 6990 : rpgfb(:, jset), rac, rbc, mab)
993 :
994 : ! Contraction step
995 35709 : 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 20970 : 0.0_dp, work(1, 1), SIZE(work, 1))
1000 :
1001 28953 : 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 11901 : 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 9069 : 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 64 : CALL neighbor_list_iterator_release(nl_iterator)
1022 :
1023 : ! Release work storage
1024 64 : DEALLOCATE (mab, basis_set_list)
1025 64 : DEALLOCATE (work)
1026 256 : DO i = 1, nm
1027 256 : NULLIFY (mint(i)%block)
1028 : END DO
1029 64 : DEALLOCATE (mint)
1030 :
1031 64 : CALL timestop(handle)
1032 :
1033 192 : 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 13456 : 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 53824 : TYPE(dbcsr_p_type), DIMENSION(2, 1) :: op_sm_set
1057 : REAL(KIND=dp), DIMENSION(3, 1) :: kvec_batch
1058 :
1059 13456 : CALL timeset(routineN, handle)
1060 :
1061 13456 : op_sm_set(1, 1)%matrix => cosmat
1062 13456 : op_sm_set(2, 1)%matrix => sinmat
1063 13456 : kvec_batch(:, 1) = kvec(:)
1064 :
1065 : CALL build_exp_ikr_matrix(qs_env, op_sm_set, kvec_batch, &
1066 15752 : sab_orb_external=sab_orb_external, basis_type=basis_type)
1067 :
1068 13456 : CALL timestop(handle)
1069 :
1070 13456 : 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 : !> \param ordered use all ordered atom pairs for a non-Hermitian cross-k overlap
1080 : !> \param ordered_kpoints image-cell mapping for the ordered neighbour list
1081 : ! **************************************************************************************************
1082 116 : SUBROUTINE build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type, ordered, ordered_kpoints)
1083 :
1084 : TYPE(qs_environment_type), POINTER :: qs_env
1085 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: cosmat, sinmat
1086 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: kvec
1087 : CHARACTER(len=*), OPTIONAL :: basis_type
1088 : LOGICAL, INTENT(IN), OPTIONAL :: ordered
1089 : TYPE(kpoint_type), POINTER, OPTIONAL :: ordered_kpoints
1090 :
1091 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_berry_kpoint_matrix'
1092 :
1093 : CHARACTER :: matrix_symmetry
1094 : INTEGER :: handle, i, iatom, ic, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, &
1095 : ldsa, ldsb, ldwork, natom, ncoa, ncob, nimg, nkind, nseta, nsetb, sgfa, sgfb
1096 : INTEGER, DIMENSION(3) :: icell
1097 116 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
1098 116 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1099 : LOGICAL :: found, full_pairs, use_cell_mapping
1100 116 : REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1101 : REAL(KIND=dp) :: dab
1102 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb
1103 : TYPE(cell_type), POINTER :: cell
1104 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1105 : TYPE(dft_control_type), POINTER :: dft_control
1106 116 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1107 : TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
1108 : TYPE(kpoint_type), POINTER :: kpoints
1109 : TYPE(neighbor_list_iterator_p_type), &
1110 116 : DIMENSION(:), POINTER :: nl_iterator
1111 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1112 116 : POINTER :: sab_orb
1113 116 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1114 116 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1115 : TYPE(qs_kind_type), POINTER :: qs_kind
1116 : TYPE(qs_ks_env_type), POINTER :: ks_env
1117 :
1118 116 : CALL timeset(routineN, handle)
1119 :
1120 116 : full_pairs = .FALSE.
1121 116 : IF (PRESENT(ordered)) full_pairs = ordered
1122 116 : matrix_symmetry = dbcsr_type_symmetric
1123 116 : IF (full_pairs) matrix_symmetry = dbcsr_type_no_symmetry
1124 :
1125 : CALL get_qs_env(qs_env, &
1126 : ks_env=ks_env, &
1127 116 : dft_control=dft_control)
1128 116 : nimg = dft_control%nimages
1129 116 : IF (nimg > 1) THEN
1130 116 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1131 116 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1132 116 : use_cell_mapping = .TRUE.
1133 : ELSE
1134 : use_cell_mapping = .FALSE.
1135 : END IF
1136 :
1137 : CALL get_qs_env(qs_env=qs_env, &
1138 : qs_kind_set=qs_kind_set, &
1139 : particle_set=particle_set, cell=cell, &
1140 116 : sab_orb=sab_orb)
1141 :
1142 116 : IF (full_pairs) THEN
1143 20 : CALL get_qs_env(qs_env, sab_all=sab_orb)
1144 20 : IF (.NOT. ASSOCIATED(sab_orb)) CPABORT("Ordered Berry integrals require k-point neighbour lists.")
1145 20 : CPASSERT(PRESENT(ordered_kpoints))
1146 20 : CALL get_kpoint_info(ordered_kpoints, cell_to_index=cell_to_index)
1147 20 : nimg = SIZE(ordered_kpoints%index_to_cell, 2)
1148 20 : use_cell_mapping = .TRUE.
1149 : END IF
1150 :
1151 116 : nkind = SIZE(qs_kind_set)
1152 116 : natom = SIZE(particle_set)
1153 464 : ALLOCATE (basis_set_list(nkind))
1154 232 : DO ikind = 1, nkind
1155 116 : qs_kind => qs_kind_set(ikind)
1156 232 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=basis_type)
1157 232 : IF (ASSOCIATED(basis_set)) THEN
1158 116 : basis_set_list(ikind)%gto_basis_set => basis_set
1159 : ELSE
1160 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
1161 : END IF
1162 : END DO
1163 :
1164 348 : ALLOCATE (row_blk_sizes(natom))
1165 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1166 116 : basis=basis_set_list)
1167 116 : CALL get_ks_env(ks_env, dbcsr_dist=dbcsr_dist)
1168 : ! (re)allocate matrix sets
1169 116 : CALL dbcsr_allocate_matrix_set(sinmat, 1, nimg)
1170 116 : CALL dbcsr_allocate_matrix_set(cosmat, 1, nimg)
1171 17582 : DO i = 1, nimg
1172 : ! sin
1173 17466 : ALLOCATE (sinmat(1, i)%matrix)
1174 : CALL dbcsr_create(matrix=sinmat(1, i)%matrix, &
1175 : name="SINMAT", &
1176 : dist=dbcsr_dist, matrix_type=matrix_symmetry, &
1177 17466 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
1178 17466 : CALL cp_dbcsr_alloc_block_from_nbl(sinmat(1, i)%matrix, sab_orb)
1179 17466 : CALL dbcsr_set(sinmat(1, i)%matrix, 0.0_dp)
1180 : ! cos
1181 17466 : ALLOCATE (cosmat(1, i)%matrix)
1182 : CALL dbcsr_create(matrix=cosmat(1, i)%matrix, &
1183 : name="COSMAT", &
1184 : dist=dbcsr_dist, matrix_type=matrix_symmetry, &
1185 17466 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
1186 17466 : CALL cp_dbcsr_alloc_block_from_nbl(cosmat(1, i)%matrix, sab_orb)
1187 17582 : CALL dbcsr_set(cosmat(1, i)%matrix, 0.0_dp)
1188 : END DO
1189 :
1190 116 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
1191 116 : ldab = ldwork
1192 464 : ALLOCATE (cosab(ldab, ldab))
1193 348 : ALLOCATE (sinab(ldab, ldab))
1194 348 : ALLOCATE (work(ldwork, ldwork))
1195 :
1196 116 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1197 98908 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1198 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1199 98792 : iatom=iatom, jatom=jatom, r=rab, cell=icell)
1200 98792 : basis_set_a => basis_set_list(ikind)%gto_basis_set
1201 98792 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
1202 98792 : basis_set_b => basis_set_list(jkind)%gto_basis_set
1203 98792 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
1204 : ASSOCIATE ( &
1205 : ! basis ikind
1206 : first_sgfa => basis_set_a%first_sgf, &
1207 : la_max => basis_set_a%lmax, &
1208 : la_min => basis_set_a%lmin, &
1209 : npgfa => basis_set_a%npgf, &
1210 : nsgfa => basis_set_a%nsgf_set, &
1211 : rpgfa => basis_set_a%pgf_radius, &
1212 : set_radius_a => basis_set_a%set_radius, &
1213 : sphi_a => basis_set_a%sphi, &
1214 : zeta => basis_set_a%zet, &
1215 : ! basis jkind, &
1216 : first_sgfb => basis_set_b%first_sgf, &
1217 : lb_max => basis_set_b%lmax, &
1218 : lb_min => basis_set_b%lmin, &
1219 : npgfb => basis_set_b%npgf, &
1220 : nsgfb => basis_set_b%nsgf_set, &
1221 : rpgfb => basis_set_b%pgf_radius, &
1222 : set_radius_b => basis_set_b%set_radius, &
1223 : sphi_b => basis_set_b%sphi, &
1224 197584 : zetb => basis_set_b%zet)
1225 :
1226 98792 : nseta = basis_set_a%nset
1227 98792 : nsetb = basis_set_b%nset
1228 :
1229 98792 : ldsa = SIZE(sphi_a, 1)
1230 98792 : ldsb = SIZE(sphi_b, 1)
1231 :
1232 42450 : IF (iatom <= jatom .OR. full_pairs) THEN
1233 56552 : irow = iatom
1234 56552 : icol = jatom
1235 : ELSE
1236 42240 : irow = jatom
1237 42240 : icol = iatom
1238 : END IF
1239 :
1240 98792 : IF (use_cell_mapping) THEN
1241 98792 : ic = cell_to_index(icell(1), icell(2), icell(3))
1242 98792 : CPASSERT(ic > 0)
1243 : ELSE
1244 : ic = 1
1245 : END IF
1246 :
1247 98792 : NULLIFY (sblock)
1248 : CALL dbcsr_get_block_p(matrix=sinmat(1, ic)%matrix, &
1249 98792 : row=irow, col=icol, BLOCK=sblock, found=found)
1250 98792 : CPASSERT(found)
1251 98792 : NULLIFY (cblock)
1252 : CALL dbcsr_get_block_p(matrix=cosmat(1, ic)%matrix, &
1253 98792 : row=irow, col=icol, BLOCK=cblock, found=found)
1254 98792 : CPASSERT(found)
1255 :
1256 98792 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
1257 395168 : rb(:) = ra + rab
1258 98792 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1259 :
1260 300038 : DO iset = 1, nseta
1261 :
1262 102454 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1263 102454 : sgfa = first_sgfa(1, iset)
1264 :
1265 311024 : DO jset = 1, nsetb
1266 :
1267 109778 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
1268 :
1269 106986 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1270 106986 : sgfb = first_sgfb(1, jset)
1271 :
1272 : ! Calculate the primitive integrals
1273 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1274 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1275 106986 : ra, rb, kvec, cosab, sinab)
1276 : CALL contract_cossin(cblock, sblock, &
1277 : iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1278 : jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1279 212232 : cosab, sinab, ldab, work, ldwork, ordered=full_pairs)
1280 :
1281 : END DO
1282 : END DO
1283 : END ASSOCIATE
1284 : END DO
1285 116 : CALL neighbor_list_iterator_release(nl_iterator)
1286 :
1287 116 : DEALLOCATE (cosab)
1288 116 : DEALLOCATE (sinab)
1289 116 : DEALLOCATE (work)
1290 116 : DEALLOCATE (basis_set_list)
1291 116 : DEALLOCATE (row_blk_sizes)
1292 :
1293 116 : CALL timestop(handle)
1294 :
1295 232 : END SUBROUTINE build_berry_kpoint_matrix
1296 :
1297 : ! **************************************************************************************************
1298 : !> \brief ...
1299 : !> \param qs_env ...
1300 : !> \param magnetic ...
1301 : !> \param nmoments ...
1302 : !> \param reference ...
1303 : !> \param ref_point ...
1304 : !> \param unit_number ...
1305 : ! **************************************************************************************************
1306 478 : SUBROUTINE qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number)
1307 :
1308 : TYPE(qs_environment_type), POINTER :: qs_env
1309 : LOGICAL, INTENT(IN) :: magnetic
1310 : INTEGER, INTENT(IN) :: nmoments, reference
1311 : REAL(dp), DIMENSION(:), POINTER :: ref_point
1312 : INTEGER, INTENT(IN) :: unit_number
1313 :
1314 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_moment_berry_phase'
1315 :
1316 478 : CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
1317 : CHARACTER(LEN=default_string_length) :: description
1318 : COMPLEX(dp) :: xphase(3), zdet, zdeta, zi(3), &
1319 : zij(3, 3), zijk(3, 3, 3), &
1320 : zijkl(3, 3, 3, 3), zphase(3), zz
1321 : INTEGER :: handle, i, ia, idim, ikind, ispin, ix, &
1322 : iy, iz, j, k, l, nao, nm, nmo, nmom, &
1323 : nmotot, tmp_dim
1324 : LOGICAL :: floating, ghost, uniform
1325 : REAL(dp) :: charge, ci(3), cij(3, 3), dd, occ, trace
1326 478 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom
1327 478 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
1328 : REAL(dp), DIMENSION(3) :: kvec, qq, rcc, ria
1329 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1330 : TYPE(cell_type), POINTER :: cell
1331 478 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat
1332 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
1333 478 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: opvec
1334 478 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: op_fm_set
1335 478 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
1336 : TYPE(cp_fm_type), POINTER :: mo_coeff
1337 : TYPE(cp_result_type), POINTER :: results
1338 478 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
1339 : TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1340 : TYPE(dft_control_type), POINTER :: dft_control
1341 : TYPE(distribution_1d_type), POINTER :: local_particles
1342 478 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1343 : TYPE(mp_para_env_type), POINTER :: para_env
1344 478 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1345 478 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1346 : TYPE(qs_rho_type), POINTER :: rho
1347 : TYPE(rt_prop_type), POINTER :: rtp
1348 :
1349 0 : CPASSERT(ASSOCIATED(qs_env))
1350 :
1351 478 : IF (ASSOCIATED(qs_env%ls_scf_env)) THEN
1352 0 : IF (unit_number > 0) WRITE (unit_number, *) "Periodic moment calculation not implemented in linear scaling code"
1353 0 : RETURN
1354 : END IF
1355 :
1356 478 : CALL timeset(routineN, handle)
1357 :
1358 : ! restrict maximum moment available
1359 478 : nmom = MIN(nmoments, 2)
1360 :
1361 478 : nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
1362 : ! rmom(:,1)=electronic
1363 : ! rmom(:,2)=nuclear
1364 : ! rmom(:,1)=total
1365 2390 : ALLOCATE (rmom(nm + 1, 3))
1366 1434 : ALLOCATE (rlab(nm + 1))
1367 478 : rmom = 0.0_dp
1368 2390 : rlab = ""
1369 478 : IF (magnetic) THEN
1370 0 : nm = 3
1371 0 : ALLOCATE (mmom(nm))
1372 0 : mmom = 0._dp
1373 : END IF
1374 :
1375 478 : NULLIFY (dft_control, rho, cell, particle_set, results, para_env, &
1376 478 : local_particles, matrix_s, mos, rho_ao)
1377 :
1378 : CALL get_qs_env(qs_env, &
1379 : dft_control=dft_control, &
1380 : rho=rho, &
1381 : cell=cell, &
1382 : results=results, &
1383 : particle_set=particle_set, &
1384 : qs_kind_set=qs_kind_set, &
1385 : para_env=para_env, &
1386 : local_particles=local_particles, &
1387 : matrix_s=matrix_s, &
1388 478 : mos=mos)
1389 :
1390 478 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1391 :
1392 478 : NULLIFY (cosmat, sinmat)
1393 478 : ALLOCATE (cosmat, sinmat)
1394 478 : CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
1395 478 : CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
1396 :
1397 2934 : ALLOCATE (op_fm_set(2, dft_control%nspins))
1398 1934 : ALLOCATE (opvec(dft_control%nspins))
1399 1934 : ALLOCATE (eigrmat(dft_control%nspins))
1400 478 : nmotot = 0
1401 978 : DO ispin = 1, dft_control%nspins
1402 500 : NULLIFY (tmp_fm_struct, mo_coeff)
1403 500 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
1404 500 : nmotot = nmotot + nmo
1405 500 : CALL cp_fm_create(opvec(ispin), mo_coeff%matrix_struct)
1406 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
1407 500 : ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
1408 1500 : DO i = 1, SIZE(op_fm_set, 1)
1409 1500 : CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
1410 : END DO
1411 500 : CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
1412 1478 : CALL cp_fm_struct_release(tmp_fm_struct)
1413 : END DO
1414 :
1415 : ! occupation
1416 978 : DO ispin = 1, dft_control%nspins
1417 500 : CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
1418 978 : IF (.NOT. uniform) THEN
1419 0 : CPWARN("Berry phase moments for non uniform MOs' occupation numbers not implemented")
1420 : END IF
1421 : END DO
1422 :
1423 : ! reference point
1424 478 : CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
1425 1912 : rcc = pbc(rcc, cell)
1426 :
1427 : ! label
1428 1912 : DO l = 1, nm
1429 1434 : ix = indco(1, l + 1)
1430 1434 : iy = indco(2, l + 1)
1431 1434 : iz = indco(3, l + 1)
1432 1912 : CALL set_label(rlab(l + 1), ix, iy, iz)
1433 : END DO
1434 :
1435 : ! nuclear contribution
1436 1932 : DO ia = 1, SIZE(particle_set)
1437 1454 : atomic_kind => particle_set(ia)%atomic_kind
1438 1454 : CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1439 1454 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1440 1932 : IF (.NOT. ghost .AND. .NOT. floating) THEN
1441 1454 : rmom(1, 2) = rmom(1, 2) - charge
1442 : END IF
1443 : END DO
1444 7648 : ria = twopi*MATMUL(cell%h_inv, rcc)
1445 1912 : zphase = CMPLX(COS(ria), SIN(ria), dp)**rmom(1, 2)
1446 :
1447 478 : zi = 0._dp
1448 478 : zij = 0._dp
1449 : zijk = 0._dp
1450 : zijkl = 0._dp
1451 :
1452 956 : DO l = 1, nmom
1453 478 : SELECT CASE (l)
1454 : CASE (1)
1455 : ! Dipole
1456 1912 : zi(:) = CMPLX(1._dp, 0._dp, dp)
1457 1932 : DO ia = 1, SIZE(particle_set)
1458 1454 : atomic_kind => particle_set(ia)%atomic_kind
1459 1454 : CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1460 1454 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1461 1932 : IF (.NOT. ghost .AND. .NOT. floating) THEN
1462 5816 : ria = particle_set(ia)%r
1463 5816 : ria = pbc(ria, cell)
1464 5816 : DO i = 1, 3
1465 17448 : kvec(:) = twopi*cell%h_inv(i, :)
1466 17448 : dd = SUM(kvec(:)*ria(:))
1467 4362 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
1468 5816 : zi(i) = zi(i)*zdeta
1469 : END DO
1470 : END IF
1471 : END DO
1472 1912 : zi = zi*zphase
1473 1912 : ci = AIMAG(LOG(zi))/twopi
1474 1912 : qq = AIMAG(LOG(zi))
1475 7648 : rmom(2:4, 2) = MATMUL(cell%hmat, ci)
1476 : CASE (2)
1477 : ! Quadrupole
1478 0 : CPABORT("Berry phase moments bigger than 1 not implemented")
1479 0 : zij(:, :) = CMPLX(1._dp, 0._dp, dp)
1480 0 : DO ia = 1, SIZE(particle_set)
1481 0 : atomic_kind => particle_set(ia)%atomic_kind
1482 0 : CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1483 0 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
1484 0 : ria = particle_set(ia)%r
1485 0 : ria = pbc(ria, cell)
1486 0 : DO i = 1, 3
1487 0 : DO j = i, 3
1488 0 : kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1489 0 : dd = SUM(kvec(:)*ria(:))
1490 0 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
1491 0 : zij(i, j) = zij(i, j)*zdeta
1492 0 : zij(j, i) = zij(i, j)
1493 : END DO
1494 : END DO
1495 : END DO
1496 0 : DO i = 1, 3
1497 0 : DO j = 1, 3
1498 0 : zij(i, j) = zij(i, j)*zphase(i)*zphase(j)
1499 0 : zz = zij(i, j)/zi(i)/zi(j)
1500 0 : cij(i, j) = AIMAG(LOG(zz))/twopi
1501 : END DO
1502 : END DO
1503 0 : cij = 0.5_dp*cij/twopi/twopi
1504 0 : cij = MATMUL(MATMUL(cell%hmat, cij), TRANSPOSE(cell%hmat))
1505 0 : DO k = 4, 9
1506 0 : ix = indco(1, k + 1)
1507 0 : iy = indco(2, k + 1)
1508 0 : iz = indco(3, k + 1)
1509 0 : IF (ix == 0) THEN
1510 0 : rmom(k + 1, 2) = cij(iy, iz)
1511 0 : ELSE IF (iy == 0) THEN
1512 0 : rmom(k + 1, 2) = cij(ix, iz)
1513 0 : ELSE IF (iz == 0) THEN
1514 0 : rmom(k + 1, 2) = cij(ix, iy)
1515 : END IF
1516 : END DO
1517 : CASE (3)
1518 : ! Octapole
1519 0 : CPABORT("Berry phase moments bigger than 2 not implemented")
1520 : CASE (4)
1521 : ! Hexadecapole
1522 0 : CPABORT("Berry phase moments bigger than 3 not implemented")
1523 : CASE DEFAULT
1524 478 : CPABORT("Berry phase moments bigger than 4 not implemented")
1525 : END SELECT
1526 : END DO
1527 :
1528 : ! electronic contribution
1529 :
1530 7648 : ria = twopi*REAL(nmotot, dp)*occ*MATMUL(cell%h_inv, rcc)
1531 1912 : xphase = CMPLX(COS(ria), SIN(ria), dp)
1532 :
1533 : ! charge
1534 478 : trace = 0.0_dp
1535 978 : DO ispin = 1, dft_control%nspins
1536 500 : CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
1537 978 : rmom(1, 1) = rmom(1, 1) + trace
1538 : END DO
1539 :
1540 478 : zi = 0._dp
1541 478 : zij = 0._dp
1542 : zijk = 0._dp
1543 : zijkl = 0._dp
1544 :
1545 956 : DO l = 1, nmom
1546 478 : SELECT CASE (l)
1547 : CASE (1)
1548 : ! Dipole
1549 1912 : DO i = 1, 3
1550 5736 : kvec(:) = twopi*cell%h_inv(i, :)
1551 1434 : CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1552 1434 : IF (qs_env%run_rtp) THEN
1553 48 : CALL get_qs_env(qs_env, rtp=rtp)
1554 48 : CALL get_rtp(rtp, mos_new=mos_new)
1555 48 : CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1556 : ELSE
1557 1386 : CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1558 : END IF
1559 1434 : zdet = CMPLX(1._dp, 0._dp, dp)
1560 2934 : DO ispin = 1, dft_control%nspins
1561 1500 : CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
1562 8322 : DO idim = 1, tmp_dim
1563 : eigrmat(ispin)%local_data(:, idim) = &
1564 : CMPLX(op_fm_set(1, ispin)%local_data(:, idim), &
1565 43968 : -op_fm_set(2, ispin)%local_data(:, idim), dp)
1566 : END DO
1567 : ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
1568 1500 : CALL cp_cfm_det(eigrmat(ispin), zdeta)
1569 1500 : zdet = zdet*zdeta
1570 4434 : IF (dft_control%nspins == 1) THEN
1571 1368 : zdet = zdet*zdeta
1572 : END IF
1573 : END DO
1574 1912 : zi(i) = zdet
1575 : END DO
1576 1912 : zi = zi*xphase
1577 : CASE (2)
1578 : ! Quadrupole
1579 0 : CPABORT("Berry phase moments bigger than 1 not implemented")
1580 0 : DO i = 1, 3
1581 0 : DO j = i, 3
1582 0 : kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1583 0 : CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1584 0 : IF (qs_env%run_rtp) THEN
1585 0 : CALL get_qs_env(qs_env, rtp=rtp)
1586 0 : CALL get_rtp(rtp, mos_new=mos_new)
1587 0 : CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1588 : ELSE
1589 0 : CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1590 : END IF
1591 0 : zdet = CMPLX(1._dp, 0._dp, dp)
1592 0 : DO ispin = 1, dft_control%nspins
1593 0 : CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
1594 0 : DO idim = 1, tmp_dim
1595 : eigrmat(ispin)%local_data(:, idim) = &
1596 : CMPLX(op_fm_set(1, ispin)%local_data(:, idim), &
1597 0 : -op_fm_set(2, ispin)%local_data(:, idim), dp)
1598 : END DO
1599 : ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
1600 0 : CALL cp_cfm_det(eigrmat(ispin), zdeta)
1601 0 : zdet = zdet*zdeta
1602 0 : IF (dft_control%nspins == 1) THEN
1603 0 : zdet = zdet*zdeta
1604 : END IF
1605 : END DO
1606 0 : zij(i, j) = zdet*xphase(i)*xphase(j)
1607 0 : zij(j, i) = zdet*xphase(i)*xphase(j)
1608 : END DO
1609 : END DO
1610 : CASE (3)
1611 : ! Octapole
1612 0 : CPABORT("Berry phase moments bigger than 2 not implemented")
1613 : CASE (4)
1614 : ! Hexadecapole
1615 0 : CPABORT("Berry phase moments bigger than 3 not implemented")
1616 : CASE DEFAULT
1617 478 : CPABORT("Berry phase moments bigger than 4 not implemented")
1618 : END SELECT
1619 : END DO
1620 956 : DO l = 1, nmom
1621 478 : SELECT CASE (l)
1622 : CASE (1)
1623 : ! Dipole (apply periodic (2 Pi) boundary conditions)
1624 1912 : ci = AIMAG(LOG(zi))
1625 1912 : DO i = 1, 3
1626 1434 : IF (qq(i) + ci(i) > pi) ci(i) = ci(i) - twopi
1627 1912 : IF (qq(i) + ci(i) < -pi) ci(i) = ci(i) + twopi
1628 : END DO
1629 9082 : rmom(2:4, 1) = MATMUL(cell%hmat, ci)/twopi
1630 : CASE (2)
1631 : ! Quadrupole
1632 0 : CPABORT("Berry phase moments bigger than 1 not implemented")
1633 0 : DO i = 1, 3
1634 0 : DO j = 1, 3
1635 0 : zz = zij(i, j)/zi(i)/zi(j)
1636 0 : cij(i, j) = AIMAG(LOG(zz))/twopi
1637 : END DO
1638 : END DO
1639 0 : cij = 0.5_dp*cij/twopi/twopi
1640 0 : cij = MATMUL(MATMUL(cell%hmat, cij), TRANSPOSE(cell%hmat))
1641 0 : DO k = 4, 9
1642 0 : ix = indco(1, k + 1)
1643 0 : iy = indco(2, k + 1)
1644 0 : iz = indco(3, k + 1)
1645 0 : IF (ix == 0) THEN
1646 0 : rmom(k + 1, 1) = cij(iy, iz)
1647 0 : ELSE IF (iy == 0) THEN
1648 0 : rmom(k + 1, 1) = cij(ix, iz)
1649 0 : ELSE IF (iz == 0) THEN
1650 0 : rmom(k + 1, 1) = cij(ix, iy)
1651 : END IF
1652 : END DO
1653 : CASE (3)
1654 : ! Octapole
1655 0 : CPABORT("Berry phase moments bigger than 2 not implemented")
1656 : CASE (4)
1657 : ! Hexadecapole
1658 0 : CPABORT("Berry phase moments bigger than 3 not implemented")
1659 : CASE DEFAULT
1660 478 : CPABORT("Berry phase moments bigger than 4 not implemented")
1661 : END SELECT
1662 : END DO
1663 :
1664 2390 : rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
1665 478 : description = "[DIPOLE]"
1666 478 : CALL cp_results_erase(results=results, description=description)
1667 : CALL put_results(results=results, description=description, &
1668 478 : values=rmom(2:4, 3))
1669 478 : IF (magnetic) THEN
1670 0 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.TRUE., mmom=mmom)
1671 : ELSE
1672 478 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.TRUE.)
1673 : END IF
1674 :
1675 478 : DEALLOCATE (rmom)
1676 478 : DEALLOCATE (rlab)
1677 478 : IF (magnetic) THEN
1678 0 : DEALLOCATE (mmom)
1679 : END IF
1680 :
1681 478 : CALL dbcsr_deallocate_matrix(cosmat)
1682 478 : CALL dbcsr_deallocate_matrix(sinmat)
1683 :
1684 478 : CALL cp_fm_release(opvec)
1685 478 : CALL cp_fm_release(op_fm_set)
1686 978 : DO ispin = 1, dft_control%nspins
1687 978 : CALL cp_cfm_release(eigrmat(ispin))
1688 : END DO
1689 478 : DEALLOCATE (eigrmat)
1690 :
1691 478 : CALL timestop(handle)
1692 :
1693 1434 : END SUBROUTINE qs_moment_berry_phase
1694 :
1695 : ! **************************************************************************************************
1696 : !> \brief ...
1697 : !> \param cosmat ...
1698 : !> \param sinmat ...
1699 : !> \param mos ...
1700 : !> \param op_fm_set ...
1701 : !> \param opvec ...
1702 : ! **************************************************************************************************
1703 1386 : SUBROUTINE op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1704 :
1705 : TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1706 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1707 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
1708 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: opvec
1709 :
1710 : INTEGER :: i, nao, nmo
1711 : TYPE(cp_fm_type), POINTER :: mo_coeff
1712 :
1713 2814 : DO i = 1, SIZE(op_fm_set, 2) ! spin
1714 1428 : CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
1715 1428 : CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff, opvec(i), ncol=nmo)
1716 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
1717 1428 : op_fm_set(1, i))
1718 1428 : CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff, opvec(i), ncol=nmo)
1719 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
1720 4242 : op_fm_set(2, i))
1721 : END DO
1722 :
1723 1386 : END SUBROUTINE op_orbbas
1724 :
1725 : ! **************************************************************************************************
1726 : !> \brief ...
1727 : !> \param cosmat ...
1728 : !> \param sinmat ...
1729 : !> \param mos ...
1730 : !> \param op_fm_set ...
1731 : !> \param mos_new ...
1732 : ! **************************************************************************************************
1733 48 : SUBROUTINE op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1734 :
1735 : TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1736 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1737 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
1738 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
1739 :
1740 : INTEGER :: i, icol, lcol, nao, newdim, nmo
1741 : LOGICAL :: double_col, double_row
1742 : TYPE(cp_fm_struct_type), POINTER :: newstruct, newstruct1
1743 : TYPE(cp_fm_type) :: work, work1, work2
1744 : TYPE(cp_fm_type), POINTER :: mo_coeff
1745 :
1746 120 : DO i = 1, SIZE(op_fm_set, 2) ! spin
1747 72 : CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
1748 72 : CALL cp_fm_get_info(mos_new(2*i), ncol_local=lcol, ncol_global=nmo)
1749 72 : double_col = .TRUE.
1750 72 : double_row = .FALSE.
1751 : CALL cp_fm_struct_double(newstruct, &
1752 : mos_new(2*i)%matrix_struct, &
1753 : mos_new(2*i)%matrix_struct%context, &
1754 : double_col, &
1755 72 : double_row)
1756 :
1757 72 : CALL cp_fm_create(work, matrix_struct=newstruct)
1758 72 : CALL cp_fm_create(work1, matrix_struct=newstruct)
1759 72 : CALL cp_fm_create(work2, matrix_struct=newstruct)
1760 72 : CALL cp_fm_get_info(work, ncol_global=newdim)
1761 :
1762 72 : CALL cp_fm_set_all(work, 0.0_dp, 0.0_dp)
1763 336 : DO icol = 1, lcol
1764 3300 : work%local_data(:, icol) = mos_new(2*i - 1)%local_data(:, icol)
1765 3372 : work%local_data(:, icol + lcol) = mos_new(2*i)%local_data(:, icol)
1766 : END DO
1767 :
1768 72 : CALL cp_dbcsr_sm_fm_multiply(cosmat, work, work1, ncol=newdim)
1769 72 : CALL cp_dbcsr_sm_fm_multiply(sinmat, work, work2, ncol=newdim)
1770 :
1771 336 : DO icol = 1, lcol
1772 3300 : work%local_data(:, icol) = work1%local_data(:, icol) - work2%local_data(:, icol + lcol)
1773 3372 : work%local_data(:, icol + lcol) = work1%local_data(:, icol + lcol) + work2%local_data(:, icol)
1774 : END DO
1775 :
1776 72 : CALL cp_fm_release(work1)
1777 72 : CALL cp_fm_release(work2)
1778 :
1779 : CALL cp_fm_struct_double(newstruct1, &
1780 : op_fm_set(1, i)%matrix_struct, &
1781 : op_fm_set(1, i)%matrix_struct%context, &
1782 : double_col, &
1783 72 : double_row)
1784 :
1785 72 : CALL cp_fm_create(work1, matrix_struct=newstruct1)
1786 :
1787 : CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i - 1), &
1788 72 : work, 0.0_dp, work1)
1789 :
1790 336 : DO icol = 1, lcol
1791 756 : op_fm_set(1, i)%local_data(:, icol) = work1%local_data(:, icol)
1792 828 : op_fm_set(2, i)%local_data(:, icol) = work1%local_data(:, icol + lcol)
1793 : END DO
1794 :
1795 : CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i), &
1796 72 : work, 0.0_dp, work1)
1797 :
1798 336 : DO icol = 1, lcol
1799 : op_fm_set(1, i)%local_data(:, icol) = &
1800 756 : op_fm_set(1, i)%local_data(:, icol) + work1%local_data(:, icol + lcol)
1801 : op_fm_set(2, i)%local_data(:, icol) = &
1802 828 : op_fm_set(2, i)%local_data(:, icol) - work1%local_data(:, icol)
1803 : END DO
1804 :
1805 72 : CALL cp_fm_release(work)
1806 72 : CALL cp_fm_release(work1)
1807 72 : CALL cp_fm_struct_release(newstruct)
1808 336 : CALL cp_fm_struct_release(newstruct1)
1809 :
1810 : END DO
1811 :
1812 48 : END SUBROUTINE op_orbbas_rtp
1813 :
1814 : ! **************************************************************************************************
1815 : !> \brief ...
1816 : !> \param qs_env ...
1817 : !> \param magnetic ...
1818 : !> \param nmoments ...
1819 : !> \param reference ...
1820 : !> \param ref_point ...
1821 : !> \param unit_number ...
1822 : !> \param vel_reprs ...
1823 : !> \param com_nl ...
1824 : ! **************************************************************************************************
1825 1082 : SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
1826 :
1827 : TYPE(qs_environment_type), POINTER :: qs_env
1828 : LOGICAL, INTENT(IN) :: magnetic
1829 : INTEGER, INTENT(IN) :: nmoments, reference
1830 : REAL(dp), DIMENSION(:), INTENT(IN), POINTER :: ref_point
1831 : INTEGER, INTENT(IN) :: unit_number
1832 : LOGICAL, INTENT(IN), OPTIONAL :: vel_reprs, com_nl
1833 :
1834 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_moment_locop'
1835 :
1836 1082 : CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
1837 : CHARACTER(LEN=default_string_length) :: description
1838 : INTEGER :: akind, handle, i, ia, iatom, idir, &
1839 : ikind, ispin, ix, iy, iz, l, nm, nmom, &
1840 : order
1841 : LOGICAL :: my_com_nl, my_velreprs
1842 : REAL(dp) :: charge, dd, strace, trace
1843 1082 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rrv_vrr, &
1844 1082 : nlcom_rv, nlcom_rvr, nlcom_rxrv, &
1845 1082 : qupole_der, rmom_vel
1846 1082 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
1847 : REAL(dp), DIMENSION(3) :: rcc, ria
1848 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1849 : TYPE(cell_type), POINTER :: cell
1850 : TYPE(cp_result_type), POINTER :: results
1851 1082 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom, matrix_s, moments, momentum, &
1852 1082 : rho_ao
1853 1082 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_der
1854 : TYPE(dbcsr_type), POINTER :: tmp_ao
1855 : TYPE(dft_control_type), POINTER :: dft_control
1856 : TYPE(distribution_1d_type), POINTER :: local_particles
1857 : TYPE(mp_para_env_type), POINTER :: para_env
1858 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1859 1082 : POINTER :: sab_all, sab_orb
1860 1082 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1861 1082 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1862 : TYPE(qs_rho_type), POINTER :: rho
1863 :
1864 0 : CPASSERT(ASSOCIATED(qs_env))
1865 :
1866 1082 : CALL timeset(routineN, handle)
1867 :
1868 1082 : my_velreprs = .FALSE.
1869 1082 : IF (PRESENT(vel_reprs)) my_velreprs = vel_reprs
1870 1082 : IF (PRESENT(com_nl)) my_com_nl = com_nl
1871 1082 : IF (my_velreprs) CALL cite_reference(Mattiat2019)
1872 :
1873 : ! reference point
1874 1082 : CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
1875 :
1876 : ! only allow for moments up to maxl set by basis
1877 1082 : nmom = MIN(nmoments, current_maxl)
1878 : ! electronic contribution
1879 1082 : NULLIFY (dft_control, rho, cell, particle_set, qs_kind_set, results, para_env, matrix_s, rho_ao, sab_all, sab_orb)
1880 : CALL get_qs_env(qs_env, &
1881 : dft_control=dft_control, &
1882 : rho=rho, &
1883 : cell=cell, &
1884 : results=results, &
1885 : particle_set=particle_set, &
1886 : qs_kind_set=qs_kind_set, &
1887 : para_env=para_env, &
1888 : matrix_s=matrix_s, &
1889 : sab_all=sab_all, &
1890 1082 : sab_orb=sab_orb)
1891 :
1892 1082 : IF (my_com_nl) THEN
1893 40 : IF ((nmom >= 1) .AND. my_velreprs) THEN
1894 40 : ALLOCATE (nlcom_rv(3))
1895 40 : nlcom_rv(:) = 0._dp
1896 : END IF
1897 40 : IF ((nmom >= 2) .AND. my_velreprs) THEN
1898 40 : ALLOCATE (nlcom_rrv(6))
1899 40 : nlcom_rrv(:) = 0._dp
1900 40 : ALLOCATE (nlcom_rvr(6))
1901 40 : nlcom_rvr(:) = 0._dp
1902 40 : ALLOCATE (nlcom_rrv_vrr(6))
1903 40 : nlcom_rrv_vrr(:) = 0._dp
1904 : END IF
1905 40 : IF (magnetic) THEN
1906 18 : ALLOCATE (nlcom_rxrv(3))
1907 18 : nlcom_rxrv = 0._dp
1908 : END IF
1909 : ! Calculate non local correction terms
1910 40 : CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr, rcc)
1911 : END IF
1912 :
1913 1082 : NULLIFY (moments)
1914 1082 : nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
1915 1082 : CALL dbcsr_allocate_matrix_set(moments, nm)
1916 4630 : DO i = 1, nm
1917 3548 : ALLOCATE (moments(i)%matrix)
1918 3548 : IF (my_velreprs .AND. (nmom >= 2)) THEN
1919 : CALL dbcsr_create(moments(i)%matrix, template=matrix_s(1)%matrix, &
1920 360 : matrix_type=dbcsr_type_symmetric)
1921 360 : CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_orb)
1922 : ELSE
1923 3188 : CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
1924 : END IF
1925 4630 : CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
1926 : END DO
1927 :
1928 : ! calculate derivatives if quadrupole in vel. reprs. is requested
1929 1082 : IF (my_velreprs .AND. (nmom >= 2)) THEN
1930 40 : NULLIFY (moments_der)
1931 40 : CALL dbcsr_allocate_matrix_set(moments_der, 3, 3)
1932 160 : DO i = 1, 3 ! x, y, z
1933 520 : DO idir = 1, 3 ! d/dx, d/dy, d/dz
1934 360 : CALL dbcsr_init_p(moments_der(i, idir)%matrix)
1935 : CALL dbcsr_create(moments_der(i, idir)%matrix, template=matrix_s(1)%matrix, &
1936 360 : matrix_type=dbcsr_type_antisymmetric)
1937 360 : CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_orb)
1938 480 : CALL dbcsr_set(moments_der(i, idir)%matrix, 0.0_dp)
1939 : END DO
1940 : END DO
1941 40 : CALL build_local_moments_der_matrix(qs_env, moments_der, 1, 2, ref_point=rcc, moments=moments)
1942 : ELSE
1943 1042 : CALL build_local_moment_matrix(qs_env, moments, nmom, ref_point=rcc)
1944 : END IF
1945 :
1946 1082 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1947 :
1948 4328 : ALLOCATE (rmom(nm + 1, 3))
1949 3246 : ALLOCATE (rlab(nm + 1))
1950 1082 : rmom = 0.0_dp
1951 5712 : rlab = ""
1952 :
1953 1082 : IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
1954 : ! Allocate matrix to store the matrix product to be traced (dbcsr_dot only works for products of
1955 : ! symmetric matrices)
1956 42 : NULLIFY (tmp_ao)
1957 42 : CALL dbcsr_init_p(tmp_ao)
1958 42 : CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
1959 42 : CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
1960 42 : CALL dbcsr_set(tmp_ao, 0.0_dp)
1961 : END IF
1962 :
1963 1082 : trace = 0.0_dp
1964 2242 : DO ispin = 1, dft_control%nspins
1965 1160 : CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
1966 2242 : rmom(1, 1) = rmom(1, 1) + trace
1967 : END DO
1968 :
1969 4630 : DO i = 1, SIZE(moments)
1970 3548 : strace = 0._dp
1971 7330 : DO ispin = 1, dft_control%nspins
1972 3782 : IF (my_velreprs .AND. nmoments >= 2) THEN
1973 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, moments(i)%matrix, &
1974 360 : 0.0_dp, tmp_ao)
1975 360 : CALL dbcsr_trace(tmp_ao, trace)
1976 : ELSE
1977 3422 : CALL dbcsr_dot(rho_ao(ispin)%matrix, moments(i)%matrix, trace)
1978 : END IF
1979 7330 : strace = strace + trace
1980 : END DO
1981 4630 : rmom(i + 1, 1) = strace
1982 : END DO
1983 :
1984 1082 : CALL dbcsr_deallocate_matrix_set(moments)
1985 :
1986 : ! nuclear contribution
1987 : CALL get_qs_env(qs_env=qs_env, &
1988 1082 : local_particles=local_particles)
1989 3400 : DO ikind = 1, SIZE(local_particles%n_el)
1990 5178 : DO ia = 1, local_particles%n_el(ikind)
1991 1778 : iatom = local_particles%list(ikind)%array(ia)
1992 : ! fold atomic positions back into unit cell
1993 14224 : ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
1994 7112 : ria = ria - rcc
1995 1778 : atomic_kind => particle_set(iatom)%atomic_kind
1996 1778 : CALL get_atomic_kind(atomic_kind, kind_number=akind)
1997 1778 : CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
1998 1778 : rmom(1, 2) = rmom(1, 2) - charge
1999 10393 : DO l = 1, nm
2000 6297 : ix = indco(1, l + 1)
2001 6297 : iy = indco(2, l + 1)
2002 6297 : iz = indco(3, l + 1)
2003 6297 : dd = 1._dp
2004 6297 : IF (ix > 0) dd = dd*ria(1)**ix
2005 6297 : IF (iy > 0) dd = dd*ria(2)**iy
2006 6297 : IF (iz > 0) dd = dd*ria(3)**iz
2007 6297 : rmom(l + 1, 2) = rmom(l + 1, 2) - charge*dd
2008 8075 : CALL set_label(rlab(l + 1), ix, iy, iz)
2009 : END DO
2010 : END DO
2011 : END DO
2012 1082 : CALL para_env%sum(rmom(:, 2))
2013 18218 : rmom(:, :) = -rmom(:, :)
2014 5712 : rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
2015 :
2016 : ! magnetic moments
2017 1082 : IF (magnetic) THEN
2018 20 : NULLIFY (magmom)
2019 20 : CALL dbcsr_allocate_matrix_set(magmom, 3)
2020 80 : DO i = 1, SIZE(magmom)
2021 60 : CALL dbcsr_init_p(magmom(i)%matrix)
2022 : CALL dbcsr_create(magmom(i)%matrix, template=matrix_s(1)%matrix, &
2023 60 : matrix_type=dbcsr_type_antisymmetric)
2024 60 : CALL cp_dbcsr_alloc_block_from_nbl(magmom(i)%matrix, sab_orb)
2025 80 : CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
2026 : END DO
2027 :
2028 20 : CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
2029 :
2030 60 : ALLOCATE (mmom(SIZE(magmom)))
2031 20 : mmom(:) = 0.0_dp
2032 20 : IF (qs_env%run_rtp) THEN
2033 : ! get imaginary part of the density in rho_ao (the real part is not needed since the trace of the product
2034 : ! of a symmetric (REAL(rho_ao)) and an anti-symmetric (L_AO) matrix is zero)
2035 : ! There may be other cases, where the imaginary part of the density is relevant
2036 12 : NULLIFY (rho_ao)
2037 12 : CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2038 : END IF
2039 : ! if the density is purely real this is an expensive way to calculate zero
2040 80 : DO i = 1, SIZE(magmom)
2041 60 : strace = 0._dp
2042 120 : DO ispin = 1, dft_control%nspins
2043 60 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2044 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, magmom(i)%matrix, &
2045 60 : 0.0_dp, tmp_ao)
2046 60 : CALL dbcsr_trace(tmp_ao, trace)
2047 120 : strace = strace + trace
2048 : END DO
2049 80 : mmom(i) = strace
2050 : END DO
2051 :
2052 20 : CALL dbcsr_deallocate_matrix_set(magmom)
2053 : END IF
2054 :
2055 : ! velocity representations
2056 1082 : IF (my_velreprs) THEN
2057 120 : ALLOCATE (rmom_vel(nm))
2058 40 : rmom_vel = 0.0_dp
2059 :
2060 120 : DO order = 1, nmom
2061 40 : SELECT CASE (order)
2062 :
2063 : CASE (1) ! expectation value of momentum
2064 40 : NULLIFY (momentum)
2065 40 : CALL dbcsr_allocate_matrix_set(momentum, 3)
2066 160 : DO i = 1, 3
2067 120 : CALL dbcsr_init_p(momentum(i)%matrix)
2068 : CALL dbcsr_create(momentum(i)%matrix, template=matrix_s(1)%matrix, &
2069 120 : matrix_type=dbcsr_type_antisymmetric)
2070 120 : CALL cp_dbcsr_alloc_block_from_nbl(momentum(i)%matrix, sab_orb)
2071 160 : CALL dbcsr_set(momentum(i)%matrix, 0.0_dp)
2072 : END DO
2073 40 : CALL build_lin_mom_matrix(qs_env, momentum)
2074 :
2075 : ! imaginary part of the density for RTP, real part gives 0 since momentum is antisymmetric
2076 40 : IF (qs_env%run_rtp) THEN
2077 30 : NULLIFY (rho_ao)
2078 30 : CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2079 120 : DO idir = 1, SIZE(momentum)
2080 90 : strace = 0._dp
2081 180 : DO ispin = 1, dft_control%nspins
2082 90 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2083 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, momentum(idir)%matrix, &
2084 90 : 0.0_dp, tmp_ao)
2085 90 : CALL dbcsr_trace(tmp_ao, trace)
2086 180 : strace = strace + trace
2087 : END DO
2088 120 : rmom_vel(idir) = rmom_vel(idir) + strace
2089 : END DO
2090 : END IF
2091 :
2092 40 : CALL dbcsr_deallocate_matrix_set(momentum)
2093 :
2094 : CASE (2) ! expectation value of quadrupole moment in vel. reprs.
2095 40 : ALLOCATE (qupole_der(9)) ! will contain the expectation values of r_\alpha * d/d r_\beta
2096 40 : qupole_der = 0._dp
2097 :
2098 40 : NULLIFY (rho_ao)
2099 40 : CALL qs_rho_get(rho, rho_ao=rho_ao)
2100 :
2101 : ! Calculate expectation value over real part
2102 40 : trace = 0._dp
2103 160 : DO i = 1, 3
2104 520 : DO idir = 1, 3
2105 360 : strace = 0._dp
2106 720 : DO ispin = 1, dft_control%nspins
2107 360 : CALL dbcsr_set(tmp_ao, 0._dp)
2108 360 : CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2109 360 : CALL dbcsr_trace(tmp_ao, trace)
2110 720 : strace = strace + trace
2111 : END DO
2112 480 : qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2113 : END DO
2114 : END DO
2115 :
2116 40 : IF (qs_env%run_rtp) THEN
2117 30 : NULLIFY (rho_ao)
2118 30 : CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2119 :
2120 : ! Calculate expectation value over imaginary part
2121 30 : trace = 0._dp
2122 120 : DO i = 1, 3
2123 390 : DO idir = 1, 3
2124 270 : strace = 0._dp
2125 540 : DO ispin = 1, dft_control%nspins
2126 270 : CALL dbcsr_set(tmp_ao, 0._dp)
2127 270 : CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2128 270 : CALL dbcsr_trace(tmp_ao, trace)
2129 540 : strace = strace + trace
2130 : END DO
2131 360 : qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2132 : END DO
2133 : END DO
2134 : END IF
2135 :
2136 : ! calculate vel. reprs. of quadrupole moment from derivatives
2137 40 : rmom_vel(4) = -2*qupole_der(1) - rmom(1, 1)
2138 40 : rmom_vel(5) = -qupole_der(2) - qupole_der(4)
2139 40 : rmom_vel(6) = -qupole_der(3) - qupole_der(7)
2140 40 : rmom_vel(7) = -2*qupole_der(5) - rmom(1, 1)
2141 40 : rmom_vel(8) = -qupole_der(6) - qupole_der(8)
2142 40 : rmom_vel(9) = -2*qupole_der(9) - rmom(1, 1)
2143 :
2144 120 : DEALLOCATE (qupole_der)
2145 : CASE DEFAULT
2146 : END SELECT
2147 : END DO
2148 : END IF
2149 :
2150 1082 : IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2151 42 : CALL dbcsr_deallocate_matrix(tmp_ao)
2152 : END IF
2153 1082 : IF (my_velreprs .AND. (nmoments >= 2)) THEN
2154 40 : CALL dbcsr_deallocate_matrix_set(moments_der)
2155 : END IF
2156 :
2157 1082 : description = "[DIPOLE]"
2158 1082 : CALL cp_results_erase(results=results, description=description)
2159 : CALL put_results(results=results, description=description, &
2160 1082 : values=rmom(2:4, 3))
2161 :
2162 1082 : IF (magnetic .AND. my_velreprs) THEN
2163 18 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., mmom=mmom, rmom_vel=rmom_vel)
2164 1064 : ELSE IF (magnetic .AND. .NOT. my_velreprs) THEN
2165 2 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., mmom=mmom)
2166 1062 : ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2167 22 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., rmom_vel=rmom_vel)
2168 : ELSE
2169 1040 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE.)
2170 : END IF
2171 :
2172 1082 : IF (my_com_nl) THEN
2173 40 : IF (magnetic) THEN
2174 72 : mmom(:) = nlcom_rxrv(:)
2175 : END IF
2176 40 : IF (my_velreprs .AND. (nmom >= 1)) THEN
2177 40 : DEALLOCATE (rmom_vel)
2178 40 : ALLOCATE (rmom_vel(21))
2179 160 : rmom_vel(1:3) = nlcom_rv
2180 : END IF
2181 40 : IF (my_velreprs .AND. (nmom >= 2)) THEN
2182 280 : rmom_vel(4:9) = nlcom_rrv
2183 280 : rmom_vel(10:15) = nlcom_rvr
2184 280 : rmom_vel(16:21) = nlcom_rrv_vrr
2185 : END IF
2186 40 : IF (magnetic .AND. .NOT. my_velreprs) THEN
2187 0 : CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom)
2188 40 : ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2189 22 : CALL print_moments_nl(unit_number, nmom, rlab, rmom_vel=rmom_vel)
2190 18 : ELSE IF (my_velreprs .AND. magnetic) THEN
2191 18 : CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom, rmom_vel=rmom_vel)
2192 : END IF
2193 :
2194 : END IF
2195 :
2196 : IF (my_com_nl) THEN
2197 40 : IF (nmom >= 1 .AND. my_velreprs) DEALLOCATE (nlcom_rv)
2198 40 : IF (nmom >= 2 .AND. my_velreprs) THEN
2199 40 : DEALLOCATE (nlcom_rrv)
2200 40 : DEALLOCATE (nlcom_rvr)
2201 40 : DEALLOCATE (nlcom_rrv_vrr)
2202 : END IF
2203 40 : IF (magnetic) DEALLOCATE (nlcom_rxrv)
2204 : END IF
2205 :
2206 1082 : DEALLOCATE (rmom)
2207 1082 : DEALLOCATE (rlab)
2208 1082 : IF (magnetic) THEN
2209 20 : DEALLOCATE (mmom)
2210 : END IF
2211 1082 : IF (my_velreprs) THEN
2212 40 : DEALLOCATE (rmom_vel)
2213 : END IF
2214 :
2215 1082 : CALL timestop(handle)
2216 :
2217 2164 : END SUBROUTINE qs_moment_locop
2218 :
2219 : ! **************************************************************************************************
2220 : !> \brief ...
2221 : !> \param label ...
2222 : !> \param ix ...
2223 : !> \param iy ...
2224 : !> \param iz ...
2225 : ! **************************************************************************************************
2226 8415 : SUBROUTINE set_label(label, ix, iy, iz)
2227 : CHARACTER(LEN=*), INTENT(OUT) :: label
2228 : INTEGER, INTENT(IN) :: ix, iy, iz
2229 :
2230 : INTEGER :: i
2231 :
2232 8415 : label = ""
2233 11581 : DO i = 1, ix
2234 11581 : WRITE (label(i:), "(A1)") "X"
2235 : END DO
2236 11581 : DO i = ix + 1, ix + iy
2237 11581 : WRITE (label(i:), "(A1)") "Y"
2238 : END DO
2239 11581 : DO i = ix + iy + 1, ix + iy + iz
2240 11581 : WRITE (label(i:), "(A1)") "Z"
2241 : END DO
2242 :
2243 8415 : END SUBROUTINE set_label
2244 :
2245 : ! **************************************************************************************************
2246 : !> \brief ...
2247 : !> \param unit_number ...
2248 : !> \param nmom ...
2249 : !> \param rmom ...
2250 : !> \param rlab ...
2251 : !> \param rcc ...
2252 : !> \param cell ...
2253 : !> \param periodic ...
2254 : !> \param mmom ...
2255 : !> \param rmom_vel ...
2256 : ! **************************************************************************************************
2257 1578 : SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmom, rmom_vel)
2258 : INTEGER, INTENT(IN) :: unit_number, nmom
2259 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: rmom
2260 : CHARACTER(LEN=8), DIMENSION(:) :: rlab
2261 : REAL(dp), DIMENSION(3), INTENT(IN) :: rcc
2262 : TYPE(cell_type), POINTER :: cell
2263 : LOGICAL :: periodic
2264 : REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2265 :
2266 : INTEGER :: i, i0, i1, j, l
2267 : REAL(dp) :: dd
2268 :
2269 1578 : IF (unit_number > 0) THEN
2270 2438 : DO l = 0, nmom
2271 805 : SELECT CASE (l)
2272 : CASE (0)
2273 805 : WRITE (unit_number, "(T3,A,T33,3F16.8)") "Reference Point [Bohr]", rcc
2274 805 : WRITE (unit_number, "(T3,A)") "Charges"
2275 : WRITE (unit_number, "(T5,A,T18,F14.8,T36,A,T42,F14.8,T60,A,T67,F14.8)") &
2276 805 : "Electronic=", rmom(1, 1), "Core=", rmom(1, 2), "Total=", rmom(1, 3)
2277 : CASE (1)
2278 805 : IF (periodic) THEN
2279 246 : WRITE (unit_number, "(T3,A)") "Dipole vectors are based on the periodic (Berry phase) operator."
2280 246 : WRITE (unit_number, "(T3,A)") "They are defined modulo integer multiples of the cell matrix [Debye]."
2281 984 : WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[X] [", cell%hmat(1, :)*debye, "] [i]"
2282 984 : WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Y]=[", cell%hmat(2, :)*debye, "]*[j]"
2283 984 : WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Z] [", cell%hmat(3, :)*debye, "] [k]"
2284 : ELSE
2285 559 : WRITE (unit_number, "(T3,A)") "Dipoles are based on the traditional operator."
2286 : END IF
2287 3220 : dd = SQRT(SUM(rmom(2:4, 3)**2))*debye
2288 805 : WRITE (unit_number, "(T3,A)") "Dipole moment [Debye]"
2289 : WRITE (unit_number, "(T5,3(A,A,E15.7,1X),T60,A,T68,F13.7)") &
2290 3220 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye, i=2, 4), "Total=", dd
2291 : CASE (2)
2292 21 : WRITE (unit_number, "(T3,A)") "Quadrupole moment [Debye*Angstrom]"
2293 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2294 84 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=5, 7)
2295 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2296 84 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=8, 10)
2297 : CASE (3)
2298 1 : WRITE (unit_number, "(T3,A)") "Octapole moment [Debye*Angstrom**2]"
2299 : WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2300 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=11, 14)
2301 : WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2302 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=15, 18)
2303 : WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2304 3 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=19, 20)
2305 : CASE (4)
2306 1 : WRITE (unit_number, "(T3,A)") "Hexadecapole moment [Debye*Angstrom**3]"
2307 : WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2308 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=21, 24)
2309 : WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2310 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=25, 28)
2311 : WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2312 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=29, 32)
2313 : WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2314 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=32, 35)
2315 : CASE DEFAULT
2316 0 : WRITE (unit_number, "(T3,A,A,I2)") "Higher moment [Debye*Angstrom**(L-1)]", &
2317 0 : " L=", l
2318 0 : i0 = (6 + 11*(l - 1) + 6*(l - 1)**2 + (l - 1)**3)/6
2319 0 : i1 = (6 + 11*l + 6*l**2 + l**3)/6 - 1
2320 0 : dd = debye/(bohr)**(l - 1)
2321 1633 : DO i = i0, i1, 3
2322 : WRITE (unit_number, "(T18,3(A,A,F14.8,4X))") &
2323 0 : (TRIM(rlab(j + 1)), "=", rmom(j + 1, 3)*dd, j=i, MIN(i1, i + 2))
2324 : END DO
2325 : END SELECT
2326 : END DO
2327 805 : IF (PRESENT(mmom)) THEN
2328 28 : IF (nmom >= 1) THEN
2329 112 : dd = SQRT(SUM(mmom(1:3)**2))
2330 28 : WRITE (unit_number, "(T3,A)") "Orbital angular momentum [a. u.]"
2331 : WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2332 112 : (TRIM(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2333 : END IF
2334 : END IF
2335 805 : IF (PRESENT(rmom_vel)) THEN
2336 96 : DO l = 1, nmom
2337 38 : SELECT CASE (l)
2338 : CASE (1)
2339 152 : dd = SQRT(SUM(rmom_vel(1:3)**2))
2340 38 : WRITE (unit_number, "(T3,A)") "Expectation value of momentum operator [a. u.]"
2341 : WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2342 152 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2343 : CASE (2)
2344 20 : WRITE (unit_number, "(T3,A)") "Expectation value of quadrupole operator in vel. repr. [a. u.]"
2345 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2346 80 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2347 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2348 138 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2349 : CASE DEFAULT
2350 : END SELECT
2351 : END DO
2352 : END IF
2353 : END IF
2354 :
2355 1578 : END SUBROUTINE print_moments
2356 :
2357 : ! **************************************************************************************************
2358 : !> \brief ...
2359 : !> \param unit_number ...
2360 : !> \param nmom ...
2361 : !> \param rlab ...
2362 : !> \param mmom ...
2363 : !> \param rmom_vel ...
2364 : ! **************************************************************************************************
2365 58 : SUBROUTINE print_moments_nl(unit_number, nmom, rlab, mmom, rmom_vel)
2366 : INTEGER, INTENT(IN) :: unit_number, nmom
2367 : CHARACTER(LEN=8), DIMENSION(:) :: rlab
2368 : REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2369 :
2370 : INTEGER :: i, l
2371 : REAL(dp) :: dd
2372 :
2373 58 : IF (unit_number > 0) THEN
2374 38 : IF (PRESENT(mmom)) THEN
2375 27 : IF (nmom >= 1) THEN
2376 108 : dd = SQRT(SUM(mmom(1:3)**2))
2377 27 : WRITE (unit_number, "(T3,A)") "Expectation value of rx[r,V_nl] [a. u.]"
2378 : WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2379 108 : (TRIM(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2380 : END IF
2381 : END IF
2382 38 : IF (PRESENT(rmom_vel)) THEN
2383 96 : DO l = 1, nmom
2384 38 : SELECT CASE (l)
2385 : CASE (1)
2386 152 : dd = SQRT(SUM(rmom_vel(1:3)**2))
2387 38 : WRITE (unit_number, "(T3,A)") "Expectation value of [r,V_nl] [a. u.]"
2388 : WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2389 152 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2390 : CASE (2)
2391 20 : WRITE (unit_number, "(T3,A)") "Expectation value of [rr,V_nl] [a. u.]"
2392 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2393 80 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2394 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2395 80 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2396 20 : WRITE (unit_number, "(T3,A)") "Expectation value of r x V_nl x r [a. u.]"
2397 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2398 80 : (TRIM(rlab(i + 1 - 6)), "=", rmom_vel(i), i=10, 12)
2399 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2400 80 : (TRIM(rlab(i + 1 - 6)), "=", rmom_vel(i), i=13, 15)
2401 20 : WRITE (unit_number, "(T3,A)") "Expectation value of r x r x V_nl + V_nl x r x r [a. u.]"
2402 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2403 80 : (TRIM(rlab(i + 1 - 12)), "=", rmom_vel(i), i=16, 18)
2404 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2405 138 : (TRIM(rlab(i + 1 - 12)), "=", rmom_vel(i), i=19, 21)
2406 : CASE DEFAULT
2407 : END SELECT
2408 : END DO
2409 : END IF
2410 : END IF
2411 :
2412 58 : END SUBROUTINE print_moments_nl
2413 :
2414 : ! **************************************************************************************************
2415 : !> \brief Calculate the expectation value of operators related to non-local potential:
2416 : !> [r, Vnl], noted rv
2417 : !> r x [r,Vnl], noted rxrv
2418 : !> [rr,Vnl], noted rrv
2419 : !> r x Vnl x r, noted rvr
2420 : !> r x r x Vnl + Vnl x r x r, noted rrv_vrr
2421 : !> Note that the 3 first operator are commutator while the 2 last
2422 : !> are not. For reading clarity the same notation is used for all 5
2423 : !> operators.
2424 : !> \param qs_env ...
2425 : !> \param nlcom_rv ...
2426 : !> \param nlcom_rxrv ...
2427 : !> \param nlcom_rrv ...
2428 : !> \param nlcom_rvr ...
2429 : !> \param nlcom_rrv_vrr ...
2430 : !> \param ref_point ...
2431 : ! **************************************************************************************************
2432 76 : SUBROUTINE calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, &
2433 : nlcom_rrv_vrr, ref_point)
2434 :
2435 : TYPE(qs_environment_type), POINTER :: qs_env
2436 : REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: nlcom_rv, nlcom_rxrv, nlcom_rrv, &
2437 : nlcom_rvr, nlcom_rrv_vrr
2438 : REAL(dp), DIMENSION(3) :: ref_point
2439 :
2440 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_commutator_nl_terms'
2441 :
2442 : INTEGER :: handle, ind, ispin
2443 : LOGICAL :: calc_rrv, calc_rrv_vrr, calc_rv, &
2444 : calc_rvr, calc_rxrv
2445 : REAL(dp) :: eps_ppnl, strace, trace
2446 : TYPE(cell_type), POINTER :: cell
2447 76 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rrv, matrix_rrv_vrr, matrix_rv, &
2448 76 : matrix_rvr, matrix_rxrv, matrix_s, &
2449 76 : rho_ao
2450 : TYPE(dbcsr_type), POINTER :: tmp_ao
2451 : TYPE(dft_control_type), POINTER :: dft_control
2452 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2453 76 : POINTER :: sab_all, sab_orb, sap_ppnl
2454 76 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2455 76 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2456 : TYPE(qs_rho_type), POINTER :: rho
2457 :
2458 76 : CALL timeset(routineN, handle)
2459 :
2460 76 : calc_rv = .FALSE.
2461 76 : calc_rxrv = .FALSE.
2462 76 : calc_rrv = .FALSE.
2463 76 : calc_rvr = .FALSE.
2464 76 : calc_rrv_vrr = .FALSE.
2465 :
2466 : ! rv, rxrv and rrv are commutator matrices: anti-symmetric.
2467 : ! The real part of the density matrix rho_ao is symmetric so that
2468 : ! the expectation value of real density matrix is zero. Hence, if
2469 : ! the density matrix is real, no need to compute these quantities.
2470 : ! This is not the case for rvr and rrv_vrr which are symmetric.
2471 :
2472 76 : IF (ALLOCATED(nlcom_rv)) THEN
2473 76 : nlcom_rv(:) = 0._dp
2474 76 : IF (qs_env%run_rtp) calc_rv = .TRUE.
2475 : END IF
2476 76 : IF (ALLOCATED(nlcom_rxrv)) THEN
2477 54 : nlcom_rxrv(:) = 0._dp
2478 54 : IF (qs_env%run_rtp) calc_rxrv = .TRUE.
2479 : END IF
2480 76 : IF (ALLOCATED(nlcom_rrv)) THEN
2481 40 : nlcom_rrv(:) = 0._dp
2482 40 : IF (qs_env%run_rtp) calc_rrv = .TRUE.
2483 : END IF
2484 76 : IF (ALLOCATED(nlcom_rvr)) THEN
2485 40 : nlcom_rvr(:) = 0._dp
2486 40 : calc_rvr = .TRUE.
2487 : END IF
2488 76 : IF (ALLOCATED(nlcom_rrv_vrr)) THEN
2489 40 : nlcom_rrv_vrr(:) = 0._dp
2490 40 : calc_rrv_vrr = .TRUE.
2491 : END IF
2492 :
2493 76 : IF (.NOT. (calc_rv .OR. calc_rrv .OR. calc_rxrv .OR. calc_rvr .OR. calc_rrv_vrr)) THEN
2494 12 : CALL timestop(handle)
2495 12 : RETURN
2496 : END IF
2497 :
2498 64 : NULLIFY (cell, matrix_s, particle_set, qs_kind_set, rho, sab_all, sab_orb, sap_ppnl)
2499 : CALL get_qs_env(qs_env, &
2500 : cell=cell, &
2501 : dft_control=dft_control, &
2502 : matrix_s=matrix_s, &
2503 : particle_set=particle_set, &
2504 : qs_kind_set=qs_kind_set, &
2505 : rho=rho, &
2506 : sab_orb=sab_orb, &
2507 : sab_all=sab_all, &
2508 64 : sap_ppnl=sap_ppnl)
2509 :
2510 64 : eps_ppnl = dft_control%qs_control%eps_ppnl
2511 :
2512 : ! Allocate storage
2513 64 : NULLIFY (matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr)
2514 64 : IF (calc_rv) THEN
2515 54 : CALL dbcsr_allocate_matrix_set(matrix_rv, 3)
2516 216 : DO ind = 1, 3
2517 162 : CALL dbcsr_init_p(matrix_rv(ind)%matrix)
2518 : CALL dbcsr_create(matrix_rv(ind)%matrix, template=matrix_s(1)%matrix, &
2519 162 : matrix_type=dbcsr_type_antisymmetric)
2520 162 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(ind)%matrix, sab_orb)
2521 216 : CALL dbcsr_set(matrix_rv(ind)%matrix, 0._dp)
2522 : END DO
2523 : END IF
2524 :
2525 64 : IF (calc_rxrv) THEN
2526 36 : CALL dbcsr_allocate_matrix_set(matrix_rxrv, 3)
2527 144 : DO ind = 1, 3
2528 108 : CALL dbcsr_init_p(matrix_rxrv(ind)%matrix)
2529 : CALL dbcsr_create(matrix_rxrv(ind)%matrix, template=matrix_s(1)%matrix, &
2530 108 : matrix_type=dbcsr_type_antisymmetric)
2531 108 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_rxrv(ind)%matrix, sab_orb)
2532 144 : CALL dbcsr_set(matrix_rxrv(ind)%matrix, 0._dp)
2533 : END DO
2534 : END IF
2535 :
2536 64 : IF (calc_rrv) THEN
2537 30 : CALL dbcsr_allocate_matrix_set(matrix_rrv, 6)
2538 210 : DO ind = 1, 6
2539 180 : CALL dbcsr_init_p(matrix_rrv(ind)%matrix)
2540 : CALL dbcsr_create(matrix_rrv(ind)%matrix, template=matrix_s(1)%matrix, &
2541 180 : matrix_type=dbcsr_type_antisymmetric)
2542 180 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv(ind)%matrix, sab_orb)
2543 210 : CALL dbcsr_set(matrix_rrv(ind)%matrix, 0._dp)
2544 : END DO
2545 : END IF
2546 :
2547 64 : IF (calc_rvr) THEN
2548 40 : CALL dbcsr_allocate_matrix_set(matrix_rvr, 6)
2549 280 : DO ind = 1, 6
2550 240 : CALL dbcsr_init_p(matrix_rvr(ind)%matrix)
2551 : CALL dbcsr_create(matrix_rvr(ind)%matrix, template=matrix_s(1)%matrix, &
2552 240 : matrix_type=dbcsr_type_symmetric)
2553 240 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_rvr(ind)%matrix, sab_orb)
2554 280 : CALL dbcsr_set(matrix_rvr(ind)%matrix, 0._dp)
2555 : END DO
2556 : END IF
2557 64 : IF (calc_rrv_vrr) THEN
2558 40 : CALL dbcsr_allocate_matrix_set(matrix_rrv_vrr, 6)
2559 280 : DO ind = 1, 6
2560 240 : CALL dbcsr_init_p(matrix_rrv_vrr(ind)%matrix)
2561 : CALL dbcsr_create(matrix_rrv_vrr(ind)%matrix, template=matrix_s(1)%matrix, &
2562 240 : matrix_type=dbcsr_type_symmetric)
2563 240 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv_vrr(ind)%matrix, sab_orb)
2564 280 : CALL dbcsr_set(matrix_rrv_vrr(ind)%matrix, 0._dp)
2565 : END DO
2566 : END IF
2567 :
2568 : ! calculate evaluation of operators in AO basis set
2569 : CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
2570 : matrix_rxrv=matrix_rxrv, matrix_rrv=matrix_rrv, matrix_rvr=matrix_rvr, &
2571 64 : matrix_rrv_vrr=matrix_rrv_vrr, ref_point=ref_point)
2572 :
2573 : ! Calculate expectation values
2574 : ! Real part
2575 64 : NULLIFY (tmp_ao)
2576 64 : CALL dbcsr_init_p(tmp_ao)
2577 64 : CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
2578 64 : CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
2579 64 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2580 :
2581 64 : IF (calc_rvr .OR. calc_rrv_vrr) THEN
2582 40 : NULLIFY (rho_ao)
2583 40 : CALL qs_rho_get(rho, rho_ao=rho_ao)
2584 :
2585 40 : IF (calc_rvr) THEN
2586 : trace = 0._dp
2587 280 : DO ind = 1, SIZE(matrix_rvr)
2588 240 : strace = 0._dp
2589 480 : DO ispin = 1, dft_control%nspins
2590 240 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2591 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rvr(ind)%matrix, &
2592 240 : 0.0_dp, tmp_ao)
2593 240 : CALL dbcsr_trace(tmp_ao, trace)
2594 480 : strace = strace + trace
2595 : END DO
2596 280 : nlcom_rvr(ind) = nlcom_rvr(ind) + strace
2597 : END DO
2598 : END IF
2599 :
2600 40 : IF (calc_rrv_vrr) THEN
2601 : trace = 0._dp
2602 280 : DO ind = 1, SIZE(matrix_rrv_vrr)
2603 240 : strace = 0._dp
2604 480 : DO ispin = 1, dft_control%nspins
2605 240 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2606 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv_vrr(ind)%matrix, &
2607 240 : 0.0_dp, tmp_ao)
2608 240 : CALL dbcsr_trace(tmp_ao, trace)
2609 480 : strace = strace + trace
2610 : END DO
2611 280 : nlcom_rrv_vrr(ind) = nlcom_rrv_vrr(ind) + strace
2612 : END DO
2613 : END IF
2614 : END IF
2615 :
2616 : ! imagninary part of the density matrix
2617 64 : NULLIFY (rho_ao)
2618 64 : CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2619 :
2620 64 : IF (calc_rv) THEN
2621 : trace = 0._dp
2622 216 : DO ind = 1, SIZE(matrix_rv)
2623 162 : strace = 0._dp
2624 324 : DO ispin = 1, dft_control%nspins
2625 162 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2626 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rv(ind)%matrix, &
2627 162 : 0.0_dp, tmp_ao)
2628 162 : CALL dbcsr_trace(tmp_ao, trace)
2629 324 : strace = strace + trace
2630 : END DO
2631 216 : nlcom_rv(ind) = nlcom_rv(ind) + strace
2632 : END DO
2633 : END IF
2634 :
2635 64 : IF (calc_rrv) THEN
2636 : trace = 0._dp
2637 210 : DO ind = 1, SIZE(matrix_rrv)
2638 180 : strace = 0._dp
2639 360 : DO ispin = 1, dft_control%nspins
2640 180 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2641 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv(ind)%matrix, &
2642 180 : 0.0_dp, tmp_ao)
2643 180 : CALL dbcsr_trace(tmp_ao, trace)
2644 360 : strace = strace + trace
2645 : END DO
2646 210 : nlcom_rrv(ind) = nlcom_rrv(ind) + strace
2647 : END DO
2648 : END IF
2649 :
2650 64 : IF (calc_rxrv) THEN
2651 : trace = 0._dp
2652 144 : DO ind = 1, SIZE(matrix_rxrv)
2653 108 : strace = 0._dp
2654 216 : DO ispin = 1, dft_control%nspins
2655 108 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2656 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rxrv(ind)%matrix, &
2657 108 : 0.0_dp, tmp_ao)
2658 108 : CALL dbcsr_trace(tmp_ao, trace)
2659 216 : strace = strace + trace
2660 : END DO
2661 144 : nlcom_rxrv(ind) = nlcom_rxrv(ind) + strace
2662 : END DO
2663 : END IF
2664 64 : CALL dbcsr_deallocate_matrix(tmp_ao)
2665 64 : IF (calc_rv) CALL dbcsr_deallocate_matrix_set(matrix_rv)
2666 64 : IF (calc_rxrv) CALL dbcsr_deallocate_matrix_set(matrix_rxrv)
2667 64 : IF (calc_rrv) CALL dbcsr_deallocate_matrix_set(matrix_rrv)
2668 64 : IF (calc_rvr) CALL dbcsr_deallocate_matrix_set(matrix_rvr)
2669 64 : IF (calc_rrv_vrr) CALL dbcsr_deallocate_matrix_set(matrix_rrv_vrr)
2670 :
2671 64 : CALL timestop(handle)
2672 76 : END SUBROUTINE calculate_commutator_nl_terms
2673 :
2674 : ! **************************************************************************************************
2675 : !> \brief Get list of kpoints from input to compute dipole moment elements
2676 : !> \param qs_env ...
2677 : !> \param xkp ...
2678 : !> \param special_pnts ...
2679 : !> \author Shridhar Shanbhag
2680 : ! **************************************************************************************************
2681 10 : SUBROUTINE get_xkp_for_dipole_calc(qs_env, xkp, special_pnts)
2682 : TYPE(qs_environment_type), POINTER :: qs_env
2683 : TYPE(section_vals_type), POINTER :: kpnts, kpset
2684 : REAL(KIND=dp), DIMENSION(3, 3) :: cart_hmat, hmat
2685 :
2686 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_xkp_for_dipole_calc'
2687 :
2688 : CHARACTER(LEN=default_string_length) :: ustr
2689 : TYPE(kpoint_type), POINTER :: kpoint_work
2690 : TYPE(cell_type), POINTER :: cell
2691 : CHARACTER(LEN=default_string_length), &
2692 10 : DIMENSION(:), POINTER :: strptr
2693 : CHARACTER(LEN=default_string_length), &
2694 10 : DIMENSION(:), POINTER :: special_pnts, spname
2695 : CHARACTER(LEN=max_line_length) :: error_message
2696 : INTEGER :: handle, i, ik, ikk, ip, &
2697 : n_ptr, npline, nkp
2698 : LOGICAL :: explicit_kpnts, explicit_kpset
2699 10 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: kspecial, xkp
2700 : REAL(KIND=dp), DIMENSION(3) :: kpptr
2701 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2702 :
2703 10 : CALL timeset(routineN, handle)
2704 10 : kpset => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINT_SET")
2705 10 : kpnts => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINTS")
2706 10 : CALL section_vals_get(kpset, explicit=explicit_kpset)
2707 10 : CALL section_vals_get(kpnts, explicit=explicit_kpnts)
2708 10 : IF (explicit_kpset .AND. explicit_kpnts) then
2709 0 : CPABORT("Both KPOINT_SET and KPOINTS present in MOMENTS section")
2710 : end if
2711 :
2712 10 : IF (explicit_kpset) THEN
2713 4 : CALL get_qs_env(qs_env, cell=cell)
2714 4 : CALL get_cell(cell, h=hmat)
2715 4 : cart_hmat(:, :) = hmat(:, :)
2716 4 : IF (cell%input_cell_canonicalized) cart_hmat(:, :) = cell%input_hmat(:, :)
2717 4 : CALL section_vals_val_get(kpset, "NPOINTS", i_val=npline)
2718 4 : CALL section_vals_val_get(kpset, "UNITS", c_val=ustr)
2719 4 : CALL uppercase(ustr)
2720 4 : CALL section_vals_val_get(kpset, "SPECIAL_POINT", n_rep_val=n_ptr)
2721 4 : CPASSERT(n_ptr > 0)
2722 12 : ALLOCATE (kspecial(3, n_ptr))
2723 12 : ALLOCATE (spname(n_ptr))
2724 8 : DO ip = 1, n_ptr
2725 4 : CALL section_vals_val_get(kpset, "SPECIAL_POINT", i_rep_val=ip, c_vals=strptr)
2726 4 : IF (SIZE(strptr(:), 1) == 4) THEN
2727 2 : spname(ip) = strptr(1)
2728 8 : DO i = 1, 3
2729 6 : CALL read_float_object(strptr(i + 1), kpptr(i), error_message)
2730 8 : IF (LEN_TRIM(error_message) > 0) CPABORT(TRIM(error_message))
2731 : END DO
2732 2 : ELSE IF (SIZE(strptr(:), 1) == 3) THEN
2733 2 : spname(ip) = "not specified"
2734 8 : DO i = 1, 3
2735 6 : CALL read_float_object(strptr(i), kpptr(i), error_message)
2736 8 : IF (LEN_TRIM(error_message) > 0) CPABORT(TRIM(error_message))
2737 : END DO
2738 : ELSE
2739 0 : CPABORT("Input SPECIAL_POINT invalid")
2740 : END IF
2741 4 : SELECT CASE (ustr)
2742 : CASE ("B_VECTOR")
2743 16 : kspecial(1:3, ip) = kpptr(1:3)
2744 : CASE ("CART_ANGSTROM")
2745 : kspecial(1:3, ip) = (kpptr(1)*cart_hmat(1, 1:3) + &
2746 : kpptr(2)*cart_hmat(2, 1:3) + &
2747 0 : kpptr(3)*cart_hmat(3, 1:3))/twopi*angstrom
2748 : CASE ("CART_BOHR")
2749 : kspecial(1:3, ip) = (kpptr(1)*cart_hmat(1, 1:3) + &
2750 : kpptr(2)*cart_hmat(2, 1:3) + &
2751 0 : kpptr(3)*cart_hmat(3, 1:3))/twopi
2752 : CASE DEFAULT
2753 4 : CPABORT("Unknown unit <"//TRIM(ustr)//"> specified for k-point definition")
2754 : END SELECT
2755 : END DO
2756 4 : nkp = (n_ptr - 1)*npline + 1
2757 4 : CPASSERT(nkp >= 1)
2758 :
2759 : ! Initialize environment and calculate MOs
2760 12 : ALLOCATE (xkp(3, nkp))
2761 12 : ALLOCATE (special_pnts(nkp))
2762 8 : special_pnts(:) = ""
2763 16 : xkp(1:3, 1) = kspecial(1:3, 1)
2764 4 : ikk = 1
2765 4 : special_pnts(ikk) = spname(1)
2766 4 : DO ik = 2, n_ptr
2767 0 : DO ip = 1, npline
2768 0 : ikk = ikk + 1
2769 : xkp(1:3, ikk) = kspecial(1:3, ik - 1) + &
2770 : REAL(ip, KIND=dp)/REAL(npline, KIND=dp)* &
2771 0 : (kspecial(1:3, ik) - kspecial(1:3, ik - 1))
2772 : END DO
2773 4 : special_pnts(ikk) = spname(ik)
2774 : END DO
2775 12 : DEALLOCATE (spname, kspecial)
2776 6 : ELSE IF (explicit_kpnts) THEN
2777 2 : CALL get_qs_env(qs_env, particle_set=particle_set, cell=cell)
2778 2 : CALL get_cell(cell, h=hmat)
2779 2 : NULLIFY (kpoint_work)
2780 2 : CALL kpoint_create(kpoint_work)
2781 2 : CALL read_kpoint_section(kpoint_work, kpnts, hmat, cell)
2782 2 : CALL kpoint_initialize(kpoint_work, particle_set, cell)
2783 2 : nkp = kpoint_work%nkp
2784 6 : ALLOCATE (xkp(3, nkp))
2785 6 : ALLOCATE (special_pnts(nkp))
2786 4 : special_pnts(:) = ""
2787 10 : xkp(1:3, :) = kpoint_work%xkp(1:3, :)
2788 2 : CALL kpoint_release(kpoint_work)
2789 : ELSE
2790 : ! use k-point mesh from DFT calculation
2791 4 : CALL get_qs_env(qs_env, kpoints=kpoint_work)
2792 4 : nkp = kpoint_work%nkp
2793 4 : nkp = kpoint_work%nkp
2794 12 : ALLOCATE (xkp(3, nkp))
2795 12 : ALLOCATE (special_pnts(nkp))
2796 252 : special_pnts(:) = ""
2797 996 : xkp(1:3, :) = kpoint_work%xkp(1:3, :)
2798 : END IF
2799 10 : CALL timestop(handle)
2800 :
2801 10 : END SUBROUTINE get_xkp_for_dipole_calc
2802 : ! **************************************************************************************************
2803 : !> \brief Calculate local moment matrix for a periodic system for all image cells
2804 : !> \param qs_env ...
2805 : !> \param moments_rs_img ...
2806 : !> \param rcc ...
2807 : !> \author Shridhar Shanbhag
2808 : ! **************************************************************************************************
2809 8 : SUBROUTINE build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc)
2810 :
2811 : TYPE(qs_environment_type), POINTER :: qs_env
2812 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_rs_img
2813 : REAL(KIND=dp), DIMENSION(3), OPTIONAL :: rcc
2814 :
2815 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_local_moment_matrix_rs_img'
2816 :
2817 : INTEGER :: handle, i_dir, iatom, ic, ikind, iset, j, jatom, jkind, jset, &
2818 : ldsa, ldsb, ldwork, ncoa, ncob, nimg, nkind, nseta, nsetb, nsize, sgfa, sgfb
2819 : INTEGER, DIMENSION(3) :: icell
2820 8 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
2821 8 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2822 : LOGICAL :: found
2823 : REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
2824 8 : REAL(dp), DIMENSION(:, :), POINTER :: dblock, work
2825 8 : REAL(dp), DIMENSION(:, :, :), POINTER :: dipab
2826 : REAL(KIND=dp) :: dab
2827 : TYPE(cell_type), POINTER :: cell
2828 8 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
2829 : TYPE(dft_control_type), POINTER :: dft_control
2830 8 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
2831 : TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
2832 : TYPE(kpoint_type), POINTER :: kpoints_all
2833 : TYPE(mp_para_env_type), POINTER :: para_env
2834 : TYPE(neighbor_list_iterator_p_type), &
2835 8 : DIMENSION(:), POINTER :: nl_iterator
2836 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2837 8 : POINTER :: sab_all
2838 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2839 8 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2840 : TYPE(qs_kind_type), POINTER :: qs_kind
2841 :
2842 8 : CALL timeset(routineN, handle)
2843 :
2844 : CALL get_qs_env(qs_env=qs_env, &
2845 : dft_control=dft_control, &
2846 : qs_kind_set=qs_kind_set, &
2847 : matrix_ks_kp=matrix_ks_kp, &
2848 : particle_set=particle_set, &
2849 : cell=cell, &
2850 : para_env=para_env, &
2851 8 : sab_all=sab_all)
2852 :
2853 8 : NULLIFY (kpoints_all)
2854 8 : CALL kpoint_create(kpoints_all)
2855 8 : CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, nimg)
2856 :
2857 8 : nkind = SIZE(qs_kind_set)
2858 36 : ALLOCATE (basis_set_list(nkind))
2859 20 : DO ikind = 1, nkind
2860 12 : qs_kind => qs_kind_set(ikind)
2861 12 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set)
2862 20 : IF (ASSOCIATED(basis_set)) THEN
2863 12 : basis_set_list(ikind)%gto_basis_set => basis_set
2864 : ELSE
2865 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
2866 : END IF
2867 : END DO
2868 :
2869 8 : rc(:) = 0._dp
2870 8 : IF (PRESENT(rcc)) rc(:) = rcc(:)
2871 :
2872 8 : CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_list)
2873 8 : CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
2874 8 : nsize = SIZE(index_to_cell, 2)
2875 8 : CPASSERT(SIZE(moments_rs_img, 2) == nsize)
2876 32 : DO i_dir = 1, 3
2877 2360 : DO j = 1, nsize
2878 2328 : ALLOCATE (moments_rs_img(i_dir, j)%matrix)
2879 : CALL dbcsr_create(matrix=moments_rs_img(i_dir, j)%matrix, &
2880 : template=matrix_ks_kp(1, 1)%matrix, &
2881 : matrix_type=dbcsr_type_no_symmetry, &
2882 2328 : name="DIPMAT")
2883 2328 : CALL cp_dbcsr_alloc_block_from_nbl(moments_rs_img(i_dir, j)%matrix, sab_all)
2884 2352 : CALL dbcsr_set(moments_rs_img(i_dir, j)%matrix, 0.0_dp)
2885 : END DO
2886 : END DO
2887 :
2888 8 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
2889 40 : ALLOCATE (dipab(ldwork, ldwork, 3))
2890 32 : ALLOCATE (work(ldwork, ldwork))
2891 :
2892 8 : CALL neighbor_list_iterator_create(nl_iterator, sab_all)
2893 2100 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2894 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2895 2092 : iatom=iatom, jatom=jatom, r=rab, cell=icell)
2896 :
2897 2092 : basis_set_a => basis_set_list(ikind)%gto_basis_set
2898 2092 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
2899 2092 : basis_set_b => basis_set_list(jkind)%gto_basis_set
2900 2092 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
2901 : ASSOCIATE ( &
2902 : ! basis ikind
2903 : first_sgfa => basis_set_a%first_sgf, &
2904 : la_max => basis_set_a%lmax, &
2905 : la_min => basis_set_a%lmin, &
2906 : npgfa => basis_set_a%npgf, &
2907 : nsgfa => basis_set_a%nsgf_set, &
2908 : rpgfa => basis_set_a%pgf_radius, &
2909 : set_radius_a => basis_set_a%set_radius, &
2910 : sphi_a => basis_set_a%sphi, &
2911 : zeta => basis_set_a%zet, &
2912 : ! basis jkind, &
2913 : first_sgfb => basis_set_b%first_sgf, &
2914 : lb_max => basis_set_b%lmax, &
2915 : lb_min => basis_set_b%lmin, &
2916 : npgfb => basis_set_b%npgf, &
2917 : nsgfb => basis_set_b%nsgf_set, &
2918 : rpgfb => basis_set_b%pgf_radius, &
2919 : set_radius_b => basis_set_b%set_radius, &
2920 : sphi_b => basis_set_b%sphi, &
2921 : zetb => basis_set_b%zet)
2922 :
2923 2092 : nseta = basis_set_a%nset
2924 2092 : nsetb = basis_set_b%nset
2925 :
2926 2092 : ldsa = SIZE(sphi_a, 1)
2927 2092 : ldsb = SIZE(sphi_b, 1)
2928 :
2929 2092 : NULLIFY (dblock)
2930 :
2931 2092 : ra = pbc(particle_set(iatom)%r(:), cell)
2932 8368 : rb(:) = ra(:) + rab(:)
2933 8368 : rac = ra - rc
2934 8368 : rbc = rb - rc
2935 8368 : dab = norm2(rab)
2936 :
2937 2092 : ic = cell_to_index(icell(1), icell(2), icell(3))
2938 :
2939 6276 : DO iset = 1, nseta
2940 :
2941 2092 : ncoa = npgfa(iset)*ncoset(la_max(iset))
2942 2092 : sgfa = first_sgfa(1, iset)
2943 :
2944 6276 : DO jset = 1, nsetb
2945 :
2946 2092 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
2947 :
2948 2092 : ncob = npgfb(jset)*ncoset(lb_max(jset))
2949 2092 : sgfb = first_sgfb(1, jset)
2950 :
2951 : CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
2952 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), 1, &
2953 2092 : rac, rbc, dipab)
2954 10460 : DO i_dir = 1, 3
2955 : CALL dbcsr_get_block_p(matrix=moments_rs_img(i_dir, ic)%matrix, &
2956 6276 : row=iatom, col=jatom, BLOCK=dblock, found=found)
2957 6276 : CPASSERT(found)
2958 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
2959 : 1.0_dp, dipab(1, 1, i_dir), ldwork, &
2960 6276 : sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
2961 :
2962 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
2963 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
2964 14644 : work(1, 1), ldwork, 1.0_dp, dblock(1, 1), SIZE(dblock, 1))
2965 : END DO
2966 : END DO
2967 : END DO
2968 : END ASSOCIATE
2969 : END DO
2970 8 : CALL neighbor_list_iterator_release(nl_iterator)
2971 8 : CALL kpoint_release(kpoints_all)
2972 8 : DEALLOCATE (dipab, work, basis_set_list)
2973 8 : CALL timestop(handle)
2974 :
2975 24 : END SUBROUTINE build_local_moment_matrix_rs_img
2976 :
2977 : ! **************************************************************************************************
2978 : !> \brief Calculates the dipole moments and berry curvature for periodic systems for kpoints
2979 : !> \param qs_env ...
2980 : !> \param xkp list of kpoints
2981 : !> \param dipole ...
2982 : !> \param rcc coordinates about which to calculate the dipole
2983 : !> \param berry_c berry curvature calculated using Ω^γ_n = Σ_m 2*Im[d^α_nm (d^β_mn)*]
2984 : !> \param do_parallel option to distribute the result in dipole across
2985 : !> different MPI ranks
2986 : !> \author Shridhar Shanbhag
2987 : ! **************************************************************************************************
2988 8 : SUBROUTINE qs_moment_kpoints_deep(qs_env, xkp, dipole, rcc, berry_c, do_parallel)
2989 : TYPE(qs_environment_type), POINTER :: qs_env
2990 : LOGICAL, OPTIONAL :: do_parallel
2991 : LOGICAL :: my_do_parallel, calc_bc
2992 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_moment_kpoints_deep'
2993 : COMPLEX(KIND=dp) :: phase, tmp_max
2994 8 : COMPLEX(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: C_k, H_k, S_k, D_k, CDC, C_dH_C, &
2995 8 : C_dS_C, dH_dk_i, dS_dk_i
2996 8 : COMPLEX(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dip
2997 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
2998 : ALLOCATABLE :: dipole
2999 : INTEGER :: handle, i_dir, ikp, nkp, &
3000 : n_img_scf, n_img_all, nao, &
3001 : num_pe, num_copy, mepos, n, m, mu, &
3002 : ispin, nspin
3003 : INTEGER, DIMENSION(3) :: periodic
3004 8 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_all
3005 8 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_all
3006 : REAL(KIND=dp), DIMENSION(3), OPTIONAL :: rcc
3007 : REAL(KIND=dp), DIMENSION(3) :: my_rcc
3008 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
3009 8 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: eigenvals
3010 8 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: bc, xkp
3011 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
3012 : ALLOCATABLE, OPTIONAL :: berry_c
3013 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
3014 8 : ALLOCATABLE :: D_rs, H_rs, S_rs
3015 : TYPE(cell_type), POINTER :: cell
3016 8 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_rs_img, matrix_ks_kp, &
3017 8 : matrix_s_kp
3018 : TYPE(dft_control_type), POINTER :: dft_control
3019 : TYPE(kpoint_type), POINTER :: kpoints_all, kpoints_scf
3020 8 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3021 : TYPE(mp_para_env_type), POINTER :: para_env
3022 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3023 8 : POINTER :: sab_all
3024 :
3025 8 : CALL timeset(routineN, handle)
3026 8 : calc_bc = PRESENT(berry_c)
3027 8 : my_do_parallel = .FALSE.
3028 8 : my_rcc = 0.0_dp
3029 8 : IF (PRESENT(do_parallel)) my_do_parallel = do_parallel
3030 8 : IF (PRESENT(rcc)) my_rcc = rcc
3031 :
3032 : CALL get_qs_env(qs_env, &
3033 : matrix_ks_kp=matrix_ks_kp, &
3034 : matrix_s_kp=matrix_s_kp, &
3035 : sab_all=sab_all, &
3036 : cell=cell, &
3037 : kpoints=kpoints_scf, &
3038 : para_env=para_env, &
3039 : dft_control=dft_control, &
3040 8 : mos=mos)
3041 :
3042 8 : CALL get_mo_set(mo_set=mos(1), nao=nao)
3043 8 : CALL get_cell(cell=cell, h=hmat, periodic=periodic)
3044 8 : nspin = SIZE(matrix_ks_kp, 1)
3045 8 : nkp = SIZE(xkp, 2)
3046 :
3047 : ! create kpoint environment kpoints_all which contains all neighbor cells R
3048 : ! without considering any lattice symmetry
3049 8 : NULLIFY (kpoints_all)
3050 8 : CALL kpoint_create(kpoints_all)
3051 8 : CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, n_img_scf)
3052 :
3053 : CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index_all, &
3054 8 : index_to_cell=index_to_cell_all)
3055 8 : n_img_all = SIZE(index_to_cell_all, 2)
3056 :
3057 8 : NULLIFY (moments_rs_img)
3058 8 : CALL dbcsr_allocate_matrix_set(moments_rs_img, 3, n_img_all)
3059 : ! D_μ,ν = <φ_μ|r|φ_ν>
3060 8 : CALL build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc=my_rcc)
3061 :
3062 80 : ALLOCATE (S_rs(1, nao, nao, n_img_all), H_rs(nspin, nao, nao, n_img_all), source=0.0_dp)
3063 40 : ALLOCATE (D_rs(3, nao, nao, n_img_all), source=0.0_dp)
3064 :
3065 : ! Convert real-space dbcsr matrices into arrays
3066 8 : CALL replicate_rs_matrices(matrix_s_kp, kpoints_scf, S_rs, cell_to_index_all)
3067 8 : CALL replicate_rs_matrices(matrix_ks_kp, kpoints_scf, H_rs, cell_to_index_all)
3068 8 : CALL replicate_rs_matrices(moments_rs_img, kpoints_all, D_rs, cell_to_index_all)
3069 :
3070 8 : mepos = 0
3071 8 : num_pe = 1
3072 8 : num_copy = nkp
3073 8 : IF (my_do_parallel) THEN
3074 8 : mepos = para_env%mepos
3075 8 : num_pe = para_env%num_pe
3076 8 : num_copy = CEILING(REAL(nkp)/num_pe)
3077 : END IF
3078 :
3079 56 : ALLOCATE (dipole(nspin, num_copy, 3, nao, nao), source=z_zero)
3080 32 : IF (calc_bc) ALLOCATE (berry_c(nspin, num_copy, 3, nao), source=0.0_dp)
3081 :
3082 : !$OMP PARALLEL DEFAULT(NONE) PRIVATE(ikp, S_k, H_k, eigenvals, C_k, ispin, n, m, &
3083 : !$OMP i_dir, dS_dk_i, dH_dk_i, D_k, dip, bc, C_dS_C, C_dH_C, CDC, tmp_max, phase) &
3084 : !$OMP SHARED(num_pe, mepos, dipole, berry_c, nao, nspin, periodic, &
3085 8 : !$OMP nkp, xkp, S_rs, H_rs, D_rs, index_to_cell_all, hmat, calc_bc)
3086 : ALLOCATE (dS_dk_i(nao, nao), C_dS_C(nao, nao), dH_dk_i(nao, nao), C_dH_C(nao, nao), source=z_zero)
3087 : ALLOCATE (CDC(nao, nao), dip(3, nao, nao), S_k(nao, nao), H_k(nao, nao), source=z_zero)
3088 : ALLOCATE (C_k(nao, nao), D_k(nao, nao), source=z_zero)
3089 : ALLOCATE (eigenvals(nao), source=0.0_dp)
3090 : IF (calc_bc) ALLOCATE (bc(3, nao), source=0.0_dp)
3091 : !$OMP DO COLLAPSE(2)
3092 : DO ispin = 1, nspin
3093 : DO ikp = 1, nkp
3094 : IF (MOD(ikp - 1, num_pe) /= mepos) CYCLE
3095 :
3096 : ! S^R -> S(k), H^R -> H(k)
3097 : S_k = 0
3098 : H_k = 0
3099 : CALL rs_to_kp(S_rs(1, :, :, :), S_k, index_to_cell_all, xkp(:, ikp))
3100 : CALL rs_to_kp(H_rs(ispin, :, :, :), H_k, index_to_cell_all, xkp(:, ikp))
3101 :
3102 : ! Diagonalize H(k)C(k) = S(k)C(k)ε(k)
3103 : CALL geeig_right(H_k, S_k, eigenvals, C_k)
3104 :
3105 : ! To have a smooth complex phase of C(k) as function of k, for every n, we force
3106 : ! the largest C_μ,n(k) to be real.
3107 : ! This is important to have a continuous dipole moment d_nm(k) as a function of k
3108 : DO n = 1, nao
3109 : tmp_max = C_k(1, n)
3110 : DO mu = 1, nao
3111 : IF (ABS(C_k(mu, n)) < ABS(tmp_max)) CYCLE
3112 : tmp_max = C_k(mu, n)
3113 : END DO
3114 : phase = tmp_max/ABS(tmp_max)
3115 : C_k(:, n) = C_k(:, n)/phase
3116 : END DO
3117 :
3118 : DO i_dir = 1, 3 ! d^x, d^y, d^z
3119 :
3120 : IF (periodic(i_dir) == 0) CYCLE
3121 : ! ∇ S(k) = Σ_R iR S^R e^(ikR), ∇ H(k) = Σ_R iR H^R e^(ikR)
3122 : CALL rs_to_kp(S_rs(1, :, :, :), dS_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
3123 : CALL rs_to_kp(H_rs(ispin, :, :, :), dH_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
3124 :
3125 : ! Σ_R D^R e^(ikR) = D(k), D_μ,ν = <φ_μ|r|φ_ν>
3126 : CALL rs_to_kp(D_rs(i_dir, :, :, :), D_k(:, :), index_to_cell_all, xkp(:, ikp))
3127 :
3128 : ! Basis transform to Kohn-Sham basis: (C^H) ∇ S C, (C^H) ∇ H C, (C^H) D C
3129 : CALL gemm_square(C_k, 'C', dS_dk_i, 'N', C_k, 'N', C_dS_C)
3130 : CALL gemm_square(C_k, 'C', dH_dk_i, 'N', C_k, 'N', C_dH_C)
3131 : CALL gemm_square(C_k, 'C', D_k, 'N', C_k, 'N', CDC)
3132 :
3133 : ! Compute the dipole
3134 : ! d_nm (k) = - i/(ε(n)-ε(m)) [ (C^H)(dH(k)/dk)C ]_nm
3135 : ! + i ε(n)/(ε(n)-ε(m)) [ (C^H)(dS(k)/dk)C ]_nm + [ (C^H)D(k)C ]_nm
3136 : DO n = 1, nao
3137 : DO m = 1, nao
3138 : IF (n == m) CYCLE ! diagonal elements would need to be computed from
3139 : ! a numerical k-derivative which is not implemented
3140 : dip(i_dir, n, m) = -gaussi*C_dH_C(n, m)/(eigenvals(n) - eigenvals(m)) &
3141 : + gaussi*eigenvals(n)*C_dS_C(n, m)/(eigenvals(n) - eigenvals(m)) &
3142 : + CDC(n, m)
3143 : END DO
3144 : END DO
3145 : END DO
3146 : ! Compute the Berry curvature from the dipoles
3147 : ! Ω^γ_n = Σ_m 2*Im[d^α_nm d^β_mn], where, α, β, γ belong to {x, y, z}
3148 : IF (calc_bc) THEN
3149 : bc = 0.0_dp
3150 : DO i_dir = 1, 3
3151 : DO n = 1, nao
3152 : DO m = 1, nao
3153 : IF (n == m) CYCLE
3154 : bc(i_dir, n) = bc(i_dir, n) &
3155 : + 2*AIMAG(dip(1 + MOD(i_dir, 3), n, m)*dip(1 + MOD(i_dir + 1, 3), m, n))
3156 : END DO
3157 : END DO
3158 : END DO
3159 : END IF
3160 : ! Store the dipoles and berry curvature for each MPI rank
3161 : dipole(ispin, CEILING(REAL(ikp)/num_pe), :, :, :) = dip(:, :, :)
3162 : IF (calc_bc) berry_c(ispin, CEILING(REAL(ikp)/num_pe), :, :) = bc(:, :)
3163 : END DO
3164 : END DO
3165 : !$OMP END DO
3166 : DEALLOCATE (dS_dk_i, C_dS_C, dH_dk_i, C_dH_C, CDC, dip, S_k, H_k, C_k, D_k, eigenvals)
3167 : IF (calc_bc) DEALLOCATE (bc)
3168 : !$OMP END PARALLEL
3169 8 : DEALLOCATE (S_rs, H_rs, D_rs)
3170 8 : CALL dbcsr_deallocate_matrix_set(moments_rs_img)
3171 8 : CALL kpoint_release(kpoints_all)
3172 8 : CALL timestop(handle)
3173 32 : END SUBROUTINE qs_moment_kpoints_deep
3174 :
3175 : ! **************************************************************************************************
3176 : !> \brief Calculates interband k-point dipoles in the existing SCF MO basis.
3177 : !> \param qs_env ...
3178 : !> \param dipole ...
3179 : !> \param rcc retained for interface compatibility; interband dipoles are origin independent
3180 : !> \param nmo_spin_out number of SCF MOs available for each spin
3181 : ! **************************************************************************************************
3182 6 : SUBROUTINE qs_moment_kpoints_scf_mos(qs_env, dipole, rcc, nmo_spin_out)
3183 : TYPE(qs_environment_type), POINTER :: qs_env
3184 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
3185 : ALLOCATABLE :: dipole
3186 : REAL(KIND=dp), DIMENSION(3), OPTIONAL :: rcc
3187 : INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT), &
3188 : OPTIONAL :: nmo_spin_out
3189 :
3190 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_moment_kpoints_scf_mos'
3191 :
3192 : INTEGER :: handle, i_dir, ikp, ikp_local, ispin, &
3193 : m, n, nao, nkp, nmo, nspin
3194 6 : INTEGER, DIMENSION(:), ALLOCATABLE :: nmo_spin
3195 : INTEGER, DIMENSION(2) :: kp_range
3196 6 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3197 : LOGICAL :: my_kpgrp
3198 : REAL(KIND=dp), PARAMETER :: eps_degenerate = 1.0E-10_dp
3199 : REAL(KIND=dp) :: cimag, creal, energy_diff
3200 6 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: eigenvalues_kp
3201 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvals
3202 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_all
3203 : TYPE(cp_fm_struct_type), POINTER :: moment_struct
3204 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
3205 : TYPE(cp_fm_type) :: fm_dummy, fm_tmp, mo_coeff_im_global, &
3206 : mo_coeff_re_global, moment_im, &
3207 : moment_re
3208 : TYPE(cp_fm_type), POINTER :: mo_coeff_im, mo_coeff_re
3209 6 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: overlap_deriv
3210 : TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix
3211 : TYPE(dft_control_type), POINTER :: dft_control
3212 6 : TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env
3213 : TYPE(kpoint_env_type), POINTER :: kp
3214 : TYPE(kpoint_type), POINTER :: kpoints_scf
3215 6 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp
3216 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_kp
3217 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3218 6 : POINTER :: sab_kp, sab_orb
3219 : TYPE(qs_ks_env_type), POINTER :: ks_env
3220 :
3221 6 : CALL timeset(routineN, handle)
3222 :
3223 6 : NULLIFY (blacs_env_all, cell_to_index, cmatrix, dft_control, eigenvals, fm_struct, kp, &
3224 6 : kp_env, kpoints_scf, ks_env, mo_coeff_im, mo_coeff_re, moment_struct, mos_kp, &
3225 6 : overlap_deriv, para_env, para_env_kp, rmatrix, sab_kp, sab_orb)
3226 : IF (PRESENT(rcc)) THEN
3227 : MARK_USED(rcc)
3228 : END IF
3229 :
3230 : CALL get_qs_env(qs_env, dft_control=dft_control, kpoints=kpoints_scf, ks_env=ks_env, &
3231 6 : para_env=para_env, sab_orb=sab_orb)
3232 6 : CPASSERT(ASSOCIATED(dft_control))
3233 6 : CPASSERT(ASSOCIATED(kpoints_scf))
3234 6 : CPASSERT(ASSOCIATED(ks_env))
3235 6 : CPASSERT(ASSOCIATED(para_env))
3236 6 : CPASSERT(ASSOCIATED(sab_orb))
3237 :
3238 : CALL get_kpoint_info(kpoints_scf, nkp=nkp, kp_range=kp_range, kp_env=kp_env, &
3239 : para_env_kp=para_env_kp, blacs_env_all=blacs_env_all, &
3240 6 : cell_to_index=cell_to_index, sab_nl=sab_kp)
3241 6 : IF (kp_range(2) >= kp_range(1)) THEN
3242 6 : CPASSERT(ASSOCIATED(kp_env))
3243 : END IF
3244 6 : CPASSERT(ASSOCIATED(para_env_kp))
3245 6 : CPASSERT(ASSOCIATED(blacs_env_all))
3246 6 : CPASSERT(ASSOCIATED(cell_to_index))
3247 6 : CPASSERT(ASSOCIATED(sab_kp))
3248 :
3249 : CALL build_overlap_matrix(ks_env, matrixkp_s=overlap_deriv, nderivative=1, &
3250 : basis_type_a="ORB", basis_type_b="ORB", sab_nl=sab_orb, &
3251 6 : ext_kpoints=kpoints_scf)
3252 :
3253 6 : nspin = dft_control%nspins
3254 6 : CALL dbcsr_get_info(overlap_deriv(1, 1)%matrix, nfullrows_total=nao)
3255 18 : ALLOCATE (nmo_spin(nspin), source=0)
3256 6 : IF (kp_range(2) >= kp_range(1)) THEN
3257 6 : kp => kp_env(1)%kpoint_env
3258 6 : mos_kp => kp%mos
3259 6 : CPASSERT(ASSOCIATED(mos_kp))
3260 12 : DO ispin = 1, nspin
3261 12 : CALL get_mo_set(mos_kp(1, ispin), nmo=nmo_spin(ispin))
3262 : END DO
3263 : END IF
3264 6 : CALL para_env%max(nmo_spin)
3265 48 : ALLOCATE (dipole(nspin, nkp, 3, MAXVAL(nmo_spin), MAXVAL(nmo_spin)), source=z_zero)
3266 6 : IF (PRESENT(nmo_spin_out)) THEN
3267 8 : ALLOCATE (nmo_spin_out(nspin))
3268 8 : nmo_spin_out(:) = nmo_spin(:)
3269 : END IF
3270 :
3271 6 : ALLOCATE (rmatrix, cmatrix)
3272 : CALL dbcsr_create(rmatrix, template=overlap_deriv(1, 1)%matrix, &
3273 6 : matrix_type=dbcsr_type_antisymmetric)
3274 : CALL dbcsr_create(cmatrix, template=overlap_deriv(1, 1)%matrix, &
3275 6 : matrix_type=dbcsr_type_symmetric)
3276 6 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_kp)
3277 6 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_kp)
3278 :
3279 258 : DO ikp = 1, nkp
3280 252 : my_kpgrp = (ikp >= kp_range(1) .AND. ikp <= kp_range(2))
3281 : IF (my_kpgrp) THEN
3282 234 : ikp_local = ikp - kp_range(1) + 1
3283 234 : kp => kp_env(ikp_local)%kpoint_env
3284 234 : mos_kp => kp%mos
3285 : ELSE
3286 252 : NULLIFY (kp, mos_kp)
3287 : END IF
3288 510 : DO ispin = 1, nspin
3289 252 : nmo = nmo_spin(ispin)
3290 756 : ALLOCATE (eigenvalues_kp(nmo), source=0.0_dp)
3291 :
3292 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nmo, &
3293 252 : para_env=para_env, context=blacs_env_all)
3294 252 : CALL cp_fm_create(mo_coeff_re_global, fm_struct)
3295 252 : CALL cp_fm_create(mo_coeff_im_global, fm_struct)
3296 252 : CALL cp_fm_create(fm_tmp, fm_struct)
3297 252 : CALL cp_fm_struct_release(fm_struct)
3298 : CALL cp_fm_struct_create(moment_struct, nrow_global=nmo, ncol_global=nmo, &
3299 252 : para_env=para_env, context=blacs_env_all)
3300 252 : CALL cp_fm_create(moment_re, moment_struct)
3301 252 : CALL cp_fm_create(moment_im, moment_struct)
3302 252 : CALL cp_fm_struct_release(moment_struct)
3303 :
3304 252 : IF (my_kpgrp) THEN
3305 234 : CALL get_mo_set(mos_kp(1, ispin), eigenvalues=eigenvals, mo_coeff=mo_coeff_re)
3306 234 : CALL get_mo_set(mos_kp(2, ispin), mo_coeff=mo_coeff_im)
3307 234 : CPASSERT(ASSOCIATED(eigenvals))
3308 234 : CPASSERT(ASSOCIATED(mo_coeff_re))
3309 234 : CPASSERT(ASSOCIATED(mo_coeff_im))
3310 772 : IF (para_env_kp%is_source()) eigenvalues_kp(1:nmo) = eigenvals(1:nmo)
3311 234 : CALL cp_fm_copy_general(mo_coeff_re, mo_coeff_re_global, para_env)
3312 234 : CALL cp_fm_copy_general(mo_coeff_im, mo_coeff_im_global, para_env)
3313 : ELSE
3314 18 : CALL cp_fm_copy_general(fm_dummy, mo_coeff_re_global, para_env)
3315 18 : CALL cp_fm_copy_general(fm_dummy, mo_coeff_im_global, para_env)
3316 : END IF
3317 252 : CALL para_env%sum(eigenvalues_kp)
3318 :
3319 1008 : DO i_dir = 1, 3
3320 756 : CALL dbcsr_set(rmatrix, 0.0_dp)
3321 756 : CALL dbcsr_set(cmatrix, 0.0_dp)
3322 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=overlap_deriv, &
3323 : ispin=i_dir + 1, xkp=kpoints_scf%xkp(:, ikp), &
3324 756 : cell_to_index=cell_to_index, sab_nl=sab_kp)
3325 :
3326 : ! Project the complex AO derivative operator as C^H A C. The
3327 : ! off-diagonal length-gauge dipoles follow from the energy-gap relation.
3328 756 : CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_re_global, fm_tmp, nmo)
3329 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3330 756 : 1.0_dp, mo_coeff_re_global, fm_tmp, 0.0_dp, moment_re)
3331 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3332 756 : -1.0_dp, mo_coeff_im_global, fm_tmp, 0.0_dp, moment_im)
3333 :
3334 756 : CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_im_global, fm_tmp, nmo)
3335 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3336 756 : 1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
3337 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3338 756 : 1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)
3339 :
3340 756 : CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_re_global, fm_tmp, nmo)
3341 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3342 756 : 1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
3343 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3344 756 : 1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)
3345 :
3346 756 : CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_im_global, fm_tmp, nmo)
3347 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3348 756 : -1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_re)
3349 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
3350 756 : 1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_im)
3351 :
3352 4236 : DO n = 1, nmo
3353 18108 : DO m = 1, nmo
3354 14124 : IF (n == m) CYCLE
3355 10896 : energy_diff = eigenvalues_kp(m) - eigenvalues_kp(n)
3356 10896 : IF (ABS(energy_diff) <= eps_degenerate) CYCLE
3357 10872 : CALL cp_fm_get_element(moment_re, m, n, creal)
3358 10872 : CALL cp_fm_get_element(moment_im, m, n, cimag)
3359 14100 : IF (para_env%is_source()) then
3360 8688 : dipole(ispin, ikp, i_dir, n, m) = CMPLX(creal, cimag, KIND=dp)/energy_diff
3361 : end if
3362 : END DO
3363 : END DO
3364 : END DO
3365 252 : CALL cp_fm_release(mo_coeff_im_global)
3366 252 : CALL cp_fm_release(mo_coeff_re_global)
3367 252 : CALL cp_fm_release(moment_im)
3368 252 : CALL cp_fm_release(moment_re)
3369 252 : CALL cp_fm_release(fm_tmp)
3370 1008 : DEALLOCATE (eigenvalues_kp)
3371 : END DO
3372 : END DO
3373 :
3374 12 : DO ispin = 1, nspin
3375 264 : DO ikp = 1, nkp
3376 1014 : DO i_dir = 1, 3
3377 35712 : CALL para_env%sum(dipole(ispin, ikp, i_dir, :, :))
3378 : END DO
3379 : END DO
3380 : END DO
3381 :
3382 6 : CALL dbcsr_deallocate_matrix(cmatrix)
3383 6 : CALL dbcsr_deallocate_matrix(rmatrix)
3384 6 : CALL dbcsr_deallocate_matrix_set(overlap_deriv)
3385 6 : DEALLOCATE (nmo_spin)
3386 6 : CALL timestop(handle)
3387 :
3388 18 : END SUBROUTINE qs_moment_kpoints_scf_mos
3389 :
3390 : ! **************************************************************************************************
3391 : !> \brief Calculate and print dipole moment elements d_nm(k) for k-point calculations
3392 : !> \param qs_env ...
3393 : !> \param nmoments ...
3394 : !> \param reference ...
3395 : !> \param ref_point ...
3396 : !> \param max_nmo ...
3397 : !> \param unit_number ...
3398 : !> \author Shridhar Shanbhag
3399 : ! **************************************************************************************************
3400 10 : SUBROUTINE qs_moment_kpoints(qs_env, nmoments, reference, ref_point, max_nmo, unit_number)
3401 : TYPE(qs_environment_type), POINTER :: qs_env
3402 : INTEGER, INTENT(IN) :: nmoments, reference, max_nmo
3403 : REAL(dp), DIMENSION(:), INTENT(IN), POINTER :: ref_point
3404 : INTEGER, INTENT(IN) :: unit_number
3405 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_moment_kpoints'
3406 10 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
3407 10 : COMPLEX(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dipole_to_print
3408 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
3409 10 : ALLOCATABLE :: dipole
3410 : INTEGER :: handle, i_dir, ikp, nmo_dim, nkp, nao, &
3411 : num_pe, mepos, n, m, &
3412 : ispin, nspin, nmin, nmax, homo
3413 10 : INTEGER, DIMENSION(:), ALLOCATABLE :: nmo_spin_scf
3414 : LOGICAL :: explicit_kpnts, explicit_kpset, use_scf_mos
3415 : REAL(KIND=dp), DIMENSION(3) :: rcc
3416 10 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: xkp
3417 10 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: bc_to_print
3418 10 : REAL(KIND=dp), DIMENSION(:, :, :, :), ALLOCATABLE :: berry_c
3419 10 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3420 : TYPE(mp_para_env_type), POINTER :: para_env
3421 : TYPE(section_vals_type), POINTER :: kpnts, kpset
3422 : CHARACTER(LEN=default_string_length), &
3423 10 : DIMENSION(:), POINTER :: special_pnts
3424 :
3425 10 : CALL timeset(routineN, handle)
3426 :
3427 10 : IF (nmoments > 1) CPABORT("KPOINT quadrupole and higher moments not implemented.")
3428 10 : IF (max_nmo < 0) CPABORT("Negative maximum number of molecular orbitals max_nmo provided.")
3429 :
3430 : CALL get_qs_env(qs_env, &
3431 : para_env=para_env, &
3432 : matrix_ks_kp=matrix_ks_kp, &
3433 10 : mos=mos)
3434 :
3435 10 : CALL get_mo_set(mo_set=mos(1), nao=nao)
3436 10 : CALL get_xkp_for_dipole_calc(qs_env, xkp, special_pnts)
3437 10 : CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
3438 10 : nspin = SIZE(matrix_ks_kp, 1)
3439 10 : nkp = SIZE(xkp, 2)
3440 :
3441 10 : kpset => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINT_SET")
3442 10 : kpnts => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINTS")
3443 10 : CALL section_vals_get(kpset, explicit=explicit_kpset)
3444 10 : CALL section_vals_get(kpnts, explicit=explicit_kpnts)
3445 10 : use_scf_mos = .NOT. explicit_kpset .AND. .NOT. explicit_kpnts
3446 :
3447 10 : IF (unit_number > 0) WRITE (unit_number, FMT="(/,T2,A)") &
3448 5 : '!-----------------------------------------------------------------------------!'
3449 10 : IF (unit_number > 0) WRITE (unit_number, "(T22,A)") "Periodic Dipole Matrix Elements"
3450 :
3451 10 : IF (use_scf_mos) THEN
3452 4 : CALL qs_moment_kpoints_scf_mos(qs_env, dipole, rcc, nmo_spin_scf)
3453 4 : nmo_dim = SIZE(dipole, 4)
3454 24 : ALLOCATE (berry_c(nspin, nkp, 3, nmo_dim), source=0.0_dp)
3455 8 : DO ispin = 1, nspin
3456 256 : DO ikp = 1, nkp
3457 996 : DO i_dir = 1, 3
3458 4160 : DO n = 1, nmo_dim
3459 17736 : DO m = 1, nmo_dim
3460 13824 : IF (n == m) CYCLE
3461 : berry_c(ispin, ikp, i_dir, n) = berry_c(ispin, ikp, i_dir, n) &
3462 : + 2*AIMAG(dipole(ispin, ikp, 1 + MOD(i_dir, 3), n, m)* &
3463 16992 : dipole(ispin, ikp, 1 + MOD(i_dir + 1, 3), m, n))
3464 : END DO
3465 : END DO
3466 : END DO
3467 : END DO
3468 : END DO
3469 : ELSE
3470 : CALL qs_moment_kpoints_deep(qs_env, &
3471 : xkp, &
3472 : dipole, &
3473 : rcc, &
3474 : berry_c, &
3475 6 : do_parallel=.TRUE.)
3476 6 : nmo_dim = nao
3477 : END IF
3478 :
3479 40 : ALLOCATE (dipole_to_print(3, nmo_dim, nmo_dim), source=z_zero)
3480 30 : ALLOCATE (bc_to_print(3, nmo_dim), source=0.0_dp)
3481 :
3482 10 : mepos = para_env%mepos
3483 10 : num_pe = para_env%num_pe
3484 :
3485 264 : DO ikp = 1, nkp
3486 520 : DO ispin = 1, nspin
3487 256 : CALL get_mo_set(mo_set=mos(ispin), homo=homo)
3488 256 : nmin = max(1, homo - (max_nmo - 1)/2)
3489 256 : nmax = min(nao, homo + max_nmo/2)
3490 256 : IF (max_nmo == 0) THEN
3491 0 : nmin = 1
3492 0 : nmax = nao
3493 : END IF
3494 256 : IF (use_scf_mos) THEN
3495 248 : nmax = min(nmax, nmo_spin_scf(ispin))
3496 : END IF
3497 256 : dipole_to_print = 0.0_dp
3498 256 : bc_to_print = 0.0_dp
3499 256 : IF (use_scf_mos) THEN
3500 19736 : dipole_to_print(:, :, :) = dipole(ispin, ikp, :, :, :)
3501 4472 : bc_to_print(:, :) = berry_c(ispin, ikp, :, :)
3502 8 : ELSE IF (mod(ikp - 1, num_pe) == mepos) THEN
3503 87268 : dipole_to_print(:, :, :) = dipole(ispin, CEILING(REAL(ikp)/num_pe), :, :, :)
3504 900 : bc_to_print(:, :) = berry_c(ispin, CEILING(REAL(ikp)/num_pe), :, :)
3505 : END IF
3506 256 : IF (.NOT. use_scf_mos) THEN
3507 8 : CALL para_env%sum(dipole_to_print)
3508 8 : CALL para_env%sum(bc_to_print)
3509 : END IF
3510 766 : IF (unit_number > 0) THEN
3511 128 : IF (special_pnts(ikp) /= "") WRITE (unit_number, "(/,2X,A,A)") &
3512 3 : "Special point: ", ADJUSTL(TRIM(special_pnts(ikp)))
3513 : WRITE (unit_number, "(/,1X,A,I3,1X,3(A,1F12.6))") &
3514 128 : "Kpoint:", ikp, ", kx:", xkp(1, ikp), ", ky:", xkp(2, ikp), ", kz:", xkp(3, ikp)
3515 128 : IF (nspin > 1) WRITE (unit_number, "(/,2X,A,I2)") "Open Shell System. Spin:", ispin
3516 : WRITE (unit_number, "(2X,A)") " kp n m Re(dx_nm) Im(dx_nm) &
3517 128 : & Re(dy_nm) Im(dy_nm) Re(dz_nm) Im(dz_nm)"
3518 676 : DO n = nmin, nmax
3519 3116 : DO m = nmin, nmax
3520 2440 : IF (n == m) CYCLE
3521 2988 : WRITE (unit_number, "(2X,I4,2I4,6(G11.3))") ikp, n, m, dipole_to_print(1:3, n, m)
3522 : END DO
3523 : END DO
3524 128 : WRITE (unit_number, "(/,1X,A)") "Berry Curvature"
3525 128 : WRITE (unit_number, "(2X,A)") " kp n YZ ZX XY"
3526 676 : DO n = nmin, nmax
3527 : WRITE (unit_number, "(2X,2I5,3(1X,G11.3))") &
3528 676 : ikp, n, bc_to_print(1, n), bc_to_print(2, n), bc_to_print(3, n)
3529 : END DO
3530 : END IF
3531 : END DO
3532 : END DO
3533 10 : DEALLOCATE (dipole_to_print, bc_to_print, berry_c, dipole)
3534 10 : IF (ALLOCATED(nmo_spin_scf)) DEALLOCATE (nmo_spin_scf)
3535 10 : DEALLOCATE (special_pnts, xkp)
3536 :
3537 10 : CALL timestop(handle)
3538 :
3539 40 : END SUBROUTINE qs_moment_kpoints
3540 :
3541 : END MODULE qs_moments
|
