Fix gradient evaluation when one of coordinates is close to another (#…

…1208) * Do not skip overlap integrals close to zero Signed-off-by: Igor S. Gerasimov <[email protected]> * Avoid division of zero by zero Signed-off-by: Igor S. Gerasimov <[email protected]> * Update tests Signed-off-by: Igor S. Gerasimov <[email protected]> * Remove unused matlist2 and nmat2 Signed-off-by: Igor S. Gerasimov <[email protected]> --------- Signed-off-by: Igor S. Gerasimov <[email protected]> Co-authored-by: Thomas Froitzheim <[email protected]>
grimme-lab · Mar 6, 2025 · fc512bb · fc512bb
1 parent 04e5303
commit fc512bb
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 29 deletions.
diff --git a/src/scf_module.F90 b/src/scf_module.F90
@@ -173,7 +173,6 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, &
    integer :: ndp,nqp
 
    integer,allocatable :: matlist (:,:)
-   integer,allocatable :: matlist2(:,:)
    integer,allocatable :: xblist(:,:)
    real(wp),allocatable :: sqrab(:)
    real(wp),allocatable :: dcndr(:,:,:)
@@ -198,7 +197,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, &
 
    integer :: ich ! file handle
    integer :: npr,ii,jj,kk,i,j,k,m,iat,jat,mi,jter,atj,kkk,mj,mm
-   integer :: ishell,jshell,np,ia,ndimv,l,nmat,nmat2
+   integer :: ishell,jshell,np,ia,ndimv,l,nmat
    integer :: ll,i1,i2,nn,ati,nxb,lin,startpdiag
    integer :: is,js,gav
 
@@ -326,8 +325,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, &
    &        X(basis%nao,basis%nao), &
    &        ves(basis%nshell), &
    &        zsh(basis%nshell),&
-   &        matlist (2,basis%nao*(basis%nao+1)/2), &
-   &        matlist2(2,basis%nao*(basis%nao+1)/2-basis%nao))
+   &        matlist (2,basis%nao*(basis%nao+1)/2))
    allocate(selfEnergy(maxval(xtbData%nshell), mol%n))
    allocate(dSEdcn(maxval(xtbData%nshell), mol%n))
 
@@ -587,24 +585,13 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, &
    ! ------------------------------------------------------------------------
    ! prepare matrix indices
    nmat =0
-   nmat2=0
    do ii=1,basis%nao
       iat=basis%aoat2(ii)
       do jj=1,ii-1
          jat=basis%aoat2(jj)
-         if(abs(S(jj,ii)).lt.neglect) then
-            S(jj,ii)=0.0_wp
-            S(ii,jj)=0.0_wp
-            cycle
-         endif
          nmat=nmat+1
          matlist(1,nmat)=ii
          matlist(2,nmat)=jj
-         if(iat.ne.jat)then
-            nmat2=nmat2+1
-            matlist2(1,nmat2)=ii
-            matlist2(2,nmat2)=jj
-         endif
       enddo
       ! CB: moved this here so j/i indices from matlist come in a reasonable order
       ! also setup CN dep. stuff
@@ -713,17 +700,21 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, &
    if (mol%npbc == 0) then
       allocate(H(basis%nao, basis%nao))
       H(:, :) = 0.0_wp
-      do m = 1, nmat2
-         i = matlist2(1,m)
-         j = matlist2(2,m)
-         k = j+i*(i-1)/2
-         !ishell = ao2sh(i)
-         !jshell = ao2sh(j)
-         ! SCC terms
-         !eh1 = autoev*(shellShift(ishell) + shellShift(jshell))
-         !H1 = -S(j,i)*eh1*0.5_wp
-         H(j,i) = H0(k)*evtoau/S(j,i)
-         H(i,j) = H(j,i)
+      do i = 1, basis%nao
+         do j = 1, i-1
+            k = j+i*(i-1)/2
+            !ishell = ao2sh(i)
+            !jshell = ao2sh(j)
+            ! SCC terms
+            !eh1 = autoev*(shellShift(ishell) + shellShift(jshell))
+            !H1 = -S(j,i)*eh1*0.5_wp
+            if (H0(k) .eq. 0.0_wp) then
+               H(j,i) = 0.0_wp
+            else
+               H(j,i) = H0(k)*evtoau/S(j,i)
+            end if
+            H(i,j) = H(j,i)
+         end do
       enddo
       call build_dSDQH0_noreset(xtbData%nShell, xtbData%hamiltonian, selfEnergy, &
          & dSEdcn, intcut, mol%n, basis%nao, basis%nbf, mol%at, mol%xyz, &

diff --git a/test/unit/test_gfn2.f90 b/test/unit/test_gfn2.f90
@@ -299,11 +299,11 @@ subroutine test_gfn2gbsa_api(error)
 
    call check_(error, hl_gap, 3.408607724814_wp, thr=1e-5_wp)
    call check_(error, energy,-22.002501380096_wp, thr=thr)
-   call check_(error, norm2(gradient),0.19441977481008E-01_wp, thr=thr)
+   call check_(error, norm2(gradient),0.19366866479022E-01_wp, thr=thr)
    call check_(error, gradient(1,1), .9038674439127e-02_wp, thr=thr)
    call check_(error, gradient(3,2),-.1394693523214e-02_wp, thr=thr)
    call check_(error, gradient(3,11),-gradient(3,10), thr=thr)
-   call check_(error, gradient(1,8),0.22890674680144E-02_wp, thr=thr)
+   call check_(error, gradient(1,8),0.17878350605942E-02_wp, thr=thr)
 
 end subroutine test_gfn2gbsa_api
 
@@ -1080,7 +1080,7 @@ subroutine test_gfn2_wbo(error)
    ! check  scc !
    call check_(error, energy, -33.314958144107_wp, thr=thr)
    call check_(error, norm2(gradient), 0.018945181484_wp, thr=thr)
-   call check_(error, hl_gap, 1.805948277212_wp, thr=thr)
+   call check_(error, hl_gap, 1.805948321662_wp, thr=thr)
 
 
    ! check wbo !