ncoord_erf Subroutine

private subroutine ncoord_erf(mol, trans, cutoff, rcov, en, cn)

Arguments

Type IntentOptional Attributes Name
type(structure_type), intent(in) :: mol

Molecular structure data

real(kind=wp), intent(in) :: trans(:,:)

Lattice points

real(kind=wp), intent(in) :: cutoff

Real space cutoff

real(kind=wp), intent(in) :: rcov(:)

Covalent radius

real(kind=wp), intent(in) :: en(:)

Electronegativity

real(kind=wp), intent(out) :: cn(:)

Error function coordination number.


Source Code

subroutine ncoord_erf(mol, trans, cutoff, rcov, en, cn)

   !> 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(:)

   !> Electronegativity
   real(wp), intent(in) :: en(:)

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

   integer :: iat, jat, izp, jzp, itr
   real(wp) :: r2, r1, rc, rij(3), countf, cutoff2, den

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

   !$omp parallel do default(none) reduction(+:cn) &
   !$omp shared(mol, trans, cutoff2, rcov, en) &
   !$omp private(jat, itr, izp, jzp, r2, rij, r1, rc, countf, den)
   do iat = 1, mol%nat
      izp = mol%id(iat)
      do jat = 1, iat
         jzp = mol%id(jat)
         den = k4*exp(-(abs(en(izp)-en(jzp)) + k5)**2/k6)

         do itr = 1, size(trans, dim=2)
            rij = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, itr))
            r2 = sum(rij**2)
            if (r2 > cutoff2 .or. r2 < 1.0e-12_wp) cycle
            r1 = sqrt(r2)

            rc = rcov(izp) + rcov(jzp)

            countf = den*erf_count(kcn, r1, rc)

            cn(iat) = cn(iat) + countf
            if (iat /= jat) then
               cn(jat) = cn(jat) + countf
            end if

         end do
      end do
   end do

end subroutine ncoord_erf