Remove error container from ceh_singlepoint. Add consistent timing. (t…

…blite#192)

app/driver_guess.f90

@@ -152,7 +152,7 @@ subroutine guess_main(config, error)
       case("eeq")
          call eeq_guess(mol, qat, dpat)
       case("ceh")
-         call ceh_singlepoint(ctx, calc_ceh, mol, error, wfn_ceh, config%accuracy, config%verbosity)
+         call ceh_singlepoint(ctx, calc_ceh, mol, wfn_ceh, config%accuracy, config%verbosity)
          if (ctx%failed()) then
             call fatal(ctx, "CEH singlepoint calculation failed")
             do while(ctx%failed())

app/driver_run.f90

@@ -188,7 +188,7 @@ subroutine run_main(config, error)
    case("eeq")
       call eeq_guess(mol, calc, wfn)
    case("ceh")
-      call ceh_singlepoint(ctx, calc_ceh, mol, error, wfn_ceh, config%accuracy, config%verbosity)
+      call ceh_singlepoint(ctx, calc_ceh, mol, wfn_ceh, config%accuracy, config%verbosity)
       if (ctx%failed()) then
          call fatal(ctx, "CEH singlepoint calculation failed")
          do while(ctx%failed())

src/tblite/ceh/singlepoint.f90

@@ -61,40 +61,38 @@ module tblite_ceh_singlepoint
 
 
    !> Run the CEH calculation (equivalent to xtb_singlepoint)
-   subroutine ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy, verbosity)
+   subroutine ceh_singlepoint(ctx, calc, mol, wfn, accuracy, verbosity)
       !> Calculation context
       type(context_type), intent(inout) :: ctx
       !> CEH calculator
       type(xtb_calculator), intent(inout) :: calc
       !> Molecular structure data
       type(structure_type), intent(in)  :: mol
-      !> Error container
-      type(error_type), allocatable, intent(out) :: error
       !> Wavefunction data
       type(wavefunction_type), intent(inout) :: wfn
       !> Accuracy for computation
       real(wp), intent(in) :: accuracy
       !> Verbosity level of output
       integer, intent(in), optional :: verbosity
 
-      !> Molecular dipole moment
+      ! Molecular dipole moment
       real(wp) :: dipole(3)
-      !> Integral container
+      ! Integral container
       type(integral_type) :: ints
-      !> Electronic solver
+      ! Electronic solver
       class(solver_type), allocatable :: solver
-      !> Adjacency list
+      ! Adjacency list
       type(adjacency_list) :: list
-      !> Potential type
+      ! Potential type
       type(potential_type) :: pot
-      !> Restart data for interaction containers and coulomb 
+      ! Restart data for interaction containers and coulomb 
       type(container_cache) :: icache, ccache
-      !> Timer
+      ! Timer
       type(timer_type) :: timer
-      real(wp) :: ttime
-
+      ! Error container
+      type(error_type), allocatable :: error
+
       logical :: grad
-
       real(wp) :: elec_entropy
       real(wp) :: nel, cutoff
       real(wp), allocatable :: tmp(:)
@@ -106,7 +104,7 @@ subroutine ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy, verbosity)
       ! self energy related arrays
       real(wp), allocatable :: selfenergy(:), dsedcn(:), dsedcn_en(:), lattr(:, :)
 
-      call timer%push("wall time CEH")
+      call timer%push("total CEH")
 
       if (present(verbosity)) then
          prlevel = verbosity
@@ -115,7 +113,7 @@ subroutine ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy, verbosity)
       end if
 
       if (prlevel > 1) then
-         call ctx%message("CEH guess")
+         call ctx%message("CEH singlepoint")
       endif
       ! Gradient logical as future starting point (not implemented yet)
       ! Entry point could either be (i) modified wavefunction type (including derivatives),
@@ -133,6 +131,7 @@ subroutine ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy, verbosity)
       end if
       call get_alpha_beta_occupation(wfn%nocc, wfn%nuhf, wfn%nel(1), wfn%nel(2))
 
+      call timer%push("hamiltonian")
       ! calculate coordination number (CN) and the EN-weighted coordination number
       if (allocated(calc%ncoord)) then
          allocate(cn(mol%nat))
@@ -168,6 +167,7 @@ subroutine ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy, verbosity)
       ints%quadrupole = 0.0_wp
       call get_hamiltonian(mol, lattr, list, calc%bas, calc%h0, selfenergy, &
       & ints%overlap, ints%overlap_diat, ints%dipole, ints%hamiltonian)
+      call timer%pop
 
       ! Get initial potential for external fields and Coulomb
       call new_potential(pot, mol, calc%bas, wfn%nspin)
@@ -194,6 +194,7 @@ subroutine ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy, verbosity)
       ! Add effective Hamiltonian to wavefunction
       call add_pot_to_h1(calc%bas, ints, pot, wfn%coeff)
 
+      call timer%push("diagonalization")
       ! Solve the effective Hamiltonian
       call ctx%new_solver(solver, calc%bas%nao)
 
@@ -202,6 +203,7 @@ subroutine ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy, verbosity)
       if (allocated(error)) then
          call ctx%set_error(error)
       end if
