! This file is part of s-dftd3. ! SPDX-Identifier: LGPL-3.0-or-later ! ! s-dftd3 is free software: you can redistribute it and/or modify it under ! the terms of the GNU Lesser General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! ! s-dftd3 is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public License ! along with s-dftd3. If not, see <https://www.gnu.org/licenses/>. module dftd3_gcp use dftd3_cutoff, only : realspace_cutoff, get_lattice_points use dftd3_gcp_param, only : gcp_param, get_gcp_param use dftd3_partition, only : work_partition, owns_pair use mctc_env, only : wp use mctc_io, only : structure_type implicit none private public :: gcp_param, get_gcp_param, get_geometric_counterpoise public :: get_geometric_counterpoise_hessian !> Geometric counterpoise correction interface get_geometric_counterpoise module procedure get_geometric_counterpoise module procedure get_geometric_counterpoise_atomic end interface get_geometric_counterpoise !> Exponents of the short-range basis incompleteness correction real(wp), parameter :: rexp_base = 0.75_wp, zexp_base = 1.5_wp real(wp), parameter :: rexp_srb = -1.0_wp, zexp_srb = 0.5_wp contains !> Geometric counterpoise correction subroutine get_geometric_counterpoise(mol, param, cutoff, energy, gradient, sigma, partition) !> Molecular structure data class(structure_type), intent(in) :: mol !> Geometric counterpoise parameters type(gcp_param), intent(in) :: param !> Realspace cutoffs type(realspace_cutoff), intent(in) :: cutoff !> Counter-poise energy real(wp), intent(inout) :: energy !> Counter-poise gradient real(wp), intent(inout), optional :: gradient(:, :) !> Counter-poise virial real(wp), intent(inout), optional :: sigma(:, :) !> Work partition of the atom pairs, absent selects the complete work type(work_partition), intent(in), optional :: partition real(wp), allocatable :: energies(:) allocate(energies(mol%nat), source=0.0_wp) call get_geometric_counterpoise_atomic(mol, param, cutoff, energies, gradient, sigma, & & partition) energy = energy + sum(energies) end subroutine get_geometric_counterpoise !> Geometric counterpoise correction with atom-resolved energies subroutine get_geometric_counterpoise_atomic(mol, param, cutoff, energies, gradient, sigma, & & partition) !> Molecular structure data class(structure_type), intent(in) :: mol !> Geometric counterpoise parameters type(gcp_param), intent(in) :: param !> Realspace cutoffs type(realspace_cutoff), intent(in) :: cutoff !> Dispersion energy real(wp), intent(inout) :: energies(:) !> Dispersion gradient real(wp), intent(inout), optional :: gradient(:, :) !> Dispersion virial real(wp), intent(inout), optional :: sigma(:, :) !> Work partition of the atom pairs, absent selects the complete work type(work_partition), intent(in), optional :: partition real(wp), allocatable :: lattr(:, :) logical :: grad grad = present(gradient) .and. present(sigma) if (allocated(param%emiss) .and. allocated(param%slater) .and. allocated(param%xv)) then call get_lattice_points(mol%periodic, mol%lattice, cutoff%gcp, lattr) if (grad) then call gcp_deriv(mol, lattr, cutoff%gcp, param%zeff, param%emiss, param%slater, & & param%xv, param%rvdw, param%sigma, param%alpha, param%beta, param%damp, & & param%dmp_scal, param%dmp_exp, energies, gradient, sigma, partition) else call gcp_energy(mol, lattr, cutoff%gcp, param%zeff, param%emiss, param%slater, & & param%xv, param%rvdw, param%sigma, param%alpha, param%beta, param%damp, & & param%dmp_scal, param%dmp_exp, energies, partition) end if end if if (param%srb .or. param%base) then call get_lattice_points(mol%periodic, mol%lattice, cutoff%srb, lattr) end if if (param%srb) then if (grad) then call srb_deriv(mol, lattr, cutoff%srb, mol%num, param%rvdw_srb, param%rscal, param%qscal, & & rexp_srb, zexp_srb, energies, gradient, sigma, partition) else call srb_energy(mol, lattr, cutoff%srb, mol%num, param%rvdw_srb, param%rscal, param%qscal, & & rexp_srb, zexp_srb, energies, partition) end if end if if (param%base) then if (grad) then call srb_deriv(mol, lattr, cutoff%srb, param%zeff, param%rvdw, param%rscal, param%qscal, & & rexp_base, zexp_base, energies, gradient, sigma, partition) else call srb_energy(mol, lattr, cutoff%srb, param%zeff, param%rvdw, param%rscal, param%qscal, & & rexp_base, zexp_base, energies, partition) end if end if end subroutine get_geometric_counterpoise_atomic !> Analytical second derivatives of the geometric counterpoise correction !> with respect to the Cartesian coordinates subroutine get_geometric_counterpoise_hessian(mol, param, cutoff, hessian, partition) !> Molecular structure data class(structure_type), intent(in) :: mol !> Geometric counterpoise parameters type(gcp_param), intent(in) :: param !> Realspace cutoffs type(realspace_cutoff), intent(in) :: cutoff !> Counter-poise hessian real(wp), intent(out), contiguous :: hessian(:, :) !> Work partition of the atom pairs, absent selects the complete work type(work_partition), intent(in), optional :: partition real(wp), allocatable :: lattr(:, :) hessian(:, :) = 0.0_wp if (allocated(param%emiss) .and. allocated(param%slater) .and. allocated(param%xv)) then call get_lattice_points(mol%periodic, mol%lattice, cutoff%gcp, lattr) call gcp_hessian(mol, lattr, cutoff%gcp, param%zeff, param%emiss, param%slater, & & param%xv, param%rvdw, param%sigma, param%alpha, param%beta, param%damp, & & param%dmp_scal, param%dmp_exp, hessian, partition) end if if (param%srb .or. param%base) then call get_lattice_points(mol%periodic, mol%lattice, cutoff%srb, lattr) end if if (param%srb) then call srb_hessian(mol, lattr, cutoff%srb, mol%num, param%rvdw_srb, param%rscal, & & param%qscal, rexp_srb, zexp_srb, hessian, partition) end if if (param%base) then call srb_hessian(mol, lattr, cutoff%srb, param%zeff, param%rvdw, param%rscal, & & param%qscal, rexp_base, zexp_base, hessian, partition) end if end subroutine get_geometric_counterpoise_hessian !> Geometric counterpoise correction subroutine gcp_energy(mol, trans, cutoff, iz, emiss, slater, xv, rvdw, escal, alpha, beta, & & damp, dmp_scal, dmp_exp, energies, partition) !> Molecular structure data type(structure_type), intent(in) :: mol !> Translation vectors real(wp), intent(in) :: trans(:, :) !> Distance cutoff real(wp), intent(in) :: cutoff !> Effective nuclear charges integer, intent(in) :: iz(:) !> Basis set superposition error per atom real(wp), intent(in) :: emiss(:) !> Slater exponents real(wp), intent(in) :: slater(:) !> Number of virtual orbitals real(wp), intent(in) :: xv(:) !> Van der Waals radii real(wp), intent(in) :: rvdw(:, :) !> Scaling factor real(wp), intent(in) :: escal !> Exponential factor real(wp), intent(in) :: alpha !> Power factor real(wp), intent(in) :: beta !> Damping flag logical, intent(in) :: damp !> Damping scaling factor real(wp), intent(in) :: dmp_scal !> Damping exponent real(wp), intent(in) :: dmp_exp !> Atom-resolved energy real(wp), intent(inout) :: energies(:) !> Work partition of the atom pairs, absent selects the complete work type(work_partition), intent(in), optional :: partition integer :: iat, jat, jtr, izp, jzp real(wp) :: argv real(wp) :: xvi, xvj real(wp) :: r1, rscal, rscalexp, r0 real(wp) :: sij, expv, bsse, emi, emj real(wp) :: vec(3) real(wp) :: dampval, grd_dmp real(wp) :: dE !$omp parallel do default(none) & !$omp reduction(+:energies) & !$omp shared(mol, iz, xv, emiss, rvdw, trans, cutoff, alpha, beta, & !$omp& dmp_scal, dmp_exp, escal, slater, damp, partition) & !$omp private(izp, jat, jzp, xvi, xvj, emi, emj, r0, jtr, vec, r1, sij, & !$omp& expv, bsse, argv, dE, dampval, grd_dmp, rscal, rscalexp) do iat = 1, mol%nat izp = mol%id(iat) xvi = merge(1.0_wp / sqrt(xv(izp)), 0.0_wp, xv(izp) >= 0.5_wp) ! the BSSE due to atom jat, Loop over all j atoms do jat = 1, iat if (.not.owns_pair(partition, iat, jat)) cycle jzp = mol%id(jat) xvj = merge(1.0_wp / sqrt(xv(jzp)), 0.0_wp, xv(jzp) >= 0.5_wp) emi = emiss(izp)*xvj*escal emj = emiss(jzp)*xvi*escal r0 = rvdw(izp, jzp) do jtr = 1, size(trans, 2) vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, jtr)) r1 = norm2(vec) if(r1 > cutoff .or. r1 < epsilon(1.0_wp)) cycle ! calulate slater overlap sij call ssovl(r1, izp, jzp, iz, slater(izp), slater(jzp), sij) ! evaluate gcp central expression argv = -alpha*r1**beta expv = exp(argv) bsse = expv/sqrt(sij) if (damp) then rscal = r1/r0 rscalexp = dmp_scal*rscal**dmp_exp dampval = (1.0_wp-1.0_wp/(1.0_wp+rscalexp)) else dampval = 1.0_wp end if dE = bsse*dampval energies(iat) = energies(iat) + emi*dE if (iat /= jat) then energies(jat) = energies(jat) + emj*dE end if end do end do end do end subroutine gcp_energy !> Geometric counterpoise correction subroutine gcp_deriv(mol, trans, cutoff, iz, emiss, slater, xv, rvdw, escal, alpha, beta, & & damp, dmp_scal, dmp_exp, energies, gradient, sigma, partition) !> Molecular structure data type(structure_type), intent(in) :: mol !> Translation vectors real(wp), intent(in) :: trans(:, :) !> Distance cutoff real(wp), intent(in) :: cutoff !> Effective nuclear charges integer, intent(in) :: iz(:) !> Basis set superposition error per atom real(wp), intent(in) :: emiss(:) !> Slater exponents real(wp), intent(in) :: slater(:) !> Number of virtual orbitals real(wp), intent(in) :: xv(:) !> Van der Waals radii real(wp), intent(in) :: rvdw(:, :) !> Scaling factor real(wp), intent(in) :: escal !> Exponential factor real(wp), intent(in) :: alpha !> Power factor real(wp), intent(in) :: beta !> Damping flag logical, intent(in) :: damp !> Damping scaling factor real(wp), intent(in) :: dmp_scal !> Damping exponent real(wp), intent(in) :: dmp_exp !> Atom-resolved energy real(wp), intent(inout) :: energies(:) !> Molecular gradient real(wp), intent(inout) :: gradient(:, :) !> Virial real(wp), intent(inout) :: sigma(:, :) !> Work partition of the atom pairs, absent selects the complete work type(work_partition), intent(in), optional :: partition integer :: iat, jat, jtr, izp, jzp real(wp) :: argv, argd, ovlpd real(wp) :: xvi, xvj real(wp) :: r1, rscal, rscalexp, rscalexpm1, r0 real(wp) :: sij, expv, expd, bsse, emi, emj, emij real(wp) :: vec(3), gs(3), gij real(wp) :: dampval, grd_dmp real(wp) :: dE, dG(3), dS(3, 3) !$omp parallel do default(none) & !$omp reduction(+:energies, gradient, sigma) & !$omp shared(mol, iz, xv, emiss, rvdw, trans, cutoff, alpha, beta, & !$omp& dmp_scal, dmp_exp, escal, slater, damp, partition) & !$omp private(izp, jat, jzp, xvi, xvj, emi, emj, r0, jtr, vec, r1, sij, gij, emij, & !$omp& expv, bsse, argv, dE, dampval, grd_dmp, dG, dS, expd, argd, ovlpd, gs, & !$omp& rscal, rscalexp, rscalexpm1) do iat = 1, mol%nat izp = mol%id(iat) xvi = merge(1.0_wp / sqrt(xv(izp)), 0.0_wp, xv(izp) >= 0.5_wp) ! the BSSE due to atom jat, Loop over all j atoms do jat = 1, iat if (.not.owns_pair(partition, iat, jat)) cycle jzp = mol%id(jat) xvj = merge(1.0_wp / sqrt(xv(jzp)), 0.0_wp, xv(jzp) >= 0.5_wp) emi = emiss(izp)*xvj*escal emj = emiss(jzp)*xvi*escal emij = emi + emj r0 = rvdw(izp, jzp) do jtr = 1, size(trans, 2) vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, jtr)) r1 = norm2(vec) if(r1 > cutoff .or. r1 < epsilon(1.0_wp)) cycle ! calulate slater overlap sij call ssovl(r1, izp, jzp, iz, slater(izp), slater(jzp), sij) ! evaluate gcp central expression argv = -alpha*r1**beta expv = exp(argv) bsse = expv/sqrt(sij) if (damp) then rscal = r1/r0 rscalexp = dmp_scal*rscal**dmp_exp dampval = (1.0_wp-1.0_wp/(1.0_wp+rscalexp)) else dampval = 1.0_wp end if call gsovl(r1, izp, jzp, iz, slater(izp), slater(jzp), gij) gs(:) = gij*vec expd = exp(-alpha*r1**beta)*(-0.5_wp) argd = 2d0*alpha*beta*r1**beta*sij/r1 ovlpd = r1*sij**1.5_wp if(damp) then rscalexpm1 = rscal**(dmp_exp-1) grd_dmp = dmp_scal*dmp_exp*rscalexpm1/r0 grd_dmp = grd_dmp/(rscalexp+1.0_wp)**2 end if dE = bsse*dampval dG = expd*(argd*vec + gs)/ovlpd*emij if(damp) then dG = dG*dampval+bsse*grd_dmp*(vec/r1)*emij end if dS = spread(dG, 1, 3) * spread(vec, 2, 3) * 0.5_wp energies(iat) = energies(iat) + emi*dE if (iat /= jat) then energies(jat) = energies(jat) + emj*dE end if sigma(:, :) = sigma(:, :) + dS if (iat /= jat) then gradient(:, iat) = gradient(:, iat) + dG gradient(:, jat) = gradient(:, jat) - dG sigma(:, :) = sigma(:, :) + dS end if end do end do end do end subroutine gcp_deriv !> Short-range bond length correction for HF-3c subroutine srb_energy(mol, trans, cutoff, iz, r0ab, rscal, qscal, rexp, zexp, energies, & & partition) !> Molecular structure data type(structure_type), intent(in) :: mol !> Translation vectors real(wp), intent(in) :: trans(:, :) !> Distance cutoff real(wp), intent(in) :: cutoff !> Effective nuclear charges integer, intent(in) :: iz(:) !> Van der Waals radii real(wp), intent(in) :: r0ab(:, :) !> Radii scaling factor real(wp), intent(in) :: rscal !> Prefactor for the SRB potential real(wp), intent(in) :: qscal !> Exponent for radii real(wp), intent(in) :: rexp !> Exponent for charges real(wp), intent(in) :: zexp !> Atom-resolved energy real(wp), intent(inout) :: energies(:) !> Work partition of the atom pairs, absent selects the complete work type(work_partition), intent(in), optional :: partition real(wp) :: fi, fj, ff, r1, expt real(wp) :: r0, vec(3), dE integer :: iat, jat, jtr, izp, jzp !$omp parallel do default(none) & !$omp reduction(+:energies) & !$omp shared(mol, trans, cutoff, iz, r0ab, rexp, zexp, rscal, qscal, partition) & !$omp private(izp, jat, jzp, r0, vec, r1, fi, fj, ff, expt, dE) do iat = 1, mol%nat izp = mol%id(iat) fi = real(iz(izp), wp) do jat = 1, iat if (.not.owns_pair(partition, iat, jat)) cycle jzp = mol%id(jat) r0 = rscal*r0ab(izp, jzp)**rexp fj = real(iz(jzp), wp) ff = -(fi*fj)**zexp do jtr = 1, size(trans, 2) vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, jtr)) r1 = norm2(vec) if(r1 > cutoff .or. r1 < epsilon(1.0_wp)) cycle expt = exp(-r0*r1) dE = qscal*ff*expt*0.5_wp energies(iat) = energies(iat) + dE if (iat /= jat) then energies(jat) = energies(jat) + dE end if end do end do end do end subroutine srb_energy !> Short-range bond length correction for HF-3c subroutine srb_deriv(mol, trans, cutoff, iz, r0ab, rscal, qscal, rexp, zexp, energies, & & gradient, sigma, partition) !> Molecular structure data type(structure_type), intent(in) :: mol !> Translation vectors real(wp), intent(in) :: trans(:, :) !> Distance cutoff real(wp), intent(in) :: cutoff !> Effective nuclear charges integer, intent(in) :: iz(:) !> Van der Waals radii real(wp), intent(in) :: r0ab(:, :) !> Radii scaling factor real(wp), intent(in) :: rscal !> Prefactor for the SRB potential real(wp), intent(in) :: qscal !> Exponent for radii real(wp), intent(in) :: rexp !> Exponent for charges real(wp), intent(in) :: zexp !> Atom-resolved energy real(wp), intent(inout) :: energies(:) !> Molecular gradient real(wp), intent(inout) :: gradient(:, :) !> Molecular virial real(wp), intent(inout) :: sigma(:, :) !> Work partition of the atom pairs, absent selects the complete work type(work_partition), intent(in), optional :: partition real(wp) :: fi, fj, ff, rf, r1, expt real(wp) :: r0, vec(3), dE, dG(3), dS(3, 3) integer :: iat, jat, jtr, izp, jzp !$omp parallel do default(none) & !$omp reduction(+:energies, gradient, sigma) & !$omp shared(mol, trans, cutoff, iz, r0ab, rexp, zexp, rscal, qscal, partition) & !$omp private(izp, jat, jzp, r0, vec, r1, fi, fj, ff, expt, rf, dE, dG, dS) do iat = 1, mol%nat izp = mol%id(iat) fi = real(iz(izp), wp) do jat = 1, iat if (.not.owns_pair(partition, iat, jat)) cycle jzp = mol%id(jat) r0 = rscal*r0ab(izp, jzp)**rexp fj = real(iz(jzp), wp) ff = -(fi*fj)**zexp do jtr = 1, size(trans, 2) vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, jtr)) r1 = norm2(vec) if(r1 > cutoff .or. r1 < epsilon(1.0_wp)) cycle expt = exp(-r0*r1) rf = qscal/r1 dE = qscal*ff*expt*0.5_wp dG(:) = -ff*r0*vec*expt*rf dS(:, :) = spread(dG, 1, 3) * spread(vec, 2, 3) * 0.5_wp energies(iat) = energies(iat) + dE if (iat /= jat) then energies(jat) = energies(jat) + dE end if sigma(:, :) = sigma + dS if (iat /= jat) then gradient(:, iat) = gradient(:, iat) + dG gradient(:, jat) = gradient(:, jat) - dG sigma(:, :) = sigma + dS end if end do end do end do end subroutine srb_deriv !> Second derivatives of the geometric counterpoise correction subroutine gcp_hessian(mol, trans, cutoff, iz, emiss, slater, xv, rvdw, escal, alpha, beta, & & damp, dmp_scal, dmp_exp, hessian, partition) !> Molecular structure data type(structure_type), intent(in) :: mol !> Translation vectors real(wp), intent(in) :: trans(:, :) !> Distance cutoff real(wp), intent(in) :: cutoff !> Effective nuclear charges integer, intent(in) :: iz(:) !> Basis set superposition error per atom real(wp), intent(in) :: emiss(:) !> Slater exponents real(wp), intent(in) :: slater(:) !> Number of virtual orbitals real(wp), intent(in) :: xv(:) !> Van der Waals radii real(wp), intent(in) :: rvdw(:, :) !> Scaling factor real(wp), intent(in) :: escal !> Exponential factor real(wp), intent(in) :: alpha !> Power factor real(wp), intent(in) :: beta !> Damping flag logical, intent(in) :: damp !> Damping scaling factor real(wp), intent(in) :: dmp_scal !> Damping exponent real(wp), intent(in) :: dmp_exp !> Second derivative of the energy w.r.t. the Cartesian coordinates real(wp), intent(inout), contiguous :: hessian(:, :) !> Work partition of the atom pairs, absent selects the complete work type(work_partition), intent(in), optional :: partition integer :: iat, jat, jtr, izp, jzp real(wp) :: xvi, xvj, emi, emj, emij real(wp) :: r0, r1, vec(3), rscal, rscalexp, dmpq1, dmpq2 real(wp) :: sij, gij, hij, srt, ovl0, ovl1, ovl2 real(wp) :: expv, expv1, expv2, arg1, arg2, rbeta real(wp) :: bsse, bsse1, bsse2 real(wp) :: dampval, damp1, damp2 real(wp) :: dE1, dE2 do iat = 1, mol%nat izp = mol%id(iat) xvi = merge(1.0_wp / sqrt(xv(izp)), 0.0_wp, xv(izp) >= 0.5_wp) do jat = 1, iat - 1 if (.not.owns_pair(partition, iat, jat)) cycle jzp = mol%id(jat) xvj = merge(1.0_wp / sqrt(xv(jzp)), 0.0_wp, xv(jzp) >= 0.5_wp) emi = emiss(izp)*xvj*escal emj = emiss(jzp)*xvi*escal emij = emi + emj r0 = rvdw(izp, jzp) do jtr = 1, size(trans, 2) vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, jtr)) r1 = norm2(vec) if(r1 > cutoff .or. r1 < epsilon(1.0_wp)) cycle call dsovl(r1, izp, jzp, iz, slater(izp), slater(jzp), sij, gij, hij) ! derivatives of the inverse square root of the overlap srt = sqrt(sij) ovl0 = 1.0_wp/srt ovl1 = -0.5_wp*gij/(sij*srt) ovl2 = (0.75_wp*gij*gij/sij - 0.5_wp*hij)/(sij*srt) ! derivatives of the exponential prefactor rbeta = r1**beta arg1 = alpha*beta*rbeta/r1 arg2 = arg1*(beta-1.0_wp)/r1 expv = exp(-alpha*rbeta) expv1 = -arg1*expv expv2 = (arg1*arg1 - arg2)*expv bsse = expv*ovl0 bsse1 = expv1*ovl0 + expv*ovl1 bsse2 = expv2*ovl0 + 2.0_wp*expv1*ovl1 + expv*ovl2 if (damp) then rscal = r1/r0 rscalexp = dmp_scal*rscal**dmp_exp dmpq1 = dmp_exp*rscalexp/(rscal*r0) dmpq2 = dmpq1*(dmp_exp-1.0_wp)/(rscal*r0) dampval = rscalexp/(1.0_wp+rscalexp) damp1 = dmpq1/(1.0_wp+rscalexp)**2 damp2 = dmpq2/(1.0_wp+rscalexp)**2 & & - 2.0_wp*dmpq1*dmpq1/(1.0_wp+rscalexp)**3 else dampval = 1.0_wp damp1 = 0.0_wp damp2 = 0.0_wp end if dE1 = emij*(bsse1*dampval + bsse*damp1) dE2 = emij*(bsse2*dampval + 2.0_wp*bsse1*damp1 + bsse*damp2) call add_pair_hessian(hessian, iat, jat, vec, r1, dE1, dE2) end do end do end do end subroutine gcp_hessian !> Second derivatives of the short-range bond length correction for HF-3c subroutine srb_hessian(mol, trans, cutoff, iz, r0ab, rscal, qscal, rexp, zexp, hessian, & & partition) !> Molecular structure data type(structure_type), intent(in) :: mol !> Translation vectors real(wp), intent(in) :: trans(:, :) !> Distance cutoff real(wp), intent(in) :: cutoff !> Effective nuclear charges integer, intent(in) :: iz(:) !> Van der Waals radii real(wp), intent(in) :: r0ab(:, :) !> Radii scaling factor real(wp), intent(in) :: rscal !> Prefactor for the SRB potential real(wp), intent(in) :: qscal !> Exponent for radii real(wp), intent(in) :: rexp !> Exponent for charges real(wp), intent(in) :: zexp !> Second derivative of the energy w.r.t. the Cartesian coordinates real(wp), intent(inout), contiguous :: hessian(:, :) !> Work partition of the atom pairs, absent selects the complete work type(work_partition), intent(in), optional :: partition real(wp) :: fi, fj, ff, r1, expt real(wp) :: r0, vec(3), dE1, dE2 integer :: iat, jat, jtr, izp, jzp do iat = 1, mol%nat izp = mol%id(iat) fi = real(iz(izp), wp) do jat = 1, iat - 1 if (.not.owns_pair(partition, iat, jat)) cycle jzp = mol%id(jat) r0 = rscal*r0ab(izp, jzp)**rexp fj = real(iz(jzp), wp) ff = -(fi*fj)**zexp do jtr = 1, size(trans, 2) vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, jtr)) r1 = norm2(vec) if(r1 > cutoff .or. r1 < epsilon(1.0_wp)) cycle expt = exp(-r0*r1) dE1 = -qscal*ff*r0*expt dE2 = qscal*ff*r0*r0*expt call add_pair_hessian(hessian, iat, jat, vec, r1, dE1, dE2) end do end do end do end subroutine srb_hessian !> Distribute the second derivatives of a radial pair potential to the hessian pure subroutine add_pair_hessian(hessian, iat, jat, vec, r1, dE1, dE2) !> Second derivative of the energy w.r.t. the Cartesian coordinates real(wp), intent(inout), contiguous :: hessian(:, :) !> Atom indices of the interacting pair integer, intent(in) :: iat, jat !> Distance vector and its norm real(wp), intent(in) :: vec(3), r1 !> First and second derivative of the pair potential w.r.t. the distance real(wp), intent(in) :: dE1, dE2 integer :: ic, jc, ii, jj real(wp) :: hblk(3, 3), fr, fc fr = dE1/r1 fc = (dE2 - fr)/(r1*r1) do ic = 1, 3 do jc = 1, 3 hblk(jc, ic) = fc*vec(ic)*vec(jc) end do hblk(ic, ic) = hblk(ic, ic) + fr end do ii = 3*(iat - 1) jj = 3*(jat - 1) do ic = 1, 3 do jc = 1, 3 hessian(jj+jc, jj+ic) = hessian(jj+jc, jj+ic) + hblk(jc, ic) hessian(ii+jc, ii+ic) = hessian(ii+jc, ii+ic) + hblk(jc, ic) hessian(jj+jc, ii+ic) = hessian(jj+jc, ii+ic) - hblk(jc, ic) hessian(ii+jc, jj+ic) = hessian(ii+jc, jj+ic) - hblk(jc, ic) end do end do end subroutine add_pair_hessian !****************************************************************************** !* calculates the s-type overlap integral over 1s, 2s and 3s slater functions !* added support for 3s functions !* ovl = overlap integral !* za = slater exponent atom A !* zb = slater exponent atom B !* R = distance between atom A and B !* Inspired by mopac7.0 !****************************************************************************** subroutine ssovl(r, iat, jat, iz, xza, xzb, ovl) real(wp), intent(in) :: r integer, intent(in) :: iat, jat integer, intent(in) :: iz(:) real(wp), intent(in) :: xza, xzb real(wp), intent(out) :: ovl integer :: ii, shell(72) logical :: debug real(wp) :: za, zb, ax, bx, norm, R05 integer :: na, nb real(wp) :: xx data shell/ & ! h, he 1, 1 & ! li-ne , 2, 2, 2, 2, 2, 2, 2, 2, & ! na-ar 3, 3, 3, 3, 3, 3, 3, 3, & ! 4s, 5s will be treated as 3s ! k-rn , no f-elements 54*3/ ! ... za = xza zb = xzb na = iz(iat) nb = iz(jat) debug = .false. !debug = .true. ! ii selects kind of ovl by multiplying the shell ! kind <1s|1s> <2s|1s> <2s|2s> <1s|3s> <2s|3s> <3s|3s> ! case: 1 2 4 3 6 9 ! ii = shell(na)*shell(nb) if(debug) write(*, *) "shell", ii R05 = R*0.5 ax = (za+zb)*R05 bx = (zb-za)*R05 ! same elements if(abs(za-zb) < 0.1) then select case (ii) case (1) ovl = 0.25d0*sqrt((za*zb*R*R)**3)*(A2(ax)*Bint(bx, 0)-Bint(bx, 2)*A0(ax)) case (2) ovl = SQRT(1.0_wp/3.0_wp) if(shell(na) < shell(nb)) then ! <1s|2s> norm = SQRT((ZA**3)*(ZB**5))*(R**4)*0.125_wp ovl = ovl*norm*(A3(ax)*Bint(bx, 0)-Bint(bx, 3)*A0(ax)+A2(ax)*Bint(bx, 1)-Bint(bx, 2)*A1(ax)) else ! switch za/zb to get <2s|1s> xx = za za = zb zb = xx ax = (za+zb)*R05 bx = (zb-za)*R05 norm = SQRT((ZA**3)*(ZB**5))*(R**4)*0.125_wp ovl = ovl*norm*(A3(ax)*Bint(bx, 0)-Bint(bx, 3)*A0(ax)+A2(ax)*Bint(bx, 1)-Bint(bx, 2)*A1(ax)) end if case (4) norm = SQRT((ZA*ZB)**5)*(R**5)*0.0625d0 ovl = norm* (A4(ax)*Bint(bx, 0)+Bint(bx, 4)*A0(ax)-2.0d0*A2(ax)*Bint(bx, 2))*(1d0/3d0) case(3) if(shell(na) < shell(nb)) then norm = SQRT((ZA**3)*(ZB**7)/7.5_wp)*(R**5)*0.0625_wp ovl = norm*(A4(ax)*Bint(bx, 0)-Bint(bx, 4)*A0(ax)+2.0d0*(A3(ax)*Bint(bx, 1)-Bint(bx, 3)*A1(ax)))/sqrt(3.0d0) else xx = za za = zb zb = xx ax = (za+zb)*R05 bx = (zb-za)*R05 norm = SQRT((ZA**3)*(ZB**7)/7.5_wp)*(R**5)*0.0625_wp ovl = norm*(A4(ax)*Bint(bx, 0)-Bint(bx, 4)*A0(ax)+2.0d0*(A3(ax)*Bint(bx, 1)-Bint(bx, 3)*A1(ax)))/sqrt(3.0d0) end if case(6) if(shell(na) < shell(nb)) then norm = SQRT((za**5)*(zb**7)/7.5_wp)*(R**6)*0.03125_wp ovl = norm*(A5(ax)*Bint(bx, 0)+A4(ax)*Bint(bx, 1) & & -2d0*(A3(ax)*Bint(bx, 2)+A2(ax)*Bint(bx, 3)) & & +A1(ax)*Bint(bx, 4)+A0(ax)*Bint(bx, 5))/3.0d0 else xx = za za = zb zb = xx ax = (za+zb)*R05 bx = (zb-za)*R05 norm = SQRT((za**5)*(zb**7)/7.5_wp)*(R**6)*0.03125_wp ovl = norm*(A5(ax)*Bint(bx, 0)+A4(ax)*Bint(bx, 1) & & -2d0*(A3(ax)*Bint(bx, 2)+A2(ax)*Bint(bx, 3)) & & +A1(ax)*Bint(bx, 4)+A0(ax)*Bint(bx, 5))/3.0d0 end if case(9) norm = sqrt((ZA*ZB*R*R)**7)/480.0d0 ovl = norm*(A6(ax)*Bint(bx, 0)-3.0d0*(A4(ax)*Bint(bx, 2) & & -A2(ax)*Bint(bx, 4))-A0(ax)*Bint(bx, 6))/3.0_wp case default error stop "invalid shell combination in ssovl" end select else ! different elements select case (ii) case (1) norm = 0.25d0*sqrt((za*zb*R*R)**3) ovl = (A2(ax)*B0(bx)-B2(bx)*A0(ax))*norm case (2) ovl = SQRT(1.0_wp/3.0_wp) if(shell(na) < shell(nb)) then ! <1s|2s> norm = SQRT((ZA**3)*(ZB**5))*(R**4)*0.125_wp ovl = ovl*norm*(A3(ax)*B0(bx)-B3(bx)*A0(ax)+A2(ax)*B1(bx)-B2(bx)*A1(ax)) else ! switch za/zb to get <2s|1s> xx = za za = zb zb = xx ax = (za+zb)*R05 bx = (zb-za)*R05 norm = SQRT((ZA**3)*(ZB**5))*(R**4)*0.125_wp ovl = ovl*norm*(A3(ax)*B0(bx)-B3(bx)*A0(ax)+A2(ax)*B1(bx)-B2(bx)*A1(ax)) end if case (4) ! <2s|2s> norm = SQRT((ZA*ZB)**5)*(R**5)*0.0625_wp ovl = norm* (A4(ax)*B0(bx)+B4(bx)*A0(ax)-2.0_wp*A2(ax)*B2(bx))*(1d0/3d0) case(3) ! <1s|3s> + <3s|1s> if(shell(na) < shell(nb)) then norm = SQRT((ZA**3)*(ZB**7)/7.5_wp)*(R**5)*0.0625_wp ovl = norm*(A4(ax)*B0(bx)-B4(bx)*A0(ax)+2.0d0*(A3(ax)*B1(bx)-B3(bx)*A1(ax)))/sqrt(3.0d0) else xx = za za = zb zb = xx ax = (za+zb)*R05 bx = (zb-za)*R05 norm = SQRT((ZA**3)*(ZB**7)/7.5_wp)*(R**5)*0.0625_wp ovl = norm*(A4(ax)*B0(bx)-B4(bx)*A0(ax)+2.0d0*(A3(ax)*B1(bx)-B3(bx)*A1(ax)))/sqrt(3.0d0) end if case(6) ! <2s|3s> + <3s|2s> if(shell(na) < shell(nb)) then norm = SQRT((za**5)*(zb**7)/7.5_wp)*(R**6)*0.03125_wp ovl = norm*(A5(ax)*B0(bx)+A4(ax)*B1(bx)-2d0*(A3(ax)*B2(bx)+A2(ax)*B3(bx))+A1(ax)*B4(bx)+A0(ax)*B5(bx))/3.0d0 else xx = za za = zb zb = xx ax = (za+zb)*R05 bx = (zb-za)*R05 norm = SQRT((za**5)*(zb**7)/7.5_wp)*(R**6)*0.03125_wp ovl = norm*(A5(ax)*B0(bx)+A4(ax)*B1(bx)-2.0_wp*(A3(ax)*B2(bx)+A2(ax)*B3(bx))+A1(ax)*B4(bx)+A0(ax)*B5(bx))/3.0d0 end if case(9) ! <3s|3> norm = sqrt((ZA*ZB*R*R)**7)/1440.0d0 ovl = norm*(A6(ax)*B0(bx)-3.0d0*(A4(ax)*B2(bx)-A2(ax)*B4(bx))-A0(ax)*B6(bx)) case default error stop "invalid shell combination in ssovl" end select end if end subroutine ssovl !**************************************** !* A(x) auxiliary integrals * !* Quantenchemie - Ein Lehrgang Vol 5 * !* p. 570 eq. 11.4.14 * !**************************************** real(wp) pure function A0(x) ! Hilfsintegral A_0 implicit none real(wp), intent(in) :: x A0 = exp(-x)/x return end function A0 real(wp) pure function A1(x) ! Hilfsintegral A_1 implicit none real(wp), intent(in) :: x A1 = ((1+x)*exp(-x))/(x**2) return end function A1 real(wp) pure function A2(x) ! Hilfsintegral A_2 implicit none real(wp), intent(in) :: x A2 = ((2d0+2d0*x+x**2)*exp(-x))/x**3 return end function A2 real(wp) pure function A3(x) ! Hilfsintegral A_3 implicit none real(wp), intent(in) :: x real(wp) :: xx real(wp) :: x2, x3, x4 x2 = x*x x3 = x2*x x4 = x3*x xx = (6d0+6d0*x+3d0*x2+x3) A3 = (xx*exp(-x))/x4 return end function A3 real(wp) pure function A4(x) ! Hilfsintegral A_4 implicit none real(wp), intent(in) :: x real(wp) :: xx real(wp) :: x2, x3, x4, x5 x2 = x*x x3 = x2*x x4 = x3*x x5 = x4*x xx = (24d0+24d0*x+12d0*x2+4d0*x3+x4) A4 = (xx*exp(-x))/x5 return end function A4 real(wp) pure function A5(x) ! Hilfsintegral A_5 implicit none real(wp), intent(in) :: x real(wp) :: xx real(wp) :: x2, x3, x4, x5, x6 x2 = x*x x3 = x2*x x4 = x3*x x5 = x4*x x6 = x5*x xx = (120d0+120d0*x+60d0*x2+20d0*x3+5d0*x4+x5) A5 = (xx*exp(-x))/x6 return end function A5 real(wp) pure function A6(x) ! Hilfsintegral A_6 implicit none real(wp), intent(in) :: x real(wp) :: xx real(wp) :: x2, x3, x4, x5, x6, x7 x2 = x*x x3 = x2*x x4 = x3*x x5 = x4*x x6 = x5*x x7 = x6*x xx = (720d0+720d0*x+360d0*x2+120d0*x3+30d0*x4+6d0*x5+x6) A6 = (xx*exp(-x))/x7 return end function A6 !************************************** !* B(x) auxiliary integrals * !* Quantenchemie - Ein Lehrgang Vol 5 * !* p. 570 eq. 11.4.14b * !************************************** real(wp) pure function B0(x) real(wp), intent(in) :: x B0 = (exp(x)-exp(-x))/x end function B0 real(wp) pure function B1(x) real(wp), intent(in) :: x real(wp) :: x2, x3 x2 = x*x x3 = x2*x B1 = ((1.0_wp-x)*exp(x)-(1.0_wp+x)*exp(-x))/x2 end function B1 real(wp) pure function B2(x) real(wp), intent(in) :: x real(wp) :: x2, x3 x2 = x*x x3 = x2*x B2 = (((2.0_wp-2*x+x2)*exp(x)) - ((2.0_wp+2.0_wp*x+x2)*exp(-x)))/x3 end function B2 real(wp) pure function B3(x) real(wp), intent(in) :: x real(wp) :: xx, yy real(wp) :: x2, x3, x4 x2 = x*x x3 = x2*x x4 = x3*x xx = (6.0_wp-6.0_wp*x+3.0_wp*x2-x3)*exp(x)/x4 yy = (6.0_wp+6.0_wp*x+3.0_wp*x2+x3)*exp(-x)/x4 B3 = xx-yy end function B3 real(wp) pure function B4(x) real(wp), intent(in) :: x real(wp) :: xx, yy real(wp) :: x2, x3, x4, x5 x2 = x*x x3 = x2*x x4 = x3*x x5 = x4*x xx = (24.0_wp-24.0_wp*x+12.0_wp*x2-4.0_wp*x3+x4)*exp(x)/x5 yy = (24.0_wp+24.0_wp*x+12.0_wp*x2+4.0_wp*x3+x4)*exp(-x)/x5 B4 = xx-yy end function B4 real(wp) pure function B5(x) real(wp), intent(in) :: x real(wp) :: xx, yy real(wp) :: x2, x3, x4, x5, x6 x2 = x*x x3 = x2*x x4 = x3*x x5 = x4*x x6 = x5*x xx = (120.0_wp-120*x+60*x2-20*x3+5*x4-x5)*exp(x)/x6 yy = (120.0_wp+120*x+60*x2+20*x3+5*x4+x5)*exp(-x)/x6 B5 = xx-yy end function B5 real(wp) function B6(x) real(wp), intent(in) :: x real(wp) :: x2, x3, x4, x5, x6, x7, yy, xx x2 = x*x x3 = x2*x x4 = x3*x x5 = x4*x x6 = x5*x x7 = x6*x xx = (720.0_wp - 720.0_wp*x+ 360.0_wp*x2 - 120.0_wp*x3 + 30.0_wp*x4 - 6.0_wp*x5 + x6)*exp(x)/x7 yy = (720.0_wp + 720.0_wp*x + 360.0_wp*x2 + 120.0_wp*x3 + 30.0_wp*x4 + 6.0_wp*x5 + x6)*exp(-x)/x7 B6 = xx-yy end function B6 real(wp) function bint(x, k) ! calculates B_k(x) ! general summation formula ! 'infinite' sum is numerically unstable. 12 terms seem ! accurate enough implicit none real(wp), intent(in) :: x real(wp) :: xx, yy integer, intent(in) :: k integer :: i bint = 0 if(abs(x)<1e-6) then do i = 0, k bint = (1.0d0+(-1d0)**i)/(dble(i)+1.0d0) end do return end if do i = 0, 12 xx = 1d0-((-1d0)**(k+i+1)) yy = dble(fact(i))*dble((k+i+1)) bint = bint+xx/yy*(-x)**i end do end function bint ! faculty function integer(wp) function fact(N) implicit none integer, intent(in) :: n integer :: j fact = 1 do j = 2, n fact = fact*j end do return end function fact subroutine gsovl(r, iat, jat, iz, xza, xzb, g) ! GRADIENT ! calculates the s-type overlap integral over 1s, 2s, 3s slater functions ! ovl = overlap integral ! za = slater exponent atom A ! zb = slater exponent atom B ! R = distance between atom A and B real(wp), intent(in) :: r integer, intent(in) :: iat, jat integer, intent(in) :: iz(:) real(wp), intent(in) :: xza, xzb real(wp), intent(out) :: g integer :: ii, shell(72) logical :: debug real(wp) :: ax, bx, R05, za, zb integer :: na, nb data shell/ & ! h, he 1, 1 & ! li-ne , 2, 2, 2, 2, 2, 2, 2, 2, & ! na-ar 3, 3, 3, 3, 3, 3, 3, 3, & ! 4s, 5s will be treated as 3s ! k-rn , no f-elements 54*3/ ! ... real(wp) :: Fa, Fb !--------------------- set exponents --------------------------------------- real(wp) :: xx real(wp) :: s0, s2 logical :: lsame za = xza zb = xzb na = iz(iat) nb = iz(jat) !---------------------------------------------------------------------------- debug = .false. !debug = .true. ! ii selects kind of ovl by multiplying the shell ! kind <1s|1s> <2s|1s> <2s|2s> <1s|3s> <2s|3s> <3s|3s> ! case: 1 2 4 3 6 9 ! ii = shell(na)*shell(nb) if(debug) write(*, *) "gshell", ii R05 = R*0.5_wp ax = (za+zb)*R05 Fa = (za+zb) bx = (zb-za)*R05 Fb = (zb-za) lsame = .false. ! ! Nearly equal exponents: the closed-form expressions below are singular for ! zb -> za, evaluate the auxiliary integrals with the stable series instead if(abs(za-zb) < 0.1) then call dsovl(r, iat, jat, iz, xza, xzb, s0, g, s2) else ! different elements select case (ii) case (1) call g1s1s(za, zb, Fa, Fb, R, g, lsame) case (2) ! <1s|2s> if(shell(na) < shell(nb)) then call g2s1s(za, zb, Fa, Fb, R, g, .false., lsame) else xx = za za = zb zb = xx call g2s1s(za, zb, Fa, Fb, R, g, .true., lsame) end if case (4) ! <2s|2s> call g2s2s(za, zb, Fa, Fb, R, g, lsame) case(3) ! <1s|3s> + <3s|1s> if(shell(na) < shell(nb)) then call g1s3s(za, zb, Fa, Fb, R, g, .false., lsame) else xx = za za = zb zb = xx call g1s3s(za, zb, Fa, Fb, R, g, .true., lsame) end if case(6) ! <2s|3s> + <3s|2s> if(shell(na) < shell(nb)) then call g2s3s(za, zb, Fa, Fb, R, g, .false., lsame) else xx = za za = zb zb = xx call g2s3s(za, zb, Fa, Fb, R, g, .true., lsame) end if case(9) ! <3s|3> call g3s3s(za, zb, Fa, Fb, R, g, lsame) case default error stop "invalid shell combination in gsovl" end select end if end subroutine gsovl !> Second derivative of the s-type overlap integral with respect to the distance. subroutine hsovl(r, iat, jat, iz, xza, xzb, h) real(wp), intent(in) :: r integer, intent(in) :: iat, jat integer, intent(in) :: iz(:) real(wp), intent(in) :: xza, xzb real(wp), intent(out) :: h real(wp) :: s0, s1 call dsovl(r, iat, jat, iz, xza, xzb, s0, s1, h) end subroutine hsovl !> Value and first two derivatives of the s-type overlap integral. !> !> The overlap is written as norm * R**m * sum_t w_t * A_p(ax) * B_q(bx) with !> ax = (za+zb)*R/2 and bx = (zb-za)*R/2. Since dA_k/dx = -A_(k+1) and !> dB_k/dx = -B_(k+1) the derivatives follow from shifted auxiliary integrals. subroutine dsovl(r, iat, jat, iz, xza, xzb, s0, s1, s2) real(wp), intent(in) :: r integer, intent(in) :: iat, jat integer, intent(in) :: iz(:) real(wp), intent(in) :: xza, xzb real(wp), intent(out) :: s0, s1, s2 integer :: ii, shell(72) integer :: na, nb, nterm, m, it, ip, iq, k integer :: pa(6), qb(6) real(wp) :: wt(6) real(wp) :: za, zb, xx, ax, bx, ha, hb, cnorm real(wp) :: f0, f1, f2, av(0:8), bv(0:8) logical :: lsame data shell/ & ! h, he 1, 1 & ! li-ne , 2, 2, 2, 2, 2, 2, 2, 2, & ! na-ar 3, 3, 3, 3, 3, 3, 3, 3, & ! 4s, 5s will be treated as 3s ! k-rn , no f-elements 54*3/ za = xza zb = xzb na = iz(iat) nb = iz(jat) lsame = abs(za-zb) < 0.1_wp ! ii selects kind of ovl by multiplying the shell ! kind <1s|1s> <2s|1s> <2s|2s> <1s|3s> <2s|3s> <3s|3s> ! case: 1 2 4 3 6 9 ii = shell(na)*shell(nb) select case(ii) case(1) m = 3 nterm = 2 wt(:nterm) = [1.0_wp, -1.0_wp] pa(:nterm) = [2, 0] qb(:nterm) = [0, 2] cnorm = 0.25_wp*sqrt((za*zb)**3) case(2) if (shell(na) >= shell(nb)) then xx = za za = zb zb = xx end if m = 4 nterm = 4 wt(:nterm) = [1.0_wp, -1.0_wp, 1.0_wp, -1.0_wp] pa(:nterm) = [3, 0, 2, 1] qb(:nterm) = [0, 3, 1, 2] cnorm = sqrt(1.0_wp/3.0_wp)*sqrt((za**3)*(zb**5))*0.125_wp case(3) if (shell(na) >= shell(nb)) then xx = za za = zb zb = xx end if m = 5 nterm = 4 wt(:nterm) = [1.0_wp, -1.0_wp, 2.0_wp, -2.0_wp] pa(:nterm) = [4, 0, 3, 1] qb(:nterm) = [0, 4, 1, 3] cnorm = sqrt((za**3)*(zb**7)/7.5_wp)*0.0625_wp/sqrt(3.0_wp) case(4) m = 5 nterm = 3 wt(:nterm) = [1.0_wp, 1.0_wp, -2.0_wp] pa(:nterm) = [4, 0, 2] qb(:nterm) = [0, 4, 2] cnorm = sqrt((za*zb)**5)*0.0625_wp/3.0_wp case(6) if (shell(na) >= shell(nb)) then xx = za za = zb zb = xx end if m = 6 nterm = 6 wt(:nterm) = [1.0_wp, 1.0_wp, -2.0_wp, -2.0_wp, 1.0_wp, 1.0_wp] pa(:nterm) = [5, 4, 3, 2, 1, 0] qb(:nterm) = [0, 1, 2, 3, 4, 5] cnorm = sqrt((za**5)*(zb**7)/7.5_wp)*0.03125_wp/3.0_wp case(9) m = 7 nterm = 4 wt(:nterm) = [1.0_wp, -3.0_wp, 3.0_wp, -1.0_wp] pa(:nterm) = [6, 4, 2, 0] qb(:nterm) = [0, 2, 4, 6] cnorm = sqrt((za*zb)**7)/1440.0_wp case default error stop "invalid shell combination in dsovl" end select ha = 0.5_wp*(za + zb) hb = 0.5_wp*(zb - za) ax = ha*r bx = hb*r call aaux_all(ax, av) if (lsame) then call bint_all(bx, bv) else call baux_all(bx, bv) end if f0 = 0.0_wp f1 = 0.0_wp f2 = 0.0_wp do it = 1, nterm ip = pa(it) iq = qb(it) f0 = f0 + wt(it)*av(ip)*bv(iq) f1 = f1 - wt(it)*(ha*av(ip+1)*bv(iq) + hb*av(ip)*bv(iq+1)) f2 = f2 + wt(it)*(ha*ha*av(ip+2)*bv(iq) & & + 2.0_wp*ha*hb*av(ip+1)*bv(iq+1) + hb*hb*av(ip)*bv(iq+2)) end do s0 = cnorm*r**m*f0 s1 = cnorm*(real(m, wp)*r**(m-1)*f0 + r**m*f1) s2 = cnorm*(real(m*(m-1), wp)*r**(m-2)*f0 + 2.0_wp*real(m, wp)*r**(m-1)*f1 & & + r**m*f2) end subroutine dsovl !> All auxiliary integrals A_k(x) from the recursion x*A_k = k*A_(k-1) + exp(-x) pure subroutine aaux_all(x, a) real(wp), intent(in) :: x real(wp), intent(out) :: a(0:) real(wp) :: ex, rx integer :: k ex = exp(-x) rx = 1.0_wp/x a(0) = ex*rx do k = 1, ubound(a, 1) a(k) = (real(k, wp)*a(k-1) + ex)*rx end do end subroutine aaux_all !> All auxiliary integrals B_k(x), sharing the two exponentials between orders pure subroutine baux_all(x, b) real(wp), intent(in) :: x real(wp), intent(out) :: b(0:) real(wp) :: ep, em, rx, term, sp, sm, sgn integer :: j, k ep = exp(x) em = exp(-x) rx = 1.0_wp/x do k = 0, ubound(b, 1) term = rx sgn = merge(1.0_wp, -1.0_wp, mod(k, 2) == 0) sp = sgn*term sm = term do j = 1, k term = term*real(k-j+1, wp)*rx sgn = -sgn sp = sp + sgn*term sm = sm + term end do b(k) = ep*sp - em*sm end do end subroutine baux_all !> All auxiliary integrals B_k(x) from the series expansion used for small x, !> the terms (-x)**i/i! are shared between the orders pure subroutine bint_all(x, b) real(wp), intent(in) :: x real(wp), intent(out) :: b(0:) integer, parameter :: nterm = 12 real(wp) :: pw(0:nterm), acc integer :: i, k if (abs(x) < 1e-6_wp) then b(:) = 0.0_wp do k = 0, ubound(b, 1), 2 b(k) = 2.0_wp/(real(k, wp) + 1.0_wp) end do return end if pw(0) = 1.0_wp do i = 1, nterm pw(i) = pw(i-1)*(-x)/real(i, wp) end do ! only terms with k+i even contribute, they enter with a factor of two do k = 0, ubound(b, 1) acc = 0.0_wp do i = mod(k, 2), nterm, 2 acc = acc + pw(i)/real(k+i+1, wp) end do b(k) = 2.0_wp*acc end do end subroutine bint_all !------------------------------------------------------------- ! Maple was used to find the analy. derivatives of ! the slater integrals (expressions over A, B aux. integrals) ! Optimized fortran code by maple with some human-corrections !------------------------------------------------------------- subroutine g1s1s(za, zb, Fa, Fb, R, g, sameElement) ! slater overlap derv. ! derivative of explicit integral expression ! using maple implicit real(wp) (t) real(wp), intent(in) :: za, zb, Fa, Fb real(wp), intent(in) :: R real(wp), intent(out) :: g logical, intent(in) :: sameElement if(sameElement) then t1 = za ** 2 t3 = zb ** 2 t5 = t1 * za * t3 * zb t6 = R ** 2 t7 = t6 ** 2 t10 = Fa * R t14 = exp(-0.5_wp * t10) t17 = sqrt(t5 * t7 * t6) g = -(1.0_wp/3.0_wp) * t5 * t7 / Fa * (0.2D1 + t10) * t14 / t17 else t1 = za ** 2 t3 = zb ** 2 t5 = t1 * za * t3 * zb t6 = Fb ** 2 t7 = Fb * R t8 = 0.5_wp * t7 t9 = exp(t8) t12 = exp(-t8) t15 = t6 * Fa t22 = Fa ** 2 t23 = t22 * t9 t27 = t22 * t12 t31 = t6 * Fb t32 = R * t31 t37 = t22 * Fa t38 = R * t37 t43 = R ** 2 t44 = t43 * t31 t51 = t43 * t37 t56 = 0.4D1 * t6 * t9 - 0.4D1 * t6 * t12 + 0.2D1 * t15 * R * t9 - & 0.2D1 * t15 * R * t12 - 0.4D1 * t23 + 0.2D1 * t23 * t7 + 0.4D1 * t27 & + 0.2D1 * t27 * t7 - 0.2D1 * t32 * t9 - 0.2D1 * t32 * t12 - 0.2D1 * t38 * & t9 + 0.2D1 * t38 * t12 - 0.1D1 * t44 * Fa * t9 - 0.1D1 & * t44 * Fa * t12 + t51 * t9 * Fb + t51 * t12 * Fb t61 = exp(-0.5_wp * Fa * R) t62 = t43 ** 2 t65 = sqrt(t5 * t62 * t43) g = -0.2D1 * t5 * R * t56 * t61 / t65 / t31 / t37 end if end subroutine g1s1s subroutine g2s1s(za, zb, Fa, Fb, R, g, switch, lsame) ! slater overlap derv. ! derivative of explicit integral expression ! using maple implicit real(wp) (t) real(wp), intent(in) :: za, zb, Fa real(wp), intent(inout) :: Fb real(wp), intent(in) :: R real(wp), intent(out) :: g real(wp) :: norm logical, intent(in) :: switch logical, intent(in) :: lsame norm = (1d0/24d0)*sqrt(za**3*zb**5*3d0) if(switch) then Fb = -Fb end if if(lsame) then t1 = Fa * R t3 = exp(-0.5000000000_wp * t1) t6 = Fa ** 2 t7 = R ** 2 g = -0.1000000000D-8 * R * t3 * (0.5333333380D10 + 0.2666666670D10 & * t1 + 0.1333333333D10 * t6 * t7) / t6 g = g*norm else t3 = exp(-0.5000000000_wp * Fa * R) t4 = Fa ** 2 t5 = t4 * Fa t6 = Fb * R t7 = 0.5000000000_wp * t6 t8 = exp(t7) t9 = t5 * t8 t11 = Fb ** 2 t12 = t11 * Fa t15 = exp(-t7) t18 = t4 ** 2 t19 = R * t18 t22 = t11 ** 2 t29 = Fb * t4 t36 = R ** 2 t37 = t36 * t18 t44 = t36 * R t48 = -0.12D2 * t9 + 0.4D1 * t12 * t8 - 0.4D1 * t12 * t15 & - 0.6D1 * t19 * t8 - 0.6D1 * t22 * t8 * R - 0.6D1 * t22 * t15 & * R + 0.4D1 * t29 * t15 - 0.4D1 * t29 * t8 + 0.6D1 * t19 * t15 & + 0.2D1 * t37 * t8 * Fb + 0.4D1 * t37 * t15 * Fb + t44 * t18 * t15 * t11 t49 = t5 * t15 t51 = t11 * Fb t58 = t51 * Fa t59 = R * t8 t76 = t36 * t15 t79 = t22 * Fa t87 = 0.12D2 * t49 - 0.12D2 * t51 * t15 - 0.1D1 * t22 * t4 * t15 * t44 & + 0.4D1 * t58 * t59 - 0.8D1 * t58 * R * t15 + 0.4D1 * t9 * t6 + 0.8D1 * & t49 * t6 + 0.2D1 * t49 * t11 * t36 + 0.4D1 * t11 * t4 * t59 - 0.2D1 * t51 & * t4 * t76 - 0.2D1 * t79 * t36 * t8 - 0.4D1 * t79 * t76 + 0.12D2 * t51 * t8 g = -0.16D2 * t3 * (t48 + t87) / t36 / t22 / t18 g = g*norm end if end subroutine g2s1s subroutine g2s2s(za, zb, Fa, Fb, R, g, SameElement) ! slater overlap derv. ! derivative of explicit integral expression ! using maple implicit real(wp) (t) real(wp), intent(in) :: za, zb, Fa, Fb real(wp), intent(in) :: R real(wp), intent(out) :: g real(wp) :: norm logical, intent(in) :: SameElement norm = 1d0/(16d0*3d0)*SQRT((ZA*ZB)**5) if(SameElement) then t2 = R ** 2 t5 = Fa ** 2 t9 = t5 * Fa t10 = t2 ** 2 t16 = exp(-Fa * R / 0.2D1) g = (-0.4266666666D2 * R - 0.2133333333D2 * Fa * t2 - 0.2133333333D1 & * t5 * t2 * R - 0.1066666666D1 * t9 * t10) * t16 / t9 g = g*norm else t1 = R ** 2 t3 = 0.3840000000D3 * t1 * Fb t4 = t1 * R t5 = Fb ** 2 t7 = 0.6400000000D2 * t4 * t5 t8 = 0.7680000000D3 * R t10 = Fa ** 2 t11 = t10 ** 2 t12 = t11 * Fa t14 = Fb * R t15 = 0.768000000D3 * t14 t17 = 0.1280000000D3 * t5 * t1 t21 = 0.256000000D3 * t5 * R t22 = t5 * Fb t24 = 0.1280000000D3 * t22 * t1 t26 = t10 * Fa t28 = t5 ** 2 t30 = 0.1280000000D3 * t1 * t28 t32 = 0.256000000D3 * t22 * R t33 = 0.512000000D3 * t5 t34 = t28 * Fb t36 = 0.6400000000D2 * t4 * t34 t40 = 0.768000000D3 * t28 * R t42 = 0.3840000000D3 * t1 * t34 t45 = 0.1536000000D4 * t28 t47 = 0.7680000000D3 * t34 * R t51 = exp(-0.5_wp * Fa * R) t53 = 0.5_wp * t14 t54 = exp(-t53) t68 = exp(t53) g = (((t3 + t7 + t8) * t12 + (0.1536000000D4 + t15 + t17) * t11 + & (-t21 - t24) * t26 + (t30 - t32 - t33 + t36) * t10 + (t40 + t42) *& Fa + t45 + t47) * t51 * t54 + ((t3 - t8 - t7) * t12 + (-0.1536000000D4 & + t15 - t17) * t11 + (-t24 + t21) * t26 + (-t30 + t33 - t32 & + t36) * t10 + (-t40 + t42) * Fa + t47 - t45) * t51 * t68) / t1 / & t12 / t34 g = g*norm end if end subroutine g2s2s subroutine g1s3s(za, zb, Fa, Fb, R, g, switch, lsame) ! slater overlap derv. ! derivative of explicit integral expression ! using maple implicit real(wp) (t) real(wp), intent(in) :: za, zb, Fa real(wp), intent(inout) :: Fb real(wp), intent(in) :: R real(wp), intent(out) :: g real(wp) :: norm logical, intent(in) :: switch logical, intent(in) :: lsame if(switch) Fb = -Fb norm = SQRT((ZA**3)*(ZB**7)/7.5_wp)/(16d0*sqrt(3d0)) if(lsame) then t1 = Fa * R t3 = exp(-0.5000000000_wp * t1) t4 = R ** 2 g = -0.1600000000D1 * t3 * t4 * R * (0.2D1 + t1) / Fa g = g*norm else t3 = exp(-0.5000_wp * Fa * R) t4 = Fb ** 2 t5 = t4 ** 2 t6 = t5 * Fb t7 = t6 * Fa t8 = R ** 2 t9 = Fb * R t10 = 0.50_wp * t9 t11 = exp(t10) t15 = exp(-t10) t16 = t8 * t15 t19 = Fa ** 2 t21 = t8 * R t22 = t21 * t15 t25 = t19 * Fa t27 = t8 ** 2 t31 = t19 ** 2 t32 = t31 * Fa t33 = t8 * t32 t45 = t4 * Fb t48 = t31 * t15 t55 = t4 * t25 t56 = t11 * R t59 = t15 * R t62 = t5 * Fa t73 = -0.6D1 * t7 * t8 * t11 - 0.18D2 * t7 * t16 - 0.6D1 * t6 * t19 & * t22 - 0.1D1 * t6 * t25 * t27 * t15 + 0.6D1 * t33 * t11 * Fb + 0.18D2 & * t33 * t15 * Fb + 0.6D1 * t21 * t32 * t15 * t4 + t27 * t32* t15 * t45 & + 0.2D1 * t48 * t45 * t21 + 0.12D2 * t48 * t4 * t8 + 0.12D2 * t55 * t56 & + 0.12D2 * t55 * t59 + 0.12D2 * t62 * t56 - 0.36D2 * t62 * t59 - 0.12D2 & * t5 * t19 * t16 - 0.2D1 * t5 * t25 * t22 t74 = t31 * t11 t79 = t45 * t19 t92 = R * t32 t95 = t45 * Fa t100 = Fb * t25 t111 = 0.12D2 * t74 * t9 + 0.36D2 * t48 * t9 + 0.12D2 * t79 * t56 - 0.12D2 & * t79 * t59 + 0.48D2 * t5 * t11 - 0.24D2 * t6 * t11 * R - 0.24D2 * t6 * t15 & * R - 0.24D2 * t92 * t11 + 0.24D2 * t95 * t11 - 0.24D2 * t95 * t15 + 0.24D2 & * t100 * t15 + 0.24D2 * t92 * t15 - 0.24D2 * t100 * t11 - 0.48D2 * t5 * t15 & - 0.48D2 * t74 + 0.48D2 * t48 g = -0.32D2 * t3 * (t73 + t111) / t8 / t6 / t32 g = g*norm end if end subroutine g1s3s subroutine g2s3s(za, zb, Fa, Fb, R, g, switch, lsame) ! slater overlap derv. ! derivative of explicit integral expression ! using maple implicit real(wp) (t) real(wp), intent(in) :: za, zb, Fa real(wp), intent(inout) :: Fb real(wp), intent(in) :: R real(wp), intent(out) :: g real(wp) :: norm logical, intent(in) :: switch logical, intent(in) :: lsame norm = sqrt((za**5)*(zb**7)/7.5_wp)/96.0d0 if(switch) Fb = -Fb if(lsame) then t1 = Fa * R t3 = exp(-0.5000000000_wp * t1) t6 = Fa ** 2 t7 = R ** 2 t14 = t6 ** 2 t15 = t7 ** 2 g = -0.2000000000D-8 * R * t3 * (0.1280000000D12 + 0.6400000000D11 & * t1 + 0.1280000000D11 * t6 * t7 + 0.1066666670D10 * t6 * Fa * t7 & * R + 0.533333333D9 * t14 * t15) / t14 g = g*norm else t3 = exp(-0.5_wp * Fa * R) t4 = Fb ** 2 t5 = t4 ** 2 t6 = Fa ** 2 t7 = t6 * Fa t8 = t5 * t7 t9 = R ** 2 t11 = 0.50_wp * Fb * R t12 = exp(t11) t13 = t9 * t12 t16 = t6 ** 2 t17 = t16 * Fa t18 = exp(-t11) t21 = t5 * Fb t28 = t9 * t18 t32 = t9 * R t33 = t32 * t18 t36 = t5 * t4 t38 = t9 ** 2 t39 = t38 * t18 t41 = t21 * Fa t42 = R * t12 t45 = t16 * t6 t46 = t4 * Fb t49 = t46 * t16 t52 = -0.6D1 * t8 * t13 + 0.120D3 * t17 * t18 + 0.120D3 * t21 * t18 & - 0.120D3 * t17 * t12 - 0.120D3 * t21 * t12 - 0.6D1 * t8 * t28 - 0.2D1 & * t5 * t16 * t33 + t36 * t7 * t39 - 0.48D2 * t41 * t42 + t45 * t46 * t39 - 0.6D1 * t49 * t13 t54 = R * t18 t60 = t46 * t6 t63 = Fb * t16 t66 = t5 * Fa t69 = t4 * t7 t72 = t36 * t6 t75 = t32 * t12 t78 = Fb * t9 t84 = Fb * t17 t87 = -0.24D2 * t46 * t7 * t54 - 0.24D2 * t5 * t6 * t42 + 0.24D2 *& t60 * t12 + 0.24D2 * t63 * t18 - 0.24D2 * t66 * t12 - 0.24D2 * t69 & * t18 + 0.9D1 * t72 * t33 + 0.3D1 * t72 * t75 + 0.24D2 * t78 * t45 & * t12 - 0.6D1 * t49 * t28 + 0.48D2 * t84 * t42 t102 = t21 * t6 t105 = t4 * t17 t113 = t45 * t4 t118 = 0.72D2 * t84 * t54 + 0.72D2 * t41 * t54 + 0.36D2 * t78 * t45 & * t18 + 0.2D1 * t46 * t17 * t33 + 0.24D2 * t4 * t16 * t42 - 0.6D1 & * t102 * t13 - 0.6D1 * t105 * t13 + 0.18D2 * t105 * t28 + 0.2D1 & * t21 * t7 * t33 - 0.3D1 * t113 * t75 + 0.9D1 * t113 * t33 t121 = t36 * Fa t130 = R * t45 t145 = 0.18D2 * t102 * t28 + 0.24D2 * t121 * t13 + 0.36D2 * t121 * & t28 - 0.24D2 * t60 * t18 - 0.24D2 * t63 * t12 + 0.60D2 * t130 * t18 & + 0.60D2 * t36 * t18 * R + 0.24D2 * t69 * t12 + 0.60D2 * t36 * t12 * & R - 0.60D2 * t130 * t12 + 0.24D2 * t66 * t18 g = 0.128D3 * t3 * (t52 + t87 + t118 + t145) / t9 / t36 / t45 g = g*norm end if end subroutine g2s3s subroutine g3s3s(za, zb, Fa, Fb, R, g, SameElement) ! slater overlap derv. ! derivative of explicit integral expression ! using maple implicit real(wp) (t) real(wp), intent(in) :: za, zb, Fa, Fb real(wp), intent(in) :: R real(wp), intent(out) :: g real(wp) :: norm logical, intent(in) :: SameElement norm = sqrt((ZA*ZB)**7)/1440.0d0 if(SameElement) then t1 = Fa * R t3 = exp(-0.5000000000_wp * t1) t5 = Fa ** 2 t6 = t5 ** 2 t7 = t6 * Fa t8 = R ** 2 t9 = t8 ** 2 g = -0.2000000000D-8 * t3 * R * (0.457142857D9 * t7 * t9 & * R + 0.7680000000D12 * t1 + 0.1536000000D12 * t5 * t8 & + 0.1280000000D11 * t5 * Fa * t8 * R + 0.914285715D9 * t6 * t9 + 0.1536000000D13) / t7 g = g*norm else t3 = exp(-0.5000000000_wp * Fa * R) t4 = Fa ** 2 t5 = t4 ** 2 t6 = t5 * t4 t7 = Fb * R t8 = 0.5000000000_wp * t7 t9 = exp(-t8) t10 = t6 * t9 t13 = Fb ** 2 t14 = t13 * Fb t15 = t13 ** 2 t16 = t15 * t14 t17 = R ** 2 t18 = t17 * R t19 = t16 * t18 t23 = exp(t8) t24 = t6 * t23 t27 = t5 * Fa t28 = t27 * t13 t29 = R * t23 t32 = t6 * t13 t33 = t17 * t9 t36 = t15 * Fb t37 = t4 * t36 t38 = t9 * R t43 = t17 * t23 t46 = t4 * Fa t47 = t5 * t46 t48 = t47 * t18 t52 = t47 * t17 t65 = 0.120D3 * t10 * t7 - 0.12D2 * t19 * t4 * t9 + 0.120D3 & * t24 * t7 + 0.24D2 * t28 * t29 + 0.24D2 * t32 * t33 + 0.24D2 * t37 & * t38 - 0.24D2 * t28 * t38 - 0.24D2 * t32 * t43 - 0.12D2 * t48 * t13 & * t23 + 0.60D2 * t52 * t23 * Fb + 0.12D2 * t48 * t13 * t9 + 0.60D2 & * t52 * t9 * Fb - 0.12D2 * t19 * t4 * t23 t66 = t17 ** 2 t67 = t16 * t66 t74 = t27 * t14 t77 = t6 * t14 t78 = t18 * t23 t81 = t46 * t15 t86 = t27 * t15 t89 = t5 * t36 t90 = t18 * t9 t97 = t46 * t36 t104 = -0.1D1 * t67 * t46 * t9 - 0.1D1 * t67 * t46 * t23 - 0.12D2 & * t74 * t43 + 0.2D1 * t77 * t78 - 0.24D2 * t81 * t29 + 0.24D2 * t81 & * t38 + 0.2D1 * t86 * t78 + 0.2D1 * t89 * t90 - 0.2D1 * t86 * t90 & + 0.24D2 * t37 * t29 + 0.12D2 * t97 * t33 + 0.2D1 * t89 * t78 - 0.12D2 * t74 * t33 t108 = t5 * t14 t111 = t15 * t13 t112 = t111 * t4 t117 = t111 * t46 t122 = t111 * Fa t129 = t4 * t15 t132 = t47 * R t139 = 0.2D1 * t77 * t90 - 0.24D2 * t108 * t38 + 0.24D2 * t112 * t43 & - 0.24D2 * t112 * t33 + 0.2D1 * t117 * t78 - 0.2D1 * t117 * t90 + 0.120D3 & * t122 * t29 - 0.120D3 * t122 * t38 + 0.12D2 * t97 * t43 - 0.48D2 * t129 & * t23 + 0.120D3 * t132 * t9 - 0.120D3 * t132 * t23 + 0.240D3 * t111 * t23 t140 = t47 * t66 t145 = t16 * R t150 = t16 * t17 t160 = t5 * t13 t170 = t140 * t14 * t23 + t140 * t14 * t9 - 0.120D3 * t145 * t9 - 0.24D2 & * t108 * t29 - 0.60D2 * t150 * Fa * t23 - 0.240D3 * t111 * t9 - 0.240D3 & * t24 + 0.240D3 * t10 + 0.48D2 * t129 * t9 - 0.48D2 * t160 * t9 + 0.48D2 & * t160 * t23 - 0.120D3 * t145 * t23 - 0.60D2 * t150 * Fa * t9 g = -0.768D3 * t3 * (t65 + t104 + t139 + t170) / t17 / t47 / t16 g = g*norm end if end subroutine g3s3s end module dftd3_gcp