diff --git a/src/atom_fit.F b/src/atom_fit.F index fcebbfe510..2ee078efa0 100644 --- a/src/atom_fit.F +++ b/src/atom_fit.F @@ -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) diff --git a/src/atom_output.F b/src/atom_output.F index 4c29b6ac79..bcfb220ee3 100644 --- a/src/atom_output.F +++ b/src/atom_output.F @@ -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 diff --git a/src/atom_pseudo.F b/src/atom_pseudo.F index 856638fb51..8ee6a4dc21 100644 --- a/src/atom_pseudo.F +++ b/src/atom_pseudo.F @@ -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 diff --git a/src/input_cp2k_atom.F b/src/input_cp2k_atom.F index 3e7607ca97..751755f87a 100644 --- a/src/input_cp2k_atom.F +++ b/src/input_cp2k_atom.F @@ -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 diff --git a/tests/ATOM/regtest-pseudo/H-ECPGTH.inp b/tests/ATOM/regtest-pseudo/H-ECPGTH.inp new file mode 100644 index 0000000000..8f16d03655 --- /dev/null +++ b/tests/ATOM/regtest-pseudo/H-ECPGTH.inp @@ -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 + >H_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 diff --git a/tests/ATOM/regtest-pseudo/TEST_FILES.toml b/tests/ATOM/regtest-pseudo/TEST_FILES.toml index cab59403de..e8e50e4deb 100644 --- a/tests/ATOM/regtest-pseudo/TEST_FILES.toml +++ b/tests/ATOM/regtest-pseudo/TEST_FILES.toml @@ -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