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

Add grand canonical SCF for slab models (#5467)

Browse source 
Co-authored-by: Ole Schütt <ole@schuett.name>
This commit is contained in:
Ziwei Chai 2026-07-01 18:23:41 +08:00 • committed by GitHub
parent e6548e2930
commit 23ff2abb53
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
12 changed files with 369 additions and 72 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/cp_control_types.F
View file

@ -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
! **************************************************************************************************

2
src/gce_methods.F
View file

@ -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)

36
src/input_cp2k_scf.F
View file

@ -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.)

37
src/qs_environment.F
View file

@ -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)

150
src/qs_mo_occupation.F
View file

@ -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

26
src/qs_scf.F
View file

@ -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

50
src/qs_scf_diagonalization.F
View file

@ -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

54
src/qs_scf_output.F
View file

@ -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

3
src/scf_control_types.F
View file

@ -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

77
tests/QS/regtest-gcscf/Li_gcscf.inp Normal file
View file

@ -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

2
tests/QS/regtest-gcscf/TEST_FILES.toml Normal file
View file

@ -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}]

1
tests/TEST_DIRS
View file

@ -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