diff --git a/src/tblite/basis/q-vszp.f90 b/src/tblite/basis/q-vszp.f90 index 9454368e..42cdd05e 100644 --- a/src/tblite/basis/q-vszp.f90 +++ b/src/tblite/basis/q-vszp.f90 @@ -31,7 +31,7 @@ module tblite_basis_qvszp use mctc_io, only : structure_type use tblite_basis_type, only : basis_type, integral_cutoff use tblite_basis_type, only : cgto_type - use tblite_integral_overlap, only : overlap_cgto, overlap_grad_cgto, msao + use tblite_integral_overlap, only : overlap_cgto, overlap_grad_decontracted_cgto, msao use tblite_ncoord, only : ncoord_type, new_ncoord implicit none private @@ -243,13 +243,13 @@ module tblite_basis_qvszp contains !> Create a new basis set - subroutine new_basis(self, mol, nsh_id, cgto, acc) + subroutine new_basis(self, mol, nshell, cgto, acc) !> Instance of the basis set data type(qvszp_basis_type), intent(out) :: self !> Molecular structure data type(structure_type), intent(in) :: mol !> Number of shells per species - integer, intent(in) :: nsh_id(:) + integer, intent(in) :: nshell(:) !> Contracted Gaussian basis functions for each shell and species class(qvszp_cgto_type), intent(in) :: cgto(:, :) !> Calculation accuracy @@ -258,12 +258,12 @@ subroutine new_basis(self, mol, nsh_id, cgto, acc) integer :: iat, isp, ish, iao, ii real(wp) :: min_alpha - self%nsh_id = nsh_id + self%nsh_id = nshell self%cgto = cgto self%intcut = integral_cutoff(acc) ! Make count of shells for each atom - self%nsh_at = nsh_id(mol%id) + self%nsh_at = nshell(mol%id) ! Create mapping between atoms and shells self%nsh = sum(self%nsh_at) @@ -302,7 +302,7 @@ subroutine new_basis(self, mol, nsh_id, cgto, acc) ii = 0 do iat = 1, mol%nat isp = mol%id(iat) - do ish = 1, nsh_id(isp) + do ish = 1, nshell(isp) self%iao_sh(ish+self%ish_at(iat)) = ii ii = ii + 2*cgto(ish, iat)%ang + 1 end do @@ -311,7 +311,7 @@ subroutine new_basis(self, mol, nsh_id, cgto, acc) min_alpha = huge(acc) do iat = 1, mol%nat isp = mol%id(iat) - do ish = 1, nsh_id(isp) + do ish = 1, nshell(isp) self%maxl = max(self%maxl, cgto(ish, iat)%ang) min_alpha = min(min_alpha, minval(cgto(ish, iat)%alpha(:cgto(ish, iat)%nprim))) end do @@ -319,6 +319,16 @@ subroutine new_basis(self, mol, nsh_id, cgto, acc) self%min_alpha = min_alpha + ! Allocate the coefficient derivatives + do iat = 1, mol%nat + isp = mol%id(iat) + do ish = 1, nshell(isp) + ! Allocate memory for the gradient of the coefficients + allocate(self%cgto(ish, iat)%dcoeffdr(3, mol%nat, maxg), & + & self%cgto(ish, iat)%dcoeffdL(3, 3, maxg), source=0.0_wp) + end do + end do + ! Create the coordination number for the environment specific scaling call new_ncoord(self%ncoord, mol, cn_type="erf", & & kcn=kcn, rcov=qvszp_cov_radii(mol%num)) @@ -367,10 +377,6 @@ subroutine add_qvszp_cgtos(self, mol, nsh_id, norm) self(ish, iat)%k1 = p_k1(izp) self(ish, iat)%k2 = p_k2(izp) self(ish, iat)%k3 = p_k3(izp) - - ! Allocate memory for the gradient of the coefficients - allocate(self(ish, iat)%dcoeffdr(3, mol%nat, maxg), self(ish, iat)%dcoeffdL(3, 3, maxg), source=0.0_wp) - end do end do @@ -380,7 +386,7 @@ subroutine add_qvszp_cgtos(self, mol, nsh_id, norm) isp = mol%id(iat) izp = mol%num(isp) do ish = 1, nsh_id(isp) - call normalize_cgto(self(ish, iat), .false.) + call normalize_cgto(self(ish, iat), mol%nat, .false.) end do end do end if @@ -415,7 +421,6 @@ subroutine scale_cgto(self, qat, cn, norm, expscal, dqatdr, dqatdL, dcndr, dcndL grad = .false. if(present(dqatdr) .and. present(dcndr) .and. present(dqatdL) .and. present(dcndL)) then - write(*,*) "in grad", 3, size(dqatdr, 2) allocate(dqeffdr(3, size(dqatdr, 2)), dqeffdL(3, 3)) grad = .true. end if @@ -445,11 +450,6 @@ subroutine scale_cgto(self, qat, cn, norm, expscal, dqatdr, dqatdL, dcndr, dcndL end if end do - !do ipr = 1, 3 - ! write(*,*) "ic", ipr - ! call write_2d_matrix(self%dcoeffdr(ipr, :, :), "dcoeffdr") - ! !write(*,*) self%dcoeffdr(:, :, :) - !end do ! Optionally scale the exponent if(present(expscal)) then do ipr = 1, maxg @@ -459,57 +459,54 @@ subroutine scale_cgto(self, qat, cn, norm, expscal, dqatdr, dqatdL, dcndr, dcndL ! Normalize the CGTOs if requested if(norm) then - call normalize_cgto(self, grad) + call normalize_cgto(self, size(dqatdr, 2), grad) end if - !do ipr = 1, 3 - ! write(*,*) "ic", ipr - ! call write_2d_matrix(self%dcoeffdr(ipr, :, :), "dcoeffdr last") - ! !write(*,*) self%dcoeffdr(:, :, :) - !end do end subroutine scale_cgto - subroutine normalize_cgto(self, grad) + subroutine normalize_cgto(self, nat, grad) !> Instance of the cgto data type(qvszp_cgto_type), intent(inout) :: self + !> Number of atoms in the molecule + integer, intent(in) :: nat !> Flag to normalize also the coefficient gradient logical, intent(in) :: grad - integer :: ipr, ic - real(wp) :: normalization, dnormalization(3), r2, vec(3) - real(wp) :: overlap(msao(max_shell), msao(max_shell)), doverlap(3, msao(max_shell), msao(max_shell)) + integer :: ipr + real(wp) :: normalization, r2, vec(3) + real(wp) :: overlap(msao(max_shell), msao(max_shell)) + real(wp), allocatable :: dnormalization(:, :), doverlap(:, :, :, :) r2 = 0.0_wp vec = 0.0_wp overlap = 0.0_wp - doverlap = 0.0_wp ! No screening has to be applied since the gaussians are on the same center if(grad) then + allocate(dnormalization(3, nat), doverlap(3, nat, msao(max_shell), msao(max_shell)), source=0.0_wp) ! The single center overlap has a derivative since it is not normalized at this point - call overlap_grad_cgto(self, self, r2, vec, 100.0_wp, overlap, doverlap) + call overlap_grad_decontracted_cgto(nat, self, self, r2, vec, 100.0_wp, overlap, doverlap) else call overlap_cgto(self, self, r2, vec, 100.0_wp, overlap) end if - do ic = 1, 3 - call write_2d_matrix(doverlap(ic, :, :), "doverlap") - end do - ! Normalization constant since all diagonal entries are equivalent normalization = 1 / sqrt(overlap(1, 1)) - dnormalization = -0.5_wp * doverlap(:, 1, 1) / sqrt(overlap(1, 1))**3 - - ! Normalize the coefficients and their derivatives - do ipr = 1, self%nprim - ! Derivative of the normalization w.r.t the position of the atom at which the CGTO is centered - self%dcoeffdr(:, self%atomidx, ipr) = self%dcoeffdr(:, self%atomidx, ipr) + self%coeff(ipr) * dnormalization - ! self%dcoeffdL(:,ipr) = self%dcoeffdL(:,ipr) + self%coeff(ipr) * dnormalization + if(grad) then + dnormalization = -0.5_wp * doverlap(:, :, 1, 1) / sqrt(overlap(1, 1))**3 + + do ipr = 1, self%nprim + ! Normalization of the coefficient derivative + self%dcoeffdr(:, :, ipr) = self%dcoeffdr(:, :, ipr) * normalization + ! self%dcoeffdL(:,ipr) = self%dcoeffdL(:,ipr) * normalization - ! Normalization of the coefficient derivative - self%dcoeffdr(:, :, ipr) = self%dcoeffdr(:, :, ipr) * normalization - ! self%dcoeffdL(:,ipr) = self%dcoeffdL(:,ipr) * normalization + ! Derivative of the normalization constant + self%dcoeffdr(:, :, ipr) = self%dcoeffdr(:, :, ipr) + self%coeff(ipr) * dnormalization + ! self%dcoeffdL(:,ipr) = self%dcoeffdL(:,ipr) + self%coeff(ipr) * dnormalization + end do + end if + do ipr = 1, self%nprim ! Normalization of the coefficient self%coeff(ipr) = self%coeff(ipr) * normalization end do diff --git a/src/tblite/integral/overlap.f90 b/src/tblite/integral/overlap.f90 index af24c83b..99553bbb 100644 --- a/src/tblite/integral/overlap.f90 +++ b/src/tblite/integral/overlap.f90 @@ -33,7 +33,7 @@ module tblite_integral_overlap implicit none private - public :: overlap_cgto, overlap_cgto_diat, overlap_grad_cgto, overlap_grad_cgto_diat + public :: overlap_cgto, overlap_cgto_diat, overlap_grad_cgto, overlap_grad_cgto_diat, overlap_grad_decontracted_cgto public :: get_overlap, get_overlap_gradient public :: maxl, msao, smap @@ -462,81 +462,6 @@ pure subroutine overlap_cgto_diat(cgtoj, cgtoi, r2, vec, intcut, & end subroutine overlap_cgto_diat - -! !pure -! subroutine overlap_grad_cgto(cgtoj, cgtoi, r2, vec, intcut, overlap, doverlap) -! !> Description of contracted Gaussian function on center j -! class(cgto_type), intent(in) :: cgtoj -! !> Description of contracted Gaussian function on center i -! class(cgto_type), intent(in) :: cgtoi -! !> Square distance between center i and j -! real(wp), intent(in) :: r2 -! !> Distance vector between center i and j, ri - rj -! real(wp), intent(in) :: vec(3) -! !> Maximum value of integral prefactor to consider -! real(wp), intent(in) :: intcut -! !> Overlap integrals for the given pair i and j -! real(wp), intent(out) :: overlap(msao(cgtoj%ang), msao(cgtoi%ang)) -! !> Overlap integral gradient for the given pair i and j -! real(wp), intent(out) :: doverlap(3, msao(cgtoj%ang), msao(cgtoi%ang)) - -! integer :: ip, jp, mli, mlj, l -! real(wp) :: eab, oab, est, s1d(0:maxl2), rpi(3), rpj(3), cc, val, grad(3), pre, dcc(3) -! real(wp) :: s3d(mlao(cgtoj%ang), mlao(cgtoi%ang)) -! real(wp) :: ds3d(3, mlao(cgtoj%ang), mlao(cgtoi%ang)) -! real(wp) :: ds3d_delta(3, mlao(cgtoj%ang), mlao(cgtoi%ang)) -! real(wp) :: temp1(3, msao(cgtoj%ang), msao(cgtoi%ang)) -! real(wp) :: temp2(3, msao(cgtoj%ang), msao(cgtoi%ang)) - - -! s3d(:, :) = 0.0_wp -! ds3d(:, :, :) = 0.0_wp -! ds3d_delta(:, :, :) = 0.0_wp -! temp1(:, :, :) = 0.0_wp -! temp2(:, :, :) = 0.0_wp - -! do ip = 1, cgtoi%nprim -! do jp = 1, cgtoj%nprim -! eab = cgtoi%alpha(ip) + cgtoj%alpha(jp) -! oab = 1.0_wp/eab -! est = cgtoi%alpha(ip) * cgtoj%alpha(jp) * r2 * oab -! if (est > intcut) cycle -! pre = exp(-est) * sqrtpi3*sqrt(oab)**3 -! rpi = -vec * cgtoj%alpha(jp) * oab -! rpj = +vec * cgtoi%alpha(ip) * oab -! do l = 0, cgtoi%ang + cgtoj%ang + 1 -! s1d(l) = overlap_1d(l, eab) -! end do -! cc = cgtoi%coeff(ip) * cgtoj%coeff(jp) * pre -! !dcc = (cgtoi%dcoeffdr(:,ip) * cgtoj%coeff(jp) + cgtoi%coeff(ip) * cgtoj%dcoeffdr(:,jp)) * pre -! do mli = 1, mlao(cgtoi%ang) -! do mlj = 1, mlao(cgtoj%ang) -! call overlap_grad_3d(rpj, rpi, cgtoj%alpha(jp), cgtoi%alpha(ip), & -! & lx(:, mlj+lmap(cgtoj%ang)), lx(:, mli+lmap(cgtoi%ang)), & -! & s1d, val, grad) -! write(*,*) "dcc*val", dcc*val, dcc, val -! s3d(mlj, mli) = s3d(mlj, mli) + cc*val -! ds3d(:, mlj, mli) = ds3d(:, mlj, mli) + cc*grad -! !ds3d_delta(:, mlj, mli) = ds3d_delta(:, mlj, mli) + dcc*val -! end do -! end do -! end do -! end do - -! call transform0(cgtoj%ang, cgtoi%ang, s3d, overlap) - -! call transform1(cgtoj%ang, cgtoi%ang, ds3d, temp1) -! call transform1(cgtoj%ang, cgtoi%ang, ds3d_delta, temp2) -! call write_2d_matrix(temp1(3, :, :), "temp1") -! call write_2d_matrix(temp2(3, :, :), "temp2") -! !call transform1(cgtoj%ang, cgtoi%ang, ds3d, doverlap) -! !call transform1(cgtoj%ang, cgtoi%ang, ds3d_delta, doverlap) - -! doverlap = temp1 !+ temp2 -! call write_2d_matrix(doverlap(3, :, :), "doverlap") - -! end subroutine overlap_grad_cgto - ! pure subroutine overlap_grad_cgto(cgtoj, cgtoi, r2, vec, intcut, overlap, doverlap) !> Description of contracted Gaussian function on center j @@ -555,74 +480,102 @@ subroutine overlap_grad_cgto(cgtoj, cgtoi, r2, vec, intcut, overlap, doverlap) real(wp), intent(out) :: doverlap(3, msao(cgtoj%ang), msao(cgtoi%ang)) integer :: ip, jp, mli, mlj, l - real(wp) :: eab, oab, est, s1d(0:maxl2), rpi(3), rpj(3), cc, val, grad(3), pre, dcc(3) + real(wp) :: eab, oab, est, s1d(0:maxl2), rpi(3), rpj(3), cc, val, grad(3), pre real(wp) :: s3d(mlao(cgtoj%ang), mlao(cgtoi%ang)) real(wp) :: ds3d(3, mlao(cgtoj%ang), mlao(cgtoi%ang)) s3d(:, :) = 0.0_wp ds3d(:, :, :) = 0.0_wp - if(allocated(cgtoi%dcoeffdr) .and. allocated(cgtoj%dcoeffdr)) then - write(*,*) "doing the gradient" - do ip = 1, cgtoi%nprim - do jp = 1, cgtoj%nprim - eab = cgtoi%alpha(ip) + cgtoj%alpha(jp) - oab = 1.0_wp/eab - est = cgtoi%alpha(ip) * cgtoj%alpha(jp) * r2 * oab - if (est > intcut) cycle - pre = exp(-est) * sqrtpi3*sqrt(oab)**3 - rpi = -vec * cgtoj%alpha(jp) * oab - rpj = +vec * cgtoi%alpha(ip) * oab - do l = 0, cgtoi%ang + cgtoj%ang + 1 - s1d(l) = overlap_1d(l, eab) - end do - cc = cgtoi%coeff(ip) * cgtoj%coeff(jp) * pre - write(*,*) "atom idx", cgtoi%atomidx, cgtoj%atomidx - dcc = (cgtoi%dcoeffdr(:,cgtoi%atomidx,ip) * cgtoj%coeff(jp) & - & + cgtoi%coeff(ip) * cgtoj%dcoeffdr(:,cgtoj%atomidx,jp)) * pre - do mli = 1, mlao(cgtoi%ang) - do mlj = 1, mlao(cgtoj%ang) - call overlap_grad_3d(rpj, rpi, cgtoj%alpha(jp), cgtoi%alpha(ip), & - & lx(:, mlj+lmap(cgtoj%ang)), lx(:, mli+lmap(cgtoi%ang)), & - & s1d, val, grad) - s3d(mlj, mli) = s3d(mlj, mli) + cc*val - ds3d(:, mlj, mli) = ds3d(:, mlj, mli) + cc*grad + dcc*val - end do - end do + do ip = 1, cgtoi%nprim + do jp = 1, cgtoj%nprim + eab = cgtoi%alpha(ip) + cgtoj%alpha(jp) + oab = 1.0_wp/eab + est = cgtoi%alpha(ip) * cgtoj%alpha(jp) * r2 * oab + if (est > intcut) cycle + pre = exp(-est) * sqrtpi3*sqrt(oab)**3 + rpi = -vec * cgtoj%alpha(jp) * oab + rpj = +vec * cgtoi%alpha(ip) * oab + do l = 0, cgtoi%ang + cgtoj%ang + 1 + s1d(l) = overlap_1d(l, eab) end do - end do - else - do ip = 1, cgtoi%nprim - do jp = 1, cgtoj%nprim - eab = cgtoi%alpha(ip) + cgtoj%alpha(jp) - oab = 1.0_wp/eab - est = cgtoi%alpha(ip) * cgtoj%alpha(jp) * r2 * oab - if (est > intcut) cycle - pre = exp(-est) * sqrtpi3*sqrt(oab)**3 - rpi = -vec * cgtoj%alpha(jp) * oab - rpj = +vec * cgtoi%alpha(ip) * oab - do l = 0, cgtoi%ang + cgtoj%ang + 1 - s1d(l) = overlap_1d(l, eab) - end do - cc = cgtoi%coeff(ip) * cgtoj%coeff(jp) * pre - do mli = 1, mlao(cgtoi%ang) - do mlj = 1, mlao(cgtoj%ang) - call overlap_grad_3d(rpj, rpi, cgtoj%alpha(jp), cgtoi%alpha(ip), & - & lx(:, mlj+lmap(cgtoj%ang)), lx(:, mli+lmap(cgtoi%ang)), & - & s1d, val, grad) - s3d(mlj, mli) = s3d(mlj, mli) + cc*val - ds3d(:, mlj, mli) = ds3d(:, mlj, mli) + cc*grad - end do + cc = cgtoi%coeff(ip) * cgtoj%coeff(jp) * pre + do mli = 1, mlao(cgtoi%ang) + do mlj = 1, mlao(cgtoj%ang) + call overlap_grad_3d(rpj, rpi, cgtoj%alpha(jp), cgtoi%alpha(ip), & + & lx(:, mlj+lmap(cgtoj%ang)), lx(:, mli+lmap(cgtoi%ang)), & + & s1d, val, grad) + s3d(mlj, mli) = s3d(mlj, mli) + cc*val + ds3d(:, mlj, mli) = ds3d(:, mlj, mli) + cc*grad end do end do end do - end if + end do call transform0(cgtoj%ang, cgtoi%ang, s3d, overlap) call transform1(cgtoj%ang, cgtoi%ang, ds3d, doverlap) end subroutine overlap_grad_cgto +! pure +subroutine overlap_grad_decontracted_cgto(nat, cgtoj, cgtoi, r2, vec, intcut, overlap, doverlap) + !> Number of atoms in the molecule + integer, intent(in) :: nat + !> Description of contracted Gaussian function on center j + class(cgto_type), intent(in) :: cgtoj + !> Description of contracted Gaussian function on center i + class(cgto_type), intent(in) :: cgtoi + !> Square distance between center i and j + real(wp), intent(in) :: r2 + !> Distance vector between center i and j, ri - rj + real(wp), intent(in) :: vec(3) + !> Maximum value of integral prefactor to consider + real(wp), intent(in) :: intcut + !> Overlap integrals for the given pair i and j + real(wp), intent(out) :: overlap(msao(cgtoj%ang), msao(cgtoi%ang)) + !> Overlap integral gradient for the given pair i and j + real(wp), intent(out) :: doverlap(3, nat, msao(cgtoj%ang), msao(cgtoi%ang)) + + integer :: ip, jp, mli, mlj, l + real(wp) :: eab, oab, est, s1d(0:maxl2), rpi(3), rpj(3), cc, val, grad(3), pre, dcc(3,nat) + real(wp) :: s3d(mlao(cgtoj%ang), mlao(cgtoi%ang)) + real(wp) :: ds3d(3, nat, mlao(cgtoj%ang), mlao(cgtoi%ang)) + + s3d(:, :) = 0.0_wp + ds3d(:, :, :, :) = 0.0_wp + + do ip = 1, cgtoi%nprim + do jp = 1, cgtoj%nprim + eab = cgtoi%alpha(ip) + cgtoj%alpha(jp) + oab = 1.0_wp/eab + est = cgtoi%alpha(ip) * cgtoj%alpha(jp) * r2 * oab + if (est > intcut) cycle + pre = exp(-est) * sqrtpi3*sqrt(oab)**3 + rpi = -vec * cgtoj%alpha(jp) * oab + rpj = +vec * cgtoi%alpha(ip) * oab + do l = 0, cgtoi%ang + cgtoj%ang + 1 + s1d(l) = overlap_1d(l, eab) + end do + cc = cgtoi%coeff(ip) * cgtoj%coeff(jp) * pre + dcc = (cgtoi%dcoeffdr(:,:,ip) * cgtoj%coeff(jp) & + & + cgtoi%coeff(ip) * cgtoj%dcoeffdr(:,:,jp)) * pre + do mli = 1, mlao(cgtoi%ang) + do mlj = 1, mlao(cgtoj%ang) + call overlap_3d(rpj, rpi, cgtoj%alpha(jp), cgtoi%alpha(ip), & + & lx(:, mlj+lmap(cgtoj%ang)), lx(:, mli+lmap(cgtoi%ang)), & + & s1d, val) + s3d(mlj, mli) = s3d(mlj, mli) + cc*val + ds3d(:, :, mlj, mli) = ds3d(:, :, mlj, mli) + dcc*val + end do + end do + end do + end do + + call transform0(cgtoj%ang, cgtoi%ang, s3d, overlap) + call transform2(cgtoj%ang, cgtoi%ang, ds3d, doverlap) + +end subroutine overlap_grad_decontracted_cgto + pure subroutine overlap_grad_cgto_diat(cgtoj, cgtoi, r2, vec, intcut, & & ksig, kpi, kdel, overlap, doverlap, overlap_diat, doverlap_diat) diff --git a/test/unit/main.f90 b/test/unit/main.f90 index 4c947dbb..32257167 100644 --- a/test/unit/main.f90 +++ b/test/unit/main.f90 @@ -68,7 +68,7 @@ program tester !new_testsuite("double-dictionary", collect_double_dictionary), & !new_testsuite("post-processing", collect_post_processing), & !new_testsuite("slater-expansion", collect_slater_expansion), & - new_testsuite("qvszp", collect_qvszp) & + new_testsuite("q-vszp", collect_qvszp) & !new_testsuite("cgto-ortho", collect_cgto_ortho), & !new_testsuite("integral-overlap", collect_integral_overlap), & !new_testsuite("integral-multipole", collect_integral_multipole), & diff --git a/test/unit/test_q_vszp.f90 b/test/unit/test_q_vszp.f90 index 20382c03..02b83aeb 100644 --- a/test/unit/test_q_vszp.f90 +++ b/test/unit/test_q_vszp.f90 @@ -56,15 +56,15 @@ subroutine collect_qvszp(testsuite) type(unittest_type), allocatable, intent(out) :: testsuite(:) testsuite = [ & - !new_unittest("norm-qvszp-mb01", test_qvszp_norm_cgto_mb01), & - !new_unittest("norm-qvszp-accl3", test_qvszp_norm_cgto_accl3), & - !new_unittest("overlap-qvszp-lih", test_qvszp_overlap_lih), & - !new_unittest("overlap-qvszp-sih4", test_qvszp_overlap_sih4), & - !new_unittest("overlap-qvszp-accl3", test_qvszp_overlap_accl3), & + new_unittest("norm-qvszp-mb01", test_qvszp_norm_cgto_mb01), & + new_unittest("norm-qvszp-accl3", test_qvszp_norm_cgto_accl3), & + new_unittest("overlap-qvszp-lih", test_qvszp_overlap_lih), & + new_unittest("overlap-qvszp-sih4", test_qvszp_overlap_sih4), & + new_unittest("overlap-qvszp-accl3", test_qvszp_overlap_accl3), & new_unittest("coefficient-grad-qvszp-lih", test_qvszp_coefficient_grad_lih), & new_unittest("coefficient-grad-qvszp-sih4", test_qvszp_coefficient_grad_sih4), & - new_unittest("coefficient-grad-qvszp-accl3", test_qvszp_coefficient_grad_accl3), & - new_unittest("overlap-diat-grad-qvszp-lih", test_qvszp_overlap_diat_grad_lih) & + new_unittest("coefficient-grad-qvszp-accl3", test_qvszp_coefficient_grad_accl3) & + !new_unittest("overlap-diat-grad-qvszp-lih", test_qvszp_overlap_diat_grad_lih) & & ] end subroutine collect_qvszp @@ -311,7 +311,7 @@ subroutine test_qvszp_coefficient_grad(mol, error) ! Numerical gradient of the hamiltonian matrix do jat = 1, mol%nat - jzp = mol%id(iat) + jzp = mol%id(jat) do jsh = 1, basr%nsh_id(jzp) ! Upper triangular matrix and diagonal numdr(ic, iat, jat, jsh, :) = & @@ -342,24 +342,18 @@ subroutine test_qvszp_coefficient_grad(mol, error) do iat = 1, mol%nat do jat = 1, mol%nat jzp = mol%id(jat) - !write(*,*) "dqatnumdr", dqatnumdr(ic, iat, :) do jsh = 1, basr%nsh_id(jzp) - do jpr = 1, basr%cgto(jsh,jat)%nprim - write(*,*) "CURRENT DIMENSION", ic, iat, jat, jsh, jpr - write(*,*) "numdr", numdr(ic, iat, jat, jsh, jpr) - write(*,*) "analytical", basr%cgto(jsh,jat)%dcoeffdr(ic,iat,jpr) call check(error, numdr(ic, iat, jat, jsh, jpr), basr%cgto(jsh,jat)%dcoeffdr(ic,iat,jpr), thr=thr) if (allocated(error)) then call test_failed(error, "Coefficient derivative does not match") - !exit num + exit num end if end do end do end do end do end do num - write(*,*) "almost at the end" end subroutine test_qvszp_coefficient_grad