Moved ljmatrix to representations, rewrote code

2 files changed, 102 insertions, 176 deletions

diff --git a/ml_exp/lj_matrix.py b/ml_exp/lj_matrix.py
deleted file mode 100644
index 5e3e89678..000000000
--- a/ml_exp/lj_matrix.py
+++ /dev/null
@@ -1,176 +0,0 @@
-"""MIT License
-
-Copyright (c) 2019 David Luevano Alvarado
-
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in all
-copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
-"""
-import math
-import numpy as np
-from numpy.linalg import eig
-
-
-def lj_matrix(mol_data,
-              nc_data,
-              diag_value=None,
-              sigma=1.0,
-              epsilon=1.0,
-              max_len=25,
-              as_eig=True,
-              bohr_radius_units=False):
-    """
-    Creates the Lennard-Jones Matrix from the molecule data given.
-    mol_data: molecule data, matrix of atom coordinates.
-    nc_data: nuclear charge data, array of atom data.
-    diag_value: if special diagonal value is to be used.
-    sigma: sigma value.
-    epsilon: epsilon value.
-    max_len: maximum amount of atoms in molecule.
-    as_eig: if data should be returned as matrix or array of eigenvalues.
-    bohr_radius_units: if units should be in bohr's radius units.
-    """
-    if bohr_radius_units:
-        conversion_rate = 0.52917721067
-    else:
-        conversion_rate = 1
-
-    mol_n = len(mol_data)
-    mol_nr = range(mol_n)
-
-    if not mol_n == len(nc_data):
-        print(''.join(['Error. Molecule matrix dimension is different ',
-                       'than the nuclear charge array dimension.']))
-    else:
-        if max_len < mol_n:
-            print(''.join(['Error. Molecule matrix dimension (mol_n) is ',
-                           'greater than max_len. Using mol_n.']))
-            max_len = None
-
-        if max_len:
-            lj = np.zeros((max_len, max_len))
-            ml_r = range(max_len)
-
-            # Actual calculation of the coulomb matrix.
-            for i in ml_r:
-                if i < mol_n:
-                    x_i = mol_data[i, 0]
-                    y_i = mol_data[i, 1]
-                    z_i = mol_data[i, 2]
-                    Z_i = nc_data[i]
-                else:
-                    break
-
-                for j in ml_r:
-                    if j < mol_n:
-                        x_j = mol_data[j, 0]
-                        y_j = mol_data[j, 1]
-                        z_j = mol_data[j, 2]
-
-                        x = (x_i-x_j)**2
-                        y = (y_i-y_j)**2
-                        z = (z_i-z_j)**2
-
-                        if i == j:
-                            if diag_value is None:
-                                lj[i, j] = (0.5*Z_i**2.4)
-                            else:
-                                lj[i, j] = diag_value
-                        else:
-                            # Calculations are done after i==j is checked
-                            # so no division by zero is done.
-
-                            # A little play with r exponents
-                            # so no square root is calculated.
-                            # Conversion factor is included in r^2.
-
-                            # 1/r^2
-                            r_2 = sigma**2/(conversion_rate**2*(x + y + z))
-
-                            r_6 = math.pow(r_2, 3)
-                            r_12 = math.pow(r_6, 2)
-                            lj[i, j] = (4*epsilon*(r_12 - r_6))
-                    else:
-                        break
-
-            # Now the value will be returned.
-            if as_eig:
-                lj_sorted = np.sort(eig(lj)[0])[::-1]
-                # Thanks to SO for the following lines of code.
-                # https://stackoverflow.com/a/43011036
-
-                # Keep zeros at the end.
-                mask = lj_sorted != 0.
-                f_mask = mask.sum(0, keepdims=1) >\
-                    np.arange(lj_sorted.shape[0]-1, -1, -1)
-
-                f_mask = f_mask[::-1]
-                lj_sorted[f_mask] = lj_sorted[mask]
-                lj_sorted[~f_mask] = 0.
-
-                return lj_sorted
-
-            else:
-                return lj
-
-        else:
-            lj_temp = []
-            # Actual calculation of the coulomb matrix.
-            for i in mol_nr:
-                x_i = mol_data[i, 0]
-                y_i = mol_data[i, 1]
-                z_i = mol_data[i, 2]
-                Z_i = nc_data[i]
-
-                lj_row = []
-                for j in mol_nr:
-                    x_j = mol_data[j, 0]
-                    y_j = mol_data[j, 1]
-                    z_j = mol_data[j, 2]
-
-                    x = (x_i-x_j)**2
-                    y = (y_i-y_j)**2
-                    z = (z_i-z_j)**2
-
-                    if i == j:
-                        if not diag_value:
-                            lj_row.append(0.5*Z_i**2.4)
-                        else:
-                            lj_row.append(diag_value)
-                    else:
-                        # Calculations are done after i==j is checked
-                        # so no division by zero is done.
-
-                        # A little play with r exponents
-                        # so no square root is calculated.
-                        # Conversion factor is included in r^2.
-
-                        # 1/r^2
-                        r_2 = sigma**2/(conversion_rate**2*(x + y + z))
-
-                        r_6 = math.pow(r_2, 3)
-                        r_12 = math.pow(r_6, 2)
-                        lj_row.append(4*epsilon*(r_12 - r_6))
-
-                lj_temp.append(np.array(lj_row))
-
-            lj = np.array(lj_temp)
-            # Now the value will be returned.
-            if as_eig:
-                return np.sort(eig(lj)[0])[::-1]
-            else:
-                return lj
diff --git a/ml_exp/representations.py b/ml_exp/representations.py
index c503c310d..b4edef8a5 100644
--- a/ml_exp/representations.py
+++ b/ml_exp/representations.py
@@ -102,3 +102,105 @@ def coulomb_matrix(coords,
         return cm_sorted
     else:
         return cm
+
+
+def lennard_jones_matrix(coords,
+                         nc,
+                         diag_value=None,
+                         sigma=1.0,
+                         epsilon=1.0,
+                         size=23,
+                         as_eig=True,
+                         bhor_ru=False):
+    """
+    Creates the Lennard-Jones Matrix from the molecule data given.
+    coords: compound coordinates.
+    nc: nuclear charge data.
+    diag_value: if special diagonal value is to be used.
+    sigma: sigma value.
+    epsilon: epsilon value.
+    size: compound size.
+    as_eig: if the representation should be as the eigenvalues.
+    bhor_ru: if radius units should be in bohr's radius units.
+    """
+    if bhor_ru:
+        cr = 0.52917721067
+    else:
+        cr = 1
+
+    n = coords.shape[0]
+
+    if not n == nc.shape[0]:
+        raise ValueError('Compound size is different than the nuclear charge\
+                         size. Arrays are not of the right shape.')
+
+    if size < n:
+        print('Error. Compound size (n) is greater han (size). Using (n)\
+              instead of (size).')
+        size = n
+
+    nr = range(size)
+    lj = np.zeros((size, size), dtype=float)
+
+    # Actual calculation of the lennard-jones matrix.
+    for i in nr:
+        if i < n:
+            x_i = coords[i, 0]
+            y_i = coords[i, 1]
+            z_i = coords[i, 2]
+            Z_i = nc[i]
+        else:
+            break
+
+        for j in nr:
+            if j < n:
+                x_j = coords[j, 0]
+                y_j = coords[j, 1]
+                z_j = coords[j, 2]
+                # Never really used because when needed it is equal to Z_i.
+                # Z_j = nc[j]
+
+                x = (x_i-x_j)**2
+                y = (y_i-y_j)**2
+                z = (z_i-z_j)**2
+
+                if i == j:
+                    if diag_value is None:
+                        lj[i, j] = (0.5*Z_i**2.4)
+                    else:
+                        lj[i, j] = diag_value
+                else:
+                    # Calculations are done after i==j is checked
+                    # so no division by zero is done.
+
+                    # A little play with r exponents
+                    # so no square root is calculated.
+                    # Conversion factor is included in r^2.
+
+                    # 1/r^2
+                    r_2 = sigma**2/(cr**2*(x + y + z))
+
+                    r_6 = math.pow(r_2, 3)
+                    r_12 = math.pow(r_6, 2)
+                    lj[i, j] = (4*epsilon*(r_12 - r_6))
+            else:
+                break
+
+    # Now the value will be returned.
+    if as_eig:
+        lj_sorted = np.sort(np.linalg.eig(lj)[0])[::-1]
+        # Thanks to SO for the following lines of code.
+        # https://stackoverflow.com/a/43011036
+
+        # Keep zeros at the end.
+        mask = lj_sorted != 0.
+        f_mask = mask.sum(0, keepdims=1) >\
+            np.arange(lj_sorted.shape[0]-1, -1, -1)
+
+        f_mask = f_mask[::-1]
+        lj_sorted[f_mask] = lj_sorted[mask]
+        lj_sorted[~f_mask] = 0.
+
+        return lj_sorted
+    else:
+        return lj