Examples

See the Quickstart guide for a fast example with the ASE EMT Calculator.

Important: This page only contains tips on how to couple certain ASE calculators to AMS. It is provided for information purposes only. The information may become outdated and is not guaranteed to work for you on your system. SCM does not provide support for how to install different ASE calculators into AMS, nor how to set up the input for them.

Warning

When you manually install packages into the AMS Python environment, you may break the SCM-supported ML Potential packages, for example by installing incompatible versions of dependencies. If this happens, it is easiest to remove the AMS Python virtual environment completely and reinstall the ML Potential packages from the package manager.

Overview of external methods/programs

Method/program

AMS Engine

Type

Website

AIMNet2

ML Potential

ML Potential

Github

ALIGNN-FF

ASE

ML potential

Github

ANI (TorchANI)

ML potential

ML potential

Github

CHGNet

ASE

ML potential

Github

CP2K

ASE

DFT

cp2k.org

DeePMD-kit

ASE

ML potential

Github

DFTpy

ASE

orbital-free DFT

Gitlab

EMT

ASE

Force field

GPAW

ASE

DFT

website

M3GNet

ML potential

ML potential

Github

NequIP

ML potential

ML potential

Github

Open Catalyst Project

ASE

ML potential

opencatalystproject.org

Quantum ESPRESSO

Quantum ESPRESSO

DFT

quantum-espresso.org

sGDML

ML potential

ML potential

Github

VASP

External

DFT

vasp.at

xTB (GFN2-xTB)

ASE

semi-empirical

Github

AIMNet2

See the ML Potential documentation.

ALIGNN-FF

Tested with: AMS2023.101, Ubuntu Linux 22.04, July 19 2023

  • SCM → Packages, select CUDA or CPU under ML options

  • Install TorchANI (will install PyTorch)

  • Install DGL into the AMS python environment using instructions from https://www.dgl.ai/pages/start.html. For example (depending on CUDA version):

amspython -m pip install  dgl -f https://data.dgl.ai/wheels/cu117/repo.html
amspython -m pip install  dglgo -f https://data.dgl.ai/wheels-test/repo.html
  • If there’s an error about the pydantic version, then do amspython -m pip uninstall pydantic followed by amspython -m pip install pydantic==1.8.0

  • Install ALIGNN: amspython -m pip install alignn

  • Place the following file in an easily accessible place, for example /home/user/alignn_ams.py:

from alignn.ff.ff import AlignnAtomwiseCalculator, default_path
import numpy as np


class MyCalculator(AlignnAtomwiseCalculator):
    def calculate(self, atoms, properties, system_changes):
        if sum(atoms.get_pbc()) == 0:
            delta = np.max(atoms.get_positions()) - np.min(atoms.get_positions())
            atoms.set_cell([delta + 10, delta + 10, delta + 10])

        return super().calculate(atoms, properties, system_changes)


def get_calculator(path: str = None):
    if path is None:
        path = default_path()

    return MyCalculator(path=path)

Specify Type File and File /path/to/alignn_ams.py in the Engine ASE block:

Listing 1 alignn_ams.run
#!/bin/sh

$AMSBIN/ams -n 1 <<eor
Task GeometryOptimization

System
    Atoms
        O 0. 0. 0.
        H 1. 0. 0.
        H -0.7 0.7 0.
    End
End

Engine ASE
    Type File
    File /path/to/alignn_ams.py
EndEngine
eor

ANI (TorchANI)

See the ML Potential documentation.

CHGNet

Tested with: AMS2023.101, Ubuntu Linux 22.04, July 5 2023

  • SCM → Packages, select CUDA or CPU under ML options

  • Install TorchANI (will install PyTorch) and pymatgen (required by CHGNet)

  • In a terminal: $AMSBIN/amspython -m pip install chgnet==0.1.4

  • Test that the following does not report any error:

$AMSBIN/amspython -c "from chgnet.model.dynamics import CHGNetCalculator"
  • Download chgnet_ams.py and place it in an easily accessible place, for example /home/user/chgnet_ams.py.

Listing 2 chgnet_ams.py
from __future__ import annotations
from chgnet.model.dynamics import CHGNetCalculator
import numpy as np
from ase import Atoms
from pymatgen.io.ase import AseAtomsAdaptor
from ase.calculators.calculator import Calculator, all_changes, all_properties


