Implement atomic masses

src/mctc/data.f90

@@ -22,9 +22,10 @@
 module mctc_data
    use mctc_data_atomicrad, only : get_atomic_rad
    use mctc_data_covrad, only : get_covalent_rad
+   use mctc_data_mass, only : get_atomic_mass
    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
+   public :: get_atomic_rad, get_covalent_rad, get_pauling_en, get_vdw_rad, get_atomic_mass
 end module mctc_data

src/mctc/data/CMakeLists.txt

@@ -18,6 +18,7 @@ list(
   APPEND srcs
   "${dir}/atomicrad.f90"
   "${dir}/covrad.f90"
+  "${dir}/mass.f90"
   "${dir}/paulingen.f90"
   "${dir}/vdwrad.f90"
 )

src/mctc/data/mass.f90

@@ -0,0 +1,108 @@
+! 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/mass.f90
+!> Provides atomic masses for all elements
+
+! atomic masses of the elements in isotop wighted form.
+
+! Source: 
+
+module mctc_data_mass
+    use mctc_env, only : wp
+    use mctc_io_convert, only : gmoltoau
+    use mctc_io_symbols, only : to_number
+    implicit none
+    private
+
+    public :: get_atomic_mass
+
+    !> Get atomic mass for a species
+    interface get_atomic_mass
+        module procedure :: get_atomic_mass_symbol
+        module procedure :: get_atomic_mass_number
+    end interface get_atomic_mass
+
+    integer, parameter :: max_elem = 118
+
+    !>  Isotope-averaged atomic masses in atomic units (in g/mol) from
+    !> https://github.com/grimme-lab/xtb/blob/main/src/param/atomicmass.f90
+    !> commit 7618f60
+    real(wp), parameter :: atomic_masses(max_elem) = gmoltoau * [ &
+      &   1.00794075_wp,   4.00260193_wp,   6.94003660_wp,   9.01218307_wp,&
+      &  10.81102805_wp,  12.01073590_wp,  14.00670321_wp,  15.99940492_wp,&
+      &  18.99840316_wp,  20.18004638_wp,  22.98976928_wp,  24.30505162_wp,&
+      &  26.98153853_wp,  28.08549871_wp,  30.97376200_wp,  32.06478741_wp,&
+      &  35.45293758_wp,  39.94779856_wp,  39.09830091_wp,  40.07802251_wp,&
+      &  44.95590828_wp,  47.86674496_wp,  50.94146504_wp,  51.99613176_wp,&
+      &  54.93804391_wp,  55.84514443_wp,  58.93319429_wp,  58.69334711_wp,&
+      &  63.54603995_wp,  65.37778253_wp,  69.72306607_wp,  72.62755016_wp,&
+      &  74.92159457_wp,  78.95938856_wp,  79.90352778_wp,  83.79800000_wp,&
+      &  85.46766360_wp,  87.61664447_wp,  88.90584030_wp,  91.22364160_wp,&
+      &  92.90637300_wp,  95.95978854_wp,  97.90721240_wp, 101.06494014_wp,&
+      & 102.90549800_wp, 106.41532751_wp, 107.86814963_wp, 112.41155782_wp,&
+      & 114.81808663_wp, 118.71011259_wp, 121.75978367_wp, 127.60312648_wp,&
+      & 126.90447190_wp, 131.29276145_wp, 132.90545196_wp, 137.32689163_wp,&
+      & 138.90546887_wp, 140.11573074_wp, 140.90765760_wp, 144.24159603_wp,&
+      & 144.91275590_wp, 150.36635571_wp, 151.96437813_wp, 157.25213065_wp,&
+      & 158.92535470_wp, 162.49947282_wp, 164.93032880_wp, 167.25908265_wp,&
+      & 168.93421790_wp, 173.05415017_wp, 174.96681496_wp, 178.48497872_wp,&
+      & 180.94787564_wp, 183.84177755_wp, 186.20670455_wp, 190.22485963_wp,&
+      & 192.21605165_wp, 195.08445686_wp, 196.96656879_wp, 200.59916703_wp,&
+      & 204.38341284_wp, 207.21690806_wp, 208.98039910_wp, 208.98243080_wp,&
+      & 209.98714790_wp, 222.01757820_wp, 223.01973600_wp, 226.02541030_wp,&
+      & 227.02775230_wp, 232.03805580_wp, 231.03588420_wp, 238.02891046_wp,&
+      & 237.04817360_wp, 244.06420530_wp, 243.06138130_wp, 247.07035410_wp,&
+      & 247.07030730_wp, 251.07958860_wp, 252.08298000_wp, 257.09510610_wp,&
+      & 258.09843150_wp, 259.10103000_wp, 262.10961000_wp, 267.12179000_wp,&
+      & 269.12791000_wp, 271.13393000_wp, 270.13336000_wp, 276.14846000_wp,&
+      & 276.15159000_wp, 280.16131000_wp, 282.16912000_wp, 284.17416000_wp,&
+      & 284.17873000_wp, 289.19042000_wp, 288.19274000_wp, 293.20449000_wp,&
+      & 292.20746000_wp, 294.21392000_wp]
+
+
+contains
+
+!> Get atomic mass for species with a given symbol
+elemental function get_atomic_mass_symbol(symbol) result(mass)
+
+   !> Element symbol
+   character(len=*), intent(in) :: symbol
+
+   !> atomic mass
+   real(wp) :: mass
+
+   mass = get_atomic_mass(to_number(symbol))
+
+end function get_atomic_mass_symbol
+
+
+!> Get atomic mass for species with a given atomic number
+elemental function get_atomic_mass_number(number) result(mass)
+
+   !> atomic number
+   integer, intent(in) :: number
+
+   !> atomic mass
+   real(wp) :: mass
+
+   if (number > 0 .and. number <= size(atomic_masses)) then
+      mass = atomic_masses(number)
+   else
+      mass = 0.0_wp
+   end if
+
+end function get_atomic_mass_number
+
+end module mctc_data_mass

