From 2de90c3493f93ab5519e4c813f6f1ead4f3a97fc Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Tue, 3 Dec 2024 11:43:06 +0100 Subject: [PATCH] EHT initial guess (#3790) --- src/CMakeLists.txt | 1 + src/dm_ls_scf.F | 3 + src/dm_ls_scf_create.F | 9 + src/input_constants.F | 3 +- src/input_cp2k_scf.F | 13 +- src/mixed_cdft_utils.F | 4 +- src/qs_eht_guess.F | 238 ++++++++++++++++++++++++++ src/qs_environment.F | 22 ++- src/qs_gamma2kp.F | 2 +- src/qs_initial_guess.F | 27 +-- src/qs_kind_types.F | 29 ++-- src/qs_subsys_methods.F | 12 +- tests/QS/regtest-eht-guess/TEST_FILES | 9 + tests/QS/regtest-eht-guess/h2o_1.inp | 46 +++++ tests/QS/regtest-eht-guess/h2o_2.inp | 46 +++++ tests/TEST_DIRS | 1 + 16 files changed, 418 insertions(+), 47 deletions(-) create mode 100644 src/qs_eht_guess.F create mode 100644 tests/QS/regtest-eht-guess/TEST_FILES create mode 100644 tests/QS/regtest-eht-guess/h2o_1.inp create mode 100644 tests/QS/regtest-eht-guess/h2o_2.inp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index dd6c95777a..f7a24a2a44 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -591,6 +591,7 @@ list( qs_harmonics_atom.F qs_hash_table_functions.F qs_initial_guess.F + qs_eht_guess.F qs_integral_utils.F qs_integrate_potential.F qs_integrate_potential_product.F diff --git a/src/dm_ls_scf.F b/src/dm_ls_scf.F index 948b7f5195..354509521a 100644 --- a/src/dm_ls_scf.F +++ b/src/dm_ls_scf.F @@ -28,6 +28,7 @@ MODULE dm_ls_scf cp_logger_get_default_unit_nr,& cp_logger_type USE dm_ls_chebyshev, ONLY: compute_chebyshev + USE dm_ls_scf_create, ONLY: ls_scf_create USE dm_ls_scf_curvy, ONLY: deallocate_curvy_data,& dm_ls_curvy_optimization USE dm_ls_scf_methods, ONLY: apply_matrix_preconditioner,& @@ -116,6 +117,8 @@ SUBROUTINE ls_scf(qs_env, nonscf) do_scf = .TRUE. IF (PRESENT(nonscf)) do_scf = .NOT. nonscf + ! Moved here from qs_environment to remove dependencies + CALL ls_scf_create(qs_env) CALL get_qs_env(qs_env, ls_scf_env=ls_scf_env) IF (do_scf) THEN diff --git a/src/dm_ls_scf_create.F b/src/dm_ls_scf_create.F index fcad2d2261..91d7d11888 100644 --- a/src/dm_ls_scf_create.F +++ b/src/dm_ls_scf_create.F @@ -92,6 +92,15 @@ SUBROUTINE ls_scf_create(qs_env) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(section_vals_type), POINTER :: input, mixing_section + NULLIFY (ls_scf_env) + CALL get_qs_env(qs_env, ls_scf_env=ls_scf_env) + ! we cannot rebuild ls_scf_env for each new optimization. It seems it holds some values + ! that need to be save between calls? + IF (ASSOCIATED(ls_scf_env)) RETURN +! IF(ASSOCIATED(ls_scf_env)) THEN +! CALL ls_scf_release(ls_scf_env) +! END IF + CALL timeset(routineN, handle) CALL cite_reference(VandeVondele2012) diff --git a/src/input_constants.F b/src/input_constants.F index 7c50b13c71..9198b3cd72 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -208,7 +208,8 @@ MODULE input_constants mopac_guess = 7, & sparse_guess = 8, & almo_scf_guess = 9, & - molecular_guess = 10 + molecular_guess = 10, & + eht_guess = 11 ! Projection formulas for the maximum overlap method INTEGER, PARAMETER, PUBLIC :: momproj_norm = 0, & diff --git a/src/input_cp2k_scf.F b/src/input_cp2k_scf.F index c06d7abada..32ee4e0942 100644 --- a/src/input_cp2k_scf.F +++ b/src/input_cp2k_scf.F @@ -33,10 +33,10 @@ MODULE input_cp2k_scf cdft_alpha_constraint, cdft_beta_constraint, cdft_charge_constraint, & cdft_magnetization_constraint, cholesky_dbcsr, cholesky_inverse, cholesky_off, & cholesky_reduce, cholesky_restore, core_guess, diag_block_davidson, diag_block_krylov, & - diag_filter_matrix, diag_ot, diag_standard, gaussian, general_roks, high_spin_roks, & - history_guess, jacobian_fd1, jacobian_fd1_backward, jacobian_fd1_central, jacobian_fd2, & - jacobian_fd2_backward, ls_2pnt, ls_3pnt, ls_gold, ls_none, mopac_guess, no_guess, & - numerical, ot_algo_irac, ot_algo_taylor_or_diag, ot_chol_irac, ot_lwdn_irac, & + diag_filter_matrix, diag_ot, diag_standard, eht_guess, gaussian, general_roks, & + high_spin_roks, history_guess, jacobian_fd1, jacobian_fd1_backward, jacobian_fd1_central, & + jacobian_fd2, jacobian_fd2_backward, ls_2pnt, ls_3pnt, ls_gold, ls_none, mopac_guess, & + no_guess, numerical, ot_algo_irac, ot_algo_taylor_or_diag, ot_chol_irac, ot_lwdn_irac, & ot_mini_broyden, ot_mini_cg, ot_mini_diis, ot_mini_sd, ot_poly_irac, ot_precond_full_all, & ot_precond_full_kinetic, ot_precond_full_single, ot_precond_full_single_inverse, & ot_precond_none, ot_precond_s_inverse, ot_precond_solver_default, & @@ -207,17 +207,18 @@ SUBROUTINE create_scf_section(section) description="Change the initial guess for the wavefunction.", & usage="SCF_GUESS RESTART", default_i_val=atomic_guess, & enum_c_vals=s2a("ATOMIC", "RESTART", "RANDOM", "CORE", & - "HISTORY_RESTART", "MOPAC", "SPARSE", "NONE"), & + "HISTORY_RESTART", "MOPAC", "EHT", "SPARSE", "NONE"), & enum_desc=s2a("Generate an atomic density using the atomic code", & "Use the RESTART file as an initial guess (and ATOMIC if not present).", & "Use random wavefunction coefficients.", & "Diagonalize the core hamiltonian for an initial guess.", & "Extrapolated from previous RESTART files.", & "Use same guess as MOPAC for semi-empirical methods or a simple diagonal density matrix for other methods", & + "Use the EHT (gfn0-xTB) code to generate an initial wavefunction.", & "Generate a sparse wavefunction using the atomic code (for OT based methods)", & "Skip initial guess (only for non-self consistent methods)."), & enum_i_vals=(/atomic_guess, restart_guess, random_guess, core_guess, & - history_guess, mopac_guess, sparse_guess, no_guess/)) + history_guess, mopac_guess, eht_guess, sparse_guess, no_guess/)) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) diff --git a/src/mixed_cdft_utils.F b/src/mixed_cdft_utils.F index 9aedb3368f..f6375ea5e4 100644 --- a/src/mixed_cdft_utils.F +++ b/src/mixed_cdft_utils.F @@ -803,7 +803,7 @@ SUBROUTINE mixed_cdft_init_structures(force_env, force_env_qs, mixed_env, mixed_ NULLIFY (qs_kind_set) CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set) CALL create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, & - force_env%para_env, force_env_section) + force_env%para_env, force_env_section, silent=.FALSE.) mixed_cdft%qs_kind_set => qs_kind_set DEALLOCATE (i_force_eval) CALL section_vals_release(force_env_section) @@ -990,7 +990,7 @@ SUBROUTINE mixed_cdft_init_structures(force_env, force_env_qs, mixed_env, mixed_ NULLIFY (qs_kind_set) CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set) CALL create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, & - force_env%para_env, force_env_section) + force_env%para_env, force_env_section, silent=.FALSE.) mixed_cdft%qs_kind_set => qs_kind_set DEALLOCATE (i_force_eval) CALL section_vals_release(force_env_section) diff --git a/src/qs_eht_guess.F b/src/qs_eht_guess.F new file mode 100644 index 0000000000..9ed1b813a7 --- /dev/null +++ b/src/qs_eht_guess.F @@ -0,0 +1,238 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Generate an initial guess (dm and orb) from EHT calculation +! ************************************************************************************************** +MODULE qs_eht_guess + USE basis_set_types, ONLY: gto_basis_set_p_type + USE cp_blacs_env, ONLY: cp_blacs_env_type + USE cp_dbcsr_api, ONLY: dbcsr_create,& + dbcsr_desymmetrize,& + dbcsr_get_info,& + dbcsr_p_type,& + dbcsr_release,& + dbcsr_type,& + dbcsr_type_no_symmetry + USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& + cp_dbcsr_sm_fm_multiply,& + dbcsr_deallocate_matrix_set + USE cp_fm_basic_linalg, ONLY: cp_fm_invert + USE cp_fm_diag, ONLY: cp_fm_geeig + USE cp_fm_struct, ONLY: cp_fm_struct_create,& + cp_fm_struct_release,& + cp_fm_struct_type + USE cp_fm_types, ONLY: cp_fm_create,& + cp_fm_release,& + cp_fm_to_fm,& + cp_fm_type + USE cp_subsys_types, ONLY: cp_subsys_type + USE input_constants, ONLY: do_method_xtb + USE input_section_types, ONLY: section_vals_duplicate,& + section_vals_get_subs_vals,& + section_vals_release,& + section_vals_type,& + section_vals_val_set + USE kinds, ONLY: dp + USE message_passing, ONLY: mp_para_env_type + USE parallel_gemm_api, ONLY: parallel_gemm + USE qs_energy_init, ONLY: qs_energies_init + USE qs_environment, ONLY: qs_init + USE qs_environment_methods, ONLY: qs_env_rebuild_pw_env + USE qs_environment_types, ONLY: get_qs_env,& + qs_env_create,& + qs_env_release,& + qs_environment_type + USE qs_integral_utils, ONLY: basis_set_list_setup + USE qs_kind_types, ONLY: qs_kind_type + USE qs_ks_types, ONLY: qs_ks_env_type + USE qs_mo_occupation, ONLY: set_mo_occupation + USE qs_mo_types, ONLY: get_mo_set,& + mo_set_type + USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type + USE qs_overlap, ONLY: build_overlap_matrix_simple + USE xtb_ks_matrix, ONLY: build_xtb_ks_matrix +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_eht_guess' + + PUBLIC :: calculate_eht_guess + +! ************************************************************************************************** + +CONTAINS + +! ************************************************************************************************** +!> \brief EHT MO guess calclulation +!> \param qs_env ... +!> \param mo_array ... +! ************************************************************************************************** + SUBROUTINE calculate_eht_guess(qs_env, mo_array) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array + + CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_eht_guess' + + INTEGER :: handle, ispin, nao, nbas, neeht, neorb, & + nkind, nmo, nspins, zero + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues + REAL(KIND=dp), DIMENSION(:), POINTER :: eigval + TYPE(cp_blacs_env_type), POINTER :: blacs_env + TYPE(cp_fm_struct_type), POINTER :: mstruct_ee, mstruct_oe, mstruct_oo + TYPE(cp_fm_type) :: fmksmat, fmorb, fmscr, fmsmat, fmvec, & + fmwork, sfull, sinv + TYPE(cp_fm_type), POINTER :: mo_coeff + TYPE(cp_subsys_type), POINTER :: cp_subsys + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, matrix_s, matrix_t, smat + TYPE(dbcsr_type) :: tempmat, tmat + TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_nl + TYPE(qs_environment_type), POINTER :: eht_env + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(section_vals_type), POINTER :: dft_section, eht_force_env_section, & + force_env_section, qs_section, & + subsys_section, xtb_section + + CALL timeset(routineN, handle) + + NULLIFY (subsys_section) + CALL get_qs_env(qs_env, & + ks_env=ks_env, & + para_env=para_env, & + input=force_env_section, & + cp_subsys=cp_subsys) + NULLIFY (eht_force_env_section) + CALL section_vals_duplicate(force_env_section, eht_force_env_section) + dft_section => section_vals_get_subs_vals(eht_force_env_section, "DFT") + qs_section => section_vals_get_subs_vals(dft_section, "QS") + CALL section_vals_val_set(qs_section, "METHOD", i_val=do_method_xtb) + xtb_section => section_vals_get_subs_vals(qs_section, "xTB") + zero = 0 + CALL section_vals_val_set(xtb_section, "GFN_TYPE", i_val=zero) + ! + ALLOCATE (eht_env) + CALL qs_env_create(eht_env) + CALL qs_init(eht_env, para_env, cp_subsys=cp_subsys, & + force_env_section=eht_force_env_section, & + subsys_section=subsys_section, & + use_motion_section=.FALSE., silent=.TRUE.) + ! + CALL get_qs_env(qs_env, nelectron_total=neorb) + CALL get_qs_env(eht_env, nelectron_total=neeht) + IF (neorb /= neeht) THEN + CPWARN("EHT has different number of electrons than calculation method.") + CPABORT("EHT Initial Guess") + END IF + ! + CALL qs_env_rebuild_pw_env(eht_env) + CALL qs_energies_init(eht_env, calc_forces=.FALSE.) + CALL build_xtb_ks_matrix(eht_env, .FALSE., .FALSE.) + ! + CALL get_qs_env(eht_env, & + matrix_s=smat, matrix_ks=ksmat) + nspins = SIZE(ksmat, 1) + CALL get_qs_env(eht_env, para_env=para_env, blacs_env=blacs_env) + CALL dbcsr_get_info(smat(1)%matrix, nfullrows_total=nao) + CALL cp_fm_struct_create(fmstruct=mstruct_ee, context=blacs_env, & + nrow_global=nao, ncol_global=nao) + CALL cp_fm_create(fmksmat, mstruct_ee) + CALL cp_fm_create(fmsmat, mstruct_ee) + CALL cp_fm_create(fmvec, mstruct_ee) + CALL cp_fm_create(fmwork, mstruct_ee) + ALLOCATE (eigenvalues(nao)) + + ! DBCSR matrix + CALL dbcsr_create(tempmat, template=smat(1)%matrix, matrix_type=dbcsr_type_no_symmetry) + + ! transfer to FM + CALL dbcsr_desymmetrize(smat(1)%matrix, tempmat) + CALL copy_dbcsr_to_fm(tempmat, fmsmat) + + !SINV of origianl basis + CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env) + CALL get_qs_env(qs_env, matrix_s=matrix_s) + CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nbas) + CALL dbcsr_create(tmat, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry) + CALL cp_fm_struct_create(fmstruct=mstruct_oo, context=blacs_env, & + nrow_global=nbas, ncol_global=nbas) + CALL cp_fm_create(sfull, mstruct_oo) + CALL cp_fm_create(sinv, mstruct_oo) + CALL dbcsr_desymmetrize(matrix_s(1)%matrix, tmat) + CALL copy_dbcsr_to_fm(tmat, sfull) + CALL cp_fm_invert(sfull, sinv) + CALL dbcsr_release(tmat) + CALL cp_fm_release(sfull) + !TMAT(bas1, bas2) + CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, sab_all=sab_nl, nkind=nkind) + IF (.NOT. ASSOCIATED(sab_nl)) THEN + CPWARN("Full neighborlist not available for this method. EHT initial guess not possible.") + CPABORT("EHT Initial Guess") + END IF + ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind)) + CALL basis_set_list_setup(basis_set_list_a, "ORB", qs_kind_set) + CALL get_qs_env(eht_env, qs_kind_set=qs_kind_set) + CALL basis_set_list_setup(basis_set_list_b, "ORB", qs_kind_set) + ! + NULLIFY (matrix_t) + CALL build_overlap_matrix_simple(ks_env, matrix_t, & + basis_set_list_a, basis_set_list_b, sab_nl) + DEALLOCATE (basis_set_list_a, basis_set_list_b) + + ! KS matrix is not spin dependent! + CALL dbcsr_desymmetrize(ksmat(1)%matrix, tempmat) + CALL copy_dbcsr_to_fm(tempmat, fmksmat) + ! diagonalize + CALL cp_fm_geeig(fmksmat, fmsmat, fmvec, eigenvalues, fmwork) + ! Sinv*T*d + CALL cp_fm_struct_create(fmstruct=mstruct_oe, context=blacs_env, & + nrow_global=nbas, ncol_global=nao) + CALL cp_fm_create(fmscr, mstruct_oe) + CALL cp_fm_create(fmorb, mstruct_oe) + CALL cp_dbcsr_sm_fm_multiply(matrix_t(1)%matrix, fmvec, fmscr, ncol=nao) + CALL parallel_gemm('N', 'N', nbas, nao, nbas, 1.0_dp, sinv, fmscr, 0.0_dp, fmorb) + ! + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff, nmo=nmo) + CALL cp_fm_to_fm(fmorb, mo_coeff, nmo, 1, 1) + NULLIFY (eigval) + CALL get_mo_set(mo_set=mo_array(ispin), eigenvalues=eigval) + IF (ASSOCIATED(eigval)) THEN + eigval(1:nmo) = eigenvalues(1:nmo) + END IF + END DO + CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear) + + DEALLOCATE (eigenvalues) + CALL dbcsr_release(tempmat) + CALL dbcsr_deallocate_matrix_set(matrix_t) + CALL cp_fm_release(fmksmat) + CALL cp_fm_release(fmsmat) + CALL cp_fm_release(fmvec) + CALL cp_fm_release(fmwork) + CALL cp_fm_release(fmscr) + CALL cp_fm_release(fmorb) + CALL cp_fm_release(sinv) + CALL cp_fm_struct_release(mstruct_ee) + CALL cp_fm_struct_release(mstruct_oe) + CALL cp_fm_struct_release(mstruct_oo) + ! + CALL qs_env_release(eht_env) + DEALLOCATE (eht_env) + CALL section_vals_release(eht_force_env_section) + + CALL timestop(handle) + + END SUBROUTINE calculate_eht_guess + +END MODULE qs_eht_guess diff --git a/src/qs_environment.F b/src/qs_environment.F index 53cc00d40c..c3bf2b4098 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -58,7 +58,6 @@ MODULE qs_environment USE distribution_1d_types, ONLY: distribution_1d_release,& distribution_1d_type USE distribution_methods, ONLY: distribute_molecules_1d - USE dm_ls_scf_create, ONLY: ls_scf_create USE ec_env_types, ONLY: energy_correction_type USE ec_environment, ONLY: ec_env_create,& ec_write_input @@ -237,10 +236,12 @@ MODULE qs_environment !> \param force_env_section ... !> \param subsys_section ... !> \param use_motion_section ... +!> \param silent ... !> \author Creation (22.05.2000,MK) ! ************************************************************************************************** SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_env, cell, cell_ref, & - qmmm, qmmm_env_qm, force_env_section, subsys_section, use_motion_section) + qmmm, qmmm_env_qm, force_env_section, subsys_section, & + use_motion_section, silent) TYPE(qs_environment_type), POINTER :: qs_env TYPE(mp_para_env_type), POINTER :: para_env @@ -253,12 +254,13 @@ SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_en TYPE(qmmm_env_qm_type), OPTIONAL, POINTER :: qmmm_env_qm TYPE(section_vals_type), POINTER :: force_env_section, subsys_section LOGICAL, INTENT(IN) :: use_motion_section + LOGICAL, INTENT(IN), OPTIONAL :: silent CHARACTER(LEN=default_string_length) :: basis_type INTEGER :: ikind, method_id, nelectron_total, & nkind, nkp_grid(3) LOGICAL :: do_admm_rpa, do_ec_hfx, do_et, do_exx, do_hfx, do_kpoints, is_identical, is_semi, & - mp2_present, my_qmmm, qmmm_decoupl, same_except_frac, silent, use_ref_cell + mp2_present, my_qmmm, qmmm_decoupl, same_except_frac, use_ref_cell REAL(KIND=dp), DIMENSION(:, :), POINTER :: rtmat TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: my_cell, my_cell_ref @@ -324,7 +326,7 @@ SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_en use_motion_section=use_motion_section, & root_section=root_section, & cp_subsys=cp_subsys, cell=cell, cell_ref=cell_ref, & - elkind=is_semi) + elkind=is_semi, silent=silent) ALLOCATE (ks_env) CALL qs_ks_env_create(ks_env) @@ -355,12 +357,10 @@ SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_en ! kpoints IF (PRESENT(kpoint_env)) THEN - silent = .TRUE. kpoints => kpoint_env CALL set_qs_env(qs_env=qs_env, kpoints=kpoints) CALL kpoint_initialize(kpoints, particle_set, my_cell) ELSE - silent = .FALSE. NULLIFY (kpoints) CALL kpoint_create(kpoints) CALL set_qs_env(qs_env=qs_env, kpoints=kpoints) @@ -483,10 +483,6 @@ SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_en CALL transport_env_create(qs_env) END IF - IF (dft_control%qs_control%do_ls_scf) THEN - CALL ls_scf_create(qs_env) - END IF - CALL get_qs_env(qs_env, harris_env=harris_env) IF (qs_env%harris_method) THEN ! initialize the Harris input density and potential integrals @@ -835,9 +831,11 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD") print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION") IF (gfn_type == 0) THEN - CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, pset="EEQ") + CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, & + silent=silent, pset="EEQ") ELSE - CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat) + CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, & + silent=silent) END IF ALLOCATE (ewald_pw) CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section) diff --git a/src/qs_gamma2kp.F b/src/qs_gamma2kp.F index d84caeadc6..d9cbadc24e 100644 --- a/src/qs_gamma2kp.F +++ b/src/qs_gamma2kp.F @@ -107,7 +107,7 @@ SUBROUTINE create_kp_from_gamma(qs_env, qs_env_kp, with_xc_terms) CALL qs_env_create(qs_env_kp) CALL qs_init(qs_env_kp, para_env, cp_subsys=cp_subsys, kpoint_env=kpoint, & force_env_section=force_env_section, subsys_section=subsys_section, & - use_motion_section=.FALSE.) + use_motion_section=.FALSE., silent=.TRUE.) CALL get_qs_env(qs_env_kp, scf_control=scf_control) scf_control%density_guess = atomic_guess diff --git a/src/qs_initial_guess.F b/src/qs_initial_guess.F index 2194094aca..dfbdeeea51 100644 --- a/src/qs_initial_guess.F +++ b/src/qs_initial_guess.F @@ -46,14 +46,9 @@ MODULE qs_initial_guess gth_potential_type,& sgp_potential_type USE hfx_types, ONLY: hfx_type - USE input_constants, ONLY: atomic_guess,& - core_guess,& - history_guess,& - mopac_guess,& - no_guess,& - random_guess,& - restart_guess,& - sparse_guess + USE input_constants, ONLY: & + atomic_guess, core_guess, eht_guess, history_guess, mopac_guess, no_guess, random_guess, & + restart_guess, sparse_guess USE input_cp2k_hfx, ONLY: ri_mo USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type,& @@ -68,6 +63,7 @@ MODULE qs_initial_guess 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_eht_guess, ONLY: calculate_eht_guess USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_kind_types, ONLY: get_qs_kind,& @@ -275,6 +271,7 @@ SUBROUTINE calculate_first_density_matrix(scf_env, qs_env) (density_guess == atomic_guess) .OR. & (density_guess == core_guess) .OR. & (density_guess == mopac_guess) .OR. & + (density_guess == eht_guess) .OR. & (density_guess == sparse_guess) .OR. & (((density_guess == restart_guess) .OR. & (density_guess == history_guess)) .AND. & @@ -1032,7 +1029,8 @@ SUBROUTINE calculate_first_density_matrix(scf_env, qs_env) IF (density_guess == mopac_guess) THEN CALL calculate_mopac_dm(p_rmpv, s_sparse(1)%matrix, has_unit_metric, dft_control, & - particle_set, atomic_kind_set, qs_kind_set, nspin, nelectron_spin, para_env) + particle_set, atomic_kind_set, qs_kind_set, & + nspin, nelectron_spin, para_env) DO ispin = 1, nspin ! The orbital transformation method (OT) requires not only an @@ -1075,7 +1073,16 @@ SUBROUTINE calculate_first_density_matrix(scf_env, qs_env) did_guess = .TRUE. END IF -! switch_surf_dip [SGh] + ! + ! EHT guess (gfn0-xTB) + IF (density_guess == eht_guess) THEN + CALL calculate_eht_guess(qs_env, mo_array) + DO ispin = 1, nspin + CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix) + END DO + did_guess = .TRUE. + END IF + ! switch_surf_dip [SGh] IF (dft_control%switch_surf_dip) THEN DO ispin = 1, nspin CALL reassign_allocated_mos(mos_last_converged(ispin), & diff --git a/src/qs_kind_types.F b/src/qs_kind_types.F index 0eb6eef42a..200d153bd8 100644 --- a/src/qs_kind_types.F +++ b/src/qs_kind_types.F @@ -1493,12 +1493,14 @@ END SUBROUTINE init_gapw_nlcc !> \param force_env_section ... !> \param no_fail ... !> \param method_id ... +!> \param silent ... !> \par History !> - Creation (09.02.2002,MK) !> - 20.09.2002,gt: adapted for POL/KG use (elp_potential) !> - 05.03.2010: split elp_potential into fist_potential and kg_potential ! ************************************************************************************************** - SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_fail, method_id) + SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, & + no_fail, method_id, silent) TYPE(qs_kind_type), INTENT(INOUT) :: qs_kind TYPE(section_vals_type), POINTER :: kind_section @@ -1506,6 +1508,7 @@ SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_f TYPE(section_vals_type), POINTER :: force_env_section LOGICAL, INTENT(IN) :: no_fail INTEGER, INTENT(IN) :: method_id + LOGICAL, INTENT(IN) :: silent CHARACTER(LEN=*), PARAMETER :: routineN = 'read_qs_kind' INTEGER, PARAMETER :: maxbas = 20 @@ -2117,13 +2120,13 @@ SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_f zeff=qs_kind%se_parameter%zeff, zeff_correction=zeff_correction) qs_kind%se_parameter%zeff = qs_kind%se_parameter%zeff - zeff_correction - check = ((potential_name /= '') .OR. explicit_potential) + check = ((potential_name /= '') .OR. explicit_potential) .AND. .NOT. silent IF (check) & CALL cp_warn(__LOCATION__, & "Information provided in the input file regarding POTENTIAL for KIND <"// & TRIM(qs_kind%name)//"> will be ignored!") - check = ((k_rep > 0) .OR. explicit_basis) + check = ((k_rep > 0) .OR. explicit_basis) .AND. .NOT. silent IF (check) & CALL cp_warn(__LOCATION__, & "Information provided in the input file regarding BASIS for KIND <"// & @@ -2150,13 +2153,13 @@ SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_f zeff_correction=0.0_dp) END IF - check = ((potential_name /= '') .OR. explicit_potential) + check = ((potential_name /= '') .OR. explicit_potential) .AND. .NOT. silent IF (check) & CALL cp_warn(__LOCATION__, & "Information provided in the input file regarding POTENTIAL for KIND <"// & TRIM(qs_kind%name)//"> will be ignored!") - check = ((k_rep > 0) .OR. explicit_basis) + check = ((k_rep > 0) .OR. explicit_basis) .AND. .NOT. silent IF (check) & CALL cp_warn(__LOCATION__, & "Information provided in the input file regarding BASIS for KIND <"// & @@ -2183,13 +2186,13 @@ SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_f zeff_correction=0.0_dp) END IF - check = ((potential_name /= '') .OR. explicit_potential) + check = ((potential_name /= '') .OR. explicit_potential) .AND. .NOT. silent IF (check) & CALL cp_warn(__LOCATION__, & "Information provided in the input file regarding POTENTIAL for KIND <"// & TRIM(qs_kind%name)//"> will be ignored!") - check = ((k_rep > 0) .OR. explicit_basis) + check = ((k_rep > 0) .OR. explicit_basis) .AND. .NOT. silent IF (check) & CALL cp_warn(__LOCATION__, & "Information provided in the input file regarding BASIS for KIND <"// & @@ -2387,7 +2390,7 @@ SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_f CPABORT("ECPs are only well tested in their semi-local form") CALL get_qs_kind(qs_kind, basis_set=orb_basis_set) CALL sgp_construction(sgp_pot=sgppot, ecp_pot=ecppot, orb_basis=orb_basis_set, error=error) - IF (iounit > 0) THEN + IF (iounit > 0 .AND. .NOT. silent) THEN WRITE (iounit, "(/,T2,'PP Transformation for ',A)") TRIM(ecppot%pname) IF (sgppot%has_local) THEN WRITE (iounit, "(T8,'Accuracy for local part:',T41,F10.3,'%',T61,F20.12)") error(4), error(1) @@ -2447,7 +2450,7 @@ SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_f zeff=upfpot%zion, z=z, has_nlcc=upfpot%core_correction) ! convert pp CALL sgp_construction(sgp_pot=sgppot, upf_pot=upfpot, error=error) - IF (iounit > 0) THEN + IF (iounit > 0 .AND. .NOT. silent) THEN WRITE (iounit, "(/,T2,'PP Transformation for ',A)") TRIM(upfpot%pname) IF (sgppot%has_local) THEN WRITE (iounit, "(T8,'Accuracy for local part:',T61,F20.12)") error(1) @@ -2617,14 +2620,17 @@ FUNCTION parse_valence_electrons(string) RESULT(n) !> \param kind_section ... !> \param para_env ... !> \param force_env_section ... +!> \param silent ... ! ************************************************************************************************** - SUBROUTINE create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, para_env, force_env_section) + SUBROUTINE create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, para_env, & + force_env_section, silent) TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(section_vals_type), POINTER :: kind_section TYPE(mp_para_env_type), POINTER :: para_env TYPE(section_vals_type), POINTER :: force_env_section + LOGICAL, INTENT(IN) :: silent CHARACTER(len=*), PARAMETER :: routineN = 'create_qs_kind_set' @@ -2664,7 +2670,8 @@ SUBROUTINE create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, para_e qs_kind_set(ikind)%name = atomic_kind_set(ikind)%name qs_kind_set(ikind)%element_symbol = atomic_kind_set(ikind)%element_symbol qs_kind_set(ikind)%natom = atomic_kind_set(ikind)%natom - CALL read_qs_kind(qs_kind_set(ikind), kind_section, para_env, force_env_section, no_fail, qs_method) + CALL read_qs_kind(qs_kind_set(ikind), kind_section, para_env, force_env_section, & + no_fail, qs_method, silent) END DO CALL timestop(handle) diff --git a/src/qs_subsys_methods.F b/src/qs_subsys_methods.F index a1683d8922..27a8e04d68 100644 --- a/src/qs_subsys_methods.F +++ b/src/qs_subsys_methods.F @@ -66,9 +66,10 @@ MODULE qs_subsys_methods !> \param cell ... !> \param cell_ref ... !> \param elkind ... +!> \param silent ... ! ************************************************************************************************** SUBROUTINE qs_subsys_create(subsys, para_env, root_section, force_env_section, subsys_section, & - use_motion_section, cp_subsys, cell, cell_ref, elkind) + use_motion_section, cp_subsys, cell, cell_ref, elkind, silent) TYPE(qs_subsys_type), INTENT(OUT) :: subsys TYPE(mp_para_env_type), POINTER :: para_env TYPE(section_vals_type), OPTIONAL, POINTER :: root_section @@ -76,9 +77,9 @@ SUBROUTINE qs_subsys_create(subsys, para_env, root_section, force_env_section, s LOGICAL, INTENT(IN) :: use_motion_section TYPE(cp_subsys_type), OPTIONAL, POINTER :: cp_subsys TYPE(cell_type), OPTIONAL, POINTER :: cell, cell_ref - LOGICAL, INTENT(IN), OPTIONAL :: elkind + LOGICAL, INTENT(IN), OPTIONAL :: elkind, silent - LOGICAL :: use_ref_cell + LOGICAL :: be_silent, use_ref_cell TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: my_cell, my_cell_ref TYPE(cp_subsys_type), POINTER :: my_cp_subsys @@ -87,6 +88,8 @@ SUBROUTINE qs_subsys_create(subsys, para_env, root_section, force_env_section, s NULLIFY (atomic_kind_set, qs_kind_set, cell_section, kind_section, my_cell, my_cell_ref, my_cp_subsys) + be_silent = .FALSE. + IF (PRESENT(silent)) be_silent = silent ! create cp_subsys IF (PRESENT(cp_subsys)) THEN my_cp_subsys => cp_subsys @@ -125,7 +128,7 @@ SUBROUTINE qs_subsys_create(subsys, para_env, root_section, force_env_section, s CALL cp_subsys_get(my_cp_subsys, atomic_kind_set=atomic_kind_set) kind_section => section_vals_get_subs_vals(subsys_section, "KIND") CALL create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, & - para_env, force_env_section) + para_env, force_env_section, be_silent) CALL num_ao_el_per_molecule(my_cp_subsys%molecule_kinds%els, & qs_kind_set) @@ -139,6 +142,7 @@ SUBROUTINE qs_subsys_create(subsys, para_env, root_section, force_env_section, s IF (.NOT. PRESENT(cell)) CALL cell_release(my_cell) IF (.NOT. PRESENT(cell_ref)) CALL cell_release(my_cell_ref) IF (.NOT. PRESENT(cp_subsys)) CALL cp_subsys_release(my_cp_subsys) + END SUBROUTINE qs_subsys_create ! ************************************************************************************************** diff --git a/tests/QS/regtest-eht-guess/TEST_FILES b/tests/QS/regtest-eht-guess/TEST_FILES new file mode 100644 index 0000000000..dfc607a4fb --- /dev/null +++ b/tests/QS/regtest-eht-guess/TEST_FILES @@ -0,0 +1,9 @@ +# runs are executed in the same order as in this file +# the second field tells which test should be run in order to compare with the last available output +# e.g. 0 means do not compare anything, running is enough +# 1 compares the last total energy in the file +# for details see cp2k/tools/do_regtest +#QS +h2o_1.inp 1 1e-12 -16.83566845657052 +h2o_2.inp 1 1e-12 -16.41572764902628 +#EOF diff --git a/tests/QS/regtest-eht-guess/h2o_1.inp b/tests/QS/regtest-eht-guess/h2o_1.inp new file mode 100644 index 0000000000..59736cfb1d --- /dev/null +++ b/tests/QS/regtest-eht-guess/h2o_1.inp @@ -0,0 +1,46 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT H2O + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.0E-8 + &END QS + &SCF + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 1 + SCF_GUESS EHT + &END SCF + &XC + &XC_FUNCTIONAL Pade + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + &END CELL + &COORD + O 0.000000 0.000000 -0.065587 + H 0.000000 -0.757136 0.520545 + H 0.000000 0.757136 0.520545 + &END COORD + &KIND H + BASIS_SET DZV-GTH-PADE + POTENTIAL GTH-PADE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH-PADE + POTENTIAL GTH-PADE-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-eht-guess/h2o_2.inp b/tests/QS/regtest-eht-guess/h2o_2.inp new file mode 100644 index 0000000000..d27cd442fb --- /dev/null +++ b/tests/QS/regtest-eht-guess/h2o_2.inp @@ -0,0 +1,46 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT H2O + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH + POTENTIAL_FILE_NAME POTENTIAL_UZH + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.0E-8 + &END QS + &SCF + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 1 + SCF_GUESS EHT + &END SCF + &XC + &XC_FUNCTIONAL Pade + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + &END CELL + &COORD + O 0.000000 0.000000 -0.065587 + H 0.000000 -0.757136 0.520545 + H 0.000000 0.757136 0.520545 + &END COORD + &KIND H + BASIS_SET TZV2P-MOLOPT-GGA-GTH-q1 + POTENTIAL GTH-GGA-q1 + &END KIND + &KIND O + BASIS_SET TZV2P-MOLOPT-GGA-GTH-q6 + POTENTIAL GTH-GGA-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index a2aff14f7f..23e5e40379 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -372,3 +372,4 @@ Fist/regtest-deepmd deepmd QS/regtest-elpa elpa mpiranks==1||mpiranks%2==0 Fist/regtest-plumed2 plumed2 Fist/regtest-quip quip +QS/regtest-eht-guess libdftd4