Skip to content

Commit

Permalink
More refactoring, adding efield output (#2742)
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Apr 20, 2023
1 parent 056b934 commit e08f4c8
Show file tree
Hide file tree
Showing 7 changed files with 218 additions and 63 deletions.
1 change: 1 addition & 0 deletions src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -267,6 +267,7 @@ list(
input_restart_rng.F
integration_grid_types.F
iao_analysis.F
iao_types.F
ipi_driver.F
iterate_matrix.F
kg_correction.F
Expand Down
87 changes: 28 additions & 59 deletions src/iao_analysis.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,9 +59,8 @@ MODULE iao_analysis
dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
dbcsr_replicate_all, dbcsr_reserve_diag_blocks, dbcsr_set, dbcsr_trace, dbcsr_type, &
dbcsr_type_no_symmetry
USE iao_types, ONLY: iao_env_type
USE input_section_types, ONLY: section_get_ivals,&
section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: default_path_length,&
Expand Down Expand Up @@ -130,31 +129,33 @@ MODULE iao_analysis
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param input_section ...
!> \param iao_env ...
!> \param unit_nr ...
! **************************************************************************************************
SUBROUTINE iao_wfn_analysis(qs_env, input_section, unit_nr)
SUBROUTINE iao_wfn_analysis(qs_env, iao_env, unit_nr)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(section_vals_type), POINTER :: input_section
TYPE(iao_env_type), INTENT(IN) :: iao_env
INTEGER, INTENT(IN) :: unit_nr

CHARACTER(len=*), PARAMETER :: routineN = 'iao_wfn_analysis'

CHARACTER(LEN=2) :: element_symbol
CHARACTER(LEN=default_string_length) :: bname
INTEGER :: handle, i, iatom, ikind, isgf, ispin, lb, lmax_oce, mao, nao, natom, nbas_oce, &
ne, nimages, nkind, norb, nref, nsgf, nspin, nvec, sweeps, ub, z
INTEGER :: handle, i, iatom, ikind, isgf, ispin, &
lb, mao, nao, natom, ne, nimages, &
nkind, norb, nref, nsgf, nspin, nvec, &
sweeps, ub, z
INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf, nsgfat
INTEGER, DIMENSION(0:3) :: elcon
INTEGER, DIMENSION(2) :: nocc
INTEGER, DIMENSION(:), POINTER :: econf
LOGICAL :: append_oce, cubes_iao, cubes_ibo, do_bondorbitals, do_center, do_charges, do_oce, &
explicit, molden_iao, molden_ibo, pos_periodic, uniform_occupation
REAL(KIND=dp) :: eps_occ, eps_svd, trace, zval
LOGICAL :: cubes_iao, cubes_ibo, do_bondorbitals, &
molden_iao, molden_ibo, &
uniform_occupation
REAL(KIND=dp) :: trace, zval
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mcharge
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: moments
REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers
TYPE(cell_type), POINTER :: cell
TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: smat_kind
TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:, :) :: oce_atom
TYPE(cp_fm_struct_type), POINTER :: fm_struct
Expand All @@ -180,11 +181,10 @@ SUBROUTINE iao_wfn_analysis(qs_env, input_section, unit_nr)
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(section_vals_type), POINTER :: iao_cubes_section, iao_molden_section, &
ibo_cc_section, ibo_cubes_section, &
ibo_molden_section, subsection
ibo_molden_section

! only do IAO analysis if explicitly requested
CALL section_vals_get(input_section, explicit=explicit)
IF (.NOT. explicit) RETURN
CPASSERT(iao_env%do_iao)

! k-points?
CALL get_qs_env(qs_env, dft_control=dft_control)
Expand All @@ -208,43 +208,12 @@ SUBROUTINE iao_wfn_analysis(qs_env, input_section, unit_nr)
END IF
CALL cite_reference(Knizia2013)

! input options
CALL section_vals_val_get(input_section, "EPS_SVD", r_val=eps_svd)
CALL section_vals_val_get(input_section, "EPS_OCC", r_val=eps_occ)
CALL section_vals_val_get(input_section, "ATOMIC_CHARGES", l_val=do_charges)
iao_molden_section => section_vals_get_subs_vals(input_section, "IAO_MOLDEN")
CALL section_vals_get(iao_molden_section, explicit=molden_iao)
iao_cubes_section => section_vals_get_subs_vals(input_section, "IAO_CUBES")
CALL section_vals_get(iao_cubes_section, explicit=cubes_iao)
subsection => section_vals_get_subs_vals(input_section, "ONE_CENTER_EXPANSION")
CALL section_vals_get(subsection, explicit=do_oce)
IF (do_oce) THEN
subsection => section_vals_get_subs_vals(input_section, "ONE_CENTER_EXPANSION")
CALL section_vals_val_get(subsection, "LMAX", i_val=lmax_oce)
CALL section_vals_val_get(subsection, "NBAS", i_val=nbas_oce)
CALL section_vals_val_get(subsection, "APPEND", l_val=append_oce)
END IF
subsection => section_vals_get_subs_vals(input_section, "BOND_ORBITALS")
CALL section_vals_get(subsection, explicit=do_bondorbitals)
IF (do_bondorbitals) THEN
subsection => section_vals_get_subs_vals(input_section, "BOND_ORBITALS")
ibo_molden_section => section_vals_get_subs_vals(subsection, "IBO_MOLDEN")
CALL section_vals_get(ibo_molden_section, explicit=molden_ibo)
ibo_cubes_section => section_vals_get_subs_vals(subsection, "IBO_CUBES")
CALL section_vals_get(ibo_cubes_section, explicit=cubes_ibo)
ibo_cc_section => section_vals_get_subs_vals(subsection, "CHARGE_CENTER")
CALL section_vals_get(ibo_cc_section, explicit=do_center)
IF (do_center) THEN
CALL get_qs_env(qs_env, cell=cell)
CALL section_vals_val_get(ibo_cc_section, "POSITION_OPERATOR_BERRY", &
l_val=pos_periodic, explicit=explicit)
IF (.NOT. explicit) THEN
! set default according to cell periodicity
pos_periodic = .TRUE.
IF (ALL(cell%perd == 0)) pos_periodic = .FALSE.
END IF
END IF
END IF
! input sections
iao_molden_section => iao_env%iao_molden_section
iao_cubes_section => iao_env%iao_cubes_section
ibo_molden_section => iao_env%ibo_molden_section
ibo_cubes_section => iao_env%ibo_cubes_section
ibo_cc_section => iao_env%ibo_cc_section

! check for or generate reference basis
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
Expand Down Expand Up @@ -318,7 +287,7 @@ SUBROUTINE iao_wfn_analysis(qs_env, input_section, unit_nr)
CALL cp_fm_create(p_orb_ref, fm_struct)
CALL cp_fm_struct_release(fm_struct)
!
CALL iao_projectors(smat, sref(1)%matrix, sro(1)%matrix, p_orb_ref, p_ref_orb, eps_svd)
CALL iao_projectors(smat, sref(1)%matrix, sro(1)%matrix, p_orb_ref, p_ref_orb, iao_env%eps_svd)

! occupied orbitals
nocc(1:2) = 0
Expand All @@ -330,7 +299,7 @@ SUBROUTINE iao_wfn_analysis(qs_env, input_section, unit_nr)
ELSE
nvec = 0
DO i = 1, norb
IF (occupation_numbers(i) > eps_occ) THEN
IF (occupation_numbers(i) > iao_env%eps_occ) THEN
nvec = i
ELSE
EXIT
Expand Down Expand Up @@ -375,7 +344,7 @@ SUBROUTINE iao_wfn_analysis(qs_env, input_section, unit_nr)
CALL cp_fm_release(smo)
END DO
!
IF (do_charges) THEN
IF (iao_env%do_charges) THEN
! population analysis
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(/,T2,A)') 'Intrinsic AO Population Analysis '
Expand Down Expand Up @@ -405,7 +374,7 @@ SUBROUTINE iao_wfn_analysis(qs_env, input_section, unit_nr)
DEALLOCATE (mcharge)
END IF
!
IF (do_oce) THEN
IF (iao_env%do_oce) THEN
! generate OCE basis
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A)') "IAO One-Center Expansion: OCE Basis Set"
Expand All @@ -415,7 +384,7 @@ SUBROUTINE iao_wfn_analysis(qs_env, input_section, unit_nr)
orbbasis => orb_basis_set_list(ikind)%gto_basis_set
CPASSERT(ASSOCIATED(orbbasis))
NULLIFY (ocebasis)
CALL create_oce_basis(ocebasis, orbbasis, lmax_oce, nbas_oce)
CALL create_oce_basis(ocebasis, orbbasis, iao_env%lmax_oce, iao_env%nbas_oce)
CALL init_orb_basis_set(ocebasis)
CALL init_interaction_radii_orb_basis(ocebasis, dft_control%qs_control%eps_pgf_orb)
oce_basis_set_list(ikind)%gto_basis_set => ocebasis
Expand Down Expand Up @@ -557,10 +526,10 @@ SUBROUTINE iao_wfn_analysis(qs_env, input_section, unit_nr)
END DO
DEALLOCATE (cloc)
!
IF (do_center) THEN
IF (iao_env%do_center) THEN
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A)') ' Calculate Localized Orbital Centers and Spread'
IF (pos_periodic) THEN
IF (iao_env%pos_periodic) THEN
WRITE (unit_nr, '(T2,A)') ' Use Berry Phase Position Operator'
ELSE
WRITE (unit_nr, '(T2,A)') ' Use Local Position Operator'
Expand All @@ -570,7 +539,7 @@ SUBROUTINE iao_wfn_analysis(qs_env, input_section, unit_nr)
! x y z m2(1) m2(2)
ALLOCATE (moments(5, nvec, nspin))
moments = 0.0_dp
IF (pos_periodic) THEN
IF (iao_env%pos_periodic) THEN
CALL center_spread_berry(qs_env, c_loc_orb, moments)
ELSE
CALL center_spread_loc(qs_env, c_loc_orb, moments)
Expand Down
161 changes: 161 additions & 0 deletions src/iao_types.F
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,161 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2023 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief Calculate ntrinsic atomic orbitals and analyze wavefunctions
!> \par History
!> 03.2023 created [JGH]
!> \author JGH
! **************************************************************************************************
MODULE iao_types
USE cell_types, ONLY: cell_type
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
#include "./base/base_uses.f90"

