Trouble loading geometry

[wladerer]

Hi everyone, my goal is relatively simple:

Plot the atomic drift due to a DFT geometry optimization

I can load in all the data I need to calculate the displacement vectors for each atom using ASE. What I cannot figure out is creating the appropriate geometry object.

I tried the VASP IO module

geom = carSile('CONTCAR').get_geometry()
geom.plot()

Which did not work.

Any clues on how to use the io modules to instantiate a geometry object ?

[zerothi]

Hi, generally one should trust sisl to parse the file correctly, I think it should work with:

import sisl
geom = sisl.get_sile('CONTCAR').read_geometry()

there are no get_geometry methods, for details please have a look at the documentation here

Thanks!

[pfebrer]

Hi @wladerer, could you print geom and geom.xyz?

[pfebrer]

Also, what happens if you do geom.plot(axes="xy")?

Maybe you can upload the CONTCAR here so that I can try, it will probably be faster :)

[wladerer]

geom.plot(axes="xy") worked, I think! It opens a tab in my browser and produces an image. I'll send the file anyway. It is a surface slab.

contcar.txt

For printing geom and geom.xyz

Geometry{na: 36, no: 36,
 Atoms{species: 2,
  Atom{Sn, Z: 50, mass(au): 118.71000, maxR: -1.00000,
   Orbital{R: -1.00000, q0: 0.0}
  }: 9,
  Atom{Pt, Z: 78, mass(au): 195.08400, maxR: -1.00000,
   Orbital{R: -1.00000, q0: 0.0}
  }: 27,
 },
 maxR: -1.00000,
 Lattice{nsc: [1 1 1],
  origin=0.0000, 0.0000, 0.0000,
  A=[4.0293, 0.0000, 0.0000],
  B=[0.0000, 4.0293, 0.0000],
  C=[0.0000, 0.0000, 48.3513],
  bc=[Periodic,
      Periodic,
      Periodic]
 }
}
[[0.00000000e+00 0.00000000e+00 9.00020381e+00]
 [0.00000000e+00 0.00000000e+00 1.30951560e+01]
 [0.00000000e+00 0.00000000e+00 1.71244347e+01]
 [0.00000000e+00 0.00000000e+00 2.11537135e+01]
 [0.00000000e+00 0.00000000e+00 2.51829923e+01]
 [0.00000000e+00 0.00000000e+00 2.92122710e+01]
 [0.00000000e+00 0.00000000e+00 3.32415498e+01]
 [0.00000000e+00 0.00000000e+00 3.73078952e+01]
 [0.00000000e+00 0.00000000e+00 4.15209128e+01]
 [2.01463938e+00 2.01463938e+00 9.10018342e+00]
 [2.01463938e+00 0.00000000e+00 7.13402416e+00]
 [3.00000000e-16 2.01463938e+00 7.13402416e+00]
 [2.01463938e+00 2.01463938e+00 1.30951560e+01]
 [2.01463938e+00 0.00000000e+00 1.10770654e+01]
 [3.00000000e-16 2.01463938e+00 1.10770654e+01]
 [2.01463938e+00 2.01463938e+00 1.71244347e+01]
 [2.01463938e+00 0.00000000e+00 1.51097953e+01]
 [3.00000000e-16 2.01463938e+00 1.51097953e+01]
 [2.01463938e+00 2.01463938e+00 2.11537135e+01]
 [2.01463938e+00 0.00000000e+00 1.91390741e+01]
 [3.00000000e-16 2.01463938e+00 1.91390741e+01]
 [2.01463938e+00 2.01463938e+00 2.51829923e+01]
 [2.01463938e+00 0.00000000e+00 2.31683529e+01]
 [3.00000000e-16 2.01463938e+00 2.31683529e+01]
 [2.01463938e+00 2.01463938e+00 2.92122710e+01]
 [2.01463938e+00 0.00000000e+00 2.71976316e+01]
 [3.00000000e-16 2.01463938e+00 2.71976316e+01]
 [2.01463938e+00 2.01463938e+00 3.32415498e+01]
 [2.01463938e+00 0.00000000e+00 3.12269104e+01]
 [3.00000000e-16 2.01463938e+00 3.12269104e+01]
 [2.01463938e+00 2.01463938e+00 3.72729057e+01]
 [2.01463938e+00 0.00000000e+00 3.52561891e+01]
 [3.00000000e-16 2.01463938e+00 3.52561891e+01]
 [2.01463938e+00 2.01463938e+00 4.12367387e+01]
 [2.01463938e+00 0.00000000e+00 3.93089580e+01]
 [3.00000000e-16 2.01463938e+00 3.93089580e+01]]

