Skip to content

Commit

Permalink
Atom PP code: NLCC options added (#507)
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Aug 9, 2019
1 parent 95014e7 commit f4b239e
Show file tree
Hide file tree
Showing 4 changed files with 248 additions and 17 deletions.
183 changes: 166 additions & 17 deletions src/atom_fit.F
View file
Expand Up @@ -19,14 +19,20 @@ MODULE atom_fit
atom_write_pseudo_param
USE atom_types, ONLY: &
GTO_BASIS, STO_BASIS, atom_basis_type, atom_gthpot_type, atom_integrals, atom_p_type, &
atom_potential_type, atom_type, create_opgrid, lmat, opgrid_type, release_opgrid, set_atom
atom_potential_type, atom_type, create_opgrid, create_opmat, lmat, opgrid_type, &
opmat_type, release_opgrid, release_opmat, set_atom
USE atom_utils, ONLY: &
atom_completeness, atom_condnumber, atom_consistent_method, atom_denmat, atom_density, &
atom_orbital_charge, atom_orbital_max, atom_orbital_nodes, atom_wfnr0, get_maxn_occ, &
integrate_grid
atom_completeness, atom_condnumber, atom_consistent_method, atom_core_density, &
atom_denmat, atom_density, atom_orbital_charge, atom_orbital_max, atom_orbital_nodes, &
atom_wfnr0, get_maxn_occ, integrate_grid
USE cp_files, ONLY: close_file,&
open_file
USE input_constants, ONLY: do_analytic
USE input_constants, ONLY: do_analytic,&
do_rhf_atom,&
do_rks_atom,&
do_rohf_atom,&
do_uhf_atom,&
do_uks_atom
USE input_section_types, ONLY: section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
Expand Down Expand Up @@ -495,7 +501,7 @@ SUBROUTINE atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)
n10, np, nre, nreinit, ntarget
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: xtob
INTEGER, DIMENSION(0:lmat) :: ncore
LOGICAL :: explicit
LOGICAL :: explicit, noopt_nlcc, preopt_nlcc
REAL(KIND=dp) :: afun, charge, de, deig, drho, dx, eig, fopt, oc, pchg, peig, pv, rcm, rcov, &
rmax, semicore_level, step_size_scaling, t_ener, t_psir0, t_semi, t_valence, t_virt, &
w_ener, w_node, w_psir0, w_semi, w_valence, w_virt, wtot
Expand Down Expand Up @@ -532,6 +538,9 @@ SUBROUTINE atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)

CALL section_vals_val_get(powell_section, "SEMICORE_LEVEL", r_val=semicore_level)

CALL section_vals_val_get(powell_section, "NOOPT_NLCC", l_val=noopt_nlcc)
CALL section_vals_val_get(powell_section, "PREOPT_NLCC", l_val=preopt_nlcc)

n = SIZE(atom_info, 1)
m = SIZE(atom_info, 2)
ALLOCATE (wem(n, m))
Expand Down Expand Up @@ -572,9 +581,16 @@ SUBROUTINE atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)
CALL open_file(file_name="POWELL_RESULT", file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
END IF

IF (ppot%gth_pot%nlcc) THEN
CALL opt_nlcc_param(atom_info, atom_refs, ppot%gth_pot, iunit, preopt_nlcc)
ELSE
noopt_nlcc = .TRUE.
preopt_nlcc = .FALSE.
END IF

ALLOCATE (xi(200))
!decide here what to optimize
CALL get_pseudo_param(xi, ostate%nvar, ppot%gth_pot)
CALL get_pseudo_param(xi, ostate%nvar, ppot%gth_pot, noopt_nlcc)
ALLOCATE (x(ostate%nvar))
x(1:ostate%nvar) = xi(1:ostate%nvar)
DEALLOCATE (xi)
Expand Down Expand Up @@ -602,7 +618,7 @@ SUBROUTINE atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)
IF (atom%state%multiplicity == -1) THEN
! no spin polarization
atom%weight = wem(i, j)
ncore = get_maxn_occ(atom_info(i, j)%atom%state%core)
ncore = get_maxn_occ(atom%state%core)
DO l = 0, lmat
np = atom%state%maxn_calc(l)
DO k = 1, np
Expand Down Expand Up @@ -773,7 +789,7 @@ SUBROUTINE atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)
ALLOCATE (xi(ostate%nvar))
DO nre = 1, nreinit
xi(:) = x(:)
CALL put_pseudo_param(x, ppot%gth_pot)
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, .TRUE.)
IF (nre == 1) THEN
WRITE (iunit, '(/," POWELL| Initial errors of target values")')
Expand Down Expand Up @@ -910,7 +926,7 @@ SUBROUTINE atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)
DO