IMPLICIT NONE
PRIVATE

PUBLIC :: iao_env_type, iao_read_input, iao_set_default

! **************************************************************************************************
TYPE iao_env_type
LOGICAL :: do_iao
!
REAL(KIND=dp) :: eps_svd
REAL(KIND=dp) :: eps_occ
! chages
LOGICAL :: do_charges
! one-center expansion
LOGICAL :: do_oce
INTEGER :: lmax_oce
INTEGER :: nbas_oce
LOGICAL :: append_oce
! Bond orbitals
LOGICAL :: do_bondorbitals
! Wannier centers
LOGICAL :: do_center
LOGICAL :: pos_periodic
! Molden
LOGICAL :: molden_iao
LOGICAL :: molden_ibo
! CUBE files
LOGICAL :: cubes_iao
LOGICAL :: cubes_ibo
! Input sections
TYPE(section_vals_type), POINTER :: iao_cubes_section, iao_molden_section, &
ibo_cubes_section, ibo_molden_section, &
ibo_cc_section
END TYPE iao_env_type
! **************************************************************************************************

CONTAINS

! **************************************************************************************************
!> \brief ...
!> \param iao_env ...
! **************************************************************************************************
SUBROUTINE iao_set_default(iao_env)
TYPE(iao_env_type), INTENT(INOUT) :: iao_env

