index | modules | gitlab

# Source code for ase.calculators.dmol

    """This module defines an ASE interface to DMol3.

    Contacts
    --------
    Adam Arvidsson <adam.arvidsson@chalmers.se>
    Erik Fransson  <erikfr@chalmers.se>
    Anders Hellman <anders.hellman@chalmers.se>

    DMol3 environment variables
    ----------------------------
    DMOL_COMMAND should point to the RunDmol script and specify the number of cores
    to prallelize over

    export DMOL_COMMAND="./RunDmol.sh -np 16"

    Example
    --------
    >>> from ase.build import bulk
    >>> from ase.calculators import DMol3

    >>> atoms = bulk('Al','fcc')
    >>> calc = DMol3()
    >>> atoms.calc = calc
    >>> print 'Potential energy %5.5f eV' % atoms.get_potential_energy()

    DMol3 calculator functionality
    -------------------------------
    This calculator does support all the functionality in DMol3.

    Firstly this calculator is limited to only handling either fully
    periodic structures (pbc = [1,1,1]) or non periodic structures (pbc=[0,0,0]).

    Internal relaxations are not supported by the calculator,
    only support for energy and forces is implemented.

    Reading eigenvalues and kpts are supported.
    Be careful with kpts and their directions (see internal coordinates below).

    Outputting the full electron density or specific bands to .grd files can be
    achieved with the plot command. The .grd files can be converted to the cube
    format using grd_to_cube().

    DMol3 internal coordinates
    ---------------------------
    DMol3 may change the atomic positions / cell vectors in order to satisfy
    certain criterion ( e.g. molecule symmetry axis along z ). Specifically this
    happens when using Symmetry on/auto. This means the forces read from .grad
    will be in a different coordinates system compared to the atoms object used.
    To solve this the rotation matrix that converts the dmol coordinate system
    to the ase coordinate system is found and applied to the forces.

    For non periodic structures (pbc=[0,0,0]) the rotation matrix can be directly
    parsed from the .rot file.
    For fully periodic structures the rotation matrix is found by reading the
    cell vectors and positions used by dmol and then solving the matrix problem
    DMol_atoms * rot_mat = ase_atoms

    DMol3 files
    ------------
    The supported DMol3 file formats are:

    car    structure file - Angstrom and cellpar description of cell.
    incoor structure file - Bohr and cellvector describption of cell.
                            Note: incoor file not used if car file present.
    outmol outfile from DMol3 - atomic units (Bohr and Hartree)
    grad   outfile for forces from DMol3 - forces in Hartree/Bohr
    grd    outfile for orbitals from DMol3 - cellpar in Angstrom

    """

    import os
    import re

    import numpy as np

    from ase import Atoms
    from ase.calculators.calculator import FileIOCalculator, Parameters, ReadError
    from ase.io import read
    from ase.io.dmol import write_dmol_car, write_dmol_incoor
    from ase.units import Bohr, Hartree

    class DMol3(FileIOCalculator):
        """ DMol3 calculator object. """

        implemented_properties = ['energy', 'forces']
        default_parameters = {'functional': 'pbe',
                              'symmetry': 'on'}
        discard_results_on_any_change = True

        def __init__(self, restart=None,
                     ignore_bad_restart_file=FileIOCalculator._deprecated,
                     label='dmol_calc/tmp', atoms=None,
                     command=None, **kwargs):
            """ Construct DMol3 calculator. """

            if command is None:
                if 'DMOL_COMMAND' in self.cfg:
                    command = self.cfg['DMOL_COMMAND'] + ' PREFIX > PREFIX.out'

            super().__init__(restart, ignore_bad_restart_file,
                             label, atoms, command=command,
                             **kwargs)

            # tracks if DMol transformed coordinate system
            self.internal_transformation = False

        def write_input(self, atoms, properties=None, system_changes=None):

            if not np.all(atoms.pbc) and np.any(atoms.pbc):
                raise RuntimeError('PBC must be all true or all false')

            self.clean()   # Remove files from old run
            self.internal_transformation = False
            self.ase_positions = atoms.positions.copy()
            self.ase_cell = atoms.cell.copy()

            FileIOCalculator.write_input(self, atoms, properties, system_changes)

            if np.all(atoms.pbc):
                write_dmol_incoor(self.label + '.incoor', atoms)
            elif not np.any(atoms.pbc):
                write_dmol_car(self.label + '.car', atoms)

            self.write_input_file()
            self.parameters.write(self.label + '.parameters.ase')

        def write_input_file(self):
            """ Writes the input file. """
            with open(self.label + '.input', 'w') as fd:
                self._write_input_file(fd)

        def _write_input_file(self, fd):
            fd.write('%-32s %s\n' % ('calculate', 'gradient'))

            # if no key about eigs
            fd.write('%-32s %s\n' % ('print', 'eigval_last_it'))

            for key, value in self.parameters.items():
                if isinstance(value, str):
                    fd.write('%-32s %s\n' % (key, value))
                elif isinstance(value, (list, tuple)):
                    for val in value:
                        fd.write('%-32s %s\n' % (key, val))
                else:
                    fd.write('%-32s %r\n' % (key, value))

        def read(self, label):
            FileIOCalculator.read(self, label)
            geometry = self.label + '.car'
            output = self.label + '.outmol'
            force = self.label + '.grad'

            for filename in [force, output, geometry]:
                if not os.path.isfile(filename):
                    raise ReadError

            self.atoms = read(geometry)
            self.parameters = Parameters.read(self.label + 'parameters.ase')
            self.read_results()

        def read_results(self):
            finished, message = self.finished_successfully()
            if not finished:
                raise RuntimeError('DMol3 run failed, see outmol file for'
                                   ' more info\n\n%s' % message)

            self.find_dmol_transformation()
            self.read_energy()
            self.read_forces()

        def finished_successfully(self):
            """ Reads outmol file and checks if job completed or failed.

            Returns
            -------
            finished (bool): True if job completed, False if something went wrong
            message (str): If job failed message contains parsed errors, else empty

            """
            finished = False
            message = ""
            for line in self._outmol_lines():
                if line.rfind('Message: DMol3 job finished successfully') > -1:
                    finished = True
                if line.startswith('Error'):
                    message += line
            return finished, message

        def find_dmol_transformation(self, tol=1e-4):
            """Finds rotation matrix that takes us from DMol internal
            coordinates to ase coordinates.

            For pbc = [False, False, False]  the rotation matrix is parsed from
            the .rot file, if this file doesnt exist no rotation is needed.

            For pbc = [True, True, True] the Dmol internal cell vectors and
            positions are parsed and compared to self.ase_cell self.ase_positions.
            The rotation matrix can then be found by a call to the helper
            function find_transformation(atoms1, atoms2)

            If a rotation matrix is needed then self.internal_transformation is
            set to True and the rotation matrix is stored in self.rotation_matrix

            Parameters
            ----------
            tol (float): tolerance for check if positions and cell are the same
            """

            if np.all(self.atoms.pbc):  # [True, True, True]
                dmol_atoms = self.read_atoms_from_outmol()
                if (np.linalg.norm(self.atoms.positions - dmol_atoms.positions) <
                        tol) and (np.linalg.norm(self.atoms.cell -
                                                 dmol_atoms.cell) < tol):
                    self.internal_transformation = False
                else:
                    R, err = find_transformation(dmol_atoms, self.atoms)
                    if abs(np.linalg.det(R) - 1.0) > tol:
                        raise RuntimeError('Error: transformation matrix does'
                                           ' not have determinant 1.0')
                    if err < tol:
                        self.internal_transformation = True
                        self.rotation_matrix = R
                    else:
                        raise RuntimeError('Error: Could not find dmol'
                                           ' coordinate transformation')
            elif not np.any(self.atoms.pbc):  # [False,False,False]
                try:
                    data = np.loadtxt(self.label + '.rot')
                except OSError:
                    self.internal_transformation = False
                else:
                    self.internal_transformation = True
                    self.rotation_matrix = data[1:].transpose()

        def read_atoms_from_outmol(self):
            """ Reads atomic positions and cell from outmol file and returns atoms
            object.

            If no cell vectors are found in outmol the cell is set to np.eye(3) and
            pbc 000.

            Formatting for cell in outmol :
             translation vector [a0]    1    5.1    0.0    5.1
             translation vector [a0]    2    5.1    5.1    0.0
             translation vector [a0]    3    0.0    5.1    5.1

            Formatting for positions in outmol:
            df              ATOMIC  COORDINATES (au)
            df            x          y          z
            df   Si     0.0   0.0   0.0
            df   Si     1.3   3.5   2.2
            df  binding energy      -0.2309046Ha

            Returns
            -------
            atoms (Atoms object): read atoms object
            """

            lines = self._outmol_lines()
            found_cell = False
            cell = np.zeros((3, 3))
            symbols = []
            positions = []
            pattern_translation_vectors = re.compile(r'\s+translation\s+vector')
            pattern_atomic_coordinates = re.compile(r'df\s+ATOMIC\s+COORDINATES')

            for i, line in enumerate(lines):
                if pattern_translation_vectors.match(line):
                    cell[int(line.split()[3]) - 1, :] = \
                        np.array([float(x) for x in line.split()[-3:]])
                    found_cell = True
                if pattern_atomic_coordinates.match(line):
                    for ind, j in enumerate(range(i + 2, i + 2 + len(self.atoms))):
                        flds = lines[j].split()
                        symbols.append(flds[1])
                        positions.append(flds[2:5])
            atoms = Atoms(symbols=symbols, positions=positions, cell=cell)
            atoms.positions *= Bohr
            atoms.cell *= Bohr

            if found_cell:
                atoms.pbc = [True, True, True]
                atoms.wrap()
            else:
                atoms.pbc = [False, False, False]
            return atoms

        def read_energy(self):
            """ Find and return last occurrence of Ef in outmole file. """
            energy_regex = re.compile(r'^Ef\s+(\S+)Ha')
            found = False
            for line in self._outmol_lines():
                match = energy_regex.match(line)
                if match:
                    energy = float(match.group(1))
                    found = True
            if not found:
                raise RuntimeError('Could not read energy from outmol')
            self.results['energy'] = energy * Hartree

        def read_forces(self):
            """ Read forces from .grad file. Applies self.rotation_matrix if
            self.internal_transformation is True. """
            with open(self.label + '.grad') as fd:
                lines = fd.readlines()

            forces = []
            for i, line in enumerate(lines):
                if line.startswith('$gradients'):
                    for j in range(i + 1, i + 1 + len(self.atoms)):
                        # force = - grad(Epot)
                        forces.append(np.array(
                            [-float(x) for x in lines[j].split()[1:4]]))

            forces = np.array(forces) * Hartree / Bohr
            if self.internal_transformation:
                forces = np.dot(forces, self.rotation_matrix)
            self.results['forces'] = forces

        def get_eigenvalues(self, kpt=0, spin=0):
            return self.read_eigenvalues(kpt, spin, 'eigenvalues')

        def get_occupations(self, kpt=0, spin=0):
            return self.read_eigenvalues(kpt, spin, 'occupations')

        def get_k_point_weights(self):
            return self.read_kpts(mode='k_point_weights')

        def get_bz_k_points(self):
            raise NotImplementedError

        def get_ibz_k_points(self):
            return self.read_kpts(mode='ibz_k_points')

        def get_spin_polarized(self):
            return self.read_spin_polarized()

        def get_fermi_level(self):
            return self.read_fermi()

        def get_energy_contributions(self):
            return self.read_energy_contributions()

        def get_xc_functional(self):
            return self.parameters['functional']

        def read_eigenvalues(self, kpt=0, spin=0, mode='eigenvalues'):
            """Reads eigenvalues from .outmol file.

            This function splits into two situations:
            1. We have no kpts just the raw eigenvalues ( Gamma point )
            2. We have eigenvalues for each k-point

            If calculation is spin_restricted then all eigenvalues
            will be returned no matter what spin parameter is set to.

            If calculation has no kpts then all eigenvalues
            will be returned no matter what kpts parameter is set to.

            Note DMol does usually NOT print all unoccupied eigenvalues.
            Meaning number of eigenvalues for different kpts can vary.
            """

            assert mode in ['eigenvalues', 'occupations']
            lines = self._outmol_lines()
            pattern_kpts = re.compile(r'Eigenvalues for kvector\s+%d' % (kpt + 1))
            for n, line in enumerate(lines):

                # 1. We have no kpts
                if line.split() == ['state', 'eigenvalue', 'occupation']:
                    spin_key = '+'
                    if self.get_spin_polarized():
                        if spin == 1:
                            spin_key = '-'
                    val_index = -2
                    if mode == 'occupations':
                        val_index = -1
                    values = []
                    m = n + 3
                    while True:
                        if lines[m].strip() == '':
                            break
                        flds = lines[m].split()
                        if flds[1] == spin_key:
                            values.append(float(flds[val_index]))
                        m += 1
                    return np.array(values)

                # 2. We have kpts
                if pattern_kpts.match(line):
                    val_index = 3
                    if self.get_spin_polarized():
                        if spin == 1:
                            val_index = 6
                    if mode == 'occupations':
                        val_index += 1
                    values = []
                    m = n + 2
                    while True:
                        if lines[m].strip() == '':
                            break
                        values.append(float(lines[m].split()[val_index]))
                        m += 1
                    return np.array(values)
            return None

        def _outmol_lines(self):
            with open(self.label + '.outmol') as fd:
                return fd.readlines()

        def read_kpts(self, mode='ibz_k_points'):
            """ Returns list of kpts coordinates or kpts weights.  """

            assert mode in ['ibz_k_points', 'k_point_weights']
            lines = self._outmol_ines()

            values = []
            for n, line in enumerate(lines):
                if line.startswith('Eigenvalues for kvector'):
                    if mode == 'ibz_k_points':
                        values.append([float(k_i)
                                       for k_i in lines[n].split()[4:7]])
                    if mode == 'k_point_weights':
                        values.append(float(lines[n].split()[8]))
            if values == []:
                return None
            return values

        def read_spin_polarized(self):
            """Reads, from outmol file, if calculation is spin polarized."""

            lines = self._outmol_lines()
            for n, line in enumerate(lines):
                if line.rfind('Calculation is Spin_restricted') > -1:
                    return False
                if line.rfind('Calculation is Spin_unrestricted') > -1:
                    return True
            raise OSError('Could not read spin restriction from outmol')

        def read_fermi(self):
            """Reads the Fermi level.

            Example line in outmol:
            Fermi Energy:           -0.225556 Ha     -6.138 eV   xyz text
            """
            lines = self._outmol_lines()
            pattern_fermi = re.compile(r'Fermi Energy:\s+(\S+)\s+Ha')
            for line in lines:
                m = pattern_fermi.match(line)
                if m:
                    return float(m.group(1)) * Hartree
            return None

        def read_energy_contributions(self):
            """Reads the different energy contributions."""

            lines = self._outmol_lines()
            energies = {}
            for n, line in enumerate(lines):
                if line.startswith('Energy components'):
                    m = n + 1
                    while lines[m].strip() != '':
                        energies[lines[m].split('=')[0].strip()] = \
                            float(re.findall(
                                r"[-+]?\d*\.\d+|\d+", lines[m])[0]) * Hartree
                        m += 1
            return energies

        def clean(self):
            """ Cleanup after dmol calculation

            Only removes dmol files in self.directory,
            does not remove the directory itself
            """
            file_extensions = ['basis', 'car', 'err', 'grad', 'input', 'inatm',
                               'incoor', 'kpoints', 'monitor', 'occup', 'outmol',
                               'outatom', 'rot', 'sdf', 'sym', 'tpotl', 'tpdensk',
                               'torder', 'out', 'parameters.ase']
            files_to_clean = ['DMol3.log', 'stdouterr.txt', 'mpd.hosts']

            files = [os.path.join(self.directory, f) for f in files_to_clean]
            files += [''.join((self.label, '.', ext)) for ext in file_extensions]

            for f in files:
                try:
                    os.remove(f)
                except OSError:
                    pass

    # Helper functions
    # ------------------

    def find_transformation(atoms1, atoms2, verbose=False, only_cell=False):
        """ Solves Ax = B where A and B are cell and positions from atoms objects.

        Uses numpys least square solver to solve the problem Ax = B where A and
        B are cell vectors and positions for atoms1 and atoms2 respectively.

        Parameters
        ----------
        atoms1 (Atoms object): First atoms object (A)
        atoms2 (Atoms object): Second atoms object (B)
        verbose (bool): If True prints for each i A[i], B[i], Ax[i]
        only_cell (bool): If True only cell in used, otherwise cell and positions.

        Returns
        -------
        x (np.array((3,3))): Least square solution to Ax = B
        error (float): The error calculated as np.linalg.norm(Ax-b)

        """

        if only_cell:
            N = 3
        elif len(atoms1) != len(atoms2):
            raise RuntimeError('Atoms object must be of same length')
        else:
            N = len(atoms1) + 3

        # Setup matrices A and B
        A = np.zeros((N, 3))
        B = np.zeros((N, 3))
        A[0:3, :] = atoms1.cell
        B[0:3, :] = atoms2.cell
        if not only_cell:
            A[3:, :] = atoms1.positions
            B[3:, :] = atoms2.positions

        # Solve least square problem Ax = B
        lstsq_fit = np.linalg.lstsq(A, B, rcond=-1)
        x = lstsq_fit[0]
        error = np.linalg.norm(np.dot(A, x) - B)

        # Print comparison between A, B and Ax
        if verbose:
            print('%17s %33s %35s %24s' % ('A', 'B', 'Ax', '|Ax-b|'))
            for a, b in zip(A, B):
                ax = np.dot(a, x)
                loss = np.linalg.norm(ax - b)
                print('(', end='')
                for a_i in a:
                    print('%8.5f' % a_i, end='')
                print(')   (', end='')
                for b_i in b:
                    print('%8.5f ' % b_i, end='')
                print(')   (', end='')
                for ax_i in ax:
                    print('%8.5f ' % ax_i, end='')
                print(')   %8.5f' % loss)

        return x, error

    def grd_to_file(atoms, grd_file, new_file):
        """ Reads grd_file and converts data to cube format and writes to
        cube_file.

        Note: content of grd_file and atoms object are assumed to match with the
              same orientation.

        Parameters
        -----------
        atoms (Atoms object): atoms object grd_file data is for
        grd_file (str): filename of .grd file
        new_file (str): filename to write grd-data to, must be ASE format
                        that supports data argument
        """
        from ase.io import write

        atoms_copy = atoms.copy()
        data, cell, origin = read_grd(grd_file)
        atoms_copy.cell = cell
        atoms_copy.positions += origin
        write(new_file, atoms_copy, data=data)

    def read_grd(filename):
        """ Reads .grd file

        Notes
        -----
        origin_xyz is offset with half a grid point in all directions to be
            compatible with the cube format
        Periodic systems is not guaranteed to be oriented correctly
        """
        from ase.geometry.cell import cellpar_to_cell

        with open(filename) as fd:
            lines = fd.readlines()

        cell_data = np.array([float(fld) for fld in lines[2].split()])
        cell = cellpar_to_cell(cell_data)
        grid = [int(fld) + 1 for fld in lines[3].split()]
        data = np.empty(grid)

        origin_data = [int(fld) for fld in lines[4].split()[1:]]
        origin_xyz = cell[0] * (-float(origin_data[0]) - 0.5) / (grid[0] - 1) + \
            cell[1] * (-float(origin_data[2]) - 0.5) / (grid[1] - 1) + \
            cell[2] * (-float(origin_data[4]) - 0.5) / (grid[2] - 1)

        # Fastest index describes which index ( x or y ) varies fastest
        # 1: x  , 3: y
        fastest_index = int(lines[4].split()[0])
        assert fastest_index in [1, 3]
        if fastest_index == 3:
            grid[0], grid[1] = grid[1], grid[0]

        dummy_counter = 5
        for i in range(grid[2]):
            for j in range(grid[1]):
                for k in range(grid[0]):  # Fastest index
                    if fastest_index == 1:
                        data[k, j, i] = float(lines[dummy_counter])
                    elif fastest_index == 3:
                        data[j, k, i] = float(lines[dummy_counter])
                    dummy_counter += 1

        return data, cell, origin_xyz

    if __name__ == '__main__':
        from ase.build import molecule

        atoms = molecule('H2')
        calc = DMol3()
        atoms.calc = calc
        # ~ 60 sec calculation
        print('Potential energy %5.5f eV' % atoms.get_potential_energy())