From 08b7a3110f24d21bd9a13d9e33a01d230536f1a0 Mon Sep 17 00:00:00 2001 From: Giacomo Becatti Date: Mon, 9 Dec 2024 17:09:12 +0100 Subject: [PATCH 1/6] Started work on EAM parser --- .../lodispp/__pycache__/pp_io.cpython-312.pyc | Bin 13488 -> 16025 bytes snow/lodispp/physics.py | 8 +- snow/lodispp/pp_io.py | 69 +++++++++++++++++- 3 files changed, 73 insertions(+), 4 deletions(-) diff --git a/snow/lodispp/__pycache__/pp_io.cpython-312.pyc b/snow/lodispp/__pycache__/pp_io.cpython-312.pyc index 8d04bef38f76fe7e3e6320e5701c0b40b95d3262..722fbd27ae02b7ed83ddddcf1ccf267bb96e5965 100644 GIT binary patch delta 5886 zcmd5=Yiu0V6`t9d*{65?Sijef$MGXxKN2SngoZpYi4#bO69Q?#K$l@>Y;Ur=v)q}* z4`#KY<%e5BYip{;l$1Ix2qP+y);|jLM@vnNkbNAuh zJPPTr&T7wnopbKJ=YHqjbMM*63#s&bUay;h?;n*%N6cP^`2-8O&z{X(`72~jF)}m6 z$gFIOvN79`jm10{<%W2g3PS=-?L&5%R~)jSIEQ)J{wzP_DkwXk?3SIh-KOxz($%@X zR#T|HLG7@u4p=~?3sjU8sJH>|S7>m92CwXSmKpM;Wv^WFEIU*>!3BMjN5oxfT$dto zO^n17NnKQvIwr%(cvz39P)(Gmq${$htB{rCNPI-pbdn6~Nup@oVnWrGxE_(B;&3FY zXriKr*NYc%YC+Z#9!w^p%An=IHsG}EL}>uQc$i(CJmZ2q2mdR7=FH}%H8bytyulvh z%x7u>tN^rr_)(r=Qe28p2`Rho$cY-<>sd7^PDds5Gx&ThUr;OLYX#ktmB3<=a44Px6B3J$Nt$N4U}FhUGgwZ{>Jtg*D!w;F zMpa9YX~M-r#Nr?ieTBkF=zdS8eu!3$s@RPDJRK+-A-Fwgejw(O&u{+QcYj)uVx$Es z(Al;3V9PkcTy0?7zNwLEZP^uAVH$;=+41ug_2(*DmMU77E83$4BfJ1^7*W)Cdac20ROlvYkZ_}tS|{DsQ8+2C?z>y+!$ zj;>z}oV_EfAYGtgx#QrHvu#S4N=}Do`{sm?$||OW^S-hv?mXr`_1B*E&mMWPe(uQ^ znwI^+%M9=8nc8*UA2{c4UGlds``gd?yO;dk%l`FKyDwDL&u%myy3E+TJ=5$*WmTt2 zW=c*UncF}_`x9WG90=H+_?n+(61OMM)c zx{5mW7xYmf7Z0|wsI6FEmAS;V?Q+SEdL~s}z*j97>#r$lE7o5Nqc4TgH#;(V>dZPm zVAk?Je=%PG_{s|S8q5^%8Q1Zx0lxAAz9#cyK2RzajGzgu_&I;Td`a+{_X^E49TV1> zzq7l{SB1~nO7oOG2=pCFKeGotm9OVwhEDU0y$RZOI)XK;wcR@I9rTkc&9vhvTV)<_ zHUi&2^LdAP+*!vtCY`3o*-zUW%)e!G$T{H3Njc2lxH6oNxU!txlMsHSHp|t-1lu`ZJt$W~rwwP@k2nFs_p`4|=N70fRSLxxo;CHX5x!*BEU;n~WBqqR|et z*{A`!)~E&AV$=g|H3C4}j7FgC#u}jO3=wFD(F`tOMF*beOCWG&+qgqg(Di z&A-Avz_c(M8I2u-2;BU{6G-<+qp*sZ!&x6y3eHk-n8492iHG;zS8y`pVbm?cVFaNq z#Hc6|SHg;>;ptp>C}*PI8y(+;GU*ACo}`76H$Bd8;WR%aL7V__Hl0#FuFh+w5jGz3 zX;J)Ir~K|4gdo?L3nP)-G=^0|yDqXFohUt22Q5~%SjA#@yDfI$B1Xu>X;znf-a;_; zvRt`vL9-l0NkpZvLNL-M7%7`Q-r5#YihP_ho(SHS2^Dx6CJ1x2!pt{%o6{J-XjyQ@#qo(rx)`FQ1J*_&PbqOR zs)i*D!>)vYiPpps25+Tw=x44>rjIJjJ zd;?JhL$gdU8dg;TS39`qK^bL4bvV=Igr+MocbbYFmLi%`I3ZfQ2kbVe>U-2=Tpl1q zU5!o$)r^(KBrqK|qrDJ-*H{dP^Gzw5REotPP-BXyphoeru~FzHe*r}uxLmaVf$^|H z1K|ABj>G|CN*HggR2=3MRpH)>Vq?+-K^E(m41VageuLBVB>JV3x zeOLjPH-=(2gXc1AhjeAwbOk2dtONpq?}?nmTNAd_~;X$XCsUH&TsB&;y`J}~?N zIbQ;}onThF&7Q{l`FofBx0~N?{A|m9_JXfzu5Il0yJG zz9NN30pD;X!KH^A58_WJ(lM~33cF+nGR4&~F#>M{@#uu8B@+pi0DkgJkt{aGMirun zdK$XS#QRhubs*es!lMdY1>nH|_hVIypwoyEKuBa9pBYApnjC>`5Z7u{g+9}HL)g+i+%8qhv>c%%=%nz11?yo z3|1ZH(U$%IB6jdKCTO#4>wC#_h#K64MUkvE7g{#j>RI#eEzRj5)R9d{=oU?JJ_g00 zm+sT67$+y7*5b5e?4~d$&!C zl@8S6aiY{QG&~dX{(yr=ULPJFA;4THral!>$d93t2mER>^^2=(7cRE8asoWa$&6Xi z@q;G!;2!F|me7yLdhkoOh2RA*>c$oqRpJ?A*ZLzty5RY%e-AC1)EV}!f$2YXpi`o)itqk}ho(&F z6m*fF;>7nML36@f}YV{_*wPct>3iz-KakY54o<7mETsQIq8E*#(V3Y#vo9HEe`hC?CC4i5?N zOi4IYDgZze6PbZO4HaY@YdP?Ci^nGs`Z7S{0?U!LA>EWI0_j@O9?1-dU4cJC@%p86 zTTGRcQDp~t8|vTz2?q2d$gOZJ%l?IVho73?#Gi6--t} z**!Maj+K~VEv74n_^s@YS#E`aWKNiS;?-kkj-6fm>&{^5bkL+kEvQww>)EZ^79912{D3`~Uy| delta 3668 zcmb7HS!^4}8Qvke6c2G_-J~d6OY0&f>WJjTmg2+}+@wj6DDsgyRc)|q?n>Htd0=)W zOX?Dg0T+F7E)0saFG*XYsEZVcTfI04P@qpO&=gSNhZJi+wLlBQ4{>b+Em9!u|Id;V zNx3Pq6u;S-`LCJz{(q<+m9FdMKSUxS0UmbgmE~Wi1mUmPsXm^{%Y#LDc~?*cl^DVz zscSup9@V3HmprQXvTxDH=Wq;V#YJ&P)^}N0^s6FZlj;Wyr~$RzH4NHk;E9goH zf*n%pbt0e>Mx9ZqMoLlxb(UM|=7Q|QaJ;47PAkiV2iX_gZ=^ZSnwur|L~~<=o!~9P z`XY(Wfa>oNN?kQ}yVL-C12nHROJWq{q#C(cja&@m{?Qz757tE(2aOigIObaJtv(E_jLkl`_7OK_F-$UM{Q@fb~}D{qBZ9AZ^AnF&2foy zyK)7F2y)1KNC1wrkyr?K6HM6eHS!ttwZcEXmiuYzTm!Pzj( zrtGn!INx2-j54i#upF!OTKC%FTEi8sPJ5!_Z*s@qV_)hoz2nc-jM@#F7$*yMv_cH7 zL~Emj+QwrTaV@<~uH5wrx}f&hK1S5IEduu11AzPOLBIocKVY9d1Ss3xfC)Pac+l291NMHvLHhvU5cntRuxuxwN!kYiQ}!_6h@E5vHf+QmwbSY-`>10ynx1zEb%^Q^ z%^_cR$TN2^o9~G1Lf2arX?4aWH+W4E4}))%H~Gye8GUl zp?dWhtMwo!^Q*by#&}*;C{;FKd2`ru%PUF|7M;(aQ8O)1r{s&ewZivPh7GOTJv2_6 z@hI`gXb1aj@apjqcy+vb-f{v}vH`^`*cUp1Qd0=i2xBPeSBgMP)d5PGmgCJSMLN#@ zI5bJlvoD6!xzi|k7@-Nki~v-Z#9Mk6MqE#s^^Hx63I9Emd=I$tx^OSbt`48@J$JMH zarUp_yN6#Ox23)xk6e58+6&uHd}CXF?xytocHsHj;m-R45uds(MSpng+L3D`*N@*I z*UjzNlQ*TK+kvB7J*frn)VV7CfG{Rfj?*t!SkpC^gg_5%-J(j?s_>7+XJh3G3Fhc> z1HUpwAXP2Fsw_f&?VA0HBr+TLQ{E*I=MS8Al{gjp$<{_uGv%cMl{F>np7kTyW9}u+ zsVZNVFo%FJdzjNxq32%ah-_+nE{^<}|Amn9IG*uIdJ)eK>lQHUEPHeOG0!Ig`|Wt5 zjK8C32LjHztjoY#DngNQbOk|2cn!eunya}-@r_o_X%mGVFxAFhql(}WAGmYFUuCCmWiWriOo{XPhN3|~{_?0gvA z`pv|ER}87}VSk@`s~lEJ4J4CaoO|h|g;TPsSjzukS#47CLRWXLUxj{Ux4t#H{C^Cr zIok9N+Zt$C7v!G$GrvA9LoH&NnP^!~uBYXVw7i*WBw@8&%3n_8q+()))bq=x44j5i zq%A{nPd$p^&+^jch;vx(PJd4NX&Fz_09;T1SzYFWHV z9I`|YvxA2whIqjXb7C1vE>|?^4?y#8@HIy{tL{aG_Rg)j0|;-wvL>jfV1 z-z4-q(B-!aPh!636z$LwzytUS+|B2Kqq%~*YG_Z>t1t%k1@jodJue~Tme9ukg7@ diff --git a/snow/lodispp/physics.py b/snow/lodispp/physics.py index 9cadb26..1501643 100644 --- a/snow/lodispp/physics.py +++ b/snow/lodispp/physics.py @@ -1,10 +1,12 @@ -import numpy as np +""" +Contains functions to ocmpute physical properties for the system, such as pressure, potential energy etc. +""" +import numpy as np from snow.lodispp.pp_io import read_rgl from tqdm import tqdm - -def properties(coords: np.ndarray, elements: np.ndarray, pot_file: str, dist_mat: np.ndarray, l_pressure: bool): +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. diff --git a/snow/lodispp/pp_io.py b/snow/lodispp/pp_io.py index 891afaf..682eb07 100644 --- a/snow/lodispp/pp_io.py +++ b/snow/lodispp/pp_io.py @@ -1,9 +1,76 @@ +""" 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: + comment = 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 = [] + + for rho in range(int(n_rho) // 5): + line = pot_file.readline().split() + for i in range(5): + + F_rho.append(float(line[i])) + + for r in range(int(n_r) // 5): + line = pot_file.readline().split() + for i in range(5): + + Z_r.append(float(line[i])) + + for r 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, Z_r, rho_r + + + + + +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() From eeaacdf56525ef5df246fb6d72847abd02cb2c6e Mon Sep 17 00:00:00 2001 From: Giacomo Becatti Date: Mon, 9 Dec 2024 17:19:57 +0100 Subject: [PATCH 2/6] Parser reads parameters --- snow/lodispp/pp_io.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/snow/lodispp/pp_io.py b/snow/lodispp/pp_io.py index 682eb07..1866f72 100644 --- a/snow/lodispp/pp_io.py +++ b/snow/lodispp/pp_io.py @@ -27,31 +27,38 @@ def read_eam(filepot: str) -> dict: lat_param = float(lat_param) n_rho, d_rho, n_r, d_r, r_cut = map(float, pot_file.readline().split()) - + print(d_r) F_rho = [] Z_r = [] rho_r = [] + r = np.arange(0, (n_rho - 1) * d_r, d_r) + print(r) for rho in range(int(n_rho) // 5): line = pot_file.readline().split() for i in range(5): F_rho.append(float(line[i])) - for r in range(int(n_r) // 5): + for _ in range(int(n_r) // 5): line = pot_file.readline().split() for i in range(5): Z_r.append(float(line[i])) - for r in range(int(n_r) // 5): + 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, Z_r, rho_r + return { + "F_rho": F_rho, + "Z_r": Z_r, + "rho_r": rho_r, + "r": r + } From 829c7810218382fe09f0a583add0046ee06bad9d Mon Sep 17 00:00:00 2001 From: Giacomo Becatti Date: Mon, 9 Dec 2024 17:22:28 +0100 Subject: [PATCH 3/6] Fixed typo in r --- snow/lodispp/pp_io.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/snow/lodispp/pp_io.py b/snow/lodispp/pp_io.py index 1866f72..5b5caa9 100644 --- a/snow/lodispp/pp_io.py +++ b/snow/lodispp/pp_io.py @@ -32,9 +32,11 @@ def read_eam(filepot: str) -> dict: Z_r = [] rho_r = [] - r = np.arange(0, (n_rho - 1) * d_r, d_r) + r = np.arange(0, (n_r - 1) * d_r, d_r) + rho = np.arange(0, (n_rho - 1) * d_rho, d_rho) + print(r) - for rho in range(int(n_rho) // 5): + for _ in range(int(n_rho) // 5): line = pot_file.readline().split() for i in range(5): @@ -57,7 +59,9 @@ def read_eam(filepot: str) -> dict: "F_rho": F_rho, "Z_r": Z_r, "rho_r": rho_r, - "r": r + "r": r, + "rho": rho, + "cut_off": r_cut } From 00fe0bafb89418988b68b8b72b21d8210f7e0129 Mon Sep 17 00:00:00 2001 From: Giacomo Becatti Date: Mon, 9 Dec 2024 17:36:28 +0100 Subject: [PATCH 4/6] Removed debug print in pp_io and added computation of pairwise parameters given coords of pair --- snow/lodispp/physics.py | 37 +++++++++++++++++++++++++++++++++++-- snow/lodispp/pp_io.py | 2 -- 2 files changed, 35 insertions(+), 4 deletions(-) diff --git a/snow/lodispp/physics.py b/snow/lodispp/physics.py index 1501643..41e7220 100644 --- a/snow/lodispp/physics.py +++ b/snow/lodispp/physics.py @@ -3,9 +3,9 @@ """ import numpy as np -from snow.lodispp.pp_io import read_rgl +from snow.lodispp.pp_io import read_rgl, read_eam from tqdm import tqdm - +from scipy.spatial.distance import pdist 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. @@ -131,3 +131,36 @@ def properties_rgl(coords: np.ndarray, elements: np.ndarray, pot_file: str, dist +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 = min(range(len(potential["r"])), key=lambda i: abs(potential["r"][i] - r_ij)) + + rho_r = potential["rho_r"][idx_closer_r] + phi_r = potential["Z_r"][idx_closer_r] + + idx_closer_rho = min(range(len(potential["rho_r"])), key=lambda i: abs(potential["rho_r"][i] - rho_r)) + + F_rho = potential["F_rho"][idx_closer_rho] + + return F_rho, phi_r + + + \ No newline at end of file diff --git a/snow/lodispp/pp_io.py b/snow/lodispp/pp_io.py index 5b5caa9..679f047 100644 --- a/snow/lodispp/pp_io.py +++ b/snow/lodispp/pp_io.py @@ -27,7 +27,6 @@ def read_eam(filepot: str) -> dict: lat_param = float(lat_param) n_rho, d_rho, n_r, d_r, r_cut = map(float, pot_file.readline().split()) - print(d_r) F_rho = [] Z_r = [] rho_r = [] @@ -35,7 +34,6 @@ def read_eam(filepot: str) -> dict: r = np.arange(0, (n_r - 1) * d_r, d_r) rho = np.arange(0, (n_rho - 1) * d_rho, d_rho) - print(r) for _ in range(int(n_rho) // 5): line = pot_file.readline().split() for i in range(5): From 1f6b57d7f3abbce22fd18eedcd4d2f8a707df893 Mon Sep 17 00:00:00 2001 From: Giacomo Becatti Date: Mon, 9 Dec 2024 18:28:22 +0100 Subject: [PATCH 5/6] Implemented (to check) EAM potential energy per atom computation --- snow/lodispp/physics.py | 45 ++++++++++++++++++++++++++++++++++++++--- snow/lodispp/pp_io.py | 4 ++-- 2 files changed, 44 insertions(+), 5 deletions(-) diff --git a/snow/lodispp/physics.py b/snow/lodispp/physics.py index 41e7220..3d9168c 100644 --- a/snow/lodispp/physics.py +++ b/snow/lodispp/physics.py @@ -6,6 +6,7 @@ from snow.lodispp.pp_io import read_rgl, read_eam from tqdm import tqdm 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. @@ -151,16 +152,54 @@ def pair_energy_eam(pot_file: str, coords: np.ndarray) -> float: potential = read_eam(pot_file) - idx_closer_r = min(range(len(potential["r"])), key=lambda i: abs(potential["r"][i] - r_ij)) + idx_closer_r = np.searchsorted(potential["r"], r_ij) rho_r = potential["rho_r"][idx_closer_r] phi_r = potential["Z_r"][idx_closer_r] - idx_closer_rho = min(range(len(potential["rho_r"])), key=lambda i: abs(potential["rho_r"][i] - rho_r)) + idx_closer_rho = np.searchsorted(potential["r"], rho_r) F_rho = potential["F_rho"][idx_closer_rho] - return F_rho, phi_r + 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 679f047..043fbfb 100644 --- a/snow/lodispp/pp_io.py +++ b/snow/lodispp/pp_io.py @@ -12,7 +12,7 @@ def read_eam(filepot: str) -> dict: ---------- filepot : str _description_ - + Returns ------- dict @@ -20,7 +20,7 @@ def read_eam(filepot: str) -> dict: """ with open(filepot, "r") as pot_file: - comment = pot_file.readline() + _ = pot_file.readline() atomic_number, mass, lat_param, lat_type = pot_file.readline().split() atomic_number = int(atomic_number) mass = float(mass) From e394693153dc021f871f858b12f6fbda6a1cd7dd Mon Sep 17 00:00:00 2001 From: Giacomo Becatti Date: Tue, 10 Dec 2024 10:43:56 +0100 Subject: [PATCH 6/6] Fixed issues in filename --- .../__pycache__/physics.cpython-312.pyc | Bin 0 -> 7049 bytes .../lodispp/__pycache__/pp_io.cpython-312.pyc | Bin 16025 -> 16276 bytes .../lodispp/__pycache__/utils.cpython-312.pyc | Bin 14166 -> 14419 bytes snow/lodispp/physics.py | 11 +++++++---- snow/lodispp/pp_io.py | 2 +- snow/lodispp/utils.py | 5 +++-- 6 files changed, 11 insertions(+), 7 deletions(-) create mode 100644 snow/lodispp/__pycache__/physics.cpython-312.pyc diff --git a/snow/lodispp/__pycache__/physics.cpython-312.pyc b/snow/lodispp/__pycache__/physics.cpython-312.pyc new file mode 100644 index 0000000000000000000000000000000000000000..34cc083b60e0525ce626d427a57bc07f78067218 GIT binary patch literal 7049 zcmb_gU2q%Mb>0OQe@lP_0fGbx5?YdyOjs0clKQdfADR|LNi-#^R;;*;+|3HRBtZdK zaCa$*$Rb5&It?9r5+KS%U^-JnGo2|s@=U0w@U5J_)SdK!YpPRpD^EF7z4#@8O!p;E zJ@+ma03}LMCfS+2d-mLO&OP_szwe&kcsw)#&#y+lJAXGs5dVvQ%m>%N&n`g{#BG8n zW(l6;Epc+z0$*#~nxJN>gl*PFB0d$jC+xF!Oxxm)gmc!3X?xt2pl4~I9fD=C-85eR z^3JF5U4omZpHkN(c;;x{{kPO?Gmqq+X4@mQ%-Q*7kn=uO2I+hDZwGW+K(}?p67i<5 z(Wk_u!o`v@Gq;?KDlri#MP$TiBDJgtOlo08jzzgRlaj=gASp2c1VxEa76e9Kkrg2^ z$jHml1%{JBM3CiWNf=~Oq9P=f7^n$JL7HD-1SL8|e~sdEA!5;;lECq-G#{^00+)!8 zI;AY}=xa;yFOD5wm?UzR4}Sty+&r1_z}pOO3%p)l(-3jqorb+B z)6p=RUJ!ch0bLK$Rk5`uLhV5PGozmgK3kh_2K|=0K90>^3!Nb0txNDUieaWpwX5!% zmYcXwYM0uHG}AC6qSmMd_ZeAkt?Q{lQzChdscunO3N2@D2Lbg==wJ8DEY} zT0sIwP7OrCmus+QB`5HFs9hY*Z;}?0Nd6AT7;|w0Ls-&0tM|7t1E)44&M{RffE9_WV%&fk z0*9alm@^ESM?Dyl6{+UMMNU~jpap@F_e2@wnKR2Ng=@eek9Qn7Z^DKB-I@fqxMWmd z08mJ=n~nXd_ORN!I@SO~S4A;i^J^;3%_F}6KUD|UL?tqW9=Y(Ipe#$tJz>XuWjPfW zYIX*ybKyx~HE7O)n&Al6k`SH&2<1S-3{m?Up`Q8o@$-M<>7G*}f_Jgx{C)sWr-V4m z=@g0&rmudDfY8^-rtXR^D=b7tQs&_op31mRrm8ftsgp~(bxDrtB&R!&sF37!%gy8P z9)ovOw{SpnKy#y2rM9=Kgx4gTCQDLI%Pf4<*Q&3*j$!oeD z7IrFj4MNO9X-6gdGAPR_5rj|=Y)ebHBM{1wsz3sqMS_VDlNdFLVsT4q!P;d^<$6wElGAqlMh{?dEV{ zNo(%So_Wyh#jG3K%?IGG)qgvf4;Eh7ZtcrXKk&8P9>@qYrZ2$KT6NRqRCjfq(SH&tea*1I;QwBb{m6WZ_;SROv7G3Rnq zh2i|U9pp!GxJc!Hn4`*_{lzJ*^GFGb(VTt9-%%JVUeiJ&U;0OOLI;Y6*WM_NYN4Ur zs(7sQ;`*h}UfMXUy*Q!u zPi#(V{b#qk&*e__ z96nUCYvG~&sU5@RaObZT_+O^=z>nF9~*pr+`R zxf_M0d|Gq&uAbQmhF9s_@JF8ih8mbC4lVrB`iYI+doYP!Uwx-Qete-aZ$X-H^as0t{44bxE zJ5I2R@7pTAH)mXc?;WZGqi<*zaIBN-f-V}s0$?D*V%6jiC(B|q2H=oqZo~ko0WRZW z5`ugfC_WDceGcGvdg6Qo-b1(tBcW%Z*Hb8$RE4eJjdDXYJi@MFC!dYd=Heo!)KJ>R z$CF;DgN7A!2xCa5lM(=x1(B6>>nxm-|L|9l-+2FZ@Gado#R9tugm%Uwl!V4#LK8|t zN4f5nAsH53kVOgN7@-Y7gYb<`!EwK+yC6hH&P{j6_?v7rF2cFaO42*X#0(O&r50k) zaC2kzf+(GbIvE|*b*LbnM)ECan5vkFg&dvq9pL^5N#BM7PIQKFx8~YEexuN;HHEU1 z552)caNFCLJ@?QbxILa9FTPs3bmz_WW83~$vNJpGKyF@hhYPCa9?YKJq5Zkj8r@ZR zU84uGr^@bNVH9eMe$9QbsA%pZ+0zdL9k(t~ZPBIqlWBC# zg>Au+{ufqWP9!*Kg?*lt{;c)N@`8+!DBN^eyzZW|FYFb>l25y7?Q^N!t6g8a_&y`f z%}FR45sPuVxD0ffso$sYo@ZvsYu8th#y1c&3t;i06)p9l0xxhhL#2LH96|*;1Ej-H z=r(-Y#_KNPHrhjGbjJnbaw1(sqR0iE#z}HEj_z!B6c!;Vy^Ja-AqiI!8?J`i8&032 z4a0``i!!Cn z_RSaWw@%kF`~0>ioSniTeegbY)Bx!nw_e)z^nF7(Y@H@X2dqctKm3NET;bKz<(4)$ zk<$p9o?SazoGbmso!EN&hV5R@)}it3?zeKT2f-KCgPW22!I_+`>}@O9@}4|Zae+1l zWUKbor60R0Ezm&4OSpa6^LU?#_;d%$@?w-_AwtWE)QYj$`p_L`sCvg4N~PGCh%cSS zaG~+h+vR91wK629An(G(htORaWef1Ot?$)`i*{U=&>+LT(@3D`R>)uI6v`Gz=$@q8 z@L7UKS~3ET9mJ#Yd*|y3kzbArZ%T_G3%eL0|2-5HD@l@%EEba5brPiie~96)h`?9G zz*od@#j-+@M++mLj;)PV2t=>@e%*Di%b+)9{PcK-hXK~)Dz literal 0 HcmV?d00001 diff --git a/snow/lodispp/__pycache__/pp_io.cpython-312.pyc b/snow/lodispp/__pycache__/pp_io.cpython-312.pyc index 722fbd27ae02b7ed83ddddcf1ccf267bb96e5965..9fcde9dabbc96d62c011c28115c29c627d67f116 100644 GIT binary patch delta 1816 zcmb_cOKclO7@qNZ*LLh(zv4J0xK5fEMoOBf4TYvDp)_d>Y7hjZmY22iZj!D12(z0e zu^hPtAtfSdQl_Ym;1I;6RfPmtZ~zHaTtGsVh;qPMajcMfK)4hMDMG;f;}Y5f7sOcp zcIKaN{`vmL?#AM!lJc!AHy}2$u2WOF{*RSq>=<`eHt;Z}c_k8a6g`?yV2olzK?J-e z02`r2za`TtCB_)8y9_sHxOW)=D@eB-v52$aKElkfj8~(FL`G5DQ+~;xcum@=({I0Q(VPI79mE@9A z@|OZ7uO_`GO`vwvk4!x0!00sF$7m%qFPuhmIN>V%L|+p>bhffV-|vah-y<>|lOl%+ zc)JA|yv2%u#uJh)j;6?T-WFa=5nE)l6d^9QWeZ?+rDv^FesVI=Ky@kMb>wDjNhNA- zN~fz*n**olXHv48K;elOzy|>5iOj%uW{jL}W5W!y6>ZxSrbDh5`nuV|&HkpJ$vKC9 zeGoM$OH=FSwx^|r;$lyEzUpiI-WS{Q#kPI%B}~6^O+H@@gqFmr5?B)W?50myOD-qZ z)pGkL-ss+(*oa?g-}dz`;UAR1nz3w@ad~9(iH+oy$W}}5N8+~f$dYg~*mTJ(575Hq zvlmO-!EXAedxWmIzkK`p z)9l8mQf!-cdLiHS!VJ)#J3koo!#V^IU|@@;MT{9Pkvs#|C;->q4EwB*jL`_jFul>$ zhv(=YUCHOj1F-4VGkR9fS*Gnx=kr7}a;l}9q!ZX301q=Tp&_6lcL%vg^!!4BtTDn3 zdb!(B66J|)xvhr$o?*Oss-irkODzz<_*Il}*pA*ley<75+}cXI=!M>X+)l6e?lXFr zi3|Y1CgD*puvtPT_eV@#XI5J-SJ~7SCyAQX z>zQz>RzAyPHsN*H4HZ-ZSK$?6!KSkLvxZLIWk$L3$^NQP?;ftVhrV+lUgC1O$Oh9= ziDk|i)-`3e$eUcs3^8&Rfr{0mI>}~#>%IV1u2=?@X3EslmXXidc#oagB((ZBYsuW43MK{Qku#FUV$y1Yv`l9s_LJUGbVf0TK~yz zJ+~fyM+!Qr!cG{KuLr-wFL6NAq=C#?3ew2Jce;L{rK7^(9bPcgbS-S&|T7h;S%{f13C;kGSGOZc_ delta 1611 zcmb_cPi)&%7`J1)iR&h5lR-<`tZ52egS0LiDzIMM$pPVBZgyHl!QO8C1J3X%|oT3anzi%m$q%)i&kA*XlaG~rtmg^HYd~}d1lR6a3RC-zms<0KiNobHuwrf$*hA`Xi zB}qJWC;20MBHzy*BWSw9s=Orp)g!dVmS^hThT7c}(#qbo`OQS?YNG##ME^$Op5-t( z9hn~ZsXei*Y{uG`rA=J@lI(gXxpw01)WuU5dN-1}KLuIcyF5XDS4T-z{W5PQ+;?Xczm6B1gKtm!<0S&nosu%#q87h0Hz^BTB1!WK7q~4VlGvwp0?Zacp9!EHY z@BqRv!u<#jA{<7TKsbW%Fv29l6ha@lcGq*evTR8ss78mG68P+Kq;WFU6%pk|HCvRm zdy$L9h6Cicfe#)+`GAQ64m?#nAM6SoL&M_;aRi=9JNj9_7$b(pDEV$rUaXRvdyYOq zMV&^hl7$WTs8%SK!L(h&w>;=WZU$k1frkaj?7U!yTT%ELK~h`q>UAWK&E?H5zzi+$Nzj+x*90y=@`JjYhK zDBPYeEHNvUysGoBDZ`5xBFGp+Xi!DN3c||_sOnZ7J?MSJ)u+L5$jQ7m`3m|Je6+Yu zDjR;;;g#^Po;zv~lpw~t#!8?&<Gz5%{B zXZuCoB`TFHH|Xa-lHny}DiyF@AFzu-u8*?St!!4&$`waZrKfvV!MDpU6{opkwi}P3 zDtdJ?G`zc$w@u@kONQfAJou0WIoK!+f7LP@XvG6mcn@PAL%Qj#G;tIrbqrpoUdHCKFW4)DoOSYWb6w<7!MR zQK**A@O|O`;w9{A0KJ!O;FPT^5=lWMl}`}40m$($YH;h+hH zXKlm*MS{WC;qoo;rtn$8GJ(tOTi}bZ3x?e{NIz`#RFG%jc~6C}o#W1VC%u#6>Fpe6 zkU1BOpe+Qmo@x?-k3AJ-opx^Spq-m6nDo5PbKInR(z6Ra&kdSCdh#|mOOpB-Nbr5} zl7J4m12_y9%=sp!3i&Ow7UXJ2%^<;aVeL%!I&;QzoDl3Ajg5M&2+GugwS z19k^PqzzsUZsS>PhIuLICA~^J;kZMZu2Chc8PZ6~$V%yiq>PM686`yx>nEEYHKf;4 z*v%Mh;W?%IR@viq^_1sI_VA?Gx7J;t)2%}m#3 z>9X2rCWZ$36kMs=#UDlGC+2rm2l&DLj8uRi7t#%(BFyyfW;q`V$I8eFx|SMFuEk=Z zt$;apBd8kF$3}-}8;sXB@h8ysmAP1(BIFqS(Q-gsB1e`u*_JDaZlztI)*T_gfmzpB z%!t;(SgSVz_v@t1Rw&1CT*qrRPLDl`V^>gz5t{2ONgo`pZ{*`>yAR{_S4hMZ8_p8A z6fwzFvs!_)}|ak6th~y z6#Jl|^)NXOXIh&D0~LSbFb=^-t>5quJe;%0d}FI@BZF{l+e!W$+FHTceuIp`O8XxE zJpMlmfsXwIJJ>OB%45V+nTe>8QK+F&2U^$Oo)DARgGbvLNC-K{KX)j{Cq~lvOde(y z>zwUn5<7@%Jz$$hS!R@lrz|GzftJp&a1lj(12|(9Ug!*#)vnt8B(fSQB*9g`L@FP4 slSWdR?6rN%4-Cg#NWnq`)N4mh4d18sx{i3fqOb(537Ta~CZI z1shC=VxXBmtWRQ}+SEkNeW}qX>I>1Bl%l3tON=k3O-!@WVyZsqnM+!cC~=bc?z!hX zbIzG>=H9jVcTw9{R;vlYc;@TgVcCDq=EgJUe8;8{F^*azwtx*Ik~xb4S;S3Dvp8TU zi0DSeh(4eP1`#}p0+^WPY+{k|h?i&j=O9)&huGxYAtT9}6#~xx7duVt5bWADIQv0j z2gK%)e3AnVx=8`ag=an~+;ugyHKZrbSus%XFLCT8Y?pRTFeH#XC|fjvr6M+uYxQDJ zG#rbnD|lY8NjPO1jKzYHP=ba}v+O#TI-;M{VK@6JQxVHSp}nl%^Z|aC+0DNm>_zBh zOTv6#w<+N=-m;INzC*xS|AX@bXJb!W3UM|&Whud>>`hByevN?U4-052W6E+yL}^(!l&1NSEoGYrWwt zNH($l4U}Q^E#m9BMaj>C+5vJsw>UY*I+i+DbIZ?n{31eMrWI5YyJOi?*Jd@`7zvYZk%o7S;bSG{-;*n~RavSH&Hv&jRam3`#) zoB1kQ%At(ib(gtXxygaSR2~PIL3fSnzWi?lzZRPTu5LI&Ye-SXBXV#&8Xg;$B@5*h zXeI|u#Y$FP+;-qFyi4Y7T{?L{GeF9yM5Dpcu#%wrxUWw0CFDpv2n4NXmy27(PHJ4Mu z#ia61+{`A+kKtSF%koM)CpFQx_9!vI=V{T>8-@p?aG6ijt_Qh=5csyXWD6x?k)U=I z12L+k^qdNY>GSM_uTmti~yOGE)uPP4YXI%Ilp!JR-Ns?(LSJV@pZSfgUHk_3e}X+4qiJ5 zelcl3yIq~%90XHx>}>J`cvFV>$xsxOp=CUrE0+=ng7;OUZLF>43Ea)*YpP5MAaBD| zj{Y4^c&?Zcg(54{hL3!Xpe(5~4ITkrIbSWAMXz zk~7*>*S-{fLAM0pCrNp$l=nh;1DfzO2drLA<=>4hkr){tl@HN3zzJd_ float: idx_closer_r = np.searchsorted(potential["r"], r_ij) rho_r = potential["rho_r"][idx_closer_r] - phi_r = potential["Z_r"][idx_closer_r] - - idx_closer_rho = np.searchsorted(potential["r"], rho_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] diff --git a/snow/lodispp/pp_io.py b/snow/lodispp/pp_io.py index 043fbfb..94f5f62 100644 --- a/snow/lodispp/pp_io.py +++ b/snow/lodispp/pp_io.py @@ -236,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