!iao
iao_env%do_iao = .FALSE.
iao_env%eps_svd = 0.0_dp
iao_env%eps_occ = 0.0_dp
! charges
iao_env%do_charges = .FALSE.
! one-center expansion
iao_env%do_oce = .FALSE.
iao_env%lmax_oce = 3
iao_env%nbas_oce = 10
iao_env%append_oce = .FALSE.
! Bond orbitals
iao_env%do_bondorbitals = .FALSE.
! Wannier centers
iao_env%do_center = .FALSE.
iao_env%pos_periodic = .FALSE.
! i/o
iao_env%molden_iao = .FALSE.
iao_env%molden_ibo = .FALSE.
iao_env%cubes_iao = .FALSE.
iao_env%cubes_ibo = .FALSE.
! Input sections
NULLIFY (iao_env%iao_cubes_section, iao_env%iao_molden_section)
NULLIFY (iao_env%ibo_cubes_section, iao_env%ibo_molden_section)
NULLIFY (iao_env%ibo_cc_section)

END SUBROUTINE iao_set_default

! **************************************************************************************************

! **************************************************************************************************
!> \brief ...
!> \param iao_env ...
!> \param iao_section ...
!> \param cell ...
! **************************************************************************************************
SUBROUTINE iao_read_input(iao_env, iao_section, cell)
TYPE(iao_env_type), INTENT(INOUT) :: iao_env
TYPE(section_vals_type), POINTER :: iao_section
TYPE(cell_type), OPTIONAL :: cell

LOGICAL :: explicit, iao_explicit
TYPE(section_vals_type), POINTER :: subsection

CALL iao_set_default(iao_env)