class AMSCHGNetCalculator(CHGNetCalculator):
    """
    Override the CHGNetCalculator calculate() method.

    Modifications made:

    * add a lattice for non-periodic systems with a gap of at least 10 angstrom

    * only calculate the stress tensor if requested

    The original code can be found at https://github.com/CederGroupHub/chgnet
    """

    def calculate(
        self,
        atoms: Atoms | None = None,
        properties: list | None = None,
        system_changes: list | None = None,
    ) -> None:
        """Calculate various properties of the atoms using CHGNet.

        Args:
            atoms (Atoms | None): The atoms object to calculate properties for.
            properties (list | None): The properties to calculate.
                Default is all properties.
            system_changes (list | None): The changes made to the system.
                Default is all changes.
        """

        if sum(atoms.get_pbc()) == 0:
            delta = np.max(atoms.get_positions()) - np.min(atoms.get_positions())
            atoms.set_cell([delta + 10, delta + 10, delta + 10])

        properties = properties or all_properties
        system_changes = system_changes or all_changes
        Calculator.calculate(
            self,
            atoms=atoms,
            properties=properties,
            system_changes=system_changes,
        )

        # Run CHGNet
        structure = AseAtomsAdaptor.get_structure(atoms)
        graph = self.model.graph_converter(structure)
        task = "e"
        if "forces" in properties:
            task += "f"
        if "stress" in properties:
            task += "s"
        if "magmoms" in properties:
            task += "m"
        model_prediction = self.model.predict_graph(graph.to(self.device), task=task)

        # Convert Result
        factor = 1 if not self.model.is_intensive else structure.composition.num_atoms

        self.results = dict()
        self.results["energy"] = model_prediction["e"] * factor
        self.results["free_energy"] = model_prediction["e"] * factor
        if "f" in model_prediction:
            self.results["forces"] = model_prediction["f"]
        if "s" in model_prediction:
            self.results["stress"] = model_prediction["s"] * self.stress_weight
        if "m" in model_prediction:
            self.results["magmoms"] = model_prediction["m"]


def get_calculator(use_device: str = "cpu", **kwargs):
    return AMSCHGNetCalculator(use_device=use_device, **kwargs)
  • Run AMS with the ASE engine and specify File /path/to/chgnet_ams.py (the path must be absolute, not relative):

Listing 3 chgnet_ams.run
#!/bin/sh

$AMSBIN/ams -n 1 <<eor
Task GeometryOptimization

System
    Atoms
        O 0. 0. 0.
        H 1. 0. 0.
        H -0.7 0.7 0.
    End
End

Engine ASE
    File /path/to/chgnet_ams.py
    Arguments
        use_device = "cpu"
    End
EndEngine
eor

If you have a CUDA-enabled GPU, you can set use_device = "cuda" to run on the GPU instead.

CP2K

Tested with: AMS2023.101, Ubuntu Linux 22.04, July 5 2023

The below works for single-node calculation but fails for multinode (MPI parallelization).

  • Install CP2K. In a terminal: sudo apt install cp2k

  • Run AMS with the CP2K ASE calculator cp2k_ams.run:

Listing 4 cp2k_ams.run
#!/bin/sh

export SCM_DISABLE_MPI=1

# set OMP_NUM_THREADS to the appropriate number below

$AMSBIN/ams -n 1 <<eor

Task GeometryOptimization

System
  Atoms
              O       2.8115000409       2.5498605363       2.0000000000 
              H       2.0000000000       2.0000000000       2.0000000000 
              H       3.6254485609       2.0005857872       2.0000000000 
  End
  Lattice
         5.6254485609     0.0000000000     0.0000000000
         0.0000000000     4.5498605363     0.0000000000
         0.0000000000     0.0000000000     4.0000000000
  End
End

