diff --git a/src/almo_scf.F b/src/almo_scf.F index 41dd342a71..f39651c0db 100644 --- a/src/almo_scf.F +++ b/src/almo_scf.F @@ -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 diff --git a/src/dm_ls_scf_qs.F b/src/dm_ls_scf_qs.F index 364c176926..3d46f23049 100644 --- a/src/dm_ls_scf_qs.F +++ b/src/dm_ls_scf_qs.F @@ -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,& diff --git a/src/energy_corrections.F b/src/energy_corrections.F index f19963282a..796aec36e2 100644 --- a/src/energy_corrections.F +++ b/src/energy_corrections.F @@ -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,& diff --git a/src/fm/cp_fm_basic_linalg.F b/src/fm/cp_fm_basic_linalg.F index fa6b17b993..b17045ef06 100644 --- a/src/fm/cp_fm_basic_linalg.F +++ b/src/fm/cp_fm_basic_linalg.F @@ -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 diff --git a/src/input_constants.F b/src/input_constants.F index cb83417f97..e54a035db4 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -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, & diff --git a/src/input_cp2k_loc.F b/src/input_cp2k_loc.F index d34d7825a5..80176bcec3 100644 --- a/src/input_cp2k_loc.F +++ b/src/input_cp2k_loc.F @@ -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 ", & + 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) diff --git a/src/qs_atomic_block.F b/src/qs_atomic_block.F new file mode 100644 index 0000000000..1a9f5b6de3 --- /dev/null +++ b/src/qs_atomic_block.F @@ -0,0 +1,186 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2022 CP2K developers group ! +! ! +! 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 diff --git a/src/qs_initial_guess.F b/src/qs_initial_guess.F index 89f3bdeab5..678155b0fe 100644 --- a/src/qs_initial_guess.F +++ b/src/qs_initial_guess.F @@ -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 ... diff --git a/src/qs_loc_methods.F b/src/qs_loc_methods.F index 3ed6d4c8a4..ff1538f5e8 100644 --- a/src/qs_loc_methods.F +++ b/src/qs_loc_methods.F @@ -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 diff --git a/src/qs_loc_states.F b/src/qs_loc_states.F index 6bf73c91e9..30d3dd46e4 100644 --- a/src/qs_loc_states.F +++ b/src/qs_loc_states.F @@ -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 diff --git a/src/qs_loc_types.F b/src/qs_loc_types.F index 2f13aac7ba..f88049b15b 100644 --- a/src/qs_loc_types.F +++ b/src/qs_loc_types.F @@ -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. diff --git a/src/qs_loc_utils.F b/src/qs_loc_utils.F index c56fd015ea..6cad2418f4 100644 --- a/src/qs_loc_utils.F +++ b/src/qs_loc_utils.F @@ -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) diff --git a/src/qs_localization_methods.F b/src/qs_localization_methods.F index 52000ad3ca..47a9a4f06d 100644 --- a/src/qs_localization_methods.F +++ b/src/qs_localization_methods.F @@ -20,16 +20,15 @@ MODULE qs_localization_methods USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,& cp_cfm_gemm,& - cp_cfm_schur_product + cp_cfm_scale,& + cp_cfm_scale_and_add,& + cp_cfm_schur_product,& + cp_cfm_trace USE cp_cfm_diag, ONLY: cp_cfm_heevd - USE cp_cfm_types, ONLY: cp_cfm_create,& - cp_cfm_get_element,& - cp_cfm_get_info,& - cp_cfm_p_type,& - cp_cfm_release,& - cp_cfm_set_all,& - cp_cfm_to_cfm,& - cp_cfm_type + USE cp_cfm_types, ONLY: & + cp_cfm_create, cp_cfm_get_element, cp_cfm_get_info, cp_cfm_get_submatrix, cp_cfm_p_type, & + cp_cfm_release, cp_cfm_set_all, cp_cfm_set_submatrix, cp_cfm_to_cfm, cp_cfm_to_fm, & + cp_cfm_type, cp_fm_to_cfm USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply USE cp_external_control, ONLY: external_control USE cp_fm_basic_linalg, ONLY: cp_fm_frobenius_norm,& @@ -45,14 +44,15 @@ MODULE qs_localization_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_maxabsrownorm, & - cp_fm_maxabsval, cp_fm_p_type, cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, & - cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type + cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_init_random, & + cp_fm_maxabsrownorm, cp_fm_maxabsval, cp_fm_p_type, cp_fm_release, cp_fm_set_all, & + cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type USE cp_gemm_interface, ONLY: cp_gemm USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,& cp_logger_get_default_unit_nr USE cp_para_types, ONLY: cp_para_env_type USE dbcsr_api, ONLY: dbcsr_p_type + USE kahan_sum, ONLY: accurate_sum USE kinds, ONLY: dp USE machine, ONLY: m_flush,& m_walltime @@ -66,11 +66,13 @@ MODULE qs_localization_methods mp_sendrecv,& mp_sum,& mp_sync + USE qs_mo_methods, ONLY: make_basis_simple #include "./base/base_uses.f90" IMPLICIT NONE PUBLIC :: initialize_weights, crazy_rotations, & - direct_mini, rotate_orbitals, approx_l1_norm_sd, jacobi_rotations, scdm_qrfact, zij_matrix + direct_mini, rotate_orbitals, approx_l1_norm_sd, jacobi_rotations, scdm_qrfact, zij_matrix, & + jacobi_cg_edf_ls CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_localization_methods' @@ -265,7 +267,7 @@ CONTAINS REAL(KIND=dp), DIMENSION(3, 3) :: metric metric = 0.0_dp - CALL dgemm('T', 'N', 3, 3, 3, 1._dp, cell%hmat, 3, cell%hmat, 3, 0.0_dp, metric, 3) + CALL dgemm('T', 'N', 3, 3, 3, 1._dp, cell%hmat(:, :), 3, cell%hmat(:, :), 3, 0.0_dp, metric(:, :), 3) weights(1) = METRIC(1, 1) - METRIC(1, 2) - METRIC(1, 3) weights(2) = METRIC(2, 2) - METRIC(1, 2) - METRIC(2, 3) @@ -292,8 +294,8 @@ CONTAINS !> \par History !> \author Joost VandeVondele (02.2010) ! ************************************************************************************************** - SUBROUTINE jacobi_rotations(weights, zij, vectors, para_env, max_iter, eps_localization, & - sweeps, out_each, target_time, start_time) + SUBROUTINE jacobi_rotations(weights, zij, vectors, para_env, max_iter, & + eps_localization, sweeps, out_each, target_time, start_time) REAL(KIND=dp), INTENT(IN) :: weights(:) TYPE(cp_fm_p_type), INTENT(INOUT) :: ZIJ(:, :) @@ -325,7 +327,8 @@ CONTAINS !> \param sweeps ... !> \param out_each ... ! ************************************************************************************************** - SUBROUTINE jacobi_rotations_serial(weights, zij, vectors, max_iter, eps_localization, sweeps, out_each) + SUBROUTINE jacobi_rotations_serial(weights, zij, vectors, max_iter, eps_localization, sweeps, & + out_each) REAL(KIND=dp), INTENT(IN) :: weights(:) TYPE(cp_fm_p_type), INTENT(INOUT) :: ZIJ(:, :) TYPE(cp_fm_type), POINTER :: vectors @@ -366,16 +369,19 @@ CONTAINS CALL cp_fm_get_info(rmat, nrow_global=nstate) tolerance = 1.0e10_dp + sweeps = 0 unit_nr = -1 IF (rmat%matrix_struct%para_env%mepos .EQ. rmat%matrix_struct%para_env%source) THEN unit_nr = cp_logger_get_default_unit_nr() WRITE (unit_nr, '(T4,A )') " Localization by iterative Jacobi rotation" END IF + ! do jacobi sweeps until converged DO WHILE (tolerance >= eps_localization .AND. sweeps < max_iter) sweeps = sweeps + 1 t1 = m_walltime() + DO istate = 1, nstate DO jstate = istate + 1, nstate DO idim = 1, dim2 @@ -387,16 +393,20 @@ CONTAINS st = SIN(theta) ct = COS(theta) CALL rotate_zij(istate, jstate, st, ct, c_zij) + CALL rotate_rmat(istate, jstate, st, ct, c_rmat) END DO END DO + CALL check_tolerance(c_zij, weights, tolerance) + t2 = m_walltime() IF (unit_nr > 0 .AND. MODULO(sweeps, out_each) == 0) THEN WRITE (unit_nr, '(T4,A,I7,T30,A,E12.4,T60,A,F8.3)') & "Iteration:", sweeps, "Tolerance:", tolerance, "Time:", t2 - t1 CALL m_flush(unit_nr) END IF + END DO DO idim = 1, dim2 @@ -415,6 +425,823 @@ CONTAINS END SUBROUTINE jacobi_rotations_serial ! ************************************************************************************************** +!> \brief very similar to jacobi_rotations_serial with some extra output options +!> \param weights ... +!> \param c_zij ... +!> \param max_iter ... +!> \param c_rmat ... +!> \param eps_localization ... +!> \param tol_out ... +!> \param jsweeps ... +!> \param out_each ... +!> \param c_zij_out ... +!> \param grad_final ... +! ************************************************************************************************** + SUBROUTINE jacobi_rotations_serial_1(weights, c_zij, max_iter, c_rmat, eps_localization, & + tol_out, jsweeps, out_each, c_zij_out, grad_final) + REAL(KIND=dp), INTENT(IN) :: weights(:) + TYPE(cp_cfm_p_type), INTENT(IN), POINTER :: c_zij(:) + INTEGER, INTENT(IN) :: max_iter + TYPE(cp_cfm_type), INTENT(INOUT), POINTER :: c_rmat + REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_localization + REAL(KIND=dp), INTENT(OUT), OPTIONAL :: tol_out + INTEGER, INTENT(OUT), OPTIONAL :: jsweeps + INTEGER, INTENT(IN), OPTIONAL :: out_each + TYPE(cp_cfm_p_type), INTENT(OUT), OPTIONAL, & + POINTER :: c_zij_out(:) + TYPE(cp_fm_type), INTENT(OUT), OPTIONAL, POINTER :: grad_final + + CHARACTER(len=*), PARAMETER :: routineN = 'jacobi_rotations_serial_1', & + routineP = moduleN//':'//routineN + + COMPLEX(KIND=dp) :: mzii + COMPLEX(KIND=dp), POINTER :: mii(:), mij(:), mjj(:) + INTEGER :: dim2, handle, idim, istate, jstate, & + nstate, sweeps, unit_nr + REAL(KIND=dp) :: alpha, avg_spread_ii, ct, spread_ii, st, & + sum_spread_ii, t1, t2, theta, tolerance + TYPE(cp_cfm_p_type), POINTER :: c_zij_local(:) + TYPE(cp_cfm_type), POINTER :: c_rmat_local + + CALL timeset(routineN, handle) + + dim2 = SIZE(c_zij) + NULLIFY (c_zij_local, c_rmat_local) + NULLIFY (mii, mij, mjj) + ALLOCATE (mii(dim2), mij(dim2), mjj(dim2)) + + ALLOCATE (c_zij_local(dim2)) + CALL cp_cfm_create(c_rmat_local, c_rmat%matrix_struct) + CALL cp_cfm_set_all(c_rmat_local, (0.0_dp, 0.0_dp), (1.0_dp, 0.0_dp)) + DO idim = 1, dim2 + NULLIFY (c_zij_local(idim)%matrix) + CALL cp_cfm_create(c_zij_local(idim)%matrix, c_zij(idim)%matrix%matrix_struct) + c_zij_local(idim)%matrix%local_data = c_zij(idim)%matrix%local_data + END DO + + CALL cp_cfm_get_info(c_rmat_local, nrow_global=nstate) + tolerance = 1.0e10_dp + + IF (PRESENT(grad_final)) CALL cp_fm_set_all(grad_final, 0.0_dp) + + sweeps = 0 + IF (PRESENT(out_each)) THEN + unit_nr = -1 + IF (c_rmat_local%matrix_struct%para_env%mepos .EQ. c_rmat_local%matrix_struct%para_env%source) THEN + unit_nr = cp_logger_get_default_unit_nr() + END IF + alpha = 0.0_dp + DO idim = 1, dim2 + alpha = alpha + weights(idim) + END DO + END IF + + ! do jacobi sweeps until converged + DO WHILE (sweeps < max_iter) + sweeps = sweeps + 1 + IF (PRESENT(eps_localization)) THEN + IF (tolerance < eps_localization) EXIT + END IF + IF (PRESENT(out_each)) t1 = m_walltime() + + DO istate = 1, nstate + DO jstate = istate + 1, nstate + DO idim = 1, dim2 + CALL cp_cfm_get_element(c_zij_local(idim)%matrix, istate, istate, mii(idim)) + CALL cp_cfm_get_element(c_zij_local(idim)%matrix, istate, jstate, mij(idim)) + CALL cp_cfm_get_element(c_zij_local(idim)%matrix, jstate, jstate, mjj(idim)) + END DO + CALL get_angle(mii, mjj, mij, weights, theta) + st = SIN(theta) + ct = COS(theta) + CALL rotate_zij(istate, jstate, st, ct, c_zij_local) + + CALL rotate_rmat(istate, jstate, st, ct, c_rmat_local) + END DO + END DO + + IF (PRESENT(grad_final)) THEN + CALL check_tolerance(c_zij_local, weights, tolerance, grad=grad_final) + ELSE + CALL check_tolerance(c_zij_local, weights, tolerance) + END IF + IF (PRESENT(tol_out)) tol_out = tolerance + + IF (PRESENT(out_each)) THEN + t2 = m_walltime() + IF (unit_nr > 0 .AND. MODULO(sweeps, out_each) == 0) THEN + sum_spread_ii = 0.0_dp + DO istate = 1, nstate + spread_ii = 0.0_dp + DO idim = 1, dim2 + CALL cp_cfm_get_element(c_zij_local(idim)%matrix, istate, istate, mzii) + spread_ii = spread_ii + weights(idim)* & + ABS(mzii)**2/twopi/twopi + END DO + sum_spread_ii = sum_spread_ii + spread_ii + END DO + sum_spread_ii = alpha*nstate/twopi/twopi - sum_spread_ii + avg_spread_ii = sum_spread_ii/nstate + WRITE (unit_nr, '(T4,A,T26,A,T48,A,T64,A)') & + "Iteration", "Avg. Spread_ii", "Tolerance", "Time" + WRITE (unit_nr, '(T4,I7,T20,F20.10,T45,E12.4,T60,F8.3)') & + sweeps, avg_spread_ii, tolerance, t2 - t1 + CALL m_flush(unit_nr) + END IF + IF (PRESENT(jsweeps)) jsweeps = sweeps + END IF + + END DO + + IF (PRESENT(c_zij_out)) THEN + DO idim = 1, dim2 + CALL cp_cfm_to_cfm(c_zij_local(idim)%matrix, c_zij_out(idim)%matrix) + END DO + END IF + CALL cp_cfm_to_cfm(c_rmat_local, c_rmat) + + DEALLOCATE (mii, mij, mjj) + DO idim = 1, dim2 + CALL cp_cfm_release(c_zij_local(idim)%matrix) + END DO + DEALLOCATE (c_zij_local) + CALL cp_cfm_release(c_rmat_local) + + CALL timestop(handle) + + END SUBROUTINE jacobi_rotations_serial_1 +! ************************************************************************************************** +!> \brief combine jacobi rotations (serial) and conjugate gradient with golden section line search +!> for partially occupied wannier functions +!> \param para_env ... +!> \param weights ... +!> \param zij ... +!> \param vectors ... +!> \param max_iter ... +!> \param eps_localization ... +!> \param iter ... +!> \param out_each ... +!> \param nextra ... +!> \param do_cg ... +!> \param nmo ... +!> \param vectors_2 ... +!> \param mos_guess ... +! ************************************************************************************************** + SUBROUTINE jacobi_cg_edf_ls(para_env, weights, zij, vectors, max_iter, eps_localization, & + iter, out_each, nextra, do_cg, nmo, vectors_2, mos_guess) + TYPE(cp_para_env_type), POINTER :: para_env + REAL(KIND=dp), INTENT(IN) :: weights(:) + TYPE(cp_fm_p_type), INTENT(INOUT) :: zij(:, :) + TYPE(cp_fm_type), POINTER :: vectors + INTEGER, INTENT(IN) :: max_iter + REAL(KIND=dp), INTENT(IN) :: eps_localization + INTEGER :: iter + INTEGER, INTENT(IN) :: out_each, nextra + LOGICAL, INTENT(IN) :: do_cg + INTEGER, INTENT(IN), OPTIONAL :: nmo + TYPE(cp_fm_type), INTENT(IN), OPTIONAL, POINTER :: vectors_2, mos_guess + + CHARACTER(len=*), PARAMETER :: routineN = 'jacobi_cg_edf_ls', & + routineP = moduleN//':'//routineN + COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), & + czero = (0.0_dp, 0.0_dp) + REAL(KIND=dp), PARAMETER :: gold_sec = 0.3819_dp + + COMPLEX(KIND=dp) :: cnorm2_Gct, cnorm2_Gct_cross, mzii + COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_cmat + COMPLEX(KIND=dp), DIMENSION(:), POINTER :: arr_zii + COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: matrix_zii + INTEGER :: dim2, handle, icinit, idim, istate, line_search_count, line_searches, lsl, lsm, & + lsr, miniter, nao, ndummy, nocc, norextra, northo, nstate, unit_nr + INTEGER, DIMENSION(1) :: iloc + LOGICAL :: do_cinit_mo, do_cinit_random, & + do_U_guess_mo, new_direction + REAL(KIND=dp) :: alpha, avg_spread_ii, beta, beta_pr, ds, ds_min, mintol, norm, norm2_Gct, & + norm2_Gct_cross, norm2_old, spread_ii, spread_sum, sum_spread_ii, t1, tol, tolc, weight + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sum_spread + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_mat, tmp_mat_1 + REAL(KIND=dp), DIMENSION(50) :: energy, pos + REAL(KIND=dp), DIMENSION(:), POINTER :: tmp_arr + TYPE(cp_blacs_env_type), POINTER :: context + TYPE(cp_cfm_p_type), DIMENSION(:), POINTER :: c_zij, zij_0 + TYPE(cp_cfm_type), POINTER :: c_tilde, ctrans_lambda, Gct_old, & + grad_ctilde, skc, tmp_cfm, tmp_cfm_1, & + tmp_cfm_2, U, UL, V, VL, zdiag + TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct + TYPE(cp_fm_type), POINTER :: id_nextra, matrix_U, matrix_V, & + matrix_V_all, rmat, tmp_fm, vectors_all + + CALL timeset(routineN, handle) + + dim2 = SIZE(zij, 2) + NULLIFY (context) + NULLIFY (rmat, c_zij, matrix_zii, arr_zii) + NULLIFY (U, UL, matrix_U) + NULLIFY (tmp_fm_struct) + NULLIFY (tmp_fm) + NULLIFY (tmp_cfm, tmp_cfm_1, tmp_cfm_2) + NULLIFY (tmp_arr) + + ALLOCATE (c_zij(dim2)) + + CALL cp_fm_get_info(zij(1, 1)%matrix, nrow_global=nstate) + + ALLOCATE (sum_spread(nstate)) + ALLOCATE (matrix_zii(nstate, dim2)) + matrix_zii = czero + sum_spread = 0.0_dp + + alpha = 0.0_dp + DO idim = 1, dim2 + alpha = alpha + weights(idim) + NULLIFY (c_zij(idim)%matrix) + CALL cp_cfm_create(c_zij(idim)%matrix, zij(1, 1)%matrix%matrix_struct) + c_zij(idim)%matrix%local_data = CMPLX(zij(1, idim)%matrix%local_data, & + zij(2, idim)%matrix%local_data, dp) + END DO + + NULLIFY (zij_0) + ALLOCATE (zij_0(dim2)) + + CALL cp_cfm_create(U, zij(1, 1)%matrix%matrix_struct) + CALL cp_fm_create(matrix_U, zij(1, 1)%matrix%matrix_struct) + + CALL cp_cfm_set_all(U, czero, cone) + CALL cp_fm_set_all(matrix_U, 0.0_dp, 1.0_dp) + + CALL cp_fm_get_info(vectors, nrow_global=nao) + IF (nextra > 0) THEN + IF (PRESENT(mos_guess)) THEN + do_cinit_random = .FALSE. + do_cinit_mo = .TRUE. + CALL cp_fm_get_info(mos_guess, ncol_global=ndummy) + ELSE + do_cinit_random = .TRUE. + do_cinit_mo = .FALSE. + ndummy = nstate + END IF + + IF (do_cinit_random) THEN + icinit = 1 + do_U_guess_mo = .FALSE. + ELSEIF (do_cinit_mo) THEN + icinit = 2 + do_U_guess_mo = .TRUE. + END IF + + nocc = nstate - nextra + northo = nmo - nocc + norextra = nmo - nstate + CALL cp_fm_struct_get(zij(1, 1)%matrix%matrix_struct, context=context) + + NULLIFY (c_tilde, V, VL, grad_ctilde, id_nextra, zdiag, tmp_cfm, & + tmp_cfm_1, ctrans_lambda, tmp_cfm_2, skc, Gct_old, & + vectors_all) + NULLIFY (matrix_V, matrix_V_all) + + ALLOCATE (tmp_cmat(nstate, nstate)) + CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=nmo, & + para_env=para_env, context=context) + DO idim = 1, dim2 + NULLIFY (zij_0(idim)%matrix) + CALL cp_cfm_create(zij_0(idim)%matrix, tmp_fm_struct) + CALL cp_cfm_set_all(zij_0(idim)%matrix, czero, cone) + CALL cp_cfm_get_submatrix(c_zij(idim)%matrix, tmp_cmat) + CALL cp_cfm_set_submatrix(zij_0(idim)%matrix, tmp_cmat) + END DO + CALL cp_fm_struct_release(tmp_fm_struct) + DEALLOCATE (tmp_cmat) + + CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=nstate, & + para_env=para_env, context=context) + CALL cp_cfm_create(V, tmp_fm_struct) + CALL cp_fm_create(matrix_V, tmp_fm_struct) + CALL cp_cfm_create(zdiag, tmp_fm_struct) + CALL cp_fm_create(rmat, tmp_fm_struct) + CALL cp_fm_struct_release(tmp_fm_struct) + CALL cp_cfm_set_all(V, czero, cone) + CALL cp_fm_set_all(matrix_V, 0.0_dp, 1.0_dp) + + CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=ndummy, & + para_env=para_env, context=context) + CALL cp_fm_create(matrix_V_all, tmp_fm_struct) + CALL cp_fm_struct_release(tmp_fm_struct) + CALL cp_fm_set_all(matrix_V_all, 0._dp, 1._dp) + + ALLOCATE (arr_zii(nstate)) + + CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=northo, ncol_global=nextra, & + para_env=para_env, context=context) + CALL cp_cfm_create(c_tilde, tmp_fm_struct) + CALL cp_cfm_create(grad_ctilde, tmp_fm_struct) + CALL cp_cfm_create(Gct_old, tmp_fm_struct) + CALL cp_cfm_create(skc, tmp_fm_struct) + CALL cp_fm_struct_release(tmp_fm_struct) + CALL cp_cfm_set_all(c_tilde, czero) + CALL cp_cfm_set_all(Gct_old, czero) + CALL cp_cfm_set_all(skc, czero) + + CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=northo, ncol_global=nstate, & + para_env=para_env, context=context) + CALL cp_cfm_create(VL, tmp_fm_struct) + CALL cp_cfm_set_all(VL, czero) + CALL cp_fm_struct_release(tmp_fm_struct) + + CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nextra, ncol_global=nextra, & + para_env=para_env, context=context) + CALL cp_fm_create(id_nextra, tmp_fm_struct) + CALL cp_cfm_create(ctrans_lambda, tmp_fm_struct) + CALL cp_fm_struct_release(tmp_fm_struct) + CALL cp_cfm_set_all(ctrans_lambda, czero) + CALL cp_fm_set_all(id_nextra, 0.0_dp, 1.0_dp) + + CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nextra, ncol_global=nstate, & + para_env=para_env, context=context) + CALL cp_cfm_create(UL, tmp_fm_struct) + CALL cp_fm_struct_release(tmp_fm_struct) + CALL cp_cfm_set_all(UL, czero) + + CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, ncol_global=nmo, & + para_env=para_env, context=context) + CALL cp_fm_create(vectors_all, tmp_fm_struct) + CALL cp_fm_struct_release(tmp_fm_struct) + ALLOCATE (tmp_mat(nao, nstate)) + CALL cp_fm_get_submatrix(vectors, tmp_mat) + CALL cp_fm_set_submatrix(vectors_all, tmp_mat, 1, 1, nao, nstate) + DEALLOCATE (tmp_mat) + ALLOCATE (tmp_mat(nao, norextra)) + CALL cp_fm_get_submatrix(vectors_2, tmp_mat) + CALL cp_fm_set_submatrix(vectors_all, tmp_mat, 1, nstate + 1, nao, norextra) + DEALLOCATE (tmp_mat) + + ! initialize c_tilde + SELECT CASE (icinit) + CASE (1) ! random coefficients + WRITE (*, *) "RANDOM INITIAL GUESS FOR C" + CALL cp_fm_create(tmp_fm, c_tilde%matrix_struct) + CALL cp_fm_init_random(tmp_fm, nextra) + CALL make_basis_simple(tmp_fm, nextra) + c_tilde%local_data = tmp_fm%local_data + CALL cp_fm_release(tmp_fm) + ALLOCATE (tmp_cmat(northo, nextra)) + CALL cp_cfm_get_submatrix(c_tilde, tmp_cmat) + CALL cp_cfm_set_submatrix(V, tmp_cmat, nocc + 1, nocc + 1, northo, nextra) + DEALLOCATE (tmp_cmat) + CASE (2) ! MO based coeffs + CALL cp_gemm("T", "N", nmo, ndummy, nao, 1.0_dp, vectors_all, mos_guess, 0.0_dp, matrix_V_all) + ALLOCATE (tmp_arr(nmo)) + ALLOCATE (tmp_mat(nmo, ndummy)) + ALLOCATE (tmp_mat_1(nmo, nstate)) + ! normalize matrix_V_all + CALL cp_fm_get_submatrix(matrix_V_all, tmp_mat) + DO istate = 1, ndummy + tmp_arr(:) = tmp_mat(:, istate) + norm = norm2(tmp_arr) + tmp_arr(:) = tmp_arr(:)/norm + tmp_mat(:, istate) = tmp_arr(:) + END DO + CALL cp_fm_set_submatrix(matrix_V_all, tmp_mat) + CALL cp_fm_get_submatrix(matrix_V_all, tmp_mat_1, 1, 1, nmo, nstate) + CALL cp_fm_set_submatrix(matrix_V, tmp_mat_1) + DEALLOCATE (tmp_arr, tmp_mat, tmp_mat_1) + CALL cp_fm_to_cfm(msourcer=matrix_V, mtarget=V) + ALLOCATE (tmp_mat(northo, ndummy)) + ALLOCATE (tmp_mat_1(northo, nextra)) + CALL cp_fm_get_submatrix(matrix_V_all, tmp_mat, nocc + 1, 1, northo, ndummy) + ALLOCATE (tmp_arr(ndummy)) + tmp_arr = 0.0_dp + DO istate = 1, ndummy + tmp_arr(istate) = norm2(tmp_mat(:, istate)) + END DO + ! find edfs + DO istate = 1, nextra + iloc = MAXLOC(tmp_arr) + tmp_mat_1(:, istate) = tmp_mat(:, iloc(1)) + tmp_arr(iloc(1)) = 0.0_dp + END DO + + DEALLOCATE (tmp_arr, tmp_mat) + + CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=northo, ncol_global=nextra, & + para_env=para_env, context=context) + CALL cp_fm_create(tmp_fm, tmp_fm_struct) + CALL cp_fm_struct_release(tmp_fm_struct) + CALL cp_fm_set_submatrix(tmp_fm, tmp_mat_1) + DEALLOCATE (tmp_mat_1) + CALL make_basis_simple(tmp_fm, nextra) + CALL cp_fm_to_cfm(msourcer=tmp_fm, mtarget=c_tilde) + CALL cp_fm_release(tmp_fm) + ! initialize U + IF (do_U_guess_mo) THEN + ALLOCATE (tmp_cmat(nocc, nstate)) + CALL cp_cfm_get_submatrix(V, tmp_cmat, 1, 1, nocc, nstate) + CALL cp_cfm_set_submatrix(U, tmp_cmat, 1, 1, nocc, nstate) + DEALLOCATE (tmp_cmat) + ALLOCATE (tmp_cmat(northo, nstate)) + CALL cp_cfm_get_submatrix(V, tmp_cmat, nocc + 1, 1, northo, nstate) + CALL cp_cfm_set_submatrix(VL, tmp_cmat, 1, 1, northo, nstate) + DEALLOCATE (tmp_cmat) + CALL cp_cfm_gemm("C", "N", nextra, nstate, northo, cone, c_tilde, VL, czero, UL) + ALLOCATE (tmp_cmat(nextra, nstate)) + CALL cp_cfm_get_submatrix(UL, tmp_cmat, 1, 1, nextra, nstate) + CALL cp_cfm_set_submatrix(U, tmp_cmat, nocc + 1, 1, nextra, nstate) + DEALLOCATE (tmp_cmat) + CALL cp_fm_create(tmp_fm, U%matrix_struct) + tmp_fm%local_data = REAL(U%local_data, KIND=dp) + CALL make_basis_simple(tmp_fm, nstate) + CALL cp_fm_to_cfm(msourcer=tmp_fm, mtarget=U) + CALL cp_fm_release(tmp_fm) + CALL cp_cfm_to_fm(U, matrix_U) + END IF + ! reevaluate V + ALLOCATE (tmp_cmat(nocc, nstate)) + CALL cp_cfm_get_submatrix(U, tmp_cmat, 1, 1, nocc, nstate) + CALL cp_cfm_set_submatrix(V, tmp_cmat, 1, 1, nocc, nstate) + DEALLOCATE (tmp_cmat) + ALLOCATE (tmp_cmat(nextra, nstate)) + CALL cp_cfm_get_submatrix(U, tmp_cmat, nocc + 1, 1, nextra, nstate) + CALL cp_cfm_set_submatrix(UL, tmp_cmat, 1, 1, nextra, nstate) + DEALLOCATE (tmp_cmat) + CALL cp_cfm_gemm("N", "N", northo, nstate, nextra, cone, c_tilde, UL, czero, VL) + ALLOCATE (tmp_cmat(northo, nstate)) + CALL cp_cfm_get_submatrix(VL, tmp_cmat) + CALL cp_cfm_set_submatrix(V, tmp_cmat, nocc + 1, 1, northo, nstate) + DEALLOCATE (tmp_cmat) + END SELECT + ELSE + DO idim = 1, dim2 + NULLIFY (zij_0(idim)%matrix) + CALL cp_cfm_create(zij_0(idim)%matrix, zij(1, 1)%matrix%matrix_struct) + CALL cp_cfm_to_cfm(c_zij(idim)%matrix, zij_0(idim)%matrix) + END DO + CALL cp_fm_create(rmat, zij(1, 1)%matrix%matrix_struct) + CALL cp_fm_set_all(rmat, 0._dp, 1._dp) + END IF + + unit_nr = -1 + IF (rmat%matrix_struct%para_env%mepos .EQ. rmat%matrix_struct%para_env%source) THEN + unit_nr = cp_logger_get_default_unit_nr() + WRITE (unit_nr, '(T4,A )') " Localization by combined Jacobi rotations and Non-Linear Conjugate Gradient" + END IF + + norm2_old = 1.0E30_dp + ds_min = 1.0_dp + new_direction = .TRUE. + iter = 0 + line_searches = 0 + line_search_count = 0 + tol = 1.0E+20_dp + mintol = 1.0E+10_dp + miniter = 0 + + !IF (nextra > 0) WRITE(*,*) 'random_guess, MO_guess, U_guess, conjugate_gradient: ', & + ! do_cinit_random, do_cinit_mo, do_U_guess_mo, do_cg + + ! do conjugate gradient until converged + DO WHILE (iter < max_iter) + iter = iter + 1 + !WRITE(*,*) 'iter = ', iter + t1 = m_walltime() + + IF (iter > 1) THEN + ! comput U + CALL cp_cfm_create(tmp_cfm, zij(1, 1)%matrix%matrix_struct) + CALL cp_cfm_create(tmp_cfm_2, zij(1, 1)%matrix%matrix_struct) + IF (para_env%num_pe == 1) THEN + CALL jacobi_rotations_serial_1(weights, c_zij, 1, tmp_cfm_2, tol_out=tol) + ELSE + CALL jacobi_rot_para_1(weights, c_zij, para_env, 1, tmp_cfm_2, tol_out=tol) + END IF + CALL cp_cfm_gemm('N', 'N', nstate, nstate, nstate, cone, U, tmp_cfm_2, czero, tmp_cfm) + CALL cp_cfm_to_cfm(tmp_cfm, U) + CALL cp_cfm_release(tmp_cfm) + CALL cp_cfm_release(tmp_cfm_2) + END IF + + IF (nextra > 0) THEN + ALLOCATE (tmp_cmat(nextra, nstate)) + CALL cp_cfm_get_submatrix(U, tmp_cmat, nocc + 1, 1, nextra, nstate) + CALL cp_cfm_set_submatrix(UL, tmp_cmat) + DEALLOCATE (tmp_cmat) + IF (iter > 1) THEN + ! orthonormalize c_tilde + CALL cp_fm_create(tmp_fm, c_tilde%matrix_struct) + tmp_fm%local_data = REAL(c_tilde%local_data, KIND=dp) + CALL make_basis_simple(tmp_fm, nextra) + CALL cp_fm_to_cfm(msourcer=tmp_fm, mtarget=c_tilde) + CALL cp_fm_release(tmp_fm) + + ALLOCATE (tmp_cmat(nocc, nstate)) + CALL cp_cfm_get_submatrix(U, tmp_cmat, 1, 1, nocc, nstate) + CALL cp_cfm_set_submatrix(V, tmp_cmat, 1, 1, nocc, nstate) + DEALLOCATE (tmp_cmat) + CALL cp_cfm_gemm("N", "N", northo, nstate, nextra, cone, c_tilde, UL, czero, VL) + ALLOCATE (tmp_cmat(northo, nstate)) + CALL cp_cfm_get_submatrix(VL, tmp_cmat) + CALL cp_cfm_set_submatrix(V, tmp_cmat, nocc + 1, 1, northo, nstate) + DEALLOCATE (tmp_cmat) + END IF + + ! reset if new_direction + IF (new_direction .AND. MOD(line_searches, 20) .EQ. 5) THEN + CALL cp_cfm_set_all(skc, czero) + CALL cp_cfm_set_all(Gct_old, czero) + norm2_old = 1.0E30_dp + END IF + + CALL cp_cfm_create(tmp_cfm, V%matrix_struct) + CALL cp_cfm_to_cfm(V, tmp_cfm) + CALL cp_cfm_create(tmp_cfm_1, V%matrix_struct) + ndummy = nmo + ELSE + CALL cp_cfm_create(tmp_cfm, zij(1, 1)%matrix%matrix_struct) + CALL cp_cfm_to_cfm(U, tmp_cfm) + CALL cp_cfm_create(tmp_cfm_1, zij(1, 1)%matrix%matrix_struct) + ndummy = nstate + END IF + ! update z_ij + DO idim = 1, dim2 + ! 'tmp_cfm_1 = zij_0*tmp_cfm' + CALL cp_cfm_gemm("N", "N", ndummy, nstate, ndummy, cone, zij_0(idim)%matrix, & + tmp_cfm, czero, tmp_cfm_1) + ! 'c_zij = tmp_cfm_dagg*tmp_cfm_1' + CALL cp_cfm_gemm("C", "N", nstate, nstate, ndummy, cone, tmp_cfm, tmp_cfm_1, & + czero, c_zij(idim)%matrix) + END DO + CALL cp_cfm_release(tmp_cfm) + CALL cp_cfm_release(tmp_cfm_1) + ! compute spread + DO istate = 1, nstate + spread_ii = 0.0_dp + DO idim = 1, dim2 + CALL cp_cfm_get_element(c_zij(idim)%matrix, istate, istate, mzii) + spread_ii = spread_ii + weights(idim)* & + ABS(mzii)**2/twopi/twopi + matrix_zii(istate, idim) = mzii + END DO + !WRITE(*,*) 'spread_ii', spread_ii + sum_spread(istate) = spread_ii + END DO + CALL mp_sum(spread_ii, c_zij(1)%matrix%matrix_struct%para_env%group) + spread_sum = accurate_sum(sum_spread) + + IF (nextra > 0) THEN + ! update c_tilde + CALL cp_cfm_set_all(zdiag, czero) + CALL cp_cfm_set_all(grad_ctilde, czero) + CALL cp_cfm_create(tmp_cfm, V%matrix_struct) + CALL cp_cfm_set_all(tmp_cfm, czero) + CALL cp_cfm_create(tmp_cfm_1, V%matrix_struct) + CALL cp_cfm_set_all(tmp_cfm_1, czero) + ALLOCATE (tmp_cmat(northo, nstate)) + DO idim = 1, dim2 + weight = weights(idim) + arr_zii = matrix_zii(:, idim) + ! tmp_cfm = zij_0*V + CALL cp_cfm_gemm("N", "N", nmo, nstate, nmo, cone, & + zij_0(idim)%matrix, V, czero, tmp_cfm) + ! tmp_cfm = tmp_cfm*diag_zij_dagg + CALL cp_cfm_column_scale(tmp_cfm, CONJG(arr_zii)) + ! tmp_cfm_1 = tmp_cfm*U_dagg + CALL cp_cfm_gemm("N", "C", nmo, nstate, nstate, cone, tmp_cfm, & + U, czero, tmp_cfm_1) + CALL cp_cfm_scale(weight, tmp_cfm_1) + ! zdiag = zdiag + tmp_cfm_1' + CALL cp_cfm_scale_and_add(cone, zdiag, cone, tmp_cfm_1) + + ! tmp_cfm = zij_0_dagg*V + CALL cp_cfm_gemm("C", "N", nmo, nstate, nmo, cone, & + zij_0(idim)%matrix, V, czero, tmp_cfm) + + ! tmp_cfm = tmp_cfm*diag_zij + CALL cp_cfm_column_scale(tmp_cfm, arr_zii) + ! tmp_cfm_1 = tmp_cfm*U_dagg + CALL cp_cfm_gemm("N", "C", nmo, nstate, nstate, cone, tmp_cfm, & + U, czero, tmp_cfm_1) + CALL cp_cfm_scale(weight, tmp_cfm_1) + ! zdiag = zdiag + tmp_cfm_1' + CALL cp_cfm_scale_and_add(cone, zdiag, cone, tmp_cfm_1) + END DO ! idim + CALL cp_cfm_release(tmp_cfm) + CALL cp_cfm_release(tmp_cfm_1) + DEALLOCATE (tmp_cmat) + ALLOCATE (tmp_cmat(northo, nextra)) + CALL cp_cfm_get_submatrix(zdiag, tmp_cmat, nocc + 1, nocc + 1, & + northo, nextra, .FALSE.) + ! 'grad_ctilde' + CALL cp_cfm_set_submatrix(grad_ctilde, tmp_cmat) + DEALLOCATE (tmp_cmat) + ! ctrans_lambda = c_tilde_dagg*grad_ctilde + CALL cp_cfm_gemm("C", "N", nextra, nextra, northo, cone, c_tilde, grad_ctilde, czero, ctrans_lambda) + !WRITE(*,*) "norm(ctrans_lambda) = ", cp_cfm_norm(ctrans_lambda, "F") + ! 'grad_ctilde = - c_tilde*ctrans_lambda + grad_ctilde' + CALL cp_cfm_gemm("N", "N", northo, nextra, nextra, -cone, c_tilde, ctrans_lambda, cone, grad_ctilde) + END IF ! nextra > 0 + + ! tolerance + IF (nextra > 0) THEN + tolc = 0.0_dp + CALL cp_fm_create(tmp_fm, grad_ctilde%matrix_struct) + CALL cp_cfm_to_fm(grad_ctilde, tmp_fm) + CALL cp_fm_maxabsval(tmp_fm, tolc) + CALL cp_fm_release(tmp_fm) + !WRITE(*,*) 'tolc = ', tolc + tol = tol + tolc + END IF + !WRITE(*,*) 'tol = ', tol + + IF (nextra > 0) THEN + !WRITE(*,*) 'new_direction: ', new_direction + IF (new_direction) THEN + line_searches = line_searches + 1 + IF (mintol > tol) THEN + mintol = tol + miniter = iter + END IF + + IF (unit_nr > 0 .AND. MODULO(iter, out_each) == 0) THEN + sum_spread_ii = alpha*nstate/twopi/twopi - spread_sum + avg_spread_ii = sum_spread_ii/nstate + WRITE (unit_nr, '(T4,A,T26,A,T48,A)') & + "Iteration", "Avg. Spread_ii", "Tolerance" + WRITE (unit_nr, '(T4,I7,T20,F20.10,T45,E12.4)') & + iter, avg_spread_ii, tol + CALL m_flush(unit_nr) + END IF + IF (tol < eps_localization) EXIT + + IF (do_cg) THEN + cnorm2_Gct = czero + cnorm2_Gct_cross = czero + CALL cp_cfm_trace(grad_ctilde, Gct_old, cnorm2_Gct_cross) + norm2_Gct_cross = REAL(cnorm2_Gct_cross, KIND=dp) + Gct_old%local_data = grad_ctilde%local_data + CALL cp_cfm_trace(grad_ctilde, Gct_old, cnorm2_Gct) + norm2_Gct = REAL(cnorm2_Gct, KIND=dp) + ! compute beta_pr + beta_pr = (norm2_Gct - norm2_Gct_cross)/norm2_old + norm2_old = norm2_Gct + beta = MAX(0.0_dp, beta_pr) + !WRITE(*,*) 'beta = ', beta + ! compute skc / ska = beta * skc / ska + grad_ctilde / G + CALL cp_cfm_scale(beta, skc) + CALL cp_cfm_scale_and_add(cone, skc, cone, Gct_old) + CALL cp_cfm_trace(skc, Gct_old, cnorm2_Gct_cross) + norm2_Gct_cross = REAL(cnorm2_Gct_cross, KIND=dp) + IF (norm2_Gct_cross .LE. 0.0_dp) THEN ! back to steepest ascent + CALL cp_cfm_scale_and_add(czero, skc, cone, Gct_old) + END IF + ELSE + CALL cp_cfm_scale_and_add(czero, skc, cone, grad_ctilde) + END IF + line_search_count = 0 + END IF + + line_search_count = line_search_count + 1 + !WRITE(*,*) 'line_search_count = ', line_search_count + energy(line_search_count) = spread_sum + + ! gold line search + new_direction = .FALSE. + IF (line_search_count .EQ. 1) THEN + lsl = 1 + lsr = 0 + lsm = 1 + pos(1) = 0.0_dp + pos(2) = ds_min/gold_sec + ds = pos(2) + ELSE + IF (line_search_count .EQ. 50) THEN + IF (ABS(energy(line_search_count) - energy(line_search_count - 1)) < 1.0E-4_dp) THEN + CPWARN("Line search failed to converge properly") + ds_min = 0.1_dp + new_direction = .TRUE. + ds = pos(line_search_count) + line_search_count = 0 + ELSE + CPABORT("No. of line searches exceeds 50") + END IF + ELSE + IF (lsr .EQ. 0) THEN + IF (energy(line_search_count - 1) .GT. energy(line_search_count)) THEN + lsr = line_search_count + pos(line_search_count + 1) = pos(lsm) + (pos(lsr) - pos(lsm))*gold_sec + ELSE + lsl = lsm + lsm = line_search_count + pos(line_search_count + 1) = pos(line_search_count)/gold_sec + END IF + ELSE + IF (pos(line_search_count) .LT. pos(lsm)) THEN + IF (energy(line_search_count) .GT. energy(lsm)) THEN + lsr = lsm + lsm = line_search_count + ELSE + lsl = line_search_count + END IF + ELSE + IF (energy(line_search_count) .GT. energy(lsm)) THEN + lsl = lsm + lsm = line_search_count + ELSE + lsr = line_search_count + END IF + END IF + IF (pos(lsr) - pos(lsm) .GT. pos(lsm) - pos(lsl)) THEN + pos(line_search_count + 1) = pos(lsm) + gold_sec*(pos(lsr) - pos(lsm)) + ELSE + pos(line_search_count + 1) = pos(lsl) + gold_sec*(pos(lsm) - pos(lsl)) + END IF + IF ((pos(lsr) - pos(lsl)) .LT. 1.0E-3_dp*pos(lsr)) THEN + new_direction = .TRUE. + END IF + END IF ! lsr .eq. 0 + END IF ! line_search_count .eq. 50 + ! now go to the suggested point + ds = pos(line_search_count + 1) - pos(line_search_count) + !WRITE(*,*) 'lsl, lsr, lsm, ds = ', lsl, lsr, lsm, ds + IF ((ABS(ds) < 1.0E-10_dp) .AND. (lsl == 1)) THEN + new_direction = .TRUE. + ds_min = 0.5_dp/alpha + ELSEIF (ABS(ds) > 10.0_dp) THEN + new_direction = .TRUE. + ds_min = 0.5_dp/alpha + ELSE + ds_min = pos(line_search_count + 1) + END IF + END IF ! first step + ! 'c_tilde = c_tilde + d*skc' + CALL cp_cfm_scale(ds, skc) + CALL cp_cfm_scale_and_add(cone, c_tilde, cone, skc) + ELSE + IF (mintol > tol) THEN + mintol = tol + miniter = iter + END IF + IF (unit_nr > 0 .AND. MODULO(iter, out_each) == 0) THEN + sum_spread_ii = alpha*nstate/twopi/twopi - spread_sum + avg_spread_ii = sum_spread_ii/nstate + WRITE (unit_nr, '(T4,A,T26,A,T48,A)') & + "Iteration", "Avg. Spread_ii", "Tolerance" + WRITE (unit_nr, '(T4,I7,T20,F20.10,T45,E12.4)') & + iter, avg_spread_ii, tol + CALL m_flush(unit_nr) + END IF + IF (tol < eps_localization) EXIT + END IF ! nextra > 0 + + END DO ! iteration + + IF ((unit_nr > 0) .AND. (iter == max_iter)) THEN + WRITE (unit_nr, '(T4,A,T4,A)') "Min. Itr.", "Min. Tol." + WRITE (unit_nr, '(T4,I7,T4,E12.4)') miniter, mintol + CALL m_flush(unit_nr) + END IF + + CALL cp_cfm_to_fm(U, matrix_U) + + IF (nextra > 0) THEN + rmat%local_data = REAL(V%local_data, KIND=dp) + CALL rotate_orbitals_edf(rmat, vectors_all, vectors) + + CALL cp_cfm_release(c_tilde) + CALL cp_cfm_release(grad_ctilde) + CALL cp_cfm_release(Gct_old) + CALL cp_cfm_release(skc) + CALL cp_cfm_release(UL) + CALL cp_cfm_release(zdiag) + CALL cp_cfm_release(ctrans_lambda) + CALL cp_fm_release(id_nextra) + CALL cp_fm_release(vectors_all) + CALL cp_cfm_release(V) + CALL cp_fm_release(matrix_V) + CALL cp_fm_release(matrix_V_all) + CALL cp_cfm_release(VL) + DEALLOCATE (arr_zii) + ELSE + rmat%local_data = matrix_U%local_data + CALL rotate_orbitals(rmat, vectors) + END IF + DO idim = 1, dim2 + CALL cp_cfm_release(zij_0(idim)%matrix) + END DO + DEALLOCATE (zij_0) + + DO idim = 1, dim2 + zij(1, idim)%matrix%local_data = REAL(c_zij(idim)%matrix%local_data, dp) + zij(2, idim)%matrix%local_data = AIMAG(c_zij(idim)%matrix%local_data) + CALL cp_cfm_release(c_zij(idim)%matrix) + END DO + DEALLOCATE (c_zij) + CALL cp_fm_release(rmat) + CALL cp_cfm_release(U) + CALL cp_fm_release(matrix_U) + DEALLOCATE (matrix_zii, sum_spread) + + CALL timestop(handle) + + END SUBROUTINE jacobi_cg_edf_ls +! ************************************************************************************************** !> \brief ... !> \param istate ... !> \param jstate ... @@ -491,11 +1318,14 @@ CONTAINS !> \param mij ... !> \param weights ... !> \param theta ... +!> \param grad_ij ... +!> \param step ... ! ************************************************************************************************** - SUBROUTINE get_angle(mii, mjj, mij, weights, theta) + SUBROUTINE get_angle(mii, mjj, mij, weights, theta, grad_ij, step) COMPLEX(KIND=dp), POINTER :: mii(:), mjj(:), mij(:) REAL(KIND=dp), INTENT(IN) :: weights(:) REAL(KIND=dp), INTENT(OUT) :: theta + REAL(KIND=dp), INTENT(IN), OPTIONAL :: grad_ij, step COMPLEX(KIND=dp) :: z11, z12, z22 INTEGER :: dim_m, idim @@ -521,6 +1351,7 @@ CONTAINS ELSE theta = 0.25_dp*pi END IF + IF (PRESENT(grad_ij)) theta = theta + step*grad_ij ! Check second derivative info d2 = a12*SIN(4._dp*theta) - b12*COS(4._dp*theta) IF (d2 <= 0._dp) THEN ! go to the maximum, not the minimum @@ -536,11 +1367,13 @@ CONTAINS !> \param zij ... !> \param weights ... !> \param tolerance ... +!> \param grad ... ! ************************************************************************************************** - SUBROUTINE check_tolerance(zij, weights, tolerance) + SUBROUTINE check_tolerance(zij, weights, tolerance, grad) TYPE(cp_cfm_p_type) :: zij(:) REAL(KIND=dp), INTENT(IN) :: weights(:) REAL(KIND=dp), INTENT(OUT) :: tolerance + TYPE(cp_fm_type), INTENT(OUT), OPTIONAL, POINTER :: grad CHARACTER(len=*), PARAMETER :: routineN = 'check_tolerance' @@ -556,6 +1389,7 @@ CONTAINS CALL cp_fm_set_all(force, 0._dp) CALL grad_at_0(zij, weights, force) CALL cp_fm_maxabsval(force, tolerance) + IF (PRESENT(grad)) CALL cp_fm_to_fm(force, grad) CALL cp_fm_release(force) CALL timestop(handle) @@ -581,6 +1415,27 @@ CONTAINS END SUBROUTINE rotate_orbitals ! ************************************************************************************************** !> \brief ... +!> \param rmat ... +!> \param vec_all ... +!> \param vectors ... +! ************************************************************************************************** + SUBROUTINE rotate_orbitals_edf(rmat, vec_all, vectors) + TYPE(cp_fm_type), POINTER :: rmat, vec_all, vectors + + INTEGER :: k, l, n + TYPE(cp_fm_type), POINTER :: wf + + NULLIFY (wf) + CALL cp_fm_create(wf, vectors%matrix_struct) + CALL cp_fm_get_info(vec_all, nrow_global=n, ncol_global=k) + CALL cp_fm_get_info(rmat, ncol_global=l) + + CALL cp_gemm("N", "N", n, l, k, 1.0_dp, vec_all, rmat, 0.0_dp, wf) + CALL cp_fm_to_fm(wf, vectors) + CALL cp_fm_release(wf) + END SUBROUTINE rotate_orbitals_edf +! ************************************************************************************************** +!> \brief ... !> \param diag ... !> \param weights ... !> \param matrix ... @@ -1361,71 +2216,36 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'jacobi_rot_para' - COMPLEX(KIND=dp) :: zi, zj - COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: c_array_me, c_array_partner - COMPLEX(KIND=dp), POINTER :: mii(:), mij(:), mjj(:) - INTEGER :: dim2, handle, i, idim, ii, ik, il1, il2, il_recv, il_recv_partner, ilow1, ilow2, & - ip, ip_has_i, ip_partner, ip_recv_from, ip_recv_partner, ipair, iperm, istat, istate, & - iu1, iu2, iup1, iup2, j, jj, jstate, k, kk, n1, n2, nblock, nblock_max, npair, nperm, & - ns_me, ns_partner, ns_recv_from, ns_recv_partner, nstate, output_unit - INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount, rdispl - INTEGER, ALLOCATABLE, DIMENSION(:, :) :: list_pair, ns_bound - LOGICAL :: should_stop - REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gmat, rmat_loc, rmat_recv, rmat_send, & - rotmat, z_ij_loc_im, z_ij_loc_re - REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: rmat_recv_all - REAL(KIND=dp) :: ct, func, gmax, grad, ri, rj, st, t1, & - t2, theta, tolerance, xlow, xstate, & - xup, zc, zr + INTEGER :: dim2, handle, i, idim, ii, ilow1, ip, j, & + nblock, nblock_max, ns_me, nstate, & + output_unit + INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ns_bound + REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rotmat, z_ij_loc_im, z_ij_loc_re + REAL(KIND=dp) :: xstate TYPE(cp_fm_type), POINTER :: rmat - TYPE(set_c_1d_type), DIMENSION(:), POINTER :: zdiag_all, zdiag_me - TYPE(set_c_2d_type), DIMENSION(:), POINTER :: cz_ij_loc, xyz_mix, xyz_mix_ns + TYPE(set_c_2d_type), DIMENSION(:), POINTER :: cz_ij_loc CALL timeset(routineN, handle) output_unit = cp_logger_get_default_io_unit() - NULLIFY (rmat, cz_ij_loc, zdiag_all, zdiag_me) - NULLIFY (xyz_mix, xyz_mix_ns) - NULLIFY (mii, mij, mjj) + NULLIFY (rmat, cz_ij_loc) dim2 = SIZE(zij, 2) - ALLOCATE (mii(dim2), mij(dim2), mjj(dim2)) CALL cp_fm_create(rmat, zij(1, 1)%matrix%matrix_struct) CALL cp_fm_set_all(rmat, 0._dp, 1._dp) CALL cp_fm_get_info(rmat, nrow_global=nstate) - ALLOCATE (rcount(para_env%num_pe), STAT=istat) - ALLOCATE (rdispl(para_env%num_pe), STAT=istat) - - tolerance = 1.0e10_dp - sweeps = 0 - - ! number of processor pairs and number of permutations - npair = (para_env%num_pe + 1)/2 - nperm = para_env%num_pe - MOD(para_env%num_pe + 1, 2) - ALLOCATE (list_pair(2, npair)) - ! Distribution of the states (XXXXX safe against more pe than states ??? XXXXX) xstate = REAL(nstate, dp)/REAL(para_env%num_pe, dp) - nblock_max = 0 ALLOCATE (ns_bound(0:para_env%num_pe - 1, 2)) - Xlow = 0.0D0 - Xup = 0.0D0 DO ip = 1, para_env%num_pe - xup = xlow + xstate - ns_bound(ip - 1, 1) = NINT(xlow) + 1 - ns_bound(ip - 1, 2) = NINT(xup) - IF (NINT(xup) .GT. nstate) THEN - ns_bound(ip - 1, 2) = nstate - END IF - IF (NINT(xlow) .GT. nstate) THEN - ns_bound(ip - 1, 1) = nstate + 1 - END IF - xlow = xup + ns_bound(ip - 1, 1) = MIN(nstate, NINT(xstate*(ip - 1))) + 1 + ns_bound(ip - 1, 2) = MIN(nstate, NINT(xstate*ip)) END DO + nblock_max = 0 DO ip = 0, para_env%num_pe - 1 nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1 nblock_max = MAX(nblock_max, nblock) @@ -1441,9 +2261,8 @@ CONTAINS CALL cp_fm_get_submatrix(zij(1, idim)%matrix, z_ij_loc_re, 1, ns_bound(ip, 1), nstate, nblock) CALL cp_fm_get_submatrix(zij(2, idim)%matrix, z_ij_loc_im, 1, ns_bound(ip, 1), nstate, nblock) IF (para_env%mepos == ip) THEN - ns_me = nblock - ALLOCATE (cz_ij_loc(idim)%c_array(nstate, ns_me)) - DO i = 1, ns_me + ALLOCATE (cz_ij_loc(idim)%c_array(nstate, nblock)) + DO i = 1, nblock DO j = 1, nstate cz_ij_loc(idim)%c_array(j, i) = CMPLX(z_ij_loc_re(j, i), z_ij_loc_im(j, i), dp) END DO @@ -1454,8 +2273,258 @@ CONTAINS DEALLOCATE (z_ij_loc_re) DEALLOCATE (z_ij_loc_im) - ! initialize rotation matrix ALLOCATE (rotmat(nstate, 2*nblock_max)) + + CALL jacobi_rot_para_core(weights, para_env, max_iter, sweeps, out_each, dim2, nstate, nblock_max, ns_bound, & + cz_ij_loc, rotmat, output_unit, eps_localization=eps_localization, & + target_time=target_time, start_time=start_time) + + ilow1 = ns_bound(para_env%mepos, 1) + ns_me = ns_bound(para_env%mepos, 2) - ns_bound(para_env%mepos, 1) + 1 + ALLOCATE (z_ij_loc_re(nstate, nblock_max)) + ALLOCATE (z_ij_loc_im(nstate, nblock_max)) + DO idim = 1, dim2 + DO ip = 0, para_env%num_pe - 1 + z_ij_loc_re = 0.0_dp + z_ij_loc_im = 0.0_dp + nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1 + IF (ip == para_env%mepos) THEN + ns_me = nblock + DO i = 1, ns_me + ii = ilow1 + i - 1 + DO j = 1, nstate + z_ij_loc_re(j, i) = REAL(cz_ij_loc(idim)%c_array(j, i), dp) + z_ij_loc_im(j, i) = AIMAG(cz_ij_loc(idim)%c_array(j, i)) + END DO + END DO + END IF + CALL mp_bcast(z_ij_loc_re, ip, para_env%group) + CALL mp_bcast(z_ij_loc_im, ip, para_env%group) + CALL cp_fm_set_submatrix(zij(1, idim)%matrix, z_ij_loc_re, 1, ns_bound(ip, 1), nstate, nblock) + CALL cp_fm_set_submatrix(zij(2, idim)%matrix, z_ij_loc_im, 1, ns_bound(ip, 1), nstate, nblock) + END DO ! ip + END DO + + DO ip = 0, para_env%num_pe - 1 + z_ij_loc_re = 0.0_dp + nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1 + IF (ip == para_env%mepos) THEN + ns_me = nblock + DO i = 1, ns_me + ii = ilow1 + i - 1 + DO j = 1, nstate + z_ij_loc_re(j, i) = rotmat(j, i) + END DO + END DO + END IF + CALL mp_bcast(z_ij_loc_re, ip, para_env%group) + CALL cp_fm_set_submatrix(rmat, z_ij_loc_re, 1, ns_bound(ip, 1), nstate, nblock) + END DO + + DEALLOCATE (z_ij_loc_re) + DEALLOCATE (z_ij_loc_im) + DO idim = 1, dim2 + DEALLOCATE (cz_ij_loc(idim)%c_array) + END DO + DEALLOCATE (cz_ij_loc) + + CALL mp_sync(para_env%group) + CALL rotate_orbitals(rmat, vectors) + CALL cp_fm_release(rmat) + + DEALLOCATE (rotmat) + DEALLOCATE (ns_bound) + + CALL timestop(handle) + + END SUBROUTINE jacobi_rot_para + +! ************************************************************************************************** +!> \brief almost identical to 'jacobi_rot_para' but with different inout variables +!> \param weights ... +!> \param czij ... +!> \param para_env ... +!> \param max_iter ... +!> \param rmat ... +!> \param tol_out ... +!> \author Soumya Ghosh (08/21) +! ************************************************************************************************** + SUBROUTINE jacobi_rot_para_1(weights, czij, para_env, max_iter, rmat, tol_out) + + REAL(KIND=dp), INTENT(IN) :: weights(:) + TYPE(cp_cfm_p_type), INTENT(IN), POINTER :: czij(:) + TYPE(cp_para_env_type), POINTER :: para_env + INTEGER, INTENT(IN) :: max_iter + TYPE(cp_cfm_type), INTENT(INOUT), POINTER :: rmat + REAL(dp), INTENT(OUT), OPTIONAL :: tol_out + + CHARACTER(len=*), PARAMETER :: routineN = 'jacobi_rot_para_1' + + COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: czij_array + INTEGER :: dim2, handle, i, idim, ii, ilow1, ip, j, & + nblock, nblock_max, ns_me, nstate, & + sweeps + INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ns_bound + REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rotmat, z_ij_loc_re + REAL(KIND=dp) :: xstate + TYPE(set_c_2d_type), DIMENSION(:), POINTER :: cz_ij_loc + + CALL timeset(routineN, handle) + + dim2 = SIZE(czij) + + CALL cp_cfm_create(rmat, czij(1)%matrix%matrix_struct) + CALL cp_cfm_set_all(rmat, CMPLX(0._dp, 0._dp, dp), CMPLX(1._dp, 0._dp, dp)) + + CALL cp_cfm_get_info(rmat, nrow_global=nstate) + + ! Distribution of the states (XXXXX safe against more pe than states ??? XXXXX) + xstate = REAL(nstate, dp)/REAL(para_env%num_pe, dp) + ALLOCATE (ns_bound(0:para_env%num_pe - 1, 2)) + DO ip = 1, para_env%num_pe + ns_bound(ip - 1, 1) = MIN(nstate, NINT(xstate*(ip - 1))) + 1 + ns_bound(ip - 1, 2) = MIN(nstate, NINT(xstate*ip)) + END DO + nblock_max = 0 + DO ip = 0, para_env%num_pe - 1 + nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1 + nblock_max = MAX(nblock_max, nblock) + END DO + + ! otbtain local part of the matrix (could be made faster, but is likely irrelevant). + ALLOCATE (czij_array(nstate, nblock_max)) + ALLOCATE (cz_ij_loc(dim2)) + DO idim = 1, dim2 + DO ip = 0, para_env%num_pe - 1 + nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1 + ! cfm --> allocatable + CALL cp_cfm_get_submatrix(czij(idim)%matrix, czij_array, 1, ns_bound(ip, 1), nstate, nblock) + IF (para_env%mepos == ip) THEN + ns_me = nblock + ALLOCATE (cz_ij_loc(idim)%c_array(nstate, ns_me)) + DO i = 1, ns_me + DO j = 1, nstate + cz_ij_loc(idim)%c_array(j, i) = czij_array(j, i) + END DO + END DO + END IF + END DO ! ip + END DO + DEALLOCATE (czij_array) + + ALLOCATE (rotmat(nstate, 2*nblock_max)) + + CALL jacobi_rot_para_core(weights, para_env, max_iter, sweeps, 1, dim2, nstate, nblock_max, ns_bound, & + cz_ij_loc, rotmat, 0, tol_out=tol_out) + + ilow1 = ns_bound(para_env%mepos, 1) + ns_me = ns_bound(para_env%mepos, 2) - ns_bound(para_env%mepos, 1) + 1 + ALLOCATE (z_ij_loc_re(nstate, nblock_max)) + + DO ip = 0, para_env%num_pe - 1 + z_ij_loc_re = 0.0_dp + nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1 + IF (ip == para_env%mepos) THEN + ns_me = nblock + DO i = 1, ns_me + ii = ilow1 + i - 1 + DO j = 1, nstate + z_ij_loc_re(j, i) = rotmat(j, i) + END DO + END DO + END IF + CALL mp_bcast(z_ij_loc_re, ip, para_env%group) + CALL cp_cfm_set_submatrix(rmat, CMPLX(z_ij_loc_re, 0.0_dp, dp), 1, ns_bound(ip, 1), nstate, nblock) + END DO + + DEALLOCATE (z_ij_loc_re) + DO idim = 1, dim2 + DEALLOCATE (cz_ij_loc(idim)%c_array) + END DO + DEALLOCATE (cz_ij_loc) + + CALL mp_sync(para_env%group) + + DEALLOCATE (rotmat) + DEALLOCATE (ns_bound) + + CALL timestop(handle) + + END SUBROUTINE jacobi_rot_para_1 + +! ************************************************************************************************** +!> \brief Parallel algorithm for jacobi rotations +!> \param weights ... +!> \param para_env ... +!> \param max_iter ... +!> \param sweeps ... +!> \param out_each ... +!> \param dim2 ... +!> \param nstate ... +!> \param nblock_max ... +!> \param ns_bound ... +!> \param cz_ij_loc ... +!> \param rotmat ... +!> \param output_unit ... +!> \param tol_out ... +!> \param eps_localization ... +!> \param target_time ... +!> \param start_time ... +!> \par History +!> split out to reuse with different input types +!> \author HF (05.2022) +! ************************************************************************************************** + SUBROUTINE jacobi_rot_para_core(weights, para_env, max_iter, sweeps, out_each, dim2, nstate, nblock_max, & + ns_bound, cz_ij_loc, rotmat, output_unit, tol_out, eps_localization, target_time, start_time) + + REAL(KIND=dp), INTENT(IN) :: weights(:) + TYPE(cp_para_env_type), POINTER :: para_env + INTEGER, INTENT(IN) :: max_iter + INTEGER, INTENT(OUT) :: sweeps + INTEGER, INTENT(IN) :: out_each, dim2, nstate, nblock_max + INTEGER, DIMENSION(0:, :), INTENT(IN) :: ns_bound + TYPE(set_c_2d_type), DIMENSION(:), POINTER :: cz_ij_loc + REAL(dp), DIMENSION(:, :), INTENT(OUT) :: rotmat + INTEGER, INTENT(IN) :: output_unit + REAL(dp), INTENT(OUT), OPTIONAL :: tol_out + REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_localization + REAL(dp), OPTIONAL :: target_time, start_time + + COMPLEX(KIND=dp) :: zi, zj + COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: c_array_me, c_array_partner + COMPLEX(KIND=dp), POINTER :: mii(:), mij(:), mjj(:) + INTEGER :: i, idim, ii, ik, il1, il2, il_recv, il_recv_partner, ilow1, ilow2, ip, ip_has_i, & + ip_partner, ip_recv_from, ip_recv_partner, ipair, iperm, istat, istate, iu1, iu2, iup1, & + iup2, j, jj, jstate, k, kk, lsweep, n1, n2, npair, nperm, ns_me, ns_partner, & + ns_recv_from, ns_recv_partner + INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount, rdispl + INTEGER, ALLOCATABLE, DIMENSION(:, :) :: list_pair + LOGICAL :: should_stop + REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gmat, rmat_loc, rmat_recv, rmat_send + REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: rmat_recv_all + REAL(KIND=dp) :: ct, func, gmax, grad, ri, rj, st, t1, & + t2, theta, tolerance, zc, zr + TYPE(set_c_1d_type), DIMENSION(:), POINTER :: zdiag_all, zdiag_me + TYPE(set_c_2d_type), DIMENSION(:), POINTER :: xyz_mix, xyz_mix_ns + + NULLIFY (zdiag_all, zdiag_me) + NULLIFY (xyz_mix, xyz_mix_ns) + NULLIFY (mii, mij, mjj) + + ALLOCATE (mii(dim2), mij(dim2), mjj(dim2)) + + ALLOCATE (rcount(para_env%num_pe), STAT=istat) + ALLOCATE (rdispl(para_env%num_pe), STAT=istat) + + tolerance = 1.0e10_dp + sweeps = 0 + + ! number of processor pairs and number of permutations + npair = (para_env%num_pe + 1)/2 + nperm = para_env%num_pe - MOD(para_env%num_pe + 1, 2) + ALLOCATE (list_pair(2, npair)) + + ! initialize rotation matrix rotmat = 0.0_dp DO i = ns_bound(para_env%mepos, 1), ns_bound(para_env%mepos, 2) ii = i - ns_bound(para_env%mepos, 1) + 1 @@ -1493,7 +2562,8 @@ CONTAINS WRITE (output_unit, '(T20,A12,T32, A22,T60, A12,A8 )') "Iteration", "Functional", "Tolerance", " Time " END IF - DO sweeps = 1, max_iter + 1 + DO lsweep = 1, max_iter + 1 + sweeps = lsweep IF (sweeps == max_iter + 1) THEN IF (output_unit > 0) THEN WRITE (output_unit, *) ' LOCALIZATION! loop did not converge within the maximum number of iterations.' @@ -1615,12 +2685,21 @@ CONTAINS rotmat(1:nstate, k:k + n2 - 1), ip_partner, para_env%group) IF (ilow1 < ilow2) THEN - CALL dgemm("N", "N", nstate, n1, n2, 1.0_dp, rotmat(1, k), nstate, rmat_loc(1 + n1, 1), n1 + n2, 0.0_dp, gmat, nstate) - CALL dgemm("N", "N", nstate, n1, n1, 1.0_dp, rotmat(1, 1), nstate, rmat_loc(1, 1), n1 + n2, 1.0_dp, gmat, nstate) + ! no longer compiles in official sdgb: + !CALL dgemm("N", "N", nstate, n1, n2, 1.0_dp, rotmat(1, k), nstate, rmat_loc(1 + n1, 1), n1 + n2, 0.0_dp, gmat, nstate) + ! probably inefficient: + CALL dgemm("N", "N", nstate, n1, n2, 1.0_dp, rotmat(1:, k:), nstate, rmat_loc(1 + n1:, 1:n1), & + n2, 0.0_dp, gmat(:, :), nstate) + CALL dgemm("N", "N", nstate, n1, n1, 1.0_dp, rotmat(1:, 1:), nstate, rmat_loc(1:, 1:), & + n1 + n2, 1.0_dp, gmat(:, :), nstate) ELSE - CALL dgemm("N", "N", nstate, n1, n2, 1.0_dp, rotmat(1, k), nstate, rmat_loc(1, n2 + 1), n1 + n2, 0.0_dp, gmat, nstate) - CALL dgemm("N", "N", nstate, n1, n1, 1.0_dp, rotmat(1, 1), nstate, & - rmat_loc(n2 + 1, n2 + 1), n1 + n2, 1.0_dp, gmat, nstate) + CALL dgemm("N", "N", nstate, n1, n2, 1.0_dp, rotmat(1:, k:), nstate, & + rmat_loc(1:, n2 + 1:), n1 + n2, 0.0_dp, gmat(:, :), nstate) + ! no longer compiles in official sdgb: + !CALL dgemm("N", "N", nstate, n1, n1, 1.0_dp, rotmat(1, 1), nstate, rmat_loc(n2 + 1, n2 + 1), n1 + n2, 1.0_dp, gmat, nstate) + ! probably inefficient: + CALL dgemm("N", "N", nstate, n1, n1, 1.0_dp, rotmat(1:, 1:), nstate, rmat_loc(n2 + 1:, n2 + 1:), & + n1, 1.0_dp, gmat(:, :), nstate) END IF CALL dcopy(nstate*n1, gmat(1, 1), 1, rotmat(1, 1), 1) @@ -1814,6 +2893,7 @@ CONTAINS CALL mp_max(gmax, para_env%group) tolerance = gmax + IF (PRESENT(tol_out)) tol_out = tolerance func = 0.0_dp DO i = ns_bound(para_env%mepos, 1), ns_bound(para_env%mepos, 2) @@ -1831,11 +2911,15 @@ CONTAINS WRITE (output_unit, '(T20,I12,T35,F20.10,T60,E12.4,F8.3)') sweeps, func, tolerance, t2 - t1 CALL m_flush(output_unit) END IF - IF (tolerance < eps_localization) EXIT - CALL external_control(should_stop, "LOC", target_time=target_time, start_time=start_time) - IF (should_stop) EXIT + IF (PRESENT(eps_localization)) THEN + IF (tolerance < eps_localization) EXIT + END IF + IF (PRESENT(target_time) .AND. PRESENT(start_time)) THEN + CALL external_control(should_stop, "LOC", target_time=target_time, start_time=start_time) + IF (should_stop) EXIT + END IF - END DO ! sweeps + END DO ! lsweep ! buffer for message passing DEALLOCATE (rmat_recv_all) @@ -1864,66 +2948,9 @@ CONTAINS DEALLOCATE (mii) DEALLOCATE (mij) DEALLOCATE (mjj) + DEALLOCATE (list_pair) - ilow1 = ns_bound(para_env%mepos, 1) - ns_me = ns_bound(para_env%mepos, 2) - ns_bound(para_env%mepos, 1) + 1 - ALLOCATE (z_ij_loc_re(nstate, nblock_max)) - ALLOCATE (z_ij_loc_im(nstate, nblock_max)) - DO idim = 1, dim2 - DO ip = 0, para_env%num_pe - 1 - z_ij_loc_re = 0.0_dp - z_ij_loc_im = 0.0_dp - nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1 - IF (ip == para_env%mepos) THEN - ns_me = nblock - DO i = 1, ns_me - ii = ilow1 + i - 1 - DO j = 1, nstate - z_ij_loc_re(j, i) = REAL(cz_ij_loc(idim)%c_array(j, i), dp) - z_ij_loc_im(j, i) = AIMAG(cz_ij_loc(idim)%c_array(j, i)) - END DO - END DO - END IF - CALL mp_bcast(z_ij_loc_re, ip, para_env%group) - CALL mp_bcast(z_ij_loc_im, ip, para_env%group) - CALL cp_fm_set_submatrix(zij(1, idim)%matrix, z_ij_loc_re, 1, ns_bound(ip, 1), nstate, nblock) - CALL cp_fm_set_submatrix(zij(2, idim)%matrix, z_ij_loc_im, 1, ns_bound(ip, 1), nstate, nblock) - END DO ! ip - END DO - - DO ip = 0, para_env%num_pe - 1 - z_ij_loc_re = 0.0_dp - nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1 - IF (ip == para_env%mepos) THEN - ns_me = nblock - DO i = 1, ns_me - ii = ilow1 + i - 1 - DO j = 1, nstate - z_ij_loc_re(j, i) = rotmat(j, i) - END DO - END DO - END IF - CALL mp_bcast(z_ij_loc_re, ip, para_env%group) - CALL cp_fm_set_submatrix(rmat, z_ij_loc_re, 1, ns_bound(ip, 1), nstate, nblock) - END DO - - DEALLOCATE (z_ij_loc_re) - DEALLOCATE (z_ij_loc_im) - DO idim = 1, dim2 - DEALLOCATE (cz_ij_loc(idim)%c_array) - END DO - DEALLOCATE (cz_ij_loc) - - CALL mp_sync(para_env%group) - CALL rotate_orbitals(rmat, vectors) - CALL cp_fm_release(rmat) - - DEALLOCATE (rotmat) - DEALLOCATE (ns_bound, list_pair) - - CALL timestop(handle) - - END SUBROUTINE jacobi_rot_para + END SUBROUTINE jacobi_rot_para_core ! ************************************************************************************************** !> \brief ... diff --git a/src/qs_scf_post_gpw.F b/src/qs_scf_post_gpw.F index 9587d1674e..c0d8ec630f 100644 --- a/src/qs_scf_post_gpw.F +++ b/src/qs_scf_post_gpw.F @@ -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 diff --git a/tests/QS/regtest-loc_powf/TEST_FILES b/tests/QS/regtest-loc_powf/TEST_FILES new file mode 100644 index 0000000000..5a2475a6b8 --- /dev/null +++ b/tests/QS/regtest-loc_powf/TEST_FILES @@ -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 diff --git a/tests/QS/regtest-loc_powf/pos_4Au_H2O_4Ne.xyz b/tests/QS/regtest-loc_powf/pos_4Au_H2O_4Ne.xyz new file mode 100644 index 0000000000..42e9ae46f2 --- /dev/null +++ b/tests/QS/regtest-loc_powf/pos_4Au_H2O_4Ne.xyz @@ -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 diff --git a/tests/QS/regtest-loc_powf/run.inp b/tests/QS/regtest-loc_powf/run.inp new file mode 100644 index 0000000000..313b1e9d05 --- /dev/null +++ b/tests/QS/regtest-loc_powf/run.inp @@ -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 diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 7419ade25f..c2c7ce3529 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -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 diff --git a/tests/TEST_TYPES b/tests/TEST_TYPES index 8af43134a3..1b51d3b17a 100644 --- a/tests/TEST_TYPES +++ b/tests/TEST_TYPES @@ -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