Source code for ccpn.util.isotopes

This file contains NMR isotope related functionalities
# Licence, Reference and Credits
__copyright__ = "Copyright (C) CCPN project ( 2014 - 2021"
__credits__ = ("Ed Brooksbank, Joanna Fox, Victoria A Higman, Luca Mureddu, Eliza Płoskoń",
               "Timothy J Ragan, Brian O Smith, Gary S Thompson & Geerten W Vuister")
__licence__ = ("CCPN licence. See")
__reference__ = ("Skinner, S.P., Fogh, R.H., Boucher, W., Ragan, T.J., Mureddu, L.G., & Vuister, G.W.",
                 "CcpNmr AnalysisAssign: a flexible platform for integrated NMR analysis",
                 "J.Biomol.Nmr (2016), 66, 111-124,")
# Last code modification
__modifiedBy__ = "$modifiedBy: Geerten Vuister $"
__dateModified__ = "$dateModified: 2021-12-23 11:27:19 +0000 (Thu, December 23, 2021) $"
__version__ = "$Revision: 3.0.4 $"
# Created
__author__ = "$Author: gvuister $"
__date__ = "$Date: 2017-04-07 10:28:48 +0000 (Fri, April 07, 2017) $"
# Start of code

import math
from collections import Counter, OrderedDict

from ccpn.util.traits.CcpNmrJson import CcpNmrJson, CcpnJsonDirectoryABC
from ccpn.util.traits.CcpNmrTraits import Unicode, Int, Float, Bool

from ccpn.framework.PathsAndUrls import ccpnConfigPath
from ccpn.util.Path import aPath
from ccpn.util.decorators import singleton

# Definitions for findMostLikelyFieldFromSpectrometerFrequencies
# The default isotopes checked:
defaultIsotopes = '1H 13C 15N 19F 31P'.split()
# The default fields (defined in MHz; checked in range Field.validRange):
fieldsInMHz = (300.0, 360.0, 400.0, 500.0, 600.0, 700.0, 750.0, 800.0, 850.0, 900.0, 950.0,
               1000.0, 1100.0, 1200.0)


