Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Factory routines for potentials used e.g. by pao_param_exp and pao_ml
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE pao_potentials
13 : USE ai_overlap, ONLY: overlap_aab
14 : USE ao_util, ONLY: exp_radius
15 : USE atomic_kind_types, ONLY: get_atomic_kind
16 : USE basis_set_types, ONLY: gto_basis_set_type
17 : USE cell_types, ONLY: cell_type,&
18 : pbc
19 : USE kinds, ONLY: dp
20 : USE mathconstants, ONLY: gamma1
21 : USE mathlib, ONLY: multinomial
22 : USE orbital_pointers, ONLY: indco,&
23 : ncoset,&
24 : orbital_pointers_maxl => current_maxl
25 : USE particle_types, ONLY: particle_type
26 : USE qs_environment_types, ONLY: get_qs_env,&
27 : qs_environment_type
28 : USE qs_kind_types, ONLY: get_qs_kind,&
29 : qs_kind_type
30 : #include "./base/base_uses.f90"
31 :
32 : IMPLICIT NONE
33 :
34 : PRIVATE
35 :
36 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_potentials'
37 :
38 : PUBLIC :: pao_guess_initial_potential, pao_calc_gaussian
39 :
40 : CONTAINS
41 :
42 : ! **************************************************************************************************
43 : !> \brief Makes an educated guess for the initial potential based on positions of neighboring atoms
44 : !> \param qs_env ...
45 : !> \param iatom ...
46 : !> \param block_V ...
47 : ! **************************************************************************************************
48 76 : SUBROUTINE pao_guess_initial_potential(qs_env, iatom, block_V)
49 : TYPE(qs_environment_type), POINTER :: qs_env
50 : INTEGER, INTENT(IN) :: iatom
51 : REAL(dp), DIMENSION(:, :), INTENT(OUT) :: block_V
52 :
53 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_guess_initial_potential'
54 :
55 : INTEGER :: handle, ikind, jatom, natoms
56 : REAL(dp), DIMENSION(3) :: Ra, Rab, Rb
57 : TYPE(cell_type), POINTER :: cell
58 : TYPE(gto_basis_set_type), POINTER :: basis_set
59 76 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
60 76 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
61 :
62 76 : CALL timeset(routineN, handle)
63 :
64 : CALL get_qs_env(qs_env, &
65 : cell=cell, &
66 : particle_set=particle_set, &
67 : qs_kind_set=qs_kind_set, &
68 76 : natom=natoms)
69 :
70 76 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
71 76 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set)
72 :
73 : ! construct matrix block_V from neighboring atoms
74 4508 : block_V = 0.0_dp
75 270 : DO jatom = 1, natoms
76 194 : IF (jatom == iatom) CYCLE
77 472 : Ra = particle_set(iatom)%r
78 472 : Rb = particle_set(jatom)%r
79 118 : Rab = pbc(ra, rb, cell)
80 270 : CALL pao_calc_gaussian(basis_set, block_V, Rab=Rab, lpot=0, beta=1.0_dp, weight=-1.0_dp)
81 : END DO
82 :
83 76 : CALL timestop(handle)
84 76 : END SUBROUTINE pao_guess_initial_potential
85 :
86 : ! **************************************************************************************************
87 : !> \brief Calculates potential term of the form r**lpot * Exp(-beta*r**2)
88 : !> One needs to call init_orbital_pointers(lpot) before calling pao_calc_gaussian().
89 : !> \param basis_set ...
90 : !> \param block_V potential term that is returned
91 : !> \param block_D derivative of potential term wrt to Rab
92 : !> \param Rab ...
93 : !> \param lpot polynomial prefactor, r**lpot
94 : !> \param beta exponent of the Gaussian
95 : !> \param weight ...
96 : !> \param min_shell ...
97 : !> \param max_shell ...
98 : !> \param min_l ...
99 : !> \param max_l ...
100 : ! **************************************************************************************************
101 478 : SUBROUTINE pao_calc_gaussian(basis_set, block_V, block_D, Rab, lpot, beta, weight, min_shell, max_shell, min_l, max_l)
102 : TYPE(gto_basis_set_type), POINTER :: basis_set
103 : REAL(dp), DIMENSION(:, :), INTENT(OUT), OPTIONAL :: block_V
104 : REAL(dp), DIMENSION(:, :, :), INTENT(OUT), &
105 : OPTIONAL :: block_D
106 : REAL(dp), DIMENSION(3) :: Rab
107 : INTEGER, INTENT(IN) :: lpot
108 : REAL(dp), INTENT(IN) :: beta, weight
109 : INTEGER, INTENT(IN), OPTIONAL :: min_shell, max_shell, min_l, max_l
110 :
111 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_gaussian'
112 :
113 : INTEGER :: handle, i, ic, iset, ishell, ishell_abs, jset, jshell, jshell_abs, la1_max, &
114 : la1_min, la2_max, la2_min, lb_max, lb_min, N, na1, na2, nb, ncfga1, ncfga2, ncfgb, &
115 : npgfa1, npgfa2, npgfb
116 : REAL(dp) :: coeff, norm2
117 478 : REAL(dp), DIMENSION(:), POINTER :: rpgfa1, rpgfa2, rpgfb, zeta1, zeta2, zetb
118 478 : REAL(dp), DIMENSION(:, :), POINTER :: new_block_V, sab
119 478 : REAL(dp), DIMENSION(:, :, :), POINTER :: dab, new_block_D, saab
120 478 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: daab
121 :
122 478 : CALL timeset(routineN, handle)
123 :
124 478 : CPASSERT(PRESENT(block_V) .NEQV. PRESENT(block_D)) ! just to keep the code simpler
125 478 : CPASSERT(PRESENT(min_shell) .EQV. PRESENT(max_shell))
126 478 : CPASSERT(PRESENT(min_l) .EQV. PRESENT(max_l))
127 :
128 478 : CPASSERT(MOD(lpot, 2) == 0) ! otherwise it's not rotationally invariant
129 478 : CPASSERT(orbital_pointers_maxl >= lpot) ! can't call init_orbital_pointers here, it's not thread-safe
130 :
131 478 : N = basis_set%nsgf ! primary basis-size
132 :
133 478 : IF (PRESENT(block_V)) THEN
134 468 : CPASSERT(SIZE(block_V, 1) == N .AND. SIZE(block_V, 2) == N)
135 1872 : ALLOCATE (new_block_V(N, N))
136 26084 : new_block_V = 0.0_dp
137 : END IF
138 478 : IF (PRESENT(block_D)) THEN
139 10 : CPASSERT(SIZE(block_D, 1) == N .AND. SIZE(block_D, 2) == N .AND. SIZE(block_D, 3) == 3)
140 50 : ALLOCATE (new_block_D(N, N, 3))
141 940 : new_block_D = 0.0_dp
142 : END IF
143 :
144 : ! setup description of potential
145 478 : lb_min = lpot
146 478 : lb_max = lpot
147 478 : ncfgb = ncoset(lb_max) - ncoset(lb_min - 1)
148 478 : npgfb = 1 ! number of exponents
149 478 : nb = npgfb*ncfgb
150 :
151 : ! initialize exponents
152 478 : ALLOCATE (rpgfb(npgfb), zetb(npgfb))
153 478 : rpgfb(1) = exp_radius(0, beta, 1.0E-12_dp, 1.0_dp) ! TODO get the EPS parameter from somewhere / precompute this elsewhere
154 478 : zetb(1) = beta
155 :
156 : ! loop over all set/shell combination and fill block_V
157 958 : DO iset = 1, basis_set%nset
158 1442 : DO jset = 1, basis_set%nset
159 2560 : DO ishell = 1, basis_set%nshell(iset)
160 7612 : DO jshell = 1, basis_set%nshell(jset)
161 5532 : IF (PRESENT(min_shell) .AND. PRESENT(max_shell)) THEN
162 0 : ishell_abs = SUM(basis_set%nshell(1:iset - 1)) + ishell
163 0 : jshell_abs = SUM(basis_set%nshell(1:jset - 1)) + jshell
164 0 : IF (MIN(ishell_abs, jshell_abs) /= min_shell) CYCLE
165 0 : IF (MAX(ishell_abs, jshell_abs) /= max_shell) CYCLE
166 : END IF
167 5532 : IF (PRESENT(min_l) .AND. PRESENT(min_l)) THEN
168 3036 : IF (MIN(basis_set%l(ishell, iset), basis_set%l(jshell, jset)) /= min_l) CYCLE
169 1436 : IF (MAX(basis_set%l(ishell, iset), basis_set%l(jshell, jset)) /= max_l) CYCLE
170 : END IF
171 :
172 : ! setup iset
173 3038 : la1_max = basis_set%l(ishell, iset)
174 3038 : la1_min = basis_set%l(ishell, iset)
175 3038 : npgfa1 = basis_set%npgf(iset)
176 3038 : ncfga1 = ncoset(la1_max) - ncoset(la1_min - 1)
177 3038 : na1 = npgfa1*ncfga1
178 3038 : zeta1 => basis_set%zet(:, iset)
179 3038 : rpgfa1 => basis_set%pgf_radius(:, iset)
180 :
181 : ! setup jset
182 3038 : la2_max = basis_set%l(jshell, jset)
183 3038 : la2_min = basis_set%l(jshell, jset)
184 3038 : npgfa2 = basis_set%npgf(jset)
185 3038 : ncfga2 = ncoset(la2_max) - ncoset(la2_min - 1)
186 3038 : na2 = npgfa2*ncfga2
187 3038 : zeta2 => basis_set%zet(:, jset)
188 3038 : rpgfa2 => basis_set%pgf_radius(:, jset)
189 :
190 : ! calculate integrals
191 3038 : IF (PRESENT(block_V)) THEN
192 14740 : ALLOCATE (saab(na1, na2, nb))
193 899278 : saab = 0.0_dp
194 : CALL overlap_aab(la1_max=la1_max, la1_min=la1_min, npgfa1=npgfa1, rpgfa1=rpgfa1, zeta1=zeta1, &
195 : la2_max=la2_max, la2_min=la2_min, npgfa2=npgfa2, rpgfa2=rpgfa2, zeta2=zeta2, &
196 : lb_max=lb_max, lb_min=lb_min, npgfb=npgfb, rpgfb=rpgfb, zetb=zetb, &
197 2948 : rab=Rab, saab=saab)
198 : END IF
199 :
200 3038 : IF (PRESENT(block_D)) THEN
201 540 : ALLOCATE (daab(na1, na2, nb, 3))
202 40530 : daab = 0.0_dp
203 : CALL overlap_aab(la1_max=la1_max, la1_min=la1_min, npgfa1=npgfa1, rpgfa1=rpgfa1, zeta1=zeta1, &
204 : la2_max=la2_max, la2_min=la2_min, npgfa2=npgfa2, rpgfa2=rpgfa2, zeta2=zeta2, &
205 : lb_max=lb_max, lb_min=lb_min, npgfb=npgfb, rpgfb=rpgfb, zetb=zetb, &
206 90 : rab=Rab, daab=daab)
207 : END IF
208 :
209 : ! sum potential terms: POW(x**2 + y**2 + z**2, lpot/2)
210 3038 : IF (PRESENT(block_V)) THEN
211 11792 : ALLOCATE (sab(na1, na2))
212 415060 : sab = 0.0_dp
213 11166 : DO ic = 1, ncfgb
214 32872 : coeff = multinomial(lpot/2, indco(:, ncoset(lpot - 1) + ic)/2)
215 1787390 : sab = sab + coeff*saab(:, :, ic)
216 : END DO
217 : CALL my_contract(sab=sab, block=new_block_V, basis_set=basis_set, &
218 2948 : iset=iset, ishell=ishell, jset=jset, jshell=jshell)
219 2948 : DEALLOCATE (sab, saab)
220 : END IF
221 :
222 4634 : IF (PRESENT(block_D)) THEN
223 450 : ALLOCATE (dab(na1, na2, 3))
224 40260 : dab = 0.0_dp
225 180 : DO ic = 1, ncfgb
226 360 : coeff = multinomial(lpot/2, indco(:, ncoset(lpot - 1) + ic)/2)
227 80520 : dab = dab + coeff*daab(:, :, ic, :)
228 : END DO
229 360 : DO i = 1, 3
230 : CALL my_contract(sab=dab(:, :, i), block=new_block_D(:, :, i), basis_set=basis_set, &
231 360 : iset=iset, ishell=ishell, jset=jset, jshell=jshell)
232 : END DO
233 90 : DEALLOCATE (dab, daab)
234 : END IF
235 : END DO
236 : END DO
237 : END DO
238 : END DO
239 :
240 478 : DEALLOCATE (rpgfb, zetb)
241 :
242 : ! post-processing
243 478 : norm2 = (2.0_dp*beta)**(-0.5_dp - lpot)*gamma1(lpot)
244 478 : IF (PRESENT(block_V)) THEN
245 26084 : block_V = block_V + weight*new_block_V/SQRT(norm2)
246 468 : DEALLOCATE (new_block_V)
247 51700 : block_V = 0.5_dp*(block_V + TRANSPOSE(block_V)) ! symmetrize
248 : END IF
249 :
250 478 : IF (PRESENT(block_D)) THEN
251 940 : block_D = block_D + weight*new_block_D/SQRT(norm2)
252 10 : DEALLOCATE (new_block_D)
253 40 : DO i = 1, 3
254 1840 : block_D(:, :, i) = 0.5_dp*(block_D(:, :, i) + TRANSPOSE(block_D(:, :, i))) ! symmetrize
255 : END DO
256 : END IF
257 :
258 478 : CALL timestop(handle)
259 478 : END SUBROUTINE pao_calc_gaussian
260 :
261 : ! **************************************************************************************************
262 : !> \brief Helper routine, contracts a basis block
263 : !> \param sab ...
264 : !> \param block ...
265 : !> \param basis_set ...
266 : !> \param iset ...
267 : !> \param ishell ...
268 : !> \param jset ...
269 : !> \param jshell ...
270 : ! **************************************************************************************************
271 3218 : SUBROUTINE my_contract(sab, block, basis_set, iset, ishell, jset, jshell)
272 : REAL(dp), DIMENSION(:, :), INTENT(IN), TARGET :: sab
273 : REAL(dp), DIMENSION(:, :), INTENT(OUT), TARGET :: block
274 : TYPE(gto_basis_set_type), POINTER :: basis_set
275 : INTEGER, INTENT(IN) :: iset, ishell, jset, jshell
276 :
277 : INTEGER :: a, b, c, d, ipgf, jpgf, l1, l2, n1, n2, &
278 : nn1, nn2, sgfa1, sgfa2, sgla1, sgla2
279 3218 : REAL(dp), DIMENSION(:, :), POINTER :: S, T1, T2, V
280 :
281 : ! first and last indices of given shell in block.
282 : ! This matrix is in the contracted spherical basis.
283 3218 : sgfa1 = basis_set%first_sgf(ishell, iset)
284 3218 : sgla1 = basis_set%last_sgf(ishell, iset)
285 3218 : sgfa2 = basis_set%first_sgf(jshell, jset)
286 3218 : sgla2 = basis_set%last_sgf(jshell, jset)
287 :
288 : ! prepare the result matrix
289 3218 : V => block(sgfa1:sgla1, sgfa2:sgla2)
290 :
291 : ! Calculate strides of sphi matrix.
292 : ! This matrix is in the uncontraced cartesian basis.
293 : ! It contains all shells of the set.
294 : ! Its index runs over all primitive gaussians of the set
295 : ! and then for each gaussian over all configurations of *the entire set*. (0->lmax)
296 3218 : nn1 = ncoset(basis_set%lmax(iset))
297 3218 : nn2 = ncoset(basis_set%lmax(jset))
298 :
299 : ! Calculate strides of sab matrix
300 : ! This matrix is also in the uncontraced cartensian basis,
301 : ! however it contains only a single shell.
302 : ! Its index runs over all primitive gaussians of the set
303 : ! and then for each gaussian over all configrations of *the given shell*.
304 3218 : l1 = basis_set%l(ishell, iset)
305 3218 : l2 = basis_set%l(jshell, jset)
306 3218 : n1 = ncoset(l1) - ncoset(l1 - 1)
307 3218 : n2 = ncoset(l2) - ncoset(l2 - 1)
308 :
309 21688 : DO ipgf = 1, basis_set%npgf(iset)
310 130794 : DO jpgf = 1, basis_set%npgf(jset)
311 : ! prepare first trafo-matrix
312 109106 : a = (ipgf - 1)*nn1 + ncoset(l1 - 1) + 1
313 109106 : T1 => basis_set%sphi(a:a + n1 - 1, sgfa1:sgla1)
314 :
315 : ! prepare second trafo-matrix
316 109106 : b = (jpgf - 1)*nn2 + ncoset(l2 - 1) + 1
317 109106 : T2 => basis_set%sphi(b:b + n2 - 1, sgfa2:sgla2)
318 :
319 : ! prepare SAB matrix
320 109106 : c = (ipgf - 1)*n1 + 1
321 109106 : d = (jpgf - 1)*n2 + 1
322 109106 : S => sab(c:c + n1 - 1, d:d + n2 - 1)
323 :
324 : ! do the transformation
325 10895956 : V = V + MATMUL(TRANSPOSE(T1), MATMUL(S, T2))
326 : END DO
327 : END DO
328 :
329 3218 : END SUBROUTINE my_contract
330 :
331 109106 : END MODULE pao_potentials
|
