How to use the ase.units function in ase

def test_apply_strain(self):
            calc = TersoffScr(**Tersoff_PRB_39_5566_Si_C__Scr)
            timestep = 1.0*units.fs

            atoms = ase.io.read('cryst_rot_mod.xyz')
            atoms.set_calculator(calc)

            # constraints
            top = atoms.positions[:, 1].max()
            bottom = atoms.positions[:, 1].min()
            fixed_mask = ((abs(atoms.positions[:, 1] - top) < 1.0) |
                          (abs(atoms.positions[:, 1] - bottom) < 1.0))
            fix_atoms = FixAtoms(mask=fixed_mask)

            # strain
            orig_height = (atoms.positions[:, 1].max() - atoms.positions[:, 1].min())
            delta_strain = timestep*1e-5*(1/units.fs)
            rigid_constraints = False
            strain_atoms = ConstantStrainRate(orig_height, delta_strain)

):
        Calculator.__init__(self, **kwargs)

        self.model = model

        self.atoms_converter = AtomsConverter(
            environment_provider=environment_provider,
            collect_triples=collect_triples,
            device=device,
        )

        self.model_energy = energy
        self.model_forces = forces

        # Convert to ASE internal units
        self.energy_units = MDUnits.parse_mdunit(energy_units) * units.Ha
        self.forces_units = MDUnits.parse_mdunit(forces_units) * units.Ha / units.Bohr

Pourbaix M (1966)
        Atlas of electrochemical equilibria in aqueous solutions.
        No. v. 1 in Atlas of Electrochemical Equilibria in Aqueous Solutions.
        Pergamon Press, New York.

    Returns list of (name, energy) tuples.
    """

    if isinstance(symbols, str):
        symbols = Formula(symbols).count().keys()
    if len(_solvated) == 0:
        for line in _aqueous.splitlines():
            energy, formula = line.split(',')
            name = formula + '(aq)'
            count, charge, aq = parse_formula(name)
            energy = float(energy) * 0.001 * units.kcal / units.mol
            _solvated.append((name, count, charge, aq, energy))
    references = []
    for name, count, charge, aq, energy in _solvated:
        for symbol in count:
            if symbol not in 'HO' and symbol not in symbols:
                break
        else:
            references.append((name, energy))
    return references

def get_zero_point_energy(self, freq=None):
        if freq is None:
            return 0.5 * self.hnu.real.sum()
        else:
            s = 0.01 * units._e / units._c / units._hplanck
            return 0.5 * freq.real.sum() / s

spring constant (eV/A^2)
        m : float
            mass (grams/mole or AMU)
        T : float
            temperature (K)
        method : str
            method for free energy computation, classical or QM.

        Returns
        --------
        float
            free energy of the Einstein crystal (eV/atom)
        """
        assert method in ['classical', 'QM']

        hbar = units._hbar * units.J  # eV/s
        m = m / units.kg              # mass kg
        k = k * units.m**2 / units.J  # spring constant J/m2
        omega = np.sqrt(k / m)        # angular frequency 1/s

        if method == 'classical':
            F_einstein = 3 * units.kB * T * np.log(hbar * omega / (units.kB * T))
        elif method == 'QM':
            log_factor = np.log(1.0 - np.exp(-hbar * omega / (units.kB * T)))
            F_einstein = 3 * units.kB * T * log_factor + 1.5 * hbar * omega

        return F_einstein

import numpy as np
from ase.calculators.calculator import Calculator
from ase import units

k_c = units.Hartree * units.Bohr


class AtomicCounterIon(Calculator):
    implemented_properties = ['energy', 'forces']

    def __init__(self, charge, epsilon, sigma, sites_per_mol=1,
                 rc=7.0, width=1.0):
        """ Counter Ion Calculator.

        A very simple, nonbonded (Coulumb and LJ)
        interaction calculator meant for single atom ions
        to charge neutralize systems (and nothing else)...
        """
        self.rc = rc
        self.width = width
        self.sites_per_mol = sites_per_mol