Engine ASE
  Type Import
  Import ase.calculators.cp2k.CP2K
  # see the ASE CP2K documentation for details about the arguments
  Arguments
     command = "OMP_NUM_THREADS=2 cp2k_shell"   # set OMP_NUM_THREADS here
     cutoff = 4000 # eV
     stress_tensor = False   # set stress_tensor here  (defaults to True)
     xc = "PBE"
     inp = """
     &FORCE_EVAL
       &DFT
         &KPOINTS
           SCHEME GAMMA
         &END KPOINTS
         &SCF
           ADDED_MOS 10
           &SMEAR
             METHOD FERMI_DIRAC
             ELECTRONIC_TEMPERATURE [K] 500.0
           &END SMEAR
         &END SCF
       &END DFT
     &END FORCE_EVAL
    """
  End
EndEngine

eor

DeePMD-kit

Tested with: AMS2023.104, Ubuntu Linux 22.04, 14 Nov 2023

  • Install deepmd-kit into the AMS python environment. This will place the db binary in a location like /home/user/.scm/python/AMS2023.1.venv/bin

  • Either train your own model or download one from the AIS Square

  • If you download a model you may need to convert it using db convert-from, see the DeepMD-kit documentation

Then specify Type Import and specify the path to the model (.pb) file in the Arguments:

Listing 5 deepmd-kit_ams.run
#!/bin/sh
NSCM=1 $AMSBIN/ams <<EOF
system
  Atoms
              O      -0.0008161597       0.3663784285      -0.0000000000 
              H      -0.8123162006      -0.1834821079      -0.0000000000 
              H       0.8131323603      -0.1828963206       0.0000000000 
  End
End

Task GeometryOptimization

Engine ASE
  Arguments
     
    model = '/absolute/path/to/graph.pb'

  End
  Import deepmd.calculator.DP
  Type Import
EndEngine
EOF

DFTpy

Tested with: AMS2023.101, Ubuntu Linux 22.04, August 3 2023

$AMSBIN/amspython -m pip install dftpy
$AMSBIN/amspython -m pip install pylibxc2
git clone https://gitlab.com/pavanello-research-group/local-pseudopotentials

Set the environment variable to the path to the pseudopotentials, for example

export DFTPY_DATA_PATH=`readlink -f local-pseudopotentials/BLPS/LDA/reci`

An ASE Calculator with some settings for Al is defined in dftpy_calculator.py:

#!/usr/bin/env amspython
import os
from dftpy.config import DefaultOption, OptionFormat
from dftpy.api.api4ase import DFTpyCalculator
from ase.calculators.calculator import Calculator


class MyCalculator(Calculator):
    implemented_properties = ["energy", "forces", "stress"]

    def __init__(self, config, **kwargs):
        Calculator.__init__(self, **kwargs)
        self.dftpy_calculator = DFTpyCalculator(config=config)

    def calculate(self, atoms, properties=None, system_changes=None):
        super().calculate(atoms, properties, system_changes)
        self.results = dict()
        self.results["energy"] = self.dftpy_calculator.get_potential_energy(atoms)
        self.results["forces"] = self.dftpy_calculator.get_forces(atoms)
        self.results["stress"] = self.dftpy_calculator.get_stress(atoms)


def get_calculator():
    config = DefaultOption()
    config["PATH"]["pppath"] = os.environ.get(
        "DFTPY_DATA_PATH",
        "/home/hellstrom/software/local-pseudopotentials/BLPS/LDA/reci",
    )
    config["PP"]["Al"] = "al.lda.recpot"
    config["OPT"]["method"] = "TN"
    config["KEDF"]["kedf"] = "WT"
    config["JOB"]["calctype"] = "Energy Force"
    config = OptionFormat(config)
    calc = MyCalculator(config=config)

    return calc

This file is referenced inside the Engine ASE block in the input to AMS:

#!/bin/sh
export SCM_DISABLE_MPI=1

$AMSBIN/ams <<EOF
Engine ASE
  File /home/hellstrom/dftpy_calculator.py # change this!
  Type File
EndEngine

MolecularDynamics
  InitialVelocities
    Temperature 1000
    Type Random
  End
  NSteps 20
  Thermostat
    Tau 100.0
    Temperature 1000
    Type NHC
  End
  Barostat
    Type MTK
    Pressure 1.0
    Tau 10000
  End
  TimeStep 1.0
  Trajectory
    SamplingFreq 1
  End
End

Task MolecularDynamics

