q-vSZP basis with 3 x nat size of the coefficient derivative (current…

…ly still with a double free error).

src/tblite/basis/CMakeLists.txt

@@ -21,6 +21,7 @@ list(
   "${dir}/ortho.f90"
   "${dir}/slater.f90"
   "${dir}/type.f90"
+  "${dir}/q-vszp.f90"
 )
 
 set(srcs "${srcs}" PARENT_SCOPE)

src/tblite/basis/meson.build

@@ -18,4 +18,5 @@ srcs += files(
   'ortho.f90',
   'slater.f90',
   'type.f90',
+  'q-vszp.f90',
 )

src/tblite/basis/type.f90

@@ -24,7 +24,7 @@ module tblite_basis_type
    implicit none
    private
 
-   public :: new_basis, get_cutoff
+   public :: new_basis, get_cutoff, integral_cutoff, maxg
 
    !> Maximum contraction length of basis functions.
    !> The limit is chosen as twice the maximum size returned by the STO-NG expansion
@@ -36,11 +36,21 @@ module tblite_basis_type
       integer :: ang = -1
       !> Contraction length of this basis function
       integer :: nprim = 0
+      !> Index of the atom at which the CGTO is centered,
+      !> for use in the coefficient gradient (default is not specified)
+      integer :: atomidx = -1 
       !> Exponent of the primitive Gaussian functions
       real(wp) :: alpha(maxg) = 0.0_wp
       !> Contraction coefficients of the primitive Gaussian functions,
       !> might contain normalization
       real(wp) :: coeff(maxg) = 0.0_wp
+      !> Derivative of the Contraction coefficient w.r.t. the atomic positions
+      real(wp), allocatable :: dcoeffdr(:,:,:)
+      !> Derivative of the Contraction coefficient w.r.t. L
+      real(wp), allocatable :: dcoeffdL(:,:,:)
+   contains 
+      !> Scales the coefficient of the cgto 
+      procedure :: scale_cgto
    end type cgto_type
 
    !> Collection of information regarding the basis set of a system
@@ -73,7 +83,10 @@ module tblite_basis_type
       !> Mapping from shells to the respective atom
       integer, allocatable :: sh2at(:)
       !> Contracted Gaussian basis functions forming the basis set
-      type(cgto_type), allocatable :: cgto(:, :)
+      class(cgto_type), allocatable :: cgto(:, :)
+   contains
+      !> Scales the coefficient of the basis
+      procedure :: scale_basis
    end type basis_type
 
    !> Get optimal real space cutoff for integral evaluation
@@ -162,10 +175,76 @@ subroutine new_basis(self, mol, nshell, cgto, acc)
 
 end subroutine new_basis
 
