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

PP optimization: Reference can be also a PP (#4635)

Browse source
This commit is contained in:
Juerg Hutter 2025-12-27 16:51:21 +01:00 • committed by GitHub
parent e5ce795793
commit 1971565b89
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
6 changed files with 207 additions and 22 deletions
Download patch file Download diff file Expand all files Collapse all files

10
src/atom_fit.F
View file

@ -630,7 +630,7 @@ CONTAINS
IF (atom%state%multiplicity == -1) THEN
! no spin polarization
atom%weight = wem(i, j)
ncore = get_maxn_occ(atom%state%core)
ncore = get_maxn_occ(atom%state%core) - get_maxn_occ(atom_refs(i, j)%atom%state%core)
DO l = 0, lmat
np = atom%state%maxn_calc(l)
DO k = 1, np
@ -805,7 +805,7 @@ CONTAINS
CALL pseudo_fit(atom_info, wfn_guess, ppot, ostate%f, wtot, pval, dener, w_ener, t_ener, .TRUE.)
IF (nre == 1) THEN
WRITE (iunit, '(/," POWELL| Initial errors of target values")')
afun = ostate%f*wtot
afun = ostate%f*wtot + 1.0E-8_dp
DO i = 1, SIZE(atom_info, 1)
DO j = 1, SIZE(atom_info, 2)
atom => atom_info(i, j)%atom
@ -954,7 +954,7 @@ CONTAINS
WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp), fopt
CALL put_pseudo_param(ostate%xopt, ppot%gth_pot, noopt_nlcc)
CALL atom_write_pseudo_param(ppot%gth_pot)
CALL atom_write_pseudo_param(ppot%gth_pot, fopt=fopt)
END IF
IF (fopt < 1.e-12_dp) EXIT
@ -972,7 +972,7 @@ CONTAINS
ostate%state = 8
CALL powell_optimize(ostate%nvar, x, ostate)
CALL put_pseudo_param(x, ppot%gth_pot, noopt_nlcc)
CALL atom_write_pseudo_param(ppot%gth_pot)
CALL atom_write_pseudo_param(ppot%gth_pot, fopt=fopt)
! dx < SQRT(ostate%rhoend)
IF ((dx*dx) < ostate%rhoend) EXIT
ostate%rhobeg = step_size_scaling*ostate%rhobeg
@ -985,7 +985,7 @@ CONTAINS
CALL put_pseudo_param(x, ppot%gth_pot, noopt_nlcc)
CALL pseudo_fit(atom_info, wfn_guess, ppot, ostate%f, wtot, pval, dener, w_ener, t_ener, .FALSE.)
afun = wtot*ostate%f
afun = wtot*ostate%f + 1.0E-8_dp
WRITE (iunit, '(/," POWELL| Final errors of target values")')
DO i = 1, SIZE(atom_info, 1)

10
src/atom_output.F
View file