system
  Atoms
             Al       0.0000000000       0.0000000000       0.0000000000 
             Al       0.0000000000       2.1200000000       2.1200000000 
             Al       2.1200000000       0.0000000000       2.1200000000 
             Al       2.1200000000       2.1200000000       0.0000000000 
             Al       0.0000000000       0.0000000000       4.2400000000 
             Al       0.0000000000       2.1200000000       6.3600000000 
             Al       2.1200000000       0.0000000000       6.3600000000 
             Al       2.1200000000       2.1200000000       4.2400000000 
             Al       0.0000000000       4.2400000000       0.0000000000 
             Al       0.0000000000       6.3600000000       2.1200000000 
             Al       2.1200000000       4.2400000000       2.1200000000 
             Al       2.1200000000       6.3600000000       0.0000000000 
             Al       0.0000000000       4.2400000000       4.2400000000 
             Al       0.0000000000       6.3600000000       6.3600000000 
             Al       2.1200000000       4.2400000000       6.3600000000 
             Al       2.1200000000       6.3600000000       4.2400000000 
             Al       4.2400000000       0.0000000000       0.0000000000 
             Al       4.2400000000       2.1200000000       2.1200000000 
             Al       6.3600000000       0.0000000000       2.1200000000 
             Al       6.3600000000       2.1200000000       0.0000000000 
             Al       4.2400000000       0.0000000000       4.2400000000 
             Al       4.2400000000       2.1200000000       6.3600000000 
             Al       6.3600000000       0.0000000000       6.3600000000 
             Al       6.3600000000       2.1200000000       4.2400000000 
             Al       4.2400000000       4.2400000000       0.0000000000 
             Al       4.2400000000       6.3600000000       2.1200000000 
             Al       6.3600000000       4.2400000000       2.1200000000 
             Al       6.3600000000       6.3600000000       0.0000000000 
             Al       4.2400000000       4.2400000000       4.2400000000 
             Al       4.2400000000       6.3600000000       6.3600000000 
             Al       6.3600000000       4.2400000000       6.3600000000 
             Al       6.3600000000       6.3600000000       4.2400000000 
  End
  Lattice
         8.4800000000     0.0000000000     0.0000000000
         0.0000000000     8.4800000000     0.0000000000
         0.0000000000     0.0000000000     8.4800000000
  End
End

EOF

EMT

See the Quickstart guide.

GPAW

Tested with: AMS2023.102, Ubuntu Linux 22.04, August 1 2023

GPAW is a planewave density functional theory code.

export C_INCLUDE_PATH=$AMSBIN/python3.8/include/python3.8/
amspython -m pip install gpaw
# VENV_BIN is something like /home/user/.scm/python/AMS2023.1.venv/bin
VENV_BIN=$(dirname $(amspython -c "import sys; print(sys.executable)"))
# set TARGET_DIR appropriately
TARGET_DIR=/home/user/gpaw
# Download and install the PAW dataset
amspython $VENV_BIN/gpaw install-data $TARGET_DIR

Follow the instructions from the install-data command.

  • Download GPAW_calculator.py and place it in an easily accessible place, for example /home/user/GPAW_calculator.py.

Listing 6 GPAW_calculator.py
import numpy
import ase
import gpaw


class ASEGPAWCalculator(ase.calculators.calculator.Calculator):
    def __init__(
        self,
        pbc=[True, True, True],
        cancel_total_force=False,
        charge=0,
        name="atoms",
        **gpaw_kwargs,
    ):
        self.name = name
        self.counter = 1
        self.pbc = pbc
        self.cancel = cancel_total_force
        self._kwargs = gpaw_kwargs
        ase.calculators.calculator.Calculator.__init__(self)
        self._setup_gpaw(charge)

    def calculate(self, atoms=None, properties=None, system_changes=None):
        atoms.center()
        atoms.set_pbc(self.pbc)
        super().calculate(atoms, properties, system_changes)

        self.gpaw.calculate(atoms, properties, system_changes)
        self.results = self.gpaw.results
        # remove total force on the molecule
        if self.cancel:
            molecule_force = self.results["forces"].sum(axis=0)
            per_atom_force = molecule_force / self.results["forces"].shape[0]
            self.results["forces"] -= per_atom_force

    def _setup_gpaw(self, charge):
        self.charge = charge
        txt = self.name
        if self.counter > 1:
            txt = txt + f"_{self.counter}"
        txt = txt + ".txt"
        self.gpaw = gpaw.GPAW(txt=txt, **self._kwargs)
        self.counter += 1

    @property
    def implemented_properties(self):
        return self.gpaw.implemented_properties