+!> Scale the basis set contraction based on the atomic environment and charge
+subroutine scale_basis(self, mol, qat, cn, norm, expscal, dqatdr, dqatdL, dcndr, dcndL)
+   !> Instance of the basis set data
+   class(basis_type), intent(inout) :: self
+   !> Molecular structure data
+   type(structure_type), intent(in) :: mol
+   !> Atomic charges for the charge scaling of the basis set 
+   real(wp), intent(in) :: qat(:)
+   !> Coordination number
+   real(wp), intent(in) :: cn(:)
+   !> Include normalization in contraction coefficients
+   logical, intent(in) :: norm   
+   !> Exponent scaling factor
+   real(wp), intent(in), optional :: expscal
+   !> Derivative of the specific atomic charge w.r.t. the positions
+   real(wp), intent(in), optional :: dqatdr(:, :, :)
+   !> Derivative of the specific atomic charge w.r.t. the lattice vectors
+   real(wp), intent(in), optional :: dqatdL(:, :, :)
+   !> Derivative of the specific coordination number w.r.t. the positions
+   real(wp), intent(in), optional :: dcndr(:, :, :)
+   !> Derivative of the specific coordination number w.r.t. the lattice vectors
+   real(wp), intent(in), optional :: dcndL(:, :, :)
+
+   integer :: iat, ish
+
+   ! Scale the coefficients for all atoms
+   if(present(dqatdr) .and. present(dqatdL) .and. present(dcndr) .and. present(dcndL)) then
+      do iat = 1, mol%nat
+         do ish = 1, self%nsh_at(iat)
+            call self%cgto(ish, iat)%scale_cgto(qat(iat), cn(iat), norm, expscal, &
+               & dqatdr(:, :, iat), dqatdL(:, :, iat), dcndr(:, :, iat), dcndL(:, :, iat))
+         end do
+      end do
+   else 
+      do iat = 1, mol%nat
+         do ish = 1, self%nsh_at(iat)
+            call self%cgto(ish, iat)%scale_cgto(qat(iat), cn(iat), norm, expscal)
+         end do
+      end do
+   end if
+
+end subroutine scale_basis
+
+!> Dummy procedure to scale the contraction coefficients of the CGTOs
+subroutine scale_cgto(self, qat, cn, norm, expscal, dqatdr, dqatdL, dcndr, dcndL)
+   !> Instance of the basis set data
+   class(cgto_type(*)), intent(inout) :: self
+   !> Atomic charges for the charge scaling of the basis set 
+   real(wp), intent(in) :: qat
+   !> Coordination number
+   real(wp), intent(in) :: cn
+   !> Include normalization in contraction coefficients
+   logical, intent(in) :: norm
+   !> Exponent scaling factor
+   real(wp), intent(in), optional :: expscal
+   !> Derivative of the specific atomic charge w.r.t. the positions (not spin-resolved)
+   real(wp), intent(in), optional :: dqatdr(:, :)
+   !> Derivative of the specific atomic charge w.r.t. the lattice vectors (not spin-resolved)
+   real(wp), intent(in), optional :: dqatdL(:, :)
+   !> Derivative of the specific coordination number w.r.t. the positions
+   real(wp), intent(in), optional :: dcndr(:, :)
+   !> Derivative of the specific coordination number w.r.t. the lattice vectors
+   real(wp), intent(in), optional :: dcndL(:, :)
+
+end
+
 !> Determine required real space cutoff for the basis set
 pure function get_cutoff(self, acc) result(cutoff)
    !> Instance of the basis set data
-   type(basis_type), intent(in) :: self
+   class(basis_type), intent(in) :: self
    !> Accuracy for the integral cutoff
    real(wp), intent(in), optional :: acc
    !> Required realspace cutoff

src/tblite/ceh/ceh.f90

@@ -880,9 +880,6 @@ module tblite_ceh_ceh
    !> Angular momentum-specific scaling factors for H0 
    real(wp), parameter   :: kll(1:4) = [0.6379_wp, 0.9517_wp, 1.18_wp, 2.84_wp]
 
-   !> Conversion constant
-   real(wp), parameter   :: kt = 3.166808578545117e-06_wp
-
    !> Specification of the CEH hamiltonian
    type, extends(tb_h0spec) :: ceh_h0spec
    contains

src/tblite/ceh/h0.f90

@@ -129,7 +129,7 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, &
 
       integer  :: itr, img, inl, ii, jj, is, js, jzp, izp, nao
       integer  :: iat, ish, jat, jsh, k, iao, jao, ij, iaosh, jaosh
-      real(wp) :: hij, rr, r2, vec(3), dtmpj(3)
+      real(wp) :: hij, r2, vec(3), dtmpj(3)
       real(wp), allocatable :: stmp(:), dtmpi(:, :), block_overlap(:,:)
       integer :: kl, l
 
@@ -145,7 +145,7 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, &
       !$omp parallel do schedule(runtime) default(none) &
       !$omp shared(mol, bas, trans, list, overlap, overlap_diat, dpint, hamiltonian, h0, selfenergy) &
       !$omp private(iat, jat, izp, jzp, itr, is, js, ish, jsh, ii, jj, iao, jao, nao, ij, iaosh, jaosh, k) &
