(git:4733e3f)
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
73 USE kinds, ONLY: default_path_length,&
75 dp
76 USE machine, ONLY: m_walltime
77 USE mathconstants, ONLY: twopi
78 USE mathlib, ONLY: invmat_symm,&
79 jacobi
83 USE orbital_pointers, ONLY: ncoset
88 USE pw_env_types, ONLY: pw_env_get,&
91 USE pw_types, ONLY: pw_c1d_gs_type,&
97 USE qs_kind_types, ONLY: get_qs_kind,&
99 USE qs_ks_types, ONLY: get_ks_env,&
106 USE qs_mo_types, ONLY: allocate_mo_set,&
108 get_mo_set,&
117 USE qs_subsys_types, ONLY: qs_subsys_get,&
119#include "./base/base_uses.f90"
120
121 IMPLICIT NONE
122 PRIVATE
123
124 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'iao_analysis'
125
127
129 MODULE PROCEDURE iao_calculate_dmat_diag, & ! diagonal occupation matrix
130 iao_calculate_dmat_full ! full occupation matrix
131 END INTERFACE
132
133! **************************************************************************************************
134
135CONTAINS
136
137! **************************************************************************************************
138!> \brief ...
139!> \param qs_env ...
140!> \param iao_env ...
141!> \param unit_nr ...
142!> \param c_iao_coef ...
143!> \param mos ...
144!> \param bond_centers ...
145! **************************************************************************************************
146 SUBROUTINE iao_wfn_analysis(qs_env, iao_env, unit_nr, c_iao_coef, mos, bond_centers)
147 TYPE(qs_environment_type), POINTER :: qs_env
148 TYPE(iao_env_type), INTENT(IN) :: iao_env
149 INTEGER, INTENT(IN) :: unit_nr
150 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
151 OPTIONAL :: c_iao_coef
152 TYPE(mo_set_type), DIMENSION(:), OPTIONAL, POINTER :: mos
153 REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL :: bond_centers
154
155 CHARACTER(len=*), PARAMETER :: routinen = 'iao_wfn_analysis'
156
157 CHARACTER(LEN=2) :: element_symbol
158 CHARACTER(LEN=default_string_length) :: bname
159 INTEGER :: handle, i, iatom, ikind, isgf, ispin, &
160 nao, natom, nimages, nkind, no, norb, &
161 nref, ns, nsgf, nspin, nvec, nx, order
162 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf, nsgfat
163 INTEGER, DIMENSION(2) :: nocc
164 LOGICAL :: cubes_iao, cubes_ibo, molden_iao, &
165 molden_ibo, uniform_occupation
166 REAL(kind=dp) :: fin, fout, t1, t2, trace, zval
167 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fdiag
168 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mcharge
169 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: moments
170 REAL(kind=dp), DIMENSION(:), POINTER :: occupation_numbers
171 TYPE(cell_type), POINTER :: cell
172 TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: smat_kind
173 TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:, :) :: oce_atom
174 TYPE(cp_fm_struct_type), POINTER :: fm_struct
175 TYPE(cp_fm_type) :: p_orb_ref, p_ref_orb, smo
176 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: c_loc_orb, ciao, cloc, cvec, iao_coef
177 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: zij_atom
178 TYPE(cp_fm_type), POINTER :: mo_coeff
179 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: sref, sro, sxo
180 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
181 TYPE(dbcsr_type) :: dmat
182 TYPE(dbcsr_type), POINTER :: smat
183 TYPE(dft_control_type), POINTER :: dft_control
184 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: oce_basis_set_list, orb_basis_set_list, &
185 ref_basis_set_list
186 TYPE(gto_basis_set_type), POINTER :: ocebasis, orbbasis, refbasis
187 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mo_iao, mo_loc
188 TYPE(mo_set_type), DIMENSION(:), POINTER :: my_mos
189 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
190 POINTER :: sro_list, srr_list, sxo_list
191 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
192 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
193 TYPE(qs_kind_type), POINTER :: qs_kind
194 TYPE(qs_ks_env_type), POINTER :: ks_env
195 TYPE(section_vals_type), POINTER :: iao_cubes_section, iao_molden_section, &
196 ibo_cc_section, ibo_cubes_section, &
197 ibo_molden_section
198
199 ! only do IAO analysis if explicitly requested
200 cpassert(iao_env%do_iao)
201
202 ! k-points?
203 CALL get_qs_env(qs_env, dft_control=dft_control)
204 nspin = dft_control%nspins
205 nimages = dft_control%nimages
206 IF (nimages > 1) THEN
207 IF (unit_nr > 0) THEN
208 WRITE (unit=unit_nr, fmt="(T2,A)") &
209 "K-Points: Intrinsic Atomic Orbitals Analysis not available."
210 END IF
211 END IF
212 IF (nimages > 1) RETURN
213
214 CALL timeset(routinen, handle)
215
216 IF (unit_nr > 0) THEN
217 WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
218 WRITE (unit=unit_nr, fmt="(T24,A)") "INTRINSIC ATOMIC ORBITALS ANALYSIS"
219 WRITE (unit=unit_nr, fmt="(T13,A)") "G. Knizia, J. Chem. Theory Comput. 9, 4834-4843 (2013)"
220 WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!'
221 END IF
222 CALL cite_reference(knizia2013)
223
224 ! input sections
225 iao_molden_section => iao_env%iao_molden_section
226 iao_cubes_section => iao_env%iao_cubes_section
227 ibo_molden_section => iao_env%ibo_molden_section
228 ibo_cubes_section => iao_env%ibo_cubes_section
229 ibo_cc_section => iao_env%ibo_cc_section
230
231 !
232 molden_iao = .false.
233 IF (ASSOCIATED(iao_molden_section)) THEN
234 CALL section_vals_get(iao_molden_section, explicit=molden_iao)
235 END IF
236 cubes_iao = .false.
237 IF (ASSOCIATED(iao_cubes_section)) THEN
238 CALL section_vals_get(iao_cubes_section, explicit=cubes_iao)
239 END IF
240 molden_ibo = .false.
241 IF (ASSOCIATED(ibo_molden_section)) THEN
242 CALL section_vals_get(ibo_molden_section, explicit=molden_ibo)
243 END IF
244 cubes_ibo = .false.
245 IF (ASSOCIATED(ibo_cubes_section)) THEN
246 CALL section_vals_get(ibo_cubes_section, explicit=cubes_ibo)
247 END IF
248
249 ! check for or generate reference basis
250 CALL create_minbas_set(qs_env, unit_nr)
251
252 ! overlap matrices
253 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s)
254 smat => matrix_s(1, 1)%matrix
255 !
256 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
257 nkind = SIZE(qs_kind_set)
258 ALLOCATE (ref_basis_set_list(nkind), orb_basis_set_list(nkind))
259 DO ikind = 1, nkind
260 qs_kind => qs_kind_set(ikind)
261 NULLIFY (ref_basis_set_list(ikind)%gto_basis_set)
262 NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
263 NULLIFY (refbasis, orbbasis)
264 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orbbasis, basis_type="ORB")
265 IF (ASSOCIATED(orbbasis)) orb_basis_set_list(ikind)%gto_basis_set => orbbasis
266 CALL get_qs_kind(qs_kind=qs_kind, basis_set=refbasis, basis_type="MIN")
267 IF (ASSOCIATED(refbasis)) ref_basis_set_list(ikind)%gto_basis_set => refbasis
268 END DO
269
270 ! neighbor lists
271 NULLIFY (srr_list, sro_list)
272 CALL setup_neighbor_list(srr_list, ref_basis_set_list, qs_env=qs_env)
273 CALL setup_neighbor_list(sro_list, ref_basis_set_list, orb_basis_set_list, qs_env=qs_env)
274
275 ! Srr and Sro overlap matrices
276 NULLIFY (sref, sro)
277 CALL get_qs_env(qs_env, ks_env=ks_env)
278 CALL build_overlap_matrix_simple(ks_env, sref, &
279 ref_basis_set_list, ref_basis_set_list, srr_list)
280 CALL build_overlap_matrix_simple(ks_env, sro, &
281 ref_basis_set_list, orb_basis_set_list, sro_list)
282 !
283 IF (PRESENT(mos)) THEN
284 my_mos => mos
285 ELSE
286 CALL get_qs_env(qs_env, mos=my_mos)
287 END IF
288 CALL get_mo_set(my_mos(1), mo_coeff=mo_coeff)
289 CALL dbcsr_get_info(sro(1)%matrix, nfullrows_total=nref, nfullcols_total=nao)
290
291 ! Projectors
292 NULLIFY (fm_struct)
293 CALL cp_fm_struct_create(fm_struct, context=mo_coeff%matrix_struct%context, nrow_global=nref, &
294 ncol_global=nao, para_env=mo_coeff%matrix_struct%para_env)
295 CALL cp_fm_create(p_ref_orb, fm_struct)
296 CALL cp_fm_struct_release(fm_struct)
297 NULLIFY (fm_struct)
298 CALL cp_fm_struct_create(fm_struct, context=mo_coeff%matrix_struct%context, nrow_global=nao, &
299 ncol_global=nref, para_env=mo_coeff%matrix_struct%para_env)
300 CALL cp_fm_create(p_orb_ref, fm_struct)
301 CALL cp_fm_struct_release(fm_struct)
302 !
303 CALL iao_projectors(smat, sref(1)%matrix, sro(1)%matrix, p_orb_ref, p_ref_orb, iao_env%eps_svd)
304
305 ! occupied orbitals, orbitals considered for IAO generation
306 nocc(1:2) = 0
307 DO ispin = 1, nspin
308 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nmo=norb, uniform_occupation=uniform_occupation, &
309 occupation_numbers=occupation_numbers)
310 IF (uniform_occupation) THEN
311 nvec = norb
312 ELSE
313 nvec = 0
314 DO i = 1, norb
315 IF (occupation_numbers(i) > iao_env%eps_occ) THEN
316 nvec = i
317 ELSE
318 EXIT
319 END IF
320 END DO
321 END IF
322 nocc(ispin) = nvec
323 END DO
324 ! generate IAOs
325 ALLOCATE (iao_coef(nspin), cvec(nspin))
326 DO ispin = 1, nspin
327 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff)
328 nvec = nocc(ispin)
329 IF (nvec > 0) THEN
330 NULLIFY (fm_struct)
331 CALL cp_fm_struct_create(fm_struct, ncol_global=nvec, template_fmstruct=mo_coeff%matrix_struct)
332 CALL cp_fm_create(cvec(ispin), fm_struct)
333 CALL cp_fm_struct_release(fm_struct)
334 CALL cp_fm_to_fm(mo_coeff, cvec(ispin), nvec)
335 !
336 NULLIFY (fm_struct)
337 CALL cp_fm_struct_create(fm_struct, ncol_global=nref, template_fmstruct=mo_coeff%matrix_struct)
338 CALL cp_fm_create(iao_coef(ispin), fm_struct)
339 CALL cp_fm_struct_release(fm_struct)
340 !
341 CALL intrinsic_ao_calc(smat, p_orb_ref, p_ref_orb, cvec(ispin), iao_coef(ispin))
342 END IF
343 END DO
344 !
345 IF (iao_env%do_charges) THEN
346 ! MOs in IAO basis
347 ALLOCATE (ciao(nspin))
348 DO ispin = 1, nspin
349 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=norb)
350 CALL cp_fm_create(smo, mo_coeff%matrix_struct)
351 NULLIFY (fm_struct)
352 CALL cp_fm_struct_create(fm_struct, nrow_global=nref, template_fmstruct=mo_coeff%matrix_struct)
353 CALL cp_fm_create(ciao(ispin), fm_struct)
354 CALL cp_fm_struct_release(fm_struct)
355 ! CIAO = A(T)*S*C
356 CALL cp_dbcsr_sm_fm_multiply(smat, mo_coeff, smo, ncol=norb)
357 CALL parallel_gemm('T', 'N', nref, norb, nao, 1.0_dp, iao_coef(ispin), smo, 0.0_dp, ciao(ispin))
358 CALL cp_fm_release(smo)
359 END DO
360 !
361 ! population analysis
362 IF (unit_nr > 0) THEN
363 WRITE (unit_nr, '(/,T2,A)') 'Intrinsic AO Population Analysis '
364 END IF
365 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
366 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
367 ALLOCATE (mcharge(natom, nspin))
368 CALL dbcsr_get_info(sref(1)%matrix, nfullrows_total=nref)
369 DO ispin = 1, nspin
370 CALL get_mo_set(my_mos(ispin), &
371 uniform_occupation=uniform_occupation, &
372 occupation_numbers=occupation_numbers)
373 ! diagonal block matrix
374 CALL dbcsr_create(dmat, template=sref(1)%matrix)
376 !
377 CALL iao_calculate_dmat(ciao(ispin), dmat, occupation_numbers, uniform_occupation)
378 CALL dbcsr_trace(dmat, trace)
379 IF (unit_nr > 0) THEN
380 WRITE (unit_nr, '(T2,A,I2,T66,F15.4)') 'Number of Electrons: Trace(Piao) Spin ', ispin, trace
381 END IF
382 CALL iao_charges(dmat, mcharge(:, ispin))
383 CALL dbcsr_release(dmat)
384 END DO
385 CALL print_atomic_charges(particle_set, qs_kind_set, unit_nr, "Intrinsic Atomic Orbital Charges", &
386 electronic_charges=mcharge)
387 DEALLOCATE (mcharge)
388 CALL cp_fm_release(ciao)
389 END IF
390 !
391 ! Deallocate the neighbor list structure
392 CALL release_neighbor_list_sets(srr_list)
393 CALL release_neighbor_list_sets(sro_list)
394 IF (ASSOCIATED(sref)) CALL dbcsr_deallocate_matrix_set(sref)
395 IF (ASSOCIATED(sro)) CALL dbcsr_deallocate_matrix_set(sro)
396 CALL cp_fm_release(p_ref_orb)
397 CALL cp_fm_release(p_orb_ref)
398 CALL cp_fm_release(cvec)
399 DEALLOCATE (orb_basis_set_list)
400 !
401 IF (iao_env%do_oce) THEN
402 ! generate OCE basis
403 IF (unit_nr > 0) THEN
404 WRITE (unit_nr, '(T2,A)') "IAO One-Center Expansion: OCE Basis Set"
405 END IF
406 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
407 nkind = SIZE(qs_kind_set)
408 ALLOCATE (oce_basis_set_list(nkind), orb_basis_set_list(nkind))
409 DO ikind = 1, nkind
410 qs_kind => qs_kind_set(ikind)
411 NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
412 NULLIFY (orbbasis)
413 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orbbasis, basis_type="ORB")
414 IF (ASSOCIATED(orbbasis)) orb_basis_set_list(ikind)%gto_basis_set => orbbasis
415 END DO
416 DO ikind = 1, nkind
417 orbbasis => orb_basis_set_list(ikind)%gto_basis_set
418 cpassert(ASSOCIATED(orbbasis))
419 NULLIFY (ocebasis)
420 CALL create_oce_basis(ocebasis, orbbasis, iao_env%lmax_oce, iao_env%nbas_oce)
421 CALL init_orb_basis_set(ocebasis)
422 CALL init_interaction_radii_orb_basis(ocebasis, dft_control%qs_control%eps_pgf_orb)
423 oce_basis_set_list(ikind)%gto_basis_set => ocebasis
424 IF (unit_nr > 0) THEN
425 qs_kind => qs_kind_set(ikind)
426 CALL get_qs_kind(qs_kind, zeff=zval, element_symbol=element_symbol)
427 CALL get_gto_basis_set(ocebasis, name=bname, nsgf=nsgf)
428 WRITE (unit_nr, '(T2,A,A,T14,A,I4,T40,A,A30)') "Kind: ", element_symbol, "NBasFun: ", nsgf, &
429 "OCE Basis: ", adjustl(trim(bname))
430 END IF
431 END DO
432 IF (unit_nr > 0) WRITE (unit_nr, *)
433 ! OCE basis overlap
434 ALLOCATE (smat_kind(nkind))
435 DO ikind = 1, nkind
436 ocebasis => oce_basis_set_list(ikind)%gto_basis_set
437 CALL get_gto_basis_set(gto_basis_set=ocebasis, nsgf=nsgf)
438 ALLOCATE (smat_kind(ikind)%array(nsgf, nsgf))
439 END DO
440 CALL oce_overlap_matrix(smat_kind, oce_basis_set_list)
441 ! overlap between IAO OCE basis and orbital basis
442 NULLIFY (sxo, sxo_list)
443 CALL setup_neighbor_list(sxo_list, oce_basis_set_list, orb_basis_set_list, qs_env=qs_env)
444 CALL get_qs_env(qs_env, ks_env=ks_env)
445 CALL build_overlap_matrix_simple(ks_env, sxo, &
446 oce_basis_set_list, orb_basis_set_list, sxo_list)
447 !
448 ! One Center Expansion of IAOs
449 CALL get_qs_env(qs_env=qs_env, natom=natom)
450 ALLOCATE (oce_atom(natom, nspin))
451 CALL iao_oce_expansion(qs_env, oce_basis_set_list, smat_kind, sxo, iao_coef, oce_atom, unit_nr)
452 !
453 DO ikind = 1, nkind
454 ocebasis => oce_basis_set_list(ikind)%gto_basis_set
455 CALL deallocate_gto_basis_set(ocebasis)
456 END DO
457 DEALLOCATE (oce_basis_set_list, orb_basis_set_list)
458 !
459 CALL release_neighbor_list_sets(sxo_list)
460 IF (ASSOCIATED(sxo)) CALL dbcsr_deallocate_matrix_set(sxo)
461 DO ikind = 1, nkind
462 DEALLOCATE (smat_kind(ikind)%array)
463 END DO
464 DEALLOCATE (smat_kind)
465 DO ispin = 1, nspin
466 DO iatom = 1, natom
467 DEALLOCATE (oce_atom(iatom, ispin)%array)
468 END DO
469 END DO
470 DEALLOCATE (oce_atom)
471 END IF
472 !
473 IF (molden_iao) THEN
474 ! Print basis functions: molden file
475 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
476 CALL get_qs_env(qs_env, cell=cell)
477 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
478 IF (unit_nr > 0) THEN
479 WRITE (unit_nr, '(T2,A)') ' Write IAO in MOLDEN Format'
480 END IF
481 ALLOCATE (mo_iao(nspin))
482 DO ispin = 1, nspin
483 CALL get_mo_set(my_mos(ispin), nao=nao)
484 CALL allocate_mo_set(mo_iao(ispin), nao, nref, nref, 0.0_dp, 1.0_dp, 0.0_dp)
485 CALL init_mo_set(mo_iao(ispin), fm_ref=iao_coef(ispin), name="iao_set")
486 CALL cp_fm_to_fm(iao_coef(ispin), mo_iao(ispin)%mo_coeff)
487 CALL set_mo_set(mo_iao(ispin), homo=nref)
488 mo_iao(ispin)%occupation_numbers = 1.0_dp
489 END DO
490 ! Print IAO basis functions: molden format
491 CALL write_mos_molden(mo_iao, qs_kind_set, particle_set, iao_molden_section, cell=cell, qs_env=qs_env)
492 DO ispin = 1, nspin
493 CALL deallocate_mo_set(mo_iao(ispin))
494 END DO
495 DEALLOCATE (mo_iao)
496 END IF
497 IF (cubes_iao) THEN
498 IF (unit_nr > 0) THEN
499 WRITE (unit_nr, '(T2,A)') ' Write IAO as CUBE Files'
500 END IF
501 ! Print basis functions: cube file
502 CALL print_iao_cubes(qs_env, iao_cubes_section, iao_coef, ref_basis_set_list)
503 END IF
504 !
505 ! Intrinsic bond orbitals
506 IF (iao_env%do_bondorbitals) THEN
507 IF (unit_nr > 0) THEN
508 WRITE (unit_nr, '(/,T2,A)') 'Intrinsic Bond Orbital Generation'
509 END IF
510 ! MOs in IAO basis -> ciao
511 ALLOCATE (ciao(nspin))
512 DO ispin = 1, nspin
513 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=norb)
514 CALL cp_fm_create(smo, mo_coeff%matrix_struct)
515 NULLIFY (fm_struct)
516 CALL cp_fm_struct_create(fm_struct, nrow_global=nref, template_fmstruct=mo_coeff%matrix_struct)
517 CALL cp_fm_create(ciao(ispin), fm_struct)
518 CALL cp_fm_struct_release(fm_struct)
519 ! CIAO = A(T)*S*C
520 CALL cp_dbcsr_sm_fm_multiply(smat, mo_coeff, smo, ncol=norb)
521 CALL parallel_gemm('T', 'N', nref, norb, nao, 1.0_dp, iao_coef(ispin), smo, 0.0_dp, ciao(ispin))
522 CALL cp_fm_release(smo)
523 END DO
524 !
525 ! localize occupied orbitals nocc(ispin), see IAO generation
526 !
527 ALLOCATE (cloc(nspin), c_loc_orb(nspin))
528 DO ispin = 1, nspin
529 NULLIFY (fm_struct)
530 CALL cp_fm_struct_create(fm_struct, ncol_global=nocc(ispin), &
531 template_fmstruct=ciao(ispin)%matrix_struct)
532 CALL cp_fm_create(cloc(ispin), fm_struct)
533 CALL cp_fm_struct_release(fm_struct)
534 CALL cp_fm_to_fm(ciao(ispin), cloc(ispin), ncol=nocc(ispin))
535 END DO
536 CALL cp_fm_release(ciao)
537 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
538 ALLOCATE (first_sgf(natom), last_sgf(natom), nsgfat(natom))
539 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
540 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf, &
541 nsgf=nsgfat, basis=ref_basis_set_list)
542
543 IF (iao_env%loc_operator /= do_iaoloc_pm2) THEN
544 cpabort("IAO localization operator NYA")
545 END IF
546 IF (iao_env%eloc_function /= do_iaoloc_enone) THEN
547 cpabort("IAO energy weight localization NYA")
548 END IF
549
550 DO ispin = 1, nspin
551 nvec = nocc(ispin)
552 IF (nvec > 0) THEN
553 ALLOCATE (zij_atom(natom, 1), fdiag(nvec))
554 NULLIFY (fm_struct)
555 CALL cp_fm_struct_create(fm_struct, nrow_global=nvec, ncol_global=nvec, &
556 template_fmstruct=cloc(ispin)%matrix_struct)
557 DO iatom = 1, natom
558 CALL cp_fm_create(zij_atom(iatom, 1), fm_struct)
559 isgf = first_sgf(iatom)
560 CALL parallel_gemm('T', 'N', nvec, nvec, nsgfat(iatom), 1.0_dp, cloc(ispin), cloc(ispin), &
561 0.0_dp, zij_atom(iatom, 1), &
562 a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf)
563 END DO
564 CALL cp_fm_struct_release(fm_struct)
565 !
566 t1 = m_walltime()
567 order = 4
568 fin = 0.0_dp
569 DO iatom = 1, natom
570 CALL cp_fm_get_diag(zij_atom(iatom, 1), fdiag)
571 fin = fin + sum(fdiag**order)
572 END DO
573 fin = fin**(1._dp/order)
574 ! QR localization
575 CALL scdm_qrfact(cloc(ispin))
576 !
577 DO iatom = 1, natom
578 isgf = first_sgf(iatom)
579 CALL parallel_gemm('T', 'N', nvec, nvec, nsgfat(iatom), 1.0_dp, cloc(ispin), cloc(ispin), &
580 0.0_dp, zij_atom(iatom, 1), &
581 a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf)
582 END DO
583 fout = 0.0_dp
584 DO iatom = 1, natom
585 CALL cp_fm_get_diag(zij_atom(iatom, 1), fdiag)
586 fout = fout + sum(fdiag**order)
587 END DO
588 fout = fout**(1._dp/order)
589 DEALLOCATE (fdiag)
590 t2 = m_walltime()
591 IF (unit_nr > 0) THEN
592 WRITE (unit_nr, '(T2,A,F14.8,A,F14.8,A,F8.2)') &
593 'SCDM pre-localization: fin=', fin, " fout=", fout, " Time=", t2 - t1
594 END IF
595 !
596 CALL pm_localization(zij_atom, cloc(ispin), order, 1.e-12_dp, unit_nr)
597 !
598 CALL cp_fm_release(zij_atom)
599 CALL get_mo_set(my_mos(ispin), nao=nao)
600 NULLIFY (fm_struct)
601 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, &
602 template_fmstruct=cloc(ispin)%matrix_struct)
603 CALL cp_fm_create(c_loc_orb(ispin), fm_struct)
604 CALL cp_fm_struct_release(fm_struct)
605 CALL parallel_gemm('N', 'N', nao, nvec, nref, 1.0_dp, iao_coef(ispin), cloc(ispin), &
606 0.0_dp, c_loc_orb(ispin))
607 END IF
608 END DO
609 DEALLOCATE (first_sgf, last_sgf, nsgfat)
610 CALL cp_fm_release(cloc)
611 !
612 IF (iao_env%do_center) THEN
613 IF (unit_nr > 0) THEN
614 WRITE (unit_nr, '(T2,A)') ' Calculate Localized Orbital Centers and Spread'
615 IF (iao_env%pos_periodic) THEN
616 WRITE (unit_nr, '(T2,A)') ' Use Berry Phase Position Operator'
617 ELSE
618 WRITE (unit_nr, '(T2,A)') ' Use Local Position Operator'
619 END IF
620 END IF
621 nvec = maxval(nocc)
622 ! x y z m2(1) m2(2)
623 ALLOCATE (moments(5, nvec, nspin))
624 moments = 0.0_dp
625 IF (iao_env%pos_periodic) THEN
626 CALL center_spread_berry(qs_env, c_loc_orb, moments)
627 ELSE
628 CALL center_spread_loc(qs_env, c_loc_orb, moments)
629 END IF
630 !
631 IF (ASSOCIATED(ibo_cc_section)) THEN
632 CALL print_center_spread(moments, nocc, ibo_cc_section)
633 END IF
634 !
635 IF (PRESENT(bond_centers)) THEN
636 nx = SIZE(bond_centers, 1)
637 no = SIZE(bond_centers, 2)
638 ns = SIZE(bond_centers, 3)
639 cpassert(no >= nvec)
640 cpassert(ns == nspin)
641 cpassert(nx >= 3)
642 nx = min(nx, 5)
643 bond_centers(1:nx, 1:nvec, 1:ns) = moments(1:nx, 1:nvec, 1:ns)
644 END IF
645 DEALLOCATE (moments)
646 END IF
647 !
648 IF (molden_ibo) THEN
649 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
650 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
651 IF (unit_nr > 0) THEN
652 WRITE (unit_nr, '(T2,A)') ' Write IBO in MOLDEN Format'
653 END IF
654 ALLOCATE (mo_loc(nspin))
655 DO ispin = 1, nspin
656 CALL get_mo_set(my_mos(ispin), nao=nao)
657 nvec = nocc(ispin)
658 CALL allocate_mo_set(mo_loc(ispin), nao, nvec, nvec, 0.0_dp, 1.0_dp, 0.0_dp)
659 CALL init_mo_set(mo_loc(ispin), fm_ref=c_loc_orb(ispin), name="ibo_orb")
660 CALL cp_fm_to_fm(c_loc_orb(ispin), mo_loc(ispin)%mo_coeff)
661 CALL set_mo_set(mo_loc(ispin), homo=nvec)
662 mo_loc(ispin)%occupation_numbers = 1.0_dp
663 END DO
664 ! Print IBO functions: molden format
665 CALL write_mos_molden(mo_loc, qs_kind_set, particle_set, ibo_molden_section, cell=cell, qs_env=qs_env)
666 DO ispin = 1, nspin
667 CALL deallocate_mo_set(mo_loc(ispin))
668 END DO
669 DEALLOCATE (mo_loc)
670 END IF
671 !
672 IF (cubes_ibo) THEN
673 IF (unit_nr > 0) THEN
674 WRITE (unit_nr, '(T2,A)') ' Write IBO on CUBE files'
675 END IF
676 ! Print localized orbital functions: cube file
677 CALL print_ibo_cubes(qs_env, ibo_cubes_section, c_loc_orb)
678 END IF
679 !
680 IF (PRESENT(mos) .AND. ALLOCATED(c_loc_orb)) THEN
681 ! c_loc_orb -> mos
682 DO ispin = 1, nspin
683 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
684 nvec = nocc(ispin)
685 CALL cp_fm_to_fm(c_loc_orb(ispin), mo_coeff, ncol=nvec)
686 END DO
687 END IF
688 CALL cp_fm_release(c_loc_orb)
689 END IF
690 DEALLOCATE (ref_basis_set_list)
691
692 IF (PRESENT(c_iao_coef)) THEN
693 CALL cp_fm_release(c_iao_coef)
694 ALLOCATE (c_iao_coef(nspin))
695 DO ispin = 1, nspin
696 CALL cp_fm_create(c_iao_coef(ispin), iao_coef(ispin)%matrix_struct)
697 CALL cp_fm_to_fm(iao_coef(ispin), c_iao_coef(ispin))
698 END DO
699 END IF
700 CALL cp_fm_release(iao_coef)
701
702 IF (unit_nr > 0) THEN
703 WRITE (unit_nr, '(T2,A)') &
704 '!----------------------------END OF IAO ANALYSIS------------------------------!'
705 END IF
706
707 CALL timestop(handle)
708
709 END SUBROUTINE iao_wfn_analysis
710
711! **************************************************************************************************
712!> \brief Computes projector matrices for ref to orb basis and reverse
713!> \param smat ...
714!> \param sref ...
715!> \param s_r_o ...
716!> \param p_o_r ...
717!> \param p_r_o ...
718!> \param eps_svd ...
719! **************************************************************************************************
720 SUBROUTINE iao_projectors(smat, sref, s_r_o, p_o_r, p_r_o, eps_svd)
721 TYPE(dbcsr_type), INTENT(INOUT) :: smat, sref, s_r_o
722 TYPE(cp_fm_type), INTENT(INOUT) :: p_o_r, p_r_o
723 REAL(kind=dp), INTENT(IN) :: eps_svd
724
725 CHARACTER(len=*), PARAMETER :: routinen = 'iao_projectors'
726
727 INTEGER :: handle, norb, nref
728 TYPE(cp_fm_struct_type), POINTER :: fm_struct
729 TYPE(cp_fm_type) :: fm_inv, fm_s, fm_sro
730
731 CALL timeset(routinen, handle)
732
733 CALL cp_fm_get_info(p_r_o, nrow_global=nref, ncol_global=norb)
734 CALL cp_fm_create(fm_sro, p_r_o%matrix_struct)
735 CALL copy_dbcsr_to_fm(s_r_o, fm_sro)
736
737 NULLIFY (fm_struct)
738 CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=norb, &
739 template_fmstruct=p_r_o%matrix_struct)
740 CALL cp_fm_create(fm_s, fm_struct, name="smat")
741 CALL cp_fm_create(fm_inv, fm_struct, name="sinv")
742 CALL cp_fm_struct_release(fm_struct)
743 CALL copy_dbcsr_to_fm(smat, fm_s)
744 CALL cp_fm_invert(fm_s, fm_inv, eps_svd=eps_svd)
745 CALL parallel_gemm('N', 'T', norb, nref, norb, 1.0_dp, fm_inv, fm_sro, 0.0_dp, p_o_r)
746 CALL cp_fm_release(fm_s)
747 CALL cp_fm_release(fm_inv)
748
749 NULLIFY (fm_struct)
750 CALL cp_fm_struct_create(fm_struct, nrow_global=nref, ncol_global=nref, &
751 template_fmstruct=p_r_o%matrix_struct)
752 CALL cp_fm_create(fm_s, fm_struct, name="sref")
753 CALL cp_fm_create(fm_inv, fm_struct, name="sinv")
754 CALL cp_fm_struct_release(fm_struct)
755 CALL copy_dbcsr_to_fm(sref, fm_s)
756 CALL cp_fm_invert(fm_s, fm_inv, eps_svd=eps_svd)
757 CALL parallel_gemm('N', 'N', nref, norb, nref, 1.0_dp, fm_inv, fm_sro, 0.0_dp, p_r_o)
758 CALL cp_fm_release(fm_s)
759 CALL cp_fm_release(fm_inv)
760
761 CALL cp_fm_release(fm_sro)
762
763 CALL timestop(handle)
764
765 END SUBROUTINE iao_projectors
766
767! **************************************************************************************************
768!> \brief Computes intrinsic orbitals for a given MO vector set
769!> \param smat ...
770!> \param p_o_r ...
771!> \param p_r_o ...
772!> \param cvec ...
773!> \param avec ...
774! **************************************************************************************************
775 SUBROUTINE intrinsic_ao_calc(smat, p_o_r, p_r_o, cvec, avec)
776 TYPE(dbcsr_type), INTENT(INOUT) :: smat
777 TYPE(cp_fm_type), INTENT(INOUT) :: p_o_r, p_r_o, cvec, avec
778
779 CHARACTER(len=*), PARAMETER :: routinen = 'intrinsic_ao_calc'
780
781 INTEGER :: handle, nao, norb, nref
782 TYPE(cp_fm_struct_type), POINTER :: fm_struct
783 TYPE(cp_fm_type) :: ctvec, pc, sc, sct, vec1
784
785 CALL timeset(routinen, handle)
786
787 ! number of orbitals
788 CALL cp_fm_get_info(cvec, ncol_global=norb)
789 ! basis set sizes
790 CALL cp_fm_get_info(p_r_o, nrow_global=nref, ncol_global=nao)
791 ! temp array
792 NULLIFY (fm_struct)
793 CALL cp_fm_struct_create(fm_struct, nrow_global=nref, ncol_global=norb, &
794 template_fmstruct=cvec%matrix_struct)
795 CALL cp_fm_create(pc, fm_struct)
796 CALL cp_fm_struct_release(fm_struct)
797 ! CT = orth(Por.Pro.C)
798 CALL parallel_gemm('N', 'N', nref, norb, nao, 1.0_dp, p_r_o, cvec, 0.0_dp, pc)
799 CALL cp_fm_create(ctvec, cvec%matrix_struct)
800 CALL parallel_gemm('N', 'N', nao, norb, nref, 1.0_dp, p_o_r, pc, 0.0_dp, ctvec)
801 CALL cp_fm_release(pc)
802 CALL make_basis_lowdin(ctvec, norb, smat)
803 ! S*C and S*CT
804 CALL cp_fm_create(sc, cvec%matrix_struct)
805 CALL cp_dbcsr_sm_fm_multiply(smat, cvec, sc, ncol=norb)
806 CALL cp_fm_create(sct, cvec%matrix_struct)
807 CALL cp_dbcsr_sm_fm_multiply(smat, ctvec, sct, ncol=norb)
808 CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=nref, &
809 template_fmstruct=cvec%matrix_struct)
810 CALL cp_fm_create(pc, fm_struct)
811 CALL cp_fm_struct_release(fm_struct)
812 ! V1 = (CT*SCT(T))Por
813 CALL parallel_gemm('T', 'N', norb, nref, nao, 1.0_dp, sct, p_o_r, 0.0_dp, pc)
814 NULLIFY (fm_struct)
815 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nref, &
816 template_fmstruct=cvec%matrix_struct)
817 CALL cp_fm_create(vec1, fm_struct)
818 CALL cp_fm_struct_release(fm_struct)
819 CALL parallel_gemm('N', 'N', nao, nref, norb, 1.0_dp, ctvec, pc, 0.0_dp, vec1)
820 ! A = C*SC(T)*V1
821 CALL parallel_gemm('T', 'N', norb, nref, nao, 1.0_dp, sc, vec1, 0.0_dp, pc)
822 CALL parallel_gemm('N', 'N', nao, nref, norb, 1.0_dp, cvec, pc, 0.0_dp, avec)
823 ! V1 = Por - V1
824 CALL cp_fm_geadd(1.0_dp, 'N', p_o_r, -1.0_dp, vec1)
825 ! A = A + V1
826 CALL cp_fm_geadd(1.0_dp, 'N', vec1, 1.0_dp, avec)
827 ! A = A - C*SC(T)*V1
828 CALL parallel_gemm('T', 'N', norb, nref, nao, 1.0_dp, sc, vec1, 0.0_dp, pc)
829 CALL parallel_gemm('N', 'N', nao, nref, norb, -1.0_dp, cvec, pc, 1.0_dp, avec)
830 ! A = orth(A)
831 CALL make_basis_lowdin(avec, nref, smat)
832
833 ! clean up
834 CALL cp_fm_release(pc)
835 CALL cp_fm_release(vec1)
836 CALL cp_fm_release(sc)
837 CALL cp_fm_release(sct)
838 CALL cp_fm_release(ctvec)
839
840 CALL timestop(handle)
841
842 END SUBROUTINE intrinsic_ao_calc
843
844! **************************************************************************************************
845!> \brief Calculate the density matrix from fm vectors using occupation numbers
846!> \param cvec ...
847!> \param density_matrix ...
848!> \param occupation ...
849!> \param uniform_occupation ...
850! **************************************************************************************************
851 SUBROUTINE iao_calculate_dmat_diag(cvec, density_matrix, occupation, uniform_occupation)
852
853 TYPE(cp_fm_type), INTENT(IN) :: cvec
854 TYPE(dbcsr_type) :: density_matrix
855 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
856 LOGICAL, INTENT(IN) :: uniform_occupation
857
858 CHARACTER(len=*), PARAMETER :: routineN = 'iao_calculate_dmat_diag'
859
860 INTEGER :: handle, ncol
861 REAL(KIND=dp) :: alpha
862 TYPE(cp_fm_type) :: fm_tmp
863
864 CALL timeset(routinen, handle)
865
866 CALL dbcsr_set(density_matrix, 0.0_dp)
867
868 CALL cp_fm_get_info(cvec, ncol_global=ncol)
869 IF (.NOT. uniform_occupation) THEN
870 CALL cp_fm_create(fm_tmp, cvec%matrix_struct)
871 CALL cp_fm_to_fm(cvec, fm_tmp)
872 CALL cp_fm_column_scale(fm_tmp, occupation(1:ncol))
873 alpha = 1.0_dp
874 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
875 matrix_v=cvec, matrix_g=fm_tmp, &
876 ncol=ncol, alpha=alpha)
877 CALL cp_fm_release(fm_tmp)
878 ELSE
879 alpha = occupation(1)
880 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
881 matrix_v=cvec, ncol=ncol, alpha=alpha)
882 END IF
883
884 CALL timestop(handle)
885
886 END SUBROUTINE iao_calculate_dmat_diag
887
888! **************************************************************************************************
889!> \brief Calculate the density matrix from fm vectors using an occupation matrix
890!> \param cvec ...
891!> \param density_matrix ...
892!> \param weight ...
893!> \param occmat ...
894! **************************************************************************************************
895 SUBROUTINE iao_calculate_dmat_full(cvec, density_matrix, weight, occmat)
896
897 TYPE(cp_fm_type), INTENT(IN) :: cvec
898 TYPE(dbcsr_type) :: density_matrix
899 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: weight
900 TYPE(cp_fm_type), INTENT(IN) :: occmat
901
902 CHARACTER(len=*), PARAMETER :: routinen = 'iao_calculate_dmat_full'
903
904 INTEGER :: handle, ic, jc, ncol
905 REAL(kind=dp) :: alpha
906 TYPE(cp_fm_struct_type), POINTER :: fm_struct
907 TYPE(cp_fm_type) :: fm1, fm2
908
909 CALL timeset(routinen, handle)
910
911 CALL dbcsr_set(density_matrix, 0.0_dp)
912
913 CALL cp_fm_struct_create(fm_struct, ncol_global=1, &
914 template_fmstruct=cvec%matrix_struct)
915 CALL cp_fm_create(fm1, fm_struct)
916 CALL cp_fm_create(fm2, fm_struct)
917 CALL cp_fm_struct_release(fm_struct)
918
919 CALL cp_fm_get_info(cvec, ncol_global=ncol)
920 DO ic = 1, ncol
921 CALL cp_fm_to_fm(cvec, fm1, ncol=1, source_start=ic, target_start=1)
922 DO jc = 1, ncol
923 CALL cp_fm_to_fm(cvec, fm2, ncol=1, source_start=jc, target_start=1)
924 CALL cp_fm_get_element(occmat, ic, jc, alpha)
925 alpha = weight(ic)*alpha
926 CALL cp_dbcsr_plus_fm_fm_t(density_matrix, fm1, fm2, ncol=1, &
927 alpha=alpha, symmetry_mode=1)
928 END DO
929 END DO
930 CALL cp_fm_release(fm1)
931 CALL cp_fm_release(fm2)
932
933 CALL timestop(handle)
934
935 END SUBROUTINE iao_calculate_dmat_full
936
937! **************************************************************************************************
938!> \brief compute the atomic charges (orthogonal basis)
939!> \param p_matrix ...
940!> \param charges ...
941! **************************************************************************************************
942 SUBROUTINE iao_charges(p_matrix, charges)
943 TYPE(dbcsr_type) :: p_matrix
944 REAL(kind=dp), DIMENSION(:) :: charges
945
946 INTEGER :: i, iblock_col, iblock_row
947 REAL(kind=dp) :: trace
948 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block
949 TYPE(dbcsr_iterator_type) :: iter
950 TYPE(mp_comm_type) :: group
951
952 charges = 0.0_dp
953
954 CALL dbcsr_iterator_start(iter, p_matrix)
955 DO WHILE (dbcsr_iterator_blocks_left(iter))
956 NULLIFY (p_block)
957 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, p_block)
958 cpassert(iblock_row == iblock_col)
959 trace = 0.0_dp
960 DO i = 1, SIZE(p_block, 1)
961 trace = trace + p_block(i, i)
962 END DO
963 charges(iblock_row) = trace
964 END DO
965 CALL dbcsr_iterator_stop(iter)
966
967 CALL dbcsr_get_info(p_matrix, group=group)
968 CALL group%sum(charges)
969
970 END SUBROUTINE iao_charges
971
972! **************************************************************************************************
973!> \brief Computes and prints the Cube Files for the intrinsic atomic orbitals
974!> \param qs_env ...
975!> \param print_section ...
976!> \param iao_coef ...
977!> \param basis_set_list ...
978! **************************************************************************************************
979 SUBROUTINE print_iao_cubes(qs_env, print_section, iao_coef, basis_set_list)
980 TYPE(qs_environment_type), POINTER :: qs_env
981 TYPE(section_vals_type), POINTER :: print_section
982 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: iao_coef
983 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
984
985 CHARACTER(LEN=default_path_length) :: filename, title
986 INTEGER :: i, i_rep, ispin, ivec, iw, j, n_rep, &
987 nat, natom, norb, nspins, nstart
988 INTEGER, DIMENSION(:), POINTER :: atom_list, blk_sizes, first_bas, stride
989 LOGICAL :: mpi_io
990 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
991 TYPE(cell_type), POINTER :: cell
992 TYPE(cp_logger_type), POINTER :: logger
993 TYPE(dft_control_type), POINTER :: dft_control
994 TYPE(particle_list_type), POINTER :: particles
995 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
996 TYPE(pw_c1d_gs_type) :: wf_g
997 TYPE(pw_env_type), POINTER :: pw_env
998 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
999 TYPE(pw_r3d_rs_type) :: wf_r
1000 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1001 TYPE(qs_subsys_type), POINTER :: subsys
1002
1003 logger => cp_get_default_logger()
1004 stride => section_get_ivals(print_section, "STRIDE")
1005 CALL section_vals_val_get(print_section, "ATOM_LIST", n_rep_val=n_rep)
1006
1007 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1008 subsys=subsys, cell=cell, particle_set=particle_set, pw_env=pw_env, dft_control=dft_control)
1009 CALL qs_subsys_get(subsys, particles=particles)
1010
1011 CALL get_qs_env(qs_env=qs_env, natom=natom)
1012 ALLOCATE (blk_sizes(natom), first_bas(0:natom))
1013 CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_sizes, basis=basis_set_list)
1014 first_bas(0) = 0
1015 DO i = 1, natom
1016 first_bas(i) = first_bas(i - 1) + blk_sizes(i)
1017 END DO
1018
1019 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1020 CALL auxbas_pw_pool%create_pw(wf_r)
1021 CALL auxbas_pw_pool%create_pw(wf_g)
1022
1023 nspins = SIZE(iao_coef)
1024 nstart = min(1, n_rep)
1025
1026 DO ispin = 1, nspins
1027 DO i_rep = nstart, n_rep
1028 CALL cp_fm_get_info(iao_coef(ispin), ncol_global=norb)
1029 IF (i_rep == 0) THEN
1030 nat = natom
1031 ELSE
1032 CALL section_vals_val_get(print_section, "ATOM_LIST", i_rep_val=i_rep, i_vals=atom_list)
1033 nat = SIZE(atom_list)
1034 END IF
1035 DO i = 1, nat
1036 IF (i_rep == 0) THEN
1037 j = i
1038 ELSE
1039 j = atom_list(i)
1040 END IF
1041 cpassert(j >= 1 .AND. j <= natom)
1042 DO ivec = first_bas(j - 1) + 1, first_bas(j)
1043 WRITE (filename, '(a4,I5.5,a1,I1.1)') "IAO_", ivec, "_", ispin
1044 WRITE (title, *) "Intrinsic Atomic Orbitals ", ivec, " atom ", j, " spin ", ispin
1045 mpi_io = .true.
1046 iw = cp_print_key_unit_nr(logger, print_section, "", extension=".cube", &
1047 middle_name=trim(filename), file_position="REWIND", log_filename=.false., &
1048 mpi_io=mpi_io)
1049 CALL calculate_wavefunction(iao_coef(ispin), ivec, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1050 cell, dft_control, particle_set, pw_env)
1051 CALL cp_pw_to_cube(wf_r, iw, title, particles=particles, stride=stride, mpi_io=mpi_io)
1052 CALL cp_print_key_finished_output(iw, logger, print_section, "", mpi_io=mpi_io)
1053 END DO
1054 END DO
1055 END DO
1056 END DO
1057 DEALLOCATE (blk_sizes, first_bas)
1058 CALL auxbas_pw_pool%give_back_pw(wf_r)
1059 CALL auxbas_pw_pool%give_back_pw(wf_g)
1060
1061 END SUBROUTINE print_iao_cubes
1062
1063! **************************************************************************************************
1064!> \brief Computes and prints the Cube Files for the intrinsic bond orbitals
1065!> \param qs_env ...
1066!> \param print_section ...
1067!> \param ibo_coef ...
1068! **************************************************************************************************
1069 SUBROUTINE print_ibo_cubes(qs_env, print_section, ibo_coef)
1070 TYPE(qs_environment_type), POINTER :: qs_env
1071 TYPE(section_vals_type), POINTER :: print_section
1072 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: ibo_coef
1073
1074 CHARACTER(LEN=default_path_length) :: filename, title
1075 INTEGER :: i, i_rep, ispin, iw, j, n_rep, norb, &
1076 nspins, nstart, nstate
1077 INTEGER, DIMENSION(:), POINTER :: state_list, stride
1078 LOGICAL :: mpi_io
1079 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1080 TYPE(cell_type), POINTER :: cell
1081 TYPE(cp_logger_type), POINTER :: logger
1082 TYPE(dft_control_type), POINTER :: dft_control
1083 TYPE(particle_list_type), POINTER :: particles
1084 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1085 TYPE(pw_c1d_gs_type) :: wf_g
1086 TYPE(pw_env_type), POINTER :: pw_env
1087 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1088 TYPE(pw_r3d_rs_type) :: wf_r
1089 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1090 TYPE(qs_subsys_type), POINTER :: subsys
1091
1092 cpassert(ASSOCIATED(print_section))
1093
1094 logger => cp_get_default_logger()
1095 stride => section_get_ivals(print_section, "STRIDE")
1096 CALL section_vals_val_get(print_section, "STATE_LIST", n_rep_val=n_rep)
1097
1098 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1099 subsys=subsys, cell=cell, particle_set=particle_set, pw_env=pw_env, dft_control=dft_control)
1100 CALL qs_subsys_get(subsys, particles=particles)
1101
1102 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1103 CALL auxbas_pw_pool%create_pw(wf_r)
1104 CALL auxbas_pw_pool%create_pw(wf_g)
1105
1106 nspins = SIZE(ibo_coef)
1107 nstart = min(1, n_rep)
1108
1109 DO ispin = 1, nspins
1110 DO i_rep = nstart, n_rep
1111 CALL cp_fm_get_info(ibo_coef(ispin), ncol_global=norb)
1112 IF (i_rep == 0) THEN
1113 nstate = norb
1114 ELSE
1115 CALL section_vals_val_get(print_section, "STATE_LIST", i_rep_val=i_rep, i_vals=state_list)
1116 nstate = SIZE(state_list)
1117 END IF
1118 DO i = 1, nstate
1119 IF (i_rep == 0) THEN
1120 j = i
1121 ELSE
1122 j = state_list(i)
1123 END IF
1124 cpassert(j >= 1 .AND. j <= norb)
1125 WRITE (filename, '(a4,I5.5,a1,I1.1)') "IBO_", j, "_", ispin
1126 WRITE (title, *) "Intrinsic Atomic Orbitals ", j, " spin ", ispin
1127 mpi_io = .true.
1128 iw = cp_print_key_unit_nr(logger, print_section, "", extension=".cube", &
1129 middle_name=trim(filename), file_position="REWIND", log_filename=.false., &
1130 mpi_io=mpi_io)
1131 CALL calculate_wavefunction(ibo_coef(ispin), j, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1132 cell, dft_control, particle_set, pw_env)
1133 CALL cp_pw_to_cube(wf_r, iw, title, particles=particles, stride=stride, mpi_io=mpi_io)
1134 CALL cp_print_key_finished_output(iw, logger, print_section, "", mpi_io=mpi_io)
1135 END DO
1136 END DO
1137 END DO
1138
1139 CALL auxbas_pw_pool%give_back_pw(wf_r)
1140 CALL auxbas_pw_pool%give_back_pw(wf_g)
1141
1142 END SUBROUTINE print_ibo_cubes
1143
1144! **************************************************************************************************
1145!> \brief prints charge center and spreads for all orbitals
1146!> \param moments ...
1147!> \param nocc ...
1148!> \param print_section ...
1149!> \param iounit ...
1150! **************************************************************************************************
1151 SUBROUTINE print_center_spread(moments, nocc, print_section, iounit)
1152 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: moments
1153 INTEGER, DIMENSION(:), INTENT(IN) :: nocc
1154 TYPE(section_vals_type), OPTIONAL, POINTER :: print_section
1155 INTEGER, INTENT(IN), OPTIONAL :: iounit
1156
1157 CHARACTER(LEN=default_path_length) :: filename
1158 INTEGER :: is, ispin, iw, nspin
1159 TYPE(cp_logger_type), POINTER :: logger
1160
1161 logger => cp_get_default_logger()
1162 nspin = SIZE(moments, 3)
1163
1164 IF (PRESENT(print_section)) THEN
1165 WRITE (filename, '(A18,I1.1)') "IBO_CENTERS_SPREAD"
1166 iw = cp_print_key_unit_nr(logger, print_section, "", extension=".csp", &
1167 middle_name=trim(filename), file_position="REWIND", log_filename=.false.)
1168 ELSEIF (PRESENT(iounit)) THEN
1169 iw = iounit
1170 ELSE
1171 iw = -1
1172 END IF
1173 IF (iw > 0) THEN
1174 DO ispin = 1, nspin
1175 WRITE (iw, "(T2,A,i1)") "Intrinsic Bond Orbitals: Centers and Spread for Spin ", ispin
1176 WRITE (iw, "(A7,T30,A6,T68,A7)") "State", "Center", "Spreads"
1177 DO is = 1, nocc(ispin)
1178 WRITE (iw, "(i7,3F15.8,8X,2F10.5)") is, moments(:, is, ispin)
1179 END DO
1180 END DO
1181 END IF
1182 IF (PRESENT(print_section)) THEN
1183 CALL cp_print_key_finished_output(iw, logger, print_section, "")
1184 END IF
1185
1186 END SUBROUTINE print_center_spread
1187
1188! **************************************************************************************************
1189!> \brief ...
1190!> \param qs_env ...
1191!> \param c_loc_orb ...
1192!> \param moments ...
1193! **************************************************************************************************
1194 SUBROUTINE center_spread_loc(qs_env, c_loc_orb, moments)
1195 TYPE(qs_environment_type), POINTER :: qs_env
1196 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: c_loc_orb
1197 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: moments
1198
1199 CHARACTER(len=*), PARAMETER :: routinen = 'center_spread_loc'
1200 INTEGER, DIMENSION(6), PARAMETER :: list = [1, 2, 3, 4, 7, 9]
1201
1202 INTEGER :: handle, i, iop, iorb, ispin, nao, norb, &
1203 nspin
1204 REAL(kind=dp) :: xmii
1205 REAL(kind=dp), DIMENSION(3) :: rpoint
1206 TYPE(cell_type), POINTER :: cell
1207 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1208 TYPE(cp_fm_type) :: ccmat, ocvec
1209 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat
1210 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1211 TYPE(dbcsr_type), POINTER :: omat, smat
1212
1213 CALL timeset(routinen, handle)
1214
1215 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
1216 smat => matrix_s(1, 1)%matrix
1217 rpoint = 0.0_dp
1218 nspin = SIZE(c_loc_orb)
1219 moments = 0.0_dp
1220
1221 ALLOCATE (dipmat(9))
1222 DO i = 1, 9
1223 ALLOCATE (dipmat(i)%matrix)
1224 CALL dbcsr_copy(dipmat(i)%matrix, smat)
1225 CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
1226 END DO
1227
1228 CALL build_local_moment_matrix(qs_env, dipmat, 2, rpoint)
1229
1230 DO ispin = 1, nspin
1231 CALL cp_fm_get_info(c_loc_orb(ispin), nrow_global=nao, ncol_global=norb)
1232 CALL cp_fm_create(ocvec, c_loc_orb(ispin)%matrix_struct)
1233 NULLIFY (fm_struct)
1234 CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=norb, &
1235 template_fmstruct=c_loc_orb(ispin)%matrix_struct)
1236 CALL cp_fm_create(ccmat, fm_struct)
1237 CALL cp_fm_struct_release(fm_struct)
1238 DO i = 1, 6
1239 iop = list(i)
1240 omat => dipmat(iop)%matrix
1241 CALL cp_dbcsr_sm_fm_multiply(omat, c_loc_orb(ispin), ocvec, ncol=norb)
1242 CALL parallel_gemm('T', 'N', norb, norb, nao, 1.0_dp, c_loc_orb(ispin), ocvec, &
1243 0.0_dp, ccmat)
1244 DO iorb = 1, norb
1245 CALL cp_fm_get_element(ccmat, iorb, iorb, xmii)
1246 IF (iop <= 3) THEN
1247 moments(iop, iorb, ispin) = moments(iop, iorb, ispin) + xmii
1248 moments(4, iorb, ispin) = moments(4, iorb, ispin) - xmii**2
1249 ELSE
1250 moments(4, iorb, ispin) = moments(4, iorb, ispin) + xmii
1251 END IF
1252 END DO
1253 END DO
1254 CALL cp_fm_release(ocvec)
1255 CALL cp_fm_release(ccmat)
1256 END DO
1257
1258 DO i = 1, 9
1259 CALL dbcsr_release(dipmat(i)%matrix)
1260 DEALLOCATE (dipmat(i)%matrix)
1261 END DO
1262 DEALLOCATE (dipmat)
1263
1264 CALL get_qs_env(qs_env=qs_env, cell=cell)
1265 DO ispin = 1, nspin
1266 DO iorb = 1, norb
1267 moments(1:3, iorb, ispin) = pbc(moments(1:3, iorb, ispin), cell)
1268 END DO
1269 END DO
1270
1271 CALL timestop(handle)
1272
1273 END SUBROUTINE center_spread_loc
1274
1275! **************************************************************************************************
1276!> \brief ...
1277!> \param qs_env ...
1278!> \param c_loc_orb ...
1279!> \param moments ...
1280! **************************************************************************************************
1281 SUBROUTINE center_spread_berry(qs_env, c_loc_orb, moments)
1282 TYPE(qs_environment_type), POINTER :: qs_env
1283 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: c_loc_orb
1284 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: moments
1285
1286 CHARACTER(len=*), PARAMETER :: routinen = 'center_spread_berry'
1287
1288 COMPLEX(KIND=dp) :: z
1289 INTEGER :: dim_op, handle, i, idir, ispin, istate, &
1290 j, jdir, nao, norb, nspin
1291 REAL(dp), DIMENSION(3) :: c, cpbc
1292 REAL(dp), DIMENSION(6) :: weights
1293 REAL(kind=dp) :: imagpart, realpart, spread_i, spread_ii
1294 TYPE(cell_type), POINTER :: cell
1295 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1296 TYPE(cp_fm_type) :: opvec
1297 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: zij
1298 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1299 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
1300
1301 CALL timeset(routinen, handle)
1302
1303 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, cell=cell)
1304
1305 IF (cell%orthorhombic) THEN
1306 dim_op = 3
1307 ELSE
1308 dim_op = 6
1309 END IF
1310 ALLOCATE (op_sm_set(2, dim_op))
1311 DO i = 1, dim_op
1312 DO j = 1, 2
1313 NULLIFY (op_sm_set(j, i)%matrix)
1314 ALLOCATE (op_sm_set(j, i)%matrix)
1315 CALL dbcsr_copy(op_sm_set(j, i)%matrix, matrix_s(1)%matrix)
1316 CALL dbcsr_set(op_sm_set(j, i)%matrix, 0.0_dp)
1317 END DO
1318 END DO
1319 CALL initialize_weights(cell, weights)
1320
1321 CALL compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
1322
1323 nspin = SIZE(c_loc_orb, 1)
1324 DO ispin = 1, nspin
1325 CALL cp_fm_get_info(c_loc_orb(ispin), nrow_global=nao, ncol_global=norb)
1326 CALL cp_fm_create(opvec, c_loc_orb(ispin)%matrix_struct)
1327 NULLIFY (fm_struct)
1328 CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=norb, &
1329 template_fmstruct=c_loc_orb(ispin)%matrix_struct)
1330 ALLOCATE (zij(2, dim_op))
1331
1332 ! Compute zij here
1333 DO i = 1, dim_op
1334 DO j = 1, 2
1335 CALL cp_fm_create(zij(j, i), fm_struct)
1336 CALL cp_fm_set_all(zij(j, i), 0.0_dp)
1337 CALL cp_dbcsr_sm_fm_multiply(op_sm_set(j, i)%matrix, c_loc_orb(ispin), opvec, ncol=norb)
1338 CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, c_loc_orb(ispin), opvec, 0.0_dp, zij(j, i))
1339 END DO
1340 END DO
1341
1342 CALL cp_fm_struct_release(fm_struct)
1343 CALL cp_fm_release(opvec)
1344
1345 DO istate = 1, norb
1346 c = 0.0_dp
1347 spread_i = 0.0_dp
1348 spread_ii = 0.0_dp
1349 DO jdir = 1, dim_op
1350 CALL cp_fm_get_element(zij(1, jdir), istate, istate, realpart)
1351 CALL cp_fm_get_element(zij(2, jdir), istate, istate, imagpart)
1352 z = cmplx(realpart, imagpart, dp)
1353 spread_i = spread_i - weights(jdir)* &
1354 log(realpart*realpart + imagpart*imagpart)/twopi/twopi
1355 spread_ii = spread_ii + weights(jdir)* &
1356 (1.0_dp - (realpart*realpart + imagpart*imagpart))/twopi/twopi
1357 IF (jdir < 4) THEN
1358 DO idir = 1, 3
1359 c(idir) = c(idir) + &
1360 (cell%hmat(idir, jdir)/twopi)*aimag(log(z))
1361 END DO
1362 END IF
1363 END DO
1364 cpbc = pbc(c, cell)
1365 moments(1:3, istate, ispin) = cpbc(1:3)
1366 moments(4, istate, ispin) = spread_i
1367 moments(5, istate, ispin) = spread_ii
1368 END DO
1369
1370 CALL cp_fm_release(zij)
1371
1372 END DO
1373
1374 DO i = 1, dim_op
1375 DO j = 1, 2
1376 CALL dbcsr_release(op_sm_set(j, i)%matrix)
1377 DEALLOCATE (op_sm_set(j, i)%matrix)
1378 END DO
1379 END DO
1380 DEALLOCATE (op_sm_set)
1381
1382 CALL timestop(handle)
1383
1384 END SUBROUTINE center_spread_berry
1385
1386! **************************************************************************************************
1387!> \brief ...
1388!> \param qs_env ...
1389!> \param oce_basis_set_list ...
1390!> \param smat_kind ...
1391!> \param sxo ...
1392!> \param iao_coef ...
1393!> \param oce_atom ...
1394!> \param unit_nr ...
1395! **************************************************************************************************
1396 SUBROUTINE iao_oce_expansion(qs_env, oce_basis_set_list, smat_kind, sxo, iao_coef, oce_atom, unit_nr)
1397 TYPE(qs_environment_type), POINTER :: qs_env
1398 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: oce_basis_set_list
1399 TYPE(cp_2d_r_p_type), DIMENSION(:) :: smat_kind
1400 TYPE(dbcsr_p_type), DIMENSION(:) :: sxo
1401 TYPE(cp_fm_type), DIMENSION(:) :: iao_coef
1402 TYPE(cp_2d_r_p_type), DIMENSION(:, :) :: oce_atom
1403 INTEGER, INTENT(IN) :: unit_nr
1404
1405 CHARACTER(len=*), PARAMETER :: routinen = 'iao_oce_expansion'
1406
1407 INTEGER :: handle, i, i1, i2, iao, iatom, ikind, &
1408 iset, ishell, ispin, l, m, maxl, n, &
1409 natom, nkind, noce, ns, nset, nsgf, &
1410 nspin
1411 INTEGER, DIMENSION(:), POINTER :: iao_blk_sizes, nshell, oce_blk_sizes, &
1412 orb_blk_sizes
1413 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, last_sgf, lval
1414 LOGICAL :: found
1415 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ev, vector
1416 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat, filter, oce_comp, prol, vec
1417 REAL(kind=dp), DIMENSION(:, :), POINTER :: ablock
1418 TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: sinv_kind
1419 TYPE(dbcsr_distribution_type) :: dbcsr_dist
1420 TYPE(dbcsr_type) :: iao_vec, sx_vec
1421 TYPE(gto_basis_set_type), POINTER :: oce_basis
1422 TYPE(mp_comm_type) :: group
1423 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1424 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1425 TYPE(qs_ks_env_type), POINTER :: ks_env
1426
1427 CALL timeset(routinen, handle)
1428
1429 ! basis sets: block sizes
1430 CALL dbcsr_get_info(sxo(1)%matrix, row_blk_size=oce_blk_sizes, col_blk_size=orb_blk_sizes, &
1431 distribution=dbcsr_dist)
1432 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, natom=natom)
1433 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, qs_kind_set=qs_kind_set)
1434 ALLOCATE (iao_blk_sizes(natom))
1435 DO iatom = 1, natom
1436 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1437 CALL get_qs_kind(qs_kind_set(ikind), nsgf=ns, basis_type="MIN")
1438 iao_blk_sizes(iatom) = ns
1439 END DO
1440
1441 CALL dbcsr_create(iao_vec, name="IAO_VEC", dist=dbcsr_dist, &
1442 matrix_type=dbcsr_type_no_symmetry, row_blk_size=orb_blk_sizes, &
1443 col_blk_size=iao_blk_sizes)
1444 CALL dbcsr_create(sx_vec, name="SX_VEC", dist=dbcsr_dist, &
1445 matrix_type=dbcsr_type_no_symmetry, row_blk_size=oce_blk_sizes, &
1446 col_blk_size=iao_blk_sizes)
1447 CALL dbcsr_reserve_diag_blocks(matrix=sx_vec)
1448
1449 nkind = SIZE(smat_kind)
1450 ALLOCATE (sinv_kind(nkind))
1451 DO ikind = 1, nkind
1452 noce = SIZE(smat_kind(ikind)%array, 1)
1453 ALLOCATE (sinv_kind(ikind)%array(noce, noce))
1454 sinv_kind(ikind)%array = smat_kind(ikind)%array
1455 CALL invmat_symm(sinv_kind(ikind)%array)
1456 END DO
1457 CALL dbcsr_get_info(iao_vec, group=group)
1458
1459 nspin = SIZE(iao_coef, 1)
1460 ALLOCATE (oce_comp(natom, nspin))
1461 oce_comp = 0.0_dp
1462 DO ispin = 1, nspin
1463 CALL copy_fm_to_dbcsr(iao_coef(ispin), iao_vec)
1464 CALL dbcsr_multiply("N", "N", 1.0_dp, sxo(1)%matrix, iao_vec, 0.0_dp, sx_vec, &
1465 retain_sparsity=.true.)
1466 DO iatom = 1, natom
1467 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1468 CALL dbcsr_get_block_p(matrix=sx_vec, &
1469 row=iatom, col=iatom, block=ablock, found=found)
1470 IF (found) THEN
1471 n = SIZE(ablock, 1)
1472 m = SIZE(ablock, 2)
1473 ablock = matmul(sinv_kind(ikind)%array, ablock)
1474 ALLOCATE (amat(n, m), bmat(m, m), ev(m), vec(m, m))
1475 amat(1:n, 1:m) = matmul(smat_kind(ikind)%array(1:n, 1:n), ablock(1:n, 1:m))
1476 bmat(1:m, 1:m) = matmul(transpose(ablock(1:n, 1:m)), amat(1:n, 1:m))
1477 CALL jacobi(bmat, ev, vec)
1478 oce_comp(iatom, ispin) = sum(ev(1:m))/real(m, kind=dp)
1479 DO i = 1, m
1480 ev(i) = 1._dp/sqrt(ev(i))
1481 bmat(1:m, i) = vec(1:m, i)*ev(i)
1482 END DO
1483 bmat(1:m, 1:m) = matmul(bmat(1:m, 1:m), transpose(vec(1:m, 1:m)))
1484 ablock(1:n, 1:m) = matmul(ablock(1:n, 1:m), bmat(1:m, 1:m))
1485 DEALLOCATE (amat, bmat, ev, vec)
1486 END IF
1487 END DO
1488 CALL dbcsr_replicate_all(sx_vec)
1489 DO iatom = 1, natom
1490 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1491 CALL dbcsr_get_block_p(matrix=sx_vec, &
1492 row=iatom, col=iatom, block=ablock, found=found)
1493 cpassert(found)
1494 n = SIZE(ablock, 1)
1495 m = SIZE(ablock, 2)
1496 ALLOCATE (oce_atom(iatom, ispin)%array(n, m))
1497 oce_atom(iatom, ispin)%array = ablock
1498 END DO
1499 END DO
1500
1501 CALL group%sum(oce_comp)
1502 IF (unit_nr > 0) THEN
1503 DO iatom = 1, natom
1504 WRITE (unit_nr, "(T4,A,I6,T30,A,2F12.4)") "Atom:", iatom, " Completeness: ", oce_comp(iatom, 1:nspin)
1505 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1506 oce_basis => oce_basis_set_list(ikind)%gto_basis_set
1507 CALL get_gto_basis_set(oce_basis, nset=nset, nshell=nshell, nsgf=nsgf, maxl=maxl, &
1508 l=lval, first_sgf=first_sgf, last_sgf=last_sgf)
1509 ALLOCATE (filter(nsgf, 0:maxl), vector(nsgf), prol(0:maxl, nspin))
1510 filter = 0.0_dp
1511 DO iset = 1, nset
1512 DO ishell = 1, nshell(iset)
1513 l = lval(ishell, iset)
1514 i1 = first_sgf(ishell, iset)
1515 i2 = last_sgf(ishell, iset)
1516 filter(i1:i2, l) = 1.0_dp
1517 END DO
1518 END DO
1519 !
1520 n = SIZE(oce_atom(iatom, 1)%array, 1)
1521 m = SIZE(oce_atom(iatom, 1)%array, 2)
1522 cpassert(n == nsgf)
1523 DO iao = 1, m
1524 prol = 0.0_dp
1525 DO ispin = 1, nspin
1526 DO l = 0, maxl
1527 vector(1:n) = oce_atom(iatom, ispin)%array(1:n, iao)*filter(1:n, l)
1528 prol(l, ispin) = sum(matmul(smat_kind(ikind)%array(1:n, 1:n), vector(1:n))*vector(1:n))
1529 END DO
1530 END DO
1531 WRITE (unit_nr, "(T4,I3,T15,A,T39,(6F7.4))") iao, " l-contributions:", (sum(prol(l, :)), l=0, maxl)
1532 END DO
1533 DEALLOCATE (filter, vector, prol)
1534 END DO
1535 WRITE (unit_nr, *)
1536 END IF
1537
1538 DEALLOCATE (oce_comp)
1539 DO ikind = 1, nkind
1540 DEALLOCATE (sinv_kind(ikind)%array)
1541 END DO
1542 DEALLOCATE (sinv_kind)
1543 DEALLOCATE (iao_blk_sizes)
1544 CALL dbcsr_release(iao_vec)
1545 CALL dbcsr_release(sx_vec)
1546
1547 CALL timestop(handle)
1548
1549 END SUBROUTINE iao_oce_expansion
1550
1551! **************************************************************************************************
1552!> \brief ...
1553!> \param smat_kind ...
1554!> \param basis_set_list ...
1555! **************************************************************************************************
1556 SUBROUTINE oce_overlap_matrix(smat_kind, basis_set_list)
1557 TYPE(cp_2d_r_p_type), DIMENSION(:) :: smat_kind
1558 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1559
1560 CHARACTER(len=*), PARAMETER :: routinen = 'oce_overlap_matrix'
1561
1562 INTEGER :: handle, ikind, iset, jset, ldsab, m1, &
1563 m2, n1, n2, ncoa, ncob, nkind, nseta, &
1564 sgfa, sgfb
1565 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
1566 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
1567 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: oint, owork
1568 REAL(kind=dp), DIMENSION(3) :: rab
1569 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, scon_a, smat, zeta
1570 TYPE(gto_basis_set_type), POINTER :: basis_set
1571
1572 CALL timeset(routinen, handle)
1573
1574 rab(1:3) = 0.0_dp
1575
1576 nkind = SIZE(smat_kind)
1577 ldsab = 0
1578 DO ikind = 1, nkind
1579 basis_set => basis_set_list(ikind)%gto_basis_set
1580 cpassert(ASSOCIATED(basis_set))
1581 CALL get_gto_basis_set(gto_basis_set=basis_set, maxco=m1, nsgf=m2)
1582 ldsab = max(m1, m2, ldsab)
1583 END DO
1584
1585 ALLOCATE (oint(ldsab, ldsab), owork(ldsab, ldsab))
1586
1587 DO ikind = 1, nkind
1588 basis_set => basis_set_list(ikind)%gto_basis_set
1589 cpassert(ASSOCIATED(basis_set))
1590 smat => smat_kind(ikind)%array
1591 smat = 0.0_dp
1592 ! basis ikind
1593 first_sgfa => basis_set%first_sgf
1594 la_max => basis_set%lmax
1595 la_min => basis_set%lmin
1596 npgfa => basis_set%npgf
1597 nseta = basis_set%nset
1598 nsgfa => basis_set%nsgf_set
1599 rpgfa => basis_set%pgf_radius
1600 scon_a => basis_set%scon
1601 zeta => basis_set%zet
1602
1603 DO iset = 1, nseta
1604
1605 ncoa = npgfa(iset)*ncoset(la_max(iset))
1606 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1607 sgfa = first_sgfa(1, iset)
1608
1609 DO jset = 1, nseta
1610
1611 ncob = npgfa(jset)*ncoset(la_max(jset))
1612 n2 = npgfa(jset)*(ncoset(la_max(jset)) - ncoset(la_min(jset) - 1))
1613 sgfb = first_sgfa(1, jset)
1614
1615 ! calculate integrals and derivatives
1616 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1617 la_max(jset), la_min(jset), npgfa(jset), rpgfa(:, jset), zeta(:, jset), &
1618 rab, sab=oint(:, :))
1619 ! Contraction
1620 CALL contraction(oint(:, :), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1621 cb=scon_a(:, sgfb:), nb=n2, mb=nsgfa(jset), fscale=1.0_dp, trans=.false.)
1622 CALL block_add("IN", owork, nsgfa(iset), nsgfa(jset), smat, &
1623 sgfa, sgfb, trans=.false.)
1624
1625 END DO
1626 END DO
1627
1628 END DO
1629 DEALLOCATE (oint, owork)
1630
1631 CALL timestop(handle)
1632
1633 END SUBROUTINE oce_overlap_matrix
1634
1635! **************************************************************************************************
1636!> \brief 2 by 2 rotation optimization for the Pipek-Mezey operator
1637!> \param zij_fm_set ...
1638!> \param vectors ...
1639!> \param order ...
1640!> \param accuracy ...
1641!> \param unit_nr ...
1642!> \par History
1643!> 05-2005 created
1644!> 10-2023 adapted from jacobi_rotation_pipek [JGH]
1645!> \author MI
1646! **************************************************************************************************
1647 SUBROUTINE pm_localization(zij_fm_set, vectors, order, accuracy, unit_nr)
1648
1649 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: zij_fm_set
1650 TYPE(cp_fm_type), INTENT(IN) :: vectors
1651 INTEGER, INTENT(IN) :: order
1652 REAL(kind=dp), INTENT(IN) :: accuracy
1653 INTEGER, INTENT(IN) :: unit_nr
1654
1655 INTEGER, PARAMETER :: max_sweeps = 250
1656
1657 INTEGER :: iatom, istate, jstate, natom, nstate, &
1658 sweeps
1659 REAL(kind=dp) :: aij, bij, ct, ftarget, mii, mij, mjj, &
1660 st, t1, t2, theta, tolerance, tt
1661 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fdiag
1662 TYPE(cp_fm_type) :: rmat
1663
1664 CALL cp_fm_create(rmat, zij_fm_set(1, 1)%matrix_struct)
1665 CALL cp_fm_set_all(rmat, 0.0_dp, 1.0_dp)
1666
1667 CALL cp_fm_get_info(rmat, nrow_global=nstate)
1668 ALLOCATE (fdiag(nstate))
1669 tolerance = 1.0e10_dp
1670 sweeps = 0
1671 natom = SIZE(zij_fm_set, 1)
1672 IF (unit_nr > 0) THEN
1673 WRITE (unit_nr, '(A,T30,A,T45,A,T63,A,T77,A)') &
1674 " Jacobi Localization ", "Sweep", "Functional", "Gradient", "Time"
1675 END IF
1676 ! do jacobi sweeps until converged
1677 DO sweeps = 1, max_sweeps
1678 t1 = m_walltime()
1679 tolerance = 0.0_dp
1680 DO istate = 1, nstate
1681 DO jstate = istate + 1, nstate
1682 aij = 0.0_dp
1683 bij = 0.0_dp
1684 DO iatom = 1, natom
1685 CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, istate, mii)
1686 CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, jstate, mij)
1687 CALL cp_fm_get_element(zij_fm_set(iatom, 1), jstate, jstate, mjj)
1688 IF (order == 2) THEN
1689 aij = aij + 4._dp*mij**2 - (mii - mjj)**2
1690 bij = bij + 4._dp*mij*(mii - mjj)
1691 ELSEIF (order == 4) THEN
1692 aij = aij - mii**4 - mjj**4 + 6._dp*(mii**2 + mjj**2)*mij**2 + &
1693 mii**3*mjj + mii*mjj**3
1694 bij = bij + 4._dp*mij*(mii**3 - mjj**3)
1695 ELSE
1696 cpabort("PM order")
1697 END IF
1698 END DO
1699 IF ((aij**2 + bij**2) < 1.e-10_dp) cycle
1700 tolerance = tolerance + bij**2
1701 theta = 0.25_dp*atan2(bij, -aij)
1702 ct = cos(theta)
1703 st = sin(theta)
1704
1705 CALL cp_fm_rot_cols(rmat, istate, jstate, ct, st)
1706
1707 DO iatom = 1, natom
1708 CALL cp_fm_rot_cols(zij_fm_set(iatom, 1), istate, jstate, ct, st)
1709 CALL cp_fm_rot_rows(zij_fm_set(iatom, 1), istate, jstate, ct, st)
1710 END DO
1711 END DO
1712 END DO
1713 ftarget = 0.0_dp
1714 DO iatom = 1, natom
1715 CALL cp_fm_get_diag(zij_fm_set(iatom, 1), fdiag)
1716 ftarget = ftarget + sum(fdiag**order)
1717 END DO
1718 ftarget = ftarget**(1._dp/order)
1719 tolerance = sqrt(tolerance)
1720 t2 = m_walltime()
1721 tt = t2 - t1
1722 IF (unit_nr > 0) THEN
1723 WRITE (unit_nr, '(T31,I4,T39,F16.8,T55,G16.8,T71,F10.2)') sweeps, ftarget, tolerance, tt
1724 END IF
1725 IF (tolerance < accuracy) EXIT
1726 END DO
1727 DEALLOCATE (fdiag)
1728
1729 CALL rotate_orbitals(rmat, vectors)
1730 CALL cp_fm_release(rmat)
1731
1732 END SUBROUTINE pm_localization
1733! **************************************************************************************************
1734
1735END 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)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:681
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:414
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)
...
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...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
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_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)
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_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
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
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
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:153
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:1476
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)
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, 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, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, 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_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 Berry operator for periodic systems used to define the spread of the MOS Here the matrix...
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 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)
Get the components of a MO set data structure.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
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)
...
Definition qs_moments.F:583
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)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:574
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 ...