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_expikr,&
16 : expikr
17 : USE ai_overlap, ONLY: overlap_ab
18 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
19 : gto_basis_set_type
20 : USE block_p_types, ONLY: block_p_type
21 : USE cell_types, ONLY: cell_type,&
22 : pbc
23 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
24 : dbcsr_get_matrix_type,&
25 : dbcsr_has_symmetry,&
26 : dbcsr_p_type,&
27 : dbcsr_set,&
28 : dbcsr_type_antisymmetric
29 : USE kinds, ONLY: dp
30 : USE orbital_pointers, ONLY: init_orbital_pointers,&
31 : ncoset
32 : USE particle_types, ONLY: particle_type
33 : USE qs_environment_types, ONLY: get_qs_env,&
34 : qs_environment_type
35 : USE qs_kind_types, ONLY: get_qs_kind,&
36 : get_qs_kind_set,&
37 : qs_kind_type
38 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
39 : neighbor_list_iterate,&
40 : neighbor_list_iterator_create,&
41 : neighbor_list_iterator_p_type,&
42 : neighbor_list_iterator_release,&
43 : neighbor_list_set_p_type
44 : #include "./base/base_uses.f90"
45 :
46 : IMPLICIT NONE
47 : PRIVATE
48 :
49 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_operators_ao'
50 :
51 : ! *** Public subroutines ***
52 :
53 : PUBLIC :: build_exp_ikr_matrix
54 : PUBLIC :: build_lin_mom_matrix, build_ang_mom_matrix
55 :
56 : CONTAINS
57 :
58 : ! **************************************************************************************************
59 : !> \brief Build real and imaginary AO matrices for exp(i*k*r).
60 : !> \param qs_env ...
61 : !> \param op_sm_set op_sm_set(1, i) is the cosine matrix and op_sm_set(2, i) is the sine matrix
62 : !> for kvec(:, i)
63 : !> \param kvec Cartesian wave vectors, one column for each output matrix pair
64 : !> \param sab_orb_external optional neighbor list used instead of the default orbital list
65 : !> \param basis_type optional basis-set name
66 : !> \param force_periodic temporarily use all three periodic directions
67 : !> \param cell_external optional cell used for coordinate folding and periodicity
68 : !> \author CP2K developers
69 : !> \note The cosine and sine matrices must be preallocated with matching symmetric
70 : !> DBCSR topology. The matrices are overwritten from zero and only canonical
71 : !> atom-pair blocks are written.
72 : ! **************************************************************************************************
73 13972 : SUBROUTINE build_exp_ikr_matrix(qs_env, op_sm_set, kvec, sab_orb_external, basis_type, force_periodic, cell_external)
74 :
75 : TYPE(qs_environment_type), POINTER :: qs_env
76 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: op_sm_set
77 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: kvec
78 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
79 : OPTIONAL, POINTER :: sab_orb_external
80 : CHARACTER(LEN=*), OPTIONAL :: basis_type
81 : LOGICAL, OPTIONAL :: force_periodic
82 : TYPE(cell_type), OPTIONAL, POINTER :: cell_external
83 :
84 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_exp_ikr_matrix'
85 :
86 : COMPLEX(KIND=dp), ALLOCATABLE :: expab(:, :)
87 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
88 : ldsa, ldsb, ldwork, ncoa, ncob, nkind, nkvec, nseta, nsetb, reim, sgfa, sgfb
89 : INTEGER, DIMENSION(3) :: perd0
90 13972 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
91 13972 : npgfb, nsgfa, nsgfb
92 13972 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
93 : LOGICAL :: found, my_force_periodic, new_atom_b
94 : REAL(KIND=dp) :: dab
95 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb
96 13972 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
97 13972 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, &
98 13972 : zeta, zetb
99 13972 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: op_cossin
100 : TYPE(cell_type), POINTER :: cell
101 13972 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
102 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
103 : TYPE(neighbor_list_iterator_p_type), &
104 13972 : DIMENSION(:), POINTER :: nl_iterator
105 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
106 13972 : POINTER :: sab_orb
107 13972 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
108 13972 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
109 : TYPE(qs_kind_type), POINTER :: qs_kind
110 :
111 13972 : CALL timeset(routineN, handle)
112 :
113 13972 : CPASSERT(SIZE(kvec, 1) == 3)
114 13972 : nkvec = SIZE(kvec, 2)
115 13972 : CPASSERT(nkvec > 0)
116 13972 : CPASSERT(SIZE(op_sm_set, 1) == 2)
117 13972 : CPASSERT(SIZE(op_sm_set, 2) == nkvec)
118 28968 : DO i = 1, nkvec
119 58960 : DO reim = 1, 2
120 29992 : CPASSERT(ASSOCIATED(op_sm_set(reim, i)%matrix))
121 29992 : CPASSERT(dbcsr_has_symmetry(op_sm_set(reim, i)%matrix))
122 44988 : CALL dbcsr_set(op_sm_set(reim, i)%matrix, 0.0_dp)
123 : END DO
124 : END DO
125 :
126 13972 : NULLIFY (qs_kind, qs_kind_set, particle_set, sab_orb, cell)
127 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
128 13972 : particle_set=particle_set, cell=cell, sab_orb=sab_orb)
129 13972 : CPASSERT(ASSOCIATED(cell))
130 13972 : CPASSERT(ASSOCIATED(qs_kind_set))
131 13972 : CPASSERT(ASSOCIATED(particle_set))
132 13972 : IF (PRESENT(cell_external)) THEN
133 494 : CPASSERT(ASSOCIATED(cell_external))
134 494 : cell => cell_external
135 : END IF
136 :
137 13972 : IF (PRESENT(sab_orb_external)) THEN
138 11160 : CPASSERT(ASSOCIATED(sab_orb_external))
139 11160 : sab_orb => sab_orb_external
140 : END IF
141 13972 : CPASSERT(ASSOCIATED(sab_orb))
142 :
143 13972 : my_force_periodic = .FALSE.
144 13972 : IF (PRESENT(force_periodic)) my_force_periodic = force_periodic
145 55888 : perd0(:) = -1
146 13972 : IF (my_force_periodic) THEN
147 1976 : perd0(:) = cell%perd(:)
148 1976 : cell%perd(:) = 1
149 : END IF
150 :
151 13972 : nkind = SIZE(qs_kind_set)
152 16784 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork, basis_type=basis_type)
153 97804 : ALLOCATE (expab(ldwork, ldwork), work(ldwork, ldwork))
154 :
155 86904 : ALLOCATE (op_cossin(2, nkvec))
156 :
157 69720 : ALLOCATE (basis_set_list(nkind))
158 41776 : DO ikind = 1, nkind
159 27804 : qs_kind => qs_kind_set(ikind)
160 27804 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
161 41776 : IF (ASSOCIATED(basis_set_a)) THEN
162 27804 : basis_set_list(ikind)%gto_basis_set => basis_set_a
163 : ELSE
164 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
165 : END IF
166 : END DO
167 :
168 13972 : last_jatom = 0
169 13972 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
170 310101 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
171 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
172 296129 : iatom=iatom, jatom=jatom, r=rab)
173 296129 : basis_set_a => basis_set_list(ikind)%gto_basis_set
174 296129 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
175 296129 : basis_set_b => basis_set_list(jkind)%gto_basis_set
176 296129 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
177 :
178 296129 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
179 1184516 : rb(:) = ra(:) + rab(:)
180 :
181 296129 : first_sgfa => basis_set_a%first_sgf
182 296129 : la_max => basis_set_a%lmax
183 296129 : la_min => basis_set_a%lmin
184 296129 : npgfa => basis_set_a%npgf
185 296129 : nsgfa => basis_set_a%nsgf_set
186 296129 : rpgfa => basis_set_a%pgf_radius
187 296129 : set_radius_a => basis_set_a%set_radius
188 296129 : sphi_a => basis_set_a%sphi
189 296129 : zeta => basis_set_a%zet
190 :
191 296129 : first_sgfb => basis_set_b%first_sgf
192 296129 : lb_max => basis_set_b%lmax
193 296129 : lb_min => basis_set_b%lmin
194 296129 : npgfb => basis_set_b%npgf
195 296129 : nsgfb => basis_set_b%nsgf_set
196 296129 : rpgfb => basis_set_b%pgf_radius
197 296129 : set_radius_b => basis_set_b%set_radius
198 296129 : sphi_b => basis_set_b%sphi
199 296129 : zetb => basis_set_b%zet
200 :
201 296129 : nseta = basis_set_a%nset
202 296129 : nsetb = basis_set_b%nset
203 296129 : ldsa = SIZE(sphi_a, 1)
204 296129 : ldsb = SIZE(sphi_b, 1)
205 296129 : IF (inode == 1) last_jatom = 0
206 :
207 296129 : IF (jatom /= last_jatom) THEN
208 : new_atom_b = .TRUE.
209 : last_jatom = jatom
210 : ELSE
211 : new_atom_b = .FALSE.
212 : END IF
213 :
214 : IF (new_atom_b) THEN
215 62032 : IF (iatom <= jatom) THEN
216 38092 : irow = iatom
217 38092 : icol = jatom
218 : ELSE
219 23940 : irow = jatom
220 23940 : icol = iatom
221 : END IF
222 :
223 162360 : DO i = 1, nkvec
224 300984 : DO reim = 1, 2
225 200656 : NULLIFY (op_cossin(reim, i)%block)
226 : CALL dbcsr_get_block_p(matrix=op_sm_set(reim, i)%matrix, &
227 300984 : row=irow, col=icol, block=op_cossin(reim, i)%block, found=found)
228 : END DO
229 162360 : IF (ASSOCIATED(op_cossin(1, i)%block) .NEQV. ASSOCIATED(op_cossin(2, i)%block)) THEN
230 0 : CPABORT("cosine and sine blocks should have the same topology")
231 : END IF
232 : END DO
233 : END IF
234 :
235 1184516 : dab = NORM2(rab)
236 1106157 : DO iset = 1, nseta
237 796056 : ncoa = npgfa(iset)*ncoset(la_max(iset))
238 796056 : sgfa = first_sgfa(1, iset)
239 4146786 : DO jset = 1, nsetb
240 3054601 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
241 1095147 : ncob = npgfb(jset)*ncoset(lb_max(jset))
242 1095147 : sgfb = first_sgfb(1, jset)
243 :
244 3297566 : DO i = 1, nkvec
245 1406363 : IF (.NOT. ASSOCIATED(op_cossin(1, i)%block)) CYCLE
246 : CALL expikr(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), &
247 : la_min(iset), lb_max(jset), npgfb(jset), zetb(:, jset), &
248 1406363 : rpgfb(:, jset), lb_min(jset), ra, rb, kvec(:, i), expab)
249 : CALL contract_expikr(op_cossin(1, i)%block, op_cossin(2, i)%block, &
250 : iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
251 : jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
252 4460964 : expab, work, ldwork)
253 : END DO
254 : END DO
255 : END DO
256 : END DO
257 13972 : CALL neighbor_list_iterator_release(nl_iterator)
258 :
259 13972 : DEALLOCATE (op_cossin, expab, work, basis_set_list)
260 :
261 15454 : IF (my_force_periodic) cell%perd(:) = perd0(:)
262 :
263 13972 : CALL timestop(handle)
264 13972 : END SUBROUTINE build_exp_ikr_matrix
265 :
266 : ! **************************************************************************************************
267 : !> \brief Calculation of the linear momentum matrix <mu|∂|nu> over
268 : !> Cartesian Gaussian functions.
269 : !> \param qs_env ...
270 : !> \param matrix ...
271 : !> \param minimum_image take into account only the first neighbors in the lists
272 : !> \date 27.02.2009
273 : !> \author VW
274 : !> \version 1.0
275 : ! **************************************************************************************************
276 1246 : SUBROUTINE build_lin_mom_matrix(qs_env, matrix, minimum_image)
277 :
278 : TYPE(qs_environment_type), POINTER :: qs_env
279 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
280 : LOGICAL, INTENT(IN), OPTIONAL :: minimum_image
281 :
282 : CHARACTER(len=*), PARAMETER :: routineN = 'build_lin_mom_matrix'
283 :
284 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
285 : ldai, maxco, maxlgto, maxsgf, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
286 1246 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nsgfa, &
287 1246 : nsgfb
288 1246 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
289 : LOGICAL :: do_symmetric, found, my_minimum_image, &
290 : new_atom_b
291 : REAL(KIND=dp) :: alpha, dab, Lxo2, Lyo2, Lzo2
292 1246 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: rr_work
293 1246 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
294 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
295 : TARGET :: intab
296 : REAL(KIND=dp), DIMENSION(3) :: rab
297 1246 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
298 1246 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
299 1246 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
300 : TYPE(cell_type), POINTER :: cell
301 1246 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
302 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
303 : TYPE(neighbor_list_iterator_p_type), &
304 1246 : DIMENSION(:), POINTER :: nl_iterator
305 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
306 1246 : POINTER :: sab_nl
307 1246 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
308 : TYPE(qs_kind_type), POINTER :: qs_kind
309 :
310 1246 : CALL timeset(routineN, handle)
311 :
312 1246 : NULLIFY (cell, sab_nl, qs_kind_set)
313 :
314 : CALL get_qs_env(qs_env=qs_env, &
315 : qs_kind_set=qs_kind_set, &
316 1246 : cell=cell)
317 :
318 1246 : nkind = SIZE(qs_kind_set)
319 :
320 1246 : my_minimum_image = .FALSE.
321 1246 : IF (PRESENT(minimum_image)) THEN
322 44 : my_minimum_image = minimum_image
323 176 : Lxo2 = NORM2(cell%hmat(:, 1))/2.0_dp
324 176 : Lyo2 = NORM2(cell%hmat(:, 2))/2.0_dp
325 176 : Lzo2 = NORM2(cell%hmat(:, 3))/2.0_dp
326 : END IF
327 :
328 : ! Take into account the symmetry of the input matrix
329 1246 : do_symmetric = dbcsr_has_symmetry(matrix(1)%matrix)
330 1246 : IF (do_symmetric) THEN
331 1244 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
332 : ELSE
333 2 : CALL get_qs_env(qs_env=qs_env, sab_all=sab_nl)
334 : END IF
335 : ! *** Allocate work storage ***
336 :
337 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
338 : maxco=maxco, &
339 : maxlgto=maxlgto, &
340 1246 : maxsgf=maxsgf)
341 :
342 1246 : ldai = ncoset(maxlgto + 1)
343 1246 : CALL init_orbital_pointers(ldai)
344 :
345 16198 : ALLOCATE (rr_work(ldai*ldai*3), intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
346 1246 : rr_work(:) = 0.0_dp
347 1246 : intab(:, :, :) = 0.0_dp
348 1246 : work(:, :) = 0.0_dp
349 :
350 5892 : ALLOCATE (basis_set_list(nkind))
351 3400 : DO ikind = 1, nkind
352 2154 : qs_kind => qs_kind_set(ikind)
353 2154 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
354 3400 : IF (ASSOCIATED(basis_set_a)) THEN
355 2154 : basis_set_list(ikind)%gto_basis_set => basis_set_a
356 : ELSE
357 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
358 : END IF
359 : END DO
360 1246 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
361 72446 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
362 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
363 71200 : iatom=iatom, jatom=jatom, r=rab)
364 71200 : basis_set_a => basis_set_list(ikind)%gto_basis_set
365 71200 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
366 71200 : basis_set_b => basis_set_list(jkind)%gto_basis_set
367 71200 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
368 : ! basis ikind
369 71200 : first_sgfa => basis_set_a%first_sgf
370 71200 : la_max => basis_set_a%lmax
371 71200 : npgfa => basis_set_a%npgf
372 71200 : nseta = basis_set_a%nset
373 71200 : nsgfa => basis_set_a%nsgf_set
374 71200 : rpgfa => basis_set_a%pgf_radius
375 71200 : set_radius_a => basis_set_a%set_radius
376 71200 : sphi_a => basis_set_a%sphi
377 71200 : zeta => basis_set_a%zet
378 : ! basis jkind
379 71200 : first_sgfb => basis_set_b%first_sgf
380 71200 : lb_max => basis_set_b%lmax
381 71200 : npgfb => basis_set_b%npgf
382 71200 : nsetb = basis_set_b%nset
383 71200 : nsgfb => basis_set_b%nsgf_set
384 71200 : rpgfb => basis_set_b%pgf_radius
385 71200 : set_radius_b => basis_set_b%set_radius
386 71200 : sphi_b => basis_set_b%sphi
387 71200 : zetb => basis_set_b%zet
388 :
389 71200 : IF (inode == 1) last_jatom = 0
390 :
391 71200 : IF (my_minimum_image) THEN
392 26557 : IF (ANY(ABS(rab(:)) > [Lxo2, Lyo2, Lzo2])) CYCLE
393 : END IF
394 :
395 67736 : IF (jatom /= last_jatom) THEN
396 : new_atom_b = .TRUE.
397 : last_jatom = jatom
398 : ELSE
399 : new_atom_b = .FALSE.
400 : END IF
401 :
402 : IF (new_atom_b) THEN
403 13819 : alpha = 1.0_dp
404 13819 : IF (do_symmetric) THEN
405 13810 : IF (iatom <= jatom) THEN
406 7583 : irow = iatom
407 7583 : icol = jatom
408 : ELSE
409 6227 : irow = jatom
410 6227 : icol = iatom
411 6227 : IF (dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_antisymmetric) THEN
412 12454 : alpha = -1.0_dp
413 : END IF
414 : END IF
415 : ELSE
416 9 : irow = iatom
417 9 : icol = jatom
418 : END IF
419 :
420 55276 : DO i = 1, 3
421 41457 : NULLIFY (integral(i)%block)
422 : CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
423 41457 : row=irow, col=icol, BLOCK=integral(i)%block, found=found)
424 55276 : CPASSERT(ASSOCIATED(INTEGRAL(i)%block))
425 : END DO
426 : END IF
427 :
428 270944 : dab = NORM2(rab)
429 :
430 213078 : DO iset = 1, nseta
431 :
432 144096 : ncoa = npgfa(iset)*ncoset(la_max(iset))
433 144096 : sgfa = first_sgfa(1, iset)
434 :
435 547953 : DO jset = 1, nsetb
436 :
437 336121 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
438 :
439 144107 : ncob = npgfb(jset)*ncoset(lb_max(jset))
440 144107 : sgfb = first_sgfb(1, jset)
441 :
442 : ! *** Calculate the primitive fermi contact integrals ***
443 :
444 : ! Keep the full primitive layout used by the contraction below.
445 : CALL overlap_ab(la_max(iset), 0, npgfa(iset), &
446 : rpgfa(:, iset), zeta(:, iset), &
447 : lb_max(jset), 0, npgfb(jset), &
448 : rpgfb(:, jset), zetb(:, jset), &
449 144107 : rab, dab=intab, rr_work=rr_work)
450 :
451 : ! *** Contraction step ***
452 :
453 720524 : DO i = 1, 3
454 :
455 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
456 : 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
457 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
458 432321 : 0.0_dp, work(1, 1), SIZE(work, 1))
459 :
460 768442 : IF (do_symmetric) THEN
461 432294 : IF (iatom <= jatom) THEN
462 :
463 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
464 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
465 : work(1, 1), SIZE(work, 1), &
466 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
467 266226 : SIZE(integral(i)%block, 1))
468 :
469 : ELSE
470 :
471 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
472 : alpha, work(1, 1), SIZE(work, 1), &
473 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
474 : 1.0_dp, integral(i)%block(sgfb, sgfa), &
475 166068 : SIZE(integral(i)%block, 1))
476 :
477 : END IF
478 : ELSE
479 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
480 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
481 : work(1, 1), SIZE(work, 1), &
482 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
483 27 : SIZE(integral(i)%block, 1))
484 : END IF
485 :
486 : END DO
487 :
488 : END DO
489 :
490 : END DO
491 :
492 : END DO
493 1246 : CALL neighbor_list_iterator_release(nl_iterator)
494 :
495 : ! *** Release work storage ***
496 :
497 1246 : DEALLOCATE (intab, rr_work, work, integral, basis_set_list)
498 :
499 1246 : CALL timestop(handle)
500 :
501 3738 : END SUBROUTINE build_lin_mom_matrix
502 :
503 : ! **************************************************************************************************
504 : !> \brief Calculation of the angular momentum matrix over
505 : !> Cartesian Gaussian functions.
506 : !> \param qs_env ...
507 : !> \param matrix ...
508 : !> \param rc ...
509 : !> \date 27.02.2009
510 : !> \author VW
511 : !> \version 1.0
512 : ! **************************************************************************************************
513 :
514 1194 : SUBROUTINE build_ang_mom_matrix(qs_env, matrix, rc)
515 :
516 : TYPE(qs_environment_type), POINTER :: qs_env
517 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
518 : REAL(dp), DIMENSION(:), INTENT(IN) :: rc
519 :
520 : CHARACTER(len=*), PARAMETER :: routineN = 'build_ang_mom_matrix'
521 :
522 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
523 : maxco, maxsgf, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
524 1194 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
525 1194 : npgfb, nsgfa, nsgfb
526 1194 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
527 : LOGICAL :: found, new_atom_b
528 : REAL(KIND=dp) :: dab
529 1194 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
530 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab
531 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rbc
532 1194 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
533 1194 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
534 1194 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
535 : TYPE(cell_type), POINTER :: cell
536 1194 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
537 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
538 : TYPE(neighbor_list_iterator_p_type), &
539 1194 : DIMENSION(:), POINTER :: nl_iterator
540 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
541 1194 : POINTER :: sab_all
542 1194 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
543 1194 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
544 : TYPE(qs_kind_type), POINTER :: qs_kind
545 :
546 1194 : CALL timeset(routineN, handle)
547 :
548 : CALL get_qs_env(qs_env=qs_env, &
549 : qs_kind_set=qs_kind_set, &
550 : particle_set=particle_set, &
551 : sab_all=sab_all, &
552 1194 : cell=cell)
553 :
554 1194 : nkind = SIZE(qs_kind_set)
555 :
556 : ! *** Allocate work storage ***
557 :
558 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
559 : maxco=maxco, &
560 1194 : maxsgf=maxsgf)
561 :
562 13134 : ALLOCATE (intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
563 1194 : intab(:, :, :) = 0.0_dp
564 1194 : work(:, :) = 0.0_dp
565 :
566 5540 : ALLOCATE (basis_set_list(nkind))
567 3152 : DO ikind = 1, nkind
568 1958 : qs_kind => qs_kind_set(ikind)
569 1958 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
570 3152 : IF (ASSOCIATED(basis_set_a)) THEN
571 1958 : basis_set_list(ikind)%gto_basis_set => basis_set_a
572 : ELSE
573 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
574 : END IF
575 : END DO
576 1194 : CALL neighbor_list_iterator_create(nl_iterator, sab_all)
577 90709 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
578 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
579 89515 : iatom=iatom, jatom=jatom, r=rab)
580 89515 : basis_set_a => basis_set_list(ikind)%gto_basis_set
581 89515 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
582 89515 : basis_set_b => basis_set_list(jkind)%gto_basis_set
583 89515 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
584 89515 : ra = pbc(particle_set(iatom)%r, cell)
585 : ! basis ikind
586 89515 : first_sgfa => basis_set_a%first_sgf
587 89515 : la_max => basis_set_a%lmax
588 89515 : la_min => basis_set_a%lmin
589 89515 : npgfa => basis_set_a%npgf
590 89515 : nseta = basis_set_a%nset
591 89515 : nsgfa => basis_set_a%nsgf_set
592 89515 : rpgfa => basis_set_a%pgf_radius
593 89515 : set_radius_a => basis_set_a%set_radius
594 89515 : sphi_a => basis_set_a%sphi
595 89515 : zeta => basis_set_a%zet
596 : ! basis jkind
597 89515 : first_sgfb => basis_set_b%first_sgf
598 89515 : lb_max => basis_set_b%lmax
599 89515 : lb_min => basis_set_b%lmin
600 89515 : npgfb => basis_set_b%npgf
601 89515 : nsetb = basis_set_b%nset
602 89515 : nsgfb => basis_set_b%nsgf_set
603 89515 : rpgfb => basis_set_b%pgf_radius
604 89515 : set_radius_b => basis_set_b%set_radius
605 89515 : sphi_b => basis_set_b%sphi
606 89515 : zetb => basis_set_b%zet
607 :
608 89515 : IF (inode == 1) last_jatom = 0
609 :
610 89515 : IF (jatom /= last_jatom) THEN
611 : new_atom_b = .TRUE.
612 : last_jatom = jatom
613 : ELSE
614 : new_atom_b = .FALSE.
615 : END IF
616 :
617 : IF (new_atom_b) THEN
618 : !IF (iatom <= jatom) THEN
619 5795 : irow = iatom
620 5795 : icol = jatom
621 : !ELSE
622 : ! irow = jatom
623 : ! icol = iatom
624 : !END IF
625 :
626 23180 : DO i = 1, 3
627 17385 : NULLIFY (INTEGRAL(i)%block)
628 : CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
629 17385 : row=irow, col=icol, BLOCK=INTEGRAL(i)%block, found=found)
630 23180 : CPASSERT(found .AND. ASSOCIATED(INTEGRAL(i)%block))
631 : END DO
632 : END IF
633 :
634 358060 : dab = NORM2(rab)
635 :
636 299354 : DO iset = 1, nseta
637 :
638 208645 : ncoa = npgfa(iset)*ncoset(la_max(iset))
639 208645 : sgfa = first_sgfa(1, iset)
640 :
641 807795 : DO jset = 1, nsetb
642 :
643 509635 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
644 :
645 : !IF(PRESENT(wancen)) THEN
646 : ! rc = wancen
647 193525 : rac = pbc(rc, ra, cell)
648 774100 : rbc = rac + rab
649 : !ELSE
650 : ! rc(1:3) = rb(1:3)
651 : ! rac(1:3) = -rab(1:3)
652 : ! rbc(1:3) = 0.0_dp
653 : !ENDIF
654 :
655 193525 : ncob = npgfb(jset)*ncoset(lb_max(jset))
656 193525 : sgfb = first_sgfb(1, jset)
657 :
658 : ! *** Calculate the primitive angular momentum integrals ***
659 :
660 : CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
661 193525 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), rac, rbc, intab)
662 :
663 : ! *** Contraction step ***
664 :
665 982745 : DO i = 1, 3
666 :
667 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
668 : 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
669 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
670 580575 : 0.0_dp, work(1, 1), SIZE(work, 1))
671 :
672 : !IF (iatom <= jatom) THEN
673 :
674 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
675 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
676 : work(1, 1), SIZE(work, 1), &
677 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
678 1090210 : SIZE(integral(i)%block, 1))
679 :
680 : !ELSE
681 : !
682 : ! CALL dgemm("T","N",nsgfb(jset),nsgfa(iset),ncoa,&
683 : ! -1.0_dp,work(1,1),SIZE(work,1),&
684 : ! sphi_a(1,sgfa),SIZE(sphi_a,1),&
685 : ! 1.0_dp,integral(i)%block(sgfb,sgfa),&
686 : ! SIZE(integral(i)%block,1))
687 : !
688 : !ENDIF
689 :
690 : END DO
691 :
692 : END DO
693 :
694 : END DO
695 :
696 : END DO
697 1194 : CALL neighbor_list_iterator_release(nl_iterator)
698 :
699 : ! *** Release work storage ***
700 :
701 1194 : DEALLOCATE (intab, work, integral, basis_set_list)
702 :
703 : ! *** Print the spin orbit matrix, if requested ***
704 :
705 : !IF (BTEST(cp_print_key_should_output(logger%iter_info,&
706 : ! qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM"),cp_p_file)) THEN
707 : ! iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM",&
708 : ! extension=".Log")
709 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(1)%matrix,4,6,qs_env,para_env,output_unit=iw)
710 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(2)%matrix,4,6,qs_env,para_env,output_unit=iw)
711 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(3)%matrix,4,6,qs_env,para_env,output_unit=iw)
712 : ! CALL cp_print_key_finished_output(iw,logger,qs_env%input,&
713 : ! "DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM")
714 : !END IF
715 :
716 1194 : CALL timestop(handle)
717 :
718 2388 : END SUBROUTINE build_ang_mom_matrix
719 :
720 : ! **************************************************************************************************
721 : END MODULE qs_operators_ao
|