-      !$omp private(r2, vec, stmp, block_overlap, dtmpi, dtmpj, hij, rr, inl, img)
+      !$omp private(r2, vec, stmp, block_overlap, dtmpi, dtmpj, hij, inl, img)
       do iat = 1, mol%nat
          izp = mol%id(iat)
          is = bas%ish_at(iat)
@@ -157,7 +157,6 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, &
             js = bas%ish_at(jat)
             vec(:) = mol%xyz(:, iat) - mol%xyz(:, jat) - trans(:, itr)
             r2 = vec(1)**2 + vec(2)**2 + vec(3)**2
-            rr = sqrt(sqrt(r2) / (h0%rad(jzp) + h0%rad(izp)))
 
             ! Get the overlap and dipole integrals for the current diatomic pair
             block_overlap = 0.0_wp
@@ -258,13 +257,12 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, &
       !$omp parallel do schedule(runtime) default(none) &
       !$omp shared(mol, bas, trans, overlap, overlap_diat, dpint, hamiltonian, h0, selfenergy) &
       !$omp private(iat, izp, itr, is, ish, jsh, ii, jj, iao, jao, nao, ij) &
-      !$omp private(r2, vec, stmp, dtmpi, hij, rr)
+      !$omp private(r2, vec, stmp, dtmpi, hij)
       do iat = 1, mol%nat
          izp = mol%id(iat)
          is = bas%ish_at(iat)
          vec(:) = 0.0_wp
          r2 = 0.0_wp
-         rr = sqrt(sqrt(r2) / (h0%rad(izp) + h0%rad(izp)))
          do ish = 1, bas%nsh_id(izp)
             ii = bas%iao_sh(is+ish)
             do jsh = 1, bas%nsh_id(izp)
@@ -275,9 +273,12 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, &
                do iao = 1, msao(bas%cgto(ish, izp)%ang)
                   do jao = 1, nao
                      ij = jao + nao*(iao-1)
+
+                     !$omp atomic
                      overlap(jj+jao, ii+iao) = overlap(jj+jao, ii+iao) &
                         + stmp(ij)
 
+                     !$omp atomic
                      overlap_diat(jj+jao, ii+iao) = overlap_diat(jj+jao, ii+iao) &
                         + stmp(ij)
 
