(git:aeba166)
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-2024 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_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 calc_forces = .false.
95 IF (PRESENT(pab) .AND. PRESENT(force_a) .AND. PRESENT(force_b)) calc_forces = .true.
96
97 IF (calc_forces) THEN
98
99 !Note: warning against numerical stability of libgrpp gradients. The day the library becomes
100 ! stable, this routine can be used immediatly as is, and the warning removed.
101 CALL cp_warn(__location__, &
102 "ECP gradients calculated with the libgrpp library are, to this date, not numerically stable. "// &
103 "Please use the reference routine 'libgrpp_local_forces_ref' instead.")
104
105 !there is a weird feature of libgrpp gradients, which is such that the gradient is calculated
106 !for a point in space, and not with respect to an atomic center. For example, if atoms A and
107 !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
108 !Because we want the forces on centers A and B seprately, we need a case study on atomic positions
109 !We always calculate the gradient wrt to atomic position of A and B, and we scale accordingly
110
111 mindist = 1.0e-6_dp
112 !If ra != rb != rc
113 IF (dab >= mindist .AND. dbc >= mindist .AND. dac >= mindist) THEN
114 fac_a = 1.0_dp
115 fac_b = 1.0_dp
116
117 !If ra = rb, but ra != rc
118 ELSE IF (dab < mindist .AND. dac >= mindist) THEN
119 fac_a = 0.5_dp
120 fac_b = 0.5_dp
121
122 !IF ra != rb but ra = rc
123 ELSE IF (dab >= mindist .AND. dac < mindist) THEN
124 fac_a = 0.5_dp
125 fac_b = 1.0_dp
126
127 !IF ra != rb but rb = rc
128 ELSE IF (dab >= mindist .AND. dbc < mindist) THEN
129 fac_a = 1.0_dp
130 fac_b = 0.5_dp
131
132 !If all atoms the same --> no force
133 ELSE
134 calc_forces = .false.
135 END IF
136 END IF
137
138 !libgrpp requires absolute positions, not relative ones
139 ra(:) = 0.0_dp
140 rb(:) = rab(:)
141 rc(:) = rac(:)
142
143 ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
144 IF (calc_forces) THEN
145 ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
146 ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
147 ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
148 END IF
149
150 DO ipgf = 1, npgfa
151 IF (rpgfa(ipgf) + rpgfc < dac) cycle
152 zeti = zeta(ipgf)
153 a_start = (ipgf - 1)*ncoset(la_max_set)
154
155 DO jpgf = 1, npgfb
156 IF (rpgfb(jpgf) + rpgfc < dbc) cycle
157 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
158 zetj = zetb(jpgf)
159 b_start = (jpgf - 1)*ncoset(lb_max_set)
160
161 DO li = la_min_set, la_max_set
162 a_offset = a_start + ncoset(li - 1)
163 ncoa = nco(li)
164 prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
165 expi = 0.25_dp*real(2*li + 3, dp)
166 normi = 1.0_dp/(prefi*zeti**expi)
167
168 DO lj = lb_min_set, lb_max_set
169 b_offset = b_start + ncoset(lj - 1)
170 ncob = nco(lj)
171 prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
172 expj = 0.25_dp*real(2*lj + 3, dp)
173 normj = 1.0_dp/(prefj*zetj**expj)
174
175 tmp(1:ncoa*ncob) = 0.0_dp
176 !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
177 !the 1/norm coefficients for PGFi and PGFj
178 CALL libgrpp_type1_integrals(ra, li, 1, [normi], [zeti], &
179 rb, lj, 1, [normj], [zetj], &
180 rc, [npot_ecp], nrpot_ecp, &
181 coeffs_ecp, alpha_ecp, tmp)
182
183 !note: tmp array is in C row-major ordering
184 DO j = 1, ncob
185 DO i = 1, ncoa
186 vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
187 END DO
188 END DO
189
190 IF (calc_forces) THEN
191 tmpx(1:ncoa*ncob) = 0.0_dp
192 tmpy(1:ncoa*ncob) = 0.0_dp
193 tmpz(1:ncoa*ncob) = 0.0_dp
194
195 !force wrt to atomic position A
196 CALL libgrpp_type1_integrals_gradient(ra, li, 1, [normi], [zeti], &
197 rb, lj, 1, [normj], [zetj], &
198 rc, [npot_ecp], nrpot_ecp, &
199 coeffs_ecp, alpha_ecp, ra, &
200 tmpx, tmpy, tmpz)
201
202 !note: tmp array is in C row-major ordering
203 !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
204 DO j = 1, ncob
205 DO i = 1, ncoa
206 force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
207 force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
208 force_a(3) = force_a(3) + fac_a*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
209 END DO
210 END DO
211
212 tmpx(1:ncoa*ncob) = 0.0_dp
213 tmpy(1:ncoa*ncob) = 0.0_dp
214 tmpz(1:ncoa*ncob) = 0.0_dp
215
216 !force wrt to atomic position B
217 CALL libgrpp_type1_integrals_gradient(ra, li, 1, [normi], [zeti], &
218 rb, lj, 1, [normj], [zetj], &
219 rc, [npot_ecp], nrpot_ecp, &
220 coeffs_ecp, alpha_ecp, rb, &
221 tmpx, tmpy, tmpz)
222
223 !note: tmp array is in C row-major ordering
224 !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
225 DO j = 1, ncob
226 DO i = 1, ncoa
227 force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
228 force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
229 force_b(3) = force_b(3) + fac_b*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
230 END DO
231 END DO
232 END IF
233
234 END DO !lj
235 END DO !li
236
237 END DO !jpgf
238 END DO !ipgf
239#else
240
241 mark_used(la_max_set)
242 mark_used(la_min_set)
243 mark_used(npgfa)
244 mark_used(rpgfa)
245 mark_used(zeta)
246 mark_used(lb_max_set)
247 mark_used(lb_min_set)
248 mark_used(npgfb)
249 mark_used(rpgfb)
250 mark_used(zetb)
251 mark_used(npot_ecp)
252 mark_used(alpha_ecp)
253 mark_used(coeffs_ecp)
254 mark_used(nrpot_ecp)
255 mark_used(rpgfc)
256 mark_used(rab)
257 mark_used(dab)
258 mark_used(rac)
259 mark_used(dac)
260 mark_used(dbc)
261 mark_used(vab)
262 mark_used(pab)
263 mark_used(force_a)
264 mark_used(force_b)
265
266 cpabort("Please compile CP2K with libgrpp support for calculations with ECPs")
267#endif
268
269 END SUBROUTINE libgrpp_local_integrals
270
271! **************************************************************************************************
272!> \brief Semi-local ECP integrals using libgrpp.
273!> \param la_max_set ...
274!> \param la_min_set ...
275!> \param npgfa ...
276!> \param rpgfa ...
277!> \param zeta ...
278!> \param lb_max_set ...
279!> \param lb_min_set ...
280!> \param npgfb ...
281!> \param rpgfb ...
282!> \param zetb ...
283!> \param lmax_ecp ...
284!> \param npot_ecp ...
285!> \param alpha_ecp ...
286!> \param coeffs_ecp ...
287!> \param nrpot_ecp ...
288!> \param rpgfc ...
289!> \param rab ...
290!> \param dab ...
291!> \param rac ...
292!> \param dac ...
293!> \param dbc ...
294!> \param vab ...
295!> \param pab ...
296!> \param force_a ...
297!> \param force_b ...
298! **************************************************************************************************
299 SUBROUTINE libgrpp_semilocal_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
300 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
301 lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
302 rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
303
304 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
305 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
306 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
307 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
308 INTEGER, INTENT(IN) :: lmax_ecp
309 INTEGER, DIMENSION(0:10), INTENT(IN) :: npot_ecp
310 REAL(kind=dp), DIMENSION(1:15, 0:10), INTENT(IN) :: alpha_ecp, coeffs_ecp
311 INTEGER, DIMENSION(1:15, 0:10), INTENT(IN) :: nrpot_ecp
312 REAL(kind=dp), INTENT(IN) :: rpgfc
313 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
314 REAL(kind=dp), INTENT(IN) :: dab
315 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
316 REAL(kind=dp), INTENT(IN) :: dac
317 REAL(kind=dp), INTENT(IN) :: dbc
318 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
319 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
320 OPTIONAL :: pab
321 REAL(kind=dp), DIMENSION(3), INTENT(INOUT), &
322 OPTIONAL :: force_a, force_b
323
324#if defined(__LIBGRPP)
325 INTEGER :: a_offset, a_start, b_offset, b_start, i, &
326 ipgf, j, jpgf, li, lj, lk, ncoa, ncob
327 LOGICAL :: calc_forces
328 REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
329 zeti, zetj, mindist, fac_a, fac_b
330 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmpx, tmpz, tmpy
331 REAL(dp), DIMENSION(3) :: ra, rb, rc
332
333 calc_forces = .false.
334 IF (PRESENT(pab) .AND. PRESENT(force_a) .AND. PRESENT(force_b)) calc_forces = .true.
335
336 IF (calc_forces) THEN
337
338 !Note: warning against numerical stability of libgrpp gradients. The day the library becomes
339 ! stable, this routine can be used immediatly as is, and the warning removed.
340 CALL cp_warn(__location__, &
341 "ECP gradients calculated with the libgrpp library are, to this date, not numerically stable. "// &
342 "Please use the reference routine 'libgrpp_semilocal_forces_ref' instead.")
343
344 !there is a weird feature of libgrpp gradients, which is such that the gradient is calculated
345 !for a point in space, and not with respect to an atomic center. For example, if atoms A and
346 !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
347 !Because we want the forces on centers A and B seprately, we need a case study on atomic positions
348 !We always calculate the gradient wrt to atomic position of A and B, and we scale accordingly
349
350 mindist = 1.0e-6_dp
351 !If ra != rb != rc
352 IF (dab >= mindist .AND. dbc >= mindist .AND. dac >= mindist) THEN
353 fac_a = 1.0_dp
354 fac_b = 1.0_dp
355
356 !If ra = rb, but ra != rc
357 ELSE IF (dab < mindist .AND. dac >= mindist) THEN
358 fac_a = 0.5_dp
359 fac_b = 0.5_dp
360
361 !IF ra != rb but ra = rc
362 ELSE IF (dab >= mindist .AND. dac < mindist) THEN
363 fac_a = 0.5_dp
364 fac_b = 1.0_dp
365
366 !IF ra != rb but rb = rc
367 ELSE IF (dab >= mindist .AND. dbc < mindist) THEN
368 fac_a = 1.0_dp
369 fac_b = 0.5_dp
370
371 !If all atoms the same --> no force
372 ELSE
373 calc_forces = .false.
374 END IF
375 END IF
376
377 !libgrpp requires absolute positions, not relative ones
378 ra(:) = 0.0_dp
379 rb(:) = rab(:)
380 rc(:) = rac(:)
381
382 ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
383 IF (calc_forces) THEN
384 ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
385 ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
386 ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
387 END IF
388
389 DO ipgf = 1, npgfa
390 IF (rpgfa(ipgf) + rpgfc < dac) cycle
391 zeti = zeta(ipgf)
392 a_start = (ipgf - 1)*ncoset(la_max_set)
393
394 DO jpgf = 1, npgfb
395 IF (rpgfb(jpgf) + rpgfc < dbc) cycle
396 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
397 zetj = zetb(jpgf)
398 b_start = (jpgf - 1)*ncoset(lb_max_set)
399
400 DO li = la_min_set, la_max_set
401 a_offset = a_start + ncoset(li - 1)
402 ncoa = nco(li)
403 prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
404 expi = 0.25_dp*real(2*li + 3, dp)
405 normi = 1.0_dp/(prefi*zeti**expi)
406
407 DO lj = lb_min_set, lb_max_set
408 b_offset = b_start + ncoset(lj - 1)
409 ncob = nco(lj)
410 prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
411 expj = 0.25_dp*real(2*lj + 3, dp)
412 normj = 1.0_dp/(prefj*zetj**expj)
413
414 !Loop over ECP angular momentum
415 DO lk = 0, lmax_ecp
416 tmp(1:ncoa*ncob) = 0.0_dp
417 !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
418 !the 1/norm coefficients for PGFi and PGFj
419 CALL libgrpp_type2_integrals(ra, li, 1, [normi], [zeti], &
420 rb, lj, 1, [normj], [zetj], &
421 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
422 coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
423
424 !note: tmp array is in C row-major ordering
425 DO j = 1, ncob
426 DO i = 1, ncoa
427 vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
428 END DO
429 END DO
430
431 IF (calc_forces) THEN
432
433 tmpx(1:ncoa*ncob) = 0.0_dp
434 tmpy(1:ncoa*ncob) = 0.0_dp
435 tmpz(1:ncoa*ncob) = 0.0_dp
436
437 !force wrt to atomic position A
438 CALL libgrpp_type2_integrals_gradient(ra, li, 1, [normi], [zeti], &
439 rb, lj, 1, [normj], [zetj], &
440 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
441 coeffs_ecp(:, lk), alpha_ecp(:, lk), ra, &
442 tmpx, tmpy, tmpz)
443
444 !note: tmp array is in C row-major ordering
445 !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
446 DO j = 1, ncob
447 DO i = 1, ncoa
448 force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
449 force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
450 force_a(3) = force_a(3) + fac_a*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
451 END DO
452 END DO
453
454 tmpx(1:ncoa*ncob) = 0.0_dp
455 tmpy(1:ncoa*ncob) = 0.0_dp
456 tmpz(1:ncoa*ncob) = 0.0_dp
457
458 !force wrt to atomic position B
459 CALL libgrpp_type2_integrals_gradient(ra, li, 1, [normi], [zeti], &
460 rb, lj, 1, [normj], [zetj], &
461 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
462 coeffs_ecp(:, lk), alpha_ecp(:, lk), rb, &
463 tmpx, tmpy, tmpz)
464 !note: tmp array is in C row-major ordering
465 !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
466 DO j = 1, ncob
467 DO i = 1, ncoa
468 force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
469 force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
470 force_b(3) = force_b(3) + fac_b*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
471 END DO
472 END DO
473
474 END IF !calc_forces
475
476 END DO !lk
477
478 END DO !lj
479 END DO !li
480
481 END DO !jpgf
482 END DO !ipgf
483
484#else
485
486 mark_used(la_max_set)
487 mark_used(la_min_set)
488 mark_used(npgfa)
489 mark_used(rpgfa)
490 mark_used(zeta)
491 mark_used(lb_max_set)
492 mark_used(lb_min_set)
493 mark_used(npgfb)
494 mark_used(rpgfb)
495 mark_used(zetb)
496 mark_used(lmax_ecp)
497 mark_used(npot_ecp)
498 mark_used(alpha_ecp)
499 mark_used(coeffs_ecp)
500 mark_used(nrpot_ecp)
501 mark_used(rpgfc)
502 mark_used(rab)
503 mark_used(dab)
504 mark_used(rac)
505 mark_used(dac)
506 mark_used(dbc)
507 mark_used(vab)
508 mark_used(pab)
509 mark_used(force_a)
510 mark_used(force_b)
511
512 cpabort("Please compile CP2K with libgrpp support for calculations with ECPs")
513#endif
514
515 END SUBROUTINE libgrpp_semilocal_integrals
516
517! **************************************************************************************************
518!> \brief Reference local ECP force routine using l+-1 integrals. No call is made to the numerically
519!> unstable gradient routine of libgrpp. Calculates both the integrals and the forces.
520!> \param la_max_set ...
521!> \param la_min_set ...
522!> \param npgfa ...
523!> \param rpgfa ...
524!> \param zeta ...
525!> \param lb_max_set ...
526!> \param lb_min_set ...
527!> \param npgfb ...
528!> \param rpgfb ...
529!> \param zetb ...
530!> \param npot_ecp ...
531!> \param alpha_ecp ...
532!> \param coeffs_ecp ...
533!> \param nrpot_ecp ...
534!> \param rpgfc ...
535!> \param rab ...
536!> \param dab ...
537!> \param rac ...
538!> \param dac ...
539!> \param dbc ...
540!> \param vab ...
541!> \param pab ...
542!> \param force_a ...
543!> \param force_b ...
544!> \note: this is a reference routine, which has no reason to be used once the libgrpp gradients
545!> become numerically stable
546! **************************************************************************************************
547 SUBROUTINE libgrpp_local_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
548 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
549 npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
550 rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
551
552 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
553 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
554 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
555 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
556 INTEGER, INTENT(IN) :: npot_ecp
557 REAL(kind=dp), DIMENSION(1:npot_ecp), INTENT(IN) :: alpha_ecp, coeffs_ecp
558 INTEGER, DIMENSION(1:npot_ecp), INTENT(IN) :: nrpot_ecp
559 REAL(kind=dp), INTENT(IN) :: rpgfc
560 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
561 REAL(kind=dp), INTENT(IN) :: dab
562 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
563 REAL(kind=dp), INTENT(IN) :: dac
564 REAL(kind=dp), INTENT(IN) :: dbc
565 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
566 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pab
567 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: force_a, force_b
568
569#if defined(__LIBGRPP)
570 INTEGER :: a_offset, a_start, b_offset, b_start, i, &
571 ipgf, j, jpgf, li, lj, ncoa, ncob, a_offset_f, &
572 b_offset_f, a_start_f, b_start_f
573 REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
574 zeti, zetj
575 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
576 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vab_f, tmpx, tmpy, tmpz
577 REAL(dp), DIMENSION(3) :: ra, rb, rc
578
579 !Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
580 ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
581 vab_f(:, :) = 0.0_dp
582
583 !libgrpp requires absolute positions, not relative ones
584 ra(:) = 0.0_dp
585 rb(:) = rab(:)
586 rc(:) = rac(:)
587
588 ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
589
590 DO ipgf = 1, npgfa
591 IF (rpgfa(ipgf) + rpgfc < dac) cycle
592 zeti = zeta(ipgf)
593 a_start = (ipgf - 1)*ncoset(la_max_set)
594 a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
595
596 DO jpgf = 1, npgfb
597 IF (rpgfb(jpgf) + rpgfc < dbc) cycle
598 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
599 zetj = zetb(jpgf)
600 b_start = (jpgf - 1)*ncoset(lb_max_set)
601 b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
602
603 DO li = max(0, la_min_set - 1), la_max_set + 1
604 a_offset = a_start + ncoset(li - 1)
605 a_offset_f = a_start_f + ncoset(li - 1)
606 ncoa = nco(li)
607 prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
608 expi = 0.25_dp*real(2*li + 3, dp)
609 normi = 1.0_dp/(prefi*zeti**expi)
610
611 DO lj = max(0, lb_min_set - 1), lb_max_set + 1
612 b_offset = b_start + ncoset(lj - 1)
613 b_offset_f = b_start_f + ncoset(lj - 1)
614 ncob = nco(lj)
615 prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
616 expj = 0.25_dp*real(2*lj + 3, dp)
617 normj = 1.0_dp/(prefj*zetj**expj)
618
619 tmp(1:ncoa*ncob) = 0.0_dp
620 !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
621 !the 1/norm coefficients for PGFi and PGFj
622 CALL libgrpp_type1_integrals(ra, li, 1, [normi], [zeti], &
623 rb, lj, 1, [normj], [zetj], &
624 rc, [npot_ecp], nrpot_ecp, &
625 coeffs_ecp, alpha_ecp, tmp)
626
627 !the l+-1 integrals for gradient calculation
628 DO j = 1, ncob
629 DO i = 1, ncoa
630 vab_f(a_offset_f + i, b_offset_f + j) = &
631 vab_f(a_offset_f + i, b_offset_f + j) + tmp((i - 1)*ncob + j)
632 END DO
633 END DO
634
635 !the actual integrals
636 IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
637 DO j = 1, ncob
638 DO i = 1, ncoa
639 vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
640 END DO
641 END DO
642 END IF
643
644 END DO !lj
645 END DO !li
646
647 END DO !jpgf
648 END DO !ipgf
649
650 ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
651 ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
652 ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
653
654 !Derivative wrt to center A
655 tmpx(:, :) = 0.0_dp
656 tmpy(:, :) = 0.0_dp
657 tmpz(:, :) = 0.0_dp
658 CALL dabdr(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, rpgfb, lb_min_set, &
659 dab, vab_f, tmpx, tmpy, tmpz)
660 DO j = 1, npgfb*ncoset(lb_max_set)
661 DO i = 1, npgfa*ncoset(la_max_set)
662 force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
663 force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
664 force_a(3) = force_a(3) + tmpz(i, j)*pab(i, j)
665 END DO
666 END DO
667
668 !Derivative wrt to center B
669 tmpx(:, :) = 0.0_dp
670 tmpy(:, :) = 0.0_dp
671 tmpz(:, :) = 0.0_dp
672 CALL adbdr(la_max_set, npgfa, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
673 dab, vab_f, tmpx, tmpy, tmpz)
674 DO j = 1, npgfb*ncoset(lb_max_set)
675 DO i = 1, npgfa*ncoset(la_max_set)
676 force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
677 force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
678 force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
679 END DO
680 END DO
681 DEALLOCATE (tmpx, tmpy, tmpz)
682
683#else
684
685 mark_used(la_max_set)
686 mark_used(la_min_set)
687 mark_used(npgfa)
688 mark_used(rpgfa)
689 mark_used(zeta)
690 mark_used(lb_max_set)
691 mark_used(lb_min_set)
692 mark_used(npgfb)
693 mark_used(rpgfb)
694 mark_used(zetb)
695 mark_used(npot_ecp)
696 mark_used(alpha_ecp)
697 mark_used(coeffs_ecp)
698 mark_used(nrpot_ecp)
699 mark_used(rpgfc)
700 mark_used(rab)
701 mark_used(dab)
702 mark_used(rac)
703 mark_used(dac)
704 mark_used(dbc)
705 mark_used(pab)
706 mark_used(vab)
707 mark_used(force_a)
708 mark_used(force_b)
709
710 cpabort("Please compile CP2K with libgrpp support for calculations with ECPs")
711#endif
712
713 END SUBROUTINE libgrpp_local_forces_ref
714
715! **************************************************************************************************
716!> \brief Reference semi-local ECP forces using l+-1 integrals. No call is made to the numerically
717!> unstable gradient routine of libgrpp. Calculates both the integrals and the forces.
718!> \param la_max_set ...
719!> \param la_min_set ...
720!> \param npgfa ...
721!> \param rpgfa ...
722!> \param zeta ...
723!> \param lb_max_set ...
724!> \param lb_min_set ...
725!> \param npgfb ...
726!> \param rpgfb ...
727!> \param zetb ...
728!> \param lmax_ecp ...
729!> \param npot_ecp ...
730!> \param alpha_ecp ...
731!> \param coeffs_ecp ...
732!> \param nrpot_ecp ...
733!> \param rpgfc ...
734!> \param rab ...
735!> \param dab ...
736!> \param rac ...
737!> \param dac ...
738!> \param dbc ...
739!> \param vab ...
740!> \param pab ...
741!> \param force_a ...
742!> \param force_b ...
743!> \note: this is a reference routine, which has no reason to be used once the libgrpp gradients
744!> become numerically stable
745! **************************************************************************************************
746 SUBROUTINE libgrpp_semilocal_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
747 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
748 lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
749 rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
750
751 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
752 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
753 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
754 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
755 INTEGER, INTENT(IN) :: lmax_ecp
756 INTEGER, DIMENSION(0:10), INTENT(IN) :: npot_ecp
757 REAL(kind=dp), DIMENSION(1:15, 0:10), INTENT(IN) :: alpha_ecp, coeffs_ecp
758 INTEGER, DIMENSION(1:15, 0:10), INTENT(IN) :: nrpot_ecp
759 REAL(kind=dp), INTENT(IN) :: rpgfc
760 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
761 REAL(kind=dp), INTENT(IN) :: dab
762 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
763 REAL(kind=dp), INTENT(IN) :: dac
764 REAL(kind=dp), INTENT(IN) :: dbc
765 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
766 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pab
767 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: force_a, force_b
768
769#if defined(__LIBGRPP)
770 INTEGER :: a_offset, a_start, b_offset, b_start, i, &
771 ipgf, j, jpgf, li, lj, lk, ncoa, ncob, &
772 a_start_f, b_start_f, a_offset_f, b_offset_f
773 REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
774 zeti, zetj
775 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
776 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vab_f, tmpx, tmpy, tmpz
777 REAL(dp), DIMENSION(3) :: ra, rb, rc
778
779 !Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
780 ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
781 vab_f(:, :) = 0.0_dp
782
783 !libgrpp requires absolute positions, not relative ones
784 ra(:) = 0.0_dp
785 rb(:) = rab(:)
786 rc(:) = rac(:)
787
788 ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
789
790 DO ipgf = 1, npgfa
791 IF (rpgfa(ipgf) + rpgfc < dac) cycle
792 zeti = zeta(ipgf)
793 a_start = (ipgf - 1)*ncoset(la_max_set)
794 a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
795
796 DO jpgf = 1, npgfb
797 IF (rpgfb(jpgf) + rpgfc < dbc) cycle
798 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
799 zetj = zetb(jpgf)
800 b_start = (jpgf - 1)*ncoset(lb_max_set)
801 b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
802
803 DO li = max(0, la_min_set - 1), la_max_set + 1
804 a_offset = a_start + ncoset(li - 1)
805 a_offset_f = a_start_f + ncoset(li - 1)
806 ncoa = nco(li)
807 prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
808 expi = 0.25_dp*real(2*li + 3, dp)
809 normi = 1.0_dp/(prefi*zeti**expi)
810
811 DO lj = max(0, lb_min_set - 1), lb_max_set + 1
812 b_offset = b_start + ncoset(lj - 1)
813 b_offset_f = b_start_f + ncoset(lj - 1)
814 ncob = nco(lj)
815 prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
816 expj = 0.25_dp*real(2*lj + 3, dp)
817 normj = 1.0_dp/(prefj*zetj**expj)
818
819 !Loop over ECP angular momentum
820 DO lk = 0, lmax_ecp
821 tmp(1:ncoa*ncob) = 0.0_dp
822 !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
823 !the 1/norm coefficients for PGFi and PGFj
824 CALL libgrpp_type2_integrals(ra, li, 1, [normi], [zeti], &
825 rb, lj, 1, [normj], [zetj], &
826 rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
827 coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
828
829 !the l+-1 integrals for gradient calculation
830 DO j = 1, ncob
831 DO i = 1, ncoa
832 vab_f(a_offset_f + i, b_offset_f + j) = &
833 vab_f(a_offset_f + i, b_offset_f + j) + tmp((i - 1)*ncob + j)
834 END DO
835 END DO
836
837 !the actual integrals
838 IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
839 DO j = 1, ncob
840 DO i = 1, ncoa
841 vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
842 END DO
843 END DO
844 END IF
845
846 END DO !lk
847
848 END DO !lj
849 END DO !li
850
851 END DO !jpgf
852 END DO !ipgf
853
854 ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
855 ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
856 ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
857
858 !Derivative wrt to center A
859 tmpx(:, :) = 0.0_dp
860 tmpy(:, :) = 0.0_dp
861 tmpz(:, :) = 0.0_dp
862 CALL dabdr(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, rpgfb, lb_min_set, &
863 0.0_dp, vab_f, tmpx, tmpy, tmpz)
864 DO j = 1, npgfb*ncoset(lb_max_set)
865 DO i = 1, npgfa*ncoset(la_max_set)
866 force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
867 force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
868 force_a(3) = force_a(3) + tmpz(i, j)*pab(i, j)
869 END DO
870 END DO
871
872 !Derivative wrt to center B
873 tmpx(:, :) = 0.0_dp
874 tmpy(:, :) = 0.0_dp
875 tmpz(:, :) = 0.0_dp
876 CALL adbdr(la_max_set, npgfa, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
877 0.0_dp, vab_f, tmpx, tmpy, tmpz)
878 DO j = 1, npgfb*ncoset(lb_max_set)
879 DO i = 1, npgfa*ncoset(la_max_set)
880 force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
881 force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
882 force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
883 END DO
884 END DO
885 DEALLOCATE (tmpx, tmpy, tmpz)
886
887#else
888
889 mark_used(la_max_set)
890 mark_used(la_min_set)
891 mark_used(npgfa)
892 mark_used(rpgfa)
893 mark_used(zeta)
894 mark_used(lb_max_set)
895 mark_used(lb_min_set)
896 mark_used(npgfb)
897 mark_used(rpgfb)
898 mark_used(zetb)
899 mark_used(lmax_ecp)
900 mark_used(npot_ecp)
901 mark_used(alpha_ecp)
902 mark_used(coeffs_ecp)
903 mark_used(nrpot_ecp)
904 mark_used(rpgfc)
905 mark_used(rab)
906 mark_used(dab)
907 mark_used(rac)
908 mark_used(dac)
909 mark_used(dbc)
910 mark_used(pab)
911 mark_used(vab)
912 mark_used(force_a)
913 mark_used(force_b)
914
915 cpabort("Please compile CP2K with libgrpp support for calculations with ECPs")
916#endif
917
918 END SUBROUTINE libgrpp_semilocal_forces_ref
919
920END MODULE libgrpp_integrals
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:303
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:551
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:750
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