Skip to content

Commit

Permalink
Fix for integral transformation and diatomic frame overlap (tblite#186)
Browse files Browse the repository at this point in the history 
* Further preparation for the new CEH version. Fix d-orbital ordering, fix diatomic-frame derivative, update Hamiltonian construction in CEH

* Update also the meson mstore version
  • Loading branch information
@thfroitzheim
thfroitzheim authored Jul 30, 2024
1 parent d2f08e8 commit 79d8f9a
Show file tree
Hide file tree
Showing 11 changed files with 2,986 additions and 1,258 deletions.
96 changes: 60 additions & 36 deletions src/tblite/ceh/h0.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,9 +25,9 @@ module tblite_ceh_h0
use tblite_xtb_spec, only : tb_h0spec
use tblite_xtb_h0, only : tb_hamiltonian
use tblite_integral_dipole, only: get_dipole_integrals, dipole_cgto, &
& dipole_cgto_diat_scal, maxl, msao
use tblite_integral_diat_trafo, only: relvec
& maxl, msao, smap
use tblite_adjlist, only : adjacency_list
use tblite_integral_diat_trafo, only: diat_trafo

implicit none
private
Expand Down Expand Up @@ -105,15 +105,15 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, &
real(wp), intent(out) :: overlap(:, :)
!> Dipole moment integral matrix
real(wp), intent(out) :: dpint(:, :, :)
!> Diatomic frame scaled overlap integral matrix
!> Diatomic frame-scaled overlap integral matrix
real(wp), intent(out) :: overlap_diat(:, :)
!> Effective Hamiltonian
real(wp), intent(out) :: hamiltonian(:, :)

integer :: itr, img, inl, ii, jj, is, js, jzp, izp, nao
integer :: iat, ish, jat, jsh, k, iao, jao, ij
real(wp) :: hij, rr, r2, vec(3), vec_diat_trafo(3), dtmpj(3)
real(wp), allocatable :: stmp(:), dtmpi(:, :), stmp_diat(:)
integer :: iat, ish, jat, jsh, k, iao, jao, ij, iaosh, jaosh
real(wp) :: hij, rr, r2, vec(3), dtmpj(3)
real(wp), allocatable :: stmp(:), dtmpi(:, :), block_overlap(:,:)
integer :: kl, l

overlap(:, :) = 0.0_wp
Expand All @@ -122,12 +122,13 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, &
hamiltonian(:, :) = 0.0_wp

! Allocate temporary matrices for overlap, diatomic frame scaled overlap and dipole moment integrals
allocate(stmp(msao(bas%maxl)**2), stmp_diat(msao(bas%maxl)**2), dtmpi(3, msao(bas%maxl)**2))
allocate(stmp(msao(bas%maxl)**2), dtmpi(3, msao(bas%maxl)**2), &
& block_overlap(smap(bas%maxl+1),smap(bas%maxl+1)))

!$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, k) &
!$omp private(r2, vec, vec_diat_trafo, stmp, stmp_diat, dtmpi, dtmpj, hij, rr, inl, img)
!$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)
do iat = 1, mol%nat
izp = mol%id(iat)
is = bas%ish_at(iat)
Expand All @@ -140,19 +141,20 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, &
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)))
call relvec(vec, sqrt(r2), vec_diat_trafo)

! Get the overlap and dipole integrals for the current diatomic pair
block_overlap = 0.0_wp
do ish = 1, bas%nsh_id(izp)
ii = bas%iao_sh(is+ish)
iaosh = smap(ish-1) ! Offset for the block overlap matrix
do jsh = 1, bas%nsh_id(jzp)
jj = bas%iao_sh(js+jsh)
jaosh = smap(jsh-1) ! Offset for the block overlap matrix

call dipole_cgto_diat_scal(bas%cgto(jsh,jzp), bas%cgto(ish,izp), r2, vec, &
bas%intcut, vec_diat_trafo, h0%ksig(izp,jzp), h0%kpi(izp,jzp), h0%kdel(izp,jzp), &
& stmp, stmp_diat, dtmpi)

hij = 0.5_wp * (selfenergy(is+ish) + selfenergy(js+jsh)) &
* h0%hscale(jsh, ish, jzp, izp)
call dipole_cgto(bas%cgto(jsh,jzp), bas%cgto(ish,izp), &
& r2, vec, bas%intcut, stmp, dtmpi)

! Store the overlap and dipole matrix
nao = msao(bas%cgto(jsh, jzp)%ang)
do iao = 1, msao(bas%cgto(ish, izp)%ang)
do jao = 1, nao
Expand All @@ -166,41 +168,70 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, &
+ stmp(ij)

