Pass angular moment instead array offset to integral evaluator

grimme-lab · Jan 25, 2021 · 0bc838a · 0bc838a
1 parent c173684
commit 0bc838a
Show file tree

Hide file tree

Showing 4 changed files with 41 additions and 33 deletions.
diff --git a/src/intgrad.f90 b/src/intgrad.f90
@@ -871,14 +871,14 @@ pure subroutine build_ds_ints(a,b,e,alpi,alpj,la,lb,t,v,g)
 
 end subroutine build_ds_ints
 
-pure subroutine get_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,intcut, &
+pure subroutine get_overlap(icao,jcao,naoi,naoj,ishtyp,jshtyp,ri,rj,point,intcut, &
       &                nprim,primcount,alp,cont,sint)
    integer, intent(in)  :: icao
    integer, intent(in)  :: jcao
    integer, intent(in)  :: naoi
    integer, intent(in)  :: naoj
-   integer, intent(in)  :: iptyp
-   integer, intent(in)  :: jptyp
+   integer, intent(in)  :: ishtyp
+   integer, intent(in)  :: jshtyp
    real(wp),intent(in)  :: ri(3)
    real(wp),intent(in)  :: rj(3)
    real(wp),intent(in)  :: point(3)
@@ -890,13 +890,15 @@ pure subroutine get_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,intcut,
    real(wp),intent(in)  :: alp(:)
    real(wp),intent(in)  :: cont(:)
 
-   integer  :: ip,iprim,mli,jp,jprim,mlj,k
+   integer  :: ip,iprim,mli,jp,jprim,mlj,k,iptyp,jptyp
    real(wp) :: rij(3),rij2,alpi,alpj,ci,cj,cc,kab,rp(3),t(0:8)
    real(wp) :: ab,est,saw(10)
 
    real(wp),parameter :: max_r2 = 2000.0_wp
 
    sint = 0.0_wp
+   iptyp = itt(ishtyp)
+   jptyp = itt(jshtyp)
 
    rij = ri - rj
    rij2 = rij(1)**2 + rij(2)**2 + rij(3)**2
@@ -916,7 +918,7 @@ pure subroutine get_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,intcut,
          if(est.gt.intcut) cycle
          call build_kab(ri,alpi,rj,alpj,ab,kab)
          rp = gpcenter(alpi,ri,alpj,rj)
-         do k = 0, 8 ! maxval(ij)
+         do k = 0, ishtyp + jshtyp
             t(k) = olapp(k, alpi+alpj)
          end do
          ! now compute integrals
@@ -938,14 +940,14 @@ pure subroutine get_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,intcut,
 
 end subroutine get_overlap
 
-pure subroutine get_grad_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,intcut, &
+pure subroutine get_grad_overlap(icao,jcao,naoi,naoj,ishtyp,jshtyp,ri,rj,point,intcut, &
       &                     nprim,primcount,alp,cont,sdq,sdqg)
    integer, intent(in)  :: icao
    integer, intent(in)  :: jcao
    integer, intent(in)  :: naoi
    integer, intent(in)  :: naoj
-   integer, intent(in)  :: iptyp
-   integer, intent(in)  :: jptyp
+   integer, intent(in)  :: ishtyp
+   integer, intent(in)  :: jshtyp
    real(wp),intent(in)  :: ri(3)
    real(wp),intent(in)  :: rj(3)
    real(wp),intent(in)  :: point(3)
@@ -958,14 +960,16 @@ pure subroutine get_grad_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,int
    real(wp),intent(in)  :: alp(:)
    real(wp),intent(in)  :: cont(:)
 
-   integer  :: ip,iprim,mli,jp,jprim,mlj,k
+   integer  :: ip,iprim,mli,jp,jprim,mlj,k,iptyp,jptyp
    real(wp) :: rij(3),rij2,alpi,alpj,ci,cj,cc,kab,rp(3),t(0:8)
    real(wp) :: ab,est,saw,sawg(3)
 
    real(wp),parameter :: max_r2 = 2000.0_wp
 
    sdqg = 0.0_wp
    sdq  = 0.0_wp
+   iptyp = itt(ishtyp)
+   jptyp = itt(jshtyp)
 
    rij = ri - rj
    rij2 = rij(1)**2 + rij(2)**2 + rij(3)**2
