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