Skip to content

Commit

Permalink
Enable sirius support of dft+U
Browse files Browse the repository at this point in the history
  • Loading branch information
@mtaillefumier @mkrack
mtaillefumier authored and mkrack committed Sep 19, 2022
1 parent e711b3a commit 86aeba6
Show file tree
Hide file tree
Showing 3 changed files with 232 additions and 10 deletions.
85 changes: 85 additions & 0 deletions src/input_cp2k_subsys.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2270,6 +2270,19 @@ SUBROUTINE create_dft_plus_u_section(section)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

#if defined(__SIRIUS)
CALL keyword_create(keyword, __LOCATION__, &
name="N", &
description="principal quantum number of the"// &
"orbitals to which the correction is applied. Ignored unless pwdft is used for the calculations", &
repeats=.FALSE., &
n_var=1, &
type_of_var=integer_t, &
default_i_val=-1, &
usage="N 2")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, &
name="U_MINUS_J", &
variants=(/"U_EFF"/), &
Expand All @@ -2283,6 +2296,78 @@ SUBROUTINE create_dft_plus_u_section(section)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, &
name="U", &
description="U parameter in the theory of Dudarev et al. Ignored unless pwdft is used", &
repeats=.FALSE., &
n_var=1, &
type_of_var=real_t, &
default_r_val=0.0_dp, &
unit_str="au_e", &
usage="U [eV] 1.4")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, &
name="J", &
description="J parameter in the theory of Dudarev et al. Ignored unless pwdft is used", &
repeats=.FALSE., &
n_var=1, &
type_of_var=real_t, &
default_r_val=0.0_dp, &
unit_str="au_e", &
usage="J [eV] 1.4")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, &
name="alpha", &
description="alpha parameter in the theory of Dudarev et al. Ignored unless pwdft is used", &
repeats=.FALSE., &
n_var=1, &
type_of_var=real_t, &
default_r_val=0.0_dp, &
unit_str="au_e", &
usage="alpha [eV] 1.4")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, &
name="beta", &
description="beta parameter in the theory of Dudarev et al. Ignored unless pwdft is used", &
repeats=.FALSE., &
n_var=1, &
type_of_var=real_t, &
default_r_val=0.0_dp, &
unit_str="au_e", &
usage="beta [eV] 1.4")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, &
name="J0", &
description="J0 parameter in the theory of Dudarev et al. Ignored unless pwdft is used", &
repeats=.FALSE., &
n_var=1, &
type_of_var=real_t, &
default_r_val=0.0_dp, &
unit_str="au_e", &
usage="J0 [eV] 1.4")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, &
name="occupation", &
description="number of electrons in the hubbard shell. Ignored unless pwdft is used", &
repeats=.FALSE., &
n_var=1, &
type_of_var=real_t, &
default_r_val=0.0_dp, &
usage="occupation 6")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
#endif

