From a2463470b639a743189deece071654ed03791c3c Mon Sep 17 00:00:00 2001 From: David Luevano Alvarado <55825613+luevano@users.noreply.github.com> Date: Fri, 21 Feb 2020 19:00:34 -0700 Subject: Moved ljmatrix to representations, rewrote code --- ml_exp/lj_matrix.py | 176 ---------------------------------------------- ml_exp/representations.py | 102 +++++++++++++++++++++++++++ 2 files changed, 102 insertions(+), 176 deletions(-) delete mode 100644 ml_exp/lj_matrix.py 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 -- cgit v1.2.3-70-g09d2