Electrostatics and second version parameters for CEH (tblite#190)

Inclusion of approximate electrostatics based on the EN-weighted CN as charges. Update of the parameters and tests to the second version of CEH.

app/cli.f90

@@ -75,7 +75,7 @@ module tblite_cli
       !> Electronic temperature
       real(wp) :: etemp = 300.0_wp
       !> Electronic temperature for the guess (currently only CEH)
-      real(wp) :: etemp_guess = 5000.0_wp
+      real(wp) :: etemp_guess = 4000.0_wp
       !> Electric field
       real(wp), allocatable :: efield(:)
       !> Spin polarization
@@ -107,7 +107,7 @@ module tblite_cli
       !> File for output of JSON dump
       character(len=:), allocatable :: json_output
       !> Electronic temperature for the guess (currently only CEH)
-      real(wp) :: etemp_guess = 5000.0_wp
+      real(wp) :: etemp_guess = 4000.0_wp
       !> Electric field
       real(wp), allocatable :: efield(:)
       !> Algorithm for electronic solver

src/tblite/ceh/singlepoint.f90

@@ -33,12 +33,13 @@ module tblite_ceh_singlepoint
    use tblite_wavefunction_mulliken, only: get_mulliken_shell_charges, &
    & get_mulliken_atomic_multipoles
    use tblite_scf_iterator, only: get_density, get_qat_from_qsh
-   use tblite_scf, only: new_potential, potential_type ! Potential for external field
+   use tblite_scf, only: new_potential, potential_type 
    use tblite_container, only : container_cache
    use tblite_scf_potential, only: add_pot_to_h1
    use tblite_scf_solver, only : solver_type
    use tblite_blas, only : gemv
    use tblite_ceh_h0, only : get_hamiltonian, get_scaled_selfenergy, get_occupation
+   use tblite_ceh_ceh, only : get_effective_qat
    use tblite_xtb_spec, only : tb_h0spec 
    use tblite_xtb_calculator, only : xtb_calculator
    use tblite_timer, only : timer_type, format_time
@@ -54,7 +55,7 @@ module tblite_ceh_singlepoint
    character(len=25), parameter :: &
       label_cutoff = "integral cutoff", &
       label_charges = "CEH atomic charges", &
-      label_dipole = "CEH molecular dipole moment / a.u."
+      label_dipole = "CEH mol. dip. mom. / a.u."
 
 contains
 
@@ -86,8 +87,8 @@ subroutine ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy, verbosity)
       type(adjacency_list) :: list
       !> Potential type
       type(potential_type) :: pot
-      !> Restart data for interaction containers
-      type(container_cache) :: icache
+      !> Restart data for interaction containers and coulomb 
+      type(container_cache) :: icache, ccache
       !> Timer
       type(timer_type) :: timer
       real(wp) :: ttime
@@ -168,14 +169,27 @@ subroutine ceh_singlepoint(ctx, calc, mol, error, wfn, accuracy, verbosity)
       call get_hamiltonian(mol, lattr, list, calc%bas, calc%h0, selfenergy, &
       & ints%overlap, ints%overlap_diat, ints%dipole, ints%hamiltonian)
 
-      ! Get initial potential
+      ! Get initial potential for external fields and Coulomb
       call new_potential(pot, mol, calc%bas, wfn%nspin)
       ! Set potential to zero
       call pot%reset
+      ! Add potential due to external field
       if (allocated(calc%interactions)) then
+         call timer%push("interactions")
          call calc%interactions%update(mol, icache)
          call calc%interactions%get_potential(mol, icache, wfn, pot)
+         call timer%pop
       endif
+      ! Add potential due to Coulomb
+      if (allocated(calc%coulomb)) then
+         call timer%push("coulomb")
+         ! Use electronegativity-weighted CN as 0th-order charge guess
+         call get_effective_qat(mol, calc%bas, cn_en, wfn%qat)
+
+         call calc%coulomb%update(mol, ccache)
+         call calc%coulomb%get_potential(mol, ccache, wfn, pot)
+         call timer%pop
+      end if
 
       ! Add effective Hamiltonian to wavefunction
       call add_pot_to_h1(calc%bas, ints, pot, wfn%coeff)

src/tblite/ncoord/erf.f90

@@ -40,7 +40,7 @@ module tblite_ncoord_erf
    end type erf_ncoord_type
 
    !> Steepness of counting function (CEH)
-   real(wp), parameter :: default_kcn = 2.60_wp
+   real(wp), parameter :: default_kcn = 3.15_wp
 
    real(wp), parameter :: default_cutoff = 25.0_wp
 