@ -757,14 +757,16 @@ CONTAINS
!> \brief Print GTH pseudo-potential parameters.
!> \param gthpot pseudo-potential
!> \param iunit output file unit
!> \param fopt ...
!> \par History
!> * 09.2012 created [Juerg Hutter]
!> \note The pseudo-potential is written into the 'iunit' file unit or as the file 'GTH-PARAMETER'
!> inside the current working directory if the I/O unit is not given explicitly.
! **************************************************************************************************
SUBROUTINE atom_write_pseudo_param(gthpot, iunit)
SUBROUTINE atom_write_pseudo_param(gthpot, iunit, fopt)
TYPE(atom_gthpot_type), INTENT(INOUT) :: gthpot
INTEGER, INTENT(IN), OPTIONAL :: iunit
REAL(KIND=dp), INTENT(IN), OPTIONAL :: fopt
INTEGER :: i, iw, j, k, n
@ -773,7 +775,11 @@ CONTAINS
ELSE
CALL open_file(file_name="GTH-PARAMETER", file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
END IF
WRITE (iw, '(A)') TRIM(ADJUSTL(gthpot%symbol))//" "//TRIM(ADJUSTL(gthpot%pname))
IF (PRESENT(fopt)) THEN
WRITE (iw, '(A,F30.8)') "# "//TRIM(ADJUSTL(gthpot%symbol)), fopt
ELSE
WRITE (iw, '(A)') TRIM(ADJUSTL(gthpot%symbol))//" "//TRIM(ADJUSTL(gthpot%pname))
END IF
WRITE (iw, '(4I5)') gthpot%econf(0:3)
WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%rc, gthpot%ncl, (gthpot%cl(i), i=1, gthpot%ncl)
IF (gthpot%lpotextended) THEN

98
src/atom_pseudo.F
View file

@ -44,6 +44,7 @@ MODULE atom_pseudo
USE periodic_table, ONLY: nelem,&
ptable
USE physcon, ONLY: bohr
USE string_utilities, ONLY: uppercase
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -67,7 +68,7 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'atom_pseudo_opt'
CHARACTER(LEN=2) :: elem
CHARACTER(LEN=2) :: elem, refm
CHARACTER(LEN=default_string_length), &
DIMENSION(:), POINTER :: tmpstringlist
INTEGER :: ads, do_eric, do_erie, handle, i, im, &
@ -76,7 +77,8 @@ CONTAINS
zval, zz
INTEGER, DIMENSION(0:lmat) :: maxn
INTEGER, DIMENSION(:), POINTER :: cn
LOGICAL :: do_gh, eri_c, eri_e, graph, pp_calc
LOGICAL :: do_gh, do_ppref, eri_c, eri_e, explicit, &
graph, pp_calc
REAL(KIND=dp) :: ne, nm
REAL(KIND=dp), DIMENSION(0:lmat, 10) :: pocc
TYPE(atom_basis_type), POINTER :: ae_basis, pp_basis
@ -85,7 +87,7 @@ CONTAINS
TYPE(atom_orbitals), POINTER :: orbitals
TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info, atom_refs
TYPE(atom_potential_type), POINTER :: ae_pot, p_pot
TYPE(atom_state), POINTER :: state, statepp
TYPE(atom_state), POINTER :: state, statepp, stateref
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: basis_section, method_section, &
opt_section, potential_section, &
@ -121,17 +123,29 @@ CONTAINS
CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%PROGRAM_BANNER")
iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%BASIS_SET", extension=".log")
IF (iw > 0) THEN
CALL atom_print_basis(ae_basis, iw, " All Electron Basis")
CALL atom_print_basis(ae_basis, iw, " All Electron/Reference Basis")
CALL atom_print_basis(pp_basis, iw, " Pseudopotential Basis")
END IF
CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%BASIS_SET")
! check on the reference method
CALL section_vals_val_get(atom_section, "REFERENCE%METHOD", c_val=refm)
CALL uppercase(refm)
do_ppref = (refm == "PP")
! read and setup information on the pseudopotential
NULLIFY (potential_section)
potential_section => section_vals_get_subs_vals(atom_section, "POTENTIAL")
ALLOCATE (ae_pot, p_pot)
CALL init_atom_potential(p_pot, potential_section, zval)
CALL init_atom_potential(ae_pot, potential_section, -1)
IF (do_ppref) THEN
potential_section => section_vals_get_subs_vals(atom_section, "REFERENCE%POTENTIAL")
CALL section_vals_get(potential_section, explicit=explicit)
CPASSERT(explicit)
CALL init_atom_potential(ae_pot, potential_section, zval)
ELSE
CALL init_atom_potential(ae_pot, potential_section, -1)
END IF
IF (.NOT. p_pot%confinement .AND. .NOT. ae_pot%confinement) THEN
!set default confinement potential
p_pot%confinement = .TRUE.
@ -208,7 +222,7 @@ CONTAINS
atom_info(in, im)%atom%xc_section => xc_section
atom_refs(in, im)%atom%xc_section => xc_section
ALLOCATE (state, statepp)
ALLOCATE (state, statepp, stateref)
! get the electronic configuration
CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", i_rep_val=in, &
@ -217,7 +231,7 @@ CONTAINS
pp_calc = INDEX(tmpstringlist(1), "CORE") /= 0
CPASSERT(.NOT. pp_calc)
! set occupations
! set occupations for the full electronic configuration
CALL atom_set_occupation(tmpstringlist, state%occ, state%occupation, state%multiplicity)
state%maxl_occ = get_maxl_occ(state%occ)
state%maxn_occ = get_maxn_occ(state%occ)
@ -232,7 +246,6 @@ CONTAINS
state%maxn_calc(k) = MIN(state%maxn_calc(k), ae_basis%nbas(k))
END DO
state%core = 0._dp
CALL set_atom(atom_refs(in, im)%atom, zcore=zval, pp_calc=.FALSE.)
IF (state%multiplicity /= -1) THEN
! set alpha and beta occupations
@ -258,6 +271,52 @@ CONTAINS
END DO
END IF
! set occupations for reference calculation
IF (do_ppref) THEN
CALL section_vals_val_get(atom_section, "REFERENCE%CORE", explicit=explicit)
IF (explicit) THEN
CALL section_vals_val_get(atom_section, "REFERENCE%CORE", c_vals=tmpstringlist)
ELSE
CALL section_vals_val_get(atom_section, "CORE", c_vals=tmpstringlist)
END IF
CALL atom_set_occupation(tmpstringlist, stateref%core, pocc)
zcore = zval - NINT(SUM(stateref%core))
CALL set_atom(atom_refs(in, im)%atom, zcore=zcore, pp_calc=.TRUE.)
stateref%occ = state%occ - stateref%core
stateref%occupation = 0._dp
DO l = 0, lmat
k = 0
DO i = 1, 10
IF (stateref%occ(l, i) /= 0._dp) THEN
k = k + 1
stateref%occupation(l, k) = state%occ(l, i)
IF (state%multiplicity /= -1) THEN
stateref%occa(l, k) = state%occa(l, i) - stateref%core(l, i)/2
stateref%occb(l, k) = state%occb(l, i) - stateref%core(l, i)/2
END IF
END IF
END DO
END DO
stateref%maxl_occ = get_maxl_occ(stateref%occ)
stateref%maxn_occ = get_maxn_occ(stateref%occ)
stateref%maxl_calc = state%maxl_calc
stateref%maxn_calc = 0
maxn = get_maxn_occ(stateref%core)
DO k = 0, stateref%maxl_calc
stateref%maxn_calc(k) = state%maxn_calc(k) - maxn(k)
stateref%maxn_calc(k) = MIN(stateref%maxn_calc(k), ae_basis%nbas(k))
END DO
stateref%multiplicity = state%multiplicity
ELSE
stateref = state
CALL set_atom(atom_refs(in, im)%atom, zcore=zval, pp_calc=.FALSE.)
END IF
! set occupations for pseudopotential calculation
CALL section_vals_val_get(atom_section, "CORE", c_vals=tmpstringlist)
CALL atom_set_occupation(tmpstringlist, statepp%core, pocc)
@ -294,7 +353,11 @@ CONTAINS
CALL section_vals_val_get(method_section, "METHOD_TYPE", i_val=method, i_rep_section=im)
CALL section_vals_val_get(method_section, "RELATIVISTIC", i_val=reltyp, i_rep_section=im)
CALL set_atom(atom_info(in, im)%atom, method_type=method)
CALL set_atom(atom_refs(in, im)%atom, method_type=method, relativistic=reltyp)
IF (do_ppref) THEN
CALL set_atom(atom_refs(in, im)%atom, method_type=method)
ELSE
CALL set_atom(atom_refs(in, im)%atom, method_type=method, relativistic=reltyp)
END IF
! calculate integrals: pseudopotential basis
! general integrals
@ -314,12 +377,16 @@ CONTAINS
eri_coulomb=eri_c, eri_exchange=eri_e)
! potential
CALL atom_ppint_setup(ae_int, ae_basis, potential=ae_pot)
! relativistic correction terms
CALL atom_relint_setup(ae_int, ae_basis, reltyp, zcore=REAL(zval, dp))
IF (do_ppref) THEN
NULLIFY (ae_int%tzora, ae_int%hdkh)
ELSE
! relativistic correction terms
CALL atom_relint_setup(ae_int, ae_basis, reltyp, zcore=REAL(zval, dp))
END IF
!
CALL set_atom(atom_refs(in, im)%atom, basis=ae_basis, integrals=ae_int, potential=ae_pot)
state%maxn_calc(:) = MIN(state%maxn_calc(:), ae_basis%nbas(:))
CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
stateref%maxn_calc(:) = MIN(stateref%maxn_calc(:), ae_basis%nbas(:))
CPASSERT(ALL(stateref%maxn_calc(:) >= stateref%maxn_occ))
CALL set_atom(atom_info(in, im)%atom, coulomb_integral_type=do_eric, &
exchange_integral_type=do_erie)
@ -337,9 +404,9 @@ CONTAINS
CALL create_atom_orbs(orbitals, mb, mo)
CALL set_atom(atom_info(in, im)%atom, orbitals=orbitals)
CALL set_atom(atom_refs(in, im)%atom, state=state)
CALL set_atom(atom_refs(in, im)%atom, state=stateref)
NULLIFY (orbitals)
mo = MAXVAL(state%maxn_calc)
mo = MAXVAL(stateref%maxn_calc)
mb = MAXVAL(atom_refs(in, im)%atom%basis%nbas)
CALL create_atom_orbs(orbitals, mb, mo)
CALL set_atom(atom_refs(in, im)%atom, orbitals=orbitals)
@ -354,6 +421,7 @@ CONTAINS
CALL calculate_atom(atom_refs(in, im)%atom, iw)
CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SCF_INFO")
END IF
DEALLOCATE (state)
END DO
END DO

