59 USE dbcsr_api,
ONLY: &
60 dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_get_info, &
61 dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
62 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
63 dbcsr_replicate_all, dbcsr_reserve_diag_blocks, dbcsr_set, dbcsr_trace, dbcsr_type, &
64 dbcsr_type_no_symmetry
118#include "./base/base_uses.f90"
123 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'iao_analysis'
128 MODULE PROCEDURE iao_calculate_dmat_diag, &
129 iao_calculate_dmat_full
148 INTEGER,
INTENT(IN) :: unit_nr
149 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
150 OPTIONAL :: c_iao_coef
151 TYPE(
mo_set_type),
DIMENSION(:),
OPTIONAL,
POINTER :: mos
152 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL :: bond_centers
154 CHARACTER(len=*),
PARAMETER :: routinen =
'iao_wfn_analysis'
156 CHARACTER(LEN=2) :: element_symbol
157 CHARACTER(LEN=default_string_length) :: bname
158 INTEGER :: handle, i, iatom, ikind, isgf, ispin, &
159 nao, natom, nimages, nkind, no, norb, &
160 nref, ns, nsgf, nspin, nvec, nx, order
161 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: first_sgf, last_sgf, nsgfat
162 INTEGER,
DIMENSION(2) :: nocc
163 LOGICAL :: cubes_iao, cubes_ibo, molden_iao, &
164 molden_ibo, uniform_occupation
165 REAL(kind=
dp) :: fin, fout, t1, t2, trace, zval
166 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: fdiag
167 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: mcharge
168 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: moments
169 REAL(kind=
dp),
DIMENSION(:),
POINTER :: occupation_numbers
174 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: c_loc_orb, ciao, cloc, cvec, iao_coef
175 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:, :) :: zij_atom
177 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: sref, sro, sxo
178 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_s
179 TYPE(dbcsr_type) :: dmat
180 TYPE(dbcsr_type),
POINTER :: smat
185 TYPE(
mo_set_type),
ALLOCATABLE,
DIMENSION(:) :: mo_iao, mo_loc
188 POINTER :: sro_list, srr_list, sxo_list
190 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
194 ibo_cc_section, ibo_cubes_section, &
198 cpassert(iao_env%do_iao)
201 CALL get_qs_env(qs_env, dft_control=dft_control)
202 nspin = dft_control%nspins
203 nimages = dft_control%nimages
204 IF (nimages > 1)
THEN
205 IF (unit_nr > 0)
THEN
206 WRITE (unit=unit_nr, fmt=
"(T2,A)") &
207 "K-Points: Intrinsic Atomic Orbitals Analysis not available."
210 IF (nimages > 1)
RETURN
212 CALL timeset(routinen, handle)
214 IF (unit_nr > 0)
THEN
215 WRITE (unit_nr,
'(/,T2,A)')
'!-----------------------------------------------------------------------------!'
216 WRITE (unit=unit_nr, fmt=
"(T24,A)")
"INTRINSIC ATOMIC ORBITALS ANALYSIS"
217 WRITE (unit=unit_nr, fmt=
"(T13,A)")
"G. Knizia, J. Chem. Theory Comput. 9, 4834-4843 (2013)"
218 WRITE (unit_nr,
'(T2,A)')
'!-----------------------------------------------------------------------------!'
223 iao_molden_section => iao_env%iao_molden_section
224 iao_cubes_section => iao_env%iao_cubes_section
225 ibo_molden_section => iao_env%ibo_molden_section
226 ibo_cubes_section => iao_env%ibo_cubes_section
227 ibo_cc_section => iao_env%ibo_cc_section
231 IF (
ASSOCIATED(iao_molden_section))
THEN
235 IF (
ASSOCIATED(iao_cubes_section))
THEN
239 IF (
ASSOCIATED(ibo_molden_section))
THEN
243 IF (
ASSOCIATED(ibo_cubes_section))
THEN
252 smat => matrix_s(1, 1)%matrix
254 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
255 nkind =
SIZE(qs_kind_set)
256 ALLOCATE (ref_basis_set_list(nkind), orb_basis_set_list(nkind))
258 qs_kind => qs_kind_set(ikind)
259 NULLIFY (ref_basis_set_list(ikind)%gto_basis_set)
260 NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
261 NULLIFY (refbasis, orbbasis)
262 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orbbasis, basis_type=
"ORB")
263 IF (
ASSOCIATED(orbbasis)) orb_basis_set_list(ikind)%gto_basis_set => orbbasis
264 CALL get_qs_kind(qs_kind=qs_kind, basis_set=refbasis, basis_type=
"MIN")
265 IF (
ASSOCIATED(refbasis)) ref_basis_set_list(ikind)%gto_basis_set => refbasis
269 NULLIFY (srr_list, sro_list)
277 ref_basis_set_list, ref_basis_set_list, srr_list)
279 ref_basis_set_list, orb_basis_set_list, sro_list)
281 IF (
PRESENT(mos))
THEN
287 CALL dbcsr_get_info(sro(1)%matrix, nfullrows_total=nref, nfullcols_total=nao)
291 CALL cp_fm_struct_create(fm_struct, context=mo_coeff%matrix_struct%context, nrow_global=nref, &
292 ncol_global=nao, para_env=mo_coeff%matrix_struct%para_env)
296 CALL cp_fm_struct_create(fm_struct, context=mo_coeff%matrix_struct%context, nrow_global=nao, &
297 ncol_global=nref, para_env=mo_coeff%matrix_struct%para_env)
301 CALL iao_projectors(smat, sref(1)%matrix, sro(1)%matrix, p_orb_ref, p_ref_orb, iao_env%eps_svd)
306 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nmo=norb, uniform_occupation=uniform_occupation, &
307 occupation_numbers=occupation_numbers)
308 IF (uniform_occupation)
THEN
313 IF (occupation_numbers(i) > iao_env%eps_occ)
THEN
323 ALLOCATE (iao_coef(nspin), cvec(nspin))
325 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff)
329 CALL cp_fm_struct_create(fm_struct, ncol_global=nvec, template_fmstruct=mo_coeff%matrix_struct)
335 CALL cp_fm_struct_create(fm_struct, ncol_global=nref, template_fmstruct=mo_coeff%matrix_struct)
339 CALL intrinsic_ao_calc(smat, p_orb_ref, p_ref_orb, cvec(ispin), iao_coef(ispin))
343 IF (iao_env%do_charges)
THEN
345 ALLOCATE (ciao(nspin))
347 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=norb)
350 CALL cp_fm_struct_create(fm_struct, nrow_global=nref, template_fmstruct=mo_coeff%matrix_struct)
355 CALL parallel_gemm(
'T',
'N', nref, norb, nao, 1.0_dp, iao_coef(ispin), smo, 0.0_dp, ciao(ispin))
360 IF (unit_nr > 0)
THEN
361 WRITE (unit_nr,
'(/,T2,A)')
'Intrinsic AO Population Analysis '
363 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
364 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
365 ALLOCATE (mcharge(natom, nspin))
366 CALL dbcsr_get_info(sref(1)%matrix, nfullrows_total=nref)
369 uniform_occupation=uniform_occupation, &
370 occupation_numbers=occupation_numbers)
372 CALL dbcsr_create(dmat, template=sref(1)%matrix)
373 CALL dbcsr_reserve_diag_blocks(dmat)
376 CALL dbcsr_trace(dmat, trace)
377 IF (unit_nr > 0)
THEN
378 WRITE (unit_nr,
'(T2,A,I2,T66,F15.4)')
'Number of Electrons: Trace(Piao) Spin ', ispin, trace
381 CALL dbcsr_release(dmat)
383 CALL print_atomic_charges(particle_set, qs_kind_set, unit_nr,
"Intrinsic Atomic Orbital Charges", &
384 electronic_charges=mcharge)
397 DEALLOCATE (orb_basis_set_list)
399 IF (iao_env%do_oce)
THEN
401 IF (unit_nr > 0)
THEN
402 WRITE (unit_nr,
'(T2,A)')
"IAO One-Center Expansion: OCE Basis Set"
404 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
405 nkind =
SIZE(qs_kind_set)
406 ALLOCATE (oce_basis_set_list(nkind), orb_basis_set_list(nkind))
408 qs_kind => qs_kind_set(ikind)
409 NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
411 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orbbasis, basis_type=
"ORB")
412 IF (
ASSOCIATED(orbbasis)) orb_basis_set_list(ikind)%gto_basis_set => orbbasis
415 orbbasis => orb_basis_set_list(ikind)%gto_basis_set
416 cpassert(
ASSOCIATED(orbbasis))
418 CALL create_oce_basis(ocebasis, orbbasis, iao_env%lmax_oce, iao_env%nbas_oce)
421 oce_basis_set_list(ikind)%gto_basis_set => ocebasis
422 IF (unit_nr > 0)
THEN
423 qs_kind => qs_kind_set(ikind)
424 CALL get_qs_kind(qs_kind, zeff=zval, element_symbol=element_symbol)
426 WRITE (unit_nr,
'(T2,A,A,T14,A,I4,T40,A,A30)')
"Kind: ", element_symbol,
"NBasFun: ", nsgf, &
427 "OCE Basis: ", adjustl(trim(bname))
430 IF (unit_nr > 0)
WRITE (unit_nr, *)
432 ALLOCATE (smat_kind(nkind))
434 ocebasis => oce_basis_set_list(ikind)%gto_basis_set
436 ALLOCATE (smat_kind(ikind)%array(nsgf, nsgf))
438 CALL oce_overlap_matrix(smat_kind, oce_basis_set_list)
440 NULLIFY (sxo, sxo_list)
444 oce_basis_set_list, orb_basis_set_list, sxo_list)
448 ALLOCATE (oce_atom(natom, nspin))
449 CALL iao_oce_expansion(qs_env, oce_basis_set_list, smat_kind, sxo, iao_coef, oce_atom, unit_nr)
452 ocebasis => oce_basis_set_list(ikind)%gto_basis_set
455 DEALLOCATE (oce_basis_set_list, orb_basis_set_list)
460 DEALLOCATE (smat_kind(ikind)%array)
462 DEALLOCATE (smat_kind)
465 DEALLOCATE (oce_atom(iatom, ispin)%array)
468 DEALLOCATE (oce_atom)
473 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
474 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
475 IF (unit_nr > 0)
THEN
476 WRITE (unit_nr,
'(T2,A)')
' Write IAO in MOLDEN Format'
478 ALLOCATE (mo_iao(nspin))
481 CALL allocate_mo_set(mo_iao(ispin), nao, nref, nref, 0.0_dp, 1.0_dp, 0.0_dp)
482 CALL init_mo_set(mo_iao(ispin), fm_ref=iao_coef(ispin), name=
"iao_set")
483 CALL cp_fm_to_fm(iao_coef(ispin), mo_iao(ispin)%mo_coeff)
485 mo_iao(ispin)%occupation_numbers = 1.0_dp
488 CALL write_mos_molden(mo_iao, qs_kind_set, particle_set, iao_molden_section)
495 IF (unit_nr > 0)
THEN
496 WRITE (unit_nr,
'(T2,A)')
' Write IAO as CUBE Files'
499 CALL print_iao_cubes(qs_env, iao_cubes_section, iao_coef, ref_basis_set_list)
503 IF (iao_env%do_bondorbitals)
THEN
504 IF (unit_nr > 0)
THEN
505 WRITE (unit_nr,
'(/,T2,A)')
'Intrinsic Bond Orbital Generation'
508 ALLOCATE (ciao(nspin))
510 CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=norb)
513 CALL cp_fm_struct_create(fm_struct, nrow_global=nref, template_fmstruct=mo_coeff%matrix_struct)
518 CALL parallel_gemm(
'T',
'N', nref, norb, nao, 1.0_dp, iao_coef(ispin), smo, 0.0_dp, ciao(ispin))
524 ALLOCATE (cloc(nspin), c_loc_orb(nspin))
528 template_fmstruct=ciao(ispin)%matrix_struct)
531 CALL cp_fm_to_fm(ciao(ispin), cloc(ispin), ncol=nocc(ispin))
534 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
535 ALLOCATE (first_sgf(natom), last_sgf(natom), nsgfat(natom))
536 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
537 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf, &
538 nsgf=nsgfat, basis=ref_basis_set_list)
541 cpabort(
"IAO localization operator NYA")
544 cpabort(
"IAO energy weight localization NYA")
550 ALLOCATE (zij_atom(natom, 1), fdiag(nvec))
553 template_fmstruct=cloc(ispin)%matrix_struct)
556 isgf = first_sgf(iatom)
557 CALL parallel_gemm(
'T',
'N', nvec, nvec, nsgfat(iatom), 1.0_dp, cloc(ispin), cloc(ispin), &
558 0.0_dp, zij_atom(iatom, 1), &
559 a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf)
568 fin = fin + sum(fdiag**order)
570 fin = fin**(1._dp/order)
575 isgf = first_sgf(iatom)
576 CALL parallel_gemm(
'T',
'N', nvec, nvec, nsgfat(iatom), 1.0_dp, cloc(ispin), cloc(ispin), &
577 0.0_dp, zij_atom(iatom, 1), &
578 a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf)
583 fout = fout + sum(fdiag**order)
585 fout = fout**(1._dp/order)
588 IF (unit_nr > 0)
THEN
589 WRITE (unit_nr,
'(T2,A,F14.8,A,F14.8,A,F8.2)') &
590 'SCDM pre-localization: fin=', fin,
" fout=", fout,
" Time=", t2 - t1
593 CALL pm_localization(zij_atom, cloc(ispin), order, 1.e-12_dp, unit_nr)
599 template_fmstruct=cloc(ispin)%matrix_struct)
602 CALL parallel_gemm(
'N',
'N', nao, nvec, nref, 1.0_dp, iao_coef(ispin), cloc(ispin), &
603 0.0_dp, c_loc_orb(ispin))
606 DEALLOCATE (first_sgf, last_sgf, nsgfat)
609 IF (iao_env%do_center)
THEN
610 IF (unit_nr > 0)
THEN
611 WRITE (unit_nr,
'(T2,A)')
' Calculate Localized Orbital Centers and Spread'
612 IF (iao_env%pos_periodic)
THEN
613 WRITE (unit_nr,
'(T2,A)')
' Use Berry Phase Position Operator'
615 WRITE (unit_nr,
'(T2,A)')
' Use Local Position Operator'
620 ALLOCATE (moments(5, nvec, nspin))
622 IF (iao_env%pos_periodic)
THEN
623 CALL center_spread_berry(qs_env, c_loc_orb, moments)
625 CALL center_spread_loc(qs_env, c_loc_orb, moments)
628 IF (
ASSOCIATED(ibo_cc_section))
THEN
632 IF (
PRESENT(bond_centers))
THEN
633 nx =
SIZE(bond_centers, 1)
634 no =
SIZE(bond_centers, 2)
635 ns =
SIZE(bond_centers, 3)
637 cpassert(ns == nspin)
640 bond_centers(1:nx, 1:nvec, 1:ns) = moments(1:nx, 1:nvec, 1:ns)
646 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
647 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
648 IF (unit_nr > 0)
THEN
649 WRITE (unit_nr,
'(T2,A)')
' Write IBO in MOLDEN Format'
651 ALLOCATE (mo_loc(nspin))
655 CALL allocate_mo_set(mo_loc(ispin), nao, nvec, nvec, 0.0_dp, 1.0_dp, 0.0_dp)
656 CALL init_mo_set(mo_loc(ispin), fm_ref=c_loc_orb(ispin), name=
"ibo_orb")
657 CALL cp_fm_to_fm(c_loc_orb(ispin), mo_loc(ispin)%mo_coeff)
659 mo_loc(ispin)%occupation_numbers = 1.0_dp
662 CALL write_mos_molden(mo_loc, qs_kind_set, particle_set, ibo_molden_section)
670 IF (unit_nr > 0)
THEN
671 WRITE (unit_nr,
'(T2,A)')
' Write IBO on CUBE files'
674 CALL print_ibo_cubes(qs_env, ibo_cubes_section, c_loc_orb)
677 IF (
PRESENT(mos) .AND.
ALLOCATED(c_loc_orb))
THEN
680 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
682 CALL cp_fm_to_fm(c_loc_orb(ispin), mo_coeff, ncol=nvec)
687 DEALLOCATE (ref_basis_set_list)
689 IF (
PRESENT(c_iao_coef))
THEN
691 ALLOCATE (c_iao_coef(nspin))
693 CALL cp_fm_create(c_iao_coef(ispin), iao_coef(ispin)%matrix_struct)
694 CALL cp_fm_to_fm(iao_coef(ispin), c_iao_coef(ispin))
699 IF (unit_nr > 0)
THEN
700 WRITE (unit_nr,
'(T2,A)') &
701 '!----------------------------END OF IAO ANALYSIS------------------------------!'
704 CALL timestop(handle)
717 SUBROUTINE iao_projectors(smat, sref, s_r_o, p_o_r, p_r_o, eps_svd)
718 TYPE(dbcsr_type),
INTENT(INOUT) :: smat, sref, s_r_o
719 TYPE(
cp_fm_type),
INTENT(INOUT) :: p_o_r, p_r_o
720 REAL(kind=
dp),
INTENT(IN) :: eps_svd
722 CHARACTER(len=*),
PARAMETER :: routinen =
'iao_projectors'
724 INTEGER :: handle, norb, nref
728 CALL timeset(routinen, handle)
736 template_fmstruct=p_r_o%matrix_struct)
742 CALL parallel_gemm(
'N',
'T', norb, nref, norb, 1.0_dp, fm_inv, fm_sro, 0.0_dp, p_o_r)
748 template_fmstruct=p_r_o%matrix_struct)
754 CALL parallel_gemm(
'N',
'N', nref, norb, nref, 1.0_dp, fm_inv, fm_sro, 0.0_dp, p_r_o)
760 CALL timestop(handle)
762 END SUBROUTINE iao_projectors
772 SUBROUTINE intrinsic_ao_calc(smat, p_o_r, p_r_o, cvec, avec)
773 TYPE(dbcsr_type),
INTENT(INOUT) :: smat
774 TYPE(
cp_fm_type),
INTENT(INOUT) :: p_o_r, p_r_o, cvec, avec
776 CHARACTER(len=*),
PARAMETER :: routinen =
'intrinsic_ao_calc'
778 INTEGER :: handle, nao, norb, nref
782 CALL timeset(routinen, handle)
791 template_fmstruct=cvec%matrix_struct)
795 CALL parallel_gemm(
'N',
'N', nref, norb, nao, 1.0_dp, p_r_o, cvec, 0.0_dp, pc)
797 CALL parallel_gemm(
'N',
'N', nao, norb, nref, 1.0_dp, p_o_r, pc, 0.0_dp, ctvec)
806 template_fmstruct=cvec%matrix_struct)
810 CALL parallel_gemm(
'T',
'N', norb, nref, nao, 1.0_dp, sct, p_o_r, 0.0_dp, pc)
813 template_fmstruct=cvec%matrix_struct)
816 CALL parallel_gemm(
'N',
'N', nao, nref, norb, 1.0_dp, ctvec, pc, 0.0_dp, vec1)
818 CALL parallel_gemm(
'T',
'N', norb, nref, nao, 1.0_dp, sc, vec1, 0.0_dp, pc)
819 CALL parallel_gemm(
'N',
'N', nao, nref, norb, 1.0_dp, cvec, pc, 0.0_dp, avec)
821 CALL cp_fm_geadd(1.0_dp,
'N', p_o_r, -1.0_dp, vec1)
825 CALL parallel_gemm(
'T',
'N', norb, nref, nao, 1.0_dp, sc, vec1, 0.0_dp, pc)
826 CALL parallel_gemm(
'N',
'N', nao, nref, norb, -1.0_dp, cvec, pc, 1.0_dp, avec)
837 CALL timestop(handle)
839 END SUBROUTINE intrinsic_ao_calc
848 SUBROUTINE iao_calculate_dmat_diag(cvec, density_matrix, occupation, uniform_occupation)
851 TYPE(dbcsr_type) :: density_matrix
852 REAL(KIND=
dp),
DIMENSION(:),
INTENT(IN) :: occupation
853 LOGICAL,
INTENT(IN) :: uniform_occupation
855 CHARACTER(len=*),
PARAMETER :: routineN =
'iao_calculate_dmat_diag'
857 INTEGER :: handle, ncol
858 REAL(KIND=
dp) :: alpha
861 CALL timeset(routinen, handle)
863 CALL dbcsr_set(density_matrix, 0.0_dp)
866 IF (.NOT. uniform_occupation)
THEN
872 matrix_v=cvec, matrix_g=fm_tmp, &
873 ncol=ncol, alpha=alpha)
876 alpha = occupation(1)
878 matrix_v=cvec, ncol=ncol, alpha=alpha)
881 CALL timestop(handle)
883 END SUBROUTINE iao_calculate_dmat_diag
892 SUBROUTINE iao_calculate_dmat_full(cvec, density_matrix, weight, occmat)
895 TYPE(dbcsr_type) :: density_matrix
896 REAL(KIND=
dp),
DIMENSION(:),
INTENT(IN) :: weight
899 CHARACTER(len=*),
PARAMETER :: routinen =
'iao_calculate_dmat_full'
901 INTEGER :: handle, ic, jc, ncol
902 REAL(kind=
dp) :: alpha
906 CALL timeset(routinen, handle)
908 CALL dbcsr_set(density_matrix, 0.0_dp)
911 template_fmstruct=cvec%matrix_struct)
918 CALL cp_fm_to_fm(cvec, fm1, ncol=1, source_start=ic, target_start=1)
920 CALL cp_fm_to_fm(cvec, fm2, ncol=1, source_start=jc, target_start=1)
922 alpha = weight(ic)*alpha
924 alpha=alpha, symmetry_mode=1)
930 CALL timestop(handle)
932 END SUBROUTINE iao_calculate_dmat_full
940 TYPE(dbcsr_type) :: p_matrix
941 REAL(kind=
dp),
DIMENSION(:) :: charges
943 INTEGER :: blk, group_handle, i, iblock_col, &
945 REAL(kind=
dp) :: trace
946 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: p_block
947 TYPE(dbcsr_iterator_type) :: iter
952 CALL dbcsr_iterator_start(iter, p_matrix)
953 DO WHILE (dbcsr_iterator_blocks_left(iter))
955 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, p_block, blk)
956 cpassert(iblock_row == iblock_col)
958 DO i = 1,
SIZE(p_block, 1)
959 trace = trace + p_block(i, i)
961 charges(iblock_row) = trace
963 CALL dbcsr_iterator_stop(iter)
965 CALL dbcsr_get_info(p_matrix, group=group_handle)
966 CALL group%set_handle(group_handle)
968 CALL group%sum(charges)
979 SUBROUTINE print_iao_cubes(qs_env, print_section, iao_coef, basis_set_list)
982 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: iao_coef
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
1000 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
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)
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)
1016 first_bas(i) = first_bas(i - 1) + blk_sizes(i)
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)
1023 nspins =
SIZE(iao_coef)
1024 nstart = min(1, n_rep)
1026 DO ispin = 1, nspins
1027 DO i_rep = nstart, n_rep
1029 IF (i_rep == 0)
THEN
1033 nat =
SIZE(atom_list)
1036 IF (i_rep == 0)
THEN
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
1047 middle_name=trim(filename), file_position=
"REWIND", log_filename=.false., &
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)
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)
1061 END SUBROUTINE print_iao_cubes
1069 SUBROUTINE print_ibo_cubes(qs_env, print_section, ibo_coef)
1072 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: ibo_coef
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
1089 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1092 cpassert(
ASSOCIATED(print_section))
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)
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)
1106 nspins =
SIZE(ibo_coef)
1107 nstart = min(1, n_rep)
1109 DO ispin = 1, nspins
1110 DO i_rep = nstart, n_rep
1112 IF (i_rep == 0)
THEN
1116 nstate =
SIZE(state_list)
1119 IF (i_rep == 0)
THEN
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
1129 middle_name=trim(filename), file_position=
"REWIND", log_filename=.false., &
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)
1139 CALL auxbas_pw_pool%give_back_pw(wf_r)
1140 CALL auxbas_pw_pool%give_back_pw(wf_g)
1142 END SUBROUTINE print_ibo_cubes
1152 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(IN) :: moments
1153 INTEGER,
DIMENSION(:),
INTENT(IN) :: nocc
1155 INTEGER,
INTENT(IN),
OPTIONAL :: iounit
1157 CHARACTER(LEN=default_path_length) :: filename
1158 INTEGER :: is, ispin, iw, nspin
1162 nspin =
SIZE(moments, 3)
1164 IF (
PRESENT(print_section))
THEN
1165 WRITE (filename,
'(A18,I1.1)')
"IBO_CENTERS_SPREAD"
1167 middle_name=trim(filename), file_position=
"REWIND", log_filename=.false.)
1168 ELSEIF (
PRESENT(iounit))
THEN
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)
1182 IF (
PRESENT(print_section))
THEN
1194 SUBROUTINE center_spread_loc(qs_env, c_loc_orb, moments)
1196 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: c_loc_orb
1197 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: moments
1199 CHARACTER(len=*),
PARAMETER :: routinen =
'center_spread_loc'
1200 INTEGER,
DIMENSION(6),
PARAMETER ::
list = (/1, 2, 3, 4, 7, 9/)
1202 INTEGER :: handle, i, iop, iorb, ispin, nao, norb, &
1204 REAL(kind=
dp) :: xmii
1205 REAL(kind=
dp),
DIMENSION(3) :: rpoint
1209 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: dipmat
1210 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_s
1211 TYPE(dbcsr_type),
POINTER :: omat, smat
1213 CALL timeset(routinen, handle)
1215 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
1216 smat => matrix_s(1, 1)%matrix
1218 nspin =
SIZE(c_loc_orb)
1221 ALLOCATE (dipmat(9))
1223 ALLOCATE (dipmat(i)%matrix)
1224 CALL dbcsr_copy(dipmat(i)%matrix, smat)
1225 CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
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)
1235 template_fmstruct=c_loc_orb(ispin)%matrix_struct)
1240 omat => dipmat(iop)%matrix
1242 CALL parallel_gemm(
'T',
'N', norb, norb, nao, 1.0_dp, c_loc_orb(ispin), ocvec, &
1247 moments(iop, iorb, ispin) = moments(iop, iorb, ispin) + xmii
1248 moments(4, iorb, ispin) = moments(4, iorb, ispin) - xmii**2
1250 moments(4, iorb, ispin) = moments(4, iorb, ispin) + xmii
1259 CALL dbcsr_release(dipmat(i)%matrix)
1260 DEALLOCATE (dipmat(i)%matrix)
1267 moments(1:3, iorb, ispin) =
pbc(moments(1:3, iorb, ispin), cell)
1271 CALL timestop(handle)
1273 END SUBROUTINE center_spread_loc
1281 SUBROUTINE center_spread_berry(qs_env, c_loc_orb, moments)
1283 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: c_loc_orb
1284 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: moments
1286 CHARACTER(len=*),
PARAMETER :: routinen =
'center_spread_berry'
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
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
1301 CALL timeset(routinen, handle)
1303 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, cell=cell)
1305 IF (cell%orthorhombic)
THEN
1310 ALLOCATE (op_sm_set(2, dim_op))
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)
1323 nspin =
SIZE(c_loc_orb, 1)
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)
1329 template_fmstruct=c_loc_orb(ispin)%matrix_struct)
1330 ALLOCATE (zij(2, dim_op))
1338 CALL parallel_gemm(
"T",
"N", norb, norb, nao, 1.0_dp, c_loc_orb(ispin), opvec, 0.0_dp, zij(j, i))
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
1359 c(idir) = c(idir) + &
1360 (cell%hmat(idir, jdir)/
twopi)*aimag(log(z))
1365 moments(1:3, istate, ispin) = cpbc(1:3)
1366 moments(4, istate, ispin) = spread_i
1367 moments(5, istate, ispin) = spread_ii
1376 CALL dbcsr_release(op_sm_set(j, i)%matrix)
1377 DEALLOCATE (op_sm_set(j, i)%matrix)
1380 DEALLOCATE (op_sm_set)
1382 CALL timestop(handle)
1384 END SUBROUTINE center_spread_berry
1396 SUBROUTINE iao_oce_expansion(qs_env, oce_basis_set_list, smat_kind, sxo, iao_coef, oce_atom, unit_nr)
1400 TYPE(dbcsr_p_type),
DIMENSION(:) :: sxo
1403 INTEGER,
INTENT(IN) :: unit_nr
1405 CHARACTER(len=*),
PARAMETER :: routinen =
'iao_oce_expansion'
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, para_group_handle
1411 INTEGER,
DIMENSION(:),
POINTER :: iao_blk_sizes, nshell, oce_blk_sizes, &
1413 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgf, last_sgf, lval
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
1419 TYPE(dbcsr_distribution_type) :: dbcsr_dist
1420 TYPE(dbcsr_type) :: iao_vec, sx_vec
1424 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1427 CALL timeset(routinen, handle)
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))
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
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, nze=0)
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, nze=0)
1447 CALL dbcsr_reserve_diag_blocks(matrix=sx_vec)
1449 nkind =
SIZE(smat_kind)
1450 ALLOCATE (sinv_kind(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
1457 CALL dbcsr_get_info(iao_vec, group=para_group_handle)
1458 CALL para_type%set_handle(para_group_handle)
1460 nspin =
SIZE(iao_coef, 1)
1461 ALLOCATE (oce_comp(natom, nspin))
1465 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, sxo(1)%matrix, iao_vec, 0.0_dp, sx_vec, &
1466 retain_sparsity=.true.)
1468 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1469 CALL dbcsr_get_block_p(matrix=sx_vec, &
1470 row=iatom, col=iatom, block=ablock, found=found)
1474 ablock = matmul(sinv_kind(ikind)%array, ablock)
1475 ALLOCATE (amat(n, m), bmat(m, m), ev(m), vec(m, m))
1476 amat(1:n, 1:m) = matmul(smat_kind(ikind)%array(1:n, 1:n), ablock(1:n, 1:m))
1477 bmat(1:m, 1:m) = matmul(transpose(ablock(1:n, 1:m)), amat(1:n, 1:m))
1478 CALL jacobi(bmat, ev, vec)
1479 oce_comp(iatom, ispin) = sum(ev(1:m))/real(m, kind=
dp)
1481 ev(i) = 1._dp/sqrt(ev(i))
1482 bmat(1:m, i) = vec(1:m, i)*ev(i)
1484 bmat(1:m, 1:m) = matmul(bmat(1:m, 1:m), transpose(vec(1:m, 1:m)))
1485 ablock(1:n, 1:m) = matmul(ablock(1:n, 1:m), bmat(1:m, 1:m))
1486 DEALLOCATE (amat, bmat, ev, vec)
1489 CALL dbcsr_replicate_all(sx_vec)
1491 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1492 CALL dbcsr_get_block_p(matrix=sx_vec, &
1493 row=iatom, col=iatom, block=ablock, found=found)
1497 ALLOCATE (oce_atom(iatom, ispin)%array(n, m))
1498 oce_atom(iatom, ispin)%array = ablock
1502 CALL para_type%sum(oce_comp)
1503 IF (unit_nr > 0)
THEN
1505 WRITE (unit_nr,
"(T4,A,I6,T30,A,2F12.4)")
"Atom:", iatom,
" Completeness: ", oce_comp(iatom, 1:nspin)
1506 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1507 oce_basis => oce_basis_set_list(ikind)%gto_basis_set
1508 CALL get_gto_basis_set(oce_basis, nset=nset, nshell=nshell, nsgf=nsgf, maxl=maxl, &
1509 l=lval, first_sgf=first_sgf, last_sgf=last_sgf)
1510 ALLOCATE (filter(nsgf, 0:maxl), vector(nsgf), prol(0:maxl, nspin))
1513 DO ishell = 1, nshell(iset)
1514 l = lval(ishell, iset)
1515 i1 = first_sgf(ishell, iset)
1516 i2 = last_sgf(ishell, iset)
1517 filter(i1:i2, l) = 1.0_dp
1521 n =
SIZE(oce_atom(iatom, 1)%array, 1)
1522 m =
SIZE(oce_atom(iatom, 1)%array, 2)
1528 vector(1:n) = oce_atom(iatom, ispin)%array(1:n, iao)*filter(1:n, l)
1529 prol(l, ispin) = sum(matmul(smat_kind(ikind)%array(1:n, 1:n), vector(1:n))*vector(1:n))
1532 WRITE (unit_nr,
"(T4,I3,T15,A,T39,(6F7.4))") iao,
" l-contributions:", (sum(prol(l, :)), l=0, maxl)
1534 DEALLOCATE (filter, vector, prol)
1539 DEALLOCATE (oce_comp)
1541 DEALLOCATE (sinv_kind(ikind)%array)
1543 DEALLOCATE (sinv_kind)
1544 DEALLOCATE (iao_blk_sizes)
1545 CALL dbcsr_release(iao_vec)
1546 CALL dbcsr_release(sx_vec)
1548 CALL timestop(handle)
1550 END SUBROUTINE iao_oce_expansion
1557 SUBROUTINE oce_overlap_matrix(smat_kind, basis_set_list)
1561 CHARACTER(len=*),
PARAMETER :: routinen =
'oce_overlap_matrix'
1563 INTEGER :: handle, ikind, iset, jset, ldsab, m1, &
1564 m2, n1, n2, ncoa, ncob, nkind, nseta, &
1566 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, npgfa, nsgfa
1567 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa
1568 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: oint, owork
1569 REAL(kind=
dp),
DIMENSION(3) :: rab
1570 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: rpgfa, scon_a, smat, zeta
1573 CALL timeset(routinen, handle)
1577 nkind =
SIZE(smat_kind)
1580 basis_set => basis_set_list(ikind)%gto_basis_set
1581 cpassert(
ASSOCIATED(basis_set))
1583 ldsab = max(m1, m2, ldsab)
1586 ALLOCATE (oint(ldsab, ldsab), owork(ldsab, ldsab))
1589 basis_set => basis_set_list(ikind)%gto_basis_set
1590 cpassert(
ASSOCIATED(basis_set))
1591 smat => smat_kind(ikind)%array
1594 first_sgfa => basis_set%first_sgf
1595 la_max => basis_set%lmax
1596 la_min => basis_set%lmin
1597 npgfa => basis_set%npgf
1598 nseta = basis_set%nset
1599 nsgfa => basis_set%nsgf_set
1600 rpgfa => basis_set%pgf_radius
1601 scon_a => basis_set%scon
1602 zeta => basis_set%zet
1606 ncoa = npgfa(iset)*
ncoset(la_max(iset))
1607 n1 = npgfa(iset)*(
ncoset(la_max(iset)) -
ncoset(la_min(iset) - 1))
1608 sgfa = first_sgfa(1, iset)
1612 ncob = npgfa(jset)*
ncoset(la_max(jset))
1613 n2 = npgfa(jset)*(
ncoset(la_max(jset)) -
ncoset(la_min(jset) - 1))
1614 sgfb = first_sgfa(1, jset)
1617 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1618 la_max(jset), la_min(jset), npgfa(jset), rpgfa(:, jset), zeta(:, jset), &
1619 rab, sab=oint(:, :))
1621 CALL contraction(oint(:, :), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1622 cb=scon_a(:, sgfb:), nb=n2, mb=nsgfa(jset), fscale=1.0_dp, trans=.false.)
1623 CALL block_add(
"IN", owork, nsgfa(iset), nsgfa(jset), smat, &
1624 sgfa, sgfb, trans=.false.)
1630 DEALLOCATE (oint, owork)
1632 CALL timestop(handle)
1634 END SUBROUTINE oce_overlap_matrix
1648 SUBROUTINE pm_localization(zij_fm_set, vectors, order, accuracy, unit_nr)
1650 TYPE(
cp_fm_type),
DIMENSION(:, :),
INTENT(IN) :: zij_fm_set
1652 INTEGER,
INTENT(IN) :: order
1653 REAL(kind=
dp),
INTENT(IN) :: accuracy
1654 INTEGER,
INTENT(IN) :: unit_nr
1656 INTEGER,
PARAMETER :: max_sweeps = 250
1658 INTEGER :: iatom, istate, jstate, natom, nstate, &
1660 REAL(kind=
dp) :: aij, bij, ct, ftarget, mii, mij, mjj, &
1661 st, t1, t2, theta, tolerance, tt
1662 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: fdiag
1665 CALL cp_fm_create(rmat, zij_fm_set(1, 1)%matrix_struct)
1669 ALLOCATE (fdiag(nstate))
1670 tolerance = 1.0e10_dp
1672 natom =
SIZE(zij_fm_set, 1)
1673 IF (unit_nr > 0)
THEN
1674 WRITE (unit_nr,
'(A,T30,A,T45,A,T63,A,T77,A)') &
1675 " Jacobi Localization ",
"Sweep",
"Functional",
"Gradient",
"Time"
1678 DO sweeps = 1, max_sweeps
1681 DO istate = 1, nstate
1682 DO jstate = istate + 1, nstate
1689 IF (order == 2)
THEN
1690 aij = aij + 4._dp*mij**2 - (mii - mjj)**2
1691 bij = bij + 4._dp*mij*(mii - mjj)
1692 ELSEIF (order == 4)
THEN
1693 aij = aij - mii**4 - mjj**4 + 6._dp*(mii**2 + mjj**2)*mij**2 + &
1694 mii**3*mjj + mii*mjj**3
1695 bij = bij + 4._dp*mij*(mii**3 - mjj**3)
1700 IF ((aij**2 + bij**2) < 1.e-10_dp) cycle
1701 tolerance = tolerance + bij**2
1702 theta = 0.25_dp*atan2(bij, -aij)
1709 CALL cp_fm_rot_cols(zij_fm_set(iatom, 1), istate, jstate, ct, st)
1710 CALL cp_fm_rot_rows(zij_fm_set(iatom, 1), istate, jstate, ct, st)
1717 ftarget = ftarget + sum(fdiag**order)
1719 ftarget = ftarget**(1._dp/order)
1720 tolerance = sqrt(tolerance)
1723 IF (unit_nr > 0)
THEN
1724 WRITE (unit_nr,
'(T31,I4,T39,F16.8,T55,G16.8,T71,F10.2)') sweeps, ftarget, tolerance, tt
1726 IF (tolerance < accuracy)
EXIT
1733 END SUBROUTINE pm_localization
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
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....
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.
subroutine, public create_oce_basis(oce_basis, orb_basis, lmax_oce, nbas_oce)
...
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.
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)
...
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.
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...
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
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, use_sp)
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, stride, 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.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Machine interface based on Fortran 2003 and POSIX.
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition of mathematical constants and functions.
real(kind=dp), parameter, public twopi
Collection of simple mathematical functions and subroutines.
subroutine, public jacobi(a, d, v)
Jacobi matrix diagonalization. The eigenvalues are returned in vector d and the eigenvectors are retu...
subroutine, public invmat_symm(a, cholesky_triangle)
returns inverse of real symmetric, positive definite matrix
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)
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, external_vector)
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, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
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, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, 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, 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_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, 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.
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>
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
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.
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.
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.
represent a pointer to a 2d array
keeps the information about the structure of a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represent a list of objects
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 ...