diff --git a/src/mctc/CMakeLists.txt b/src/mctc/CMakeLists.txt index 89bb8407..b60c58da 100644 --- a/src/mctc/CMakeLists.txt +++ b/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" ) diff --git a/src/mctc/cutoff.f90 b/src/mctc/cutoff.f90 new file mode 100644 index 00000000..df1fad2a --- /dev/null +++ b/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 diff --git a/src/mctc/data.f90 b/src/mctc/data.f90 new file mode 100644 index 00000000..ac671c92 --- /dev/null +++ b/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 diff --git a/src/mctc/data/CMakeLists.txt b/src/mctc/data/CMakeLists.txt new file mode 100644 index 00000000..4a268349 --- /dev/null +++ b/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) diff --git a/src/mctc/data/atomicrad.f90 b/src/mctc/data/atomicrad.f90 new file mode 100644 index 00000000..f72b1373 --- /dev/null +++ b/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 diff --git a/src/mctc/data/covrad.f90 b/src/mctc/data/covrad.f90 new file mode 100644 index 00000000..b3cb9cef --- /dev/null +++ b/src/mctc/data/covrad.f90 @@ -0,0 +1,105 @@ +! 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/covrad.f90 +!> Provides covalent radii for all elements + +!> Covalent radii (taken from Pyykko and Atsumi, Chem. Eur. J. 15, 2009, +!> 188-197), values for metals decreased by 10 % +module mctc_data_covrad + use mctc_env, only : wp + use mctc_io_convert, only : aatoau + use mctc_io_symbols, only : to_number + implicit none + private + + public :: get_covalent_rad + + + !> Covalent radii for coordination number + interface get_covalent_rad + module procedure :: get_covalent_rad_num + module procedure :: get_covalent_rad_sym + end interface get_covalent_rad + + + integer, parameter :: max_elem = 118 + + !> Covalent radii (taken from Pyykko and Atsumi, Chem. Eur. J. 15, 2009, + !> 188-197), values for metals decreased by 10 % + real(wp), parameter :: covalent_rad_2009(max_elem) = aatoau * [ & + & 0.32_wp,0.46_wp, & ! H,He + & 1.20_wp,0.94_wp,0.77_wp,0.75_wp,0.71_wp,0.63_wp,0.64_wp,0.67_wp, & ! Li-Ne + & 1.40_wp,1.25_wp,1.13_wp,1.04_wp,1.10_wp,1.02_wp,0.99_wp,0.96_wp, & ! Na-Ar + & 1.76_wp,1.54_wp, & ! K,Ca + & 1.33_wp,1.22_wp,1.21_wp,1.10_wp,1.07_wp, & ! Sc- + & 1.04_wp,1.00_wp,0.99_wp,1.01_wp,1.09_wp, & ! -Zn + & 1.12_wp,1.09_wp,1.15_wp,1.10_wp,1.14_wp,1.17_wp, & ! Ga-Kr + & 1.89_wp,1.67_wp, & ! Rb,Sr + & 1.47_wp,1.39_wp,1.32_wp,1.24_wp,1.15_wp, & ! Y- + & 1.13_wp,1.13_wp,1.08_wp,1.15_wp,1.23_wp, & ! -Cd + & 1.28_wp,1.26_wp,1.26_wp,1.23_wp,1.32_wp,1.31_wp, & ! In-Xe + & 2.09_wp,1.76_wp, & ! Cs,Ba + & 1.62_wp,1.47_wp,1.58_wp,1.57_wp,1.56_wp,1.55_wp,1.51_wp, & ! La-Eu + & 1.52_wp,1.51_wp,1.50_wp,1.49_wp,1.49_wp,1.48_wp,1.53_wp, & ! Gd-Yb + & 1.46_wp,1.37_wp,1.31_wp,1.23_wp,1.18_wp, & ! Lu- + & 1.16_wp,1.11_wp,1.12_wp,1.13_wp,1.32_wp, & ! -Hg + & 1.30_wp,1.30_wp,1.36_wp,1.31_wp,1.38_wp,1.42_wp, & ! Tl-Rn + & 2.01_wp,1.81_wp, & ! Fr,Ra + & 1.67_wp,1.58_wp,1.52_wp,1.53_wp,1.54_wp,1.55_wp,1.49_wp, & ! Ac-Am + & 1.49_wp,1.51_wp,1.51_wp,1.48_wp,1.50_wp,1.56_wp,1.58_wp, & ! Cm-No + & 1.45_wp,1.41_wp,1.34_wp,1.29_wp,1.27_wp, & ! Lr- + & 1.21_wp,1.16_wp,1.15_wp,1.09_wp,1.22_wp, & ! -Cn + & 1.36_wp,1.43_wp,1.46_wp,1.58_wp,1.48_wp,1.57_wp ] ! Nh-Og + + !> D3 covalent radii used to construct the coordination number + real(wp), parameter :: covalent_rad_d3(max_elem) = & + & 4.0_wp / 3.0_wp * covalent_rad_2009 + +contains + + +!> Get covalent radius for a given element symbol +elemental function get_covalent_rad_sym(sym) result(rad) + + !> Element symbol + character(len=*), intent(in) :: sym + + !> Covalent radius + real(wp) :: rad + + rad = get_covalent_rad(to_number(sym)) + +end function get_covalent_rad_sym + + +!> Get covalent radius for a given atomic number +elemental function get_covalent_rad_num(num) result(rad) + + !> Atomic number + integer, intent(in) :: num + + !> Covalent radius + real(wp) :: rad + + if (num > 0 .and. num <= size(covalent_rad_d3)) then + rad = covalent_rad_d3(num) + else + rad = 0.0_wp + end if + +end function get_covalent_rad_num + + +end module mctc_data_covrad diff --git a/src/mctc/data/meson.build b/src/mctc/data/meson.build new file mode 100644 index 00000000..7781eaf0 --- /dev/null +++ b/src/mctc/data/meson.build @@ -0,0 +1,20 @@ +# 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. + +srcs += files( + 'atomicrad.f90', + 'covrad.f90', + 'paulingen.f90', + 'vdwrad.f90' +) diff --git a/src/mctc/data/paulingen.f90 b/src/mctc/data/paulingen.f90 new file mode 100644 index 00000000..445c9d4a --- /dev/null +++ b/src/mctc/data/paulingen.f90 @@ -0,0 +1,102 @@ +! 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/paulingen.f90 +!> Provides electronegativities for all elements + +!> Pauling electronegativities +module mctc_data_paulingen + use mctc_env, only : wp + use mctc_io_symbols, only : to_number + implicit none + private + + public :: get_pauling_en, pauling_en + + + !> Get electronegativity for a species + interface get_pauling_en + module procedure :: get_pauling_en_symbol + module procedure :: get_pauling_en_number + end interface get_pauling_en + + + integer, parameter :: max_elem = 118 + + !> Pauling electronegativities, used for the covalent coordination number. + real(wp), parameter :: pauling_en(max_elem) = [ & + & 2.20_wp,3.00_wp, & ! H,He + & 0.98_wp,1.57_wp,2.04_wp,2.55_wp,3.04_wp,3.44_wp,3.98_wp,4.50_wp, & ! Li-Ne + & 0.93_wp,1.31_wp,1.61_wp,1.90_wp,2.19_wp,2.58_wp,3.16_wp,3.50_wp, & ! Na-Ar + & 0.82_wp,1.00_wp, & ! K,Ca + & 1.36_wp,1.54_wp,1.63_wp,1.66_wp,1.55_wp, & ! Sc- + & 1.83_wp,1.88_wp,1.91_wp,1.90_wp,1.65_wp, & ! -Zn + & 1.81_wp,2.01_wp,2.18_wp,2.55_wp,2.96_wp,3.00_wp, & ! Ga-Kr + & 0.82_wp,0.95_wp, & ! Rb,Sr + & 1.22_wp,1.33_wp,1.60_wp,2.16_wp,1.90_wp, & ! Y- + & 2.20_wp,2.28_wp,2.20_wp,1.93_wp,1.69_wp, & ! -Cd + & 1.78_wp,1.96_wp,2.05_wp,2.10_wp,2.66_wp,2.60_wp, & ! In-Xe + & 0.79_wp,0.89_wp, & ! Cs,Ba + & 1.10_wp,1.12_wp,1.13_wp,1.14_wp,1.15_wp,1.17_wp,1.18_wp, & ! La-Eu + & 1.20_wp,1.21_wp,1.22_wp,1.23_wp,1.24_wp,1.25_wp,1.26_wp, & ! Gd-Yb + & 1.27_wp,1.30_wp,1.50_wp,2.36_wp,1.90_wp, & ! Lu- + & 2.20_wp,2.20_wp,2.28_wp,2.54_wp,2.00_wp, & ! -Hg + & 1.62_wp,2.33_wp,2.02_wp,2.00_wp,2.20_wp,2.20_wp, & ! Tl-Rn + + & 0.79_wp,0.90_wp, & ! Fr,Ra + & 1.10_wp,1.30_wp,1.50_wp,1.38_wp,1.36_wp,1.28_wp,1.30_wp, & ! Ac-Am + & 1.30_wp,1.30_wp,1.30_wp,1.30_wp,1.30_wp,1.30_wp,1.30_wp, & ! Cm-No + & 1.30_wp, & ! Lr + ! only dummies below + & 1.50_wp,1.50_wp,1.50_wp,1.50_wp, & ! Rf-Bh + & 1.50_wp,1.50_wp,1.50_wp,1.50_wp,1.50_wp, & ! Hs-Cn + & 1.50_wp,1.50_wp,1.50_wp,1.50_wp,1.50_wp,1.50_wp ] ! Nh-Og + + +contains + + +!> Get electronegativity for species with a given symbol +elemental function get_pauling_en_symbol(symbol) result(en) + + !> Element symbol + character(len=*), intent(in) :: symbol + + !> atomic EN + real(wp) :: en + + en = get_pauling_en(to_number(symbol)) + +end function get_pauling_en_symbol + + +!> Get electronegativity for species with a given atomic number +elemental function get_pauling_en_number(number) result(en) + + !> Atomic number + integer, intent(in) :: number + + !> atomic EN + real(wp) :: en + + if (number > 0 .and. number <= size(pauling_en, dim=1)) then + en = pauling_en(number) + else + en = -1.0_wp + end if + +end function get_pauling_en_number + + +end module mctc_data_paulingen diff --git a/src/mctc/data/vdwrad.f90 b/src/mctc/data/vdwrad.f90 new file mode 100644 index 00000000..3f29d229 --- /dev/null +++ b/src/mctc/data/vdwrad.f90 @@ -0,0 +1,1213 @@ +! 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/vdwrad.f90 +!> Provides van der Waals radii for all elements and pairs of elements + +!> Pairwise van der Waals or cutoff radii based on the distance +!> where the first-order PBE0/def2-QZVP energy falls below a cutoff value +!> defined by: +!> S. Grimme, J. Antony, S. Ehrlich, H. Krieg, +!> J. Chem. Phys. 132, 154104 (2010), DOI: 10.1063/1.3382344 +!> extended by: +!> L. Trombach, S. Ehlert, S. Grimme, P Schwerdtfeger, J.-M. Mewes, +!> PCCP, 2019, 21,18048, DOI: 10.1039/c9cp02455g +!> L. Wittmann, I. Gordiy, M. Friede, B. Helmich-Paris, S. Grimme, A. Hansen, M. Bursch, +!> PCCP, 2024, 26, 21379, DOI: 10.1039/d4cp01514b +module mctc_data_vdwrad + use mctc_env, only : wp + use mctc_io_convert, only : aatoau + use mctc_io_symbols, only : to_number + use mctc_data_covrad, only : get_covalent_rad + implicit none + private + + public :: get_vdw_rad + + + interface get_vdw_rad + module procedure :: get_vdw_rad_sym + module procedure :: get_vdw_rad_pair_sym + module procedure :: get_vdw_rad_num + module procedure :: get_vdw_rad_pair_num + end interface get_vdw_rad + + + integer, parameter :: max_elem = 103 + + real(wp), parameter :: vdwrad(max_elem*(1+max_elem)/2) = aatoau * [ & + 2.1823_wp, 1.8547_wp, 1.7347_wp, 2.9086_wp, 2.5732_wp, & + 3.4956_wp, 2.3550_wp, 2.5095_wp, 2.9802_wp, 3.0982_wp, & + 2.5141_wp, 2.3917_wp, 2.9977_wp, 2.9484_wp, 3.2160_wp, & + 2.4492_wp, 2.2527_wp, 3.1933_wp, 3.0214_wp, 2.9531_wp, & + 2.9103_wp, 2.3667_wp, 2.1328_wp, 2.8784_wp, 2.7660_wp, & + 2.7776_wp, 2.7063_wp, 2.6225_wp, 2.1768_wp, 2.0625_wp, & + 2.6395_wp, 2.6648_wp, 2.6482_wp, 2.5697_wp, 2.4846_wp, & + 2.4817_wp, 2.0646_wp, 1.9891_wp, 2.5086_wp, 2.6908_wp, & + 2.6233_wp, 2.4770_wp, 2.3885_wp, 2.3511_wp, 2.2996_wp, & + 1.9892_wp, 1.9251_wp, 2.4190_wp, 2.5473_wp, 2.4994_wp, & + 2.4091_wp, 2.3176_wp, 2.2571_wp, 2.1946_wp, 2.1374_wp, & + 2.9898_wp, 2.6397_wp, 3.6031_wp, 3.1219_wp, 3.7620_wp, & + 3.2485_wp, 2.9357_wp, 2.7093_wp, 2.5781_wp, 2.4839_wp, & + 3.7082_wp, 2.5129_wp, 2.7321_wp, 3.1052_wp, 3.2962_wp, & + 3.1331_wp, 3.2000_wp, 2.9586_wp, 3.0822_wp, 2.8582_wp, & + 2.7120_wp, 3.2570_wp, 3.4839_wp, 2.8766_wp, 2.7427_wp, & + 3.2776_wp, 3.2363_wp, 3.5929_wp, 3.2826_wp, 3.0911_wp, & + 2.9369_wp, 2.9030_wp, 2.7789_wp, 3.3921_wp, 3.3970_wp, & + 4.0106_wp, 2.8884_wp, 2.6605_wp, 3.7513_wp, 3.1613_wp, & + 3.3605_wp, 3.3325_wp, 3.0991_wp, 2.9297_wp, 2.8674_wp, & + 2.7571_wp, 3.8129_wp, 3.3266_wp, 3.7105_wp, 3.7917_wp, & + 2.8304_wp, 2.5538_wp, 3.3932_wp, 3.1193_wp, 3.1866_wp, & + 3.1245_wp, 3.0465_wp, 2.8727_wp, 2.7664_wp, 2.6926_wp, & + 3.4608_wp, 3.2984_wp, 3.5142_wp, 3.5418_wp, 3.5017_wp, & + 2.6190_wp, 2.4797_wp, 3.1331_wp, 3.0540_wp, 3.0651_wp, & + 2.9879_wp, 2.9054_wp, 2.8805_wp, 2.7330_wp, 2.6331_wp, & + 3.2096_wp, 3.5668_wp, 3.3684_wp, 3.3686_wp, 3.3180_wp, & + 3.3107_wp, 2.4757_wp, 2.4019_wp, 2.9789_wp, 3.1468_wp, & + 2.9768_wp, 2.8848_wp, 2.7952_wp, 2.7457_wp, 2.6881_wp, & + 2.5728_wp, 3.0574_wp, 3.3264_wp, 3.3562_wp, 3.2529_wp, & + 3.1916_wp, 3.1523_wp, 3.1046_wp, 2.3725_wp, 2.3289_wp, & + 2.8760_wp, 2.9804_wp, 2.9093_wp, 2.8040_wp, 2.7071_wp, & + 2.6386_wp, 2.5720_wp, 2.5139_wp, 2.9517_wp, 3.1606_wp, & + 3.2085_wp, 3.1692_wp, 3.0982_wp, 3.0352_wp, 2.9730_wp, & + 2.9148_wp, 3.2147_wp, 2.8315_wp, 3.8724_wp, 3.4621_wp, & + 3.8823_wp, 3.3760_wp, 3.0746_wp, 2.8817_wp, 2.7552_wp, & + 2.6605_wp, 3.9740_wp, 3.6192_wp, 3.6569_wp, 3.9586_wp, & + 3.6188_wp, 3.3917_wp, 3.2479_wp, 3.1434_wp, 4.2411_wp, & + 2.7597_wp, 3.0588_wp, 3.3474_wp, 3.6214_wp, 3.4353_wp, & + 3.4729_wp, 3.2487_wp, 3.3200_wp, 3.0914_wp, 2.9403_wp, & + 3.4972_wp, 3.7993_wp, 3.6773_wp, 3.8678_wp, 3.5808_wp, & + 3.8243_wp, 3.5826_wp, 3.4156_wp, 3.8765_wp, 4.1035_wp, & + 2.7361_wp, 2.9765_wp, 3.2475_wp, 3.5004_wp, 3.4185_wp, & + 3.4378_wp, 3.2084_wp, 3.2787_wp, 3.0604_wp, 2.9187_wp, & + 3.4037_wp, 3.6759_wp, 3.6586_wp, 3.8327_wp, 3.5372_wp, & + 3.7665_wp, 3.5310_wp, 3.3700_wp, 3.7788_wp, 3.9804_wp, & + 3.8903_wp, 2.6832_wp, 2.9060_wp, 3.2613_wp, 3.4359_wp, & + 3.3538_wp, 3.3860_wp, 3.1550_wp, 3.2300_wp, 3.0133_wp, & + 2.8736_wp, 3.4024_wp, 3.6142_wp, 3.5979_wp, 3.5295_wp, & + 3.4834_wp, 3.7140_wp, 3.4782_wp, 3.3170_wp, 3.7434_wp, & + 3.9623_wp, 3.8181_wp, 3.7642_wp, 2.6379_wp, 2.8494_wp, & + 3.1840_wp, 3.4225_wp, 3.2771_wp, 3.3401_wp, 3.1072_wp, & + 3.1885_wp, 2.9714_wp, 2.8319_wp, 3.3315_wp, 3.5979_wp, & + 3.5256_wp, 3.4980_wp, 3.4376_wp, 3.6714_wp, 3.4346_wp, & + 3.2723_wp, 3.6859_wp, 3.8985_wp, 3.7918_wp, 3.7372_wp, & + 3.7211_wp, 2.9230_wp, 2.6223_wp, 3.4161_wp, 2.8999_wp, & + 3.0557_wp, 3.3308_wp, 3.0555_wp, 2.8508_wp, 2.7385_wp, & + 2.6640_wp, 3.5263_wp, 3.0277_wp, 3.2990_wp, 3.7721_wp, & + 3.5017_wp, 3.2751_wp, 3.1368_wp, 3.0435_wp, 3.7873_wp, & + 3.2858_wp, 3.2140_wp, 3.1727_wp, 3.2178_wp, 3.4414_wp, & + 2.5490_wp, 2.7623_wp, 3.0991_wp, 3.3252_wp, 3.1836_wp, & + 3.2428_wp, 3.0259_wp, 3.1225_wp, 2.9032_wp, 2.7621_wp, & + 3.2490_wp, 3.5110_wp, 3.4429_wp, 3.3845_wp, 3.3574_wp, & + 3.6045_wp, 3.3658_wp, 3.2013_wp, 3.6110_wp, 3.8241_wp, & + 3.7090_wp, 3.6496_wp, 3.6333_wp, 3.0896_wp, 3.5462_wp, & + 2.4926_wp, 2.7136_wp, 3.0693_wp, 3.2699_wp, 3.1272_wp, & + 3.1893_wp, 2.9658_wp, 3.0972_wp, 2.8778_wp, 2.7358_wp, & + 3.2206_wp, 3.4566_wp, 3.3896_wp, 3.3257_wp, 3.2946_wp, & + 3.5693_wp, 3.3312_wp, 3.1670_wp, 3.5805_wp, 3.7711_wp, & + 3.6536_wp, 3.5927_wp, 3.5775_wp, 3.0411_wp, 3.4885_wp, & + 3.4421_wp, 2.4667_wp, 2.6709_wp, 3.0575_wp, 3.2357_wp, & + 3.0908_wp, 3.1537_wp, 2.9235_wp, 3.0669_wp, 2.8476_wp, & + 2.7054_wp, 3.2064_wp, 3.4519_wp, 3.3593_wp, 3.2921_wp, & + 3.2577_wp, 3.2161_wp, 3.2982_wp, 3.1339_wp, 3.5606_wp, & + 3.7582_wp, 3.6432_wp, 3.5833_wp, 3.5691_wp, 3.0161_wp, & + 3.4812_wp, 3.4339_wp, 3.4327_wp, 2.4515_wp, 2.6338_wp, & + 3.0511_wp, 3.2229_wp, 3.0630_wp, 3.1265_wp, 2.8909_wp, & + 3.0253_wp, 2.8184_wp, 2.6764_wp, 3.1968_wp, 3.4114_wp, & + 3.3492_wp, 3.2691_wp, 3.2320_wp, 3.1786_wp, 3.2680_wp, & + 3.1036_wp, 3.5453_wp, 3.7259_wp, 3.6090_wp, 3.5473_wp, & + 3.5327_wp, 3.0018_wp, 3.4413_wp, 3.3907_wp, 3.3593_wp, & + 3.3462_wp, 2.4413_wp, 2.6006_wp, 3.0540_wp, 3.1987_wp, & + 3.0490_wp, 3.1058_wp, 2.8643_wp, 2.9948_wp, 2.7908_wp, & + 2.6491_wp, 3.1950_wp, 3.3922_wp, 3.3316_wp, 3.2585_wp, & + 3.2136_wp, 3.1516_wp, 3.2364_wp, 3.0752_wp, 3.5368_wp, & + 3.7117_wp, 3.5941_wp, 3.5313_wp, 3.5164_wp, 2.9962_wp, & + 3.4225_wp, 3.3699_wp, 3.3370_wp, 3.3234_wp, 3.3008_wp, & + 2.4318_wp, 2.5729_wp, 3.0416_wp, 3.1639_wp, 3.0196_wp, & + 3.0843_wp, 2.8413_wp, 2.7436_wp, 2.7608_wp, 2.6271_wp, & + 3.1811_wp, 3.3591_wp, 3.3045_wp, 3.2349_wp, 3.1942_wp, & + 3.1291_wp, 3.2111_wp, 3.0534_wp, 3.5189_wp, 3.6809_wp, & + 3.5635_wp, 3.5001_wp, 3.4854_wp, 2.9857_wp, 3.3897_wp, & + 3.3363_wp, 3.3027_wp, 3.2890_wp, 3.2655_wp, 3.2309_wp, & + 2.8502_wp, 2.6934_wp, 3.2467_wp, 3.1921_wp, 3.5663_wp, & + 3.2541_wp, 3.0571_wp, 2.9048_wp, 2.8657_wp, 2.7438_wp, & + 3.3547_wp, 3.3510_wp, 3.9837_wp, 3.6871_wp, 3.4862_wp, & + 3.3389_wp, 3.2413_wp, 3.1708_wp, 3.6096_wp, 3.6280_wp, & + 3.6860_wp, 3.5568_wp, 3.4836_wp, 3.2868_wp, 3.3994_wp, & + 3.3476_wp, 3.3170_wp, 3.2950_wp, 3.2874_wp, 3.2606_wp, & + 3.9579_wp, 2.9226_wp, 2.6838_wp, 3.7867_wp, 3.1732_wp, & + 3.3872_wp, 3.3643_wp, 3.1267_wp, 2.9541_wp, 2.8505_wp, & + 2.7781_wp, 3.8475_wp, 3.3336_wp, 3.7359_wp, 3.8266_wp, & + 3.5733_wp, 3.3959_wp, 3.2775_wp, 3.1915_wp, 3.9878_wp, & + 3.8816_wp, 3.5810_wp, 3.5364_wp, 3.5060_wp, 3.8097_wp, & + 3.3925_wp, 3.3348_wp, 3.3019_wp, 3.2796_wp, 3.2662_wp, & + 3.2464_wp, 3.7136_wp, 3.8619_wp, 2.9140_wp, 2.6271_wp, & + 3.4771_wp, 3.1774_wp, 3.2560_wp, 3.1970_wp, 3.1207_wp, & + 2.9406_wp, 2.8322_wp, 2.7571_wp, 3.5455_wp, 3.3514_wp, & + 3.5837_wp, 3.6177_wp, 3.5816_wp, 3.3902_wp, 3.2604_wp, & + 3.1652_wp, 3.7037_wp, 3.6283_wp, 3.5858_wp, 3.5330_wp, & + 3.4884_wp, 3.5789_wp, 3.4094_wp, 3.3473_wp, 3.3118_wp, & + 3.2876_wp, 3.2707_wp, 3.2521_wp, 3.5570_wp, 3.6496_wp, & + 3.6625_wp, 2.7300_wp, 2.5870_wp, 3.2471_wp, 3.1487_wp, & + 3.1667_wp, 3.0914_wp, 3.0107_wp, 2.9812_wp, 2.8300_wp, & + 2.7284_wp, 3.3259_wp, 3.3182_wp, 3.4707_wp, 3.4748_wp, & + 3.4279_wp, 3.4182_wp, 3.2547_wp, 3.1353_wp, 3.5116_wp, & + 3.9432_wp, 3.8828_wp, 3.8303_wp, 3.7880_wp, 3.3760_wp, & + 3.7218_wp, 3.3408_wp, 3.3059_wp, 3.2698_wp, 3.2446_wp, & + 3.2229_wp, 3.4422_wp, 3.5023_wp, 3.5009_wp, 3.5268_wp, & + 2.6026_wp, 2.5355_wp, 3.1129_wp, 3.2863_wp, 3.1029_wp, & + 3.0108_wp, 2.9227_wp, 2.8694_wp, 2.8109_wp, 2.6929_wp, & + 3.1958_wp, 3.4670_wp, 3.4018_wp, 3.3805_wp, 3.3218_wp, & + 3.2815_wp, 3.2346_wp, 3.0994_wp, 3.3937_wp, 3.7266_wp, & + 3.6697_wp, 3.6164_wp, 3.5730_wp, 3.2522_wp, 3.5051_wp, & + 3.4686_wp, 3.4355_wp, 3.4084_wp, 3.3748_wp, 3.3496_wp, & + 3.3692_wp, 3.4052_wp, 3.3910_wp, 3.3849_wp, 3.3662_wp, & + 2.5087_wp, 2.4814_wp, 3.0239_wp, 3.1312_wp, 3.0535_wp, & + 2.9457_wp, 2.8496_wp, 2.7780_wp, 2.7828_wp, 2.6532_wp, & + 3.1063_wp, 3.3143_wp, 3.3549_wp, 3.3120_wp, 3.2421_wp, & + 3.1787_wp, 3.1176_wp, 3.0613_wp, 3.3082_wp, 3.5755_wp, & + 3.5222_wp, 3.4678_wp, 3.4231_wp, 3.1684_wp, 3.3528_wp, & + 3.3162_wp, 3.2827_wp, 3.2527_wp, 3.2308_wp, 3.2029_wp, & + 3.3173_wp, 3.3343_wp, 3.3092_wp, 3.2795_wp, 3.2452_wp, & + 3.2096_wp, 3.2893_wp, 2.8991_wp, 4.0388_wp, 3.6100_wp, & + 3.9388_wp, 3.4475_wp, 3.1590_wp, 2.9812_wp, 2.8586_wp, & + 2.7683_wp, 4.1428_wp, 3.7911_wp, 3.8225_wp, 4.0372_wp, & + 3.7059_wp, 3.4935_wp, 3.3529_wp, 3.2492_wp, 4.4352_wp, & + 4.0826_wp, 3.9733_wp, 3.9254_wp, 3.8646_wp, 3.9315_wp, & + 3.7837_wp, 3.7465_wp, 3.7211_wp, 3.7012_wp, 3.6893_wp, & + 3.6676_wp, 3.7736_wp, 4.0660_wp, 3.7926_wp, 3.6158_wp, & + 3.5017_wp, 3.4166_wp, 4.6176_wp, 2.8786_wp, 3.1658_wp, & + 3.5823_wp, 3.7689_wp, 3.5762_wp, 3.5789_wp, 3.3552_wp, & + 3.4004_wp, 3.1722_wp, 3.0212_wp, 3.7241_wp, 3.9604_wp, & + 3.8500_wp, 3.9844_wp, 3.7035_wp, 3.9161_wp, 3.6751_wp, & + 3.5075_wp, 4.1151_wp, 4.2877_wp, 4.1579_wp, 4.1247_wp, & + 4.0617_wp, 3.4874_wp, 3.9848_wp, 3.9280_wp, 3.9079_wp, & + 3.8751_wp, 3.8604_wp, 3.8277_wp, 3.8002_wp, 3.9981_wp, & + 3.7544_wp, 4.0371_wp, 3.8225_wp, 3.6718_wp, 4.3092_wp, & + 4.4764_wp, 2.8997_wp, 3.0953_wp, 3.4524_wp, 3.6107_wp, & + 3.6062_wp, 3.5783_wp, 3.3463_wp, 3.3855_wp, 3.1746_wp, & + 3.0381_wp, 3.6019_wp, 3.7938_wp, 3.8697_wp, 3.9781_wp, & + 3.6877_wp, 3.8736_wp, 3.6451_wp, 3.4890_wp, 3.9858_wp, & + 4.1179_wp, 4.0430_wp, 3.9563_wp, 3.9182_wp, 3.4002_wp, & + 3.8310_wp, 3.7716_wp, 3.7543_wp, 3.7203_wp, 3.7053_wp, & + 3.6742_wp, 3.8318_wp, 3.7631_wp, 3.7392_wp, 3.9892_wp, & + 3.7832_wp, 3.6406_wp, 4.1701_wp, 4.3016_wp, 4.2196_wp, & + 2.8535_wp, 3.0167_wp, 3.3978_wp, 3.5363_wp, 3.5393_wp, & + 3.5301_wp, 3.2960_wp, 3.3352_wp, 3.1287_wp, 2.9967_wp, & + 3.6659_wp, 3.7239_wp, 3.8070_wp, 3.7165_wp, 3.6368_wp, & + 3.8162_wp, 3.5885_wp, 3.4336_wp, 3.9829_wp, 4.0529_wp, & + 3.9584_wp, 3.9025_wp, 3.8607_wp, 3.3673_wp, 3.7658_wp, & + 3.7035_wp, 3.6866_wp, 3.6504_wp, 3.6339_wp, 3.6024_wp, & + 3.7708_wp, 3.7283_wp, 3.6896_wp, 3.9315_wp, 3.7250_wp, & + 3.5819_wp, 4.1457_wp, 4.2280_wp, 4.1130_wp, 4.0597_wp, & + 3.0905_wp, 2.7998_wp, 3.6448_wp, 3.0739_wp, 3.2996_wp, & + 3.5262_wp, 3.2559_wp, 3.0518_wp, 2.9394_wp, 2.8658_wp, & + 3.7514_wp, 3.2295_wp, 3.5643_wp, 3.7808_wp, 3.6931_wp, & + 3.4723_wp, 3.3357_wp, 3.2429_wp, 4.0280_wp, 3.5589_wp, & + 3.4636_wp, 3.4994_wp, 3.4309_wp, 3.6177_wp, 3.2946_wp, & + 3.2376_wp, 3.2050_wp, 3.1847_wp, 3.1715_wp, 3.1599_wp, & + 3.5555_wp, 3.8111_wp, 3.7693_wp, 3.5718_wp, 3.4498_wp, & + 3.3662_wp, 4.1608_wp, 3.7417_wp, 3.6536_wp, 3.6154_wp, & + 3.8596_wp, 3.0301_wp, 2.7312_wp, 3.5821_wp, 3.0473_wp, & + 3.2137_wp, 3.4679_wp, 3.1975_wp, 2.9969_wp, 2.8847_wp, & + 2.8110_wp, 3.6931_wp, 3.2076_wp, 3.4943_wp, 3.5956_wp, & + 3.6379_wp, 3.4190_wp, 3.2808_wp, 3.1860_wp, 3.9850_wp, & + 3.5105_wp, 3.4330_wp, 3.3797_wp, 3.4155_wp, 3.6033_wp, & + 3.2737_wp, 3.2145_wp, 3.1807_wp, 3.1596_wp, 3.1461_wp, & + 3.1337_wp, 3.4812_wp, 3.6251_wp, 3.7152_wp, 3.5201_wp, & + 3.3966_wp, 3.3107_wp, 4.1128_wp, 3.6899_wp, 3.6082_wp, & + 3.5604_wp, 3.7834_wp, 3.7543_wp, 2.9189_wp, 2.6777_wp, & + 3.4925_wp, 2.9648_wp, 3.1216_wp, 3.2940_wp, 3.0975_wp, & + 2.9757_wp, 2.8493_wp, 2.7638_wp, 3.6085_wp, 3.1214_wp, & + 3.4006_wp, 3.4793_wp, 3.5147_wp, 3.3806_wp, 3.2356_wp, & + 3.1335_wp, 3.9144_wp, 3.4183_wp, 3.3369_wp, 3.2803_wp, & + 3.2679_wp, 3.4871_wp, 3.1714_wp, 3.1521_wp, 3.1101_wp, & + 3.0843_wp, 3.0670_wp, 3.0539_wp, 3.3890_wp, 3.5086_wp, & + 3.5895_wp, 3.4783_wp, 3.3484_wp, 3.2559_wp, 4.0422_wp, & + 3.5967_wp, 3.5113_wp, 3.4576_wp, 3.6594_wp, 3.6313_wp, & + 3.5690_wp, 2.8578_wp, 2.6334_wp, 3.4673_wp, 2.9245_wp, & + 3.0732_wp, 3.2435_wp, 3.0338_wp, 2.9462_wp, 2.8143_wp, & + 2.7240_wp, 3.5832_wp, 3.0789_wp, 3.3617_wp, 3.4246_wp, & + 3.4505_wp, 3.3443_wp, 3.1964_wp, 3.0913_wp, 3.8921_wp, & + 3.3713_wp, 3.2873_wp, 3.2281_wp, 3.2165_wp, 3.4386_wp, & + 3.1164_wp, 3.1220_wp, 3.0761_wp, 3.0480_wp, 3.0295_wp, & + 3.0155_wp, 3.3495_wp, 3.4543_wp, 3.5260_wp, 3.4413_wp, & + 3.3085_wp, 3.2134_wp, 4.0170_wp, 3.5464_wp, 3.4587_wp, & + 3.4006_wp, 3.6027_wp, 3.5730_wp, 3.4945_wp, 3.4623_wp, & + 2.8240_wp, 2.5960_wp, 3.4635_wp, 2.9032_wp, 3.0431_wp, & + 3.2115_wp, 2.9892_wp, 2.9148_wp, 2.7801_wp, 2.6873_wp, & + 3.5776_wp, 3.0568_wp, 3.3433_wp, 3.3949_wp, 3.4132_wp, & + 3.3116_wp, 3.1616_wp, 3.0548_wp, 3.8859_wp, 3.3719_wp, & + 3.2917_wp, 3.2345_wp, 3.2274_wp, 3.4171_wp, 3.1293_wp, & + 3.0567_wp, 3.0565_wp, 3.0274_wp, 3.0087_wp, 2.9939_wp, & + 3.3293_wp, 3.4249_wp, 3.4902_wp, 3.4091_wp, 3.2744_wp, & + 3.1776_wp, 4.0078_wp, 3.5374_wp, 3.4537_wp, 3.3956_wp, & + 3.5747_wp, 3.5430_wp, 3.4522_wp, 3.4160_wp, 3.3975_wp, & + 2.8004_wp, 2.5621_wp, 3.4617_wp, 2.9154_wp, 3.0203_wp, & + 3.1875_wp, 2.9548_wp, 2.8038_wp, 2.7472_wp, 2.6530_wp, & + 3.5736_wp, 3.0584_wp, 3.3304_wp, 3.3748_wp, 3.3871_wp, & + 3.2028_wp, 3.1296_wp, 3.0214_wp, 3.8796_wp, 3.3337_wp, & + 3.2492_wp, 3.1883_wp, 3.1802_wp, 3.4050_wp, 3.0756_wp, & + 3.0478_wp, 3.0322_wp, 3.0323_wp, 3.0163_wp, 3.0019_wp, & + 3.3145_wp, 3.4050_wp, 3.4656_wp, 3.3021_wp, 3.2433_wp, & + 3.1453_wp, 3.9991_wp, 3.5017_wp, 3.4141_wp, 3.3520_wp, & + 3.5583_wp, 3.5251_wp, 3.4243_wp, 3.3851_wp, 3.3662_wp, & + 3.3525_wp, 2.7846_wp, 2.5324_wp, 3.4652_wp, 2.8759_wp, & + 3.0051_wp, 3.1692_wp, 2.9273_wp, 2.7615_wp, 2.7164_wp, & + 2.6212_wp, 3.5744_wp, 3.0275_wp, 3.3249_wp, 3.3627_wp, & + 3.3686_wp, 3.1669_wp, 3.0584_wp, 2.9915_wp, 3.8773_wp, & + 3.3099_wp, 3.2231_wp, 3.1600_wp, 3.1520_wp, 3.4023_wp, & + 3.0426_wp, 3.0099_wp, 2.9920_wp, 2.9809_wp, 2.9800_wp, & + 2.9646_wp, 3.3068_wp, 3.3930_wp, 3.4486_wp, 3.2682_wp, & + 3.1729_wp, 3.1168_wp, 3.9952_wp, 3.4796_wp, 3.3901_wp, & + 3.3255_wp, 3.5530_wp, 3.5183_wp, 3.4097_wp, 3.3683_wp, & + 3.3492_wp, 3.3360_wp, 3.3308_wp, 2.5424_wp, 2.6601_wp, & + 3.2555_wp, 3.2807_wp, 3.1384_wp, 3.1737_wp, 2.9397_wp, & + 2.8429_wp, 2.8492_wp, 2.7225_wp, 3.3875_wp, 3.4910_wp, & + 3.4520_wp, 3.3608_wp, 3.3036_wp, 3.2345_wp, 3.2999_wp, & + 3.1487_wp, 3.7409_wp, 3.8392_wp, 3.7148_wp, 3.6439_wp, & + 3.6182_wp, 3.1753_wp, 3.5210_wp, 3.4639_wp, 3.4265_wp, & + 3.4075_wp, 3.3828_wp, 3.3474_wp, 3.4071_wp, 3.3754_wp, & + 3.3646_wp, 3.3308_wp, 3.4393_wp, 3.2993_wp, 3.8768_wp, & + 3.9891_wp, 3.8310_wp, 3.7483_wp, 3.3417_wp, 3.3019_wp, & + 3.2250_wp, 3.1832_wp, 3.1578_wp, 3.1564_wp, 3.1224_wp, & + 3.4620_wp, 2.9743_wp, 2.8058_wp, 3.4830_wp, 3.3474_wp, & + 3.6863_wp, 3.3617_wp, 3.1608_wp, 3.0069_wp, 2.9640_wp, & + 2.8427_wp, 3.5885_wp, 3.5219_wp, 4.1314_wp, 3.8120_wp, & + 3.6015_wp, 3.4502_wp, 3.3498_wp, 3.2777_wp, 3.8635_wp, & + 3.8232_wp, 3.8486_wp, 3.7215_wp, 3.6487_wp, 3.4724_wp, & + 3.5627_wp, 3.5087_wp, 3.4757_wp, 3.4517_wp, 3.4423_wp, & + 3.4139_wp, 4.1028_wp, 3.8388_wp, 3.6745_wp, 3.5562_wp, & + 3.4806_wp, 3.4272_wp, 4.0182_wp, 3.9991_wp, 4.0007_wp, & + 3.9282_wp, 3.7238_wp, 3.6498_wp, 3.5605_wp, 3.5211_wp, & + 3.5009_wp, 3.4859_wp, 3.4785_wp, 3.5621_wp, 4.2623_wp, & + 3.0775_wp, 2.8275_wp, 4.0181_wp, 3.3385_wp, 3.5379_wp, & + 3.5036_wp, 3.2589_wp, 3.0804_wp, 3.0094_wp, 2.9003_wp, & + 4.0869_wp, 3.5088_wp, 3.9105_wp, 3.9833_wp, 3.7176_wp, & + 3.5323_wp, 3.4102_wp, 3.3227_wp, 4.2702_wp, 4.0888_wp, & + 3.7560_wp, 3.7687_wp, 3.6681_wp, 3.6405_wp, 3.5569_wp, & + 3.4990_wp, 3.4659_wp, 3.4433_wp, 3.4330_wp, 3.4092_wp, & + 3.8867_wp, 4.0190_wp, 3.7961_wp, 3.6412_wp, 3.5405_wp, & + 3.4681_wp, 4.3538_wp, 4.2136_wp, 3.9381_wp, 3.8912_wp, & + 3.9681_wp, 3.7909_wp, 3.6774_wp, 3.6262_wp, 3.5999_wp, & + 3.5823_wp, 3.5727_wp, 3.5419_wp, 4.0245_wp, 4.1874_wp, & + 3.0893_wp, 2.7917_wp, 3.7262_wp, 3.3518_wp, 3.4241_wp, & + 3.5433_wp, 3.2773_wp, 3.0890_wp, 2.9775_wp, 2.9010_wp, & + 3.8048_wp, 3.5362_wp, 3.7746_wp, 3.7911_wp, 3.7511_wp, & + 3.5495_wp, 3.4149_wp, 3.3177_wp, 4.0129_wp, 3.8370_wp, & + 3.7739_wp, 3.7125_wp, 3.7152_wp, 3.7701_wp, 3.5813_wp, & + 3.5187_wp, 3.4835_wp, 3.4595_wp, 3.4439_wp, 3.4242_wp, & + 3.7476_wp, 3.8239_wp, 3.8346_wp, 3.6627_wp, 3.5479_wp, & + 3.4639_wp, 4.1026_wp, 3.9733_wp, 3.9292_wp, 3.8667_wp, & + 3.9513_wp, 3.8959_wp, 3.7698_wp, 3.7089_wp, 3.6765_wp, & + 3.6548_wp, 3.6409_wp, 3.5398_wp, 3.8759_wp, 3.9804_wp, & + 4.0150_wp, 2.9091_wp, 2.7638_wp, 3.5066_wp, 3.3377_wp, & + 3.3481_wp, 3.2633_wp, 3.1810_wp, 3.1428_wp, 2.9872_wp, & + 2.8837_wp, 3.5929_wp, 3.5183_wp, 3.6729_wp, 3.6596_wp, & + 3.6082_wp, 3.5927_wp, 3.4224_wp, 3.2997_wp, 3.8190_wp, & + 4.1865_wp, 4.1114_wp, 4.0540_wp, 3.6325_wp, 3.5697_wp, & + 3.5561_wp, 3.5259_wp, 3.4901_wp, 3.4552_wp, 3.4315_wp, & + 3.4091_wp, 3.6438_wp, 3.6879_wp, 3.6832_wp, 3.7043_wp, & + 3.5557_wp, 3.4466_wp, 3.9203_wp, 4.2919_wp, 4.2196_wp, & + 4.1542_wp, 3.7573_wp, 3.7039_wp, 3.6546_wp, 3.6151_wp, & + 3.5293_wp, 3.4849_wp, 3.4552_wp, 3.5192_wp, 3.7673_wp, & + 3.8359_wp, 3.8525_wp, 3.8901_wp, 2.7806_wp, 2.7209_wp, & + 3.3812_wp, 3.4958_wp, 3.2913_wp, 3.1888_wp, 3.0990_wp, & + 3.0394_wp, 2.9789_wp, 2.8582_wp, 3.4716_wp, 3.6883_wp, & + 3.6105_wp, 3.5704_wp, 3.5059_wp, 3.4619_wp, 3.4138_wp, & + 3.2742_wp, 3.7080_wp, 3.9773_wp, 3.9010_wp, 3.8409_wp, & + 3.7944_wp, 3.4465_wp, 3.7235_wp, 3.6808_wp, 3.6453_wp, & + 3.6168_wp, 3.5844_wp, 3.5576_wp, 3.5772_wp, 3.5959_wp, & + 3.5768_wp, 3.5678_wp, 3.5486_wp, 3.4228_wp, 3.8107_wp, & + 4.0866_wp, 4.0169_wp, 3.9476_wp, 3.6358_wp, 3.5800_wp, & + 3.5260_wp, 3.4838_wp, 3.4501_wp, 3.4204_wp, 3.3553_wp, & + 3.6487_wp, 3.6973_wp, 3.7398_wp, 3.7405_wp, 3.7459_wp, & + 3.7380_wp, 2.6848_wp, 2.6740_wp, 3.2925_wp, 3.3386_wp, & + 3.2473_wp, 3.1284_wp, 3.0301_wp, 2.9531_wp, 2.9602_wp, & + 2.8272_wp, 3.3830_wp, 3.5358_wp, 3.5672_wp, 3.5049_wp, & + 3.4284_wp, 3.3621_wp, 3.3001_wp, 3.2451_wp, 3.6209_wp, & + 3.8299_wp, 3.7543_wp, 3.6920_wp, 3.6436_wp, 3.3598_wp, & + 3.5701_wp, 3.5266_wp, 3.4904_wp, 3.4590_wp, 3.4364_wp, & + 3.4077_wp, 3.5287_wp, 3.5280_wp, 3.4969_wp, 3.4650_wp, & + 3.4304_wp, 3.3963_wp, 3.7229_wp, 3.9402_wp, 3.8753_wp, & + 3.8035_wp, 3.5499_wp, 3.4913_wp, 3.4319_wp, 3.3873_wp, & + 3.3520_wp, 3.3209_wp, 3.2948_wp, 3.5052_wp, 3.6465_wp, & + 3.6696_wp, 3.6577_wp, 3.6388_wp, 3.6142_wp, 3.5889_wp, & + 3.3968_wp, 3.0122_wp, 4.2241_wp, 3.7887_wp, 4.0049_wp, & + 3.5384_wp, 3.2698_wp, 3.1083_wp, 2.9917_wp, 2.9057_wp, & + 4.3340_wp, 3.9900_wp, 4.6588_wp, 4.1278_wp, 3.8125_wp, & + 3.6189_wp, 3.4851_wp, 3.3859_wp, 4.6531_wp, 4.3134_wp, & + 4.2258_wp, 4.1309_wp, 4.0692_wp, 4.0944_wp, 3.9850_wp, & + 3.9416_wp, 3.9112_wp, 3.8873_wp, 3.8736_wp, 3.8473_wp, & + 4.6027_wp, 4.1538_wp, 3.8994_wp, 3.7419_wp, 3.6356_wp, & + 3.5548_wp, 4.8353_wp, 4.5413_wp, 4.3891_wp, 4.3416_wp, & + 4.3243_wp, 4.2753_wp, 4.2053_wp, 4.1790_wp, 4.1685_wp, & + 4.1585_wp, 4.1536_wp, 4.0579_wp, 4.1980_wp, 4.4564_wp, & + 4.2192_wp, 4.0528_wp, 3.9489_wp, 3.8642_wp, 5.0567_wp, & + 3.0630_wp, 3.3271_wp, 4.0432_wp, 4.0046_wp, 4.1555_wp, & + 3.7426_wp, 3.5130_wp, 3.5174_wp, 3.2884_wp, 3.1378_wp, & + 4.1894_wp, 4.2321_wp, 4.1725_wp, 4.1833_wp, 3.8929_wp, & + 4.0544_wp, 3.8118_wp, 3.6414_wp, 4.6373_wp, 4.6268_wp, & + 4.4750_wp, 4.4134_wp, 4.3458_wp, 3.8582_wp, 4.2583_wp, & + 4.1898_wp, 4.1562_wp, 4.1191_wp, 4.1069_wp, 4.0639_wp, & + 4.1257_wp, 4.1974_wp, 3.9532_wp, 4.1794_wp, 3.9660_wp, & + 3.8130_wp, 4.8160_wp, 4.8272_wp, 4.6294_wp, 4.5840_wp, & + 4.0770_wp, 4.0088_wp, 3.9103_wp, 3.8536_wp, 3.8324_wp, & + 3.7995_wp, 3.7826_wp, 4.2294_wp, 4.3380_wp, 4.4352_wp, & + 4.1933_wp, 4.4580_wp, 4.2554_wp, 4.1072_wp, 5.0454_wp, & + 5.1814_wp, 3.0632_wp, 3.2662_wp, 3.6432_wp, 3.8088_wp, & + 3.7910_wp, 3.7381_wp, 3.5093_wp, 3.5155_wp, 3.3047_wp, & + 3.1681_wp, 3.7871_wp, 3.9924_wp, 4.0637_wp, 4.1382_wp, & + 3.8591_wp, 4.0164_wp, 3.7878_wp, 3.6316_wp, 4.1741_wp, & + 4.3166_wp, 4.2395_wp, 4.1831_wp, 4.1107_wp, 3.5857_wp, & + 4.0270_wp, 3.9676_wp, 3.9463_wp, 3.9150_wp, 3.9021_wp, & + 3.8708_wp, 4.0240_wp, 4.1551_wp, 3.9108_wp, 4.1337_wp, & + 3.9289_wp, 3.7873_wp, 4.3666_wp, 4.5080_wp, 4.4232_wp, & + 4.3155_wp, 3.8461_wp, 3.8007_wp, 3.6991_wp, 3.6447_wp, & + 3.6308_wp, 3.5959_wp, 3.5749_wp, 4.0359_wp, 4.3124_wp, & + 4.3539_wp, 4.1122_wp, 4.3772_wp, 4.1785_wp, 4.0386_wp, & + 4.7004_wp, 4.8604_wp, 4.6261_wp, 2.9455_wp, 3.2470_wp, & + 3.6108_wp, 3.8522_wp, 3.6625_wp, 3.6598_wp, 3.4411_wp, & + 3.4660_wp, 3.2415_wp, 3.0944_wp, 3.7514_wp, 4.0397_wp, & + 3.9231_wp, 4.0561_wp, 3.7860_wp, 3.9845_wp, 3.7454_wp, & + 3.5802_wp, 4.1366_wp, 4.3581_wp, 4.2351_wp, 4.2011_wp, & + 4.1402_wp, 3.5381_wp, 4.0653_wp, 4.0093_wp, 3.9883_wp, & + 3.9570_wp, 3.9429_wp, 3.9112_wp, 3.8728_wp, 4.0682_wp, & + 3.8351_wp, 4.1054_wp, 3.8928_wp, 3.7445_wp, 4.3415_wp, & + 4.5497_wp, 4.3833_wp, 4.3122_wp, 3.8051_wp, 3.7583_wp, & + 3.6622_wp, 3.6108_wp, 3.5971_wp, 3.5628_wp, 3.5408_wp, & + 4.0780_wp, 4.0727_wp, 4.2836_wp, 4.0553_wp, 4.3647_wp, & + 4.1622_wp, 4.0178_wp, 4.5802_wp, 4.9125_wp, 4.5861_wp, & + 4.6201_wp, 2.9244_wp, 3.2241_wp, 3.5848_wp, 3.8293_wp, & + 3.6395_wp, 3.6400_wp, 3.4204_wp, 3.4499_wp, 3.2253_wp, & + 3.0779_wp, 3.7257_wp, 4.0170_wp, 3.9003_wp, 4.0372_wp, & + 3.7653_wp, 3.9672_wp, 3.7283_wp, 3.5630_wp, 4.1092_wp, & + 4.3347_wp, 4.2117_wp, 4.1793_wp, 4.1179_wp, 3.5139_wp, & + 4.0426_wp, 3.9867_wp, 3.9661_wp, 3.9345_wp, 3.9200_wp, & + 3.8883_wp, 3.8498_wp, 4.0496_wp, 3.8145_wp, 4.0881_wp, & + 3.8756_wp, 3.7271_wp, 4.3128_wp, 4.5242_wp, 4.3578_wp, & + 4.2870_wp, 3.7796_wp, 3.7318_wp, 3.6364_wp, 3.5854_wp, & + 3.5726_wp, 3.5378_wp, 3.5155_wp, 4.0527_wp, 4.0478_wp, & + 4.2630_wp, 4.0322_wp, 4.3449_wp, 4.1421_wp, 3.9975_wp, & + 4.5499_wp, 4.8825_wp, 4.5601_wp, 4.5950_wp, 4.5702_wp, & + 2.9046_wp, 3.2044_wp, 3.5621_wp, 3.8078_wp, 3.6185_wp, & + 3.6220_wp, 3.4019_wp, 3.4359_wp, 3.2110_wp, 3.0635_wp, & + 3.7037_wp, 3.9958_wp, 3.8792_wp, 4.0194_wp, 3.7460_wp, & + 3.9517_wp, 3.7128_wp, 3.5474_wp, 4.0872_wp, 4.3138_wp, & + 4.1906_wp, 4.1593_wp, 4.0973_wp, 3.4919_wp, 4.0216_wp, & + 3.9657_wp, 3.9454_wp, 3.9134_wp, 3.8986_wp, 3.8669_wp, & + 3.8289_wp, 4.0323_wp, 3.7954_wp, 4.0725_wp, 3.8598_wp, & + 3.7113_wp, 4.2896_wp, 4.5021_wp, 4.3325_wp, 4.2645_wp, & + 3.7571_wp, 3.7083_wp, 3.6136_wp, 3.5628_wp, 3.5507_wp, & + 3.5155_wp, 3.4929_wp, 4.0297_wp, 4.0234_wp, 4.2442_wp, & + 4.0112_wp, 4.3274_wp, 4.1240_wp, 3.9793_wp, 4.5257_wp, & + 4.8568_wp, 4.5353_wp, 4.5733_wp, 4.5485_wp, 4.5271_wp, & + 2.8878_wp, 3.1890_wp, 3.5412_wp, 3.7908_wp, 3.5974_wp, & + 3.6078_wp, 3.3871_wp, 3.4243_wp, 3.1992_wp, 3.0513_wp, & + 3.6831_wp, 3.9784_wp, 3.8579_wp, 4.0049_wp, 3.7304_wp, & + 3.9392_wp, 3.7002_wp, 3.5347_wp, 4.0657_wp, 4.2955_wp, & + 4.1705_wp, 4.1424_wp, 4.0800_wp, 3.4717_wp, 4.0043_wp, & + 3.9485_wp, 3.9286_wp, 3.8965_wp, 3.8815_wp, 3.8500_wp, & + 3.8073_wp, 4.0180_wp, 3.7796_wp, 4.0598_wp, 3.8470_wp, & + 3.6983_wp, 4.2678_wp, 4.4830_wp, 4.3132_wp, 4.2444_wp, & + 3.7370_wp, 3.6876_wp, 3.5935_wp, 3.5428_wp, 3.5314_wp, & + 3.4958_wp, 3.4730_wp, 4.0117_wp, 4.0043_wp, 4.2287_wp, & + 3.9939_wp, 4.3134_wp, 4.1096_wp, 3.9646_wp, 4.5032_wp, & + 4.8356_wp, 4.5156_wp, 4.5544_wp, 4.5297_wp, 4.5083_wp, & + 4.4896_wp, 2.8709_wp, 3.1737_wp, 3.5199_wp, 3.7734_wp, & + 3.5802_wp, 3.5934_wp, 3.3724_wp, 3.4128_wp, 3.1877_wp, & + 3.0396_wp, 3.6624_wp, 3.9608_wp, 3.8397_wp, 3.9893_wp, & + 3.7145_wp, 3.9266_wp, 3.6877_wp, 3.5222_wp, 4.0448_wp, & + 4.2771_wp, 4.1523_wp, 4.1247_wp, 4.0626_wp, 3.4530_wp, & + 3.9866_wp, 3.9310_wp, 3.9115_wp, 3.8792_wp, 3.8641_wp, & + 3.8326_wp, 3.7892_wp, 4.0025_wp, 3.7636_wp, 4.0471_wp, & + 3.8343_wp, 3.6854_wp, 4.2464_wp, 4.4635_wp, 4.2939_wp, & + 4.2252_wp, 3.7169_wp, 3.6675_wp, 3.5739_wp, 3.5235_wp, & + 3.5126_wp, 3.4768_wp, 3.4537_wp, 3.9932_wp, 3.9854_wp, & + 4.2123_wp, 3.9765_wp, 4.2992_wp, 4.0951_wp, 3.9500_wp, & + 4.4811_wp, 4.8135_wp, 4.4959_wp, 4.5351_wp, 4.5105_wp, & + 4.4891_wp, 4.4705_wp, 4.4515_wp, 2.8568_wp, 3.1608_wp, & + 3.5050_wp, 3.7598_wp, 3.5665_wp, 3.5803_wp, 3.3601_wp, & + 3.4031_wp, 3.1779_wp, 3.0296_wp, 3.6479_wp, 3.9471_wp, & + 3.8262_wp, 3.9773_wp, 3.7015_wp, 3.9162_wp, 3.6771_wp, & + 3.5115_wp, 4.0306_wp, 4.2634_wp, 4.1385_wp, 4.1116_wp, & + 4.0489_wp, 3.4366_wp, 3.9732_wp, 3.9176_wp, 3.8983_wp, & + 3.8659_wp, 3.8507_wp, 3.8191_wp, 3.7757_wp, 3.9907_wp, & + 3.7506_wp, 4.0365_wp, 3.8235_wp, 3.6745_wp, 4.2314_wp, & + 4.4490_wp, 4.2792_wp, 4.2105_wp, 3.7003_wp, 3.6510_wp, & + 3.5578_wp, 3.5075_wp, 3.4971_wp, 3.4609_wp, 3.4377_wp, & + 3.9788_wp, 3.9712_wp, 4.1997_wp, 3.9624_wp, 4.2877_wp, & + 4.0831_wp, 3.9378_wp, 4.4655_wp, 4.7974_wp, 4.4813_wp, & + 4.5209_wp, 4.4964_wp, 4.4750_wp, 4.4565_wp, 4.4375_wp, & + 4.4234_wp, 2.6798_wp, 3.0151_wp, 3.2586_wp, 3.5292_wp, & + 3.5391_wp, 3.4902_wp, 3.2887_wp, 3.3322_wp, 3.1228_wp, & + 2.9888_wp, 3.4012_wp, 3.7145_wp, 3.7830_wp, 3.6665_wp, & + 3.5898_wp, 3.8077_wp, 3.5810_wp, 3.4265_wp, 3.7726_wp, & + 4.0307_wp, 3.9763_wp, 3.8890_wp, 3.8489_wp, 3.2706_wp, & + 3.7595_wp, 3.6984_wp, 3.6772_wp, 3.6428_wp, 3.6243_wp, & + 3.5951_wp, 3.7497_wp, 3.6775_wp, 3.6364_wp, 3.9203_wp, & + 3.7157_wp, 3.5746_wp, 3.9494_wp, 4.2076_wp, 4.1563_wp, & + 4.0508_wp, 3.5329_wp, 3.4780_wp, 3.3731_wp, 3.3126_wp, & + 3.2846_wp, 3.2426_wp, 3.2135_wp, 3.7491_wp, 3.9006_wp, & + 3.8332_wp, 3.8029_wp, 4.1436_wp, 3.9407_wp, 3.7998_wp, & + 4.1663_wp, 4.5309_wp, 4.3481_wp, 4.2911_wp, 4.2671_wp, & + 4.2415_wp, 4.2230_wp, 4.2047_wp, 4.1908_wp, 4.1243_wp, & + 2.5189_wp, 2.9703_wp, 3.3063_wp, 3.6235_wp, 3.4517_wp, & + 3.3989_wp, 3.2107_wp, 3.2434_wp, 3.0094_wp, 2.8580_wp, & + 3.4253_wp, 3.8157_wp, 3.7258_wp, 3.6132_wp, 3.5297_wp, & + 3.7566_wp, 3.5095_wp, 3.3368_wp, 3.7890_wp, 4.1298_wp, & + 4.0190_wp, 3.9573_wp, 3.9237_wp, 3.2677_wp, 3.8480_wp, & + 3.8157_wp, 3.7656_wp, 3.7317_wp, 3.7126_wp, 3.6814_wp, & + 3.6793_wp, 3.6218_wp, 3.5788_wp, 3.8763_wp, 3.6572_wp, & + 3.5022_wp, 3.9737_wp, 4.3255_wp, 4.1828_wp, 4.1158_wp, & + 3.5078_wp, 3.4595_wp, 3.3600_wp, 3.3088_wp, 3.2575_wp, & + 3.2164_wp, 3.1856_wp, 3.8522_wp, 3.8665_wp, 3.8075_wp, & + 3.7772_wp, 4.1391_wp, 3.9296_wp, 3.7772_wp, 4.2134_wp, & + 4.7308_wp, 4.3787_wp, 4.3894_wp, 4.3649_wp, 4.3441_wp, & + 4.3257_wp, 4.3073_wp, 4.2941_wp, 4.1252_wp, 4.2427_wp, & + 3.0481_wp, 2.9584_wp, 3.6919_wp, 3.5990_wp, 3.8881_wp, & + 3.4209_wp, 3.1606_wp, 3.1938_wp, 2.9975_wp, 2.8646_wp, & + 3.8138_wp, 3.7935_wp, 3.7081_wp, 3.9155_wp, 3.5910_wp, & + 3.4808_wp, 3.4886_wp, 3.3397_wp, 4.1336_wp, 4.1122_wp, & + 3.9888_wp, 3.9543_wp, 3.8917_wp, 3.5894_wp, 3.8131_wp, & + 3.7635_wp, 3.7419_wp, 3.7071_wp, 3.6880_wp, 3.6574_wp, & + 3.6546_wp, 3.9375_wp, 3.6579_wp, 3.5870_wp, 3.6361_wp, & + 3.5039_wp, 4.3149_wp, 4.2978_wp, 4.1321_wp, 4.1298_wp, & + 3.8164_wp, 3.7680_wp, 3.7154_wp, 3.6858_wp, 3.6709_wp, & + 3.6666_wp, 3.6517_wp, 3.8174_wp, 3.8608_wp, 4.1805_wp, & + 3.9102_wp, 3.8394_wp, 3.8968_wp, 3.7673_wp, 4.5274_wp, & + 4.6682_wp, 4.3344_wp, 4.3639_wp, 4.3384_wp, 4.3162_wp, & + 4.2972_wp, 4.2779_wp, 4.2636_wp, 4.0253_wp, 4.1168_wp, & + 4.1541_wp, 2.8136_wp, 3.0951_wp, 3.4635_wp, 3.6875_wp, & + 3.4987_wp, 3.5183_wp, 3.2937_wp, 3.3580_wp, 3.1325_wp, & + 2.9832_wp, 3.6078_wp, 3.8757_wp, 3.7616_wp, 3.9222_wp, & + 3.6370_wp, 3.8647_wp, 3.6256_wp, 3.4595_wp, 3.9874_wp, & + 4.1938_wp, 4.0679_wp, 4.0430_wp, 3.9781_wp, 3.3886_wp, & + 3.9008_wp, 3.8463_wp, 3.8288_wp, 3.7950_wp, 3.7790_wp, & + 3.7472_wp, 3.7117_wp, 3.9371_wp, 3.6873_wp, 3.9846_wp, & + 3.7709_wp, 3.6210_wp, 4.1812_wp, 4.3750_wp, 4.2044_wp, & + 4.1340_wp, 3.6459_wp, 3.5929_wp, 3.5036_wp, 3.4577_wp, & + 3.4528_wp, 3.4146_wp, 3.3904_wp, 3.9014_wp, 3.9031_wp, & + 4.1443_wp, 3.8961_wp, 4.2295_wp, 4.0227_wp, 3.8763_wp, & + 4.4086_wp, 4.7097_wp, 4.4064_wp, 4.4488_wp, 4.4243_wp, & + 4.4029_wp, 4.3842_wp, 4.3655_wp, 4.3514_wp, 4.1162_wp, & + 4.2205_wp, 4.1953_wp, 4.2794_wp, 2.8032_wp, 3.0805_wp, & + 3.4519_wp, 3.6700_wp, 3.4827_wp, 3.5050_wp, 3.2799_wp, & + 3.3482_wp, 3.1233_wp, 2.9747_wp, 3.5971_wp, 3.8586_wp, & + 3.7461_wp, 3.9100_wp, 3.6228_wp, 3.8535_wp, 3.6147_wp, & + 3.4490_wp, 3.9764_wp, 4.1773_wp, 4.0511_wp, 4.0270_wp, & + 3.9614_wp, 3.3754_wp, 3.8836_wp, 3.8291_wp, 3.8121_wp, & + 3.7780_wp, 3.7619_wp, 3.7300_wp, 3.6965_wp, 3.9253_wp, & + 3.6734_wp, 3.9733_wp, 3.7597_wp, 3.6099_wp, 4.1683_wp, & + 4.3572_wp, 4.1862_wp, 4.1153_wp, 3.6312_wp, 3.5772_wp, & + 3.4881_wp, 3.4429_wp, 3.4395_wp, 3.4009_wp, 3.3766_wp, & + 3.8827_wp, 3.8868_wp, 4.1316_wp, 3.8807_wp, 4.2164_wp, & + 4.0092_wp, 3.8627_wp, 4.3936_wp, 4.6871_wp, 4.3882_wp, & + 4.4316_wp, 4.4073_wp, 4.3858_wp, 4.3672_wp, 4.3485_wp, & + 4.3344_wp, 4.0984_wp, 4.2036_wp, 4.1791_wp, 4.2622_wp, & + 4.2450_wp, 2.7967_wp, 3.0689_wp, 3.4445_wp, 3.6581_wp, & + 3.4717_wp, 3.4951_wp, 3.2694_wp, 3.3397_wp, 3.1147_wp, & + 2.9661_wp, 3.5898_wp, 3.8468_wp, 3.7358_wp, 3.9014_wp, & + 3.6129_wp, 3.8443_wp, 3.6054_wp, 3.4396_wp, 3.9683_wp, & + 4.1656_wp, 4.0394_wp, 4.0158_wp, 3.9498_wp, 3.3677_wp, & + 3.8718_wp, 3.8164_wp, 3.8005_wp, 3.7662_wp, 3.7500_wp, & + 3.7181_wp, 3.6863_wp, 3.9170_wp, 3.6637_wp, 3.9641_wp, & + 3.7503_wp, 3.6004_wp, 4.1590_wp, 4.3448_wp, 4.1739_wp, & + 4.1029_wp, 3.6224_wp, 3.5677_wp, 3.4785_wp, 3.4314_wp, & + 3.4313_wp, 3.3923_wp, 3.3680_wp, 3.8698_wp, 3.8758_wp, & + 4.1229_wp, 3.8704_wp, 4.2063_wp, 3.9987_wp, 3.8519_wp, & + 4.3832_wp, 4.6728_wp, 4.3759_wp, 4.4195_wp, 4.3952_wp, & + 4.3737_wp, 4.3551_wp, 4.3364_wp, 4.3223_wp, 4.0861_wp, & + 4.1911_wp, 4.1676_wp, 4.2501_wp, 4.2329_wp, 4.2208_wp, & + 2.7897_wp, 3.0636_wp, 3.4344_wp, 3.6480_wp, 3.4626_wp, & + 3.4892_wp, 3.2626_wp, 3.3344_wp, 3.1088_wp, 2.9597_wp, & + 3.5804_wp, 3.8359_wp, 3.7251_wp, 3.8940_wp, 3.6047_wp, & + 3.8375_wp, 3.5990_wp, 3.4329_wp, 3.9597_wp, 4.1542_wp, & + 4.0278_wp, 4.0048_wp, 3.9390_wp, 3.3571_wp, 3.8608_wp, & + 3.8056_wp, 3.7899_wp, 3.7560_wp, 3.7400_wp, 3.7081_wp, & + 3.6758_wp, 3.9095_wp, 3.6552_wp, 3.9572_wp, 3.7436_wp, & + 3.5933_wp, 4.1508_wp, 4.3337_wp, 4.1624_wp, 4.0916_wp, & + 3.6126_wp, 3.5582_wp, 3.4684_wp, 3.4212_wp, 3.4207_wp, & + 3.3829_wp, 3.3586_wp, 3.8604_wp, 3.8658_wp, 4.1156_wp, & + 3.8620_wp, 4.1994_wp, 3.9917_wp, 3.8446_wp, 4.3750_wp, & + 4.6617_wp, 4.3644_wp, 4.4083_wp, 4.3840_wp, 4.3625_wp, & + 4.3439_wp, 4.3253_wp, 4.3112_wp, 4.0745_wp, 4.1807_wp, & + 4.1578_wp, 4.2390_wp, 4.2218_wp, 4.2097_wp, 4.1986_wp, & + 2.8395_wp, 3.0081_wp, 3.3171_wp, 3.4878_wp, 3.5360_wp, & + 3.5145_wp, 3.2809_wp, 3.3307_wp, 3.1260_wp, 2.9940_wp, & + 3.4741_wp, 3.6675_wp, 3.7832_wp, 3.6787_wp, 3.6156_wp, & + 3.8041_wp, 3.5813_wp, 3.4301_wp, 3.8480_wp, 3.9849_wp, & + 3.9314_wp, 3.8405_wp, 3.8029_wp, 3.2962_wp, 3.7104_wp, & + 3.6515_wp, 3.6378_wp, 3.6020_wp, 3.5849_wp, 3.5550_wp, & + 3.7494_wp, 3.6893_wp, 3.6666_wp, 3.9170_wp, 3.7150_wp, & + 3.5760_wp, 4.0268_wp, 4.1596_wp, 4.1107_wp, 3.9995_wp, & + 3.5574_wp, 3.5103_wp, 3.4163_wp, 3.3655_wp, 3.3677_wp, & + 3.3243_wp, 3.2975_wp, 3.7071_wp, 3.9047_wp, 3.8514_wp, & + 3.8422_wp, 3.8022_wp, 3.9323_wp, 3.7932_wp, 4.2343_wp, & + 4.4583_wp, 4.3115_wp, 4.2457_wp, 4.2213_wp, 4.1945_wp, & + 4.1756_wp, 4.1569_wp, 4.1424_wp, 4.0620_wp, 4.0494_wp, & + 3.9953_wp, 4.0694_wp, 4.0516_wp, 4.0396_wp, 4.0280_wp, & + 4.0130_wp, 2.9007_wp, 2.9674_wp, 3.8174_wp, 3.5856_wp, & + 3.6486_wp, 3.5339_wp, 3.2832_wp, 3.3154_wp, 3.1144_wp, & + 2.9866_wp, 3.9618_wp, 3.8430_wp, 3.9980_wp, 3.8134_wp, & + 3.6652_wp, 3.7985_wp, 3.5756_wp, 3.4207_wp, 4.4061_wp, & + 4.2817_wp, 4.1477_wp, 4.0616_wp, 3.9979_wp, 3.6492_wp, & + 3.8833_wp, 3.8027_wp, 3.7660_wp, 3.7183_wp, 3.6954_wp, & + 3.6525_wp, 3.9669_wp, 3.8371_wp, 3.7325_wp, 3.9160_wp, & + 3.7156_wp, 3.5714_wp, 4.6036_wp, 4.4620_wp, 4.3092_wp, & + 4.2122_wp, 3.8478_wp, 3.7572_wp, 3.6597_wp, 3.5969_wp, & + 3.5575_wp, 3.5386_wp, 3.5153_wp, 3.7818_wp, 4.1335_wp, & + 4.0153_wp, 3.9177_wp, 3.8603_wp, 3.9365_wp, 3.7906_wp, & + 4.7936_wp, 4.7410_wp, 4.5461_wp, 4.5662_wp, 4.5340_wp, & + 4.5059_wp, 4.4832_wp, 4.4604_wp, 4.4429_wp, 4.2346_wp, & + 4.4204_wp, 4.3119_wp, 4.3450_wp, 4.3193_wp, 4.3035_wp, & + 4.2933_wp, 4.1582_wp, 4.2450_wp, 2.8559_wp, 2.9050_wp, & + 3.8325_wp, 3.5442_wp, 3.5077_wp, 3.4905_wp, 3.2396_wp, & + 3.2720_wp, 3.0726_wp, 2.9467_wp, 3.9644_wp, 3.8050_wp, & + 3.8981_wp, 3.7762_wp, 3.6216_wp, 3.7531_wp, 3.5297_wp, & + 3.3742_wp, 4.3814_wp, 4.2818_wp, 4.1026_wp, 4.0294_wp, & + 3.9640_wp, 3.6208_wp, 3.8464_wp, 3.7648_wp, 3.7281_wp, & + 3.6790_wp, 3.6542_wp, 3.6117_wp, 3.8650_wp, 3.8010_wp, & + 3.6894_wp, 3.8713_wp, 3.6699_wp, 3.5244_wp, 4.5151_wp, & + 4.4517_wp, 4.2538_wp, 4.1483_wp, 3.8641_wp, 3.7244_wp, & + 3.6243_wp, 3.5589_wp, 3.5172_wp, 3.4973_wp, 3.4715_wp, & + 3.7340_wp, 4.0316_wp, 3.9958_wp, 3.8687_wp, 3.8115_wp, & + 3.8862_wp, 3.7379_wp, 4.7091_wp, 4.7156_wp, 4.5199_wp, & + 4.5542_wp, 4.5230_wp, 4.4959_wp, 4.4750_wp, 4.4529_wp, & + 4.4361_wp, 4.1774_wp, 4.3774_wp, 4.2963_wp, 4.3406_wp, & + 4.3159_wp, 4.3006_wp, 4.2910_wp, 4.1008_wp, 4.1568_wp, & + 4.0980_wp, 2.8110_wp, 2.8520_wp, 3.7480_wp, 3.5105_wp, & + 3.4346_wp, 3.3461_wp, 3.1971_wp, 3.2326_wp, 3.0329_wp, & + 2.9070_wp, 3.8823_wp, 3.7928_wp, 3.8264_wp, 3.7006_wp, & + 3.5797_wp, 3.7141_wp, 3.4894_wp, 3.3326_wp, 4.3048_wp, & + 4.2217_wp, 4.0786_wp, 3.9900_wp, 3.9357_wp, 3.6331_wp, & + 3.8333_wp, 3.7317_wp, 3.6957_wp, 3.6460_wp, 3.6197_wp, & + 3.5779_wp, 3.7909_wp, 3.7257_wp, 3.6476_wp, 3.5729_wp, & + 3.6304_wp, 3.4834_wp, 4.4368_wp, 4.3921_wp, 4.2207_wp, & + 4.1133_wp, 3.8067_wp, 3.7421_wp, 3.6140_wp, 3.5491_wp, & + 3.5077_wp, 3.4887_wp, 3.4623_wp, 3.6956_wp, 3.9568_wp, & + 3.8976_wp, 3.8240_wp, 3.7684_wp, 3.8451_wp, 3.6949_wp, & + 4.6318_wp, 4.6559_wp, 4.4533_wp, 4.4956_wp, 4.4641_wp, & + 4.4366_wp, 4.4155_wp, 4.3936_wp, 4.3764_wp, 4.1302_wp, & + 4.3398_wp, 4.2283_wp, 4.2796_wp, 4.2547_wp, 4.2391_wp, & + 4.2296_wp, 4.0699_wp, 4.1083_wp, 4.0319_wp, 3.9855_wp, & + 2.7676_wp, 2.8078_wp, 3.6725_wp, 3.4804_wp, 3.3775_wp, & + 3.2411_wp, 3.1581_wp, 3.1983_wp, 2.9973_wp, 2.8705_wp, & + 3.8070_wp, 3.7392_wp, 3.7668_wp, 3.6263_wp, 3.5402_wp, & + 3.6807_wp, 3.4545_wp, 3.2962_wp, 4.2283_wp, 4.1698_wp, & + 4.0240_wp, 3.9341_wp, 3.8711_wp, 3.5489_wp, 3.7798_wp, & + 3.7000_wp, 3.6654_wp, 3.6154_wp, 3.5882_wp, 3.5472_wp, & + 3.7289_wp, 3.6510_wp, 3.6078_wp, 3.5355_wp, 3.5963_wp, & + 3.4480_wp, 4.3587_wp, 4.3390_wp, 4.1635_wp, 4.0536_wp, & + 3.7193_wp, 3.6529_wp, 3.5512_wp, 3.4837_wp, 3.4400_wp, & + 3.4191_wp, 3.3891_wp, 3.6622_wp, 3.8934_wp, 3.8235_wp, & + 3.7823_wp, 3.7292_wp, 3.8106_wp, 3.6589_wp, 4.5535_wp, & + 4.6013_wp, 4.3961_wp, 4.4423_wp, 4.4109_wp, 4.3835_wp, & + 4.3625_wp, 4.3407_wp, 4.3237_wp, 4.0863_wp, 4.2835_wp, & + 4.1675_wp, 4.2272_wp, 4.2025_wp, 4.1869_wp, 4.1774_wp, & + 4.0126_wp, 4.0460_wp, 3.9815_wp, 3.9340_wp, 3.8955_wp, & + 2.6912_wp, 2.7604_wp, 3.6037_wp, 3.4194_wp, 3.3094_wp, & + 3.1710_wp, 3.0862_wp, 3.1789_wp, 2.9738_wp, 2.8427_wp, & + 3.7378_wp, 3.6742_wp, 3.6928_wp, 3.5512_wp, 3.4614_wp, & + 3.4087_wp, 3.4201_wp, 3.2607_wp, 4.1527_wp, 4.0977_wp, & + 3.9523_wp, 3.8628_wp, 3.8002_wp, 3.4759_wp, 3.7102_wp, & + 3.6466_wp, 3.6106_wp, 3.5580_wp, 3.5282_wp, 3.4878_wp, & + 3.6547_wp, 3.5763_wp, 3.5289_wp, 3.5086_wp, 3.5593_wp, & + 3.4099_wp, 4.2788_wp, 4.2624_wp, 4.0873_wp, 3.9770_wp, & + 3.6407_wp, 3.5743_wp, 3.5178_wp, 3.4753_wp, 3.3931_wp, & + 3.3694_wp, 3.3339_wp, 3.6002_wp, 3.8164_wp, 3.7478_wp, & + 3.7028_wp, 3.6952_wp, 3.7669_wp, 3.6137_wp, 4.4698_wp, & + 4.5488_wp, 4.3168_wp, 4.3646_wp, 4.3338_wp, 4.3067_wp, & + 4.2860_wp, 4.2645_wp, 4.2478_wp, 4.0067_wp, 4.2349_wp, & + 4.0958_wp, 4.1543_wp, 4.1302_wp, 4.1141_wp, 4.1048_wp, & + 3.9410_wp, 3.9595_wp, 3.8941_wp, 3.8465_wp, 3.8089_wp, & + 3.7490_wp, 2.7895_wp, 2.5849_wp, 3.6484_wp, 3.0162_wp, & + 3.1267_wp, 3.2125_wp, 3.0043_wp, 2.9572_wp, 2.8197_wp, & + 2.7261_wp, 3.7701_wp, 3.2446_wp, 3.5239_wp, 3.4696_wp, & + 3.4261_wp, 3.3508_wp, 3.1968_wp, 3.0848_wp, 4.1496_wp, & + 3.6598_wp, 3.5111_wp, 3.4199_wp, 3.3809_wp, 3.5382_wp, & + 3.2572_wp, 3.2100_wp, 3.1917_wp, 3.1519_wp, 3.1198_wp, & + 3.1005_wp, 3.5071_wp, 3.5086_wp, 3.5073_wp, 3.4509_wp, & + 3.3120_wp, 3.2082_wp, 4.2611_wp, 3.8117_wp, 3.6988_wp, & + 3.5646_wp, 3.6925_wp, 3.6295_wp, 3.5383_wp, 3.4910_wp, & + 3.4625_wp, 3.4233_wp, 3.4007_wp, 3.2329_wp, 3.6723_wp, & + 3.6845_wp, 3.6876_wp, 3.6197_wp, 3.4799_wp, 3.3737_wp, & + 4.4341_wp, 4.0525_wp, 3.9011_wp, 3.8945_wp, 3.8635_wp, & + 3.8368_wp, 3.8153_wp, 3.7936_wp, 3.7758_wp, 3.4944_wp, & + 3.4873_wp, 3.9040_wp, 3.7110_wp, 3.6922_wp, 3.6799_wp, & + 3.6724_wp, 3.5622_wp, 3.6081_wp, 3.5426_wp, 3.4922_wp, & + 3.4498_wp, 3.3984_wp, 3.4456_wp, 2.7522_wp, 2.5524_wp, & + 3.5742_wp, 2.9508_wp, 3.0751_wp, 3.0158_wp, 2.9644_wp, & + 2.8338_wp, 2.7891_wp, 2.6933_wp, 3.6926_wp, 3.1814_wp, & + 3.4528_wp, 3.4186_wp, 3.3836_wp, 3.2213_wp, 3.1626_wp, & + 3.0507_wp, 4.0548_wp, 3.5312_wp, 3.4244_wp, 3.3409_wp, & + 3.2810_wp, 3.4782_wp, 3.1905_wp, 3.1494_wp, 3.1221_wp, & + 3.1128_wp, 3.0853_wp, 3.0384_wp, 3.4366_wp, 3.4562_wp, & + 3.4638_wp, 3.3211_wp, 3.2762_wp, 3.1730_wp, 4.1632_wp, & + 3.6825_wp, 3.5822_wp, 3.4870_wp, 3.6325_wp, 3.5740_wp, & + 3.4733_wp, 3.4247_wp, 3.3969_wp, 3.3764_wp, 3.3525_wp, & + 3.1984_wp, 3.5989_wp, 3.6299_wp, 3.6433_wp, 3.4937_wp, & + 3.4417_wp, 3.3365_wp, 4.3304_wp, 3.9242_wp, 3.7793_wp, & + 3.7623_wp, 3.7327_wp, 3.7071_wp, 3.6860_wp, 3.6650_wp, & + 3.6476_wp, 3.3849_wp, 3.3534_wp, 3.8216_wp, 3.5870_wp, & + 3.5695_wp, 3.5584_wp, 3.5508_wp, 3.4856_wp, 3.5523_wp, & + 3.4934_wp, 3.4464_wp, 3.4055_wp, 3.3551_wp, 3.3888_wp, & + 3.3525_wp, 2.7202_wp, 2.5183_wp, 3.4947_wp, 2.8731_wp, & + 3.0198_wp, 3.1457_wp, 2.9276_wp, 2.7826_wp, 2.7574_wp, & + 2.6606_wp, 3.6090_wp, 3.0581_wp, 3.3747_wp, 3.3677_wp, & + 3.3450_wp, 3.1651_wp, 3.1259_wp, 3.0147_wp, 3.9498_wp, & + 3.3857_wp, 3.2917_wp, 3.2154_wp, 3.1604_wp, 3.4174_wp, & + 3.0735_wp, 3.0342_wp, 3.0096_wp, 3.0136_wp, 2.9855_wp, & + 2.9680_wp, 3.3604_wp, 3.4037_wp, 3.4243_wp, 3.2633_wp, & + 3.1810_wp, 3.1351_wp, 4.0557_wp, 3.5368_wp, 3.4526_wp, & + 3.3699_wp, 3.5707_wp, 3.5184_wp, 3.4085_wp, 3.3595_wp, & + 3.3333_wp, 3.3143_wp, 3.3041_wp, 3.1094_wp, 3.5193_wp, & + 3.5745_wp, 3.6025_wp, 3.4338_wp, 3.3448_wp, 3.2952_wp, & + 4.2158_wp, 3.7802_wp, 3.6431_wp, 3.6129_wp, 3.5853_wp, & + 3.5610_wp, 3.5406_wp, 3.5204_wp, 3.5036_wp, 3.2679_wp, & + 3.2162_wp, 3.7068_wp, 3.4483_wp, 3.4323_wp, 3.4221_wp, & + 3.4138_wp, 3.3652_wp, 3.4576_wp, 3.4053_wp, 3.3618_wp, & + 3.3224_wp, 3.2711_wp, 3.3326_wp, 3.2950_wp, 3.2564_wp, & + 2.5315_wp, 2.6104_wp, 3.2734_wp, 3.2299_wp, 3.1090_wp, & + 2.9942_wp, 2.9159_wp, 2.8324_wp, 2.8350_wp, 2.7216_wp, & + 3.3994_wp, 3.4475_wp, 3.4354_wp, 3.3438_wp, 3.2807_wp, & + 3.2169_wp, 3.2677_wp, 3.1296_wp, 3.7493_wp, 3.8075_wp, & + 3.6846_wp, 3.6104_wp, 3.5577_wp, 3.2052_wp, 3.4803_wp, & + 3.4236_wp, 3.3845_wp, 3.3640_wp, 3.3365_wp, 3.3010_wp, & + 3.3938_wp, 3.3624_wp, 3.3440_wp, 3.3132_wp, 3.4035_wp, & + 3.2754_wp, 3.8701_wp, 3.9523_wp, 3.8018_wp, 3.7149_wp, & + 3.3673_wp, 3.3199_wp, 3.2483_wp, 3.2069_wp, 3.1793_wp, & + 3.1558_wp, 3.1395_wp, 3.4097_wp, 3.5410_wp, 3.5228_wp, & + 3.5116_wp, 3.4921_wp, 3.4781_wp, 3.4690_wp, 4.0420_wp, & + 4.1759_wp, 4.0078_wp, 4.0450_wp, 4.0189_wp, 3.9952_wp, & + 3.9770_wp, 3.9583_wp, 3.9434_wp, 3.7217_wp, 3.8228_wp, & + 3.7826_wp, 3.8640_wp, 3.8446_wp, 3.8314_wp, 3.8225_wp, & + 3.6817_wp, 3.7068_wp, 3.6555_wp, 3.6159_wp, 3.5831_wp, & + 3.5257_wp, 3.2133_wp, 3.1689_wp, 3.1196_wp, 3.3599_wp, & + 2.9852_wp, 2.7881_wp, 3.5284_wp, 3.3493_wp, 3.6958_wp, & + 3.3642_wp, 3.1568_wp, 3.0055_wp, 2.9558_wp, 2.8393_wp, & + 3.6287_wp, 3.5283_wp, 4.1511_wp, 3.8259_wp, 3.6066_wp, & + 3.4527_wp, 3.3480_wp, 3.2713_wp, 3.9037_wp, 3.8361_wp, & + 3.8579_wp, 3.7311_wp, 3.6575_wp, 3.5176_wp, 3.5693_wp, & + 3.5157_wp, 3.4814_wp, 3.4559_wp, 3.4445_wp, 3.4160_wp, & + 4.1231_wp, 3.8543_wp, 3.6816_wp, 3.5602_wp, 3.4798_wp, & + 3.4208_wp, 4.0542_wp, 4.0139_wp, 4.0165_wp, 3.9412_wp, & + 3.7698_wp, 3.6915_wp, 3.6043_wp, 3.5639_wp, 3.5416_wp, & + 3.5247_wp, 3.5153_wp, 3.5654_wp, 4.2862_wp, 4.0437_wp, & + 3.8871_wp, 3.7741_wp, 3.6985_wp, 3.6413_wp, 4.2345_wp, & + 4.3663_wp, 4.3257_wp, 4.0869_wp, 4.0612_wp, 4.0364_wp, & + 4.0170_wp, 3.9978_wp, 3.9834_wp, 3.9137_wp, 3.8825_wp, & + 3.8758_wp, 3.9143_wp, 3.8976_wp, 3.8864_wp, 3.8768_wp, & + 3.9190_wp, 4.1613_wp, 4.0566_wp, 3.9784_wp, 3.9116_wp, & + 3.8326_wp, 3.7122_wp, 3.6378_wp, 3.5576_wp, 3.5457_wp, & + 4.3127_wp, 3.1160_wp, 2.8482_wp, 4.0739_wp, 3.3599_wp, & + 3.5698_wp, 3.5366_wp, 3.2854_wp, 3.1039_wp, 2.9953_wp, & + 2.9192_wp, 4.1432_wp, 3.5320_wp, 3.9478_wp, 4.0231_wp, & + 3.7509_wp, 3.5604_wp, 3.4340_wp, 3.3426_wp, 4.3328_wp, & + 3.8288_wp, 3.7822_wp, 3.7909_wp, 3.6907_wp, 3.6864_wp, & + 3.5793_wp, 3.5221_wp, 3.4883_wp, 3.4649_wp, 3.4514_wp, & + 3.4301_wp, 3.9256_wp, 4.0596_wp, 3.8307_wp, 3.6702_wp, & + 3.5651_wp, 3.4884_wp, 4.4182_wp, 4.2516_wp, 3.9687_wp, & + 3.9186_wp, 3.9485_wp, 3.8370_wp, 3.7255_wp, 3.6744_wp, & + 3.6476_wp, 3.6295_wp, 3.6193_wp, 3.5659_wp, 4.0663_wp, & + 4.2309_wp, 4.0183_wp, 3.8680_wp, 3.7672_wp, 3.6923_wp, & + 4.5240_wp, 4.4834_wp, 4.1570_wp, 4.3204_wp, 4.2993_wp, & + 4.2804_wp, 4.2647_wp, 4.2481_wp, 4.2354_wp, 3.8626_wp, & + 3.8448_wp, 4.2267_wp, 4.1799_wp, 4.1670_wp, 3.8738_wp, & + 3.8643_wp, 3.8796_wp, 4.0575_wp, 4.0354_wp, 3.9365_wp, & + 3.8611_wp, 3.7847_wp, 3.7388_wp, 3.6826_wp, 3.6251_wp, & + 3.5492_wp, 4.0889_wp, 4.2764_wp, 3.1416_wp, 2.8325_wp, & + 3.7735_wp, 3.3787_wp, 3.4632_wp, 3.5923_wp, 3.3214_wp, & + 3.1285_wp, 3.0147_wp, 2.9366_wp, 3.8527_wp, 3.5602_wp, & + 3.8131_wp, 3.8349_wp, 3.7995_wp, 3.5919_wp, 3.4539_wp, & + 3.3540_wp, 4.0654_wp, 3.8603_wp, 3.7972_wp, 3.7358_wp, & + 3.7392_wp, 3.8157_wp, 3.6055_wp, 3.5438_wp, 3.5089_wp, & + 3.4853_wp, 3.4698_wp, 3.4508_wp, 3.7882_wp, 3.8682_wp, & + 3.8837_wp, 3.7055_wp, 3.5870_wp, 3.5000_wp, 4.1573_wp, & + 4.0005_wp, 3.9568_wp, 3.8936_wp, 3.9990_wp, 3.9433_wp, & + 3.8172_wp, 3.7566_wp, 3.7246_wp, 3.7033_wp, 3.6900_wp, & + 3.5697_wp, 3.9183_wp, 4.0262_wp, 4.0659_wp, 3.8969_wp, & + 3.7809_wp, 3.6949_wp, 4.2765_wp, 4.2312_wp, 4.1401_wp, & + 4.0815_wp, 4.0580_wp, 4.0369_wp, 4.0194_wp, 4.0017_wp, & + 3.9874_wp, 3.8312_wp, 3.8120_wp, 3.9454_wp, 3.9210_wp, & + 3.9055_wp, 3.8951_wp, 3.8866_wp, 3.8689_wp, 3.9603_wp, & + 3.9109_wp, 3.9122_wp, 3.8233_wp, 3.7438_wp, 3.7436_wp, & + 3.6981_wp, 3.6555_wp, 3.5452_wp, 3.9327_wp, 4.0658_wp, & + 4.1175_wp, 2.9664_wp, 2.8209_wp, 3.5547_wp, 3.3796_wp, & + 3.3985_wp, 3.3164_wp, 3.2364_wp, 3.1956_wp, 3.0370_wp, & + 2.9313_wp, 3.6425_wp, 3.5565_wp, 3.7209_wp, 3.7108_wp, & + 3.6639_wp, 3.6484_wp, 3.4745_wp, 3.3492_wp, 3.8755_wp, & + 4.2457_wp, 3.7758_wp, 3.7161_wp, 3.6693_wp, 3.6155_wp, & + 3.5941_wp, 3.5643_wp, 3.5292_wp, 3.4950_wp, 3.4720_wp, & + 3.4503_wp, 3.6936_wp, 3.7392_wp, 3.7388_wp, 3.7602_wp, & + 3.6078_wp, 3.4960_wp, 3.9800_wp, 4.3518_wp, 4.2802_wp, & + 3.8580_wp, 3.8056_wp, 3.7527_wp, 3.7019_wp, 3.6615_wp, & + 3.5768_wp, 3.5330_wp, 3.5038_wp, 3.5639_wp, 3.8192_wp, & + 3.8883_wp, 3.9092_wp, 3.9478_wp, 3.7995_wp, 3.6896_wp, & + 4.1165_wp, 4.5232_wp, 4.4357_wp, 4.4226_wp, 4.4031_wp, & + 4.3860_wp, 4.3721_wp, 4.3580_wp, 4.3466_wp, 4.2036_wp, & + 4.2037_wp, 3.8867_wp, 4.2895_wp, 4.2766_wp, 4.2662_wp, & + 4.2598_wp, 3.8408_wp, 3.9169_wp, 3.8681_wp, 3.8250_wp, & + 3.7855_wp, 3.7501_wp, 3.6753_wp, 3.5499_wp, 3.4872_wp, & + 3.5401_wp, 3.8288_wp, 3.9217_wp, 3.9538_wp, 4.0054_wp, & + 2.8388_wp, 2.7890_wp, 3.4329_wp, 3.5593_wp, 3.3488_wp, & + 3.2486_wp, 3.1615_wp, 3.1000_wp, 3.0394_wp, 2.9165_wp, & + 3.5267_wp, 3.7479_wp, 3.6650_wp, 3.6263_wp, 3.5658_wp, & + 3.5224_wp, 3.4762_wp, 3.3342_wp, 3.7738_wp, 4.0333_wp, & + 3.9568_wp, 3.8975_wp, 3.8521_wp, 3.4929_wp, 3.7830_wp, & + 3.7409_wp, 3.7062_wp, 3.6786_wp, 3.6471_wp, 3.6208_wp, & + 3.6337_wp, 3.6519_wp, 3.6363_wp, 3.6278_wp, 3.6110_wp, & + 3.4825_wp, 3.8795_wp, 4.1448_wp, 4.0736_wp, 4.0045_wp, & + 3.6843_wp, 3.6291_wp, 3.5741_wp, 3.5312_wp, 3.4974_wp, & + 3.4472_wp, 3.4034_wp, 3.7131_wp, 3.7557_wp, 3.7966_wp, & + 3.8005_wp, 3.8068_wp, 3.8015_wp, 3.6747_wp, 4.0222_wp, & + 4.3207_wp, 4.2347_wp, 4.2191_wp, 4.1990_wp, 4.1811_wp, & + 4.1666_wp, 4.1521_wp, 4.1401_wp, 3.9970_wp, 3.9943_wp, & + 3.9592_wp, 4.0800_wp, 4.0664_wp, 4.0559_wp, 4.0488_wp, & + 3.9882_wp, 4.0035_wp, 3.9539_wp, 3.9138_wp, 3.8798_wp, & + 3.8355_wp, 3.5359_wp, 3.4954_wp, 3.3962_wp, 3.5339_wp, & + 3.7595_wp, 3.8250_wp, 3.8408_wp, 3.8600_wp, 3.8644_wp, & + 2.7412_wp, 2.7489_wp, 3.3374_wp, 3.3950_wp, 3.3076_wp, & + 3.1910_wp, 3.0961_wp, 3.0175_wp, 3.0280_wp, 2.8929_wp, & + 3.4328_wp, 3.5883_wp, 3.6227_wp, 3.5616_wp, 3.4894_wp, & + 3.4241_wp, 3.3641_wp, 3.3120_wp, 3.6815_wp, 3.8789_wp, & + 3.8031_wp, 3.7413_wp, 3.6939_wp, 3.4010_wp, 3.6225_wp, & + 3.5797_wp, 3.5443_wp, 3.5139_wp, 3.4923_wp, 3.4642_wp, & + 3.5860_wp, 3.5849_wp, 3.5570_wp, 3.5257_wp, 3.4936_wp, & + 3.4628_wp, 3.7874_wp, 3.9916_wp, 3.9249_wp, 3.8530_wp, & + 3.5932_wp, 3.5355_wp, 3.4757_wp, 3.4306_wp, 3.3953_wp, & + 3.3646_wp, 3.3390_wp, 3.5637_wp, 3.7053_wp, 3.7266_wp, & + 3.7177_wp, 3.6996_wp, 3.6775_wp, 3.6558_wp, 3.9331_wp, & + 4.1655_wp, 4.0879_wp, 4.0681_wp, 4.0479_wp, 4.0299_wp, & + 4.0152_wp, 4.0006_wp, 3.9883_wp, 3.8500_wp, 3.8359_wp, & + 3.8249_wp, 3.9269_wp, 3.9133_wp, 3.9025_wp, 3.8948_wp, & + 3.8422_wp, 3.8509_wp, 3.7990_wp, 3.7570_wp, 3.7219_wp, & + 3.6762_wp, 3.4260_wp, 3.3866_wp, 3.3425_wp, 3.5294_wp, & + 3.7022_wp, 3.7497_wp, 3.7542_wp, 3.7494_wp, 3.7370_wp, & + 3.7216_wp, 3.3286_wp, 3.5286_wp, 4.5857_wp, 4.2143_wp, & + 3.9714_wp, 3.9428_wp, 3.8857_wp, 3.7714_wp, 3.7857_wp, & + 3.8286_wp, 4.8714_wp, 4.6571_wp, 4.4857_wp, 4.3571_wp, & + 4.4429_wp, 4.3286_wp, 4.2857_wp, 4.2428_wp, 5.3857_wp, & + 5.0714_wp, 4.7714_wp, 4.6143_wp, 4.6000_wp, 4.4429_wp, & + 4.4000_wp, 4.3571_wp, 4.3000_wp, 4.2857_wp, 4.3143_wp, & + 4.4286_wp, 4.4714_wp, 4.4286_wp, 4.5143_wp, 4.4429_wp, & + 4.5000_wp, 4.5429_wp, 5.5714_wp, 5.2572_wp, 4.9714_wp, & + 4.8572_wp, 4.7571_wp, 4.6428_wp, 4.5143_wp, 4.4857_wp, & + 4.4857_wp, 4.4143_wp, 4.5143_wp, 4.6286_wp, 4.7000_wp, & + 4.6714_wp, 4.6714_wp, 4.6286_wp, 4.7571_wp, 4.7429_wp, & + 5.8571_wp, 5.3857_wp, 5.1857_wp, 4.9714_wp, 5.1286_wp, & + 5.1143_wp, 5.1000_wp, 5.0857_wp, 5.0286_wp, 5.0428_wp, & + 5.0286_wp, 5.0143_wp, 5.0000_wp, 5.0000_wp, 4.9857_wp, & + 5.0571_wp, 4.9571_wp, 4.8286_wp, 4.7429_wp, 4.6286_wp, & + 4.5572_wp, 4.5286_wp, 4.4572_wp, 4.4714_wp, 4.4857_wp, & + 4.7571_wp, 4.7286_wp, 4.7286_wp, 4.8143_wp, 4.7429_wp, & + 4.8429_wp, 4.9000_wp, 5.7428_wp, 3.0429_wp, 3.2429_wp, & + 4.3000_wp, 3.9286_wp, 3.6857_wp, 3.6571_wp, 3.6000_wp, & + 3.4857_wp, 3.5000_wp, 3.5428_wp, 4.5857_wp, 4.3714_wp, & + 4.2000_wp, 4.0714_wp, 4.1572_wp, 4.0429_wp, 4.0000_wp, & + 3.9571_wp, 5.1000_wp, 4.7857_wp, 4.4857_wp, 4.3286_wp, & + 4.3143_wp, 4.1572_wp, 4.1143_wp, 4.0714_wp, 4.0143_wp, & + 4.0000_wp, 4.0286_wp, 4.1429_wp, 4.1857_wp, 4.1429_wp, & + 4.2285_wp, 4.1572_wp, 4.2143_wp, 4.2571_wp, 5.2857_wp, & + 4.9714_wp, 4.6857_wp, 4.5714_wp, 4.4714_wp, 4.3571_wp, & + 4.2285_wp, 4.2000_wp, 4.2000_wp, 4.1286_wp, 4.2285_wp, & + 4.3429_wp, 4.4143_wp, 4.3857_wp, 4.3857_wp, 4.3429_wp, & + 4.4714_wp, 4.4572_wp, 5.5714_wp, 5.1000_wp, 4.9000_wp, & + 4.6857_wp, 4.8429_wp, 4.8286_wp, 4.8143_wp, 4.8000_wp, & + 4.7429_wp, 4.7571_wp, 4.7429_wp, 4.7286_wp, 4.7143_wp, & + 4.7143_wp, 4.7000_wp, 4.7714_wp, 4.6714_wp, 4.5429_wp, & + 4.4572_wp, 4.3429_wp, 4.2714_wp, 4.2428_wp, 4.1715_wp, & + 4.1857_wp, 4.2000_wp, 4.4714_wp, 4.4429_wp, 4.4429_wp, & + 4.5286_wp, 4.4572_wp, 4.5572_wp, 4.6143_wp, 5.4571_wp, & + 5.1714_wp, 2.8428_wp, 3.0429_wp, 4.1000_wp, 3.7286_wp, & + 3.4857_wp, 3.4572_wp, 3.4000_wp, 3.2857_wp, 3.3000_wp, & + 3.3429_wp, 4.3857_wp, 4.1715_wp, 4.0000_wp, 3.8714_wp, & + 3.9571_wp, 3.8428_wp, 3.8000_wp, 3.7572_wp, 4.9000_wp, & + 4.5857_wp, 4.2857_wp, 4.1286_wp, 4.1143_wp, 3.9571_wp, & + 3.9143_wp, 3.8714_wp, 3.8143_wp, 3.8000_wp, 3.8286_wp, & + 3.9428_wp, 3.9857_wp, 3.9428_wp, 4.0286_wp, 3.9571_wp, & + 4.0143_wp, 4.0571_wp, 5.0857_wp, 4.7714_wp, 4.4857_wp, & + 4.3714_wp, 4.2714_wp, 4.1572_wp, 4.0286_wp, 4.0000_wp, & + 4.0000_wp, 3.9286_wp, 4.0286_wp, 4.1429_wp, 4.2143_wp, & + 4.1857_wp, 4.1857_wp, 4.1429_wp, 4.2714_wp, 4.2571_wp, & + 5.3714_wp, 4.9000_wp, 4.7000_wp, 4.4857_wp, 4.6428_wp, & + 4.6286_wp, 4.6143_wp, 4.6000_wp, 4.5429_wp, 4.5572_wp, & + 4.5429_wp, 4.5286_wp, 4.5143_wp, 4.5143_wp, 4.5000_wp, & + 4.5714_wp, 4.4714_wp, 4.3429_wp, 4.2571_wp, 4.1429_wp, & + 4.0714_wp, 4.0429_wp, 3.9714_wp, 3.9857_wp, 4.0000_wp, & + 4.2714_wp, 4.2428_wp, 4.2428_wp, 4.3286_wp, 4.2571_wp, & + 4.3571_wp, 4.4143_wp, 5.2572_wp, 4.9714_wp, 4.7714_wp, & + 2.7143_wp, 2.9143_wp, 3.9714_wp, 3.6000_wp, 3.3572_wp, & + 3.3286_wp, 3.2714_wp, 3.1571_wp, 3.1714_wp, 3.2143_wp, & + 4.2571_wp, 4.0429_wp, 3.8714_wp, 3.7429_wp, 3.8286_wp, & + 3.7143_wp, 3.6714_wp, 3.6286_wp, 4.7714_wp, 4.4572_wp, & + 4.1572_wp, 4.0000_wp, 3.9857_wp, 3.8286_wp, 3.7857_wp, & + 3.7429_wp, 3.6857_wp, 3.6714_wp, 3.7000_wp, 3.8143_wp, & + 3.8571_wp, 3.8143_wp, 3.9000_wp, 3.8286_wp, 3.8857_wp, & + 3.9286_wp, 4.9571_wp, 4.6428_wp, 4.3571_wp, 4.2428_wp, & + 4.1429_wp, 4.0286_wp, 3.9000_wp, 3.8714_wp, 3.8714_wp, & + 3.8000_wp, 3.9000_wp, 4.0143_wp, 4.0857_wp, 4.0571_wp, & + 4.0571_wp, 4.0143_wp, 4.1429_wp, 4.1286_wp, 5.2429_wp, & + 4.7714_wp, 4.5714_wp, 4.3571_wp, 4.5143_wp, 4.5000_wp, & + 4.4857_wp, 4.4714_wp, 4.4143_wp, 4.4286_wp, 4.4143_wp, & + 4.4000_wp, 4.3857_wp, 4.3857_wp, 4.3714_wp, 4.4429_wp, & + 4.3429_wp, 4.2143_wp, 4.1286_wp, 4.0143_wp, 3.9428_wp, & + 3.9143_wp, 3.8428_wp, 3.8571_wp, 3.8714_wp, 4.1429_wp, & + 4.1143_wp, 4.1143_wp, 4.2000_wp, 4.1286_wp, 4.2285_wp, & + 4.2857_wp, 5.1286_wp, 4.8429_wp, 4.6428_wp, 4.5143_wp, & + 2.6286_wp, 2.8286_wp, 3.8857_wp, 3.5143_wp, 3.2714_wp, & + 3.2429_wp, 3.1857_wp, 3.0715_wp, 3.0857_wp, 3.1285_wp, & + 4.1715_wp, 3.9571_wp, 3.7857_wp, 3.6571_wp, 3.7429_wp, & + 3.6286_wp, 3.5857_wp, 3.5428_wp, 4.6857_wp, 4.3714_wp, & + 4.0714_wp, 3.9143_wp, 3.9000_wp, 3.7429_wp, 3.7000_wp, & + 3.6571_wp, 3.6000_wp, 3.5857_wp, 3.6143_wp, 3.7286_wp, & + 3.7714_wp, 3.7286_wp, 3.8143_wp, 3.7429_wp, 3.8000_wp, & + 3.8428_wp, 4.8714_wp, 4.5572_wp, 4.2714_wp, 4.1572_wp, & + 4.0571_wp, 3.9428_wp, 3.8143_wp, 3.7857_wp, 3.7857_wp, & + 3.7143_wp, 3.8143_wp, 3.9286_wp, 4.0000_wp, 3.9714_wp, & + 3.9714_wp, 3.9286_wp, 4.0571_wp, 4.0429_wp, 5.1571_wp, & + 4.6857_wp, 4.4857_wp, 4.2714_wp, 4.4286_wp, 4.4143_wp, & + 4.4000_wp, 4.3857_wp, 4.3286_wp, 4.3429_wp, 4.3286_wp, & + 4.3143_wp, 4.3000_wp, 4.3000_wp, 4.2857_wp, 4.3571_wp, & + 4.2571_wp, 4.1286_wp, 4.0429_wp, 3.9286_wp, 3.8571_wp, & + 3.8286_wp, 3.7572_wp, 3.7714_wp, 3.7857_wp, 4.0571_wp, & + 4.0286_wp, 4.0286_wp, 4.1143_wp, 4.0429_wp, 4.1429_wp, & + 4.2000_wp, 5.0428_wp, 4.7571_wp, 4.5572_wp, 4.4286_wp, & + 4.3429_wp, 2.6429_wp, 2.8428_wp, 3.9000_wp, 3.5286_wp, & + 3.2857_wp, 3.2571_wp, 3.2000_wp, 3.0857_wp, 3.1000_wp, & + 3.1428_wp, 4.1857_wp, 3.9714_wp, 3.8000_wp, 3.6714_wp, & + 3.7572_wp, 3.6429_wp, 3.6000_wp, 3.5571_wp, 4.7000_wp, & + 4.3857_wp, 4.0857_wp, 3.9286_wp, 3.9143_wp, 3.7572_wp, & + 3.7143_wp, 3.6714_wp, 3.6143_wp, 3.6000_wp, 3.6286_wp, & + 3.7429_wp, 3.7857_wp, 3.7429_wp, 3.8286_wp, 3.7572_wp, & + 3.8143_wp, 3.8571_wp, 4.8857_wp, 4.5714_wp, 4.2857_wp, & + 4.1715_wp, 4.0714_wp, 3.9571_wp, 3.8286_wp, 3.8000_wp, & + 3.8000_wp, 3.7286_wp, 3.8286_wp, 3.9428_wp, 4.0143_wp, & + 3.9857_wp, 3.9857_wp, 3.9428_wp, 4.0714_wp, 4.0571_wp, & + 5.1714_wp, 4.7000_wp, 4.5000_wp, 4.2857_wp, 4.4429_wp, & + 4.4286_wp, 4.4143_wp, 4.4000_wp, 4.3429_wp, 4.3571_wp, & + 4.3429_wp, 4.3286_wp, 4.3143_wp, 4.3143_wp, 4.3000_wp, & + 4.3714_wp, 4.2714_wp, 4.1429_wp, 4.0571_wp, 3.9428_wp, & + 3.8714_wp, 3.8428_wp, 3.7714_wp, 3.7857_wp, 3.8000_wp, & + 4.0714_wp, 4.0429_wp, 4.0429_wp, 4.1286_wp, 4.0571_wp, & + 4.1572_wp, 4.2143_wp, 5.0571_wp, 4.7714_wp, 4.5714_wp, & + 4.4429_wp, 4.3571_wp, 4.3714_wp, 2.6572_wp, 2.8571_wp, & + 3.9143_wp, 3.5428_wp, 3.3000_wp, 3.2714_wp, 3.2143_wp, & + 3.1000_wp, 3.1143_wp, 3.1571_wp, 4.2000_wp, 3.9857_wp, & + 3.8143_wp, 3.6857_wp, 3.7714_wp, 3.6571_wp, 3.6143_wp, & + 3.5714_wp, 4.7143_wp, 4.4000_wp, 4.1000_wp, 3.9428_wp, & + 3.9286_wp, 3.7714_wp, 3.7286_wp, 3.6857_wp, 3.6286_wp, & + 3.6143_wp, 3.6429_wp, 3.7572_wp, 3.8000_wp, 3.7572_wp, & + 3.8428_wp, 3.7714_wp, 3.8286_wp, 3.8714_wp, 4.9000_wp, & + 4.5857_wp, 4.3000_wp, 4.1857_wp, 4.0857_wp, 3.9714_wp, & + 3.8428_wp, 3.8143_wp, 3.8143_wp, 3.7429_wp, 3.8428_wp, & + 3.9571_wp, 4.0286_wp, 4.0000_wp, 4.0000_wp, 3.9571_wp, & + 4.0857_wp, 4.0714_wp, 5.1857_wp, 4.7143_wp, 4.5143_wp, & + 4.3000_wp, 4.4572_wp, 4.4429_wp, 4.4286_wp, 4.4143_wp, & + 4.3571_wp, 4.3714_wp, 4.3571_wp, 4.3429_wp, 4.3286_wp, & + 4.3286_wp, 4.3143_wp, 4.3857_wp, 4.2857_wp, 4.1572_wp, & + 4.0714_wp, 3.9571_wp, 3.8857_wp, 3.8571_wp, 3.7857_wp, & + 3.8000_wp, 3.8143_wp, 4.0857_wp, 4.0571_wp, 4.0571_wp, & + 4.1429_wp, 4.0714_wp, 4.1715_wp, 4.2285_wp, 5.0714_wp, & + 4.7857_wp, 4.5857_wp, 4.4572_wp, 4.3714_wp, 4.3857_wp, & + 4.4000_wp, 2.6714_wp, 2.8714_wp, 3.9286_wp, 3.5571_wp, & + 3.3143_wp, 3.2857_wp, 3.2286_wp, 3.1143_wp, 3.1285_wp, & + 3.1714_wp, 4.2143_wp, 4.0000_wp, 3.8286_wp, 3.7000_wp, & + 3.7857_wp, 3.6714_wp, 3.6286_wp, 3.5857_wp, 4.7286_wp, & + 4.4143_wp, 4.1143_wp, 3.9571_wp, 3.9428_wp, 3.7857_wp, & + 3.7429_wp, 3.7000_wp, 3.6429_wp, 3.6286_wp, 3.6571_wp, & + 3.7714_wp, 3.8143_wp, 3.7714_wp, 3.8571_wp, 3.7857_wp, & + 3.8428_wp, 3.8857_wp, 4.9143_wp, 4.6000_wp, 4.3143_wp, & + 4.2000_wp, 4.1000_wp, 3.9857_wp, 3.8571_wp, 3.8286_wp, & + 3.8286_wp, 3.7572_wp, 3.8571_wp, 3.9714_wp, 4.0429_wp, & + 4.0143_wp, 4.0143_wp, 3.9714_wp, 4.1000_wp, 4.0857_wp, & + 5.2000_wp, 4.7286_wp, 4.5286_wp, 4.3143_wp, 4.4714_wp, & + 4.4572_wp, 4.4429_wp, 4.4286_wp, 4.3714_wp, 4.3857_wp, & + 4.3714_wp, 4.3571_wp, 4.3429_wp, 4.3429_wp, 4.3286_wp, & + 4.4000_wp, 4.3000_wp, 4.1715_wp, 4.0857_wp, 3.9714_wp, & + 3.9000_wp, 3.8714_wp, 3.8000_wp, 3.8143_wp, 3.8286_wp, & + 4.1000_wp, 4.0714_wp, 4.0714_wp, 4.1572_wp, 4.0857_wp, & + 4.1857_wp, 4.2428_wp, 5.0857_wp, 4.8000_wp, 4.6000_wp, & + 4.4714_wp, 4.3857_wp, 4.4000_wp, 4.4143_wp, 4.4286_wp, & + 2.5857_wp, 2.7857_wp, 3.8428_wp, 3.4714_wp, 3.2286_wp, & + 3.2000_wp, 3.1428_wp, 3.0286_wp, 3.0429_wp, 3.0857_wp, & + 4.1286_wp, 3.9143_wp, 3.7429_wp, 3.6143_wp, 3.7000_wp, & + 3.5857_wp, 3.5428_wp, 3.5000_wp, 4.6428_wp, 4.3286_wp, & + 4.0286_wp, 3.8714_wp, 3.8571_wp, 3.7000_wp, 3.6571_wp, & + 3.6143_wp, 3.5571_wp, 3.5428_wp, 3.5714_wp, 3.6857_wp, & + 3.7286_wp, 3.6857_wp, 3.7714_wp, 3.7000_wp, 3.7572_wp, & + 3.8000_wp, 4.8286_wp, 4.5143_wp, 4.2285_wp, 4.1143_wp, & + 4.0143_wp, 3.9000_wp, 3.7714_wp, 3.7429_wp, 3.7429_wp, & + 3.6714_wp, 3.7714_wp, 3.8857_wp, 3.9571_wp, 3.9286_wp, & + 3.9286_wp, 3.8857_wp, 4.0143_wp, 4.0000_wp, 5.1143_wp, & + 4.6428_wp, 4.4429_wp, 4.2285_wp, 4.3857_wp, 4.3714_wp, & + 4.3571_wp, 4.3429_wp, 4.2857_wp, 4.3000_wp, 4.2857_wp, & + 4.2714_wp, 4.2571_wp, 4.2571_wp, 4.2428_wp, 4.3143_wp, & + 4.2143_wp, 4.0857_wp, 4.0000_wp, 3.8857_wp, 3.8143_wp, & + 3.7857_wp, 3.7143_wp, 3.7286_wp, 3.7429_wp, 4.0143_wp, & + 3.9857_wp, 3.9857_wp, 4.0714_wp, 4.0000_wp, 4.1000_wp, & + 4.1572_wp, 5.0000_wp, 4.7143_wp, 4.5143_wp, 4.3857_wp, & + 4.3000_wp, 4.3143_wp, 4.3286_wp, 4.3429_wp, 4.2571_wp, & + 2.5857_wp, 2.7857_wp, 3.8428_wp, 3.4714_wp, 3.2286_wp, & + 3.2000_wp, 3.1428_wp, 3.0286_wp, 3.0429_wp, 3.0857_wp, & + 4.1286_wp, 3.9143_wp, 3.7429_wp, 3.6143_wp, 3.7000_wp, & + 3.5857_wp, 3.5428_wp, 3.5000_wp, 4.6428_wp, 4.3286_wp, & + 4.0286_wp, 3.8714_wp, 3.8571_wp, 3.7000_wp, 3.6571_wp, & + 3.6143_wp, 3.5571_wp, 3.5428_wp, 3.5714_wp, 3.6857_wp, & + 3.7286_wp, 3.6857_wp, 3.7714_wp, 3.7000_wp, 3.7572_wp, & + 3.8000_wp, 4.8286_wp, 4.5143_wp, 4.2285_wp, 4.1143_wp, & + 4.0143_wp, 3.9000_wp, 3.7714_wp, 3.7429_wp, 3.7429_wp, & + 3.6714_wp, 3.7714_wp, 3.8857_wp, 3.9571_wp, 3.9286_wp, & + 3.9286_wp, 3.8857_wp, 4.0143_wp, 4.0000_wp, 5.1143_wp, & + 4.6428_wp, 4.4429_wp, 4.2285_wp, 4.3857_wp, 4.3714_wp, & + 4.3571_wp, 4.3429_wp, 4.2857_wp, 4.3000_wp, 4.2857_wp, & + 4.2714_wp, 4.2571_wp, 4.2571_wp, 4.2428_wp, 4.3143_wp, & + 4.2143_wp, 4.0857_wp, 4.0000_wp, 3.8857_wp, 3.8143_wp, & + 3.7857_wp, 3.7143_wp, 3.7286_wp, 3.7429_wp, 4.0143_wp, & + 3.9857_wp, 3.9857_wp, 4.0714_wp, 4.0000_wp, 4.1000_wp, & + 4.1572_wp, 5.0000_wp, 4.7143_wp, 4.5143_wp, 4.3857_wp, & + 4.3000_wp, 4.3143_wp, 4.3286_wp, 4.3429_wp, 4.2571_wp, & + 4.2571_wp, 2.6143_wp, 2.8143_wp, 3.8714_wp, 3.5000_wp, & + 3.2571_wp, 3.2286_wp, 3.1714_wp, 3.0572_wp, 3.0715_wp, & + 3.1143_wp, 4.1572_wp, 3.9428_wp, 3.7714_wp, 3.6429_wp, & + 3.7286_wp, 3.6143_wp, 3.5714_wp, 3.5286_wp, 4.6714_wp, & + 4.3571_wp, 4.0571_wp, 3.9000_wp, 3.8857_wp, 3.7286_wp, & + 3.6857_wp, 3.6429_wp, 3.5857_wp, 3.5714_wp, 3.6000_wp, & + 3.7143_wp, 3.7572_wp, 3.7143_wp, 3.8000_wp, 3.7286_wp, & + 3.7857_wp, 3.8286_wp, 4.8572_wp, 4.5429_wp, 4.2571_wp, & + 4.1429_wp, 4.0429_wp, 3.9286_wp, 3.8000_wp, 3.7714_wp, & + 3.7714_wp, 3.7000_wp, 3.8000_wp, 3.9143_wp, 3.9857_wp, & + 3.9571_wp, 3.9571_wp, 3.9143_wp, 4.0429_wp, 4.0286_wp, & + 5.1429_wp, 4.6714_wp, 4.4714_wp, 4.2571_wp, 4.4143_wp, & + 4.4000_wp, 4.3857_wp, 4.3714_wp, 4.3143_wp, 4.3286_wp, & + 4.3143_wp, 4.3000_wp, 4.2857_wp, 4.2857_wp, 4.2714_wp, & + 4.3429_wp, 4.2428_wp, 4.1143_wp, 4.0286_wp, 3.9143_wp, & + 3.8428_wp, 3.8143_wp, 3.7429_wp, 3.7572_wp, 3.7714_wp, & + 4.0429_wp, 4.0143_wp, 4.0143_wp, 4.1000_wp, 4.0286_wp, & + 4.1286_wp, 4.1857_wp, 5.0286_wp, 4.7429_wp, 4.5429_wp, & + 4.4143_wp, 4.3286_wp, 4.3429_wp, 4.3571_wp, 4.3714_wp, & + 4.2857_wp, 4.2857_wp, 4.3143_wp, 2.6143_wp, 2.8143_wp, & + 3.8714_wp, 3.5000_wp, 3.2571_wp, 3.2286_wp, 3.1714_wp, & + 3.0572_wp, 3.0715_wp, 3.1143_wp, 4.1572_wp, 3.9428_wp, & + 3.7714_wp, 3.6429_wp, 3.7286_wp, 3.6143_wp, 3.5714_wp, & + 3.5286_wp, 4.6714_wp, 4.3571_wp, 4.0571_wp, 3.9000_wp, & + 3.8857_wp, 3.7286_wp, 3.6857_wp, 3.6429_wp, 3.5857_wp, & + 3.5714_wp, 3.6000_wp, 3.7143_wp, 3.7572_wp, 3.7143_wp, & + 3.8000_wp, 3.7286_wp, 3.7857_wp, 3.8286_wp, 4.8572_wp, & + 4.5429_wp, 4.2571_wp, 4.1429_wp, 4.0429_wp, 3.9286_wp, & + 3.8000_wp, 3.7714_wp, 3.7714_wp, 3.7000_wp, 3.8000_wp, & + 3.9143_wp, 3.9857_wp, 3.9571_wp, 3.9571_wp, 3.9143_wp, & + 4.0429_wp, 4.0286_wp, 5.1429_wp, 4.6714_wp, 4.4714_wp, & + 4.2571_wp, 4.4143_wp, 4.4000_wp, 4.3857_wp, 4.3714_wp, & + 4.3143_wp, 4.3286_wp, 4.3143_wp, 4.3000_wp, 4.2857_wp, & + 4.2857_wp, 4.2714_wp, 4.3429_wp, 4.2428_wp, 4.1143_wp, & + 4.0286_wp, 3.9143_wp, 3.8428_wp, 3.8143_wp, 3.7429_wp, & + 3.7572_wp, 3.7714_wp, 4.0429_wp, 4.0143_wp, 4.0143_wp, & + 4.1000_wp, 4.0286_wp, 4.1286_wp, 4.1857_wp, 5.0286_wp, & + 4.7429_wp, 4.5429_wp, 4.4143_wp, 4.3286_wp, 4.3429_wp, & + 4.3571_wp, 4.3714_wp, 4.2857_wp, 4.2857_wp, 4.3143_wp, & + 4.3143_wp, 2.5714_wp, 2.7714_wp, 3.8286_wp, 3.4572_wp, & + 3.2143_wp, 3.1857_wp, 3.1285_wp, 3.0143_wp, 3.0286_wp, & + 3.0715_wp, 4.1143_wp, 3.9000_wp, 3.7286_wp, 3.6000_wp, & + 3.6857_wp, 3.5714_wp, 3.5286_wp, 3.4857_wp, 4.6286_wp, & + 4.3143_wp, 4.0143_wp, 3.8571_wp, 3.8428_wp, 3.6857_wp, & + 3.6429_wp, 3.6000_wp, 3.5428_wp, 3.5286_wp, 3.5571_wp, & + 3.6714_wp, 3.7143_wp, 3.6714_wp, 3.7572_wp, 3.6857_wp, & + 3.7429_wp, 3.7857_wp, 4.8143_wp, 4.5000_wp, 4.2143_wp, & + 4.1000_wp, 4.0000_wp, 3.8857_wp, 3.7572_wp, 3.7286_wp, & + 3.7286_wp, 3.6571_wp, 3.7572_wp, 3.8714_wp, 3.9428_wp, & + 3.9143_wp, 3.9143_wp, 3.8714_wp, 4.0000_wp, 3.9857_wp, & + 5.1000_wp, 4.6286_wp, 4.4286_wp, 4.2143_wp, 4.3714_wp, & + 4.3571_wp, 4.3429_wp, 4.3286_wp, 4.2714_wp, 4.2857_wp, & + 4.2714_wp, 4.2571_wp, 4.2428_wp, 4.2428_wp, 4.2285_wp, & + 4.3000_wp, 4.2000_wp, 4.0714_wp, 3.9857_wp, 3.8714_wp, & + 3.8000_wp, 3.7714_wp, 3.7000_wp, 3.7143_wp, 3.7286_wp, & + 4.0000_wp, 3.9714_wp, 3.9714_wp, 4.0571_wp, 3.9857_wp, & + 4.0857_wp, 4.1429_wp, 4.9857_wp, 4.7000_wp, 4.5000_wp, & + 4.3714_wp, 4.2857_wp, 4.3000_wp, 4.3143_wp, 4.3286_wp, & + 4.2428_wp, 4.2428_wp, 4.2714_wp, 4.2714_wp, 4.2285_wp, & + 2.6000_wp, 2.8000_wp, 3.8571_wp, 3.4857_wp, 3.2429_wp, & + 3.2143_wp, 3.1571_wp, 3.0429_wp, 3.0572_wp, 3.1000_wp, & + 4.1429_wp, 3.9286_wp, 3.7572_wp, 3.6286_wp, 3.7143_wp, & + 3.6000_wp, 3.5571_wp, 3.5143_wp, 4.6571_wp, 4.3429_wp, & + 4.0429_wp, 3.8857_wp, 3.8714_wp, 3.7143_wp, 3.6714_wp, & + 3.6286_wp, 3.5714_wp, 3.5571_wp, 3.5857_wp, 3.7000_wp, & + 3.7429_wp, 3.7000_wp, 3.7857_wp, 3.7143_wp, 3.7714_wp, & + 3.8143_wp, 4.8429_wp, 4.5286_wp, 4.2428_wp, 4.1286_wp, & + 4.0286_wp, 3.9143_wp, 3.7857_wp, 3.7572_wp, 3.7572_wp, & + 3.6857_wp, 3.7857_wp, 3.9000_wp, 3.9714_wp, 3.9428_wp, & + 3.9428_wp, 3.9000_wp, 4.0286_wp, 4.0143_wp, 5.1286_wp, & + 4.6571_wp, 4.4572_wp, 4.2428_wp, 4.4000_wp, 4.3857_wp, & + 4.3714_wp, 4.3571_wp, 4.3000_wp, 4.3143_wp, 4.3000_wp, & + 4.2857_wp, 4.2714_wp, 4.2714_wp, 4.2571_wp, 4.3286_wp, & + 4.2285_wp, 4.1000_wp, 4.0143_wp, 3.9000_wp, 3.8286_wp, & + 3.8000_wp, 3.7286_wp, 3.7429_wp, 3.7572_wp, 4.0286_wp, & + 4.0000_wp, 4.0000_wp, 4.0857_wp, 4.0143_wp, 4.1143_wp, & + 4.1715_wp, 5.0143_wp, 4.7286_wp, 4.5286_wp, 4.4000_wp, & + 4.3143_wp, 4.3286_wp, 4.3429_wp, 4.3571_wp, 4.2714_wp, & + 4.2714_wp, 4.3000_wp, 4.3000_wp, 4.2571_wp, 4.2857_wp, & + 2.6857_wp, 2.8857_wp, 3.9428_wp, 3.5714_wp, 3.3286_wp, & + 3.3000_wp, 3.2429_wp, 3.1285_wp, 3.1428_wp, 3.1857_wp, & + 4.2285_wp, 4.0143_wp, 3.8428_wp, 3.7143_wp, 3.8000_wp, & + 3.6857_wp, 3.6429_wp, 3.6000_wp, 4.7429_wp, 4.4286_wp, & + 4.1286_wp, 3.9714_wp, 3.9571_wp, 3.8000_wp, 3.7572_wp, & + 3.7143_wp, 3.6571_wp, 3.6429_wp, 3.6714_wp, 3.7857_wp, & + 3.8286_wp, 3.7857_wp, 3.8714_wp, 3.8000_wp, 3.8571_wp, & + 3.9000_wp, 4.9285_wp, 4.6143_wp, 4.3286_wp, 4.2143_wp, & + 4.1143_wp, 4.0000_wp, 3.8714_wp, 3.8428_wp, 3.8428_wp, & + 3.7714_wp, 3.8714_wp, 3.9857_wp, 4.0571_wp, 4.0286_wp, & + 4.0286_wp, 3.9857_wp, 4.1143_wp, 4.1000_wp, 5.2143_wp, & + 4.7429_wp, 4.5429_wp, 4.3286_wp, 4.4857_wp, 4.4714_wp, & + 4.4572_wp, 4.4429_wp, 4.3857_wp, 4.4000_wp, 4.3857_wp, & + 4.3714_wp, 4.3571_wp, 4.3571_wp, 4.3429_wp, 4.4143_wp, & + 4.3143_wp, 4.1857_wp, 4.1000_wp, 3.9857_wp, 3.9143_wp, & + 3.8857_wp, 3.8143_wp, 3.8286_wp, 3.8428_wp, 4.1143_wp, & + 4.0857_wp, 4.0857_wp, 4.1715_wp, 4.1000_wp, 4.2000_wp, & + 4.2571_wp, 5.1000_wp, 4.8143_wp, 4.6143_wp, 4.4857_wp, & + 4.4000_wp, 4.4143_wp, 4.4286_wp, 4.4429_wp, 4.3571_wp, & + 4.3571_wp, 4.3857_wp, 4.3857_wp, 4.3429_wp, 4.3714_wp, & + 4.4572_wp, 2.7143_wp, 2.9143_wp, 3.9714_wp, 3.6000_wp, & + 3.3572_wp, 3.3286_wp, 3.2714_wp, 3.1571_wp, 3.1714_wp, & + 3.2143_wp, 4.2571_wp, 4.0429_wp, 3.8714_wp, 3.7429_wp, & + 3.8286_wp, 3.7143_wp, 3.6714_wp, 3.6286_wp, 4.7714_wp, & + 4.4572_wp, 4.1572_wp, 4.0000_wp, 3.9857_wp, 3.8286_wp, & + 3.7857_wp, 3.7429_wp, 3.6857_wp, 3.6714_wp, 3.7000_wp, & + 3.8143_wp, 3.8571_wp, 3.8143_wp, 3.9000_wp, 3.8286_wp, & + 3.8857_wp, 3.9286_wp, 4.9571_wp, 4.6428_wp, 4.3571_wp, & + 4.2428_wp, 4.1429_wp, 4.0286_wp, 3.9000_wp, 3.8714_wp, & + 3.8714_wp, 3.8000_wp, 3.9000_wp, 4.0143_wp, 4.0857_wp, & + 4.0571_wp, 4.0571_wp, 4.0143_wp, 4.1429_wp, 4.1286_wp, & + 5.2429_wp, 4.7714_wp, 4.5714_wp, 4.3571_wp, 4.5143_wp, & + 4.5000_wp, 4.4857_wp, 4.4714_wp, 4.4143_wp, 4.4286_wp, & + 4.4143_wp, 4.4000_wp, 4.3857_wp, 4.3857_wp, 4.3714_wp, & + 4.4429_wp, 4.3429_wp, 4.2143_wp, 4.1286_wp, 4.0143_wp, & + 3.9428_wp, 3.9143_wp, 3.8428_wp, 3.8571_wp, 3.8714_wp, & + 4.1429_wp, 4.1143_wp, 4.1143_wp, 4.2000_wp, 4.1286_wp, & + 4.2285_wp, 4.2857_wp, 5.1286_wp, 4.8429_wp, 4.6428_wp, & + 4.5143_wp, 4.4286_wp, 4.4429_wp, 4.4572_wp, 4.4714_wp, & + 4.3857_wp, 4.3857_wp, 4.4143_wp, 4.4143_wp, 4.3714_wp, & + 4.4000_wp, 4.4857_wp, 4.5143_wp, 2.5286_wp, 2.7286_wp, & + 3.7857_wp, 3.4143_wp, 3.1714_wp, 3.1428_wp, 3.0857_wp, & + 2.9714_wp, 2.9857_wp, 3.0286_wp, 4.0714_wp, 3.8571_wp, & + 3.6857_wp, 3.5571_wp, 3.6429_wp, 3.5286_wp, 3.4857_wp, & + 3.4429_wp, 4.5857_wp, 4.2714_wp, 3.9714_wp, 3.8143_wp, & + 3.8000_wp, 3.6429_wp, 3.6000_wp, 3.5571_wp, 3.5000_wp, & + 3.4857_wp, 3.5143_wp, 3.6286_wp, 3.6714_wp, 3.6286_wp, & + 3.7143_wp, 3.6429_wp, 3.7000_wp, 3.7429_wp, 4.7714_wp, & + 4.4572_wp, 4.1715_wp, 4.0571_wp, 3.9571_wp, 3.8428_wp, & + 3.7143_wp, 3.6857_wp, 3.6857_wp, 3.6143_wp, 3.7143_wp, & + 3.8286_wp, 3.9000_wp, 3.8714_wp, 3.8714_wp, 3.8286_wp, & + 3.9571_wp, 3.9428_wp, 5.0571_wp, 4.5857_wp, 4.3857_wp, & + 4.1715_wp, 4.3286_wp, 4.3143_wp, 4.3000_wp, 4.2857_wp, & + 4.2285_wp, 4.2428_wp, 4.2285_wp, 4.2143_wp, 4.2000_wp, & + 4.2000_wp, 4.1857_wp, 4.2571_wp, 4.1572_wp, 4.0286_wp, & + 3.9428_wp, 3.8286_wp, 3.7572_wp, 3.7286_wp, 3.6571_wp, & + 3.6714_wp, 3.6857_wp, 3.9571_wp, 3.9286_wp, 3.9286_wp, & + 4.0143_wp, 3.9428_wp, 4.0429_wp, 4.1000_wp, 4.9428_wp, & + 4.6571_wp, 4.4572_wp, 4.3286_wp, 4.2428_wp, 4.2571_wp, & + 4.2714_wp, 4.2857_wp, 4.2000_wp, 4.2000_wp, 4.2285_wp, & + 4.2285_wp, 4.1857_wp, 4.2143_wp, 4.3000_wp, 4.3286_wp, & + 4.1429_wp & + ] + + +contains + + + !> Get van-der-Waals radius for a single element symbol +elemental function get_vdw_rad_sym(sym) result(rad) + + !> Element symbol + character(len=*), intent(in) :: sym + + !> Van-der-Waals radius + real(wp) :: rad + + rad = get_vdw_rad(to_number(sym)) + +end function get_vdw_rad_sym + + +!> Get van-der-Waals radius for a pair of element symbols +elemental function get_vdw_rad_pair_sym(sym1, sym2) result(rad) + + !> Element symbol + character(len=*), intent(in) :: sym1 + + !> Element symbol + character(len=*), intent(in) :: sym2 + + !> Van-der-Waals radius + real(wp) :: rad + + rad = get_vdw_rad(to_number(sym1), to_number(sym2)) + +end function get_vdw_rad_pair_sym + + +!> Get van-der-Waals radius for a single atomic number +elemental function get_vdw_rad_num(num) result(rad) + + !> Atomic number + integer, intent(in) :: num + + !> Van-der-Waals radius + real(wp) :: rad + + rad = 0.5_wp * get_vdw_rad(num, num) + +end function get_vdw_rad_num + + +!> Get van-der-Waals radius for a pair of atomic numbers +elemental function get_vdw_rad_pair_num(num1, num2) result(rad) + + !> Atomic number + integer, intent(in) :: num1 + + !> Atomic number + integer, intent(in) :: num2 + + !> Van-der-Waals radius + real(wp) :: rad + + ! Magic number from ratio of van-der-Waals and covalent radii over element 1-94. + ! + ! L. Trombach, S. Ehlert, S. Grimme, P Schwerdtfeger, J.-M. Mewes, PCCP, 2019. + ! DOI: 10.1039/c9cp02455g + real(wp), parameter :: cov2vdw = 1.45_wp + + if (num1 > 0 .and. num1 <= max_elem .and. num2 > 0 .and. num2 <= max_elem) then + if (num1 > num2) then + rad = vdwrad(num2 + num1*(num1-1)/2) + else + rad = vdwrad(num1 + num2*(num2-1)/2) + end if + else + ! best estimate + if (num1 > 0 .and. num1 <= max_elem) then + rad = 0.5_wp * vdwrad(num1*(num1+1)/2) + cov2vdw * get_covalent_rad(num2) + else if (num2 > 0 .and. num2 <= max_elem) then + rad = 0.5_wp * vdwrad(num2*(num2+1)/2) + cov2vdw * get_covalent_rad(num1) + else + rad = cov2vdw * (get_covalent_rad(num1) + get_covalent_rad(num2)) + end if + end if + +end function get_vdw_rad_pair_num + + +end module mctc_data_vdwrad + + diff --git a/src/mctc/meson.build b/src/mctc/meson.build index d1e35b02..03cd35b6 100644 --- a/src/mctc/meson.build +++ b/src/mctc/meson.build @@ -12,11 +12,16 @@ # See the License for the specific language governing permissions and # limitations under the License. +subdir('data') subdir('env') subdir('io') +subdir('ncoord') srcs += files( + 'cutoff.f90', + 'data.f90', 'env.f90', 'io.f90', + 'ncoord.f90', 'version.F90', ) diff --git a/src/mctc/ncoord.f90 b/src/mctc/ncoord.f90 new file mode 100644 index 00000000..b45176aa --- /dev/null +++ b/src/mctc/ncoord.f90 @@ -0,0 +1,102 @@ +! 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/ncoord +!> Contains the implementation for the coordination number evaluators. + +!> @file mctc/ncoord.f90 +!> Reexports the coordination number evaluation modules. + +!> Proxy module to expose coordination number containers +module mctc_ncoord + use mctc_env, only : wp + use mctc_io, only : structure_type + use mctc_ncoord_type, only : ncoord_type + use mctc_ncoord_dexp, only : dexp_ncoord_type, new_dexp_ncoord + use mctc_ncoord_exp, only : exp_ncoord_type, new_exp_ncoord + use mctc_ncoord_erf, only : erf_ncoord_type, new_erf_ncoord + use mctc_ncoord_erf_en, only : erf_en_ncoord_type, new_erf_en_ncoord + use mctc_ncoord_erf_dftd4, only : erf_dftd4_ncoord_type, new_erf_dftd4_ncoord + implicit none + private + + public :: ncoord_type, new_ncoord + +contains + +!> Create a new generic coordination number container +subroutine new_ncoord(self, mol, cn_type, kcn, cutoff, rcov, en, cut, norm_exp) + !> Instance of the coordination number container + class(ncoord_type), allocatable, intent(out) :: self + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Coordination number type + character(len=*), intent(in) :: cn_type + !> Optional steepness of counting function + real(wp), intent(in), optional :: kcn + !> Optional real space cutoff + real(wp), intent(in), optional :: cutoff + !> Optional set of covalent radii to be used in CN + real(wp), intent(in), optional :: rcov(:) + !> Optional set of electronegativity to be use din CN + real(wp), intent(in), optional :: en(:) + !> Optional cutoff for the maximum coordination number + real(wp), intent(in), optional :: cut + !> Optional exponent of the distance normalization + real(wp), intent(in), optional :: norm_exp + + select case(cn_type) + case("exp") + block + type(exp_ncoord_type), allocatable :: tmp + allocate(tmp) + call new_exp_ncoord(tmp, mol, kcn=kcn, cutoff=cutoff, rcov=rcov, cut=cut) + call move_alloc(tmp, self) + end block + case("dexp") + block + type(dexp_ncoord_type), allocatable :: tmp + allocate(tmp) + call new_dexp_ncoord(tmp, mol, cutoff=cutoff, rcov=rcov, cut=cut) + call move_alloc(tmp, self) + end block + case("erf") + block + type(erf_ncoord_type), allocatable :: tmp + allocate(tmp) + call new_erf_ncoord(tmp, mol, kcn=kcn, cutoff=cutoff, & + & rcov=rcov, cut=cut, norm_exp=norm_exp) + call move_alloc(tmp, self) + end block + case("erf_en") + block + type(erf_en_ncoord_type), allocatable :: tmp + allocate(tmp) + call new_erf_en_ncoord(tmp, mol, kcn=kcn, cutoff=cutoff, & + & rcov=rcov, en=en, cut=cut, norm_exp=norm_exp) + call move_alloc(tmp, self) + end block + case("dftd4") + block + type(erf_dftd4_ncoord_type), allocatable :: tmp + allocate(tmp) + call new_erf_dftd4_ncoord(tmp, mol, kcn=kcn, cutoff=cutoff, & + & rcov=rcov, en=en, cut=cut, norm_exp=norm_exp) + call move_alloc(tmp, self) + end block + end select + +end subroutine new_ncoord + +end module mctc_ncoord diff --git a/src/mctc/ncoord/CMakeLists.txt b/src/mctc/ncoord/CMakeLists.txt new file mode 100644 index 00000000..356c1ba1 --- /dev/null +++ b/src/mctc/ncoord/CMakeLists.txt @@ -0,0 +1,27 @@ +# 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. + +add_subdirectory("erf") + +set(dir "${CMAKE_CURRENT_SOURCE_DIR}") + +list( + APPEND srcs + "${dir}/dexp.f90" + "${dir}/exp.f90" + "${dir}/erf.f90" + "${dir}/type.f90" +) + +set(srcs "${srcs}" PARENT_SCOPE) diff --git a/src/mctc/ncoord/dexp.f90 b/src/mctc/ncoord/dexp.f90 new file mode 100644 index 00000000..80cccefc --- /dev/null +++ b/src/mctc/ncoord/dexp.f90 @@ -0,0 +1,162 @@ +! 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/ncoord/dexp.f90 +!> Provides a coordination number implementation with double exponential counting function + +!> Coordination number implementation using a double exponential counting function as in GFN2-xTB +module mctc_ncoord_dexp + use mctc_env, only : wp + use mctc_io, only : structure_type + use mctc_data_covrad, only : get_covalent_rad + use mctc_ncoord_type, only : ncoord_type + implicit none + private + + public :: new_dexp_ncoord + + + !> Coordination number evaluator + type, public, extends(ncoord_type) :: dexp_ncoord_type + real(wp), allocatable :: rcov(:) + contains + !> Evaluates the dexp counting function + procedure :: ncoord_count + !> Evaluates the derivative of the dexp counting function + procedure :: ncoord_dcount + end type dexp_ncoord_type + + !> Steepness of first counting function + real(wp),parameter :: ka = 10.0_wp + !> Steepness of second counting function + real(wp),parameter :: kb = 20.0_wp + !> Offset of the second counting function + real(wp),parameter :: r_shift = 2.0_wp + + real(wp), parameter :: default_cutoff = 25.0_wp + + +contains + + +subroutine new_dexp_ncoord(self, mol, cutoff, rcov, cut) + !> Coordination number container + type(dexp_ncoord_type), intent(out) :: self + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Real space cutoff + real(wp), intent(in), optional :: cutoff + !> Covalent radii + real(wp), intent(in), optional :: rcov(:) + !> Cutoff for the maximum coordination number + real(wp), intent(in), optional :: cut + + if (present(cutoff)) then + self%cutoff = cutoff + else + self%cutoff = default_cutoff + end if + + allocate(self%rcov(mol%nid)) + if (present(rcov)) then + self%rcov(:) = rcov + else + self%rcov(:) = get_covalent_rad(mol%num) + end if + + self%directed_factor = 1.0_wp + + if (present(cut)) then + self%cut = cut + else + ! Negative value deactivates the cutoff + self%cut = -1.0_wp + end if + +end subroutine new_dexp_ncoord + + +!> Double-exponential counting function for coordination number contributions. +elemental function ncoord_count(self, izp, jzp, r) result(count) + !> Coordination number container + class(dexp_ncoord_type), intent(in) :: self + !> Atom i index + integer, intent(in) :: izp + !> Atom j index + integer, intent(in) :: jzp + !> Current distance. + real(wp), intent(in) :: r + + real(wp) :: rc, count + + rc = self%rcov(izp) + self%rcov(jzp) + + count = exp_count(ka, r, rc) * exp_count(kb, r, rc + r_shift) + +end function ncoord_count + +!> Derivative of the double-exponential counting function w.r.t. the distance. +elemental function ncoord_dcount(self, izp, jzp, r) result(count) + !> Coordination number container + class(dexp_ncoord_type), intent(in) :: self + !> Atom i index + integer, intent(in) :: izp + !> Atom j index + integer, intent(in) :: jzp + !> Current distance. + real(wp), intent(in) :: r + + real(wp) :: rc, count + + rc = self%rcov(izp) + self%rcov(jzp) + + count = (exp_dcount(ka, r, rc) * exp_count(kb, r, rc + r_shift) & + & + exp_count(ka, r, rc) * exp_dcount(kb, r, rc + r_shift)) + +end function ncoord_dcount + + +!> Mono-exponential counting function for coordination number contributions. +elemental function exp_count(k, r, r0) result(count) + !> Steepness of the counting function. + real(wp), intent(in) :: k + !> Current distance. + real(wp), intent(in) :: r + !> Cutoff radius. + real(wp), intent(in) :: r0 + + real(wp) :: count + + count = 1.0_wp/(1.0_wp+exp(-k*(r0/r-1.0_wp))) + +end function exp_count + +!> Derivative of the mono-exponential counting function w.r.t. the distance. +elemental function exp_dcount(k, r, r0) result(count) + !> Steepness of the counting function. + real(wp), intent(in) :: k + !> Current distance. + real(wp), intent(in) :: r + !> Cutoff radius. + real(wp), intent(in) :: r0 + + real(wp) :: count + real(wp) :: expterm + + expterm = exp(-k*(r0/r-1._wp)) + count = (-k*r0*expterm)/(r**2*((expterm+1._wp)**2)) + +end function exp_dcount + +end module mctc_ncoord_dexp diff --git a/src/mctc/ncoord/erf.f90 b/src/mctc/ncoord/erf.f90 new file mode 100644 index 00000000..dd123f5e --- /dev/null +++ b/src/mctc/ncoord/erf.f90 @@ -0,0 +1,150 @@ +! 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/ncoord/erf +!> Provides variations (EN) of the error function based coordination number + +!> @file mctc/ncoord/erf.f90 +!> Provides error function based coordination number + +!> Coordination number implementation with single error function +module mctc_ncoord_erf + use mctc_env, only : wp + use mctc_io, only : structure_type + use mctc_io_constants, only : pi + use mctc_data_covrad, only : get_covalent_rad + use mctc_ncoord_type, only : ncoord_type + implicit none + private + + public :: new_erf_ncoord + + !> Coordination number evaluator + type, public, extends(ncoord_type) :: erf_ncoord_type + real(wp), allocatable :: rcov(:) + !> Exponent of the distance normalization + real(wp) :: norm_exp + contains + !> Evaluates the error counting function + procedure :: ncoord_count + !> Evaluates the derivative of the error counting function + procedure :: ncoord_dcount + end type erf_ncoord_type + + !> Steepness of counting function (CEH) + real(wp), parameter :: default_kcn = 3.15_wp + + !> Exponent of distance normalization + real(wp), parameter :: default_norm_exp = 1.0_wp + + real(wp), parameter :: default_cutoff = 25.0_wp + +contains + + + subroutine new_erf_ncoord(self, mol, kcn, cutoff, rcov, cut, norm_exp) + !> Coordination number container + type(erf_ncoord_type), intent(out) :: self + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Steepness of counting function + real(wp), optional :: kcn + !> Real space cutoff + real(wp), intent(in), optional :: cutoff + !> Covalent radii + real(wp), intent(in), optional :: rcov(:) + !> Cutoff for the maximum coordination number + real(wp), intent(in), optional :: cut + !> Exponent of the distance normalization + real(wp), intent(in), optional :: norm_exp + + if(present(kcn)) then + self%kcn = kcn + else + self%kcn = default_kcn + end if + + if (present(cutoff)) then + self%cutoff = cutoff + else + self%cutoff = default_cutoff + end if + + allocate(self%rcov(mol%nid)) + if (present(rcov)) then + self%rcov(:) = rcov + else + self%rcov(:) = get_covalent_rad(mol%num) + end if + + self%directed_factor = 1.0_wp + + if (present(cut)) then + self%cut = cut + else + ! Negative value deactivates the cutoff + self%cut = -1.0_wp + end if + + if (present(norm_exp)) then + self%norm_exp = norm_exp + else + self%norm_exp = default_norm_exp + end if + + end subroutine new_erf_ncoord + + + !> Error counting function for coordination number contributions. + elemental function ncoord_count(self, izp, jzp, r) result(count) + !> Coordination number container + class(erf_ncoord_type), intent(in) :: self + !> Atom i index + integer, intent(in) :: izp + !> Atom j index + integer, intent(in) :: jzp + !> Current distance. + real(wp), intent(in) :: r + + real(wp) :: rc, count + + rc = (self%rcov(izp) + self%rcov(jzp)) + + count = 0.5_wp * (1.0_wp + erf(-self%kcn*(r-rc)/rc**self%norm_exp)) + + end function ncoord_count + + !> Derivative of the error counting function w.r.t. the distance. + elemental function ncoord_dcount(self, izp, jzp, r) result(count) + !> Coordination number container + class(erf_ncoord_type), intent(in) :: self + !> Atom i index + integer, intent(in) :: izp + !> Atom j index + integer, intent(in) :: jzp + !> Current distance. + real(wp), intent(in) :: r + + real(wp), parameter :: sqrtpi = sqrt(pi) + real(wp) :: rc, exponent, expterm, count + + rc = self%rcov(izp) + self%rcov(jzp) + + exponent = self%kcn*(r-rc)/rc**self%norm_exp + expterm = exp(-exponent**2) + count = -(self%kcn*expterm)/(sqrtpi*rc**self%norm_exp) + + end function ncoord_dcount + +end module mctc_ncoord_erf diff --git a/src/mctc/ncoord/erf/CMakeLists.txt b/src/mctc/ncoord/erf/CMakeLists.txt new file mode 100644 index 00000000..948715c2 --- /dev/null +++ b/src/mctc/ncoord/erf/CMakeLists.txt @@ -0,0 +1,23 @@ +# 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}/dftd4.f90" + "${dir}/en.f90" +) + +set(srcs "${srcs}" PARENT_SCOPE) diff --git a/src/mctc/ncoord/erf/dftd4.f90 b/src/mctc/ncoord/erf/dftd4.f90 new file mode 100644 index 00000000..fd761432 --- /dev/null +++ b/src/mctc/ncoord/erf/dftd4.f90 @@ -0,0 +1,135 @@ +! 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/ncoord/erf_en.f90 +!> Provides an implementation for the CN as used dftd4 + +!> Coordination number implementation with single error function and EN-weighting for dftd4 +module mctc_ncoord_erf_dftd4 + use mctc_env, only : wp + use mctc_io, only : structure_type + use mctc_data_covrad, only : get_covalent_rad + use mctc_data_paulingen, only : get_pauling_en + use mctc_ncoord_erf, only : erf_ncoord_type + implicit none + private + + public :: new_erf_dftd4_ncoord + + !> Coordination number evaluator + type, public, extends(erf_ncoord_type) :: erf_dftd4_ncoord_type + real(wp), allocatable :: en(:) + contains + !> Evaluates pairwise electronegativity factor + procedure :: get_en_factor + end type erf_dftd4_ncoord_type + + !> Steepness of counting function + real(wp), parameter :: default_kcn = 7.5_wp + + !> Exponent of distance normalization + real(wp), parameter :: default_norm_exp = 1.0_wp + + real(wp), parameter :: default_cutoff = 25.0_wp + + !> Parameter for electronegativity scaling + real(wp), parameter :: k4 = 4.10451_wp + + !> Parameter for electronegativity scaling + real(wp), parameter :: k5 = 19.08857_wp + + !> Parameter for electronegativity scaling + real(wp), parameter :: k6 = 2*11.28174_wp**2 + +contains + + + subroutine new_erf_dftd4_ncoord(self, mol, kcn, cutoff, rcov, en, cut, norm_exp) + !> Coordination number container + type(erf_dftd4_ncoord_type), intent(out) :: self + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Steepness of counting function + real(wp), optional :: kcn + !> Real space cutoff + real(wp), intent(in), optional :: cutoff + !> Covalent radii + real(wp), intent(in), optional :: rcov(:) + !> Electronegativity + real(wp), intent(in), optional :: en(:) + !> Cutoff for the maximum coordination number + real(wp), intent(in), optional :: cut + !> Exponent of the distance normalization + real(wp), intent(in), optional :: norm_exp + + if(present(kcn)) then + self%kcn = kcn + else + self%kcn = default_kcn + end if + + if (present(cutoff)) then + self%cutoff = cutoff + else + self%cutoff = default_cutoff + end if + + allocate(self%rcov(mol%nid)) + if (present(rcov)) then + self%rcov(:) = rcov + else + self%rcov(:) = get_covalent_rad(mol%num) + end if + + allocate(self%en(mol%nid)) + if (present(en)) then + self%en(:) = en + else + self%en(:) = get_pauling_en(mol%num) + end if + + self%directed_factor = 1.0_wp + + if (present(cut)) then + self%cut = cut + else + ! Negative value deactivates the cutoff + self%cut = -1.0_wp + end if + + if (present(norm_exp)) then + self%norm_exp = norm_exp + else + self%norm_exp = default_norm_exp + end if + + end subroutine new_erf_dftd4_ncoord + + + !> Evaluates pairwise electronegativity factor if non applies + elemental function get_en_factor(self, izp, jzp) result(en_factor) + !> Coordination number container + class(erf_dftd4_ncoord_type), intent(in) :: self + !> Atom i index + integer, intent(in) :: izp + !> Atom j index + integer, intent(in) :: jzp + + real(wp) :: en_factor + + en_factor = k4*exp(-(abs(self%en(izp)-self%en(jzp)) + k5)**2/k6) + + end function get_en_factor + +end module mctc_ncoord_erf_dftd4 diff --git a/src/mctc/ncoord/erf/en.f90 b/src/mctc/ncoord/erf/en.f90 new file mode 100644 index 00000000..4a41c77f --- /dev/null +++ b/src/mctc/ncoord/erf/en.f90 @@ -0,0 +1,128 @@ +! 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/ncoord/erf/en.f90 +!> Provides an implementation for the electronegativity-weighted CN + +!> Coordination number implementation with single error function and EN-weighting. +module mctc_ncoord_erf_en + use mctc_env, only : wp + use mctc_io, only : structure_type + use mctc_data_covrad, only : get_covalent_rad + use mctc_data_paulingen, only : get_pauling_en + use mctc_ncoord_erf, only : erf_ncoord_type + implicit none + private + + public :: new_erf_en_ncoord + + !> Coordination number evaluator + type, public, extends(erf_ncoord_type) :: erf_en_ncoord_type + real(wp), allocatable :: en(:) + contains + !> Evaluates pairwise electronegativity factor + procedure :: get_en_factor + end type erf_en_ncoord_type + + !> Steepness of counting function (CEH) + real(wp), parameter :: default_kcn = 2.65_wp + + !> Exponent of distance normalization + real(wp), parameter :: default_norm_exp = 1.0_wp + + real(wp), parameter :: default_cutoff = 25.0_wp + +contains + + + subroutine new_erf_en_ncoord(self, mol, kcn, cutoff, rcov, en, cut, norm_exp) + !> Coordination number container + type(erf_en_ncoord_type), intent(out) :: self + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Steepness of counting function + real(wp), optional :: kcn + !> Real space cutoff + real(wp), intent(in), optional :: cutoff + !> Covalent radii + real(wp), intent(in), optional :: rcov(:) + !> Electronegativity (normalized to F) + real(wp), intent(in), optional :: en(:) + !> Cutoff for the maximum coordination number + real(wp), intent(in), optional :: cut + !> Exponent of the distance normalization + real(wp), intent(in), optional :: norm_exp + + if(present(kcn)) then + self%kcn = kcn + else + self%kcn = default_kcn + end if + + if (present(cutoff)) then + self%cutoff = cutoff + else + self%cutoff = default_cutoff + end if + + allocate(self%rcov(mol%nid)) + if (present(rcov)) then + self%rcov(:) = rcov + else + self%rcov(:) = get_covalent_rad(mol%num) + end if + + allocate(self%en(mol%nid)) + if (present(en)) then + self%en(:) = en + else + self%en(:) = get_pauling_en(mol%num) + end if + + ! CN is directed due to the EN contribution + ! i.e. added to higher EN and removed from lower EN species + self%directed_factor = -1.0_wp + + if (present(cut)) then + self%cut = cut + else + ! Negative value deactivates the cutoff + self%cut = -1.0_wp + end if + + if (present(norm_exp)) then + self%norm_exp = norm_exp + else + self%norm_exp = default_norm_exp + end if + + end subroutine new_erf_en_ncoord + + + !> Evaluates pairwise electronegativity factor if non applies + elemental function get_en_factor(self, izp, jzp) result(en_factor) + !> Coordination number container + class(erf_en_ncoord_type), intent(in) :: self + !> Atom i index + integer, intent(in) :: izp + !> Atom j index + integer, intent(in) :: jzp + + real(wp) :: en_factor + + en_factor = self%en(jzp) - self%en(izp) + + end function get_en_factor + +end module mctc_ncoord_erf_en diff --git a/src/mctc/ncoord/erf/meson.build b/src/mctc/ncoord/erf/meson.build new file mode 100644 index 00000000..7c30c5f7 --- /dev/null +++ b/src/mctc/ncoord/erf/meson.build @@ -0,0 +1,18 @@ +# 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. + +srcs += files( + 'dftd4.f90', + 'en.f90', +) diff --git a/src/mctc/ncoord/exp.f90 b/src/mctc/ncoord/exp.f90 new file mode 100644 index 00000000..c1848f06 --- /dev/null +++ b/src/mctc/ncoord/exp.f90 @@ -0,0 +1,131 @@ +! 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/ncoord/exp.f90 +!> Provides a coordination number implementation with exponential counting function + +!> Coordination number implementation using an exponential counting function as in dftd3. +module mctc_ncoord_exp + use mctc_env, only : wp + use mctc_io, only : structure_type + use mctc_data_covrad, only : get_covalent_rad + use mctc_ncoord_type, only : ncoord_type + implicit none + private + + public :: new_exp_ncoord + + !> Coordination number evaluator + type, public, extends(ncoord_type) :: exp_ncoord_type + real(wp), allocatable :: rcov(:) + contains + !> Evaluates the exponential counting function + procedure :: ncoord_count + !> Evaluates the derivative of the exponential counting function + procedure :: ncoord_dcount + end type exp_ncoord_type + + !> Steepness of counting function + real(wp), parameter :: default_kcn = 16.0_wp + + real(wp), parameter :: default_cutoff = 25.0_wp + +contains + + + subroutine new_exp_ncoord(self, mol, kcn, cutoff, rcov, cut) + !> Coordination number container + type(exp_ncoord_type), intent(out) :: self + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Steepness of counting function + real(wp), optional :: kcn + !> Real space cutoff + real(wp), intent(in), optional :: cutoff + !> Covalent radii + real(wp), intent(in), optional :: rcov(:) + !> Cutoff for the maximum coordination number + real(wp), intent(in), optional :: cut + + if(present(kcn)) then + self%kcn = kcn + else + self%kcn = default_kcn + end if + + if (present(cutoff)) then + self%cutoff = cutoff + else + self%cutoff = default_cutoff + end if + + allocate(self%rcov(mol%nid)) + if (present(rcov)) then + self%rcov(:) = rcov + else + self%rcov(:) = get_covalent_rad(mol%num) + end if + + self%directed_factor = 1.0_wp + + if (present(cut)) then + self%cut = cut + else + ! Negative value deactivates the cutoff + self%cut = -1.0_wp + end if + + end subroutine new_exp_ncoord + + + !> Exponential counting function for coordination number contributions. + elemental function ncoord_count(self, izp, jzp, r) result(count) + !> Coordination number container + class(exp_ncoord_type), intent(in) :: self + !> Atom i index + integer, intent(in) :: izp + !> Atom j index + integer, intent(in) :: jzp + !> Current distance. + real(wp), intent(in) :: r + + real(wp) :: rc, count + + rc = self%rcov(izp) + self%rcov(jzp) + + count =1.0_wp/(1.0_wp+exp(-self%kcn*(rc/r-1.0_wp))) + + end function ncoord_count + + !> Derivative of the exponential counting function w.r.t. the distance. + elemental function ncoord_dcount(self, izp, jzp, r) result(count) + !> Coordination number container + class(exp_ncoord_type), intent(in) :: self + !> Atom i index + integer, intent(in) :: izp + !> Atom j index + integer, intent(in) :: jzp + !> Current distance. + real(wp), intent(in) :: r + + real(wp) :: rc, expterm, count + + rc = self%rcov(izp) + self%rcov(jzp) + + expterm = exp(-self%kcn*(rc/r-1._wp)) + count = (-self%kcn*rc*expterm)/(r**2*((expterm+1._wp)**2)) + + end function ncoord_dcount + +end module mctc_ncoord_exp diff --git a/src/mctc/ncoord/meson.build b/src/mctc/ncoord/meson.build new file mode 100644 index 00000000..895008ae --- /dev/null +++ b/src/mctc/ncoord/meson.build @@ -0,0 +1,22 @@ +# 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. + +subdir('erf') + +srcs += files( + 'dexp.f90', + 'exp.f90', + 'erf.f90', + 'type.f90', +) diff --git a/src/mctc/ncoord/type.f90 b/src/mctc/ncoord/type.f90 new file mode 100644 index 00000000..72bb62e2 --- /dev/null +++ b/src/mctc/ncoord/type.f90 @@ -0,0 +1,386 @@ +! 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/ncoord/type.f90 +!> Provides a coordination number evalulator base class + +!> Declaration of base class for coordination number evalulations +module mctc_ncoord_type + use mctc_env, only : wp + use mctc_io, only : structure_type + use mctc_cutoff, only : get_lattice_points + + implicit none + private + + !> Abstract base class for coordination number evaluator + type, public, abstract :: ncoord_type + !> Radial cutoff for the coordination number + real(wp) :: cutoff + !> Steepness of counting function + real(wp) :: kcn + !> Factor determining whether the CN is evaluated with direction + !> if +1 the CN contribution is added equally to both partners + !> if -1 (i.e. with the EN-dep.) it is added to one and subtracted from the other + real(wp) :: directed_factor + !> Cutoff for the maximum coordination number (negative value, no cutoff) + real(wp) :: cut = -1.0_wp + contains + !> Obtains lattice information and calls get_coordination number + procedure :: get_cn + !> Decides whether the energy or gradient is calculated + procedure :: get_coordination_number + !> Evaluates the CN from the specific counting function + procedure :: ncoord + !> Evaluates derivative of the CN from the specific counting function + procedure :: ncoord_d + !> Evaluates pairwise electronegativity factor + procedure :: get_en_factor + !> Add CN derivative of an arbitrary function + procedure :: add_coordination_number_derivs + !> Evaluates the counting function (exp, dexp, erf, ...) + procedure(ncoord_count), deferred :: ncoord_count + !> Evaluates the derivative of the counting function (exp, dexp, erf, ...) + procedure(ncoord_dcount), deferred :: ncoord_dcount + end type ncoord_type + + abstract interface + + !> Abstract counting function + elemental function ncoord_count(self, izp, jzp, r) result(count) + import :: ncoord_type, wp + !> Instance of coordination number container + class(ncoord_type), intent(in) :: self + !> Atom i index + integer, intent(in) :: izp + !> Atom j index + integer, intent(in) :: jzp + !> Current distance. + real(wp), intent(in) :: r + + real(wp) :: count + end function ncoord_count + + !> Abstract derivative of the counting function w.r.t. the distance. + elemental function ncoord_dcount(self, izp, jzp, r) result(count) + import :: ncoord_type, wp + !> Instance of coordination number container + class(ncoord_type), intent(in) :: self + !> Atom i index + integer, intent(in) :: izp + !> Atom j index + integer, intent(in) :: jzp + !> Current distance. + real(wp), intent(in) :: r + + real(wp) :: count + end function ncoord_dcount + + end interface + +contains + + !> Wrapper for CN using the CN cutoff for the lattice + subroutine get_cn(self, mol, cn, dcndr, dcndL) + !> Coordination number container + class(ncoord_type), intent(in) :: self + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Error function coordination number. + real(wp), intent(out) :: cn(:) + !> Derivative of the CN with respect to the Cartesian coordinates. + real(wp), intent(out), optional :: dcndr(:, :, :) + !> Derivative of the CN with respect to strain deformations. + real(wp), intent(out), optional :: dcndL(:, :, :) + + real(wp), allocatable :: lattr(:, :) + + call get_lattice_points(mol%periodic, mol%lattice, self%cutoff, lattr) + call get_coordination_number(self, mol, lattr, cn, dcndr, dcndL) + end subroutine get_cn + + !> Geometric fractional coordination number + subroutine get_coordination_number(self, mol, trans, cn, dcndr, dcndL) + + !> Coordination number container + class(ncoord_type), intent(in) :: self + + !> Molecular structure data + type(structure_type), intent(in) :: mol + + !> Lattice points + real(wp), intent(in) :: trans(:, :) + + !> Error function coordination number. + real(wp), intent(out) :: cn(:) + + !> Derivative of the CN with respect to the Cartesian coordinates. + real(wp), intent(out), optional :: dcndr(:, :, :) + + !> Derivative of the CN with respect to strain deformations. + real(wp), intent(out), optional :: dcndL(:, :, :) + + if (present(dcndr) .and. present(dcndL)) then + call ncoord_d(self, mol, trans, cn, dcndr, dcndL) + else + call ncoord(self, mol, trans, cn) + end if + + if (self%cut > 0.0_wp) then + call cut_coordination_number(self%cut, cn, dcndr, dcndL) + end if + + end subroutine get_coordination_number + + + subroutine ncoord(self, mol, trans, cn) + !> Coordination number container + class(ncoord_type), intent(in) :: self + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Lattice points + real(wp), intent(in) :: trans(:, :) + !> Error function coordination number. + real(wp), intent(out) :: cn(:) + + integer :: iat, jat, izp, jzp, itr + real(wp) :: r2, r1, rij(3), countf, cutoff2, den + + cn(:) = 0.0_wp + cutoff2 = self%cutoff**2 + + !$omp parallel do schedule(runtime) default(none) reduction(+:cn) & + !$omp shared(self, mol, trans, cutoff2) & + !$omp private(jat, itr, izp, jzp, r2, rij, r1, den, countf) + do iat = 1, mol%nat + izp = mol%id(iat) + do jat = 1, iat + jzp = mol%id(jat) + den = self%get_en_factor(izp, jzp) + + 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) + + countf = den * self%ncoord_count(izp, jzp, r1) + + cn(iat) = cn(iat) + countf + if (iat /= jat) then + cn(jat) = cn(jat) + countf * self%directed_factor + end if + + end do + end do + end do + + end subroutine ncoord + + subroutine ncoord_d(self, mol, trans, cn, dcndr, dcndL) + !> Coordination number container + class(ncoord_type), intent(in) :: self + !> Molecular structure data + type(structure_type), intent(in) :: mol + !> Lattice points + real(wp), intent(in) :: trans(:, :) + !> Error function coordination number. + real(wp), intent(out) :: cn(:) + !> Derivative of the CN with respect to the Cartesian coordinates. + real(wp), intent(out) :: dcndr(:, :, :) + !> Derivative of the CN with respect to strain deformations. + real(wp), intent(out) :: dcndL(:, :, :) + + integer :: iat, jat, izp, jzp, itr + real(wp) :: r2, r1, rij(3), countf, countd(3), sigma(3, 3), cutoff2, den + + cn(:) = 0.0_wp + dcndr(:, :, :) = 0.0_wp + dcndL(:, :, :) = 0.0_wp + cutoff2 = self%cutoff**2 + + !$omp parallel do schedule(runtime) default(none) & + !$omp reduction(+:cn, dcndr, dcndL) shared(mol, trans, cutoff2) & + !$omp shared(self) & + !$omp private(jat, itr, izp, jzp, r2, rij, r1, den, countf, countd, sigma) + do iat = 1, mol%nat + izp = mol%id(iat) + do jat = 1, iat + jzp = mol%id(jat) + den = self%get_en_factor(izp, jzp) + + 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) + + countf = den * self%ncoord_count(izp, jzp, r1) + countd = den * self%ncoord_dcount(izp, jzp, r1) * rij/r1 + + cn(iat) = cn(iat) + countf + if (iat /= jat) then + cn(jat) = cn(jat) + countf * self%directed_factor + end if + + dcndr(:, iat, iat) = dcndr(:, iat, iat) + countd + dcndr(:, jat, jat) = dcndr(:, jat, jat) - countd * self%directed_factor + dcndr(:, iat, jat) = dcndr(:, iat, jat) + countd * self%directed_factor + dcndr(:, jat, iat) = dcndr(:, jat, iat) - countd + + sigma = spread(countd, 1, 3) * spread(rij, 2, 3) + + dcndL(:, :, iat) = dcndL(:, :, iat) + sigma + if (iat /= jat) then + dcndL(:, :, jat) = dcndL(:, :, jat) + sigma * self%directed_factor + end if + + end do + end do + end do + + end subroutine ncoord_d + + + subroutine add_coordination_number_derivs(self, mol, trans, dEdcn, gradient, sigma) + + !> Coordination number container + class(ncoord_type), intent(in) :: self + + !> Molecular structure data + type(structure_type), intent(in) :: mol + + !> Lattice points + real(wp), intent(in) :: trans(:, :) + + !> 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(:, :) + + integer :: iat, jat, izp, jzp, itr + real(wp) :: r2, r1, rc, rij(3), countf, countd(3), ds(3, 3), cutoff2, den + + cutoff2 = self%cutoff**2 + + !$omp parallel do schedule(runtime) default(none) & + !$omp reduction(+:gradient, sigma) & + !$omp shared(self, mol, trans, cutoff2, dEdcn) & + !$omp private(iat, jat, itr, izp, jzp, r2, rij, r1, rc, countd, ds, den) + do iat = 1, mol%nat + izp = mol%id(iat) + do jat = 1, iat + jzp = mol%id(jat) + den = self%get_en_factor(izp, jzp) + + 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) + + countd = den * self%ncoord_dcount(izp, jzp, r1) * rij/r1 + + gradient(:, iat) = gradient(:, iat) + countd & + & * (dEdcn(iat) + dEdcn(jat) * self%directed_factor) + gradient(:, jat) = gradient(:, jat) - countd & + & * (dEdcn(iat) + dEdcn(jat) * self%directed_factor) + + ds = spread(countd, 1, 3) * spread(rij, 2, 3) + + sigma(:, :) = sigma(:, :) & + & + ds * (dEdcn(iat) + & + & merge(dEdcn(jat) * self%directed_factor, 0.0_wp, jat /= iat)) + end do + end do + end do + + end subroutine add_coordination_number_derivs + + + !> Evaluates pairwise electronegativity factor if non applies + elemental function get_en_factor(self, izp, jzp) result(en_factor) + !> Coordination number container + class(ncoord_type), intent(in) :: self + !> Atom i index + integer, intent(in) :: izp + !> Atom j index + integer, intent(in) :: jzp + + real(wp) :: en_factor + + en_factor = 1.0_wp + + end function get_en_factor + + + !> Cutoff function for large coordination numbers + pure subroutine cut_coordination_number(cn_max, cn, dcndr, dcndL) + + !> Maximum CN (not strictly obeyed) + real(wp), intent(in) :: cn_max + + !> On input coordination number, on output modified CN + real(wp), intent(inout) :: cn(:) + + !> On input derivative of CN w.r.t. cartesian coordinates, + !> on output derivative of modified CN + real(wp), intent(inout), optional :: dcndr(:, :, :) + + !> On input derivative of CN w.r.t. strain deformation, + !> on output derivative of modified CN + real(wp), intent(inout), optional :: dcndL(:, :, :) + + real(wp) :: dcnpdcn + integer :: iat + + if (present(dcndL)) then + do iat = 1, size(cn) + dcnpdcn = dlog_cn_cut(cn(iat), cn_max) + dcndL(:, :, iat) = dcnpdcn*dcndL(:, :, iat) + enddo + endif + + if (present(dcndr)) then + do iat = 1, size(cn) + dcnpdcn = dlog_cn_cut(cn(iat), cn_max) + dcndr(:, :, iat) = dcnpdcn*dcndr(:, :, iat) + enddo + endif + + do iat = 1, size(cn) + cn(iat) = log_cn_cut(cn(iat), cn_max) + enddo + + end subroutine cut_coordination_number + + elemental function log_cn_cut(cn, cnmax) result(cnp) + real(wp), intent(in) :: cn + real(wp), intent(in) :: cnmax + real(wp) :: cnp + cnp = log(1.0_wp + exp(cnmax)) - log(1.0_wp + exp(cnmax - cn)) + end function log_cn_cut + + elemental function dlog_cn_cut(cn, cnmax) result(dcnpdcn) + real(wp), intent(in) :: cn + real(wp), intent(in) :: cnmax + real(wp) :: dcnpdcn + dcnpdcn = exp(cnmax)/(exp(cnmax) + exp(cn)) + end function dlog_cn_cut + +end module mctc_ncoord_type diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1509148b..f2d577b5 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -40,6 +40,7 @@ set( "write-turbomole" "write-vasp" "write-xyz" + "ncoord" ) set( test-srcs diff --git a/test/main.f90 b/test/main.f90 index b89bbb08..f2e64e09 100644 --- a/test/main.f90 +++ b/test/main.f90 @@ -43,6 +43,7 @@ program tester use test_write_turbomole, only : collect_write_turbomole use test_write_vasp, only : collect_write_vasp use test_write_xyz, only : collect_write_xyz + use test_ncoord, only : collect_ncoord implicit none integer :: stat, is character(len=:), allocatable :: suite_name, test_name @@ -76,7 +77,8 @@ program tester & new_testsuite("write-qchem", collect_write_qchem), & & new_testsuite("write-turbomole", collect_write_turbomole), & & new_testsuite("write-vasp", collect_write_vasp), & - & new_testsuite("write-xyz", collect_write_xyz) & + & new_testsuite("write-xyz", collect_write_xyz), & + & new_testsuite("ncoord", collect_ncoord) & & ] call get_argument(1, suite_name) diff --git a/test/meson.build b/test/meson.build index b2acfda0..51c3cb4a 100644 --- a/test/meson.build +++ b/test/meson.build @@ -38,6 +38,7 @@ tests = [ 'write-turbomole', 'write-vasp', 'write-xyz', + 'ncoord', ] test_srcs = files( diff --git a/test/test_ncoord.f90 b/test/test_ncoord.f90 new file mode 100644 index 00000000..4cba76cb --- /dev/null +++ b/test/test_ncoord.f90 @@ -0,0 +1,1714 @@ +! 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. + +module test_ncoord + use mctc_env, only : wp + use mctc_env_testing, only : new_unittest, unittest_type, error_type, & + & test_failed + use mctc_io_structure, only : structure_type + use testsuite_structure, only : get_structure + use mctc_cutoff, only : get_lattice_points + use mctc_data_covrad, only : get_covalent_rad + use mctc_data_paulingen, only : get_pauling_en + use mctc_ncoord_dexp, only : dexp_ncoord_type, new_dexp_ncoord + use mctc_ncoord_exp, only : exp_ncoord_type, new_exp_ncoord + use mctc_ncoord_erf, only : erf_ncoord_type, new_erf_ncoord + use mctc_ncoord_erf_en, only : erf_en_ncoord_type, new_erf_en_ncoord + use mctc_ncoord_erf_dftd4, only : erf_dftd4_ncoord_type, new_erf_dftd4_ncoord + use mctc_ncoord_type, only : ncoord_type + use mctc_ncoord, only : new_ncoord + implicit none + private + + public :: collect_ncoord + + real(wp), parameter :: thr = 100*epsilon(1.0_wp) + real(wp), parameter :: thr2 = sqrt(epsilon(1.0_wp)) + +contains + + + !> Collect all exported unit tests + subroutine collect_ncoord(testsuite) + + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: testsuite(:) + + testsuite = [ & + & new_unittest("cn-mb01_dexp", test_cn_mb01_dexp), & + & new_unittest("cn-mb02_dexp", test_cn_mb02_dexp), & + & new_unittest("cn-mb03_dexp", test_cn_mb03_dexp), & + & new_unittest("cn-acetic_dexp", test_cn_acetic_dexp), & + & new_unittest("dcndr-mb04_dexp", test_dcndr_mb04_dexp), & + & new_unittest("dcndr-mb05_dexp", test_dcndr_mb05_dexp), & + & new_unittest("dcndr-ammonia_dexp", test_dcndr_ammonia_dexp), & + & new_unittest("dcndL-mb06_dexp", test_dcndL_mb06_dexp), & + & new_unittest("dcndL-mb07_dexp", test_dcndL_mb07_dexp), & + & new_unittest("dcndL-antracene_dexp", test_dcndL_anthracene_dexp), & + & new_unittest("cn-mb01_exp", test_cn_mb01_exp), & + & new_unittest("cn-mb02_exp", test_cn_mb02_exp), & + & new_unittest("cn-mb03_exp", test_cn_mb03_exp), & + & new_unittest("cn-acetic_exp", test_cn_acetic_exp), & + & new_unittest("dcndr-mb04_exp", test_dcndr_mb04_exp), & + & new_unittest("dcndr-mb05_exp", test_dcndr_mb05_exp), & + & new_unittest("dcndr-ammonia_exp", test_dcndr_ammonia_exp), & + & new_unittest("dcndL-mb06_exp", test_dcndL_mb06_exp), & + & new_unittest("dcndL-mb07_exp", test_dcndL_mb07_exp), & + & new_unittest("dcndL-antracene_exp", test_dcndL_anthracene_exp), & + & new_unittest("cn-mb01_erf", test_cn_mb01_erf), & + & new_unittest("cn-mb01_cutoff_erf", test_cn_cutoff_mb01_erf), & + & new_unittest("cn-mb02_erf", test_cn_mb02_erf), & + & new_unittest("cn-mb03_erf", test_cn_mb03_erf), & + & new_unittest("cn-acetic_erf", test_cn_acetic_erf), & + & new_unittest("dcndr-mb04_erf", test_dcndr_mb04_erf), & + & new_unittest("dcndr-mb05_erf", test_dcndr_mb05_erf), & + & new_unittest("dcndr-ammonia_erf", test_dcndr_ammonia_erf), & + & new_unittest("dcndL-mb06_erf", test_dcndL_mb06_erf), & + & new_unittest("dcndL-mb07_erf", test_dcndL_mb07_erf), & + & new_unittest("dcndL-antracene_erf", test_dcndL_anthracene_erf), & + & new_unittest("cn-mb01_erf_en", test_cn_mb01_erf_en), & + & new_unittest("cn-mb02_erf_en", test_cn_mb02_erf_en), & + & new_unittest("cn-mb03_erf_en", test_cn_mb03_erf_en), & + & new_unittest("cn-acetic_erf_en", test_cn_acetic_erf_en), & + & new_unittest("dcndr-mb04_erf_en", test_dcndr_mb04_erf_en), & + & new_unittest("dcndr-mb05_erf_en", test_dcndr_mb05_erf_en), & + & new_unittest("dcndr-ammonia_erf_en", test_dcndr_ammonia_erf_en), & + & new_unittest("dcndL-mb06_erf_en", test_dcndL_mb06_erf_en), & + & new_unittest("dcndL-mb07_erf_en", test_dcndL_mb07_erf_en), & + & new_unittest("dcndL-antracene_erf_en", test_dcndL_anthracene_erf_en), & + & new_unittest("cn-mb01_erf_dftd4", test_cn_mb01_erf_dftd4), & + & new_unittest("dfdcn-mb01_erf_dftd4", test_dfdcn_mb01_erf_dftd4), & + & new_unittest("cn-mb02_erf_dftd4", test_cn_mb02_erf_dftd4), & + & new_unittest("cn-mb03_erf_dftd4", test_cn_mb03_erf_dftd4), & + & new_unittest("cn-acetic_erf_dftd4", test_cn_acetic_erf_dftd4), & + & new_unittest("dcndr-mb04_erf_dftd4", test_dcndr_mb04_erf_dftd4), & + & new_unittest("dcndr-mb05_erf_dftd4", test_dcndr_mb05_erf_dftd4), & + & new_unittest("dcndr-ammonia_erf_dftd4", test_dcndr_ammonia_erf_dftd4), & + & new_unittest("dcndL-mb06_erf_dftd4", test_dcndL_mb06_erf_dftd4), & + & new_unittest("dcndL-mb07_erf_dftd4", test_dcndL_mb07_erf_dftd4), & + & new_unittest("dcndL-antracene_erf_dftd4", test_dcndL_anthracene_erf_dftd4) & + & ] + + end subroutine collect_ncoord + + + subroutine test_cn_gen(error, mol, ncoord, ref) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + !> Molecular structure data + type(structure_type) :: mol + + !> Coordination number type + class(ncoord_type) :: ncoord + + !> Reference CNs + real(wp), intent(in) :: ref(:) + real(wp), allocatable :: cn(:) + real(wp), allocatable :: lattr(:, :) + + allocate(cn(mol%nat)) + + call get_lattice_points(mol%periodic, mol%lattice, ncoord%cutoff, lattr) + call ncoord%get_coordination_number(mol, lattr, cn) + + if (any(abs(cn - ref) > thr)) then + call test_failed(error, "Coordination numbers do not match") + print'(3es21.14)', cn + end if + + end subroutine test_cn_gen + + + subroutine test_numgrad(error, mol, ncoord) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + !> Coordination number type + class(ncoord_type) :: ncoord + + integer :: iat, ic + real(wp), allocatable :: cn(:), cnr(:), cnl(:) + real(wp), allocatable :: dcndr(:, :, :), dcndL(:, :, :) + real(wp), allocatable :: numdr(:, :, :) + real(wp), allocatable :: lattr(:, :) + real(wp), parameter :: step = 1.0e-6_wp + + allocate(cn(mol%nat), cnr(mol%nat), cnl(mol%nat), & + & dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & + & numdr(3, mol%nat, mol%nat)) + + call get_lattice_points(mol%periodic, mol%lattice, ncoord%cutoff, lattr) + + do iat = 1, mol%nat + do ic = 1, 3 + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + call ncoord%get_coordination_number(mol, lattr, cnr) + mol%xyz(ic, iat) = mol%xyz(ic, iat) - 2*step + call ncoord%get_coordination_number(mol, lattr, cnl) + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + numdr(ic, iat, :) = 0.5_wp*(cnr - cnl)/step + end do + end do + + call ncoord%get_coordination_number(mol, lattr, cn, dcndr, dcndL) + + if (any(abs(dcndr - numdr) > thr2)) then + call test_failed(error, "Derivative of coordination number does not match") + end if + + end subroutine test_numgrad + + + subroutine test_numsigma(error, mol, ncoord) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + !> Coordination number type + class(ncoord_type) :: ncoord + + integer :: ic, jc + real(wp) :: eps(3, 3) + real(wp), allocatable :: cn(:), cnr(:), cnl(:), xyz(:, :) + real(wp), allocatable :: dcndr(:, :, :), dcndL(:, :, :) + real(wp), allocatable :: numdL(:, :, :) + real(wp), allocatable :: lattr(:, :), trans(:, :) + real(wp), parameter :: unity(3, 3) = reshape(& + & [1, 0, 0, 0, 1, 0, 0, 0, 1], shape(unity)) + real(wp), parameter :: step = 1.0e-6_wp + + allocate(cn(mol%nat), cnr(mol%nat), cnl(mol%nat), & + & dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), xyz(3, mol%nat), & + & numdL(3, 3, mol%nat)) + call get_lattice_points(mol%periodic, mol%lattice, ncoord%cutoff, lattr) + + eps(:, :) = unity + xyz(:, :) = mol%xyz + trans = lattr + do ic = 1, 3 + do jc = 1, 3 + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = matmul(eps, xyz) + lattr(:, :) = matmul(eps, trans) + call ncoord%get_coordination_number(mol, lattr, cnr) + eps(jc, ic) = eps(jc, ic) - 2*step + mol%xyz(:, :) = matmul(eps, xyz) + lattr(:, :) = matmul(eps, trans) + call ncoord%get_coordination_number(mol, lattr, cnl) + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = xyz + lattr(:, :) = trans + numdL(jc, ic, :) = 0.5_wp*(cnr - cnl)/step + end do + end do + + call ncoord%get_coordination_number(mol, lattr, cn, dcndr, dcndL) + + if (any(abs(dcndL - numdL) > 10.0_wp * thr2)) then + call test_failed(error, "Derivative of coordination number does not match") + end if + + end subroutine test_numsigma + + + !> ---------------------------------------------------- + !> Tests for double-exponential coordination number + !> ---------------------------------------------------- + subroutine test_cn_mb01_dexp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(ncoord_type), allocatable :: dexp_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: cn(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref(16) = [& + & 4.11453659059991E+0_wp, 9.32058998762811E-1_wp, 2.03554597140311E+0_wp, & + & 1.42227835389358E+0_wp, 1.12812426574031E+0_wp, 1.05491602558828E+0_wp, & + & 1.52709064704269E+0_wp, 1.95070367247232E+0_wp, 3.83759889196540E+0_wp, & + & 1.09388314007182E+0_wp, 1.07090773695340E+0_wp, 2.00285254082830E+0_wp, & + & 4.36400837813955E+0_wp, 3.83469860546080E+0_wp, 3.91542517673963E+0_wp, & + & 5.58571682419960E+0_wp] + + call get_structure(mol, "mindless01") + + allocate(rcov(mol%nid), cn(mol%nat)) + rcov(:) = get_covalent_rad(mol%num) + ! Test also the external interface + call new_ncoord(dexp_ncoord, mol, "dexp", cutoff=cutoff, rcov=rcov) + call dexp_ncoord%get_cn(mol, cn) + + if (any(abs(cn - ref) > thr)) then + call test_failed(error, "Coordination numbers do not match") + print'(3es21.14)', cn + end if + + end subroutine test_cn_mb01_dexp + + + subroutine test_cn_mb02_dexp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(dexp_ncoord_type) :: dexp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref(16) = [& + & 9.16720394154780E-1_wp, 3.81678481096046E+0_wp, 3.58669504034282E+0_wp, & + & 2.90017371823910E+0_wp, 5.09896950232009E+0_wp, 1.13810142820063E+0_wp, & + & 9.35735928561309E-1_wp, 9.43324181507483E-1_wp, 4.86704300506033E+0_wp, & + & 1.22024405676881E+0_wp, 3.86661472534511E+0_wp, 4.02897853934353E+0_wp, & + & 1.96869999324996E+0_wp, 9.04495366214972E-1_wp, 1.49977675358986E+0_wp, & + & 2.04937665719480E+0_wp] + + call get_structure(mol, "mindless02") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_dexp_ncoord(dexp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_cn_gen(error, mol, dexp_ncoord, ref) + + end subroutine test_cn_mb02_dexp + + + subroutine test_cn_mb03_dexp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(dexp_ncoord_type) :: dexp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref(16) = [& + & 4.03003984525826E+0_wp, 2.88404829383272E+0_wp, 1.12030447706437E+0_wp, & + & 4.89326653038565E+0_wp, 6.97177319120287E+0_wp, 4.43722632690247E+0_wp, & + & 4.74866591812568E+0_wp, 1.30338583081493E+0_wp, 1.06130190313902E+0_wp, & + & 1.04048373066340E+0_wp, 1.92270060977307E+0_wp, 2.98737904327655E+0_wp, & + & 4.69081207436918E+0_wp, 1.26260125612717E+0_wp, 3.36368006785498E+0_wp, & + & 1.35743886389476E+0_wp] + + call get_structure(mol, "mindless03") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_dexp_ncoord(dexp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_cn_gen(error, mol, dexp_ncoord, ref) + + end subroutine test_cn_mb03_dexp + + + subroutine test_cn_acetic_dexp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(dexp_ncoord_type) :: dexp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref(32) = [& + & 1.29670485827103E+0_wp, 2.12821173792238E+0_wp, 1.29689958901536E+0_wp, & + & 2.12819492380168E+0_wp, 2.12820788692351E+0_wp, 1.29730145110972E+0_wp, & + & 1.29709466512645E+0_wp, 2.12822682780515E+0_wp, 3.12090335292641E+0_wp, & + & 4.03264857560058E+0_wp, 3.12090536181014E+0_wp, 4.03267505042477E+0_wp, & + & 4.03270042180352E+0_wp, 3.12093082635011E+0_wp, 3.12100563126373E+0_wp, & + & 4.03263801987003E+0_wp, 1.00224480975784E+0_wp, 1.11740483416248E+0_wp, & + & 1.00446562107381E+0_wp, 1.00138537125699E+0_wp, 1.00139516076697E+0_wp, & + & 1.00452681395494E+0_wp, 1.00222513359088E+0_wp, 1.11751807844285E+0_wp, & + & 1.00457316822347E+0_wp, 1.00217171056793E+0_wp, 1.11725423574458E+0_wp, & + & 1.00146807045654E+0_wp, 1.00452614353119E+0_wp, 1.00217802827893E+0_wp, & + & 1.00141181936828E+0_wp, 1.11720394408805E+0_wp] + + call get_structure(mol, "x02") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_dexp_ncoord(dexp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_cn_gen(error, mol, dexp_ncoord, ref) + + end subroutine test_cn_acetic_dexp + + + subroutine test_dcndr_mb04_dexp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(dexp_ncoord_type) :: dexp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless04") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_dexp_ncoord(dexp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numgrad(error, mol, dexp_ncoord) + + end subroutine test_dcndr_mb04_dexp + + + subroutine test_dcndr_mb05_dexp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(dexp_ncoord_type) :: dexp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless05") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_dexp_ncoord(dexp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numgrad(error, mol, dexp_ncoord) + + end subroutine test_dcndr_mb05_dexp + + + subroutine test_dcndr_ammonia_dexp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(dexp_ncoord_type) :: dexp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "x04") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_dexp_ncoord(dexp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numgrad(error, mol, dexp_ncoord) + + end subroutine test_dcndr_ammonia_dexp + + + subroutine test_dcndL_mb06_dexp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(dexp_ncoord_type) :: dexp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless06") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_dexp_ncoord(dexp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numsigma(error, mol, dexp_ncoord) + + end subroutine test_dcndL_mb06_dexp + + + subroutine test_dcndL_mb07_dexp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(dexp_ncoord_type) :: dexp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless07") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_dexp_ncoord(dexp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numsigma(error, mol, dexp_ncoord) + + end subroutine test_dcndL_mb07_dexp + + + subroutine test_dcndL_anthracene_dexp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(dexp_ncoord_type) :: dexp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "x05") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_dexp_ncoord(dexp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numsigma(error, mol, dexp_ncoord) + + end subroutine test_dcndL_anthracene_dexp + + + !> ---------------------------------------------------- + !> Tests for mono-exponential coordination number + !> ---------------------------------------------------- + subroutine test_cn_mb01_exp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(ncoord_type), allocatable :: exp_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: cn(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref(16) = [& + & 4.15066368951397E+0_wp, 9.78868026389781E-1_wp, 2.01080985633859E+0_wp, & + & 1.47865697827818E+0_wp, 1.03577822442117E+0_wp, 1.01206994314781E+0_wp, & + & 1.50329777127401E+0_wp, 1.99858468272609E+0_wp, 3.89181927539324E+0_wp, & + & 1.04323373360740E+0_wp, 1.01526584450636E+0_wp, 1.99315213227354E+0_wp, & + & 4.63526560889683E+0_wp, 3.87312260639335E+0_wp, 3.99316800677884E+0_wp, & + & 5.45068226903888E+0_wp] + + call get_structure(mol, "mindless01") + + allocate(rcov(mol%nid), cn(mol%nat)) + rcov(:) = get_covalent_rad(mol%num) + + ! Test also the external interface + call new_ncoord(exp_ncoord, mol, "exp", cutoff=cutoff, rcov=rcov) + call exp_ncoord%get_cn(mol, cn) + + if (any(abs(cn - ref) > thr)) then + call test_failed(error, "Coordination numbers do not match") + print'(3es21.14)', cn + end if + + end subroutine test_cn_mb01_exp + + + subroutine test_cn_mb02_exp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(exp_ncoord_type) :: exp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref(16) = [& + & 9.68434565947196E-1_wp, 3.94488220883276E+0_wp, 3.82701677880409E+0_wp, & + & 2.99161201243234E+0_wp, 5.46904892971914E+0_wp, 1.06438740698832E+0_wp, & + & 9.77627896999762E-1_wp, 9.81715643929621E-1_wp, 5.06702924169705E+0_wp, & + & 1.08324093335500E+0_wp, 4.00161384251868E+0_wp, 4.03067393311321E+0_wp, & + & 2.00249301823990E+0_wp, 9.73399178780401E-1_wp, 1.66868575900646E+0_wp, & + & 2.03273936244268E+0_wp] + + call get_structure(mol, "mindless02") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_exp_ncoord(exp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_cn_gen(error, mol, exp_ncoord, ref) + + end subroutine test_cn_mb02_exp + + + subroutine test_cn_mb03_exp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(exp_ncoord_type) :: exp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref(16) = [& + & 4.04992403668766E+0_wp, 2.98487291056010E+0_wp, 1.03236609075831E+0_wp, & + & 4.86782876076800E+0_wp, 7.48154833549122E+0_wp, 4.76588477343835E+0_wp, & + & 4.92432613481557E+0_wp, 1.19761963808520E+0_wp, 1.01809037129368E+0_wp, & + & 1.01042713594272E+0_wp, 1.99186789992505E+0_wp, 3.03157225205500E+0_wp, & + & 4.89702969622395E+0_wp, 1.11663432026701E+0_wp, 3.52903544775683E+0_wp, & + & 1.33563958333433E+0_wp] + + call get_structure(mol, "mindless03") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_exp_ncoord(exp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_cn_gen(error, mol, exp_ncoord, ref) + + end subroutine test_cn_mb03_exp + + + subroutine test_cn_acetic_exp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(exp_ncoord_type) :: exp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref(32) = [& + & 1.08784433872404E+0_wp, 2.04822026776264E+0_wp, 1.08792777906791E+0_wp, & + & 2.04820799570728E+0_wp, 2.04821951599101E+0_wp, 1.08811939753571E+0_wp, & + & 1.08802819063717E+0_wp, 2.04823763269824E+0_wp, 3.03590173052949E+0_wp, & + & 4.04097072127914E+0_wp, 3.03590425814990E+0_wp, 4.04098287934520E+0_wp, & + & 4.04097935270722E+0_wp, 3.03590447694319E+0_wp, 3.03593501518861E+0_wp, & + & 4.04096875806481E+0_wp, 1.00135680977725E+0_wp, 1.03683422933766E+0_wp, & + & 1.00190020398973E+0_wp, 1.00114101857689E+0_wp, 1.00114344249343E+0_wp, & + & 1.00191213773267E+0_wp, 1.00135278116294E+0_wp, 1.03690135098048E+0_wp, & + & 1.00192295753089E+0_wp, 1.00133927917777E+0_wp, 1.03675374108303E+0_wp, & + & 1.00115868342640E+0_wp, 1.00191299668631E+0_wp, 1.00134039650187E+0_wp, & + & 1.00114761354986E+0_wp, 1.03670619005681E+0_wp] + + call get_structure(mol, "x02") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_exp_ncoord(exp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_cn_gen(error, mol, exp_ncoord, ref) + + end subroutine test_cn_acetic_exp + + + subroutine test_dcndr_mb04_exp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(exp_ncoord_type) :: exp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless04") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_exp_ncoord(exp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numgrad(error, mol, exp_ncoord) + + end subroutine test_dcndr_mb04_exp + + + subroutine test_dcndr_mb05_exp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(exp_ncoord_type) :: exp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless05") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_exp_ncoord(exp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numgrad(error, mol, exp_ncoord) + + end subroutine test_dcndr_mb05_exp + + + subroutine test_dcndr_ammonia_exp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(exp_ncoord_type) :: exp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "x04") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_exp_ncoord(exp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numgrad(error, mol, exp_ncoord) + + end subroutine test_dcndr_ammonia_exp + + + subroutine test_dcndL_mb06_exp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(exp_ncoord_type) :: exp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless06") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_exp_ncoord(exp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numsigma(error, mol, exp_ncoord) + + end subroutine test_dcndL_mb06_exp + + + subroutine test_dcndL_mb07_exp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(exp_ncoord_type) :: exp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless07") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_exp_ncoord(exp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numsigma(error, mol, exp_ncoord) + + end subroutine test_dcndL_mb07_exp + + + subroutine test_dcndL_anthracene_exp(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(exp_ncoord_type) :: exp_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "x05") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_exp_ncoord(exp_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numsigma(error, mol, exp_ncoord) + + end subroutine test_dcndL_anthracene_exp + + + !> ---------------------------------------------------- + !> Tests for error-function based coordination number + !> using the Pyykko covalent radii and Pauling EN + !> ---------------------------------------------------- + subroutine test_cn_mb01_erf(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(ncoord_type), allocatable :: erf_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: cn(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: kcn = 2.60_wp + real(wp), parameter :: ref(16) = [& + & +4.02932551939856E+00_wp, +7.85103341426644E-01_wp, +1.81014102249862E+00_wp, & + & +1.27570668822915E+00_wp, +1.07819362139425E+00_wp, +9.25571697114998E-01_wp, & + & +1.46154464565793E+00_wp, +1.68476840890299E+00_wp, +3.36480649393304E+00_wp, & + & +1.01246571460105E+00_wp, +9.41673850989885E-01_wp, +1.86611893141095E+00_wp, & + & +3.82207104984211E+00_wp, +3.43775046556670E+00_wp, +3.43498366851429E+00_wp, & + & +5.23834442955556E+00_wp] + + call get_structure(mol, "mindless01") + + allocate(rcov(mol%nid), cn(mol%nat)) + rcov(:) = get_covalent_rad(mol%num) + + ! Test also the external interface + call new_ncoord(erf_ncoord, mol, "erf", kcn=kcn, cutoff=cutoff, rcov=rcov) + call erf_ncoord%get_cn(mol, cn) + + if (any(abs(cn - ref) > thr)) then + call test_failed(error, "Coordination numbers do not match") + print'(3es21.14)', cn + end if + + end subroutine test_cn_mb01_erf + + subroutine test_cn_cutoff_mb01_erf(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(ncoord_type), allocatable :: erf_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: cn(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: kcn = 7.5_wp + real(wp), parameter :: cn_max = 8.0_wp + real(wp), parameter :: ref(16) = [& + & 4.01821724338495E+0_wp, 9.72247109465167E-1_wp, 1.98487635195765E+0_wp, & + & 1.47199898291333E+0_wp, 9.96978323172462E-1_wp, 9.96288834201663E-1_wp, & + & 1.45078799827281E+0_wp, 1.99055055374771E+0_wp, 3.83042334184661E+0_wp, & + & 1.00185131973923E+0_wp, 9.96142041368325E-1_wp, 1.92309124625511E+0_wp, & + & 4.58700117126191E+0_wp, 3.80489139049207E+0_wp, 3.94005009430351E+0_wp, & + & 5.27144183788523E+0_wp] + + call get_structure(mol, "mindless01") + + allocate(rcov(mol%nid), cn(mol%nat)) + rcov(:) = get_covalent_rad(mol%num) + + ! Test also the external interface + call new_ncoord(erf_ncoord, mol, "erf", kcn=kcn, cutoff=cutoff, rcov=rcov, cut=cn_max) + call erf_ncoord%get_cn(mol, cn) + + if (any(abs(cn - ref) > thr)) then + call test_failed(error, "Coordination numbers do not match") + print'(3es21.14)', cn + end if + + end subroutine test_cn_cutoff_mb01_erf + + + subroutine test_cn_mb02_erf(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_ncoord_type) :: erf_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: kcn = 2.60_wp + real(wp), parameter :: norm_exp = 0.8_wp + real(wp), parameter :: ref(16) = [& + & 7.82883150684095E-01_wp, 3.47414568086695E+00_wp, 3.21490904495080E+00_wp, & + & 2.57204039419665E+00_wp, 4.70139897809430E+00_wp, 1.01600663641776E+00_wp, & + & 8.02988129457653E-01_wp, 8.05575211220716E-01_wp, 4.53312550965630E+00_wp, & + & 1.07341286735542E+00_wp, 3.48131954350286E+00_wp, 3.75872922171681E+00_wp, & + & 1.83133763195600E+00_wp, 8.02241219305798E-01_wp, 1.33431699912730E+00_wp, & + & 1.85405528497185E+00_wp] + + call get_structure(mol, "mindless02") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_erf_ncoord(erf_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov, norm_exp=norm_exp) + call test_cn_gen(error, mol, erf_ncoord, ref) + + end subroutine test_cn_mb02_erf + + + subroutine test_cn_mb03_erf(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_ncoord_type) :: erf_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: kcn = 2.60_wp + real(wp), parameter :: ref(16) = [& + & +3.65604858453391E+00_wp, 2.47444042434446E+00_wp, 1.04729945244460E+00_wp, & + & +4.69028664318902E+00_wp, 6.13331280895969E+00_wp, 3.96706985488549E+00_wp, & + & +4.19093957903296E+00_wp, 1.24758706681051E+00_wp, 9.56572882868544E-01_wp, & + & +9.15465255206712E-01_wp, 1.65382740011051E+00_wp, 2.71619299069121E+00_wp, & + & +4.13366036611672E+00_wp, 1.25908852786395E+00_wp, 3.12414929301946E+00_wp, & + & +1.24550167348772E+00_wp] + + call get_structure(mol, "mindless03") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_erf_ncoord(erf_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov) + call test_cn_gen(error, mol, erf_ncoord, ref) + + end subroutine test_cn_mb03_erf + + + subroutine test_cn_acetic_erf(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_ncoord_type) :: erf_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: kcn = 2.60_wp + real(wp), parameter :: ref(32) = [& + & 1.32759422973205E+00_wp, 1.94146232311076E+00_wp, 1.32785842732676E+00_wp, & + & 1.94145796402605E+00_wp, 1.94144168404587E+00_wp, 1.32846648409150E+00_wp, & + & 1.32817586952618E+00_wp, 1.94146973878954E+00_wp, 2.77909305549586E+00_wp, & + & 3.54216537574710E+00_wp, 2.77910670281296E+00_wp, 3.54220807172730E+00_wp, & + & 3.54223833711160E+00_wp, 2.77910422212159E+00_wp, 2.77925865658396E+00_wp, & + & 3.54214599507453E+00_wp, 8.43683239278435E-01_wp, 1.03961906362830E+00_wp, & + & 8.45225661192366E-01_wp, 8.42052731344644E-01_wp, 8.42069624420118E-01_wp, & + & 8.45292837568377E-01_wp, 8.43651700598541E-01_wp, 1.03980036426124E+00_wp, & + & 8.45340647438754E-01_wp, 8.43551963743648E-01_wp, 1.03938961293610E+00_wp, & + & 8.42194599506036E-01_wp, 8.45292255082151E-01_wp, 8.43583848813426E-01_wp, & + & 8.42107560458205E-01_wp, 1.03929415502029E+00_wp] + + call get_structure(mol, "x02") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_erf_ncoord(erf_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov) + call test_cn_gen(error, mol, erf_ncoord, ref) + + end subroutine test_cn_acetic_erf + + + subroutine test_dcndr_mb04_erf(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_ncoord_type) :: erf_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless04") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_erf_ncoord(erf_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numgrad(error, mol, erf_ncoord) + + end subroutine test_dcndr_mb04_erf + + + subroutine test_dcndr_mb05_erf(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_ncoord_type) :: erf_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: norm_exp = 0.8_wp + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless05") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_erf_ncoord(erf_ncoord, mol, cutoff=cutoff, & + & rcov=rcov, norm_exp=norm_exp) + call test_numgrad(error, mol, erf_ncoord) + + end subroutine test_dcndr_mb05_erf + + + subroutine test_dcndr_ammonia_erf(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_ncoord_type) :: erf_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "x04") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_erf_ncoord(erf_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numgrad(error, mol, erf_ncoord) + + end subroutine test_dcndr_ammonia_erf + + + subroutine test_dcndL_mb06_erf(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_ncoord_type) :: erf_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless06") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_erf_ncoord(erf_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numsigma(error, mol, erf_ncoord) + + end subroutine test_dcndL_mb06_erf + + + subroutine test_dcndL_mb07_erf(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_ncoord_type) :: erf_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless07") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_erf_ncoord(erf_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numsigma(error, mol, erf_ncoord) + + end subroutine test_dcndL_mb07_erf + + + subroutine test_dcndL_anthracene_erf(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_ncoord_type) :: erf_ncoord + real(wp), allocatable :: rcov(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "x05") + + allocate(rcov(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + + call new_erf_ncoord(erf_ncoord, mol, cutoff=cutoff, rcov=rcov) + call test_numsigma(error, mol, erf_ncoord) + + end subroutine test_dcndL_anthracene_erf + + + !> ---------------------------------------------------- + !> Tests for electronegativity-weighted + !> error-function-based coordination number + !> using the Pyykko covalent radii and Pauling EN + !> ---------------------------------------------------- + subroutine test_cn_mb01_erf_en(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(ncoord_type), allocatable :: erf_en_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + real(wp), allocatable :: cn(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: kcn = 2.60_wp + real(wp), parameter :: ref(16) = [& + & +6.49835771218084E+00_wp, -9.24780426380943E-02_wp, -3.19404023806046E+00_wp, & + & -7.26463405313327E-01_wp, -1.92171552979195E+00_wp, +6.71703342153316E-01_wp, & + & -6.47861832477202E-01_wp, -2.56263447944906E+00_wp, -3.28660571935921E+00_wp, & + & +8.86404306801211E-01_wp, +6.65084753274741E-01_wp, -2.51172025807030E+00_wp, & + & +4.24336914659746E-01_wp, +2.98914644571775E+00_wp, -2.72992153136905E+00_wp, & + & +5.53840756174106E+00_wp] + + call get_structure(mol, "mindless01") + + allocate(rcov(mol%nid), en(mol%nid), cn(mol%nat)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + ! Test also the external interface + call new_ncoord(erf_en_ncoord, mol, "erf_en", & + & kcn=kcn, cutoff=cutoff, rcov=rcov, en=en) + call erf_en_ncoord%get_cn(mol, cn) + + if (any(abs(cn - ref) > thr)) then + call test_failed(error, "Coordination numbers do not match") + print'(3es21.14)', cn + end if + + end subroutine test_cn_mb01_erf_en + + + subroutine test_cn_mb02_erf_en(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_en_ncoord_type) :: erf_en_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: kcn = 2.60_wp + real(wp), parameter :: norm_exp = 0.8_wp + real(wp), parameter :: ref(16) = [& + & -1.20766490409683E-01_wp, -3.16431597148754E+00_wp, -5.05004208597497E-02_wp, & + & -4.25983254262925E+00_wp, 5.12905009209427E+00_wp, 8.28195323146111E-01_wp, & + & -1.20642619114040E-01_wp, -1.30015328757524E-01_wp, 1.04018711784936E+00_wp, & + & 8.14999806238091E-01_wp, 6.97386143793777E-01_wp, 6.83247561402713E+00_wp, & + & -5.18157773756499E+00_wp, -2.40202426638317E-01_wp, -3.02851981491491E-01_wp, & + & -1.77158857819614E+00_wp] + + call get_structure(mol, "mindless02") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_en_ncoord(erf_en_ncoord, mol, kcn=kcn, cutoff=cutoff, & + & rcov=rcov, en=en, norm_exp=norm_exp) + call test_cn_gen(error, mol, erf_en_ncoord, ref) + + end subroutine test_cn_mb02_erf_en + + + subroutine test_cn_mb03_erf_en(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_en_ncoord_type) :: erf_en_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: kcn = 2.60_wp + real(wp), parameter :: ref(16) = [& + & -1.57129310680311E+00_wp, -5.12538648600132E+00_wp, +3.44404010675507E-02_wp, & + & +5.53888057043349E+00_wp, +6.08403921871331E+00_wp, +1.70484401191915E+00_wp, & + & -3.05797643039751E+00_wp, -1.77339255176255E-01_wp, +1.77101353984179E-01_wp, & + & +3.11611982445210E-01_wp, -4.72001776406229E+00_wp, -1.75545403486359E+00_wp, & + & -2.87572032143436E+00_wp, +4.45164095433560E-01_wp, +5.12260227693958E+00_wp, & + & -1.35496512197608E-01_wp] + + call get_structure(mol, "mindless03") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_en_ncoord(erf_en_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov, en=en) + call test_cn_gen(error, mol, erf_en_ncoord, ref) + + end subroutine test_cn_mb03_erf_en + + + subroutine test_cn_acetic_erf_en(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_en_ncoord_type) :: erf_en_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: kcn = 2.60_wp + real(wp), parameter :: ref(32) = [& + & -1.10960425359405E+00_wp, -1.87307032436238E+00_wp,-1.10988726298238E+00_wp, & + & -1.87300941691613E+00_wp, -1.87307851037536E+00_wp,-1.11055533892111E+00_wp, & + & -1.11022640841129E+00_wp, -1.87314182780441E+00_wp, 1.48844491507249E+00_wp, & + & -5.92197262250447E-01_wp, 1.48846528050850E+00_wp,-5.92187160672111E-01_wp, & + & -5.92267911690038E-01_wp, 1.48840420948914E+00_wp, 1.48839664548042E+00_wp, & + & -5.92201115944998E-01_wp, 2.95283098472335E-01_wp, 1.20097019987756E+00_wp, & + & 2.96054614475136E-01_wp, 2.94737814418179E-01_wp, 2.94743337206434E-01_wp, & + & 2.96081399455567E-01_wp, 2.95272376357621E-01_wp, 1.20121563992199E+00_wp, & + & 2.96099722915932E-01_wp, 2.95237799928904E-01_wp, 1.20068087168607E+00_wp, & + & 2.94786022309938E-01_wp, 2.96081562273217E-01_wp, 2.95248938198809E-01_wp, & + & 2.94755972287658E-01_wp, 1.20046637358882E+00_wp] + + call get_structure(mol, "x02") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_en_ncoord(erf_en_ncoord, mol, kcn=kcn, cutoff=cutoff, rcov=rcov, en=en) + call test_cn_gen(error, mol, erf_en_ncoord, ref) + + end subroutine test_cn_acetic_erf_en + + + subroutine test_dcndr_mb04_erf_en(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_en_ncoord_type) :: erf_en_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless04") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_en_ncoord(erf_en_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_numgrad(error, mol, erf_en_ncoord) + + end subroutine test_dcndr_mb04_erf_en + + + subroutine test_dcndr_mb05_erf_en(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_en_ncoord_type) :: erf_en_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: norm_exp = 0.8_wp + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless05") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_en_ncoord(erf_en_ncoord, mol, cutoff=cutoff, & + & rcov=rcov, en=en, norm_exp=norm_exp) + call test_numgrad(error, mol, erf_en_ncoord) + + end subroutine test_dcndr_mb05_erf_en + + + subroutine test_dcndr_ammonia_erf_en(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_en_ncoord_type) :: erf_en_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "x04") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_en_ncoord(erf_en_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_numgrad(error, mol, erf_en_ncoord) + + end subroutine test_dcndr_ammonia_erf_en + + + subroutine test_dcndL_mb06_erf_en(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_en_ncoord_type) :: erf_en_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless06") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_en_ncoord(erf_en_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_numsigma(error, mol, erf_en_ncoord) + + end subroutine test_dcndL_mb06_erf_en + + + subroutine test_dcndL_mb07_erf_en(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_en_ncoord_type) :: erf_en_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless07") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_en_ncoord(erf_en_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_numsigma(error, mol, erf_en_ncoord) + + end subroutine test_dcndL_mb07_erf_en + + + subroutine test_dcndL_anthracene_erf_en(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_en_ncoord_type) :: erf_en_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "x05") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_en_ncoord(erf_en_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_numsigma(error, mol, erf_en_ncoord) + + end subroutine test_dcndL_anthracene_erf_en + + + !> ---------------------------------------------------- + !> Tests for DFT-D4 coordination number + !> using the Pyykko covalent radii and Pauling EN + !> ---------------------------------------------------- + subroutine test_cn_mb01_erf_dftd4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(ncoord_type), allocatable :: erf_dftd4_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + real(wp), allocatable :: cn(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref(16) = [& + & 3.07349677110402E+0_wp, 9.31461605116103E-1_wp, 1.43709439375839E+0_wp, & + & 1.33309431581960E+0_wp, 7.20743527030337E-1_wp, 8.59659004770982E-1_wp, & + & 1.35782158177921E+0_wp, 1.53940006996025E+0_wp, 3.19400368195259E+0_wp, & + & 8.12162111631342E-1_wp, 8.59533443784854E-1_wp, 1.53347108155587E+0_wp, & + & 4.23314989525721E+0_wp, 3.03048504567396E+0_wp, 3.45229319488306E+0_wp, & + & 4.28478289652264E+0_wp] + + call get_structure(mol, "mindless01") + + allocate(rcov(mol%nid), en(mol%nid), cn(mol%nat)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + ! Test also the external interface + call new_ncoord(erf_dftd4_ncoord, mol, "dftd4", & + & cutoff=cutoff, rcov=rcov, en=en) + call erf_dftd4_ncoord%get_cn(mol, cn) + + if (any(abs(cn - ref) > thr)) then + call test_failed(error, "Coordination numbers do not match") + print'(3es21.14)', cn + end if + + end subroutine test_cn_mb01_erf_dftd4 + + subroutine test_dfdcn_mb01_erf_dftd4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(ncoord_type), allocatable :: erf_dftd4_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + real(wp) :: gradient(3, 16), sigma(3,3), dEdcn(16) + real(wp), allocatable :: cn(:) + real(wp), allocatable :: lattr(:, :) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref_gradient(3, 16) = reshape([ & + & 2.3515337861584880_wp, -3.0453827536691351_wp, 0.5392536733971742_wp, & + & -0.4086227335489838_wp, 0.0501074408829545_wp, 0.2137313744331271_wp, & + & 0.0249303991608369_wp, -0.0373339882416463_wp, -0.0756057413224991_wp, & + & -0.5436172455452489_wp, 0.4933244472902760_wp, 1.2499429094951524_wp, & + & 0.0187062301843386_wp, 0.0079863446319139_wp, -0.0342050060044000_wp, & + & 0.0541178064781164_wp, 0.0314514945594538_wp, 0.0256701683856457_wp, & + & -1.1131746078453495_wp, -1.6288054931512879_wp, 0.4146928754871603_wp, & + & 0.0887505649437715_wp, 0.0350320577127468_wp, -0.0407677557149725_wp, & + & -0.0499761354250785_wp, 0.1548143128637227_wp, -0.2724294261044626_wp, & + & -0.2140027394072616_wp, -0.0253209254791824_wp, -0.0023845891460968_wp, & + & 0.0091192286107324_wp, -0.0054889861418698_wp, -0.0708601381531837_wp, & + & 0.1710608463902656_wp, 0.5486016273258228_wp, 0.0456750463319099_wp, & + & -1.4496220571725840_wp, 0.2848532829040934_wp, 2.1631929338961609_wp, & + & 0.2437819446177886_wp, 0.0864993947414499_wp, -0.8074962768633803_wp, & + & 0.0760143426254694_wp, 0.2932453158431059_wp, -0.1384011976765554_wp, & + & 0.7410003697746987_wp, 2.7564164279275825_wp, -3.2100088504407793_wp], & + & shape(ref_gradient)) + real(wp), parameter :: ref_sigma(3, 3) = reshape([ & + &-14.102412014163017_wp, 5.058593018200457_wp, 7.2112237247416280_wp, & + & 5.058593018200459_wp,-18.111874478203490_wp, 6.8763096498221916_wp, & + & 7.211223724741628_wp, 6.876309649822192_wp, -22.636906090720050_wp], & + & shape(ref_sigma)) + + dEdcn = 1.0_wp + gradient = 0.0_wp + sigma = 0.0_wp + + call get_structure(mol, "mindless01") + + allocate(rcov(mol%nid), en(mol%nid), cn(mol%nat)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + ! Test also the external interface + call new_ncoord(erf_dftd4_ncoord, mol, "dftd4", & + & cutoff=cutoff, rcov=rcov, en=en) + call get_lattice_points(mol%periodic, mol%lattice, cutoff, lattr) + call erf_dftd4_ncoord%add_coordination_number_derivs(mol, lattr, & + & dEdcn, gradient, sigma) + + if (any(abs(gradient - ref_gradient) > thr)) then + call test_failed(error, "Coordination number gradient does not match") + print'(3es21.14)', gradient - ref_gradient + end if + + if (any(abs(sigma - ref_sigma) > thr)) then + call test_failed(error, "Coordination numbers sigma does not match") + print'(3es21.14)', sigma - ref_sigma + end if + + end subroutine test_dfdcn_mb01_erf_dftd4 + + subroutine test_cn_mb02_erf_dftd4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_dftd4_ncoord_type) :: erf_dftd4_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: norm_exp = 0.8_wp + real(wp), parameter :: ref(16) = [& + & 9.42718287075170E-1_wp, 3.36848356833576E+0_wp, 3.63463906109694E+0_wp, & + & 2.27670250096768E+0_wp, 4.71481260163654E+0_wp, 8.10661654309892E-1_wp, & + & 9.47358382643435E-1_wp, 9.48316042740831E-1_wp, 4.64606177656574E+0_wp, & + & 8.10384366749581E-1_wp, 3.74392916782297E+0_wp, 2.93620819522544E+0_wp, & + & 1.24273507745740E+0_wp, 9.30562514461567E-1_wp, 1.67284106029936E+0_wp, & + & 1.68882066428832E+0_wp] + + call get_structure(mol, "mindless02") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_dftd4_ncoord(erf_dftd4_ncoord, mol, cutoff=cutoff, & + & rcov=rcov, en=en, norm_exp=norm_exp) + call test_cn_gen(error, mol, erf_dftd4_ncoord, ref) + + end subroutine test_cn_mb02_erf_dftd4 + + + subroutine test_cn_mb03_erf_dftd4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_dftd4_ncoord_type) :: erf_dftd4_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref(16) = [& + & 3.70329575672631E+0_wp, 2.11476582851627E+0_wp, 9.23682826697708E-1_wp, & + & 3.88999767055963E+0_wp, 6.17490081567969E+0_wp, 4.10595506858888E+0_wp, & + & 4.22916534938705E+0_wp, 1.04287687415719E+0_wp, 9.23686985143960E-1_wp, & + & 9.24487848931390E-1_wp, 1.22927636800717E+0_wp, 2.59853001985457E+0_wp, & + & 4.30015470969650E+0_wp, 8.29081650895103E-1_wp, 2.65239010637793E+0_wp, & + & 1.19840431336618E+0_wp] + + call get_structure(mol, "mindless03") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_dftd4_ncoord(erf_dftd4_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_cn_gen(error, mol, erf_dftd4_ncoord, ref) + + end subroutine test_cn_mb03_erf_dftd4 + + + subroutine test_cn_acetic_erf_dftd4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_dftd4_ncoord_type) :: erf_dftd4_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + real(wp), parameter :: ref(32) = [& + & 8.57927937964499E-1_wp, 1.65578811302979E+0_wp, 8.57942519815764E-01_wp, & + & 1.65577032315393E+0_wp, 1.65579750428472E+0_wp, 8.57976292566367E-01_wp, & + & 8.57959806712906E-1_wp, 1.65580595790859E+0_wp, 2.68757316963446E+00_wp, & + & 3.75387937835847E+0_wp, 2.68757320525778E+0_wp, 3.75388279785624E+00_wp, & + & 3.75389467927099E+0_wp, 2.68757350196962E+0_wp, 2.68757458902025E+00_wp, & + & 3.75387626586114E+0_wp, 9.24588078703746E-1_wp, 8.00786883779774E-01_wp, & + & 9.25187896946826E-1_wp, 9.24714572053233E-1_wp, 9.24714022896347E-01_wp, & + & 9.25191607736594E-1_wp, 9.24588475035834E-1_wp, 8.00785878187927E-01_wp, & + & 9.25198590606433E-1_wp, 9.24591610342463E-1_wp, 8.00781165120561E-01_wp, & + & 9.24717032271544E-1_wp, 9.25190368617526E-1_wp, 9.24586232339746E-01_wp, & + & 9.24711233829591E-1_wp, 8.00775373869118E-1_wp] + + call get_structure(mol, "x02") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_dftd4_ncoord(erf_dftd4_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_cn_gen(error, mol, erf_dftd4_ncoord, ref) + + end subroutine test_cn_acetic_erf_dftd4 + + + subroutine test_dcndr_mb04_erf_dftd4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_dftd4_ncoord_type) :: erf_dftd4_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless04") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_dftd4_ncoord(erf_dftd4_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_numgrad(error, mol, erf_dftd4_ncoord) + + end subroutine test_dcndr_mb04_erf_dftd4 + + + subroutine test_dcndr_mb05_erf_dftd4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_dftd4_ncoord_type) :: erf_dftd4_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless05") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_dftd4_ncoord(erf_dftd4_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_numgrad(error, mol, erf_dftd4_ncoord) + + end subroutine test_dcndr_mb05_erf_dftd4 + + + subroutine test_dcndr_ammonia_erf_dftd4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_dftd4_ncoord_type) :: erf_dftd4_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "x04") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_dftd4_ncoord(erf_dftd4_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_numgrad(error, mol, erf_dftd4_ncoord) + + end subroutine test_dcndr_ammonia_erf_dftd4 + + + subroutine test_dcndL_mb06_erf_dftd4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_dftd4_ncoord_type) :: erf_dftd4_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless06") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_dftd4_ncoord(erf_dftd4_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_numsigma(error, mol, erf_dftd4_ncoord) + + end subroutine test_dcndL_mb06_erf_dftd4 + + + subroutine test_dcndL_mb07_erf_dftd4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_dftd4_ncoord_type) :: erf_dftd4_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "mindless07") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_dftd4_ncoord(erf_dftd4_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_numsigma(error, mol, erf_dftd4_ncoord) + + end subroutine test_dcndL_mb07_erf_dftd4 + + + subroutine test_dcndL_anthracene_erf_dftd4(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + type(erf_dftd4_ncoord_type) :: erf_dftd4_ncoord + real(wp), allocatable :: rcov(:) + real(wp), allocatable :: en(:) + + real(wp), parameter :: cutoff = 30.0_wp + + call get_structure(mol, "x05") + + allocate(rcov(mol%nid), en(mol%nid)) + rcov(:) = get_covalent_rad(mol%num) + en(:) = get_pauling_en(mol%num) + + call new_erf_dftd4_ncoord(erf_dftd4_ncoord, mol, cutoff=cutoff, rcov=rcov, en=en) + call test_numsigma(error, mol, erf_dftd4_ncoord) + + end subroutine test_dcndL_anthracene_erf_dftd4 + + +end module test_ncoord diff --git a/test/testsuite_structure.f90 b/test/testsuite_structure.f90 index dd642459..95b50fab 100644 --- a/test/testsuite_structure.f90 +++ b/test/testsuite_structure.f90 @@ -532,65 +532,65 @@ end subroutine x04 subroutine x05(self) type(structure_type), intent(out) :: self integer, parameter :: nat = 48 - character(len=*), parameter :: sym(nat) = [character(len=4) ::& - & "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", & - & "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", & + character(len=*), parameter :: sym(nat) = [character(len=4)::& & "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", & - & "H", "H", "H", "H", "H", "H"] + & "H", "H", "H", "H", "H", "H", "C", "C", "C", "C", "C", "C", "C", "C", & + & "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", & + & "C", "C", "C", "C", "C", "C"] real(wp), parameter :: xyz(3, nat) = reshape([& - & 9.36760291026725_wp, 7.80244140446007_wp, 2.94420668811823_wp, & - & -2.60635659089003_wp, 11.13901918927945_wp, 16.83643328201347_wp, & - & 7.22248107756189_wp, 8.57559920837242_wp, 7.03672928400636_wp, & - & -1.76362632094559_wp, 0.64297750017893_wp, 14.40186664953440_wp, & - & 5.83962338438617_wp, 10.13957405739921_wp, 8.57608868644902_wp, & - & -0.31075648460643_wp, 10.41940511029398_wp, 12.90511986423551_wp, & - & 0.36884840129967_wp, 7.95718616550666_wp, 13.74008757124138_wp, & - & 5.13899142512075_wp, 1.27780458134152_wp, 7.74094984443455_wp, & - & -0.43135113991668_wp, 7.08361972732341_wp, 16.04852770237523_wp, & - & 5.85948725126428_wp, 2.12703296273982_wp, 5.39640022648409_wp, & - & 7.98202471665911_wp, 9.37739667910611_wp, 4.56827791982260_wp, & - & 2.58521074631042_wp, 6.23620255015096_wp, 0.63970266218240_wp, & - & -7.37287492489702_wp, 6.83831475498402_wp, 13.72263180036320_wp, & - & -0.33592003395561_wp, 6.47879071615157_wp, 7.13632985901710_wp, & - & -9.65567668851609_wp, 2.93064163769937_wp, 14.45885460740137_wp, & - & 0.44488747834642_wp, 5.60046986687546_wp, 4.70073641648637_wp, & - &-10.35051844686700_wp, 3.77330440377895_wp, 16.93226888683481_wp, & - & -2.67740619557892_wp, 2.59670085855714_wp, 7.83370501910091_wp, & - & -8.03719513476919_wp, 7.73259684150050_wp, 16.06666801328784_wp, & - & -1.95562894546284_wp, 1.67604311050164_wp, 5.51722154256241_wp, & - & 1.91027980799140_wp, 7.09482655347090_wp, 3.06297438409324_wp, & - & -8.19632876596279_wp, 4.41220669444970_wp, 12.90802915938188_wp, & - & 8.06439765037426_wp, 1.37696801270714_wp, 1.32184680649995_wp, & - & 0.30124800389322_wp, 2.27996651956055_wp, 1.45618778825843_wp, & - & 10.17812550339424_wp, 8.63514254729744_wp, 0.55738672304122_wp, & - & 7.30609985362735_wp, 0.55954890552544_wp, 3.73485042789423_wp, & - & -1.86811945051845_wp, 5.02597852648671_wp, 8.64608290497032_wp, & - & -0.39194685290862_wp, 3.14311853344512_wp, 3.87021821970436_wp, & - & 0.34646325134892_wp, 11.08026825220325_wp, 11.06935462688116_wp, & - & 5.28103743231810_wp, 9.50927498680126_wp, 10.45446654094714_wp, & - & -2.26789230598577_wp, 2.53398791223687_wp, 13.76284852738644_wp, & - & 7.76197910212852_wp, 6.69602202299020_wp, 7.67831443128346_wp, & - & 9.83472737946160_wp, 5.89603575648236_wp, 3.55909477405225_wp, & - & 5.35897816853694_wp, 4.02359018773944_wp, 4.77723376533483_wp, & - & 7.54765344185304_wp, 3.26594082001099_wp, 0.68608024951558_wp, & - & -6.21783671028862_wp, 7.98605223285624_wp, 12.46513175710043_wp, & - & -3.91122464008637_wp, 1.48801747180495_wp, 9.05064606532323_wp, & - & -2.61184696118149_wp, 11.15215041991691_wp, 4.87991477050055_wp, & - & -0.35396280221428_wp, 0.43333061103608_wp, 0.82127690631710_wp, & - & 2.51885417779770_wp, 8.96364971376210_wp, 3.67683565997560_wp, & - & 5.60295691186475_wp, 1.08117572257201_wp, 13.82633961558057_wp, & - & 0.27625154557950_wp, 8.34365186719871_wp, 7.75926129035577_wp, & - & -7.67466304247654_wp, 3.74319313352409_wp, 11.03238946502150_wp, & - & -2.49374222683799_wp, 5.73280097562707_wp, 10.47517387698890_wp, & - & 1.55529102490321_wp, 6.77673381130491_wp, 12.54402499606942_wp, & - & 4.03617550083984_wp, 2.48938700817517_wp, 8.98510135702578_wp, & - & -7.43995626641283_wp, 9.60481600971518_wp, 16.67385502383442_wp, & - & 0.10707944729661_wp, 5.20506134431823_wp, 16.69319327980730_wp],& + & -3.66228842747680_wp, 0.96035860621450_wp, 8.20178752989411_wp, & + & -0.50644649048699_wp, 6.57586752834559_wp, 8.97052794903630_wp, & + & 7.42888971418823_wp, 10.27084821061875_wp, 8.96429185419075_wp, & + & 4.27323674976949_wp, 4.65515031591658_wp, 8.19574040761964_wp, & + & -0.89043875491593_wp, 3.60219514987424_wp, 5.52140058173462_wp, & + & -3.27829616304786_wp, 9.21770407200532_wp, 11.65091489719579_wp, & + & 4.65722901419843_wp, 7.62901166695902_wp, 11.64467880235024_wp, & + & 7.04489744975929_wp, 2.01350274482793_wp, 5.51516448688908_wp, & + & 1.42220756992726_wp, 4.25377257494854_wp, 1.50176502235078_wp, & + & -5.59094248789105_wp, 9.86928149707963_wp, 15.67055045657963_wp, & + & 2.34439371678417_wp, 6.97743424188471_wp, 15.66431436173409_wp, & + & 9.35773274717355_wp, 1.36173634718255_wp, 1.49571790007631_wp, & + & 12.09556735694052_wp, 6.36535208416555_wp, 2.47270609254561_wp, & + & -0.39344089298280_wp, 0.74984316203446_wp, 14.68732616926479_wp, & + & -8.32877709765801_wp, 4.86566576009662_wp, 14.69337329153926_wp, & + & 4.16004217969423_wp, 10.48117468222771_wp, 2.47894218739116_wp, & + & 10.71852423150072_wp, 7.94780839436631_wp, 6.63822847680108_wp, & + & 0.98360223245700_wp, 2.33229947223522_wp, 10.52161481243824_wp, & + & -6.95192294478929_wp, 3.28320944989586_wp, 10.52785090728379_wp, & + & 2.78318802682550_wp, 8.89871837202695_wp, 6.64446457164663_wp, & + & -3.08044188112999_wp, 0.28175810347616_wp, 6.34569893677352_wp, & + & -1.08848200940487_wp, 5.89726702560724_wp, 10.82642756958582_wp, & + & 6.84704316784142_wp, 10.94944871335710_wp, 10.82038044731135_wp, & + & 4.85508329611630_wp, 5.33375081865493_wp, 6.33965181449905_wp, & + & -1.53256755143637_wp, 1.75838977387366_wp, 4.86320911667262_wp, & + & -2.63616736652742_wp, 7.37389869600475_wp, 12.30910636225780_wp, & + & 5.29916883814780_wp, 9.47281704295959_wp, 12.30305923998332_wp, & + & 6.40295762580992_wp, 3.85711914825742_wp, 4.85697302182707_wp, & + & -0.75003213460554_wp, 0.92350895485444_wp, 2.41525843093813_wp, & + & -3.41889175592932_wp, 6.53920684955661_wp, 14.75705704799229_wp, & + & 4.51663342131697_wp, 10.30750888940773_wp, 14.75100992571782_wp, & + & 7.18549304264075_wp, 4.69199996727664_wp, 2.40902233609258_wp, & + & 0.79519657909300_wp, 2.40297521381811_wp, 0.84527431042847_wp, & + & -4.96393149705678_wp, 8.01867310852027_wp, 16.32704116850195_wp, & + & 2.97159368018951_wp, 8.82804263044407_wp, 16.32080507365639_wp, & + & 8.73053278376821_wp, 3.21253370831298_wp, 0.83922718815400_wp, & + & 14.30333390483585_wp, 9.69108036255379_wp, 1.53218960629421_wp, & + & -2.60120744087813_wp, 4.07557144042271_wp, 15.62784265551618_wp, & + &-10.53673261812442_wp, 1.53993748170838_wp, 15.63388977779065_wp, & + & 6.36799770016063_wp, 7.15544640383947_wp, 1.53823672856869_wp, & + & 12.71294074664985_wp, 8.21350382930090_wp, 3.13392111874485_wp, & + & -1.01081428269213_wp, 2.59799490716981_wp, 14.02592217049447_wp, & + & -8.94615048736734_wp, 3.01751401496127_wp, 14.03215826534002_wp, & + & 4.77741556940355_wp, 8.63302293709236_wp, 3.14015721359040_wp, & + & 11.95497176405906_wp, 9.10224183107713_wp, 5.45903963327914_wp, & + & -0.25265632753026_wp, 3.48673290894605_wp, 11.70099262853126_wp, & + & -8.18818150477655_wp, 2.12896498575612_wp, 11.70703975080573_wp, & + & 4.01944658681277_wp, 7.74447390788721_wp, 5.46508675555361_wp],& & shape(xyz)) real(wp), parameter :: lattice(3, 3) = reshape([& - & 15.90091151_wp, 0.00000000_wp, 0.00000000_wp, & - & 0.00000000_wp, 11.32002641_wp, 0.00000000_wp, & - &-12.11389533_wp, 0.00000000_wp, 17.11350086_wp],& + & 15.87093489225165_wp, 0.00000000000000_wp, -0.01227036698519_wp, & + & 0.00000000000000_wp, 11.23112234609398_wp, 0.00000000000000_wp, & + &-12.10426425260664_wp, 0.00000000000000_wp, 17.17840096263681_wp],& & shape(lattice)) call new(self, sym, xyz, lattice=lattice) end subroutine x05