summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorDavid Luevano Alvarado <55825613+luevano@users.noreply.github.com>2020-02-21 19:00:34 -0700
committerDavid Luevano Alvarado <55825613+luevano@users.noreply.github.com>2020-02-21 19:00:34 -0700
commita2463470b639a743189deece071654ed03791c3c (patch)
treedbeaecde71e3fe1827102f3f30e044ba000ec24b
parent8a108f8757027750f825c6ef9dd55adb4c09def0 (diff)
Moved ljmatrix to representations, rewrote code
-rw-r--r--ml_exp/lj_matrix.py176
-rw-r--r--ml_exp/representations.py102
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