comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Locpowf (#2154)

Browse source 
implementing partially occupied wannier states [sghosh]
This commit is contained in:
hforbert 2022-06-07 16:29:57 +02:00 • committed by GitHub
parent aaa3bad4b0
commit 278fffb80f
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
19 changed files with 2222 additions and 387 deletions
Download patch file Download diff file Expand all files Collapse all files

4
src/almo_scf.F
View file

@ -87,10 +87,10 @@ MODULE almo_scf
USE mscfg_types, ONLY: get_matrix_from_submatrices,&
molecular_scf_guess_env_type
USE particle_types, ONLY: particle_type
USE qs_atomic_block, ONLY: calculate_atomic_block_dm
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_initial_guess, ONLY: calculate_atomic_block_dm,&
calculate_mopac_dm
USE qs_initial_guess, ONLY: calculate_mopac_dm
USE qs_kind_types, ONLY: qs_kind_type
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_p_type

4
src/dm_ls_scf_qs.F
View file

@ -50,6 +50,7 @@ MODULE dm_ls_scf_qs
REALSPACE,&
RECIPROCALSPACE,&
pw_p_type
USE qs_atomic_block, ONLY: calculate_atomic_block_dm
USE qs_collocate_density, ONLY: calculate_rho_elec
USE qs_density_mixing_types, ONLY: direct_mixing_nr,&
gspace_mixing_nr
@ -57,8 +58,7 @@ MODULE dm_ls_scf_qs
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_gspace_mixing, ONLY: gspace_mixing
USE qs_initial_guess, ONLY: calculate_atomic_block_dm,&
calculate_mopac_dm
USE qs_initial_guess, ONLY: calculate_mopac_dm
USE qs_kind_types, ONLY: qs_kind_type
USE qs_ks_methods, ONLY: qs_ks_update_qs_env
USE qs_ks_types, ONLY: qs_ks_did_change,&

2
src/energy_corrections.F
View file

@ -126,6 +126,7 @@ MODULE energy_corrections
RECIPROCALSPACE,&
pw_p_type,&
pw_type
USE qs_atomic_block, ONLY: calculate_atomic_block_dm
USE qs_collocate_density, ONLY: calculate_rho_elec
USE qs_core_energies, ONLY: calculate_ecore_overlap,&
calculate_ptrace
@ -139,7 +140,6 @@ MODULE energy_corrections
USE qs_force_types, ONLY: qs_force_type,&
total_qs_force,&
zero_qs_force
USE qs_initial_guess, ONLY: calculate_atomic_block_dm
USE qs_integrate_potential, ONLY: integrate_v_core_rspace,&
integrate_v_rspace
USE qs_kind_types, ONLY: get_qs_kind,&

140
src/fm/cp_fm_basic_linalg.F
View file

@ -58,7 +58,8 @@ MODULE cp_fm_basic_linalg
cp_fm_qr_factorization, & ! compute the QR factorization of a rectangular matrix
cp_fm_solve, & ! solves the equation A*B=C A and C are input
cp_fm_pdgeqpf, & ! compute a QR factorization with column pivoting of a M-by-N distributed matrix
cp_fm_pdorgqr ! generates an M-by-N as first N columns of a product of K elementary reflectors
cp_fm_pdorgqr, & ! generates an M-by-N as first N columns of a product of K elementary reflectors
cp_fm_Gram_Schmidt_orthonorm ! Gram-Schmidt orthonormalization of columns of a full matrix
REAL(KIND=dp), EXTERNAL :: dlange, pdlange, pdlatra
REAL(KIND=sp), EXTERNAL :: slange, pslange, pslatra
@ -2285,4 +2286,141 @@ CONTAINS
END SUBROUTINE cp_fm_pdorgqr
! **************************************************************************************************
!> \brief Orthonormalizes selected rows and columns of a full matrix, matrix_a
!> \param matrix_a ...
!> \param B ...
!> \param nrows number of rows of matrix_a, optional, defaults to size(matrix_a,1)
!> \param ncols number of columns of matrix_a, optional, defaults to size(matrix_a, 2)
!> \param start_row starting index of rows, optional, defaults to 1
!> \param start_col starting index of columns, optional, defaults to 1
!> \param do_norm ...
!> \param do_print ...
! **************************************************************************************************
SUBROUTINE cp_fm_Gram_Schmidt_orthonorm(matrix_a, B, nrows, ncols, start_row, start_col, &
do_norm, do_print)
TYPE(cp_fm_type), INTENT(IN), POINTER :: matrix_a
REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: B
INTEGER, INTENT(IN), OPTIONAL :: nrows, ncols, start_row, start_col
LOGICAL, INTENT(IN), OPTIONAL :: do_norm, do_print
CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_Gram_Schmidt_orthonorm', &
routineP = moduleN//':'//routineN
INTEGER :: end_col_global, end_col_local, end_row_global, end_row_local, handle, i, j, &
j_col, ncol_global, ncol_local, nrow_global, nrow_local, start_col_global, &
start_col_local, start_row_global, start_row_local, this_col
INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
LOGICAL :: my_do_norm, my_do_print
REAL(KIND=dp) :: norm
REAL(kind=dp), DIMENSION(:, :), POINTER :: a
CALL timeset(routineN, handle)
CPASSERT(ASSOCIATED(matrix_a))
CPASSERT(matrix_a%ref_count .GT. 0)
my_do_norm = .TRUE.
IF (PRESENT(do_norm)) my_do_norm = do_norm
my_do_print = .FALSE.
IF (PRESENT(do_print) .AND. (my_do_norm)) my_do_print = do_print
IF (SIZE(B) /= 0) THEN
IF (PRESENT(nrows)) THEN
nrow_global = nrows
ELSE
nrow_global = SIZE(B, 1)
END IF
IF (PRESENT(ncols)) THEN
ncol_global = ncols
ELSE
ncol_global = SIZE(B, 2)
END IF
IF (PRESENT(start_row)) THEN
start_row_global = start_row
ELSE
start_row_global = 1
END IF
IF (PRESENT(start_col)) THEN
start_col_global = start_col
ELSE
start_col_global = 1
END IF
end_row_global = start_row_global + nrow_global - 1
end_col_global = start_col_global + ncol_global - 1
CALL cp_fm_get_info(matrix=matrix_a, &
nrow_global=nrow_global, ncol_global=ncol_global, &
nrow_local=nrow_local, ncol_local=ncol_local, &
row_indices=row_indices, col_indices=col_indices)
IF (end_row_global > nrow_global) THEN
end_row_global = nrow_global
END IF
IF (end_col_global > ncol_global) THEN
end_col_global = ncol_global
END IF
! find out row/column indices of locally stored matrix elements that
! needs to be copied.
! Arrays row_indices and col_indices are assumed to be sorted in
! ascending order
DO start_row_local = 1, nrow_local
IF (row_indices(start_row_local) >= start_row_global) EXIT
END DO
DO end_row_local = start_row_local, nrow_local
IF (row_indices(end_row_local) > end_row_global) EXIT
END DO
end_row_local = end_row_local - 1
DO start_col_local = 1, ncol_local
IF (col_indices(start_col_local) >= start_col_global) EXIT
END DO
DO end_col_local = start_col_local, ncol_local
IF (col_indices(end_col_local) > end_col_global) EXIT
END DO
end_col_local = end_col_local - 1
a => matrix_a%local_data
this_col = col_indices(start_col_local) - start_col_global + 1
B(:, this_col) = a(:, start_col_local)
IF (my_do_norm) THEN
norm = SQRT(accurate_dot_product(B(:, this_col), B(:, this_col)))
B(:, this_col) = B(:, this_col)/norm
IF (my_do_print) WRITE (*, '(I3,F8.3)') this_col, norm
END IF
DO i = start_col_local + 1, end_col_local
this_col = col_indices(i) - start_col_global + 1
B(:, this_col) = a(:, i)
DO j = start_col_local, i - 1
j_col = col_indices(j) - start_col_global + 1
B(:, this_col) = B(:, this_col) - &
accurate_dot_product(B(:, j_col), B(:, this_col))* &
B(:, j_col)/accurate_dot_product(B(:, j_col), B(:, j_col))
END DO
IF (my_do_norm) THEN
norm = SQRT(accurate_dot_product(B(:, this_col), B(:, this_col)))
B(:, this_col) = B(:, this_col)/norm
IF (my_do_print) WRITE (*, '(I3,F8.3)') this_col, norm
END IF
END DO
CALL mp_sum(B, matrix_a%matrix_struct%para_env%group)
END IF
CALL timestop(handle)
END SUBROUTINE cp_fm_Gram_Schmidt_orthonorm
END MODULE cp_fm_basic_linalg

16
src/input_constants.F
View file

@ -448,7 +448,15 @@ MODULE input_constants
do_loc_crazy = 2, &
do_loc_direct = 3, &
do_loc_l1_norm_sd = 4, &
do_loc_scdm = 5
do_loc_scdm = 5, &
do_loc_gapo = 6
INTEGER, PARAMETER, PUBLIC :: do_loc_cpo_atomic = 0, &
do_loc_cpo_restart = 1, &
do_loc_cpo_random = 2
INTEGER, PARAMETER, PUBLIC :: do_loc_cpo_space_wan = 0, &
do_loc_cpo_space_nmo = 1
INTEGER, PARAMETER, PUBLIC :: do_loc_min = 0, &
do_loc_max = 1, &
@ -460,11 +468,13 @@ MODULE input_constants
state_loc_range = 1, &
state_loc_list = 2, &
energy_loc_range = 3, &
state_loc_none = 4
state_loc_none = 4, &
state_loc_mixed = 5
INTEGER, PARAMETER, PUBLIC :: do_loc_homo = 0, &
do_loc_lumo = 1, &
do_loc_both = 2
do_loc_both = 2, &
do_loc_mixed = 3
INTEGER, PARAMETER, PUBLIC :: orb_s = 0, &
orb_px = 1, &

63
src/input_cp2k_loc.F
View file

@ -12,9 +12,10 @@ MODULE input_cp2k_loc
high_print_level,&
low_print_level
USE input_constants, ONLY: &
do_loc_both, do_loc_crazy, do_loc_direct, do_loc_homo, do_loc_jacobi, do_loc_l1_norm_sd, &
do_loc_lumo, do_loc_max, do_loc_min, do_loc_none, do_loc_scdm, op_loc_berry, op_loc_boys, &
op_loc_pipek
do_loc_both, do_loc_cpo_atomic, do_loc_cpo_random, do_loc_cpo_restart, &
do_loc_cpo_space_nmo, do_loc_cpo_space_wan, do_loc_crazy, do_loc_direct, do_loc_gapo, &
do_loc_homo, do_loc_jacobi, do_loc_l1_norm_sd, do_loc_lumo, do_loc_max, do_loc_min, &
do_loc_mixed, do_loc_none, do_loc_scdm, op_loc_berry, op_loc_boys, op_loc_pipek
USE input_cp2k_mm, ONLY: create_dipoles_section
USE input_cp2k_motion_print, ONLY: add_format_keyword
USE input_keyword_types, ONLY: keyword_create,&
@ -134,22 +135,51 @@ CONTAINS
CALL keyword_create( &
keyword, __LOCATION__, name="METHOD", &
description="Method of optimization if any", &
usage="METHOD (JACOBI|CRAZY|DIRECT|L1SD|SCDM|NONE)", &
enum_c_vals=s2a("NONE", "JACOBI", "CRAZY", "L1SD", "DIRECT", "SCDM"), &
usage="METHOD (JACOBI|CRAZY|DIRECT|GAPO|L1SD|SCDM|NONE)", &
enum_c_vals=s2a("NONE", "JACOBI", "CRAZY", "GAPO", "L1SD", "DIRECT", "SCDM"), &
enum_i_vals=(/do_loc_none, &
do_loc_jacobi, &
do_loc_crazy, &
do_loc_gapo, &
do_loc_l1_norm_sd, &
do_loc_direct, do_loc_scdm/), &
enum_desc=s2a("No localization is applied", &
"Using 2 x 2 rotations of the orbitals, slow but robust", &
"A new fast method is applied, might be slightly less robust than jacobi, but usually much faster", &
"Gradient ascent for partially occupied wannier functions", &
"Steepest descent minimization of an approximate l1 norm", &
"Using a direct minimisation approacha", "Use QR factorization"), &
default_i_val=do_loc_jacobi)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="CPO_GUESS", &
description="Initial guess for coefficients if METHOD GAPO is used", &
usage="CPO_GUESS (ATOMIC|RESTART|RANDOM)", &
enum_c_vals=s2a("ATOMIC", "RESTART", "RANDOM"), &
enum_i_vals=(/do_loc_cpo_atomic, do_loc_cpo_restart, do_loc_cpo_random/), &
default_i_val=do_loc_cpo_atomic)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="CPO_GUESS_SPACE", &
description="Orbital space from which initial guess for coefficients is determined "// &
"if METHOD GAPO and CPO_GUESS ATOMIC are employed", &
usage="CPO_GUESS_SPACE (WAN|ALL)", &
enum_c_vals=s2a("WAN", "ALL"), &
enum_i_vals=(/do_loc_cpo_space_wan, do_loc_cpo_space_nmo/), &
default_i_val=do_loc_cpo_space_wan)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="CG_PO", &
description="Use conjugate gradient in conjunction with METHOD GAPO. If FALSE, "// &
" steepest descent is used instead.", &
usage="CG_PO", default_l_val=.TRUE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="JACOBI_FALLBACK", &
description="Use Jacobi method in case no convergence was achieved"// &
" by using the crazy rotations method.", &
@ -181,6 +211,14 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="LOCMIXD_RESTART_FILE_NAME", &
description="File name where to read the MOS from"// &
"which to restart the localization procedure for MIXED states", &
usage="LOCMIXD_RESTART_FILE_NAME <FILENAME>", &
type_of_var=lchar_t)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="LOCLUMO_RESTART_FILE_NAME", &
description="File name where to read the MOS from"// &
"which to restart the localization procedure for unoccupied states", &
@ -218,11 +256,20 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="NEXTRA", &
description="Number of orbitals above fully occupied MOs to be localized, "// &
"up to now only valid in combination with GPW. "// &
"This keyword has to be present for STATES MIXED option. "// &
"Otherwise, only the fully occupied MOs are localized.", &
usage="NEXTRA 5", default_i_val=0)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="STATES", &
description="Which states to localize, LUMO up to now only available in GPW", &
usage="STATES (HOMO|LUMO|ALL)", &
enum_c_vals=s2a("OCCUPIED", "UNOCCUPIED", "ALL"), &
enum_i_vals=(/do_loc_homo, do_loc_lumo, do_loc_both/), &
usage="STATES (HOMO|LUMO|MIXED|ALL)", &
enum_c_vals=s2a("OCCUPIED", "UNOCCUPIED", "MIXED", "ALL"), &
enum_i_vals=(/do_loc_homo, do_loc_lumo, do_loc_mixed, do_loc_both/), &
default_i_val=do_loc_homo)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

