SMTG-Bham / sumo

Heavyweight plotting tools for ab initio calculations
https://smtg-bham.github.io/sumo/
MIT License
192 stars 78 forks source link

Combination of certain sites of elements is not support in SUMO? #75

Open Nokimann opened 5 years ago

Nokimann commented 5 years ago

Let say the structure is

site / element 1 / A 2 / A 3 / B 4 / B 5 / C 6 / C

Could I extract figure and data as divided two parts: sum(A.1,B.3,C.5) & sum(A.2,B.4,C.6) on the command line? Or, should I create my own script manually?

ajjackson commented 5 years ago

We've discussed this idea in the past and basically we would need to

So it would be good to do but is quite a bit of work to do well. We're interested but it won't happen soon :-( I think your best bet is to run dosplot using something like --atoms A.1,B.1,C.1, copy the .dat files to another name, re-run dosplot with --atoms A.2,B.2,C.2, and then make your own plot summing data from the .dat files.

Nokimann commented 5 years ago

On your answer, 3 / B, 4 / B, 5 / C, and 6 / C correspond to B.1, B.2, C.1, and C.2, respectively. Right?

What about using either '{}', '[]', '()', or '<>' ? Also, I suggest a range of sites, e.g., 1-5 is same as 1, 2, 3, 4, and 5 because someone might consider large structures (~100 atoms)

And we may sum the individual densities object with reduce(add_densities, densities). Here add_densities is a function from pymatgen and 'densities' is a list of densities

ajjackson commented 5 years ago

Yes, I think the counting restarts at one for each element type.

If you are comfortable working with the Pymatgen objects then you might as well write a script that uses Pymatgen directly to get a pymatgen.electronic_structure.dos you are happy with. In principle you could feed this back into a sumo.plotting.dos_plotter.SDOSPlotter to benefit from some of Sumo's effort on plotting outputs.

Most of those brackets have special meanings in Bash which would be likely to trip up users. Certainly ranges should be supported. The challenge is to enable selection, grouping, splitting and plotting over multiple axes by any combination of atoms, elements, orbital and spin in notation that is concise, legible and unambiguous. It's not a simple problem, but it beats what we have done so far which is to individually code logic for all the arrangements that we encountered routinely.

Nokimann commented 5 years ago

I did change some codes in dos_plotter.py, dosplot.py and dos.py to create the summations of atoms.

@ajjackson This is my command line, and I hope this is helpful to give SUMO an idea of language for users in the future. sumo-dosplot --bundle-atoms Mo.1,S.1_2/Mo.2,S.3_4

'_' indicates the range of sites '/' means a partition

By the way, could you advise me with a guide about modifying sumo-bandplot? I'm appreciated for your previous guide.

Below are some parts that I changed. Sorry for a messy organization. I think the important parts are '_bundle_atoms' and 'get_nk_pdos' functions

In short, I replaced --orbitals, --elements, and --atoms with one --bundle-atoms flag. I created a new type function: bundleatoms in getparser. I modified a function: _loaddos with a new _get_nkpdos for creating the summation of atoms I removed _sortorbitals in _writefiles and _dos_plotdata

In dosplot.py

def dosplot(filename=None, code='vasp', prefix=None, directory=None,
            bundle_atoms=None, subplot=False,
            shift=True, total_only=False, plot_total=True, legend_on=True,
            legend_frame_on=False, legend_cutoff=3., gaussian=None, height=6.,
            width=8., xmin=-6., xmax=6., num_columns=2, colours=None, yscale=1,
            xlabel='Energy (eV)', ylabel='Arb. units',
            style=None, no_base_style=False,
            image_format='pdf', dpi=400, plt=None, fonts=None):

    if code.lower() == 'vasp':
        if not filename:
            if os.path.exists('vasprun.xml'):
                filename = 'vasprun.xml'
            elif os.path.exists('vasprun.xml.gz'):
                filename = 'vasprun.xml.gz'
            else:
                logging.error('ERROR: No vasprun.xml found!')
                sys.exit()

        dos, pdos = load_dos(filename, bundle_atoms, gaussian,
                             total_only)
    ...

def _bundle_atoms(bundle_atoms_string):
    list_atoms = []
    for bundle in bundle_atoms_string.split('/'):
        atoms = {}
        for split in bundle.split(','):
            sites = split.split('.')
            el = sites.pop(0)
            for i, s in list(enumerate(sites)):
                if '_' in s:
                    start, before_end = tuple(map(int, s.split('_')))
                    sites.pop(i)
                    sites.extend(range(start, before_end + 1))
            sites = sorted(list(map(int, sites)))
            atoms[el] = np.array(sites) - 1
        list_atoms.append(atoms)
    print(list_atoms)

def _get_parser():
    ...
    parser.add_argument('-ba', '--bundle-atoms', type=_bundle_atoms, metavar='BA',
                        help=('atoms to include (e.g. "Mo.1,S.1_2/Mo.2,S.3_4")'))
    ...

In dos.py

def load_dos(vasprun, bundle_atoms=None,
    ...
    pdos = {}
    if not total_only:
        pdos = get_nk_pdos(dos, bundle_atoms=bundle_atoms)
    return dos, pdos

def get_nk_pdos(dos, bundle_atoms=None):
    from collections import defaultdict
    from functools import reduce

    from pymatgen.electronic_structure.dos import Dos, add_densities

    d = defaultdict(lambda : np.arange(0))
    for s in dos.structure.sites:
        d[s.specie] = np.append(d[s.specie], s)

    pdos = defaultdict(dict)
    for i, bundle in enumerate(bundle_atoms):
        bundle_densities = []
        for el, sites in bundle.items():
            map_densities = map(lambda s: dos.get_site_dos(s).densities,
                                d[get_el_sp(el)][sites])
            bundle_densities.extend(map_densities)
        densities = reduce(add_densities, bundle_densities)
        bundle_dos = Dos(dos.efermi, dos.energies, densities)
        info = ', '.join('{}-{}'.format(el, sites)
                         for el, sites in bundle.items())
        pdos['sum ' + str(i)][info] = bundle_dos
    return pdos

def write_files(dos, pdos, prefix=None, directory=None, zero_to_efermi=True):
    ...
    spin = len(dos.densities)
    for el, el_pdos in pdos.items():
        header = ['energy']
        pdos_data = [eners]
        for orb in el_pdos: # <-- just remove sort_orbitals function
            for spin, sign, label in sdata:
                header.append('{}{}'.format(orb, label))
                pdos_data.append(el_pdos[orb].densities[spin] * sign)
        pdos_data = np.stack(pdos_data, axis=1)
    ...

In dos_plotter.py

class SDOSPlotter(object):
    ...
    def dos_plot_data(self, yscale=1, xmin=-6., xmax=6., colours=None,
                      plot_total=True, legend_cutoff=3, subplot=False,
                      zero_to_efermi=True, cache=None):
        ...
        for el, el_pdos in pdos.items():
            el_lines = []
            for orb in el_pdos: # <-- I just remove sort_orbitals function
                dmax = max([max(d[mask])
                            for d in el_pdos[orb].densities.values()])
                ymax = dmax if dmax > ymax else ymax
                label = None if dmax < cutoff else '{} ({})'.format(el, orb)
                colour, cache = get_cached_colour(el, orb, colours,
                                                  cache=cache)
                el_lines.append({'label': label, 'alpha': 0.25,
                                 'colour': colour,
                                 'dens': el_pdos[orb].densities})
            if subplot:
                lines.append(el_lines)
            else:
                lines[0].extend(el_lines)
        ...
ajjackson commented 5 years ago

I'm less familiar with the orbital projected band structures and I know there are some challenges there. @utf may be able to advise on a workaround for your specific goal.

Your changes to dosplot look to be along the right general lines, but I would prefer not to bring them into Sumo until we are able to dedicate a bit more time to the problem. For now I suggest you work with these modifications that suit your needs. You may find it helpful to maintain your own fork of Sumo so you can merge in other updates while keeping these changes?

I appreciate this contribution and it seems to be a good start on the syntax problem.

Nokimann commented 4 years ago

Hope to be added