# AMS looks for "get_calculator"
get_calculator = ASEGPAWCalculator
  • Run AMS with the ASE engine and specify File /path/to/GPAW_calculator.py (the path must be absolute, not relative):

Listing 7 GPAW_ams.run
AMS_JOBNAME=gpaw $AMSBIN/ams -n 1 <<EOF
properties
  gradients yes
End

system
  Atoms
   H      4.630000    5.000000    5.000000
   H      5.370000    5.000000    5.000000
  End
  Lattice
    10.0 0.0 0.0
    0.0 10.0 0.0
    0.0 0.0 10.0
  End
End

task GeometryOptimization

GeometryOptimization
  Method Quasi-Newton
  Convergence
    Gradients 0.00019446905
  End
End

Engine ASE
  File /path/to/GPAW_calculator.py
EndEngine

EOF

GPAW always requires a lattice defined in the AMS system since they are part of the basis set definition for planewaves. For non-periodic systems you can turn off periodic boundary conditions in GPAW by specifying the following block in the ASE Engine:

Arguments
    pbc = [False, False, False]
End

M3GNet

See the ML Potential documentation.

NequIP

See the ML Potential documentation.

Open Catalyst Project

Tested with: AMS2023.102, Ubuntu Linux 22.04, August 1 2023

The Open Catalyst Project (OCP) has several pre-trained machine learning potentials for catalytic systems. They write “The aim is to use AI to model and discover new catalysts for use in renewable energy storage to help in addressing climate change”.

  • To install OCP, first install all ML Potential packages through AMSPackages.

  • Next create a directory for the OCP source, for example /home/user/programs/ocp.

  • Obtain the OCP source code and go into the directory:

    git clone https://github.com/Open-Catalyst-Project/ocp.git /home/user/programs/ocp
    cd /home/user/programs/ocp
    
  • Install ocp and the requirements. This requires suitable compilers to be installed, see the respective documentation of each package:

    amspython -m pip install -e .
    amspython -m pip install torch-geometric
    export CPPFLAGS=-I$AMSBIN/python3.8/include/python3.8/
    amspython -m pip install torch-scatter #takes a while
    amspython -m pip install torch-sparse #also takes a while
    amspython -m pip install mbdb
    amspython -m pip install orjson
    amspython -m pip install wandb
    amspython -m pip install lmdb
    amspython -m pip install numba
    amspython -m pip install e3nn
    
  • Obtain a checkpoint file from the OCP project, which contains a model architecture and pre-trained parameters. For example:

    wget https://dl.fbaipublicfiles.com/opencatalystproject/models/2022_09/oc22/s2ef/gndt_oc22_all_s2ef.pt /home/user/Projects/ocp/gndt_oc22_all_s2ef.pt
    
  • Run AMS with the OCP ASE Calculator and specify the absolute path to the checkpoint file ocp_ams.run:

    Listing 8 ocp_ams.run
    #!/bin/sh
    
    $AMSBIN/ams -n 1 <<eor
    Task GeometryOptimization
    
    System
        Atoms
            O 0. 0. 0.
            H 1. 0. 0.
            H -0.7 0.7 0.
        End
    End
    
    Engine ASE
        Type Import
        Import ocpmodels.common.relaxation.ase_utils.OCPCalculator
        !absolute path
        Arguments
            checkpoint_path = "/home/user/Projects/ocp/gndt_oc22_all_s2ef.pt"
            cpu = True
        End
    EndEngine
    eor
    

Important: OCP models typically only predict energies and forces but no stress tensor. When the stress tensor is required, e.g. for lattice optimizations, the stress tensor is computed using finite differences which substantially slows down calculations.

Quantum ESPRESSO

See the Quantum ESPRESSO documentation.

sGDML

See the ML Potential documentation.

VASP

See the VASP via AMS documentation.

xTB (GFN2-xTB) (method 2)

Tested with: AMS2023.101, Ubuntu Linux 22.04, July 12 2023

Follow the instructions to install xtb-python into a separate anaconda environment. See the xtb-python documentation for details.

