! 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 mctc_env, only : wp use mctc_io, only : structure_type use dftd3_gcp_param, only : gcp_param, get_gcp_param use dftd3_cutoff, only : realspace_cutoff, get_lattice_points implicit none private public :: gcp_param, get_gcp_param, get_geometric_counterpoise !> Geometric counterpoise correction interface get_geometric_counterpoise module procedure get_geometric_counterpoise module procedure get_geometric_counterpoise_atomic end interface get_geometric_counterpoise contains !> Geometric counterpoise correction subroutine get_geometric_counterpoise(mol, param, cutoff, energy, gradient, sigma) !> 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(:, :) real(wp), allocatable :: energies(:) allocate(energies(mol%nat), source=0.0_wp) call get_geometric_counterpoise_atomic(mol, param, cutoff, energies, gradient, sigma) 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) !> 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(:, :) real(wp), allocatable :: lattr(:, :) real(wp), parameter :: rexp_base = 0.75_wp, zexp_base = 1.5_wp, rexp_srb = -1.0_wp, zexp_srb = 0.5_wp 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) 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) 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) else call srb_energy(mol, lattr, cutoff%srb, mol%num, param%rvdw_srb, param%rscal, param%qscal, & & rexp_srb, zexp_srb, energies) 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) else call srb_energy(mol, lattr, cutoff%srb, param%zeff, param%rvdw, param%rscal, param%qscal, & & rexp_base, zexp_base, energies) end if end if end subroutine get_geometric_counterpoise_atomic !> Geometric counterpoise correction subroutine gcp_energy(mol, trans, cutoff, iz, emiss, slater, xv, rvdw, escal, alpha, beta, & & damp, dmp_scal, dmp_exp, energies) !> 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(:) 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 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 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 endif 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) !> 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(:, :) 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) 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 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 endif 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 endif dE = bsse*dampval dG = expd*(argd*vec + gs)/ovlpd*emij if(damp) then dG = dG*dampval+bsse*grd_dmp*(vec/r1)*emij endif 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) !> 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(:) real(wp) :: fi, fj, ff, r1, expt real(wp) :: r0, vec(3), dE integer iat, jat, jtr, izp, jzp do iat = 1, mol%nat izp = mol%id(iat) fi = real(iz(izp), wp) do jat = 1, iat 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) !> 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(:, :) 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 do iat = 1, mol%nat izp = mol%id(iat) fi = real(iz(izp), wp) do jat = 1, iat 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 !****************************************************************************** !* 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) integer ii, shell(72) logical debug real(wp) za, zb, R, ovl, ax, bx, norm, R05 integer na, nb real(wp) Bxx0, Bxx1, Bxx2, xx, Bxx4, Bxx6 real(wp) Bxx3, Bxx5 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) xza, xzb integer iat, jat, iz(*) 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._wp/3._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)) endif 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.d0*(A3(ax)*Bint(bx, 1)-Bint(bx, 3)*A1(ax)))/sqrt(3.d0) 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.d0*(A3(ax)*Bint(bx, 1)-Bint(bx, 3)*A1(ax)))/sqrt(3.d0) endif 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.d0 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.d0 endif case(9) norm = sqrt((ZA*ZB*R*R)**7)/480.d0 ovl = norm*(A6(ax)*Bint(bx, 0)-3.d0*(A4(ax)*Bint(bx, 2) & & -A2(ax)*Bint(bx, 4))-A0(ax)*Bint(bx, 6))/3._wp 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._wp/3._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)) endif 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.d0*(A3(ax)*B1(bx)-B3(bx)*A1(ax)))/sqrt(3.d0) 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.d0*(A3(ax)*B1(bx)-B3(bx)*A1(ax)))/sqrt(3.d0) endif 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.d0 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.d0 endif case(9) ! <3s|3> norm = sqrt((ZA*ZB*R*R)**7)/1440.d0 ovl = norm*(A6(ax)*B0(bx)-3.d0*(A4(ax)*B2(bx)-A2(ax)*B4(bx))-A0(ax)*B6(bx)) end select endif 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 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 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 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 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 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 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 !************************************** !* 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 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 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 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 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 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 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 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).lt.1e-6) then do i = 0, k bint = (1.d0+(-1d0)**i)/(dble(i)+1.d0) end do return endif do i = 0, 12 xx = 1d0-((-1d0)**(k+i+1)) yy = dble(fact(i))*dble((k+i+1)) bint = bint+xx/yy*(-x)**i enddo end function bint ! faculty function integer(wp) function fact(N) implicit none integer j, n fact = 1 do j = 2, n fact = fact*j enddo return end 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 integer ii, shell(72) logical debug real(wp) ax, bx, R05, za, zb, R 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) g, Fa, Fb !--------------------- set exponents --------------------------------------- real(wp) xza, xzb real(wp) xx integer iat, jat, iz(*) 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. ! ! same elements if(abs(za-zb) < 0.1) then lsame = .true. ! call arguments: gtype(exponents, argumentDeriv., distance, gradient, (Switch shell), sameElement) select case (ii) case (1) call g1s1s(za, zb, Fa, Fb, R, g, lsame) case (2) 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) call g2s2s(za, zb, Fa, Fb, R, g, lsame) case(3) 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) 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) call g3s3s(za, zb, Fa, Fb, R, g, lsame) end select else ! different elements lsame = .false. 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) end select end if end subroutine gsovl !------------------------------------------------------------- ! 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) za, zb, Fa, Fb real(wp) g, R logical 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) za, zb, Fa, Fb real(wp) g, R, norm logical switch logical lsame norm = (1d0/24d0)*sqrt(za**3*zb**5*3d0) if(switch) then Fb = -Fb endif 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) za, zb, Fa, Fb real(wp) g, R, norm logical 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 endif 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) za, zb, Fa, Fb real(wp) g, R, norm logical switch logical 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 endif 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) za, zb, Fa, Fb real(wp) g, R, norm logical switch logical lsame norm = sqrt((za**5)*(zb**7)/7.5_wp)/96.d0 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 endif 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) za, zb, Fa, Fb real(wp) g, R, norm logical SameElement norm = sqrt((ZA*ZB)**7)/1440.d0 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 endif end subroutine g3s3s end module dftd3_gcp