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 : MODULE atom_set_basis
10 : USE ai_onecenter, ONLY: sg_overlap
11 : USE atom_types, ONLY: CGTO_BASIS,&
12 : GTO_BASIS,&
13 : atom_basis_type,&
14 : lmat
15 : USE basis_set_types, ONLY: get_gto_basis_set,&
16 : gto_basis_set_type
17 : USE input_constants, ONLY: do_gapw_log
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: dfac,&
20 : twopi
21 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
22 : create_grid_atom,&
23 : grid_atom_type
24 : #include "./base/base_uses.f90"
25 :
26 : IMPLICIT NONE
27 :
28 : PRIVATE
29 :
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_set_basis'
31 :
32 : INTEGER, PARAMETER :: nua = 40, nup = 20
33 : REAL(KIND=dp), DIMENSION(nua), PARAMETER :: ugbs = [0.007299_dp, 0.013705_dp, 0.025733_dp, &
34 : 0.048316_dp, 0.090718_dp, 0.170333_dp, 0.319819_dp, 0.600496_dp, 1.127497_dp, 2.117000_dp, &
35 : 3.974902_dp, 7.463317_dp, 14.013204_dp, 26.311339_dp, 49.402449_dp, 92.758561_dp, &
36 : 174.164456_dp, 327.013024_dp, 614.003114_dp, 1152.858743_dp, 2164.619772_dp, &
37 : 4064.312984_dp, 7631.197056_dp, 14328.416324_dp, 26903.186074_dp, 50513.706789_dp, &
38 : 94845.070265_dp, 178082.107320_dp, 334368.848683_dp, 627814.487663_dp, 1178791.123851_dp, &
39 : 2213310.684886_dp, 4155735.557141_dp, 7802853.046713_dp, 14650719.428954_dp, &
40 : 27508345.793637_dp, 51649961.080194_dp, 96978513.342764_dp, 182087882.613702_dp, &
41 : 341890134.751331_dp]
42 :
43 : PUBLIC :: set_kind_basis_atomic
44 :
45 : ! **************************************************************************************************
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief ...
51 : !> \param basis ...
52 : !> \param orb_basis_set ...
53 : !> \param has_pp ...
54 : !> \param agrid ...
55 : !> \param cp2k_norm ...
56 : ! **************************************************************************************************
57 8640 : SUBROUTINE set_kind_basis_atomic(basis, orb_basis_set, has_pp, agrid, cp2k_norm)
58 : TYPE(atom_basis_type), INTENT(INOUT) :: basis
59 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
60 : LOGICAL, INTENT(IN) :: has_pp
61 : TYPE(grid_atom_type), OPTIONAL :: agrid
62 : LOGICAL, INTENT(IN), OPTIONAL :: cp2k_norm
63 :
64 : INTEGER :: i, ii, ipgf, j, k, l, m, ngp, nj, nr, &
65 : ns, nset, nsgf, nu, quadtype
66 8640 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
67 8640 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, last_sgf, ls
68 : LOGICAL :: has_basis, set_norm
69 : REAL(KIND=dp) :: al, an, cn, ear, en, rk
70 8640 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
71 8640 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
72 : TYPE(grid_atom_type), POINTER :: grid
73 :
74 8640 : IF (ASSOCIATED(orb_basis_set)) THEN
75 : has_basis = .TRUE.
76 : ELSE
77 6 : has_basis = .FALSE.
78 : END IF
79 :
80 8640 : IF (PRESENT(cp2k_norm)) THEN
81 0 : set_norm = cp2k_norm
82 : ELSE
83 : set_norm = .FALSE.
84 : END IF
85 :
86 8640 : NULLIFY (grid)
87 8640 : IF (PRESENT(agrid)) THEN
88 6 : ngp = agrid%nr
89 6 : quadtype = agrid%quadrature
90 : ELSE
91 8634 : ngp = 400
92 8634 : quadtype = do_gapw_log
93 : END IF
94 8640 : CALL allocate_grid_atom(grid)
95 8640 : CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
96 8640 : grid%nr = ngp
97 8640 : basis%grid => grid
98 :
99 8640 : NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
100 :
101 8640 : IF (has_basis) THEN
102 : ! fill in the basis data structures
103 8634 : basis%basis_type = CGTO_BASIS
104 8634 : basis%eps_eig = 1.e-12_dp
105 : CALL get_gto_basis_set(orb_basis_set, &
106 : nset=nset, nshell=nshell, npgf=npgf, lmin=lmin, lmax=lmax, &
107 : l=ls, nsgf=nsgf, zet=zet, gcc=gcc, &
108 8634 : first_sgf=first_sgf, last_sgf=last_sgf)
109 60438 : basis%nprim = 0
110 60438 : basis%nbas = 0
111 26252 : DO i = 1, nset
112 41829 : DO j = lmin(i), MIN(lmax(i), lmat)
113 41829 : basis%nprim(j) = basis%nprim(j) + npgf(i)
114 : END DO
115 60434 : DO j = 1, nshell(i)
116 34182 : l = ls(j, i)
117 51800 : IF (l <= lmat) THEN
118 34182 : basis%nbas(l) = basis%nbas(l) + 1
119 34182 : k = basis%nbas(l)
120 : END IF
121 : END DO
122 : END DO
123 :
124 60438 : nj = MAXVAL(basis%nprim)
125 60438 : ns = MAXVAL(basis%nbas)
126 25896 : ALLOCATE (basis%am(nj, 0:lmat))
127 307122 : basis%am = 0._dp
128 43158 : ALLOCATE (basis%cm(nj, ns, 0:lmat))
129 702342 : basis%cm = 0._dp
130 60438 : DO j = 0, lmat
131 51804 : nj = 0
132 51804 : ns = 0
133 51804 : cn = 2.0_dp**(j + 2)/SQRT(dfac(2*j + 1))/twopi**0.25_dp
134 51804 : en = (2*j + 3)*0.25_dp
135 166146 : DO i = 1, nset
136 157512 : IF (j >= lmin(i) .AND. j <= lmax(i)) THEN
137 99413 : DO ipgf = 1, npgf(i)
138 99413 : basis%am(nj + ipgf, j) = zet(ipgf, i)
139 : END DO
140 82447 : DO ii = 1, nshell(i)
141 82447 : IF (ls(ii, i) == j) THEN
142 34182 : ns = ns + 1
143 34182 : IF (set_norm) THEN
144 0 : DO ipgf = 1, npgf(i)
145 0 : an = cn*zet(ipgf, i)**en
146 0 : basis%cm(nj + ipgf, ns, j) = an*gcc(ipgf, ii, i)
147 : END DO
148 : ELSE
149 157522 : DO ipgf = 1, npgf(i)
150 157522 : basis%cm(nj + ipgf, ns, j) = gcc(ipgf, ii, i)
151 : END DO
152 : END IF
153 : END IF
154 : END DO
155 24211 : nj = nj + npgf(i)
156 : END IF
157 : END DO
158 : END DO
159 : ! Normalization
160 17268 : IF (set_norm) THEN
161 0 : CALL normalize_basis_cp2k(basis)
162 : END IF
163 : ELSE
164 : ! use default basis
165 6 : IF (has_pp) THEN
166 : nu = nup
167 : ELSE
168 0 : nu = nua
169 : END IF
170 6 : basis%geometrical = .FALSE.
171 6 : basis%aval = 0._dp
172 6 : basis%cval = 0._dp
173 42 : basis%start = 0
174 6 : basis%eps_eig = 1.e-12_dp
175 :
176 6 : basis%basis_type = GTO_BASIS
177 42 : basis%nbas = nu
178 42 : basis%nprim = nu
179 12 : ALLOCATE (basis%am(nu, 0:lmat))
180 42 : DO i = 0, lmat
181 762 : basis%am(1:nu, i) = ugbs(1:nu)
182 : END DO
183 : END IF
184 :
185 : ! initialize basis function on a radial grid
186 8640 : nr = basis%grid%nr
187 60480 : m = MAXVAL(basis%nbas)
188 43194 : ALLOCATE (basis%bf(nr, m, 0:lmat))
189 25914 : ALLOCATE (basis%dbf(nr, m, 0:lmat))
190 25914 : ALLOCATE (basis%ddbf(nr, m, 0:lmat))
191 :
192 42806208 : basis%bf = 0._dp
193 42806208 : basis%dbf = 0._dp
194 42806208 : basis%ddbf = 0._dp
195 60480 : DO l = 0, lmat
196 136402 : DO i = 1, basis%nprim(l)
197 75922 : al = basis%am(i, l)
198 127762 : IF (basis%basis_type == GTO_BASIS) THEN
199 669600 : DO k = 1, nr
200 668880 : rk = basis%grid%rad(k)
201 668880 : ear = EXP(-al*basis%grid%rad(k)**2)
202 668880 : basis%bf(k, i, l) = rk**l*ear
203 668880 : basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
204 : basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - &
205 669600 : 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
206 : END DO
207 75202 : ELSEIF (basis%basis_type == CGTO_BASIS) THEN
208 30156002 : DO k = 1, nr
209 30080800 : rk = basis%grid%rad(k)
210 30080800 : ear = EXP(-al*basis%grid%rad(k)**2)
211 90036802 : DO j = 1, basis%nbas(l)
212 59880800 : basis%bf(k, j, l) = basis%bf(k, j, l) + rk**l*ear*basis%cm(i, j, l)
213 : basis%dbf(k, j, l) = basis%dbf(k, j, l) &
214 59880800 : + (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis%cm(i, j, l)
215 : basis%ddbf(k, j, l) = basis%ddbf(k, j, l) + &
216 : (REAL(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + &
217 89961600 : 4._dp*al*rk**(l + 2))*ear*basis%cm(i, j, l)
218 : END DO
219 : END DO
220 : ELSE
221 0 : CPABORT('Atom basis type?')
222 : END IF
223 : END DO
224 : END DO
225 :
226 8640 : END SUBROUTINE set_kind_basis_atomic
227 :
228 : ! **************************************************************************************************
229 : !> \brief ...
230 : !> \param basis ...
231 : ! **************************************************************************************************
232 0 : SUBROUTINE normalize_basis_cp2k(basis)
233 : TYPE(atom_basis_type), INTENT(INOUT) :: basis
234 :
235 : INTEGER :: ii, l, n, np
236 : REAL(KIND=dp) :: fnorm
237 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: smat
238 :
239 0 : DO l = 0, lmat
240 0 : n = basis%nbas(l)
241 0 : np = basis%nprim(l)
242 0 : IF (n > 0) THEN
243 0 : ALLOCATE (smat(np, np))
244 0 : CALL sg_overlap(smat, l, basis%am(1:np, l), basis%am(1:np, l))
245 0 : DO ii = 1, basis%nbas(l)
246 0 : fnorm = DOT_PRODUCT(basis%cm(1:np, ii, l), MATMUL(smat, basis%cm(1:np, ii, l)))
247 0 : fnorm = 1._dp/SQRT(fnorm)
248 0 : basis%cm(1:np, ii, l) = fnorm*basis%cm(1:np, ii, l)
249 : END DO
250 0 : DEALLOCATE (smat)
251 : END IF
252 : END DO
253 :
254 0 : END SUBROUTINE normalize_basis_cp2k
255 :
256 : ! **************************************************************************************************
257 :
258 0 : END MODULE atom_set_basis
|
