weight_references_hessian Subroutine

subroutine weight_references_hessian(self, mol, cn, q, gwvec, gwdcn, gwd2cn)

   !> 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 and its CN derivatives
   real(wp), intent(out) :: gwvec(:, :, :), gwdcn(:, :, :), gwd2cn(:, :, :)

   integer :: iat, izp, iref, igw
   real(wp) :: norm, dnorm, d2norm, gw, expw, expd, expdd
   real(wp) :: gwk, dgwk, d2gwk, wf, dcn, zi, gi, zscale, maxcn
   real(wp), parameter :: eps_norm = tiny(1.0_wp)**0.5_wp

   gwvec(:, :, :) = 0.0_wp
   gwdcn(:, :, :) = 0.0_wp
   gwd2cn(:, :, :) = 0.0_wp

   !$omp parallel do default(none) schedule(runtime) &
   !$omp shared(gwvec, gwdcn, gwd2cn, mol, self, cn, q) &
   !$omp private(iat, izp, iref, igw, norm, dnorm, d2norm, gw, expw, expd, expdd, &
   !$omp& gwk, dgwk, d2gwk, wf, dcn, zi, gi, zscale, 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
      d2norm = 0.0_wp
      do iref = 1, self%ref(izp)
         do igw = 1, self%ngw(iref, izp)
            wf = igw * self%wf
            dcn = self%cn(iref, izp) - cn(iat)
            gw = weight_cn(wf, cn(iat), self%cn(iref, izp))
            norm = norm + gw
            dnorm = dnorm + 2.0_wp*wf*dcn*gw
            d2norm = d2norm + (4.0_wp*wf*wf*dcn*dcn - 2.0_wp*wf)*gw
         end do
      end do

      if (abs(norm) > eps_norm) then
         norm = 1.0_wp / norm
      else
         norm = 0.0_wp
      end if

      do iref = 1, self%ref(izp)
         expw = 0.0_wp
         expd = 0.0_wp
         expdd = 0.0_wp
         do igw = 1, self%ngw(iref, izp)
            wf = igw * self%wf
            dcn = self%cn(iref, izp) - cn(iat)
            gw = weight_cn(wf, cn(iat), self%cn(iref, izp))
            expw = expw + gw
            expd = expd + 2.0_wp*wf*dcn*gw
            expdd = expdd + (4.0_wp*wf*wf*dcn*dcn - 2.0_wp*wf)*gw
         end do

         gwk = expw*norm
         if (is_exceptional(gwk) .or. norm == 0.0_wp) 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

         dgwk = norm*(expd - expw*dnorm*norm)
         if (is_exceptional(dgwk) .or. norm == 0.0_wp) dgwk = 0.0_wp

         d2gwk = expdd*norm - 2.0_wp*expd*dnorm*norm**2 &
            & - expw*d2norm*norm**2 + 2.0_wp*expw*dnorm**2*norm**3
         if (is_exceptional(d2gwk) .or. norm == 0.0_wp) d2gwk = 0.0_wp

         zscale = zeta(self%ga, gi, self%q(iref, izp)+zi, q(iat)+zi)
         gwvec(iref, iat, 1) = gwk*zscale
         gwdcn(iref, iat, 1) = dgwk*zscale
         gwd2cn(iref, iat, 1) = d2gwk*zscale
      end do
   end do

end subroutine weight_references_hessian