tools/i-pi/ipi/utils/units.py

   1 """Contains fundamental constants in atomic units.
   2
   3 Copyright (C) 2013, Joshua More and Michele Ceriotti
   4
   5 This program is free software: you can redistribute it and/or modify
   6 it under the terms of the GNU General Public License as published by
   7 the Free Software Foundation, either version 3 of the License, or
   8 (at your option) any later version.
   9
  10 This program is distributed in the hope that it will be useful,
  11 but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  13 GNU General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License
  16 along with this program. If not, see <http.//www.gnu.org/licenses/>.
  17
  18
  19 Classes:
  20    Constants: Class whose members are fundamental constants.
  21    Elements: Class which contains the mass of different elements
  22    Units: Class which contains the methods needed to transform
  23       between different systems of units.
  24 """
  25
  26 import re
  27 from ipi.utils.messages import verbosity, info
  28
  29 __all__ = ['Constants', 'Elements', 'unit_to_internal', 'unit_to_user']
  30
  31
  32 class Constants:
  33    """Class whose members are fundamental constants.
  34
  35    Attributes:
  36       kb: Boltzmann constant.
  37       hbar: Reduced Planck's constant.
  38       amu: Atomic mass unit.
  39    """
  40
  41    kb = 1.0
  42    hbar = 1.0
  43    amu = 1822.8885
  44
  45
  46 class Elements(dict):
  47    """Class which contains the mass of different elements.
  48
  49    Attributes:
  50       mass_list: A dictionary containing the masses of different elements.
  51          Has the form {"label": Mass in a.m.u.}. Note that the generic "X"
  52          label is assumed to be an electron.
  53    """
  54
  55    mass_list={
  56       "X"   :    1.0000/Constants.amu,
  57       "H"   :   1.00794,
  58       "D"   :    2.0141,
  59       "Z"   :  1.382943, #an interpolated H-D atom, based on y=1/sqrt(m) scaling
  60       "H2"  :    2.0160,
  61       "He"  :  4.002602,
  62       "Li"  :    6.9410,
  63       "Be"  :  9.012182,
  64       "B"   :    10.811,
  65       "C"   :   12.0107,
  66       "N"   :  14.00674,
  67       "O"   :   15.9994,
  68       "F"   : 18.998403,
  69       "Ne"  :   20.1797,
  70       "Na"  : 22.989770,
  71       "Mg"  :   24.3050,
  72       "Al"  : 26.981538,
  73       "Si"  :   28.0855,
  74       "P"   : 30.973761,
  75       "S"   :    32.066,
  76       "Cl"  :   35.4527,
  77       "Ar"  :   39.9480,
  78       "K"   :   39.0983,
  79       "Ca"  :    40.078,
  80       "Sc"  : 44.955910,
  81       "Ti"  :    47.867,
  82       "V"   :   50.9415,
  83       "Cr"  :   51.9961,
  84       "Mn"  : 54.938049,
  85       "Fe"  :    55.845,
  86       "Co"  : 58.933200,
  87       "Ni"  :   58.6934,
  88       "Cu"  :    63.546,
  89       "Zn"  :     65.39,
  90       "Ga"  :    69.723,
  91       "Ge"  :     72.61,
  92       "As"  :  74.92160,
  93       "Se"  :     78.96,
  94       "Br"  :    79.904,
  95       "Kr"  :     83.80,
  96       "Rb"  :   85.4678,
  97       "Sr"  :     87.62,
  98       "Y"   :  88.90585,
  99       "Zr"  :    91.224,
 100       "Nb"  :  92.90638,
 101       "Mo"  :     95.94,
 102       "Tc"  :        98,
 103       "Ru"  :    101.07,
 104       "Rh"  : 102.90550,
 105       "Pd"  :    106.42,
 106       "Ag"  :  107.8682,
 107       "Cd"  :   112.411,
 108       "In"  :   114.818,
 109       "Sn"  :   118.710,
 110       "Sb"  :   121.760,
 111       "Te"  :    127.60,
 112       "I"   : 126.90447,
 113       "Xe"  :    131.29,
 114       "Cs"  : 132.90545,
 115       "Ba"  :   137.327,
 116       "La"  :  138.9055,
 117       "Ce"  :   140.166,
 118       "Pr"  : 140.90765,
 119       "Nd"  :    144.24,
 120       "Pm"  :       145,
 121       "Sm"  :    150.36,
 122       "Eu"  :   151.964,
 123       "Gd"  :    157.25,
 124       "Tb"  : 158.92534,
 125       "Dy"  :    162.50,
 126       "Ho"  : 164.93032,
 127       "Er"  :    167.26,
 128       "Tm"  : 168.93241,
 129       "Yb"  :    173.04,
 130       "Lu"  :   174.967,
 131       "Hf"  :    178.49,
 132       "Ta"  :  180.9479,
 133       "W"   :    183.84,
 134       "Re"  :   186.207,
 135       "Os"  :    190.23,
 136       "Ir"  :   192.217,
 137       "Pt"  :   195.078,
 138       "Au"  : 196.96655,
 139       "Hg"  :    200.59,
 140       "Tl"  :  204.3833,
 141       "Pb"  :     207.2,
 142       "Bi"  : 208.98038,
 143       "Po"  :       209,
 144       "At"  :       210,
 145       "Rn"  :       222,
 146       "Fr"  :       223,
 147       "Ra"  :       226,
 148       "Ac"  :       227,
 149       "Th"  :  232.0381,
 150       "Pa"  : 231.03588,
 151       "U"   :  238.0289,
 152       "Np"  :       237,
 153       "Pu"  :       244,
 154       "Am"  :       243,
 155       "Cm"  :       247,
 156       "Bk"  :       247,
 157       "Cf"  :       251,
 158       "Es"  :       252,
 159       "Fm"  :       257,
 160       "Md"  :       258,
 161       "No"  :       259,
 162       "Lr"  :       262,
 163       "Rf"  :       267,
 164       "Db"  :       268,
 165       "Sg"  :       269,
 166       "Bh"  :       270,
 167       "Hs"  :       269,
 168       "Mt"  :       278,
 169       "Ds"  :       281,
 170       "Rg"  :       280,
 171       "Cn"  :       285,
 172       "Uut" :       286,
 173       "Fl"  :       289,
 174       "Uup" :       288,
 175       "Lv"  :       293,
 176       "Uus" :       294,
 177       "Uuo" :       294
 178    }
 179
 180    @classmethod
 181    def mass(cls, label):
 182       """Function to access the mass_list attribute.
 183
 184       Note that this does not require an instance of the Elements class to be
 185       created, as this is a class method. Therefore using Elements.mass(label)
 186       will give the mass of the element with the atomic symbol given by label.
 187
 188       Args:
 189          label: The atomic symbol of the atom whose mass is required.
 190
 191       Returns:
 192          A float giving the mass of the atom with atomic symbol label.
 193       """
 194
 195       try:
 196          return cls.mass_list[label]*Constants.amu
 197       except KeyError:
 198          info("Unknown element given, you must specify the mass", verbosity.low)
 199          return -1.0
 200
 201 # these are the conversion FROM the unit stated to internal (atomic) units
 202 UnitMap = {
 203    "undefined": {
 204       ""             : 1.00
 205       },
 206    "energy":   {
 207       ""             : 1.00,
 208       "atomic_unit"  : 1.00,
 209       "electronvolt" : 0.036749326,
 210       "j/mol"        : 0.00000038087989,
 211       "cal/mol"      : 0.0000015946679,
 212       "kelvin"       : 3.1668152e-06
 213       },
 214    "temperature":   {
 215       ""             : 1.00,
 216       "atomic_unit"  : 1.00,
 217       "kelvin"       : 3.1668152e-06
 218       },
 219    "time":     {
 220       ""             : 1.00,
 221       "atomic_unit"  : 1.00,
 222       "second"       : 4.1341373e+16
 223       },
 224    "frequency" :   {   # NB Internally, ANGULAR frequencies are used.
 225       ""             : 1.00,
 226       "atomic_unit"  : 1.00,
 227       "inversecm"    : 4.5563353e-06,
 228       "hertz*rad"    : 2.4188843e-17,
 229       "hertz"        : 1.5198298e-16
 230       },
 231    "ms-momentum" :   {   # TODO fill up units here (mass-scaled momentum)
 232       ""             : 1.00,
 233       "atomic_unit"  : 1.00
 234       },
 235    "length" :     {
 236       ""             : 1.00,
 237       "atomic_unit"  : 1.00,
 238       "angstrom"     : 1.8897261,
 239       "meter"        : 1.8897261e+10
 240       },
 241    "volume" :     {
 242       ""             : 1.00,
 243       "atomic_unit"  : 1.00,
 244       "angstrom3"    : 6.748334231,
 245       },
 246    "velocity":    {
 247       ""            : 1.00,
 248       "atomic_unit" : 1.00,
 249       "m/s"         : 4.5710289e-7
 250       },
 251    "momentum":    {
 252       ""             : 1.00,
 253       "atomic_unit"  : 1.00
 254       },
 255    "mass":        {
 256       ""             : 1.00,
 257       "atomic_unit"  : 1.00,
 258       "dalton"       : 1.00*Constants.amu,
 259       "electronmass" : 1.00
 260       },
 261    "pressure" :     {
 262       ""             : 1.00,
 263       "atomic_unit"  : 1.00,
 264       "bar"          : 3.398827377e-9,
 265       "atmosphere"   : 3.44386184e-9,
 266       "pascal"       : 3.398827377e-14
 267       },
 268    "density" : {
 269       ""             : 1.00,
 270       "atomic_unit"  : 1.00,
 271       "g/cm3"        : 162.67263
 272       },
 273     "force" : {
 274       ""             : 1.00,
 275       "atomic_unit"  : 1.00,
 276       "newton"       : 12137805
 277       }
 278
 279 }
 280
 281 # a list of magnitude prefixes
 282 UnitPrefix = {
 283    "" : 1.0,
 284    "yotta" : 1e24, "zetta" : 1e21, "exa" : 1e18, "peta" : 1e15,
 285    "tera" : 1e12, "giga" : 1e9, "mega" : 1e6, "kilo" : 1e3,
 286    "milli" : 1e-3, "micro" : 1e-6, "nano" : 1e-9, "pico" : 1e-12,
 287    "femto" : 1e-15, "atto" : 1e-18, "zepto" : 1e-21, "yocto" : 1e-24
 288 }
 289
 290 # builds a RE to match prefix and split out the base unit
 291 UnitPrefixRE = ""
 292 for key in UnitPrefix:
 293    UnitPrefixRE = UnitPrefixRE + key + "|"
 294 UnitPrefixRE = " *(" + UnitPrefixRE[1:] + ")(.*) *"
 295 UnitPrefixRE = re.compile(UnitPrefixRE)
 296
 297 ########################################################################
 298 #  Atomic units are used EVERYWHERE internally. In order to quickly    #
 299 #  interface with any "outside" unit, we set up a simple conversion    #
 300 #  library.                                                            #
 301 ########################################################################
 302
 303 def unit_to_internal(family, unit, number):
 304    """Converts a number of given dimensions and units into internal units.
 305
 306    Args:
 307       family: The dimensionality of the number.
 308       unit: The units 'number' is originally in.
 309       number: The value of the parameter in the units 'unit'.
 310
 311    Returns:
 312       The number in internal units.
 313
 314    Raises:
 315       ValueError: Raised if the user specified units aren't given in the
 316          UnitMap dictionary.
 317       IndexError: Raised if the programmer specified dimensionality for the
 318          parameter isn't in UnitMap. Shouldn't happen, for obvious reasons.
 319       TypeError: Raised if the prefix is correct, but the base unit is not, in
 320          the user specified unit string.
 321    """
 322
 323    if not (family == "number" or family in UnitMap):
 324       raise IndexError(family + " is an undefined units kind.")
 325    if family == "number":
 326       return number
 327
 328
 329    if unit == "":
 330       prefix = ""
 331       base = ""
 332    else:
 333       m = UnitPrefixRE.match(unit);
 334       if m is None:
 335          raise ValueError("Unit " + unit + " is not structured with a prefix+base syntax.")
 336       prefix = m.group(1)
 337       base = m.group(2)
 338
 339    if not prefix in UnitPrefix:
 340       raise TypeError(prefix + " is not a valid unit prefix.")
 341    if not base in UnitMap[family]:
 342       raise TypeError(base + " is an undefined unit for kind " + family + ".")
 343
 344    return number*UnitMap[family][base]*UnitPrefix[prefix]
 345
 346 def unit_to_user(family, unit, number):
 347    """Converts a number of given dimensions from internal to user units.
 348
 349    Args:
 350       family: The dimensionality of the number.
 351       unit: The units 'number' should be changed to.
 352       number: The value of the parameter in internal units.
 353
 354    Returns:
 355       The number in the user specified units
 356    """
 357
 358    return number/unit_to_internal(family, unit, 1.0)