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 Rotationally invariant parametrization of Fock matrix.
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE pao_linpot_rotinv
13 : USE ai_overlap, ONLY: overlap_aab
14 : USE atomic_kind_types, ONLY: get_atomic_kind
15 : USE basis_set_types, ONLY: gto_basis_set_type
16 : USE cell_types, ONLY: cell_type,&
17 : pbc
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: gamma1
20 : USE mathlib, ONLY: multinomial
21 : USE orbital_pointers, ONLY: indco,&
22 : ncoset
23 : USE particle_types, ONLY: particle_type
24 : USE qs_environment_types, ONLY: get_qs_env,&
25 : qs_environment_type
26 : USE qs_kind_types, ONLY: get_qs_kind,&
27 : pao_potential_type,&
28 : qs_kind_type
29 : #include "./base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 :
33 : PRIVATE
34 :
35 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_linpot_rotinv'
36 :
37 : PUBLIC :: linpot_rotinv_count_terms, linpot_rotinv_calc_terms, linpot_rotinv_calc_forces
38 :
39 : CONTAINS
40 :
41 : ! **************************************************************************************************
42 : !> \brief Count number of terms for given atomic kind
43 : !> \param qs_env ...
44 : !> \param ikind ...
45 : !> \param nterms ...
46 : ! **************************************************************************************************
47 538 : SUBROUTINE linpot_rotinv_count_terms(qs_env, ikind, nterms)
48 : TYPE(qs_environment_type), POINTER :: qs_env
49 : INTEGER, INTENT(IN) :: ikind
50 : INTEGER, INTENT(OUT) :: nterms
51 :
52 : CHARACTER(len=*), PARAMETER :: routineN = 'linpot_rotinv_count_terms'
53 :
54 : INTEGER :: handle, ipot, iset, ishell, ishell_abs, &
55 : lmax, lmin, lpot, max_shell, &
56 : min_shell, npots, nshells, pot_maxl
57 538 : INTEGER, ALLOCATABLE, DIMENSION(:) :: shell_l
58 : TYPE(gto_basis_set_type), POINTER :: basis_set
59 538 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
60 538 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
61 :
62 538 : CALL timeset(routineN, handle)
63 :
64 538 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
65 538 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, pao_potentials=pao_potentials)
66 :
67 1132 : nshells = SUM(basis_set%nshell)
68 538 : npots = SIZE(pao_potentials)
69 :
70 538 : CPWARN_IF(npots == 0, "Found no PAO_POTENTIAL section")
71 :
72 : ! fill shell_l
73 1614 : ALLOCATE (shell_l(nshells))
74 1132 : DO iset = 1, basis_set%nset
75 2832 : DO ishell = 1, basis_set%nshell(iset)
76 1756 : ishell_abs = SUM(basis_set%nshell(1:iset - 1)) + ishell
77 2294 : shell_l(ishell_abs) = basis_set%l(ishell, iset)
78 : END DO
79 : END DO
80 :
81 538 : nterms = 0
82 :
83 : ! terms sensing neighboring atoms
84 1081 : DO ipot = 1, npots
85 543 : pot_maxl = pao_potentials(ipot)%maxl ! maxl is taken from central atom
86 543 : IF (pot_maxl < 0) THEN
87 0 : CPABORT("ROTINV parametrization requires non-negative PAO_POTENTIAL%MAXL")
88 : END IF
89 543 : IF (MOD(pot_maxl, 2) /= 0) THEN
90 0 : CPABORT("ROTINV parametrization requires even-numbered PAO_POTENTIAL%MAXL")
91 : END IF
92 2796 : DO max_shell = 1, nshells
93 5875 : DO min_shell = 1, max_shell
94 9903 : DO lpot = 0, pot_maxl, 2
95 4571 : lmin = shell_l(min_shell)
96 4571 : lmax = shell_l(max_shell)
97 4571 : IF (lmin == 0 .AND. lmax == 0) CYCLE ! covered by central terms
98 8188 : nterms = nterms + 1
99 : END DO
100 : END DO
101 : END DO
102 : END DO
103 :
104 : ! spherical symmetric terms on central atom
105 2238 : DO max_shell = 1, nshells
106 5825 : DO min_shell = 1, max_shell
107 3587 : IF (shell_l(min_shell) /= shell_l(max_shell)) CYCLE ! need quadratic block
108 5287 : nterms = nterms + 1
109 : END DO
110 : END DO
111 :
112 538 : CALL timestop(handle)
113 :
114 1076 : END SUBROUTINE linpot_rotinv_count_terms
115 :
116 : ! **************************************************************************************************
117 : !> \brief Calculate all potential terms of the rotinv parametrization
118 : !> \param qs_env ...
119 : !> \param iatom ...
120 : !> \param V_blocks ...
121 : ! **************************************************************************************************
122 195 : SUBROUTINE linpot_rotinv_calc_terms(qs_env, iatom, V_blocks)
123 : TYPE(qs_environment_type), POINTER :: qs_env
124 : INTEGER, INTENT(IN) :: iatom
125 : REAL(dp), DIMENSION(:, :, :), INTENT(OUT), TARGET :: V_blocks
126 :
127 : CHARACTER(len=*), PARAMETER :: routineN = 'linpot_rotinv_calc_terms'
128 :
129 : INTEGER :: handle, i, ic, ikind, ipot, iset, ishell, ishell_abs, jatom, jkind, jset, jshell, &
130 : jshell_abs, kterm, la1_max, la1_min, la2_max, la2_min, lb_max, lb_min, lpot, N, na1, na2, &
131 : natoms, nb, ncfga1, ncfga2, ncfgb, npgfa1, npgfa2, npgfb, npots, pot_maxl, sgfa1, sgfa2, &
132 : sgla1, sgla2
133 : REAL(dp) :: coeff, norm2, pot_beta, pot_weight, &
134 : rpgfa_max, tab
135 : REAL(dp), DIMENSION(3) :: Ra, Rab, Rb
136 195 : REAL(dp), DIMENSION(:), POINTER :: rpgfa1, rpgfa2, rpgfb, zeta1, zeta2, zetb
137 195 : REAL(dp), DIMENSION(:, :), POINTER :: T1, T2, V12, V21
138 195 : REAL(dp), DIMENSION(:, :, :), POINTER :: block_V_full, saab, saal
139 : TYPE(cell_type), POINTER :: cell
140 : TYPE(gto_basis_set_type), POINTER :: basis_set
141 195 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: ipao_potentials, jpao_potentials
142 195 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
143 195 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
144 :
145 195 : CALL timeset(routineN, handle)
146 :
147 : CALL get_qs_env(qs_env, &
148 : natom=natoms, &
149 : cell=cell, &
150 : particle_set=particle_set, &
151 195 : qs_kind_set=qs_kind_set)
152 :
153 195 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
154 195 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, pao_potentials=ipao_potentials)
155 195 : npots = SIZE(ipao_potentials)
156 195 : N = basis_set%nsgf ! primary basis-size
157 195 : CPASSERT(SIZE(V_blocks, 1) == N .AND. SIZE(V_blocks, 2) == N)
158 195 : kterm = 0 ! init counter
159 :
160 391 : DO ipot = 1, npots
161 196 : pot_maxl = ipao_potentials(ipot)%maxl ! taken from central atom
162 :
163 : ! setup description of potential
164 196 : lb_min = 0
165 196 : lb_max = pot_maxl
166 196 : ncfgb = ncoset(lb_max) - ncoset(lb_min - 1)
167 196 : npgfb = 1 ! number of exponents
168 196 : nb = npgfb*ncfgb
169 196 : ALLOCATE (rpgfb(npgfb), zetb(npgfb))
170 :
171 : ! build block_V_full
172 980 : ALLOCATE (block_V_full(N, N, pot_maxl/2 + 1))
173 8116 : block_V_full = 0.0_dp
174 :
175 418 : DO iset = 1, basis_set%nset
176 666 : DO jset = 1, iset
177 :
178 : ! setup iset
179 248 : la1_max = basis_set%lmax(iset)
180 248 : la1_min = basis_set%lmin(iset)
181 248 : npgfa1 = basis_set%npgf(iset)
182 248 : ncfga1 = ncoset(la1_max) - ncoset(la1_min - 1)
183 248 : na1 = npgfa1*ncfga1
184 248 : zeta1 => basis_set%zet(:, iset)
185 248 : rpgfa1 => basis_set%pgf_radius(:, iset)
186 :
187 : ! setup jset
188 248 : la2_max = basis_set%lmax(jset)
189 248 : la2_min = basis_set%lmin(jset)
190 248 : npgfa2 = basis_set%npgf(jset)
191 248 : ncfga2 = ncoset(la2_max) - ncoset(la2_min - 1)
192 248 : na2 = npgfa2*ncfga2
193 248 : zeta2 => basis_set%zet(:, jset)
194 248 : rpgfa2 => basis_set%pgf_radius(:, jset)
195 :
196 : ! radius of most diffuse basis-function
197 3296 : rpgfa_max = MAX(MAXVAL(rpgfa1), MAXVAL(rpgfa2))
198 :
199 : ! allocate space for integrals
200 2232 : ALLOCATE (saab(na1, na2, nb), saal(na1, na2, pot_maxl/2 + 1))
201 150372 : saal = 0.0_dp
202 :
203 : ! loop over neighbors
204 754 : DO jatom = 1, natoms
205 506 : IF (jatom == iatom) CYCLE ! no self-interaction
206 258 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
207 258 : CALL get_qs_kind(qs_kind_set(jkind), pao_potentials=jpao_potentials)
208 258 : IF (SIZE(jpao_potentials) /= npots) THEN
209 0 : CPABORT("Not all KINDs have the same number of PAO_POTENTIAL sections")
210 : END IF
211 :
212 : ! initialize exponents
213 258 : pot_weight = jpao_potentials(ipot)%weight ! taken from remote atom
214 258 : pot_beta = jpao_potentials(ipot)%beta ! taken from remote atom
215 258 : rpgfb(1) = jpao_potentials(ipot)%beta_radius ! taken from remote atom
216 258 : zetb(1) = pot_beta
217 :
218 : ! calculate direction
219 1032 : Ra = particle_set(iatom)%r
220 1032 : Rb = particle_set(jatom)%r
221 258 : Rab = pbc(ra, rb, cell)
222 :
223 : ! distance screening
224 1032 : tab = SQRT(SUM(Rab**2))
225 258 : IF (rpgfa_max + rpgfb(1) < tab) CYCLE
226 :
227 : ! calculate actual integrals
228 689244 : saab = 0.0_dp
229 : CALL overlap_aab(la1_max=la1_max, la1_min=la1_min, npgfa1=npgfa1, rpgfa1=rpgfa1, zeta1=zeta1, &
230 : la2_max=la2_max, la2_min=la2_min, npgfa2=npgfa2, rpgfa2=rpgfa2, zeta2=zeta2, &
231 : lb_max=lb_max, lb_min=lb_min, npgfb=npgfb, rpgfb=rpgfb, zetb=zetb, &
232 258 : rab=Rab, saab=saab)
233 :
234 : ! sum neighbor contributions according to remote atom's weight and normalization
235 1058 : DO lpot = 0, pot_maxl, 2
236 294 : norm2 = (2.0_dp*pot_beta)**(-0.5_dp - lpot)*gamma1(lpot)
237 : ! sum potential terms: POW(x**2 + y**2 + z**2, lpot/2)
238 1436 : DO ic = ncoset(lpot - 1) + 1, ncoset(lpot)
239 2544 : coeff = multinomial(lpot/2, indco(:, ic)/2)
240 959082 : saal(:, :, lpot/2 + 1) = saal(:, :, lpot/2 + 1) + saab(:, :, ic)*coeff*pot_weight/SQRT(norm2)
241 : END DO
242 : END DO
243 : END DO ! jatom
244 :
245 : ! find bounds of set-pair and setup transformation matrices
246 248 : sgfa1 = basis_set%first_sgf(1, iset)
247 248 : sgla1 = sgfa1 + basis_set%nsgf_set(iset) - 1
248 248 : sgfa2 = basis_set%first_sgf(1, jset)
249 248 : sgla2 = sgfa2 + basis_set%nsgf_set(jset) - 1
250 248 : T1 => basis_set%scon(1:na1, sgfa1:sgla1)
251 248 : T2 => basis_set%scon(1:na2, sgfa2:sgla2)
252 :
253 : ! transform into primary basis
254 248 : DO lpot = 0, pot_maxl, 2
255 268 : V12 => block_V_full(sgfa1:sgla1, sgfa2:sgla2, lpot/2 + 1)
256 268 : V21 => block_V_full(sgfa2:sgla2, sgfa1:sgla1, lpot/2 + 1)
257 1735432 : V12 = MATMUL(TRANSPOSE(T1), MATMUL(saal(:, :, lpot/2 + 1), T2))
258 15364 : V21 = TRANSPOSE(V12)
259 : END DO
260 470 : DEALLOCATE (saab, saal)
261 : END DO ! jset
262 : END DO ! iset
263 196 : DEALLOCATE (rpgfb, zetb)
264 :
265 : ! block_V_full is ready -------------------------------------------------------------------
266 : ! split the full blocks into shell-pair sub-blocks
267 418 : DO iset = 1, basis_set%nset
268 666 : DO jset = 1, iset
269 1114 : DO ishell = 1, basis_set%nshell(iset)
270 2792 : DO jshell = 1, basis_set%nshell(jset)
271 1900 : IF (basis_set%l(ishell, iset) == 0 .AND. basis_set%l(jshell, jset) == 0) CYCLE ! covered by central terms
272 1090 : ishell_abs = SUM(basis_set%nshell(1:iset - 1)) + ishell
273 1012 : jshell_abs = SUM(basis_set%nshell(1:jset - 1)) + jshell
274 986 : IF (ishell_abs < jshell_abs) CYCLE
275 :
276 : ! find bounds of shell-pair
277 632 : sgfa1 = basis_set%first_sgf(ishell, iset)
278 632 : sgla1 = basis_set%last_sgf(ishell, iset)
279 632 : sgfa2 = basis_set%first_sgf(jshell, jset)
280 632 : sgla2 = basis_set%last_sgf(jshell, jset)
281 :
282 1276 : DO lpot = 0, pot_maxl, 2
283 728 : kterm = kterm + 1
284 34760 : V_blocks(:, :, kterm) = 0.0_dp
285 10896 : V_blocks(sgfa1:sgla1, sgfa2:sgla2, kterm) = block_V_full(sgfa1:sgla1, sgfa2:sgla2, lpot/2 + 1)
286 14188 : V_blocks(sgfa2:sgla2, sgfa1:sgla1, kterm) = block_V_full(sgfa2:sgla2, sgfa1:sgla1, lpot/2 + 1)
287 : END DO ! lpot
288 : END DO ! jshell
289 : END DO ! ishell
290 : END DO ! jset
291 : END DO ! iset
292 391 : DEALLOCATE (block_V_full)
293 : END DO ! ipot
294 :
295 : ! terms on central atom ----------------------------------------------------------------------
296 :
297 416 : DO iset = 1, basis_set%nset
298 663 : DO jset = 1, iset
299 1109 : DO ishell = 1, basis_set%nshell(iset)
300 2779 : DO jshell = 1, basis_set%nshell(jset)
301 1891 : IF (basis_set%l(ishell, iset) /= basis_set%l(jshell, jset)) CYCLE ! need quadratic block
302 1139 : ishell_abs = SUM(basis_set%nshell(1:iset - 1)) + ishell
303 1139 : jshell_abs = SUM(basis_set%nshell(1:jset - 1)) + jshell
304 1113 : IF (ishell_abs < jshell_abs) CYCLE
305 864 : kterm = kterm + 1
306 864 : sgfa1 = basis_set%first_sgf(ishell, iset)
307 864 : sgla1 = basis_set%last_sgf(ishell, iset)
308 864 : sgfa2 = basis_set%first_sgf(jshell, jset)
309 864 : sgla2 = basis_set%last_sgf(jshell, jset)
310 864 : CPASSERT((sgla1 - sgfa1) == (sgla2 - sgfa2)) ! should be a quadratic block
311 31096 : V_blocks(:, :, kterm) = 0.0_dp
312 2134 : DO i = 1, sgla1 - sgfa1 + 1 ! set diagonal of sub-block
313 1270 : V_blocks(sgfa1 - 1 + i, sgfa2 - 1 + i, kterm) = 1.0_dp
314 2134 : V_blocks(sgfa2 - 1 + i, sgfa1 - 1 + i, kterm) = 1.0_dp
315 : END DO
316 31096 : norm2 = SUM(V_blocks(:, :, kterm)**2)
317 32764 : V_blocks(:, :, kterm) = V_blocks(:, :, kterm)/SQRT(norm2) ! normalize
318 : END DO ! jshell
319 : END DO ! ishell
320 : END DO ! jset
321 : END DO ! iset
322 :
323 195 : CPASSERT(SIZE(V_blocks, 3) == kterm) ! ensure we generated all terms
324 :
325 195 : CALL timestop(handle)
326 195 : END SUBROUTINE linpot_rotinv_calc_terms
327 :
328 : ! **************************************************************************************************
329 : !> \brief Calculate force contribution from rotinv parametrization
330 : !> \param qs_env ...
331 : !> \param iatom ...
332 : !> \param M_blocks ...
333 : !> \param forces ...
334 : ! **************************************************************************************************
335 32 : SUBROUTINE linpot_rotinv_calc_forces(qs_env, iatom, M_blocks, forces)
336 : TYPE(qs_environment_type), POINTER :: qs_env
337 : INTEGER, INTENT(IN) :: iatom
338 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: M_blocks
339 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: forces
340 :
341 : CHARACTER(len=*), PARAMETER :: routineN = 'linpot_rotinv_calc_forces'
342 :
343 : INTEGER :: handle, i, ic, ikind, ipot, iset, ishell, ishell_abs, jatom, jkind, jset, jshell, &
344 : jshell_abs, kterm, la1_max, la1_min, la2_max, la2_min, lb_max, lb_min, lpot, N, na1, na2, &
345 : natoms, nb, ncfga1, ncfga2, ncfgb, npgfa1, npgfa2, npgfb, npots, nshells, pot_maxl, &
346 : sgfa1, sgfa2, sgla1, sgla2
347 : REAL(dp) :: coeff, f, norm2, pot_beta, pot_weight, &
348 : rpgfa_max, tab
349 : REAL(dp), DIMENSION(3) :: Ra, Rab, Rb
350 32 : REAL(dp), DIMENSION(:), POINTER :: rpgfa1, rpgfa2, rpgfb, zeta1, zeta2, zetb
351 32 : REAL(dp), DIMENSION(:, :), POINTER :: block_D, T1, T2
352 32 : REAL(dp), DIMENSION(:, :, :), POINTER :: block_M_full, dab
353 32 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: daab
354 : TYPE(cell_type), POINTER :: cell
355 : TYPE(gto_basis_set_type), POINTER :: basis_set
356 32 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: ipao_potentials, jpao_potentials
357 32 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
358 32 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
359 :
360 32 : CALL timeset(routineN, handle)
361 :
362 : CALL get_qs_env(qs_env, &
363 : natom=natoms, &
364 : cell=cell, &
365 : particle_set=particle_set, &
366 32 : qs_kind_set=qs_kind_set)
367 :
368 32 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
369 32 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, pao_potentials=ipao_potentials)
370 32 : npots = SIZE(ipao_potentials)
371 66 : nshells = SUM(basis_set%nshell)
372 32 : N = basis_set%nsgf ! primary basis-size
373 32 : CPASSERT(SIZE(M_blocks, 1) == N .AND. SIZE(M_blocks, 2) == N)
374 32 : kterm = 0 ! init counter
375 128 : ALLOCATE (block_D(N, N))
376 :
377 64 : DO ipot = 1, npots
378 32 : pot_maxl = ipao_potentials(ipot)%maxl ! taken from central atom
379 :
380 : ! build block_M_full
381 160 : ALLOCATE (block_M_full(N, N, pot_maxl/2 + 1))
382 1048 : block_M_full = 0.0_dp
383 66 : DO iset = 1, basis_set%nset
384 102 : DO jset = 1, iset
385 170 : DO ishell = 1, basis_set%nshell(iset)
386 432 : DO jshell = 1, basis_set%nshell(jset)
387 296 : IF (basis_set%l(ishell, iset) == 0 .AND. basis_set%l(jshell, jset) == 0) CYCLE ! covered by central terms
388 166 : ishell_abs = SUM(basis_set%nshell(1:iset - 1)) + ishell
389 160 : jshell_abs = SUM(basis_set%nshell(1:jset - 1)) + jshell
390 158 : IF (ishell_abs < jshell_abs) CYCLE
391 : ! find bounds of shell-pair
392 98 : sgfa1 = basis_set%first_sgf(ishell, iset)
393 98 : sgla1 = basis_set%last_sgf(ishell, iset)
394 98 : sgfa2 = basis_set%first_sgf(jshell, jset)
395 98 : sgla2 = basis_set%last_sgf(jshell, jset)
396 296 : DO lpot = 0, pot_maxl, 2
397 98 : kterm = kterm + 1
398 746 : block_M_full(sgfa1:sgla1, sgfa2:sgla2, lpot/2 + 1) = M_blocks(sgfa1:sgla1, sgfa2:sgla2, kterm)
399 1174 : block_M_full(sgfa2:sgla2, sgfa1:sgla1, lpot/2 + 1) = M_blocks(sgfa2:sgla2, sgfa1:sgla1, kterm)
400 : END DO ! lpot
401 : END DO ! jshell
402 : END DO ! ishell
403 : END DO ! jset
404 : END DO ! iset
405 :
406 : ! setup description of potential
407 32 : lb_min = 0
408 32 : lb_max = pot_maxl
409 32 : ncfgb = ncoset(lb_max) - ncoset(lb_min - 1)
410 32 : npgfb = 1 ! number of exponents
411 32 : nb = npgfb*ncfgb
412 32 : ALLOCATE (rpgfb(npgfb), zetb(npgfb))
413 :
414 66 : DO iset = 1, basis_set%nset
415 102 : DO jset = 1, iset
416 :
417 : ! setup iset
418 36 : la1_max = basis_set%lmax(iset)
419 36 : la1_min = basis_set%lmin(iset)
420 36 : npgfa1 = basis_set%npgf(iset)
421 36 : ncfga1 = ncoset(la1_max) - ncoset(la1_min - 1)
422 36 : na1 = npgfa1*ncfga1
423 36 : zeta1 => basis_set%zet(:, iset)
424 36 : rpgfa1 => basis_set%pgf_radius(:, iset)
425 :
426 : ! setup jset
427 36 : la2_max = basis_set%lmax(jset)
428 36 : la2_min = basis_set%lmin(jset)
429 36 : npgfa2 = basis_set%npgf(jset)
430 36 : ncfga2 = ncoset(la2_max) - ncoset(la2_min - 1)
431 36 : na2 = npgfa2*ncfga2
432 36 : zeta2 => basis_set%zet(:, jset)
433 36 : rpgfa2 => basis_set%pgf_radius(:, jset)
434 :
435 : ! radius of most diffuse basis-function
436 460 : rpgfa_max = MAX(MAXVAL(rpgfa1), MAXVAL(rpgfa2))
437 :
438 : ! find bounds of set-pair and setup transformation matrices
439 36 : sgfa1 = basis_set%first_sgf(1, iset)
440 36 : sgla1 = sgfa1 + basis_set%nsgf_set(iset) - 1
441 36 : sgfa2 = basis_set%first_sgf(1, jset)
442 36 : sgla2 = sgfa2 + basis_set%nsgf_set(jset) - 1
443 36 : T1 => basis_set%scon(1:na1, sgfa1:sgla1)
444 36 : T2 => basis_set%scon(1:na2, sgfa2:sgla2)
445 :
446 : ! allocate space for integrals
447 360 : ALLOCATE (daab(na1, na2, nb, 3), dab(na1, na2, 3))
448 :
449 : ! loop over neighbors
450 110 : DO jatom = 1, natoms
451 74 : IF (jatom == iatom) CYCLE ! no self-interaction
452 38 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
453 38 : CALL get_qs_kind(qs_kind_set(jkind), pao_potentials=jpao_potentials)
454 38 : IF (SIZE(jpao_potentials) /= npots) THEN
455 0 : CPABORT("Not all KINDs have the same number of PAO_POTENTIAL sections")
456 : END IF
457 :
458 : ! initialize exponents
459 38 : pot_weight = jpao_potentials(ipot)%weight ! taken from remote atom
460 38 : pot_beta = jpao_potentials(ipot)%beta ! taken from remote atom
461 38 : rpgfb(1) = jpao_potentials(ipot)%beta_radius ! taken from remote atom
462 38 : zetb(1) = pot_beta
463 :
464 : ! calculate direction
465 152 : Ra = particle_set(iatom)%r
466 152 : Rb = particle_set(jatom)%r
467 38 : Rab = pbc(ra, rb, cell)
468 :
469 : ! distance screening
470 152 : tab = SQRT(SUM(Rab**2))
471 38 : IF (rpgfa_max + rpgfb(1) < tab) CYCLE
472 :
473 : ! calculate actual integrals
474 59774 : daab = 0.0_dp
475 : CALL overlap_aab(la1_max=la1_max, la1_min=la1_min, npgfa1=npgfa1, rpgfa1=rpgfa1, zeta1=zeta1, &
476 : la2_max=la2_max, la2_min=la2_min, npgfa2=npgfa2, rpgfa2=rpgfa2, zeta2=zeta2, &
477 : lb_max=lb_max, lb_min=lb_min, npgfb=npgfb, rpgfb=rpgfb, zetb=zetb, &
478 38 : rab=Rab, daab=daab)
479 :
480 : ! sum neighbor contributions according to remote atom's weight and normalization
481 150 : DO lpot = 0, pot_maxl, 2
482 : ! sum potential terms: POW(x**2 + y**2 + z**2, lpot/2)
483 59660 : dab = 0.0_dp
484 76 : DO ic = ncoset(lpot - 1) + 1, ncoset(lpot)
485 38 : norm2 = (2.0_dp*pot_beta)**(-0.5_dp - lpot)*gamma1(lpot)
486 152 : coeff = multinomial(lpot/2, indco(:, ic)/2)
487 119320 : dab = dab + coeff*daab(:, :, ic, :)*pot_weight/SQRT(norm2)
488 : END DO
489 226 : DO i = 1, 3
490 : ! transform into primary basis
491 3750 : block_D = 0.0_dp
492 728652 : block_D(sgfa1:sgla1, sgfa2:sgla2) = MATMUL(TRANSPOSE(T1), MATMUL(dab(:, :, i), T2))
493 6306 : block_D(sgfa2:sgla2, sgfa1:sgla1) = TRANSPOSE(block_D(sgfa1:sgla1, sgfa2:sgla2))
494 : ! calculate and add forces
495 3750 : f = SUM(block_M_full(:, :, lpot/2 + 1)*block_D)
496 114 : forces(iatom, i) = forces(iatom, i) - f
497 152 : forces(jatom, i) = forces(jatom, i) + f
498 : END DO
499 : END DO ! lpot
500 : END DO ! jatom
501 70 : DEALLOCATE (dab, daab)
502 : END DO ! jset
503 : END DO ! iset
504 64 : DEALLOCATE (rpgfb, zetb, block_M_full)
505 : END DO ! ipot
506 32 : DEALLOCATE (block_D)
507 :
508 32 : CALL timestop(handle)
509 64 : END SUBROUTINE linpot_rotinv_calc_forces
510 :
511 382 : END MODULE pao_linpot_rotinv
|
