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 : !> \par History
10 : !> none
11 : !> \author MI (08.01.2004)
12 : ! **************************************************************************************************
13 : MODULE soft_basis_set
14 :
15 : USE ao_util, ONLY: exp_radius
16 : USE basis_set_types, ONLY: copy_gto_basis_set,&
17 : deallocate_gto_basis_set,&
18 : get_gto_basis_set,&
19 : gto_basis_set_type,&
20 : init_cphi_and_sphi
21 : USE kinds, ONLY: default_string_length,&
22 : dp
23 : USE memory_utilities, ONLY: reallocate
24 : USE orbital_pointers, ONLY: indco,&
25 : nco,&
26 : ncoset,&
27 : nso
28 : USE orbital_symbols, ONLY: cgf_symbol,&
29 : sgf_symbol
30 : #include "../base/base_uses.f90"
31 :
32 : IMPLICIT NONE
33 :
34 : PRIVATE
35 :
36 : ! *** Global parameters (only in this module)
37 :
38 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'soft_basis_set'
39 :
40 : INTEGER, PARAMETER :: max_name_length = 60
41 :
42 : ! *** Public subroutines ***
43 :
44 : PUBLIC :: create_soft_basis
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief create the soft basis from a GTO basis
50 : !> \param orb_basis ...
51 : !> \param soft_basis ...
52 : !> \param eps_fit ...
53 : !> \param rc ...
54 : !> \param paw_atom ...
55 : !> \param paw_type_forced ...
56 : !> \param gpw_r3d_rs_type_forced ...
57 : !> \version 1.0
58 : ! **************************************************************************************************
59 1932 : SUBROUTINE create_soft_basis(orb_basis, soft_basis, eps_fit, rc, paw_atom, &
60 : paw_type_forced, gpw_r3d_rs_type_forced)
61 :
62 : TYPE(gto_basis_set_type), POINTER :: orb_basis, soft_basis
63 : REAL(dp), INTENT(IN) :: eps_fit, rc
64 : LOGICAL, INTENT(OUT) :: paw_atom
65 : LOGICAL, INTENT(IN) :: paw_type_forced, gpw_r3d_rs_type_forced
66 :
67 : CHARACTER(LEN=default_string_length) :: bsname
68 : INTEGER :: ico, ipgf, ipgf_s, iset, iset_s, ishell, lshell, lshell_old, m, maxco, maxpgf, &
69 : maxpgf_s, maxshell, maxshell_s, ncgf, nset, nset_s, nsgf
70 1932 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iset_s2h
71 1932 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
72 1932 : INTEGER, DIMENSION(:, :), POINTER :: l, n
73 : LOGICAL :: my_gpw_r3d_rs_type_forced
74 : REAL(KIND=dp) :: minzet, radius
75 1932 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
76 1932 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
77 :
78 1932 : NULLIFY (gcc, l, lmax, lmin, n, npgf, nshell, zet)
79 1932 : paw_atom = .FALSE.
80 1932 : my_gpw_r3d_rs_type_forced = gpw_r3d_rs_type_forced
81 1932 : IF (paw_type_forced) THEN
82 18 : paw_atom = .TRUE.
83 : my_gpw_r3d_rs_type_forced = .FALSE.
84 : END IF
85 1914 : IF (.NOT. my_gpw_r3d_rs_type_forced) THEN
86 : CALL get_gto_basis_set(gto_basis_set=orb_basis, name=bsname, &
87 1904 : maxpgf=maxpgf, maxshell=maxshell, nset=nset)
88 :
89 1904 : soft_basis%name = TRIM(bsname)//"_soft"
90 :
91 1904 : CALL reallocate(npgf, 1, nset)
92 1904 : CALL reallocate(nshell, 1, nset)
93 1904 : CALL reallocate(lmax, 1, nset)
94 1904 : CALL reallocate(lmin, 1, nset)
95 :
96 1904 : CALL reallocate(n, 1, maxshell, 1, nset)
97 1904 : CALL reallocate(l, 1, maxshell, 1, nset)
98 :
99 1904 : CALL reallocate(zet, 1, maxpgf, 1, nset)
100 1904 : CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
101 :
102 5712 : ALLOCATE (iset_s2h(nset))
103 :
104 7416 : iset_s2h = 0
105 1904 : iset_s = 0
106 1904 : maxpgf_s = 0
107 1904 : maxshell_s = 0
108 :
109 7416 : DO iset = 1, nset ! iset
110 5512 : minzet = orb_basis%zet(orb_basis%npgf(iset), iset)
111 13338 : DO ipgf = orb_basis%npgf(iset) - 1, 1, -1
112 13338 : IF (orb_basis%zet(ipgf, iset) < minzet) THEN
113 42 : minzet = orb_basis%zet(ipgf, iset)
114 : END IF
115 : END DO
116 5512 : radius = exp_radius(orb_basis%lmax(iset), minzet, eps_fit, 1.0_dp)
117 :
118 : ! The soft basis contains this set
119 5512 : iset_s = iset_s + 1
120 5512 : nshell(iset_s) = orb_basis%nshell(iset)
121 5512 : lmax(iset_s) = orb_basis%lmax(iset)
122 5512 : lmin(iset_s) = orb_basis%lmin(iset)
123 :
124 5512 : iset_s2h(iset_s) = iset
125 :
126 14408 : DO ishell = 1, nshell(iset_s)
127 8896 : n(ishell, iset_s) = orb_basis%n(ishell, iset)
128 14408 : l(ishell, iset_s) = orb_basis%l(ishell, iset)
129 : END DO
130 :
131 5512 : IF (nshell(iset_s) > maxshell_s) THEN
132 2342 : maxshell_s = nshell(iset_s)
133 : END IF
134 :
135 5512 : IF (radius < rc) THEN
136 : ! The soft basis does not contain this set
137 : ! For the moment I keep the set as a dummy set
138 : ! with no exponents, in order to have the right number of contractions
139 : ! In a second time it can be taken away, by creating a pointer
140 : ! which connects the remaining sets to the right contraction index
141 486 : paw_atom = .TRUE.
142 486 : npgf(iset_s) = 0
143 486 : CYCLE
144 : END IF
145 :
146 5026 : ipgf_s = 0
147 16300 : DO ipgf = 1, orb_basis%npgf(iset) ! ipgf
148 11274 : IF (orb_basis%zet(ipgf, iset) > 100.0_dp) THEN
149 : ! The soft basis does not contain this exponent
150 678 : paw_atom = .TRUE.
151 678 : CYCLE
152 : END IF
153 :
154 : radius = exp_radius(orb_basis%lmax(iset), orb_basis%zet(ipgf, iset), &
155 10596 : eps_fit, 1.0_dp)
156 10596 : IF (radius < rc) THEN
157 : ! The soft basis does not contain this exponent
158 2266 : paw_atom = .TRUE.
159 2266 : CYCLE
160 : END IF
161 :
162 : ! The soft basis contains this exponent
163 8330 : ipgf_s = ipgf_s + 1
164 8330 : zet(ipgf_s, iset_s) = orb_basis%zet(ipgf, iset)
165 :
166 8330 : lshell_old = orb_basis%l(1, iset)
167 8330 : radius = exp_radius(lshell_old, zet(ipgf_s, iset_s), eps_fit, 1.0_dp)
168 :
169 31470 : DO ishell = 1, nshell(iset_s)
170 18114 : lshell = orb_basis%l(ishell, iset)
171 18114 : IF (lshell == lshell_old) THEN
172 : ELSE
173 3346 : lshell_old = lshell
174 3346 : radius = exp_radius(lshell_old, zet(ipgf_s, iset_s), eps_fit, 1.0_dp)
175 : END IF
176 26444 : IF (radius < rc) THEN
177 4 : gcc(ipgf_s, ishell, iset_s) = 0.0_dp
178 4 : paw_atom = .TRUE.
179 : ELSE
180 18110 : gcc(ipgf_s, ishell, iset_s) = orb_basis%gcc(ipgf, ishell, iset)
181 : END IF
182 : END DO
183 : END DO
184 5026 : npgf(iset_s) = ipgf_s
185 6930 : IF (ipgf_s > maxpgf_s) THEN
186 2132 : maxpgf_s = ipgf_s
187 : END IF
188 : END DO
189 1904 : nset_s = iset_s
190 1904 : IF (paw_atom) THEN
191 1620 : soft_basis%nset = nset_s
192 1620 : CALL reallocate(soft_basis%lmax, 1, nset_s)
193 1620 : CALL reallocate(soft_basis%lmin, 1, nset_s)
194 1620 : CALL reallocate(soft_basis%npgf, 1, nset_s)
195 1620 : CALL reallocate(soft_basis%nshell, 1, nset_s)
196 1620 : CALL reallocate(soft_basis%n, 1, maxshell_s, 1, nset_s)
197 1620 : CALL reallocate(soft_basis%l, 1, maxshell_s, 1, nset_s)
198 1620 : CALL reallocate(soft_basis%zet, 1, maxpgf_s, 1, nset_s)
199 1620 : CALL reallocate(soft_basis%gcc, 1, maxpgf_s, 1, maxshell_s, 1, nset_s)
200 :
201 : ! *** Copy the basis set information into the data structure ***
202 :
203 6654 : DO iset = 1, nset_s
204 5034 : soft_basis%lmax(iset) = lmax(iset)
205 5034 : soft_basis%lmin(iset) = lmin(iset)
206 5034 : soft_basis%npgf(iset) = npgf(iset)
207 5034 : soft_basis%nshell(iset) = nshell(iset)
208 13014 : DO ishell = 1, nshell(iset)
209 7980 : soft_basis%n(ishell, iset) = n(ishell, iset)
210 7980 : soft_basis%l(ishell, iset) = l(ishell, iset)
211 28552 : DO ipgf = 1, npgf(iset)
212 23518 : soft_basis%gcc(ipgf, ishell, iset) = gcc(ipgf, ishell, iset)
213 : END DO
214 : END DO
215 13992 : DO ipgf = 1, npgf(iset)
216 12372 : soft_basis%zet(ipgf, iset) = zet(ipgf, iset)
217 : END DO
218 : END DO
219 :
220 : ! *** Initialise the depending soft_basis pointers ***
221 1620 : CALL reallocate(soft_basis%set_radius, 1, nset_s)
222 1620 : CALL reallocate(soft_basis%pgf_radius, 1, maxpgf_s, 1, nset_s)
223 1620 : CALL reallocate(soft_basis%first_cgf, 1, maxshell_s, 1, nset_s)
224 1620 : CALL reallocate(soft_basis%first_sgf, 1, maxshell_s, 1, nset_s)
225 1620 : CALL reallocate(soft_basis%last_cgf, 1, maxshell_s, 1, nset_s)
226 1620 : CALL reallocate(soft_basis%last_sgf, 1, maxshell_s, 1, nset_s)
227 1620 : CALL reallocate(soft_basis%ncgf_set, 1, nset_s)
228 1620 : CALL reallocate(soft_basis%nsgf_set, 1, nset_s)
229 :
230 1620 : maxco = 0
231 1620 : ncgf = 0
232 1620 : nsgf = 0
233 :
234 6654 : DO iset = 1, nset_s
235 5034 : soft_basis%ncgf_set(iset) = 0
236 5034 : soft_basis%nsgf_set(iset) = 0
237 13014 : DO ishell = 1, nshell(iset)
238 7980 : lshell = soft_basis%l(ishell, iset)
239 7980 : soft_basis%first_cgf(ishell, iset) = ncgf + 1
240 7980 : ncgf = ncgf + nco(lshell)
241 7980 : soft_basis%last_cgf(ishell, iset) = ncgf
242 : soft_basis%ncgf_set(iset) = &
243 7980 : soft_basis%ncgf_set(iset) + nco(lshell)
244 7980 : soft_basis%first_sgf(ishell, iset) = nsgf + 1
245 7980 : nsgf = nsgf + nso(lshell)
246 7980 : soft_basis%last_sgf(ishell, iset) = nsgf
247 : soft_basis%nsgf_set(iset) = &
248 13014 : soft_basis%nsgf_set(iset) + nso(lshell)
249 : END DO
250 6654 : maxco = MAX(maxco, npgf(iset)*ncoset(lmax(iset)))
251 : END DO
252 1620 : soft_basis%ncgf = ncgf
253 1620 : soft_basis%nsgf = nsgf
254 :
255 1620 : CALL reallocate(soft_basis%cphi, 1, maxco, 1, ncgf)
256 1620 : CALL reallocate(soft_basis%sphi, 1, maxco, 1, nsgf)
257 1620 : CALL reallocate(soft_basis%scon, 1, maxco, 1, nsgf)
258 1620 : CALL reallocate(soft_basis%lx, 1, ncgf)
259 1620 : CALL reallocate(soft_basis%ly, 1, ncgf)
260 1620 : CALL reallocate(soft_basis%lz, 1, ncgf)
261 1620 : CALL reallocate(soft_basis%m, 1, nsgf)
262 1620 : CALL reallocate(soft_basis%norm_cgf, 1, ncgf)
263 4860 : ALLOCATE (soft_basis%cgf_symbol(ncgf))
264 4860 : ALLOCATE (soft_basis%sgf_symbol(nsgf))
265 :
266 1620 : ncgf = 0
267 1620 : nsgf = 0
268 :
269 6654 : DO iset = 1, nset_s
270 14634 : DO ishell = 1, nshell(iset)
271 7980 : lshell = soft_basis%l(ishell, iset)
272 26838 : DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
273 18858 : ncgf = ncgf + 1
274 18858 : soft_basis%lx(ncgf) = indco(1, ico)
275 18858 : soft_basis%ly(ncgf) = indco(2, ico)
276 18858 : soft_basis%lz(ncgf) = indco(3, ico)
277 : soft_basis%cgf_symbol(ncgf) = &
278 : cgf_symbol(n(ishell, iset), (/soft_basis%lx(ncgf), &
279 : soft_basis%ly(ncgf), &
280 83412 : soft_basis%lz(ncgf)/))
281 : END DO
282 30698 : DO m = -lshell, lshell
283 17684 : nsgf = nsgf + 1
284 17684 : soft_basis%m(nsgf) = m
285 : soft_basis%sgf_symbol(nsgf) = &
286 25664 : sgf_symbol(n(ishell, iset), lshell, m)
287 : END DO
288 : END DO
289 : END DO
290 :
291 : ! *** Normalization factor of the contracted Gaussians ***
292 1620 : soft_basis%norm_type = orb_basis%norm_type
293 39336 : soft_basis%norm_cgf = orb_basis%norm_cgf
294 : ! *** Initialize the transformation matrices ***
295 1620 : CALL init_cphi_and_sphi(soft_basis)
296 : END IF
297 :
298 3808 : DEALLOCATE (gcc, l, lmax, lmin, n, npgf, nshell, zet, iset_s2h)
299 : END IF
300 :
301 1932 : IF (.NOT. paw_atom) THEN
302 312 : CALL deallocate_gto_basis_set(soft_basis)
303 312 : CALL copy_gto_basis_set(orb_basis, soft_basis)
304 312 : CALL get_gto_basis_set(gto_basis_set=orb_basis, name=bsname)
305 312 : soft_basis%name = TRIM(bsname)//"_soft"
306 : END IF
307 :
308 1932 : END SUBROUTINE create_soft_basis
309 :
310 : END MODULE soft_basis_set
|
