Test matgl with LAMMPS

.ci_support/environment-matgl.yml

@@ -1,4 +1,8 @@
 channels:
 - conda-forge
 dependencies:
+- lammps =2023.11.21
 - matgl =0.9.2
+- pandas =2.2.0
+- pylammpsmpi =0.2.11
+- jinja2 =3.1.3

.github/workflows/unittests_matgl.yml

@@ -33,4 +33,5 @@ jobs:
         python -m unittest tests/test_evcurve_ase_matgl.py
         python -m unittest tests/test_phonons_ase_matgl.py
         python -m unittest tests/test_quasiharmonic_ase_matgl.py
-        python -m unittest tests/test_ase_md_matgl.py
+        python -m unittest tests/test_ase_md_matgl.py
+        python -m unittest tests/test_evcurve_lammps_matgl.py

tests/static/lammps/potential_LAMMPS/M3GNET/matgl_driver.py

@@ -0,0 +1,127 @@
+"""
+Copyright (c) 2023, AdvanceSoft Corp.
+
+This source code is licensed under the GNU General Public License Version 2
+found in the LICENSE file in the root directory of this source tree.
+"""
+
+from ase import Atoms
+from ase.calculators.mixing import SumCalculator
+
+import matgl
+from matgl.ext.ase import M3GNetCalculator
+
+def m3gnet_initialize(model_name = None, dftd3 = False):
+    """
+    Initialize GNNP of M3GNet.
+    Args:
+        model_name (str): name of model for GNNP.
+        dftd3 (bool): to add correction of DFT-D3.
+    Returns:
+        cutoff: cutoff radius.
+    """
+
+    # Create M3GNetCalculator, that is pre-trained
+    global myCalculator
+
+    if model_name is not None:
+        myPotential = matgl.load_model(model_name)
+    else:
+        myPotential = matgl.load_model("M3GNet-MP-2021.2.8-PES")
+
+    myCalculator = M3GNetCalculator(
+        potential      = myPotential,
+        compute_stress = True,
+        stress_weight  = 1.0
+    )
+
+    # Add DFT-D3 to calculator without three-body term
+    global m3gnetCalculator
+    global dftd3Calculator
+
+    m3gnetCalculator = myCalculator
+    dftd3Calculator  = None
+
+    if dftd3:
+        from dftd3.ase import DFTD3
+        #from torch_dftd.torch_dftd3_calculator import TorchDFTD3Calculator
+
+        dftd3Calculator = DFTD3(
+            method  = "PBE",
+            damping = "d3zero",
+            s9      = 0.0
+        )
+        #dftd3Calculator = TorchDFTD3Calculator(
+        #    xc      = "pbe",
+        #    damping = "zero",
+        #    abc     = False
+        #)
+
+        myCalculator = SumCalculator([m3gnetCalculator, dftd3Calculator])
+
+    # Atoms object of ASE, that is empty here
+    global myAtoms
+
+    myAtoms = None
+
+    return myPotential.model.cutoff
+
+def m3gnet_get_energy_forces_stress(cell, atomic_numbers, positions):
+    """
+    Predict total energy, atomic forces and stress w/ pre-trained GNNP of M3GNet.
+    Args:
+        cell: lattice vectors in angstroms.
+        atomic_numbers: atomic numbers for all atoms.
+        positions: xyz coordinates for all atoms in angstroms.
+    Returns:
+        energy:  total energy.
+        forcces: atomic forces.
+        stress:  stress tensor (Voigt order).
+    """
+
+    # Initialize Atoms
+    global myAtoms
+    global myCalculator
+
+    if myAtoms is not None and len(myAtoms.numbers) != len(atomic_numbers):
+        myAtoms = None
+
+    if myAtoms is None:
+        myAtoms = Atoms(
+            numbers   = atomic_numbers,
+            positions = positions,
+            cell      = cell,
+            pbc       = [True, True, True]
+        )
+
+        myAtoms.calc = myCalculator
+
+    else:
+        myAtoms.set_cell(cell)
+        myAtoms.set_atomic_numbers(atomic_numbers)
+        myAtoms.set_positions(positions)
+
+    # Predicting energy, forces and stress
+    energy = myAtoms.get_potential_energy()
+    forces = myAtoms.get_forces().tolist()
+
+    global m3gnetCalculator
+    global dftd3Calculator
+
+    if dftd3Calculator is None:
+        stress = myAtoms.get_stress().tolist()
+    else:
+        # to avoid the bug of SumCalculator
+        myAtoms.calc = m3gnetCalculator
+        stress1 = myAtoms.get_stress()
+
+        myAtoms.calc = dftd3Calculator
+        stress2 = myAtoms.get_stress()
+
+        stress = stress1 + stress2
+        stress = stress.tolist()
+
+        myAtoms.calc = myCalculator
+
+    return energy, forces, stress
+