!$omp atomic
overlap_diat(jj+jao, ii+iao) = overlap_diat(jj+jao, ii+iao) &
+ stmp_diat(ij)
block_overlap(jaosh+jao, iaosh+iao) = block_overlap(jaosh+jao, iaosh+iao) &
+ stmp(ij)

do k = 1, 3
!$omp atomic
dpint(k, jj+jao, ii+iao) = dpint(k, jj+jao, ii+iao) &
+ dtmpi(k, ij)
end do

!$omp atomic
hamiltonian(jj+jao, ii+iao) = hamiltonian(jj+jao, ii+iao) &
+ stmp_diat(ij) * hij

if (iat /= jat) then
!$omp atomic
overlap(ii+iao, jj+jao) = overlap(ii+iao, jj+jao) &
+ stmp(ij)

!$omp atomic
overlap_diat(ii+iao, jj+jao) = overlap_diat(ii+iao, jj+jao) &
+ stmp_diat(ij)

do k = 1, 3
!$omp atomic
dpint(k, ii+iao, jj+jao) = dpint(k, ii+iao, jj+jao) &
+ dtmpj(k)
end do

end if
end do
end do
end do
end do

! Diatomic frame transformation and scaling of the overlap
call diat_trafo(block_overlap, vec, h0%ksig(izp,jzp), h0%kpi(izp,jzp), h0%kdel(izp,jzp), &
& bas%nsh_at(jat)-1, bas%nsh_at(iat)-1)

! Setup the Hamiltonian and store the diatomic frame scaled overlap.
do ish = 1, bas%nsh_id(izp)
ii = bas%iao_sh(is+ish)
iaosh = smap(ish-1) ! Offset for the block overlap matrix
do jsh = 1, bas%nsh_id(jzp)
jj = bas%iao_sh(js+jsh)
jaosh = smap(jsh-1) ! Offset for the block overlap matrix

hij = 0.5_wp * h0%hscale(jsh, ish, jzp, izp) * (selfenergy(is+ish) + selfenergy(js+jsh))

nao = msao(bas%cgto(jsh, jzp)%ang)
do iao = 1, msao(bas%cgto(ish, izp)%ang)
do jao = 1, nao
ij = jao + nao*(iao-1)

!$omp atomic
overlap_diat(jj+jao, ii+iao) = overlap_diat(jj+jao, ii+iao) &
+ block_overlap(jaosh+jao, iaosh+iao)

!$omp atomic
hamiltonian(jj+jao, ii+iao) = hamiltonian(jj+jao, ii+iao) &
+ block_overlap(jaosh+jao, iaosh+iao) * hij

if (iat /= jat) then
!$omp atomic
overlap_diat(ii+iao, jj+jao) = overlap_diat(ii+iao, jj+jao) &
+ block_overlap(jaosh+jao, iaosh+iao)

!$omp atomic
hamiltonian(ii+iao, jj+jao) = hamiltonian(ii+iao, jj+jao) &
+ stmp_diat(ij) * hij
+ block_overlap(jaosh+jao, iaosh+iao) * hij
end if
end do
end do

end do
end do

Expand All @@ -223,10 +254,6 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, &
jj = bas%iao_sh(is+jsh)
call dipole_cgto(bas%cgto(jsh,izp), bas%cgto(ish,izp), &
& r2, vec, bas%intcut, stmp, dtmpi)

hij = 0.5_wp * (selfenergy(is+ish) + selfenergy(is+jsh)) &
& * h0%hscale(ish, jsh, izp, izp)

nao = msao(bas%cgto(jsh, izp)%ang)
do iao = 1, msao(bas%cgto(ish, izp)%ang)
do jao = 1, nao
Expand All @@ -239,13 +266,10 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, &

dpint(:, jj+jao, ii+iao) = dpint(:, jj+jao, ii+iao) &
+ dtmpi(:, ij)

hamiltonian(jj+jao, ii+iao) = hamiltonian(jj+jao, ii+iao) &
+ stmp(ij) * hij
end do
end do
end do
! diagonal term (AO(i) == AO(j))
! diagonal term (AO(i) == AO(j)) as on-site off-diagonal the hamiltonian is zero
do iao = 1, msao(bas%cgto(ish, izp)%ang)
hamiltonian(ii+iao, ii+iao) = selfenergy(is + ish)
enddo
Expand Down
Loading

0 comments on commit 79d8f9a

Please sign in to comment.