186
src/qs_atomic_block.F Normal file
View file

@ -0,0 +1,186 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2022 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Routine to return block diagonal density matrix. Blocks correspond to the atomic densities
!> \par History
!> 2006.03 Moved here from qs_scf.F [Joost VandeVondele]
!> 2022.05 split from qs_initial_guess.F to break circular dependency [Harald Forbert]
! **************************************************************************************************
MODULE qs_atomic_block
USE atom_kind_orbitals, ONLY: calculate_atomic_orbitals
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind_set
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: &
dbcsr_add_on_diag, dbcsr_dot, dbcsr_get_info, dbcsr_iterator_blocks_left, &
dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type
USE kinds, ONLY: dp
USE message_passing, ONLY: mp_sum
USE particle_types, ONLY: particle_type
USE qs_kind_types, ONLY: qs_kind_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_atomic_block'
PUBLIC :: calculate_atomic_block_dm
TYPE atom_matrix_type
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mat
END TYPE atom_matrix_type
CONTAINS
! **************************************************************************************************
!> \brief returns a block diagonal density matrix. Blocks correspond to the atomic densities.
!> \param pmatrix ...
!> \param matrix_s ...
!> \param particle_set ...
!> \param atomic_kind_set ...
!> \param qs_kind_set ...
!> \param nspin ...
!> \param nelectron_spin ...
!> \param ounit ...
!> \param para_env ...
! **************************************************************************************************
SUBROUTINE calculate_atomic_block_dm(pmatrix, matrix_s, particle_set, atomic_kind_set, &
qs_kind_set, nspin, nelectron_spin, ounit, para_env)
TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: pmatrix
TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
INTEGER, INTENT(IN) :: nspin
INTEGER, DIMENSION(:), INTENT(IN) :: nelectron_spin
INTEGER, INTENT(IN) :: ounit
TYPE(cp_para_env_type) :: para_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atomic_block_dm'
INTEGER :: blk, group, handle, icol, ikind, irow, &
ispin, natom, nc, nkind, nocc(2)
INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nok
REAL(dp), DIMENSION(:, :), POINTER :: pdata
REAL(KIND=dp) :: rds, rscale, trps1
TYPE(atom_matrix_type), ALLOCATABLE, DIMENSION(:) :: pmat
TYPE(atomic_kind_type), POINTER :: atomic_kind
TYPE(dbcsr_iterator_type) :: iter
TYPE(dbcsr_type), POINTER :: matrix_p
TYPE(qs_kind_type), POINTER :: qs_kind
CALL timeset(routineN, handle)
natom = SIZE(particle_set)
nkind = SIZE(atomic_kind_set)
ALLOCATE (kind_of(natom))
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
ALLOCATE (pmat(nkind))
ALLOCATE (nok(2, nkind))
! precompute the atomic blocks corresponding to spherical atoms
DO ikind = 1, nkind
atomic_kind => atomic_kind_set(ikind)
qs_kind => qs_kind_set(ikind)
NULLIFY (pmat(ikind)%mat)
IF (ounit > 0) THEN
WRITE (UNIT=ounit, FMT="(/,T2,A)") &
"Guess for atomic kind: "//TRIM(atomic_kind%name)
END IF
CALL calculate_atomic_orbitals(atomic_kind, qs_kind, iunit=ounit, &
pmat=pmat(ikind)%mat, nocc=nocc)
nok(1:2, ikind) = nocc(1:2)
END DO
rscale = 1.0_dp
IF (nspin == 2) rscale = 0.5_dp
DO ispin = 1, nspin
IF ((ounit > 0) .AND. (nspin > 1)) THEN
WRITE (UNIT=ounit, FMT="(/,T2,A,I0)") "Spin ", ispin
END IF
matrix_p => pmatrix(ispin)%matrix
CALL dbcsr_set(matrix_p, 0.0_dp)
nocc(ispin) = 0
CALL dbcsr_iterator_start(iter, matrix_p)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, irow, icol, pdata, blk)
ikind = kind_of(irow)
IF (icol .EQ. irow) THEN
IF (ispin == 1) THEN
pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale + &
pmat(ikind)%mat(:, :, 2)*rscale
ELSE
pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale - &
pmat(ikind)%mat(:, :, 2)*rscale
END IF
nocc(ispin) = nocc(ispin) + nok(ispin, ikind)
END IF
END DO
CALL dbcsr_iterator_stop(iter)
CALL dbcsr_dot(matrix_p, matrix_s, trps1)
rds = 0.0_dp
! could be a ghost-atoms-only simulation
IF (nelectron_spin(ispin) > 0) THEN
rds = REAL(nelectron_spin(ispin), dp)/trps1
END IF
CALL dbcsr_scale(matrix_p, rds)
IF (ounit > 0) THEN
IF (nspin > 1) THEN
WRITE (UNIT=ounit, FMT="(T2,A,I1)") &
"Re-scaling the density matrix to get the right number of electrons for spin ", ispin
ELSE
WRITE (UNIT=ounit, FMT="(T2,A)") &
"Re-scaling the density matrix to get the right number of electrons"
END IF
WRITE (ounit, '(T19,A,T44,A,T67,A)') "# Electrons", "Trace(P)", "Scaling factor"
WRITE (ounit, '(T20,I10,T40,F12.3,T67,F14.3)') nelectron_spin(ispin), trps1, rds
END IF
IF (nspin > 1) THEN
group = para_env%group
CALL mp_sum(nocc, group)
IF (nelectron_spin(ispin) > nocc(ispin)) THEN
rds = 0.99_dp
CALL dbcsr_scale(matrix_p, rds)
rds = (1.0_dp - rds)*nelectron_spin(ispin)
CALL dbcsr_get_info(matrix_p, nfullcols_total=nc)
rds = rds/REAL(nc, KIND=dp)
CALL dbcsr_add_on_diag(matrix_p, rds)
IF (ounit > 0) THEN
WRITE (UNIT=ounit, FMT="(T4,A,/,T4,A,T59,F20.12)") &
"More MOs than initial guess orbitals detected", &
"Add constant to diagonal elements ", rds
END IF
END IF
END IF
END DO
DO ikind = 1, nkind
IF (ASSOCIATED(pmat(ikind)%mat)) THEN
DEALLOCATE (pmat(ikind)%mat)
END IF
END DO
DEALLOCATE (pmat)
DEALLOCATE (kind_of, nok)
CALL timestop(handle)
END SUBROUTINE calculate_atomic_block_dm
END MODULE qs_atomic_block

156
src/qs_initial_guess.F
View file

