(git:34ef472)
lri_environment_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 Calculates integral matrices for LRIGPW method
10 !> lri : local resolution of the identity
11 !> \par History
12 !> created JGH [08.2012]
13 !> Dorothea Golze [02.2014] (1) extended, re-structured, cleaned
14 !> (2) heavily debugged
15 !> \authors JGH
16 !> Dorothea Golze
17 ! **************************************************************************************************
18 MODULE lri_environment_methods
19  USE atomic_kind_types, ONLY: atomic_kind_type,&
20  get_atomic_kind,&
21  get_atomic_kind_set
22  USE basis_set_types, ONLY: gto_basis_set_type
23  USE cell_types, ONLY: cell_type
24  USE core_ppl, ONLY: build_core_ppl_ri
25  USE cp_control_types, ONLY: dft_control_type
26  USE cp_log_handling, ONLY: cp_get_default_logger,&
27  cp_logger_type
28  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
29  cp_print_key_unit_nr
30  USE dbcsr_api, ONLY: dbcsr_create,&
31  dbcsr_dot,&
32  dbcsr_get_block_diag,&
33  dbcsr_get_block_p,&
34  dbcsr_p_type,&
35  dbcsr_release,&
36  dbcsr_replicate_all,&
37  dbcsr_type
38  USE generic_os_integrals, ONLY: int_overlap_aabb_os
39  USE input_constants, ONLY: do_lri_inv,&
40  do_lri_inv_auto,&
41  do_lri_pseudoinv_diag,&
42  do_lri_pseudoinv_svd
43  USE input_section_types, ONLY: section_vals_type
44  USE kinds, ONLY: dp
45  USE lri_compression, ONLY: lri_comp,&
46  lri_cont_mem,&
47  lri_decomp_i
48  USE lri_environment_types, ONLY: &
49  allocate_lri_coefs, allocate_lri_ints, allocate_lri_ints_rho, allocate_lri_ppl_ints, &
50  allocate_lri_rhos, deallocate_lri_ints, deallocate_lri_ints_rho, deallocate_lri_ppl_ints, &
51  lri_density_create, lri_density_release, lri_density_type, lri_environment_type, &
52  lri_int_rho_type, lri_int_type, lri_kind_type, lri_list_type, lri_rhoab_type
53  USE lri_integrals, ONLY: allocate_int_type,&
54  deallocate_int_type,&
55  int_type,&
56  lri_int2
57  USE mathlib, ONLY: get_pseudo_inverse_diag,&
58  get_pseudo_inverse_svd,&
59  invmat_symm
60  USE message_passing, ONLY: mp_para_env_type
61  USE particle_types, ONLY: particle_type
62  USE pw_types, ONLY: pw_c1d_gs_type,&
63  pw_r3d_rs_type
64  USE qs_collocate_density, ONLY: calculate_lri_rho_elec
65  USE qs_environment_types, ONLY: get_qs_env,&
66  qs_environment_type
67  USE qs_force_types, ONLY: qs_force_type
68  USE qs_kind_types, ONLY: qs_kind_type
69  USE qs_neighbor_list_types, ONLY: get_iterator_info,&
70  neighbor_list_iterate,&
71  neighbor_list_iterator_create,&
72  neighbor_list_iterator_p_type,&
73  neighbor_list_iterator_release,&
74  neighbor_list_set_p_type
75  USE qs_rho_types, ONLY: qs_rho_get,&
76  qs_rho_set,&
77  qs_rho_type
78  USE virial_types, ONLY: virial_type
79 
80 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
81 #include "./base/base_uses.f90"
82 
83  IMPLICIT NONE
84 
85  PRIVATE
86 
87 ! **************************************************************************************************
88 
89  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_environment_methods'
90 
91  PUBLIC :: build_lri_matrices, calculate_lri_densities, &
92  calculate_lri_integrals, calculate_avec_lri, &
93  v_int_ppl_update, v_int_ppl_energy, lri_kg_rho_update, lri_print_stat
94 
95 ! **************************************************************************************************
96 
97 CONTAINS
98 
99 ! **************************************************************************************************
100 !> \brief creates and initializes an lri_env
101 !> \param lri_env the lri_environment you want to create
102 !> \param qs_env ...
103 ! **************************************************************************************************
104  SUBROUTINE build_lri_matrices(lri_env, qs_env)
105 
106  TYPE(lri_environment_type), POINTER :: lri_env
107  TYPE(qs_environment_type), POINTER :: qs_env
108 
109  ! calculate the integrals needed to do the local (2-center) expansion
110  ! of the (pair) densities
111  CALL calculate_lri_integrals(lri_env, qs_env)
112 
113  ! calculate integrals for local pp (if used in LRI)
114  IF (lri_env%ppl_ri) THEN
115  CALL calculate_lri_ppl_integrals(lri_env, qs_env, .false.)
116  END IF
117 
118  END SUBROUTINE build_lri_matrices
119 
120 ! **************************************************************************************************
121 !> \brief calculates integrals needed for the LRI density fitting,
122 !> integrals are calculated once, before the SCF starts
123 !> \param lri_env ...
124 !> \param qs_env ...
125 ! **************************************************************************************************
126  SUBROUTINE calculate_lri_integrals(lri_env, qs_env)
127 
128  TYPE(lri_environment_type), POINTER :: lri_env
129  TYPE(qs_environment_type), POINTER :: qs_env
130 
131  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_lri_integrals'
132 
133  INTEGER :: handle, iac, iatom, ikind, ilist, jatom, &
134  jkind, jneighbor, mepos, nba, nbb, &
135  nfa, nfb, nkind, nlist, nn, nneighbor, &
136  nthread
137  LOGICAL :: e1c, use_virial
138  REAL(kind=dp) :: cmem, cpff, cpsr, cptt, dab
139  REAL(kind=dp), DIMENSION(3) :: rab
140  TYPE(cell_type), POINTER :: cell
141  TYPE(dft_control_type), POINTER :: dft_control
142  TYPE(gto_basis_set_type), POINTER :: fbasa, fbasb, obasa, obasb
143  TYPE(int_type) :: lriint
144  TYPE(lri_int_type), POINTER :: lrii
145  TYPE(lri_list_type), POINTER :: lri_ints
146  TYPE(neighbor_list_iterator_p_type), &
147  DIMENSION(:), POINTER :: nl_iterator
148  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
149  POINTER :: soo_list
150  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
151  TYPE(virial_type), POINTER :: virial
152 
153  CALL timeset(routinen, handle)
154  NULLIFY (cell, dft_control, fbasa, fbasb, lrii, lri_ints, nl_iterator, &
155  obasa, obasb, particle_set, soo_list, virial)
156 
157  lri_env%stat%pairs_tt = 0.0_dp
158  lri_env%stat%pairs_sr = 0.0_dp
159  lri_env%stat%pairs_ff = 0.0_dp
160  lri_env%stat%overlap_error = 0.0_dp
161  lri_env%stat%abai_mem = 0.0_dp
162 
163  IF (ASSOCIATED(lri_env%soo_list)) THEN
164  soo_list => lri_env%soo_list
165 
166  CALL get_qs_env(qs_env=qs_env, cell=cell, dft_control=dft_control, &
167  nkind=nkind, particle_set=particle_set, virial=virial)
168  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
169  nthread = 1
170 !$ nthread = omp_get_max_threads()
171 
172  IF (ASSOCIATED(lri_env%lri_ints)) THEN
173  CALL deallocate_lri_ints(lri_env%lri_ints)
174  END IF
175 
176  ! allocate matrices storing the LRI integrals
177  CALL allocate_lri_ints(lri_env, lri_env%lri_ints, nkind)
178  lri_ints => lri_env%lri_ints
179 
180  CALL neighbor_list_iterator_create(nl_iterator, soo_list, nthread=nthread)
181 !$OMP PARALLEL DEFAULT(NONE)&
182 !$OMP SHARED (nthread,nl_iterator,lri_env,lri_ints,nkind,use_virial)&
183 !$OMP PRIVATE (mepos,ikind,jkind,iatom,jatom,nlist,ilist,nneighbor,jneighbor,rab,dab,&
184 !$OMP e1c,iac,obasa,obasb,fbasa,fbasb,lrii,lriint,nba,nbb,nfa,nfb,nn,cptt,cpsr,cpff,cmem)
185 
186  mepos = 0
187 !$ mepos = omp_get_thread_num()
188 
189  cptt = 0.0_dp
190  cpsr = 0.0_dp
191  cpff = 0.0_dp
192  cmem = 0.0_dp
193  DO WHILE (neighbor_list_iterate(nl_iterator, mepos) == 0)
194 
195  CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
196  nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor, &
197  iatom=iatom, jatom=jatom, r=rab)
198  iac = ikind + nkind*(jkind - 1)
199  dab = sqrt(sum(rab*rab))
200 
201  obasa => lri_env%orb_basis(ikind)%gto_basis_set
202  obasb => lri_env%orb_basis(jkind)%gto_basis_set
203  fbasa => lri_env%ri_basis(ikind)%gto_basis_set
204  fbasb => lri_env%ri_basis(jkind)%gto_basis_set
205 
206  IF (.NOT. ASSOCIATED(obasa)) cycle
207  IF (.NOT. ASSOCIATED(obasb)) cycle
208 
209  lrii => lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
210 
211  ! exact 1 center approximation
212  e1c = .false.
213  IF (iatom == jatom .AND. dab < lri_env%delta) e1c = .true.
214  ! forces: not every pair is giving contributions
215  ! no forces for self-pair aa
216  IF (iatom == jatom .AND. dab < lri_env%delta) THEN
217  lrii%calc_force_pair = .false.
218  ELSE
219  ! forces for periodic self-pair aa' required for virial
220  IF (iatom == jatom .AND. .NOT. use_virial) THEN
221  lrii%calc_force_pair = .false.
222  ELSE
223  lrii%calc_force_pair = .true.
224  END IF
225  END IF
226 
227  IF (e1c .AND. lri_env%exact_1c_terms) THEN
228  ! nothing to do
229  ELSE
230  cptt = cptt + 1.0_dp
231 
232  nba = obasa%nsgf
233  nbb = obasb%nsgf
234  nfa = fbasa%nsgf
235  nfb = fbasb%nsgf
236 
237  lrii%nba = nba
238  lrii%nbb = nbb
239  lrii%nfa = nfa
240  lrii%nfb = nfb
241 
242  ! full method is used
243  ! calculate integrals (fa,fb), (a,b), (a,b,fa) and (a,b,fb)
244  CALL allocate_int_type(lriint=lriint, &
245  nba=nba, nbb=nbb, nfa=nfa, nfb=nfb, skip_sab=(.NOT. lrii%lrisr))
246  CALL lri_int2(lri_env, lrii, lriint, rab, obasa, obasb, fbasa, fbasb, &
247  iatom, jatom, ikind, jkind)
248  ! store abaint/abbint in compressed form
249  IF (e1c) THEN
250  CALL lri_comp(lriint%abaint, lrii%abascr, lrii%cabai)
251  cmem = cmem + lri_cont_mem(lrii%cabai)
252  ELSE
253  CALL lri_comp(lriint%abaint, lrii%abascr, lrii%cabai)
254  cmem = cmem + lri_cont_mem(lrii%cabai)
255  CALL lri_comp(lriint%abbint, lrii%abbscr, lrii%cabbi)
256  cmem = cmem + lri_cont_mem(lrii%cabbi)
257  END IF
258  ! store overlap
259  lrii%soo(1:nba, 1:nbb) = lriint%sooint(1:nba, 1:nbb)
260 
261  ! Full LRI method
262  IF (lrii%lrisr) THEN
263  cpsr = cpsr + 1.0_dp
264  ! construct and invert S matrix
265  ! calculate Sinv*n and n*Sinv*n
266  IF (e1c) THEN
267  lrii%sinv(1:nfa, 1:nfa) = lri_env%bas_prop(ikind)%ri_ovlp_inv(1:nfa, 1:nfa)
268  lrii%n(1:nfa) = lri_env%bas_prop(ikind)%int_fbas(1:nfa)
269  CALL dgemv("N", nfa, nfa, 1.0_dp, lrii%sinv(1, 1), nfa, &
270  lrii%n(1), 1, 0.0_dp, lrii%sn, 1)
271  lrii%nsn = sum(lrii%sn(1:nfa)*lrii%n(1:nfa))
272  ELSE
273  nn = nfa + nfb
274  CALL inverse_lri_overlap(lri_env, lrii%sinv, lri_env%bas_prop(ikind)%ri_ovlp, &
275  lri_env%bas_prop(jkind)%ri_ovlp, lriint%sabint)
276  lrii%n(1:nfa) = lri_env%bas_prop(ikind)%int_fbas(1:nfa)
277  lrii%n(nfa + 1:nn) = lri_env%bas_prop(jkind)%int_fbas(1:nfb)
278  CALL dgemv("N", nn, nn, 1.0_dp, lrii%sinv(1, 1), nn, &
279  lrii%n(1), 1, 0.0_dp, lrii%sn, 1)
280  lrii%nsn = sum(lrii%sn(1:nn)*lrii%n(1:nn))
281  END IF
282  IF (lri_env%store_integrals) THEN
283  IF (ALLOCATED(lrii%sab)) DEALLOCATE (lrii%sab)
284  ALLOCATE (lrii%sab(nfa, nfb))
285  lrii%sab(1:nfa, 1:nfb) = lriint%sabint(1:nfa, 1:nfb)
286  END IF
287  END IF
288 
289  ! Distant Pair methods
290  IF (lrii%lriff) THEN
291  cpff = cpff + 1.0_dp
292  cpassert(.NOT. e1c)
293  cpassert(.NOT. lri_env%store_integrals)
294  ! calculate Sinv*n and n*Sinv*n for A and B centers
295  lrii%na(1:nfa) = lri_env%bas_prop(ikind)%int_fbas(1:nfa)
296  lrii%nb(1:nfb) = lri_env%bas_prop(jkind)%int_fbas(1:nfb)
297  CALL dgemv("N", nfa, nfa, 1.0_dp, lrii%asinv(1, 1), nfa, &
298  lrii%na(1), 1, 0.0_dp, lrii%sna, 1)
299  lrii%nsna = sum(lrii%sna(1:nfa)*lrii%na(1:nfa))
300  CALL dgemv("N", nfb, nfb, 1.0_dp, lrii%bsinv(1, 1), nfb, &
301  lrii%nb(1), 1, 0.0_dp, lrii%snb, 1)
302  lrii%nsnb = sum(lrii%snb(1:nfb)*lrii%nb(1:nfb))
303  END IF
304 
305  CALL deallocate_int_type(lriint=lriint)
306 
307  END IF
308 
309  END DO
310 !$OMP CRITICAL(UPDATE)
311  lri_env%stat%pairs_tt = lri_env%stat%pairs_tt + cptt
312  lri_env%stat%pairs_sr = lri_env%stat%pairs_sr + cpsr
313  lri_env%stat%pairs_ff = lri_env%stat%pairs_ff + cpff
314  lri_env%stat%abai_mem = lri_env%stat%abai_mem + cmem
315 !$OMP END CRITICAL(UPDATE)
316 
317 !$OMP END PARALLEL
318 
319  CALL neighbor_list_iterator_release(nl_iterator)
320 
321  IF (lri_env%debug) THEN
322  CALL output_debug_info(lri_env, qs_env, lri_ints, soo_list)
323  END IF
324 
325  END IF
326 
327  CALL timestop(handle)
328 
329  END SUBROUTINE calculate_lri_integrals
330 
331 ! **************************************************************************************************
332 !> \brief ...
333 !> \param rho_struct ...
334 !> \param qs_env ...
335 !> \param lri_env ...
336 !> \param lri_density ...
337 !> \param atomlist ...
338 ! **************************************************************************************************
339  SUBROUTINE lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist)
340 
341  TYPE(qs_rho_type), POINTER :: rho_struct
342  TYPE(qs_environment_type), POINTER :: qs_env
343  TYPE(lri_environment_type), POINTER :: lri_env
344  TYPE(lri_density_type), POINTER :: lri_density
345  INTEGER, DIMENSION(:), INTENT(IN) :: atomlist
346 
347  CHARACTER(LEN=*), PARAMETER :: routinen = 'lri_kg_rho_update'
348 
349  INTEGER :: handle, ispin, nspins
350  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
351  TYPE(dbcsr_type) :: pmat_diag
352  TYPE(dft_control_type), POINTER :: dft_control
353  TYPE(mp_para_env_type), POINTER :: para_env
354  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
355  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
356 
357  CALL timeset(routinen, handle)
358 
359  cpassert(ASSOCIATED(rho_struct))
360 
361  CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env)
362 
363  CALL qs_rho_get(rho_struct, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
364 
365  nspins = dft_control%nspins
366 
367  cpassert(.NOT. dft_control%use_kinetic_energy_density)
368  cpassert(.NOT. lri_env%exact_1c_terms)
369 
370  DO ispin = 1, nspins
371  CALL calculate_lri_rho_elec(rho_g(ispin), rho_r(ispin), qs_env, &
372  lri_density%lri_coefs(ispin)%lri_kinds, tot_rho_r(ispin), &
373  "LRI_AUX", lri_env%exact_1c_terms, pmat=pmat_diag, atomlist=atomlist)
374  END DO
375 
376  CALL qs_rho_set(rho_struct, rho_r_valid=.true., rho_g_valid=.true.)
377 
378  CALL timestop(handle)
379 
380  END SUBROUTINE lri_kg_rho_update
381 
382 ! **************************************************************************************************
383 !> \brief calculates integrals needed for the LRI ppl method,
384 !> integrals are calculated once, before the SCF starts
385 !> For forces we assume integrals are available and will not be updated
386 !> \param lri_env ...
387 !> \param qs_env ...
388 !> \param calculate_forces ...
389 ! **************************************************************************************************
390  SUBROUTINE calculate_lri_ppl_integrals(lri_env, qs_env, calculate_forces)
391 
392  TYPE(lri_environment_type), POINTER :: lri_env
393  TYPE(qs_environment_type), POINTER :: qs_env
394  LOGICAL, INTENT(IN) :: calculate_forces
395 
396  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_lri_ppl_integrals'
397 
398  INTEGER :: handle, ikind, ispin, na, nb, nkind, &
399  nspin
400  LOGICAL :: use_virial
401  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
402  TYPE(lri_density_type), POINTER :: lri_density
403  TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_ppl_coef
404  TYPE(mp_para_env_type), POINTER :: para_env
405  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
406  POINTER :: sac_lri
407  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
408  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
409  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
410  TYPE(virial_type), POINTER :: virial
411 
412  CALL timeset(routinen, handle)
413 
414  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, &
415  particle_set=particle_set, atomic_kind_set=atomic_kind_set)
416  IF (calculate_forces) THEN
417  cpassert(ASSOCIATED(lri_env%lri_ppl_ints))
418  CALL get_qs_env(qs_env, lri_density=lri_density)
419  nspin = SIZE(lri_density%lri_coefs, 1)
420  ELSE
421  IF (ASSOCIATED(lri_env%lri_ppl_ints)) THEN
422  CALL deallocate_lri_ppl_ints(lri_env%lri_ppl_ints)
423  END IF
424  CALL allocate_lri_ppl_ints(lri_env, lri_env%lri_ppl_ints, atomic_kind_set)
425  END IF
426  !
427  CALL get_qs_env(qs_env, sac_lri=sac_lri, force=force, virial=virial, para_env=para_env)
428  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
429  use_virial = use_virial .AND. calculate_forces
430  !
431  CALL get_qs_env(qs_env, nkind=nkind)
432  ALLOCATE (lri_ppl_coef(nkind))
433  DO ikind = 1, nkind
434  na = SIZE(lri_env%lri_ppl_ints%lri_ppl(ikind)%v_int, 1)
435  nb = SIZE(lri_env%lri_ppl_ints%lri_ppl(ikind)%v_int, 2)
436  NULLIFY (lri_ppl_coef(ikind)%acoef)
437  NULLIFY (lri_ppl_coef(ikind)%v_int)
438  NULLIFY (lri_ppl_coef(ikind)%v_dadr)
439  NULLIFY (lri_ppl_coef(ikind)%v_dfdr)
440  ALLOCATE (lri_ppl_coef(ikind)%v_int(na, nb))
441  lri_ppl_coef(ikind)%v_int = 0.0_dp
442  IF (calculate_forces) THEN
443  ALLOCATE (lri_ppl_coef(ikind)%acoef(na, nb))
444  lri_ppl_coef(ikind)%acoef = 0.0_dp
445  DO ispin = 1, nspin
446  lri_ppl_coef(ikind)%acoef(1:na, 1:nb) = lri_ppl_coef(ikind)%acoef(1:na, 1:nb) + &
447  lri_density%lri_coefs(ispin)%lri_kinds(ikind)%acoef(1:na, 1:nb)
448  END DO
449  END IF
450  END DO
451  !
452  CALL build_core_ppl_ri(lri_ppl_coef, force, virial, calculate_forces, use_virial, &
453  qs_kind_set, atomic_kind_set, particle_set, sac_lri, &
454  "LRI_AUX")
455  !
456  IF (.NOT. calculate_forces) THEN
457  DO ikind = 1, nkind
458  CALL para_env%sum(lri_ppl_coef(ikind)%v_int)
459  lri_env%lri_ppl_ints%lri_ppl(ikind)%v_int = lri_ppl_coef(ikind)%v_int
460  END DO
461  END IF
462  !
463  DO ikind = 1, nkind
464  IF (ASSOCIATED(lri_ppl_coef(ikind)%acoef)) DEALLOCATE (lri_ppl_coef(ikind)%acoef)
465  IF (ASSOCIATED(lri_ppl_coef(ikind)%v_int)) DEALLOCATE (lri_ppl_coef(ikind)%v_int)
466  IF (ASSOCIATED(lri_ppl_coef(ikind)%v_dadr)) DEALLOCATE (lri_ppl_coef(ikind)%v_dadr)
467  IF (ASSOCIATED(lri_ppl_coef(ikind)%v_dfdr)) DEALLOCATE (lri_ppl_coef(ikind)%v_dfdr)
468  END DO
469  DEALLOCATE (lri_ppl_coef)
470 
471  CALL timestop(handle)
472 
473  END SUBROUTINE calculate_lri_ppl_integrals
474 
475 ! **************************************************************************************************
476 !> \brief calculates overlap integrals (aabb) of the orbital basis set,
477 !> required for LRI basis set optimization
478 !> \param lri_env ...
479 !> \param qs_env ...
480 ! **************************************************************************************************
481  SUBROUTINE calculate_lri_overlap_aabb(lri_env, qs_env)
482 
483  TYPE(lri_environment_type), POINTER :: lri_env
484  TYPE(qs_environment_type), POINTER :: qs_env
485 
486  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_lri_overlap_aabb'
487 
488  INTEGER :: handle, iac, iatom, ikind, ilist, jatom, &
489  jkind, jneighbor, nba, nbb, nkind, &
490  nlist, nneighbor
491  REAL(kind=dp) :: dab
492  REAL(kind=dp), DIMENSION(3) :: rab
493  TYPE(cell_type), POINTER :: cell
494  TYPE(gto_basis_set_type), POINTER :: obasa, obasb
495  TYPE(lri_int_rho_type), POINTER :: lriir
496  TYPE(lri_list_type), POINTER :: lri_ints_rho
497  TYPE(neighbor_list_iterator_p_type), &
498  DIMENSION(:), POINTER :: nl_iterator
499  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
500  POINTER :: soo_list
501  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
502 
503  CALL timeset(routinen, handle)
504  NULLIFY (cell, lriir, lri_ints_rho, nl_iterator, obasa, obasb, &
505  particle_set, soo_list)
506 
507  IF (ASSOCIATED(lri_env%soo_list)) THEN
508  soo_list => lri_env%soo_list
509 
510  CALL get_qs_env(qs_env=qs_env, nkind=nkind, particle_set=particle_set, &
511  cell=cell)
512 
513  IF (ASSOCIATED(lri_env%lri_ints_rho)) THEN
514  CALL deallocate_lri_ints_rho(lri_env%lri_ints_rho)
515  END IF
516 
517  CALL allocate_lri_ints_rho(lri_env, lri_env%lri_ints_rho, nkind)
518  lri_ints_rho => lri_env%lri_ints_rho
519 
520  CALL neighbor_list_iterator_create(nl_iterator, soo_list)
521  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
522 
523  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
524  nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor, &
525  iatom=iatom, jatom=jatom, r=rab)
526 
527  iac = ikind + nkind*(jkind - 1)
528  dab = sqrt(sum(rab*rab))
529 
530  obasa => lri_env%orb_basis(ikind)%gto_basis_set
531  obasb => lri_env%orb_basis(jkind)%gto_basis_set
532  IF (.NOT. ASSOCIATED(obasa)) cycle
533  IF (.NOT. ASSOCIATED(obasb)) cycle
534 
535  lriir => lri_ints_rho%lri_atom(iac)%lri_node(ilist)%lri_int_rho(jneighbor)
536 
537  nba = obasa%nsgf
538  nbb = obasb%nsgf
539 
540  ! calculate integrals (aa,bb)
541  CALL int_overlap_aabb_os(lriir%soaabb, obasa, obasb, rab, lri_env%debug, &
542  lriir%dmax_aabb)
543 
544  END DO
545 
546  CALL neighbor_list_iterator_release(nl_iterator)
547 
548  END IF
549 
550  CALL timestop(handle)
551 
552  END SUBROUTINE calculate_lri_overlap_aabb
553 
554 ! **************************************************************************************************
555 !> \brief performs the fitting of the density and distributes the fitted
556 !> density on the grid
557 !> \param lri_env the lri environment
558 !> \param lri_density ...
559 !> \param qs_env ...
560 !> \param pmatrix ...
561 !> \param cell_to_index ...
562 !> \param lri_rho_struct ...
563 !> \param atomic_kind_set ...
564 !> \param para_env ...
565 !> \param response_density ...
566 ! **************************************************************************************************
567  SUBROUTINE calculate_lri_densities(lri_env, lri_density, qs_env, pmatrix, cell_to_index, &
568  lri_rho_struct, atomic_kind_set, para_env, response_density)
569 
570  TYPE(lri_environment_type), POINTER :: lri_env
571  TYPE(lri_density_type), POINTER :: lri_density
572  TYPE(qs_environment_type), POINTER :: qs_env
573  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmatrix
574  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
575  TYPE(qs_rho_type), INTENT(INOUT) :: lri_rho_struct
576  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
577  TYPE(mp_para_env_type), POINTER :: para_env
578  LOGICAL, INTENT(IN) :: response_density
579 
580  CALL calculate_avec_lri(lri_env, lri_density, pmatrix, cell_to_index, response_density)
581  CALL distribute_lri_density_on_the_grid(lri_env, lri_density, qs_env, &
582  lri_rho_struct, atomic_kind_set, para_env, &
583  response_density)
584 
585  END SUBROUTINE calculate_lri_densities
586 
587 ! **************************************************************************************************
588 !> \brief performs the fitting of the density; solves the linear system of
589 !> equations; yield the expansion coefficients avec
590 !> \param lri_env the lri environment
591 !> lri_density the environment for the fitting
592 !> pmatrix density matrix
593 !> \param lri_density ...
594 !> \param pmatrix ...
595 !> \param cell_to_index ...
596 !> \param response_density ...
597 ! **************************************************************************************************
598  SUBROUTINE calculate_avec_lri(lri_env, lri_density, pmatrix, cell_to_index, response_density)
599 
600  TYPE(lri_environment_type), POINTER :: lri_env
601  TYPE(lri_density_type), POINTER :: lri_density
602  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmatrix
603  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
604  LOGICAL, INTENT(IN), OPTIONAL :: response_density
605 
606  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_avec_lri'
607 
608  INTEGER :: handle, i, iac, iatom, ic, ikind, ilist, ispin, jatom, jkind, jneighbor, mepos, &
609  nba, nbb, nfa, nfb, nkind, nlist, nn, nneighbor, nspin, nthread
610  INTEGER, DIMENSION(3) :: cell
611  LOGICAL :: found, lresponse, trans, use_cell_mapping
612  REAL(kind=dp) :: dab, rab(3), threshold
613  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: m
614  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: int3, paa, pab, pbb
615  REAL(kind=dp), DIMENSION(:, :), POINTER :: pbij
616  TYPE(dbcsr_type), POINTER :: pmat
617  TYPE(lri_int_type), POINTER :: lrii
618  TYPE(lri_list_type), POINTER :: lri_rho
619  TYPE(lri_rhoab_type), POINTER :: lrho
620  TYPE(neighbor_list_iterator_p_type), &
621  DIMENSION(:), POINTER :: nl_iterator
622  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
623  POINTER :: soo_list
624 
625  CALL timeset(routinen, handle)
626  NULLIFY (lrii, lri_rho, nl_iterator, pbij, pmat, soo_list)
627 
628  lresponse = .false.
629  IF (PRESENT(response_density)) lresponse = response_density
630 
631  IF (ASSOCIATED(lri_env%soo_list)) THEN
632  soo_list => lri_env%soo_list
633 
634  nspin = SIZE(pmatrix, 1)
635  use_cell_mapping = (SIZE(pmatrix, 2) > 1)
636  nkind = lri_env%lri_ints%nkind
637  nthread = 1
638 !$ nthread = omp_get_max_threads()
639 
640  IF (ASSOCIATED(lri_density)) THEN
641  CALL lri_density_release(lri_density)
642  ELSE
643  ALLOCATE (lri_density)
644  END IF
645  CALL lri_density_create(lri_density)
646  lri_density%nspin = nspin
647 
648  ! allocate structure lri_rhos and vectors tvec and avec
649  CALL allocate_lri_rhos(lri_env, lri_density%lri_rhos, nspin, nkind)
650 
651  DO ispin = 1, nspin
652  lri_rho => lri_density%lri_rhos(ispin)%lri_list
653 
654  CALL neighbor_list_iterator_create(nl_iterator, soo_list, nthread=nthread)
655 !$OMP PARALLEL DEFAULT(NONE)&
656 !$OMP SHARED (nthread,nl_iterator,lri_env,lri_rho,pmatrix,nkind,cell_to_index,use_cell_mapping,ispin,&
657 !$OMP lresponse)&
658 !$OMP PRIVATE (mepos,ikind,jkind,iatom,jatom,nlist,ilist,nneighbor,jneighbor,rab,dab,iac,&
659 !$OMP trans,found,pmat,pbij,pab,paa,pbb,int3,lrho,lrii,nba,nbb,nfa,nfb,nn,threshold,i,m,cell,ic)
660 
661  mepos = 0
662 !$ mepos = omp_get_thread_num()
663  DO WHILE (neighbor_list_iterate(nl_iterator, mepos) == 0)
664  CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, &
665  jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor, &
666  r=rab, cell=cell)
667 
668  iac = ikind + nkind*(jkind - 1)
669  dab = sqrt(sum(rab*rab))
670 
671  IF (.NOT. ASSOCIATED(lri_env%lri_ints%lri_atom(iac)%lri_node)) cycle
672  IF (iatom == jatom .AND. dab < lri_env%delta .AND. lri_env%exact_1c_terms) cycle
673 
674  IF (use_cell_mapping) THEN
675  ic = cell_to_index(cell(1), cell(2), cell(3))
676  cpassert(ic > 0)
677  ELSE
678  ic = 1
679  END IF
680  pmat => pmatrix(ispin, ic)%matrix
681  ! get the density matrix Pab
682  ! this needs to be response density for LRI-TDDFT
683  NULLIFY (pbij)
684  IF (iatom <= jatom) THEN
685  CALL dbcsr_get_block_p(matrix=pmat, row=iatom, col=jatom, block=pbij, found=found)
686  trans = .false.
687  ELSE
688  CALL dbcsr_get_block_p(matrix=pmat, row=jatom, col=iatom, block=pbij, found=found)
689  trans = .true.
690  END IF
691  cpassert(found)
692 
693  lrho => lri_rho%lri_atom(iac)%lri_node(ilist)%lri_rhoab(jneighbor)
694  lrii => lri_env%lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
695 
696  nba = lrii%nba
697  nbb = lrii%nbb
698  nfa = lrii%nfa
699  nfb = lrii%nfb
700 
701  nn = nfa + nfb
702 
703  ! compute tvec = SUM_ab Pab *(a,b,x) and charge constraint
704  ALLOCATE (pab(nba, nbb))
705  IF (trans) THEN
706  pab(1:nba, 1:nbb) = transpose(pbij(1:nbb, 1:nba))
707  ELSE
708  pab(1:nba, 1:nbb) = pbij(1:nba, 1:nbb)
709  END IF
710 
711  ALLOCATE (int3(nba, nbb))
712  IF (lrii%lrisr) THEN
713  ! full LRI method
714  lrho%charge = sum(pab(1:nba, 1:nbb)*lrii%soo(1:nba, 1:nbb))
715  DO i = 1, nfa
716  CALL lri_decomp_i(int3, lrii%cabai, i)
717  lrho%tvec(i) = sum(pab(1:nba, 1:nbb)*int3(1:nba, 1:nbb))
718  END DO
719  IF (dab > lri_env%delta) THEN
720  DO i = 1, nfb
721  CALL lri_decomp_i(int3, lrii%cabbi, i)
722  lrho%tvec(nfa + i) = sum(pab(1:nba, 1:nbb)*int3(1:nba, 1:nbb))
723  END DO
724  END IF
725  !
726  IF (iatom == jatom .AND. dab < lri_env%delta) THEN
727  lrho%nst = sum(lrho%tvec(1:nfa)*lrii%sn(1:nfa))
728  ELSE
729  lrho%nst = sum(lrho%tvec(1:nn)*lrii%sn(1:nn))
730  END IF
731  lrho%lambda = (lrho%charge - lrho%nst)/lrii%nsn
732  !
733  ! solve the linear system of equations
734  ALLOCATE (m(nn))
735  m = 0._dp
736  IF (iatom == jatom .AND. dab < lri_env%delta) THEN
737  m(1:nfa) = lrho%tvec(1:nfa) + lrho%lambda*lrii%n(1:nfa)
738  CALL dgemv("N", nfa, nfa, 1.0_dp, lrii%sinv(1, 1), nfa, &
739  m(1), 1, 0.0_dp, lrho%avec, 1)
740  ELSE
741  m(1:nn) = lrho%tvec(1:nn) + lrho%lambda*lrii%n(1:nn)
742  CALL dgemv("N", nn, nn, 1.0_dp, lrii%sinv(1, 1), nn, &
743  m(1), 1, 0.0_dp, lrho%avec, 1)
744  END IF
745  DEALLOCATE (m)
746  END IF
747 
748  IF (lrii%lriff) THEN
749  ! distant pair approximations
750  ALLOCATE (paa(nba, nbb), pbb(nba, nbb))
751  paa(1:nba, 1:nbb) = pab(1:nba, 1:nbb)*lri_env%wmat(ikind, jkind)%mat(1:nba, 1:nbb)
752  pbb(1:nba, 1:nbb) = pab(1:nba, 1:nbb)*(1._dp - lri_env%wmat(ikind, jkind)%mat(1:nba, 1:nbb))
753  !
754  threshold = lri_env%eps_o3_int/max(sum(abs(paa(1:nba, 1:nbb))), 1.0e-14_dp)
755  lrho%chargea = sum(paa(1:nba, 1:nbb)*lrii%soo(1:nba, 1:nbb))
756  DO i = 1, nfa
757  IF (lrii%abascr(i) > threshold) THEN
758  CALL lri_decomp_i(int3, lrii%cabai, i)
759  lrho%tveca(i) = sum(paa(1:nba, 1:nbb)*int3(1:nba, 1:nbb))
760  ELSE
761  lrho%tveca(i) = 0.0_dp
762  END IF
763  END DO
764  threshold = lri_env%eps_o3_int/max(sum(abs(pbb(1:nba, 1:nbb))), 1.0e-14_dp)
765  lrho%chargeb = sum(pbb(1:nba, 1:nbb)*lrii%soo(1:nba, 1:nbb))
766  DO i = 1, nfb
767  IF (lrii%abbscr(i) > threshold) THEN
768  CALL lri_decomp_i(int3, lrii%cabbi, i)
769  lrho%tvecb(i) = sum(pbb(1:nba, 1:nbb)*int3(1:nba, 1:nbb))
770  ELSE
771  lrho%tvecb(i) = 0.0_dp
772  END IF
773  END DO
774  !
775  lrho%nsta = sum(lrho%tveca(1:nfa)*lrii%sna(1:nfa))
776  lrho%nstb = sum(lrho%tvecb(1:nfb)*lrii%snb(1:nfb))
777  lrho%lambdaa = (lrho%chargea - lrho%nsta)/lrii%nsna
778  lrho%lambdab = (lrho%chargeb - lrho%nstb)/lrii%nsnb
779  ! solve the linear system of equations
780  ALLOCATE (m(nfa))
781  m(1:nfa) = lrho%tveca(1:nfa) + lrho%lambdaa*lrii%na(1:nfa)
782  CALL dgemv("N", nfa, nfa, 1.0_dp, lrii%asinv(1, 1), nfa, &
783  m(1), 1, 0.0_dp, lrho%aveca, 1)
784  DEALLOCATE (m)
785  ALLOCATE (m(nfb))
786  m(1:nfb) = lrho%tvecb(1:nfb) + lrho%lambdab*lrii%nb(1:nfb)
787  CALL dgemv("N", nfb, nfb, 1.0_dp, lrii%bsinv(1, 1), nfb, &
788  m(1), 1, 0.0_dp, lrho%avecb, 1)
789  DEALLOCATE (m)
790  !
791  DEALLOCATE (paa, pbb)
792  END IF
793 
794  DEALLOCATE (pab, int3)
795 
796  END DO
797 !$OMP END PARALLEL
798  CALL neighbor_list_iterator_release(nl_iterator)
799 
800  END DO
801 
802  END IF
803 
804  CALL timestop(handle)
805 
806  END SUBROUTINE calculate_avec_lri
807 
808 ! **************************************************************************************************
809 !> \brief sums up avec and distributes the fitted density on the grid
810 !> \param lri_env the lri environment
811 !> \param lri_density ...
812 !> \param qs_env ...
813 !> \param lri_rho_struct ...
814 !> \param atomic_kind_set ...
815 !> \param para_env ...
816 !> \param response_density ...
817 ! **************************************************************************************************
818  SUBROUTINE distribute_lri_density_on_the_grid(lri_env, lri_density, qs_env, &
819  lri_rho_struct, atomic_kind_set, para_env, &
820  response_density)
821 
822  TYPE(lri_environment_type), POINTER :: lri_env
823  TYPE(lri_density_type), POINTER :: lri_density
824  TYPE(qs_environment_type), POINTER :: qs_env
825  TYPE(qs_rho_type), INTENT(INOUT) :: lri_rho_struct
826  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
827  TYPE(mp_para_env_type), POINTER :: para_env
828  LOGICAL, INTENT(IN) :: response_density
829 
830  CHARACTER(LEN=*), PARAMETER :: routinen = 'distribute_lri_density_on_the_grid'
831 
832  INTEGER :: atom_a, atom_b, handle, iac, iatom, &
833  ikind, ilist, ispin, jatom, jkind, &
834  jneighbor, natom, nfa, nfb, nkind, &
835  nspin
836  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
837  REAL(kind=dp) :: dab, fw, rab(3), str
838  REAL(kind=dp), DIMENSION(:), POINTER :: aci, acj, tot_rho_r
839  TYPE(atomic_kind_type), POINTER :: atomic_kind
840  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
841  TYPE(dbcsr_type) :: pmat_diag
842  TYPE(lri_int_type), POINTER :: lrii
843  TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
844  TYPE(lri_list_type), POINTER :: lri_rho
845  TYPE(lri_rhoab_type), POINTER :: lrho
846  TYPE(neighbor_list_iterator_p_type), &
847  DIMENSION(:), POINTER :: nl_iterator
848  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
849  POINTER :: soo_list
850  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
851  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
852  TYPE(qs_rho_type), POINTER :: rho
853 
854  CALL timeset(routinen, handle)
855  NULLIFY (aci, acj, atomic_kind, lri_coef, lri_rho, &
856  nl_iterator, soo_list, rho_r, rho_g, tot_rho_r)
857 
858  IF (ASSOCIATED(lri_env%soo_list)) THEN
859  soo_list => lri_env%soo_list
860 
861  lri_env%stat%rho_tt = 0.0_dp
862  lri_env%stat%rho_sr = 0.0_dp
863  lri_env%stat%rho_ff = 0.0_dp
864  lri_env%stat%rho_1c = 0.0_dp
865 
866  nspin = lri_density%nspin
867  nkind = lri_env%lri_ints%nkind
868 
869  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
870  atom_of_kind=atom_of_kind)
871 
872  ! allocate the arrays to hold RI expansion coefficients lri_coefs
873  CALL allocate_lri_coefs(lri_env, lri_density, atomic_kind_set)
874  DO ispin = 1, nspin
875 
876  lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
877  lri_rho => lri_density%lri_rhos(ispin)%lri_list
878 
879  ! sum up expansion coefficients
880  CALL neighbor_list_iterator_create(nl_iterator, soo_list)
881  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
882  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
883  iatom=iatom, jatom=jatom, ilist=ilist, inode=jneighbor, r=rab)
884  dab = sqrt(sum(rab*rab))
885  atom_a = atom_of_kind(iatom)
886  atom_b = atom_of_kind(jatom)
887  aci => lri_coef(ikind)%acoef(atom_a, :)
888  acj => lri_coef(jkind)%acoef(atom_b, :)
889  iac = ikind + nkind*(jkind - 1)
890  lrho => lri_rho%lri_atom(iac)%lri_node(ilist)%lri_rhoab(jneighbor)
891  lrii => lri_env%lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
892  nfa = lrho%nfa
893  nfb = lrho%nfb
894 
895  IF (lrii%lrisr) THEN
896  IF (iatom == jatom .AND. dab < lri_env%delta) THEN
897  !self pair aa
898  IF (.NOT. lri_env%exact_1c_terms) THEN
899  aci(1:nfa) = aci(1:nfa) + lrho%avec(1:nfa)
900  lri_env%stat%rho_sr = lri_env%stat%rho_sr + sum(lrho%avec(:)*lrii%n(:))
901  END IF
902  ELSE
903  IF (iatom == jatom) THEN
904  !periodic self pair aa'
905  fw = lrii%wsr
906  ELSE
907  !pairs ab
908  fw = 2.0_dp*lrii%wsr
909  END IF
910  aci(1:nfa) = aci(1:nfa) + fw*lrho%avec(1:nfa)
911  acj(1:nfb) = acj(1:nfb) + fw*lrho%avec(nfa + 1:nfa + nfb)
912  lri_env%stat%rho_sr = lri_env%stat%rho_sr + fw*sum(lrho%avec(:)*lrii%n(:))
913  END IF
914  END IF
915  !
916  IF (lrii%lriff) THEN
917  IF (iatom == jatom) THEN
918  fw = lrii%wff
919  ELSE
920  fw = 2.0_dp*lrii%wff
921  END IF
922  aci(1:nfa) = aci(1:nfa) + fw*lrho%aveca(1:nfa)
923  acj(1:nfb) = acj(1:nfb) + fw*lrho%avecb(1:nfb)
924  lri_env%stat%rho_sr = lri_env%stat%rho_sr + fw*sum(lrho%aveca(:)*lrii%na(:))
925  lri_env%stat%rho_sr = lri_env%stat%rho_sr + fw*sum(lrho%avecb(:)*lrii%nb(:))
926  END IF
927 
928  END DO
929  CALL neighbor_list_iterator_release(nl_iterator)
930 
931  ! replicate the acoef information
932  DO ikind = 1, nkind
933  atomic_kind => atomic_kind_set(ikind)
934  CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom)
935  DO iatom = 1, natom
936  aci => lri_coef(ikind)%acoef(iatom, :)
937  CALL para_env%sum(aci)
938  END DO
939  END DO
940 
941  END DO
942 
943  !distribute fitted density on the grid
944  CALL qs_rho_get(lri_rho_struct, rho_r=rho_r, rho_g=rho_g, &
945  tot_rho_r=tot_rho_r)
946  DO ispin = 1, nspin
947  IF (lri_env%exact_1c_terms) THEN
948  ! 1 center terms (if requested)
949  CALL get_qs_env(qs_env, rho=rho, matrix_s_kp=matrix_s)
950  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
951  CALL dbcsr_create(pmat_diag, name="Block diagonal density", template=matrix_p(1, 1)%matrix)
952  CALL dbcsr_get_block_diag(matrix_p(ispin, 1)%matrix, pmat_diag)
953  str = 0.0_dp
954  CALL dbcsr_dot(matrix_s(1, 1)%matrix, pmat_diag, str)
955  lri_env%stat%rho_1c = lri_env%stat%rho_1c + str
956  CALL dbcsr_replicate_all(pmat_diag)
957  END IF
958  !
959  IF (.NOT. response_density) THEN
960  CALL calculate_lri_rho_elec(rho_g(ispin), &
961  rho_r(ispin), qs_env, &
962  lri_density%lri_coefs(ispin)%lri_kinds, tot_rho_r(ispin), &
963  "LRI_AUX", lri_env%exact_1c_terms, pmat=pmat_diag)
964  ELSE
965  CALL calculate_lri_rho_elec(rho_g(ispin), &
966  rho_r(ispin), qs_env, &
967  lri_density%lri_coefs(ispin)%lri_kinds, tot_rho_r(ispin), &
968  "P_LRI_AUX", lri_env%exact_1c_terms, pmat=pmat_diag)
969  END IF
970  lri_env%stat%rho_tt = lri_env%stat%rho_tt + tot_rho_r(ispin)
971  !
972  IF (lri_env%exact_1c_terms) CALL dbcsr_release(pmat_diag)
973  END DO
974 
975  END IF
976 
977  CALL timestop(handle)
978 
979  END SUBROUTINE distribute_lri_density_on_the_grid
980 
981 ! **************************************************************************************************
982 !> \brief get inverse or pseudoinverse of lri overlap matrix for aux basis set
983 !> \param lri_env lri environment
984 !> \param sinv on exit its inverse
985 !> \param sa ...
986 !> \param sb ...
987 !> \param sab ...
988 ! **************************************************************************************************
989  SUBROUTINE inverse_lri_overlap(lri_env, sinv, sa, sb, sab)
990 
991  TYPE(lri_environment_type) :: lri_env
992  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: sinv
993  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: sa, sb, sab
994 
995  CHARACTER(LEN=*), PARAMETER :: routinen = 'inverse_lri_overlap'
996 
997  INTEGER :: handle, info, n, nfa, nfb, nn
998  INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
999  REAL(kind=dp) :: anorm, delta, rcond, rskip
1000  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
1001  REAL(kind=dp), DIMENSION(:, :), POINTER :: s
1002  REAL(kind=dp), EXTERNAL :: dlange
1003 
1004  CALL timeset(routinen, handle)
1005 
1006  nfa = SIZE(sa, 1)
1007  nfb = SIZE(sb, 1)
1008  nn = nfa + nfb
1009  n = SIZE(sinv, 1)
1010  cpassert(n == nn)
1011  ALLOCATE (s(n, n))
1012  s(1:nfa, 1:nfa) = sa(1:nfa, 1:nfa)
1013  s(1:nfa, nfa + 1:nn) = sab(1:nfa, 1:nfb)
1014  s(nfa + 1:nn, 1:nfa) = transpose(sab(1:nfa, 1:nfb))
1015  s(nfa + 1:nn, nfa + 1:nn) = sb(1:nfb, 1:nfb)
1016 
1017  rskip = 1.e-8_dp ! parameter for pseudo inverse
1018 
1019  sinv(:, :) = s
1020  SELECT CASE (lri_env%lri_overlap_inv)
1021  CASE (do_lri_inv)
1022  CALL invmat_symm(sinv)
1023  CASE (do_lri_pseudoinv_svd)
1024  CALL get_pseudo_inverse_svd(s, sinv, rskip)
1025  CASE (do_lri_pseudoinv_diag)
1026  CALL get_pseudo_inverse_diag(s, sinv, rskip)
1027  CASE (do_lri_inv_auto)
1028  ! decide whether to calculate inverse or pseudoinverse
1029  ALLOCATE (work(3*n), iwork(n))
1030  ! norm of matrix
1031  anorm = dlange('1', n, n, sinv, n, work)
1032  ! Cholesky factorization (fails if matrix not positive definite)
1033  CALL dpotrf('U', n, sinv, n, info)
1034  IF (info == 0) THEN
1035  ! condition number
1036  CALL dpocon('U', n, sinv, n, anorm, rcond, work, iwork, info)
1037  IF (info /= 0) THEN
1038  cpabort("DPOCON failed")
1039  END IF
1040  IF (log(1._dp/rcond) > lri_env%cond_max) THEN
1041  CALL get_pseudo_inverse_svd(s, sinv, rskip)
1042  ELSE
1043  CALL invmat_symm(sinv, "U")
1044  END IF
1045  ELSE
1046  CALL get_pseudo_inverse_svd(s, sinv, rskip)
1047  END IF
1048  DEALLOCATE (work, iwork)
1049  CASE DEFAULT
1050  cpabort("No initialization for LRI overlap available?")
1051  END SELECT
1052 
1053  delta = inv_test(s, sinv)
1054 !$OMP CRITICAL(sum_critical)
1055  lri_env%stat%overlap_error = max(delta, lri_env%stat%overlap_error)
1056 !$OMP END CRITICAL(sum_critical)
1057 
1058  DEALLOCATE (s)
1059 
1060  CALL timestop(handle)
1061 
1062  END SUBROUTINE inverse_lri_overlap
1063 
1064 ! **************************************************************************************************
1065 !> \brief ...
1066 !> \param amat ...
1067 !> \param ainv ...
1068 !> \return ...
1069 ! **************************************************************************************************
1070  FUNCTION inv_test(amat, ainv) RESULT(delta)
1071  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: amat, ainv
1072  REAL(kind=dp) :: delta
1073 
1074  INTEGER :: i, n
1075  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1076 
1077  n = SIZE(amat, 1)
1078  ALLOCATE (work(n, n))
1079  work(1:n, 1:n) = matmul(amat(1:n, 1:n), ainv(1:n, 1:n))
1080  DO i = 1, n
1081  work(i, i) = work(i, i) - 1.0_dp
1082  END DO
1083  delta = maxval(abs(work))
1084  DEALLOCATE (work)
1085  END FUNCTION inv_test
1086 
1087 ! **************************************************************************************************
1088 !> \brief debug output
1089 !> \param lri_env ...
1090 !> \param qs_env ...
1091 !> \param lri_ints ...
1092 !> \param soo_list ...
1093 ! **************************************************************************************************
1094  SUBROUTINE output_debug_info(lri_env, qs_env, lri_ints, soo_list)
1095 
1096  TYPE(lri_environment_type), POINTER :: lri_env
1097  TYPE(qs_environment_type), POINTER :: qs_env
1098  TYPE(lri_list_type), POINTER :: lri_ints
1099  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1100  POINTER :: soo_list
1101 
1102  CHARACTER(LEN=*), PARAMETER :: routinen = 'output_debug_info'
1103 
1104  INTEGER :: handle, iac, ikind, ilist, iunit, jkind, &
1105  jneighbor, nkind
1106  REAL(kind=dp) :: dmax_aabb, dmax_ab, dmax_aba, dmax_abb, &
1107  dmax_oo
1108  TYPE(cp_logger_type), POINTER :: logger
1109  TYPE(dft_control_type), POINTER :: dft_control
1110  TYPE(lri_int_rho_type), POINTER :: lriir
1111  TYPE(lri_int_type), POINTER :: lrii
1112  TYPE(mp_para_env_type), POINTER :: para_env
1113  TYPE(neighbor_list_iterator_p_type), &
1114  DIMENSION(:), POINTER :: nl_iterator
1115  TYPE(section_vals_type), POINTER :: input
1116 
1117  CALL timeset(routinen, handle)
1118  NULLIFY (input, logger, lrii, lriir, nl_iterator, para_env)
1119  CALL get_qs_env(qs_env, dft_control=dft_control, input=input, nkind=nkind, &
1120  para_env=para_env)
1121  dmax_ab = 0._dp
1122  dmax_oo = 0._dp
1123  dmax_aba = 0._dp
1124  dmax_abb = 0._dp
1125  dmax_aabb = 0._dp
1126 
1127  CALL neighbor_list_iterator_create(nl_iterator, soo_list)
1128  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1129 
1130  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
1131  ilist=ilist, inode=jneighbor)
1132 
1133  iac = ikind + nkind*(jkind - 1)
1134  lrii => lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
1135 
1136  dmax_ab = max(dmax_ab, lrii%dmax_ab)
1137  dmax_oo = max(dmax_oo, lrii%dmax_oo)
1138  dmax_aba = max(dmax_aba, lrii%dmax_aba)
1139  dmax_abb = max(dmax_abb, lrii%dmax_abb)
1140 
1141  IF (dft_control%qs_control%lri_optbas) THEN
1142  lriir => lri_env%lri_ints_rho%lri_atom(iac)%lri_node(ilist)%lri_int_rho(jneighbor)
1143  dmax_aabb = max(dmax_aabb, lriir%dmax_aabb)
1144  END IF
1145 
1146  END DO
1147 
1148  CALL neighbor_list_iterator_release(nl_iterator)
1149  CALL para_env%max(dmax_ab)
1150  CALL para_env%max(dmax_oo)
1151  CALL para_env%max(dmax_aba)
1152  CALL para_env%max(dmax_abb)
1153  CALL para_env%max(dmax_aabb)
1154 
1155  logger => cp_get_default_logger()
1156  iunit = cp_print_key_unit_nr(logger, input, "PRINT%PROGRAM_RUN_INFO", &
1157  extension=".lridebug")
1158 
1159  IF (iunit > 0) THEN
1160  WRITE (iunit, fmt="(/,T2,A)") "DEBUG INFO FOR LRI INTEGRALS"
1161  WRITE (iunit, fmt="(T2,A,T69,ES12.5)") "Maximal deviation of integrals "// &
1162  "[ai|bi]; fit basis", dmax_ab
1163  WRITE (iunit, fmt="(T2,A,T69,ES12.5)") "Maximal deviation of integrals "// &
1164  "[a|b]; orbital basis", dmax_oo
1165  WRITE (iunit, fmt="(T2,A,T69,ES12.5)") "Maximal deviation of integrals "// &
1166  "[a|b|ai]", dmax_aba
1167  WRITE (iunit, fmt="(T2,A,T69,ES12.5)") "Maximal deviation of integrals "// &
1168  "[a|b|bi]", dmax_abb
1169  IF (dft_control%qs_control%lri_optbas) THEN
1170  WRITE (iunit, fmt="(T2,A,T69,ES12.5,/)") "Maximal deviation of integrals "// &
1171  "[aa|bb]; orbital basis", &
1172  dmax_aabb
1173  END IF
1174  END IF
1175 
1176  CALL cp_print_key_finished_output(iunit, logger, input, &
1177  "PRINT%PROGRAM_RUN_INFO")
1178  CALL timestop(handle)
1179 
1180  END SUBROUTINE output_debug_info
1181 
1182 ! **************************************************************************************************
1183 !> \brief ...
1184 !> \param qs_env ...
1185 !> \param lri_v_int ...
1186 !> \param calculate_forces ...
1187 ! **************************************************************************************************
1188  SUBROUTINE v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1189  TYPE(qs_environment_type), POINTER :: qs_env
1190  TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
1191  LOGICAL, INTENT(IN) :: calculate_forces
1192 
1193  INTEGER :: ikind, natom, nfa, nkind
1194  REAL(kind=dp), DIMENSION(:, :), POINTER :: v_int
1195  TYPE(lri_environment_type), POINTER :: lri_env
1196 
1197  CALL get_qs_env(qs_env, lri_env=lri_env, nkind=nkind)
1198 
1199  DO ikind = 1, nkind
1200  natom = SIZE(lri_v_int(ikind)%v_int, 1)
1201  nfa = SIZE(lri_v_int(ikind)%v_int, 2)
1202  v_int => lri_env%lri_ppl_ints%lri_ppl(ikind)%v_int
1203  cpassert(SIZE(v_int, 1) == natom)
1204  cpassert(SIZE(v_int, 2) == nfa)
1205  lri_v_int(ikind)%v_int(:, :) = lri_v_int(ikind)%v_int(:, :) + v_int(:, :)
1206  END DO
1207 
1208  IF (calculate_forces) THEN
1209  CALL calculate_lri_ppl_integrals(lri_env, qs_env, calculate_forces)
1210  END IF
1211 
1212  END SUBROUTINE v_int_ppl_update
1213 
1214 ! **************************************************************************************************
1215 !> \brief ...
1216 !> \param qs_env ...
1217 !> \param lri_v_int ...
1218 !> \param ecore_ppl_ri ...
1219 ! **************************************************************************************************
1220  SUBROUTINE v_int_ppl_energy(qs_env, lri_v_int, ecore_ppl_ri)
1221  TYPE(qs_environment_type), POINTER :: qs_env
1222  TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
1223  REAL(kind=dp), INTENT(INOUT) :: ecore_ppl_ri
1224 
1225  INTEGER :: ikind, natom, nfa, nkind
1226  REAL(kind=dp), DIMENSION(:, :), POINTER :: v_int
1227  TYPE(lri_environment_type), POINTER :: lri_env
1228 
1229  CALL get_qs_env(qs_env, lri_env=lri_env, nkind=nkind)
1230  DO ikind = 1, nkind
1231  natom = SIZE(lri_v_int(ikind)%v_int, 1)
1232  nfa = SIZE(lri_v_int(ikind)%v_int, 2)
1233  v_int => lri_env%lri_ppl_ints%lri_ppl(ikind)%v_int
1234  cpassert(SIZE(v_int, 1) == natom)
1235  cpassert(SIZE(v_int, 2) == nfa)
1236  ecore_ppl_ri = ecore_ppl_ri + sum(v_int(:, :)*lri_v_int(ikind)%acoef(:, :))
1237  END DO
1238 
1239  END SUBROUTINE v_int_ppl_energy
1240 
1241 ! **************************************************************************************************
1242 !> \brief ...
1243 !> \param qs_env ...
1244 !> \param ltddfpt ...
1245 !> \param tddfpt_lri_env ...
1246 ! **************************************************************************************************
1247  SUBROUTINE lri_print_stat(qs_env, ltddfpt, tddfpt_lri_env)
1248 
1249  TYPE(qs_environment_type), POINTER :: qs_env
1250  LOGICAL, OPTIONAL :: ltddfpt
1251  TYPE(lri_environment_type), OPTIONAL, POINTER :: tddfpt_lri_env
1252 
1253  INTEGER :: iunit
1254  LOGICAL :: my_ltddfpt
1255  REAL(kind=dp) :: abai_mem, coef_mem, oint_mem, overlap_error, pairs_ff, pairs_sr, pairs_tt, &
1256  ppli_mem, ppx, rho_1c, rho_ff, rho_sr, rho_tt, rhos_mem
1257  TYPE(cp_logger_type), POINTER :: logger
1258  TYPE(lri_environment_type), POINTER :: lri_env
1259  TYPE(mp_para_env_type), POINTER :: para_env
1260  TYPE(section_vals_type), POINTER :: input
1261 
1262  my_ltddfpt = .false.
1263  IF (PRESENT(ltddfpt)) my_ltddfpt = ltddfpt
1264 
1265  IF (.NOT. my_ltddfpt) THEN
1266  CALL get_qs_env(qs_env, lri_env=lri_env, input=input, para_env=para_env)
1267  ELSE
1268  CALL get_qs_env(qs_env, input=input, para_env=para_env)
1269  NULLIFY (lri_env)
1270  lri_env => tddfpt_lri_env
1271  END IF
1272 
1273  IF (lri_env%statistics) THEN
1274  logger => cp_get_default_logger()
1275  iunit = cp_print_key_unit_nr(logger, input, "PRINT%PROGRAM_RUN_INFO", extension=".Log")
1276  pairs_tt = lri_env%stat%pairs_tt
1277  CALL para_env%sum(pairs_tt)
1278  pairs_sr = lri_env%stat%pairs_sr
1279  CALL para_env%sum(pairs_sr)
1280  pairs_ff = lri_env%stat%pairs_ff
1281  CALL para_env%sum(pairs_ff)
1282  overlap_error = lri_env%stat%overlap_error
1283  CALL para_env%sum(overlap_error)
1284  rho_tt = -lri_env%stat%rho_tt
1285  rho_sr = lri_env%stat%rho_sr
1286  CALL para_env%sum(rho_sr)
1287  rho_ff = lri_env%stat%rho_ff
1288  CALL para_env%sum(rho_ff)
1289  IF (lri_env%exact_1c_terms) THEN
1290  rho_1c = lri_env%stat%rho_1c
1291  ELSE
1292  rho_1c = 0.0_dp
1293  END IF
1294  ppx = 8.e-6_dp
1295  coef_mem = lri_env%stat%coef_mem*ppx
1296  CALL para_env%sum(coef_mem)
1297  oint_mem = lri_env%stat%oint_mem*ppx
1298  CALL para_env%sum(oint_mem)
1299  rhos_mem = lri_env%stat%rhos_mem*ppx
1300  CALL para_env%sum(rhos_mem)
1301  abai_mem = lri_env%stat%abai_mem*ppx
1302  CALL para_env%sum(abai_mem)
1303  IF (lri_env%ppl_ri) THEN
1304  ppli_mem = lri_env%stat%ppli_mem*ppx
1305  CALL para_env%sum(ppli_mem)
1306  ELSE
1307  ppli_mem = 0.0_dp
1308  END IF
1309  IF (iunit > 0) THEN
1310  !
1311  WRITE (iunit, fmt="(/,T2,A,A,A)") "********************************", &
1312  " LRI STATISTICS ", "*******************************"
1313  !
1314  WRITE (iunit, fmt="(T4,A,T63,F16.0)") "Total number of atom pairs", pairs_tt
1315  ppx = pairs_sr/pairs_tt*100._dp
1316  WRITE (iunit, fmt="(T4,A,T52,A,F5.1,A,T63,F16.0)") "Near field atom pairs", &
1317  "[", ppx, "%]", pairs_sr
1318  ppx = pairs_ff/pairs_tt*100._dp
1319  WRITE (iunit, fmt="(T4,A,T52,A,F5.1,A,T63,F16.0)") "Far field atom pairs", &
1320  "[", ppx, "%]", pairs_ff
1321  !
1322  WRITE (iunit, fmt="(T4,A,T63,G16.8)") "Max. error in pair overlap inversion", overlap_error
1323  !
1324  WRITE (iunit, fmt="(T4,A,T63,F16.2)") "Total charge approximated", rho_tt
1325  ppx = rho_sr/rho_tt*100._dp
1326  WRITE (iunit, fmt="(T4,A,T52,A,F5.1,A,T63,F16.2)") "Near field charge", &
1327  "[", ppx, "%]", rho_sr
1328  ppx = rho_ff/rho_tt*100._dp
1329  WRITE (iunit, fmt="(T4,A,T52,A,F5.1,A,T63,F16.2)") "Far field charge", &
1330  "[", ppx, "%]", rho_ff
1331  IF (lri_env%exact_1c_terms) THEN
1332  ppx = rho_1c/rho_tt*100._dp
1333  WRITE (iunit, fmt="(T4,A,T52,A,F5.1,A,T63,F16.2)") "One center charge", &
1334  "[", ppx, "%]", rho_1c
1335  END IF
1336  !
1337  WRITE (iunit, fmt="(T4,A,T63,F9.0,A)") "Max. memory/task for expansion coeficients", coef_mem, " Mbytes"
1338  WRITE (iunit, fmt="(T4,A,T63,F9.0,A)") "Max. memory/task for overlap matrices", oint_mem, " Mbytes"
1339  WRITE (iunit, fmt="(T4,A,T63,F9.0,A)") "Max. memory/task for density expansions", rhos_mem, " Mbytes"
1340  WRITE (iunit, fmt="(T4,A,T63,F9.0,A)") "Max. memory/task for aba/abb integrals", abai_mem, " Mbytes"
1341  IF (lri_env%ppl_ri) THEN
1342  WRITE (iunit, fmt="(T4,A,T63,F9.0,A)") "Max. memory/task for ppl integrals", ppli_mem, " Mbytes"
1343  END IF
1344  !
1345  WRITE (iunit, fmt="(T2,A,A)") "********************************", &
1346  "***********************************************"
1347  !
1348  END IF
1349 
1350  CALL cp_print_key_finished_output(iunit, logger, input, "PRINT%PROGRAM_RUN_INFO")
1351  END IF
1352 
1353  END SUBROUTINE lri_print_stat
1354 
1355 END MODULE lri_environment_methods
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition: atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
basis_set_types
Definition: basis_set_types.F:15
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
core_ppl
Calculation of the local pseudopotential contribution to the core Hamiltonian <a|V(local)|b> = <a|Sum...
Definition: core_ppl.F:18
core_ppl::build_core_ppl_ri
subroutine, public build_core_ppl_ri(lri_ppl_coef, force, virial, calculate_forces, use_virial, qs_kind_set, atomic_kind_set, particle_set, sac_ppl, basis_type)
...
Definition: core_ppl.F:743
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
generic_os_integrals
Calculation of contracted, spherical Gaussian integrals using the (OS) integral scheme....
Definition: generic_os_integrals.F:18
generic_os_integrals::int_overlap_aabb_os
subroutine, public int_overlap_aabb_os(saabb, oba, obb, rab, debug, dmax)
calculate overlap integrals (aa,bb)
Definition: generic_os_integrals.F:824
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_lri_inv_auto
integer, parameter, public do_lri_inv_auto
Definition: input_constants.F:1102
input_constants::do_lri_pseudoinv_svd
integer, parameter, public do_lri_pseudoinv_svd
Definition: input_constants.F:1102
input_constants::do_lri_inv
integer, parameter, public do_lri_inv
Definition: input_constants.F:1102
input_constants::do_lri_pseudoinv_diag
integer, parameter, public do_lri_pseudoinv_diag
Definition: input_constants.F:1102
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
lri_compression
integral compression (fix point accuracy)
Definition: lri_compression.F:14
lri_compression::lri_decomp_i
subroutine, public lri_decomp_i(aval, cont, ival)
...
Definition: lri_compression.F:135
lri_compression::lri_cont_mem
real(kind=dp) function, public lri_cont_mem(cont)
...
Definition: lri_compression.F:107
lri_compression::lri_comp
subroutine, public lri_comp(aval, amax, cont)
...
Definition: lri_compression.F:42
lri_environment_methods
Calculates integral matrices for LRIGPW method lri : local resolution of the identity.
Definition: lri_environment_methods.F:18
lri_environment_methods::v_int_ppl_update
subroutine, public v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
...
Definition: lri_environment_methods.F:1189
lri_environment_methods::calculate_lri_integrals
subroutine, public calculate_lri_integrals(lri_env, qs_env)
calculates integrals needed for the LRI density fitting, integrals are calculated once,...
Definition: lri_environment_methods.F:127
lri_environment_methods::lri_print_stat
subroutine, public lri_print_stat(qs_env, ltddfpt, tddfpt_lri_env)
...
Definition: lri_environment_methods.F:1248
lri_environment_methods::lri_kg_rho_update
subroutine, public lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist)
...
Definition: lri_environment_methods.F:340
lri_environment_methods::calculate_lri_densities
subroutine, public calculate_lri_densities(lri_env, lri_density, qs_env, pmatrix, cell_to_index, lri_rho_struct, atomic_kind_set, para_env, response_density)
performs the fitting of the density and distributes the fitted density on the grid
Definition: lri_environment_methods.F:569
lri_environment_methods::calculate_avec_lri
subroutine, public calculate_avec_lri(lri_env, lri_density, pmatrix, cell_to_index, response_density)
performs the fitting of the density; solves the linear system of equations; yield the expansion coeff...
Definition: lri_environment_methods.F:599
lri_environment_methods::build_lri_matrices
subroutine, public build_lri_matrices(lri_env, qs_env)
creates and initializes an lri_env
Definition: lri_environment_methods.F:105
lri_environment_methods::v_int_ppl_energy
subroutine, public v_int_ppl_energy(qs_env, lri_v_int, ecore_ppl_ri)
...
Definition: lri_environment_methods.F:1221
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition: lri_environment_types.F:18
lri_environment_types::allocate_lri_ints_rho
subroutine, public allocate_lri_ints_rho(lri_env, lri_ints_rho, nkind)
allocate lri_ints_rho, storing integral for the exact density
Definition: lri_environment_types.F:786
lri_environment_types::deallocate_lri_ppl_ints
subroutine, public deallocate_lri_ppl_ints(lri_ppl_ints)
deallocates the given lri_ppl_ints
Definition: lri_environment_types.F:1213
lri_environment_types::allocate_lri_ppl_ints
subroutine, public allocate_lri_ppl_ints(lri_env, lri_ppl_ints, atomic_kind_set)
allocate lri_ppl_ints, matrices that store LRI integrals
Definition: lri_environment_types.F:753
lri_environment_types::allocate_lri_coefs
subroutine, public allocate_lri_coefs(lri_env, lri_density, atomic_kind_set)
creates and initializes lri_coefs
Definition: lri_environment_types.F:974
lri_environment_types::lri_density_release
subroutine, public lri_density_release(lri_density)
releases the given lri_density
Definition: lri_environment_types.F:552
lri_environment_types::allocate_lri_ints
subroutine, public allocate_lri_ints(lri_env, lri_ints, nkind)
allocate lri_ints, matrices that store LRI integrals
Definition: lri_environment_types.F:567
lri_environment_types::deallocate_lri_ints
subroutine, public deallocate_lri_ints(lri_ints)
deallocates the given lri_ints
Definition: lri_environment_types.F:1106
lri_environment_types::deallocate_lri_ints_rho
subroutine, public deallocate_lri_ints_rho(lri_ints_rho)
deallocates the given lri_ints_rho
Definition: lri_environment_types.F:1237
lri_environment_types::lri_density_create
subroutine, public lri_density_create(lri_density)
creates and initializes an lri_density environment
Definition: lri_environment_types.F:537
lri_environment_types::allocate_lri_rhos
subroutine, public allocate_lri_rhos(lri_env, lri_rhos, nspin, nkind)
creates and initializes lri_rhos
Definition: lri_environment_types.F:862
lri_integrals
Calculates integrals for LRIGPW method lri : local resolution of the identity.
Definition: lri_integrals.F:19
lri_integrals::deallocate_int_type
subroutine, public deallocate_int_type(lriint, lridint)
...
Definition: lri_integrals.F:671
lri_integrals::allocate_int_type
subroutine, public allocate_int_type(lriint, lridint, nba, nbb, nfa, nfb, skip_sab, skip_soo, skip_aba, skip_abb, skip_dsab, skip_dsoo, skip_daba, skip_dabb)
...
Definition: lri_integrals.F:588
lri_integrals::lri_int2
subroutine, public lri_int2(lri_env, lrii, lriint, rab, obasa, obasb, fbasa, fbasb, iatom, jatom, ikind, jkind)
calcuates the lri integrals using solid harmonic Gaussians
Definition: lri_integrals.F:113
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::get_pseudo_inverse_svd
subroutine, public get_pseudo_inverse_svd(a, a_pinverse, rskip, determinant, sval)
returns the pseudoinverse of a real, square matrix using singular value decomposition
Definition: mathlib.F:949
mathlib::get_pseudo_inverse_diag
subroutine, public get_pseudo_inverse_diag(a, a_pinverse, rskip)
returns the pseudoinverse of a real, symmetric and positive definite matrix using diagonalization.
Definition: mathlib.F:1030
mathlib::invmat_symm
subroutine, public invmat_symm(a, cholesky_triangle)
returns inverse of real symmetric, positive definite matrix
Definition: mathlib.F:574
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
pw_types
Definition: pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition: qs_collocate_density.F:38
qs_collocate_density::calculate_lri_rho_elec
subroutine, public calculate_lri_rho_elec(lri_rho_g, lri_rho_r, qs_env, lri_coef, total_rho, basis_type, exact_1c_terms, pmat, atomlist)
Collocates the fitted lri density on a grid.
Definition: qs_collocate_density.F:517
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_force_types
Definition: qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition: qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition: qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition: qs_neighbor_list_types.F:549
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_set
subroutine, public qs_rho_set(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
...
Definition: qs_rho_types.F:308
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
virial_types
Definition: virial_types.F:13