43
src/input_cp2k_atom.F
View file

@ -168,6 +168,10 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL create_atom_reference_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_atom_print_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
@ -198,6 +202,45 @@ CONTAINS
END SUBROUTINE create_atom_section
! **************************************************************************************************
!> \brief Create the atom section for reference method
!> \param section the section to create
!> \author jgh
! **************************************************************************************************
SUBROUTINE create_atom_reference_section(section)
TYPE(section_type), POINTER :: section
TYPE(keyword_type), POINTER :: keyword
TYPE(section_type), POINTER :: subsection
CPASSERT(.NOT. ASSOCIATED(section))
CALL section_create(section, __LOCATION__, name="REFERENCE", &
description="Section to specify a reference method for PP optimization.", &
n_keywords=1, n_subsections=1, repeats=.FALSE.)
NULLIFY (keyword)
CALL keyword_create(keyword, __LOCATION__, name="METHOD", &
description="METHOD to be used as reference: [AE|PP] ", &
usage="METHOD PP", &
type_of_var=char_t, default_c_val="AE", n_var=-1)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="CORE", &
description="Specifies the core electrons for a pseudopotential", &
usage="CORE 1s2 ... or CORE [Ne] or CORE none for 0 electron cores", repeats=.FALSE., &
n_var=-1, type_of_var=char_t)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
NULLIFY (subsection)
CALL create_potential_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
END SUBROUTINE create_atom_reference_section
! **************************************************************************************************
!> \brief Create the print atom section
!> \param section the section to create