@ -38,11 +38,11 @@ MODULE qs_initial_guess
cp_print_key_unit_nr
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: &
dbcsr_add_on_diag, dbcsr_checksum, dbcsr_copy, dbcsr_dot, dbcsr_filter, dbcsr_get_diag, &
dbcsr_get_info, dbcsr_get_num_blocks, dbcsr_get_occupation, dbcsr_iterator_blocks_left, &
dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
dbcsr_multiply, dbcsr_nfullrows_total, dbcsr_p_type, dbcsr_release, dbcsr_scale, &
dbcsr_set, dbcsr_set_diag, dbcsr_type, dbcsr_verify_matrix
dbcsr_checksum, dbcsr_copy, dbcsr_dot, dbcsr_filter, dbcsr_get_diag, dbcsr_get_num_blocks, &
dbcsr_get_occupation, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
dbcsr_nfullrows_total, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, &
dbcsr_set_diag, dbcsr_type, dbcsr_verify_matrix
USE external_potential_types, ONLY: all_potential_type,&
gth_potential_type,&
sgp_potential_type
@ -66,6 +66,7 @@ MODULE qs_initial_guess
mp_sum
USE particle_methods, ONLY: get_particle_set
USE particle_types, ONLY: particle_type
USE qs_atomic_block, ONLY: calculate_atomic_block_dm
USE qs_density_matrices, ONLY: calculate_density_matrix
USE qs_dftb_utils, ONLY: get_dftb_atom_param
USE qs_environment_types, ONLY: get_qs_env,&
@ -107,7 +108,7 @@ MODULE qs_initial_guess
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_initial_guess'
PUBLIC :: calculate_first_density_matrix, calculate_atomic_block_dm, calculate_mopac_dm
PUBLIC :: calculate_first_density_matrix, calculate_mopac_dm
PUBLIC :: calculate_atomic_fock_matrix
TYPE atom_matrix_type
@ -1069,149 +1070,6 @@ CONTAINS
END SUBROUTINE calculate_first_density_matrix
! **************************************************************************************************
!> \brief returns a block diagonal density matrix. Blocks correspond to the atomic densities.
!> \param pmatrix ...
!> \param matrix_s ...
!> \param particle_set ...
!> \param atomic_kind_set ...
!> \param qs_kind_set ...
!> \param nspin ...
!> \param nelectron_spin ...
!> \param ounit ...
!> \param para_env ...
! **************************************************************************************************
SUBROUTINE calculate_atomic_block_dm(pmatrix, matrix_s, particle_set, atomic_kind_set, &
qs_kind_set, nspin, nelectron_spin, ounit, para_env)
TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: pmatrix
TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
INTEGER, INTENT(IN) :: nspin
INTEGER, DIMENSION(:), INTENT(IN) :: nelectron_spin
INTEGER, INTENT(IN) :: ounit
TYPE(cp_para_env_type) :: para_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atomic_block_dm'
INTEGER :: blk, group, handle, icol, ikind, irow, &
ispin, natom, nc, nkind, nocc(2)
INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nok
REAL(dp), DIMENSION(:, :), POINTER :: pdata
REAL(KIND=dp) :: rds, rscale, trps1
TYPE(atom_matrix_type), ALLOCATABLE, DIMENSION(:) :: pmat
TYPE(atomic_kind_type), POINTER :: atomic_kind
TYPE(dbcsr_iterator_type) :: iter
TYPE(dbcsr_type), POINTER :: matrix_p
TYPE(qs_kind_type), POINTER :: qs_kind
CALL timeset(routineN, handle)
natom = SIZE(particle_set)
nkind = SIZE(atomic_kind_set)
ALLOCATE (kind_of(natom))
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
ALLOCATE (pmat(nkind))
ALLOCATE (nok(2, nkind))
! precompute the atomic blocks corresponding to spherical atoms
DO ikind = 1, nkind
atomic_kind => atomic_kind_set(ikind)
qs_kind => qs_kind_set(ikind)
NULLIFY (pmat(ikind)%mat)
IF (ounit > 0) THEN
WRITE (UNIT=ounit, FMT="(/,T2,A)") &
"Guess for atomic kind: "//TRIM(atomic_kind%name)
END IF
CALL calculate_atomic_orbitals(atomic_kind, qs_kind, iunit=ounit, &
pmat=pmat(ikind)%mat, nocc=nocc)
nok(1:2, ikind) = nocc(1:2)
END DO
rscale = 1.0_dp
IF (nspin == 2) rscale = 0.5_dp
DO ispin = 1, nspin
IF ((ounit > 0) .AND. (nspin > 1)) THEN
WRITE (UNIT=ounit, FMT="(/,T2,A,I0)") "Spin ", ispin
END IF
matrix_p => pmatrix(ispin)%matrix
CALL dbcsr_set(matrix_p, 0.0_dp)
nocc(ispin) = 0
CALL dbcsr_iterator_start(iter, matrix_p)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, irow, icol, pdata, blk)
ikind = kind_of(irow)
IF (icol .EQ. irow) THEN
IF (ispin == 1) THEN
pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale + &
pmat(ikind)%mat(:, :, 2)*rscale
ELSE
pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale - &
pmat(ikind)%mat(:, :, 2)*rscale
END IF
nocc(ispin) = nocc(ispin) + nok(ispin, ikind)
END IF
END DO
CALL dbcsr_iterator_stop(iter)
CALL dbcsr_dot(matrix_p, matrix_s, trps1)
rds = 0.0_dp
! could be a ghost-atoms-only simulation
IF (nelectron_spin(ispin) > 0) THEN
rds = REAL(nelectron_spin(ispin), dp)/trps1
END IF
CALL dbcsr_scale(matrix_p, rds)
IF (ounit > 0) THEN
IF (nspin > 1) THEN
WRITE (UNIT=ounit, FMT="(T2,A,I1)") &
"Re-scaling the density matrix to get the right number of electrons for spin ", ispin
ELSE
WRITE (UNIT=ounit, FMT="(T2,A)") &
"Re-scaling the density matrix to get the right number of electrons"
END IF
WRITE (ounit, '(T19,A,T44,A,T67,A)') "# Electrons", "Trace(P)", "Scaling factor"
WRITE (ounit, '(T20,I10,T40,F12.3,T67,F14.3)') nelectron_spin(ispin), trps1, rds
END IF
IF (nspin > 1) THEN
group = para_env%group
CALL mp_sum(nocc, group)
IF (nelectron_spin(ispin) > nocc(ispin)) THEN
rds = 0.99_dp
CALL dbcsr_scale(matrix_p, rds)
rds = (1.0_dp - rds)*nelectron_spin(ispin)
CALL dbcsr_get_info(matrix_p, nfullcols_total=nc)
rds = rds/REAL(nc, KIND=dp)
CALL dbcsr_add_on_diag(matrix_p, rds)
IF (ounit > 0) THEN
WRITE (UNIT=ounit, FMT="(T4,A,/,T4,A,T59,F20.12)") &
"More MOs than initial guess orbitals detected", &
"Add constant to diagonal elements ", rds
END IF
END IF
END IF
END DO
DO ikind = 1, nkind
IF (ASSOCIATED(pmat(ikind)%mat)) THEN
DEALLOCATE (pmat(ikind)%mat)
END IF
END DO
DEALLOCATE (pmat)
DEALLOCATE (kind_of, nok)
CALL timestop(handle)
END SUBROUTINE calculate_atomic_block_dm
! **************************************************************************************************
!> \brief returns a block diagonal fock matrix.
!> \param matrix_f ...

266
src/qs_loc_methods.F
View file