H[r] /= 2 * self.delta
                r += 1
        H += H.copy().T
        self.H = H
        m = self.atoms.get_masses()
        if 0 in [m[index] for index in self.indices]:
            raise RuntimeError('Zero mass encountered in one or more of '
                               'the vibrated atoms. Use Atoms.set_masses()'
                               ' to set all masses to non-zero values.')

        self.im = np.repeat(m[self.indices]**-0.5, 3)
        omega2, modes = np.linalg.eigh(self.im[:, None] * H * self.im)
        self.modes = modes.T.copy()

        # Conversion factor:
        s = units._hbar * 1e10 / sqrt(units._e * units._amu)
        self.hnu = s * omega2.astype(complex)**0.5

def read(self, method='standard', direction='central', inter = True):
        self.method = method.lower()
        self.direction = direction.lower()
        assert self.method in ['standard', 'frederiksen']
        if direction != 'central':
            raise NotImplementedError(
                'Only central difference is implemented at the moment.')

        # Get "static" dipole moment polarizability and forces
        name = '%s.eq.pckl' % self.name
        [forces_zero, dipole_zero, freq_zero,
            noninPol_zero, pol_zero] = pickle.load(open(name, "rb"))
        self.dipole_zero = (sum(dipole_zero**2)**0.5) / units.Debye
        self.force_zero = max([sum((forces_zero[j])**2)**0.5
                               for j in self.indices])
        self.noninPol_zero = noninPol_zero * (units.Bohr)**3  # Ang**3
        self.pol_zero = pol_zero * (units.Bohr)**3  # Ang**3

        ndof = 3 * len(self.indices)
        H = np.empty((ndof, ndof))
        dpdx = np.empty((ndof, 3))
        dadx = np.empty((ndof, self.pol_zero.shape[0], 3, 3), dtype=complex)

        r = 0
        for a in self.indices:
            for i in 'xyz':
                name = '%s.%d%s' % (self.name, a, i)
                [fminus, dminus, frminus, noninpminus, pminus] = pickle.load(
                    open(name + '-.pckl', "rb"))

def __init__(self, atomsx, express=np.zeros((3,3)), fixstrain=np.ones((3,3))):
        """box relaxation
        atomsx: an Atoms object
        express: external pressure, a 3*3 lower triangular matrix in the unit of GPa;
                 define positive values as compression
        fixstrain: 3*3 matrix as express. 0 fixes strain at the corresponding direction
        """
        self.atomsx = atomsx 
        self.express= express * units.GPa
        if express[0][1]**2+express[0][2]**2+express[1][2]**2 > 1e-5:
           express[0][1] = 0
           express[0][2] = 0
           express[1][2] = 0
           print("warning: xy, xz, yz components of the external pressure will be set to zero")
        self.fixstrain = fixstrain

        cell       = atomsx.get_cell()
        vol        = atomsx.get_volume()
        self.natom = atomsx.get_number_of_atoms()
        avglen     = (vol/self.natom)**(1.0/3.0)
        self.jacobian = avglen * self.natom**0.5
        Atoms.__init__(self,atomsx)

import logging
import os

import numpy as np

from ase import Atoms, units
from tqdm import tqdm

import schnetpack as spk
from schnetpack import Properties

logging.basicConfig(level=os.environ.get("LOGLEVEL", "INFO"))

# Conversion from ppm to atomic units. Alpha is the fine structure constant and 1e6 are the ppm
ppm2au = 1.0 / (units.alpha ** 2 * 1e6)


class OrcaParserException(Exception):
    pass


class OrcaParser:
    main_properties = [
        Properties.energy,
        Properties.forces,
        Properties.dipole_moment,
        Properties.polarizability,
        Properties.shielding,
    ]
    hessian_properties = [
        Properties.hessian,