Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 : USE commutator_rpnl, ONLY: build_com_mom_nl
32 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_lu_decompose
33 : USE cp_cfm_types, ONLY: cp_cfm_create,&
34 : cp_cfm_get_info,&
35 : cp_cfm_release,&
36 : cp_cfm_type
37 : USE cp_control_types, ONLY: dft_control_type
38 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
39 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
40 : dbcsr_allocate_matrix_set,&
41 : dbcsr_deallocate_matrix_set
42 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
43 : cp_fm_struct_double,&
44 : cp_fm_struct_release,&
45 : cp_fm_struct_type
46 : USE cp_fm_types, ONLY: cp_fm_create,&
47 : cp_fm_get_info,&
48 : cp_fm_release,&
49 : cp_fm_set_all,&
50 : cp_fm_type
51 : USE cp_result_methods, ONLY: cp_results_erase,&
52 : put_results
53 : USE cp_result_types, ONLY: cp_result_type
54 : USE dbcsr_api, ONLY: &
55 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_type, dbcsr_dot, &
56 : dbcsr_get_block_p, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_set, dbcsr_trace, &
57 : dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
58 : USE distribution_1d_types, ONLY: distribution_1d_type
59 : USE kinds, ONLY: default_string_length,&
60 : dp
61 : USE kpoint_types, ONLY: get_kpoint_info,&
62 : kpoint_type
63 : USE mathconstants, ONLY: pi,&
64 : twopi
65 : USE message_passing, ONLY: mp_para_env_type
66 : USE moments_utils, ONLY: get_reference_point
67 : USE orbital_pointers, ONLY: current_maxl,&
68 : indco,&
69 : ncoset
70 : USE parallel_gemm_api, ONLY: parallel_gemm
71 : USE particle_methods, ONLY: get_particle_set
72 : USE particle_types, ONLY: particle_type
73 : USE physcon, ONLY: bohr,&
74 : debye
75 : USE qs_environment_types, ONLY: get_qs_env,&
76 : qs_environment_type
77 : USE qs_kind_types, ONLY: get_qs_kind,&
78 : get_qs_kind_set,&
79 : qs_kind_type
80 : USE qs_ks_types, ONLY: get_ks_env,&
81 : qs_ks_env_type
82 : USE qs_mo_types, ONLY: get_mo_set,&
83 : mo_set_type
84 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
85 : neighbor_list_iterate,&
86 : neighbor_list_iterator_create,&
87 : neighbor_list_iterator_p_type,&
88 : neighbor_list_iterator_release,&
89 : neighbor_list_set_p_type
90 : USE qs_operators_ao, ONLY: build_lin_mom_matrix
91 : USE qs_rho_types, ONLY: qs_rho_get,&
92 : qs_rho_type
93 : USE rt_propagation_types, ONLY: get_rtp,&
94 : rt_prop_type
95 : #include "./base/base_uses.f90"
96 :
97 : IMPLICIT NONE
98 :
99 : PRIVATE
100 :
101 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_moments'
102 :
103 : ! Public subroutines
104 : PUBLIC :: build_berry_moment_matrix, build_local_moment_matrix
105 : PUBLIC :: build_berry_kpoint_matrix, build_local_magmom_matrix
106 : PUBLIC :: qs_moment_berry_phase, qs_moment_locop
107 : PUBLIC :: dipole_deriv_ao
108 : PUBLIC :: build_local_moments_der_matrix
109 : PUBLIC :: build_dsdv_moments
110 : PUBLIC :: dipole_velocity_deriv
111 :
112 : CONTAINS
113 :
114 : ! *****************************************************************************
115 : !> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
116 : !> to the basis function on the right
117 : !> difdip(beta, alpha) = < mu | r_beta | ∂_alpha nu > * (mu - nu)
118 : !> \param qs_env ...
119 : !> \param difdip ...
120 : !> \param order The order of the derivative (1 for dipole moment)
121 : !> \param lambda The atom on which we take the derivative
122 : !> \param rc ...
123 : !> \author Edward Ditler
124 : ! **************************************************************************************************
125 6 : SUBROUTINE dipole_velocity_deriv(qs_env, difdip, order, lambda, rc)
126 : TYPE(qs_environment_type), POINTER :: qs_env
127 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: difdip
128 : INTEGER, INTENT(IN) :: order, lambda
129 : REAL(KIND=dp), DIMENSION(3) :: rc
130 :
131 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dipole_velocity_deriv'
132 :
133 : INTEGER :: handle, i, iatom, icol, idir, ikind, inode, irow, iset, j, jatom, jkind, jset, &
134 : last_jatom, lda, ldab, ldb, M_dim, maxsgf, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
135 : sgfb
136 : LOGICAL :: found
137 : REAL(dp) :: dab
138 : REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
139 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
140 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
141 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab, difmab2
142 6 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mint, mint2
143 : TYPE(cell_type), POINTER :: cell
144 6 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
145 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
146 : TYPE(neighbor_list_iterator_p_type), &
147 6 : DIMENSION(:), POINTER :: nl_iterator
148 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
149 6 : POINTER :: sab_all
150 6 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
151 6 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
152 : TYPE(qs_kind_type), POINTER :: qs_kind
153 :
154 6 : CALL timeset(routineN, handle)
155 :
156 6 : NULLIFY (cell, particle_set, qs_kind_set, sab_all)
157 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, &
158 6 : qs_kind_set=qs_kind_set, sab_all=sab_all)
159 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
160 6 : maxco=ldab, maxsgf=maxsgf)
161 :
162 6 : nkind = SIZE(qs_kind_set)
163 6 : natom = SIZE(particle_set)
164 :
165 6 : M_dim = ncoset(order) - 1
166 :
167 30 : ALLOCATE (basis_set_list(nkind))
168 :
169 30 : ALLOCATE (mab(ldab, ldab, M_dim))
170 36 : ALLOCATE (difmab2(ldab, ldab, M_dim, 3))
171 24 : ALLOCATE (work(ldab, maxsgf))
172 78 : ALLOCATE (mint(3, 3))
173 78 : ALLOCATE (mint2(3, 3))
174 :
175 4920 : mab(1:ldab, 1:ldab, 1:M_dim) = 0.0_dp
176 14766 : difmab2(1:ldab, 1:ldab, 1:M_dim, 1:3) = 0.0_dp
177 414 : work(1:ldab, 1:maxsgf) = 0.0_dp
178 :
179 24 : DO i = 1, 3
180 78 : DO j = 1, 3
181 54 : NULLIFY (mint(i, j)%block)
182 72 : NULLIFY (mint2(i, j)%block)
183 : END DO
184 : END DO
185 :
186 : ! Set the basis_set_list(nkind) to point to the corresponding basis sets
187 18 : DO ikind = 1, nkind
188 12 : qs_kind => qs_kind_set(ikind)
189 12 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
190 18 : IF (ASSOCIATED(basis_set_a)) THEN
191 12 : basis_set_list(ikind)%gto_basis_set => basis_set_a
192 : ELSE
193 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
194 : END IF
195 : END DO
196 :
197 6 : CALL neighbor_list_iterator_create(nl_iterator, sab_all)
198 33 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
199 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
200 27 : iatom=iatom, jatom=jatom, r=rab)
201 :
202 27 : basis_set_a => basis_set_list(ikind)%gto_basis_set
203 27 : basis_set_b => basis_set_list(jkind)%gto_basis_set
204 27 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
205 27 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
206 :
207 : ASSOCIATE ( &
208 : ! basis ikind
209 : first_sgfa => basis_set_a%first_sgf, &
210 : la_max => basis_set_a%lmax, &
211 : la_min => basis_set_a%lmin, &
212 : npgfa => basis_set_a%npgf, &
213 : nsgfa => basis_set_a%nsgf_set, &
214 : rpgfa => basis_set_a%pgf_radius, &
215 : set_radius_a => basis_set_a%set_radius, &
216 : sphi_a => basis_set_a%sphi, &
217 : zeta => basis_set_a%zet, &
218 : ! basis jkind, &
219 : first_sgfb => basis_set_b%first_sgf, &
220 : lb_max => basis_set_b%lmax, &
221 : lb_min => basis_set_b%lmin, &
222 : npgfb => basis_set_b%npgf, &
223 : nsgfb => basis_set_b%nsgf_set, &
224 : rpgfb => basis_set_b%pgf_radius, &
225 : set_radius_b => basis_set_b%set_radius, &
226 : sphi_b => basis_set_b%sphi, &
227 : zetb => basis_set_b%zet)
228 :
229 27 : nseta = basis_set_a%nset
230 27 : nsetb = basis_set_b%nset
231 :
232 18 : IF (inode == 1) last_jatom = 0
233 :
234 : ! this guarantees minimum image convention
235 : ! anything else would not make sense
236 27 : IF (jatom == last_jatom) THEN
237 : CYCLE
238 : END IF
239 :
240 27 : last_jatom = jatom
241 :
242 27 : irow = iatom
243 27 : icol = jatom
244 :
245 108 : DO i = 1, 3
246 351 : DO j = 1, 3
247 243 : NULLIFY (mint(i, j)%block)
248 : CALL dbcsr_get_block_p(matrix=difdip(i, j)%matrix, &
249 : row=irow, col=icol, BLOCK=mint(i, j)%block, &
250 243 : found=found)
251 243 : CPASSERT(found)
252 2025 : mint(i, j)%block = 0._dp
253 : END DO
254 : END DO
255 :
256 : ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
257 27 : ra = pbc(particle_set(iatom)%r(:), cell)
258 108 : rb(:) = ra(:) + rab(:)
259 27 : rac = pbc(rc, ra, cell)
260 108 : rbc = rac + rab
261 108 : dab = norm2(rab)
262 :
263 81 : DO iset = 1, nseta
264 27 : ncoa = npgfa(iset)*ncoset(la_max(iset))
265 27 : sgfa = first_sgfa(1, iset)
266 81 : DO jset = 1, nsetb
267 27 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
268 27 : ncob = npgfb(jset)*ncoset(lb_max(jset))
269 27 : sgfb = first_sgfb(1, jset)
270 27 : ldab = MAX(ncoa, ncob)
271 27 : lda = ncoset(la_max(iset))*npgfa(iset)
272 27 : ldb = ncoset(lb_max(jset))*npgfb(jset)
273 162 : ALLOCATE (difmab(lda, ldb, M_dim, 3))
274 :
275 : ! Calculate integral difmab(beta, alpha) = (a|r_beta|db_alpha)
276 : ! difmab(beta, alpha) = < a | r_beta | ∂_alpha b >
277 : ! difmab(j, idir) = < a | r_j | ∂_idir b >
278 : CALL diff_momop_velocity(la_max(iset), npgfa(iset), zeta(:, iset), &
279 : rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
280 : zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
281 27 : difmab, lambda=lambda, iatom=iatom, jatom=jatom)
282 :
283 : ! *** Contraction step ***
284 :
285 108 : DO idir = 1, 3 ! derivative of AO function
286 351 : DO j = 1, 3 ! position operator r_j
287 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
288 : 1.0_dp, difmab(1, 1, j, idir), SIZE(difmab, 1), &
289 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
290 243 : 0.0_dp, work(1, 1), SIZE(work, 1))
291 :
292 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
293 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
294 : work(1, 1), SIZE(work, 1), &
295 : 1.0_dp, mint(j, idir)%block(sgfa, sgfb), &
296 324 : SIZE(mint(j, idir)%block, 1))
297 : END DO !j
298 : END DO !idir
299 54 : DEALLOCATE (difmab)
300 : END DO !jset
301 : END DO !iset
302 : END ASSOCIATE
303 : END DO!iterator
304 :
305 6 : CALL neighbor_list_iterator_release(nl_iterator)
306 :
307 24 : DO i = 1, 3
308 78 : DO j = 1, 3
309 72 : NULLIFY (mint(i, j)%block)
310 : END DO
311 : END DO
312 :
313 6 : DEALLOCATE (mab, difmab2, basis_set_list, work, mint, mint2)
314 :
315 6 : CALL timestop(handle)
316 18 : END SUBROUTINE dipole_velocity_deriv
317 :
318 : ! **************************************************************************************************
319 : !> \brief Builds the moments for the derivative of the overlap with respect to nuclear velocities
320 : !> \param qs_env ...
321 : !> \param moments ...
322 : !> \param nmoments ...
323 : !> \param ref_point ...
324 : !> \param ref_points ...
325 : !> \param basis_type ...
326 : !> \author Edward Ditler
327 : ! **************************************************************************************************
328 6 : SUBROUTINE build_dsdv_moments(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
329 :
330 : TYPE(qs_environment_type), POINTER :: qs_env
331 : TYPE(dbcsr_p_type), DIMENSION(:) :: moments
332 : INTEGER, INTENT(IN) :: nmoments
333 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
334 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
335 : OPTIONAL :: ref_points
336 : CHARACTER(len=*), OPTIONAL :: basis_type
337 :
338 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_dsdv_moments'
339 :
340 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
341 : maxco, maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
342 : INTEGER, DIMENSION(3) :: image_cell
343 : LOGICAL :: found
344 : REAL(KIND=dp) :: dab
345 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
346 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
347 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
348 6 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
349 : TYPE(cell_type), POINTER :: cell
350 6 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
351 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
352 : TYPE(neighbor_list_iterator_p_type), &
353 6 : DIMENSION(:), POINTER :: nl_iterator
354 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
355 6 : POINTER :: sab_orb
356 6 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
357 6 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
358 : TYPE(qs_kind_type), POINTER :: qs_kind
359 :
360 6 : IF (nmoments < 1) RETURN
361 :
362 6 : CALL timeset(routineN, handle)
363 :
364 6 : NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
365 :
366 6 : nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
367 6 : CPASSERT(SIZE(moments) >= nm)
368 :
369 6 : NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
370 : CALL get_qs_env(qs_env=qs_env, &
371 : qs_kind_set=qs_kind_set, &
372 : particle_set=particle_set, cell=cell, &
373 6 : sab_orb=sab_orb)
374 :
375 6 : nkind = SIZE(qs_kind_set)
376 6 : natom = SIZE(particle_set)
377 :
378 : ! Allocate work storage
379 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
380 : maxco=maxco, maxsgf=maxsgf, &
381 18 : basis_type=basis_type)
382 :
383 30 : ALLOCATE (mab(maxco, maxco, nm))
384 4920 : mab(:, :, :) = 0.0_dp
385 :
386 24 : ALLOCATE (work(maxco, maxsgf))
387 414 : work(:, :) = 0.0_dp
388 :
389 36 : ALLOCATE (mint(nm))
390 24 : DO i = 1, nm
391 24 : NULLIFY (mint(i)%block)
392 : END DO
393 :
394 30 : ALLOCATE (basis_set_list(nkind))
395 18 : DO ikind = 1, nkind
396 12 : qs_kind => qs_kind_set(ikind)
397 12 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
398 18 : IF (ASSOCIATED(basis_set_a)) THEN
399 12 : basis_set_list(ikind)%gto_basis_set => basis_set_a
400 : ELSE
401 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
402 : END IF
403 : END DO
404 6 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
405 24 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
406 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
407 18 : iatom=iatom, jatom=jatom, r=rab, cell=image_cell)
408 18 : basis_set_a => basis_set_list(ikind)%gto_basis_set
409 18 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
410 18 : basis_set_b => basis_set_list(jkind)%gto_basis_set
411 18 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
412 : ! basis ikind
413 : ASSOCIATE ( &
414 : first_sgfa => basis_set_a%first_sgf, &
415 : la_max => basis_set_a%lmax, &
416 : la_min => basis_set_a%lmin, &
417 : npgfa => basis_set_a%npgf, &
418 : nsgfa => basis_set_a%nsgf_set, &
419 : rpgfa => basis_set_a%pgf_radius, &
420 : set_radius_a => basis_set_a%set_radius, &
421 : sphi_a => basis_set_a%sphi, &
422 : zeta => basis_set_a%zet, &
423 : ! basis jkind, &
424 : first_sgfb => basis_set_b%first_sgf, &
425 : lb_max => basis_set_b%lmax, &
426 : lb_min => basis_set_b%lmin, &
427 : npgfb => basis_set_b%npgf, &
428 : nsgfb => basis_set_b%nsgf_set, &
429 : rpgfb => basis_set_b%pgf_radius, &
430 : set_radius_b => basis_set_b%set_radius, &
431 : sphi_b => basis_set_b%sphi, &
432 : zetb => basis_set_b%zet)
433 :
434 18 : nseta = basis_set_a%nset
435 18 : nsetb = basis_set_b%nset
436 :
437 15 : IF (inode == 1) last_jatom = 0
438 :
439 : ! this guarantees minimum image convention
440 : ! anything else would not make sense
441 18 : IF (jatom == last_jatom) THEN
442 : CYCLE
443 : END IF
444 :
445 18 : last_jatom = jatom
446 :
447 18 : IF (iatom <= jatom) THEN
448 12 : irow = iatom
449 12 : icol = jatom
450 : ELSE
451 6 : irow = jatom
452 6 : icol = iatom
453 : END IF
454 :
455 72 : DO i = 1, nm
456 54 : NULLIFY (mint(i)%block)
457 : CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
458 54 : row=irow, col=icol, BLOCK=mint(i)%block, found=found)
459 396 : mint(i)%block = 0._dp
460 : END DO
461 :
462 : ! fold atomic position back into unit cell
463 18 : IF (PRESENT(ref_points)) THEN
464 0 : rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
465 18 : ELSE IF (PRESENT(ref_point)) THEN
466 72 : rc(:) = ref_point(:)
467 : ELSE
468 0 : rc(:) = 0._dp
469 : END IF
470 : ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
471 : ! by folding around the center, such screwing can be avoided for a proper choice of center.
472 : ! we dont use PBC at this point
473 :
474 : ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
475 72 : ra(:) = particle_set(iatom)%r(:)
476 72 : rb(:) = ra(:) + rab(:)
477 18 : rac = pbc(rc, ra, cell)
478 72 : rbc = rac + rab
479 :
480 72 : dab = NORM2(rab)
481 :
482 54 : DO iset = 1, nseta
483 :
484 18 : ncoa = npgfa(iset)*ncoset(la_max(iset))
485 18 : sgfa = first_sgfa(1, iset)
486 :
487 54 : DO jset = 1, nsetb
488 :
489 18 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
490 :
491 18 : ncob = npgfb(jset)*ncoset(lb_max(jset))
492 18 : sgfb = first_sgfb(1, jset)
493 :
494 : ! Calculate the primitive integrals
495 : CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
496 : rpgfa(:, iset), la_min(iset), &
497 : lb_max(jset), npgfb(jset), zetb(:, jset), &
498 18 : rpgfb(:, jset), nmoments, rac, rbc, mab)
499 :
500 : ! Contraction step
501 90 : DO i = 1, nm
502 :
503 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
504 : 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
505 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
506 54 : 0.0_dp, work(1, 1), SIZE(work, 1))
507 :
508 72 : IF (iatom <= jatom) THEN
509 :
510 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
511 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
512 : work(1, 1), SIZE(work, 1), &
513 : 1.0_dp, mint(i)%block(sgfa, sgfb), &
514 36 : SIZE(mint(i)%block, 1))
515 :
516 : ELSE
517 :
518 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
519 : 1.0_dp, work(1, 1), SIZE(work, 1), &
520 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
521 : 1.0_dp, mint(i)%block(sgfb, sgfa), &
522 18 : SIZE(mint(i)%block, 1))
523 :
524 : END IF
525 :
526 : END DO
527 :
528 : END DO
529 : END DO
530 : END ASSOCIATE
531 :
532 : END DO ! iterator
533 :
534 6 : CALL neighbor_list_iterator_release(nl_iterator)
535 :
536 : ! Release work storage
537 6 : DEALLOCATE (mab, basis_set_list)
538 6 : DEALLOCATE (work)
539 24 : DO i = 1, nm
540 24 : NULLIFY (mint(i)%block)
541 : END DO
542 6 : DEALLOCATE (mint)
543 :
544 6 : CALL timestop(handle)
545 :
546 12 : END SUBROUTINE build_dsdv_moments
547 :
548 : ! **************************************************************************************************
549 : !> \brief ...
550 : !> \param qs_env ...
551 : !> \param moments ...
552 : !> \param nmoments ...
553 : !> \param ref_point ...
554 : !> \param ref_points ...
555 : !> \param basis_type ...
556 : ! **************************************************************************************************
557 2482 : SUBROUTINE build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
558 :
559 : TYPE(qs_environment_type), POINTER :: qs_env
560 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moments
561 : INTEGER, INTENT(IN) :: nmoments
562 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
563 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
564 : OPTIONAL :: ref_points
565 : CHARACTER(len=*), OPTIONAL :: basis_type
566 :
567 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_local_moment_matrix'
568 :
569 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
570 : maxco, maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
571 : LOGICAL :: found
572 : REAL(KIND=dp) :: dab
573 2482 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
574 2482 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
575 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
576 2482 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
577 : TYPE(cell_type), POINTER :: cell
578 2482 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
579 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
580 : TYPE(neighbor_list_iterator_p_type), &
581 2482 : DIMENSION(:), POINTER :: nl_iterator
582 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
583 2482 : POINTER :: sab_orb
584 2482 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
585 2482 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
586 : TYPE(qs_kind_type), POINTER :: qs_kind
587 :
588 2482 : IF (nmoments < 1) RETURN
589 :
590 2482 : CALL timeset(routineN, handle)
591 :
592 2482 : NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
593 :
594 2482 : nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
595 2482 : CPASSERT(SIZE(moments) >= nm)
596 :
597 2482 : NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
598 : CALL get_qs_env(qs_env=qs_env, &
599 : qs_kind_set=qs_kind_set, &
600 : particle_set=particle_set, cell=cell, &
601 2482 : sab_orb=sab_orb)
602 :
603 2482 : nkind = SIZE(qs_kind_set)
604 2482 : natom = SIZE(particle_set)
605 :
606 : ! Allocate work storage
607 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
608 : maxco=maxco, maxsgf=maxsgf, &
609 7234 : basis_type=basis_type)
610 :
611 12410 : ALLOCATE (mab(maxco, maxco, nm))
612 2831314 : mab(:, :, :) = 0.0_dp
613 :
614 9928 : ALLOCATE (work(maxco, maxsgf))
615 456014 : work(:, :) = 0.0_dp
616 :
617 15014 : ALLOCATE (mint(nm))
618 10050 : DO i = 1, nm
619 10050 : NULLIFY (mint(i)%block)
620 : END DO
621 :
622 12510 : ALLOCATE (basis_set_list(nkind))
623 7546 : DO ikind = 1, nkind
624 5064 : qs_kind => qs_kind_set(ikind)
625 5064 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
626 7546 : IF (ASSOCIATED(basis_set_a)) THEN
627 5064 : basis_set_list(ikind)%gto_basis_set => basis_set_a
628 : ELSE
629 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
630 : END IF
631 : END DO
632 2482 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
633 30270 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
634 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
635 27788 : iatom=iatom, jatom=jatom, r=rab)
636 27788 : basis_set_a => basis_set_list(ikind)%gto_basis_set
637 27788 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
638 27788 : basis_set_b => basis_set_list(jkind)%gto_basis_set
639 27788 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
640 : ASSOCIATE ( &
641 : ! basis ikind
642 : first_sgfa => basis_set_a%first_sgf, &
643 : la_max => basis_set_a%lmax, &
644 : la_min => basis_set_a%lmin, &
645 : npgfa => basis_set_a%npgf, &
646 : nsgfa => basis_set_a%nsgf_set, &
647 : rpgfa => basis_set_a%pgf_radius, &
648 : set_radius_a => basis_set_a%set_radius, &
649 : sphi_a => basis_set_a%sphi, &
650 : zeta => basis_set_a%zet, &
651 : ! basis jkind, &
652 : first_sgfb => basis_set_b%first_sgf, &
653 : lb_max => basis_set_b%lmax, &
654 : lb_min => basis_set_b%lmin, &
655 : npgfb => basis_set_b%npgf, &
656 : nsgfb => basis_set_b%nsgf_set, &
657 : rpgfb => basis_set_b%pgf_radius, &
658 : set_radius_b => basis_set_b%set_radius, &
659 : sphi_b => basis_set_b%sphi, &
660 : zetb => basis_set_b%zet)
661 :
662 27788 : nseta = basis_set_a%nset
663 27788 : nsetb = basis_set_b%nset
664 :
665 6478 : IF (inode == 1) last_jatom = 0
666 :
667 : ! this guarantees minimum image convention
668 : ! anything else would not make sense
669 27788 : IF (jatom == last_jatom) THEN
670 : CYCLE
671 : END IF
672 :
673 7708 : last_jatom = jatom
674 :
675 7708 : IF (iatom <= jatom) THEN
676 5015 : irow = iatom
677 5015 : icol = jatom
678 : ELSE
679 2693 : irow = jatom
680 2693 : icol = iatom
681 : END IF
682 :
683 31426 : DO i = 1, nm
684 23718 : NULLIFY (mint(i)%block)
685 : CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
686 23718 : row=irow, col=icol, BLOCK=mint(i)%block, found=found)
687 1268324 : mint(i)%block = 0._dp
688 : END DO
689 :
690 : ! fold atomic position back into unit cell
691 7708 : IF (PRESENT(ref_points)) THEN
692 0 : rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
693 7708 : ELSE IF (PRESENT(ref_point)) THEN
694 28352 : rc(:) = ref_point(:)
695 : ELSE
696 620 : rc(:) = 0._dp
697 : END IF
698 : ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
699 : ! by folding around the center, such screwing can be avoided for a proper choice of center.
700 61664 : ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
701 61664 : rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
702 : ! we dont use PBC at this point
703 30832 : rab(:) = ra(:) - rb(:)
704 30832 : rac(:) = ra(:) - rc(:)
705 30832 : rbc(:) = rb(:) - rc(:)
706 30832 : dab = NORM2(rab)
707 :
708 47624 : DO iset = 1, nseta
709 :
710 12128 : ncoa = npgfa(iset)*ncoset(la_max(iset))
711 12128 : sgfa = first_sgfa(1, iset)
712 :
713 40763 : DO jset = 1, nsetb
714 :
715 20927 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
716 :
717 20846 : ncob = npgfb(jset)*ncoset(lb_max(jset))
718 20846 : sgfb = first_sgfb(1, jset)
719 :
720 : ! Calculate the primitive integrals
721 : CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
722 : rpgfa(:, iset), la_min(iset), &
723 : lb_max(jset), npgfb(jset), zetb(:, jset), &
724 20846 : rpgfb(:, jset), nmoments, rac, rbc, mab)
725 :
726 : ! Contraction step
727 97708 : DO i = 1, nm
728 :
729 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
730 : 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
731 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
732 64734 : 0.0_dp, work(1, 1), SIZE(work, 1))
733 :
734 85661 : IF (iatom <= jatom) THEN
735 :
736 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
737 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
738 : work(1, 1), SIZE(work, 1), &
739 : 1.0_dp, mint(i)%block(sgfa, sgfb), &
740 42451 : SIZE(mint(i)%block, 1))
741 :
742 : ELSE
743 :
744 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
745 : 1.0_dp, work(1, 1), SIZE(work, 1), &
746 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
747 : 1.0_dp, mint(i)%block(sgfb, sgfa), &
748 22283 : SIZE(mint(i)%block, 1))
749 :
750 : END IF
751 :
752 : END DO
753 :
754 : END DO
755 : END DO
756 : END ASSOCIATE
757 :
758 : END DO
759 2482 : CALL neighbor_list_iterator_release(nl_iterator)
760 :
761 : ! Release work storage
762 2482 : DEALLOCATE (mab, basis_set_list)
763 2482 : DEALLOCATE (work)
764 10050 : DO i = 1, nm
765 10050 : NULLIFY (mint(i)%block)
766 : END DO
767 2482 : DEALLOCATE (mint)
768 :
769 2482 : CALL timestop(handle)
770 :
771 4964 : END SUBROUTINE build_local_moment_matrix
772 :
773 : ! **************************************************************************************************
774 : !> \brief Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b >
775 : !> Optionally stores the multipole moments themselves for free.
776 : !> Note that the multipole moments are symmetric while their derivatives are anti-symmetric
777 : !> Only first derivatives are performed, e. g. x d/dy
778 : !> \param qs_env ...
779 : !> \param moments_der will contain the derivatives of the multipole moments
780 : !> \param nmoments_der order of the moments with derivatives
781 : !> \param nmoments order of the multipole moments (no derivatives, same output as
782 : !> build_local_moment_matrix, needs moments as arguments to store results)
783 : !> \param ref_point ...
784 : !> \param moments contains the multipole moments, optionally for free, up to order nmoments
785 : !> \note
786 : !> Adapted from rRc_xyz_der_ao in qs_operators_ao
787 : ! **************************************************************************************************
788 30 : SUBROUTINE build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, &
789 30 : ref_point, moments)
790 : TYPE(qs_environment_type), POINTER :: qs_env
791 : TYPE(dbcsr_p_type), DIMENSION(:, :), &
792 : INTENT(INOUT), POINTER :: moments_der
793 : INTEGER, INTENT(IN) :: nmoments_der, nmoments
794 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
795 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
796 : OPTIONAL, POINTER :: moments
797 :
798 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_local_moments_der_matrix'
799 :
800 : INTEGER :: dimders, handle, i, iatom, icol, ider, ii, ikind, inode, ipgf, irow, iset, j, &
801 : jatom, jkind, jpgf, jset, last_jatom, M_dim, maxco, maxsgf, na, nb, ncoa, ncob, nda, ndb, &
802 : nders, nkind, nm, nmom_build, nseta, nsetb, sgfa, sgfb
803 : LOGICAL :: found
804 : REAL(KIND=dp) :: dab
805 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
806 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
807 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
808 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
809 30 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab_tmp
810 30 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mom_block
811 30 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mom_block_der
812 : TYPE(cell_type), POINTER :: cell
813 30 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
814 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
815 : TYPE(neighbor_list_iterator_p_type), &
816 30 : DIMENSION(:), POINTER :: nl_iterator
817 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
818 30 : POINTER :: sab_orb
819 30 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
820 30 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
821 : TYPE(qs_kind_type), POINTER :: qs_kind
822 :
823 30 : nmom_build = MAX(nmoments, nmoments_der) ! build moments up to order nmom_buiod
824 30 : IF (nmom_build < 1) RETURN
825 :
826 30 : CALL timeset(routineN, handle)
827 :
828 30 : nders = 1 ! only first order derivatives
829 30 : dimders = ncoset(nders) - 1
830 :
831 30 : NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
832 : CALL get_qs_env(qs_env=qs_env, &
833 : qs_kind_set=qs_kind_set, &
834 : particle_set=particle_set, &
835 : cell=cell, &
836 30 : sab_orb=sab_orb)
837 :
838 30 : nkind = SIZE(qs_kind_set)
839 :
840 : ! Work storage
841 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
842 30 : maxco=maxco, maxsgf=maxsgf)
843 :
844 30 : IF (nmoments > 0) THEN
845 28 : CPASSERT(PRESENT(moments))
846 28 : nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
847 28 : CPASSERT(SIZE(moments) == nm)
848 : ! storage for integrals
849 140 : ALLOCATE (mab(maxco, maxco, nm))
850 : ! blocks
851 272620 : mab(:, :, :) = 0.0_dp
852 336 : ALLOCATE (mom_block(nm))
853 280 : DO i = 1, nm
854 280 : NULLIFY (mom_block(i)%block)
855 : END DO
856 : END IF
857 :
858 30 : IF (nmoments_der > 0) THEN
859 30 : M_dim = ncoset(nmoments_der) - 1
860 30 : CPASSERT(SIZE(moments_der, dim=1) == M_dim)
861 30 : CPASSERT(SIZE(moments_der, dim=2) == dimders)
862 : ! storage for integrals
863 180 : ALLOCATE (difmab(maxco, maxco, M_dim, dimders))
864 287454 : difmab(:, :, :, :) = 0.0_dp
865 : ! blocks
866 516 : ALLOCATE (mom_block_der(M_dim, dimders))
867 132 : DO i = 1, M_dim
868 438 : DO ider = 1, dimders
869 408 : NULLIFY (mom_block_der(i, ider)%block)
870 : END DO
871 : END DO
872 : END IF
873 :
874 120 : ALLOCATE (work(maxco, maxsgf))
875 6254 : work(:, :) = 0.0_dp
876 :
877 30 : NULLIFY (basis_set_a, basis_set_b, basis_set_list)
878 30 : NULLIFY (qs_kind)
879 128 : ALLOCATE (basis_set_list(nkind))
880 68 : DO ikind = 1, nkind
881 38 : qs_kind => qs_kind_set(ikind)
882 38 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
883 68 : IF (ASSOCIATED(basis_set_a)) THEN
884 38 : basis_set_list(ikind)%gto_basis_set => basis_set_a
885 : ELSE
886 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
887 : END IF
888 : END DO
889 :
890 : ! Calculate derivatives looping over neighbour list
891 30 : NULLIFY (nl_iterator)
892 30 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
893 213 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
894 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
895 183 : iatom=iatom, jatom=jatom, r=rab)
896 183 : basis_set_a => basis_set_list(ikind)%gto_basis_set
897 183 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
898 183 : basis_set_b => basis_set_list(jkind)%gto_basis_set
899 183 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
900 : ASSOCIATE ( &
901 : ! basis ikind
902 : first_sgfa => basis_set_a%first_sgf, &
903 : la_max => basis_set_a%lmax, &
904 : la_min => basis_set_a%lmin, &
905 : npgfa => basis_set_a%npgf, &
906 : nsgfa => basis_set_a%nsgf_set, &
907 : rpgfa => basis_set_a%pgf_radius, &
908 : set_radius_a => basis_set_a%set_radius, &
909 : sphi_a => basis_set_a%sphi, &
910 : zeta => basis_set_a%zet, &
911 : ! basis jkind, &
912 : first_sgfb => basis_set_b%first_sgf, &
913 : lb_max => basis_set_b%lmax, &
914 : lb_min => basis_set_b%lmin, &
915 : npgfb => basis_set_b%npgf, &
916 : nsgfb => basis_set_b%nsgf_set, &
917 : rpgfb => basis_set_b%pgf_radius, &
918 : set_radius_b => basis_set_b%set_radius, &
919 : sphi_b => basis_set_b%sphi, &
920 : zetb => basis_set_b%zet)
921 :
922 183 : nseta = basis_set_a%nset
923 183 : nsetb = basis_set_b%nset
924 :
925 : ! reference point
926 183 : IF (PRESENT(ref_point)) THEN
927 732 : rc(:) = ref_point(:)
928 : ELSE
929 0 : rc(:) = 0._dp
930 : END IF
931 : ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
932 : ! by folding around the center, such screwing can be avoided for a proper choice of center.
933 1464 : ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
934 1464 : rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
935 : ! we dont use PBC at this point
936 732 : rab(:) = ra(:) - rb(:)
937 732 : rac(:) = ra(:) - rc(:)
938 732 : rbc(:) = rb(:) - rc(:)
939 732 : dab = NORM2(rab)
940 :
941 : ! get blocks
942 183 : IF (inode == 1) last_jatom = 0
943 :
944 183 : IF (jatom == last_jatom) THEN
945 : CYCLE
946 : END IF
947 :
948 57 : last_jatom = jatom
949 :
950 57 : IF (iatom <= jatom) THEN
951 38 : irow = iatom
952 38 : icol = jatom
953 : ELSE
954 19 : irow = jatom
955 19 : icol = iatom
956 : END IF
957 :
958 57 : IF (nmoments > 0) THEN
959 510 : DO i = 1, nm
960 459 : NULLIFY (mom_block(i)%block)
961 : ! get block from pre calculated overlap matrix
962 : CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
963 459 : row=irow, col=icol, BLOCK=mom_block(i)%block, found=found)
964 459 : CPASSERT(found .AND. ASSOCIATED(mom_block(i)%block))
965 21003 : mom_block(i)%block = 0._dp
966 : END DO
967 : END IF
968 57 : IF (nmoments_der > 0) THEN
969 264 : DO i = 1, M_dim
970 885 : DO ider = 1, dimders
971 621 : NULLIFY (mom_block_der(i, ider)%block)
972 : CALL dbcsr_get_block_p(matrix=moments_der(i, ider)%matrix, &
973 : row=irow, col=icol, &
974 : block=mom_block_der(i, ider)%block, &
975 621 : found=found)
976 621 : CPASSERT(found .AND. ASSOCIATED(mom_block_der(i, ider)%block))
977 22455 : mom_block_der(i, ider)%block = 0._dp
978 : END DO
979 : END DO
980 : END IF
981 :
982 330 : DO iset = 1, nseta
983 :
984 90 : ncoa = npgfa(iset)*ncoset(la_max(iset))
985 90 : sgfa = first_sgfa(1, iset)
986 :
987 303 : DO jset = 1, nsetb
988 :
989 156 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
990 :
991 156 : ncob = npgfb(jset)*ncoset(lb_max(jset))
992 156 : sgfb = first_sgfb(1, jset)
993 :
994 156 : NULLIFY (mab_tmp)
995 : ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
996 780 : npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(nmom_build) - 1))
997 :
998 : ! Calculate the primitive integrals (need l+1 for derivatives)
999 : CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), &
1000 : rpgfa(:, iset), la_min(iset), &
1001 : lb_max(jset) + 1, npgfb(jset), zetb(:, jset), &
1002 156 : rpgfb(:, jset), nmom_build, rac, rbc, mab_tmp)
1003 :
1004 156 : IF (nmoments_der > 0) THEN
1005 : CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), &
1006 : rpgfa(:, iset), la_min(iset), &
1007 : lb_max(jset), npgfb(jset), zetb(:, jset), &
1008 : rpgfb(:, jset), lb_min(jset), &
1009 156 : nmoments_der, rac, rbc, difmab, mab_ext=mab_tmp)
1010 : END IF
1011 :
1012 156 : IF (nmoments > 0) THEN
1013 : ! copy subset of mab_tmp (l+1) to mab (l)
1014 830400 : mab = 0.0_dp
1015 1500 : DO ii = 1, nm
1016 1350 : na = 0
1017 1350 : nda = 0
1018 5604 : DO ipgf = 1, npgfa(iset)
1019 4104 : nb = 0
1020 4104 : ndb = 0
1021 19467 : DO jpgf = 1, npgfb(jset)
1022 74871 : DO j = 1, ncoset(lb_max(jset))
1023 395523 : DO i = 1, ncoset(la_max(iset))
1024 380160 : mab(i + na, j + nb, ii) = mab_tmp(i + nda, j + ndb, ii)
1025 : END DO ! i
1026 : END DO ! j
1027 15363 : nb = nb + ncoset(lb_max(jset))
1028 19467 : ndb = ndb + ncoset(lb_max(jset) + 1)
1029 : END DO ! jpgf
1030 4104 : na = na + ncoset(la_max(iset))
1031 5454 : nda = nda + ncoset(la_max(iset) + 1)
1032 : END DO ! ipgf
1033 : END DO
1034 : ! Contraction step
1035 1500 : DO i = 1, nm
1036 :
1037 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1038 : 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
1039 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1040 1350 : 0.0_dp, work(1, 1), SIZE(work, 1))
1041 :
1042 1500 : IF (iatom <= jatom) THEN
1043 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1044 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1045 : work(1, 1), SIZE(work, 1), &
1046 : 1.0_dp, mom_block(i)%block(sgfa, sgfb), &
1047 900 : SIZE(mom_block(i)%block, 1))
1048 : ELSE
1049 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1050 : 1.0_dp, work(1, 1), SIZE(work, 1), &
1051 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1052 : 1.0_dp, mom_block(i)%block(sgfb, sgfa), &
1053 450 : SIZE(mom_block(i)%block, 1))
1054 : END IF
1055 : END DO
1056 : END IF
1057 :
1058 156 : IF (nmoments_der > 0) THEN
1059 660 : DO i = 1, M_dim
1060 2172 : DO ider = 1, dimders
1061 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1062 : 1.0_dp, difmab(1, 1, i, ider), SIZE(difmab, 1), &
1063 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1064 1512 : 0._dp, work(1, 1), SIZE(work, 1))
1065 :
1066 2016 : IF (iatom <= jatom) THEN
1067 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1068 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1069 : work(1, 1), SIZE(work, 1), &
1070 : 1._dp, mom_block_der(i, ider)%block(sgfa, sgfb), &
1071 1008 : SIZE(mom_block_der(i, ider)%block, 1))
1072 : ELSE
1073 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1074 : -1.0_dp, work(1, 1), SIZE(work, 1), &
1075 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1076 : 1.0_dp, mom_block_der(i, ider)%block(sgfb, sgfa), &
1077 504 : SIZE(mom_block_der(i, ider)%block, 1))
1078 : END IF
1079 : END DO
1080 : END DO
1081 : END IF
1082 246 : DEALLOCATE (mab_tmp)
1083 : END DO
1084 : END DO
1085 : END ASSOCIATE
1086 : END DO
1087 30 : CALL neighbor_list_iterator_release(nl_iterator)
1088 :
1089 : ! deallocations
1090 30 : DEALLOCATE (basis_set_list)
1091 30 : DEALLOCATE (work)
1092 30 : IF (nmoments > 0) THEN
1093 28 : DEALLOCATE (mab)
1094 280 : DO i = 1, nm
1095 280 : NULLIFY (mom_block(i)%block)
1096 : END DO
1097 28 : DEALLOCATE (mom_block)
1098 : END IF
1099 30 : IF (nmoments_der > 0) THEN
1100 30 : DEALLOCATE (difmab)
1101 132 : DO i = 1, M_dim
1102 438 : DO ider = 1, dimders
1103 408 : NULLIFY (mom_block_der(i, ider)%block)
1104 : END DO
1105 : END DO
1106 30 : DEALLOCATE (mom_block_der)
1107 : END IF
1108 :
1109 30 : CALL timestop(handle)
1110 :
1111 90 : END SUBROUTINE build_local_moments_der_matrix
1112 :
1113 : ! **************************************************************************************************
1114 : !> \brief ...
1115 : !> \param qs_env ...
1116 : !> \param magmom ...
1117 : !> \param nmoments ...
1118 : !> \param ref_point ...
1119 : !> \param ref_points ...
1120 : !> \param basis_type ...
1121 : ! **************************************************************************************************
1122 16 : SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
1123 :
1124 : TYPE(qs_environment_type), POINTER :: qs_env
1125 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom
1126 : INTEGER, INTENT(IN) :: nmoments
1127 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
1128 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
1129 : OPTIONAL :: ref_points
1130 : CHARACTER(len=*), OPTIONAL :: basis_type
1131 :
1132 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_local_magmom_matrix'
1133 :
1134 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, maxco, &
1135 : maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
1136 : LOGICAL :: found
1137 : REAL(KIND=dp) :: dab
1138 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1139 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
1140 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
1141 16 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
1142 : TYPE(cell_type), POINTER :: cell
1143 16 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1144 16 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1145 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1146 : TYPE(neighbor_list_iterator_p_type), &
1147 16 : DIMENSION(:), POINTER :: nl_iterator
1148 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1149 16 : POINTER :: sab_orb
1150 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1151 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1152 : TYPE(qs_kind_type), POINTER :: qs_kind
1153 :
1154 16 : IF (nmoments < 1) RETURN
1155 :
1156 16 : CALL timeset(routineN, handle)
1157 :
1158 16 : NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
1159 16 : NULLIFY (matrix_s)
1160 :
1161 16 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
1162 :
1163 : ! magnetic dipoles/angular moments only
1164 16 : nm = 3
1165 :
1166 16 : NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
1167 : CALL get_qs_env(qs_env=qs_env, &
1168 : qs_kind_set=qs_kind_set, &
1169 : particle_set=particle_set, cell=cell, &
1170 16 : sab_orb=sab_orb)
1171 :
1172 16 : nkind = SIZE(qs_kind_set)
1173 16 : natom = SIZE(particle_set)
1174 :
1175 : ! Allocate work storage
1176 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1177 16 : maxco=maxco, maxsgf=maxsgf)
1178 :
1179 80 : ALLOCATE (mab(maxco, maxco, nm))
1180 210436 : mab(:, :, :) = 0.0_dp
1181 :
1182 64 : ALLOCATE (work(maxco, maxsgf))
1183 13380 : work(:, :) = 0.0_dp
1184 :
1185 64 : ALLOCATE (mint(nm))
1186 64 : DO i = 1, nm
1187 64 : NULLIFY (mint(i)%block)
1188 : END DO
1189 :
1190 80 : ALLOCATE (basis_set_list(nkind))
1191 48 : DO ikind = 1, nkind
1192 32 : qs_kind => qs_kind_set(ikind)
1193 96 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1194 48 : IF (ASSOCIATED(basis_set_a)) THEN
1195 32 : basis_set_list(ikind)%gto_basis_set => basis_set_a
1196 : ELSE
1197 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
1198 : END IF
1199 : END DO
1200 16 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1201 175 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1202 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1203 159 : iatom=iatom, jatom=jatom, r=rab)
1204 159 : basis_set_a => basis_set_list(ikind)%gto_basis_set
1205 159 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
1206 159 : basis_set_b => basis_set_list(jkind)%gto_basis_set
1207 159 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
1208 : ASSOCIATE ( &
1209 : ! basis ikind
1210 : first_sgfa => basis_set_a%first_sgf, &
1211 : la_max => basis_set_a%lmax, &
1212 : la_min => basis_set_a%lmin, &
1213 : npgfa => basis_set_a%npgf, &
1214 : nsgfa => basis_set_a%nsgf_set, &
1215 : rpgfa => basis_set_a%pgf_radius, &
1216 : set_radius_a => basis_set_a%set_radius, &
1217 : sphi_a => basis_set_a%sphi, &
1218 : zeta => basis_set_a%zet, &
1219 : ! basis jkind, &
1220 : first_sgfb => basis_set_b%first_sgf, &
1221 : lb_max => basis_set_b%lmax, &
1222 : lb_min => basis_set_b%lmin, &
1223 : npgfb => basis_set_b%npgf, &
1224 : nsgfb => basis_set_b%nsgf_set, &
1225 : rpgfb => basis_set_b%pgf_radius, &
1226 : set_radius_b => basis_set_b%set_radius, &
1227 : sphi_b => basis_set_b%sphi, &
1228 : zetb => basis_set_b%zet)
1229 :
1230 159 : nseta = basis_set_a%nset
1231 159 : nsetb = basis_set_b%nset
1232 :
1233 159 : IF (iatom <= jatom) THEN
1234 106 : irow = iatom
1235 106 : icol = jatom
1236 : ELSE
1237 53 : irow = jatom
1238 53 : icol = iatom
1239 : END IF
1240 :
1241 636 : DO i = 1, nm
1242 477 : NULLIFY (mint(i)%block)
1243 : CALL dbcsr_get_block_p(matrix=magmom(i)%matrix, &
1244 477 : row=irow, col=icol, BLOCK=mint(i)%block, found=found)
1245 33123 : mint(i)%block = 0._dp
1246 1113 : CPASSERT(ASSOCIATED(mint(i)%block))
1247 : END DO
1248 :
1249 : ! fold atomic position back into unit cell
1250 159 : IF (PRESENT(ref_points)) THEN
1251 0 : rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
1252 159 : ELSE IF (PRESENT(ref_point)) THEN
1253 636 : rc(:) = ref_point(:)
1254 : ELSE
1255 0 : rc(:) = 0._dp
1256 : END IF
1257 : ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
1258 : ! by folding around the center, such screwing can be avoided for a proper choice of center.
1259 1272 : ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
1260 1272 : rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
1261 : ! we dont use PBC at this point
1262 636 : rab(:) = ra(:) - rb(:)
1263 636 : rac(:) = ra(:) - rc(:)
1264 636 : rbc(:) = rb(:) - rc(:)
1265 636 : dab = NORM2(rab)
1266 :
1267 594 : DO iset = 1, nseta
1268 :
1269 276 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1270 276 : sgfa = first_sgfa(1, iset)
1271 :
1272 945 : DO jset = 1, nsetb
1273 :
1274 510 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
1275 :
1276 510 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1277 510 : sgfb = first_sgfb(1, jset)
1278 :
1279 : ! Calculate the primitive integrals
1280 : CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), &
1281 : rpgfa(:, iset), la_min(iset), &
1282 : lb_max(jset), npgfb(jset), zetb(:, jset), &
1283 510 : rpgfb(:, jset), rac, rbc, mab)
1284 :
1285 : ! Contraction step
1286 2316 : DO i = 1, nm
1287 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1288 : 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
1289 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1290 1530 : 0.0_dp, work(1, 1), SIZE(work, 1))
1291 :
1292 2040 : IF (iatom <= jatom) THEN
1293 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1294 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1295 : work(1, 1), SIZE(work, 1), &
1296 : 1.0_dp, mint(i)%block(sgfa, sgfb), &
1297 1020 : SIZE(mint(i)%block, 1))
1298 : ELSE
1299 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1300 : -1.0_dp, work(1, 1), SIZE(work, 1), &
1301 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1302 : 1.0_dp, mint(i)%block(sgfb, sgfa), &
1303 510 : SIZE(mint(i)%block, 1))
1304 : END IF
1305 :
1306 : END DO
1307 :
1308 : END DO
1309 : END DO
1310 : END ASSOCIATE
1311 : END DO
1312 16 : CALL neighbor_list_iterator_release(nl_iterator)
1313 :
1314 : ! Release work storage
1315 16 : DEALLOCATE (mab, basis_set_list)
1316 16 : DEALLOCATE (work)
1317 64 : DO i = 1, nm
1318 64 : NULLIFY (mint(i)%block)
1319 : END DO
1320 16 : DEALLOCATE (mint)
1321 :
1322 16 : CALL timestop(handle)
1323 :
1324 48 : END SUBROUTINE build_local_magmom_matrix
1325 :
1326 : ! **************************************************************************************************
1327 : !> \brief ...
1328 : !> \param qs_env ...
1329 : !> \param cosmat ...
1330 : !> \param sinmat ...
1331 : !> \param kvec ...
1332 : !> \param sab_orb_external ...
1333 : !> \param basis_type ...
1334 : ! **************************************************************************************************
1335 8660 : SUBROUTINE build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
1336 :
1337 : TYPE(qs_environment_type), POINTER :: qs_env
1338 : TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1339 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: kvec
1340 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1341 : OPTIONAL, POINTER :: sab_orb_external
1342 : CHARACTER(len=*), OPTIONAL :: basis_type
1343 :
1344 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_berry_moment_matrix'
1345 :
1346 : INTEGER :: handle, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, ldsa, &
1347 : ldsb, ldwork, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
1348 : LOGICAL :: found
1349 8660 : REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1350 : REAL(KIND=dp) :: dab
1351 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb
1352 : TYPE(cell_type), POINTER :: cell
1353 8660 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1354 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1355 : TYPE(neighbor_list_iterator_p_type), &
1356 8660 : DIMENSION(:), POINTER :: nl_iterator
1357 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1358 8660 : POINTER :: sab_orb
1359 8660 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1360 8660 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1361 : TYPE(qs_kind_type), POINTER :: qs_kind
1362 :
1363 8660 : CALL timeset(routineN, handle)
1364 :
1365 8660 : NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
1366 : CALL get_qs_env(qs_env=qs_env, &
1367 : qs_kind_set=qs_kind_set, &
1368 : particle_set=particle_set, cell=cell, &
1369 8660 : sab_orb=sab_orb)
1370 :
1371 8660 : IF (PRESENT(sab_orb_external)) sab_orb => sab_orb_external
1372 :
1373 8660 : CALL dbcsr_set(sinmat, 0.0_dp)
1374 8660 : CALL dbcsr_set(cosmat, 0.0_dp)
1375 :
1376 12876 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork, basis_type=basis_type)
1377 8660 : ldab = ldwork
1378 34640 : ALLOCATE (cosab(ldab, ldab))
1379 34640 : ALLOCATE (sinab(ldab, ldab))
1380 34640 : ALLOCATE (work(ldwork, ldwork))
1381 :
1382 8660 : nkind = SIZE(qs_kind_set)
1383 8660 : natom = SIZE(particle_set)
1384 :
1385 43360 : ALLOCATE (basis_set_list(nkind))
1386 26040 : DO ikind = 1, nkind
1387 17380 : qs_kind => qs_kind_set(ikind)
1388 17380 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1389 26040 : IF (ASSOCIATED(basis_set_a)) THEN
1390 17380 : basis_set_list(ikind)%gto_basis_set => basis_set_a
1391 : ELSE
1392 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
1393 : END IF
1394 : END DO
1395 8660 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1396 385932 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1397 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1398 377272 : iatom=iatom, jatom=jatom, r=rab)
1399 377272 : basis_set_a => basis_set_list(ikind)%gto_basis_set
1400 377272 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
1401 377272 : basis_set_b => basis_set_list(jkind)%gto_basis_set
1402 377272 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
1403 : ASSOCIATE ( &
1404 : ! basis ikind
1405 : first_sgfa => basis_set_a%first_sgf, &
1406 : la_max => basis_set_a%lmax, &
1407 : la_min => basis_set_a%lmin, &
1408 : npgfa => basis_set_a%npgf, &
1409 : nsgfa => basis_set_a%nsgf_set, &
1410 : rpgfa => basis_set_a%pgf_radius, &
1411 : set_radius_a => basis_set_a%set_radius, &
1412 : sphi_a => basis_set_a%sphi, &
1413 : zeta => basis_set_a%zet, &
1414 : ! basis jkind, &
1415 : first_sgfb => basis_set_b%first_sgf, &
1416 : lb_max => basis_set_b%lmax, &
1417 : lb_min => basis_set_b%lmin, &
1418 : npgfb => basis_set_b%npgf, &
1419 : nsgfb => basis_set_b%nsgf_set, &
1420 : rpgfb => basis_set_b%pgf_radius, &
1421 : set_radius_b => basis_set_b%set_radius, &
1422 : sphi_b => basis_set_b%sphi, &
1423 : zetb => basis_set_b%zet)
1424 :
1425 377272 : nseta = basis_set_a%nset
1426 377272 : nsetb = basis_set_b%nset
1427 :
1428 377272 : ldsa = SIZE(sphi_a, 1)
1429 377272 : ldsb = SIZE(sphi_b, 1)
1430 :
1431 377272 : IF (iatom <= jatom) THEN
1432 212258 : irow = iatom
1433 212258 : icol = jatom
1434 : ELSE
1435 165014 : irow = jatom
1436 165014 : icol = iatom
1437 : END IF
1438 :
1439 377272 : NULLIFY (cblock)
1440 : CALL dbcsr_get_block_p(matrix=cosmat, &
1441 377272 : row=irow, col=icol, BLOCK=cblock, found=found)
1442 377272 : NULLIFY (sblock)
1443 : CALL dbcsr_get_block_p(matrix=sinmat, &
1444 377272 : row=irow, col=icol, BLOCK=sblock, found=found)
1445 377272 : IF (ASSOCIATED(cblock) .AND. .NOT. ASSOCIATED(sblock) .OR. &
1446 377272 : .NOT. ASSOCIATED(cblock) .AND. ASSOCIATED(sblock)) THEN
1447 0 : CPABORT("")
1448 : END IF
1449 :
1450 1131816 : IF (ASSOCIATED(cblock) .AND. ASSOCIATED(sblock)) THEN
1451 :
1452 377272 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
1453 1509088 : rb(:) = ra + rab
1454 377272 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1455 :
1456 1317470 : DO iset = 1, nseta
1457 :
1458 940198 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1459 940198 : sgfa = first_sgfa(1, iset)
1460 :
1461 4583145 : DO jset = 1, nsetb
1462 :
1463 3265675 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
1464 :
1465 1462057 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1466 1462057 : sgfb = first_sgfb(1, jset)
1467 :
1468 : ! Calculate the primitive integrals
1469 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1470 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1471 1462057 : ra, rb, kvec, cosab, sinab)
1472 : CALL contract_cossin(cblock, sblock, &
1473 : iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1474 : jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1475 4205873 : cosab, sinab, ldab, work, ldwork)
1476 :
1477 : END DO
1478 : END DO
1479 :
1480 : END IF
1481 : END ASSOCIATE
1482 : END DO
1483 8660 : CALL neighbor_list_iterator_release(nl_iterator)
1484 :
1485 8660 : DEALLOCATE (cosab)
1486 8660 : DEALLOCATE (sinab)
1487 8660 : DEALLOCATE (work)
1488 8660 : DEALLOCATE (basis_set_list)
1489 :
1490 8660 : CALL timestop(handle)
1491 :
1492 8660 : END SUBROUTINE build_berry_moment_matrix
1493 :
1494 : ! **************************************************************************************************
1495 : !> \brief ...
1496 : !> \param qs_env ...
1497 : !> \param cosmat ...
1498 : !> \param sinmat ...
1499 : !> \param kvec ...
1500 : !> \param basis_type ...
1501 : ! **************************************************************************************************
1502 0 : SUBROUTINE build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
1503 :
1504 : TYPE(qs_environment_type), POINTER :: qs_env
1505 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: cosmat, sinmat
1506 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: kvec
1507 : CHARACTER(len=*), OPTIONAL :: basis_type
1508 :
1509 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_berry_kpoint_matrix'
1510 :
1511 : INTEGER :: handle, i, iatom, ic, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, &
1512 : ldsa, ldsb, ldwork, natom, ncoa, ncob, nimg, nkind, nseta, nsetb, sgfa, sgfb
1513 : INTEGER, DIMENSION(3) :: icell
1514 0 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
1515 0 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1516 : LOGICAL :: found, use_cell_mapping
1517 0 : REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1518 : REAL(KIND=dp) :: dab
1519 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb
1520 : TYPE(cell_type), POINTER :: cell
1521 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1522 : TYPE(dft_control_type), POINTER :: dft_control
1523 0 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1524 : TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
1525 : TYPE(kpoint_type), POINTER :: kpoints
1526 : TYPE(neighbor_list_iterator_p_type), &
1527 0 : DIMENSION(:), POINTER :: nl_iterator
1528 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1529 0 : POINTER :: sab_orb
1530 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1531 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1532 : TYPE(qs_kind_type), POINTER :: qs_kind
1533 : TYPE(qs_ks_env_type), POINTER :: ks_env
1534 :
1535 0 : CALL timeset(routineN, handle)
1536 :
1537 : CALL get_qs_env(qs_env, &
1538 : ks_env=ks_env, &
1539 0 : dft_control=dft_control)
1540 0 : nimg = dft_control%nimages
1541 0 : IF (nimg > 1) THEN
1542 0 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1543 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1544 0 : use_cell_mapping = .TRUE.
1545 : ELSE
1546 : use_cell_mapping = .FALSE.
1547 : END IF
1548 :
1549 : CALL get_qs_env(qs_env=qs_env, &
1550 : qs_kind_set=qs_kind_set, &
1551 : particle_set=particle_set, cell=cell, &
1552 0 : sab_orb=sab_orb)
1553 :
1554 0 : nkind = SIZE(qs_kind_set)
1555 0 : natom = SIZE(particle_set)
1556 0 : ALLOCATE (basis_set_list(nkind))
1557 0 : DO ikind = 1, nkind
1558 0 : qs_kind => qs_kind_set(ikind)
1559 0 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=basis_type)
1560 0 : IF (ASSOCIATED(basis_set)) THEN
1561 0 : basis_set_list(ikind)%gto_basis_set => basis_set
1562 : ELSE
1563 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
1564 : END IF
1565 : END DO
1566 :
1567 0 : ALLOCATE (row_blk_sizes(natom))
1568 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1569 0 : basis=basis_set_list)
1570 0 : CALL get_ks_env(ks_env, dbcsr_dist=dbcsr_dist)
1571 : ! (re)allocate matrix sets
1572 0 : CALL dbcsr_allocate_matrix_set(sinmat, 1, nimg)
1573 0 : CALL dbcsr_allocate_matrix_set(cosmat, 1, nimg)
1574 0 : DO i = 1, nimg
1575 : ! sin
1576 0 : ALLOCATE (sinmat(1, i)%matrix)
1577 : CALL dbcsr_create(matrix=sinmat(1, i)%matrix, &
1578 : name="SINMAT", &
1579 : dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1580 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
1581 0 : nze=0)
1582 0 : CALL cp_dbcsr_alloc_block_from_nbl(sinmat(1, i)%matrix, sab_orb)
1583 0 : CALL dbcsr_set(sinmat(1, i)%matrix, 0.0_dp)
1584 : ! cos
1585 0 : ALLOCATE (cosmat(1, i)%matrix)
1586 : CALL dbcsr_create(matrix=cosmat(1, i)%matrix, &
1587 : name="COSMAT", &
1588 : dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1589 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
1590 0 : nze=0)
1591 0 : CALL cp_dbcsr_alloc_block_from_nbl(cosmat(1, i)%matrix, sab_orb)
1592 0 : CALL dbcsr_set(cosmat(1, i)%matrix, 0.0_dp)
1593 : END DO
1594 :
1595 0 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
1596 0 : ldab = ldwork
1597 0 : ALLOCATE (cosab(ldab, ldab))
1598 0 : ALLOCATE (sinab(ldab, ldab))
1599 0 : ALLOCATE (work(ldwork, ldwork))
1600 :
1601 0 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1602 0 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1603 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1604 0 : iatom=iatom, jatom=jatom, r=rab, cell=icell)
1605 0 : basis_set_a => basis_set_list(ikind)%gto_basis_set
1606 0 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
1607 0 : basis_set_b => basis_set_list(jkind)%gto_basis_set
1608 0 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
1609 : ASSOCIATE ( &
1610 : ! basis ikind
1611 : first_sgfa => basis_set_a%first_sgf, &
1612 : la_max => basis_set_a%lmax, &
1613 : la_min => basis_set_a%lmin, &
1614 : npgfa => basis_set_a%npgf, &
1615 : nsgfa => basis_set_a%nsgf_set, &
1616 : rpgfa => basis_set_a%pgf_radius, &
1617 : set_radius_a => basis_set_a%set_radius, &
1618 : sphi_a => basis_set_a%sphi, &
1619 : zeta => basis_set_a%zet, &
1620 : ! basis jkind, &
1621 : first_sgfb => basis_set_b%first_sgf, &
1622 : lb_max => basis_set_b%lmax, &
1623 : lb_min => basis_set_b%lmin, &
1624 : npgfb => basis_set_b%npgf, &
1625 : nsgfb => basis_set_b%nsgf_set, &
1626 : rpgfb => basis_set_b%pgf_radius, &
1627 : set_radius_b => basis_set_b%set_radius, &
1628 : sphi_b => basis_set_b%sphi, &
1629 0 : zetb => basis_set_b%zet)
1630 :
1631 0 : nseta = basis_set_a%nset
1632 0 : nsetb = basis_set_b%nset
1633 :
1634 0 : ldsa = SIZE(sphi_a, 1)
1635 0 : ldsb = SIZE(sphi_b, 1)
1636 :
1637 0 : IF (iatom <= jatom) THEN
1638 0 : irow = iatom
1639 0 : icol = jatom
1640 : ELSE
1641 0 : irow = jatom
1642 0 : icol = iatom
1643 : END IF
1644 :
1645 0 : IF (use_cell_mapping) THEN
1646 0 : ic = cell_to_index(icell(1), icell(2), icell(3))
1647 0 : CPASSERT(ic > 0)
1648 : ELSE
1649 : ic = 1
1650 : END IF
1651 :
1652 0 : NULLIFY (sblock)
1653 : CALL dbcsr_get_block_p(matrix=sinmat(1, ic)%matrix, &
1654 0 : row=irow, col=icol, BLOCK=sblock, found=found)
1655 0 : CPASSERT(found)
1656 0 : NULLIFY (cblock)
1657 : CALL dbcsr_get_block_p(matrix=cosmat(1, ic)%matrix, &
1658 0 : row=irow, col=icol, BLOCK=cblock, found=found)
1659 0 : CPASSERT(found)
1660 :
1661 0 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
1662 0 : rb(:) = ra + rab
1663 0 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1664 :
1665 0 : DO iset = 1, nseta
1666 :
1667 0 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1668 0 : sgfa = first_sgfa(1, iset)
1669 :
1670 0 : DO jset = 1, nsetb
1671 :
1672 0 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
1673 :
1674 0 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1675 0 : sgfb = first_sgfb(1, jset)
1676 :
1677 : ! Calculate the primitive integrals
1678 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1679 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1680 0 : ra, rb, kvec, cosab, sinab)
1681 : CALL contract_cossin(cblock, sblock, &
1682 : iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1683 : jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1684 0 : cosab, sinab, ldab, work, ldwork)
1685 :
1686 : END DO
1687 : END DO
1688 : END ASSOCIATE
1689 : END DO
1690 0 : CALL neighbor_list_iterator_release(nl_iterator)
1691 :
1692 0 : DEALLOCATE (cosab)
1693 0 : DEALLOCATE (sinab)
1694 0 : DEALLOCATE (work)
1695 0 : DEALLOCATE (basis_set_list)
1696 0 : DEALLOCATE (row_blk_sizes)
1697 :
1698 0 : CALL timestop(handle)
1699 :
1700 0 : END SUBROUTINE build_berry_kpoint_matrix
1701 :
1702 : ! **************************************************************************************************
1703 : !> \brief ...
1704 : !> \param qs_env ...
1705 : !> \param magnetic ...
1706 : !> \param nmoments ...
1707 : !> \param reference ...
1708 : !> \param ref_point ...
1709 : !> \param unit_number ...
1710 : ! **************************************************************************************************
1711 340 : SUBROUTINE qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number)
1712 :
1713 : TYPE(qs_environment_type), POINTER :: qs_env
1714 : LOGICAL, INTENT(IN) :: magnetic
1715 : INTEGER, INTENT(IN) :: nmoments, reference
1716 : REAL(dp), DIMENSION(:), POINTER :: ref_point
1717 : INTEGER, INTENT(IN) :: unit_number
1718 :
1719 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_moment_berry_phase'
1720 :
1721 340 : CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
1722 : CHARACTER(LEN=default_string_length) :: description
1723 : COMPLEX(dp) :: xphase(3), zdet, zdeta, zi(3), &
1724 : zij(3, 3), zijk(3, 3, 3), &
1725 : zijkl(3, 3, 3, 3), zphase(3), zz
1726 : INTEGER :: handle, i, ia, idim, ikind, ispin, ix, &
1727 : iy, iz, j, k, l, nao, nm, nmo, nmom, &
1728 : nmotot, tmp_dim
1729 : LOGICAL :: floating, ghost, uniform
1730 : REAL(dp) :: charge, ci(3), cij(3, 3), dd, occ, trace
1731 340 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom
1732 340 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
1733 : REAL(dp), DIMENSION(3) :: kvec, qq, rcc, ria
1734 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1735 : TYPE(cell_type), POINTER :: cell
1736 340 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat
1737 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
1738 340 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: opvec
1739 340 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: op_fm_set
1740 340 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
1741 : TYPE(cp_fm_type), POINTER :: mo_coeff
1742 : TYPE(cp_result_type), POINTER :: results
1743 340 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
1744 : TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1745 : TYPE(dft_control_type), POINTER :: dft_control
1746 : TYPE(distribution_1d_type), POINTER :: local_particles
1747 340 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1748 : TYPE(mp_para_env_type), POINTER :: para_env
1749 340 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1750 340 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1751 : TYPE(qs_rho_type), POINTER :: rho
1752 : TYPE(rt_prop_type), POINTER :: rtp
1753 :
1754 0 : CPASSERT(ASSOCIATED(qs_env))
1755 :
1756 340 : IF (ASSOCIATED(qs_env%ls_scf_env)) THEN
1757 0 : IF (unit_number > 0) WRITE (unit_number, *) "Periodic moment calculation not implemented in linear scaling code"
1758 0 : RETURN
1759 : END IF
1760 :
1761 340 : CALL timeset(routineN, handle)
1762 :
1763 : ! restrict maximum moment available
1764 340 : nmom = MIN(nmoments, 2)
1765 :
1766 340 : nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
1767 : ! rmom(:,1)=electronic
1768 : ! rmom(:,2)=nuclear
1769 : ! rmom(:,1)=total
1770 1700 : ALLOCATE (rmom(nm + 1, 3))
1771 1020 : ALLOCATE (rlab(nm + 1))
1772 5440 : rmom = 0.0_dp
1773 1700 : rlab = ""
1774 340 : IF (magnetic) THEN
1775 0 : nm = 3
1776 0 : ALLOCATE (mmom(nm))
1777 0 : mmom = 0._dp
1778 : END IF
1779 :
1780 340 : NULLIFY (dft_control, rho, cell, particle_set, results, para_env, &
1781 340 : local_particles, matrix_s, mos, rho_ao)
1782 :
1783 : CALL get_qs_env(qs_env, &
1784 : dft_control=dft_control, &
1785 : rho=rho, &
1786 : cell=cell, &
1787 : results=results, &
1788 : particle_set=particle_set, &
1789 : qs_kind_set=qs_kind_set, &
1790 : para_env=para_env, &
1791 : local_particles=local_particles, &
1792 : matrix_s=matrix_s, &
1793 340 : mos=mos)
1794 :
1795 340 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1796 :
1797 340 : NULLIFY (cosmat, sinmat)
1798 340 : ALLOCATE (cosmat, sinmat)
1799 340 : CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
1800 340 : CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
1801 340 : CALL dbcsr_set(cosmat, 0.0_dp)
1802 340 : CALL dbcsr_set(sinmat, 0.0_dp)
1803 :
1804 2106 : ALLOCATE (op_fm_set(2, dft_control%nspins))
1805 1382 : ALLOCATE (opvec(dft_control%nspins))
1806 1382 : ALLOCATE (eigrmat(dft_control%nspins))
1807 340 : nmotot = 0
1808 702 : DO ispin = 1, dft_control%nspins
1809 362 : NULLIFY (tmp_fm_struct, mo_coeff)
1810 362 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
1811 362 : nmotot = nmotot + nmo
1812 362 : CALL cp_fm_create(opvec(ispin), mo_coeff%matrix_struct)
1813 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
1814 362 : ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
1815 1086 : DO i = 1, SIZE(op_fm_set, 1)
1816 1086 : CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
1817 : END DO
1818 362 : CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
1819 1064 : CALL cp_fm_struct_release(tmp_fm_struct)
1820 : END DO
1821 :
1822 : ! occupation
1823 702 : DO ispin = 1, dft_control%nspins
1824 362 : CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
1825 702 : IF (.NOT. uniform) THEN
1826 0 : CPWARN("Berry phase moments for non uniform MOs' occupation numbers not implemented")
1827 : END IF
1828 : END DO
1829 :
1830 : ! reference point
1831 340 : CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
1832 1360 : rcc = pbc(rcc, cell)
1833 :
1834 : ! label
1835 1360 : DO l = 1, nm
1836 1020 : ix = indco(1, l + 1)
1837 1020 : iy = indco(2, l + 1)
1838 1020 : iz = indco(3, l + 1)
1839 1360 : CALL set_label(rlab(l + 1), ix, iy, iz)
1840 : END DO
1841 :
1842 : ! nuclear contribution
1843 1380 : DO ia = 1, SIZE(particle_set)
1844 1040 : atomic_kind => particle_set(ia)%atomic_kind
1845 1040 : CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1846 1040 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1847 1380 : IF (.NOT. ghost .AND. .NOT. floating) THEN
1848 1040 : rmom(1, 2) = rmom(1, 2) - charge
1849 : END IF
1850 : END DO
1851 5440 : ria = twopi*MATMUL(cell%h_inv, rcc)
1852 1360 : zphase = CMPLX(COS(ria), SIN(ria), dp)**rmom(1, 2)
1853 :
1854 340 : zi = 0._dp
1855 340 : zij = 0._dp
1856 : zijk = 0._dp
1857 : zijkl = 0._dp
1858 :
1859 680 : DO l = 1, nmom
1860 340 : SELECT CASE (l)
1861 : CASE (1)
1862 : ! Dipole
1863 1360 : zi(:) = CMPLX(1._dp, 0._dp, dp)
1864 1380 : DO ia = 1, SIZE(particle_set)
1865 1040 : atomic_kind => particle_set(ia)%atomic_kind
1866 1040 : CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1867 1040 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1868 1380 : IF (.NOT. ghost .AND. .NOT. floating) THEN
1869 4160 : ria = particle_set(ia)%r
1870 4160 : ria = pbc(ria, cell)
1871 4160 : DO i = 1, 3
1872 12480 : kvec(:) = twopi*cell%h_inv(i, :)
1873 12480 : dd = SUM(kvec(:)*ria(:))
1874 3120 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
1875 4160 : zi(i) = zi(i)*zdeta
1876 : END DO
1877 : END IF
1878 : END DO
1879 1360 : zi = zi*zphase
1880 1360 : ci = AIMAG(LOG(zi))/twopi
1881 1360 : qq = AIMAG(LOG(zi))
1882 5440 : rmom(2:4, 2) = MATMUL(cell%hmat, ci)
1883 : CASE (2)
1884 : ! Quadrupole
1885 0 : CPABORT("Berry phase moments bigger than 1 not implemented")
1886 0 : zij(:, :) = CMPLX(1._dp, 0._dp, dp)
1887 0 : DO ia = 1, SIZE(particle_set)
1888 0 : atomic_kind => particle_set(ia)%atomic_kind
1889 0 : CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1890 0 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
1891 0 : ria = particle_set(ia)%r
1892 0 : ria = pbc(ria, cell)
1893 0 : DO i = 1, 3
1894 0 : DO j = i, 3
1895 0 : kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1896 0 : dd = SUM(kvec(:)*ria(:))
1897 0 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
1898 0 : zij(i, j) = zij(i, j)*zdeta
1899 0 : zij(j, i) = zij(i, j)
1900 : END DO
1901 : END DO
1902 : END DO
1903 0 : DO i = 1, 3
1904 0 : DO j = 1, 3
1905 0 : zij(i, j) = zij(i, j)*zphase(i)*zphase(j)
1906 0 : zz = zij(i, j)/zi(i)/zi(j)
1907 0 : cij(i, j) = AIMAG(LOG(zz))/twopi
1908 : END DO
1909 : END DO
1910 0 : cij = 0.5_dp*cij/twopi/twopi
1911 0 : cij = MATMUL(MATMUL(cell%hmat, cij), TRANSPOSE(cell%hmat))
1912 0 : DO k = 4, 9
1913 0 : ix = indco(1, k + 1)
1914 0 : iy = indco(2, k + 1)
1915 0 : iz = indco(3, k + 1)
1916 0 : IF (ix == 0) THEN
1917 0 : rmom(k + 1, 2) = cij(iy, iz)
1918 0 : ELSE IF (iy == 0) THEN
1919 0 : rmom(k + 1, 2) = cij(ix, iz)
1920 0 : ELSE IF (iz == 0) THEN
1921 0 : rmom(k + 1, 2) = cij(ix, iy)
1922 : END IF
1923 : END DO
1924 : CASE (3)
1925 : ! Octapole
1926 0 : CPABORT("Berry phase moments bigger than 2 not implemented")
1927 : CASE (4)
1928 : ! Hexadecapole
1929 0 : CPABORT("Berry phase moments bigger than 3 not implemented")
1930 : CASE DEFAULT
1931 340 : CPABORT("Berry phase moments bigger than 4 not implemented")
1932 : END SELECT
1933 : END DO
1934 :
1935 : ! electronic contribution
1936 :
1937 5440 : ria = twopi*REAL(nmotot, dp)*occ*MATMUL(cell%h_inv, rcc)
1938 1360 : xphase = CMPLX(COS(ria), SIN(ria), dp)
1939 :
1940 : ! charge
1941 340 : trace = 0.0_dp
1942 702 : DO ispin = 1, dft_control%nspins
1943 362 : CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
1944 702 : rmom(1, 1) = rmom(1, 1) + trace
1945 : END DO
1946 :
1947 340 : zi = 0._dp
1948 340 : zij = 0._dp
1949 : zijk = 0._dp
1950 : zijkl = 0._dp
1951 :
1952 680 : DO l = 1, nmom
1953 340 : SELECT CASE (l)
1954 : CASE (1)
1955 : ! Dipole
1956 1360 : DO i = 1, 3
1957 4080 : kvec(:) = twopi*cell%h_inv(i, :)
1958 1020 : CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1959 1020 : IF (qs_env%run_rtp) THEN
1960 48 : CALL get_qs_env(qs_env, rtp=rtp)
1961 48 : CALL get_rtp(rtp, mos_new=mos_new)
1962 48 : CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1963 : ELSE
1964 972 : CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1965 : END IF
1966 1020 : zdet = CMPLX(1._dp, 0._dp, dp)
1967 2106 : DO ispin = 1, dft_control%nspins
1968 1086 : CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
1969 6252 : DO idim = 1, tmp_dim
1970 : eigrmat(ispin)%local_data(:, idim) = &
1971 : CMPLX(op_fm_set(1, ispin)%local_data(:, idim), &
1972 38586 : -op_fm_set(2, ispin)%local_data(:, idim), dp)
1973 : END DO
1974 1086 : CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
1975 1086 : zdet = zdet*zdeta
1976 3192 : IF (dft_control%nspins == 1) THEN
1977 954 : zdet = zdet*zdeta
1978 : END IF
1979 : END DO
1980 1360 : zi(i) = zdet
1981 : END DO
1982 1360 : zi = zi*xphase
1983 : CASE (2)
1984 : ! Quadrupole
1985 0 : CPABORT("Berry phase moments bigger than 1 not implemented")
1986 0 : DO i = 1, 3
1987 0 : DO j = i, 3
1988 0 : kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1989 0 : CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1990 0 : IF (qs_env%run_rtp) THEN
1991 0 : CALL get_qs_env(qs_env, rtp=rtp)
1992 0 : CALL get_rtp(rtp, mos_new=mos_new)
1993 0 : CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1994 : ELSE
1995 0 : CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1996 : END IF
1997 0 : zdet = CMPLX(1._dp, 0._dp, dp)
1998 0 : DO ispin = 1, dft_control%nspins
1999 0 : CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
2000 0 : DO idim = 1, tmp_dim
2001 : eigrmat(ispin)%local_data(:, idim) = &
2002 : CMPLX(op_fm_set(1, ispin)%local_data(:, idim), &
2003 0 : -op_fm_set(2, ispin)%local_data(:, idim), dp)
2004 : END DO
2005 0 : CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
2006 0 : zdet = zdet*zdeta
2007 0 : IF (dft_control%nspins == 1) THEN
2008 0 : zdet = zdet*zdeta
2009 : END IF
2010 : END DO
2011 0 : zij(i, j) = zdet*xphase(i)*xphase(j)
2012 0 : zij(j, i) = zdet*xphase(i)*xphase(j)
2013 : END DO
2014 : END DO
2015 : CASE (3)
2016 : ! Octapole
2017 0 : CPABORT("Berry phase moments bigger than 2 not implemented")
2018 : CASE (4)
2019 : ! Hexadecapole
2020 0 : CPABORT("Berry phase moments bigger than 3 not implemented")
2021 : CASE DEFAULT
2022 340 : CPABORT("Berry phase moments bigger than 4 not implemented")
2023 : END SELECT
2024 : END DO
2025 680 : DO l = 1, nmom
2026 340 : SELECT CASE (l)
2027 : CASE (1)
2028 : ! Dipole (apply periodic (2 Pi) boundary conditions)
2029 1360 : ci = AIMAG(LOG(zi))
2030 1360 : DO i = 1, 3
2031 1020 : IF (qq(i) + ci(i) > pi) ci(i) = ci(i) - twopi
2032 1360 : IF (qq(i) + ci(i) < -pi) ci(i) = ci(i) + twopi
2033 : END DO
2034 6460 : rmom(2:4, 1) = MATMUL(cell%hmat, ci)/twopi
2035 : CASE (2)
2036 : ! Quadrupole
2037 0 : CPABORT("Berry phase moments bigger than 1 not implemented")
2038 0 : DO i = 1, 3
2039 0 : DO j = 1, 3
2040 0 : zz = zij(i, j)/zi(i)/zi(j)
2041 0 : cij(i, j) = AIMAG(LOG(zz))/twopi
2042 : END DO
2043 : END DO
2044 0 : cij = 0.5_dp*cij/twopi/twopi
2045 0 : cij = MATMUL(MATMUL(cell%hmat, cij), TRANSPOSE(cell%hmat))
2046 0 : DO k = 4, 9
2047 0 : ix = indco(1, k + 1)
2048 0 : iy = indco(2, k + 1)
2049 0 : iz = indco(3, k + 1)
2050 0 : IF (ix == 0) THEN
2051 0 : rmom(k + 1, 1) = cij(iy, iz)
2052 0 : ELSE IF (iy == 0) THEN
2053 0 : rmom(k + 1, 1) = cij(ix, iz)
2054 0 : ELSE IF (iz == 0) THEN
2055 0 : rmom(k + 1, 1) = cij(ix, iy)
2056 : END IF
2057 : END DO
2058 : CASE (3)
2059 : ! Octapole
2060 0 : CPABORT("Berry phase moments bigger than 2 not implemented")
2061 : CASE (4)
2062 : ! Hexadecapole
2063 0 : CPABORT("Berry phase moments bigger than 3 not implemented")
2064 : CASE DEFAULT
2065 340 : CPABORT("Berry phase moments bigger than 4 not implemented")
2066 : END SELECT
2067 : END DO
2068 :
2069 1700 : rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
2070 340 : description = "[DIPOLE]"
2071 340 : CALL cp_results_erase(results=results, description=description)
2072 : CALL put_results(results=results, description=description, &
2073 340 : values=rmom(2:4, 3))
2074 340 : IF (magnetic) THEN
2075 0 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.TRUE., mmom=mmom)
2076 : ELSE
2077 340 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.TRUE.)
2078 : END IF
2079 :
2080 340 : DEALLOCATE (rmom)
2081 340 : DEALLOCATE (rlab)
2082 340 : IF (magnetic) THEN
2083 0 : DEALLOCATE (mmom)
2084 : END IF
2085 :
2086 340 : CALL dbcsr_deallocate_matrix(cosmat)
2087 340 : CALL dbcsr_deallocate_matrix(sinmat)
2088 :
2089 340 : CALL cp_fm_release(opvec)
2090 340 : CALL cp_fm_release(op_fm_set)
2091 702 : DO ispin = 1, dft_control%nspins
2092 702 : CALL cp_cfm_release(eigrmat(ispin))
2093 : END DO
2094 340 : DEALLOCATE (eigrmat)
2095 :
2096 340 : CALL timestop(handle)
2097 :
2098 1020 : END SUBROUTINE qs_moment_berry_phase
2099 :
2100 : ! **************************************************************************************************
2101 : !> \brief ...
2102 : !> \param cosmat ...
2103 : !> \param sinmat ...
2104 : !> \param mos ...
2105 : !> \param op_fm_set ...
2106 : !> \param opvec ...
2107 : ! **************************************************************************************************
2108 972 : SUBROUTINE op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
2109 :
2110 : TYPE(dbcsr_type), POINTER :: cosmat, sinmat
2111 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2112 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
2113 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: opvec
2114 :
2115 : INTEGER :: i, nao, nmo
2116 : TYPE(cp_fm_type), POINTER :: mo_coeff
2117 :
2118 1986 : DO i = 1, SIZE(op_fm_set, 2) ! spin
2119 1014 : CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
2120 1014 : CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff, opvec(i), ncol=nmo)
2121 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
2122 1014 : op_fm_set(1, i))
2123 1014 : CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff, opvec(i), ncol=nmo)
2124 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
2125 3000 : op_fm_set(2, i))
2126 : END DO
2127 :
2128 972 : END SUBROUTINE op_orbbas
2129 :
2130 : ! **************************************************************************************************
2131 : !> \brief ...
2132 : !> \param cosmat ...
2133 : !> \param sinmat ...
2134 : !> \param mos ...
2135 : !> \param op_fm_set ...
2136 : !> \param mos_new ...
2137 : ! **************************************************************************************************
2138 48 : SUBROUTINE op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
2139 :
2140 : TYPE(dbcsr_type), POINTER :: cosmat, sinmat
2141 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2142 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
2143 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
2144 :
2145 : INTEGER :: i, icol, lcol, nao, newdim, nmo
2146 : LOGICAL :: double_col, double_row
2147 : TYPE(cp_fm_struct_type), POINTER :: newstruct, newstruct1
2148 : TYPE(cp_fm_type) :: work, work1, work2
2149 : TYPE(cp_fm_type), POINTER :: mo_coeff
2150 :
2151 120 : DO i = 1, SIZE(op_fm_set, 2) ! spin
2152 72 : CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
2153 72 : CALL cp_fm_get_info(mos_new(2*i), ncol_local=lcol, ncol_global=nmo)
2154 72 : double_col = .TRUE.
2155 72 : double_row = .FALSE.
2156 : CALL cp_fm_struct_double(newstruct, &
2157 : mos_new(2*i)%matrix_struct, &
2158 : mos_new(2*i)%matrix_struct%context, &
2159 : double_col, &
2160 72 : double_row)
2161 :
2162 72 : CALL cp_fm_create(work, matrix_struct=newstruct)
2163 72 : CALL cp_fm_create(work1, matrix_struct=newstruct)
2164 72 : CALL cp_fm_create(work2, matrix_struct=newstruct)
2165 72 : CALL cp_fm_get_info(work, ncol_global=newdim)
2166 :
2167 72 : CALL cp_fm_set_all(work, 0.0_dp, 0.0_dp)
2168 336 : DO icol = 1, lcol
2169 3300 : work%local_data(:, icol) = mos_new(2*i - 1)%local_data(:, icol)
2170 3372 : work%local_data(:, icol + lcol) = mos_new(2*i)%local_data(:, icol)
2171 : END DO
2172 :
2173 72 : CALL cp_dbcsr_sm_fm_multiply(cosmat, work, work1, ncol=newdim)
2174 72 : CALL cp_dbcsr_sm_fm_multiply(sinmat, work, work2, ncol=newdim)
2175 :
2176 336 : DO icol = 1, lcol
2177 3300 : work%local_data(:, icol) = work1%local_data(:, icol) - work2%local_data(:, icol + lcol)
2178 3372 : work%local_data(:, icol + lcol) = work1%local_data(:, icol + lcol) + work2%local_data(:, icol)
2179 : END DO
2180 :
2181 72 : CALL cp_fm_release(work1)
2182 72 : CALL cp_fm_release(work2)
2183 :
2184 : CALL cp_fm_struct_double(newstruct1, &
2185 : op_fm_set(1, i)%matrix_struct, &
2186 : op_fm_set(1, i)%matrix_struct%context, &
2187 : double_col, &
2188 72 : double_row)
2189 :
2190 72 : CALL cp_fm_create(work1, matrix_struct=newstruct1)
2191 :
2192 : CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i - 1), &
2193 72 : work, 0.0_dp, work1)
2194 :
2195 336 : DO icol = 1, lcol
2196 756 : op_fm_set(1, i)%local_data(:, icol) = work1%local_data(:, icol)
2197 828 : op_fm_set(2, i)%local_data(:, icol) = work1%local_data(:, icol + lcol)
2198 : END DO
2199 :
2200 : CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i), &
2201 72 : work, 0.0_dp, work1)
2202 :
2203 336 : DO icol = 1, lcol
2204 : op_fm_set(1, i)%local_data(:, icol) = &
2205 756 : op_fm_set(1, i)%local_data(:, icol) + work1%local_data(:, icol + lcol)
2206 : op_fm_set(2, i)%local_data(:, icol) = &
2207 828 : op_fm_set(2, i)%local_data(:, icol) - work1%local_data(:, icol)
2208 : END DO
2209 :
2210 72 : CALL cp_fm_release(work)
2211 72 : CALL cp_fm_release(work1)
2212 72 : CALL cp_fm_struct_release(newstruct)
2213 264 : CALL cp_fm_struct_release(newstruct1)
2214 :
2215 : END DO
2216 :
2217 48 : END SUBROUTINE op_orbbas_rtp
2218 :
2219 : ! **************************************************************************************************
2220 : !> \brief ...
2221 : !> \param qs_env ...
2222 : !> \param magnetic ...
2223 : !> \param nmoments ...
2224 : !> \param reference ...
2225 : !> \param ref_point ...
2226 : !> \param unit_number ...
2227 : !> \param vel_reprs ...
2228 : !> \param com_nl ...
2229 : ! **************************************************************************************************
2230 778 : SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
2231 :
2232 : TYPE(qs_environment_type), POINTER :: qs_env
2233 : LOGICAL, INTENT(IN) :: magnetic
2234 : INTEGER, INTENT(IN) :: nmoments, reference
2235 : REAL(dp), DIMENSION(:), INTENT(IN), POINTER :: ref_point
2236 : INTEGER, INTENT(IN) :: unit_number
2237 : LOGICAL, INTENT(IN), OPTIONAL :: vel_reprs, com_nl
2238 :
2239 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_moment_locop'
2240 :
2241 778 : CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
2242 : CHARACTER(LEN=default_string_length) :: description
2243 : INTEGER :: akind, handle, i, ia, iatom, idir, &
2244 : ikind, ispin, ix, iy, iz, l, nm, nmom, &
2245 : order
2246 : LOGICAL :: my_com_nl, my_velreprs
2247 : REAL(dp) :: charge, dd, strace, trace
2248 778 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rrv_vrr, &
2249 778 : nlcom_rv, nlcom_rvr, nlcom_rxrv, &
2250 778 : qupole_der, rmom_vel
2251 778 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
2252 : REAL(dp), DIMENSION(3) :: rcc, ria
2253 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2254 : TYPE(cell_type), POINTER :: cell
2255 : TYPE(cp_result_type), POINTER :: results
2256 778 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom, matrix_s, moments, momentum, &
2257 778 : rho_ao
2258 778 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_der
2259 : TYPE(dbcsr_type), POINTER :: tmp_ao
2260 : TYPE(dft_control_type), POINTER :: dft_control
2261 : TYPE(distribution_1d_type), POINTER :: local_particles
2262 : TYPE(mp_para_env_type), POINTER :: para_env
2263 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2264 778 : POINTER :: sab_all, sab_orb
2265 778 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2266 778 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2267 : TYPE(qs_rho_type), POINTER :: rho
2268 :
2269 0 : CPASSERT(ASSOCIATED(qs_env))
2270 :
2271 778 : CALL timeset(routineN, handle)
2272 :
2273 778 : my_velreprs = .FALSE.
2274 778 : IF (PRESENT(vel_reprs)) my_velreprs = vel_reprs
2275 778 : IF (PRESENT(com_nl)) my_com_nl = com_nl
2276 778 : IF (my_velreprs) CALL cite_reference(Mattiat2019)
2277 :
2278 : ! reference point
2279 778 : CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
2280 :
2281 : ! only allow for moments up to maxl set by basis
2282 778 : nmom = MIN(nmoments, current_maxl)
2283 : ! electronic contribution
2284 778 : NULLIFY (dft_control, rho, cell, particle_set, qs_kind_set, results, para_env, matrix_s, rho_ao, sab_all, sab_orb)
2285 : CALL get_qs_env(qs_env, &
2286 : dft_control=dft_control, &
2287 : rho=rho, &
2288 : cell=cell, &
2289 : results=results, &
2290 : particle_set=particle_set, &
2291 : qs_kind_set=qs_kind_set, &
2292 : para_env=para_env, &
2293 : matrix_s=matrix_s, &
2294 : sab_all=sab_all, &
2295 778 : sab_orb=sab_orb)
2296 :
2297 778 : IF (my_com_nl) THEN
2298 28 : IF ((nmom >= 1) .AND. my_velreprs) THEN
2299 28 : ALLOCATE (nlcom_rv(3))
2300 112 : nlcom_rv(:) = 0._dp
2301 : END IF
2302 28 : IF ((nmom >= 2) .AND. my_velreprs) THEN
2303 28 : ALLOCATE (nlcom_rrv(6))
2304 196 : nlcom_rrv(:) = 0._dp
2305 28 : ALLOCATE (nlcom_rvr(6))
2306 196 : nlcom_rvr(:) = 0._dp
2307 28 : ALLOCATE (nlcom_rrv_vrr(6))
2308 196 : nlcom_rrv_vrr(:) = 0._dp
2309 : END IF
2310 28 : IF (magnetic) THEN
2311 6 : ALLOCATE (nlcom_rxrv(3))
2312 24 : nlcom_rxrv = 0._dp
2313 : END IF
2314 : ! Calculate non local correction terms
2315 28 : CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr, rcc)
2316 : END IF
2317 :
2318 778 : NULLIFY (moments)
2319 778 : nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
2320 778 : CALL dbcsr_allocate_matrix_set(moments, nm)
2321 3342 : DO i = 1, nm
2322 2564 : ALLOCATE (moments(i)%matrix)
2323 2564 : IF (my_velreprs .AND. (nmom >= 2)) THEN
2324 : CALL dbcsr_create(moments(i)%matrix, template=matrix_s(1)%matrix, &
2325 252 : matrix_type=dbcsr_type_symmetric)
2326 252 : CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_orb)
2327 : ELSE
2328 2312 : CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
2329 : END IF
2330 3342 : CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
2331 : END DO
2332 :
2333 : ! calculate derivatives if quadrupole in vel. reprs. is requested
2334 778 : IF (my_velreprs .AND. (nmom >= 2)) THEN
2335 28 : NULLIFY (moments_der)
2336 28 : CALL dbcsr_allocate_matrix_set(moments_der, 3, 3)
2337 112 : DO i = 1, 3 ! x, y, z
2338 364 : DO idir = 1, 3 ! d/dx, d/dy, d/dz
2339 252 : CALL dbcsr_init_p(moments_der(i, idir)%matrix)
2340 : CALL dbcsr_create(moments_der(i, idir)%matrix, template=matrix_s(1)%matrix, &
2341 252 : matrix_type=dbcsr_type_antisymmetric)
2342 252 : CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_orb)
2343 336 : CALL dbcsr_set(moments_der(i, idir)%matrix, 0.0_dp)
2344 : END DO
2345 : END DO
2346 28 : CALL build_local_moments_der_matrix(qs_env, moments_der, 1, 2, ref_point=rcc, moments=moments)
2347 : ELSE
2348 750 : CALL build_local_moment_matrix(qs_env, moments, nmom, ref_point=rcc)
2349 : END IF
2350 :
2351 778 : CALL qs_rho_get(rho, rho_ao=rho_ao)
2352 :
2353 3112 : ALLOCATE (rmom(nm + 1, 3))
2354 2334 : ALLOCATE (rlab(nm + 1))
2355 13138 : rmom = 0.0_dp
2356 4120 : rlab = ""
2357 :
2358 778 : IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2359 : ! Allocate matrix to store the matrix product to be traced (dbcsr_dot only works for products of
2360 : ! symmetric matrices)
2361 30 : NULLIFY (tmp_ao)
2362 30 : CALL dbcsr_init_p(tmp_ao)
2363 30 : CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
2364 30 : CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
2365 30 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2366 : END IF
2367 :
2368 778 : trace = 0.0_dp
2369 1626 : DO ispin = 1, dft_control%nspins
2370 848 : CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
2371 1626 : rmom(1, 1) = rmom(1, 1) + trace
2372 : END DO
2373 :
2374 3342 : DO i = 1, SIZE(moments)
2375 2564 : strace = 0._dp
2376 5338 : DO ispin = 1, dft_control%nspins
2377 2774 : IF (my_velreprs .AND. nmoments >= 2) THEN
2378 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, moments(i)%matrix, &
2379 252 : 0.0_dp, tmp_ao)
2380 252 : CALL dbcsr_trace(tmp_ao, trace)
2381 : ELSE
2382 2522 : CALL dbcsr_dot(rho_ao(ispin)%matrix, moments(i)%matrix, trace)
2383 : END IF
2384 5338 : strace = strace + trace
2385 : END DO
2386 3342 : rmom(i + 1, 1) = strace
2387 : END DO
2388 :
2389 778 : CALL dbcsr_deallocate_matrix_set(moments)
2390 :
2391 : ! nuclear contribution
2392 : CALL get_qs_env(qs_env=qs_env, &
2393 778 : local_particles=local_particles)
2394 2434 : DO ikind = 1, SIZE(local_particles%n_el)
2395 3645 : DO ia = 1, local_particles%n_el(ikind)
2396 1211 : iatom = local_particles%list(ikind)%array(ia)
2397 : ! fold atomic positions back into unit cell
2398 9688 : ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
2399 4844 : ria = ria - rcc
2400 1211 : atomic_kind => particle_set(iatom)%atomic_kind
2401 1211 : CALL get_atomic_kind(atomic_kind, kind_number=akind)
2402 1211 : CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
2403 1211 : rmom(1, 2) = rmom(1, 2) - charge
2404 6779 : DO l = 1, nm
2405 3912 : ix = indco(1, l + 1)
2406 3912 : iy = indco(2, l + 1)
2407 3912 : iz = indco(3, l + 1)
2408 3912 : dd = 1._dp
2409 3912 : IF (ix > 0) dd = dd*ria(1)**ix
2410 3912 : IF (iy > 0) dd = dd*ria(2)**iy
2411 3912 : IF (iz > 0) dd = dd*ria(3)**iz
2412 3912 : rmom(l + 1, 2) = rmom(l + 1, 2) - charge*dd
2413 5123 : CALL set_label(rlab(l + 1), ix, iy, iz)
2414 : END DO
2415 : END DO
2416 : END DO
2417 778 : CALL para_env%sum(rmom(:, 2))
2418 13138 : rmom(:, :) = -rmom(:, :)
2419 4120 : rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
2420 :
2421 : ! magnetic moments
2422 778 : IF (magnetic) THEN
2423 8 : NULLIFY (magmom)
2424 8 : CALL dbcsr_allocate_matrix_set(magmom, 3)
2425 32 : DO i = 1, SIZE(magmom)
2426 24 : CALL dbcsr_init_p(magmom(i)%matrix)
2427 : CALL dbcsr_create(magmom(i)%matrix, template=matrix_s(1)%matrix, &
2428 24 : matrix_type=dbcsr_type_antisymmetric)
2429 24 : CALL cp_dbcsr_alloc_block_from_nbl(magmom(i)%matrix, sab_orb)
2430 32 : CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
2431 : END DO
2432 :
2433 8 : CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
2434 :
2435 24 : ALLOCATE (mmom(SIZE(magmom)))
2436 32 : mmom(:) = 0.0_dp
2437 8 : IF (qs_env%run_rtp) THEN
2438 : ! get imaginary part of the density in rho_ao (the real part is not needed since the trace of the product
2439 : ! of a symmetric (REAL(rho_ao)) and an anti-symmetric (L_AO) matrix is zero)
2440 : ! There may be other cases, where the imaginary part of the density is relevant
2441 4 : NULLIFY (rho_ao)
2442 4 : CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2443 : END IF
2444 : ! if the density is purely real this is an expensive way to calculate zero
2445 32 : DO i = 1, SIZE(magmom)
2446 24 : strace = 0._dp
2447 48 : DO ispin = 1, dft_control%nspins
2448 24 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2449 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, magmom(i)%matrix, &
2450 24 : 0.0_dp, tmp_ao)
2451 24 : CALL dbcsr_trace(tmp_ao, trace)
2452 48 : strace = strace + trace
2453 : END DO
2454 32 : mmom(i) = strace
2455 : END DO
2456 :
2457 8 : CALL dbcsr_deallocate_matrix_set(magmom)
2458 : END IF
2459 :
2460 : ! velocity representations
2461 778 : IF (my_velreprs) THEN
2462 84 : ALLOCATE (rmom_vel(nm))
2463 280 : rmom_vel = 0.0_dp
2464 :
2465 84 : DO order = 1, nmom
2466 28 : SELECT CASE (order)
2467 :
2468 : CASE (1) ! expectation value of momentum
2469 28 : NULLIFY (momentum)
2470 28 : CALL dbcsr_allocate_matrix_set(momentum, 3)
2471 112 : DO i = 1, 3
2472 84 : CALL dbcsr_init_p(momentum(i)%matrix)
2473 : CALL dbcsr_create(momentum(i)%matrix, template=matrix_s(1)%matrix, &
2474 84 : matrix_type=dbcsr_type_antisymmetric)
2475 84 : CALL cp_dbcsr_alloc_block_from_nbl(momentum(i)%matrix, sab_orb)
2476 112 : CALL dbcsr_set(momentum(i)%matrix, 0.0_dp)
2477 : END DO
2478 28 : CALL build_lin_mom_matrix(qs_env, momentum)
2479 :
2480 : ! imaginary part of the density for RTP, real part gives 0 since momentum is antisymmetric
2481 28 : IF (qs_env%run_rtp) THEN
2482 22 : NULLIFY (rho_ao)
2483 22 : CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2484 88 : DO idir = 1, SIZE(momentum)
2485 66 : strace = 0._dp
2486 132 : DO ispin = 1, dft_control%nspins
2487 66 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2488 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, momentum(idir)%matrix, &
2489 66 : 0.0_dp, tmp_ao)
2490 66 : CALL dbcsr_trace(tmp_ao, trace)
2491 132 : strace = strace + trace
2492 : END DO
2493 88 : rmom_vel(idir) = rmom_vel(idir) + strace
2494 : END DO
2495 : END IF
2496 :
2497 28 : CALL dbcsr_deallocate_matrix_set(momentum)
2498 :
2499 : CASE (2) ! expectation value of quadrupole moment in vel. reprs.
2500 28 : ALLOCATE (qupole_der(9)) ! will contain the expectation values of r_\alpha * d/d r_\beta
2501 280 : qupole_der = 0._dp
2502 :
2503 28 : NULLIFY (rho_ao)
2504 28 : CALL qs_rho_get(rho, rho_ao=rho_ao)
2505 :
2506 : ! Calculate expectation value over real part
2507 28 : trace = 0._dp
2508 112 : DO i = 1, 3
2509 364 : DO idir = 1, 3
2510 252 : strace = 0._dp
2511 504 : DO ispin = 1, dft_control%nspins
2512 252 : CALL dbcsr_set(tmp_ao, 0._dp)
2513 252 : CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2514 252 : CALL dbcsr_trace(tmp_ao, trace)
2515 504 : strace = strace + trace
2516 : END DO
2517 336 : qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2518 : END DO
2519 : END DO
2520 :
2521 28 : IF (qs_env%run_rtp) THEN
2522 22 : NULLIFY (rho_ao)
2523 22 : CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2524 :
2525 : ! Calculate expectation value over imaginary part
2526 22 : trace = 0._dp
2527 88 : DO i = 1, 3
2528 286 : DO idir = 1, 3
2529 198 : strace = 0._dp
2530 396 : DO ispin = 1, dft_control%nspins
2531 198 : CALL dbcsr_set(tmp_ao, 0._dp)
2532 198 : CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2533 198 : CALL dbcsr_trace(tmp_ao, trace)
2534 396 : strace = strace + trace
2535 : END DO
2536 264 : qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2537 : END DO
2538 : END DO
2539 : END IF
2540 :
2541 : ! calculate vel. reprs. of quadrupole moment from derivatives
2542 28 : rmom_vel(4) = -2*qupole_der(1) - rmom(1, 1)
2543 28 : rmom_vel(5) = -qupole_der(2) - qupole_der(4)
2544 28 : rmom_vel(6) = -qupole_der(3) - qupole_der(7)
2545 28 : rmom_vel(7) = -2*qupole_der(5) - rmom(1, 1)
2546 28 : rmom_vel(8) = -qupole_der(6) - qupole_der(8)
2547 28 : rmom_vel(9) = -2*qupole_der(9) - rmom(1, 1)
2548 :
2549 84 : DEALLOCATE (qupole_der)
2550 : CASE DEFAULT
2551 : END SELECT
2552 : END DO
2553 : END IF
2554 :
2555 778 : IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2556 30 : CALL dbcsr_deallocate_matrix(tmp_ao)
2557 : END IF
2558 778 : IF (my_velreprs .AND. (nmoments >= 2)) THEN
2559 28 : CALL dbcsr_deallocate_matrix_set(moments_der)
2560 : END IF
2561 :
2562 778 : description = "[DIPOLE]"
2563 778 : CALL cp_results_erase(results=results, description=description)
2564 : CALL put_results(results=results, description=description, &
2565 778 : values=rmom(2:4, 3))
2566 :
2567 778 : IF (magnetic .AND. my_velreprs) THEN
2568 6 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., mmom=mmom, rmom_vel=rmom_vel)
2569 772 : ELSE IF (magnetic .AND. .NOT. my_velreprs) THEN
2570 2 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., mmom=mmom)
2571 770 : ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2572 22 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., rmom_vel=rmom_vel)
2573 : ELSE
2574 748 : CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE.)
2575 : END IF
2576 :
2577 778 : IF (my_com_nl) THEN
2578 28 : IF (magnetic) THEN
2579 24 : mmom(:) = nlcom_rxrv(:)
2580 : END IF
2581 28 : IF (my_velreprs .AND. (nmom >= 1)) THEN
2582 28 : DEALLOCATE (rmom_vel)
2583 28 : ALLOCATE (rmom_vel(21))
2584 112 : rmom_vel(1:3) = nlcom_rv
2585 : END IF
2586 28 : IF (my_velreprs .AND. (nmom >= 2)) THEN
2587 196 : rmom_vel(4:9) = nlcom_rrv
2588 196 : rmom_vel(10:15) = nlcom_rvr
2589 196 : rmom_vel(16:21) = nlcom_rrv_vrr
2590 : END IF
2591 28 : IF (magnetic .AND. .NOT. my_velreprs) THEN
2592 0 : CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom)
2593 28 : ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2594 22 : CALL print_moments_nl(unit_number, nmom, rlab, rmom_vel=rmom_vel)
2595 6 : ELSE IF (my_velreprs .AND. magnetic) THEN
2596 6 : CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom, rmom_vel=rmom_vel)
2597 : END IF
2598 :
2599 : END IF
2600 :
2601 778 : IF (my_com_nl) THEN
2602 28 : IF (nmom >= 1 .AND. my_velreprs) DEALLOCATE (nlcom_rv)
2603 28 : IF (nmom >= 2 .AND. my_velreprs) THEN
2604 28 : DEALLOCATE (nlcom_rrv)
2605 28 : DEALLOCATE (nlcom_rvr)
2606 28 : DEALLOCATE (nlcom_rrv_vrr)
2607 : END IF
2608 28 : IF (magnetic) DEALLOCATE (nlcom_rxrv)
2609 : END IF
2610 :
2611 778 : DEALLOCATE (rmom)
2612 778 : DEALLOCATE (rlab)
2613 778 : IF (magnetic) THEN
2614 8 : DEALLOCATE (mmom)
2615 : END IF
2616 778 : IF (my_velreprs) THEN
2617 28 : DEALLOCATE (rmom_vel)
2618 : END IF
2619 :
2620 778 : CALL timestop(handle)
2621 :
2622 1556 : END SUBROUTINE qs_moment_locop
2623 :
2624 : ! **************************************************************************************************
2625 : !> \brief ...
2626 : !> \param label ...
2627 : !> \param ix ...
2628 : !> \param iy ...
2629 : !> \param iz ...
2630 : ! **************************************************************************************************
2631 4932 : SUBROUTINE set_label(label, ix, iy, iz)
2632 : CHARACTER(LEN=*), INTENT(OUT) :: label
2633 : INTEGER, INTENT(IN) :: ix, iy, iz
2634 :
2635 : INTEGER :: i
2636 :
2637 4932 : label = ""
2638 6709 : DO i = 1, ix
2639 6709 : WRITE (label(i:), "(A1)") "X"
2640 : END DO
2641 6709 : DO i = ix + 1, ix + iy
2642 6709 : WRITE (label(i:), "(A1)") "Y"
2643 : END DO
2644 6709 : DO i = ix + iy + 1, ix + iy + iz
2645 6709 : WRITE (label(i:), "(A1)") "Z"
2646 : END DO
2647 :
2648 4932 : END SUBROUTINE set_label
2649 :
2650 : ! **************************************************************************************************
2651 : !> \brief ...
2652 : !> \param unit_number ...
2653 : !> \param nmom ...
2654 : !> \param rmom ...
2655 : !> \param rlab ...
2656 : !> \param rcc ...
2657 : !> \param cell ...
2658 : !> \param periodic ...
2659 : !> \param mmom ...
2660 : !> \param rmom_vel ...
2661 : ! **************************************************************************************************
2662 1118 : SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmom, rmom_vel)
2663 : INTEGER, INTENT(IN) :: unit_number, nmom
2664 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: rmom
2665 : CHARACTER(LEN=8), DIMENSION(:) :: rlab
2666 : REAL(dp), DIMENSION(3), INTENT(IN) :: rcc
2667 : TYPE(cell_type), POINTER :: cell
2668 : LOGICAL :: periodic
2669 : REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2670 :
2671 : INTEGER :: i, i0, i1, j, l
2672 : REAL(dp) :: dd
2673 :
2674 1118 : IF (unit_number > 0) THEN
2675 1724 : DO l = 0, nmom
2676 569 : SELECT CASE (l)
2677 : CASE (0)
2678 569 : WRITE (unit_number, "(T3,A,T33,3F16.8)") "Reference Point [Bohr]", rcc
2679 569 : WRITE (unit_number, "(T3,A)") "Charges"
2680 : WRITE (unit_number, "(T5,A,T18,F14.8,T36,A,T42,F14.8,T60,A,T67,F14.8)") &
2681 569 : "Electronic=", rmom(1, 1), "Core=", rmom(1, 2), "Total=", rmom(1, 3)
2682 : CASE (1)
2683 569 : IF (periodic) THEN
2684 180 : WRITE (unit_number, "(T3,A)") "Dipole vectors are based on the periodic (Berry phase) operator."
2685 180 : WRITE (unit_number, "(T3,A)") "They are defined modulo integer multiples of the cell matrix [Debye]."
2686 720 : WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[X] [", cell%hmat(1, :)*debye, "] [i]"
2687 720 : WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Y]=[", cell%hmat(2, :)*debye, "]*[j]"
2688 720 : WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Z] [", cell%hmat(3, :)*debye, "] [k]"
2689 : ELSE
2690 389 : WRITE (unit_number, "(T3,A)") "Dipoles are based on the traditional operator."
2691 : END IF
2692 2276 : dd = SQRT(SUM(rmom(2:4, 3)**2))*debye
2693 569 : WRITE (unit_number, "(T3,A)") "Dipole moment [Debye]"
2694 : WRITE (unit_number, "(T5,3(A,A,E15.7,1X),T60,A,T68,F13.7)") &
2695 2276 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye, i=2, 4), "Total=", dd
2696 : CASE (2)
2697 15 : WRITE (unit_number, "(T3,A)") "Quadrupole moment [Debye*Angstrom]"
2698 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2699 60 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=5, 7)
2700 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2701 60 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=8, 10)
2702 : CASE (3)
2703 1 : WRITE (unit_number, "(T3,A)") "Octapole moment [Debye*Angstrom**2]"
2704 : WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2705 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=11, 14)
2706 : WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2707 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=15, 18)
2708 : WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2709 3 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=19, 20)
2710 : CASE (4)
2711 1 : WRITE (unit_number, "(T3,A)") "Hexadecapole moment [Debye*Angstrom**3]"
2712 : WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2713 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=21, 24)
2714 : WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2715 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=25, 28)
2716 : WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2717 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=29, 32)
2718 : WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2719 5 : (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=32, 35)
2720 : CASE DEFAULT
2721 0 : WRITE (unit_number, "(T3,A,A,I2)") "Higher moment [Debye*Angstrom**(L-1)]", &
2722 0 : " L=", l
2723 0 : i0 = (6 + 11*(l - 1) + 6*(l - 1)**2 + (l - 1)**3)/6
2724 0 : i1 = (6 + 11*l + 6*l**2 + l**3)/6 - 1
2725 0 : dd = debye/(bohr)**(l - 1)
2726 1155 : DO i = i0, i1, 3
2727 : WRITE (unit_number, "(T18,3(A,A,F14.8,4X))") &
2728 0 : (TRIM(rlab(j + 1)), "=", rmom(j + 1, 3)*dd, j=i, MIN(i1, i + 2))
2729 : END DO
2730 : END SELECT
2731 : END DO
2732 569 : IF (PRESENT(mmom)) THEN
2733 4 : IF (nmom >= 1) THEN
2734 16 : dd = SQRT(SUM(mmom(1:3)**2))
2735 4 : WRITE (unit_number, "(T3,A)") "Orbital angular momentum [a. u.]"
2736 : WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2737 16 : (TRIM(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2738 : END IF
2739 : END IF
2740 569 : IF (PRESENT(rmom_vel)) THEN
2741 42 : DO l = 1, nmom
2742 14 : SELECT CASE (l)
2743 : CASE (1)
2744 56 : dd = SQRT(SUM(rmom_vel(1:3)**2))
2745 14 : WRITE (unit_number, "(T3,A)") "Expectation value of momentum operator [a. u.]"
2746 : WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2747 56 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2748 : CASE (2)
2749 14 : WRITE (unit_number, "(T3,A)") "Expectation value of quadrupole operator in vel. repr. [a. u.]"
2750 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2751 56 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2752 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2753 84 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2754 : CASE DEFAULT
2755 : END SELECT
2756 : END DO
2757 : END IF
2758 : END IF
2759 :
2760 1118 : END SUBROUTINE print_moments
2761 :
2762 : ! **************************************************************************************************
2763 : !> \brief ...
2764 : !> \param unit_number ...
2765 : !> \param nmom ...
2766 : !> \param rlab ...
2767 : !> \param mmom ...
2768 : !> \param rmom_vel ...
2769 : ! **************************************************************************************************
2770 28 : SUBROUTINE print_moments_nl(unit_number, nmom, rlab, mmom, rmom_vel)
2771 : INTEGER, INTENT(IN) :: unit_number, nmom
2772 : CHARACTER(LEN=8), DIMENSION(:) :: rlab
2773 : REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2774 :
2775 : INTEGER :: i, l
2776 : REAL(dp) :: dd
2777 :
2778 28 : IF (unit_number > 0) THEN
2779 14 : IF (PRESENT(mmom)) THEN
2780 3 : IF (nmom >= 1) THEN
2781 12 : dd = SQRT(SUM(mmom(1:3)**2))
2782 3 : WRITE (unit_number, "(T3,A)") "Expectation value of rx[r,V_nl] [a. u.]"
2783 : WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2784 12 : (TRIM(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2785 : END IF
2786 : END IF
2787 14 : IF (PRESENT(rmom_vel)) THEN
2788 42 : DO l = 1, nmom
2789 14 : SELECT CASE (l)
2790 : CASE (1)
2791 56 : dd = SQRT(SUM(rmom_vel(1:3)**2))
2792 14 : WRITE (unit_number, "(T3,A)") "Expectation value of [r,V_nl] [a. u.]"
2793 : WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2794 56 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2795 : CASE (2)
2796 14 : WRITE (unit_number, "(T3,A)") "Expectation value of [rr,V_nl] [a. u.]"
2797 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2798 56 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2799 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2800 56 : (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2801 14 : WRITE (unit_number, "(T3,A)") "Expectation value of r x V_nl x r [a. u.]"
2802 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2803 56 : (TRIM(rlab(i + 1 - 6)), "=", rmom_vel(i), i=10, 12)
2804 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2805 56 : (TRIM(rlab(i + 1 - 6)), "=", rmom_vel(i), i=13, 15)
2806 14 : WRITE (unit_number, "(T3,A)") "Expectation value of r x r x V_nl + V_nl x r x r [a. u.]"
2807 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2808 56 : (TRIM(rlab(i + 1 - 12)), "=", rmom_vel(i), i=16, 18)
2809 : WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2810 84 : (TRIM(rlab(i + 1 - 12)), "=", rmom_vel(i), i=19, 21)
2811 : CASE DEFAULT
2812 : END SELECT
2813 : END DO
2814 : END IF
2815 : END IF
2816 :
2817 28 : END SUBROUTINE print_moments_nl
2818 :
2819 : ! **************************************************************************************************
2820 : !> \brief Calculate the expectation value of operators related to non-local potential:
2821 : !> [r, Vnl], noted rv
2822 : !> r x [r,Vnl], noted rxrv
2823 : !> [rr,Vnl], noted rrv
2824 : !> r x Vnl x r, noted rvr
2825 : !> r x r x Vnl + Vnl x r x r, noted rrv_vrr
2826 : !> Note that the 3 first operator are commutator while the 2 last
2827 : !> are not. For reading clarity the same notation is used for all 5
2828 : !> operators.
2829 : !> \param qs_env ...
2830 : !> \param nlcom_rv ...
2831 : !> \param nlcom_rxrv ...
2832 : !> \param nlcom_rrv ...
2833 : !> \param nlcom_rvr ...
2834 : !> \param nlcom_rrv_vrr ...
2835 : !> \param ref_point ...
2836 : ! **************************************************************************************************
2837 28 : SUBROUTINE calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, &
2838 : nlcom_rrv_vrr, ref_point)
2839 :
2840 : TYPE(qs_environment_type), POINTER :: qs_env
2841 : REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: nlcom_rv, nlcom_rxrv, nlcom_rrv, &
2842 : nlcom_rvr, nlcom_rrv_vrr
2843 : REAL(dp), DIMENSION(3) :: ref_point
2844 :
2845 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_commutator_nl_terms'
2846 :
2847 : INTEGER :: handle, ind, ispin
2848 : LOGICAL :: calc_rrv, calc_rrv_vrr, calc_rv, &
2849 : calc_rvr, calc_rxrv
2850 : REAL(dp) :: eps_ppnl, strace, trace
2851 : TYPE(cell_type), POINTER :: cell
2852 28 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rrv, matrix_rrv_vrr, matrix_rv, &
2853 28 : matrix_rvr, matrix_rxrv, matrix_s, &
2854 28 : rho_ao
2855 : TYPE(dbcsr_type), POINTER :: tmp_ao
2856 : TYPE(dft_control_type), POINTER :: dft_control
2857 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2858 28 : POINTER :: sab_all, sab_orb, sap_ppnl
2859 28 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2860 28 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2861 : TYPE(qs_rho_type), POINTER :: rho
2862 :
2863 28 : CALL timeset(routineN, handle)
2864 :
2865 28 : calc_rv = .FALSE.
2866 28 : calc_rxrv = .FALSE.
2867 28 : calc_rrv = .FALSE.
2868 28 : calc_rvr = .FALSE.
2869 28 : calc_rrv_vrr = .FALSE.
2870 :
2871 : ! rv, rxrv and rrv are commutator matrices: anti-symmetric.
2872 : ! The real part of the density matrix rho_ao is symmetric so that
2873 : ! the expectation value of real density matrix is zero. Hence, if
2874 : ! the density matrix is real, no need to compute these quantities.
2875 : ! This is not the case for rvr and rrv_vrr which are symmetric.
2876 :
2877 28 : IF (ALLOCATED(nlcom_rv)) THEN
2878 112 : nlcom_rv(:) = 0._dp
2879 28 : IF (qs_env%run_rtp) calc_rv = .TRUE.
2880 : END IF
2881 28 : IF (ALLOCATED(nlcom_rxrv)) THEN
2882 24 : nlcom_rxrv(:) = 0._dp
2883 6 : IF (qs_env%run_rtp) calc_rxrv = .TRUE.
2884 : END IF
2885 28 : IF (ALLOCATED(nlcom_rrv)) THEN
2886 196 : nlcom_rrv(:) = 0._dp
2887 28 : IF (qs_env%run_rtp) calc_rrv = .TRUE.
2888 : END IF
2889 28 : IF (ALLOCATED(nlcom_rvr)) THEN
2890 196 : nlcom_rvr(:) = 0._dp
2891 : calc_rvr = .TRUE.
2892 : END IF
2893 28 : IF (ALLOCATED(nlcom_rrv_vrr)) THEN
2894 196 : nlcom_rrv_vrr(:) = 0._dp
2895 : calc_rrv_vrr = .TRUE.
2896 : END IF
2897 :
2898 28 : IF (.NOT. (calc_rv .OR. calc_rrv .OR. calc_rxrv &
2899 0 : .OR. calc_rvr .OR. calc_rrv_vrr)) RETURN
2900 :
2901 28 : NULLIFY (cell, matrix_s, particle_set, qs_kind_set, rho, sab_all, sab_orb, sap_ppnl)
2902 : CALL get_qs_env(qs_env, &
2903 : cell=cell, &
2904 : dft_control=dft_control, &
2905 : matrix_s=matrix_s, &
2906 : particle_set=particle_set, &
2907 : qs_kind_set=qs_kind_set, &
2908 : rho=rho, &
2909 : sab_orb=sab_orb, &
2910 : sab_all=sab_all, &
2911 28 : sap_ppnl=sap_ppnl)
2912 :
2913 28 : eps_ppnl = dft_control%qs_control%eps_ppnl
2914 :
2915 : ! Allocate storage
2916 28 : NULLIFY (matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr)
2917 28 : IF (calc_rv) THEN
2918 22 : CALL dbcsr_allocate_matrix_set(matrix_rv, 3)
2919 88 : DO ind = 1, 3
2920 66 : CALL dbcsr_init_p(matrix_rv(ind)%matrix)
2921 : CALL dbcsr_create(matrix_rv(ind)%matrix, template=matrix_s(1)%matrix, &
2922 66 : matrix_type=dbcsr_type_antisymmetric)
2923 66 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(ind)%matrix, sab_orb)
2924 88 : CALL dbcsr_set(matrix_rv(ind)%matrix, 0._dp)
2925 : END DO
2926 : END IF
2927 :
2928 28 : IF (calc_rxrv) THEN
2929 4 : CALL dbcsr_allocate_matrix_set(matrix_rxrv, 3)
2930 16 : DO ind = 1, 3
2931 12 : CALL dbcsr_init_p(matrix_rxrv(ind)%matrix)
2932 : CALL dbcsr_create(matrix_rxrv(ind)%matrix, template=matrix_s(1)%matrix, &
2933 12 : matrix_type=dbcsr_type_antisymmetric)
2934 12 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_rxrv(ind)%matrix, sab_orb)
2935 16 : CALL dbcsr_set(matrix_rxrv(ind)%matrix, 0._dp)
2936 : END DO
2937 : END IF
2938 :
2939 28 : IF (calc_rrv) THEN
2940 22 : CALL dbcsr_allocate_matrix_set(matrix_rrv, 6)
2941 154 : DO ind = 1, 6
2942 132 : CALL dbcsr_init_p(matrix_rrv(ind)%matrix)
2943 : CALL dbcsr_create(matrix_rrv(ind)%matrix, template=matrix_s(1)%matrix, &
2944 132 : matrix_type=dbcsr_type_antisymmetric)
2945 132 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv(ind)%matrix, sab_orb)
2946 154 : CALL dbcsr_set(matrix_rrv(ind)%matrix, 0._dp)
2947 : END DO
2948 : END IF
2949 :
2950 28 : IF (calc_rvr) THEN
2951 28 : CALL dbcsr_allocate_matrix_set(matrix_rvr, 6)
2952 196 : DO ind = 1, 6
2953 168 : CALL dbcsr_init_p(matrix_rvr(ind)%matrix)
2954 : CALL dbcsr_create(matrix_rvr(ind)%matrix, template=matrix_s(1)%matrix, &
2955 168 : matrix_type=dbcsr_type_symmetric)
2956 168 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_rvr(ind)%matrix, sab_orb)
2957 196 : CALL dbcsr_set(matrix_rvr(ind)%matrix, 0._dp)
2958 : END DO
2959 : END IF
2960 28 : IF (calc_rrv_vrr) THEN
2961 28 : CALL dbcsr_allocate_matrix_set(matrix_rrv_vrr, 6)
2962 196 : DO ind = 1, 6
2963 168 : CALL dbcsr_init_p(matrix_rrv_vrr(ind)%matrix)
2964 : CALL dbcsr_create(matrix_rrv_vrr(ind)%matrix, template=matrix_s(1)%matrix, &
2965 168 : matrix_type=dbcsr_type_symmetric)
2966 168 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv_vrr(ind)%matrix, sab_orb)
2967 196 : CALL dbcsr_set(matrix_rrv_vrr(ind)%matrix, 0._dp)
2968 : END DO
2969 : END IF
2970 :
2971 : ! calculate evaluation of operators in AO basis set
2972 : CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
2973 : matrix_rxrv=matrix_rxrv, matrix_rrv=matrix_rrv, matrix_rvr=matrix_rvr, &
2974 28 : matrix_rrv_vrr=matrix_rrv_vrr, ref_point=ref_point)
2975 :
2976 : ! Calculate expectation values
2977 : ! Real part
2978 28 : NULLIFY (tmp_ao)
2979 28 : CALL dbcsr_init_p(tmp_ao)
2980 28 : CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
2981 28 : CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
2982 28 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2983 :
2984 28 : IF (calc_rvr .OR. calc_rrv_vrr) THEN
2985 28 : NULLIFY (rho_ao)
2986 28 : CALL qs_rho_get(rho, rho_ao=rho_ao)
2987 :
2988 28 : IF (calc_rvr) THEN
2989 28 : trace = 0._dp
2990 196 : DO ind = 1, SIZE(matrix_rvr)
2991 168 : strace = 0._dp
2992 336 : DO ispin = 1, dft_control%nspins
2993 168 : CALL dbcsr_set(tmp_ao, 0.0_dp)
2994 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rvr(ind)%matrix, &
2995 168 : 0.0_dp, tmp_ao)
2996 168 : CALL dbcsr_trace(tmp_ao, trace)
2997 336 : strace = strace + trace
2998 : END DO
2999 196 : nlcom_rvr(ind) = nlcom_rvr(ind) + strace
3000 : END DO
3001 : END IF
3002 :
3003 28 : IF (calc_rrv_vrr) THEN
3004 28 : trace = 0._dp
3005 196 : DO ind = 1, SIZE(matrix_rrv_vrr)
3006 168 : strace = 0._dp
3007 336 : DO ispin = 1, dft_control%nspins
3008 168 : CALL dbcsr_set(tmp_ao, 0.0_dp)
3009 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv_vrr(ind)%matrix, &
3010 168 : 0.0_dp, tmp_ao)
3011 168 : CALL dbcsr_trace(tmp_ao, trace)
3012 336 : strace = strace + trace
3013 : END DO
3014 196 : nlcom_rrv_vrr(ind) = nlcom_rrv_vrr(ind) + strace
3015 : END DO
3016 : END IF
3017 : END IF
3018 :
3019 : ! imagninary part of the density matrix
3020 28 : NULLIFY (rho_ao)
3021 28 : CALL qs_rho_get(rho, rho_ao_im=rho_ao)
3022 :
3023 28 : IF (calc_rv) THEN
3024 22 : trace = 0._dp
3025 88 : DO ind = 1, SIZE(matrix_rv)
3026 66 : strace = 0._dp
3027 132 : DO ispin = 1, dft_control%nspins
3028 66 : CALL dbcsr_set(tmp_ao, 0.0_dp)
3029 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rv(ind)%matrix, &
3030 66 : 0.0_dp, tmp_ao)
3031 66 : CALL dbcsr_trace(tmp_ao, trace)
3032 132 : strace = strace + trace
3033 : END DO
3034 88 : nlcom_rv(ind) = nlcom_rv(ind) + strace
3035 : END DO
3036 : END IF
3037 :
3038 28 : IF (calc_rrv) THEN
3039 22 : trace = 0._dp
3040 154 : DO ind = 1, SIZE(matrix_rrv)
3041 132 : strace = 0._dp
3042 264 : DO ispin = 1, dft_control%nspins
3043 132 : CALL dbcsr_set(tmp_ao, 0.0_dp)
3044 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv(ind)%matrix, &
3045 132 : 0.0_dp, tmp_ao)
3046 132 : CALL dbcsr_trace(tmp_ao, trace)
3047 264 : strace = strace + trace
3048 : END DO
3049 154 : nlcom_rrv(ind) = nlcom_rrv(ind) + strace
3050 : END DO
3051 : END IF
3052 :
3053 28 : IF (calc_rxrv) THEN
3054 4 : trace = 0._dp
3055 16 : DO ind = 1, SIZE(matrix_rxrv)
3056 12 : strace = 0._dp
3057 24 : DO ispin = 1, dft_control%nspins
3058 12 : CALL dbcsr_set(tmp_ao, 0.0_dp)
3059 : CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rxrv(ind)%matrix, &
3060 12 : 0.0_dp, tmp_ao)
3061 12 : CALL dbcsr_trace(tmp_ao, trace)
3062 24 : strace = strace + trace
3063 : END DO
3064 16 : nlcom_rxrv(ind) = nlcom_rxrv(ind) + strace
3065 : END DO
3066 : END IF
3067 28 : CALL dbcsr_deallocate_matrix(tmp_ao)
3068 28 : IF (calc_rv) CALL dbcsr_deallocate_matrix_set(matrix_rv)
3069 28 : IF (calc_rxrv) CALL dbcsr_deallocate_matrix_set(matrix_rxrv)
3070 28 : IF (calc_rrv) CALL dbcsr_deallocate_matrix_set(matrix_rrv)
3071 28 : IF (calc_rvr) CALL dbcsr_deallocate_matrix_set(matrix_rvr)
3072 28 : IF (calc_rrv_vrr) CALL dbcsr_deallocate_matrix_set(matrix_rrv_vrr)
3073 :
3074 28 : CALL timestop(handle)
3075 28 : END SUBROUTINE calculate_commutator_nl_terms
3076 :
3077 : ! *****************************************************************************
3078 : !> \brief ...
3079 : !> \param qs_env ...
3080 : !> \param difdip ...
3081 : !> \param deltaR ...
3082 : !> \param order ...
3083 : !> \param rcc ...
3084 : !> \note calculate matrix elements <a|r_beta |db/dR_alpha > + <da/dR_alpha | r_beta | b >
3085 : !> be aware: < a | r_beta| db/dR_alpha > = - < da/dR_alpha | r_beta | b > only valid
3086 : !> if alpha .neq.beta
3087 : !> if alpha=beta: < a | r_beta| db/dR_alpha > = - < da/dR_alpha | r_beta | b > - < a | b >
3088 : !> modified from qs_efield_mo_derivatives
3089 : !> SL July 2015
3090 : ! **************************************************************************************************
3091 468 : SUBROUTINE dipole_deriv_ao(qs_env, difdip, deltaR, order, rcc)
3092 : TYPE(qs_environment_type), POINTER :: qs_env
3093 : TYPE(dbcsr_p_type), DIMENSION(:, :), &
3094 : INTENT(INOUT), POINTER :: difdip
3095 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
3096 : POINTER :: deltaR
3097 : INTEGER, INTENT(IN) :: order
3098 : REAL(KIND=dp), DIMENSION(3), OPTIONAL :: rcc
3099 :
3100 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dipole_deriv_ao'
3101 :
3102 : INTEGER :: handle, i, iatom, icol, idir, ikind, inode, irow, iset, j, jatom, jkind, jset, &
3103 : last_jatom, lda, ldab, ldb, M_dim, maxsgf, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
3104 : sgfb
3105 468 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
3106 468 : npgfb, nsgfa, nsgfb
3107 468 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
3108 : LOGICAL :: found
3109 : REAL(dp) :: dab
3110 : REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
3111 468 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
3112 468 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
3113 468 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
3114 468 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
3115 468 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab
3116 468 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: difmab2
3117 468 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mint, mint2
3118 : TYPE(cell_type), POINTER :: cell
3119 468 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
3120 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
3121 : TYPE(neighbor_list_iterator_p_type), &
3122 468 : DIMENSION(:), POINTER :: nl_iterator
3123 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3124 468 : POINTER :: sab_all, sab_orb
3125 468 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3126 468 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3127 : TYPE(qs_kind_type), POINTER :: qs_kind
3128 :
3129 468 : CALL timeset(routineN, handle)
3130 :
3131 468 : NULLIFY (cell, particle_set, qs_kind_set, sab_orb, sab_all)
3132 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, &
3133 468 : qs_kind_set=qs_kind_set, sab_orb=sab_orb, sab_all=sab_all)
3134 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
3135 468 : maxco=ldab, maxsgf=maxsgf)
3136 :
3137 468 : nkind = SIZE(qs_kind_set)
3138 468 : natom = SIZE(particle_set)
3139 :
3140 468 : M_dim = ncoset(order) - 1
3141 :
3142 468 : IF (PRESENT(rcc)) THEN
3143 468 : rc = rcc
3144 : ELSE
3145 0 : rc = 0._dp
3146 : END IF
3147 :
3148 2340 : ALLOCATE (basis_set_list(nkind))
3149 :
3150 2340 : ALLOCATE (mab(ldab, ldab, M_dim))
3151 2808 : ALLOCATE (difmab2(ldab, ldab, M_dim, 3))
3152 1872 : ALLOCATE (work(ldab, maxsgf))
3153 6084 : ALLOCATE (mint(3, 3))
3154 6084 : ALLOCATE (mint2(3, 3))
3155 :
3156 383760 : mab(1:ldab, 1:ldab, 1:M_dim) = 0.0_dp
3157 1151748 : difmab2(1:ldab, 1:ldab, 1:M_dim, 1:3) = 0.0_dp
3158 32292 : work(1:ldab, 1:maxsgf) = 0.0_dp
3159 :
3160 1872 : DO i = 1, 3
3161 6084 : DO j = 1, 3
3162 4212 : NULLIFY (mint(i, j)%block)
3163 5616 : NULLIFY (mint2(i, j)%block)
3164 : END DO
3165 : END DO
3166 :
3167 : ! Set the basis_set_list(nkind) to point to the corresponding basis sets
3168 1404 : DO ikind = 1, nkind
3169 936 : qs_kind => qs_kind_set(ikind)
3170 936 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
3171 1404 : IF (ASSOCIATED(basis_set_a)) THEN
3172 936 : basis_set_list(ikind)%gto_basis_set => basis_set_a
3173 : ELSE
3174 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
3175 : END IF
3176 : END DO
3177 :
3178 468 : CALL neighbor_list_iterator_create(nl_iterator, sab_all)
3179 20586 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3180 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
3181 20118 : iatom=iatom, jatom=jatom, r=rab)
3182 :
3183 20118 : basis_set_a => basis_set_list(ikind)%gto_basis_set
3184 20118 : basis_set_b => basis_set_list(jkind)%gto_basis_set
3185 20118 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
3186 20118 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
3187 :
3188 : ! basis ikind
3189 20118 : first_sgfa => basis_set_a%first_sgf
3190 20118 : la_max => basis_set_a%lmax
3191 20118 : la_min => basis_set_a%lmin
3192 20118 : npgfa => basis_set_a%npgf
3193 20118 : nseta = basis_set_a%nset
3194 20118 : nsgfa => basis_set_a%nsgf_set
3195 20118 : rpgfa => basis_set_a%pgf_radius
3196 20118 : set_radius_a => basis_set_a%set_radius
3197 20118 : sphi_a => basis_set_a%sphi
3198 20118 : zeta => basis_set_a%zet
3199 : ! basis jkind
3200 20118 : first_sgfb => basis_set_b%first_sgf
3201 20118 : lb_max => basis_set_b%lmax
3202 20118 : lb_min => basis_set_b%lmin
3203 20118 : npgfb => basis_set_b%npgf
3204 20118 : nsetb = basis_set_b%nset
3205 20118 : nsgfb => basis_set_b%nsgf_set
3206 20118 : rpgfb => basis_set_b%pgf_radius
3207 20118 : set_radius_b => basis_set_b%set_radius
3208 20118 : sphi_b => basis_set_b%sphi
3209 20118 : zetb => basis_set_b%zet
3210 :
3211 20118 : IF (inode == 1) last_jatom = 0
3212 :
3213 : ! this guarentees minimum image convention
3214 : ! anything else would not make sense
3215 20118 : IF (jatom == last_jatom) THEN
3216 : CYCLE
3217 : END IF
3218 :
3219 2106 : last_jatom = jatom
3220 :
3221 2106 : irow = iatom
3222 2106 : icol = jatom
3223 :
3224 8424 : DO i = 1, 3
3225 27378 : DO j = 1, 3
3226 18954 : NULLIFY (mint(i, j)%block)
3227 : CALL dbcsr_get_block_p(matrix=difdip(i, j)%matrix, &
3228 : row=irow, col=icol, BLOCK=mint(i, j)%block, &
3229 18954 : found=found)
3230 25272 : CPASSERT(found)
3231 : END DO
3232 : END DO
3233 :
3234 8424 : ra(:) = particle_set(iatom)%r(:)
3235 8424 : rb(:) = particle_set(jatom)%r(:)
3236 2106 : rab(:) = pbc(rb, ra, cell)
3237 8424 : rac(:) = pbc(ra - rc, cell)
3238 8424 : rbc(:) = pbc(rb - rc, cell)
3239 2106 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
3240 :
3241 4680 : DO iset = 1, nseta
3242 2106 : ncoa = npgfa(iset)*ncoset(la_max(iset))
3243 2106 : sgfa = first_sgfa(1, iset)
3244 24330 : DO jset = 1, nsetb
3245 2106 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
3246 2106 : ncob = npgfb(jset)*ncoset(lb_max(jset))
3247 2106 : sgfb = first_sgfb(1, jset)
3248 2106 : ldab = MAX(ncoa, ncob)
3249 2106 : lda = ncoset(la_max(iset))*npgfa(iset)
3250 2106 : ldb = ncoset(lb_max(jset))*npgfb(jset)
3251 12636 : ALLOCATE (difmab(lda, ldb, M_dim, 3))
3252 :
3253 : ! Calculate integral (da|r|b)
3254 : CALL diff_momop2(la_max(iset), npgfa(iset), zeta(:, iset), &
3255 : rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
3256 : zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
3257 2106 : difmab, deltaR=deltaR, iatom=iatom, jatom=jatom)
3258 :
3259 : ! *** Contraction step ***
3260 :
3261 8424 : DO idir = 1, 3 ! derivative of AO function
3262 27378 : DO j = 1, 3 ! position operator r_j
3263 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
3264 : 1.0_dp, difmab(1, 1, j, idir), SIZE(difmab, 1), &
3265 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
3266 18954 : 0.0_dp, work(1, 1), SIZE(work, 1))
3267 :
3268 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
3269 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
3270 : work(1, 1), SIZE(work, 1), &
3271 : 1.0_dp, mint(j, idir)%block(sgfa, sgfb), &
3272 25272 : SIZE(mint(j, idir)%block, 1))
3273 : END DO !j
3274 : END DO !idir
3275 4212 : DEALLOCATE (difmab)
3276 : END DO !jset
3277 : END DO !iset
3278 : END DO!iterator
3279 :
3280 468 : CALL neighbor_list_iterator_release(nl_iterator)
3281 :
3282 1872 : DO i = 1, 3
3283 6084 : DO j = 1, 3
3284 5616 : NULLIFY (mint(i, j)%block)
3285 : END DO
3286 : END DO
3287 :
3288 468 : DEALLOCATE (mab, difmab2, basis_set_list, work, mint, mint2)
3289 :
3290 468 : CALL timestop(handle)
3291 1404 : END SUBROUTINE dipole_deriv_ao
3292 :
3293 : END MODULE qs_moments
|
