Moving the coordination number to mctc-lib

src/mctc/CMakeLists.txt

@@ -12,15 +12,20 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+add_subdirectory("data")
 add_subdirectory("env")
 add_subdirectory("io")
+add_subdirectory("ncoord")
 
 set(dir "${CMAKE_CURRENT_SOURCE_DIR}")
 
 list(
   APPEND srcs
+  "${dir}/cutoff.f90"
+  "${dir}/data.f90"
   "${dir}/env.f90"
   "${dir}/io.f90"
+  "${dir}/ncoord.f90"
   "${dir}/version.F90"
 )
 

src/mctc/cutoff.f90

@@ -0,0 +1,198 @@
+! This file is part of mctc-lib.
+!
+! Licensed under the Apache License, Version 2.0 (the "License");
+! you may not use this file except in compliance with the License.
+! You may obtain a copy of the License at
+!
+!     http://www.apache.org/licenses/LICENSE-2.0
+!
+! Unless required by applicable law or agreed to in writing, software
+! distributed under the License is distributed on an "AS IS" BASIS,
+! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+! See the License for the specific language governing permissions and
+! limitations under the License.
+
+!> @file mctc/cutoff.f90
+!> Provides realspace cutoff and lattice generation
+
+!> Realspace cutoff and lattice point generator utilities
+module mctc_cutoff
+   use mctc_env, only : wp
+   use mctc_io_math, only : matinv_3x3
+   implicit none
+   private
+
+   public :: get_lattice_points, wrap_to_central_cell
+
+
+   interface get_lattice_points
+      module procedure :: get_lattice_points_cutoff
+      module procedure :: get_lattice_points_rep_3d
+   end interface get_lattice_points
+
+
+contains
+
+
+!> Generate lattice points from repeatitions
+subroutine get_lattice_points_rep_3d(lat, rep, origin, trans)
+
+   !> Lattice vectors
+   real(wp), intent(in) :: lat(:, :)
+
+   !> Repeatitions of lattice points to generate
+   integer, intent(in) :: rep(:)
+
+   !> Include the origin in the generated lattice points
+   logical, intent(in) :: origin
+
+   !> Generated lattice points
+   real(wp), allocatable, intent(out) :: trans(:, :)
+
+   integer :: itr, ix, iy, iz, jx, jy, jz
+
+   itr = 0
+   if (origin) then
+      allocate(trans(3, product(2*rep+1)))
+      do ix = 0, rep(1)
+         do iy = 0, rep(2)
+            do iz = 0, rep(3)
+               do jx = 1, merge(-1, 1, ix > 0), -2
+                  do jy = 1, merge(-1, 1, iy > 0), -2
+                     do jz = 1, merge(-1, 1, iz > 0), -2
+                        itr = itr + 1
+                        trans(:, itr) = lat(:, 1)*ix*jx &
+                           & + lat(:, 2)*iy*jy + lat(:, 3)*iz*jz
+                     end do
+                  end do
+               end do
+            end do
+         end do
+      end do
+   else
+      allocate(trans(3, product(2*rep+1)-1))
+      do ix = 0, rep(1)
+         do iy = 0, rep(2)
+            do iz = 0, rep(3)
+               if (ix == 0 .and. iy == 0 .and. iz == 0) cycle
+               do jx = 1, merge(-1, 1, ix > 0), -2
+                  do jy = 1, merge(-1, 1, iy > 0), -2
+                     do jz = 1, merge(-1, 1, iz > 0), -2
+                        itr = itr + 1
+                        trans(:, itr) = lat(:, 1)*ix*jx &
+                           & + lat(:, 2)*iy*jy + lat(:, 3)*iz*jz
+                     end do
+                  end do
+               end do
+            end do
+         end do
+      end do
+   end if
+
+end subroutine get_lattice_points_rep_3d
+
+
+!> Create lattice points within a given cutoff
+subroutine get_lattice_points_cutoff(periodic, lat, rthr, trans)
+
+   !> Periodic dimensions
+   logical, intent(in) :: periodic(:)
+
+   !> Real space cutoff
+   real(wp), intent(in) :: rthr
+
+   !> Lattice parameters
+   real(wp), intent(in) :: lat(:, :)
+
+   !> Generated lattice points
+   real(wp), allocatable, intent(out) :: trans(:, :)
+
+   integer :: rep(3)
+
+   if (.not.any(periodic)) then
+      allocate(trans(3, 1))
+      trans(:, :) = 0.0_wp
+   else
+      call get_translations(lat, rthr, rep)
+      call get_lattice_points(lat, rep, .true., trans)
+   end if
+
+end subroutine get_lattice_points_cutoff
+
+
+!> Generate a supercell based on a realspace cutoff, this subroutine
+!> doesn't know anything about the convergence behaviour of the
+!> associated property.
+pure subroutine get_translations(lat, rthr, rep)
+   real(wp), intent(in) :: rthr
+   real(wp), intent(in) :: lat(3, 3)
+   integer, intent(out) :: rep(3)
+   real(wp) :: normx(3), normy(3), normz(3)
+   real(wp) :: cos10, cos21, cos32
+
+   ! find normal to the plane...
+   call crossproduct(lat(:, 2), lat(:, 3), normx)
+   call crossproduct(lat(:, 3), lat(:, 1), normy)
+   call crossproduct(lat(:, 1), lat(:, 2), normz)
+   ! ...normalize it...
+   normx = normx/norm2(normx)
+   normy = normy/norm2(normy)
+   normz = normz/norm2(normz)
+   ! cos angles between normals and lattice vectors
+   cos10 = sum(normx*lat(:, 1))
+   cos21 = sum(normy*lat(:, 2))
+   cos32 = sum(normz*lat(:, 3))
+   rep(1) = ceiling(abs(rthr/cos10))
+   rep(2) = ceiling(abs(rthr/cos21))
+   rep(3) = ceiling(abs(rthr/cos32))
+
+contains
+
+   pure subroutine crossproduct(a, b, c)
+      real(wp), intent(in) :: a(3)
+      real(wp), intent(in) :: b(3)
+      real(wp), intent(out) :: c(3)
+      c(1)=a(2)*b(3)-b(2)*a(3)
+      c(2)=a(3)*b(1)-b(3)*a(1)
+      c(3)=a(1)*b(2)-b(1)*a(2)
+   end subroutine crossproduct
+
+end subroutine get_translations
+
+
+subroutine wrap_to_central_cell(xyz, lattice, periodic)
+   real(wp), intent(inout) :: xyz(:, :)
+   real(wp), intent(in) :: lattice(:, :)
+   logical, intent(in) :: periodic(:)
+   real(wp) :: invlat(3, 3), vec(3)
+   integer :: iat, idir
+
+   if (.not.any(periodic)) return
+
+   invlat = matinv_3x3(lattice)
+   do iat = 1, size(xyz, 2)
+      vec(:) = matmul(invlat, xyz(:, iat))
+      vec(:) = shift_back_abc(vec)
+      xyz(:, iat) = matmul(lattice, vec)
+   end do
+
+end subroutine wrap_to_central_cell
+
+
+elemental function shift_back_abc(in) result(out)
+   !> fractional coordinate in (-∞,+∞)
+   real(wp),intent(in) :: in
+   !> fractional coordinate in [0,1)
+   real(wp) :: out
+   real(wp),parameter :: p_pbc_eps = 1.0e-14_wp
+   out = in
+   if(in < (0.0_wp - p_pbc_eps)) &
+      out = in + real(ceiling(-in),wp)
+   if(in > (1.0_wp + p_pbc_eps)) &
+      out = in - real(floor  ( in),wp)
+   if (abs(in - 1.0_wp) < p_pbc_eps) &
+      out = in - 1.0_wp
+end function shift_back_abc
+
+
+end module mctc_cutoff