@ -25,8 +25,11 @@ MODULE qs_loc_methods
USE cell_types, ONLY: cell_type,&
pbc
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
cp_dbcsr_sm_fm_multiply
cp_dbcsr_sm_fm_multiply,&
dbcsr_allocate_matrix_set,&
dbcsr_deallocate_matrix_set
USE cp_fm_basic_linalg, ONLY: cp_fm_schur_product
USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
fm_pool_create_fm
@ -34,8 +37,8 @@ MODULE qs_loc_methods
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: &
cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_p_type, &
cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type
cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_init_random, &
cp_fm_p_type, cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type
USE cp_gemm_interface, ONLY: cp_gemm
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type,&
@ -49,13 +52,17 @@ MODULE qs_loc_methods
USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
USE cp_units, ONLY: cp_unit_from_cp2k
USE dbcsr_api, ONLY: dbcsr_copy,&
dbcsr_create,&
dbcsr_deallocate_matrix,&
dbcsr_p_type,&
dbcsr_set
dbcsr_set,&
dbcsr_type,&
dbcsr_type_symmetric
USE input_constants, ONLY: &
do_loc_crazy, do_loc_direct, do_loc_jacobi, do_loc_l1_norm_sd, do_loc_none, do_loc_scdm, &
dump_dcd, dump_dcd_aligned_cell, dump_xmol, op_loc_berry, op_loc_boys, op_loc_pipek, &
state_loc_list
do_loc_cpo_atomic, do_loc_cpo_random, do_loc_cpo_restart, do_loc_cpo_space_nmo, &
do_loc_cpo_space_wan, do_loc_crazy, do_loc_direct, do_loc_gapo, do_loc_jacobi, &
do_loc_l1_norm_sd, do_loc_none, do_loc_scdm, dump_dcd, dump_dcd_aligned_cell, dump_xmol, &
op_loc_berry, op_loc_boys, op_loc_pipek, state_loc_list
USE input_section_types, ONLY: section_get_ival,&
section_get_ivals,&
section_get_lval,&
@ -89,6 +96,7 @@ MODULE qs_loc_methods
REALSPACE,&
RECIPROCALSPACE,&
pw_p_type
USE qs_atomic_block, ONLY: calculate_atomic_block_dm
USE qs_collocate_density, ONLY: calculate_wavefunction
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
@ -100,13 +108,17 @@ MODULE qs_loc_methods
USE qs_localization_methods, ONLY: approx_l1_norm_sd,&
crazy_rotations,&
direct_mini,&
jacobi_cg_edf_ls,&
jacobi_rotations,&
scdm_qrfact,&
zij_matrix
USE qs_matrix_pools, ONLY: mpools_get
USE qs_mo_methods, ONLY: make_basis_simple,&
make_basis_sm
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_p_type
USE qs_moments, ONLY: build_local_moment_matrix
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE qs_subsys_types, ONLY: qs_subsys_get,&
qs_subsys_type
USE string_utilities, ONLY: xstring
@ -144,6 +156,10 @@ CONTAINS
!> \param weights ...
!> \param ispin ...
!> \param print_loc_section ...
!> \param nextra ...
!> \param nmo ...
!> \param vectors_2 ...
!> \param guess_mos ...
!> \par History
!> 04.2005 created [MI]
!> \author MI
@ -153,7 +169,8 @@ CONTAINS
!> The file for the centers and the spreads have a xyz format
! **************************************************************************************************
SUBROUTINE optimize_loc_berry(method, qs_loc_env, vectors, op_sm_set, &
zij_fm_set, para_env, cell, weights, ispin, print_loc_section)
zij_fm_set, para_env, cell, weights, ispin, print_loc_section, &
nextra, nmo, vectors_2, guess_mos)
INTEGER, INTENT(IN) :: method
TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env
@ -165,20 +182,20 @@ CONTAINS
REAL(dp), DIMENSION(:) :: weights
INTEGER, INTENT(IN) :: ispin
TYPE(section_vals_type), POINTER :: print_loc_section
INTEGER, INTENT(IN), OPTIONAL :: nextra, nmo
TYPE(cp_fm_type), INTENT(IN), OPTIONAL, POINTER :: vectors_2, guess_mos
CHARACTER(len=*), PARAMETER :: routineN = 'optimize_loc_berry'
INTEGER :: handle, max_iter, nao, nmoloc, out_each, &
output_unit, sweeps
LOGICAL :: converged, crazy_use_diag, &
do_jacobi_refinement
do_jacobi_refinement, my_do_mixed
REAL(dp) :: crazy_scale, eps_localization, &
max_crazy_angle, start_time, &
target_time
TYPE(cp_logger_type), POINTER :: logger
! INTEGER :: i,j
CALL timeset(routineN, handle)
logger => cp_get_default_logger()
output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", &
@ -198,6 +215,7 @@ CONTAINS
target_time = qs_loc_env%target_time
start_time = qs_loc_env%start_time
do_jacobi_refinement = qs_loc_env%localized_wfn_control%jacobi_refinement
my_do_mixed = qs_loc_env%localized_wfn_control%do_mixed
CALL centers_spreads_berry(qs_loc_env, zij_fm_set, nmoloc, cell, weights, &
ispin, print_loc_section, only_initial_out=.TRUE.)
SELECT CASE (method)
@ -205,6 +223,28 @@ CONTAINS
CALL jacobi_rotations(weights, zij_fm_set, vectors, para_env, max_iter=max_iter, &
eps_localization=eps_localization, sweeps=sweeps, &
out_each=out_each, target_time=target_time, start_time=start_time)
CASE (do_loc_gapo)
IF (my_do_mixed) THEN
IF (nextra > 0) THEN
IF (PRESENT(guess_mos)) THEN
CALL jacobi_cg_edf_ls(para_env, weights, zij_fm_set, vectors, max_iter, &
eps_localization, sweeps, out_each, nextra, &
qs_loc_env%localized_wfn_control%do_cg_po, &
nmo=nmo, vectors_2=vectors_2, mos_guess=guess_mos)
ELSE
CALL jacobi_cg_edf_ls(para_env, weights, zij_fm_set, vectors, max_iter, &
eps_localization, sweeps, out_each, nextra, &
qs_loc_env%localized_wfn_control%do_cg_po, &
nmo=nmo, vectors_2=vectors_2)
END IF
ELSE
CALL jacobi_cg_edf_ls(para_env, weights, zij_fm_set, vectors, max_iter, &
eps_localization, sweeps, out_each, 0, &
qs_loc_env%localized_wfn_control%do_cg_po)
END IF
ELSE
CPABORT("GAPO works only with STATES MIXED")
END IF
CASE (do_loc_scdm)
! Decomposition
CALL scdm_qrfact(vectors)
@ -361,6 +401,8 @@ CONTAINS
SELECT CASE (method)
CASE (do_loc_jacobi)
CALL jacobi_rotation_pipek(zij_fm_set, vectors, sweeps)
CASE (do_loc_gapo)
CPABORT("GAPO and Pipek not implemented.")
CASE (do_loc_crazy)
CPABORT("Crazy and Pipek not implemented.")
CASE (do_loc_l1_norm_sd)
@ -418,8 +460,8 @@ CONTAINS
INTEGER :: idir, istate, jdir, nstates, &
output_unit, unit_out_s
LOGICAL :: my_only_init
REAL(dp) :: spread_i, spread_ii, sum_spread_i, &
sum_spread_ii
REAL(dp) :: avg_spread_ii, spread_i, spread_ii, &
sum_spread_i, sum_spread_ii
REAL(dp), DIMENSION(3) :: c, c2, cpbc
REAL(dp), DIMENSION(:, :), POINTER :: centers
REAL(KIND=dp) :: imagpart, realpart
@ -447,6 +489,7 @@ CONTAINS
CPASSERT(SIZE(centers, 2) == nmoloc)
sum_spread_i = 0.0_dp
sum_spread_ii = 0.0_dp
avg_spread_ii = 0.0_dp
DO istate = 1, nmoloc
c = 0.0_dp
c2 = 0.0_dp
@ -475,18 +518,21 @@ CONTAINS
sum_spread_ii = sum_spread_ii + centers(5, istate)
IF (unit_out_s > 0 .AND. .NOT. my_only_init) WRITE (unit_out_s, '(I6,2F16.8)') istate, centers(4:5, istate)
END DO
avg_spread_ii = sum_spread_ii/REAL(nmoloc, KIND=dp)
! Print of wannier centers
print_key => section_vals_get_subs_vals(print_loc_section, "WANNIER_CENTERS")
IF (.NOT. my_only_init) CALL print_wannier_centers(qs_loc_env, print_key, centers, logger, ispin)
IF (output_unit > 0) THEN
WRITE (output_unit, '(T4, A, 2x, A26, A26)') " Spread Functional ", "sum_in -w_i ln(|z_in|^2)", &
"sum_in w_i(1-|z_in|^2)"
WRITE (output_unit, '(T4, A, 2x, 2A26,/,T23, A28)') " Spread Functional ", "sum_in -w_i ln(|z_in|^2)", &
"sum_in w_i(1-|z_in|^2)", "sum_in w_i(1-|z_in|^2)/n"
IF (my_only_init) THEN
WRITE (output_unit, '(T4,A,T38,2F20.10)') " Initial Spread (Berry) : ", sum_spread_i, sum_spread_ii
WRITE (output_unit, '(T4,A,T38,2F20.10,/,T38,F20.10)') " Initial Spread (Berry) : ", &
sum_spread_i, sum_spread_ii, avg_spread_ii
ELSE
WRITE (output_unit, '(T4,A,T38,2F20.10)') " Total Spread (Berry) : ", sum_spread_i, sum_spread_ii
WRITE (output_unit, '(T4,A,T38,2F20.10,/,T38,F20.10)') " Total Spread (Berry) : ", &
sum_spread_i, sum_spread_ii, avg_spread_ii
END IF
END IF
@ -610,25 +656,35 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'qs_loc_driver'
CHARACTER(LEN=default_string_length) :: my_pos
INTEGER :: dim_op, handle, i, imo, imoloc, ir, &
ispin, istate, j, jstate, l_spin, lb, &
loc_method, n_rep, nao, ncubes, nmo, &
nmosub, s_spin, ub
INTEGER :: dim_op, handle, i, imo, imoloc, ir, ispin, istate, j, jstate, l_spin, lb, &
loc_method, n_rep, nao, ncubes, ndummy, nextra, ngextra, nguess, nmo, nmosub, norextra, &
s_spin, ub
INTEGER, DIMENSION(2) :: nelectron_spin
INTEGER, DIMENSION(:), POINTER :: bounds, list, list_cubes
LOGICAL :: append_cube, list_cubes_setup
LOGICAL :: append_cube, do_ortho, has_unit_metric, &
list_cubes_setup, my_guess_atomic, &
my_guess_wan
LOGICAL, SAVE :: first_time = .TRUE.
REAL(dp), DIMENSION(6) :: weights
REAL(KIND=dp), DIMENSION(:, :), POINTER :: centers, vecbuffer
REAL(KIND=dp), DIMENSION(:, :), POINTER :: centers, tmp_mat, vecbuffer
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: moloc_coeff
TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: op_fm_set
TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(cp_fm_type), POINTER :: mo_coeff, mos_guess, tmp_fm, tmp_fm_1, &
vectors_2
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, op_sm_set
TYPE(dbcsr_type), POINTER :: refmatrix, tmatrix
TYPE(dft_control_type), POINTER :: dft_control
TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(section_vals_type), POINTER :: print_key
CALL timeset(routineN, handle)
@ -652,6 +708,53 @@ CONTAINS
s_spin = myspin
l_spin = myspin
END IF
my_guess_atomic = .FALSE.
! SGh-wan: if atomic guess and do_mixed and nextra > 0
! read CPO_GUESS; CASE ATOMIC / RESTART / RANDOM (0/1/2)
! read CPO_GUESS_SPACE if CASE ATOMIC; CASE ALL / WAN
nextra = localized_wfn_control%nextra
IF (nextra > 0) THEN
my_guess_atomic = .TRUE.
my_guess_wan = .FALSE.
do_ortho = .TRUE.
SELECT CASE (localized_wfn_control%coeff_po_guess)
CASE (do_loc_cpo_atomic)
my_guess_atomic = .TRUE.
NULLIFY (atomic_kind_set, qs_kind_set, particle_set, matrix_s_kp, sab_orb, p_rmpv, &
refmatrix, tmatrix)
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set, &
matrix_s_kp=matrix_s_kp, &
has_unit_metric=has_unit_metric, &
nelectron_spin=nelectron_spin, &
sab_orb=sab_orb)
refmatrix => matrix_s_kp(1, 1)%matrix
! create p_rmpv
CALL dbcsr_allocate_matrix_set(p_rmpv, dft_control%nspins)
DO ispin = 1, dft_control%nspins
ALLOCATE (p_rmpv(ispin)%matrix)
tmatrix => p_rmpv(ispin)%matrix
CALL dbcsr_create(matrix=tmatrix, template=refmatrix, &
matrix_type=dbcsr_type_symmetric, nze=0)
CALL cp_dbcsr_alloc_block_from_nbl(tmatrix, sab_orb)
CALL dbcsr_set(tmatrix, 0.0_dp)
END DO
CALL calculate_atomic_block_dm(p_rmpv, refmatrix, &
particle_set, atomic_kind_set, qs_kind_set, &
dft_control%nspins, nelectron_spin, 0, para_env)
CASE (do_loc_cpo_restart)
my_guess_atomic = .FALSE.
my_guess_wan = .TRUE.
CASE (do_loc_cpo_random)
my_guess_atomic = .FALSE.
END SELECT
END IF
DO ispin = s_spin, l_spin
@ -685,8 +788,116 @@ CONTAINS
END DO
CALL cp_fm_struct_release(tmp_fm_struct)
CALL optimize_loc_berry(loc_method, qs_loc_env, moloc_coeff(ispin)%matrix, op_sm_set, &
op_fm_set, para_env, cell, weights, ispin, print_loc_section)
IF (localized_wfn_control%do_mixed) THEN
IF (nextra > 0) THEN
NULLIFY (vectors_2, mos_guess, tmp_fm, tmp_fm_1)
norextra = nmo - nmosub
CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff)
CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
ncol_global=norextra, para_env=para_env, context=mo_coeff%matrix_struct%context)
CALL cp_fm_create(vectors_2, tmp_fm_struct)
CALL cp_fm_struct_release(tmp_fm_struct)
ALLOCATE (tmp_mat(nao, norextra))
CALL cp_fm_get_submatrix(mo_coeff, tmp_mat, 1, nmosub + 1)
CALL cp_fm_set_submatrix(vectors_2, tmp_mat)
DEALLOCATE (tmp_mat)
! if guess "atomic" generate MOs based on atomic densities and
! pass on to optimize_loc_berry
IF (my_guess_atomic .OR. my_guess_wan) THEN
SELECT CASE (localized_wfn_control%coeff_po_guess_mo_space)
CASE (do_loc_cpo_space_wan)
ndummy = nmosub
CASE (do_loc_cpo_space_nmo)
ndummy = nmo
do_ortho = .FALSE.
END SELECT
CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
ncol_global=ndummy, para_env=para_env, &
context=mo_coeff%matrix_struct%context)
CALL cp_fm_create(mos_guess, tmp_fm_struct)
CALL cp_fm_set_all(mos_guess, 0.0_dp)
IF (my_guess_atomic) THEN
CALL cp_fm_create(tmp_fm, tmp_fm_struct)
CALL cp_fm_create(tmp_fm_1, tmp_fm_struct)
CALL cp_fm_set_all(tmp_fm, 0.0_dp)
CALL cp_fm_set_all(tmp_fm_1, 0.0_dp)
CALL cp_fm_init_random(tmp_fm, ndummy)
IF (has_unit_metric) THEN
CALL cp_fm_to_fm(tmp_fm, tmp_fm_1)
ELSE
! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc))
CALL cp_dbcsr_sm_fm_multiply(refmatrix, tmp_fm, tmp_fm_1, ndummy)
END IF
CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, tmp_fm_1, mos_guess, ndummy)
CALL cp_fm_release(tmp_fm)
CALL cp_fm_release(tmp_fm_1)
CALL cp_fm_struct_release(tmp_fm_struct)
ELSEIF (my_guess_wan) THEN
nguess = localized_wfn_control%nguess(ispin)
ALLOCATE (tmp_mat(nao, nguess))
CALL cp_fm_get_submatrix(moloc_coeff(ispin)%matrix, tmp_mat, 1, 1, nao, nguess)
CALL cp_fm_set_submatrix(mos_guess, tmp_mat, 1, 1, nao, nguess)
DEALLOCATE (tmp_mat)
ngextra = nmosub - nguess
!WRITE(*,*) 'nguess, ngextra = ', nguess, ngextra
CALL cp_fm_struct_release(tmp_fm_struct)
IF (ngextra > 0) THEN
CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
ncol_global=ngextra, para_env=para_env, &
context=mo_coeff%matrix_struct%context)
CALL cp_fm_create(tmp_fm, tmp_fm_struct)
CALL cp_fm_init_random(tmp_fm, ngextra)
ALLOCATE (tmp_mat(nao, ngextra))
CALL cp_fm_get_submatrix(tmp_fm, tmp_mat, 1, 1, nao, ngextra)
CALL cp_fm_set_submatrix(mos_guess, tmp_mat, 1, nguess + 1, nao, ngextra)
DEALLOCATE (tmp_mat)
CALL cp_fm_release(tmp_fm)
CALL cp_fm_struct_release(tmp_fm_struct)
ELSE
do_ortho = .FALSE.
END IF
ALLOCATE (tmp_mat(nao, nmosub))
CALL cp_fm_get_submatrix(mo_coeff, tmp_mat, 1, 1, nao, nmosub)
CALL cp_fm_set_submatrix(moloc_coeff(ispin)%matrix, tmp_mat)
DEALLOCATE (tmp_mat)
END IF
IF (do_ortho) THEN
IF ((my_guess_atomic) .OR. (my_guess_wan)) THEN
!! and ortho the result
IF (has_unit_metric) THEN
CALL make_basis_simple(mos_guess, ndummy)
ELSE
CALL make_basis_sm(mos_guess, ndummy, refmatrix)
END IF
END IF
END IF
CALL optimize_loc_berry(loc_method, qs_loc_env, moloc_coeff(ispin)%matrix, op_sm_set, &
op_fm_set, para_env, cell, weights, ispin, print_loc_section, &
nextra=nextra, nmo=nmo, vectors_2=vectors_2, &
guess_mos=mos_guess)
CALL cp_fm_release(mos_guess)
ELSE
CALL optimize_loc_berry(loc_method, qs_loc_env, moloc_coeff(ispin)%matrix, op_sm_set, &
op_fm_set, para_env, cell, weights, ispin, print_loc_section, &
nextra=nextra, nmo=nmo, vectors_2=vectors_2)
END IF
CALL cp_fm_release(vectors_2)
ELSE
CALL optimize_loc_berry(loc_method, qs_loc_env, moloc_coeff(ispin)%matrix, op_sm_set, &
op_fm_set, para_env, cell, weights, ispin, print_loc_section, nextra=0)
END IF
ELSE
CALL optimize_loc_berry(loc_method, qs_loc_env, moloc_coeff(ispin)%matrix, op_sm_set, &
op_fm_set, para_env, cell, weights, ispin, print_loc_section)
END IF
! Here we dealloctate op_fm_set
IF (ASSOCIATED(op_fm_set)) THEN
@ -825,6 +1036,7 @@ CONTAINS
DEALLOCATE (list_cubes)
END IF
END DO ! ispin
IF (my_guess_atomic) CALL dbcsr_deallocate_matrix_set(p_rmpv)
first_time = .FALSE.
CALL timestop(handle)
END SUBROUTINE qs_loc_driver

