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