@@ -316,18 +317,18 @@ subroutine get_hamiltonian_gradient(mol, trans, list, bas, h0, selfenergy, &
       real(wp), intent(in) :: dsedL(:, :, :)
       !> Density dependent potential shifts on the Hamiltonian
       type(potential_type), intent(in) :: pot
-      !> Derivative of the electronic energy w.r.t. coordinate displacements
+      !> Derivative of the overlap matrix w.r.t. coordinates
       real(wp), intent(inout) :: doverlap(:, :, :)
-      !> Derivative of the electronic energy w.r.t. coordinate displacements
+      !> Derivative of the diatomic frame-scaled overlap matrix w.r.t. coordinates
       real(wp), intent(inout) :: doverlap_diat(:, :, :)
-      !> Derivative of the electronic energy w.r.t. coordinate displacements
+      !> Derivative of the Hamiltonian matrix w.r.t. coordinates
       real(wp), intent(inout) :: dh0dr(:, :, :)
-      !> Derivative of the electronic energy w.r.t. the lattice vector
+      !> Derivative of the Hamiltonian matrix w.r.t. the lattice vector
       real(wp), intent(inout) :: dh0dL(:, :, :)
 
       integer :: iat, jat, izp, jzp, itr, img, inl, spin, nspin
       integer :: ish, jsh, is, js, ii, jj, iao, jao, nao, ij, iaosh, jaosh
-      real(wp) :: rr, r2, vec(3), hij, vij
+      real(wp) :: r2, vec(3), hij, vij
       real(wp), allocatable :: stmp(:), dstmp(:, :)
       real(wp), allocatable :: block_overlap(:,:), block_doverlap(:,:,:)
       real(wp), allocatable :: dhijdr(:,:), dhijdL(:,:), dvijdr(:,:), dvijdL(:,:)
@@ -344,7 +345,7 @@ subroutine get_hamiltonian_gradient(mol, trans, list, bas, h0, selfenergy, &
       !$omp parallel do schedule(runtime) default(none) reduction(+: doverlap, doverlap_diat, dh0dr, dh0dL) &
       !$omp shared(mol, bas, trans, h0, selfenergy, dsedr, dsedL, pot, list) &
       !$omp private(iat, jat, izp, jzp, itr, is, js, ish, jsh, ii, jj, iao, jao, nao, ij, iaosh, jaosh, &
-      !$omp& r2, vec, stmp, dstmp, block_overlap, block_doverlap, hij, vij, dhijdr, dhijdL, dhvjdr, dvijdL, rr, inl, img)
+      !$omp& r2, vec, stmp, dstmp, block_overlap, block_doverlap, hij, vij, dhijdr, dhijdL, dvijdr, dvijdL, inl, img)
       do iat = 1, mol%nat
          izp = mol%id(iat)
          is = bas%ish_at(iat)
@@ -357,7 +358,6 @@ subroutine get_hamiltonian_gradient(mol, trans, list, bas, h0, selfenergy, &
             if (iat == jat) cycle
             vec(:) = mol%xyz(:, iat) - mol%xyz(:, jat) - trans(:, itr)
             r2 = vec(1)**2 + vec(2)**2 + vec(3)**2
-            rr = sqrt(sqrt(r2) / (h0%rad(jzp) + h0%rad(izp)))
 
             ! Setup the potential intermediate for the current atom pair
             vij = - 0.5_wp * (pot%vat(iat,1) + pot%vat(jat,1))
@@ -373,7 +373,7 @@ subroutine get_hamiltonian_gradient(mol, trans, list, bas, h0, selfenergy, &
                do jsh = 1, bas%nsh_id(jzp)
                   jj = bas%iao_sh(js+jsh)
                   jaosh = smap(jsh-1) ! Offset for the block overlap matrix
-
+                  
                   call overlap_grad_cgto(bas%cgto(jsh,jzp), bas%cgto(ish,izp), r2, vec, &
                   & bas%intcut, stmp, dstmp)
 
@@ -387,21 +387,16 @@ subroutine get_hamiltonian_gradient(mol, trans, list, bas, h0, selfenergy, &
                         block_overlap(jaosh+jao, iaosh+iao) = block_overlap(jaosh+jao, iaosh+iao) &
                            + stmp(ij)
 
-                        !$omp atomic
                         block_doverlap(:, jaosh+jao, iaosh+iao) = block_doverlap(:, jaosh+jao, iaosh+iao) &
                            + dstmp(:, ij)
 
-                        !$omp atomic
                         doverlap(:, jj+jao, ii+iao) = doverlap(:, jj+jao, ii+iao) &
                            & + dstmp(:,ij)
-                        !$omp atomic
                         doverlap(:, ii+iao, jj+jao) = doverlap(:, ii+iao, jj+jao) &
                            & - dstmp(:,ij)
 
-                        !$omp atomic
                         dh0dr(:, jj+jao, ii+iao) = dh0dr(:, jj+jao, ii+iao) &
                            & + dstmp(:, ij) * vij + stmp(ij) * dvijdr(:,iat)
-                        !$omp atomic
                         dh0dr(:, ii+iao, jj+jao) = dh0dr(:, ii+iao, jj+jao) &
                            & - dstmp(:, ij) * vij - stmp(ij) * dvijdr(:,iat)
 
@@ -457,8 +452,8 @@ subroutine get_hamiltonian_gradient(mol, trans, list, bas, h0, selfenergy, &
       end do
 
       !$omp parallel do schedule(runtime) default(none) reduction(+:dh0dr) &
-      !$omp shared(mol, bas, dsedr) &
-      !$omp private(iat, izp, jzp, is, ish, ii, iao)
+      !$omp shared(mol, bas, dsedr, pot) &
+      !$omp private(iat, izp, jzp, is, ish, ii, iao, dvijdr, dvijdL)
       do iat = 1, mol%nat
          izp = mol%id(iat)
          is = bas%ish_at(iat)