How to use the ase.build.molecule function in ase

local_a = desc.create(atoms).toarray()
            re = RematchKernel()
            envkernel = re.compute_envkernel(local_a, local_a)
            glosim = re.rematch(envkernel, gamma=0.01)

            if (0.99 < glosim < 1.01):
                pass
            else:
                is_pass = False

        # Check randomly a few combinations of molecules, just for no errors
        for molname1 in g2.names:
            for molname2 in g2.names:
                if np.random.rand(1, 1) < 0.99:
                    continue
                atoms1 = molecule(molname1)
                atoms2 = molecule(molname2)

                local_a = desc.create(atoms1).toarray()
                local_b = desc.create(atoms2).toarray()
                if len(local_a) == 0 or len(local_b) == 0:
                    continue

                envkernel = re.compute_envkernel(local_a, local_b)
                glosim = re.rematch(envkernel, gamma=0.01)
        self.assertTrue(is_pass)
        return is_pass

def test_parallel_sparse(self):
        """Tests creating sparse output parallelly.
        """
        # Test indices
        samples = [molecule("CO"), molecule("N2O")]
        desc = copy.deepcopy(default_desc_k2)
        desc.species = ["C", "O", "N"]
        desc.sparse = True
        n_features = desc.get_number_of_features()

        # Multiple systems, serial job
        output = desc.create(
            system=samples,
            positions=[[0], [0, 1]],
            n_jobs=1,
        ).toarray()
        assumed = np.empty((3, n_features))
        assumed[0, :] = desc.create(samples[0], [0]).toarray()
        assumed[1, :] = desc.create(samples[1], [0]).toarray()
        assumed[2, :] = desc.create(samples[1], [1]).toarray()
        self.assertTrue(np.allclose(output, assumed))