Because xTB runs in a separate python environment, the communication with AMS happens through files. Running xTB in this way introduces significant overhead, as the program is restarted for every evaluation. This also means that the SCF runs from scratch every time. The overhead becomes less noticeable for larger molecules.

If you receive a Segmentation fault error message, try to run ulimit -s unlimited. See also the xtb documentation.

  • Download xtb_calculator.py and place it in an easily accessible place, for example /home/user/xtb_calculator.py. IMPORTANT: Modify the path to the conda environment where you have installed xtb-python!

Listing 10 xtb_calculator.py
import os
from ase.calculators.calculator import Calculator, all_changes
import subprocess
import pickle

""" 
Example for using xtb (GFN2-xTB) from a different python environment together with AMS 

IMPORTANT: Set the xtb_python_environment variable to the conda python executable where xtb is installed!

Ideally xtb-python would be installed into the AMS python environment, then one
could just access the XTB calculator directly without having to go the
roundabout way via files. But it's easier to install xtb into a separate conda environment.

To install xtb and xtb-python into a conda environment, do:
conda config --add channels conda-forge
conda install xtb
conda install xtb-python
conda install ase

"""

xtb_python_environment = "/home/hellstrom/anaconda3/bin/python3"


class MyCalculator(Calculator):
    """this class is used inside the AMS python environment"""

    implemented_properties = ["energy", "forces", "charges"]

    def __init__(self, python: str, file: str):
        Calculator.__init__(self)
        self.python = python
        self.file = file

    def calculate(self, atoms=None, properties=["energy"], system_changes=all_changes):
        with open("atoms.pkl", "wb") as f:
            pickle.dump(atoms, f)
        subprocess.run([self.python, self.file])
        with open("results.pkl", "rb") as f:
            self.results = pickle.load(f)
        if os.path.exists("atoms.pkl"):
            os.remove("atoms.pkl")
        if os.path.exists("results.pkl"):
            os.remove("results.pkl")


def get_calculator(**kwargs):
    """this function is imported in the amspython environment"""
    return MyCalculator(python=os.path.abspath(xtb_python_environment), file=os.path.abspath(__file__))


def run_xtb_and_pickle():
    """this is run in the python environment where xTB is installed"""
    from xtb.ase.calculator import XTB

    with open("atoms.pkl", "rb") as f:
        atoms = pickle.load(f)
    atoms.calc = XTB(method="GFN2-xTB")
    atoms.get_forces()
    with open("results.pkl", "wb") as f:
        pickle.dump(atoms.calc.results, f)


if __name__ == "__main__":
    """this is run in the python environment where xTB is installed"""
    run_xtb_and_pickle()
  • Run AMS with the ASE engine and specify File /path/to/xtb_calculator.py (the path must be absolute, not relative):

Listing 11 run_xtb.py
#!/usr/bin/env amspython
import scm.plams as plams
from os.path import abspath

"""
Example showing how to run the external xtb program using GFN2-xTB.

NOTE 1: To instead run GFN1-xTB, you can use the DFTB engine inside AMS.

NOTE 2: You must modify the path below to give the correct path to the xtb_calculator.py file

NOTE 3: To run this script, use the command 
$AMSBIN/amspython run_xtb.py
"""


def main():
    plams.init()
    molecule = plams.from_smiles("O")
    s = plams.Settings()
    s.input.ams.Task = "GeometryOptimization"
    s.input.ASE.File = abspath("/home/hellstrom/xtb_calculator.py")
    s.runscript.nproc = 1  # make sure to run AMS in serial!
    job = plams.AMSJob(settings=s, molecule=molecule, name="opt_GFN2-xTB")
    job.run()

    optimized_molecule = job.results.get_main_molecule()

    plams.log(f"Energy: {job.results.get_energy(unit='eV'):.3f} eV")

    # access atomic charges
    # the charges are stored in "Other%charges" on ase.rkf
    plams.log("Atomic charges:")
    try:
        charges = job.results.get_charges()
        for at, charge in zip(optimized_molecule, charges):
            plams.log(f"{at.symbol} {charge:.3f}")
    except KeyError:
        plams.log("Couldn't get charges")


if __name__ == "__main__":
    main()

The above shows a PLAMS script but you can also use a normal .run file. See the CHGNet example and modify accordingly.