diff --git a/src/aobasis/basis_set_types.F b/src/aobasis/basis_set_types.F index f56232b2ab..57d240a4b4 100644 --- a/src/aobasis/basis_set_types.F +++ b/src/aobasis/basis_set_types.F @@ -123,6 +123,7 @@ MODULE basis_set_types write_orb_basis_set PUBLIC :: allocate_sto_basis_set, & + read_sto_basis_set, & create_gto_from_sto_basis, & deallocate_sto_basis_set, & set_sto_basis_set, & @@ -2223,6 +2224,179 @@ CONTAINS END SUBROUTINE set_sto_basis_set +! ************************************************************************************************** +!> \brief ... +!> \param element_symbol ... +!> \param basis_set_name ... +!> \param sto_basis_set ... +!> \param para_env ... +!> \param dft_section ... +! ************************************************************************************************** + SUBROUTINE read_sto_basis_set(element_symbol, basis_set_name, sto_basis_set, para_env, dft_section) + + ! Read a Slater-type orbital (STO) basis set from the database file. + + CHARACTER(LEN=*), INTENT(IN) :: element_symbol, basis_set_name + TYPE(sto_basis_set_type), POINTER :: sto_basis_set + TYPE(cp_para_env_type), POINTER :: para_env + TYPE(section_vals_type), POINTER :: dft_section + + CHARACTER(len=*), PARAMETER :: routineN = 'read_sto_basis_set', & + routineP = moduleN//':'//routineN + + CHARACTER(LEN=10) :: nlsym + CHARACTER(LEN=2) :: lsym + CHARACTER(LEN=240) :: line + CHARACTER(LEN=242) :: line2 + CHARACTER(len=default_path_length) :: basis_set_file_name, tmp + CHARACTER(LEN=default_path_length), DIMENSION(:), & + POINTER :: cbasis + CHARACTER(LEN=LEN(basis_set_name)) :: bsname + CHARACTER(LEN=LEN(basis_set_name)+2) :: bsname2 + CHARACTER(LEN=LEN(element_symbol)) :: symbol + CHARACTER(LEN=LEN(element_symbol)+2) :: symbol2 + INTEGER :: ibasis, irep, iset, nbasis, nq, nset, & + strlen1, strlen2 + LOGICAL :: basis_found, found, match + REAL(KIND=dp) :: zet + TYPE(cp_parser_type), POINTER :: parser + + line = "" + line2 = "" + symbol = "" + symbol2 = "" + bsname = "" + bsname2 = "" + + nbasis = 1 + + sto_basis_set%name = basis_set_name + CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", & + n_rep_val=nbasis) + ALLOCATE (cbasis(nbasis)) + DO ibasis = 1, nbasis + CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", & + i_rep_val=ibasis, c_val=cbasis(ibasis)) + NULLIFY (parser) + basis_set_file_name = cbasis(ibasis) + tmp = basis_set_file_name + CALL uppercase(tmp) + END DO + + ! Search for the requested basis set in the basis set file + ! until the basis set is found or the end of file is reached + + basis_found = .FALSE. + basis_loop: DO ibasis = 1, nbasis + IF (basis_found) EXIT basis_loop + NULLIFY (parser) + basis_set_file_name = cbasis(ibasis) + CALL parser_create(parser, basis_set_file_name, para_env=para_env) + + bsname = basis_set_name + symbol = element_symbol + irep = 0 + + tmp = basis_set_name + CALL uppercase(tmp) + + IF (tmp .NE. "NONE") THEN + search_loop: DO + + CALL parser_search_string(parser, TRIM(bsname), .TRUE., found, line) + IF (found) THEN + CALL uppercase(symbol) + CALL uppercase(bsname) + + match = .FALSE. + CALL uppercase(line) + ! Check both the element symbol and the basis set name + line2 = " "//line//" " + symbol2 = " "//TRIM(symbol)//" " + bsname2 = " "//TRIM(bsname)//" " + strlen1 = LEN_TRIM(symbol2)+1 + strlen2 = LEN_TRIM(bsname2)+1 + + IF ((INDEX(line2, symbol2(:strlen1)) > 0) .AND. & + (INDEX(line2, bsname2(:strlen2)) > 0)) match = .TRUE. + IF (match) THEN + ! Read the basis set information + CALL parser_get_object(parser, nset, newline=.TRUE.) + sto_basis_set%nshell = nset + + CALL reallocate(sto_basis_set%nq, 1, nset) + CALL reallocate(sto_basis_set%lq, 1, nset) + CALL reallocate(sto_basis_set%zet, 1, nset) + ALLOCATE (sto_basis_set%symbol(nset)) + + DO iset = 1, nset + CALL parser_get_object(parser, nq, newline=.TRUE.) + CALL parser_get_object(parser, lsym) + CALL parser_get_object(parser, zet) + sto_basis_set%nq(iset) = nq + sto_basis_set%zet(iset) = zet + WRITE (nlsym, "(I2,A)") nq, TRIM(lsym) + sto_basis_set%symbol(iset) = TRIM(nlsym) + SELECT CASE (TRIM (lsym)) + CASE ("S", "s") + sto_basis_set%lq(iset) = 0 + CASE ("P", "p") + sto_basis_set%lq(iset) = 1 + CASE ("D", "d") + sto_basis_set%lq(iset) = 2 + CASE ("F", "f") + sto_basis_set%lq(iset) = 3 + CASE ("G", "g") + sto_basis_set%lq(iset) = 4 + CASE ("H", "h") + sto_basis_set%lq(iset) = 5 + CASE ("I", "i", "J", "j") + sto_basis_set%lq(iset) = 6 + CASE ("K", "k") + sto_basis_set%lq(iset) = 7 + CASE ("L", "l") + sto_basis_set%lq(iset) = 8 + CASE ("M", "m") + sto_basis_set%lq(iset) = 9 + CASE DEFAULT + CALL cp_abort(__LOCATION__, & + "The requested basis set <"//TRIM(bsname)// & + "> for element <"//TRIM(symbol)//"> has an ivalid component: ") + END SELECT + END DO + + basis_found = .TRUE. + EXIT search_loop + END IF + ELSE + EXIT search_loop + END IF + END DO search_loop + ELSE + match = .FALSE. + ENDIF + + CALL parser_release(parser) + + END DO basis_loop + + IF (tmp .NE. "NONE") THEN + IF (.NOT. basis_found) THEN + basis_set_file_name = "" + DO ibasis = 1, nbasis + basis_set_file_name = TRIM(basis_set_file_name)//"<"//TRIM(cbasis(ibasis))//"> " + END DO + CALL cp_abort(__LOCATION__, & + "The requested basis set <"//TRIM(bsname)// & + "> for element <"//TRIM(symbol)//"> was not "// & + "found in the basis set files "// & + TRIM(basis_set_file_name)) + END IF + END IF + DEALLOCATE (cbasis) + + END SUBROUTINE read_sto_basis_set + ! ************************************************************************************************** !> \brief ... !> \param sto_basis_set ... diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 8fab9685c3..b89eb02db2 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -2961,6 +2961,12 @@ CONTAINS CALL keyword_release(keyword) ! Integers + CALL keyword_create(keyword, name="STO_NG", & + description="Order of Gaussian type expansion of Slater orbital basis sets.", & + usage="STO_NG", default_i_val=6) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, name="LMAXN1", & variants=(/"LMAXRHO1"/), & description="GAPW : max L number for expansion of the atomic densities in spherical gaussians", & diff --git a/src/input_cp2k_subsys.F b/src/input_cp2k_subsys.F index 7fefdc640c..641de1453b 100644 --- a/src/input_cp2k_subsys.F +++ b/src/input_cp2k_subsys.F @@ -1064,8 +1064,11 @@ CONTAINS CALL keyword_release(keyword) CALL keyword_create(keyword, name="BASIS_SET", & - description="The primary Gaussian basis set (NONE implies no basis used, meaningful with GHOST)", & - usage="BASIS_SET [type] DZVP", type_of_var=char_t, default_c_vals=(/" ", " "/), & + description="The primary Gaussian basis set (NONE implies no basis used, meaningful with GHOST)"// & + "Defaults are set for TYPE {ORB} and FORM {GTO}. Possible values for TYPE are {ORB, AUX, RI_AUX, LRI, ...}. "// & + "Possible values for FORM are {GTO, STO}. Where STO results in a GTO expansion of a Slater type basis."// & + "If a value for FORM is given, also TYPE has to be set explicitely.", & + usage="BASIS_SET [type] [form] DZVP", type_of_var=char_t, default_c_vals=(/" ", " ", " "/), & citations=(/VandeVondele2005a, VandeVondele2007/), & repeats=.TRUE., n_var=-1) CALL section_add_keyword(section, keyword) diff --git a/src/qs_kind_types.F b/src/qs_kind_types.F index c0b4498cc0..dc59096ffc 100644 --- a/src/qs_kind_types.F +++ b/src/qs_kind_types.F @@ -35,9 +35,11 @@ MODULE qs_kind_types remove_basis_from_container,& remove_basis_set_container USE basis_set_types, ONLY: & - allocate_gto_basis_set, combine_basis_sets, deallocate_gto_basis_set, get_gto_basis_set, & - gto_basis_set_type, init_aux_basis_set, init_orb_basis_set, read_gto_basis_set, & - write_gto_basis_set, write_orb_basis_set + allocate_gto_basis_set, allocate_sto_basis_set, combine_basis_sets, & + create_gto_from_sto_basis, deallocate_gto_basis_set, deallocate_sto_basis_set, & + get_gto_basis_set, gto_basis_set_type, init_aux_basis_set, init_orb_basis_set, & + read_gto_basis_set, read_sto_basis_set, sto_basis_set_type, write_gto_basis_set, & + write_orb_basis_set USE cp_control_types, ONLY: dft_control_type,& qs_control_type USE cp_log_handling, ONLY: cp_get_default_logger,& @@ -1372,7 +1374,8 @@ CONTAINS CHARACTER(LEN=default_string_length), & DIMENSION(:), POINTER :: tmpstringlist CHARACTER(LEN=default_string_length), & - DIMENSION(maxbas) :: basis_set_name, basis_set_type + DIMENSION(maxbas) :: basis_set_form, basis_set_name, & + basis_set_type INTEGER :: handle, i, i_rep, ipaodesc, ipaopot, ipos, j, jj, k_rep, l, m, n_rep, nb_rep, & nexp, ngauss, nlcc, nloc, nnl, norbitals, npaodesc, npaopot, nppnl, nspin, z INTEGER, DIMENSION(:), POINTER :: add_el, elec_conf, orbitals @@ -1390,6 +1393,7 @@ CONTAINS TYPE(section_vals_type), POINTER :: basis_section, bs_section, dft_plus_u_section, & dft_section, enforce_occupation_section, kgpot_section, pao_desc_section, & pao_pot_section, potential_section, spin_section + TYPE(sto_basis_set_type), POINTER :: sto_basis_set CALL timeset(routineN, handle) @@ -1398,6 +1402,7 @@ CONTAINS update_input = .TRUE. basis_set_name(:) = "" basis_set_type(:) = "" + basis_set_form(:) = "" aux_basis_set_name = "" ri_aux_basis_set_name = "" orb_basis_set_name = "" @@ -1494,15 +1499,26 @@ CONTAINS CALL section_vals_val_get(kind_section, i_rep_section=k_rep, & keyword_name="BASIS_SET", i_rep_val=i, c_vals=tmpstringlist) IF (SIZE(tmpstringlist) == 1) THEN - ! default is orbital type + ! default is orbital type and GTO basis_set_type(i) = "ORB" + basis_set_form(i) = "GTO" basis_set_name(i) = tmpstringlist(1) ELSEIF (SIZE(tmpstringlist) == 2) THEN + ! default is GTO basis_set_type(i) = tmpstringlist(1) + basis_set_form(i) = "GTO" basis_set_name(i) = tmpstringlist(2) + ELSEIF (SIZE(tmpstringlist) == 3) THEN + basis_set_type(i) = tmpstringlist(1) + basis_set_form(i) = tmpstringlist(2) + basis_set_name(i) = tmpstringlist(3) ELSE CPABORT("") END IF + ! check that we have a valid basis set form + IF (basis_set_form(i) /= "GTO" .AND. basis_set_form(i) /= "STO") THEN + CPABORT("BASIS_SET_FORM invalid") + END IF END DO ! Deprecated keywords @@ -1540,6 +1556,7 @@ CONTAINS nb_rep = nb_rep+1 CPASSERT(nb_rep <= maxbas) basis_set_type(nb_rep) = "RI_AUX" + basis_set_form(nb_rep) = "GTO" basis_set_name(nb_rep) = ri_aux_basis_set_name END IF ! end of deprecated input @@ -1975,13 +1992,30 @@ CONTAINS CASE DEFAULT + ! set nauss for STO expansion + CALL section_vals_val_get(dft_section, "QS%STO_NG", i_val=ngauss) ! Allocate and initialise the basis set data set structure ! first external basis sets DO i = 1, nb_rep - NULLIFY (tmp_basis_set) - CALL allocate_gto_basis_set(tmp_basis_set) - CALL read_gto_basis_set(qs_kind%element_symbol, basis_set_name(i), & - tmp_basis_set, para_env, dft_section) + SELECT CASE (basis_set_form (i)) + CASE ("GTO") + NULLIFY (tmp_basis_set) + CALL allocate_gto_basis_set(tmp_basis_set) + CALL read_gto_basis_set(qs_kind%element_symbol, basis_set_name(i), & + tmp_basis_set, para_env, dft_section) + CASE ("STO") + NULLIFY (sto_basis_set) + CALL allocate_sto_basis_set(sto_basis_set) + CALL read_sto_basis_set(qs_kind%element_symbol, basis_set_name(i), & + sto_basis_set, para_env, dft_section) + NULLIFY (tmp_basis_set) + CALL create_gto_from_sto_basis(sto_basis_set, tmp_basis_set, ngauss) + CALL deallocate_sto_basis_set(sto_basis_set) + CASE DEFAULT + CALL cp_abort(__LOCATION__, & + "Invalid basis set form "//TRIM(basis_set_form(i))// & + "for atomic kind <"//TRIM(qs_kind%name)//">") + END SELECT tmp = basis_set_type(i) CALL uppercase(tmp) CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, tmp) diff --git a/tests/QS/regtest-sto/BASIS_STO b/tests/QS/regtest-sto/BASIS_STO new file mode 100644 index 0000000000..827205c233 --- /dev/null +++ b/tests/QS/regtest-sto/BASIS_STO @@ -0,0 +1,44 @@ +H DZ_AE + 2 + 1 S 0.760 + 1 S 1.280 + +H QZ_AE + 8 + 1 S 3.300 + 1 S 2.000 + 1 S 1.400 + 1 S 1.000 + 1 S 0.710 + 2 P 2.000 + 2 P 1.000 + 3 D 2.500 + 3 D 1.500 + +O DZ_AE + 6 + 1 S 9.80 + 1 S 6.80 + 2 S 1.70 + 2 S 2.82 + 2 P 3.06 + 2 P 1.30 + +O QZ_AE + 16 + 1 S 14.800 + 1 S 8.750 + 1 S 6.650 + 2 S 7.100 + 2 S 3.200 + 2 S 2.050 + 2 S 1.500 + 2 S 0.750 + 2 P 5.700 + 2 P 3.050 + 2 P 1.650 + 2 P 1.000 + 3 D 2.500 + 3 D 1.500 + 4 F 4.000 + 4 F 2.000 diff --git a/tests/QS/regtest-sto/H2O_t1.inp b/tests/QS/regtest-sto/H2O_t1.inp new file mode 100644 index 0000000000..94d370870f --- /dev/null +++ b/tests/QS/regtest-sto/H2O_t1.inp @@ -0,0 +1,50 @@ + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME ./BASIS_STO + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 250 + NGRIDS 5 + &END MGRID + &QS + METHOD GAPW + STO_NG 4 + &END QS + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &OT ON + MINIMIZER DIIS + PRECONDITIONER FULL_ALL + &END OT + &END SCF + &XC + &XC_FUNCTIONAL Pade + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + &COORD + O 0.000000 0.000000 -0.065587 + H 0.000000 -0.757136 0.520545 + H 0.000000 0.757136 0.520545 + &END COORD + &KIND H + BASIS_SET ORB STO DZ_AE + POTENTIAL ALL + &END KIND + &KIND O + BASIS_SET ORB STO DZ_AE + POTENTIAL ALL + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT sto1 + PRINT_LEVEL MEDIUM +&END GLOBAL diff --git a/tests/QS/regtest-sto/TEST_FILES b/tests/QS/regtest-sto/TEST_FILES new file mode 100644 index 0000000000..8e308d37ac --- /dev/null +++ b/tests/QS/regtest-sto/TEST_FILES @@ -0,0 +1,8 @@ +# runs are executed in the same order as in this file +# the second field tells which test should be run in order to compare with the last available output +# e.g. 0 means do not compare anything, running is enough +# 1 compares the last total energy in the file +# for details see cp2k/tools/do_regtest +# tests GAPW STO +H2O_t1.inp 1 1e-13 -75.659689099442119 +#EOF diff --git a/tests/QS/regtest-sto/TEST_FILES_RESET b/tests/QS/regtest-sto/TEST_FILES_RESET new file mode 100644 index 0000000000..792d600548 --- /dev/null +++ b/tests/QS/regtest-sto/TEST_FILES_RESET @@ -0,0 +1 @@ +#