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