WIP: ASE calculator interface

examples/ase_api_example.py

@@ -0,0 +1,89 @@
+"""
+API example modified for an ASE calculator.
+
+The result is meaningless, apart from showing that the code runs, because of using a Morse potential.
+
+"""
+import ase.io
+import numpy as np
+from ase.calculators.morse import MorsePotential
+
+from pygsm.coordinate_systems import DelocalizedInternalCoordinates
+from pygsm.level_of_theories.ase import ASELoT
+from pygsm.optimizers import eigenvector_follow
+from pygsm.potential_energy_surfaces import PES
+from pygsm.utilities import elements, manage_xyz, nifty
+from pygsm.wrappers import Molecule
+
+
+def main(geom):
+    nifty.printcool(" Building the LOT")
+    lot = ASELoT.from_options(MorsePotential(), geom=geom)
+
+    nifty.printcool(" Building the PES")
+    pes = PES.from_options(
+        lot=lot,
+        ad_idx=0,
+        multiplicity=1,
+    )
+
+    nifty.printcool("Building the topology")
+    atom_symbols = manage_xyz.get_atoms(geom)
+    ELEMENT_TABLE = elements.ElementData()
+    atoms = [ELEMENT_TABLE.from_symbol(atom) for atom in atom_symbols]
+    # top = Topology.build_topology(
+    #     xyz,
+    #     atoms,
+    # )
+
+    # nifty.printcool("Building Primitive Internal Coordinates")
+    # p1 = PrimitiveInternalCoordinates.from_options(
+    #     xyz=xyz,
+    #     atoms=atoms,
+    #     addtr=False,  # Add TRIC
+    #     topology=top,
+    # )
+
+    nifty.printcool("Building Delocalized Internal Coordinates")
+    coord_obj1 = DelocalizedInternalCoordinates.from_options(
+        xyz=xyz,
+        atoms=atoms,
+        addtr=False,  # Add TRIC
+    )
+
+    nifty.printcool("Building Molecule")
+    reactant = Molecule.from_options(
+        geom=geom,
+        PES=pes,
+        coord_obj=coord_obj1,
+        Form_Hessian=True,
+    )
+
+    nifty.printcool("Creating optimizer")
+    optimizer = eigenvector_follow.from_options(Linesearch='backtrack', OPTTHRESH=0.0005, DMAX=0.5, abs_max_step=0.5,
+                                                conv_Ediff=0.5)
+
+    nifty.printcool("initial energy is {:5.4f} kcal/mol".format(reactant.energy))
+    geoms, energies = optimizer.optimize(
+        molecule=reactant,
+        refE=reactant.energy,
+        opt_steps=500,
+        verbose=True,
+    )
+
+    nifty.printcool("Final energy is {:5.4f}".format(reactant.energy))
+    manage_xyz.write_xyz('minimized.xyz', geoms[-1], energies[-1], scale=1.)
+
+
+if __name__ == '__main__':
+    diels_adler = ase.io.read("diels_alder.xyz", ":")
+    xyz = diels_adler[0].positions
+
+    # this is a hack
+    geom = np.column_stack([diels_adler[0].symbols, xyz]).tolist()
+    for i in range(len(geom)):
+        for j in [1, 2, 3]:
+            geom[i][j] = float(geom[i][j])
+    # --------------------------
+
+    main(geom)

pygsm/level_of_theories/ase.py

@@ -0,0 +1,103 @@
+"""
+Level Of Theory for ASE calculators
+https://gitlab.com/ase/ase
+
+Written by Tamas K. Stenczel in 2021
+"""
+try:
+    from ase import Atoms
+    from ase.calculators.calculator import Calculator
+    from ase.data import atomic_numbers
+    from ase import units
+except ModuleNotFoundError:
+    print("ASE not installed, ASE-based calculators will not work")
+
+from .base_lot import Lot
+
+
+class ASELoT(Lot):
+    """
+    Warning:
+        multiplicity is not implemented, the calculator ignores it
+    """
+
+    def __init__(self, calculator: Calculator, options):
+        super(ASELoT, self).__init__(options)
+
+        self.ase_calculator = calculator
+
+    @classmethod
+    def from_options(cls, calculator: Calculator, **kwargs):
+        """ Returns an instance of this class with default options updated from values in kwargs"""
+        return cls(calculator, cls.default_options().set_values(kwargs))
+
+    @classmethod
+    def copy(cls, lot, options, copy_wavefunction=True):
+        assert isinstance(lot, ASELoT)
+        return cls(lot.ase_calculator, lot.options.copy().set_values(options))
+
+    def run(self, geom, mult, ad_idx, runtype='gradient'):
+        # run ASE
+        self.run_ase_atoms(xyz_to_ase(geom), mult, ad_idx, runtype)
+
+    def run_ase_atoms(self, atoms: Atoms, mult, ad_idx, runtype='gradient'):
+        # set the calculator
+        atoms.set_calculator(self.ase_calculator)
+
+        # perform gradient calculation if needed
+        if runtype == "gradient":
+            self._Gradients[(mult, ad_idx)] = self.Gradient(- atoms.get_forces() / units.Ha * units.Bohr,
+                                                            'Hartree/Bohr')
+        elif runtype == "energy":
+            pass
+        else:
+            raise NotImplementedError(f"Run type {runtype} is not implemented in the ASE calculator interface")
+
+        # energy is always calculated -> cached if force calculation was done
+        self._Energies[(mult, ad_idx)] = self.Energy(atoms.get_potential_energy() / units.Ha, 'Hartree')
+
+        # write E to scratch
+        self.write_E_to_file()
+
+        self.hasRanForCurrentCoords = True
+
+
+def xyz_to_ase(xyz):
+    """
+
+    Parameters
+    ----------
+    xyz : np.ndarray, shape=(N, 4)
+
+
+    Returns
+    -------
+    atoms : ase.Atoms
+        ASE's Atoms object
+
+    """
+
+    # compatible with list-of-list as well
+    numbers = [atomic_numbers[x[0]] for x in xyz]
+    pos = [x[1:4] for x in xyz]
+    return geom_to_ase(numbers, pos)
+
+
+def geom_to_ase(numbers, positions, **kwargs):
+    """Geometry to ASE atoms object
+
+    Parameters
+    ----------
+    numbers : array_like, shape=(N_atoms,)
+        atomic numbers
+    positions : array_like, shape=(N_atoms,3)
+        positions of atoms in Angstroms
+    kwargs
+
+    Returns
+    -------
+    atoms : ase.Atoms
+        ASE's Atoms object
+    """
+
+    return Atoms(numbers=numbers, positions=positions, **kwargs)