Line data Source code
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 :
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_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 96 : SUBROUTINE libgrpp_local_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
62 96 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
63 96 : npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
64 192 : 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 96 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmpx, tmpy, tmpz
92 : REAL(dp), DIMENSION(3) :: ra, rb, rc
93 :
94 96 : calc_forces = .FALSE.
95 96 : 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 0 : "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 0 : mindist = 1.0E-6_dp
112 : !If ra != rb != rc
113 0 : 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 0 : 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 0 : 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 0 : 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 0 : calc_forces = .FALSE.
135 : END IF
136 : END IF
137 :
138 : !libgrpp requires absolute positions, not relative ones
139 96 : ra(:) = 0.0_dp
140 96 : rb(:) = rab(:)
141 96 : rc(:) = rac(:)
142 :
143 288 : ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
144 96 : IF (calc_forces) THEN
145 0 : ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
146 0 : ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
147 0 : ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
148 : END IF
149 :
150 384 : DO ipgf = 1, npgfa
151 288 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
152 288 : zeti = zeta(ipgf)
153 288 : a_start = (ipgf - 1)*ncoset(la_max_set)
154 :
155 1248 : DO jpgf = 1, npgfb
156 864 : IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
157 864 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
158 864 : zetj = zetb(jpgf)
159 864 : b_start = (jpgf - 1)*ncoset(lb_max_set)
160 :
161 2016 : DO li = la_min_set, la_max_set
162 864 : a_offset = a_start + ncoset(li - 1)
163 864 : ncoa = nco(li)
164 864 : prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
165 864 : expi = 0.25_dp*REAL(2*li + 3, dp)
166 864 : normi = 1.0_dp/(prefi*zeti**expi)
167 :
168 2592 : DO lj = lb_min_set, lb_max_set
169 864 : b_offset = b_start + ncoset(lj - 1)
170 864 : ncob = nco(lj)
171 864 : prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
172 864 : expj = 0.25_dp*REAL(2*lj + 3, dp)
173 864 : normj = 1.0_dp/(prefj*zetj**expj)
174 :
175 4320 : 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 5184 : coeffs_ecp, alpha_ecp, tmp)
182 :
183 : !note: tmp array is in C row-major ordering
184 2592 : DO j = 1, ncob
185 6048 : DO i = 1, ncoa
186 5184 : 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 1728 : IF (calc_forces) THEN
191 0 : tmpx(1:ncoa*ncob) = 0.0_dp
192 0 : tmpy(1:ncoa*ncob) = 0.0_dp
193 0 : 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 0 : 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 0 : DO j = 1, ncob
205 0 : DO i = 1, ncoa
206 0 : force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
207 0 : force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
208 0 : 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 0 : tmpx(1:ncoa*ncob) = 0.0_dp
213 0 : tmpy(1:ncoa*ncob) = 0.0_dp
214 0 : 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 0 : 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 0 : DO j = 1, ncob
226 0 : DO i = 1, ncoa
227 0 : force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
228 0 : force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
229 0 : 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 96 : 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 9368 : SUBROUTINE libgrpp_semilocal_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
300 9368 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
301 : lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
302 18736 : 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 9368 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmpx, tmpz, tmpy
331 : REAL(dp), DIMENSION(3) :: ra, rb, rc
332 :
333 9368 : calc_forces = .FALSE.
334 9368 : 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 0 : "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 0 : mindist = 1.0E-6_dp
351 : !If ra != rb != rc
352 0 : 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 0 : 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 0 : 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 0 : 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 0 : calc_forces = .FALSE.
374 : END IF
375 : END IF
376 :
377 : !libgrpp requires absolute positions, not relative ones
378 9368 : ra(:) = 0.0_dp
379 9368 : rb(:) = rab(:)
380 9368 : rc(:) = rac(:)
381 :
382 28104 : ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
383 9368 : IF (calc_forces) THEN
384 0 : ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
385 0 : ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
386 0 : ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
387 : END IF
388 :
389 22865 : DO ipgf = 1, npgfa
390 13497 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
391 11075 : zeti = zeta(ipgf)
392 11075 : a_start = (ipgf - 1)*ncoset(la_max_set)
393 :
394 37277 : DO jpgf = 1, npgfb
395 16834 : IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
396 14100 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
397 13188 : zetj = zetb(jpgf)
398 13188 : b_start = (jpgf - 1)*ncoset(lb_max_set)
399 :
400 39873 : DO li = la_min_set, la_max_set
401 13188 : a_offset = a_start + ncoset(li - 1)
402 13188 : ncoa = nco(li)
403 13188 : prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
404 13188 : expi = 0.25_dp*REAL(2*li + 3, dp)
405 13188 : normi = 1.0_dp/(prefi*zeti**expi)
406 :
407 43210 : DO lj = lb_min_set, lb_max_set
408 13188 : b_offset = b_start + ncoset(lj - 1)
409 13188 : ncob = nco(lj)
410 13188 : prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
411 13188 : expj = 0.25_dp*REAL(2*lj + 3, dp)
412 13188 : normj = 1.0_dp/(prefj*zetj**expj)
413 :
414 : !Loop over ECP angular momentum
415 84268 : DO lk = 0, lmax_ecp
416 500450 : 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 347352 : coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
423 :
424 : !note: tmp array is in C row-major ordering
425 218184 : DO j = 1, ncob
426 660742 : DO i = 1, ncoa
427 602850 : 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 71080 : IF (calc_forces) THEN
432 :
433 0 : tmpx(1:ncoa*ncob) = 0.0_dp
434 0 : tmpy(1:ncoa*ncob) = 0.0_dp
435 0 : 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 0 : 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 0 : DO j = 1, ncob
447 0 : DO i = 1, ncoa
448 0 : force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
449 0 : force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
450 0 : 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 0 : tmpx(1:ncoa*ncob) = 0.0_dp
455 0 : tmpy(1:ncoa*ncob) = 0.0_dp
456 0 : 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 0 : 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 0 : DO j = 1, ncob
467 0 : DO i = 1, ncoa
468 0 : force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
469 0 : force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
470 0 : 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 9368 : 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 24 : SUBROUTINE libgrpp_local_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
548 24 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
549 24 : npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
550 24 : 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 24 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
576 24 : 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 96 : ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
581 11112 : vab_f(:, :) = 0.0_dp
582 :
583 : !libgrpp requires absolute positions, not relative ones
584 24 : ra(:) = 0.0_dp
585 24 : rb(:) = rab(:)
586 24 : rc(:) = rac(:)
587 :
588 72 : ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
589 :
590 96 : DO ipgf = 1, npgfa
591 72 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
592 72 : zeti = zeta(ipgf)
593 72 : a_start = (ipgf - 1)*ncoset(la_max_set)
594 72 : a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
595 :
596 312 : DO jpgf = 1, npgfb
597 216 : IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
598 216 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
599 216 : zetj = zetb(jpgf)
600 216 : b_start = (jpgf - 1)*ncoset(lb_max_set)
601 216 : b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
602 :
603 828 : DO li = MAX(0, la_min_set - 1), la_max_set + 1
604 540 : a_offset = a_start + ncoset(li - 1)
605 540 : a_offset_f = a_start_f + ncoset(li - 1)
606 540 : ncoa = nco(li)
607 540 : prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
608 540 : expi = 0.25_dp*REAL(2*li + 3, dp)
609 540 : normi = 1.0_dp/(prefi*zeti**expi)
610 :
611 2106 : DO lj = MAX(0, lb_min_set - 1), lb_max_set + 1
612 1350 : b_offset = b_start + ncoset(lj - 1)
613 1350 : b_offset_f = b_start_f + ncoset(lj - 1)
614 1350 : ncob = nco(lj)
615 1350 : prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
616 1350 : expj = 0.25_dp*REAL(2*lj + 3, dp)
617 1350 : normj = 1.0_dp/(prefj*zetj**expj)
618 :
619 11934 : 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 8100 : coeffs_ecp, alpha_ecp, tmp)
626 :
627 : !the l+-1 integrals for gradient calculation
628 5130 : DO j = 1, ncob
629 15714 : DO i = 1, ncoa
630 : vab_f(a_offset_f + i, b_offset_f + j) = &
631 14364 : 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 1890 : IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
637 648 : DO j = 1, ncob
638 1512 : DO i = 1, ncoa
639 1296 : 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 96 : ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
651 72 : ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
652 72 : ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
653 :
654 : !Derivative wrt to center A
655 1554 : tmpx(:, :) = 0.0_dp
656 1554 : tmpy(:, :) = 0.0_dp
657 1554 : 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 24 : dab, vab_f, tmpx, tmpy, tmpz)
660 204 : DO j = 1, npgfb*ncoset(lb_max_set)
661 1554 : DO i = 1, npgfa*ncoset(la_max_set)
662 1350 : force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
663 1350 : force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
664 1530 : 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 1554 : tmpx(:, :) = 0.0_dp
670 1554 : tmpy(:, :) = 0.0_dp
671 1554 : 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 24 : dab, vab_f, tmpx, tmpy, tmpz)
674 204 : DO j = 1, npgfb*ncoset(lb_max_set)
675 1554 : DO i = 1, npgfa*ncoset(la_max_set)
676 1350 : force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
677 1350 : force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
678 1530 : force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
679 : END DO
680 : END DO
681 24 : 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 24 : 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 2874 : SUBROUTINE libgrpp_semilocal_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
747 2874 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
748 : lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
749 2874 : 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 2874 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
776 2874 : 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 11496 : ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
781 576650 : vab_f(:, :) = 0.0_dp
782 :
783 : !libgrpp requires absolute positions, not relative ones
784 2874 : ra(:) = 0.0_dp
785 2874 : rb(:) = rab(:)
786 2874 : rc(:) = rac(:)
787 :
788 8622 : ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
789 :
790 7384 : DO ipgf = 1, npgfa
791 4510 : IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
792 3728 : zeti = zeta(ipgf)
793 3728 : a_start = (ipgf - 1)*ncoset(la_max_set)
794 3728 : a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
795 :
796 12826 : DO jpgf = 1, npgfb
797 6224 : IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
798 5290 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
799 4817 : zetj = zetb(jpgf)
800 4817 : b_start = (jpgf - 1)*ncoset(lb_max_set)
801 4817 : b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
802 :
803 21639 : DO li = MAX(0, la_min_set - 1), la_max_set + 1
804 12312 : a_offset = a_start + ncoset(li - 1)
805 12312 : a_offset_f = a_start_f + ncoset(li - 1)
806 12312 : ncoa = nco(li)
807 12312 : prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
808 12312 : expi = 0.25_dp*REAL(2*li + 3, dp)
809 12312 : normi = 1.0_dp/(prefi*zeti**expi)
810 :
811 50252 : DO lj = MAX(0, lb_min_set - 1), lb_max_set + 1
812 31716 : b_offset = b_start + ncoset(lj - 1)
813 31716 : b_offset_f = b_start_f + ncoset(lj - 1)
814 31716 : ncob = nco(lj)
815 31716 : prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
816 31716 : expj = 0.25_dp*REAL(2*lj + 3, dp)
817 31716 : normj = 1.0_dp/(prefj*zetj**expj)
818 :
819 : !Loop over ECP angular momentum
820 175473 : DO lk = 0, lmax_ecp
821 1652083 : 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 788670 : coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
828 :
829 : !the l+-1 integrals for gradient calculation
830 580299 : DO j = 1, ncob
831 2100937 : DO i = 1, ncoa
832 : vab_f(a_offset_f + i, b_offset_f + j) = &
833 1969492 : 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 163161 : IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
839 70924 : DO j = 1, ncob
840 203815 : DO i = 1, ncoa
841 184188 : 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 11496 : ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
855 8622 : ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
856 8622 : ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
857 :
858 : !Derivative wrt to center A
859 102609 : tmpx(:, :) = 0.0_dp
860 102609 : tmpy(:, :) = 0.0_dp
861 102609 : 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 2874 : 0.0_dp, vab_f, tmpx, tmpy, tmpz)
864 18477 : DO j = 1, npgfb*ncoset(lb_max_set)
865 102609 : DO i = 1, npgfa*ncoset(la_max_set)
866 84132 : force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
867 84132 : force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
868 99735 : 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 102609 : tmpx(:, :) = 0.0_dp
874 102609 : tmpy(:, :) = 0.0_dp
875 102609 : 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 2874 : 0.0_dp, vab_f, tmpx, tmpy, tmpz)
878 18477 : DO j = 1, npgfb*ncoset(lb_max_set)
879 102609 : DO i = 1, npgfa*ncoset(la_max_set)
880 84132 : force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
881 84132 : force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
882 99735 : force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
883 : END DO
884 : END DO
885 2874 : 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 2874 : END SUBROUTINE libgrpp_semilocal_forces_ref
919 :
920 : END MODULE libgrpp_integrals
|
