Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \par History
10 : !> created 07.2005
11 : !> \author MI (07.2005)
12 : ! **************************************************************************************************
13 : MODULE qs_operators_ao
14 : USE ai_moments, ONLY: diff_momop,&
15 : diffop,&
16 : moment
17 : USE ai_os_rr, ONLY: os_rr_ovlp
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_type_antisymmetric,&
28 : dbcsr_type_no_symmetry
29 : USE cp_log_handling, ONLY: cp_get_default_logger,&
30 : cp_logger_type
31 : USE kinds, ONLY: default_string_length,&
32 : dp
33 : USE mathconstants, ONLY: pi
34 : USE message_passing, ONLY: mp_para_env_type
35 : USE orbital_pointers, ONLY: coset,&
36 : init_orbital_pointers,&
37 : ncoset
38 : USE particle_types, ONLY: particle_type
39 : USE qs_environment_types, ONLY: get_qs_env,&
40 : qs_environment_type
41 : USE qs_kind_types, ONLY: get_qs_kind,&
42 : get_qs_kind_set,&
43 : qs_kind_type
44 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
45 : neighbor_list_iterate,&
46 : neighbor_list_iterator_create,&
47 : neighbor_list_iterator_p_type,&
48 : neighbor_list_iterator_release,&
49 : neighbor_list_set_p_type
50 : #include "./base/base_uses.f90"
51 :
52 : IMPLICIT NONE
53 : PRIVATE
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_operators_ao'
56 :
57 : ! *** Public subroutines ***
58 :
59 : PUBLIC :: p_xyz_ao, rRc_xyz_ao, rRc_xyz_der_ao
60 : PUBLIC :: build_lin_mom_matrix, build_ang_mom_matrix
61 :
62 : CONTAINS
63 :
64 : ! **************************************************************************************************
65 : !> \brief Calculation of the linear momentum matrix <mu|∂|nu> over
66 : !> Cartesian Gaussian functions.
67 : !> \param qs_env ...
68 : !> \param matrix ...
69 : !> \date 27.02.2009
70 : !> \author VW
71 : !> \version 1.0
72 : ! **************************************************************************************************
73 930 : SUBROUTINE build_lin_mom_matrix(qs_env, matrix)
74 :
75 : TYPE(qs_environment_type), POINTER :: qs_env
76 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
77 :
78 : CHARACTER(len=*), PARAMETER :: routineN = 'build_lin_mom_matrix'
79 :
80 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
81 : ldai, maxco, maxlgto, maxsgf, natom, ncoa, ncob, neighbor_list_id, nkind, nseta, nsetb, &
82 : sgfa, sgfb
83 930 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
84 930 : npgfb, nsgfa, nsgfb
85 930 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
86 : LOGICAL :: do_symmetric, found, new_atom_b
87 : REAL(KIND=dp) :: dab, rab2
88 930 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
89 930 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab, rr_work
90 : REAL(KIND=dp), DIMENSION(3) :: rab
91 930 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
92 930 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
93 930 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
94 : TYPE(cell_type), POINTER :: cell
95 : TYPE(cp_logger_type), POINTER :: logger
96 930 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
97 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
98 : TYPE(mp_para_env_type), POINTER :: para_env
99 : TYPE(neighbor_list_iterator_p_type), &
100 930 : DIMENSION(:), POINTER :: nl_iterator
101 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
102 930 : POINTER :: sab_nl
103 930 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
104 930 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
105 : TYPE(qs_kind_type), POINTER :: qs_kind
106 :
107 930 : CALL timeset(routineN, handle)
108 :
109 930 : NULLIFY (cell, sab_nl, qs_kind_set, particle_set, para_env)
110 930 : NULLIFY (logger)
111 :
112 930 : logger => cp_get_default_logger()
113 :
114 : CALL get_qs_env(qs_env=qs_env, &
115 : qs_kind_set=qs_kind_set, &
116 : particle_set=particle_set, &
117 : neighbor_list_id=neighbor_list_id, &
118 : para_env=para_env, &
119 930 : cell=cell)
120 :
121 930 : nkind = SIZE(qs_kind_set)
122 930 : natom = SIZE(particle_set)
123 :
124 : ! Take into account the symmetry of the input matrix
125 930 : do_symmetric = dbcsr_has_symmetry(matrix(1)%matrix)
126 930 : IF (do_symmetric) THEN
127 928 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
128 : ELSE
129 2 : CALL get_qs_env(qs_env=qs_env, sab_all=sab_nl)
130 : END IF
131 : ! *** Allocate work storage ***
132 :
133 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
134 : maxco=maxco, &
135 : maxlgto=maxlgto, &
136 930 : maxsgf=maxsgf)
137 :
138 930 : ldai = ncoset(maxlgto + 1)
139 930 : CALL init_orbital_pointers(ldai)
140 :
141 13950 : ALLOCATE (rr_work(ldai, ldai, 3), intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
142 789900 : rr_work(:, :, :) = 0.0_dp
143 649476 : intab(:, :, :) = 0.0_dp
144 149076 : work(:, :) = 0.0_dp
145 :
146 4402 : ALLOCATE (basis_set_list(nkind))
147 2542 : DO ikind = 1, nkind
148 1612 : qs_kind => qs_kind_set(ikind)
149 1612 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
150 2542 : IF (ASSOCIATED(basis_set_a)) THEN
151 1612 : basis_set_list(ikind)%gto_basis_set => basis_set_a
152 : ELSE
153 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
154 : END IF
155 : END DO
156 930 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
157 53301 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
158 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
159 52371 : iatom=iatom, jatom=jatom, r=rab)
160 52371 : basis_set_a => basis_set_list(ikind)%gto_basis_set
161 52371 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
162 52371 : basis_set_b => basis_set_list(jkind)%gto_basis_set
163 52371 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
164 : ! basis ikind
165 52371 : first_sgfa => basis_set_a%first_sgf
166 52371 : la_max => basis_set_a%lmax
167 52371 : la_min => basis_set_a%lmin
168 52371 : npgfa => basis_set_a%npgf
169 52371 : nseta = basis_set_a%nset
170 52371 : nsgfa => basis_set_a%nsgf_set
171 52371 : rpgfa => basis_set_a%pgf_radius
172 52371 : set_radius_a => basis_set_a%set_radius
173 52371 : sphi_a => basis_set_a%sphi
174 52371 : zeta => basis_set_a%zet
175 : ! basis jkind
176 52371 : first_sgfb => basis_set_b%first_sgf
177 52371 : lb_max => basis_set_b%lmax
178 52371 : lb_min => basis_set_b%lmin
179 52371 : npgfb => basis_set_b%npgf
180 52371 : nsetb = basis_set_b%nset
181 52371 : nsgfb => basis_set_b%nsgf_set
182 52371 : rpgfb => basis_set_b%pgf_radius
183 52371 : set_radius_b => basis_set_b%set_radius
184 52371 : sphi_b => basis_set_b%sphi
185 52371 : zetb => basis_set_b%zet
186 :
187 52371 : IF (inode == 1) last_jatom = 0
188 52371 : IF (jatom /= last_jatom) THEN
189 : new_atom_b = .TRUE.
190 : last_jatom = jatom
191 : ELSE
192 : new_atom_b = .FALSE.
193 : END IF
194 :
195 : IF (new_atom_b) THEN
196 3208 : IF (do_symmetric) THEN
197 3199 : IF (iatom <= jatom) THEN
198 2014 : irow = iatom
199 2014 : icol = jatom
200 : ELSE
201 1185 : irow = jatom
202 1185 : icol = iatom
203 : END IF
204 : ELSE
205 9 : irow = iatom
206 9 : icol = jatom
207 : END IF
208 :
209 12832 : DO i = 1, 3
210 9624 : NULLIFY (integral(i)%block)
211 : CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
212 9624 : row=irow, col=icol, BLOCK=integral(i)%block, found=found)
213 12832 : CPASSERT(ASSOCIATED(INTEGRAL(i)%block))
214 : END DO
215 : END IF
216 :
217 52371 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
218 52371 : dab = SQRT(rab2)
219 :
220 170837 : DO iset = 1, nseta
221 :
222 117536 : ncoa = npgfa(iset)*ncoset(la_max(iset))
223 117536 : sgfa = first_sgfa(1, iset)
224 :
225 444053 : DO jset = 1, nsetb
226 :
227 274146 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
228 :
229 109994 : ncob = npgfb(jset)*ncoset(lb_max(jset))
230 109994 : sgfb = first_sgfb(1, jset)
231 :
232 : ! *** Calculate the primitive fermi contact integrals ***
233 :
234 : CALL lin_mom(la_max(iset), la_min(iset), npgfa(iset), &
235 : rpgfa(:, iset), zeta(:, iset), &
236 : lb_max(jset), lb_min(jset), npgfb(jset), &
237 : rpgfb(:, jset), zetb(:, jset), &
238 109994 : rab, intab, SIZE(rr_work, 1), rr_work)
239 :
240 : ! *** Contraction step ***
241 :
242 557512 : DO i = 1, 3
243 :
244 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
245 : 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
246 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
247 329982 : 0.0_dp, work(1, 1), SIZE(work, 1))
248 :
249 604128 : IF (do_symmetric) THEN
250 329955 : IF (iatom <= jatom) THEN
251 :
252 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
253 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
254 : work(1, 1), SIZE(work, 1), &
255 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
256 207738 : SIZE(integral(i)%block, 1))
257 :
258 : ELSE
259 :
260 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
261 : -1.0_dp, work(1, 1), SIZE(work, 1), &
262 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
263 : 1.0_dp, integral(i)%block(sgfb, sgfa), &
264 122217 : SIZE(integral(i)%block, 1))
265 :
266 : END IF
267 : ELSE
268 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
269 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
270 : work(1, 1), SIZE(work, 1), &
271 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
272 27 : SIZE(integral(i)%block, 1))
273 : END IF
274 :
275 : END DO
276 :
277 : END DO
278 :
279 : END DO
280 :
281 : END DO
282 930 : CALL neighbor_list_iterator_release(nl_iterator)
283 :
284 : ! *** Release work storage ***
285 :
286 930 : DEALLOCATE (intab, work, integral, basis_set_list)
287 :
288 : ! *** Print the spin orbit matrix, if requested ***
289 :
290 : !IF (BTEST(cp_print_key_should_output(logger%iter_info,&
291 : ! qs_env%input,"DFT%PRINT%AO_MATRICES/LINEAR_MOMENTUM"),cp_p_file)) THEN
292 : ! iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/LINEA_MOMENTUM",&
293 : ! extension=".Log")
294 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(1)%matrix,4,6,qs_env,para_env,output_unit=iw)
295 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(2)%matrix,4,6,qs_env,para_env,output_unit=iw)
296 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(3)%matrix,4,6,qs_env,para_env,output_unit=iw)
297 : ! CALL cp_print_key_finished_output(iw,logger,qs_env%input,&
298 : ! "DFT%PRINT%AO_MATRICES/LINEAR_MOMENTUM")
299 : !END IF
300 :
301 930 : CALL timestop(handle)
302 :
303 2790 : END SUBROUTINE build_lin_mom_matrix
304 :
305 : ! **************************************************************************************************
306 : !> \brief Calculation of the primitive paramagnetic spin orbit integrals over
307 : !> Cartesian Gaussian-type functions.
308 : !> \param la_max ...
309 : !> \param la_min ...
310 : !> \param npgfa ...
311 : !> \param rpgfa ...
312 : !> \param zeta ...
313 : !> \param lb_max ...
314 : !> \param lb_min ...
315 : !> \param npgfb ...
316 : !> \param rpgfb ...
317 : !> \param zetb ...
318 : !> \param rab ...
319 : !> \param intab ...
320 : !> \param ldrr ...
321 : !> \param rr ...
322 : !> \date 02.03.2009
323 : !> \author VW
324 : !> \version 1.0
325 : ! **************************************************************************************************
326 219988 : SUBROUTINE lin_mom(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, &
327 109994 : rab, intab, ldrr, rr)
328 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
329 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
330 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
331 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
332 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
333 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: intab
334 : INTEGER, INTENT(IN) :: ldrr
335 : REAL(dp), DIMENSION(0:ldrr-1, 0:ldrr-1, 3), &
336 : INTENT(INOUT) :: rr
337 :
338 : INTEGER :: ax, ay, az, bx, by, bz, coa, cob, i, &
339 : ipgf, j, jpgf, la, lb, ma, mb, na, nb
340 : REAL(dp) :: dab, dumx, dumy, dumz, f0, rab2, xhi, zet
341 : REAL(dp), DIMENSION(3) :: rap, rbp
342 :
343 : ! *** Calculate the distance of the centers a and c ***
344 :
345 109994 : rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
346 109994 : dab = SQRT(rab2)
347 :
348 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
349 :
350 109994 : na = 0
351 :
352 354886 : DO ipgf = 1, npgfa
353 :
354 244892 : nb = 0
355 :
356 873734 : DO jpgf = 1, npgfb
357 :
358 : ! *** Screening ***
359 :
360 628842 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
361 1508461 : DO j = nb + 1, nb + ncoset(lb_max)
362 5166579 : DO i = na + 1, na + ncoset(la_max)
363 3658118 : intab(i, j, 1) = 0.0_dp
364 3658118 : intab(i, j, 2) = 0.0_dp
365 4763494 : intab(i, j, 3) = 0.0_dp
366 : END DO
367 : END DO
368 403085 : nb = nb + ncoset(lb_max)
369 403085 : CYCLE
370 : END IF
371 :
372 : ! *** Calculate some prefactors ***
373 225757 : zet = zeta(ipgf) + zetb(jpgf)
374 225757 : xhi = zeta(ipgf)*zetb(jpgf)/zet
375 903028 : rap = zetb(jpgf)*rab/zet
376 903028 : rbp = -zeta(ipgf)*rab/zet
377 :
378 225757 : f0 = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
379 :
380 : ! *** Calculate the recurrence relation ***
381 :
382 225757 : CALL os_rr_ovlp(rap, la_max + 1, rbp, lb_max, zet, ldrr, rr)
383 :
384 : ! *** Calculate the primitive linear momentum integrals ***
385 538526 : DO lb = lb_min, lb_max
386 972225 : DO bx = 0, lb
387 1313015 : DO by = 0, lb - bx
388 566547 : bz = lb - bx - by
389 566547 : cob = coset(bx, by, bz)
390 566547 : mb = nb + cob
391 1862245 : DO la = la_min, la_max
392 2688722 : DO ax = 0, la
393 3824563 : DO ay = 0, la - ax
394 1702388 : az = la - ax - ay
395 1702388 : coa = coset(ax, ay, az)
396 1702388 : ma = na + coa
397 : !
398 : !
399 : ! (a|p_x|b) = 2*a*(a+1x|b) - N_x(a)*(a-1_x|b)
400 1702388 : dumx = 2.0_dp*zeta(ipgf)*rr(ax + 1, bx, 1)
401 1702388 : IF (ax .GT. 0) dumx = dumx - REAL(ax, dp)*rr(ax - 1, bx, 1)
402 1702388 : intab(ma, mb, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
403 : !
404 : ! (a|p_y|b)
405 1702388 : dumy = 2.0_dp*zeta(ipgf)*rr(ay + 1, by, 2)
406 1702388 : IF (ay .GT. 0) dumy = dumy - REAL(ay, dp)*rr(ay - 1, by, 2)
407 1702388 : intab(ma, mb, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
408 : !
409 : ! (a|p_z|b)
410 1702388 : dumz = 2.0_dp*zeta(ipgf)*rr(az + 1, bz, 3)
411 1702388 : IF (az .GT. 0) dumz = dumz - REAL(az, dp)*rr(az - 1, bz, 3)
412 2962564 : intab(ma, mb, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
413 : !
414 : END DO
415 : END DO
416 : END DO !la
417 :
418 : END DO
419 : END DO
420 : END DO !lb
421 :
422 470649 : nb = nb + ncoset(lb_max)
423 :
424 : END DO
425 :
426 354886 : na = na + ncoset(la_max)
427 :
428 : END DO
429 :
430 109994 : END SUBROUTINE lin_mom
431 :
432 : ! **************************************************************************************************
433 : !> \brief Calculation of the angular momentum matrix over
434 : !> Cartesian Gaussian functions.
435 : !> \param qs_env ...
436 : !> \param matrix ...
437 : !> \param rc ...
438 : !> \date 27.02.2009
439 : !> \author VW
440 : !> \version 1.0
441 : ! **************************************************************************************************
442 :
443 1250 : SUBROUTINE build_ang_mom_matrix(qs_env, matrix, rc)
444 :
445 : TYPE(qs_environment_type), POINTER :: qs_env
446 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
447 : REAL(dp), DIMENSION(:), INTENT(IN) :: rc
448 :
449 : CHARACTER(len=*), PARAMETER :: routineN = 'build_ang_mom_matrix'
450 :
451 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
452 : ldai, maxco, maxlgto, maxsgf, natom, ncoa, ncob, neighbor_list_id, nkind, nseta, nsetb, &
453 : sgfa, sgfb
454 1250 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
455 1250 : npgfb, nsgfa, nsgfb
456 1250 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
457 : LOGICAL :: found, new_atom_b
458 : REAL(KIND=dp) :: dab, rab2
459 1250 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
460 1250 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab, rr_work
461 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rbc
462 1250 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
463 1250 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
464 1250 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
465 : TYPE(cell_type), POINTER :: cell
466 : TYPE(cp_logger_type), POINTER :: logger
467 1250 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
468 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
469 : TYPE(mp_para_env_type), POINTER :: para_env
470 : TYPE(neighbor_list_iterator_p_type), &
471 1250 : DIMENSION(:), POINTER :: nl_iterator
472 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
473 1250 : POINTER :: sab_all
474 1250 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
475 1250 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
476 : TYPE(qs_kind_type), POINTER :: qs_kind
477 :
478 1250 : CALL timeset(routineN, handle)
479 :
480 1250 : NULLIFY (cell, sab_all, qs_kind_set, particle_set, para_env)
481 1250 : NULLIFY (logger)
482 :
483 1250 : logger => cp_get_default_logger()
484 :
485 : CALL get_qs_env(qs_env=qs_env, &
486 : qs_kind_set=qs_kind_set, &
487 : particle_set=particle_set, &
488 : neighbor_list_id=neighbor_list_id, &
489 : para_env=para_env, &
490 : sab_all=sab_all, &
491 1250 : cell=cell)
492 :
493 1250 : nkind = SIZE(qs_kind_set)
494 1250 : natom = SIZE(particle_set)
495 :
496 : ! *** Allocate work storage ***
497 :
498 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
499 : maxco=maxco, &
500 : maxlgto=maxlgto, &
501 1250 : maxsgf=maxsgf)
502 :
503 1250 : ldai = ncoset(maxlgto + 1)
504 1250 : CALL init_orbital_pointers(ldai)
505 :
506 18750 : ALLOCATE (rr_work(ldai, ldai, 3), intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
507 1130660 : rr_work(:, :, :) = 0.0_dp
508 673460 : intab(:, :, :) = 0.0_dp
509 201912 : work(:, :) = 0.0_dp
510 :
511 5796 : ALLOCATE (basis_set_list(nkind))
512 3296 : DO ikind = 1, nkind
513 2046 : qs_kind => qs_kind_set(ikind)
514 2046 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
515 3296 : IF (ASSOCIATED(basis_set_a)) THEN
516 2046 : basis_set_list(ikind)%gto_basis_set => basis_set_a
517 : ELSE
518 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
519 : END IF
520 : END DO
521 1250 : CALL neighbor_list_iterator_create(nl_iterator, sab_all)
522 94981 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
523 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
524 93731 : iatom=iatom, jatom=jatom, r=rab)
525 93731 : basis_set_a => basis_set_list(ikind)%gto_basis_set
526 93731 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
527 93731 : basis_set_b => basis_set_list(jkind)%gto_basis_set
528 93731 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
529 93731 : ra = pbc(particle_set(iatom)%r, cell)
530 : ! basis ikind
531 93731 : first_sgfa => basis_set_a%first_sgf
532 93731 : la_max => basis_set_a%lmax
533 93731 : la_min => basis_set_a%lmin
534 93731 : npgfa => basis_set_a%npgf
535 93731 : nseta = basis_set_a%nset
536 93731 : nsgfa => basis_set_a%nsgf_set
537 93731 : rpgfa => basis_set_a%pgf_radius
538 93731 : set_radius_a => basis_set_a%set_radius
539 93731 : sphi_a => basis_set_a%sphi
540 93731 : zeta => basis_set_a%zet
541 : ! basis jkind
542 93731 : first_sgfb => basis_set_b%first_sgf
543 93731 : lb_max => basis_set_b%lmax
544 93731 : lb_min => basis_set_b%lmin
545 93731 : npgfb => basis_set_b%npgf
546 93731 : nsetb = basis_set_b%nset
547 93731 : nsgfb => basis_set_b%nsgf_set
548 93731 : rpgfb => basis_set_b%pgf_radius
549 93731 : set_radius_b => basis_set_b%set_radius
550 93731 : sphi_b => basis_set_b%sphi
551 93731 : zetb => basis_set_b%zet
552 :
553 93731 : IF (inode == 1) last_jatom = 0
554 :
555 93731 : IF (jatom /= last_jatom) THEN
556 : new_atom_b = .TRUE.
557 : last_jatom = jatom
558 : ELSE
559 : new_atom_b = .FALSE.
560 : END IF
561 :
562 : IF (new_atom_b) THEN
563 : !IF (iatom <= jatom) THEN
564 5987 : irow = iatom
565 5987 : icol = jatom
566 : !ELSE
567 : ! irow = jatom
568 : ! icol = iatom
569 : !END IF
570 :
571 23948 : DO i = 1, 3
572 17961 : NULLIFY (INTEGRAL(i)%block)
573 : CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
574 17961 : row=irow, col=icol, BLOCK=INTEGRAL(i)%block, found=found)
575 23948 : CPASSERT(ASSOCIATED(INTEGRAL(i)%block))
576 : END DO
577 : END IF
578 :
579 93731 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
580 93731 : dab = SQRT(rab2)
581 :
582 312058 : DO iset = 1, nseta
583 :
584 217077 : ncoa = npgfa(iset)*ncoset(la_max(iset))
585 217077 : sgfa = first_sgfa(1, iset)
586 :
587 837307 : DO jset = 1, nsetb
588 :
589 526499 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
590 :
591 : !IF(PRESENT(wancen)) THEN
592 : ! rc = wancen
593 201893 : rac = pbc(rc, ra, cell)
594 807572 : rbc = rac + rab
595 : !ELSE
596 : ! rc(1:3) = rb(1:3)
597 : ! rac(1:3) = -rab(1:3)
598 : ! rbc(1:3) = 0.0_dp
599 : !ENDIF
600 :
601 201893 : ncob = npgfb(jset)*ncoset(lb_max(jset))
602 201893 : sgfb = first_sgfb(1, jset)
603 :
604 : ! *** Calculate the primitive angular momentum integrals ***
605 :
606 : CALL ang_mom(la_max(iset), la_min(iset), npgfa(iset), &
607 : rpgfa(:, iset), zeta(:, iset), &
608 : lb_max(jset), lb_min(jset), npgfb(jset), &
609 : rpgfb(:, jset), zetb(:, jset), &
610 201893 : rab, rac, intab, SIZE(rr_work, 1), rr_work)
611 :
612 : ! *** Contraction step ***
613 :
614 1024649 : DO i = 1, 3
615 :
616 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
617 : 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
618 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
619 605679 : 0.0_dp, work(1, 1), SIZE(work, 1))
620 :
621 : !IF (iatom <= jatom) THEN
622 :
623 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
624 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
625 : work(1, 1), SIZE(work, 1), &
626 : 1.0_dp, integral(i)%block(sgfa, sgfb), &
627 1132178 : SIZE(integral(i)%block, 1))
628 :
629 : !ELSE
630 : !
631 : ! CALL dgemm("T","N",nsgfb(jset),nsgfa(iset),ncoa,&
632 : ! -1.0_dp,work(1,1),SIZE(work,1),&
633 : ! sphi_a(1,sgfa),SIZE(sphi_a,1),&
634 : ! 1.0_dp,integral(i)%block(sgfb,sgfa),&
635 : ! SIZE(integral(i)%block,1))
636 : !
637 : !ENDIF
638 :
639 : END DO
640 :
641 : END DO
642 :
643 : END DO
644 :
645 : END DO
646 1250 : CALL neighbor_list_iterator_release(nl_iterator)
647 :
648 : ! *** Release work storage ***
649 :
650 1250 : DEALLOCATE (intab, work, integral, basis_set_list)
651 :
652 : ! *** Print the spin orbit matrix, if requested ***
653 :
654 : !IF (BTEST(cp_print_key_should_output(logger%iter_info,&
655 : ! qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM"),cp_p_file)) THEN
656 : ! iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM",&
657 : ! extension=".Log")
658 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(1)%matrix,4,6,qs_env,para_env,output_unit=iw)
659 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(2)%matrix,4,6,qs_env,para_env,output_unit=iw)
660 : ! CALL cp_dbcsr_write_sparse_matrix(matrix(3)%matrix,4,6,qs_env,para_env,output_unit=iw)
661 : ! CALL cp_print_key_finished_output(iw,logger,qs_env%input,&
662 : ! "DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM")
663 : !END IF
664 :
665 1250 : CALL timestop(handle)
666 :
667 3750 : END SUBROUTINE build_ang_mom_matrix
668 :
669 : ! **************************************************************************************************
670 : !> \brief Calculation of the primitive paramagnetic spin orbit integrals over
671 : !> Cartesian Gaussian-type functions.
672 : !> \param la_max ...
673 : !> \param la_min ...
674 : !> \param npgfa ...
675 : !> \param rpgfa ...
676 : !> \param zeta ...
677 : !> \param lb_max ...
678 : !> \param lb_min ...
679 : !> \param npgfb ...
680 : !> \param rpgfb ...
681 : !> \param zetb ...
682 : !> \param rab ...
683 : !> \param rac ...
684 : !> \param intab ...
685 : !> \param ldrr ...
686 : !> \param rr ...
687 : !> \date 02.03.2009
688 : !> \author VW
689 : !> \version 1.0
690 : ! **************************************************************************************************
691 403786 : SUBROUTINE ang_mom(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, &
692 201893 : rab, rac, intab, ldrr, rr)
693 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
694 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
695 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
696 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
697 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab, rac
698 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: intab
699 : INTEGER, INTENT(IN) :: ldrr
700 : REAL(dp), DIMENSION(0:ldrr-1, 0:ldrr-1, 3), &
701 : INTENT(INOUT) :: rr
702 :
703 : INTEGER :: ax, ay, az, bx, by, bz, coa, cob, i, &
704 : ipgf, j, jpgf, la, lb, ma, mb, na, nb
705 : REAL(dp) :: dab, dumx, dumy, dumz, f0, rab2, xhi, zet
706 : REAL(dp), DIMENSION(3) :: rap, rbp
707 :
708 : ! *** Calculate the distance of the centers a and c ***
709 :
710 201893 : rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
711 201893 : dab = SQRT(rab2)
712 :
713 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
714 :
715 201893 : na = 0
716 :
717 676304 : DO ipgf = 1, npgfa
718 :
719 474411 : nb = 0
720 :
721 1742562 : DO jpgf = 1, npgfb
722 :
723 : ! *** Screening ***
724 :
725 1268151 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
726 3336529 : DO j = nb + 1, nb + ncoset(lb_max)
727 12446889 : DO i = na + 1, na + ncoset(la_max)
728 9110360 : intab(i, j, 1) = 0.0_dp
729 9110360 : intab(i, j, 2) = 0.0_dp
730 11620957 : intab(i, j, 3) = 0.0_dp
731 : END DO
732 : END DO
733 825932 : nb = nb + ncoset(lb_max)
734 825932 : CYCLE
735 : END IF
736 :
737 : ! *** Calculate some prefactors ***
738 442219 : zet = zeta(ipgf) + zetb(jpgf)
739 442219 : xhi = zeta(ipgf)*zetb(jpgf)/zet
740 1768876 : rap = zetb(jpgf)*rab/zet
741 1768876 : rbp = -zeta(ipgf)*rab/zet
742 :
743 442219 : f0 = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
744 :
745 : ! *** Calculate the recurrence relation ***
746 :
747 442219 : CALL os_rr_ovlp(rap, la_max + 1, rbp, lb_max, zet, ldrr, rr)
748 :
749 : ! *** Calculate the primitive Fermi contact integrals ***
750 :
751 1087122 : DO lb = lb_min, lb_max
752 2011607 : DO bx = 0, lb
753 2802912 : DO by = 0, lb - bx
754 1233524 : bz = lb - bx - by
755 1233524 : cob = coset(bx, by, bz)
756 1233524 : mb = nb + cob
757 4112005 : DO la = la_min, la_max
758 6109494 : DO ax = 0, la
759 8873886 : DO ay = 0, la - ax
760 3997916 : az = la - ax - ay
761 3997916 : coa = coset(ax, ay, az)
762 3997916 : ma = na + coa
763 : !
764 3997916 : dumx = -2.0_dp*zeta(ipgf)*rr(ax + 1, bx, 1)
765 3997916 : dumy = -2.0_dp*zeta(ipgf)*rr(ay + 1, by, 2)
766 3997916 : dumz = -2.0_dp*zeta(ipgf)*rr(az + 1, bz, 3)
767 3997916 : IF (ax .GT. 0) dumx = dumx + REAL(ax, dp)*rr(ax - 1, bx, 1)
768 3997916 : IF (ay .GT. 0) dumy = dumy + REAL(ay, dp)*rr(ay - 1, by, 2)
769 3997916 : IF (az .GT. 0) dumz = dumz + REAL(az, dp)*rr(az - 1, bz, 3)
770 : !
771 : ! (a|l_z|b)
772 : intab(ma, mb, 1) = -f0*rr(ax, bx, 1)*( &
773 : & (rr(ay + 1, by, 2) + rac(2)*rr(ay, by, 2))*dumz &
774 3997916 : & - (rr(az + 1, bz, 3) + rac(3)*rr(az, bz, 3))*dumy)
775 : !
776 : ! (a|l_y|b)
777 : intab(ma, mb, 2) = -f0*rr(ay, by, 2)*( &
778 : & (rr(az + 1, bz, 3) + rac(3)*rr(az, bz, 3))*dumx &
779 3997916 : & - (rr(ax + 1, bx, 1) + rac(1)*rr(ax, bx, 1))*dumz)
780 : !
781 : ! (a|l_z|b)
782 : intab(ma, mb, 3) = -f0*rr(az, bz, 3)*( &
783 : & (rr(ax + 1, bx, 1) + rac(1)*rr(ax, bx, 1))*dumy &
784 6919890 : & - (rr(ay + 1, by, 2) + rac(2)*rr(ay, by, 2))*dumx)
785 : !
786 : END DO
787 : END DO
788 : END DO !la
789 :
790 : END DO
791 : END DO
792 : END DO !lb
793 :
794 916630 : nb = nb + ncoset(lb_max)
795 :
796 : END DO
797 :
798 676304 : na = na + ncoset(la_max)
799 :
800 : END DO
801 :
802 201893 : END SUBROUTINE ang_mom
803 :
804 : ! **************************************************************************************************
805 : !> \brief Calculation of the components of the dipole operator in the velocity form
806 : !> The elements of the sparse matrices are the integrals in the
807 : !> basis functions
808 : !> \param op matrix representation of the p operator
809 : !> calculated in terms of the contracted basis functions
810 : !> \param qs_env environment for the lists and the basis sets
811 : !> \param minimum_image take into account only the first neighbors in the lists
812 : !> \par History
813 : !> 06.2005 created [MI]
814 : !> \author MI
815 : ! **************************************************************************************************
816 :
817 74 : SUBROUTINE p_xyz_ao(op, qs_env, minimum_image)
818 :
819 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op
820 : TYPE(qs_environment_type), POINTER :: qs_env
821 : LOGICAL, INTENT(IN), OPTIONAL :: minimum_image
822 :
823 : CHARACTER(len=*), PARAMETER :: routineN = 'p_xyz_ao'
824 :
825 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
826 : ldab, ldsa, ldsb, ldwork, maxl, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
827 74 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
828 74 : npgfb, nsgfa, nsgfb
829 74 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
830 : LOGICAL :: found, my_minimum_image, new_atom_b
831 : REAL(KIND=dp) :: alpha, dab, Lxo2, Lyo2, Lzo2, rab2
832 : REAL(KIND=dp), DIMENSION(3) :: ra, rab
833 74 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
834 74 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, &
835 74 : zeta, zetb
836 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: difab
837 74 : TYPE(block_p_type), DIMENSION(:), POINTER :: op_dip
838 : TYPE(cell_type), POINTER :: cell
839 74 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
840 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
841 : TYPE(neighbor_list_iterator_p_type), &
842 74 : DIMENSION(:), POINTER :: nl_iterator
843 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
844 74 : POINTER :: sab_orb
845 74 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
846 74 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
847 : TYPE(qs_kind_type), POINTER :: qs_kind
848 :
849 74 : CALL timeset(routineN, handle)
850 :
851 74 : NULLIFY (qs_kind, qs_kind_set)
852 74 : NULLIFY (cell, particle_set)
853 74 : NULLIFY (sab_orb)
854 74 : NULLIFY (difab, op_dip, work)
855 74 : NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb)
856 74 : NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
857 :
858 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
859 : cell=cell, particle_set=particle_set, &
860 74 : sab_orb=sab_orb)
861 :
862 74 : nkind = SIZE(qs_kind_set)
863 :
864 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
865 74 : maxco=ldwork, maxlgto=maxl)
866 :
867 74 : my_minimum_image = .FALSE.
868 74 : IF (PRESENT(minimum_image)) THEN
869 32 : my_minimum_image = minimum_image
870 128 : Lxo2 = SQRT(SUM(cell%hmat(:, 1)**2))/2.0_dp
871 128 : Lyo2 = SQRT(SUM(cell%hmat(:, 2)**2))/2.0_dp
872 128 : Lzo2 = SQRT(SUM(cell%hmat(:, 3)**2))/2.0_dp
873 : END IF
874 :
875 74 : ldab = ldwork
876 :
877 370 : ALLOCATE (difab(ldab, ldab, 3))
878 73952 : difab(1:ldab, 1:ldab, 1:3) = 0.0_dp
879 296 : ALLOCATE (work(ldwork, ldwork))
880 24626 : work(1:ldwork, 1:ldwork) = 0.0_dp
881 296 : ALLOCATE (op_dip(3))
882 :
883 296 : DO i = 1, 3
884 296 : NULLIFY (op_dip(i)%block)
885 : END DO
886 :
887 342 : ALLOCATE (basis_set_list(nkind))
888 194 : DO ikind = 1, nkind
889 120 : qs_kind => qs_kind_set(ikind)
890 120 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
891 194 : IF (ASSOCIATED(basis_set_a)) THEN
892 120 : basis_set_list(ikind)%gto_basis_set => basis_set_a
893 : ELSE
894 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
895 : END IF
896 : END DO
897 74 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
898 9841 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
899 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
900 9767 : iatom=iatom, jatom=jatom, r=rab)
901 9767 : basis_set_a => basis_set_list(ikind)%gto_basis_set
902 9767 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
903 9767 : basis_set_b => basis_set_list(jkind)%gto_basis_set
904 9767 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
905 9767 : ra = pbc(particle_set(iatom)%r, cell)
906 : ! basis ikind
907 9767 : first_sgfa => basis_set_a%first_sgf
908 9767 : la_max => basis_set_a%lmax
909 9767 : la_min => basis_set_a%lmin
910 9767 : npgfa => basis_set_a%npgf
911 9767 : nseta = basis_set_a%nset
912 9767 : nsgfa => basis_set_a%nsgf_set
913 9767 : rpgfa => basis_set_a%pgf_radius
914 9767 : set_radius_a => basis_set_a%set_radius
915 9767 : sphi_a => basis_set_a%sphi
916 9767 : zeta => basis_set_a%zet
917 : ! basis jkind
918 9767 : first_sgfb => basis_set_b%first_sgf
919 9767 : lb_max => basis_set_b%lmax
920 9767 : lb_min => basis_set_b%lmin
921 9767 : npgfb => basis_set_b%npgf
922 9767 : nsetb = basis_set_b%nset
923 9767 : nsgfb => basis_set_b%nsgf_set
924 9767 : rpgfb => basis_set_b%pgf_radius
925 9767 : set_radius_b => basis_set_b%set_radius
926 9767 : sphi_b => basis_set_b%sphi
927 9767 : zetb => basis_set_b%zet
928 :
929 9767 : IF (inode == 1) THEN
930 363 : last_jatom = 0
931 363 : alpha = 1.0_dp
932 : END IF
933 9767 : ldsa = SIZE(sphi_a, 1)
934 9767 : ldsb = SIZE(sphi_b, 1)
935 :
936 9767 : IF (my_minimum_image) THEN
937 8345 : IF (ABS(rab(1)) > Lxo2 .OR. ABS(rab(2)) > Lyo2 .OR. ABS(rab(3)) > Lzo2) CYCLE
938 : END IF
939 :
940 6303 : IF (jatom /= last_jatom) THEN
941 : new_atom_b = .TRUE.
942 : last_jatom = jatom
943 : ELSE
944 : new_atom_b = .FALSE.
945 : END IF
946 :
947 : IF (new_atom_b) THEN
948 4921 : IF (iatom <= jatom) THEN
949 2532 : irow = iatom
950 2532 : icol = jatom
951 2532 : alpha = 1.0_dp
952 : ELSE
953 2389 : irow = jatom
954 2389 : icol = iatom
955 2389 : IF (dbcsr_get_matrix_type(op(1)%matrix) == dbcsr_type_antisymmetric) THEN
956 : !IF(op(1)%matrix%symmetry=="antisymmetric") THEN
957 4778 : alpha = -1.0_dp
958 : END IF
959 : END IF
960 :
961 19684 : DO i = 1, 3
962 14763 : NULLIFY (op_dip(i)%block)
963 : CALL dbcsr_get_block_p(matrix=op(i)%matrix, &
964 14763 : row=irow, col=icol, block=op_dip(i)%block, found=found)
965 19684 : CPASSERT(ASSOCIATED(op_dip(i)%block))
966 : END DO
967 : END IF ! new_atom_b
968 6303 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
969 6303 : dab = SQRT(rab2)
970 :
971 23726 : DO iset = 1, nseta
972 :
973 17349 : ncoa = npgfa(iset)*ncoset(la_max(iset))
974 17349 : sgfa = first_sgfa(1, iset)
975 :
976 78981 : DO jset = 1, nsetb
977 :
978 51865 : ncob = npgfb(jset)*ncoset(lb_max(jset))
979 51865 : sgfb = first_sgfb(1, jset)
980 :
981 69214 : IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN
982 :
983 : ! *** Calculate the primitive overlap integrals ***
984 : CALL diffop(la_max(iset), npgfa(iset), zeta(:, iset), &
985 : rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
986 22848 : zetb(:, jset), rpgfb(:, jset), lb_min(jset), rab, difab)
987 :
988 : ! *** Contraction ***
989 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
990 : alpha, difab(1, 1, 1), ldab, sphi_b(1, sgfb), ldsb, &
991 22848 : 0.0_dp, work(1, 1), ldwork)
992 22848 : IF (iatom <= jatom) THEN
993 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
994 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
995 : work(1, 1), ldwork, &
996 : 1.0_dp, op_dip(1)%block(sgfa, sgfb), &
997 12802 : SIZE(op_dip(1)%block, 1))
998 :
999 : ELSE
1000 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1001 : 1.0_dp, work(1, 1), ldwork, &
1002 : sphi_a(1, sgfa), ldsa, &
1003 : 1.0_dp, op_dip(1)%block(sgfb, sgfa), &
1004 10046 : SIZE(op_dip(1)%block, 1))
1005 :
1006 : END IF
1007 :
1008 : ! *** Contraction ***
1009 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1010 : alpha, difab(1, 1, 2), ldab, sphi_b(1, sgfb), ldsb, &
1011 22848 : 0.0_dp, work(1, 1), ldwork)
1012 22848 : IF (iatom <= jatom) THEN
1013 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1014 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
1015 : work(1, 1), ldwork, &
1016 : 1.0_dp, op_dip(2)%block(sgfa, sgfb), &
1017 12802 : SIZE(op_dip(2)%block, 1))
1018 : ELSE
1019 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1020 : 1.0_dp, work(1, 1), ldwork, &
1021 : sphi_a(1, sgfa), ldsa, &
1022 : 1.0_dp, op_dip(2)%block(sgfb, sgfa), &
1023 10046 : SIZE(op_dip(2)%block, 1))
1024 : END IF
1025 :
1026 : ! *** Contraction ***
1027 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1028 : alpha, difab(1, 1, 3), ldab, sphi_b(1, sgfb), ldsb, &
1029 22848 : 0.0_dp, work(1, 1), ldwork)
1030 22848 : IF (iatom <= jatom) THEN
1031 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1032 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
1033 : work(1, 1), ldwork, &
1034 : 1.0_dp, op_dip(3)%block(sgfa, sgfb), &
1035 12802 : SIZE(op_dip(3)%block, 1))
1036 : ELSE
1037 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1038 : 1.0_dp, work(1, 1), ldwork, &
1039 : sphi_a(1, sgfa), ldsa, &
1040 : 1.0_dp, op_dip(3)%block(sgfb, sgfa), &
1041 10046 : SIZE(op_dip(3)%block, 1))
1042 : END IF
1043 : END IF ! >= dab
1044 :
1045 : END DO ! jset
1046 :
1047 : END DO ! iset
1048 :
1049 : END DO
1050 74 : CALL neighbor_list_iterator_release(nl_iterator)
1051 :
1052 296 : DO i = 1, 3
1053 296 : NULLIFY (op_dip(i)%block)
1054 : END DO
1055 74 : DEALLOCATE (op_dip)
1056 :
1057 74 : DEALLOCATE (difab, work, basis_set_list)
1058 :
1059 74 : CALL timestop(handle)
1060 :
1061 148 : END SUBROUTINE p_xyz_ao
1062 :
1063 : ! **************************************************************************************************
1064 : !> \brief Calculation of the components of the dipole operator in the length form
1065 : !> by taking the relative position operator r-Rc, with respect a reference point Rc
1066 : !> Probably it does not work for PBC, or maybe yes if the wfn are
1067 : !> sufficiently localized
1068 : !> The elements of the sparse matrices are the integrals in the
1069 : !> basis functions
1070 : !> \param op matrix representation of the p operator
1071 : !> calculated in terms of the contracted basis functions
1072 : !> \param qs_env environment for the lists and the basis sets
1073 : !> \param rc reference vector position
1074 : !> \param order maximum order of the momentum, for the dipole order = 1, order = -2 for quad only
1075 : !> \param minimum_image take into account only the first neighbors in the lists
1076 : !> \param soft ...
1077 : !> \par History
1078 : !> 03.2006 created [MI]
1079 : !> 06.2019 added quarupole only option (A.Bussy)
1080 : !> \author MI
1081 : ! **************************************************************************************************
1082 :
1083 40 : SUBROUTINE rRc_xyz_ao(op, qs_env, rc, order, minimum_image, soft)
1084 :
1085 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op
1086 : TYPE(qs_environment_type), POINTER :: qs_env
1087 : REAL(dp) :: Rc(3)
1088 : INTEGER, INTENT(IN) :: order
1089 : LOGICAL, INTENT(IN), OPTIONAL :: minimum_image, soft
1090 :
1091 : CHARACTER(len=*), PARAMETER :: routineN = 'rRc_xyz_ao'
1092 :
1093 : CHARACTER(LEN=default_string_length) :: basis_type
1094 : INTEGER :: handle, iatom, icol, ikind, imom, inode, irow, iset, jatom, jkind, jset, &
1095 : last_jatom, ldab, ldsa, ldsb, ldwork, M_dim, maxl, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
1096 : sgfb, smom
1097 40 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, npgfa, npgfb, &
1098 40 : nsgfa, nsgfb
1099 40 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1100 : LOGICAL :: found, my_minimum_image, my_soft, &
1101 : new_atom_b
1102 : REAL(KIND=dp) :: dab, Lxo2, Lyo2, Lzo2, rab2
1103 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
1104 40 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
1105 40 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, &
1106 40 : zeta, zetb
1107 40 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab
1108 40 : TYPE(block_p_type), DIMENSION(:), POINTER :: op_dip
1109 : TYPE(cell_type), POINTER :: cell
1110 40 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1111 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1112 : TYPE(neighbor_list_iterator_p_type), &
1113 40 : DIMENSION(:), POINTER :: nl_iterator
1114 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1115 40 : POINTER :: sab_orb
1116 40 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1117 40 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1118 : TYPE(qs_kind_type), POINTER :: qs_kind
1119 :
1120 40 : CALL timeset(routineN, handle)
1121 :
1122 40 : NULLIFY (qs_kind, qs_kind_set)
1123 40 : NULLIFY (cell, particle_set)
1124 40 : NULLIFY (sab_orb)
1125 40 : NULLIFY (mab, op_dip, work)
1126 40 : NULLIFY (la_max, la_min, lb_max, npgfa, npgfb, nsgfa, nsgfb)
1127 40 : NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
1128 :
1129 40 : my_soft = .FALSE.
1130 40 : IF (PRESENT(soft)) my_soft = soft
1131 14 : IF (my_soft) THEN
1132 0 : basis_type = "ORB_SOFT"
1133 : ELSE
1134 40 : basis_type = "ORB"
1135 : END IF
1136 :
1137 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
1138 40 : cell=cell, particle_set=particle_set, sab_orb=sab_orb)
1139 :
1140 40 : nkind = SIZE(qs_kind_set)
1141 :
1142 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1143 40 : maxco=ldwork, maxlgto=maxl)
1144 :
1145 40 : my_minimum_image = .FALSE.
1146 40 : IF (PRESENT(minimum_image)) THEN
1147 38 : my_minimum_image = minimum_image
1148 152 : Lxo2 = SQRT(SUM(cell%hmat(:, 1)**2))/2.0_dp
1149 152 : Lyo2 = SQRT(SUM(cell%hmat(:, 2)**2))/2.0_dp
1150 152 : Lzo2 = SQRT(SUM(cell%hmat(:, 3)**2))/2.0_dp
1151 : END IF
1152 :
1153 40 : ldab = ldwork
1154 :
1155 40 : smom = 1
1156 40 : IF (order == -2) smom = 4
1157 40 : M_dim = ncoset(ABS(order)) - 1
1158 40 : CPASSERT(M_dim <= SIZE(op, 1))
1159 :
1160 200 : ALLOCATE (mab(ldab, ldab, 1:M_dim))
1161 29872 : mab(1:ldab, 1:ldab, 1:M_dim) = 0.0_dp
1162 160 : ALLOCATE (work(ldwork, ldwork))
1163 9944 : work(1:ldwork, 1:ldwork) = 0.0_dp
1164 240 : ALLOCATE (op_dip(smom:M_dim))
1165 :
1166 160 : DO imom = smom, M_dim
1167 160 : NULLIFY (op_dip(imom)%block)
1168 : END DO
1169 :
1170 188 : ALLOCATE (basis_set_list(nkind))
1171 108 : DO ikind = 1, nkind
1172 68 : qs_kind => qs_kind_set(ikind)
1173 68 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1174 108 : IF (ASSOCIATED(basis_set_a)) THEN
1175 68 : basis_set_list(ikind)%gto_basis_set => basis_set_a
1176 : ELSE
1177 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
1178 : END IF
1179 : END DO
1180 40 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1181 157 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1182 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1183 117 : iatom=iatom, jatom=jatom, r=rab)
1184 117 : basis_set_a => basis_set_list(ikind)%gto_basis_set
1185 117 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
1186 117 : basis_set_b => basis_set_list(jkind)%gto_basis_set
1187 117 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
1188 117 : ra = pbc(particle_set(iatom)%r, cell)
1189 : ! basis ikind
1190 117 : first_sgfa => basis_set_a%first_sgf
1191 117 : la_max => basis_set_a%lmax
1192 117 : la_min => basis_set_a%lmin
1193 117 : npgfa => basis_set_a%npgf
1194 117 : nseta = basis_set_a%nset
1195 117 : nsgfa => basis_set_a%nsgf_set
1196 117 : rpgfa => basis_set_a%pgf_radius
1197 117 : set_radius_a => basis_set_a%set_radius
1198 117 : sphi_a => basis_set_a%sphi
1199 117 : zeta => basis_set_a%zet
1200 : ! basis jkind
1201 117 : first_sgfb => basis_set_b%first_sgf
1202 117 : lb_max => basis_set_b%lmax
1203 117 : npgfb => basis_set_b%npgf
1204 117 : nsetb = basis_set_b%nset
1205 117 : nsgfb => basis_set_b%nsgf_set
1206 117 : rpgfb => basis_set_b%pgf_radius
1207 117 : set_radius_b => basis_set_b%set_radius
1208 117 : sphi_b => basis_set_b%sphi
1209 117 : zetb => basis_set_b%zet
1210 :
1211 117 : ldsa = SIZE(sphi_a, 1)
1212 117 : ldsb = SIZE(sphi_b, 1)
1213 117 : IF (inode == 1) last_jatom = 0
1214 :
1215 117 : IF (my_minimum_image) THEN
1216 91 : IF (ABS(rab(1)) > Lxo2 .OR. ABS(rab(2)) > Lyo2 .OR. ABS(rab(3)) > Lzo2) CYCLE
1217 : END IF
1218 :
1219 460 : rb = rab + ra
1220 :
1221 115 : IF (jatom /= last_jatom) THEN
1222 : new_atom_b = .TRUE.
1223 : last_jatom = jatom
1224 : ELSE
1225 : new_atom_b = .FALSE.
1226 : END IF
1227 :
1228 : IF (new_atom_b) THEN
1229 97 : IF (iatom <= jatom) THEN
1230 65 : irow = iatom
1231 65 : icol = jatom
1232 : ELSE
1233 32 : irow = jatom
1234 32 : icol = iatom
1235 : END IF
1236 :
1237 388 : DO imom = smom, M_dim
1238 291 : NULLIFY (op_dip(imom)%block)
1239 : CALL dbcsr_get_block_p(matrix=op(imom)%matrix, &
1240 291 : row=irow, col=icol, block=op_dip(imom)%block, found=found)
1241 388 : CPASSERT(ASSOCIATED(op_dip(imom)%block))
1242 : END DO ! imom
1243 : END IF ! new_atom_b
1244 :
1245 115 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1246 115 : dab = SQRT(rab2)
1247 :
1248 503 : DO iset = 1, nseta
1249 :
1250 348 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1251 348 : sgfa = first_sgfa(1, iset)
1252 :
1253 1621 : DO jset = 1, nsetb
1254 :
1255 1156 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1256 1156 : sgfb = first_sgfb(1, jset)
1257 :
1258 1504 : IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN
1259 :
1260 1036 : rac = pbc(rc, ra, cell)
1261 1036 : rbc = pbc(rc, rb, cell)
1262 :
1263 : ! *** Calculate the primitive overlap integrals ***
1264 : CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
1265 : rpgfa(:, iset), la_min(iset), &
1266 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), &
1267 1036 : ABS(order), rac, rbc, mab)
1268 :
1269 4144 : DO imom = smom, M_dim
1270 : ! *** Contraction ***
1271 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1272 : 1.0_dp, mab(1, 1, imom), ldab, sphi_b(1, sgfb), ldsb, &
1273 3108 : 0.0_dp, work(1, 1), ldwork)
1274 4144 : IF (iatom <= jatom) THEN
1275 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1276 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
1277 : work(1, 1), ldwork, &
1278 : 1.0_dp, op_dip(imom)%block(sgfa, sgfb), &
1279 2385 : SIZE(op_dip(imom)%block, 1))
1280 : ELSE
1281 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1282 : 1.0_dp, work(1, 1), ldwork, &
1283 : sphi_a(1, sgfa), ldsa, &
1284 : 1.0_dp, op_dip(imom)%block(sgfb, sgfa), &
1285 723 : SIZE(op_dip(imom)%block, 1))
1286 : END IF
1287 :
1288 : END DO ! imom
1289 : END IF ! >= dab
1290 :
1291 : END DO ! jset
1292 :
1293 : END DO ! iset
1294 :
1295 : END DO
1296 40 : CALL neighbor_list_iterator_release(nl_iterator)
1297 :
1298 160 : DO imom = smom, M_dim
1299 160 : NULLIFY (op_dip(imom)%block)
1300 : END DO
1301 40 : DEALLOCATE (op_dip)
1302 :
1303 40 : DEALLOCATE (mab, work, basis_set_list)
1304 :
1305 40 : CALL timestop(handle)
1306 :
1307 80 : END SUBROUTINE rRc_xyz_ao
1308 :
1309 : ! **************************************************************************************************
1310 : !> \brief Calculation of the multipole operators integrals
1311 : !> and of its derivatives of the type
1312 : !> [\mu | op | d(\nu)/dR(\nu)]-[d(\mu)/dR(\mu)| op | \nu]
1313 : !> by taking the relative position operator r-Rc, with respect a reference point Rc
1314 : !> The derivative are with respect to the primitive position,
1315 : !> The multipole operator is symmetric and if it does not depend on R(\mu) or R(\nu)
1316 : !> therefore [\mu | op | d(\nu)/dR(\nu)] = -[d(\mu)/dR(\mu)| op | \nu]
1317 : !> [\mu|op|d(\nu)/dR]-[d(\mu)/dR|op|\nu]=2[\mu|op|d(\nu)/dR]
1318 : !> When it is not the case a correction term is needed
1319 : !>
1320 : !> The momentum operator [\mu|M|\nu] is symmetric, the number of components is
1321 : !> determined by the order: 3 for order 1 (x,y,x), 9 for order 2(xx,xy,xz,yy,yz,zz)
1322 : !> The derivative of the type [\mu | op | d(\nu)/dR_i(\nu)], where
1323 : !> i indicates the cartesian direction, is antisymmetric only when
1324 : !> the no component M =(r_i) or (r_i r_j) is in the same cartesian
1325 : !> direction of the derivative, indeed
1326 : !> d([\mu|M|\nu])/dr_i = [d(\mu)/dr_i|M|\nu] + [\mu|M|d(\nu)/dr_i] + [\mu |d(M)/dr_i|\nu]
1327 : !> 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]
1328 : !> Therefore we cannot use an antisymmetric matrix
1329 : !>
1330 : !> The same holds for the derivative with respect to the electronic position r
1331 : !> taking into account that [\mu|op|d(\nu)/dR] = -[\mu|op|d(\nu)/dr]
1332 : !> \param op matrix representation of the p operator
1333 : !> calculated in terms of the contracted basis functions
1334 : !> \param op_der ...
1335 : !> \param qs_env environment for the lists and the basis sets
1336 : !> \param rc reference vector position
1337 : !> \param order maximum order of the momentum, for the dipole order = 1
1338 : !> \param minimum_image take into account only the first neighbors in the lists
1339 : !> \param soft ...
1340 : !> \par History
1341 : !> 03.2006 created [MI]
1342 : !> \author MI
1343 : !> \note
1344 : !> Probably it does not work for PBC, or maybe yes if the wfn are
1345 : !> sufficiently localized
1346 : !> The elements of the sparse matrices are the integrals in the
1347 : !> basis functions
1348 : ! **************************************************************************************************
1349 3750 : SUBROUTINE rRc_xyz_der_ao(op, op_der, qs_env, rc, order, minimum_image, soft)
1350 :
1351 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op
1352 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_der
1353 : TYPE(qs_environment_type), POINTER :: qs_env
1354 : REAL(dp) :: Rc(3)
1355 : INTEGER, INTENT(IN) :: order
1356 : LOGICAL, INTENT(IN), OPTIONAL :: minimum_image, soft
1357 :
1358 : CHARACTER(len=*), PARAMETER :: routineN = 'rRc_xyz_der_ao'
1359 :
1360 : CHARACTER(LEN=default_string_length) :: basis_type
1361 : INTEGER :: handle, i, iatom, icol, idir, ikind, imom, inode, ipgf, irow, iset, j, jatom, &
1362 : jkind, jpgf, jset, last_jatom, lda_min, ldab, ldb_min, ldsa, ldsb, ldwork, M_dim, maxl, &
1363 : na, nb, ncoa, ncob, nda, ndb, nkind, nseta, nsetb, sgfa, sgfb
1364 3750 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1365 3750 : npgfb, nsgfa, nsgfb
1366 3750 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1367 : LOGICAL :: my_minimum_image, my_soft, new_atom_b, &
1368 : op_der_found, op_found
1369 : REAL(KIND=dp) :: alpha, alpha_der, dab, Lxo2, Lyo2, Lzo2, &
1370 : rab2
1371 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
1372 3750 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
1373 3750 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, &
1374 3750 : zeta, zetb
1375 3750 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab, mab_tmp
1376 3750 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: difmab
1377 3750 : TYPE(block_p_type), DIMENSION(:), POINTER :: op_dip
1378 3750 : TYPE(block_p_type), DIMENSION(:, :), POINTER :: op_dip_der
1379 : TYPE(cell_type), POINTER :: cell
1380 3750 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1381 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1382 : TYPE(neighbor_list_iterator_p_type), &
1383 3750 : DIMENSION(:), POINTER :: nl_iterator
1384 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1385 3750 : POINTER :: sab_all
1386 3750 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1387 3750 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1388 : TYPE(qs_kind_type), POINTER :: qs_kind
1389 :
1390 3750 : CALL timeset(routineN, handle)
1391 :
1392 3750 : CPASSERT(ASSOCIATED(op))
1393 3750 : CPASSERT(ASSOCIATED(op_der))
1394 : !IF(.NOT.op_sm_der(1,1)%matrix%symmetry=="none") THEN
1395 3750 : IF (.NOT. dbcsr_get_matrix_type(op_der(1, 1)%matrix) .EQ. dbcsr_type_no_symmetry) THEN
1396 3750 : CPABORT("")
1397 : END IF
1398 :
1399 3750 : NULLIFY (qs_kind, qs_kind_set)
1400 3750 : NULLIFY (cell, particle_set)
1401 3750 : NULLIFY (sab_all)
1402 3750 : NULLIFY (difmab, mab, mab_tmp)
1403 3750 : NULLIFY (op_dip, op_dip_der, work)
1404 3750 : NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb)
1405 3750 : NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
1406 :
1407 3750 : my_soft = .FALSE.
1408 3750 : IF (PRESENT(soft)) my_soft = soft
1409 3750 : IF (my_soft) THEN
1410 2022 : basis_type = "ORB_SOFT"
1411 : ELSE
1412 1728 : basis_type = "ORB"
1413 : END IF
1414 :
1415 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
1416 : cell=cell, particle_set=particle_set, &
1417 3750 : sab_all=sab_all)
1418 :
1419 3750 : nkind = SIZE(qs_kind_set)
1420 :
1421 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1422 3750 : maxco=ldwork, maxlgto=maxl)
1423 :
1424 3750 : my_minimum_image = .FALSE.
1425 3750 : IF (PRESENT(minimum_image)) THEN
1426 3750 : my_minimum_image = minimum_image
1427 15000 : Lxo2 = SQRT(SUM(cell%hmat(:, 1)**2))/2.0_dp
1428 15000 : Lyo2 = SQRT(SUM(cell%hmat(:, 2)**2))/2.0_dp
1429 15000 : Lzo2 = SQRT(SUM(cell%hmat(:, 3)**2))/2.0_dp
1430 : END IF
1431 :
1432 3750 : ldab = ldwork
1433 :
1434 3750 : M_dim = ncoset(order) - 1
1435 3750 : CPASSERT(M_dim <= SIZE(op, 1))
1436 :
1437 18750 : ALLOCATE (mab(ldab, ldab, M_dim))
1438 6053640 : mab(1:ldab, 1:ldab, 1:M_dim) = 0.0_dp
1439 22500 : ALLOCATE (difmab(ldab, ldab, M_dim, 3))
1440 18164670 : difmab(1:ldab, 1:ldab, 1:M_dim, 1:3) = 0.0_dp
1441 :
1442 15000 : ALLOCATE (work(ldwork, ldwork))
1443 672210 : work(1:ldwork, 1:ldwork) = 0.0_dp
1444 45000 : ALLOCATE (op_dip(M_dim))
1445 123750 : ALLOCATE (op_dip_der(M_dim, 3))
1446 :
1447 37500 : DO imom = 1, M_dim
1448 33750 : NULLIFY (op_dip(imom)%block)
1449 138750 : DO i = 1, 3
1450 135000 : NULLIFY (op_dip_der(imom, i)%block)
1451 : END DO
1452 : END DO
1453 :
1454 17388 : ALLOCATE (basis_set_list(nkind))
1455 9888 : DO ikind = 1, nkind
1456 6138 : qs_kind => qs_kind_set(ikind)
1457 6138 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1458 9888 : IF (ASSOCIATED(basis_set_a)) THEN
1459 6138 : basis_set_list(ikind)%gto_basis_set => basis_set_a
1460 : ELSE
1461 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
1462 : END IF
1463 : END DO
1464 3750 : CALL neighbor_list_iterator_create(nl_iterator, sab_all)
1465 284943 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1466 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1467 281193 : iatom=iatom, jatom=jatom, r=rab)
1468 281193 : basis_set_a => basis_set_list(ikind)%gto_basis_set
1469 281193 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
1470 281193 : basis_set_b => basis_set_list(jkind)%gto_basis_set
1471 281193 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
1472 281193 : ra = pbc(particle_set(iatom)%r, cell)
1473 : ! basis ikind
1474 281193 : first_sgfa => basis_set_a%first_sgf
1475 281193 : la_max => basis_set_a%lmax
1476 281193 : la_min => basis_set_a%lmin
1477 281193 : npgfa => basis_set_a%npgf
1478 281193 : nseta = basis_set_a%nset
1479 281193 : nsgfa => basis_set_a%nsgf_set
1480 281193 : rpgfa => basis_set_a%pgf_radius
1481 281193 : set_radius_a => basis_set_a%set_radius
1482 281193 : sphi_a => basis_set_a%sphi
1483 281193 : zeta => basis_set_a%zet
1484 : ! basis jkind
1485 281193 : first_sgfb => basis_set_b%first_sgf
1486 281193 : lb_max => basis_set_b%lmax
1487 281193 : lb_min => basis_set_b%lmin
1488 281193 : npgfb => basis_set_b%npgf
1489 281193 : nsetb = basis_set_b%nset
1490 281193 : nsgfb => basis_set_b%nsgf_set
1491 281193 : rpgfb => basis_set_b%pgf_radius
1492 281193 : set_radius_b => basis_set_b%set_radius
1493 281193 : sphi_b => basis_set_b%sphi
1494 281193 : zetb => basis_set_b%zet
1495 :
1496 281193 : ldsa = SIZE(sphi_a, 1)
1497 281193 : IF (ldsa .EQ. 0) CYCLE
1498 281172 : ldsb = SIZE(sphi_b, 1)
1499 281172 : IF (ldsb .EQ. 0) CYCLE
1500 281172 : IF (inode == 1) last_jatom = 0
1501 :
1502 281172 : IF (my_minimum_image) THEN
1503 0 : IF (ABS(rab(1)) > Lxo2 .OR. ABS(rab(2)) > Lyo2 .OR. ABS(rab(3)) > Lzo2) CYCLE
1504 : END IF
1505 :
1506 1124688 : rb = rab + ra
1507 :
1508 281172 : IF (jatom /= last_jatom) THEN
1509 : new_atom_b = .TRUE.
1510 : last_jatom = jatom
1511 : ELSE
1512 : new_atom_b = .FALSE.
1513 : END IF
1514 :
1515 : IF (new_atom_b) THEN
1516 17958 : irow = iatom
1517 17958 : icol = jatom
1518 17958 : alpha_der = 2.0_dp
1519 :
1520 179580 : DO imom = 1, M_dim
1521 161622 : NULLIFY (op_dip(imom)%block)
1522 : CALL dbcsr_get_block_p(matrix=op(imom)%matrix, &
1523 : row=irow, col=icol, &
1524 : block=op_dip(imom)%block, &
1525 161622 : found=op_found)
1526 161622 : CPASSERT(op_found .AND. ASSOCIATED(op_dip(imom)%block))
1527 826068 : DO idir = 1, 3
1528 484866 : NULLIFY (op_dip_der(imom, idir)%block)
1529 : CALL dbcsr_get_block_p(matrix=op_der(imom, idir)%matrix, &
1530 : row=irow, col=icol, &
1531 : block=op_dip_der(imom, idir)%block, &
1532 484866 : found=op_der_found)
1533 646488 : CPASSERT(op_der_found .AND. ASSOCIATED(op_dip_der(imom, idir)%block))
1534 : END DO ! idir
1535 : END DO ! imom
1536 : END IF ! new_atom_b
1537 :
1538 281172 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1539 281172 : dab = SQRT(rab2)
1540 :
1541 936090 : DO iset = 1, nseta
1542 :
1543 651168 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1544 651168 : sgfa = first_sgfa(1, iset)
1545 :
1546 2511669 : DO jset = 1, nsetb
1547 :
1548 1579308 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1549 1579308 : sgfb = first_sgfb(1, jset)
1550 :
1551 2230476 : IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN
1552 :
1553 596982 : rac = pbc(rc, ra, cell)
1554 2387928 : rbc = rac + rab
1555 : ! rac = pbc(rc,ra,cell)
1556 : ! rbc = pbc(rc,rb,cell)
1557 :
1558 : ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
1559 2958078 : npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(order) - 1))
1560 :
1561 596982 : lda_min = MAX(0, la_min(iset) - 1)
1562 596982 : ldb_min = MAX(0, lb_min(jset) - 1)
1563 : ! *** Calculate the primitive overlap integrals ***
1564 : CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), &
1565 : rpgfa(:, iset), lda_min, &
1566 : lb_max(jset) + 1, npgfb(jset), zetb(:, jset), rpgfb(:, jset), &
1567 596982 : order, rac, rbc, mab_tmp)
1568 :
1569 : ! *** Calculate the derivatives
1570 : CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), &
1571 : rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
1572 : zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
1573 596982 : difmab, mab_ext=mab_tmp)
1574 :
1575 : ! Contract and copy in the sparse matrix
1576 874316388 : mab = 0.0_dp
1577 5969820 : DO imom = 1, M_dim
1578 5372838 : na = 0
1579 5372838 : nda = 0
1580 17098830 : DO ipgf = 1, npgfa(iset)
1581 11725992 : nb = 0
1582 11725992 : ndb = 0
1583 42864039 : DO jpgf = 1, npgfb(jset)
1584 127844055 : DO j = 1, ncoset(lb_max(jset))
1585 486284148 : DO i = 1, ncoset(la_max(iset))
1586 455146101 : mab(i + na, j + nb, imom) = mab_tmp(i + nda, j + ndb, imom)
1587 : END DO ! i
1588 : END DO ! j
1589 31138047 : nb = nb + ncoset(lb_max(jset))
1590 42864039 : ndb = ndb + ncoset(lb_max(jset) + 1)
1591 : END DO ! jpgf
1592 11725992 : na = na + ncoset(la_max(iset))
1593 17098830 : nda = nda + ncoset(la_max(iset) + 1)
1594 : END DO ! ipgf
1595 :
1596 : ! *** Contraction ***
1597 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1598 : 1.0_dp, mab(1, 1, imom), ldab, sphi_b(1, sgfb), ldsb, &
1599 5372838 : 0.0_dp, work(1, 1), ldwork)
1600 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1601 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
1602 : work(1, 1), ldwork, &
1603 : 1.0_dp, op_dip(imom)%block(sgfa, sgfb), &
1604 5372838 : SIZE(op_dip(imom)%block, 1))
1605 :
1606 5372838 : alpha = -1.0_dp !-alpha_der
1607 22088334 : DO idir = 1, 3
1608 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1609 : alpha, difmab(1, 1, imom, idir), ldab, sphi_b(1, sgfb), ldsb, &
1610 16118514 : 0.0_dp, work(1, 1), ldwork)
1611 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1612 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
1613 : work(1, 1), ldwork, &
1614 : 1.0_dp, op_dip_der(imom, idir)%block(sgfa, sgfb), &
1615 21491352 : SIZE(op_dip_der(imom, idir)%block, 1))
1616 :
1617 : END DO ! idir
1618 :
1619 : END DO ! imom
1620 :
1621 596982 : DEALLOCATE (mab_tmp)
1622 : END IF ! >= dab
1623 :
1624 : END DO ! jset
1625 :
1626 : END DO ! iset
1627 :
1628 : END DO
1629 3750 : CALL neighbor_list_iterator_release(nl_iterator)
1630 :
1631 15000 : DO i = 1, 3
1632 15000 : NULLIFY (op_dip(i)%block)
1633 : END DO
1634 3750 : DEALLOCATE (op_dip, op_dip_der)
1635 :
1636 3750 : DEALLOCATE (mab, difmab, work, basis_set_list)
1637 :
1638 3750 : CALL timestop(handle)
1639 :
1640 7500 : END SUBROUTINE rRc_xyz_der_ao
1641 :
1642 : END MODULE qs_operators_ao
1643 :
|