67
tests/ATOM/regtest-pseudo/H-ECPGTH.inp Normal file
View file

@ -0,0 +1,67 @@
&GLOBAL
PROGRAM_NAME ATOM
&END GLOBAL
&ATOM
CORE 1s0
ELECTRON_CONFIGURATION 1s1
ELEMENT H
MAX_ANGULAR_MOMENTUM 1
RUN_TYPE PSEUDOPOTENTIAL_OPTIMIZATION
&AE_BASIS
&END AE_BASIS
&METHOD
METHOD_TYPE KOHN-SHAM
&XC
&XC_FUNCTIONAL PBE0
&END XC_FUNCTIONAL
&END XC
&END METHOD
&OPTIMIZATION
EPS_SCF 1.e-8
&END OPTIMIZATION
&POTENTIAL
CONFINEMENT 200. 4.0 12.0
CONFINEMENT_TYPE BARRIER
PSEUDO_TYPE GTH
&GTH_POTENTIAL
# H 0.00000000
1 0 0 0
0.19557818168233 2 -4.13775404933047 0.68536512348064
0
&END GTH_POTENTIAL
&END POTENTIAL
&POWELL
ACCURACY 1.e-14
MAX_FUN 50
MAX_INIT 1
STEP_SIZE 0.002
TARGET_POT_SEMICORE [eV] 0.0010
TARGET_POT_VALENCE [eV] 0.0001
TARGET_POT_VIRTUAL [eV] 0.0010
WEIGHT_POT_NODE 10.0
WEIGHT_POT_SEMICORE 2.0
WEIGHT_POT_VALENCE 5.0
WEIGHT_POT_VIRTUAL 1.0
WEIGHT_PSIR0 0.0
&END POWELL
&PP_BASIS
&END PP_BASIS
&REFERENCE
METHOD PP
&POTENTIAL
CONFINEMENT 200. 4.0 12.0
CONFINEMENT_TYPE BARRIER
PSEUDO_TYPE ECP
&ECP
H nelec 0
H ul
1 21.24359508259891 1.00000000000000
3 21.24359508259891 21.24359508259891
2 21.77696655044365 -10.85192405303825
H S
2 1.000000000000000 0.00000000000000
&END ECP
&END POTENTIAL
&END REFERENCE
&END ATOM

1
tests/ATOM/regtest-pseudo/TEST_FILES.toml
View file

@ -23,4 +23,5 @@
"C_grb2.inp" = [{matcher="M035", tol=1e-12, ref=-5.360945665362}]
"Li_NLCC.inp" = [{matcher="M061", tol=1e-12, ref=0.5527955347}]
"Au.inp" = [{matcher="M035", tol=1e-10, ref=-136.583322713061}]
"H-ECPGTH.inp" = [{matcher="M056", tol=1e-08, ref=0.0000000000}]
#EOF