9
src/qs_loc_states.F
View file

@ -82,7 +82,7 @@ CONTAINS
INTEGER :: handle, ispin, mystate, ns, output_unit
INTEGER, DIMENSION(:), POINTER :: lstates, marked_states_spin
LOGICAL :: do_homo
LOGICAL :: do_homo, do_mixed
REAL(KIND=dp), DIMENSION(:, :), POINTER :: scenter
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, matrix_s
@ -101,6 +101,7 @@ CONTAINS
output_unit = cp_logger_get_default_io_unit(logger)
loc_print_section => section_vals_get_subs_vals(loc_section, "PRINT")
do_homo = qs_loc_env%localized_wfn_control%do_homo
do_mixed = qs_loc_env%localized_wfn_control%do_mixed
IF (BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, &
"WANNIER_STATES"), cp_p_file)) THEN
CALL get_qs_env(qs_env=qs_env, WannierCentres=wc)
@ -122,7 +123,11 @@ CONTAINS
IF (output_unit > 0 .AND. do_homo) WRITE (output_unit, "(/,T2,A,I3)") &
"LOCALIZATION| Computing localization properties "// &
"for OCCUPIED ORBITALS. Spin:", ispin
IF (output_unit > 0 .AND. (.NOT. do_homo)) WRITE (output_unit, "(/,T2,A,I3)") &
IF (output_unit > 0 .AND. do_mixed) WRITE (output_unit, "(/,T2,A,/,T16,A,I3)") &
"LOCALIZATION| Computing localization properties for OCCUPIED, ", &
"PARTIALLY OCCUPIED and UNOCCUPIED ORBITALS. Spin:", ispin
IF (output_unit > 0 .AND. (.NOT. do_homo) .AND. (.NOT. do_mixed)) &
WRITE (output_unit, "(/,T2,A,I3)") &
"LOCALIZATION| Computing localization properties "// &
"for UNOCCUPIED ORBITALS. Spin:", ispin

7
src/qs_loc_types.F
View file

@ -127,11 +127,13 @@ MODULE qs_loc_types
INTEGER :: min_or_max
INTEGER :: localization_method
INTEGER :: operator_type
INTEGER, DIMENSION(2) :: nloc_states
INTEGER, DIMENSION(2) :: nloc_states, nguess
INTEGER :: set_of_states
INTEGER, DIMENSION(2, 2) :: lu_bound_states
INTEGER :: max_iter
INTEGER :: out_each
INTEGER :: nextra
INTEGER :: coeff_po_guess, coeff_po_guess_mo_space
REAL(KIND=dp) :: eps_localization
REAL(KIND=dp) :: max_crazy_angle
REAL(KIND=dp) :: crazy_scale
@ -142,6 +144,7 @@ MODULE qs_loc_types
LOGICAL :: print_centers
LOGICAL :: print_spreads
LOGICAL :: do_homo
LOGICAL :: do_mixed, do_cg_po
LOGICAL :: loc_restart
LOGICAL :: use_history
INTEGER, POINTER, DIMENSION(:, :) :: loc_states
@ -298,6 +301,8 @@ CONTAINS
localized_wfn_control%ref_count = 1
localized_wfn_control%nloc_states = 0
localized_wfn_control%nextra = 0
localized_wfn_control%nguess = 0
localized_wfn_control%lu_bound_states = 0
localized_wfn_control%lu_ene_bound = 0.0_dp
localized_wfn_control%print_cubes = .FALSE.

196
src/qs_loc_utils.F
View file

