(git:ed6f26b)
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-2025 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 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
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 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 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 max_val_temp = max_val
96 DO ipgf = 1, npgfa
97 DO jpgf = 1, npgfb
98 r1 = max(0.0_dp, tmp_r_1(jpgf, ipgf)%x(1)*rab2 + tmp_r_1(jpgf, ipgf)%x(2))
99 DO la = la_min, la_max
100 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 max_val_temp, potential_parameter, r1, r1, p_work)
106
107 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 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 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)
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 INTEGER, DIMENSION(:), POINTER :: nsgfb, nsgfc
156 REAL(dp) :: pmax_tmp
157
158 CALL timeset(routinen, handle)
159
160 DO i = 1, SIZE(pmax_set)
161 pmax_set(i)%p_kind = 0.0_dp
162 END DO
163
164 pmax_atom = log_zero
165
166 nimg = SIZE(full_density_alpha, 2)
167
168 DO jatom = 1, natom
169 jkind = kind_of(jatom)
170 nsetb = basis_parameter(jkind)%nset
171 nsgfb => basis_parameter(jkind)%nsgf
172
173 DO katom = 1, natom
174 IF (.NOT. is_assoc_atomic_block(jatom, katom) >= 1) cycle
175 kkind = kind_of(katom)
176 IF (kkind < jkind) cycle
177 nsetc = basis_parameter(kkind)%nset
178 nsgfc => basis_parameter(kkind)%nsgf
179 act_atomic_block_offset = atomic_block_offset(katom, jatom)
180 DO jset = 1, nsetb
181 DO kset = 1, nsetc
182 IF (katom >= jatom) THEN
183 pmax_tmp = 0.0_dp
184 act_set_offset = set_offset(kset, jset, kkind, jkind)
185 DO img = 1, nimg
186 i = act_set_offset + act_atomic_block_offset - 1
187 DO mc = 1, nsgfc(kset)
188 DO mb = 1, nsgfb(jset)
189 pmax_tmp = max(pmax_tmp, abs(full_density_alpha(i, img)))
190 IF (ASSOCIATED(full_density_beta)) THEN
191 pmax_tmp = max(pmax_tmp, abs(full_density_beta(i, img)))
192 END IF
193 i = i + 1
194 END DO
195 END DO
196 END DO
197 IF (pmax_tmp == 0.0_dp) THEN
198 pmax_tmp = log_zero
199 ELSE
200 pmax_tmp = log10(pmax_tmp)
201 END IF
202 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 map_atom_to_kind_atom(katom)) = pmax_tmp
207 ELSE
208 pmax_tmp = 0.0_dp
209 act_set_offset = set_offset(jset, kset, jkind, kkind)
210 DO img = 1, nimg
211 DO mc = 1, nsgfc(kset)
212 i = act_set_offset + act_atomic_block_offset - 1 + mc - 1
213 DO mb = 1, nsgfb(jset)
214 pmax_tmp = max(pmax_tmp, abs(full_density_alpha(i, img)))
215 IF (ASSOCIATED(full_density_beta)) THEN
216 pmax_tmp = max(pmax_tmp, abs(full_density_beta(i, img)))
217 END IF
218 i = i + nsgfc(kset)
219 END DO
220 END DO
221 END DO
222 IF (pmax_tmp == 0.0_dp) THEN
223 pmax_tmp = log_zero
224 ELSE
225 pmax_tmp = log10(pmax_tmp)
226 END IF
227 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 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 DO jatom = 1, natom
239 jkind = kind_of(jatom)
240 nsetb = basis_parameter(jkind)%nset
241 DO katom = 1, natom
242 IF (.NOT. is_assoc_atomic_block(jatom, katom) >= 1) cycle
243 kkind = kind_of(katom)
244 IF (kkind < jkind) cycle
245 nsetc = basis_parameter(kkind)%nset
246 pmax_tmp = log_zero
247 DO jset = 1, nsetb
248 DO kset = 1, nsetc
249 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 map_atom_to_kind_atom(katom)), pmax_tmp)
254 END DO
255 END DO
256 pmax_atom(jatom, katom) = pmax_tmp
257 pmax_atom(katom, jatom) = pmax_tmp
258 END DO
259 END DO
260
261 CALL timestop(handle)
262
263 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 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 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
323 npgfb, nsgfa, nsgfb
324 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 REAL(dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
328 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 POINTER :: tmp_r_1
333 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
334 TYPE(qs_kind_type), POINTER :: qs_kind
335
336!$OMP MASTER
337 CALL timeset(routinen, handle)
338!$OMP END MASTER
339
340 CALL get_qs_env(qs_env=qs_env, &
341 qs_kind_set=qs_kind_set)
342
343 nkind = SIZE(qs_kind_set, 1)
344
345!$OMP MASTER
346 ALLOCATE (coeffs_pgf(max_pgf, max_pgf, max_set, max_set, nkind, nkind))
347
348 DO ikind = 1, nkind
349 DO jkind = 1, nkind
350 DO iset = 1, max_set
351 DO jset = 1, max_set
352 DO ipgf = 1, max_pgf
353 DO jpgf = 1, max_pgf
354 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!$OMP BARRIER
364 ra = 0.0_dp
365 rb = 0.0_dp
366 DO ikind = 1, nkind
367 NULLIFY (qs_kind, orb_basis)
368 qs_kind => qs_kind_set(ikind)
369 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
370 NULLIFY (la_max, la_min, npgfa, zeta)
371
372 la_max => basis_parameter(ikind)%lmax
373 la_min => basis_parameter(ikind)%lmin
374 npgfa => basis_parameter(ikind)%npgf
375 nseta = basis_parameter(ikind)%nset
376 zeta => basis_parameter(ikind)%zet
377 set_radius_a => basis_parameter(ikind)%set_radius
378 first_sgfa => basis_parameter(ikind)%first_sgf
379 sphi_a => basis_parameter(ikind)%sphi
380 nsgfa => basis_parameter(ikind)%nsgf
381 rpgfa => basis_parameter(ikind)%pgf_radius
382
383 DO jkind = 1, nkind
384 NULLIFY (qs_kind, orb_basis)
385 qs_kind => qs_kind_set(jkind)
386 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
387 NULLIFY (lb_max, lb_min, npgfb, zetb)
388
389 lb_max => basis_parameter(jkind)%lmax
390 lb_min => basis_parameter(jkind)%lmin
391 npgfb => basis_parameter(jkind)%npgf
392 nsetb = basis_parameter(jkind)%nset
393 zetb => basis_parameter(jkind)%zet
394 set_radius_b => basis_parameter(jkind)%set_radius
395 first_sgfb => basis_parameter(jkind)%first_sgf
396 sphi_b => basis_parameter(jkind)%sphi
397 nsgfb => basis_parameter(jkind)%nsgf
398 rpgfb => basis_parameter(jkind)%pgf_radius
399
400 DO iset = 1, nseta
401 ncoa = npgfa(iset)*ncoset(la_max(iset))
402 sgfa = first_sgfa(1, iset)
403 max_contraction_a = maxval((/(sum(abs(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/))
404 DO jset = 1, nsetb
405 ncob = npgfb(jset)*ncoset(lb_max(jset))
406 sgfb = first_sgfb(1, jset)
407 max_contraction_b = maxval((/(sum(abs(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/))
408 radius = set_radius_a(iset) + set_radius_b(jset)
409 DO ipgf = 1, npgfa(iset)
410 DO jpgf = 1, npgfb(jset)
411 radius = rpgfa(ipgf, iset) + rpgfb(jpgf, jset)
412 DO i = i_thread, 100, n_threads
413 rb(1) = 0.0_dp + real(i, dp)*0.01_dp*radius
414 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 radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(2))
417 DO la = la_min(iset), la_max(iset)
418 DO lb = lb_min(jset), lb_max(jset)
419 !Build primitives
420 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 max_val_temp, potential_parameter, r1, r1, p_work)
425 max_val = max(max_val, max_val_temp)
426 END DO !lb
427 END DO !la
428 max_val = sqrt(max_val)
429 max_val = max_val*max_contraction_a*max_contraction_b
430 DATA(1, i) = rb(1)
431 IF (max_val == 0.0_dp) THEN
432 DATA(2, i) = powell_min_log
433 ELSE
434 DATA(2, i) = log10((max_val))
435 END IF
436 END DO
437!$OMP BARRIER
438!$OMP MASTER
439 CALL optimize_it(DATA, x, powell_min_log)
440 coeffs_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x = x
441!$OMP END MASTER
442!$OMP BARRIER
443
444 END DO
445 END DO
446 END DO
447 END DO
448 END DO
449 END DO
450
451!$OMP MASTER
452 ALLOCATE (coeffs_set(max_set, max_set, nkind, nkind))
453
454 DO ikind = 1, nkind
455 DO jkind = 1, nkind
456 DO iset = 1, max_set
457 DO jset = 1, max_set
458 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!$OMP BARRIER
465 ra = 0.0_dp
466 rb = 0.0_dp
467 DO ikind = 1, nkind
468 NULLIFY (qs_kind, orb_basis)
469 qs_kind => qs_kind_set(ikind)
470 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
471 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 la_max => basis_parameter(ikind)%lmax
483 la_min => basis_parameter(ikind)%lmin
484 npgfa => basis_parameter(ikind)%npgf
485 nseta = basis_parameter(ikind)%nset
486 zeta => basis_parameter(ikind)%zet
487 set_radius_a => basis_parameter(ikind)%set_radius
488 first_sgfa => basis_parameter(ikind)%first_sgf
489 sphi_a => basis_parameter(ikind)%sphi
490 nsgfa => basis_parameter(ikind)%nsgf
491
492 DO jkind = 1, nkind
493 NULLIFY (qs_kind, orb_basis)
494 qs_kind => qs_kind_set(jkind)
495 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
496 NULLIFY (lb_max, lb_min, npgfb, zetb)
497
498 lb_max => basis_parameter(jkind)%lmax
499 lb_min => basis_parameter(jkind)%lmin
500 npgfb => basis_parameter(jkind)%npgf
501 nsetb = basis_parameter(jkind)%nset
502 zetb => basis_parameter(jkind)%zet
503 set_radius_b => basis_parameter(jkind)%set_radius
504 first_sgfb => basis_parameter(jkind)%first_sgf
505 sphi_b => basis_parameter(jkind)%sphi
506 nsgfb => basis_parameter(jkind)%nsgf
507
508 DO iset = 1, nseta
509 ncoa = npgfa(iset)*ncoset(la_max(iset))
510 sgfa = first_sgfa(1, iset)
511 max_contraction_a = maxval((/(sum(abs(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/))
512 DO jset = 1, nsetb
513 ncob = npgfb(jset)*ncoset(lb_max(jset))
514 sgfb = first_sgfb(1, jset)
515 max_contraction_b = maxval((/(sum(abs(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/))
516 radius = set_radius_a(iset) + set_radius_b(jset)
517 tmp_r_1 => radii_pgf(:, :, jset, iset, jkind, ikind)
518 DO i = i_thread, 100, n_threads
519 rb(1) = 0.0_dp + real(i, dp)*0.01_dp*radius
520 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 max_val, potential_parameter, tmp_r_1, rb(1)**2, p_work)
526 max_val = sqrt(max_val)
527 max_val = max_val*max_contraction_a*max_contraction_b
528 DATA(1, i) = rb(1)
529 IF (max_val == 0.0_dp) THEN
530 DATA(2, i) = powell_min_log
531 ELSE
532 DATA(2, i) = log10((max_val))
533 END IF
534 END DO
535!$OMP BARRIER
536!$OMP MASTER
537 CALL optimize_it(DATA, x, powell_min_log)
538 coeffs_set(jset, iset, jkind, ikind)%x = x
539!$OMP END MASTER
540!$OMP BARRIER
541 END DO
542 END DO
543
544 END DO
545 END DO
546
547 ! ** now kinds
548!$OMP MASTER
549 ALLOCATE (coeffs_kind(nkind, nkind))
550
551 DO ikind = 1, nkind
552 DO jkind = 1, nkind
553 coeffs_kind(jkind, ikind)%x(:) = 0.0_dp
554 END DO
555 END DO
556!$OMP END MASTER
557 ra = 0.0_dp
558 rb = 0.0_dp
559 DO ikind = 1, nkind
560 NULLIFY (qs_kind, orb_basis)
561 qs_kind => qs_kind_set(ikind)
562 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
563 NULLIFY (la_max, la_min, npgfa, zeta)
564
565 la_max => basis_parameter(ikind)%lmax
566 la_min => basis_parameter(ikind)%lmin
567 npgfa => basis_parameter(ikind)%npgf
568 nseta = basis_parameter(ikind)%nset
569 zeta => basis_parameter(ikind)%zet
570 kind_radius_a = basis_parameter(ikind)%kind_radius
571 first_sgfa => basis_parameter(ikind)%first_sgf
572 sphi_a => basis_parameter(ikind)%sphi
573 nsgfa => basis_parameter(ikind)%nsgf
574
575 DO jkind = 1, nkind
576 NULLIFY (qs_kind, orb_basis)
577 qs_kind => qs_kind_set(jkind)
578 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
579 NULLIFY (lb_max, lb_min, npgfb, zetb)
580
581 lb_max => basis_parameter(jkind)%lmax
582 lb_min => basis_parameter(jkind)%lmin
583 npgfb => basis_parameter(jkind)%npgf
584 nsetb = basis_parameter(jkind)%nset
585 zetb => basis_parameter(jkind)%zet
586 kind_radius_b = basis_parameter(jkind)%kind_radius
587 first_sgfb => basis_parameter(jkind)%first_sgf
588 sphi_b => basis_parameter(jkind)%sphi
589 nsgfb => basis_parameter(jkind)%nsgf
590
591 radius = kind_radius_a + kind_radius_b
592 DO iset = 1, nseta
593 ncoa = npgfa(iset)*ncoset(la_max(iset))
594 sgfa = first_sgfa(1, iset)
595 max_contraction_a = maxval((/(sum(abs(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/))
596 DO jset = 1, nsetb
597 ncob = npgfb(jset)*ncoset(lb_max(jset))
598 sgfb = first_sgfb(1, jset)
599 max_contraction_b = maxval((/(sum(abs(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/))
600 DO i = i_thread, 100, n_threads
601 rb(1) = 0.0_dp + real(i, dp)*0.01_dp*radius
602 max_val = 0.0_dp
603 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 max_val, potential_parameter, tmp_r_1, rb(1)**2, p_work)
609 DATA(1, i) = rb(1)
610 max_val = sqrt(max_val)
611 max_val = max_val*max_contraction_a*max_contraction_b
612 IF (max_val == 0.0_dp) THEN
613 DATA(2, i) = max(powell_min_log, DATA(2, i))
614 ELSE
615 DATA(2, i) = max(log10(max_val), DATA(2, i))
616 END IF
617 END DO
618 END DO
619 END DO
620!$OMP BARRIER
621!$OMP MASTER
622 CALL optimize_it(DATA, x, powell_min_log)
623 coeffs_kind(jkind, ikind)%x = x
624!$OMP END MASTER
625!$OMP BARRIER
626 END DO
627 END DO
628
629!$OMP MASTER
630 CALL timestop(handle)
631!$OMP END MASTER
632
633 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 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 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
675 npgfb, nsgfa, nsgfb
676 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 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 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
683 TYPE(qs_kind_type), POINTER :: qs_kind
684
685!$OMP MASTER
686 CALL timeset(routinen, handle)
687!$OMP END MASTER
688 CALL get_qs_env(qs_env=qs_env, &
689 qs_kind_set=qs_kind_set)
690
691 nkind = SIZE(qs_kind_set, 1)
692 ra = 0.0_dp
693 rb = 0.0_dp
694!$OMP MASTER
695 ALLOCATE (radii_pgf(max_pgf, max_pgf, max_set, max_set, nkind, nkind))
696 DO ikind = 1, nkind
697 DO jkind = 1, nkind
698 DO iset = 1, max_set
699 DO jset = 1, max_set
700 DO ipgf = 1, max_pgf
701 DO jpgf = 1, max_pgf
702 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 DATA = 0.0_dp
711!$OMP END MASTER
712!$OMP BARRIER
713
714 DO ikind = 1, nkind
715 NULLIFY (qs_kind, orb_basis)
716 qs_kind => qs_kind_set(ikind)
717 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
718 NULLIFY (la_max, la_min, npgfa, zeta)
719
720 la_max => basis_parameter(ikind)%lmax
721 la_min => basis_parameter(ikind)%lmin
722 npgfa => basis_parameter(ikind)%npgf
723 nseta = basis_parameter(ikind)%nset
724 zeta => basis_parameter(ikind)%zet
725 first_sgfa => basis_parameter(ikind)%first_sgf
726 sphi_a => basis_parameter(ikind)%sphi
727 nsgfa => basis_parameter(ikind)%nsgf
728 rpgfa => basis_parameter(ikind)%pgf_radius
729
730 DO jkind = 1, nkind
731 NULLIFY (qs_kind, orb_basis)
732 qs_kind => qs_kind_set(jkind)
733 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
734 NULLIFY (lb_max, lb_min, npgfb, zetb)
735
736 lb_max => basis_parameter(jkind)%lmax
737 lb_min => basis_parameter(jkind)%lmin
738 npgfb => basis_parameter(jkind)%npgf
739 nsetb = basis_parameter(jkind)%nset
740 zetb => basis_parameter(jkind)%zet
741 first_sgfb => basis_parameter(jkind)%first_sgf
742 sphi_b => basis_parameter(jkind)%sphi
743 nsgfb => basis_parameter(jkind)%nsgf
744 rpgfb => basis_parameter(jkind)%pgf_radius
745
746 DO iset = 1, nseta
747 ncoa = npgfa(iset)*ncoset(la_max(iset))
748 sgfa = first_sgfa(1, iset)
749 max_contraction_a = maxval((/(sum(abs(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/))
750 DO jset = 1, nsetb
751 ncob = npgfb(jset)*ncoset(lb_max(jset))
752 sgfb = first_sgfb(1, jset)
753 max_contraction_b = maxval((/(sum(abs(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/))
754 DO ipgf = 1, npgfa(iset)
755 DO jpgf = 1, npgfb(jset)
756 radius = rpgfa(ipgf, iset) + rpgfb(jpgf, jset)
757 DO i = i_thread, 100, n_threads
758 rb(1) = 0.0_dp + 0.01_dp*radius*i
759 r_max = 0.0_dp
760 DO la = la_min(iset), la_max(iset)
761 DO lb = lb_min(jset), lb_max(jset)
762 zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
763 ff = zetb(jpgf, jset)/zetp
764 rab = 0.0_dp
765 rab(1) = rb(1)
766 rab2 = rb(1)**2
767 prefactor = exp(-zeta(ipgf, iset)*ff*rab2)
768 rap(:) = ff*rab(:)
769 rp(:) = ra(:) + rap(:)
770 rb(:) = ra(:) + rab(:)
771 cutoff = 1.0_dp
772 r1 = exp_radius_very_extended( &
773 la, la, lb, lb, ra=ra, rb=rb, rp=rp, &
774 zetp=zetp, eps=eps_schwarz, prefactor=prefactor, cutoff=cutoff, epsabs=1.0e-12_dp)
775 r_max = max(r_max, r1)
776 END DO
777 END DO
778 DATA(1, i) = rb(1)
779 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!$OMP BARRIER
784!$OMP MASTER
785 CALL optimize_it(DATA, x, 0.0_dp)
786 radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x = x
787!$OMP END MASTER
788!$OMP BARRIER
789 END DO !jpgf
790 END DO !ipgf
791 END DO
792 END DO
793 END DO
794 END DO
795!$OMP MASTER
796 CALL timestop(handle)
797!$OMP END MASTER
798 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 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 x(1) = 0.0_dp
827 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 DO k = 0, 4, 4
834
835 small_weight = 1.0_dp
836 large_weight = small_weight*(10.0_dp**k)
837
838 ! init opt run
839 opt_state%state = 0
840 opt_state%nvar = 2
841 opt_state%iprint = 3
842 opt_state%unit = default_output_unit
843 opt_state%maxfun = 100000
844 opt_state%rhobeg = 0.1_dp
845 opt_state%rhoend = 0.000001_dp
846
847 DO
848
849 ! compute function
850 IF (opt_state%state == 2) THEN
851 opt_state%f = 0.0_dp
852 DO i = 0, 100
853 f = x(1)*DATA(1, i)**2 + x(2)
854 IF (f > DATA(2, i)) THEN
855 weight = small_weight
856 ELSE
857 weight = large_weight
858 END IF
859 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 IF (opt_state%state == -1) EXIT
864 CALL powell_optimize(opt_state%nvar, x, opt_state)
865 END DO
866
867 ! dealloc mem
868 opt_state%state = 8
869 CALL powell_optimize(opt_state%nvar, x, opt_state)
870
871 END DO
872
873 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 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 min_ij = min(i, j)
896 get_1d_idx = min_ij*n + max(i, j) - (min_ij - 1)*min_ij/2 - n
897
898 END FUNCTION get_1d_idx
899
900END 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:209
basis_set_types
Definition basis_set_types.F:15
hfx_libint_interface
Interface to the Libint-Library.
Definition hfx_libint_interface.F:15
hfx_libint_interface::evaluate_eri_screen
subroutine, public evaluate_eri_screen(libint, 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, p_work)
Evaluates max-abs values of electron repulsion integrals for a primitive quartet.
Definition hfx_libint_interface.F:1300
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, lib, potential_parameter, coeffs_set, coeffs_kind, coeffs_pgf, radii_pgf, max_set, max_pgf, n_threads, i_thread, p_work)
calculates screening functions for schwarz screening
Definition hfx_screening_methods.F:304
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)
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:660
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:15
hfx_types::log_zero
real(dp), parameter, public log_zero
Definition hfx_types.F:118
hfx_types::powell_min_log
real(dp), parameter, public powell_min_log
Definition hfx_types.F:119
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
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition libint_wrapper.F:15
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:53
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:55
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, 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, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, 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, 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, init_u_ramping_each_scf, reltmat, ghost, 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:454
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
hfx_types::hfx_basis_type
Definition hfx_types.F:263
hfx_types::hfx_p_kind
Definition hfx_types.F:296
hfx_types::hfx_potential_type
Definition hfx_types.F:138
hfx_types::hfx_screen_coeff_type
Definition hfx_types.F:291
libint_wrapper::cp_libint_t
Definition libint_wrapper.F:69
powell::opt_state_type
Definition powell.F:18
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171