IF (ostate%state == 2) THEN
CALL put_pseudo_param(x, ppot%gth_pot)
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.)
fopt = MIN(fopt, ostate%f)
END IF
Expand All @@ -925,7 +941,7 @@ SUBROUTINE atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)
IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0 .AND. ostate%nf > 2) THEN
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)
CALL put_pseudo_param(ostate%xopt, ppot%gth_pot, noopt_nlcc)
CALL atom_write_pseudo_param(ppot%gth_pot)
END IF

Expand All @@ -943,7 +959,7 @@ SUBROUTINE atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)
END IF
ostate%state = 8
CALL powell_optimize(ostate%nvar, x, ostate)
CALL put_pseudo_param(x, ppot%gth_pot)
CALL put_pseudo_param(x, ppot%gth_pot, noopt_nlcc)
CALL atom_write_pseudo_param(ppot%gth_pot)
! dx < SQRT(ostate%rhoend)
IF ((dx*dx) < ostate%rhoend) EXIT
Expand All @@ -955,7 +971,7 @@ SUBROUTINE atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)
WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') ostate%fopt

CALL put_pseudo_param(x, ppot%gth_pot)
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

Expand Down Expand Up @@ -1144,6 +1160,135 @@ SUBROUTINE atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)
END IF

END SUBROUTINE atom_fit_pseudo

! **************************************************************************************************
!> \brief Fit NLCC density on core densities
!> \param atom_info ...
!> \param atom_refs ...
!> \param gthpot ...
!> \param iunit ...
!> \param preopt_nlcc ...
! **************************************************************************************************
SUBROUTINE opt_nlcc_param(atom_info, atom_refs, gthpot, iunit, preopt_nlcc)
TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info, atom_refs
TYPE(atom_gthpot_type), INTENT(inout) :: gthpot
INTEGER, INTENT(in) :: iunit
LOGICAL, INTENT(in) :: preopt_nlcc

CHARACTER(len=*), PARAMETER :: routineN = 'opt_nlcc_param', routineP = moduleN//':'//routineN

INTEGER :: i, im, j, k, method
REAL(KIND=dp) :: rcov, zcore, zrc, zrch
TYPE(atom_type), POINTER :: aref, atom
TYPE(opgrid_type), POINTER :: corden, den, den1, den2, density
TYPE(opmat_type), POINTER :: denmat, dma, dmb

CPASSERT(gthpot%nlcc)

IF (iunit > 0) THEN
WRITE (iunit, '(/," Core density information")')
WRITE (iunit, '(A,T37,A,T55,A,T75,A)') " State Density:", "Full", "Rcov/2", "Rcov/4"
END IF

rcov = ptable(atom_refs(1, 1)%atom%z)%covalent_radius*bohr
atom => atom_info(1, 1)%atom
NULLIFY (density)
CALL create_opgrid(density, atom%basis%grid)
density%op = 0.0_dp
im = 0
DO i = 1, SIZE(atom_info, 1)
DO j = 1, SIZE(atom_info, 2)
atom => atom_info(i, j)%atom
aref => atom_refs(i, j)%atom
IF (atom_consistent_method(atom%method_type, atom%state%multiplicity)) THEN
NULLIFY (den, denmat)
CALL create_opgrid(den, atom%basis%grid)
CALL create_opmat(denmat, atom%basis%nbas)

method = atom%method_type
SELECT CASE (method)
CASE (do_rks_atom, do_rhf_atom)
CALL atom_denmat(denmat%op, aref%orbitals%wfn, &
atom%basis%nbas, atom%state%core, &
aref%state%maxl_occ, aref%state%maxn_occ)
CASE (do_uks_atom, do_uhf_atom)
NULLIFY (dma, dmb)
CALL create_opmat(dma, atom%basis%nbas)
CALL create_opmat(dmb, atom%basis%nbas)
CALL atom_denmat(dma%op, aref%orbitals%wfna, &
atom%basis%nbas, atom%state%core, &
aref%state%maxl_occ, aref%state%maxn_occ)
CALL atom_denmat(dmb%op, aref%orbitals%wfnb, &
atom%basis%nbas, atom%state%core, &
aref%state%maxl_occ, aref%state%maxn_occ)
denmat%op = 0.5_dp*(dma%op+dmb%op)
CALL release_opmat(dma)
CALL release_opmat(dmb)
CASE (do_rohf_atom)
CPABORT("")
CASE DEFAULT
CPABORT("")
END SELECT

