Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \par History
10 : !> created 07.2005
11 : !> \author MI (07.2005)
12 : ! **************************************************************************************************
13 : MODULE qs_operators_ao
14 : USE ai_angmom, ONLY: angmom
15 : USE ai_moments, ONLY: contract_cossin,&
16 : cossin,&
17 : diff_momop,&
18 : moment
19 : USE ai_overlap, ONLY: overlap_ab
20 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
21 : gto_basis_set_type
22 : USE block_p_types, ONLY: block_p_type
23 : USE cell_types, ONLY: cell_type,&
24 : pbc
25 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
26 : dbcsr_get_matrix_type,&
27 : dbcsr_has_symmetry,&
28 : dbcsr_p_type,&
29 : dbcsr_set,&
30 : dbcsr_type_antisymmetric,&
31 : dbcsr_type_no_symmetry
32 : USE kinds, ONLY: default_string_length,&
33 : dp
34 : USE orbital_pointers, ONLY: init_orbital_pointers,&
35 : ncoset
36 : USE particle_types, ONLY: particle_type
37 : USE qs_environment_types, ONLY: get_qs_env,&
38 : qs_environment_type
39 : USE qs_kind_types, ONLY: get_qs_kind,&
40 : get_qs_kind_set,&
41 : qs_kind_type
42 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
43 : neighbor_list_iterate,&
44 : neighbor_list_iterator_create,&
45 : neighbor_list_iterator_p_type,&
46 : neighbor_list_iterator_release,&
47 : neighbor_list_set_p_type
48 : #include "./base/base_uses.f90"
49 :
50 : IMPLICIT NONE
51 : PRIVATE
52 :
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_operators_ao'
54 :
55 : ! *** Public subroutines ***
56 :
57 : PUBLIC :: build_exp_ikr_matrix
58 : PUBLIC :: rRc_xyz_ao, rRc_xyz_der_ao
59 : PUBLIC :: build_lin_mom_matrix, build_ang_mom_matrix
60 :
61 : CONTAINS
62 :
63 : ! **************************************************************************************************
64 : !> \brief Build real and imaginary AO matrices for exp(i*k*r).
65 : !> \param qs_env ...
66 : !> \param op_sm_set op_sm_set(1, i) is the cosine matrix and op_sm_set(2, i) is the sine matrix
67 : !> for kvec(:, i)
68 : !> \param kvec Cartesian wave vectors, one column for each output matrix pair
69 : !> \param sab_orb_external optional neighbor list used instead of the default orbital list
70 : !> \param basis_type optional basis-set name
71 : !> \param force_periodic temporarily use all three periodic directions
72 : !> \param cell_external optional cell used for coordinate folding and periodicity
73 : !> \author CP2K developers
74 : !> \note The cosine and sine matrices must be preallocated with matching symmetric
75 : !> DBCSR topology. The matrices are overwritten from zero and only canonical
76 : !> atom-pair blocks are written.
77 : ! **************************************************************************************************
78 13960 : SUBROUTINE build_exp_ikr_matrix(qs_env, op_sm_set, kvec, sab_orb_external, basis_type, force_periodic, cell_external)
79 :
80 : TYPE(qs_environment_type), POINTER :: qs_env
81 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: op_sm_set
82 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: kvec
83 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
84 : OPTIONAL, POINTER :: sab_orb_external
85 : CHARACTER(LEN=*), OPTIONAL :: basis_type
86 : LOGICAL, OPTIONAL :: force_periodic
87 : TYPE(cell_type), OPTIONAL, POINTER :: cell_external
88 :
89 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_exp_ikr_matrix'
90 :
91 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
92 : ldsa, ldsb, ldwork, ncoa, ncob, nkind, nkvec, nseta, nsetb, reim, sgfa, sgfb
93 : INTEGER, DIMENSION(3) :: perd0
94 13960 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
95 13960 : npgfb, nsgfa, nsgfb
96 13960 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
97 : LOGICAL :: found, my_force_periodic, new_atom_b
98 : REAL(KIND=dp) :: dab
99 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb
100 13960 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
101 13960 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cosab, rpgfa, rpgfb, sinab, sphi_a, &
102 13960 : sphi_b, work, zeta, zetb
103 13960 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: op_cossin
104 : TYPE(cell_type), POINTER :: cell
105 13960 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
106 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
107 : TYPE(neighbor_list_iterator_p_type), &
108 13960 : DIMENSION(:), POINTER :: nl_iterator
109 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
110 13960 : POINTER :: sab_orb
111 13960 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
112 13960 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
113 : TYPE(qs_kind_type), POINTER :: qs_kind
114 :
115 13960 : CALL timeset(routineN, handle)
116 :
117 13960 : CPASSERT(SIZE(kvec, 1) == 3)
118 13960 : nkvec = SIZE(kvec, 2)
119 13960 : CPASSERT(nkvec > 0)
120 13960 : CPASSERT(SIZE(op_sm_set, 1) == 2)
121 13960 : CPASSERT(SIZE(op_sm_set, 2) == nkvec)
122 28964 : DO i = 1, nkvec
123 58972 : DO reim = 1, 2
124 30008 : CPASSERT(ASSOCIATED(op_sm_set(reim, i)%matrix))
125 30008 : CPASSERT(dbcsr_has_symmetry(op_sm_set(reim, i)%matrix))
126 45012 : CALL dbcsr_set(op_sm_set(reim, i)%matrix, 0.0_dp)
127 : END DO
128 : END DO
129 :
130 13960 : NULLIFY (qs_kind, qs_kind_set, particle_set, sab_orb, cell)
131 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
132 13960 : particle_set=particle_set, cell=cell, sab_orb=sab_orb)
133 13960 : CPASSERT(ASSOCIATED(cell))
134 13960 : CPASSERT(ASSOCIATED(qs_kind_set))
135 13960 : CPASSERT(ASSOCIATED(particle_set))
136 13960 : IF (PRESENT(cell_external)) THEN
137 504 : CPASSERT(ASSOCIATED(cell_external))
138 504 : cell => cell_external
139 : END IF
140 :
141 13960 : IF (PRESENT(sab_orb_external)) THEN
142 11160 : CPASSERT(ASSOCIATED(sab_orb_external))
143 11160 : sab_orb => sab_orb_external
144 : END IF
145 13960 : CPASSERT(ASSOCIATED(sab_orb))
146 :
147 13960 : my_force_periodic = .FALSE.
148 13960 : IF (PRESENT(force_periodic)) my_force_periodic = force_periodic
149 504 : IF (my_force_periodic) THEN
150 2016 : perd0(:) = cell%perd(:)
151 2016 : cell%perd(:) = 1
152 : END IF
153 :
154 13960 : nkind = SIZE(qs_kind_set)
155 16760 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork, basis_type=basis_type)
156 111680 : ALLOCATE (cosab(ldwork, ldwork), sinab(ldwork, ldwork), work(ldwork, ldwork))
157 :
158 86892 : ALLOCATE (op_cossin(2, nkvec))
159 :
160 69662 : ALLOCATE (basis_set_list(nkind))
161 41742 : DO ikind = 1, nkind
162 27782 : qs_kind => qs_kind_set(ikind)
163 27782 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
164 41742 : IF (ASSOCIATED(basis_set_a)) THEN
165 27782 : basis_set_list(ikind)%gto_basis_set => basis_set_a
166 : ELSE
167 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
168 : END IF
169 : END DO
170 :
171 13960 : last_jatom = 0
172 13960 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
173 309559 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
174 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
175 295599 : iatom=iatom, jatom=jatom, r=rab)
176 295599 : basis_set_a => basis_set_list(ikind)%gto_basis_set
177 295599 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
178 295599 : basis_set_b => basis_set_list(jkind)%gto_basis_set
179 295599 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
180 :
181 295599 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
182 1182396 : rb(:) = ra(:) + rab(:)
183 :
184 295599 : first_sgfa => basis_set_a%first_sgf
185 295599 : la_max => basis_set_a%lmax
186 295599 : la_min => basis_set_a%lmin
187 295599 : npgfa => basis_set_a%npgf
188 295599 : nsgfa => basis_set_a%nsgf_set
189 295599 : rpgfa => basis_set_a%pgf_radius
190 295599 : set_radius_a => basis_set_a%set_radius
191 295599 : sphi_a => basis_set_a%sphi
192 295599 : zeta => basis_set_a%zet
193 :
194 295599 : first_sgfb => basis_set_b%first_sgf
195 295599 : lb_max => basis_set_b%lmax
196 295599 : lb_min => basis_set_b%lmin
197 295599 : npgfb => basis_set_b%npgf
198 295599 : nsgfb => basis_set_b%nsgf_set
199 295599 : rpgfb => basis_set_b%pgf_radius
200 295599 : set_radius_b => basis_set_b%set_radius
201 295599 : sphi_b => basis_set_b%sphi
202 295599 : zetb => basis_set_b%zet
203 :
204 295599 : nseta = basis_set_a%nset
205 295599 : nsetb = basis_set_b%nset
206 295599 : ldsa = SIZE(sphi_a, 1)
207 295599 : ldsb = SIZE(sphi_b, 1)
208 295599 : IF (inode == 1) last_jatom = 0
209 :
210 295599 : IF (jatom /= last_jatom) THEN
211 : new_atom_b = .TRUE.
212 : last_jatom = jatom
213 : ELSE
214 : new_atom_b = .FALSE.
215 : END IF
216 :
217 : IF (new_atom_b) THEN
218 62003 : IF (iatom <= jatom) THEN
219 38072 : irow = iatom
220 38072 : icol = jatom
221 : ELSE
222 23931 : irow = jatom
223 23931 : icol = iatom
224 : END IF
225 :
226 162356 : DO i = 1, nkvec
227 301059 : DO reim = 1, 2
228 200706 : NULLIFY (op_cossin(reim, i)%block)
229 : CALL dbcsr_get_block_p(matrix=op_sm_set(reim, i)%matrix, &
230 301059 : row=irow, col=icol, block=op_cossin(reim, i)%block, found=found)
231 : END DO
232 162356 : IF (ASSOCIATED(op_cossin(1, i)%block) .NEQV. ASSOCIATED(op_cossin(2, i)%block)) THEN
233 0 : CPABORT("cosine and sine blocks should have the same topology")
234 : END IF
235 : END DO
236 : END IF
237 :
238 1182396 : dab = NORM2(rab)
239 1104555 : DO iset = 1, nseta
240 794996 : ncoa = npgfa(iset)*ncoset(la_max(iset))
241 794996 : sgfa = first_sgfa(1, iset)
242 4143076 : DO jset = 1, nsetb
243 3052481 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
244 1094156 : ncob = npgfb(jset)*ncoset(lb_max(jset))
245 1094156 : sgfb = first_sgfb(1, jset)
246 :
247 3296636 : DO i = 1, nkvec
248 1407484 : IF (.NOT. ASSOCIATED(op_cossin(1, i)%block)) CYCLE
249 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), &
250 : la_min(iset), lb_max(jset), npgfb(jset), zetb(:, jset), &
251 1407484 : rpgfb(:, jset), lb_min(jset), ra, rb, kvec(:, i), cosab, sinab)
252 : CALL contract_cossin(op_cossin(1, i)%block, op_cossin(2, i)%block, &
253 : iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
254 : jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
255 4459965 : cosab, sinab, ldwork, work, ldwork)
256 : END DO
257 : END DO
258 : END DO
259 : END DO
260 13960 : CALL neighbor_list_iterator_release(nl_iterator)
261 :
262 13960 : DEALLOCATE (op_cossin, cosab, sinab, work, basis_set_list)
263 :
264 15472 : IF (my_force_periodic) cell%perd(:) = perd0(:)
265 :
266 13960 : CALL timestop(handle)
267 13960 : END SUBROUTINE build_exp_ikr_matrix
268 :
269 : ! **************************************************************************************************
270 : !> \brief Calculation of the linear momentum matrix <mu|∂|nu> over
271 : !> Cartesian Gaussian functions.
272 : !> \param qs_env ...
273 : !> \param matrix ...
274 : !> \param minimum_image take into account only the first neighbors in the lists
275 : !> \date 27.02.2009
276 : !> \author VW
277 : !> \version 1.0
278 : ! **************************************************************************************************
279 1296 : SUBROUTINE build_lin_mom_matrix(qs_env, matrix, minimum_image)
280 :
281 : TYPE(qs_environment_type), POINTER :: qs_env
282 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
283 : LOGICAL, INTENT(IN), OPTIONAL :: minimum_image
284 :
285 : CHARACTER(len=*), PARAMETER :: routineN = 'build_lin_mom_matrix'
286 :
287 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
288 : ldai, maxco, maxlgto, maxsgf, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
289 1296 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nsgfa, &
290 1296 : nsgfb
291 1296 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
292 : LOGICAL :: do_symmetric, found, my_minimum_image, &
293 : new_atom_b
294 : REAL(KIND=dp) :: alpha, dab, Lxo2, Lyo2, Lzo2
295 1296 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: rr_work
296 1296 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
297 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
298 : TARGET :: intab
299 : REAL(KIND=dp), DIMENSION(3) :: rab
300 1296 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
301 1296 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
302 1296 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
303 : TYPE(cell_type), POINTER :: cell
304 1296 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
305 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
306 : TYPE(neighbor_list_iterator_p_type), &
307 1296 : DIMENSION(:), POINTER :: nl_iterator
308 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
309 1296 : POINTER :: sab_nl
310 1296 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
311 : TYPE(qs_kind_type), POINTER :: qs_kind
312 :
313 1296 : CALL timeset(routineN, handle)
314 :
315 1296 : NULLIFY (cell, sab_nl, qs_kind_set)
316 :
317 : CALL get_qs_env(qs_env=qs_env, &
318 : qs_kind_set=qs_kind_set, &
319 1296 : cell=cell)
320 :
321 1296 : nkind = SIZE(qs_kind_set)
322 :
323 1296 : my_minimum_image = .FALSE.
324 1296 : IF (PRESENT(minimum_image)) THEN
325 44 : my_minimum_image = minimum_image
326 176 : Lxo2 = NORM2(cell%hmat(:, 1))/2.0_dp
327 176 : Lyo2 = NORM2(cell%hmat(:, 2))/2.0_dp
328 176 : Lzo2 = NORM2(cell%hmat(:, 3))/2.0_dp
329 : END IF
330 :
331 : ! Take into account the symmetry of the input matrix
332 1296 : do_symmetric = dbcsr_has_symmetry(matrix(1)%matrix)
333 1296 : IF (do_symmetric) THEN
334 1294 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
335 : ELSE
336 2 : CALL get_qs_env(qs_env=qs_env, sab_all=sab_nl)
337 : END IF
338 : ! *** Allocate work storage ***
339 :
340 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
341 : maxco=maxco, &
342 : maxlgto=maxlgto, &
343 1296 : maxsgf=maxsgf)
344 :
345 1296 : ldai = ncoset(maxlgto + 1)
346 1296 : CALL init_orbital_pointers(ldai)
347 :
348 16848 : ALLOCATE (rr_work(ldai*ldai*3), intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
349 1296 : rr_work(:) = 0.0_dp
350 1296 : intab(:, :, :) = 0.0_dp
351 1296 : work(:, :) = 0.0_dp
352 :
353 6132 : ALLOCATE (basis_set_list(nkind))
354 3540 : DO ikind = 1, nkind
355 2244 : qs_kind => qs_kind_set(ikind)
356 2244 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
357 3540 : IF (ASSOCIATED(basis_set_a)) THEN
358 2244 : basis_set_list(ikind)%gto_basis_set => basis_set_a
359 : ELSE
360 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
361 : END IF
362 : END DO
363 1296 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
364 75246 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
365 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
366 73950 : iatom=iatom, jatom=jatom, r=rab)
367 73950 : basis_set_a => basis_set_list(ikind)%gto_basis_set
368 73950 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
369 73950 : basis_set_b => basis_set_list(jkind)%gto_basis_set
370 73950 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
371 : ! basis ikind
372 73950 : first_sgfa => basis_set_a%first_sgf
373 73950 : la_max => basis_set_a%lmax
374 73950 : npgfa => basis_set_a%npgf
375 73950 : nseta = basis_set_a%nset
376 73950 : nsgfa => basis_set_a%nsgf_set
377 73950 : rpgfa => basis_set_a%pgf_radius
378 73950 : set_radius_a => basis_set_a%set_radius
379 73950 : sphi_a => basis_set_a%sphi
380 73950 : zeta => basis_set_a%zet
381 : ! basis jkind
382 73950 : first_sgfb => basis_set_b%first_sgf
383 73950 : lb_max => basis_set_b%lmax
384 73950 : npgfb => basis_set_b%npgf
385 73950 : nsetb = basis_set_b%nset
386 73950 : nsgfb => basis_set_b%nsgf_set
387 73950 : rpgfb => basis_set_b%pgf_radius
388 73950 : set_radius_b => basis_set_b%set_radius
389 73950 : sphi_b => basis_set_b%sphi
390 73950 : zetb => basis_set_b%zet
391 :
392 73950 : IF (inode == 1) last_jatom = 0
393 :
394 73950 : IF (my_minimum_image) THEN
395 26557 : IF (ANY(ABS(rab(:)) > [Lxo2, Lyo2, Lzo2])) CYCLE
396 : END IF
397 :
398 70486 : IF (jatom /= last_jatom) THEN
399 : new_atom_b = .TRUE.
400 : last_jatom = jatom
401 : ELSE
402 : new_atom_b = .FALSE.
403 : END IF
404 :
405 : IF (new_atom_b) THEN
406 13954 : alpha = 1.0_dp
407 13954 : IF (do_symmetric) THEN
408 13945 : IF (iatom <= jatom) THEN
409 7673 : irow = iatom
410 7673 : icol = jatom
411 : ELSE
412 6272 : irow = jatom
413 6272 : icol = iatom
414 6272 : IF (dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_antisymmetric) THEN
415 12544 : alpha = -1.0_dp
416 : END IF
417 : END IF
418 : ELSE
419 9 : irow = iatom
420 9 : icol = jatom
421 : END IF
422 :
423 55816 : DO i = 1, 3
424 41862 : NULLIFY (integral(i)%block)
425 : CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
426 41862 : row=irow, col=icol, BLOCK=integral(i)%block, found=found)
427 55816 : CPASSERT(ASSOCIATED(INTEGRAL(i)%block))
428 : END DO
429 : END IF
430 :
431 281944 : dab = NORM2(rab)
432 :
433 221378 : DO iset = 1, nseta
434 :
435 149596 : ncoa = npgfa(iset)*ncoset(la_max(iset))
436 149596 : sgfa = first_sgfa(1, iset)
437 :
438 567203 : DO jset = 1, nsetb
439 :
440 347121 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
441 :
442 149387 : ncob = npgfb(jset)*ncoset(lb_max(jset))
443 149387 : sgfb = first_sgfb(1, jset)
444 :
445 : ! *** Calculate the primitive fermi contact integrals ***
446 :
447 : ! Keep the full primitive layout used by the contraction below.
448 : CALL overlap_ab(la_max(iset), 0, npgfa(iset), &
449 : rpgfa(:, iset), zeta(:, iset), &
450 : lb_max(jset), 0, npgfb(jset), &
451 : rpgfb(:, jset), zetb(:, jset), &
452 149387 : rab, dab=intab, rr_work=rr_work)
453 :
454 : ! *** Contraction step ***
455 :
456 747144 : DO i = 1, 3
457 :
458 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
459 : 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
460 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
461 448161 : 0.0_dp, work(1, 1), SIZE(work, 1))
462 :
463 795282 : IF (do_symmetric) THEN
464 448134 : IF (iatom <= jatom) THEN
465 :
466 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
467 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
468 : work(1, 1), SIZE(work, 1), &
469 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
470 276486 : SIZE(integral(i)%block, 1))
471 :
472 : ELSE
473 :
474 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
475 : alpha, work(1, 1), SIZE(work, 1), &
476 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
477 : 1.0_dp, integral(i)%block(sgfb, sgfa), &
478 171648 : SIZE(integral(i)%block, 1))
479 :
480 : END IF
481 : ELSE
482 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
483 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
484 : work(1, 1), SIZE(work, 1), &
485 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
486 27 : SIZE(integral(i)%block, 1))
487 : END IF
488 :
489 : END DO
490 :
491 : END DO
492 :
493 : END DO
494 :
495 : END DO
496 1296 : CALL neighbor_list_iterator_release(nl_iterator)
497 :
498 : ! *** Release work storage ***
499 :
500 1296 : DEALLOCATE (intab, rr_work, work, integral, basis_set_list)
501 :
502 1296 : CALL timestop(handle)
503 :
504 3888 : END SUBROUTINE build_lin_mom_matrix
505 :
506 : ! **************************************************************************************************
507 : !> \brief Calculation of the angular momentum matrix over
508 : !> Cartesian Gaussian functions.
509 : !> \param qs_env ...
510 : !> \param matrix ...
511 : !> \param rc ...
512 : !> \date 27.02.2009
513 : !> \author VW
514 : !> \version 1.0
515 : ! **************************************************************************************************
516 :
517 1250 : SUBROUTINE build_ang_mom_matrix(qs_env, matrix, rc)
518 :
519 : TYPE(qs_environment_type), POINTER :: qs_env
520 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
521 : REAL(dp), DIMENSION(:), INTENT(IN) :: rc
522 :
523 : CHARACTER(len=*), PARAMETER :: routineN = 'build_ang_mom_matrix'
524 :
525 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
526 : maxco, maxsgf, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
527 1250 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
528 1250 : npgfb, nsgfa, nsgfb
529 1250 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
530 : LOGICAL :: found, new_atom_b
531 : REAL(KIND=dp) :: dab
532 1250 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
533 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab
534 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rbc
535 1250 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
536 1250 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
537 1250 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
538 : TYPE(cell_type), POINTER :: cell
539 1250 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
540 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
541 : TYPE(neighbor_list_iterator_p_type), &
542 1250 : DIMENSION(:), POINTER :: nl_iterator
543 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
544 1250 : POINTER :: sab_all
545 1250 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
546 1250 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
547 : TYPE(qs_kind_type), POINTER :: qs_kind
548 :
549 1250 : CALL timeset(routineN, handle)
550 :
551 : CALL get_qs_env(qs_env=qs_env, &
552 : qs_kind_set=qs_kind_set, &
553 : particle_set=particle_set, &
554 : sab_all=sab_all, &
555 1250 : cell=cell)
556 :
557 1250 : nkind = SIZE(qs_kind_set)
558 :
559 : ! *** Allocate work storage ***
560 :
561 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
562 : maxco=maxco, &
563 1250 : maxsgf=maxsgf)
564 :
565 13750 : ALLOCATE (intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
566 1250 : intab(:, :, :) = 0.0_dp
567 1250 : work(:, :) = 0.0_dp
568 :
569 5796 : ALLOCATE (basis_set_list(nkind))
570 3296 : DO ikind = 1, nkind
571 2046 : qs_kind => qs_kind_set(ikind)
572 2046 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
573 3296 : IF (ASSOCIATED(basis_set_a)) THEN
574 2046 : basis_set_list(ikind)%gto_basis_set => basis_set_a
575 : ELSE
576 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
577 : END IF
578 : END DO
579 1250 : CALL neighbor_list_iterator_create(nl_iterator, sab_all)
580 94981 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
581 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
582 93731 : iatom=iatom, jatom=jatom, r=rab)
583 93731 : basis_set_a => basis_set_list(ikind)%gto_basis_set
584 93731 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
585 93731 : basis_set_b => basis_set_list(jkind)%gto_basis_set
586 93731 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
587 93731 : ra = pbc(particle_set(iatom)%r, cell)
588 : ! basis ikind
589 93731 : first_sgfa => basis_set_a%first_sgf
590 93731 : la_max => basis_set_a%lmax
591 93731 : la_min => basis_set_a%lmin
592 93731 : npgfa => basis_set_a%npgf
593 93731 : nseta = basis_set_a%nset
594 93731 : nsgfa => basis_set_a%nsgf_set
595 93731 : rpgfa => basis_set_a%pgf_radius
596 93731 : set_radius_a => basis_set_a%set_radius
597 93731 : sphi_a => basis_set_a%sphi
598 93731 : zeta => basis_set_a%zet
599 : ! basis jkind
600 93731 : first_sgfb => basis_set_b%first_sgf
601 93731 : lb_max => basis_set_b%lmax
602 93731 : lb_min => basis_set_b%lmin
603 93731 : npgfb => basis_set_b%npgf
604 93731 : nsetb = basis_set_b%nset
605 93731 : nsgfb => basis_set_b%nsgf_set
606 93731 : rpgfb => basis_set_b%pgf_radius
607 93731 : set_radius_b => basis_set_b%set_radius
608 93731 : sphi_b => basis_set_b%sphi
609 93731 : zetb => basis_set_b%zet
610 :
611 93731 : IF (inode == 1) last_jatom = 0
612 :
613 93731 : IF (jatom /= last_jatom) THEN
614 : new_atom_b = .TRUE.
615 : last_jatom = jatom
616 : ELSE
617 : new_atom_b = .FALSE.
618 : END IF
619 :
620 : IF (new_atom_b) THEN
621 : !IF (iatom <= jatom) THEN
622 5987 : irow = iatom
623 5987 : icol = jatom
624 : !ELSE
625 : ! irow = jatom
626 : ! icol = iatom
627 : !END IF
628 :
629 23948 : DO i = 1, 3
630 17961 : NULLIFY (INTEGRAL(i)%block)
631 : CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
632 17961 : row=irow, col=icol, BLOCK=INTEGRAL(i)%block, found=found)
633 23948 : CPASSERT(ASSOCIATED(INTEGRAL(i)%block))
634 : END DO
635 : END IF
636 :
637 374924 : dab = NORM2(rab)
638 :
639 312058 : DO iset = 1, nseta
640 :
641 217077 : ncoa = npgfa(iset)*ncoset(la_max(iset))
642 217077 : sgfa = first_sgfa(1, iset)
643 :
644 837307 : DO jset = 1, nsetb
645 :
646 526499 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
647 :
648 : !IF(PRESENT(wancen)) THEN
649 : ! rc = wancen
650 201893 : rac = pbc(rc, ra, cell)
651 807572 : rbc = rac + rab
652 : !ELSE
653 : ! rc(1:3) = rb(1:3)
654 : ! rac(1:3) = -rab(1:3)
655 : ! rbc(1:3) = 0.0_dp
656 : !ENDIF
657 :
658 201893 : ncob = npgfb(jset)*ncoset(lb_max(jset))
659 201893 : sgfb = first_sgfb(1, jset)
660 :
661 : ! *** Calculate the primitive angular momentum integrals ***
662 :
663 : CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
664 201893 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), rac, rbc, intab)
665 :
666 : ! *** Contraction step ***
667 :
668 1024649 : DO i = 1, 3
669 :
670 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
671 : 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
672 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
673 605679 : 0.0_dp, work(1, 1), SIZE(work, 1))
674 :
675 : !IF (iatom <= jatom) THEN
676 :
677 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
678 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
679 : work(1, 1), SIZE(work, 1), &
680 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
681 1132178 : SIZE(integral(i)%block, 1))
682 :
683 : !ELSE
684 : !
685 : ! CALL dgemm("T","N",nsgfb(jset),nsgfa(iset),ncoa,&
686 : ! -1.0_dp,work(1,1),SIZE(work,1),&
687 : ! sphi_a(1,sgfa),SIZE(sphi_a,1),&
688 : ! 1.0_dp,integral(i)%block(sgfb,sgfa),&
689 : ! SIZE(integral(i)%block,1))
690 : !
691 : !ENDIF
692 :
693 : END DO
694 :
695 : END DO
696 :
697 : END DO
698 :
699 : END DO
700 1250 : CALL neighbor_list_iterator_release(nl_iterator)
701 :
702 : ! *** Release work storage ***
703 :
704 1250 : DEALLOCATE (intab, work, integral, basis_set_list)
705 :
706 : ! *** Print the spin orbit matrix, if requested ***
707 :
708 : !IF (BTEST(cp_print_key_should_output(logger%iter_info,&
709 : ! qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM"),cp_p_file)) THEN
710 : ! iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM",&
711 : ! extension=".Log")
712 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(1)%matrix,4,6,qs_env,para_env,output_unit=iw)
713 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(2)%matrix,4,6,qs_env,para_env,output_unit=iw)
714 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(3)%matrix,4,6,qs_env,para_env,output_unit=iw)
715 : ! CALL cp_print_key_finished_output(iw,logger,qs_env%input,&
716 : ! "DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM")
717 : !END IF
718 :
719 1250 : CALL timestop(handle)
720 :
721 2500 : END SUBROUTINE build_ang_mom_matrix
722 :
723 : ! **************************************************************************************************
724 : !> \brief Calculation of the components of the dipole operator in the length form
725 : !> by taking the relative position operator r-Rc, with respect a reference point Rc
726 : !> Probably it does not work for PBC, or maybe yes if the wfn are
727 : !> sufficiently localized
728 : !> The elements of the sparse matrices are the integrals in the
729 : !> basis functions
730 : !> \param op matrix representation of the p operator
731 : !> calculated in terms of the contracted basis functions
732 : !> \param qs_env environment for the lists and the basis sets
733 : !> \param rc reference vector position
734 : !> \param order maximum order of the momentum, for the dipole order = 1, order = -2 for quad only
735 : !> \param minimum_image take into account only the first neighbors in the lists
736 : !> \param soft ...
737 : !> \par History
738 : !> 03.2006 created [MI]
739 : !> 06.2019 added quarupole only option (A.Bussy)
740 : !> \author MI
741 : ! **************************************************************************************************
742 :
743 60 : SUBROUTINE rRc_xyz_ao(op, qs_env, rc, order, minimum_image, soft)
744 :
745 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op
746 : TYPE(qs_environment_type), POINTER :: qs_env
747 : REAL(dp) :: Rc(3)
748 : INTEGER, INTENT(IN) :: order
749 : LOGICAL, INTENT(IN), OPTIONAL :: minimum_image, soft
750 :
751 : CHARACTER(len=*), PARAMETER :: routineN = 'rRc_xyz_ao'
752 :
753 : CHARACTER(LEN=default_string_length) :: basis_type
754 : INTEGER :: handle, iatom, icol, ikind, imom, inode, irow, iset, jatom, jkind, jset, &
755 : last_jatom, ldab, ldsa, ldsb, ldwork, M_dim, maxl, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
756 : sgfb, smom
757 60 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, npgfa, npgfb, &
758 60 : nsgfa, nsgfb
759 60 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
760 : LOGICAL :: found, my_minimum_image, my_soft, &
761 : new_atom_b
762 : REAL(KIND=dp) :: dab, Lxo2, Lyo2, Lzo2, rab2
763 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
764 60 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
765 60 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, &
766 60 : zeta, zetb
767 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab
768 60 : TYPE(block_p_type), DIMENSION(:), POINTER :: op_dip
769 : TYPE(cell_type), POINTER :: cell
770 60 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
771 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
772 : TYPE(neighbor_list_iterator_p_type), &
773 60 : DIMENSION(:), POINTER :: nl_iterator
774 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
775 60 : POINTER :: sab_orb
776 60 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
777 60 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
778 : TYPE(qs_kind_type), POINTER :: qs_kind
779 :
780 60 : CALL timeset(routineN, handle)
781 :
782 60 : NULLIFY (qs_kind, qs_kind_set)
783 60 : NULLIFY (cell, particle_set)
784 60 : NULLIFY (sab_orb)
785 60 : NULLIFY (mab, op_dip, work)
786 60 : NULLIFY (la_max, la_min, lb_max, npgfa, npgfb, nsgfa, nsgfb)
787 60 : NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
788 :
789 60 : my_soft = .FALSE.
790 60 : IF (PRESENT(soft)) my_soft = soft
791 28 : IF (my_soft) THEN
792 0 : basis_type = "ORB_SOFT"
793 : ELSE
794 60 : basis_type = "ORB"
795 : END IF
796 :
797 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
798 60 : cell=cell, particle_set=particle_set, sab_orb=sab_orb)
799 :
800 60 : nkind = SIZE(qs_kind_set)
801 :
802 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
803 60 : maxco=ldwork, maxlgto=maxl)
804 :
805 60 : my_minimum_image = .FALSE.
806 60 : IF (PRESENT(minimum_image)) THEN
807 58 : my_minimum_image = minimum_image
808 232 : Lxo2 = SQRT(SUM(cell%hmat(:, 1)**2))/2.0_dp
809 232 : Lyo2 = SQRT(SUM(cell%hmat(:, 2)**2))/2.0_dp
810 232 : Lzo2 = SQRT(SUM(cell%hmat(:, 3)**2))/2.0_dp
811 : END IF
812 :
813 60 : ldab = ldwork
814 :
815 60 : smom = 1
816 60 : IF (order == -2) smom = 4
817 60 : M_dim = ncoset(ABS(order)) - 1
818 60 : CPASSERT(M_dim <= SIZE(op, 1))
819 :
820 300 : ALLOCATE (mab(ldab, ldab, 1:M_dim))
821 46032 : mab(1:ldab, 1:ldab, 1:M_dim) = 0.0_dp
822 240 : ALLOCATE (work(ldwork, ldwork))
823 15324 : work(1:ldwork, 1:ldwork) = 0.0_dp
824 360 : ALLOCATE (op_dip(smom:M_dim))
825 :
826 240 : DO imom = smom, M_dim
827 240 : NULLIFY (op_dip(imom)%block)
828 : END DO
829 :
830 282 : ALLOCATE (basis_set_list(nkind))
831 162 : DO ikind = 1, nkind
832 102 : qs_kind => qs_kind_set(ikind)
833 102 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
834 162 : IF (ASSOCIATED(basis_set_a)) THEN
835 102 : basis_set_list(ikind)%gto_basis_set => basis_set_a
836 : ELSE
837 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
838 : END IF
839 : END DO
840 60 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
841 462 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
842 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
843 402 : iatom=iatom, jatom=jatom, r=rab)
844 402 : basis_set_a => basis_set_list(ikind)%gto_basis_set
845 402 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
846 402 : basis_set_b => basis_set_list(jkind)%gto_basis_set
847 402 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
848 402 : ra = pbc(particle_set(iatom)%r, cell)
849 : ! basis ikind
850 402 : first_sgfa => basis_set_a%first_sgf
851 402 : la_max => basis_set_a%lmax
852 402 : la_min => basis_set_a%lmin
853 402 : npgfa => basis_set_a%npgf
854 402 : nseta = basis_set_a%nset
855 402 : nsgfa => basis_set_a%nsgf_set
856 402 : rpgfa => basis_set_a%pgf_radius
857 402 : set_radius_a => basis_set_a%set_radius
858 402 : sphi_a => basis_set_a%sphi
859 402 : zeta => basis_set_a%zet
860 : ! basis jkind
861 402 : first_sgfb => basis_set_b%first_sgf
862 402 : lb_max => basis_set_b%lmax
863 402 : npgfb => basis_set_b%npgf
864 402 : nsetb = basis_set_b%nset
865 402 : nsgfb => basis_set_b%nsgf_set
866 402 : rpgfb => basis_set_b%pgf_radius
867 402 : set_radius_b => basis_set_b%set_radius
868 402 : sphi_b => basis_set_b%sphi
869 402 : zetb => basis_set_b%zet
870 :
871 402 : ldsa = SIZE(sphi_a, 1)
872 402 : ldsb = SIZE(sphi_b, 1)
873 402 : IF (inode == 1) last_jatom = 0
874 :
875 402 : IF (my_minimum_image) THEN
876 121 : IF (ABS(rab(1)) > Lxo2 .OR. ABS(rab(2)) > Lyo2 .OR. ABS(rab(3)) > Lzo2) CYCLE
877 : END IF
878 :
879 1600 : rb = rab + ra
880 :
881 400 : IF (jatom /= last_jatom) THEN
882 : new_atom_b = .TRUE.
883 : last_jatom = jatom
884 : ELSE
885 : new_atom_b = .FALSE.
886 : END IF
887 :
888 : IF (new_atom_b) THEN
889 151 : IF (iatom <= jatom) THEN
890 100 : irow = iatom
891 100 : icol = jatom
892 : ELSE
893 51 : irow = jatom
894 51 : icol = iatom
895 : END IF
896 :
897 604 : DO imom = smom, M_dim
898 453 : NULLIFY (op_dip(imom)%block)
899 : CALL dbcsr_get_block_p(matrix=op(imom)%matrix, &
900 453 : row=irow, col=icol, block=op_dip(imom)%block, found=found)
901 604 : CPASSERT(ASSOCIATED(op_dip(imom)%block))
902 : END DO ! imom
903 : END IF ! new_atom_b
904 :
905 400 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
906 400 : dab = SQRT(rab2)
907 :
908 1412 : DO iset = 1, nseta
909 :
910 952 : ncoa = npgfa(iset)*ncoset(la_max(iset))
911 952 : sgfa = first_sgfa(1, iset)
912 :
913 3856 : DO jset = 1, nsetb
914 :
915 2502 : ncob = npgfb(jset)*ncoset(lb_max(jset))
916 2502 : sgfb = first_sgfb(1, jset)
917 :
918 3454 : IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN
919 :
920 1921 : rac = pbc(rc, ra, cell)
921 1921 : rbc = pbc(rc, rb, cell)
922 :
923 : ! *** Calculate the primitive overlap integrals ***
924 : CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
925 : rpgfa(:, iset), la_min(iset), &
926 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), &
927 1921 : ABS(order), rac, rbc, mab)
928 :
929 7684 : DO imom = smom, M_dim
930 : ! *** Contraction ***
931 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
932 : 1.0_dp, mab(1, 1, imom), ldab, sphi_b(1, sgfb), ldsb, &
933 5763 : 0.0_dp, work(1, 1), ldwork)
934 7684 : IF (iatom <= jatom) THEN
935 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
936 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
937 : work(1, 1), ldwork, &
938 : 1.0_dp, op_dip(imom)%block(sgfa, sgfb), &
939 4215 : SIZE(op_dip(imom)%block, 1))
940 : ELSE
941 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
942 : 1.0_dp, work(1, 1), ldwork, &
943 : sphi_a(1, sgfa), ldsa, &
944 : 1.0_dp, op_dip(imom)%block(sgfb, sgfa), &
945 1548 : SIZE(op_dip(imom)%block, 1))
946 : END IF
947 :
948 : END DO ! imom
949 : END IF ! >= dab
950 :
951 : END DO ! jset
952 :
953 : END DO ! iset
954 :
955 : END DO
956 60 : CALL neighbor_list_iterator_release(nl_iterator)
957 :
958 240 : DO imom = smom, M_dim
959 240 : NULLIFY (op_dip(imom)%block)
960 : END DO
961 60 : DEALLOCATE (op_dip)
962 :
963 60 : DEALLOCATE (mab, work, basis_set_list)
964 :
965 60 : CALL timestop(handle)
966 :
967 120 : END SUBROUTINE rRc_xyz_ao
968 :
969 : ! **************************************************************************************************
970 : !> \brief Calculation of the multipole operators integrals
971 : !> and of its derivatives of the type
972 : !> [\mu | op | d(\nu)/dR(\nu)]-[d(\mu)/dR(\mu)| op | \nu]
973 : !> by taking the relative position operator r-Rc, with respect a reference point Rc
974 : !> The derivative are with respect to the primitive position,
975 : !> The multipole operator is symmetric and if it does not depend on R(\mu) or R(\nu)
976 : !> therefore [\mu | op | d(\nu)/dR(\nu)] = -[d(\mu)/dR(\mu)| op | \nu]
977 : !> [\mu|op|d(\nu)/dR]-[d(\mu)/dR|op|\nu]=2[\mu|op|d(\nu)/dR]
978 : !> When it is not the case a correction term is needed
979 : !>
980 : !> The momentum operator [\mu|M|\nu] is symmetric, the number of components is
981 : !> determined by the order: 3 for order 1 (x,y,x), 9 for order 2(xx,xy,xz,yy,yz,zz)
982 : !> The derivative of the type [\mu | op | d(\nu)/dR_i(\nu)], where
983 : !> i indicates the cartesian direction, is antisymmetric only when
984 : !> the no component M =(r_i) or (r_i r_j) is in the same cartesian
985 : !> direction of the derivative, indeed
986 : !> d([\mu|M|\nu])/dr_i = [d(\mu)/dr_i|M|\nu] + [\mu|M|d(\nu)/dr_i] + [\mu |d(M)/dr_i|\nu]
987 : !> d([\mu|M|\nu])/dr_i = -[d(\mu)/dR_i(\mu)|M|\nu] -[\mu|M|d(\nu)/dR_i(\nu)] + [\mu |d(M)/dr_i|\nu]
988 : !> Therefore we cannot use an antisymmetric matrix
989 : !>
990 : !> The same holds for the derivative with respect to the electronic position r
991 : !> taking into account that [\mu|op|d(\nu)/dR] = -[\mu|op|d(\nu)/dr]
992 : !> \param op matrix representation of the p operator
993 : !> calculated in terms of the contracted basis functions
994 : !> \param op_der ...
995 : !> \param qs_env environment for the lists and the basis sets
996 : !> \param rc reference vector position
997 : !> \param order maximum order of the momentum, for the dipole order = 1
998 : !> \param minimum_image take into account only the first neighbors in the lists
999 : !> \param soft ...
1000 : !> \par History
1001 : !> 03.2006 created [MI]
1002 : !> \author MI
1003 : !> \note
1004 : !> Probably it does not work for PBC, or maybe yes if the wfn are
1005 : !> sufficiently localized
1006 : !> The elements of the sparse matrices are the integrals in the
1007 : !> basis functions
1008 : ! **************************************************************************************************
1009 3750 : SUBROUTINE rRc_xyz_der_ao(op, op_der, qs_env, rc, order, minimum_image, soft)
1010 :
1011 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op
1012 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_der
1013 : TYPE(qs_environment_type), POINTER :: qs_env
1014 : REAL(dp) :: Rc(3)
1015 : INTEGER, INTENT(IN) :: order
1016 : LOGICAL, INTENT(IN), OPTIONAL :: minimum_image, soft
1017 :
1018 : CHARACTER(len=*), PARAMETER :: routineN = 'rRc_xyz_der_ao'
1019 :
1020 : CHARACTER(LEN=default_string_length) :: basis_type
1021 : INTEGER :: handle, i, iatom, icol, idir, ikind, imom, inode, ipgf, irow, iset, j, jatom, &
1022 : jkind, jpgf, jset, last_jatom, lda_min, ldab, ldb_min, ldsa, ldsb, ldwork, M_dim, maxl, &
1023 : na, nb, ncoa, ncob, nda, ndb, nkind, nseta, nsetb, sgfa, sgfb
1024 3750 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1025 3750 : npgfb, nsgfa, nsgfb
1026 3750 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1027 : LOGICAL :: my_minimum_image, my_soft, new_atom_b, &
1028 : op_der_found, op_found
1029 : REAL(KIND=dp) :: alpha, alpha_der, dab, Lxo2, Lyo2, Lzo2, &
1030 : rab2
1031 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
1032 3750 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
1033 3750 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, &
1034 3750 : zeta, zetb
1035 3750 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab, mab_tmp
1036 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: difmab
1037 3750 : TYPE(block_p_type), DIMENSION(:), POINTER :: op_dip
1038 3750 : TYPE(block_p_type), DIMENSION(:, :), POINTER :: op_dip_der
1039 : TYPE(cell_type), POINTER :: cell
1040 3750 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1041 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1042 : TYPE(neighbor_list_iterator_p_type), &
1043 3750 : DIMENSION(:), POINTER :: nl_iterator
1044 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1045 3750 : POINTER :: sab_all
1046 3750 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1047 3750 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1048 : TYPE(qs_kind_type), POINTER :: qs_kind
1049 :
1050 3750 : CALL timeset(routineN, handle)
1051 :
1052 3750 : CPASSERT(ASSOCIATED(op))
1053 3750 : CPASSERT(ASSOCIATED(op_der))
1054 : !IF(.NOT.op_sm_der(1,1)%matrix%symmetry=="none") THEN
1055 3750 : CPASSERT(dbcsr_get_matrix_type(op_der(1, 1)%matrix) == dbcsr_type_no_symmetry)
1056 :
1057 3750 : NULLIFY (qs_kind, qs_kind_set)
1058 3750 : NULLIFY (cell, particle_set)
1059 3750 : NULLIFY (sab_all)
1060 3750 : NULLIFY (difmab, mab, mab_tmp)
1061 3750 : NULLIFY (op_dip, op_dip_der, work)
1062 3750 : NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb)
1063 3750 : NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
1064 :
1065 3750 : my_soft = .FALSE.
1066 3750 : IF (PRESENT(soft)) my_soft = soft
1067 3750 : IF (my_soft) THEN
1068 2022 : basis_type = "ORB_SOFT"
1069 : ELSE
1070 1728 : basis_type = "ORB"
1071 : END IF
1072 :
1073 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
1074 : cell=cell, particle_set=particle_set, &
1075 3750 : sab_all=sab_all)
1076 :
1077 3750 : nkind = SIZE(qs_kind_set)
1078 :
1079 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1080 3750 : maxco=ldwork, maxlgto=maxl)
1081 :
1082 3750 : my_minimum_image = .FALSE.
1083 3750 : IF (PRESENT(minimum_image)) THEN
1084 3750 : my_minimum_image = minimum_image
1085 15000 : Lxo2 = SQRT(SUM(cell%hmat(:, 1)**2))/2.0_dp
1086 15000 : Lyo2 = SQRT(SUM(cell%hmat(:, 2)**2))/2.0_dp
1087 15000 : Lzo2 = SQRT(SUM(cell%hmat(:, 3)**2))/2.0_dp
1088 : END IF
1089 :
1090 3750 : ldab = ldwork
1091 :
1092 3750 : M_dim = ncoset(order) - 1
1093 3750 : CPASSERT(M_dim <= SIZE(op, 1))
1094 :
1095 18750 : ALLOCATE (mab(ldab, ldab, M_dim))
1096 6053640 : mab(1:ldab, 1:ldab, 1:M_dim) = 0.0_dp
1097 22500 : ALLOCATE (difmab(ldab, ldab, M_dim, 3))
1098 18164670 : difmab(1:ldab, 1:ldab, 1:M_dim, 1:3) = 0.0_dp
1099 :
1100 15000 : ALLOCATE (work(ldwork, ldwork))
1101 672210 : work(1:ldwork, 1:ldwork) = 0.0_dp
1102 45000 : ALLOCATE (op_dip(M_dim))
1103 123750 : ALLOCATE (op_dip_der(M_dim, 3))
1104 :
1105 37500 : DO imom = 1, M_dim
1106 33750 : NULLIFY (op_dip(imom)%block)
1107 138750 : DO i = 1, 3
1108 135000 : NULLIFY (op_dip_der(imom, i)%block)
1109 : END DO
1110 : END DO
1111 :
1112 17388 : ALLOCATE (basis_set_list(nkind))
1113 9888 : DO ikind = 1, nkind
1114 6138 : qs_kind => qs_kind_set(ikind)
1115 6138 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1116 9888 : IF (ASSOCIATED(basis_set_a)) THEN
1117 6138 : basis_set_list(ikind)%gto_basis_set => basis_set_a
1118 : ELSE
1119 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
1120 : END IF
1121 : END DO
1122 3750 : CALL neighbor_list_iterator_create(nl_iterator, sab_all)
1123 284943 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1124 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1125 281193 : iatom=iatom, jatom=jatom, r=rab)
1126 281193 : basis_set_a => basis_set_list(ikind)%gto_basis_set
1127 281193 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
1128 281193 : basis_set_b => basis_set_list(jkind)%gto_basis_set
1129 281193 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
1130 281193 : ra = pbc(particle_set(iatom)%r, cell)
1131 : ! basis ikind
1132 281193 : first_sgfa => basis_set_a%first_sgf
1133 281193 : la_max => basis_set_a%lmax
1134 281193 : la_min => basis_set_a%lmin
1135 281193 : npgfa => basis_set_a%npgf
1136 281193 : nseta = basis_set_a%nset
1137 281193 : nsgfa => basis_set_a%nsgf_set
1138 281193 : rpgfa => basis_set_a%pgf_radius
1139 281193 : set_radius_a => basis_set_a%set_radius
1140 281193 : sphi_a => basis_set_a%sphi
1141 281193 : zeta => basis_set_a%zet
1142 : ! basis jkind
1143 281193 : first_sgfb => basis_set_b%first_sgf
1144 281193 : lb_max => basis_set_b%lmax
1145 281193 : lb_min => basis_set_b%lmin
1146 281193 : npgfb => basis_set_b%npgf
1147 281193 : nsetb = basis_set_b%nset
1148 281193 : nsgfb => basis_set_b%nsgf_set
1149 281193 : rpgfb => basis_set_b%pgf_radius
1150 281193 : set_radius_b => basis_set_b%set_radius
1151 281193 : sphi_b => basis_set_b%sphi
1152 281193 : zetb => basis_set_b%zet
1153 :
1154 281193 : ldsa = SIZE(sphi_a, 1)
1155 281193 : IF (ldsa == 0) CYCLE
1156 281172 : ldsb = SIZE(sphi_b, 1)
1157 281172 : IF (ldsb == 0) CYCLE
1158 281172 : IF (inode == 1) last_jatom = 0
1159 :
1160 281172 : IF (my_minimum_image) THEN
1161 0 : IF (ABS(rab(1)) > Lxo2 .OR. ABS(rab(2)) > Lyo2 .OR. ABS(rab(3)) > Lzo2) CYCLE
1162 : END IF
1163 :
1164 1124688 : rb = rab + ra
1165 :
1166 281172 : IF (jatom /= last_jatom) THEN
1167 : new_atom_b = .TRUE.
1168 : last_jatom = jatom
1169 : ELSE
1170 : new_atom_b = .FALSE.
1171 : END IF
1172 :
1173 : IF (new_atom_b) THEN
1174 17958 : irow = iatom
1175 17958 : icol = jatom
1176 17958 : alpha_der = 2.0_dp
1177 :
1178 179580 : DO imom = 1, M_dim
1179 161622 : NULLIFY (op_dip(imom)%block)
1180 : CALL dbcsr_get_block_p(matrix=op(imom)%matrix, &
1181 : row=irow, col=icol, &
1182 : block=op_dip(imom)%block, &
1183 161622 : found=op_found)
1184 161622 : CPASSERT(op_found .AND. ASSOCIATED(op_dip(imom)%block))
1185 826068 : DO idir = 1, 3
1186 484866 : NULLIFY (op_dip_der(imom, idir)%block)
1187 : CALL dbcsr_get_block_p(matrix=op_der(imom, idir)%matrix, &
1188 : row=irow, col=icol, &
1189 : block=op_dip_der(imom, idir)%block, &
1190 484866 : found=op_der_found)
1191 646488 : CPASSERT(op_der_found .AND. ASSOCIATED(op_dip_der(imom, idir)%block))
1192 : END DO ! idir
1193 : END DO ! imom
1194 : END IF ! new_atom_b
1195 :
1196 281172 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1197 281172 : dab = SQRT(rab2)
1198 :
1199 936090 : DO iset = 1, nseta
1200 :
1201 651168 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1202 651168 : sgfa = first_sgfa(1, iset)
1203 :
1204 2511669 : DO jset = 1, nsetb
1205 :
1206 1579308 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1207 1579308 : sgfb = first_sgfb(1, jset)
1208 :
1209 2230476 : IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN
1210 :
1211 596982 : rac = pbc(rc, ra, cell)
1212 2387928 : rbc = rac + rab
1213 : ! rac = pbc(rc,ra,cell)
1214 : ! rbc = pbc(rc,rb,cell)
1215 :
1216 : ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
1217 2958078 : npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(order) - 1))
1218 :
1219 596982 : lda_min = MAX(0, la_min(iset) - 1)
1220 596982 : ldb_min = MAX(0, lb_min(jset) - 1)
1221 : ! *** Calculate the primitive overlap integrals ***
1222 : CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), &
1223 : rpgfa(:, iset), lda_min, &
1224 : lb_max(jset) + 1, npgfb(jset), zetb(:, jset), rpgfb(:, jset), &
1225 596982 : order, rac, rbc, mab_tmp)
1226 :
1227 : ! *** Calculate the derivatives
1228 : CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), &
1229 : rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
1230 : zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
1231 596982 : difmab, mab_ext=mab_tmp)
1232 :
1233 : ! Contract and copy in the sparse matrix
1234 874316388 : mab = 0.0_dp
1235 5969820 : DO imom = 1, M_dim
1236 5372838 : na = 0
1237 5372838 : nda = 0
1238 17098830 : DO ipgf = 1, npgfa(iset)
1239 11725992 : nb = 0
1240 11725992 : ndb = 0
1241 42864039 : DO jpgf = 1, npgfb(jset)
1242 127844055 : DO j = 1, ncoset(lb_max(jset))
1243 486284148 : DO i = 1, ncoset(la_max(iset))
1244 455146101 : mab(i + na, j + nb, imom) = mab_tmp(i + nda, j + ndb, imom)
1245 : END DO ! i
1246 : END DO ! j
1247 31138047 : nb = nb + ncoset(lb_max(jset))
1248 42864039 : ndb = ndb + ncoset(lb_max(jset) + 1)
1249 : END DO ! jpgf
1250 11725992 : na = na + ncoset(la_max(iset))
1251 17098830 : nda = nda + ncoset(la_max(iset) + 1)
1252 : END DO ! ipgf
1253 :
1254 : ! *** Contraction ***
1255 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1256 : 1.0_dp, mab(1, 1, imom), ldab, sphi_b(1, sgfb), ldsb, &
1257 5372838 : 0.0_dp, work(1, 1), ldwork)
1258 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1259 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
1260 : work(1, 1), ldwork, &
1261 : 1.0_dp, op_dip(imom)%block(sgfa, sgfb), &
1262 5372838 : SIZE(op_dip(imom)%block, 1))
1263 :
1264 5372838 : alpha = -1.0_dp !-alpha_der
1265 22088334 : DO idir = 1, 3
1266 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1267 : alpha, difmab(1, 1, imom, idir), ldab, sphi_b(1, sgfb), ldsb, &
1268 16118514 : 0.0_dp, work(1, 1), ldwork)
1269 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1270 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
1271 : work(1, 1), ldwork, &
1272 : 1.0_dp, op_dip_der(imom, idir)%block(sgfa, sgfb), &
1273 21491352 : SIZE(op_dip_der(imom, idir)%block, 1))
1274 :
1275 : END DO ! idir
1276 :
1277 : END DO ! imom
1278 :
1279 596982 : DEALLOCATE (mab_tmp)
1280 : END IF ! >= dab
1281 :
1282 : END DO ! jset
1283 :
1284 : END DO ! iset
1285 :
1286 : END DO
1287 3750 : CALL neighbor_list_iterator_release(nl_iterator)
1288 :
1289 15000 : DO i = 1, 3
1290 15000 : NULLIFY (op_dip(i)%block)
1291 : END DO
1292 3750 : DEALLOCATE (op_dip, op_dip_der)
1293 :
1294 3750 : DEALLOCATE (mab, difmab, work, basis_set_list)
1295 :
1296 3750 : CALL timestop(handle)
1297 :
1298 7500 : END SUBROUTINE rRc_xyz_der_ao
1299 :
1300 : END MODULE qs_operators_ao
|
