Started work on EAM

snow/lodispp/physics.py

@@ -1,10 +1,13 @@
-import numpy as np
+"""
+Contains functions to ocmpute physical properties for the system, such as pressure, potential energy etc.
+"""
 
-from snow.lodispp.pp_io import read_rgl
+import numpy as np
+from snow.lodispp.pp_io import read_rgl, read_eam
 from tqdm import tqdm
-
-
-def properties(coords: np.ndarray, elements: np.ndarray, pot_file: str, dist_mat: np.ndarray, l_pressure: bool):
+from scipy.spatial.distance import pdist
+from snow.lodispp.utils import nearest_neighbours
+def properties_rgl(coords: np.ndarray, elements: np.ndarray, pot_file: str, dist_mat: np.ndarray, l_pressure: bool):
     """Calculate energy, density and pressure for a given atomic configuration given an interatomic potential parameter file.
 
 
@@ -125,7 +128,81 @@ def properties(coords: np.ndarray, elements: np.ndarray, pot_file: str, dist_mat
         if l_pressure:
             press_atoms[i] = (presr_i + presb_i) / at_vol[0]
 
-    return density, pot_energy, press_atoms
+    return density, -pot_energy, press_atoms
 
 
 
+def pair_energy_eam(pot_file: str, coords: np.ndarray) -> float:
+    """Computes the energy associated with a pair of atoms
+
+    Parameters
+    ----------
+    pot_file : str
+        Path to the EAM potential file
+    coords : np.ndarray
+        Coordinates of two atoms of a pair
+
+    Returns
+    -------
+    float
+        Potential energy of a pair
+    """
+
+    r_ij = pdist(coords)[0]
+
+    potential = read_eam(pot_file)
+
+    idx_closer_r = np.searchsorted(potential["r"], r_ij)
+
+    rho_r = potential["rho_r"][idx_closer_r]
+
+    Z_r = potential["Z_r"][idx_closer_r]
+    phi_r = 27.2 * 0.529 * Z_r * Z_r / r_ij
+
+    idx_closer_rho = np.searchsorted(potential["rho_r"], rho_r)
+
+
+    F_rho = potential["F_rho"][idx_closer_rho]
+
+    return F_rho + 0.5 * phi_r
+
+def energy_eam(coords: np.ndarray, pot_file: str) -> np.ndarray:
+    """_summary_
+
+    Parameters
+    ----------
+    coords : np.ndarray
+        _description_
+    pot_file : str
+        _description_
+
+    Returns
+    -------
+    np.ndarray
+        _description_
+    """
+
+    N_atoms = np.shape(coords)[0]
+
+    potential = read_eam(pot_file)
+    cut_off = potential["cut_off"]
+
+    neigh_list = nearest_neighbours(1, coords = coords, cut_off = cut_off)
+
+    pot_en = np.zeros(N_atoms)
+    print("Computing energies for each atom")
+    for i in tqdm(range (N_atoms)):
+        en_i = 0
+        for neigh in neigh_list[i]:
+            if neigh != i:
+                coords_ij = np.array([coords[i], coords[neigh]])
+                en_ij = pair_energy_eam(pot_file=pot_file, coords=coords_ij)
+                en_i += en_ij
+        pot_en[i] = en_i
+    return pot_en
+
+
+
+
+
+

snow/lodispp/pp_io.py

@@ -1,9 +1,85 @@
+""" Contains input output functions related to reading structures, potential files etc. """
 from typing import Tuple
 import numpy as np
 import os
 import inspect
 
-def read_rgl(filepot: str):
+
+def read_eam(filepot: str) -> dict:
+    """Reads an EAM (Embedded Atom Model) potential and returns a dictionary with all factors and parameters
+
+    Parameters
+    ----------
+    filepot : str
+        _description_
+
+    Returns
+    -------
+    dict
+        _description_
+    """
+
+    with open(filepot, "r") as pot_file:
+        _ = pot_file.readline()
+        atomic_number, mass, lat_param, lat_type = pot_file.readline().split()
+        atomic_number = int(atomic_number)
+        mass = float(mass)
+        lat_param = float(lat_param)
+
+        n_rho, d_rho, n_r, d_r, r_cut = map(float, pot_file.readline().split())
+        F_rho = []
+        Z_r = []
+        rho_r = []
+
+        r = np.arange(0, (n_r - 1) * d_r, d_r)
+        rho = np.arange(0, (n_rho - 1) * d_rho, d_rho)
+
+        for _ in range(int(n_rho) // 5):
+            line = pot_file.readline().split()
+            for i in range(5):
+
+                F_rho.append(float(line[i]))
+
+        for _ in range(int(n_r) // 5):
+            line = pot_file.readline().split()
+            for i in range(5):
+
+                Z_r.append(float(line[i]))
+
+        for _ in range(int(n_r) // 5):
+            line = pot_file.readline().split()
+            for i in range(5):
+
+                rho_r.append(float(line[i]))
+
+
+    return {
+        "F_rho": F_rho,
+        "Z_r": Z_r,
+        "rho_r": rho_r,
+        "r": r,
+        "rho": rho,
+        "cut_off": r_cut
+    }
+
+
+
+
+
+def read_rgl(filepot: str) -> dict:
+    """Reads the parameters from a RGL potential file, computes the necessary factors and outputs a dictionary with the necessary 
+    factors and parameters
+
+    Parameters
+    ----------
+    filepot : str
+        Path to the potential parameter file
+
+    Returns
+    -------
+    dict
+        Dictionary with all the parameters and factors for the RGL potential
+    """
     with open(filepot, 'r') as file:
         lines = file.readlines()
 
@@ -160,7 +236,7 @@ def read_xyz(file_path: str) -> Tuple[np.ndarray, np.ndarray]:
         filepath = os.path.join(script_dir, file_path)
 
         # Open the file
-        with open(filepath, "r") as xyz_file:
+        with open(file_path, "r") as xyz_file:
             # Number of atoms
             n_atoms = int(xyz_file.readline().strip())
 

snow/lodispp/utils.py

@@ -103,8 +103,9 @@ def pddf_calculator(index_frame, coords, bin_precision = None, bin_count = None)
     for i in range(n_bins_int):
         for j in  range(n_atoms):
             for k in range(n_atoms):
-                if (dist_mat[j,k] < bin_precision * i and dist_mat[j,k] >= (bin_precision * (i - 1))):
-                    dist_count[i] += 1
+                if k != j:
+                    if (dist_mat[j,k] < bin_precision * i and dist_mat[j,k] >= (bin_precision * (i - 1))):
+                        dist_count[i] += 1
         dist[i] = bin_precision * i
 
     return dist, dist_count