diff --git a/jarvis/__init__.py b/jarvis/__init__.py index afffb341..c6964e76 100644 --- a/jarvis/__init__.py +++ b/jarvis/__init__.py @@ -1,6 +1,6 @@ """Version number.""" -__version__ = "2024.4.10" +__version__ = "2024.4.20" import os diff --git a/jarvis/analysis/defects/surface.py b/jarvis/analysis/defects/surface.py index 6290e6e4..224e9f0f 100644 --- a/jarvis/analysis/defects/surface.py +++ b/jarvis/analysis/defects/surface.py @@ -35,10 +35,11 @@ def __init__( atoms=None, indices=[0, 0, 1], layers=3, - thickness=25, + thickness=None, vacuum=18.0, tol=1e-10, from_conventional_structure=True, + use_thickness_c=True, ): """Initialize the class. @@ -67,6 +68,7 @@ def __init__( self.vacuum = vacuum self.layers = layers self.thickness = thickness + self.use_thickness_c = use_thickness_c # Note thickness overwrites layers def to_dict(self): @@ -147,7 +149,26 @@ def make_surface(self): cartesian=True, ) if self.thickness is not None and (self.thickness) > 0: - self.layers = int(self.thickness / new_atoms.lattice.c) + 1 + if not self.use_thickness_c: + new_lat = new_atoms.lattice_mat # lat_lengths() + a1 = new_lat[0] + a2 = new_lat[1] + a3 = new_lat[2] + new_lat = np.array( + [ + a1, + a2, + np.cross(a1, a2) + * np.dot(a3, np.cross(a1, a2)) + / norm(np.cross(a1, a2)) ** 2, + ] + ) + + a3 = new_lat[2] + self.layers = int(self.thickness / np.linalg.norm(a3)) + 1 + else: + self.layers = int(self.thickness / new_atoms.lattice.c) + 1 + # print("self.layers", self.layers) # dims=get_supercell_dims(new_atoms,enforce_c_size=self.thickness) # print ('dims=',dims,self.layers) # surf_atoms = new_atoms.make_supercell_matrix([1, 1, dims[2]]) diff --git a/jarvis/core/atoms.py b/jarvis/core/atoms.py index 722e7734..53633158 100644 --- a/jarvis/core/atoms.py +++ b/jarvis/core/atoms.py @@ -20,10 +20,26 @@ import string import datetime from collections import defaultdict +from sklearn.metrics import mean_absolute_error +import zipfile +import json amu_gm = 1.66054e-24 ang_cm = 1e-8 +with zipfile.ZipFile( + str( + os.path.join( + os.path.dirname(__file__), "mineral_name_prototype.json.zip" + ) + ), + "r", +) as z: + # Open the specific JSON file within the zip + with z.open("mineral_name_prototype.json") as file: + # Load the JSON data from the file + mineral_json_file = json.load(file) + class Atoms(object): """ @@ -649,10 +665,10 @@ def check_polar(self): dn = dn + Specie(element).Z polar = False if up != dn: - print("Seems like a polar materials.") + # print("Seems like a polar materials.") polar = True if up == dn: - print("Non-polar") + # print("Non-polar") polar = False return polar @@ -1039,6 +1055,106 @@ def density(self): ) return den + def plot_atoms( + self=None, + colors=[], + sizes=[], + cutoff=1.9, + opacity=0.5, + bond_width=2, + filename=None, + ): + """Plot atoms using plotly.""" + import plotly.graph_objects as go + + fig = go.Figure() + if not colors: + colors = ["blue", "green", "red"] + + unique_elements = self.uniq_species + if len(unique_elements) > len(colors): + raise ValueError("Provide more colors.") + color_map = {} + size_map = {} + for ii, i in enumerate(unique_elements): + color_map[i] = colors[ii] + size_map[i] = Specie(i).Z * 2 + cart_coords = self.cart_coords + elements = self.elements + atoms_arr = [] + + for ii, i in enumerate(cart_coords): + atoms_arr.append( + [ + i[0], + i[1], + i[2], + color_map[elements[ii]], + size_map[elements[ii]], + ] + ) + # atoms = [ + # (0, 0, 0, 'red'), # Atom 1 + # (1, 1, 1, 'blue'), # Atom 2 + # (2, 0, 1, 'green') # Atom 3 + # ] + + # Create a scatter plot for the 3D points + trace1 = go.Scatter3d( + x=[atom[0] for atom in atoms_arr], + y=[atom[1] for atom in atoms_arr], + z=[atom[2] for atom in atoms_arr], + mode="markers", + marker=dict( + size=[atom[4] for atom in atoms_arr], # Marker size + color=[atom[3] for atom in atoms_arr], # Marker color + opacity=opacity, + ), + ) + fig.add_trace(trace1) + + # Update plot layout + fig.update_layout( + title="3D Atom Coordinates", + scene=dict( + xaxis_title="X Coordinates", + yaxis_title="Y Coordinates", + zaxis_title="Z Coordinates", + ), + margin=dict(l=0, r=0, b=0, t=0), # Tight layout + ) + if bond_width is not None: + nbs = self.get_all_neighbors(r=5) + bonds = [] + for i in nbs: + for j in i: + if j[2] <= cutoff: + bonds.append([j[0], j[1]]) + for bond in bonds: + # print(bond) + # Extract coordinates of the first and second atom in each bond + x_coords = [atoms_arr[bond[0]][0], atoms_arr[bond[1]][0]] + y_coords = [atoms_arr[bond[0]][1], atoms_arr[bond[1]][1]] + z_coords = [atoms_arr[bond[0]][2], atoms_arr[bond[1]][2]] + + fig.add_trace( + go.Scatter3d( + x=x_coords, + y=y_coords, + z=z_coords, + mode="lines", + line=dict(color="grey", width=bond_width), + marker=dict( + size=0.1 + ), # Small marker size to make lines prominent + ) + ) + # Show the plot + if filename is not None: + fig.write_image(filename) + else: + fig.show() + @property def atomic_numbers(self): """Get list of atomic numbers of atoms in the atoms object.""" @@ -1168,6 +1284,105 @@ def packing_fraction(self): pf = np.array([4 * np.pi * total_rad / (3 * self.volume)]) return round(pf[0], 5) + def get_alignn_feats( + self, + model_name="jv_formation_energy_peratom_alignn", + max_neighbors=12, + neighbor_strategy="k-nearest", + use_canonize=True, + atom_features="cgcnn", + line_graph=True, + cutoff=8, + model="", + ): + """Get ALIGNN features.""" + + def get_val(model, g, lg): + activation = {} + + def getActivation(name): + # the hook signature + def hook(model, input, output): + activation[name] = output.detach() + + return hook + + h = model.readout.register_forward_hook(getActivation("readout")) + out = model([g, lg]) + del out + h.remove() + return activation["readout"][0] + + from alignn.graphs import Graph + from alignn.pretrained import get_figshare_model + + g, lg = Graph.atom_dgl_multigraph( + self, + cutoff=cutoff, + atom_features=atom_features, + max_neighbors=max_neighbors, + neighbor_strategy=neighbor_strategy, + compute_line_graph=line_graph, + use_canonize=use_canonize, + ) + if model == "": + model = get_figshare_model( + model_name="jv_formation_energy_peratom_alignn" + ) + h = get_val(model, g, lg) + return h + + def get_mineral_prototype_name( + self, prim=True, include_c_over_a=False, digits=3 + ): + from jarvis.analysis.structure.spacegroup import Spacegroup3D + + spg = Spacegroup3D(self) + number = spg.space_group_number + cnv_atoms = spg.conventional_standard_structure + n_conv = cnv_atoms.num_atoms + # hall_number=str(spg._dataset['hall_number']) + wyc = "".join(list((sorted(set(spg._dataset["wyckoffs"]))))) + name = ( + (self.composition.prototype_new) + + "_" + + str(number) + + "_" + + str(wyc) + + "_" + + str(n_conv) + ) + # if include_com: + if include_c_over_a: + # print(cnv_atoms) + abc = cnv_atoms.lattice.abc + ca = round(abc[2] / abc[0], digits) + # com_positions = "_".join(map(str,self.get_center_of_mass())) + # round(np.sum(spg._dataset['std_positions']),round_digit) + name += "_" + str(ca) + # name+="_"+str(com_positions) + return name + + def get_minaral_name(self, model=""): + """Get mineral prototype.""" + mae = np.inf + feats = self.get_alignn_feats(model=model) + nm = self.get_mineral_prototype_name() + if nm in mineral_json_file: + for i in mineral_json_file[nm]: + maem = mean_absolute_error(i[1], feats) + if maem < mae: + mae = maem + name = i[0] + else: + for i, j in mineral_json_file.items(): + for k in j: + maem = mean_absolute_error(k[1], feats) + if maem < mae: + mae = maem + name = k[0] + return name + def lattice_points_in_supercell(self, supercell_matrix): """ Adapted from Pymatgen. @@ -1236,6 +1451,7 @@ def describe( """Describe for NLP applications.""" from jarvis.analysis.diffraction.xrd import XRD + min_name = self.get_minaral_name() if include_spg: from jarvis.analysis.structure.spacegroup import Spacegroup3D @@ -1264,6 +1480,7 @@ def describe( fracs[i] = round(j, 3) info = {} chem_info = { + "mineral_name": min_name, "atomic_formula": self.composition.reduced_formula, "prototype": self.composition.prototype, "molecular_weight": round(self.composition.weight / 2, 2), @@ -1372,8 +1589,30 @@ def describe( + struct_info["wyckoff"] + "." ) + if min_name is not None: + line3 = ( + chem_info["atomic_formula"] + + " is " + + min_name + + "-derived structured and" + + " crystallizes in the " + + struct_info["crystal_system"] + + " " + + str(struct_info["spg_symbol"]) + + " spacegroup." + ) + else: + line3 = ( + chem_info["atomic_formula"] + + " crystallizes in the " + + struct_info["crystal_system"] + + " " + + str(struct_info["spg_symbol"]) + + " spacegroup." + ) info["desc_1"] = line1 info["desc_2"] = line2 + info["desc_3"] = line3 return info def make_supercell_matrix(self, scaling_matrix): diff --git a/jarvis/core/composition.py b/jarvis/core/composition.py index 4965d561..223adb11 100644 --- a/jarvis/core/composition.py +++ b/jarvis/core/composition.py @@ -1,4 +1,5 @@ """Modules handling chemical composition.""" + # from math import gcd import string from jarvis.core.specie import Specie @@ -76,9 +77,36 @@ def prototype(self): reduced, repeat = self.reduce() N = 0 for specie, count in reduced.items(): - proto = proto + str(all_upper[N]) + str(round(count, 3)) - N = N + 1 - return proto.replace("1", "") + if count != 1: + proto = proto + str(all_upper[N]) + str(round(count, 3)) + N = N + 1 + else: + proto = proto + str(all_upper[N]) + N = N + 1 + return proto # .replace("1", "") + + @property + def prototype_new(self): + """Get chemical prototypes such as A, AB etc.""" + proto = "" + all_upper = string.ascii_uppercase + # print('reduce',self.reduce()) + reduced, repeat = self.reduce() + items = sorted(list(reduced.values()), reverse=True) + # print('items',items) + N = 0 + for c in items: + if c == 1: + proto = proto + str(all_upper[N]) + N = N + 1 + else: + proto = proto + str(all_upper[N]) + str(round(c, 3)) + N = N + 1 + + # for specie, count in reduced.items(): + # proto = proto + str(all_upper[N]) + str(round(count, 3)) + # N = N + 1 + return proto # .replace("1", "") def to_dict(self): """Return dictionary format.""" diff --git a/jarvis/core/mineral_name_prototype.json.zip b/jarvis/core/mineral_name_prototype.json.zip new file mode 100644 index 00000000..56bbffcb Binary files /dev/null and b/jarvis/core/mineral_name_prototype.json.zip differ diff --git a/jarvis/core/utils.py b/jarvis/core/utils.py index a777f3dd..7fde5013 100644 --- a/jarvis/core/utils.py +++ b/jarvis/core/utils.py @@ -249,16 +249,13 @@ def operate_affine(cart_coord=[], affine_matrix=[]): def gaussian(x, sigma): """Get Gaussian profile.""" - return np.exp(-(x ** 2) / (2 * sigma ** 2)) + return np.exp(-(x**2) / (2 * sigma**2)) def lorentzian2(x, gamma): """Get Lorentziann profile.""" return ( - gamma - / 2 - / (np.pi * (x ** 2 + (gamma / 2) ** 2)) - / (2 / (np.pi * gamma)) + gamma / 2 / (np.pi * (x**2 + (gamma / 2) ** 2)) / (2 / (np.pi * gamma)) ) @@ -275,7 +272,14 @@ def digitize_array(values=[], max_len=10): def bond_angle( - dist1, dist2, bondx1, bondx2, bondy1, bondy2, bondz1, bondz2, + dist1, + dist2, + bondx1, + bondx2, + bondy1, + bondy2, + bondz1, + bondz2, ): """Get an angle.""" nm = dist1 * dist2 @@ -324,7 +328,7 @@ def volumetric_grid_reshape(data=[], final_grid=[50, 50, 50]): def cos_formula(a, b, c): """Get angle between three edges for oblique triangles.""" - res = (a ** 2 + b ** 2 - c ** 2) / (2 * a * b) + res = (a**2 + b**2 - c**2) / (2 * a * b) res = -1.0 if res < -1.0 else res res = 1.0 if res > 1.0 else res return np.arccos(res) diff --git a/jarvis/io/vasp/outputs.py b/jarvis/io/vasp/outputs.py index ed15e329..5f639d04 100644 --- a/jarvis/io/vasp/outputs.py +++ b/jarvis/io/vasp/outputs.py @@ -199,7 +199,6 @@ def read_file(self, lines=""): def chg_set(self, text, start, end, volume, ng): """Return CHGCAR sets.""" - lines_0 = text[start:end] # tmp = np.empty((ng)) # p = np.fromstring('\n'.join(lines_0), sep=' ') @@ -231,7 +230,6 @@ def vac_potential( plot=True, ): """Calculate vacuum potential used in work-function calculation.""" - if use_ase: from ase.calculators.vasp import VaspChargeDensity @@ -514,6 +512,80 @@ def efermi(self): efermi.append(float(i.split()[2])) return efermi[-1] + def all_structures(self, elements=[]): + """Get all structure snapshots.""" + force_pattern = "TOTAL-FORCE (eV/Angst)" + if not elements: + try: + atoms = Atoms.from_poscar( + self.filename.replace("OUTCAR", "POSCAR") + ).elements + except Exception: + print("Check POSCAR exsist, or pass elements.") + pass + blocks = [] + for ii, i in enumerate(self.data): + if " free energy TOTEN =" in i: + blocks.append(i) + if " free energy ML TOTEN =" in i: + blocks.append(i) + nions = self.nions + atoms_array = [] + energy_array = [] + force_array = [] + stress_array = [] + for ii, i in enumerate(self.data): + if "in kB" in i: + stress = [ + float(j) for j in self.data[ii].split("in kB")[1].split() + ] + stress_array.append(stress) + if "VOLUME and BASIS-vectors are now :" in i: + tmp1 = [float(j) for j in self.data[ii + 5].split()] + tmp2 = [float(j) for j in self.data[ii + 6].split()] + tmp3 = [float(j) for j in self.data[ii + 7].split()] + lat_mat = [ + [tmp1[0], tmp1[1], tmp1[2]], + [tmp2[0], tmp2[1], tmp2[2]], + [tmp3[0], tmp3[1], tmp3[2]], + ] + if " free energy TOTEN =" in i: + energy = ( + self.data[ii] + .split(" free energy TOTEN =")[1] + .split("eV")[0] + ) + energy_array.append(energy) + if " free energy ML TOTEN =" in i: + energy = ( + self.data[ii] + .split(" free energy TOTEN =")[1] + .split("eV")[0] + ) + energy_array.append(energy) + + if force_pattern in i: + coords = [] + forces = [] + for j in range(nions): + tmp = [float(k) for k in self.data[ii + 2 + j].split()] + coords.append([tmp[0], tmp[1], tmp[2]]) + forces.append([tmp[3], tmp[4], tmp[5]]) + # print(tmp) + atoms = Atoms( + lattice_mat=lat_mat, + coords=coords, + elements=elements, + cartesian=True, + ) + + atoms_array.append(atoms) + if len(blocks) != len(atoms_array): + print( + "WARNING: check OUTCAR parser", len(blocks), len(atoms_array) + ) + return atoms_array, energy_array, force_array, stress_array + @property def all_band_energies(self): """Get all band energies.""" @@ -2377,10 +2449,3 @@ def parse_raman_dat( info["indices"] = indices info["intensity"] = intensity return info - - -""" -kp='/users/knc6/Software/Devs/jarvis/jarvis/examples/vasp/SiOptb88/MAIN-RELAX-bulk@mp_149/KPOINT' -kpt=Kpoints(filename=kp) -print (kpt) -""" diff --git a/jarvis/tests/testfiles/core/test_atoms.py b/jarvis/tests/testfiles/core/test_atoms.py index 0d40f4d0..d044e1c3 100644 --- a/jarvis/tests/testfiles/core/test_atoms.py +++ b/jarvis/tests/testfiles/core/test_atoms.py @@ -96,7 +96,6 @@ def test_from_cif(): def test_basic_atoms(): - Si = Atoms( lattice_mat=FIXTURES["lattice_mat"], coords=FIXTURES["coords"], @@ -112,6 +111,7 @@ def test_basic_atoms(): print(opt.from_optimade(opt_info)) polar = Si.check_polar + # prot = Si.get_prototype_name() Si.props = ["a", "a"] vac_pad = VacuumPadding(Si) den_2d = round(vac_pad.get_effective_2d_slab().density, 2) @@ -119,7 +119,7 @@ def test_basic_atoms(): den_lll_red = round(Si.get_lll_reduced_structure().density, 2) strng = Si.get_string() scell_nat_old = Si.make_supercell_old([2, 2, 2]).num_atoms - descr = Si.describe() + # descr = Si.describe() scell_nat = Si.make_supercell([2, 2, 2]).num_atoms scell_nat2 = Si.make_supercell_matrix( [[2, 0, 0], [0, 2, 0], [0, 0, 2]] @@ -266,4 +266,5 @@ def test_remove_sites_by_indices(): assert Si2_supercell_without_two_atoms.num_atoms == 14 +# test_basic_atoms() # test_remove_sites_by_indices() diff --git a/jarvis/tests/testfiles/core/test_composition.py b/jarvis/tests/testfiles/core/test_composition.py index 9571536b..bb185687 100644 --- a/jarvis/tests/testfiles/core/test_composition.py +++ b/jarvis/tests/testfiles/core/test_composition.py @@ -6,12 +6,14 @@ def test_comp(): cc = Composition(comp) assert ( - cc.prototype, + cc.prototype_new, cc.formula, cc.reduced_formula, round(cc.weight, 4), cc.to_dict(), - ) == ("AB2", "Li2O4", "LiO2", 77.8796, comp) + ) == ("A2B", "Li2O4", "LiO2", 77.8796, comp) + # ) == ("AB2", "Li2O4", "LiO2", 77.8796, comp) + assert cc.prototype == "AB2" c = Composition.from_string("Al2O3Al5Co6O1") td = c.to_dict() fd = Composition.from_dict(td) diff --git a/jarvis/tests/testfiles/io/vasp/test_outputs.py b/jarvis/tests/testfiles/io/vasp/test_outputs.py index 03e3696f..df885a3a 100644 --- a/jarvis/tests/testfiles/io/vasp/test_outputs.py +++ b/jarvis/tests/testfiles/io/vasp/test_outputs.py @@ -196,7 +196,7 @@ "OUTCAR", ) ) -#TODO: +# TODO: # wder = Waveder( # os.path.join( # os.path.dirname(__file__), @@ -266,16 +266,16 @@ def test_locpot(): (2, 56, 56, 56), ) vac = loc.vac_potential()[0] - #assert round(vac, 2) == round(7.62302803577618, 2) + # assert round(vac, 2) == round(7.62302803577618, 2) - #td = loc.to_dict() - #fd = Locpot.from_dict(td) + # td = loc.to_dict() + # fd = Locpot.from_dict(td) vac = loc.vac_potential(direction="Y")[0] - #assert round(vac, 2) == round(7.62302803577618, 2) + # assert round(vac, 2) == round(7.62302803577618, 2) vac = loc.vac_potential(direction="Z")[0] - #assert round(vac, 2) == round(7.62302803577618, 2) + # assert round(vac, 2) == round(7.62302803577618, 2) def test_vrun(): @@ -329,6 +329,7 @@ def test_out(): out_efg = Outcar( os.path.join(os.path.dirname(__file__), "OUTCAR.EFG-JVASP-12148") ) + print(out_efg.all_structures()) rl, imag = opt_out.freq_dielectric_tensor() nedos = opt_out.nedos out_efg_raw = Outcar( diff --git a/setup.py b/setup.py index 4a03d57a..bf1d1c99 100644 --- a/setup.py +++ b/setup.py @@ -12,7 +12,7 @@ setup( name="jarvis-tools", - version="2024.4.10", + version="2024.4.20", long_description=long_d, install_requires=[ "numpy>=1.20.1", @@ -34,6 +34,7 @@ "magpie.json", "element_charge.json", "atom_init.json", + "mineral_name_prototype.json.zip", ], "jarvis.tasks.lammps.templates": [ "displace.mod",