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