(git:f7b9015)
Loading...
Searching...
No Matches
hfx_screening_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
14MODULE 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
48CONTAINS
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 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 max_val_temp = max_val
94 DO ipgf = 1, npgfa
95 DO jpgf = 1, npgfb
96 r1 = max(0.0_dp, tmp_r_1(jpgf, ipgf)%x(1)*rab2 + tmp_r_1(jpgf, ipgf)%x(2))
97 DO la = la_min, la_max
98 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 max_val_temp, potential_parameter, r1, r1)
104
105 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 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 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 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 INTEGER, DIMENSION(:), POINTER :: nsgfb, nsgfc
157 REAL(dp) :: pmax_tmp
158
159 CALL timeset(routinen, handle)
160
161 DO i = 1, SIZE(pmax_set)
162 pmax_set(i)%p_kind = 0.0_dp
163 END DO
164
165 pmax_atom = log_zero
166
167 nimg = SIZE(full_density_alpha, 2)
168
169 DO jatom = 1, natom
170 jkind = kind_of(jatom)
171 nsetb = basis_parameter(jkind)%nset
172 nsgfb => basis_parameter(jkind)%nsgf
173
174 DO katom = 1, natom
175 IF (.NOT. is_assoc_atomic_block(jatom, katom) >= 1) cycle
176 kkind = kind_of(katom)
177 IF (kkind < jkind) cycle
178 nsetc = basis_parameter(kkind)%nset
179 nsgfc => basis_parameter(kkind)%nsgf
180 act_atomic_block_offset = atomic_block_offset(katom, jatom)
181 DO jset = 1, nsetb
182 DO kset = 1, nsetc
183 IF (katom >= jatom) THEN
184 pmax_tmp = 0.0_dp
185 act_set_offset = set_offset(kset, jset, kkind, jkind)
186 DO img = 1, nimg
187 i = act_set_offset + act_atomic_block_offset - 1
188 DO mc = 1, nsgfc(kset)
189 DO mb = 1, nsgfb(jset)
190 pmax_tmp = max(pmax_tmp, abs(full_density_alpha(i, img)))
191 IF (ASSOCIATED(full_density_beta)) THEN
192 pmax_tmp = max(pmax_tmp, abs(full_density_beta(i, img)))
193 END IF
194 i = i + 1
195 END DO
196 END DO
197 END DO
198 IF (pmax_tmp == 0.0_dp) THEN
199 pmax_tmp = log_zero
200 ELSE
201 pmax_tmp = log10(real(nspins, dp)*pmax_tmp)
202 END IF
203 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 map_atom_to_kind_atom(katom)) = pmax_tmp
208 ELSE
209 pmax_tmp = 0.0_dp
210 act_set_offset = set_offset(jset, kset, jkind, kkind)
211 DO img = 1, nimg
212 DO mc = 1, nsgfc(kset)
213 i = act_set_offset + act_atomic_block_offset - 1 + mc - 1
214 DO mb = 1, nsgfb(jset)
215 pmax_tmp = max(pmax_tmp, abs(full_density_alpha(i, img)))
216 IF (ASSOCIATED(full_density_beta)) THEN
217 pmax_tmp = max(pmax_tmp, abs(full_density_beta(i, img)))
218 END IF
219 i = i + nsgfc(kset)
220 END DO
221 END DO
222 END DO
223 IF (pmax_tmp == 0.0_dp) THEN
224 pmax_tmp = log_zero
225 ELSE
226 pmax_tmp = log10(real(nspins, dp)*pmax_tmp)
227 END IF
228 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 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 DO jatom = 1, natom
240 jkind = kind_of(jatom)
241 nsetb = basis_parameter(jkind)%nset
242 DO katom = 1, natom
243 IF (.NOT. is_assoc_atomic_block(jatom, katom) >= 1) cycle
244 kkind = kind_of(katom)
245 IF (kkind < jkind) cycle
246 nsetc = basis_parameter(kkind)%nset
247 pmax_tmp = log_zero
248 DO jset = 1, nsetb
249 DO kset = 1, nsetc
250 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 map_atom_to_kind_atom(katom)), pmax_tmp)
255 END DO
256 END DO
257 pmax_atom(jatom, katom) = pmax_tmp
258 pmax_atom(katom, jatom) = pmax_tmp
259 END DO
260 END DO
261
262 CALL timestop(handle)
263
264 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 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 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
321 npgfb, nsgfa, nsgfb
322 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 REAL(dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
326 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 POINTER :: tmp_r_1
331 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
332 TYPE(qs_kind_type), POINTER :: qs_kind
333
334!$OMP MASTER
335 CALL timeset(routinen, handle)
336!$OMP END MASTER
337
338 CALL get_qs_env(qs_env=qs_env, &
339 qs_kind_set=qs_kind_set)
340
341 nkind = SIZE(qs_kind_set, 1)
342
343!$OMP MASTER
344 ALLOCATE (coeffs_pgf(max_pgf, max_pgf, max_set, max_set, nkind, nkind))
345
346 DO ikind = 1, nkind
347 DO jkind = 1, nkind
348 DO iset = 1, max_set
349 DO jset = 1, max_set
350 DO ipgf = 1, max_pgf
351 DO jpgf = 1, max_pgf
352 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!$OMP BARRIER
362 ra = 0.0_dp
363 rb = 0.0_dp
364 DO ikind = 1, nkind
365 NULLIFY (qs_kind, orb_basis)
366 qs_kind => qs_kind_set(ikind)
367 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
368 NULLIFY (la_max, la_min, npgfa, zeta)
369
370 la_max => basis_parameter(ikind)%lmax
371 la_min => basis_parameter(ikind)%lmin
372 npgfa => basis_parameter(ikind)%npgf
373 nseta = basis_parameter(ikind)%nset
374 zeta => basis_parameter(ikind)%zet
375 set_radius_a => basis_parameter(ikind)%set_radius
376 first_sgfa => basis_parameter(ikind)%first_sgf
377 sphi_a => basis_parameter(ikind)%sphi
378 nsgfa => basis_parameter(ikind)%nsgf
379 rpgfa => basis_parameter(ikind)%pgf_radius
380
381 DO jkind = 1, nkind
382 NULLIFY (qs_kind, orb_basis)
383 qs_kind => qs_kind_set(jkind)
384 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
385 NULLIFY (lb_max, lb_min, npgfb, zetb)
386
387 lb_max => basis_parameter(jkind)%lmax
388 lb_min => basis_parameter(jkind)%lmin
389 npgfb => basis_parameter(jkind)%npgf
390 nsetb = basis_parameter(jkind)%nset
391 zetb => basis_parameter(jkind)%zet
392 set_radius_b => basis_parameter(jkind)%set_radius
393 first_sgfb => basis_parameter(jkind)%first_sgf
394 sphi_b => basis_parameter(jkind)%sphi
395 nsgfb => basis_parameter(jkind)%nsgf
396 rpgfb => basis_parameter(jkind)%pgf_radius
397
398 DO iset = 1, nseta
399 ncoa = npgfa(iset)*ncoset(la_max(iset))
400 sgfa = first_sgfa(1, iset)
401 max_contraction_a = maxval([(sum(abs(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)])
402 DO jset = 1, nsetb
403 ncob = npgfb(jset)*ncoset(lb_max(jset))
404 sgfb = first_sgfb(1, jset)
405 max_contraction_b = maxval([(sum(abs(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)])
406 radius = set_radius_a(iset) + set_radius_b(jset)
407 DO ipgf = 1, npgfa(iset)
408 DO jpgf = 1, npgfb(jset)
409 radius = rpgfa(ipgf, iset) + rpgfb(jpgf, jset)
410 DO i = i_thread, 100, n_threads
411 rb(1) = 0.0_dp + real(i, dp)*0.01_dp*radius
412 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 radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(2))
415 DO la = la_min(iset), la_max(iset)
416 DO lb = lb_min(jset), lb_max(jset)
417 !Build primitives
418 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 max_val_temp, potential_parameter, r1, r1)
423 max_val = max(max_val, max_val_temp)
424 END DO !lb
425 END DO !la
426 max_val = sqrt(max_val)
427 max_val = max_val*max_contraction_a*max_contraction_b
428 DATA(1, i) = rb(1)
429 IF (max_val == 0.0_dp) THEN
430 DATA(2, i) = powell_min_log
431 ELSE
432 DATA(2, i) = log10((max_val))
433 END IF
434 END DO
435!$OMP BARRIER
436!$OMP MASTER
437 CALL optimize_it(DATA, x, powell_min_log)
438 coeffs_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x = x
439!$OMP END MASTER
440!$OMP BARRIER
441
442 END DO
443 END DO
444 END DO
445 END DO
446 END DO
447 END DO
448
449!$OMP MASTER
450 ALLOCATE (coeffs_set(max_set, max_set, nkind, nkind))
451
452 DO ikind = 1, nkind
453 DO jkind = 1, nkind
454 DO iset = 1, max_set
455 DO jset = 1, max_set
456 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!$OMP BARRIER
463 ra = 0.0_dp
464 rb = 0.0_dp
465 DO ikind = 1, nkind
466 NULLIFY (qs_kind, orb_basis)
467 qs_kind => qs_kind_set(ikind)
468 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
469 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 la_max => basis_parameter(ikind)%lmax
481 la_min => basis_parameter(ikind)%lmin
482 npgfa => basis_parameter(ikind)%npgf
483 nseta = basis_parameter(ikind)%nset
484 zeta => basis_parameter(ikind)%zet
485 set_radius_a => basis_parameter(ikind)%set_radius
486 first_sgfa => basis_parameter(ikind)%first_sgf
487 sphi_a => basis_parameter(ikind)%sphi
488 nsgfa => basis_parameter(ikind)%nsgf
489
490 DO jkind = 1, nkind
491 NULLIFY (qs_kind, orb_basis)
492 qs_kind => qs_kind_set(jkind)
493 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
494 NULLIFY (lb_max, lb_min, npgfb, zetb)
495
496 lb_max => basis_parameter(jkind)%lmax
497 lb_min => basis_parameter(jkind)%lmin
498 npgfb => basis_parameter(jkind)%npgf
499 nsetb = basis_parameter(jkind)%nset
500 zetb => basis_parameter(jkind)%zet
501 set_radius_b => basis_parameter(jkind)%set_radius
502 first_sgfb => basis_parameter(jkind)%first_sgf
503 sphi_b => basis_parameter(jkind)%sphi
504 nsgfb => basis_parameter(jkind)%nsgf
505
506 DO iset = 1, nseta
507 ncoa = npgfa(iset)*ncoset(la_max(iset))
508 sgfa = first_sgfa(1, iset)
509 max_contraction_a = maxval([(sum(abs(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)])
510 DO jset = 1, nsetb
511 ncob = npgfb(jset)*ncoset(lb_max(jset))
512 sgfb = first_sgfb(1, jset)
513 max_contraction_b = maxval([(sum(abs(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)])
514 radius = set_radius_a(iset) + set_radius_b(jset)
515 tmp_r_1 => radii_pgf(:, :, jset, iset, jkind, ikind)
516 DO i = i_thread, 100, n_threads
517 rb(1) = 0.0_dp + real(i, dp)*0.01_dp*radius
518 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 max_val, potential_parameter, tmp_r_1, rb(1)**2)
524 max_val = sqrt(max_val)
525 max_val = max_val*max_contraction_a*max_contraction_b
526 DATA(1, i) = rb(1)
527 IF (max_val == 0.0_dp) THEN
528 DATA(2, i) = powell_min_log
529 ELSE
530 DATA(2, i) = log10((max_val))
531 END IF
532 END DO
533!$OMP BARRIER
534!$OMP MASTER
535 CALL optimize_it(DATA, x, powell_min_log)
536 coeffs_set(jset, iset, jkind, ikind)%x = x
537!$OMP END MASTER
538!$OMP BARRIER
539 END DO
540 END DO
541
542 END DO
543 END DO
544
545 ! ** now kinds
546!$OMP MASTER
547 ALLOCATE (coeffs_kind(nkind, nkind))
548
549 DO ikind = 1, nkind
550 DO jkind = 1, nkind
551 coeffs_kind(jkind, ikind)%x(:) = 0.0_dp
552 END DO
553 END DO
554!$OMP END MASTER
555 ra = 0.0_dp
556 rb = 0.0_dp
557 DO ikind = 1, nkind
558 NULLIFY (qs_kind, orb_basis)
559 qs_kind => qs_kind_set(ikind)
560 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
561 NULLIFY (la_max, la_min, npgfa, zeta)
562
563 la_max => basis_parameter(ikind)%lmax
564 la_min => basis_parameter(ikind)%lmin
565 npgfa => basis_parameter(ikind)%npgf
566 nseta = basis_parameter(ikind)%nset
567 zeta => basis_parameter(ikind)%zet
568 kind_radius_a = basis_parameter(ikind)%kind_radius
569 first_sgfa => basis_parameter(ikind)%first_sgf
570 sphi_a => basis_parameter(ikind)%sphi
571 nsgfa => basis_parameter(ikind)%nsgf
572
573 DO jkind = 1, nkind
574 NULLIFY (qs_kind, orb_basis)
575 qs_kind => qs_kind_set(jkind)
576 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
577 NULLIFY (lb_max, lb_min, npgfb, zetb)
578
579 lb_max => basis_parameter(jkind)%lmax
580 lb_min => basis_parameter(jkind)%lmin
581 npgfb => basis_parameter(jkind)%npgf
582 nsetb = basis_parameter(jkind)%nset
583 zetb => basis_parameter(jkind)%zet
584 kind_radius_b = basis_parameter(jkind)%kind_radius
585 first_sgfb => basis_parameter(jkind)%first_sgf
586 sphi_b => basis_parameter(jkind)%sphi
587 nsgfb => basis_parameter(jkind)%nsgf
588
589 radius = kind_radius_a + kind_radius_b
590 DO iset = 1, nseta
591 ncoa = npgfa(iset)*ncoset(la_max(iset))
592 sgfa = first_sgfa(1, iset)
593 max_contraction_a = maxval([(sum(abs(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)])
594 DO jset = 1, nsetb
595 ncob = npgfb(jset)*ncoset(lb_max(jset))
596 sgfb = first_sgfb(1, jset)
597 max_contraction_b = maxval([(sum(abs(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)])
598 DO i = i_thread, 100, n_threads
599 rb(1) = 0.0_dp + real(i, dp)*0.01_dp*radius
600 max_val = 0.0_dp
601 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 max_val, potential_parameter, tmp_r_1, rb(1)**2)
607 DATA(1, i) = rb(1)
608 max_val = sqrt(max_val)
609 max_val = max_val*max_contraction_a*max_contraction_b
610 IF (max_val == 0.0_dp) THEN
611 DATA(2, i) = max(powell_min_log, DATA(2, i))
612 ELSE
613 DATA(2, i) = max(log10(max_val), DATA(2, i))
614 END IF
615 END DO
616 END DO
617 END DO
618!$OMP BARRIER
619!$OMP MASTER
620 CALL optimize_it(DATA, x, powell_min_log)
621 coeffs_kind(jkind, ikind)%x = x
622!$OMP END MASTER
623!$OMP BARRIER
624 END DO
625 END DO
626
627!$OMP MASTER
628 CALL timestop(handle)
629!$OMP END MASTER
630
631 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 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 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
673 npgfb, nsgfa, nsgfb
674 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 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 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
681 TYPE(qs_kind_type), POINTER :: qs_kind
682
683!$OMP MASTER
684 CALL timeset(routinen, handle)
685!$OMP END MASTER
686 CALL get_qs_env(qs_env=qs_env, &
687 qs_kind_set=qs_kind_set)
688
689 nkind = SIZE(qs_kind_set, 1)
690 ra = 0.0_dp
691 rb = 0.0_dp
692!$OMP MASTER
693 ALLOCATE (radii_pgf(max_pgf, max_pgf, max_set, max_set, nkind, nkind))
694 DO ikind = 1, nkind
695 DO jkind = 1, nkind
696 DO iset = 1, max_set
697 DO jset = 1, max_set
698 DO ipgf = 1, max_pgf
699 DO jpgf = 1, max_pgf
700 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 DATA = 0.0_dp
709!$OMP END MASTER
710!$OMP BARRIER
711
712 DO ikind = 1, nkind
713 NULLIFY (qs_kind, orb_basis)
714 qs_kind => qs_kind_set(ikind)
715 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
716 NULLIFY (la_max, la_min, npgfa, zeta)
717
718 la_max => basis_parameter(ikind)%lmax
719 la_min => basis_parameter(ikind)%lmin
720 npgfa => basis_parameter(ikind)%npgf
721 nseta = basis_parameter(ikind)%nset
722 zeta => basis_parameter(ikind)%zet
723 first_sgfa => basis_parameter(ikind)%first_sgf
724 sphi_a => basis_parameter(ikind)%sphi
725 nsgfa => basis_parameter(ikind)%nsgf
726 rpgfa => basis_parameter(ikind)%pgf_radius
727
728 DO jkind = 1, nkind
729 NULLIFY (qs_kind, orb_basis)
730 qs_kind => qs_kind_set(jkind)
731 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
732 NULLIFY (lb_max, lb_min, npgfb, zetb)
733
734 lb_max => basis_parameter(jkind)%lmax
735 lb_min => basis_parameter(jkind)%lmin
736 npgfb => basis_parameter(jkind)%npgf
737 nsetb = basis_parameter(jkind)%nset
738 zetb => basis_parameter(jkind)%zet
739 first_sgfb => basis_parameter(jkind)%first_sgf
740 sphi_b => basis_parameter(jkind)%sphi
741 nsgfb => basis_parameter(jkind)%nsgf
742 rpgfb => basis_parameter(jkind)%pgf_radius
743
744 DO iset = 1, nseta
745 ncoa = npgfa(iset)*ncoset(la_max(iset))
746 sgfa = first_sgfa(1, iset)
747 max_contraction_a = maxval([(sum(abs(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)])
748 DO jset = 1, nsetb
749 ncob = npgfb(jset)*ncoset(lb_max(jset))
750 sgfb = first_sgfb(1, jset)
751 max_contraction_b = maxval([(sum(abs(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)])
752 DO ipgf = 1, npgfa(iset)
753 DO jpgf = 1, npgfb(jset)
754 radius = rpgfa(ipgf, iset) + rpgfb(jpgf, jset)
755 DO i = i_thread, 100, n_threads
756 rb(1) = 0.0_dp + 0.01_dp*radius*i
757 r_max = 0.0_dp
758 DO la = la_min(iset), la_max(iset)
759 DO lb = lb_min(jset), lb_max(jset)
760 zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
761 ff = zetb(jpgf, jset)/zetp
762 rab = 0.0_dp
763 rab(1) = rb(1)
764 rab2 = rb(1)**2
765 prefactor = exp(-zeta(ipgf, iset)*ff*rab2)
766 rap(:) = ff*rab(:)
767 rp(:) = ra(:) + rap(:)
768 rb(:) = ra(:) + rab(:)
769 cutoff = 1.0_dp
770 r1 = exp_radius_very_extended( &
771 la, la, lb, lb, ra=ra, rb=rb, rp=rp, &
772 zetp=zetp, eps=eps_schwarz, prefactor=prefactor, cutoff=cutoff, epsabs=1.0e-12_dp)
773 r_max = max(r_max, r1)
774 END DO
775 END DO
776 DATA(1, i) = rb(1)
777 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!$OMP BARRIER
782!$OMP MASTER
783 CALL optimize_it(DATA, x, 0.0_dp)
784 radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x = x
785!$OMP END MASTER
786!$OMP BARRIER
787 END DO !jpgf
788 END DO !ipgf
789 END DO
790 END DO
791 END DO
792 END DO
793!$OMP MASTER
794 CALL timestop(handle)
795!$OMP END MASTER
796 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 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 x(1) = 0.0_dp
825 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 DO k = 0, 4, 4
832
833 small_weight = 1.0_dp
834 large_weight = small_weight*(10.0_dp**k)
835
836 ! init opt run
837 opt_state%state = 0
838 opt_state%nvar = 2
839 opt_state%iprint = 3
840 opt_state%unit = default_output_unit
841 opt_state%maxfun = 100000
842 opt_state%rhobeg = 0.1_dp
843 opt_state%rhoend = 0.000001_dp
844
845 DO
846
847 ! compute function
848 IF (opt_state%state == 2) THEN
849 opt_state%f = 0.0_dp
850 DO i = 0, 100
851 f = x(1)*DATA(1, i)**2 + x(2)
852 IF (f > DATA(2, i)) THEN
853 weight = small_weight
854 ELSE
855 weight = large_weight
856 END IF
857 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 IF (opt_state%state == -1) EXIT
862 CALL powell_optimize(opt_state%nvar, x, opt_state)
863 END DO
864
865 ! dealloc mem
866 opt_state%state = 8
867 CALL powell_optimize(opt_state%nvar, x, opt_state)
868
869 END DO
870
871 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 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 min_ij = min(i, j)
894 get_1d_idx = min_ij*n + max(i, j) - (min_ij - 1)*min_ij/2 - n
895
896 END FUNCTION get_1d_idx
897
898END MODULE hfx_screening_methods
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius_very_extended
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Definition ao_util.F:208
basis_set_types
Definition basis_set_types.F:15
four_center_integral_interface
Four-center integral evaluation with an encapsulated Libint backend.
Definition four_center_integral_interface.F:15
four_center_integral_interface::evaluate_eri_screen
subroutine, public evaluate_eri_screen(context, a, b, c, d, zeta_a, zeta_b, zeta_c, zeta_d, n_a, n_b, n_c, n_d, max_val, potential_parameter, r1, r2)
Evaluates max-abs values of electron repulsion integrals for a primitive quartet.
Definition four_center_integral_interface.F:433
hfx_screening_methods
Several screening methods used in HFX calcualtions.
Definition hfx_screening_methods.F:14
hfx_screening_methods::calc_screening_functions
subroutine, public calc_screening_functions(qs_env, basis_parameter, context, potential_parameter, coeffs_set, coeffs_kind, coeffs_pgf, radii_pgf, max_set, max_pgf, n_threads, i_thread)
calculates screening functions for schwarz screening
Definition hfx_screening_methods.F:303
hfx_screening_methods::update_pmax_mat
subroutine, public update_pmax_mat(pmax_set, map_atom_to_kind_atom, set_offset, atomic_block_offset, pmax_atom, full_density_alpha, full_density_beta, natom, kind_of, basis_parameter, nkind, is_assoc_atomic_block, nspins)
updates the maximum of the density matrix in compressed form for screening purposes
Definition hfx_screening_methods.F:138
hfx_screening_methods::calc_pair_dist_radii
subroutine, public calc_pair_dist_radii(qs_env, basis_parameter, radii_pgf, max_set, max_pgf, eps_schwarz, n_threads, i_thread)
calculates radius functions for longrange screening
Definition hfx_screening_methods.F:658
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:16
hfx_types::log_zero
real(dp), parameter, public log_zero
Definition hfx_types.F:125
hfx_types::powell_min_log
real(dp), parameter, public powell_min_log
Definition hfx_types.F:126
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::default_output_unit
integer, parameter, public default_output_unit
Definition machine.F:48
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
powell
Definition powell.F:9
powell::powell_optimize
subroutine, public powell_optimize(n, x, optstate)
...
Definition powell.F:52
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:483
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
four_center_integral_interface::four_center_integral_context_type
Context owning the energy and first-derivative backend workspaces.
Definition four_center_integral_interface.F:48
hfx_types::hfx_basis_type
Definition hfx_types.F:260
hfx_types::hfx_p_kind
Definition hfx_types.F:289
powell::opt_state_type
Definition powell.F:18
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190