Raise error for too heavy elements (tblite#188)

Signed-off-by: Marcel Müller <[email protected]>
Co-authored-by: Sebastian Ehlert <[email protected]>

app/driver_guess.f90

@@ -115,7 +115,8 @@ subroutine guess_main(config, error)
       method = "ceh"
       if (allocated(config%method)) method = config%method
       if (method == "ceh") then
-         call new_ceh_calculator(calc_ceh, mol)
+         call new_ceh_calculator(calc_ceh, mol, error)
+         if (allocated(error)) return
          call new_wavefunction(wfn_ceh, mol%nat, calc_ceh%bas%nsh, calc_ceh%bas%nao, 1, config%etemp_guess * kt)
       end if
 

app/driver_run.f90

@@ -142,11 +142,11 @@ subroutine run_main(config, error)
       case default
          call fatal_error(error, "Unknown method '"//method//"' requested")
       case("gfn2")
-         call new_gfn2_calculator(calc, mol)
+         call new_gfn2_calculator(calc, mol, error)
       case("gfn1")
-         call new_gfn1_calculator(calc, mol)
+         call new_gfn1_calculator(calc, mol, error)
       case("ipea1")
-         call new_ipea1_calculator(calc, mol)
+         call new_ipea1_calculator(calc, mol, error)
       end select
    end if
    if (allocated(error)) return
@@ -156,7 +156,8 @@ subroutine run_main(config, error)
    call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, nspin, config%etemp * kt)
 
    if (config%guess == "ceh") then
-      call new_ceh_calculator(calc_ceh, mol)
+      call new_ceh_calculator(calc_ceh, mol, error)
+      if (allocated(error)) return
       call new_wavefunction(wfn_ceh, mol%nat, calc_ceh%bas%nsh, calc_ceh%bas%nao, 1, config%etemp_guess * kt)
       if (config%grad) then
          call ctx%message("WARNING: CEH gradient not yet implemented. Stopping.")

src/tblite/api/calculator.f90

@@ -88,6 +88,7 @@ function new_gfn2_calculator_api(vctx, vmol) result(vcalc) &
    type(vp_structure), pointer :: mol
    type(c_ptr) :: vcalc
    type(vp_calculator), pointer :: calc
+   type(error_type), allocatable :: error
 
    if (debug) print '("[Info]", 1x, a)', "new_gfn2_calculator"
 
@@ -99,8 +100,14 @@ function new_gfn2_calculator_api(vctx, vmol) result(vcalc) &
    call c_f_pointer(vmol, mol)
 
    allocate(calc)
-   call new_gfn2_calculator(calc%ptr, mol%ptr)
-   vcalc = c_loc(calc)
+   call new_gfn2_calculator(calc%ptr, mol%ptr, error)
+   if (allocated(error)) then
+      deallocate(calc)
+      call ctx%ptr%set_error(error)
+      return
+   else
+      vcalc = c_loc(calc)
+   end if
 
 end function new_gfn2_calculator_api
 
@@ -113,6 +120,7 @@ function new_ipea1_calculator_api(vctx, vmol) result(vcalc) &
    type(vp_structure), pointer :: mol
    type(c_ptr) :: vcalc
    type(vp_calculator), pointer :: calc
+   type(error_type), allocatable :: error
 
    if (debug) print '("[Info]", 1x, a)', "new_ipea1_calculator"
 
@@ -124,8 +132,14 @@ function new_ipea1_calculator_api(vctx, vmol) result(vcalc) &
    call c_f_pointer(vmol, mol)
 
    allocate(calc)
-   call new_ipea1_calculator(calc%ptr, mol%ptr)
-   vcalc = c_loc(calc)
+   call new_ipea1_calculator(calc%ptr, mol%ptr, error)
+   if (allocated(error)) then
+      deallocate(calc)
+      call ctx%ptr%set_error(error)
+      return
+   else
+      vcalc = c_loc(calc)
+   end if
 
 end function new_ipea1_calculator_api
 
@@ -138,6 +152,7 @@ function new_gfn1_calculator_api(vctx, vmol) result(vcalc) &
    type(vp_structure), pointer :: mol
    type(c_ptr) :: vcalc
    type(vp_calculator), pointer :: calc
+   type(error_type), allocatable :: error
 
    if (debug) print '("[Info]", 1x, a)', "new_gfn1_calculator"
 
@@ -149,8 +164,14 @@ function new_gfn1_calculator_api(vctx, vmol) result(vcalc) &
    call c_f_pointer(vmol, mol)
 
    allocate(calc)
-   call new_gfn1_calculator(calc%ptr, mol%ptr)
-   vcalc = c_loc(calc)
+   call new_gfn1_calculator(calc%ptr, mol%ptr, error)
+   if (allocated(error)) then
+      deallocate(calc)
+      call ctx%ptr%set_error(error)
+      return
+   else
+      vcalc = c_loc(calc)
+   end if
 
 end function new_gfn1_calculator_api
 
@@ -184,6 +205,8 @@ function new_xtb_calculator_api(vctx, vmol, vparam) result(vcalc) &
    call new_xtb_calculator(calc%ptr, mol%ptr, param%ptr, error)
    if (allocated(error)) then
       deallocate(calc)
+      call ctx%ptr%set_error(error)
+      return
    else
       vcalc = c_loc(calc)
    end if

src/tblite/ceh/ceh.f90

@@ -18,7 +18,7 @@
 !> Contains the specification of the Charge Extended Hückel (CEH) method.
 
 module tblite_ceh_ceh
-   use mctc_env, only : error_type, wp
+   use mctc_env, only : wp, error_type, fatal_error
    use mctc_io, only: structure_type
    use tblite_basis_ortho, only : orthogonalize
    use tblite_basis_slater, only : slater_to_gauss
@@ -438,10 +438,18 @@ module tblite_ceh_ceh
 contains
 
 
-   subroutine new_ceh_calculator(calc,mol)
+   subroutine new_ceh_calculator(calc, mol, error)
       !> Instance of the CEH evaluator
       type(xtb_calculator), intent(out) :: calc
       type(structure_type), intent(in)  :: mol
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+
+      ! Check if all atoms of mol%nat are supported (Z <= 86)
+      if (any(mol%num > max_elem)) then
+         call fatal_error(error, "No support for elements with Z >" // format_string(max_elem, '(i0)') // ".")
+         return
+      end if
 
       call add_ceh_basis(calc, mol)
       call add_ncoord(calc, mol)

src/tblite/xtb/gfn1.f90

@@ -19,7 +19,7 @@
 
 !> Implementation of the GFN1-xTB Hamiltonian to parametrize an xTB calculator.
 module tblite_xtb_gfn1
-   use mctc_env, only : wp
+   use mctc_env, only : wp, error_type, fatal_error
    use mctc_io, only : structure_type
    use mctc_io_symbols, only : to_symbol
    use tblite_basis_ortho, only : orthogonalize
@@ -37,6 +37,7 @@ module tblite_xtb_gfn1
    use tblite_xtb_calculator, only : xtb_calculator
    use tblite_xtb_h0, only : new_hamiltonian
    use tblite_xtb_spec, only : tb_h0spec
+   use tblite_output_format, only : format_string
    implicit none
    private
 
@@ -515,11 +516,19 @@ module tblite_xtb_gfn1
 contains
 
 
-subroutine new_gfn1_calculator(calc, mol)
+subroutine new_gfn1_calculator(calc, mol, error)
    !> Instance of the xTB evaluator
    type(xtb_calculator), intent(out) :: calc
    !> Molecular structure data
    type(structure_type), intent(in) :: mol
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   ! Check if all atoms of mol%nat are supported (Z <= 86)
+   if (any(mol%num > max_elem)) then
+      call fatal_error(error, "No support for elements with Z >" // format_string(max_elem, '(i0)') // ".")
+      return
+   end if
 
    call add_basis(calc, mol)
    call add_ncoord(calc, mol)

src/tblite/xtb/gfn2.f90

@@ -19,7 +19,7 @@
 
 !> Implementation of the GFN2-xTB Hamiltonian to parametrize an xTB calculator.
 module tblite_xtb_gfn2
-   use mctc_env, only : wp
+   use mctc_env, only : wp, error_type, fatal_error
    use mctc_io, only : structure_type
    use mctc_io_symbols, only : to_symbol
    use tblite_basis_type, only : basis_type, new_basis, cgto_type
@@ -36,6 +36,7 @@ module tblite_xtb_gfn2
    use tblite_xtb_calculator, only : xtb_calculator
    use tblite_xtb_h0, only : new_hamiltonian
    use tblite_xtb_spec, only : tb_h0spec
+   use tblite_output_format, only : format_string
    implicit none
    private
 
@@ -566,11 +567,19 @@ module tblite_xtb_gfn2
 contains
 
 
-subroutine new_gfn2_calculator(calc, mol)
+subroutine new_gfn2_calculator(calc, mol, error)
    !> Instance of the xTB evaluator
    type(xtb_calculator), intent(out) :: calc
    !> Molecular structure data
    type(structure_type), intent(in) :: mol
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   ! Check if all atoms of mol%nat are supported (Z <= 86)
+   if (any(mol%num > max_elem)) then
+      call fatal_error(error, "No support for elements with Z >" // format_string(max_elem, '(i0)') // ".")
+      return
+   end if
 
    call add_basis(calc, mol)
    call add_ncoord(calc, mol)

src/tblite/xtb/ipea1.f90

@@ -19,7 +19,7 @@
 
 !> Implementation of the IPEA1-xTB Hamiltonian to parametrize an xTB calculator.
 module tblite_xtb_ipea1
-   use mctc_env, only : wp
+   use mctc_env, only : wp, error_type, fatal_error
    use mctc_io, only : structure_type
    use mctc_io_symbols, only : to_symbol
    use tblite_basis_ortho, only : orthogonalize
@@ -37,6 +37,7 @@ module tblite_xtb_ipea1
    use tblite_xtb_calculator, only : xtb_calculator
    use tblite_xtb_h0, only : new_hamiltonian
    use tblite_xtb_spec, only : tb_h0spec
+   use tblite_output_format, only : format_string
    implicit none
    private
 
@@ -525,11 +526,19 @@ module tblite_xtb_ipea1
 contains
 
 
-subroutine new_ipea1_calculator(calc, mol)
+subroutine new_ipea1_calculator(calc, mol, error)
    !> Instance of the xTB evaluator
    type(xtb_calculator), intent(out) :: calc
    !> Molecular structure data
    type(structure_type), intent(in) :: mol
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   ! Check if all atoms of mol%nat are supported (Z <= 86)
+   if (any(mol%num > max_elem)) then
+      call fatal_error(error, "No support for elements with Z >" // format_string(max_elem, '(i0)') // ".")
+      return
+   end if
 
    call add_basis(calc, mol)
    call add_ncoord(calc, mol)

test/unit/test_ceh.f90

@@ -338,7 +338,8 @@ subroutine test_q_gen(error, mol, ref)
       integer :: i
       allocate(cn(mol%nat))
 
-      call new_ceh_calculator(calc, mol)
+      call new_ceh_calculator(calc, mol, error)
+      if (allocated(error)) return
       call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt)
       ctx%verbosity = 0
       call ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy)
@@ -1038,7 +1039,8 @@ subroutine test_q_ef_chrg_mb01(error)
 
       call get_structure(mol, "MB16-43", "01")
       mol%charge = 2.0_wp
-      call new_ceh_calculator(calc, mol)
+      call new_ceh_calculator(calc, mol, error)
+      if (allocated(error)) return
       call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt)
       cont = electric_field(efield)
       call calc%push_back(cont)
@@ -1071,7 +1073,8 @@ subroutine test_d_mb01(error)
 
       call get_structure(mol, "MB16-43", "01")
 
-      call new_ceh_calculator(calc, mol)
+      call new_ceh_calculator(calc, mol, error)
+      if (allocated(error)) return
       call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt)
       ctx%verbosity = 0
       call ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy)
