get_partitioned_coordination_number Subroutine

subroutine get_partitioned_coordination_number(mol, trans, cutoff, rcov, cn, 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(:)

   !> Error function coordination number.
   real(wp), intent(out) :: cn(:)

   !> 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, countf
   real(wp), allocatable :: cn_local(:)

   cn(:) = 0.0_wp
   cutoff2 = cutoff*cutoff

   !$omp parallel default(none) &
   !$omp shared(mol, trans, cutoff2, rcov, cn, partition) &
   !$omp private(iat, jat, itr, izp, jzp, vec, r2, r1, rc, countf, cn_local)
   allocate(cn_local(size(cn)), 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)

            countf = 1.0_wp/(1.0_wp + exp(-default_kcn*(rc/r1 - 1.0_wp)))

            cn_local(iat) = cn_local(iat) + countf
            if (iat /= jat) cn_local(jat) = cn_local(jat) + countf
         end do
      end do
   end do
   !$omp end do
   !$omp critical (get_partitioned_coordination_number_)
   cn(:) = cn(:) + cn_local(:)
   !$omp end critical (get_partitioned_coordination_number_)
   deallocate(cn_local)
   !$omp end parallel

end subroutine get_partitioned_coordination_number