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