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