Import numpy as np in Arkane statmech

arkane/statmech.py

@@ -36,7 +36,7 @@
 
 import os.path
 import math
-import numpy
+import numpy as np
 import logging
 
 from rdkit.Chem import GetPeriodicTable
@@ -120,8 +120,8 @@ def load(self):
                     angles.append(float(tokens[0]) / angleFactor)
                     energies.append(float(tokens[1]) / energyFactor)
 
-        angles = numpy.array(angles)
-        energies = numpy.array(energies)
+        angles = np.array(angles)
+        energies = np.array(energies)
         energies -= energies[0]
 
         return angles, energies
@@ -489,8 +489,8 @@ def load(self, pdep=False):
                     fourierRotor = HinderedRotor(inertia=(inertia, "amu*angstrom^2"), symmetry=symmetry)
                     fourierRotor.fitFourierPotentialToData(angle, v_list)
 
-                    Vlist_cosine = numpy.zeros_like(angle)
-                    Vlist_fourier = numpy.zeros_like(angle)
+                    Vlist_cosine = np.zeros_like(angle)
+                    Vlist_fourier = np.zeros_like(angle)
                     for i in range(angle.shape[0]):
                         Vlist_cosine[i] = cosineRotor.getPotential(angle[i])
                         Vlist_fourier[i] = fourierRotor.getPotential(angle[i])
@@ -504,9 +504,9 @@ def load(self, pdep=False):
                         rotorCount += 1
                         conformer.modes.append(rotor)
                     elif fit == 'best':