[docs]class IsotopeRecord(CcpNmrJson): """Class to store isotope information """ classVersion = VERSION saveAllTraitsToJson = True isotopeCode = Unicode(allow_none=False, default_value='??').tag(info='The isotope code identifier, e.g. 1H, 13C, etc') elementNumber = Int(allow_none=False, default_value=0).tag(info='The element number') massNumber = Int(allow_none=False, default_value=0).tag(info='The mass number') isRadioactive = Bool(allow_none=True, default_value=None).tag(info='Flag identifying if isotope is radioactive') symbol = Unicode(allow_none=False, default_value='?').tag(info='The symbol identifying the isotope; eg. H for 1H, C for 13C, etc') name = Unicode(allow_none=True, default_value=None).tag(info='The name of the isotop; e.g. hydrogen') spin = Float(allow_none=False, default_value=0.0).tag(info='The spin of the isotope') gFactor = Float(allow_none=False, default_value=0.0).tag(info='The g-factor of the isotope') abundance = Float(allow_none=True, default_value=None).tag(info='The abundance of the isotope') quadrupoleMoment = Float(allow_none=True, default_value=None).tag(info='The quadrupoleMoment of the isotope') # to allow them to be sorted def __le__(self, other): return self.massNumber <= other.massNumber def __lt__(self, other): return self.massNumber < other.massNumber def __str__(self): return '<IsotopeRecord %s>' % self.isotopeCode def __repr__(self): return 'IsotopeRecord(%s)' % ', '.join(['%s=%r'%(k,v) for k,v in self.items()])
#end class # register the class for restoring IsotopeRecord.register()
[docs]def getIsotopeRecords(): """Coveniance funtions to return the istopeRecords Instance """ return IsotopeRecords()
[docs]@singleton class IsotopeRecords(CcpnJsonDirectoryABC): """Singleton class to contain all isotopeRecords as (isotopeCode, IsotopeRecord) (key,value) pairs """ attributeName = 'isotopeCode' directory = aPath(ccpnConfigPath) / 'isotopeRecords' sorted = True
[docs] def isotopesWithSpin(self, minValue=0.5, maxValue=None): """Return list of all isotopeRecords that have spin between minValue and maxValue (inclusive) """ if maxValue is None: maxValue = 1000.0 # Some rediculus high number result = [] for record in self.values(): if minValue <= record.spin <= maxValue: result.append(record) return result
[docs] def isotopesWithSpinHalf(self): return self.isotopesWithSpin(minValue=0.5, maxValue=0.5)
#end class # create an instance isotopeRecords = getIsotopeRecords()
[docs]class Nucleus(str): """ Class that behaves as a string, but holds the information about a particular nucleus and has usefull (!?) conversion methods """ muDivHbar = 7.62259328547 # MHz/T (from def __init__(self, isotopeCode): super().__init__() self.isotopeCode = isotopeCode self.isotopeRecord = isotopeRecords.get(self) @property def gamma(self): """return gamma in MHz/T; NB the isotopeRecord for 1H has the shielded g-value; most likely also for the other nuclei (if available) """ if not self.isotopeRecord: raise RuntimeError('Undefined isotopeRecord for %s' % self) gFac = self.isotopeRecord.gFactor return gFac * self.muDivHbar
[docs] def frequencyAtField(self, field): """Return absolute frequency (MHz) at field (T) """ if not self.isotopeRecord: raise RuntimeError('Undefined isotopeRecord for %s' % self) freq = self.gamma * field return math.fabs(freq)
[docs] def fieldAtFrequency(self, frequency): """Return field (T) at absolute frequency (MHz) """ if not self.isotopeRecord: raise RuntimeError('Undefined isotopeRecord for %s' % self) field = frequency / self.gamma return math.fabs(field)
@property def axisCode(self): """Return axis code string """ if not self.isotopeRecord: raise RuntimeError('Undefined isotopeRecord for %s' % self) return self.isotopeRecord.symbol def __str__(self): return '<Nucleus %s>' % super().__str__()
[docs]class Field(dict): """Class to hold the frequencies of various nuclei for a given field """ validRange = (-1.0, 0.5) # Variation in spectrometer frequences (MHz); accommodate Agilent/Varian and Bruker def __init__(self, field): super().__init__() self.field = field # add the default nuclei and their frequencies that will be tried for isotope in defaultIsotopes: self.addNucleus(isotope)
[docs] def addNucleus(self, isotopeCode): """Add a Nucleus instance defined by isotopeCode """ if isotopeCode not in isotopeRecords: raise ValueError('Undefined isotope "%s"' % isotopeCode) nuc = Nucleus(isotopeCode) self[nuc] = nuc.frequencyAtField(self.field)
@property def nuclei(self): """return a tuple of Nucleus instances for this field """ return tuple([n for n, f in self.items()]) @property def frequencies(self): """return a tuple of frequencies for this field """ return tuple([f for n, f in self.items()])
[docs] def isInRange(self, nucleus, freq): """return True in freq is withing validRange """ if nucleus not in self: raise ValueError('frequency of nucleus "%s" is not defined' % nucleus) lower, upper = (self[nucleus] + self.validRange[0], self[nucleus] + self.validRange[1]) #print('>', lower, freq, upper) return (lower <= freq <= upper)
[docs] def findNucleus(self, freq): """Find nucleus corresponding with freq; return None if not found """ for nucleus in self.nuclei: if self.isInRange(nucleus, freq): return nucleus return None
def __le__(self, other): return self.field <= other.field def __lt__(self, other): return self.field < other.field def __str__(self): return '<Field %.4fT>' % self.field def __repr__(self): return 'Field(%f)' % self.field
[docs]def findMostLikelyFieldFromSpectrometerFrequencies(spectrometerFrequencies): """Return the most likely field corresponding to spectrometerFrequencies, or None if it cannot be determined """ # list of fields to check field500 = 11.7433073413584267 # 500.0 MHz reference; tweaked in 14+ th digit to minimize rounding fieldList = [Field(freq / 500.0 * field500) for freq in fieldsInMHz] # list of Field instances fieldsFound = [] for sf in spectrometerFrequencies: for indx, field in enumerate(fieldList): nucleus = field.findNucleus(sf) if nucleus is not None: fieldsFound.append(indx) #print('Found>>', field.field, sf, nucleus) counter = Counter(fieldsFound) mostCommon = counter.most_common() #print('mostCommon', mostCommon) if len(mostCommon) > 0: return fieldList[mostCommon[0][0]] else: return None
[docs]def findNucleiFromSpectrometerFrequencies(spectrometerFrequencies): """return a tuple of Nuclei instances corresponding to spectrometerFrequencies or None if not identified """ theField = findMostLikelyFieldFromSpectrometerFrequencies(spectrometerFrequencies) if theField is None: return None nuclei = [theField.findNucleus(sf) for sf in spectrometerFrequencies] return tuple(nuclei)
[docs]def isotopeCode2Nucleus(isotopeCode): """:return the nucleus symbol for isotopeCode or None if not found """ if not isotopeCode: return None record = isotopeRecords.get(isotopeCode) if record is None: return None else: return record.symbol.upper()
#================================================================================================= # This routines were copied/adapted from util.Common #================================================================================================= # A translation dict DEFAULT_ISOTOPE_DICT = OrderedDict(( ('H', '1H'), ('D', '2H'), ('B', '11B'), ('C', '13C'), ('N', '15N'), ('O', '17O'), ('F', '19F'), ('P', '31P'), ('S', '33S'), ('K', '39K'), ('V', '51V'), ('Y', '89Y'), ('I', '127I'), ('W', '183W'), ('U', '235U'), ('HE', '3He'), ('LI', '7Li'), ('BE', '9Be'), ('NE', '21Ne'), ('NA', '23Na'), ('MG', '25Mg'), ('AL', '27Al'), ('SI', '29Si'), ('CL', '35Cl'), ('AR', '40Ar'), ('CA', '43Ca'), ('SC', '45Sc'), ('TI', '47Ti'), ('CR', '53Cr'), ('MN', '55Mn'), ('FE', '57Fe'), ('CO', '59Co'), ('NI', '61Ni'), ('CU', '63Cu'), ('ZN', '67Zn'), ('GA', '69Ga'), ('GE', '73Ge'), ('AS', '75As'), ('SE', '77Se'), ('BR', '79Br'), ('KR', '83Kr'), ('RB', '85Rb'), ('SR', '87Sr'), ('ZR', '91Zr'), ('NB', '93Nb'), ('MO', '95Mo'), ('TC', '99Tc'), ('RU', '99Ru'), ('RH', '103Rh'), ('PD', '105Pd'), ('AG', '107Ag'), ('CD', '111Cd'), ('IN', '115In'), ('SN', '119Sn'), ('SB', '121Sb'), ('TE', '125Te'), ('XE', '129Xe'), ('CS', '133Cs'), ('BA', '137Ba'), ('LA', '139La'), ('CE', '140Ce'), ('PR', '141Pr'), ('ND', '143Nd'), ('PM', '147Pm'), ('SM', '144Sm'), ('EU', '153Eu'), ('GD', '157Gd'), ('TB', '159Tb'), ('DY', '163Dy'), ('HO', '165Ho'), ('ER', '167Er'), ('TM', '169Tm'), ('YB', '171Yb'), ('LU', '175Lu'), ('HF', '177Hf'), ('TA', '181Ta'), ('RE', '187Re'), ('OS', '187Os'), ('IR', '193Ir'), ('PT', '195Pt'), ('AU', '197Au'), ('HG', '199Hg'), ('TL', '205Tl'), ('PB', '207Pb'), ('BI', '209Bi'), ('PO', '209Po'), ('AC', '227Ac'), ('TH', '232Th'), ('NP', '237Np'), ('PU', '239Pu'), ('AM', '243Am'), ))
[docs]def checkIsotope(isotopeCode) -> str: """Convert isotopeCode string to most probable isotope code - defaulting to '1H' This function is intended for external format isotope specifications; mostly in SpectrumSource parameter reading routines """ defaultIsotope = '1H' if not isotopeCode: raise ValueError('Invalid isotopeCode (%r)' % isotopeCode) isotopeCode = isotopeCode.strip().upper() if isotopeCode in isotopeRecords: # Superfluous but should speed things up return isotopeCode for ic in isotopeRecords.keys(): # NB checking this first means that e.g. 'H13C' returns '13C' rather than '1H' if ic.upper() in isotopeCode: return ic # NB order of checking means that e.g. 'CA' returns Calcium rather than Carbon result = (DEFAULT_ISOTOPE_DICT.get(isotopeCode[:2]) or DEFAULT_ISOTOPE_DICT.get(isotopeCode[0])) if result is None: result = defaultIsotope return result
[docs]def name2IsotopeCode(name=None): """Get standard isotope code matching atom name or axisCode string """ if not name: return None result = DEFAULT_ISOTOPE_DICT.get(name[0]) if result is None: if name[0].isdigit(): ss = name.title() for key in isotopeRecords: if ss.startswith(key): if name[:len(key)].isupper(): result = key break else: result = DEFAULT_ISOTOPE_DICT.get(name[:2]) # return result
[docs]def name2ElementSymbol(name): """Get standard element symbol matching name or axisCode NB, the first letter takes precedence, so e.g. 'CD' returns 'C' (carbon) rather than 'CD' (Cadmium)""" # NB, We deliberately do NOT use 'value in Constants.DEFAULT_ISOTOPE_DICT' # We want to avoid elements that are in the dict but have value None. if not name: return None elif DEFAULT_ISOTOPE_DICT.get(name[0]) is not None: return name[0] elif DEFAULT_ISOTOPE_DICT.get(name[:2]) is not None: return name[:2] elif name[0].isdigit(): ss = name.title() for key, record in isotopeRecords.items(): if ss.startswith(key): if name[:len(key)].isupper(): return record.symbol.upper() break # return None