+      call timer%pop
 
       ! Get charges and dipole moment from density and integrals
       call get_mulliken_shell_charges(calc%bas, ints%overlap, wfn%density, wfn%n0sh, &
@@ -214,7 +216,27 @@ subroutine ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy, verbosity)
       dipole(:) = tmp + sum(wfn%dpat(:, :, 1), 2)
 
       call timer%pop
-      ttime = timer%get("wall time CEH")
+
+      block
+         integer :: it
+         real(wp) :: ttime, stime
+         character(len=*), parameter :: label(*) = [character(len=20):: &
+            & "coulomb", "hamiltonian", "diagonalization"]
+         if (prlevel > 1) then
+            ttime = timer%get("total CEH")
+            call ctx%message(" total CEH:"//repeat(" ", 16)//format_time(ttime))
+         end if
+         if (prlevel > 2) then
+            do it = 1, size(label)
+               stime = timer%get(label(it))
+               if (stime <= epsilon(0.0_wp)) cycle
+               call ctx%message(" - "//label(it)//format_time(stime) &
+                  & //" ("//format_string(int(stime/ttime*100), '(i3)')//"%)")
+            end do
+            call ctx%message("")
+         end if
+      end block
+
 
    end subroutine ceh_singlepoint
 

src/tblite/coulomb/charge/effective.f90

@@ -71,7 +71,7 @@ end function average_interface
    real(wp), parameter :: sqrtpi = sqrt(pi)
    real(wp), parameter :: eps = sqrt(epsilon(0.0_wp))
    real(wp), parameter :: conv = eps
-   character(len=*), parameter :: label = "isotropic Klopman-Ohno electrostatics"
+   character(len=*), parameter :: label = "isotropic Klopman-Ohno-Mataga-Nishimoto electrostatics"
 
 contains
 

test/unit/test_ceh.f90

@@ -475,7 +475,11 @@ subroutine test_q_gen(error, mol, ref)
       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)
+      call ceh_singlepoint(ctx, calc, mol, wfn, accuracy)
+      if (ctx%failed()) then
+         call ctx%get_error(error)
+         return
+      end if
 
       do i = 1, mol%nat
          call check(error, wfn%qat(i,1), ref(i), thr=1e-6_wp)
@@ -1257,7 +1261,11 @@ subroutine test_q_ef_chrg_mb01(error)
       cont = electric_field(efield)
       call calc%push_back(cont)
       ctx%verbosity = 0
-      call ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy)
+      call ceh_singlepoint(ctx, calc, mol, wfn, accuracy)
+      if (ctx%failed()) then
+         call ctx%get_error(error)
+         return
+      end if
 
       do i = 1, mol%nat
          call check(error, wfn%qat(i,1), ref(i), thr=5e-6_wp, message="Calculated charge&
@@ -1289,7 +1297,11 @@ subroutine test_d_mb01(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)
+      call ceh_singlepoint(ctx, calc, mol, wfn, accuracy)
+      if (ctx%failed()) then
+         call ctx%get_error(error)
+         return
+      end if
       tmp = 0.0_wp
       dipole = 0.0_wp
       call gemv(mol%xyz, wfn%qat(:, 1), tmp)
@@ -1335,7 +1347,11 @@ subroutine test_d_field_mb04(error)
       call calc%push_back(cont)
 
       ctx%verbosity = 0
-      call ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy)
+      call ceh_singlepoint(ctx, calc, mol, wfn, accuracy)
+      if (ctx%failed()) then
+         call ctx%get_error(error)
+         return
+      end if
       tmp = 0.0_wp
       dipole = 0.0_wp
       call gemv(mol%xyz, wfn%qat(:, 1), tmp)
@@ -1384,7 +1400,11 @@ subroutine test_d_hcn(error)
       call new_wavefunction(wfn1, mol1%nat, calc1%bas%nsh, calc1%bas%nao, 1, kt)
       cont1 = electric_field(efield)
       call calc1%push_back(cont1)
-      call ceh_singlepoint(ctx, calc1, mol1, error, wfn1, accuracy)
+      call ceh_singlepoint(ctx, calc1, mol1, wfn1, accuracy)
+      if (ctx%failed()) then
+         call ctx%get_error(error)
+         return
+      end if
       tmp = 0.0_wp
       dip1 = 0.0_wp
       call gemv(mol1%xyz, wfn1%qat(:, 1), tmp)
@@ -1397,7 +1417,11 @@ subroutine test_d_hcn(error)
       call new_wavefunction(wfn2, mol2%nat, calc2%bas%nsh, calc2%bas%nao, 1, kt)
       cont2 = electric_field(efield)
       call calc2%push_back(cont2)
-      call ceh_singlepoint(ctx, calc2, mol2, error, wfn2, accuracy)
+      call ceh_singlepoint(ctx, calc2, mol2, wfn2, accuracy)
+      if (ctx%failed()) then
+         call ctx%get_error(error)
+         return
+      end if
       tmp = 0.0_wp
       dip2 = 0.0_wp
       call gemv(mol2%xyz, wfn2%qat(:, 1), tmp)