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 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 hfx_libint_interface, ONLY: evaluate_eri_screen
18 : USE hfx_types, ONLY: hfx_basis_type,&
19 : hfx_p_kind,&
20 : hfx_potential_type,&
21 : hfx_screen_coeff_type,&
22 : log_zero,&
23 : powell_min_log
24 : USE kinds, ONLY: dp,&
25 : int_8
26 : USE libint_wrapper, ONLY: cp_libint_t
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 lib ...
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 contains info for libint
66 : !> \param tmp_R_1 pgf_based screening coefficients
67 : !> \param rab2 Squared Distance of centers ab
68 : !> \param p_work ...
69 : !> \par History
70 : !> 03.2007 created [Manuel Guidon]
71 : !> 02.2009 refactored [Manuel Guidon]
72 : !> \author Manuel Guidon
73 : ! **************************************************************************************************
74 8491272 : SUBROUTINE screen4(lib, ra, rb, zeta, zetb, &
75 : la_min, la_max, lb_min, lb_max, &
76 : npgfa, npgfb, &
77 : max_val, potential_parameter, tmp_R_1, &
78 : rab2, p_work)
79 :
80 : TYPE(cp_libint_t) :: lib
81 : REAL(dp), INTENT(IN) :: ra(3), rb(3)
82 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, zetb
83 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, npgfa, &
84 : npgfb
85 : REAL(dp), INTENT(INOUT) :: max_val
86 : TYPE(hfx_potential_type) :: potential_parameter
87 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
88 : POINTER :: tmp_R_1
89 : REAL(dp) :: rab2
90 : REAL(dp), DIMENSION(:), POINTER :: p_work
91 :
92 : INTEGER :: ipgf, jpgf, la, lb
93 : REAL(dp) :: max_val_temp, R1
94 :
95 8491272 : max_val_temp = max_val
96 22994872 : DO ipgf = 1, npgfa
97 52865016 : DO jpgf = 1, npgfb
98 29870144 : R1 = MAX(0.0_dp, tmp_R_1(jpgf, ipgf)%x(1)*rab2 + tmp_R_1(jpgf, ipgf)%x(2))
99 83896256 : DO la = la_min, la_max
100 124202932 : DO lb = lb_min, lb_max
101 : !Build primitives
102 : CALL evaluate_eri_screen(lib, ra, rb, ra, rb, &
103 : zeta(ipgf), zetb(jpgf), zeta(ipgf), zetb(jpgf), &
104 : la, lb, la, lb, &
105 54810276 : max_val_temp, potential_parameter, R1, R1, p_work)
106 :
107 94332788 : max_val = MAX(max_val, max_val_temp)
108 : END DO !lb
109 : END DO !la
110 : END DO !jpgf
111 : END DO !ipgf
112 8491272 : END SUBROUTINE screen4
113 :
114 : ! **************************************************************************************************
115 : !> \brief updates the maximum of the density matrix in compressed form for screening purposes
116 : !> \param pmax_set buffer to store matrix
117 : !> \param map_atom_to_kind_atom ...
118 : !> \param set_offset ...
119 : !> \param atomic_block_offset ...
120 : !> \param pmax_atom ...
121 : !> \param full_density_alpha ...
122 : !> \param full_density_beta ...
123 : !> \param natom ...
124 : !> \param kind_of ...
125 : !> \param basis_parameter ...
126 : !> \param nkind ...
127 : !> \param is_assoc_atomic_block ...
128 : !> \par History
129 : !> 09.2007 created [Manuel Guidon]
130 : !> \author Manuel Guidon
131 : !> \note
132 : !> - updates for each pair of shells the maximum absolute value of p
133 : ! **************************************************************************************************
134 1412 : 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 706 : nkind, is_assoc_atomic_block)
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 :
151 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_pmax_mat'
152 :
153 : INTEGER :: act_atomic_block_offset, act_set_offset, handle, i, img, jatom, jkind, jset, &
154 : katom, kind_kind_idx, kkind, kset, mb, mc, nimg, nsetb, nsetc
155 706 : INTEGER, DIMENSION(:), POINTER :: nsgfb, nsgfc
156 : REAL(dp) :: pmax_tmp
157 :
158 706 : CALL timeset(routineN, handle)
159 :
160 2728 : DO i = 1, SIZE(pmax_set)
161 70266 : pmax_set(i)%p_kind = 0.0_dp
162 : END DO
163 :
164 10630 : pmax_atom = log_zero
165 :
166 706 : nimg = SIZE(full_density_alpha, 2)
167 :
168 2934 : DO jatom = 1, natom
169 2228 : jkind = kind_of(jatom)
170 2228 : nsetb = basis_parameter(jkind)%nset
171 2228 : nsgfb => basis_parameter(jkind)%nsgf
172 :
173 10630 : DO katom = 1, natom
174 7696 : IF (.NOT. is_assoc_atomic_block(jatom, katom) >= 1) CYCLE
175 7648 : kkind = kind_of(katom)
176 7648 : IF (kkind < jkind) CYCLE
177 6126 : nsetc = basis_parameter(kkind)%nset
178 6126 : nsgfc => basis_parameter(kkind)%nsgf
179 6126 : act_atomic_block_offset = atomic_block_offset(katom, jatom)
180 23796 : DO jset = 1, nsetb
181 66528 : DO kset = 1, nsetc
182 58832 : IF (katom >= jatom) THEN
183 37262 : pmax_tmp = 0.0_dp
184 37262 : act_set_offset = set_offset(kset, jset, kkind, jkind)
185 74524 : DO img = 1, nimg
186 37262 : i = act_set_offset + act_atomic_block_offset - 1
187 144200 : DO mc = 1, nsgfc(kset)
188 298624 : DO mb = 1, nsgfb(jset)
189 191686 : pmax_tmp = MAX(pmax_tmp, ABS(full_density_alpha(i, img)))
190 191686 : IF (ASSOCIATED(full_density_beta)) THEN
191 40636 : pmax_tmp = MAX(pmax_tmp, ABS(full_density_beta(i, img)))
192 : END IF
193 261362 : i = i + 1
194 : END DO
195 : END DO
196 : END DO
197 37262 : IF (pmax_tmp == 0.0_dp) THEN
198 : pmax_tmp = log_zero
199 : ELSE
200 37036 : pmax_tmp = LOG10(pmax_tmp)
201 : END IF
202 37262 : kind_kind_idx = INT(get_1D_idx(jkind, kkind, INT(nkind, int_8)))
203 : pmax_set(kind_kind_idx)%p_kind(jset, &
204 : kset, &
205 : map_atom_to_kind_atom(jatom), &
206 37262 : map_atom_to_kind_atom(katom)) = pmax_tmp
207 : ELSE
208 6128 : pmax_tmp = 0.0_dp
209 6128 : act_set_offset = set_offset(jset, kset, jkind, kkind)
210 12256 : DO img = 1, nimg
211 23698 : DO mc = 1, nsgfc(kset)
212 11442 : i = act_set_offset + act_atomic_block_offset - 1 + mc - 1
213 43572 : DO mb = 1, nsgfb(jset)
214 26002 : pmax_tmp = MAX(pmax_tmp, ABS(full_density_alpha(i, img)))
215 26002 : IF (ASSOCIATED(full_density_beta)) THEN
216 2842 : pmax_tmp = MAX(pmax_tmp, ABS(full_density_beta(i, img)))
217 : END IF
218 37444 : i = i + nsgfc(kset)
219 : END DO
220 : END DO
221 : END DO
222 6128 : IF (pmax_tmp == 0.0_dp) THEN
223 : pmax_tmp = log_zero
224 : ELSE
225 6102 : pmax_tmp = LOG10(pmax_tmp)
226 : END IF
227 6128 : kind_kind_idx = INT(get_1D_idx(jkind, kkind, INT(nkind, int_8)))
228 : pmax_set(kind_kind_idx)%p_kind(jset, &
229 : kset, &
230 : map_atom_to_kind_atom(jatom), &
231 6128 : map_atom_to_kind_atom(katom)) = pmax_tmp
232 : END IF
233 : END DO
234 : END DO
235 : END DO
236 : END DO
237 :
238 2934 : DO jatom = 1, natom
239 2228 : jkind = kind_of(jatom)
240 2228 : nsetb = basis_parameter(jkind)%nset
241 10630 : DO katom = 1, natom
242 7696 : IF (.NOT. is_assoc_atomic_block(jatom, katom) >= 1) CYCLE
243 7648 : kkind = kind_of(katom)
244 7648 : IF (kkind < jkind) CYCLE
245 6126 : nsetc = basis_parameter(kkind)%nset
246 6126 : pmax_tmp = log_zero
247 21568 : DO jset = 1, nsetb
248 64958 : DO kset = 1, nsetc
249 43390 : kind_kind_idx = INT(get_1D_idx(jkind, kkind, INT(nkind, int_8)))
250 : pmax_tmp = MAX(pmax_set(kind_kind_idx)%p_kind(jset, &
251 : kset, &
252 : map_atom_to_kind_atom(jatom), &
253 58832 : map_atom_to_kind_atom(katom)), pmax_tmp)
254 : END DO
255 : END DO
256 6126 : pmax_atom(jatom, katom) = pmax_tmp
257 9924 : pmax_atom(katom, jatom) = pmax_tmp
258 : END DO
259 : END DO
260 :
261 706 : CALL timestop(handle)
262 :
263 706 : END SUBROUTINE update_pmax_mat
264 :
265 : ! **************************************************************************************************
266 : !> \brief calculates screening functions for schwarz screening
267 : !> \param qs_env qs_env
268 : !> \param basis_parameter ...
269 : !> \param lib structure to libint
270 : !> \param potential_parameter contains infos on potential
271 : !> \param coeffs_set set based coefficients
272 : !> \param coeffs_kind kind based coefficients
273 : !> \param coeffs_pgf pgf based coefficients
274 : !> \param radii_pgf coefficients for long-range screening
275 : !> \param max_set Maximum Number of basis set sets in the system
276 : !> \param max_pgf Maximum Number of basis set pgfs in the system
277 : !> \param n_threads ...
278 : !> \param i_thread Thread ID of current task
279 : !> \param p_work ...
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 1216 : SUBROUTINE calc_screening_functions(qs_env, basis_parameter, lib, potential_parameter, &
301 : coeffs_set, coeffs_kind, coeffs_pgf, radii_pgf, &
302 : max_set, max_pgf, n_threads, i_thread, &
303 : p_work)
304 : TYPE(qs_environment_type), POINTER :: qs_env
305 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
306 : TYPE(cp_libint_t) :: lib
307 : TYPE(hfx_potential_type) :: potential_parameter
308 : TYPE(hfx_screen_coeff_type), &
309 : DIMENSION(:, :, :, :), POINTER :: coeffs_set
310 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
311 : POINTER :: coeffs_kind
312 : TYPE(hfx_screen_coeff_type), &
313 : DIMENSION(:, :, :, :, :, :), POINTER :: coeffs_pgf, radii_pgf
314 : INTEGER, INTENT(IN) :: max_set, max_pgf, n_threads, i_thread
315 : REAL(dp), DIMENSION(:), POINTER :: p_work
316 :
317 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_screening_functions'
318 :
319 : INTEGER :: handle, i, ikind, ipgf, iset, jkind, &
320 : jpgf, jset, la, lb, ncoa, ncob, nkind, &
321 : nseta, nsetb, sgfa, sgfb
322 1216 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
323 1216 : npgfb, nsgfa, nsgfb
324 1216 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
325 : REAL(dp) :: kind_radius_a, kind_radius_b, max_contraction_a, max_contraction_b, max_val, &
326 : max_val_temp, R1, ra(3), radius, rb(3), x(2)
327 1216 : REAL(dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
328 1216 : REAL(dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
329 : REAL(dp), SAVE :: DATA(2, 0:100)
330 : TYPE(gto_basis_set_type), POINTER :: orb_basis
331 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
332 1216 : POINTER :: tmp_R_1
333 1216 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
334 : TYPE(qs_kind_type), POINTER :: qs_kind
335 :
336 2432 : !$OMP MASTER
337 1216 : CALL timeset(routineN, handle)
338 : !$OMP END MASTER
339 :
340 : CALL get_qs_env(qs_env=qs_env, &
341 1216 : qs_kind_set=qs_kind_set)
342 :
343 1216 : nkind = SIZE(qs_kind_set, 1)
344 :
345 1216 : !$OMP MASTER
346 1190272 : ALLOCATE (coeffs_pgf(max_pgf, max_pgf, max_set, max_set, nkind, nkind))
347 :
348 3440 : DO ikind = 1, nkind
349 8740 : DO jkind = 1, nkind
350 22274 : DO iset = 1, max_set
351 79648 : DO jset = 1, max_set
352 264988 : DO ipgf = 1, max_pgf
353 1155838 : DO jpgf = 1, max_pgf
354 2907440 : coeffs_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(:) = 0.0_dp
355 : END DO
356 : END DO
357 : END DO
358 : END DO
359 : END DO
360 : END DO
361 :
362 : !$OMP END MASTER
363 1216 : !$OMP BARRIER
364 1216 : ra = 0.0_dp
365 1216 : rb = 0.0_dp
366 3440 : DO ikind = 1, nkind
367 2224 : NULLIFY (qs_kind, orb_basis)
368 2224 : qs_kind => qs_kind_set(ikind)
369 2224 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
370 2224 : NULLIFY (la_max, la_min, npgfa, zeta)
371 :
372 2224 : la_max => basis_parameter(ikind)%lmax
373 2224 : la_min => basis_parameter(ikind)%lmin
374 2224 : npgfa => basis_parameter(ikind)%npgf
375 2224 : nseta = basis_parameter(ikind)%nset
376 2224 : zeta => basis_parameter(ikind)%zet
377 2224 : set_radius_a => basis_parameter(ikind)%set_radius
378 2224 : first_sgfa => basis_parameter(ikind)%first_sgf
379 2224 : sphi_a => basis_parameter(ikind)%sphi
380 2224 : nsgfa => basis_parameter(ikind)%nsgf
381 2224 : rpgfa => basis_parameter(ikind)%pgf_radius
382 :
383 8740 : DO jkind = 1, nkind
384 5300 : NULLIFY (qs_kind, orb_basis)
385 5300 : qs_kind => qs_kind_set(jkind)
386 5300 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
387 5300 : NULLIFY (lb_max, lb_min, npgfb, zetb)
388 :
389 5300 : lb_max => basis_parameter(jkind)%lmax
390 5300 : lb_min => basis_parameter(jkind)%lmin
391 5300 : npgfb => basis_parameter(jkind)%npgf
392 5300 : nsetb = basis_parameter(jkind)%nset
393 5300 : zetb => basis_parameter(jkind)%zet
394 5300 : set_radius_b => basis_parameter(jkind)%set_radius
395 5300 : first_sgfb => basis_parameter(jkind)%first_sgf
396 5300 : sphi_b => basis_parameter(jkind)%sphi
397 5300 : nsgfb => basis_parameter(jkind)%nsgf
398 5300 : rpgfb => basis_parameter(jkind)%pgf_radius
399 :
400 20386 : DO iset = 1, nseta
401 12862 : ncoa = npgfa(iset)*ncoset(la_max(iset))
402 12862 : sgfa = first_sgfa(1, iset)
403 518186 : max_contraction_a = MAXVAL((/(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/))
404 60198 : DO jset = 1, nsetb
405 42036 : ncob = npgfb(jset)*ncoset(lb_max(jset))
406 42036 : sgfb = first_sgfb(1, jset)
407 1160944 : max_contraction_b = MAXVAL((/(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/))
408 42036 : radius = set_radius_a(iset) + set_radius_b(jset)
409 126698 : DO ipgf = 1, npgfa(iset)
410 261708 : DO jpgf = 1, npgfb(jset)
411 147872 : radius = rpgfa(ipgf, iset) + rpgfb(jpgf, jset)
412 15082944 : DO i = i_thread, 100, n_threads
413 14935072 : rb(1) = 0.0_dp + REAL(i, dp)*0.01_dp*radius
414 14935072 : max_val = 0.0_dp
415 : R1 = MAX(0.0_dp, radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(1)*rb(1)**2 + &
416 14935072 : radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(2))
417 34696328 : DO la = la_min(iset), la_max(iset)
418 62101466 : DO lb = lb_min(jset), lb_max(jset)
419 : !Build primitives
420 27405138 : max_val_temp = 0.0_dp
421 : CALL evaluate_eri_screen(lib, ra, rb, ra, rb, &
422 : zeta(ipgf, iset), zetb(jpgf, jset), zeta(ipgf, iset), zetb(jpgf, jset), &
423 : la, lb, la, lb, &
424 27405138 : max_val_temp, potential_parameter, R1, R1, p_work)
425 47166394 : max_val = MAX(max_val, max_val_temp)
426 : END DO !lb
427 : END DO !la
428 14935072 : max_val = SQRT(max_val)
429 14935072 : max_val = max_val*max_contraction_a*max_contraction_b
430 14935072 : DATA(1, i) = rb(1)
431 15082944 : IF (max_val == 0.0_dp) THEN
432 0 : DATA(2, i) = powell_min_log
433 : ELSE
434 14935072 : DATA(2, i) = LOG10((max_val))
435 : END IF
436 : END DO
437 147872 : !$OMP BARRIER
438 147872 : !$OMP MASTER
439 147872 : CALL optimize_it(DATA, x, powell_min_log)
440 443616 : coeffs_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x = x
441 : !$OMP END MASTER
442 219672 : !$OMP BARRIER
443 :
444 : END DO
445 : END DO
446 : END DO
447 : END DO
448 : END DO
449 : END DO
450 :
451 1216 : !$OMP MASTER
452 91600 : ALLOCATE (coeffs_set(max_set, max_set, nkind, nkind))
453 :
454 3440 : DO ikind = 1, nkind
455 8740 : DO jkind = 1, nkind
456 22274 : DO iset = 1, max_set
457 79648 : DO jset = 1, max_set
458 193544 : coeffs_set(jset, iset, jkind, ikind)%x(:) = 0.0_dp
459 : END DO
460 : END DO
461 : END DO
462 : END DO
463 : !$OMP END MASTER
464 1216 : !$OMP BARRIER
465 1216 : ra = 0.0_dp
466 1216 : rb = 0.0_dp
467 3440 : DO ikind = 1, nkind
468 2224 : NULLIFY (qs_kind, orb_basis)
469 2224 : qs_kind => qs_kind_set(ikind)
470 2224 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
471 2224 : NULLIFY (la_max, la_min, npgfa, zeta)
472 : ! CALL get_gto_basis_set(gto_basis_set=orb_basis,&
473 : ! lmax=la_max,&
474 : ! lmin=la_min,&
475 : ! npgf=npgfa,&
476 : ! nset=nseta,&
477 : ! zet=zeta,&
478 : ! set_radius=set_radius_a,&
479 : ! first_sgf=first_sgfa,&
480 : ! sphi=sphi_a,&
481 : ! nsgf_set=nsgfa)
482 2224 : la_max => basis_parameter(ikind)%lmax
483 2224 : la_min => basis_parameter(ikind)%lmin
484 2224 : npgfa => basis_parameter(ikind)%npgf
485 2224 : nseta = basis_parameter(ikind)%nset
486 2224 : zeta => basis_parameter(ikind)%zet
487 2224 : set_radius_a => basis_parameter(ikind)%set_radius
488 2224 : first_sgfa => basis_parameter(ikind)%first_sgf
489 2224 : sphi_a => basis_parameter(ikind)%sphi
490 2224 : nsgfa => basis_parameter(ikind)%nsgf
491 :
492 8740 : DO jkind = 1, nkind
493 5300 : NULLIFY (qs_kind, orb_basis)
494 5300 : qs_kind => qs_kind_set(jkind)
495 5300 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
496 5300 : NULLIFY (lb_max, lb_min, npgfb, zetb)
497 :
498 5300 : lb_max => basis_parameter(jkind)%lmax
499 5300 : lb_min => basis_parameter(jkind)%lmin
500 5300 : npgfb => basis_parameter(jkind)%npgf
501 5300 : nsetb = basis_parameter(jkind)%nset
502 5300 : zetb => basis_parameter(jkind)%zet
503 5300 : set_radius_b => basis_parameter(jkind)%set_radius
504 5300 : first_sgfb => basis_parameter(jkind)%first_sgf
505 5300 : sphi_b => basis_parameter(jkind)%sphi
506 5300 : nsgfb => basis_parameter(jkind)%nsgf
507 :
508 20386 : DO iset = 1, nseta
509 12862 : ncoa = npgfa(iset)*ncoset(la_max(iset))
510 12862 : sgfa = first_sgfa(1, iset)
511 518186 : max_contraction_a = MAXVAL((/(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/))
512 60198 : DO jset = 1, nsetb
513 42036 : ncob = npgfb(jset)*ncoset(lb_max(jset))
514 42036 : sgfb = first_sgfb(1, jset)
515 1160944 : max_contraction_b = MAXVAL((/(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/))
516 42036 : radius = set_radius_a(iset) + set_radius_b(jset)
517 42036 : tmp_R_1 => radii_pgf(:, :, jset, iset, jkind, ikind)
518 4287672 : DO i = i_thread, 100, n_threads
519 4245636 : rb(1) = 0.0_dp + REAL(i, dp)*0.01_dp*radius
520 4245636 : max_val = 0.0_dp
521 : CALL screen4(lib, ra, rb, &
522 : zeta(:, iset), zetb(:, jset), &
523 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
524 : npgfa(iset), npgfb(jset), &
525 4245636 : max_val, potential_parameter, tmp_R_1, rb(1)**2, p_work)
526 4245636 : max_val = SQRT(max_val)
527 4245636 : max_val = max_val*max_contraction_a*max_contraction_b
528 4245636 : DATA(1, i) = rb(1)
529 4287672 : IF (max_val == 0.0_dp) THEN
530 0 : DATA(2, i) = powell_min_log
531 : ELSE
532 4245636 : DATA(2, i) = LOG10((max_val))
533 : END IF
534 : END DO
535 42036 : !$OMP BARRIER
536 42036 : !$OMP MASTER
537 42036 : CALL optimize_it(DATA, x, powell_min_log)
538 126108 : coeffs_set(jset, iset, jkind, ikind)%x = x
539 : !$OMP END MASTER
540 54898 : !$OMP BARRIER
541 : END DO
542 : END DO
543 :
544 : END DO
545 : END DO
546 :
547 : ! ** now kinds
548 1216 : !$OMP MASTER
549 14820 : ALLOCATE (coeffs_kind(nkind, nkind))
550 :
551 3440 : DO ikind = 1, nkind
552 8740 : DO jkind = 1, nkind
553 18124 : coeffs_kind(jkind, ikind)%x(:) = 0.0_dp
554 : END DO
555 : END DO
556 : !$OMP END MASTER
557 1216 : ra = 0.0_dp
558 1216 : rb = 0.0_dp
559 3440 : DO ikind = 1, nkind
560 2224 : NULLIFY (qs_kind, orb_basis)
561 2224 : qs_kind => qs_kind_set(ikind)
562 2224 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
563 2224 : NULLIFY (la_max, la_min, npgfa, zeta)
564 :
565 2224 : la_max => basis_parameter(ikind)%lmax
566 2224 : la_min => basis_parameter(ikind)%lmin
567 2224 : npgfa => basis_parameter(ikind)%npgf
568 2224 : nseta = basis_parameter(ikind)%nset
569 2224 : zeta => basis_parameter(ikind)%zet
570 2224 : kind_radius_a = basis_parameter(ikind)%kind_radius
571 2224 : first_sgfa => basis_parameter(ikind)%first_sgf
572 2224 : sphi_a => basis_parameter(ikind)%sphi
573 2224 : nsgfa => basis_parameter(ikind)%nsgf
574 :
575 8740 : DO jkind = 1, nkind
576 5300 : NULLIFY (qs_kind, orb_basis)
577 5300 : qs_kind => qs_kind_set(jkind)
578 5300 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
579 5300 : NULLIFY (lb_max, lb_min, npgfb, zetb)
580 :
581 5300 : lb_max => basis_parameter(jkind)%lmax
582 5300 : lb_min => basis_parameter(jkind)%lmin
583 5300 : npgfb => basis_parameter(jkind)%npgf
584 5300 : nsetb = basis_parameter(jkind)%nset
585 5300 : zetb => basis_parameter(jkind)%zet
586 5300 : kind_radius_b = basis_parameter(jkind)%kind_radius
587 5300 : first_sgfb => basis_parameter(jkind)%first_sgf
588 5300 : sphi_b => basis_parameter(jkind)%sphi
589 5300 : nsgfb => basis_parameter(jkind)%nsgf
590 :
591 5300 : radius = kind_radius_a + kind_radius_b
592 18162 : DO iset = 1, nseta
593 12862 : ncoa = npgfa(iset)*ncoset(la_max(iset))
594 12862 : sgfa = first_sgfa(1, iset)
595 518186 : max_contraction_a = MAXVAL((/(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/))
596 60198 : DO jset = 1, nsetb
597 42036 : ncob = npgfb(jset)*ncoset(lb_max(jset))
598 42036 : sgfb = first_sgfb(1, jset)
599 1160944 : max_contraction_b = MAXVAL((/(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/))
600 4300534 : DO i = i_thread, 100, n_threads
601 4245636 : rb(1) = 0.0_dp + REAL(i, dp)*0.01_dp*radius
602 4245636 : max_val = 0.0_dp
603 4245636 : tmp_R_1 => radii_pgf(:, :, jset, iset, jkind, ikind)
604 : CALL screen4(lib, ra, rb, &
605 : zeta(:, iset), zetb(:, jset), &
606 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
607 : npgfa(iset), npgfb(jset), &
608 4245636 : max_val, potential_parameter, tmp_R_1, rb(1)**2, p_work)
609 4245636 : DATA(1, i) = rb(1)
610 4245636 : max_val = SQRT(max_val)
611 4245636 : max_val = max_val*max_contraction_a*max_contraction_b
612 4287672 : IF (max_val == 0.0_dp) THEN
613 143036 : DATA(2, i) = MAX(powell_min_log, DATA(2, i))
614 : ELSE
615 4102600 : DATA(2, i) = MAX(LOG10(max_val), DATA(2, i))
616 : END IF
617 : END DO
618 : END DO
619 : END DO
620 5300 : !$OMP BARRIER
621 5300 : !$OMP MASTER
622 5300 : CALL optimize_it(DATA, x, powell_min_log)
623 15900 : coeffs_kind(jkind, ikind)%x = x
624 : !$OMP END MASTER
625 7524 : !$OMP BARRIER
626 : END DO
627 : END DO
628 :
629 1216 : !$OMP MASTER
630 1216 : CALL timestop(handle)
631 : !$OMP END MASTER
632 :
633 1216 : END SUBROUTINE calc_screening_functions
634 :
635 : ! **************************************************************************************************
636 : !> \brief calculates radius functions for longrange screening
637 : !> \param qs_env qs_env
638 : !> \param basis_parameter ...
639 : !> \param radii_pgf pgf based coefficients
640 : !> \param max_set Maximum Number of basis set sets in the system
641 : !> \param max_pgf Maximum Number of basis set pgfs in the system
642 : !> \param eps_schwarz ...
643 : !> \param n_threads ...
644 : !> \param i_thread Thread ID of current task
645 : !> \par History
646 : !> 02.2009 created [Manuel Guidon]
647 : !> \author Manuel Guidon
648 : !> \note
649 : !> This routine calculates the pair-distribution radius of a product
650 : !> Gaussian and uses the powell optimiztion routines in order to fit
651 : !> the results in the following form:
652 : !>
653 : !> (ab| = (ab(Rab) = c2*Rab^2 + c0
654 : !>
655 : ! **************************************************************************************************
656 :
657 1216 : SUBROUTINE calc_pair_dist_radii(qs_env, basis_parameter, &
658 : radii_pgf, max_set, max_pgf, eps_schwarz, &
659 : n_threads, i_thread)
660 :
661 : TYPE(qs_environment_type), POINTER :: qs_env
662 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
663 : TYPE(hfx_screen_coeff_type), &
664 : DIMENSION(:, :, :, :, :, :), POINTER :: radii_pgf
665 : INTEGER, INTENT(IN) :: max_set, max_pgf
666 : REAL(dp) :: eps_schwarz
667 : INTEGER, INTENT(IN) :: n_threads, i_thread
668 :
669 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_pair_dist_radii'
670 :
671 : INTEGER :: handle, i, ikind, ipgf, iset, jkind, &
672 : jpgf, jset, la, lb, ncoa, ncob, nkind, &
673 : nseta, nsetb, sgfa, sgfb
674 1216 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
675 1216 : npgfb, nsgfa, nsgfb
676 1216 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
677 : REAL(dp) :: cutoff, ff, max_contraction_a, max_contraction_b, prefactor, R1, R_max, ra(3), &
678 : rab(3), rab2, radius, rap(3), rb(3), rp(3), x(2), zetp
679 1216 : REAL(dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
680 : REAL(dp), SAVE :: DATA(2, 0:100)
681 : TYPE(gto_basis_set_type), POINTER :: orb_basis
682 1216 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
683 : TYPE(qs_kind_type), POINTER :: qs_kind
684 :
685 2432 : !$OMP MASTER
686 1216 : CALL timeset(routineN, handle)
687 : !$OMP END MASTER
688 : CALL get_qs_env(qs_env=qs_env, &
689 1216 : qs_kind_set=qs_kind_set)
690 :
691 1216 : nkind = SIZE(qs_kind_set, 1)
692 1216 : ra = 0.0_dp
693 1216 : rb = 0.0_dp
694 1216 : !$OMP MASTER
695 1190272 : ALLOCATE (radii_pgf(max_pgf, max_pgf, max_set, max_set, nkind, nkind))
696 3440 : DO ikind = 1, nkind
697 8740 : DO jkind = 1, nkind
698 22274 : DO iset = 1, max_set
699 79648 : DO jset = 1, max_set
700 264988 : DO ipgf = 1, max_pgf
701 1155838 : DO jpgf = 1, max_pgf
702 2907440 : radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(:) = 0.0_dp
703 : END DO
704 : END DO
705 : END DO
706 : END DO
707 : END DO
708 : END DO
709 :
710 1216 : DATA = 0.0_dp
711 : !$OMP END MASTER
712 1216 : !$OMP BARRIER
713 :
714 3440 : DO ikind = 1, nkind
715 2224 : NULLIFY (qs_kind, orb_basis)
716 2224 : qs_kind => qs_kind_set(ikind)
717 2224 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
718 2224 : NULLIFY (la_max, la_min, npgfa, zeta)
719 :
720 2224 : la_max => basis_parameter(ikind)%lmax
721 2224 : la_min => basis_parameter(ikind)%lmin
722 2224 : npgfa => basis_parameter(ikind)%npgf
723 2224 : nseta = basis_parameter(ikind)%nset
724 2224 : zeta => basis_parameter(ikind)%zet
725 2224 : first_sgfa => basis_parameter(ikind)%first_sgf
726 2224 : sphi_a => basis_parameter(ikind)%sphi
727 2224 : nsgfa => basis_parameter(ikind)%nsgf
728 2224 : rpgfa => basis_parameter(ikind)%pgf_radius
729 :
730 8740 : DO jkind = 1, nkind
731 5300 : NULLIFY (qs_kind, orb_basis)
732 5300 : qs_kind => qs_kind_set(jkind)
733 5300 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
734 5300 : NULLIFY (lb_max, lb_min, npgfb, zetb)
735 :
736 5300 : lb_max => basis_parameter(jkind)%lmax
737 5300 : lb_min => basis_parameter(jkind)%lmin
738 5300 : npgfb => basis_parameter(jkind)%npgf
739 5300 : nsetb = basis_parameter(jkind)%nset
740 5300 : zetb => basis_parameter(jkind)%zet
741 5300 : first_sgfb => basis_parameter(jkind)%first_sgf
742 5300 : sphi_b => basis_parameter(jkind)%sphi
743 5300 : nsgfb => basis_parameter(jkind)%nsgf
744 5300 : rpgfb => basis_parameter(jkind)%pgf_radius
745 :
746 20386 : DO iset = 1, nseta
747 12862 : ncoa = npgfa(iset)*ncoset(la_max(iset))
748 12862 : sgfa = first_sgfa(1, iset)
749 518186 : max_contraction_a = MAXVAL((/(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/))
750 60198 : DO jset = 1, nsetb
751 42036 : ncob = npgfb(jset)*ncoset(lb_max(jset))
752 42036 : sgfb = first_sgfb(1, jset)
753 1160944 : max_contraction_b = MAXVAL((/(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/))
754 126698 : DO ipgf = 1, npgfa(iset)
755 261708 : DO jpgf = 1, npgfb(jset)
756 147872 : radius = rpgfa(ipgf, iset) + rpgfb(jpgf, jset)
757 147872 : DO i = i_thread, 100, n_threads
758 14935072 : rb(1) = 0.0_dp + 0.01_dp*radius*i
759 14935072 : R_max = 0.0_dp
760 34696328 : DO la = la_min(iset), la_max(iset)
761 62101466 : DO lb = lb_min(jset), lb_max(jset)
762 27405138 : zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
763 27405138 : ff = zetb(jpgf, jset)/zetp
764 27405138 : rab = 0.0_dp
765 27405138 : rab(1) = rb(1)
766 27405138 : rab2 = rb(1)**2
767 27405138 : prefactor = EXP(-zeta(ipgf, iset)*ff*rab2)
768 109620552 : rap(:) = ff*rab(:)
769 109620552 : rp(:) = ra(:) + rap(:)
770 109620552 : rb(:) = ra(:) + rab(:)
771 27405138 : cutoff = 1.0_dp
772 : R1 = exp_radius_very_extended( &
773 : la, la, lb, lb, ra=ra, rb=rb, rp=rp, &
774 27405138 : zetp=zetp, eps=eps_schwarz, prefactor=prefactor, cutoff=cutoff, epsabs=1.0E-12_dp)
775 47166394 : R_max = MAX(R_max, R1)
776 : END DO
777 : END DO
778 14935072 : DATA(1, i) = rb(1)
779 14935072 : DATA(2, i) = R_max
780 : END DO
781 : ! the radius can not be negative, we take that into account in the code as well by using a MAX
782 : ! the functional form we use for fitting does not seem particularly accurate
783 147872 : !$OMP BARRIER
784 147872 : !$OMP MASTER
785 147872 : CALL optimize_it(DATA, x, 0.0_dp)
786 443616 : radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x = x
787 : !$OMP END MASTER
788 219672 : !$OMP BARRIER
789 : END DO !jpgf
790 : END DO !ipgf
791 : END DO
792 : END DO
793 : END DO
794 : END DO
795 1216 : !$OMP MASTER
796 1216 : CALL timestop(handle)
797 : !$OMP END MASTER
798 1216 : END SUBROUTINE calc_pair_dist_radii
799 :
800 : !
801 : !
802 : ! little driver routine for the powell minimizer
803 : ! data is the data to fit, x is of the form (x(1)*DATA(1)**2+x(2))
804 : ! only values of DATA(2) larger than fmin are taken into account
805 : ! it constructs an approximate upper bound of the fitted function
806 : !
807 : !
808 : ! **************************************************************************************************
809 : !> \brief ...
810 : !> \param DATA ...
811 : !> \param x ...
812 : !> \param fmin ...
813 : ! **************************************************************************************************
814 343080 : SUBROUTINE optimize_it(DATA, x, fmin)
815 :
816 : REAL(KIND=dp), INTENT(IN) :: DATA(2, 0:100)
817 : REAL(KIND=dp), INTENT(OUT) :: x(2)
818 : REAL(KIND=dp), INTENT(IN) :: fmin
819 :
820 : INTEGER :: i, k
821 : REAL(KIND=dp) :: f, large_weight, small_weight, weight
822 : TYPE(opt_state_type) :: opt_state
823 :
824 : ! initial values
825 :
826 343080 : x(1) = 0.0_dp
827 343080 : x(2) = 0.0_dp
828 :
829 : ! we go in two steps, first we do the symmetric weight to get a good, unique initial guess
830 : ! we restart afterwards for the assym.
831 : ! the assym function appears to have several local minima, depending on the data to fit
832 : ! the loop over k can make the switch gradual, but there is not much need, seemingly
833 686160 : DO k = 0, 4, 4
834 :
835 686160 : small_weight = 1.0_dp
836 686160 : large_weight = small_weight*(10.0_dp**k)
837 :
838 : ! init opt run
839 686160 : opt_state%state = 0
840 686160 : opt_state%nvar = 2
841 686160 : opt_state%iprint = 3
842 686160 : opt_state%unit = default_output_unit
843 686160 : opt_state%maxfun = 100000
844 686160 : opt_state%rhobeg = 0.1_dp
845 686160 : opt_state%rhoend = 0.000001_dp
846 :
847 53303754 : DO
848 :
849 : ! compute function
850 53989914 : IF (opt_state%state == 2) THEN
851 51931434 : opt_state%f = 0.0_dp
852 5297006268 : DO i = 0, 100
853 5245074834 : f = x(1)*DATA(1, i)**2 + x(2)
854 5245074834 : IF (f > DATA(2, i)) THEN
855 : weight = small_weight
856 : ELSE
857 1134975616 : weight = large_weight
858 : END IF
859 5297006268 : IF (DATA(2, i) > fmin) opt_state%f = opt_state%f + weight*(f - DATA(2, i))**2
860 : END DO
861 : END IF
862 :
863 53989914 : IF (opt_state%state == -1) EXIT
864 53303754 : CALL powell_optimize(opt_state%nvar, x, opt_state)
865 : END DO
866 :
867 : ! dealloc mem
868 686160 : opt_state%state = 8
869 1029240 : CALL powell_optimize(opt_state%nvar, x, opt_state)
870 :
871 : END DO
872 :
873 343080 : END SUBROUTINE optimize_it
874 :
875 : ! **************************************************************************************************
876 : !> \brief Given a 2d index pair, this function returns a 1d index pair for
877 : !> a symmetric upper triangle NxN matrix
878 : !> The compiler should inline this function, therefore it appears in
879 : !> several modules
880 : !> \param i 2d index
881 : !> \param j 2d index
882 : !> \param N matrix size
883 : !> \return ...
884 : !> \par History
885 : !> 03.2009 created [Manuel Guidon]
886 : !> \author Manuel Guidon
887 : ! **************************************************************************************************
888 86780 : PURE FUNCTION get_1D_idx(i, j, N)
889 : INTEGER, INTENT(IN) :: i, j
890 : INTEGER(int_8), INTENT(IN) :: N
891 : INTEGER(int_8) :: get_1D_idx
892 :
893 : INTEGER(int_8) :: min_ij
894 :
895 86780 : min_ij = MIN(i, j)
896 86780 : get_1D_idx = min_ij*N + MAX(i, j) - (min_ij - 1)*min_ij/2 - N
897 :
898 86780 : END FUNCTION get_1D_idx
899 :
900 : END MODULE hfx_screening_methods
|
