(git:e7e05ae)
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-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  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 
900 END 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:45
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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
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, 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_r3d_rs_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_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
powell::opt_state_type
Definition: powell.F:18