weight_references_hessian Subroutine

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

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

   integer :: iat, izp, iref, igw, jat, jzp
   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, jat, jzp, 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
      do jat = 1, mol%nat
         jzp = mol%id(jat)
         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(izp, jzp)
               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(izp, jzp)
               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, jat) = gwk*zscale
            gwdcn(iref, iat, jat) = dgwk*zscale
            gwd2cn(iref, iat, jat) = d2gwk*zscale
         end do
      end do
   end do

end subroutine weight_references_hessian