CALL section_vals_get(iao_section, explicit=iao_explicit)
IF (iao_explicit) THEN
iao_env%do_iao = .TRUE.
! input options
CALL section_vals_val_get(iao_section, "EPS_SVD", r_val=iao_env%eps_svd)
CALL section_vals_val_get(iao_section, "EPS_OCC", r_val=iao_env%eps_occ)
CALL section_vals_val_get(iao_section, "ATOMIC_CHARGES", l_val=iao_env%do_charges)
iao_env%iao_molden_section => section_vals_get_subs_vals(iao_section, "IAO_MOLDEN")
CALL section_vals_get(iao_env%iao_molden_section, explicit=iao_env%molden_iao)
iao_env%iao_cubes_section => section_vals_get_subs_vals(iao_section, "IAO_CUBES")
CALL section_vals_get(iao_env%iao_cubes_section, explicit=iao_env%cubes_iao)
subsection => section_vals_get_subs_vals(iao_section, "ONE_CENTER_EXPANSION")
CALL section_vals_get(subsection, explicit=iao_env%do_oce)
IF (iao_env%do_oce) THEN
subsection => section_vals_get_subs_vals(iao_section, "ONE_CENTER_EXPANSION")
CALL section_vals_val_get(subsection, "LMAX", i_val=iao_env%lmax_oce)
CALL section_vals_val_get(subsection, "NBAS", i_val=iao_env%nbas_oce)
CALL section_vals_val_get(subsection, "APPEND", l_val=iao_env%append_oce)
END IF
subsection => section_vals_get_subs_vals(iao_section, "BOND_ORBITALS")
CALL section_vals_get(subsection, explicit=iao_env%do_bondorbitals)
IF (iao_env%do_bondorbitals) THEN
subsection => section_vals_get_subs_vals(iao_section, "BOND_ORBITALS")
iao_env%ibo_molden_section => section_vals_get_subs_vals(subsection, "IBO_MOLDEN")
CALL section_vals_get(iao_env%ibo_molden_section, explicit=iao_env%molden_ibo)
iao_env%ibo_cubes_section => section_vals_get_subs_vals(subsection, "IBO_CUBES")
CALL section_vals_get(iao_env%ibo_cubes_section, explicit=iao_env%cubes_ibo)
iao_env%ibo_cc_section => section_vals_get_subs_vals(subsection, "CHARGE_CENTER")
CALL section_vals_get(iao_env%ibo_cc_section, explicit=iao_env%do_center)
IF (iao_env%do_center) THEN
CALL section_vals_val_get(iao_env%ibo_cc_section, "POSITION_OPERATOR_BERRY", &
l_val=iao_env%pos_periodic, explicit=explicit)
IF (.NOT. explicit) THEN
! set default according to cell periodicity
iao_env%pos_periodic = .TRUE.
IF (PRESENT(cell)) THEN
IF (ALL(cell%perd == 0)) iao_env%pos_periodic = .FALSE.
END IF
END IF
END IF
END IF
END IF

END SUBROUTINE iao_read_input

! **************************************************************************************************

END MODULE iao_types
18 changes: 15 additions & 3 deletions src/qs_efield_berry.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -669,9 +669,9 @@ SUBROUTINE qs_dispfield_derivatives(qs_env, just_energy, calculate_forces)
CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_dispfield_derivatives'

COMPLEX(dp) :: zdet, zdeta, zi(3)
INTEGER :: handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, jatom, jkind, jset, &
ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, nseta, nsetb, sgfa, &
sgfb
INTEGER :: handle, i, ia, iatom, icol, idir, ikind, iodeb, irow, iset, ispin, jatom, jkind, &
jset, ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, nseta, nsetb, &
sgfa, sgfb
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
npgfb, nsgfa, nsgfb
Expand Down Expand Up @@ -1093,6 +1093,18 @@ SUBROUTINE qs_dispfield_derivatives(qs_env, just_energy, calculate_forces)
END DO
energy%efield = 0.25_dp*omega/twopi*ener_field

! debugging output
IF (para_env%is_source()) THEN
iodeb = -1
IF (iodeb > 0) THEN
WRITE (iodeb, '(A,T61,F20.10)') " Polarisation Quantum: ", 2._dp*twopi*twopi*hmat(3, 3)
WRITE (iodeb, '(A,T21,3F20.10)') " Polarisation: ", 2._dp*twopi*ci(1:3)
WRITE (iodeb, '(A,T21,3F20.10)') " Displacement: ", fieldpol(1:3)
WRITE (iodeb, '(A,T21,3F20.10)') " E-Field: ", ((fieldpol(i) - 2._dp*twopi*ci(i)), i=1, 3)
WRITE (iodeb, '(A,T61,F20.10)') " Disp Free Energy:", energy%efield
END IF
END IF

IF (.NOT. just_energy) THEN
DO i = 1, 3
di(i) = -omega*(fieldpol(i) - 2._dp*twopi*ci(i))*dfilter(i)
Expand Down

0 comments on commit e08f4c8

Please sign in to comment.