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