Thank you!

[pfebrer]

Ok, thanks! I can reproduce the problem, I'll try to fix it tomorrow.

There must be something unique in this geometry that triggers a hidden bug because I can do

import sisl
import sisl.viz

sisl.geom.graphene().write("CONTCAR")
sisl.get_sile("CONTCAR").read_geometry().plot()

And I get the 3D plot without problems 😅

[pfebrer]

This is fixed by: #681. Thanks for the report!

Apparently no one had generated a 3D plot with different species that have the same color 😅

By the way the default colors are defined here:

sisl/src/sisl/viz/data_sources/atom_data.py

Lines 68 to 87 in ab45ef4

	class AtomDefaultColors(AtomData):
	_atoms_colors = {
	"H": "#cccccc",
	"O": "red",
	"Cl": "green",
	"N": "blue",
	"C": "grey",
	"S": "yellow",
	"P": "orange",
	"Au": "gold",
	"else": "pink",
	}
	def function(self, geometry, atoms=None):
	return np.array(
	[
	self._atoms_colors.get(atom.symbol, self._atoms_colors["else"])
	for atom in geometry.sub(atoms).atoms
	]
	)

I did this at the very beginning and it does not make sense to show two different species in pink, but I had no idea what to use. Any suggestions of some default colors that I can take from some place? :)

[wladerer]

I really like the Jmol colors. I have implemented my own mini visualization suite in plotly using the Jmol colors and it works well. I think ASE uses that for their coloring scheme

def plot_atomic_drift(initial_file: str, final_file: str, output_file: str = None):
    '''Plots the atomic drift between two structures'''
    from ase.io import read
    from ase.data.colors import jmol_colors, cpk_colors
    from ase.data import covalent_radii
    import numpy as np
    import plotly.graph_objects as go

    initial = read(initial_file)
    final = read(final_file)

    #atom types 
    atomic_numbers = initial.get_atomic_numbers()
    atom_types = initial.get_chemical_symbols()
    atomic_radii = [ covalent_radii[number] for number in atomic_numbers ]
    atomic_colors = [ jmol_colors[number] for number in atomic_numbers ]

etc. etc.

[pfebrer]

Cool! So I just looked at ase's source and they take it from here: https://jmol.sourceforge.net/jscolors/#color_U

[wladerer]

I know it's not my place to suggest things, but if you expand it to rgba rather than rgb, then you have an additional handle on opacity for plotting purposes which might improve backward and forward compatibility?

[pfebrer]

Suggestions are most welcome :)

Opacity is controlled by a separate opacity key. See: https://zerothi.github.io/sisl/visualization/viz_module/showcase/GeometryPlot.html#Custom-styles-for-atoms.

This is so that it is independent of color and is easier to handle. For example you can set the color to a color name and still use opacity. It would be very hard to handle if there are multiple places where opacity can be tweaked. And also some backends accept rgba and some don't. I think blender didn't for example. So the bottomline is that it's easier to keep opacity separate 😅

[pfebrer]

Thanks for the discussion! Jmol colours are now the default, so your structure should look better by default :)

By the way, any suggestion on the documentation would be highly appreciated 👍