@@ -985,7 +989,7 @@ pure subroutine get_grad_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,int
          if(est.gt.intcut) cycle
          call build_kab(ri,alpi,rj,alpj,ab,kab)
          rp = gpcenter(alpi,ri,alpj,rj)
-         do k = 0, 8 ! maxval(ij)+1
+         do k = 0, ishtyp + jshtyp + 1
             t(k) = olapp(k, alpi+alpj)
          end do
          !--------------- compute gradient ----------
@@ -1008,14 +1012,14 @@ pure subroutine get_grad_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,int
    enddo  ! ip : loop over i prims
 end subroutine get_grad_overlap
 
-pure subroutine get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,intcut, &
+pure subroutine get_multiints(icao,jcao,naoi,naoj,ishtyp,jshtyp,ri,rj,point,intcut, &
       &                       nprim,primcount,alp,cont,ss,dd,qq)
    integer, intent(in)  :: icao
    integer, intent(in)  :: jcao
    integer, intent(in)  :: naoi
    integer, intent(in)  :: naoj
-   integer, intent(in)  :: iptyp
-   integer, intent(in)  :: jptyp
+   integer, intent(in)  :: ishtyp
+   integer, intent(in)  :: jshtyp
    real(wp),intent(in)  :: ri(3)
    real(wp),intent(in)  :: rj(3)
    real(wp),intent(in)  :: point(3)
@@ -1029,7 +1033,7 @@ pure subroutine get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,intcut
    real(wp),intent(in)  :: alp(:)
    real(wp),intent(in)  :: cont(:)
 
-   integer  :: ip,iprim,mli,jp,jprim,mlj,k
+   integer  :: ip,iprim,mli,jp,jprim,mlj,k,iptyp,jptyp
    real(wp) :: rij(3),rij2,alpi,alpj,ci,cj,cc,kab,rp(3),t(0:8)
    real(wp) :: ab,est,saw(10)
 
@@ -1038,6 +1042,8 @@ pure subroutine get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,intcut
    ss = 0.0_wp
    dd = 0.0_wp
    qq = 0.0_wp
+   iptyp = itt(ishtyp)
+   jptyp = itt(jshtyp)
 
    rij = rj - rj
    rij2 = rij(1)**2 + rij(2)**2 + rij(3)**2
@@ -1055,7 +1061,7 @@ pure subroutine get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,intcut
          if(est.gt.intcut) cycle
          call build_kab(ri,alpi,rj,alpj,ab,kab)
          rp = gpcenter(alpi,ri,alpj,rj)
-         do k = 0, 8 ! maxval(ij) + 2
+         do k = 0, ishtyp + jshtyp + 2
             t(k) = olapp(k, alpi+alpj)
          end do
          ! now compute integrals  for different components of i(e.g., px,py,pz)
@@ -1084,14 +1090,14 @@ pure subroutine get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point,intcut
 end subroutine get_multiints
 
 