@ -50,9 +50,9 @@ MODULE qs_loc_utils
dbcsr_type
USE distribution_1d_types, ONLY: distribution_1d_type
USE input_constants, ONLY: &
do_loc_crazy, do_loc_direct, do_loc_jacobi, do_loc_l1_norm_sd, do_loc_none, do_loc_scdm, &
energy_loc_range, op_loc_berry, op_loc_boys, op_loc_pipek, state_loc_all, state_loc_list, &
state_loc_none, state_loc_range
do_loc_crazy, do_loc_direct, do_loc_gapo, do_loc_jacobi, do_loc_l1_norm_sd, do_loc_none, &
do_loc_scdm, energy_loc_range, op_loc_berry, op_loc_boys, op_loc_pipek, state_loc_all, &
state_loc_list, state_loc_mixed, state_loc_none, state_loc_range
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
@ -500,7 +500,8 @@ CONTAINS
CALL cp_fm_get_info(moloc_coeff(ispin)%matrix, nrow_global=naosub, &
ncol_global=nmosub)
CPASSERT(nao == naosub)
IF (localized_wfn_control%do_homo) THEN
IF ((localized_wfn_control%do_homo) .OR. &
(localized_wfn_control%set_of_states == state_loc_mixed)) THEN
CPASSERT(nmo >= nmosub)
ELSE
CPASSERT(nao - nmo >= nmosub)
@ -521,7 +522,8 @@ CONTAINS
ELSE
mat_ptr => mo_coeff
END IF
IF (localized_wfn_control%set_of_states == state_loc_list) THEN
IF ((localized_wfn_control%set_of_states == state_loc_list) .OR. &
(localized_wfn_control%set_of_states == state_loc_mixed)) THEN
ALLOCATE (vecbuffer(1, nao))
IF (localized_wfn_control%do_homo) THEN
my_occ = occupations(localized_wfn_control%loc_states(1, ispin))
@ -906,9 +908,10 @@ CONTAINS
!> \param coeff_localized ...
!> \param do_homo ...
!> \param evals ...
!> \param do_mixed ...
! **************************************************************************************************
SUBROUTINE loc_write_restart(qs_loc_env, section, mo_array, coeff_localized, &
do_homo, evals)
do_homo, evals, do_mixed)
TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env
TYPE(section_vals_type), POINTER :: section
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mo_array
@ -916,6 +919,7 @@ CONTAINS
LOGICAL, INTENT(IN) :: do_homo
TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: evals
LOGICAL, INTENT(IN), OPTIONAL :: do_mixed
CHARACTER(LEN=*), PARAMETER :: routineN = 'loc_write_restart'
@ -923,6 +927,7 @@ CONTAINS
CHARACTER(LEN=default_string_length) :: my_middle
INTEGER :: handle, ispin, max_block, nao, nloc, &
nmo, output_unit, rst_unit
LOGICAL :: my_do_mixed
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: print_key
@ -942,8 +947,12 @@ CONTAINS
! Open file
rst_unit = -1
my_do_mixed = .FALSE.
IF (PRESENT(do_mixed)) my_do_mixed = do_mixed
IF (do_homo) THEN
my_middle = "LOC_HOMO"
ELSEIF (my_do_mixed) THEN
my_middle = "LOC_MIXED"
ELSE
my_middle = "LOC_LUMO"
END IF
@ -974,7 +983,7 @@ CONTAINS
nloc = qs_loc_env%localized_wfn_control%nloc_states(ispin)
IF (rst_unit > 0) THEN
WRITE (rst_unit) qs_loc_env%localized_wfn_control%loc_states(1:nloc, ispin)
IF (do_homo) THEN
IF (do_homo .OR. my_do_mixed) THEN
WRITE (rst_unit) nmo, &
mo_array(ispin)%mo_set%homo, &
mo_array(ispin)%mo_set%lfomo, &
@ -1012,9 +1021,10 @@ CONTAINS
!> \param do_homo ...
!> \param restart_found ...
!> \param evals ...
!> \param do_mixed ...
! **************************************************************************************************
SUBROUTINE loc_read_restart(qs_loc_env, mos, mos_localized, section, section2, para_env, &
do_homo, restart_found, evals)
do_homo, restart_found, evals, do_mixed)
TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos
@ -1025,6 +1035,7 @@ CONTAINS
LOGICAL, INTENT(INOUT) :: restart_found
TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: evals
LOGICAL, INTENT(IN), OPTIONAL :: do_mixed
CHARACTER(len=*), PARAMETER :: routineN = 'loc_read_restart'
@ -1033,7 +1044,7 @@ CONTAINS
CHARACTER(LEN=default_string_length) :: my_middle
INTEGER :: group, handle, homo_read, i, ispin, lfomo_read, max_nloc, n_rep_val, nao, &
nelectron_read, nloc, nmo, nmo_read, nspin, output_unit, rst_unit, source
LOGICAL :: file_exists
LOGICAL :: file_exists, my_do_mixed
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig_read, occ_read
REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
TYPE(cp_logger_type), POINTER :: logger
@ -1052,8 +1063,12 @@ CONTAINS
output_unit = cp_print_key_unit_nr(logger, section2, &
"PROGRAM_RUN_INFO", extension=".Log")
my_do_mixed = .FALSE.
IF (PRESENT(do_mixed)) my_do_mixed = do_mixed
IF (do_homo) THEN
fname_key = "LOCHOMO_RESTART_FILE_NAME"
ELSEIF (my_do_mixed) THEN
fname_key = "LOCMIXD_RESTART_FILE_NAME"
ELSE
fname_key = "LOCLUMO_RESTART_FILE_NAME"
IF (.NOT. PRESENT(evals)) &
@ -1069,6 +1084,8 @@ CONTAINS
print_key => section_vals_get_subs_vals(section2, "LOC_RESTART")
IF (do_homo) THEN
my_middle = "LOC_HOMO"
ELSEIF (my_do_mixed) THEN
my_middle = "LOC_MIXED"
ELSE
my_middle = "LOC_LUMO"
END IF
@ -1114,7 +1131,7 @@ CONTAINS
qs_loc_env%localized_wfn_control%loc_states = 0
DO ispin = 1, nspin
IF (do_homo) THEN
IF (do_homo .OR. do_mixed) THEN
nmo = mos(ispin)%mo_set%nmo
ELSE
nmo = SIZE(evals(ispin)%array, 1)
@ -1122,7 +1139,7 @@ CONTAINS
IF (para_env%ionode .AND. (nmo > 0)) THEN
nloc = qs_loc_env%localized_wfn_control%nloc_states(ispin)
READ (rst_unit) qs_loc_env%localized_wfn_control%loc_states(1:nloc, ispin)
IF (do_homo) THEN
IF (do_homo .OR. do_mixed) THEN
READ (rst_unit) nmo_read, homo_read, lfomo_read, nelectron_read
ALLOCATE (eig_read(nmo_read), occ_read(nmo_read))
eig_read = 0.0_dp
@ -1144,7 +1161,7 @@ CONTAINS
"the allocated MOs. The read MO set will be truncated!")
nmo = MIN(nmo, nmo_read)
IF (do_homo) THEN
IF (do_homo .OR. do_mixed) THEN
mos(ispin)%mo_set%eigenvalues(1:nmo) = eig_read(1:nmo)
mos(ispin)%mo_set%occupation_numbers(1:nmo) = occ_read(1:nmo)
DEALLOCATE (eig_read, occ_read)
@ -1154,7 +1171,7 @@ CONTAINS
END IF
END IF
IF (do_homo) THEN
IF (do_homo .OR. do_mixed) THEN
CALL mp_bcast(mos(ispin)%mo_set%eigenvalues, source, group)
CALL mp_bcast(mos(ispin)%mo_set%occupation_numbers, source, group)
ELSE
@ -1193,31 +1210,43 @@ CONTAINS
!> \param qs_loc_env ...
!> \param loc_section ...
!> \param do_homo ...
!> \param do_mixed ...
!> \param do_xas ...
!> \param nloc_xas ...
!> \param spin_xas ...
!> \par History
!> 2009 created
! **************************************************************************************************
SUBROUTINE qs_loc_control_init(qs_loc_env, loc_section, do_homo, do_xas, nloc_xas, spin_xas)
SUBROUTINE qs_loc_control_init(qs_loc_env, loc_section, do_homo, do_mixed, &
do_xas, nloc_xas, spin_xas)
TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env
TYPE(section_vals_type), POINTER :: loc_section
LOGICAL, INTENT(IN) :: do_homo
LOGICAL, INTENT(IN), OPTIONAL :: do_xas
LOGICAL, INTENT(IN), OPTIONAL :: do_mixed, do_xas
INTEGER, INTENT(IN), OPTIONAL :: nloc_xas, spin_xas
CHARACTER(len=*), PARAMETER :: routineN = 'qs_loc_control_init', &
routineP = moduleN//':'//routineN
LOGICAL :: my_do_mixed
TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
NULLIFY (localized_wfn_control)
IF (PRESENT(do_mixed)) THEN
my_do_mixed = do_mixed
ELSE
my_do_mixed = .FALSE.
END IF
CALL localized_wfn_control_create(localized_wfn_control)
CALL set_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
CALL localized_wfn_control_release(localized_wfn_control)
CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
localized_wfn_control%do_homo = do_homo
CALL read_loc_section(localized_wfn_control, loc_section, &
qs_loc_env%do_localize, do_xas, nloc_xas, spin_xas)
localized_wfn_control%do_mixed = my_do_mixed
CALL read_loc_section(localized_wfn_control, loc_section, qs_loc_env%do_localize, &
my_do_mixed, do_xas, nloc_xas, spin_xas)
END SUBROUTINE qs_loc_control_init
@ -1231,9 +1260,12 @@ CONTAINS
!> \param do_mo_cubes ...
!> \param mo_loc_history ...
!> \param evals ...
!> \param tot_zeff_corr ...
!> \param do_mixed ...
! **************************************************************************************************
SUBROUTINE qs_loc_init(qs_env, qs_loc_env, localize_section, mos_localized, &
do_homo, do_mo_cubes, mo_loc_history, evals)
do_homo, do_mo_cubes, mo_loc_history, evals, &
tot_zeff_corr, do_mixed)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env
TYPE(section_vals_type), POINTER :: localize_section
@ -1243,14 +1275,16 @@ CONTAINS
POINTER :: mo_loc_history
TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: evals
REAL(KIND=dp), INTENT(IN), OPTIONAL :: tot_zeff_corr
LOGICAL, OPTIONAL :: do_mixed
CHARACTER(len=*), PARAMETER :: routineN = 'qs_loc_init'
INTEGER :: handle, homo, i, ilast_intocc, ilow, &
ispin, iup, n_mo(2), n_mos(2), nao, &
nelectron, nmoloc(2), nspin, &
output_unit
LOGICAL :: my_do_homo, my_do_mo_cubes, restart_found
INTEGER :: handle, homo, i, ilast_intocc, ilow, ispin, iup, n_mo(2), n_mos(2), nao, &
nelectron, nextra, nmoloc(2), nocc, npocc, nspin, output_unit
LOGICAL :: my_do_homo, my_do_mixed, my_do_mo_cubes, &
restart_found
REAL(KIND=dp) :: maxocc, my_tot_zeff_corr
REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, occupation
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(cp_logger_type), POINTER :: logger
@ -1287,22 +1321,65 @@ CONTAINS
ELSE
my_do_mo_cubes = .FALSE.
END IF
IF (PRESENT(do_mixed)) THEN
my_do_mixed = do_mixed
ELSE
my_do_mixed = .FALSE.
END IF
IF (PRESENT(tot_zeff_corr)) THEN
my_tot_zeff_corr = tot_zeff_corr
ELSE
my_tot_zeff_corr = 0.0_dp
END IF
restart_found = .FALSE.
IF (qs_loc_env%do_localize) THEN
! Some setup for MOs to be localized
CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
IF (localized_wfn_control%loc_restart) THEN
IF (localized_wfn_control%nextra > 0) THEN
! currently only the occupied guess is read
my_do_homo = .FALSE.
END IF
CALL loc_read_restart(qs_loc_env, mos, mos_localized, localize_section, &
loc_print_section, para_env, my_do_homo, restart_found, evals=evals)
loc_print_section, para_env, my_do_homo, restart_found, evals=evals, &
do_mixed=my_do_mixed)
IF (output_unit > 0) WRITE (output_unit, "(/,T2,A,A)") "LOCALIZATION| ", &
" The orbitals to be localized are read from localization restart file."
nmoloc = localized_wfn_control%nloc_states
localized_wfn_control%nguess = nmoloc
IF (localized_wfn_control%nextra > 0) THEN
! reset different variables in localized_wfn_control:
! lu_bound_states, nloc_states, loc_states
localized_wfn_control%loc_restart = restart_found
localized_wfn_control%set_of_states = state_loc_mixed
DO ispin = 1, nspin
CALL get_mo_set(mos(ispin)%mo_set, homo=homo, occupation_numbers=occupation, &
maxocc=maxocc)
nextra = localized_wfn_control%nextra
nocc = homo
DO i = nocc, 1, -1
IF (maxocc - occupation(i) < localized_wfn_control%eps_occ) THEN
ilast_intocc = i
EXIT
END IF
END DO
nocc = ilast_intocc
npocc = homo - nocc
nmoloc(ispin) = nocc + nextra
localized_wfn_control%lu_bound_states(1, ispin) = 1
localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
localized_wfn_control%nloc_states(ispin) = nmoloc(ispin)
END DO
my_do_homo = .FALSE.
END IF
END IF
IF (.NOT. restart_found) THEN
nmoloc = 0
DO ispin = 1, nspin
CALL get_mo_set(mos(ispin)%mo_set, nmo=n_mo(ispin), nelectron=nelectron, homo=homo, nao=nao, &
mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues, occupation_numbers=occupation)
mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues, occupation_numbers=occupation, &
maxocc=maxocc)
! Get eigenstates (only needed if not already calculated before)
IF ((.NOT. my_do_mo_cubes) &
! .OR. section_get_ival(dft_section,"PRINT%MO_CUBES%NHOMO")==0)&
@ -1366,6 +1443,28 @@ CONTAINS
"LOCALIZATION| Spin ", ispin, " The first ", &
nmoloc(ispin), " virtual orbitals are localized,", " with energies from ", &
mo_eigenvalues(homo + 1), " to ", mo_eigenvalues(n_mo(ispin)), " [a.u.]."
ELSE IF (localized_wfn_control%set_of_states == state_loc_mixed) THEN
nextra = localized_wfn_control%nextra
nocc = homo
DO i = nocc, 1, -1
IF (maxocc - occupation(i) < localized_wfn_control%eps_occ) THEN
ilast_intocc = i
EXIT
END IF
END DO
nocc = ilast_intocc
npocc = homo - nocc
nmoloc(ispin) = nocc + nextra
localized_wfn_control%lu_bound_states(1, ispin) = 1
localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
IF (output_unit > 0) &
WRITE (output_unit, "(/,T2,A,I4,A,I6,A,/,T15,A,I6,/,T15,A,I6,/,T15,A,I6,/,T15,A,F12.6,A)") &
"LOCALIZATION| Spin ", ispin, " The first ", &
nmoloc(ispin), " orbitals are localized.", &
"Number of fully occupied MOs: ", nocc, &
"Number of partially occupied MOs: ", npocc, &
"Number of extra degrees of freedom: ", nextra, &
"Excess charge: ", my_tot_zeff_corr, " electrons"
ELSE
nmoloc(ispin) = MIN(localized_wfn_control%nloc_states(1), n_mo(ispin))
IF (output_unit > 0 .AND. my_do_homo) WRITE (output_unit, "(/,T2,A,I4,A,I6,A)") "LOCALIZATION| Spin ", ispin, &
@ -1389,11 +1488,11 @@ CONTAINS
END IF
END DO ! ispin
n_mos(:) = nao - n_mo(:)
IF (my_do_homo) n_mos = n_mo
IF (my_do_homo .OR. my_do_mixed) n_mos = n_mo
CALL set_loc_wfn_lists(localized_wfn_control, nmoloc, n_mos, nspin)
END IF
CALL set_loc_centers(localized_wfn_control, nmoloc, nspin)
IF (my_do_homo) THEN
IF (my_do_homo .OR. my_do_mixed) THEN
CALL qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, &
loc_coeff=mos_localized, mo_loc_history=mo_loc_history)
END IF
@ -1413,6 +1512,7 @@ CONTAINS
!> \param localized_wfn_control ...
!> \param loc_section ...
!> \param localize ...
!> \param do_mixed ...
!> \param do_xas ...
!> \param nloc_xas ...
!> \param spin_channel_xas ...
@ -1420,19 +1520,19 @@ CONTAINS
!> 05.2005 created [MI]
! **************************************************************************************************
SUBROUTINE read_loc_section(localized_wfn_control, loc_section, &
localize, do_xas, nloc_xas, spin_channel_xas)
localize, do_mixed, do_xas, nloc_xas, spin_channel_xas)
TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
TYPE(section_vals_type), POINTER :: loc_section
LOGICAL, INTENT(OUT) :: localize
LOGICAL, INTENT(IN), OPTIONAL :: do_xas
LOGICAL, INTENT(IN), OPTIONAL :: do_mixed, do_xas
INTEGER, INTENT(IN), OPTIONAL :: nloc_xas, spin_channel_xas
INTEGER :: i, ind, ir, n_list, n_rep, n_state, &
nline, other_spin, output_unit, &
spin_xas
nextra, nline, other_spin, &
output_unit, spin_xas
INTEGER, DIMENSION(:), POINTER :: list, loc_list
LOGICAL :: my_do_xas
LOGICAL :: my_do_mixed, my_do_xas
REAL(dp), POINTER :: ene(:)
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: loc_print_section
@ -1444,6 +1544,10 @@ CONTAINS
CPASSERT(PRESENT(nloc_xas))
END IF
IF (PRESENT(spin_channel_xas)) spin_xas = spin_channel_xas
my_do_mixed = .FALSE.
IF (PRESENT(do_mixed)) THEN
my_do_mixed = do_mixed
END IF
CPASSERT(ASSOCIATED(loc_section))
NULLIFY (logger)
logger => cp_get_default_logger()
@ -1457,6 +1561,7 @@ CONTAINS
localized_wfn_control%lu_ene_bound = 0.0_dp
localized_wfn_control%nloc_states = 0
localized_wfn_control%set_of_states = 0
localized_wfn_control%nextra = 0
n_state = 0
CALL section_vals_val_get(loc_section, "MAX_ITER", &
@ -1487,6 +1592,14 @@ CONTAINS
l_val=localized_wfn_control%loc_restart)
CALL section_vals_val_get(loc_section, "USE_HISTORY", &
l_val=localized_wfn_control%use_history)
CALL section_vals_val_get(loc_section, "NEXTRA", &
i_val=localized_wfn_control%nextra)
CALL section_vals_val_get(loc_section, "CPO_GUESS", &
i_val=localized_wfn_control%coeff_po_guess)
CALL section_vals_val_get(loc_section, "CPO_GUESS_SPACE", &
i_val=localized_wfn_control%coeff_po_guess_mo_space)
CALL section_vals_val_get(loc_section, "CG_PO", &
l_val=localized_wfn_control%do_cg_po)
IF (localized_wfn_control%do_homo) THEN
! List of States HOMO
@ -1569,6 +1682,9 @@ CONTAINS
localized_wfn_control%nloc_states(spin_xas) = nloc_xas
localized_wfn_control%lu_bound_states(1, spin_xas) = 1
localized_wfn_control%lu_bound_states(2, spin_xas) = nloc_xas
ELSE IF (my_do_mixed) THEN
localized_wfn_control%set_of_states = state_loc_mixed
nextra = localized_wfn_control%nextra
ELSE
localized_wfn_control%set_of_states = state_loc_all
END IF
@ -1622,6 +1738,9 @@ CONTAINS
WRITE (UNIT=output_unit, FMT="(T2,A,T65,/,f16.6,A,f16.6,A)") &
"LOCALIZE| Orbitals to be localized: Those with energy in the range between ", &
localized_wfn_control%lu_ene_bound(1), " and ", localized_wfn_control%lu_ene_bound(2), " a.u."
CASE (state_loc_mixed)
WRITE (UNIT=output_unit, FMT="(T2,A,I4,A)") &
"LOCALIZE| Orbitals to be localized: Occupied orbitals + ", nextra, " orbitals"
CASE DEFAULT
WRITE (UNIT=output_unit, FMT="(T2,A)") &
"LOCALIZE| Orbitals to be localized: None "
@ -1652,6 +1771,10 @@ CONTAINS
"LOCALIZE| scaling: ", localized_wfn_control%crazy_scale
WRITE (UNIT=output_unit, FMT="(T2,A,L1)") &
"LOCALIZE| use diag:", localized_wfn_control%crazy_use_diag
CASE (do_loc_gapo)
WRITE (UNIT=output_unit, FMT="(T2,A)") &
"LOCALIZE| Optimal unitary transformation generated by gradient ascent algorithm "// &
" for partially occupied wannier functions"
CASE (do_loc_direct)
WRITE (UNIT=output_unit, FMT="(T2,A)") &
"LOCALIZE| Optimal unitary transformation generated by direct algorithm"
@ -1790,6 +1913,15 @@ CONTAINS
END DO
END DO
END IF
CASE (state_loc_mixed)
! Mixed
ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
localized_wfn_control%loc_states = 0
DO ispin = 1, nspins
DO i = 1, nmoloc(ispin)
localized_wfn_control%loc_states(i, ispin) = i
END DO
END DO
END SELECT
CALL timestop(state)

