Skip to content

Commit

Permalink
Fix q-vSZP coefficient graident for full 3 x N dimensions. Workaround…
Browse files Browse the repository at this point in the history 
… for the double deallocate error.
  • Loading branch information
@thfroitzheim
thfroitzheim committed Jul 30, 2024
1 parent 0656054 commit 0283c3a
Show file tree
Hide file tree
Showing 4 changed files with 132 additions and 188 deletions.
85 changes: 41 additions & 44 deletions src/tblite/basis/q-vszp.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -311,14 +311,24 @@ 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
end do

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))
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
Loading

0 comments on commit 0283c3a

Please sign in to comment.