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

Direct unconstrained variable-metric localization of one-electron orbitals (#762)

Browse source 
* NLMO reimplemented, compiles and works for 2 test files
* NLMO refactoring to enable virtual localization in the future
* Diagonal BFGS is implemented and tested on Si. Same convergence as SD
* Adds L-BFGS: code compiles, 1st test fails
* L-BFGS minor bug resolved. Energy goes up in 1st tiny L-BFGS line search
* q: remove minus sign
* Compiled and tested refactoring NLMO objct func with virtual and compact implementation
* Modified NLMOs keyword description
* Tested code and added NLMO input file
* Fix NLMO memory leaks, enter correct reference values for NLMO tests
* Fix print bug in NLMO PCG optimizer
* Fix all valgrind errors related to uninitialized variables
* Fix unused dummy argument in almo_scf_lbfgs_types.F
This commit is contained in:
rzk1 2020-02-11 11:26:17 -05:00 • committed by GitHub
parent abecc61cbb
commit 0482296e8b
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
16 changed files with 2988 additions and 220 deletions
Download patch file Download diff file Expand all files Collapse all files

482
src/almo_scf.F
View file

@ -5,7 +5,7 @@
! **************************************************************************************************
!> \brief Routines for all ALMO-based SCF methods
!> RZK-warning marks unresolved issues
!> 'RZK-warning' marks unresolved issues
!> \par History
!> 2011.05 created [Rustam Z Khaliullin]
!> \author Rustam Z Khaliullin
@ -17,6 +17,7 @@ MODULE almo_scf
distribute_domains,&
orthogonalize_mos
USE almo_scf_optimizer, ONLY: almo_scf_block_diagonal,&
almo_scf_construct_nlmos,&
almo_scf_xalmo_eigensolver,&
almo_scf_xalmo_pcg,&
almo_scf_xalmo_trustr
@ -44,6 +45,7 @@ MODULE almo_scf
USE cp_blacs_env, ONLY: cp_blacs_env_release,&
cp_blacs_env_retain
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
USE cp_fm_types, ONLY: cp_fm_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
@ -53,9 +55,11 @@ MODULE almo_scf
cp_para_env_retain
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: &
dbcsr_add_on_diag, dbcsr_binary_read, dbcsr_checksum, dbcsr_copy, dbcsr_create, &
dbcsr_add, dbcsr_add_on_diag, dbcsr_binary_read, dbcsr_checksum, dbcsr_copy, dbcsr_create, &
dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, &
dbcsr_get_info, dbcsr_get_stored_coordinates, dbcsr_multiply, dbcsr_nblkcols_total, &
dbcsr_get_info, dbcsr_get_stored_coordinates, dbcsr_init_random, &
dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_nblkcols_total, &
dbcsr_nblkrows_total, dbcsr_p_type, dbcsr_release, dbcsr_reserve_block2d, dbcsr_scale, &
dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric, dbcsr_work_create
USE domain_submatrix_methods, ONLY: init_submatrices,&
@ -74,6 +78,7 @@ MODULE almo_scf
USE kinds, ONLY: default_path_length,&
dp
USE mathlib, ONLY: binomial
USE message_passing, ONLY: mp_sum
USE molecule_types, ONLY: get_molecule_set_info,&
molecule_type
USE mscfg_types, ONLY: get_matrix_from_submatrices,&
@ -116,7 +121,7 @@ CONTAINS
! **************************************************************************************************
SUBROUTINE almo_entry_scf(qs_env, calc_forces)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL :: calc_forces
LOGICAL, INTENT(IN) :: calc_forces
CHARACTER(len=*), PARAMETER :: routineN = 'almo_entry_scf', routineP = moduleN//':'//routineN
@ -142,6 +147,9 @@ CONTAINS
! allow electron delocalization
CALL almo_scf_delocalization(qs_env, almo_scf_env)
! construct NLMOs
CALL construct_nlmos(qs_env, almo_scf_env)
! electron correlation methods
!CALL almo_correlation_main(qs_env,almo_scf_env)
@ -167,8 +175,8 @@ CONTAINS
! **************************************************************************************************
SUBROUTINE almo_scf_init(qs_env, almo_scf_env, calc_forces)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(almo_scf_env_type) :: almo_scf_env
LOGICAL :: calc_forces
TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
LOGICAL, INTENT(IN) :: calc_forces
CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_init', routineP = moduleN//':'//routineN
@ -198,6 +206,7 @@ CONTAINS
almo_scf_env%opt_xalmo_diis%optimizer_type = optimizer_diis
almo_scf_env%opt_xalmo_pcg%optimizer_type = optimizer_pcg
almo_scf_env%opt_xalmo_trustr%optimizer_type = optimizer_trustr
almo_scf_env%opt_nlmo_pcg%optimizer_type = optimizer_pcg
almo_scf_env%opt_xalmo_newton_pcg_solver%optimizer_type = optimizer_lin_eq_pcg
! get info from the qs_env
@ -508,7 +517,7 @@ CONTAINS
! **************************************************************************************************
SUBROUTINE almo_scf_initial_guess(qs_env, almo_scf_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(almo_scf_env_type) :: almo_scf_env
TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_initial_guess', &
routineP = moduleN//':'//routineN
@ -760,7 +769,7 @@ CONTAINS
!> \author Rustam Khaliullin
! **************************************************************************************************
SUBROUTINE almo_scf_store_extrapolation_data(almo_scf_env)
TYPE(almo_scf_env_type) :: almo_scf_env
TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_store_extrapolation_data', &
routineP = moduleN//':'//routineN
@ -1259,7 +1268,7 @@ CONTAINS
! **************************************************************************************************
SUBROUTINE almo_scf_main(qs_env, almo_scf_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(almo_scf_env_type) :: almo_scf_env
TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_main', routineP = moduleN//':'//routineN
@ -1351,7 +1360,7 @@ CONTAINS
SUBROUTINE almo_scf_delocalization(qs_env, almo_scf_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(almo_scf_env_type) :: almo_scf_env
TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_delocalization', &
routineP = moduleN//':'//routineN
@ -1658,6 +1667,444 @@ CONTAINS
END SUBROUTINE almo_scf_delocalization
! **************************************************************************************************
!> \brief orbital localization
!> \param qs_env ...
!> \param almo_scf_env ...
!> \par History
!> 2018.09 created [Ziling Luo]
!> \author Ziling Luo
! **************************************************************************************************
SUBROUTINE construct_nlmos(qs_env, almo_scf_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
CHARACTER(len=*), PARAMETER :: routineN = 'construct_nlmos', &
routineP = moduleN//':'//routineN
INTEGER :: ispin
IF (almo_scf_env%construct_nlmos) THEN
DO ispin = 1, almo_scf_env%nspins
CALL orthogonalize_mos(ket=almo_scf_env%matrix_t(ispin), &
overlap=almo_scf_env%matrix_sigma(ispin), &
metric=almo_scf_env%matrix_s(1), &
retain_locality=.FALSE., &
only_normalize=.FALSE., &
nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
eps_filter=almo_scf_env%eps_filter, &
order_lanczos=almo_scf_env%order_lanczos, &
eps_lanczos=almo_scf_env%eps_lanczos, &
max_iter_lanczos=almo_scf_env%max_iter_lanczos)
ENDDO
CALL construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals=.FALSE.)
IF (almo_scf_env%opt_nlmo_pcg%opt_penalty%virtual_nlmos) THEN
CALL construct_virtuals(almo_scf_env)
CALL construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals=.TRUE.)
ENDIF
IF (almo_scf_env%opt_nlmo_pcg%opt_penalty%compactification_filter_start .GT. 0.0_dp) THEN
CALL nlmo_compactification(qs_env, almo_scf_env, almo_scf_env%matrix_t)
ENDIF
ENDIF
END SUBROUTINE construct_nlmos
! **************************************************************************************************
!> \brief Calls NLMO optimization
!> \param qs_env ...
!> \param almo_scf_env ...
!> \param virtuals ...
!> \par History
!> 2019.10 created [Ziling Luo]
!> \author Ziling Luo
! **************************************************************************************************
SUBROUTINE construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
LOGICAL, INTENT(IN) :: virtuals
CHARACTER(len=*), PARAMETER :: routineN = 'construct_nlmos_wrapper', &
routineP = moduleN//':'//routineN
REAL(KIND=dp) :: det_diff, prev_determinant
almo_scf_env%overlap_determinant = 1.0
! KEEP: initial_vol_coeff = almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength
almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength = &
-1.0_dp*almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength !NEW1
! loop over the strength of the orthogonalization penalty
prev_determinant = 10.0_dp
DO WHILE (almo_scf_env%overlap_determinant .GT. almo_scf_env%opt_nlmo_pcg%opt_penalty%final_determinant)
IF (.NOT. virtuals) THEN
CALL almo_scf_construct_nlmos(qs_env=qs_env, &
optimizer=almo_scf_env%opt_nlmo_pcg, &
matrix_s=almo_scf_env%matrix_s(1), &
matrix_mo_in=almo_scf_env%matrix_t, &
matrix_mo_out=almo_scf_env%matrix_t, &
template_matrix_sigma=almo_scf_env%matrix_sigma_inv, &
overlap_determinant=almo_scf_env%overlap_determinant, &
mat_distr_aos=almo_scf_env%mat_distr_aos, &
virtuals=virtuals, &
eps_filter=almo_scf_env%eps_filter)
ELSE
CALL almo_scf_construct_nlmos(qs_env=qs_env, &
optimizer=almo_scf_env%opt_nlmo_pcg, &
matrix_s=almo_scf_env%matrix_s(1), &
matrix_mo_in=almo_scf_env%matrix_v, &
matrix_mo_out=almo_scf_env%matrix_v, &
template_matrix_sigma=almo_scf_env%matrix_sigma_vv, &
overlap_determinant=almo_scf_env%overlap_determinant, &
mat_distr_aos=almo_scf_env%mat_distr_aos, &
virtuals=virtuals, &
eps_filter=almo_scf_env%eps_filter)
ENDIF
det_diff = prev_determinant - almo_scf_env%overlap_determinant
almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength = &
almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength/ &
ABS(almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength_dec_factor)
IF (det_diff < almo_scf_env%opt_nlmo_pcg%opt_penalty%determinant_tolerance) THEN
EXIT
ENDIF
prev_determinant = almo_scf_env%overlap_determinant
ENDDO
END SUBROUTINE construct_nlmos_wrapper
! **************************************************************************************************
!> \brief Construct virtual orbitals
!> \param almo_scf_env ...
!> \par History
!> 2019.10 created [Ziling Luo]
!> \author Ziling Luo
! **************************************************************************************************
SUBROUTINE construct_virtuals(almo_scf_env)
TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
CHARACTER(len=*), PARAMETER :: routineN = 'construct_virtuals', &
routineP = moduleN//':'//routineN
INTEGER :: ispin, n
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
TYPE(dbcsr_type) :: tempNV1, tempVOcc1, tempVOcc2, tempVV1, &
tempVV2
DO ispin = 1, almo_scf_env%nspins
CALL dbcsr_create(tempNV1, &
template=almo_scf_env%matrix_v(ispin), &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_create(tempVOcc1, &
template=almo_scf_env%matrix_vo(ispin), &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_create(tempVOcc2, &
template=almo_scf_env%matrix_vo(ispin), &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_create(tempVV1, &
template=almo_scf_env%matrix_sigma_vv(ispin), &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_create(tempVV2, &
template=almo_scf_env%matrix_sigma_vv(ispin), &
matrix_type=dbcsr_type_no_symmetry)
! Generate random virtual matrix
CALL dbcsr_init_random(almo_scf_env%matrix_v(ispin), &
keep_sparsity=.FALSE.)
! Project the orbital subspace out
CALL dbcsr_multiply("N", "N", 1.0_dp, &
almo_scf_env%matrix_s(1), &
almo_scf_env%matrix_v(ispin), &
0.0_dp, tempNV1, &
filter_eps=almo_scf_env%eps_filter)
CALL dbcsr_multiply("T", "N", 1.0_dp, &
tempNV1, &
almo_scf_env%matrix_t(ispin), &
0.0_dp, tempVOcc1, &
filter_eps=almo_scf_env%eps_filter)
CALL dbcsr_multiply("N", "N", 1.0_dp, &
tempVOcc1, &
almo_scf_env%matrix_sigma_inv(ispin), &
0.0_dp, tempVOcc2, &
filter_eps=almo_scf_env%eps_filter)
CALL dbcsr_multiply("N", "T", 1.0_dp, &
almo_scf_env%matrix_t(ispin), &
tempVOcc2, &
0.0_dp, tempNV1, &
filter_eps=almo_scf_env%eps_filter)
CALL dbcsr_add(almo_scf_env%matrix_v(ispin), tempNV1, 1.0_dp, -1.0_dp)
! compute VxV overlap
CALL dbcsr_multiply("N", "N", 1.0_dp, &
almo_scf_env%matrix_s(1), &
almo_scf_env%matrix_v(ispin), &
0.0_dp, tempNV1, &
filter_eps=almo_scf_env%eps_filter)
CALL dbcsr_multiply("T", "N", 1.0_dp, &
almo_scf_env%matrix_v(ispin), &
tempNV1, &
0.0_dp, tempVV1, &
filter_eps=almo_scf_env%eps_filter)
CALL orthogonalize_mos(ket=almo_scf_env%matrix_v(ispin), &
overlap=tempVV1, &
metric=almo_scf_env%matrix_s(1), &
retain_locality=.FALSE., &
only_normalize=.FALSE., &
nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
eps_filter=almo_scf_env%eps_filter, &
order_lanczos=almo_scf_env%order_lanczos, &
eps_lanczos=almo_scf_env%eps_lanczos, &
max_iter_lanczos=almo_scf_env%max_iter_lanczos)
! compute VxV block of the KS matrix
CALL dbcsr_multiply("N", "N", 1.0_dp, &
almo_scf_env%matrix_ks(ispin), &
almo_scf_env%matrix_v(ispin), &
0.0_dp, tempNV1, &
filter_eps=almo_scf_env%eps_filter)
CALL dbcsr_multiply("T", "N", 1.0_dp, &
almo_scf_env%matrix_v(ispin), &
tempNV1, &
0.0_dp, tempVV1, &
filter_eps=almo_scf_env%eps_filter)
CALL dbcsr_get_info(tempVV1, nfullrows_total=n)
ALLOCATE (eigenvalues(n))
CALL cp_dbcsr_syevd(tempVV1, tempVV2, &
eigenvalues, &
para_env=almo_scf_env%para_env, &
blacs_env=almo_scf_env%blacs_env)
DEALLOCATE (eigenvalues)
CALL dbcsr_multiply("N", "N", 1.0_dp, &
almo_scf_env%matrix_v(ispin), &
tempVV2, &
0.0_dp, tempNV1, &
filter_eps=almo_scf_env%eps_filter)
CALL dbcsr_copy(almo_scf_env%matrix_v(ispin), tempNV1)
CALL dbcsr_release(tempNV1)
CALL dbcsr_release(tempVOcc1)
CALL dbcsr_release(tempVOcc2)
CALL dbcsr_release(tempVV1)
CALL dbcsr_release(tempVV2)
ENDDO
END SUBROUTINE construct_virtuals
! **************************************************************************************************
!> \brief Compactify (set small blocks to zero) orbitals
!> \param qs_env ...
!> \param almo_scf_env ...
!> \param matrix ...
!> \par History
!> 2019.10 created [Ziling Luo]
!> \author Ziling Luo
! **************************************************************************************************
SUBROUTINE nlmo_compactification(qs_env, almo_scf_env, matrix)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:), &
INTENT(IN) :: matrix
CHARACTER(len=*), PARAMETER :: routineN = 'nlmo_compactification', &
routineP = moduleN//':'//routineN
INTEGER :: iblock_col, iblock_col_size, iblock_row, &
iblock_row_size, icol, irow, ispin, &
Ncols, Nrows, nspins, para_group, &
unit_nr
LOGICAL :: element_by_element
REAL(KIND=dp) :: energy, eps_local, eps_start, &
max_element, spin_factor
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occ, retained
REAL(kind=dp), DIMENSION(:, :), POINTER :: data_p
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_iterator_type) :: iter
TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_p_tmp, matrix_t_tmp
! define the output_unit
logger => cp_get_default_logger()
IF (logger%para_env%ionode) THEN
unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
ELSE
unit_nr = -1
ENDIF
nspins = SIZE(matrix)
element_by_element = .FALSE.
IF (nspins .EQ. 1) THEN
spin_factor = 2.0_dp
ELSE
spin_factor = 1.0_dp
ENDIF
ALLOCATE (matrix_t_tmp(nspins))
ALLOCATE (matrix_p_tmp(nspins))
ALLOCATE (retained(nspins))
ALLOCATE (occ(2))
DO ispin = 1, nspins
! init temporary storage
CALL dbcsr_create(matrix_t_tmp(ispin), &
template=matrix(ispin), &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_copy(matrix_t_tmp(ispin), matrix(ispin))
CALL dbcsr_create(matrix_p_tmp(ispin), &
template=almo_scf_env%matrix_p(ispin), &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_copy(matrix_p_tmp(ispin), almo_scf_env%matrix_p(ispin))
ENDDO
IF (unit_nr > 0) THEN
WRITE (unit_nr, *)
WRITE (unit_nr, '(T2,A)') &
"Energy dependence on the (block-by-block) filtering of the NLMO coefficients"
IF (unit_nr > 0) WRITE (unit_nr, '(T2,A13,A20,A20,A25)') &
"EPS filter", "Occupation Alpha", "Occupation Beta", "Energy"
ENDIF
eps_start = almo_scf_env%opt_nlmo_pcg%opt_penalty%compactification_filter_start
eps_local = MAX(eps_start, 10E-14_dp)
DO
IF (eps_local > 0.11_dp) EXIT
DO ispin = 1, nspins
retained(ispin) = 0
CALL dbcsr_work_create(matrix_t_tmp(ispin), work_mutable=.TRUE.)
CALL dbcsr_iterator_start(iter, matrix_t_tmp(ispin))
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, &
row_size=iblock_row_size, col_size=iblock_col_size)
DO icol = 1, iblock_col_size
IF (element_by_element) THEN
DO irow = 1, iblock_row_size
IF (ABS(data_p(irow, icol)) .LT. eps_local) THEN
data_p(irow, icol) = 0.0_dp
ELSE
retained(ispin) = retained(ispin) + 1
ENDIF
ENDDO
ELSE ! rows are blocked
max_element = 0.0_dp
DO irow = 1, iblock_row_size
IF (ABS(data_p(irow, icol)) .GT. max_element) THEN
max_element = ABS(data_p(irow, icol))
ENDIF
ENDDO
IF (max_element .LT. eps_local) THEN
DO irow = 1, iblock_row_size
data_p(irow, icol) = 0.0_dp
ENDDO
ELSE
retained(ispin) = retained(ispin) + iblock_row_size
ENDIF
ENDIF ! block rows?
ENDDO ! icol
ENDDO ! iterator
CALL dbcsr_iterator_stop(iter)
CALL dbcsr_finalize(matrix_t_tmp(ispin))
CALL dbcsr_filter(matrix_t_tmp(ispin), eps_local)
CALL dbcsr_get_info(matrix_t_tmp(ispin), group=para_group, &
nfullrows_total=Nrows, &
nfullcols_total=Ncols)
CALL mp_sum(retained(ispin), para_group)
!devide by the total no. elements
occ(ispin) = retained(ispin)/Nrows/Ncols
! compute the global projectors (for the density matrix)
CALL almo_scf_t_to_proj( &
t=matrix_t_tmp(ispin), &
p=matrix_p_tmp(ispin), &
eps_filter=almo_scf_env%eps_filter, &
orthog_orbs=.FALSE., &
nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
s=almo_scf_env%matrix_s(1), &
sigma=almo_scf_env%matrix_sigma(ispin), &
sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
use_guess=.FALSE., &
algorithm=almo_scf_env%sigma_inv_algorithm, &
inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
inverse_accelerator=almo_scf_env%order_lanczos, &
eps_lanczos=almo_scf_env%eps_lanczos, &
max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
para_env=almo_scf_env%para_env, &
blacs_env=almo_scf_env%blacs_env)
! compute dm from the projector(s)
CALL dbcsr_scale(matrix_p_tmp(ispin), spin_factor)
ENDDO
! the KS matrix is updated outside the spin loop
CALL almo_dm_to_almo_ks(qs_env, &
matrix_p_tmp, &
almo_scf_env%matrix_ks, &
energy, &
almo_scf_env%eps_filter, &
almo_scf_env%mat_distr_aos)
IF (nspins .LT. 2) occ(2) = occ(1)
IF (unit_nr > 0) WRITE (unit_nr, '(T2,E13.3,F20.10,F20.10,F25.15)') &
eps_local, occ(1), occ(2), energy
eps_local = 2.0_dp*eps_local
ENDDO
DO ispin = 1, nspins
CALL dbcsr_release(matrix_t_tmp(ispin))
CALL dbcsr_release(matrix_p_tmp(ispin))
ENDDO
DEALLOCATE (matrix_t_tmp)
DEALLOCATE (matrix_p_tmp)
DEALLOCATE (occ)
DEALLOCATE (retained)
END SUBROUTINE nlmo_compactification
! *****************************************************************************
!> \brief after SCF we have the final density and KS matrices compute various
!> post-scf quantities
@ -1669,7 +2116,7 @@ CONTAINS
! **************************************************************************************************
SUBROUTINE almo_scf_post(qs_env, almo_scf_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(almo_scf_env_type) :: almo_scf_env
TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_post', routineP = moduleN//':'//routineN
@ -1687,6 +2134,7 @@ CONTAINS
! orthogonalize orbitals before returning them to QS
ALLOCATE (matrix_t_processed(almo_scf_env%nspins))
!ALLOCATE (matrix_v_processed(almo_scf_env%nspins))
DO ispin = 1, almo_scf_env%nspins
@ -1697,6 +2145,13 @@ CONTAINS
CALL dbcsr_copy(matrix_t_processed(ispin), &
almo_scf_env%matrix_t(ispin))
!CALL dbcsr_create(matrix_v_processed(ispin), &
! template=almo_scf_env%matrix_v(ispin), &
! matrix_type=dbcsr_type_no_symmetry)
!CALL dbcsr_copy(matrix_v_processed(ispin), &
! almo_scf_env%matrix_v(ispin))
IF (almo_scf_env%return_orthogonalized_mos) THEN
CALL orthogonalize_mos(ket=matrix_t_processed(ispin), &
@ -1736,6 +2191,9 @@ CONTAINS
CALL copy_dbcsr_to_fm(matrix_t_processed(ispin), mo_coeff)
CALL dbcsr_release(matrix_t_processed(ispin))
ENDDO
DO ispin = 1, almo_scf_env%nspins
CALL dbcsr_release(matrix_t_processed(ispin))
ENDDO
DEALLOCATE (matrix_t_processed)
! calculate post scf properties
@ -2233,7 +2691,7 @@ CONTAINS
! **************************************************************************************************
SUBROUTINE almo_scf_clean_up(almo_scf_env)
TYPE(almo_scf_env_type) :: almo_scf_env
TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_clean_up', &
routineP = moduleN//':'//routineN

77
src/almo_scf_env_methods.F
View file

@ -15,10 +15,9 @@ MODULE almo_scf_env_methods
almo_scf_env_type
USE cp_control_types, ONLY: dft_control_type
USE input_constants, ONLY: &
almo_constraint_distance, almo_deloc_none, almo_deloc_scf, almo_deloc_x, &
almo_deloc_x_then_scf, almo_deloc_xalmo_1diag, almo_domain_layout_atomic, &
almo_domain_layout_molecular, almo_frz_crystal, almo_mat_distr_molecular, &
almo_occ_vol_penalty_none, almo_scf_diag, almo_scf_skip, almo_scf_trustr, cg_hager_zhang, &
almo_constraint_distance, almo_deloc_none, almo_deloc_xalmo_1diag, &
almo_domain_layout_atomic, almo_domain_layout_molecular, almo_frz_crystal, &
almo_mat_distr_molecular, almo_scf_diag, almo_scf_skip, almo_scf_trustr, cg_hager_zhang, &
do_bondparm_vdw, molecular_guess, tensor_orthogonal, virt_full, virt_minimal, virt_number, &
xalmo_trial_r0_out
USE input_section_types, ONLY: section_vals_get_subs_vals,&
@ -107,8 +106,9 @@ CONTAINS
INTEGER :: handle
TYPE(section_vals_type), POINTER :: almo_analysis_section, almo_opt_diis_section, &
almo_opt_pcg_section, almo_scf_section, matrix_iterate_section, penalty_section, &
xalmo_opt_newton_pcg_section, xalmo_opt_pcg_section, xalmo_opt_trustr_section
almo_opt_pcg_section, almo_scf_section, matrix_iterate_section, nlmo_opt_pcg_section, &
penalty_section, xalmo_opt_newton_pcg_section, xalmo_opt_pcg_section, &
xalmo_opt_trustr_section
CALL timeset(routineN, handle)
@ -121,12 +121,13 @@ CONTAINS
"XALMO_OPTIMIZER_PCG")
xalmo_opt_trustr_section => section_vals_get_subs_vals(almo_scf_section, &
"XALMO_OPTIMIZER_TRUSTR")
nlmo_opt_pcg_section => section_vals_get_subs_vals(almo_scf_section, &
"NLMO_OPTIMIZER_PCG")
almo_analysis_section => section_vals_get_subs_vals(almo_scf_section, "ANALYSIS")
xalmo_opt_newton_pcg_section => section_vals_get_subs_vals(xalmo_opt_pcg_section, &
"XALMO_NEWTON_PCG_SOLVER")
matrix_iterate_section => section_vals_get_subs_vals(almo_scf_section, &
"MATRIX_ITERATE")
penalty_section => section_vals_get_subs_vals(almo_scf_section, "PENALTY")
! read user input
! common ALMO options
@ -154,6 +155,8 @@ CONTAINS
almo_scf_env%xalmo_extrapolation_order = MAX(0, almo_scf_env%xalmo_extrapolation_order)
CALL section_vals_val_get(almo_scf_section, "RETURN_ORTHOGONALIZED_MOS", &
l_val=almo_scf_env%return_orthogonalized_mos)
CALL section_vals_val_get(almo_scf_section, "CONSTRUCT_NLMOS", &
l_val=almo_scf_env%construct_nlmos)
CALL section_vals_val_get(matrix_iterate_section, "EPS_LANCZOS", &
r_val=almo_scf_env%eps_lanczos)
@ -243,6 +246,49 @@ CONTAINS
CALL section_vals_val_get(xalmo_opt_pcg_section, "PRECONDITIONER", &
i_val=almo_scf_env%opt_xalmo_pcg%preconditioner)
penalty_section => section_vals_get_subs_vals(nlmo_opt_pcg_section, "PENALTY")
CALL section_vals_val_get(nlmo_opt_pcg_section, "EPS_ERROR", &
r_val=almo_scf_env%opt_nlmo_pcg%eps_error)
CALL section_vals_val_get(nlmo_opt_pcg_section, "MAX_ITER", &
i_val=almo_scf_env%opt_nlmo_pcg%max_iter)
CALL section_vals_val_get(nlmo_opt_pcg_section, "EPS_ERROR_EARLY", &
r_val=almo_scf_env%opt_nlmo_pcg%eps_error_early)
CALL section_vals_val_get(nlmo_opt_pcg_section, "MAX_ITER_EARLY", &
i_val=almo_scf_env%opt_nlmo_pcg%max_iter_early)
CALL section_vals_val_get(nlmo_opt_pcg_section, "MAX_ITER_OUTER_LOOP", &
i_val=almo_scf_env%opt_nlmo_pcg%max_iter_outer_loop)
CALL section_vals_val_get(nlmo_opt_pcg_section, "LIN_SEARCH_EPS_ERROR", &
r_val=almo_scf_env%opt_nlmo_pcg%lin_search_eps_error)
CALL section_vals_val_get(nlmo_opt_pcg_section, "LIN_SEARCH_STEP_SIZE_GUESS", &
r_val=almo_scf_env%opt_nlmo_pcg%lin_search_step_size_guess)
CALL section_vals_val_get(nlmo_opt_pcg_section, "PRECOND_FILTER_THRESHOLD", &
r_val=almo_scf_env%opt_nlmo_pcg%neglect_threshold)
CALL section_vals_val_get(nlmo_opt_pcg_section, "CONJUGATOR", &
i_val=almo_scf_env%opt_nlmo_pcg%conjugator)
CALL section_vals_val_get(nlmo_opt_pcg_section, "PRECONDITIONER", &
i_val=almo_scf_env%opt_nlmo_pcg%preconditioner)
CALL section_vals_val_get(penalty_section, &
"OPERATOR", &
i_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%operator_type)
CALL section_vals_val_get(penalty_section, &
"PENALTY_STRENGTH", &
r_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength)
CALL section_vals_val_get(penalty_section, &
"PENALTY_STRENGTH_DECREASE_FACTOR", &
r_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength_dec_factor)
CALL section_vals_val_get(penalty_section, &
"DETERMINANT_TOLERANCE", &
r_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%determinant_tolerance)
CALL section_vals_val_get(penalty_section, &
"FINAL_DETERMINANT", &
r_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%final_determinant)
CALL section_vals_val_get(penalty_section, &
"COMPACTIFICATION_FILTER_START", &
r_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%compactification_filter_start)
CALL section_vals_val_get(penalty_section, &
"VIRTUAL_NLMOS", &
l_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%virtual_nlmos)
CALL section_vals_val_get(xalmo_opt_newton_pcg_section, "EPS_ERROR", &
r_val=almo_scf_env%opt_xalmo_newton_pcg_solver%eps_error)
CALL section_vals_val_get(xalmo_opt_newton_pcg_section, "MAX_ITER", &
@ -252,13 +298,6 @@ CONTAINS
CALL section_vals_val_get(xalmo_opt_newton_pcg_section, "PRECONDITIONER", &
i_val=almo_scf_env%opt_xalmo_newton_pcg_solver%preconditioner)
CALL section_vals_val_get(penalty_section, &
"OCCUPIED_VOLUME_PENALTY_COEFF", &
r_val=almo_scf_env%penalty%occ_vol_coeff)
CALL section_vals_val_get(penalty_section, &
"OCCUPIED_VOLUME_PENALTY_METHOD", &
i_val=almo_scf_env%penalty%occ_vol_method)
CALL section_vals_val_get(almo_analysis_section, "_SECTION_PARAMETERS_", &
l_val=almo_scf_env%almo_analysis%do_analysis)
CALL section_vals_val_get(almo_analysis_section, "FROZEN_MO_ENERGY_TERM", &
@ -504,16 +543,6 @@ CONTAINS
ENDIF ! end analysis settings
! check penalty settings
IF (almo_scf_env%penalty%occ_vol_method .NE. almo_occ_vol_penalty_none) THEN
IF (almo_scf_env%deloc_method .NE. almo_deloc_x .AND. &
almo_scf_env%deloc_method .NE. almo_deloc_scf .AND. &
almo_scf_env%deloc_method .NE. almo_deloc_x_then_scf) THEN
CALL cp_abort(__LOCATION__, &
"Occupied volume penalty seems to work only with completely delocalized orbitals")
ENDIF
ENDIF ! end penalty settings
CALL timestop(handle)
END SUBROUTINE almo_scf_init_read_write_input

330
src/almo_scf_lbfgs_types.F Normal file
View file

@ -0,0 +1,330 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2020 CP2K developers group !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Limited memory BFGS
!> \par History
!> 2019.10 created [Rustam Z Khaliullin]
!> \author Rustam Z Khaliullin
! **************************************************************************************************
MODULE almo_scf_lbfgs_types
!USE cp_external_control, ONLY: external_control
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_unit_nr
!USE cp_log_handling, ONLY: cp_to_string
USE dbcsr_api, ONLY: dbcsr_add,&
dbcsr_copy,&
dbcsr_create,&
dbcsr_dot,&
dbcsr_release,&
dbcsr_scale,&
dbcsr_type
USE kinds, ONLY: dp
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'almo_scf_lbfgs_types'
PUBLIC :: lbfgs_seed, &
lbfgs_create, &
lbfgs_release, &
lbfgs_get_direction, &
lbfgs_history_type
TYPE lbfgs_history_type
INTEGER :: nstore
! istore counts the total number of action=2 pushes
! istore is designed to become more than nstore eventually
! there are two counters: the main variable and gradient
INTEGER, DIMENSION(2) :: istore
TYPE(dbcsr_type), DIMENSION(:, :, :), ALLOCATABLE :: matrix
REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: rho
END TYPE lbfgs_history_type
CONTAINS
! **************************************************************************************************
!> \brief interface subroutine to store the first variable/gradient pair
!> \param history ...
!> \param variable ...
!> \param gradient ...
! **************************************************************************************************
SUBROUTINE lbfgs_seed(history, variable, gradient)
TYPE(lbfgs_history_type), INTENT(INOUT) :: history
TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: variable, gradient
CALL lbfgs_history_push(history, variable, vartype=1, action=1)
CALL lbfgs_history_push(history, gradient, vartype=2, action=1)
END SUBROUTINE lbfgs_seed
! **************************************************************************************************
!> \brief interface subroutine to store a variable/gradient pair
!> and predict direction
!> \param history ...
!> \param variable ...
!> \param gradient ...
!> \param direction ...
! **************************************************************************************************
SUBROUTINE lbfgs_get_direction(history, variable, gradient, direction)
TYPE(lbfgs_history_type), INTENT(INOUT) :: history
TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: variable, gradient
TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: direction
! action 2 will calculate delta = (new - old)
! in the last used storage cell
CALL lbfgs_history_push(history, variable, vartype=1, action=2)
CALL lbfgs_history_push(history, gradient, vartype=2, action=2)
! compute rho for the last stored value
CALL lbfgs_history_last_rho(history)
CALL lbfgs_history_direction(history, gradient, direction)
! action 1 will seed the next storage cell
CALL lbfgs_history_push(history, variable, vartype=1, action=1)
CALL lbfgs_history_push(history, gradient, vartype=2, action=1)
END SUBROUTINE lbfgs_get_direction
! **************************************************************************************************
!> \brief create history storage for limited memory bfgs
!> \param history ...
!> \param nspins ...
!> \param nstore ...
! **************************************************************************************************
SUBROUTINE lbfgs_create(history, nspins, nstore)
TYPE(lbfgs_history_type), INTENT(INOUT) :: history
INTEGER, INTENT(IN) :: nspins, nstore
INTEGER :: nallocate
nallocate = MAX(1, nstore)
history%nstore = nallocate
history%istore(:) = 0 ! total number of action-2 pushes
ALLOCATE (history%matrix(nspins, nallocate, 2))
ALLOCATE (history%rho(nspins, nallocate))
END SUBROUTINE lbfgs_create
! **************************************************************************************************
!> \brief release the bfgs history
!> \param history ...
! **************************************************************************************************
SUBROUTINE lbfgs_release(history)
TYPE(lbfgs_history_type), INTENT(INOUT) :: history
INTEGER :: ispin, istore, ivartype
! delete history
DO ispin = 1, SIZE(history%matrix, 1)
DO ivartype = 1, 2
DO istore = 1, MIN(history%istore(ivartype) + 1, history%nstore)
!WRITE(*,*) "ZREL: ispin,istore,vartype", ispin, istore, ivartype
CALL dbcsr_release(history%matrix(ispin, istore, ivartype))
ENDDO
ENDDO
ENDDO
DEALLOCATE (history%matrix)
DEALLOCATE (history%rho)
END SUBROUTINE lbfgs_release
! **************************************************************************************************
!> \brief once all data in the last cell is stored, compute rho
!> \param history ...
! **************************************************************************************************
SUBROUTINE lbfgs_history_last_rho(history)
TYPE(lbfgs_history_type), INTENT(INOUT) :: history
INTEGER :: ispin, istore
!logger => cp_get_default_logger()
!IF (logger%para_env%mepos == logger%para_env%source) THEN
! unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
!ELSE
! unit_nr = -1
!ENDIF
DO ispin = 1, SIZE(history%matrix, 1)
istore = MOD(history%istore(1) - 1, history%nstore) + 1
CALL dbcsr_dot(history%matrix(ispin, istore, 1), &
history%matrix(ispin, istore, 2), &
history%rho(ispin, istore))
history%rho(ispin, istore) = 1.0_dp/history%rho(ispin, istore)
!IF (unit_nr > 0) THEN
! WRITE (unit_nr, *) "Rho in cell ", istore, " is computed ", history%rho(ispin, istore)
!ENDIF
ENDDO ! ispin
END SUBROUTINE lbfgs_history_last_rho
! **************************************************************************************************
!> \brief store data in history
!> vartype - which data piece to store: 1 - variable, 2 - gradient
!> operation - what to do: 1 - erase existing and store new
!> 2 - store = new - existing
!> \param history ...
!> \param matrix ...
!> \param vartype ...
!> \param action ...
! **************************************************************************************************
SUBROUTINE lbfgs_history_push(history, matrix, vartype, action)
TYPE(lbfgs_history_type), INTENT(INOUT) :: history
TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: matrix
INTEGER, INTENT(IN) :: vartype, action
INTEGER :: ispin, istore
!logger => cp_get_default_logger()
!IF (logger%para_env%mepos == logger%para_env%source) THEN
! unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
!ELSE
! unit_nr = -1
!ENDIF
! increase the counter: it moves the pointer to the next cell
! for action==1 this is a "pretend" increase; the pointer will be moved back in the end
history%istore(vartype) = history%istore(vartype) + 1
DO ispin = 1, SIZE(history%matrix, 1)
istore = MOD(history%istore(vartype) - 1, history%nstore) + 1
!IF (unit_nr > 0) THEN
! WRITE (unit_nr, *) "Action ", action, " modifying cell ", istore
!END IF
IF (history%istore(vartype) <= history%nstore .AND. &
action .EQ. 1) THEN
!WRITE(*,*) "ZCRE: ispin,istore,vartype", ispin, istore, vartype
CALL dbcsr_create(history%matrix(ispin, istore, vartype), &
template=matrix(ispin))
!IF (unit_nr > 0) THEN
! WRITE (unit_nr, *) "Creating new matrix..."
!END IF
ENDIF
IF (action .EQ. 1) THEN
CALL dbcsr_copy(history%matrix(ispin, istore, vartype), matrix(ispin))
ELSE
CALL dbcsr_add(history%matrix(ispin, istore, vartype), matrix(ispin), -1.0_dp, 1.0_dp)
ENDIF
ENDDO ! ispin
! allow the pointer to move forward only if deltas are stored (action==2)
! otherwise return the pointer to the previous value
IF (action .EQ. 1) THEN
history%istore(vartype) = history%istore(vartype) - 1
ENDIF
END SUBROUTINE lbfgs_history_push
! **************************************************************************************************
!> \brief use history data to construct dir = -Hinv.grad
!> \param history ...
!> \param gradient ...
!> \param direction ...
! **************************************************************************************************
SUBROUTINE lbfgs_history_direction(history, gradient, direction)
TYPE(lbfgs_history_type), INTENT(INOUT) :: history
TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: gradient
TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: direction
INTEGER :: ispin, istore, iterm, nterms
REAL(KIND=dp) :: beta, gammak
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha
TYPE(dbcsr_type) :: q
!logger => cp_get_default_logger()
!IF (logger%para_env%mepos == logger%para_env%source) THEN
! unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
!ELSE
! unit_nr = -1
!ENDIF
IF (history%istore(1) .NE. history%istore(2)) THEN
CPABORT("BFGS APIs are not used correctly")
ENDIF
nterms = MIN(history%istore(1), history%nstore)
!IF (unit_nr > 0) THEN
! WRITE (unit_nr, *) "L-BFGS terms used: ", nterms
!END IF
ALLOCATE (alpha(nterms))
DO ispin = 1, SIZE(history%matrix, 1)
CALL dbcsr_create(q, template=gradient(ispin))
CALL dbcsr_copy(q, gradient(ispin))
! loop over all stored items
DO iterm = 1, nterms
! location: from recent to oldest stored
istore = MOD(history%istore(1) - iterm, history%nstore) + 1
!IF (unit_nr > 0) THEN
! WRITE (unit_nr, *) "Record locator: ", istore
!END IF
CALL dbcsr_dot(history%matrix(ispin, istore, 1), q, alpha(iterm))
alpha(iterm) = history%rho(ispin, istore)*alpha(iterm)
CALL dbcsr_add(q, history%matrix(ispin, istore, 2), 1.0_dp, -alpha(iterm))
! use the most recent term to
! compute gamma_k, Nocedal (7.20) and then get H0
IF (iterm .EQ. 1) THEN
CALL dbcsr_dot(history%matrix(ispin, istore, 2), history%matrix(ispin, istore, 2), gammak)
gammak = 1.0_dp/(gammak*history%rho(ispin, istore))
!IF (unit_nr > 0) THEN
! WRITE (unit_nr, *) "Gamma_k: ", gammak
!END IF
ENDIF
ENDDO ! iterm, first loop from recent to oldest
! now q stores Nocedal's r = (gamma_k*I).q
CALL dbcsr_scale(q, gammak)
! loop over all stored items
DO iterm = nterms, 1, -1
! location: from oldest to recent stored
istore = MOD(history%istore(1) - iterm, history%nstore) + 1
CALL dbcsr_dot(history%matrix(ispin, istore, 2), q, beta)
beta = history%rho(ispin, istore)*beta
CALL dbcsr_add(q, history%matrix(ispin, istore, 1), 1.0_dp, alpha(iterm) - beta)
ENDDO ! iterm, forst loop from recent to oldest
!RZK-warning: unclear whether q should be multiplied by minus one
CALL dbcsr_scale(q, -1.0)
CALL dbcsr_copy(direction(ispin), q)
CALL dbcsr_release(q)
ENDDO !ispin
DEALLOCATE (alpha)
END SUBROUTINE lbfgs_history_direction
END MODULE almo_scf_lbfgs_types

47
src/almo_scf_methods.F
View file

@ -75,10 +75,53 @@ MODULE almo_scf_methods
distribute_domains, &
almo_scf_ks_to_ks_xx, &
construct_domain_r_down, &
xalmo_initial_guess
xalmo_initial_guess, &
fill_matrix_with_ones
CONTAINS
! **************************************************************************************************
!> \brief Fill all matrix blocks with 1.0_dp
!> \param matrix ...
!> \par History
!> 2019.09 created [Rustam Z Khaliullin]
!> \author Rustam Z Khaliullin
! **************************************************************************************************
SUBROUTINE fill_matrix_with_ones(matrix)
TYPE(dbcsr_type), INTENT(INOUT) :: matrix
INTEGER :: col, hold, iblock_col, iblock_row, &
mynode, nblkcols_tot, nblkrows_tot, row
LOGICAL :: tr
REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_new_block
TYPE(dbcsr_distribution_type) :: dist
CALL dbcsr_get_info(matrix, distribution=dist)
CALL dbcsr_distribution_get(dist, mynode=mynode)
CALL dbcsr_work_create(matrix, work_mutable=.TRUE.)
nblkrows_tot = dbcsr_nblkrows_total(matrix)
nblkcols_tot = dbcsr_nblkcols_total(matrix)
DO row = 1, nblkrows_tot
DO col = 1, nblkcols_tot
tr = .FALSE.
iblock_row = row
iblock_col = col
CALL dbcsr_get_stored_coordinates(matrix, &
iblock_row, iblock_col, hold)
IF (hold .EQ. mynode) THEN
NULLIFY (p_new_block)
CALL dbcsr_reserve_block2d(matrix, &
iblock_row, iblock_col, p_new_block)
CPASSERT(ASSOCIATED(p_new_block))
p_new_block(:, :) = 1.0_dp
ENDIF
ENDDO
ENDDO
CALL dbcsr_finalize(matrix)
END SUBROUTINE fill_matrix_with_ones
! **************************************************************************************************
!> \brief builds projected KS matrices for the overlapping domains
!> also computes the DIIS error vector as a by-product
@ -2198,7 +2241,7 @@ CONTAINS
TYPE(domain_submatrix_type), DIMENSION(:), &
INTENT(IN), OPTIONAL :: subm_s_inv, subm_s_inv_half, &
subm_s_half, subm_r_down
TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: matrix_trimmer
TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_trimmer
TYPE(dbcsr_type), INTENT(IN) :: dpattern
TYPE(domain_map_type), INTENT(IN) :: map
INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain

1865
src/almo_scf_optimizer.F
View file

File diff suppressed because it is too large Load diff

15
src/almo_scf_types.F
View file

@ -46,8 +46,11 @@ MODULE almo_scf_types
! methods that add penalty terms to the energy functional
TYPE penalty_type
REAL(KIND=dp) :: occ_vol_coeff
INTEGER :: occ_vol_method
REAL(KIND=dp) :: final_determinant, penalty_strength, &
determinant_tolerance, penalty_strength_dec_factor, &
compactification_filter_start
INTEGER :: operator_type
LOGICAL :: virtual_nlmos
END TYPE penalty_type
@ -74,6 +77,7 @@ MODULE almo_scf_types
neglect_threshold
INTEGER :: optimizer_type ! diis, pcg, etc.
TYPE(penalty_type) :: opt_penalty
INTEGER :: preconditioner, & ! preconditioner type
conjugator, & ! conjugator type
@ -207,7 +211,7 @@ MODULE almo_scf_types
LOGICAL :: s_sqrt_done
REAL(KIND=dp) :: almo_scf_energy
LOGICAL :: orthogonal_basis, fixed_mu
LOGICAL :: return_orthogonalized_mos
LOGICAL :: return_orthogonalized_mos, construct_nlmos
!! Smearing control
!! smear flag allow to retrieve eigenvalues in almo_scf with diag algorithm and create occupation-scaled ALMO orbitals
@ -239,6 +243,9 @@ MODULE almo_scf_types
! mo overlap inversion algorithm
INTEGER :: sigma_inv_algorithm
! Determinant of the ALMO overlap matrix
REAL(KIND=dp) :: overlap_determinant
! ALMO SCF delocalization control
LOGICAL :: perturbative_delocalization
INTEGER :: quencher_radius_type
@ -354,7 +361,6 @@ MODULE almo_scf_types
! Options for various subsection options collected neatly
TYPE(almo_analysis_type) :: almo_analysis
TYPE(penalty_type) :: penalty
! Options for various optimizers collected neatly
TYPE(optimizer_options_type) :: opt_block_diag_diis
@ -362,6 +368,7 @@ MODULE almo_scf_types
TYPE(optimizer_options_type) :: opt_xalmo_diis
TYPE(optimizer_options_type) :: opt_xalmo_pcg
TYPE(optimizer_options_type) :: opt_xalmo_trustr
TYPE(optimizer_options_type) :: opt_nlmo_pcg
TYPE(optimizer_options_type) :: opt_xalmo_newton_pcg_solver
TYPE(optimizer_options_type) :: opt_k_pcg

5
src/input_constants.F
View file

@ -917,8 +917,9 @@ MODULE input_constants
almo_scf_trustr = 4, &
almo_scf_skip = 0
INTEGER, PARAMETER, PUBLIC :: almo_occ_vol_penalty_none = 0, &
almo_occ_vol_penalty_lndet = 1
INTEGER, PARAMETER, PUBLIC :: penalty_type_none = 0, &
penalty_type_lndet = 1, &
penalty_type_nlmo = 2
! optimizer parameters
INTEGER, PARAMETER, PUBLIC :: cg_zero = 0, &

89
src/input_cp2k_almo.F
View file

@ -17,13 +17,13 @@ MODULE input_cp2k_almo
USE input_constants, ONLY: &
almo_deloc_none, almo_deloc_scf, almo_deloc_x, almo_deloc_x_then_scf, &
almo_deloc_xalmo_1diag, almo_deloc_xalmo_scf, almo_deloc_xalmo_x, almo_frz_crystal, &
almo_frz_none, almo_occ_vol_penalty_lndet, almo_occ_vol_penalty_none, almo_scf_diag, &
almo_scf_pcg, almo_scf_skip, almo_scf_trustr, atomic_guess, cg_dai_yuan, cg_fletcher, &
cg_fletcher_reeves, cg_hager_zhang, cg_hestenes_stiefel, cg_liu_storey, cg_polak_ribiere, &
cg_zero, molecular_guess, optimizer_diis, optimizer_lin_eq_pcg, optimizer_pcg, &
optimizer_trustr, spd_inversion_dense_cholesky, spd_inversion_ls_hotelling, &
spd_inversion_ls_taylor, trustr_cauchy, trustr_dogleg, trustr_steihaug, xalmo_prec_domain, &
xalmo_prec_full, xalmo_prec_zero, xalmo_trial_r0_out, xalmo_trial_simplex
almo_frz_none, almo_scf_diag, almo_scf_pcg, almo_scf_skip, almo_scf_trustr, atomic_guess, &
cg_dai_yuan, cg_fletcher, cg_fletcher_reeves, cg_hager_zhang, cg_hestenes_stiefel, &
cg_liu_storey, cg_polak_ribiere, cg_zero, molecular_guess, op_loc_berry, op_loc_pipek, &
optimizer_diis, optimizer_lin_eq_pcg, optimizer_pcg, optimizer_trustr, &
spd_inversion_dense_cholesky, spd_inversion_ls_hotelling, spd_inversion_ls_taylor, &
trustr_cauchy, trustr_dogleg, trustr_steihaug, xalmo_prec_domain, xalmo_prec_full, &
xalmo_prec_zero, xalmo_trial_r0_out, xalmo_trial_simplex
USE input_keyword_types, ONLY: keyword_create,&
keyword_release,&
keyword_type
@ -46,6 +46,7 @@ MODULE input_cp2k_almo
INTEGER, PARAMETER, PRIVATE :: optimizer_xalmo_pcg = 3
INTEGER, PARAMETER, PRIVATE :: optimizer_xalmo_trustr = 4
INTEGER, PARAMETER, PRIVATE :: optimizer_newton_pcg_solver = 5
INTEGER, PARAMETER, PRIVATE :: optimizer_nlmo_pcg = 6
PUBLIC :: create_almo_scf_section
@ -210,6 +211,12 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="CONSTRUCT_NLMOS", &
description="Turns on post-SCF construction of NLMOs", &
usage="CONSTRUCT_NLMOS .TRUE.", default_l_val=.FALSE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
!CALL keyword_create(keyword, __LOCATION__, name="DOMAIN_LAYOUT_MOS",&
! description="Each electron in the system is constrained to its own delocalization domain."//&
! " This keyword creates groups of electrons that share the same domain.",&
@ -533,7 +540,7 @@ CONTAINS
CALL section_release(subsection)
NULLIFY (subsection)
CALL create_penalty_section(subsection)
CALL create_optimizer_section(subsection, optimizer_nlmo_pcg)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
@ -584,12 +591,21 @@ CONTAINS
CASE (optimizer_block_diagonal_pcg)
CALL section_create(section, __LOCATION__, name="ALMO_OPTIMIZER_PCG", &
description="Controls the PCG optimization of block-diagonal ALMOs.", &
n_keywords=9, n_subsections=0, repeats=.FALSE.)
n_keywords=9, n_subsections=1, repeats=.FALSE.)
optimizer_type = optimizer_pcg
CASE (optimizer_nlmo_pcg)
CALL section_create(section, __LOCATION__, name="NLMO_OPTIMIZER_PCG", &
description="Controls the PCG optimization of nonorthogonal localized MOs.", &
n_keywords=9, n_subsections=1, repeats=.FALSE.)
optimizer_type = optimizer_pcg
NULLIFY (subsection)
CALL create_penalty_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CASE (optimizer_xalmo_pcg)
CALL section_create(section, __LOCATION__, name="XALMO_OPTIMIZER_PCG", &
description="Controls the PCG optimization of extended ALMOs.", &
n_keywords=10, n_subsections=1, repeats=.FALSE.)
n_keywords=10, n_subsections=2, repeats=.FALSE.)
NULLIFY (subsection)
CALL create_optimizer_section(subsection, optimizer_newton_pcg_solver)
CALL section_add_subsection(section, subsection)
@ -848,25 +864,54 @@ CONTAINS
CALL section_create(section, __LOCATION__, name="PENALTY", &
description="Add penalty terms to the energy functional.", &
n_keywords=2, n_subsections=0, repeats=.FALSE.)
n_keywords=3, n_subsections=0, repeats=.FALSE.)
NULLIFY (keyword)
CALL keyword_create( &
keyword, __LOCATION__, name="OCCUPIED_VOLUME_PENALTY_METHOD", &
description="Penalty that prevents nonorthogonal orbitals from becoming linear dependent.", &
usage="OCCUPIED_VOLUME_PENALTY_METHOD LNDET", &
default_i_val=almo_occ_vol_penalty_none, &
enum_c_vals=s2a("NONE", "LNDET"), &
enum_desc=s2a("Do not use penalties", &
"Use -coeff*ln(det(MO-overlap)) term."), &
enum_i_vals=(/almo_occ_vol_penalty_none, almo_occ_vol_penalty_lndet/))
keyword, __LOCATION__, name="OPERATOR", &
description="Type of opertator which defines the spread functional", &
usage="OPERATOR PIPEK", &
enum_c_vals=s2a("BERRY", "PIPEK"), &
enum_i_vals=(/op_loc_berry, op_loc_pipek/), &
default_i_val=op_loc_berry)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="OCCUPIED_VOLUME_PENALTY_COEFF", &
description="Multiplication factor that determines the strength of the penalty term.", &
usage="OCCUPIED_VOLUME_PENALTY_COEFF 10.0", default_r_val=1.0_dp)
CALL keyword_create(keyword, __LOCATION__, name="PENALTY_STRENGTH", &
description="Strength of the orthogonalization penalty", &
usage="PENALTY_STRENGTH 1.1", default_r_val=1.1_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="PENALTY_STRENGTH_DECREASE_FACTOR", &
description="Factor that decreases the strength of the orthogonalization penalty.", &
usage="PENALTY_STRENGTH_DECREASE_FACTOR 1.1", default_r_val=1.1_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="DETERMINANT_TOLERANCE", &
description="Stop the optimization of the penalty strength if the determinant of the overlap "// &
"changes less than this tolerance threshold.", &
usage="DETERMINANT_TOLERANCE 1.0E-4", default_r_val=1.0E-3_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="FINAL_DETERMINANT", &
description="The final determinant that obtained after optimization.", &
usage="FINAL_DETERMINANT 0.1", default_r_val=0.1_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="COMPACTIFICATION_FILTER_START", &
description="Set orbital coefficients with absolute value smaller than this value to zero.", &
usage="COMPACTIFICATION_FILTER_START 1.e-6", default_r_val=-1.0_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="VIRTUAL_NLMOS", &
description="Localize virtual oribtals", &
usage="VIRTUAL_NLMOS .TRUE.", default_l_val=.FALSE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

9
src/iterate_matrix.F
View file

@ -28,6 +28,7 @@ MODULE iterate_matrix
m_walltime
USE mathconstants, ONLY: ifac
USE mathlib, ONLY: abnormal_value
USE message_passing, ONLY: mp_sum
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -70,8 +71,9 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'determinant', routineP = moduleN//':'//routineN
INTEGER :: handle, i, max_iter_lanczos, nsize, &
order_lanczos, sign_iter, unit_nr
INTEGER :: group, handle, i, max_iter_lanczos, &
nsize, order_lanczos, sign_iter, &
unit_nr
INTEGER(KIND=int_8) :: flop1
INTEGER, SAVE :: recursion_depth = 0
REAL(KIND=dp) :: det0, eps_lanczos, frobnorm, maxnorm, &
@ -102,9 +104,10 @@ CONTAINS
matrix_type=dbcsr_type_no_symmetry)
! compute the product of the diagonal elements
CALL dbcsr_get_info(matrix, nfullrows_total=nsize)
CALL dbcsr_get_info(matrix, nfullrows_total=nsize, group=group)
ALLOCATE (diagonal(nsize))
CALL dbcsr_get_diag(matrix, diagonal)
CALL mp_sum(diagonal, group)
det = PRODUCT(diagonal)
! create diagonal SQRTI matrix

11
src/qs_loc_methods.F
View file

@ -302,6 +302,7 @@ CONTAINS
NULLIFY (particle_set)
! get rows and cols of the input
CALL cp_fm_get_info(vectors, nrow_global=nao, ncol_global=nmoloc)
! replicate the input kind of matrix
CALL cp_fm_create(opvec, vectors%matrix_struct)
CALL cp_fm_set_all(opvec, 0.0_dp)
@ -329,18 +330,22 @@ CONTAINS
CALL cp_fm_get_info(zij_fm_set(iatom, 1)%matrix, ncol_global=ldz)
isgf = first_sgf(iatom)
ncol = nsgf(iatom)
! multiply fmxfm, using only part of the ao : Ct x S
CALL cp_gemm('N', 'N', nao, nmoloc, nao, 1.0_dp, ov_fm, vectors, 0.0_dp, opvec, &
a_first_col=1, a_first_row=1, b_first_col=1, b_first_row=1)
CALL cp_gemm('T', 'N', nmoloc, nmoloc, ncol, 0.5_dp, vectors, opvec, &
0.0_dp, zij_fm_set(iatom, 1)%matrix, &
a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf)
CALL cp_gemm('N', 'N', nao, nmoloc, ncol, 1.0_dp, ov_fm, vectors, 0.0_dp, opvec, &
a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf)
a_first_col=isgf, a_first_row=1, b_first_col=1, b_first_row=isgf)
CALL cp_gemm('T', 'N', nmoloc, nmoloc, nao, 0.5_dp, vectors, opvec, &
1.0_dp, zij_fm_set(iatom, 1)%matrix, &
a_first_col=1, a_first_row=1, b_first_col=1, b_first_row=1)
END DO ! iatom
! And now perform the optimization and rotate the orbitals
@ -646,9 +651,12 @@ CONTAINS
END IF
DO ispin = s_spin, l_spin
CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, nmo=nmo)
loc_method = localized_wfn_control%localization_method
SELECT CASE (localized_wfn_control%operator_type)
CASE (op_loc_berry)
! Here we allocate op_fm_set with the RIGHT size for uks
NULLIFY (tmp_fm_struct, mo_coeff)
@ -691,6 +699,7 @@ CONTAINS
CPABORT("Boys localization not implemented")
CASE (op_loc_pipek)
CALL optimize_loc_pipek(qs_env, loc_method, qs_loc_env, moloc_coeff(ispin)%matrix, &
op_fm_set, ispin, print_loc_section)

76
src/qs_loc_utils.F
View file

@ -99,7 +99,7 @@ MODULE qs_loc_utils
! *** Public ***
PUBLIC :: jacobi_rotation_pipek, qs_loc_env_init, loc_write_restart, &
retain_history, qs_loc_init, &
retain_history, qs_loc_init, compute_berry_operator, &
set_loc_centers, set_loc_wfn_lists, qs_loc_control_init
CONTAINS
@ -278,7 +278,9 @@ CONTAINS
ct = COS(theta)
st = SIN(theta)
CALL rotate_rmat_real(istate, jstate, st, ct, rmat)
CALL rotate_zij_real(istate, jstate, st, ct, zij_fm_set)
END DO
END DO
@ -649,21 +651,13 @@ CONTAINS
END SUBROUTINE qs_loc_env_init
! **************************************************************************************************
!> \brief Computes the Berry operator for periodic systems
!> used to define the spread of the MOS
!> Here the matrix elements of the type <mu|cos(kr)|nu> and <mu|sin(kr)|nu>
!> are computed, where mu and nu are the contracted basis functions.
!> Namely the Berry operator is exp(ikr)
!> k is defined somewhere
!> the pair lists are exploited and sparse matrixes are constructed
!> \brief A wrapper to compute the Berry operator for periodic systems
!> \param qs_loc_env new environment for the localization calculations
!> \param qs_env the qs_env in which the qs_env lives
!> \par History
!> 04.2005 created [MI]
!> 04.2018 modified [RZK, ZL]
!> \author MI
!> \note
!> The intgrals are computed analytically using the primitives GTO
!> The contraction is performed block-wise
! **************************************************************************************************
SUBROUTINE get_berry_operator(qs_loc_env, qs_env)
TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env
@ -672,9 +666,51 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'get_berry_operator', &
routineP = moduleN//':'//routineN
INTEGER :: dim_op, handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, &
last_jatom, ldab, ldsa, ldsb, ldwork, maxl, ncoa, ncob, nkind, nrow, nseta, nsetb, sgfa, &
sgfb
INTEGER :: dim_op, handle
TYPE(cell_type), POINTER :: cell
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
CALL timeset(routineN, handle)
NULLIFY (cell, op_sm_set)
CALL get_qs_loc_env(qs_loc_env=qs_loc_env, op_sm_set=op_sm_set, &
cell=cell, dim_op=dim_op)
CALL compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
CALL timestop(handle)
END SUBROUTINE get_berry_operator
! **************************************************************************************************
!> \brief Computes the Berry operator for periodic systems
!> used to define the spread of the MOS
!> Here the matrix elements of the type <mu|cos(kr)|nu> and <mu|sin(kr)|nu>
!> are computed, where mu and nu are the contracted basis functions.
!> Namely the Berry operator is exp(ikr)
!> k is defined somewhere
!> the pair lists are exploited and sparse matrixes are constructed
!> \param qs_env the qs_env in which the qs_env lives
!> \param cell ...
!> \param op_sm_set ...
!> \param dim_op ...
!> \par History
!> 04.2005 created [MI]
!> 04.2018 wrapped old code [RZK, ZL]
!> \author MI
!> \note
!> The intgrals are computed analytically using the primitives GTO
!> The contraction is performed block-wise
! **************************************************************************************************
SUBROUTINE compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cell_type), POINTER :: cell
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
INTEGER :: dim_op
CHARACTER(len=*), PARAMETER :: routineN = 'compute_berry_operator', &
routineP = moduleN//':'//routineN
INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
ldab, ldsa, ldsb, ldwork, maxl, ncoa, ncob, nkind, nrow, nseta, nsetb, sgfa, sgfb
INTEGER, DIMENSION(3) :: perd0
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
npgfb, nsgfa, nsgfb
@ -686,8 +722,6 @@ CONTAINS
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cosab, rpgfa, rpgfb, sinab, sphi_a, &
sphi_b, work, zeta, zetb
TYPE(block_p_type), DIMENSION(:), POINTER :: op_cos, op_sin
TYPE(cell_type), POINTER :: cell
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
TYPE(neighbor_list_iterator_p_type), &
@ -699,9 +733,8 @@ CONTAINS
TYPE(qs_kind_type), POINTER :: qs_kind
CALL timeset(routineN, handle)
NULLIFY (qs_kind, qs_kind_set)
NULLIFY (cell, op_sm_set, particle_set)
NULLIFY (particle_set)
NULLIFY (sab_orb)
NULLIFY (cosab, sinab, work)
NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb)
@ -722,9 +755,6 @@ CONTAINS
ALLOCATE (work(ldwork, ldwork))
work = 0.0_dp
CALL get_qs_loc_env(qs_loc_env=qs_loc_env, op_sm_set=op_sm_set, &
cell=cell, dim_op=dim_op)
ALLOCATE (op_cos(dim_op))
ALLOCATE (op_sin(dim_op))
DO i = 1, dim_op
@ -838,7 +868,6 @@ CONTAINS
! *** Calculate the primitive overlap integrals ***
DO i = 1, dim_op
kvec(1:3) = vector_k(1:3, i)
cosab = 0.0_dp
sinab = 0.0_dp
@ -846,7 +875,6 @@ CONTAINS
la_min(iset), lb_max(jset), npgfb(jset), zetb(:, jset), &
rpgfb(:, jset), lb_min(jset), &
ra, rb, kvec, cosab, sinab)
CALL contract_cossin(op_cos(i)%block, op_sin(i)%block, &
iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
@ -876,7 +904,7 @@ CONTAINS
DEALLOCATE (cosab, sinab, work, basis_set_list)
CALL timestop(handle)
END SUBROUTINE get_berry_operator
END SUBROUTINE compute_berry_operator
! **************************************************************************************************
!> \brief ...

93
tests/QS/regtest-nlmo/Si-nlmos.inp Normal file
View file

@ -0,0 +1,93 @@
&GLOBAL
PROJECT NLMOS
RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME GTH_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 200
NGRIDS 4
&END MGRID
&QS
ALMO_SCF T
EPS_DEFAULT 1.0E-08
&END QS
&ALMO_SCF
EPS_FILTER 1.0E-09
ALMO_ALGORITHM SKIP
MO_OVERLAP_INV_ALG DENSE_CHOLESKY
DELOCALIZE_METHOD FULL_SCF
ALMO_SCF_GUESS ATOMIC
XALMO_TRIAL_WF SIMPLE
CONSTRUCT_NLMOS TRUE
&XALMO_OPTIMIZER_PCG
MAX_ITER 50
EPS_ERROR 1.0E-3
CONJUGATOR HESTENES_STIEFEL
PRECONDITIONER DEFAULT
LIN_SEARCH_EPS_ERROR 0.1
LIN_SEARCH_STEP_SIZE_GUESS 0.2
MAX_ITER_OUTER_LOOP 10
&END XALMO_OPTIMIZER_PCG
&NLMO_OPTIMIZER_PCG
MAX_ITER 200
EPS_ERROR 0.1
CONJUGATOR ZERO
PRECONDITIONER DEFAULT
LIN_SEARCH_EPS_ERROR 0.5
LIN_SEARCH_STEP_SIZE_GUESS 0.001
MAX_ITER_OUTER_LOOP 0
&PENALTY
OPERATOR PIPEK
PENALTY_STRENGTH 0.01
DETERMINANT_TOLERANCE 1.0E-8
PENALTY_STRENGTH_DECREASE_FACTOR 2.0
FINAL_DETERMINANT 0.6
COMPACTIFICATION_FILTER_START 1.0E-2
VIRTUAL_NLMOS FALSE
&END PENALTY
&END NLMO_OPTIMIZER_PCG
&END ALMO_SCF
&XC
&XC_FUNCTIONAL BLYP
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 5.430710 5.430710 5.430710
&END CELL
&COORD
SCALED T
Si 0.000000 0.000000 0.000000
Si 0.000000 0.500000 0.500000
Si 0.500000 0.000000 0.500000
Si 0.500000 0.500000 0.000000
Si 0.250000 0.250000 0.250000
Si 0.250000 0.750000 0.750000
Si 0.750000 0.250000 0.750000
Si 0.750000 0.750000 0.250000
&END COORD
&KIND Si
BASIS_SET SZV-GTH
POTENTIAL GTH-PBE-q4
&END KIND
&END SUBSYS
&END FORCE_EVAL

3
tests/QS/regtest-nlmo/TEST_FILES Normal file
View file

@ -0,0 +1,3 @@
pipek_C6H6.inp 90 1e-04 177.2020586640
Si-nlmos.inp 90 1e-04 118.5296581594
#EOF

102
tests/QS/regtest-nlmo/pipek_C6H6.inp Normal file
View file

@ -0,0 +1,102 @@
&GLOBAL
PROJECT C6H6
RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME GTH_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 200
NGRIDS 4
&END MGRID
&QS
ALMO_SCF T
EPS_DEFAULT 1.0E-8
&END QS
&ALMO_SCF
EPS_FILTER 1.0E-09
ALMO_ALGORITHM SKIP
MO_OVERLAP_INV_ALG DENSE_CHOLESKY
DELOCALIZE_METHOD FULL_SCF
ALMO_SCF_GUESS ATOMIC
XALMO_TRIAL_WF SIMPLE
CONSTRUCT_NLMOS TRUE
RETURN_ORTHOGONALIZED_MOS FALSE
&XALMO_OPTIMIZER_PCG
MAX_ITER 50
EPS_ERROR 1.0E-3
CONJUGATOR HESTENES_STIEFEL
PRECONDITIONER DEFAULT
LIN_SEARCH_EPS_ERROR 0.01
LIN_SEARCH_STEP_SIZE_GUESS 0.2
MAX_ITER_OUTER_LOOP 2
&END XALMO_OPTIMIZER_PCG
&NLMO_OPTIMIZER_PCG
MAX_ITER 20000
EPS_ERROR 0.2
CONJUGATOR ZERO
PRECONDITIONER DEFAULT
LIN_SEARCH_EPS_ERROR 0.1
LIN_SEARCH_STEP_SIZE_GUESS 0.0001
MAX_ITER_OUTER_LOOP 5
&PENALTY
OPERATOR PIPEK
PENALTY_STRENGTH 0.01
DETERMINANT_TOLERANCE 1.0E-8
PENALTY_STRENGTH_DECREASE_FACTOR 2.0
FINAL_DETERMINANT 0.5
!COMPACTIFICATION_FILTER_START 1.0E-3
VIRTUAL_NLMOS FALSE
&END PENALTY
&END NLMO_OPTIMIZER_PCG
&END ALMO_SCF
&XC
&XC_FUNCTIONAL BLYP
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 9.0 9.0 5.0
&END CELL
&TOPOLOGY
&CENTER_COORDINATES
&END
&GENERATE
CREATE_MOLECULES TRUE
&END
&END
&COORD
C 3.5560000000 4.5610000000 0.0000000000
C 4.5060000000 3.5350000000 0.0000000000
C 5.8680000000 3.8450000000 0.0000000000
C 6.2820000000 5.1790000000 0.0000000000
C 5.3320000000 6.2050000000 0.0000000000
C 3.9700000000 5.8950000000 0.0000000000
H 2.5000000000 4.3210000000 0.0000000000
H 4.1850000000 2.5000000000 0.0000000000
H 6.6040000000 3.0490000000 0.0000000000
H 7.3390000000 5.4190000000 0.0000000000
H 5.6530000000 7.2400000000 0.0000000000
H 3.2340000000 6.6910000000 0.0000000000
&END COORD
&KIND C
BASIS_SET SZV-GTH
POTENTIAL GTH-BLYP-q4
&END KIND
&KIND H
BASIS_SET SZV-GTH
POTENTIAL GTH-BLYP-q1
&END KIND
&END SUBSYS
&END FORCE_EVAL

1
tests/TEST_DIRS
View file

@ -263,3 +263,4 @@ QS/regtest-sto
QS/regtest-ri-laplace-mp2-cubic libint mpiranks==2||mpiranks==4||mpiranks==6||mpiranks==16||mpiranks==24
QS/regtest-dm-ls-scf-4
QS/regtest-almo-trustr
QS/regtest-nlmo

3
tests/TEST_TYPES
View file

@ -1,4 +1,4 @@
89
90
Total energy:!3
POTENTIAL ENERGY!4
Total energy \[eV\]:!4
@ -88,6 +88,7 @@ DIPOLE : CheckSum =!5
POLAR : CheckSum =!5
XAS excitation energy (eV): !7
Electronic density on regular grids:!7
Final localization: !3
#
# 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