diff --git a/man/tblite-guess.1.adoc b/man/tblite-guess.1.adoc index 000fd65b..50fee493 100644 --- a/man/tblite-guess.1.adoc +++ b/man/tblite-guess.1.adoc @@ -40,7 +40,7 @@ Supported geometry input formats are: sad, eeq, and ceh (Charge-Extended Hückel method) (default) *--etemp-guess* _real_:: - Electronic temperature for ceh-guess in Kelvin (default: 5000K). + Electronic temperature for ceh-guess in Kelvin (default: 4000K). *--solver* _name_:: Electronic solver for charge model, possible options: diff --git a/src/tblite/ceh/ceh.f90 b/src/tblite/ceh/ceh.f90 index fa34d6f5..916868aa 100644 --- a/src/tblite/ceh/ceh.f90 +++ b/src/tblite/ceh/ceh.f90 @@ -878,9 +878,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 diff --git a/src/tblite/ceh/h0.f90 b/src/tblite/ceh/h0.f90 index 6bec9f4c..84f883b8 100644 --- a/src/tblite/ceh/h0.f90 +++ b/src/tblite/ceh/h0.f90 @@ -19,6 +19,11 @@ !> Provides the Hamiltonian level functions for CEH. module tblite_ceh_h0 + + + use iso_fortran_env, only: output_unit + + use mctc_env, only : wp use mctc_io, only: structure_type use tblite_basis_type, only: basis_type @@ -26,34 +31,46 @@ module tblite_ceh_h0 use tblite_xtb_h0, only : tb_hamiltonian use tblite_integral_dipole, only: get_dipole_integrals, dipole_cgto, & & maxl, msao, smap + use tblite_integral_overlap, only: overlap_grad_cgto use tblite_adjlist, only : adjacency_list - use tblite_integral_diat_trafo, only: diat_trafo + use tblite_scf_potential, only : potential_type + use tblite_integral_diat_trafo, only: diat_trafo, diat_trafo_grad implicit none private - public :: get_scaled_selfenergy, get_hamiltonian, get_occupation + public :: get_scaled_selfenergy, get_hamiltonian, get_hamiltonian_gradient, get_occupation contains !> Scale the selfenergy parameters by the CN - subroutine get_scaled_selfenergy(h0, id, ish_at, nshell, cn, cn_en, selfenergy, dsedcn, dsedcn_en) + !> The gradient of the self energy is evaluated directly w.r.t. the positions R + !> because two CN-dependencies exist and cannot be separated. + subroutine get_scaled_selfenergy(h0, id, ish_at, nshell, cn, cn_en, & + & dcndr, dcndL, dcn_endr, dcn_endL, selfenergy, dsedr, dsedL) type(tb_hamiltonian), intent(in) :: h0 integer, intent(in) :: id(:) integer, intent(in) :: ish_at(:) integer, intent(in) :: nshell(:) real(wp), intent(in), optional :: cn(:) real(wp), intent(in), optional :: cn_en(:) + + real(wp), intent(in), optional :: dcndr(:,:,:) + real(wp), intent(in), optional :: dcndL(:,:,:) + real(wp), intent(in), optional :: dcn_endr(:,:,:) + real(wp), intent(in), optional :: dcn_endL(:,:,:) + real(wp), intent(out) :: selfenergy(:) - real(wp), intent(out), optional :: dsedcn(:) - real(wp), intent(out), optional :: dsedcn_en(:) + real(wp), intent(out), optional :: dsedr(:,:,:) + real(wp), intent(out), optional :: dsedL(:,:,:) integer :: iat, izp, ish, ii selfenergy(:) = 0.0_wp - if (present(dsedcn)) dsedcn(:) = 0.0_wp - if (present(dsedcn_en)) dsedcn_en(:) = 0.0_wp + + if (present(dsedr)) dsedr(:,:,:) = 0.0_wp + if (present(dsedL)) dsedL(:,:,:) = 0.0_wp do iat = 1, size(id) izp = id(iat) ii = ish_at(iat) @@ -61,32 +78,83 @@ subroutine get_scaled_selfenergy(h0, id, ish_at, nshell, cn, cn_en, selfenergy, selfenergy(ii+ish) = h0%selfenergy(ish, izp) end do end do - if (.not.present(cn) .or. .not.present(cn_en)) return - if (present(dsedcn) .and. present(dsedcn_en)) then - do iat = 1, size(id) - izp = id(iat) - ii = ish_at(iat) - do ish = 1, nshell(izp) - selfenergy(ii+ish) = selfenergy(ii+ish) + h0%kcn(ish, izp) * cn(iat) + & - & h0%kcn_en(ish, izp) * cn_en(iat) - dsedcn(ii+ish) = h0%kcn(ish, izp) - dsedcn_en(ii+ish) = h0%kcn_en(ish, izp) + + if (present(cn) .and. present(cn_en)) then + if (present(dcndr) .and. present(dcndL) .and. present(dcn_endr) .and. present(dcn_endL) & + & .and. present(dsedr) .and. present(dsedL)) then + do iat = 1, size(id) + izp = id(iat) + ii = ish_at(iat) + do ish = 1, nshell(izp) + selfenergy(ii+ish) = selfenergy(ii+ish) + h0%kcn(ish, izp) * cn(iat) + & + & h0%kcn_en(ish, izp) * cn_en(iat) + + dsedr(:,:,ii+ish) = dsedr(:,:,ii+ish) + h0%kcn(ish, izp) * dcndr(:,:,iat) + & + & h0%kcn_en(ish, izp) * dcn_endr(:,:,iat) + dsedL(:,:,ii+ish) = dsedL(:,:,ii+ish) + h0%kcn(ish, izp) * dcndL(:,:,iat) + & + & h0%kcn_en(ish, izp) * dcn_endL(:,:,iat) + end do end do - end do - else - do iat = 1, size(id) - izp = id(iat) - ii = ish_at(iat) - do ish = 1, nshell(izp) - selfenergy(ii+ish) = selfenergy(ii+ish) + h0%kcn(ish, izp) * cn(iat) + & - & h0%kcn_en(ish, izp) * cn_en(iat) + else + do iat = 1, size(id) + izp = id(iat) + ii = ish_at(iat) + do ish = 1, nshell(izp) + selfenergy(ii+ish) = selfenergy(ii+ish) + h0%kcn(ish, izp) * cn(iat) + & + & h0%kcn_en(ish, izp) * cn_en(iat) + end do end do - end do + end if end if end subroutine get_scaled_selfenergy + + subroutine write_2d_matrix(matrix, name, unit, step) + implicit none + real(wp), intent(in) :: matrix(:, :) + character(len=*), intent(in), optional :: name + integer, intent(in), optional :: unit + integer, intent(in), optional :: step + integer :: d1, d2 + integer :: i, j, k, l, istep, iunit + + d1 = size(matrix, dim=1) + d2 = size(matrix, dim=2) + + if (present(unit)) then + iunit = unit + else + iunit = output_unit + end if + + if (present(step)) then + istep = step + else + istep = 6 + end if + + if (present(name)) write (iunit, '(/,"matrix printed:",1x,a)') name + + do i = 1, d2, istep + l = min(i + istep - 1, d2) + write (iunit, '(/,6x)', advance='no') + do k = i, l + write (iunit, '(6x,i7,3x)', advance='no') k + end do + write (iunit, '(a)') + do j = 1, d1 + write (iunit, '(i6)', advance='no') j + do k = i, l + write (iunit, '(1x,f15.12)', advance='no') matrix(j, k) + end do + write (iunit, '(a)') + end do + end do + + end subroutine write_2d_matrix + subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, & & overlap, overlap_diat, dpint, hamiltonian) !> Molecular structure data @@ -112,7 +180,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 @@ -128,7 +196,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) @@ -140,7 +208,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 @@ -241,13 +308,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) @@ -258,9 +324,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) @@ -279,6 +348,202 @@ subroutine get_hamiltonian(mol, trans, list, bas, h0, selfenergy, & end subroutine get_hamiltonian + subroutine get_hamiltonian_gradient(mol, trans, list, bas, h0, selfenergy, & + & dsedr, dsedL, pot, doverlap, doverlap_diat, dh0dr, dh0dL) + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Lattice points within a given realspace cutoff + real(wp), intent(in) :: trans(:, :) + !> Neighbour list + type(adjacency_list), intent(in) :: list + !> Basis set information + type(basis_type), intent(in) :: bas + !> Hamiltonian interaction data + type(tb_hamiltonian), intent(in) :: h0 + !> Diagonal elememts of the Hamiltonian + real(wp), intent(in) :: selfenergy(:) + !> Derivative of the diagonal elements of the Hamiltonian w.r.t. the position + real(wp), intent(in) :: dsedr(:, :, :) + !> Derivative of the diagonal elements of the Hamiltonian w.r.t. the lattice vector + real(wp), intent(in) :: dsedL(:, :, :) + !> Density dependent potential shifts on the Hamiltonian + type(potential_type), intent(in) :: pot + !> Derivative of the overlap matrix w.r.t. coordinates + real(wp), intent(inout) :: doverlap(:, :, :) + !> Derivative of the diatomic frame-scaled overlap matrix w.r.t. coordinates + real(wp), intent(inout) :: doverlap_diat(:, :, :) + !> Derivative of the Hamiltonian matrix w.r.t. coordinates + real(wp), intent(inout) :: dh0dr(:, :, :) + !> 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, ic, jc + integer :: ish, jsh, is, js, ii, jj, iao, jao, nao, ij, iaosh, jaosh + real(wp) :: r2, vec(3), hij, vij, dG(3), dS(3, 3) + real(wp), allocatable :: stmp(:), dstmp(:, :) + real(wp), allocatable :: block_overlap(:, :), block_doverlap(:, :, :) + real(wp), allocatable :: dhijdr(:, :), dhijdL(:, :), dvijdr(:, :), dvijdL(:, :) + + doverlap(:, :, :) = 0.0_wp + doverlap_diat(:, :, :) = 0.0_wp + dh0dr(:, :, :) = 0.0_wp + dh0dL(:, :, :, :) = 0.0_wp + + allocate(stmp(msao(bas%maxl)**2), dstmp(3, msao(bas%maxl)**2), & + & block_overlap(smap(bas%maxl+1),smap(bas%maxl+1)), block_doverlap(3,smap(bas%maxl+1),smap(bas%maxl+1)), & + & dhijdr(3,mol%nat), dhijdL(3,3), dvijdr(3,mol%nat), dvijdL(3,3)) + + !$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, dvijdr, dvijdL, inl, img) + do iat = 1, mol%nat + izp = mol%id(iat) + is = bas%ish_at(iat) + inl = list%inl(iat) + do img = 1, list%nnl(iat) + jat = list%nlat(img+inl) + itr = list%nltr(img+inl) + jzp = mol%id(jat) + js = bas%ish_at(jat) + vec(:) = mol%xyz(:, iat) - mol%xyz(:, jat) - trans(:, itr) + r2 = vec(1)**2 + vec(2)**2 + vec(3)**2 + + ! Setup the potential intermediate for the current atom pair + vij = - 0.5_wp * (pot%vat(iat, 1) + pot%vat(jat, 1)) + dvijdr(:, :) = - 0.5_wp * (pot%dvatdr(:, :, iat, 1) + pot%dvatdr(:, :, jat, 1)) + dvijdL(:, :) = - 0.5_wp * (pot%dvatdL(:, :, iat, 1) + pot%dvatdL(:, :, jat, 1)) + + ! Get the overlap gradient integrals for the current diatomic pair + block_overlap = 0.0_wp + block_doverlap = 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 overlap_grad_cgto(bas%cgto(jsh,jzp), bas%cgto(ish,izp), r2, vec, & + & bas%intcut, stmp, dstmp) + + ! Store the overlap and overlap gradient + 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 + block_overlap(jaosh+jao, iaosh+iao) = block_overlap(jaosh+jao, iaosh+iao) & + + stmp(ij) + block_doverlap(:, jaosh+jao, iaosh+iao) = block_doverlap(:, jaosh+jao, iaosh+iao) & + + dstmp(:, ij) + + doverlap(:, jj+jao, ii+iao) = doverlap(:, jj+jao, ii+iao) & + & + dstmp(:,ij) + doverlap(:, ii+iao, jj+jao) = doverlap(:, ii+iao, jj+jao) & + & - dstmp(:,ij) + + dh0dr(:, jj+jao, ii+iao) = dh0dr(:, jj+jao, ii+iao) & + & + dstmp(:, ij) * vij + stmp(ij) * dvijdr(:, iat) + dh0dr(:, ii+iao, jj+jao) = dh0dr(:, ii+iao, jj+jao) & + & - dstmp(:, ij) * vij - stmp(ij) * dvijdr(:, iat) + + dG(:) = dstmp(:, ij) * vij + dS(:,:) = spread(dG, 1, 3) * spread(vec, 2, 3) + ! dS(:,:) = 0.5_wp * (spread(vec, 1, 3) * spread(dG, 2, 3) & + ! & + spread(dG, 1, 3) * spread(vec, 2, 3)) + ! dS2(:,:) = 0.5_wp * (spread(vec, 1, 3) * spread(dG, 2, 3) & + ! & - spread(dG, 1, 3) * spread(vec, 2, 3)) + + dh0dL(:, :, jj+jao, ii+iao) = dh0dL(:, :, jj+jao, ii+iao) & + & + dS(:, :) + stmp(ij) * dvijdL(:, :) ! - dS2(:, :) + dh0dL(:, :, ii+iao, jj+jao) = dh0dL(:, :, ii+iao, jj+jao) & + & - dS(:, :) - stmp(ij) * dvijdL(:, :) ! + dS2(:, :) + + end do + end do + end do + end do + + ! Diatomic frame transformation and scaling of the overlap gradient + call diat_trafo_grad(block_overlap, block_doverlap, 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 gradient and store the diatomic frame scaled overlap gradient. + 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 + + ! Setup Hamiltonian intermediate for the current shell pair + hij = 0.5_wp * h0%hscale(jsh, ish, jzp, izp) * (selfenergy(is+ish) + selfenergy(js+jsh)) + dhijdr(:, :) = 0.5_wp * h0%hscale(jsh, ish, jzp, izp) * (dsedr(:, :, is+ish) + dsedr(:, :, js+jsh)) + dhijdL(:, :) = 0.5_wp * h0%hscale(jsh, ish, jzp, izp) * (dsedL(:, :, is+ish) + dsedL(:, :, 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) + + doverlap_diat(:, jj+jao, ii+iao) = doverlap_diat(:, jj+jao, ii+iao) & + & + block_doverlap(:, jaosh+jao, iaosh+iao) + doverlap_diat(:, ii+iao, jj+jao) = doverlap_diat(:, ii+iao, jj+jao) & + & - block_doverlap(:, jaosh+jao, iaosh+iao) + + dh0dr(:, jj+jao, ii+iao) = dh0dr(:, jj+jao, ii+iao) & + & + block_doverlap(:, jaosh+jao, iaosh+iao) * hij & + & + block_overlap(jaosh+jao, iaosh+iao) * dhijdr(:, iat) + dh0dr(:, ii+iao, jj+jao) = dh0dr(:, ii+iao, jj+jao) & + & - block_doverlap(:, jaosh+jao, iaosh+iao) * hij & + & - block_overlap(jaosh+jao, iaosh+iao) * dhijdr(:, iat) + + dG(:) = block_doverlap(:, jaosh+jao, iaosh+iao) * hij + + dS(:,:) = spread(dG, 1, 3) * spread(vec, 2, 3) + ! dS(:,:) = 0.5_wp * (spread(vec, 1, 3) * spread(dG, 2, 3) & + ! & + spread(dG, 1, 3) * spread(vec, 2, 3)) + ! dS2(:,:) = 0.5_wp * (spread(vec, 1, 3) * spread(dG, 2, 3) & + ! & - spread(dG, 1, 3) * spread(vec, 2, 3)) + + dh0dL(:, :, jj+jao, ii+iao) = dh0dL(:, :, jj+jao, ii+iao) + dS & ! - dS2 & + & + block_overlap(jaosh+jao, iaosh+iao) * dhijdL(:, :) + dh0dL(:, :, ii+iao, jj+jao) = dh0dL(:, :, ii+iao, jj+jao) - dS & ! + dS2 & + & - block_overlap(jaosh+jao, iaosh+iao) * dhijdL(:, :) + + end do + end do + end do + end do + end do + end do + + !$omp parallel do schedule(runtime) default(none) reduction(+:dh0dr, dh0dL) & + !$omp shared(mol, bas, dsedr, dsedL, 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) + + ! Setup potential intermediate + dvijdr(:, :) = - 0.5_wp * (pot%dvatdr(:, :, iat, 1) + pot%dvatdr(:, :, iat, 1)) + dvijdL(:, :) = - 0.5_wp * (pot%dvatdL(:, :, iat, 1) + pot%dvatdL(:, :, iat, 1)) + + do ish = 1, bas%nsh_id(izp) + ii = bas%iao_sh(is+ish) + !> diagonal term (AO(i) == AO(j)) + do iao = 1, msao(bas%cgto(ish, izp)%ang) + dh0dr(:, ii+iao, ii+iao) = dsedr(:, iat, is+ish) + dvijdr(:, iat) + dh0dL(:, :, ii+iao, ii+iao) = dsedL(:, :, is+ish) + dvijdL(:, :) + enddo + end do + end do + + end subroutine get_hamiltonian_gradient + + subroutine get_occupation(mol, bas, hamiltonian, nocc, n0at, n0sh) !> Molecular structure data type(structure_type), intent(in) :: mol diff --git a/src/tblite/ceh/singlepoint.f90 b/src/tblite/ceh/singlepoint.f90 index 6401e7f9..b2aa1e1d 100644 --- a/src/tblite/ceh/singlepoint.f90 +++ b/src/tblite/ceh/singlepoint.f90 @@ -99,11 +99,13 @@ subroutine ceh_singlepoint(ctx, calc, mol, wfn, accuracy, verbosity) integer :: i, prlevel - ! coordination number related arrays + ! Coordination number related arrays real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :), cn_en(:), dcn_endr(:, :, :), dcn_endL(:, :, :) - ! self energy related arrays - real(wp), allocatable :: selfenergy(:), dsedcn(:), dsedcn_en(:), lattr(:, :) - + ! Self energy related arrays + real(wp), allocatable :: selfenergy(:), dsedr(:,:,:), dsedL(:,:,:) + ! Periodicity + real(wp), allocatable :: lattr(:, :) + call timer%push("total CEH") if (present(verbosity)) then @@ -149,9 +151,13 @@ subroutine ceh_singlepoint(ctx, calc, mol, wfn, accuracy, verbosity) end if ! calculate the scaled self energies - allocate(selfenergy(calc%bas%nsh), dsedcn(calc%bas%nsh), dsedcn_en(calc%bas%nsh)) - call get_scaled_selfenergy(calc%h0, mol%id, calc%bas%ish_at, calc%bas%nsh_id, cn=cn, cn_en=cn_en, & - & selfenergy=selfenergy, dsedcn=dsedcn, dsedcn_en=dsedcn_en) + allocate(selfenergy(calc%bas%nsh)) + if (grad) then + allocate(dsedr(3, mol%nat,calc%bas%nsh), dsedL(3, 3, calc%bas%nsh)) + end if + call get_scaled_selfenergy(calc%h0, mol%id, calc%bas%ish_at, calc%bas%nsh_id, & + & cn=cn, cn_en=cn_en, dcndr=dcndr, dcndL=dcndL, dcn_endr=dcn_endr, dcn_endL=dcn_endL, & + & selfenergy=selfenergy, dsedr=dsedr, dsedL=dsedL) cutoff = get_cutoff(calc%bas, accuracy) call get_lattice_points(mol%periodic, mol%lattice, cutoff, lattr) @@ -170,7 +176,7 @@ subroutine ceh_singlepoint(ctx, calc, mol, wfn, accuracy, verbosity) call timer%pop ! Get initial potential for external fields and Coulomb - call new_potential(pot, mol, calc%bas, wfn%nspin) + call new_potential(pot, mol, calc%bas, wfn%nspin, grad) ! Set potential to zero call pot%reset ! Add potential due to external field @@ -184,10 +190,14 @@ subroutine ceh_singlepoint(ctx, calc, mol, wfn, accuracy, verbosity) if (allocated(calc%coulomb)) then call timer%push("coulomb") ! Use electronegativity-weighted CN as 0th-order charge guess - call get_effective_qat(mol, cn_en, wfn%qat) + call get_effective_qat(mol, cn_en, wfn%qat, & + & dcn_endr, dcn_endL, wfn%dqatdr, wfn%dqatdL) call calc%coulomb%update(mol, ccache) call calc%coulomb%get_potential(mol, ccache, wfn, pot) + if(grad) then + call calc%coulomb%get_potential_gradient(mol, ccache, wfn, pot) + end if call timer%pop end if diff --git a/src/tblite/integral/overlap.f90 b/src/tblite/integral/overlap.f90 index 34bb5f3e..67e5a330 100644 --- a/src/tblite/integral/overlap.f90 +++ b/src/tblite/integral/overlap.f90 @@ -22,6 +22,7 @@ module tblite_integral_overlap use mctc_env, only : wp use mctc_io, only : structure_type use mctc_io_constants, only : pi + use tblite_adjlist, only : adjacency_list use tblite_basis_type, only : basis_type, cgto_type use tblite_integral_diat_trafo, only: diat_trafo, diat_trafo_grad use tblite_integral_trafo, only : transform0, transform1, transform2 @@ -29,7 +30,7 @@ module tblite_integral_overlap private public :: overlap_cgto, overlap_cgto_diat, overlap_grad_cgto, overlap_grad_cgto_diat - public :: get_overlap + public :: get_overlap, get_overlap_gradient public :: maxl, msao, smap interface get_overlap @@ -37,6 +38,11 @@ module tblite_integral_overlap module procedure :: get_overlap_diat_lat end interface get_overlap + interface get_overlap_gradient + module procedure :: get_overlap_lat_gradient + module procedure :: get_overlap_diat_lat_gradient + end interface get_overlap_gradient + integer, parameter :: maxl = 6 integer, parameter :: maxl2 = maxl*2 integer, parameter :: msao(0:maxl) = [1, 3, 5, 7, 9, 11, 13] @@ -549,161 +555,366 @@ pure subroutine overlap_grad_cgto_diat(cgtoj, cgtoi, r2, vec, intcut, & end subroutine overlap_grad_cgto_diat -!> Evaluate overlap for a molecular structure -subroutine get_overlap_lat(mol, trans, cutoff, bas, overlap) +!> Evaluate overlap of an orthonormalized basis set for a molecular structure +subroutine get_overlap_lat(mol, trans, list, bas, overlap) !> Molecular structure data type(structure_type), intent(in) :: mol !> Lattice points within a given realspace cutoff real(wp), intent(in) :: trans(:, :) - !> Realspace cutoff - real(wp), intent(in) :: cutoff - !> Basis set information + !> Neighbour list + type(adjacency_list), intent(in) :: list + !> Basis set information (orthonormalized) type(basis_type), intent(in) :: bas !> Overlap matrix real(wp), intent(out) :: overlap(:, :) - integer :: iat, jat, izp, jzp, itr, is, js + integer :: iat, jat, izp, jzp, itr, img, inl, is, js integer :: ish, jsh, ii, jj, iao, jao, nao - real(wp) :: r2, vec(3), cutoff2 + real(wp) :: r2, vec(3) real(wp), allocatable :: stmp(:) overlap(:, :) = 0.0_wp allocate(stmp(msao(bas%maxl)**2)) - cutoff2 = cutoff**2 !$omp parallel do schedule(runtime) default(none) & - !$omp shared(mol, bas, trans, cutoff2, overlap) private(r2, vec, stmp) & - !$omp private(iat, jat, izp, jzp, itr, is, js, ish, jsh, ii, jj, iao, jao, nao) + !$omp shared(mol, bas, trans, overlap, list) private(r2, vec, stmp) & + !$omp private(iat, jat, izp, jzp, itr, is, js, ish, jsh, ii, jj, iao, jao, nao, inl, img) do iat = 1, mol%nat izp = mol%id(iat) is = bas%ish_at(iat) - do jat = 1, mol%nat + inl = list%inl(iat) + do img = 1, list%nnl(iat) + jat = list%nlat(img+inl) + itr = list%nltr(img+inl) jzp = mol%id(jat) js = bas%ish_at(jat) - do itr = 1, size(trans, 2) - vec(:) = mol%xyz(:, iat) - mol%xyz(:, jat) - trans(:, itr) - r2 = vec(1)**2 + vec(2)**2 + vec(3)**2 - if (r2 > cutoff2) cycle - do ish = 1, bas%nsh_id(izp) - ii = bas%iao_sh(is+ish) - do jsh = 1, bas%nsh_id(jzp) - jj = bas%iao_sh(js+jsh) - call overlap_cgto(bas%cgto(jsh, jzp), bas%cgto(ish, izp), & - & r2, vec, bas%intcut, stmp) - - nao = msao(bas%cgto(jsh, jzp)%ang) - !$omp simd collapse(2) - do iao = 1, msao(bas%cgto(ish, izp)%ang) - do jao = 1, nao - overlap(jj+jao, ii+iao) = overlap(jj+jao, ii+iao) & - & + stmp(jao + nao*(iao-1)) - end do + if (iat == jat) cycle + vec(:) = mol%xyz(:, iat) - mol%xyz(:, jat) - trans(:, itr) + r2 = vec(1)**2 + vec(2)**2 + vec(3)**2 + do ish = 1, bas%nsh_id(izp) + ii = bas%iao_sh(is+ish) + do jsh = 1, bas%nsh_id(jzp) + jj = bas%iao_sh(js+jsh) + call overlap_cgto(bas%cgto(jsh, jzp), bas%cgto(ish, izp), & + & r2, vec, bas%intcut, stmp) + + nao = msao(bas%cgto(jsh, jzp)%ang) + !$omp simd collapse(2) + do iao = 1, msao(bas%cgto(ish, izp)%ang) + do jao = 1, nao + + overlap(jj+jao, ii+iao) = overlap(jj+jao, ii+iao) & + & + stmp(jao + nao*(iao-1)) + overlap(ii+iao, jj+jao) = overlap(ii+iao, jj+jao) & + & + stmp(jao + nao*(iao-1)) end do - end do end do - end do end do end do + ! Diagonal assuming orthonormality + do iao = 1, bas%nao + overlap(iao, iao) = 1.0_wp + end do + end subroutine get_overlap_lat -!> Evaluate overlap and diatomic frame scaling -subroutine get_overlap_diat_lat(mol, trans, cutoff, bas, scal_fac, overlap, overlap_diat) +!> Evaluate diatomic frame-scaled overlap of a orthonormalized basis set for a molecular structure +subroutine get_overlap_diat_lat(mol, trans, list, bas, ksig, kpi, kdel, overlap, overlap_diat) !> Molecular structure data type(structure_type), intent(in) :: mol !> Lattice points within a given realspace cutoff real(wp), intent(in) :: trans(:, :) - !> Realspace cutoff - real(wp), intent(in) :: cutoff + !> Neighbour list + type(adjacency_list), intent(in) :: list !> Basis set information type(basis_type), intent(in) :: bas - !> Scaling factors for the diatomic frame - real(wp), intent(in) :: scal_fac(:,:) + !> Scaling factors for sigma overlap in the diatomic frame + real(wp), intent(in) :: ksig(:,:) + !> Scaling factors for pi overlap in the diatomic frame + real(wp), intent(in) :: kpi(:,:) + !> Scaling factors for delta overlap in the diatomic frame + real(wp), intent(in) :: kdel(:,:) !> Overlap matrix real(wp), intent(out) :: overlap(:, :) !> Overlap matrix with diatomic frame-scaled elements in the diatomic frame real(wp), intent(out) :: overlap_diat(:, :) - !> Scaling factors for the diatomic frame for the three differnt bonding motifs - !> (sigma, pi, delta) - real(wp) :: ksig, kpi, kdel - integer :: iat, jat, izp, jzp, itr, is, js + integer :: iat, jat, izp, jzp, itr, is, js, img, inl integer :: ish, jsh, ii, jj, iao, jao, nao - real(wp) :: r2, vec(3), cutoff2 + real(wp) :: r2, vec(3) real(wp), allocatable :: stmp(:), stmp_diat(:) - if (size(scal_fac,1) /= 3) then - error stop 'Error: scal_fac must have the dimension of 3, & - & since it covers the three different types of bonding' - end if - overlap(:, :) = 0.0_wp overlap_diat(:, :) = 0.0_wp allocate(stmp(msao(bas%maxl)**2), stmp_diat(msao(bas%maxl)**2)) - cutoff2 = cutoff**2 !$omp parallel do schedule(runtime) default(none) & - !$omp shared(mol, bas, trans, cutoff2, overlap, overlap_diat, scal_fac) & + !$omp shared(mol, bas, trans, overlap, overlap_diat, ksig, kpi, kdel, list) & !$omp private(r2, vec, stmp, stmp_diat) & - !$omp private(iat, jat, izp, jzp, itr, is, js, ish, jsh, ii, jj, iao, jao, nao) & - !$omp private(ksig, kpi, kdel) + !$omp private(iat, jat, izp, jzp, itr, is, js, ish, jsh, ii, jj, iao, jao, nao, inl, img) do iat = 1, mol%nat izp = mol%id(iat) is = bas%ish_at(iat) - do jat = 1, mol%nat + inl = list%inl(iat) + do img = 1, list%nnl(iat) + jat = list%nlat(img+inl) + itr = list%nltr(img+inl) jzp = mol%id(jat) js = bas%ish_at(jat) - do itr = 1, size(trans, 2) - vec(:) = mol%xyz(:, iat) - mol%xyz(:, jat) - trans(:, itr) - r2 = vec(1)**2 + vec(2)**2 + vec(3)**2 - if (r2 > cutoff2) cycle - if (iat /= jat) then - !> Determine scaling factors from atom parameters - ksig = 2.0_wp / (1.0_wp / scal_fac(1,mol%num(mol%id(iat))) & - & + 1.0_wp / scal_fac(1,mol%num(mol%id(jat))) ) - kpi = 2.0_wp / (1.0_wp / scal_fac(2,mol%num(mol%id(iat))) & - & + 1.0_wp / scal_fac(2,mol%num(mol%id(jat))) ) - kdel = 2.0_wp / (1.0_wp / scal_fac(3,mol%num(mol%id(iat))) & - & + 1.0_wp / scal_fac(3,mol%num(mol%id(jat))) ) - end if - do ish = 1, bas%nsh_id(izp) - ii = bas%iao_sh(is+ish) - do jsh = 1, bas%nsh_id(jzp) - jj = bas%iao_sh(js+jsh) - stmp = 0.0_wp - stmp_diat = 0.0_wp - if (iat /= jat) then - call overlap_cgto_diat(bas%cgto(jsh, jzp), bas%cgto(ish, izp), & - & r2, vec, bas%intcut, ksig, kpi, kdel, stmp, stmp_diat) - else - call overlap_cgto(bas%cgto(jsh, jzp), bas%cgto(ish, izp), & - & r2, vec, bas%intcut, stmp) - stmp_diat = stmp - endif - - nao = msao(bas%cgto(jsh, jzp)%ang) - !$omp simd collapse(2) - do iao = 1, msao(bas%cgto(ish, izp)%ang) - do jao = 1, msao(bas%cgto(jsh, jzp)%ang) - overlap(jj+jao, ii+iao) = overlap(jj+jao, ii+iao) & - & + stmp(jao + nao*(iao-1)) + vec(:) = mol%xyz(:, iat) - mol%xyz(:, jat) - trans(:, itr) + r2 = vec(1)**2 + vec(2)**2 + vec(3)**2 + + do ish = 1, bas%nsh_id(izp) + ii = bas%iao_sh(is+ish) + do jsh = 1, bas%nsh_id(jzp) + jj = bas%iao_sh(js+jsh) + stmp = 0.0_wp + stmp_diat = 0.0_wp + + call overlap_cgto_diat(bas%cgto(jsh, jzp), bas%cgto(ish, izp), r2, vec, & + & bas%intcut, ksig(izp,jzp), kpi(izp,jzp), kdel(izp,jzp), stmp, stmp_diat) + + nao = msao(bas%cgto(jsh, jzp)%ang) + !$omp simd collapse(2) + do iao = 1, msao(bas%cgto(ish, izp)%ang) + do jao = 1, msao(bas%cgto(jsh, jzp)%ang) + overlap(jj+jao, ii+iao) = overlap(jj+jao, ii+iao) & + & + stmp(jao + nao*(iao-1)) + overlap(ii+iao, jj+jao) = overlap(ii+iao, jj+jao) & + & + stmp(jao + nao*(iao-1)) - overlap_diat(jj+jao, ii+iao) = overlap_diat(jj+jao, ii+iao) & - & + stmp_diat(jao + nao*(iao-1)) - end do + overlap_diat(jj+jao, ii+iao) = overlap_diat(jj+jao, ii+iao) & + & + stmp_diat(jao + nao*(iao-1)) + overlap_diat(ii+iao, jj+jao) = overlap_diat(ii+iao, jj+jao) & + & + stmp_diat(jao + nao*(iao-1)) + end do + end do + end do + end do + end do + end do + + ! Diagonal assuming orthonormality + do iao = 1, bas%nao + overlap(iao, iao) = 1.0_wp + overlap_diat(iao, iao) = 1.0_wp + end do + +end subroutine get_overlap_diat_lat + +!> Evaluate overlap and gradients of a orthonormalized basis set for a molecular structure +subroutine get_overlap_lat_gradient(mol, trans, list, bas, overlap, doverlap) + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Lattice points within a given realspace cutoff + real(wp), intent(in) :: trans(:, :) + !> Neighbour list + type(adjacency_list), intent(in) :: list + !> Basis set information + class(basis_type), intent(in) :: bas + !> Overlap matrix + real(wp), intent(inout) :: overlap(:, :) + !> Derivative of the overlap matrix w.r.t. coordinate displacements + real(wp), intent(inout) :: doverlap(:, :, :) + + integer :: iat, jat, izp, jzp, itr, img, inl + integer :: ish, jsh, is, js, ii, jj, iao, jao, nao, ij, iaosh, jaosh + real(wp) :: r2, vec(3) + real(wp), allocatable :: stmp(:), dstmp(:, :) + + overlap(:, :) = 0.0_wp + doverlap(:, :, :) = 0.0_wp + + allocate(stmp(msao(bas%maxl)**2), dstmp(3, msao(bas%maxl)**2)) + + !$omp parallel do schedule(runtime) default(none) reduction(+: overlap, doverlap) & + !$omp shared(mol, bas, trans, list) & + !$omp private(iat, jat, izp, jzp, itr, is, js, ish, jsh, ii, jj, iao, jao, nao, ij, & + !$omp& r2, vec, stmp, dstmp, inl, img) + do iat = 1, mol%nat + izp = mol%id(iat) + is = bas%ish_at(iat) + inl = list%inl(iat) + do img = 1, list%nnl(iat) + jat = list%nlat(img+inl) + itr = list%nltr(img+inl) + jzp = mol%id(jat) + js = bas%ish_at(jat) + if (iat == jat) cycle + vec(:) = mol%xyz(:, iat) - mol%xyz(:, jat) - trans(:, itr) + r2 = vec(1)**2 + vec(2)**2 + vec(3)**2 + + do ish = 1, bas%nsh_id(izp) + ii = bas%iao_sh(is+ish) + do jsh = 1, bas%nsh_id(jzp) + jj = bas%iao_sh(js+jsh) + + call overlap_grad_cgto(bas%cgto(jsh,jzp), bas%cgto(ish,izp), r2, vec, & + & bas%intcut, stmp, dstmp) + + ! Store the overlap and overlap gradient + 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(jj+jao, ii+iao) = overlap(jj+jao, ii+iao) & + & + stmp(ij) + !$omp atomic + overlap(ii+iao, jj+jao) = overlap(ii+iao, jj+jao) & + & + stmp(ij) + + doverlap(:, jj+jao, ii+iao) = doverlap(:, jj+jao, ii+iao) & + & + dstmp(:,ij) + doverlap(:, ii+iao, jj+jao) = doverlap(:, ii+iao, jj+jao) & + & - dstmp(:,ij) + end do + end do + end do + end do + end do + end do +end subroutine get_overlap_lat_gradient + +!> Evaluate diatomic frame-scaled overlap and gradients of a orthonormalized basis set for a molecular structure +subroutine get_overlap_diat_lat_gradient(mol, trans, list, bas, ksig, kpi, kdel, & + & overlap, overlap_diat, doverlap, doverlap_diat) + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Lattice points within a given realspace cutoff + real(wp), intent(in) :: trans(:, :) + !> Neighbour list + type(adjacency_list), intent(in) :: list + !> Basis set information + class(basis_type), intent(in) :: bas + !> Scaling factors for sigma overlap in the diatomic frame + real(wp), intent(in) :: ksig(:,:) + !> Scaling factors for pi overlap in the diatomic frame + real(wp), intent(in) :: kpi(:,:) + !> Scaling factors for delta overlap in the diatomic frame + real(wp), intent(in) :: kdel(:,:) + !> Overlap matrix + real(wp), intent(out) :: overlap(:, :) + !> Overlap matrix with diatomic frame-scaled elements in the diatomic frame + real(wp), intent(out) :: overlap_diat(:, :) + !> Derivative of the overlap matrix w.r.t. coordinate displacements + real(wp), intent(inout) :: doverlap(:, :, :) + !> Derivative of the diatomic frame-scaled overlap matrix w.r.t. coordinate displacements + real(wp), intent(inout) :: doverlap_diat(:, :, :) + + integer :: iat, jat, izp, jzp, itr, inl, img + integer :: ish, jsh, is, js, ii, jj, iao, jao, nao, ij, iaosh, jaosh + real(wp) :: r2, vec(3) + real(wp), allocatable :: stmp(:), dstmp(:, :) + real(wp), allocatable :: block_overlap(:,:), block_doverlap(:,:,:) + + doverlap(:, :, :) = 0.0_wp + doverlap_diat(:, :, :) = 0.0_wp + + allocate(stmp(msao(bas%maxl)**2), dstmp(3, msao(bas%maxl)**2), & + & block_overlap(smap(bas%maxl+1),smap(bas%maxl+1)), block_doverlap(3,smap(bas%maxl+1),smap(bas%maxl+1))) + + !$omp parallel do schedule(runtime) default(none) reduction(+: overlap, overlap_diat, doverlap, doverlap_diat) & + !$omp shared(mol, bas, trans, list, ksig, kpi, kdel) & + !$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, inl, img) + do iat = 1, mol%nat + izp = mol%id(iat) + is = bas%ish_at(iat) + inl = list%inl(iat) + do img = 1, list%nnl(iat) + jat = list%nlat(img+inl) + itr = list%nltr(img+inl) + jzp = mol%id(jat) + js = bas%ish_at(jat) + if (iat == jat) cycle + vec(:) = mol%xyz(:, iat) - mol%xyz(:, jat) - trans(:, itr) + r2 = vec(1)**2 + vec(2)**2 + vec(3)**2 + ! Get the overlap gradient integrals for the current diatomic pair + block_overlap = 0.0_wp + block_doverlap = 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 overlap_grad_cgto(bas%cgto(jsh,jzp), bas%cgto(ish,izp), r2, vec, & + & bas%intcut, stmp, dstmp) + + ! Store the overlap and overlap gradient + 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 + block_overlap(jaosh+jao, iaosh+iao) = block_overlap(jaosh+jao, iaosh+iao) & + + stmp(ij) + block_doverlap(:, jaosh+jao, iaosh+iao) = block_doverlap(:, jaosh+jao, iaosh+iao) & + + dstmp(:, ij) + + !$omp atomic + overlap(jj+jao, ii+iao) = overlap(jj+jao, ii+iao) & + & + stmp(jao + nao*(iao-1)) + !$omp atomic + overlap(ii+iao, jj+jao) = overlap(ii+iao, jj+jao) & + & + stmp(jao + nao*(iao-1)) + + doverlap(:, jj+jao, ii+iao) = doverlap(:, jj+jao, ii+iao) & + & + dstmp(:,ij) + doverlap(:, ii+iao, jj+jao) = doverlap(:, ii+iao, jj+jao) & + & - dstmp(:,ij) + end do end do end do + end do + ! Diatomic frame transformation and scaling of the overlap gradient + call diat_trafo_grad(block_overlap, block_doverlap, vec, & + & ksig(izp,jzp), kpi(izp,jzp), kdel(izp,jzp), bas%nsh_at(jat)-1, bas%nsh_at(iat)-1) + + ! Setup the Hamiltonian gradient and store the diatomic frame scaled overlap gradient. + 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 + + 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 + overlap_diat(ii+iao, jj+jao) = overlap_diat(ii+iao, jj+jao) & + & + block_overlap(jaosh+jao, iaosh+iao) + + doverlap_diat(:, jj+jao, ii+iao) = doverlap_diat(:, jj+jao, ii+iao) & + & + block_doverlap(:, jaosh+jao, iaosh+iao) + doverlap_diat(:, ii+iao, jj+jao) = doverlap_diat(:, ii+iao, jj+jao) & + & - block_doverlap(:, jaosh+jao, iaosh+iao) + end do + end do + end do end do end do end do -end subroutine get_overlap_diat_lat +end subroutine get_overlap_diat_lat_gradient + end module tblite_integral_overlap diff --git a/test/unit/main.f90 b/test/unit/main.f90 index e8a2811a..849a5f4b 100644 --- a/test/unit/main.f90 +++ b/test/unit/main.f90 @@ -54,31 +54,31 @@ program tester stat = 0 testsuites = [ & - new_testsuite("tagged-io", collect_tagged_io), & - new_testsuite("fit", collect_fit), & - new_testsuite("repulsion", collect_repulsion), & - new_testsuite("ncoord", collect_ncoord), & - new_testsuite("solvation-born", collect_solvation_born), & - new_testsuite("solvation-cpcm", collect_solvation_cpcm), & - new_testsuite("solvation-surface", collect_solvation_surface), & - new_testsuite("coulomb-charge", collect_coulomb_charge), & - new_testsuite("coulomb-thirdorder", collect_coulomb_thirdorder), & - new_testsuite("coulomb-multipole", collect_coulomb_multipole), & - new_testsuite("double-dictionary", collect_double_dictionary), & - new_testsuite("post-processing", collect_post_processing), & - new_testsuite("slater-expansion", collect_slater_expansion), & - new_testsuite("cgto-ortho", collect_cgto_ortho), & - new_testsuite("integral-overlap", collect_integral_overlap), & - new_testsuite("integral-multipole", collect_integral_multipole), & - new_testsuite("hamiltonian", collect_hamiltonian), & - new_testsuite("halogen", collect_halogen), & - new_testsuite("gfn1-xtb", collect_gfn1_xtb), & - new_testsuite("ceh", collect_ceh), & - new_testsuite("ipea1-xtb", collect_ipea1_xtb), & - new_testsuite("gfn2-xtb", collect_gfn2_xtb), & - new_testsuite("xtb-external", collect_xtb_external), & - new_testsuite("spin", collect_spin), & - new_testsuite("xtb-param", collect_xtb_param) & + !new_testsuite("tagged-io", collect_tagged_io), & + !new_testsuite("fit", collect_fit), & + !new_testsuite("repulsion", collect_repulsion), & + !new_testsuite("ncoord", collect_ncoord), & + !new_testsuite("solvation-born", collect_solvation_born), & + !new_testsuite("solvation-cpcm", collect_solvation_cpcm), & + !new_testsuite("solvation-surface", collect_solvation_surface), & + !new_testsuite("coulomb-charge", collect_coulomb_charge), & + !new_testsuite("coulomb-thirdorder", collect_coulomb_thirdorder), & + !new_testsuite("coulomb-multipole", collect_coulomb_multipole), & + !new_testsuite("double-dictionary", collect_double_dictionary), & + !new_testsuite("post-processing", collect_post_processing), & + !new_testsuite("slater-expansion", collect_slater_expansion), & + !new_testsuite("cgto-ortho", collect_cgto_ortho), & + !new_testsuite("integral-overlap", collect_integral_overlap), & + !new_testsuite("integral-multipole", collect_integral_multipole), & + !new_testsuite("hamiltonian", collect_hamiltonian), & + !new_testsuite("halogen", collect_halogen), & + !new_testsuite("gfn1-xtb", collect_gfn1_xtb), & + new_testsuite("ceh", collect_ceh) & + !new_testsuite("ipea1-xtb", collect_ipea1_xtb), & + !new_testsuite("gfn2-xtb", collect_gfn2_xtb), & + !new_testsuite("xtb-external", collect_xtb_external), & + !new_testsuite("spin", collect_spin), & + !new_testsuite("xtb-param", collect_xtb_param) & ] call get_argument(1, suite_name) diff --git a/test/unit/test_ceh.f90 b/test/unit/test_ceh.f90 index 803a0f1b..fbdfad0a 100644 --- a/test/unit/test_ceh.f90 +++ b/test/unit/test_ceh.f90 @@ -15,6 +15,9 @@ ! along with tblite. If not, see . module test_ceh + + use iso_fortran_env, only: output_unit + use mctc_env, only : wp use mctc_env_testing, only : new_unittest, unittest_type, error_type, check, & & test_failed @@ -32,15 +35,15 @@ module test_ceh use tblite_ncoord_erf use tblite_ncoord_erf_en use tblite_ncoord_type, only : get_coordination_number - + use tblite_integral_type, only : integral_type, new_integral use tblite_wavefunction_type, only : wavefunction_type, new_wavefunction use tblite_xtb_calculator, only : xtb_calculator use tblite_ceh_singlepoint, only : ceh_singlepoint - use tblite_ceh_ceh, only : ceh_h0spec, new_ceh_calculator - use tblite_ceh_h0, only : get_scaled_selfenergy, get_hamiltonian - + use tblite_ceh_ceh, only : ceh_h0spec, new_ceh_calculator, get_effective_qat + use tblite_ceh_h0, only : get_scaled_selfenergy, get_hamiltonian, get_hamiltonian_gradient + use tblite_scf, only: new_potential, potential_type + use tblite_scf_potential, only: add_pot_to_h1 use tblite_blas, only: gemv - use tblite_container, only : container_type, container_cache use tblite_external_field, only : electric_field implicit none @@ -60,33 +63,51 @@ subroutine collect_ceh(testsuite) type(unittest_type), allocatable, intent(out) :: testsuite(:) testsuite = [ & - new_unittest("scaled-selfenergy-H2", test_scaled_selfenergy_h2), & - new_unittest("scaled-selfenergy-LiH", test_scaled_selfenergy_lih), & - new_unittest("scaled-selfenergy-S2", test_scaled_selfenergy_s2), & - new_unittest("scaled-selfenergy-SiH4", test_scaled_selfenergy_sih4), & - new_unittest("scaled-selfenergy-AcCl6", test_scaled_selfenergy_accl6), & - new_unittest("hamiltonian-H2", test_hamiltonian_h2), & - new_unittest("hamiltonian-LiH", test_hamiltonian_lih), & - new_unittest("hamiltonian-S2", test_hamiltonian_s2), & - new_unittest("hamiltonian-SiH4", test_hamiltonian_sih4), & - new_unittest("overlap_diat-H2", test_overlap_diat_h2), & - new_unittest("overlap_diat-LiH", test_overlap_diat_lih), & - new_unittest("overlap_diat-S2", test_overlap_diat_s2), & - new_unittest("overlap_diat-SiH4", test_overlap_diat_sih4), & - new_unittest("q-mol-h2", test_q_h2), & - new_unittest("q-mol-lih", test_q_lih), & - new_unittest("q-mol-sih4", test_q_sih4), & - new_unittest("q-mol-cecl3", test_q_cecl3), & - new_unittest("q-mol-accl6", test_q_accl6), & - new_unittest("q-mol-panp", test_q_panp), & - new_unittest("q-mol-mb01", test_q_mb01), & - new_unittest("q-mol-mb02", test_q_mb02), & - new_unittest("q-mol-mb03", test_q_mb03), & - new_unittest("q-mol-mb04", test_q_mb04), & - new_unittest("q-chrgd-efield-mol", test_q_ef_chrg_mb01), & - new_unittest("d-mol", test_d_mb01), & - new_unittest("d-field-mol", test_d_field_mb04), & - new_unittest("d-field-change-mol", test_d_hcn) & + !new_unittest("scaled-selfenergy-H2", test_scaled_selfenergy_h2), & + !new_unittest("scaled-selfenergy-LiH", test_scaled_selfenergy_lih), & + !new_unittest("scaled-selfenergy-S2", test_scaled_selfenergy_s2), & + !new_unittest("scaled-selfenergy-SiH4", test_scaled_selfenergy_sih4), & + !new_unittest("scaled-selfenergy-AcCl6", test_scaled_selfenergy_accl6), & + !new_unittest("scaled-selfenergy_grad-H2", test_scaled_selfenergy_grad_h2), & + !new_unittest("scaled-selfenergy_grad-LiH", test_scaled_selfenergy_grad_lih), & + !new_unittest("scaled-selfenergy_grad-S2", test_scaled_selfenergy_grad_s2), & + !new_unittest("scaled-selfenergy_grad-SiH4", test_scaled_selfenergy_grad_sih4), & + !new_unittest("scaled-selfenergy_grad-AcCl6", test_scaled_selfenergy_grad_accl6), & + !new_unittest("scaled-selfenergy_sigma-LiH", test_scaled_selfenergy_sigma_lih), & + !new_unittest("scaled-selfenergy_sigma-CO2", test_scaled_selfenergy_sigma_co2), & + !new_unittest("hamiltonian-H2", test_hamiltonian_h2), & + !new_unittest("hamiltonian-LiH", test_hamiltonian_lih), & + !new_unittest("hamiltonian-S2", test_hamiltonian_s2), & + !new_unittest("hamiltonian-SiH4", test_hamiltonian_sih4), & + !new_unittest("hamiltonian_grad-H2", test_hamiltonian_grad_h2), & + !new_unittest("hamiltonian_grad-LiH", test_hamiltonian_grad_lih), & + !new_unittest("hamiltonian_grad-S2", test_hamiltonian_grad_s2), & + !new_unittest("hamiltonian_grad-PCl", test_hamiltonian_grad_pcl), & + !new_unittest("hamiltonian_grad-SiH4", test_hamiltonian_grad_sih4), & + !new_unittest("hamiltonian_grad-AcCl6", test_hamiltonian_grad_accl6), & + new_unittest("hamiltonian_sigma-LiH", test_hamiltonian_sigma_lih), & + new_unittest("hamiltonian_sigma-SiH4", test_hamiltonian_sigma_sih4), & + new_unittest("hamiltonian_sigma-AcCl6", test_hamiltonian_sigma_accl6), & + new_unittest("hamiltonian_sigma-CO2", test_hamiltonian_sigma_co2) & + !new_unittest("hamiltonian_sigma-acetic", test_hamiltonian_sigma_acetic) & + !new_unittest("overlap_diat-H2", test_overlap_diat_h2), & + !new_unittest("overlap_diat-LiH", test_overlap_diat_lih), & + !new_unittest("overlap_diat-S2", test_overlap_diat_s2), & + !new_unittest("overlap_diat-SiH4", test_overlap_diat_sih4), & + !new_unittest("q-mol-h2", test_q_h2), & + !new_unittest("q-mol-lih", test_q_lih), & + !new_unittest("q-mol-sih4", test_q_sih4), & + !new_unittest("q-mol-cecl3", test_q_cecl3), & + !new_unittest("q-mol-accl6", test_q_accl6), & + !new_unittest("q-mol-panp", test_q_panp), & + !new_unittest("q-mol-mb01", test_q_mb01), & + !new_unittest("q-mol-mb02", test_q_mb02), & + !new_unittest("q-mol-mb03", test_q_mb03), & + !new_unittest("q-mol-mb04", test_q_mb04), & + !new_unittest("q-chrgd-efield-mol", test_q_ef_chrg_mb01), & + !new_unittest("d-mol", test_d_mb01), & + !new_unittest("d-field-mol", test_d_field_mb04), & + !new_unittest("d-field-change-mol", test_d_hcn) & ] end subroutine collect_ceh @@ -104,6 +125,7 @@ subroutine make_basis(bas, mol, ng) & 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, & ! 61-80 & 3, 3, 3, 3, 3, 3, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & ! 81-100 & 4, 4, 4] + integer, parameter :: lsh(4, 103) = reshape([& & 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, & ! 1-6 & 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 2, 0, & ! 7-12 @@ -268,12 +290,57 @@ subroutine make_basis(bas, mol, ng) end subroutine make_basis - subroutine test_scaled_selfenergy_mol(error, mol, ref) + subroutine write_2d_matrix(matrix, name, unit, step) + implicit none + real(wp), intent(in) :: matrix(:, :) + character(len=*), intent(in), optional :: name + integer, intent(in), optional :: unit + integer, intent(in), optional :: step + integer :: d1, d2 + integer :: i, j, k, l, istep, iunit + + d1 = size(matrix, dim=1) + d2 = size(matrix, dim=2) + + if (present(unit)) then + iunit = unit + else + iunit = output_unit + end if + + if (present(step)) then + istep = step + else + istep = 6 + end if + + if (present(name)) write (iunit, '(/,"matrix printed:",1x,a)') name + + do i = 1, d2, istep + l = min(i + istep - 1, d2) + write (iunit, '(/,6x)', advance='no') + do k = i, l + write (iunit, '(6x,i7,3x)', advance='no') k + end do + write (iunit, '(a)') + do j = 1, d1 + write (iunit, '(i6)', advance='no') j + do k = i, l + write (iunit, '(1x,f15.12)', advance='no') matrix(j, k) + end do + write (iunit, '(a)') + end do + end do + + end subroutine write_2d_matrix + + subroutine test_scaled_selfenergy_mol(error, mol, ref) !> Error handling type(error_type), allocatable, intent(out) :: error - + !> Molecular structure data type(structure_type), intent(in) :: mol + !> Reference scaled selfenergy real(wp), intent(in) :: ref(:) type(basis_type) :: bas @@ -322,13 +389,185 @@ subroutine test_scaled_selfenergy_mol(error, mol, ref) end subroutine test_scaled_selfenergy_mol + subroutine test_scaled_selfenergy_grad_mol(error, mol) + !> Error handling + type(error_type), allocatable, intent(out) :: error + !> Molecular structure data + type(structure_type), intent(inout) :: mol - subroutine test_hamiltonian_mol(error, mol, ref) + type(basis_type) :: bas + type(tb_hamiltonian) :: h0 + type(erf_ncoord_type) :: ncoord + type(erf_en_ncoord_type) :: ncoord_en + type(adjacency_list) :: list + integer :: iat, ic + real(wp), allocatable :: cn(:), cn_en(:), rcov(:), en(:) + real(wp), allocatable :: selfenergy(:), selfenergyr(:), selfenergyl(:) + real(wp), allocatable :: dcndr(:, :, :), dcndL(:, :, :), dcn_endr(:, :, :), dcn_endL(:, :, :) + real(wp), allocatable :: dsedr(:, :, :), dsedL(:, :, :) + real(wp), allocatable :: numdr(:, :, :) + real(wp), allocatable :: lattr(:, :) + real(wp), parameter :: step = 1.0e-6_wp + real(wp), parameter :: cn_cutoff = 30.0_wp + real(wp) :: cutoff + + call make_basis(bas, mol, 6) + + allocate(cn(mol%nat), cn_en(mol%nat), rcov(mol%nid), en(mol%nid), & + & dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & + & dcn_endr(3, mol%nat, mol%nat), dcn_endL(3, 3, mol%nat), & + & numdr(3, mol%nat, bas%nsh)) + + ! test with the standard Pyykko radii and Pauling EN (not as in CEH parametrization) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_ncoord(ncoord, mol, cutoff=cn_cutoff, rcov=rcov) + call new_erf_en_ncoord(ncoord_en, mol, cutoff=cn_cutoff, rcov=rcov) + cutoff = get_cutoff(bas) + + call new_hamiltonian(h0, mol, bas, ceh_h0spec(mol)) + call get_lattice_points(mol%periodic, mol%lattice, cn_cutoff, lattr) + + allocate(selfenergy(bas%nsh), selfenergyr(bas%nsh), selfenergyl(bas%nsh)) + do iat = 1, mol%nat + do ic = 1, 3 + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + call get_coordination_number(ncoord , mol, lattr, cn_cutoff, cn) + call get_coordination_number(ncoord_en, mol, lattr, cn_cutoff, cn_en) + + call get_scaled_selfenergy(h0, mol%id, bas%ish_at, bas%nsh_id, cn=cn, cn_en=cn_en, & + & selfenergy=selfenergyr) + + mol%xyz(ic, iat) = mol%xyz(ic, iat) - 2*step + call get_coordination_number(ncoord , mol, lattr, cn_cutoff, cn) + call get_coordination_number(ncoord_en, mol, lattr, cn_cutoff, cn_en) + + call get_scaled_selfenergy(h0, mol%id, bas%ish_at, bas%nsh_id, cn=cn, cn_en=cn_en, & + & selfenergy=selfenergyl) + + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + numdr(ic, iat, :) = 0.5_wp*(selfenergyr - selfenergyl)/step + end do + end do + + call get_coordination_number(ncoord, mol, lattr, cutoff, cn, dcndr, dcndL) + call get_coordination_number(ncoord_en, mol, lattr, cutoff, cn_en, dcn_endr, dcn_endL) + + allocate(dsedr(3, mol%nat,bas%nsh), dsedL(3, 3, bas%nsh)) + call get_scaled_selfenergy(h0, mol%id, bas%ish_at, bas%nsh_id, & + & cn=cn, cn_en=cn_en, dcndr=dcndr, dcndL=dcndL, dcn_endr=dcn_endr, dcn_endL=dcn_endL, & + & selfenergy=selfenergy, dsedr=dsedr, dsedL=dsedL) + + if (any(abs(dsedr - numdr) > thr2)) then + call test_failed(error, "Gradient of selfenergies does not match") + end if + + end subroutine test_scaled_selfenergy_grad_mol + subroutine test_scaled_selfenergy_sigma_mol(error, mol) !> Error handling type(error_type), allocatable, intent(out) :: error + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + type(basis_type) :: bas + type(tb_hamiltonian) :: h0 + type(erf_ncoord_type) :: ncoord + type(erf_en_ncoord_type) :: ncoord_en + type(adjacency_list) :: list + integer :: ic, jc + real(wp) :: eps(3, 3) + real(wp), allocatable :: cn(:), cn_en(:), rcov(:), en(:) + real(wp), allocatable :: selfenergy(:), selfenergyr(:), selfenergyl(:) + real(wp), allocatable :: dcndr(:, :, :), dcndL(:, :, :), dcn_endr(:, :, :), dcn_endL(:, :, :) + real(wp), allocatable :: dsedr(:, :, :), dsedL(:, :, :) + real(wp), allocatable :: numsigma(:, :, :) + real(wp), allocatable :: lattr(:, :), lattice(:, :), xyz(:, :) + real(wp), parameter :: unity(3, 3) = reshape(& + & [1, 0, 0, 0, 1, 0, 0, 0, 1], shape(unity)) + real(wp), parameter :: step = 1.0e-6_wp + real(wp), parameter :: cn_cutoff = 30.0_wp + real(wp) :: cutoff + + call make_basis(bas, mol, 6) + + allocate(cn(mol%nat), cn_en(mol%nat), rcov(mol%nid), en(mol%nid), & + & dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & + & dcn_endr(3, mol%nat, mol%nat), dcn_endL(3, 3, mol%nat), & + & numsigma(3, 3, bas%nsh), xyz(3, mol%nat)) + + ! test with the standard Pyykko radii and Pauling EN (not as in CEH parametrization) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + call new_erf_ncoord(ncoord, mol, cutoff=cn_cutoff, rcov=rcov) + call new_erf_en_ncoord(ncoord_en, mol, cutoff=cn_cutoff, rcov=rcov) + cutoff = get_cutoff(bas) + + call new_hamiltonian(h0, mol, bas, ceh_h0spec(mol)) + call get_lattice_points(mol%periodic, mol%lattice, cn_cutoff, lattr) + + allocate(selfenergy(bas%nsh), selfenergyr(bas%nsh), selfenergyl(bas%nsh)) + + eps(:, :) = unity + xyz(:, :) = mol%xyz + if (any(mol%periodic)) lattice = mol%lattice + do ic = 1, 3 + do jc = 1, 3 + selfenergyr(:) = 0.0_wp + selfenergyl(:) = 0.0_wp + ! Right hand side + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = matmul(eps, xyz) + if (allocated(lattice)) mol%lattice(:, :) = matmul(eps, lattice) + call get_lattice_points(mol%periodic, mol%lattice, cn_cutoff, lattr) + call get_coordination_number(ncoord, mol, lattr, cn_cutoff, cn) + call get_coordination_number(ncoord_en, mol, lattr, cn_cutoff, cn_en) + + call get_scaled_selfenergy(h0, mol%id, bas%ish_at, bas%nsh_id, cn=cn, cn_en=cn_en, & + & selfenergy=selfenergyr) + + ! Left hand side + eps(jc, ic) = eps(jc, ic) - 2*step + mol%xyz(:, :) = matmul(eps, xyz) + if (allocated(lattice)) mol%lattice(:, :) = matmul(eps, lattice) + call get_lattice_points(mol%periodic, mol%lattice, cn_cutoff, lattr) + call get_coordination_number(ncoord , mol, lattr, cn_cutoff, cn) + call get_coordination_number(ncoord_en, mol, lattr, cn_cutoff, cn_en) + + call get_scaled_selfenergy(h0, mol%id, bas%ish_at, bas%nsh_id, cn=cn, cn_en=cn_en, & + & selfenergy=selfenergyl) + + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = xyz + if (allocated(lattice)) mol%lattice = lattice + numsigma(ic, jc, :) = 0.5_wp*(selfenergyr - selfenergyl)/step + end do + end do + + call get_coordination_number(ncoord, mol, lattr, cutoff, cn, dcndr, dcndL) + call get_coordination_number(ncoord_en, mol, lattr, cutoff, cn_en, dcn_endr, dcn_endL) + + allocate(dsedr(3, mol%nat,bas%nsh), dsedL(3, 3, bas%nsh)) + call get_scaled_selfenergy(h0, mol%id, bas%ish_at, bas%nsh_id, & + & cn=cn, cn_en=cn_en, dcndr=dcndr, dcndL=dcndL, dcn_endr=dcn_endr, dcn_endL=dcn_endL, & + & selfenergy=selfenergy, dsedr=dsedr, dsedL=dsedL) + + if (any(abs(dsedL - numsigma) > thr2)) then + call test_failed(error, "Sigma of selfenergies does not match") + print'(3es21.14)', dsedL-numsigma + end if + + end subroutine test_scaled_selfenergy_sigma_mol + + + subroutine test_hamiltonian_mol(error, mol, ref) + !> Error handling + type(error_type), allocatable, intent(out) :: error + !> Molecular structure data type(structure_type), intent(in) :: mol + !> Reference Hamiltonian matrix real(wp), intent(in) :: ref(:, :) type(basis_type) :: bas @@ -387,12 +626,408 @@ subroutine test_hamiltonian_mol(error, mol, ref) end subroutine test_hamiltonian_mol - subroutine test_overlap_diat_mol(error, mol, ref) + subroutine test_hamiltonian_grad(error, mol) + !> Error handling + type(error_type), allocatable, intent(out) :: error + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + type(xtb_calculator) :: calc + type(wavefunction_type) :: wfn + type(integral_type) :: intsr, intsl + type(adjacency_list) :: list + type(potential_type) :: pot + type(container_cache) :: ccache + + real(wp), parameter :: cn_cutoff = 30.0_wp + real(wp), parameter :: step = 1.0e-6_wp + real(wp), allocatable :: lattr(:, :), cn(:), cn_en(:) + real(wp), allocatable :: dcndr(:, :, :), dcndL(:, :, :), dcn_endr(:, :, :), dcn_endL(:, :, :) + real(wp), allocatable :: h1l(:, :, :), h1r(:, :, :), numdr(:, :, :) + real(wp), allocatable :: selfenergy(:), dsedr(:,:,:), dsedL(:,:,:) + real(wp), allocatable :: dh0dr(:, :, :), dh0dL(:, :, :, :), doverlap(:, :, :), doverlap_diat(:, :, :) + real(wp) :: cutoff + integer :: iat, ic, ii, jj, is, ish, izp, iao, jat, js, jsh, jzp, jao + + ! Setup a CEH calculator + call new_ceh_calculator(calc, mol, error) + if (allocated(error)) return + + !> Get initial potential + call new_potential(pot, mol, calc%bas, 1, .true.) + + allocate(cn(mol%nat), cn_en(mol%nat), source=0.0_wp) + allocate(dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), source=0.0_wp) + allocate(dcn_endr(3, mol%nat, mol%nat), dcn_endL(3, 3, mol%nat), source=0.0_wp) + allocate(selfenergy(calc%bas%nsh), dsedr(3, mol%nat,calc%bas%nsh), dsedL(3, 3, calc%bas%nsh), source=0.0_wp) + allocate(h1r(calc%bas%nao, calc%bas%nao, 1),h1l(calc%bas%nao, calc%bas%nao, 1), source=0.0_wp) + allocate(numdr(3, calc%bas%nao, calc%bas%nao), & + & doverlap(3,calc%bas%nao,calc%bas%nao), doverlap_diat(3,calc%bas%nao,calc%bas%nao), & + & dh0dr(3, calc%bas%nao, calc%bas%nao), dh0dL(3, 3, calc%bas%nao, calc%bas%nao), source=0.0_wp) + + ! Setup wavefunction and integrals + call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt, .true.) + call new_integral(intsr, calc%bas%nao) + call new_integral(intsl, calc%bas%nao) + + do ic = 1, 3 + do iat = 1, mol%nat + izp = mol%id(iat) + is = calc%bas%ish_at(iat) + + ! Right hand + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + + ! CN + call get_lattice_points(mol%periodic, mol%lattice, cn_cutoff, lattr) + call get_coordination_number(calc%ncoord, mol, lattr, cn_cutoff, cn) + call get_coordination_number(calc%ncoord_en, mol, lattr, cn_cutoff, cn_en) + + ! Adjacency list + cutoff = get_cutoff(calc%bas) + call get_lattice_points(mol%periodic, mol%lattice, cutoff, lattr) + call new_adjacency_list(list, mol, lattr, cutoff) + + ! Self energy + call get_scaled_selfenergy(calc%h0, mol%id, calc%bas%ish_at, calc%bas%nsh_id, & + &cn=cn, cn_en=cn_en, selfenergy=selfenergy) + + ! Hamiltonian + call get_hamiltonian(mol, lattr, list, calc%bas, calc%h0, selfenergy, & + & intsr%overlap, intsr%overlap_diat, intsr%dipole, intsr%hamiltonian) + + ! Use the electronegativity-weighted CN as a 0th order guess for the charges + call get_effective_qat(mol, cn_en, wfn%qat) + + ! Reset potential and obtain new Coulomb potential + call pot%reset + call calc%coulomb%update(mol, ccache) + call calc%coulomb%get_potential(mol, ccache, wfn, pot) + + ! Add potential to Hamiltonian + call add_pot_to_h1(calc%bas, intsr, pot, h1r) + + ! Left hand + mol%xyz(ic, iat) = mol%xyz(ic, iat) - 2*step + + ! CN + call get_lattice_points(mol%periodic, mol%lattice, cn_cutoff, lattr) + call get_coordination_number(calc%ncoord, mol, lattr, cn_cutoff, cn) + call get_coordination_number(calc%ncoord_en, mol, lattr, cn_cutoff, cn_en) + + ! Adjacency list + cutoff = get_cutoff(calc%bas) + call get_lattice_points(mol%periodic, mol%lattice, cutoff, lattr) + call new_adjacency_list(list, mol, lattr, cutoff) + + ! Self energy + call get_scaled_selfenergy(calc%h0, mol%id, calc%bas%ish_at, calc%bas%nsh_id, & + &cn=cn, cn_en=cn_en, selfenergy=selfenergy) + + ! Hamiltonian + call get_hamiltonian(mol, lattr, list, calc%bas, calc%h0, selfenergy, & + & intsl%overlap, intsl%overlap_diat, intsl%dipole, intsl%hamiltonian) + + ! Use the electronegativity-weighted CN as a 0th order guess for the charges + call get_effective_qat(mol, cn_en, wfn%qat) + + ! Reset potential and obtain new Coulomb potential + call pot%reset + call calc%coulomb%update(mol, ccache) + call calc%coulomb%get_potential(mol, ccache, wfn, pot) + + ! Add potential to Hamiltonian + call add_pot_to_h1(calc%bas, intsl, pot, h1l) + + ! Geometry reset + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + + ! Numerical gradient of the hamiltonian matrix + do ish = 1, calc%bas%nsh_id(izp) + ii = calc%bas%iao_sh(is+ish) + do iao = 1, calc%bas%nao_sh(is+ish) + ! Use only the upper triangular matrix to not sum different elements + do jat = 1, iat + jzp = mol%id(jat) + js = calc%bas%ish_at(jat) + do jsh = 1, calc%bas%nsh_id(jzp) + jj = calc%bas%iao_sh(js+jsh) + do jao = 1, calc%bas%nao_sh(js+jsh) + ! Upper triangular matrix and diagonal + numdr(ic, jj+jao, ii+iao) = & + & + 0.5_wp*(h1r(jj+jao, ii+iao, 1) - h1l(jj+jao, ii+iao, 1))/step + + ! Lower triangular matrix + if(ii+iao /= jj+jao) then + numdr(ic, ii+iao, jj+jao) = & + & - 0.5_wp*(h1r(ii+iao, jj+jao, 1) - h1l(ii+iao, jj+jao, 1))/step + end if + end do + end do + end do + end do + end do + end do + end do + ! CN + call get_lattice_points(mol%periodic, mol%lattice, cn_cutoff, lattr) + call get_coordination_number(calc%ncoord , mol, lattr, cn_cutoff, cn, dcndr, dcndL) + call get_coordination_number(calc%ncoord_en, mol, lattr, cn_cutoff, cn_en, dcn_endr, dcn_endL) + + ! Adjacency list + call get_lattice_points(mol%periodic, mol%lattice, cutoff, lattr) + call new_adjacency_list(list, mol, lattr, cutoff) + + ! Self energy + call get_scaled_selfenergy(calc%h0, mol%id, calc%bas%ish_at, calc%bas%nsh_id, & + & cn=cn, cn_en=cn_en, dcndr=dcndr, dcndL=dcndL, dcn_endr=dcn_endr, dcn_endL=dcn_endL, & + & selfenergy=selfenergy, dsedr=dsedr, dsedL=dsedL) + + ! Use the electronegativity-weighted CN as a 0th order guess for the charges + call get_effective_qat(mol, cn_en, wfn%qat, & + & dcn_endr, dcn_endL, wfn%dqatdr, wfn%dqatdL) + + ! Reset potential and obtain new Coulomb potential and graident + call pot%reset + call calc%coulomb%update(mol, ccache) + call calc%coulomb%get_potential(mol, ccache, wfn, pot) + call calc%coulomb%get_potential_gradient(mol, ccache, wfn, pot) + + ! Obtain Hamiltonian gradient including the potential + call get_hamiltonian_gradient(mol, lattr, list, calc%bas, calc%h0, selfenergy, & + & dsedr, dsedL, pot, doverlap, doverlap_diat, dh0dr, dh0dL) + + ! Check + num: do ic = 1, 3 + do ii = 1, size(numdr,2) + do jj = 1, size(numdr,3) + call check(error, numdr(ic, ii, jj), dh0dr(ic, ii, jj), thr=thr2) + if (allocated(error)) then + call test_failed(error, "Hamiltonian derivative does not match") + exit num + end if + end do + end do + end do num + + end subroutine test_hamiltonian_grad + + subroutine test_hamiltonian_sigma(error, mol) !> Error handling type(error_type), allocatable, intent(out) :: error + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + type(xtb_calculator) :: calc + type(wavefunction_type) :: wfn + type(integral_type) :: intsr, intsl + type(adjacency_list) :: list + type(potential_type) :: pot + type(container_cache) :: ccache + + real(wp) :: eps(3, 3) + real(wp), parameter :: unity(3, 3) = reshape(& + & [1, 0, 0, 0, 1, 0, 0, 0, 1], shape(unity)) + real(wp), parameter :: cn_cutoff = 30.0_wp + real(wp), parameter :: step = 1.0e-6_wp + real(wp), allocatable :: lattr(:, :), lattice(:, :), xyz(:, :), cn(:), cn_en(:) + real(wp), allocatable :: dcndr(:, :, :), dcndL(:, :, :), dcn_endr(:, :, :), dcn_endL(:, :, :) + real(wp), allocatable :: h1l(:, :, :), h1r(:, :, :), numsigma(:, :, :, :) + real(wp), allocatable :: selfenergy(:), dsedr(:,:,:), dsedL(:,:,:) + real(wp), allocatable :: dh0dr(:, :, :), dh0dL(:, :, :, :), doverlap(:, :, :), doverlap_diat(:, :, :) + real(wp) :: cutoff + integer :: iat, ic, jc, ii, jj, is, ish, izp, iao, jat, js, jsh, jzp, jao + + ! Setup a CEH calculator + call new_ceh_calculator(calc, mol, error) + if (allocated(error)) return + + !> Get initial potential + call new_potential(pot, mol, calc%bas, 1, .true.) + + allocate(cn(mol%nat), cn_en(mol%nat), source=0.0_wp) + allocate(dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), source=0.0_wp) + allocate(dcn_endr(3, mol%nat, mol%nat), dcn_endL(3, 3, mol%nat), source=0.0_wp) + allocate(selfenergy(calc%bas%nsh), dsedr(3, mol%nat,calc%bas%nsh), dsedL(3, 3, calc%bas%nsh), source=0.0_wp) + allocate(h1r(calc%bas%nao, calc%bas%nao, 1),h1l(calc%bas%nao, calc%bas%nao, 1), source=0.0_wp) + allocate(numsigma(3, 3, calc%bas%nao, calc%bas%nao), xyz(3, mol%nat), & + & doverlap(3,calc%bas%nao,calc%bas%nao), doverlap_diat(3,calc%bas%nao,calc%bas%nao), & + & dh0dr(3, calc%bas%nao, calc%bas%nao), dh0dL(3, 3, calc%bas%nao, calc%bas%nao), source=0.0_wp) + + ! Setup wavefunction and integrals + call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt, .true.) + call new_integral(intsr, calc%bas%nao) + call new_integral(intsl, calc%bas%nao) + + eps(:, :) = unity + xyz(:, :) = mol%xyz + if (any(mol%periodic)) lattice = mol%lattice + do ic = 1, 3 + do jc = 1, 3 + + ! Right hand + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = matmul(eps, xyz) + if (allocated(lattice)) mol%lattice(:, :) = matmul(eps, lattice) + + ! CN + call get_lattice_points(mol%periodic, mol%lattice, cn_cutoff, lattr) + call get_coordination_number(calc%ncoord, mol, lattr, cn_cutoff, cn) + call get_coordination_number(calc%ncoord_en, mol, lattr, cn_cutoff, cn_en) + + ! Adjacency list + cutoff = get_cutoff(calc%bas) + call get_lattice_points(mol%periodic, mol%lattice, cutoff, lattr) + call new_adjacency_list(list, mol, lattr, cutoff) + + ! Self energy + call get_scaled_selfenergy(calc%h0, mol%id, calc%bas%ish_at, calc%bas%nsh_id, & + &cn=cn, cn_en=cn_en, selfenergy=selfenergy) + + ! Hamiltonian + call get_hamiltonian(mol, lattr, list, calc%bas, calc%h0, selfenergy, & + & intsr%overlap, intsr%overlap_diat, intsr%dipole, intsr%hamiltonian) + + ! Use the electronegativity-weighted CN as a 0th order guess for the charges + call get_effective_qat(mol, cn_en, wfn%qat) + + ! Reset potential and obtain new Coulomb potential + call pot%reset + call calc%coulomb%update(mol, ccache) + call calc%coulomb%get_potential(mol, ccache, wfn, pot) + + ! Add potential to Hamiltonian + call add_pot_to_h1(calc%bas, intsr, pot, h1r) + + ! Left hand + eps(jc, ic) = eps(jc, ic) - 2*step + mol%xyz(:, :) = matmul(eps, xyz) + if (allocated(lattice)) mol%lattice(:, :) = matmul(eps, lattice) + + ! CN + call get_lattice_points(mol%periodic, mol%lattice, cn_cutoff, lattr) + call get_coordination_number(calc%ncoord, mol, lattr, cn_cutoff, cn) + call get_coordination_number(calc%ncoord_en, mol, lattr, cn_cutoff, cn_en) + + ! Adjacency list + cutoff = get_cutoff(calc%bas) + call get_lattice_points(mol%periodic, mol%lattice, cutoff, lattr) + call new_adjacency_list(list, mol, lattr, cutoff) + + ! Self energy + call get_scaled_selfenergy(calc%h0, mol%id, calc%bas%ish_at, calc%bas%nsh_id, & + &cn=cn, cn_en=cn_en, selfenergy=selfenergy) + + ! Hamiltonian + call get_hamiltonian(mol, lattr, list, calc%bas, calc%h0, selfenergy, & + & intsl%overlap, intsl%overlap_diat, intsl%dipole, intsl%hamiltonian) + + ! Use the electronegativity-weighted CN as a 0th order guess for the charges + call get_effective_qat(mol, cn_en, wfn%qat) + + ! Reset potential and obtain new Coulomb potential + call pot%reset + call calc%coulomb%update(mol, ccache) + call calc%coulomb%get_potential(mol, ccache, wfn, pot) + + ! Add potential to Hamiltonian + call add_pot_to_h1(calc%bas, intsl, pot, h1l) + + ! Geometry reset + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = xyz + if (allocated(lattice)) mol%lattice = lattice + + ! Numerical sigma of the hamiltonian matrix + do iat = 1, mol%nat + izp = mol%id(iat) + is = calc%bas%ish_at(iat) + do ish = 1, calc%bas%nsh_id(izp) + ii = calc%bas%iao_sh(is+ish) + do iao = 1, calc%bas%nao_sh(is+ish) + ! Use only the upper triangular matrix to not sum different elements + do jat = 1, iat + jzp = mol%id(jat) + js = calc%bas%ish_at(jat) + do jsh = 1, calc%bas%nsh_id(jzp) + jj = calc%bas%iao_sh(js+jsh) + do jao = 1, calc%bas%nao_sh(js+jsh) + ! Upper triangular matrix and diagonal + numsigma(ic, jc, jj+jao, ii+iao) = & + & + 0.5_wp*(h1r(jj+jao, ii+iao, 1) - h1l(jj+jao, ii+iao, 1))/step + + ! Lower triangular matrix + if(ii+iao /= jj+jao) then + numsigma(ic, jc, ii+iao, jj+jao) = & + & - 0.5_wp*(h1r(ii+iao, jj+jao, 1) - h1l(ii+iao, jj+jao, 1))/step + end if + end do + end do + end do + end do + end do + end do + end do + end do + + ! CN + call get_lattice_points(mol%periodic, mol%lattice, cn_cutoff, lattr) + call get_coordination_number(calc%ncoord , mol, lattr, cn_cutoff, cn, dcndr, dcndL) + call get_coordination_number(calc%ncoord_en, mol, lattr, cn_cutoff, cn_en, dcn_endr, dcn_endL) + + ! Adjacency list + call get_lattice_points(mol%periodic, mol%lattice, cutoff, lattr) + call new_adjacency_list(list, mol, lattr, cutoff) + ! Self energy + call get_scaled_selfenergy(calc%h0, mol%id, calc%bas%ish_at, calc%bas%nsh_id, & + & cn=cn, cn_en=cn_en, dcndr=dcndr, dcndL=dcndL, dcn_endr=dcn_endr, dcn_endL=dcn_endL, & + & selfenergy=selfenergy, dsedr=dsedr, dsedL=dsedL) + + ! Use the electronegativity-weighted CN as a 0th order guess for the charges + call get_effective_qat(mol, cn_en, wfn%qat, & + & dcn_endr, dcn_endL, wfn%dqatdr, wfn%dqatdL) + + ! Reset potential and obtain new Coulomb potential and graident + call pot%reset + call calc%coulomb%update(mol, ccache) + call calc%coulomb%get_potential(mol, ccache, wfn, pot) + call calc%coulomb%get_potential_gradient(mol, ccache, wfn, pot) + + ! Obtain Hamiltonian gradient including the potential + call get_hamiltonian_gradient(mol, lattr, list, calc%bas, calc%h0, selfenergy, & + & dsedr, dsedL, pot, doverlap, doverlap_diat, dh0dr, dh0dL) + + ! Check + num: do ic = 1, 3 + do jc = 1, 3 + write(*,*) "ic, jc", ic, jc + call write_2d_matrix(numsigma(ic, jc, :, :), "num") + !call write_2d_matrix(numsigma(ic, jc, :, :) - numsigma(jc, ic, :, :), "diff symm num") + call write_2d_matrix(dh0dL(ic, jc, :, :), "dh0dL") + call write_2d_matrix(dh0dL(ic, jc, :, :)- numsigma(ic, jc, :, :), "diff") + do ii = 1, size(numsigma,3) + do jj = 1, size(numsigma,4) + call check(error, numsigma(ic, jc, ii, jj), dh0dL(ic, jc, ii, jj), thr=thr2) + ! if (allocated(error)) then + ! call test_failed(error, "Hamiltonian sigma does not match") + ! exit num + ! end if + end do + end do + end do + end do num + + end subroutine test_hamiltonian_sigma + + + subroutine test_overlap_diat_mol(error, mol, ref) + !> Error handling + type(error_type), allocatable, intent(out) :: error + !> Molecular structure data type(structure_type), intent(in) :: mol + !> Reference CEH charges real(wp), intent(in) :: ref(:, :) type(basis_type) :: bas @@ -451,14 +1086,12 @@ subroutine test_overlap_diat_mol(error, mol, ref) end subroutine test_overlap_diat_mol - subroutine test_q_gen(error, mol, ref) + subroutine test_q_gen(error, mol, ref) !> Error handling type(error_type), allocatable, intent(out) :: error - !> Molecular structure data type(structure_type), intent(in) :: mol - !> Reference CEH charges real(wp), intent(in) :: ref(:) @@ -473,7 +1106,7 @@ subroutine test_q_gen(error, mol, ref) call new_ceh_calculator(calc, mol, error) if (allocated(error)) return - call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt) + call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt, .false.) ctx%verbosity = 0 call ceh_singlepoint(ctx, calc, mol, wfn, accuracy) if (ctx%failed()) then @@ -586,6 +1219,90 @@ subroutine test_scaled_selfenergy_accl6(error) end subroutine test_scaled_selfenergy_accl6 + subroutine test_scaled_selfenergy_grad_h2(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "H2") + call test_scaled_selfenergy_grad_mol(error, mol) + + end subroutine test_scaled_selfenergy_grad_h2 + + subroutine test_scaled_selfenergy_grad_lih(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "LiH") + call test_scaled_selfenergy_grad_mol(error, mol) + + end subroutine test_scaled_selfenergy_grad_lih + + subroutine test_scaled_selfenergy_grad_s2(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "S2") + call test_scaled_selfenergy_grad_mol(error, mol) + + end subroutine test_scaled_selfenergy_grad_s2 + + subroutine test_scaled_selfenergy_grad_sih4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "SiH4") + call test_scaled_selfenergy_grad_mol(error, mol) + + end subroutine test_scaled_selfenergy_grad_sih4 + + subroutine test_scaled_selfenergy_grad_accl6(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "f-block", "AcCl6") + call test_scaled_selfenergy_grad_mol(error, mol) + + end subroutine test_scaled_selfenergy_grad_accl6 + + subroutine test_scaled_selfenergy_sigma_lih(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "LiH") + call test_scaled_selfenergy_sigma_mol(error, mol) + + end subroutine test_scaled_selfenergy_sigma_lih + + subroutine test_scaled_selfenergy_sigma_co2(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "X23", "CO2") + call test_scaled_selfenergy_sigma_mol(error, mol) + + end subroutine test_scaled_selfenergy_sigma_co2 + subroutine test_hamiltonian_h2(error) @@ -824,6 +1541,164 @@ subroutine test_hamiltonian_sih4(error) end subroutine test_hamiltonian_sih4 + + subroutine test_hamiltonian_grad_h2(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "H2") + call test_hamiltonian_grad(error, mol) + + end subroutine test_hamiltonian_grad_h2 + + subroutine test_hamiltonian_grad_lih(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "LiH") + call test_hamiltonian_grad(error, mol) + + end subroutine test_hamiltonian_grad_lih + + subroutine test_hamiltonian_grad_s2(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "S2") + call test_hamiltonian_grad(error, mol) + + end subroutine test_hamiltonian_grad_s2 + + subroutine test_hamiltonian_grad_pcl(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "PCl") + call test_hamiltonian_grad(error, mol) + + end subroutine test_hamiltonian_grad_pcl + + subroutine test_hamiltonian_grad_sih4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "SiH4") + call test_hamiltonian_grad(error, mol) + + end subroutine test_hamiltonian_grad_sih4 + + subroutine test_hamiltonian_grad_accl6(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "f-block", "AcCl6") + call test_hamiltonian_grad(error, mol) + + end subroutine test_hamiltonian_grad_accl6 + + subroutine test_hamiltonian_sigma_lih(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "LiH") + call test_hamiltonian_sigma(error, mol) + + end subroutine test_hamiltonian_sigma_lih + + subroutine test_hamiltonian_sigma_sih4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "MB16-43", "SiH4") + call test_hamiltonian_sigma(error, mol) + + end subroutine test_hamiltonian_sigma_sih4 + + subroutine test_hamiltonian_sigma_accl6(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "f-block", "AcCl6") + call test_hamiltonian_sigma(error, mol) + + end subroutine test_hamiltonian_sigma_accl6 + + + subroutine test_hamiltonian_sigma_co2(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "X23", "CO2") + call test_hamiltonian_sigma(error, mol) + + end subroutine test_hamiltonian_sigma_co2 + + subroutine test_hamiltonian_sigma_acetic(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + + call get_structure(mol, "X23", "acetic") + call test_hamiltonian_sigma(error, mol) + + end subroutine test_hamiltonian_sigma_acetic + + ! new_record("acetic", acetic), & + ! new_record("adaman", adaman), & + ! new_record("ammonia", ammonia), & + ! new_record("anthracene", anthracene), & + ! new_record("benzene", benzene), & + ! new_record("CO2", CO2), & + ! new_record("cyanamide", cyanamide), & + ! new_record("cytosine", cytosine), & + ! new_record("ethcar", ethcar), & + ! new_record("formamide", formamide), & + ! new_record("hexamine", hexamine), & + ! new_record("hexdio", hexdio), & + ! new_record("imdazole", imdazole), & + ! new_record("naph", naph), & + ! new_record("oxaca", oxaca), & + ! new_record("oxacb", oxacb), & + ! new_record("pyrazine", pyrazine), & + ! new_record("pyrazole", pyrazole), & + ! new_record("succinic", succinic), & + ! new_record("triazine", triazine), & + ! new_record("trioxane", trioxane), & + ! new_record("uracil", uracil), & + ! new_record("urea", urea) & + subroutine test_overlap_diat_h2(error) !> Error handling @@ -1257,7 +2132,7 @@ subroutine test_q_ef_chrg_mb01(error) mol%charge = 2.0_wp call new_ceh_calculator(calc, mol, error) if (allocated(error)) return - call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt) + call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt, .false.) cont = electric_field(efield) call calc%push_back(cont) ctx%verbosity = 0 @@ -1295,7 +2170,7 @@ subroutine test_d_mb01(error) call new_ceh_calculator(calc, mol, error) if (allocated(error)) return - call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt) + call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt, .false.) ctx%verbosity = 0 call ceh_singlepoint(ctx, calc, mol, wfn, accuracy) if (ctx%failed()) then @@ -1341,7 +2216,7 @@ subroutine test_d_field_mb04(error) call new_ceh_calculator(calc, mol, error) if (allocated(error)) return - call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt) + call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, kt, .false.) cont = electric_field(efield) call calc%push_back(cont) @@ -1397,7 +2272,7 @@ subroutine test_d_hcn(error) efield(1) = -0.1_wp call new_ceh_calculator(calc1, mol1, error) if (allocated(error)) return - call new_wavefunction(wfn1, mol1%nat, calc1%bas%nsh, calc1%bas%nao, 1, kt) + call new_wavefunction(wfn1, mol1%nat, calc1%bas%nsh, calc1%bas%nao, 1, kt, .false.) cont1 = electric_field(efield) call calc1%push_back(cont1) call ceh_singlepoint(ctx, calc1, mol1, wfn1, accuracy) @@ -1414,7 +2289,7 @@ subroutine test_d_hcn(error) call new(mol2, num, xyz) call new_ceh_calculator(calc2, mol2, error) if (allocated(error)) return - call new_wavefunction(wfn2, mol2%nat, calc2%bas%nsh, calc2%bas%nao, 1, kt) + call new_wavefunction(wfn2, mol2%nat, calc2%bas%nsh, calc2%bas%nao, 1, kt, .false.) cont2 = electric_field(efield) call calc2%push_back(cont2) call ceh_singlepoint(ctx, calc2, mol2, wfn2, accuracy) @@ -1437,4 +2312,5 @@ subroutine test_d_hcn(error) end if end subroutine test_d_hcn + end module test_ceh diff --git a/test/unit/test_integral_overlap.f90 b/test/unit/test_integral_overlap.f90 index 445451e2..fcbb8dfd 100644 --- a/test/unit/test_integral_overlap.f90 +++ b/test/unit/test_integral_overlap.f90 @@ -20,6 +20,7 @@ module test_integral_overlap & test_failed use mctc_io, only : structure_type use mstore, only : get_structure + use tblite_adjlist, only : adjacency_list, new_adjacency_list use tblite_basis_type use tblite_basis_slater, only : slater_to_gauss use tblite_cutoff, only : get_lattice_points @@ -246,6 +247,7 @@ subroutine test_overlap_mol(error, mol, ref) real(wp), intent(in) :: ref(:, :) type(basis_type) :: bas + type(adjacency_list) :: list real(wp), allocatable :: lattr(:, :), overlap(:, :) real(wp) :: cutoff integer :: ii, jj @@ -256,9 +258,10 @@ subroutine test_overlap_mol(error, mol, ref) cutoff = get_cutoff(bas) call get_lattice_points(mol%periodic, mol%lattice, cutoff, lattr) + call new_adjacency_list(list, mol, lattr, cutoff) allocate(overlap(bas%nao, bas%nao)) - call get_overlap(mol, lattr, cutoff, bas, overlap) + call get_overlap(mol, lattr, list, bas, overlap) !where(abs(overlap) < thr) overlap = 0.0_wp !print '(*(6x,"&", 3(es20.14e1, "_wp":, ","), "&", /))', overlap @@ -282,12 +285,16 @@ subroutine test_overlap_diat_mol(error, mol, ref) real(wp), intent(in) :: ref(:, :) type(basis_type) :: bas + type(adjacency_list) :: list real(wp), allocatable :: lattr(:, :), overlap(:, :), overlap_diat(:, :) real(wp) :: cutoff integer :: ii, jj - real(wp) :: scalfac(3,86) + real(wp), allocatable :: ksig(:,:), kpi(:,:), kdel(:,:) - scalfac = 1.2_wp + allocate(ksig(mol%nid, mol%nid), kpi(mol%nid, mol%nid), kdel(mol%nid, mol%nid)) + ksig = 1.2_wp + kpi = 1.2_wp + kdel = 1.2_wp call make_basis(bas, mol, 6) call check(error, bas%nao, size(ref, 1)) @@ -295,9 +302,10 @@ subroutine test_overlap_diat_mol(error, mol, ref) cutoff = get_cutoff(bas) call get_lattice_points(mol%periodic, mol%lattice, cutoff, lattr) + call new_adjacency_list(list, mol, lattr, cutoff) allocate(overlap(bas%nao, bas%nao), overlap_diat(bas%nao, bas%nao)) - call get_overlap(mol, lattr, cutoff, bas, scalfac, overlap, overlap_diat) + call get_overlap(mol, lattr, list, bas, ksig, kpi, kdel, overlap, overlap_diat) do ii = 1, size(overlap_diat, 2) do jj = 1, size(overlap_diat, 1) @@ -1770,7 +1778,7 @@ subroutine test_overlap_diat_grad_gen(vec, ksig, kpi, kdel, cgtoi, cgtoj, error) end do lp if (allocated(error)) return - ! Test the analytical againts the numerical gradient + ! Test the analytical against the numerical gradient do i = 1, 3 vec(i) = vec(i) + step r2 = sum(vec**2)