From 6e9fc6819b9386a8df866c51e992a362fb9c4277 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Sun, 14 Apr 2024 17:26:21 -0500 Subject: [PATCH 1/6] Implement Fitsnap interface for SNAP and ACE --- structuretoolkit/analyse/__init__.py | 9 + structuretoolkit/analyse/fitsnap.py | 263 +++++++++++++++++++++++++++ tests/test_fitsnap.py | 67 +++++++ 3 files changed, 339 insertions(+) create mode 100644 structuretoolkit/analyse/fitsnap.py create mode 100644 tests/test_fitsnap.py diff --git a/structuretoolkit/analyse/__init__.py b/structuretoolkit/analyse/__init__.py index dd13bacc0..569bcefc7 100644 --- a/structuretoolkit/analyse/__init__.py +++ b/structuretoolkit/analyse/__init__.py @@ -29,6 +29,15 @@ get_snap_descriptor_derivatives, ) +try: + from structuretoolkit.analyse.fitsnap import ( + get_snap_descriptor_derivatives as get_snap_descriptor_derivatives_fitsnap, + get_ace_descriptor_derivatives as get_ace_descriptor_derivatives_fitsnap, + get_ace_descriptor_derivatives, + ) +except ImportError: + pass + def get_symmetry( structure, use_magmoms=False, use_elements=True, symprec=1e-5, angle_tolerance=-1.0 diff --git a/structuretoolkit/analyse/fitsnap.py b/structuretoolkit/analyse/fitsnap.py new file mode 100644 index 000000000..87b8ea0f1 --- /dev/null +++ b/structuretoolkit/analyse/fitsnap.py @@ -0,0 +1,263 @@ +from ase.atoms import Atoms +from typing import Optional, Union +import random +import numpy as np +from fitsnap3lib.fitsnap import FitSnap + + +def get_ace_descriptor_derivatives( + structure: Atoms, + atom_types: list[str], + ranks: list[int] = [1, 2, 3, 4, 5, 6], + lmax: list[int] = [1, 2, 2, 2, 1, 1], + nmax: list[int] = [22, 2, 2, 2, 1, 1], + nmaxbase: int = 22, + rcutfac: float = 4.604694451, + lambda_value: float = 3.059235105, + lmin: list[int] = [1, 1, 1, 1, 1, 1], + bzeroflag: bool = True, + cutoff: float = 10.0, +) -> np.ndarray: + """ + Calculate per atom ACE descriptors using FitSNAP https://fitsnap.github.io + + Args: + structure (ase.atoms.Atoms): atomistic structure as ASE atoms object + atom_types (list[str]): list of element types + ranks (list): + lmax (list): + nmax (list): + nmaxbase (int): + rcutfac (float): + lambda_value (float): + lmin (list): + cutoff (float): cutoff radius for the construction of the neighbor list + + Returns: + np.ndarray: Numpy array with the calculated descriptor derivatives + """ + settings = { + "ACE": { + "numTypes": len(atom_types), + "ranks": " ".join([str(r) for r in ranks]), + "lmax": " ".join([str(l) for l in lmax]), + "nmax": " ".join([str(n) for n in nmax]), + "nmaxbase": nmaxbase, + "rcutfac": rcutfac, + "lambda": lambda_value, + "type": " ".join(atom_types), + "lmin": " ".join([str(l) for l in lmin]), + "bzeroflag": True, + "bikflag": True, + }, + "CALCULATOR": { + "calculator": "LAMMPSPACE", + "energy": 1, + "force": 1, + "stress": 0, + }, + "REFERENCE": { + "units": "metal", + "atom_style": "atomic", + "pair_style": "zero " + str(cutoff), + "pair_coeff": "* *", + }, + } + fs = FitSnap(settings, comm=None, arglist=["--overwrite"]) + a, b, w = fs.calculator.process_single(_ase_scraper(data=[structure])[0]) + return a + + +def get_snap_descriptor_derivatives( + structure: Atoms, + atom_types: list[str], + twojmax: int = 6, + element_radius: list[int] = [4.0], + rcutfac: float = 1.0, + rfac0: float = 0.99363, + rmin0: float = 0.0, + bzeroflag: bool = False, + quadraticflag: bool = False, + weights: Optional[Union[list, np.ndarray]] = None, + cutoff: float = 10.0, +) -> np.ndarray: + """ + Calculate per atom SNAP descriptors using FitSNAP https://fitsnap.github.io + + Args: + structure (ase.atoms.Atoms): atomistic structure as ASE atoms object + atom_types (list[str]): list of element types + twojmax (int): band limit for bispectrum components (non-negative integer) + element_radius (list[int]): list of radii for the individual elements + rcutfac (float): scale factor applied to all cutoff radii (positive real) + rfac0 (float): parameter in distance to angle conversion (0 < rcutfac < 1) + rmin0 (float): parameter in distance to angle conversion (distance units) + bzeroflag (bool): subtract B0 + quadraticflag (bool): generate quadratic terms + weights (list/np.ndarry/None): list of neighbor weights, one for each type + cutoff (float): cutoff radius for the construction of the neighbor list + + Returns: + np.ndarray: Numpy array with the calculated descriptor derivatives + """ + if weights is None: + weights = [1.0] * len(atom_types) + settings = { + "BISPECTRUM": { + "numTypes": len(atom_types), + "twojmax": twojmax, + "rcutfac": rcutfac, + "rfac0": rfac0, + "rmin0": rmin0, + "wj": " ".join([str(w) for w in weights]), + "radelem": " ".join([str(r) for r in element_radius]), + "type": " ".join(atom_types), + "wselfallflag": 0, + "chemflag": 0, + "bzeroflag": bzeroflag, + "quadraticflag": quadraticflag, + }, + "CALCULATOR": { + "calculator": "LAMMPSSNAP", + "energy": 1, + "force": 1, + "stress": 0, + }, + "REFERENCE": { + "units": "metal", + "atom_style": "atomic", + "pair_style": "zero " + str(cutoff), + "pair_coeff": "* *", + }, + } + fs = FitSnap(settings, comm=None, arglist=["--overwrite"]) + a, b, w = fs.calculator.process_single(_ase_scraper(data=[structure])[0]) + return a + + +def _assign_validation(group_table): + """ + Given a dictionary of group info, add another key for test bools. + + Args: + group_table: Dictionary of group names. Must have keys "nconfigs" and "testing_size". + + Modifies the dictionary in place by adding another key "test_bools". + """ + + for name in group_table: + nconfigs = group_table[name]["nconfigs"] + assert "testing_size" in group_table[name] + assert group_table[name]["testing_size"] <= 1.0 + test_bools = [ + random.random() < group_table[name]["testing_size"] + for i in range(0, nconfigs) + ] + + group_table[name]["test_bools"] = test_bools + + +def _ase_scraper(data) -> list: + """ + Function to organize groups and allocate shared arrays used in Calculator. For now when using + ASE frames, we don't have groups. + + Args: + s: fitsnap instance. + data: List of ASE frames or dictionary group table containing frames. + + Returns a list of data dictionaries suitable for fitsnap descriptor calculator. + If running in parallel, this list will be distributed over procs, so that each proc will have a + portion of the list. + """ + + # Simply collate data from Atoms objects if we have a list of Atoms objects. + if type(data) == list: + # s.data = [collate_data(atoms) for atoms in data] + return [_collate_data(atoms) for atoms in data] + # If we have a dictionary, assume we are dealing with groups. + elif type(data) == dict: + _assign_validation(data) + # s.data = [] + ret = [] + for name in data: + frames = data[name]["frames"] + # Extend the fitsnap data list with this group. + # s.data.extend([collate_data(atoms, name, data[name]) for atoms in frames]) + ret.extend([_collate_data(atoms, name, data[name]) for atoms in frames]) + return ret + else: + raise Exception("Argument must be list or dictionary for ASE scraper.") + + +def _get_apre(cell): + """ + Calculate transformed ASE cell for LAMMPS calculations. Thank you Jan Janssen! + + Args: + cell: ASE atoms cell. + + Returns transformed cell as np.array which is suitable for LAMMPS. + """ + a, b, c = cell + an, bn, cn = [np.linalg.norm(v) for v in cell] + + alpha = np.arccos(np.dot(b, c) / (bn * cn)) + beta = np.arccos(np.dot(a, c) / (an * cn)) + gamma = np.arccos(np.dot(a, b) / (an * bn)) + + xhi = an + xyp = np.cos(gamma) * bn + yhi = np.sin(gamma) * bn + xzp = np.cos(beta) * cn + yzp = (bn * cn * np.cos(alpha) - xyp * xzp) / yhi + zhi = np.sqrt(cn**2 - xzp**2 - yzp**2) + + return np.array(((xhi, 0, 0), (xyp, yhi, 0), (xzp, yzp, zhi))) + + +def _collate_data(atoms, name: str = None, group_dict: dict = None) -> dict: + """ + Function to organize fitting data for FitSNAP from ASE atoms objects. + + Args: + atoms: ASE atoms object for a single configuration of atoms. + name: Optional name of this configuration. + group_dict: Optional dictionary containing group information. + + Returns a data dictionary for a single configuration. + """ + + # Transform ASE cell to be appropriate for LAMMPS. + apre = _get_apre(cell=atoms.cell) + R = np.dot(np.linalg.inv(atoms.cell), apre) + positions = np.matmul(atoms.get_positions(), R) + cell = apre.T + + # Make a data dictionary for this config. + + data = {} + data["Group"] = name # 'ASE' # TODO: Make this customizable for ASE groups. + data["File"] = None + data["Stress"] = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] + data["Positions"] = positions + data["Energy"] = 0.0 + data["AtomTypes"] = atoms.get_chemical_symbols() + data["NumAtoms"] = len(atoms) + data["Forces"] = np.array([0.0, 0.0, 0.0] * len(atoms)) + data["QMLattice"] = cell + data["test_bool"] = 0 + data["Lattice"] = cell + data["Rotation"] = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + data["Translation"] = np.zeros((len(atoms), 3)) + # Inject the weights. + if group_dict is not None: + data["eweight"] = group_dict["eweight"] if "eweight" in group_dict else 1.0 + data["fweight"] = group_dict["fweight"] if "fweight" in group_dict else 1.0 + data["vweight"] = group_dict["vweight"] if "vweight" in group_dict else 1.0 + else: + data["eweight"] = 1.0 + data["fweight"] = 1.0 + data["vweight"] = 1.0 + + return data diff --git a/tests/test_fitsnap.py b/tests/test_fitsnap.py new file mode 100644 index 000000000..812ab2b43 --- /dev/null +++ b/tests/test_fitsnap.py @@ -0,0 +1,67 @@ +from ase.build.bulk import bulk +import structuretoolkit as stk +import unittest + + +try: + from structuretoolkit.analyse.fitsnap import ( + get_ace_descriptor_derivatives, + get_snap_descriptor_derivatives, + ) + + skip_snap_test = False +except ImportError: + skip_snap_test = True + + +@unittest.skipIf( + skip_snap_test, "LAMMPS is not installed, so the SNAP tests are skipped." +) +class TestSNAP(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.structure = bulk("Cu", cubic=True) + cls.numtypes = 1 + cls.twojmax = 6 + cls.rcutfac = 1.0 + cls.rfac0 = 0.99363 + cls.rmin0 = 0.0 + cls.bzeroflag = False + cls.quadraticflag = False + cls.radelem = [4.0] + cls.type = ['Cu'] + cls.wj = [1.0] + + def test_get_snap_descriptor_derivatives(self): + n_coeff = len(stk.analyse.get_snap_descriptor_names( + twojmax=self.twojmax + )) + mat_a = get_snap_descriptor_derivatives( + structure=self.structure, + atom_types=self.type, + twojmax=self.twojmax, + element_radius=self.radelem, + rcutfac=self.rcutfac, + rfac0=self.rfac0, + rmin0=self.rmin0, + bzeroflag=self.bzeroflag, + quadraticflag=self.quadraticflag, + weights=self.wj, + cutoff=10.0, + ) + self.assertEqual(mat_a.shape, (len(self.structure) * 3 + 7, n_coeff + 1)) + + def test_get_ace_descriptor_derivatives(self): + mat_a = get_ace_descriptor_derivatives( + structure=self.structure, + atom_types=self.type, + ranks=[1, 2, 3, 4, 5, 6], + lmax=[1, 2, 2, 2, 1, 1], + nmax=[22, 2, 2, 2, 1, 1], + nmaxbase=22, + rcutfac=4.604694451, + lambda_value=3.059235105, + lmin=[1, 1, 1, 1, 1, 1], + cutoff=10.0, + ) + self.assertEqual(mat_a.shape, (16, 68)) From f36b0719fa6173d80bae7e45a7e4a7421006f820 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Sun, 14 Apr 2024 17:29:52 -0500 Subject: [PATCH 2/6] Add Fitsnap as depdenceny --- .ci_support/environment-lammps.yml | 3 ++- pyproject.toml | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/.ci_support/environment-lammps.yml b/.ci_support/environment-lammps.yml index 4a201417c..6e3ab212c 100644 --- a/.ci_support/environment-lammps.yml +++ b/.ci_support/environment-lammps.yml @@ -1,4 +1,5 @@ channels: - conda-forge dependencies: -- lammps =2023.11.21 \ No newline at end of file +- lammps =2023.11.21 +- fitsnap3 =3.1.0.2 \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index f030475d5..0aa713815 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -56,6 +56,7 @@ phonopy = [ "spglib==2.4.0", ] pyxtal = ["pyxtal==0.6.2"] +fitsnap3 = ["fitsnap3==3.1.0.2"] [tool.setuptools.packages.find] include = ["structuretoolkit*"] From b2cfe3313735adefdf285070c7439e1fc8f056d7 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Sun, 14 Apr 2024 17:33:56 -0500 Subject: [PATCH 3/6] reuse existing function --- structuretoolkit/analyse/fitsnap.py | 29 ++--------------------------- 1 file changed, 2 insertions(+), 27 deletions(-) diff --git a/structuretoolkit/analyse/fitsnap.py b/structuretoolkit/analyse/fitsnap.py index 87b8ea0f1..8990e7735 100644 --- a/structuretoolkit/analyse/fitsnap.py +++ b/structuretoolkit/analyse/fitsnap.py @@ -3,6 +3,7 @@ import random import numpy as np from fitsnap3lib.fitsnap import FitSnap +from structuretoolkit.analyse.snap import _get_lammps_compatible_cell def get_ace_descriptor_derivatives( @@ -190,32 +191,6 @@ def _ase_scraper(data) -> list: raise Exception("Argument must be list or dictionary for ASE scraper.") -def _get_apre(cell): - """ - Calculate transformed ASE cell for LAMMPS calculations. Thank you Jan Janssen! - - Args: - cell: ASE atoms cell. - - Returns transformed cell as np.array which is suitable for LAMMPS. - """ - a, b, c = cell - an, bn, cn = [np.linalg.norm(v) for v in cell] - - alpha = np.arccos(np.dot(b, c) / (bn * cn)) - beta = np.arccos(np.dot(a, c) / (an * cn)) - gamma = np.arccos(np.dot(a, b) / (an * bn)) - - xhi = an - xyp = np.cos(gamma) * bn - yhi = np.sin(gamma) * bn - xzp = np.cos(beta) * cn - yzp = (bn * cn * np.cos(alpha) - xyp * xzp) / yhi - zhi = np.sqrt(cn**2 - xzp**2 - yzp**2) - - return np.array(((xhi, 0, 0), (xyp, yhi, 0), (xzp, yzp, zhi))) - - def _collate_data(atoms, name: str = None, group_dict: dict = None) -> dict: """ Function to organize fitting data for FitSNAP from ASE atoms objects. @@ -229,7 +204,7 @@ def _collate_data(atoms, name: str = None, group_dict: dict = None) -> dict: """ # Transform ASE cell to be appropriate for LAMMPS. - apre = _get_apre(cell=atoms.cell) + apre = _get_lammps_compatible_cell(cell=atoms.cell) R = np.dot(np.linalg.inv(atoms.cell), apre) positions = np.matmul(atoms.get_positions(), R) cell = apre.T From 96b222bac21ddb168a1563168e8619845b59009a Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Sun, 14 Apr 2024 17:36:53 -0500 Subject: [PATCH 4/6] Drop Python 3.9 support --- .github/workflows/unittests.yml | 5 ----- 1 file changed, 5 deletions(-) diff --git a/.github/workflows/unittests.yml b/.github/workflows/unittests.yml index ac2f3e232..9a212c5ad 100644 --- a/.github/workflows/unittests.yml +++ b/.github/workflows/unittests.yml @@ -39,11 +39,6 @@ jobs: label: linux-64-py-3-10 prefix: /usr/share/miniconda3/envs/my-env - - operating-system: ubuntu-latest - python-version: '3.9' - label: linux-64-py-3-9 - prefix: /usr/share/miniconda3/envs/my-env - steps: - uses: actions/checkout@v4 - name: Merge conda environment From 8618e9dffc1706e7dcaa20f4e3c48d63546255b8 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Tue, 16 Apr 2024 19:34:22 -0500 Subject: [PATCH 5/6] clean up interface --- structuretoolkit/analyse/fitsnap.py | 25 ++++++++++++++++++------- tests/test_fitsnap.py | 20 +++++++++++++------- 2 files changed, 31 insertions(+), 14 deletions(-) diff --git a/structuretoolkit/analyse/fitsnap.py b/structuretoolkit/analyse/fitsnap.py index 8990e7735..d7178948b 100644 --- a/structuretoolkit/analyse/fitsnap.py +++ b/structuretoolkit/analyse/fitsnap.py @@ -9,13 +9,18 @@ def get_ace_descriptor_derivatives( structure: Atoms, atom_types: list[str], - ranks: list[int] = [1, 2, 3, 4, 5, 6], - lmax: list[int] = [1, 2, 2, 2, 1, 1], - nmax: list[int] = [22, 2, 2, 2, 1, 1], + ranks: list[int] = [1, 2, 3, 4], + lmax: list[int] = [0, 5, 2, 1], + nmax: list[int] = [22, 5, 3, 1], + mumax: int = 1, nmaxbase: int = 22, - rcutfac: float = 4.604694451, - lambda_value: float = 3.059235105, - lmin: list[int] = [1, 1, 1, 1, 1, 1], + erefs: list[float] = [0.0], + rcutfac: float = 4.5, + rcinner: float = 1.2, + drcinner: float = 0.01, + RPI_heuristic: str = "root_SO3_span", + lambda_value: float = 1.275, + lmin: list[int] = [0, 0, 1, 1], bzeroflag: bool = True, cutoff: float = 10.0, ) -> np.ndarray: @@ -43,13 +48,19 @@ def get_ace_descriptor_derivatives( "ranks": " ".join([str(r) for r in ranks]), "lmax": " ".join([str(l) for l in lmax]), "nmax": " ".join([str(n) for n in nmax]), + "mumax": mumax, "nmaxbase": nmaxbase, "rcutfac": rcutfac, + "erefs": " ".join([str(e) for e in erefs]), + "rcinner": rcinner, + "drcinner": drcinner, + "RPI_heuristic": RPI_heuristic, "lambda": lambda_value, "type": " ".join(atom_types), "lmin": " ".join([str(l) for l in lmin]), - "bzeroflag": True, + "bzeroflag": bzeroflag, "bikflag": True, + "dgradflag": True, }, "CALCULATOR": { "calculator": "LAMMPSPACE", diff --git a/tests/test_fitsnap.py b/tests/test_fitsnap.py index 812ab2b43..044ac7cd5 100644 --- a/tests/test_fitsnap.py +++ b/tests/test_fitsnap.py @@ -55,13 +55,19 @@ def test_get_ace_descriptor_derivatives(self): mat_a = get_ace_descriptor_derivatives( structure=self.structure, atom_types=self.type, - ranks=[1, 2, 3, 4, 5, 6], - lmax=[1, 2, 2, 2, 1, 1], - nmax=[22, 2, 2, 2, 1, 1], + ranks=[1, 2, 3, 4], + lmax=[0, 5, 2, 1], + nmax=[22, 5, 3, 1], + mumax=1, nmaxbase=22, - rcutfac=4.604694451, - lambda_value=3.059235105, - lmin=[1, 1, 1, 1, 1, 1], + erefs=[0.0], + rcutfac=4.5, + rcinner=1.2, + drcinner=0.01, + RPI_heuristic="root_SO3_span", + lambda_value=1.275, + lmin=[0, 0, 1, 1], + bzeroflag=True, cutoff=10.0, ) - self.assertEqual(mat_a.shape, (16, 68)) + self.assertEqual(mat_a.shape, (16, 141)) From 9ceca67da339e3aac7a04161907c94de402a556a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 18 Jul 2024 06:33:30 +0000 Subject: [PATCH 6/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- structuretoolkit/analyse/__init__.py | 8 ++++++-- structuretoolkit/analyse/fitsnap.py | 6 ++++-- tests/test_fitsnap.py | 6 ++---- 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/structuretoolkit/analyse/__init__.py b/structuretoolkit/analyse/__init__.py index 62fb9c87d..325c74da3 100644 --- a/structuretoolkit/analyse/__init__.py +++ b/structuretoolkit/analyse/__init__.py @@ -31,10 +31,14 @@ try: from structuretoolkit.analyse.fitsnap import ( - get_snap_descriptor_derivatives as get_snap_descriptor_derivatives_fitsnap, - get_ace_descriptor_derivatives as get_ace_descriptor_derivatives_fitsnap, get_ace_descriptor_derivatives, ) + from structuretoolkit.analyse.fitsnap import ( + get_ace_descriptor_derivatives as get_ace_descriptor_derivatives_fitsnap, + ) + from structuretoolkit.analyse.fitsnap import ( + get_snap_descriptor_derivatives as get_snap_descriptor_derivatives_fitsnap, + ) except ImportError: pass diff --git a/structuretoolkit/analyse/fitsnap.py b/structuretoolkit/analyse/fitsnap.py index d7178948b..f9521fb06 100644 --- a/structuretoolkit/analyse/fitsnap.py +++ b/structuretoolkit/analyse/fitsnap.py @@ -1,8 +1,10 @@ -from ase.atoms import Atoms -from typing import Optional, Union import random +from typing import Optional, Union + import numpy as np +from ase.atoms import Atoms from fitsnap3lib.fitsnap import FitSnap + from structuretoolkit.analyse.snap import _get_lammps_compatible_cell diff --git a/tests/test_fitsnap.py b/tests/test_fitsnap.py index 044ac7cd5..fea3a709c 100644 --- a/tests/test_fitsnap.py +++ b/tests/test_fitsnap.py @@ -29,13 +29,11 @@ def setUpClass(cls): cls.bzeroflag = False cls.quadraticflag = False cls.radelem = [4.0] - cls.type = ['Cu'] + cls.type = ["Cu"] cls.wj = [1.0] def test_get_snap_descriptor_derivatives(self): - n_coeff = len(stk.analyse.get_snap_descriptor_names( - twojmax=self.twojmax - )) + n_coeff = len(stk.analyse.get_snap_descriptor_names(twojmax=self.twojmax)) mat_a = get_snap_descriptor_derivatives( structure=self.structure, atom_types=self.type,