weight_references Subroutine

subroutine weight_references(self, mol, cn, q, gwvec, gwdcn, gwdq)
   !DEC$ ATTRIBUTES DLLEXPORT :: weight_references

   !> Instance of the dispersion model
   class(d4_model), intent(in) :: self

   !> Molecular structure data
   class(structure_type), intent(in) :: mol

   !> Coordination number of every atom
   real(wp), intent(in) :: cn(:)

   !> Partial charge of every atom
   real(wp), intent(in) :: q(:)

   !> weighting for the atomic reference systems
   real(wp), intent(out) :: gwvec(:, :, :)

   !> derivative of the weighting function w.r.t. the coordination number
   real(wp), intent(out), optional :: gwdcn(:, :, :)

   !> derivative of the weighting function w.r.t. the charge scaling
   real(wp), intent(out), optional :: gwdq(:, :, :)

   integer :: iat, izp, iref, igw
   real(wp) :: norm, dnorm, gw, expw, expd, gwk, dgwk, wf, zi, gi, maxcn

   if (present(gwdcn) .and. present(gwdq)) then
      gwvec(:, :, :) = 0.0_wp
      gwdcn(:, :, :) = 0.0_wp
      gwdq(:, :, :) = 0.0_wp

      !$omp parallel do default(none) schedule(runtime) &
      !$omp shared(gwvec, gwdcn, gwdq, mol, self, cn, q) private(iat, izp, iref, &
      !$omp& igw, norm, dnorm, gw, expw, expd, gwk, dgwk, wf, zi, gi, maxcn)
      do iat = 1, mol%nat
         izp = mol%id(iat)
         zi = self%zeff(izp)
         gi = self%eta(izp) * self%gc
         norm = 0.0_wp
         dnorm = 0.0_wp
         do iref = 1, self%ref(izp)
            do igw = 1, self%ngw(iref, izp)
               wf = igw * self%wf
               gw = weight_cn(wf, cn(iat), self%cn(iref, izp))
               norm = norm + gw
               dnorm = dnorm + 2*wf * (self%cn(iref, izp) - cn(iat)) * gw
            end do
         end do
         norm = 1.0_wp / norm
         do iref = 1, self%ref(izp)
            expw = 0.0_wp
            expd = 0.0_wp
            do igw = 1, self%ngw(iref, izp)
               wf = igw * self%wf
               gw = weight_cn(wf, cn(iat), self%cn(iref, izp))
               expw = expw + gw
               expd = expd + 2*wf * (self%cn(iref, izp) - cn(iat)) * gw
            end do
            gwk = expw * norm
            if (is_exceptional(gwk)) then
               maxcn = maxval(self%cn(:self%ref(izp), izp))
               if (abs(maxcn - self%cn(iref, izp)) < 1e-12_wp) then
                  gwk = 1.0_wp
               else
                  gwk = 0.0_wp
               end if
            end if
            gwvec(iref, iat, 1) = gwk * zeta(self%ga, gi, self%q(iref, izp)+zi, q(iat)+zi)
            gwdq(iref, iat, 1) = gwk * dzeta(self%ga, gi, self%q(iref, izp)+zi, q(iat)+zi)

            dgwk = norm * (expd - expw * dnorm * norm)
            if (is_exceptional(dgwk)) then
               dgwk = 0.0_wp
            end if
            gwdcn(iref, iat, 1) = dgwk * zeta(self%ga, gi, self%q(iref, izp)+zi, q(iat)+zi)
         end do
      end do

   else

      gwvec(:, :, :) = 0.0_wp

      !$omp parallel do default(none) schedule(runtime) &
      !$omp shared(gwvec, mol, self, cn, q) &
      !$omp private(iat, izp, iref, igw, norm, gw, expw, gwk, wf, zi, gi, maxcn)
      do iat = 1, mol%nat
         izp = mol%id(iat)
         zi = self%zeff(izp)
         gi = self%eta(izp) * self%gc
         norm = 0.0_wp
         do iref = 1, self%ref(izp)
            do igw = 1, self%ngw(iref, izp)
               wf = igw * self%wf
               norm = norm + weight_cn(wf, cn(iat), self%cn(iref, izp))
            end do
         end do
         norm = 1.0_wp / norm
         do iref = 1, self%ref(izp)
            expw = 0.0_wp
            do igw = 1, self%ngw(iref, izp)
               wf = igw * self%wf
               expw = expw + weight_cn(wf, cn(iat), self%cn(iref, izp))
            end do
            gwk = expw * norm
            if (is_exceptional(gwk)) then
               maxcn = maxval(self%cn(:self%ref(izp), izp))
               if (abs(maxcn - self%cn(iref, izp)) < 1e-12_wp) then
                  gwk = 1.0_wp
               else
                  gwk = 0.0_wp
               end if
            end if
            gwvec(iref, iat, 1) = gwk * zeta(self%ga, gi, self%q(iref, izp)+zi, q(iat)+zi)
         end do
      end do
   end if

end subroutine weight_references