add_coordination_number_derivs Subroutine

subroutine add_coordination_number_derivs(mol, trans, cutoff, rcov, dEdcn, gradient, &
      & sigma, partition)

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

   !> Lattice points
   real(wp), intent(in) :: trans(:, :)

   !> Real space cutoff
   real(wp), intent(in) :: cutoff

   !> Covalent radius
   real(wp), intent(in) :: rcov(:)

   !> Derivative of expression with respect to the coordination number
   real(wp), intent(in) :: dEdcn(:)

   !> Derivative of the CN with respect to the Cartesian coordinates
   real(wp), intent(inout) :: gradient(:, :)

   !> Derivative of the CN with respect to strain deformations
   real(wp), intent(inout) :: sigma(:, :)

   !> Work partition of the pair loop, absent selects the complete work
   type(work_partition), intent(in), optional :: partition

   integer :: iat, jat, izp, jzp, itr
   real(wp) :: vec(3), r2, r1, rc, cutoff2, expterm, cf, dcf, countd(3), ds(3, 3)
   real(wp), allocatable :: gradient_local(:, :), sigma_local(:, :)

   cutoff2 = cutoff*cutoff

   !$omp parallel default(none) &
   !$omp shared(mol, trans, cutoff2, rcov, dEdcn, gradient, sigma, partition) &
   !$omp private(iat, jat, itr, izp, jzp, vec, r2, r1, rc, expterm, cf, dcf) &
   !$omp private(countd, ds, gradient_local, sigma_local)
   allocate(gradient_local(size(gradient, 1), size(gradient, 2)), source=0.0_wp)
   allocate(sigma_local(size(sigma, 1), size(sigma, 2)), source=0.0_wp)
   !$omp do schedule(runtime)
   do iat = 1, mol%nat
      izp = mol%id(iat)
      do jat = 1, iat
         if (.not.owns_pair(partition, iat, jat)) cycle
         jzp = mol%id(jat)
         rc = rcov(izp) + rcov(jzp)
         do itr = 1, size(trans, 2)
            vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, itr))
            r2 = vec(1)*vec(1) + vec(2)*vec(2) + vec(3)*vec(3)
            if (r2 > cutoff2 .or. r2 < 1.0e-12_wp) cycle
            r1 = sqrt(r2)

            expterm = exp(-default_kcn*(rc/r1 - 1.0_wp))
            cf = 1.0_wp/(1.0_wp + expterm)
            dcf = -cf*(1.0_wp - cf)*default_kcn*rc/r2
            countd(:) = dcf * vec/r1

            gradient_local(:, iat) = gradient_local(:, iat) &
               & + countd*(dEdcn(iat) + dEdcn(jat))
            gradient_local(:, jat) = gradient_local(:, jat) &
               & - countd*(dEdcn(iat) + dEdcn(jat))

            ds(:, :) = spread(countd, 1, 3) * spread(vec, 2, 3)
            sigma_local(:, :) = sigma_local(:, :) &
               & + ds*(dEdcn(iat) + merge(dEdcn(jat), 0.0_wp, jat /= iat))
         end do
      end do
   end do
   !$omp end do
   !$omp critical (add_coordination_number_derivs_)
   gradient(:, :) = gradient(:, :) + gradient_local(:, :)
   sigma(:, :) = sigma(:, :) + sigma_local(:, :)
   !$omp end critical (add_coordination_number_derivs_)
   deallocate(gradient_local, sigma_local)
   !$omp end parallel

end subroutine add_coordination_number_derivs