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