@@ -1112,7 +1115,8 @@ subroutine test_d_field_mb04(error)
       efield = 0.0_wp
       efield(2) = 0.2_wp
 
-      call new_ceh_calculator(calc, mol)
+      call new_ceh_calculator(calc, mol, error)
+      if (allocated(error)) return
       call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt)
 
       cont = electric_field(efield)
@@ -1163,7 +1167,8 @@ subroutine test_d_hcn(error)
       call new(mol1, num, xyz)
       efield = 0.0_wp
       efield(1) = -0.1_wp
-      call new_ceh_calculator(calc1, mol1)
+      call new_ceh_calculator(calc1, mol1, error)
+      if (allocated(error)) return
       call new_wavefunction(wfn1, mol1%nat, calc1%bas%nsh, calc1%bas%nao, 1, kt)
       cont1 = electric_field(efield)
       call calc1%push_back(cont1)
@@ -1175,7 +1180,8 @@ subroutine test_d_hcn(error)
 
       xyz(1, :) = xyz(1, :) - 1.0_wp
       call new(mol2, num, xyz)
-      call new_ceh_calculator(calc2, mol2)
+      call new_ceh_calculator(calc2, mol2, error)
+      if (allocated(error)) return
       call new_wavefunction(wfn2, mol2%nat, calc2%bas%nsh, calc2%bas%nao, 1, kt)
       cont2 = electric_field(efield)
       call calc2%push_back(cont2)