-pure subroutine get_grad_multiint(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj, &
+pure subroutine get_grad_multiint(icao,jcao,naoi,naoj,ishtyp,jshtyp,ri,rj, &
       &                           intcut,nprim,primcount,alp,cont,sdq,sdqg)
    integer, intent(in)  :: icao
    integer, intent(in)  :: jcao
    integer, intent(in)  :: naoi
    integer, intent(in)  :: naoj
-   integer, intent(in)  :: iptyp
-   integer, intent(in)  :: jptyp
+   integer, intent(in)  :: ishtyp
+   integer, intent(in)  :: jshtyp
    real(wp),intent(in)  :: ri(3)
    real(wp),intent(in)  :: rj(3)
    real(wp),intent(in)  :: intcut
@@ -1103,14 +1109,16 @@ pure subroutine get_grad_multiint(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj, &
    real(wp),intent(in)  :: alp(:)
    real(wp),intent(in)  :: cont(:)
 
-   integer  :: ip,iprim,mli,jp,jprim,mlj,k
+   integer  :: ip,iprim,mli,jp,jprim,mlj,k,iptyp,jptyp
    real(wp) :: rij(3),rp(3),rij2,alpi,alpj,ci,cj,cc,kab,t(0:8)
    real(wp) :: ab,est,saw(10),sawg(3,10)
 
    real(wp),parameter :: max_r2 = 2000.0_wp
 
    sdqg = 0.0_wp
    sdq  = 0.0_wp
+   iptyp = itt(ishtyp)
+   jptyp = itt(jshtyp)
 
    rij = ri - rj
    rij2 = rij(1)**2 + rij(2)**2 + rij(3)**2
@@ -1131,7 +1139,7 @@ pure subroutine get_grad_multiint(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj, &
          if(est.gt.intcut) cycle
          call build_kab(ri,alpi,rj,alpj,ab,kab)
          rp = gpcenter(alpi,ri,alpj,rj)
-         do k = 0, 8 ! maxval(ij) + 3
+         do k = 0, ishtyp + jshtyp + 3
             t(k) = olapp(k, alpi+alpj)
          end do
          !--------------- compute gradient ----------
@@ -1245,7 +1253,7 @@ subroutine sdqint(nShell,angShell,nat,at,nbf,nao,xyz,intcut,caoshell,saoshell, &
                ss = 0.0_wp
                dd = 0.0_wp
                qq = 0.0_wp
-               call get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ra,rb,point, &
+               call get_multiints(icao,jcao,naoi,naoj,ishtyp,jshtyp,ra,rb,point, &
                   &               intcut,nprim,primcount,alp,cont,ss,dd,qq)
                !transform from CAO to SAO
                call dtrf2(ss,ishtyp,jshtyp)
@@ -1298,7 +1306,7 @@ subroutine sdqint(nShell,angShell,nat,at,nbf,nao,xyz,intcut,caoshell,saoshell, &
             ss = 0.0_wp
             dd = 0.0_wp
             qq = 0.0_wp
-            call get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ra,ra,point, &
+            call get_multiints(icao,jcao,naoi,naoj,ishtyp,jshtyp,ra,ra,point, &
                &               intcut,nprim,primcount,alp,cont,ss,dd,qq)
             !transform from CAO to SAO
             !call dtrf2(ss,ishtyp,jshtyp)
@@ -1493,7 +1501,7 @@ subroutine sdqint_gpu(nShell, angShell, nat, at, nbf, nao, xyz, trans, &
                   dd = 0.0_wp
                   qq = 0.0_wp
 
-                  !call get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ra,rb,point, &
+                  !call get_multiints(icao,jcao,naoi,naoj,ishtyp,jshtyp,ra,rb,point, &
                   !   &               intcut,nprim,primcount,alp,cont,ss,dd,qq)
                   !$acc loop vector private(saw,t,e) independent collapse(2)
                   do ip = 1,nprim(icao+1)
@@ -1878,7 +1886,7 @@ subroutine sdqint_gpu(nShell, angShell, nat, at, nbf, nao, xyz, trans, &
             ss = 0.0_wp
             dd2 = 0.0_wp
             qq = 0.0_wp
-            call get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ra,ra,point, &
+            call get_multiints(icao,jcao,naoi,naoj,ishtyp,jshtyp,ra,ra,point, &
                &               intcut,nprim,primcount,alp,cont,ss,dd2,qq)
             !transform from CAO to SAO
             !call dtrf2(ss,ishtyp,jshtyp)

diff --git a/src/peeq_module.f90 b/src/peeq_module.f90
@@ -944,7 +944,7 @@ subroutine ccm_build_SH0(nShell, hData, selfEnergy, nat, at, basis, nbf, nao, &
                      &               hData%atomicRad(ati),hData%atomicRad(atj),ri,rj)
 
                   ! get overlap integral
-                  call get_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point, &
+                  call get_overlap(icao,jcao,naoi,naoj,ishtyp,jshtyp,ri,rj,point, &
                      &             intcut,basis%nprim,basis%primcount, &
                      &             basis%alp,basis%cont,ss)
 
@@ -1137,7 +1137,7 @@ subroutine pbc_build_SH0(nShell, hData, selfEnergy, nat, at, basis, nbf, nao, &
                      &               hData%atomicRad(ati),hData%atomicRad(atj),ri,rj)
 
                   ! get overlap integral
-                  call get_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj,point, &
+                  call get_overlap(icao,jcao,naoi,naoj,ishtyp,jshtyp,ri,rj,point, &
                      &             intcut,basis%nprim,basis%primcount, &
                      &             basis%alp,basis%cont,ss)
 
@@ -1334,7 +1334,7 @@ subroutine ccm_build_dSH0(nShell, hData, selfEnergy, dSEdcn, dSEdq, nat, basis,
                      & hData%atomicRad(ati),hData%atomicRad(atj),rij2,ri,rj,&
                      & shpoly,dshpoly)
 
-                  call get_grad_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj, &
+                  call get_grad_overlap(icao,jcao,naoi,naoj,ishtyp,jshtyp,ri,rj, &
                      &                  point,thr,basis%nprim,basis%primcount, &
                      &                  basis%alp,basis%cont,sdq,sdqg)
 
@@ -1573,7 +1573,7 @@ subroutine pbc_build_dSH0(nShell, hData, selfEnergy, dSEdcn, dSEdq, nat, basis,
                      & hData%atomicRad(ati),hData%atomicRad(atj),rij2,ri,rj,&
                      & shpoly,dshpoly)
 
-                  call get_grad_overlap(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj, &
+                  call get_grad_overlap(icao,jcao,naoi,naoj,ishtyp,jshtyp,ri,rj, &
                      &                  point,thr,basis%nprim,basis%primcount, &
                      &                  basis%alp,basis%cont,sdq,sdqg)
 

diff --git a/src/xtb/hamiltonian.f90 b/src/xtb/hamiltonian.f90
@@ -256,7 +256,7 @@ subroutine build_SDQH0(nShell, hData, nat, at, nbf, nao, xyz, trans, selfEnergy,
                   ss = 0.0_wp
                   dd = 0.0_wp
                   qq = 0.0_wp
-                  call get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ra,rb,point, &
+                  call get_multiints(icao,jcao,naoi,naoj,ishtyp,jshtyp,ra,rb,point, &
                      &               intcut,nprim,primcount,alp,cont,ss,dd,qq)
                   !transform from CAO to SAO
                   call dtrf2(ss,ishtyp,jshtyp)
@@ -320,7 +320,7 @@ subroutine build_SDQH0(nShell, hData, nat, at, nbf, nao, xyz, trans, selfEnergy,
             ss = 0.0_wp
             dd = 0.0_wp
             qq = 0.0_wp
-            call get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ra,ra,point, &
+            call get_multiints(icao,jcao,naoi,naoj,ishtyp,jshtyp,ra,ra,point, &
                &               intcut,nprim,primcount,alp,cont,ss,dd,qq)
             !transform from CAO to SAO
             !call dtrf2(ss,ishtyp,jshtyp)
@@ -476,7 +476,7 @@ subroutine build_dSDQH0(nShell, hData, selfEnergy, dSEdcn, intcut, nat, nao, nbf
                      & shpoly,dshpoly)
 
                   sdqg = 0;sdq = 0
-                  call get_grad_multiint(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj, &
+                  call get_grad_multiint(icao,jcao,naoi,naoj,ishtyp,jshtyp,ri,rj, &
                      &                   intcut,nprim,primcount,alp,cont,sdq,sdqg)
                   tmp(1:6,1:6) = sdq(1,1:6,1:6)
                   call dtrf2(tmp,ishtyp,jshtyp)
@@ -668,7 +668,7 @@ subroutine build_dSDQH0_noreset(nShell, hData, selfEnergy, dSEdcn, intcut, &
                hav = 0.5_wp * (hii + hjj)
 
                sdqg = 0;sdq = 0
-               call get_grad_multiint(icao,jcao,naoi,naoj,iptyp,jptyp,ri,rj, &
+               call get_grad_multiint(icao,jcao,naoi,naoj,ishtyp,jshtyp,ri,rj, &
                   &                   intcut,nprim,primcount,alp,cont,sdq,sdqg)
                do k = 1,19 ! 1 S, 2-4 D, 5-10 Q, 11-13 D, 14-19 Q
                   do ixyz = 1,3

diff --git a/src/xtb/hamiltonian_gpu.f90 b/src/xtb/hamiltonian_gpu.f90
@@ -722,7 +722,7 @@ subroutine build_SDQH0_gpu(nShell, hData, nat, at, nbf, nao, xyz, trans, selfEne
             ss = 0.0_wp
             dd2 = 0.0_wp
             qq = 0.0_wp
-            call get_multiints(icao,jcao,naoi,naoj,iptyp,jptyp,ra,ra,point, &
+            call get_multiints(icao,jcao,naoi,naoj,ishtyp,jshtyp,ra,ra,point, &
                &               intcut,nprim,primcount,alp,cont,ss,dd2,qq)
             !transform from CAO to SAO
             !call dtrf2(ss,ishtyp,jshtyp)