1# -*- coding: utf-8 -*- 2# 3# Copyright (c) 2020, the cclib development team 4# 5# This file is part of cclib (http://cclib.github.io) and is distributed under 6# the terms of the BSD 3-Clause License. 7 8"""Calculate properties of nuclei based on data parsed by cclib.""" 9 10import logging 11 12import numpy as np 13 14from cclib.method.calculationmethod import Method 15from cclib.parser.utils import PeriodicTable 16from cclib.parser.utils import find_package 17 18_found_periodictable = find_package("periodictable") 19if _found_periodictable: 20 import periodictable as pt 21 22_found_scipy = find_package("scipy") 23if _found_scipy: 24 import scipy.constants 25 26 27def _check_periodictable(found_periodictable): 28 if not _found_periodictable: 29 raise ImportError("You must install `periodictable` to use this function") 30 31 32def _check_scipy(found_scipy): 33 if not _found_scipy: 34 raise ImportError("You must install `scipy` to use this function") 35 36 37def get_most_abundant_isotope(element): 38 """Given a `periodictable` element, return the most abundant 39 isotope. 40 """ 41 most_abundant_isotope = element.isotopes[0] 42 abundance = 0 43 for iso in element: 44 if iso.abundance > abundance: 45 most_abundant_isotope = iso 46 abundance = iso.abundance 47 return most_abundant_isotope 48 49 50def get_isotopic_masses(charges): 51 """Return the masses for the given nuclei, respresented by their 52 nuclear charges. 53 """ 54 _check_periodictable(_found_periodictable) 55 masses = [] 56 for charge in charges: 57 el = pt.elements[charge] 58 isotope = get_most_abundant_isotope(el) 59 mass = isotope.mass 60 masses.append(mass) 61 return np.array(masses) 62 63 64class Nuclear(Method): 65 """A container for methods pertaining to atomic nuclei.""" 66 67 def __init__(self, data, progress=None, loglevel=logging.INFO, logname="Log"): 68 69 self.required_attrs = ('natom','atomcoords','atomnos','charge') 70 71 super(Nuclear, self).__init__(data, progress, loglevel, logname) 72 73 def __str__(self): 74 """Return a string representation of the object.""" 75 return "Nuclear" 76 77 def __repr__(self): 78 """Return a representation of the object.""" 79 return "Nuclear" 80 81 def stoichiometry(self): 82 """Return the stoichemistry of the object according to the Hill system""" 83 cclib_pt = PeriodicTable() 84 elements = [cclib_pt.element[ano] for ano in self.data.atomnos] 85 counts = {el: elements.count(el) for el in set(elements)} 86 87 formula = "" 88 elcount = lambda el, c: "%s%i" % (el, c) if c > 1 else el 89 if 'C' in elements: 90 formula += elcount('C', counts['C']) 91 counts.pop('C') 92 if 'H' in elements: 93 formula += elcount('H', counts['H']) 94 counts.pop('H') 95 for el, c in sorted(counts.items()): 96 formula += elcount(el, c) 97 98 if getattr(self.data, 'charge', 0): 99 magnitude = abs(self.data.charge) 100 sign = "+" if self.data.charge > 0 else "-" 101 formula += "(%s%i)" % (sign, magnitude) 102 return formula 103 104 def repulsion_energy(self, atomcoords_index=-1): 105 """Return the nuclear repulsion energy.""" 106 nre = 0.0 107 for i in range(self.data.natom): 108 ri = self.data.atomcoords[atomcoords_index][i] 109 zi = self.data.atomnos[i] 110 for j in range(i+1, self.data.natom): 111 rj = self.data.atomcoords[0][j] 112 zj = self.data.atomnos[j] 113 d = np.linalg.norm(ri-rj) 114 nre += zi*zj/d 115 return nre 116 117 def center_of_mass(self, atomcoords_index=-1): 118 """Return the center of mass.""" 119 charges = self.data.atomnos 120 coords = self.data.atomcoords[atomcoords_index] 121 masses = get_isotopic_masses(charges) 122 123 mwc = coords * masses[:, np.newaxis] 124 numerator = np.sum(mwc, axis=0) 125 denominator = np.sum(masses) 126 127 return numerator / denominator 128 129 def moment_of_inertia_tensor(self, atomcoords_index=-1): 130 """Return the moment of inertia tensor.""" 131 charges = self.data.atomnos 132 coords = self.data.atomcoords[atomcoords_index] 133 masses = get_isotopic_masses(charges) 134 135 moi_tensor = np.empty((3, 3)) 136 137 moi_tensor[0][0] = np.sum(masses * (coords[:, 1]**2 + coords[:, 2]**2)) 138 moi_tensor[1][1] = np.sum(masses * (coords[:, 0]**2 + coords[:, 2]**2)) 139 moi_tensor[2][2] = np.sum(masses * (coords[:, 0]**2 + coords[:, 1]**2)) 140 141 moi_tensor[0][1] = np.sum(masses * coords[:, 0] * coords[:, 1]) 142 moi_tensor[0][2] = np.sum(masses * coords[:, 0] * coords[:, 2]) 143 moi_tensor[1][2] = np.sum(masses * coords[:, 1] * coords[:, 2]) 144 145 moi_tensor[1][0] = moi_tensor[0][1] 146 moi_tensor[2][0] = moi_tensor[0][2] 147 moi_tensor[2][1] = moi_tensor[1][2] 148 149 return moi_tensor 150 151 def principal_moments_of_inertia(self, units='amu_bohr_2'): 152 """Return the principal moments of inertia in 3 kinds of units: 153 1. [amu][bohr]^2 154 2. [amu][angstrom]^2 155 3. [g][cm]^2 156 and the principal axes. 157 """ 158 choices = ('amu_bohr_2', 'amu_angstrom_2', 'g_cm_2') 159 units = units.lower() 160 if units not in choices: 161 raise ValueError("Invalid units, pick one of {}".format(choices)) 162 moi_tensor = self.moment_of_inertia_tensor() 163 principal_moments, principal_axes = np.linalg.eigh(moi_tensor) 164 if units == "amu_bohr_2": 165 _check_scipy(_found_scipy) 166 bohr2ang = scipy.constants.value("atomic unit of length") / scipy.constants.angstrom 167 conv = 1 / bohr2ang ** 2 168 elif units == "amu_angstrom_2": 169 conv = 1 170 elif units == "g_cm_2": 171 _check_scipy(_found_scipy) 172 amu2g = scipy.constants.value("unified atomic mass unit") * scipy.constants.kilo 173 conv = amu2g * (scipy.constants.angstrom / scipy.constants.centi) ** 2 174 return conv * principal_moments, principal_axes 175 176 def rotational_constants(self, units='ghz'): 177 """Compute the rotational constants in 1/cm or GHz.""" 178 choices = ('invcm', 'ghz') 179 units = units.lower() 180 if units not in choices: 181 raise ValueError("Invalid units, pick one of {}".format(choices)) 182 principal_moments = self.principal_moments_of_inertia("amu_angstrom_2")[0] 183 _check_scipy(_found_scipy) 184 bohr2ang = scipy.constants.value('atomic unit of length') / scipy.constants.angstrom 185 xfamu = 1 / scipy.constants.value('electron mass in u') 186 xthz = scipy.constants.value('hartree-hertz relationship') 187 rotghz = xthz * (bohr2ang ** 2) / (2 * xfamu * scipy.constants.giga) 188 if units == 'ghz': 189 conv = rotghz 190 if units == 'invcm': 191 ghz2invcm = scipy.constants.giga * scipy.constants.centi / scipy.constants.c 192 conv = rotghz * ghz2invcm 193 return conv / principal_moments 194 195 196del find_package 197