'''Reactions using StatMech'''
        ideal_gas_param = presets['idealgas']
        self.H2O_sm = StatMech(atoms=molecule('H2O'),
                               symmetrynumber=2,
                               vib_wavenumbers=[3825.434, 3710.2642, 1582.432],
                               potentialenergy=-6.7598,
                               spin=0.,
                               **ideal_gas_param)
        self.H2_sm = StatMech(atoms=molecule('H2'),
                              symmetrynumber=2,
                              vib_wavenumbers=[4306.1793],
                              potentialenergy=-14.2209,
                              spin=0.,
                              **ideal_gas_param)
        self.O2_sm = StatMech(atoms=molecule('O2'),
                              symmetrynumber=2,
                              vib_wavenumbers=[1556.],
                              potentialenergy=-9.862407,
                              spin=1.,
                              **ideal_gas_param)
        # This is an arbitrary transition state for testing
        self.H2O_TS_sm = StatMech(atoms=molecule('H2O'),
                                  symmetrynumber=1.,
                                  vib_wavenumbers=[4000., 3900., 1600.],
                                  potentialenergy=-5.7598,
                                  spin=0.,
                                  **ideal_gas_param)
        self.rxn_sm = rxn.Reaction(
                reactants=[self.H2_sm, self.O2_sm],
                reactants_stoich=[1., 0.5],
                products=[self.H2O_sm],

def test_parallel_sparse(self):
        """Tests creating sparse output parallelly.
        """
        # Test indices
        samples = [molecule("CO"), molecule("N2O")]
        desc = copy.deepcopy(default_desc_k2)
        desc.species = ["C", "O", "N"]
        desc.sparse = True
        n_features = desc.get_number_of_features()

        # Multiple systems, serial job
        output = desc.create(
            system=samples,
            n_jobs=1,
        ).toarray()
        assumed = np.empty((2, n_features))
        assumed[0, :] = desc.create(samples[0]).toarray()
        assumed[1, :] = desc.create(samples[1]).toarray()
        self.assertTrue(np.allclose(output, assumed))

        # Multiple systems, parallel job

def test_minimize_cell(self):
        """Cell minimization.
        """
        system = molecule("H2O")
        system.set_cell([3, 3, 3])

        # Minimize in x-direction
        x_minimized_system = matid.geometry.get_minimized_cell(system, 0, 0.1)
        x_cell = x_minimized_system.get_cell()
        x_expected_cell = np.array([
            [0.1, 0., 0.],
            [0., 3., 0.],
            [0., 0., 3.]
        ])
        self.assertTrue(np.allclose(x_expected_cell, x_cell, atol=0.001, rtol=0))

        # Minimize in y-direction
        y_minimized_system = matid.geometry.get_minimized_cell(system, 1, 0.1)
        y_cell = y_minimized_system.get_cell()
        y_expected_cell = np.array([

def test_glosim_molecules(self):
        is_pass = True

        # Check if the same molecules give global similarity of around 1
        all_atomtypes = [1, 3, 4, 5, 6, 7, 8, 9, 11, 13, 14, 15, 16, 17]
        desc = SOAP(all_atomtypes, 10.0, 2, 0, periodic=False, crossover=True)
        for molname in g2.names:
            atoms = molecule(molname)
            local_a = desc.create(atoms).toarray()
            re = RematchKernel()
            envkernel = re.compute_envkernel(local_a, local_a)
            glosim = re.rematch(envkernel, gamma=0.01)

            if (0.99 < glosim < 1.01):
                pass
            else:
                is_pass = False

        # Check randomly a few combinations of molecules, just for no errors
        for molname1 in g2.names:
            for molname2 in g2.names:
                if np.random.rand(1, 1) < 0.99:
                    continue
                atoms1 = molecule(molname1)

# Undefined elements
        with self.assertRaises(ValueError):
            desc.get_location((2, 1))
        with self.assertRaises(ValueError):
            desc.get_location(("He", "H"))

        # Check that pairwise distances are not supported
        with self.assertRaises(ValueError):
            loc_oo = desc.get_location(("H", "O"))
            loc_oo = desc.get_location(("H", "C"))
            loc_oo = desc.get_location(("C", "H"))

        # Check that slices in the output are correctly empty or filled
        co2 = molecule("CO2")
        h2o = molecule("H2O")
        co2_out = desc.create(co2)
        h2o_out = desc.create(h2o)

        # H-H
        self.assertTrue(co2_out[:, loc_hh].sum() == 0)
        self.assertTrue(h2o_out[:, loc_hh].sum() != 0)

        # C-C
        self.assertTrue(co2_out[:, loc_cc].sum() != 0)
        self.assertTrue(h2o_out[:, loc_cc].sum() == 0)

        # O-O
        self.assertTrue(co2_out[:, loc_oo].sum() != 0)
        self.assertTrue(h2o_out[:, loc_oo].sum() != 0)

from ase.build import molecule
from ase.neb import NEB
from ase.calculators.emt import EMT
from ase.optimize.fire import FIRE as QuasiNewton

# Optimise molecule
initial = molecule('C2H6')
initial.set_calculator(EMT())
relax = QuasiNewton(initial)
relax.run(fmax=0.05)

# Create final state
final = initial.copy()
final.positions[2:5] = initial.positions[[3, 4, 2]]

# Generate blank images
images = [initial]

for i in range(9):
    images.append(initial.copy())

for image in images:
    image.set_calculator(EMT())

"grid": {"min": 0, "max": 5, "n": 100, "sigma": 0.1},
        "weighting": {"function": "exponential", "scale": 0.5, "cutoff": 1e-3},
    },
    k3={
        "geometry": {"function": "angle"},
        "grid": {"min": 0, "max": 180, "n": 100, "sigma": 0.1},
        "weighting": {"function": "exponential", "scale": 0.5, "cutoff": 1e-3},
    },
    periodic=False,
    normalization="l2_each",
)

# Create
from ase.build import molecule

water = molecule("H2O")

# Create MBTR output for the system
mbtr_water = lmbtr.create(water, positions=[0])

print(mbtr_water)
print(mbtr_water.shape)

# Surface sites
# Build a surface and extract different adsorption positions
from ase.build import fcc111, add_adsorbate
slab_pure = fcc111('Al', size=(2, 2, 3), vacuum=10.0)
slab_ads = slab_pure.copy()
add_adsorbate(slab_ads, 'H', 1.5, 'ontop')
ontop_pos = slab_ads.get_positions()[-1]
add_adsorbate(slab_ads, 'H', 1.5, 'bridge')
bridge_pos = slab_ads.get_positions()[-1]

# ## Create Species for Phase Diagram
# We will be considering six CO/Pt(111) configurations. The configurations have CO adsorbed in different sites and different coverages.
# 
# <img width="600" src="images/configurations.png">
# 
# First, we initialize the species as a dictionary to enable easy `Reaction` initialization.

# In[1]:


from ase.build import molecule

from pmutt.statmech import StatMech, presets

species = {
    'CO': StatMech(name='CO', atoms=molecule('CO'),
                   potentialenergy=-14.8021,
                   vib_wavenumbers=[2121.2], symmetrynumber=1, 
                   **presets['idealgas']),
    'Pt': StatMech(name='Pt', potentialenergy=-383.161235,
                   **presets['electronic']),
    'CO(S) 1/16ML fcc': StatMech(
            name='CO(S) 1/16ML fcc', potentialenergy=-399.48282843,
            vib_wavenumbers=[1731.942697, 349.970617, 322.15111,
                             319.114152, 161.45669],
            **presets['harmonic']),
    'CO(S) 1/16ML br': StatMech(
            name='CO(S) 1/16ML br', potentialenergy=-399.464095,
            vib_wavenumbers=[1831.626557, 394.436054, 388.098645,
                             373.063005, 203.887416, 52.987012],
            **presets['harmonic']),
    'CO(S) 1/16ML top': StatMech(