src/mctc/data.f90

@@ -0,0 +1,30 @@
+! This file is part of mctc-lib.
+!
+! Licensed under the Apache License, Version 2.0 (the "License");
+! you may not use this file except in compliance with the License.
+! You may obtain a copy of the License at
+!
+!     http://www.apache.org/licenses/LICENSE-2.0
+!
+! Unless required by applicable law or agreed to in writing, software
+! distributed under the License is distributed on an "AS IS" BASIS,
+! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+! See the License for the specific language governing permissions and
+! limitations under the License.
+
+!> @dir mctc/data
+!> Contains element data used for defining interactions.
+
+!> @file mctc/data.f90
+!> Reexports access to element-specific data.
+
+!> Proxy module for providing access to element data.
+module mctc_data
+   use mctc_data_atomicrad, only : get_atomic_rad
+   use mctc_data_covrad, only : get_covalent_rad
+   use mctc_data_paulingen, only : get_pauling_en
+   use mctc_data_vdwrad, only : get_vdw_rad
+   implicit none
+
+   public :: get_atomic_rad, get_covalent_rad, get_pauling_en, get_vdw_rad
+end module mctc_data

src/mctc/data/CMakeLists.txt

@@ -0,0 +1,25 @@
+# This file is part of mctc-lib.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+set(dir "${CMAKE_CURRENT_SOURCE_DIR}")
+
+list(
+  APPEND srcs
+  "${dir}/atomicrad.f90"
+  "${dir}/covrad.f90"
+  "${dir}/paulingen.f90"
+  "${dir}/vdwrad.f90"
+)
+
+set(srcs "${srcs}" PARENT_SCOPE)

src/mctc/data/atomicrad.f90

