diff --git a/src/cp_control_types.F b/src/cp_control_types.F index b114546385..d57974958f 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -612,7 +612,8 @@ MODULE cp_control_types INTEGER :: surf_normal = 0 REAL(KIND=dp) :: dist_edge = 0.0_dp, & gau_a = 0.0_dp, & - gau_c = 0.0_dp + gau_c = 0.0_dp, & + charge = 0.0_dp END TYPE pcc_control_type ! ************************************************************************************************** diff --git a/src/gce_methods.F b/src/gce_methods.F index 1d2b3547fc..c70efed5b8 100644 --- a/src/gce_methods.F +++ b/src/gce_methods.F @@ -212,6 +212,8 @@ CONTAINS tot_charge = pw_integrate_function(rho_tot_gspace, isign=-1) + pcc_env%charge = tot_charge + dz = pcc_density_r%pw_grid%dr(c) ucenter = pcc_density_r%pw_grid%bounds(2, c) - INT(pcc_env%dist_edge/dz) diff --git a/src/input_cp2k_scf.F b/src/input_cp2k_scf.F index 322bdf01e7..2162acfc57 100644 --- a/src/input_cp2k_scf.F +++ b/src/input_cp2k_scf.F @@ -13,14 +13,9 @@ !> \author fawzi ! ************************************************************************************************** MODULE input_cp2k_scf - USE bibliography, ONLY: Becke1988b,& - Holmberg2017,& - Holmberg2018,& - Schiffmann2015,& - Stewart1982,& - VandeVondele2003,& - VandeVondele2005a,& - Weber2008 + USE bibliography, ONLY: & + Becke1988b, Chai2024a, Holmberg2017, Holmberg2018, Schiffmann2015, Stewart1982, & + VandeVondele2003, VandeVondele2005a, Weber2008 USE cp_output_handling, ONLY: add_last_numeric,& cp_print_key_section_create,& high_print_level,& @@ -2048,7 +2043,7 @@ CONTAINS END SUBROUTINE create_cdft_opt_section ! ************************************************************************************************** -!> \brief Create CP2K input section for the grand canonical ensemble DFT +!> \brief Create CP2K input section for the grand canonical SCF !> \param section ... !> \date !> \author Ziwei Chai @@ -2064,7 +2059,26 @@ CONTAINS CALL section_create(section, __LOCATION__, & name="GCE", & - description="Define the settings of the grand canonical ensemble DFT", & + description="Enables grand canonical SCF calculations for "// & + "slab models with a target work function specified by the "// & + "user. The current implementation is intended for three "// & + "dimensional periodic calculations and should be used "// & + "together with the planar counter charge model. The planar "// & + "counter charge screens the net charge of the DFT slab, so "// & + "that the work function can be obtained from the planar "// & + "averaged Hartree potential in the implicit solvent region "// & + "or vacuum region. During the SCF iterations, the work "// & + "function is mixed toward the target value. "// & + "The reported GCE charges during SCF iterations are the "// & + "current net charge of the DFT system excluding the planar "// & + "counter charge. A positive value means electron deficiency. "// & + "The input CHARGE value in the DFT section defines the "// & + "reference charge state. The quantity N_e^(abs,alpha) is "// & + "evaluated as the input CHARGE value minus the reported GCE "// & + "charge. The reported GCE free energy is the sum of the DFT "// & + "total energy and the product of the target work function "// & + "and N_e^(abs,alpha).", & + citations=[Chai2024a], & n_keywords=3, & n_subsections=0, & repeats=.FALSE.) @@ -2073,7 +2087,7 @@ CONTAINS CALL keyword_create(keyword, __LOCATION__, & name="_SECTION_PARAMETERS_", & - description="Controls the activation of grand canonical ensemble DFT", & + description="Controls the activation of grand canonical SCF", & usage="&GCE ON", & default_l_val=.FALSE., & lone_keyword_l_val=.TRUE.) diff --git a/src/qs_environment.F b/src/qs_environment.F index 45d38fdb19..bf7ae4766d 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -90,16 +90,16 @@ MODULE qs_environment USE hfx_types, ONLY: compare_hfx_sections,& hfx_create USE input_constants, ONLY: & - debug_run, dispersion_d2, dispersion_d3, dispersion_d3bj, do_et_ddapc, do_method_am1, & - do_method_dftb, do_method_gapw, do_method_gapw_xc, do_method_gpw, do_method_lrigpw, & - do_method_mndo, do_method_mndod, do_method_ofgpw, do_method_pdg, do_method_pm3, & - do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rigpw, do_method_rm1, & - do_method_xtb, do_qmmm_gauss, do_qmmm_swave, general_roks, gfn1xtb, hden_atomic, & - kg_tnadd_embed_ri, linear_response_run, rel_none, rel_trans_atom, smear_fermi_dirac, & - tblite_scc_mixer_tblite, tddfpt_kernel_none, vdw_pairpot_dftd2, vdw_pairpot_dftd3, & - vdw_pairpot_dftd3bj, vdw_pairpot_dftd4, wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, & - wfi_use_prev_wf_method_nr, xc_vdw_fun_none, xc_vdw_fun_nonloc, xc_vdw_fun_pairpot, & - xtb_vdw_type_d3, xtb_vdw_type_d4, xtb_vdw_type_none + debug_run, diag_ot, dispersion_d2, dispersion_d3, dispersion_d3bj, do_et_ddapc, & + do_method_am1, do_method_dftb, do_method_gapw, do_method_gapw_xc, do_method_gpw, & + do_method_lrigpw, do_method_mndo, do_method_mndod, do_method_ofgpw, do_method_pdg, & + do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rigpw, & + do_method_rm1, do_method_xtb, do_qmmm_gauss, do_qmmm_swave, general_roks, gfn1xtb, & + hden_atomic, kg_tnadd_embed_ri, linear_response_run, rel_none, rel_trans_atom, & + smear_fermi_dirac, tblite_scc_mixer_tblite, tddfpt_kernel_none, vdw_pairpot_dftd2, & + vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, vdw_pairpot_dftd4, wfi_linear_ps_method_nr, & + wfi_linear_wf_method_nr, wfi_use_prev_wf_method_nr, xc_vdw_fun_none, xc_vdw_fun_nonloc, & + xc_vdw_fun_pairpot, xtb_vdw_type_d3, xtb_vdw_type_d4, xtb_vdw_type_none USE input_section_types, ONLY: section_get_ival,& section_get_ivals,& section_vals_get,& @@ -1603,6 +1603,23 @@ CONTAINS END IF CALL scf_c_create(scf_control) CALL scf_c_read_parameters(scf_control, dft_section) + IF (scf_control%gce%do_gce) THEN + IF (.NOT. ALL(cell%perd == 1)) THEN + CPABORT("Grand canonical SCF is only implemented for 3D periodic calculations.") + END IF + IF (.NOT. scf_control%smear%do_smear) THEN + CPABORT("Grand canonical SCF requires smearing.") + END IF + IF (scf_control%smear%method /= smear_fermi_dirac) THEN + CPABORT("Grand canonical SCF is only implemented for Fermi-Dirac way of smearing.") + END IF + IF (scf_control%use_ot .OR. .NOT. scf_control%use_diag .OR. & + scf_control%diagonalization%method == diag_ot) THEN + CALL cp_abort(__LOCATION__, & + "Grand canonical SCF requires standard diagonalization. "// & + "It is not implemented with OT.") + END IF + END IF IF (.NOT. dft_control%qs_control%do_ls_scf) THEN SELECT CASE (dft_control%qs_control%method_id) CASE (do_method_dftb) diff --git a/src/qs_mo_occupation.F b/src/qs_mo_occupation.F index 0ad635f56a..38d815e707 100644 --- a/src/qs_mo_occupation.F +++ b/src/qs_mo_occupation.F @@ -29,9 +29,11 @@ MODULE qs_mo_occupation has_uniform_occupation,& mo_set_type,& set_mo_set - USE scf_control_types, ONLY: smear_type + USE scf_control_types, ONLY: gce_type,& + smear_type USE smearing_utils, ONLY: SmearFixed,& - SmearFixedDerivMV + SmearFixedDerivMV,& + SmearOcc USE util, ONLY: sort USE xas_env_types, ONLY: get_xas_env,& xas_environment_type @@ -57,14 +59,16 @@ CONTAINS !> !> \param mo_array ... !> \param smear ... +!> \param gce ... !> \date 10.03.2011 (MI) !> \author MI !> \version 1.0 ! ************************************************************************************************** - SUBROUTINE set_mo_occupation_3(mo_array, smear) + SUBROUTINE set_mo_occupation_3(mo_array, smear, gce) TYPE(mo_set_type), DIMENSION(2), INTENT(INOUT) :: mo_array TYPE(smear_type) :: smear + TYPE(gce_type), OPTIONAL, POINTER :: gce CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mo_occupation_3' @@ -73,7 +77,7 @@ CONTAINS lfomo_a, lfomo_b, nmo_a, nmo_b, & xas_estate INTEGER, ALLOCATABLE, DIMENSION(:) :: all_index - LOGICAL :: is_large + LOGICAL :: do_gce, is_large REAL(KIND=dp) :: all_nelec, kTS, mu, nelec_a, nelec_b, & occ_estate, smear_width REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: all_eigval, all_occ @@ -96,23 +100,11 @@ CONTAINS CALL sort(all_eigval, all_nmo, all_index) - xas_estate = -1 - occ_estate = 0.0_dp - - nelec_a = 0.0_dp - nelec_b = 0.0_dp - all_nelec = 0.0_dp - nelec_a = accurate_sum(occ_a(:)) - nelec_b = accurate_sum(occ_b(:)) - all_nelec = nelec_a + nelec_b - - DO i = 1, all_nmo - IF (all_index(i) <= nmo_a) THEN - all_occ(i) = occ_a(all_index(i)) - ELSE - all_occ(i) = occ_b(all_index(i) - nmo_a) - END IF - END DO + IF (PRESENT(gce)) THEN + do_gce = gce%do_gce + ELSE + do_gce = .FALSE. + END IF SELECT CASE (smear%method) CASE (smear_fermi_dirac) @@ -131,8 +123,32 @@ CONTAINS CPABORT("set_mo_occupation_3: unsupported smearing method") END SELECT - CALL SmearFixed(all_occ, mu, kTS, all_eigval, all_nelec, & - smear_width, 1._dp, smear%method, xas_estate, occ_estate) + IF (.NOT. do_gce) THEN + xas_estate = -1 + occ_estate = 0.0_dp + + nelec_a = 0.0_dp + nelec_b = 0.0_dp + all_nelec = 0.0_dp + nelec_a = accurate_sum(occ_a(:)) + nelec_b = accurate_sum(occ_b(:)) + all_nelec = nelec_a + nelec_b + + DO i = 1, all_nmo + IF (all_index(i) <= nmo_a) THEN + all_occ(i) = occ_a(all_index(i)) + ELSE + all_occ(i) = occ_b(all_index(i) - nmo_a) + END IF + END DO + + CALL SmearFixed(all_occ, mu, kTS, all_eigval, all_nelec, & + smear_width, 1._dp, smear%method, xas_estate, occ_estate) + ELSE + gce%prev_workfunction = gce%ref_esp - mo_array(1)%mu + mu = gce%ref_esp - ((1.0_dp - gce%mixing_coef)*gce%prev_workfunction + gce%mixing_coef*gce%target_workfunction) + CALL SmearOcc(all_occ, all_nelec, kTS, all_eigval, mu, smear_width, 1._dp, smear%method) + END IF is_large = ABS(all_occ(1) - 1.0_dp) > smear%eps_fermi_dirac ! this is not a real problem, but the smearing width might be a bit large @@ -143,10 +159,11 @@ CONTAINS CALL cp_warn(__LOCATION__, & TRIM(method_label)//" smearing includes the last MO => "// & "Add more MOs for proper smearing.") - - ! check that the total electron count is accurate - is_large = (ABS(all_nelec - accurate_sum(all_occ(:))) > smear%eps_fermi_dirac*all_nelec) - CPWARN_IF(is_large, "Total number of electrons is not accurate") + IF (.NOT. do_gce) THEN + ! check that the total electron count is accurate + is_large = (ABS(all_nelec - accurate_sum(all_occ(:))) > smear%eps_fermi_dirac*all_nelec) + CPWARN_IF(is_large, "Total number of electrons is not accurate") + END IF DO i = 1, all_nmo IF (all_index(i) <= nmo_a) THEN @@ -205,6 +222,7 @@ CONTAINS !> \param eval_deriv ... !> \param tot_zeff_corr ... !> \param probe ... +!> \param gce ... !> \date 25.01.2010 (MK) !> \par History !> 10.2019 Added functionality to adjust mo occupation if the core @@ -215,7 +233,7 @@ CONTAINS !> \author Matthias Krack (MK) !> \version 1.0 ! ************************************************************************************************** - SUBROUTINE set_mo_occupation_2(mo_array, smear, eval_deriv, tot_zeff_corr, probe) + SUBROUTINE set_mo_occupation_2(mo_array, smear, eval_deriv, tot_zeff_corr, probe, gce) TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mo_array TYPE(smear_type) :: smear @@ -223,6 +241,7 @@ CONTAINS REAL(KIND=dp), OPTIONAL :: tot_zeff_corr TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, & POINTER :: probe + TYPE(gce_type), OPTIONAL, POINTER :: gce CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mo_occupation_2' @@ -240,12 +259,20 @@ CONTAINS CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe) ELSE IF (PRESENT(eval_deriv)) THEN ! Change of MO occupancy to account for CORE_CORRECTION is not yet implemented - CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv) + IF (PRESENT(gce)) THEN + CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv, gce=gce) + ELSE + CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv) + END IF ELSE IF (PRESENT(tot_zeff_corr)) THEN CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr) ELSE - CALL set_mo_occupation_1(mo_array(1), smear=smear) + IF (PRESENT(gce)) THEN + CALL set_mo_occupation_1(mo_array(1), smear=smear, gce=gce) + ELSE + CALL set_mo_occupation_1(mo_array(1), smear=smear) + END IF END IF END IF CALL timestop(handle) @@ -269,7 +296,7 @@ CONTAINS END IF IF (smear%fixed_mag_mom /= -1.0_dp) THEN CPASSERT(.NOT. (PRESENT(eval_deriv))) - CALL set_mo_occupation_3(mo_array, smear=smear) + CALL set_mo_occupation_3(mo_array, smear=smear, gce=gce) CALL timestop(handle) RETURN END IF @@ -414,6 +441,7 @@ CONTAINS !> \param xas_env ... !> \param tot_zeff_corr ... !> \param probe ... +!> \param gce ... !> \date 17.04.2002 (v1.0), 26.08.2008 (v1.1) !> \par History !> 10.2019 Added functionality to adjust mo occupation if the core @@ -424,7 +452,7 @@ CONTAINS !> \author Matthias Krack !> \version 1.1 ! ************************************************************************************************** - SUBROUTINE set_mo_occupation_1(mo_set, smear, eval_deriv, xas_env, tot_zeff_corr, probe) + SUBROUTINE set_mo_occupation_1(mo_set, smear, eval_deriv, xas_env, tot_zeff_corr, probe, gce) TYPE(mo_set_type), INTENT(INOUT) :: mo_set TYPE(smear_type), OPTIONAL :: smear @@ -433,16 +461,18 @@ CONTAINS REAL(KIND=dp), OPTIONAL :: tot_zeff_corr TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, & POINTER :: probe + TYPE(gce_type), OPTIONAL, POINTER :: gce CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mo_occupation_1' CHARACTER(LEN=20) :: method_label - INTEGER :: handle, i_first, imo, ir, irmo, nmo, & + INTEGER :: handle, i, i_first, imo, ir, irmo, nmo, & nomo, xas_estate - LOGICAL :: equal_size, is_large + LOGICAL :: do_gce, equal_size, is_large REAL(KIND=dp) :: delectron, e1, e2, edelta, edist, & - el_count, nelec, occ_estate, & - total_zeff_corr, xas_nelectron + el_count, gce_mu, my_nelec, nelec, & + occ_estate, total_zeff_corr, & + xas_nelectron REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp_v CALL timeset(routineN, handle) @@ -472,6 +502,54 @@ CONTAINS is_large = ABS(el_count - xas_nelectron) > xas_nelectron*EPSILON(el_count) CPASSERT(.NOT. is_large) ELSE + IF (PRESENT(gce)) THEN + do_gce = gce%do_gce + ELSE + do_gce = .FALSE. + END IF + ! GCE workfunction (Fermi energy) mixing + IF (do_gce) THEN + IF (smear%method /= smear_fermi_dirac) THEN + CPABORT("Grand canonical ensemble DFT SCF now only support Fermi Dirac smearing.") + END IF + IF (gce%prev_workfunction < -1000.0_dp) THEN + my_nelec = REAL(mo_set%nelectron, dp) + CALL SmearFixed(mo_set%occupation_numbers, mo_set%mu, mo_set%kTS, mo_set%eigenvalues, my_nelec, & + smear%electronic_temperature, mo_set%maxocc, smear_fermi_dirac) + gce%prev_workfunction = -501.0_dp + ELSE IF (gce%prev_workfunction < -500.0_dp) THEN + my_nelec = REAL(mo_set%nelectron, dp) + CALL SmearFixed(mo_set%occupation_numbers, mo_set%mu, mo_set%kTS, mo_set%eigenvalues, my_nelec, & + smear%electronic_temperature, mo_set%maxocc, smear_fermi_dirac) + gce%prev_workfunction = gce%ref_esp - mo_set%mu + ELSE + gce_mu = gce%ref_esp - ((1.0_dp - gce%mixing_coef)*gce%prev_workfunction & + + gce%mixing_coef*gce%target_workfunction) + CALL SmearOcc(mo_set%occupation_numbers, my_nelec, mo_set%kTS, mo_set%eigenvalues, gce_mu, & + smear%electronic_temperature, mo_set%maxocc, smear_fermi_dirac) + mo_set%mu = gce_mu + gce%prev_workfunction = gce%ref_esp - mo_set%mu + is_large = ABS(MAXVAL(mo_set%occupation_numbers) - mo_set%maxocc) > smear%eps_fermi_dirac + CPWARN_IF(is_large, "Fermi-Dirac smearing includes the first MO") + END IF + DO i = 1, SIZE(mo_set%occupation_numbers) + IF (mo_set%occupation_numbers(i) < mo_set%maxocc) THEN + mo_set%lfomo = i + EXIT + END IF + END DO + DO i = SIZE(mo_set%occupation_numbers), 1, -1 + IF (mo_set%occupation_numbers(i) > smear%eps_fermi_dirac) THEN + mo_set%homo = i + EXIT + END IF + END DO + mo_set%uniform_occupation = .FALSE. + mo_set%n_el_f = my_nelec + CALL timestop(handle) + RETURN + END IF + IF (MODULO(mo_set%nelectron, INT(mo_set%maxocc)) == 0) THEN nomo = NINT(mo_set%nelectron/mo_set%maxocc) ! Initialize MO occupations diff --git a/src/qs_scf.F b/src/qs_scf.F index 4f9fd9ce7b..30865f773b 100644 --- a/src/qs_scf.F +++ b/src/qs_scf.F @@ -97,6 +97,7 @@ MODULE qs_scf USE message_passing, ONLY: mp_comm_type,& mp_para_env_type USE particle_types, ONLY: particle_type + USE physcon, ONLY: evolt USE preconditioner, ONLY: prepare_preconditioner,& restart_preconditioner USE pw_env_types, ONLY: pw_env_get,& @@ -160,6 +161,7 @@ MODULE qs_scf qs_scf_set_loop_flags USE qs_scf_output, ONLY: qs_scf_cdft_info,& qs_scf_cdft_initial_info,& + qs_scf_gce_info,& qs_scf_loop_info,& qs_scf_loop_print,& qs_scf_outer_loop_info,& @@ -434,6 +436,16 @@ CONTAINS output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", & extension=".scfLog") + IF (scf_control%gce%do_gce .AND. output_unit > 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T2,78('-'))") + WRITE (UNIT=output_unit, FMT="(T31,A)") & + "GRAND-CANONICAL SCF" + WRITE (UNIT=output_unit, FMT="(T20,A,F12.6,A)") & + "Target work function (TWF):", & + evolt*scf_control%gce%target_workfunction, " eV" + WRITE (UNIT=output_unit, FMT="(T2,78('-'))") + END IF + IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") & "SCF WAVEFUNCTION OPTIMIZATION" @@ -579,6 +591,10 @@ CONTAINS CALL qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy) + IF (scf_control%gce%do_gce) THEN + CALL qs_scf_gce_info(output_unit, qs_env, just_energy) + END IF + IF (.NOT. just_energy) energy%tot_old = energy%total ! check for external communicator and if the intermediate energy should be sent @@ -775,8 +791,14 @@ CONTAINS probe=dft_control%probe) END IF ELSE - CALL set_mo_occupation(mo_set=mos(ispin), & - smear=scf_control%smear) + IF (.NOT. scf_control%gce%do_gce) THEN + CALL set_mo_occupation(mo_set=mos(ispin), & + smear=scf_control%smear) + ELSE + CALL set_mo_occupation(mo_set=mos(ispin), & + smear=scf_control%smear, & + gce=scf_control%gce) + END IF END IF END IF END DO diff --git a/src/qs_scf_diagonalization.F b/src/qs_scf_diagonalization.F index 1b68fa757c..69a055808c 100644 --- a/src/qs_scf_diagonalization.F +++ b/src/qs_scf_diagonalization.F @@ -404,8 +404,14 @@ CONTAINS smear=scf_control%smear, & probe=probe) ELSE - CALL set_mo_occupation(mo_array=mos, & - smear=scf_control%smear) + IF (.NOT. scf_control%gce%do_gce) THEN + CALL set_mo_occupation(mo_array=mos, & + smear=scf_control%smear) + ELSE + CALL set_mo_occupation(mo_array=mos, & + smear=scf_control%smear, & + gce=scf_control%gce) + END IF END IF END IF @@ -1175,8 +1181,14 @@ CONTAINS !rotate the mos by the eigenvectors of C'HC CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, mo_coeff) - CALL set_mo_occupation(mo_set=mos(ispin), & - smear=scf_control%smear) + IF (.NOT. scf_control%gce%do_gce) THEN + CALL set_mo_occupation(mo_set=mos(ispin), & + smear=scf_control%smear) + ELSE + CALL set_mo_occupation(mo_set=mos(ispin), & + smear=scf_control%smear, & + gce=scf_control%gce) + END IF ! does not yet handle k-points CALL calculate_density_matrix(mos(ispin), & @@ -1387,8 +1399,14 @@ CONTAINS jacobi_threshold=scf_control%diagonalization%jacobi_threshold) END DO - CALL set_mo_occupation(mo_array=mos, & - smear=scf_control%smear) + IF (.NOT. scf_control%gce%do_gce) THEN + CALL set_mo_occupation(mo_array=mos, & + smear=scf_control%smear) + ELSE + CALL set_mo_occupation(mo_array=mos, & + smear=scf_control%smear, & + gce=scf_control%gce) + END IF DO ispin = 1, nspin ! does not yet handle k-points @@ -1956,8 +1974,14 @@ CONTAINS CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "SOLVE") END IF - CALL set_mo_occupation(mo_set=mos(ispin), & - smear=scf_control%smear) + IF (.NOT. scf_control%gce%do_gce) THEN + CALL set_mo_occupation(mo_set=mos(ispin), & + smear=scf_control%smear) + ELSE + CALL set_mo_occupation(mo_set=mos(ispin), & + smear=scf_control%smear, & + gce=scf_control%gce) + END IF ! does not yet handle k-points CALL calculate_density_matrix(mos(ispin), & @@ -2070,8 +2094,14 @@ CONTAINS END IF END DO !ispin - CALL set_mo_occupation(mo_array=mos, & - smear=scf_control%smear) + IF (.NOT. scf_control%gce%do_gce) THEN + CALL set_mo_occupation(mo_array=mos, & + smear=scf_control%smear) + ELSE + CALL set_mo_occupation(mo_array=mos, & + smear=scf_control%smear, & + gce=scf_control%gce) + END IF DO ispin = 1, nspins ! does not yet handle k-points diff --git a/src/qs_scf_output.F b/src/qs_scf_output.F index 292cd37b00..4cbcb2e85c 100644 --- a/src/qs_scf_output.F +++ b/src/qs_scf_output.F @@ -97,7 +97,8 @@ MODULE qs_scf_output qs_scf_write_mos, & qs_scf_cdft_info, & qs_scf_cdft_initial_info, & - qs_scf_cdft_constraint_info + qs_scf_cdft_constraint_info, & + qs_scf_gce_info CONTAINS @@ -876,6 +877,18 @@ CONTAINS e_extrapolated = energy%total - 0.5_dp*energy%kTS WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") & "Total energy (extrapolated to T->0): ", e_extrapolated + IF (scf_control%gce%do_gce) THEN + WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & + "GCE work function [eV]: ", scf_control%gce%prev_workfunction*evolt + WRITE (UNIT=output_unit, FMT="((T3,A,T56,ES25.10))") & + "GCE WF-TWF [eV]: ", (scf_control%gce%prev_workfunction - & + scf_control%gce%target_workfunction)*evolt + WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") & + "GCE charge [e]: ", dft_control%pcc_control%charge + WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") & + "GCE free energy: ", (dft_control%pcc_control%charge + dft_control%charge) & + *scf_control%gce%prev_workfunction*evolt + END IF CASE (smear_gaussian) e_extrapolated = energy%total - 0.5_dp*energy%kTS WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") & @@ -1300,4 +1313,43 @@ CONTAINS END SUBROUTINE qs_scf_cdft_constraint_info +! ************************************************************************************************** +!> \brief Print grand canonical SCF information for the current SCF iteration. +!> \param output_unit output unit used for SCF program run information +!> \param qs_env QS environment +!> \param just_energy whether this is an energy-only step +! ************************************************************************************************** + SUBROUTINE qs_scf_gce_info(output_unit, qs_env, just_energy) + + INTEGER, INTENT(IN) :: output_unit + TYPE(qs_environment_type), POINTER :: qs_env + LOGICAL, INTENT(IN) :: just_energy + + REAL(KIND=dp) :: charge, current_wf_ev, delta_wf_ev, & + free_ener, target_wf_ev + TYPE(dft_control_type), POINTER :: dft_control + + IF (output_unit <= 0) RETURN + IF (just_energy) RETURN + + current_wf_ev = qs_env%scf_control%gce%prev_workfunction*evolt + target_wf_ev = qs_env%scf_control%gce%target_workfunction*evolt + delta_wf_ev = current_wf_ev - target_wf_ev + + CALL get_qs_env(qs_env, dft_control=dft_control) + charge = dft_control%pcc_control%charge + free_ener = (charge + dft_control%charge)*qs_env%scf_control%gce%prev_workfunction + + WRITE (UNIT=output_unit, & + FMT="(T8,A,T13,A,T24,A,T27,F6.1,A,T40,A,T56,A,T59,ES10.2,A)") & + "GCE", "WF", "=", current_wf_ev, " eV", & + "WF-TWF", "=", delta_wf_ev, " eV" + + WRITE (UNIT=output_unit, & + FMT="(T13,A,T24,A,T27,F7.3,A,T40,A,T56,A,T59,F14.10,A)") & + "Charge", "=", charge, " e", & + "GCE free energy", "=", free_ener, " a.u." + + END SUBROUTINE qs_scf_gce_info + END MODULE qs_scf_output diff --git a/src/scf_control_types.F b/src/scf_control_types.F index 8f4ccb2c01..d0adbf13a5 100644 --- a/src/scf_control_types.F +++ b/src/scf_control_types.F @@ -152,7 +152,8 @@ MODULE scf_control_types REAL(KIND=dp) :: target_workfunction = 0.16_dp, & ref_esp = 0.0_dp, & mixing_coef = 0.3_dp, & - prev_workfunction = -1001.0_dp + prev_workfunction = -1001.0_dp, & + nelec = 0.0_dp END TYPE gce_type CONTAINS diff --git a/tests/QS/regtest-gcscf/Li_gcscf.inp b/tests/QS/regtest-gcscf/Li_gcscf.inp new file mode 100644 index 0000000000..40073e7efe --- /dev/null +++ b/tests/QS/regtest-gcscf/Li_gcscf.inp @@ -0,0 +1,77 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT Li_gcscf + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT + LSD + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 100 + REL_CUTOFF 40 + &END MGRID + &PLANAR_COUNTER_CHARGE ON + DIST_EDGE 1.5 + GAU_C 1.5 + PARALLEL_PLANE XY + &END PLANAR_COUNTER_CHARGE + &QS + EPS_DEFAULT 1.0E-10 + &END QS + &SCF + ADDED_MOS 10 + CHOLESKY INVERSE + EPS_SCF 1.0E-9 + MAX_SCF 100 + &DIAGONALIZATION + ALGORITHM STANDARD + &END DIAGONALIZATION + &GCE ON + MIXING_COEF 0.3 + TARGET_WORKFUNCTION [eV] 4.5 + &END GCE + &MIXING + ALPHA 0.1 + METHOD NEW_PULAY_MIXING + NBROYDEN 5 + QK 1 + QKAPPA 0.25 + QM 0.75 + &END MIXING + &SMEAR ON + ELECTRONIC_TEMPERATURE [K] 300 + METHOD FERMI_DIRAC + &END SMEAR + &END SCF + &XC + &XC_FUNCTIONAL PADE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 9.0 6.0 10.0 + PERIODIC XYZ + &END CELL + &COORD + Li 1.5 1.5 0.9 + Li 4.5 1.5 0.9 + Li 7.5 1.5 0.9 + Li 1.5 4.5 0.9 + Li 4.5 4.5 0.9 + Li 7.5 4.5 0.9 + &END COORD + &KIND Li + BASIS_SET DZVP-MOLOPT-SR-GTH-q3 + POTENTIAL GTH-LDA-q3 + &END KIND + &TOPOLOGY + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-gcscf/TEST_FILES.toml b/tests/QS/regtest-gcscf/TEST_FILES.toml new file mode 100644 index 0000000000..6cb119fb04 --- /dev/null +++ b/tests/QS/regtest-gcscf/TEST_FILES.toml @@ -0,0 +1,2 @@ +# runs are executed in the same order as in this file +"Li_gcscf.inp" = [{matcher="E_total", tol=1e-9, ref=-44.036008353017763}] diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index a6ce63bb11..464a34eb85 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -12,6 +12,7 @@ QS/regtest-dft-d3-lib-ref libint s_dftd3 QS/regtest-dft-d3bj-lib-ref libint s_dftd3 QS/regtest-pcc QS/regtest-new-pulay-mixing +QS/regtest-gcscf QS/regtest-kubo-transport QS/regtest-gw-realspace libint !ifx QS/regtest-moments-kpoints