1303
src/qs_localization_methods.F
View file

File diff suppressed because it is too large Load diff

98
src/qs_scf_post_gpw.F
View file

@ -66,14 +66,9 @@ MODULE qs_scf_post_gpw
hirshfeld_type,&
release_hirshfeld_type,&
set_hirshfeld_info
USE input_constants, ONLY: do_loc_both,&
do_loc_homo,&
do_loc_lumo,&
ot_precond_full_all,&
radius_covalent,&
radius_user,&
ref_charge_atomic,&
ref_charge_mulliken
USE input_constants, ONLY: &
do_loc_both, do_loc_homo, do_loc_lumo, do_loc_mixed, ot_precond_full_all, radius_covalent, &
radius_user, ref_charge_atomic, ref_charge_mulliken
USE input_section_types, ONLY: section_get_ival,&
section_get_ivals,&
section_get_lval,&
@ -249,22 +244,22 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'scf_post_calculation_gpw'
INTEGER :: handle, homo, ispin, min_lumos, n_rep, &
nhomo, nlumo, nlumo_stm, nlumo_tddft, &
nlumos, nmo, output_unit, unit_nr
INTEGER :: handle, homo, ispin, min_lumos, n_rep, nchk_nmoloc, nhomo, nlumo, nlumo_stm, &
nlumo_tddft, nlumos, nmo, output_unit, unit_nr
INTEGER, DIMENSION(:, :, :), POINTER :: marked_states
LOGICAL :: check_write, compute_lumos, do_homo, do_kpoints, do_mo_cubes, do_stm, &
LOGICAL :: check_write, compute_lumos, do_homo, do_kpoints, do_mixed, do_mo_cubes, do_stm, &
do_wannier_cubes, has_homo, has_lumo, loc_explicit, loc_print_explicit, my_do_mp2, &
my_localized_wfn, p_loc, p_loc_homo, p_loc_lumo
my_localized_wfn, p_loc, p_loc_homo, p_loc_lumo, p_loc_mixed
REAL(dp) :: e_kin
REAL(KIND=dp) :: gap, homo_lumo(2, 2)
REAL(KIND=dp) :: gap, homo_lumo(2, 2), total_zeff_corr
REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
TYPE(admm_type), POINTER :: admm_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: occupied_evals, unoccupied_evals, &
unoccupied_evals_stm
TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: mixed_evals, occupied_evals, &
unoccupied_evals, unoccupied_evals_stm
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: homo_localized, lumo_localized, lumo_ptr, &
mo_loc_history, occupied_orbs, unoccupied_orbs, unoccupied_orbs_stm
mixed_localized, mixed_orbs, mo_loc_history, occupied_orbs, unoccupied_orbs, &
unoccupied_orbs_stm
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, matrix_p_mp2, matrix_s, &
@ -280,7 +275,8 @@ CONTAINS
TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env_homo, qs_loc_env_lumo
TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env_homo, qs_loc_env_lumo, &
qs_loc_env_mixed
TYPE(qs_rho_type), POINTER :: rho
TYPE(qs_scf_env_type), POINTER :: scf_env
TYPE(qs_subsys_type), POINTER :: subsys
@ -314,14 +310,16 @@ CONTAINS
mo_coeff, ks_rmpv, matrix_s, qs_loc_env_homo, qs_loc_env_lumo, scf_control, &
unoccupied_orbs, unoccupied_orbs_stm, mo_eigenvalues, unoccupied_evals, &
unoccupied_evals_stm, molecule_set, mo_derivs, &
subsys, particles, input, print_key, kinetic_m, marked_states)
NULLIFY (homo_localized, lumo_localized, lumo_ptr, rho_ao)
subsys, particles, input, print_key, kinetic_m, marked_states, &
mixed_orbs, mixed_evals, qs_loc_env_mixed)
NULLIFY (homo_localized, lumo_localized, lumo_ptr, rho_ao, mixed_localized)
has_homo = .FALSE.
has_lumo = .FALSE.
p_loc = .FALSE.
p_loc_homo = .FALSE.
p_loc_lumo = .FALSE.
p_loc_mixed = .FALSE.
CPASSERT(ASSOCIATED(scf_env))
CPASSERT(scf_env%ref_count > 0)
@ -408,10 +406,12 @@ CONTAINS
section_get_ival(localize_section, "STATES") == do_loc_both) .AND. p_loc
p_loc_lumo = (section_get_ival(localize_section, "STATES") == do_loc_lumo .OR. &
section_get_ival(localize_section, "STATES") == do_loc_both) .AND. p_loc
p_loc_mixed = (section_get_ival(localize_section, "STATES") == do_loc_mixed) .AND. p_loc
CALL section_vals_val_get(localize_section, "LIST_UNOCCUPIED", n_rep_val=n_rep)
ELSE
p_loc_homo = .FALSE.
p_loc_lumo = .FALSE.
p_loc_mixed = .FALSE.
n_rep = 0
END IF
@ -660,6 +660,63 @@ CONTAINS
END IF
END IF
IF (p_loc_mixed) THEN
IF (do_kpoints) THEN
CPWARN("Localization not implemented for k-point calculations!!")
ELSEIF (dft_control%restricted) THEN
IF (output_unit > 0) WRITE (output_unit, *) &
" Unclear how we define MOs / localization in the restricted case... skipping"
ELSE
ALLOCATE (mixed_orbs(dft_control%nspins))
ALLOCATE (mixed_evals(dft_control%nspins))
ALLOCATE (mixed_localized(dft_control%nspins))
DO ispin = 1, dft_control%nspins
CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, &
eigenvalues=mo_eigenvalues)
mixed_orbs(ispin)%matrix => mo_coeff
mixed_evals(ispin)%array => mo_eigenvalues
CALL cp_fm_create(mixed_localized(ispin)%matrix, mixed_orbs(ispin)%matrix%matrix_struct)
CALL cp_fm_to_fm(mixed_orbs(ispin)%matrix, mixed_localized(ispin)%matrix)
END DO
CALL get_qs_env(qs_env, mo_loc_history=mo_loc_history)
do_homo = .FALSE.
do_mixed = .TRUE.
total_zeff_corr = scf_env%sum_zeff_corr
CALL qs_loc_env_create(qs_loc_env_mixed)
CALL qs_loc_control_init(qs_loc_env_mixed, localize_section, do_homo=do_homo, do_mixed=do_mixed)
CALL qs_loc_init(qs_env, qs_loc_env_mixed, localize_section, mixed_localized, do_homo, &
do_mo_cubes, mo_loc_history=mo_loc_history, tot_zeff_corr=total_zeff_corr, &
do_mixed=do_mixed)
DO ispin = 1, dft_control%nspins
CALL cp_fm_get_info(mixed_localized(ispin)%matrix, ncol_global=nchk_nmoloc)
END DO
CALL get_localization_info(qs_env, qs_loc_env_mixed, localize_section, mixed_localized, &
wf_r, wf_g, particles, mixed_orbs, mixed_evals, marked_states)
!retain the homo_localized for future use
IF (qs_loc_env_mixed%localized_wfn_control%use_history) THEN
CALL retain_history(mo_loc_history, mixed_localized)
CALL set_qs_env(qs_env, mo_loc_history=mo_loc_history)
END IF
!write restart for localization of occupied orbitals
CALL loc_write_restart(qs_loc_env_mixed, loc_print_section, mos, &
mixed_localized, do_homo, do_mixed=do_mixed)
CALL cp_fm_vect_dealloc(mixed_localized)
DEALLOCATE (mixed_orbs)
DEALLOCATE (mixed_evals)
! Print Total Dipole if the localization has been performed
! Revisit the formalism later
!IF (qs_loc_env_mixed%do_localize) THEN
! CALL loc_dipole(input, dft_control, qs_loc_env_mixed, logger, qs_env)
!END IF
END IF
END IF
! Deallocate grids needed to compute wavefunctions
IF (((do_mo_cubes .OR. do_wannier_cubes) .AND. (nlumo /= 0 .OR. nhomo /= 0)) .OR. p_loc) THEN
CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_r%pw)
@ -670,6 +727,7 @@ CONTAINS
IF (.NOT. do_kpoints) THEN
IF (p_loc_homo) CALL qs_loc_env_destroy(qs_loc_env_homo)
IF (p_loc_lumo) CALL qs_loc_env_destroy(qs_loc_env_lumo)
IF (p_loc_mixed) CALL qs_loc_env_destroy(qs_loc_env_mixed)
END IF
! generate a mix of wfns, and write to a restart

