Computing density on grid from density matrix in .DM file

[rreho]

I was trying to compute the density on a grid by reading the information from a .DM file that was the output of a SIESTA calculation (including SOC) but I am running into issues. First the basis information was missing but even after manually supplying that by reading it from the ion file and attaching the proper geometry object I still get the following error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/var/folders/hx/z4tmr0ln4bdct36z3dr7ppt80000gp/T/ipykernel_24460/1110118028.py in <module>
----> 1 ChiDM.density(grid1)

~/opt/anaconda3/envs/siesta/lib/python3.10/site-packages/sisl/physics/densitymatrix.py in density(self, grid, spinor, tol, eta)
    714 
    715                     # Calculate the psi_j component
--> 716                     psi = o.psi_spher(ja_r1, ja_theta1, ja_cos_phi1, cos_phi=True)
    717 
    718                     # Now add this orbital to all components

~/opt/anaconda3/envs/siesta/lib/python3.10/site-packages/sisl/orbital.py in psi_spher(self, r, theta, phi, cos_phi)
   1101              basis function value at point `r`
   1102         """
-> 1103         return self.orb.psi_spher(r, theta, phi, self.m, cos_phi)
   1104 
   1105     def __getstate__(self):

AttributeError: 'Orbital' object has no attribute 'psi_spher'

Do you know a way to be able to do this?

[pfebrer]

Hi! How do you read the basis?

[pfebrer]

I would suggest you read the density matrix from the FDF file to avoid problems (if possible):

DM = sisl.get_sile("myfile.fdf").read_density_matrix()

[rreho]

I would suggest you read the density matrix from the FDF file to avoid problems (if possible):

DM = sisl.get_sile("myfile.fdf").read_density_matrix()

If we try this, we still get the error there is no orbital information.

Hi! How do you read the basis?

We read the basis like this:

pb_ion = sisl.io.siesta.ionxmlSileSiesta('Pb.ion.xml')

pb_atom =pb_ion.read_basis()

geom = sisl.get_sile('Pb.fdf').read_geometry()
geom_pb =  sisl.Geometry([0,0,0], atoms = pb_atom)
geom.replace(np.arange(geom.na), geom_pb)

DM = sisl.get_sile('bulk-Pb.DM').read_density_matrix(geometry = geom)

[pfebrer]

Ok, Geometry.replace returns a new geometry. So you should do:

geom = geom.replace(np.arange(geom.na), geom_pb)

But still, I don't know why it can not read the basis properly from the fdf file when you call sisl.get_sile("myfile.fdf").read_density_matrix() (without any extra arguments). That's strange.

[pfebrer]

But I think even

geom = sisl.get_sile('Pb.fdf').read_geometry()

should get you the geometry with the correct basis if the basis files are there (in the same directory as the fdf file). Could you share the output of print(geom) after this operation?

[rreho]

It turns out that the correct ion.xml files were missing from the directory. Now it works, thank you!