-                        rms_cosine = numpy.sqrt(numpy.sum((Vlist_cosine - v_list) * (Vlist_cosine - v_list)) /
+                        rms_cosine = np.sqrt(np.sum((Vlist_cosine - v_list) * (Vlist_cosine - v_list)) /
                                                 (len(v_list) - 1)) / 4184.
-                        rms_fourier = numpy.sqrt(numpy.sum((Vlist_fourier - v_list) * (Vlist_fourier - v_list)) /
+                        rms_fourier = np.sqrt(np.sum((Vlist_fourier - v_list) * (Vlist_fourier - v_list)) /
                                                  (len(v_list) - 1)) / 4184.
 
                         # Keep the rotor with the most accurate potential
@@ -525,7 +525,7 @@ def load(self, pdep=False):
                         rotorCount += 1
 
             logging.debug('    Determining frequencies from reduced force constant matrix...')
-            frequencies = numpy.array(projectRotors(conformer, F, rotors, linear, is_ts, label=self.species.label))
+            frequencies = np.array(projectRotors(conformer, F, rotors, linear, is_ts, label=self.species.label))
 
         elif len(conformer.modes) > 2:
             if len(rotors) > 0:
@@ -534,15 +534,15 @@ def load(self, pdep=False):
                              ' Gaussian to generate the force constant matrix, if running Molpro include keyword print,'
                              ' hessian')
             frequencies = conformer.modes[2].frequencies.value_si
-            rotors = numpy.array([])
+            rotors = np.array([])
         else:
             if len(rotors) > 0:
                 logging.warn('Force Constant Matrix Missing Ignoring rotors, if running Gaussian if not already'
                              ' present you need to add the keyword iop(7/33=1) in your Gaussian frequency job for'
                              ' Gaussian to generate the force constant matrix, if running Molpro include keyword print,'
                              ' hessian')
-            frequencies = numpy.array([])
-            rotors = numpy.array([])
+            frequencies = np.array([])
+            rotors = np.array([])
 
         for mode in conformer.modes:
             if isinstance(mode, HarmonicOscillator):
@@ -596,9 +596,9 @@ def plotHinderedRotor(self, angle, v_list, cosineRotor, fourierRotor, rotor, rot
         except ImportError:
             return
 
-        phi = numpy.arange(0, 6.3, 0.02, numpy.float64)
-        Vlist_cosine = numpy.zeros_like(phi)
-        Vlist_fourier = numpy.zeros_like(phi)
+        phi = np.arange(0, 6.3, 0.02, np.float64)
+        Vlist_cosine = np.zeros_like(phi)
+        Vlist_fourier = np.zeros_like(phi)
         for i in range(phi.shape[0]):
             Vlist_cosine[i] = cosineRotor.getPotential(phi[i])
             Vlist_fourier[i] = fourierRotor.getPotential(phi[i])
@@ -984,12 +984,12 @@ def is_linear(coordinates):
         return True
 
     # A tensor containing all distance vectors in the molecule
-    d = -numpy.array([c[:, numpy.newaxis] - c[numpy.newaxis, :] for c in coordinates.T])
+    d = -np.array([c[:, np.newaxis] - c[np.newaxis, :] for c in coordinates.T])
     for i in range(2, len(coordinates)):
         if i > 1:
-            u1 = d[:, 0, 1] / numpy.linalg.norm(d[:, 0, 1])  # unit vector between atoms 0 and 1
-            u2 = d[:, 1, i] / numpy.linalg.norm(d[:, 1, i])  # unit vector between atoms 1 and i
-            a = math.degrees(numpy.arccos(numpy.clip(numpy.dot(u1, u2), -1.0, 1.0)))  # angle between atoms 0, 1, i
+            u1 = d[:, 0, 1] / np.linalg.norm(d[:, 0, 1])  # unit vector between atoms 0 and 1
+            u2 = d[:, 1, i] / np.linalg.norm(d[:, 1, i])  # unit vector between atoms 1 and i
+            a = math.degrees(np.arccos(np.clip(np.dot(u1, u2), -1.0, 1.0)))  # angle between atoms 0, 1, i
             if abs(180 - a) > epsilon and abs(a) > epsilon:
                 return False
     return True
@@ -1036,22 +1036,22 @@ def projectRotors(conformer, F, rotors, linear, is_ts, label):
         coordinates[i, 1] -= ym
         coordinates[i, 2] -= zm
     # Make vector with the root of the mass in amu for each atom
-    amass = numpy.sqrt(mass / constants.amu)
+    amass = np.sqrt(mass / constants.amu)
 
     # Rotation matrix
     I = conformer.getMomentOfInertiaTensor()
-    PMoI, Ixyz = numpy.linalg.eigh(I)
+    PMoI, Ixyz = np.linalg.eigh(I)
 
     external = 6
     if linear:
         external = 5
 
-    D = numpy.zeros((Natoms * 3, external), numpy.float64)
+    D = np.zeros((Natoms * 3, external), np.float64)
 
-    P = numpy.zeros((Natoms, 3), numpy.float64)
+    P = np.zeros((Natoms, 3), np.float64)
 
     # Transform the coordinates to the principal axes
-    P = numpy.dot(coordinates, Ixyz)
+    P = np.dot(coordinates, Ixyz)
 
     for i in range(Natoms):
         # Projection vectors for translation
@@ -1075,9 +1075,9 @@ def projectRotors(conformer, F, rotors, linear, is_ts, label):
     # Make sure projection matrix is orthonormal
     import scipy.linalg
 
-    I = numpy.identity(Natoms * 3, numpy.float64)
+    I = np.identity(Natoms * 3, np.float64)
 
-    P = numpy.zeros((Natoms * 3, 3 * Natoms + external), numpy.float64)
+    P = np.zeros((Natoms * 3, 3 * Natoms + external), np.float64)
 
     P[:, 0:external] = D[:, 0:external]
     P[:, external:external + 3 * Natoms] = I[:, 0:3 * Natoms]
@@ -1088,7 +1088,7 @@ def projectRotors(conformer, F, rotors, linear, is_ts, label):
             norm += P[j, i] * P[j, i]
         for j in range(3 * Natoms):
             if norm > 1E-15:
-                P[j, i] /= numpy.sqrt(norm)
+                P[j, i] /= np.sqrt(norm)
             else:
                 P[j, i] = 0.0
         for j in range(i + 1, 3 * Natoms + external):
@@ -1110,7 +1110,7 @@ def projectRotors(conformer, F, rotors, linear, is_ts, label):
             i += 1
 
     # T is the transformation vector from cartesian to internal coordinates
-    T = numpy.zeros((Natoms * 3, 3 * Natoms - external), numpy.float64)
+    T = np.zeros((Natoms * 3, 3 * Natoms - external), np.float64)
 
     T[:, 0:3 * Natoms - external] = P[:, external:3 * Natoms]
 
@@ -1124,19 +1124,19 @@ def projectRotors(conformer, F, rotors, linear, is_ts, label):
                 for v in range(3):
                     Fm[3 * i + u, 3 * j + v] /= math.sqrt(mass[i] * mass[j])
 
-    Fint = numpy.dot(T.T, numpy.dot(Fm, T))
+    Fint = np.dot(T.T, np.dot(Fm, T))
 
     # Get eigenvalues of internal force constant matrix, V = 3N-6 * 3N-6
-    eig, V = numpy.linalg.eigh(Fint)
+    eig, V = np.linalg.eigh(Fint)
 
     logging.debug('Frequencies from internal Hessian')
     for i in range(3 * Natoms - external):
-        with numpy.warnings.catch_warnings():
-            numpy.warnings.filterwarnings('ignore', r'invalid value encountered in sqrt')
-            logging.debug(numpy.sqrt(eig[i]) / (2 * math.pi * constants.c * 100))
+        with np.warnings.catch_warnings():
+            np.warnings.filterwarnings('ignore', r'invalid value encountered in sqrt')
+            logging.debug(np.sqrt(eig[i]) / (2 * math.pi * constants.c * 100))
 
     # Now we can start thinking about projecting out the internal rotations
-    Dint = numpy.zeros((3 * Natoms, Nrotors), numpy.float64)
+    Dint = np.zeros((3 * Natoms, Nrotors), np.float64)
 
     counter = 0
     for i, rotor in enumerate(rotors):
@@ -1159,24 +1159,24 @@ def projectRotors(conformer, F, rotors, linear, is_ts, label):
             atom = j + 1
             if atom in top:
                 e31 = coordinates[atom - 1, :] - coordinates[pivot1 - 1, :]
-                Dint[3 * (atom - 1):3 * (atom - 1) + 3, counter] = numpy.cross(e31, e12) * amass[atom - 1]
+                Dint[3 * (atom - 1):3 * (atom - 1) + 3, counter] = np.cross(e31, e12) * amass[atom - 1]
             else:
                 e31 = coordinates[atom - 1, :] - coordinates[pivot2 - 1, :]
-                Dint[3 * (atom - 1):3 * (atom - 1) + 3, counter] = numpy.cross(e31, -e12) * amass[atom - 1]
+                Dint[3 * (atom - 1):3 * (atom - 1) + 3, counter] = np.cross(e31, -e12) * amass[atom - 1]
         counter += 1
 
     # Normal modes in mass weighted cartesian coordinates
-    Vmw = numpy.dot(T, V)
-    eigM = numpy.zeros((3 * Natoms - external, 3 * Natoms - external), numpy.float64)
+    Vmw = np.dot(T, V)
+    eigM = np.zeros((3 * Natoms - external, 3 * Natoms - external), np.float64)
 
     for i in range(3 * Natoms - external):
         eigM[i, i] = eig[i]
 
-    Fm = numpy.dot(Vmw, numpy.dot(eigM, Vmw.T))
+    Fm = np.dot(Vmw, np.dot(eigM, Vmw.T))
 
     # Internal rotations are not normal modes => project them on the normal modes and orthogonalize
     # Dintproj =  (3N-6) x (3N) x (3N) x (Nrotors)
-    Dintproj = numpy.dot(Vmw.T, Dint)
+    Dintproj = np.dot(Vmw.T, Dint)
 
     # Reconstruct Dint
     for i in range(Nrotors):
@@ -1191,33 +1191,33 @@ def projectRotors(conformer, F, rotors, linear, is_ts, label):
         for j in range(3 * Natoms):
             norm += Dint[j, i] * Dint[j, i]
         for j in range(3 * Natoms):
-            Dint[j, i] /= numpy.sqrt(norm)
+            Dint[j, i] /= np.sqrt(norm)
         for j in range(i + 1, Nrotors):
             proj = 0.0
             for k in range(3 * Natoms):
                 proj += Dint[k, i] * Dint[k, j]
             for k in range(3 * Natoms):
                 Dint[k, j] -= proj * Dint[k, i]
 
-    Dintproj = numpy.dot(Vmw.T, Dint)
-    Proj = numpy.dot(Dint, Dint.T)
-    I = numpy.identity(Natoms * 3, numpy.float64)
+    Dintproj = np.dot(Vmw.T, Dint)
+    Proj = np.dot(Dint, Dint.T)
+    I = np.identity(Natoms * 3, np.float64)
     Proj = I - Proj
-    Fm = numpy.dot(Proj, numpy.dot(Fm, Proj))
+    Fm = np.dot(Proj, np.dot(Fm, Proj))
     # Get eigenvalues of mass-weighted force constant matrix
-    eig, V = numpy.linalg.eigh(Fm)
+    eig, V = np.linalg.eigh(Fm)
     eig.sort()
 
     # Convert eigenvalues to vibrational frequencies in cm^-1
     # Only keep the modes that don't correspond to translation, rotation, or internal rotation
 
     logging.debug('Frequencies from projected Hessian')
     for i in range(3 * Natoms):
-        with numpy.warnings.catch_warnings():
-            numpy.warnings.filterwarnings('ignore', r'invalid value encountered in sqrt')
-            logging.debug(numpy.sqrt(eig[i]) / (2 * math.pi * constants.c * 100))
+        with np.warnings.catch_warnings():
+            np.warnings.filterwarnings('ignore', r'invalid value encountered in sqrt')
+            logging.debug(np.sqrt(eig[i]) / (2 * math.pi * constants.c * 100))
 
-    return numpy.sqrt(eig[-Nvib:]) / (2 * math.pi * constants.c * 100)
+    return np.sqrt(eig[-Nvib:]) / (2 * math.pi * constants.c * 100)
 
 
 def assign_frequency_scale_factor(model_chemistry):