(git:0e87843)
Loading...
Searching...
No Matches
libgrpp_integrals.F
Go to the documentation of this file.
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!> \brief Local and semi-local ECP integrals using the libgrpp library
10! **************************************************************************************************
11
12MODULE libgrpp_integrals
13 USE ai_cartesian_operators, ONLY: apply_cartesian_operator
14 USE kinds, ONLY: dp
15 USE libgrpp, ONLY: libgrpp_init,&
16 libgrpp_type1_integrals,&
17 libgrpp_type1_integrals_gradient,&
18 libgrpp_type2_integrals,&
19 libgrpp_type2_integrals_gradient
20 USE mathconstants, ONLY: pi
21 USE orbital_pointers, ONLY: nco,&
22 ncoset
23#include "./base/base_uses.f90"
24
25 IMPLICIT NONE
26 PRIVATE
27
28 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libgrpp_integrals'
29
30 PUBLIC :: libgrpp_semilocal_integrals, libgrpp_local_integrals, &
31 libgrpp_local_forces_ref, libgrpp_semilocal_forces_ref
32
33CONTAINS
34
35! **************************************************************************************************
36!> \brief Local ECP integrals using libgrpp
37!> \param la_max_set ...
38!> \param la_min_set ...
39!> \param npgfa ...
40!> \param rpgfa ...
41!> \param zeta ...
42!> \param lb_max_set ...
43!> \param lb_min_set ...
44!> \param npgfb ...
45!> \param rpgfb ...
46!> \param zetb ...
47!> \param npot_ecp ...
48!> \param alpha_ecp ...
49!> \param coeffs_ecp ...
50!> \param nrpot_ecp ...
51!> \param rpgfc ...
52!> \param rab ...
53!> \param dab ...
54!> \param rac ...
55!> \param dac ...
56!> \param dbc ...
57!> \param vab ...
58!> \param pab ...
59!> \param force_a ...
60!> \param force_b ...
61! **************************************************************************************************
62 SUBROUTINE libgrpp_local_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
63 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
64 npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
65 rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
66
67 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
68 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
69 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
70 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
71 INTEGER, INTENT(IN) :: npot_ecp
72 REAL(kind=dp), DIMENSION(1:npot_ecp), INTENT(IN) :: alpha_ecp, coeffs_ecp
73 INTEGER, DIMENSION(1:npot_ecp), INTENT(IN) :: nrpot_ecp
74 REAL(kind=dp), INTENT(IN) :: rpgfc
75 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
76 REAL(kind=dp), INTENT(IN) :: dab
77 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
78 REAL(kind=dp), INTENT(IN) :: dac, dbc
79 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
80 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
81 OPTIONAL :: pab
82 REAL(kind=dp), DIMENSION(3), INTENT(INOUT), &
83 OPTIONAL :: force_a, force_b
84
85 INTEGER :: a_offset, a_start, b_offset, b_start, i, &
86 ipgf, j, jpgf, li, lj, ncoa, ncob
87 LOGICAL :: calc_forces
88 REAL(dp) :: expi, expj, fac_a, fac_b, mindist, &
89 normi, normj, prefi, prefj, zeti, zetj
90 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmpx, tmpy, tmpz
91 REAL(dp), DIMENSION(3) :: ra, rb, rc
92
93 CALL libgrpp_init()
94
95 calc_forces = .false.
96 IF (PRESENT(pab) .AND. PRESENT(force_a) .AND. PRESENT(force_b)) calc_forces = .true.
97
98 IF (calc_forces) THEN
99
100 !Note: warning against numerical stability of libgrpp gradients. The day the library becomes
101 ! stable, this routine can be used immediatly as is, and the warning removed.
102 CALL cp_warn(__location__, &
103 "ECP gradients calculated with the libgrpp library are, to this date, not numerically stable. "// &
104 "Please use the reference routine 'libgrpp_local_forces_ref' instead.")
105
106 !there is a weird feature of libgrpp gradients, which is such that the gradient is calculated
107 !for a point in space, and not with respect to an atomic center. For example, if atoms A and
108 !B are the same (and C is different), then d<A | U_C | B>/dPx = d<A | U_C | B>/dAx + d<A | U_C | B>/dBx
109 !Because we want the forces on centers A and B seprately, we need a case study on atomic positions
110 !We always calculate the gradient wrt to atomic position of A and B, and we scale accordingly
111
112 mindist = 1.0e-6_dp
113 !If ra != rb != rc
114 IF (dab >= mindist .AND. dbc >= mindist .AND. dac >= mindist) THEN
115 fac_a = 1.0_dp
116 fac_b = 1.0_dp
117
118 !If ra = rb, but ra != rc
119 ELSE IF (dab < mindist .AND. dac >= mindist) THEN
120 fac_a = 0.5_dp
121 fac_b = 0.5_dp
122
123 !IF ra != rb but ra = rc
124 ELSE IF (dab >= mindist .AND. dac < mindist) THEN
125 fac_a = 0.5_dp
126 fac_b = 1.0_dp
127
128 !IF ra != rb but rb = rc
129 ELSE IF (dab >= mindist .AND. dbc < mindist) THEN
130 fac_a = 1.0_dp
131 fac_b = 0.5_dp
132
133 !If all atoms the same --> no force
134 ELSE
135 calc_forces = .false.
136 END IF
137 END IF
138
139 !libgrpp requires absolute positions, not relative ones
140 ra(:) = 0.0_dp
141 rb(:) = rab(:)
142 rc(:) = rac(:)
143
144 ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
145 IF (calc_forces) THEN
146 ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
147 ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
148 ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
149 END IF
150
151 DO ipgf = 1, npgfa
152 IF (rpgfa(ipgf) + rpgfc < dac) cycle
153 zeti = zeta(ipgf)
154 a_start = (ipgf - 1)*ncoset(la_max_set)
155
156 DO jpgf = 1, npgfb
157 IF (rpgfb(jpgf) + rpgfc < dbc) cycle
158 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
159 zetj = zetb(jpgf)
160 b_start = (jpgf - 1)*ncoset(lb_max_set)
161
162 DO li = la_min_set, la_max_set
163 a_offset = a_start + ncoset(li - 1)
164 ncoa = nco(li)
165 prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
166 expi = 0.25_dp*real(2*li + 3, dp)
167 normi = 1.0_dp/(prefi*zeti**expi)
168
169 DO lj = lb_min_set, lb_max_set
170 b_offset = b_start + ncoset(lj - 1)
171 ncob = nco(lj)
172 prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
173 expj = 0.25_dp*real(2*lj + 3, dp)
174 normj = 1.0_dp/(prefj*zetj**expj)
175
176 tmp(1:ncoa*ncob) = 0.0_dp
177 !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
178 !the 1/norm coefficients for PGFi and PGFj
179 CALL libgrpp_type1_integrals(ra, li, 1, [normi], [zeti], &
180 rb, lj, 1, [normj], [zetj], &
181 rc, [npot_ecp], nrpot_ecp, &
182 coeffs_ecp, alpha_ecp, tmp)
183
184 !note: tmp array is in C row-major ordering
185 DO j = 1, ncob
186 DO i = 1, ncoa
187 vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
188 END DO
189 END DO
190
191 IF (calc_forces) THEN
192 tmpx(1:ncoa*ncob) = 0.0_dp
193 tmpy(1:ncoa*ncob) = 0.0_dp
194 tmpz(1:ncoa*ncob) = 0.0_dp
195
196 !force wrt to atomic position A
197 CALL libgrpp_type1_integrals_gradient(ra, li, 1, [normi], [zeti], &
198 rb, lj, 1, [normj], [zetj], &
199 rc, [npot_ecp], nrpot_ecp, &
200 coeffs_ecp, alpha_ecp, ra, &
201 tmpx, tmpy, tmpz)
202
203 !note: tmp array is in C row-major ordering
204 !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
205 DO j = 1, ncob
206 DO i = 1, ncoa
207 force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
208 force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
209 force_a(3) = force_a(3) + fac_a*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
210 END DO
211 END DO
212
213 tmpx(1:ncoa*ncob) = 0.0_dp
214 tmpy(1:ncoa*ncob) = 0.0_dp
215 tmpz(1:ncoa*ncob) = 0.0_dp
216
217 !force wrt to atomic position B
218 CALL libgrpp_type1_integrals_gradient(ra, li, 1, [normi], [zeti], &
219 rb, lj, 1, [normj], [zetj], &
220 rc, [npot_ecp], nrpot_ecp, &
221 coeffs_ecp, alpha_ecp, rb, &
222 tmpx, tmpy, tmpz)
223
224 !note: tmp array is in C row-major ordering
225 !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
226 DO j = 1, ncob
227 DO i = 1, ncoa
228 force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
229 force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
230 force_b(3) = force_b(3) + fac_b*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
231 END DO
232 END DO
233 END IF
234
235 END DO !lj
236 END DO !li
237
238 END DO !jpgf
239 END DO !ipgf
240 END SUBROUTINE libgrpp_local_integrals
241
242! **************************************************************************************************
243!> \brief Semi-local ECP integrals using libgrpp.
244!> \param la_max_set ...
245!> \param la_min_set ...
246!> \param npgfa ...
247!> \param rpgfa ...
248!> \param zeta ...
249!> \param lb_max_set ...
250!> \param lb_min_set ...
251!> \param npgfb ...
252!> \param rpgfb ...
253!> \param zetb ...
254!> \param lmax_ecp ...
255!> \param npot_ecp ...
256!> \param alpha_ecp ...
257!> \param coeffs_ecp ...
258!> \param nrpot_ecp ...
259!> \param rpgfc ...
260!> \param rab ...
261!> \param dab ...
262!> \param rac ...
263!> \param dac ...
264!> \param dbc ...
265!> \param vab ...
266!> \param pab ...
267!> \param force_a ...
268!> \param force_b ...
269! **************************************************************************************************
270 SUBROUTINE libgrpp_semilocal_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
271 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
272 lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
273 rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
274
275 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
276 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
277 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
278 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
279 INTEGER, INTENT(IN) :: lmax_ecp
280 INTEGER, DIMENSION(0:10), INTENT(IN) :: npot_ecp
281 REAL(kind=dp), DIMENSION(1:15, 0:10), INTENT(IN) :: alpha_ecp, coeffs_ecp
282 INTEGER, DIMENSION(1:15, 0:10), INTENT(IN) :: nrpot_ecp
283 REAL(kind=dp), INTENT(IN) :: rpgfc
284 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
285 REAL(kind=dp), INTENT(IN) :: dab
286 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
287 REAL(kind=dp), INTENT(IN) :: dac, dbc
288 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
289 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
290 OPTIONAL :: pab
291 REAL(kind=dp), DIMENSION(3), INTENT(INOUT), &
292 OPTIONAL :: force_a, force_b
293
294 INTEGER :: a_offset, a_start, b_offset, b_start, i, &
295 ipgf, j, jpgf, li, lj, lk, ncoa, ncob
296 LOGICAL :: calc_forces
297 REAL(dp) :: expi, expj, fac_a, fac_b, mindist, &
298 normi, normj, prefi, prefj, zeti, zetj
299 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmpx, tmpy, tmpz
300 REAL(dp), DIMENSION(3) :: ra, rb, rc
301
302 CALL libgrpp_init()
303
304 calc_forces = .false.
305 IF (PRESENT(pab) .AND. PRESENT(force_a) .AND. PRESENT(force_b)) calc_forces = .true.
306
307 IF (calc_forces) THEN
308
309 !Note: warning against numerical stability of libgrpp gradients. The day the library becomes
310 ! stable, this routine can be used immediatly as is, and the warning removed.
311 CALL cp_warn(__location__, &
312 "ECP gradients calculated with the libgrpp library are, to this date, not numerically stable. "// &
313 "Please use the reference routine 'libgrpp_semilocal_forces_ref' instead.")
314
315 !there is a weird feature of libgrpp gradients, which is such that the gradient is calculated
316 !for a point in space, and not with respect to an atomic center. For example, if atoms A and
317 !B are the same (and C is different), then d<A | U_C | B>/dPx = d<A | U_C | B>/dAx + d<A | U_C | B>/dBx
318 !Because we want the forces on centers A and B seprately, we need a case study on atomic positions
319 !We always calculate the gradient wrt to atomic position of A and B, and we scale accordingly
320
321 mindist = 1.0e-6_dp
322 !If ra != rb != rc
323 IF (dab >= mindist .AND. dbc >= mindist .AND. dac >= mindist) THEN
324 fac_a = 1.0_dp
325 fac_b = 1.0_dp
326
327 !If ra = rb, but ra != rc
328 ELSE IF (dab < mindist .AND. dac >= mindist) THEN
329 fac_a = 0.5_dp
330 fac_b = 0.5_dp
331
332 !IF ra != rb but ra = rc
333 ELSE IF (dab >= mindist .AND. dac < mindist) THEN
334 fac_a = 0.5_dp
335 fac_b = 1.0_dp
336
337 !IF ra != rb but rb = rc
338 ELSE IF (dab >= mindist .AND. dbc < mindist) THEN
339 fac_a = 1.0_dp
340 fac_b = 0.5_dp
341
342 !If all atoms the same --> no force
343 ELSE
344 calc_forces = .false.
345 END IF
346 END IF
347
348 !libgrpp requires absolute positions, not relative ones
349 ra(:) = 0.0_dp
350 rb(:) = rab(:)
351 rc(:) = rac(:)
352
353 ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
354 IF (calc_forces) THEN
355 ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
356 ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
357 ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
358 END IF
359
360 DO ipgf = 1, npgfa
361 IF (rpgfa(ipgf) + rpgfc < dac) cycle
362 zeti = zeta(ipgf)
363 a_start = (ipgf - 1)*ncoset(la_max_set)
364
365 DO jpgf = 1, npgfb
366 IF (rpgfb(jpgf) + rpgfc < dbc) cycle
367 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
368 zetj = zetb(jpgf)
369 b_start = (jpgf - 1)*ncoset(lb_max_set)
370
371 DO li = la_min_set, la_max_set
372 a_offset = a_start + ncoset(li - 1)
373 ncoa = nco(li)
374 prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
375 expi = 0.25_dp*real(2*li + 3, dp)
376 normi = 1.0_dp/(prefi*zeti**expi)
377
378 DO lj = lb_min_set, lb_max_set
379 b_offset = b_start + ncoset(lj - 1)
380 ncob = nco(lj)
381 prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
382 expj = 0.25_dp*real(2*lj + 3, dp)
383 normj = 1.0_dp/(prefj*zetj**expj)
384
385 !Loop over ECP angular momentum
386 DO lk = 0, lmax_ecp
387 tmp(1:ncoa*ncob) = 0.0_dp
388 !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
389 !the 1/norm coefficients for PGFi and PGFj
390 CALL libgrpp_type2_integrals(ra, li, 1, [normi], [zeti], &
391 rb, lj, 1, [normj], [zetj], &
392 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
393 coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
394
395 !note: tmp array is in C row-major ordering
396 DO j = 1, ncob
397 DO i = 1, ncoa
398 vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
399 END DO
400 END DO
401
402 IF (calc_forces) THEN
403
404 tmpx(1:ncoa*ncob) = 0.0_dp
405 tmpy(1:ncoa*ncob) = 0.0_dp
406 tmpz(1:ncoa*ncob) = 0.0_dp
407
408 !force wrt to atomic position A
409 CALL libgrpp_type2_integrals_gradient(ra, li, 1, [normi], [zeti], &
410 rb, lj, 1, [normj], [zetj], &
411 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
412 coeffs_ecp(:, lk), alpha_ecp(:, lk), ra, &
413 tmpx, tmpy, tmpz)
414
415 !note: tmp array is in C row-major ordering
416 !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
417 DO j = 1, ncob
418 DO i = 1, ncoa
419 force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
420 force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
421 force_a(3) = force_a(3) + fac_a*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
422 END DO
423 END DO
424
425 tmpx(1:ncoa*ncob) = 0.0_dp
426 tmpy(1:ncoa*ncob) = 0.0_dp
427 tmpz(1:ncoa*ncob) = 0.0_dp
428
429 !force wrt to atomic position B
430 CALL libgrpp_type2_integrals_gradient(ra, li, 1, [normi], [zeti], &
431 rb, lj, 1, [normj], [zetj], &
432 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
433 coeffs_ecp(:, lk), alpha_ecp(:, lk), rb, &
434 tmpx, tmpy, tmpz)
435 !note: tmp array is in C row-major ordering
436 !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
437 DO j = 1, ncob
438 DO i = 1, ncoa
439 force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
440 force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
441 force_b(3) = force_b(3) + fac_b*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
442 END DO
443 END DO
444
445 END IF !calc_forces
446
447 END DO !lk
448
449 END DO !lj
450 END DO !li
451
452 END DO !jpgf
453 END DO !ipgf
454 END SUBROUTINE libgrpp_semilocal_integrals
455
456! **************************************************************************************************
457!> \brief Reference local ECP force routine using l+-1 integrals. No call is made to the numerically
458!> unstable gradient routine of libgrpp. Calculates both the integrals and the forces.
459!> \param la_max_set ...
460!> \param la_min_set ...
461!> \param npgfa ...
462!> \param rpgfa ...
463!> \param zeta ...
464!> \param lb_max_set ...
465!> \param lb_min_set ...
466!> \param npgfb ...
467!> \param rpgfb ...
468!> \param zetb ...
469!> \param npot_ecp ...
470!> \param alpha_ecp ...
471!> \param coeffs_ecp ...
472!> \param nrpot_ecp ...
473!> \param rpgfc ...
474!> \param rab ...
475!> \param dab ...
476!> \param rac ...
477!> \param dac ...
478!> \param dbc ...
479!> \param vab ...
480!> \param pab ...
481!> \param force_a ...
482!> \param force_b ...
483!> \note: this is a reference routine, which has no reason to be used once the libgrpp gradients
484!> become numerically stable
485! **************************************************************************************************
486 SUBROUTINE libgrpp_local_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
487 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
488 npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
489 rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
490
491 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
492 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
493 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
494 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
495 INTEGER, INTENT(IN) :: npot_ecp
496 REAL(kind=dp), DIMENSION(1:npot_ecp), INTENT(IN) :: alpha_ecp, coeffs_ecp
497 INTEGER, DIMENSION(1:npot_ecp), INTENT(IN) :: nrpot_ecp
498 REAL(kind=dp), INTENT(IN) :: rpgfc
499 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
500 REAL(kind=dp), INTENT(IN) :: dab
501 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
502 REAL(kind=dp), INTENT(IN) :: dac, dbc
503 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
504 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pab
505 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: force_a, force_b
506
507 INTEGER :: a_offset, a_offset_f, a_start, &
508 a_start_f, b_offset, b_offset_f, &
509 b_start, b_start_f, i, ipgf, j, jpgf, &
510 li, lj, ncoa, ncob
511 REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
512 zeti, zetj
513 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
514 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vab_f
515 REAL(dp), DIMENSION(3) :: ra, rb, rc
516
517 CALL libgrpp_init()
518
519 !Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
520 ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
521 vab_f(:, :) = 0.0_dp
522
523 !libgrpp requires absolute positions, not relative ones
524 ra(:) = 0.0_dp
525 rb(:) = rab(:)
526 rc(:) = rac(:)
527
528 ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
529
530 DO ipgf = 1, npgfa
531 IF (rpgfa(ipgf) + rpgfc < dac) cycle
532 zeti = zeta(ipgf)
533 a_start = (ipgf - 1)*ncoset(la_max_set)
534 a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
535
536 DO jpgf = 1, npgfb
537 IF (rpgfb(jpgf) + rpgfc < dbc) cycle
538 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
539 zetj = zetb(jpgf)
540 b_start = (jpgf - 1)*ncoset(lb_max_set)
541 b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
542
543 DO li = max(0, la_min_set - 1), la_max_set + 1
544 a_offset = a_start + ncoset(li - 1)
545 a_offset_f = a_start_f + ncoset(li - 1)
546 ncoa = nco(li)
547 prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
548 expi = 0.25_dp*real(2*li + 3, dp)
549 normi = 1.0_dp/(prefi*zeti**expi)
550
551 DO lj = max(0, lb_min_set - 1), lb_max_set + 1
552 b_offset = b_start + ncoset(lj - 1)
553 b_offset_f = b_start_f + ncoset(lj - 1)
554 ncob = nco(lj)
555 prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
556 expj = 0.25_dp*real(2*lj + 3, dp)
557 normj = 1.0_dp/(prefj*zetj**expj)
558
559 tmp(1:ncoa*ncob) = 0.0_dp
560 !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
561 !the 1/norm coefficients for PGFi and PGFj
562 CALL libgrpp_type1_integrals(ra, li, 1, [normi], [zeti], &
563 rb, lj, 1, [normj], [zetj], &
564 rc, [npot_ecp], nrpot_ecp, &
565 coeffs_ecp, alpha_ecp, tmp)
566
567 !the l+-1 integrals for gradient calculation
568 DO j = 1, ncob
569 DO i = 1, ncoa
570 vab_f(a_offset_f + i, b_offset_f + j) = &
571 vab_f(a_offset_f + i, b_offset_f + j) + tmp((i - 1)*ncob + j)
572 END DO
573 END DO
574
575 !the actual integrals
576 IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
577 DO j = 1, ncob
578 DO i = 1, ncoa
579 vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
580 END DO
581 END DO
582 END IF
583
584 END DO !lj
585 END DO !li
586
587 END DO !jpgf
588 END DO !ipgf
589
590 CALL contract_basis_center_forces(la_max_set, npgfa, zeta, lb_max_set, npgfb, zetb, &
591 vab_f, pab, force_a, force_b)
592 END SUBROUTINE libgrpp_local_forces_ref
593
594! **************************************************************************************************
595!> \brief Reference semi-local ECP forces using l+-1 integrals. No call is made to the numerically
596!> unstable gradient routine of libgrpp. Calculates both the integrals and the forces.
597!> \param la_max_set ...
598!> \param la_min_set ...
599!> \param npgfa ...
600!> \param rpgfa ...
601!> \param zeta ...
602!> \param lb_max_set ...
603!> \param lb_min_set ...
604!> \param npgfb ...
605!> \param rpgfb ...
606!> \param zetb ...
607!> \param lmax_ecp ...
608!> \param npot_ecp ...
609!> \param alpha_ecp ...
610!> \param coeffs_ecp ...
611!> \param nrpot_ecp ...
612!> \param rpgfc ...
613!> \param rab ...
614!> \param dab ...
615!> \param rac ...
616!> \param dac ...
617!> \param dbc ...
618!> \param vab ...
619!> \param pab ...
620!> \param force_a ...
621!> \param force_b ...
622!> \note: this is a reference routine, which has no reason to be used once the libgrpp gradients
623!> become numerically stable
624! **************************************************************************************************
625 SUBROUTINE libgrpp_semilocal_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
626 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
627 lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
628 rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
629
630 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
631 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
632 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
633 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
634 INTEGER, INTENT(IN) :: lmax_ecp
635 INTEGER, DIMENSION(0:10), INTENT(IN) :: npot_ecp
636 REAL(kind=dp), DIMENSION(1:15, 0:10), INTENT(IN) :: alpha_ecp, coeffs_ecp
637 INTEGER, DIMENSION(1:15, 0:10), INTENT(IN) :: nrpot_ecp
638 REAL(kind=dp), INTENT(IN) :: rpgfc
639 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
640 REAL(kind=dp), INTENT(IN) :: dab
641 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
642 REAL(kind=dp), INTENT(IN) :: dac, dbc
643 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
644 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pab
645 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: force_a, force_b
646
647 INTEGER :: a_offset, a_offset_f, a_start, &
648 a_start_f, b_offset, b_offset_f, &
649 b_start, b_start_f, i, ipgf, j, jpgf, &
650 li, lj, lk, ncoa, ncob
651 REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
652 zeti, zetj
653 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
654 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vab_f
655 REAL(dp), DIMENSION(3) :: ra, rb, rc
656
657 CALL libgrpp_init()
658
659 !Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
660 ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
661 vab_f(:, :) = 0.0_dp
662
663 !libgrpp requires absolute positions, not relative ones
664 ra(:) = 0.0_dp
665 rb(:) = rab(:)
666 rc(:) = rac(:)
667
668 ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
669
670 DO ipgf = 1, npgfa
671 IF (rpgfa(ipgf) + rpgfc < dac) cycle
672 zeti = zeta(ipgf)
673 a_start = (ipgf - 1)*ncoset(la_max_set)
674 a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
675
676 DO jpgf = 1, npgfb
677 IF (rpgfb(jpgf) + rpgfc < dbc) cycle
678 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
679 zetj = zetb(jpgf)
680 b_start = (jpgf - 1)*ncoset(lb_max_set)
681 b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
682
683 DO li = max(0, la_min_set - 1), la_max_set + 1
684 a_offset = a_start + ncoset(li - 1)
685 a_offset_f = a_start_f + ncoset(li - 1)
686 ncoa = nco(li)
687 prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
688 expi = 0.25_dp*real(2*li + 3, dp)
689 normi = 1.0_dp/(prefi*zeti**expi)
690
691 DO lj = max(0, lb_min_set - 1), lb_max_set + 1
692 b_offset = b_start + ncoset(lj - 1)
693 b_offset_f = b_start_f + ncoset(lj - 1)
694 ncob = nco(lj)
695 prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
696 expj = 0.25_dp*real(2*lj + 3, dp)
697 normj = 1.0_dp/(prefj*zetj**expj)
698
699 !Loop over ECP angular momentum
700 DO lk = 0, lmax_ecp
701 tmp(1:ncoa*ncob) = 0.0_dp
702 !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
703 !the 1/norm coefficients for PGFi and PGFj
704 CALL libgrpp_type2_integrals(ra, li, 1, [normi], [zeti], &
705 rb, lj, 1, [normj], [zetj], &
706 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
707 coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
708
709 !the l+-1 integrals for gradient calculation
710 DO j = 1, ncob
711 DO i = 1, ncoa
712 vab_f(a_offset_f + i, b_offset_f + j) = &
713 vab_f(a_offset_f + i, b_offset_f + j) + tmp((i - 1)*ncob + j)
714 END DO
715 END DO
716
717 !the actual integrals
718 IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
719 DO j = 1, ncob
720 DO i = 1, ncoa
721 vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
722 END DO
723 END DO
724 END IF
725
726 END DO !lk
727
728 END DO !lj
729 END DO !li
730
731 END DO !jpgf
732 END DO !ipgf
733
734 CALL contract_basis_center_forces(la_max_set, npgfa, zeta, lb_max_set, npgfb, zetb, &
735 vab_f, pab, force_a, force_b)
736 END SUBROUTINE libgrpp_semilocal_forces_ref
737
738! **************************************************************************************************
739!> \brief Contract the bra and ket basis-center derivatives of precomputed ECP integrals.
740!> \param la_max highest bra shell
741!> \param npgfa number of bra primitives
742!> \param zeta bra Gaussian exponents
743!> \param lb_max highest ket shell
744!> \param npgfb number of ket primitives
745!> \param zetb ket Gaussian exponents
746!> \param vab primitive integrals including raised shells and zero screened blocks
747!> \param pab primitive density block
748!> \param force_a accumulated bra-center force
749!> \param force_b accumulated ket-center force
750! **************************************************************************************************
751 SUBROUTINE contract_basis_center_forces(la_max, npgfa, zeta, lb_max, npgfb, zetb, &
752 vab, pab, force_a, force_b)
753 INTEGER, INTENT(IN) :: la_max, npgfa
754 REAL(kind=dp), INTENT(IN) :: zeta(:)
755 INTEGER, INTENT(IN) :: lb_max, npgfb
756 REAL(kind=dp), INTENT(IN) :: zetb(:), vab(:, :), pab(:, :)
757 REAL(kind=dp), INTENT(INOUT) :: force_a(3), force_b(3)
758
759 INTEGER :: i, j, na, nb
760 REAL(kind=dp), ALLOCATABLE :: deriv_ab(:, :, :)
761
762 na = npgfa*ncoset(la_max)
763 nb = npgfb*ncoset(lb_max)
764 ALLOCATE (deriv_ab(na, nb, 3))
765 CALL apply_cartesian_operator(la_max, 0, npgfa, lb_max, 0, npgfb, vab, deriv_ab, &
766 right=.false., derivative=.true., exponents=zeta)
767 DO j = 1, nb
768 DO i = 1, na
769 force_a = force_a + deriv_ab(i, j, :)*pab(i, j)
770 END DO
771 END DO
772 CALL apply_cartesian_operator(la_max, 0, npgfa, lb_max, 0, npgfb, vab, deriv_ab, &
773 right=.true., derivative=.true., exponents=zetb)
774 DO j = 1, nb
775 DO i = 1, na
776 force_b = force_b + deriv_ab(i, j, :)*pab(i, j)
777 END DO
778 END DO
779 END SUBROUTINE contract_basis_center_forces
780
781END MODULE libgrpp_integrals
libgrpp::libgrpp_init
Definition libgrpp.F:41
libgrpp::libgrpp_type1_integrals_gradient
Definition libgrpp.F:147
libgrpp::libgrpp_type1_integrals
Definition libgrpp.F:73
libgrpp::libgrpp_type2_integrals_gradient
Definition libgrpp.F:172
libgrpp::libgrpp_type2_integrals
Definition libgrpp.F:97
ai_cartesian_operators
Coordinate and basis-center derivative ladders for Cartesian Gaussian integrals. Derivatives act on a...
Definition ai_cartesian_operators.F:13
ai_cartesian_operators::apply_cartesian_operator
subroutine, public apply_cartesian_operator(la_max, la_min, npgfa, lb_max, lb_min, npgfb, ab, values, right, derivative, exponents, origin)
Apply a Cartesian operator on either side of a primitive integral block. The input includes one extra...
Definition ai_cartesian_operators.F:192
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
libgrpp_integrals
Local and semi-local ECP integrals using the libgrpp library.
Definition libgrpp_integrals.F:12
libgrpp_integrals::libgrpp_semilocal_integrals
subroutine, public libgrpp_semilocal_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
Semi-local ECP integrals using libgrpp.
Definition libgrpp_integrals.F:274
libgrpp_integrals::libgrpp_local_integrals
subroutine, public libgrpp_local_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
Local ECP integrals using libgrpp.
Definition libgrpp_integrals.F:66
libgrpp_integrals::libgrpp_local_forces_ref
subroutine, public libgrpp_local_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
Reference local ECP force routine using l+-1 integrals. No call is made to the numerically unstable g...
Definition libgrpp_integrals.F:490
libgrpp_integrals::libgrpp_semilocal_forces_ref
subroutine, public libgrpp_semilocal_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
Reference semi-local ECP forces using l+-1 integrals. No call is made to the numerically unstable gra...
Definition libgrpp_integrals.F:629
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42