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 : 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 13960 : 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 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
87 : ldsa, ldsb, ldwork, ncoa, ncob, nkind, nkvec, nseta, nsetb, reim, sgfa, sgfb
88 : INTEGER, DIMENSION(3) :: perd0
89 13960 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
90 13960 : npgfb, nsgfa, nsgfb
91 13960 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
92 : LOGICAL :: found, my_force_periodic, new_atom_b
93 : REAL(KIND=dp) :: dab
94 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb
95 13960 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
96 13960 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cosab, rpgfa, rpgfb, sinab, sphi_a, &
97 13960 : sphi_b, work, zeta, zetb
98 13960 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: op_cossin
99 : TYPE(cell_type), POINTER :: cell
100 13960 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
101 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
102 : TYPE(neighbor_list_iterator_p_type), &
103 13960 : DIMENSION(:), POINTER :: nl_iterator
104 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
105 13960 : POINTER :: sab_orb
106 13960 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
107 13960 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
108 : TYPE(qs_kind_type), POINTER :: qs_kind
109 :
110 13960 : CALL timeset(routineN, handle)
111 :
112 13960 : CPASSERT(SIZE(kvec, 1) == 3)
113 13960 : nkvec = SIZE(kvec, 2)
114 13960 : CPASSERT(nkvec > 0)
115 13960 : CPASSERT(SIZE(op_sm_set, 1) == 2)
116 13960 : CPASSERT(SIZE(op_sm_set, 2) == nkvec)
117 28964 : DO i = 1, nkvec
118 58972 : DO reim = 1, 2
119 30008 : CPASSERT(ASSOCIATED(op_sm_set(reim, i)%matrix))
120 30008 : CPASSERT(dbcsr_has_symmetry(op_sm_set(reim, i)%matrix))
121 45012 : CALL dbcsr_set(op_sm_set(reim, i)%matrix, 0.0_dp)
122 : END DO
123 : END DO
124 :
125 13960 : NULLIFY (qs_kind, qs_kind_set, particle_set, sab_orb, cell)
126 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
127 13960 : particle_set=particle_set, cell=cell, sab_orb=sab_orb)
128 13960 : CPASSERT(ASSOCIATED(cell))
129 13960 : CPASSERT(ASSOCIATED(qs_kind_set))
130 13960 : CPASSERT(ASSOCIATED(particle_set))
131 13960 : IF (PRESENT(cell_external)) THEN
132 504 : CPASSERT(ASSOCIATED(cell_external))
133 504 : cell => cell_external
134 : END IF
135 :
136 13960 : IF (PRESENT(sab_orb_external)) THEN
137 11160 : CPASSERT(ASSOCIATED(sab_orb_external))
138 11160 : sab_orb => sab_orb_external
139 : END IF
140 13960 : CPASSERT(ASSOCIATED(sab_orb))
141 :
142 13960 : my_force_periodic = .FALSE.
143 13960 : IF (PRESENT(force_periodic)) my_force_periodic = force_periodic
144 504 : IF (my_force_periodic) THEN
145 2016 : perd0(:) = cell%perd(:)
146 2016 : cell%perd(:) = 1
147 : END IF
148 :
149 13960 : nkind = SIZE(qs_kind_set)
150 16760 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork, basis_type=basis_type)
151 111680 : ALLOCATE (cosab(ldwork, ldwork), sinab(ldwork, ldwork), work(ldwork, ldwork))
152 :
153 86892 : ALLOCATE (op_cossin(2, nkvec))
154 :
155 69662 : ALLOCATE (basis_set_list(nkind))
156 41742 : DO ikind = 1, nkind
157 27782 : qs_kind => qs_kind_set(ikind)
158 27782 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
159 41742 : IF (ASSOCIATED(basis_set_a)) THEN
160 27782 : basis_set_list(ikind)%gto_basis_set => basis_set_a
161 : ELSE
162 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
163 : END IF
164 : END DO
165 :
166 13960 : last_jatom = 0
167 13960 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
168 309559 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
169 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
170 295599 : iatom=iatom, jatom=jatom, r=rab)
171 295599 : basis_set_a => basis_set_list(ikind)%gto_basis_set
172 295599 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
173 295599 : basis_set_b => basis_set_list(jkind)%gto_basis_set
174 295599 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
175 :
176 295599 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
177 1182396 : rb(:) = ra(:) + rab(:)
178 :
179 295599 : first_sgfa => basis_set_a%first_sgf
180 295599 : la_max => basis_set_a%lmax
181 295599 : la_min => basis_set_a%lmin
182 295599 : npgfa => basis_set_a%npgf
183 295599 : nsgfa => basis_set_a%nsgf_set
184 295599 : rpgfa => basis_set_a%pgf_radius
185 295599 : set_radius_a => basis_set_a%set_radius
186 295599 : sphi_a => basis_set_a%sphi
187 295599 : zeta => basis_set_a%zet
188 :
189 295599 : first_sgfb => basis_set_b%first_sgf
190 295599 : lb_max => basis_set_b%lmax
191 295599 : lb_min => basis_set_b%lmin
192 295599 : npgfb => basis_set_b%npgf
193 295599 : nsgfb => basis_set_b%nsgf_set
194 295599 : rpgfb => basis_set_b%pgf_radius
195 295599 : set_radius_b => basis_set_b%set_radius
196 295599 : sphi_b => basis_set_b%sphi
197 295599 : zetb => basis_set_b%zet
198 :
199 295599 : nseta = basis_set_a%nset
200 295599 : nsetb = basis_set_b%nset
201 295599 : ldsa = SIZE(sphi_a, 1)
202 295599 : ldsb = SIZE(sphi_b, 1)
203 295599 : IF (inode == 1) last_jatom = 0
204 :
205 295599 : IF (jatom /= last_jatom) THEN
206 : new_atom_b = .TRUE.
207 : last_jatom = jatom
208 : ELSE
209 : new_atom_b = .FALSE.
210 : END IF
211 :
212 : IF (new_atom_b) THEN
213 62003 : IF (iatom <= jatom) THEN
214 38072 : irow = iatom
215 38072 : icol = jatom
216 : ELSE
217 23931 : irow = jatom
218 23931 : icol = iatom
219 : END IF
220 :
221 162356 : DO i = 1, nkvec
222 301059 : DO reim = 1, 2
223 200706 : NULLIFY (op_cossin(reim, i)%block)
224 : CALL dbcsr_get_block_p(matrix=op_sm_set(reim, i)%matrix, &
225 301059 : row=irow, col=icol, block=op_cossin(reim, i)%block, found=found)
226 : END DO
227 162356 : IF (ASSOCIATED(op_cossin(1, i)%block) .NEQV. ASSOCIATED(op_cossin(2, i)%block)) THEN
228 0 : CPABORT("cosine and sine blocks should have the same topology")
229 : END IF
230 : END DO
231 : END IF
232 :
233 1182396 : dab = NORM2(rab)
234 1104555 : DO iset = 1, nseta
235 794996 : ncoa = npgfa(iset)*ncoset(la_max(iset))
236 794996 : sgfa = first_sgfa(1, iset)
237 4143076 : DO jset = 1, nsetb
238 3052481 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
239 1094156 : ncob = npgfb(jset)*ncoset(lb_max(jset))
240 1094156 : sgfb = first_sgfb(1, jset)
241 :
242 3296636 : DO i = 1, nkvec
243 1407484 : IF (.NOT. ASSOCIATED(op_cossin(1, i)%block)) CYCLE
244 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), &
245 : la_min(iset), lb_max(jset), npgfb(jset), zetb(:, jset), &
246 1407484 : rpgfb(:, jset), lb_min(jset), ra, rb, kvec(:, i), cosab, sinab)
247 : CALL contract_cossin(op_cossin(1, i)%block, op_cossin(2, i)%block, &
248 : iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
249 : jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
250 4459965 : cosab, sinab, ldwork, work, ldwork)
251 : END DO
252 : END DO
253 : END DO
254 : END DO
255 13960 : CALL neighbor_list_iterator_release(nl_iterator)
256 :
257 13960 : DEALLOCATE (op_cossin, cosab, sinab, work, basis_set_list)
258 :
259 15472 : IF (my_force_periodic) cell%perd(:) = perd0(:)
260 :
261 13960 : CALL timestop(handle)
262 13960 : END SUBROUTINE build_exp_ikr_matrix
263 :
264 : ! **************************************************************************************************
265 : !> \brief Calculation of the linear momentum matrix <mu|∂|nu> over
266 : !> Cartesian Gaussian functions.
267 : !> \param qs_env ...
268 : !> \param matrix ...
269 : !> \param minimum_image take into account only the first neighbors in the lists
270 : !> \date 27.02.2009
271 : !> \author VW
272 : !> \version 1.0
273 : ! **************************************************************************************************
274 1296 : SUBROUTINE build_lin_mom_matrix(qs_env, matrix, minimum_image)
275 :
276 : TYPE(qs_environment_type), POINTER :: qs_env
277 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
278 : LOGICAL, INTENT(IN), OPTIONAL :: minimum_image
279 :
280 : CHARACTER(len=*), PARAMETER :: routineN = 'build_lin_mom_matrix'
281 :
282 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
283 : ldai, maxco, maxlgto, maxsgf, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
284 1296 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nsgfa, &
285 1296 : nsgfb
286 1296 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
287 : LOGICAL :: do_symmetric, found, my_minimum_image, &
288 : new_atom_b
289 : REAL(KIND=dp) :: alpha, dab, Lxo2, Lyo2, Lzo2
290 1296 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: rr_work
291 1296 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
292 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
293 : TARGET :: intab
294 : REAL(KIND=dp), DIMENSION(3) :: rab
295 1296 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
296 1296 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
297 1296 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
298 : TYPE(cell_type), POINTER :: cell
299 1296 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
300 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
301 : TYPE(neighbor_list_iterator_p_type), &
302 1296 : DIMENSION(:), POINTER :: nl_iterator
303 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
304 1296 : POINTER :: sab_nl
305 1296 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
306 : TYPE(qs_kind_type), POINTER :: qs_kind
307 :
308 1296 : CALL timeset(routineN, handle)
309 :
310 1296 : NULLIFY (cell, sab_nl, qs_kind_set)
311 :
312 : CALL get_qs_env(qs_env=qs_env, &
313 : qs_kind_set=qs_kind_set, &
314 1296 : cell=cell)
315 :
316 1296 : nkind = SIZE(qs_kind_set)
317 :
318 1296 : my_minimum_image = .FALSE.
319 1296 : IF (PRESENT(minimum_image)) THEN
320 44 : my_minimum_image = minimum_image
321 176 : Lxo2 = NORM2(cell%hmat(:, 1))/2.0_dp
322 176 : Lyo2 = NORM2(cell%hmat(:, 2))/2.0_dp
323 176 : Lzo2 = NORM2(cell%hmat(:, 3))/2.0_dp
324 : END IF
325 :
326 : ! Take into account the symmetry of the input matrix
327 1296 : do_symmetric = dbcsr_has_symmetry(matrix(1)%matrix)
328 1296 : IF (do_symmetric) THEN
329 1294 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
330 : ELSE
331 2 : CALL get_qs_env(qs_env=qs_env, sab_all=sab_nl)
332 : END IF
333 : ! *** Allocate work storage ***
334 :
335 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
336 : maxco=maxco, &
337 : maxlgto=maxlgto, &
338 1296 : maxsgf=maxsgf)
339 :
340 1296 : ldai = ncoset(maxlgto + 1)
341 1296 : CALL init_orbital_pointers(ldai)
342 :
343 16848 : ALLOCATE (rr_work(ldai*ldai*3), intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
344 1296 : rr_work(:) = 0.0_dp
345 1296 : intab(:, :, :) = 0.0_dp
346 1296 : work(:, :) = 0.0_dp
347 :
348 6132 : ALLOCATE (basis_set_list(nkind))
349 3540 : DO ikind = 1, nkind
350 2244 : qs_kind => qs_kind_set(ikind)
351 2244 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
352 3540 : IF (ASSOCIATED(basis_set_a)) THEN
353 2244 : basis_set_list(ikind)%gto_basis_set => basis_set_a
354 : ELSE
355 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
356 : END IF
357 : END DO
358 1296 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
359 75246 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
360 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
361 73950 : iatom=iatom, jatom=jatom, r=rab)
362 73950 : basis_set_a => basis_set_list(ikind)%gto_basis_set
363 73950 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
364 73950 : basis_set_b => basis_set_list(jkind)%gto_basis_set
365 73950 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
366 : ! basis ikind
367 73950 : first_sgfa => basis_set_a%first_sgf
368 73950 : la_max => basis_set_a%lmax
369 73950 : npgfa => basis_set_a%npgf
370 73950 : nseta = basis_set_a%nset
371 73950 : nsgfa => basis_set_a%nsgf_set
372 73950 : rpgfa => basis_set_a%pgf_radius
373 73950 : set_radius_a => basis_set_a%set_radius
374 73950 : sphi_a => basis_set_a%sphi
375 73950 : zeta => basis_set_a%zet
376 : ! basis jkind
377 73950 : first_sgfb => basis_set_b%first_sgf
378 73950 : lb_max => basis_set_b%lmax
379 73950 : npgfb => basis_set_b%npgf
380 73950 : nsetb = basis_set_b%nset
381 73950 : nsgfb => basis_set_b%nsgf_set
382 73950 : rpgfb => basis_set_b%pgf_radius
383 73950 : set_radius_b => basis_set_b%set_radius
384 73950 : sphi_b => basis_set_b%sphi
385 73950 : zetb => basis_set_b%zet
386 :
387 73950 : IF (inode == 1) last_jatom = 0
388 :
389 73950 : IF (my_minimum_image) THEN
390 26557 : IF (ANY(ABS(rab(:)) > [Lxo2, Lyo2, Lzo2])) CYCLE
391 : END IF
392 :
393 70486 : IF (jatom /= last_jatom) THEN
394 : new_atom_b = .TRUE.
395 : last_jatom = jatom
396 : ELSE
397 : new_atom_b = .FALSE.
398 : END IF
399 :
400 : IF (new_atom_b) THEN
401 13954 : alpha = 1.0_dp
402 13954 : IF (do_symmetric) THEN
403 13945 : IF (iatom <= jatom) THEN
404 7673 : irow = iatom
405 7673 : icol = jatom
406 : ELSE
407 6272 : irow = jatom
408 6272 : icol = iatom
409 6272 : IF (dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_antisymmetric) THEN
410 12544 : alpha = -1.0_dp
411 : END IF
412 : END IF
413 : ELSE
414 9 : irow = iatom
415 9 : icol = jatom
416 : END IF
417 :
418 55816 : DO i = 1, 3
419 41862 : NULLIFY (integral(i)%block)
420 : CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
421 41862 : row=irow, col=icol, BLOCK=integral(i)%block, found=found)
422 55816 : CPASSERT(ASSOCIATED(INTEGRAL(i)%block))
423 : END DO
424 : END IF
425 :
426 281944 : dab = NORM2(rab)
427 :
428 221378 : DO iset = 1, nseta
429 :
430 149596 : ncoa = npgfa(iset)*ncoset(la_max(iset))
431 149596 : sgfa = first_sgfa(1, iset)
432 :
433 567203 : DO jset = 1, nsetb
434 :
435 347121 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
436 :
437 149387 : ncob = npgfb(jset)*ncoset(lb_max(jset))
438 149387 : sgfb = first_sgfb(1, jset)
439 :
440 : ! *** Calculate the primitive fermi contact integrals ***
441 :
442 : ! Keep the full primitive layout used by the contraction below.
443 : CALL overlap_ab(la_max(iset), 0, npgfa(iset), &
444 : rpgfa(:, iset), zeta(:, iset), &
445 : lb_max(jset), 0, npgfb(jset), &
446 : rpgfb(:, jset), zetb(:, jset), &
447 149387 : rab, dab=intab, rr_work=rr_work)
448 :
449 : ! *** Contraction step ***
450 :
451 747144 : DO i = 1, 3
452 :
453 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
454 : 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
455 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
456 448161 : 0.0_dp, work(1, 1), SIZE(work, 1))
457 :
458 795282 : IF (do_symmetric) THEN
459 448134 : IF (iatom <= jatom) THEN
460 :
461 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
462 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
463 : work(1, 1), SIZE(work, 1), &
464 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
465 276486 : SIZE(integral(i)%block, 1))
466 :
467 : ELSE
468 :
469 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
470 : alpha, work(1, 1), SIZE(work, 1), &
471 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
472 : 1.0_dp, integral(i)%block(sgfb, sgfa), &
473 171648 : SIZE(integral(i)%block, 1))
474 :
475 : END IF
476 : ELSE
477 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
478 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
479 : work(1, 1), SIZE(work, 1), &
480 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
481 27 : SIZE(integral(i)%block, 1))
482 : END IF
483 :
484 : END DO
485 :
486 : END DO
487 :
488 : END DO
489 :
490 : END DO
491 1296 : CALL neighbor_list_iterator_release(nl_iterator)
492 :
493 : ! *** Release work storage ***
494 :
495 1296 : DEALLOCATE (intab, rr_work, work, integral, basis_set_list)
496 :
497 1296 : CALL timestop(handle)
498 :
499 3888 : END SUBROUTINE build_lin_mom_matrix
500 :
501 : ! **************************************************************************************************
502 : !> \brief Calculation of the angular momentum matrix over
503 : !> Cartesian Gaussian functions.
504 : !> \param qs_env ...
505 : !> \param matrix ...
506 : !> \param rc ...
507 : !> \date 27.02.2009
508 : !> \author VW
509 : !> \version 1.0
510 : ! **************************************************************************************************
511 :
512 1250 : SUBROUTINE build_ang_mom_matrix(qs_env, matrix, rc)
513 :
514 : TYPE(qs_environment_type), POINTER :: qs_env
515 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
516 : REAL(dp), DIMENSION(:), INTENT(IN) :: rc
517 :
518 : CHARACTER(len=*), PARAMETER :: routineN = 'build_ang_mom_matrix'
519 :
520 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
521 : maxco, maxsgf, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
522 1250 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
523 1250 : npgfb, nsgfa, nsgfb
524 1250 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
525 : LOGICAL :: found, new_atom_b
526 : REAL(KIND=dp) :: dab
527 1250 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
528 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab
529 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rbc
530 1250 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
531 1250 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
532 1250 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
533 : TYPE(cell_type), POINTER :: cell
534 1250 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
535 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
536 : TYPE(neighbor_list_iterator_p_type), &
537 1250 : DIMENSION(:), POINTER :: nl_iterator
538 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
539 1250 : POINTER :: sab_all
540 1250 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
541 1250 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
542 : TYPE(qs_kind_type), POINTER :: qs_kind
543 :
544 1250 : CALL timeset(routineN, handle)
545 :
546 : CALL get_qs_env(qs_env=qs_env, &
547 : qs_kind_set=qs_kind_set, &
548 : particle_set=particle_set, &
549 : sab_all=sab_all, &
550 1250 : cell=cell)
551 :
552 1250 : nkind = SIZE(qs_kind_set)
553 :
554 : ! *** Allocate work storage ***
555 :
556 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
557 : maxco=maxco, &
558 1250 : maxsgf=maxsgf)
559 :
560 13750 : ALLOCATE (intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
561 1250 : intab(:, :, :) = 0.0_dp
562 1250 : work(:, :) = 0.0_dp
563 :
564 5796 : ALLOCATE (basis_set_list(nkind))
565 3296 : DO ikind = 1, nkind
566 2046 : qs_kind => qs_kind_set(ikind)
567 2046 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
568 3296 : IF (ASSOCIATED(basis_set_a)) THEN
569 2046 : basis_set_list(ikind)%gto_basis_set => basis_set_a
570 : ELSE
571 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
572 : END IF
573 : END DO
574 1250 : CALL neighbor_list_iterator_create(nl_iterator, sab_all)
575 94981 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
576 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
577 93731 : iatom=iatom, jatom=jatom, r=rab)
578 93731 : basis_set_a => basis_set_list(ikind)%gto_basis_set
579 93731 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
580 93731 : basis_set_b => basis_set_list(jkind)%gto_basis_set
581 93731 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
582 93731 : ra = pbc(particle_set(iatom)%r, cell)
583 : ! basis ikind
584 93731 : first_sgfa => basis_set_a%first_sgf
585 93731 : la_max => basis_set_a%lmax
586 93731 : la_min => basis_set_a%lmin
587 93731 : npgfa => basis_set_a%npgf
588 93731 : nseta = basis_set_a%nset
589 93731 : nsgfa => basis_set_a%nsgf_set
590 93731 : rpgfa => basis_set_a%pgf_radius
591 93731 : set_radius_a => basis_set_a%set_radius
592 93731 : sphi_a => basis_set_a%sphi
593 93731 : zeta => basis_set_a%zet
594 : ! basis jkind
595 93731 : first_sgfb => basis_set_b%first_sgf
596 93731 : lb_max => basis_set_b%lmax
597 93731 : lb_min => basis_set_b%lmin
598 93731 : npgfb => basis_set_b%npgf
599 93731 : nsetb = basis_set_b%nset
600 93731 : nsgfb => basis_set_b%nsgf_set
601 93731 : rpgfb => basis_set_b%pgf_radius
602 93731 : set_radius_b => basis_set_b%set_radius
603 93731 : sphi_b => basis_set_b%sphi
604 93731 : zetb => basis_set_b%zet
605 :
606 93731 : IF (inode == 1) last_jatom = 0
607 :
608 93731 : IF (jatom /= last_jatom) THEN
609 : new_atom_b = .TRUE.
610 : last_jatom = jatom
611 : ELSE
612 : new_atom_b = .FALSE.
613 : END IF
614 :
615 : IF (new_atom_b) THEN
616 : !IF (iatom <= jatom) THEN
617 5987 : irow = iatom
618 5987 : icol = jatom
619 : !ELSE
620 : ! irow = jatom
621 : ! icol = iatom
622 : !END IF
623 :
624 23948 : DO i = 1, 3
625 17961 : NULLIFY (INTEGRAL(i)%block)
626 : CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
627 17961 : row=irow, col=icol, BLOCK=INTEGRAL(i)%block, found=found)
628 23948 : CPASSERT(found .AND. ASSOCIATED(INTEGRAL(i)%block))
629 : END DO
630 : END IF
631 :
632 374924 : dab = NORM2(rab)
633 :
634 312058 : DO iset = 1, nseta
635 :
636 217077 : ncoa = npgfa(iset)*ncoset(la_max(iset))
637 217077 : sgfa = first_sgfa(1, iset)
638 :
639 837307 : DO jset = 1, nsetb
640 :
641 526499 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
642 :
643 : !IF(PRESENT(wancen)) THEN
644 : ! rc = wancen
645 201893 : rac = pbc(rc, ra, cell)
646 807572 : rbc = rac + rab
647 : !ELSE
648 : ! rc(1:3) = rb(1:3)
649 : ! rac(1:3) = -rab(1:3)
650 : ! rbc(1:3) = 0.0_dp
651 : !ENDIF
652 :
653 201893 : ncob = npgfb(jset)*ncoset(lb_max(jset))
654 201893 : sgfb = first_sgfb(1, jset)
655 :
656 : ! *** Calculate the primitive angular momentum integrals ***
657 :
658 : CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
659 201893 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), rac, rbc, intab)
660 :
661 : ! *** Contraction step ***
662 :
663 1024649 : DO i = 1, 3
664 :
665 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
666 : 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
667 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
668 605679 : 0.0_dp, work(1, 1), SIZE(work, 1))
669 :
670 : !IF (iatom <= jatom) THEN
671 :
672 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
673 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
674 : work(1, 1), SIZE(work, 1), &
675 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
676 1132178 : SIZE(integral(i)%block, 1))
677 :
678 : !ELSE
679 : !
680 : ! CALL dgemm("T","N",nsgfb(jset),nsgfa(iset),ncoa,&
681 : ! -1.0_dp,work(1,1),SIZE(work,1),&
682 : ! sphi_a(1,sgfa),SIZE(sphi_a,1),&
683 : ! 1.0_dp,integral(i)%block(sgfb,sgfa),&
684 : ! SIZE(integral(i)%block,1))
685 : !
686 : !ENDIF
687 :
688 : END DO
689 :
690 : END DO
691 :
692 : END DO
693 :
694 : END DO
695 1250 : CALL neighbor_list_iterator_release(nl_iterator)
696 :
697 : ! *** Release work storage ***
698 :
699 1250 : DEALLOCATE (intab, work, integral, basis_set_list)
700 :
701 : ! *** Print the spin orbit matrix, if requested ***
702 :
703 : !IF (BTEST(cp_print_key_should_output(logger%iter_info,&
704 : ! qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM"),cp_p_file)) THEN
705 : ! iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM",&
706 : ! extension=".Log")
707 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(1)%matrix,4,6,qs_env,para_env,output_unit=iw)
708 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(2)%matrix,4,6,qs_env,para_env,output_unit=iw)
709 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(3)%matrix,4,6,qs_env,para_env,output_unit=iw)
710 : ! CALL cp_print_key_finished_output(iw,logger,qs_env%input,&
711 : ! "DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM")
712 : !END IF
713 :
714 1250 : CALL timestop(handle)
715 :
716 2500 : END SUBROUTINE build_ang_mom_matrix
717 :
718 : ! **************************************************************************************************
719 : END MODULE qs_operators_ao
|