test/unit/test_gfn1_xtb.f90

@@ -183,7 +183,8 @@ subroutine test_e_pse(error)
    do izp = 1, 86
       if (izp == 25) cycle
       call new(mol, [izp], xyz, uhf=uhf(izp))
-      call new_gfn1_calculator(calc, mol)
+      call new_gfn1_calculator(calc, mol, error)
+      if (allocated(error)) return
       call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt)
 
       energy = 0.0_wp
@@ -254,7 +255,8 @@ subroutine test_e_pse_cation(error)
    do izp = 1, 86
       if (izp == 79) cycle  ! SCF does not converge for gold
       call new(mol, [izp], xyz, uhf=uhf(izp), charge=1.0_wp)
-      call new_gfn1_calculator(calc, mol)
+      call new_gfn1_calculator(calc, mol, error)
+      if (allocated(error)) return
       call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt)
 
       energy = 0.0_wp
@@ -326,7 +328,8 @@ subroutine test_e_pse_anion(error)
       if (izp == 2) cycle  ! Helium doesn't have enough orbitals for negative charge
       if (any(izp == [21, 22, 23, 25, 43, 57, 58, 59])) cycle  ! not converging
       call new(mol, [izp], xyz, uhf=uhf(izp), charge=-1.0_wp)
-      call new_gfn1_calculator(calc, mol)
+      call new_gfn1_calculator(calc, mol, error)
+      if (allocated(error)) return
       call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt)
 
       energy = 0.0_wp
