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