6
tests/QS/regtest-loc_powf/TEST_FILES Normal file
View file

@ -0,0 +1,6 @@
# runs are executed in the same order as in this file
# the second field tells which option from cp2k/tests/TEST_TYPES must be grepped to verify the results
# 3rd field ---> tolerance
# 4th field ---> reference result
run.inp 99 1.0e-02 94.0
#EOF

13
tests/QS/regtest-loc_powf/pos_4Au_H2O_4Ne.xyz Normal file
View file

@ -0,0 +1,13 @@
11
ABC 20.000 5.222 6.031
H 2.9511341151 1.2986179214 2.0687351408
H 4.0189849780 2.0237200504 1.0684299686
O 3.6013024600 1.1516070983 1.3233030935
Au 0.0000000000 0.0000000000 0.0000000000
Au 0.0000000000 0.0000000000 3.0155867401
Au 0.0000000000 2.6111560658 1.5080350826
Au 0.0000000000 2.6111560658 4.5231383976
Ne 10.0000000000 0.0000000000 0.0000000000
Ne 10.0000000000 0.0000000000 3.0155867401
Ne 10.0000000000 2.6111560658 1.5080350826
Ne 10.0000000000 2.6111560658 4.5231383976

136
tests/QS/regtest-loc_powf/run.inp Normal file
View file

@ -0,0 +1,136 @@
&GLOBAL
PROJECT run
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_SET
BASIS_SET_FILE_NAME GTH_BASIS_SETS
POTENTIAL_FILE_NAME GTH_POTENTIALS
CHARGE 0
MULTIPLICITY 1
SURFACE_DIPOLE_CORRECTION T
SURF_DIP_DIR X
SURF_DIP_POS 20.000
SURF_DIP_SWITCH T
CORE_CORR_DIP F
&MGRID
NGRIDS 5
CUTOFF [Ry] 300
REL_CUTOFF [Ry] 40
&END
&QS
METHOD GPW
EPS_DEFAULT 1.0E-8
&END
&POISSON
PERIODIC XYZ
&END
&SCF
SCF_GUESS ATOMIC
MAX_SCF 100
EPS_SCF 1.0E-4
ADDED_MOS 10
&DIAGONALIZATION
ALGORITHM STANDARD
&END
&MIXING
METHOD BROYDEN_MIXING
ALPHA 0.69
BETA 1.55
NBROYDEN 13
&END
&SMEAR
METHOD ENERGY_WINDOW
WINDOW_SIZE 0.0037
&END
&OUTER_SCF ! repeat the inner SCF cycle 10 times
MAX_SCF 10
EPS_SCF 1.0E-4 ! must match the above
&END
&PRINT
&RESTART OFF
&END RESTART
&RESTART_HISTORY OFF
&END RESTART_HISTORY
&END PRINT
&END SCF
&XC
&XC_FUNCTIONAL
&PBE
&END
&END XC_FUNCTIONAL
&END XC
&LOCALIZE T
MAX_ITER 200
OUT_ITER_EACH 100
EPS_LOCALIZATION 94.0
METHOD GAPO
CPO_GUESS ATOMIC
CPO_GUESS_SPACE WAN
CG_PO T
NEXTRA 2
STATES MIXED
&PRINT
&WANNIER_CUBES OFF
&END WANNIER_CUBES
&WANNIER_CENTERS ON
FILENAME =wannier.xyz
IONS+CENTERS T
FORMAT XMOL
&EACH
QS_SCF 50
&END EACH
&END WANNIER_CENTERS
&WANNIER_SPREADS OFF
&END WANNIER_SPREADS
&END PRINT
&END LOCALIZE
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 30.000 5.222 6.031
PERIODIC XYZ
&END CELL
&TOPOLOGY
COORD_FILE_NAME ./pos_4Au_H2O_4Ne.xyz
COORD_FILE_FORMAT XYZ
&CENTER_COORDINATES FALSE
&END
&END
&KIND O
ELEMENT O
BASIS_SET DZVP-GTH-q6
POTENTIAL GTH-PBE-q6
&END KIND
&KIND H
ELEMENT H
BASIS_SET DZVP-GTH-q1
POTENTIAL GTH-PBE-q1
&END KIND
&KIND Au
ELEMENT Au
BASIS_SET TZ-GTH
POTENTIAL GTH-PBE-q11
&END KIND
&KIND Ne
ELEMENT Ne
BASIS_SET DZVP-GTH-q8
POTENTIAL GTH-PBE-q8
&END KIND
&END SUBSYS
&END FORCE_EVAL

1
tests/TEST_DIRS
View file

@ -15,6 +15,7 @@ QS/regtest-double-hybrid-stress-numer-meta libxc
QS/regtest-double-hybrid-grad-numer-meta libxc
QS/regtest-double-hybrid-stress-meta libxc
QS/regtest-double-hybrid-grad-meta libxc
QS/regtest-loc_powf
QS/regtest-admm-libxc libint libxc
QS/regtest-mp2-admm-stress-numer libint
QS/regtest-mp2-admm-grad-numer libint

3
tests/TEST_TYPES
View file

@ -1,4 +1,4 @@
98
99
Total energy:!3
MD| Potential energy!5
Total energy \[eV\]:!4
@ -97,6 +97,7 @@ Fermi energy: !3
APT | 1 2 !7
r(1) !3
GW bandgap (eV) !4
Total Spread (Berry) : !6
#
# these are the tests the can be selected for regtesting.
# do regtest will grep for test_grep (first column) and look if the numeric value