refactored the Wigner-Seitz cell representation (#16)

- does not include molecule specific information anymore
- only topological information in WSC (setup only required once)
- now possible for arbitrary periodicities
- tests for 0D and 3D Wigner-Seitz cell setup

TESTSUITE/tbdef_wsc.f90

@@ -1,11 +1,90 @@
-subroutine test_wigner_seitz
+subroutine test_wigner_seitz_0d
+   use iso_fortran_env, wp => real64
+   use assertion
+   use tbdef_molecule
+   use tbdef_wsc
+   implicit none
+   type(tb_molecule) :: mol
+   type(tb_wsc) :: wsc
 
-   ipar = 0
-   do i = 1, mol%n
-      do j = 1, i
-         ipar = ipar + mol%wsc%itbl(j,i)
-      enddo
-   enddo
-   call assert_eq(ipar,mol%wsl%np)
+   real(wp),parameter :: thr = 1.0e-10_wp
+   integer, parameter :: nat = 3
+   integer, parameter :: at(nat) = [8,1,1]
+   real(wp),parameter :: xyz(3,nat) = reshape(&
+      [ 0.00000000000000_wp,  0.00000000000000_wp, -0.73578586109551_wp, &
+      & 1.44183152868459_wp,  0.00000000000000_wp,  0.36789293054775_wp, &
+      &-1.44183152868459_wp,  0.00000000000000_wp,  0.36789293054775_wp  &
+      & ],shape(xyz))
 
-end subroutine test_wigner_seitz
+   call mol%allocate(nat)
+   mol%at  = at
+   mol%xyz = xyz
+   mol%chrg = 0.0_wp
+   call mol%update
+
+   call generate_wsc(mol,wsc)
+
+   call assert(allocated(wsc%lattr))
+   call assert(allocated(wsc%itbl))
+
+   call assert(all(wsc%itbl <= 1))
+   call assert(all(wsc%lattr == 0))
+
+   call terminate(afail)
+end subroutine test_wigner_seitz_0D
+
+subroutine test_wigner_seitz_3D
+   use iso_fortran_env, wp => real64, istdout => output_unit
+   use assertion
+
+   use tbdef_molecule
+   use tbdef_wsc
+
+   use pbc_tools
+
+   implicit none
+
+   real(wp),parameter :: thr = 1.0e-10_wp
+   ! CaF2
+   integer, parameter :: nat = 3
+   integer, parameter :: at(nat) = [9,9,20]
+   real(wp),parameter :: abc(3,nat) = reshape(&
+      &[0.25_wp, 0.25_wp, 0.25_wp, &
+      & 0.75_wp, 0.75_wp, 0.75_wp, &
+      & 0.00_wp, 0.00_wp, 0.00_wp], shape(abc))
+   real(wp),parameter :: lattice(3,3) = reshape(&
+      &[5.9598811567890_wp,      2.1071361905157_wp,      3.6496669404404_wp,    &
+      & 0.0000000000000_wp,      6.3214085715472_wp,      3.6496669404404_wp,    &
+      & 0.0000000000000_wp,      0.0000000000000_wp,      7.2993338808807_wp],   &
+      & shape(lattice))
+
+   type(tb_molecule) :: mol
+   type(tb_wsc) :: wsc
+
+   call mol%allocate(nat)
+   mol%at   = at
+   mol%lattice = lattice
+   call coord_trafo(nat,lattice,abc,mol%xyz)
+   mol%npbc = 3
+   mol%pbc  = .true.
+   call mol%update
+
+   call generate_wsc(mol,wsc)
+
+   call assert(allocated(wsc%lattr))
+   call assert(allocated(wsc%itbl))
+
+   print*,wsc%itbl
+   call assert(all(abs(wsc%lattr) <= 1))
+   call assert_eq(wsc%itbl(1,1), 12)
+   call assert_eq(wsc%itbl(2,2), 12)
+   call assert_eq(wsc%itbl(3,3), 12)
+   call assert_eq(wsc%itbl(1,2), 6)
+   call assert_eq(wsc%itbl(1,3), 4)
+   call assert_eq(wsc%itbl(2,3), 4)
+   call assert_eq(wsc%itbl(1,2), wsc%itbl(2,1))
+   call assert_eq(wsc%itbl(1,3), wsc%itbl(3,1))
+   call assert_eq(wsc%itbl(2,3), wsc%itbl(3,2))
+
+   call terminate(afail)
+end subroutine test_wigner_seitz_3D

TESTSUITE/tests_peeq.f90

@@ -102,6 +102,11 @@ program peeq_tester
       case('mic');  call test_class_molecule_mic_distances
       case('axis'); call test_class_molecule_axis_trafo
       end select
+   case('tbdef_wsc')
+      select case(sec)
+      case('0d'); call test_wigner_seitz_0d
+      case('3d'); call test_wigner_seitz_3d
+      end select
    case('symmetry')
       select case(sec)
       case('water');  call test_symmetry_water

meson.build

@@ -320,6 +320,7 @@ xtb_srcs += 'xtb/c_api.f90'
 xtb_test = ['TESTSUITE/tests_peeq.f90']
 xtb_test += 'TESTSUITE/assertion.f90'
 xtb_test += 'TESTSUITE/tbdef_molecule.f90'
+xtb_test += 'TESTSUITE/tbdef_wsc.f90'
 xtb_test += 'TESTSUITE/geometry_reader.f90'
 xtb_test += 'TESTSUITE/eeq_model.f90'
 xtb_test += 'TESTSUITE/pbc_tools.f90'
@@ -414,6 +415,9 @@ test('Method: GFN2-xTB/GBSA'            ,
 test('Molecule Class: axis',xtb_test,args: ['tbdef_molecule','axis'])
 test('Molecule Class: MIC', xtb_test,args: ['tbdef_molecule','mic'])
 
+test('Wigner-Seitz Cell (0D)', xtb_test,args: ['tbdef_wsc','0d'])
+test('Wigner-Seitz Cell (3D)', xtb_test,args: ['tbdef_wsc','3d'])
+
 test('Geometry Reader: coord 3D',xtb_test,args: ['geometry_reader','coord_3d_a'])
 test('Geometry Reader: coord 3D',xtb_test,args: ['geometry_reader','coord_3d_b'])
 test('Geometry Reader: coord 2D',xtb_test,args: ['geometry_reader','coord_2d'])

xtb/eeq_model.f90

@@ -2073,7 +2073,8 @@ subroutine eeq_chrgeq_core(mol,chrgeq,cn,dcndr,dcndL,q,dqdr,dqdL, &
       do j = 1, i-1
          gamij = 1.0_wp/sqrt(alpha2(i)+alpha2(j))
          do wscAt = 1, mol%wsc%itbl(j,i)
-            riw   = mol%xyz(:,i) - mol%wsc%xyz(:,wscAt,j,i)
+            riw = mol%xyz(:,i) - mol%xyz(:,j) &
+               &  - matmul(mol%lattice,mol%wsc%lattr(:,wscAt,j,i))
             Amat(i,j) = Amat(i,j) + mol%wsc%w(j,i) * ( &
                ! reciprocal lattice sums
                + eeq_ewald_3d_rec(riw,ewaldCutR,mol%rec_lat,mol%volume,cf) &
@@ -2187,7 +2188,8 @@ subroutine eeq_chrgeq_core(mol,chrgeq,cn,dcndr,dcndL,q,dqdr,dqdL, &
          ! over WSC partner
          gamij = 1.0_wp/sqrt(alpha2(i) + alpha2(j))
          do wscAt = 1, mol%wsc%itbl(j,i)
-            riw = mol%xyz(:,i) - mol%wsc%xyz(:,wscAt,j,i)
+            riw = mol%xyz(:,i) - mol%xyz(:,j) &
+               &  - matmul(mol%lattice,mol%wsc%lattr(:,wscAt,j,i))
             call eeq_ewald_dx_3d_rec(riw,ewaldCutR,mol%rec_lat,mol%volume,cf, &
                &                     dAtmp,stmp)
             dAtmp = dAtmp*mol%wsc%w(j,i)

xtb/generate_wsc.f90

@@ -16,7 +16,7 @@
 ! along with xtb.  If not, see <https://www.gnu.org/licenses/>.
 
 !> generate a Wigner--Seitz cell from a given structure
-subroutine generate_wsc(mol,wsc,rep)
+subroutine generate_wsc(mol,wsc)
    use iso_fortran_env, wp => real64
    use tbdef_molecule
    use tbdef_wsc
@@ -25,43 +25,48 @@ subroutine generate_wsc(mol,wsc,rep)
    type(tb_molecule),intent(inout) :: mol
    !> Wigner--Seitz cell data type (might be contained in mol)
    type(tb_wsc),     intent(inout) :: wsc
-   !> images of the unit cell to consider
-   integer,          intent(in)    :: rep(3)
 ! ------------------------------------------------------------------------
 !  Variables
 ! ------------------------------------------------------------------------
+   integer  :: rep(3)
    integer  :: ii,jj,ich
    integer  :: aa,bb,cc
    integer  :: c,wc
    integer  :: minpos
    integer  :: nminpos
-   real(wp) :: t(3)
+   real(wp) :: t(3),rw(3)
    real(wp) :: mindist
    real(wp) :: nmindist
    !> overall WSC tolerance to consider atoms as WSC-images
    real(wp),parameter :: tol = 0.01_wp
-   real(wp),allocatable,dimension(:,:,:) :: txyz
+   integer, allocatable,dimension(:,:)   :: lattr
    real(wp),allocatable,dimension(:)     :: dist
    logical, allocatable,dimension(:)     :: trans
 
+   where(mol%pbc)
+      rep = 1
+   elsewhere
+      rep = 0
+   endwhere
+
 ! ------------------------------------------------------------------------
 !  allocate space for the WSC first
-   call wsc%allocate(mol%n,rep,mol%lattice)
+   call wsc%allocate(mol%n,rep)
 
 ! ------------------------------------------------------------------------
 ! initialize
 ! ------------------------------------------------------------------------
-   allocate( txyz(3,wsc%cells,mol%n) ); txyz = 0.0_wp
-   allocate( dist(wsc%cells) );     dist = 0.0_wp
-   allocate( trans(wsc%cells) );    trans = .true.
+   allocate( lattr(3,wsc%cells),      source = 0 )
+   allocate( dist(wsc%cells),         source = 0.0_wp )
+   allocate( trans(wsc%cells),        source = .true. )
 ! ------------------------------------------------------------------------
 ! Create the Wigner-Seitz Cell (WSC)
 ! ------------------------------------------------------------------------
    wsc%at  = 0
    wsc%itbl= 0
 !$omp parallel default(none) &
-!$omp private(ii,jj,wc,c,dist,trans,t) &
-!$omp shared(mol,wsc,rep,txyz) &
+!$omp private(ii,jj,wc,c,dist,trans,t,lattr,rw) &
+!$omp shared(mol,wsc,rep) &
 !$omp shared(mindist,minpos,nmindist,nminpos)
 !$omp do schedule(dynamic)
    ! Each WSC of one atom consists of n atoms
@@ -71,14 +76,16 @@ subroutine generate_wsc(mol,wsc,rep)
          ! find according neighbours
          c=0
          dist = 0.0_wp
+         lattr = 0
          do aa=-rep(1),rep(1),1
             do bb=-rep(2),rep(2),1
                do cc=-rep(3),rep(3),1
                   if ((aa.eq.0 .and. bb.eq.0 .and. cc.eq.0).and.ii.eq.jj) cycle
                   t = [aa,bb,cc]
                   c=c+1
-                  txyz(:,c,ii) = mol%xyz(:,jj) + matmul(mol%lattice,t)
-                  dist(c)=norm2(mol%xyz(:,ii)-txyz(:,c,ii))
+                  lattr(:,c) = [aa,bb,cc]
+                  rw = mol%xyz(:,jj) + matmul(mol%lattice,t)
+                  dist(c)=norm2(mol%xyz(:,ii)-rw)
                end do
             end do
          end do
@@ -90,22 +97,28 @@ subroutine generate_wsc(mol,wsc,rep)
          mindist=dist(minpos)
          trans(minpos)=.false.
          wc=1
-         wsc%xyz(:,wc,jj,ii)=txyz(:,minpos,ii)
+         wsc%lattr(:,wc,jj,ii)=lattr(:,minpos)
          ! get other images with same distance
-         find_images : do
-            nminpos=minloc(dist(:c),dim=1,mask=trans(:c))
-            nmindist=dist(nminpos)
-            if(abs(mindist-nmindist).lt.tol)then
-               trans(nminpos)=.false.
-               wc=wc+1
-               wsc%xyz(:,wc,jj,ii)=txyz(:,nminpos,ii)
-            else
-               wsc%w(jj,ii)=1.0_wp/real(wc,wp)
-               wsc%itbl(jj,ii)=wc
-               wsc%at(jj,ii)=jj
-               exit find_images
-            end if
-         end do find_images
+         if (c > 1) then
+            find_images : do
+               nminpos=minloc(dist(:c),dim=1,mask=trans(:c))
+               nmindist=dist(nminpos)
+               if(abs(mindist-nmindist).lt.tol)then
+                  trans(nminpos)=.false.
+                  wc=wc+1
+                  wsc%lattr(:,wc,jj,ii)=lattr(:,nminpos)
+               else
+                  wsc%w(jj,ii)=1.0_wp/real(wc,wp)
+                  wsc%itbl(jj,ii)=wc
+                  wsc%at(jj,ii)=jj
+                  exit find_images
+               end if
+            end do find_images
+         else
+            wsc%w(jj,ii) = 1.0_wp
+            wsc%itbl(jj,ii) = 1
+            wsc%at(jj,ii)=jj
+         endif
       end do
    end do
 !$omp enddo

xtb/peeq_module.f90

@@ -423,7 +423,8 @@ subroutine peeq &
 ! ---------------------------------------
    if (mol%npbc > 0) then
       if (ccm) then
-         call ccm_build_SH0(mol%n,mol%at,basis,nbf,nao,mol%xyz,wfn%q,cn,intcut, &
+         call ccm_build_SH0(mol%n,mol%at,basis,nbf,nao,mol%xyz,mol%lattice, &
+            &               wfn%q,cn,intcut, &
             &               param%kmagic,ken,param%alphaj,param%kcnsh, &
             &               param%xbdamp,s,h0,mol%wsc)
       else
@@ -522,7 +523,8 @@ subroutine peeq &
    call dmat(nao,tmp,wfn%C,pew)
    if (mol%npbc > 0) then
       if (ccm) then
-         call ccm_build_dSH0(mol%n,basis,intcut,nao,nbf,mol%at,mol%xyz,wfn%q,cn, &
+         call ccm_build_dSH0(mol%n,basis,intcut,nao,nbf,mol%at,mol%xyz, &
+            &                mol%lattice,wfn%q,cn, &
             &                wfn%P,Pew,g,sigma,dhdcn,dhdq,param%kmagic,ken, &
             &                param%alphaj,param%kcnsh,param%xbdamp,mol%wsc)
       else
@@ -1253,7 +1255,7 @@ end subroutine mol_build_SH0
 ! ------------------------------------------------------------------------
 !  Calculate the periodic AO overlap matrix
 ! ------------------------------------------------------------------------
-subroutine ccm_build_SH0(nat,at,basis,nbf,nao,xyz,q,cn,intcut, &
+subroutine ccm_build_SH0(nat,at,basis,nbf,nao,xyz,lattice,q,cn,intcut, &
       &                  kmagic,ken,alphaj,kcn,xbdamp,sint,h0,wsc)
 
    use tbdef_basisset
@@ -1275,6 +1277,7 @@ subroutine ccm_build_SH0(nat,at,basis,nbf,nao,xyz,q,cn,intcut, &
    integer, intent(in)  :: nao
    integer, intent(in)  :: at(nat)
    real(wp),intent(in)  :: xyz(3,nat)
+   real(wp),intent(in)  :: lattice(3,3)
    real(wp),intent(in)  :: q(nat)
    real(wp),intent(in)  :: cn(nat)
    real(wp),intent(in)  :: intcut
@@ -1302,7 +1305,7 @@ subroutine ccm_build_SH0(nat,at,basis,nbf,nao,xyz,q,cn,intcut, &
    real(wp) :: shpoly,km
    logical  :: valaoi,valaoj
 
-   real(wp)  ri(3),rj(3),f1,f2
+   real(wp)  ri(3),rj(3),t(3),f1,f2
    real(wp),parameter ::point(3) = 0.0_wp
    real(wp) dtmp(3),qtmp(6),ss(6,6),dd(3,6,6),qq(6,6,6),tmp(6,6)
    real(wp) sstmp(6,6) !TESTST
@@ -1323,9 +1326,9 @@ subroutine ccm_build_SH0(nat,at,basis,nbf,nao,xyz,q,cn,intcut, &
    !$omp&        ri,rj,icao,naoi,iptyp,jsh,jshmax,jshtyp,jcao,naoj,jptyp, &
    !$omp&        ss,saw,est,alpi,alpj,ab,iprim,jprim,ip,jp,km,shpoly, &
    !$omp&        mli,mlj,tmp,zi,zj,zetaij,enpoly,iao,jao, &
-   !$omp&        ii,jj,k,den,den2,den4,i,j,il,jl,hii,hjj,hav) &
+   !$omp&        ii,jj,k,den,den2,den4,i,j,il,jl,hii,hjj,hav,t) &
    !$omp reduction (+:sint,h0) &
-   !$omp shared(wsc,basis,at,ao_n,ao_l,xyz,intcut,nat,kqat,kcnat,cn,q,en, &
+   !$omp shared(wsc,basis,at,ao_n,ao_l,xyz,lattice,intcut,nat,kqat,kcnat,cn,q,en, &
    !$omp        ken,gam3,xbdamp,kcn,kmagic,kpair,alphaj)
    !$omp do schedule(runtime)
    do iat = 1, nat
@@ -1394,7 +1397,8 @@ subroutine ccm_build_SH0(nat,at,basis,nbf,nao,xyz,q,cn,intcut, &
 
                wscAt: do kat = 1, wsc%itbl(jat,iat)
                   ss = 0.0_wp
-                  rj = wsc%xyz(:,kat,jat,iat)
+                  t = wsc%lattr(:,kat,jat,iat)
+                  rj = xyz(:,jat) + matmul(lattice,t)
 
                   ! distance dependent polynomial
                   shpoly = rfactor(il,jl,ati,atj,ri,rj)
@@ -2018,7 +2022,7 @@ end subroutine mol_build_dSH0
 ! ------------------------------------------------------------------------
 !  Calculate the gradient resulting from a periodic AO overlap matrix
 ! ------------------------------------------------------------------------
-subroutine ccm_build_dSH0(nat,basis,thr,nao,nbf,at,xyz,q,cn,P,Pew,g,sigma, &
+subroutine ccm_build_dSH0(nat,basis,thr,nao,nbf,at,xyz,lattice,q,cn,P,Pew,g,sigma,&
       &                   dHdcn,dHdq,kmagic,ken,alphaj,kcn,xbdamp,wsc)
    use mctc_constants, only : pi
    use mctc_econv
@@ -2042,6 +2046,7 @@ subroutine ccm_build_dSH0(nat,basis,thr,nao,nbf,at,xyz,q,cn,P,Pew,g,sigma, &
    integer, intent(in)      :: nbf
    integer, intent(in) :: at(nat)
    real(wp),intent(in) :: xyz(3,nat)
+   real(wp),intent(in) :: lattice(3,3)
    real(wp),intent(in) :: q(nat)
    real(wp),intent(in) :: cn(nat)
    real(wp),intent(in) :: P(nao,nao)
@@ -2066,7 +2071,7 @@ subroutine ccm_build_dSH0(nat,basis,thr,nao,nbf,at,xyz,q,cn,P,Pew,g,sigma, &
    real(wp) tmp1,tmp2,tmp3,step,step2,step3,s00r,s00l,s00,alpj
    real(wp) skj,r1,r2,tt,t1,t2,t3,t4
    real(wp) thr2,f,ci,cc,cj,alpi,rab2,ab,est
-   real(wp) f1,f2,tmp(6,6),rij(3),ri(3),rj(3),rij2
+   real(wp) f1,f2,tmp(6,6),rij(3),ri(3),rj(3),rij2,t(3)
    real(wp) stmp,ral(3,3),rar(3,3),rbl(3,3),pre
    real(wp) dtmp,qtmp,rbr(3,3),r2l(3),r2r(3),qqa(6,6,6,3)
    real(wp)  ss(6,6,3),ddc(3,6,6,3),qqc(6,6,6,3),dda(3,6,6,3)
@@ -2089,10 +2094,10 @@ subroutine ccm_build_dSH0(nat,basis,thr,nao,nbf,at,xyz,q,cn,P,Pew,g,sigma, &
    !$omp private(iat,jat,ixyz,ati,cc,ci,rab2,atj,ish,ishtyp,valaoi,valaoj,g_xyz, &
    !$omp&        ri,rj,icao,naoi,iptyp,jsh,jshmax,jshtyp,jcao,naoj,jptyp,rij2, &
    !$omp&        sdq,sdqg,est,alpi,alpj,ab,iprim,jprim,ip,jp,km,shpoly,dshpoly, &
-   !$omp&        mli,mlj,dum,dumdum,tmp,stmp,rij,zi,zj,zetaij,enpoly,iao,jao, &
+   !$omp&        mli,mlj,dum,dumdum,tmp,stmp,rij,zi,zj,zetaij,enpoly,iao,jao,t, &
    !$omp&        ii,jj,k,den,den2,den4,i,j,il,jl,hii,hjj,hav,Pij,Hij,HPij,H0sr) &
    !$omp reduction (+:g,sigma,dhdcn,dhdq) &
-   !$omp shared(wsc,basis,at,ao_n,ao_l,xyz,thr,nat,kqat,kcnat,cn,q,en, &
+   !$omp shared(wsc,basis,at,ao_n,ao_l,xyz,lattice,thr,nat,kqat,kcnat,cn,q,en, &
    !$omp        ken,gam3,xbdamp,kcn,kmagic,kpair,alphaj,P,Pew)
    !$omp do schedule(runtime)
    do iat = 1, nat
@@ -2160,7 +2165,8 @@ subroutine ccm_build_dSH0(nat,basis,thr,nao,nbf,at,xyz,q,cn,P,Pew,g,sigma, &
                hav = 0.5_wp * km * (hii + hjj)
 
                wscAt: do kat = 1,wsc%itbl(jat,iat)
-                  rj = wsc%xyz(:,kat,jat,iat)
+                  t = wsc%lattr(:,kat,jat,iat)
+                  rj = xyz(:,jat) + matmul(lattice,t)
                   rij = ri - rj
                   rij2 = sum(rij**2)