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