@@ -358,7 +361,8 @@ subroutine test_e_mb01(error)
 
    energy = 0.0_wp
 
-   call new_gfn1_calculator(calc, mol)
+   call new_gfn1_calculator(calc, mol, error)
+   if (allocated(error)) return
    call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt)
    call xtb_singlepoint(ctx, mol, calc, wfn, acc, energy, verbosity=0)
 
@@ -405,7 +409,8 @@ subroutine test_g_mb02(error)
    gradient(:, :) = 0.0_wp
    sigma(:, :) = 0.0_wp
 
-   call new_gfn1_calculator(calc, mol)
+   call new_gfn1_calculator(calc, mol, error)
+   if (allocated(error)) return
    call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt)
    call xtb_singlepoint(ctx, mol, calc, wfn, acc, energy, gradient, sigma, 0)
 
@@ -440,7 +445,8 @@ subroutine test_s_mb03(error)
    gradient(:, :) = 0.0_wp
    sigma(:, :) = 0.0_wp
 
-   call new_gfn1_calculator(calc, mol)
+   call new_gfn1_calculator(calc, mol, error)
+   if (allocated(error)) return
    call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt)
    call xtb_singlepoint(ctx, mol, calc, wfn, acc, energy, gradient, sigma, 0)
 
@@ -475,7 +481,8 @@ subroutine test_error_mb01(error)
 
    energy = 0.0_wp
 
-   call new_gfn1_calculator(calc, mol)
+   call new_gfn1_calculator(calc, mol, error)
+   if (allocated(error)) return
    call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt)
    call xtb_singlepoint(ctx, mol, calc, wfn, acc, energy, verbosity=0)