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