src/tblite/ncoord/erf_en.f90

@@ -42,7 +42,7 @@ module tblite_ncoord_erf_en
    end type erf_en_ncoord_type
 
    !> Steepness of counting function (CEH)
-   real(wp), parameter :: default_kcn = 2.60_wp
+   real(wp), parameter :: default_kcn = 2.65_wp
 
    real(wp), parameter :: default_cutoff = 25.0_wp
 

src/tblite/wavefunction/type.f90

@@ -145,5 +145,4 @@ subroutine get_alpha_beta_occupation(nocc, nuhf, nalp, nbet)
    nbet = ntmp / 2
 end subroutine get_alpha_beta_occupation
 
-
 end module tblite_wavefunction_type

test/unit/test_ncoord.f90

@@ -752,6 +752,7 @@ subroutine test_cn_mb01_erf(error)
       real(wp), allocatable :: rcov(:)
 
       real(wp), parameter :: cutoff = 30.0_wp
+      real(wp), parameter :: kcn = 2.60_wp
       real(wp), parameter :: ref(16) = [&
       & +4.02932551939856E+00_wp, +7.85103341426644E-01_wp, +1.81014102249862E+00_wp, &
       & +1.27570668822915E+00_wp, +1.07819362139425E+00_wp, +9.25571697114998E-01_wp, &
@@ -765,7 +766,7 @@ subroutine test_cn_mb01_erf(error)
       allocate(rcov(mol%nid))
       rcov(:) = get_covalent_rad(mol%num)
 
-      call new_erf_ncoord(erf_ncoord, mol, cutoff=cutoff, rcov=rcov)
+      call new_erf_ncoord(erf_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov)
       call test_cn_gen(error, mol, erf_ncoord, cutoff, ref)
 
    end subroutine test_cn_mb01_erf
@@ -781,6 +782,7 @@ subroutine test_cn_mb02_erf(error)
       real(wp), allocatable :: rcov(:)
 
       real(wp), parameter :: cutoff = 30.0_wp
+      real(wp), parameter :: kcn = 2.60_wp
       real(wp), parameter :: ref(16) = [&
       & +7.69338455488404E-01_wp, +3.43906282200685E+00_wp, +3.09682161007250E+00_wp, &
       & +2.51887440744333E+00_wp, +4.48473180889284E+00_wp, +1.04427988033896E+00_wp, &
@@ -794,7 +796,7 @@ subroutine test_cn_mb02_erf(error)
       allocate(rcov(mol%nid))
       rcov(:) = get_covalent_rad(mol%num)
 
-      call new_erf_ncoord(erf_ncoord, mol, cutoff=cutoff, rcov=rcov)
+      call new_erf_ncoord(erf_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov)
       call test_cn_gen(error, mol, erf_ncoord, cutoff, ref)
 
    end subroutine test_cn_mb02_erf
@@ -810,6 +812,7 @@ subroutine test_cn_mb03_erf(error)
       real(wp), allocatable :: rcov(:)
 
       real(wp), parameter :: cutoff = 30.0_wp
+      real(wp), parameter :: kcn = 2.60_wp
       real(wp), parameter :: ref(16) = [&
       & +3.65604858453391E+00_wp, 2.47444042434446E+00_wp, 1.04729945244460E+00_wp, &
       & +4.69028664318902E+00_wp, 6.13331280895969E+00_wp, 3.96706985488549E+00_wp, &
@@ -823,7 +826,7 @@ subroutine test_cn_mb03_erf(error)
       allocate(rcov(mol%nid))
       rcov(:) = get_covalent_rad(mol%num)
 
-      call new_erf_ncoord(erf_ncoord, mol, cutoff=cutoff, rcov=rcov)
+      call new_erf_ncoord(erf_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov)
       call test_cn_gen(error, mol, erf_ncoord, cutoff, ref)
 
    end subroutine test_cn_mb03_erf
@@ -839,6 +842,7 @@ subroutine test_cn_acetic_erf(error)
       real(wp), allocatable :: rcov(:)
 
       real(wp), parameter :: cutoff = 30.0_wp
+      real(wp), parameter :: kcn = 2.60_wp
       real(wp), parameter :: ref(32) = [&
       & +1.04618537521867E+00_wp, 1.04621440498369E+00_wp, 1.04615269011270E+00_wp, &
       & +1.04619270680968E+00_wp, 8.39537394573134E-01_wp, 8.39512679829574E-01_wp, &
@@ -857,7 +861,7 @@ subroutine test_cn_acetic_erf(error)
       allocate(rcov(mol%nid))
       rcov(:) = get_covalent_rad(mol%num)
 
-      call new_erf_ncoord(erf_ncoord, mol, cutoff=cutoff, rcov=rcov)
+      call new_erf_ncoord(erf_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov)
       call test_cn_gen(error, mol, erf_ncoord, cutoff, ref)
 
    end subroutine test_cn_acetic_erf
@@ -1011,6 +1015,7 @@ subroutine test_cn_mb01_erf_en(error)
       real(wp), allocatable :: en(:)
 
       real(wp), parameter :: cutoff = 30.0_wp
+      real(wp), parameter :: kcn = 2.60_wp
       real(wp), parameter :: ref(16) = [&
       & +6.49835771218084E+00_wp, -9.24780426380943E-02_wp, -3.19404023806046E+00_wp, &
       & -7.26463405313327E-01_wp, -1.92171552979195E+00_wp, +6.71703342153316E-01_wp, &
@@ -1025,7 +1030,7 @@ subroutine test_cn_mb01_erf_en(error)
       rcov(:) = get_covalent_rad(mol%num)
       en(:) = get_pauling_en(mol%num)
 
-      call new_erf_en_ncoord(erf_en_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en)
+      call new_erf_en_ncoord(erf_en_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov, en=en)
       call test_cn_gen(error, mol, erf_en_ncoord, cutoff, ref)
 
    end subroutine test_cn_mb01_erf_en
@@ -1042,6 +1047,7 @@ subroutine test_cn_mb02_erf_en(error)
       real(wp), allocatable :: en(:)
 
       real(wp), parameter :: cutoff = 30.0_wp
+      real(wp), parameter :: kcn = 2.60_wp
       real(wp), parameter :: ref(16) = [&
       & -1.12234584286121E-01_wp, -2.94004224176729E+00_wp, -5.58562445598307E-02_wp, &
       & -4.14055422698917E+00_wp, +4.94118711261405E+00_wp, +6.99329422648258E-01_wp, &
@@ -1056,7 +1062,7 @@ subroutine test_cn_mb02_erf_en(error)
       rcov(:) = get_covalent_rad(mol%num)
       en(:) = get_pauling_en(mol%num)
 
-      call new_erf_en_ncoord(erf_en_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en)
+      call new_erf_en_ncoord(erf_en_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov, en=en)
       call test_cn_gen(error, mol,  erf_en_ncoord, cutoff, ref)
 
    end subroutine test_cn_mb02_erf_en
@@ -1073,6 +1079,7 @@ subroutine test_cn_mb03_erf_en(error)
       real(wp), allocatable :: en(:)
 
       real(wp), parameter :: cutoff = 30.0_wp
+      real(wp), parameter :: kcn = 2.60_wp
       real(wp), parameter :: ref(16) = [&
       & -1.57129310680311E+00_wp, -5.12538648600132E+00_wp, +3.44404010675507E-02_wp, &
       & +5.53888057043349E+00_wp, +6.08403921871331E+00_wp, +1.70484401191915E+00_wp, &
@@ -1087,7 +1094,7 @@ subroutine test_cn_mb03_erf_en(error)
       rcov(:) = get_covalent_rad(mol%num)
       en(:) = get_pauling_en(mol%num)
 
-      call new_erf_en_ncoord(erf_en_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en)
+      call new_erf_en_ncoord(erf_en_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov, en=en)
       call test_cn_gen(error, mol, erf_en_ncoord, cutoff, ref)
 
    end subroutine test_cn_mb03_erf_en
@@ -1104,6 +1111,7 @@ subroutine test_cn_acetic_erf_en(error)
       real(wp), allocatable :: en(:)
 
       real(wp), parameter :: cutoff = 30.0_wp
+      real(wp), parameter :: kcn = 2.60_wp
       real(wp), parameter :: ref(32) = [&
       & +1.21457976613492E+00_wp, +1.21461935183543E+00_wp, +1.21453529477741E+00_wp, &
       & +1.21459408863060E+00_wp, +2.94977263784175E-01_wp, +2.94968468529279E-01_wp, &
@@ -1123,7 +1131,7 @@ subroutine test_cn_acetic_erf_en(error)
       rcov(:) = get_covalent_rad(mol%num)
       en(:) = get_pauling_en(mol%num)
 
-      call new_erf_en_ncoord(erf_en_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en)
+      call new_erf_en_ncoord(erf_en_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov, en=en)
       call test_cn_gen(error, mol, erf_en_ncoord, cutoff, ref)
 
    end subroutine test_cn_acetic_erf_en