(git:d18deda)
Loading...
Searching...
No Matches
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-2025 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! **************************************************************************************************
18MODULE 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_dbcsr_api, ONLY: dbcsr_create,&
27 dbcsr_get_block_p,&
28 dbcsr_p_type,&
29 dbcsr_release,&
30 dbcsr_replicate_all,&
31 dbcsr_type
32 USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
33 dbcsr_get_block_diag
34 USE cp_log_handling, ONLY: cp_get_default_logger,&
35 cp_logger_type
36 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
37 cp_print_key_unit_nr
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
97CONTAINS
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, potrf=.false.)
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
1355END MODULE lri_environment_methods
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
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:783
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_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:679
cp_dbcsr_api::dbcsr_replicate_all
subroutine, public dbcsr_replicate_all(matrix)
...
Definition cp_dbcsr_api.F:1133
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_contrib::dbcsr_get_block_diag
subroutine, public dbcsr_get_block_diag(matrix, diag)
Copies the diagonal blocks of matrix into diag.
Definition cp_dbcsr_contrib.F:450
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:853
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:1063
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:1124
input_constants::do_lri_pseudoinv_svd
integer, parameter, public do_lri_pseudoinv_svd
Definition input_constants.F:1124
input_constants::do_lri_inv
integer, parameter, public do_lri_inv
Definition input_constants.F:1124
input_constants::do_lri_pseudoinv_diag
integer, parameter, public do_lri_pseudoinv_diag
Definition input_constants.F:1124
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:788
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:1215
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:755
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:976
lri_environment_types::lri_density_release
subroutine, public lri_density_release(lri_density)
releases the given lri_density
Definition lri_environment_types.F:554
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:569
lri_environment_types::deallocate_lri_ints
subroutine, public deallocate_lri_ints(lri_ints)
deallocates the given lri_ints
Definition lri_environment_types.F:1108
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:1239
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:539
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:864
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:938
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:1019
mathlib::invmat_symm
subroutine, public invmat_symm(a, potrf, uplo)
returns inverse of real symmetric, positive definite matrix
Definition mathlib.F:580
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:518
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
lri_environment_types::lri_density_type
Definition lri_environment_types.F:359
lri_environment_types::lri_environment_type
Definition lri_environment_types.F:271
lri_environment_types::lri_int_rho_type
Definition lri_environment_types.F:155
lri_environment_types::lri_int_type
Definition lri_environment_types.F:92
lri_environment_types::lri_kind_type
Definition lri_environment_types.F:332
lri_environment_types::lri_list_type
Definition lri_environment_types.F:174
lri_environment_types::lri_rhoab_type
Definition lri_environment_types.F:59
lri_integrals::int_type
Definition lri_integrals.F:39
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
virial_types::virial_type
Definition virial_types.F:26