diff --git a/src/input_cp2k_subsys.F b/src/input_cp2k_subsys.F index 7d7974cf5a..f03e0ac138 100644 --- a/src/input_cp2k_subsys.F +++ b/src/input_cp2k_subsys.F @@ -1199,6 +1199,12 @@ CONTAINS CALL section_add_keyword(section,keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, name="LMAX_DFTB",& + description="The maximum l-quantum number of the DFTB basis for this kind.",& + usage="LMAX_DFTB 1", default_i_val=-1) + CALL section_add_keyword(section,keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, name="MAO",& description="The number of MAOs (Modified Atomic Orbitals) for this kind.",& usage="MAO 4", default_i_val=-1) diff --git a/src/qs_dftb_matrices.F b/src/qs_dftb_matrices.F index 7d0b42d4b3..34a28f9206 100644 --- a/src/qs_dftb_matrices.F +++ b/src/qs_dftb_matrices.F @@ -312,7 +312,6 @@ CONTAINS CALL get_dftb_atom_param(dftb_kind_a,& defined=defined,lmax=lmaxi,skself=skself,& eta=eta_a,natorb=natorb_a) - IF (.NOT.defined .OR. natorb_a < 1) CYCLE CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b) CALL get_dftb_atom_param(dftb_kind_b,& diff --git a/src/qs_dftb_parameters.F b/src/qs_dftb_parameters.F index d8dc3d1550..9b4a87a1c1 100644 --- a/src/qs_dftb_parameters.F +++ b/src/qs_dftb_parameters.F @@ -261,6 +261,9 @@ CONTAINS IF ( ABS(uwork(k)) >= 1.e-12_dp ) n_urpoly = k END DO END IF +! Polynomials of length 1 are not allowed, it seems we should use spline after all +! This is creative guessing! + IF ( n_urpoly < 2 ) n_urpoly = 0 END IF CALL mp_bcast(n_urpoly,para_env%source,para_env%group) @@ -299,27 +302,37 @@ CONTAINS ! In the DFTB-Slako convention they are on orbital 10 (s-s-sigma), ! 7 (p-p-sigma) and 3 (d-d-sigma). ! - lmax=0 - DO l=0,3 - SELECT CASE (l) - CASE DEFAULT - CPABORT("") - CASE (0) - lp = 10 - CASE (1) - lp = 7 - CASE (2) - lp = 3 - CASE (3) - lp = 3 ! this is wrong but we don't allow f anyway - END SELECT - ! Technical note: In some slako files dummies are included in the - ! first matrix elements, so remove them. - IF ( (ABS(skself(l)) > 0._dp) .OR. & - (SUM(ABS(fmat(ngrd/10:ngrd,lp))) > 0._dp) ) lmax=l - END DO - ! l=2 (d) is maximum - lmax = MIN ( 2, lmax ) + ! We also allow lmax to be set in the input (in KIND) + ! + CALL get_qs_kind(qs_kind_set(ikind),lmax_dftb=lmax) + IF ( lmax < 0 ) THEN + lmax=0 + DO l=0,3 + SELECT CASE (l) + CASE DEFAULT + CPABORT("") + CASE (0) + lp = 10 + CASE (1) + lp = 7 + CASE (2) + lp = 3 + CASE (3) + lp = 3 ! this is wrong but we don't allow f anyway + END SELECT + ! Technical note: In some slako files dummies are included in the + ! first matrix elements, so remove them. + IF ( (ABS(skself(l)) > 0._dp) .OR. & + (SUM(ABS(fmat(ngrd/10:ngrd,lp))) > 0._dp) ) lmax=l + END DO + ! l=2 (d) is maximum + lmax = MIN(2,lmax) + END IF + IF ( lmax > 2 ) THEN + CALL cp_abort(__LOCATION__,"Maximum L allowed is d. "//& + "Use KIND/LMAX_DFTB to set smaller values if needed.") + END IF + ! CALL set_dftb_atom_param(dftb_parameter=dftb_atom_a,& lmax=lmax, natorb=(lmax+1)**2) @@ -461,6 +474,9 @@ CONTAINS IF ( ABS(uwork(k)) >= 1.e-12_dp ) n_urpoly = k END DO END IF +! Polynomials of length 1 are not allowed, it seems we should use spline after all +! This is creative guessing! + IF ( n_urpoly < 2 ) n_urpoly = 0 END IF CALL mp_bcast(n_urpoly,para_env%source,para_env%group) diff --git a/src/qs_kind_types.F b/src/qs_kind_types.F index 13eba389a7..98566300e4 100644 --- a/src/qs_kind_types.F +++ b/src/qs_kind_types.F @@ -138,6 +138,7 @@ MODULE qs_kind_types LOGICAL :: paw_atom = .FALSE. ! needs atomic rho1 LOGICAL :: gpw_type_forced = .FALSE. ! gpw atom even if with hard exponents LOGICAL :: ghost = .FALSE. + INTEGER :: lmax_dftb = -1 REAL(KIND = dp) :: dudq_dftb3 = 0.0_dp INTEGER, DIMENSION(:,:), POINTER :: addel => Null() INTEGER, DIMENSION(:,:), POINTER :: laddel => Null() @@ -299,6 +300,7 @@ CONTAINS !> \param zeff ... !> \param elec_conf ... !> \param mao ... +!> \param lmax_dftb ... !> \param alpha_core_charge ... !> \param ccore_charge ... !> \param core_charge ... @@ -349,7 +351,7 @@ CONTAINS basis_set, basis_type, ncgf, nsgf, & all_potential, tnadd_potential, gth_potential, & se_parameter, dftb_parameter, scptb_parameter, & - dftb3_param, zeff, elec_conf, mao,& + dftb3_param, zeff, elec_conf, mao, lmax_dftb, & alpha_core_charge, ccore_charge, core_charge, core_charge_radius,& soft_basis_set, hard_basis_set, paw_proj_set, softb, & paw_atom, hard_radius, hard0_radius, max_rad_local, & @@ -380,7 +382,7 @@ CONTAINS POINTER :: scptb_parameter REAL(KIND=dp), INTENT(OUT), OPTIONAL :: dftb3_param, zeff INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf - INTEGER, INTENT(OUT), OPTIONAL :: mao + INTEGER, INTENT(OUT), OPTIONAL :: mao, lmax_dftb REAL(KIND=dp), INTENT(OUT), OPTIONAL :: alpha_core_charge, & ccore_charge, core_charge, & core_charge_radius @@ -657,6 +659,8 @@ CONTAINS IF (PRESENT(mao)) mao = qs_kind%mao + IF (PRESENT(lmax_dftb)) lmax_dftb = qs_kind%lmax_dftb + IF (PRESENT(pao_basis_size)) pao_basis_size = qs_kind%pao_basis_size IF (PRESENT(pao_potential_maxl)) pao_potential_maxl = qs_kind%pao_potential_maxl IF (PRESENT(pao_potential_neighbors)) pao_potential_neighbors = qs_kind%pao_potential_neighbors @@ -1414,6 +1418,8 @@ CONTAINS ! DFTB3 param CALL section_vals_val_get(kind_section,i_rep_section=k_rep,& keyword_name="DFTB3_PARAM",r_val=qs_kind%dudq_dftb3) + CALL section_vals_val_get(kind_section,i_rep_section=k_rep,& + keyword_name="LMAX_DFTB",i_val=qs_kind%lmax_dftb) ! MAOS CALL section_vals_val_get(kind_section,i_rep_section=k_rep,& diff --git a/tests/DFTB/regtest-nonscc/TEST_FILES b/tests/DFTB/regtest-nonscc/TEST_FILES index ca52184f57..defaa63569 100644 --- a/tests/DFTB/regtest-nonscc/TEST_FILES +++ b/tests/DFTB/regtest-nonscc/TEST_FILES @@ -23,6 +23,7 @@ h2o-32_atprop.inp 1 6e-13 -131.08636569540289 # co2_1.inp 1 1.0E-14 -8.55586246566686 co2_2.inp 1 1.0E-14 -8.55586246566686 +co2_3.inp 1 1.0E-14 -8.55586246566686 # si_kp1.inp 1 1.0E-13 -10.063582034753917 si_kp2.inp 0 diff --git a/tests/DFTB/regtest-nonscc/co2_3.inp b/tests/DFTB/regtest-nonscc/co2_3.inp new file mode 100644 index 0000000000..1f0c53ded9 --- /dev/null +++ b/tests/DFTB/regtest-nonscc/co2_3.inp @@ -0,0 +1,49 @@ +#CPQA INCLUDE DFTB/nonscc/nonscc_parameter +#CPQA INCLUDE uff_table +#CPQA INCLUDE DFTB/nonscc/hh +&FORCE_EVAL + &DFT + &QS + METHOD DFTB + &DFTB + SELF_CONSISTENT F + &PARAMETER + PARAM_FILE_PATH DFTB/nonscc + SK_FILE O O oo + SK_FILE C C cc + SK_FILE C O co + SK_FILE O C oc + &END PARAMETER + &END DFTB + &END QS + &SCF + SCF_GUESS NONE + &MIXING + METHOD DIRECT_P_MIXING + ALPHA 1. + &END + &END SCF + &END DFT + &SUBSYS + &KIND O + LMAX_DFTB 1 + &END KIND + &KIND C + LMAX_DFTB 1 + &END KIND + &CELL + ABC 20.0 20.0 20.0 + PERIODIC NONE + &END CELL + &COORD + C 0.0 0.0 0.0 + O 0.0 0.0 +1.2 + O 0.0 0.0 -1.2 + &END COORD + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT coo + RUN_TYPE ENERGY_FORCE + PRINT_LEVEL HIGH +&END GLOBAL