tests/test_evcurve_lammps_matgl.py

@@ -0,0 +1,67 @@
+import os
+
+from ase.build import bulk
+import numpy as np
+import unittest
+
+from atomistics.workflows import EnergyVolumeCurveWorkflow, optimize_positions_and_volume
+
+
+try:
+    import pandas
+    import matgl
+    from matgl.ext.ase import M3GNetCalculator
+    from atomistics.calculators import (
+        evaluate_with_lammps, get_potential_by_name
+    )
+
+    skip_lammps_matgl_test = False
+except ImportError:
+    skip_lammps_matgl_test = True
+
+
+@unittest.skipIf(
+    skip_lammps_matgl_test, "LAMMPS and or matgl are not installed, so the corresponding tests are skipped."
+)
+class TestEvCurve(unittest.TestCase):
+    def test_calc_evcurve(self):
+        structure = bulk("Al", cubic=True)
+        df_pot_selected = pandas.DataFrame({
+            "Config": [[
+                "pair_style m3gnet " + os.path.abspath(os.path.join(__file__, "..", "static", "lammps", "potential_LAMMPS", "M3GNET")),
+                "pair_coeff * * M3GNet-MP-2021.2.8-PES Al"
+            ]],
+            "Filename": [[]],
+            "Species": [["Al"]]
+        })
+        task_dict = optimize_positions_and_volume(structure=structure)
+        result_dict = evaluate_with_lammps(
+            task_dict=task_dict,
+            potential_dataframe=df_pot_selected,
+        )
+        workflow = EnergyVolumeCurveWorkflow(
+            structure=result_dict["structure_with_optimized_positions_and_volume"],
+            num_points=11,
+            fit_type='polynomial',
+            fit_order=3,
+            vol_range=0.05,
+            axes=('x', 'y', 'z'),
+            strains=None,
+        )
+        task_dict = workflow.generate_structures()
+        result_dict = evaluate_with_lammps(
+            task_dict=task_dict,
+            potential_dataframe=df_pot_selected,
+        )
+        fit_dict = workflow.analyse_structures(output_dict=result_dict)
+        thermal_properties_dict = workflow.get_thermal_properties(
+            temperatures=[100, 1000],
+            output_keys=["temperatures", "volumes"]
+        )
+        temperatures_ev, volumes_ev = thermal_properties_dict["temperatures"], thermal_properties_dict["volumes"]
+        self.assertTrue(np.isclose(fit_dict['volume_eq'], 66.56048874824006, atol=1e-04))
+        self.assertTrue(np.isclose(fit_dict['bulkmodul_eq'], 50.96266448851179, atol=1e-02))
+        self.assertTrue(np.isclose(fit_dict['b_prime_eq'], 4.674534962000779, atol=1e-02))
+        self.assertEqual(len(temperatures_ev), 2)
+        self.assertEqual(len(volumes_ev), 2)
+        self.assertTrue(volumes_ev[0] < volumes_ev[-1])