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 Several screening methods used in HFX calcualtions
10 : !> \par History
11 : !> 04.2008 created [Manuel Guidon]
12 : !> \author Manuel Guidon
13 : ! **************************************************************************************************
14 : MODULE hfx_screening_methods
15 : USE ao_util, ONLY: exp_radius_very_extended
16 : USE basis_set_types, ONLY: gto_basis_set_type
17 : USE four_center_integral_interface, ONLY: evaluate_eri_screen,&
18 : four_center_integral_context_type
19 : USE hfx_types, ONLY: hfx_basis_type,&
20 : hfx_p_kind,&
21 : hfx_potential_type,&
22 : hfx_screen_coeff_type,&
23 : log_zero,&
24 : powell_min_log
25 : USE kinds, ONLY: dp,&
26 : int_8
27 : USE machine, ONLY: default_output_unit
28 : USE orbital_pointers, ONLY: ncoset
29 : USE powell, ONLY: opt_state_type,&
30 : powell_optimize
31 : USE qs_environment_types, ONLY: get_qs_env,&
32 : qs_environment_type
33 : USE qs_kind_types, ONLY: get_qs_kind,&
34 : qs_kind_type
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 : PRIVATE
39 :
40 : PUBLIC :: update_pmax_mat, &
41 : calc_screening_functions, &
42 : calc_pair_dist_radii
43 :
44 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_screening_methods'
45 :
46 : !***
47 :
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief calculates max values of two-electron integrals in a quartet/shell
52 : !> w.r.t. different zetas using the library lib_int
53 : !> \param context four-center integral context
54 : !> \param ra position
55 : !> \param rb position
56 : !> \param zeta zeta
57 : !> \param zetb zeta
58 : !> \param la_min angular momentum
59 : !> \param la_max angular momentum
60 : !> \param lb_min angular momentum
61 : !> \param lb_max angular momentum
62 : !> \param npgfa number of primitive cartesian gaussian in actual shell
63 : !> \param npgfb number of primitive cartesian gaussian in actual shell
64 : !> \param max_val schwarz screening value
65 : !> \param potential_parameter radial operator parameters
66 : !> \param tmp_R_1 pgf_based screening coefficients
67 : !> \param rab2 Squared Distance of centers ab
68 : !> \par History
69 : !> 03.2007 created [Manuel Guidon]
70 : !> 02.2009 refactored [Manuel Guidon]
71 : !> \author Manuel Guidon
72 : ! **************************************************************************************************
73 10106464 : SUBROUTINE screen4(context, ra, rb, zeta, zetb, &
74 : la_min, la_max, lb_min, lb_max, &
75 : npgfa, npgfb, &
76 : max_val, potential_parameter, tmp_R_1, &
77 : rab2)
78 :
79 : TYPE(four_center_integral_context_type) :: context
80 : REAL(dp), INTENT(IN) :: ra(3), rb(3)
81 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, zetb
82 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, npgfa, &
83 : npgfb
84 : REAL(dp), INTENT(INOUT) :: max_val
85 : TYPE(hfx_potential_type) :: potential_parameter
86 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
87 : POINTER :: tmp_R_1
88 : REAL(dp) :: rab2
89 :
90 : INTEGER :: ipgf, jpgf, la, lb
91 : REAL(dp) :: max_val_temp, R1
92 :
93 10106464 : max_val_temp = max_val
94 27398472 : DO ipgf = 1, npgfa
95 62889468 : DO jpgf = 1, npgfb
96 35490996 : R1 = MAX(0.0_dp, tmp_R_1(jpgf, ipgf)%x(1)*rab2 + tmp_R_1(jpgf, ipgf)%x(2))
97 99412684 : DO la = la_min, la_max
98 146375664 : DO lb = lb_min, lb_max
99 : !Build primitives
100 : CALL evaluate_eri_screen(context, ra, rb, ra, rb, &
101 : zeta(ipgf), zetb(jpgf), zeta(ipgf), zetb(jpgf), &
102 : la, lb, la, lb, &
103 64254988 : max_val_temp, potential_parameter, R1, R1)
104 :
105 110884668 : max_val = MAX(max_val, max_val_temp)
106 : END DO !lb
107 : END DO !la
108 : END DO !jpgf
109 : END DO !ipgf
110 10106464 : END SUBROUTINE screen4
111 :
112 : ! **************************************************************************************************
113 : !> \brief updates the maximum of the density matrix in compressed form for screening purposes
114 : !> \param pmax_set buffer to store matrix
115 : !> \param map_atom_to_kind_atom ...
116 : !> \param set_offset ...
117 : !> \param atomic_block_offset ...
118 : !> \param pmax_atom ...
119 : !> \param full_density_alpha ...
120 : !> \param full_density_beta ...
121 : !> \param natom ...
122 : !> \param kind_of ...
123 : !> \param basis_parameter ...
124 : !> \param nkind ...
125 : !> \param is_assoc_atomic_block ...
126 : !> \param nspins number of spin channels, used to express bounds in restricted-density units
127 : !> \par History
128 : !> 09.2007 created [Manuel Guidon]
129 : !> \author Manuel Guidon
130 : !> \note
131 : !> - stores log10(nspins*max(abs(p_alpha),abs(p_beta))) for each shell pair,
132 : !> so equivalent restricted and unrestricted closed-shell densities have the same bound
133 : ! **************************************************************************************************
134 1612 : SUBROUTINE update_pmax_mat(pmax_set, map_atom_to_kind_atom, set_offset, atomic_block_offset, &
135 : pmax_atom, full_density_alpha, full_density_beta, natom, &
136 : kind_of, basis_parameter, &
137 806 : nkind, is_assoc_atomic_block, nspins)
138 :
139 : TYPE(hfx_p_kind), DIMENSION(:), POINTER :: pmax_set
140 : INTEGER, DIMENSION(:), POINTER :: map_atom_to_kind_atom
141 : INTEGER, DIMENSION(:, :, :, :), POINTER :: set_offset
142 : INTEGER, DIMENSION(:, :), POINTER :: atomic_block_offset
143 : REAL(dp), DIMENSION(:, :), POINTER :: pmax_atom, full_density_alpha, &
144 : full_density_beta
145 : INTEGER, INTENT(IN) :: natom
146 : INTEGER :: kind_of(*)
147 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
148 : INTEGER :: nkind
149 : INTEGER, DIMENSION(:, :) :: is_assoc_atomic_block
150 : INTEGER, INTENT(IN) :: nspins
151 :
152 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_pmax_mat'
153 :
154 : INTEGER :: act_atomic_block_offset, act_set_offset, handle, i, img, jatom, jkind, jset, &
155 : katom, kind_kind_idx, kkind, kset, mb, mc, nimg, nsetb, nsetc
156 806 : INTEGER, DIMENSION(:), POINTER :: nsgfb, nsgfc
157 : REAL(dp) :: pmax_tmp
158 :
159 806 : CALL timeset(routineN, handle)
160 :
161 3124 : DO i = 1, SIZE(pmax_set)
162 78968 : pmax_set(i)%p_kind = 0.0_dp
163 : END DO
164 :
165 11838 : pmax_atom = log_zero
166 :
167 806 : nimg = SIZE(full_density_alpha, 2)
168 :
169 3316 : DO jatom = 1, natom
170 2510 : jkind = kind_of(jatom)
171 2510 : nsetb = basis_parameter(jkind)%nset
172 2510 : nsgfb => basis_parameter(jkind)%nsgf
173 :
174 11838 : DO katom = 1, natom
175 8522 : IF (.NOT. is_assoc_atomic_block(jatom, katom) >= 1) CYCLE
176 8474 : kkind = kind_of(katom)
177 8474 : IF (kkind < jkind) CYCLE
178 6768 : nsetc = basis_parameter(kkind)%nset
179 6768 : nsgfc => basis_parameter(kkind)%nsgf
180 6768 : act_atomic_block_offset = atomic_block_offset(katom, jatom)
181 26562 : DO jset = 1, nsetb
182 74770 : DO kset = 1, nsetc
183 66248 : IF (katom >= jatom) THEN
184 42190 : pmax_tmp = 0.0_dp
185 42190 : act_set_offset = set_offset(kset, jset, kkind, jkind)
186 84380 : DO img = 1, nimg
187 42190 : i = act_set_offset + act_atomic_block_offset - 1
188 164582 : DO mc = 1, nsgfc(kset)
189 347264 : DO mb = 1, nsgfb(jset)
190 224872 : pmax_tmp = MAX(pmax_tmp, ABS(full_density_alpha(i, img)))
191 224872 : IF (ASSOCIATED(full_density_beta)) THEN
192 42600 : pmax_tmp = MAX(pmax_tmp, ABS(full_density_beta(i, img)))
193 : END IF
194 305074 : i = i + 1
195 : END DO
196 : END DO
197 : END DO
198 42190 : IF (pmax_tmp == 0.0_dp) THEN
199 : pmax_tmp = log_zero
200 : ELSE
201 41964 : pmax_tmp = LOG10(REAL(nspins, dp)*pmax_tmp)
202 : END IF
203 42190 : kind_kind_idx = INT(get_1D_idx(jkind, kkind, INT(nkind, int_8)))
204 : pmax_set(kind_kind_idx)%p_kind(jset, &
205 : kset, &
206 : map_atom_to_kind_atom(jatom), &
207 42190 : map_atom_to_kind_atom(katom)) = pmax_tmp
208 : ELSE
209 6774 : pmax_tmp = 0.0_dp
210 6774 : act_set_offset = set_offset(jset, kset, jkind, kkind)
211 13548 : DO img = 1, nimg
212 26252 : DO mc = 1, nsgfc(kset)
213 12704 : i = act_set_offset + act_atomic_block_offset - 1 + mc - 1
214 48444 : DO mb = 1, nsgfb(jset)
215 28966 : pmax_tmp = MAX(pmax_tmp, ABS(full_density_alpha(i, img)))
216 28966 : IF (ASSOCIATED(full_density_beta)) THEN
217 3010 : pmax_tmp = MAX(pmax_tmp, ABS(full_density_beta(i, img)))
218 : END IF
219 41670 : i = i + nsgfc(kset)
220 : END DO
221 : END DO
222 : END DO
223 6774 : IF (pmax_tmp == 0.0_dp) THEN
224 : pmax_tmp = log_zero
225 : ELSE
226 6748 : pmax_tmp = LOG10(REAL(nspins, dp)*pmax_tmp)
227 : END IF
228 6774 : kind_kind_idx = INT(get_1D_idx(jkind, kkind, INT(nkind, int_8)))
229 : pmax_set(kind_kind_idx)%p_kind(jset, &
230 : kset, &
231 : map_atom_to_kind_atom(jatom), &
232 6774 : map_atom_to_kind_atom(katom)) = pmax_tmp
233 : END IF
234 : END DO
235 : END DO
236 : END DO
237 : END DO
238 :
239 3316 : DO jatom = 1, natom
240 2510 : jkind = kind_of(jatom)
241 2510 : nsetb = basis_parameter(jkind)%nset
242 11838 : DO katom = 1, natom
243 8522 : IF (.NOT. is_assoc_atomic_block(jatom, katom) >= 1) CYCLE
244 8474 : kkind = kind_of(katom)
245 8474 : IF (kkind < jkind) CYCLE
246 6768 : nsetc = basis_parameter(kkind)%nset
247 6768 : pmax_tmp = log_zero
248 24052 : DO jset = 1, nsetb
249 73016 : DO kset = 1, nsetc
250 48964 : kind_kind_idx = INT(get_1D_idx(jkind, kkind, INT(nkind, int_8)))
251 : pmax_tmp = MAX(pmax_set(kind_kind_idx)%p_kind(jset, &
252 : kset, &
253 : map_atom_to_kind_atom(jatom), &
254 66248 : map_atom_to_kind_atom(katom)), pmax_tmp)
255 : END DO
256 : END DO
257 6768 : pmax_atom(jatom, katom) = pmax_tmp
258 11032 : pmax_atom(katom, jatom) = pmax_tmp
259 : END DO
260 : END DO
261 :
262 806 : CALL timestop(handle)
263 :
264 806 : END SUBROUTINE update_pmax_mat
265 :
266 : ! **************************************************************************************************
267 : !> \brief calculates screening functions for schwarz screening
268 : !> \param qs_env qs_env
269 : !> \param basis_parameter ...
270 : !> \param context structure to libint
271 : !> \param potential_parameter contains infos on potential
272 : !> \param coeffs_set set based coefficients
273 : !> \param coeffs_kind kind based coefficients
274 : !> \param coeffs_pgf pgf based coefficients
275 : !> \param radii_pgf coefficients for long-range screening
276 : !> \param max_set Maximum Number of basis set sets in the system
277 : !> \param max_pgf Maximum Number of basis set pgfs in the system
278 : !> \param n_threads ...
279 : !> \param i_thread Thread ID of current task
280 : !> \par History
281 : !> 02.2009 created [Manuel Guidon]
282 : !> \author Manuel Guidon
283 : !> \note
284 : !> This routine calculates (ab|ab) for different distances Rab = |a-b|
285 : !> and uses the powell optimiztion routines in order to fit the results
286 : !> in the following form:
287 : !>
288 : !> (ab|ab) = (ab|ab)(Rab) = c2*Rab^2 + c0
289 : !>
290 : !> The missing linear term assures that the functions is monotonically
291 : !> decaying such that c2 can be used as upper bound when applying the
292 : !> Minimum Image Convention in the periodic case. Furthermore
293 : !> it seems to be a good choice to fit the logarithm of the (ab|ab)
294 : !> The fitting takes place at several levels: kind, set and pgf within
295 : !> the corresponding ranges of the prodiuct charge distributions
296 : !> Doing so, we only need arrays of size nkinds^2*2 instead of big
297 : !> screening matrices
298 : ! **************************************************************************************************
299 :
300 1490 : SUBROUTINE calc_screening_functions(qs_env, basis_parameter, context, potential_parameter, &
301 : coeffs_set, coeffs_kind, coeffs_pgf, radii_pgf, &
302 : max_set, max_pgf, n_threads, i_thread)
303 : TYPE(qs_environment_type), POINTER :: qs_env
304 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
305 : TYPE(four_center_integral_context_type) :: context
306 : TYPE(hfx_potential_type) :: potential_parameter
307 : TYPE(hfx_screen_coeff_type), &
308 : DIMENSION(:, :, :, :), POINTER :: coeffs_set
309 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
310 : POINTER :: coeffs_kind
311 : TYPE(hfx_screen_coeff_type), &
312 : DIMENSION(:, :, :, :, :, :), POINTER :: coeffs_pgf, radii_pgf
313 : INTEGER, INTENT(IN) :: max_set, max_pgf, n_threads, i_thread
314 :
315 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_screening_functions'
316 :
317 : INTEGER :: handle, i, ikind, ipgf, iset, jkind, &
318 : jpgf, jset, la, lb, ncoa, ncob, nkind, &
319 : nseta, nsetb, sgfa, sgfb
320 1490 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
321 1490 : npgfb, nsgfa, nsgfb
322 1490 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
323 : REAL(dp) :: kind_radius_a, kind_radius_b, max_contraction_a, max_contraction_b, max_val, &
324 : max_val_temp, R1, ra(3), radius, rb(3), x(2)
325 1490 : REAL(dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
326 1490 : REAL(dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
327 : REAL(dp), SAVE :: DATA(2, 0:100)
328 : TYPE(gto_basis_set_type), POINTER :: orb_basis
329 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
330 1490 : POINTER :: tmp_R_1
331 1490 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
332 : TYPE(qs_kind_type), POINTER :: qs_kind
333 :
334 2980 : !$OMP MASTER
335 1490 : CALL timeset(routineN, handle)
336 : !$OMP END MASTER
337 :
338 : CALL get_qs_env(qs_env=qs_env, &
339 1490 : qs_kind_set=qs_kind_set)
340 :
341 1490 : nkind = SIZE(qs_kind_set, 1)
342 :
343 1490 : !$OMP MASTER
344 1387872 : ALLOCATE (coeffs_pgf(max_pgf, max_pgf, max_set, max_set, nkind, nkind))
345 :
346 4184 : DO ikind = 1, nkind
347 10366 : DO jkind = 1, nkind
348 26490 : DO iset = 1, max_set
349 94686 : DO jset = 1, max_set
350 313240 : DO ipgf = 1, max_pgf
351 1346482 : DO jpgf = 1, max_pgf
352 3377304 : coeffs_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(:) = 0.0_dp
353 : END DO
354 : END DO
355 : END DO
356 : END DO
357 : END DO
358 : END DO
359 :
360 : !$OMP END MASTER
361 1490 : !$OMP BARRIER
362 1490 : ra = 0.0_dp
363 1490 : rb = 0.0_dp
364 4184 : DO ikind = 1, nkind
365 2694 : NULLIFY (qs_kind, orb_basis)
366 2694 : qs_kind => qs_kind_set(ikind)
367 2694 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
368 2694 : NULLIFY (la_max, la_min, npgfa, zeta)
369 :
370 2694 : la_max => basis_parameter(ikind)%lmax
371 2694 : la_min => basis_parameter(ikind)%lmin
372 2694 : npgfa => basis_parameter(ikind)%npgf
373 2694 : nseta = basis_parameter(ikind)%nset
374 2694 : zeta => basis_parameter(ikind)%zet
375 2694 : set_radius_a => basis_parameter(ikind)%set_radius
376 2694 : first_sgfa => basis_parameter(ikind)%first_sgf
377 2694 : sphi_a => basis_parameter(ikind)%sphi
378 2694 : nsgfa => basis_parameter(ikind)%nsgf
379 2694 : rpgfa => basis_parameter(ikind)%pgf_radius
380 :
381 10366 : DO jkind = 1, nkind
382 6182 : NULLIFY (qs_kind, orb_basis)
383 6182 : qs_kind => qs_kind_set(jkind)
384 6182 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
385 6182 : NULLIFY (lb_max, lb_min, npgfb, zetb)
386 :
387 6182 : lb_max => basis_parameter(jkind)%lmax
388 6182 : lb_min => basis_parameter(jkind)%lmin
389 6182 : npgfb => basis_parameter(jkind)%npgf
390 6182 : nsetb = basis_parameter(jkind)%nset
391 6182 : zetb => basis_parameter(jkind)%zet
392 6182 : set_radius_b => basis_parameter(jkind)%set_radius
393 6182 : first_sgfb => basis_parameter(jkind)%first_sgf
394 6182 : sphi_b => basis_parameter(jkind)%sphi
395 6182 : nsgfb => basis_parameter(jkind)%nsgf
396 6182 : rpgfb => basis_parameter(jkind)%pgf_radius
397 :
398 24180 : DO iset = 1, nseta
399 15304 : ncoa = npgfa(iset)*ncoset(la_max(iset))
400 15304 : sgfa = first_sgfa(1, iset)
401 578060 : max_contraction_a = MAXVAL([(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)])
402 71518 : DO jset = 1, nsetb
403 50032 : ncob = npgfb(jset)*ncoset(lb_max(jset))
404 50032 : sgfb = first_sgfb(1, jset)
405 1301730 : max_contraction_b = MAXVAL([(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)])
406 50032 : radius = set_radius_a(iset) + set_radius_b(jset)
407 150940 : DO ipgf = 1, npgfa(iset)
408 311334 : DO jpgf = 1, npgfb(jset)
409 175698 : radius = rpgfa(ipgf, iset) + rpgfb(jpgf, jset)
410 17921196 : DO i = i_thread, 100, n_threads
411 17745498 : rb(1) = 0.0_dp + REAL(i, dp)*0.01_dp*radius
412 17745498 : max_val = 0.0_dp
413 : R1 = MAX(0.0_dp, radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(1)*rb(1)**2 + &
414 17745498 : radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(2))
415 41060338 : DO la = la_min(iset), la_max(iset)
416 73187832 : DO lb = lb_min(jset), lb_max(jset)
417 : !Build primitives
418 32127494 : max_val_temp = 0.0_dp
419 : CALL evaluate_eri_screen(context, ra, rb, ra, rb, &
420 : zeta(ipgf, iset), zetb(jpgf, jset), zeta(ipgf, iset), zetb(jpgf, jset), &
421 : la, lb, la, lb, &
422 32127494 : max_val_temp, potential_parameter, R1, R1)
423 55442334 : max_val = MAX(max_val, max_val_temp)
424 : END DO !lb
425 : END DO !la
426 17745498 : max_val = SQRT(max_val)
427 17745498 : max_val = max_val*max_contraction_a*max_contraction_b
428 17745498 : DATA(1, i) = rb(1)
429 17921196 : IF (max_val == 0.0_dp) THEN
430 0 : DATA(2, i) = powell_min_log
431 : ELSE
432 17745498 : DATA(2, i) = LOG10((max_val))
433 : END IF
434 : END DO
435 175698 : !$OMP BARRIER
436 175698 : !$OMP MASTER
437 175698 : CALL optimize_it(DATA, x, powell_min_log)
438 527094 : coeffs_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x = x
439 : !$OMP END MASTER
440 261302 : !$OMP BARRIER
441 :
442 : END DO
443 : END DO
444 : END DO
445 : END DO
446 : END DO
447 : END DO
448 :
449 1490 : !$OMP MASTER
450 109300 : ALLOCATE (coeffs_set(max_set, max_set, nkind, nkind))
451 :
452 4184 : DO ikind = 1, nkind
453 10366 : DO jkind = 1, nkind
454 26490 : DO iset = 1, max_set
455 94686 : DO jset = 1, max_set
456 230284 : coeffs_set(jset, iset, jkind, ikind)%x(:) = 0.0_dp
457 : END DO
458 : END DO
459 : END DO
460 : END DO
461 : !$OMP END MASTER
462 1490 : !$OMP BARRIER
463 1490 : ra = 0.0_dp
464 1490 : rb = 0.0_dp
465 4184 : DO ikind = 1, nkind
466 2694 : NULLIFY (qs_kind, orb_basis)
467 2694 : qs_kind => qs_kind_set(ikind)
468 2694 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
469 2694 : NULLIFY (la_max, la_min, npgfa, zeta)
470 : ! CALL get_gto_basis_set(gto_basis_set=orb_basis,&
471 : ! lmax=la_max,&
472 : ! lmin=la_min,&
473 : ! npgf=npgfa,&
474 : ! nset=nseta,&
475 : ! zet=zeta,&
476 : ! set_radius=set_radius_a,&
477 : ! first_sgf=first_sgfa,&
478 : ! sphi=sphi_a,&
479 : ! nsgf_set=nsgfa)
480 2694 : la_max => basis_parameter(ikind)%lmax
481 2694 : la_min => basis_parameter(ikind)%lmin
482 2694 : npgfa => basis_parameter(ikind)%npgf
483 2694 : nseta = basis_parameter(ikind)%nset
484 2694 : zeta => basis_parameter(ikind)%zet
485 2694 : set_radius_a => basis_parameter(ikind)%set_radius
486 2694 : first_sgfa => basis_parameter(ikind)%first_sgf
487 2694 : sphi_a => basis_parameter(ikind)%sphi
488 2694 : nsgfa => basis_parameter(ikind)%nsgf
489 :
490 10366 : DO jkind = 1, nkind
491 6182 : NULLIFY (qs_kind, orb_basis)
492 6182 : qs_kind => qs_kind_set(jkind)
493 6182 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
494 6182 : NULLIFY (lb_max, lb_min, npgfb, zetb)
495 :
496 6182 : lb_max => basis_parameter(jkind)%lmax
497 6182 : lb_min => basis_parameter(jkind)%lmin
498 6182 : npgfb => basis_parameter(jkind)%npgf
499 6182 : nsetb = basis_parameter(jkind)%nset
500 6182 : zetb => basis_parameter(jkind)%zet
501 6182 : set_radius_b => basis_parameter(jkind)%set_radius
502 6182 : first_sgfb => basis_parameter(jkind)%first_sgf
503 6182 : sphi_b => basis_parameter(jkind)%sphi
504 6182 : nsgfb => basis_parameter(jkind)%nsgf
505 :
506 24180 : DO iset = 1, nseta
507 15304 : ncoa = npgfa(iset)*ncoset(la_max(iset))
508 15304 : sgfa = first_sgfa(1, iset)
509 578060 : max_contraction_a = MAXVAL([(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)])
510 71518 : DO jset = 1, nsetb
511 50032 : ncob = npgfb(jset)*ncoset(lb_max(jset))
512 50032 : sgfb = first_sgfb(1, jset)
513 1301730 : max_contraction_b = MAXVAL([(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)])
514 50032 : radius = set_radius_a(iset) + set_radius_b(jset)
515 50032 : tmp_R_1 => radii_pgf(:, :, jset, iset, jkind, ikind)
516 5103264 : DO i = i_thread, 100, n_threads
517 5053232 : rb(1) = 0.0_dp + REAL(i, dp)*0.01_dp*radius
518 5053232 : max_val = 0.0_dp
519 : CALL screen4(context, ra, rb, &
520 : zeta(:, iset), zetb(:, jset), &
521 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
522 : npgfa(iset), npgfb(jset), &
523 5053232 : max_val, potential_parameter, tmp_R_1, rb(1)**2)
524 5053232 : max_val = SQRT(max_val)
525 5053232 : max_val = max_val*max_contraction_a*max_contraction_b
526 5053232 : DATA(1, i) = rb(1)
527 5103264 : IF (max_val == 0.0_dp) THEN
528 0 : DATA(2, i) = powell_min_log
529 : ELSE
530 5053232 : DATA(2, i) = LOG10((max_val))
531 : END IF
532 : END DO
533 50032 : !$OMP BARRIER
534 50032 : !$OMP MASTER
535 50032 : CALL optimize_it(DATA, x, powell_min_log)
536 150096 : coeffs_set(jset, iset, jkind, ikind)%x = x
537 : !$OMP END MASTER
538 65336 : !$OMP BARRIER
539 : END DO
540 : END DO
541 :
542 : END DO
543 : END DO
544 :
545 : ! ** now kinds
546 1490 : !$OMP MASTER
547 17816 : ALLOCATE (coeffs_kind(nkind, nkind))
548 :
549 4184 : DO ikind = 1, nkind
550 10366 : DO jkind = 1, nkind
551 21240 : coeffs_kind(jkind, ikind)%x(:) = 0.0_dp
552 : END DO
553 : END DO
554 : !$OMP END MASTER
555 1490 : ra = 0.0_dp
556 1490 : rb = 0.0_dp
557 4184 : DO ikind = 1, nkind
558 2694 : NULLIFY (qs_kind, orb_basis)
559 2694 : qs_kind => qs_kind_set(ikind)
560 2694 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
561 2694 : NULLIFY (la_max, la_min, npgfa, zeta)
562 :
563 2694 : la_max => basis_parameter(ikind)%lmax
564 2694 : la_min => basis_parameter(ikind)%lmin
565 2694 : npgfa => basis_parameter(ikind)%npgf
566 2694 : nseta = basis_parameter(ikind)%nset
567 2694 : zeta => basis_parameter(ikind)%zet
568 2694 : kind_radius_a = basis_parameter(ikind)%kind_radius
569 2694 : first_sgfa => basis_parameter(ikind)%first_sgf
570 2694 : sphi_a => basis_parameter(ikind)%sphi
571 2694 : nsgfa => basis_parameter(ikind)%nsgf
572 :
573 10366 : DO jkind = 1, nkind
574 6182 : NULLIFY (qs_kind, orb_basis)
575 6182 : qs_kind => qs_kind_set(jkind)
576 6182 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
577 6182 : NULLIFY (lb_max, lb_min, npgfb, zetb)
578 :
579 6182 : lb_max => basis_parameter(jkind)%lmax
580 6182 : lb_min => basis_parameter(jkind)%lmin
581 6182 : npgfb => basis_parameter(jkind)%npgf
582 6182 : nsetb = basis_parameter(jkind)%nset
583 6182 : zetb => basis_parameter(jkind)%zet
584 6182 : kind_radius_b = basis_parameter(jkind)%kind_radius
585 6182 : first_sgfb => basis_parameter(jkind)%first_sgf
586 6182 : sphi_b => basis_parameter(jkind)%sphi
587 6182 : nsgfb => basis_parameter(jkind)%nsgf
588 :
589 6182 : radius = kind_radius_a + kind_radius_b
590 21486 : DO iset = 1, nseta
591 15304 : ncoa = npgfa(iset)*ncoset(la_max(iset))
592 15304 : sgfa = first_sgfa(1, iset)
593 578060 : max_contraction_a = MAXVAL([(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)])
594 71518 : DO jset = 1, nsetb
595 50032 : ncob = npgfb(jset)*ncoset(lb_max(jset))
596 50032 : sgfb = first_sgfb(1, jset)
597 1301730 : max_contraction_b = MAXVAL([(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)])
598 5118568 : DO i = i_thread, 100, n_threads
599 5053232 : rb(1) = 0.0_dp + REAL(i, dp)*0.01_dp*radius
600 5053232 : max_val = 0.0_dp
601 5053232 : tmp_R_1 => radii_pgf(:, :, jset, iset, jkind, ikind)
602 : CALL screen4(context, ra, rb, &
603 : zeta(:, iset), zetb(:, jset), &
604 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
605 : npgfa(iset), npgfb(jset), &
606 5053232 : max_val, potential_parameter, tmp_R_1, rb(1)**2)
607 5053232 : DATA(1, i) = rb(1)
608 5053232 : max_val = SQRT(max_val)
609 5053232 : max_val = max_val*max_contraction_a*max_contraction_b
610 5103264 : IF (max_val == 0.0_dp) THEN
611 153394 : DATA(2, i) = MAX(powell_min_log, DATA(2, i))
612 : ELSE
613 4899838 : DATA(2, i) = MAX(LOG10(max_val), DATA(2, i))
614 : END IF
615 : END DO
616 : END DO
617 : END DO
618 6182 : !$OMP BARRIER
619 6182 : !$OMP MASTER
620 6182 : CALL optimize_it(DATA, x, powell_min_log)
621 18546 : coeffs_kind(jkind, ikind)%x = x
622 : !$OMP END MASTER
623 8876 : !$OMP BARRIER
624 : END DO
625 : END DO
626 :
627 1490 : !$OMP MASTER
628 1490 : CALL timestop(handle)
629 : !$OMP END MASTER
630 :
631 1490 : END SUBROUTINE calc_screening_functions
632 :
633 : ! **************************************************************************************************
634 : !> \brief calculates radius functions for longrange screening
635 : !> \param qs_env qs_env
636 : !> \param basis_parameter ...
637 : !> \param radii_pgf pgf based coefficients
638 : !> \param max_set Maximum Number of basis set sets in the system
639 : !> \param max_pgf Maximum Number of basis set pgfs in the system
640 : !> \param eps_schwarz ...
641 : !> \param n_threads ...
642 : !> \param i_thread Thread ID of current task
643 : !> \par History
644 : !> 02.2009 created [Manuel Guidon]
645 : !> \author Manuel Guidon
646 : !> \note
647 : !> This routine calculates the pair-distribution radius of a product
648 : !> Gaussian and uses the powell optimiztion routines in order to fit
649 : !> the results in the following form:
650 : !>
651 : !> (ab| = (ab(Rab) = c2*Rab^2 + c0
652 : !>
653 : ! **************************************************************************************************
654 :
655 1490 : SUBROUTINE calc_pair_dist_radii(qs_env, basis_parameter, &
656 : radii_pgf, max_set, max_pgf, eps_schwarz, &
657 : n_threads, i_thread)
658 :
659 : TYPE(qs_environment_type), POINTER :: qs_env
660 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
661 : TYPE(hfx_screen_coeff_type), &
662 : DIMENSION(:, :, :, :, :, :), POINTER :: radii_pgf
663 : INTEGER, INTENT(IN) :: max_set, max_pgf
664 : REAL(dp) :: eps_schwarz
665 : INTEGER, INTENT(IN) :: n_threads, i_thread
666 :
667 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_pair_dist_radii'
668 :
669 : INTEGER :: handle, i, ikind, ipgf, iset, jkind, &
670 : jpgf, jset, la, lb, ncoa, ncob, nkind, &
671 : nseta, nsetb, sgfa, sgfb
672 1490 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
673 1490 : npgfb, nsgfa, nsgfb
674 1490 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
675 : REAL(dp) :: cutoff, ff, max_contraction_a, max_contraction_b, prefactor, R1, R_max, ra(3), &
676 : rab(3), rab2, radius, rap(3), rb(3), rp(3), x(2), zetp
677 1490 : REAL(dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
678 : REAL(dp), SAVE :: DATA(2, 0:100)
679 : TYPE(gto_basis_set_type), POINTER :: orb_basis
680 1490 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
681 : TYPE(qs_kind_type), POINTER :: qs_kind
682 :
683 2980 : !$OMP MASTER
684 1490 : CALL timeset(routineN, handle)
685 : !$OMP END MASTER
686 : CALL get_qs_env(qs_env=qs_env, &
687 1490 : qs_kind_set=qs_kind_set)
688 :
689 1490 : nkind = SIZE(qs_kind_set, 1)
690 1490 : ra = 0.0_dp
691 1490 : rb = 0.0_dp
692 1490 : !$OMP MASTER
693 1387872 : ALLOCATE (radii_pgf(max_pgf, max_pgf, max_set, max_set, nkind, nkind))
694 4184 : DO ikind = 1, nkind
695 10366 : DO jkind = 1, nkind
696 26490 : DO iset = 1, max_set
697 94686 : DO jset = 1, max_set
698 313240 : DO ipgf = 1, max_pgf
699 1346482 : DO jpgf = 1, max_pgf
700 3377304 : radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(:) = 0.0_dp
701 : END DO
702 : END DO
703 : END DO
704 : END DO
705 : END DO
706 : END DO
707 :
708 1490 : DATA = 0.0_dp
709 : !$OMP END MASTER
710 1490 : !$OMP BARRIER
711 :
712 4184 : DO ikind = 1, nkind
713 2694 : NULLIFY (qs_kind, orb_basis)
714 2694 : qs_kind => qs_kind_set(ikind)
715 2694 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
716 2694 : NULLIFY (la_max, la_min, npgfa, zeta)
717 :
718 2694 : la_max => basis_parameter(ikind)%lmax
719 2694 : la_min => basis_parameter(ikind)%lmin
720 2694 : npgfa => basis_parameter(ikind)%npgf
721 2694 : nseta = basis_parameter(ikind)%nset
722 2694 : zeta => basis_parameter(ikind)%zet
723 2694 : first_sgfa => basis_parameter(ikind)%first_sgf
724 2694 : sphi_a => basis_parameter(ikind)%sphi
725 2694 : nsgfa => basis_parameter(ikind)%nsgf
726 2694 : rpgfa => basis_parameter(ikind)%pgf_radius
727 :
728 10366 : DO jkind = 1, nkind
729 6182 : NULLIFY (qs_kind, orb_basis)
730 6182 : qs_kind => qs_kind_set(jkind)
731 6182 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
732 6182 : NULLIFY (lb_max, lb_min, npgfb, zetb)
733 :
734 6182 : lb_max => basis_parameter(jkind)%lmax
735 6182 : lb_min => basis_parameter(jkind)%lmin
736 6182 : npgfb => basis_parameter(jkind)%npgf
737 6182 : nsetb = basis_parameter(jkind)%nset
738 6182 : zetb => basis_parameter(jkind)%zet
739 6182 : first_sgfb => basis_parameter(jkind)%first_sgf
740 6182 : sphi_b => basis_parameter(jkind)%sphi
741 6182 : nsgfb => basis_parameter(jkind)%nsgf
742 6182 : rpgfb => basis_parameter(jkind)%pgf_radius
743 :
744 24180 : DO iset = 1, nseta
745 15304 : ncoa = npgfa(iset)*ncoset(la_max(iset))
746 15304 : sgfa = first_sgfa(1, iset)
747 578060 : max_contraction_a = MAXVAL([(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)])
748 71518 : DO jset = 1, nsetb
749 50032 : ncob = npgfb(jset)*ncoset(lb_max(jset))
750 50032 : sgfb = first_sgfb(1, jset)
751 1301730 : max_contraction_b = MAXVAL([(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)])
752 150940 : DO ipgf = 1, npgfa(iset)
753 311334 : DO jpgf = 1, npgfb(jset)
754 175698 : radius = rpgfa(ipgf, iset) + rpgfb(jpgf, jset)
755 175698 : DO i = i_thread, 100, n_threads
756 17745498 : rb(1) = 0.0_dp + 0.01_dp*radius*i
757 17745498 : R_max = 0.0_dp
758 41060338 : DO la = la_min(iset), la_max(iset)
759 73187832 : DO lb = lb_min(jset), lb_max(jset)
760 32127494 : zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
761 32127494 : ff = zetb(jpgf, jset)/zetp
762 32127494 : rab = 0.0_dp
763 32127494 : rab(1) = rb(1)
764 32127494 : rab2 = rb(1)**2
765 32127494 : prefactor = EXP(-zeta(ipgf, iset)*ff*rab2)
766 128509976 : rap(:) = ff*rab(:)
767 128509976 : rp(:) = ra(:) + rap(:)
768 128509976 : rb(:) = ra(:) + rab(:)
769 32127494 : cutoff = 1.0_dp
770 : R1 = exp_radius_very_extended( &
771 : la, la, lb, lb, ra=ra, rb=rb, rp=rp, &
772 32127494 : zetp=zetp, eps=eps_schwarz, prefactor=prefactor, cutoff=cutoff, epsabs=1.0E-12_dp)
773 55442334 : R_max = MAX(R_max, R1)
774 : END DO
775 : END DO
776 17745498 : DATA(1, i) = rb(1)
777 17745498 : DATA(2, i) = R_max
778 : END DO
779 : ! the radius can not be negative, we take that into account in the code as well by using a MAX
780 : ! the functional form we use for fitting does not seem particularly accurate
781 175698 : !$OMP BARRIER
782 175698 : !$OMP MASTER
783 175698 : CALL optimize_it(DATA, x, 0.0_dp)
784 527094 : radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x = x
785 : !$OMP END MASTER
786 261302 : !$OMP BARRIER
787 : END DO !jpgf
788 : END DO !ipgf
789 : END DO
790 : END DO
791 : END DO
792 : END DO
793 1490 : !$OMP MASTER
794 1490 : CALL timestop(handle)
795 : !$OMP END MASTER
796 1490 : END SUBROUTINE calc_pair_dist_radii
797 :
798 : !
799 : !
800 : ! little driver routine for the powell minimizer
801 : ! data is the data to fit, x is of the form (x(1)*DATA(1)**2+x(2))
802 : ! only values of DATA(2) larger than fmin are taken into account
803 : ! it constructs an approximate upper bound of the fitted function
804 : !
805 : !
806 : ! **************************************************************************************************
807 : !> \brief ...
808 : !> \param DATA ...
809 : !> \param x ...
810 : !> \param fmin ...
811 : ! **************************************************************************************************
812 407610 : SUBROUTINE optimize_it(DATA, x, fmin)
813 :
814 : REAL(KIND=dp), INTENT(IN) :: DATA(2, 0:100)
815 : REAL(KIND=dp), INTENT(OUT) :: x(2)
816 : REAL(KIND=dp), INTENT(IN) :: fmin
817 :
818 : INTEGER :: i, k
819 : REAL(KIND=dp) :: f, large_weight, small_weight, weight
820 : TYPE(opt_state_type) :: opt_state
821 :
822 : ! initial values
823 :
824 407610 : x(1) = 0.0_dp
825 407610 : x(2) = 0.0_dp
826 :
827 : ! we go in two steps, first we do the symmetric weight to get a good, unique initial guess
828 : ! we restart afterwards for the assym.
829 : ! the assym function appears to have several local minima, depending on the data to fit
830 : ! the loop over k can make the switch gradual, but there is not much need, seemingly
831 815220 : DO k = 0, 4, 4
832 :
833 815220 : small_weight = 1.0_dp
834 815220 : large_weight = small_weight*(10.0_dp**k)
835 :
836 : ! init opt run
837 815220 : opt_state%state = 0
838 815220 : opt_state%nvar = 2
839 815220 : opt_state%iprint = 3
840 815220 : opt_state%unit = default_output_unit
841 815220 : opt_state%maxfun = 100000
842 815220 : opt_state%rhobeg = 0.1_dp
843 815220 : opt_state%rhoend = 0.000001_dp
844 :
845 63207150 : DO
846 :
847 : ! compute function
848 64022370 : IF (opt_state%state == 2) THEN
849 61576710 : opt_state%f = 0.0_dp
850 6280824420 : DO i = 0, 100
851 6219247710 : f = x(1)*DATA(1, i)**2 + x(2)
852 6219247710 : IF (f > DATA(2, i)) THEN
853 : weight = small_weight
854 : ELSE
855 1317426764 : weight = large_weight
856 : END IF
857 6280824420 : IF (DATA(2, i) > fmin) opt_state%f = opt_state%f + weight*(f - DATA(2, i))**2
858 : END DO
859 : END IF
860 :
861 64022370 : IF (opt_state%state == -1) EXIT
862 63207150 : CALL powell_optimize(opt_state%nvar, x, opt_state)
863 : END DO
864 :
865 : ! dealloc mem
866 815220 : opt_state%state = 8
867 1222830 : CALL powell_optimize(opt_state%nvar, x, opt_state)
868 :
869 : END DO
870 :
871 407610 : END SUBROUTINE optimize_it
872 :
873 : ! **************************************************************************************************
874 : !> \brief Given a 2d index pair, this function returns a 1d index pair for
875 : !> a symmetric upper triangle NxN matrix
876 : !> The compiler should inline this function, therefore it appears in
877 : !> several modules
878 : !> \param i 2d index
879 : !> \param j 2d index
880 : !> \param N matrix size
881 : !> \return ...
882 : !> \par History
883 : !> 03.2009 created [Manuel Guidon]
884 : !> \author Manuel Guidon
885 : ! **************************************************************************************************
886 97928 : PURE FUNCTION get_1D_idx(i, j, N)
887 : INTEGER, INTENT(IN) :: i, j
888 : INTEGER(int_8), INTENT(IN) :: N
889 : INTEGER(int_8) :: get_1D_idx
890 :
891 : INTEGER(int_8) :: min_ij
892 :
893 97928 : min_ij = MIN(i, j)
894 97928 : get_1D_idx = min_ij*N + MAX(i, j) - (min_ij - 1)*min_ij/2 - N
895 :
896 97928 : END FUNCTION get_1D_idx
897 :
898 : END MODULE hfx_screening_methods
|
