index | modules | gitlab

# Source code for ase.calculators.demon.demon

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

    http://www.demon-software.com

    """
    import os
    import os.path as op
    import subprocess
    import shutil

    import numpy as np

    from ase.units import Bohr, Hartree
    import ase.data
    from ase.calculators.calculator import FileIOCalculator, ReadError
    from ase.calculators.calculator import Parameters, all_changes
    from ase.calculators.calculator import equal
    import ase.io
    from .demon_io import parse_xray

    m_e_to_amu = 1822.88839

    class Parameters_deMon(Parameters):
        """Parameters class for the calculator.
        Documented in Base_deMon.__init__

        The options here are the most important ones that the user needs to be
        aware of. Further options accepted by deMon can be set in the dictionary
        input_arguments.

        """
        def __init__(
                self,
                label='rundir',
                atoms=None,
                command=None,
                restart=None,
                basis_path=None,
                ignore_bad_restart_file=FileIOCalculator._deprecated,
                deMon_restart_path='.',
                title='deMon input file',
                scftype='RKS',
                forces=False,
                dipole=False,
                xc='VWN',
                guess='TB',
                print_out='MOE',
                basis={},
                ecps={},
                mcps={},
                auxis={},
                augment={},
                input_arguments=None):
            kwargs = locals()
            kwargs.pop('self')
            Parameters.__init__(self, **kwargs)

    class Demon(FileIOCalculator):
        """Calculator interface to the deMon code. """

        implemented_properties = [
            'energy',
            'forces',
            'dipole',
            'eigenvalues']

        def __init__(self, **kwargs):
            """ASE interface to the deMon code.

            The deMon2k code can be obtained from http://www.demon-software.com

            The DEMON_COMMAND environment variable must be set to run the executable, in bash it would be set along the lines of
            export DEMON_COMMAND="deMon.4.3.6.std > deMon_ase.out 2>&1"

            Parameters:

            label : str
                relative path to the run directory
            atoms : Atoms object
                the atoms object
            command  : str
                Command to run deMon. If not present the environment varable DEMON_COMMAND will be used
            restart  : str
                Relative path to ASE restart directory for parameters and atoms object and results
            basis_path  : str
                Relative path to the directory containing BASIS, AUXIS, ECPS, MCPS and AUGMENT
            ignore_bad_restart_file : bool
                Ignore broken or missing ASE restart files
                By default, it is an error if the restart
                file is missing or broken.
            deMon_restart_path  : str
                Relative path to the deMon restart dir
            title : str
                Title in the deMon input file.
            scftype : str
                Type of scf
            forces  : bool
                If True a force calculation will be enforced.
            dipole  : bool
                If True a dipole calculation will be enforced
            xc : str
                xc-functional
            guess : str
                guess for initial density and wave functions
            print_out : str | list
                Options for the printing in deMon
            basis : dict
                Definition of basis sets.
            ecps  : dict
                Definition of ECPs
            mcps  : dict
                Definition of MCPs
            auxis  : dict
                Definition of AUXIS
            augment : dict
                Definition of AUGMENT
            input_arguments : dict
                Explicitly given input arguments. The key is the input keyword
                and the value is either a str, a list of str (will be written on the same line as the keyword),
                or a list of lists of str (first list is written on the first line, the others on following lines.)

            For example usage, see the tests h2o.py and h2o_xas_xes.py in the directory ase/test/demon

            """

            parameters = Parameters_deMon(**kwargs)

            # Setup the run command
            command = parameters['command']
            if command is None:
                command = os.environ.get('DEMON_COMMAND')

            if command is None:
                mess = 'The "DEMON_COMMAND" environment is not defined.'
                raise ValueError(mess)
            else:
                parameters['command'] = command

            # Call the base class.
            FileIOCalculator.__init__(
                self,
                **parameters)

        def __getitem__(self, key):
            """Convenience method to retrieve a parameter as
            calculator[key] rather than calculator.parameters[key]

                Parameters:
                    key       : str, the name of the parameters to get.
            """
            return self.parameters[key]

        def set(self, **kwargs):
            """Set all parameters.

            Parameters:
                kwargs  : Dictionary containing the keywords for deMon
            """
            # Put in the default arguments.
            kwargs = self.default_parameters.__class__(**kwargs)

            if 'parameters' in kwargs:
                filename = kwargs.pop('parameters')
                parameters = Parameters.read(filename)
                parameters.update(kwargs)
                kwargs = parameters

            changed_parameters = {}

            for key, value in kwargs.items():
                oldvalue = self.parameters.get(key)
                if key not in self.parameters or not equal(value, oldvalue):
                    changed_parameters[key] = value
                    self.parameters[key] = value

            return changed_parameters

        def link_file(self, fromdir, todir, filename):

            if op.exists(todir + '/' + filename):
                os.remove(todir + '/' + filename)

            if op.exists(fromdir + '/' + filename):
                os.symlink(fromdir + '/' + filename,
                           todir + '/' + filename)
            else:
                raise RuntimeError(
                    "{0} doesn't exist".format(fromdir + '/' + filename))

        def calculate(self,
                      atoms=None,
                      properties=['energy'],
                      system_changes=all_changes):
            """Capture the RuntimeError from FileIOCalculator.calculate
            and add a little debug information from the deMon output.

            See base FileIocalculator for documentation.
            """

            if atoms is not None:
                self.atoms = atoms.copy()

            self.write_input(self.atoms, properties, system_changes)
            if self.command is None:
                raise RuntimeError('Please set $%s environment variable ' %
                                   ('DEMON_COMMAND') +
                                   'or supply the command keyword')
            command = self.command  # .replace('PREFIX', self.prefix)

            # basis path
            basis_path = self.parameters['basis_path']
            if basis_path is None:
                basis_path = os.environ.get('DEMON_BASIS_PATH')

            if basis_path is None:
                raise RuntimeError('Please set basis_path keyword,' +
                                   ' or the DEMON_BASIS_PATH' +
                                   ' environment variable')

            # link restart file
            value = self.parameters['guess']
            if value.upper() == 'RESTART':
                value2 = self.parameters['deMon_restart_path']

                if op.exists(self.directory + '/deMon.rst')\
                        or op.islink(self.directory + '/deMon.rst'):
                    os.remove(self.directory + '/deMon.rst')
                abspath = op.abspath(value2)

                if op.exists(abspath + '/deMon.mem') \
                        or op.islink(abspath + '/deMon.mem'):

                    shutil.copy(abspath + '/deMon.mem',
                                self.directory + '/deMon.rst')
                else:
                    raise RuntimeError(
                        "{0} doesn't exist".format(abspath + '/deMon.rst'))

            abspath = op.abspath(basis_path)

            for name in ['BASIS', 'AUXIS', 'ECPS', 'MCPS', 'FFDS']:
                self.link_file(abspath, self.directory, name)

            subprocess.check_call(command, shell=True, cwd=self.directory)

            try:
                self.read_results()
            except Exception:  # XXX Which kind of exception?
                with open(self.directory + '/deMon.out', 'r') as fd:
                    lines = fd.readlines()
                debug_lines = 10
                print('##### %d last lines of the deMon.out' % debug_lines)
                for line in lines[-20:]:
                    print(line.strip())
                print('##### end of deMon.out')
                raise RuntimeError

        def set_label(self, label):
            """Set label directory """

            self.label = label

            # in our case self.directory = self.label
            self.directory = self.label
            if self.directory == '':
                self.directory = os.curdir

        def write_input(self, atoms, properties=None, system_changes=None):
            """Write input (in)-file.
            See calculator.py for further details.

            Parameters:
                 atoms        : The Atoms object to write.
                 properties   : The properties which should be calculated.
                 system_changes : List of properties changed since last run.

            """
            # Call base calculator.
            FileIOCalculator.write_input(
                self,
                atoms=atoms,
                properties=properties,
                system_changes=system_changes)

            if system_changes is None and properties is None:
                return

            filename = self.label + '/deMon.inp'

            add_print = ''

            # Start writing the file.
            with open(filename, 'w') as fd:

                # write keyword argument keywords
                value = self.parameters['title']
                self._write_argument('TITLE', value, fd)

                fd.write('#\n')

                value = self.parameters['scftype']
                self._write_argument('SCFTYPE', value, fd)

                value = self.parameters['xc']
                self._write_argument('VXCTYPE', value, fd)

                value = self.parameters['guess']
                self._write_argument('GUESS', value, fd)

                # obtain forces through a single BOMD step
                # only if forces is in properties, or if keyword forces is True
                value = self.parameters['forces']
                if 'forces' in properties or value:

                    self._write_argument('DYNAMICS',
                                         ['INT=1', 'MAX=0', 'STEP=0'], fd)
                    self._write_argument('TRAJECTORY', 'FORCES', fd)
                    self._write_argument('VELOCITIES', 'ZERO', fd)
                    add_print = add_print + ' ' + 'MD OPT'

                # if dipole is True, enforce dipole calculation.
                # Otherwise only if asked for
                value = self.parameters['dipole']
                if 'dipole' in properties or value:
                    self._write_argument('DIPOLE', '', fd)

                # print argument, here other options could change this
                value = self.parameters['print_out']
                assert(type(value) is str)
                value = value + add_print

                if not len(value) == 0:
                    self._write_argument('PRINT', value, fd)
                    fd.write('#\n')

                # write general input arguments
                self._write_input_arguments(fd)

                fd.write('#\n')

                # write basis set, ecps, mcps, auxis, augment
                basis = self.parameters['basis']
                if 'all' not in basis:
                    basis['all'] = 'DZVP'
                self._write_basis(fd, atoms, basis, string='BASIS')

                ecps = self.parameters['ecps']
                if not len(ecps) == 0:
                    self._write_basis(fd, atoms, ecps, string='ECPS')

                mcps = self.parameters['mcps']
                if not len(mcps) == 0:
                    self._write_basis(fd, atoms, mcps, string='MCPS')

                auxis = self.parameters['auxis']
                if not len(auxis) == 0:
                    self._write_basis(fd, atoms, auxis, string='AUXIS')

                augment = self.parameters['augment']
                if not len(augment) == 0:
                    self._write_basis(fd, atoms, augment, string='AUGMENT')

                # write geometry
                self._write_atomic_coordinates(fd, atoms)

                # write xyz file for good measure.
                ase.io.write(self.label + '/deMon_atoms.xyz', self.atoms)

        def read(self, restart_path):
            """Read parameters from directory restart_path."""

            self.set_label(restart_path)

            if not op.exists(restart_path + '/deMon.inp'):
                raise ReadError('The restart_path file {0} does not exist'
                                .format(restart_path))

            self.atoms = self.deMon_inp_to_atoms(restart_path + '/deMon.inp')

            self.read_results()

        def _write_input_arguments(self, fd):
            """Write directly given input-arguments."""
            input_arguments = self.parameters['input_arguments']

            # Early return
            if input_arguments is None:
                return

            for key, value in input_arguments.items():
                self._write_argument(key, value, fd)

        def _write_argument(self, key, value, fd):
            """Write an argument to file.
             key :  a string coresponding to the input keyword
             value : the arguments, can be a string, a number or a list
             f :  and open file
            """

            # for only one argument, write on same line
            if not isinstance(value, (tuple, list)):
                line = key.upper()
                line += '    ' + str(value).upper()
                fd.write(line)
                fd.write('\n')

            # for a list, write first argument on the first line,
            # then the rest on new lines
            else:
                line = key
                if not isinstance(value[0], (tuple, list)):
                    for i in range(len(value)):
                        line += '  ' + str(value[i].upper())
                    fd.write(line)
                    fd.write('\n')
                else:
                    for i in range(len(value)):
                        for j in range(len(value[i])):
                            line += '  ' + str(value[i][j]).upper()
                        fd.write(line)
                        fd.write('\n')
                        line = ''

        def _write_atomic_coordinates(self, fd, atoms):
            """Write atomic coordinates.

            Parameters:
            - f:     An open file object.
            - atoms: An atoms object.
            """

            fd.write('#\n')
            fd.write('# Atomic coordinates\n')
            fd.write('#\n')
            fd.write('GEOMETRY CARTESIAN ANGSTROM\n')

            for i in range(len(atoms)):
                xyz = atoms.get_positions()[i]
                chem_symbol = atoms.get_chemical_symbols()[i]
                chem_symbol += str(i + 1)

                # if tag is set to 1 then we have a ghost atom,
                # set nuclear charge to 0
                if(atoms.get_tags()[i] == 1):
                    nuc_charge = str(0)
                else:
                    nuc_charge = str(atoms.get_atomic_numbers()[i])

                mass = atoms.get_masses()[i]

                line = '{0:6s}'.format(chem_symbol).rjust(10) + ' '
                line += '{0:.5f}'.format(xyz[0]).rjust(10) + ' '
                line += '{0:.5f}'.format(xyz[1]).rjust(10) + ' '
                line += '{0:.5f}'.format(xyz[2]).rjust(10) + ' '
                line += '{0:5s}'.format(nuc_charge).rjust(10) + ' '
                line += '{0:.5f}'.format(mass).rjust(10) + ' '

                fd.write(line)
                fd.write('\n')

        # routine to write basis set inormation, including ecps and auxis
        def _write_basis(self, fd, atoms, basis={}, string='BASIS'):
            """Write basis set, ECPs, AUXIS, or AUGMENT basis

            Parameters:
            - f:     An open file object.
            - atoms: An atoms object.
            - basis: A dictionary specifying the basis set
            - string: 'BASIS', 'ECP','AUXIS' or 'AUGMENT'
            """

            # basis for all atoms
            line = '{0}'.format(string).ljust(10)

            if 'all' in basis:
                default_basis = basis['all']
                line += '({0})'.format(default_basis).rjust(16)

            fd.write(line)
            fd.write('\n')

            # basis for all atomic species
            chemical_symbols = atoms.get_chemical_symbols()
            chemical_symbols_set = set(chemical_symbols)

            for i in range(chemical_symbols_set.__len__()):
                symbol = chemical_symbols_set.pop()

                if symbol in basis:
                    line = '{0}'.format(symbol).ljust(10)
                    line += '({0})'.format(basis[symbol]).rjust(16)
                    fd.write(line)
                    fd.write('\n')

            # basis for individual atoms
            for i in range(len(atoms)):

                if i in basis:
                    symbol = str(chemical_symbols[i])
                    symbol += str(i + 1)

                    line = '{0}'.format(symbol).ljust(10)
                    line += '({0})'.format(basis[i]).rjust(16)
                    fd.write(line)
                    fd.write('\n')

        # Analysis routines
        def read_results(self):
            """Read the results from output files."""
            self.read_energy()
            self.read_forces(self.atoms)
            self.read_eigenvalues()
            self.read_dipole()
            self.read_xray()

        def read_energy(self):
            """Read energy from deMon's text-output file."""
            with open(self.label + '/deMon.out', 'r') as fd:
                text = fd.read().upper()

            lines = iter(text.split('\n'))

            for line in lines:
                if line.startswith(' TOTAL ENERGY                ='):
                    self.results['energy'] = float(line.split()[-1]) * Hartree
                    break
            else:
                raise RuntimeError

        def read_forces(self, atoms):
            """Read the forces from the deMon.out file."""

            natoms = len(atoms)
            filename = self.label + '/deMon.out'

            if op.isfile(filename):
                with open(filename, 'r') as fd:
                    lines = fd.readlines()

                    # find line where the orbitals start
                    flag_found = False
                    for i in range(len(lines)):
                        if lines[i].rfind('GRADIENTS OF TIME STEP 0 IN A.U.') > -1:
                            start = i + 4
                            flag_found = True
                            break

                    if flag_found:
                        self.results['forces'] = np.zeros((natoms, 3), float)
                        for i in range(natoms):
                            line = [s for s in lines[i + start].strip().split(' ')
                                    if len(s) > 0]
                            f = -np.array([float(x) for x in line[2:5]])
                            self.results['forces'][i, :] = f * (Hartree / Bohr)

        def read_eigenvalues(self):
            """Read eigenvalues from the 'deMon.out' file."""
            assert os.access(self.label + '/deMon.out', os.F_OK)

            # Read eigenvalues
            with open(self.label + '/deMon.out', 'r') as fd:
                lines = fd.readlines()

            # try  PRINT MOE
            eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(
                lines, 'ALPHA MO ENERGIES', 6)
            eig_beta, occ_beta = self.read_eigenvalues_one_spin(
                lines, 'BETA MO ENERGIES', 6)

            # otherwise try PRINT MOS
            if len(eig_alpha) == 0 and len(eig_beta) == 0:
                eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(
                    lines, 'ALPHA MO COEFFICIENTS', 5)
                eig_beta, occ_beta = self.read_eigenvalues_one_spin(
                    lines, 'BETA MO COEFFICIENTS', 5)

            self.results['eigenvalues'] = np.array([eig_alpha, eig_beta]) * Hartree
            self.results['occupations'] = np.array([occ_alpha, occ_beta])

        def read_eigenvalues_one_spin(self, lines, string, neigs_per_line):
            """Utility method for retreiving eigenvalues after the string "string"
            with neigs_per_line eigenvlaues written per line
            """
            eig = []
            occ = []

            skip_line = False
            more_eigs = False

            # find line where the orbitals start
            for i in range(len(lines)):
                if lines[i].rfind(string) > -1:
                    ii = i
                    more_eigs = True
                    break

            while more_eigs:
                # search for two empty lines in a row preceding a line with
                # numbers
                for i in range(ii + 1, len(lines)):
                    if len(lines[i].split()) == 0 and \
                            len(lines[i + 1].split()) == 0 and \
                            len(lines[i + 2].split()) > 0:
                        ii = i + 2
                        break

                # read eigenvalues, occupations
                line = lines[ii].split()
                if len(line) < neigs_per_line:
                    # last row
                    more_eigs = False
                if line[0] != str(len(eig) + 1):
                    more_eigs = False
                    skip_line = True

                if not skip_line:
                    line = lines[ii + 1].split()
                    for l in line:
                        eig.append(float(l))
                    line = lines[ii + 3].split()
                    for l in line:
                        occ.append(float(l))
                    ii = ii + 3

            return eig, occ

        def read_dipole(self):
            """Read dipole moment."""
            dipole = np.zeros(3)
            with open(self.label + '/deMon.out', 'r') as fd:
                lines = fd.readlines()

                for i in range(len(lines)):
                    if lines[i].rfind('DIPOLE') > -1 and lines[i].rfind('XAS') == -1:
                        dipole[0] = float(lines[i + 1].split()[3])
                        dipole[1] = float(lines[i + 2].split()[3])
                        dipole[2] = float(lines[i + 3].split()[3])

                        # debye to e*Ang
                        self.results['dipole'] = dipole * 0.2081943482534

                        break

        def read_xray(self):
            """Read deMon.xry if present."""

            # try to read core IP from, .out file
            filename = self.label + '/deMon.out'
            core_IP = None
            if op.isfile(filename):
                with open(filename, 'r') as fd:
                    lines = fd.readlines()

                for i in range(len(lines)):
                    if lines[i].rfind('IONIZATION POTENTIAL') > -1:
                        core_IP = float(lines[i].split()[3])

            try:
                mode, ntrans, E_trans, osc_strength, trans_dip = parse_xray(self.label + '/deMon.xry')
            except ReadError:
                pass
            else:
                xray_results = {'xray_mode': mode,
                                'ntrans': ntrans,
                                'E_trans': E_trans,
                                'osc_strength': osc_strength,  # units?
                                'trans_dip': trans_dip,  # units?
                                'core_IP': core_IP}

                self.results['xray'] = xray_results

        def deMon_inp_to_atoms(self, filename):
            """Routine to read deMon.inp and convert it to an atoms object."""

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

            # find line where geometry starts
            for i in range(len(lines)):
                if lines[i].rfind('GEOMETRY') > -1:
                    if lines[i].rfind('ANGSTROM'):
                        coord_units = 'Ang'
                    elif lines.rfind('Bohr'):
                        coord_units = 'Bohr'
                    ii = i
                    break

            chemical_symbols = []
            xyz = []
            atomic_numbers = []
            masses = []

            for i in range(ii + 1, len(lines)):
                try:
                    line = lines[i].split()

                    if(len(line) > 0):
                        for symbol in ase.data.chemical_symbols:
                            found = None
                            if line[0].upper().rfind(symbol.upper()) > -1:
                                found = symbol
                                break

                            if found is not None:
                                chemical_symbols.append(found)
                            else:
                                break

                            xyz.append([float(line[1]), float(line[2]), float(line[3])])

                    if len(line) > 4:
                        atomic_numbers.append(int(line[4]))

                    if len(line) > 5:
                        masses.append(float(line[5]))

                except Exception:  # XXX Which kind of exception?
                    raise RuntimeError

            if coord_units == 'Bohr':
                xyz = xyz * Bohr

            natoms = len(chemical_symbols)

            # set atoms object
            atoms = ase.Atoms(symbols=chemical_symbols, positions=xyz)

            # if atomic numbers were read in, set them
            if(len(atomic_numbers) == natoms):
                atoms.set_atomic_numbers(atomic_numbers)

            # if masses were read in, set them
            if(len(masses) == natoms):
                atoms.set_masses(masses)

            return atoms