(git:f2099e5)
Loading...
Searching...
No Matches
iao_analysis.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculate intrinsic atomic orbitals and analyze wavefunctions
10!> \par History
11!> 03.2023 created [JGH]
12!> \author JGH
13! **************************************************************************************************
15 USE ai_contraction, ONLY: block_add,&
17 USE ai_overlap, ONLY: overlap_ab
27 USE bibliography, ONLY: knizia2013,&
28 cite_reference
29 USE cell_types, ONLY: cell_type,&
30 pbc
33 USE cp_dbcsr_api, ONLY: &
37 dbcsr_replicate_all, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
53 USE cp_fm_types, ONLY: cp_fm_create,&
66 USE iao_types, ONLY: iao_env_type,&
75 USE kinds, ONLY: default_path_length,&
77 dp
78 USE machine, ONLY: m_walltime
79 USE mathconstants, ONLY: twopi
80 USE mathlib, ONLY: invmat_symm,&
81 jacobi
85 USE orbital_pointers, ONLY: ncoset
90 USE pw_env_types, ONLY: pw_env_get,&
93 USE pw_types, ONLY: pw_c1d_gs_type,&
99 USE qs_kind_types, ONLY: get_qs_kind,&
101 USE qs_ks_types, ONLY: get_ks_env,&
108 USE qs_mo_types, ONLY: allocate_mo_set,&
110 get_mo_set,&
119 USE qs_subsys_types, ONLY: qs_subsys_get,&
121#include "./base/base_uses.f90"
122
123 IMPLICIT NONE
124 PRIVATE
125
126 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'iao_analysis'
127
129
131 MODULE PROCEDURE iao_calculate_dmat_diag, & ! diagonal occupation matrix
132 iao_calculate_dmat_full ! full occupation matrix
133 END INTERFACE
134
135! **************************************************************************************************
136
137CONTAINS
138
139! **************************************************************************************************
140!> \brief ...
141!> \param qs_env ...
142!> \param iao_env ...
143!> \param unit_nr ...
144!> \param c_iao_coef ...
145!> \param mos ...
146!> \param bond_centers ...
147! **************************************************************************************************
148 SUBROUTINE iao_wfn_analysis(qs_env, iao_env, unit_nr, c_iao_coef, mos, bond_centers)
149 TYPE(qs_environment_type), POINTER :: qs_env
150 TYPE(iao_env_type), INTENT(IN) :: iao_env
151 INTEGER, INTENT(IN) :: unit_nr
152 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
153 OPTIONAL :: c_iao_coef
154 TYPE(mo_set_type), DIMENSION(:), OPTIONAL, POINTER :: mos
155 REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL :: bond_centers
156
157 CHARACTER(len=*), PARAMETER :: routinen = 'iao_wfn_analysis'
158
159 CHARACTER(LEN=2) :: element_symbol
160 CHARACTER(LEN=default_string_length) :: bname
161 INTEGER :: handle, i, iatom, ikind, isgf, ispin, &
162 nao, natom, nimages, nkind, no, norb, &
163 nref, ns, nsgf, nspin, nvec, nx, order
164 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf, nsgfat
165 INTEGER, DIMENSION(2) :: nocc
166 LOGICAL :: cubes_iao, cubes_ibo, molden_iao, &
167 molden_ibo, uniform_occupation
168 REAL(kind=dp) :: fin, fout, t1, t2, trace, zval
169 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fdiag
170 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mcharge
171 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: moments
172 REAL(kind=dp), DIMENSION(:), POINTER :: occupation_numbers
173 TYPE(cell_type), POINTER :: cell
174 TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: smat_kind
175 TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:, :) :: oce_atom
176 TYPE(cp_fm_struct_type), POINTER :: fm_struct
177 TYPE(cp_fm_type) :: p_orb_ref, p_ref_orb, smo
178 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: c_loc_orb, ciao, cloc, cvec, iao_coef
179 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: zij_atom
180 TYPE(cp_fm_type), POINTER :: mo_coeff
181 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: sref, sro, sxo
182 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
183 TYPE(dbcsr_type) :: dmat
184 TYPE(dbcsr_type), POINTER :: smat
185 TYPE(dft_control_type), POINTER :: dft_control
186 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: oce_basis_set_list, orb_basis_set_list, &
187 ref_basis_set_list
188 TYPE(gto_basis_set_type), POINTER :: ocebasis, orbbasis, refbasis
189 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mo_iao, mo_loc
190 TYPE(mo_set_type), DIMENSION(:), POINTER :: my_mos
191 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
192 POINTER :: sro_list, srr_list, sxo_list
193 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
194 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
195 TYPE(qs_kind_type), POINTER :: qs_kind
196 TYPE(qs_ks_env_type), POINTER :: ks_env
197 TYPE(section_vals_type), POINTER :: iao_cubes_section, iao_molden_section, &
198 ibo_cc_section, ibo_cubes_section, &
199 ibo_molden_section
200
201 ! only do IAO analysis if explicitly requested
202 cpassert(iao_env%do_iao)
203
204 ! k-points?
205 NULLIFY (particle_set)
206 CALL get_qs_env(qs_env, dft_control=dft_control, particle_set=particle_set)
207 nspin = dft_control%nspins
208 nimages = dft_control%nimages
209 IF (nimages > 1) THEN
210 IF (unit_nr > 0) THEN
211 WRITE (unit=unit_nr, fmt="(T2,A)") &
212 "K-Points: Intrinsic Atomic Orbitals Analysis not available."
213 END IF
214 END IF
215 IF (nimages > 1) RETURN
216
217 CALL timeset(routinen, handle)
218
219 IF (unit_nr > 0) THEN
220 WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
221 WRITE (unit=unit_nr, fmt="(T24,A)") "INTRINSIC ATOMIC ORBITALS ANALYSIS"
222 WRITE (unit=unit_nr, fmt="(T13,A)") "G. Knizia, J. Chem. Theory Comput. 9, 4834-4843 (2013)"
223 WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!'
224 END IF
225 CALL cite_reference(knizia2013)
226
227 ! input sections
228 iao_molden_section => iao_env%iao_molden_section
229 iao_cubes_section => iao_env%iao_cubes_section
230 ibo_molden_section => iao_env%ibo_molden_section
231 ibo_cubes_section => iao_env%ibo_cubes_section
232 ibo_cc_section => iao_env%ibo_cc_section
233
234 !
235 molden_iao = .false.
236 IF (ASSOCIATED(iao_molden_section)) THEN
237 CALL section_vals_get(iao_molden_section, explicit=molden_iao)
238 END IF
239 cubes_iao = .false.
240 IF (ASSOCIATED(iao_cubes_section)) THEN
241 CALL section_vals_get(iao_cubes_section, explicit=cubes_iao)
242 END IF
243 molden_ibo = .false.
244 IF (ASSOCIATED(ibo_molden_section)) THEN
245 CALL section_vals_get(ibo_molden_section, explicit=molden_ibo)
246 END IF
247 cubes_ibo = .false.
248 IF (ASSOCIATED(ibo_cubes_section)) THEN
249 CALL section_vals_get(ibo_cubes_section, explicit=cubes_ibo)
250 END IF
251
252 ! check for or generate reference basis
253 CALL create_minbas_set(qs_env, unit_nr)
254
255 ! overlap matrices
256 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s)
257 smat => matrix_s(1, 1)%matrix
258
259 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
260 nkind = SIZE(qs_kind_set)
261 ALLOCATE (ref_basis_set_list(nkind), orb_basis_set_list(nkind))
262 DO ikind = 1, nkind
263 qs_kind => qs_kind_set(ikind)
264 NULLIFY (ref_basis_set_list(ikind)%gto_basis_set)
265 NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
266 NULLIFY (refbasis, orbbasis)
267 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orbbasis, basis_type="ORB")
268 IF (ASSOCIATED(orbbasis)) orb_basis_set_list(ikind)%gto_basis_set => orbbasis
269 CALL get_qs_kind(qs_kind=qs_kind, basis_set=refbasis, basis_type="MIN")
270 IF (ASSOCIATED(refbasis)) ref_basis_set_list(ikind)%gto_basis_set => refbasis
271 END DO
272
273 ! neighbor lists
274 NULLIFY (srr_list, sro_list)
275 CALL setup_neighbor_list(srr_list, ref_basis_set_list, qs_env=qs_env)
276 CALL setup_neighbor_list(sro_list, ref_basis_set_list, orb_basis_set_list, qs_env=qs_env)
277
278 ! Srr and Sro overlap matrices
279 NULLIFY (sref, sro)
280 CALL get_qs_env(qs_env, ks_env=ks_env)
281 CALL build_overlap_matrix_simple(ks_env, sref, &
282 ref_basis_set_list, ref_basis_set_list, srr_list)
283 CALL build_overlap_matrix_simple(ks_env, sro, &
284 ref_basis_set_list, orb_basis_set_list, sro_list)
285
286 IF (PRESENT(mos)) THEN
287 my_mos => mos
288 ELSE
289 CALL get_qs_env(qs_env, mos=my_mos)
290 END IF
291 CALL get_mo_set(my_mos(1), mo_coeff=mo_coeff)
292
293 IF (iao_env%do_fragments) CALL organise_printout_for_rose(qs_env, orb_basis_set_list)
294
295 IF (PRESENT(mos)) THEN
296 my_mos => mos
297 ELSE
298 CALL get_qs_env(qs_env, mos=my_mos)
299 END IF
300 CALL get_mo_set(my_mos(1), mo_coeff=mo_coeff)
301 CALL dbcsr_get_info(sro(1)%matrix, nfullrows_total=nref, nfullcols_total=nao)
302
303 ! Projectors
304 NULLIFY (fm_struct)
305 CALL cp_fm_struct_create(fm_struct, context=mo_coeff%matrix_struct%context, nrow_global=nref, &
306 ncol_global=nao, para_env=mo_coeff%matrix_struct%para_env)
307 CALL cp_fm_create(p_ref_orb, fm_struct)
308 CALL cp_fm_struct_release(fm_struct)
309 NULLIFY (fm_struct)
310 CALL cp_fm_struct_create(fm_struct, context=mo_coeff%matrix_struct%context, nrow_global=nao, &
311 ncol_global=nref, para_env=mo_coeff%matrix_struct%para_env)
312 CALL cp_fm_create(p_orb_ref, fm_struct)
313 CALL cp_fm_struct_release(fm_struct)
314 !
315 CALL iao_projectors(smat, sref(1)%matrix, sro(1)%matrix, p_orb_ref, p_ref_orb, iao_env%eps_svd)
316
317 ! occupied orbitals, orbitals considered for IAO generation
318 nocc(1:2) = 0
319 DO ispin = 1, nspin
320 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nmo=norb, uniform_occupation=uniform_occupation, &
321 occupation_numbers=occupation_numbers)
322 IF (uniform_occupation) THEN
323 nvec = norb
324 ELSE
325 nvec = 0
326 DO i = 1, norb
327 IF (occupation_numbers(i) > iao_env%eps_occ) THEN
328 nvec = i
329 ELSE
330 EXIT
331 END IF
332 END DO
333 END IF
334 nocc(ispin) = nvec
335 END DO
336 ! generate IAOs
337 ALLOCATE (iao_coef(nspin), cvec(nspin))
338 DO ispin = 1, nspin
339 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff)
340 nvec = nocc(ispin)
341 IF (nvec > 0) THEN
342 NULLIFY (fm_struct)
343 CALL cp_fm_struct_create(fm_struct, ncol_global=nvec, template_fmstruct=mo_coeff%matrix_struct)
344 CALL cp_fm_create(cvec(ispin), fm_struct)
345 CALL cp_fm_struct_release(fm_struct)
346 CALL cp_fm_to_fm(mo_coeff, cvec(ispin), nvec)
347 !
348 NULLIFY (fm_struct)
349 CALL cp_fm_struct_create(fm_struct, ncol_global=nref, template_fmstruct=mo_coeff%matrix_struct)
350 CALL cp_fm_create(iao_coef(ispin), fm_struct)
351 CALL cp_fm_struct_release(fm_struct)
352 !
353 CALL intrinsic_ao_calc(smat, p_orb_ref, p_ref_orb, cvec(ispin), iao_coef(ispin))
354 END IF
355 END DO
356 !
357 IF (iao_env%do_charges) THEN
358 ! MOs in IAO basis
359 ALLOCATE (ciao(nspin))
360 DO ispin = 1, nspin
361 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=norb)
362 CALL cp_fm_create(smo, mo_coeff%matrix_struct)
363 NULLIFY (fm_struct)
364 CALL cp_fm_struct_create(fm_struct, nrow_global=nref, template_fmstruct=mo_coeff%matrix_struct)
365 CALL cp_fm_create(ciao(ispin), fm_struct)
366 CALL cp_fm_struct_release(fm_struct)
367 ! CIAO = A(T)*S*C
368 CALL cp_dbcsr_sm_fm_multiply(smat, mo_coeff, smo, ncol=norb)
369 CALL parallel_gemm('T', 'N', nref, norb, nao, 1.0_dp, iao_coef(ispin), smo, 0.0_dp, ciao(ispin))
370 CALL cp_fm_release(smo)
371 END DO
372 !
373 ! population analysis
374 IF (unit_nr > 0) THEN
375 WRITE (unit_nr, '(/,T2,A)') 'Intrinsic AO Population Analysis '
376 END IF
377 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
378 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
379 ALLOCATE (mcharge(natom, nspin))
380 CALL dbcsr_get_info(sref(1)%matrix, nfullrows_total=nref)
381 DO ispin = 1, nspin
382 CALL get_mo_set(my_mos(ispin), &
383 uniform_occupation=uniform_occupation, &
384 occupation_numbers=occupation_numbers)
385 ! diagonal block matrix
386 CALL dbcsr_create(dmat, template=sref(1)%matrix)
388 !
389 CALL iao_calculate_dmat(ciao(ispin), dmat, occupation_numbers, uniform_occupation)
390 CALL dbcsr_trace(dmat, trace)
391 IF (unit_nr > 0) THEN
392 WRITE (unit_nr, '(T2,A,I2,T66,F15.4)') 'Number of Electrons: Trace(Piao) Spin ', ispin, trace
393 END IF
394 CALL iao_charges(dmat, mcharge(:, ispin))
395 CALL dbcsr_release(dmat)
396 END DO
397 CALL print_atomic_charges(particle_set, qs_kind_set, unit_nr, "Intrinsic Atomic Orbital Charges", &
398 electronic_charges=mcharge)
399 DEALLOCATE (mcharge)
400 CALL cp_fm_release(ciao)
401 END IF
402 !
403 ! Deallocate the neighbor list structure
404 CALL release_neighbor_list_sets(srr_list)
405 CALL release_neighbor_list_sets(sro_list)
406 IF (ASSOCIATED(sref)) CALL dbcsr_deallocate_matrix_set(sref)
407 IF (ASSOCIATED(sro)) CALL dbcsr_deallocate_matrix_set(sro)
408 CALL cp_fm_release(p_ref_orb)
409 CALL cp_fm_release(p_orb_ref)
410 CALL cp_fm_release(cvec)
411 DEALLOCATE (orb_basis_set_list)
412 !
413 IF (iao_env%do_oce) THEN
414 ! generate OCE basis
415 IF (unit_nr > 0) THEN
416 WRITE (unit_nr, '(T2,A)') "IAO One-Center Expansion: OCE Basis Set"
417 END IF
418 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
419 nkind = SIZE(qs_kind_set)
420 ALLOCATE (oce_basis_set_list(nkind), orb_basis_set_list(nkind))
421 DO ikind = 1, nkind
422 qs_kind => qs_kind_set(ikind)
423 NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
424 NULLIFY (orbbasis)
425 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orbbasis, basis_type="ORB")
426 IF (ASSOCIATED(orbbasis)) orb_basis_set_list(ikind)%gto_basis_set => orbbasis
427 END DO
428 DO ikind = 1, nkind
429 orbbasis => orb_basis_set_list(ikind)%gto_basis_set
430 cpassert(ASSOCIATED(orbbasis))
431 NULLIFY (ocebasis)
432 CALL create_oce_basis(ocebasis, orbbasis, iao_env%lmax_oce, iao_env%nbas_oce)
433 CALL init_orb_basis_set(ocebasis)
434 CALL init_interaction_radii_orb_basis(ocebasis, dft_control%qs_control%eps_pgf_orb)
435 oce_basis_set_list(ikind)%gto_basis_set => ocebasis
436 IF (unit_nr > 0) THEN
437 qs_kind => qs_kind_set(ikind)
438 CALL get_qs_kind(qs_kind, zeff=zval, element_symbol=element_symbol)
439 CALL get_gto_basis_set(ocebasis, name=bname, nsgf=nsgf)
440 WRITE (unit_nr, '(T2,A,A,T14,A,I4,T40,A,A30)') "Kind: ", element_symbol, "NBasFun: ", nsgf, &
441 "OCE Basis: ", adjustl(trim(bname))
442 END IF
443 END DO
444 IF (unit_nr > 0) WRITE (unit_nr, *)
445 ! OCE basis overlap
446 ALLOCATE (smat_kind(nkind))
447 DO ikind = 1, nkind
448 ocebasis => oce_basis_set_list(ikind)%gto_basis_set
449 CALL get_gto_basis_set(gto_basis_set=ocebasis, nsgf=nsgf)
450 ALLOCATE (smat_kind(ikind)%array(nsgf, nsgf))
451 END DO
452 CALL oce_overlap_matrix(smat_kind, oce_basis_set_list)
453 ! overlap between IAO OCE basis and orbital basis
454 NULLIFY (sxo, sxo_list)
455 CALL setup_neighbor_list(sxo_list, oce_basis_set_list, orb_basis_set_list, qs_env=qs_env)
456 CALL get_qs_env(qs_env, ks_env=ks_env)
457 CALL build_overlap_matrix_simple(ks_env, sxo, &
458 oce_basis_set_list, orb_basis_set_list, sxo_list)
459 !
460 ! One Center Expansion of IAOs
461 CALL get_qs_env(qs_env=qs_env, natom=natom)
462 ALLOCATE (oce_atom(natom, nspin))
463 CALL iao_oce_expansion(qs_env, oce_basis_set_list, smat_kind, sxo, iao_coef, oce_atom, unit_nr)
464 !
465 DO ikind = 1, nkind
466 ocebasis => oce_basis_set_list(ikind)%gto_basis_set
467 CALL deallocate_gto_basis_set(ocebasis)
468 END DO
469 DEALLOCATE (oce_basis_set_list, orb_basis_set_list)
470 !
471 CALL release_neighbor_list_sets(sxo_list)
472 IF (ASSOCIATED(sxo)) CALL dbcsr_deallocate_matrix_set(sxo)
473 DO ikind = 1, nkind
474 DEALLOCATE (smat_kind(ikind)%array)
475 END DO
476 DEALLOCATE (smat_kind)
477 DO ispin = 1, nspin
478 DO iatom = 1, natom
479 DEALLOCATE (oce_atom(iatom, ispin)%array)
480 END DO
481 END DO
482 DEALLOCATE (oce_atom)
483 END IF
484 !
485 IF (molden_iao) THEN
486 ! Print basis functions: molden file
487 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
488 CALL get_qs_env(qs_env, cell=cell)
489 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
490 IF (unit_nr > 0) THEN
491 WRITE (unit_nr, '(T2,A)') ' Write IAO in MOLDEN Format'
492 END IF
493 ALLOCATE (mo_iao(nspin))
494 DO ispin = 1, nspin
495 CALL get_mo_set(my_mos(ispin), nao=nao)
496 CALL allocate_mo_set(mo_iao(ispin), nao, nref, nref, 0.0_dp, 1.0_dp, 0.0_dp)
497 CALL init_mo_set(mo_iao(ispin), fm_ref=iao_coef(ispin), name="iao_set")
498 CALL cp_fm_to_fm(iao_coef(ispin), mo_iao(ispin)%mo_coeff)
499 CALL set_mo_set(mo_iao(ispin), homo=nref)
500 mo_iao(ispin)%occupation_numbers = 1.0_dp
501 END DO
502 ! Print IAO basis functions: molden format
503 CALL write_mos_molden(mo_iao, qs_kind_set, particle_set, iao_molden_section, cell=cell, qs_env=qs_env)
504 DO ispin = 1, nspin
505 CALL deallocate_mo_set(mo_iao(ispin))
506 END DO
507 DEALLOCATE (mo_iao)
508 END IF
509 IF (cubes_iao) THEN
510 IF (unit_nr > 0) THEN
511 WRITE (unit_nr, '(T2,A)') ' Write IAO as CUBE Files'
512 END IF
513 ! Print basis functions: cube file
514 CALL print_iao_cubes(qs_env, iao_cubes_section, iao_coef, ref_basis_set_list)
515 END IF
516 !
517 ! Intrinsic bond orbitals
518 IF (iao_env%do_bondorbitals) THEN
519 IF (unit_nr > 0) THEN
520 WRITE (unit_nr, '(/,T2,A)') 'Intrinsic Bond Orbital Generation'
521 END IF
522 ! MOs in IAO basis -> ciao
523 ALLOCATE (ciao(nspin))
524 DO ispin = 1, nspin
525 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=norb)
526 CALL cp_fm_create(smo, mo_coeff%matrix_struct)
527 NULLIFY (fm_struct)
528 CALL cp_fm_struct_create(fm_struct, nrow_global=nref, template_fmstruct=mo_coeff%matrix_struct)
529 CALL cp_fm_create(ciao(ispin), fm_struct)
530 CALL cp_fm_struct_release(fm_struct)
531 ! CIAO = A(T)*S*C
532 CALL cp_dbcsr_sm_fm_multiply(smat, mo_coeff, smo, ncol=norb)
533 CALL parallel_gemm('T', 'N', nref, norb, nao, 1.0_dp, iao_coef(ispin), smo, 0.0_dp, ciao(ispin))
534 CALL cp_fm_release(smo)
535 END DO
536 !
537 ! localize occupied orbitals nocc(ispin), see IAO generation
538 !
539 ALLOCATE (cloc(nspin), c_loc_orb(nspin))
540 DO ispin = 1, nspin
541 NULLIFY (fm_struct)
542 CALL cp_fm_struct_create(fm_struct, ncol_global=nocc(ispin), &
543 template_fmstruct=ciao(ispin)%matrix_struct)
544 CALL cp_fm_create(cloc(ispin), fm_struct)
545 CALL cp_fm_struct_release(fm_struct)
546 CALL cp_fm_to_fm(ciao(ispin), cloc(ispin), ncol=nocc(ispin))
547 END DO
548 CALL cp_fm_release(ciao)
549 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
550 ALLOCATE (first_sgf(natom), last_sgf(natom), nsgfat(natom))
551 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
552 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf, &
553 nsgf=nsgfat, basis=ref_basis_set_list)
554
555 IF (iao_env%loc_operator /= do_iaoloc_pm2) THEN
556 cpabort("IAO localization operator NYA")
557 END IF
558 IF (iao_env%eloc_function /= do_iaoloc_enone) THEN
559 cpabort("IAO energy weight localization NYA")
560 END IF
561
562 DO ispin = 1, nspin
563 nvec = nocc(ispin)
564 IF (nvec > 0) THEN
565 ALLOCATE (zij_atom(natom, 1), fdiag(nvec))
566 NULLIFY (fm_struct)
567 CALL cp_fm_struct_create(fm_struct, nrow_global=nvec, ncol_global=nvec, &
568 template_fmstruct=cloc(ispin)%matrix_struct)
569 DO iatom = 1, natom
570 CALL cp_fm_create(zij_atom(iatom, 1), fm_struct)
571 isgf = first_sgf(iatom)
572 CALL parallel_gemm('T', 'N', nvec, nvec, nsgfat(iatom), 1.0_dp, cloc(ispin), cloc(ispin), &
573 0.0_dp, zij_atom(iatom, 1), &
574 a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf)
575 END DO
576 CALL cp_fm_struct_release(fm_struct)
577 !
578 t1 = m_walltime()
579 order = 4
580 fin = 0.0_dp
581 DO iatom = 1, natom
582 CALL cp_fm_get_diag(zij_atom(iatom, 1), fdiag)
583 fin = fin + sum(fdiag**order)
584 END DO
585 fin = fin**(1._dp/order)
586 ! QR localization
587 CALL scdm_qrfact(cloc(ispin))
588 !
589 DO iatom = 1, natom
590 isgf = first_sgf(iatom)
591 CALL parallel_gemm('T', 'N', nvec, nvec, nsgfat(iatom), 1.0_dp, cloc(ispin), cloc(ispin), &
592 0.0_dp, zij_atom(iatom, 1), &
593 a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf)
594 END DO
595 fout = 0.0_dp
596 DO iatom = 1, natom
597 CALL cp_fm_get_diag(zij_atom(iatom, 1), fdiag)
598 fout = fout + sum(fdiag**order)
599 END DO
600 fout = fout**(1._dp/order)
601 DEALLOCATE (fdiag)
602 t2 = m_walltime()
603 IF (unit_nr > 0) THEN
604 WRITE (unit_nr, '(T2,A,F14.8,A,F14.8,A,F8.2)') &
605 'SCDM pre-localization: fin=', fin, " fout=", fout, " Time=", t2 - t1
606 END IF
607 !
608 CALL pm_localization(zij_atom, cloc(ispin), order, 1.e-12_dp, unit_nr)
609 !
610 CALL cp_fm_release(zij_atom)
611 CALL get_mo_set(my_mos(ispin), nao=nao)
612 NULLIFY (fm_struct)
613 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, &
614 template_fmstruct=cloc(ispin)%matrix_struct)
615 CALL cp_fm_create(c_loc_orb(ispin), fm_struct)
616 CALL cp_fm_struct_release(fm_struct)
617 CALL parallel_gemm('N', 'N', nao, nvec, nref, 1.0_dp, iao_coef(ispin), cloc(ispin), &
618 0.0_dp, c_loc_orb(ispin))
619 END IF
620 END DO
621 DEALLOCATE (first_sgf, last_sgf, nsgfat)
622 CALL cp_fm_release(cloc)
623 !
624 IF (iao_env%do_center) THEN
625 IF (unit_nr > 0) THEN
626 WRITE (unit_nr, '(T2,A)') ' Calculate Localized Orbital Centers and Spread'
627 IF (iao_env%pos_periodic) THEN
628 WRITE (unit_nr, '(T2,A)') ' Use Berry Phase Position Operator'
629 ELSE
630 WRITE (unit_nr, '(T2,A)') ' Use Local Position Operator'
631 END IF
632 END IF
633 nvec = maxval(nocc)
634 ! x y z m2(1) m2(2)
635 ALLOCATE (moments(5, nvec, nspin))
636 moments = 0.0_dp
637 IF (iao_env%pos_periodic) THEN
638 CALL center_spread_berry(qs_env, c_loc_orb, moments)
639 ELSE
640 CALL center_spread_loc(qs_env, c_loc_orb, moments)
641 END IF
642 !
643 IF (ASSOCIATED(ibo_cc_section)) THEN
644 CALL print_center_spread(moments, nocc, ibo_cc_section)
645 END IF
646
647 IF (iao_env%do_fragments) CALL print_fragment_association(qs_env, moments, unit_nr)
648
649 IF (PRESENT(bond_centers)) THEN
650 nx = SIZE(bond_centers, 1)
651 no = SIZE(bond_centers, 2)
652 ns = SIZE(bond_centers, 3)
653 cpassert(no >= nvec)
654 cpassert(ns == nspin)
655 cpassert(nx >= 3)
656 nx = min(nx, 5)
657 bond_centers(1:nx, 1:nvec, 1:ns) = moments(1:nx, 1:nvec, 1:ns)
658 END IF
659 DEALLOCATE (moments)
660 END IF
661 !
662 IF (molden_ibo) THEN
663 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
664 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
665 IF (unit_nr > 0) THEN
666 WRITE (unit_nr, '(T2,A)') ' Write IBO in MOLDEN Format'
667 END IF
668 ALLOCATE (mo_loc(nspin))
669 DO ispin = 1, nspin
670 CALL get_mo_set(my_mos(ispin), nao=nao)
671 nvec = nocc(ispin)
672 CALL allocate_mo_set(mo_loc(ispin), nao, nvec, nvec, 0.0_dp, 1.0_dp, 0.0_dp)
673 CALL init_mo_set(mo_loc(ispin), fm_ref=c_loc_orb(ispin), name="ibo_orb")
674 CALL cp_fm_to_fm(c_loc_orb(ispin), mo_loc(ispin)%mo_coeff)
675 CALL set_mo_set(mo_loc(ispin), homo=nvec)
676 mo_loc(ispin)%occupation_numbers = 1.0_dp
677 END DO
678 ! Print IBO functions: molden format
679 CALL write_mos_molden(mo_loc, qs_kind_set, particle_set, ibo_molden_section, cell=cell, qs_env=qs_env)
680 DO ispin = 1, nspin
681 CALL deallocate_mo_set(mo_loc(ispin))
682 END DO
683 DEALLOCATE (mo_loc)
684 END IF
685 !
686 IF (cubes_ibo) THEN
687 IF (unit_nr > 0) THEN
688 WRITE (unit_nr, '(T2,A)') ' Write IBO on CUBE files'
689 END IF
690 ! Print localized orbital functions: cube file
691 CALL print_ibo_cubes(qs_env, ibo_cubes_section, c_loc_orb)
692 END IF
693 !
694 IF (PRESENT(mos) .AND. ALLOCATED(c_loc_orb)) THEN
695 ! c_loc_orb -> mos
696 DO ispin = 1, nspin
697 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
698 nvec = nocc(ispin)
699 CALL cp_fm_to_fm(c_loc_orb(ispin), mo_coeff, ncol=nvec)
700 END DO
701 END IF
702 CALL cp_fm_release(c_loc_orb)
703 END IF
704 DEALLOCATE (ref_basis_set_list)
705
706 IF (PRESENT(c_iao_coef)) THEN
707 CALL cp_fm_release(c_iao_coef)
708 ALLOCATE (c_iao_coef(nspin))
709 DO ispin = 1, nspin
710 CALL cp_fm_create(c_iao_coef(ispin), iao_coef(ispin)%matrix_struct)
711 CALL cp_fm_to_fm(iao_coef(ispin), c_iao_coef(ispin))
712 END DO
713 END IF
714 CALL cp_fm_release(iao_coef)
715
716 IF (unit_nr > 0) THEN
717 WRITE (unit_nr, '(T2,A)') &
718 '!----------------------------END OF IAO ANALYSIS------------------------------!'
719 END IF
720
721 CALL timestop(handle)
722
723 END SUBROUTINE iao_wfn_analysis
724
725! **************************************************************************************************
726!> \brief Computes projector matrices for ref to orb basis and reverse
727!> \param smat ...
728!> \param sref ...
729!> \param s_r_o ...
730!> \param p_o_r ...
731!> \param p_r_o ...
732!> \param eps_svd ...
733! **************************************************************************************************
734 SUBROUTINE iao_projectors(smat, sref, s_r_o, p_o_r, p_r_o, eps_svd)
735 TYPE(dbcsr_type), INTENT(INOUT) :: smat, sref, s_r_o
736 TYPE(cp_fm_type), INTENT(INOUT) :: p_o_r, p_r_o
737 REAL(kind=dp), INTENT(IN) :: eps_svd
738
739 CHARACTER(len=*), PARAMETER :: routinen = 'iao_projectors'
740
741 INTEGER :: handle, norb, nref
742 TYPE(cp_fm_struct_type), POINTER :: fm_struct
743 TYPE(cp_fm_type) :: fm_inv, fm_s, fm_sro
744
745 CALL timeset(routinen, handle)
746
747 CALL cp_fm_get_info(p_r_o, nrow_global=nref, ncol_global=norb)
748 CALL cp_fm_create(fm_sro, p_r_o%matrix_struct)
749 CALL copy_dbcsr_to_fm(s_r_o, fm_sro)
750
751 NULLIFY (fm_struct)
752 CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=norb, &
753 template_fmstruct=p_r_o%matrix_struct)
754 CALL cp_fm_create(fm_s, fm_struct, name="smat")
755 CALL cp_fm_create(fm_inv, fm_struct, name="sinv")
756 CALL cp_fm_struct_release(fm_struct)
757 CALL copy_dbcsr_to_fm(smat, fm_s)
758 CALL cp_fm_invert(fm_s, fm_inv, eps_svd=eps_svd)
759 CALL parallel_gemm('N', 'T', norb, nref, norb, 1.0_dp, fm_inv, fm_sro, 0.0_dp, p_o_r)
760 CALL cp_fm_release(fm_s)
761 CALL cp_fm_release(fm_inv)
762
763 NULLIFY (fm_struct)
764 CALL cp_fm_struct_create(fm_struct, nrow_global=nref, ncol_global=nref, &
765 template_fmstruct=p_r_o%matrix_struct)
766 CALL cp_fm_create(fm_s, fm_struct, name="sref")
767 CALL cp_fm_create(fm_inv, fm_struct, name="sinv")
768 CALL cp_fm_struct_release(fm_struct)
769 CALL copy_dbcsr_to_fm(sref, fm_s)
770 CALL cp_fm_invert(fm_s, fm_inv, eps_svd=eps_svd)
771 CALL parallel_gemm('N', 'N', nref, norb, nref, 1.0_dp, fm_inv, fm_sro, 0.0_dp, p_r_o)
772 CALL cp_fm_release(fm_s)
773 CALL cp_fm_release(fm_inv)
774
775 CALL cp_fm_release(fm_sro)
776
777 CALL timestop(handle)
778
779 END SUBROUTINE iao_projectors
780
781! **************************************************************************************************
782!> \brief Computes intrinsic orbitals for a given MO vector set
783!> \param smat ...
784!> \param p_o_r ...
785!> \param p_r_o ...
786!> \param cvec ...
787!> \param avec ...
788! **************************************************************************************************
789 SUBROUTINE intrinsic_ao_calc(smat, p_o_r, p_r_o, cvec, avec)
790 TYPE(dbcsr_type), INTENT(INOUT) :: smat
791 TYPE(cp_fm_type), INTENT(INOUT) :: p_o_r, p_r_o, cvec, avec
792
793 CHARACTER(len=*), PARAMETER :: routinen = 'intrinsic_ao_calc'
794
795 INTEGER :: handle, nao, norb, nref
796 TYPE(cp_fm_struct_type), POINTER :: fm_struct
797 TYPE(cp_fm_type) :: ctvec, pc, sc, sct, vec1
798
799 CALL timeset(routinen, handle)
800
801 ! number of orbitals
802 CALL cp_fm_get_info(cvec, ncol_global=norb)
803 ! basis set sizes
804 CALL cp_fm_get_info(p_r_o, nrow_global=nref, ncol_global=nao)
805 ! temp array
806 NULLIFY (fm_struct)
807 CALL cp_fm_struct_create(fm_struct, nrow_global=nref, ncol_global=norb, &
808 template_fmstruct=cvec%matrix_struct)
809 CALL cp_fm_create(pc, fm_struct)
810 CALL cp_fm_struct_release(fm_struct)
811 ! CT = orth(Por.Pro.C)
812 CALL parallel_gemm('N', 'N', nref, norb, nao, 1.0_dp, p_r_o, cvec, 0.0_dp, pc)
813 CALL cp_fm_create(ctvec, cvec%matrix_struct)
814 CALL parallel_gemm('N', 'N', nao, norb, nref, 1.0_dp, p_o_r, pc, 0.0_dp, ctvec)
815 CALL cp_fm_release(pc)
816 CALL make_basis_lowdin(ctvec, norb, smat)
817 ! S*C and S*CT
818 CALL cp_fm_create(sc, cvec%matrix_struct)
819 CALL cp_dbcsr_sm_fm_multiply(smat, cvec, sc, ncol=norb)
820 CALL cp_fm_create(sct, cvec%matrix_struct)
821 CALL cp_dbcsr_sm_fm_multiply(smat, ctvec, sct, ncol=norb)
822 CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=nref, &
823 template_fmstruct=cvec%matrix_struct)
824 CALL cp_fm_create(pc, fm_struct)
825 CALL cp_fm_struct_release(fm_struct)
826 ! V1 = (CT*SCT(T))Por
827 CALL parallel_gemm('T', 'N', norb, nref, nao, 1.0_dp, sct, p_o_r, 0.0_dp, pc)
828 NULLIFY (fm_struct)
829 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nref, &
830 template_fmstruct=cvec%matrix_struct)
831 CALL cp_fm_create(vec1, fm_struct)
832 CALL cp_fm_struct_release(fm_struct)
833 CALL parallel_gemm('N', 'N', nao, nref, norb, 1.0_dp, ctvec, pc, 0.0_dp, vec1)
834 ! A = C*SC(T)*V1
835 CALL parallel_gemm('T', 'N', norb, nref, nao, 1.0_dp, sc, vec1, 0.0_dp, pc)
836 CALL parallel_gemm('N', 'N', nao, nref, norb, 1.0_dp, cvec, pc, 0.0_dp, avec)
837 ! V1 = Por - V1
838 CALL cp_fm_geadd(1.0_dp, 'N', p_o_r, -1.0_dp, vec1)
839 ! A = A + V1
840 CALL cp_fm_geadd(1.0_dp, 'N', vec1, 1.0_dp, avec)
841 ! A = A - C*SC(T)*V1
842 CALL parallel_gemm('T', 'N', norb, nref, nao, 1.0_dp, sc, vec1, 0.0_dp, pc)
843 CALL parallel_gemm('N', 'N', nao, nref, norb, -1.0_dp, cvec, pc, 1.0_dp, avec)
844 ! A = orth(A)
845 CALL make_basis_lowdin(avec, nref, smat)
846
847 ! clean up
848 CALL cp_fm_release(pc)
849 CALL cp_fm_release(vec1)
850 CALL cp_fm_release(sc)
851 CALL cp_fm_release(sct)
852 CALL cp_fm_release(ctvec)
853
854 CALL timestop(handle)
855
856 END SUBROUTINE intrinsic_ao_calc
857
858! **************************************************************************************************
859!> \brief Calculate the density matrix from fm vectors using occupation numbers
860!> \param cvec ...
861!> \param density_matrix ...
862!> \param occupation ...
863!> \param uniform_occupation ...
864! **************************************************************************************************
865 SUBROUTINE iao_calculate_dmat_diag(cvec, density_matrix, occupation, uniform_occupation)
866
867 TYPE(cp_fm_type), INTENT(IN) :: cvec
868 TYPE(dbcsr_type) :: density_matrix
869 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
870 LOGICAL, INTENT(IN) :: uniform_occupation
871
872 CHARACTER(len=*), PARAMETER :: routineN = 'iao_calculate_dmat_diag'
873
874 INTEGER :: handle, ncol
875 REAL(KIND=dp) :: alpha
876 TYPE(cp_fm_type) :: fm_tmp
877
878 CALL timeset(routinen, handle)
879
880 CALL dbcsr_set(density_matrix, 0.0_dp)
881
882 CALL cp_fm_get_info(cvec, ncol_global=ncol)
883 IF (.NOT. uniform_occupation) THEN
884 CALL cp_fm_create(fm_tmp, cvec%matrix_struct)
885 CALL cp_fm_to_fm(cvec, fm_tmp)
886 CALL cp_fm_column_scale(fm_tmp, occupation(1:ncol))
887 alpha = 1.0_dp
888 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
889 matrix_v=cvec, matrix_g=fm_tmp, &
890 ncol=ncol, alpha=alpha)
891 CALL cp_fm_release(fm_tmp)
892 ELSE
893 alpha = occupation(1)
894 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
895 matrix_v=cvec, ncol=ncol, alpha=alpha)
896 END IF
897
898 CALL timestop(handle)
899
900 END SUBROUTINE iao_calculate_dmat_diag
901
902! **************************************************************************************************
903!> \brief Calculate the density matrix from fm vectors using an occupation matrix
904!> \param cvec ...
905!> \param density_matrix ...
906!> \param weight ...
907!> \param occmat ...
908! **************************************************************************************************
909 SUBROUTINE iao_calculate_dmat_full(cvec, density_matrix, weight, occmat)
910
911 TYPE(cp_fm_type), INTENT(IN) :: cvec
912 TYPE(dbcsr_type) :: density_matrix
913 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: weight
914 TYPE(cp_fm_type), INTENT(IN) :: occmat
915
916 CHARACTER(len=*), PARAMETER :: routinen = 'iao_calculate_dmat_full'
917
918 INTEGER :: handle, ic, jc, ncol
919 REAL(kind=dp) :: alpha
920 TYPE(cp_fm_struct_type), POINTER :: fm_struct
921 TYPE(cp_fm_type) :: fm1, fm2
922
923 CALL timeset(routinen, handle)
924
925 CALL dbcsr_set(density_matrix, 0.0_dp)
926
927 CALL cp_fm_struct_create(fm_struct, ncol_global=1, &
928 template_fmstruct=cvec%matrix_struct)
929 CALL cp_fm_create(fm1, fm_struct)
930 CALL cp_fm_create(fm2, fm_struct)
931 CALL cp_fm_struct_release(fm_struct)
932
933 CALL cp_fm_get_info(cvec, ncol_global=ncol)
934 DO ic = 1, ncol
935 CALL cp_fm_to_fm(cvec, fm1, ncol=1, source_start=ic, target_start=1)
936 DO jc = 1, ncol
937 CALL cp_fm_to_fm(cvec, fm2, ncol=1, source_start=jc, target_start=1)
938 CALL cp_fm_get_element(occmat, ic, jc, alpha)
939 alpha = weight(ic)*alpha
940 CALL cp_dbcsr_plus_fm_fm_t(density_matrix, fm1, fm2, ncol=1, &
941 alpha=alpha, symmetry_mode=1)
942 END DO
943 END DO
944 CALL cp_fm_release(fm1)
945 CALL cp_fm_release(fm2)
946
947 CALL timestop(handle)
948
949 END SUBROUTINE iao_calculate_dmat_full
950
951! **************************************************************************************************
952!> \brief compute the atomic charges (orthogonal basis)
953!> \param p_matrix ...
954!> \param charges ...
955! **************************************************************************************************
956 SUBROUTINE iao_charges(p_matrix, charges)
957 TYPE(dbcsr_type) :: p_matrix
958 REAL(kind=dp), DIMENSION(:) :: charges
959
960 INTEGER :: i, iblock_col, iblock_row
961 REAL(kind=dp) :: trace
962 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block
963 TYPE(dbcsr_iterator_type) :: iter
964 TYPE(mp_comm_type) :: group
965
966 charges = 0.0_dp
967
968 CALL dbcsr_iterator_start(iter, p_matrix)
969 DO WHILE (dbcsr_iterator_blocks_left(iter))
970 NULLIFY (p_block)
971 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, p_block)
972 cpassert(iblock_row == iblock_col)
973 trace = 0.0_dp
974 DO i = 1, SIZE(p_block, 1)
975 trace = trace + p_block(i, i)
976 END DO
977 charges(iblock_row) = trace
978 END DO
979 CALL dbcsr_iterator_stop(iter)
980
981 CALL dbcsr_get_info(p_matrix, group=group)
982 CALL group%sum(charges)
983
984 END SUBROUTINE iao_charges
985
986! **************************************************************************************************
987!> \brief Computes and prints the Cube Files for the intrinsic atomic orbitals
988!> \param qs_env ...
989!> \param print_section ...
990!> \param iao_coef ...
991!> \param basis_set_list ...
992! **************************************************************************************************
993 SUBROUTINE print_iao_cubes(qs_env, print_section, iao_coef, basis_set_list)
994 TYPE(qs_environment_type), POINTER :: qs_env
995 TYPE(section_vals_type), POINTER :: print_section
996 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: iao_coef
997 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
998
999 CHARACTER(LEN=default_path_length) :: filename, title
1000 INTEGER :: i, i_rep, ispin, ivec, iw, j, n_rep, &
1001 nat, natom, norb, nspins, nstart
1002 INTEGER, DIMENSION(:), POINTER :: atom_list, blk_sizes, first_bas, stride
1003 LOGICAL :: mpi_io
1004 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1005 TYPE(cell_type), POINTER :: cell
1006 TYPE(cp_logger_type), POINTER :: logger
1007 TYPE(dft_control_type), POINTER :: dft_control
1008 TYPE(particle_list_type), POINTER :: particles
1009 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1010 TYPE(pw_c1d_gs_type) :: wf_g
1011 TYPE(pw_env_type), POINTER :: pw_env
1012 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1013 TYPE(pw_r3d_rs_type) :: wf_r
1014 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1015 TYPE(qs_subsys_type), POINTER :: subsys
1016
1017 logger => cp_get_default_logger()
1018 stride => section_get_ivals(print_section, "STRIDE")
1019 CALL section_vals_val_get(print_section, "ATOM_LIST", n_rep_val=n_rep)
1020
1021 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1022 subsys=subsys, cell=cell, particle_set=particle_set, pw_env=pw_env, dft_control=dft_control)
1023 CALL qs_subsys_get(subsys, particles=particles)
1024
1025 CALL get_qs_env(qs_env=qs_env, natom=natom)
1026 ALLOCATE (blk_sizes(natom), first_bas(0:natom))
1027 CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_sizes, basis=basis_set_list)
1028 first_bas(0) = 0
1029 DO i = 1, natom
1030 first_bas(i) = first_bas(i - 1) + blk_sizes(i)
1031 END DO
1032
1033 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1034 CALL auxbas_pw_pool%create_pw(wf_r)
1035 CALL auxbas_pw_pool%create_pw(wf_g)
1036
1037 nspins = SIZE(iao_coef)
1038 nstart = min(1, n_rep)
1039
1040 DO ispin = 1, nspins
1041 DO i_rep = nstart, n_rep
1042 CALL cp_fm_get_info(iao_coef(ispin), ncol_global=norb)
1043 IF (i_rep == 0) THEN
1044 nat = natom
1045 ELSE
1046 CALL section_vals_val_get(print_section, "ATOM_LIST", i_rep_val=i_rep, i_vals=atom_list)
1047 nat = SIZE(atom_list)
1048 END IF
1049 DO i = 1, nat
1050 IF (i_rep == 0) THEN
1051 j = i
1052 ELSE
1053 j = atom_list(i)
1054 END IF
1055 cpassert(j >= 1 .AND. j <= natom)
1056 DO ivec = first_bas(j - 1) + 1, first_bas(j)
1057 WRITE (filename, '(a4,I5.5,a1,I1.1)') "IAO_", ivec, "_", ispin
1058 WRITE (title, *) "Intrinsic Atomic Orbitals ", ivec, " atom ", j, " spin ", ispin
1059 mpi_io = .true.
1060 iw = cp_print_key_unit_nr(logger, print_section, "", extension=".cube", &
1061 middle_name=trim(filename), file_position="REWIND", log_filename=.false., &
1062 mpi_io=mpi_io)
1063 CALL calculate_wavefunction(iao_coef(ispin), ivec, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1064 cell, dft_control, particle_set, pw_env)
1065 CALL cp_pw_to_cube(wf_r, iw, title, particles=particles, stride=stride, mpi_io=mpi_io)
1066 CALL cp_print_key_finished_output(iw, logger, print_section, "", mpi_io=mpi_io)
1067 END DO
1068 END DO
1069 END DO
1070 END DO
1071 DEALLOCATE (blk_sizes, first_bas)
1072 CALL auxbas_pw_pool%give_back_pw(wf_r)
1073 CALL auxbas_pw_pool%give_back_pw(wf_g)
1074
1075 END SUBROUTINE print_iao_cubes
1076
1077! **************************************************************************************************
1078!> \brief Computes and prints the Cube Files for the intrinsic bond orbitals
1079!> \param qs_env ...
1080!> \param print_section ...
1081!> \param ibo_coef ...
1082! **************************************************************************************************
1083 SUBROUTINE print_ibo_cubes(qs_env, print_section, ibo_coef)
1084 TYPE(qs_environment_type), POINTER :: qs_env
1085 TYPE(section_vals_type), POINTER :: print_section
1086 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: ibo_coef
1087
1088 CHARACTER(LEN=default_path_length) :: filename, title
1089 INTEGER :: i, i_rep, ispin, iw, j, n_rep, norb, &
1090 nspins, nstart, nstate
1091 INTEGER, DIMENSION(:), POINTER :: state_list, stride
1092 LOGICAL :: mpi_io
1093 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1094 TYPE(cell_type), POINTER :: cell
1095 TYPE(cp_logger_type), POINTER :: logger
1096 TYPE(dft_control_type), POINTER :: dft_control
1097 TYPE(particle_list_type), POINTER :: particles
1098 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1099 TYPE(pw_c1d_gs_type) :: wf_g
1100 TYPE(pw_env_type), POINTER :: pw_env
1101 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1102 TYPE(pw_r3d_rs_type) :: wf_r
1103 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1104 TYPE(qs_subsys_type), POINTER :: subsys
1105
1106 cpassert(ASSOCIATED(print_section))
1107
1108 logger => cp_get_default_logger()
1109 stride => section_get_ivals(print_section, "STRIDE")
1110 CALL section_vals_val_get(print_section, "STATE_LIST", n_rep_val=n_rep)
1111
1112 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1113 subsys=subsys, cell=cell, particle_set=particle_set, pw_env=pw_env, dft_control=dft_control)
1114 CALL qs_subsys_get(subsys, particles=particles)
1115
1116 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1117 CALL auxbas_pw_pool%create_pw(wf_r)
1118 CALL auxbas_pw_pool%create_pw(wf_g)
1119
1120 nspins = SIZE(ibo_coef)
1121 nstart = min(1, n_rep)
1122
1123 DO ispin = 1, nspins
1124 DO i_rep = nstart, n_rep
1125 CALL cp_fm_get_info(ibo_coef(ispin), ncol_global=norb)
1126 IF (i_rep == 0) THEN
1127 nstate = norb
1128 ELSE
1129 CALL section_vals_val_get(print_section, "STATE_LIST", i_rep_val=i_rep, i_vals=state_list)
1130 nstate = SIZE(state_list)
1131 END IF
1132 DO i = 1, nstate
1133 IF (i_rep == 0) THEN
1134 j = i
1135 ELSE
1136 j = state_list(i)
1137 END IF
1138 cpassert(j >= 1 .AND. j <= norb)
1139 WRITE (filename, '(a4,I5.5,a1,I1.1)') "IBO_", j, "_", ispin
1140 WRITE (title, *) "Intrinsic Atomic Orbitals ", j, " spin ", ispin
1141 mpi_io = .true.
1142 iw = cp_print_key_unit_nr(logger, print_section, "", extension=".cube", &
1143 middle_name=trim(filename), file_position="REWIND", log_filename=.false., &
1144 mpi_io=mpi_io)
1145 CALL calculate_wavefunction(ibo_coef(ispin), j, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1146 cell, dft_control, particle_set, pw_env)
1147 CALL cp_pw_to_cube(wf_r, iw, title, particles=particles, stride=stride, mpi_io=mpi_io)
1148 CALL cp_print_key_finished_output(iw, logger, print_section, "", mpi_io=mpi_io)
1149 END DO
1150 END DO
1151 END DO
1152
1153 CALL auxbas_pw_pool%give_back_pw(wf_r)
1154 CALL auxbas_pw_pool%give_back_pw(wf_g)
1155
1156 END SUBROUTINE print_ibo_cubes
1157
1158! **************************************************************************************************
1159!> \brief prints charge center and spreads for all orbitals
1160!> \param moments ...
1161!> \param nocc ...
1162!> \param print_section ...
1163!> \param iounit ...
1164! **************************************************************************************************
1165 SUBROUTINE print_center_spread(moments, nocc, print_section, iounit)
1166 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: moments
1167 INTEGER, DIMENSION(:), INTENT(IN) :: nocc
1168 TYPE(section_vals_type), OPTIONAL, POINTER :: print_section
1169 INTEGER, INTENT(IN), OPTIONAL :: iounit
1170
1171 CHARACTER(LEN=default_path_length) :: filename
1172 INTEGER :: is, ispin, iw, nspin
1173 TYPE(cp_logger_type), POINTER :: logger
1174
1175 logger => cp_get_default_logger()
1176 nspin = SIZE(moments, 3)
1177
1178 IF (PRESENT(print_section)) THEN
1179 WRITE (filename, '(A18,I1.1)') "IBO_CENTERS_SPREAD"
1180 iw = cp_print_key_unit_nr(logger, print_section, "", extension=".csp", &
1181 middle_name=trim(filename), file_position="REWIND", log_filename=.false.)
1182 ELSE IF (PRESENT(iounit)) THEN
1183 iw = iounit
1184 ELSE
1185 iw = -1
1186 END IF
1187 IF (iw > 0) THEN
1188 DO ispin = 1, nspin
1189 WRITE (iw, "(T2,A,i1)") "Intrinsic Bond Orbitals: Centers and Spread for Spin ", ispin
1190 WRITE (iw, "(A7,T30,A6,T68,A7)") "State", "Center", "Spreads"
1191 DO is = 1, nocc(ispin)
1192 WRITE (iw, "(i7,3F15.8,8X,2F10.5)") is, moments(:, is, ispin)
1193 END DO
1194 END DO
1195 END IF
1196 IF (PRESENT(print_section)) THEN
1197 CALL cp_print_key_finished_output(iw, logger, print_section, "")
1198 END IF
1199
1200 END SUBROUTINE print_center_spread
1201
1202! **************************************************************************************************
1203!> \brief ...
1204!> \param qs_env ...
1205!> \param c_loc_orb ...
1206!> \param moments ...
1207! **************************************************************************************************
1208 SUBROUTINE center_spread_loc(qs_env, c_loc_orb, moments)
1209 TYPE(qs_environment_type), POINTER :: qs_env
1210 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: c_loc_orb
1211 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: moments
1212
1213 CHARACTER(len=*), PARAMETER :: routinen = 'center_spread_loc'
1214 INTEGER, DIMENSION(6), PARAMETER :: list = [1, 2, 3, 4, 7, 9]
1215
1216 INTEGER :: handle, i, iop, iorb, ispin, nao, norb, &
1217 nspin
1218 REAL(kind=dp) :: xmii
1219 REAL(kind=dp), DIMENSION(3) :: rpoint
1220 TYPE(cell_type), POINTER :: cell
1221 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1222 TYPE(cp_fm_type) :: ccmat, ocvec
1223 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat
1224 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1225 TYPE(dbcsr_type), POINTER :: omat, smat
1226
1227 CALL timeset(routinen, handle)
1228
1229 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
1230 smat => matrix_s(1, 1)%matrix
1231 rpoint = 0.0_dp
1232 nspin = SIZE(c_loc_orb)
1233 moments = 0.0_dp
1234
1235 ALLOCATE (dipmat(9))
1236 DO i = 1, 9
1237 ALLOCATE (dipmat(i)%matrix)
1238 CALL dbcsr_copy(dipmat(i)%matrix, smat)
1239 CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
1240 END DO
1241
1242 CALL build_local_moment_matrix(qs_env, dipmat, 2, rpoint)
1243
1244 DO ispin = 1, nspin
1245 CALL cp_fm_get_info(c_loc_orb(ispin), nrow_global=nao, ncol_global=norb)
1246 CALL cp_fm_create(ocvec, c_loc_orb(ispin)%matrix_struct)
1247 NULLIFY (fm_struct)
1248 CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=norb, &
1249 template_fmstruct=c_loc_orb(ispin)%matrix_struct)
1250 CALL cp_fm_create(ccmat, fm_struct)
1251 CALL cp_fm_struct_release(fm_struct)
1252 DO i = 1, 6
1253 iop = list(i)
1254 omat => dipmat(iop)%matrix
1255 CALL cp_dbcsr_sm_fm_multiply(omat, c_loc_orb(ispin), ocvec, ncol=norb)
1256 CALL parallel_gemm('T', 'N', norb, norb, nao, 1.0_dp, c_loc_orb(ispin), ocvec, &
1257 0.0_dp, ccmat)
1258 DO iorb = 1, norb
1259 CALL cp_fm_get_element(ccmat, iorb, iorb, xmii)
1260 IF (iop <= 3) THEN
1261 moments(iop, iorb, ispin) = moments(iop, iorb, ispin) + xmii
1262 moments(4, iorb, ispin) = moments(4, iorb, ispin) - xmii**2
1263 ELSE
1264 moments(4, iorb, ispin) = moments(4, iorb, ispin) + xmii
1265 END IF
1266 END DO
1267 END DO
1268 CALL cp_fm_release(ocvec)
1269 CALL cp_fm_release(ccmat)
1270 END DO
1271
1272 DO i = 1, 9
1273 CALL dbcsr_release(dipmat(i)%matrix)
1274 DEALLOCATE (dipmat(i)%matrix)
1275 END DO
1276 DEALLOCATE (dipmat)
1277
1278 CALL get_qs_env(qs_env=qs_env, cell=cell)
1279 DO ispin = 1, nspin
1280 DO iorb = 1, norb
1281 moments(1:3, iorb, ispin) = pbc(moments(1:3, iorb, ispin), cell)
1282 END DO
1283 END DO
1284
1285 CALL timestop(handle)
1286
1287 END SUBROUTINE center_spread_loc
1288
1289! **************************************************************************************************
1290!> \brief ...
1291!> \param qs_env ...
1292!> \param c_loc_orb ...
1293!> \param moments ...
1294! **************************************************************************************************
1295 SUBROUTINE center_spread_berry(qs_env, c_loc_orb, moments)
1296 TYPE(qs_environment_type), POINTER :: qs_env
1297 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: c_loc_orb
1298 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: moments
1299
1300 CHARACTER(len=*), PARAMETER :: routinen = 'center_spread_berry'
1301
1302 COMPLEX(KIND=dp) :: z
1303 INTEGER :: dim_op, handle, i, idir, ispin, istate, &
1304 j, jdir, nao, norb, nspin
1305 REAL(dp), DIMENSION(3) :: c, cpbc
1306 REAL(dp), DIMENSION(6) :: weights
1307 REAL(kind=dp) :: imagpart, realpart, spread_i, spread_ii
1308 TYPE(cell_type), POINTER :: cell
1309 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1310 TYPE(cp_fm_type) :: opvec
1311 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: zij
1312 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1313 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
1314
1315 CALL timeset(routinen, handle)
1316
1317 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, cell=cell)
1318
1319 IF (cell%orthorhombic) THEN
1320 dim_op = 3
1321 ELSE
1322 dim_op = 6
1323 END IF
1324 ALLOCATE (op_sm_set(2, dim_op))
1325 DO i = 1, dim_op
1326 DO j = 1, 2
1327 NULLIFY (op_sm_set(j, i)%matrix)
1328 ALLOCATE (op_sm_set(j, i)%matrix)
1329 CALL dbcsr_copy(op_sm_set(j, i)%matrix, matrix_s(1)%matrix)
1330 END DO
1331 END DO
1332 CALL initialize_weights(cell, weights)
1333
1334 CALL compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
1335
1336 nspin = SIZE(c_loc_orb, 1)
1337 DO ispin = 1, nspin
1338 CALL cp_fm_get_info(c_loc_orb(ispin), nrow_global=nao, ncol_global=norb)
1339 CALL cp_fm_create(opvec, c_loc_orb(ispin)%matrix_struct)
1340 NULLIFY (fm_struct)
1341 CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=norb, &
1342 template_fmstruct=c_loc_orb(ispin)%matrix_struct)
1343 ALLOCATE (zij(2, dim_op))
1344
1345 ! Compute zij here
1346 DO i = 1, dim_op
1347 DO j = 1, 2
1348 CALL cp_fm_create(zij(j, i), fm_struct)
1349 CALL cp_fm_set_all(zij(j, i), 0.0_dp)
1350 CALL cp_dbcsr_sm_fm_multiply(op_sm_set(j, i)%matrix, c_loc_orb(ispin), opvec, ncol=norb)
1351 CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, c_loc_orb(ispin), opvec, 0.0_dp, zij(j, i))
1352 END DO
1353 END DO
1354
1355 CALL cp_fm_struct_release(fm_struct)
1356 CALL cp_fm_release(opvec)
1357
1358 DO istate = 1, norb
1359 c = 0.0_dp
1360 spread_i = 0.0_dp
1361 spread_ii = 0.0_dp
1362 DO jdir = 1, dim_op
1363 CALL cp_fm_get_element(zij(1, jdir), istate, istate, realpart)
1364 CALL cp_fm_get_element(zij(2, jdir), istate, istate, imagpart)
1365 z = cmplx(realpart, imagpart, dp)
1366 spread_i = spread_i - weights(jdir)* &
1367 log(realpart*realpart + imagpart*imagpart)/twopi/twopi
1368 spread_ii = spread_ii + weights(jdir)* &
1369 (1.0_dp - (realpart*realpart + imagpart*imagpart))/twopi/twopi
1370 IF (jdir < 4) THEN
1371 DO idir = 1, 3
1372 c(idir) = c(idir) + &
1373 (cell%hmat(idir, jdir)/twopi)*aimag(log(z))
1374 END DO
1375 END IF
1376 END DO
1377 cpbc = pbc(c, cell)
1378 moments(1:3, istate, ispin) = cpbc(1:3)
1379 moments(4, istate, ispin) = spread_i
1380 moments(5, istate, ispin) = spread_ii
1381 END DO
1382
1383 CALL cp_fm_release(zij)
1384
1385 END DO
1386
1387 DO i = 1, dim_op
1388 DO j = 1, 2
1389 CALL dbcsr_release(op_sm_set(j, i)%matrix)
1390 DEALLOCATE (op_sm_set(j, i)%matrix)
1391 END DO
1392 END DO
1393 DEALLOCATE (op_sm_set)
1394
1395 CALL timestop(handle)
1396
1397 END SUBROUTINE center_spread_berry
1398
1399! **************************************************************************************************
1400!> \brief ...
1401!> \param qs_env ...
1402!> \param oce_basis_set_list ...
1403!> \param smat_kind ...
1404!> \param sxo ...
1405!> \param iao_coef ...
1406!> \param oce_atom ...
1407!> \param unit_nr ...
1408! **************************************************************************************************
1409 SUBROUTINE iao_oce_expansion(qs_env, oce_basis_set_list, smat_kind, sxo, iao_coef, oce_atom, unit_nr)
1410 TYPE(qs_environment_type), POINTER :: qs_env
1411 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: oce_basis_set_list
1412 TYPE(cp_2d_r_p_type), DIMENSION(:) :: smat_kind
1413 TYPE(dbcsr_p_type), DIMENSION(:) :: sxo
1414 TYPE(cp_fm_type), DIMENSION(:) :: iao_coef
1415 TYPE(cp_2d_r_p_type), DIMENSION(:, :) :: oce_atom
1416 INTEGER, INTENT(IN) :: unit_nr
1417
1418 CHARACTER(len=*), PARAMETER :: routinen = 'iao_oce_expansion'
1419
1420 INTEGER :: handle, i, i1, i2, iao, iatom, ikind, &
1421 iset, ishell, ispin, l, m, maxl, n, &
1422 natom, nkind, noce, ns, nset, nsgf, &
1423 nspin
1424 INTEGER, DIMENSION(:), POINTER :: iao_blk_sizes, nshell, oce_blk_sizes, &
1425 orb_blk_sizes
1426 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, last_sgf, lval
1427 LOGICAL :: found
1428 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ev, vector
1429 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat, filter, oce_comp, prol, vec
1430 REAL(kind=dp), DIMENSION(:, :), POINTER :: ablock
1431 TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: sinv_kind
1432 TYPE(dbcsr_distribution_type) :: dbcsr_dist
1433 TYPE(dbcsr_type) :: iao_vec, sx_vec
1434 TYPE(gto_basis_set_type), POINTER :: oce_basis
1435 TYPE(mp_comm_type) :: group
1436 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1437 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1438 TYPE(qs_ks_env_type), POINTER :: ks_env
1439
1440 CALL timeset(routinen, handle)
1441
1442 ! basis sets: block sizes
1443 CALL dbcsr_get_info(sxo(1)%matrix, row_blk_size=oce_blk_sizes, col_blk_size=orb_blk_sizes, &
1444 distribution=dbcsr_dist)
1445 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, natom=natom)
1446 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, qs_kind_set=qs_kind_set)
1447 ALLOCATE (iao_blk_sizes(natom))
1448 DO iatom = 1, natom
1449 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1450 CALL get_qs_kind(qs_kind_set(ikind), nsgf=ns, basis_type="MIN")
1451 iao_blk_sizes(iatom) = ns
1452 END DO
1453
1454 CALL dbcsr_create(iao_vec, name="IAO_VEC", dist=dbcsr_dist, &
1455 matrix_type=dbcsr_type_no_symmetry, row_blk_size=orb_blk_sizes, &
1456 col_blk_size=iao_blk_sizes)
1457 CALL dbcsr_create(sx_vec, name="SX_VEC", dist=dbcsr_dist, &
1458 matrix_type=dbcsr_type_no_symmetry, row_blk_size=oce_blk_sizes, &
1459 col_blk_size=iao_blk_sizes)
1460 CALL dbcsr_reserve_diag_blocks(matrix=sx_vec)
1461
1462 nkind = SIZE(smat_kind)
1463 ALLOCATE (sinv_kind(nkind))
1464 DO ikind = 1, nkind
1465 noce = SIZE(smat_kind(ikind)%array, 1)
1466 ALLOCATE (sinv_kind(ikind)%array(noce, noce))
1467 sinv_kind(ikind)%array = smat_kind(ikind)%array
1468 CALL invmat_symm(sinv_kind(ikind)%array)
1469 END DO
1470 CALL dbcsr_get_info(iao_vec, group=group)
1471
1472 nspin = SIZE(iao_coef, 1)
1473 ALLOCATE (oce_comp(natom, nspin))
1474 oce_comp = 0.0_dp
1475 DO ispin = 1, nspin
1476 CALL copy_fm_to_dbcsr(iao_coef(ispin), iao_vec)
1477 CALL dbcsr_multiply("N", "N", 1.0_dp, sxo(1)%matrix, iao_vec, 0.0_dp, sx_vec, &
1478 retain_sparsity=.true.)
1479 DO iatom = 1, natom
1480 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1481 CALL dbcsr_get_block_p(matrix=sx_vec, &
1482 row=iatom, col=iatom, block=ablock, found=found)
1483 IF (found) THEN
1484 n = SIZE(ablock, 1)
1485 m = SIZE(ablock, 2)
1486 ablock = matmul(sinv_kind(ikind)%array, ablock)
1487 ALLOCATE (amat(n, m), bmat(m, m), ev(m), vec(m, m))
1488 amat(1:n, 1:m) = matmul(smat_kind(ikind)%array(1:n, 1:n), ablock(1:n, 1:m))
1489 bmat(1:m, 1:m) = matmul(transpose(ablock(1:n, 1:m)), amat(1:n, 1:m))
1490 CALL jacobi(bmat, ev, vec)
1491 oce_comp(iatom, ispin) = sum(ev(1:m))/real(m, kind=dp)
1492 DO i = 1, m
1493 ev(i) = 1._dp/sqrt(ev(i))
1494 bmat(1:m, i) = vec(1:m, i)*ev(i)
1495 END DO
1496 bmat(1:m, 1:m) = matmul(bmat(1:m, 1:m), transpose(vec(1:m, 1:m)))
1497 ablock(1:n, 1:m) = matmul(ablock(1:n, 1:m), bmat(1:m, 1:m))
1498 DEALLOCATE (amat, bmat, ev, vec)
1499 END IF
1500 END DO
1501 CALL dbcsr_replicate_all(sx_vec)
1502 DO iatom = 1, natom
1503 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1504 CALL dbcsr_get_block_p(matrix=sx_vec, &
1505 row=iatom, col=iatom, block=ablock, found=found)
1506 cpassert(found)
1507 n = SIZE(ablock, 1)
1508 m = SIZE(ablock, 2)
1509 ALLOCATE (oce_atom(iatom, ispin)%array(n, m))
1510 oce_atom(iatom, ispin)%array = ablock
1511 END DO
1512 END DO
1513
1514 CALL group%sum(oce_comp)
1515 IF (unit_nr > 0) THEN
1516 DO iatom = 1, natom
1517 WRITE (unit_nr, "(T4,A,I6,T30,A,2F12.4)") "Atom:", iatom, " Completeness: ", oce_comp(iatom, 1:nspin)
1518 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1519 oce_basis => oce_basis_set_list(ikind)%gto_basis_set
1520 CALL get_gto_basis_set(oce_basis, nset=nset, nshell=nshell, nsgf=nsgf, maxl=maxl, &
1521 l=lval, first_sgf=first_sgf, last_sgf=last_sgf)
1522 ALLOCATE (filter(nsgf, 0:maxl), vector(nsgf), prol(0:maxl, nspin))
1523 filter = 0.0_dp
1524 DO iset = 1, nset
1525 DO ishell = 1, nshell(iset)
1526 l = lval(ishell, iset)
1527 i1 = first_sgf(ishell, iset)
1528 i2 = last_sgf(ishell, iset)
1529 filter(i1:i2, l) = 1.0_dp
1530 END DO
1531 END DO
1532 !
1533 n = SIZE(oce_atom(iatom, 1)%array, 1)
1534 m = SIZE(oce_atom(iatom, 1)%array, 2)
1535 cpassert(n == nsgf)
1536 DO iao = 1, m
1537 prol = 0.0_dp
1538 DO ispin = 1, nspin
1539 DO l = 0, maxl
1540 vector(1:n) = oce_atom(iatom, ispin)%array(1:n, iao)*filter(1:n, l)
1541 prol(l, ispin) = sum(matmul(smat_kind(ikind)%array(1:n, 1:n), vector(1:n))*vector(1:n))
1542 END DO
1543 END DO
1544 WRITE (unit_nr, "(T4,I3,T15,A,T39,(6F7.4))") iao, " l-contributions:", (sum(prol(l, :)), l=0, maxl)
1545 END DO
1546 DEALLOCATE (filter, vector, prol)
1547 END DO
1548 WRITE (unit_nr, *)
1549 END IF
1550
1551 DEALLOCATE (oce_comp)
1552 DO ikind = 1, nkind
1553 DEALLOCATE (sinv_kind(ikind)%array)
1554 END DO
1555 DEALLOCATE (sinv_kind)
1556 DEALLOCATE (iao_blk_sizes)
1557 CALL dbcsr_release(iao_vec)
1558 CALL dbcsr_release(sx_vec)
1559
1560 CALL timestop(handle)
1561
1562 END SUBROUTINE iao_oce_expansion
1563
1564! **************************************************************************************************
1565!> \brief ...
1566!> \param smat_kind ...
1567!> \param basis_set_list ...
1568! **************************************************************************************************
1569 SUBROUTINE oce_overlap_matrix(smat_kind, basis_set_list)
1570 TYPE(cp_2d_r_p_type), DIMENSION(:) :: smat_kind
1571 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1572
1573 CHARACTER(len=*), PARAMETER :: routinen = 'oce_overlap_matrix'
1574
1575 INTEGER :: handle, ikind, iset, jset, ldsab, m1, &
1576 m2, n1, n2, ncoa, ncob, nkind, nseta, &
1577 sgfa, sgfb
1578 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
1579 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
1580 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: oint, owork
1581 REAL(kind=dp), DIMENSION(3) :: rab
1582 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, scon_a, smat, zeta
1583 TYPE(gto_basis_set_type), POINTER :: basis_set
1584
1585 CALL timeset(routinen, handle)
1586
1587 rab(1:3) = 0.0_dp
1588
1589 nkind = SIZE(smat_kind)
1590 ldsab = 0
1591 DO ikind = 1, nkind
1592 basis_set => basis_set_list(ikind)%gto_basis_set
1593 cpassert(ASSOCIATED(basis_set))
1594 CALL get_gto_basis_set(gto_basis_set=basis_set, maxco=m1, nsgf=m2)
1595 ldsab = max(m1, m2, ldsab)
1596 END DO
1597
1598 ALLOCATE (oint(ldsab, ldsab), owork(ldsab, ldsab))
1599
1600 DO ikind = 1, nkind
1601 basis_set => basis_set_list(ikind)%gto_basis_set
1602 cpassert(ASSOCIATED(basis_set))
1603 smat => smat_kind(ikind)%array
1604 smat = 0.0_dp
1605 ! basis ikind
1606 first_sgfa => basis_set%first_sgf
1607 la_max => basis_set%lmax
1608 la_min => basis_set%lmin
1609 npgfa => basis_set%npgf
1610 nseta = basis_set%nset
1611 nsgfa => basis_set%nsgf_set
1612 rpgfa => basis_set%pgf_radius
1613 scon_a => basis_set%scon
1614 zeta => basis_set%zet
1615
1616 DO iset = 1, nseta
1617
1618 ncoa = npgfa(iset)*ncoset(la_max(iset))
1619 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1620 sgfa = first_sgfa(1, iset)
1621
1622 DO jset = 1, nseta
1623
1624 ncob = npgfa(jset)*ncoset(la_max(jset))
1625 n2 = npgfa(jset)*(ncoset(la_max(jset)) - ncoset(la_min(jset) - 1))
1626 sgfb = first_sgfa(1, jset)
1627
1628 ! calculate integrals and derivatives
1629 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1630 la_max(jset), la_min(jset), npgfa(jset), rpgfa(:, jset), zeta(:, jset), &
1631 rab, sab=oint(:, :))
1632 ! Contraction
1633 CALL contraction(oint(:, :), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1634 cb=scon_a(:, sgfb:), nb=n2, mb=nsgfa(jset), fscale=1.0_dp, trans=.false.)
1635 CALL block_add("IN", owork, nsgfa(iset), nsgfa(jset), smat, &
1636 sgfa, sgfb, trans=.false.)
1637
1638 END DO
1639 END DO
1640
1641 END DO
1642 DEALLOCATE (oint, owork)
1643
1644 CALL timestop(handle)
1645
1646 END SUBROUTINE oce_overlap_matrix
1647
1648! **************************************************************************************************
1649!> \brief 2 by 2 rotation optimization for the Pipek-Mezey operator
1650!> \param zij_fm_set ...
1651!> \param vectors ...
1652!> \param order ...
1653!> \param accuracy ...
1654!> \param unit_nr ...
1655!> \par History
1656!> 05-2005 created
1657!> 10-2023 adapted from jacobi_rotation_pipek [JGH]
1658!> \author MI
1659! **************************************************************************************************
1660 SUBROUTINE pm_localization(zij_fm_set, vectors, order, accuracy, unit_nr)
1661
1662 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: zij_fm_set
1663 TYPE(cp_fm_type), INTENT(IN) :: vectors
1664 INTEGER, INTENT(IN) :: order
1665 REAL(kind=dp), INTENT(IN) :: accuracy
1666 INTEGER, INTENT(IN) :: unit_nr
1667
1668 INTEGER, PARAMETER :: max_sweeps = 250
1669
1670 INTEGER :: iatom, istate, jstate, natom, nstate, &
1671 sweeps
1672 REAL(kind=dp) :: aij, bij, ct, ftarget, mii, mij, mjj, &
1673 st, t1, t2, theta, tolerance, tt
1674 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fdiag
1675 TYPE(cp_fm_type) :: rmat
1676
1677 CALL cp_fm_create(rmat, zij_fm_set(1, 1)%matrix_struct)
1678 CALL cp_fm_set_all(rmat, 0.0_dp, 1.0_dp)
1679
1680 CALL cp_fm_get_info(rmat, nrow_global=nstate)
1681 ALLOCATE (fdiag(nstate))
1682 tolerance = 1.0e10_dp
1683 sweeps = 0
1684 natom = SIZE(zij_fm_set, 1)
1685 IF (unit_nr > 0) THEN
1686 WRITE (unit_nr, '(A,T30,A,T45,A,T63,A,T77,A)') &
1687 " Jacobi Localization ", "Sweep", "Functional", "Gradient", "Time"
1688 END IF
1689 ! do jacobi sweeps until converged
1690 DO sweeps = 1, max_sweeps
1691 t1 = m_walltime()
1692 tolerance = 0.0_dp
1693 DO istate = 1, nstate
1694 DO jstate = istate + 1, nstate
1695 aij = 0.0_dp
1696 bij = 0.0_dp
1697 DO iatom = 1, natom
1698 CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, istate, mii)
1699 CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, jstate, mij)
1700 CALL cp_fm_get_element(zij_fm_set(iatom, 1), jstate, jstate, mjj)
1701 IF (order == 2) THEN
1702 aij = aij + 4._dp*mij**2 - (mii - mjj)**2
1703 bij = bij + 4._dp*mij*(mii - mjj)
1704 ELSE IF (order == 4) THEN
1705 aij = aij - mii**4 - mjj**4 + 6._dp*(mii**2 + mjj**2)*mij**2 + &
1706 mii**3*mjj + mii*mjj**3
1707 bij = bij + 4._dp*mij*(mii**3 - mjj**3)
1708 ELSE
1709 cpabort("PM order")
1710 END IF
1711 END DO
1712 IF ((aij**2 + bij**2) < 1.e-10_dp) cycle
1713 tolerance = tolerance + bij**2
1714 theta = 0.25_dp*atan2(bij, -aij)
1715 ct = cos(theta)
1716 st = sin(theta)
1717
1718 CALL cp_fm_rot_cols(rmat, istate, jstate, ct, st)
1719
1720 DO iatom = 1, natom
1721 CALL cp_fm_rot_cols(zij_fm_set(iatom, 1), istate, jstate, ct, st)
1722 CALL cp_fm_rot_rows(zij_fm_set(iatom, 1), istate, jstate, ct, st)
1723 END DO
1724 END DO
1725 END DO
1726 ftarget = 0.0_dp
1727 DO iatom = 1, natom
1728 CALL cp_fm_get_diag(zij_fm_set(iatom, 1), fdiag)
1729 ftarget = ftarget + sum(fdiag**order)
1730 END DO
1731 ftarget = ftarget**(1._dp/order)
1732 tolerance = sqrt(tolerance)
1733 t2 = m_walltime()
1734 tt = t2 - t1
1735 IF (unit_nr > 0) THEN
1736 WRITE (unit_nr, '(T31,I4,T39,F16.8,T55,G16.8,T71,F10.2)') sweeps, ftarget, tolerance, tt
1737 END IF
1738 IF (tolerance < accuracy) EXIT
1739 END DO
1740 DEALLOCATE (fdiag)
1741
1742 CALL rotate_orbitals(rmat, vectors)
1743 CALL cp_fm_release(rmat)
1744
1745 END SUBROUTINE pm_localization
1746! **************************************************************************************************
1747
1748END MODULE iao_analysis
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab, rr_work)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:273
simple routine to print charges for all atomic charge methods (currently mulliken,...
subroutine, public print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges, atomic_charges)
generates a unified output format for atomic charges
Define the atomic kind types and their sub types.
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.
Automatic generation of auxiliary basis sets of different kind.
Definition auto_basis.F:15
subroutine, public create_oce_basis(oce_basis, orb_basis, lmax_oce, nbas_oce)
...
Definition auto_basis.F:432
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
subroutine, public deallocate_gto_basis_set(gto_basis_set)
...
subroutine, public init_orb_basis_set(gto_basis_set)
Initialise a Gaussian-type orbital (GTO) basis set data set.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public knizia2013
Handles all functions related to the CELL.
Definition cell_types.F:15
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_replicate_all(matrix)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_trace(matrix, trace)
Computes the trace of the given matrix, also known as the sum of its diagonal elements.
subroutine, public dbcsr_reserve_diag_blocks(matrix)
Reserves all diagonal blocks.
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_rot_rows(matrix, irow, jrow, cs, sn)
Applies a planar rotation defined by cs and sn to the i'th and j'th rows.
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
subroutine, public cp_fm_rot_cols(matrix, icol, jcol, cs, sn)
Applies a planar rotation defined by cs and sn to the i'th and j'th columnns.
subroutine, public cp_fm_geadd(alpha, trans, matrix_a, beta, matrix_b)
interface to BLACS geadd: matrix_b = beta*matrix_b + alpha*opt(matrix_a) where opt(matrix_a) can be e...
subroutine, public cp_fm_invert(matrix_a, matrix_inverse, det_a, eps_svd, eigval)
Inverts a cp_fm_type matrix, optionally returning the determinant of the input matrix.
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
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)
...
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,...
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io)
...
Calculate intrinsic atomic orbitals and analyze wavefunctions.
subroutine, public iao_wfn_analysis(qs_env, iao_env, unit_nr, c_iao_coef, mos, bond_centers)
...
subroutine, public print_center_spread(moments, nocc, print_section, iounit)
prints charge center and spreads for all orbitals
subroutine, public iao_charges(p_matrix, charges)
compute the atomic charges (orthogonal basis)
Calculate ntrinsic atomic orbitals and analyze wavefunctions.
Definition iao_types.F:14
subroutine, public organise_printout_for_rose(qs_env, orb_basis_set_list)
Prints fragmented overlap matrices for interface with ROSE.
Definition iao_types.F:210
subroutine, public print_fragment_association(qs_env, moments, unit_nr)
Associates IBO centers to fragments.
Definition iao_types.F:332
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_iaoloc_enone
integer, parameter, public do_iaoloc_pm2
objects that represent the structure of input sections and the data contained in an input section
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
integer, parameter, public default_path_length
Definition kinds.F:58
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
Definition of mathematical constants and functions.
real(kind=dp), parameter, public twopi
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
subroutine, public jacobi(a, d, v)
Jacobi matrix diagonalization. The eigenvalues are returned in vector d and the eigenvectors are retu...
Definition mathlib.F:1605
subroutine, public invmat_symm(a, potrf, uplo)
returns inverse of real symmetric, positive definite matrix
Definition mathlib.F:588
Interface to the message passing library MPI.
generate or use from input minimal basis set
subroutine, public create_minbas_set(qs_env, unit_nr, basis_type, primitive)
...
Functions handling the MOLDEN format. Split from mode_selective.
subroutine, public write_mos_molden(mos, qs_kind_set, particle_set, print_section, cell, unoccupied_orbs, unoccupied_evals, qs_env, calc_energies)
Write out the MOs in molden format for visualisation.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
basic linear algebra operations for full matrixes
represent a simple array based list of the given type
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Define the data structure for the particle information.
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type)
maps a given wavefunction on the grid
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Calculate the interaction radii for the operator matrix calculation.
subroutine, public init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short)
...
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Some utilities for the construction of the localization environment.
subroutine, public compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
Computes the reciprocal-space operators used by Berry localization. The operators contain the contrac...
Localization methods such as 2x2 Jacobi rotations Steepest Decents Conjugate Gradient.
subroutine, public rotate_orbitals(rmat, vectors)
...
subroutine, public initialize_weights(cell, weights)
...
subroutine, public scdm_qrfact(vectors)
...
collects routines that perform operations directly related to MOs
subroutine, public make_basis_lowdin(vmatrix, ncol, matrix_s)
return a set of S orthonormal vectors (C^T S C == 1) where a Loedwin transformation is applied to kee...
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public set_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, uniform_occupation, kts, mu, flexible_electron_count)
Set the components of a MO set data structure.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name, counter, complex_coeff)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count, cmo_coeff)
Get the components of a MO set data structure.
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>.
Definition qs_moments.F:14
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type, all_images, minimum_image, neighbor_image, first_component)
...
Definition qs_moments.F:166
Define the neighbor list data types and the corresponding functionality.
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Generate the atomic neighbor lists.
subroutine, public setup_neighbor_list(ab_list, basis_set_a, basis_set_b, qs_env, mic, symmetric, molecular, operator_type)
Build a neighborlist.
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
subroutine, public build_overlap_matrix_simple(ks_env, matrix_s, basis_set_list_a, basis_set_list_b, sab_nl, lcart)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:575
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
represent a pointer to a 2d array
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
contained for different pw related things
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...