im = im+1
CALL atom_density(den%op, denmat%op, atom%basis, aref%state%maxl_occ, typ="RHO")
density%op = density%op+den%op
zcore = integrate_grid(den%op, atom%basis%grid)
zcore = fourpi*zcore
NULLIFY (den1, den2)
CALL create_opgrid(den1, atom%basis%grid)
CALL create_opgrid(den2, atom%basis%grid)
den1%op = 0.0_dp
den2%op = 0.0_dp
DO k = 1, atom%basis%grid%nr
IF (atom%basis%grid%rad(k) < 0.50_dp*rcov) THEN
den1%op(k) = den%op(k)
END IF
IF (atom%basis%grid%rad(k) < 0.25_dp*rcov) THEN
den2%op(k) = den%op(k)
END IF
END DO
zrc = integrate_grid(den1%op, atom%basis%grid)
zrc = fourpi*zrc
zrch = integrate_grid(den2%op, atom%basis%grid)
zrch = fourpi*zrch
CALL release_opgrid(den1)
CALL release_opgrid(den2)
CALL release_opgrid(den)
CALL release_opmat(denmat)
IF (iunit > 0) THEN
WRITE (iunit, '(2I5,T31,F10.4,T51,F10.4,T71,F10.4)') i, j, zcore, zrc, zrch
END IF
END IF
END DO
END DO
density%op = density%op/REAL(im, KIND=dp)
!
IF (preopt_nlcc) THEN
gthpot%nexp_nlcc = 1
gthpot%nct_nlcc = 0
gthpot%cval_nlcc = 0._dp
gthpot%alpha_nlcc = 0._dp
gthpot%nct_nlcc(1) = 1
gthpot%cval_nlcc(1, 1) = 1.0_dp
gthpot%alpha_nlcc(1) = gthpot%rc
NULLIFY (corden)
CALL create_opgrid(corden, atom%basis%grid)
CALL atom_core_density(corden%op, atom%potential, typ="RHO", rr=atom%basis%grid%rad)
DO k = 1, atom%basis%grid%nr
IF (atom%basis%grid%rad(k) > 0.25_dp*rcov) THEN
corden%op(k) = 0.0_dp
END IF
END DO
zrc = integrate_grid(corden%op, atom%basis%grid)
zrc = fourpi*zrc
gthpot%cval_nlcc(1, 1) = zrch/zrc
CALL release_opgrid(corden)
END IF
CALL release_opgrid(density)

END SUBROUTINE opt_nlcc_param

! **************************************************************************************************
!> \brief Low level routine to fit an atomic electron density.
!> \param density electron density on an atomic radial grid 'atom%basis%grid'
Expand Down Expand Up @@ -1552,11 +1697,13 @@ END FUNCTION get_error_value
!> \param pvec ...
!> \param nval ...
!> \param gthpot ...
!> \param noopt_nlcc ...
! **************************************************************************************************
SUBROUTINE get_pseudo_param(pvec, nval, gthpot)
SUBROUTINE get_pseudo_param(pvec, nval, gthpot, noopt_nlcc)
REAL(KIND=dp), DIMENSION(:), INTENT(out) :: pvec
INTEGER, INTENT(out) :: nval
TYPE(atom_gthpot_type), INTENT(in) :: gthpot
LOGICAL, INTENT(in) :: noopt_nlcc

CHARACTER(len=*), PARAMETER :: routineN = 'get_pseudo_param', &
routineP = moduleN//':'//routineN
Expand Down Expand Up @@ -1592,7 +1739,7 @@ SUBROUTINE get_pseudo_param(pvec, nval, gthpot)
END DO
END DO
END IF
IF (gthpot%nlcc) THEN
IF (gthpot%nlcc .AND. (.NOT. noopt_nlcc)) THEN
DO j = 1, gthpot%nexp_nlcc
ival = ival+1
pvec(ival) = rcpro(-1, gthpot%alpha_nlcc(j))
Expand Down Expand Up @@ -1628,10 +1775,12 @@ END SUBROUTINE get_pseudo_param
!> \brief ...
!> \param pvec ...
!> \param gthpot ...
!> \param noopt_nlcc ...
! **************************************************************************************************
SUBROUTINE put_pseudo_param(pvec, gthpot)
SUBROUTINE put_pseudo_param(pvec, gthpot, noopt_nlcc)
REAL(KIND=dp), DIMENSION(:), INTENT(in) :: pvec
TYPE(atom_gthpot_type), INTENT(inout) :: gthpot
LOGICAL, INTENT(in) :: noopt_nlcc