@@ -0,0 +1,97 @@
+! This file is part of mctc-lib.
+!
+! Licensed under the Apache License, Version 2.0 (the "License");
+! you may not use this file except in compliance with the License.
+! You may obtain a copy of the License at
+!
+!     http://www.apache.org/licenses/LICENSE-2.0
+!
+! Unless required by applicable law or agreed to in writing, software
+! distributed under the License is distributed on an "AS IS" BASIS,
+! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+! See the License for the specific language governing permissions and
+! limitations under the License.
+
+!> @file mctc/data/atomicrad.f90
+!> Provides atomic radii for all elements
+
+!> Atomic radii of the elements
+!>
+!> M. Mantina, R. Valero, C. J. Cramer, and D. G. Truhlar,
+!> in CRC Handbook of Chemistry and Physics, 91st Edition (2010-2011),
+!> edited by W. M. Haynes (CRC Press, Boca Raton, FL, 2010), pages 9-49-9-50;
+!> corrected Nov. 17, 2010 for the 92nd edition.
+module mctc_data_atomicrad
+   use mctc_env, only : wp
+   use mctc_io_convert, only : aatoau
+   use mctc_io_symbols, only : to_number
+   implicit none
+   private
+
+   public :: get_atomic_rad, atomic_rad
+
+
+   !> Get atomic radius for a species
+   interface get_atomic_rad
+      module procedure :: get_atomic_rad_symbol
+      module procedure :: get_atomic_rad_number
+   end interface get_atomic_rad
+
+
+   integer, parameter :: max_elem = 118
+
+   !> Atomic radii
+   real(wp), parameter :: atomic_rad(max_elem) = aatoau * [ &
+      & 0.32_wp, 0.37_wp, 1.30_wp, 0.99_wp, 0.84_wp, 0.75_wp, 0.71_wp, 0.64_wp, &
+      & 0.60_wp, 0.62_wp, 1.60_wp, 1.40_wp, 1.24_wp, 1.14_wp, 1.09_wp, 1.04_wp, &
+      & 1.00_wp, 1.01_wp, 2.00_wp, 1.74_wp, 1.59_wp, 1.48_wp, 1.44_wp, 1.30_wp, &
+      & 1.29_wp, 1.24_wp, 1.18_wp, 1.17_wp, 1.22_wp, 1.20_wp, 1.23_wp, 1.20_wp, &
+      & 1.20_wp, 1.18_wp, 1.17_wp, 1.16_wp, 2.15_wp, 1.90_wp, 1.76_wp, 1.64_wp, &
+      & 1.56_wp, 1.46_wp, 1.38_wp, 1.36_wp, 1.34_wp, 1.30_wp, 1.36_wp, 1.40_wp, &
+      & 1.42_wp, 1.40_wp, 1.40_wp, 1.37_wp, 1.36_wp, 1.36_wp, 2.38_wp, 2.06_wp, &
+      & 1.94_wp, 1.84_wp, 1.90_wp, 1.88_wp, 1.86_wp, 1.85_wp, 1.83_wp, 1.82_wp, &
+      & 1.81_wp, 1.80_wp, 1.79_wp, 1.77_wp, 1.77_wp, 1.78_wp, 1.74_wp, 1.64_wp, &
+      & 1.58_wp, 1.50_wp, 1.41_wp, 1.36_wp, 1.32_wp, 1.30_wp, 1.30_wp, 1.32_wp, &
+      & 1.44_wp, 1.45_wp, 1.50_wp, 1.42_wp, 1.48_wp, 1.46_wp, 2.42_wp, 2.11_wp, &
+      & 2.01_wp, 1.90_wp, 1.84_wp, 1.83_wp, 1.80_wp, 1.80_wp, 1.73_wp, 1.68_wp, &
+      & 1.68_wp, 1.68_wp, 1.65_wp, 1.67_wp, 1.73_wp, 1.76_wp, 1.61_wp, 1.57_wp, &
+      & 1.49_wp, 1.43_wp, 1.41_wp, 1.34_wp, 1.29_wp, 1.28_wp, 1.21_wp, 1.22_wp, &
+      & 1.36_wp, 1.43_wp, 1.62_wp, 1.75_wp, 1.65_wp, 1.57_wp]
+
+
+contains
+
+
+!> Get atomic radius for species with a given symbol
+elemental function get_atomic_rad_symbol(symbol) result(radius)
+
+   !> Element symbol
+   character(len=*), intent(in) :: symbol
+
+   !> atomic radius
+   real(wp) :: radius
+
+   radius = get_atomic_rad(to_number(symbol))
+
+end function get_atomic_rad_symbol
+
+
+!> Get atomic radius for species with a given atomic number
+elemental function get_atomic_rad_number(number) result(radius)
+
+   !> Atomic number
+   integer, intent(in) :: number
+
+   !> atomic radius
+   real(wp) :: radius
+
+   if (number > 0 .and. number <= size(atomic_rad, dim=1)) then
+      radius = atomic_rad(number)
+   else
+      radius = -1.0_wp
+   end if
+
+end function get_atomic_rad_number
+
+
+end module mctc_data_atomicrad