Line data Source code
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 :
12 : MODULE 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 :
32 : CONTAINS
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 0 : SUBROUTINE libgrpp_local_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
62 0 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
63 0 : npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
64 0 : 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 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmpx, tmpy, tmpz
92 : REAL(dp), DIMENSION(3) :: ra, rb, rc
93 :
94 0 : CALL libgrpp_init()
95 :
96 0 : calc_forces = .FALSE.
97 0 : 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 0 : "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 0 : mindist = 1.0E-6_dp
114 : !If ra != rb != rc
115 0 : 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 0 : 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 0 : 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 0 : 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 0 : calc_forces = .FALSE.
137 : END IF
138 : END IF
139 :
140 : !libgrpp requires absolute positions, not relative ones
141 0 : ra(:) = 0.0_dp
142 0 : rb(:) = rab(:)
143 0 : rc(:) = rac(:)
144 :
145 0 : ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
146 0 : IF (calc_forces) THEN
147 0 : ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
148 0 : ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
149 0 : ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
150 : END IF
151 :
152 0 : DO ipgf = 1, npgfa
153 0 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
154 0 : zeti = zeta(ipgf)
155 0 : a_start = (ipgf - 1)*ncoset(la_max_set)
156 :
157 0 : DO jpgf = 1, npgfb
158 0 : IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
159 0 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
160 0 : zetj = zetb(jpgf)
161 0 : b_start = (jpgf - 1)*ncoset(lb_max_set)
162 :
163 0 : DO li = la_min_set, la_max_set
164 0 : a_offset = a_start + ncoset(li - 1)
165 0 : ncoa = nco(li)
166 0 : prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
167 0 : expi = 0.25_dp*REAL(2*li + 3, dp)
168 0 : normi = 1.0_dp/(prefi*zeti**expi)
169 :
170 0 : DO lj = lb_min_set, lb_max_set
171 0 : b_offset = b_start + ncoset(lj - 1)
172 0 : ncob = nco(lj)
173 0 : prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
174 0 : expj = 0.25_dp*REAL(2*lj + 3, dp)
175 0 : normj = 1.0_dp/(prefj*zetj**expj)
176 :
177 0 : 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 0 : coeffs_ecp, alpha_ecp, tmp)
184 :
185 : !note: tmp array is in C row-major ordering
186 0 : DO j = 1, ncob
187 0 : DO i = 1, ncoa
188 0 : 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 0 : IF (calc_forces) THEN
193 0 : tmpx(1:ncoa*ncob) = 0.0_dp
194 0 : tmpy(1:ncoa*ncob) = 0.0_dp
195 0 : 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 0 : 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 0 : DO j = 1, ncob
207 0 : DO i = 1, ncoa
208 0 : force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
209 0 : force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
210 0 : 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 0 : tmpx(1:ncoa*ncob) = 0.0_dp
215 0 : tmpy(1:ncoa*ncob) = 0.0_dp
216 0 : 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 0 : 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 0 : DO j = 1, ncob
228 0 : DO i = 1, ncoa
229 0 : force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
230 0 : force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
231 0 : 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 0 : 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 8980 : SUBROUTINE libgrpp_semilocal_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
302 8980 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
303 : lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
304 17960 : 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 8980 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmpx, tmpz, tmpy
333 : REAL(dp), DIMENSION(3) :: ra, rb, rc
334 :
335 8980 : CALL libgrpp_init()
336 :
337 8980 : calc_forces = .FALSE.
338 8980 : 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 0 : "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 0 : mindist = 1.0E-6_dp
355 : !If ra != rb != rc
356 0 : 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 0 : 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 0 : 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 0 : 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 0 : calc_forces = .FALSE.
378 : END IF
379 : END IF
380 :
381 : !libgrpp requires absolute positions, not relative ones
382 8980 : ra(:) = 0.0_dp
383 8980 : rb(:) = rab(:)
384 8980 : rc(:) = rac(:)
385 :
386 26940 : ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
387 8980 : IF (calc_forces) THEN
388 0 : ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
389 0 : ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
390 0 : ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
391 : END IF
392 :
393 21618 : DO ipgf = 1, npgfa
394 12638 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
395 10216 : zeti = zeta(ipgf)
396 10216 : a_start = (ipgf - 1)*ncoset(la_max_set)
397 :
398 34094 : DO jpgf = 1, npgfb
399 14898 : IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
400 12164 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
401 11380 : zetj = zetb(jpgf)
402 11380 : b_start = (jpgf - 1)*ncoset(lb_max_set)
403 :
404 35558 : DO li = la_min_set, la_max_set
405 11540 : a_offset = a_start + ncoset(li - 1)
406 11540 : ncoa = nco(li)
407 11540 : prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
408 11540 : expi = 0.25_dp*REAL(2*li + 3, dp)
409 11540 : normi = 1.0_dp/(prefi*zeti**expi)
410 :
411 38394 : DO lj = lb_min_set, lb_max_set
412 11956 : b_offset = b_start + ncoset(lj - 1)
413 11956 : ncob = nco(lj)
414 11956 : prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
415 11956 : expj = 0.25_dp*REAL(2*lj + 3, dp)
416 11956 : normj = 1.0_dp/(prefj*zetj**expj)
417 :
418 : !Loop over ECP angular momentum
419 76096 : DO lk = 0, lmax_ecp
420 466322 : 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 315600 : coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
427 :
428 : !note: tmp array is in C row-major ordering
429 200097 : DO j = 1, ncob
430 613819 : DO i = 1, ncoa
431 561219 : 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 64556 : IF (calc_forces) THEN
436 :
437 0 : tmpx(1:ncoa*ncob) = 0.0_dp
438 0 : tmpy(1:ncoa*ncob) = 0.0_dp
439 0 : 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 0 : 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 0 : DO j = 1, ncob
451 0 : DO i = 1, ncoa
452 0 : force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
453 0 : force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
454 0 : 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 0 : tmpx(1:ncoa*ncob) = 0.0_dp
459 0 : tmpy(1:ncoa*ncob) = 0.0_dp
460 0 : 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 0 : 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 0 : DO j = 1, ncob
471 0 : DO i = 1, ncoa
472 0 : force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
473 0 : force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
474 0 : 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 8980 : 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 0 : SUBROUTINE libgrpp_local_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
552 0 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
553 0 : npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
554 0 : 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 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
580 0 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vab_f, tmpx, tmpy, tmpz
581 : REAL(dp), DIMENSION(3) :: ra, rb, rc
582 :
583 0 : CALL libgrpp_init()
584 :
585 : !Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
586 0 : ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
587 0 : vab_f(:, :) = 0.0_dp
588 :
589 : !libgrpp requires absolute positions, not relative ones
590 0 : ra(:) = 0.0_dp
591 0 : rb(:) = rab(:)
592 0 : rc(:) = rac(:)
593 :
594 0 : ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
595 :
596 0 : DO ipgf = 1, npgfa
597 0 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
598 0 : zeti = zeta(ipgf)
599 0 : a_start = (ipgf - 1)*ncoset(la_max_set)
600 0 : a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
601 :
602 0 : DO jpgf = 1, npgfb
603 0 : IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
604 0 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
605 0 : zetj = zetb(jpgf)
606 0 : b_start = (jpgf - 1)*ncoset(lb_max_set)
607 0 : b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
608 :
609 0 : DO li = MAX(0, la_min_set - 1), la_max_set + 1
610 0 : a_offset = a_start + ncoset(li - 1)
611 0 : a_offset_f = a_start_f + ncoset(li - 1)
612 0 : ncoa = nco(li)
613 0 : prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
614 0 : expi = 0.25_dp*REAL(2*li + 3, dp)
615 0 : normi = 1.0_dp/(prefi*zeti**expi)
616 :
617 0 : DO lj = MAX(0, lb_min_set - 1), lb_max_set + 1
618 0 : b_offset = b_start + ncoset(lj - 1)
619 0 : b_offset_f = b_start_f + ncoset(lj - 1)
620 0 : ncob = nco(lj)
621 0 : prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
622 0 : expj = 0.25_dp*REAL(2*lj + 3, dp)
623 0 : normj = 1.0_dp/(prefj*zetj**expj)
624 :
625 0 : 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 0 : coeffs_ecp, alpha_ecp, tmp)
632 :
633 : !the l+-1 integrals for gradient calculation
634 0 : DO j = 1, ncob
635 0 : DO i = 1, ncoa
636 : vab_f(a_offset_f + i, b_offset_f + j) = &
637 0 : 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 0 : IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
643 0 : DO j = 1, ncob
644 0 : DO i = 1, ncoa
645 0 : 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 0 : ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
657 0 : ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
658 0 : ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
659 :
660 : !Derivative wrt to center A
661 0 : tmpx(:, :) = 0.0_dp
662 0 : tmpy(:, :) = 0.0_dp
663 0 : 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 0 : dab, vab_f, tmpx, tmpy, tmpz)
666 0 : DO j = 1, npgfb*ncoset(lb_max_set)
667 0 : DO i = 1, npgfa*ncoset(la_max_set)
668 0 : force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
669 0 : force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
670 0 : 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 0 : tmpx(:, :) = 0.0_dp
676 0 : tmpy(:, :) = 0.0_dp
677 0 : 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 0 : dab, vab_f, tmpx, tmpy, tmpz)
680 0 : DO j = 1, npgfb*ncoset(lb_max_set)
681 0 : DO i = 1, npgfa*ncoset(la_max_set)
682 0 : force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
683 0 : force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
684 0 : force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
685 : END DO
686 : END DO
687 0 : 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 0 : 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 2400 : SUBROUTINE libgrpp_semilocal_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
753 2400 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
754 : lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
755 2400 : 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 2400 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
782 2400 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vab_f, tmpx, tmpy, tmpz
783 : REAL(dp), DIMENSION(3) :: ra, rb, rc
784 :
785 2400 : CALL libgrpp_init()
786 :
787 : !Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
788 9600 : ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
789 427592 : vab_f(:, :) = 0.0_dp
790 :
791 : !libgrpp requires absolute positions, not relative ones
792 2400 : ra(:) = 0.0_dp
793 2400 : rb(:) = rab(:)
794 2400 : rc(:) = rac(:)
795 :
796 7200 : ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
797 :
798 5976 : DO ipgf = 1, npgfa
799 3576 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
800 2794 : zeti = zeta(ipgf)
801 2794 : a_start = (ipgf - 1)*ncoset(la_max_set)
802 2794 : a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
803 :
804 9492 : DO jpgf = 1, npgfb
805 4298 : IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
806 3364 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
807 2980 : zetj = zetb(jpgf)
808 2980 : b_start = (jpgf - 1)*ncoset(lb_max_set)
809 2980 : b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
810 :
811 14377 : DO li = MAX(0, la_min_set - 1), la_max_set + 1
812 7821 : a_offset = a_start + ncoset(li - 1)
813 7821 : a_offset_f = a_start_f + ncoset(li - 1)
814 7821 : ncoa = nco(li)
815 7821 : prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
816 7821 : expi = 0.25_dp*REAL(2*li + 3, dp)
817 7821 : normi = 1.0_dp/(prefi*zeti**expi)
818 :
819 32773 : DO lj = MAX(0, lb_min_set - 1), lb_max_set + 1
820 20654 : b_offset = b_start + ncoset(lj - 1)
821 20654 : b_offset_f = b_start_f + ncoset(lj - 1)
822 20654 : ncob = nco(lj)
823 20654 : prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
824 20654 : expj = 0.25_dp*REAL(2*lj + 3, dp)
825 20654 : normj = 1.0_dp/(prefj*zetj**expj)
826 :
827 : !Loop over ECP angular momentum
828 125507 : DO lk = 0, lmax_ecp
829 1210945 : 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 582192 : coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
836 :
837 : !the l+-1 integrals for gradient calculation
838 425907 : DO j = 1, ncob
839 1539820 : DO i = 1, ncoa
840 : vab_f(a_offset_f + i, b_offset_f + j) = &
841 1442788 : 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 117686 : IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
847 51400 : DO j = 1, ncob
848 149614 : DO i = 1, ncoa
849 135465 : 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 9600 : ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
863 7200 : ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
864 7200 : ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
865 :
866 : !Derivative wrt to center A
867 75627 : tmpx(:, :) = 0.0_dp
868 75627 : tmpy(:, :) = 0.0_dp
869 75627 : 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 2400 : 0.0_dp, vab_f, tmpx, tmpy, tmpz)
872 14583 : DO j = 1, npgfb*ncoset(lb_max_set)
873 75627 : DO i = 1, npgfa*ncoset(la_max_set)
874 61044 : force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
875 61044 : force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
876 73227 : 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 75627 : tmpx(:, :) = 0.0_dp
882 75627 : tmpy(:, :) = 0.0_dp
883 75627 : 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 2400 : 0.0_dp, vab_f, tmpx, tmpy, tmpz)
886 14583 : DO j = 1, npgfb*ncoset(lb_max_set)
887 75627 : DO i = 1, npgfa*ncoset(la_max_set)
888 61044 : force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
889 61044 : force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
890 73227 : force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
891 : END DO
892 : END DO
893 2400 : 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 2400 : END SUBROUTINE libgrpp_semilocal_forces_ref
927 :
928 : END MODULE libgrpp_integrals
|
