Coverage for cclib/parser/utils.py : 73%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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.
8"""Utilities often used by cclib parsers and scripts"""
10import importlib
11import sys
12from itertools import accumulate
14import numpy
15import periodictable
18def find_package(package):
19 """Check if a package exists without importing it.
21 Derived from https://stackoverflow.com/a/14050282
22 """
23 module_spec = importlib.util.find_spec(package)
24 return module_spec is not None and module_spec.loader is not None
27_found_scipy = find_package("scipy")
28if _found_scipy:
29 import scipy.spatial
32def symmetrize(m, use_triangle='lower'):
33 """Symmetrize a square NumPy array by reflecting one triangular
34 section across the diagonal to the other.
35 """
37 if use_triangle not in ('lower', 'upper'):
38 raise ValueError
39 if not len(m.shape) == 2:
40 raise ValueError
41 if not (m.shape[0] == m.shape[1]):
42 raise ValueError
44 dim = m.shape[0]
46 lower_indices = numpy.tril_indices(dim, k=-1)
47 upper_indices = numpy.triu_indices(dim, k=1)
49 ms = m.copy()
51 if use_triangle == 'lower':
52 ms[upper_indices] = ms[lower_indices]
53 if use_triangle == 'upper':
54 ms[lower_indices] = ms[upper_indices]
56 return ms
59_BUILTIN_FLOAT = float
62def float(number):
63 """Convert a string to a float.
65 This method should perform certain checks that are specific to cclib,
66 including avoiding the problem with Ds instead of Es in scientific notation.
67 Another point is converting string signifying numerical problems (*****)
68 to something we can manage (Numpy's NaN).
69 """
71 if list(set(number)) == ['*']:
72 return numpy.nan
74 return _BUILTIN_FLOAT(number.replace("D", "E"))
77def convertor(value, fromunits, tounits):
78 """Convert from one set of units to another.
80 Sources:
81 NIST 2010 CODATA (http://physics.nist.gov/cuu/Constants/index.html)
82 Documentation of GAMESS-US or other programs as noted
83 """
85 _convertor = {
87 "time_au_to_fs": lambda x: x * 0.02418884,
88 "fs_to_time_au": lambda x: x / 0.02418884,
90 "Angstrom_to_bohr": lambda x: x * 1.8897261245,
91 "bohr_to_Angstrom": lambda x: x * 0.5291772109,
93 "wavenumber_to_eV": lambda x: x / 8065.54429,
94 "wavenumber_to_hartree": lambda x: x / 219474.6313708,
95 "wavenumber_to_kcal/mol": lambda x: x / 349.7550112,
96 "wavenumber_to_kJ/mol": lambda x: x / 83.5934722814,
97 "wavenumber_to_nm": lambda x: 1e7 / x,
98 "wavenumber_to_Hz": lambda x: x * 29.9792458,
100 "eV_to_wavenumber": lambda x: x * 8065.54429,
101 "eV_to_hartree": lambda x: x / 27.21138505,
102 "eV_to_kcal/mol": lambda x: x * 23.060548867,
103 "eV_to_kJ/mol": lambda x: x * 96.4853364596,
105 "hartree_to_wavenumber": lambda x: x * 219474.6313708,
106 "hartree_to_eV": lambda x: x * 27.21138505,
107 "hartree_to_kcal/mol": lambda x: x * 627.50947414,
108 "hartree_to_kJ/mol": lambda x: x * 2625.4996398,
110 "kcal/mol_to_wavenumber": lambda x: x * 349.7550112,
111 "kcal/mol_to_eV": lambda x: x / 23.060548867,
112 "kcal/mol_to_hartree": lambda x: x / 627.50947414,
113 "kcal/mol_to_kJ/mol": lambda x: x * 4.184,
115 "kJ/mol_to_wavenumber": lambda x: x * 83.5934722814,
116 "kJ/mol_to_eV": lambda x: x / 96.4853364596,
117 "kJ/mol_to_hartree": lambda x: x / 2625.49963978,
118 "kJ/mol_to_kcal/mol": lambda x: x / 4.184,
119 "nm_to_wavenumber": lambda x: 1e7 / x,
121 # Taken from GAMESS docs, "Further information",
122 # "Molecular Properties and Conversion Factors"
123 "Debye^2/amu-Angstrom^2_to_km/mol": lambda x: x * 42.255,
125 # Conversion for charges and multipole moments.
126 "e_to_coulomb": lambda x: x * 1.602176565 * 1e-19,
127 "e_to_statcoulomb": lambda x: x * 4.80320425 * 1e-10,
128 "coulomb_to_e": lambda x: x * 0.6241509343 * 1e19,
129 "statcoulomb_to_e": lambda x: x * 0.2081943527 * 1e10,
130 "ebohr_to_Debye": lambda x: x * 2.5417462300,
131 "ebohr2_to_Buckingham": lambda x: x * 1.3450341749,
132 "ebohr2_to_Debye.ang": lambda x: x * 1.3450341749,
133 "ebohr3_to_Debye.ang2": lambda x: x * 0.7117614302,
134 "ebohr4_to_Debye.ang3": lambda x: x * 0.3766479268,
135 "ebohr5_to_Debye.ang4": lambda x: x * 0.1993134985,
137 "hartree/bohr2_to_mDyne/angstrom": lambda x: x * 8.23872350 / 0.5291772109
138 }
140 return _convertor["%s_to_%s" % (fromunits, tounits)](value)
142def _get_rmat_from_vecs(a, b):
143 """Get rotation matrix from two 3D vectors, a and b
144 Args:
145 a (np.ndaray): 3d vector with shape (3,0)
146 b (np.ndaray): 3d vector with shape (3,0)
147 Returns:
148 np.ndarray
149 """
150 a_ = (a / numpy.linalg.norm(a, 2))
151 b_ = (b / numpy.linalg.norm(b, 2))
152 v = numpy.cross(a_, b_)
153 s = numpy.linalg.norm(v, 2)
154 c = numpy.dot(a_, b_)
155 # skew-symmetric cross product of v
156 vx = numpy.array([[0, -v[2], v[1]],
157 [v[2], 0, -v[0]],
158 [-v[1], v[0], 0]])
159 rmat = numpy.identity(3) + vx + numpy.matmul(vx, vx) * ((1-c)/s**2)
160 return rmat
162def get_rotation(a, b):
163 """Get rotation part for transforming a to b, where a and b are same positions with different orientations
164 If one atom positions, i.e (1,3) shape array, are given, it returns identify transformation
166 Args:
167 a (np.ndarray): positions with shape(N,3)
168 b (np.ndarray): positions with shape(N,3)
169 Returns:
170 A scipy.spatial.transform.Rotation object
171 """
172 if not _found_scipy:
173 raise ImportError("You must install `scipy` to use this function")
175 assert a.shape == b.shape
176 if a.shape[0] == 1:
177 return scipy.spatial.transform.Rotation.from_euler('xyz', [0,0,0])
178 # remove translation part
179 a_ = a - a[0]
180 b_ = b - b[0]
181 if hasattr(scipy.spatial.transform.Rotation, "align_vectors"):
182 r, _ = scipy.spatial.transform.Rotation.align_vectors(b_, a_)
183 else:
184 if numpy.linalg.matrix_rank(a_) == 1:
185 # in the case of linear molecule, e.g. O2, C2H2
186 idx = numpy.argmax(numpy.linalg.norm(a_, ord=2, axis=1))
187 rmat = _get_rmat_from_vecs(a_[idx], b_[idx])
188 r = scipy.spatial.transform.Rotation.from_dcm(rmat)
189 else:
190 # scipy.spatial.transform.Rotation.match_vectors has bug
191 # Kabsch Algorithm
192 cov = numpy.dot(b_.T, a_)
193 V, S, W = numpy.linalg.svd(cov)
194 if ((numpy.linalg.det(V) * numpy.linalg.det(W))< 0.0):
195 S[-1] = -S[-1]
196 V[:,-1] = -V[:,-1]
197 rmat = numpy.dot(V, W)
198 r = scipy.spatial.transform.Rotation.from_dcm(rmat)
199 return r
201class PeriodicTable:
202 """Allows conversion between element name and atomic no."""
204 def __init__(self):
205 self.element = [None]
206 self.number = {}
208 for e in periodictable.elements:
209 if e.symbol != 'n':
210 self.element.append(e.symbol)
211 self.number[e.symbol] = e.number
214class WidthSplitter:
215 """Split a line based not on a character, but a given number of field
216 widths.
217 """
219 def __init__(self, widths):
220 self.start_indices = [0] + list(accumulate(widths))[:-1]
221 self.end_indices = list(accumulate(widths))
223 def split(self, line, truncate=True):
224 """Split the given line using the field widths passed in on class
225 initialization.
226 """
227 elements = [line[start:end].strip()
228 for (start, end) in zip(self.start_indices, self.end_indices)]
229 # Handle lines that contain fewer fields than specified in the
230 # widths; they are added as empty strings, so remove them.
231 if truncate:
232 while len(elements) and elements[-1] == '':
233 elements.pop()
234 return elements