src/mctc/data/meson.build

@@ -15,6 +15,7 @@
 srcs += files(
   'atomicrad.f90',
   'covrad.f90',
+  'mass.f90',
   'paulingen.f90',
   'vdwrad.f90'
 )

test/test_data.f90

@@ -16,9 +16,10 @@ module test_data
    use mctc_env, only : wp
    use mctc_env_testing, only : new_unittest, unittest_type, error_type, &
    & test_failed, check
+   use mctc_io_convert, only : gmoltoau
    use mctc_io_structure, only : structure_type
    use testsuite_structure, only : get_structure
-   use mctc_data, only : get_atomic_rad, get_covalent_rad, get_pauling_en, get_vdw_rad
+   use mctc_data, only : get_atomic_rad, get_covalent_rad, get_atomic_mass, get_pauling_en, get_vdw_rad
    implicit none
    private
 
@@ -44,7 +45,9 @@ subroutine collect_data(testsuite)
       & new_unittest("pauling_en", test_pauling_en), &
       & new_unittest("pauling_en_mb03", test_pauling_en_mb03), &
       & new_unittest("vdw_rad", test_vdw_rad), &
-      & new_unittest("vdw_rad_mb04", test_vdw_rad_mb04) &
+      & new_unittest("vdw_rad_mb04", test_vdw_rad_mb04), &
+      & new_unittest("atomic_mass", test_atomic_mass), &
+      & new_unittest("atomic_mass_mb05", test_atomic_mass_mb05) &
       & ]
 
    end subroutine collect_data
@@ -253,5 +256,49 @@ subroutine test_vdw_rad_mb04(error)
 
    end subroutine test_vdw_rad_mb04
 
+   subroutine test_atomic_mass(error)
+
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+
+      call check(error, get_atomic_mass("C"), get_atomic_mass(6))
+      if (allocated(error)) return
+      call check(error, get_atomic_mass("Am"), get_atomic_mass(95))
+      if (allocated(error)) return
+      call check(error, get_atomic_mass("Og"), get_atomic_mass(118))
+      if (allocated(error)) return
+      call check(error, get_atomic_mass("X"), get_atomic_mass(-1))
+
+   end subroutine test_atomic_mass  
+
+
+   subroutine test_atomic_mass_mb05(error)
+
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+
+      type(structure_type) :: mol
+      real(wp), allocatable :: mass_sym(:), mass_num(:)
+
+      real(wp), parameter :: ref(6) = gmoltoau * [&
+      & 1.08110280500000E+1_wp, 3.09737620000000E+1_wp, 1.0079407500000E+0_wp, &
+      & 3.54529375800000E+1_wp, 1.40067032100000E+1_wp, 2.80854987100000E+1_wp]
+
+      call get_structure(mol, "mindless05")
+
+      mass_sym = get_atomic_mass(mol%sym)
+      mass_num = get_atomic_mass(mol%num)
+
+      if (any(abs(mass_sym - mass_num) > thr) .or. any(abs(mass_sym - ref) > thr)) then
+         call test_failed(error, "Atomic masses do not match")
+         print'(3es21.14)', mass_sym
+         print'("---")'
+         print'(3es21.14)', mass_num
+         print'("---")'
+         print'(3es21.14)', ref
+      end if
+
+   end subroutine test_atomic_mass_mb05
+
 
 end module test_data