weight_references Subroutine

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

   !> Instance of the dispersion model
   class(d4s_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(:)

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

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

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

   integer :: iat, izp, iref, igw, jat, jzp
   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) &
      !$omp private(iat, izp, iref, igw, zi, gi, jat, jzp) &
      !$omp private(norm, dnorm, gw, expw, expd, gwk, dgwk, wf, maxcn)      
      do iat = 1, mol%nat
         izp = mol%id(iat)
         zi = self%zeff(izp)
         gi = self%eta(izp) * self%gc

         do jat = 1, mol%nat
            jzp = mol%id(jat)

            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(izp, jzp)
                  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(izp, jzp)
                  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, jat) = gwk * zeta(self%ga, gi, self%q(iref, izp)+zi, q(iat)+zi)
               gwdq(iref, iat, jat) = 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, jat) = dgwk * zeta(self%ga, gi, self%q(iref, izp)+zi, q(iat)+zi)
            end do

         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, zi, gi, jat, jzp) &
      !$omp private(norm, gw, expw, gwk, wf, maxcn)
      do iat = 1, mol%nat
         izp = mol%id(iat)
         zi = self%zeff(izp)
         gi = self%eta(izp) * self%gc

         do jat = 1, mol%nat
            jzp = mol%id(jat)

            norm = 0.0_wp
            do iref = 1, self%ref(izp)
               do igw = 1, self%ngw(iref, izp)
                  wf = igw * self%wf(izp, jzp)
                  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(izp, jzp)
                  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, jat) = gwk * zeta(self%ga, gi, self%q(iref, izp)+zi, q(iat)+zi)
            end do

         end do 
      end do
   end if

end subroutine weight_references