diff --git a/snow/lodispp/__pycache__/physics.cpython-312.pyc b/snow/lodispp/__pycache__/physics.cpython-312.pyc new file mode 100644 index 0000000..34cc083 Binary files /dev/null and b/snow/lodispp/__pycache__/physics.cpython-312.pyc differ diff --git a/snow/lodispp/__pycache__/pp_io.cpython-312.pyc b/snow/lodispp/__pycache__/pp_io.cpython-312.pyc index 8d04bef..9fcde9d 100644 Binary files a/snow/lodispp/__pycache__/pp_io.cpython-312.pyc and b/snow/lodispp/__pycache__/pp_io.cpython-312.pyc differ diff --git a/snow/lodispp/__pycache__/utils.cpython-312.pyc b/snow/lodispp/__pycache__/utils.cpython-312.pyc index b2fdfda..c8a2fdd 100644 Binary files a/snow/lodispp/__pycache__/utils.cpython-312.pyc and b/snow/lodispp/__pycache__/utils.cpython-312.pyc differ diff --git a/snow/lodispp/physics.py b/snow/lodispp/physics.py index 9cadb26..75f7250 100644 --- a/snow/lodispp/physics.py +++ b/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 + + + + + + \ No newline at end of file diff --git a/snow/lodispp/pp_io.py b/snow/lodispp/pp_io.py index 891afaf..94f5f62 100644 --- a/snow/lodispp/pp_io.py +++ b/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()) diff --git a/snow/lodispp/utils.py b/snow/lodispp/utils.py index 073eb37..ba2a1f5 100644 --- a/snow/lodispp/utils.py +++ b/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