1# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding: utf-8 -*- 2# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 3# 4# MDAnalysis --- https://www.mdanalysis.org 5# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors 6# (see the file AUTHORS for the full list of names) 7# 8# Released under the GNU Public Licence, v2 or any higher version 9# 10# Please cite your use of MDAnalysis in published work: 11# 12# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, 13# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. 14# MDAnalysis: A Python package for the rapid analysis of molecular dynamics 15# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th 16# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. 17# 18# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. 19# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. 20# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 21# 22 23r""" 24Constants and unit conversion --- :mod:`MDAnalysis.units` 25=============================================================== 26 27The base units of MDAnalysis trajectories are the **Å** (**ångström**) for 28**length** and **ps** (**pico second**) for **time**. By default, all positions 29are in Å and all times are in ps, regardless of how the MD code stored 30trajectory data. By default, MDAnalysis converts automatically to the 31MDAnalysis units when reading trajectories and converts back when writing. This 32makes it possible to write scripts that can be agnostic of the specifics of how 33a particular MD code stores trajectory data. Other base units are listed in the 34table on :ref:`table-baseunits`. 35 36.. _table-baseunits: 37 38.. Table:: Base units in MDAnalysis 39 40 =========== ============== =============================================== 41 quantity unit SI units 42 =========== ============== =============================================== 43 length Å :math:`10^{-10}` m 44 time ps :math:`10^{-12}` s 45 energy kJ/mol :math:`1.66053892103219 \times 10^{-21}` J 46 charge :math:`e` :math:`1.602176565 \times 10^{-19}` As 47 force kJ/(mol·Å) :math:`1.66053892103219 \times 10^{-11}` J/m 48 speed Å/ps :math:`100` m/s 49 =========== ============== =============================================== 50 51Implementation notes 52-------------------- 53 54All conversions with :func:`convert` are carried out in a simple fashion: the 55conversion factor :math:`f_{b,b'}` from the base unit :math:`b` to another unit 56:math:`b'` is precomputed and stored (see :ref:`Data`). The numerical value of 57a quantity in unit :math:`b` is :math:`X/b` (e.g. for :math:`X = 581.23\,\mathrm{ps}` the numerical value is :math:`X/\mathrm{ps} = 591.23`). [#funits]_ 60 61The new numerical value :math:`X'/b'` of the quantity (in units of :math:`b'`) 62is then 63 64.. math:: 65 66 X'/b' = f_{b,b'} X/b 67 68The function :func:`get_conversion_factor` returns the appropriate factor 69:math:`f_{b,b'}`. 70 71Conversion between different units is always carried out via the base unit as 72an intermediate step:: 73 74 x is in u1: from u1 to b: x' = x / factor[u1] 75 from b to u2: x'' = x' * factor[u2] 76 so f[u1,u2] = factor[u2]/factor[u1] 77 78 79Conversions 80~~~~~~~~~~~ 81 82Examples for how to calculate some of the conversion factors that are 83hard-coded in :mod:`~MDAnalysis.units` (see :ref:`Data`). 84 85density: 86 Base unit is :math:`\mathrm{Å}^{-3}`:: 87 88 n/x = n/A**3 * densityUnit_factor[x] 89 90 Example for how to calculate the conversion factor 91 :math:`f_{\mathrm{Å}^{-3},\mathrm{nm}^{-3}}` from :math:`\mathrm{Å^-3}` to 92 :math:`\mathrm{nm}^{-3}`: 93 94 .. math:: 95 96 f_{\mathrm{Å}^{-3},\mathrm{nm}^{-3}} 97 = \frac{1\,\mathrm{nm}^{-3}}{1\,\mathrm{Å}^{-3}} 98 = \frac{(10\,\mathrm{Å})^{-3}}{1\,\mathrm{Å}^{-3}} 99 = 10^{-3} 100 101concentration: 102 Example for how to convert the conversion factor to Molar (mol/l):: 103 104 factor = 1 A**-3 / (N_Avogadro * (10**-9 dm)**-3) 105 106 relative to a density rho0 in g/cm^3:: 107 108 M(H2O) = 18 g/mol Molar mass of water 109 110 factor = 1/(1e-24 * N_Avogadro / M(H2O)) 111 112 from `rho/rho0 = n/(N_A * M**-1) / rho0` where `[n] = 1/Volume`, `[rho] = mass/Volume` 113 114 115Note 116---- 117In the future me might move towards using the Quantities_ package or 118:mod:`scipy.constants`. 119 120 121.. _Quantities: http://packages.python.org/quantities/ 122 123Functions 124--------- 125 126.. autofunction:: get_conversion_factor 127.. autofunction:: convert 128 129.. _Data: 130 131Data 132---- 133 134.. autodata:: constants 135.. autodata:: lengthUnit_factor 136.. autodata:: N_Avogadro 137.. autodata:: water 138.. autodata:: densityUnit_factor 139.. autodata:: timeUnit_factor 140.. autodata:: speedUnit_factor 141.. autodata:: forceUnit_factor 142.. autodata:: chargeUnit_factor 143.. autodata:: conversion_factor 144.. autodata:: unit_types 145 146 147References and footnotes 148------------------------ 149 150.. [Jorgensen1998] W. Jorgensen, C. Jenson, J Comp Chem 19 (1998), 1179-1186 151 152.. _AKMA: http://www.charmm.org/documentation/c37b1/usage.html#%20AKMA 153.. _electron charge: http://physics.nist.gov/cgi-bin/cuu/Value?e 154.. _`Avogadro's constant`: http://physics.nist.gov/cgi-bin/cuu/Value?na 155 156.. Rubric:: Footnotes 157 158.. [#funits] One can also consider the conversion factor to carry 159 units :math:`b'/b`, in which case the conversion formula would 160 become 161 162 .. math:: 163 164 X' = f_{b,b'} X 165 166""" 167 168from __future__ import unicode_literals, division 169 170#: `Avogadro's constant`_ in mol**-1. 171#: 172#: .. deprecated:: 0.9.0 173#: Use *N_Avogadro* in dict :data:`constants` instead. 174#: 175#: .. versionchanged:: 0.9.0 176#: Use CODATA 2010 value, which differs from the previously used value 177#: 6.02214179e+23 mol**-1 by -5.00000000e+16 mol**-1. 178N_Avogadro = 6.02214129e+23 # mol**-1 179 180# 181# NOTE: Whenever a constant is added to the constants dict, you also 182# MUST add an appropriate entry to 183# test_units:TestConstants.constants_reference ! 184 185#: Values of physical constants are taken from `CODATA 2010 at NIST`_. The 186#: thermochemical calorie is defined in the `ISO 80000-5:2007`_ standard 187#: and is also listed in the `NIST Guide to SI: Appendix B.8: Factors for Units`_. 188#: 189#: .. _`CODATA 2010 at NIST`: 190#: http://physics.nist.gov/cuu/Constants/ 191#: .. _`ISO 80000-5:2007`: 192#: http://www.iso.org/iso/catalogue_detail?csnumber=31890 193#: .. _`NIST Guide to SI: Appendix B.8: Factors for Units`: 194#: http://physics.nist.gov/Pubs/SP811/appenB8.html#C 195#: 196#: .. versionadded:: 0.9.0 197constants = { 198 'N_Avogadro': 6.02214129e+23, # mol**-1 199 'elementary_charge': 1.602176565e-19, # As 200 'calorie': 4.184, # J 201} 202 203#: The basic unit of *length* in MDAnalysis is the Angstrom. 204#: Conversion factors between the base unit and other lengthUnits *x* are stored. 205#: Conversions follow `L/x = L/Angstrom * lengthUnit_factor[x]`. 206#: *x* can be *nm*/*nanometer* or *fm*. 207lengthUnit_factor = { 208 'Angstrom': 1.0, 'A': 1.0, 'angstrom': 1.0, 209 u'\u212b': 1.0, # Unicode and UTF-8 encoded symbol for angstroms 210 'nm': 1.0 / 10, 'nanometer': 1.0 / 10, 211 'pm': 1e2, 'picometer': 1e2, 212 'fm': 1e5, 'femtometer': 1e5, 213} 214 215 216#: water density values at T=298K, P=1atm [Jorgensen1998]_ 217#: ======== ========= 218#: model g cm**-3 219#: ======== ========= 220#: SPC 0.985(1) 221#: TIP3P 1.002(1) 222#: TIP4P 1.001(1) 223#: exp 0.997 224#: ======== ========= 225#: 226#: and molar mass 18.016 g mol**-1. 227water = { 228 'exp': 0.997, 'SPC': 0.985, 'TIP3P': 1.002, 'TIP4P': 1.001, # in g cm**-3 229 'MolarMass': 18.016, # in g mol**-1 230} 231 232#: The basic unit for *densities* is Angstroem**(-3), i.e. 233#: the volume per molecule in A**3. Especially for water 234#: it can be convenient to measure the density relative to bulk, and 235#: hence a number of values are pre-stored in :data:`water`. 236densityUnit_factor = { 237 'Angstrom^{-3}': 1 / 1.0, 'A^{-3}': 1 / 1.0, 238 '\u212b^{-3}': 1 / 1.0, 239 'nm^{-3}': 1 / 1e-3, 'nanometer^{-3}': 1 / 1e-3, 240 'Molar': 1 / (1e-27 * constants['N_Avogadro']), 241 'SPC': 1 / (1e-24 * constants['N_Avogadro'] * water['SPC'] / water['MolarMass']), 242 'TIP3P': 1 / (1e-24 * constants['N_Avogadro'] * water['TIP3P'] / water['MolarMass']), 243 'TIP4P': 1 / (1e-24 * constants['N_Avogadro'] * water['TIP4P'] / water['MolarMass']), 244 'water': 1 / (1e-24 * constants['N_Avogadro'] * water['exp'] / water['MolarMass']), 245} 246 247 248#: For *time*, the basic unit is ps; in particular CHARMM's 249#: 1 AKMA_ time unit = 4.888821E-14 sec is supported. 250timeUnit_factor = { 251 'ps': 1.0, 'picosecond': 1.0, # 1/1.0 252 'fs': 1e3, 'femtosecond': 1e3, # 1/1e-3, 253 'ns': 1e-3, 'nanosecond': 1e-3, # 1/1e3, 254 'second': 1e-12, 'sec': 1e-12, 's': 1e-12, # 1/1e12, 255 'AKMA': 1 / 4.888821e-2, 256} 257# getting the factor f: 1200ps * f = 1.2 ns ==> f = 1/1000 ns/ps 258 259#: For *speed*, the basic unit is Angstrom/ps. 260speedUnit_factor = { 261 'Angstrom/ps': 1.0, 'A/ps': 1.0, '\u212b/ps': 1.0, 262 'Angstrom/picosecond': 1.0, 263 'angstrom/picosecond': 1.0, # 1 264 'Angstrom/fs': 1.0 * 1e3, 265 'Angstrom/femtosecond': 1.0 * 1e3, 266 'angstrom/femtosecond': 1.0 * 1e3, 267 'angstrom/fs': 1.0 * 1e3, 268 'Angstrom/AKMA': 4.888821e-2, 269 'nm/ps': 0.1, 'nanometer/ps': 0.1, 'nanometer/picosecond': 0.1, # 1/10 270 'nm/ns': 0.1 / 1e-3, 271 'pm/ps': 1e2, 272 'm/s': 1e-10 / 1e-12, 273} 274# (TODO: build this combinatorically from lengthUnit and timeUnit) 275 276#: *Energy* is measured in kJ/mol. 277energyUnit_factor = { 278 'kJ/mol': 1.0, 279 'kcal/mol': 1/constants['calorie'], 280 'J': 1e3/constants['N_Avogadro'], 281 'eV': 1e3/(constants['N_Avogadro'] * constants['elementary_charge']), 282 } 283 284#: For *force* the basic unit is kJ/(mol*Angstrom). 285forceUnit_factor = { 286 'kJ/(mol*Angstrom)': 1.0, 'kJ/(mol*A)': 1.0, 287 'kJ/(mol*\u212b)': 1.0, 288 'kJ/(mol*nm)': 10.0, 289 'Newton': 1e13/constants['N_Avogadro'], 290 'N': 1e13/constants['N_Avogadro'], 291 'J/m': 1e13/constants['N_Avogadro'], 292 'kcal/(mol*Angstrom)': 1/constants['calorie'], 293} 294# (TODO: build this combinatorically from lengthUnit and energyUnit) 295 296#: *Charge* is measured in multiples of the `electron charge`_ *e*, with the value 297#: *elementary_charge* in :data:`constants`. 298#: The `conversion factor to Amber charge units`_ is 18.2223. 299#: 300#: .. _`conversion factor to Amber charge units`: http://ambermd.org/formats.html#parm 301#: 302#: .. versionchanged:: 0.9.0 303#: Use CODATA 2010 value for *elementary charge*, which differs from the previously used value 304#: *e* = 1.602176487 x 10**(-19) C by 7.8000000e-27 C. 305chargeUnit_factor = { 306 'e': 1.0, 307 'Amber': 18.2223, # http://ambermd.org/formats.html#parm 308 'C': constants['elementary_charge'], 'As': constants['elementary_charge'], 309} 310 311#: :data:`conversion_factor` is used by :func:`get_conversion_factor`: 312#: Note: any observable with a unit (i.e. one with an entry in 313#: the :attr:`unit` attribute) needs an entry in :data:`conversion_factor` 314conversion_factor = { 315 'length': lengthUnit_factor, 316 'density': densityUnit_factor, 317 'time': timeUnit_factor, 318 'charge': chargeUnit_factor, 319 'speed': speedUnit_factor, 320 'force': forceUnit_factor, 321 'energy': energyUnit_factor, 322} 323 324#: Generated lookup table (dict): returns the type of unit for a known input unit. 325#: Note: Any unit must be *unique* because this dict is used to guess the 326#: unit type. 327unit_types = {} 328for utype, ufactor in conversion_factor.items(): 329 for unit in ufactor.keys(): 330 assert not unit in unit_types # see comment! 331 unit_types[unit] = utype 332 333 334def get_conversion_factor(unit_type, u1, u2): 335 """generate the conversion factor u1 -> u2 by using the base unit as an intermediate 336 337 f[u1 -> u2] = factor[u2]/factor[u1] 338 339 Conversion of X (in u1) to X' (in u2): 340 341 X' = conversion_factor * X 342 """ 343 # x is in u1: from u1 to b: x' = x / factor[u1] 344 # from b to u2: x'' = x' * factor[u2] 345 # so f[u1,u2] = factor[u2]/factor[u1] 346 return conversion_factor[unit_type][u2] / conversion_factor[unit_type][u1] 347 348 349def convert(x, u1, u2): 350 """Convert value *x* in unit *u1* to new value in *u2*. 351 352 Returns 353 ------- 354 float 355 Converted value. 356 357 Raises 358 ------ 359 ValueError 360 The units are not known or if one attempts to convert between 361 incompatible units. 362 """ 363 try: 364 ut1 = unit_types[u1] 365 except KeyError: 366 raise ValueError("unit '{0}' not recognized.\n" 367 "It must be one of {1}.".format(u1, ", ".join(unit_types))) 368 try: 369 ut2 = unit_types[u2] 370 except KeyError: 371 raise ValueError("unit '{0}' not recognized.\n" 372 "It must be one of {1}.".format(u2, ", ".join(unit_types))) 373 if ut1 != ut2: 374 raise ValueError("Cannot convert between unit types " 375 "{0} --> {1}".format(u1, u2)) 376 return x * get_conversion_factor(ut1, u1, u2) 377