1"""Contains fundamental constants in atomic units.
2
3Copyright (C) 2013, Joshua More and Michele Ceriotti
4
5This program is free software: you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation, either version 3 of the License, or
8(at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13GNU General Public License for more details.
14
15You should have received a copy of the GNU General Public License
16along with this program. If not, see <http.//www.gnu.org/licenses/>.
17
18
19Classes:
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
26import re
27from ipi.utils.messages import verbosity, info
28
29__all__ = ['Constants', 'Elements', 'unit_to_internal', 'unit_to_user']
30
31
32class 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
46class 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
202UnitMap = {
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
282UnitPrefix = {
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
291UnitPrefixRE = ""
292for key in UnitPrefix:
293   UnitPrefixRE = UnitPrefixRE + key + "|"
294UnitPrefixRE = " *(" + UnitPrefixRE[1:] + ")(.*) *"
295UnitPrefixRE = 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
303def 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
346def 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)
359