CHARACTER(len=*), PARAMETER :: routineN = 'put_pseudo_param', &
routineP = moduleN//':'//routineN
Expand Down Expand Up @@ -1665,7 +1814,7 @@ SUBROUTINE put_pseudo_param(pvec, gthpot)
END DO
END DO
END IF
IF (gthpot%nlcc) THEN
IF (gthpot%nlcc .AND. (.NOT. noopt_nlcc)) THEN
DO j = 1, gthpot%nexp_nlcc
ival = ival+1
gthpot%alpha_nlcc(j) = rcpro(1, pvec(ival))
Expand Down
16 changes: 16 additions & 0 deletions src/input_cp2k_atom.F
View file
Expand Up @@ -1223,6 +1223,22 @@ SUBROUTINE create_powell_section(section)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="NOOPT_NLCC", &
description="Don't optimize NLCC parameters.", &
usage="NOOPT_NLCC T", &
type_of_var=logical_t, &
default_l_val=.FALSE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="PREOPT_NLCC", &
description="Optimize NLCC parameters by fitting core charge density.", &
usage="PREOPT_NLCC T", &
type_of_var=logical_t, &
default_l_val=.FALSE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

END SUBROUTINE create_powell_section

! **************************************************************************************************
Expand Down
65 changes: 65 additions & 0 deletions tests/ATOM/regtest-pseudo/Li_NLCC.inp
View file
@@ -0,0 +1,65 @@
&GLOBAL
PROGRAM_NAME ATOM
&END GLOBAL
&ATOM
ELEMENT Li

RUN_TYPE PSEUDOPOTENTIAL_OPTIMIZATION

ELECTRON_CONFIGURATION [He] 2s1
CORE [He]
MAX_ANGULAR_MOMENTUM 2

&METHOD
METHOD_TYPE KOHN-SHAM
RELATIVISTIC DKH(3)
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END METHOD
&OPTIMIZATION
EPS_SCF 1.e-8
&END

&AE_BASIS
BASIS_TYPE GEOMETRICAL_GTO
&END AE_BASIS
&PP_BASIS
BASIS_TYPE GEOMETRICAL_GTO
&END PP_BASIS
&POTENTIAL
PSEUDO_TYPE GTH
&GTH_POTENTIAL
1 0 0 0
0.62945358642262 2 -1.59552429936639 0.40867157960381
NLCC 1
0.62741295302047 1 22.95670274799660
2
0.70551855332551 1 2.20970642447945
1.12054973659941 1 0.08625487394423
&END
CONFINEMENT_TYPE BARRIER
CONFINEMENT 200. 4.0 12.0
&END POTENTIAL

&POWELL
ACCURACY 1.e-14
STEP_SIZE 0.1
MAX_INIT 1
MAX_FUN 20
STEP_SIZE_SCALING 0.90
WEIGHT_PSIR0 0.0
TARGET_POT_SEMICORE [eV] 0.0300
TARGET_POT_VALENCE [eV] 0.0030
TARGET_POT_VIRTUAL [eV] 0.0300
WEIGHT_POT_NODE 10.0
WEIGHT_POT_SEMICORE 2.0
WEIGHT_POT_VALENCE 5.0
WEIGHT_POT_VIRTUAL 1.0
SEMICORE_LEVEL [eV] 15.0
PREOPT_NLCC T
NOOPT_NLCC T
&END
&END ATOM

1 change: 1 addition & 0 deletions tests/ATOM/regtest-pseudo/TEST_FILES
View file
Expand Up @@ -16,4 +16,5 @@ C_basis1.inp 35 1.e-12
C_basis2.inp 35 1.e-12 -37.728835366125
C_basis3.inp 35 1.e-12 -13.923782735900
C_grb.inp 35 1.e-12 -5.360945665362
Li_NLCC.inp 61 1.e-12 0.5527955347
#EOF

0 comments on commit f4b239e

Please sign in to comment.