CALL keyword_create(keyword, __LOCATION__, &
name="U_RAMPING", &
description="Increase the effective U parameter stepwise using the specified "// &
Expand Down
110 changes: 105 additions & 5 deletions src/qs_kind_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,12 +123,19 @@ MODULE qs_kind_types
! **************************************************************************************************
TYPE dft_plus_u_type
INTEGER :: l = -1
INTEGER :: n = -1
INTEGER :: max_scf = -1
REAL(KIND=dp) :: eps_u_ramping = 0.0_dp
REAL(KIND=dp) :: eps_scf = HUGE(0.0_dp)
REAL(KIND=dp) :: u_minus_j_target
REAL(KIND=dp) :: u_minus_j = 0.0_dp
REAL(KIND=dp) :: u_ramping = 0.0_dp
REAL(KIND=dp) :: U = 0.0_dp
REAL(KIND=dp) :: J = 0.0_dp
REAL(KIND=dp) :: alpha = 0.0_dp
REAL(KIND=dp) :: beta = 0.0_dp
REAL(KIND=dp) :: J0 = 0.0_dp
REAL(KIND=dp) :: occupation = -1.0_dp
INTEGER, DIMENSION(:), POINTER :: orbitals => Null()
LOGICAL :: init_u_ramping_each_scf = .FALSE.
LOGICAL :: smear = .FALSE.
Expand Down Expand Up @@ -403,7 +410,14 @@ END SUBROUTINE deallocate_qs_kind_set
!> \param lmax_rho0 ...
!> \param dft_plus_u_atom ...
!> \param l_of_dft_plus_u ...
!> \param n_of_dft_plus_u ...
!> \param u_minus_j ...
!> \param U_of_dft_plus_u ...
!> \param J_of_dft_plus_u ...
!> \param alpha_of_dft_plus_u ...
!> \param beta_of_dft_plus_u ...
!> \param J0_of_dft_plus_u ...
!> \param occupation_of_dft_plus_u ...
!> \param dispersion ...
!> \param bs_occupation ...
!> \param magnetization ...
Expand Down Expand Up @@ -440,7 +454,9 @@ SUBROUTINE get_qs_kind(qs_kind, &
covalent_radius, vdw_radius, &
gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, &
ngrid_ang, ngrid_rad, lmax_rho0, &
dft_plus_u_atom, l_of_dft_plus_u, u_minus_j, dispersion, &
dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, &
u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, &
alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, &
bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, &
max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, &
init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, &
Expand Down Expand Up @@ -476,8 +492,9 @@ SUBROUTINE get_qs_kind(qs_kind, &
TYPE(grid_atom_type), OPTIONAL, POINTER :: grid_atom
INTEGER, INTENT(OUT), OPTIONAL :: ngrid_ang, ngrid_rad, lmax_rho0
LOGICAL, INTENT(OUT), OPTIONAL :: dft_plus_u_atom
INTEGER, INTENT(OUT), OPTIONAL :: l_of_dft_plus_u
REAL(KIND=dp), INTENT(OUT), OPTIONAL :: u_minus_j
INTEGER, INTENT(OUT), OPTIONAL :: l_of_dft_plus_u, n_of_dft_plus_u
REAL(KIND=dp), INTENT(OUT), OPTIONAL :: u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, &
alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u
TYPE(qs_atom_dispersion_type), OPTIONAL, POINTER :: dispersion
LOGICAL, INTENT(OUT), OPTIONAL :: bs_occupation
REAL(KIND=dp), INTENT(OUT), OPTIONAL :: magnetization
Expand Down Expand Up @@ -683,6 +700,13 @@ SUBROUTINE get_qs_kind(qs_kind, &
l_of_dft_plus_u = -1
END IF
END IF
IF (PRESENT(n_of_dft_plus_u)) THEN
IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
n_of_dft_plus_u = qs_kind%dft_plus_u%n
ELSE
n_of_dft_plus_u = -1
END IF
END IF
IF (PRESENT(u_minus_j)) THEN
IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
u_minus_j = qs_kind%dft_plus_u%u_minus_j
Expand All @@ -697,6 +721,49 @@ SUBROUTINE get_qs_kind(qs_kind, &
u_minus_j_target = 0.0_dp
END IF
END IF
IF (PRESENT(U_of_dft_plus_u)) THEN
IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
U_of_dft_plus_u = qs_kind%dft_plus_u%U
ELSE
U_of_dft_plus_u = 0.0_dp
END IF
END IF
IF (PRESENT(J_of_dft_plus_u)) THEN
IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
J_of_dft_plus_u = qs_kind%dft_plus_u%J
ELSE
J_of_dft_plus_u = 0.0_dp
END IF
END IF
IF (PRESENT(alpha_of_dft_plus_u)) THEN
IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
alpha_of_dft_plus_u = qs_kind%dft_plus_u%alpha
ELSE
alpha_of_dft_plus_u = 0.0_dp
END IF
END IF
IF (PRESENT(beta_of_dft_plus_u)) THEN
IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
beta_of_dft_plus_u = qs_kind%dft_plus_u%beta
ELSE
beta_of_dft_plus_u = 0.0_dp
END IF
END IF
IF (PRESENT(J0_of_dft_plus_u)) THEN
IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
J0_of_dft_plus_u = qs_kind%dft_plus_u%J0
ELSE
J0_of_dft_plus_u = 0.0_dp
END IF
END IF
IF (PRESENT(occupation_of_dft_plus_u)) THEN
IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
occupation_of_dft_plus_u = qs_kind%dft_plus_u%occupation
ELSE
occupation_of_dft_plus_u = -1.0_dp
END IF
END IF

IF (PRESENT(init_u_ramping_each_scf)) THEN
IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
init_u_ramping_each_scf = qs_kind%dft_plus_u%init_u_ramping_each_scf
Expand Down Expand Up @@ -1441,7 +1508,7 @@ SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_f
DIMENSION(maxbas) :: basis_set_form, basis_set_name, &
basis_set_type
INTEGER :: handle, i, i_rep, iounit, ipaodesc, ipaopot, ipos, j, jj, k_rep, l, m, n_rep, &
nb_rep, nexp, ngauss, nlcc, nloc, nnl, norbitals, npaodesc, npaopot, nppnl, nspin, z
nb_rep, nexp, ngauss, nlcc, nloc, nnl, norbitals, npaodesc, npaopot, nppnl, nspin, nu, z
INTEGER, DIMENSION(:), POINTER :: add_el, elec_conf, orbitals
LOGICAL :: check, explicit, explicit_basis, explicit_kgpot, explicit_potential, nobasis, &
section_enabled, subsection_enabled, update_input
Expand Down Expand Up @@ -1786,7 +1853,6 @@ SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_f
spin_section => section_vals_get_subs_vals(bs_section, "BETA")
CALL section_vals_get(spin_section, explicit=explicit)
IF (explicit) THEN
NULLIFY (add_el)
CALL section_vals_val_get(spin_section, &
keyword_name="NEL", i_vals=add_el)
Expand Down Expand Up @@ -1832,6 +1898,40 @@ SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_f
qs_kind%dft_plus_u%l = l
IF (l < 0) &
CPABORT("DFT+U| Invalid orbital angular momentum quantum number specified: l < 0")
#if defined(__SIRIUS)
CALL section_vals_val_get(dft_plus_u_section, &
keyword_name="N", &
i_val=nu)
qs_kind%dft_plus_u%n = nu
IF (nu < 0) &
CPABORT("DFT+U| the principal quantum number n can not be negative")
CALL section_vals_val_get(dft_plus_u_section, &
keyword_name="U", &
r_val=qs_kind%dft_plus_u%U)
CALL section_vals_val_get(dft_plus_u_section, &
keyword_name="J", &
r_val=qs_kind%dft_plus_u%J)
CALL section_vals_val_get(dft_plus_u_section, &
keyword_name="alpha", &
r_val=qs_kind%dft_plus_u%alpha)
CALL section_vals_val_get(dft_plus_u_section, &
keyword_name="beta", &
r_val=qs_kind%dft_plus_u%beta)
CALL section_vals_val_get(dft_plus_u_section, &
keyword_name="J0", &
r_val=qs_kind%dft_plus_u%J0)
CALL section_vals_val_get(dft_plus_u_section, &
keyword_name="occupation", &
r_val=qs_kind%dft_plus_u%occupation)
#endif
CALL section_vals_val_get(dft_plus_u_section, &
keyword_name="U_MINUS_J", &
r_val=qs_kind%dft_plus_u%u_minus_j_target)
Expand Down
47 changes: 42 additions & 5 deletions src/sirius_interface.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,17 +130,18 @@ SUBROUTINE cp_sirius_create_env(pwdft_env)
CHARACTER(len=default_string_length) :: label
INTEGER :: i, iatom, ibeta, ifun, ikind, iwf, j, l, &
n, natom, nbeta, nkind, nmesh, &
num_mag_dims, sirius_mpi_comm, vdw_func
num_mag_dims, sirius_mpi_comm, vdw_func, nu, lu
INTEGER, DIMENSION(:), POINTER :: mpi_grid_dims
INTEGER(KIND=C_INT), DIMENSION(3) :: k_grid, k_shift
INTEGER, DIMENSION(:), POINTER :: kk
LOGICAL :: up, use_ref_cell
LOGICAL(4) :: use_so, use_symmetry
LOGICAL(4) :: use_so, use_symmetry, dft_plus_u_atom
REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:) :: fun
REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:, :) :: dion
REAL(KIND=C_DOUBLE), DIMENSION(3) :: a1, a2, a3, v1, v2
REAL(KIND=dp) :: al, angle1, angle2, cval, focc, &
magnetization, mass, pf, rl, zeff
magnetization, mass, pf, rl, zeff, alpha_u, beta_u, &
J0_u, J_u, U_u, occ_u, u_minus_J
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, ef, fe, locpot, rc, rp
REAL(KIND=dp), DIMENSION(3) :: vr, vs
REAL(KIND=dp), DIMENSION(:), POINTER :: density
Expand Down Expand Up @@ -285,10 +286,12 @@ SUBROUTINE cp_sirius_create_env(pwdft_env)
CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
NULLIFY (upf_pot, gth_potential)
CALL get_qs_kind(qs_kind_set(ikind), upf_potential=upf_pot, gth_potential=gth_potential)

IF (ASSOCIATED(upf_pot)) THEN
CALL sirius_add_atom_type(sctx, label, fname=upf_pot%filename, &
symbol=element_symbol, &
mass=REAL(mass/massunit, KIND=C_DOUBLE))

ELSEIF (ASSOCIATED(gth_potential)) THEN
!
NULLIFY (atom_grid)
Expand Down Expand Up @@ -393,6 +396,7 @@ SUBROUTINE cp_sirius_create_env(pwdft_env)
DO iwf = 1, SIZE(wavefunction, 2)
focc = wfninfo(1, iwf)
l = NINT(wfninfo(2, iwf))
nu = NINT(wfninfo(3, iwf))
IF (up) THEN
fun(1:nmesh) = wavefunction(1:nmesh, iwf)*rp(i)
ELSE
Expand All @@ -402,7 +406,7 @@ SUBROUTINE cp_sirius_create_env(pwdft_env)
END IF
CALL sirius_add_atom_type_radial_function(sctx, &
label, "ps_atomic_wf", &
fun(1), nmesh, l=l, occ=REAL(focc, KIND=C_DOUBLE), n=-1)
fun(1), nmesh, l=l, occ=REAL(focc, KIND=C_DOUBLE), n=nu)
END DO

! set total charge density of a free atom (to compute initial rho(r))
Expand All @@ -416,6 +420,39 @@ SUBROUTINE cp_sirius_create_env(pwdft_env)
CALL sirius_add_atom_type_radial_function(sctx, label, "ps_rho_total", &
fun(1), nmesh)

CALL get_qs_kind(qs_kind_set(ikind), &
dft_plus_u_atom=dft_plus_u_atom, &
l_of_dft_plus_u=lu, &
n_of_dft_plus_u=nu, &
u_minus_j_target=u_minus_j, &
U_of_dft_plus_u=U_u, &
J_of_dft_plus_u=J_u, &
alpha_of_dft_plus_u=alpha_u, &
beta_of_dft_plus_u=beta_u, &
J0_of_dft_plus_u=J0_u, &
occupation_of_dft_plus_u=occ_u)

IF (dft_plus_u_atom) THEN
IF (nu < 1) THEN
CALL cpabort("CP2K/SIRIUS (hubbard): principal quantum number not specified")
END IF

IF (lu < 0) THEN
CALL cpabort("CP2K/SIRIUS (hubbard): l can not be negative.")
END IF

IF (occ_u < 0.0) THEN
CALL cpabort("CP2K/SIRIUS (hubbard): the occupation number can not be negative.")
END IF

IF (ABS(u_minus_j) < 1e-8) THEN
CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
NINT(occ_u), U_u, J_u, alpha_u, beta_u, J0_u)
ELSE
CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
NINT(occ_u), u_minus_j, 0.0_dp, alpha_u, beta_u, J0_u)
END IF
END IF
IF (ASSOCIATED(density)) DEALLOCATE (density)
IF (ASSOCIATED(wavefunction)) DEALLOCATE (wavefunction)
IF (ASSOCIATED(wfninfo)) DEALLOCATE (wfninfo)
Expand All @@ -426,7 +463,7 @@ SUBROUTINE cp_sirius_create_env(pwdft_env)
!
ELSE
CALL cp_abort(__LOCATION__, &
'CP2K/SIRIUS: atomic kind needs UPF or GTH potential definition')
"CP2K/SIRIUS: atomic kind needs UPF or GTH potential definition")
END IF
END DO

Expand Down

0 comments on commit 86aeba6

Please sign in to comment.