(git:ed6f26b)
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 kinds, ONLY: dp
14 USE mathconstants, ONLY: pi
15 USE ai_derivatives, ONLY: dabdr_noscreen, adbdr, dabdr
16 USE orbital_pointers, ONLY: nco, &
17 ncoset
18#if defined(__LIBGRPP)
19 USE libgrpp, ONLY: libgrpp_init, libgrpp_type1_integrals, libgrpp_type2_integrals, &
20 libgrpp_type1_integrals_gradient, libgrpp_type2_integrals_gradient
21#endif
22#include "./base/base_uses.f90"
23
24 IMPLICIT NONE
25 PRIVATE
26
27 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libgrpp_integrals'
28
29 PUBLIC :: libgrpp_semilocal_integrals, libgrpp_local_integrals, &
30 libgrpp_local_forces_ref, libgrpp_semilocal_forces_ref
31
32CONTAINS
33
34! **************************************************************************************************
35!> \brief Local ECP integrals using libgrpp
36!> \param la_max_set ...
37!> \param la_min_set ...
38!> \param npgfa ...
39!> \param rpgfa ...
40!> \param zeta ...
41!> \param lb_max_set ...
42!> \param lb_min_set ...
43!> \param npgfb ...
44!> \param rpgfb ...
45!> \param zetb ...
46!> \param npot_ecp ...
47!> \param alpha_ecp ...
48!> \param coeffs_ecp ...
49!> \param nrpot_ecp ...
50!> \param rpgfc ...
51!> \param rab ...
52!> \param dab ...
53!> \param rac ...
54!> \param dac ...
55!> \param dbc ...
56!> \param vab ...
57!> \param pab ...
58!> \param force_a ...
59!> \param force_b ...
60! **************************************************************************************************
61 SUBROUTINE libgrpp_local_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
62 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
63 npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
64 rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
65
66 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
67 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
68 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
69 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
70 INTEGER, INTENT(IN) :: npot_ecp
71 REAL(kind=dp), DIMENSION(1:npot_ecp), INTENT(IN) :: alpha_ecp, coeffs_ecp
72 INTEGER, DIMENSION(1:npot_ecp), INTENT(IN) :: nrpot_ecp
73 REAL(kind=dp), INTENT(IN) :: rpgfc
74 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
75 REAL(kind=dp), INTENT(IN) :: dab
76 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
77 REAL(kind=dp), INTENT(IN) :: dac
78 REAL(kind=dp), INTENT(IN) :: 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#if defined(__LIBGRPP)
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, normi, normj, prefi, prefj, &
90 zeti, zetj, mindist, fac_a, fac_b
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#else
242
243 mark_used(la_max_set)
244 mark_used(la_min_set)
245 mark_used(npgfa)
246 mark_used(rpgfa)
247 mark_used(zeta)
248 mark_used(lb_max_set)
249 mark_used(lb_min_set)
250 mark_used(npgfb)
251 mark_used(rpgfb)
252 mark_used(zetb)
253 mark_used(npot_ecp)
254 mark_used(alpha_ecp)
255 mark_used(coeffs_ecp)
256 mark_used(nrpot_ecp)
257 mark_used(rpgfc)
258 mark_used(rab)
259 mark_used(dab)
260 mark_used(rac)
261 mark_used(dac)
262 mark_used(dbc)
263 mark_used(vab)
264 mark_used(pab)
265 mark_used(force_a)
266 mark_used(force_b)
267
268 cpabort("Please compile CP2K with libgrpp support for calculations with ECPs")
269#endif
270
271 END SUBROUTINE libgrpp_local_integrals
272
273! **************************************************************************************************
274!> \brief Semi-local ECP integrals using libgrpp.
275!> \param la_max_set ...
276!> \param la_min_set ...
277!> \param npgfa ...
278!> \param rpgfa ...
279!> \param zeta ...
280!> \param lb_max_set ...
281!> \param lb_min_set ...
282!> \param npgfb ...
283!> \param rpgfb ...
284!> \param zetb ...
285!> \param lmax_ecp ...
286!> \param npot_ecp ...
287!> \param alpha_ecp ...
288!> \param coeffs_ecp ...
289!> \param nrpot_ecp ...
290!> \param rpgfc ...
291!> \param rab ...
292!> \param dab ...
293!> \param rac ...
294!> \param dac ...
295!> \param dbc ...
296!> \param vab ...
297!> \param pab ...
298!> \param force_a ...
299!> \param force_b ...
300! **************************************************************************************************
301 SUBROUTINE libgrpp_semilocal_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
302 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
303 lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
304 rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
305
306 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
307 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
308 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
309 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
310 INTEGER, INTENT(IN) :: lmax_ecp
311 INTEGER, DIMENSION(0:10), INTENT(IN) :: npot_ecp
312 REAL(kind=dp), DIMENSION(1:15, 0:10), INTENT(IN) :: alpha_ecp, coeffs_ecp
313 INTEGER, DIMENSION(1:15, 0:10), INTENT(IN) :: nrpot_ecp
314 REAL(kind=dp), INTENT(IN) :: rpgfc
315 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
316 REAL(kind=dp), INTENT(IN) :: dab
317 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
318 REAL(kind=dp), INTENT(IN) :: dac
319 REAL(kind=dp), INTENT(IN) :: dbc
320 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
321 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
322 OPTIONAL :: pab
323 REAL(kind=dp), DIMENSION(3), INTENT(INOUT), &
324 OPTIONAL :: force_a, force_b
325
326#if defined(__LIBGRPP)
327 INTEGER :: a_offset, a_start, b_offset, b_start, i, &
328 ipgf, j, jpgf, li, lj, lk, ncoa, ncob
329 LOGICAL :: calc_forces
330 REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
331 zeti, zetj, mindist, fac_a, fac_b
332 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmpx, tmpz, tmpy
333 REAL(dp), DIMENSION(3) :: ra, rb, rc
334
335 CALL libgrpp_init()
336
337 calc_forces = .false.
338 IF (PRESENT(pab) .AND. PRESENT(force_a) .AND. PRESENT(force_b)) calc_forces = .true.
339
340 IF (calc_forces) THEN
341
342 !Note: warning against numerical stability of libgrpp gradients. The day the library becomes
343 ! stable, this routine can be used immediatly as is, and the warning removed.
344 CALL cp_warn(__location__, &
345 "ECP gradients calculated with the libgrpp library are, to this date, not numerically stable. "// &
346 "Please use the reference routine 'libgrpp_semilocal_forces_ref' instead.")
347
348 !there is a weird feature of libgrpp gradients, which is such that the gradient is calculated
349 !for a point in space, and not with respect to an atomic center. For example, if atoms A and
350 !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
351 !Because we want the forces on centers A and B seprately, we need a case study on atomic positions
352 !We always calculate the gradient wrt to atomic position of A and B, and we scale accordingly
353
354 mindist = 1.0e-6_dp
355 !If ra != rb != rc
356 IF (dab >= mindist .AND. dbc >= mindist .AND. dac >= mindist) THEN
357 fac_a = 1.0_dp
358 fac_b = 1.0_dp
359
360 !If ra = rb, but ra != rc
361 ELSE IF (dab < mindist .AND. dac >= mindist) THEN
362 fac_a = 0.5_dp
363 fac_b = 0.5_dp
364
365 !IF ra != rb but ra = rc
366 ELSE IF (dab >= mindist .AND. dac < mindist) THEN
367 fac_a = 0.5_dp
368 fac_b = 1.0_dp
369
370 !IF ra != rb but rb = rc
371 ELSE IF (dab >= mindist .AND. dbc < mindist) THEN
372 fac_a = 1.0_dp
373 fac_b = 0.5_dp
374
375 !If all atoms the same --> no force
376 ELSE
377 calc_forces = .false.
378 END IF
379 END IF
380
381 !libgrpp requires absolute positions, not relative ones
382 ra(:) = 0.0_dp
383 rb(:) = rab(:)
384 rc(:) = rac(:)
385
386 ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
387 IF (calc_forces) THEN
388 ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
389 ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
390 ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
391 END IF
392
393 DO ipgf = 1, npgfa
394 IF (rpgfa(ipgf) + rpgfc < dac) cycle
395 zeti = zeta(ipgf)
396 a_start = (ipgf - 1)*ncoset(la_max_set)
397
398 DO jpgf = 1, npgfb
399 IF (rpgfb(jpgf) + rpgfc < dbc) cycle
400 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
401 zetj = zetb(jpgf)
402 b_start = (jpgf - 1)*ncoset(lb_max_set)
403
404 DO li = la_min_set, la_max_set
405 a_offset = a_start + ncoset(li - 1)
406 ncoa = nco(li)
407 prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
408 expi = 0.25_dp*real(2*li + 3, dp)
409 normi = 1.0_dp/(prefi*zeti**expi)
410
411 DO lj = lb_min_set, lb_max_set
412 b_offset = b_start + ncoset(lj - 1)
413 ncob = nco(lj)
414 prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
415 expj = 0.25_dp*real(2*lj + 3, dp)
416 normj = 1.0_dp/(prefj*zetj**expj)
417
418 !Loop over ECP angular momentum
419 DO lk = 0, lmax_ecp
420 tmp(1:ncoa*ncob) = 0.0_dp
421 !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
422 !the 1/norm coefficients for PGFi and PGFj
423 CALL libgrpp_type2_integrals(ra, li, 1, [normi], [zeti], &
424 rb, lj, 1, [normj], [zetj], &
425 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
426 coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
427
428 !note: tmp array is in C row-major ordering
429 DO j = 1, ncob
430 DO i = 1, ncoa
431 vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
432 END DO
433 END DO
434
435 IF (calc_forces) THEN
436
437 tmpx(1:ncoa*ncob) = 0.0_dp
438 tmpy(1:ncoa*ncob) = 0.0_dp
439 tmpz(1:ncoa*ncob) = 0.0_dp
440
441 !force wrt to atomic position A
442 CALL libgrpp_type2_integrals_gradient(ra, li, 1, [normi], [zeti], &
443 rb, lj, 1, [normj], [zetj], &
444 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
445 coeffs_ecp(:, lk), alpha_ecp(:, lk), ra, &
446 tmpx, tmpy, tmpz)
447
448 !note: tmp array is in C row-major ordering
449 !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
450 DO j = 1, ncob
451 DO i = 1, ncoa
452 force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
453 force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
454 force_a(3) = force_a(3) + fac_a*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
455 END DO
456 END DO
457
458 tmpx(1:ncoa*ncob) = 0.0_dp
459 tmpy(1:ncoa*ncob) = 0.0_dp
460 tmpz(1:ncoa*ncob) = 0.0_dp
461
462 !force wrt to atomic position B
463 CALL libgrpp_type2_integrals_gradient(ra, li, 1, [normi], [zeti], &
464 rb, lj, 1, [normj], [zetj], &
465 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
466 coeffs_ecp(:, lk), alpha_ecp(:, lk), rb, &
467 tmpx, tmpy, tmpz)
468 !note: tmp array is in C row-major ordering
469 !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
470 DO j = 1, ncob
471 DO i = 1, ncoa
472 force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
473 force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
474 force_b(3) = force_b(3) + fac_b*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
475 END DO
476 END DO
477
478 END IF !calc_forces
479
480 END DO !lk
481
482 END DO !lj
483 END DO !li
484
485 END DO !jpgf
486 END DO !ipgf
487
488#else
489
490 mark_used(la_max_set)
491 mark_used(la_min_set)
492 mark_used(npgfa)
493 mark_used(rpgfa)
494 mark_used(zeta)
495 mark_used(lb_max_set)
496 mark_used(lb_min_set)
497 mark_used(npgfb)
498 mark_used(rpgfb)
499 mark_used(zetb)
500 mark_used(lmax_ecp)
501 mark_used(npot_ecp)
502 mark_used(alpha_ecp)
503 mark_used(coeffs_ecp)
504 mark_used(nrpot_ecp)
505 mark_used(rpgfc)
506 mark_used(rab)
507 mark_used(dab)
508 mark_used(rac)
509 mark_used(dac)
510 mark_used(dbc)
511 mark_used(vab)
512 mark_used(pab)
513 mark_used(force_a)
514 mark_used(force_b)
515
516 cpabort("Please compile CP2K with libgrpp support for calculations with ECPs")
517#endif
518
519 END SUBROUTINE libgrpp_semilocal_integrals
520
521! **************************************************************************************************
522!> \brief Reference local ECP force routine using l+-1 integrals. No call is made to the numerically
523!> unstable gradient routine of libgrpp. Calculates both the integrals and the forces.
524!> \param la_max_set ...
525!> \param la_min_set ...
526!> \param npgfa ...
527!> \param rpgfa ...
528!> \param zeta ...
529!> \param lb_max_set ...
530!> \param lb_min_set ...
531!> \param npgfb ...
532!> \param rpgfb ...
533!> \param zetb ...
534!> \param npot_ecp ...
535!> \param alpha_ecp ...
536!> \param coeffs_ecp ...
537!> \param nrpot_ecp ...
538!> \param rpgfc ...
539!> \param rab ...
540!> \param dab ...
541!> \param rac ...
542!> \param dac ...
543!> \param dbc ...
544!> \param vab ...
545!> \param pab ...
546!> \param force_a ...
547!> \param force_b ...
548!> \note: this is a reference routine, which has no reason to be used once the libgrpp gradients
549!> become numerically stable
550! **************************************************************************************************
551 SUBROUTINE libgrpp_local_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
552 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
553 npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
554 rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
555
556 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
557 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
558 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
559 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
560 INTEGER, INTENT(IN) :: npot_ecp
561 REAL(kind=dp), DIMENSION(1:npot_ecp), INTENT(IN) :: alpha_ecp, coeffs_ecp
562 INTEGER, DIMENSION(1:npot_ecp), INTENT(IN) :: nrpot_ecp
563 REAL(kind=dp), INTENT(IN) :: rpgfc
564 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
565 REAL(kind=dp), INTENT(IN) :: dab
566 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
567 REAL(kind=dp), INTENT(IN) :: dac
568 REAL(kind=dp), INTENT(IN) :: dbc
569 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
570 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pab
571 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: force_a, force_b
572
573#if defined(__LIBGRPP)
574 INTEGER :: a_offset, a_start, b_offset, b_start, i, &
575 ipgf, j, jpgf, li, lj, ncoa, ncob, a_offset_f, &
576 b_offset_f, a_start_f, b_start_f
577 REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
578 zeti, zetj
579 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
580 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vab_f, tmpx, tmpy, tmpz
581 REAL(dp), DIMENSION(3) :: ra, rb, rc
582
583 CALL libgrpp_init()
584
585 !Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
586 ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
587 vab_f(:, :) = 0.0_dp
588
589 !libgrpp requires absolute positions, not relative ones
590 ra(:) = 0.0_dp
591 rb(:) = rab(:)
592 rc(:) = rac(:)
593
594 ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
595
596 DO ipgf = 1, npgfa
597 IF (rpgfa(ipgf) + rpgfc < dac) cycle
598 zeti = zeta(ipgf)
599 a_start = (ipgf - 1)*ncoset(la_max_set)
600 a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
601
602 DO jpgf = 1, npgfb
603 IF (rpgfb(jpgf) + rpgfc < dbc) cycle
604 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
605 zetj = zetb(jpgf)
606 b_start = (jpgf - 1)*ncoset(lb_max_set)
607 b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
608
609 DO li = max(0, la_min_set - 1), la_max_set + 1
610 a_offset = a_start + ncoset(li - 1)
611 a_offset_f = a_start_f + ncoset(li - 1)
612 ncoa = nco(li)
613 prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
614 expi = 0.25_dp*real(2*li + 3, dp)
615 normi = 1.0_dp/(prefi*zeti**expi)
616
617 DO lj = max(0, lb_min_set - 1), lb_max_set + 1
618 b_offset = b_start + ncoset(lj - 1)
619 b_offset_f = b_start_f + ncoset(lj - 1)
620 ncob = nco(lj)
621 prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
622 expj = 0.25_dp*real(2*lj + 3, dp)
623 normj = 1.0_dp/(prefj*zetj**expj)
624
625 tmp(1:ncoa*ncob) = 0.0_dp
626 !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
627 !the 1/norm coefficients for PGFi and PGFj
628 CALL libgrpp_type1_integrals(ra, li, 1, [normi], [zeti], &
629 rb, lj, 1, [normj], [zetj], &
630 rc, [npot_ecp], nrpot_ecp, &
631 coeffs_ecp, alpha_ecp, tmp)
632
633 !the l+-1 integrals for gradient calculation
634 DO j = 1, ncob
635 DO i = 1, ncoa
636 vab_f(a_offset_f + i, b_offset_f + j) = &
637 vab_f(a_offset_f + i, b_offset_f + j) + tmp((i - 1)*ncob + j)
638 END DO
639 END DO
640
641 !the actual integrals
642 IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
643 DO j = 1, ncob
644 DO i = 1, ncoa
645 vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
646 END DO
647 END DO
648 END IF
649
650 END DO !lj
651 END DO !li
652
653 END DO !jpgf
654 END DO !ipgf
655
656 ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
657 ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
658 ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
659
660 !Derivative wrt to center A
661 tmpx(:, :) = 0.0_dp
662 tmpy(:, :) = 0.0_dp
663 tmpz(:, :) = 0.0_dp
664 CALL dabdr(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, rpgfb, lb_min_set, &
665 dab, vab_f, tmpx, tmpy, tmpz)
666 DO j = 1, npgfb*ncoset(lb_max_set)
667 DO i = 1, npgfa*ncoset(la_max_set)
668 force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
669 force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
670 force_a(3) = force_a(3) + tmpz(i, j)*pab(i, j)
671 END DO
672 END DO
673
674 !Derivative wrt to center B
675 tmpx(:, :) = 0.0_dp
676 tmpy(:, :) = 0.0_dp
677 tmpz(:, :) = 0.0_dp
678 CALL adbdr(la_max_set, npgfa, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
679 dab, vab_f, tmpx, tmpy, tmpz)
680 DO j = 1, npgfb*ncoset(lb_max_set)
681 DO i = 1, npgfa*ncoset(la_max_set)
682 force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
683 force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
684 force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
685 END DO
686 END DO
687 DEALLOCATE (tmpx, tmpy, tmpz)
688
689#else
690
691 mark_used(la_max_set)
692 mark_used(la_min_set)
693 mark_used(npgfa)
694 mark_used(rpgfa)
695 mark_used(zeta)
696 mark_used(lb_max_set)
697 mark_used(lb_min_set)
698 mark_used(npgfb)
699 mark_used(rpgfb)
700 mark_used(zetb)
701 mark_used(npot_ecp)
702 mark_used(alpha_ecp)
703 mark_used(coeffs_ecp)
704 mark_used(nrpot_ecp)
705 mark_used(rpgfc)
706 mark_used(rab)
707 mark_used(dab)
708 mark_used(rac)
709 mark_used(dac)
710 mark_used(dbc)
711 mark_used(pab)
712 mark_used(vab)
713 mark_used(force_a)
714 mark_used(force_b)
715
716 cpabort("Please compile CP2K with libgrpp support for calculations with ECPs")
717#endif
718
719 END SUBROUTINE libgrpp_local_forces_ref
720
721! **************************************************************************************************
722!> \brief Reference semi-local ECP forces using l+-1 integrals. No call is made to the numerically
723!> unstable gradient routine of libgrpp. Calculates both the integrals and the forces.
724!> \param la_max_set ...
725!> \param la_min_set ...
726!> \param npgfa ...
727!> \param rpgfa ...
728!> \param zeta ...
729!> \param lb_max_set ...
730!> \param lb_min_set ...
731!> \param npgfb ...
732!> \param rpgfb ...
733!> \param zetb ...
734!> \param lmax_ecp ...
735!> \param npot_ecp ...
736!> \param alpha_ecp ...
737!> \param coeffs_ecp ...
738!> \param nrpot_ecp ...
739!> \param rpgfc ...
740!> \param rab ...
741!> \param dab ...
742!> \param rac ...
743!> \param dac ...
744!> \param dbc ...
745!> \param vab ...
746!> \param pab ...
747!> \param force_a ...
748!> \param force_b ...
749!> \note: this is a reference routine, which has no reason to be used once the libgrpp gradients
750!> become numerically stable
751! **************************************************************************************************
752 SUBROUTINE libgrpp_semilocal_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
753 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
754 lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
755 rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
756
757 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
758 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
759 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
760 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
761 INTEGER, INTENT(IN) :: lmax_ecp
762 INTEGER, DIMENSION(0:10), INTENT(IN) :: npot_ecp
763 REAL(kind=dp), DIMENSION(1:15, 0:10), INTENT(IN) :: alpha_ecp, coeffs_ecp
764 INTEGER, DIMENSION(1:15, 0:10), INTENT(IN) :: nrpot_ecp
765 REAL(kind=dp), INTENT(IN) :: rpgfc
766 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
767 REAL(kind=dp), INTENT(IN) :: dab
768 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
769 REAL(kind=dp), INTENT(IN) :: dac
770 REAL(kind=dp), INTENT(IN) :: dbc
771 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
772 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pab
773 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: force_a, force_b
774
775#if defined(__LIBGRPP)
776 INTEGER :: a_offset, a_start, b_offset, b_start, i, &
777 ipgf, j, jpgf, li, lj, lk, ncoa, ncob, &
778 a_start_f, b_start_f, a_offset_f, b_offset_f
779 REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
780 zeti, zetj
781 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
782 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vab_f, tmpx, tmpy, tmpz
783 REAL(dp), DIMENSION(3) :: ra, rb, rc
784
785 CALL libgrpp_init()
786
787 !Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
788 ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
789 vab_f(:, :) = 0.0_dp
790
791 !libgrpp requires absolute positions, not relative ones
792 ra(:) = 0.0_dp
793 rb(:) = rab(:)
794 rc(:) = rac(:)
795
796 ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
797
798 DO ipgf = 1, npgfa
799 IF (rpgfa(ipgf) + rpgfc < dac) cycle
800 zeti = zeta(ipgf)
801 a_start = (ipgf - 1)*ncoset(la_max_set)
802 a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
803
804 DO jpgf = 1, npgfb
805 IF (rpgfb(jpgf) + rpgfc < dbc) cycle
806 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
807 zetj = zetb(jpgf)
808 b_start = (jpgf - 1)*ncoset(lb_max_set)
809 b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
810
811 DO li = max(0, la_min_set - 1), la_max_set + 1
812 a_offset = a_start + ncoset(li - 1)
813 a_offset_f = a_start_f + ncoset(li - 1)
814 ncoa = nco(li)
815 prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
816 expi = 0.25_dp*real(2*li + 3, dp)
817 normi = 1.0_dp/(prefi*zeti**expi)
818
819 DO lj = max(0, lb_min_set - 1), lb_max_set + 1
820 b_offset = b_start + ncoset(lj - 1)
821 b_offset_f = b_start_f + ncoset(lj - 1)
822 ncob = nco(lj)
823 prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
824 expj = 0.25_dp*real(2*lj + 3, dp)
825 normj = 1.0_dp/(prefj*zetj**expj)
826
827 !Loop over ECP angular momentum
828 DO lk = 0, lmax_ecp
829 tmp(1:ncoa*ncob) = 0.0_dp
830 !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
831 !the 1/norm coefficients for PGFi and PGFj
832 CALL libgrpp_type2_integrals(ra, li, 1, [normi], [zeti], &
833 rb, lj, 1, [normj], [zetj], &
834 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
835 coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
836
837 !the l+-1 integrals for gradient calculation
838 DO j = 1, ncob
839 DO i = 1, ncoa
840 vab_f(a_offset_f + i, b_offset_f + j) = &
841 vab_f(a_offset_f + i, b_offset_f + j) + tmp((i - 1)*ncob + j)
842 END DO
843 END DO
844
845 !the actual integrals
846 IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
847 DO j = 1, ncob
848 DO i = 1, ncoa
849 vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
850 END DO
851 END DO
852 END IF
853
854 END DO !lk
855
856 END DO !lj
857 END DO !li
858
859 END DO !jpgf
860 END DO !ipgf
861
862 ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
863 ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
864 ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
865
866 !Derivative wrt to center A
867 tmpx(:, :) = 0.0_dp
868 tmpy(:, :) = 0.0_dp
869 tmpz(:, :) = 0.0_dp
870 CALL dabdr(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, rpgfb, lb_min_set, &
871 0.0_dp, vab_f, tmpx, tmpy, tmpz)
872 DO j = 1, npgfb*ncoset(lb_max_set)
873 DO i = 1, npgfa*ncoset(la_max_set)
874 force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
875 force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
876 force_a(3) = force_a(3) + tmpz(i, j)*pab(i, j)
877 END DO
878 END DO
879
880 !Derivative wrt to center B
881 tmpx(:, :) = 0.0_dp
882 tmpy(:, :) = 0.0_dp
883 tmpz(:, :) = 0.0_dp
884 CALL adbdr(la_max_set, npgfa, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
885 0.0_dp, vab_f, tmpx, tmpy, tmpz)
886 DO j = 1, npgfb*ncoset(lb_max_set)
887 DO i = 1, npgfa*ncoset(la_max_set)
888 force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
889 force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
890 force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
891 END DO
892 END DO
893 DEALLOCATE (tmpx, tmpy, tmpz)
894
895#else
896
897 mark_used(la_max_set)
898 mark_used(la_min_set)
899 mark_used(npgfa)
900 mark_used(rpgfa)
901 mark_used(zeta)
902 mark_used(lb_max_set)
903 mark_used(lb_min_set)
904 mark_used(npgfb)
905 mark_used(rpgfb)
906 mark_used(zetb)
907 mark_used(lmax_ecp)
908 mark_used(npot_ecp)
909 mark_used(alpha_ecp)
910 mark_used(coeffs_ecp)
911 mark_used(nrpot_ecp)
912 mark_used(rpgfc)
913 mark_used(rab)
914 mark_used(dab)
915 mark_used(rac)
916 mark_used(dac)
917 mark_used(dbc)
918 mark_used(pab)
919 mark_used(vab)
920 mark_used(force_a)
921 mark_used(force_b)
922
923 cpabort("Please compile CP2K with libgrpp support for calculations with ECPs")
924#endif
925
926 END SUBROUTINE libgrpp_semilocal_forces_ref
927
928END 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::dabdr_noscreen
subroutine, public dabdr_noscreen(la_max, npgfa, zeta, lb_max, npgfb, 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:366
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:305
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:65
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:555
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:756
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